【PDPTW】python调用guribo求解PDPTW问题（Li & Lim‘s benchmark）之二

panbaoran913

已于 2023-01-02 13:06:08 修改

阅读量1.1k

点赞数 3

分类专栏：车辆路径规划问题文章标签： python guribo PDPTW

于 2023-01-01 15:48:43 首次发布

本文链接：https://blog.csdn.net/panbaoran913/article/details/128512553

版权

车辆路径规划问题专栏收录该内容

21 篇文章 15 订阅

订阅专栏

原文连接：知乎《使用Python调用Gurobi求解PDPTW问题（Li & Lim’s benchmark）》

分析文章：《[PDPTW]解读Python调用Gurobi求解PDPTW问题（Li & Lim‘s benchmark)之一》

修改的代码在本篇。在这里使编码更加简洁。

全部代码在Gitee :https://gitee.com/pan-baoran/pdptw_gurobi_lc101

文章目录

修改

以及修改公示约束（8）与代码不符合的问题。
约束(8)
添加depot的时间窗口约束,添加节点的时间窗口约束.这里的V包含depot和节点.depot的id为0.
$\forall i \in V, \forall k \in K$ ,
$e_i \leq B_i^q \leq l_i$

    # 约束（8）depot和customer的时间窗约束
    for k in Pro.Vehicles.keys():
        # 两个depot的时间窗约束
        model.addConstr(b[0, k] >= Pro.TimeWindow[0][0])
        model.addConstr(b[0, k] <= Pro.TimeWindow[0][1])
        model.addConstr(b[node_num - 1, k] >= Pro.TimeWindow[node_num - 1][0])
        model.addConstr(b[node_num - 1, k] <= Pro.TimeWindow[node_num - 1][1])
        model.addConstr(b[1, k] >= Pro.TimeWindow[1][0])
        #节点的时间窗口约束
        for i in range(node_num):
            model.addConstr(b[i,k] >= Pro.TimeWindow[i][0])
            model.addConstr(b[i,k] <= Pro.TimeWindow[i][1])

当depot不分为depot_start和depot_end的时候，修改为：

    for k in Pro.Vehicles.keys():
        for i in range(node_num):
            model.addConstr(b[i,k] >= Pro.TimeWindow[i][0])
            model.addConstr(b[i,k] <= Pro.TimeWindow[i][1])

现在代码可以分为两篇：utlis.py和test.py

utlis.py

class Pro():用于存储本问题的相关信息
def read_txt_data(DataPath,Pro)：读取车辆信息，并将结果放入Pro
def calculate_euclid_distance(x1, y1, x2, y2):定义了两个节点的距离公示
def Construct_DisTime_Matrix(Pro)：计算任意两个节点的距离矩阵（与车辆无关）和时间矩阵（含车辆）
def build_model(Pro)：建立PDPTW的问题模型


from gurobipy import *
import numpy as np
import pandas as pd
import math
import random
import time
class Pro():
    def __init__(self):
        self.Vehicles={} #用于存放车辆，适用于单depot的情况
        #self.Depots={}#用于存放depot,使用于多depot的情况。{d1:{车辆集1}，d2:{车辆集2}，...}
        self.Locations={}
        self.Demand={}
        self.TimeWindow={}
        self.ServiceTime={}
        self.Request={}
        self.node_num=None
        self.veh_num=None
        self.pair_num=None #成对约束的数量
        self.EarliestTime=None##depot的服务的最早时间
        self.LatestTime=None #depot的服务的最晚时间
        self.DistanceMatrix={}#用于存放节点(包括depot和需求节点)彼此之间的距离
        self.LongestDistance=None #最大节点间的距离
        self.TimeMatrix={}
