python调用SCIP求解1-PDVRP （MTZ约束消除子环路）

1. 1-PDVRP 问题定义

One-Commodity M-M Pickup and Delivery Vehicle Routing Problem (1-PDVRP) ，a new variant of classical VRP. It differs from the general VRP with pickup and delivery in that the commodity provided by the pickup customers is the same as that needed by delivery customers. This is an important characteristic which distinguishes the 1-PDVRP from general vehicle routing problem with pickup and delivery (VRPPD)

2. 1-PDVRP

2.1 数学模型

数学模型见参考文献[1]

2.1 消除子环路（Miller-Tucker-Zemlin剪枝法）

一种经典方法是使用Miller-Tucker-Zemlin (MTZ) 列举变量技巧，在线性规划模型中加入额外的变量和约束来防止子环路的出现。具体做法是为每辆车引入一串顺序变量$ (u_i)（其中 (i) 代表节点）$，代表车辆访问节点的顺序。然后添加以下约束：
$$u_j-u_i+|N|\cdot x_{ij}\le |N|-1\text{对于所有},(i,j)\in A,i\ne j$$

这里，$(N)$ 是节点集，$(A)$ 是边集，$(x_{ij})$ 是决策变量，表示车辆从节点 $(i)$ 到节点 $(j)$ 是否行驶，$(|N|)$ 是节点总数。这个约束确保如果从$(i)$ 到 $(j)$ 有行驶$(x_{ij}=1)$，那么 $(j)$ 节点被访问的顺序必须在 $(i)$ 之后，除非它们属于同一个子环路，这通过限制 $(u_j-u_i)$ 的值来避免。

3. python调用SCIP求解1-PDVRP （MTZ约束消除子环路）

3.1 代码

import os
import re
import math
import pandas as pd
import numpy as np
import itertools
from collections import defaultdict
import pyscipopt as opt
import matplotlib.pyplot as plt
import networkx as nx


from pyscipopt import Model, Conshdlr, quicksum, SCIP_RESULT


