A Power Flow in Matlab

Simulation and Analysis of Power Systems and Power Electronics in MATLAB
Montreal, September 17, 2001

A Power Flow in Matlab

This section describes the development and extension of a basic power flow program in Matlab. To the extent feasible, complex vector operations are used, even for the calculation of the Jacobian matrix[1].

The problem

The main problem of interest in power system analysis is the solution of the so-called power flow (or load flow) problem. The assumptions for the basic power flow problem are:

· A passive linear network of admittances describes the transmission grid. The nodal admittance matrix Y-bus (denoted by Y) characterizes the behavior of the nodal voltages to nodal current injections. The components of Y are complex numbers of the form . There are n nodes in the network. The dimension of Y is n by n.

· Every node can have active and reactive generation and active and reactive demand, each represented as a complex number. The complete set of generation injections for the network is denoted by the vector Sg, and the demands by the vector Sd, each a vector of dimension n.

· Every node (or bus) is classified as either a PV or a PQ node. For a PV node, the voltage magnitude and the net power injection are fixed, while the generation Q is a variable. For a PQ node, both the active and reactive net power injections are fixed (either positive or negative).

· One node (usually a PV node) is designated as the angle reference for all voltages. At least one node (usually the same PV node) is designated as a node where the net injection P is not specified but rather determined as a variable. Thus, in this node both voltage magnitude and angle are known. The node is called the slack node. When multiple nodes are designated as slack nodes, some means for resolving the uncertainty as to which node picks up the “slack” is necessary.

The following complex equations must be satisfied at every node in the network, no exceptions:

where denotes element-wise product of two vectors (using Matlab notation), and * as a superscript denotes the complex conjugate. The dot denotes the ordinary dot product of a matrix times a vector.

Since at least one voltage V is known (and possibly many magnitudes of several other voltages), this means that only a subset of these equations is solved as a simultaneous set. In particular, the equations that are selected for a simultaneous solution are:

· The real and reactive parts of these equations for all PQ nodes.

· The real-only part of the above equations for all PV nodes excluding the slack.

Likewise, the only variables of immediate interest are:

· The voltage magnitudes and angles (or real and imaginary parts of the voltages) for all PQ nodes.

· The voltage angles at all PV nodes.

All other variables (such as Q at PV nodes or P and Q at the slack node) can be determined after the fact from the above equation.

There are situations where the explicit rectangular form of the power flow equations is more desirable than the concise complex vectorized version above. For these situations, the following is the set of power flow equations of interest:

where gij and boj are the real and imaginary components of the entries of the nodal admittance matrix, di is the angle of Vi, ai denotes all neighbors or node i excluding i itself, and denotes the neighbors of i including i itself.

In matrix form, these equations can be expressed as:

and from here the Jacobian can be determined from:

Building Y-bus (vectorized version)

The following vectorized code builds the Y-bus matrix:

function Ybus=bldybus(S)
% bldybus.m Build the y-bus matrices
Ybus=sparse(S.Bus.n,S.Bus.n);
YL=[];
if isfield(S.Bus,'G'),
YL=S.Bus.G;
end
if isfield(S.Bus,'B'),
YL=YL+i*S.Bus.B;
end
if ~isempty(YL),
Ybus=Ybus+diag(sparse(YL));
end
n=S.Bus.n;
Y11=(1./S.Branch.Z+(i*S.Branch.B)/2).*S.Branch.Status;
Ybus=Ybus+sparse(S.Branch.From,S.Branch.From,Y11,n,n);
Y12=(-1./(S.Branch.TAP.*S.Branch.Z)).*S.Branch.Status;
Ybus=Ybus+sparse(S.Branch.From,S.Branch.To,Y12,n,n);
Y21=(-1./(conj(S.Branch.TAP).*S.Branch.Z)).*S.Branch.Status;
Ybus=Ybus+sparse(S.Branch.To,S.Branch.From,Y21,n,n);
Y22=(1./((abs(S.Branch.TAP).^2).*S.Branch.Z)+(i*S.Branch.B)/2).*S.Branch.Status;
Ybus=Ybus+sparse(S.Branch.To,S.Branch.To,Y22,n,n);
if isfield(S.Branch,'YI'),
Ybus=Ybus+sparse(S.Branch.From,S.Branch.From,S.Branch.YI,n,n);
end
if isfield(S.Branch,'YJ'),
Ybus=Ybus+sparse(S.Branch.To,S.Branch.To,S.Branch.YJ,n,n);
end
if isfield(S,'Shunt'),
Ybus=Ybus+sparse(S.Shunt.I,S.Shunt.I,i*S.Shunt.BINIT,n,n);
end;

