Python編程實現(xiàn)蟻群算法詳解

更新時間：2017年11月13日 17:19:59 作者：金良山莊

這篇文章主要介紹了Python編程實現(xiàn)蟻群算法詳解，涉及螞蟻算法的簡介，主要原理及公式，以及Python中的實現(xiàn)代碼，具有一定參考價值，需要的朋友可以了解下。

簡介

蟻群算法(ant colony optimization, ACO)，又稱螞蟻算法，是一種用來在圖中尋找優(yōu)化路徑的機率型算法。它由Marco Dorigo于1992年在他的博士論文中提出，其靈感來源于螞蟻在尋找食物過程中發(fā)現(xiàn)路徑的行為。蟻群算法是一種模擬進化算法，初步的研究表明該算法具有許多優(yōu)良的性質(zhì)。針對PID控制器參數(shù)優(yōu)化設(shè)計問題，將蟻群算法設(shè)計的結(jié)果與遺傳算法設(shè)計的結(jié)果進行了比較，數(shù)值仿真結(jié)果表明，蟻群算法具有一種新的模擬進化優(yōu)化方法的有效性和應用價值。
定義

各個螞蟻在沒有事先告訴他們食物在什么地方的前提下開始尋找食物。當一只找到食物以后，它會向環(huán)境釋放一種揮發(fā)性分泌物pheromone (稱為信息素,該物質(zhì)隨著時間的推移會逐漸揮發(fā)消失，信息素濃度的大小表征路徑的遠近)來實現(xiàn)的，吸引其他的螞蟻過來，這樣越來越多的螞蟻會找到食物。有些螞蟻并沒有像其它螞蟻一樣總重復同樣的路，他們會另辟蹊徑，如果另開辟的道路比原來的其他道路更短，那么，漸漸地，更多的螞蟻被吸引到這條較短的路上來。最后，經(jīng)過一段時間運行，可能會出現(xiàn)一條最短的路徑被大多數(shù)螞蟻重復著。

解決的問題

三維地形中，給出起點和重點，找到其最優(yōu)路徑。

作圖源碼：

from mpl_toolkits.mplot3d import proj3d
from mpl_toolkits.mplot3d import Axes3D
import numpy as np

height3d = np.array([[2000,1400,800,650,500,750,1000,950,900,800,700,900,1100,1050,1000,1150,1300,1250,1200,1350,1500],          [1100,900,700,625,550,825,1100,1150,1200,925,650,750,850,950,1050,1175,1300,1350,1400,1425,1450],          [200,400,600,600,600,900,1200,1350,1500,1050,600,600,600,850,1100,1200,1300,1450,1600,1500,1400],          [450,500,550,575,600,725,850,875,900,750,600,600,600,725,850,900,950,1150,1350,1400,1450],          [700,600,500,550,600,550,500,400,300,450,600,600,600,600,600,600,600,850,1100,1300,1500],          [500,525,550,575,600,575,550,450,350,475,600,650,700,650,600,600,600,725,850,1150,1450],          [300,450,600,600,600,600,600,500,400,500,600,700,800,700,600,600,600,600,600,1000,1400],          [550,525,500,550,600,875,1150,900,650,725,800,700,600,875,1150,1175,1200,975,750,875,1000],          [800,600,400,500,600,1150,1700,1300,900,950,1000,700,400,1050,1700,1750,1800,1350,900,750,600],          [650,600,550,625,700,1175,1650,1275,900,1100,1300,1275,1250,1475,1700,1525,1350,1200,1050,950,850],          [500,600,700,750,800,1200,1600,1250,900,1250,1600,1850,2100,1900,1700,1300,900,1050,1200,1150,1100],          [400,375,350,600,850,1200,1550,1250,950,1225,1500,1750,2000,1950,1900,1475,1050,975,900,1175,1450],          [300,150,0,450,900,1200,1500,1250,1000,1200,1400,1650,1900,2000,2100,1650,1200,900,600,1200,1800],          [600,575,550,750,950,1275,1600,1450,1300,1300,1300,1525,1750,1625,1500,1450,1400,1125,850,1200,1550],          [900,1000,1100,1050,1000,1350,1700,1650,1600,1400,1200,1400,1600,1250,900,1250,1600,1350,1100,1200,1300],          [750,850,950,900,850,1000,1150,1175,1200,1300,1400,1325,1250,1125,1000,1150,1300,1075,850,975,1100],          [600,700,800,750,700,650,600,700,800,1200,1600,1250,900,1000,1100,1050,1000,800,600,750,900],          [750,775,800,725,650,700,750,775,800,1000,1200,1025,850,975,1100,950,800,900,1000,1050,1100],          [900,850,800,700,600,750,900,850,800,800,800,800,800,950,1100,850,600,1000,1400,1350,1300],          [750,800,850,850,850,850,850,825,800,750,700,775,850,1000,1150,875,600,925,1250,1100,950],          [600,750,900,1000,1100,950,800,800,800,700,600,750,900,1050,1200,900,600,850,1100,850,600]])

