牛顿-拉夫逊法潮流计算.doc

%形成节点导纳矩阵%
G(1,1)=0.0157;B(1,1)=-0.1189;
G(1,2)=-0.0086;B(1,2)=0.0649;
G(1,3)=-0.0071;B(1,3)=0.0541;
G(1,4)=0;B(1,4)=0;
G(1,5)=0;B(1,5)=0;

G(2,1)=-0.0086;B(2,1)=0.0649;
G(2,2)=0.0111;B(2,2)=-0.0927;
G(2,3)=0;B(2,3)=0;
G(2,4)=0;B(2,4)=0;
G(2,5)=-0.0026;B(2,5)=0.0278;

G(3,1)=-0.0071;B(3,1)=0.0541;
G(3,2)=0;B(3,2)=0;
G(3,3)=0.0092;B(3,3)=-0.0763;
G(3,4)=0;B(3,4)=0;
G(3,5)=-0.0021;B(3,5)=0.0223;

G(4,1)=0;B(4,1)=0;
G(4,2)=0;B(4,2)=0;
G(4,3)=0;B(4,3)=0;
G(4,4)=0.0003;B(4,4)=-0.0127;
G(4,5)=-0.0003;B(4,5)=0.0127;

G(5,1)=0;B(5,1)=0;
G(5,2)=-0.0026;B(5,2)=0.0278;
G(5,3)=-0.0021;B(5,3)=0.0223;
G(5,4)=-0.0003;B(5,4)=0.0127
G(5,5)=0.0050;B(5,5)=-0.0628
Y=G+j*B;
%设定节点起始计算电压%
delt(1)=0;delt(2)=0;delt(3)=0;delt(4)=0;
u(1)=220;u(2)=220;u(3)=220;u(4)=220;
%在最大负荷下时的各节点注入功率%
p(1)=-60.55;q(1)=-39.29;p(2)=-150.83;q(2)=-101.46;
p(3)=-90.55;q(3)=-53.81;p(4)=117;q(4)=80;

%或者在最小负荷下时的各节点注入功率%
p(1)=-40.37;q(1)=-25.45;p(2)=-120.58;q(2)=-75.9;
p(3)=-60.38;q(3)=-34.08;p(4)=117;q(4)=80;
%置迭代次数k=0%
k=0;precision=1;
N1=4;(4个PQ节点)
%应用式（5.3）计算%
while precision>0.145
delt(5)=0;u(5)=220;
for m=1:N1
for n=1:N1+1

pt(n)=u(m)*u(n)*(G(m,n)*cos(delt(m)-delt(n))+B(m,n)*sin(delt(m)-delt(n)));
qt(n)=u(m)*u(n)*(G(m,n)*sin(delt(m)-delt(n))-B(m,n)*cos(delt(m)-delt(n)));
end
pp(m)=p(m)-sum(pt);qq(m)=q(m)-sum(qt);
end
%计算雅可比矩阵%
for m=1:N1
for n=1:N1+1
h0(n)=u(m)*u(n)*(G(m,n)*sin(delt(m)-delt(n))-B(m,n)*cos(delt(m)-delt(n)));
n0(n)=-u(m)*u(n)*(G(m,n)*cos(delt(m)-delt(n))+B(m,n)*sin(delt(m)-delt(n)));
j0(n)=-u(m)*u(n)*(G(m,n)*cos(delt(m)-delt(n))+B(m,n)*sin(delt(m)-delt(n)));
L0(n)=-u(m)*u(n)*(G(m,n)*sin(delt(m)-delt(n))-B(m,n)*cos(delt(m)-delt(n)));
end
H(m,m)=sum(h0)-u(m)^2*(G(m,m)*sin(delt(m)-delt(m))-B(m,m)*cos(delt(m)-delt(m)));
N(m,m)=sum(n0)-2*u(m)^2*G(m,m)+u(m)^2*(G(m,m)*cos(delt(m)-delt(m))+B(m,m)*sin(delt(m)-delt(m)));
J(m,m)=sum(j0)+u(m)^2*(G(m,m)*cos(delt(m)-delt(m))+B(m,m)*sin(delt(m)-delt(m)));
L(m,m)=sum(L0)+2*u(m)^2*B(m,m)+u(m)^2*(G(m,m)*sin(delt(m)-delt(m))-B(m,m)*cos(delt(m)-delt(m)));
End
%当m=n时,形成雅可比矩阵的对角形式,见公式5.5%
for m=1:N1
JJ(2*m-1,2*m-1)=H(m,m); JJ(2*m-1,2*m)=N(m,m);
JJ(2*m,2*m-1)=J(m,m); JJ(2*m,2*m)=L(m,m);
end
for m=1:N1
for n=1:N1
if m==n
else
%当m≠n时,形成雅可比矩阵的非对角形式,见公式5.6%
H(m,n)=-u(m)*u(n)*(G(m,n)*sin(delt(m)-delt(n))-B(m,n)*cos(delt(m)-delt(n)));
J(m,n)=u(m)*u(n)*(G(m,n)*cos(delt(m)-delt(n))+B(m,n)*sin(delt(m)-delt(n)));
N(m,n)=-J(m,n);L(m,n)=H(m,n);
JJ(2*m-1,2*n-1)=H(m,n);JJ(2*m-1,2*n)=N(m,n);
JJ(2*m,2*n-1)=J(m,n); JJ(2*m,2*n)=L(m,n);
end
end
end
%得到%
for m=1:N1
PP(2*m-1)=pp(m); PP(2*m)=qq(m);
end
%解修正方程式,计算电压修正量%
uu=-JJ\PP';precision=max(abs(uu));
%当,计算各节点电压新值%
for n=1:N1
delt(n)=delt(n)+uu(2*n-1);
u(n)=u(n)+uu(2*n);
%增大迭代次数%
k=k+1;
end
%输出每次迭代的电压新值,直到结束%
k-1,u,delt
end
%计算平衡节点的功率,并输出%
for n=1:N1+1
U(n)=u(n)*(cos(delt(n))+j*sin(delt(n)));
end
for m=1:N1+1
I(m)=Y(5,m)*U(m);
end
S5=U(5)*sum(conj(I))
%计算各线路的功率,并输出%
for m=1:N1+1
for n=1:N1+1
S(m,n)=U(m)*(conj(U(m))-conj(U(n)))*conj(-Y(m,n));
end
end
S