万字长文带你了解蚁群算法及求解复杂约束问题【源码实现】

蚁群算法

蚁群算法是一种源于大自然生物世界的新的仿生进化算法，由意大利学者M. Dorigo, V. Maniezzo和A. Colorni等人于20世纪90年代初期通过模拟自然界中蚂蚁集体寻径行为而提出的一种基于种群的启发式随机搜索算法"。蚂蚁有能力在没有任何提示的情形下找到从巢穴到食物源的最短路径，并且能随环境的变化，适应性地搜索新的路径，产生新的选择。其根本原因是蚂蚁在寻找食物时,能在其走过的路径.上释放一种特殊的分泌物一信息素 (也称外激素)，随着时间的推移该物质会逐渐挥发，后来的蚂蚁选择该路径的概率与当时这条路径上信息素的强度成正比。当一条路径上通过的蚂蚁越来越多时，其留下的信息素也越来越多，后来蚂蚁选择该路径的概率也就越高，从而更增加了该路径上的信息素强度。而强度大的信息素会吸引更多的蚂蚁，从而形成一种正反馈机制。通过这种正反馈机制，蚂蚁最终可以发现最短路径。
蚁群算法是对自然界蚂蚁的寻径方式进行模拟而得出的一种仿生算法。蚂蚁在运动过程中，能够在它所经过的路径.上留下信息素进行信息传递，而且蚂蚁在运动过程中能够感知这种物质，并以此来指导自己的运动方向。因此，由大量蚂蚁组成的蚁群的集体行为便表现出一种信息正反馈现象:某一路径上走过的蚂蚁越多，则后来者选择该路径的概率就越大。

真实蚁群觅食过程

为了说明蚁群算法的原理，先简要介绍一下蚂蚁搜寻食物的具体过程。在自然界中，蚁群在寻找食物时，它们总能找到一- 条从食物到巢穴之间的最优路径。这是因为蚂蚁在寻找路径时会在路径上释放出一种特殊的信息素。蚁群算法的信息交互主要是通过信息素来完成的。蚂蚁在运动过程中，能够感知这种物质的存在和强度。初始阶段，环境中没有信息素的遗留，蚂蚁寻找事物完全是随机选择路径，随后寻找该事物源的过程中就会受到先前蚂蚁所残留的信息素的影响，其表现为蚂蚁在选择路径时趋向于选择信息素浓度高的路径。同时，信息素是一种挥发性化学物，会随着时间的推移而慢慢地消逝。如果每只蚂蚁在单位距离留下的信息素相同，那对于较短路径上残留的信息素浓度就相对较高，这被后来的蚂蚁选择的概率就大,从而导致这条短路径上走的蚂蚁就越多。而经过的蚂蚁越多，.该路径.上残留的信息素就将更多，这样使得整个蚂蚁的集体行为构成了信息素的正反馈过程，最终整个蚁群会找出最优路径。

若蚂蚁从A点出发，速度相同，食物在D点，则它可能随机选择路线ABD
或ACD。假设初始时每条路线分配一只蚂蚁，每个时间单位行走一步。图所示为经过8个时间单位时的情形:走路线ABD的蚂蚁到达终点;而走路线ACD的蚂蚁刚好走到C点，为一半路程。
图2表示从开始算起，经过16个时间单位时的情形:走路线ABD的蚂蚁
到达终点后得到食物又返回了起点A，而走路线ACD的蚂蚁刚好走到D点。
假设蚂蚁每经过一-处所留下的信息素为1个单位，则经过32个时间单位后，
所有开始一起出发的蚂蚁都经过不同路径从D点取得了食物。此时ABD的路线往返了2趟，每一处的信息素为4个单位:而ACD的路线往返了一趟，每一处
的信息素为2个单位，其比值为2: 1。
寻找食物的过程继续进行，则按信息素的指导，蚁群在ABD路线上增派一只蚂蚁(共2只)，而ACD路线上仍然为一只蚂蚁。再经过32个时间单位后，两条线路.上的信息素单位积累为12和4，比值为3: 1。
若按以上规则继续,蚁群在ABD路线上再增派-一只蚂蚁(共3只),而ACD
路线.上仍然为一只蚂蚁。再经过32个时间单位后，两条线路上的信息素单位积累为24和6，比值为4: 1。
若继续进行，则按信息素的指导，最终所有的蚂蚁都会放弃ACD路线，而选择ABD路线。这也就是前面所提到的正反馈效应。
信息素的更新方式有两种: 一是挥发，也就是所有路径上的信息素以一定的比率减少，模拟自然蚁群的信息素随时间挥发的过程:二是增强，给评价值“好”(有蚂蚁走过)的边增加信息素。

