vrp python vrp Python包

1. 介绍

基于NetworkX包（操作图）和cspy包（启发式算法）开发，虽然性能不是很好，但是功能强大，易于上手，可以用来学习column generation求解VRP问题的写法。

支持如下类型的问题：
the Capacitated VRP (CVRP),
the CVRP with resource constraints,
the CVRP with time windows (CVRPTW),
the CVRP with simultaneous distribution and collection (CVRPSDC),
the CVRP with heterogeneous fleet (HFCVRP).

1.1 性能

性能可以参考官方文档的这张图，和google ortools进行的对比：

1.2 简单例子

from networkx import DiGraph
G = DiGraph()
for v in [1, 2, 3, 4, 5]:
    G.add_edge("Source", v, cost=10)
    G.add_edge(v, "Sink", cost=10)
G.add_edge(1, 2, cost=10)
G.add_edge(2, 3, cost=10)
G.add_edge(3, 4, cost=15)
G.add_edge(4, 5, cost=10)
from vrpy import VehicleRoutingProblem
prob = VehicleRoutingProblem(G)
prob.num_stops = 3
prob.solve(time_limit=10)

结果为

1.3 加上capacity约束

for v in G.nodes():
    if v not in ["Source", "Sink"]:
        G.nodes[v]["demand"] = 5
prob.load_capacity = 10
prob.solve()
prob.best_routes

1.4 加上时间约束

for (u, v) in G.edges():
    G.edges[u,v]["time"] = 20
G.edges[4,5]["time"] = 25
prob.duration = 60
prob.solve()
prob.best_value

1.5 加上time window约束

time_windows = {1: (5, 100), 2: (5, 20), 3: (5, 100), 4: (5, 100), 5: (5, 100),"Sink":(5,110),"Source":(5,100)}
for v in G.nodes():
    G.nodes[v]["lower"] = time_windows[v][0]
    G.nodes[v]["upper"] = time_windows[v][1]
    if v not in ["Source","Sink"]:
        G.nodes[v]["service_time"] = 1        
from vrpy import VehicleRoutingProblem
prob = VehicleRoutingProblem(G)
prob.time_windows = True
prob.num_stops = 3
prob.solve()
prob.best_routes

1.6 混合fleet

两辆车的 load_capacity 分别是5和20；
fixed_cost 分别是0和5；
travel costs，第二辆车比第一辆车贵1；

from networkx import DiGraph
from vrpy import VehicleRoutingProblem
G = DiGraph()
for v in [1, 2, 3, 4, 5]:
    G.add_edge("Source", v, cost=[10, 11])
    G.add_edge(v, "Sink", cost=[10, 11])
    G.nodes[v]["demand"] = 5
G.add_edge(1, 2, cost=[10, 11])
G.add_edge(2, 3, cost=[10, 11])
G.add_edge(3, 4, cost=[15, 16])
G.add_edge(4, 5, cost=[10, 11])
prob=VehicleRoutingProblem(G, mixed_fleet=True, fixed_cost=[0, 5], load_capacity=[5, 20])
prob.solve()
prob.best_routes

1.7 综合例子

import networkx as nx
from vrpy import VehicleRoutingProblem

# Create graph
G = nx.DiGraph()
for v in [1, 2, 3, 4, 5]:
       G.add_edge("Source", v, cost=10, time=20)
   G.add_edge(v, "Sink", cost=10, time=20)
   G.nodes[v]["demand"] = 5
   G.nodes[v]["upper"] = 100
   G.nodes[v]["lower"] = 5
   G.nodes[v]["service_time"] = 1
G.nodes[2]["upper"] = 20
G.nodes["Sink"]["upper"] = 110
G.nodes["Source"]["upper"] = 100
G.add_edge(1, 2, cost=10, time=20)
G.add_edge(2, 3, cost=10, time=20)
G.add_edge(3, 4, cost=15, time=20)
G.add_edge(4, 5, cost=10, time=25)

# Create vrp
prob = VehicleRoutingProblem(G, num_stops=3, load_capacity=10, duration=64, time_windows=True)

# Solve and display solution
prob.solve()
print(prob.best_routes)
print(prob.best_value)

下面是ortools的一个例子：

from networkx import DiGraph, from_numpy_matrix, relabel_nodes, set_node_attributes
from numpy import array