fig = figure()
ax = Axes3D(fig)
X = np.arange(21)
Y = np.arange(21)
X, Y = np.meshgrid(X, Y)
Z = -20*np.exp(-0.2*np.sqrt(np.sqrt(((X-10)**2+(Y-10)**2)/2)))+20+np.e-np.exp((np.cos(2*np.pi*X)+np.sin(2*np.pi*Y))/2)
ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap='cool')
ax.set_xlabel('X axis')
ax.set_ylabel('Y axis')
ax.set_zlabel('Z')
ax.set_title('3D map')


point0 = [0,9,Z[0][9]] 
point1 = [20,7,Z[20][7]]

ax.plot([point0[0]],[point0[1]],[point0[2]],'r',marker = u'o',markersize = 15)
ax.plot([point1[0]],[point1[1]],[point1[2]],'r',marker = u'o',markersize = 15)

x0,y0,_ = proj3d.proj_transform(point0[0],point0[1],point0[2], ax.get_proj())
x1,y1,_ = proj3d.proj_transform(point1[0],point1[1],point1[2], ax.get_proj())

label = pylab.annotate(
  "start", 
  xy = (x0, y0), xytext = (-20, 20),
  textcoords = 'offset points', ha = 'right', va = 'bottom',
  bbox = dict(boxstyle = 'round,pad=0.5', fc = 'yellow', alpha = 1),
  arrowprops = dict(arrowstyle = '->', connectionstyle = 'arc3,rad=0'),fontsize=15)
label2 = pylab.annotate(
  "end", 
  xy = (x1, y1), xytext = (-20, 20),
  textcoords = 'offset points', ha = 'right', va = 'bottom',
  bbox = dict(boxstyle = 'round,pad=0.5', fc = 'yellow', alpha = 1),
  arrowprops = dict(arrowstyle = '->', connectionstyle = 'arc3,rad=0'),fontsize=15)
def update_position(e):
  x2, y2, _ = proj3d.proj_transform(point0[0],point0[1],point0[2],ax.get_proj())
  label.xy = x2,y2
  label.update_positions(fig.canvas.renderer)

  x1,y1,_ = proj3d.proj_transform(point1[0],point1[1],point1[2],ax.get_proj())
  label2.xy = x1,y1
  label2.update_positions(fig.canvas.renderer)
  fig.canvas.draw()

fig.canvas.mpl_connect('button_release_event', update_position)

基本原理

螞蟻k根據(jù)各個城市間鏈接路徑上的信息素濃度決定其下一個訪問城市，設(shè)P^k_ij(t)表示t時刻螞蟻k從城市i轉(zhuǎn)移到矩陣j的概率，其計算公式為

計算完城市間的轉(zhuǎn)移概率后，采用與遺傳算法中一樣的輪盤賭方法選擇下一個待訪問的城市。

當所有的螞蟻完成一次循環(huán)后，各個城市間鏈接路徑上的信息素濃度需進行更新，計算公式為

