Post

带时间窗的取送货路径规划问题PDPTW

Posted Updated
By ZoeyWang
23 min read
带时间窗的取送货路径规划问题PDPTW

1.问题描述

带时间窗的取送货路径规划问题（pickup and delivery problem with time window，PDPTW）是路径规划问题的一个重要子类，在供应链领域比较常见，例如某物流企业需要指定交易未来某段时间内车辆的运输调度计划，在满足货物运输要求前提下最小化运输成本。

（1）对车辆，已知起始停靠点，结束停靠点，最早开始时间，最晚结束时间，最大行驶距离，速度，最大载量，匹配关系等信息

（2）对订单，已知货物量，提货地点，提货服务时间，提货时间窗，送货地点，送货服务时间，送货时间窗，匹配关系等信息

所有车辆从其对应的起始停靠点，依次取送货，最后回到结束停靠点，假设有V1，V2，V3三辆车，六个订单A，B，C，D，E，F。运输路径用图表示如下：

2.数据分析

数据有3张表：地点信息表，车辆信息表，订单信息表

（1）地点信息表：所有订单可能取送货的地点，包括停靠点ID，纬度和经度

（2）车辆信息表：可用车辆的信息，包括每辆车的起始和结束点，服务时间，时间和距离限制，运输速度，单位成本，匹配关系等

（3）订单信息表：要取送货的订单集合，包括每个订单的取送货时间窗，服务时间，取送货地点，匹配关系等

3.数学模型

对于数学模型，首先要明确想要模型求解出什么信息？有一些是目标的决策变量，有一些是辅助的决策变量。在该问题中，我们想知道每辆车负责运输那些订单，这些订单的运输线路是那样的？那么最直观的建模思路是：（1）用一个0-1变量来表示是否运输该订单（2）用另一个0-1变量来表示在选中的订单中的运输边。但是（2）实际上已经包括（1），因为选中了某条运输边，就代表这个运输边包含的订单也被选中；并且在Gurobi建模中决策变量不能作为索引，即不能通过（1）的信息来约束（2）。

换成另一个建模思路：对每辆车，获取匹配关系对应的订单的取送货候选点集合，该车辆只会经过该候选点集合组成的边。如果候选点集合是真实的地点，那么会出现一个真实地点比如L1会有多个订单都经过，即L1在多辆车的候选点集合中，决策变量只会体现出是否经过L1，但是不清楚经过L1取哪个订单，因此候选点集合不能是真实的地点，而是虚拟的地点，即假设每个订单的取送货点都是独立的，但是计算距离时采用虚拟点对应真实点的经纬度。

对每辆车：虚拟点包括起始虚拟点，结束虚拟点，取送货虚拟点。

4.代码实现

以下代码建模部分用pulp库实现，可支持调用多种求解器（cplex, gurobi等）进行求解，完整代码在github：Routing-Problem/PDPTW

（1）数据预处理System：对数据进行预处理，获取建模所需的数据信息

PDPTW问题系统包括3个基本对象：地点Location，订单Order，车辆Vehicle，依次读取数据表创建这些对象。接着基于对象信息提取所需数据：（1）所有的虚拟点集合，依次用数字编号，每个节点包含元组信息（订单号，取/送货真实点，取/送货量），奇数为取货，货量为正数，偶数为送货，货量为负数；（2）每辆车可以匹配的订单号集合，可经过虚拟点集合

1
2
3
4
5
6
7
8
9
10
11
12
13
14

# V1
['O1', 'O2', 'O3', 'O5', 'O7', 'O9', 'O10']
[1, 2, 3, 4, 5, 6, 9, 10, 11, 12, 13, 14, 15, 16]

# V2
['O1', 'O2', 'O3', 'O5', 'O7', 'O9', 'O10']
[1, 2, 3, 4, 5, 6, 9, 10, 11, 12, 13, 14, 15, 16]

# V3
['O1', 'O2', 'O3', 'O5', 'O7', 'O9', 'O10', 'O4']
[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16]