蚁群算法流程

以TSP问题为例

我自己画的下面算例的流程图

还有些相关术语，自己见代码吧。这个代码简单。这个代码都看不懂，我劝你放弃挣扎，躺平。

复杂约束算例1

求函数f(x, y)=20(x^2- y^2) -(1-y)^2 -3(1 +y)^2+ 0.3的最小值，其中x的取值范围为[-5，5], y的取值范围为[-5, 5]。这是一个有多个局部极值的函数.

matlab版代码

%%%%%%%%%%%%%%%%%%%%蚁群算法求函数极值%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%初始化%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;               %清除所有变量
close all;               %清图
clc;                     %清屏
m=20;                    %蚂蚁个数
G_max=200;               %最大迭代次数
Rho=0.9;                 %信息素蒸发系数
P0=0.2;                  %转移概率常数
XMAX= 5;                 %搜索变量x最大值
XMIN= -5;                %搜索变量x最小值
YMAX= 5;                 %搜索变量y最大值
YMIN= -5;                %搜索变量y最小值
%%%%%%%%%%%%%%%%%随机设置蚂蚁初始位置%%%%%%%%%%%%%%%%%%%%%%
for i=1:mX(i,1)=(XMIN+(XMAX-XMIN)*rand);X(i,2)=(YMIN+(YMAX-YMIN)*rand);Tau(i)=func(X(i,1),X(i,2));
end
step=0.1;                %局部搜索步长
for NC=1:G_maxlamda=1/NC;[Tau_best,BestIndex]=min(Tau);%%%%%%%%%%%%%%%%%%计算状态转移概率%%%%%%%%%%%%%%%%%%%%for i=1:mP(NC,i)=(Tau(BestIndex)-Tau(i))/Tau(BestIndex);end%%%%%%%%%%%%%%%%%%%%%%位置更新%%%%%%%%%%%%%%%%%%%%%%%%for i=1:m%%%%%%%%%%%%%%%%%局部搜索%%%%%%%%%%%%%%%%%%%%%%if P(NC,i)<P0temp1=X(i,1)+(2*rand-1)*step*lamda;temp2=X(i,2)+(2*rand-1)*step*lamda;else%%%%%%%%%%%%%%%%全局搜索%%%%%%%%%%%%%%%%%%%%%%%temp1=X(i,1)+(XMAX-XMIN)*(rand-0.5);temp2=X(i,2)+(YMAX-YMIN)*(rand-0.5);end%%%%%%%%%%%%%%%%%%%%%边界处理%%%%%%%%%%%%%%%%%%%%%%%if temp1<XMINtemp1=XMIN;endif temp1>XMAXtemp1=XMAX;endif temp2<YMINtemp2=YMIN;endif temp2>YMAXtemp2=YMAX;end%%%%%%%%%%%%%%%%%%蚂蚁判断是否移动%%%%%%%%%%%%%%%%%%if func(temp1,temp2)<func(X(i,1),X(i,2))X(i,1)=temp1;X(i,2)=temp2;endend%%%%%%%%%%%%%%%%%%%%%%%更新信息素%%%%%%%%%%%%%%%%%%%%%%%for i=1:mTau(i)=(1-Rho)*Tau(i)+func(X(i,1),X(i,2));end[value,index]=min(Tau);trace(NC)=func(X(index,1),X(index,2));
end
[min_value,min_index]=min(Tau);
minX=X(min_index,1);                           %最优变量
minY=X(min_index,2) ;                        %最优变量
minValue=func(X(min_index,1),X(min_index,2));  %最优值
figure
plot(trace)
xlabel('搜索次数');
ylabel('适应度值');
title('适应度进化曲线')
%%%%%%%%%%%%%%%%%%%%%%%适应度函数%%%%%%%%%%%%%%%%%%%%%%%
function value=func(x,y)
value =20*(x^2-y^2)^2-(1-y)^2-3*(1+y)^2+0.3;
end

