FILO：超大规模车辆路径问题的快速启发式求解

---

本文仅对FILO算法进行简单介绍，同时演示如何利用FILO的算法库求解基于.vrp文件的VRP问题，并对比与具有商用许可的Gurobi的求解效率和效果对比，而算法的具体实现细节以及理论部分，读者可以通过参考资料链接自行了解。

1. 背景引入

车辆路径问题（Vehicle routing problem, VRP） 是经典的组合优化问题，旨在为一组车辆设计最优的配送路径，以满足特定的需求和约束条件。该问题通常涉及多个配送点、车辆的容量限制以及服务时间窗口等因素。目标是最小化总运输成本（如距离、时间或费用），同时确保所有客户的需求都得到满足。车辆路径问题广泛应用于物流、配送、公共交通等领域。而带容量约束的车辆路径问题（Capacitated vehicle routing problem, CVRP） 是在VRP的基础上，为每个客户节点增加一个对应的需求量，同时车辆有一个装车容量，以此限制车辆一趟能满足的需求量，如下图所示，起始点$(x_0,y_0)$处，有无穷多的容量为$Q$的车辆，以及有多个配送点，配送节点$i$的需求量为$q_i$，以最小化总的配送成本（距离）为目标，得到一组车辆的配送路径方案，这是一个典型的CVRP问题示例。

车辆路径是一类NP-Hard问题，以CVRP为例，如果不固定问题的最大求解时间，许多前沿的算法的计算复杂度表现出指数级的增长趋势。

因此，Luca Accorsi和Daniele Vigo决定设计一个快速、高效、可扩展的启发式方法，如下。

2. FILO库的使用步骤

FILO(Fast-ILS Localized Optimization，快速迭代局部优化) 从局部搜索加速技术、解空间剪支技术、参数优化等方面进行精心设计，能够在普通的算力系统下求解超大规模的CVRP问题（尽管超大规模问题的解可能质量较差），例如可以在10GB内存每分钟执行上万次的迭代。其中关键的局部搜索引擎包括：

1. 分层随机变邻域下降（Hierarchical Randomized Variable Neighborhood Descent, HRVND）作为邻域勘探技术；
2. 静态移动描述符（Static Move Descriptors, SMD）用于加速局部搜索的迭代过程，将遍历邻域算子的做法替换为有组织的算子检查；
3. Granular Neighborhoods和Selective Vertex Caching作为解空间的剪支技术，限制局部搜索只会转移到更有希望的空间进行搜索。

这里我们介绍如何使用FILO2（Github地址），这是比FILO性能更好的版本，能够以更短的时间求解更大规模的算例。下文将介绍FILO2库如何求解.vrp文件格式的VRP问题，而对算法感兴趣的读者可以阅读相关的报告¹²及论文³。

（1）获取项目代码

通过git工具将FILO代码及数据下载到本地。整个项目文件约$3G$，其中“result”文件夹的结果数据约占$1.5G$。

git clone git@github.com:acco93/filo2.git

（2）编译构建项目

下载得到一个“filo2”文件夹，由于filo2是由C++实现的，需要通过cmake和make进行编译和构建项目（这是C++项目使用最广泛的构建系统），因此需要在电脑上安装好 cmake（相关安装教程），安装后通过cmake --version检查是否正确安装。

安装好编译器后，依次执行下面的代码：

cd filo2
mkdir build
cd build
cmake .. -DCMAKE_BUILD_TYPE=Release -DENABLE_VERBOSE=1
make -j

前三行代码表示在“filo2”文件夹下创建“build”文件夹，并把工作目录切换到“build”当中。通过cmake命令配置构建环境，这个过程中会自动创建Makefile文件，其中包含了如何如何编译源代码以及链接对象文件生成最终可执行文件的信息，然后使用make命令根据Makefile文件中的规则来编译和构建项目。

其中，cmake后面接的是..，表示从上一层目录当中“filo2”寻找定义项目如何构建的配置文件（通常是名为 CMakeLists.txt 的文件，这是cmake的核心），-DCMAKE_BUILD_TYPE=Release 这个选项设置构建类型为 Release；-DENABLE_VERBOSE=1 表示在程序执行时打印过程信息。

（3）使用命令行调用可执行文件

