【数学建模】实验二【五】Python实现蚁群算法

关于蚁群算法

这个网上的教程已经很多很详细了，比如

现代优化算法（五）： 蚁群算法

这里摘取部分，不详述

应用举例

这里继续研究 1.2 中的问题

例 已知敌方 100 个目标的经度、纬度如表 1 所示。


我方有一个基地，经度和纬度为（70,40）。假设我方飞机的速度为 1000 公里/小时。 我方派一架飞机从基地出发，侦察完敌方所有目标，再返回原来的基地。在敌方每一目 标点的侦察时间不计，求该架飞机所花费的时间（假设我方飞机巡航时间可以充分长）。

具体算法分析详见前面推的文章

数据集下载

链接：https://pan.baidu.com/s/1AmLRAZ4ypbUkUOijWmbCvQ
提取码：ygtl

Python代码

我们编写如下的 Python 程序如下：

from math import pi, sin, cos, acos

import matplotlib.pyplot as plt
import pandas as pd
import numpy as np

# # 读取数据
sj = pd.read_csv('../sj.csv').values
x = sj[..., 0:8:2].reshape(100, 1)
y = sj[..., 1:9:2].reshape(100, 1)
sj = np.concatenate([x, y], axis=1)
# 设置起始结束的经纬度为(70，40)
d1 = np.array([70, 40])
sj = np.vstack((d1, sj, d1))
# 单位化成弧度
sj = sj * pi / 180
# 初始化距离矩阵d
d = np.zeros((102, 102))
for i in range(101):
	for j in range(i + 1, 102):
		temp = cos(sj[i, 0] - sj[j, 0]) * cos(sj[i, 1]) * cos(sj[j, 1]) + sin(sj[i, 1]) * sin(sj[j, 1])
		d[i, j] = 6370 * acos(temp)
d = d + d.T

n = 102  # 地点数量
m = 75  # 蚂蚁数量
alpha = 1  # 信息素重要程度因子
beta = 5  # 启发函数重要程度因子
vol = 0.1  # 信息素发挥（volatilization）因子
Q = 1  # 常系数
Heu_F = 1 / d  # 距离倒数
Tau = np.ones((n, n))  # 信息素矩阵
Table = np.zeros((m, n)).astype(int)  # 路径记录表（蚂蚁*地点）
iter_max = 500  # 最大迭代次数
Route_best = np.zeros((iter_max, n)).astype(int)  # 每一次迭代后的最佳路线
Length_best = np.zeros(iter_max)  # 各代最佳路径长度
Length_ave = np.zeros(iter_max)  # 各代平均路径长度

for iter in range(iter_max):
	# 随机各蚂蚁起始点
	Table[:, 0] = np.random.choice(n, size=m)
	# 构建解空间
	areas = np.arange(n)
	# 逐个蚂蚁选择
	for i in range(m):
		# 逐个地点路径选择
		for j in range(1, n):
			# Table是路径记录表，每次j之前的地点都是被选择过的
			has_visited = Table[i, :j]  # 已访问的地点集合（禁忌表）
			# areas中有元素属于禁忌表中元素的时候取True，取反后变成False，产生的是逻辑数组
			allow_index = ~np.in1d(areas, has_visited)  # 剩下待访问的地点集合
			allow = areas[allow_index]
			P = np.zeros(allow.shape)
			# 计算地点间的转移概率
			for k in range(allow.size):
				tempi, tempj = has_visited[-1], allow[k]
				#  Heu_F是距离的到数据帧，所以说越远概率越低
				P[k] = (Tau[tempi, tempj] ** alpha) * (Heu_F[tempi, tempj] ** beta)
			P = P / sum(P)
			# 轮盘赌法选择下一个访问地点
			Pc = np.cumsum(P)  # 累加函数，把前几个累加到1
			target_index = np.where(Pc >= np.random.random(1))[0][0]
			target = allow[target_index]
			Table[i, j] = target
	# 首先记录上一次迭代之后的最佳路线放在第一个位置（类似一个学习的效果）
	if iter >= 1:
		Table[0] = Route_best[iter - 1]
	# 记录各个蚂蚁的路径距离
	Length = np.zeros(m)
	for i in range(m):
		Route = Table[i]  # 取出一条路径
		# for循环最后是一个点-》最后一个点之间的距离
		for j in range(n - 1):
			Length[i] = Length[i] + d[Route[j], Route[j + 1]]
		Length[i] = Length[i] + d[Route[n - 1], Route[0]]
	# 找到最短距离的值的大小和位置
	min_Length, min_index = np.min(Length), np.argmin(Length)
	# 此时迭代的路线最小值记录
	Length_best[iter] = min_Length
	Route_best[iter] = Table[min_index]
	# 此时迭代的路线平均值记录
	Length_ave[iter] = np.mean(Length)
	# 更新信息素
	Delta_Tau = np.zeros((n, n))
	# 逐个蚂蚁计算
	for i in range(m):
		# 逐个地点计算
		for j in range(n - 1):
			# 这里可以理解成每只蚂蚁带着等量的信息素，均匀撒在路线上
			Delta_Tau[Table[i, j], Table[i, j + 1]] += (Q / Length[i])
		# 路线最后一个点和起点的信息素
		Delta_Tau[Table[i, n - 1], Table[i, 0]] += (Q / Length[i])
	Tau = (1 - vol) * Tau + Delta_Tau  # 信息素挥发一部分再加上增加的
	# 清空路线图
	Table = np.zeros((m, n)).astype(int)
# 找到每次迭代之后记录下来中最小值里面的最小值
Shortest_Length, index = np.min(Length_best), np.argmin(Length_best)
Shortest_Route = Route_best[index]
print("Shortest_Route={}, Shortest_Length={}".format(Shortest_Route, Shortest_Length))
# 画出路径图
for i in range(-1,101):
	dx, dy = [sj[Shortest_Route[i]][0], sj[Shortest_Route[i + 1]][0]], [sj[Shortest_Route[i]][1], sj[Shortest_Route[i + 1]][1]]
	plt.plot(dx, dy)
plt.show()

结果截图

计算结果为 43 小时左右（取决于各参数的选择）。其中的一个巡航路径如图 5 所示。