# Distance matrix
DISTANCES = [
[0,548,776,696,582,274,502,194,308,194,536,502,388,354,468,776,662,0], # from Source
[0,0,684,308,194,502,730,354,696,742,1084,594,480,674,1016,868,1210,548],
[0,684,0,992,878,502,274,810,468,742,400,1278,1164,1130,788,1552,754,776],
[0,308,992,0,114,650,878,502,844,890,1232,514,628,822,1164,560,1358,696],
[0,194,878,114,0,536,764,388,730,776,1118,400,514,708,1050,674,1244,582],
[0,502,502,650,536,0,228,308,194,240,582,776,662,628,514,1050,708,274],
[0,730,274,878,764,228,0,536,194,468,354,1004,890,856,514,1278,480,502],
[0,354,810,502,388,308,536,0,342,388,730,468,354,320,662,742,856,194],
[0,696,468,844,730,194,194,342,0,274,388,810,696,662,320,1084,514,308],
[0,742,742,890,776,240,468,388,274,0,342,536,422,388,274,810,468,194],
[0,1084,400,1232,1118,582,354,730,388,342,0,878,764,730,388,1152,354,536],
[0,594,1278,514,400,776,1004,468,810,536,878,0,114,308,650,274,844,502],
[0,480,1164,628,514,662,890,354,696,422,764,114,0,194,536,388,730,388],
[0,674,1130,822,708,628,856,320,662,388,730,308,194,0,342,422,536,354],
[0,1016,788,1164,1050,514,514,662,320,274,388,650,536,342,0,764,194,468],
[0,868,1552,560,674,1050,1278,742,1084,810,1152,274,388,422,764,0,798,776],
[0,1210,754,1358,1244,708,480,856,514,468,354,844,730,536,194,798,0,662],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0], # from Sink
]

# Demands (key: node, value: amount)
DEMAND = {1: 1, 2: 1, 3: 2, 4: 4, 5: 2, 6: 4, 7: 8, 8: 8, 9: 1, 10: 2, 11: 1, 12: 2, 13: 4, 14: 4, 15: 8, 16: 8}

# The matrix is transformed into a DiGraph
A = array(DISTANCES, dtype=[("cost", int)])
G = from_numpy_matrix(A, create_using=nx.DiGraph())

# The demands are stored as node attributes
set_node_attributes(G, values=DEMAND, name="demand")

# The depot is relabeled as Source and Sink
G = relabel_nodes(G, {0: "Source", 17: "Sink"})

>>> from vrpy import VehicleRoutingProblem
>>> prob = VehicleRoutingProblem(G, load_capacity=15)
>>> prob.solve()
>>> prob.best_value
6208.0
>>> prob.best_routes
{1: ['Source', 12, 11, 15, 13, 'Sink'], 2: ['Source', 1, 3, 4, 7, 'Sink'], 3: ['Source', 5, 2, 6, 8, 'Sink'], 4: ['Source', 14, 16, 10, 9, 'Sink']}
>>> prob.best_routes_load
{1: 15, 2: 15, 3: 15, 4: 15}

2. API参考

2.1 VehicleRoutingProblem类

下面是系统参数，包含了所有必须的数据

G (DiGraph) – The underlying network.

num_stops (int, optional) – Maximum number of stops. Defaults to None.

load_capacity (list, optional) – Maximum load per vehicle. Each item of the list points to a different capacity. Defaults to None.

duration (int, optional) – Maximum duration of route. Defaults to None.

time_windows (bool, optional) – True if time windows on vertices. Defaults to False.

pickup_delivery (bool, optional) – True if pickup and delivery constraints. Defaults to False.

distribution_collection (bool, optional) – True if distribution and collection are simultaneously enforced. Defaults to False.

drop_penalty (int, optional) – Value of penalty if node is dropped. Defaults to None.

fixed_cost (int, optional) – Fixed cost per vehicle. Defaults to 0.

num_vehicles (int, optional) – Maximum number of vehicles available. Defaults to None (in this case num_vehicles is unbounded).

periodic (int, optional) – Time span if vertices are to be visited periodically. Defaults to None.

mixed_fleet (bool, optional) – True if heterogeneous fleet. Defaults to False.

minimize_global_span (bool, optional) – True if global span is minimized (instead of total cost). Defaults to False.

2.2 solve方法

包含了求解所需要的配置

initial_routes (list, optional) – List of routes (ordered list of nodes). Feasible solution for first iteration. Defaults to None.

preassignments (list, optional) – List of preassigned routes (ordered list of nodes). If the route contains Source and Sink nodes, it is locked, otherwise it may be extended. Defaults to None.

pricing_strategy (str, optional) –

Strategy used for solving the sub problem. Options available :

”Exact”: the subproblem is solved exactly;

”BestEdges1”: some edges are removed;

”BestEdges2”: some edges are removed (with a different strategy);

”BestPaths”: some edges are removed (with a different strategy);

