Python實現蟻群演演算法

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.3
    return 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] #初始化x
    X[i,1]=np.random.uniform(YMIN,YMAX,1)[0] #初始化y
    Tau[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 =XMIN
        if temp1 > XMAX:
            temp1 =XMAX
        if temp2 < XMIN:
            temp2 =XMIN
        if temp2 > XMAX:
            temp2 =XMAX
 
 
        #==判斷螞蟻是否移動(選更優===
        if func(temp1, temp2) < func(X[i, 0], X[i, 1]):
            X[i, 0] = temp1
            X[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()

最優變數x 5.0最優變數y 5.0
最優目標函數值 -123.7

#====================匯入相關庫=============================
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()