Consider two dams in the same valley. These two dams are connected together, as we suppose that the water turbined by the first dam is an input of the second dam. The goal of this notebook is to show how to use the Stochastic Dual Dynamic Programming (SDDP) algorithm to find how to manage optimally these two dams.
Assume that we have a water inflow $W_t$ which arrives in the first dam between $t$ and $t+1$. This dam turbines a quantity $U_t^1$ of water, and spill a quantity $S_t^1$.
Its dynamic is:
$$ X_{t+1}^1 = X_{t}^1 + W^t - U_t^1 - S_t^1 $$The turbined flow arrives in the second dam, which turbines a quantity $U_t^2$ and spills a quantity $S_t^2$. So its dynamic is:
$$ X_{t+1}^2 = X_{t}^2 + U_t^1 + S_t^1 - U_t^2 - S_t^2 $$Thus, we could define the state: $$ X_t = (X_t^1, X_t^2)$$
and the control: $$U_t = (U_t^1, U_t^2, S_t^1, S_t^2) $$
The two turbines produce a quantity of electricity proportionnal to the flow turbined, and this electricity is sold into the market at a price $c_t$. So we gain at each timestep:
$$ C_t(X_t, U_t) = c_t (U_t^1 + U_t^2) $$Here, we suppose that costs are negative, as we sell electricity onto the network.
We want to maximize our expected gains, so we minimize the following quantity:
$$ J = \mathbb{E} \left[ \sum_{i=1}^{T_f} C_t(X_t, U_t) \right]$$First, we need to import some modules:
using JuMP, Clp, StochDynamicProgramming, PyPlot
INFO: Recompiling stale cache file /home/fpacaud/.julia/lib/v0.5/PyPlot.ji for module PyPlot.
JuMP is the julia module for Mathematical Programming, Clp the module calling a linear solver (can be replaced by CPLEX or Gurobi), StochDynamicProgramming is a SDDP module in Julia, PyPlot is used here to plot the results
Then, we define the constants of this problem:
# Number of timesteps (as we manage the dams over a year, it is equal to the number of weeks):
TF = 52
# Capacity of dams:
VOLUME_MAX = 100
VOLUME_MIN = 0
# Specify the maximum flow that could be turnined:
CONTROL_MAX = round(Int, .4/7. * VOLUME_MAX) + 1
CONTROL_MIN = 0
# Some statistics about aleas (water inflow):
W_MAX = round(Int, .5/7. * VOLUME_MAX)
W_MIN = 0
DW = 1
T0 = 1
# Define aleas' space:
N_ALEAS = Int(round(Int, (W_MAX - W_MIN) / DW + 1))
ALEAS = linspace(W_MIN, W_MAX, N_ALEAS);
Now, we generate a random process to simulate the evolution of electricity prices over a year:
COST = -66*2.7*(1 + .5*(rand(TF) - .5));
We could plot the evolution of prices with matplotlib:
plot(COST, color="k", lw=2)
xlim(0, 53)
xlabel("Week")
ylabel("Cost")
PyObject <matplotlib.text.Text object at 0x7f1aa19ef710>
We could now define the dynamic of our system:
# Define dynamic of the dam:
function dynamic(t, x, u, w)
return [x[1] - u[1] - u[3] + w[1], x[2] - u[2] - u[4] + u[1] + u[3]]
end
dynamic (generic function with 1 method)
and the cost at time $t$:
# Define cost corresponding to each timestep:
function cost_t(t, x, u, w)
return COST[t] * (u[1] + u[2])
end
cost_t (generic function with 1 method)
Now, we build a function to simulate the evolution of water inflow over one year. We aim to have less water in summer than in winter.
"""Build aleas probabilities for each month."""
function build_aleas()
aleas = zeros(N_ALEAS, TF)
# take into account seasonality effects:
unorm_prob = linspace(1, N_ALEAS, N_ALEAS)
proba1 = unorm_prob / sum(unorm_prob)
proba2 = proba1[N_ALEAS:-1:1]
for t in 1:TF
aleas[:, t] = (1 - sin(pi*t/TF)) * proba1 + sin(pi*t/TF) * proba2
end
return aleas
end
"""Build an admissible scenario for water inflow."""
function build_scenarios(n_scenarios::Int64, probabilities)
scenarios = zeros(n_scenarios, TF)
for scen in 1:n_scenarios
for t in 1:TF
Pcum = cumsum(probabilities[:, t])
n_random = rand()
prob = findfirst(x -> x > n_random, Pcum)
scenarios[scen, t] = prob
end
end
return scenarios
end
build_scenarios
We could test our generator with one scenario:
scenario = build_scenarios(1, build_aleas());
plot(scenario', lw=2, color="blue")
xlabel("Week")
ylabel("Water inflow")
xlim(0, 53)
(0,53)
To use these scenarios in SDDP, we must use a discrete distribution for each timestep. The following function generates n_scenarios
and returns a vector of NoiseLaw
corresponding to the evolution of aleas distribution along the time:
"""Build probability distribution at each timestep.
Return a Vector{NoiseLaw}"""
function generate_probability_laws(n_scenarios)
aleas = build_scenarios(n_scenarios, build_aleas())
laws = Vector{NoiseLaw}(TF-1)
# uniform probabilities:
proba = 1/n_scenarios*ones(n_scenarios)
for t=1:TF-1
laws[t] = NoiseLaw(aleas[:, t], proba)
end
return laws
end
generate_probability_laws
We generate 10 scenarios and fit a probability distribution at each timestep:
N_SCENARIO = 10
aleas = generate_probability_laws(10);
We define the bounds over the state and the control:
x_bounds = [(VOLUME_MIN, VOLUME_MAX), (VOLUME_MIN, VOLUME_MAX)];
u_bounds = [(CONTROL_MIN, CONTROL_MAX), (CONTROL_MIN, CONTROL_MAX), (0, 100), (0, 100)];
and the initial position $X_0$:
x0 = [50, 50]
2-element Array{Int64,1}: 50 50
We build an instance of LinearDynamicLinearCostSPModel
to translate our problem in SDDP:
model = LinearSPModel(TF, # number of timestep
u_bounds, # control bounds
x0, # initial state
cost_t, # cost function
dynamic, # dynamic function
aleas);
We add bounds to state:
set_state_bounds(model, x_bounds)
2-element Array{Tuple{Int64,Int64},1}: (0,100) (0,100)
We define the parameters of the algorithm:
# LP solver:
solver = ClpSolver()
# Set tolerance at 0.1% for convergence:
EPSILON = 0.001
# Maximum iterations:
MAX_ITER = 20
params = SDDPparameters(solver,
passnumber=10,
gap=EPSILON,
compute_ub=5,
montecarlo_in_iter=100,
max_iterations=MAX_ITER);
# this function build a SDDPInterface object, and run SDDP onto it. Then, it renders `sddp`
sddp = @time solve_SDDP(model, params, 1);
SDDP Interface initialized Initialize cuts Starting SDDP iterations Compute upper-bound with 100 scenarios... Compute upper-bound with 100 scenarios... Compute upper-bound with 100 scenarios... Compute upper-bound with 100 scenarios... Compute final upper-bound with 1000 scenarios... ############################################################ SDDP CONVERGENCE - Exact lower bound: -1.0147e+05 [Gap < -0.49%] - Estimation of upper-bound: -1.0126e+05 - Upper-bound's s.t.d: 4.6818e+03 - Confidence interval (95%): [-1.0155e+05 , -1.0097e+05] ############################################################ 70.956267 seconds (16.29 M allocations: 1.676 GB, 0.79% gc time)
# We plot evolution of lower-bound:
plot(sddp.stats.upper_bounds, lw=2, c="darkred", lw=2)
plot(sddp.stats.lower_bounds, lw=2, c="darkblue", lw=2)
grid()
xlabel("Number of iterations")
ylabel("Lower bound");
# and evolution of execution time:
plot(sddp.stats.exectime, lw=2, c="k")
grid()
xlabel("Number of iterations")
ylabel("Time (s)");
The algorithm returns the bellman functions (V
) and a vector of JuMP.Model
used to approximate these functions with linear cuts.
We suppose given a scenario of inflows:
alea_year = Array([7.0 7.0 8.0 3.0 1.0 1.0 3.0 4.0 3.0 2.0 6.0 5.0 2.0 6.0 4.0 7.0 3.0 4.0 1.0 1.0 6.0 2.0 2.0 8.0 3.0 7.0 3.0 1.0 4.0 2.0 4.0 1.0 3.0 2.0 8.0 1.0 5.0 5.0 2.0 1.0 6.0 7.0 5.0 1.0 7.0 7.0 7.0 4.0 3.0 2.0 8.0 7.0])
1×52 Array{Float64,2}: 7.0 7.0 8.0 3.0 1.0 1.0 3.0 4.0 … 7.0 7.0 4.0 3.0 2.0 8.0 7.0
We store this scenario as a 3D array, so it could be used to compute a forward-pass:
aleas = zeros(52, 1, 1)
aleas[:, 1, 1] = alea_year;
We have only one scenario, so we set the forwardPassNumber equal to 1:
Find the optimal control with a forward simulation:
costs, stocks = StochDynamicProgramming.simulate(sddp, 100);
SDDP_COST = costs[1]
println("SDDP cost: ", SDDP_COST)
SDDP cost: -99495.47365243685
And the optimal solution is:
plot(stocks[:, :, 1])
plot(stocks[:, :, 2])
xlabel("Week")
ylabel("Stocks")
grid()
ylim(0, 100)
xlim(0, 53);
To check the results given by SDDP, we solve the deterministic problem with JuMP:
m = Model(solver=solver)
@variable(m, VOLUME_MIN <= x1[1:(TF+1)] <= VOLUME_MAX)
@variable(m, VOLUME_MIN <= x2[1:(TF+1)] <= VOLUME_MAX)
@variable(m, CONTROL_MIN <= u1[1:TF] <= CONTROL_MAX)
@variable(m, CONTROL_MIN <= u2[1:TF] <= CONTROL_MAX)
@variable(m, u3[1:TF] >= 0)
@variable(m, u4[1:TF] >= 0)
@objective(m, Min, sum(COST[i]*(u1[i] + u2[i]) for i = 1:TF))
for i in 1:TF
@constraint(m, x1[i+1] - x1[i] + u1[i] + u3[i] - alea_year[i] == 0)
@constraint(m, x2[i+1] - x2[i] + u2[i] + u4[i] - u1[i] - u3[i] == 0)
end
@constraint(m, x1[1] == x0[1])
@constraint(m, x2[1] == x0[2])
status = solve(m)
println(status)
LP_COST = getobjectivevalue(m)
println("LP value: ", LP_COST)
Optimal LP value: -104297.05688368101
And we plot the evolution of the stocks:
plot(getvalue(x1))
plot(getvalue(x2))
xlabel("Week")
ylabel("Stocks")
grid()
ylim(0, 100)
xlim(0, 53);
The solution given by the solver is more optimistic, as it assumes that the future is known in advance.
If we consider the costs, we have a discrepancy between the two solutions :
abs((LP_COST - SDDP_COST)/LP_COST)
0.046037571669919765