# 虚拟点集合，共8个订单，因此16个虚拟点
{1: ('O1', 'L1', 20), 2: ('O1', 'L2', -20), 3: ('O2', 'L1', 4), 4: ('O2', 'L3', -4), 5: ('O3', 'L1', 10), 6: ('O3', 'L4', -10), 7: ('O4', 'L1', 10), 8: ('O4', 'L5', -10), 9: ('O5', 'L1', 10), 10: ('O5', 'L6', -10), 11: ('O7', 'L1', 10), 12: ('O7', 'L8', -10), 13: ('O9', 'L1', 10), 14: ('O9', 'L10', -10), 15: ('O10', 'L2', 2), 16: ('O10', 'L1', -2)}

（2）求解模型Solver：基于System类的数据信息，开始构建模型并求解

  • 添加边变量：对每辆车k，添加0-1变量，表示是否经过边(i,j)，边(i,j)来自车辆k可经过虚拟点集合，除了虚拟点，车k还要经过起始点和结束点
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

def add_edge_var(self):
    # 变量x_kij: 车辆k是否经过边ij，ij是带order维度的虚拟扩展点，以及车辆的起始点
    for veh_k, veh_obj in self.system.vehicle_obj_dict.items():
        for node_i in veh_obj.alter_node_list:
            for node_j in veh_obj.alter_node_list:
                if node_i != node_j:
                    self.x_vars[(veh_k, node_i, node_j)] = LpVariable(name=f'x_{veh_k}_{node_i}_{node_j}',
                                                                      lowBound=0, upBound=1, cat=LpBinary)
    for veh_k, veh_obj in self.system.vehicle_obj_dict.items():
        node_origin, node_des = 0, self.system.dummy_node_count + 1
        for node in veh_obj.alter_node_list:
            self.x_vars[(veh_k, node_origin, node)] = LpVariable(name=f'{veh_k}_{node_origin}_{node}',
                                                                 lowBound=0, upBound=1, cat=LpBinary)
            self.x_vars[(veh_k, node, node_des)] = LpVariable(name=f'{veh_k}_{node}_{node_des}',
                                                              lowBound=0, upBound=1, cat=LpBinary)
  • 添加基本约束：在添加完边变量后，对其有如下约束：（1）保证每个货物都被配送，即对每个取货虚拟点，一定存在某辆车k从该点出发（2）保证取货后要有对应的送货，即车辆k从某个订单的取货虚拟点出发，也要从该订单的送货虚拟点出发（3）路径平衡约束，即每辆车k一定要从起始虚拟点出发，回到结束虚拟点，取送虚拟点的进出边相等。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43

def add_basic_cons(self):
    # (1) 保证每个货物都被配送：对所有虚拟点中的每个pick点，要求必须存在(pick,j)
    for pick_node in self.system.dummy_node_dict.keys():
        if pick_node % 2 != 0:
            lhs_origin = 0
            for key in self.x_vars.keys():
                k, i, j = key[0], key[1], key[2]
                if i == pick_node:
                    lhs_origin += self.x_vars[key]
            self.model += (lhs_origin == 1, f'{pick_node}_pick_cons')

    # (3) 保证取货后要有对应的送货：车辆k经过pick点也一定也要经过对应的del点
    for veh_k, veh_obj in self.system.vehicle_obj_dict.items():
        for node in veh_obj.alter_node_list:
            if node % 2 != 0:
                pick_node, del_node = node, node + 1
                lhs_origin, rhs = 0, 0
                for key in self.x_vars.keys():
                    k, i, j = key[0], key[1], key[2]
                    if k == veh_k and i == pick_node:
                        lhs_origin += self.x_vars[key]
                    if k == veh_k and i == del_node:
                        rhs += self.x_vars[key]
                self.model += (lhs_origin == rhs, f'{veh_k}_{pick_node}_{del_node}_pair_cons')

    # (4) 路径平衡约束：每辆车一定经过从起点出发和回到终点，其他虚拟点要保证有进有出
    for veh_k, veh_obj in self.system.vehicle_obj_dict.items():
        lhs_origin, lhs_des = 0, 0
        for node in veh_obj.alter_node_list:
            lsh, rhs = 0, 0
            for key in self.x_vars.keys():
                k, i, j = key[0], key[1], key[2]
                if k == veh_k and j == node and i == 0:
                    lhs_origin += self.x_vars[key]
                if k == veh_k and i == node and j == (self.system.dummy_node_count + 1):
                    lhs_des += self.x_vars[key]
                if k == veh_k and i == node:
                    lsh += self.x_vars[key]
                if k == veh_k and j == node:
                    rhs += self.x_vars[key]
            self.model += (lsh == rhs, f'{veh_k}_{node}_inout_cons')
        self.model += (lhs_origin == 1, f'{veh_k}_origin_cons')
        self.model += (lhs_des == 1, f'{veh_k}_des_cons')
  • 添加节点变量：上述约束和变量没有反应取送货信息，因此还需要对节点添加一些辅助信息，（1）每辆车到达某虚拟点的载货量（2）每辆车到达虚拟点的时间（3）每辆车在虚拟点的等待时间
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