其中，Δτkij表示第k只螞蟻在城市i與城市j連接路徑上釋放的信息素濃度；Δτij表示所有螞蟻在城市i與城市j連接路徑上釋放的信息素濃度之和。

螞蟻釋放信息素的模型

程序代碼：

import numpy as np
import matplotlib.pyplot as plt
%pylab
coordinates = np.array([[565.0,575.0],[25.0,185.0],[345.0,750.0],[945.0,685.0],[845.0,655.0],
            [880.0,660.0],[25.0,230.0],[525.0,1000.0],[580.0,1175.0],[650.0,1130.0],
            [1605.0,620.0],[1220.0,580.0],[1465.0,200.0],[1530.0, 5.0],[845.0,680.0],
            [725.0,370.0],[145.0,665.0],[415.0,635.0],[510.0,875.0],[560.0,365.0],
            [300.0,465.0],[520.0,585.0],[480.0,415.0],[835.0,625.0],[975.0,580.0],
            [1215.0,245.0],[1320.0,315.0],[1250.0,400.0],[660.0,180.0],[410.0,250.0],
            [420.0,555.0],[575.0,665.0],[1150.0,1160.0],[700.0,580.0],[685.0,595.0],
            [685.0,610.0],[770.0,610.0],[795.0,645.0],[720.0,635.0],[760.0,650.0],
            [475.0,960.0],[95.0,260.0],[875.0,920.0],[700.0,500.0],[555.0,815.0],
            [830.0,485.0],[1170.0, 65.0],[830.0,610.0],[605.0,625.0],[595.0,360.0],
            [1340.0,725.0],[1740.0,245.0]])
def getdistmat(coordinates):
  num = coordinates.shape[0]
  distmat = np.zeros((52,52))
  for i in range(num):
    for j in range(i,num):
      distmat[i][j] = distmat[j][i]=np.linalg.norm(coordinates[i]-coordinates[j])
  return distmat
distmat = getdistmat(coordinates)
numant = 40 #螞蟻個數(shù)
numcity = coordinates.shape[0] #城市個數(shù)
alpha = 1  #信息素重要程度因子
beta = 5  #啟發(fā)函數(shù)重要程度因子
rho = 0.1  #信息素的揮發(fā)速度
Q = 1
iter = 0
itermax = 250
etatable = 1.0/(distmat+np.diag([1e10]*numcity)) #啟發(fā)函數(shù)矩陣，表示螞蟻從城市i轉(zhuǎn)移到矩陣j的期望程度
pheromonetable = np.ones((numcity,numcity)) # 信息素矩陣
pathtable = np.zeros((numant,numcity)).astype(int) #路徑記錄表
distmat = getdistmat(coordinates) #城市的距離矩陣
lengthaver = np.zeros(itermax) #各代路徑的平均長度
lengthbest = np.zeros(itermax) #各代及其之前遇到的最佳路徑長度
pathbest = np.zeros((itermax,numcity)) # 各代及其之前遇到的最佳路徑長度
while iter < itermax:
  # 隨機產(chǎn)生各個螞蟻的起點城市
  if numant <= numcity:#城市數(shù)比螞蟻數(shù)多
    pathtable[:,0] = np.random.permutation(range(0,numcity))[:numant]
  else: #螞蟻數(shù)比城市數(shù)多，需要補足
    pathtable[:numcity,0] = np.random.permutation(range(0,numcity))[:]
    pathtable[numcity:,0] = np.random.permutation(range(0,numcity))[:numant-numcity]
  length = np.zeros(numant) #計算各個螞蟻的路徑距離
  for i in range(numant):
    visiting = pathtable[i,0] # 當前所在的城市
    #visited = set() #已訪問過的城市，防止重復
    #visited.add(visiting) #增加元素
    unvisited = set(range(numcity))#未訪問的城市
    unvisited.remove(visiting) #刪除元素
    for j in range(1,numcity):#循環(huán)numcity-1次，訪問剩余的numcity-1個城市
      #每次用輪盤法選擇下一個要訪問的城市
      listunvisited = list(unvisited)
      probtrans = np.zeros(len(listunvisited))
      for k in range(len(listunvisited)):
        probtrans[k] = np.power(pheromonetable[visiting][listunvisited[k]],alpha)\
            *np.power(etatable[visiting][listunvisited[k]],alpha)
      cumsumprobtrans = (probtrans/sum(probtrans)).cumsum()
      cumsumprobtrans -= np.random.rand()
      k = listunvisited[find(cumsumprobtrans>0)[0]] #下一個要訪問的城市
      pathtable[i,j] = k
      unvisited.remove(k)
      #visited.add(k)
      length[i] += distmat[visiting][k]
      visiting = k
    length[i] += distmat[visiting][pathtable[i,0]] #螞蟻的路徑距離包括最后一個城市和第一個城市的距離
  #print length
  # 包含所有螞蟻的一個迭代結(jié)束后，統(tǒng)計本次迭代的若干統(tǒng)計參數(shù)
  lengthaver[iter] = length.mean()
  if iter == 0:
    lengthbest[iter] = length.min()
    pathbest[iter] = pathtable[length.argmin()].copy()   
  else:
    if length.min() > lengthbest[iter-1]:
      lengthbest[iter] = lengthbest[iter-1]
      pathbest[iter] = pathbest[iter-1].copy()
    else:
      lengthbest[iter] = length.min()
      pathbest[iter] = pathtable[length.argmin()].copy()  
  # 更新信息素
  changepheromonetable = np.zeros((numcity,numcity))
  for i in range(numant):
    for j in range(numcity-1):
      changepheromonetable[pathtable[i,j]][pathtable[i,j+1]] += Q/distmat[pathtable[i,j]][pathtable[i,j+1]]
    changepheromonetable[pathtable[i,j+1]][pathtable[i,0]] += Q/distmat[pathtable[i,j+1]][pathtable[i,0]]
  pheromonetable = (1-rho)*pheromonetable + changepheromonetable
  iter += 1 #迭代次數(shù)指示器+1
  #觀察程序執(zhí)行進度，該功能是非必須的
  if (iter-1)%20==0: 
    print iter-1
