Basically the dynamic method that samples a fixed set of prices (ala ribs) on the backwards pass in addition to the price sampled on the forward pass.

@lkapelevich this weekend I intend to squash the git history of SDDP.jl to remove the code that was split into the pro version out of the repo. Because you've forked it, I have no idea what that will do to your copy.

What's the ETA on the changes you want to bring in for the Lagrangean stuff?

At the moment it's horrible (because I know nothing about JS).

https://github.com/odow/SDDP.jl/blob/master/src/visualizer/simulation.js

Bonus points for making it easy to embed into other sites such as @EyobZewdie's Part IV report.

Would be nice to have this re-implemented.

Concretely, the dynamics function would change from

foo(price::T, noise::S, stage::Int, markov_state::Int) where T where S

to

foo(price::T, noise::S) where T where S

and you would create stage-dependent value functions via the functional approach

function buildvaluefunction(stage::Int, markov_state::Int)
    PriceInterpolation(
        dynamics = (price, noise) -> price + noise,
        min_price = -t,
        max_price = 2t,
        lipschitz_constant = 100 * i
    )
end
m=SDDPModel(
# other arguments
value_function = buildvaluefunction
) do sp, t, i
# subproblem definition
end

A common pattern is to define a set of noises, converge the model, and then simulate using some out-of-sample observations.

The current approach is to dump the cuts to file when solving the model, then rebuild a different model with the new observations, load the cuts and simulate.

This should be improved.

Should be resolved along with #29

SDDP.jl needs a better name. Suggestions that have been floated

• ELISE.jl - Extensible LIbrary for SDDP Experiments
• SMOKIE.jl - Stochastic Multistage Optimisation Kit Including Extensibility
• ELMO (S).jl - Extensible Library for Multistage Optimisation (with Stochasticity)
• JESS.jl - Julia's Extensible SDDP Solver
• JuDO.jl - Julia for Dynamic Optimization
• SMORE.jl - Stochastic Multi-stage Optimization Research Environment
• JAST.jl - Just Another SDDP Toolbox
• JEMSO.jl - JuMP Extension for Multi-stage Stochastic Optimisation
• Kokako.jl - for no reason other than Kokako's are cool birds
• Mohua.jl - Multistage Optimisation H?? Uncertainty A??
• Moa.jl - Multistage Optimisation Application

At this stage I'm leaning towards Moa.jlKōkako.jl... It riffs nicely with WEKA and Pajarito.jl

@EyobZewdie reports that the cut output file becomes corrupt when using the async solver.

Presumably this is because each process writes cuts to the same file. I'd hoped the OS would put a lock on a file being edited but apparently not.

Running the posm_price.jl example without the line

t2 = deepcopy(t)

results in an error ERROR: type UnionAll has no field name.

I cannot reproduce this yet.

Potentially related
JuliaLang/julia#21700
JuliaLang/julia#20804
JuliaLang/julia#19456

Provided the value function is still convex w.r.t. the state variables, and we can obtain valid subgradients, there is nothing to say that the controls have to be convex.

In a previous iteration of SDDP.jl, we implemented this by solving the MIP's on the forward pass, then on the backwards pass, we solve the MIP, fix the integer variables and then solve the LP relaxation. This doesn't necessarily provide valid subgradients but it will create a policy that might be good enough.

This is particularly useful if you want to do things like use SOSII to model things.

Leaving some code that may be potentially useful.

nonzero(x) = abs(x) > 1e-10

function fixvariable!(m::Model, idx, value)
    setupperbound(Variable(m, idx), value)
    setlowerbound(Variable(m, idx), value)
end

