StructDualDynProg: Problem 5.2, 2 stages
This tutorial was generated using Literate.jl. Download the source as a .jl
file. Download the source as a .ipynb
file.
This example comes from StochasticDualDynamicProgramming.jl
using SDDP, HiGHS, Test
function test_prob52_2stages()
model = SDDP.LinearPolicyGraph(;
stages = 2,
lower_bound = 0.0,
optimizer = HiGHS.Optimizer,
) do subproblem, stage
# ========== Problem data ==========
n = 4
m = 3
i_c = [16, 5, 32, 2]
C = [25, 80, 6.5, 160]
T = [8760, 7000, 1500] / 8760
D2 = [diff([0, 3919, 7329, 10315]) diff([0, 7086, 9004, 11169])]
p2 = [0.9, 0.1]
# ========== State Variables ==========
@variable(subproblem, x[i = 1:n] >= 0, SDDP.State, initial_value = 0.0)
# ========== Variables ==========
@variables(subproblem, begin
y[1:n, 1:m] >= 0
v[1:n] >= 0
penalty >= 0
rhs_noise[1:m] # Dummy variable for RHS noise term.
end)
# ========== Constraints ==========
@constraints(
subproblem,
begin
[i = 1:n], x[i].out == x[i].in + v[i]
[i = 1:n], sum(y[i, :]) <= x[i].in
[j = 1:m], sum(y[:, j]) + penalty >= rhs_noise[j]
end
)
if stage == 2
# No investment in last stage.
@constraint(subproblem, sum(v) == 0)
end
# ========== Uncertainty ==========
if stage != 1 # no uncertainty in first stage
SDDP.parameterize(subproblem, 1:size(D2, 2), p2) do ω
for j in 1:m
JuMP.fix(rhs_noise[j], D2[j, ω])
end
end
end
# ========== Stage objective ==========
@stageobjective(subproblem, i_c' * v + C' * y * T + 1e6 * penalty)
return
end
SDDP.train(model; log_frequency = 10)
@test SDDP.calculate_bound(model) ≈ 340315.52 atol = 0.1
return
end
test_prob52_2stages()
-------------------------------------------------------------------
SDDP.jl (c) Oscar Dowson and contributors, 2017-24
-------------------------------------------------------------------
problem
nodes : 2
state variables : 4
scenarios : 2.00000e+00
existing cuts : false
options
solver : serial mode
risk measure : SDDP.Expectation()
sampling scheme : SDDP.InSampleMonteCarlo
subproblem structure
VariableRef : [29, 29]
AffExpr in MOI.EqualTo{Float64} : [4, 5]
AffExpr in MOI.GreaterThan{Float64} : [3, 3]
AffExpr in MOI.LessThan{Float64} : [4, 4]
VariableRef in MOI.GreaterThan{Float64} : [22, 22]
VariableRef in MOI.LessThan{Float64} : [1, 1]
numerical stability report
matrix range [1e+00, 1e+00]
objective range [1e+00, 1e+06]
bounds range [0e+00, 0e+00]
rhs range [0e+00, 0e+00]
-------------------------------------------------------------------
iteration simulation bound time (s) solves pid
-------------------------------------------------------------------
10 3.455904e+05 3.147347e+05 8.118153e-03 54 1
20 3.336455e+05 3.402383e+05 1.420712e-02 104 1
30 3.337559e+05 3.403155e+05 2.133512e-02 158 1
40 3.337559e+05 3.403155e+05 2.817106e-02 208 1
48 3.337559e+05 3.403155e+05 3.657508e-02 248 1
-------------------------------------------------------------------
status : simulation_stopping
total time (s) : 3.657508e-02
total solves : 248
best bound : 3.403155e+05
simulation ci : 1.351676e+08 ± 1.785770e+08
numeric issues : 0
-------------------------------------------------------------------