Python实现动态规划Labeling算法求解SPPRC问题

本文中的课件来自清华大学深圳国际研究生院，物流与交通学部戚铭尧教授《物流地理信息系统》课程。

SPPRC问题

  带资源约束的最短路径问题（shortest path problem with resource constraints）是一个众所周知的NP-Hard问题。 除了作为网络问题直接应用外，SPPRC还用作列生成解决方案方法的基础，用于解决车辆路径规划问题和人员排班问题等。
　　考虑一个有向图$G=(N,A)$， N = { v 1 , v 2 , . . . , v i , . . . , v n } N=\{v_1,v_2,...,v_i,...,v_n\} N={v1​,v2​,...,vi​,...,vn​}表示节点的集合，并且 A = { ( i , j ) ∣ v i ∈ N , v j ∈ N , i ≠ j } A=\{(i,j)|v_i \in N,v_j \in N ,i \ne j\} A={(i,j)∣vi​∈N,vj​∈N,i​=j}表示弧的集合。对于每一段弧$(i,j)\in A$都有一个非负的权重(vrptw问题中可能为负，需要作特殊处理)，$c_{ij}$和$t_{ij}$，表示通过这段弧的成本和资源消耗。
　　SPPRC问题包括找到从起始节点$v_s\in N$到结束节点$v_t\in N$的一条路径$P$，使该路径的总成本最小化，但不超过最大资源消耗$T$。即使只存在一种资源，SPPRC也是一个NP-Hard。

Labeling算法

  直接求解SPPRC问题是比较困难的，故通常采用一种动态规划的方法：$labelcorrectionalgorithm$算法。
通常我们考虑有$L$种资源，可能包括时间、最大装载重量、最大装载体积等。
  对于每一条从起始节点$v_s$扩展到点$v_i$的路径$X_{pi}$都有一个标签$(R_i,C_i)$与之相关联。$R_i\equiv (T_i^1,T_i^2,...,T_i^L)$是路径使用的每L个资源的数量；$C_i$是$cost$。
  此外，我们还需要设置一些优超规则，称为$dominancerule$，以帮助减少不必要的扩展。令$X_{pi}$和$X_{pi}^{\ast }$是从$v_s$扩展到$v_i$的两条不同的路径，与之关联的标签分别为$(R_i,C_i)$和$(R_i^{\ast },C_i^{\ast })$。路径$X_{pi}$可以$dominates$路径$X_{pi}^{\ast }$当且仅当：
$$C_i\le C_i^{\ast },T_i^l\le T_i^{l\ast },\forall l\in (1,...,L),(R_i,C_i)\ne (R_i^{\ast },C_i^{\ast })$$

下面介绍一个简单的小例子：

一个有向图如上图所示，我们可以看到每一段弧上都有一个非负的权重，分别表示通过该弧的$time$与$cost$。首先初始化出发节点$p$的标签，设置为$[0,0]$，随后令其向可达节点$v_2$和$v_3$进行扩展，则分别设置一个新的标签$[2,2]$。再由$v_2$向$v_3$和$d$进行扩展，则$v_3$得到一个新的标签$[3,0]$，$d$得到一个新的标签$[4,4]$，同理继续扩展到$d$，生成一个新的标签$[5,2]$。然后回过头由$v_3$扩展到$v_2$，$v_2$得到一个新的标签$[4,1]$。我们可以发现没有继续向$d$扩展了，因为其生成的标签为$[6,2]$，与$d$已有的标签$[5,2]$相比满足了$dominancerule$的要求，不必再扩展。

Python编程实现

考虑这样一个有向图。

首先在python中对该图进行定义：

import numpy as np
import networkx as nx
from time import *
# 点中包含时间窗属性
Nodes = {'s': (0, 0)
         , '1': (6, 14)
         , '2': (9, 12)
         , '3': (8, 12)
         , 't': (9, 15)
        }
# 弧的属性包括travel_time与distance
Arcs = {('s','1'): (8, 3)
        ,('s','2'): (5, 5)
        ,('s','3'): (12, 2)
        ,('1','t'): (4, 7)
        ,('2','t'): (2, 6)
        ,('3','t'): (4, 3)
       }