def read_txt_data(DataPath,Pro):
    # 读取车辆信息u
    VehiclesInfo = pd.read_table(DataPath, nrows=1, names=['K', 'C', 'S'])

    for i in range(int(VehiclesInfo.iloc[0, 0])):
        Pro.Vehicles[i] = [VehiclesInfo.iloc[0, 1], VehiclesInfo.iloc[0, 2]]  # 键i=车辆的序号，值为[车辆容量、速度]

    Pro.veh_num=len(Pro.Vehicles)
    # 读取Depot和任务信息
    ColumnNames = ['TaskNo', 'X', 'Y', 'Demand', 'ET', 'LT', 'ST', 'PI', 'DI']
    # [任务标号,x坐标,y坐标,需求任务]
    TaskData = pd.read_table(DataPath, skiprows=[0], names=ColumnNames, index_col='TaskNo')

    # 提取Depot和取送货点（Customer）的位置坐标 Locations
    Pro.node_num = TaskData.shape[0]
    for i in range(Pro.node_num):
        if i not in Pro.Locations.keys():
            Pro.Locations[i] = [TaskData.iloc[i, 0], TaskData.iloc[i, 1]]  # 键为depot或客户编号，值为相应的坐标（x，y）

    # 提取Depot和取送货点的Demand
    for i in range(Pro.node_num):
        if i not in Pro.Demand.keys():
            Pro.Demand[i] = TaskData.iloc[i, 2]

    # 提取Depot和取送货点的最早和最晚取送货时间及时间窗

    Pro.EarliestTime = TaskData.sort_values(by='ET').iloc[0, 3]
    Pro.LatestTime = TaskData.sort_values(by='LT', ascending=False).iloc[0, 4]
    for i in range(Pro.node_num):
        if i not in Pro.TimeWindow.keys():
            Pro.TimeWindow[i] = [TaskData.iloc[i, 3], TaskData.iloc[i, 4]]
    # 提取Depot和取送货点的服务时间ServiceTime
    for i in range(Pro.node_num):
        if i not in Pro.ServiceTime.keys():
            Pro.ServiceTime[i] = TaskData.iloc[i, 5]

    # 提取运输Request
    # 对于取货任务，PICKUP索引为0，而同一行的DELIVERY索引给出相应送货任务的索引
    count = 0  # 记录运输需求的数量
    for i in range(1, Pro.node_num - 1):
        if TaskData.iloc[i, 6] == 0:
            Pro.Request[count] = [i, TaskData.iloc[i, 7]]  # 将取送货点组合在一起，键为count，值为[取货点，送货点]
            count += 1
    Pro.pair_num=count


def calculate_euclid_distance(x1, y1, x2, y2):
    """  计算Euclid距离.

    具体描述：计算两点(x1, y1), (x2, y2)之间的Euclid距离.
    """
    EuclidDistance = math.sqrt((x2-x1)**2 + (y2-y1)**2)
    return EuclidDistance

def Construct_DisTime_Matrix(Pro):
    #获得两个节点之间距离和时间距离
    Locations=Pro.Locations
    for i in Locations.keys():
        for j in Locations.keys():
            if i != j:
                Pro.DistanceMatrix[i, j] = calculate_euclid_distance(Locations[i][0], Locations[i][1],
                                                                    Locations[j][0], Locations[j][1])
    key = max(Pro.DistanceMatrix, key=Pro.DistanceMatrix.get)
    Pro.LongestDistance = Pro.DistanceMatrix[key]
    for k in range(Pro.veh_num):
        for i in Locations.keys():
            for j in Locations.keys():
                if i != j:
                    Pro.TimeMatrix[i, j, k] = Pro.DistanceMatrix[i, j] / Pro.Vehicles[k][1]



