vrp 节约算法 c++_列生成算法求解车辆路径问题的简单实现

准备

车辆路径问题的变种有很多种，这次我们考虑较为简单的CVRP (capacitated vehicle routing problem)。可参考

SysOpt：(FOCUS)有容量限制的 VRP 模型的分支切割法的实现 by python3 & Gurobi​zhuanlan.zhihu.com

列生成法的简单框架

• Step A 准备能够得到可行解的变量
• Step B 在Step A中准备的变量，用主问题求解 (线性规划问题) 得到最优解
• Step C 使用最优解的对偶信息定义子问题，求解
• Step C-1 如果得到的最优值是负的，从最优解中生成新变量(列)加入到主问题中，Go to Step B
• Step C-2 如果得到的最优值是非负，中止

本次使用的主问题和子问题可参考如下

俽澔：基于基本最短路径列生成的车辆路径问题​zhuanlan.zhihu.com

注: 上述链接的问题带有时间窗，大家可以自动忽略关于时间的约束

数据

本次使用的数据可参考

master苏：ortools系列：容量约束VRP问题​zhuanlan.zhihu.com

计算环境

我使用的计算环境如下

• OS: Microsoft Windows 10 Pro
• CPU: Intel(R) Core(TM) i7-6950X
• Memory: 64.0 GB
• Solver: Gurobi Optimizer version 8.0.1

结果及改良点

针对上述问题，使用我的服务器大约需要 98s 左右。如果你的machine性能不好时间可能更长。即便如此，想要高速化我个人感觉有两点可再下下功夫

• 使用好的初始解，这次用的是随机生成比如，可以用启发式算法如

SysOpt：对于车辆路径问题的节约里程法的简单实现​zhuanlan.zhihu.com

• 子问题求解的高速化，这次使用的子问题使用模型求解，可以用动态规划等，可参考

俽澔：基于基本最短路径列生成的车辆路径问题​zhuanlan.zhihu.com

Code

####################################
#Content: CG algorithm for solving CVRP
#Author: SysOpt
#Date: 2020.03.12
#Code reference: Mr. Samarthmistry and Prof. Kubo
####################################
from collections import defaultdict
from numpy import *
from gurobipy import *
from itertools import chain, combinations
import time
import numpy as np
import networkx
distance = np.array([
    [0, 548, 776, 696, 582, 274, 502, 194, 308, 194, 536, 502, 388, 354, 468, 776, 662],
    [548, 0, 684, 308, 194, 502, 730, 354, 696, 742, 1084, 594, 480, 674, 1016, 868, 1210],
    [776, 684, 0, 992, 878, 502, 274, 810, 468, 742, 400, 1278, 1164, 1130, 788, 1552, 754],
    [696, 308, 992, 0, 114, 650, 878, 502, 844, 890, 1232, 514, 628, 822, 1164, 560, 1358],
    [582, 194, 878, 114, 0, 536, 764, 388, 730, 776, 1118, 400, 514, 708, 1050, 674, 1244],
    [274, 502, 502, 650, 536, 0, 228, 308, 194, 240, 582, 776, 662, 628, 514, 1050, 708],
    [502, 730, 274, 878, 764, 228, 0, 536, 194, 468, 354, 1004, 890, 856, 514, 1278, 480],
    [194, 354, 810, 502, 388, 308, 536, 0, 342, 388, 730, 468, 354, 320, 662, 742, 856],
    [308, 696, 468, 844, 730, 194, 194, 342, 0, 274, 388, 810, 696, 662, 320, 1084, 514],
    [194, 742, 742, 890, 776, 240, 468, 388, 274, 0, 342, 536, 422, 388, 274, 810, 468],
    [536, 1084, 400, 1232, 1118, 582, 354, 730, 388, 342, 0, 878, 764, 730, 388, 1152, 354],
    [502, 594, 1278, 514, 400, 776, 1004, 468, 810, 536, 878, 0, 114, 308, 650, 274, 844],
    [388, 480, 1164, 628, 514, 662, 890, 354, 696, 422, 764, 114, 0, 194, 536, 388, 730],
    [354, 674, 1130, 822, 708, 628, 856, 320, 662, 388, 730, 308, 194, 0, 342, 422, 536],
    [468, 1016, 788, 1164, 1050, 514, 514, 662, 320, 274, 388, 650, 536, 342, 0, 764, 194],
    [776, 868, 1552, 560, 674, 1050, 1278, 742, 1084, 810, 1152, 274, 388, 422, 764, 0, 798],
    [662, 1210, 754, 1358, 1244, 708, 480, 856, 514, 468, 354, 844, 730, 536, 194, 798, 0]
])

demand = np.array([0, 1, 1, 2, 4, 2, 4, 8, 8, 1, 2, 1, 2, 4, 4, 8, 8])
capacity = 13

numofvehicles=6

numnodes=17

def generate_random_solution(numnodes, distance, capacity, demand, numofvehicles):
    customers = list(range(1,numnodes))
    routes = defaultdict(list)
    
    remaining_capacity = ones(numofvehicles, dtype=int) * capacity
       
    for vehicle in range(numofvehicles):
      
        # Try to feasibly add customers to the vehicle
        for id in customers:
            q = demand[id]
            # If there is remaining capacity, or it is the last vehicle
            if q <= remaining_capacity[vehicle]: #or vehicle == (numofvehicles - 1):
                routes[vehicle].append(id)
                remaining_capacity[vehicle] -= q

        # Remove from the list the customers actually added
        for id in routes[vehicle]:
            customers.remove(id)
                  
        # Add the depot to the start and end of of the route
        routes[vehicle].insert(0, 0)
        routes[vehicle].append(0)
                           
    return routes
