基于蚁群算法(ACO)的函数寻优代码详解

前言

蚁群算法与遗传算法一并属于启发式算法，其原理有一定的相似性。
蚁群算法的仿生原理可以这样举例：在不远处的地上有一块奶糖，这时候你用手放个蚂蚁在地上，在无其他因素影响的情况下，这只蚂蚁会爬到奶糖那边干饭，蚂蚁从你手里出去到奶糖之间走过的路径就是一个解。
这时候问题来了，如果你想知道从自己在的位置到奶糖的最短路径（最优解）是什么，该怎么办？显然光凭那一只蚂蚁太为难人了。于是你会想到，在地上放一大把蚂蚁，这些蚂蚁在无其他因素影响的情况下，应该都会爬过去干饭，这个时候观察一下他们的路径就好了。
现在又来个问题，怎么去观察哪个路径最短呢。幸亏蚂蚁在觅食的时候会在路上留下信息素（荷尔蒙），当路径很短的时候，路径上的信息素会残留的更多，会吸引其他的蚂蚁也过来走。这样，信息素浓度最高的路径就成了蚂蚁最多的路径，也就是最短路径（最优解）。
由上可见，蚁群算法是非常适合解决TSP问题（旅行商问题）的。实际上，原始的蚁群算法就是用来解决TSP问题，不过稍加改进之后可以像遗传算法一样解决函数寻优问题。

函数寻优

问题本身很简单，四个变量x1,x2,x3,x4组成方程y = x1平方+x2平方+x3平方+x4平方，四个变量取值范围都是[1,30]，求最大值。

源代码

代码扒至网络，作为学习资料，并加上了很多注释（原先代码就几句注释，难以快速上手）。

import numpy as np
import matplotlib.pyplot as plt
#每一个个体就是一个蚂蚁
#一个种群就是一个蚁群class ACO:def __init__(self, parameters):"""Ant Colony Optimizationparameter: 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]) #加入变量约束范围，加入两个列表，bound成为2维的self.pop_x = np.zeros((self.pop_size, self.var_num))  # 所有蚂蚁的位置，pop_x是2维的self.g_best = np.zeros((1, self.var_num))  # 全局蚂蚁最优的位置，大小等于一个个体# 初始化第0代初始全局最优解temp = -1for i in range(self.pop_size):for j in range(self.var_num):#对2维列表pop_x做遍历，每一行可以视作一个个体，每个个体里面含一组完整的解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]#gbest为最优位置列表更新最大值，即最优个体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 ** 4return ydef update_operator(self, gen, t, t_max):"""更新算子：根据概率更新下一时刻的位置和信息素，挑选最好的位置保存gen是当前的代数，t是信息素列表，t_max是当前信息素列表中信息素最大的值每个个体都对应一个信息素量，当信息素相对少时该个体便大概率进行行动，迭代次数多了之后个体的优度整体提升"""rou = 0.8   # 信息素挥发系数Q = 1       # 信息释放总量，蚂蚁们工作循环一次释放的信息总量lamda = 1 / gen #lamda随着代数增加而减小，用于局部搜索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) * lamdaelse:#进行全局搜索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, tdef 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 = 200#迭代代数popsize = 100#蚁群的个体数量low = [1, 1, 1, 1]#四个变量的下界up = [30, 30, 30, 30]#四个变量的上界parameters = [NGEN, popsize, low, up]#做成参数列表aco = ACO(parameters)#参数代入aco.main()#运行

代码解析

首先明确一些代码中的基础概念：
在蚁群算法中一个个体就是一个蚂蚁，等效于遗传算法中的个体（染色体）。
在蚁群算法中，一个蚂蚁身上一般不只有一个解，像该例子中蚂蚁身上有四个解（x1、x2、x3、x4），可以等效为遗传算法中染色体上有四个基因。

初始化蚁群

该程序使用一个类对整个算法进行了封装，方便进行移植使用。关于python的模块化设计与对象化编程这里不加赘述。
构造函数（不知道这里称呼为构造函数是否合适，套用C++的说法）里有两个参数，第一个self是对象自身，不需解释；第二个参数parameters显然是一个列表，使用列表管理多个参数是很高效的方式。

def __init__(self, parameters)

在类中的属性中记录了迭代代数、蚁群大小、一个蚂蚁身上有几个变量、变量的约束范围，这些属性是依靠函数传参得到的。此外还有两个重要的列表：pop_x列表记录了当前蚁群内所有蚂蚁的数据，本质是个二维数组，行数是蚁群中蚂蚁的数量，列数是一个蚂蚁的长度（四个变量）。
g_best列表是一个蚂蚁的长度，专门用来记录当前“最强的蚂蚁”，也就是当前的最优解。

# 初始化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]) #加入变量约束范围，加入两个列表，bound成为2维的self.pop_x = np.zeros((self.pop_size, self.var_num))  # 所有蚂蚁的位置，pop_x是2维的self.g_best = np.zeros((1, self.var_num))  # 全局蚂蚁最优的位置，大小等于一个个体

初始化蚁群，也就是构建蚁群。构建的方式和遗传算法中的实数编码法是一样的，给每个蚂蚁的每个解分配边界范围内的一个随机值。当一个蚂蚁拥有初始解后，立刻对其进行适应度评估，对蚁群评估完毕后，将适应度最高的个体存起来，并称其为最优个体。

 # 初始化第0代初始全局最优解temp = -1for i in range(self.pop_size):for j in range(self.var_num):#对2维列表pop_x做遍历，每一行可以视作一个个体，每个个体里面含一组完整的解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]#gbest为最优位置列表更新最大值，即最优个体temp = fit#同时刷新最大值的比较标准

