A Power Flow in Matlab

A Power Flow in Matlab © 1999, 2000, 2001 Fernando L. Alvarado 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 [pic]. 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: [pic] where [pic] denotes element-wise product of two vectors (using Matlab
notation), and * as a superscript denotes the complex conjugate. The dot
[pic]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, ?i is the angle of Vi, ?i 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: [pic] and from here the Jacobian can be determined from: [pic] 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.S
tatus;
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 nargin0.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 p-code 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 p-code 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