python版代码

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
# @Author: yudengwu(余登武）
# @Date  : 2021/5/26
#@email:1344732766@qq.com
import numpy as np
from tqdm import tqdm#进度条设置
import matplotlib.pyplot as plt
import matplotlib as mpl
import matplotlib; matplotlib.use('TkAgg')
mpl.rcParams['font.sans-serif'] = ['SimHei']  # 指定默认字体
mpl.rcParams['axes.unicode_minus'] = False  # 解决保存图像是负号'-'显示为方块的问题#============蚁群算法求函数极值================#=======适应度函数===
def func(x,y):value = 20*np.power(x*x-y*y,2)-np.power(1-y,2)-3*np.power(1+y,2)+0.3return value
#=======初始化参数====
m=20                   #蚂蚁个数
G_max=200              #最大迭代次数
Rho=0.9                #信息素蒸发系数
P0=0.2                 #转移概率常数
XMAX= 5               #搜索变量x最大值
XMIN= -5              #搜索变量x最小值
YMAX= 5               #搜索变量y最大值
YMIN= -5              #搜索变量y最小值
X=np.zeros(shape=(m,2)) #蚁群 shape=(20, 2)
Tau=np.zeros(shape=(m,)) #信息素
P=np.zeros(shape=(G_max,m)) #状态转移矩阵
fitneess_value_list=[] #迭代记录最优目标函数值
#==随机设置蚂蚁初始位置==
for i in range(m):#遍历每一个蚂蚁X[i,0]=np.random.uniform(XMIN,XMAX,1)[0] #初始化xX[i,1]=np.random.uniform(YMIN,YMAX,1)[0] #初始化yTau[i]=func(X[i,0],X[i,1])step=0.1;                #局部搜索步长
for NC in range(G_max):#遍历每一代lamda=1/(NC+1)BestIndex=np.argmin(Tau) #最优索引Tau_best=Tau[BestIndex] #最优信息素#计算状态转移概率for i in range(m):#遍历每一个蚂蚁P[NC,i]=np.abs((Tau_best-Tau[i]))/np.abs(Tau_best)+0.01 #即例最优信息素的距离#=======位置更新==========for i in range(m):  # 遍历每一个蚂蚁#===局部搜索===if P[NC,i]<P0:temp1 = X[i, 0] + (2 * np.random.random() - 1) * step * lamda # x(2 * np.random.random() - 1) 转换到【-1,1】区间temp2 = X[i,1] + (2 * np.random.random() - 1) * step * lamda #y#===全局搜索===else:temp1 = X[i, 0] + (XMAX - XMIN) * (np.random.random() - 0.5)temp2 = X[i, 0] + (YMAX - YMIN) * (np.random.random() - 0.5)#=====边界处理=====if temp1 < XMIN:temp1 =XMINif temp1 > XMAX:temp1 =XMAXif temp2 < XMIN:temp2 =XMINif temp2 > XMAX:temp2 =XMAX##判断蚂蚁是否移动(选更优）if func(temp1, temp2) < func(X[i, 0], X[i, 1]):X[i, 0] = temp1X[i, 1]= temp2#=====更新信息素========for i in range(m):  # 遍历每一个蚂蚁Tau[i] = (1 - Rho) * Tau[i] + func(X[i, 0], X[i, 1]) #(1 - Rho) * Tau[i] 信息蒸发后保留的index=np.argmin(Tau)#最小值索引value=Tau[index]#最小值fitneess_value_list.append(func(X[index,0],X[index,1])) #记录最优目标函数值#打印结果
min_index=np.argmin(Tau)#最优值索引
minX=X[min_index,0]  #最优变量x
minY=X[min_index,1]  #最优变量y
minValue=func(X[min_index,0],X[min_index,1])  #最优目标函数值print('最优变量x',minX,end='')
print('最优变量y',minY,end='\n')
print('最优目标函数值',minValue)plt.plot(fitneess_value_list,label='迭代曲线')
plt.legend()
plt.show()