The main code

The code to solve the power flow problem is:

function S=mpflow(S)
% Save as mpflow.m file
if nargin<1, S=readcf; end
S.Ybus=bldybus(S); % Build the Y-bus matrix
S.snet= sparse(S.Machine.BusRef,1,S.Machine.MW+j*S.Machine.MVAR,S.Bus.n,1)...
-sparse(S.Load.BusRef,1,S.Load.MW+j*S.Load.MVAR,S.Bus.n,1);
S=pfsolve(S);
if S.mismatch>0.001,
disp('Power flow failed to converge');
else
disp('Power flow solved');
end;

This code is virtually self-explanatory:

· The data are read from the common format file using the readcf.m function. This code is listed in the appendix. The result of reading this information is stored in several tables: The Z table contains branch data, the L table contains load data and the A table contains generator data. In addition, the slack bus(es), generation bus(es) and load bus(es) are identified.

· The appropriate matrices are constructed using the bldybus.m function, also available as a pcode file.

· Some data conversion is performed (the solver assumes that net injections are used and does not distinguish between the positive power delivered to loads and the negative power delivered by generators).

· The solver is called. It returns the complex voltages at every bus, and also the mismatch or maximum error.

· The user is given an indication of completion.

The main routine within this code if pfsolve.m. While we have only provided the pcode for this routine, it is useful to examine a simplified version of this routine:

function S=pfsolve(S)
% This routine solves the power flow problem using a Newton-Raphson
% iteration with full Jacobian update at every iteration.
% Originally developed by Professor Christopher DeMarco, 1997
% and expanded by Fernando Alvarado 1998, 1999
%
% Compute initial mismatch
MAXITER=15; n=S.Bus.n;
vbus=S.Bus.Voltages;
S.mismatch = Inf; % pick off largest component of relevant mismatch
% Newton-Raphson Iteration
S.itcnt = 0;
while S.mismatch > .0001
fullmiss = pfmiss(S);
rmiss=[real(fullmiss(S.PVlist));real(fullmiss(S.PQlist)); ...
imag(fullmiss(S.PQlist))];
S.mismatch = max(abs(rmiss));
[dsdd, dsdv] = pflowjac(S.Ybus,vbus);
% dsdd and dsdv are composed of all (complex) partial derivatives
% rjac is a selection from these
rjac = [ ...
real(dsdd(S.PVlist,S.PVlist)) real(dsdd(S.PVlist,S.PQlist)) ...
real(dsdv(S.PVlist,S.PQlist)); ...
real(dsdd(S.PQlist,S.PVlist)) real(dsdd(S.PQlist,S.PQlist)) ...
real(dsdv(S.PQlist,S.PQlist)); ...
imag(dsdd(S.PQlist,S.PVlist)) imag(dsdd(S.PQlist,S.PQlist)) ...
imag(dsdv(S.PQlist,S.PQlist)) ];
x=[angle(vbus(S.PVlist))*pi/180;
angle(vbus(S.PQlist))*pi/180;
abs(vbus(S.PQlist))];
dx = -rjac\rmiss; % Here is the actual update
del=angle(vbus);
vmag=abs(vbus);
n1=length(S.PVlist); n2=n1+length(S.PQlist); n3=length(dx);
del(S.PVlist)=del(S.PVlist)+dx(1:n1);
del(S.PQlist)=del(S.PQlist)+dx((n1+1):n2);
vmag(S.PQlist)=vmag(S.PQlist)+dx((n2+1):n3);
vbus=vmag.*exp(j*del);
S.Bus.Voltages=vbus;
S.itcnt=S.itcnt+1;
if (S.itcnt > MAXITER)|((S.itcnt>3)&(S.mismatch>100)),
S.itcnt=Inf;
break;
end;
end
if (S.itcnt<=MAXITER),
S.Bus.Voltages=vbus;
disp([int2str(S.itcnt) ' iterations']);
end

