1
The Power Flow Equations
1.0 The Admittance Matrix
Current injections at a bus are analogous to power injections. The
student may have already been introduced to them in the form of
current sources at a node. Current injections may be either positive
(into the bus) or negative (out of the bus). Unlike current flowing
through a branch (and thus is a branch quantity), a current injection
is a nodal quantity. The admittance matrix, a fundamental network
analysis tool that we shall use heavily, relates current injections at
a bus to the bus voltages. Thus, the admittance matrix relates nodal
quantities. We motivate these ideas by introducing a simple
example. We assume that all electrical variables in this document
are given in the per-unit system.
Fig. 1 shows a network represented in a hybrid fashion using one-
line diagram representation for the nodes (buses 1-4) and circuit
representation for the branches connecting the nodes and the
branches to ground. The branches connecting the nodes represent
lines. The branches to ground represent any shunt elements at the
buses, including the charging capacitance at either end of the line.
All branches are denoted with their admittance values y
ij
for a
branch connecting bus i to bus j and y
i
for a shunt element at bus i.
The current injections at each bus i are denoted by I
i
.
y
4
y
3
y
1
y
2
I
4
I
3
I
2
I
1
4
3
2
1
y
34
y
23
y
13
y
12
Fig. 1: Network for Motivating Admittance Matrix
2
Kirchoff’s Current Law (KCL) requires that each of the current
injections be equal to the sum of the currents flowing out of the
bus and into the lines connecting the bus to other buses, or to the
ground. Therefore, recalling Ohm’s Law, I=V/z=Vy, the current
injected into bus 1 may be written as:
I
1
=(V
1
-V
2
)y
12
+ (V
1
-V
3
)y
13
+ V
1
y
1
(1)
To be complete, we may also consider that bus 1 is “connected” to
bus 4 through an infinite impedance, which implies that the
corresponding admittance y
14
is zero. The advantage to doing this
is that it allows us to consider that bus 1 could be connected to any
bus in the network. Then, we have:
I
1
=(V
1
-V
2
)y
12
+ (V
1
-V
3
)y
13
+ (V
1
-V
4
)y
14
+ V
1
y
1
(2)
Note that the current contribution of the term containing y
14
is zero
since y
14
is zero. Rearranging eq. 2, we have:
I
1
= V
1
(
y
1
+ y
12
+ y
13
+ y
14
) + V
2
(-y
12
)+ V
3
(-y
13
) + V
4
(-y
14
) (3)
Similarly, we may develop the current injections at buses 2, 3, and
4 as:
I
2
= V
1
(-y
21
) + V
2
(
y
2
+ y
21
+ y
23
+ y
24
) + V
3
(-y
23
) + V
4
(-y
24
) (4)
I
3
= V
1
(-y
31
)+ V
2
(-y
32
) + V
3
(
y
3
+ y
31
+ y
32
+ y
34
) + V
4
(-y
34
)
I
4
= V
1
(-y
41
)+ V
2
(-y
42
) + V
3
(-y
34
)+ V
4
(
y
4
+ y
41
+ y
42
+ y
43
)
where we recognize that the admittance of the circuit from bus k to
bus i is the same as the admittance from bus i to bus k, i.e., y
ki
=y
ik
From eqs. (3) and (4), we see that the current injections are linear
functions of the nodal voltages. Therefore, we may write these
equations in a more compact form using matrices according to:
3
4
3
2
1
4342414434241
3434323133231
2423242321221
1413121413121
4
3
2
1
V
V
V
V
yyyyyyy
yyyyyyy
yyyyyyy
yyyyyyy
I
I
I
I
(5)
The matrix containing the network admittances in eq. (5) is the
admittance matrix, also known as the Y-bus, and denoted as:
4342414434241
3434323133231
2423242321221
1413121413121
yyyyyyy
yyyyyyy
yyyyyyy
yyyyyyy
Y
(6)
Denoting the element in row i, column j, as Y
ij
, we rewrite eq. (6)
as:
44434241
34333231
24232221
14131211
YYYY
YYYY
YYYY
YYYY
Y
(7)
where the terms Y
ij
are not admittances but rather elements of the
admittance matrix. Therefore, eq. (5) becomes:
4
3
2
1
44434241
34333231
24232221
14131211
4
3
2
1
V
V
V
V
YYYY
YYYY
YYYY
YYYY
I
I
I
I
(8)
By using eq. (7) and (8), and defining the vectors V and I, we may
write eq. (8) in compact form according to:
4
4
3
2
1
4
3
2
1
,
I
I
I
I
I
V
V
V
V
V
VYI
(9)
We make several observations about the admittance matrix given
in eqs. (6) and (7). These observations hold true for any linear
network of any size.
1. The matrix is symmetric, i.e., Y
ij
=Y
ji
.
2. A diagonal element Y
ii
is obtained as the sum of admittances for
all branches connected to bus i, including the shunt branch, i.e.,
N
ikk
ikiii
yyY
,1
, where we emphasize once again that y
ik
is
non-zero only when there exists a physical connection between
buses i and k.
3. The off-diagonal elements are the negative of the admittances
connecting buses i and j, i.e., Y
ij
=-y
ji
.
These observations enable us to formulate the admittance matrix
very quickly from the network based on visual inspection. The
following example will clarify.
5
Example 1
Consider the network given in Fig. 2, where the numbers indicate
admittances.
j0.4
j0.3
j0.1
j0.2
I
4
I
3
I
2
I
1
4
3
2
1
2-j3
2-j5
1-j4
2-j4
Fig. 2: Circuit for Example 1
The admittance matrix is given by inspection as:
6.223200
327.1155241
0528.8442
041429.73
44434241
34333231
24232221
14131211
jj
jjjj
jjj
jjj
YYYY
YYYY
YYYY
YYYY
Y
2.0 The power flow equations
Define the net complex power injection into a bus as S
k
=S
gk
-S
dk
. In
this section, we desire to derive an expression for this quantity in
terms of network voltages and admittances. We begin by
reminding the reader that all quantities are assumed to be in per
unit, so we may utilize single-phase power relations. Drawing on
the familiar relation for complex power, we may express S
k
as:
6
S
k
=V
k
I
k
*
(10)
From eq. (8), we see that the current injection into any bus k may
be expressed as
j
N
j
kjk
VYI
1
(11)
where, again, we emphasize that the Y
kj
terms are admittance
matrix elements and not admittances. Substitution of eq. (11) into
eq. (10) yields:
*
1
*
*
1
j
N
j
kjkj
N
j
kjkk
VYVVYVS
(12)
Recall that V
k
is a phasor, having magnitude and angle, so that
V
k
=|V
k
|
k
. Also, Y
kj
, being a function of admittances, is therefore
generally complex, and we define G
kj
and B
kj
as the real and
imaginary parts of the admittance matrix element Y
kj
, respectively,
so that Y
kj
=G
kj
+jB
kj
. Then we may rewrite eq. (12) as
N
j
kjkjjkjk
N
j
kjkjjjkk
jj
N
j
kjkjkk
jj
N
j
kjkjkkj
N
j
kjkk
jBGVV
jBGVV
VjBGV
VjBGVVYVS
1
1
1
*
1
*
*
1
*
)()(
)(
)(
)(
(13)
7
Recall, from the Euler relation, that a phasor may be expressed as
complex function of sinusoids, i.e., V=|V|=|V|{cos+jsin}.
With this, we may rewrite eq. (13) as
N
j
kjkjjkjkjk
N
j
kjkjjkjkk
jBGjVV
jBGVVS
1
1
)()sin()cos(
)()(
(14)
If we now perform the algebraic multiplication of the two terms
inside the parentheses of eq. (14), and then collect real and
imaginary parts, and recall that S
k
=P
k
+jQ
k
, we can express eq. (14)
as two equations, one for the real part, P
k
, and one for the
imaginary part, Q
k
, according to:
N
j
jkkjjkkjjkk
N
j
jkkjjkkjjkk
BGVVQ
BGVVP
1
1
)cos()sin(
)sin()cos(
(15)
The two equations of (15) are called the power flow equations, and
they form the fundamental building block from which we attack
the power flow problem.
3.0 Solving the power flow problem
The standard power flow problem is as follows:
Given that for each bus (node) in the network, we know 2 out of
the following 4 variables: P
k
, Q
k
, |V
k
|, θ
k
, so that for each bus, there
are two equations available – those of eq. (15) above, and there are
two unknown variables. Thus the power flow problem is to solve
8
the power flow equations of (15) for the remaining two variables
per bus.
This problem is one where we are required to solve simultaneous
nonlinear equations. Because most power systems are very large
interconnections, with many buses, the number of power flow
equations (and thus the number of unknowns) is very large. For
example, a model of the eastern interconnection in the US can have
50,000 buses.
The approach to solving the power flow problem is to use an
iterative algorithm. The Newton-Raphson algorithm is the most
commonly used algorithm in commercial power flow programs.
Starting with a reasonable guess at the solution (where the
“solution” is a numerical value of all of the unknown variables),
this algorithm checks to see how close the solution is, and then if it
is not close enough, updates the solution in a direction that is sure
to improve it, and then repeats the check. This process continues
until the check is satisfied. Usually, this process requires 5-20
iterations to converge to a satisfactory solution. For large
networks, it is computationally intensive.
In this class, we are very interested in optimization methods for
finding maximum surplus solutions to the problem of how to
dispatch the generation. So far, we have dealt with problems where
all generation and load was considered to be at the same bus (node)
and were thus able to ignore the network. But in reality, generation
and load are located at various buses, and the transportation
mechanism for moving electrical energy from supply to
consumption is the transmission network. If there are losses or
constraints in the transmission network (which there are), these
will influence how supply can be allocated, and the most
economically desirable solutions may not be feasible.
9
To account for the network in the economic optimization problems
we have posed, we must account for the equations that correspond
to the network. These are the power flow equations. This can be
done, and problem that results is referred to as the optimal power
flow (OPF). However, because the power flow equations are
nonlinear, the OPF requires a nonlinear optimization method for its
solution. Nonlinear optimization (usually called nonlinear
programming instead of linear programming) is a rich, interesting,
and highly applicable area. You can take entire courses on this
subject (see, for example, IE 631).
But we do not have time in this class to learn nonlinear
programming methods. And fortunately, since we have learned
linear programming, we do not need to do so if we can convert our
nonlinear problem into a linear one.
In our problem, where we desire to maximize social surplus, the
objective function may be nonlinear. It is possible to convert a
nonlinear objective function into a linear objective function using
piecewise linear approximations. This method is very effective in
approximating an objective function that is a separable function
(can be separated into components), since each component is a
function of only one variable. In other words, we are able to apply
piecewise linear approximation to each individual utility or cost
function because each one is a function of only one variable. Fig. 3
illustrates a piecewise linear approximation of a cost function.
P
i,min
P
i3
P
i2
P
i1
C
i
$/hr
P
i
(MW)
Fig. 3
10
The power flow equations are functions of many variables, and it is
very complex to see how to apply piecewise linear approximations,
since a piecewise linear approximation for any one variable will
depend on the value of the other variables.
Thus, we seek another method of converting our nonlinear power
flow equations into linear equations.
4.0 Approximations to the power flow equations
Let’s reconsider the power flow equations:
N
j
jkkjjkkjjkk
N
j
jkkjjkkjjkk
BGVVQ
BGVVP
1
1
)cos()sin(
)sin()cos(
(15)
We will make use of three practical observations regarding high
voltage electric transmission systems.
Observation 1: The resistance of transmission circuits is
significantly less than the reactance. Usually, it is the case that the
x/r ratio is between 2 and 10. So any given transmission circuit
with impedance of z=r-jx will have an admittance of
jbg
xr
jx
xr
r
xr
jxr
jxr
jxr
jxrjxrz
y
2222
22
111
(16)
From eq. (16), we see that
22
xr
r
g
and
22
xr
x
b
(17)
If r is very small compared to x, then we observe that g will be
very small compared to b, and it is reasonable to approximate eqs.
(17) as
11
0g
and
x
b
1
(18)
Now, if g=0, then the real part of all of the Y-bus elements will
also be zero, that is, g=0G=0.
Applying this conclusion to the power flow equations of eq. (15):
N
j
jkkjjkk
N
j
jkkjjkk
BVVQ
BVVP
1
1
)cos(
)sin(
(19)
Observation 2: For most typical operating conditions, the
difference in angles of the voltage phasors at two buses k and j
connected by a circuit, which is θ
k
-θ
j
for buses k and j, is less than
10-15 degrees. It is extremely rare to ever see such angular
separation exceed 30 degrees. Thus, we say that the angular
separation across any transmission circuit is “small.”
Consider that, in eqs. (19), the angular separation across a
transmission circuit, θ
k
-θ
j
, appears as the argument of the
trigonometric functions sine and cosine. What do these functions
look like for small angles? We can answer this question by
recalling that the sine and cosine functions represent the vertical
and horizontal components of a unit (length=1) vector making an
angle δ=θ
k
-θ
j
with the positive x-axis, as illustrated in Fig. 3.
12
cosδ
δ
sinδ
Fig. 3: Trig functions of a small angle
In Fig. 3, it is clear that as the angle δ=θ
k
-θ
j
gets smaller and
smaller, the cosine function approaches 1.0.
One might be tempted to accept the approximation that the sine
function goes to zero. This it does, as the angle goes to zero. But
an even better approximation is that the sine of a small angle is the
angle itself (when the angle is given in radians). This can be
observed in Fig. 3 from the fact that the vertical line, representing
the sine, is almost the same length as the indicated radial distance
along the circle, which is the angle (when measured in radians).
Applying these conclusions from observation 2 to eqs. (19):
N
j
kjjkk
N
j
jkkjjkk
BVVQ
BVVP
1
1
)(
(20)
Note that we have made significant progress at this point, in
relation to obtaining linear power flow equations, since we have
13
eliminated the trigonometric terms. However, we still have product
terms in the voltage variables, and so we are not done yet. Our next
and last observation will take care of these product terms.
Before we do that, however, let’s investigate the expressions of eq.
(20) a little.
Recall that the quantity B
kj
is not actually a susceptance but rather
an element in the Y-bus matrix.
If k≠j, then B
kj
=-b
kj
, i.e., the Y-bus element in row k column j is
the negative of the susceptance of the circuit connecting bus k to
bus j.
If k=j, then
N
kjj
kjkkk
bbB
,1
Reactive power flow:
The reactive power flow equation of eqs. (20) may be rewritten by
pulling out the k=j term from the summation.
N
kj
j
kjjkkkk
N
j
kjjkk
BVVBVBVVQ
,1
2
1
Then substitute susceptances for the Y-bus elements:
N
kj
j
kjjk
N
kjj
kjkkk
N
kj
j
kjjk
N
kjj
kjkkk
N
kj
j
kjjk
N
kjj
kjkkk
bVVbVbV
bVVbVbV
bVVbbVQ
,1,1
22
,1,1
22
,1,1
2
14
Now bring all the terms in the two summations under a single
summation.
N
kjj
kjjkkjkkkk
bVVbVbVQ
,1
22
Factor out the |V
k
| and the –b
kj
in the summation:
N
kjj
jkkjkkkk
VVbVbVQ
,1
2
(21)
Note because we defined the circuit admittance between buses k
and j as y
kj
=g
kj
+jb
kj
, and because all circuits have series elements
that are inductive, the numerical value of b
kj
is negative. Thus, we
can rewrite eq. (21) as
N
kjj
jkkjkkkk
VVbVbVQ
,1
2
(22)
So there are two main terms in eq. (22).
The first term corresponds to the reactive power supplied (if a
capacitor) or consumed (if an inductor) by the shunt susceptance
modeled at bus k.
The second term corresponds to the reactive power flowing on
the circuits connected to bus k. Only these circuits will have
nonzero b
kj
. One sees, then, that each circuit will have per-unit
reactive flow in proportion to (a) the bus k voltage magnitude
and (b) the difference in per-unit voltages at the circuit’s
terminating buses. The direction of flow will be from the higher
voltage bus to the lower voltage bus.
Real power flow: Now consider the real power flow equation from
eqs. (20), and, as with the reactive power flow equation, let’s pull
out the j=k term. Thus,
N
kj
j
jkkjjkkkkkk
N
j
jkkjjkk
BVVBVBVVP
,1
2
1
)()()(
Here, we see that the first term is zero, so that:
15
N
kj
j
jkkjjkk
BVVP
,1
)(
(23)
Some comments about this expression:
There is no “first term” corresponding to shunt elements as
there was for the reactive power equation. The reason for this is
that, because we assumed that R=0 for the entire network, there
can be no shunt resistive element in our model. This actually
conforms to reality since we never connect a resistive shunt in
the transmission system (this would be equivalent to a giant
heater!). The only place where we do actually see an effect
which should be modeled as a resistor to ground is in
transformers the core loss is so modeled. But the value of this
resistance tends to be extremely large, implying the
corresponding conductance (G) is extremely small, and it is
very reasonable to assume it is zero.
Therefore the term that we see in eq. (23) represents the real
power flow on the circuits connected to bus k. One sees, then,
that each circuit will have per-unit real power flow in
proportion to (a) the bus k and j voltage magnitudes and (b) the
angular difference across the circuit. Furthermore, recalling that
B
kj
=-b
kj
, and also that all transmission circuits have series
elements that are inductive, the numerical value of b
kj
is
negative, implying that the numerical value of B
kj
is positive.
Therefore, the direction of flow will be from the bus with the
larger angle to the bus with the smaller angle.
Observation 3: In the per-unit system, the numerical values of
voltage magnitudes |V
k
| and |V
j
| are very close to 1.0. Typical
range under most operating conditions is 0.95 to 1.05. Let’s
consider the implications of this fact in terms of the expressions for
reactive and real power flow eqs. (22) and (23), repeated here for
convenience:
16
N
kjj
jkkjkkkk
VVbVbVQ
,1
2
N
kj
j
jkkjjkk
BVVP
,1
)(
Given that 0.95<|V
k
| and |V
j
|<1.05, then we incur little error in the
above expressions if we assume |V
k
|=|V
j
|=1.0 everywhere that they
occur as a multiplying factor. We cannot make this approximation,
however, where they occur as a difference, in the reactive power
equation, because the difference of two numbers close to 1.0 can
range significantly. For example, 1.05-0.95=0.1, but 1.01-
1.0=0.01, an order of magnitude difference.
Making this approximation results in:
N
kjj
jkkjkk
VVbbQ
,1
(24)
N
kj
j
jkkjk
BP
,1
)(
(25)
With these equations, we can narrow our statements about power
flow.
Reactive power flow across circuits is determined by the
difference in the voltage phasor magnitudes between the
terminating buses.
Real power flow across circuits is determined by the difference
in voltage phasor angles between the terminating buses.
Finally, it is interesting to note that the disparity between the
maximum reactive power flow and the maximum real power flow
across a circuit.
The reactive power flow equation is proportional to the circuit
susceptance and the difference in voltage phasor magnitudes.
17
The maximum difference in voltage phasor magnitudes will be
on the order of 1.05-0.95=0.1.
The real power flow equation is proportional to the circuit
susceptance and the difference in voltage phasor angles. The
maximum difference in voltage phasor angles will be, in
radians, about 0.52 (which corresponds to 30 degrees).
We see from these last two bullets that real power flow across
circuits tends to be significantly larger than reactive power flow,
i.e., usually, we will see that
kjkj
QP
This conclusion is consistent with operational experience, which is
characterized by an old operator’s saying: “Vars don’t travel.”
5.0 Real vs. Reactive Power Flow
Recall that our original intent was to represent the network in our
optimization problem because of our concern that network
constraints may limit the ability to most economically dispatch the
generation. There are actually several different causes of network
constraints, but here, we will limit our interest to the type of
constraint that is most common in most networks, and that is
circuit overload.
Circuit overload is caused by high current magnitude. When the
current magnitude exceeds a given threshold for a particular circuit
(called the rating), we say that overload has occurred.
In the per-unit system, we recall that
kjkkjkjkj
IVjQPS
where V
k
is the bus k nodal voltage phasor and I
kj
is the phasor of
the current flowing from bus k to bus j. Thus, we have that:
18
k
kjkj
kj
V
jQP
I
Taking the magnitude (since that is what determines circuit
overload), we have:
k
kjkj
kj
V
QP
I
22
Given our conclusion on the previous page that generally,
P
kj
>>Q
kj
, we may approximate the above expression according to:
k
kj
k
kj
kj
V
P
V
P
I
2
and if |V
k
|≈1.0, then we have that
kjkj
PI
Thus, in assessing circuit overload, it is reasonable to look at real
power flows only. As a result of this conclusion, we will build into
our optimization formulation only the real power flow equations,
i.e., eq. (25).
6.0 The DC Power Flow – an example
Let’s study the real power flow expression given in eq. (25).
N
kj
j
jkkjk
BP
,1
)(
19
It is worthwhile to perform a simple example to illustrate use of
this expression.
Consider the 4-bus network given in Fig. 4. All 5 lines have the
same admittance, and this admittance has no real part indicating
we are assuming R=0 for this network. The real power values for
each of the three generators and each of the two loads are given.
All numerical quantities are given in per-unit.
The problem is to compute the real power flows on all circuits.
y
13
=-j10
y
14
=-j10
y
34
=-j10
y
23
=-j10
y
12
=-j10
P
g1
=2pu
P
d3
=4pu
P
d2
=1pu
1
2
3
4
P
g2
=2pu
P
g4
=1pu
Fig. 4: Four-bus network used in example
We first write down eq. (25) for each bus, beginning with bus 1.
414114313113212112
4114311321121
)()()(
BBBBBB
BBBP
Collecting terms in the same variables results in:
41431321211413121
BBBBBBP
(26)
Repeating the process for bus 2:
20
424224323223121221
4224322312212
)()()(
BBBBBB
BBBP
Again, collecting terms in the same variables results in:
42432322423211212
BBBBBBP
(27)
Repeating eqs. (26) and (27), together with the relations for buses
3 and 4, yields:
41431321211413121
BBBBBBP
42432322423211212
BBBBBBP
43433432312321313
BBBBBBP
44342413432421414
BBBBBBP
We can write these equations in matrix form, according to:
4
3
2
1
434241434241
343432313231
242324232121
141312141312
4
3
2
1
BBBBBB
BBBBBB
BBBBBB
BBBBBB
P
P
P
P
(28)
Remember, the left-hand-side vector is the injections, which is the
generation less the demand.
To get the matrix, it is helpful to first write down the Y-bus:
44434241
34333231
24232221
14131211
BBBB
BBBB
BBBB
BBBB
jY
4342414434241
3434323133231
2423242321221
1413121413121
bbbbbbb
bbbbbbb
bbbbbbb
bbbbbbb
j
21
2010010
10301010
0102010
10101030
jY
So we readily observe here that, for example, B
11
=-30, B
12
=10,
B
13
=10, and B
14
=10, and it is similar for the other three rows.
So using the Y-bus values, we can express eq. (28) as:
4
3
2
1
2010010
10301010
0102010
10101030
1
4
1
2
(29)
(Observe that we omit the j). Then, the angles are given by:
1
4
1
2
2010010
10301010
0102010
10101030
1
4
3
2
1
(30)
However, when we evaluate the above expression by taking the
inverse of the indicated matrix, we find it is singular, i.e., it does
not have an inverse. The problem here is that we have a
dependency in the 4 equations, implying that one of the equations
may be obtained from the other three. For example, if we add the
bottom three rows and then multiply by -1, we get the top row (in
terms of the injection vector, this is just saying that the sum of the
generation must equal the demand).
This dependency occurs because all four angles are not
independent. What is important are the angular differences. Thus,
we are free to choose any one of them as the reference, with a
fixed value of 0 degrees. This angle is then no longer a variable (it
22
is 0 degrees), and, referring to eq. (29), the corresponding column
in the matrix may be eliminated, since those are the numbers that
get multiplied by this 0 degree angle.
To fix this problem, we need to eliminate one of the equations and
one of the variables (by setting the variable to zero).
We choose to eliminate the first equation, and set the first variable,
θ
1
, to zero (which means we are choosing θ
1
as the reference).
This results in:
025.0
15.0
025.0
1
4
1
20100
103010
01020
1
4
3
2
(31)
The solution on the right-hand-side gives the angles on the bus
voltage phasors at buses, 2, 3, and 4.
However, the problem statement requires us to compute the power
flows on the lines (this is usually the information needed by
operational and planning engineers as they study the power
system).
We can get the power flows easily by employing just one term
from the summation in eq. (25), which gives the flow across circuit
k-j according to:
)(
jkkjkj
BP
(32)
We utilize the Y-bus elements together with the bus angles given
by eq. (31) to make these calculations, as follows:
25.0)025.00(10)(
211212
BP
5.1)15.00(10)(
311313
BP
23
25.0)025.00(10)(
411414
BP
25.1)15.0025.0(10)(
322323
BP
25.1)025.015.0(10)(
433434
BP
These computed flows are illustrated in Fig. 5. The power flowing
into a bus equals the power flowing out of that bus.
P
13
=1.5
P
14
=0.25
P
43
=1.25
P
23
=1.25
P
12
=0.25
P
g1
=2pu
P
d3
=4pu
P
d2
=1pu
1
2
3
4
P
g2
=2pu
P
g4
=1pu
Fig. 5: Computed flows for four-bus network used in example
7.0 The DC Power Flow – Generalization
We desire to generalize the above procedure.
We assume that we are given the network with the following
information:
Total number of buses is N, total number of branches is M.
Bus number 1 identified as the reference
Real power injections at all buses except bus 1
Network topology
Admittances for all branches.
24
The DC power flow equations, based on eq. (25) are given in
matrix form as
'BP
(33)
where
P is the vector of nodal injections for buses 2, …, N
θ is the vector of nodal phase angles for buses 2,…N
B’ is the “B-prime” matrix. Generalization of its development
requires a few comments.
Development of the B’ matrix:
Compare the matrix of eq. (28) with the Y-bus matrix, all repeated
here for convenience:
4
3
2
1
434241434241
343432313231
242324232121
141312141312
4
3
2
1
BBBBBB
BBBBBB
BBBBBB
BBBBBB
P
P
P
P
44434241
34333231
24232221
14131211
BBBB
BBBB
BBBB
BBBB
jY
4342414434241
3434323133231
2423242321221
1413121413121
bbbbbbb
bbbbbbb
bbbbbbb
bbbbbbb
j
From the above, we can develop a procedure to obtain the B’
matrix from the Y-bus, as follows:
1. Remove the “j” from the Y-bus.
2. Replace diagonal element B’
kk
with the sum of the non-diagonal
elements in row k. Alternatively, subtract b
k
(shunt term) from
B
kk
, & multiply by -1 (if there is no b
k
, then just multiple by -1).
3. Multiply all off-diagonals by -1.
4. Remove row 1 and column 1.
If there are
no b
k
, then
steps 2-3
simplify to
“multiply
Y-bus by -1
25
Comparison of the numerical values of the Y-bus with the
numerical values of the B’ matrix for our example will confirm the
above procedure:
2010010
10301010
0102010
10101030
jY
20100
103010
01020
'B
Another way to remember the B’ matrix is to observe that, since its
non-diagonal elements are the negative of the Y-bus matrix, the B’
non-diagonal elements are susceptances. However, one must be
careful to note that the B’ matrix element in position row k,
column j is the susceptance of the branch connecting buses k+1
and j+1, since the B’ matrix does not have a column or row
corresponding to bus 1.
Question: Why are shunt terms excluded in the B’ matrix? That is,
why does excluding them not affect real power calculations?
Although eq. (33) provides the ability to compute the angles, it
does not provide the line flows. A systematic method of computing
the line flows is:
)( ADP
B
(34)
where:
P
B
is the vector of branch flows. It has dimension of M x 1.
Branches are ordered arbitrarily, but whatever order is chosen
must also be used in D and A.
θ is (as before) the vector of nodal phase angles for buses 2,…N
26
D is an M x M matrix having non-diagonal elements of zeros;
the diagonal element in position row k, column k contains the
negative of the susceptance of the k
th
branch.
A is the M x N-1 node-arc incidence matrix. It is also called the
adjacency matrix, or the connection matrix. Its development
requires a few comments.
Development of the node-arc incidence matrix:
This matrix is well known in any discipline that has reason to
structure its problems using a network of nodes and “arcs” (or
branches or edges). Any type of transportation engineering is
typical of such a discipline.
The node-arc incidence matrix contains a number of rows equal to
the number of arcs and a number of columns equal to the number
of nodes.
Element (k,j) of A is 1 if the k
th
branch begins at node j, -1 if the
k
th
branch terminates at node j, and 0 otherwise.
A branch is said to “begin” at node j if the power flowing across
branch k is defined positive for a direction from node j to the other
node.
A branch is said to “terminate” at node j if the power flowing
across branch k is defined positive for a direction to node j from
the other node.
Note that matrix A is of dimension M x N-1, i.e., it has only N-1
columns. This is because we do not form a column with the
reference bus, in order to conform to the vector θ, which is of
dimension (N-1) x 1. This works because the angle being excluded,
θ
1
, is zero.
27
We can illustrate development of the node-arc incidence matrix for
our example system. Consider numbering the branches as given in
Fig. 6. Positive direction of flow is as given by the indicated
arrows. The numbers in the circles are bus (node) numbers. The
numbers next to each branch are branch numbers.
5
1
4
3
2
P
g1
=2pu
P
d3
=4pu
P
d2
=1pu
1
2
3
4
P
g2
=2pu
P
g4
=1pu
Fig. 6: Branches numbering for development of incidence matrix
Therefore, the node-arc incidence matrix is given as
number branch
5
4
3
2
1
01-0
11-0
01-1
001-
1-00
A
4 3 2
number node
28
The D-matrix is formed by placing the negative of the susceptance
of each branch along the diagonal of an M x M matrix, where
M=5.
100000
010000
001000
000100
000010
D
Combining A, D, and θ based on eq. (34) yields
)( ADP
B
4
3
2
100000
010000
001000
000100
000010
01-0
11-0
01-1
001-
1-00
3
43
32
2
4
3
43
32
2
4
5
4
3
2
1
10
)(10
)(10
10
10
100000
010000
001000
000100
000010
B
B
B
B
B
P
P
P
P
P
Plugging in the solution for θ that we obtained, which was:
29
025.0
15.0
025.0
4
3
2
We find that the above evaluates to
5.1
25.1
25.1
25.0
25.0
5
4
3
2
1
B
B
B
B
B
P
P
P
P
P
This solution, obtained systematically, in a way that can be
efficiently programmed, agrees with the solution we obtained
manually and is displayed in Fig. 5, repeated here for convenience.
P
13
=1.5
P
14
=0.25
P
43
=1.25
P
23
=1.25
P
12
=0.25
P
g1
=2pu
P
d3
=4pu
P
d2
=1pu
1
2
3
4
P
g2
=2pu
P
g4
=1pu