”Hyper”: choose from the above using a hyper_heuristic (see hyper_heuristic.py);

Defaults to “BestEdges1”.

cspy (bool, optional) – True if cspy is used for subproblem. Defaults to True.

exact (bool, optional) – True if only cspy’s exact algorithm is used to generate columns. Otherwise, heuristics will be used until they produce +ve reduced cost columns, after which the exact algorithm is used. Defaults to True.

time_limit (int, optional) – Maximum number of seconds allowed for solving (for finding columns). Defaults to None.

solver (str, optional) – Solver used. Three options available: “cbc”, “cplex”, “gurobi”. Using “cplex” or “gurobi” requires installation. Not available by default. Additionally, “gurobi” requires pulp to be installed from source. Defaults to “cbc”, available by default.

dive (bool, optional) – True if diving heuristic is used. Defaults to False.

greedy (bool, optional) – True if randomized greedy algorithm is used to generate extra columns. Only valid for capacity constraints, time constraints, num stops constraints. Defaults to False.

max_iter (int, optional) – maximum number of iterations for the column generation procedure.

run_exact (int, optional) – if a pricing strategy is selected, this parameter controls the number of iterations after which the exact algorithm is run. Defaults to 1.

3. 列生成算法介绍

首先要理解对偶变量和marginal cost的意义

使用set covering对VRP问题建模，$\Omega$为可行解的集合：

举个只有3个点的例子，距离分别是0-1(3)，0-2(4)，1-2(5)。

首先枚举所有可能的路径（注意这里是有顺序的，在time windows等约束条件下，访问的先后顺序是有影响的）：

$\Omega$={0-1-0(6); 0-2-0(8); 0-1-2-0(12);0-2-1-0(12)}

则建模为：
$\min 6\lambda_1+8\lambda_2+12\lambda_{3}+12\lambda_{4}$
s.t.
$\lambda_1+\lambda_3+\lambda_4=1$（对应点1）
$\lambda_2+\lambda_3+\lambda_4=1$（对应点2）
另外再加上0-1变量约束。
每一个约束对应一个对偶变量称为$\pi_v$，变量$\lambda_i$的marginal cost为：$\hat{c}_i=c_i-\Sigma\pi_va_{v,i}$。

3.1 限制主问题

实际求解的时候，是用启发式算法生成$\Omega$的一个子集，称为“限制主问题”，然后不断求解子问题进行补充，称为列生成算法。以上面的问题为例，假设我们的初始问题为：
$\min 6\lambda_1+8\lambda_2$
s.t.
$\lambda_1=1$（对应点1）
$\lambda_2=1$（对应点2）

我们再生成一个$\lambda_3$对应的轨迹
$\min 6\lambda_1+8\lambda_2+12\lambda_{3}$
s.t.
$\lambda_1+\lambda_3=1$（对应点1）
$\lambda_2+\lambda_3=1$（对应点2）
初始限制主问题的解为$(1,1,0)$，对偶解为$(6,8)$。$\lambda_1$和$\lambda_2$是基变量，残差为0；$\lambda_3$为非基变量，残差为$12-6*1-8*1=-2<0$，因此可以把$\lambda_3$转入进基变量。
剩下的问题就是，如何不断迭代生成新的轨迹。

3.2 子问题

3.2.1 精确算法

生成新的“列”也是一个线性规划问题，可以以最小化残差为目标来进行求解。令$x_{uv}$为边$u-v$是否连接，则子问题是如下最短路问题：

生成$\lambda_3$路径的子问题为：

$\min 3x_{01}+3x_{10}+4x_{02}+4x_{20}+5x_{12}+5x_{21}-6*(x_{01}+x_{10}+x_{12}+x_{21})-8*(x_{02}+x_{20}+x_{12}+x_{21})$
$x_{01}+x_{21}=x_{10}+x_{12}=1$
$x_{02}+x_{12}=x_{20}+x_{21}=1$
$x_{02}+x_{01}=x_{10}+x_{20}=1$
使用julia进行求解：

using JuMP
using Cbc
prob = Model(Cbc.Optimizer)
@variable(prob, x[i=0:2, j=0:2], Bin)
for i = 0:2
    @constraint(prob, sum(x[i,j] for j=0:2) == 1)
    @constraint(prob, sum(x[j,i] for j=0:2) == 1)
end
@objective(prob,Min,3x[0,1]+3x[1,0]+4x[0,2]+4x[2,0]+5x[1,2]+5x[2,1]-6*(x[0,1]+x[1,0]+x[1,2]+x[2,1])-8*(x[0,2]+x[2,0]+x[1,2]+x[2,1]))

optimize!(prob)
# Extract the values of x
x_val = value.(x)

