Skip to content

Commit

Permalink
Refactor reserves (#562)
Browse files Browse the repository at this point in the history 
This fixes #573, where Reserves + Num_VRE_Bins > 1 would lead to the
VRE_NO_POWER_OUT resources being able to contribute an infinite amount
to the reserves.

This fixes #574, where VRE power output constraints were being applied
length(VRE) times to all VRE resources, when Reserves = 1.

This fixed #566, where some constraints in storage_symmetric were redundant.
---------

Co-authored-by: Luca Bonaldo <[email protected]>
  • Loading branch information
@cfe316 @lbonaldo
cfe316 and lbonaldo authored Nov 18, 2023
1 parent a7ae6ba commit 2ef54dc
Show file tree
Hide file tree
Showing 9 changed files with 282 additions and 492 deletions.
29 changes: 28 additions & 1 deletion src/model/expression_manipulation.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,6 +47,33 @@ function fill_with_const!(arr::Array{GenericAffExpr{C,T}, dims}, con::Real) wher
return nothing
end

###### ###### ###### ###### ###### ######
# Create an expression from some first-dimension indices of a 2D variable array,
# where all of the 2nd-dimension indices are kept
###### ###### ###### ###### ###### ######
#
function extract_time_series_to_expression(var::Matrix{VariableRef},
set::AbstractVector{Int})
TIME_DIM = 2
time_range = 1:size(var)[TIME_DIM]

aff_exprs_data = AffExpr.(0, var[set, :] .=> 1)
new_axes = (set, time_range)
expr = JuMP.Containers.DenseAxisArray(aff_exprs_data, new_axes...)
return expr
end

function extract_time_series_to_expression(var::JuMP.Containers.DenseAxisArray{VariableRef, 2, Tuple{X, Base.OneTo{Int64}}, Y},
set::AbstractVector{Int}) where {X, Y}
TIME_DIM = 2
time_range = var.axes[TIME_DIM]

aff_exprs = AffExpr.(0, var[set, :] .=> 1)
new_axes = (set, time_range)
expr = JuMP.Containers.DenseAxisArray(aff_exprs.data, new_axes...)
return expr
end

###### ###### ###### ###### ###### ######
# Element-wise addition of one expression into another
# Both arrays must have the same dimensions
Expand Down Expand Up @@ -141,4 +168,4 @@ end

function sum_expression(expr::AbstractArray{AbstractJuMPScalar, dims}) where {dims}
return _sum_expression(expr)
end
end
131 changes: 62 additions & 69 deletions src/model/resources/curtailable_variable_renewable/curtailable_variable_renewable.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,24 +47,24 @@ function curtailable_variable_renewable!(EP::Model, inputs::Dict, setup::Dict)
add_similar_to_expression!(EP[:eCapResMarBalance], eCapResMarBalanceVRE)
end

### Constratints ###
# For resource for which we are modeling hourly power output
for y in VRE_POWER_OUT
# Define the set of generator indices corresponding to the different sites (or bins) of a particular VRE technology (E.g. wind or solar) in a particular zone.
# For example the wind resource in a particular region could be include three types of bins corresponding to different sites with unique interconnection, hourly capacity factor and maximim available capacity limits.
VRE_BINS = intersect(dfGen[dfGen[!,:R_ID].>=y,:R_ID], dfGen[dfGen[!,:R_ID].<=y+dfGen[y,:Num_VRE_Bins]-1,:R_ID])

# Constraints on contribution to regulation and reserves
if Reserves == 1
curtailable_variable_renewable_reserves!(EP, inputs)
else
# Maximum power generated per hour by renewable generators must be less than
# sum of product of hourly capacity factor for each bin times its the bin installed capacity
# Note: inequality constraint allows curtailment of output below maximum level.
@constraint(EP, [t=1:T], EP[:vP][y,t] <= sum(inputs["pP_Max"][yy,t]*EP[:eTotalCap][yy] for yy in VRE_BINS))
end

end
### Constraints ###
if Reserves == 1
# Constraints on power output and contribution to regulation and reserves
curtailable_variable_renewable_reserves!(EP, inputs)
remove_reserves_for_binned_vre_resources!(EP, inputs)
else
# For resource for which we are modeling hourly power output
for y in VRE_POWER_OUT
# Define the set of generator indices corresponding to the different sites (or bins) of a particular VRE technology (E.g. wind or solar) in a particular zone.
# For example the wind resource in a particular region could be include three types of bins corresponding to different sites with unique interconnection, hourly capacity factor and maximim available capacity limits.
VRE_BINS = intersect(dfGen[dfGen[!,:R_ID].>=y,:R_ID], dfGen[dfGen[!,:R_ID].<=y+dfGen[y,:Num_VRE_Bins]-1,:R_ID])

# Maximum power generated per hour by renewable generators must be less than
# sum of product of hourly capacity factor for each bin times its the bin installed capacity
# Note: inequality constraint allows curtailment of output below maximum level.
@constraint(EP, [t=1:T], EP[:vP][y,t] <= sum(inputs["pP_Max"][yy,t]*EP[:eTotalCap][yy] for yy in VRE_BINS))
end
end

# Set power variables for all bins that are not being modeled for hourly output to be zero
for y in VRE_NO_POWER_OUT
Expand Down Expand Up @@ -102,60 +102,53 @@ The amount of frequency regulation and operating reserves procured in each time
```
"""
function curtailable_variable_renewable_reserves!(EP::Model, inputs::Dict)

dfGen = inputs["dfGen"]
T = inputs["T"]

VRE_POWER_OUT = intersect(dfGen[dfGen.Num_VRE_Bins.>=1,:R_ID], inputs["VRE"])

for y in VRE_POWER_OUT
# Define the set of generator indices corresponding to the different sites (or bins) of a particular VRE technology (E.g. wind or solar) in a particular zone.
# For example the wind resource in a particular region could be include three types of bins corresponding to different sites with unique interconnection, hourly capacity factor and maximim available capacity limits.
VRE_BINS = intersect(dfGen[dfGen[!,:R_ID].>=y,:R_ID], dfGen[dfGen[!,:R_ID].<=y+dfGen[y,:Num_VRE_Bins]-1,:R_ID])

if y in inputs["REG"] && y in inputs["RSV"] # Resource eligible for regulation and spinning reserves
@constraints(EP, begin
# For VRE, reserve contributions must be less than the specified percentage of the capacity factor for the hour
[t=1:T], EP[:vREG][y,t] <= dfGen[y,:Reg_Max]*sum(inputs["pP_Max"][yy,t]*EP[:eTotalCap][yy] for yy in VRE_BINS)
[t=1:T], EP[:vRSV][y,t] <= dfGen[y,:Rsv_Max]*sum(inputs["pP_Max"][yy,t]*EP[:eTotalCap][yy] for yy in VRE_BINS)

# Power generated and regulation reserve contributions down per hour must be greater than zero
[t=1:T], EP[:vP][y,t]-EP[:vREG][y,t] >= 0

# Power generated and reserve contributions up per hour by renewable generators must be less than
# hourly capacity factor. Note: inequality constraint allows curtailment of output below maximum level.
[t=1:T], EP[:vP][y,t]+EP[:vREG][y,t]+EP[:vRSV][y,t] <= sum(inputs["pP_Max"][yy,t]*EP[:eTotalCap][yy] for yy in VRE_BINS)
end)
elseif y in inputs["REG"] # Resource only eligible for regulation reserves
@constraints(EP, begin
# For VRE, reserve contributions must be less than the specified percentage of the capacity factor for the hour
[t=1:T], EP[:vREG][y,t] <= dfGen[y,:Reg_Max]*sum(inputs["pP_Max"][yy,t]*EP[:eTotalCap][yy] for yy in VRE_BINS)

# Power generated and regulation reserve contributions down per hour must be greater than zero
[t=1:T], EP[:vP][y,t]-EP[:vREG][y,t] >= 0

# Power generated and reserve contributions up per hour by renewable generators must be less than
# hourly capacity factor. Note: inequality constraint allows curtailment of output below maximum level.
[t=1:T], EP[:vP][y,t]+EP[:vREG][y,t] <= sum(inputs["pP_Max"][yy,t]*EP[:eTotalCap][yy] for yy in VRE_BINS)
end)

elseif y in inputs["RSV"] # Resource only eligible for spinning reserves - only available in up, no down spinning reserves

@constraints(EP, begin
# For VRE, reserve contributions must be less than the specified percentage of the capacity factor for the hour
[t=1:T], EP[:vRSV][y,t] <= dfGen[y,:Rsv_Max]*sum(inputs["pP_Max"][yy,t]*EP[:eTotalCap][yy] for yy in VRE_BINS)

# Power generated and reserve contributions up per hour by renewable generators must be less than
# hourly capacity factor. Note: inequality constraint allows curtailment of output below maximum level.
[t=1:T], EP[:vP][y,t]+EP[:vRSV][y,t] <= sum(inputs["pP_Max"][yy,t]*EP[:eTotalCap][yy] for yy in VRE_BINS)
end)
else # Resource not eligible for reserves
# Maximum power generated per hour by renewable generators must be less than
# sum of product of hourly capacity factor for each bin times its the bin installed capacity
# Note: inequality constraint allows curtailment of output below maximum level.
@constraint(EP, [t=1:T], EP[:vP][y,t] <= sum(inputs["pP_Max"][yy,t]*EP[:eTotalCap][yy] for yy in VRE_BINS))
end
end
VRE = inputs["VRE"]
VRE_POWER_OUT = intersect(dfGen[dfGen.Num_VRE_Bins.>=1,:R_ID], VRE)
REG = intersect(VRE_POWER_OUT, inputs["REG"])
RSV = intersect(VRE_POWER_OUT, inputs["RSV"])

eTotalCap = EP[:eTotalCap]
vP = EP[:vP]
vREG = EP[:vREG]
vRSV = EP[:vRSV]
hourly_capacity_factor(y, t) = inputs["pP_Max"][y, t]
reg_max(y) = dfGen[y, :Reg_Max]
rsv_max(y) = dfGen[y, :Rsv_Max]

hourly_capacity(y, t) = hourly_capacity_factor(y, t) * eTotalCap[y]
resources_in_bin(y) = UnitRange(y, y + dfGen[y, :Num_VRE_Bins] - 1)
hourly_bin_capacity(y, t) = sum(hourly_capacity(yy, t) for yy in resources_in_bin(y))

@constraint(EP, [y in REG, t in 1:T], vREG[y, t] <= reg_max(y) * hourly_bin_capacity(y, t))
@constraint(EP, [y in RSV, t in 1:T], vRSV[y, t] <= rsv_max(y) * hourly_bin_capacity(y, t))

expr = extract_time_series_to_expression(vP, VRE_POWER_OUT)
add_similar_to_expression!(expr[REG, :], -vREG[REG, :])
@constraint(EP, [y in VRE_POWER_OUT, t in 1:T], expr[y, t] >= 0)

expr = extract_time_series_to_expression(vP, VRE_POWER_OUT)
add_similar_to_expression!(expr[REG, :], +vREG[REG, :])
add_similar_to_expression!(expr[RSV, :], +vRSV[RSV, :])
@constraint(EP, [y in VRE_POWER_OUT, t in 1:T], expr[y, t] <= hourly_bin_capacity(y, t))
end

function remove_reserves_for_binned_vre_resources!(EP::Model, inputs::Dict)
dfGen = inputs["dfGen"]

VRE = inputs["VRE"]
VRE_POWER_OUT = intersect(dfGen[dfGen.Num_VRE_Bins.>=1,:R_ID], VRE)
REG = inputs["REG"]
RSV = inputs["RSV"]

VRE_NO_POWER_OUT = setdiff(VRE, VRE_POWER_OUT)

for y in intersect(VRE_NO_POWER_OUT, REG)
fix.(EP[:vREG][y,:], 0.0, force=true)
end
for y in intersect(VRE_NO_POWER_OUT, RSV)
fix.(EP[:vRSV][y,:], 0.0, force=true)
end
end
55 changes: 21 additions & 34 deletions src/model/resources/hydro/hydro_res.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -185,45 +185,32 @@ function hydro_res_reserves!(EP::Model, inputs::Dict)
T = inputs["T"] # Number of time steps (hours)

HYDRO_RES = inputs["HYDRO_RES"]
REG = inputs["REG"]
RSV = inputs["RSV"]

HYDRO_RES_REG_RSV = intersect(HYDRO_RES, inputs["REG"], inputs["RSV"]) # Set of reservoir hydro resources with both regulation and spinning reserves
HYDRO_RES_REG = intersect(HYDRO_RES, REG) # Set of reservoir hydro resources with regulation reserves
HYDRO_RES_RSV = intersect(HYDRO_RES, RSV) # Set of reservoir hydro resources with spinning reserves

HYDRO_RES_REG = intersect(HYDRO_RES, inputs["REG"]) # Set of reservoir hydro resources with regulation reserves
HYDRO_RES_RSV = intersect(HYDRO_RES, inputs["RSV"]) # Set of reservoir hydro resources with spinning reserves
vP = EP[:vP]
vREG = EP[:vREG]
vRSV = EP[:vRSV]
eTotalCap = EP[:eTotalCap]
reg_max(y) = dfGen[y, :Reg_Max]
rsv_max(y) = dfGen[y, :Rsv_Max]

HYDRO_RES_REG_ONLY = setdiff(HYDRO_RES_REG, HYDRO_RES_RSV) # Set of reservoir hydro resources only with regulation reserves
HYDRO_RES_RSV_ONLY = setdiff(HYDRO_RES_RSV, HYDRO_RES_REG) # Set of reservoir hydro resources only with spinning reserves
max_up_reserves_lhs = extract_time_series_to_expression(vP, HYDRO_RES)
max_dn_reserves_lhs = extract_time_series_to_expression(vP, HYDRO_RES)

if !isempty(HYDRO_RES_REG_RSV)
@constraints(EP, begin
# Maximum storage contribution to reserves is a specified fraction of installed capacity
cRegulation[y in HYDRO_RES_REG_RSV, t in 1:T], EP[:vREG][y,t] <= dfGen[y,:Reg_Max]*EP[:eTotalCap][y]
cReserve[y in HYDRO_RES_REG_RSV, t in 1:T], EP[:vRSV][y,t] <= dfGen[y,:Rsv_Max]*EP[:eTotalCap][y]
# Maximum discharging rate and contribution to reserves up must be less than power rating
cMaxReservesUp[y in HYDRO_RES_REG_RSV, t in 1:T], EP[:vP][y,t]+EP[:vREG][y,t]+EP[:vRSV][y,t] <= EP[:eTotalCap][y]
# Maximum discharging rate and contribution to regulation down must be greater than zero
cMaxReservesDown[y in HYDRO_RES_REG_RSV, t in 1:T], EP[:vP][y,t]-EP[:vREG][y,t] >= 0
end)
end
S = HYDRO_RES_REG
add_similar_to_expression!(max_up_reserves_lhs[S, :], vREG[S, :])
add_similar_to_expression!(max_dn_reserves_lhs[S, :], -vREG[S, :])

if !isempty(HYDRO_RES_REG_ONLY)
@constraints(EP, begin
# Maximum storage contribution to reserves is a specified fraction of installed capacity
cRegulation[y in HYDRO_RES_REG_ONLY, t in 1:T], EP[:vREG][y,t] <= dfGen[y,:Reg_Max]*EP[:eTotalCap][y]
# Maximum discharging rate and contribution to reserves up must be less than power rating
cMaxReservesUp[y in HYDRO_RES_REG_ONLY, t in 1:T], EP[:vP][y,t]+EP[:vREG][y,t] <= EP[:eTotalCap][y]
# Maximum discharging rate and contribution to regulation down must be greater than zero
cMaxReservesDown[y in HYDRO_RES_REG_ONLY, t in 1:T], EP[:vP][y,t]-EP[:vREG][y,t] >= 0
end)
end
S = HYDRO_RES_RSV
add_similar_to_expression!(max_up_reserves_lhs[S, :], vRSV[S, :])

if !isempty(HYDRO_RES_RSV_ONLY)
@constraints(EP, begin
# Maximum storage contribution to reserves is a specified fraction of installed capacity
cReserve[y in HYDRO_RES_RSV_ONLY, t in 1:T], EP[:vRSV][y,t] <= dfGen[y,:Rsv_Max]*EP[:eTotalCap][y]
# Maximum discharging rate and contribution to reserves up must be less than power rating
cMaxReservesUp[y in HYDRO_RES_RSV_ONLY, t in 1:T], EP[:vP][y,t]+EP[:vRSV][y,t] <= EP[:eTotalCap][y]
end)
end
@constraint(EP, [y in HYDRO_RES, t in 1:T], max_up_reserves_lhs[y, t] <= eTotalCap[y])
@constraint(EP, [y in HYDRO_RES, t in 1:T], max_dn_reserves_lhs[y, t] >= 0)

@constraint(EP, [y in HYDRO_RES_REG, t in 1:T], vREG[y, t] <= reg_max(y) * eTotalCap[y])
@constraint(EP, [y in HYDRO_RES_RSV, t in 1:T], vRSV[y, t] <= rsv_max(y) * eTotalCap[y])
end
Loading

0 comments on commit 2ef54dc

Please sign in to comment.