VRP问题求解工具
Software Tools for Vehicle Routing
# 求解效率、灵活度各具特色 #
文章前言 |SUMMER
《INFORMS Journal on Computing》的软件工具专刊与第十二届 DIMACS 车辆路径问题实现挑战赛以及后续的车辆路径计算方法研讨会一同发布。这些活动将研究人员聚集在一起，通过计算实验评估解决各种车辆路径问题（VRP）变体的最新技术水平。这一举措以及这两个活动都是向David S. Johnson的遗产致敬，他是计算算法领域的先驱人物，也是 DIMACS 算法实现挑战赛的创始人。
2021 - 2022 年挑战赛考虑了七种经典和新兴的 VRP 变体，包括容量限制车辆路径问题、带时间窗的车辆路径问题、库存路径问题、拆分配送、电动汽车路径问题、弧路径问题和动态网约车问题。来自 36 个参赛团队的 47 份求解器提交成果凸显了开源软件工具在解决这些困难优化问题上实现卓越计算性能的价值。
在这一成功的基础上，后续研讨会强调了车辆路径问题的计算方法和算法，整合了 DIMACS 实现挑战赛、EURO Meets NeurIPS 竞赛以及亚马逊最后一英里挑战赛的见解。研讨会展示了表现出色的团队的报告。
认识到这些活动带来的进步，我们为这一特别版块发布了征稿通知，旨在突出车辆路径研究社区有影响力的开源软件和数据贡献。
四篇被接受的论文展示了最先进的求解器软件包和算法。
Part.01
PyVRP
Wouda 等人（2024 年）论文中的 PyVRP Python 软件包实现了一种混合遗传搜索求解器，专注于带时间窗的车辆路径问题，将 Python 的灵活性与高性能的 C++ 组件相结合。
from pyvrp import Modelfrom pyvrp.plotting import plot_solutionfrom pyvrp.stop import MaxRuntimeCOORDS = [(456, 320),(228, 0),(912, 0),(0, 80),(114, 80),(570, 160),(798, 160),(342, 240),(684, 240),(570, 400),(912, 400),(114, 480),(228, 480),(342, 560),(684, 560),(0, 640),(798, 640),]DEMANDS = [0, 1, 1, 2, 4, 2, 4, 8, 8, 1, 2, 1, 2, 4, 4, 8, 8]m = Model()m.add_vehicle_type(4, capacity=15)depot = m.add_depot(x=COORDS[0][0], y=COORDS[0][1])clients = [m.add_client(x=COORDS[idx][0], y=COORDS[idx][1], delivery=DEMANDS[idx])for idx in range(1, len(COORDS))]locations = [depot] + clientsfor frm in locations:for to in locations:distance = abs(frm.x - to.x) + abs(frm.y - to.y) # Manhattanm.add_edge(frm, to, distance=distance)res = m.solve(stop=MaxRuntime(1), display=True)_, ax = plt.subplots(figsize=(8, 8))plot_solution(res.best, m.data(), ax=ax)
Part.02
RoutingBlocks
Klein 和 Schiffer（2024）的 “RoutingBlocks” 是一个多功能的 Python 工具包，旨在通过优化的数据结构和问题抽象，简化具有中间站点（例如电动汽车路径问题）的车辆路径问题算法的开发。
import randomimport sysfrom itertools import productfrom pathlib import Pathfrom typing import Listimport routingblocksimport routingblocks_cvrp as cvrpimport vrplibdef read_instance(instance_name: str, basedir: Path = Path('instances/')):instance_file = basedir / instance_nameif not instance_file.exists():basedir.mkdir(parents=True, exist_ok=True)# Download the CVRP problem instance if it does not existvrplib.download_instance(instance_name, str(instance_file))# Load the CVRP problem instancereturn vrplib.read_instance(instance_file)def create_cvrp_instance(instance):# Create CVRP verticesn = len(instance['demand'])vertices = [cvrp.create_cvrp_vertex(i, str(i), False, i == instance['depot'][0], cvrp.CVRPVertexData(instance['demand'][i]))for i in range(n)]# Create CVRP arcsarcs = [[cvrp.create_cvrp_arc(cvrp.CVRPArcData(instance['edge_weight'][i][j])) for j in range(n)] for i in range(n)]return routingblocks.Instance(vertices, arcs, len(vertices))def main(instance_name: str):instance = read_instance(instance_name)cpp_instance = create_cvrp_instance(instance)evaluation = cvrp.CVRPEvaluation(instance['capacity'])max_demand = instance['demand'].max()max_dist = instance['edge_weight'].max()# Set the penalty factor for load violationsevaluation.overload_penalty_factor = max_dist / max_demand# Create a simple solution by applying local search to a random solutionrandgen = routingblocks.Random()random_solution = generate_random_solution(evaluation, cpp_instance, randgen)optimized_solution = optimize_solution(evaluation, cpp_instance, random_solution)print(optimized_solution)print("Cost:", optimized_solution.cost, optimized_solution.cost_components)def optimize_solution(evaluation: routingblocks.Evaluation, instance: routingblocks.Instance,solution: routingblocks.Solution):# Create a local search solversolver = routingblocks.LocalSearch(instance, evaluation, evaluation, routingblocks.BestImprovementPivotingRule())# Create some operators# Create the arc set - by default all arcs are includedfull_arc_set = routingblocks.ArcSet(len(instance))operators = [routingblocks.operators.SwapOperator_0_1(instance, full_arc_set),routingblocks.operators.SwapOperator_1_1(instance, full_arc_set),routingblocks.operators.InterRouteTwoOptOperator(instance, full_arc_set)]# Optimize the solution (inplace)solver.optimize(solution, operators)return solutiondef distribute_randomly(sequence, num_subsequences: int, randgen=random.Random()) -> List[List]:subsequences = [[] for _ in range(num_subsequences)]for item in sequence:subsequences[randgen.randint(0, len(subsequences) - 1)].append(item)return subsequencesdef generate_random_solution(evaluation: routingblocks.Evaluation, instance: routingblocks.Instance,random: routingblocks.Random):customers = [x.vertex_id for x in instance.customers]sol = routingblocks.Solution(evaluation, instance,[routingblocks.create_route(evaluation, instance, r) for r indistribute_randomly(customers, instance.fleet_size,random)])return solif __name__ == '__main__':if len(sys.argv) != 2:print(f'Usage: {sys.argv[0]} <instance_name>')sys.exit(1)main(sys.argv[1])
Part.03
VRPSolverEasy
Errami等人（2023）的 “VRPSolverEasy” 提供了一个 Python 库，它与 “VRPSolver”（用于复杂车辆路径问题的先进分支切割定价求解器）相连接。通过这个软件包，用户可以根据熟悉的元素（如仓库、客户、连接和车辆类型）来定义问题，而无需或只需很少的混合整数规划建模知识。
import VRPSolverEasy as vrpseimport mathdef compute_euclidean_distance(x_i, x_j, y_i, y_j):"""compute the euclidean distance between 2 points from graph"""return round(math.sqrt((x_i - x_j)**2 +(y_i - y_j)**2), 3)# Datacost_per_distance = 10begin_time = 0end_time = 5000nb_point = 7# Map with names and coordinatescoordinates = {"Wisconsin, USA": (44.50, -89.50),"West Virginia, USA": (39.000000, -80.500000),"Vermont, USA": (44.000000, -72.699997),"Texas, the USA": (31.000000, -100.000000),"South Dakota, the US": (44.500000, -100.000000),"Rhode Island, the US": (41.742325, -71.742332),"Oregon, the US": (44.000000, -120.500000)}# Demands of pointsdemands = [0, 500, 300, 600, 658, 741, 436]# Initialisationmodel = vrpse.Model()# Add vehicle typemodel.add_vehicle_type(id=1,start_point_id=0,end_point_id=0,name="VEH1",capacity=1100,max_number=6,var_cost_dist=cost_per_distance,tw_end=5000)# Add depotmodel.add_depot(id=0, name="D1", tw_begin=0, tw_end=5000)coordinates_keys = list(coordinates.keys())# Add customersfor i in range(1, nb_point):model.add_customer(id=i,name=coordinates_keys[i],demand=demands[i],tw_begin=begin_time,tw_end=end_time)# Add linkscoordinates_values = list(coordinates.values())for i in range(0, 7):for j in range(i + 1, 7):dist = compute_euclidean_distance(coordinates_values[i][0],coordinates_values[j][0],coordinates_values[i][1],coordinates_values[j][1])model.add_link(start_point_id=i,end_point_id=j,distance=dist,time=dist)model.solve()model.export()if model.solution.is_defined():print(model.solution)
Part.04
AILS-II
Maximo 等人（2024）展示了 AILS-II，这是一种自适应迭代局部搜索元启发式算法的 Java 实现，能够有效地解决大规模容量限制车辆路径问题实例。
package SearchMethod;import java.lang.reflect.InvocationTargetException;import java.text.DecimalFormat;import java.util.HashMap;import java.util.Random;import Auxiliary.Distance;import Data.Instance;import DiversityControl.DistAdjustment;import DiversityControl.OmegaAdjustment;import DiversityControl.AcceptanceCriterion;import DiversityControl.IdealDist;import Improvement.LocalSearch;import Improvement.IntraLocalSearch;import Improvement.FeasibilityPhase;import Perturbation.InsertionHeuristic;import Perturbation.Perturbation;import Solution.Solution;public class AILSII{//----------Problema------------Solution solution,referenceSolution,bestSolution;Instance instance;Distance pairwiseDistance;double bestF=Double.MAX_VALUE;double executionMaximumLimit;double optimal;//----------caculoLimiar------------int numIterUpdate;//----------Metricas------------int iterator,iteratorMF;long first,ini;double timeAF,totalTime,time;Random rand=new Random();HashMap<String,OmegaAdjustment>omegaSetup=new HashMap<String,OmegaAdjustment>();double distanceLS;Perturbation[] pertubOperators;Perturbation selectedPerturbation;FeasibilityPhase feasibilityOperator;ConstructSolution constructSolution;LocalSearch localSearch;InsertionHeuristic insertionHeuristic;IntraLocalSearch intraLocalSearch;AcceptanceCriterion acceptanceCriterion;// ----------Mare------------DistAdjustment distAdjustment;// ---------Print----------boolean print=true;IdealDist idealDist;double epsilon;DecimalFormat deci=new DecimalFormat("0.0000");StoppingCriterionType stoppingCriterionType;public AILSII(Instance instance,InputParameters reader){this.instance=instance;Config config=reader.getConfig();this.optimal=reader.getBest();this.executionMaximumLimit=reader.getTimeLimit();this.epsilon=config.getEpsilon();this.stoppingCriterionType=config.getStoppingCriterionType();this.idealDist=new IdealDist();this.solution =new Solution(instance,config);this.referenceSolution =new Solution(instance,config);this.bestSolution =new Solution(instance,config);this.numIterUpdate=config.getGamma();this.pairwiseDistance=new Distance();this.pertubOperators=new Perturbation[config.getPerturbation().length];this.distAdjustment=new DistAdjustment( idealDist, config, executionMaximumLimit);this.intraLocalSearch=new IntraLocalSearch(instance,config);this.localSearch=new LocalSearch(instance,config,intraLocalSearch);this.feasibilityOperator=new FeasibilityPhase(instance,config,intraLocalSearch);this.constructSolution=new ConstructSolution(instance,config);OmegaAdjustment newOmegaAdjustment;for (int i = 0; i < config.getPerturbation().length; i++){newOmegaAdjustment=new OmegaAdjustment(config.getPerturbation()[i], config,instance.getSize(),idealDist);omegaSetup.put(config.getPerturbation()[i]+"", newOmegaAdjustment);}this.acceptanceCriterion=new AcceptanceCriterion(instance,config,executionMaximumLimit);try{for (int i = 0; i < pertubOperators.length; i++){this.pertubOperators[i]=(Perturbation) Class.forName("Perturbation."+config.getPerturbation()[i]).getConstructor(Instance.class,Config.class,HashMap.class,IntraLocalSearch.class).newInstance(instance,config,omegaSetup,intraLocalSearch);}} catch (InstantiationException | IllegalAccessException | IllegalArgumentException| InvocationTargetException | NoSuchMethodException | SecurityException| ClassNotFoundException e) {e.printStackTrace();}}public void search(){iterator=0;first=System.currentTimeMillis();referenceSolution.numRoutes=instance.getMinNumberRoutes();constructSolution.construct(referenceSolution);feasibilityOperator.makeFeasible(referenceSolution);localSearch.localSearch(referenceSolution,true);bestSolution.clone(referenceSolution);while(!stoppingCriterion()){iterator++;solution.clone(referenceSolution);selectedPerturbation=pertubOperators[rand.nextInt(pertubOperators.length)];selectedPerturbation.applyPerturbation(solution);feasibilityOperator.makeFeasible(solution);localSearch.localSearch(solution,true);distanceLS=pairwiseDistance.pairwiseSolutionDistance(solution,referenceSolution);evaluateSolution();distAdjustment.distAdjustment();selectedPerturbation.getChosenOmega().setDistance(distanceLS);//updateif(acceptanceCriterion.acceptSolution(solution))referenceSolution.clone(solution);}totalTime=(double)(System.currentTimeMillis()-first)/1000;}public void evaluateSolution(){if((solution.f-bestF)<-epsilon){bestF=solution.f;bestSolution.clone(solution);iteratorMF=iterator;timeAF=(double)(System.currentTimeMillis()-first)/1000;if(print){System.out.println("solution quality: "+bestF+" gap: "+deci.format(getGap())+"%"+" K: "+solution.numRoutes+" iteration: "+iterator+" eta: "+deci.format(acceptanceCriterion.getEta())+" omega: "+deci.format(selectedPerturbation.omega)+" time: "+timeAF);}}}private boolean stoppingCriterion(){switch(stoppingCriterionType){case Iteration: if(bestF<=optimal||executionMaximumLimit<=iterator)return true;break;case Time: if(bestF<=optimal||executionMaximumLimit<(System.currentTimeMillis()-first)/1000)return true;break;}return false;}public static void main(String[] args){InputParameters reader=new InputParameters();reader.readingInput(args);Instance instance=new Instance(reader);AILSII ailsII=new AILSII(instance,reader);ailsII.search();}public Solution getBestSolution() {return bestSolution;}public double getBestF() {return bestF;}public double getGap(){return 100*((bestF-optimal)/optimal);}public boolean isPrint() {return print;}public void setPrint(boolean print) {this.print = print;}public Solution getSolution() {return solution;}public int getIterator() {return iterator;}public String printOmegas(){String str="";for (int i = 0; i < pertubOperators.length; i++){str+="\n"+omegaSetup.get(this.pertubOperators[i].perturbationType+""+referenceSolution.numRoutes);}return str;}public Perturbation[] getPertubOperators() {return pertubOperators;}public double getTotalTime() {return totalTime;}public double getTimePerIteration(){return totalTime/iterator;}public double getTimeAF() {return timeAF;}public int getIteratorMF() {return iteratorMF;}public double getConvergenceIteration(){return (double)iteratorMF/iterator;}public double convergenceTime(){return (double)timeAF/totalTime;}}
参考文献 | References
Errami N, Queiroga E, Sadykov R, Uchoa E (2024) Vrpsolvereasy: A python library for the exact solution of a rich vehicle routing problem. INFORMS J. Comput. 36(4):956–964.
Klein PS, Schiffer M (2024) Routingblocks: An open-source python package for vehicle routing problems with intermediate stops. INFORMS J. Comput. 36(4):966–973.
Maximo V, Cordeau JF, Nascimento M (2024) AILS-II: An adaptive iterated local search heuristic for the large-scale capacitated vehicle routing problem. INFORMS J. Comput. 36(4):974–986.
Wouda NA, Lan L, Kool W (2024) PyVRP: A high-performance VRP solver package. INFORMS J. Comput. 36(4):943ߝ955.