结果为：

包里面的代码如下：

def _formulate(self):
        # minimize reduced cost
        self.prob += pulp.lpSum(
            [
                self.sub_G.edges[i, j]["weight"] * self.x[(i, j)]
                for (i, j) in self.sub_G.edges()
            ]
        )
        # flow balance
        for v in self.sub_G.nodes():
            if v not in ["Source", "Sink"]:
                in_flow = pulp.lpSum(
                    [self.x[(i, v)] for i in self.sub_G.predecessors(v)]
                )
                out_flow = pulp.lpSum(
                    [self.x[(v, j)] for j in self.sub_G.successors(v)]
                )
                self.prob += in_flow == out_flow, "flow_balance_%s" % v

        # Start at Source and end at Sink
        self.prob += (
            pulp.lpSum([self.x[("Source", v)] for v in self.sub_G.successors("Source")])
            == 1,
            "start_at_source",
        )
        self.prob += (
            pulp.lpSum([self.x[(u, "Sink")] for u in self.sub_G.predecessors("Sink")])
            == 1,
            "end_at_sink",
        )
        # Forbid route Source-Sink
        self.prob += (
            pulp.lpSum([self.x[(i, j)] for (i, j) in self.sub_G.edges()]) >= 2,
            "at_least_1_stop",
        )

        # Problem specific constraints
        if self.time_windows:
            self._add_time_windows()
        if self.num_stops:
            self._add_max_stops()
        if self.load_capacity:
            self._add_max_load()
        if self.duration:
            self._add_max_duration()
        self.elementarity = False
        if negative_edge_cycle(self.G):
            logger.debug("negative cycle found")
            self._add_elementarity()
            self.elementarity = True
        if self.pickup_delivery:
            self._add_pickup_delivery()
        if self.distribution_collection:
            self._add_distribution_collection()

3.2.2 简化算法

使用edge remove方法移除一部分边，减少子问题求解的计算量。vrpy中有三种remove方法：
1）移除条件$c_{ij}>\alpha \pi_{max}$，代码如下，有一个入参alpha。逐步增大alpha，直到找到一个合适的路径。

def remove_edges_1(self, alpha):
    """
    Removes edges based on criteria described here :
    https://pubsonline.informs.org/doi/10.1287/trsc.1050.0118

    Edges for which [cost > alpha x largest dual value] are removed,
    where 0 < alpha < 1 is a parameter.
    """
    self.sub_G = self.G.copy()
    largest_dual = max(
        self.duals[v] for v in self.duals if v != "upper_bound_vehicles"
    )

    for (u, v) in self.G.edges():
        if self.G.edges[u, v]["cost"][self.vehicle_type] > alpha * largest_dual:
            self.sub_G.remove_edge(u, v)
    # If pruning the graph disconnects the source and the sink,
    # do not solve the subproblem.
    try:
        if not has_path(self.sub_G, "Source", "Sink"):
            self.run_subsolve = False
    except NetworkXException:
        self.run_subsolve = False

2）移除reduce cost最大的一部分边（比例为ratio）。边的reduce cost由边的cost减去那一列的dual values*那一列是否包含这条边 得到，计算逻辑如下：

def add_reduced_cost_attribute(self):
     """Substracts the dual values to compute reduced cost on each edge."""
     for edge in self.G.edges(data=True):
         edge[2]["weight"] = edge[2]["cost"][self.vehicle_type]
         if self.route:
             edge[2]["weight"] *= -self.duals[
                 "makespan_%s" % self.route.graph["name"]
             ]
         for v in self.duals:
             if edge[0] == v:
                 edge[2]["weight"] -= self.duals[v]
     if "upper_bound_vehicles" in self.duals:
         for v in self.G.successors("Source"):
             self.G.edges["Source", v]["weight"] -= self.duals[
                 "upper_bound_vehicles"
             ][self.vehicle_type]

移除边的代码为：

def remove_edges_2(self, ratio):
    """
    Removes edges based on criteria described here :
    https://www.sciencedirect.com/science/article/abs/pii/S0377221717306045

    Edges are sorted by non decreasing reduced cost, and only
    the K|E| ones with lowest reduced cost are kept, where K is a parameter (ratio).
    """
    self.sub_G = self.G.copy()
    # Sort the edges by non decreasing reduced cost
    reduced_cost = {}
    for (u, v) in self.G.edges():
        if u != "Source" and v != "Sink":
            reduced_cost[(u, v)] = self.G.edges[u, v]["weight"]
    sorted_edges = sorted(reduced_cost, key=reduced_cost.get)
    # Keep the best ones
    limit = int(ratio * len(sorted_edges))
    self.sub_G.remove_edges_from(sorted_edges[limit:])
    # If pruning the graph disconnects the source and the sink,
    # do not solve the subproblem.
    try:
        if not has_path(self.sub_G, "Source", "Sink"):
            self.run_subsolve = False
    except NetworkXException:
        self.run_subsolve = False

