运筹系列93：VRP精确算法

1. 基础版本

定义$x_{ijk}$为边$ij$是否由车辆$k$去运输。如果有时间窗约束的话，再加上一个变量$c_{ik}$即可，表示第k辆车到达节点i时的时间点。
第一类客户流量约束，要求每个点都有1个入度和1个出度，并且需要是同一辆车；
第二类车辆容量约束，要求每辆车的总访问边数为容量+1，并且仅有一个从depot出去的边

using JuMP, HiGHS

K = 3
N = 13 #34
Q = 4 #13
q = ones(Int,N);q[1]=0
CVRP = Model(HiGHS.Optimizer)
set_silent(CVRP)
@variable(CVRP,x[1:N,1:N,1:K],Bin)
for i in 2:N
    @constraint(CVRP, sum(x[i,:,:]) == 1)
    @constraint(CVRP, sum(x[:,i,:]) == 1)
    for k in 1:K;@constraint(CVRP, sum(x[i,:,k])==sum(x[:,i,k]));end
end
for k in 1:K
    @constraint(CVRP, sum(x[1,2:N,k]) == 1)
    @constraint(CVRP, sum(x[:,:,k])==Q+1)
    for i in 1:N, j in 1:N;@constraint(CVRP, x[i,j,k]+x[j,i,k]<=1);end
end
@objective(CVRP,Min, sum(x[i,j,k]*distmat[i,j] for i=1:N,j=1:N,k in 1:K))
@time optimize!(CVRP)
draw_vrp(sum(value.(x),dims = 3))

速度比后面的MTZ快一些，但是整体还是慢。

2. MTZ模型

MTZ是Miller-Tucker-Zemlin inequalities的缩写。除了定义是否用到边$x_{ij}$外，还需要定义一个$u_i$用来表示此时车辆的当前载货量。注意这里x变量需要定义为有向。
这里定义为pickup问题，代码为：

using JuMP, HiGHS

k = 3 # number of vehicles
N = 11 # number of points, 0 as depot
Q = 4 # vehicle capacity
q = ones(Int,N);q[1]=0 # demand
CVRP = Model(HiGHS.Optimizer)
set_silent(CVRP)
@variable(CVRP,x[1:N,1:N],Bin)
@variable(CVRP,u[1:N],lower_bound = 0, upper_bound = Q)
# 约束1：出度约束
for i in 2:N
    @constraint(CVRP, sum(x[i,1:i-1]) + sum(x[i,i+1:N]) == 1)
    @constraint(CVRP, sum(x[1:i-1,i]) + sum(x[i+1:N,i]) == 1)
end
@constraint(CVRP, sum(x[1,1:N]) == k)
# 约束2：流量约束。若存在i->j，则u_j-u_i==q_j；否则u_j-q_j和u_i没有关系。此外，需要有u_j-q_j>=0
for i=2:N, j=[2:i-1;i+1:N]
    @constraint(CVRP,u[i] - u[j] + Q*x[i,j] <= Q-q[j])
end
for i in 2:N
    @constraint(CVRP,q[i] <= u[i] <= Q)
end

@objective(CVRP,Min, sum(x[i,j]*distmat[i,j] for i=1:N,j=1:N))
@time optimize!(CVRP)

MTZ的求解速度不快，10个点3辆车都需要3秒左右时间。

3. 分支定界法（行生成）

使用Two-index vehicle flow formulations。按照tsp的方式使用行生成法速度极慢（cut的效率太低），因此考虑使用branch-and-cut直接求解。需要cut的主要有2个：1）容量约束；2）subtour约束。如下例子：

using TSPLIB,JuMP, HiGHS, Distances
N = 13
Q = 4
k = 3
m = Model(HiGHS.Optimizer)
set_silent(m)
@variable(m, x[1:N,1:N]>=0,Bin)
@objective(m, Min, sum(x[i,j]*distmat[i,j] for i=1:N,j=1:N))
@constraint(m, x[1,1] == 0)
@constraint(m, sum(x[1,j] for j in 2:N) == k)
@constraint(m, sum(x[j,1] for j in 2:N) == k)
for i=2:N 
    for j in 1:N;@constraint(m, x[i,j]+x[j,i] <= 1);end
    @constraint(m, sum(x[i,j] for j in 1:N) == 1)
    @constraint(m, sum(x[j,i] for j in 1:N) == 1)
end
optimize!(m)
draw_vrp(x)

接下来定义寻找tour的函数，以及branch and cut的代码：

# find all subtours
function tours(x)
    g = JuMP.value.(x)
    # 第一步，找到所有从1出发的tour
    abnormal_paths = []
    paths = []
    path = [1]
    left = collect(1:N)
    while true
        v, idx = findmax(g[path[end],left])
        if v==0
            break
        else
            g[left[idx],path[end]]=0
            g[path[end],left[idx]]=0
            push!(path,left[idx])
        end
        if path[end]==1
            if length(path)>Q+2;push!(abnormal_paths,path);end
            push!(paths,path)
            path = [1]
            setdiff!(left,path[2:end-1])
        end
    end
    # 第二步，找到所有孤立的环(subtour)
    left = collect(1:N)
    for path in paths;setdiff!(left,path);end
    while length(left)>0   
        path = [left[1]]
        while true
            v, idx = findmax(g[path[end],left])
            if idx == 1
                break
            else
                g[left[idx],path[end]]=0
                g[path[end],left[idx]]=0
                push!(path,left[idx])
            end
        end
        setdiff!(left,path)
        push!(paths,path)
        push!(abnormal_paths,path)
    end
    return paths,abnormal_paths
