Hydro-thermal scheduling
Problem Description
In a hydro-thermal problem, the agent controls a hydro-electric generator and reservoir. Each time period, they need to choose a generation quantity from thermal g_t, and hydro g_h, in order to meet demand w_d, which is a stagewise-independent random variable. The state variable, x, is the quantity of water in the reservoir at the start of each time period, and it has a minimum level of 5 units and a maximum level of 15 units. We assume that there are 10 units of water in the reservoir at the start of time, so that x_0 = 10. The state-variable is connected through time by the water balance constraint: x.out = x.in - g_h - s + w_i, where x.out is the quantity of water at the end of the time period, x.in is the quantity of water at the start of the time period, s is the quantity of water spilled from the reservoir, and w_i is a stagewise-independent random variable that represents the inflow into the reservoir during the time period.
We assume that there are three stages, t=1, 2, 3, representing summer-fall, winter, and spring, and that we are solving this problem in an infinite-horizon setting with a discount factor of 0.95.
In each stage, the agent incurs the cost of spillage, plus the cost of thermal generation. We assume that the cost of thermal generation is dependent on the stage t = 1, 2, 3, and that in each stage, w is drawn from the set (w_i, w_d) = {(0, 7.5), (3, 5), (10, 2.5)} with equal probability.
Importing packages
For this example, in addition to SDDP, we need HiGHS as a solver and Statisitics to compute the mean of our simulations.
using HiGHS
using SDDP
using StatisticsConstructing the policy graph
There are three stages in our problem, so we construct a linear policy graph with three stages using SDDP.LinearGraph:
graph = SDDP.LinearGraph(3)Root
0
Nodes
1
2
3
Arcs
0 => 1 w.p. 1.0
1 => 2 w.p. 1.0
2 => 3 w.p. 1.0
Then, because we want to solve an infinite-horizon problem, we add an additional edge between node 3 and node 1 with probability 0.95:
SDDP.add_edge(graph, 3 => 1, 0.95)Constructing the model
Much of the macro code (i.e., lines starting with @) in the first part of the following should be familiar to users of JuMP.
Inside the do-end block, sp is a standard JuMP model, and t is an index for the state variable that will be called with t = 1, 2, 3.
The state variable x, constructed by passing the SDDP.State tag to @variable is actually a Julia struct with two fields: x.in and x.out corresponding to the incoming and outgoing state variables respectively. Both x.in and x.out are standard JuMP variables. The initial_value keyword provides the value of the state variable in the root node (i.e., x_0).
Compared to a JuMP model, one key difference is that we use @stageobjective instead of @objective. The SDDP.parameterize function takes a list of supports for w and parameterizes the JuMP model sp by setting the right-hand sides of the appropriate constraints (note how the constraints initially have a right-hand side of 0). By default, it is assumed that the realizations have uniform probability, but a probability mass vector can also be provided.
model = SDDP.PolicyGraph(
graph,
sense = :Min,
lower_bound = 0.0,
optimizer = HiGHS.Optimizer,
) do sp, t
@variable(sp, 5 <= x <= 15, SDDP.State, initial_value = 10)
@variable(sp, g_t >= 0)
@variable(sp, g_h >= 0)
@variable(sp, s >= 0)
@constraint(sp, balance, x.out - x.in + g_h + s == 0)
@constraint(sp, demand, g_h + g_t == 0)
@stageobjective(sp, s + t * g_t)
SDDP.parameterize(sp, [[0, 7.5], [3, 5], [10, 2.5]]) do w
set_normalized_rhs(balance, w[1])
return set_normalized_rhs(demand, w[2])
end
endA policy graph with 3 nodes.
Node indices: 1, 2, 3
Training the policy
Once a model has been constructed, the next step is to train the policy. This can be achieved using SDDP.train. There are many options that can be passed, but iteration_limit terminates the training after the prescribed number of SDDP iterations.
SDDP.train(model, iteration_limit = 100)------------------------------------------------------------------------------
SDDP.jl (c) Oscar Dowson and SDDP.jl contributors, 2017-23
Problem
Nodes : 3
State variables : 1
Scenarios : Inf
Existing cuts : false
Subproblem structure : (min, max)
Variables : (6, 6)
VariableRef in MOI.LessThan{Float64} : (1, 1)
VariableRef in MOI.GreaterThan{Float64} : (5, 5)
AffExpr in MOI.EqualTo{Float64} : (2, 2)
Options
Solver : serial mode
Risk measure : SDDP.Expectation()
Sampling scheme : SDDP.InSampleMonteCarlo
Numerical stability report
Non-zero Matrix range [1e+00, 1e+00]
Non-zero Objective range [1e+00, 3e+00]
Non-zero Bounds range [5e+00, 2e+01]
Non-zero RHS range [2e+00, 1e+01]
No problems detected
Iteration Simulation Bound Time (s) Proc. ID # Solves
1 1.940000e+02 8.617986e+01 9.489083e-02 1 267
2 8.850000e+01 1.066277e+02 1.105609e-01 1 402
3 1.760000e+02 1.299277e+02 1.275361e-01 1 549
4 3.540083e+02 1.546951e+02 1.611218e-01 1 816
5 2.685000e+02 1.730631e+02 1.975241e-01 1 1107
6 3.385817e+02 1.847377e+02 2.467339e-01 1 1470
7 2.776220e+01 1.879896e+02 2.558401e-01 1 1533
8 1.321895e+02 1.925415e+02 2.697918e-01 1 1632
9 2.205322e+02 2.010833e+02 3.010809e-01 1 1851
10 1.996429e+02 2.046703e+02 3.282509e-01 1 2034
11 8.023055e+01 2.070885e+02 3.410950e-01 1 2121
12 6.942059e+01 2.084161e+02 3.507810e-01 1 2184
13 5.672949e+02 2.187836e+02 4.360650e-01 1 2775
14 1.742083e+02 2.204837e+02 4.682219e-01 1 2982
15 8.638027e+02 2.275930e+02 6.091249e-01 1 3873
16 2.500000e+00 2.277098e+02 6.117699e-01 1 3888
17 2.893484e+02 2.285684e+02 6.591361e-01 1 4167
18 3.064593e+02 2.300139e+02 7.093339e-01 1 4458
19 7.853518e+01 2.302053e+02 7.285879e-01 1 4569
20 7.974752e+02 2.321294e+02 8.652229e-01 1 5352
21 2.088776e+02 2.322546e+02 8.928258e-01 1 5511
22 7.597159e+02 2.337308e+02 1.036667e+00 1 6294
23 2.371941e+02 2.339254e+02 1.072508e+00 1 6489
24 2.834876e+01 2.339254e+02 1.074971e+00 1 6504
25 9.054597e+02 2.349435e+02 1.317811e+00 1 7491
26 3.066430e+01 2.349839e+02 1.327241e+00 1 7542
27 6.414651e+02 2.353374e+02 1.462371e+00 1 8265
28 4.871049e+02 2.354549e+02 1.544904e+00 1 8688
29 1.047948e+02 2.355252e+02 1.568495e+00 1 8811
30 1.687157e+02 2.355556e+02 1.607776e+00 1 8994
31 7.261254e+02 2.357417e+02 1.775061e+00 1 9789
32 4.492668e+01 2.357959e+02 1.789669e+00 1 9852
33 1.026583e+02 2.357982e+02 1.814028e+00 1 9963
34 5.574223e+02 2.359182e+02 1.936833e+00 1 10530
35 2.729798e+02 2.359679e+02 2.000950e+00 1 10821
36 3.670739e+02 2.359917e+02 2.105078e+00 1 11268
37 1.064541e+02 2.360178e+02 2.132552e+00 1 11391
38 2.668830e+02 2.360337e+02 2.202269e+00 1 11682
39 2.650000e+01 2.360346e+02 2.214082e+00 1 11733
40 1.905546e+01 2.360598e+02 2.223415e+00 1 11772
41 3.435508e+02 2.360827e+02 2.312203e+00 1 12159
42 1.020289e+02 2.361113e+02 2.333026e+00 1 12246
43 1.443556e+02 2.361180e+02 2.366890e+00 1 12393
44 5.542970e+02 2.361749e+02 2.489658e+00 1 12900
45 2.100219e+02 2.361927e+02 2.549145e+00 1 13131
46 1.785834e+02 2.362168e+02 2.597539e+00 1 13314
47 1.013430e+02 2.362230e+02 2.621315e+00 1 13401
48 7.203920e+01 2.362319e+02 2.647166e+00 1 13500
49 3.640509e+02 2.362434e+02 2.751302e+00 1 13875
50 2.757174e+02 2.362652e+02 2.820198e+00 1 14130
51 3.457205e+02 2.362718e+02 2.922911e+00 1 14493
52 3.015133e+02 2.363003e+02 3.028077e+00 1 14868
53 1.157184e+02 2.363062e+02 3.077174e+00 1 15051
54 8.404997e+01 2.363069e+02 3.101297e+00 1 15138
55 2.530183e+02 2.363173e+02 3.154777e+00 1 15345
56 5.651378e+01 2.363207e+02 3.170855e+00 1 15408
57 8.599815e+01 2.363252e+02 3.201237e+00 1 15519
58 3.007672e+02 2.363312e+02 3.297075e+00 1 15870
59 1.914912e+02 2.363317e+02 3.342249e+00 1 16041
60 3.203913e+01 2.363360e+02 3.356054e+00 1 16092
61 6.600653e+01 2.363382e+02 3.375784e+00 1 16167
62 2.461688e+02 2.363422e+02 3.477498e+00 1 16506
63 1.850595e+01 2.363446e+02 3.481664e+00 1 16521
64 2.070102e+02 2.363465e+02 3.549432e+00 1 16752
65 2.005412e+00 2.363480e+02 3.553320e+00 1 16767
66 1.720059e+02 2.363534e+02 3.613916e+00 1 16986
67 1.579869e+02 2.363557e+02 3.672982e+00 1 17193
68 2.460232e+02 2.363641e+02 3.736540e+00 1 17424
69 2.234751e+02 2.363649e+02 3.787713e+00 1 17583
70 1.099001e+01 2.363678e+02 3.800498e+00 1 17622
71 2.329982e+02 2.363699e+02 3.870344e+00 1 17841
72 3.005747e+02 2.363764e+02 3.954277e+00 1 18120
73 7.794549e+01 2.363778e+02 3.974398e+00 1 18183
74 2.720713e+02 2.363821e+02 4.061216e+00 1 18450
75 1.035026e+02 2.363821e+02 4.076381e+00 1 18501
76 7.799956e+01 2.363825e+02 4.087912e+00 1 18540
77 4.301657e+01 2.363834e+02 4.103778e+00 1 18591
78 5.850102e+01 2.363839e+02 4.122771e+00 1 18654
79 9.099923e+01 2.363850e+02 4.165399e+00 1 18789
80 3.451472e+01 2.363875e+02 4.181850e+00 1 18840
81 7.499980e+01 2.363885e+02 4.209658e+00 1 18927
82 2.925173e+02 2.363931e+02 4.320508e+00 1 19266
83 7.548398e+01 2.363935e+02 4.348609e+00 1 19353
84 2.500000e+01 2.363938e+02 4.356973e+00 1 19380
85 3.225046e+02 2.363970e+02 4.422204e+00 1 19599
86 3.505412e+02 2.363999e+02 4.551562e+00 1 20010
87 4.190096e+02 2.364045e+02 4.677771e+00 1 20409
88 4.355059e+01 2.364052e+02 4.698049e+00 1 20472
89 1.530011e+02 2.364072e+02 4.752772e+00 1 20643
90 2.534959e+02 2.364083e+02 4.834090e+00 1 20886
91 3.390081e+02 2.364109e+02 4.933929e+00 1 21189
92 1.552565e+01 2.364110e+02 4.942851e+00 1 21216
93 1.275182e+02 2.364112e+02 5.004943e+00 1 21399
94 2.105008e+02 2.364128e+02 5.077047e+00 1 21606
95 8.435292e+02 2.364143e+02 5.330378e+00 1 22317
96 4.405284e+02 2.364174e+02 5.491970e+00 1 22752
97 1.535465e+02 2.364176e+02 5.541896e+00 1 22887
98 9.850103e+01 2.364177e+02 5.570416e+00 1 22962
99 1.069949e+02 2.364188e+02 5.615007e+00 1 23085
100 8.750134e+01 2.364192e+02 5.641679e+00 1 23160
Terminating training
Status : iteration_limit
Total time (s) : 5.641679e+00
Total solves : 23160
Best bound : 2.364192e+02
Simulation CI : 2.247979e+02 ± 4.010270e+01
------------------------------------------------------------------------------Simulating the policy
After training, we can simulate the policy using SDDP.simulate.
sims = SDDP.simulate(model, 100, [:g_t])
mu = round(mean([s[1][:g_t] for s in sims]), digits = 2)
println("On average, $(mu) units of thermal are used in the first stage.")On average, 2.55 units of thermal are used in the first stage.Extracting the water values
Finally, we can use SDDP.ValueFunction and SDDP.evaluate to obtain and evaluate the value function at different points in the state-space. Note that since we are minimizing, the price has a negative sign: each additional unit of water leads to a decrease in the the expected long-run cost.
V = SDDP.ValueFunction(model[1])
cost, price = SDDP.evaluate(V, x = 10)(233.52670428215484, Dict(:x => -0.6430290959743883))This page was generated using Literate.jl.