Power Flow Problem#
Overview#
A typical nonlinear equation is the power flow equations. Power flow equations are algebraic equations enforcing the power balancing at each bus. Intuitively, the equations require that, at each bus, the power injection (generation minus load) is equal to the power that leaves the bus through connected lines.
Using the complex power definition, the power flow equations are as simple as:
where \(P_i\) and \(Q_i\) are the real and reactive power injection at bus \(i\), \(V_i\) and \(V_j\) are the voltage at bus \(i\) and \(j\), \(Y\) is the admittance matrix, and
If bus k has no injection, the corresponding \(P_k\) and \(Q_k\) are zero.
For a system with \(N\) buses, there are \(2N\) real-valued equations. There has to be \(2N\) unkonwns to solve the equations.
As we have learned from power system analysis, there are three types of buses, PQ, PV and Slack. The unknowns are:
Bus Type |
Knowns |
Unknowns |
|---|---|---|
PQ |
\(P_i, Q_i\) |
\(V_i, \theta_i\) |
PV |
\(P_i, V_i\) |
\(Q_i, \theta_i\) |
Slack |
\(V_i, \theta_i\) |
\(P_i, Q_i\) |
Example: 3-bus System#
We consider a three-bus system from [1] (Example 3.11). I have taken the liberty to use a 0-based index to align with Python code.
It is a three bus sytem with the following bus parameters:
Bus |
Type |
V |
Pgen |
Qgen |
Pload |
Qload |
|---|---|---|---|---|---|---|
0 |
Slack |
1.02 |
– |
– |
0.0 |
0.0 |
1 |
PV |
1.00 |
0.5 |
– |
0.0 |
0.0 |
2 |
PQ |
– |
0 |
0 |
1.2 |
0.5 |
And the line parameters are given as
From |
To |
R |
X |
B |
|---|---|---|---|---|
0 |
1 |
0.02 |
0.3 |
0.15 |
0 |
2 |
0.01 |
0.1 |
0.1 |
1 |
2 |
0.01 |
0.1 |
0.1 |
First, let’s make the admittance matrix.
import numpy as np
from scipy.optimize import fsolve
Y = np.zeros((3, 3), dtype=np.complex128)
# Diagonal elements
Y[0, 0] = 1 / (0.02 + 0.3j) + 1 / (0.01 + 0.1j) + 1j * 0.15 / 2 + 1j * 0.1 / 2
Y[1, 1] = 1 / (0.02 + 0.3j) + 1 / (0.01 + 0.1j) + 1j * 0.15 / 2 + 1j * 0.1 / 2
Y[2, 2] = 1 / (0.01 + 0.1j) + 1 / (0.01 + 0.1j) + 1j * 0.1 / 2 + 1j * 0.1 / 2
# Off-diagonal elements
Y[0, 1] -= 1 / (0.02 + 0.3j)
Y[1, 0] -= 1 / (0.02 + 0.3j)
Y[0, 2] -= 1 / (0.01 + 0.1j)
Y[2, 0] -= 1 / (0.01 + 0.1j)
Y[1, 2] -= 1 / (0.01 + 0.1j)
Y[2, 1] -= 1 / (0.01 + 0.1j)
print(Y)
[[ 1.21133795-13.09457417j -0.22123894 +3.31858407j
-0.99009901 +9.9009901j ]
[-0.22123894 +3.31858407j 1.21133795-13.09457417j
-0.99009901 +9.9009901j ]
[-0.99009901 +9.9009901j -0.99009901 +9.9009901j
1.98019802-19.7019802j ]]
We can verify that the result matches the textbook:
# print Y in a polar form (using degrees)
print("Y magnitudes:")
print(np.abs(Y))
print("Y angles:")
print(np.angle(Y, deg=True))
Y magnitudes:
[[13.15048335 3.32595053 9.9503719 ]
[ 3.32595053 13.15048335 9.9503719 ]
[ 9.9503719 9.9503719 19.80124259]]
Y angles:
[[-84.71478916 93.81407483 95.71059314]
[ 93.81407483 -84.71478916 95.71059314]
[ 95.71059314 95.71059314 -84.26061502]]
We can move on and write the power flow equations. To conform with the standard form of f(x) = 0, the power injections are subtracted from both sides.
TODO: briefly explain the signs
def pf_3bus(x):
# we can make the unknowns to be
# x = [P0, Q0, theta1, Q1, theta2, V2]
# where theta is in radians and V is in per unit
V = np.array([1.02, 1.00, x[5]])
theta = np.array([0.0, x[2], x[4]])
Vc = V * np.exp(1j * theta)
# calculate the power injection into the network
# note that power leaves the bus into the network
S = np.diag(Vc) @ np.conj(Y @ Vc)
# power leaving each bus via the lines
# minus power injection at each bus shall equal to 0
P = np.real(S) - np.array([x[0], 0.5, -1.2])
Q = np.imag(S) - np.array([x[1], x[3], -0.5])
return np.concatenate([P, Q])
initial_guess = [0, 0, 0, 0, 0, 1.0]
sol, infodict, ier, mesg = fsolve(pf_3bus, initial_guess, full_output=True)
print(sol)
[ 0.70867913 0.28056964 -0.01009172 -0.04462179 -0.06351473 0.98158465]
This result is exactly the same as the one in the textbook.
Notes on the formulation#
If you still have some impressions from the Power System Analysis course, you may remember the lengthy derivations of the power flow equations and its Jacobian matrix. In fact, this formulation comes much simpler than the textbook ones, so it can be interesting to understand why both formulations work.
You may also have many questions about the formulation, such as
Why is the P&Q at a slack bus and the Q at a PV bus included as variables, resulting in six variables for three buses?
Following up on the previous question, why do we include two equations at the slack bus and one reactive power equation at the PV bus?
Why are the variables not ordered in the typical way, namely, unknown phase angles and then unknown voltage magnitudes, yet it solves?
Following up on the previous question, why is that ordering the convention (as written in textbooks and implemented in many programs)? Is that still relevant?
These are left as exercises.
Larger Systems#
The implementation above is elegent due to the very few lines of code. However, without providing the Jacobian matrix, computational speed challenges are expected.
We will use the ANDES software to load the IEEE 300-bus system.
import urllib.request
import os
import tempfile
# URL of the file
url = "https://raw.githubusercontent.com/MATPOWER/matpower/refs/heads/master/data/case14.m"
# Create a temporary file
with tempfile.NamedTemporaryFile(delete=False, suffix='.m') as temp_file:
# Download the file
urllib.request.urlretrieve(url, temp_file.name)
print(f"File downloaded and saved to {temp_file.name}")
# Verify if the file exists
if os.path.exists(temp_file.name):
print(f"File size: {os.path.getsize(temp_file.name)} bytes")
else:
print("Download failed")
File downloaded and saved to /var/folders/__/n5kx_m_s0tbg6n5qd7rh51700000gn/T/tmp61i9hwa5.m
File size: 4597 bytes
import andes
from kvxopt import matrix
ss = andes.load(temp_file.name)
Y = matrix(ss.Line.build_y())
nbus = ss.Bus.n
P = np.zeros(nbus)
Q = np.zeros(nbus)
V = np.ones(nbus)
theta = np.zeros(nbus)
Generating code for 1 models on 12 processes.
pq_loc = ss.Bus.idx2uid(ss.PQ.bus.v)
pv_loc = ss.Bus.idx2uid(ss.PV.bus.v)
slack_loc = ss.Bus.idx2uid(ss.Slack.bus.v)
gen_loc = np.concatenate([pv_loc, slack_loc])
non_gen_loc = np.setdiff1d(np.arange(ss.Bus.n), gen_loc)
non_slack_loc = np.setdiff1d(np.arange(ss.Bus.n), slack_loc)
shunt_loc = ss.Bus.idx2uid(ss.Shunt.bus.v)
def pf_large(x, ss, P, Q, V, theta):
# PQ
P[pq_loc] -= ss.PQ.p0.v
Q[pq_loc] -= ss.PQ.q0.v
# PV
V[pv_loc] = ss.PV.v0.v
P[pv_loc] += ss.PV.p0.v
# Slack
theta[slack_loc] = ss.Slack.a0.v
V[slack_loc] = ss.Slack.v0.v
# retrieve unknowns from `x`
# unknowns in `x` are grouped by theta and then V
# note that the `non_gen_loc` can have random order
# no sorting is needed, because x will be in the same order
theta[non_slack_loc] = x[:len(non_slack_loc)]
V[non_gen_loc] = x[len(non_slack_loc):]
# shunt elements
P[shunt_loc] += ss.Shunt.g.v * V[shunt_loc] ** 2
Q[shunt_loc] += ss.Shunt.b.v * V[shunt_loc] ** 2
Vc = V * np.exp(1j * theta)
# calculate the power injection into the network
# note that power leaves the bus into the network
S = np.diag(Vc) @ np.conj(Y @ Vc)
# power leaving each bus via the lines
# minus power injection at each bus shall equal to 0
Pmis = np.real(S) - P
Qmis = np.imag(S) - Q
return np.concatenate([Pmis[non_slack_loc],
Qmis[non_gen_loc]])
initial_guess = np.concatenate([ss.Bus.a0.v[non_slack_loc],
ss.Bus.v0.v[non_gen_loc]])
sol = fsolve(pf_large, initial_guess, args=(ss, P, Q, V, theta))
/var/folders/__/n5kx_m_s0tbg6n5qd7rh51700000gn/T/ipykernel_49347/1904602603.py:4: RuntimeWarning: The iteration is not making good progress, as measured by the
improvement from the last ten iterations.
sol = fsolve(pf_large, initial_guess, args=(ss, P, Q, V, theta))
V
array([1.06 , 1.045 , 1.01 , 1.01919406, 1.02024525,
1.07 , 1.06202177, 1.09 , 1.05594395, 1.05095431,
1.05697622, 1.05499798, 1.04999583, 1.0359623 ])
theta
array([ 0. , -0.08690277, -0.22195922, -0.18031046, -0.15328149,
-0.2482051 , -0.23335178, -0.23317861, -0.2607525 , -0.26354768,
-0.25814474, -0.26303971, -0.264607 , -0.27995707])
from andes.linsolvers.scipy import spmatrix_to_csc
Y_spmatrix = ss.Line.build_y()
Ycsc = spmatrix_to_csc(Y_spmatrix)
print(Ycsc)
<Compressed Sparse Column sparse matrix of dtype 'complex128'
with 54 stored elements and shape (14, 14)>
Coords Values
(0, 0) (6.025029055768224-19.498070205514384j)
(1, 0) (-4.999131600798035+15.263086523179553j)
(4, 0) (-1.025897454970189+4.234983682334831j)
(0, 1) (-4.999131600798035+15.263086523179553j)
(1, 1) (9.521323610814779-30.354715398779067j)
(2, 1) (-1.1350191923073958+4.781863151757718j)
(3, 1) (-1.686033150614943+5.115838325872083j)
(4, 1) (-1.7011396670944048+5.193927397969713j)
(1, 2) (-1.1350191923073958+4.781863151757718j)
(2, 2) (3.1209949022329564-9.850680129351638j)
(3, 2) (-1.9859757099255606+5.0688169775939205j)
(1, 3) (-1.686033150614943+5.115838325872083j)
(2, 3) (-1.9859757099255606+5.0688169775939205j)
(3, 3) (10.512989522036175-38.343131738471556j)
(4, 3) (-6.840980661495672+21.578553981691588j)
(6, 3) 4.889512660317341j
(8, 3) 1.8554995578159006j
(0, 4) (-1.025897454970189+4.234983682334831j)
(1, 4) (-1.7011396670944048+5.193927397969713j)
(3, 4) (-6.840980661495672+21.578553981691588j)
(4, 4) (9.568017783560265-34.97540411445229j)
(5, 4) 4.257445335253384j
(4, 5) 4.257445335253384j
(5, 5) (6.579923407466222-17.34073280991911j)
(10, 5) (-1.9550285631772604+4.0940743442404415j)
: :
(7, 6) 5.676979846721544j
(8, 6) 9.09008271975275j
(6, 7) 5.676979846721544j
(7, 7) -5.676979846721544j
(3, 8) 1.8554995578159006j
(6, 8) 9.09008271975275j
(8, 8) (5.326055039467359-24.28250637526788j)
(9, 8) (-3.902049552447428+10.365394127060915j)
(13, 8) (-1.4240054870199312+3.0290504569306034j)
(8, 9) (-3.902049552447428+10.365394127060915j)
(9, 9) (5.782934306147828-14.768337876521436j)
(10, 9) (-1.8808847537003996+4.402943749460521j)
(5, 10) (-1.9550285631772604+4.0940743442404415j)
(9, 10) (-1.8808847537003996+4.402943749460521j)
(10, 10) (3.83591331687766-8.497018093700962j)
(5, 11) (-1.525967440450974+3.1759639650294003j)
(11, 11) (4.014992027272893-5.427938591201612j)
(12, 11) (-2.4890245868219187+2.251974626172212j)
(5, 12) (-3.0989274038379877+6.102755448193116j)
(11, 12) (-2.4890245868219187+2.251974626172212j)
(12, 12) (6.724946148466233-10.66969354947068j)
(13, 12) (-1.1369941578063267+2.314963475105352j)
(8, 13) (-1.4240054870199312+3.0290504569306034j)
(12, 13) (-1.1369941578063267+2.314963475105352j)
(13, 13) (2.560999644826258-5.344013932035955j)
References#
Mariesa L Crow. Computational methods for electric power systems. CRC Press, 3 edition, 2016.