粒子群算法

通过上一节我们得知遗传算法是一种智能的随机搜索算法，胡乱随机的过程中隐隐透着一丝章法，这是借鉴了大自然的生存法则。而今天登场的粒子群算法也源于日常生活。

当地上掉落一坨蜂蜜时，我们希望找到它的方位，怎么办呢？于是乎我们放出100只蚂蚁，让他们在这片广袤无垠的土地上随机搜索。但并不是完全随机，当有一只蚂蚁嗅到了极高浓度的蜂蜜气息时，它会对其他蚂蚁大喊：“Come here！“，于是乎所有蚂蚁会往它的方位靠拢。然而每一只蚂蚁也是一个单独的个体，它们会有自己的搜索记忆，对于自己搜索路径上出现过最高浓度的方位也一直念念不忘，因此每一只蚂蚁还会向自己的历史最强方位靠拢。由此每个蚂蚁的运动方式都是通过运动惯性+个体认知+群体协作决定的。

---

Ω

基础的算法流程非常简单易懂：

1. 有必要的话需对优化问题的可行解进行编码

2. 初始化各粒子的位置和速度，每个位置坐标对应一个可行解

3. 指定粒子适应度(蜂蜜气息浓度)计算方式，计算各个粒子的适应值

4. 与个体最高适应值进行比较，更高则进行替换，并记录当前位置

5. 与群体最高适应值进行比较，更高则进行替换，并记录当前位置

6. 更新当前粒子速度

   其中w为惯性系数，表示维持当前速度的趋势因子；均为每次更新时产生的随机数；则为对个体最优、群体最优方位的学习因子；分别为个体最优位置坐标、群体最优位置坐标

7. 更新当前粒子方位

那么在算法开始前，我们往往需要指定【粒子个数、最大迭代次数、学习因子*2、惯性系数】，另外由于可行解域往往是一块有限的区域，因此我们不能让速度过大，需要设定一个速度上限，而这个上限往往通过方位的最大绝对值进行设定，即

注意算法中的v，x都可以是n维向量，下面请欣赏我随便应用的一个案例。

---

ƺ

这个案例是Gurobi求解器官方给出的，一个建塔问题，官方文档见这。

简单地说，就是设计一个通讯塔的选址方案，每个塔的搭建费用和覆盖的地区不同，不同地区所涵盖的人口也不同，由于资金有预算限制，因此我们需要在预算范围内求出能够覆盖最大人数的搭建方案。

∻ 导入数据

原题中只有6座塔，那么意味着只有64种情况，我30个粒子就已经占了一大半，这还迭代个毛球，于是乎我把塔数改成了10，同时更改了个别塔所涵盖的地区。下面是基本信息的导入：

import numpy as np
import random as rd

# 各地区人口数量
pop = np.array([523, 690, 420, 1010, 1200, 850, 400, 1008, 950])
# 建塔费用
cost = np.array([4.2, 6.1, 5.2, 5.5, 4.8, 9.2, 7.6, 5.8, 8.7, 4.9])
tn = len(cost)
# 资金预算
budget = 20
# 每座塔涵盖的地区
cover = np.array([[False] * len(pop)] * tn)
cover[0][[0, 5]] = True
cover[1][[8]] = True
cover[2][[3, 4]] = True
cover[3][[2, 6]] = True
cover[4][[6, 8]] = True
cover[5][[0]] = True
cover[6][[1, 5]] = True
cover[7][[2]] = True
cover[8][[3, 7]] = True
cover[9][[7]] = True

∻ 特殊化处理

在本题中没有直接的坐标，我们需要对每个方案进行编码，由于每座塔只有建/不建两种情况，那么一个10位二进制数即可表示一种方案

但是这里有最大资金约束，意味着并不是随便一个10位二进制数就是可行解。这里我采取的方法是，初始粒子生成的都是可行解，但对于搜索过程中产生的不可行解，我们将其适应值设为0，从而对搜索过程不会造成太大影响。

