采用Python+Gurobi求解带有集群回程需求的VRPCB问题

1. 模型

1.1 VRPB问题介绍

带有回程需求的VRP问题（VRP with Backhauls,VRPB)最初由 Deif, . Bodin, LD. 于1984年出。相关学者将其进一步细分为四类子问题：
第一类VRPCB：

• 第一类VRPB(VRPCB)：
  – 客户只能是linehaul客户和backhaul客户中的一类
  – 车辆在访问backhaul客户集群前，必须先向linehaul客户集群交付货物
  – 客户只能被访问一次

• 第二类VRPB（VRPBM)
  – 不考虑客户集群限制，车辆可以交叉服务linehaul客户和backhaul客户
  – 客户只能被访问一次

• 第三类VRPB（VRPPD）
  – 客户可以同时是linehaul客户和backhaul客户
  – 客户可以2次被访问，即：先访问几个客户交付货物，以便部分清空车辆装载。然后，再返回访问该客户提货

• 第四类VRPB（VRPPD)
  – 客户可以同时是linehaul客户和backhaul客户
  – 客户只能被访问一次

1.2 数学模型

这里暂时只研究VRPCB问题，借用 Toth, P., Vigo, D. 学者于1997给出的整数线性规模。

1.2.1 模型参数

1.2.2 数学模型

1.2.3 模型分解

上述模型在复现时有一定困难，尤其是 约束（6）和约束（7）中的 $$\mathcal{S},\mathcal{F}$$
这里根据 帖子（其实有一点小bug）提供的分解模型进行复现，模型如下：

2. 数据结构

原网址暂时想不起来了，回头再补。这里直接贴上代码读取的文件结构，各列依次是：节点id，节点横坐标，节点纵坐标，linehaul客户需求和backhaul客户需求

3. Gurobi源码

import copy
import csv
import math
import xlsxwriter
import matplotlib.pyplot as plt
from gurobipy import GRB,Model,quicksum

# 读取文件
def read_data(filename):
    Q = {}
    L = []
    B = []
    V = []
    XY = {}
    depot = None
    # 读取网络节点、需求
    with open(filename, 'r') as f:
        node_reader = csv.DictReader(f)
        for row in node_reader:
            if float(row['linehaul']) == 0 and  float(row['backhaul']) == 0:
                depot = row['id']
                Q[depot] = 0
            elif float(row['linehaul']) > 0 and  float(row['backhaul']) == 0:
                L.append(row['id'])
                Q[row['id']] = float(row['linehaul'])
            elif float(row['linehaul']) == 0 and  float(row['backhaul']) > 0:
                B.append(row['id'])
                Q[row['id']] = float(row['backhaul'])
            V.append(row['id'])
            XY[row['id']] = (float(row['x_coord']), float(row['y_coord']))
    # 计算网络弧
    L0 = L + [depot]
    B0 = B + [depot]
    AL = [ (i, j) for i in L0 for j in L if i != j ]
    AB = [ (i, j) for i in B for j in B0 if i != j]
    AC = [ (i, j) for i in L for j in B0 if i != j]
    # 计算网络弧距离
    Cost = {}
    for i in V:
        x1, y1 = XY[i][0], XY[i][0]
        for j in V:
            x2, y2 = XY[j][0], XY[j][0]
            Cost[i,j] = math.sqrt( (x1-x2)**2 + (y1-y2)**2 )
    return depot,L,L0,B,B0,AL,AB,AC,Q,Cost,XY
# 提取结果，形成车辆路径
def extract_routes(depot,L,B,B0,X,Y,Z):
    L = copy.deepcopy(L)
    B = copy.deepcopy(B)
    B0 = copy.deepcopy(B0)
    route_list = []
    V = []
    while len(L):
        # 提取 派送阶段路径
        route = [depot]
        cur_node = depot
        for j in L:
            if X[depot, j].x > 0:
                cur_node = j
                route.append(j)
                L.remove(j)
                break
        stop = True
        while len(L) > 0:
            for j in L:
                if X[cur_node, j].x > 0:
                    cur_node = j
                    route.append(j)
                    if j != depot:
                        L.remove(j)
                    stop = False
                    break
            if stop:
                break
            else:
                stop = True
        # 提取 取货阶段路径
        for j in B:
            if Z[cur_node,j].x > 0:
                cur_node = j
                route.append(j)
                B0.remove(j)
                break
        if cur_node in B:
            while cur_node != depot:
                for j in B0:
                    if Y[cur_node, j].x > 0:
                        cur_node = j
                        route.append(j)
                        if j != depot:
                            B0.remove(j)
                        break
        else:
            route.append(depot)
        route_list.append(route)
        V.extend(route[1:-1])
    print(len(V))
    return route_list
