Generation expansion
This tutorial was generated using Literate.jl. Download the source as a .jl file. Download the source as a .ipynb file.
using SDDP
import HiGHS
import Test
function 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(ω -> JuMP.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
return
end
generation_expansion(SDDP.ContinuousConicDuality())
generation_expansion(SDDP.LagrangianDuality())-------------------------------------------------------------------
SDDP.jl (c) Oscar Dowson and contributors, 2017-23
-------------------------------------------------------------------
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.078257e+06 1.970447e+00 920 1
20 5.494568e+05 2.078257e+06 2.485726e+00 1340 1
30 4.985879e+04 2.078257e+06 4.347530e+00 2260 1
40 3.799447e+06 2.078257e+06 4.915712e+00 2680 1
50 1.049867e+06 2.078257e+06 6.895854e+00 3600 1
60 3.985191e+04 2.078257e+06 7.484916e+00 4020 1
-------------------------------------------------------------------
status : simulation_stopping
total time (s) : 7.484916e+00
total solves : 4020
best bound : 2.078257e+06
simulation ci : 2.031697e+06 ± 3.922745e+05
numeric issues : 0
-------------------------------------------------------------------
-------------------------------------------------------------------
SDDP.jl (c) Oscar Dowson and contributors, 2017-23
-------------------------------------------------------------------
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
-------------------------------------------------------------------
10L 4.986663e+04 2.079056e+06 3.561376e+00 920 1
20L 3.799878e+06 2.079267e+06 6.280585e+00 1340 1
30L 3.003923e+04 2.079457e+06 1.024400e+01 2260 1
40L 5.549882e+06 2.079457e+06 1.292256e+01 2680 1
50L 2.799466e+06 2.079457e+06 1.692570e+01 3600 1
60L 3.549880e+06 2.079457e+06 1.990280e+01 4020 1
-------------------------------------------------------------------
status : simulation_stopping
total time (s) : 1.990280e+01
total solves : 4020
best bound : 2.079457e+06
simulation ci : 2.352197e+06 ± 5.377546e+05
numeric issues : 0
-------------------------------------------------------------------