Asset management with modifications

This tutorial was generated using Literate.jl. Download the source as a .jl file. Download the source as a .ipynb file.

A modified version of the Asset Management Problem Taken from the book J.R. Birge, F. Louveaux, Introduction to Stochastic Programming, Springer Series in Operations Research and Financial Engineering, Springer New York, New York, NY, 2011

using SDDP, HiGHS, Test

function asset_management_stagewise(; cut_type)
    w_s = [1.25, 1.06]
    w_b = [1.14, 1.12]
    Phi = [-1, 5]
    Psi = [0.02, 0.0]

    model = SDDP.MarkovianPolicyGraph(
        sense = :Max,
        transition_matrices = Array{Float64,2}[
            [1.0]',
            [0.5 0.5],
            [0.5 0.5; 0.5 0.5],
            [0.5 0.5; 0.5 0.5],
        ],
        upper_bound = 1000.0,
        optimizer = HiGHS.Optimizer,
    ) do subproblem, node
        t, i = node
        @variable(subproblem, xs >= 0, SDDP.State, initial_value = 0)
        @variable(subproblem, xb >= 0, SDDP.State, initial_value = 0)
        if t == 1
            @constraint(subproblem, xs.out + xb.out == 55 + xs.in + xb.in)
            @stageobjective(subproblem, 0)
        elseif t == 2 || t == 3
            @variable(subproblem, phi)
            @constraint(
                subproblem,
                w_s[i] * xs.in + w_b[i] * xb.in + phi == xs.out + xb.out
            )
            SDDP.parameterize(subproblem, [1, 2], [0.6, 0.4]) do ω
                JuMP.fix(phi, Phi[ω])
                @stageobjective(subproblem, Psi[ω] * xs.out)
            end
        else
            @variable(subproblem, u >= 0)
            @variable(subproblem, v >= 0)
            @constraint(
                subproblem,
                w_s[i] * xs.in + w_b[i] * xb.in + u - v == 80,
            )
            @stageobjective(subproblem, -4u + v)
        end
    end
    SDDP.train(
        model;
        cut_type = cut_type,
        log_frequency = 10,
        risk_measure = (node) -> begin
            if node[1] != 3
                SDDP.Expectation()
            else
                SDDP.EAVaR(lambda = 0.5, beta = 0.5)
            end
        end,
    )
    @test SDDP.calculate_bound(model) ≈ 1.278 atol = 1e-3
    return
end

asset_management_stagewise(cut_type = SDDP.SINGLE_CUT)

asset_management_stagewise(cut_type = SDDP.MULTI_CUT)
-------------------------------------------------------------------
         SDDP.jl (c) Oscar Dowson and contributors, 2017-23
-------------------------------------------------------------------
problem
  nodes           : 7
  state variables : 2
  scenarios       : 3.20000e+01
  existing cuts   : false
options
  solver          : serial mode
  risk measure    : #4
  sampling scheme : SDDP.InSampleMonteCarlo
subproblem structure
  VariableRef                             : [5, 7]
  AffExpr in MOI.EqualTo{Float64}         : [1, 1]
  VariableRef in MOI.GreaterThan{Float64} : [2, 5]
  VariableRef in MOI.LessThan{Float64}    : [1, 1]
numerical stability report
  matrix range     [1e+00, 1e+00]
  objective range  [2e-02, 4e+00]
  bounds range     [1e+03, 1e+03]
  rhs range        [6e+01, 8e+01]
-------------------------------------------------------------------
 iteration    simulation      bound        time (s)     solves  pid
-------------------------------------------------------------------
        10   4.375000e+00  6.798204e+00  1.585679e-01       278   1
        20   1.230957e+01  1.358825e+00  1.765730e-01       428   1
        30   8.859026e+00  1.278410e+00  2.106240e-01       706   1
        40  -2.315795e+01  1.278410e+00  2.374361e-01       856   1
        49   3.014193e+01  1.278410e+00  2.627370e-01       991   1
-------------------------------------------------------------------
status         : simulation_stopping
total time (s) : 2.627370e-01
total solves   : 991
best bound     :  1.278410e+00
simulation ci  : -1.755629e+00 ± 5.526921e+00
numeric issues : 0
-------------------------------------------------------------------

-------------------------------------------------------------------
         SDDP.jl (c) Oscar Dowson and contributors, 2017-23
-------------------------------------------------------------------
problem
  nodes           : 7
  state variables : 2
  scenarios       : 3.20000e+01
  existing cuts   : false
options
  solver          : serial mode
  risk measure    : #4
  sampling scheme : SDDP.InSampleMonteCarlo
subproblem structure
  VariableRef                             : [5, 7]
  AffExpr in MOI.EqualTo{Float64}         : [1, 1]
  VariableRef in MOI.GreaterThan{Float64} : [2, 5]
  VariableRef in MOI.LessThan{Float64}    : [1, 1]
numerical stability report
  matrix range     [1e+00, 1e+00]
  objective range  [2e-02, 4e+00]
  bounds range     [1e+03, 1e+03]
  rhs range        [6e+01, 8e+01]
-------------------------------------------------------------------
 iteration    simulation      bound        time (s)     solves  pid
-------------------------------------------------------------------
        10  -2.454000e+00  1.647363e+00  3.884983e-02       278   1
        20  -3.575118e+00  1.278410e+00  6.837082e-02       428   1
        30  -5.003795e+01  1.278410e+00  1.220548e-01       706   1
        40   6.835609e+00  1.278410e+00  1.679418e-01       856   1
-------------------------------------------------------------------
status         : simulation_stopping
total time (s) : 1.679418e-01
total solves   : 856
best bound     :  1.278410e+00
simulation ci  :  4.369345e+00 ± 4.780393e+00
numeric issues : 0
-------------------------------------------------------------------