# 绘制车辆路径
def draw_routes(route_list,XY,L,B):
    for route in route_list:
        path_x = []
        path_y = []
        for n in route:
            path_x.append(XY[n][0])
            path_y.append(XY[n][1])
        plt.plot(path_x, path_y, linewidth=0.5, ms=5,color='black')
        linehual_point_x = [XY[n][0] for n in L]
        linehual_point_y = [XY[n][1] for n in L]
        backhual_point_x = [XY[n][0] for n in B]
        backhual_point_y = [XY[n][1] for n in B]
        plt.scatter(linehual_point_x, linehual_point_y, marker='s', c='b', s=5)
        plt.scatter(backhual_point_x, backhual_point_y, marker='o', c='r', s=5)
        plt.show()
# 保存结果
def save_file(route_list,total_cost,Cost):
    wb = xlsxwriter.Workbook('路径方案.xlsx')
    ws = wb.add_worksheet()
    ws.write(0,0,'总费用')
    ws.write(0,1,total_cost)
    ws.write(1,0,'车辆')
    ws.write(1,1,'路径')
    ws.write(1,2,'距离')
    for row,route in enumerate(route_list):
        route_str = [str(i) for i in route]
        dist = sum(Cost[route[i], route[i + 1]] for i in range(len(route) - 1))
        ws.write(row + 2, 0, f'{row + 1}')
        ws.write(row+2,1,'-'.join(route_str))
        ws.write(row + 2, 2, dist)
        row += 1
    wb.close()
# 建模和求解
def solve_model(depot,L,L0,B,B0,AL,AB,AC,Q,Cost,K,CAP,XY):
    """

    :param depot:车场id
    :param L:linehaul节点集合
    :param B:backhaul节点集合
    :param AL:linehaul节点衔接弧集合
    :param AB:backhaul节点衔接弧集合
    :param AC:linehaul节点和backhaul节点衔接弧集合
    :param Q:节点需求集合
    :param Cost:网络弧费用
    :return:
    """
    model = Model()
    # 添加变量
    X = model.addVars(AL,vtype=GRB.BINARY,name='X[i,j]')
    Y = model.addVars(AB,vtype=GRB.BINARY,name='Y[i,j]')
    Z = model.addVars(AC,vtype=GRB.BINARY,name='Z[i,j]')
    U1 = model.addVars(L0,vtype=GRB.CONTINUOUS,name='U[i]')
    U2 = model.addVars(B0, vtype=GRB.CONTINUOUS, name='U[i]')
    # 目标函数
    obj = (quicksum(X[i,j]*Cost[i,j] for i,j in AL) + quicksum(Y[i,j]*Cost[i,j] for i,j in AB) +
           quicksum(Z[i,j]*Cost[i,j] for i,j in AC))
    model.setObjective(obj,GRB.MINIMIZE)
    # linebaul 相关约束
    model.addConstr( quicksum(X[depot,j] for j in L) == K ) # 车辆数约束
    model.addConstrs( (quicksum(X[i,j] for i in L0 if i != j) == 1 for j in L) ) # 派送需求约束
    model.addConstrs( (U1[i] - U1[j] + CAP*X[i,j] <= CAP - Q[j] for i,j in AL) ) # 破圈约束
    # backbaul 相关约束
    model.addConstr( quicksum([Y[j,depot] for j in B]) + quicksum(Z[i,depot] for i in L) == K ) # 车辆数约束
    model.addConstrs( (quicksum(Y[i,j] for j in B0 if i != j) == 1 for i in B) ) # 取货需求约束
    model.addConstrs( (U2[i] - U2[j] + CAP*Y[i,j] <= CAP - Q[j] for i,j in AB) ) # 破圈约束
    # connection 相关约束
    model.addConstr( quicksum(Z[i,j] for i,j in AC) == K ) # 车辆数约束
    model.addConstrs( quicksum(X[i,j] for j in L if i != j ) + quicksum(Z[i,j] for j in B0) <= 1 for i in L ) # 接续约束
    model.addConstrs( quicksum(Y[i,j] for i in B if i != j) + quicksum(Z[i,j] for i in L) <= 1 for j in B) # 接续约束
    # 模型求解
    model.Params.TimeLimit = 300  # 规模较大时可设置求解时间限制
    model.optimize()
    # 判断求解状态
    if model.status == GRB.Status.OPTIMAL or model.status == GRB.Status.TIME_LIMIT:
        route_list = extract_routes(depot,L,B,B0,X,Y,Z)
        draw_routes(route_list, XY, L, B)
        save_file(route_list, model.objVal,Cost)
if __name__=='__main__':
    filename=r'demand-X-n120-50-k3.csv'
    depot,L,L0,B,B0,AL,AB,AC,Q,Cost, XY = read_data(filename)
    solve_model(depot=depot,L=L,L0=L0,B=B,B0=B0,AL=AL,AB=AB,AC=AC,Q=Q,Cost=Cost,K=3,CAP=21,XY=XY)

4. 求解结果

参考

1. https://blog.csdn.net/qq_44149045/article/details/128941919
2. Koc,Cagri,Laporte,et al.Vehicle routing with backhauls: Review and research perspectives[J].Computers, 2018.
3. Parragh, S.N., Doerner, K.F. & Hartl, R.F. A survey on pickup and delivery problems . Journal für Betriebswirtschaft 58, 21–51 (2008). https://doi.org/10.1007/s11301-008-0033-7