def add_node_var(self):
    # 变量q_ki：车辆k到达点i的载货量
    # 变量a_ki：车辆k到达点i的时间
    # 变量w_ki：车辆k在点i的等待时间
    for veh_k, veh_obj in self.system.vehicle_obj_dict.items():
        for node_i in veh_obj.alter_node_list:
            self.q_vars[(veh_k, node_i)] = LpVariable(name=f'q_{veh_k}_{node_i}', lowBound=0, cat=LpContinuous)
            self.a_vars[(veh_k, node_i)] = LpVariable(name=f'a_{veh_k}_{node_i}', lowBound=0, cat=LpContinuous)
            self.w_vars[(veh_k, node_i)] = LpVariable(name=f'w_{veh_k}_{node_i}', lowBound=0, cat=LpContinuous)
        node_origin, node_des = 0, self.system.dummy_node_count + 1
        self.q_vars[(veh_k, node_origin)] = LpVariable(name=f'q_{veh_k}_{node_origin}', lowBound=0, cat=LpContinuous)
        self.a_vars[(veh_k, node_origin)] = LpVariable(name=f'a_{veh_k}_{node_origin}', lowBound=0, cat=LpContinuous)
        self.w_vars[(veh_k, node_origin)] = LpVariable(name=f'w_{veh_k}_{node_origin}', lowBound=0, cat=LpContinuous)

        self.q_vars[(veh_k, node_des)] = LpVariable(name=f'q_{veh_k}_{node_des}', lowBound=0, cat=LpContinuous)
        self.a_vars[(veh_k, node_des)] = LpVariable(name=f'a_{veh_k}_{node_des}', lowBound=0, cat=LpContinuous)
        self.w_vars[(veh_k, node_des)] = LpVariable(name=f'w_{veh_k}_{node_des}', lowBound=0, cat=LpContinuous)
  • 添加负载约束：（1）每辆车经过某个虚拟点，其负载量要加上该节点的取/送货量，因此如果经过边(i,j)，那么到达i的负载量+节点j的负载量=到达j的负载量；（2）每辆车到达任何虚拟点，其负载量都不能超过该车的最大载量
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