def build_model(Pro):
    # 建立pdptw问题的模型.
    # 创建模型
    model = Model()  # gurobipy中的Model()
    # 变量下标存放列表
    node_num=Pro.node_num
    ijk=[(i,j,k) for i in range(node_num) for j in range(node_num) if i!=j for k in range(Pro.veh_num)]
    ik=[(i,k) for i in range(node_num) for k in range(Pro.veh_num)]
    # 添加决策变量x_ijk，意为车辆k经过弧（i，j）
    # 添加状态变量q_ik,意为车辆k离开节点i的存储量
    # 添加状态变量b_ik,意为车辆k到达节点i的时间
    x = model.addVars(ijk, vtype=GRB.BINARY, name='x')  # 变量x_{ijk},GRB.BINARY表示{0,1}变量
    q = model.addVars(ik, vtype=GRB.INTEGER, name='q')  # 变量q_{ik},GRB.INTEGER为整数变量
    b = model.addVars(ik, vtype=GRB.CONTINUOUS, name='b')  # 变量b_{ik},GRB.CONTINUOUS为连续变量
    model.update()#将变量跟新加入，不添加该行的话，查看model的变量为空
    ## ================== VRPTW===================

    # 约束（1）every vertex has to be served exactly once，除了depot外，每个节点只被服务一次
    for i in range(1, node_num - 1):
        model.addConstr(x.sum(i, '*', '*') == 1,name="eq1-i%s"%i)  # 星号*的用法

    # 约束(2-3) guarantee that every vehicle starts at the depot and returns to the depot
    # at the end of its route. 每辆车必须从depot出发，最后回到depot
    # 约束(4) Flow conservation 流平衡约束
    for k in Pro.Vehicles.keys():
        model.addConstr(x.sum(0, '*', k) == 1,name="eq2-k%s"%k)  # 车辆k从depot出发
        model.addConstr(x.sum('*', node_num - 1, k) == 1,name="eq3-k%s"%k)  # 车辆k回到depot
        for i in range(1, node_num - 1):  # 车辆k在P和D的流平衡约束
            model.addConstr(x.sum(i, '*', k) == x.sum('*', i, k),name="eq4-k%s-i%s"%(k,i))

    # 约束 (5) Time variables are used to eliminate subtours,用时间变量来消除子回路约束
    # 需要用大M法来线性化该约束，M定义为 2*(LatestTime+LongestDistance)
    vars=model.getVars()
    for i,j,k in ijk:
        if i !=j:
            model.addConstr(b[j,k]+2*(1-x[i,j,k])*(Pro.LatestTime+Pro.LongestDistance) >=
            b[i,k]+Pro.ServiceTime[i]+Pro.TimeMatrix[i,j,k],name="eq5-k%s-i%s-j%s"%(k,i,j))

    # 约束(6-7) guarantee that a vehicle’s capacity is not exceeded throughout its tour，载货量平衡与车辆载量约束
    # 约束(6) 载货量平衡约束，需要用大M法来线性化该约束，M定义为 100*车辆最大载量

    for i,j,k in ijk:
        if i!=j:
            model.addConstr(q[j,k]+(1-x[i,j,k])*(100*Pro.Vehicles[k][0])>=q[i,k]+Pro.Demand[j],name="eq6-k%s-i%s-j%s"%(k,i,j))

    # 约束（7）车辆载量约束
    for i,k in ik:
        model.addConstr(q[i,k] >=0, name="eq7.1-k%s-i%s"%(k,i))
        model.addConstr(q[i,k] >=Pro.Demand[i],name="eq7.2-k%s-i%s"%(k,i))
        model.addConstr(q[i,k] <= Pro.Vehicles[k][0],name="eq7.3-k%s-i%s"%(k,i))
        model.addConstr(q[i,k] <=Pro.Vehicles[k][0]+Pro.Demand[i],name="eq7.4-k%s-i%s"%(k,i))

    # 约束（8）depot和customer的时间窗约束
    for k in Pro.Vehicles.keys():
        # 两个depot的时间窗约束
        model.addConstr(b[0, k] >= Pro.TimeWindow[0][0])
        model.addConstr(b[0, k] <= Pro.TimeWindow[0][1])
        model.addConstr(b[node_num - 1, k] >= Pro.TimeWindow[node_num - 1][0])
        model.addConstr(b[node_num - 1, k] <= Pro.TimeWindow[node_num - 1][1])
        #节点的时间窗口约束
        for i in range(node_num):
            model.addConstr(b[i,k] >= Pro.TimeWindow[i][0])
            model.addConstr(b[i,k] <= Pro.TimeWindow[i][1])
    # 旧代码的写法
    # for k in Pro.Vehicles.keys():
    #     # 两个depot的时间窗约束
    #
    #     model.addConstr(b[0, k] >= Pro.TimeWindow[0][0])
    #     model.addConstr(b[0, k] <= Pro.TimeWindow[0][1])
    #     model.addConstr(b[node_num - 1, k] >= Pro.TimeWindow[node_num - 1][0])
    #     model.addConstr(b[node_num - 1, k] <= Pro.TimeWindow[node_num - 1][1])
    #     Request=Pro.Request
    #     for r in range(Pro.pair_num):
    #         # 运输请求i两个节点的左时间窗
    #         model.addConstr(b[Request[r][0], k] >= Pro.TimeWindow[Request[r][0]][0])
    #         model.addConstr(b[Request[r][1], k] >= Pro.TimeWindow[Request[r][1]][0])
    #         # 运输请求i两个节点的右时间窗
    #         model.addConstr(b[Request[r][0], k] <= Pro.TimeWindow[Request[r][0]][1])
    #         model.addConstr(b[Request[r][1], k] <= Pro.TimeWindow[Request[r][1]][1])

    ## ================== PDPTW===================

    # 约束（9）both origin and destination of a request must be served by the same vehicle
    # 保证取货后要有对应的送货，取货和送货由同一辆车完成
    for p in range(Pro.pair_num):
        for k in Pro.Vehicles.keys():
            model.addConstr(x.sum(Pro.Request[p][0], '*', k) == x.sum('*', Pro.Request[p][1], k),name="eq9-p%s"%p)

    # 约束 (10) delivery can only occur after pickup,先取后送货约束
    for p in range(Pro.pair_num):
        for k in Pro.Vehicles.keys():
            model.addConstr(b[Pro.Request[p][0], k] <= b[Pro.Request[p][1], k],name="eq10-p%s-k%s"%(p,k))

    ## 设置目标函数
    # 目标函数（3）:最小化车辆总的行驶距离
    DistanceCost = 0
    for k in Pro.Vehicles.keys():
        for i in range(node_num):
            for j in range(node_num):
                if i != j:
                    DistanceCost += (Pro.DistanceMatrix[i, j] * x[i, j, k])
    model.setObjective(DistanceCost, GRB.MINIMIZE)
    model.update()

    model.__data = x, b, q

    return model

