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  5.343089e-01       920   1
        20   5.494568e+05  2.078257e+06  7.377729e-01      1340   1
        30   4.985879e+04  2.078257e+06  1.293831e+00      2260   1
        40   3.799447e+06  2.078257e+06  1.499373e+00      2680   1
        50   1.049867e+06  2.078257e+06  2.069469e+00      3600   1
        60   3.985191e+04  2.078257e+06  2.278581e+00      4020   1
-------------------------------------------------------------------
status         : simulation_stopping
total time (s) : 2.278581e+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.079119e+06  9.651542e-01       920   1
        20L  3.799878e+06  2.079330e+06  1.691671e+00      1340   1
        30L  3.003923e+04  2.079457e+06  2.824466e+00      2260   1
        40L  5.549882e+06  2.079457e+06  3.636492e+00      2680   1
        50L  2.799466e+06  2.079457e+06  4.829626e+00      3600   1
        60L  3.549880e+06  2.079457e+06  5.610949e+00      4020   1
-------------------------------------------------------------------
status         : simulation_stopping
total time (s) : 5.610949e+00
total solves   : 4020
best bound     :  2.079457e+06
simulation ci  :  2.352204e+06 ± 5.377531e+05
numeric issues : 0
-------------------------------------------------------------------