Power Flow Solution
Power Flow Solution
Power Flow Solution
Use the following data to solve the power flow problem using Newton-Raphson method. Take bus 1
as the slack bus and base MVA as 100.
Build the system with the following data in PowerWorld software. Present the power flow result
from PowerWorld software, using Newton-Raphson method. How does it compare with the
computed power flow solution?
Is the generator at the slack bus producing or consuming reactive power? What can you do to make
sure that it does not consume any reactive power? Verify your solution with the software, and show
the results.
Data for the 5 bus system
Bus data
Number
Nom kV
PU Volt
Angle (Deg)
Load MW
Load Mvar
Gen
MW
Gen
Mvar
B Shunt
Mvar
15
394.8
114.18
345
700
180
345
1.05
80
40
520
337.35
345
345
Line data
From Number
To Number
Transformer
0.0015
0.02
Line
0.009
0.1
1.72
Line
0.0045
0.05
0.88
Transformer
0.00075
0.01
Line
0.00225
0.025
0.44
Acknowledgement:
We would like to thank Joe H. Chow for code for finding inverse of matrix.
CODES USED:
1.Calc.m
function [delP,delQ,P,Q,conv_flag] = ...
calc(nbus,V,ang,Y,Pg,Qg,Pl,Ql,sw_bno,g_bno,tol)
% Syntax: [delP,delQ,P,Q,conv_flag] =
%
calc(nbus,V,ang,Y,Pg,Qg,Pl,Ql,sw_bno,g_bno,tol)
%
% Purpose: calculates power mismatch and checks convergence
%
also determines the values of P and Q based on the
%
supplied values of voltage magnitude and angle
% Version: 2.0 eliminates do loop
% Input:
nbus
- total number of buses
%
bus_type -load_bus(3), gen_bus(2), swing_bus(1)
%
V
- magnitude of bus voltage
%
ang
- angle(rad) of bus voltage
%
Y
- admittance matrix
%
Pg
- real power of generation
%
Qg
- reactive power of generation
%
Pl
- real power of load
%
Ql
- reactive power of load
%
sw_bno - a vector having zeros at all swing_bus locations ones
otherwise
%
g_bno - a vector having zeros at all generator bus locations ones
otherwise
%
tol
- a tolerance of computational error
%
% Output: delP
- real power mismatch
%
delQ
- reactive power mismatch
%
P
- calculated real power
%
Q
- calculated reactive power
%
conv_flag - 0, converged
%
1, not yet converged
%
% ************************************************************
jay = sqrt(-1);
swing_bus = 1;
gen_bus = 2;
load_bus = 3;
% voltage in rectangular coordinate
V_rect = V.*exp(jay*ang);
% bus current injection
cur_inj = Y*V_rect;
% power output based on voltages
S = V_rect.*conj(cur_inj);
P = real(S); Q = imag(S);
delP = Pg - Pl - P;
delQ = Qg - Ql - Q;
% zero out mismatches on swing bus and generation bus
delP=delP.*sw_bno;
delQ=delQ.*sw_bno;
delQ=delQ.*g_bno;
% total mismatch
[pmis,ip]=max(abs(delP));
[qmis,iq]=max(abs(delQ));
mism = pmis+qmis;
ifmism>tol,
conv_flag = 1;
else
conv_flag = 0;
end
return
2. chq_lim
function f = chq_lim(qg_max,qg_min)
% function for detecting generator vars outside limit
% sets Qg to zero if limit exceded, sets Ql to negative of limit
% sets bus_type to 3, and recalculates ang_red and volt_red
% changes generatorbus_type to type 3
% recalculates the generator index
% inputs: qg_max and qg_min are the last two clumns of the bus matrix
3.data5bus
bus = [
1 1.0 0
2 1.0 0
3 1.05 0
4 1.00 0
5 1.00 0
];
3.948
0.0
5.2
0.0
0.0
1.1418
0.0
3.3735
0.0
0.0
0.0
7
0.8
0.0
0.0
0.0
1.8
0.4
0.0
0.0
0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
1
3
2
3
3
99.0
0.0
99.0
0.0
0.0
-99
0
-99
0
0
15
345
345
345
345
1.1
1.1
1.1
1.1
1.1
.9;
.9;
.9;
.9;
.9;
line = [
5 1 0.0015
0.02
0.0
1 0.0 1.1 0.9 0.05;
4 2 0.009
0.1
1.72
0 0.0 0
0.0 0.0;
5 2 0.0045
0.05
0.88
0 0.0 0
0.0 0.0;
3 4 0.00075
0.01
0.0
1
0.0 1.1 0.9 0.05;
5 4 0.00225
0.025 0.44
0
0.0 0.0 0.0 0.00;
];
4.form_jac
function [Jac11,Jac12,Jac21,Jac22]=form_jac(V,ang,Y,ang_red,volt_red)
% Syntax: [Jac] = form_jac(V,ang,Y,ang_red,volt_red)
%
[Jac11,Jac12,Jac21,Jac22] = form_jac(V,ang,Y,...
%
ang_red,volt_red)
% Purpose: form the Jacobian matrix using sparse matrix techniques
% Input:
V
- magnitude of bus voltage
%
ang
- angle(rad) of bus voltage
%
Y
- admittance matrix
%
ang_red - matrix to eliminate swing bus voltage magnitude and
angle
%
entries
%
volt_red - matrix to eliminate generator bus voltage magnitude
%
entries
% Output: Jac
- jacobian matrix
%
Jac11,Jac12,Jac21,Jac22 - submatrices of
%
jacobian matrix
% ***********************************************************
jay = sqrt(-1);
exp_ang = exp(jay*ang);
% Voltage rectangular coordinates
V_rect = V.*exp_ang;
CV_rect=conj(V_rect);
Y_con = conj(Y);
%vector of conjugate currents
i_c=Y_con*CV_rect;
% complex power vector
S=V_rect.*i_c;
S=sparse(diag(S));
Vdia=sparse(diag(V_rect));
CVdia=conj(Vdia);
Vmag=sparse(diag(abs(V)));
S1=Vdia*Y_con*CVdia;
t1=((S+S1)/Vmag)*volt_red';
t2=(S-S1)*ang_red';
J11=-ang_red*imag(t2);
J12=ang_red*real(t1);
J21=volt_red*real(t2);
J22=volt_red*imag(t1);
ifnargout> 3
Jac11 = J11; clear J11
Jac12 = J12; clear J12
Jac21 = J21; clear J21
Jac22 = J22; clear J22
else
end
end
if ~isempty(min_v_idx)
% change on load taps to correct voltage
% assumes that bus to be corrected is the to bus
for fb = 1 : length(min_v_idx)
chk_fb = find(line(:,2) == bus(min_v_idx(fb),1));
if ~isempty(chk_fb)
% freeze dc taps
if ~isempty(n_dcl)
forkt = 1:2*n_dcl
dc_chk = find(chk_fb==ac_line(kt));
if ~isempty(dc_chk);chk_fb(dc_chk)=[];end
end
end
end
if ~isempty(chk_fb)
if line(chk_fb,8)~=0
% can change tap
disp('voltage low changing tap on line');disp(chk_fb)
vm2 = volt_min(min_v_idx(fb));
verror = vm2 - V(min_v_idx(fb));
% voltage too low tap must be reduced
ifverror/vm2<line(chk_fb,10)
% reduce tap by one increment
tap = line(chk_fb,6)-line(chk_fb,10);
else
tap = line(chk_fb,6) - verror./vm2;
tap_set = fix( (tap-line(chk_fb,9))/line(chk_fb,10));
tap = tap_set*line(chk_fb,10) + line(chk_fb,9);
end
tap = min(line(chk_fb,8),max(tap, line(chk_fb,9)));
disp('taps reset to');tap
% reset tap in line data
line(chk_fb,6) = tap;
end
end
end
end
end
7.loadflow.m
function [bus_sol,line_sol,line_flow] = ...
loadflow(bus,line,tol,iter_max,acc,display,flag)
% Syntax:
[bus_sol,line_sol,line_flow] =
% loadflow(bus,line,tol,iter_max,acc,display,flag)
% Purpose:
solve the load-flow equations of power systems
%
modified to eliminate do loops and improve the use
%
sparse matices
% (c) Copyright 1991 Joe H. Chow - All Rights Reserved
% Input:
%
%
%
%
%
%
%
%
%
bus
line
tol
iter_max
acc
display
flag
bus data
line data
tolerance for convergence
maximum number of iterations
acceleration factor
'y', generate load-flow study report
else, no load-flow study report
- 1, form new Jacobian every iteration
2, form new Jacobian every other
iteration
% Output:
bus_sol
- bus solution (see report for the
%
solution format)
%
line_sol - modified line matrix
%
line_flow - line flow solution (see report)
% Algorithm: Newton-Raphson method using the polar form of
%
the equations for P(real power) and Q(reactive power).
% Calls:
Y_sparse, calc, form_jacchq_lim
% ***********************************************************
globalbus_int
globalQgbus_typeg_bnoPQV_noPQ_noang_redvolt_red
global Q Ql
globalgen_chg_idx
globalac_linen_dcl
tt = clock;
% start the total time clock
jay = sqrt(-1);
load_bus = 3;
gen_bus = 2;
swing_bus = 1;
if exist('flag') == 0
flag = 1;
end
lf_flag = 1;
% set solution defaults
ifisempty(tol);tol = 1e-9;end
ifisempty(iter_max);iter_max = 30;end
ifisempty(acc);acc = 1.0; end;
ifisempty(display);display = 'n';end;
if flag <1 | flag > 2
error('LOADFLOW: flag not recognized')
end
[nlinenlc] = size(line);
% number of lines and no of line cols
[nbusncol] = size(bus);
% number of buses and number of col
% set defaults
% bus data defaults
ifncol<15
% set generator var limits
ifncol<12
bus(:,11) = 9999*ones(nbus,1);
bus(:,12) = -9999*ones(nbus,1);
end
ifncol<13;bus(:,13) = ones(nbus,1);end
bus(:,14) = 1.5*ones(nbus,1);
bus(:,15) = 0.5*ones(nbus,1);
volt_min = bus(:,15);
volt_max = bus(:,14);
else
volt_min = bus(:,15);
volt_max = bus(:,14);
end
no_vmin_idx = find(volt_min==0);
if ~isempty(no_vmin_idx)
volt_min(no_vmin_idx) = 0.5*ones(length(no_vmin_idx),1);
end
no_vmax_idx = find(volt_max==0);
if ~isempty(no_vmax_idx)
volt_max(no_vmax_idx) = 1.5*ones(length(no_vmax_idx),1);
end
no_mxv = find(bus(:,11)==0);
no_mnv = find(bus(:,12)==0);
if ~isempty(no_mxv);bus(no_mxv,11)=9999*ones(length(no_mxv),1);end
if ~isempty(no_mnv);bus(no_mnv,12) = -9999*ones(length(no_mnv),1);end
no_vrate = find(bus(:,13)==0);
if ~isempty(no_vrate);bus(no_vrate,13) = ones(length(no_vrate),1);end
tap_it = 0;
tap_it_max = 10;
no_taps = 0;
% line data defaults, sets all tap ranges to zero - this fixes taps
ifnlc< 10
line(:,7:10) = zeros(nline,4);
no_taps = 1;
% disable tap changing
end
% outer loop for on-load tap changers
mm_chk=1;
while (tap_it<tap_it_max&mm_chk)
tap_it = tap_it+1;
% build admittance matrix Y
[Y,nSW,nPV,nPQ,SB] = y_sparse(bus,line);
% process bus data
bus_no = bus(:,1);
V = bus(:,2);
ang = bus(:,3)*pi/180;
Pg = bus(:,4);
Qg = bus(:,5);
Pl = bus(:,6);
Ql = bus(:,7);
Gb = bus(:,8);
Bb = bus(:,9);
bus_type = round(bus(:,10));
qg_max = bus(:,11);
qg_min = bus(:,12);
sw_bno=ones(nbus,1);
g_bno=sw_bno;
% set up index for Jacobian calculation
%% form PQV_no and PQ_no
bus_zeros=zeros(nbus,1);
swing_index=find(bus_type==1);
sw_bno(swing_index)=bus_zeros(swing_index);
PQV_no=find(bus_type>=2);
PQ_no=find(bus_type==3);
gen_index=find(bus_type==2);
g_bno(gen_index)=bus_zeros(gen_index);
%sw_bno is a vector having ones everywhere but the swing bus locations
%g_bno is a vector having ones everywhere but the generator bus
locations
% construct sparse angle reduction matrix
il = length(PQV_no);
ii = [1:1:il]';
ang_red = sparse(ii,PQV_no,ones(il,1),il,nbus);
% construct sparse voltage reduction matrix
il = length(PQ_no);
ii = [1:1:il]';
volt_red = sparse(ii,PQ_no,ones(il,1),il,nbus);
iter = 0;
[delP,delQ,P,Q,conv_flag] =...
calc(nbus,V,ang,Y,Pg,Qg,Pl,Ql,sw_bno,g_bno,tol);
st = clock;
% start the iteration time clock
%% start iteration process for main Newton_Raphson solution
while (conv_flag == 1 &iter<iter_max)
iter = iter + 1;
% Form the Jacobean matrix
clearJac
Jac=form_jac(V,ang,Y,ang_red,volt_red);
% reduced real and reactive power mismatch vectors
red_delP = ang_red*delP;
red_delQ = volt_red*delQ;
cleardelPdelQ
% solve for voltage magnitude and phase angle increments
temp = Jac\[red_delP; red_delQ];
% expand solution vectors to all buses
delAng = ang_red'*temp(1:length(PQV_no),:);
delV = volt_red'*temp(length(PQV_no)+1:length(PQV_no)+length(PQ_no),:);
% update voltage magnitude and phase angle
V = V + acc*delV;
V = max(V,volt_min); % voltage higher than minimum
V = min(V,volt_max); % voltage lower than maximum
ang = ang + acc*delAng;
% calculate the power mismatch and check convergence
[delP,delQ,P,Q,conv_flag] =...
calc(nbus,V,ang,Y,Pg,Qg,Pl,Ql,sw_bno,g_bno,tol);
% check if Qg is outside limits
gen_index=find(bus_type==2);
Qg(gen_index) = Q(gen_index) + Ql(gen_index);
lim_flag = chq_lim(qg_max,qg_min);
iflim_flag == 1;
disp('Qg at var limit');
end
end
ifiter == iter_max
imstr = int2str(iter_max);
disp(['inner ac load flow failed to converge after ', imstr,' iterations'])
tistr = int2str(tap_it);
disp(['at tap iteration number ' tistr])
else
disp('inner load flow iterations')
disp(iter)
end
ifno_taps == 0
lftap
else
mm_chk = 0;
end
end
iftap_it>= tap_it_max
titstr = int2str(tap_it_max);
disp(['tap iteration failed to converge after',titstr,' iterations'])
else
disp(' tap iterations ')
disp(tap_it)
end
ste = clock;
% end the iteration time clock
vmx_idx = find(V==volt_max);
vmn_idx = find(V==volt_min);
if ~isempty(vmx_idx)
disp('voltages at')
bus(vmx_idx,1)'
disp('are at the max limit')
end
if ~isempty(vmn_idx)
disp('voltages at')
bus(vmn_idx,1)'
disp('are at the min limit');
end
gen_index=find(bus_type==2);
load_index = find(bus_type==3);
Pg(gen_index) = P(gen_index) + Pl(gen_index);
Qg(gen_index) = Q(gen_index) + Ql(gen_index);
gend_idx = find((bus(:,10)==2)&(bus_type~=2));
if ~isempty(gend_idx)
disp('the following generators are at their var limits')
disp('
bus#
Qg')
disp([bus(gend_idx,1) Q(gend_idx)])
Qlg = Ql(gend_idx)-bus(gend_idx,7);% the generator var part of the load
Qg(gend_idx)=Qg(gend_idx)-Qlg;% restore the generator vars
Ql(gend_idx)=bus(gend_idx,7);% restore the original load vars
end
Pl(load_index) = Pg(load_index) - P(load_index);
Ql(load_index) = Qg(load_index) - Q(load_index);
Pg(SB) = P(SB) + Pl(SB); Qg(SB) = Q(SB) + Ql(SB);
VV = V.*exp(jay*ang); % solution voltage
% calculate the line flows and power losses
tap_index = find(abs(line(:,6))>0);
tap_ratio = ones(nline,1);
tap_ratio(tap_index)=line(tap_index,6);
phase_shift(:,1) = line(:,7);
tps = tap_ratio.*exp(jay*phase_shift*pi/180);
from_bus = line(:,1);
from_int = bus_int(round(from_bus));
to_bus = line(:,2);
to_int = bus_int(round(to_bus));
r = line(:,3);
rx = line(:,4);
chrg = line(:,5);
z = r + jay*rx;
y = ones(nline,1)./z;
MW_s = VV(from_int).*conj((VV(from_int) - tps.*VV(to_int)).*y ...
+ VV(from_int).*(jay*chrg/2))./(tps.*conj(tps));
P_s = real(MW_s);
% active power sent out by from_bus
% to to_bus
Q_s = imag(MW_s);
% reactive power sent out by
% from_bus to to_bus
MW_r = VV(to_int).*conj((VV(to_int) ...
- VV(from_int)./tps).*y ...
+ VV(to_int).*(jay*chrg/2));
P_r = real(MW_r);
% active power received by to_bus
% from from_bus
Q_r = imag(MW_r);
% reactive power received by
% to_bus from from_bus
iline = [1:1:nline]';
line_ffrom = [ilinefrom_busto_bus P_s Q_s];
line_fto
= [ilineto_busfrom_busP_rQ_r];
% keyboard
P_loss = sum(P_s) + sum(P_r) ;
Q_loss = sum(Q_s) + sum(Q_r) ;
bus_sol=[bus_no V ang*180/pi PgQg Pl Ql Gb Bb...
bus_typeqg_maxqg_minbus(:,13) volt_maxvolt_min];
line_sol = line;
line_flow(1:nline, :) =[ilinefrom_busto_bus P_s Q_s];
line_flow(1+nline:2*nline,:) = [ilineto_busfrom_busP_rQ_r];
% Give warning of non-convergence
ifconv_flag == 1
disp('ac load flow failed to converge')
error('stop')
end
% display results
if display == 'y',
clc
disp('
LOAD-FLOW STUDY')
disp('
REPORT OF POWER FLOW CALCULATIONS ')
disp(' ')
disp(date)
fprintf('SWING BUS
: BUS %g \n', SB)
fprintf('NUMBER OF ITERATIONS
: %g \n', iter)
fprintf('SOLUTION TIME
: %g sec.\n',etime(ste,st))
fprintf('TOTAL TIME
: %g sec.\n',etime(clock,tt))
fprintf('TOTAL REAL POWER LOSSES
: %g.\n',P_loss)
fprintf('TOTAL REACTIVE POWER LOSSES: %g.\n\n',Q_loss)
ifconv_flag == 0,
disp('
GENERATION
LOAD')
disp('
BUS
VOLTS
ANGLE
REAL REACTIVE
REAL
REACTIVE ')
disp(bus_sol(:,1:7))
disp('
LINE FLOWS
')
disp('
LINE FROM BUS
TO BUS
REAL REACTIVE
')
disp(line_ffrom)
disp(line_fto)
end
end; %
ifiter>iter_max,
disp('Note: Solution did not converge in %g iterations.\n', iter_max)
lf_flag = 0
end
return
8.y_sparse.m
function
[Y,nSW,nPV,nPQ,SB] = y_sparse(bus,line)
% Syntax:
[Y,nSW,nPV,nPQ,SB] = y_sparse(bus,line)
%
% Purpose:
build sparse admittance matrix Y from the line data
%
% Input:
bus - bus data
%
line - line data
%
% Output:
Y
- admittance matrix
%
nSW - total number of swing buses
%
nPV - total number generator buses
%
nPQ - total number of load buses
%
SB
- internal bus numbers of swing bus
% ************************************************************
globalbus_int
jay = sqrt(-1);
swing_bus = 1;
gen_bus = 2;
load_bus = 3;
nline = length(line(:,1));
nbus = length(bus(:,1));
r=zeros(nline,1);
rx=zeros(nline,1);
chrg=zeros(nline,1);
z=zeros(nline,1);
y=zeros(nline,1);
% number of lines
% number of buses
Y = sparse(1,1,0,nbus,nbus);
% set up internal bus numbers for second indexing of buses
busmax = max(bus(:,1));
bus_int = zeros(busmax,1);
ibus = [1:1:nbus]';
fori = 1:nbus
bus_int(round(bus(i,1))) = i;
end
% process line data and build admittance matrix Y
r = line(:,3);
rx = line(:,4);
chrg =jay*sparse(diag( 0.5*line(:,5)));
z = r + jay*rx;
% line impedance
y = sparse(diag(ones(nline,1)./z));
% determine connection matrices including tap changers and phase shifters
from_bus = round(line(:,1));
from_int = bus_int(from_bus);
to_bus = round(line(:,2));
to_int = bus_int(to_bus);
tap_index = find(abs(line(:,6))>0);
tap=ones(nline,1);
tap(tap_index)=1. ./line(tap_index,6);
phase_shift = line(:,7);
tap = tap.*exp(-jay*phase_shift*pi/180);
% sparse matrix formulation
iline = [1:1:nline]';
C_from = sparse(from_int,iline,tap,nbus,nline,nline);
C_to = sparse(to_int,iline,ones(nline,1),nbus,nline,nline);
C_line = C_from - C_to;
% form Y matrix from primative line ys and connection matrices
Y=C_from*chrg*C_from' + C_to*chrg*C_to' ;
Y = Y + C_line*y*C_line';
Gb = bus(:,8);
% bus conductance
Bb = bus(:,9);
% bus susceptance
% add diagonal shunt admittances
Y = Y + sparse(ibus,ibus,Gb+jay*Bb,nbus,nbus);
ifnargout> 1
% count buses of different types
nSW = 0;
nPV = 0;
nPQ = 0;
bus_type=round(bus(:,10));
load_index=find(bus_type==3);
gen_index=find(bus_type==2);
SB=find(bus_type==1);
nSW=length(SB);
nPV=length(gen_index);
nPQ=length(load_index);
end
return
Solution by matlab code (Newton Rapson Method):
CONCLUSIONS:
On adding a shunt capacitance (which will produce 21MVar) compensates for
the power consumed by Generator at the slack Bus. A capacitor connected in
parallel produces reactive power thus, the reactive power demand of
generator at 1 is met by shunt capacitive compensation at bus 1. By adding
this reactive power of the generator becomes zero.