Integer variables when constructing a PolicyGraph

[hideakiv]

Thank you all for maintaining this package! I have a question.

I am trying to write a package involving mixed-integer variables and am testing if I can use the interface in SDDP.jl to construct a multistage model. I set up the control variables such that @variable(subproblem, x, Int) but when I write_subproblem_to_file, the variables are x_ free. Does this mean that the integer constraint is lost?

[odow]

Do you have a minimal working example demonstrating the problem?

Did you write the subproblems to file during a train? We relax the integrality during a train call so we can get dual information: https://github.com/odow/SDDP.jl/blob/a41baa39912ad11ccd3088f8c20906b639d754d0/src/algorithm.jl#L994-L996

I've been meaning to change this so we only relax on the backward pass, and retain integrality on the forward pass (#431).

When stimulating, we use the integer policy.

I see you have: https://github.com/hideakiv/DRMSMIP.jl

SDDP.jl can already solve distributionally robust multistage stochastic integer problems. Can you explain what's different about DRMSMIP?

[hideakiv]

Thank you for the answer! We are focusing on problems that involve general mixed-integer state variables, which could be challenging for SDDiP algorithm. I will construct a minimal example and post it again.

[odow]

SDDP.jl works with mixed-integer state variables. (Although it is not guaranteed to converge to the optimal policy for some problems.)

Re SDDiP: https://github.com/odow/SDDP.jl/issues/246.

Are you working with Kibaek? I'd be interested in a call to discuss in more detail.

[odow]

Here's an example of an SMIP with L1-Wasserstein:

using SDDP, GLPK, Test

model = SDDP.LinearPolicyGraph(
    stages = 3,
    lower_bound = 0.0,
    optimizer = GLPK.Optimizer,
) do sp, t
    DEMAND = [[100.0], [100.0, 300.0], [100.0, 300.0]]
    @variable(sp, 0 <= x <= 100, Int, SDDP.State, initial_value = 0)
    @variable(sp, 0 <= u <= 200, Int)
    @variable(sp, v >= 0, Int)
    @constraint(sp, demand, x.out == x.in + u + v)
    @stageobjective(sp, 100u + 300v + 50x.out)
    SDDP.parameterize(ω -> JuMP.set_normalized_rhs(demand, -ω), sp, DEMAND[t])
end

SDDP.train(
    model;
    iteration_limit = 20,
    risk_measure = SDDP.Wasserstein(GLPK.Optimizer; alpha=50.0) do x, y
        return abs(x.term - y.term)
    end,
)

Here's what I see when I write out a subproblem:

julia> model = SDDP.LinearPolicyGraph(
           stages = 3,
           lower_bound = 0.0,
           optimizer = GLPK.Optimizer,
       ) do sp, t
           DEMAND = [[100.0], [100.0, 300.0], [100.0, 300.0]]
           @variable(sp, 0 <= x <= 100, Int, SDDP.State, initial_value = 0)
           @variable(sp, 0 <= u <= 200, Int)
           @variable(sp, v >= 0, Int)
           @constraint(sp, demand, x.out == x.in + u + v)
           @stageobjective(sp, 100u + 300v + 50x.out)
           SDDP.parameterize(ω -> JuMP.set_normalized_rhs(demand, -ω), sp, DEMAND[t])
       end
A policy graph with 3 nodes.
 Node indices: 1, 2, 3

julia> JuMP.write_to_file(model[1].subproblem, "test.lp")

shell> cat test.lp
minimize
obj: 
subject to
demand: 1 x_out - 1 u - 1 v - 1 x_in = 0
Bounds
x_out <= 100
u <= 200
x_out >= 0
u >= 0
v >= 0
x4 >= 0
x_in free
General
x_out
u
v
End

Not that the x_in is free, but that's because it's just a place-holder for the incoming state variable.

Also note that the objective hasn't been set yet. You need to call parameterize. But this will also add the cost-to-go variable:

julia> SDDP.parameterize(model[1], 100.0)

julia> JuMP.write_to_file(model[1].subproblem, "test.lp")