def add_load_cons(self):
    # (5) 负载平衡约束：如果车辆k经过边(i,j)，那么到达i的负载量+节点j的负载量=到达j的负载量
    for key in self.x_vars.keys():
        k, i, j = key[0], key[1], key[2]
        veh_obj = self.system.vehicle_obj_dict[k]
        big_M = 2 * veh_obj.max_load
        lhs = self.q_vars[(k, j)]
        if i == 0:
            rhs = self.q_vars[(k, i)] + (self.x_vars[key] - 1) * big_M
        else:
            load_i = self.system.dummy_node_dict[i][2]
            rhs = self.q_vars[(k, i)] + load_i + (self.x_vars[key] - 1) * big_M
        # 为了将约束线性化，引入bigM，必须要用不等式，导致该约束限制不了环路的出现
        self.model += (lhs >= rhs, f'{key}_load_balance_cons')

    # (6) 负载最大约束：车辆k到达节点i的负载量不能大于最大载量
    for key in self.q_vars.keys():
        k, i = key[0], key[1]
        veh_obj = self.system.vehicle_obj_dict[k]
        self.model += (self.q_vars[key] <= veh_obj.max_load, f'{k}_{i}_max_load_cons')
  • 添加时间约束：（1）时间窗约束：每辆车在虚拟点开始服务的时间要在时间窗内（2）时间关系约束：如果车辆k经过边(i,j)，那么需要满足到达i的时间+等待时间+服务时间+运输时间不能超过到达j的时间，通过时间关系约束可以限制环路
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44

def add_time_cons(self):
    # (7) 时间平衡约束：如果车辆k经过ij，那么到达j的时间=到达i的时间+在i等待的时间+服务i的时间+运输耗时
    # 可以限制环路
    for key in self.x_vars.keys():
        k, i, j = key[0], key[1], key[2]
        veh_obj = self.system.vehicle_obj_dict[k]
        big_M = 2 * veh_obj.end
        if j == self.system.dummy_node_count + 1:
            loca_j_name = veh_obj.des
        else:
            loca_j_name = self.system.dummy_node_dict[j][1]
        if i == 0:
            loca_i_name = veh_obj.origin
            service_i = 0
        else:
            loca_i_name = self.system.dummy_node_dict[i][1]
            order_id = self.system.dummy_node_dict[i][0]
            order_obj = self.system.order_obj_dict[order_id]
            service_i = order_obj.pick_service
        trans_time = round((self.system.dist_matrix[(loca_i_name, loca_j_name)] / veh_obj.speed) * 3600)
        lhs = self.a_vars[(k, i)] + self.w_vars[(k, i)] + trans_time + (self.x_vars[key] - 1) * big_M + service_i
        self.model += (lhs <= self.a_vars[(k, j)], f'{k}_{i}_{j}_time_balance_cons')

    # (8)(9)(10) 时间窗约束
    # 车辆k在节点i的等待时间=max(0, 节点i可开始服务的时间-到达节点i的时间）
    # 节点i最早开始服务 <= 在节点i的等待时间 + 到达节点i的时间 <= 节点i最晚开始服务
    # 起点的到达时间>=车辆的最早开始时间，终点的到达时间<=车辆的最晚时间
    # 约束2可以实现约束13，去掉约束1也能避免线性化操作
    for key in self.a_vars.keys():
        k, i = key[0], key[1]
        veh_obj = self.system.vehicle_obj_dict[k]
        lhs = self.a_vars[key] + self.w_vars[key]
        if i == 0 or i == self.system.dummy_node_count + 1:
            self.model += (lhs <= veh_obj.end, f'{k}_{i}_window_end_cons')
            self.model += (lhs >= veh_obj.start, f'{k}_{i}_window_start_cons')
        else:
            order_id = self.system.dummy_node_dict[i][0]
            order_obj = self.system.order_obj_dict[order_id]
            if i % 2 != 0:
                self.model += (lhs <= order_obj.pick_end, f'{k}_{i}_window_end_cons')
                self.model += (lhs >= order_obj.pick_start, f'{k}_{i}_window_start_cons')
            else:
                self.model += (lhs <= order_obj.del_end, f'{k}_{i}_window_end_cons')
                self.model += (lhs >= order_obj.del_start, f'{k}_{i}_window_start_cons')
  • 添加顺序约束：添加基本约束后可以得到带环路的图，添加负载约束只能部分限制取送顺序，即不超过最大载量即可，添加时间约束后可以消除环路，得到每辆车从起始虚拟点到结束虚拟点的路径，但是存在先送货再取货的情况，因此需要限制车辆到达取货虚拟点和送货虚拟点的顺序，可以通过时间来约束，到达取货虚拟点时间+等待时间+服务时间+取送运输时间不能大于到达送货虚拟点时间
