Optimization Under Uncertainty: Load and PV/Wind

src/REopt.jl

@@ -145,6 +145,7 @@ include("core/electric_heater.jl")
 include("core/cst_ssc.jl")
 include("core/cst.jl")
 include("core/ashp.jl")
+include("core/uncertainty.jl")
 include("core/scenario.jl")
 include("core/bau_scenario.jl")
 include("core/reopt_inputs.jl")

src/constraints/battery_degradation.jl

@@ -29,7 +29,7 @@ function constrain_degradation_variables(m, p; b="ElectricStorage")
     end
 
     @constraint(m, [d in days],
-        m[:Eavg][d] == 1/ts_per_day * sum(m[:dvStoredEnergy][b, ts] for ts in ts0[d]:tsF[d])
+        m[:Eavg][d] == 1/ts_per_day * sum(sum(p.scenario_probabilities[s] * m[:dvStoredEnergy][s, b, ts] for s in 1:p.n_scenarios) for ts in ts0[d]:tsF[d])
     )
 
     @constraint(m, [d in days],
@@ -38,13 +38,13 @@ function constrain_degradation_variables(m, p; b="ElectricStorage")
 
     # Power in equals power into storage from grid or local production
     @constraint(m, [ts in p.time_steps],
-    sum(m[:dvSegmentChargePower][ts, j] for j in 1:J) == sum(
-        m[:dvProductionToStorage][b, t, ts] for t in p.techs.elec) + m[:dvGridToStorage][b, ts]
+    sum(m[:dvSegmentChargePower][ts, j] for j in 1:J) == sum(p.scenario_probabilities[s] * (sum(
+        m[:dvProductionToStorage][s, b, t, ts] for t in p.techs.elec) + m[:dvGridToStorage][s, b, ts]) for s in 1:p.n_scenarios)
     )
 
     # Power out equals power discharged from storage to any destination
     @constraint(m, [ts in p.time_steps],
-    sum(m[:dvSegmentDischargePower][ts, j] for j in 1:J) == m[:dvDischargeFromStorage][b, ts]);
+    sum(m[:dvSegmentDischargePower][ts, j] for j in 1:J) == sum(p.scenario_probabilities[s] * m[:dvDischargeFromStorage][s, b, ts] for s in 1:p.n_scenarios));
 
     # Balance charging with daily e_plus, here is only collect all power across the day, so don't need to times efficiency
     @constraint(m, [d in days, j in 1:J], m[:Eplus_sum][d, j] == sum(m[:dvSegmentChargePower][ts0[d]:tsF[d], j])*p.hours_per_time_step)

src/constraints/chp_constraints.jl

@@ -9,32 +9,33 @@ function add_chp_fuel_burn_constraints(m, p; _n="")
 
     # Fuel cost
     m[:TotalCHPFuelCosts] = @expression(m, 
-        sum(p.pwf_fuel[t] * m[:dvFuelUsage][t, ts] * p.fuel_cost_per_kwh[t][ts] for t in p.techs.chp, ts in p.time_steps)
+        sum(p.scenario_probabilities[s] * p.pwf_fuel[t] * m[:dvFuelUsage][s, t, ts] * p.fuel_cost_per_kwh[t][ts] 
+            for s in 1:p.n_scenarios, t in p.techs.chp, ts in p.time_steps)
     )      
     # Conditionally add dvFuelBurnYIntercept if coefficient p.FuelBurnYIntRate is greater than ~zero
     if abs(fuel_burn_intercept) > 1.0E-7
         dv = "dvFuelBurnYIntercept"*_n
         m[Symbol(dv)] = @variable(m, [p.techs.chp, p.time_steps], base_name=dv)
 
         #Constraint (1c1): Total Fuel burn for CHP **with** y-intercept fuel burn and supplementary firing
