Geatpy求解TSP问题

import geatpy as ea
import numpy as np
from math import sin, asin, cos, radians, fabs, sqrt
import matplotlib.pyplot as plt

# TSP坐标点
places=[[117.000923, 36.675807], [115.48333, 38.03333], [125.35, 43.88333], [127.63333, 47.75], [123.38333, 41.8], 
        [111.670801, 41.818311], [87.68333, 43.76667], [103.73333, 36.03333], [106.26667, 37.46667], [112.53333, 37.86667], 
        [108.95, 34.26667], [113.65, 34.76667], [117.283042, 31.86119], [119.78333, 32.05], [120.2, 30.26667], [118.3, 26.08333],
        [113.23333, 23.16667], [115.9, 28.68333], [110.35, 20.01667], [108.320004, 22.82402], [106.71667, 26.56667], [113.0, 28.21667],
        [114.298572, 30.584355], [104.06667, 30.66667], [102.73333, 25.05], [91.0, 30.6], [96.75, 36.56667], [117.2, 39.13333], 
        [121.55333, 31.2], [106.45, 29.56667], [116.41667, 39.91667], [121.3, 25.03], [114.1, 22.2], [113.5, 22.2]]
#print(len(places))
def get_distance(position1,position2):
    '''
        获取相邻两个城市的距离
        position1:第一个地区的经纬度，为列表形式
    
        position2:第二个地区的经纬度，为列表形式
    
        传入两个地区的地理坐标，返回二者的距离
    '''
    lng1,lat1,lng2,lat2 = (position1[0],position1[1],position2[0],position2[1])
    lng1, lat1, lng2, lat2 = map(radians, [float(lng1), float(lat1), float(lng2), float(lat2)]) # 经纬度转换成弧度
    dlon=lng2-lng1
    dlat=lat2-lat1
    a=sin(dlat/2)**2 + cos(lat1) * cos(lat2) * sin(dlon/2)**2
    distance=2*asin(sqrt(a))*6371*1000 # 地球平均半径，6371km
    distance=round(distance/1000,3)
    return distance

# 定义问题
class MyProblem(ea.Problem):  # 继承Problem父类
    def __init__(self):
        name = 'MyProblem'  # 初始化name（函数名称，可以随意设置）
        M = 1  # 初始化M（目标维数）
        maxormins = [1]  # 初始化maxormins（目标最小最大化标记列表，1：最小化该目标；-1：最大化该目标）
        Dim = len(places) - 1 # 初始化Dim（决策变量维数）
        varTypes = [1] * Dim  # 初始化varTypes（决策变量的类型，元素为0表示对应的变量是连续的；1表示是离散的）
        lb = [1] * Dim  # 决策变量下界
        ub = [Dim] * Dim  # 决策变量上界
        lbin = [1] * Dim  # 决策变量下边界（0表示不包含该变量的下边界，1表示包含）
        ubin = [1] * Dim  # 决策变量上边界（0表示不包含该变量的上边界，1表示包含）
        # 调用父类构造方法完成实例化
        ea.Problem.__init__(self, name, M, maxormins, Dim, varTypes, lb, ub, lbin, ubin)
        # 新增一个属性存储旅行地坐标
        self.places = np.array(places)

    def aimFunc(self, pop):  # 目标函数
        x = pop.Phen.copy()  # 得到决策变量矩阵
        # 添加最后回到出发地
        X = np.hstack([x, x[:, [0]]]).astype(int)
        ObjV = []  # 存储所有种群个体对应的总路程
        for i in range(pop.sizes):
            journey = self.places[X[i], :]  # 按既定顺序到达的地点坐标
            distance = 0.0
            for j in range(len(journey)-1):
                distance+=get_distance(journey[j],journey[j+1])
            #distance = np.sum(np.sqrt(np.sum(np.diff(journey.T) ** 2, 0)))  # 计算总路程
            #print("dis",distance)
            ObjV.append(distance)
        pop.ObjV = np.array([ObjV]).T


"""调用模板求解"""
if __name__ == '__main__':
    """================================实例化问题对象============================"""
    problem = MyProblem()  # 生成问题对象
    """==================================种群设置=============================="""
    Encoding = 'P'  # 编码方式，采用排列编码
    NIND = 100  # 种群规模
    Field = ea.crtfld(Encoding, problem.varTypes, problem.ranges, problem.borders)  # 创建区域描述器
    population = ea.Population(Encoding, Field, NIND)  # 实例化种群对象（此时种群还没被初始化，仅仅是完成种群对象的实例化）
    """================================算法参数设置============================="""
    myAlgorithm = ea.soea_SEGA_templet(problem, population)  # 实例化一个算法模板对象
    myAlgorithm.MAXGEN = 500 # 最大进化代数
    myAlgorithm.mutOper.Pm = 0.8  # 变异概率
    myAlgorithm.logTras = 0  # 设置每隔多少代记录日志，若设置成0则表示不记录日志
    myAlgorithm.verbose = True  # 设置是否打印输出日志信息
    myAlgorithm.drawing = 1  # 设置绘图方式（0：不绘图；1：绘制结果图；2：绘制目标空间过程动画；3：绘制决策空间过程动画）
    """===========================调用算法模板进行种群进化========================"""
    [BestIndi, population] = myAlgorithm.run()  # 执行算法模板，得到最优个体以及最后一代种群
    BestIndi.save()  # 把最优个体的信息保存到文件中
    """==================================输出结果=============================="""
    print('评价次数：%s' % myAlgorithm.evalsNum)
    print('时间已过 %s 秒' % myAlgorithm.passTime)
    if BestIndi.sizes != 0:
        print('最短路程为：%s' % BestIndi.ObjV[0][0])
        print('最佳路线为：')
        best_journey = np.hstack([0, BestIndi.Phen[0, :], 0])
        for i in range(len(best_journey)):
            print(int(best_journey[i]), end=' ')
        print()
        # 绘图
        plt.figure()
        plt.plot(problem.places[best_journey.astype(int), 0], problem.places[best_journey.astype(int), 1], c='black')
        plt.plot(problem.places[best_journey.astype(int), 0], problem.places[best_journey.astype(int), 1], 'o',
                 c='black')
        for i in range(len(best_journey)):
            plt.text(problem.places[int(best_journey[i]), 0], problem.places[int(best_journey[i]), 1],
                     chr(int(best_journey[i]) + 65), fontsize=20)
        plt.grid(True)
        plt.xlabel('x坐标')
        plt.ylabel('y坐标')
        plt.savefig('roadmap.svg', dpi=600, bbox_inches='tight')
    else:
        print('没找到可行解。')

运行结果：

import folium
path = []
for i in best_journey:
    path.append(places[i])
m1 = folium.Map(location=[34.26667,108.95000], width=1000, height=500,zoom_start=4,
               #tiles='http://webst02.is.autonavi.com/appmaptile?style=6&x={x}&y={y}&z={z}', # 高德卫星图
               #tiles='http://webrd02.is.autonavi.com/appmaptile?lang=zh_cn&size=1&scale=1&style=7&x={x}&y={y}&z={z}', # 高德街道图
               #attr='default'
               tiles='OpenStreetMap')
for i in range(len(path)-1):
    begin = path[i][::-1]
    end = path[i+1][::-1]
    coordinates=[begin,end] 
    aline=folium.PolyLine(locations=coordinates,weight=2,color = 'red') 
    m1.add_child(aline) 
begin = end
end = path[0][::-1]
coordinates=[begin,end] 
aline=folium.PolyLine(locations=coordinates,weight=2,color = 'red') 
m1.add_child(aline)
m1

运行结果: