Matlab自带的粒子群函数

Matlab自带的粒子群函数

学习参考： 清风数模

particleswarm()

**注意：**这个函数求的是最小值的，如果目标函数是求最大值则需要添加负号从而转换为求最小值。

Matlab中particleswarm函数采用的是自适应的领域模式(lbest)，搜索前期使用领域模式。后期使用全局模式

自适应的体现

如果适应度开始停滞是，粒子群搜索会从领域模式向全局模式转换，一旦适应度开始下降，则又恢复到领域模式，一面陷入局部最优。当适应度的停滞次数足够大时，惯性系数开始逐渐变小，从而有利于局部搜索。

matlab自带的粒子群算法particleswarm如何实现过程

预设参数的选取

1. 粒子个数SwarmSize：默认设置为：min{100,10*nvars}，nvars是变量个数
2. 惯性权重InertiaRange：默认设置的范围为：[0.1,1.1]，**注意：**在迭代过程中惯性权重会采取自适应措施，随着迭代过程不断调整。
3. 个体学习因子SelfAdjustmentWeight：默认设置为1.49（和压缩因子的系数几乎相同）
4. 社会学习因子SocialAdjustmentWeight：默认设置为1.49（和压缩因子的系数几乎相同）
5. 邻域内粒子的比例MinNeighborsFraction：默认设置为: 0.25，由于采取的是邻域模式，因此定义了一个“邻域最少粒子数目": minNeighborhoodSize = max{2,(粒子数目*邻域内粒子的比例)的整数部分}，在迭代开始后，每个粒子会有一个邻域，初始时邻域内的粒子个数(记为Q)就等于“邻域最少粒子数目”，后续邻域内的粒子个数Q会自适应调整。变为全局模式是把Q变大

变量初始化和适应度的计算

1. 速度初始化：最大速度就是上界和下界的差额：vmax=ub-lb；v=-vmax+2*vmax.*rand(n,narvs);
2. 位置初始化：将每个粒子的位置均匀分布在上下界约束内
3. 计算每个粒子的适应度：适应度仍设置为我们要优化的目标函数，由于该函数求解的是最小值问题，因此，最优解应为适应度最小即目标函数越小的解。
4. 初始化个体最优位置：循环还没有开始所以这里的个体最优位置就是我们初始化得到的位置
5. 初始化所有粒子的最优位置：每个粒子的适应度我们都已经求出来了，所以我们只需要找到适应度最低的那个粒子，并记其位置为所有粒子的最优位置.

更新粒子的速度和位置

在每次迭代中，我们要分别更新每一个粒子的信息。例如:对于现在要更新的粒子i ,我们要进行以下几个操作：

1. 随机生成粒子i的邻域，邻域内一共Q个粒子(包含粒子i本身)，并找到这Q个粒子中位置最佳的那个粒子，此时它的目标函数值最小，记其位置为lbest,
2. 更新粒子i的速度，公式为: v = w*v + c1*u1*(pbest-x) + c2*u2*(Ibest-x);(注:这里省略了一些上下标)，这个速度公式和基本粒子群算法最大的不同在于：这里的群体信息交流部分使用的是邻域内的最优位置，而不是整个粒子群的最优位置；
3. 更新粒子i的位置，公式为: x=x+v，这个和基本粒子群算法一样;
4. 修正位置和速度:如果粒子i的位置超过了约束，就将其位置修改到边界处;另外如果这个粒子的位置在边界处，我们还需要查看其速度是否超过了最大速度，如果超过的话将这个速度变为0。(注意，如果是多元函数的话可能只有某个分量超过了约束，我们的修改只需要针对这个分量即可)
5. 计算粒子i的适应度，如果小于其历史最佳的适应度，就更新粒子i的历史最佳位置为现在的位置;另外还需要判断粒子i的适应度是否要小于所有粒子迄今为止找到的最小适应度，如果小的话需要更新所有的粒子的最佳位置为粒子i的位置。

自适应调整参数

假设在第d次迭代过程中，所有的粒子的信息都已经更新好了，那么在开始下一次的迭代之前，需要更新模型中的参数，这里就体现了自适应过程。
规则如下：记此时所有粒子的最小适应度为a，上一次迭代完成后所有粒子的最小适应度为b。比较a和b的相对大小，如果a<b，则记flag = 1;否则记flag = 0.
如果: flag = 0，那么做下面的操作:

1. 更新c=c+1﹔这里的c表示“停滞次数计数器"，在开始迭代前就初始为0
2. 更新Q= min{Q+ minNeighborhoodSize, SwarmSize}; Q：邻域内的粒子个数；minNeighborhoodSize:邻域最少粒子数目； SwarmSize:粒子的总数

如果: flag = 1，那么做下面的操作:

1. 更新Q= minNeighborhoodSize
2. 更新c = max{c-1,0}
3. 判断c的大小，如果c<2，则更新w = 2*w；如果c >5，则更新w = w/2；这里的w是惯性权重，如果计算的结果超过了惯性权重的上界或低于下界都需要将其进行调整到边界处。

自动退出迭代循环（参考matlab官方文档说明）

