实验7 粒子群优化算法求解tsp问题[通俗易懂]

大家好，又见面了，我是你们的朋友全栈君。如果您正在找激活码,请点击查看最新教程,关注关注公众号 “全栈程序员社区” 获取激活教程,可能之前旧版本教程已经失效.最新Idea2022.1教程亲测有效,一键激活。

Jetbrains全系列IDE稳定放心使用

传送门（所有的实验都使用python实现）

实验1 BP神经网络实验

实验2 som网实验

实验3 hopfield实现八皇后问题

实验4 模糊搜索算法预测薄冰厚度

实验5 遗传算法求解tsp问题

实验6 蚁群算法求解tsp问题

实验7 粒子群优化算法求解tsp问题

实验8 分布估计算法求解背包问题

实验9 模拟退火算法求解背包问题

实验10 禁忌搜索算法求解tsp问题

一、实验目的

理解并使用粒子群优化算法

二、实验内容

实现基于粒子群优化算法的旅行商路线寻找

三、实验环境

使用Python3.0 在 eclipse进行编辑

四、实验步骤

1、输入介绍：

    城市总数目为14，采用以下城市坐标，坐标取整数。

2、初始解设定：

    给每一个粒子赋予随机解，和空交换序（即速度为0）。

3、计算每个粒子的下个位置：

（1）首先计算粒子当前位置与局部最优解的差，结果为一个交换序ss1，并以概率u1保留其中的交换子。同理计算粒子当前位置与全局最优解的差，以概率u2保存在交换序ss2。

（2）其次合并粒子当前速度speed，交换序ss1，交换序ss2三个交换序，以合并结果更新粒子速度

（3）最后将速度作用在粒子当前位置

4、计算粒子函数适应值：

    求出粒子函数适应值，并更新局部最优解与全局最优解

5.终止条件

    若全局最优不满足条件回到步骤3，否则打印结果，结束迭代。

运行截图：

路线随机选取 距离48.7

100个粒子迭代100次   距离 40.4

100个粒子迭代 600次   距离32.3

五、总结

    多次实验之后发现测试组数据的14个城市，所能达到的最优解gbest = 32.3；算法的效率受到粒子个数的影响十分明显。

通过表格可以看到随着粒子个数的减少,迭代次数逐渐增加，当粒子个数小于40时，迭代次数已经超过一万也还没有找到最优解。说明粒子的个数越多寻找的效率越高，但是每次迭代时间与粒子个数成正比，所以并非粒子越多越好，个数取50是最佳。

Python源码

#coding:gbk
import random
import math
import matplotlib.pyplot as plt
global n,m,u1,u2;    #n个粒子， m个城市 u1,u2为交换子保留概率
n=100;m=14;          
gbest = 10000;    gbestway= [0.0]*m;       #gbest记录全局最优适应值，gbestway记录全局最优解
pbestway =  [[0]*(m) for i in range(n)];   pbest = [0.0]*n;   #pbest记录个体最优适应值，gbestway记录个体最优解
road = [[0]*(m) for i in range(n)]   #开辟n*m的数组记录每个粒子当前位置
class no:                       #该类表示每个点的坐标
    def __init__(self,x,y):
        self.x=x;
        self.y=y;
p=[];
class so:                       #该类表示交换子
    def __init__(self,x,y):
        self.x=x;
        self.y=y;
ss=[so]*100;     global co;        #暂存减法结果            
speed = [[so]*(100) for i in range(n)]     #每个粒子的速度，即交换序
spco= [0.0]*n;                             #记录交换序的个数
def draw(t):              #该函数用于描绘路线图
    x=[0]*(m+1);y=[0]*(m+1);
    for i in range(m):
        x[i] =p[t[i]].x;
        y[i] =p[t[i]].y;
    x[m] =p[t[0]].x;
    y[m] =p[t[0]].y;
    plt.plot(x,y,color='r',marker='*' ); 
    plt.show();
def  mycol():                           #城市坐标输入
        p.append(no( 16 , 96 ));
        p.append(no( 16 , 94 )); p.append(no( 20 , 92 )); p.append(no( 22 , 93 )); p.append(no( 25 , 97 )); p.append(no( 22 , 96 )); p.append(no( 20 , 97 ));
        p.append(no( 17 , 96 )); p.append(no( 16 , 97 )); p.append(no( 14 , 98 )); p.append(no( 17, 97 )); p.append(no( 21 , 95 )); p.append(no( 19 , 97 ));
        p.append(no( 20 , 94 )); 