def get_routecst_ksit(route,numnodes):
    
    ksit  = [0] * (numnodes+1)
    routecst = 0
    routelength = len(route)
    lastcustomer = 0
    for i,cust in enumerate(route[1:]):
      ksit[cust] +=1      
      routecst += distance[cust][lastcustomer]
      lastcustomer = cust
          
    return routecst,ksit


def callback(model, where):
    if where == GRB.callback.MIPSOL:
        solution = np.array(model.cbGetSolution(model._vars)).reshape(numnodes+1, numnodes+1)
        solution = np.where(solution>0.5, 1, 0) 
        G_sol = networkx.DiGraph(solution)
        
        for route in networkx.simple_cycles(G_sol):
            route.append(route[0])
            subtour = route
            expr = LinExpr()
            for i in range(len(subtour)-1):
                expr += model._vars[(numnodes+1)*subtour[i] + subtour[i+1]]
            model.cbLazy(expr <= len(subtour)-2)

def transform(demand, distance,numnodes):
    dem=np.append(demand,demand[0])
    dist = np.pad(distance, ((0,1), (0,1)), mode='constant', constant_values=9999999)
    for i in range(1, numnodes):
        dist[i,numnodes] = distance[i,0]
        dist[i,0] = 9999999
    dist[numnodes,numnodes] = 0
    return dem, dist

routes=generate_random_solution(numnodes, distance, capacity, demand, numofvehicles)


a = defaultdict(list)
cost = defaultdict(float)
for r in range(len(routes)):
    cost[r],a[r]=get_routecst_ksit(routes[r],numnodes)
initial=len(routes)

dem, dist = transform(demand, distance,numnodes)
# Model
master = Model("MASTER")

# Create decision variables 
y = {}
for r in range(initial):
    y[r] = master.addVar(vtype=GRB.BINARY,name='y_%s' % (r))

master.update()
# constraints
constr = {}
for i in range(numnodes):
    constr[i] = master.addConstr(
      quicksum(a[r][i] * y[r] for r in range(initial)) >= 1)
               
vehicle = master.addConstr(
    quicksum(y[r] for r in range(initial)) <=  numofvehicles)

master.setObjective(quicksum(y[r]*cost[r] for r in range(initial)), GRB.MINIMIZE)

master.update()

EPS = 1.e-6
start_time=time.time()
print("--- CG is start! Good luck! ----")
while True:
    relax=master.relax()                
    relax.optimize()
    pi = [c.Pi for c in relax.getConstrs()]

    subprob=Model("sub")
    x = subprob.addVars(numnodes+1, numnodes+1, vtype=GRB.BINARY, name='x')
    subprob.update()

    subprob.addConstr(quicksum(dem[i]*x.sum(i,'*') for i in range(numnodes+1)) <= capacity, name='capacity')
    
    subprob.addConstr((quicksum(x[0,j] for j in range(1, numnodes+1)) == 1), name='out')

    subprob.addConstrs((x.sum('*',k) == x.sum(k,'*') for k in range(1, numnodes)), name='flow')

    subprob.addConstr((quicksum(x[i,numnodes] for i in range(numnodes)) == 1), name='in')
    
    subprob.addConstrs((x[i,i] == 0 for i in range(numnodes+1)), name='self')

    subprob.addConstrs((x[i,j] + x[j,i]<= 1 for i in range(numnodes) for j in range(numnodes)), name='cycle')
   
    subprob.setObjective(quicksum((dist[i,j]-pi[i])*x[i,j] for i in range(numnodes+1) for j in range(numnodes+1)), GRB.MINIMIZE)
    
    subprob.params.DualReductions = 0
    subprob.params.LazyConstraints = 1
    subprob.params.Threads = 8
    subprob.params.MipFocus=1
    subprob.params.OutputFlag=0
    subprob.update()
    subprob._vars=subprob.getVars()
    subprob.optimize(callback)
    
    if subprob.ObjVal >= 0:
        print('---The CG algorithm is over!---n')
        print("--- The computational time is :%s s ---" % (time.time() - start_time))
        break
    
    solution = np.array(subprob.X).reshape(numnodes+1, numnodes+1)
    solution = np.where(solution>0.5, 1, 0)
    G_sol = networkx.DiGraph(solution)
    routesimple = list(networkx.shortest_path(G_sol, source=0,target=numnodes))
    routesimple.pop()
    routesimple.append(0)
    
    count = len(routes)
    K=count
    
    flag_new_route=True
    for route in routes.values():
        if routesimple == route:
            flag_new_route=False
            break

    if flag_new_route:
        print("---New route---")
        print(routesimple)
    
        routes[K]=routesimple
        cost[K],a[K]=get_routecst_ksit(routes[K],numnodes)
        col=Column()
        for i in range(numnodes):
            col.addTerms(a[K][i],constr[i])
        y[K] = master.addVar(vtype=GRB.BINARY,obj=cost[K], column=col,name='y_%s' % (K))
        master.update()
        master.params.OutputFlag=0

        #master.write("MP" + str(iter) + ".lp")
master.optimize()
for var in master.getVars():
    if var.X > 0.5:
        print('%g %s' % (var.X, var.varName))

大家可以拓展为带有时间窗的车辆路径问题，期待！另外大家可以对比割平面方法，对于小规模问题割平面方法应该优于列生成，但是也有cut一直加，停不下来的情况，具体可参考

SysOpt：(FOCUS)有容量限制的 VRP 模型的分支切割法的实现 by python3 & Gurobi​zhuanlan.zhihu.com

如果有什么疑问或者建议，可以给我发邮件。➡ zhaoyou728 atoutlook.com