∻ 参数设置

接下来设置各个参数：

# 参数设定
pn = 20           # 粒子数量
iter = 500        # 迭代次数
maxx = 2**tn-1    # x_max(全为1)
r_vx = 0.5        # x_max/v_max
maxv = r_vx*maxx  # v_max
w = 0.5           # 惯性系数
c1 = 5            # 个体学习因子
c2 = 3            # 群体学习因子

∻ 解码函数

在计算下一步坐标时我们不直接对一个十维向量进行操作，毕竟每一位就两个选择也不能存在小数。我先将二进制数转换为十进制数，在进行坐标的计算，即这里的坐标和速度都是1维的，那么我们就需要一个将十进制数转换为十位二进制列表的解码函数：

def decode(p):
    b = list(bin(p)[2:])
    a = ['0'] * (tn - len(b))
    a.extend(b)
    a = list(map(int, a))
    return a

转换时别忘了把前置0补全。

∻ 粒子群初始化

初始化起点坐标和速度，这里构造的每个粒子都是可行解：

# 初始化粒子群
x = []
v = []
for i in range(pn):
    a = [1] * tn
    while sum(cost * a) > budget:
        c = rd.randint(0, maxx)
        a = decode(c)
    x.append(c)
    v.append(r_vx * c)

∻ 适应度函数

很显然本问中的适应度函数即为涵盖的人口总数，我们只需判断每个地区是否有塔覆盖即可，对于不可行解则直接返回0。

⚠️在else里会出大问题，由于idx为空，所以和pop相乘的也是个None。

# 适应度函数
def adpt(p):
    a = decode(p)
    if sum(cost * a) > budget or p == 0:
        return 0
    else:
        idx = list(np.where(np.array(a) == 1)[0])
        return sum(list(map(bool, sum(cover[idx]))) * pop)

∻ 迭代搜索

这里注意当速度和坐标超出上下限时需要及时调整。

# 粒子群搜索
maxp = 0
optx = 0
maxp_idv = list(map(adpt, x))
optx_idv = x
record_p = []
record_iter = []
for i in range(iter):
    for j in range(pn):
        v[j] = w * v[j] + c1 * rd.random() * (optx - x[j]) + c2 * rd.random() * (optx_idv[j] - x[j])
        if abs(v[j]) > maxv:
            v[j] = maxv if v[j] > 0 else -maxv
        x[j] = round(x[j] + v[j])
        if x[j] < 0:
            x[j] = 0
        elif x[j] > maxx:
            x[j] = maxx
        new_adpt = adpt(x[j])
        if new_adpt > maxp_idv[j]:
            maxp_idv[j] = new_adpt
            optx_idv[j] = x[j]
            if new_adpt > maxp:
                maxp = new_adpt
                optx = x[j]
                record_p.append(new_adpt)
                record_iter.append(i)

print("最大覆盖人数：\t", maxp)
print("最佳建塔计划：\t", decode(optx))

最终求得的结果：

这与Gurobi给出的最优解结果一致：

∻ 绘制搜索过程

import matplotlib.pyplot as plt

# 设置字体
plt.rc('font', family='Times New Roman')
# 设置图像的像素
plt.rcParams['figure.dpi'] = 150
# 设置字体的颜色
plt.rcParams['text.color'] = 'black'

plt.plot(record_iter, record_p, color='b', linestyle='-.', marker='*', label='Survived')
plt.ylabel('max Population', fontsize=13)
plt.xlabel('iteration num', fontsize=13)
plt.xticks(ticks=record_iter)
plt.show()

用matplotlib库我们可以绘制出最大人口数的变化过程，由于运行速度很快，所以我多运行了几次：

咳咳，虽然吧算法鲁棒性不高，不过就以上四种情况来看，五百轮迭代是足够了，最坏也就408次。只能说调参真的是一门技术活，我觉得我这里的参数肯定不是最优的，或许用智能算法搜索智能算法的参数会是个好主意(禁止套娃\doge)。