import math
from dataclasses import dataclass
from plot import draw_prices, draw_schedule
from report import report_schedule
@dataclass
class Period:
price: float
inflow: float
@dataclass
class Turbine:
efficiency: float
min_discharge: float
max_discharge: float
min_production: float
max_production: float
@dataclass
class Plant:
periods: list[Period]
turbines: list[Turbine]
min_volume: float = 5e6
max_volume: float = 30e6
initial_volume: float = 25e6
final_volume: float = 25e6
a_head: float = 1e-7
b_head: float = 50.0
a_tail_q: float = 0.005
a_tail_s: float = 0.001
b_tail: float = 30.0
@property
def num_periods(self):
return len(self.periods)
@property
def num_turbines(self):
return len(self.turbines)
@dataclass
class Schedule:
revenue: float
discharge: list[list[float]]
production: list[list[float]]
volume: list[float]
headwater: list[float]
tailwater: list[float]
net_head: list[float]Hydro Unit Commitment
Schedule the water discharge of several turbines sharing one reservoir to maximize the revenue from selling generated hydropower.
Introduction
Given a set of production units, the unit commitment problem consists in determining the status (on/off) of each unit at each period as well as the quantity of electricity to produce, in order to maximize the profit or satisfy demands while minimizing operational costs.
Here, we focus on hydropower generation.
The powerhouse contains a turbine and a generator. The water goes from the upstream reservoir to the downstream river through the penstock, passing by the powerhouse. The flow of water through the blades of a turbine generates electricity which is then transferred from the generator to its final destination through the long distance power lines. The water that passes through the penstock is used to generate power, but water could also be released from the reservoir directly to the downstream river. This water is said to be spilled, as it is not used for production, but rather to avoid an overflow of the reservoir.
Image source: https://commons.wikimedia.org/wiki/File:Hydroelectric_dam-es.svg
How much electricity the water yields depends on how far it falls before reaching the turbine. That drop is measured between two water levels, the headwater of the upstream reservoir and the tailwater of the downstream river. Their difference is the net head, the effective height of the water column driving the turbine. Both levels move as the plant operates, the reservoir dropping as water leaves it and the river rising as discharged and spilled water arrive, so the net head changes from one period to the next.
The quantity of electricity \(p_{j, t}\) produced by turbine \(j\) at period \(t\) is a nonlinear function of the net head \(h^\text{net}_{t}\) and the water discharge \(q_{j,t}\):
\[ p_{j,t} = \eta_j \cdot q_{j,t} \cdot h^\text{net}_{t} \]
where \(\eta_j\) is called the efficiency of turbine \(j\), a coefficient that also carries the density of water, gravity and the conversion to megawatts.
That product of two variables is what makes the problem nonlinear, and it is also what ties the turbines together. They draw on one reservoir, so they share one net head, and water one of them discharges lowers what the others can produce later.
Problem description
Here, we consider a problem composed of several turbines and one reservoir. The objective is to maximize the revenue from selling generated power.
Input:
\(H\) periods; for each period \(t = 1, \dots, H\)
- a price \(R_t\) for selling a unit of electricity
- an inflow \(\Omega_t\)
\(N\) turbines; for each turbine \(j = 1, \dots, N\):
- a minimum discharge \(Q^\text{min}_j\)
- a maximum discharge \(Q^\text{max}_j\)
- a minimum production \(P^\text{min}_j\)
- a maximum production \(P^\text{max}_j\)
- an efficiency \(\eta_j\)
a minimum volume \(V^\text{min}\) for the reservoir
a maximum volume \(V^\text{max}\) for the reservoir
an initial volume \(V^\text{init}\) for the reservoir
a minimum final volume \(V^\text{final}\) for the reservoir
a function \(F^\text{head}(v) = A^\text{head} v + B^\text{head}\) that returns the headwater depending on the volume in the reservoir
a function \(F^\text{tail}(q, s) = A^\text{tail}_q q + A^\text{tail}_s s + B^\text{tail}\) that returns the tailwater depending on the water discharge and the spilled water
Problem: find the quantity of water to send through each turbine at each period such that:
The volume in the reservoir always stays in \([V^\text{min}, V^\text{max}]\)
The initial volume is equal to \(V^\text{init}\)
The final volume is greater than \(V^\text{final}\)
If a turbine is used
- the water discharge in this turbine belongs to \([Q^\text{min}, Q^\text{max}]\); otherwise, it is null.
- the production of this turbine belongs to \([P^\text{min}, P^\text{max}]\); otherwise, it is null.
Objective: maximize the revenue
Mixed-integer nonlinear program
Variables:
- \(x_{j, t} \in \{ 0, 1 \}\), \(j = 1, \dots, N\), \(t = 1, \dots, H\): \(x_{j, t} = 1\) iff turbine \(j\) is on at period \(t\), otherwise \(0\).
- \(q_{j, t} \in [0, Q_j^\text{max}]\), \(j = 1, \dots, N\), \(t = 1, \dots, H\): water discharge in turbine \(j\) at period \(t\).
- \(q^\text{tot}_t \in \mathbb{R}^+\), \(t = 1, \dots, H\): total water discharge at period \(t\).
- \(s_t \in \mathbb{R}^+\), \(t = 1, \dots, H\): water spilled at period \(t\).
- \(p_{j, t} \in [0, P_j^\text{max}]\), \(j = 1, \dots, N\), \(t = 1, \dots, H\): production of turbine \(j\) at period \(t\)
- \(v_t \in [V^\text{min}, V^\text{max}]\), \(t = 1, \dots, H\): volume in the reservoir at period \(t\)
- \(h^\text{head}_t \in \mathbb{R}\), \(t = 1, \dots, H\): head water at period \(t\)
- \(h^\text{tail}_t \in \mathbb{R}\), \(t = 1, \dots, H\): tail water at period \(t\)
- \(h^\text{net}_t \in \mathbb{R}^+\), \(t = 1, \dots, H\): net head at period \(t\)
Objective: maximize the revenue
\[ \max \sum_{t = 1}^H \sum_{j = 1}^N R_t p_{j, t} \]
Constraints:
- Initial and final volume:
\[ v_1 = V^\text{init} \]
\[ v_H \ge V^\text{final} \]
- Water balance: volume at \(t\) = volume at \(t - 1\) + inflow at \(t\) - consumption at \(t\). A period is one hour and the flows are per second, hence the factor 3600.
\[ \forall t = 2, \dots, H \qquad v_t = v_{t - 1} + 3600 \cdot \Omega_t - 3600 \cdot q^\text{tot}_t - 3600 \cdot s_t \]
- Water discharge limits
\[ \forall j = 1, \dots, N \quad \forall t = 1, \dots, H \qquad q_{j, t} \ge Q_j^\text{min} x_{j, t} \]
\[ \forall j = 1, \dots, N \quad \forall t = 1, \dots, H \qquad q_{j, t} \le Q_j^\text{max} x_{j, t} \]
- Production limits
\[ \forall j = 1, \dots, N \quad \forall t = 1, \dots, H \qquad p_{j, t} \ge P_j^\text{min} x_{j, t} \]
\[ \forall j = 1, \dots, N \quad \forall t = 1, \dots, H \qquad p_{j, t} \le P_j^\text{max} x_{j, t} \]
Both sides are gated by \(x_{j, t}\), so an idle turbine is held at zero discharge and zero production, and a running one is kept between its bounds.
- Total water discharge
\[ \forall t = 1, \dots, H \qquad q_t^\text{tot} = \sum_{j = 1}^N q_{j, t} \]
- Headwater
\[ \forall t = 1, \dots, H \qquad h_t^\text{head} = A^\text{head} v_t + B^\text{head} \]
- Tailwater
\[ \forall t = 1, \dots, H \qquad h_t^\text{tail} = A^\text{tail}_q q^\text{tot}_t + A^\text{tail}_s s_t + B^\text{tail} \]
- Net head
\[ \forall t = 1, \dots, H \qquad h_t^\text{net} = h_t^\text{head} - h_t^\text{tail} \]
- Hydropower production function
\[ \forall j = 1, \dots, N \quad \forall t = 1, \dots, H \qquad p_{j, t} = \eta_j \cdot q_{j,t} \cdot h^\text{net}_{t} \]
The model has the following properties:
- Continuous and binary variables
- Quadratic structures
- Non-convex
- Therefore, it is a non-convex MIQCQP
Input data
A Period carries what the horizon offers at one hour, the price a unit of electricity sells for and the inflow arriving at the reservoir. A Turbine carries its efficiency and the bounds it has to respect while it runs. A Plant holds those two lists plus everything that belongs to the reservoir rather than to any one turbine, namely the volume it must stay within, where it starts and where it has to end, and the coefficients of the two level functions.
A Schedule carries the revenue it earns, the discharge and production of every turbine at every period, and the four series the reservoir follows, namely volume, headwater, tailwater and net head.
The plant is generated rather than read from a file. Ten turbines draw on one reservoir over 168 hourly periods, a week of operation. Prices follow a daily cycle with noise on top, and the inflows are drawn flat. Each language draws with its own generator, so the plant scheduled here is not the one the Julia page schedules, and the two schedules cannot be compared hour by hour.
import random
def generate_plant(num_periods, num_turbines):
draw = random.Random(num_periods * num_turbines)
turbines = []
for _ in range(num_turbines):
min_production = draw.uniform(0, 10)
max_production = min_production + draw.uniform(0, 40)
min_discharge = draw.triangular(0, 100, 10)
max_discharge = draw.triangular(min_discharge, 100, max(min_discharge, 90))
turbines.append(
Turbine(
draw.uniform(0.5, 0.9) * 1e-6 * 1000 * 9.81,
min_discharge,
max_discharge,
min_production,
max_production,
)
)
periods = []
for t in range(num_periods):
inflow = draw.uniform(10, 100) * num_turbines
cycle = 40 * (1 + math.cos(t / 24 * 2 * math.pi)) + 10
periods.append(Period(cycle + draw.uniform(-10, 10), inflow))
return Plant(periods, turbines)
plant = generate_plant(168, 10)Every figure is a stack of panels sharing one time axis, so a period reads straight down from one panel to the next. The discharge panel is a step area per turbine stacked into the total, since a turbine holds its discharge for the whole hour rather than sliding to the next value, and every other panel is a line.
Prices are the input the schedule reacts to, so they come first. The cycle is daily, and the noise on top is what keeps the schedule from being the same day repeated seven times.
draw_prices(plant)Model implementation
import pyomo.environ as pyo
SOLVER_NAME = "knitroampl"maximize_revenue(plant) builds the model, solves it with Knitro, and returns a Schedule. The only nonlinear part is the production function, a product of two variables, so every line of the model is a Constraint written as it reads. The turbine limits are constraints rather than variable bounds, because each side is gated by the on/off variable, so an idle turbine is pinned to zero.
Proving optimality on a nonconvex problem of this size takes far longer than reaching a good schedule, so the branch and bound stops at the first integer feasible point.
def maximize_revenue(plant, *, mip_terminate=1, num_threads=1):
model = pyo.ConcreteModel()
n, h = plant.num_turbines, plant.num_periods
periods, turbines = plant.periods, plant.turbines
model.N = pyo.RangeSet(0, n - 1)
model.H = pyo.RangeSet(0, h - 1)
model.x = pyo.Var(model.N, model.H, within=pyo.Binary)
model.q = pyo.Var(model.N, model.H, within=pyo.NonNegativeReals)
model.p = pyo.Var(model.N, model.H, within=pyo.NonNegativeReals)
model.qtot = pyo.Var(model.H, within=pyo.NonNegativeReals)
model.s = pyo.Var(model.H, within=pyo.NonNegativeReals)
model.v = pyo.Var(model.H, bounds=(plant.min_volume, plant.max_volume))
model.h_head = pyo.Var(model.H)
model.h_tail = pyo.Var(model.H)
model.h_net = pyo.Var(model.H, within=pyo.NonNegativeReals)
model.objective = pyo.Objective(
expr=pyo.quicksum(
periods[t].price * model.p[j, t] for j in model.N for t in model.H
),
sense=pyo.maximize,
)
def water_balance_rule(model, t):
return model.v[t] == model.v[t - 1] + 3600 * (
periods[t].inflow - model.qtot[t] - model.s[t]
)
model.water_balance = pyo.Constraint(
pyo.RangeSet(1, h - 1), rule=water_balance_rule
)
def min_production_rule(model, t, j):
return model.p[j, t] >= turbines[j].min_production * model.x[j, t]
def max_production_rule(model, t, j):
return model.p[j, t] <= turbines[j].max_production * model.x[j, t]
def min_discharge_rule(model, t, j):
return model.q[j, t] >= turbines[j].min_discharge * model.x[j, t]
def max_discharge_rule(model, t, j):
return model.q[j, t] <= turbines[j].max_discharge * model.x[j, t]
model.min_production = pyo.Constraint(model.H, model.N, rule=min_production_rule)
model.max_production = pyo.Constraint(model.H, model.N, rule=max_production_rule)
model.min_discharge = pyo.Constraint(model.H, model.N, rule=min_discharge_rule)
model.max_discharge = pyo.Constraint(model.H, model.N, rule=max_discharge_rule)
def total_discharge_rule(model, t):
return model.qtot[t] == pyo.quicksum(model.q[j, t] for j in model.N)
def headwater_rule(model, t):
return model.h_head[t] == plant.a_head * model.v[t] + plant.b_head
def tailwater_rule(model, t):
return model.h_tail[t] == (
plant.a_tail_q * model.qtot[t] + plant.a_tail_s * model.s[t] + plant.b_tail
)
def net_head_rule(model, t):
return model.h_net[t] == model.h_head[t] - model.h_tail[t]
model.total_discharge = pyo.Constraint(model.H, rule=total_discharge_rule)
model.headwater = pyo.Constraint(model.H, rule=headwater_rule)
model.tailwater = pyo.Constraint(model.H, rule=tailwater_rule)
model.net_head = pyo.Constraint(model.H, rule=net_head_rule)
def production_rule(model, t, j):
return model.p[j, t] == turbines[j].efficiency * model.q[j, t] * model.h_net[t]
model.production = pyo.Constraint(model.H, model.N, rule=production_rule)
model.initial_volume = pyo.Constraint(expr=model.v[0] == plant.initial_volume)
model.final_volume = pyo.Constraint(expr=model.v[h - 1] >= plant.final_volume)
model.initial_spill = pyo.Constraint(expr=model.s[0] == 0)
def initial_discharge_rule(model, j):
return model.q[j, 0] == 0
model.initial_discharge = pyo.Constraint(model.N, rule=initial_discharge_rule)
solver = pyo.SolverFactory(SOLVER_NAME)
solver.options["mip_terminate"] = mip_terminate
solver.options["numthreads"] = num_threads
solver.solve(model, tee=True)
revenue = model.objective()
return Schedule(
revenue,
[[model.q[j, t]() for t in model.H] for j in model.N],
[[model.p[j, t]() for t in model.H] for j in model.N],
[model.v[t]() for t in model.H],
[model.h_head[t]() for t in model.H],
[model.h_tail[t]() for t in model.H],
[model.h_net[t]() for t in model.H],
)schedule = maximize_revenue(plant)Artelys Knitro 16.0.0: mip_terminate=1
numthreads=1
=======================================
Commercial License
Artelys Knitro 16.0.0
=======================================
Knitro changing mip_method from AUTO to 1.
No start point provided -- Knitro computing one.
Knitro presolve eliminated 537 variables (9%) and 568 constraints (6%) in 0.02s.
concurrent_evals 0
datacheck 0
feastol 1e-06
feastol_abs 1e-06
findiff_numthreads 1
hessian_no_f 1
hessopt 1
mip_terminate 1
numthreads 1
opttol 1e-06
opttol_abs 0.001
Knitro changing mip_root_nlpalg from AUTO to 1.
Knitro changing mip_node_nlpalg from AUTO to 1.
Knitro changing mip_branchrule from AUTO to 2.
Knitro changing mip_selectrule from AUTO to 2.
Knitro changing mip_mir from AUTO to 2.
Knitro changing mip_clique from AUTO to 0.
Knitro changing mip_zerohalf from AUTO to 0.
Knitro changing mip_liftproject from AUTO to 0.
Knitro changing mip_knapsack from AUTO to 2.
Knitro changing mip_gomory from AUTO to 0.
Knitro changing mip_cut_flowcover from AUTO to 2.
Knitro changing mip_cut_probing from AUTO to 1.
Knitro changing mip_rounding from AUTO to 3.
Knitro changing mip_heuristic_strategy from AUTO to 1.
Knitro changing mip_heuristic_feaspump from AUTO to 1.
Knitro changing mip_heuristic_misqp from AUTO to 0.
Knitro changing mip_heuristic_mpec from AUTO to 1.
Knitro changing mip_heuristic_diving from AUTO to 1926.
Knitro changing mip_heuristic_fixpropagate from AUTO to 62.
Knitro changing mip_heuristic_lns from AUTO to 0.
Knitro changing mip_heuristic_localsearch from AUTO to 1.
Knitro changing mip_pseudoinit from AUTO to 1.
WARNING: Problem appears to have nonlinear equalities and be non-convex.
The Knitro mixed integer solver is designed for convex problems.
For non-convex problems it is only a heuristic, and the reported
bounds and optimality claims cannot be verified.
Problem Characteristics | Presolved
-----------------------
Problem type: MIQCQP
Objective: maximize / linear
Number of variables: 6048 | 5511
bounds: lower upper range | lower upper range
3864 0 1848 | 3673 0 1836
free fixed | free fixed
336 0 | 2 0
cont. binary integer | cont. binary integer
4368 1680 0 | 3841 1670 0
Number of constraints: 9252 | 8684
eq. ineq. range | eq. ineq. range
linear: 851 6721 0 | 334 6680 0
quadratic: 1680 0 0 | 1670 0 0
Number of nonzeros:
objective Jacobian Hessian | objective Jacobian Hessian
linear: 1680 18993 | 1670 17867
quadratic: 0 3360 1680 | 0 3340 1670
total: 1680 22353 1680 | 1670 21207 1670
Knitro using Branch and Bound method with 1 thread.
Initial points
--------------
No initial point provided for the root node relaxation.
No primal point provided for the MIP.
Coefficient range:
linear objective: [3e+00, 1e+02] | [3e+00, 1e+02]
linear constraints: [1e-07, 4e+03] | [1e-06, 9e+01]
quadratic objective: [0e+00, 0e+00] | [0e+00, 0e+00]
quadratic constraints: [5e-03, 9e-03] | [5e-03, 9e-03]
variable bounds: [1e+00, 3e+07] | [1e+00, 3e+07]
constraint bounds: [3e+01, 2e+07] | [2e+01, 8e+03]
Root node relaxation
--------------------
Iter Objective Feasibility Optimality Time
error error (secs)
---- --------- ----------- ---------- ------
0 247761. 827.401 97.4654 0.249
1 249238. 827.397 12.5380 0.272
2 253050. 827.383 12.5380 0.286
3 263643. 827.338 12.5380 0.299
4 287276. 827.187 14.7175 0.312
5 336336. 826.672 18.1836 0.325
6 411702. 824.881 34.2399 0.338
7 524584. 816.759 60.2453 0.351
8 588576. 793.798 149.537 0.366
9 575258. 667.432 901.344 0.384
10 596826. 272.947 692.022 0.409
11 602935. 129.260 251.746 0.428
12 605967. 53.7016 414.947 0.448
13 595911. 23.4134 176.803 0.468
14 530045. 0.831240 482.595 0.489
15 589952. 1.96940e-02 208.573 0.511
16 662970. 2.10447e-02 52.8708 0.532
17 686642. 5.15266e-02 28.2246 0.550
18 700192. 0.178270 16.1639 0.566
19 721048. 0.124058 3698.36 0.593
20 731591. 5.62418e-02 853.742 0.619
30 739831. 8.60304e-08 1.18735e-04 0.841
Tree search
-----------
Nodes Best solution Best bound Gap Time
Expl | Unexpl value value (secs)
--------------- ------------- ---------- --- ------
1 2 739673. FCRD inf 2.544
EXIT: Terminating at first integer feasible point.
HINT: The problem may be a non-convex mixed-integer problem. Set
mip_multistart=1 to enable a mixed-integer multistart heuristic,
which may improve the chances of finding the global solution.
Final Statistics for MIP
------------------------
Final objective value = 7.39673482859438285e+05
Final bound value = +inf
Final optimality gap (abs / rel) = inf / inf
# of root cutting plane rounds = 1
# of restarts = 0
# of nodes processed = 1 (0.898s)
# of strong branching evaluations = 0 (0.000s)
# of function evaluations = 0 (0.000s)
# of gradient evaluations = 0 (0.000s)
# of hessian evaluations = 0 (0.000s)
# of hessian-vector evaluations = 0
# of subproblems processed = 2 (2.455s)
Total program time (secs) = 2.54436 (3.624 CPU time)
Time spent in evaluations (secs) = 0.00000
Cuts statistics (gen / add)
---------------------------
Knapsack cuts = 0 / 0
Mixed-integer rounding cuts = 0 / 0
Flow-cover cuts = 0 / 0
Probing cuts = 0 / 0
Heuristics statistics (calls / successes / time)
------------------------------------------------
Feasibility pump = 0 / 0 / 0.000s
Rounding heuristic = 1 / 1 / 1.566s
MPEC heuristic = 0 / 0 / 0.000s
Local search heuristic = 0 / 0 / 0.000s
===========================================================================
WARNING: Loading a SolverResults object with a warning status into
model.name="unknown";
- termination condition: maxIterations
- message from solver: Knitro 16.0.0\x3a MIP\x3a Integer feasible point
found.; objective 739673.4828594383; optimality gap Infinity; 1 nodes; 2
subproblem solvesShow the full outputHide the full output
Output visualization
What the schedule earns, and which turbines earn it:
report_schedule(plant, schedule)Revenue 739,673 over 168 hours
turbine hours on production revenue share
(MWh)
-----------------------------------------------
1 147 1,377 76,320 10.3%
2 157 2,187 115,938 15.7%
3 129 1,097 66,678 9.0%
4 0 0 0 0.0%
5 136 1,686 98,531 13.3%
6 144 1,686 94,640 12.8%
7 118 1,121 71,424 9.7%
8 155 1,971 104,212 14.1%
9 147 2,027 111,929 15.1%
10 0 0 0 0.0%
Reservoir 25.0M at the start, 25.0M at the end, 20.5M at its lowest
Net head 18.7m to 22.9m, so the same water is worth 1.22 times more at the top than at the bottomA turbine whose minimum production is more than its maximum discharge can reach at the heads the reservoir offers can never be switched on at all, so the model leaves it off for the whole week and it stands at zero hours in the table.
The schedule reads in four panels. They show the discharge of every turbine stacked into the total, the volume left in the reservoir, the production summed over the turbines, and the three water levels. The turbines run through the dear hours and idle through the cheap ones, refilling the reservoir while the water is worth less than it will be later. The levels follow. The headwater falls with the volume, the tailwater rises with what is discharged and spilled, and the net head, the gap between them, is what every unit of water is worth.
draw_schedule(plant, schedule)