如果是第一种方式退出迭代的话，我们可以将粒子群算法得到的解作为初始值，继续调用其他的函数来进行混合求解，例如我们熟悉的fmincon函数

修改函数的参数

使用CODE

%% 求解函数y = x1^2+x2^2-x1*x2-10*x1-4*x2+60在[-15,15]内的最小值（最小值为8）
narvs = 2; % 变量个数
x_lb = [-15 -15]; % x的下界(长度等于变量的个数，每个变量对应一个下界约束)，没有下界约束就-inf
x_ub = [15 15]; % x的上界，自变量的上界
[x,fval,exitflag,output] = particleswarm(@Obj_fun2, narvs, x_lb, x_ub)%重要的是前两个变量   
% x是最佳位置，fval最优值，exitflag退出标志，output一个结构体返回运行过程的信息

%% 绘制最佳的函数值随迭代次数的变化图
options = optimoptions('particleswarm','PlotFcn','pswplotbestf')   
[x,fval] = particleswarm(@Obj_fun3,narvs,x_lb,x_ub,options)

%% 展示函数的迭代过程 Ctrl+C会中断迭代过程，主要看第一列和第三列
options = optimoptions('particleswarm','Display','iter');
[x,fval] = particleswarm(@Obj_fun3,narvs,x_lb,x_ub,options)

%% 修改粒子数量，默认的是：min(100,10*nvars)
options = optimoptions('particleswarm','SwarmSize',50);
[x,fval] = particleswarm(@Obj_fun3,narvs,x_lb,x_ub,options)

%% 在粒子群算法结束后继续调用其他函数进行混合求解（单词hybrid  n.混合物合成物; adj.混合的; 杂种的;） 
options = optimoptions('particleswarm','HybridFcn',@fmincon);
[x,fval] = particleswarm(@Obj_fun3,narvs,x_lb,x_ub,options)

%% 惯性权重的变化范围，默认的是0.1-1.1
options = optimoptions('particleswarm','InertiaRange',[0.2 1.2]);
[x,fval] = particleswarm(@Obj_fun3,narvs,x_lb,x_ub,options)

%% 个体学习因子，默认的是1.49（压缩因子）
options = optimoptions('particleswarm','SelfAdjustmentWeight',2);
[x,fval] = particleswarm(@Obj_fun3,narvs,x_lb,x_ub,options)

%% 社会学习因子，默认的是1.49（压缩因子）
options = optimoptions('particleswarm','SocialAdjustmentWeight',2);
[x,fval] = particleswarm(@Obj_fun3,narvs,x_lb,x_ub,options)

%% 最大的迭代次数，默认的是200*nvars
options = optimoptions('particleswarm','MaxIterations',10000);
[x,fval] = particleswarm(@Obj_fun3,narvs,x_lb,x_ub,options)

%% 领域内粒子的比例 MinNeighborsFraction，默认是0.25 
options = optimoptions('particleswarm','MinNeighborsFraction',0.2);
[x,fval] = particleswarm(@Obj_fun3,narvs,x_lb,x_ub,options)

%% 函数容忍度FunctionTolerance, 默认1e-6, 用于控制自动退出迭代的参数
options = optimoptions('particleswarm','FunctionTolerance',1e-8);
[x,fval] = particleswarm(@Obj_fun3,narvs,x_lb,x_ub,options)

%% 最大停滞迭代数MaxStallIterations, 默认20, 用于控制自动退出迭代的参数
options = optimoptions('particleswarm','MaxStallIterations',50);
[x,fval] = particleswarm(@Obj_fun3,narvs,x_lb,x_ub,options)

%% 不考虑计算时间，同时修改三个控制迭代退出的参数
tic
options = optimoptions('particleswarm','FunctionTolerance',1e-12,'MaxStallIterations',100,'MaxIterations',100000);
[x,fval] = particleswarm(@Obj_fun3,narvs,x_lb,x_ub,options)
toc

%% 在粒子群结束后调用其他函数进行混合求解
tic
options = optimoptions('particleswarm','FunctionTolerance',1e-12,'MaxStallIterations',50,'MaxIterations',20000,'HybridFcn',@fmincon);
[x,fval] = particleswarm(@Obj_fun3,narvs,x_lb,x_ub,options)
toc

对于其他优化问题

如果需要使用粒子群算法解决非线性约束问题的

1. 要在每一次生成或者粒子移动的时候要判断是否满足约束条件
2. 在适应度函数(目标函数)上构造惩罚项

如果需要使用粒子群算法解决组合优化问题（最好用模拟退火）

1. 需要对原来粒子群算法的运动公式和运算符号进行重新定义

如果需要使用粒子群算法解决非线性约束问题的

1. 要在每一次生成或者粒子移动的时候要判断是否满足约束条件
2. 在适应度函数(目标函数)上构造惩罚项

如果需要使用粒子群算法解决组合优化问题（最好用模拟退火）

1. 需要对原来粒子群算法的运动公式和运算符号进行重新定义
2. 利用其他的智能算法和粒子群算法结合，进行混合求解