%本程序的功能是用牛顿——拉夫逊法进行潮流计算
clear;
n=6;
nl=6;
isb=1;
pr=0.0001;
B1=[1 2 (0.065+0.208i)*27/484/2 0.00000364*27*484*2 1 0;
2 3 (0.065+0.208i)*17/484/2 0.00000364*17*484*2 1 0;
2 4 (0.065+0.208i)*6/484/2 0.00000364*6*484*2 1 0;
3 6 (0.065+0.208i)*18/484/2 0.00000364*18*484*2 1 0;
4 5 (0.065+0.208i)*5/484/2 0.00000364*5*484*2 1 0;
5 6 (0.065+0.208i)*24/484/2 0.00000364*24*484*2 1 0];
B2=[0 0 1.0455 1.05 0 1;
0 (250+120i)/100 1 1.05 0 2;
0 (330+190i)/100 1 1.05 0 2;
0 (250+80i)/100 1 1.05 0 2;
0 (400+200i)/100 1 1.05 0 2;
(1+0.6197i)*251.1/100 0 1.05 1.05 0 3];
%?B1 矩阵: 1、 支路首端号; 2、 末端号; 3、 支路阻抗; 4、 支路对地电纳
%?????????5、 支路的变比; 6、 支路首端处于 K 侧为 1, 1 侧为 0
%?B2 矩阵: 1、 该节点发电机功率; 2、 该节点负荷功率; 3、 节点电压初始值
%?????????4、 PV 节点电压 V 的给定值; 5、 节点所接的无功补偿设备的容量
%?????????6、 节点分类标号: 1 为平衡节点(应为 1 号节点); 2 为 PQ 节点;
%????????????3 为 PV 节点;
isn_1=input('是否在进行n‐1校验,是为1,否为0:');
if isn_1==0
isge=input('请输入发电机出力状态,满载为 1,轻载为 0:');
if isge==0,B2(6,1)=(1+0.6197i)*251.1/100;
else
if isge==1,B2(6,1)=(1+0.6197i)*753.3/100;
end
end
else
if isn_1==1,B2(6,1)=(1+0.6197i)*753.3/100;
nn_1=input('请输入n‐1校验断线序号,从 1‐6:');
B1(nn_1,3)=B1(nn_1,3)*2;
B1(nn_1,4)=B1(nn_1,4)/2;
end
end
Y=zeros(n);e=zeros(1,n);f=zeros(1,n);V=zeros(1,n);sida=zeros(1,n);S1=zeros(nl);
%?%?%‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐
for i=1:nl %支路数
if B1(i,6)==0 %左节点处于 1 侧
p=B1(i,1);q=B1(i,2);
else %左节点处于 K 侧
p=B1(i,2);q=B1(i,1);
end
Y(p,q)=Y(p,q)-1./(B1(i,3)*B1(i,5)); %非对角元
Y(q,p)=Y(p,q); %非对角元
Y(q,q)=Y(q,q)+1./(B1(i,3)*B1(i,5)^2)+B1(i,4)./2; %对角元 K 侧
Y(p,p)=Y(p,p)+1./B1(i,3)+B1(i,4)./2; %对角元 1 侧
end
%求导纳矩阵
disp('导纳矩阵?Y=');
disp(Y)
%‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐
G=real(Y);B=imag(Y); %分解出导纳阵的实部和虚部
for i=1:n %给定各节点初始电压的实部和虚部
e(i)=real(B2(i,3));
f(i)=imag(B2(i,3));
V(i)=B2(i