粒子群算法
- 求解函数的最值问题
启发式算法，数模中常称为智能优化算法。
粒子群算法时智能算法中的一种，全称为粒子群优化算法（Particle Swarm Optimization, PSO）。通过模拟鸟群觅食行为而发展起来的一种基于群体协作的搜索算法。
核心：利用群体中的个体对信息的共享使整个群体的运动在问题求解空间中产生从无序到有序的演化过程
一群鸟在搜索食物：
- 所有的鸟都不知道食物在哪
- 它们知道自己的当前位置距离食物有多远
- 它们知道离食物最近的鸟的位置
这只鸟历史最佳位置pbest(d)
当前位置x(d)
上一步位置x(d-1)
当前距离食物最近的鸟的位置gbest(d)
这只鸟第d步所在的位置+第d-1步的速度*运动的时间
x(d) = x(d-1) + v(d-1)*t（每一步运动的时间t一般取1）
这只鸟第d步的速度=上一步自身的速度惯性+自我认知部分+社会认知部分
v(d) = wv(d-1) + c1r1*(pbest(d)-x(d)) + c2r2(gbest(d)-x(d))（三个部分的和）
c1:个体学习因子，也称为个体加速因子
c2:社会学习因子，也称为社会加速因子
r1，r2: [0,1]上随机数
w：称为惯性权重
粒子：优化问题的候选解
位置：候选解所在的位置
速度：候选解移动的速度
适应度：评价粒子优劣的值，一般设置为目标函数值
个体最佳位置：单个粒子迄今为止找到的最佳位置
群体最佳位置：所有粒子迄今为止找到的最佳位置
其中r1，r2是[0,1]上的随机数
求解函数的最值问题
启发式搜索——利用中间信息来改进搜索策略则称为启发式算法
盲目搜索——在搜索过程中获取的中间信息不用来改进策略
求函数在[-3,3]内的最大值
可以用[x,fval]=fmincon(@Obj_fun1, x0, A, b, Aeq, bee, x_lb, x_ub)（数学规划方法详见5种数学规划模型）
使用粒子群算法求解：
写出目标函数Obj_fun.m
function y = Obj_fun1(x)
%如果调用fmincon函数，则需要添加负号改为求最小值
y = -(11*sin(x) + 7*cos(5*x))
end设置粒子群算法的参数code1.m
n = 10; %粒子数量
narvs = 1; %变量个数（函数中有几个自变量）
c1 = 2; %每个粒子的个体学习因子
c2 = 2; %每个粒子的社会学习因子
w = 0.9; %惯性权重
K = 50; %迭代权重
vmax = 1.2; %粒子的最大速度
x_lb = -3; %x的下界
x_ub = 3; %x的上界- 增加粒子数量会增加我们找到更好结果的可能性，但会增加运算的时间
- c1，c2，w这三个量有很多改进空间
- 迭代的次数k不一定越大越好，如果现在已经找到最优解了，再增加迭代次数是没有意义的
- 这里出现了粒子的最大速度，是为了保证下一步的位置离当前位置别太远，一般取自变量可行域的10%至20%
- x的上界和下界是保证粒子不会飞出定义域
初始化粒子的位置和速度
x = zeros(n, narvs);
for I = 1:narvs
%随机初始化粒子所在的位置在定义域内
x(: , i) = x_lb(i) + (x_ub(i) - x_lb(i))*rand(n, 1);
end
%随机初始化粒子的速度，设置为[-vmax, vmax]
v = -vmax + 2*vmax.* rand(n, narvs);- 写成下标的形式，保证程序的兼容性，“适配”多元函数
- 行向量可以和矩阵点乘（2017版及后Matlab支持）
- 行向量可以和矩阵相加（2017版及后Matlab支持）
计算适应度
fit = zeros(n,1); %初始化这n个粒子的适应度全为0
for I = 1 : n %循环整个粒子群，计算每一个粒子的适应度
fit(i) = Obj_fun1(x(i, :)); %调用Obj_fun1函数来计算适应度（这里写成x(i , : )主要是为了和以后遇到的多元函数互通）
end
best = x; %初始化这n个粒子迄今为止找到的最佳位置（是一个n*narvs的向量）
ind = find(fit == max(fit), 1); %找到适应度最大的那个粒子的下标
best = x(ind, : ); %定义所有粒子迄今为止找到的最佳位置（是一个1*narvs的向量）
function y = Obj_fun1(x)
y = 11*sin(x) + 7*cos(5*x);
end- 这里的适应度实际上就是我们的目标函数
- 这里可以直接计算出pbest和gbest，在后面将用于计算粒子和速度以更新粒子的位置
循环体：更新粒子速度和位置
for d = 1 : K %开始迭代，一共迭代K次
for i = 1 : n %依次更新第i个粒子的速度与位置
%更新第i个粒子的速度
v(i, :) = w*v(i, :) + c1*rand(1)*(pbest(i,:) - x(i, :)) + c2*rand(1)*(best - x(i, :));
%如果粒子的速度超过了最大速度限制，就对其进行调整
for j = 1 : narvs
if v(i,j) < -vmax(j)
v(i,j) =- vmax(j);
else if v(i,j) > vmax(j)
v(i,j) = vmax(j)
end
end
x(i, :) = x(i, :) + v(i,:); %更新第I个粒子的位置
%如果粒子的位置超出了定义域，就对其进行调整
for j = 1 : narvs
if x(i,j) < x_lb(j)
x(i,j) = x_lb(j);
else if x(i,j) > x_ub(j)
x(i,j) = x_ub(j);
end
end
fit(i) = Obj_fun1(x(i,:)); %重新计算第i个粒子的适应度
%如果第i个粒子的适应度大于这个粒子迄今为止的最佳位置对应的适应度
if fit(i) > Obj_fun1(pbest(i,:))
pbest(i,:) = x(i,:); %那就更新第i个粒子迄今为止找到的最佳位置
end
%如果第i个粒子的适应度大于所有的粒子迄今为止找到的最佳位置对应的适应度
if Obj_fun1(pbest(i,:)) > Obj_fun1(gbest)
gbest = pbest(i,:); %那就更新所有粒子迄今为止找到的最佳位置
end
end
end