1
2
3
4
5
6
7
8
9
10
11
12
13

def add_seq_cons(self):
    # (13) 保证货物先取后送
    for key in self.a_vars.keys():
        k, i = key[0], key[1]
        veh_obj = self.system.vehicle_obj_dict[k]
        if i % 2 != 0 and i != 0 and i != (self.system.dummy_node_count + 1):
            node_pick, node_del = i, i + 1
            order_id, loca_i_name, _ = self.system.dummy_node_dict[node_pick]
            loca_j_name = self.system.dummy_node_dict[node_del][1]
            order_obj = self.system.order_obj_dict[order_id]
            trans_time = round((self.system.dist_matrix[(loca_i_name, loca_j_name)] / veh_obj.speed) * 3600)
            lhs = self.a_vars[(k, node_pick)] + self.w_vars[(k, node_pick)] + order_obj.pick_service + trans_time
            self.model += (lhs <= self.a_vars[(k, node_del)], f'{k}_{i}_{i + 1}_seq_cons')
  • 添加限制约束：上述约束添加后已经可以得到满足要求的车辆调度计划，还可以添加额外的限制约束，（1）最大运输距离限制：计算每辆车的运输总距离，不能超过最大值（2）最大运输时间限制：计算每辆车到达结束虚拟点的时间，不能超过最大时间
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

def add_limit_cons(self):
    # (11) 最大运输距离限制：车辆k经过的ij总运输距离不能超过车辆最大距离
    for veh_k, veh_obj in self.system.vehicle_obj_dict.items():
        lhs = 0
        for key in self.x_vars.keys():
            k, i, j = key[0], key[1], key[2]
            if k == veh_k:
                if j == self.system.dummy_node_count + 1:
                    loca_j_name = veh_obj.des
                else:
                    loca_j_name = self.system.dummy_node_dict[j][1]
                if i == 0:
                    loca_i_name = veh_obj.origin
                else:
                    loca_i_name = self.system.dummy_node_dict[i][1]
                lhs += self.x_vars[key] * self.system.dist_matrix[(loca_i_name, loca_j_name)]
        self.model += (lhs <= veh_obj.max_distance, f'{veh_k}_distance_limit_cons')

    # (12) 最大运输时间限制：车辆k到达终点的时间-车辆从起点出发的时间不能超过最大时长
    for veh_k, veh_obj in self.system.vehicle_obj_dict.items():
        time_delta = self.a_vars[(veh_k, self.system.dummy_node_count + 1)] - self.a_vars[(veh_k, 0)]
        self.model += (time_delta <= (veh_obj.max_time * 3600), f'{veh_k}_time_limit_cons')
  • 添加目标函数：用每辆车的单位运输成本乘上运输距离即可
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

def add_trans_cost_objs(self):
    # 对于变量kij，找到ij的距离和k的单位运输成本，即可计算总运输成本
    for key in self.x_vars.keys():
        k, i, j = key[0], key[1], key[2]
        veh_obj = self.system.vehicle_obj_dict[k]
        if i == 0:
            loca_i_name = veh_obj.origin
        else:
            loca_i_name = self.system.dummy_node_dict[i][1]
        if j == len(self.system.dummy_node_dict) + 1:
            loca_j_name = veh_obj.des
        else:
            loca_j_name = self.system.dummy_node_dict[j][1]
        if loca_i_name != loca_j_name:
            dist = self.system.dist_matrix[(loca_i_name, loca_j_name)]
            self.objs += self.x_vars[key] * veh_obj.unit_cost * dist
    self.model += (self.objs, f'minimize_transport_cost')

5.求解结果

求解参数设置：NoRelHeurTime是一种启发式算法，可以提前找到一个可行解，对于大规模混合整数线性规划问题，提前拿到第一个可行解是非常非常有助于整个Branch and Cut的搜索过程的

1
2

