网重构中，很重要的一个约束条件为配网应随时保持开环、辐射的状态： 配电网系统是属于闭环设计但是开环运行的系统，因此，在开关的开闭过程中，随时保持配电网的开环状态时很重要。Mendoza等利用图论，尤其是基本环向量的方法，不仅仅有效地保持了网络的辐射性，还大大地节省了搜寻时间和空间。还有文献对如何利用基本环向量去创造可行解的规则进行了详细阐述，文 中定义了基本节点，外节点，内节点，环向量，基本支路向量以及禁止的组向量 避免不可行解。本文中，应用了一个更加简单但同样有效的基于基本环向量的 方法保持配电网的开环，下面将进行详细阐述：
一、理论分析
本文中算例采用的是IEEE 33节点系统，因此以IEEE 33节点系统为例说明，配网中保持配网辐射状态的方法，IEEE 33节点系统有33个节点，32个分段开关和5个联络开关，为了使得网络保持辐射运行状态，具体的步骤如下：
1．将系统中的各个支路进行标号，0代表支路上的开关是断开的，1代表支 路上的开关是闭合的； 2．利用图论决定基本的环路，环路的个数与联络开关的个数是一样的，对33 节点系统而言即5个，如下图所示。
3.然后，根据基本的环路列出基本环矩阵M。 M的行数即基本环的个数， 列数为支路数最多的那个环的支路数，其他小一些的环中没有这么多支路的空余 的地方值都填充为0。本文中，为了更简便，每条支路，包括环共有的支路，都只会包含在一个支路环里，这样就省略了节点分类以及还需创建其他支路的麻烦。 编码的规则是从图的最上面的支路编码到最下面的支路，从最左边的支路编码到最右边的支路。公共的支路只会被包含在首次将其编码进去的环内。通过上面的 编码即得到IEEE 33节点的基本环矩阵M：
4．从矩阵M的每一行中选取一个非零元素，在每一次的网络重构中，每个环路中的都有且只有一条支路(非零元素)能被选为断开，每个环中被选中的断开的支路开关状态改为0，其余支路的开关状态都保持为1。
通过上面的步骤，在重构过程中没有非可行解的产生，因此整个重构过程的步骤大大简化。
二、程序介绍
上述理论非常适合应用于智能算法的配电网重构中，首先需要通过配电网参数形成M矩阵，具体代码如下：
1.基本环矩阵M的matlab代码
%产生回路编号矩阵M b=32; %支路数，开断开关数 n=33; %节点数 LL=5; %联络开关数 Sb=10; %MW Vb=12.66; %KV Zb=Vb^2/Sb; %ohm %第一列存支路号，第二列存首节点号，第三列存尾节点号，第四列存支路自阻抗(ohm),第五列存尾节点给定功率(MW\Mvar) %其中：33-37支路为5个联络开关 Z=[1 0 1 0.0922+i*0.047 0.1000+i*0.0600 2 1 2 0.4930+i*0.2511 0.0900+i*0.0400 3 2 3 0.3660+i*0.1864 0.1200+i*0.0800 4 3 4 0.3811+i*0.1941 0.0600+i*0.0300 5 4 5 0.8190+i*0.7070 0.0600+i*0.0200 6 5 6 0.1872+i*0.6188 0.2000+i*0.1000 7 6 7 0.7114+i*0.2351 0.2000+i*0.1000 8 7 8 1.0300+i*0.7400 0.0600+i*0.0200 9 8 9 1.0440+i*0.7400 0.0600+i*0.0200 10 9 10 0.1966+i*0.0650 0.0450+i*0.0300 11 10 11 0.3744+i*0.1238 0.0600+i*0.0350 12 11 12 1.4680+i*1.1550 0.0600+i*0.0350 13 12 13 0.5416+i*0.7129 0.1200+i*0.0800 14 13 14 0.5910+i*0.5260 0.0600+i*0.0100 15 14 15 0.7463+i*0.5450 0.0600+i*0.0200 16 15 16 1.2890+i*1.7210 0.0600+i*0.0200 17 16 17 0.3720+i*0.5740 0.0900+i*0.0400 18 1 18 0.1640+i*0.1565 0.0900+i*0.0400 19 18 19 1.5042+i*1.3554 0.0900+i*0.0400 20 19 20 0.4095+i*0.4784 0.0900+i*0.0400 21 20 21 0.7089+i*0.9373 0.0900+i*0.0400 22 2 22 0.4512+i*0.3083 0.0900+i*0.0500 23 22 23 0.8980+i*0.7091 0.4200+i*0.2000 24 23 24 0.8960+i*0.7011 0.4200+i*0.2000 25 5 25 0.2030+i*0.1034 0.0600+i*0.0250 26 25 26 0.2842+i*0.1447 0.0600+i*0.0250 27 26 27 1.0590+i*0.9337 0.0600+i*0.0200 28 27 28 0.8042+i*0.7006 0.1200+i*0.0700 29 28 29 0.5075+i*0.2585 0.2000+i*0.6000 30 29 30 0.9744+i*0.9630 0.1500+i*0.0700 31 30 31 0.3105+i*0.3619 0.2100+i*0.1000 32 31 32 0.3410+i*0.5362 0.0600+i*0.0400]; C=[33 7 20 2.0+i*2.0 0 %联络开关矩阵C 34 8 14 2.0+i*2.0 0 35 11 21 2.0+i*2.0 0 36 17 32 0.5+i*0.5 0 37 24 28 0.5+i*0.5 0 ]; % 定义节点到母节点路径矩阵A A=zeros(b); %对矩阵A进行初始化 for i=1:b A(i,1)=i; %将路径矩阵A的第一列储存为支路号 end for j=1:b %利用Z矩阵中的首末节点列向量进行路径搜索 n1=j; k=2; while(j~=0) %一直搜寻到母节点0结束 for m=1:b if(Z(m,3)==j) %在末节点定位j A(n1,k)=Z(m,2); j=Z(m,2); %将首节点赋给j,然后再在末节点定位j k=k+1; end end end end %遍历A矩阵形成单回路矩阵H for i=1:LL jd=max(intersect(A(C(i,2),:),A(C(i,3),:))); %取回路交叉点jd,保证是单回路 h=1; while(A(C(i,2),h)~=jd) M(i,h)=A(C(i,2),h); h=h+1; end h1=1; while(A(C(i,3),h1)~=jd) M(i,h)=A(C(i,3),h1); h=h+1; h1=h1+1; end M(i,h)=C(i,1); end
2.智能算法重构代码
如粒子群算法中，变量可以选取M矩阵每行中的非零元素，通过迭代得到最佳的重构后网络，重构子函数如下所示。
%该为智能算法程序计算子函数，加入不可行解判断程序（在指定节点处加入DG） b=32; k=1; n=33; LL=5; %联络开关数 z Sb=10; %MW Vb=12.66; %KV Zb=Vb^2/Sb; %ohm check=1; checkhl=1; checkgd=1; H=[ 2 3 4 5 6 7 33 20 19 18 0 0 0 0 0 0 0 0 0 0 0 9 10 11 12 13 14 34 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 3 4 5 6 7 8 9 10 11 35 21 20 19 18 0 0 0 0 0 0 6 7 8 9 10 11 12 13 14 15 16 17 36 32 31 30 29 28 27 26 25 3 4 5 25 26 27 28 37 24 23 22 0 0 0 0 0 0 0 0 0 0];%由函数matrixH生成 for i1=1:LL a(1,i1)=H(i1,Swarm1(1,i1)); end %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) BranchM(:,4)=BranchM(:,4)/Zb; %阻抗标幺化 NodeM(:,2)=NodeM(:,2)/Sb; %功率标幺化 %下面进行分层前推回代法潮流计算 V=ones(n,1); %节点电压 J=zeros(n,1); %支路电流 k=1; %记录迭代次数 V0=zeros(n,1); t=0; while(max(abs(V-V0))>1e-3) %判断收敛性（收敛精度设为1e-6） V0=V; %记录上一次迭代的电压值 %1、回代求支路电流矩阵J I=conj(NodeM(:,2)./V); %节点注入电流 for i1=h:-1:2 ss=find(LayerM(:,i1)~=0); for i2=1:length(ss) if min(LayerM(i2,i1)==NU) J(LayerM(i2,i1))=I(LayerM(i2,i1)); else sumJ=0; gx=find(NU==LayerM(i2,i1)); %求得节点LayerM(i2,i1)的下层节点矩阵（大小不定） for i3=1:length(gx) sumJ=sumJ+J(gx(1,i3)); end J(LayerM(i2,i1))=I(LayerM(i2,i1))+sumJ; end end end %2、前推求节点电压矩阵V Z=zeros(n,1); %支路阻抗 for i1=2:h ss=find(LayerM(:,i1)~=0); for i2=1:length(ss) m1=NU(1,LayerM(i2,i1)); %首节点 n1=LayerM(i2,i1); %尾节点 for i3=1:b %从BranchM中搜寻得到各支路阻抗矩阵 if (BranchM(i3,2)==m1&&BranchM(i3,3)==n1)||(BranchM(i3,2)==n1&&BranchM(i3,3)==m1) Z(n1,1)=BranchM(i3,4); end end V(n1,1)=V(m1,1)-Z(n1,1)*J(n1,1); end end k=k+1; %迭代次数k if k>20 guan=10000; break end end %% Vby=abs(V); V=V*Vb; %反标幺 Vm=abs(V); %节点电压幅值 %———————————————————————————————— %节点电压上下限约束（排除a=[7 9 2 12 3];这种状况） %if min(Vm)<Vb*0.9||min(Vm)>Vb %guan=10000; %end %———————————————————————————————— Va=angle(V); %节点电压弧度 %下面计算有功网损fploss Sloss=zeros(n,1); Z=zeros(n,1); %支路阻抗 for i1=2:h %！！！此for循环程序段是为了搜寻网络各支路的首末节点与支路阻抗，应对网络拓扑变化，造成的潮流流向变化 ss=find(LayerM(:,i1)~=0); for i2=1:length(ss) m1=NU(1,LayerM(i2,i1)); %首节点 n1=LayerM(i2,i1); %尾节点 for i3=1:b %从BranchM中搜寻得到各支路阻抗矩阵 if (BranchM(i3,2)==m1&&BranchM(i3,3)==n1)||(BranchM(i3,2)==n1&&BranchM(i3,3)==m1) Z(n1,1)=BranchM(i3,4); end end Sloss(n1)=(abs(V(m1)-V(n1)))^2/conj(Z(n1,1)*Zb)*1000; %各支路的功率损耗(单位为kW) end end Ploss=real(Sloss); %有功损耗 Qloss=imag(Sloss); %无功损耗 fPloss=sum(Ploss); %系统的总有功损耗 guan=fPloss; check=0; end checkhl=0; end if k>20 guan=10000; end
完整粒子群算法的重构程序代码链接：