function fixmodel!(m::Model)
    !m.internalModelLoaded && error("Unable to fix an unsolved model")
    # Follow what Gurobi does for Special Ordered Sets
    if length(m.sosconstr) > 0
        for sos in m.sosconstr
            cols = Int[v.col for v in sos.terms]
            if sos.sostype == :SOS1
                for c in cols
                    nonzero(m.colVal[c]) || fixvariable!(m, c, 0.)
                end
            elseif sos.sostype == :SOS2
                colorder = sortperm(sos.weights)
                for (i, sosidx) in enumerate(colorder)
                    c = cols[sosidx]
                    nonzero(m.colVal[c]) || fixvariable!(m, c, 0.)
                end
            else
                error("Unable to create fixed model with SOS of type $(sos.sostype)")
            end
        end
        m.sosconstr = JuMP.SOSConstraint[]
        # Do this for now since no efficient way to update/remove SOS constraints
        m.internalModelLoaded = false
    end
    # Fix all binary and integer variables to their current value
    for (i, v) in enumerate(m.colCat)
        if v == :Bin || v == :Int
            setcategory(Variable(m, i), :Cont)
            fixvariable!(m, i, m.colVal[i])
        end
    end
end

function solvenoncontinuous!(sp::Model, scenario::Int)
    lb = copy(sp.colLower)
    ub = copy(sp.colUpper)
    cat = copy(sp.colCat)
    sos = copy(sp.sosconstr)

    fixmodel!(sp)
    status = solve(sp)
    @assert status == :Optimal

    storecontinuous!(sp, scenario)

    sp.colLower = lb
    sp.colUpper = ub
    sp.colCat = cat
    if length(sos) > 0
        sp.sosconstr = sos
        sp.internalModelLoaded = false
    end
end

function SOSII!(m::JuMP.Model, f::Function, x, lb::Float64, ub::Float64, n::Int64=10)
    #   This function adds a SOS of Type II to the JuMP model
    #      Σλ = 1         # Convexity
    #      a∙λ = z        # x value
    #      f(z) = f(a)∙λ  # approximated function value
    #      where λ is a SOS of Type II
    xx = linspace(lb, ub, n) # Set x value
    fx = f(xx)               # Estimate function value
    SOSII!(m, x, xx, fx)
end

function SOSII!(m::JuMP.Model, x, xx::AbstractVector, fx::AbstractVector)
    n = length(xx)
    @assert length(fx) == n

    y = @variable(m)                                    # Estimated function value variable
    λ = @variable(m, [1:n], lowerbound=0, upperbound=1) # SOS variable set
    @constraints(m, begin
        sum(λ) == 1          # Convexity
        x      == dot(xx, λ)
        y      == dot(fx, λ)
    end)
    addSOS2(m, [i*λ[i] for i =1:n])
    y  # Return estimated function value variable
end

SOSII!{T<:Real}(m::JuMP.Model, x, xx::Vector{NTuple{2, T}}) = SOSII!(m, x, [y[1] for y in xx], [y[2] for y in xx])

function bilinear(x::Variable, y::Variable, n::Int=10)
    # This macro is used to replace a bilinear term in a JuMP constraint with the expansion
    #   xy = 0.25 * [(x+y)^2 - (x-y)^2]
    #     where (x+y)^2 and (x-y)^2 are approximated with SOS2
    m = x.m

    # (x + y)^2 SOS2
    v1 = SOSII!(m, (a) -> a.^2, x + y, getlowerbound(x) + getlowerbound(y), getupperbound(x) + getupperbound(y), n)

    # (x + y)^2 SOS2
    v2 = SOSII!(m, (a) -> a.^2, x - y, getlowerbound(x) - getupperbound(y), getupperbound(x) - getlowerbound(y), n)

    # Return expression to be placed in constraint
    (0.25 * (v1 - v2))
end

Big things that aren't really tested (from CodeCov):

