电力系统分析潮流计算大作业（源程序及实验报告）
源程序如下：
采用直角坐标系的牛顿-拉夫逊迭代
function chaoliujisuan()
m=3; %m=PQ节点个数
v=1;%v=PV节点个数
P=[-0.8055 -0.18 0]; %P=PQ节点的P值
Q=[-0.5320 -0.12 0]; %Q=PQ节点的Q值
PP=[0.5];%PP=PV节点的P值
V=[1.0];%V=PV节点的U值
E=[1 1 1 1.0 1.0]'; %E=PQ,PV,Vθ节点e的初值
F=[0 0 0 0 0]'; %F=PQ,PV,Vθ节点f的初值
G=[
6.3110 -3.5587 -2.7523 0 0;
-3.5587 8.5587 -5 0 0;
-2.7523 -5 7.7523 0 0;
0 0 0 0 0;
0 0 0 0 0
];
B=[
-20.4022 11.3879 9.1743 0 0;
11.3879 -31.00937 15 4.9889 0;
9.1743 15 -28.7757 0 4.9889;
0 4.9889 0 5.2493 0;
0 0 4.9889 0 -5.2493
];
Y=G+j*B;
X=[]; %X=△X
n=m+v+1;%总的节点数
FX=ones(2*n-2,1);%F(x)矩阵
F1=zeros(n-1,n-1);%F(x)导数矩阵
a=0;%记录迭代次数
EF=zeros(n-1,n-1);%最后的节点电压矩阵
while max(FX)>=10^(-5)
for i=1:m %PQ节点
FX(i)=P(i);%△P
FX(n+i-1)=Q(i);%△Q
for w=1:n
FX(i)=
FX(i)-E(i)*G(i,w)*E(w)+E(i)*B(i,w)*F(w)-F(i)*G(i,w)*F(w)-F(i)*B(i,w)*E(w); %△P
FX(n+i-1)=FX(n+i-1)-F(i)*G(i,w)*E(w)+F(i)*B(i,w)*F(w)+E(i)*G(i,w)*F(w)+E(i)*B(i ,w)*E(w); %△Q
end
end
for i=m+1:n-1 %PV节点
FX(i)=PP(i-m);%△P
FX(n+i-1)=V(i-m)^2-E(i)^2-F(i)^2; %△Q
for w=1:n
FX(i)=
FX(i)-E(i)*G(i,w)*E(w)+E(i)*B(i,w)*F(w)-F(i)*G(i,w)*F(w)-F(i)*B(i,w)*E(w);%△P end
end
for i=1:m %PQ节点
for w=1:n-1
if i~=w
F1(i,w)=-(G(i,w)*E(i)+B(i,w)*F(i));
F1(i,n+w-1)=B(i,w)*E(i)-G(i,w)*F(i);
F1(n+i-1,w)=B(i,w)*E(i)-G(i,w)*F(i);
F1(n+i-1,n+w-1)=G(i,w)*E(i)+B(i,w)*F(i);
else
F1(i,w)=-G(i,i)*E(i)-B(i,i)*F(i);
F1(i,n+w-1)=B(i,i)*E(i)-G(i,i)*F(i);
F1(n+i-1,w)=B(i,i)*E(i)-G(i,i)*F(i);
F1(n+i-1,n+w-1)=G(i,i)*E(i)+B(i,i)*F(i);
for k=1:n
F1(i,w)=F1(i,w)-G(i,k)*E(k)+B(i,k)*F(k);
F1(i,n+w-1)= F1(i,n+w-1)-G(i,k)*F(k)-B(i,k)*E(k);
F1(n+i-1,w)=F1(n+i-1,w)+G(i,k)*F(k)+B(i,k)*E(k);
F1(n+i-1,n+w-1)=F1(n+i-1,n+w-1)-G(i,k)*E(k)+B(i,k)*F(k);
end
end
end
end
for i=m+1:n-1 %PV节点
for w=1:n-1
if i~=w
F1(i,w)=-(G(i,w)*E(i)+B(i,w)*F(i));
F1(i,n+w-1)=B(i,w)*E(i)-G(i,w)*F(i);
F1(n+i-1,w)=0;
F1(n+i-1,n+w-1)=0;
else
F1(i,w)=-G(i,i)*E(i)-B(i,i)*F(i);
F1(i,n+w-1)=B(i,i)*E(i)-G(i,i)*F(i);
F1(n+i-1,w)=-2*E(i);
F1(n+i-1,n+w-1)=-2*F(i);
for k=1:n
F1(i,w)=F1(i,w)-G(i,k)*E(k)+B(i,k)*F(k);
F1(i,n+w-1)= F1(i,n+w-1)-G(i,k)*F(k)-B(i,k)*E(k);
end
end
end
end
X=inv(F1)*(-FX) ;
for i=1:n-1
E(i)=E(i)+X(i);
F(i)=F(i)+X(n+i-1);
end
a=a+1;
fprintf('第%d次迭代后的节点电压分别为:\n',a);
disp(E+j*F);
fprintf('第%d次迭代后功率偏差△P △Q电压偏差△V的平方分别为:\n',a); disp(FX);
end
disp('收敛后的节点电压用极坐标表示为:');
EF=E+j*F;
for i=1:n-1
fprintf('%d号节点电压的幅值为:',i)
disp(abs(EF(i)));
fprintf('%d号节点电压的相角度数为',i)
disp(angle(EF(i))*180/pi);
end
PPH=0;
for i=1:n
PPH=PPH+EF(n)*conj(Y(n,i))*conj(EF(i));
end
fprintf('平衡节点的功率');
disp(PPH);
运行结果：