-        @constraint(m, CHPFuelBurnCon[t in p.techs.chp, ts in p.time_steps],
-            m[Symbol("dvFuelUsage"*_n)][t,ts]  == p.hours_per_time_step * (
+        @constraint(m, CHPFuelBurnCon[s in 1:p.n_scenarios, t in p.techs.chp, ts in p.time_steps],
+            m[Symbol("dvFuelUsage"*_n)][s,t,ts]  == p.hours_per_time_step * (
                 m[Symbol("dvFuelBurnYIntercept"*_n)][t,ts] +
-                p.production_factor[t,ts] * fuel_burn_slope * m[Symbol("dvRatedProduction"*_n)][t,ts] +
-                m[Symbol("dvSupplementaryThermalProduction"*_n)][t,ts] / p.s.chp.supplementary_firing_efficiency
+                p.production_factor_by_scenario[s][t][ts] * fuel_burn_slope * m[Symbol("dvRatedProduction"*_n)][s,t,ts] +
+                m[Symbol("dvSupplementaryThermalProduction"*_n)][s,t,ts] / p.s.chp.supplementary_firing_efficiency
             )
         )
         #Constraint (1d): Y-intercept fuel burn for CHP
-        @constraint(m, CHPFuelBurnYIntCon[t in p.techs.chp, ts in p.time_steps],
+        @constraint(m, CHPFuelBurnYIntCon[s in 1:p.n_scenarios, t in p.techs.chp, ts in p.time_steps],
                     fuel_burn_intercept * m[Symbol("dvSize"*_n)][t] - p.s.chp.max_kw * 
-                    (1-m[Symbol("binCHPIsOnInTS"*_n)][t,ts])  <= m[Symbol("dvFuelBurnYIntercept"*_n)][t,ts]
+                    (1-m[Symbol("binCHPIsOnInTS"*_n)][s,t,ts])  <= m[Symbol("dvFuelBurnYIntercept"*_n)][t,ts]
                     )
     else
         #Constraint (1c2): Total Fuel burn for CHP **without** y-intercept fuel burn
-        @constraint(m, CHPFuelBurnConLinear[t in p.techs.chp, ts in p.time_steps],
-            m[Symbol("dvFuelUsage"*_n)][t,ts]  == p.hours_per_time_step * (
-                p.production_factor[t,ts] * fuel_burn_slope * m[Symbol("dvRatedProduction"*_n)][t,ts] +
-                m[Symbol("dvSupplementaryThermalProduction"*_n)][t,ts] / p.s.chp.supplementary_firing_efficiency
+        @constraint(m, CHPFuelBurnConLinear[s in 1:p.n_scenarios, t in p.techs.chp, ts in p.time_steps],
+            m[Symbol("dvFuelUsage"*_n)][s,t,ts]  == p.hours_per_time_step * (
+                p.production_factor_by_scenario[s][t][ts] * fuel_burn_slope * m[Symbol("dvRatedProduction"*_n)][s,t,ts] +
+                m[Symbol("dvSupplementaryThermalProduction"*_n)][s,t,ts] / p.s.chp.supplementary_firing_efficiency
             )
         )
     end 
@@ -58,28 +59,28 @@ function add_chp_thermal_production_constraints(m, p; _n="")
             m[Symbol("dvHeatingProductionYIntercept"*_n)][t,ts] <= thermal_prod_intercept * m[Symbol("dvSize"*_n)][t]
         )
         # Constraint (2a-2): Upper Bounds on Thermal Production Y-Intercept
-        @constraint(m, CHPYInt2a2Con[t in p.techs.chp, ts in p.time_steps],
+        @constraint(m, CHPYInt2a2Con[s in 1:p.n_scenarios, t in p.techs.chp, ts in p.time_steps],
             m[Symbol("dvHeatingProductionYIntercept"*_n)][t,ts] <= thermal_prod_intercept * p.s.chp.max_kw 
-            * m[Symbol("binCHPIsOnInTS"*_n)][t,ts]
+            * m[Symbol("binCHPIsOnInTS"*_n)][s,t,ts]
         )
         #Constraint (2b): Lower Bounds on Thermal Production Y-Intercept
-        @constraint(m, CHPYInt2bCon[t in p.techs.chp, ts in p.time_steps],
+        @constraint(m, CHPYInt2bCon[s in 1:p.n_scenarios, t in p.techs.chp, ts in p.time_steps],
             m[Symbol("dvHeatingProductionYIntercept"*_n)][t,ts] >= thermal_prod_intercept * m[Symbol("dvSize"*_n)][t] 
-            - thermal_prod_intercept * p.s.chp.max_kw * (1 - m[Symbol("binCHPIsOnInTS"*_n)][t,ts])
+            - thermal_prod_intercept * p.s.chp.max_kw * (1 - m[Symbol("binCHPIsOnInTS"*_n)][s,t,ts])
         )
         # Constraint (2c): Thermal Production of CHP
         # Note: p.HotWaterAmbientFactor[t,ts] * p.HotWaterThermalFactor[t,ts] removed from this but present in math
-        @constraint(m, CHPThermalProductionCon[t in p.techs.chp, ts in p.time_steps],
-        sum(m[Symbol("dvHeatingProduction"*_n)][t,q,ts] for q in p.heating_loads) ==
-            thermal_prod_slope * p.production_factor[t,ts] * m[Symbol("dvRatedProduction"*_n)][t,ts] 
+        @constraint(m, CHPThermalProductionCon[s in 1:p.n_scenarios, t in p.techs.chp, ts in p.time_steps],
+        sum(m[Symbol("dvHeatingProduction"*_n)][s,t,q,ts] for q in p.heating_loads) ==
+            thermal_prod_slope * p.production_factor_by_scenario[s][t][ts] * m[Symbol("dvRatedProduction"*_n)][s,t,ts] 
             + m[Symbol("dvHeatingProductionYIntercept"*_n)][t,ts] +
-            m[Symbol("dvSupplementaryThermalProduction"*_n)][t,ts]
+            m[Symbol("dvSupplementaryThermalProduction"*_n)][s,t,ts]
         )
     else
-        @constraint(m, CHPThermalProductionConLinear[t in p.techs.chp, ts in p.time_steps],
-            sum(m[Symbol("dvHeatingProduction"*_n)][t,q,ts] for q in p.heating_loads) ==
-            thermal_prod_slope * p.production_factor[t,ts] * m[Symbol("dvRatedProduction"*_n)][t,ts] +
-            m[Symbol("dvSupplementaryThermalProduction"*_n)][t,ts]
+        @constraint(m, CHPThermalProductionConLinear[s in 1:p.n_scenarios, t in p.techs.chp, ts in p.time_steps],
+            sum(m[Symbol("dvHeatingProduction"*_n)][s,t,q,ts] for q in p.heating_loads) ==
+            thermal_prod_slope * p.production_factor_by_scenario[s][t][ts] * m[Symbol("dvRatedProduction"*_n)][s,t,ts] +
+            m[Symbol("dvSupplementaryThermalProduction"*_n)][s,t,ts]
         )        
     end
 
@@ -98,39 +99,39 @@ function add_chp_supplementary_firing_constraints(m, p; _n="")
     thermal_prod_slope = (thermal_prod_full_load - thermal_prod_half_load) / (1.0 - 0.5)  # [kWt/kWe]
 
     # Constrain upper limit of dvSupplementaryThermalProduction, using auxiliary variable for (size * useSupplementaryFiring)
-    @constraint(m, CHPSupplementaryFireCon[t in p.techs.chp, ts in p.time_steps],
-                m[Symbol("dvSupplementaryThermalProduction"*_n)][t,ts] <=
-                (p.s.chp.supplementary_firing_max_steam_ratio - 1.0) * p.production_factor[t,ts] * (thermal_prod_slope * m[Symbol("dvSupplementaryFiringSize"*_n)][t] + m[Symbol("dvHeatingProductionYIntercept"*_n)][t,ts])
+    @constraint(m, CHPSupplementaryFireCon[s in 1:p.n_scenarios, t in p.techs.chp, ts in p.time_steps],
+                m[Symbol("dvSupplementaryThermalProduction"*_n)][s,t,ts] <=
+                (p.s.chp.supplementary_firing_max_steam_ratio - 1.0) * p.production_factor_by_scenario[s][t][ts] * (thermal_prod_slope * m[Symbol("dvSupplementaryFiringSize"*_n)][t] + m[Symbol("dvHeatingProductionYIntercept"*_n)][t,ts])
                 )
     if solver_is_compatible_with_indicator_constraints(p.s.settings.solver_name)
         # Constrain lower limit of 0 if CHP tech is off
-        @constraint(m, NoCHPSupplementaryFireOffCon[t in p.techs.chp, ts in p.time_steps],
-                !m[Symbol("binCHPIsOnInTS"*_n)][t,ts] => {m[Symbol("dvSupplementaryThermalProduction"*_n)][t,ts] <= 0.0}
+        @constraint(m, NoCHPSupplementaryFireOffCon[s in 1:p.n_scenarios, t in p.techs.chp, ts in p.time_steps],
+                !m[Symbol("binCHPIsOnInTS"*_n)][s,t,ts] => {m[Symbol("dvSupplementaryThermalProduction"*_n)][s,t,ts] <= 0.0}
                 )
     else
         #There's no upper bound specified for the CHP supplementary firing, so assume the entire heat load as a reasonable maximum that wouldn't be exceeded (but might not be the best possible value). 
         max_supplementary_firing_size = maximum(p.s.dhw_load.loads_kw .+ p.s.space_heating_load.loads_kw)
-        @constraint(m, NoCHPSupplementaryFireOffCon[t in p.techs.chp, ts in p.time_steps],
-                m[Symbol("dvSupplementaryThermalProduction"*_n)][t,ts] <= (p.s.chp.supplementary_firing_max_steam_ratio - 1.0) * p.production_factor[t,ts] * (thermal_prod_slope * max_supplementary_firing_size + m[Symbol("dvHeatingProductionYIntercept"*_n)][t,ts])
+        @constraint(m, NoCHPSupplementaryFireOffCon[s in 1:p.n_scenarios, t in p.techs.chp, ts in p.time_steps],
+                m[Symbol("dvSupplementaryThermalProduction"*_n)][s,t,ts] <= (p.s.chp.supplementary_firing_max_steam_ratio - 1.0) * p.production_factor_by_scenario[s][t][ts] * (thermal_prod_slope * max_supplementary_firing_size + m[Symbol("dvHeatingProductionYIntercept"*_n)][t,ts])
                 )
     end
 end
 
 function add_binCHPIsOnInTS_constraints(m, p; _n="")
     # Note, min_turn_down_fraction for CHP is only enforced in p.time_steps_with_grid
-    @constraint(m, [t in p.techs.chp, ts in p.time_steps_with_grid],
-        m[Symbol("dvRatedProduction"*_n)][t, ts] <= p.s.chp.max_kw * m[Symbol("binCHPIsOnInTS"*_n)][t, ts]
+    @constraint(m, [s in 1:p.n_scenarios, t in p.techs.chp, ts in p.time_steps_with_grid],
+        m[Symbol("dvRatedProduction"*_n)][s, t, ts] <= p.s.chp.max_kw * m[Symbol("binCHPIsOnInTS"*_n)][s, t, ts]
     )
-    @constraint(m, [t in p.techs.chp, ts in p.time_steps_with_grid],
-        p.s.chp.min_turn_down_fraction * m[Symbol("dvSize"*_n)][t] - m[Symbol("dvRatedProduction"*_n)][t, ts] <=
-        p.s.chp.max_kw * (1 - m[Symbol("binCHPIsOnInTS"*_n)][t, ts])
+    @constraint(m, [s in 1:p.n_scenarios, t in p.techs.chp, ts in p.time_steps_with_grid],
+        p.s.chp.min_turn_down_fraction * m[Symbol("dvSize"*_n)][t] - m[Symbol("dvRatedProduction"*_n)][s, t, ts] <=
+        p.s.chp.max_kw * (1 - m[Symbol("binCHPIsOnInTS"*_n)][s, t, ts])
     )
 end
 
 
 function add_chp_rated_prod_constraint(m, p; _n="")
-    @constraint(m, [t in p.techs.chp, ts in p.time_steps],
-        m[Symbol("dvSize"*_n)][t] >= m[Symbol("dvRatedProduction"*_n)][t, ts]
+    @constraint(m, [s in 1:p.n_scenarios, t in p.techs.chp, ts in p.time_steps],
+        m[Symbol("dvSize"*_n)][t] >= m[Symbol("dvRatedProduction"*_n)][s, t, ts]
     )
 end
 
@@ -146,9 +147,9 @@ function add_chp_hourly_om_charges(m, p; _n="")
     m[Symbol(dv)] = @variable(m, [p.techs.chp, p.time_steps], base_name=dv, lower_bound=0)
 
     #Constraint CHP-hourly-om-a: om per hour, per time step >= per_unit_size_cost * size for when on, >= zero when off
-	@constraint(m, CHPHourlyOMBySizeA[t in p.techs.chp, ts in p.time_steps],
+	@constraint(m, CHPHourlyOMBySizeA[s in 1:p.n_scenarios, t in p.techs.chp, ts in p.time_steps],
         p.s.chp.om_cost_per_hr_per_kw_rated * m[Symbol("dvSize"*_n)][t] -
-        p.s.chp.max_kw * p.s.chp.om_cost_per_hr_per_kw_rated * (1-m[Symbol("binCHPIsOnInTS"*_n)][t,ts])
+        p.s.chp.max_kw * p.s.chp.om_cost_per_hr_per_kw_rated * (1-m[Symbol("binCHPIsOnInTS"*_n)][s,t,ts])
             <= m[Symbol("dvOMByHourBySizeCHP"*_n)][t, ts]
     )
 	#Constraint CHP-hourly-om-b: om per hour, per time step <= per_unit_size_cost * size for each hour
@@ -157,8 +158,8 @@ function add_chp_hourly_om_charges(m, p; _n="")
             >= m[Symbol("dvOMByHourBySizeCHP"*_n)][t, ts]
     )
 	#Constraint CHP-hourly-om-c: om per hour, per time step <= zero when off, <= per_unit_size_cost*max_size
-	@constraint(m, CHPHourlyOMBySizeC[t in p.techs.chp, ts in p.time_steps],
-        p.s.chp.max_kw * p.s.chp.om_cost_per_hr_per_kw_rated * m[Symbol("binCHPIsOnInTS"*_n)][t,ts]
+	@constraint(m, CHPHourlyOMBySizeC[s in 1:p.n_scenarios, t in p.techs.chp, ts in p.time_steps],
+        p.s.chp.max_kw * p.s.chp.om_cost_per_hr_per_kw_rated * m[Symbol("binCHPIsOnInTS"*_n)][s,t,ts]
             >= m[Symbol("dvOMByHourBySizeCHP"*_n)][t, ts]
     )
 
@@ -179,14 +180,14 @@ function add_chp_constraints(m, p; _n="")
     @warn """Adding binary variable to model CHP. 
                 Some solvers are very slow with integer variables"""
     @variables m begin
-        binCHPIsOnInTS[p.techs.chp, p.time_steps], Bin  # 1 If technology t is operating in time step; 0 otherwise
+        binCHPIsOnInTS[1:p.n_scenarios, p.techs.chp, p.time_steps], Bin  # 1 If technology t is operating in time step in scenario s; 0 otherwise
     end    
 
     m[:TotalHourlyCHPOMCosts] = 0
     m[:TotalCHPFuelCosts] = 0
     m[:TotalCHPPerUnitProdOMCosts] = @expression(m, p.third_party_factor * p.pwf_om *
-        sum(p.s.chp.om_cost_per_kwh * p.hours_per_time_step *
-        m[:dvRatedProduction][t, ts] for t in p.techs.chp, ts in p.time_steps)
+        sum(p.scenario_probabilities[s] * p.s.chp.om_cost_per_kwh * p.hours_per_time_step *
+        m[:dvRatedProduction][s, t, ts] for s in 1:p.n_scenarios, t in p.techs.chp, ts in p.time_steps)
     )
 
     if p.s.chp.om_cost_per_hr_per_kw_rated > 1.0E-7
@@ -203,8 +204,8 @@ function add_chp_constraints(m, p; _n="")
     else
         for t in p.techs.chp
             fix(m[Symbol("dvSupplementaryFiringSize"*_n)][t], 0.0, force=true)
-            for ts in p.time_steps
-                fix(m[Symbol("dvSupplementaryThermalProduction"*_n)][t,ts], 0.0, force=true)
+            for s in 1:p.n_scenarios, ts in p.time_steps
+                fix(m[Symbol("dvSupplementaryThermalProduction"*_n)][s,t,ts], 0.0, force=true)
             end
         end
     end