The following comments apply to this code:

The implementation of the idea of using a single vectorized complex computation to obtain a complex Jacobian matrix from which the real Jacobian is extracted originated with Professor Christopher DeMarco of the University of Wisconsin in Madison.
The input voltages are used as the default value of the output voltages in case the routine terminates prematurely due to convergence or other difficulties.
The norm of the error (mismatch) is used to determine convergence tolerance is the infinite norm (largest error).
The tolerance for this error is fixed at 0.0001, a suitable value for most applications. Normal termination is based on the mismatch being smaller that this value.
Although most computations are done in complex mode, it is still necessary to have available vectors with the polar representation of the voltages.
This is where the complex version of the complete Jacobian is computed.
This is where the real portion of the Jacobian of interest to us is extracted.
This is the main computational step: the solution of the simultaneous sparse linear equations.
This is where the correction is added to the previous solution value.
The mismatch fullmiss is computed for all buses, but only the portion of interest to us rmiss is extracted.
The solution also can terminate by exceeding the maximum allowed number of iterations,.

Extending the power flow

The code above simply solves the basic power flow with no limits and finds the voltages and angles at every node. After this is done, there are generally a number of additional computations of interest. Some of these are described next.

Flows and overloads

The real objective of solving a power flow is quite often to find the individual line flows in the network, and most particularly, to find whether there are any overloaded lines under the presumed conditions. While developing the code for doing this computation using for and while loops is straightforward. This section illustrates how to write this code in a vectorized manner. The following code determines and displays every flow in every line in the network, identifies any overloads, and displays a list of overloaded lines.

function [Sflow,iOvld]=flows(S)
% Calculate all the line flows and report overloads
Vf=S.Bus.Voltages;
Sflow=(Vf(S.Branch.From,1)-Vf(S.Branch.To,1))./S.Branch.Z+ ...
Vf(S.Branch.From,1).*(j*S.Branch.B/2);
iOvld=find((S.Branch.Status.*abs(Sflow))>S.Branch.RateValue);

The following comments apply to this code:

The computations of all flows is done is one single vectorized line of code.
An even better way to do this is to use the concept of the flow admittance matrix Yflow (not illustrated).
The find function finds all overloads in a vectorized manner.

Penalty factors

The following code illustrates the construction of the loss penalty factors in accordance with the Transposed Jacobian method (Alvarado 1979):

