前言：之前分享的采用智能算法的配电网优化重构和故障重构程序，主要是采用罚函数的方式确保网络不出现环网和孤岛，但是在多目标pareto求解中，采用罚函数的方式大大降低了程序可用性，因此主要的难点就在于此，本次和大家看一下多目标重构的编程方法。
一、主要理论
以IEEE 33节点系统为例说明，配网中保持配网辐射状态的方法，IEEE 33节点系统有33个节点，32个分段开关和5个联络开关，为了使得网络保持辐射运行状态，具体的步骤如下：
1．将系统中的各个支路进行标号，0代表支路上的开关是断开的，1代表支 路上的开关是闭合的； 2．利用图论决定基本的环路，环路的个数与联络开关的个数是一样的，对33 节点系统而言即5个，如下图所示。
3.然后，根据基本的环路列出基本环矩阵M。 M的行数即基本环的个数， 列数为支路数最多的那个环的支路数，其他小一些的环中没有这么多支路的空余 的地方值都填充为0。本文中，为了更简便，每条支路，包括环共有的支路，都只会包含在一个支路环里，这样就省略了节点分类以及还需创建其他支路的麻烦。 编码的规则是从图的最上面的支路编码到最下面的支路，从最左边的支路编码到最右边的支路。公共的支路只会被包含在首次将其编码进去的环内。通过上面的编码即得到IEEE 33节点的基本环矩阵M：
4．从矩阵M的每一行中选取一个非零元素，在每一次的网络重构中，每个环路中的都有且只有一条支路(非零元素)能被选为断开，每个环中被选中的断开的支路开关状态改为0，其余支路的开关状态都保持为1。
通过上面的步骤，在重构过程中没有非可行解的产生，因此整个重构过程的步骤大大简化。
该算例的目标函数分别是网损、电压偏差和开关开断成本。
二、主要程序内容
通过引入kxx函数改变传统罚函数方式，在目标函数中增加方案校验和修正方式，确保满足配网联络性和辐射性特点。
al=0.9,0.88,0.89,0.91,0.92,0.94,0.96,1.02,1.04,1.06,1.08,1.05,1.02,0.98,0.96,0.95,0.98,0.99,1.02,1.05,1.07,1.05,1.03,0.96];%不同时间负荷比例 pw=50.*[0.359165 0.378539 0.397914 0.546945 0.518629 0.497765 0.459016 0.318927 0.253353 0.4307 0.490313 0.561848 0.582712 0.603577 0.643815 0.663189 0.718331 0.631893 0.530551 0.460507 0.38152 0.371088 0.359165 0.338301]; pv=50.*[0,0,0,0,0,0.00670613838163443,0.0740392771149517,0.240411546113239,0.450742969667699,0.699735172142443,0.855341269188865,0.963572441035704,1,0.917690448475439,0.500560277957363,0.463527102400363,0.341920844346065,0.191942376214751,0.136370650849977,0,0,0,0,0]; for it=1:24 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; NodeM(:,2)=al(it).*NodeM(:,2); NodeM(7,2)=NodeM(7,2)-pw(it)/1000;%加入风电 NodeM(9,2)=NodeM(9,2)-pw(it)*2/1000;%加入风电 NodeM(18,2)=NodeM(18,2)-pv(it)/1000;%加入光伏 NodeM(25,2)=NodeM(25,2)-pv(it)/1000;%加入光伏 Dim=5; lb=ones(1,Dim); ub=[10 7 15 21 11]; % a1=rand(sizepop, Dim); % a2=repmat(ub-lb,sizepop,1); % a3=repmat(lb,sizepop,1); % pop = round(a1.* a2+ a3); % for i1=1:LL % a(1,i1)=H(i1,Swarm1(1,(it-1)*5+i1)); % state(a(1,i1),it)=0; % end [checkhl,checkgd]=kxx(Swarm1(1,(it-1)*5+1:it*5)); while checkhl==0 || checkgd==0 Swarm1(1,(it-1)*5+1:it*5)=round(rand(1,5).* (ub-lb)+ lb); [checkhl,checkgd]=kxx(Swarm1(1,(it-1)*5+1:it*5)); end for i1=1:LL a(1,i1)=H(i1,Swarm1(1,(it-1)*5+i1)); state(a(1,i1),it)=0; 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
三、程序结果
程序运行过程解集实时优化图
可私信获取程序链接~
