配电网粒子群算法实例

算例：
        


一、故障信息的数学表示
在上图中K表示断路器，每一个断路器上均有一个FTU装置，可以反馈断路器开关是否过流，用表示上传的故障信息，反映的是各分段开关处是否流过故障电流有故障电流为1,否则为0）。即：



因为FTU上传的信息可分为有故障信息及无故障信息两类,对于分段区间来讲也只能是有故障及无故障两种情况,所以我们可以用二进制编码规则对配电网故障定位问题进行数学建模。以上图所示辐射状配电网为例,系统拥有12个分段开关,我们可以用一串12位的二进制代码表示FTU的上传信息,作为程序的输入,1代表对应的开关有过流信息,0代表对应的开关无过流信息。同时用另一串12位的二进制代码作为程序的输出,代表对应馈线区间发生故障,代表无故障。
二、二进制PSO算法的配电网故障定位原理
应用二进制PSO算法求解配电网故障定位问题,粒子的位置代表配电网中馈线区段的状态,粒子的维数代表配电网的馈线区段总数。每一馈线区段存在0和1两种状态,0表示正常状态,1表示故障状态,馈线区段的状态为待求量。因此,N段馈线区段的状态求解就转化成N维粒子群优化求解,每个粒子的N维位置都表示为配电网N段馈线区段的潜在状态。每次迭代过程中,通过评价函数评价各粒子位置优劣,更新粒子的当前最优位置和全体粒子的最优位置,进而更新粒子的速度和位置,直到满足程序终止条件为止。最终得出的粒子群的全局最优位置就是所求的各馈线区段的实际状态。
三、构造配电网故障定位评价函数

    基于待求各馈线区段实际状态下所对应的信息应与实际上传故障信息偏差最小的原则，构造如下评价函数：



表达式的值为每个潜在解对应的适应度值,值越小表示解越优良,因此评价函数应取极小值。
式中：

为第j个开关FTU上传的故障信息,取值为1认为该开关流过了故障电流,为0则未流过；

为各开关节点的期望状态,若该开关流过了故障电流,则期望状态为1,相反,期望状态为0,表达式为各区段状态的函数。N为配电网中馈线区段的总数，

为配电网中各设备状态,为1表明设备故障,取0则设备正常。

表示一个权系数乘以故障设备数,。是根据故障诊断理论中“最小集”概念设置的权重系数,取值范围介于与1之间,表明故障区间数越少解越优,避免出现误诊断,本文中权系数取0.5。
四、求期望函数



针对不同的网络，期望函数的表达式不同，无法用统一的格式表示出来。仅以上述算例为例来说明如何求取期望函数。
在辐射型配电网中,任何区间发生故障,必将导致其上游开关流过故障过电流。这是求取期望函数的原则。假设图一的馈线区段6处发生故障,那么分段开关K1,K2,K3 ,K6,都将有过电流若馈线区段11处发生故障,则K1,K2,K3 ,K6,K10都将有过电流。由此类推我们可得到各个分段开关的期望函数:



用粒子群算法求解配电网故障定位评价函数的极小值即可。

Matlab源代码：

    clc;

    clear all;

    global Ij;

    global Vmax;                       %允许误差

    MaxNum=120;                    %粒子最大迭代次数

    D=12; %粒子的维数                    

    Particlesize=30;                    %粒子群规模

    c1=2.1;                            %每个粒子的个体学习因子，也称为加速常数

    c2=2.1;                            %每个粒子的社会学习因子，也称为加速常数

    w=1;                           %惯性因子

    Vmax=4;                        %粒子的最大飞翔速度

    S_b=randint(Particlesize,D);     %粒子所在的位置

    V=Vmax*rand(Particlesize,D);         %粒子的飞翔速度

     Ij=[1 1 1 1 0 1 1 1 0 0 1 1];%FTU 反馈的值这里就是输入，请修改这里，看结果

    fitness=F(S_b');%这里注意下，要转置

    %inline定义的适应度函数会使程序运行速度大大降低

    for i=1:Particlesize

            f(i)=F(S_b(i,:));

    end

    personalbest_x=S_b;%这里是种群里的每个粒子都认为自己是最好的

    personalbest_faval=f;

    [globalbest_favali]=min(personalbest_faval);%全局最优

    globalbest_x=personalbest_x(i,:);

    k=1;

    Start=0;

    while k<=MaxNum

        Remember=globalbest_faval;

        for i=1:Particlesize

                f(i)=F(S_b(i,:));%计算适应度函数值

            if f(i)<personalbest_faval(i) %判断当前位置是否是历史上最佳位置

                personalbest_faval(i)=f(i);

                personalbest_x(i,:)=S_b(i,:);

            end

        end%这里就是把整个种群中的所有粒子都检查了一遍，保证当前最优

        [globalbest_favali]=min(personalbest_faval);%更新种群最优

        globalbest_x=personalbest_x(i,:);

        for i=1:Particlesize %更新粒子群里每个个体的最新位置

           V(i,:)=w*V(i,:)+c1*rand*(personalbest_x(i,:)-S_b(i,:))...

               +c2*rand*(globalbest_x-S_b(i,:));%更新了速度

           for j=1:length(Ij)

               if rand(1)<Sigmoid(V(i,j))%这里体现了速度决定位置

                   S_b(i,j)=1;

               else

                    S_b(i,j)=0;

               end

           end

        end

    if globalbest_faval-Remember==0

        Start=Start+1;

        if Start==20

            disp(globalbest_x);

            break;

        end

    end

        k=k+1;

end

F.m

function fitness=F(SB)

global Ij;

IE=zeros(1,length(Ij));

IE(1)=SB(1)|SB(2)|SB(3)|SB(4)|SB(5)|SB(6)|SB(7)|SB(8)|SB(9)|SB(10)|SB(11)|SB(12);

IE(2)=SB(2)|SB(3)|SB(4)|SB(5)|SB(6)|SB(7)|SB(8)|SB(9)|SB(10)|SB(11)|SB(12);

IE(3)=SB(3)|SB(4)|SB(5)|SB(6)|SB(7)|SB(8)|SB(9)|SB(10)|SB(11)|SB(12);

IE(4)=SB(4)|SB(5);

IE(5)=SB(5);

IE(6)=SB(6)|SB(7)|SB(8)|SB(9)|SB(10)|SB(11)|SB(12);

IE(7)=SB(7)|SB(8)|SB(9)|SB(10);

IE(8)=SB(8)|SB(9)|SB(10);

IE(9)=SB(9);

IE(10)=SB(10);

IE(11)=SB(11)|SB(12);

IE(12)=SB(12);

fitness=sum(abs(Ij-IE))+0.5*sum(abs(SB));

end

Sigmoid.m

function AN=Sigmoid(V)

global Vmax;

    if V>Vmax

        AN=0.98;

    elseif V>-Vmax

       AN=1/(1+exp(-V));

    else

        AN=-0.98;

    end

end