Function beta=penalty(S)
% Construct the penalty factors
[dSdd,dSdv]=pflowjac(S.Ybus,S.Bus.Voltages);
Jac=[real(dSdd) real(dSdv); imag(dSdd) imag(dSdv)];
[n1,n2]=size(Jac);
Jacr=Jac(2:n1,2:n2); % Assuming the slack bus is bus 1
rhs=-Jac(1,2:n2)';
beta=[1; Jacr'\rhs];
disp('Beta');
disp(full(reshape(beta,n2/2,2)));

The following comments apply to this code:

The same Jacobian computation routine as in the power flow is used. In fact, it is possible to save this Jacobian re-computation and just use the last Jacobian from the solution itself.
As before, real and imaginary parts of this complex Jacobian are extracted. Unlike the previous case, however, all entries of the Jacobian are used. This assumes that no buses are holding voltage constant.
The slack bus row and column are then excluded. In this simple version, the slack bus must be node 1.
The right hand side corresponds to the row of the omitted slack bus active power.
The solution involves the use of the transposed reduced Jacobian.
The reshape function permits a more elegant display of the computed beta inverse penalty factors.

Using the power flow

We now illustrate the use of the above code to solve several simple power flow cases. Consider the case illustrated below, from Wood and Wollenberg:

The common format data for this case is (save it as 6bus.cf):

Wood and Wollenberg test case
BUS DATA FOLLOWS
1 BUS1 1 1 3 1.0500 0.00 0.00 0.00 100.00 .00 100.00
2 BUS2 1 1 2 1.0500 0.00 0.00 0.00 50.00 .00 100.00
3 BUS3 1 1 2 1.0700 0.00 0.00 0.00 60.00 .00 100.00
4 BUS4 1 1 0 1.0000 0.00 70.00 70.00 .00
5 BUS5 1 1 0 1.0000 0.00 70.00 70.00 .00
6 BUS6 1 1 0 1.0000 0.00 70.00 70.00 .00
-999
BRANCH DATA FOLLOWS
1 2 1 1 1 0 .100000 .200000 .020000 100
1 4 1 1 1 0 .050000 .200000 .020000 100
1 5 1 1 1 0 .080000 .300000 .030000 100
2 3 1 1 1 0 .050000 .250000 .030000 100
2 4 1 1 1 0 .050000 .100000 .010000 100
2 5 1 1 1 0 .100000 .300000 .020000 100
2 6 1 1 1 0 .070000 .200000 .025000 100
3 5 1 1 1 0 .120000 .260000 .025000 100
3 6 1 1 1 0 .020000 .100000 .010000 100
4 5 1 1 1 0 .200000 .400000 .040000 100
5 6 1 1 1 0 .100000 .300000 .030000 100
-999
INTERCHANGE DATA FOLLOWS
-9
TIE LINE DATA FOLLOWS
-999
LOSS ZONE DATA FOLLOWS
-99
END OF DATA

Because of width restrictions, the “desired voltage” field generally supplied with PV type buses has been excluded in the above data. However, the desired voltage magnitude for these buses is listed in the voltage magnitude field. Execution of the above power flow with this data as input yields the following:

» fname='6bus';
» mpflow;
mpflow (c)1997 F. Alvarado and C. DeMarco
2 iterations
Power flow solved
» mis
mis =
1.4496e-006
» VV
VV =
(1,1) 1.0500
(2,1) 1.0478- 0.0682i
(3,1) 1.0669- 0.0809i
(4,1) 0.9838- 0.0719i
(5,1) 0.9756- 0.0892i
(6,1) 0.9960- 0.1042i
» rect2pol(VV)
ans =
(1,1) 1.0500
(2,1) 1.0500- 3.7248i
(3,1) 1.0700- 4.3340i
(4,1) 0.9864- 4.1788i
(5,1) 0.9797- 5.2240i
(6,1) 1.0014- 5.9705i

The following comments apply:

The output from the code is minimal. The final mismatch and the VV output had to be manually extracted after the solution ended.
The final mismatch is quite small.
The rect2pol function gives the output as complex numbers, but the interpretation is that the real part of the numbers is the magnitude of VV and the imaginary part is the angle in degrees.

If, in addition, we execute the flows.m script:

» flows
Flows:
1 2 0.2773 0.1381 -1.0000 0.3098 1.0000
1 4 0.4161 -0.2164 -1.0000 0.4691 1.0000
1 5 0.3393 -0.1418 -1.0000 0.3678 1.0000
2 3 0.0349 0.0991 -1.0000 0.1051 1.0000
2 4 0.2856 -0.4919 -1.0000 0.5688 1.0000
2 5 0.1358 -0.1851 -1.0000 0.2296 1.0000
2 6 0.2417 -0.1615 -1.0000 0.2907 1.0000
3 5 0.1611 -0.2641 -1.0000 0.3094 1.0000
3 6 0.3609 -0.6319 -1.0000 0.7277 1.0000
4 5 0.0443 0.0206 -1.0000 0.0488 1.0000
5 6 0.0258 0.0908 -1.0000 0.0944 1.0000
Overloaded lines:
None

The columns here correspond to: the terminal nodes, the active and reactive flows, the minimum limit, the flow, and the maximum limit. There are no overloaded lines in this case.

And if, in addition, we execute the penalty.m script we obtain:

» penalty
Beta
1.0000 -0.1590
1.0154 -0.1777
1.0188 -0.1792
1.0289 -0.1460
1.0331 -0.1504
1.0264 -0.1702

The first column entries are the inverses of the penalty factors for the active powers at every node in the network. The second column entries are the corresponding factors for the reactive powers.