test.py

from gurobipy import *
import time
from utlis import Pro, read_txt_data, Construct_DisTime_Matrix, build_model

Pro=Pro()
start = time.time()
# DataPath="lc101.txt"
DataPath="smallcase.txt"
read_txt_data(DataPath,Pro)
Construct_DisTime_Matrix(Pro)

model=build_model(Pro)
model.params.TimeLimit = 1000    # 设定模型求解时间
model.write('PDPTW_%s.lp' % DataPath)
model.write('PDPTW_%s.mps' % DataPath)
model.setParam(GRB.Param.LogFile, 'PDPTW_%s.log' % DataPath)
model.optimize()
x,b,q = model.__data

end = time.time()
print('程序运行时间：')
print(end-start)

程序运行时间：
0.13767194747924805
查看结果

Optimize a model with 2588 rows, 1100 columns and 9060 nonzeros
Variable types: 100 continuous, 1000 integer (900 binary)
Coefficient statistics:
Matrix range [1e+00, 1e+04]
Objective range [1e+01, 4e+01]
Bounds range [1e+00, 1e+00]
RHS range [1e+00, 1e+04]

Presolve removed 740 rows and 50 columns
Presolve time: 0.03s
Presolved: 1848 rows, 1050 columns, 11628 nonzeros
Variable types: 100 continuous, 950 integer (850 binary)

Root relaxation: objective 8.000000e+01, 60 iterations, 0.00 seconds (0.00 work units)

Cutting planes:
Learned: 1
Gomory: 1
Implied bound: 39
MIR: 13
StrongCG: 3
RLT: 5
Relax-and-lift: 3

Explored 68 nodes (2136 simplex iterations) in 0.60 seconds (0.36 work units)
Thread count was 8 (of 8 available processors)

输出结果

# 输出方式2
for k in Pro.Vehicles.keys():
    for i in range(Pro.node_num):
        for j in range(Pro.node_num):
            if i != j:
                if (x[i,j,k].X- 0.5 > 0):
                    print('tour for vehicle %s:' % k)
                    print('%s-->%s' % (i,j))

