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

1.软件版本

MATLAB2013b

2.本算法理论知识

传统的PSO如下：

 第一，引入惯性因子：

 对应的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.操作步骤与仿真结论

A06-05