Solving Power Flow Problems with a Matlab Implementation

of the Power System Applications Data Dictionary

Fernando L. Alvarado
alvarado@engr.wisc.edu
ECE Department, The University of Wisconsin, Madison, Wisconsin 53705

Abstract 2. The Data Dictionary

This paper implements a power flow application and The PSADD is an attempt by IEEE to extend and
variations using the IEEE Power System Application renovate the old standby of the power industry, the
Data Dictionary within a Matlab environment. It IEEE Common Format [1], which served the industry
describes a number of useful data and implementation well for many years. The Data Dictionary is meant
techniques for a variety of applications. The primarily as a definition of terms used for analytical
techniques include the use of very compact and applications within power system models. It is
efficient code for the computation of Power Transfer organized around the notion of objects. Within this
Distribution Factors. Power Transfer Distribution paper, we consider objects to be structures within the
Factors are an important element of present and Matlab environment. Each object contains a number
proposed congestion management strategies for of attributes or fields. As such, each object constitutes
power systems, particularly when these systems must a relational database. The following is a list of all data
operate in a deregulated environment. Keywords: types defined in the dictionary. Boldface indicates
databases, IEEE Common Format, data exchange. data types that are used in the present paper. Some of
the data types are defined but declared obsolescent
• Analog Data
• Area
1. Introduction • Bus Data
This paper provides a brief overview of the major • Capacitor Data: obsolete, see shunt device data.
elements of the IEEE Power System Application Data • Company Data
Dictionary (PSADD). It utilizes PSADD elements • Connectivity Node Data
toward the creation of flexible general-purpose power • Data Exchange Data
flow environments that have great expandability and • Equivalent AC Series Device Data
flexibility. The PSADD is intended primarily as a • Equivalent Injection/Shunt Device Data
means for exchanging application data among power • Generator Data: obsolete, see machine data.
system analytical applications of various types. The • Interchange Control Data
paper then uses an environment closely based on
• Line Data
PSADD concepts toward the creation of sample
applications in Matlab. The two main applications • Load Data
considered are: a simple “vanilla” power flow and a • Machine (Synchronous) Data
rapid and efficient computation of Power Transfer • Miscellaneous Application Data
Distribution Factors (PTDF). The paper is intended as • Ownership Data
a means of documenting efforts conducted by PSerc1 • Rating Data
toward the creation of a unified environment for • Reactor Data: obsolete, see shunt device data.
power flow analysis, as a means for illustrating and • Set Data
popularizing the Data Dictionary itself, and as a • Shunt (Nonsynchronous) Device Data
means for illustrating Matlab vectorized computation • Station Data
techniques and methods. The paper contains much • Status Data
actual code because it is intended for use in its • Switching Device Data
electronic form, where users are expected to be able to
• TieSwitch Data
“cut and paste” examples to a Matlab environment
• System
directly from the text of the paper.
• Terminal
1
The Power Systems Engineering Research Center, an • Transaction Data
NSF and industry sponsored university consortium. • Transformer Data
 abbreviation of one of the root
Each field of each object type has a type: integer, text, entities). Some possibilities are:
real or date. We only consider text or real types. In BS=Bus LN=Line
some cases, two alternative types are permitted in the TR=Xformer SW=Switch
MA=Machine LD=Load
PSADD for certain fields. We do not allow this. Only
SH=Shunt
a small subset of all available PSADD fields is used in Analog.DeviceRef Device assignment I
this paper. Fields or data types mentioned but not (Bus.Number, Line.Number) of
used are indicated with an asterisk*. T denotes Text value on device for device type.
(string) data, R denotes real data, I denotes integer Analog.Value Numeric value R
data. The Transformer and Load data types have
subtypes not described here. The next important category is the Bus type:
Name Definition
The PSADD integrates two viewpoints: a bus-oriented Bus.Number Bus number I
viewpoint common to analytical studies, and an Bus.Name Bus name T
equipment-oriented viewpoint common to measured Bus.BasekV Base voltage RMS R
values and system operation. The Miscellaneous Bus.Status I
Application Data describes this information. Here is a Bus.Vmag Magnitude of bus voltage. Solved R
subset of the Miscellaneous Application Data table: value.
Name Description Bus.Vangle Phase angle of bus voltage. R
Misc.Title User defined case title T Solved value
Misc.BaseMVA MVA base R Bus.MWMismatch Net Mw not assigned to R
generator, load or shunt
Misc.ModelConnect BM or Bus.Model T
ivityType BS or BusSection* Bus.MvarMismatc Net Mvar not assigned to R
BO or Both* h generator, load or shunt
Misc.ContingencyG RF = Use response factors T Bus.TypeCode *L or Load (no control) T
enMWResponse PC = Proportional to capacity *CB or ControlToBand
PO = Proportional to output *CT or ControlToTarget
*S = Swing (or reference)
Misc.DeviceRefMet Method (i.e., Number or Name) T
*I = Isolated
hod used as for referencing a device.
Bus.VoltTarget Voltage magnitude setpoint R
Bus.LowVoltLim Voltage control band lower limit R
Four other tables of generic interest are those for Bus.UpVoltLim Voltage control band upper limit R
Analog Data, Status Data, Ratings Data and Set Data. Bus.ControlPriorit T= transformer T
The Set Data is useful for grouping information into yCode S= Synchronous machine
Zones, Areas, Companies, Locations and almost any N =nonsynchronous shunt
other user-desired grouping. It is not used here. The
Ratings Data specifies rating types. This paper The nodal injections types used in this paper include
assumes all flow ratings are in MVA, thus no such Machine, Load and Shunt Element, defined next:
data is required. Status Data is used to denote the Name Description
status of a device or component. We use 1 to denote in Machine.Number Unique Number I
service, active or closed, 0 otherwise. Thus, we do not Machine.BusRef Bus.Number to which I
need Status Data as a separate data type. Analog Data machine is connected
is used to specify “metered” information (as opposed Machine.Id Identifier T
to “parameters” of components). Rather than entering Machine.Status I
this information as part of each device type, it is Machine.MW MW output, solved value R
entered as “analog data” with a pointer to the device Machine.Mvar MVAR output, solved value R
Machine.ControlMode PQ = hold P and Q constant I
type where it is used. Here is a subset of the
PV = adjust Q to hold V
definition of Analog Data: NL = adjust Q, no limits
Name Description SW = case swing Machine
Analog.Number Unique number I Machine.ControlledBu Bus (Bus.Number) whose I
Analog.Type VOLT=voltage magnitude T sRef voltage is controlled
FLOW=Power or current Machine.RatingId Identification of ratings T
MW = Active Power Machine.RatingValue A vector of possible generator R
MVAR = Reactive Power maximum MW values
TAP = Transformer Tap Machine.RatingMinim A vector of possible generator R
ANGL = Voltage Angle um(+) minimum MW values
Analog.AssignDev Device type associated with T Machine.MinQOutput Lower limit for reactive R
ice value (complete spelling or power
Machine.MaxQOutput Upper limit for reactive power R Shunt.BlockSuscept Number of discrete values I
Machine.MinOperatin Machine minimum operating R ValueNo available. Zero for continuous
gVolt voltage adjustment
Machine.MaxOperatin Machine maximum operating R Shunt.BlockSuscept Susceptance of each block A R
gVolt voltage ance vector.
Machine.PostContRes Percent of Machine MW R Shunt.MW Magnitude of Shunt output MW. R
pFact response to system MW A solved value
imbalance. Values are Shunt.Mvar Magnitude of Shunt output R
normalized for each island Mvar. A solved value
(+) Not in the current version of the Data Dictionary.
Two series PSADD elements are implemented here:
Name Description
Load.Number Unique Number I Name Definitions
Load.BusRef Bus.Number of bus to which load I Line.Number Unique number I
is connected Line.FromBusRef Bus.Number of "from" bus I
Load.Id Identifier T Line.ToBusRef Bus.Number of "to" bus I
Line.Circuit Alphanumeric circuit ID T
Load.Status I Line.Status I
Load.ModelType* CONST = constant MW and T Line.SectionR Series resistance R
MVAR Line.SectionX Series reactance R
VFUNC = voltage depend Line.ShuntConducta Total uniformly shunt R
INDMOT=induction motor nce conductance of this section
Load.ModelMinVo Minimum model voltage for this R Line.ShuntSuscepta Total shunt susceptance of this R
ltage* models nce section
Load.ModelMaxV Maximum voltage for this load R Line.FromMW Line MW flow at From end. A R
oltage* model solved value
Load.MinMW* Minimum real power R Line.FromMvar Line Mvar flow at From end. A R
Load.MaxMW* Maximum real power R solved value.
Load.ReferenceM Reference real power at one per R Line.ToMW Line MW flow at To end. A R
W* unit (for scaling) solved value
Load.MinMvar* Minimum reactive power R Line.ToMvar Line Mvar flow at To end. A R
consumed solved value
Load.MaxMvar* Maximum reactive power R Line.TypeId Type of equipment type (LINE, T
Load.ReferenceMv Reference reactive power at one R CAPACITOR*)
ar* per unit (for scaling) Line.Length* Length R
Load.Mvar Magnitude of Load Mvar. Solved R Line.LossAssignBus Bus.Number to which losses I
value Ref should be assigned
Load.Mw Magnitude of Load MW. Solved R Line.FromLumpedC Conductance of lumped shunt R
value ond element at "from" bus
Line.FromLumpedS Susceptance of lumped shunt R
Name Description uscept element at "from" bus
Shunt.Number Unique Number I Line.ToLumpedCon Conductance of lumped shunt R
Shunt.BusRef Bus.Number to which shunt is I d element at "to" bus
connected Line.ToLumpedSuc Susceptance of lumped shunt R
Shunt.Id Identifier I ept element at "to" bus
Shunt.Status I Line.RateId ID of line rating value. A vector T
Shunt.BaseMVA* MVA of this device R of text strings
Shunt.BasekV* Base voltage of this device R Line.RateValue Rating value (a vector) R
Shunt.Type Identifier of this device T
(Capacitor, Reactor) Name Definition
Shunt.ControlMode NONE = fixed shunt value I Xformer.Number Unique Number I
VOLT = hold voltage between Xformer.PrimaryBu Bus.Number to which winding I
limits sRef awinding A is connected
MVAR=hold MVAR flow* Xformer.Secondary Bus.Number to which winding I
Shunt.ControlBusRe Bus.Number at which voltage is I BusRef B is connected
f to be controlled Xformer.Circuit Identifier T
Shunt.MaxRegVolta Value at which shunt will R Xformer.Type FIXED Fix ratio and angle T
ge change out to lower voltage VOLT Fix angle var ratio
Shunt.MinRegVolta Value at which shunt will R MVAR Fix angle var ratio
ge change out to raise voltage MW – Fix ratio var angle
Xformer.TapRatio Tap ratio value R these types as a single structure with vector structure
Xformer.TapAngle Tap angle value R entries, or to represent each data type entry as a
Xformer.TapRatioM Minimum Tap ratio value of the R structure and to represent the entire data type as a
in transformer vector of structures. Either answer is correct and will
Xformer.TapRatioM Maximum Tap ratio value of the R work, but because of the vectorized nature of Matlab
ax transformer
computations, our experiments indicate that it is quite
Xformer.TapRatioSt Tap ratio step size value of the R
ep transformer
efficient to represent data types as structures of
Xformer.TapAngle Minimum Tap angle value of the R vectors. Furthermore, for the sake all consistence
Min transformer almost all vectors will be column vectors. When it is
Xformer.TapAngle Maximum Tap angle value of R necessary to represent sets (actually vectors) of
Max the transformer variable-length vectors, this is done by using the
Xformer.TapAngleS Tap angle step size value of the R notion of cells, which are vectors of arbitrary entry
tep transformer types. This is the convention adopted in the data
Xformer.Status I dictionary.
Xformer.BaseMVA * R
Xformer.FromMw Mw flow From bus. Solved R Start with a definition of all buses along with, loads,
value machines and shunts connected to them.
Xformer.FromMvar Mvar flow From bus. Solved R % Initialization of Bus Data Type
value Bus.Number=[]; Bus.Name=’’;
Xformer.ToMw Mw flow To bus. Solved value R Bus.BasekV=[]; Bus.Status=[];
Xformer.ToMvar Mvar flow To bus. Solved value R Bus.Vmag=[]; Bus.Vangle=[];
Xformer.LossAssign Bus.Number to which losses I Bus.MwMismatch=[]; Bus.MvarMismatch=[];
Bus.TypeCode=’’; Bus.VoltTarget=[];
BusRef should be assigned
Bus.LowVoltLim=[]; Bus.UpVoltLim=[];
Xformer.Pri-SecR s.c. resistance winding A to B R
Xformer.Pri-SecX s.c. reactance winding A to B R % Initialization of Machine Data Type
Xformer.Pri- s.c. zero sequence resistance R Machine.Number=[]; Machine.BusRef=[];
SecZeroR* winding A to B Machine.Id=’’; Machine.Status=[];
Xformer.Pri- s.c. zero sequence reactance R Machine.Mw=[]; Machine.Mvar=[];
Machine.ControlMode=’’;
SecZeroX* winding A to B Machine.ControlledBusRef=[];
Machine.MinQOutput=[]; Machine.MaxQOutput=[];
3. Rendition in Matlab Machine.MinOperatingVolt=[];
Machine.MaxOperatingVolt=[];
The implementation of the Data Dictionary in Matlab Machine.PostContRespFact=[];
is done by means of sparse arrays, structures and cells
% Initialization of Load Data Type
of types real and character. For the sake of clarity, all Load.Number=[]; Load.Id=’’;
concepts will be illustrated by means of an example. Load.BusRef=[]; Load.Status=[];
The example in question corresponds to the 6-bus Load.Mvar=[]; Load.Mw=[];
system from Wood and Wollenberg. The Common
% Initialization of Shunt Data Type
Format File (CFF) describing this data is illustrated in Shunt.Number=[]; Shunt.BusRef=[];
the appendix. The CFF is not sufficient to describe all Shunt.Id=’’; Shunt.Status=[];
the necessary information. Thus, additional Shunt.Type=’’; Shunt.ControlMode=[];
Shunt.ControlBusRef=’’;
assumptions will be made as needed. Shunt.MaxRegVoltage=[];
Shunt.MinRegVoltage=[];
To implement the data dictionary, one needs to first Shunt.BlockSusceptValueNo=[];
establish the necessary structures for all the base data Shunt.BlockSusceptance=[];
Shunt.Mw=[]; Shunt.Mvar=[];
types and the corresponding fields (or attributes) of
each data type. Some data types (such as Misc)
contain only scalar or simple text attributes. The Likewise, all structures for lines can be initialized:
Line.Number=[]; Line.FromBusRef=[];
complete definition for the Misc data type subset is: Line.ToBusRef=[]; Line.Circuit=’’;
Misc.Title=’This is a sample case’; Line.Status=[];
Misc.BaseMVA=100; Misc.ConnectivityType=’BM’; Line.SectionR=[]; Line.SectionX=[];
Line.ShuntConductance=[];
Misc.ContingencyGenMWResponse=’PC’; Line.ShuntSusceptance=[];
Misc.DeviceRefMethod=’Number’; Line.FromMw=[]; Line.FromMvar=[];
Line.ToMw=[]; Line.ToMvar=[];
Line.TypeId=’’; Line.LossAssignBusRef=[];
Most other data types include vector or even sets of Line.FromLumpedCond=[];
vectors as their data elements. The question arises Line.FromLumpedSuscept=[];
immediately in Matlab whether it is better to represent
Line.ToLumpedCond=[];Line.ToLumpedSuscept=[]; Load.BusRef=[4;5;6];
Line.RateId=’’; Line.RateValue=[]; Load.Id=strvcat(’1’,’1’,’1’);
Load.Status=[1;1;1];
Since no transformers are used in the example, details % Shunt data specification
of transformers are not illustrated. nsh=1; Shunt.Number=[1]; Shunt.BusRef=5;
Shunt.Id=’1’; Shunt.Status=1;
Shunt.Type=’Capacitor’;
At this point, the table values can be loaded for the Shunt.ControlMode=’NONE’;
example at hand more or less directly from the Data Shunt.ControlBusRef=5;
Dictionary. Some assumptions are made: Shunt.MaxRegVoltage=1.25*ones(nsh,1);
• All machine minimum powers are 10% of the Shunt.MinRegVoltage=0.85*ones(nsh,1);
Shunt.BlockSusceptValueNo={2}; % CELL ARRAY
machine maximum powers. Shunt.BlockSusceptValue={[0.1,0.1]}; % ”
• All machine maximum powers are 200 MW (or
2pu). % Line data specification
nL=11; Line.Number=(1:nL)’;
• We prefer to work in per-unit everywhere, even Line.FromBusRef=[1;1;1;2;2;2;2;3;3;4;5];
when the Data Dictionary calls for Mw or Mvar. Line.ToBusRef=[2;4;5;3;4;5;6;5;6;5;2];
• Minimum voltages anywhere are specified at 0.85 Line.Circuit=strvcat(’1’,’1’,’1’,’1’,...
’1’,’1’,’1’,’1’,’1’,’1’,’1’);
pu, maximum voltages at 1.2 pu. Line.SectionR=[.1;.05;.008;.05;.05; ...
• All shunt susceptances are assumed to be .1;.07;.12;.02;.02;.1];
adjustable in exactly two equal-sized steps. Since Line.SectionX=[.2;.2;.3;.25;.1;.3; ...
.2;.26;.1;.4;.3];
there are no susceptances in this example, a Line.ShuntConductance=zeros(nL,1);
capacitive susceptance with a value of 0.02 pu is Line.ShuntSusceptance=[.02;.02;.03; ...
assumed at bus 5. .02;.01;.02;.025;.025;.01;.04;.03];
Line.TypeId=char(ones(nL,1)*double(’Line’));
• Normally, data types would be constructed using Line.LossAssignBusRef=Line.FromBusRef;
a function that reads raw data from a database, Line.FromLumpedCond=zeros(nL,1);
from a common format data file, or some other Line.FromLumpedSuscept=zeros(nL,1);
available source. Here, data is constructed Line.ToLumpedCond=zeros(nL,1);
Line.ToLumpedSuscept=zeros(nL,1);
directly for the example at hand, without resorting ratid=strvcat(’NORMAL’,’EMERG’);
to a reading and translation function. ratval=[1.5;2];
for k=1:nL,
% Bus data specification Line.RateId{k,1}=ratid;
n=6; Bus.Number=(1:n)’; Line.RateValue{k,1}=ratval;
Bus.Name=strvcat(’BUS1’,’BUS2’,... end
’BUS3’,’BUS4’,’BUS5’,’BUS6’);
Bus.BasekV=138*ones(n,1); The specification of metered parameters for all data
Bus.Status=ones(n,1);
Bus.Vmag=[1.05;1.05;1.07;1;1;1]; types is done by means of the Analog data type. We
Bus.Vangle=zeros(n,1); consider first the specification of data values for all
Bus.TypeCode=... loads, followed by the specification of all machines,
strvcat(’S’,’CT’,’CT’,’L’,’L’,’L’); all shunts, all lines and all transformers (this example
Bus.VoltTarget=Bus.Vmag;
Bus.LowVoltLim=0.85*ones(n,1); has no transformers). Ordinarily this would also be
Bus.UpVoltLim=1.2*ones(n,1); accomplished by an input routine that would translate
appropriate database entries, common format file
% Machine data specification
ng=3; Machine.Number=(1:ng)’;
information and/or other such available data into the
Machine.Id=strvcat(Machine.Id,’1’,’1’,’1’); Data Dictionary structures. Here, the data is specified
Machine.Status=ones(ng,1); directly for the example of interest, starting with an
Machine.ControlMode=strvcat(’SW’,’PV’,’PV’); initialization of all Analog data types:
Machine.BusRef=[1;2;3];
%Initialization of all Analog data
Machine.ControlledBusRef=[1;2;3];
Analog.Number=[]; Analog.Type=’’;
Machine.RatingId=...
Analog.AssignDevice=’’;
strvcat(’NORMAL’,’NORMAL’,’NORMAL’);
Analog.DeviceRef=[]; Analog.Value=[];
Machine.RatingValue=[1.5;1.5;1.5];
Machine.RatingMinimum=[0.15;0.15;0.15];
% Specification of all Load operating Analog
Machine.MinQOutput=[-0.25;-0.25;-0.25];
data
Machine.MaxQOutput=[0.75;0.75;0.75];
for k=1:nl,
Machine.MinOperatingVolt=0.85*ones(ng,1);
na=length(Analog.Number);
Machine.MaxOperatingVolt=1.2*ones(ng,1);
Analog.Number=[Analog.Number;na+1;na+2];
Machine.PostContRespFact=[1;1;1];
Analog.Type=strvcat(Analog.Type,’MW’,’MVAR’);
Analog.AssignDevice=...
% Load data specification
strvcat(Analog.AssignDevice,’Load’,’Load’);
nl=3; Load.Number=(1:3)’;
Analog.DeviceRef=[Analog.DeviceRef;k;k];
end ept;
Analog.Value=[Analog.Value;0.7;0.7;... % Shunt devices at bus not considered here
0.7;0.7;0.7;0.7]; PQlist=strmatch(’L’,Bus.TypeCode);
PVlist=strmatch(’CT’,Bus.TypeCode);
% Specification of all Machine operating SlackList=strmatch(’S’,Bus.TypeCode);
Analog data GenList=union(SlackList,PVlist);
for k=1:ng, NonSlack=union(PQlist,PVlist);
na=length(Analog.Number);
Analog.Number=[Analog.Number;na+1];
Analog.Type=strvcat(Analog.Type,’MW’); Next we construct the Y-bus matrix:
Analog.AssignDevice = ... n=length(Bus.Number);
strvcat(Analog.AssignDevice,’Machine’); Ybus=sparse([],[],[],n,n);
end Ybus=Ybus+sparse(Line.FromBusRef,Line.FromBus
Analog.DeviceRef=[Analog.DeviceRef;1;2;3]; Ref,1./Line.Z+Line.ShuntY2,n,n);
Analog.Value=[Analog.Value;0;0.5;0.6]; Ybus=Ybus+sparse(Line.FromBusRef,Line.ToBusRe
f,-1./Line.Z,n,n);
Ybus=Ybus+sparse(Line.ToBusRef,Line.FromBusRe
4. A “Vanilla” Newton Power Flow f,-1./Line.Z,n,n);
Ybus=Ybus+sparse(Line.ToBusRef,Line.ToBusRef,
The above structures can be used directly in the 1./Line.Z+Line.ShuntY2,n,n);
implementation of a basic Newton power flow. Ybus=Ybus+sparse(Line.FromBusRef,Line.FromBus
Convenient computation in Matlab requires a few Ref,Line.YI,n,n);
Ybus=Ybus+sparse(Line.ToBusRef,Line.ToBusRef,
additional ideas: Line.YJ,n,n);
• Most data types used for computation should be
complex. Thus, a few data types are added to the We then illustrate the code to construct the complete
definition of the Bus data type (and other data power flow Jacobian for any system (Note: this highly
types as needed). vectorized version of the code is due to Christopher
• An Internal data type is defined and used to DeMarco of the University of Wisconsin):
contain all important intermediate matrix function [dSdd,dSdv]=pflowjac(Yb,Vb)
variables. Ib=Yb*Vb;
• An overall encapsulation of all data types into a dSdd=j*diag(conj(Ib).*Vb)-
j*diag(Vb)*conj(Yb)*diag(conj(Vb));
single master data types is done to be able to refer dSdv=diag(conj(Ib).*(Vb./abs(Vb)))
to all the variables and data for a single study by +diag(Vb)*conj(Yb)*diag(conj(Vb)./abs(Vb));
means of a single variable. This is a structure of
structures. Finally, the complete master vanilla Newton solver is:
Vmag=Bus.Vmag; Vang=Bus.Vangle*pi/180;
The first step is to convert pertinent Data Dictionary for iter=1:10,
V=Vmag.*exp(j*Vang);
structures into Matlab structures suitable for [dSdd,dSdv]=pflowjac(Ybus,V);
computation. Also, it is extremely useful to associate misvect=V.*conj(Ybus*V)+Bus.SL-Bus.PM;
any analog measurement data from the Analog data rmiss=[real(misvect(NonSlack)); ...
types and associate it with the pertinent data type in a imag(misvect(PQlist))];
mismatch=max(abs(rmiss));
manner suitable for computation. The following code if mismatch<0.0001, break; end
performs these tasks and further encapsulates all rjac=[real(dSdd(NonSlack,NonSlack)) ...
information: real(dSdv(NonSlack,PQlist)); ...
j=sqrt(-1); imag(dSdd(PQlist,NonSlack)) ...
kPL=intersect(strmatch(’Load’,Analog.AssignDe imag(dSdv(PQlist,PQlist))];
vice),strmatch(’MW’,Analog.Type)); dx=rjac\rmiss;
kQL=intersect(strmatch(’Load’,Analog.AssignDe dxAng=dx(1:length(NonSlack));
vice),strmatch(’MVAR’,Analog.Type)); dxMag=dx(length(NonSlack)+1:length(dx));
kPM=intersect(strmatch(’Machine’,Analog.Assig Vang(NonSlack)=Vang(NonSlack)-dxAng;
nDevice),strmatch(’MW’,Analog.Type)); Vmag(PQlist)=Vmag(PQlist)-dxMag;
Bus.SL=zeros(n,1); Bus.PM=zeros(n,1); end
iPL=Load.BusRef(Analog.DeviceRef(kPL)); Bus.Vmag=Vmag; Bus.Vangle=Vang*180/pi;
iQL=Load.BusRef(Analog.DeviceRef(kQL));
iPM=Machine.BusRef(Analog.DeviceRef(kPM));
Bus.SL(iPL)=Bus.SL(iPL)+Analog.Value(kPL);
To run the example in this paper, “cut and paste”
Bus.SL(iQL)=Bus.SL(iQL)+j*Analog.Value(kQL); every Matlab code segment to a Matlab environment.
Bus.PM(iPM)=Bus.PM(iPM)+Analog.Value(kPM); Store function code segments in m-files.
Line.Z=Line.SectionR+j*Line.SectionX;
Line.ShuntY2=(Line.ShuntConductance+j*Line.Sh
untSusceptance)/2; Extrapolation of these techniques to considerably
Line.YI=Line.FromLumpedCond+j*Line.FromLumped more complex cases and problems is immediate. We
Suscept; illustrate only one such extension.
Line.YJ=Line.ToLumpedCond+j*Line.ToLumpedSusc
5. Power Transfer Distribution Factors At this point, the Matlab environment contains already
a wealth of useful information for answering a large
As an example of the usefulness of this approach to
number of questions about the system.
practical power system computations, we illustrate the
code to compute Power Transfer Distribution Factors
for an increase in active power demand at any location 6. Conclusions
in the network, relative to the power being supplied at Power System Application Data Dictionary data types
the slack generator. By having the PTDFs for any two have been used as the basis for a direct
(or more) locations in the grid, the PTDF for any implementation of a simple power flow in Matlab.
bilateral (or multilateral) transaction not involving the The implementation is remarkably efficient and easy
slack generator can also be readily found, although to extend to any desired additional purpose: new
this is not done here. For the example above, the methods can be tested with ease, new functions and
PTDF code requires the determination of a Jacobian objectives implemented, all while maintaining the
matrix that relates the flows at either end of a line to ability to communicate between and among diverse
changes in voltage magnitudes and angles. The applications in a standardized manner.
vectorized code to construct such matrix for the flow
at the “from” end of every line in the system is given Acknowledgements
next. Here Vb is the vector of complex nodal voltages
Professor Christopher DeMarco made available the
and Vs refers to the sending end voltages of the lines
key Jacobian routine to this author. Joann Staron has
of interest. Yf is either YflowI or YflowJ, the been for years instrumental in developing the PSADD.
from or to end “flow” admittance matrices. Thanks are also due to Scott Greene and Taiyou Yong.
function [YflowI,YflowJ]=bldyflow(Bus,Line)
n=length(Bus.Number);
nb=length(Line.FromBusRef); nbe=(1:nb)'; References
YflowI=sparse(nbe,Line.FromBusRef,1./Line.Z+i
*Line.ShuntY2+Line.YI,nb,n); [1] IEEE Committee Report, “Common format for
YflowI=YflowI+sparse(nbe,Line.ToBusRef,- exchange of solved load flow data,” IEEE
1./Line.Z,nb,n); Transactions on Power Apparatus and Systems,
YflowJ=sparse(nbe,Line.FromBusRef,-
1./Line.Z,nb,n); Volume PAS-92, Number 6, pages 1916-1925,
YflowJ=YflowJ+sparse(nbe,Line.ToBusRef,1./Lin November/December, 1973.
e.Z+i*Line.ShuntY2+Line.YJ,nb,n);

function [dFdd,dFdv]=flowjac(Yf,Line,Vb,Vs)
nb=length(Line.FromBusRef); n=length(Vb);
If=Yflow*Vb;
dFdd=j*sparse((1:nb)',Line.FromBusRef,
(conj(If).*Vs),nb,n)-
j*diag(Vs)*conj(Yf)*diag(conj(Vb));
dFdv=sparse((1:nb)',Line.FromBusRef,
(conj(If).*(Vs./abs(Vs))),nb,n)
+diag(Vs)*conj(Yf)*diag(conj(Vb)./abs(Vb));

This code can be used to determine the sensitivity of

flows on any line to any injections at a given operating
point. Sample code to do this calculation for our
example for an injection at bus 5 would be:
[YflowI,YflowJ]=bldyflow(Bus,Line);
[dFdd,dFdv]=flowjac(YflowI,Line,V,V(Line.From
BusRef));
[dSdd,dSdv]=pflowjac(Ybus,V);
rJ=[real(dSdd(NonSlack,NonSlack)),real(dSdv(N
onSlack,PQlist));imag(dSdd(PQlist,NonSlack)),
imag(dSdv(PQlist,PQlist))];
rJf=[real(dFdd(:,NonSlack)),real(dFdv(:,PQlis
t));imag(dFdd(:,NonSlack)),imag(dFdv(:,PQlist
))];
e=zeros(n,1); e(5)=1;
er=[real(e(NonSlack));imag(e(PQlist))];
dFP5=rJf*(rJ\er)