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 = 10
    pi = np.pi
    x = 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] - 6
    ee += max(0, e1)
    """计算第二个约束的惩罚项"""
    e2 = 3 * X[0] - 2 * X[1] - 5
    ee += 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_e
        else:
            return child,child_fitness,child_e
    # 规则2，如果child违反约束而parent没有违反约束，则取parent
    if parent_e < 0.00001 and child_e  >= 0.00001:
        return parent,parent_fitness,parent_e
    # 规则3，如果parent违反约束而child没有违反约束，则取child
    if 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_e
    else:
        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]  # 初始化x
        X[i, 1] =np.random.uniform(YMIN, YMAX, 1)[0]  # 初始化y
        Tau[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]  # 初始化x
        if (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]=pbesti
    return 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 Tau
 
#=============主函数======================
def main():
    X,Tau=initialization() #初始化蚂蚁群X 和信息素 Tau
    for 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]  # 最优变量x
    minY = X[min_index, 1]  # 最优变量y
    minValue = 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()