执行完上面的代码之后，会在“build”文件夹中生成 filo2 可执行文件，需要传入运行数据。如下命令所示，传入“…/instances/X/”文件夹下的X-n129-k18.vrp文件，回车即可求解该例子。

./filo2 ../instances/X/X-n129-k18.vrp

在“instances“文件夹当中有三类数据集，其中，“X”是标准规模数据集；“B”是大规模数据集；“I”是超大规模数据集。

3. 与Gurobi的求解效率对比

这一部分，我们需要对实例数据的输入文件.vrp和输出文件.sol做出简单的解析，并通过实验对比FILO2与具有商业版License的Gurobi的求解性能对比。

CVRP的输入和输出文件的更详细解析可以参考文章：Python调用vrplib库解析具有VRPLIB标准格式的CVRP实例数据，其中提到可以通过vrplib进行自动化解析。

（1）输入文件的内容解析

CVRP算例的标准输入文件以.vrp为后缀（有的数据集以.txt为后缀），其中的数据内容的格式及含义是固定的，以X-n129-k18.vrp文件为例，具体的内容含义如下：

1. NAME: 实例文件名，值为X-n129-k18；
2. COMMENT：数据注释说明；
3. TYPE：数据的问题类型，CVRP类型；
4. DIMENSION：节点数，值为$29$，该值与文件名对应；
5. EDGE_WEIGHT_TYPE：基于节点坐标距离取权值的方式（EUC_2D：双精度的欧几里得距离，不进行四舍五入）；
6. CAPACITY：车辆的容量，值为$39$；
7. NODE_COORD_SECTION：存储节点坐标信息，每行数据的三个值分别表示节点的索引信息、$x$坐标、$y$坐标；
8. DEMAND_SECTION：存储节点的需求信息，没行数据的两个值分别表示节点的索引信息、节点的需求量；
9. DEPOT_SECTION：存储 Depot 节点信息，下述例子表示索引值为$1$的是 Depot；
10. EOF：固定写法，标记数据文件到达结尾。

NAME : 	X-n129-k18	
COMMENT : 	"Generated by Uchoa, Pecin, Pessoa, Poggi, Subramanian, and Vidal (2013)"	
TYPE : 	CVRP	
DIMENSION : 	129	
EDGE_WEIGHT_TYPE : 	EUC_2D	
CAPACITY : 	39	
NODE_COORD_SECTION		
1	0	0
2	51	570
...
129	610	408
DEMAND_SECTION		
1	0	
2	9	
...
129	3	
DEPOT_SECTION		
	1	
	-1	
EOF		

（2）输出文件的内容解析

每个实例的输出文件包含.out和.sol，以X-n129-k18.vrp的求解结果为例。

其中的，X-n129-k18.vrp.out的内容为计算结果总览，如下所示仅包含两个值，首个值$29066$表示成本（距离）目标值，第二个值$64$表示求解时间，单位是$s$。

29066	64

另一个文件X-n129-k18.vrp_seed-0.vrp.sol的内容为结果方案，记录了全部的结果路线。

Route #1: 15 66 89 103 75 91
Route #2: 3
Route #3: 83 109 51 120 69 97 31 99
...
Route #18: 110 17 101 58 40 55 41 61 85 11 126
Cost 29066.000000

（3）对比129个节点的CVRP示例

上述问题是简单的CVRP问题，由$29$个节点以及小容量车型组成。FILO2通过$64s$将问题目标求解到$29066$，利用有商业版License的Gurobi求解相同的算例，通过vrplib解析X-n129-k18.vrp数据集，以双下标变量进行建模，具体的实现代码如下：

import vrplib
from gurobipy import Env, GRB, Model, quicksum
import os
import time


vrplib.download_instance("X-n129-k18", os.path.dirname(os.path.abspath(__file__)))
instance = vrplib.read_instance("X-n129-k18.vrp")

# 连接远程服务器
env = Env(empty=True)
server_address = "Gurobi服务器地址"
env.setParam(GRB.Param.ComputeServer, server_address)
env.setParam(GRB.Param.ServerPassword, "账号密码")
env.start()
m = Model(env=env)
m.Params.MIPGap = 0.001     # 0.001
m.Params.outputFlag = 1    # 0

