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.
# install dependencies and select solver
%pip install -q amplpy pandas matplotlib numpy scipy
SOLVER = "highs"
from amplpy import AMPL, ampl_notebook
ampl = ampl_notebook(
modules=["highs"], # modules to install
license_uuid="default", # license to use
) # instantiate AMPL object and register magics
# Load packages
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
def get_firststagecost(model):
return model.get_value("sum{i in V: is_generator[i]==1} c[i]*p[i]")
def get_optimal_production_levels(model):
q = model.get_data("{i in CG} p[i]")
q = q.to_list()
q = [i[1] for i in q]
return [np.round(i, 2) for i in q]
def get_participation_factors(model):
q = [i[1] for i in model.get_data("{i in CG} alpha[i]")]
return [np.round(i, 2) for i in q]
def get_averagerecoursecost(model, df, T, abstotalimbalances):
q_alpha = [i[1] for i in model.var["alpha"].get_values()]
averagerecoursecost = (
1
/ T
* sum(
sum(
2 * df.iloc[i]["c"] * q_alpha[i] * abstotalimbalances[t]
for i in df.index.values
)
for t in range(T)
)
)
return averagerecoursecost
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 \(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, shall 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 $\( \begin{align*} \begin{array}{lllll} Q(\bar{p},\Delta) := &\min \quad & \sum_{i \in \mathcal{G}^{\text{coal}} \cup \mathcal{G}^{\text{gas}}} \hat{c}_i(r_i) \\ &\text{s.t.} \quad & p_i = \bar{p}_i + r_i & \forall \, i \in \mathcal{G}^{\text{coal}} \cup \mathcal{G}^{\text{gas}} \\ && p_{i}^{\min } \leq p_{i} \leq p_{i}^{\max } & \forall \, i \in \mathcal{G}^{\text{coal}} \cup \mathcal{G}^{\text{gas}} \\ && \sum_{j: (i, j) \in E} f_{ij} - \sum_{j: (j, i) \in E} f_{ji} = p_i - d_i & \forall \, i \in \mathcal{G}^{\text{coal}} \cup \mathcal{G}^{\text{gas}} \\ && \sum_{j: (i, j) \in E} f_{ij} - \sum_{j: (j, i) \in E} f_{ji} = \bar{p}_i + \Delta_i - d_i & \forall \, i \in \mathcal{G}^{\text{wind}} \cup \mathcal{G}^{\text{solar}}\\ && \sum_{j: (i, j) \in E} f_{ij} - \sum_{j: (j, i) \in E} f_{ji} = \bar{p}_i - d_i & \forall \, i \in V \setminus (\mathcal{G}^{\text{wind}} \cup \mathcal{G}^{\text{solar}} \cup \mathcal{G}^{\text{coal}} \cup \mathcal{G}^{\text{gas}}) \\ && f_{ij} = b_{ij}(\theta_i - \theta_j), & \forall \, (i, j) \in E \\ && -f_{ij}^{\max} \leq f_{ij} \leq f_{ij}^{\max} & \forall \, (i, j) \in E\\ && \theta_i \in \mathbb{R} & \forall \, i \in V \\ && f_{ij} \in \mathbb{R} & \forall \, (i, j) \in E\\ && r_i \in \mathbb{R} & \forall \, i \in \mathcal{G}^{\text{coal}} \cup \mathcal{G}^{\text{gas}}, \end{array} \end{align*} \)$
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, remembering that among \(\bar{p}_i\)’s, only those with \(i \in \mathcal{G}^{\text{coal}} \cup \mathcal{G}^{\text{gas}}\) are actual decision variables while \(\bar{p}_i\), \(i \in \mathcal{G}^{\text{solar}} \cup \mathcal{G}^{\text{wind}}\) are known parameters (think of them as wind and solar power forecasts).
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}} 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 \(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 $\( \begin{align*} \begin{array}{llll} \min \quad & \sum_{i \in \mathcal{G}^{\text{coal}} \cup \mathcal{G}^{\text{gas}}} c_i(\bar{p}_i) + \mathbb{E} Q(\bar{p}, \alpha, \Delta) \\ \text{s.t.} \quad & \sum_{i \in \mathcal{G}^{\text{coal}} \cup \mathcal{G}^{\text{gas}}} \alpha_i = 1 \\ & p_{i}^{\min } \leq \bar{p}_{i} \leq p_{i}^{\max } & \forall i \in \mathcal{G}^{\text{coal}} \cup \mathcal{G}^{\text{gas}} \\ & \alpha_i \geq 0 & \forall i \in \mathcal{G}^{\text{coal}} \cup \mathcal{G}^{\text{gas}} \end{array} \end{align*} \)\( where the second stage problem is \)\( \begin{align*} \begin{array}{lllll} Q(\bar{p},\alpha,\Delta) := &\min \quad & \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) \\ &\text{s.t.} \quad & p_{i}^{\min } \leq p_i - \alpha_i \sum_{j \in \mathcal{G}^{\text{wind}} \cup \mathcal{G}^{\text{solar}}} \Delta_j \leq p_{i}^{\max } & \forall i \in \mathcal{G}^{\text{coal}} \cup \mathcal{G}^{\text{gas}} \\ && \sum_{j: (i, j) \in E} f_{ij} - \sum_{j: (j, i) \in E} f_{ji} = \bar{p}_i - \alpha_i \sum_{j \in \mathcal{G}^{\text{wind}} \cup \mathcal{G}^{\text{solar}}} \Delta_j - d_i & \forall \, i \in \mathcal{G}^{\text{coal}} \cup \mathcal{G}^{\text{gas}} \\ && \sum_{j: (i, j) \in E} f_{ij} - \sum_{j: (j, i) \in E} f_{ji} = \bar{p}_i + \Delta_i - d_i & \forall \, i \in \mathcal{G}^{\text{wind}} \cup \mathcal{G}^{\text{solar}}\\ && \sum_{j: (i, j) \in E} f_{ij} - \sum_{j: (j, i) \in E} f_{ji} = \bar{p}_i - d_i & \forall \, i \in V \setminus (\mathcal{G}^{\text{wind}} \cup \mathcal{G}^{\text{solar}} \cup \mathcal{G}^{\text{coal}} \cup \mathcal{G}^{\text{gas}}) \\ && f_{ij} = b_{ij}(\theta_i - \theta_j), & \forall \, (i, j) \in E \\ && -f_{ij}^{\max} \leq f_{ij} \leq f_{ij}^{\max} & \forall (i, j) \in E\\ && \theta_i \in \mathbb{R} & \forall i \in V \\ && f_{ij} \in \mathbb{R} & \forall (i, j) \in E. \end{array} \end{align*} \)$
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: $\( \begin{align*} c_i(x) := c_i x, \quad \hat{c}_i(x) = 2 c_i |x|, \end{align*} \)\( 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.
AMPL 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.
# Download the data
base_url = (
"https://raw.githubusercontent.com/ampl/mo-book.ampl.com/dev/notebooks/10/"
)
nodes_df = pd.read_csv(base_url + "nodes.csv").set_index("node_id")
edges_df = pd.read_csv(base_url + "edges.csv").set_index(["node_id1", "node_id2"])
# Replace 'na' by an empty string
nodes_df.fillna("", inplace=True)
# Rename column acording to the model
nodes_df.rename(columns={"c_var": "c"}, inplace=True)
# Remove unused column
nodes_df.drop(["c_fixed"], axis=1, inplace=True)
network = {"nodes": nodes_df, "edges": edges_df}
SAA implementation#
The following cell is the AMPL 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.
%%writefile opf_ldr.mod
# Define sets
set T;
set V;
set E within V cross V;
set SWH;
set CG;
set NG;
# Define parameters
param uniformparticipationfactors;
param p_min{V} >= 0;
param p_max{i in V} >= p_min[i];
param is_generator{V};
param energy_type{V} symbolic;
param c{V};
param d{V};
param abs_total_imbalance{T} >= 0;
param total_imbalance{T};
param nT := card(T);
param b{E};
param f_max{E};
param imbalances{T, V};
# Declare decision variables
var p{i in V} >= p_min[i], <= p_max[i];
var r{V, T} >= 0;
var alpha{V} >= 0, <= 1;
var theta{V, T};
var f{E, T};
# Declare objective function including the recourse actions
minimize obj:
sum{i in V: is_generator[i] == 1} c[i] * p[i] +
1 / nT * sum{t in T, i in V: energy_type[i] == "coal" or energy_type[i] == "gas"}
2 * c[i] * alpha[i] * abs_total_imbalance[t];
# Declare constraints
s.t. windsolarhydro_nopartecipationfactors {i in SWH}: alpha[i] = 0;
s.t. load_nopartecipationfactors {i in NG}: alpha[i] = 0;
s.t. sum_one: sum{i in V} alpha[i] = 1;
s.t. uniformparticipationfactors_check {i in CG}:
uniformparticipationfactors == 1 ==> alpha[i] = 1 / card(CG);
# Second-stage production levels must also meet generator limits
s.t. power_withrecourse {i in V, t in T}:
r[i, t] = p[i] - alpha[i] * total_imbalance[t];
s.t. generation_upper_bound_withrecourse {i in CG, t in T}:
r[i, t] <= p_max[i];
s.t. generation_lower_bound_withrecourse {i in CG, t in T}:
r[i, t] >= p_min[i];
# Expressions for outgoing and incoming flows
s.t. flow_conservation {i in V, t in T}:
sum{(j, i) in E} f[j, i, t] -
sum{(i, j) in E} f[i, j, t]
== r[i, t] + imbalances[t, i] - d[i];
s.t. susceptance {(i, j) in E, t in T}:
f[i, j, t] = b[i, j] * (theta[i, t] - theta[j, t]);
s.t. flows_upper_bound {(i, j) in E, t in T}:
f[i, j, t] <= f_max[i, j];
s.t. flows_lower_bound {(i, j) in E, t in T}:
-f[i, j, t] <= f_max[i, j];
Overwriting opf_ldr.mod
# Define an OPF problem with recourse actions for the conventional generators based on participation factors
def OPF_participationfactors(
network,
imbalances,
totalimbalances,
abstotalimbalances,
uniformparticipationfactors=False,
):
nimbalances = len(imbalances)
temp = {}
for i, j in enumerate(imbalances):
for k, v in j.items():
temp[i, k] = v
imbalances = temp
nodes_df = network["nodes"]
edges_df = network["edges"]
ampl = AMPL()
ampl.read("opf_ldr.mod")
# load data in pandas data frames
ampl.set_data(nodes_df, "V")
ampl.set_data(edges_df, "E")
# data for remaining sets
ampl.set["T"] = list(range(nimbalances))
ampl.set["SWH"] = nodes_df[
nodes_df["energy_type"].isin(["wind", "solar", "hydro"])
].index
ampl.set["CG"] = nodes_df[nodes_df["energy_type"].isin(["coal", "gas"])].index
ampl.set["NG"] = nodes_df[nodes_df["energy_type"].eq("")].index
# data for parameters
ampl.param["abs_total_imbalance"] = abstotalimbalances
ampl.param["total_imbalance"] = totalimbalances
ampl.param["uniformparticipationfactors"] = 1 if uniformparticipationfactors else 0
ampl.param["imbalances"] = imbalances
ampl.option["solver"] = "highs"
ampl.option["presolve"] = 0 # or change presolve_eps
ampl.get_output("solve;")
return ampl
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.index.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,
)
firststagecost = get_firststagecost(m)
print(f"First-stage energy production cost = {firststagecost:.2f}")
print(
f"The optimal production levels for the conventional generators are: {get_optimal_production_levels(m)}"
)
First-stage energy production cost = 40385.16
The optimal production levels for the conventional generators are: [135.69, 227.38, 235.31, 371.35, 0, 0, 298.39, 336.95, 0, 0]
It is key to understand how this solution would perform if there were perturbations in renewable energy production.
First of all, it is not guaranteed that it is possible to find a feasible solution in such a case. Following an approach similar to that of the OPF problem with wind curtailment, one can solve for the remaining variables when keeping the initial solution \(\bar{p}\) fixed and check that with uniform participation factors, it is not possible to find a feasible flow in any of the scenarios we consider. If instead we allow for non-uniform participation factors, then this is not possible in 13 out of 100 scenarios.
Putting the feasibility issues aside for a moment, let us check how much extra cost would there be with uniform participation factors would there be on average across our scenarios. We can calculate this by taking the total imbalance and computing the cost of the recourse action to cover it, assuming that every coal and gas generator adjusts its production proportionally to the optimal participation factor previously obtained.
averagerecoursecost = get_averagerecoursecost(m, nodes_df, T, abstotalimbalances)
averagetotalcost = m.obj["obj"].value() + averagerecoursecost
print(
f"Average energy production cost due to recourse actions = {averagerecoursecost:.2f} (but including many infeasible scenarios!)"
)
print(
f"Average total cost = {averagetotalcost:.2f} (but including many infeasible scenarios!)"
)
Average energy production cost due to recourse actions = 5354.25 (but including many infeasible scenarios!)
Average total cost = 45739.41 (but including many infeasible scenarios!)
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 the ten conventional generators in \(\in \mathcal{G}^{\text{coal}} \cup \mathcal{G}^{\text{gas}}\).
m = OPF_participationfactors(
network,
imbalances,
totalimbalances,
abstotalimbalances,
uniformparticipationfactors=True,
)
print(
"The optimal production levels for the conventional generators are",
get_optimal_production_levels(m),
)
q = [i[1] for i in m.get_data("{i in CG} alpha[i]")]
print(
"The participation factors for the conventional generators are",
[np.round(i, 2) for i in q],
)
firststagecost = get_firststagecost(m)
averagerecoursecost = get_averagerecoursecost(m, nodes_df, T, abstotalimbalances)
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.obj['obj'].value():.2f}")
The optimal production levels for the conventional generators are [127.0, 205.65, 213.57, 349.62, 27.39, 27.39, 276.66, 323.04, 27.39, 27.39]
The participation factors for the conventional generators are [0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1]
First-stage energy production cost = 41719.19
Average energy production cost due to recourse actions = 5354.25
Total cost = 47073.44
We see that the total average production cost is slightly higher than in the first “perfect forecast” nominal scenario, but the benefit of this newly obtained solution is that we are sure that in all scenarios we have a feasible power flow dispatch.
We argue that using this solution should be preferable: even if the ‘nominal’ solution had a slightly lower average production cost, we need to factor in how costly will it be for the network operator when the solution becomes infeasible. Having an infeasible network configuration means there is a risk of cascading failures and/or blackout, which besides possibly damaging the infrastructure is dramatically more expensive from a financial and societal perspective, and having a 13% chance of this happening is just unaffordable.
Next, we also optimize the participation factors \(\alpha_i\)’s jointly with the initial power levels, to see if we can achieve a reduction in the average total cost. We can do this using the same function as before but changing the argument uniformparticipationfactors to False.
m = OPF_participationfactors(
network,
imbalances,
totalimbalances,
abstotalimbalances,
uniformparticipationfactors=False,
)
print(
"The optimal production levels for the conventional generators are",
get_optimal_production_levels(m),
)
print(
"The optimal participation factors for the conventional generators are",
get_participation_factors(m),
)
firststagecost = get_firststagecost(m)
averagerecoursecost = get_averagerecoursecost(m, nodes_df, T, abstotalimbalances)
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.obj['obj'].value():.2f}")
The optimal production levels for the conventional generators are [134.75, 227.38, 235.31, 371.35, 0, 0, 298.39, 281.66, 0, 56.24]
The optimal participation factors for the conventional generators are [0.06, 0, -0.0, -0.0, 0, 0, 0, 0.73, 0, 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.