end
    
paths,abnormal_paths = tours(x)
while length(abnormal_paths) > 0
    for path in paths
        s = setdiff(path,1)
        sn = setdiff(2:N,s)
        @constraint(m, sum(x[i,j] for i in s, j in setdiff(1:N,s)) >= ceil(length(s)/Q))
        @constraint(m, sum(x[i,j] for i in sn, j in setdiff(1:N,sn)) >= ceil(length(sn)/Q))
    end
    optimize!(m)
    paths,abnormal_paths = tours(x)
end
draw_vrp(x)

4. set-partitioning方法

方法很直观，把所有的子路径用TSP问题求解（使用Concorde库），然后用set-partitioning方法选择最合适的几条路线组合成VRP的结果。

using JuMP, HiGHS, Combinatorics, Concorde 

k = 3
N = 13 #34
Q = 4 #13

function getRoutes(k,N,Q)
    Qm = N-1-(k-1)*Q
    route_dists = Dict()
    # 求解所有子路径的最优解
    for q in Qm:Q
        for c in combinations(2:N,q)
            c_index_tour,c_tour_length = Concorde.solve_tsp(floor.(Int,distmat[[1;c],[1;c]].*100)) 
            c_tour = [1;c][c_index_tour]
            route_dists[c_tour] = c_tour_length
        end
    end
    route_dists
end
@time route_dists = getRoutes(k,N,Q);
CVRP = Model(HiGHS.Optimizer)
set_silent(CVRP)
routes = collect(keys(route_dists))
route_dists = collect(values(route_dists))
rn = length(routes)
@variable(CVRP,x[1:rn],Bin)
@objective(CVRP,Min, sum(x[i]*route_dists[i] for i in 1:rn))
a = zeros(Int,rn,N)
for i in 1:rn,j in routes[i];a[i,j]=1;end
for j in 2:N;@constraint(CVRP,sum(x[i]*a[i,j] for i in 1:rn)==1);end
@constraint(CVRP,sum(x)==k)
@time optimize!(CVRP)
rs = routes[findall(x->x>0.1,value.(x))]
plt.figure(figsize=(10,7)) 
plt.scatter(pos[1:N,1],pos[1:N,2],c="red")
for i in 1:N;plt.text(pos[i,1], pos[i,2], i);end
for r in rs
    l = pos[[r;1],:]
    PyPlot.plot(l[:,1], l[:,2], color="b")
end

5. 列生成方法

接下来设计基于set-partitioning的列生成算法。首先基于启发式算法给出一些tsp子路径。这里模型变量注意要改成线性问题，不然无法给出对偶解。

# 用启发式算法先生成一些路径
routes = [[2,3,4,5],[4,5,6,7],[6,7,8,9],[2,3,4,6]]
route_dists = [Concorde.solve_tsp(floor.(Int,distmat[[1;r],[1;r]].*100))[2]/100 for r in routes]
rn = length(routes)
cg = Model(HiGHS.Optimizer)
set_silent(cg)
@variable(cg,1>=x[1:rn]>=0)#,Bin)
@objective(cg,Min, sum(x[i]*route_dists[i] for i in 1:rn))
a = zeros(Int,rn,N)
for i in 1:rn,j in routes[i];a[i,j]=1;end
cons = []
for j in 2:N;con = @constraint(cg,sum(x[i]*a[i,j] for i in 1:rn)==1);push!(cons,con);end
@constraint(cg,sum(x)==k)
@time optimize!(cg)

接下来我们建立子问题来生成新的路径，针对上面的主问题，每一个路径都对应一个残差：cx-πx（这里π是对偶变量的解），我们以最小化残差为目标，如果结果小于0，则说明可以将新的路径转入基变量，完整代码如下：

