配网重构中，很重要的一个约束条件为配网应随时保持开环、辐射的状态： 配电网系统是属于闭环设计但是开环运行的系统，因此，在开关的开闭过程中，随时保持配电网的开环状态时很重要。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);
end2.智能算法重构代码
如粒子群算法中，变量可以选取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
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.7320+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 %要保证联络开关的首节点必须是T=3节点，编程决定的。（如果两个节点T=2时就不做要求了）
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 %节点参数矩阵(电源节点负荷为0)
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
% 10 0.4200+i*0.2000
11 0.0450+i*0.0300
12 0.0600+i*0.0350
13 0.0600+i*0.0350
14 0.1200+i*0.0800
% 14 0.4200+i*0.2000
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 ];
% NodeM(:,2)=NodeM(:,2)+SDG;
%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
https://download.csdn.net/download/zhangxd212489/73526211