牛顿法潮流计算Matlab程序

Temp1=Temp1+Um(k)*(G(L,k)*cos(delta)+B(L,k)*sin(delta));
Temp2=Temp2+Um(k)*(G(L,k)*sin(delta)-B(L,k)*cos(delta));
end
P(L)=Um(L)*Temp1;
Q(L)=Um(L)*Temp2;
Correct=Jacobi\deltaPQ;
% Part 3.4: Correcting voltage for PQ nodes and angles for PV nodes;
5-j*15 1.667-j*5 0];
Y=[y(1,2)+y(1,3) -y(1,2) -y(1,3);
-y(1,2) y(1,2)+y(2,3) -y(2,3);
-y(1,3) -y(2,3) y(1,3)+y(2,3)];
G=tive matrix;
Temp1=G(2,3)*cos(delta)+B(2,3)*sin(delta);
N23=-Um(2)*Um(3)*Temp1;
N33=-P(3)-G(3,3)*Um(3)*Um(3);
% J matrix;
J32=-N23;
J33=-P(3)+G(3,3)*Um(3)^2;
% L matrix;
% This program is to calculate power flow using Newton-Raphonson iterative
% method. Developed by Wang shishan, 2010-10-27.
for k=2:3
sita(k)=sita(k)-Correct(k-1);
end
Um(3)=Um(3)-Correct(3)*Um(3);
% Part 3.5: Checking Q2?
Q2_G=SB*Q(2)-Qout(2);
if Q2_G>Q2_range(2)
Q(2)=Q2_range(2)/SB;
elseif Q2_G<Q2_range(1)
Q(2)=Q2_range(1)/SB;
else
disp('Reactive power at 2th node in a norm range');
end
end % END-iteration.
% Part 4: Postprocess----Power flowing in every branch;
B=imag(Y); % Image part, susceptance matrix;
% Part 2: Initialing voltage, power value;
SB=100 ; % MVA, base power;
Pin=[0;20;0]; % Input active power at 2th node;
for k=(L+1):3
delta=sita(L)-sita(k);
Temp1=G(L,k)*sin(delta)-B(L,k)*cos(delta);
H(L,k)=-Um(L)*Um(k)*Temp1;
H(k,L)=H(L,k);
end
end
% N matrix;
delta=sita(2)-sita(3);
% Part 3.1: Solving error vectors of active and reactive powers;
for L=2:3 % L: No. of node;
Temp1=0;
Temp2=0;
for k=1:3
delta=sita(L)-sita(k);
L33=-Q(3)+B(3,3)*Um(3)^2;
% Creating Jacobi matrix;
Jacobi=[H(2,2) H(2,3) N23;
H(3,2) H(3,3) N33;
J32 J33 L33]
% Part 3.3: Solving the correcting equation;
Q2_range=[0 35]; % Range of min and max for Q2;
U=[1.05;1.03;1.0];
sita=angle(U); % Phase angle of voltage;
Um=abs(U); % Amplititude of voltage;
Nmax=10;
TL=1E-06; % Tolerance level;
% Part 3: Iterative process;
for Ite=1:Nmax % Iterative vary;
end
deltaP=(Pin-Pout)/SB-P';
deltaQ=(Qin-Qout)/SB-Q';
deltaPQ=[deltaP(2) deltaP(3) deltaQ(3)]';
if max(deltaPQ)<TL
disp('Iteration is sucessful.')
for k=1:Nn
U(k)=Um(k)*exp(j*sita(k));
end
for I=1:Nn
for J=1:Nn
S(I,J)=U(I)*conj((U(I)-U(J))*y(I,J));
end
end
%THE END,
Ite
Um
sita
break;
end
% Part 3.2: Creating Jacobi elements coreponding matrix J;
% H matrix;
for L=2:3
H(L,L)=Q(L)+B(L,L)*Um(L)*Um(L);
% Part 1: Creating admittance matrix Y;
clear;
Nn=3; % Number of nodes;
y=[0 1.25-j*3.75 5-j*15
1.25-j*3.75 0 1.667-j*5
Pout=[0;50;60]; % Output active power at 2th node;
Qin=[0;0;0]; % Inpouring reactive power;
Qout=[0; -20;25]; % Outing reactive power;