全部初始化代码：

def __init__(self, parameters):"""Ant Colony Optimizationparameter: 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]) #加入变量约束范围，加入两个列表，bound成为2维的self.pop_x = np.zeros((self.pop_size, self.var_num))  # 所有蚂蚁的位置，pop_x是2维的self.g_best = np.zeros((1, self.var_num))  # 全局蚂蚁最优的位置，大小等于一个个体# 初始化第0代初始全局最优解temp = -1for i in range(self.pop_size):for j in range(self.var_num):#对2维列表pop_x做遍历，每一行可以视作一个个体，每个个体里面含一组完整的解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]#gbest为最优位置列表更新最大值，即最优个体temp = fit#同时刷新最大值的比较标准

蚁群算法核心

update_operator函数是蚁群算法的核心，也是蚁群算法基本原理的体现。
和遗传算法不同，蚁群算法不存在淘汰机制，因此适应度并不是蚁群算法里的“硬通货”。蚁群算法里的硬通货是每个蚂蚁的“信息素”，信息素决定了每个蚂蚁的采取行动概率。
首先熟悉下蚁群算法核心算子中的几个参数。目前尚未接触自适应蚁群算法，但是从经验中可以得知这些参数可以自适应调整的。
rou是信息素挥发系数，每个蚂蚁的信息素是会挥发的，不会一直保存。信息素挥发系数用来参与计算下一个时刻蚂蚁的信息素值。
Q是每次迭代信息释放总量，也就是所有蚂蚁产生的信息素总和，这里设置为1是为了方便计算。
lamda是个变系数，用于给局部搜索做参数，显然代数越多，该系数越小，这样可以使得当代数足够多时增强解的收敛性。
pi是一个列表，该列表记录了当前每个蚂蚁的转移功率，也就是采取行动概率。
t是一个列表，该列表记录了当前每个蚂蚁的信息素值。
tmax记录了当前蚁群的最高信息素值。

         rou = 0.8   # 信息素挥发系数Q = 1       # 信息释放总量，蚂蚁们工作循环一次释放的信息总量lamda = 1 / gen #lamda随着代数增加而减小，用于局部搜索pi = np.zeros(self.pop_size)#概率表，存储每个个体的转移概率

接下来的代码是每次迭代中的核心操作。
首先，计算每一个蚂蚁的行动概率，一个蚂蚁的信息素越少，说明其优度较低，需要提高行动概率来改进优度。

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) * lamdaelse:#进行全局搜索self.pop_x[i][j] = self.pop_x[i][j] + np.random.uniform(-1, 1) * (self.bound[1][j] - self.bound[0][j]) / 2

搜索是蚁群算法中唯一更新变量值的方式，当搜索结束后蚁群内所有变量的值就不会产生变化。
一次搜索完毕后需要对全部蚂蚁更新信息素的值，同时更新蚁群中的最优解。

         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)#对信息素序列进行检索得到最大值

核心代码如下：

 def update_operator(self, gen, t, t_max):"""更新算子：根据概率更新下一时刻的位置和信息素，挑选最好的位置保存gen是当前的代数，t是信息素列表，t_max是当前信息素列表中信息素最大的值每个个体都对应一个信息素量，当信息素相对少时该个体便大概率进行行动，迭代次数多了之后个体的优度整体提升"""rou = 0.8   # 信息素挥发系数Q = 1       # 信息释放总量，蚂蚁们工作循环一次释放的信息总量lamda = 1 / gen #lamda随着代数增加而减小，用于局部搜索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) * lamdaelse:#进行全局搜索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

模块主程序

将蚁群算法的核心算子加入到迭代程序中就成了蚁群算法的主程序，popobj列表记录各代的最优质值，best记录当代最大值所处的位置。利用popobj可以观察蚁群算法的迭代收敛情况。
第一代迭代之前没有信息素，但是有适应度，因此第一代的时候直接使用适应度当做信息素，之后的迭代再使用信息素作为计算用的参数。
在主程序中使用了map函数，该函数在python写的算法中非常常见。
map函数的第一个参数是一个函数的名字，第二个参数一般是数据存储结构（列表），函数的返回值就是第二个参数中的所有变量按照第一个参数（函数）运算过后的值。

     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 ----")

主程序

python本质上是个脚本语言，不需要函数入口这种东西，也就是没有其他语言的main()函数。
而且每个python程序文件本质上都是模块，因此需要区分正在运行的哪个文件是脚本，哪个文件只是当模块。
因此在程序中增加if __name__ == '__main__':以声明：以下的代码只有在当前文件为脚本的时候能够执行；当前文件作为模块被其他文件调用时，以下代码绝不会执行。
主程序中主要设定几个超参数，并调用类中的主程序，可以看到模块化编程的主程序较为简洁。

if __name__ == '__main__':NGEN = 200#迭代代数popsize = 100#蚁群的个体数量low = [1, 1, 1, 1]#四个变量的下界up = [30, 30, 30, 30]#四个变量的上界parameters = [NGEN, popsize, low, up]#做成参数列表aco = ACO(parameters)#参数代入aco.main()#运行

运行结果

下图为蚁群内100个蚂蚁，200次迭代的试验图像结果（跑了几次挑的最好的）。
可以看到在100代左右蚁群就开始收敛，算法的收敛性目测比遗传算法更强。