Documentation improvements

[odow]

Overview

I want to make a big overhaul of the documentation. At the moment, we have a lot of pointless tutorial/examples.

We should instead create a few, more useful ones, along the lines of:

• Two-stage (news vendor)
  • [docs] add two-stage newsvendor tutorial #557
• Three-stage (airline problem)
  • Add airline scheduling tutorial #529
• MDP
  • Can a policy graph model MDPs? #505
  • The simples dynamic problem that you can solve in sddp.jl. #513
• Finance (something with objective states or a Markov chain)
  • Add Markovian discretization of finance models #428
  • Document how to add auto-regressive stochastic process #443
• Optimal control (Finance) [docs] add example_reservoir.jl #572
  • Add examples of Stochastic Optimal Control #429
  • Simple Infinite Horizon Discrete Time non-stochastic problem #528
• Hydro-thermal scheduling [docs] add example_reservoir.jl #572

Other

@bernardokp had these suggestions

• <img width="564" alt="image" src="https://user-images.githubusercontent.com/8177701/153268333-857997e4-11b7-4401-9499-fd3b73309404.png">
• 
• 
• 

[odow]

From #516

Doubts from the Vector auto-regressive model,
model = SDDP.LinearPolicyGraph(
    stages = 3,
    sense = :Min,
    lower_bound = 0.0,
    optimizer = GLPK.Optimizer,
) do sp, t
    @variable(sp, 0 <= x <= 200, SDDP.State, initial_value = 200)
    @variable(sp, g_t >= 0)
    @variable(sp, g_h >= 0)
    @variable(sp, s >= 0)
    @constraint(sp, g_h + g_t == 150)
    c = [50, 100, 150]
    @stageobjective(sp, c[t] * g_t)
    # =========================================================================
    # New stuff below Here
    # Add inflow as a state
    @variable(sp, inflow[1:2], SDDP.State, initial_value = 50.0)
    # Add the random variable as a control variable
    @variable(sp, ε[1:2])

    #why inflow[1:2] and ε[1:2] is written in this way and not inflow[2:1] and ε[2:1]?
    # Needed to transpose? Seems to me that the problem set up is 2x1 


This is a Julia thing. 2:1 means "the range from 2 to 1 in steps of +1", which is an empty set:

julia> x = 2:1
2:1

julia> collect(x)
Int64[]
 
You could write instead

julia> y = 2:-1:1
2:-1:1

julia> collect(y)
2-element Vector{Int64}:
 2
 1

which is "the range from 2 to 1 in steps of -1" and this is the set you expect. Also, inflow and ε are vectors. So they only have one dimension. That is, they aren't 2x1 matrices, but length-2 vectors:

julia> x = [1, 2]
2-element Vector{Int64}:
 1
 2

julia> size(x)
(2,)

julia> y = reshape(x, 2, 1)
2×1 Matrix{Int64}:
 1
 2

julia> size(y)
(2, 1)

julia> x == y
false

Conceptually in maths these are same object, but computationally they behave slightly differently. 


    # The equation describing our statistical model
    A = [0.8 0.2; 0.2 0.8]
    inflow_in = [inflow[i].in for i in 1:2]
    inflow_out = [inflow[i].in for i in 1:2]

    #Why does inflow_in and inflow_out have the same structure?

This is a typo! It should be inflow[i].out
 

    @constraint(sp, inflow_out .== A * inflow_in .+ ε)

    #why element by element and not regular matrix addition?

I don't understand the difference.
 
    # Where is b?


You could add inflow_out .== A * inflow_in .+ b .+ ε. 
 
    # Even though in the equation there is a constant b, the equation here and in the entire set up there is no mention of the "intercept" b
    # b is not included in the matrix A as well, that could be one possibility. 

I just didn't write it. You could also think about including it in the ε, not the A matrix.
 

    # The new water balance constraint using the state variable
    @constraint(sp, x.out == x.in - g_h - s + inflow[1].out + inflow[2].out)

    #all variables in the same unit, right? MWH, I think. 

Ahhh. We don't talk about units very much in this example. But yes, they should all have the same units. I guess here it is "cubic meters"? Or perhaps "MWh equivalent volume"?
 
    # Here we have a line equation, not a matrix 

