Solve the hydro-thermal power system planning problem: periodical SDDP¶
The hydro-thermal power system planning problem is periodical with a period of 12. In this tutorial, we use periodical SDDP to solve the problem for infinite horizon.
The syntax of periodical SDDP is very similar to the classical SDDP. We will only highlight the differences. More details can be found at http://www.optimization-online.org/DB_FILE/2019/09/7367.pdf
[1]:
import pandas
import numpy
import gurobipy
from msppy.msp import MSLP
from msppy.solver import PSDDP
from msppy.evaluation import EvaluationTrue, Evaluation
import sys
gamma = numpy.array(pandas.read_csv(
"./data/gamma.csv",
names=[0,1,2,3],
index_col=0,
skiprows=1,
))
sigma = [
numpy.array(pandas.read_csv(
"./data/sigma_{}.csv".format(i),
names=[0,1,2,3],
index_col=0,
skiprows=1,
)) for i in range(12)
]
exp_mu = numpy.array(pandas.read_csv(
"./data/exp_mu.csv",
names=[0,1,2,3],
index_col=0,
skiprows=1,
))
hydro_ = pandas.read_csv("./data/hydro.csv", index_col=0)
demand = pandas.read_csv("./data/demand.csv", index_col=0)
deficit_ = pandas.read_csv("./data/deficit.csv", index_col=0)
exchange_ub = pandas.read_csv("./data/exchange.csv", index_col=0)
exchange_cost = pandas.read_csv("./data/exchange_cost.csv", index_col=0)
thermal_ = [pandas.read_csv("./data/thermal_{}.csv".format(i),
index_col=0) for i in range(4)]
stored_initial = hydro_['INITIAL'][:4]
inflow_initial = hydro_['INITIAL'][4:8]
def sampler(t):
def inner(random_state):
noise = numpy.exp(
random_state.multivariate_normal(mean=[0]*4, cov=sigma[t%12]))
coef = [None]*4
rhs = [None]*4
for i in range(4):
coef[i] = -noise[i]*gamma[t%12][i]*exp_mu[t%12][i]/exp_mu[(t-1)%12][i]
rhs[i] = noise[i]*(1-gamma[t%12][i])*exp_mu[t%12][i]
return (coef+rhs)
return inner
Build the true problem and make discretization¶
[2]:
HydroThermal = MSLP(T=13, bound=0, discount=0.9906)
Academic license - for non-commercial use only
Academic license - for non-commercial use only
Academic license - for non-commercial use only
Academic license - for non-commercial use only
Academic license - for non-commercial use only
Academic license - for non-commercial use only
Academic license - for non-commercial use only
Academic license - for non-commercial use only
Academic license - for non-commercial use only
Academic license - for non-commercial use only
Academic license - for non-commercial use only
Academic license - for non-commercial use only
Academic license - for non-commercial use only
Periodical SDDP algorithm solves the problem for a single period plus an initial stage. In this case, the number of stages to consider is 13, set by \(T=13\).
[3]:
for t in range(13):
m = HydroThermal[t]
stored_now,stored_past = m.addStateVars(4, ub=hydro_['UB'][:4], name="stored")
inflow_now,inflow_past = m.addStateVars(4, name="inflow")
spill = m.addVars(4, obj=0.001, name="spill")
hydro = m.addVars(4, ub=hydro_['UB'][-4:], name="hydro")
deficit = m.addVars(
[(i,j) for i in range(4) for j in range(4)],
ub = [
demand.iloc[t%12][i] * deficit_['DEPTH'][j]
for i in range(4) for j in range(4)
],
obj = [
deficit_['OBJ'][j]
for i in range(4) for j in range(4)
],
name = "deficit")
thermal = [None] * 4
for i in range(4):
thermal[i] = m.addVars(
len(thermal_[i]),
ub=thermal_[i]['UB'],
lb=thermal_[i]['LB'],
obj=thermal_[i]['OBJ'],
name="thermal_{}".format(i)
)
exchange = m.addVars(5,5, obj=exchange_cost.values.flatten(),
ub=exchange_ub.values.flatten(), name="exchange")
thermal_sum = m.addVars(4, name="thermal_sum")
m.addConstrs(thermal_sum[i] ==
gurobipy.quicksum(thermal[i].values()) for i in range(4))
for i in range(4):
m.addConstr(
thermal_sum[i]
+ gurobipy.quicksum(deficit[(i,j)] for j in range(4))
+ hydro[i]
- gurobipy.quicksum(exchange[(i,j)] for j in range(5))
+ gurobipy.quicksum(exchange[(j,i)] for j in range(5))
== demand.iloc[t%12][i],
name = 'demand',
)
m.addConstr(
gurobipy.quicksum(exchange[(j,4)] for j in range(5))
- gurobipy.quicksum(exchange[(4,j)] for j in range(5))
== 0
)
m.addConstrs(
stored_now[i] + spill[i] + hydro[i] - stored_past[i] == inflow_now[i]
for i in range(4)
)
if t == 0:
m.addConstrs(stored_past[i] == stored_initial[i] for i in range(4))
m.addConstrs(inflow_now[i] == inflow_initial[i] for i in range(4))
else:
TS = m.addConstrs(inflow_now[i] + inflow_past[i] == 0 for i in range(4))
m.add_continuous_uncertainty(
uncertainty=sampler(t),
locations=(
[(TS[i],inflow_past[i]) for i in range(4)]
+ [TS[i] for i in range(4)]
),
)
HydroThermal.discretize(n_samples=100, random_state=888)
Solve the problem¶
We now call PSDDP solver to run the periodical SDDP for 20 iterations.
Backward passes of the periodical SDDP generates cuts for the first 13 stages.
Forward passes of the periodical SDDP generates trial points. Trial points can be just obtained from solving these 13 stages. They can also be obtained from later stages (since the problem is periodical). It is often found solving more stages makes the algorithm converge faster. Here we set \(\textrm{forward_T}=120\), meaning that trial points are generated from the first 120 stages.
[4]:
HT_psddp = PSDDP(HydroThermal)
HT_psddp.solve(max_iterations=10, forward_T=120)
----------------------------------------------------------------
SDDP Solver, Lingquan Ding
----------------------------------------------------------------
Iteration Bound Value Time
----------------------------------------------------------------
1 245082.919600 1734139943.585172 0.426226
2 275523.302075 396460564.856688 0.773433
3 3605638.658380 343399590.174994 0.456558
4 7472497.857329 591001246.721298 0.357867
5 7472497.857329 175229142.726806 0.365556
6 8608579.050937 455731628.093480 0.337117
7 18568417.721664 316499896.612110 0.328819
8 18568417.721664 252349984.539381 0.382036
9 18568417.721664 1034789138.913040 0.362973
10 29642949.251290 483878406.759454 0.357292
----------------------------------------------------------------
Time: 4.147876501083374 seconds
Algorithm stops since iteration:10 has reached
Evaluate the obtained policy¶
The obtained policy is implementable feasible for any finite number of stages. We can for example, set \(\textrm{query}_T=60}\) as below to evaluate the policy for the first 60 stages.
[5]:
result = Evaluation(HydroThermal)
result.run(n_simulations=10, query_T=60)
result.CI
[5]:
(199193170.50824797, 406535057.5869812)