3）不考虑资源约束，去除路径最短的k个路线，代码如下：

def remove_edges_3(self, beta):
     """
     Heuristic pruning:
     1. Normalize weights in interval [-1,1]
     2. Set all negative weights to 0
     3. Compute beta shortest paths (beta is a paramater)
        https://networkx.github.io/documentation/networkx-1.10/reference/generated/networkx.algorithms.simple_paths.shortest_simple_paths.html
     4. Remove all edges that do not belong to these paths
     """
     # Normalize weights
     max_weight = max([self.G.edges[i, j]["weight"] for (i, j) in self.G.edges()])
     min_weight = min(self.G.edges[i, j]["weight"] for (i, j) in self.G.edges())
     for edge in self.G.edges(data=True):
         edge[2]["pos_weight"] = (
             -max_weight - min_weight + 2 * edge[2]["weight"]
         ) / (max_weight - min_weight)
         edge[2]["pos_weight"] = max(0, edge[2]["pos_weight"])
     # Compute beta shortest paths
     best_paths = list(
         islice(
             shortest_simple_paths(self.G, "Source", "Sink", weight="pos_weight"),
             beta,
         )
     )
     # Store these paths as a list of DiGraphs
     best_paths_list = []
     for path in best_paths:
         H = DiGraph()
         add_path(H, path)
         best_paths_list.append(H)
     # Merge the paths into one graph
     induced_graph = compose_all(best_paths_list)
     # Create subgraph induced by the edges of this graph
     self.sub_G = self.G.edge_subgraph(induced_graph.edges()).copy()

3.2.3 使用cspy

参考这篇文章：
。
cspy是一个有约束条件下生成最短路的一个算法包，它是好几个列生成开源算法包中用于列生成算法的基础算法包，主要是用于生成新的列。最开始用于公交司机调度问题，后面被广泛用于交通相关的问题中。
cspy是基于networkx的，算法清单包括：

• BiDirectional: Bidirectional and monodirectional algorithms
• Tabu Heuristic Search
• Greedy Elimination Procedure
• GRASP
• Particle Swarm Optimization with combined Local and Global Expanding Neighbourhood Topology (PSOLGENT)

3.2.4 使用greedy算法

思路及其简单：

We can attempt to find negative-reduced-cost elementary path with a simple greedy algorithm that builds a path starting in σ and then proceeds by choosing one random arc among the K1 outgoing arcs of least reduced cost that do not close a cycle, until τ is reached. The same procedure can also be applied backward, starting in τ and ending at σ. After the path is constructed, we can check that the capacity constraint has not been violated by a single pass over the list of visited ports; if the path is not capacity-feasible, or if its reduced cost is non-negative, it is discarded. The algorithm can be applied in both directions many times,
as its computational time is very small, thereby increasing the chances that a feasible path of negative reduced cost is found.

代码如下：

def solve(self, n_runs=5):
        """The forward and backwards search are run."""
        more_routes = False
        # The forward search is run n_runs times
        for _ in range(n_runs):
            self._initialize_run()
            self.run_forward()
            if self._new_node and self._weight < 0:
                logger.debug("negative column %s" % self._weight)
                more_routes = True
                self._add_new_route()
        # The backwards search is run n_runs times
        for _ in range(n_runs):
            self._initialize_run()
            self.run_backwards()
            if self._new_node and self._weight < 0:
                logger.debug("negative column %s" % self._weight)
                more_routes = True
                self._add_new_route()
        return self.routes, more_routes

    def run_forward(self):
        """
        A path starting from Source is randomly greedily extended
        until Sink is reached.
        The procedure aborts if path becomes infeasible.
        """
        self._current_path = ["Source"]
        extend = True
        new_node = True
        while extend and new_node:
            new_node = self._get_next_node()
            extend = self._update(forward=True)

3.3 代码简单分析

文件夹如下

VehicleRoutingProblem类中包含一个重要的_solve函数，

其中最重要的是column generation函数，迭代生成一系列新的列（find column）

里面迭代代用find columns函数

1. 使用pulp求解master problem（master.lp文件保存在同文件夹下），得到每一个约束对应的dual和relaxed_cost
2. 根据dual求解sub problem，默认使用的是bestEdge1算法
3. 另外再使用greedy算法再生成一批columns，然后更新限制主问题