Multi-Stage Instantiation

[rtwalker]

I ran into an error with instantiating multistage models.

Taking the inventory example from the multistage chapter of Birge and Louveaux:

inventory_model = @stochastic_model begin
    @stage 1 begin
        # number of units produced in regular time work during month 1
        @variable(model, x₁ >= 0)
        # number of units produced in overtime work during month 1
        @variable(model, w₁ >= 0)
        # number of units stored at month 1
        @decision(model, y₁ >= 0)
        # minimize total cost = production cost + holding cost
        @objective(model, Min, x₁ + 3*w₁ + 0.5*y₁)
        # constraint on production capacity
        @constraint(model, x₁ <= 2)
        # must meet already known month 1 demand
        @constraint(model, x₁ + w₁ - y₁ == 1)
    end

    @stage 2 begin
        @uncertain d
        # number of units produced in regular time work during month 2
        @variable(model, x₂ >= 0)
        # number of units produced in overtime work during month 2
        @variable(model, w₂ >= 0)
        # number of units stored at month 2
        @decision(model, y₂ >= 0)
        # minimize total cost = production cost + holding cost
        @objective(model, Min, x₂ + 3*w₂ + 0.5*y₂)
        # constraint on production capacity
        @constraint(model, x₂ <= 2)
        # must meet month 2 demand
        @constraint(model, x₂ + w₂ - y₂ == d - y₁)
    end

    @stage 3 begin
        @uncertain d
        # number of units produced in regular time work during month 3
        @variable(model, x₃ >= 0)
        # number of units produced in overtime work during month 3
        @variable(model, w₃ >= 0)
        # number of units stored at month 3
        @decision(model, y₃ >= 0)
        # minimize total cost = production cost + holding cost
        @objective(model, Min, x₃ + 3*w₃ + 0.5*y₃)
        # constraint on production capacity
        @constraint(model, x₃ <= 2)
        # must meet month 3 demand
        @constraint(model, x₃ + w₃ - y₃ == d - y₂)
    end
end

s₂₁ = Scenario(d = 1.0, probability = 0.5)
s₂₂ = Scenario(d = 3.0, probability = 0.5)

s₃₁ = Scenario(d = 1.0, probability = 0.25)
s₃₂ = Scenario(d = 3.0, probability = 0.25)
s₃₃ = Scenario(d = 1.0, probability = 0.25)
s₃₄ = Scenario(d = 3.0, probability = 0.25)

Instantiating the DEP gives the following error:

