% 6.061/6.690 Spring Term 2011
% Problem 7.1: Simple-Minded Load Flow Example
% First, impedances
Z1=.005+j*.1;
Z2=.01+j*.1;
Z3=.005+j*.15;
Z4=.001+j*.05;
Z5=.005+j*.1;
Z6=.005+j*.2;
Z7=.01+j*.3;
Z8=.005+j*.05;
%
% Enter the node-incidence Matrix Here: It should have as
% many rows as there are buses and as many columns as
% there are lines
NI=[];
% This is the vector of "known" voltage magnitudes
VNM = [0 0 1 1 0 0]';
% And the vector of known voltage angles
VNA = [0 0 0 0 0 0]';
% and this is the "key" to which are actually known
KNM = [0 0 1 1 0 0]';
KNA = [0 0 0 1 0 0]';
% and which are to be manipulated by the system
KUM = 1 - KNM;
KUA = 1 - KNA;
% Here are the known loads (positive is INTO network
% Use zeros for unknowns
P=[2 -2 1 0 -1 0]';
Q=[0 -.5 0 0 -.5 0]';
% and here are the corresponding vectors to indicate
% which elements should be checked in error checking
PC = [1 1 1 0 1 1]';
QC = [1 1 0 0 1 1]';
Check = KNM + KNA + PC + QC;
% Unknown P and Q vectors
PU = 1 - PC;
QU = 1 - QC;
fprintf('Here is the line admittance matrix:\n');
Y=[1/Z1 0 0 0 0 0 0 0;
0 1/Z2 0 0 0 0 0 0;
0 0 1/Z3 0 0 0 0 0;
0 0 0 1/Z4 0 0 0 0;
0 0 0 0 1/Z5 0 0 0;
0 0 0 0 0 1/Z6 0 0;
0 0 0 0 0 0 1/Z7 0;
0 0 0 0 0 0 0 1/Z8]
% Construct Node-Admittance Matrix
fprintf('And here is the bus admittance matrix\n')
YN=NI*Y*NI'
% Now: here are some starting voltage magnitudes and angles
VM = [1 1 1 1 1 1]';
VA = [.0965 .146 .00713 .0261 0 0]';
% Here starts a loop
Error = 1;
Tol=1e-16;
N = length(VNM);
% Construct a candidate voltage from what we have so far
VMAG = VNM .* KNM + VM .* KUM;
VANG = VNA .* KNA + VA .* KUA;
V = VMAG .* exp(j .* VANG);
% and calculate power to start
I = (YN*V);
PI = real(V .* conj(I));
QI = imag(V .* conj(I));
%pause
while(Error>Tol);
for i=1:N, % Run through all of the buses
% What we do depends on what bus!
if (KUM(i) == 1) & (KUA(i) == 1), % don't know voltage magnitude or angle
pvc= (P(i)-j*Q(i))/conj(V(i));
for n=1:N,
if n ~=i, pvc = pvc - (YN(i,n) * V(n)); end
end
V(i) = pvc/YN(i,i);
elseif (KUM(i) == 0) & (KUA(i) == 1), % know magnitude but not angle
% first must generate an estimate for Q
Qn = imag(V(i) * conj(YN(i,:)*V));
pvc= (P(i)-j*Qn)/conj(V(i));
for n=1:N,
if n ~=i, pvc = pvc - (YN(i,n) * V(n)); end
end
pv=pvc/YN(i,i);
V(i) = VNM(i) * exp(j*angle(pv));
end % probably should have more cases
end % one shot through voltage list: check error
% Now calculate currents indicated by this voltage expression
I = (YN*V);
% For error checking purposes, compute indicated power
PI = real(V .* conj(I));
QI = imag(V .* conj(I));
% Now we find out how close we are to desired conditions
PERR = (P-PI) .* PC;
QERR = (Q-QI) .* QC;
Error = sum(abs(PERR) .^2 + abs(QERR) .^2);
end
fprintf('Here are the voltages\n')
V
fprintf('Magnitudes\n')
abs(V)
fprintf('Real Power\n')
PI
fprintf('Reactive Power\n')
QI
fprintf('Line Voltages are\n')
Vline = NI'*V
fprintf('Line Currents are\n')
Iline = Y*Vline