def get_value(t):        #计算粒子的函数适应值
    ans = 0.0;
    for i in range(1,m):     #两点距离公式
        ans += math.sqrt((p[t[i]].x-p[t[i-1]].x) *(p[t[i]].x-p[t[i-1]].x)  +(p[t[i]].y-p[t[i-1]].y) *(p[t[i]].y-p[t[i-1]].y));
    ans +=  math.sqrt((p[t[0]].x-p[t[m-1]].x) * (p[t[0]].x-p[t[m-1]].x)  +(p[t[0]].y-p[t[m-1]].y) *(p[t[0]].y-p[t[m-1]].y));
    return ans;
def cop(a,b,le):     #把b数组的值赋值a数组
    for i in range(le):
        a[i]=b[i]
def rand(g):                 # 随机解生成函数
    vis = [0]*m
    for i in range(m):
        vis[i]=0;
    on= 0
    while on<m:
        te = random.randint(0,m-1);
        if(vis[te]==0):
            vis[te]=1;
            g[on]=te;
            on+=1;
def find (g,ob):         #在数组g中寻找，值为ob的下标，并返回下标
    for i in range(m):
        if(g[i]==ob):
            return i;
def  cat(a, c,u):   #求解a减去解b的交换序结果  将其保存在ss序列 u为保留概率
    global co;
    co = 0;
    b = [0]*m;
    for i in range(m):
        b[i]=c[i]
    ob=0;
    for i in range(m):
        if(a[i]==b[i]):continue;
        ob=find(b,a[i])
        if(random.random() < u): 
            ss[co]=so(i,ob)
            co+=1
        b[i],b[ob]=b[ob],b[i]
def add(g,sv,le):                  #解g加上长度为le的交换序sv
    a=0;b=0;
    for i in range(le):
        a=sv[i].x;b=sv[i].y;
        g[a],g[b]=g[b], g[a];
        
def init():            #初始化函数
    global gbest
    for i in range(n):
        spco[i]=0;
        rand(road[i]);
        cop(pbestway[i],road[i],m); #将局部最优设为初始化的随机解
        pbest[i]= get_value(road[i]);
        if(pbest[i] < gbest):
            gbest = pbest[i];
            cop(gbestway,pbestway[i],m);
def update(i,r):   #更新i个体的函数适应值
    global gbest;
    te = get_value(r);        #计算适应值
    if(te < pbest[i]):         #个体最优更新
        pbest[i] =te;cop(pbestway[i],r,m);
    if(te < gbest):              #全局最优更新
            gbest = te;cop(gbestway,r,m);
    
def slove():    #执行函数
    global co,gbest,u1,u2
    t1 = [0]*m;t2=[0]*m;
    for i in range(m):
        t1[i]=i
        t2[i]=i;    
    u1 = 0.6;     #个体最优交换子保留概率
    u2 = 0.8     #全局最优交换子保留概率
    for i in range(n):
        for k in range(m):t1[k]=k;t2[k]=k; #构造两个解t1，t2求基本交换序
        add(t1,speed[i],spco[i]);
        cat(pbestway[i],road[i],u1);       #与个体最优解相减
        add(t1,ss,co);
        cat(gbestway,road[i],u2);            #与全局最优解相减
        add(t1,ss,co);
        cat(t1,t2,1)             #求出基本交换序
        for j in range(co): speed[i][j].x = ss[j].x;speed[i][j].y = ss[j].y    #更新个体速度;
        spco[i]=co;
        add(road[i],speed[i],spco[i]);  #将速度作用到当前位置
        update(i,road[i]);   #更新函数适应值
        
mycol();            #数据输入    
init();                 #数据初始化
for i in range(10000):      #控制迭代次数
    slove();
    if(gbest<32.5):  print('迭代次数',i); break;    #达到最优解提前退出
print(round(gbest,3));       #打印最优解距离保留三位小数
draw(gbestway);                      #画图描绘路线
print(gbestway);                     #打印路线，以序列表示