routes = [[2,3,4,5],[4,5,6,7],[6,7,8,9],[2,3,4,6]]
flag = true
while flag
    route_dists = [Concorde.solve_tsp(floor.(Int,distmat[[1;r],[1;r]].*100))[2]/100 for r in routes]
    rn = length(routes)
    cg = Model(HiGHS.Optimizer)
    set_silent(cg)
    @variable(cg,1>=x[1:rn]>=0)#,Bin)
    @objective(cg,Min, sum(x[i]*route_dists[i] for i in 1:rn))
    a = zeros(Int,rn,N)
    for i in 1:rn,j in routes[i];a[i,j]=1;end
    cons = []
    for j in 2:N;con = @constraint(cg,sum(x[i]*a[i,j] for i in 1:rn)==1);push!(cons,con);end
    optimize!(cg)
    sp = Model(HiGHS.Optimizer)
    set_silent(sp)
    @variable(sp,1>=y[1:N,1:N]>=0,Bin)
    p = [0.0]
    for c in cons;push!(p,dual(c));end
    @objective(sp, Min, sum(distmat[i,j]*y[i,j] for i in 1:N, j in 1:N)-sum(p[i]*sum(y[i,j] for j in 1:N) for i in 2:N))
    for i in 1:N
        for j in 1:N;@constraint(sp,y[i,j]+y[j,i]<=1);end
        @constraint(sp,sum(y[i,j] for j in 1:N)==sum(y[j,i] for j in 1:N))
    end;
    @constraint(sp,sum(y)==Q+1)
    @constraint(sp,sum(y[1,:])==1)
    @constraint(sp,sum(y[:,1])==1)
    optimize!(sp)
    flag = (objective_value(sp)<-0.01)
    yv = value.(y)
    route = [1]
    while length(route)<=Q
        push!(route,findfirst(x->x>=0.5,yv[route[end],:]))
    end
    route = sort(route[2:end])
    if objective_value(sp)<0
        println("残差为",objective_value(sp),", 转入",route)
    else
        println("残差为",objective_value(sp),", 结束")
    end
    push!(routes,route)
end

结果如下：

残差为-38.831449600051535, 转入[2, 5, 7, 9]
残差为-39.150033402407914, 转入[3, 5, 6, 9]
残差为-38.81684590554077, 转入[5, 6, 7, 8]
残差为-49.543219451392396, 转入[3, 4, 7, 8]
残差为-24.007403774319204, 转入[4, 5, 8, 9]
残差为-32.952086371511896, 转入[2, 3, 8, 9]
残差为-16.53067003901679, 转入[2, 3, 7, 9]
残差为-13.61155193986133, 转入[2, 5, 8, 9]
残差为-10.133530272379385, 转入[3, 4, 5, 6]
残差为-8.572327787226246, 转入[2, 5, 6, 9]
残差为-7.603418916493028, 转入[5, 6, 7, 9]
残差为-7.600187922969602, 转入[5, 6, 8, 9]
残差为-5.700835963904262, 转入[2, 3, 4, 8]
残差为-13.765527089229668, 转入[2, 3, 4, 9]
残差为-1.4118420231930586, 转入[3, 4, 5, 9]
残差为-7.4652916644221925, 转入[5, 7, 8, 9]
残差为-7.5561902615201175, 转入[2, 3, 4, 7]
残差为-4.869792462350416, 转入[3, 4, 7, 9]
残差为-1.337641019305987, 转入[3, 5, 7, 8]
残差为-7.712871227881135, 转入[2, 4, 7, 8]
残差为-4.741848715534764, 转入[2, 4, 5, 6]
残差为-5.479199809684055, 转入[2, 3, 5, 6]
残差为-2.4061894903914096, 转入[2, 7, 8, 9]
残差为-1.413493799681662, 转入[2, 4, 5, 9]
残差为-4.5133099650613415, 转入[4, 7, 8, 9]
残差为-1.2752916644218404, 转入[5, 7, 8, 9]
残差为-1.0853635570299593, 转入[2, 3, 4, 6]
残差为-0.5040970280641338, 转入[2, 3, 7, 8]
残差为0.023151284460664213, 结束

此时主问题的解即为最优解。（如果时间有限，我们也可以在中间停下，也是可行解。）

5. 关于测试数据

测试案例可参考 http://vrp.atd-lab.inf.puc-rio.br/index.php/en/。
我们这里用的数据为：

pos = [121.472641	31.231707
123.429092	41.796768
125.324501	43.886841
126.642464	45.756966
116.405289	39.904987
117.190186	39.125595
111.75199	40.84149
106.23248	38.48644
112.549248	37.857014
114.502464	38.045475
117.000923	36.675808
113.665413	34.757977
108.948021	34.263161
114.298569	30.584354
118.76741	32.041546
117.283043	31.861191
112.982277	28.19409
115.892151	28.676493
120.15358	30.287458
119.306236	26.075302
113.28064	23.125177
121.520076	25.030724
110.19989	20.04422
108.320007	22.82402
106.504959	29.533155
102.71225	25.040609
106.713478	26.578342
104.065735	30.659462
103.83417	36.06138
101.77782	36.61729
91.1145	29.64415
87.61688	43.82663
114.16546	22.27534
113.54913	22.19875];

function dis(i,j)
    A = pos[i,:];B = pos[j,:]
    sqrt(sum((A-B).^2))
end

function drawTree(t,n)
    plt.figure(figsize=(10,7)) 
    plt.scatter(pos[1:n,1],pos[1:n,2],c="red")
    for i in 1:n;plt.text(pos[i,1], pos[i,2], i);end
    for i in 1:length(t)
        l = pos[collect(t[i]),:]
        PyPlot.plot(l[:,1], l[:,2], color="b")
    end
end

function draw_vrp(x)
    xv = value.(x)
    t = []
    for i in 1:size(x)[1],j in 1:size(x)[1]
        if xv[i,j]>0.1;push!(t,(i,j));end
    end
    drawTree(t,size(x)[1])
end