Refactor reserves

src/model/expression_manipulation.jl

@@ -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
@@ -141,4 +168,4 @@ end
 
 function sum_expression(expr::AbstractArray{AbstractJuMPScalar, dims}) where {dims}
     return _sum_expression(expr)
-end
+end

src/model/resources/curtailable_variable_renewable/curtailable_variable_renewable.jl

@@ -47,24 +47,24 @@
 		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
@@ -102,60 +102,53 @@
 ```
 """
 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

src/model/resources/hydro/hydro_res.jl

@@ -185,45 +185,32 @@
 	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