# 设置CVRP参数
num_customers = instance["dimension"] - 1   # 客户数量
capacity = instance["capacity"]             # 车辆容量
depot = 0                                   # 起点（仓库）编号
demands = instance["demand"]
distances = instance["edge_weight"]
num_vehicles = int(sum(demands)/capacity) + 5                           # 车辆数量

t1 = time.time()
# 定义决策变量
x = {(i, j): m.addVar(vtype=GRB.BINARY, name=f"x_{i}{j}")
     for i in range(num_customers + 1) for j in range(num_customers + 1)}
u = {i: m.addVar(lb=0, ub=capacity, vtype=GRB.CONTINUOUS, name=f"u_{i}")
     for i in range(num_customers + 1)}    # 车辆到达i的剩余容量（累计需求）

# 目标函数：最小化总距离
m.setObjective(quicksum(distances[i][j] * x[i, j]
                        for i in range(num_customers + 1) for j in range(num_customers + 1) if i != j
                        ), GRB.MINIMIZE)

# 约束条件
# 1. 所有车辆由中心出发，最后回到中心
m.addConstr(quicksum(x[0, j] for j in range(1, num_customers + 1)) ==
            quicksum(x[j, 0] for j in range(1, num_customers + 1)), name=f"flow_balance")

# 2. 车辆动态平衡，每个客户点只被服务1次
for j in range(1, num_customers + 1):
    m.addConstr(quicksum(x[i, j] for i in range(num_customers + 1) if i != j) == 1, name=f"once_visit")
    m.addConstr(quicksum(x[i, j] for i in range(num_customers + 1) if i != j) ==
                    quicksum(x[j, i] for i in range(num_customers + 1) if i != j), name=f"flow_balance")

# 3. 车辆容量约束（MTZ约束，消除子环路）
for i in range(num_customers + 1):
    for j in range(num_customers + 1):
        if i != j and j != depot:
            m.addConstr(u[i] + demands[i] <= u[j] + capacity * (1 - x[i, j]), name=f"capacity_{i}{j}")

t2 = time.time()
print("建模时长", t2 - t1)
# 求解模型
m.optimize()
t3 = time.time()
print("求解时长：", t3 - t2)

# 输出结果
if m.status == GRB.OPTIMAL:
    print("Optimal solution found!")
    for i in range(num_customers + 1):
        for j in range(num_customers + 1):
            if x[i, j].x > 0.5:
                print(f"Edge ({i}, {j}) with distance {distances[i][j]}")
else:
    print("No optimal solution found.")

然而，CVRP问题具有极强的对称性（下界收敛会非常慢）的特点，通过Gurobi求解到$753s$时，可行解的目标值只下降到$32177.2066$，而FILO2通过$64s$将问题目标求解到$29066$。从这个实验结果看，FILO2的在大规模CVRP问题上的表现非常好，能在可接受的时间范围内给出一个不错的解。

当然有人会说，子环路消除约束的复杂度非常大，通过处理成惰性约束可以增加求解效率。但随着问题规模的增大，启发式的求解效率将越来越有优势，尽管无法说明其具体理论最优值的（最小）Gap。

由于小编用的是服务器版的Gurobi，由于通信开销限制了实时回调的使用，因此，如果感兴趣的话读者可以自行在本地测试回调函数，具体可以参考文章：利用Gurobi的Callback方法在求解过程中添加惰性约束。

4. 总结

经过小编的多组其他数据集的实验，FILO2在大规模CVRP问题上的优势非常显著，还可以探索不同的随机种子，以生成更优的结果；对于具有标准数据格式的CVRP问题而言，FILO2的使用方式也极其简单。

总体而言，FILO2作为大规模的CVRP的启发式工具，相比线性规划求解器能在可接受的时间范围内，获得一个较优的解，如果想获取解的准确Gap，可以将FILO2得到的较优解作为规划求解器的初始解。

---

1. FILO slides: https://github.com/acco93/filo/blob/master/docs/slides.pdf ↩︎

2. FILO Research Report: https://github.com/acco93/filo/blob/master/docs/report.pdf ↩︎

3. Luca Accorsi, Daniele Vigo (2021) A Fast and Scalable Heuristic for the Solution of Large-Scale Capacitated Vehicle Routing Problems. Transportation Science 55(4):832-856 ↩︎