# 做出平均路徑長度和最優(yōu)路徑長度    
fig,axes = plt.subplots(nrows=2,ncols=1,figsize=(12,10))
axes[0].plot(lengthaver,'k',marker = u'')
axes[0].set_title('Average Length')
axes[0].set_xlabel(u'iteration')
axes[1].plot(lengthbest,'k',marker = u'')
axes[1].set_title('Best Length')
axes[1].set_xlabel(u'iteration')
fig.savefig('Average_Best.png',dpi=500,bbox_inches='tight')
plt.close()
#作出找到的最優(yōu)路徑圖
bestpath = pathbest[-1]
plt.plot(coordinates[:,0],coordinates[:,1],'r.',marker=u'$\cdot$')
plt.xlim([-100,2000])
plt.ylim([-100,1500])
for i in range(numcity-1):#
  m,n = bestpath[i],bestpath[i+1]
  print m,n
  plt.plot([coordinates[m][0],coordinates[n][0]],[coordinates[m][1],coordinates[n][1]],'k')
plt.plot([coordinates[bestpath[0]][0],coordinates[n][0]],[coordinates[bestpath[0]][1],coordinates[n][1]],'b')
ax=plt.gca()
ax.set_title("Best Path")
ax.set_xlabel('X axis')
ax.set_ylabel('Y_axis')
plt.savefig('Best Path.png',dpi=500,bbox_inches='tight')
plt.close()

總結(jié)

以上就是本文關(guān)于Python編程實現(xiàn)蟻群算法詳解的全部內(nèi)容，希望對大家有所幫助。感興趣的朋友可以繼續(xù)參閱本站：python實現(xiàn)圖片處理和特征提取詳解、python圖像常規(guī)操作、python先序遍歷二叉樹問題等，有什么問題可以隨時留言，小編會及時回復大家的。感謝朋友們對本站的支持！

您可能感興趣的文章: