10.3 Optimal power flow problem with recourse actions#
In this notebook we illustrate an application of the idea of linear decision rules to a two-stage optimal power flow problem in which the power of the conventional generators has to adapt automatically to balance out energy surplus/shortage due to changes in renewable resources input.
We will be working with a modified version of the energy dispatch problem and the OPF problem with wind curtailment. Please refer to the corresponding notebooks for the basics of power networks.
Preamble: Install Pyomo and a solver#
The following cell sets and verifies a global SOLVER for the notebook. If run on Google Colab, the cell installs Pyomo and the HiGHS solver, while, if run elsewhere, it assumes Pyomo and HiGHS have been previously installed. It then sets to use HiGHS as solver via the appsi module and a test is performed to verify that it is available. The solver interface is stored in a global object SOLVER for later use.
import sys
if 'google.colab' in sys.modules:
%pip install pyomo >/dev/null 2>/dev/null
%pip install highspy >/dev/null 2>/dev/null
solver = 'appsi_highs'
import pyomo.environ as pyo
SOLVER = pyo.SolverFactory(solver)
assert SOLVER.available(), f"Solver {solver} is not available."
Problem description#
We consider a variant of the Optimal Power Flow problem in which each conventional generator \(i\) commits in advance to produce a specific amount \(\bar{p}_i\) of energy as determined by the OPF problem assuming the renewable energy production from all solar panels and wind turbines will be equal to the forecasted one, also denoted as \(\bar{p}_j\). The realized renewable energy output of generator \(j \in \mathcal{G}^{\text{wind}} \cup \mathcal{G}^{\text{solar}}\), however, may deviate from their forecast values, by an amount \(\Delta_j\), and result in a power production of
Then, the conventional generators need to take a recourse action to make sure that the network is balanced, i.e., that the total energy production equals the total energy demand. This means that the problem has a two-stage structure:
first, the ‘nominal’ energy generation levels are set for the coal and gas units
then, the actual renewable energy output of the wind/solar generators is observed
power generation levels of the coal and gas units need to be adapted.
If we were optimizing for the average-case total cost, then a proper two-stage formulation of our problem would be as follows:
where the second stage problem is
where \(c_i(.)\) and \(\hat{c}_i(.)\) are the cost functions related to the pre-committed \(\bar{p}_i\) and changed \(r_i\) amounts of energy, with first-stage decisions \(\bar{p}_i\) for coal, gas and hydro generators. Wind and solar forecasts are fixed parameters. Hydro production is chosen within its bounds and remains unchanged across scenarios; in these data it has zero variable cost. At load nodes \(\bar{p}_i=0\). Thus the final balance row in each recourse formulation covers both hydro and load nodes.
In the OPF problem with wind curtailment we solved a similar problem using the SAA method. More specifically, we generated a number of scenarios for the unknown parameters and combined the second-stage problems related to each of these scenarios with the single first-stage problem obtaining one large MILO problem. In this notebook, we are going to use a slightly different approach, namely implementing a specific type of linear decision rules. For each conventional generator, we set a recourse action, that is real-time adjustment of its power production, based on the realization of the renewable energy. More specifically, each conventional generator \(i \in \mathcal{G}^{\text{coal}} \cup \mathcal{G}^{\text{gas}}\) has a participation factor \(\alpha_i \geq 0\) which determines to which extent that generator responds to the total imbalance \(\sum_j \Delta_j\). Specifically, the power production after the recourse action at the conventional generator \(i\) is denoted by \(p_i\) and is given by
In the notation of the problem above, we picked the recourse action \(r_i\) to be a function of the uncertain parameters \(\Delta\)’s rather than a decision variable, since
The participation factor \(\alpha_i \in [0,1]\) indicates the fraction of the power imbalance that generator \(i\) needs to help compensate. To ensure that the power balance is satisfied, we need to have \(\sum_{i \in \mathcal{G}^{\text{coal}} \cup \mathcal{G}^{\text{gas}}} \alpha_i = 1\). Indeed, in this case, assuming the power was balanced in the first stage, i.e., \(\sum_{i \in \mathcal{G}} \bar{p}_i - \sum_{i \in V} d_i =0\), then the net power balance after the second stage is
The participation factors \(\alpha_i\)’s do not have to be equal for all generators and in fact, they can be optimized jointly together with the initial power levels \(\bar{p}_i\). Since the energy produced as recourse action is more expensive, we account for this by adding to the objective function the cost term \(\sum_{i \in \mathcal{G}^{\text{coal}} \cup \mathcal{G}^{\text{gas}}} \hat{c}_i(-\alpha_i \sum_{j \in \mathcal{G}^{\text{wind}} \cup \mathcal{G}^{\text{solar}}} \Delta_j)\) for some cost functions \(\hat{c}_i(.)\).
The resulting two-stage stochastic problem is
where the second stage problem is
The only remaining question is where to take the values for \(\Delta_j\)’s from. One way of doing things would be to construct an uncertainty set for \(\Delta_j\) and make sure that the inequality constraints hold for all realizations of \(\Delta_j\) using the technique of robust counterparts/adversarial approach from Chapter 8. This would be a solid approach if we were optimizing for the worst-case value of the objective function, instead of an expectation.
However, since in this particular application we optimize for the expected value it makes sense for us to resort to the SAA method. More specifically, we sample \(T\) realizations of the renewable fluctuations \(\Delta\)’s, denoted as \(\{\Delta^s\}_{s=1,\dots,T}\), and we approximate the expectation through an empirical average across all those samples while enforcing that the constraints hold for every such realization. In this way, the resulting problem we actually solve is:
To recover a linear problem, we assume that the energy generation costs are modeled as:
where \(c_i\) is the unit production cost of the \(i\)-th generator. This structure for the \(\hat{c}_i(\cdot)\) functions means we assume that any real-time adjustment in the energy dispatch of a generator is twice as costly as a pre-scheduled unit of energy generated there.
The adjustment in the recourse cost has a minus sign because it compensates the renewable imbalance. For the absolute-value cost used here, the sign does not change the numerical value; a general asymmetric adjustment cost would depend on it.
Pyomo solution#
Data imports#
Importantly for our problem, the variable costs of producing a unit of energy per generator are stored in their corresponding ‘c_var’ attributes. For the rest, we are using the same data elements as in the OPF problem with wind curtailment.
The CSVs bundled in notebooks/10/data are used when running from the notebook folder or the repository root. For a standalone notebook, such as one opened in Colab, the same data are read from a fixed repository snapshot. Keep nodes.csv and edges.csv together when supplying local files.
# Load packages
import pyomo.environ as pyo
from IPython.display import Markdown, HTML
import numpy as np
import numpy.random as rnd
import pandas as pd
import matplotlib.pyplot as plt
from ast import literal_eval as make_tuple
import networkx as nx
import time
# Use the bundled data when available; otherwise download the same snapshot.
from pathlib import Path
data_files = ("nodes.csv", "edges.csv")
data_dir = next(
(folder for folder in (Path("data"), Path("notebooks/10/data"))
if any((folder / name).is_file() for name in data_files)),
None,
)
if data_dir is not None:
if not all((data_dir / name).is_file() for name in data_files):
raise FileNotFoundError(f"Keep nodes.csv and edges.csv together in {data_dir}.")
data_source = str(data_dir)
else:
data_source = (
"https://raw.githubusercontent.com/mobook/MO-book/"
"e745704ca19c0eef476beb32bb24dfba80244043/notebooks/10/data"
)
nodes_df = pd.read_csv(f"{data_source}/nodes.csv", index_col=0)
edges_df = pd.read_csv(f"{data_source}/edges.csv", index_col=0)
# Rework the nodes and edges data into a form that is convenient for NetworkX to work with
nodes = nodes_df.set_index("node_id").T.to_dict()
edges = edges_df.set_index(edges_df["edge_id"].apply(make_tuple)).T.to_dict()
I = {"nodes": nodes, "edges": edges}
# Initialize a network for demonstration purposes
network = I
SAA implementation#
The following cell is the Pyomo function implementation of our final optimization problem above, with an optional argument indicating if the participation factors should be taken as uniform across all conventional generators. For ease of accounting for the second-stage costs \(\hat{c}_i(\cdot)\), we assume that the function takes as an argument not only a list of scenarios for the renewable energy generation, but also the corresponding total and absolute values of the imbalances in the network.
In the code, model.p denotes the scheduled production \(\bar p\). model.r[i,t] stores \(\bar p_i-\alpha_i\sum_j\Delta_j^t\), which is the realized production for a conventional generator, rather than the adjustment \(r_i\) in the mathematical notation. Wind and solar deviations are added separately in the node balances. The displayed balances use outgoing minus incoming flow, whereas the code uses incoming minus outgoing flow. Negating all flows and phase angles maps one convention to the other and leaves the production decisions and costs unchanged.
# OPF with recourse actions for conventional generators based on participation factors
def OPF_participationfactors(
network,
imbalances,
totalimbalances,
abstotalimbalances,
uniformparticipationfactors=False,
):
# Define a model
model = pyo.ConcreteModel("OPF with participation factors")
# Define sets
model.T = pyo.Set(initialize=range(len(imbalances)))
model.V = pyo.Set(initialize=network["nodes"].keys())
model.E = pyo.Set(initialize=network["edges"].keys())
model.SWH = pyo.Set(
initialize=[
i
for i, data in network["nodes"].items()
if data["energy_type"] in ["wind", "solar", "hydro"]
]
)
model.CG = pyo.Set(
initialize=[
i
for i, data in network["nodes"].items()
if data["energy_type"] in ["coal", "gas"]
]
)
model.NG = pyo.Set(
initialize=[
i for i, data in network["nodes"].items() if pd.isna(data["energy_type"])
]
)
# Declare decision variables
model.p = pyo.Var(model.V, domain=pyo.NonNegativeReals)
model.r = pyo.Var(model.V, model.T, domain=pyo.NonNegativeReals)
model.alpha = pyo.Var(model.V, domain=pyo.NonNegativeReals)
model.theta = pyo.Var(model.V, model.T, domain=pyo.Reals)
model.f = pyo.Var(model.E, model.T, domain=pyo.Reals)
model.abs_total_imbalance = pyo.Param(
model.T, domain=pyo.NonNegativeReals, initialize=abstotalimbalances
)
model.total_imbalance = pyo.Param(
model.T, domain=pyo.Reals, initialize=totalimbalances
)
# Declare objective function including the recourse actions
model.objective = pyo.Objective(
expr=sum(
data["c_var"] * model.p[i]
for i, data in network["nodes"].items()
if data["is_generator"]
)
+ 1
/ len(model.T)
* sum(
sum(
2 * data["c_var"] * model.alpha[i] * model.abs_total_imbalance[t]
for i, data in network["nodes"].items()
if data["energy_type"] in ["coal", "gas"]
)
for t in model.T
),
sense=pyo.minimize,
)
# Declare constraints
# First-stage production levels must meet generator limits
model.generation_upper_bound = pyo.Constraint(
model.V, rule=lambda m, i: m.p[i] <= network["nodes"][i]["p_max"]
)
model.generation_lower_bound = pyo.Constraint(
model.V, rule=lambda m, i: network["nodes"][i]["p_min"] <= m.p[i]
)
# Wind, solar and hydro generators have zero participation factors
model.windsolarhydro_nopartecipationfactors = pyo.Constraint(
model.SWH, rule=lambda m, i: m.alpha[i] == 0
)
# Load nodes have zero participation factors
model.load_nopartecipationfactors = pyo.Constraint(
model.NG, rule=lambda m, i: m.alpha[i] == 0
)
# Participation factors must sum to one
model.sum_one = pyo.Constraint(expr=sum(model.alpha[i] for i in model.V) == 1)
if uniformparticipationfactors:
# Participation factors must be equal
model.equal_participationfactors = pyo.Constraint(
model.CG, rule=lambda m, i: m.alpha[i] == 1 / len(model.CG)
)
# Second-stage production levels must also meet generator limits
model.power_withrecourse = pyo.Constraint(
model.V,
model.T,
rule=lambda m, i, t: m.r[i, t] == m.p[i] - m.alpha[i] * m.total_imbalance[t],
)
model.generation_upper_bound_withrecourse = pyo.Constraint(
model.CG,
model.T,
rule=lambda m, i, t: m.r[i, t] <= network["nodes"][i]["p_max"],
)
model.generation_lower_bound_withrecourse = pyo.Constraint(
model.CG,
model.T,
rule=lambda m, i, t: network["nodes"][i]["p_min"] <= m.r[i, t],
)
# Expressions for outgoing and incoming flows
model.outgoing_flow = pyo.Expression(
model.V,
model.T,
rule=lambda m, i, t: sum(m.f[i, j, t] for j in model.V if (i, j) in model.E),
)
model.incoming_flow = pyo.Expression(
model.V,
model.T,
rule=lambda m, i, t: sum(m.f[j, i, t] for j in model.V if (j, i) in model.E),
)
# Net power production at each node after recourse actions
model.flow_conservation = pyo.Constraint(
model.V,
model.T,
rule=lambda m, i, t: m.incoming_flow[i, t] - m.outgoing_flow[i, t]
== m.r[i, t] + imbalances[t][i] - network["nodes"][i]["d"],
)
model.susceptance = pyo.Constraint(
model.E,
model.T,
rule=lambda m, i, j, t: m.f[(i, j), t]
== network["edges"][(i, j)]["b"] * (m.theta[i, t] - m.theta[j, t]),
)
model.flows_upper_bound = pyo.Constraint(
model.E,
model.T,
rule=lambda m, i, j, t: m.f[(i, j), t] <= network["edges"][(i, j)]["f_max"],
)
model.flows_lower_bound = pyo.Constraint(
model.E,
model.T,
rule=lambda m, i, j, t: -m.f[(i, j), t] <= network["edges"][(i, j)]["f_max"],
)
return model
Scenario generation#
Next, we generate \(T=100\) scenarios in which the wind and solar production deviate from the forecasted values. Such deviations, named imbalances, are generated uniformly at random assuming the realized wind or solar power is between 0.5 and 1.5 times the forecasted value.
For ease of calculations, for each scenario, we define a separate data structure with the total energy imbalance and the total absolute imbalance.
# Define the set of nodes with possible deviations from forecasts, i.e. those with either a wind or a solar generator
SW = {48, 53, 58, 60, 64, 65}
SW_df = nodes_df[nodes_df["node_id"].isin(SW)]
# Define the number of scenarios and the random seed
T = 100
seed = 0
rng = np.random.default_rng(seed)
# Imbalances are generated uniformly at random assuming the realized
# wind or solar power is between 0.5 and 1.5 times the forecasted value
imbalances = [
{
i: (
rng.uniform(-nodes_df["p_min"][i] / 2, nodes_df["p_min"][i] / 2)
if i in SW
else 0
)
for i in nodes_df.index
}
for t in range(T)
]
totalimbalances = {t: sum(imbalances[t].values()) for t in range(len(imbalances))}
abstotalimbalances = {t: abs(totalimbalances[t]) for t in range(len(totalimbalances))}
Perfect forecast case (no imbalances)#
We first solve the optimization model in the case where the forecast for solar and wind power are perfect, meaning there is no imbalance. In this case, the recourse actions are not needed and the second stage part of the problem is trivial.
# Define trivial arrays for the case of perfect forecast and no need of recourse actions
zeroimbalances = [{i: 0 for i in nodes_df.index}]
zerototalimbalances = {0: sum(zeroimbalances[0].values())}
zeroabstotalimbalances = {0: abs(zerototalimbalances[0])}
# Solve the model
m = OPF_participationfactors(
network,
zeroimbalances,
zerototalimbalances,
zeroabstotalimbalances,
uniformparticipationfactors=True,
)
results = SOLVER.solve(m)
pyo.assert_optimal_termination(results)
firststagecost = sum(
data["c_var"] * m.p[i].value
for i, data in network["nodes"].items()
if data["is_generator"]
)
print(f"First-stage energy production cost = {firststagecost:.2f}")
print(
f"The optimal production levels for the conventional generators are: {[np.round(m.p[i].value, 2) for i in m.CG]}"
)
First-stage energy production cost = 40385.16
The optimal production levels for the conventional generators are: [np.float64(135.69), np.float64(227.38), np.float64(235.31), np.float64(371.35), np.float64(-0.0), np.float64(-0.0), np.float64(298.39), np.float64(336.95), np.float64(-0.0), np.float64(-0.0)]
It is key to understand how the nominal solution would perform if renewable production differed from the forecast. Balancing total supply and demand does not by itself ensure a feasible dispatch: generator and line limits must also hold.
The next calculation evaluates the adjustment cost for the nominal production plan with uniform participation factors. It is a cost-only estimate; it does not solve the scenario flow problems or establish feasibility. A feasibility count would also need to specify whether participation factors are fixed in advance or reoptimized separately after each scenario is observed.
averagerecoursecost = (
1
/ T
* sum(
sum(
2 * data["c_var"] * m.alpha[i].value * abstotalimbalances[t]
for i, data in network["nodes"].items()
)
for t in range(T)
)
)
averagetotalcost = m.objective() + 1 / T * sum(
sum(
2 * data["c_var"] * m.alpha[i].value * abstotalimbalances[t]
for i, data in network["nodes"].items()
)
for t in range(T)
)
print(
f"Average energy production cost due to recourse actions = {averagerecoursecost:.2f} (cost-only estimate; network feasibility not checked)"
)
print(
f"Average total cost = {averagetotalcost:.2f} (cost-only estimate; network feasibility not checked)"
)
Average energy production cost due to recourse actions = 5354.25 (cost-only estimate; network feasibility not checked)
Average total cost = 45739.41 (cost-only estimate; network feasibility not checked)
Stochastic case (nonzero imbalances)#
If we assume that the forecast for solar and wind power is not perfect, then the total energy imbalance in the network will be nonzero in each scenario. The resulting average total cost of energy production would be much higher than the deterministic scenario. This is intuitive because recourse actions are needed to cover the imbalance, and the recourse actions are much more expensive than the first-stage production-level decisions.
We now solve the two-stage stochastic optimization model that accounts for the fluctuations of solar and wind power from their forecasts, using the 100 generated scenarios. In this case, the recourse actions are still needed, but we assume fixed uniform participation factors equal to \(0.1\) for all ten conventional generators in \(\mathcal{G}^{\mathrm{coal}}\cup\mathcal{G}^{\mathrm{gas}}\).
m = OPF_participationfactors(
network,
imbalances,
totalimbalances,
abstotalimbalances,
uniformparticipationfactors=True,
)
results = SOLVER.solve(m)
pyo.assert_optimal_termination(results)
print(
"The optimal production levels for the conventional generators are",
[np.round(m.p[i].value, 2) for i in m.CG],
)
print(
"The participation factors for the conventional generators are",
[np.round(m.alpha[i].value, 2) for i in m.CG],
)
firststagecost = sum(
data["c_var"] * m.p[i].value
for i, data in network["nodes"].items()
if data["is_generator"]
)
averagerecoursecost = (
1
/ T
* sum(
sum(
2 * data["c_var"] * m.alpha[i].value * m.abs_total_imbalance[t]
for i, data in network["nodes"].items()
)
for t in m.T
)
)
print(f"\nFirst-stage energy production cost = {firststagecost:.2f}")
print(
f"Average energy production cost due to recourse actions = {averagerecoursecost:.2f}"
)
print(f"Total cost = {m.objective():.2f}")
The optimal production levels for the conventional generators are [np.float64(127.0), np.float64(205.65), np.float64(213.57), np.float64(349.62), np.float64(27.39), np.float64(27.39), np.float64(276.66), np.float64(323.04), np.float64(27.39), np.float64(27.39)]
The participation factors for the conventional generators are [np.float64(0.1), np.float64(0.1), np.float64(0.1), np.float64(0.1), np.float64(0.1), np.float64(0.1), np.float64(0.1), np.float64(0.1), np.float64(0.1), np.float64(0.1)]
First-stage energy production cost = 41719.19
Average energy production cost due to recourse actions = 5354.25
Total cost = 47073.44
The stochastic model chooses scheduled production so that a feasible dispatch exists for every sampled scenario under the uniform participation rule. Its objective includes the average cost of the required adjustments.
This feasibility statement concerns the sampled network model. It neither guarantees feasibility for every future realization nor estimates a probability of cascading failures or blackouts.
Next, we optimize the participation factors jointly with the scheduled production, keeping one common set of factors across all scenarios. We use the same function with uniformparticipationfactors=False.
m = OPF_participationfactors(
network,
imbalances,
totalimbalances,
abstotalimbalances,
uniformparticipationfactors=False,
)
results = SOLVER.solve(m)
pyo.assert_optimal_termination(results)
print(
"The optimal production levels for the conventional generators are",
[np.round(m.p[i].value, 2) for i in m.CG],
)
print(
"The optimal participation factors for the conventional generators are",
[np.round(m.alpha[i].value, 2) for i in m.CG],
)
firststagecost = sum(
data["c_var"] * m.p[i].value
for i, data in network["nodes"].items()
if data["is_generator"]
)
averagerecoursecost = (
1
/ T
* sum(
sum(
2 * data["c_var"] * m.alpha[i].value * m.abs_total_imbalance[t]
for i, data in network["nodes"].items()
)
for t in m.T
)
)
print(f"\nFirst-stage energy production cost = {firststagecost:.2f}")
print(
f"Average energy production cost due to recourse actions = {averagerecoursecost:.2f}"
)
print(f"Total cost = {m.objective():.2f}")
The optimal production levels for the conventional generators are [np.float64(134.75), np.float64(227.38), np.float64(235.31), np.float64(371.35), np.float64(0.0), np.float64(-0.0), np.float64(298.39), np.float64(281.66), np.float64(-0.0), np.float64(56.24)]
The optimal participation factors for the conventional generators are [np.float64(0.06), np.float64(-0.0), np.float64(-0.0), np.float64(0.0), np.float64(0.0), np.float64(0.0), np.float64(-0.0), np.float64(0.73), np.float64(0.0), np.float64(0.21)]
First-stage energy production cost = 40446.13
Average energy production cost due to recourse actions = 5969.05
Total cost = 46415.18
This energy dispatch is about 1.4% cheaper than the solution with uniform participation factors. It might seem like a small difference but in view of the high volumes of energy produced and consumed, say at the national level, this makes a huge difference.