Python求解CVRP问题

一.问题描述

考虑带容量限制的VRP问题，即CVRP。

1.1 假设

• 配送中心：只有一个配送中心；
• 配送方式：任何车辆配送货物结束后需要返回配送中心；
• 车型：只考虑一种车型，且容量约束相同；
• 需求不可拆分：客户需求只能由一辆车满足，即任一客户需求小于车辆容量；
• 车辆充足：不限制车辆数量，即配送车辆需求均能满足；

1.2 要求

• 优化目标：最小化车辆启动成本和车辆行驶成本之和；
• 约束条件：客户访问约束，容积约束；
• 已知信息：配送中心位置、客户点位置、客户点需求、车辆最大容积、车辆启动成本、车辆单位距离行驶成本；（数据如有需要请后台私信我）

二.遗传算法

遗传算法是一种模仿自然界生物进化机制发展的随机全局搜索和优化方法，本质是一种高效、并行、全局搜索的方法，能在搜索过程中自动获取和积累有关搜索空间的知识，自适应地控制搜索过程以求得近优解或最优解。

2.1 基本概念

• 个体（individual）：指染色体带有特征的实体；
• 种群（population）：个体的集合，该集合内个体数称为种群的大小；
• 编码(coding)：DNA中遗传信息在一个长链上按一定的模式排列，遗传编码可看作从表现型到基因型的映射；
• 适应度(fitness)：度量某个物种对于生存环境的适应程度；
• 选择(selection)：以一定的概率从种群中选择若干个个体，一般选择过程是一种基于适应度的轮盘赌过程；
• 解码(decoding)：基因型到表现型的映射。
• 交叉(crossover)：两个染色体的某一相同位置处DNA被切断，前后两串分别交叉组合形成两个新的染色体，也称基因重组或杂交；
• 变异(mutation)：复制时可能（很小的概率）产生某些复制差错，变异产生新的染色体，表现出新的性状；

2.2 算法流程

开始循环直至找到满意的解。

1. 评估每条染色体所对应个体的适应度；
2. 遵照适应度越高，选择概率越大的原则，从种群中选择两个个体作为父方和母方；
3. 抽取父母双方的染色体，进行染色体交叉，产生子代；
4. 对子代的染色体进行变异。
5. 重复2，3，4步骤，直到新种群的产生。

达到终止条件结束循环。

三.Python代码实现

step.0 导入相应模块

import pandas as pd
import math
import random
import numpy as np
import copy
import xlsxwriter
import matplotlib.pyplot as plt

step.1 定义解、网络节点、全局参数的类

#数据结构：网络节点
class Sol():
    def __init__(self):
        self.nodes_seq=None # 解的编码
        self.obj=None # 目标函数
        self.fit=None  # 适应度
        self.routes=None # 解的解码
# 数据结构：网络节点
class Node():
    def __init__(self):
        self.id=0 # 节点id
        self.name=''  # 节点名称，可选
        self.seq_no=0 # 节点映射id
        self.x_coord=0 # 节点平面横坐标
        self.y_coord=0  # 节点平面纵坐标
        self.demand=0 # 节点需求
# 数据结构：全局参数
class Model():
    def __init__(self):
        self.best_sol=None # 全局最优解
        self.node_list=[]  # 需求节点集合
        self.sol_list=[]  # 解的集合
        self.node_seq_no_list=[]  #需求节点映射id集合
        self.depot=None # 车场节点
        self.number_of_nodes=0 # 需求节点数量
        self.opt_type=0 # 优化目标类型
        self.vehicle_cap=0 # 车辆最大容量
        self.pc=0.5 # 交叉概率
        self.pm=0.2 # 变异概率
        self.n_select=80 # 种群选择数量
        self.popsize=100 # 种群规模

step.2 处理表格数据

def readXlsxFile(filepath,model):
    # 建议在xlsx文件中，第一行为表头，其中: x_coord,y_coord,demand是必须项；车辆基地数据放在表头下首行
    node_seq_no =-1 #车辆基地的seq_no值为-1,剩余需求节点的seq_no 依次编号为 0,1,2,...
    df = pd.read_excel('file_path')
    for i in range(df.shape[0]):
        node=Node()
        node.id=node_seq_no
        node.seq_no=node_seq_no
        node.x_coord= df['x_coord'][i]
        node.y_coord= df['y_coord'][i]
        node.demand=df['demand'][i]
        if df['demand'][i] == 0:
            model.depot=node
        else:
            model.node_list.append(node)
            model.node_seq_no_list.append(node_seq_no)
        try:
            node.name=df['name'][i]
        except:
            pass
        try:
            node.id=df['id'][i]
        except:
            pass
        node_seq_no=node_seq_no+1
    model.number_of_nodes=len(model.node_list)

step.3 构造初始解

def genInitialSol(model):
    nodes_seq=copy.deepcopy(model.node_seq_no_list)
    for i in range(model.popsize):
        seed=int(random.randint(0,10))
        random.seed(seed)
        random.shuffle(nodes_seq)
        sol=Sol()
        sol.nodes_seq=copy.deepcopy(nodes_seq)
        model.sol_list.append(sol)