There is only one state, so x is a scalar. It isn't a matrix.
 

    # Assume we have some empirical residuals:
    Ω₁ = [-10.0, 0.1, 9.6]
    Ω₂ = [-10.0, 0.1, 9.6]

    # It's ok to have empirical residuals, but it would not be the case for them to be normal and i.i.d: 
    # E(ε_t)=0, E(ε_t ε_t^T)= \Sigma_ε, and E(ε_t ε_t^T)= 0 for s \neq t

SDDP requires a finite discrete sample space. You could sample a set from the normal distribution, but you can't use the true normal distribution.
 

    Ω = [(ω₁, ω₂) for ω₁ in Ω₁ for ω₂ in Ω₂]
    SDDP.parameterize(sp, Ω) do ω
        JuMP.fix(ε[1], ω[1])
        JuMP.fix(ε[2], ω[2])
        return
    end
end

And, from the objective state, some concerns:

# fuel_cost[t] = ω * fuel_cost[t-1]
# at the beginning of the problem it is mentioned that the problem evolves as auto-regressive process
#However, in the problem "fuel_cost[t] = ω * fuel_cost[t-1]" ω is random
# In general, the error is random and the Coefficient matrix is fixed to give the deterministic part of the forecast
# What am I missing? 

You can think of this as a log-normal AR(1) model.

log(fuel_cost[t]) = log(w * fuel_cost[t-1])
log(fuel_cost[t]) = log(fuel_cost[t-1]) + log(w)

So now the errors are additive in the log space. It's just a modeling trick.
 

model = SDDP.LinearPolicyGraph(
    stages = 3,
    sense = :Min,
    lower_bound = 0.0,
    optimizer = GLPK.Optimizer,
) do subproblem, t
    @variable(subproblem, 0 <= volume <= 200, SDDP.State, initial_value = 200)
    @variables(subproblem, begin
        thermal_generation >= 0
        hydro_generation >= 0
        hydro_spill >= 0
        inflow
    end)

    # inflow is defined as variable, but is missing the inequality companion (no big deal)

This is a free variable, with bounds of (-inf, inf).
 

    @constraints(
        subproblem,
        begin
            volume.out == volume.in + inflow - hydro_generation - hydro_spill
            demand_constraint, thermal_generation + hydro_generation == 150.0
        end
    )

    SDDP.add_objective_state(
        subproblem,
        initial_value = (50.0, 50.0),
        lipschitz = (10_000.0, 10_000.0),
        lower_bound = (50.0, 50.0),
        upper_bound = (150.0, 150.0),
    ) do fuel_cost, ω
        fuel_cost′ = fuel_cost[1] + 0.5 * (fuel_cost[1] - fuel_cost[2]) + ω.fuel
        return (fuel_cost′, fuel_cost[1])
    end

    #what does the prime (') mean in fuel_cost′? (transpose or outgoing variable)

It's just a name. It doesn't have a special meaning.
 
    #why 0.5 * (fuel_cost[1] - fuel_cost[2])? 0.5 for average?

It's just a weight. Ideally you would fit this parameter to data.
 
    # What is the real meaning of ω.fuel? the randomness applied to fuel?

ω is a realization from our sample space which we define below. So ω.fuel is the "f" component in the next line
 

    Ω = [
        (fuel = f, inflow = w) for f in [-10.0, -5.0, 5.0, 10.0] for
        w in [0.0, 50.0, 100.0]
    ]

    # why fuel can be negative? (negative price, loss)

It's not the fuel price, it's the ω in the stochastic process for fuel.
 

    SDDP.parameterize(subproblem, Ω) do ω
        (fuel_cost, fuel_cost_old) = SDDP.objective_state(subproblem)
        @stageobjective(subproblem, fuel_cost * thermal_generation)
        return JuMP.fix(inflow, ω.inflow)
    end
end

    #why fule_cost, and fuel_cost_old (what is the role of "old", same as .in and .out?)

fuel_cost is fuel_cost[t]
fued_cost_old is fuel_cost[t-1]
 
    #So, I need to understand the prime, otherwise, I'm missing something. 

I don't understand the queston.
 
    # Here we have (inflow, ω.inflow), how to conciliate ω.inflow and ω.fuel? meaning?

ω is a realization from our sample space. This fixes inflow to the value ω.inflow
 
    
    #What parameterize and JuMP.fix actually do? 

Parameterize updates a JuMP model with the realization.
fix sets the lower and upper bounds on the variable to a value