www.pudn.com > IEEE33DG.rar > fitness_cgfcPQV.asv, change:2015-04-17,size:9314b


function guan=fitness_cgfcPQV(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(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)节点 
%风力发电机数据 
node=18;    %风机加入的节点号 
Pe=0.6;     %单台风机的额定有功 
NodeM(node,2)=NodeM(node,2)-Pe;      %添加DG有功 
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; 
init=NodeM(node,2); 
while(max(abs(V-V0))>1e-3) %判断收敛性(收敛精度设为1e-6) 
NodeM(node,2)=init; 
V0=V;                      %记录上一次迭代的电压值 
SN=r2*((abs(V(node,1))*Ue)^2-sqrt((abs(V(node,1))*Ue)^4-4*x3^2*Pe^2))/(2*Pe*x3^2); 
Q=([r2^2+x3*(xm+x3)*SN^2]*Pe/(SN*r2*xm))/Sb; 
 
cos=Pe/sqrt(Pe^2+(Q*Sb)^2); 
 
Qc=Pe*(sqrt(1/(cos^2)-1)-sqrt(1/(cosb^2)-1))/Sb;     %并联电容器组无功补偿 
nc=fix(Qc*Sb/Qunit); 
Qz=nc*Qunit*(abs(V(node,1))^2)/Sb; 
 
cosc=Pe/sqrt(Pe^2+((Q-Qz)*Sb)^2); 
if cosc<0.9 
    nc=nc+1; 
else 
    nc=nc-1; 
end 
Qz=nc*Qunit*(abs(V(node,1))^2)/Sb; 
cosc=Pe/sqrt(Pe^2+((Q-Qz)*Sb)^2); 
while(cosc<0.9) 
nc=nc+1; 
Qz=nc*Qunit*(abs(V(node,1))^2)/Sb; 
cosc=Pe/sqrt(Pe^2+((Q-Qz)*Sb)^2); 
end 
NodeM(node,2)=NodeM(node,2)+i*(Q-Qz); 
 
%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>30 
   guan=10000; 
   disp('k>30') 
   jy=1; 
   break 
end 
end 
%输出端并联电容器补偿的参数 
%disp('风机吸收的无功功率') 
%Q 
%disp('按照初始功率因数核算的补偿无功') 
%Qc 
%disp('投切的并联电容器容量') 
%Qz 
%disp('投切的并联电容器组数') 
%n  
%disp('风机功率因数') 
%cosc 
 
V1=V                %标幺 
V=V*Vb;            %反标幺 
Vm=abs(V);       %节点电压幅值 
%———————————————————————————————— 
%节点电压上下限约束(排除a=[7 9 2 12 3];这种状况) 
if min(Vm)<Vb*0.9||max(Vm)>Vb 
    guan=10000; 
    disp('节点电压越限') 
    jy=1; 
end 
%———————————————————————————————— 
Va=angle(V);     %节点电压弧度 
while(jy==0) 
%下面计算有功网损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; 
jy=1; 
end 
check=0; 
end 
checkhl=0; 
 end