NoRelHeurTime = 400  # 控制启发式算法的时间限制
timeLimit = 600  # 控制求解器求解时间限制

1000s后找到gap为45%的最优解，与官方给出代码运行效果基本一致，求解的调度计划也一致

1
2
3
4
5
6
7

Explored 664725 nodes (14855641 simplex iterations) in 600.01 seconds (680.73 work units)
Thread count was 8 (of 8 available processors)

Solution count 10: 846.682 846.682 846.682 ... 876.259

Time limit reached
Best objective 8.466817613171e+02, best bound 4.627762720459e+02, gap 45.3424%

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

车辆id,停靠点,订单id,到达时间,离开时间,地点类型,到达载货量,载货变动量
V1,L1,,2017-01-06 07:00:00,2017-01-06 07:00:00,起点,0.0,
V1,L1,O6,2017-01-06 07:00:00,2017-01-06 08:00:00,取货,0.0,20.0
V1,L7,O6,2017-01-06 08:12:26,2017-01-06 09:12:26,送货,20.0,-20.0
V1,L1,O5,2017-01-06 09:24:52,2017-01-06 10:24:52,取货,0.0,10.0
V1,L1,O2,2017-01-06 10:24:52,2017-01-06 11:24:52,取货,10.0,4.0
V1,L1,O7,2017-01-06 11:24:52,2017-01-06 12:24:52,取货,14.0,10.0
V1,L6,O5,2017-01-06 12:44:53,2017-01-06 13:44:53,送货,24.0,-10.0
V1,L8,O7,2017-01-06 13:49:52,2017-01-06 14:49:52,送货,14.0,-10.0
V1,L3,O2,2017-01-06 15:00:00,2017-01-06 16:00:00,送货,4.0,-4.0
V1,L1,,2017-01-06 16:26:36,2017-01-06 16:26:36,终点,0.0,
V2,L1,,2017-01-06 07:00:00,2017-01-06 07:00:00,起点,0.0,
V2,L1,O1,2017-01-06 07:00:07,2017-01-06 08:00:07,取货,0.0,20.0
V2,L2,O1,2017-01-06 08:15:40,2017-01-06 09:15:40,送货,20.0,-20.0
V2,L2,O10,2017-01-06 09:15:40,2017-01-06 10:15:40,取货,0.0,2.0
V2,L1,O10,2017-01-06 10:31:13,2017-01-06 11:31:13,送货,2.0,-2.0
V2,L1,O9,2017-01-06 11:31:13,2017-01-06 12:31:13,取货,0.0,10.0
V2,L1,O3,2017-01-06 12:31:13,2017-01-06 13:31:13,取货,10.0,10.0
V2,L10,O9,2017-01-06 13:56:10,2017-01-06 14:56:10,送货,20.0,-10.0
V2,L4,O3,2017-01-06 15:00:00,2017-01-06 16:00:00,送货,10.0,-10.0
V2,L1,,2017-01-06 16:26:36,2017-01-06 16:26:36,终点,0.0,
V3,L1,,2017-01-06 07:00:00,2017-01-06 07:00:00,起点,0.0,
V3,L1,O4,2017-01-06 07:00:00,2017-01-06 08:00:00,取货,0.0,10.0
V3,L5,O4,2017-01-06 08:44:45,2017-01-06 09:44:45,送货,10.0,-10.0
V3,L1,O8,2017-01-06 10:29:30,2017-01-06 11:29:30,取货,0.0,10.0
V3,L9,O8,2017-01-06 11:52:29,2017-01-06 12:52:29,送货,10.0,-10.0
V3,L1,,2017-01-06 16:44:45,2017-01-06 16:44:45,终点,0.0,

PDPTW在如此小规模问题下都需要较长求解时间，求解时间是比较长的，在实际应用中，车辆调度次数比较频繁，约束更复杂，因此往往采用启发式算法来实现，能比较快的找到一个次优解。

运筹优化, GurobiCase
OR
This post is licensed under CC BY 4.0 by the author.