33节点配网重构.rar_33节点_IEEE 33_IEEE 配网_ieee_节点重构

function guan=fitness_cgfcPQVmulti(Swarm1) %加入不可行解判断程序（多个节点处加入分布式风电节点DG） b=32; n=33; LL=5; %联络开关数 Sb=10; %MW Vb=12.66; %KV Zb=Vb^2/Sb; %ohm check=1; checkhl=1; checkgd=1; jy=0; H=[7 6 5 4 3 2 20 19 18 33 0 0 0 0 0 0 0 0 0 0 0 14 13 12 11 10 9 34 0 0 0 0 0 0 0 0 0 0 0 0 0 0 11 10 9 8 7 6 5 4 3 2 21 20 19 18 35 0 0 0 0 0 0 17 16 15 14 13 12 11 10 9 8 7 6 32 31 30 29 28 27 26 25 36 24 23 22 28 27 26 25 5 4 3 37 0 0 0 0 0 0 0 0 0 0];%由函数matrixH生成 for i1=1:LL a(1,i1)=H(i1,Swarm1(1,i1)); end %a=[7 9 2 12 3] %狭长网络，电压越限的情况 %a=[7 14 9 31 37] %重构后的最优网络 %a=[33 34 35 36 37]; %原始网络 BranchM=[1 1 2 0.0922+i*0.047 %支路参数矩阵 2 2 3 0.4930+i*0.2511 3 3 4 0.3660+i*0.1864 4 4 5 0.3811+i*0.1941 5 5 6 0.8190+i*0.7070 6 6 7 0.1872+i*0.6188 7 7 8 0.7114+i*0.2351 8 8 9 1.0300+i*0.7400 9 9 10 1.0440+i*0.7400 10 10 11 0.1966+i*0.0650 11 11 12 0.3744+i*0.1238 12 12 13 1.4680+i*1.1550 13 13 14 0.5416+i*0.7129 14 14 15 0.5910+i*0.5260 15 15 16 0.7463+i*0.5450 16 16 17 1.2890+i*1.7210 17 17 18 0.3720+i*0.5740 18 2 19 0.1640+i*0.1565 19 19 20 1.5042+i*1.3554 20 20 21 0.4095+i*0.4784 21 21 22 0.7089+i*0.9373 22 3 23 0.4512+i*0.3083 23 23 24 0.8980+i*0.7091 24 24 25 0.8960+i*0.7011 25 6 26 0.2030+i*0.1034 26 26 27 0.2842+i*0.1447 27 27 28 1.0590+i*0.9337 28 28 29 0.8042+i*0.7006 29 29 30 0.5075+i*0.2585 30 30 31 0.9744+i*0.9630 31 31 32 0.3105+i*0.3619 32 32 33 0.3410+i*0.5362 33 8 21 2.0+i*2.0 34 9 15 2.0+i*2.0 35 12 22 2.0+i*2.0 36 18 33 0.5+i*0.5 37 25 29 0.5+i*0.5 ]; NodeM=[1 0 2 0.1000+i*0.0600 3 0.0900+i*0.0400 4 0.1200+i*0.0800 5 0.0600+i*0.0300 6 0.0600+i*0.0200 7 0.2000+i*0.1000 8 0.2000+i*0.1000 9 0.0600+i*0.0200 10 0.0600+i*0.0200 11 0.0450+i*0.0300 12 0.0600+i*0.0350 13 0.0600+i*0.0350 14 0.1200+i*0.0800 15 0.0600+i*0.0100 16 0.0600+i*0.0200 17 0.0600+i*0.0200 18 0.0900+i*0.0400 19 0.0900+i*0.0400 20 0.0900+i*0.0400 21 0.0900+i*0.0400 22 0.0900+i*0.0400 23 0.0900+i*0.0500 24 0.4200+i*0.2000 25 0.4200+i*0.2000 26 0.0600+i*0.0250 27 0.0600+i*0.0250 28 0.0600+i*0.0200 29 0.1200+i*0.0700 30 0.2000+i*0.6000 31 0.1500+i*0.0700 32 0.2100+i*0.1000 33 0.0600+i*0.0400 ]; %1、判断是否形成环路，F为支路环路关联矩阵(行表示回路，列表示断开开关，若任意两行相同，则表示形成了环路) F=zeros(5); for i1=1:LL %回路 for i2=1:LL %断开开关 if max(a(1,i2)==H(i1,:)) F(i1,i2)=1; end end end for i1=1:LL-1 for i2=i1+1:LL if min(F(:,i1)==F(:,i2)) checkhl=0; disp('出现环路')%出现环路时 guan=10000; end end end for i1=1:LL %按照断开开关矩阵，剔除Z矩阵中的断开支路 j=i1-1; for i2=1:b+LL-j if BranchM(i2,1)==a(1,i1) BranchM(i2,:)=[]; break end end end NodeN=zeros(n); %节点-节点关联矩阵A for i1=1:b NodeN(BranchM(i1,2),BranchM(i1,3))=1; NodeN(BranchM(i1,3),BranchM(i1,2))=1; end LayerM=[1]; %节点分层矩阵，电源节点号记“1” NU=zeros(1,n); %上层节点矩阵（有33列的行矩阵） while(checkhl==1) %以下用循环求取矩阵LayerM和NU while(checkgd==1) h=1; while(min(NU(2:33)~=0)==0) %NU矩阵的2-最后都有上层节点了，表示循环结束了 m=max(find(LayerM(:,h))); %m为矩阵LayerM第h列非零元素的个数 k=1; for i1=1:m g=LayerM(i1,h); %LayerM的第i1行第h列元素 ss=find(NodeN(g,:)==1); for i2=1:length(ss) if LayerM~=ss(1,i2) %排除相同节点 LayerM(k,h+1)=ss(1,i2); NU(1,ss(1,i2))=g; k=k+1; %k表示第h层含有的节点数 end end end h=h+1; %h表示网络分层的层数 if length(LayerM(1,:))==h-1 %如果网络分层矩阵没有搜索到下层节点，说明形成了断点，后边网络形成了孤岛，与电源节点没有连通回路 checkgd=0; disp('形成孤岛') guan=10000; check=0; break %结束循环 end end if min(NU(2:33)~=0) %若解可行，已经计算完LayerM，则让其跳出最外层while循环 checkgd=0; disp('可行解') end end while(check==1) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %添加两个风机pq（v）节点 %风力发电机数据 node1=18; %风机加入的节点号 node2=23; Pe=0.6; %单台风机的额定有功 NodeM(node1,2)=NodeM(node1,2)-Pe; %添加DG有功 NodeM(node2,2)=NodeM(node2,2)-Pe; BranchM(:,4)=BranchM(:,4)/Zb; %阻抗标幺化 NodeM(:,2)=NodeM(:,2)/Sb; %功率标幺化 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% x1=0.00453; %pu定子电抗 x2=0.149152; %pu转子电抗 x3=0.199904; %定子电抗与转子电抗之和 r2=0.00486; %pu转子电阻 xm=2.205952; %pu励磁电抗 Pe=0.6; %单台风机的额定有功 Qunit=0.04; %单组并联电容器容量 cosb=0.9; %设定风机的初始功率因数为0.9 Ue=0.69; %机端电压 %下面进行分层前推回代法潮流计算 V=ones(n,1); %节点电压 J=zeros(n,1); %支路电流 k=1; %记录迭代次数 V0=zeros(n,1); t=0; init1=NodeM(node1,2); init2=NodeM(node2,2); while(max(abs(V-V0))>1e-3) %判断收敛性（收敛精度设为1e-6） NodeM(node1,2)=init1; NodeM(node2,2)=init2; V0=V; %记录上一次迭代的电压值 %第一个风机节点 SN1=r2*((abs(V(node1,1))*Ue)^2-sqrt((abs(V(node1,1))*Ue)^4-4*x3^2*Pe^2))/(2*Pe*x3^2); Q1=([r2^2+x3*(xm+x3)*SN1^2]*Pe/(SN1*r2*xm))/Sb; cos1=Pe/sqrt(Pe^2+(Q1*Sb)^2); Qc1=Pe*(sqrt(1/(cos1^2)-1)-sqrt(1/(cosb^2)-1))/Sb; %并联电容器组无功补偿(nc为并联电容器的组数) nc1=fix(Qc1*Sb/Qunit); Qz1=nc1*Qunit*(abs(V(node1,1))^2)/Sb; cosc1=Pe/sqrt(Pe^2+((Q1-Qz1)*Sb)^2); %补偿后风机的功率因数（要求在0.9-1.0之间） %第二个风机节点 SN2=r2*((abs(V(node2,1))*Ue)^2-sqrt((abs(V(node2,1))*Ue)^4-4*x3^2*Pe^2))/(2*Pe*x3^2); Q2=([r2^2+x3*(xm+x3)*SN1^2]*Pe/(SN2*r2*xm))/Sb; cos2=Pe/sqrt(Pe^2+(Q2*Sb)^2); Qc2=Pe*(sqrt(1/(cos2^2)-1)-sqrt(1/(cosb^2)-1))/Sb; %并联电容器组无功补偿(nc为并联电容器的组数) nc2=fix(Qc2*Sb/Qunit); Qz2=nc2*Qunit*(abs(V(node2,1))^2)/Sb; cosc2=Pe/sqrt(Pe^2+((Q2-Qz2)*Sb)^2); %补偿后风机的功率因数（要求在0.9-1.0之间） %第一个接入风机的节点无功补偿 if cosc1<0.9 nc1=nc1+1; else nc1=nc1-1; end Qz1=nc1*Qunit*(abs(V(node1,1))^2)/Sb; cosc1=Pe/sqrt(Pe^2+((Q1-Qz1)*Sb)^2); while(cosc1<0.9) nc1=nc1+1; Qz1=nc1*Qunit*(abs(V(node1,1))^2)/Sb; cosc1=Pe/sqrt(Pe^2+((Q1-Qz1)*Sb)^2); if nc1>20 disp('网络变为狭长网络，末端功率因数过低，故终止进行无功补偿') break end end NodeM(node1,2)=NodeM(node1,2)+i*(Q1-Qz1); %第一个接入风机的节点无功补偿 if cosc2<0.9 nc2=nc2+1; else nc2=nc2-1; end Qz2=nc2*Qunit*(abs(V(node2,1))^2)/Sb; cosc2=Pe/sqrt(Pe^2+((Q2-Qz2)*Sb)^2); while(cosc2<0.9) nc2=nc2+1; Qz2=nc2*Qunit*(abs(V(node2,1))^2)/Sb; cosc2=Pe/sqrt(Pe^2+((Q2-Qz2)*Sb)^2); if nc2>20 disp('网络变为狭长网络，末端功率因数过低，故终止进行无功补偿') break end end NodeM(node2,2)=NodeM(node2,2)+i*(Q2-Qz2); %1、回代求支路电流矩阵J