采用Split思想对TSP序列进行切割，得到可行车辆路径

def splitRoutes(nodes_seq,model):
    """
    采用简单的分割方法：按顺序依次检查路径的容量约束，在超出车辆容量限制的位置插入车场。
    例如某TSP解为：[1,2,3,4,5,6,7,8,9,10],累计需求为：[10,20,30,40,50,60,70,80,90,10]，车辆容量为：30，则应在3,6,9节点后插入车场，
    即得到：[0,1,2,3,0,4,5,6,0,7,8,9,0,10,0]
    """
    num_vehicle = 0
    vehicle_routes = []
    route = []
    remained_cap = model.vehicle_cap
    for node_no in nodes_seq:
        if remained_cap - model.node_list[node_no].demand >= 0:
            route.append(node_no)
            remained_cap = remained_cap - model.node_list[node_no].demand
        else:
            vehicle_routes.append(route)
            route = [node_no]
            num_vehicle = num_vehicle + 1
            remained_cap =model.vehicle_cap - model.node_list[node_no].demand
    vehicle_routes.append(route)
    return num_vehicle,vehicle_routes

计算所切割车辆路径的总行驶距离

def calDistance(route,model):
    distance=0
    depot=model.depot
    for i in range(len(route)-1):
        from_node=model.node_list[route[i]]
        to_node=model.node_list[route[i+1]]
        distance+=math.sqrt((from_node.x_coord-to_node.x_coord)**2+(from_node.y_coord-to_node.y_coord)**2)
    first_node=model.node_list[route[0]]
    last_node=model.node_list[route[-1]]
    distance+=math.sqrt((depot.x_coord-first_node.x_coord)**2+(depot.y_coord-first_node.y_coord)**2)
    distance+=math.sqrt((depot.x_coord-last_node.x_coord)**2+(depot.y_coord - last_node.y_coord)**2)
    return distance

step.4 计算适应度

def calFit(model):
    #calculate fit value：fit=Objmax-obj
    Objmax=-float('inf')
    best_sol=Sol()#record the local best solution
    best_sol.obj=float('inf')

    for sol in model.sol_list:
        nodes_seq=sol.nodes_seq
        num_vehicle, vehicle_routes = splitRoutes(nodes_seq, model)
        if model.opt_type==0:
            sol.obj=num_vehicle
            sol.routes=vehicle_routes
            if sol.obj>Objmax:
                Objmax=sol.obj
            if sol.obj<best_sol.obj:
                best_sol=copy.deepcopy(sol)
        else:
            distance=0
            for route in vehicle_routes:
                distance+=calDistance(route,model)
            sol.obj=distance
            sol.routes=vehicle_routes
            if sol.obj>Objmax:
                Objmax=sol.obj
            if sol.obj < best_sol.obj:
                best_sol = copy.deepcopy(sol)
    #calculate fit value
    for sol in model.sol_list:
        sol.fit=Objmax-sol.obj
    #update the global best solution
    if best_sol.obj<model.best_sol.obj:
        model.best_sol=best_sol

step.5 选择父代

def selectSol(model):
    sol_list=copy.deepcopy(model.sol_list)
    model.sol_list=[]
    for i in range(model.n_select):
        f1_index=random.randint(0,len(sol_list)-1)
        f2_index=random.randint(0,len(sol_list)-1)
        f1_fit=sol_list[f1_index].fit
        f2_fit=sol_list[f2_index].fit
        if f1_fit<f2_fit:
            model.sol_list.append(sol_list[f2_index])
        else:
            model.sol_list.append(sol_list[f1_index])

step.6 对父代染色体进行交叉

def crossSol(model):
    sol_list=copy.deepcopy(model.sol_list)
    model.sol_list=[]
    while True:
        f1_index = random.randint(0, len(sol_list) - 1)
        f2_index = random.randint(0, len(sol_list) - 1)
        if f1_index!=f2_index:
            f1 = copy.deepcopy(sol_list[f1_index])
            f2 = copy.deepcopy(sol_list[f2_index])
            if random.random() <= model.pc:
                cro1_index=int(random.randint(0,model.number_of_nodes-1))
                cro2_index=int(random.randint(cro1_index,model.number_of_nodes-1))
                new_c1_f = []
                new_c1_m=f1.nodes_seq[cro1_index:cro2_index+1]
                new_c1_b = []
                new_c2_f = []
                new_c2_m=f2.nodes_seq[cro1_index:cro2_index+1]
                new_c2_b = []
                for index in range(model.number_of_nodes):
                    if len(new_c1_f)<cro1_index:
                        if f2.nodes_seq[index] not in new_c1_m:
                            new_c1_f.append(f2.nodes_seq[index])
                    else:
                        if f2.nodes_seq[index] not in new_c1_m:
                            new_c1_b.append(f2.nodes_seq[index])
                for index in range(model.number_of_nodes):
                    if len(new_c2_f)<cro1_index:
                        if f1.nodes_seq[index] not in new_c2_m:
                            new_c2_f.append(f1.nodes_seq[index])
                    else:
                        if f1.nodes_seq[index] not in new_c2_m:
                            new_c2_b.append(f1.nodes_seq[index])
                new_c1=copy.deepcopy(new_c1_f)
                new_c1.extend(new_c1_m)
                new_c1.extend(new_c1_b)
                f1.nodes_seq=new_c1
                new_c2=copy.deepcopy(new_c2_f)
                new_c2.extend(new_c2_m)
                new_c2.extend(new_c2_b)
                f2.nodes_seq=new_c2
                model.sol_list.append(copy.deepcopy(f1))
                model.sol_list.append(copy.deepcopy(f2))
            else:
                model.sol_list.append(copy.deepcopy(f1))
                model.sol_list.append(copy.deepcopy(f2))
            if len(model.sol_list)>model.popsize:
                break

