function [line_flow, bus]=loadlflow(line, bus, tolerance, max_iter, min_condition)
% [line_flow, bus_out]=loadlflow(line, bus, tolerance, max_iter, min_conditiion)
% Solve a load flow for a given electrical system (buses and lines)
%
% Usage:
%
% [line_flow, bus_out]=loadflow(line, bus)
% [line_flow, bus_out]=loadflow(line, bus, tolerance)
% [line_flow, bus_out]=loadflow(line, bus, tolerance, max_iter)
% [line_flow, bus_out]=loadflow(line, bus, tolerance, max_iter, min_condition)
%
% If un-specified defaults are:
% tolerance = 1e-5
% max_iter = 5
% min_condition = 1e-5
%
% Inputs:
% line / Line Matrix:
% Column Name Meaning
% 1 from ID of bus connected to one end
% 2 to ID of bus on the other end
% 3 resistance Resistance of the line in per unit (R)
% 4 reactance reactance of the line in per unit (X)
% 5 susceptance Susceptance of the line in per unit (line charging)
% Models the overall ground capacitative effects of
% a medium length line (50-150 miles). (B)
% Modeled as a pi-connected line See [1]
% 6 tap Transformer tap ratio (Step up and down ratio)
% 7 tapAngle Transformer phase shift in degrees
% 8 Pmax Max Real Power (Magnitude)
% 9 Smax Max Apparent Power (Magnitude)
% 10 factsSkip Skip line in FACTS Device placement
% 11 outageSkip Skip line outages (never outage)
%
% bus / Bus Matrix:
% Column Name Meaning
% 1 number ID of the bus (Sequential, with no Gaps 1:1:N)
% 2 voltage Voltage on the bus (initial)
% 3 theta Voltage Angle in degrees (initial)
% 4 Pgen P generated on the bus
% 5 Qgen Q generated on the bus
% 6 Pload P load on the bus
% 7 Qload Q load on the bus
% 8 G Shunt Conductance ("resistance" to ground)
% 9 B Shunt Susceptance ("reactance" to ground)
% 10 type Type of bus Matlab: 1=Swing, 2=Generator, 3=Normal
% Swing: theta and V are known (The reference bus)
% (P & Q are unknown)
% Generator: V and P are known
% (Q & theta are unknown)
% Normal/Load: P & Q are known
% (V & theta are unknown)
%
% tolerance: Convergence tolerance for power mis-matches
%
% max_iter: Maximum number of iterations before failure
%
% min_condition: Minimum acceptable reciprical-condition number estimate on Jacobian
%
% Return Values:
% line_flow / Line Flows:
% Column Name Meaning
% 1 from ID of bus connected to one end
% 2 to ID of bus on the other end
% 3 P_from2to Real Power going from "from bus" to "to bus"
% 4 P_to2from Real Power going from "to bus" to "from bus"
% 5 Q_from2to Reactive Power going from "from bus" to "to bus"
% 6 Q_to2from Reactive Power going from "to bus" to "from bus"
% 7 S_mag Maximum apparent power magnitued
% (max(|S(from->to)|, |S(to->From)|))
% Error Conditions:
% Message ID Error Message
% 'loadflow:ill_condition' 'Ill conditioned Jacobian'
% 'loadflow:max_iterations' 'Max Iterations Exceeded'
%
% Errors can be caught with "try/catch" blocks
% Ex:
% try
% [line_flow, bus_out]=loadflow(line,bus);
% ... % Analyze / record results
% catch
% err=lasterror;
% if err.identifier == 'loadflow:ill_condition'
% disp('Ill Conditioned Jacobian. Disregarding results.');
% end
% end
%
% Bus Model:
%
% ------------------------- <- The Bus (Voltage, Theta, P, Q)
% | | ^ |
% B - G \ | |
% - / | |
% | \ | |
% | / | |
% | | | V
% Gnd Gnd Pgen Pload
% Qgen Qload
% Each bus has a Voltage, Phase Angle (theta), Real Power, and Reactive power
% as well as a bus type. There are 3 types of buses:
% Type 1 - The "swing bus"
% The swing bus is assumed to be the "reference" for the phase angle.
% (I.e. it's theta is 0 by definition). This is always a form of
% Generator bus, and usually chosen as the generator with the largest
% output in the system. Since it's a generator, it's Voltage is also
% assumed to be known. Available real and reactive power are unknown.
% Type of bus Matlab: 1=Swing, 2=Generator, 3=Normal
% (theta and V are known / fixed, P & Q are unknown / solved for.
% Initial P and Q are "guesses" of correct values)
%
% Type 2 - The "PV buses" / Generator buses
% The PV buses represent all the remaining generators in the system
% Since they are generating and contributing power, the voltage and
% available power are assumed to be known and the phase angle and available
% reactive power are unknown (V and P are known / fixed, Q & theta are
% unknown / solved for. Initial P and Q are "guesses" of correct values)
%
% Type 3 - The "PQ buses" / Load buses
% The load buses represent the normal consumers of power (the loads).
% At the loads the real and reactive power needs are known, but the
% bus voltage and phase angle are unknown (P & Q are known / fixed,
% V & theta are unknown / solved for. Initial V and theta are "guesses"
% of correct values). At these buses, Pload is greater than Pgen
%
% B - The "susceptance" to ground. (Reciprical of reactance)
% This represents the losses due to the interaction of the ground
% and the AC signal on the bus.
%
% G - The "conductance" to ground. (Reciprical of resistance)
% This represents the losses due to AC signal to ground.
%
% Pgen / Qgen - These represent the real and reactive power being
% injected into the bus.
%
% Pload / Qload - These represent the real and reactive power being
% consumed at the bus.
%
% Line Model:
%
% tap @ tapAngle
% | |
% ---> 3---------/\/\/\/\-^^^^^^^-------------------|
% | 3 | R X | |
% |----------3 - - |
% | 3 - B/2 B/2 - |
% | 3 | | To Bus
% | | | |
% From Bus Gnd Gnd Gnd
%
% R - Resistance of the line in per unit
%
% X - Reactance of the line in per unit
%
% B - The susceptance (1/capacitance) of the line to ground (shunt)
% This is a pi-modeled line where the total capacitive effect
% to ground is modeled by a capacitance of half the total
% effective capacitance is placed on each side of the line.
% This is an accurate for a medium length lines (50-150 miles).
%
% tap - Multiplier to step up (or down) voltage for a transformer
%
% tapAngle - The phase shift imposed by a phase shifting transformer.
%
% References:
% 1. Power System Analysis: Short-Circuit Load Flow and Harmonics
% J.C. Das. Dekker. 2002.
% 2. Computational Methods for Electric Power Systems
% Mariesa Crow. CRC Press. 2003
% 3. Computer Analysis of Power System
% J. Arrillaga and C.P. Arnold. Wiley. 1994
% Copyright 2005 William M. Siever All Rights Reserved
% I. Initialize constants:
% A. Minimum allowed Jacobian reciporical condition number estimate
if nargin<5
min_condition=1e-5;
end
% B. Maximum allowed iterations
if nargin<4
max_iter=5;
end
% C. Tolerance of P and Q mismatch
if nargin<3
tolerance=1e-5;
end
% II. Do some simple data conversions
% A. Convert the bus angle to radians
bus(:,3)=bus(:,3)*pi/180;
% B. Convert the line angle to radians
line(:,7)=line(:,7)*pi/180;
% C. If transformer tap is 0, make it a 1
line(:,6)=line(:,6)+(line(:,6)==0);
% III. Setup for Newton-Raphson Iteration:
% A. Get the admittance Matrix
[y_mag, y_ang]=ymat(line,bus);
% B. Find the buses that will have a deriviative with respect to
% phase angle (Everything but the reference buses)
angleMap=find(bus(:,10)~=1);
% C. Find the buses that will have a deriviative with respect to
% voltage (the load buses)
vMap=find(bus(:,10)==3);
% D. Find the buses with generators
genMap=find(bus(:,10)~=3);
% E. Find the swing buses
swingMap=find(bus(:,10)==1);
% IV. Iteration
for i=1:max_iter
% A. Compute the Jacobian:
[jac, deltaP, deltaQ]=jmat(bus, y_mag, y_ang);
% B. Find complete P mismatch:
% Pmismatch = P(generation) - P(load) - deltaP
Pmismatch=bus(:,4)-bus(:,6)-deltaP;
% C. Find complete Q mismatch:
% Qmismatch = Q(generation) - Q(load) - deltaQ
Qmismatch=bus(:,5)-bus(:,7)-deltaQ;
% Update the Qgen on all generator buses
bus(genMap,5)=bus(genMap,5) + -Qmismatch(genMap);
% Update the Pgen on the swing bus
bus(swingMap,4)=bus(swingMap,4) + -Pmismatch(swingMap);
% D. Determine if current tolerance is acceptable
if all(abs(Pmismatch) < tolerance) & ...
all(abs(Qmismatch) < tolerance)
% A. Get From Voltages and Angles
Vf = bus(line(:,1),2);
Vf_ang = bus(line(:,1),3);
% B. Get To Voltages and Angles
Vt = bus(line(:,2),2);
Vt_ang = bus(line(:,2),3);
% C. Compute the From Transformer phase shift voltage multiplier
Vf_tps = Vf./line(:,6);
% D. Get each line's Admittance Magnitude and Angle
y_line_mag=(1./sqrt(line(:,3).^2+line(:,4).^2));
y_line_ang=-atan2(line(:,4),line(:,3));
% E. Compute the Sin and Cos terms of the Admittance Angles
cos_y_line_ang = cos(-y_line_ang);
sin_y_line_ang = sin(-y_line_ang);
% F. Compute the from and to angles
Vf_delta_ang = Vf_ang - line(:,7) - Vt_ang - y_line_ang;
Vt_delta_ang = Vt_ang + line(:,7) - Vf_ang - y_line_ang;
% G. Compute the real power flow
Pflow_f2t=Vf_tps.*y_line_mag.* ...
(Vf_tps.*cos_y_line_ang - Vt.*cos(Vf_delta_ang));
Pflow_t2f=Vt.*y_line_mag.* ...
(Vt.*cos_y_line_ang - Vf_tps.*cos(Vt_delta_ang));
% H. Compute the reactive power flow
Qflow_f2t=Vf_tps.*(y_line_mag.* ...
(Vf_tps.*sin_y_line_ang - Vt.*sin(Vf_delta_ang)) - ...
Vf_tps.*line(:,5)./2);
Qflow_t2f=Vt.*(y_line_mag.* ...
(Vt.*sin_y_line_ang - Vf_tps.*sin(Vt_delta_ang)) - ...
Vt.*line(:,5)./2);
% I. Consolidate the Real and Reactive power flows
Pflow=[Pflow_f2t Pflow_t2f];
Qflow=[Qflow_f2t Qflow_t2f];
% J. Compute the active power flow
Sflow=max(sqrt(Pflow.^2+Qflow.^2)')';
% K. Prepare the matrix of all lines and flows
line_flow=[line(:,1:2) Pflow Qflow Sflow];
% M. Convert the bus angle to degrees
bus(:,3)=bus(:,3)*180/pi;
% N. Convert the line angle to degrees
line(:,7)=line(:,7)*180/pi;
% Loadflow complete - return
return; % Only valid exit point
end
% E. Perform LU Decomposition & verify Jacobian condition
% (This uses an estimate of the reciprocal condition number)
[l,u,p,q] = lu(jac);
% [l,u,p] = lu(jac); %% BSIEVER: For Matlab Versions prior to 7
condition=full(min(abs(diag(u)))/max(abs(diag(u))));
if condition