www.pudn.com > IEEE33DG.rar > powflow_guanDG.m, change:2015-04-17,size:8163b


clc 
clear 
b=32; 
n=33; 
LL=5;           %联络开关数   
Sb=10;          %MW 
Vb=12.66;       %KV 
Zb=Vb^2/Sb;     %ohm 
check=1; 
checkhl=1; 
checkgd=1; 
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=[6 13 8 31 37];  %可行解 
%a=[6 13 8 25 3];   %不可行解 
a=[33 34 35 36 37]; %原始网络 
%a=[6 10 14 17 28];    %孤岛不可行解 
%a=[7    14     9    32    28]; 
%a=[7 14 9 32 37];  %139.4731kw文献 基于十进制编码的差分进化算法在配电网重构中的应用 
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      %要保证联络开关的首节点必须是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 
       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 ]; 
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
%添加DG的数据(powsys上的论文Study of Reconfiguration for the Distribution System With Distributed Generators结果相同) 
PDG=zeros(n,1); 
QDG=zeros(n,1); 
PDG(3,1)=0.05;  QDG(3,1)=PDG(3,1)*tan(acos(0.8)); 
PDG(6,1)=0.10;   QDG(6,1)=PDG(6,1)*tan(acos(0.8));    
PDG(24,1)=0.20;  QDG(24,1)=PDG(24,1)*tan(acos(0.9)); 
PDG(29,1)=0.10;  QDG(29,1)=PDG(29,1)*tan(acos(1)); 
SDG=-PDG-i*QDG; 
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-6) %判断收敛性(收敛精度设为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  
end 
V=V*Vb;         %反标幺 
Vm=abs(V);      %节点电压幅值 
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