Air conditioning
This tutorial was generated using Literate.jl. Download the source as a .jl file. Download the source as a .ipynb file.
Taken from Anthony Papavasiliou's notes on SDDP This is a variation of the problem that first appears in the book Introduction to Stochastic Programming by Birge and Louveaux, 1997, Springer-Verlag, New York, on page 237, Example 1. For a rescaled problem, they reported an optimal value of 6.25 with a first-stage solution of x1 = 2 (production)and y1 = 1 (store production). On this variation, without rescaling, it would be equivalent to 62500, 200 and 100, respectively.
Consider the following problem
- Produce air conditioners for 3 months
- 200 units/month at 100 $/unit
- Overtime costs 300 $/unit
- Known demand of 100 units for period 1
- Equally likely demand, 100 or 300 units, for periods 2, 3
- Storage cost is 50 $/unit
- All demand must be met
The known optimal solution is $62,500
using SDDP, HiGHS, Testfunction air_conditioning_model(duality_handler) model = SDDP.LinearPolicyGraph(; stages = 3, lower_bound = 0.0, optimizer = HiGHS.Optimizer, ) do sp, stage @variable( sp, 0 <= stored_production <= 100, Int, SDDP.State, initial_value = 0 ) @variable(sp, 0 <= production <= 200, Int) @variable(sp, overtime >= 0, Int) @variable(sp, demand) DEMAND = [[100.0], [100.0, 300.0], [100.0, 300.0]] SDDP.parameterize(ω -> fix(demand, ω), sp, DEMAND[stage]) @constraint( sp, stored_production.out == stored_production.in + production + overtime - demand ) @stageobjective( sp, 100 * production + 300 * overtime + 50 * stored_production.out ) end SDDP.train(model; duality_handler = duality_handler) lb = SDDP.calculate_bound(model) println("Lower bound is: $lb") @test isapprox(lb, 62_500.0, atol = 0.1) sims = SDDP.simulate(model, 1, [:production, :stored_production]) x1 = sims[1][1][:production] y1 = sims[1][1][:stored_production].out @test isapprox(x1, 200, atol = 0.1) @test isapprox(y1, 100, atol = 0.1) println( "With first stage solutions $(x1) (production) and $(y1) (stored_production).", ) returnendfor duality_handler in [SDDP.LagrangianDuality(), SDDP.ContinuousConicDuality()] air_conditioning_model(duality_handler)end-------------------------------------------------------------------
SDDP.jl (c) Oscar Dowson and contributors, 2017-26
-------------------------------------------------------------------
problem
nodes : 3
state variables : 1
scenarios : 4.00000e+00
existing cuts : false
options
solver : serial mode
risk measure : SDDP.Expectation()
sampling scheme : SDDP.InSampleMonteCarlo
subproblem structure
VariableRef : [6, 6]
AffExpr in MOI.EqualTo{Float64} : [1, 1]
VariableRef in MOI.GreaterThan{Float64} : [4, 4]
VariableRef in MOI.Integer : [3, 3]
VariableRef in MOI.LessThan{Float64} : [2, 3]
numerical stability report
matrix range [1e+00, 1e+00]
objective range [1e+00, 3e+02]
bounds range [1e+02, 2e+02]
rhs range [0e+00, 0e+00]
-------------------------------------------------------------------
iteration simulation bound time (s) solves pid
-------------------------------------------------------------------
1L 1.100000e+05 6.250000e+04 4.911990e-01 8 1
20L 9.500000e+04 6.250000e+04 1.090166e+00 172 1
-------------------------------------------------------------------
status : simulation_stopping
total time (s) : 1.090166e+00
total solves : 172
best bound : 6.250000e+04
numeric issues : 0
-------------------------------------------------------------------
Lower bound is: 62500.0
With first stage solutions 200.0 (production) and 100.0 (stored_production).
-------------------------------------------------------------------
SDDP.jl (c) Oscar Dowson and contributors, 2017-26
-------------------------------------------------------------------
problem
nodes : 3
state variables : 1
scenarios : 4.00000e+00
existing cuts : false
options
solver : serial mode
risk measure : SDDP.Expectation()
sampling scheme : SDDP.InSampleMonteCarlo
subproblem structure
VariableRef : [6, 6]
AffExpr in MOI.EqualTo{Float64} : [1, 1]
VariableRef in MOI.GreaterThan{Float64} : [4, 4]
VariableRef in MOI.Integer : [3, 3]
VariableRef in MOI.LessThan{Float64} : [2, 3]
numerical stability report
matrix range [1e+00, 1e+00]
objective range [1e+00, 3e+02]
bounds range [1e+02, 2e+02]
rhs range [0e+00, 0e+00]
-------------------------------------------------------------------
iteration simulation bound time (s) solves pid
-------------------------------------------------------------------
1 3.000000e+04 6.250000e+04 3.396034e-03 8 1
20 5.500000e+04 6.250000e+04 2.640460e-01 172 1
-------------------------------------------------------------------
status : simulation_stopping
total time (s) : 2.640460e-01
total solves : 172
best bound : 6.250000e+04
numeric issues : 0
-------------------------------------------------------------------
Lower bound is: 62500.0
With first stage solutions 200.0 (production) and 100.0 (stored_production).