SDDP.jl
SDDP.jl copied to clipboard
Documentation feedback
Some feedback from a user in my emails (search "Suggestions for JuMP tutorial" 28 February).
https://odow.github.io/SDDP.jl/stable/tutorial/basic/07_arma/ https://odow.github.io/SDDP.jl/stable/tutorial/advanced/11_objective_states/
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