struct Bus
number::Int
min_voltage::Float64
max_voltage::Float64
active_demand::Float64
reactive_demand::Float64
end
struct Branch
from_bus::Int
to_bus::Int
y_kk::ComplexF64
y_km::ComplexF64
y_mm::ComplexF64
max_flow::Float64
ratio::Float64
phase_shift::Float64
min_angle::Float64
max_angle::Float64
end
struct Generator
bus::Int
min_active_power::Float64
max_active_power::Float64
min_reactive_power::Float64
max_reactive_power::Float64
cost::Vector{Float64}
end
struct Network
buses::Vector{Bus}
branches::Vector{Branch}
generators::Vector{Generator}
base_power::Float64
positions::Vector{Tuple{Float64,Float64}}
name::String
index::Dict{Int,Int}
end
function Network(buses, branches, generators, base_power, positions, name)
return Network(
buses,
branches,
generators,
base_power,
positions,
name,
Dict(bus.number => k for (k, bus) in enumerate(buses)),
)
end
num_buses(network::Network) = length(network.buses)
num_branches(network::Network) = length(network.branches)
num_generators(network::Network) = length(network.generators)
struct Dispatch
cost::Float64
v::Vector{Float64}
theta::Vector{Float64}
p_km::Vector{Float64}
q_km::Vector{Float64}
p_mk::Vector{Float64}
q_mk::Vector{Float64}
p_g::Vector{Float64}
q_g::Vector{Float64}
endAC 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.
const DATA_DIR = joinpath("shared", "data")
function _matrix(text, name)
opening = last(findfirst(Regex("mpc\\.$name\\s*=\\s*\\["), text))
closing = findnext(']', text, opening)
rows = Vector{Float64}[]
for line in split(text[(opening + 1):(closing - 1)], ';')
stripped = strip(replace(line, r"%.*" => ""))
isempty(stripped) && continue
push!(rows, [parse(Float64, field) for field in split(stripped)])
end
return rows
end
function _scalar(text, name)
found = match(Regex("mpc\\.$name\\s*=\\s*([0-9.eE+-]+)"), text)
return parse(Float64, found.captures[1])
end
function read_matpower(name, positions, label)
text = read(joinpath(DATA_DIR, name), String)
base = _scalar(text, "baseMVA")
buses = [
Bus(Int(row[1]), row[13], row[12], row[3] / base, row[4] / base) for
row in _matrix(text, "bus")
]
branches = Branch[]
for row in _matrix(text, "branch")
r, x, b = row[3], row[4], row[5]
y = complex(r / (r^2 + x^2), -x / (r^2 + x^2))
shunt = complex(0.0, b / 2)
tau = row[9] == 0 ? 1.0 : row[9]
push!(
branches,
Branch(
Int(row[1]),
Int(row[2]),
(y + shunt) / tau^2,
-y / tau,
y + shunt,
row[6] / base,
tau,
deg2rad(row[10]),
deg2rad(row[12]),
deg2rad(row[13]),
),
)
end
gen = _matrix(text, "gen")
gencost = _matrix(text, "gencost")
degree = Int(gencost[1][4]) - 1
generators = [
Generator(
Int(gen[i][1]),
gen[i][10] / base,
gen[i][9] / base,
gen[i][5] / base,
gen[i][4] / base,
[gencost[i][5 + j] * base^(degree - j) for j in 0:degree],
) for i in eachindex(gen)
]
return Network(buses, branches, generators, base, positions, label)
endThe 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.0, 48.0),
(72.0, 8.0),
(12.0, 8.0),
(42.0, 38.0),
(22.0, 30.0),
(22.0, 18.0),
(42.0, 12.0),
(62.0, 18.0),
(62.0, 30.0),
],
"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
function without_branches(network::Network, removed)
kept = [
branch for
branch in network.branches if !((branch.from_bus, branch.to_bus) in removed)
]
return Network(
network.buses,
kept,
network.generators,
network.base_power,
network.positions,
network.name,
)
endBoth 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
using JuMP
using 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 the same @constraint macro takes all of them.
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.
function minimize_cost(network::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]
model = Model(KNITRO.Optimizer)
@variable(
model,
buses[k].min_voltage <= v[k in eachindex(buses)] <= buses[k].max_voltage,
start = 1.0
)
@variable(model, theta[eachindex(buses)], start = 0.0)
@variable(model, p_km[eachindex(branches)])
@variable(model, q_km[eachindex(branches)])
@variable(model, p_mk[eachindex(branches)])
@variable(model, q_mk[eachindex(branches)])
@variable(
model,
generators[i].min_active_power <=
p_g[i in eachindex(generators)] <=
generators[i].max_active_power
)
@variable(
model,
generators[i].min_reactive_power <=
q_g[i in eachindex(generators)] <=
generators[i].max_reactive_power
)
@objective(
model,
Min,
sum(
c * p_g[i]^(j - 1) for (i, generator) in enumerate(generators) for
(j, c) in enumerate(reverse(generator.cost))
)
)
@constraint(model, theta[1] == 0)
for (km, branch) in enumerate(branches)
k, m = ends[km]
G_kk, B_kk = real(branch.y_kk), imag(branch.y_kk)
G_km, B_km = real(branch.y_km), imag(branch.y_km)
G_mm, B_mm = real(branch.y_mm), imag(branch.y_mm)
shifted = theta[k] - theta[m] - branch.phase_shift
cos_km = v[k] * v[m] * cos(shifted)
sin_km = v[k] * v[m] * sin(shifted)
@constraint(model, p_km[km] == G_kk * v[k]^2 + G_km * cos_km + B_km * sin_km)
@constraint(model, q_km[km] == -B_kk * v[k]^2 - B_km * cos_km + G_km * sin_km)
@constraint(model, p_mk[km] == G_mm * v[m]^2 + G_km * cos_km - B_km * sin_km)
@constraint(model, q_mk[km] == -B_mm * v[m]^2 - B_km * cos_km - G_km * sin_km)
if branch.max_flow > 0
@constraint(model, p_km[km]^2 + q_km[km]^2 <= branch.max_flow^2)
@constraint(model, p_mk[km]^2 + q_mk[km]^2 <= branch.max_flow^2)
end
@constraint(model, branch.min_angle <= theta[k] - theta[m] <= branch.max_angle)
end
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 = @expression(model, sum(generated[i] for i in supplying))
leaving_flow = @expression(model, sum(forward[km] for km in leaving))
arriving_flow = @expression(model, sum(backward[km] for km in arriving))
@constraint(model, supplied - demand == leaving_flow + arriving_flow)
end
end
optimize!(model)
variables = (v, theta, p_km, q_km, p_mk, q_mk, p_g, q_g)
return Dispatch(objective_value(model), (value.(x) for x in variables)...)
endOutput 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.
Knitro presolve eliminated 1 variable (2%) and 3 constraints (3%) in 0.00s.
datacheck 0
feastol 1e-06
feastol_abs 0.001
hessian_no_f 1
mip_numthreads 1
ms_numthreads 1
numthreads 1
opttol 1e-06
opttol_abs 0.001
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 constraints: 82 | 88
eq. ineq. range | eq. ineq. range
linear: 19 0 9 | 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 61 | 3 74
quadratic: 3 36 39 | 3 36 39
nonlinear: 0 180 63 | 0 180 59
total: 3 277 102 | 3 286 98
Coefficient range:
linear objective: [1e+02, 5e+02] | [3e+00, 2e+01]
linear constraints: [1e+00, 1e+00] | [1e+00, 1e+00]
quadratic objective: [9e+02, 1e+03] | [3e+01, 4e+01]
quadratic constraints: [1e+00, 1e+00] | [1e+00, 1e+00]
variable bounds: [1e-01, 3e+00] | [1e-01, 6e+00]
constraint bounds: [3e-01, 9e+00] | [3e-01, 9e+00]
Knitro using the Interior-Point/Barrier Direct algorithm.
Iter Objective FeasError OptError ||Step|| Time
-------- -------------- --------- --------- --------- --------
0 5.642018e+03 1.25e+00
5 5.296686e+03 3.35e-08 1.29e-05 1.01e-04 14.62
EXIT: Locally optimal solution found.
Final Statistics
----------------
Final objective value = 5.29668619828880e+03
Final feasibility error (abs / rel) = 3.35e-08 / 2.68e-08
Final optimality error (abs / rel) = 1.29e-05 / 1.51e-07
# of iterations = 5
# of CG iterations = 0
# of function evaluations = 10
# of gradient evaluations = 8
# of Hessian evaluations = 5
Total program time (secs) = 14.62381 ( 14.613 CPU time)
Time spent in evaluations (secs) = 11.64504
================================================================================
Show the full outputHide the full output
draw_dispatch(case9, dispatch_9)report_dispatch(case9, dispatch_9)9-bus network: $5,296.69 per hour
generator bus active reactive of range cost share
(MW) (MVAr) ($/h)
------------------------------------------------------------------
0 1 89.8 13.0 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.0718 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
const 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.
Knitro presolve eliminated 1 variable (0%) and 5 constraints (1%) in 0.00s.
datacheck 0
feastol 1e-06
feastol_abs 0.001
hessian_no_f 1
mip_numthreads 1
ms_numthreads 1
numthreads 1
opttol 1e-06
opttol_abs 0.001
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 constraints: 401 | 442
eq. ineq. range | eq. ineq. range
linear: 79 0 46 | 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 297 | 10 380
quadratic: 10 184 194 | 10 184 194
nonlinear: 0 920 301 | 0 920 295
total: 10 1401 495 | 10 1476 489
Coefficient range:
linear objective: [3e+01, 3e+01] | [1e+01, 1e+01]
linear constraints: [1e+00, 1e+00] | [1e+00, 1e+00]
quadratic objective: [1e+02, 1e+02] | [4e+01, 4e+01]
quadratic constraints: [1e+00, 1e+00] | [1e+00, 1e+00]
variable bounds: [9e-01, 1e+01] | [9e-01, 1e+01]
constraint bounds: [2e-02, 3e+02] | [2e-02, 3e+02]
Knitro using the Interior-Point/Barrier Direct algorithm.
Iter Objective FeasError OptError ||Step|| Time
-------- -------------- --------- --------- --------- --------
0 1.279100e+03 1.00e+01
8 4.186418e+04 3.76e-08 6.06e-04 1.67e-03 0.02
EXIT: Locally optimal solution found.
Final Statistics
----------------
Final objective value = 4.18641777354737e+04
Final feasibility error (abs / rel) = 3.76e-08 / 3.40e-09
Final optimality error (abs / rel) = 6.06e-04 / 9.80e-07
# of iterations = 8
# of CG iterations = 0
# of function evaluations = 13
# of gradient evaluations = 11
# of Hessian evaluations = 8
Total program time (secs) = 0.01757 ( 0.018 CPU time)
Time spent in evaluations (secs) = 0.00598
================================================================================
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 19Switching 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.
const 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.
Knitro presolve eliminated 3 variables (1%) and 8 constraints (2%) in 0.00s.
datacheck 0
feastol 1e-06
feastol_abs 0.001
hessian_no_f 1
mip_numthreads 1
ms_numthreads 1
numthreads 1
opttol 1e-06
opttol_abs 0.001
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 constraints: 387 | 423
eq. ineq. range | eq. ineq. range
linear: 79 0 44 | 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 285 | 10 362
quadratic: 10 176 186 | 10 174 184
nonlinear: 0 880 293 | 0 880 287
total: 10 1341 479 | 10 1406 471
Coefficient range:
linear objective: [3e+01, 3e+01] | [1e+01, 1e+01]
linear constraints: [1e+00, 1e+00] | [1e+00, 1e+00]
quadratic objective: [1e+02, 1e+02] | [4e+01, 4e+01]
quadratic constraints: [1e+00, 1e+00] | [1e+00, 1e+00]
variable bounds: [9e-01, 1e+01] | [9e-01, 1e+01]
constraint bounds: [2e-02, 3e+02] | [2e-02, 3e+02]
Knitro using the Interior-Point/Barrier Direct algorithm.
Iter Objective FeasError OptError ||Step|| Time
-------- -------------- --------- --------- --------- --------
0 1.279100e+03 1.00e+01
10 4.191498e+04 4.97e-01 7.24e+01 1.80e+00 0.02
20 4.638528e+04 3.29e-02 1.37e+02 2.46e-01 0.04
30 4.668746e+04 1.34e-03 5.58e-01 4.99e-03 0.06
40 4.718586e+04 1.13e-02 1.70e+08 4.57e+00 0.08
50 4.666578e+04 6.68e-03 9.78e+02 4.17e-01 0.10
60 4.670029e+04 1.72e-02 9.78e+02 6.57e-01 0.12
70 4.677971e+04 1.85e-02 9.78e+02 7.16e-07 0.14
80 4.652233e+04 1.71e-02 9.77e+02 3.49e-02 0.16
90 4.670337e+04 1.76e-02 9.78e+02 6.46e-01 0.18
100 4.678074e+04 1.85e-02 9.78e+02 3.64e-04 0.20
110 4.652233e+04 1.71e-02 9.77e+02 3.51e-02 0.22
120 4.670337e+04 1.76e-02 9.78e+02 6.46e-01 0.24
130 4.678074e+04 1.85e-02 9.78e+02 3.64e-04 0.25
140 4.652233e+04 1.71e-02 9.77e+02 3.51e-02 0.28
150 4.670337e+04 1.76e-02 9.78e+02 6.46e-01 0.29
160 4.678074e+04 1.85e-02 9.78e+02 3.64e-04 0.31
162 4.678172e+04 1.85e-02 9.78e+02 3.08e-05 0.32
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.67817223142079e+04
Final feasibility error (abs / rel) = 1.85e-02 / 1.67e-03
Final optimality error (abs / rel) = 9.78e+02 / 1.00e+00
# of iterations = 162
# of CG iterations = 5
# of function evaluations = 221
# of gradient evaluations = 177
# of Hessian evaluations = 163
Total program time (secs) = 0.31623 ( 0.316 CPU time)
Time spent in evaluations (secs) = 0.10902
================================================================================
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.72 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.1 5.5%
6 36 349.7 72.2 60% 1,328.3 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).