clear;close;clc;
Branch=[1 2 0.0192+0.0575j 0.0264j 1 0
1 3 0.0452+0.1852j 0.0204j 1 0
2 4 0.057+0.1737j 0.0184j 1 0
3 4 0.0132+0.0379j 0.0042j 1 0
2 5 0.0472+0.1983j 0.0209j 1 0
2 6 0.0581+0.1763j 0.0187j 1 0
4 6 0.0119+0.0414j 0.0045j 1 0
5 7 0.046+0.116j 0.0102j 1 0
6 7 0.0267+0.082j 0.0085j 1 0
6 8 0.012+0.042j 0.0045j 1 0
6 9 0.208j 0 1.0155 1
6 10 0.556j 0 0.9629 1
9 11 0.208j 0 1 0
9 10 0.11j 0 1 0
4 12 0.256j 0 1.0129 1
12 13 0.14j 0 1 0
12 14 0.1231+0.2559j 0 1 0
12 15 0.0662+0.1304j 0 1 0
12 16 0.0945+0.1987j 0 1 0
14 15 0.221+0.1997j 0 1 0
16 17 0.0824+0.1932j 0 1 0 %第三问中断开即去掉这一行;第三问中采用双回路,则将支路对地阻抗变为之前的1/2
15 18 0.107+0.2185j 0 1 0
18 19 0.0639+0.1292j 0 1 0
19 20 0.034+0.068j 0 1 0
10 20 0.0936+0.209j 0 1 0
10 17 0.0324+0.0845j 0 1 0
10 21 0.0348+0.0749j 0 1 0
10 22 0.0727+0.1799j 0 1 0
21 22 0.0116+0.0236j 0 1 0
15 23 0.1+0.202j 0 1 0
22 24 0.115+0.179j 0 1 0
23 24 0.132+0.27j 0 1 0
24 25 0.1885+0.3292j 0 1 0
25 26 0.2554+0.38j 0 1 0
25 27 0.1093+0.2087 0 1 0
28 27 0.396j 0 0.9581 1
27 29 0.2198+0.4153j 0 1 0
27 30 0.3202+0.6027j 0 1 0
29 30 0.2399+0.4533j 0 1 0
8 28 0.0636+0.2j 0.0214j 1 0
6 28 0.0169+0.0599j 0.0065j 1 0
]
%Branch矩阵:1、支路首端号;2、支路末端号;3、支路阻抗;4、支路对地导纳;5、变压器变比k;6、支路首端处于K侧为1,1侧为0(是否含有变压器)
Y=zeros(30); %节点导纳矩阵
for i=1:41 %因为有41条支路(第三问中断开16、17之间支路后改为40)
if Branch(i,6)==0 %不含变压器的支路
j=Branch(i,1);
k=Branch(i,2);
Y(j,k)=Y(j,k)-1/Branch(i,3);%互导纳
Y(k,j)=Y(j,k);
Y(j,j)=Y(j,j)+1/Branch(i,3)+Branch(i,4);%自导纳
Y(k,k)=Y(k,k)+1/Branch(i,3)+Branch(i,4);
else %含变压器的支路
j=Branch(i,1);
k=Branch(i,2);
Y(j,k)=Y(j,k)-Branch(i,5)/Branch(i,3);
Y(k,j)=Y(j,k);
Y(j,j)=Y(j,j)+1/Branch(i,3);
Y(k,k)=Y(k,k)+Branch(i,5)^2/Branch(i,3);
end
end %计算每个节点的自导纳和互导纳
disp('节点导纳矩阵');
Y
G=real(Y);B=imag(Y); %G为节点导纳矩阵的实部,B为其虚部
V=[1.05;1.0338;1;1;1.0058;
1;1;1.023;1;1;
1.0913;1;1.0883;1;1;
1;1;1;1;1;
1;1;1;1;1;
1;1;1;1;1]; %给定节点电压的初始计算值
e=real(V);f=imag(V); %e为节点电压实部,f为其虚部
disp('节点电压的实部:')
e
disp('节点电压的虚部:')
f
disp('节点注入有功功率:')
Ps=[0;0.3586;-0.024;-0.076;-0.6964;
0;-0.228;0.05;0;-0.058;
0.1793;-0.112;0.1691;-0.062;-0.082;
-0.035;-0.09;-0.032;-0.095;-0.022; %增加20%后为-0.0308
-0.175;0;-0.032;-0.087;0;
-0.035;0;0;-0.024;-0.106]
disp('节点注入无功功率:')
Qs=[0;-0.127;-0.012;-0.016;-0.19;
0;-0.109;-0.3;0;-0.02;
0;-0.075;0;-0.016;-0.025;
-0.018;-0.058;-0.009;-0.034;-0.007; %增加30%前为-0.0091
-0.112;0;-0.016;-0.067;0;
-0.023;0;0;-0.009;-0.019]
%由各节点电压向量(状态变量)可得各节点注入功率:
P=e.*(G*e-B*f)+f.*(G*f+B*e);
Q=f.*(G*e-B*f)-e.*(G*f+B*e);
del_W=[
0.3586-P(2);1.0338^2-e(2)^2-f(2)^2; %1节点为平衡节点故不参与计算
-0.024-P(3);-0.012-Q(3);
-0.076-P(4);-0.012-Q(4);
-0.6964-P(5);1.0058^2-e(5)^2-f(5)^2;
0-P(6);0-Q(6);
-0.228-P(7);-0.109-Q(7);
0.05-P(8);1.023^2-e(8)^2-f(8)^2;
0-P(9);0-Q(9);
-0.058-P(10);-0.02-Q(10);
0.1793-P(11);1.0913^2-e(11)^2-f(11)^2;
-0.112-P(12);-0.075-Q(12);
0.1691-P(13);1.0883^2-e(13)^2-f(13)^2;
-0.062-P(14);-0.016-Q(14);
-0.082-P(15);-0.025-Q(15);
-0.035-P(16);-0.018-Q(16);
-0.09-P(17);-0.058-Q(17);
-0.032-P(18);-0.009-Q(18);
-0.095-P(19);-0.034-Q(19);
-0.022-P(20);-0.007-Q(20);
-0.175-P(21);-0.112-Q(21);
0-P(22);0-Q(22);
-0.032-P(23);-0.016-Q(23);
-0.087-P(24);-0.067-Q(24);
0-P(25);0-Q(25);
-0.035-P(26);-0.023-Q(26);
0-P(27);0-Q(27);
0-P(28);0-Q(28);
-0.024-P(29);-0.009-Q(29);
-0.106-P(30);-0.019-Q(30);
] %del_W是失配功率
n=0;%记录迭代次数
while (any(del_W>1e-5))&(n<5) %循环结束的条件:失配功率<1e-5或迭代达到5次
n=n+1;
disp('迭代次数:');disp(n)
%--------------------------------------------------------
%修正方程式:
% [0.3586-P(2) ] [del_e(2)]
% [1.0338^2-e(2)^2-f(2)^2 ] [del_f(2)]
% ?????? + J * ??????
% [-0.106-P(30) ] [del_e(30)]
% [-0.019-Q(30) ] [del_f(30)]
%雅克比矩阵为:J=[
% P1_e1 P1_f1 P1_e2 P1_f2 P1_e3 P1_e3 ???
% V1_e1 V1_f1 V1_e2 V1_f2 V1_e3 V1_f3 ???
% ??? ???
% P30_e1 P30_f1 P30_e2 P30_f2 P30_e3 P30_f3 ???
% Q30_e1 Q30_f1 Q30_e2 Q30_f2 Q30_e3 Q30_f3 ???
% ]
%--------------------------------------------------------
%-------------------求雅可比矩阵的参数-----------------------%
GeBf=G*e-B*f;%辅助计算函数
GfBe=G*f+B*e;%辅助计算函数
for i=3:4 %PQ节点3、4
for j=2:30 %计算时注意平衡节点不参与
if i==j %i=j和i!=j的计算公式不同
P_e(i,j)=-GeBf(i)-G(i,j)*e(i)-B(i,j)*f(i);
P_f(i,j)=-GfBe(i)-G(i,j)*f(i)+B(i,j)*e(i);
Q_e(i,j)=GfBe(i)+B(i,j)*e(i)-G(i,j)*f(i);
Q_f(i,j)=-GeBf(i)+G(i,j)*e(i)+B(i,j)*f(i);
else
P_e(i,j)=-G(i,j)*e(i)-B(i,j)*f(i);
Q_f(i,j)=-P_e(i,j);
P_f(i,j)=B(i,j)*e(i)-G(i,j)*f(i);
Q_e(i,j)=P_f(i,j);
end
end
end
for i=6:7 %PQ节点6、7
for j=2:30
if i==j
P_e(i,j)=-GeBf(i)-G(i,j)*e(i)-B(i,j)*f(i);
P_f(i,j)=-GfBe(i)-G(i,j)*f(i)+B(i,j)*e(i);
Q_e(i,j)=GfBe(i)+B(i,j)*e(i)-G(i,j)*f(i);
Q_f(i,j)=-GeBf(i)+G(i,j)*e(i)+B(i,j)*f(i);
else
P_e(i,j)=-G(i,j)*e(i)-B(i,j)*f(i);
Q_f(i,j)=-P_e(i,j);
P_f(i,j)=B(i,j)*e(i)-G(i,j)*f(i);
Q_e(i,j)=P_f(i,j);
end
end
end
for i=9:10 %PQ节点9、10
for j=2:30
if i==j
P_e(i,j)=-GeBf(i)-G(i,j)*e(i)-B(i,j)*f(i);
P_f(i,j)=-GfBe(i)-G(i,j)*f(i)+B(i,j)*e(i);
Q_e(i,j)=GfBe(i)+B(i,j)*e(i)-G(i,j)*f(i);
Q_f(i,j)=-GeBf(i)+G(i,j)*e(i)+B(i,j)*f(i);
else
P_e(i,j)=-G(i,j)*e(i)-B(i,j)*f(i);
Q_f(i,j)=-P_e(i,j);
P_f(i,j)=B(i,j)*e(i)-G(i,j)*f(i);
Q_e(i,j)=P_f(i,j);
end
end
end
for i=12:12 %PQ节点12
for j=2:30
if j==i
P_e(i,j)=-GeBf(i)-G(i,j)*e(i)-B(i,j)*f(i);
P_f(i,j)=-GfBe(i)-G(i,j)*f(i)+B(i,j)*e(i);
Q_e(i,j)=GfBe(i)+B(i,j)*e(i)-G(i,j)*f(i);
Q_f(i,j)=-GeBf(i)+G(i,j)*e(i)+B(i,j)*f(i);
else
P_e(i,j)=-G(i,j)*e(i)-B(i,j)*f(i);
Q_f(i,j)=-P_e(i,j);
P_f(i,j)=B(i,j)*e(i)-G(i,j)*f(i);
Q_e(i,j)=P_f(i,j);
end
end
end
for i=14:30 %PQ节点14~30
for j=2:30
if i==j
P_e(i,j)=-GeBf(i)-G(i,j)*e(i)-B(i,j)*f(i);
P_f(i,j)=-GfBe(i)-G(i,j)*f(i)+B(i,j)*e(i);
Q_e(i,j)=GfBe(i)+B(i,j)*e(i)-G(i,j)*f(i);
Q_f(i,j)=-GeBf(i)+G(i,j)*e(i)+B(i,j)*f(i);
else
P_e(i,j)=-G(i,j)*e(i)-B(i,j)*f(i);
Q_f(i,j)=-P_e(i,j);
P_f(i,j)=B(i,j)*e(i)-G(i,j)*f(i);
Q_e(i,j)=P_f(i,j);
end
end
end
%以下是PV节点求解,由于PV节点独立,故独立求解
for j=2:30
if j==2 %PV节点2
P_e(j,j)=-GeBf(j)-G(j,j)*e(j)-B(j,