简介
蚁群算法（Ant Colony Optimization，简称ACO）是由Marco Dorigo于1992年发明的一种启发式算法，是通过模拟蚁群寻找食物的过程发现最短路径的行为。和粒子群算法一样，它也属于群集智能的一种。
算法思想
蚁群觅食过程中，每只蚂蚁在所走过的路径上均会释放出一种信息素，该信息素随时间的推移逐渐挥发。因此，每条路径上的信息素同时存在正负反馈两种机制。
- 正反馈：蚂蚁每次经过该路径均会释放信息素使得该路径上的信息素浓度增加；
- 负反馈：每条路径上的信息素随时间推移会逐渐挥发。
由此，我们可以判断，在起点与终点之间，当相同数量的蚂蚁初始同时经过两条不同的路径时，路径上初始信息素的浓度是相同的；不过，当路径越短时，信息素挥发时间也越短，残留信息素浓度也将越高。随后的蚂蚁将根据路径上残留信息素浓度的大小对路径进行选择 — 浓度越高，选择概率越大。最终导致信息素浓度越高的路径上蚂蚁的选择数目越多，而更多的蚂蚁也将同时导致该路径上残留信息素浓度越高（即高者越高，低者越低）。因此，在理想情况下，整个蚁群将逐渐向信息素浓度最高的路径（即最短路径）进行转移。
核心算子
该算法主要是位置x和信息素t的更新：
一、根据初始的信息素计算概率：
二、根据概率计算下一时刻的位置：
三、根据当前信息素更新下一时刻的信息素：
根据这三步形成迭代循环过程。
当然在这里我们不再去求路径问题，将其稍微修改可以计算非线性多元函数。分别将上述的概率公式修改为：
信息素更新公式修改为：
需要注意的是，初始信息素可以用目标函数值代替，直到信息素更新一轮后，使用新的信息素计算概率。
示例
我们还是以一个简单的例子来介绍算法的流程。假设我们需要寻找函数y=x12+x22+x33+x44在[1,30]之间的最大值。我们很容易就知道，当x1=x2=x3=x4=30时，该函数能取到最大值。
构建一个叫ACO的类，这个类包括算法的所有操作方法：
|
1 2 3 4 5 6 7 8 9 10 11 12 |
class ACO: def __init__(self, parameters): pass def fitness(self, ind_var): pass def update_operator(self, gen, t, t_max): pass def main(self): pass |
使用__init__()方法初始化参数，包括自变量可取的最大值，最小值，种群大小，迭代代数；使用fitness()方法作为适应度函数评估该个体的函数值，在这里就是函数y的值；使用update_operator()方法根据概率更新蚁群下一时刻的位置和信息素，挑选出当前代蚁群中最好位置作为历史记录；使用main()方法实现整个算法的循环。
__init__()方法
__init__()方法的代码如下：
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 |
def __init__(self, parameters): # 初始化 self.NGEN = parameters[0] # 迭代的代数 self.pop_size = parameters[1] # 种群大小 self.var_num = len(parameters[2]) # 变量个数 self.bound = [] # 变量的约束范围 self.bound.append(parameters[2]) self.bound.append(parameters[3]) self.pop_x = np.zeros((self.pop_size, self.var_num)) # 所有蚂蚁的位置 self.g_best = np.zeros((1, self.var_num)) # 全局蚂蚁最优的位置 # 初始化第0代初始全局最优解 temp = -1 for i in range(self.pop_size): for j in range(self.var_num): self.pop_x[i][j] = np.random.uniform(self.bound[0][j], self.bound[1][j]) fit = self.fitness(self.pop_x[i]) if fit > temp: self.g_best = self.pop_x[i] temp = fit |
初始化方法接受传入的参数，包括最大值，最小值，种群大小和迭代代数。通过这些参数初始化产生一个蚂蚁种群，并记录当前代全局最优位置。
fitness()方法
fitness()方法用来评估该位置的函数值，代码如下：
|
1 2 3 4 5 6 7 |
def fitness(self, ind_var): x1 = ind_var[0] x2 = ind_var[1] x3 = ind_var[2] x4 = ind_var[3] y = x1**2 + x2**2 + x3**3 + x4**4 return y |
update_operator()方法
update_operator()方法里包含核心算子，根据概率更新下一时刻的位置和信息素，代码如下：
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 |
def update_operator(self, gen, t, t_max): rou = 0.8 # 信息素挥发系数 Q = 1 # 信息释放总量 lamda = 1 / gen pi = np.zeros(self.pop_size) for i in range(self.pop_size): for j in range(self.var_num): pi[i] = (t_max - t[i]) / t_max # 更新位置 if pi[i] < np.random.uniform(0, 1): self.pop_x[i][j] = self.pop_x[i][j] + np.random.uniform(-1, 1) * lamda else: self.pop_x[i][j] = self.pop_x[i][j] + np.random.uniform(-1, 1) * ( self.bound[1][j] - self.bound[0][j]) / 2 # 越界保护 if self.pop_x[i][j] < self.bound[0][j]: self.pop_x[i][j] = self.bound[0][j] if self.pop_x[i][j] > self.bound[1][j]: self.pop_x[i][j] = self.bound[1][j] # 更新t值 t[i] = (1 - rou) * t[i] + Q * self.fitness(self.pop_x[i]) # 更新全局最优值 if self.fitness(self.pop_x[i]) > self.fitness(self.g_best): self.g_best = self.pop_x[i] t_max = np.max(t) return t_max, t |
需要注意的是，在更新位置的过程中有可能超过了我们定义的定义域范围，因此需要进行越界保护，强行将位置拉进定义域内，最后更新全局最好位置。
main()方法
ACO算法所有的轮子都写好后，我们接下来将它们整合到流程中。代码实现如下：
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 |
def main(self): popobj = [] best = np.zeros((1, self.var_num))[0] for gen in range(1, self.NGEN + 1): if gen == 1: tmax, t = self.update_operator(gen, np.array(list(map(self.fitness, self.pop_x))), np.max(np.array(list(map(self.fitness, self.pop_x))))) else: tmax, t = self.update_operator(gen, t, tmax) popobj.append(self.fitness(self.g_best)) print('############ Generation {} ############'.format(str(gen))) print(self.g_best) print(self.fitness(self.g_best)) if self.fitness(self.g_best) > self.fitness(best): best = self.g_best.copy() print('最好的位置：{}'.format(best)) print('最大的函数值：{}'.format(self.fitness(best))) print("---- End of (successful) Searching ----") plt.figure() plt.title("Figure1") plt.xlabel("iterators", size=14) plt.ylabel("fitness", size=14) t = [t for t in range(1, self.NGEN + 1)] plt.plot(t, popobj, color='b', linewidth=2) plt.show() |
传入循环迭代指定代数，输出目前找到的最好的位置和函数值，并画图得到整个迭代过程。
完整代码
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 |
import numpy as np import matplotlib.pyplot as plt class ACO: def __init__(self, parameters): """ Ant Colony Optimization parameter: a list type, like [NGEN, pop_size, var_num_min, var_num_max] """ # 初始化 self.NGEN = parameters[0] # 迭代的代数 self.pop_size = parameters[1] # 种群大小 self.var_num = len(parameters[2]) # 变量个数 self.bound = [] # 变量的约束范围 self.bound.append(parameters[2]) self.bound.append(parameters[3]) self.pop_x = np.zeros((self.pop_size, self.var_num)) # 所有蚂蚁的位置 self.g_best = np.zeros((1, self.var_num)) # 全局蚂蚁最优的位置 # 初始化第0代初始全局最优解 temp = -1 for i in range(self.pop_size): for j in range(self.var_num): self.pop_x[i][j] = np.random.uniform(self.bound[0][j], self.bound[1][j]) fit = self.fitness(self.pop_x[i]) if fit > temp: self.g_best = self.pop_x[i] temp = fit def fitness(self, ind_var): """ 个体适应值计算 """ x1 = ind_var[0] x2 = ind_var[1] x3 = ind_var[2] x4 = ind_var[3] y = x1 ** 2 + x2 ** 2 + x3 ** 3 + x4 ** 4 return y def update_operator(self, gen, t, t_max): """ 更新算子：根据概率更新下一时刻的位置 """ rou = 0.8 # 信息素挥发系数 Q = 1 # 信息释放总量 lamda = 1 / gen pi = np.zeros(self.pop_size) for i in range(self.pop_size): for j in range(self.var_num): pi[i] = (t_max - t[i]) / t_max # 更新位置 if pi[i] < np.random.uniform(0, 1): self.pop_x[i][j] = self.pop_x[i][j] + np.random.uniform(-1, 1) * lamda else: self.pop_x[i][j] = self.pop_x[i][j] + np.random.uniform(-1, 1) * ( self.bound[1][j] - self.bound[0][j]) / 2 # 越界保护 if self.pop_x[i][j] < self.bound[0][j]: self.pop_x[i][j] = self.bound[0][j] if self.pop_x[i][j] > self.bound[1][j]: self.pop_x[i][j] = self.bound[1][j] # 更新t值 t[i] = (1 - rou) * t[i] + Q * self.fitness(self.pop_x[i]) # 更新全局最优值 if self.fitness(self.pop_x[i]) > self.fitness(self.g_best): self.g_best = self.pop_x[i] t_max = np.max(t) return t_max, t def main(self): popobj = [] best = np.zeros((1, self.var_num))[0] for gen in range(1, self.NGEN + 1): if gen == 1: tmax, t = self.update_operator(gen, np.array(list(map(self.fitness, self.pop_x))), np.max(np.array(list(map(self.fitness, self.pop_x))))) else: tmax, t = self.update_operator(gen, t, tmax) popobj.append(self.fitness(self.g_best)) print('############ Generation {} ############'.format(str(gen))) print(self.g_best) print(self.fitness(self.g_best)) if self.fitness(self.g_best) > self.fitness(best): best = self.g_best.copy() print('最好的位置：{}'.format(best)) print('最大的函数值：{}'.format(self.fitness(best))) print("---- End of (successful) Searching ----") plt.figure() plt.title("Figure1") plt.xlabel("iterators", size=14) plt.ylabel("fitness", size=14) t = [t for t in range(1, self.NGEN + 1)] plt.plot(t, popobj, color='b', linewidth=2) plt.show() if __name__ == '__main__': NGEN = 100 popsize = 100 low = [1, 1, 1, 1] up = [30, 30, 30, 30] parameters = [NGEN, popsize, low, up] aco = ACO(parameters) aco.main() |
在if __name__ == "__main__":语句后面，我们设定所有的参数。总共跑NGEN=100代，每代的种群大小为100。
最后得到以下结果：
迭代过程如下：
总结
我们将蚁群算法改成了可以实现非线性函数最优化搜索的方式，该算法的核心就是根据概率更新下一时刻的位置和信息素。与粒子群算法非常相似，所需的参数比遗传算法要少。本文的代码能够快速运用到其他问题中，可供借鉴和参考。至此，四大启发式算法就介绍完毕。
参考文献：