• time limit termination (covered by #41)
• user interruption termination (not sure if possible?)
• asyncronous (problems with coverage on remote workers) (improved by #43, #47)
• visualiser (covered by #42, #46)
• saving the model (covered by #41)

@jgre432 found a bug with makov states with a zero transition probability. Fix in progress

map vs pmap

Hi,

I have an issue with the show function applied in a SimulationPlot variable like is described in draft_paper.pdf. When i use show(simPl) the plot never occurs and i see in console all data structure of the variable.
Has anyone experienced this issue? How did you solve it?

Thanks!
Rodrigo Porteiro

See https://github.com/odow/SDDP.jl/blob/master/src/solve_asyncronous.jl#L257

What's the purpose of both m.storage.states and stage[...].state?
Currently visited states are saved in stage.state but m.storage.states is passed into modifyprobability!.

I don't think there are updated docs clarifying this at the moment. Just leaving this here to answer a question.

The first column is either a single value or a confidence interval from a Monte Carlo simulation. When the entry is a single value, its is the objective from a single (and the most recent) forward pass in SDDP. The single value is not an average.

To compute a statistical upper bound, you can:

1. Use the simulate function after you finished solving your model
2. If you want to get bounds while solving, use the simulation keyword argument in the solve function. You can control how frequently these will run.

There is an example of both in: https://github.com/odow/SDDP.jl/blob/master/examples/Newsvendor/newsvendor.jl.

The first time we solve the model, we set the frequency of simulations to be every 10 iterations:
https://github.com/odow/SDDP.jl/blob/master/examples/Newsvendor/newsvendor.jl#L99
So in the 10th iteration you will see a calculation of % Gap and a confidence interval.

Then there is an example using the simulate function:
https://github.com/odow/SDDP.jl/blob/master/examples/Newsvendor/newsvendor.jl#L134
And in the line below we compute the mean objective of 500 simulations.

Currently we make no claims about compatibility going forward which isn't great if we want to change existing functionality (ref #84).

Many users run into the issue where they use the @noise macro but have the uncertainty in the constraint matrix rather than the RHS vector.

Actions

• throw an error when used in the constraint matrix
• rename to @rhsnoise

Aimed at an audience with some stochastic programming / dynamic programming knowledge but no SDDP background.

Enhancement 1

Instead of rebuilding the model every time we want to change the set of cuts we include, we could try "turning off" contraints by relaxing their RHS.

Since Θ ≤ α + π'x gets re-written as -∞ ≤ Θ - π'x ≤ α, we can just store the α of each constraint, and "turn off" the constraint by replacing it with some big-M (i.e. +∞). Hopefully Gurobi is smart enough to understand that it doesn't need to consider these. This changes the model rebuild from a lot of work to a change in the RHS (fast).

Enhancement 2

Another option is to check the hash of new cuts to avoid adding duplicate cuts: i.e. ask a cut oracle if it should add the cut.

It actually turns out this is only efficient for very small models. Anecdotally, >98% of cuts in a large-scale farming model are unique.

The requests keep piling up for this one.

Waiting until JuMP/MOI is released.

As suggested by @adow031, it might be easier to specify objective bounds and transition matrices via a function that takes as input the stage and markov state instead of forming jagged arrays. For example

function foo(t::Int,i::Int)
    100 * t + i
end

SDDPModel(
    objective_bound = foo
) do sp, t, i

end

Instead of

SDDPModel(
    objective_bound = [ [101, 102, 103], [201, 202, 203], ... ]
) do sp, t, i

end

This would require adding a the following method to https://github.com/odow/SDDP.jl/blob/master/src/utilities.jl#L70

getel(T, f;:Function, t::Int, i::Int) = f(t,i)

along with some tests.

Add a test for #90

The worst observation chosen in https://github.com/odow/SDDP.jl/blob/master/src/avar_riskaversion.jl#L60-L69
may in fact have zero probability of being reached (i.e. from a markov state with zero transition probability).

Disabling one line in JuMP gives me 20-50+% speed improvement across a range of models.

See jump-dev/JuMP.jl#1110 for the discussion and jump-dev/JuMP.jl#1111 for a fix.

Profiling was done as proposed in #27

It would be useful to go through and check numerical stability in each of the subproblem (#6).

• Analyze coefficients for user-warning. Things to check
  • Large coefficients > 1e7
    log10(maximum(abs.(coefficients))) > 7
  • Small coefficients < 1e-2
    -6 < log10(minimum(abs.(coefficients))) < -2
  • High dynamic range
    log10(maximum(abs.(coefficients)) - minimum(abs.(coefficients))) > 8
• Throw away cuts that are numerically almost identical.

We could report something like Gurobi does

Optimize a model with 44 rows, 19 columns and 106 nonzeros
Coefficient statistics:
  Matrix range    [2e-05, 1e+03]
  Objective range [1e-04, 9e+02]
  Bounds range    [1e+01, 4e+05]
  RHS range       [2e+00, 6e+05]

-------------------------------------------------------------------------------
                      SDDP Solver. © Oscar Dowson, 2017.
-------------------------------------------------------------------------------
    Solver:
        Serial solver
    Model:
        Stages:         4
        States:         2
        Subproblems:    7
        Value Function: Default
        Numerical Stability Issues: none
-------------------------------------------------------------------------------

-------------------------------------------------------------------------------
                      SDDP Solver. © Oscar Dowson, 2017.
-------------------------------------------------------------------------------
    Solver:
        Serial solver
    Model:
        Stages:         4
        States:         2
        Subproblems:    7
        Value Function: Default
        Numerical Stability Issues:
            Objective Coefficients: 
                Stage 1, Markov State 1: [-1e-2, +1e13]
             Constraint Coefficients: 
                Stage 13, Markov State 3: [-1e-2, +1e13]
                Stage 14, Markov State 3: [-1e-2, +1e13]
                Stage 14, Markov State 3: [-1e-2, +1e13]
-------------------------------------------------------------------------------

This might also be useful

function getlinconstrmatrix(m)
    A = zeros(length(m.linconstr), MathProgBase.numvar(m))
    for (row, con) in enumerate(m.linconstr)
        for (coef, var) in zip(con.terms.coeffs, con.terms.vars)
            A[row, var.col] += coef
        end
    end
    A
end

Occasionally with

• Gurobi I get status: Numeric
• CPLEX I get InfOrUnBnd
• CPLEX I get CPX_OPT_INFEAS here
• Clp I get crashes, or incorrect solutions.
• Non-repeatable "Infeasible" status using Gurobi v7.0
• Non-repeatable "Unbounded" status using Gurobi v7.0

Action items:

• Isolate cause
• Figure out how to make subproblems more numerically robust

This will be significantly helped by the JuMP re-write.

Hi,

Im developing an example of hydro dams model. I want to use some historical data to generate trajectories extracting series directly from weekly data in a deterministic way. The first trajectory to use gets the first week of history and advances in that serie, the second trajectory gets the second week and so on. If historical data is shorter than trajectory i get modulus and get the initial data again. When trying to develop in this way i cant find a way in the framework to get the trajectory identifier to index the first week that i must take (if first trajectory the initial week must be week 1, if nth trajectory the initial week must be nth % size). Anyone had this issue? Any suggestions to avoid changing framework code?

Thanks,
Rodrigo

When adding a new cut using a cut selection method, check if an existing cut is scheduled to be removed. Then, instead of appending a new constraint, modify the coefficients of the existing constraint to the new values.

Fairly good evidence that this is a good thing to do as both DOASA (@DrAndyPhilpott) and PSR (@joaquimg) report very good results.

We have two options for implementation

1. Wait for MathOptInterface to drop
2. Implement a hacky route through MathProgBase

Can be implemented in conjunction with #70

Waiting until JuMP/MOI is released.

Add some examples with conic subproblems.

At one point I made a start https://github.com/odow/SDDP.jl/tree/jumpsprint/docs

Relevant docs

https://web.stanford.edu/~lcambier/fast/
https://stochdynamicprogramming.readthedocs.io/en/latest/
https://blegat.github.io/StructDualDynProg.jl/latest/
http://www.juliaopt.org/MathOptInterface.jl/latest/

Tick once addressed:

• quick start
• timer outputs
• value functions :( this is a lot of work
• static price interpolation
• bad lipschitz constants
• #65: intro to SDDP
• #113: log output
• #120: Markovian policy graphs
• #121: price interpolation and Markov states
• #132: multiple noise constraints
• DRO #140
• Conic stuff
• plotting utilities
• parallel
• more complex usages of @state
• function input for SDDPModel keywords
• writing cuts to file
• reading cuts from file

Refer to PRs updating the docs: #34, #130 #133 #135 #140 #141

To clean up the price interpolation code.

We can elaborate on this once my thesis is written, saving this here while I remember.

This is a nice way to avoid confusing the whole stage/node/do sp, t/do sp, t, i situation.

struct PolicyGraph{T}
    root::T
    subproblems::Dict{T, Vector{Tuple{T, Float64}}}
end

function LinearGraph(stages::Int)
    subproblems = Dict{Int, Vector{Tuple{Int, Float64}}}()
    subproblems[0] = Tuple{Int, Float64}[]
    for t in 1:stages
        subproblems[t] = Tuple{Int, Float64}[]
        push!(subproblems[t-1], (t, 1.0))
    end
    PolicyGraph(0, subproblems)
end

function MarkovianGraph(transition_probabilities::AbstractVector{AbstractArray{Float64, 2}})
    # check is Markov transition matrix
    i = 1
    for T in transition_probabilities
        @assert size(T, 1) == i
        i = size(T, 2)
        for j in 1:size(T, 1)
            @assert isapprox(sum(T[j,:]), 1.0)
        end
    end

    typ = Tuple{Int, Int}
    subproblems = Dict{typ, Vector{Tuple{typ, Float64}}}()
    subproblems[(0,1)] = Tuple{typ, Float64}[]
    for (t, T) in enumerate(transition_probabilities)
        for j in 1:size(T, 2)
            subproblems[(t, j)] = Tuple{typ, Float64}[]
        end
        for i in 1:size(T, 1)
            for j in 1:size(T, 2)
                push!(subproblems[(t-1, i)], ((t, j), T[i,j]))
            end
        end
    end
    PolicyGraph((0,1), subproblems)
end

SDDPModel(
    policy_graph = LinearGraph(3)
) do sp, index
    stage = index
    # stage = 1,2,3
end

SDDPModel(
    policy_graph = MarkovianGraph([ 
        [1.0]',
        [0.2 0.8],
        [0.5 0.5; 0.1 0.9]
    ])
) do sp, index
    (stage, markov_state) = index
    # (stage, markov_state) = (1,1), (2,1), (2,2), (3,1), (3,2)

end

function children(g::PolicyGraph{T}, subproblem::T) where T
    [s[2] for s in g.subproblems[subproblem]]
end

g = LinearGraph(3)
children(g, 0) # [1]
children(g, 2) # [3]
children(g, 3) # []

g = MarkovianGraph([ 
        [1.0]',
        [0.2 0.8],
        [0.5 0.5; 0.1 0.9]
    ]
children(g, (0,1)) # [(1,1)]
children(g, (2,2)) # [(3,1), (3,2)]
children(g, (3,1)) # []

For example, pass a vector, one for each stage

risk_measure = [NestedAVaR(lambda=0.5, beta=0.5), NestedAVaR(lambda=0.5, beta=0.25)]

A vector of vectors, one for each stage and markov state

risk_measure = [
    [NestedAVaR(lambda=0.5, beta=0.5), NestedAVaR(lambda=0.5, beta=0.25)],
    [Expectation(), NestedAVaR(lambda=0.25, beta=0.25)]
]

Or a function

function stagedependentrisk(stage, markov_state)
    if state == 1
        return Expectation()
    else
        return NestedAVaR(lambda=0.5, beta=markov_state / 4)
    end
end

risk_measure = stagedependentrisk

cc @Thuener

To do this, begin with random forward pass that ignores the unboundedness to avoid specifying an initial bound.

Bonus points for tying into #63

• Tidy up the multitude of addcut related methods to more explanatory
• Change methods to overload on SDDPModel{V} (#26)
• Squash previous dev commits
• Specify scenarios with a sparse b vector Revisit in new JuMP
• Add random stagewise-independent cost vectors (#25)
• Document @visualise (this sneaked into #7)
• addconstraint function rather than macro. Build directly Revisit in new JuMP

Edit: removed old links to ShareLatex.

http://www.optimization-online.org/DB_HTML/2017/12/6388.html

Feedback welcomed. @Thuener @lkapelevich

• Preprint on OptimizationOnline
• Re-write into a proper paper rather than a tutorial style
• Benchmark against DOASA

For example, forward sample with the riskadjusted_distribution of the risk measure, or uniformly.

Implement a stage dependent objective bound as suggested by @adow031

This would require using the getel functions for https://github.com/odow/SDDP.jl/blob/master/src/SDDP.jl#L182

Will look similar to #66

Shouldn't be more than a few lines with tests.

Implement an oracle that can be used to sample the subproblem or noise in the forward pass of the algorithm.

It probably takes the form of a new overloadable type that gets dispatched in the above functions.

Bonus points for solving #10 in addition to closing #29 and #79.

Potentially looks something like

abstract type SamplingOracle end

struct Default <: SamplingOracle end
function samplenoise(::Default, m::SDDPModel, sp::JuMP.Model)
    ex = ext(sp)
    noiseidx = sample(ex.noiseprobability)
    return ex.noises[noiseidx]
end

struct Uniform <: SamplingOracle end
function samplenoise(::Uniform, m::SDDPModel, sp::JuMP.Model)
    return rand(ext(sp).noises)
end

struct Deterministic{T} <: SamplingOracle 
    data::Vector{T}
end
function samplenoise(oracle::Deterministic{T}, m::SDDPModel, sp::JuMP.Model) where T
    oracle.data[ext(sp).stage]
end

# ...

solve(m, sampling_oracle=Default())

The variables status and keep_iterating here:

are overwritten on this line

@rhsnoise(sp, w=[1,2,3], 2x+y<=w)
# becomes
@rhsnoise(sp, [1,2,3], w->2x+y<=w)

This sorts out the using different noise fiasco since we can just throw new observations into the function (ref #10). It plays nicer with namespacing issues

w = 0
@rhsnoise(sp, w=[1,2,3], 2x+y<=w)
# what is w? 3!

w = 0
@rhsnoise(sp, [1,2,3], w->2x+y<=w)
# what is w? 0!

And is probably a good start towards noises in the constraint coefficients. It also slightly unifies our price interpolation stuff since we would have

@stageobjective(sp, [1,2,3], w->w*x)
# and
@stageobjective(sp, p->p*x)

But we should really have

@stageobjective(sp, [1,2,3], (p, w) -> p*x + w)

to allow for stagewise independent noise on-top-of the price states.

Nothing actionable yet, leaving this as a reminder.

In the hydro valley examples, generation only includes the last level instead of all levels. E.g.:
https://github.com/odow/SDDP.jl/blob/master/examples/visualisation.jl#L84

The whole value function abstraction needs to be re-done to streamline the multiple addcut redirections and minimize the code overlap with the new ones in development.

You should be able to set a solve hook on the backwards pass and on the forwards pass

@Thuener wondered how to incorporate a noise in the constraint matrix (i.e. LHS of the constraint).

The answer is that due to the way JuMP is structured, it is difficult to modify constraint coefficients in an efficient manner.

Users should utilize the Markov state feature to achieve random constraint coefficients. For example:

https://github.com/odow/SDDP.jl/blob/master/examples/asset_management.jl#L40

The lipschitz_constant argument to the dynamic interpolation method could take a tuple giving a different constant for each of the different price dimensions.

As suggested by @adow031

Integrate https://github.com/KristofferC/TimerOutputs.jl so users can see where the time is spent. This will also help development.

Greetings,

I'm working on an out-of-sample simulation function based on the existing 'simulate' functions. Currently it receives a list of SDDP.Noise objects (of size equal to the number of stages). I created methods to save the results, including the values of the SDDP.Noise argument in a key temporarily named ':new_noise'.

I think the name ':noise' for the storage key that represents the index for the noise in the model is unclear. Especially if you intend to implement sampling oracles as per #80. I'd like to suggest renaming the ':noise' key to ':noiseidx' or similar. What's your opinion on this?

Thanks,
Vitor.