class One_M_M_PDVRP(object):
    def __init__(self):
        self.vehicleNum  = 2         # 车辆数
        self.nodeNum     = None      # 顶点数量
        self.capacity    = 40        # 车容量
        self.cor_X       = []        # x轴坐标
        self.cor_Y       = []        # y轴坐标
        self.disMatrix   = [[]]      # 距离成本
        self.demand      = []        # 需求

    def prepare_data(self):
        data = pd.DataFrame([[  0, 250, 250, 0],
                             [  1, 387, 297, -10],
                             [  2,   5, 297,  10],
                             [  3, 355, 177,  20],
                             [  4,  78, 346, -30],
                             [  5, 286, 159,  20],
                             [  6, 322, 465, -10],
                             [  7, 393, 408,  30],
                             [  8,  89, 216, -10],
                             [  9,  76, 345,  30],
                             [ 10, 410, 285,  20]], columns=['NO', 'XCOORD', 'YCOORD', 'DEMAND'])

        self.nodeNum = len(data)
        self.cor_X = data['XCOORD'].values
        self.cor_Y = data['YCOORD'].values
        # 需求q > 0 代表 a pickup customer (顾客点提货), 需求q < 0 代表 a delivery customer(顾客点交货)
        self.demand = data['DEMAND'].values
        
        # 计算欧式距离矩阵
        self.disMatrix = np.full((self.nodeNum, self.nodeNum), 0)
        for i in range(self.nodeNum):
            for j in range(self.nodeNum):
                self.disMatrix[i][j] = math.sqrt((self.cor_X[i] - self.cor_X[j]) ** 2 + (self.cor_Y[i] - self.cor_Y[j]) ** 2)
    
    def print_vehicle_route(self, pdp_route):
        # 使用defaultdict来收集属于同一车辆的所有边
        routes_by_vehicle = defaultdict(list)
        for start, end, vehicle in pdp_route:
            routes_by_vehicle[vehicle].append((start, end))
    
        # 构建车辆及其路径的表示
        vehicles_and_paths = {}

        for vehicle_id, edges in routes_by_vehicle.items():
            routes = []
            for i in range(len(edges)):
                if i == 0:
                    routes.append(edges[0][:2])
                    edges.pop(0)
                else:
                    pre = routes[-1]
                    connected_node = pre[1] 
                    for j in range(len(edges)):
                        aft = edges[j]
                        if connected_node in aft:
                            if aft[0] == connected_node:
                                chosen_node = (aft[0], aft[1])
                            elif aft[1] == connected_node:
                                chosen_node = (aft[1], aft[0])
                            routes.append(chosen_node)
                            edges.pop(j)
                            break
                i = i+1
            
            # 解析显示: 从第一个为需求>0的点开始。点不一定从仓库点出发
            path_route = [edge[0] for edge in routes]
            route_demand = [self.demand[i] for i in path_route]
            start_idx = next((index for index, value in enumerate(route_demand) if value > 0), None)
            adj_routes = path_route[start_idx:] + path_route[:start_idx]
            
            path = " -> ".join(str(edge) for edge in adj_routes)
            vehicles_and_paths[vehicle_id] = path
            
            print(f'车辆 {vehicle_id} 路径 {path}')

    def plt_route_pic(self, vrp_route):
        plt.figure(figsize=(10, 10), dpi=200)
        # 绘制坐标的散点图
        for i in range(0, self.nodeNum):
            if i == 0:
                plt.scatter(self.cor_X[i], self.cor_Y[i], marker='^', c='r')
                plt.annotate('depot', (self.cor_X[i], self.cor_Y[i]))
            else:
                plt.scatter(self.cor_X[i], self.cor_Y[i], marker='o', c='black')
                if self.demand[i] > 0:
                    # action_demand = str((i, '+' + str(self.demand[i])))
                    action_demand = '(' + str(i) + ',+' + str(self.demand[i]) + ')'
                else:
                    action_demand = '(' + str(i) + ',' + str(self.demand[i]) + ')'
                plt.annotate(action_demand, (self.cor_X[i], self.cor_Y[i] + 5) )
        # 绘制箭头
        for idx in vrp_route:
            i,j,k = idx
            plt.arrow(self.cor_X[i], self.cor_Y[i], self.cor_X[j] - self.cor_X[i], self.cor_Y[j] - self.cor_Y[i],
                      length_includes_head=True, head_width=4.5, head_length=5.8, fc='blue', ec='blue', linewidth=1, shape='full')

        plt.xlabel('x')
        plt.ylabel('y')
        plt.title('1-PDVRP')
        plt.savefig('./1-PDVRP.png')

    def solve_pdp(self):
        model = opt.Model()
        # ==========定义变量==========
        x = {}
        f = {}
        # x_i_j_k: 0-1变量，表示车辆k经过弧（i，j）, i != j
        for i in range(self.nodeNum):
            for j in range(self.nodeNum):
                if i != j:
                    f[i, j] = model.addVar(vtype='I', name='f_' + str(i) + '_' + str(j), lb=0)
                    for k in range(self.vehicleNum):
                        x[i, j, k] = model.addVar(vtype='B', name='x_' + str(i) + '_' + str(j) + '_' + str(k))

        # ==========定义约束==========
        # 约束6.2：对于客户点: 保证每个客户点都被服务一次
        for i in range(1, self.nodeNum):
            model.addCons(opt.quicksum(x[i, j, k] for j in range(self.nodeNum) for k in range(self.vehicleNum) if i!=j) == 1,
                          name='customer_visited_' + str(i))

        # 约束6.3: 对于客户点: 进出是同一辆车，流平衡，服务之后需离开         
        for i in range(self.nodeNum):
            for k in range(self.vehicleNum):
                model.addCons(opt.quicksum(x[i, j, k] for j in range(self.nodeNum) if i!=j) ==
                              opt.quicksum(x[j, i, k] for j in range(self.nodeNum) if i!=j),
                              name='flow_' + str(k) + '_' + str(i))
        
        # 约束6.4: 限制了每个弧上的流量最多为Q 
        for i in range(0, self.nodeNum):
            for j in range(0, self.nodeNum):
                if i != j:
                    model.addCons(f[i, j] <= opt.quicksum(x[i, j, k] for k in range(self.vehicleNum)) * self.capacity)
        
        # 约束6.5: 流量守恒约束，保证了客户点的需求。流量守恒不足以确保所有路线都从depot开始和结束
        for i in range(1, self.nodeNum):
            model.addCons(opt.quicksum(f[i, j] for j in range(self.nodeNum) if i!=j) -
                          opt.quicksum(f[j, i] for j in range(self.nodeNum) if i!=j) == self.demand[i])
                    
        # 约束6.6: Miller-Tucker-Zemlin法 消除子环路
        u = {}
        for k in range(self.vehicleNum):
            for i in range(self.nodeNum):
                u[k, i] = model.addVar(vtype='I', name='u_' + str(k) + '_' + str(i), lb=0)
        
        for k in range(self.vehicleNum):
            for i in range(1, self.nodeNum):
                for j in range(1, self.nodeNum):
                    if i != j:
                        model.addCons(u[k, i] - u[k, j] + self.nodeNum * x[i, j, k] <= self.nodeNum - 1)
                    
        # 自定义约束: 每个车辆行驶线路最大距离约束（某些文献中有此约束）
        # for k in range(self.vehicleNum):
        #     model.addCons(opt.quicksum(x[i, j, k] * self.disMatrix[i, j] for i in range(self.nodeNum) for j in range(self.nodeNum) if i !=j) <= 1000)
                    
        # ==========定义目标==========
        model.setObjective(opt.quicksum(x[i, j, k] * self.disMatrix[i, j] for (i, j, k) in x))
        model.setMinimize()
        model.setParam('limits/gap', 0.001)
        model.optimize()

        # ==========输出结果==========
        if model.getStatus() != 'infeasible':
            print('model_status = ', model.getStatus())
            print('model_gap =', model.getGap())
            print('model_obj =', model.getObjVal())
    
            # model.writeProblem('pdp.lp')
            pdp_route = []
            for k in range(self.vehicleNum):
                for i in range(self.nodeNum):
                    for j in range(self.nodeNum):
                        if(i != j and model.getVal(x[i,j,k]) > 0):
                            pdp_route.append((i,j,k))      
            # print('pdp_route: ', pdp_route) 
        
            # 打印车辆路径
            self.print_vehicle_route(pdp_route)
            # 绘制结果路径
            self.plt_route_pic(pdp_route)
        else:
            print('model is infeasible')
        
if __name__ == '__main__':
    PDP = One_M_M_PDVRP()
    # 预处理数据
    PDP.prepare_data()
    # 求解模型并绘制结果
    PDP.solve_pdp()
        

3.2 求日志 & 结果

参考文献

[1] Toth P , Vigo D .Vehicle Routing: Problems, Methods, and Applications, Second Edition[J].Society for Industrial and Applied Mathematics, 2014.DOI:10.1137/1.9781611973594.ch9.
[2] Hipólito Hernández-Pérez,Juan-José Salazar-González.A branch-and-cut algorithm for a traveling salesman problem with pickup and delivery[J].Discrete Applied Mathematics, 2004, 145(1):126-139.DOI:10.1016/j.dam.2003.09.013.