% 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