step.7 子代染色体变异

def muSol(model):
    sol_list=copy.deepcopy(model.sol_list)
    model.sol_list=[]
    while True:
        f1_index = int(random.randint(0, len(sol_list) - 1))
        f1 = copy.deepcopy(sol_list[f1_index])
        m1_index=random.randint(0,model.number_of_nodes-1)
        m2_index=random.randint(0,model.number_of_nodes-1)
        if m1_index!=m2_index:
            if random.random() <= model.pm:
                node1=f1.nodes_seq[m1_index]
                f1.nodes_seq[m1_index]=f1.nodes_seq[m2_index]
                f1.nodes_seq[m2_index]=node1
                model.sol_list.append(copy.deepcopy(f1))
            else:
                model.sol_list.append(copy.deepcopy(f1))
            if len(model.sol_list)>model.popsize:
                break

step.8 输出结果

绘制目标函数收敛曲线

def plotObj(obj_list):
    plt.rcParams['font.sans-serif'] = ['SimHei'] #show chinese
    plt.rcParams['axes.unicode_minus'] = False  # Show minus sign
    plt.plot(np.arange(1,len(obj_list)+1),obj_list)
    plt.xlabel('Iterations')
    plt.ylabel('Obj Value')
    plt.grid()
    plt.xlim(1,len(obj_list)+1)
    plt.show()

绘制优化车辆路径

def plotRoutes(model):
    for route in model.best_sol.routes:
        x_coord = [model.depot.x_coord]
        y_coord = [model.depot.y_coord]
        for node_no in route:
            x_coord.append(model.node_list[node_no].x_coord)
            y_coord.append(model.node_list[node_no].y_coord)
        x_coord.append(model.depot.x_coord)
        y_coord.append(model.depot.y_coord)
        plt.plot(x_coord, y_coord, marker='s', color='b', linewidth=0.5)
    plt.show()

输出优化结果

def outPut(model):
    work=xlsxwriter.Workbook('result.xlsx')
    worksheet=work.add_worksheet()
    worksheet.write(0,0,'opt_type')
    worksheet.write(1,0,'obj')
    if model.opt_type==0:
        worksheet.write(0,1,'number of vehicles')
    else:
        worksheet.write(0, 1, 'drive distance of vehicles')
    worksheet.write(1,1,model.best_sol.obj)
    for row,route in enumerate(model.best_sol.routes):
        route.insert(0,model.depot.id)
        route.append(model.depot.id)
        worksheet.write(row+2,0,'v'+str(row+1))
        r=[str(i)for i in route]
        worksheet.write(row+2,1, '-'.join(r))
    work.close()

step.9 模型参数赋值

def run(filepath,epochs,pc,pm,popsize,n_select,v_cap,opt_type):
    """
    :param filepath: Xlsx文件路径
    :param epochs: 迭代次数
    :param pc: 交叉概率
    :param pm: 变异概率
    :param popsize: 种群规模
    :param n_select: 优秀父代保留数量
    :param v_cap: 车辆容量
    :param opt_type: 优化目标，0：最小化车辆数，1：最小化总行驶距离
    :return:
    """
    model=Model()
    model.vehicle_cap=v_cap
    model.opt_type=opt_type
    model.pc=pc
    model.pm=pm
    model.popsize=popsize
    model.n_select=n_select

    readXlsxFile(filepath,model)
    genInitialSol(model)
    history_best_obj = []
    best_sol=Sol()
    best_sol.obj=float('inf')
    model.best_sol=best_sol
    for ep in range(epochs):
        calFit(model)
        selectSol(model)
        crossSol(model)
        muSol(model)
        history_best_obj.append(model.best_sol.obj)
        print("%s/%s， best obj: %s" % (ep,epochs,model.best_sol.obj))
    plotObj(history_best_obj)
    plotRoutes(model)
    outPut(model)
if __name__=='__main__':
     file=r'filepath'
     run(filepath=file,epochs=500,pc=0.2,popsize=100,n_select=80,v_cap=80,opt_type=1)
        

输出结果如下：