# create the directed Graph
Graph = nx.DiGraph()
cnt = 0
# add nodes into the graph
for name in Nodes.keys():
    cnt += 1
    Graph.add_node(name
                   , time_window = (Nodes[name][0], Nodes[name][1])
                   , min_dis = 0  
                   )
# add edges into the graph
for key in Arcs.keys():
    Graph.add_edge(key[0], key[1]
                   , travel_time = Arcs[key][0]
                   , length = Arcs[key][1]
                   )

创建Label类

class Label:
    path = []
    travel_time = 0
    distance = 0

# dominance rule
def dominate(Q, Path_set):
    QCopy = copy.deepcopy(Q)
    PathsCopy = copy.deepcopy(Path_set)
    
    # dominate Q
    for label in QCopy:
        for another_label in Q:
            if (label.path[-1] == another_label.path[
                -1] and label.time < another_label.time and label.dis < another_label.dis):
                Q.remove(another_label)
            print('dominated path (Q) : ', another_label.path)

    # dominate Paths
    for key_1 in PathsCopy.keys():
        for key_2 in PathsCopy.keys():
            if (PathsCopy[key_1].path[-1] == PathsCopy[key_2].path[-1]
                    and PathsCopy[key_1].travel_time < PathsCopy[key_2].travel_time
                    and PathsCopy[key_1].length < PathsCopy[key_2].length
                    and (key_2 in Path_set.keys())):
                Path_set.pop(key_2)
                print('dominated path (P) : ', PathsCopy[key_1].path)
                
    return Q, Path_set
    
# labeling algorithm
def Labeling_SPPRC(Graph, org, des):
    # initial Q
    Q = []
    Path_set = {}

    # creat initial label
    label = Label()
    label.path = [org]
    label.travel_time = 0
    label.distance = 0
    Q.append(label)

    count = 0

    while(len(Q) > 0):
        count += 1
        current_path = Q.pop()

        # extend the current label
        last_node  = current_path.path[-1]
        for child in Graph.successors(last_node):
            extended_path = copy.deepcopy(current_path)
            arc = (last_node, child)

            # check the feasibility
            arrive_time = current_path.travel_time + Graph.edges[arc]['travel_time']
            time_window = Graph.nodes[child]['time_window']
            if(arrive_time >= time_window[0] and arrive_time <= time_window[1] and last_node != des):
                extended_path.path.append(child)
                extended_path.travel_time += Graph.edges[arc]['travel_time']
                extended_path.distance += Graph.edges[arc]['length']
                Q.append(extended_path)

    Path_set[count] = current_path
    # 调用dominance rule
    Q, Path_set = dominance(Q, Path_set)
    
 # filtering Paths, only keep feasible solutions
    Path_set_copy = copy.deepcopy(Path_set)
    for key in Path_set_copy.keys():
        if(Path_set[key].path[-1] != des):
            Path_set.pop(key)

    # choose optimal solution
    opt_path = {}
    min_dis = 1e6
    for key in Path_set.keys():
        if(Path_set[key].distance < min_dis):
            min_dis = Path_set[key].distance
            opt_path[1] = Path_set[key]
            
    return Graph, Q, Path_set, opt_path

调用labeling算法计算

org = 's'
des = 't'
begin_time = time()
Graph, Q, Path_set, opt_path = Labeling_SPPRC(Graph, org, des)
end_time = time()
print("计算时间： ", end_time - begin_time)
print('optimal path : ', opt_path[1].path )
print('optimal path (distance): ', opt_path[1].distance)
print('optimal path (time): ', opt_path[1].travel_time)

最终运行结果如图所示：

我们也可以尝试计算更多的点，并且加入$loadcapacity$的约束，对$Solomonc101.txt$算例的计算结果为：

对$Solomonc1\_10\_1.TXT$算例取300个客户点的计算结果为：

至于500个点往上，我的电脑就跑不出来啦。。。
感谢大家的阅读。

作者：夏旸，清华大学，工业工程系/深圳国际研究生院 (硕士在读)
刘兴禄，清华大学，清华伯克利深圳学院（博士在读）

欢迎大家关注公众号运小筹。