shell> cat test.lp
minimize
obj: 100 u + 300 v + 50 x_out + 1 x4
subject to
demand: 1 x_out - 1 u - 1 v - 1 x_in = -100
Bounds
x_out <= 100
u <= 200
x_out >= 0
u >= 0
v >= 0
x4 >= 0
x_in free
General
x_out
u
v
End

[odow]

Since it looks like you might also be interested in partially observable stuff (https://arxiv.org/pdf/1906.05988.pdf), this paper may interest you: https://www.sciencedirect.com/science/article/abs/pii/S0167637720300870

Here's SDDP.jl solving a distributional robust partially observable multistage stochastic integer program:

julia> using SDDP, GLPK, Test

julia> graph = SDDP.MarkovianGraph(
           [convert(Matrix{Float64}, [1.0]'), [0.5 0.5], [0.8 0.2; 0.2 0.8]],
       );

julia> SDDP.add_ambiguity_set(graph, [(1,1)], 0.0)

julia> for t in 2:3
           SDDP.add_ambiguity_set(graph, [(t, 1), (t, 2)], 1e3)
       end

julia> graph
Root
 (0, 1)
Nodes
 (1, 1)
 (2, 1)
 (2, 2)
 (3, 1)
 (3, 2)
Arcs
 (0, 1) => (1, 1) w.p. 1.0
 (1, 1) => (2, 1) w.p. 0.5
 (1, 1) => (2, 2) w.p. 0.5
 (2, 1) => (3, 1) w.p. 0.8
 (2, 1) => (3, 2) w.p. 0.2
 (2, 2) => (3, 1) w.p. 0.2
 (2, 2) => (3, 2) w.p. 0.8
Partition
 {
    (1, 1)
 } {
    (2, 1)
    (2, 2)
 } {
    (3, 1)
    (3, 2)
 }

julia> model = SDDP.PolicyGraph(
           graph,
           lower_bound = 0.0,
           optimizer = GLPK.Optimizer,
       ) do sp, node
           t, i = node
           DEMAND = [[100.0], [100.0, 300.0], [100.0, 300.0]]
           P = [
               [[1.0], [0.2, 0.8], [0.2, 0.8]],
               [[1.0], [0.8, 0.2], [0.8, 0.2]],
           ]
           @variable(sp, 0 <= x <= 100, Int, SDDP.State, initial_value = 0)
           @variable(sp, 0 <= u <= 200, Int)
           @variable(sp, v >= 0, Int)
           @constraint(sp, demand, x.out == x.in + u + v)
           @stageobjective(sp, 100u + 300v + 50x.out)
           SDDP.parameterize(sp, DEMAND[t], P[i][t]) do ω
               JuMP.set_normalized_rhs(demand, -ω)
           end
       end
A policy graph with 5 nodes.
 Node indices: (1, 1), (2, 1), (2, 2), (3, 1), (3, 2)

julia> SDDP.train(
           model;
           iteration_limit = 20,
           risk_measure = SDDP.Wasserstein(GLPK.Optimizer; alpha=50.0) do x, y
               return abs(x.term - y.term)
           end,
       )
------------------------------------------------------------------------------
                      SDDP.jl (c) Oscar Dowson, 2017-21

Problem
  Nodes           : 5
  State variables : 1
  Scenarios       : 1.60000e+01
  Existing cuts   : false
  Subproblem structure                      : (min, max)
    Variables                               : (6, 7)
    VariableRef in MOI.LessThan{Float64}    : (3, 5)
    AffExpr in MOI.LessThan{Float64}        : (2, 2)
    AffExpr in MOI.GreaterThan{Float64}     : (1, 2)
    AffExpr in MOI.EqualTo{Float64}         : (1, 1)
    VariableRef in MOI.GreaterThan{Float64} : (5, 6)
    VariableRef in MOI.Integer              : (3, 3)
Options
  Solver          : serial mode
  Risk measure    : SDDP.Wasserstein
  Sampling scheme : SDDP.InSampleMonteCarlo

Numerical stability report
  Non-zero Matrix range     [1e+00, 1e+00]
  Non-zero Objective range  [1e+00, 3e+02]
  Non-zero Bounds range     [1e+02, 1e+03]
  Non-zero RHS range        [1e+02, 3e+02]
No problems detected

 Iteration    Simulation       Bound         Time (s)    Proc. ID   # Solves
        1    3.000000e+04   7.030000e+04   9.419203e-02          1         12
        2    4.000000e+04   7.030000e+04   9.842896e-02          1         24
        3    5.500000e+04   7.118235e+04   1.016860e-01          1         36
        4    9.632353e+04   7.120000e+04   1.053181e-01          1         48
        5    5.500000e+04   7.120000e+04   1.089039e-01          1         60
        6    4.000000e+04   7.120000e+04   1.118469e-01          1         72
        7    4.000000e+04   7.120000e+04   1.148059e-01          1         84
        8    6.000000e+04   7.120000e+04   1.178849e-01          1         96
        9    4.000000e+04   7.120000e+04   1.208880e-01          1        108
       10    4.000000e+04   7.120000e+04   1.239669e-01          1        120
       11    9.500000e+04   7.120000e+04   1.269560e-01          1        132
       12    4.000000e+04   7.120000e+04   1.298389e-01          1        144
       13    4.000000e+04   7.120000e+04   1.328020e-01          1        156
       14    5.500000e+04   7.120000e+04   1.359589e-01          1        168
       15    4.000000e+04   7.120000e+04   1.392291e-01          1        180
       16    4.000000e+04   7.120000e+04   1.424029e-01          1        192
       17    5.500000e+04   7.120000e+04   1.453950e-01          1        204
       18    4.000000e+04   7.120000e+04   1.483951e-01          1        216
       19    6.000000e+04   7.120000e+04   1.512909e-01          1        228
       20    4.000000e+04   7.120000e+04   1.542110e-01          1        240

Terminating training
  Status         : iteration_limit
  Total time (s) : 1.542110e-01
  Total solves   : 240
  Best bound     :  7.120000e+04
  Simulation CI  :  5.006618e+04 ± 7.769264e+03
------------------------------------------------------------------------------

julia> simulations = SDDP.simulate(model, 3, [:x]);

julia> simulations[1][3]
Dict{Symbol, Any} with 7 entries:
  :bellman_term    => 0.0
  :noise_term      => 100.0
  :node_index      => (3, 2)
  :stage_objective => 0.0
  :objective_state => nothing
  :belief          => Dict((3, 2)=>0.894737, (3, 1)=>0.105263, (1, 1)=>0.0, (0, 1)=>0.0, (2, 2)=>0.0, (2, 1)=>0.0)
  :x               => State{Float64}(100.0, 0.0)

[hideakiv]

Thank you! Upon discussing with Kibaek, we decided to use other packages that would suit our needs... I really appreciate the examples, and I would definitely consider using them in my future projects!

As for the SMIP L1-Wasserstein example, does "General" indicate that variables x_out, u, and v are integers?

[odow]

does "General" indicate that variables x_out, u, and v are integers?

Yes. It's the LP file format: https://www.ibm.com/docs/en/icos/12.7.1.0?topic=representation-mip-features-in-lp-file-format

How are you proposing to solve the distributionally robust SMIP? Is it distributional robust against the end-of-horizon distribution? So PH with different weights + branch-and-bound?

we decided to use other packages that would suit our needs

The NREL folks have https://github.com/NREL/ProgressiveHedging.jl. We should try co-ordinate on packages, rather than having a plethora of solvers and modeling interfaces.

To interface with ProgressiveHedging.jl, I just need to write a function that takes a policy graph and returns a list of the JuMP models for each scenario, along with meta-data about the state variables in each stage.

Is this something you need also?

You should also take a look at @martinbiel's https://github.com/martinbiel/StochasticPrograms.jl

[odow]

Closing because there isn't anything to do on the SDDP.jl side of things here. If you want to discuss SMIP further, let's use https://github.com/odow/SDDP.jl/issues/246