复杂约束算例2

问题如下：其中a=10

python版求解

import pandas as pd
import numpy as np
from tqdm import tqdm#进度条设置
import matplotlib.pyplot as plt
import matplotlib; matplotlib.use('TkAgg')
from pylab import *
mpl.rcParams['font.sans-serif'] = ['SimHei']
mpl.rcParams['axes.unicode_minus'] = False#=======================定义一些函数==========================def calc_f(X):"""计算粒子的的适应度值，也就是目标函数值，X 的维度是 size * 2 """A = 10pi = np.pix = X[0]y = X[1]return 2 * A + x ** 2 - A * np.cos(2 * pi * x) + y ** 2 - A * np.cos(2 * pi * y)def calc_e(X):"""计算蚂蚁的惩罚项，X 的维度是 size * 2 """ee = 0"""计算第一个约束的惩罚项"""e1 = X[0] + X[1] - 6ee += max(0, e1)"""计算第二个约束的惩罚项"""e2 = 3 * X[0] - 2 * X[1] - 5ee += max(0, e2)return ee#子代和父辈之间的选择操作
def update_best(parent,parent_fitness,parent_e,child,child_fitness,child_e):"""针对不同问题，合理选择惩罚项的阈值。本例中阈值为0.00001:param parent: 父辈个体:param parent_fitness:父辈适应度值:param parent_e    ：父辈惩罚项:param child:  子代个体:param child_fitness 子代适应度值:param child_e  ：子代惩罚项:return: 父辈 和子代中较优者、适应度、惩罚项"""# 规则1，如果 parent 和 child 都没有违反约束，则取适应度小的if parent_e <= 0.00001 and child_e <= 0.00001:if parent_fitness <= child_fitness:return parent,parent_fitness,parent_eelse:return child,child_fitness,child_e# 规则2，如果child违反约束而parent没有违反约束，则取parentif parent_e < 0.00001 and child_e  >= 0.00001:return parent,parent_fitness,parent_e# 规则3，如果parent违反约束而child没有违反约束，则取childif parent_e >= 0.00001 and child_e < 0.00001:return child,child_fitness,child_e# 规则4，如果两个都违反约束，则取适应度值小的if parent_fitness <= child_fitness:return parent,parent_fitness,parent_eelse:return child,child_fitness,child_e#=======================定义一些参数==========================
m=20                   #蚂蚁个数
G_max=200              #最大迭代次数
Rho=0.9                #信息素蒸发系数
P0=0.2                 #转移概率常数
XMAX= 2               #搜索变量x最大值
XMIN= 1              #搜索变量x最小值
YMAX= 0               #搜索变量y最大值
YMIN= -1              #搜索变量y最小值
step=0.1                #局部搜索步长
P=np.zeros(shape=(G_max,m)) #状态转移矩阵
fitneess_value_list=[] #迭代记录最优目标函数值#=======================初始化蚂蚁群体位置和信息素==========================
def initialization():""":return: 初始化蚁群和初始信息素"""X = np.zeros(shape=(m, 2))  # 蚁群 shape=(20, 2)Tau = np.zeros(shape=(m,))  # 信息素for i in range(m):  # 遍历每一个蚂蚁X[i, 0] = np.random.uniform(XMIN, XMAX, 1)[0]  # 初始化xX[i, 1] =np.random.uniform(YMIN, YMAX, 1)[0]  # 初始化yTau[i] = calc_f(X[i])#计算信息素return X,Tau#位置更新
def position_update(NC,P,X):""":param NC: 当前迭代次数:param P: 状态转移矩阵:param X: 蚁群:return: 蚁群X"""lamda = 1 / (NC + 1)# =======位置更新==========for i in range(m):  # 遍历每一个蚂蚁# ===局部搜索===if P[NC, i] < P0:temp1 = X[i, 0] + (2 * np.random.random() - 1) * step * lamda  # x(2 * np.random.random() - 1) 转换到【-1,1】区间temp2 = X[i, 1] + (2 * np.random.random() - 1) * step * lamda  # y# ===全局搜索===else:temp1 = X[i, 0] + (XMAX - XMIN) * (np.random.random() - 0.5)temp2 = X[i, 0] + (YMAX - YMIN) * (np.random.random() - 0.5)# =====边界处理=====if (temp1 < XMIN) or (temp1 > XMAX):temp1 = np.random.uniform(XMIN, XMAX, 1)[0]  # 初始化xif (temp2 < YMIN) or (temp2 > YMAX):temp2 = np.random.uniform(YMIN, YMAX, 1)[0]  # 初始化y##判断蚂蚁是否移动(选更优）#子代蚂蚁children=np.array([temp1,temp2])#子代个体蚂蚁children_fit=calc_f(children) #子代目标函数值children_e=calc_e(children) #子代惩罚项parent=X[i]#父辈个体蚂蚁parent_fit=calc_f(parent)#父辈目标函数值parent_e=calc_e(parent)#父辈惩罚项pbesti, pbest_fitness, pbest_e = update_best(parent, parent_fit, parent_e, children, children_fit,children_e)X[i]=pbestireturn X#信息素更新
def Update_information(Tau,X):""":param Tau: 信息素:param X: 蚂蚁群:return: Tau信息素"""for i in range(m):  # 遍历每一个蚂蚁Tau[i] = (1 - Rho) * Tau[i] + calc_f(X[i]) #(1 - Rho) * Tau[i] 信息蒸发后保留的return Taudef main():X,Tau=initialization() #初始化蚂蚁群X 和信息素 Taufor NC in tqdm(range(G_max)):  # 遍历每一代BestIndex = np.argmin(Tau)  # 最优索引Tau_best = Tau[BestIndex]  # 最优信息素# 计算状态转移概率for i in range(m):  # 遍历每一个蚂蚁P[NC, i] = np.abs((Tau_best - Tau[i])) / np.abs(Tau_best) + 0.01  # 即离最优信息素的距离# =======位置更新==========X=position_update(NC,P,X) #X.shape=(20, 2)# =====更新信息素========Tau=Update_information(Tau, X)# =====记录最优目标函数值========index = np.argmin(Tau)  # 最小值索引value = Tau[index]  # 最小值fitneess_value_list.append(calc_f(X[index]))  # 记录最优目标函数值# 打印结果min_index = np.argmin(Tau)  # 最优值索引minX = X[min_index, 0]  # 最优变量xminY = X[min_index, 1]  # 最优变量yminValue = calc_f(X[min_index])  # 最优目标函数值print('最优变量x', minX, end='')print('最优变量y', minY, end='\n')print('最优目标函数值', minValue)print('最优变量对应的惩罚项',calc_e(X[min_index]))plt.plot(fitneess_value_list, label='迭代曲线')plt.legend()plt.show()if __name__=='__main__':main()

作者：电气-余登武。原创不易，禁止转载