|
| 1 | +# using Revise |
| 2 | +using JuMP |
| 3 | +using HiGHS |
| 4 | +# using Xpress |
| 5 | +using JSON |
| 6 | +using REopt |
| 7 | +using Logging |
| 8 | +using DotEnv |
| 9 | +# using PlotlyJS # Commented out for GitHub Actions |
| 10 | +DotEnv.load!() |
| 11 | + |
| 12 | +# ============================================================================ |
| 13 | +# SOLVER CONFIGURATION - Change solver here |
| 14 | +# ============================================================================ |
| 15 | +const USE_XPRESS = false # Set to false to use HiGHS |
| 16 | + |
| 17 | +""" |
| 18 | +Helper function to create model with configured solver |
| 19 | +""" |
| 20 | +function create_model(mip_gap::Float64=0.05; verbose::Bool=false) |
| 21 | + if USE_XPRESS |
| 22 | + m = Model(Xpress.Optimizer) |
| 23 | + set_optimizer_attribute(m, "OUTPUTLOG", verbose ? 1 : 0) |
| 24 | + set_optimizer_attribute(m, "MIPRELSTOP", mip_gap) |
| 25 | + else |
| 26 | + m = Model(HiGHS.Optimizer) |
| 27 | + set_optimizer_attribute(m, "output_flag", verbose) |
| 28 | + set_optimizer_attribute(m, "log_to_console", verbose) |
| 29 | + set_optimizer_attribute(m, "mip_rel_gap", mip_gap) |
| 30 | + end |
| 31 | + return m |
| 32 | +end |
| 33 | + |
| 34 | +println("\n" * "="^80) |
| 35 | +println("Testing Monte Carlo vs Discrete OUU Methods") |
| 36 | +println("="^80) |
| 37 | + |
| 38 | +# Suppress info/warning messages |
| 39 | +original_logger = Logging.global_logger() |
| 40 | +Logging.global_logger(Logging.SimpleLogger(stderr, Logging.Error)) |
| 41 | + |
| 42 | +# ============================================================================ |
| 43 | +# Test 0: Baseline (No Uncertainty) with BAU |
| 44 | +# ============================================================================ |
| 45 | +println("\n" * "─"^80) |
| 46 | +println("Test 0: Baseline (No Uncertainty)") |
| 47 | +println("─"^80) |
| 48 | +println("Single scenario with no uncertainty, comparing BAU vs technology optimization") |
| 49 | + |
| 50 | +scenario_baseline = JSON.parsefile("scenarios/ouu_base.json") |
| 51 | + |
| 52 | +print("Building scenario... ") |
| 53 | +s_baseline = Scenario(scenario_baseline) |
| 54 | +println("✓") |
| 55 | + |
| 56 | +print("Building REoptInputs... ") |
| 57 | +inputs_baseline = REoptInputs(s_baseline) |
| 58 | +println("✓") |
| 59 | +println(" └─ Number of scenarios: ", inputs_baseline.n_scenarios) |
| 60 | + |
| 61 | +# Create BAU and technology optimization models |
| 62 | +m_bau = create_model() |
| 63 | + |
| 64 | +m_tech = create_model() |
| 65 | + |
| 66 | +print("Building and solving models (BAU + Technology Optimization)... ") |
| 67 | +results_baseline = run_reopt([m_bau, m_tech], inputs_baseline) |
| 68 | +println("✓") |
| 69 | + |
| 70 | +if termination_status(m_bau) == MOI.OPTIMAL && termination_status(m_tech) == MOI.OPTIMAL |
| 71 | + println("\n📊 Baseline Results:") |
| 72 | + println(" BAU (Business as Usual):") |
| 73 | + println(" └─ Objective value: \$", round(objective_value(m_bau), digits=2)) |
| 74 | + println(" └─ Grid energy supplied: ", round(results_baseline["ElectricUtility"]["annual_energy_supplied_kwh_bau"], digits=1), " kWh") |
| 75 | + println("\n Technology Optimal:") |
| 76 | + println(" └─ Objective value: \$", round(objective_value(m_tech), digits=2)) |
| 77 | + println(" └─ Grid energy supplied: ", round(results_baseline["ElectricUtility"]["annual_energy_supplied_kwh"], digits=1), " kWh") |
| 78 | + println(" └─ PV size: ", round(results_baseline["PV"]["size_kw"], digits=1), " kW") |
| 79 | + println(" └─ Battery power: ", round(results_baseline["ElectricStorage"]["size_kw"], digits=1), " kW") |
| 80 | + println(" └─ Battery energy: ", round(results_baseline["ElectricStorage"]["size_kwh"], digits=1), " kWh") |
| 81 | + |
| 82 | + savings = objective_value(m_bau) - objective_value(m_tech) |
| 83 | + println("\n Technology Savings: \$", round(savings, digits=2), " (", round(savings/objective_value(m_bau)*100, digits=2), "%)") |
| 84 | +else |
| 85 | + println("\n⚠️ Baseline model status - BAU: ", termination_status(m_bau), ", Tech: ", termination_status(m_tech)) |
| 86 | +end |
| 87 | + |
| 88 | +# ============================================================================ |
| 89 | +# Test 1: Time-Invariant Method (Original Implementation) |
| 90 | +# ============================================================================ |
| 91 | +println("\n" * "─"^80) |
| 92 | +println("Test 1: Time-Invariant Method") |
| 93 | +println("─"^80) |
| 94 | +println("Creates 3 load scenarios × 3 PV scenarios = 9 total scenarios") |
| 95 | +println("Each scenario has uniform deviation across all timesteps") |
| 96 | + |
| 97 | +scenario_invariant = JSON.parsefile("scenarios/ouu_base.json") |
| 98 | +scenario_invariant["ElectricLoad"]["uncertainty"] = Dict( |
| 99 | + "enabled" => true, |
| 100 | + "method" => "time_invariant", |
| 101 | + "deviation_fractions" => [-0.1, 0.0, 0.1], |
| 102 | + "deviation_probabilities" => [0.25, 0.50, 0.25] |
| 103 | +) |
| 104 | +scenario_invariant["PV"]["production_uncertainty"] = Dict( |
| 105 | + "enabled" => true, |
| 106 | + "method" => "time_invariant", |
| 107 | + "deviation_fractions" => [-0.2, 0.0, 0.2], |
| 108 | + "deviation_probabilities" => [0.25, 0.50, 0.25] |
| 109 | +) |
| 110 | + |
| 111 | +print("Building scenario... ") |
| 112 | +s_invariant = Scenario(scenario_invariant) |
| 113 | +println("✓") |
| 114 | + |
| 115 | +print("Building REoptInputs... ") |
| 116 | +inputs_invariant = REoptInputs(s_invariant) |
| 117 | +println("✓") |
| 118 | +println(" └─ Number of scenarios: ", inputs_invariant.n_scenarios) |
| 119 | +println(" └─ Scenario probabilities: ", round.(inputs_invariant.scenario_probabilities[1:min(10, end)], digits=4)) |
| 120 | + |
| 121 | +m_invariant = create_model() |
| 122 | + |
| 123 | +m_invariant_bau = create_model() |
| 124 | + |
| 125 | +print("Building and solving models (BAU + Technology Optimization)... ") |
| 126 | +results_invariant = run_reopt([m_invariant_bau, m_invariant], inputs_invariant) |
| 127 | +println("✓") |
| 128 | + |
| 129 | +if termination_status(m_invariant) == MOI.OPTIMAL |
| 130 | + println("\n📊 Time-Invariant Method Results:") |
| 131 | + println(" Technology Optimal:") |
| 132 | + println(" └─ Objective value: \$", round(objective_value(m_invariant), digits=2)) |
| 133 | + println(" └─ Grid energy supplied: ", round(results_invariant["ElectricUtility"]["annual_energy_supplied_kwh"], digits=1), " kWh") |
| 134 | + println(" └─ PV size: ", round(results_invariant["PV"]["size_kw"], digits=1), " kW") |
| 135 | + println(" └─ Battery power: ", round(results_invariant["ElectricStorage"]["size_kw"], digits=1), " kW") |
| 136 | + println(" └─ Battery energy: ", round(results_invariant["ElectricStorage"]["size_kwh"], digits=1), " kWh") |
| 137 | + if termination_status(m_invariant_bau) == MOI.OPTIMAL |
| 138 | + println("\n BAU:") |
| 139 | + println(" └─ Objective value: \$", round(objective_value(m_invariant_bau), digits=2)) |
| 140 | + println(" └─ Grid energy supplied: ", round(results_invariant["ElectricUtility"]["annual_energy_supplied_kwh_bau"], digits=1), " kWh") |
| 141 | + savings = objective_value(m_invariant_bau) - objective_value(m_invariant) |
| 142 | + println(" └─ Technology Savings: \$", round(savings, digits=2), " (", round(savings/objective_value(m_invariant_bau)*100, digits=2), "%)") |
| 143 | + end |
| 144 | +else |
| 145 | + println("\n⚠️ Time-Invariant model status: ", termination_status(m_invariant)) |
| 146 | +end |
| 147 | + |
| 148 | +# ============================================================================ |
| 149 | +# Test 2: Monte Carlo Method (New Implementation) |
| 150 | +# ============================================================================ |
| 151 | +println("\n" * "─"^80) |
| 152 | +println("Test 2: Monte Carlo Method") |
| 153 | +println("─"^80) |
| 154 | +println("Creates 9 scenarios (3 load samples × 3 PV samples)") |
| 155 | +println("Each scenario has timestep-varying deviations sampled from distribution") |
| 156 | + |
| 157 | +scenario_mc = JSON.parsefile("scenarios/ouu_base.json") |
| 158 | +scenario_mc["ElectricLoad"]["uncertainty"] = Dict( |
| 159 | + "enabled" => true, |
| 160 | + "method" => "discrete", |
| 161 | + "deviation_fractions" => [-0.1, 0.0, 0.1], |
| 162 | + "deviation_probabilities" => [0.25, 0.50, 0.25], |
| 163 | + "n_samples" => 3 # Each sample has different deviation per timestep |
| 164 | +) |
| 165 | +scenario_mc["PV"]["production_uncertainty"] = Dict( |
| 166 | + "enabled" => true, |
| 167 | + "method" => "discrete", |
| 168 | + "deviation_fractions" => [-0.2, 0.0, 0.2], |
| 169 | + "deviation_probabilities" => [0.25, 0.50, 0.25], |
| 170 | + "n_samples" => 3 |
| 171 | +) |
| 172 | + |
| 173 | +print("Building scenario... ") |
| 174 | +s_mc = Scenario(scenario_mc) |
| 175 | +println("✓") |
| 176 | + |
| 177 | +print("Building REoptInputs... ") |
| 178 | +inputs_mc = REoptInputs(s_mc) |
| 179 | +println("✓") |
| 180 | +println(" └─ Number of scenarios: ", inputs_mc.n_scenarios) |
| 181 | +println(" └─ Scenario probabilities (first 5): ", round.(inputs_mc.scenario_probabilities[1:min(5, end)], digits=4)) |
| 182 | +println(" └─ All equal? ", all(x -> isapprox(x, inputs_mc.scenario_probabilities[1], atol=1e-10), inputs_mc.scenario_probabilities)) |
| 183 | + |
| 184 | +m_mc = create_model() |
| 185 | + |
| 186 | +m_mc_bau = create_model() |
| 187 | + |
| 188 | +print("Building and solving models (BAU + Technology Optimization)... ") |
| 189 | +results_mc = run_reopt([m_mc_bau, m_mc], inputs_mc) |
| 190 | +println("✓") |
| 191 | + |
| 192 | +if termination_status(m_mc) == MOI.OPTIMAL |
| 193 | + println("\n📊 Monte Carlo Method Results:") |
| 194 | + println(" Technology Optimal:") |
| 195 | + println(" └─ Objective value: \$", round(objective_value(m_mc), digits=2)) |
| 196 | + println(" └─ Grid energy supplied: ", round(results_mc["ElectricUtility"]["annual_energy_supplied_kwh"], digits=1), " kWh") |
| 197 | + println(" └─ PV size: ", round(results_mc["PV"]["size_kw"], digits=1), " kW") |
| 198 | + println(" └─ Battery power: ", round(results_mc["ElectricStorage"]["size_kw"], digits=1), " kW") |
| 199 | + println(" └─ Battery energy: ", round(results_mc["ElectricStorage"]["size_kwh"], digits=1), " kWh") |
| 200 | + if termination_status(m_mc_bau) == MOI.OPTIMAL |
| 201 | + println("\n BAU:") |
| 202 | + println(" └─ Objective value: \$", round(objective_value(m_mc_bau), digits=2)) |
| 203 | + println(" └─ Grid energy supplied: ", round(results_mc["ElectricUtility"]["annual_energy_supplied_kwh_bau"], digits=1), " kWh") |
| 204 | + savings = objective_value(m_mc_bau) - objective_value(m_mc) |
| 205 | + println(" └─ Technology Savings: \$", round(savings, digits=2), " (", round(savings/objective_value(m_mc_bau)*100, digits=2), "%)") |
| 206 | + end |
| 207 | +else |
| 208 | + println("\n⚠️ Monte Carlo model status: ", termination_status(m_mc)) |
| 209 | +end |
| 210 | + |
| 211 | +# ============================================================================ |
| 212 | +# Comparison |
| 213 | +# ============================================================================ |
| 214 | +if termination_status(m_invariant) == MOI.OPTIMAL && termination_status(m_mc) == MOI.OPTIMAL |
| 215 | + println("\n" * "─"^80) |
| 216 | + println("Comparison:") |
| 217 | + println("─"^80) |
| 218 | + |
| 219 | + pv_diff = results_mc["PV"]["size_kw"] - results_invariant["PV"]["size_kw"] |
| 220 | + batt_kw_diff = results_mc["ElectricStorage"]["size_kw"] - results_invariant["ElectricStorage"]["size_kw"] |
| 221 | + batt_kwh_diff = results_mc["ElectricStorage"]["size_kwh"] - results_invariant["ElectricStorage"]["size_kwh"] |
| 222 | + cost_diff = objective_value(m_mc) - objective_value(m_invariant) |
| 223 | + |
| 224 | + println(" Δ PV size: ", round(pv_diff, digits=1), " kW (", round(pv_diff/results_invariant["PV"]["size_kw"]*100, digits=1), "%)") |
| 225 | + println(" Δ Battery power: ", round(batt_kw_diff, digits=1), " kW (", round(batt_kw_diff/results_invariant["ElectricStorage"]["size_kw"]*100, digits=1), "%)") |
| 226 | + println(" Δ Battery energy: ", round(batt_kwh_diff, digits=1), " kWh (", round(batt_kwh_diff/results_invariant["ElectricStorage"]["size_kwh"]*100, digits=1), "%)") |
| 227 | + println(" Δ Cost: \$", round(cost_diff, digits=2), " (", round(cost_diff/objective_value(m_invariant)*100, digits=2), "%)") |
| 228 | + |
| 229 | + println("\n💡 Key Insight:") |
| 230 | + println(" Monte Carlo captures timestep-level uncertainty, while discrete assumes") |
| 231 | + println(" all timesteps move together. This can lead to different optimal sizing.") |
| 232 | + |
| 233 | + # ============================================================================ |
| 234 | + # Visualization: Compare scenario profiles |
| 235 | + # ============================================================================ |
| 236 | + println("\n" * "─"^80) |
| 237 | + println("Plotting scenario profiles...") |
| 238 | + println("─"^80) |
| 239 | + |
| 240 | + # Plot first 168 hours (1 week) for visibility |
| 241 | + # Plotting code commented out for GitHub Actions (PlotlyJS not available) |
| 242 | + # plot_hours = 1:min(8760, length(inputs_invariant.time_steps)) |
| 243 | + # |
| 244 | + # # Find baseline scenario (zero deviation) in time_invariant - should be scenario 2 with deviation=0.0 |
| 245 | + # baseline_scenario_id = 5 # Middle scenario with 0.0 deviation |
| 246 | + # mc_scenario_id = 2 # First Monte Carlo scenario |
| 247 | + # |
| 248 | + # # Create load comparison plot |
| 249 | + # load_trace1 = PlotlyJS.scatter( |
| 250 | + # x=collect(plot_hours), |
| 251 | + # y=inputs_discrete.loads_kw_by_scenario[baseline_scenario_id][plot_hours], |
| 252 | + # mode="lines", |
| 253 | + # name="Baseline (No Deviation)", |
| 254 | + # line=attr(width=2, color="blue") |
| 255 | + # ) |
| 256 | + # load_trace2 = PlotlyJS.scatter( |
| 257 | + # x=collect(plot_hours), |
| 258 | + # y=inputs_mc.loads_kw_by_scenario[mc_scenario_id][plot_hours], |
| 259 | + # mode="lines", |
| 260 | + # name="Discrete Monte Carlo Sample", |
| 261 | + # line=attr(width=2, color="red", dash="dash") |
| 262 | + # ) |
| 263 | + # |
| 264 | + # load_layout = PlotlyJS.Layout( |
| 265 | + # title="Load Profile Comparison (First Week)", |
| 266 | + # xaxis_title="Hour", |
| 267 | + # yaxis_title="Load (kW)", |
| 268 | + # hovermode="x unified" |
| 269 | + # ) |
| 270 | + # |
| 271 | + # load_plot = PlotlyJS.plot([load_trace1, load_trace2], load_layout) |
| 272 | + # |
| 273 | + # # Create PV production comparison plot |
| 274 | + # pv_trace1 = PlotlyJS.scatter( |
| 275 | + # x=collect(plot_hours), |
| 276 | + # y=inputs_discrete.production_factor_by_scenario[baseline_scenario_id]["PV"][plot_hours], |
| 277 | + # mode="lines", |
| 278 | + # name="Baseline (No Deviation)", |
| 279 | + # line=attr(width=2, color="blue") |
| 280 | + # ) |
| 281 | + # pv_trace2 = PlotlyJS.scatter( |
| 282 | + # x=collect(plot_hours), |
| 283 | + # y=inputs_mc.production_factor_by_scenario[mc_scenario_id]["PV"][plot_hours], |
| 284 | + # mode="lines", |
| 285 | + # name="Discrete Monte Carlo Sample", |
| 286 | + # line=attr(width=2, color="red", dash="dash") |
| 287 | + # ) |
| 288 | + # |
| 289 | + # pv_layout = PlotlyJS.Layout( |
| 290 | + # title="PV Production Profile Comparison (First Week)", |
| 291 | + # xaxis_title="Hour", |
| 292 | + # yaxis_title="Production Factor", |
| 293 | + # hovermode="x unified" |
| 294 | + # ) |
| 295 | + # |
| 296 | + # pv_plot = PlotlyJS.plot([pv_trace1, pv_trace2], pv_layout) |
| 297 | + # |
| 298 | + # # Save plots |
| 299 | + # PlotlyJS.savefig(load_plot, "load_comparison.html") |
| 300 | + # PlotlyJS.savefig(pv_plot, "pv_comparison.html") |
| 301 | + # println(" └─ Plots saved to load_comparison.html and pv_comparison.html") |
| 302 | +end |
| 303 | + |
| 304 | +# Restore logger |
| 305 | +Logging.global_logger(original_logger) |
| 306 | + |
| 307 | +println("\n" * "="^80) |
| 308 | +println("✅ Tests Complete") |
| 309 | +println("="^80) |
0 commit comments