【PSO】基于惯性因子和变异因子的PSO粒子群优化仿真

1.软件版本

MATLAB2013b

2.本算法理论知识

传统的PSO如下：

第一，引入惯性因子：

对应的MATLAB代码为：

第二，引入变异因子：

对应的MATLAB代码为：

 

3.部分核心代码

clc;
clear;
close all;

% RandStream.setDefaultStream(RandStream('mt19937ar','seed',555));
%定义粒子规模
popsize = 20;
num     = 5; 
%注意，如下两组参数，在实际中，可有配电网的公式得到，这里重点是对PSO在可靠性满足的情况下进行费用优化，故直接通过参数进行输入
%定义价格初始参数
%F = min(alpha*ASUI + beta*f);
%F           = [0.025,0.075,0.04,0.035,0.05];

%计算f = CS + CM + CL
%CS开关设备投资费用
%CS开关设备投资费用
M = 5;            %开关类型的总数
N = [2 10 3 10 24];      %每个类型的开关的安装台数,该参数的设置和M值有关
Cs= [0.0003 0.0004 0.0005 0.0006 0.0005];%每个类型的开关的单价,该参数的设置和M值有关
i = 0.6;          %贴现率
P = [2 2 2 2 2];   %第j中开关的使用寿命,该参数的设置和M值有关
for j = 1:M
    CS_sub(j) = N(j) * Cs(j) * ( ((1+i)^P(j))*i )/( ((1+i)^P(j)) - 1 );
end
CS = 0;
for j = 1:M
    CS = CS + CS_sub(j);%公式4.3
end
%CM
%CM
Cso = 0.1;%投资费用
r   = 0.3; %运行费用所在投资费用的比例
CM = Cso * r;
%CL
%CL
%系统年停电损失费用
LP = num;
T  = 4;
for i = 1:LP
    for j = 1:T
        EENS(i,j) = 0.05*N(i)*rand(); 
        C(i,j)    = 0.05*rand(); 
    end
end
for i = 1:LP
    for j = 1:T
        CL_sub(i,j) = EENS(i,j)*C(i,j);    
    end
    CL(i) = sum(CL_sub(i,:));
end
for i = 1:LP
    f(i) = CS  + CM + CL(i);
end

%ASUI
NT = 1000;
for i = 1:LP
    U(i) = round(100*rand);
    ASUI(i) = U(i)*N(i)*100 / NT/8760; 
end
alpha = 0.5;
beta  = 0.5;
for i = 1:LP
    F(i) = alpha*ASUI(i) + beta*f(i);
end

%在不同价格指数下的可靠性指标
Reliability =[0.95,0.92,0.75,0.88,0.83];
F           = F';
r0          = 0.9;
Pop_sub     = zeros(popsize,num+2);                   
New_pop_sub = zeros(popsize,num+2);    
PSO_V         = zeros(popsize,num); 
%在论文中，这两个参数定义为2.05
c1          = 2.05;%加速度常数
c2          = 2.05;%加速度常数
MaxIteration= 50;%最大迭代次数
%最大最小速度
V_max       =  1;
V_min       = -1;
%最大最小个体
P_max       = 3;
P_min       = 0 ;
%% 粒子参数初始化
%粒子参数初始化，且初始化得满足可靠性指标
[BEST,Pop_sub,New_pop_sub] = PSO_initial(Pop_sub,New_pop_sub,popsize,num,Reliability,F,r0);

%速度初始化
for i=1:popsize
    for j=1:num
        PSO_V(i,j)=rands(1);
    end
end

%% 主程序，通过PSO进行费用优化
ws = zeros(MaxIteration,1);

for k=1:MaxIteration
    k
    ws(k) = 0.9-0.5*k/MaxIteration;%惯性权重,对应论文公式4.17

    %最大值计算
    for i=1:popsize                                                      
        if Pop_sub(i,num+2) < New_pop_sub(i,num+2)  
           New_pop_sub(i,:) = Pop_sub(i,:);  
        end
        
        if New_pop_sub(i,num+2) < BEST(num+2)  
           BEST = New_pop_sub(i,:);  
        end
    end
        
    %将计算状态进行保存
    PSO_state(k,1) = k;
    
    PSO_state(k,2) = BEST(num+2);
    
    PSO_state(k,3) = min(Pop_sub(:,num+2));
    
    PSO_state(k,4) = max(Pop_sub(:,num+2));
    
    %更新位置与速度  
    for i=1:popsize
        for j=1:num
            %更新速度,对应论文公式4.15
            PSO_V(i,j)     = ws(k)*PSO_V(i,j) + c1*rand*(New_pop_sub(i,j)-Pop_sub(i,j)) + c2*rand*(BEST(j)-Pop_sub(i,j)); 
            
            %对速度进行限速
            if PSO_V(i,j) > V_max
               PSO_V(i,j) = V_max;%小于最大速度
            end
            if PSO_V(i,j) < V_min
               PSO_V(i,j) = V_min;%大于最小速度    
            end            
       
            
            %更新位置,对应论文公式4.16
            Pop_sub(i,j) = Pop_sub(i,j)+PSO_V(i,j);                           
            
            %粒子变异
            tmp{k} = Pop_sub(i,j);
            if k > 1
               %随机选择两个点作为父母点
               p = rand(1,1);            
               Pop_sub(i,j) = p*tmp{k-1} + (1-p)*tmp{k-1};
            end
            
            %由父母例子产生子例子，并带有5%的变异概率
            if rand>0.95
               %变异
               Pop_sub(i,j) = 0.3*rand;
            end           
            

            Pop_sub(i,j)=ceil(Pop_sub(i,j));
            

            if Pop_sub(i,j) > V_max
               Pop_sub(i,j) = V_max;
            end
            if Pop_sub(i,j) < V_min
               Pop_sub(i,j) = V_min;  
            end           
            
        end
    end
    
     %保证可靠性的要求，同时计算价格
     for i=1:popsize,
         
         Pop_sub(i,num+1)=1;
         for j=1:num
             Pop_sub(i,num+1)=Pop_sub(i,num+1)*(1-(1-Reliability (j))^Pop_sub(i,j));
         end
         %保证可靠性要求
         if Pop_sub(i,num+1)<r0  
             Pop_sub(i,:)=func_Reliability(Pop_sub(i,:),F,Reliability ,r0,num);  
         end   
         
         Pop_sub(i,num+2)=Pop_sub(i,1:5)*F;
     end
     pause(0.1)
end

figure;
plot(1:MaxIteration,ws,'b-*');
title('惯性权重变化曲线');
grid on;

figure
plot(PSO_state(:,1),PSO_state(:,2),'LineWidth',4);
xlabel('迭代次数');
ylabel('费用');
grid on;
 

4.操作步骤与仿真结论

5.参考文献

[1]陈发堂, 吕新慧, 鄢学锋. 配电网开关优化配置智能化算法与应用的研究[J]. 华东电力, 2009(10):4.A06-05