import cmath
import math
from dataclasses import dataclass, field
from plot import draw_dispatch, draw_layout
from report import report_dispatch, report_network
@dataclass
class Bus:
number: int
min_voltage: float
max_voltage: float
active_demand: float
reactive_demand: float
@dataclass
class Branch:
from_bus: int
to_bus: int
y_kk: complex
y_km: complex
y_mm: complex
max_flow: float
ratio: float
phase_shift: float
min_angle: float
max_angle: float
@dataclass
class Generator:
bus: int
min_active_power: float
max_active_power: float
min_reactive_power: float
max_reactive_power: float
cost: list[float]
@dataclass
class Network:
buses: list[Bus]
branches: list[Branch]
generators: list[Generator]
base_power: float
positions: list[tuple[float, float]]
name: str
index: dict[int, int] = field(init=False)
def __post_init__(self):
self.index = {bus.number: k for k, bus in enumerate(self.buses)}
@property
def num_buses(self):
return len(self.buses)
@property
def num_branches(self):
return len(self.branches)
@property
def num_generators(self):
return len(self.generators)
@dataclass
class Dispatch:
cost: float
v: list[float]
theta: list[float]
p_km: list[float]
q_km: list[float]
p_mk: list[float]
q_mk: list[float]
p_g: list[float]
q_g: list[float]AC Optimal Power Flow
Dispatch generators over a transmission network at minimum cost while respecting the nonlinear AC power-flow physics, voltage limits, and line ratings.
Introduction
An electricity transport network allows the electricity produced by producers (nuclear power plants, solar panels, wind turbines, etc.) to be delivered to consumers (electricity distributors or industrial consumers). At every moment, production must be exactly equal to consumption.
The optimal power flow (OPF) is the problem of optimizing power flows on this type of network. Most electrical grids operate on alternating current, which involves the complex impedance of the equipment; this is referred to as ACOPF (alternating current OPF).
The objective of an ACOPF problem is to minimize operational costs such that the underlying grid constraints on generation, demand, and voltage and power flow limits are satisfied.
A simple power flow model is presented here; real problems are more complex and solving them requires specialized tools like OpenReac, which itself relies on Knitro.
Every figure on this page is read the same way. A branch carries a pale band as wide as the limit on the power that may cross it, and a darker core inside it as wide as the power actually crossing, so the room a branch has left is the gap between the two. The red arrowhead points the way the active power goes, and it is the only mark that means a direction. A generator’s disc grows with what it produces, and the short dark stub on a bus is the load drawn there, with its size in MVA beside it.
Physical explanation
An electrical branch is modeled by a “\(\Pi\)-type two-port network” preceded by a transformer (which acts on the magnitude and/or phase of the voltage).
The physical parameters associated with an electrical branch \(km\), where \(k\) is the from bus and \(m\) is the to bus, are:
- a resistance \(r_{km}\)
- a reactance \(x_{km}\)
- a shunt susceptance \(b_{km}\)
- a transformer ratio \(\tau _{km}\)
- a transformer phase-shift angle \(\phi_{km}\)
The resistance \(r_{km}\) and the reactance \(x_{km}\) define what is known as the complex impedance \(z_{km}\) of the \(\Pi\) model: \(\qquad z_{km} = r_{km} + jx_{km}\)
The series admittance \(y_{km}\) is the inverse of the complex impedance: \(\qquad y_{km} = \dfrac{1}{z_{km}}\)
The admittance matrix \(Y_{km}\) takes into account both the physical parameters of the \(\Pi\)-modeling \((y_{km},~ b_{km})\) and those of the transformer \((\tau _{km},~ \phi_{km})\). It is defined as follows:
\[ Y_{km} = \begin{pmatrix} \frac{y_{km}+jb_{km}}{\tau _{km}^2} & −y_{km} \frac{1}{\tau _{km} e^{−j\phi_{km}}}\\ −y_{km} \frac{1}{\tau _{km} e^{j\phi_{km}}} & y_{km}+jb_{km} \end{pmatrix} \]
For each branch \(km\), the admittance matrix \(Y_{km}\) allows us to define the current as a function of the voltage in the following way:
\[ \begin{pmatrix} i_k\\ i_m \end{pmatrix} =Y_{km} \begin{pmatrix} v_k\\ v_m \end{pmatrix} \]
The complex power at the origin and destination of a branch is defined based on the voltage and the current. More specifically, on a branch \(km\) going from node \(k\) to node \(m\), we have:
\[ \begin{aligned} s_{km} &= v_k \cdot \overline{i_k}\\ s_{mk} &= v_m \cdot \overline{i_m} \end{aligned} \]
It is possible to express the power on the branches solely as a function of the voltage:
\[ \begin{aligned} s_{km} &= \dfrac{\overline{y_{km}} - j b_{km}}{\tau_{km}^2} |v_k|^2 - \overline{y_{km}} \dfrac{1}{\tau_{km}} e^{-j\phi_{km}} v_k \overline{v_m}\\ s_{mk} &= -\overline{y_{km}} \dfrac{1}{\tau_{km}} e^{j\phi_{km}} \overline{v_k} v_m + (\overline{y_{km}} - j b_{km}) |v_m|^2 \end{aligned} \]
Thus, the voltages allow us to express both the power and the currents on the branches. In practice, it is therefore sufficient to determine the voltages and generation powers in order to know all the physical quantities on the network.
Problem description
In what follows, in order to express the model in the nonlinear programming formalism, we will represent the power quantities in Cartesian coordinates: a complex power \(S\) is written as \(S = P + jQ\), where \(P\) is the real part of the power, called the active power, and \(Q\) is the imaginary part, called the reactive power.
Input:
- The power system is represented by a graph \((N, B)\), where nodes \(N\) stand for buses (production and/or consumption points), and arcs \(B\) stand for branches. We denote as \(G\) the set of generators (\(G\subset N\)).
- For each branch \(km\in B\):
- \(\tau_{km}>0\) and \(\phi_{km}\) the ratio and the angle in the case of transformers (for a non-transformer branch, \(\tau_{km}=1\) and \(\phi_{km}=0\))
- the admittance coefficients of branch \(km\), derived from its physical parameters. The transformer phase shift \(\phi_{km}\) is not folded into these coefficients; it is handled through \(\theta_{km}\) in the model. Note that the off-diagonal coefficient is symmetric: \(G_{km}+jB_{km} = G_{mk}+jB_{mk}\). \[G_{kk}+jB_{kk}=\dfrac{y_{km}+jb_{km}}{\tau _{km}^2}, \qquad G_{km}+jB_{km}=-\dfrac{y_{km}}{\tau _{km}}, \qquad G_{mm}+jB_{mm}=y_{km}+jb_{km}\]
- \(L_{km} \in \mathbb{R}\) the limit on the power that flows through the branch
- For each bus \(k\in N\):
- \(V_{k}^\text{min}\) and \(V_{k}^\text{max}\) the voltage limits of bus \(k\), with \(0\le V_{k}^\text{min}\le V_{k}^\text{max}\)
- \(P_{k}^d\) and \(Q_{k}^d\) the active and reactive loads (demands), \(S_{k}^d =P_{k}^d +jQ_{k}^d\)
- \(G_k\) the set of generators located at bus \(k\)
- For each generator \(i\in G\):
- \(P_{i}^\text{min}\), \(P_{i}^\text{max}\), \(Q_{i}^\text{min}\) and \(Q_{i}^\text{max}\) the active and reactive generation limits, with \(P_{i}^\text{min}\le P_{i}^\text{max}\) and \(Q_{i}^\text{min}\le Q_{i}^\text{max}\)
- \(F_{i}\) a polynomial function defining the generation cost, represented by a cost vector \([C_{i,n}, ..., C_{i,1}, C_{i,0}]\) such that \(F_i=C_{i,n}X^n + ... + C_{i,1}X + C_{i,0}\).
- We don’t consider the piecewise linear case here.
We assume that:
- there are no multiple parallel branches \(km\)
- there may be multiple generators at a given bus
- the same bus may have generators and a nonzero load (demand)
Problem: find the voltage of each bus and the generation of each generator that minimize operational costs, such that the underlying grid constraints on generation, demand, and voltage and power flow limits are satisfied.
Nonlinear model
The power flow is expressed using the un-rotated admittance coefficients defined above; the transformer phase shift \(\phi_{km}\) enters only through the shifted angle difference \(\theta_{km}\).
Variables
- \(v_k = |v_k|e^{j\theta _k}\in \mathbb{C}\), \(V_k^\text{min}\le |v_k|\le V_k^\text{max}\), \(k\in N\), the complex voltage at bus \(k\)
- \(s_{km} = p_{km} + jq_{km}, s_{mk} = p_{mk} + jq_{mk} \in \mathbb{C}\), \(km \in B\), the complex power injected into the branch at \(k\) and at \(m\) (due to losses from the Joule effect, the power varies between the source node and the destination node)
- \(s_i^g=p_i^g+jq_i^g \in \mathbb{C}\), \(P_i^\text{min}\le p_i^g\le P_i^\text{max}\) and \(Q_i^\text{min}\le q_i^g\le Q_i^\text{max}\), \(i\in G\), the generation of generator \(i\)
- Three auxiliary variables per branch, introduced so that the nonconvex voltage products appear only in their definitions and the power flow equations below remain linear in \(G\), \(B\), \(\cos_{km}\), and \(\sin_{km}\):
- \(\theta _{km} = \theta _k - \theta _m - \phi _{km}\), the shifted angle difference (this is where the transformer phase shift enters the model)
- \(\cos _{km} = |v_k| \, |v_m| \cos (\theta _{km})\)
- \(\sin _{km} = |v_k| \, |v_m| \sin (\theta _{km})\)
Objective: minimize operational costs
Each generator \(i\) has a polynomial cost \(F_i\) of its active power output. The objective is the total generation cost:
\[ \min \quad \sum_{i\in G} F_i(p_i^g) \quad = \quad \min \quad \sum_{i\in G} C_{i,n}(p_i^g)^n + ... + C_{i,1}p_i^g + C_{i,0} \]
Constraints
- AC power flow laws: splitting the admittance coefficients into their real and imaginary parts, the power at each end of a branch becomes:
\[ \begin{aligned} p_{km} &= G_{kk}|v_k|^2 + G_{km}\cos_{km} + B_{km}\sin_{km}\\ q_{km} &= -B_{kk}|v_k|^2 - B_{km}\cos_{km} + G_{km}\sin_{km}\\ p_{mk} &= G_{mm}|v_m|^2 + G_{km}\cos_{km} - B_{km}\sin_{km}\\ q_{mk} &= -B_{mm}|v_m|^2 - B_{km}\cos_{km} - G_{km}\sin_{km} \end{aligned} \]
These four equations are linear in \(\cos_{km}\) and \(\sin_{km}\). The nonconvexity of the entire model is confined to the definitions of \(\cos_{km}\) and \(\sin_{km}\) above.
- Flow balance (excluding bus shunt admittance for simplicity; the branch line-charging susceptance \(b_{km}\) is still accounted for in \(B_{kk}\) and \(B_{mm}\)): \[ \forall k \in N, \qquad \sum_{i\in G(k)} s_i^g-S_k^d = \sum_{km\in \delta^+(k)} s_{km} + \sum_{mk\in \delta^-(k)} s_{mk} \]
where \(\delta^+(k)\) is the set of branches leaving bus \(k\) and \(\delta^-(k)\) the set arriving at it. The left-hand side is the net excess of generation over load at bus \(k\). The right-hand side is the total power injected into the grid at bus \(k\).
This gives two equations, for the imaginary and the real part:
\[ \forall k \in N, \qquad \begin{aligned} \sum_{i\in G(k)} p_i^g - P_k^d &= \sum_{km\in \delta^+(k)} p_{km} + \sum_{mk\in \delta^-(k)} p_{mk}\\ \sum_{i\in G(k)} q_i^g - Q_k^d &= \sum_{km\in \delta^+(k)} q_{km} + \sum_{mk\in \delta^-(k)} q_{mk} \end{aligned} \]
- Branch limits:
\[ \forall km \in B, \qquad \begin{aligned} |s_{km}| &\le L_{km}\\ |s_{mk}| &\le L_{km} \end{aligned} \]
Since \(L_{km}\) is positive, we implement these inequalities as follows:
\[ \forall km \in B, \qquad \begin{aligned} p_{km}^2+q_{km}^2 &\le L_{km}^2\\ p_{mk}^2+q_{mk}^2 &\le L_{km}^2 \end{aligned} \]
- Phase angle difference:
\[ \forall km \in B, \qquad \theta_{km}^\text{min} \le \theta_k - \theta_m \le \theta_{km}^\text{max} \]
- Reference angle. Every constraint above uses angle differences only, so the angles are fixed only up to a common shift. Pinning one bus makes the solution unique:
\[ \theta_1 = 0 \]
Input data
A Bus is a point of the grid, carrying the band its voltage magnitude has to stay inside and the active and reactive power drawn there. A Branch is a line or a transformer between two buses, carrying the three blocks of its admittance matrix, the limit on the power that may flow through it, and the transformer’s ratio and phase shift. A Generator sits at a bus with an active and a reactive output range and a polynomial cost. A Network holds the three lists, the base power the data was scaled against, and where to draw each bus.
A Dispatch carries its cost, the voltage magnitude and angle at every bus, the active and reactive power entering each branch at both of its ends, and what every generator produces.
The two cases come from Matpower and are read straight from their .m form under shared/data. A case file is a MATLAB script, but the parts we need are four numeric matrices and one scalar, so scraping the brackets out of the text is enough.
The values in the data set are given with a baseMVA of 100. It means we have to divide all power values (active and reactive power values, given in MW and MVAr in the data set) by 100 to have values in p.u.
We have to do the same thing for the costs, which are given in $/hr and MW (or MVAr). We must therefore multiply each component of the cost vectors by the corresponding power of 100.
It is already done with voltages (with a baseKV of 345 in this set). Angles are converted from degrees to radians.
import re
from pathlib import Path
DATA_DIR = Path("shared") / "data"
def _matrix(text, name):
opening = re.search(rf"mpc\.{name}\s*=\s*\[", text).end()
closing = text.index("]", opening)
rows = []
for line in text[opening:closing].split(";"):
stripped = re.sub("%.*", "", line).strip()
if stripped:
rows.append([float(field) for field in stripped.split()])
return rows
def _scalar(text, name):
return float(re.search(rf"mpc\.{name}\s*=\s*([0-9.eE+-]+)", text).group(1))
def read_matpower(name, positions, label):
text = (DATA_DIR / name).read_text(encoding="utf-8")
base = _scalar(text, "baseMVA")
buses = [
Bus(int(row[0]), row[12], row[11], row[2] / base, row[3] / base)
for row in _matrix(text, "bus")
]
branches = []
for row in _matrix(text, "branch"):
r, x, b = row[2], row[3], row[4]
y = complex(r / (r**2 + x**2), -x / (r**2 + x**2))
shunt = complex(0.0, b / 2)
tau = row[8] if row[8] != 0 else 1.0
branches.append(
Branch(
int(row[0]),
int(row[1]),
(y + shunt) / tau**2,
-y / tau,
y + shunt,
row[5] / base,
tau,
math.radians(row[9]),
math.radians(row[11]),
math.radians(row[12]),
)
)
gen = _matrix(text, "gen")
gencost = _matrix(text, "gencost")
degree = int(gencost[0][3]) - 1
generators = [
Generator(
int(gen[i][0]),
gen[i][9] / base,
gen[i][8] / base,
gen[i][4] / base,
gen[i][3] / base,
[gencost[i][4 + j] * base ** (degree - j) for j in range(degree + 1)],
)
for i in range(len(gen))
]
return Network(buses, branches, generators, base, positions, label)The 9-bus case of Chow (1982) is the smallest one worth drawing, with three generators, six load and tie buses, and nine branches of which three are transformers. The bus coordinates are not in the case file, so they are given here.
case9 = read_matpower(
"case9.m",
[
(42, 48),
(72, 8),
(12, 8),
(42, 38),
(22, 30),
(22, 18),
(42, 12),
(62, 18),
(62, 30),
],
"9-bus network",
)Every branch is drawn twice, as a pale band as wide as the limit on the power it may carry and a darker core inside it as wide as the power actually flowing, so the room a branch has left is the gap between the two. The red arrowhead sitting on a branch points the way the active power goes, and it is the only mark in the figure that means a direction. A generator’s disc grows with what it produces, and the short dark stub on a bus is the load drawn there, with its size in MVA beside it.
Show the outage helper
def without_branches(network, removed):
kept = [
branch
for branch in network.branches
if (branch.from_bus, branch.to_bus) not in removed
]
return Network(
network.buses,
kept,
network.generators,
network.base_power,
network.positions,
network.name,
)Both cases hold every voltage inside a narrow band, and the generators are the only places power enters the grid, so the layout below is really the map of how far that power has to travel.
report_network(case9)9-bus network: 9 buses, 9 branches, 3 generators
3 of the buses draw 315 MW and 115 MVAr between them
generator bus active reactive at full
(MW) (MVAr) ($/MWh)
--------------------------------------------------------------
0 1 10 to 250 -300 to 300 33.10
1 2 10 to 300 -300 to 300 28.70
2 3 10 to 270 -300 to 300 35.32
Branch limits run from 150 to 300 MVA
Voltage magnitudes are held between 0.90 and 1.10 per unitdraw_layout(case9)Model implementation
import knitrominimize_cost(network) builds the model, solves it with Knitro, and returns a Dispatch. The flow balance, the generation limits and the angle differences are affine, and the branch limits are quadratic. The four power flow equations are the nonconvex part. Each one multiplies two voltage magnitudes by the cosine or the sine of an angle difference, and add_constraint takes them written that way.
The angles enter every constraint as differences, so pinning the first bus to zero is what makes the solution unique rather than a whole circle of them.
def minimize_cost(network):
buses, branches, generators = network.buses, network.branches, network.generators
index = network.index
ends = [(index[branch.from_bus], index[branch.to_bus]) for branch in branches]
sites = [index[generator.bus] for generator in generators]
prob = knitro.Problem()
v = [prob.add_variable(lb=bus.min_voltage, ub=bus.max_voltage) for bus in buses]
theta = [prob.add_variable() for _ in buses]
p_km = [prob.add_variable() for _ in branches]
q_km = [prob.add_variable() for _ in branches]
p_mk = [prob.add_variable() for _ in branches]
q_mk = [prob.add_variable() for _ in branches]
p_g = [
prob.add_variable(lb=generator.min_active_power, ub=generator.max_active_power)
for generator in generators
]
q_g = [
prob.add_variable(
lb=generator.min_reactive_power, ub=generator.max_reactive_power
)
for generator in generators
]
prob.add_objective(
prob.nsum(
c * p_g[i] ** j
for i, generator in enumerate(generators)
for j, c in enumerate(reversed(generator.cost))
)
)
prob.add_constraint(theta[0] == 0.0)
for km, branch in enumerate(branches):
k, m = ends[km]
G_kk, B_kk = branch.y_kk.real, branch.y_kk.imag
G_km, B_km = branch.y_km.real, branch.y_km.imag
G_mm, B_mm = branch.y_mm.real, branch.y_mm.imag
shifted = theta[k] - theta[m] - branch.phase_shift
v_k, v_m = v[k], v[m]
cos_km = v_k * v_m * prob.cos(shifted)
sin_km = v_k * v_m * prob.sin(shifted)
prob.add_constraint(p_km[km] == G_kk * v_k**2 + G_km * cos_km + B_km * sin_km)
prob.add_constraint(q_km[km] == -B_kk * v_k**2 - B_km * cos_km + G_km * sin_km)
prob.add_constraint(p_mk[km] == G_mm * v_m**2 + G_km * cos_km - B_km * sin_km)
prob.add_constraint(q_mk[km] == -B_mm * v_m**2 - B_km * cos_km - G_km * sin_km)
if branch.max_flow > 0:
prob.add_constraint(p_km[km] ** 2 + q_km[km] ** 2 <= branch.max_flow**2)
prob.add_constraint(p_mk[km] ** 2 + q_mk[km] ** 2 <= branch.max_flow**2)
prob.add_constraint(theta[k] - theta[m] >= branch.min_angle)
prob.add_constraint(theta[k] - theta[m] <= branch.max_angle)
for k, bus in enumerate(buses):
supplying = [i for i, site in enumerate(sites) if site == k]
leaving = [km for km, (start, _) in enumerate(ends) if start == k]
arriving = [km for km, (_, stop) in enumerate(ends) if stop == k]
for generated, forward, backward, demand in (
(p_g, p_km, p_mk, bus.active_demand),
(q_g, q_km, q_mk, bus.reactive_demand),
):
supplied = prob.nsum(generated[i] for i in supplying)
leaving_flow = prob.nsum(forward[km] for km in leaving)
arriving_flow = prob.nsum(backward[km] for km in arriving)
prob.add_constraint(supplied - demand == leaving_flow + arriving_flow)
prob.solve()
variables = (v, theta, p_km, q_km, p_mk, q_mk, p_g, q_g)
values = [[x.value for x in group] for group in variables]
return Dispatch(prob.get_attr(knitro.KN_ATTR_OBJ_VALUE), *values)Output visualization
Solving picks a voltage for every bus and an output for every generator. The figure below draws what that dispatch does to the grid. Each branch carries a pale band as wide as its limit and a darker core as wide as the power crossing it, so a branch running out of room is one whose core has filled its band.
dispatch_9 = minimize_cost(case9)=======================================
Commercial License
Artelys Knitro 16.0.0
=======================================
Knitro using 1 thread.
No start point provided -- Knitro computing one.
Knitro presolve eliminated 1 variable (2%) and 3 constraints (3%) in 0.00s.
concurrent_evals 0
feastol 1e-06
feastol_abs 0.001
opttol 1e-06
opttol_abs 0.001
Knitro shifted start point further inside presolved bounds (4 variables).
Knitro fixing 1 variable eliminated from the presolve.
Problem Characteristics | Presolved
-----------------------
Problem type: NLP
Objective: minimize / quadratic
Number of variables: 60 | 59
bounds: lower upper range | lower upper range
0 0 15 | 0 0 16
free fixed | free fixed
45 0 | 43 0
Number of expressions: 427 | 72
Number of constraints: 91 | 88
eq. ineq. range | eq. ineq. range
linear: 19 18 0 | 18 16 0
quadratic: 0 18 0 | 0 18 0
nonlinear: 36 0 0 | 36 0 0
Number of nonzeros:
objective Jacobian Hessian | objective Jacobian Hessian
linear: 3 115 | 3 110
quadratic: 3 66 48 | 3 66 48
nonlinear: 0 144 54 | 0 144 50
total: 3 295 102 | 3 286 98
Coefficient range:
linear objective: [1e+02, 5e+02] | [2e+00, 1e+01]
linear constraints: [1e+00, 1e+00] | [1e+00, 1e+00]
quadratic objective: [9e+02, 1e+03] | [2e+01, 2e+01]
quadratic constraints: [1e+00, 2e+01] | [1e+00, 2e+01]
variable bounds: [1e-01, 3e+00] | [1e-01, 6e+00]
constraint bounds: [2e+00, 9e+00] | [3e-01, 9e+00]
Knitro using the Interior-Point/Barrier Direct algorithm.
Iter Objective FeasError OptError ||Step|| Time
-------- -------------- --------- --------- --------- --------
0 1.330766e+04 1.67e+00
5 5.296682e+03 1.03e-05 1.08e-05 2.18e-02 0.02
EXIT: Locally optimal solution found.
Final Statistics
----------------
Final objective value = 5.29668208813826e+03
Final feasibility error (abs / rel) = 1.03e-05 / 8.52e-07
Final optimality error (abs / rel) = 1.08e-05 / 2.19e-07
# of iterations = 5
# of CG iterations = 0
# of function evaluations = 11
# of gradient evaluations = 8
# of Hessian evaluations = 5
Total program time (secs) = 0.02059 ( 0.025 CPU time)
Time spent in evaluations (secs) = 0.00029
================================================================================
Show the full outputHide the full output
draw_dispatch(case9, dispatch_9)report_dispatch(case9, dispatch_9)9-bus network: $5,296.68 per hour
generator bus active reactive of range cost share
(MW) (MVAr) ($/h)
------------------------------------------------------------------
0 1 89.8 12.9 33% 1,486.0 28.1%
1 2 134.3 0.0 43% 2,294.8 43.3%
2 3 94.2 -22.6 32% 1,515.9 28.6%
branch carrying limit loading
(MVA) (MVA)
----------------------------------------
8 to 2 134.6 250 53.9%
5 to 6 57.3 150 38.2%
1 to 4 90.7 250 36.3%
3 to 6 96.9 300 32.3%
8 to 9 72.8 250 29.1%
6 to 7 38.6 150 25.7%
7 to 8 64.0 250 25.6%
9 to 4 62.5 250 25.0%
4 to 5 35.4 250 14.2%
318.3 MW generated against 315.0 MW of demand, 3.3 MW lost in the branches
Voltages run from 1.0717 pu at bus 9 to 1.1000 pu at bus 8No generator here runs flat out, and the cheapest is not simply used first. The cost curves are quadratic, so what an economic dispatch equalises is the marginal cost of the next megawatt, and the three meet it at nearly the same figure. Nothing on the network binds either, the busiest branch running at 54 percent of its rating, so this dispatch is settled by the cost curves rather than by the grid. Where the solution does press against a bound is the voltages, and they ride at the top of their band. The same power at a higher voltage is a lower current, and current is what the branches lose to heat.
A larger network
Now, let’s solve another, larger, use case. The New England system is a standard test network an order of magnitude bigger, and eleven of its branches are transformers with an off-nominal ratio. The case file carries no coordinates, so they are given here.
Show the bus coordinates
CASE39_POSITIONS = [
(29.2, 19.6),
(31.6, 28.1),
(42.4, 30.1),
(52.3, 26.8),
(54.5, 17.9),
(63.2, 12.4),
(57.1, 5.7),
(49.0, 10.4),
(39.4, 7.5),
(80.3, 25.1),
(73.0, 15.4),
(69.7, 22.3),
(70.9, 29.5),
(60.9, 32.8),
(59.9, 41.7),
(55.0, 49.9),
(44.9, 46.6),
(44.0, 38.4),
(50.7, 58.2),
(51.2, 67.0),
(60.4, 57.4),
(68.7, 61.8),
(74.1, 56.3),
(68.5, 48.7),
(24.3, 34.7),
(24.6, 43.8),
(34.2, 47.6),
(11.3, 42.8),
(20.4, 55.5),
(22.9, 24.4),
(68.3, 5.0),
(88.0, 21.7),
(43.2, 63.2),
(49.0, 75.0),
(73.0, 69.2),
(83.4, 59.2),
(14.2, 32.1),
(8.0, 53.3),
(30.5, 11.3),
]case39 = read_matpower("case39.m", CASE39_POSITIONS, "39-bus network")report_network(case39)39-bus network: 39 buses, 46 branches, 10 generators
21 of the buses draw 6,254 MW and 1,387 MVAr between them
generator bus active reactive at full
(MW) (MVAr) ($/MWh)
--------------------------------------------------------------
0 30 0 to 1,040 140 to 400 10.70
1 31 0 to 646 -100 to 300 6.76
2 32 0 to 725 150 to 300 7.55
3 33 0 to 652 0 to 250 6.82
4 34 0 to 508 0 to 167 5.38
5 35 0 to 687 -100 to 300 7.17
6 36 0 to 580 0 to 240 6.10
7 37 0 to 564 0 to 250 5.94
8 38 0 to 865 -150 to 300 8.95
9 39 0 to 1,100 -100 to 300 11.30
Branch limits run from 480 to 1,800 MVA
Voltage magnitudes are held between 0.94 and 1.06 per unitdraw_layout(case39)dispatch_39 = minimize_cost(case39)=======================================
Commercial License
Artelys Knitro 16.0.0
=======================================
Knitro using 1 thread.
No start point provided -- Knitro computing one.
Knitro presolve eliminated 1 variable (0%) and 5 constraints (1%) in 0.00s.
concurrent_evals 0
feastol 1e-06
feastol_abs 0.001
opttol 1e-06
opttol_abs 0.001
Knitro shifted start point further inside presolved bounds (15 variables).
Knitro fixing 1 variable eliminated from the presolve.
Problem Characteristics | Presolved
-----------------------
Problem type: NLP
Objective: minimize / quadratic
Number of variables: 282 | 281
bounds: lower upper range | lower upper range
0 0 59 | 0 0 61
free fixed | free fixed
223 0 | 220 0
Number of expressions: 2221 | 354
Number of constraints: 447 | 442
eq. ineq. range | eq. ineq. range
linear: 79 92 0 | 78 88 0
quadratic: 0 92 0 | 0 92 0
nonlinear: 184 0 0 | 184 0 0
Number of nonzeros:
objective Jacobian Hessian | objective Jacobian Hessian
linear: 10 573 | 10 564
quadratic: 10 360 233 | 10 360 233
nonlinear: 0 736 262 | 0 736 256
total: 10 1493 495 | 10 1476 489
Coefficient range:
linear objective: [3e+01, 3e+01] | [2e+00, 2e+00]
linear constraints: [1e+00, 1e+00] | [2e-01, 1e+00]
quadratic objective: [1e+02, 1e+02] | [8e+00, 8e+00]
quadratic constraints: [7e-01, 4e+02] | [7e-01, 1e+02]
variable bounds: [9e-01, 1e+01] | [9e-01, 1e+01]
constraint bounds: [6e+00, 3e+02] | [2e-02, 3e+02]
Knitro using the Interior-Point/Barrier Direct algorithm.
Iter Objective FeasError OptError ||Step|| Time
-------- -------------- --------- --------- --------- --------
0 1.610205e+04 2.33e+01
9 4.186418e+04 1.07e-08 2.92e-07 1.86e-03 0.05
EXIT: Locally optimal solution found.
Final Statistics
----------------
Final objective value = 4.18641777890036e+04
Final feasibility error (abs / rel) = 1.07e-08 / 8.04e-11
Final optimality error (abs / rel) = 2.92e-07 / 2.69e-09
# of iterations = 9
# of CG iterations = 0
# of function evaluations = 15
# of gradient evaluations = 12
# of Hessian evaluations = 9
Total program time (secs) = 0.05058 ( 0.039 CPU time)
Time spent in evaluations (secs) = 0.00171
================================================================================
Show the full outputHide the full output
draw_dispatch(case39, dispatch_39)report_dispatch(case39, dispatch_39)39-bus network: $41,864.18 per hour
generator bus active reactive of range cost share
(MW) (MVAr) ($/h)
------------------------------------------------------------------
0 30 671.6 140.0 65% 4,712.1 11.3%
1 31 646.0 300.0 100% 4,367.2 10.4%
2 32 671.2 300.0 93% 4,706.0 11.2%
3 33 652.0 115.1 100% 4,446.8 10.6%
4 34 508.0 139.6 100% 2,733.2 6.5%
5 35 661.5 222.9 96% 4,573.9 10.9%
6 36 580.0 60.6 100% 3,538.2 8.5%
7 37 564.0 8.6 100% 3,350.4 8.0%
8 38 654.0 -32.8 76% 4,474.0 10.7%
9 39 689.6 82.1 63% 4,962.4 11.9%
branch carrying limit loading
(MVA) (MVA)
----------------------------------------
2 to 3 454.9 500 91.0%
6 to 11 386.3 480 80.5%
16 to 19 472.1 600 78.7%
10 to 32 700.3 900 77.8%
22 to 35 680.7 900 75.6%
2 to 30 674.5 900 74.9%
19 to 33 651.2 900 72.4%
21 to 22 626.6 900 69.6%
The 38 branches carrying less than these are not listed
6,297.8 MW generated against 6,254.2 MW of demand, 43.6 MW lost in the branches
Voltages run from 0.9968 pu at bus 20 to 1.0600 pu at bus 26Switching two branches off
Let’s see how turning off a single branch can change the whole power flow through the network, because the other branches compensate to carry all the necessary power. The table above leaves the branch from 2 to 3 at the top of the loading column. It is switched off here together with the one from 21 to 22, and the same problem is solved on what is left.
OUTAGE = [(2, 3), (21, 22)]
case39_outage = without_branches(case39, OUTAGE)dispatch_outage = minimize_cost(case39_outage)=======================================
Commercial License
Artelys Knitro 16.0.0
=======================================
Knitro using 1 thread.
No start point provided -- Knitro computing one.
Knitro presolve eliminated 3 variables (1%) and 8 constraints (2%) in 0.00s.
concurrent_evals 0
feastol 1e-06
feastol_abs 0.001
opttol 1e-06
opttol_abs 0.001
Knitro shifted start point further inside presolved bounds (17 variables).
Knitro fixing 3 variables eliminated from the presolve.
Problem Characteristics | Presolved
-----------------------
Problem type: NLP
Objective: minimize / quadratic
Number of variables: 274 | 271
bounds: lower upper range | lower upper range
0 0 59 | 0 0 61
free fixed | free fixed
215 0 | 210 0
Number of expressions: 2135 | 342
Number of constraints: 431 | 423
eq. ineq. range | eq. ineq. range
linear: 79 88 0 | 76 84 0
quadratic: 0 88 0 | 0 87 0
nonlinear: 176 0 0 | 176 0 0
Number of nonzeros:
objective Jacobian Hessian | objective Jacobian Hessian
linear: 10 549 | 10 536
quadratic: 10 344 225 | 10 342 223
nonlinear: 0 704 254 | 0 704 248
total: 10 1429 479 | 10 1406 471
Coefficient range:
linear objective: [3e+01, 3e+01] | [2e+00, 2e+00]
linear constraints: [1e+00, 1e+00] | [2e-01, 1e+00]
quadratic objective: [1e+02, 1e+02] | [8e+00, 8e+00]
quadratic constraints: [7e-01, 4e+02] | [7e-01, 1e+02]
variable bounds: [9e-01, 1e+01] | [9e-01, 1e+01]
constraint bounds: [6e+00, 3e+02] | [2e-02, 3e+02]
Knitro using the Interior-Point/Barrier Direct algorithm.
Iter Objective FeasError OptError ||Step|| Time
-------- -------------- --------- --------- --------- --------
0 1.701602e+04 2.33e+01
10 4.537491e+04 1.12e-01 1.91e+01 2.34e+00 0.03
20 4.675140e+04 1.37e-03 8.27e-02 6.41e-03 0.05
30 4.674739e+04 4.71e-03 1.74e+02 4.09e-03 0.08
40 4.678173e+04 1.85e-02 1.74e+02 1.90e-04 0.10
50 4.669712e+04 1.14e-02 1.74e+02 4.60e-01 0.12
60 4.678175e+04 1.85e-02 1.74e+02 8.99e-03 0.13
70 4.672082e+04 9.25e-03 1.74e+02 2.55e-01 0.15
80 4.678085e+04 1.85e-02 1.74e+02 5.06e-04 0.17
90 4.670535e+04 9.46e-03 1.74e+02 1.31e+00 0.19
100 4.678080e+04 1.85e-02 1.74e+02 3.30e-02 0.21
103 4.678175e+04 1.85e-02 1.74e+02 3.66e-05 0.22
EXIT: Convergence to an infeasible point. Problem appears to be locally
infeasible. If problem is believed to be feasible, try multistart
to search for feasible points, or decrease infeastol.
Final Statistics
----------------
Final objective value = 4.67817488808525e+04
Final feasibility error (abs / rel) = 1.85e-02 / 1.40e-04
Final optimality error (abs / rel) = 1.74e+02 / 9.99e-01
# of iterations = 103
# of CG iterations = 24
# of function evaluations = 141
# of gradient evaluations = 118
# of Hessian evaluations = 104
Total program time (secs) = 0.21798 ( 0.235 CPU time)
Time spent in evaluations (secs) = 0.01833
================================================================================
Show the full outputHide the full output
draw_dispatch(case39_outage, dispatch_outage, removed=OUTAGE)report_dispatch(case39_outage, dispatch_outage)39-bus network: $46,781.75 per hour
generator bus active reactive of range cost share
(MW) (MVAr) ($/h)
------------------------------------------------------------------
0 30 869.5 140.0 84% 7,820.8 16.7%
1 31 646.0 300.0 100% 4,367.2 9.3%
2 32 725.0 300.0 100% 5,473.9 11.7%
3 33 652.0 250.0 100% 4,446.8 9.5%
4 34 508.0 167.0 100% 2,733.2 5.8%
5 35 490.7 184.7 71% 2,555.2 5.5%
6 36 349.7 72.2 60% 1,328.2 2.8%
7 37 345.1 0.0 61% 1,294.6 2.8%
8 38 643.3 60.2 74% 4,331.6 9.3%
9 39 1,100.0 213.3 100% 12,430.2 26.6%
branch carrying limit loading
(MVA) (MVA)
----------------------------------------
26 to 27 600.0 600 100.0%
23 to 24 600.0 600 100.0%
25 to 26 600.0 600 100.0%
2 to 25 500.0 500 100.0%
2 to 30 869.5 900 96.6%
22 to 23 513.0 600 85.5%
16 to 19 501.5 600 83.6%
10 to 32 746.1 900 82.9%
The 36 branches carrying less than these are not listed
6,329.3 MW generated against 6,254.2 MW of demand, 75.0 MW lost in the branches
Voltages run from 0.9819 pu at bus 37 to 1.0600 pu at bus 39The power those two were carrying has to go somewhere, and the dotted lines show where it no longer can. Their neighbours take it up, four of them ending the solve pinned at exactly their limits, and the branch from 2 to 25 goes from a third of its rating to all of it. What the branches nearby cannot carry is made up further away and at a higher price. The generator at bus 39 runs up to its ceiling while the two at buses 36 and 37, both flat out before, back off by about 40 percent. Both dispatches are locally optimal points of a nonconvex problem, so read the gap between the two costs as an indication rather than a measurement.
References
- “Resolution of optimization problems for large-scale power transmission networks with semidefenite programming” - J. Sliwak (2021) HAL
- The two case files come from Matpower. The 9-bus case is based on data from J. H. Chow, ed., Time-Scale Modeling of Dynamic Networks with Applications to Power Systems, Springer-Verlag (1982), after R. P. Schulz, A. E. Turner and D. N. Ewart, “Long Term Power System Dynamics”, EPRI Report 90-7-0 (1974). The 39-bus case is the New England system of G. Bills et al., “On-Line Stability Analysis Study”, EPRI Report 90-7-0 (1970).