out:
tour for vehicle 0:
0–>9
tour for vehicle 1:
0–>9
tour for vehicle 2:
0–>9
tour for vehicle 3:
0–>9
tour for vehicle 4:
0–>9
tour for vehicle 5:
0–>9
tour for vehicle 6:
0–>9
tour for vehicle 7:
0–>9
tour for vehicle 8:
0–>9
tour for vehicle 9:
0–>1
tour for vehicle 9:
1–>2
tour for vehicle 9:
2–>3
tour for vehicle 9:
3–>4
tour for vehicle 9:
4–>5
tour for vehicle 9:
5–>6
tour for vehicle 9:
6–>7
tour for vehicle 9:
7–>8
tour for vehicle 9:
8–>9

运行DataPath=“lc101.txt”

运行该程序，会出现模型不可解的现象。运行下面的代码。检查哪个约束条件出了问题。

if model.status == GRB.Status.INFEASIBLE:
        print('Optimization was stopped with status %d' % model.status)
        # do IIS, find infeasible constraints
        model.computeIIS()
        for c in model.getConstrs():
            if c.IISConstr:
                print('%s' % c.constrName)

out:

eq1-i97
eq3-k9
eq3-k13
eq9-p50-k9
eq9-p50-k13

检查：
公式1
$\forall i \in P\cup D=\{1,2,\cdots,n,n+1,\cdots,n+n\}$
$\sum_{k\in K}\sum_{j:(i,j)\in A}x^k_{ij}=1 \tag{1}$
对于i=97的时候。
在这里插入图片描述
发现，与“smallcase.txt”最后位置少一行depot_end的一行数据.所以，导致在运行过程中，删除了97的对应节</f点。
公式3
$\forall k \in K$
$\sum_{j:(0,j)\in A}x_{0,j}^k=1$
$\sum_{i:(i,n+\tilde{n}+1)\in A} x^k_{i, n+\tilde{n}+1}=1$

for k in Vehicles.keys():
    model.addConstr(x.sum(0,'*',k) == 1)       # 车辆k从depot出发
    model.addConstr(x.sum('*',nrows-1,k) == 1)       # 车辆k回到depot

表示：对于车辆k来说，仅有一次从depot出发，并仅有一次返回depot.这里的nrows-1值得使106行。而不是depot行。

约束(9)
弧(i,j)必须由同一辆车经过(先取货再送货).假设这些弧集为特定的集合 $\tilde{A}$
$\forall k \in K$ , $\forall (i,j)\in \tilde{A}$
$\sum_{j:(i,j)\in \tilde{A}}x_{ij}^k=\sum_{i:(i,j)\in \tilde{A}}x_{ij}^k$

    for r in range(len(Request)):
        for k in range(len(Vehicles)):
            model.addConstr(x.sum(Request[r][0],'*',k)==x.sum('*',Request[r][1],k))

Request = {…, 48: [92, 93], 49: [96, 94], 50: [97, 106], 51: [98, 95], 52: [100, 99]}
检查Request中的key=50的值为：50 : [97, 106]

经查验，就是因为两个文件中depot的数量不同引起的。在”smallcase.txt"中将depot划分为depot_start和depot_end.而在“lc101.txt”中仅仅有一个depot.

要检查的地方有

数据读取的方式

#def read_txt_data(DataPath,Pro)
    count = 0  # 记录运输需求的数量
    for i in range(1, Pro.node_num ):
        if TaskData.iloc[i, 6] == 0:
            Pro.Request[count] = [i, TaskData.iloc[i, 7]]  # 将取送货点组合在一起，键为count，值为[取货点，送货点]
            count += 1
    Pro.pair_num=count

公式3：关于depot的索引修改

# def build_model(Pro):
## 约束3中，返回depot的数字id为0，而不死node_num-1
    for k in Pro.Vehicles.keys():
        model.addConstr(x.sum(0, '*', k) == 1,name="eq2-k%s"%k)  # 车辆k从depot出发
        # model.addConstr(x.sum('*', node_num - 1, k) == 1,name="eq3-k%s"%k)
        model.addConstr(x.sum('*', 0, k) == 1,name="eq3-k%s"%k)# 车辆k回到depot ###公式3
        for i in range(1, node_num - 1):  # 车辆k在P和D的流平衡约束
            model.addConstr(x.sum(i, '*', k) == x.sum('*', i, k),name="eq4-k%s-i%s"%(k,i))