Generation expansion
This tutorial was generated using Literate.jl. Download the source as a .jl file. Download the source as a .ipynb file.
using SDDPimport HiGHSimport Testfunction generation_expansion(duality_handler) build_cost = 1e4 use_cost = 4 num_units = 5 capacities = ones(num_units) demand_vals = 0.5 * [ 5 5 5 5 5 5 5 5 4 3 1 3 0 9 8 17 0 9 4 2 19 19 13 7 25 11 4 14 4 6 15 12 6 7 5 3 8 4 17 13 ] # Cost of unmet demand penalty = 5e5 # Discounting rate rho = 0.99 model = SDDP.LinearPolicyGraph(; stages = 5, lower_bound = 0.0, optimizer = HiGHS.Optimizer, ) do sp, stage @variable( sp, 0 <= invested[1:num_units] <= 1, SDDP.State, Int, initial_value = 0 ) @variables(sp, begin generation >= 0 unmet >= 0 demand end) @constraints( sp, begin # Can't un-invest investment[i in 1:num_units], invested[i].out >= invested[i].in # Generation capacity sum(capacities[i] * invested[i].out for i in 1:num_units) >= generation # Meet demand or pay a penalty unmet >= demand - sum(generation) # For fewer iterations order the units to break symmetry, units are identical (tougher numerically) [j in 1:(num_units-1)], invested[j].out <= invested[j+1].out end ) # Demand is uncertain SDDP.parameterize(ω -> fix(demand, ω), sp, demand_vals[stage, :]) @expression( sp, investment_cost, build_cost * sum(invested[i].out - invested[i].in for i in 1:num_units) ) @stageobjective( sp, (investment_cost + generation * use_cost) * rho^(stage - 1) + penalty * unmet ) end if get(ARGS, 1, "") == "--write" # Run `$ julia generation_expansion.jl --write` to update the benchmark # model directory model_dir = joinpath(@__DIR__, "..", "..", "..", "benchmarks", "models") SDDP.write_to_file( model, joinpath(model_dir, "generation_expansion.sof.json.gz"); test_scenarios = 100, ) exit(0) end SDDP.train(model; log_frequency = 10, duality_handler = duality_handler) Test.@test SDDP.calculate_bound(model) ≈ 2.078860e6 atol = 1e3 returnendgeneration_expansion(SDDP.ContinuousConicDuality())-------------------------------------------------------------------
SDDP.jl (c) Oscar Dowson and contributors, 2017-26
-------------------------------------------------------------------
problem
nodes : 5
state variables : 5
scenarios : 3.27680e+04
existing cuts : false
options
solver : serial mode
risk measure : SDDP.Expectation()
sampling scheme : SDDP.InSampleMonteCarlo
subproblem structure
VariableRef : [14, 14]
AffExpr in MOI.GreaterThan{Float64} : [7, 7]
AffExpr in MOI.LessThan{Float64} : [4, 4]
VariableRef in MOI.GreaterThan{Float64} : [8, 8]
VariableRef in MOI.Integer : [5, 5]
VariableRef in MOI.LessThan{Float64} : [5, 6]
numerical stability report
matrix range [1e+00, 1e+00]
objective range [1e+00, 5e+05]
bounds range [1e+00, 1e+00]
rhs range [0e+00, 0e+00]
-------------------------------------------------------------------
iteration simulation bound time (s) solves pid
-------------------------------------------------------------------
10 2.549668e+06 2.075239e+06 1.899990e+00 920 1
20 5.494568e+05 2.078257e+06 2.351343e+00 1340 1
30 4.985879e+04 2.078257e+06 4.134070e+00 2260 1
40 3.799447e+06 2.078257e+06 4.592905e+00 2680 1
41 1.299866e+06 2.078257e+06 4.632773e+00 2722 1
-------------------------------------------------------------------
status : simulation_stopping
total time (s) : 4.632773e+00
total solves : 2722
best bound : 2.078257e+06
numeric issues : 0
-------------------------------------------------------------------