julia> inventory_dep = instantiate(inventory_model, ([s₂₁,s₂₂], [s₃₁,s₃₂,s₃₃,s₃₄]))
ERROR: KeyError: key :y₂ not found
Stacktrace:
 [1] getindex at /Users/ryanwalker/.julia/packages/JuMP/YXK4e/src/JuMP.jl:936 [inlined]
 [2] (::var"#9#16")(::Model, ::NamedTuple{(),Tuple{}}, ::Scenario{NamedTuple{(:d,),Tuple{Float64}}}) at ./none:0
 [3] generate!(::StochasticProgram{3,Tuple{StochasticPrograms.Stage{NamedTuple{(),Tuple{}}},StochasticPrograms.Stage{NamedTuple{(),Tuple{}}},StochasticPrograms.Stage{NamedTuple{(),Tuple{}}}},DeterministicEquivalent{3,2,Tuple{Array{Scenario{NamedTuple{(:d,),Tuple{Float64}}},1},Array{Scenario{NamedTuple{(:d,),Tuple{Float64}}},1}}}}, ::DeterministicEquivalent{3,2,Tuple{Array{Scenario{NamedTuple{(:d,),Tuple{Float64}}},1},Array{Scenario{NamedTuple{(:d,),Tuple{Float64}}},1}}}, ::Int64) at /Users/ryanwalker/.julia/dev/StochasticPrograms/src/methods/deterministic_equivalent/generation.jl:44
 [4] generate!(::StochasticProgram{3,Tuple{StochasticPrograms.Stage{NamedTuple{(),Tuple{}}},StochasticPrograms.Stage{NamedTuple{(),Tuple{}}},StochasticPrograms.Stage{NamedTuple{(),Tuple{}}}},DeterministicEquivalent{3,2,Tuple{Array{Scenario{NamedTuple{(:d,),Tuple{Float64}}},1},Array{Scenario{NamedTuple{(:d,),Tuple{Float64}}},1}}}}, ::DeterministicEquivalent{3,2,Tuple{Array{Scenario{NamedTuple{(:d,),Tuple{Float64}}},1},Array{Scenario{NamedTuple{(:d,),Tuple{Float64}}},1}}}) at /Users/ryanwalker/.julia/dev/StochasticPrograms/src/methods/deterministic_equivalent/generation.jl:11
 [5] generate!(::StochasticProgram{3,Tuple{StochasticPrograms.Stage{NamedTuple{(),Tuple{}}},StochasticPrograms.Stage{NamedTuple{(),Tuple{}}},StochasticPrograms.Stage{NamedTuple{(),Tuple{}}}},DeterministicEquivalent{3,2,Tuple{Array{Scenario{NamedTuple{(:d,),Tuple{Float64}}},1},Array{Scenario{NamedTuple{(:d,),Tuple{Float64}}},1}}}}) at /Users/ryanwalker/.julia/dev/StochasticPrograms/src/methods/generation.jl:92
 [6] instantiate(::StochasticModel{3,Tuple{StageParameters{NamedTuple{(),Tuple{}}},StageParameters{NamedTuple{(),Tuple{}}},StageParameters{NamedTuple{(),Tuple{}}}}}, ::Tuple{Array{Scenario{NamedTuple{(:d,),Tuple{Float64}}},1},Array{Scenario{NamedTuple{(:d,),Tuple{Float64}}},1}}; instantiation::StochasticPrograms.UnspecifiedInstantiation, optimizer::Nothing, defer::Bool, kw::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at /Users/ryanwalker/.julia/dev/StochasticPrograms/src/methods/api.jl:93
 [7] instantiate(::StochasticModel{3,Tuple{StageParameters{NamedTuple{(),Tuple{}}},StageParameters{NamedTuple{(),Tuple{}}},StageParameters{NamedTuple{(),Tuple{}}}}}, ::Tuple{Array{Scenario{NamedTuple{(:d,),Tuple{Float64}}},1},Array{Scenario{NamedTuple{(:d,),Tuple{Float64}}},1}}) at /Users/ryanwalker/.julia/dev/StochasticPrograms/src/methods/api.jl:81
 [8] top-level scope at REPL[10]:1

The LShaped.Optimizer seems to handle the same model correctly however:

julia> inventory_lshaped = instantiate(inventory_model, ([s₂₁,s₂₂], [s₃₁,s₃₂,s₃₃,s₃₄]), optimizer = LShaped.Optimizer)
Stochastic program with:
 * Stage 1:
   * 1 decision variable
 * Stage 2:
   * 1 decision variable
   * 2 scenarios of type Scenario
 * Stage 3:
   * 4 scenarios of type Scenario
Structure: Vertical
Solver name: L-shaped with disaggregate cuts

[martinbiel]

There is no actual support for multi-stage models in StochasticPrograms yet. There is some experimental groundwork, which is likely to change. Hence, you can instantiate the model in vertical form, but the result is somewhat meaningless and you cannot solve the model.

I might open an issue with a roadmap for multi-stage support when I have more time for development again :)

[odow]

If you want to solve multistage programs, you can use https://github.com/odow/SDDP.jl.

At some point I'd like to separate the SDDP backend so it could be used as a solver for @martinbiel's user-friendly modeling layer.

[martinbiel]

That would be great :)

[rtwalker]

Hence, you can instantiate the model in vertical form, but the result is somewhat meaningless and you cannot solve the model.

I'm working on a nested L-Shaped implementation at the moment, so the instantiation part gets me very close to what I was looking for. Previously most of my trouble was coming from specifying the model, but drawing on the "user-friendly modeling layer" here has alleviated nearly all of that trouble.

In other words, the multistage functionality is good enough for me at the moment, so I'm ok to close this issue given that the above isn't really a bug for the time being.

[odow]

I'm working on a nested L-Shaped implementation at the moment

What for? Any reason to used nested Benders over SDDP ?

[rtwalker]

What for? Any reason to used nested Benders over SDDP ?

Mostly just to have as a point of comparison. I'll try to give SDDP a spin too!

[odow]

You probably already know this, but SDDP is just nested Benders with sampling. For very small problems (e.g., 3 stage), SDDP may be a little slower due to unlucky sampling. But it should only take a few seconds to solve so it's not a big deal. For larger problems, nested Benders won't work. What applications are you trying to solve?

[rtwalker]

Eventually it'll be comparing solution methods for generation expansion planning problems. In the meantime, it's just persuading others that Julia/JuMP/StochasticPrograms would be a good choice for modeling problems and implementing some new ideas.

[odow]

persuading others that Julia/JuMP/StochasticPrograms would be a good choice for modeling problems and implementing some new ideas.

This 💯

We're starting to have a nice ecosystem of tooling that makes it easy to do some crazy stuff. E.g., @andrewrosemberg was solving multistage stochastic hydro-thermal scheduling problems with SDP subproblems: https://proceedings.juliacon.org/papers/10.21105/jcon.00035

There are also tools like: https://github.com/reganbaucke/JuDGE.jl

And we're in the process of starting a file format to collect instances to make implementation of new ideas easier: https://odow.github.io/StochOptFormat/

(Edit: sorry for the noise @martinbiel)