# Recitation 10: Column generation and vehicle routing

Today we will spend some time thinking about and implementing column generation, in the particular context of vehicle routing problems. Given time, we will also think about heuristic ways to solve vehicle routing problems.


```julia
using Plots, Random, LinearAlgebra, Printf
```

## 1 - Vehicle routing problems

In an earlier recitation, we talked about the traveling salesman problem (TSP), where the goal is to find a minimum-length tour visiting a set of locations. The TSP has many "cousins", i.e., similar problems with additional complexity to represent practical considerations. One family of such problems are called vehicle routing problems (VRP), and seek to find the optimal set of routes for a *fleet* of vehicles to service a set of customers.

Some examples include:
- Open vehicle routing problem (OVRP): we are given a set of customers that must each be visited exactly once by some vehicle. All routes start from a customer location and must end at a single location, called the "depot". An example of an application of this setting is car-pooling.
- Capacitated vehicle routing problem (CVRP): in this setting, routes must begin and end at the depot, and each customer $i$ is now associated with a demand $d_i$, such that the total demand along any route cannot excced the vehicle capacity $Q$. This is a simplified version of the problem many delivery services face.

## 2 - Solving OVRP using column generation

Let's generate a problem instance, with $N$ locations uniformly sampled from a square (side length $L$)


```julia
using JuMP, Gurobi, LightGraphs

const GRB_ENV = Gurobi.Env()
```


```julia
N = 10
L = 100
```

We'll use the Euclidean metric for distances.


```julia
Random.seed!(1234)
locations = rand(N, 2) .* L
locations = vcat(locations, [L/2 L/2])
dist = [norm(locations[i, :] .- locations[j, :]) for i =1:N+1, j = 1:N+1];
```

Let's plot this example to get a sense of what it looks like!


```julia
plot(size=(500,500), xlim=(0, L), ylim=(0, L))
scatter!(locations[1:N, 1], locations[1:N, 2], label="Customers", legend=:topleft)
scatter!(locations[[N+1], 1], locations[[N+1], 2], label="Depot")
```

Now we need to formulate the problem. 

Setup:
- Let $\mathcal{R}$ designate the set of ALL possible routes starting at some customer and ending at the depot.
- Let $c_r$ designate the length of route $r\in\mathcal{R}$.
- Let $\mathcal{I}_r\subseteq[N]$ designate the set of customers visited by route $r\in\mathcal{R}$.
- For any customer $i\in[N]$, we let $\mathcal{R}_i\subseteq\mathcal{R}$ designate the subset of routes that visit customer $i$.

Then we can formulate OVRP as the following set covering problem:

$$\begin{align}
\min\quad &\sum_{r\in\mathcal{R}}c_r y_r\\
\text{s.t.}\quad & \sum_{r\in\mathcal{R}_i}y_r \ge 1 & \forall i \in [N]\\
&y_r\in\{0,1\} & \forall t\in\mathcal{R}
\end{align}$$

Note that we allow customers to be visited multiple times here. However, given a feasible solution that exhibits this behavior, we can always delete the extra visit and lower the total cost, so we won't worry too much about this.

This formulation is so simple! But of course, the complexity is "hidden" in the exponential size of the set $\mathcal{R}$. In fact, we cannot hope to formulate this problem explicitly as written.

However, most routes in the set $\mathcal{R}$ are useless! For example, consider a random route:


```julia
route = vcat(shuffle(1:N)[1:4], N+1)
```

If we plot this, we see that no one in their right mind would recommend this route!


```julia
plot(size=(500,500), xlim=(0, L), ylim=(0, L))
scatter!(locations[1:N, 1], locations[1:N, 2], label="", legend=:topleft)
scatter!(locations[[N+1], 1], locations[[N+1], 2], label="")
plot!(locations[route, 1], locations[route, 2], label="")
```

So clearly, we can get a good solution (maybe even an optimal one!) to the problem with a restricted set of routes.

Let's consider the linear relaxation of the problem, denoted by (P):

$$\begin{align}
\min\quad &\sum_{r\in\mathcal{R}}c_r y_r\\
\text{s.t.}\quad & \sum_{r\in\mathcal{R}_i}y_r \ge 1 & \forall i \in [N]\\
&y_r\ge 0 & \forall r\in\mathcal{R}
\end{align}$$

Notice that we've eliminated the upper bounds on $y_r$. It doesn't break the validity of the formulation, because from any solution with $y_r>1$, we can construct a feasible solution with $y_r'=1$ which lowers the total cost by $c_r(y_r-1)$.

Now we can write down the dual, denoted by (D):
$$\begin{align}
\max\quad & \sum_{i=1}^N p_i\\
\text{s.t.}\quad & \sum_{i\in\mathcal{I}_r}p_i \le c_r & \forall r\in\mathcal{R}\\
& p_i \ge 0&\forall i\in [N]
\end{align}$$

Assume we solve the main problem with a restricted set of routes  $\mathcal{R}^k\subset\mathcal{R}$, and we get a primal feasible solution $y^k$, and a dual solution $p^k$. The solution $y^k$ is optimal in the full primal if the corresponding dual solution $p^k$ is feasible in the full dual.

For dual feasibility, we require that $\sum_{i\in\mathcal{I}_r}p^k_i\le c_r$ for **every** route $r\in\mathcal{R}$. We can re-write this requirement as $c_r-\sum_{i\in\mathcal{I}_r}p^k_i\ge 0$. In other words, the **reduced cost** of each variable $y_r$ must be nonnegative.

So we propose the following column generation algorithm:
1. Solve (P) over a restricted set of routes $\mathcal{R}^k$, obtain a primal solution $y^k$ and a dual solution $p^k$
2. Find the route $r$ with the smallest reduced cost $\bar{c} = c_r-\sum_{i\in\mathcal{I}_r}p^k_i$
3. If $\bar{c}\ge 0$, terminate. Otherwise add the corresponding $r$ to $\mathcal{R}^k$ (obtaining $\mathcal{R}^{k+1}$) and iterate.

Let's start implementing this algorithm! First, we need an initial set of routes that will produce a feasible primal solution. An easy place to start is to just have each customer in its own route.


```julia
"Create initial set of routes and associated costs"
function initial_routes(distances::Matrix)
    n = size(distances, 1) - 1 # number of customers
    return [[i, n+1] for i in 1:n], distances[1:n, n+1]
end
```

Let's plot these initial routes!


```julia
plot(size=(500,500), xlim=(0, L), ylim=(0, L))
scatter!(locations[1:N, 1], locations[1:N, 2], label="", color="blue", legend=:topleft)
routes, costs = initial_routes(dist)
for route in routes
    plot!(locations[route, 1], locations[route, 2], label="", color="black")
end
scatter!(locations[[N+1], 1], locations[[N+1], 2], label="", color="orange")
```

Second, we need a way to compute the smallest reduced cost, i.e.
$$(\min_{r\in\mathcal{R}}c_r - \sum_{i\in\mathcal{I}_r}p_i^k)$$

This is a kind of **shortest path problem**, where the distance between nodes $i$ and $j$ is the true distance $d_{ij}$, **minus** the dual variable associated with $i$. We can formulate it as follows:

$$\begin{align}
\min\quad & \sum_{i=1}^{N+1}\sum_{j=1}^{N+1}(d_{ij} - p_i^k)x_{ij}\\
\text{s.t.}\quad & \sum_{j=1}^{N+1}x_{ji} + z_i = \sum_{j=1}^{N+1}x_{ij}&\forall i\in[N]\\
& \sum_{j=1}^{N+1}x_{ij}\le 1 & \forall i \in [N]\\
& x_{ii}=0 & \forall i \in [N]\\
& \sum_{j=1}^{N+1}x_{j,N+1} = 1\\
& \sum_{j=1}^{N+1}x_{N+1,j} = 0\\
& \sum_{(i, j) \in C} x_{ij} \le |C| - 1 & \forall C \subseteq [N]\\
& x_{ij}\in\{0,1\}, z_i\in\{0,1\}\\
\end{align}$$

The first constraint is flow conservation, the second limits the flow through each customer to 1, the third removes self-edges, the fourth ensures the route visits the depot, and the fifth ensures the depot is terminal. The sixth constraint eliminates subtours (which by construction do not include the depot).

This subproblem is not that simple! In fact, the formulation sleuths among you will probably notice that if we change the right-hand side of the fourth constraint, we are actually formulating the full OVRP. However, this no longer holds for "harder" vehicle routing problems (e.g. CVRP).

We can implement this formulation below. Note that we need to include a callback to deal with the exponential number of subtour elimination constraints.


```julia
"""
    Given the induced graph as an adjacency list (i.e., next[i] is the next node to visit after node i),
        compute all subtours.
    Return them as a list of lists of nodes in the same component
"""
function find_subtours(next::Vector{Int})
    n = length(next)
    g = DiGraph(n)
    for i = 1:n
        if next[i] != 0
            add_edge!(g, i, next[i])
        end
    end
    components = strongly_connected_components(g)
    return sort(components, by=length)
end

"Find next column to add, returns route, cost, and whether or not the reduced cost is nonnegative"
function new_route(distances::Matrix, dual_variables::Vector)
    n = size(distances, 1) - 1 # number of customers
    @assert length(dual_variables) == n
    push!(dual_variables, 0) # depot
    model = Model(() -> Gurobi.Optimizer(GRB_ENV))
    set_optimizer_attributes(model, "TimeLimit" => 10, "OutputFlag" => 0)
    @variable(model, x[i=1:(n+1), j=1:(n+1)], Bin)
    @variable(model, z[i=1:n], Bin)
    @constraint(model, flow_conservation[i = 1:n],
                sum(x[j, i] for j = 1:(n+1)) + z[i] == sum(x[i, j] for j = 1:(n+1)))
    @constraint(model, max_flow[i = 1:n],
                sum(x[i, j] for j = 1:(n+1)) <= 1)
    @constraint(model, no_self_flow[i = 1:(n+1)], x[i, i] == 0)
    @constraint(model, into_depot, sum(x[j, n+1] for j = 1:(n+1)) == 1)
    @constraint(model, out_of_depot, sum(x[n+1, j] for j = 1:(n+1)) == 0)
    @objective(model, Min, sum(x[i, j] * (distances[i, j] - dual_variables[i]) for i = 1:(n+1), j=1:(n+1)))
        
    "Define the callback function"
    function eliminate_subtours(cb_data)
        status = callback_node_status(cb_data, model)
        if status == MOI.CALLBACK_NODE_STATUS_INTEGER
            # get value of current solution
            next = zeros(Int, n + 1)
            # your code here
            
            subtours = find_subtours(next)
            # solve dual subproblems
            for subtour in subtours
                if length(subtour) == 1 || n + 1 in subtour
                    continue
                else
                    # add subtour elimination cut
                    
                    
                end
            end
        end
    end
    # set callback function and attach to model
    MOI.set(model, MOI.LazyConstraintCallback(), eliminate_subtours)
    optimize!(model)
        
    route = [argmax(value.(z))]
    while route[end] != n+1
        push!(route, argmax(value.(x[route[end], :])))
    end
    reduced_cost = sum(distances[route[i], route[i+1]] - dual_variables[route[i]] for i = 1:(length(route)-1))
    return route, sum(distances[route[i], route[i+1]] for i = 1:(length(route)-1)), reduced_cost
end
```

Whew! We implemented the subproblem, which was not easy. Now, the main problem should be a piece of cake!


```julia
"Solve OVRP relaxation using column generation"
function ovrp(distances::Matrix; T::Int = 10, verbose::Bool = false)
    n = size(distances, 1) - 1 # number of customers
    routes, costs = initial_routes(dist)
    t = 0
    while true
        t += 1
        model = Model(() -> Gurobi.Optimizer(GRB_ENV))
        set_optimizer_attributes(model, "TimeLimit" => 10, "OutputFlag" => 0)
        # 
        # define restricted main problem
        
        
        optimize!(model)
        p = [-JuMP.shadow_price(visit_customer[i]) for i=1:n]
        route, cost, reduced_cost = new_route(distances, p)
        verbose && @printf("Found column with reduced cost: %.2g\n", reduced_cost)
        if reduced_cost > -1e-6 || t > T 
            return [routes[r] for r = eachindex(routes) if value(y[r]) >= 0.1],
                   sum(costs[r] for r = eachindex(routes) if value(y[r]) >= 0.1)
        end
        push!(routes, route); push!(costs, cost)
    end
end
```

Now, let's try to solve the problem, and look at the solution we obtain!


```julia
routes, total_cost = ovrp(dist, T=100, verbose = true)
println("Initial cost: ", sum(dist[1:N, N+1]))
println("Final cost: ", total_cost)
p = plot(size=(500,500), xlim=(0, L), ylim=(0, L))
p = scatter!(locations[1:N, 1], locations[1:N, 2], label="", legend=:topleft)
for route in routes
    p = plot!(locations[route, 1], locations[route, 2], label="", color="black")
end
p = scatter!(locations[[N+1], 1], locations[[N+1], 2], label="")
```

In a relatively small number of steps, we converge to this solution, which is optimal.

This is great news! Let's try a larger problem.


```julia
N = 20
L = 100
```


```julia
Random.seed!(1234)
locations = rand(N, 2) .* L
locations = vcat(locations, [L/2 L/2])
dist = [norm(locations[i, :] .- locations[j, :]) for i =1:N+1, j = 1:N+1];
plot(size=(500,500), xlim=(0, L), ylim=(0, L))
scatter!(locations[1:N, 1], locations[1:N, 2], label="Customers", legend=:topleft)
scatter!(locations[[N+1], 1], locations[[N+1], 2], label="Depot")
```

Now let's solve it using column generation!


```julia
routes, total_cost = ovrp(dist, T=2, verbose = true)
println("Initial cost: ", sum(dist[1:N, N+1]))
println("Final cost: ", total_cost)
p = plot(size=(500,500), xlim=(0, L), ylim=(0, L))
p = scatter!(locations[1:N, 1], locations[1:N, 2], label="", legend=:topleft)
for route in routes
    p = plot!(locations[route, 1], locations[route, 2], label="", color="black")
end
p = scatter!(locations[[N+1], 1], locations[[N+1], 2], label="")
```

It looks like we quickly pivot away from the bad initial solution, but then we get stuck on the single-route solution above, which is clearly suboptimal.

This is frustrating! What is the point of adding all these columns if we just end up picking the route we added at the very beginning?

When we add columns, we are also improving the **lower bound**. So even though we make no primal progress, we are always making dual progress!

Let's modify our OVRP solver to report upper and lower bounds at each iteration.

- **Upper bound:** this part is easy, we just report the objective of our current feasible solution.
- **Lower bound:** this part is trickier. Here, we need to use the fact that in the optimal solution, we will have no more than $N$ routes, so $\sum_{r\in\mathcal{R}}y^*_r \le N$. Then, given a feasible solution with objective $\hat{c}$ and the smallest reduced cost $\bar{c}$, the most improvement we can get on $\hat{c}$ is $$\hat{c} + \sum_{r\in\mathcal{R}}\bar{c}y^*_r \ge \hat{c} + \bar{c}N$$

(remember $\bar{c}$ is negative!)


```julia
"Solve OVRP relaxation using column generation"
function ovrp_progress(distances::Matrix; T::Int = 10, verbose::Bool = false)
    n = size(distances, 1) - 1 # number of customers
    routes, costs = initial_routes(dist)
    t = 0
    upper_bounds = []
    lower_bounds = []
    while true
        t += 1
        model = Model(() -> Gurobi.Optimizer(GRB_ENV))
        set_optimizer_attributes(model, "TimeLimit" => 10, "OutputFlag" => 0)
        @variable(model, y[r = eachindex(routes)] >= 0)
        @constraint(model, visit_customer[i=1:n],
                    sum(y[r] for r = eachindex(routes) if i in routes[r]) >= 1)
        @objective(model, Min, sum(costs[r] * y[r] for r = eachindex(routes)))
        optimize!(model)
        # 
        # save upper bound

        p = [-JuMP.shadow_price(visit_customer[i]) for i=1:n]
        route, cost, reduced_cost = new_route(distances, p)
        verbose && @printf("Found column with reduced cost: %.2g\n", reduced_cost)
        #
        # save lower bound
        
        if reduced_cost > -1e-6 || t > T 
            return [routes[r] for r = eachindex(routes) if value(y[r]) >= 0.1], objective_value(model),
                   upper_bounds, lower_bounds
        end
        push!(routes, route); push!(costs, cost)
    end
end
```

Now, let's try again on our 20-city instance!


```julia
routes, total_cost, upper_bounds, lower_bounds = ovrp_progress(dist, T=10, verbose = false)
plot(0:(length(upper_bounds) - 1), [upper_bounds, lower_bounds], legend=:bottomright, label=["Upper" "Lower"])
```

Notice that unlike Benders decomposition, there are no monotonic guarantees on the lower bound, but there is a monotonic guarantee on the upper bound. Convergence is not great!

## 3 - Solving CVRP using column generation

As we mentioned earlier, OVRP is kind of "easy" and using column generation is probably overkill. But we'll see that we can get a much bigger edge on the harder problem of CVRP.

We're going to consider a super-simple version of CVRP, where each customer has a demand of 1, and the capacity of the vehicle is $Q\ge 1$. Our initial routes remain feasible. 

All we need to do is modify the subproblem so that it only finds feasible routes, i.e. routes with less than $Q$ customers, or equivalently, less than $Q$ edges.


```julia
"Find next column to add, returns route, cost, and whether or not the reduced cost is nonnegative"
function new_route_cap(distances::Matrix, dual_variables::Vector, Q::Int=length(dual_variables))
    n = size(distances, 1) - 1 # number of customers
    @assert length(dual_variables) == n
    push!(dual_variables, 0) # depot
    model = Model(() -> Gurobi.Optimizer(GRB_ENV))
    set_optimizer_attributes(model, "TimeLimit" => 10, "OutputFlag" => 0)
    @variable(model, x[i=1:(n+1), j=1:(n+1)], Bin)
    @variable(model, z[i=1:n], Bin)
    @constraint(model, flow_conservation[i = 1:n],
                sum(x[j, i] for j = 1:(n+1)) + z[i] == sum(x[i, j] for j = 1:(n+1)))
    @constraint(model, max_flow[i = 1:n],
                sum(x[i, j] for j = 1:(n+1)) <= 1)
    @constraint(model, no_self_flow[i = 1:(n+1)], x[i, i] == 0)
    @constraint(model, into_depot, sum(x[j, n+1] for j = 1:(n+1)) == 1)
    @constraint(model, out_of_depot, sum(x[n+1, j] for j = 1:(n+1)) == 0)
    @objective(model, Min, sum(x[i, j] * (distances[i, j] - dual_variables[i]) for i = 1:(n+1), j=1:(n+1)))
        
    # CAPACITY CONSTRAINT
    
     
    "Define the callback function"
    function eliminate_subtours(cb_data)
        status = callback_node_status(cb_data, model)
        if status == MOI.CALLBACK_NODE_STATUS_INTEGER
            # get value of current solution
            next = zeros(Int, n + 1)
            for i = 1:(n+1), j=1:(n+1)
                if callback_value(cb_data, x[i, j]) > 0.5
                    next[i] = j
                end
            end
            subtours = find_subtours(next)
            # solve dual subproblems
            for subtour in subtours
                if length(subtour) == 1 || n + 1 in subtour
                    continue
                else
                    cut = @build_constraint(sum(x[subtour[i], subtour[i+1]] for i = 1:(length(subtour) - 1)) +
                                            x[subtour[length(subtour)], subtour[1]] <= length(subtour) - 1)
                    MOI.submit(model, MOI.LazyConstraint(cb_data), cut)
                end
            end
        end
    end
    # set callback function and attach to model
    MOI.set(model, MOI.LazyConstraintCallback(), eliminate_subtours)
    optimize!(model)
    
    next = zeros(Int, n + 1)
    for i = 1:(n+1), j=1:(n+1)
        if value(x[i, j]) > 0.5
            next[i] = j
        end
    end
    
    route = [argmax(value.(z))]
    while route[end] != n+1
        push!(route, argmax(value.(x[route[end], :])))
    end
    reduced_cost = sum(distances[route[i], route[i+1]] - dual_variables[route[i]] for i = 1:(length(route)-1))
    return route, sum(distances[route[i], route[i+1]] for i = 1:(length(route)-1)), reduced_cost
end
```

And we correspondingly change the main problem:


```julia
"Solve CVRP relaxation using column generation"
function cvrp_progress(distances::Matrix, Q::Int; T::Int = 10, verbose::Bool = false)
    n = size(distances, 1) - 1 # number of customers
    routes, costs = initial_routes(dist)
    t = 0
    upper_bounds = []
    lower_bounds = []
    while true
        t += 1
        model = Model(() -> Gurobi.Optimizer(GRB_ENV))
        set_optimizer_attributes(model, "TimeLimit" => 10, "OutputFlag" => 0)
        @variable(model, y[r = eachindex(routes)] >= 0)
        @constraint(model, visit_customer[i=1:n],
                    sum(y[r] for r = eachindex(routes) if i in routes[r]) >= 1)
        @objective(model, Min, sum(costs[r] * y[r] for r = eachindex(routes)))
        optimize!(model)
        push!(upper_bounds, objective_value(model))
        p = [-JuMP.shadow_price(visit_customer[i]) for i=1:n]
        route, cost, reduced_cost = new_route_cap(distances, p, Q)
        verbose && @printf("Found column with reduced cost: %.2g\n", reduced_cost)
        push!(lower_bounds, objective_value(model) + n * reduced_cost)
        if reduced_cost > -1e-6 || t > T 
            return [routes[r] for r = eachindex(routes) if value(y[r]) >= 0.1], objective_value(model),
                   upper_bounds, lower_bounds
        end
        push!(routes, route); push!(costs, cost)
    end
end
```

Let's try it on our 20-city example.


```julia
Q = 4
@time routes, total_cost, upper, lower = cvrp_progress(dist, Q, T=100, verbose = true)
println("Initial cost: ", sum(dist[1:N, N+1]))
println("Final cost: ", total_cost)
```

Let's look at the progress curve:


```julia
plot(0:(length(upper) - 1), [upper, lower], legend=:bottomright, label=["Upper" "Lower"])
```

And let's observe the solution


```julia
p = plot(size=(500,500), xlim=(0, L), ylim=(0, L))
p = scatter!(locations[1:N, 1], locations[1:N, 2], label="", legend=:topleft)
for route in routes
    p = plot!(locations[route, 1], locations[route, 2], label="", color="black")
end
p = scatter!(locations[[N+1], 1], locations[[N+1], 2], label="")
```

This solution makes sense, but we do observe some funky routes. Recall that all this time, we have been solving the LP relaxation! So we are paying the price of keeping all routes with a nonzero (but possibly fractional) selection variable $y_r$.

## 4 - Restoring integrality

That was all very nice, but we still don't have an integral solution! What should we do about that?

Two approaches:
- integrate the pricing subproblem into the branch and bound tree. This is called "branch and price", it is the state of the art, but also hard to implement. We will discuss in lecture next week.
- use heuristics!

The simplest heuristic we can use is "price and branch". We simply keep adding columns to the relaxation until we no longer can, then we solve an integer program over the columns we've added.

We note that this heuristic will produce a feasible solution because of the initial routes we've chosen, but in general, it may not find a feasible solution!

Let's implement it for CVRP:


```julia
"Solve OVRP relaxation using column generation"
function cvrp_progress_integral(distances::Matrix, Q::Int; T::Int = 10, verbose::Bool = false)
    n = size(distances, 1) - 1 # number of customers
    routes, costs = initial_routes(dist)
    t = 0
    upper_bounds = []
    lower_bounds = []
    while true
        t += 1
        model = Model(() -> Gurobi.Optimizer(GRB_ENV))
        set_optimizer_attributes(model, "TimeLimit" => 10, "OutputFlag" => 0)
        @variable(model, y[r = eachindex(routes)] >= 0)
        @constraint(model, visit_customer[i=1:n],
                    sum(y[r] for r = eachindex(routes) if i in routes[r]) >= 1)
        @objective(model, Min, sum(costs[r] * y[r] for r = eachindex(routes)))
        optimize!(model)
        push!(upper_bounds, objective_value(model))
        p = [-JuMP.shadow_price(visit_customer[i]) for i=1:n]
        route, cost, reduced_cost = new_route_cap(distances, p, Q)
        verbose && @printf("Found column with reduced cost: %.2g\n", reduced_cost)
        push!(lower_bounds, objective_value(model) + n * reduced_cost)
        if reduced_cost > -1e-6 || t > T
            model = Model(() -> Gurobi.Optimizer(GRB_ENV))
            set_optimizer_attributes(model, "TimeLimit" => 10, "OutputFlag" => 0)
            # 
            # define new model with integer variables
            
            optimize!(model)
            return [routes[r] for r = eachindex(routes) if value(y[r]) >= 0.1], objective_value(model),
                   upper_bounds, lower_bounds
        end
        push!(routes, route); push!(costs, cost)
    end
end
```


```julia
Q = 4
@time routes, total_cost, upper, lower = cvrp_progress_integral(dist, Q, T=100, verbose = true)
println("Initial cost: ", sum(dist[1:N, N+1]))
println("Final cost: ", total_cost)
println("Gap: ", (total_cost - lower[end]) / lower[end])
```

We see that we terminated with a tiny gap, so this heuristic worked really well! Let's look at it plotted below:


```julia
plot(0:(length(upper) - 1), [upper, lower], legend=:bottomright, label=["Upper" "Lower"])
scatter!([length(upper) - 1], [total_cost], label="Integer solution")
```

And we can take a look at the solution:


```julia
p = plot(size=(500,500), xlim=(0, L), ylim=(0, L))
p = scatter!(locations[1:N, 1], locations[1:N, 2], label="", legend=:topleft)
for route in routes
    p = plot!(locations[route, 1], locations[route, 2], label="", color="black")
end
p = scatter!(locations[[N+1], 1], locations[[N+1], 2], label="")
```

## 5 - Large-scale neighborhood search

In the last part of this recitation, we'll talk a little more about heuristics, since this is often a great way to solve practical problems.

The one heuristic I want to talk about is large-scale neighborhood search (LNS). It has a fancy name, but it's really the simplest idea in the world:

> My MIP doesn't solve? OK, I'll just fix most of the variables. Then I'll keep changing a few at a time until I converge.

Let's apply this idea to the CVRP. Recall that its cousin the TSP has a simple formulation with a polynomial number of constraints, which doesn't scale very well. Let's write down the standard MTZ formulation for CVRP:


```julia
function cvrp_mtz(distances::Matrix, Q::Int; verbose::Bool=true)
    n = size(distances, 1) - 1 # number of customers
    model = Model(() -> Gurobi.Optimizer(GRB_ENV))
    set_optimizer_attributes(model, "TimeLimit" => 5, "OutputFlag" => ifelse(verbose, 1, 0))
    @variable(model, x[1:(n+1), 1:(n+1)], Bin)
    # Minimize total distance, but don't count distance leaving the depot
    @objective(model, Min,
               sum(dist[i, j] * x[i, j] for i = 1:n, j = 1:(n+1)))
    @constraint(model, one_successor[i=1:n],
                sum(x[i, j] for j=1:(n+1)) == 1)
    @constraint(model, flow_conservation[i=1:(n+1)],
                sum(x[i, j] for j=1:(n+1)) == sum(x[j, i] for j=1:(n+1)))
    @constraint(model, no_self_edges[i=1:(n+1)],
                x[i, i] == 0)
    # variable designates the number of customers that have been visited before this one
    @variable(model, 0 <= u[i=1:n] <= Q)
    @constraint(model, [i=1:n, j=1:n; j != i], u[j] >= u[i] + 1 - (Q) * (1 - x[i, j]))
    # 
    # warm start the solution with our initial routes, as before
    
    optimize!(model)
    x_val = value.(x)
    routes = []
    while sum(x_val) > 1e-6
        current = findfirst(a -> a > 0.5, x_val[n+1, :])
        x_val[n+1, current] = 0.0
        route = []
        while true
            push!(route, current)
            current == n+1 && break
            next = findfirst(a -> a > 0.5, x_val[current, :])
            x_val[current, next] = 0.0
            current = next
        end
        push!(routes, route)
    end
    return routes, objective_value(model)
end 
```

Let's solve our 20-customer instance.


```julia
routes, cost = cvrp_mtz(dist, 4)
println("Total cost: ", cost)
```

Looks like we did pretty well! Let's visualize the routes.


```julia
p = plot(size=(500,500), xlim=(0, L), ylim=(0, L))
p = scatter!(locations[1:N, 1], locations[1:N, 2], label="", legend=:topleft)
for route in routes
    p = plot!(locations[route, 1], locations[route, 2], label="", color="black")
end
p = scatter!(locations[[N+1], 1], locations[[N+1], 2], label="")
```

This matches our price-and-branch solution from earlier. Can we scale up?


```julia
N = 40

Random.seed!(1234)
locations = rand(N, 2) .* L
locations = vcat(locations, [L/2 L/2])
dist = [norm(locations[i, :] .- locations[j, :]) for i =1:N+1, j = 1:N+1];
plot(size=(500,500), xlim=(0, L), ylim=(0, L))
scatter!(locations[1:N, 1], locations[1:N, 2], label="Customers", legend=:topleft)
scatter!(locations[[N+1], 1], locations[[N+1], 2], label="Depot")
```


```julia
routes, total_cost = cvrp_mtz(dist, 4)
println("Total cost: ", total_cost)
```

Let's visualize:


```julia
p = plot(size=(500,500), xlim=(0, L), ylim=(0, L))
p = scatter!(locations[1:N, 1], locations[1:N, 2], label="", legend=:topleft)
for route in routes
    p = plot!(locations[route, 1], locations[route, 2], label="", color="black")
end
p = scatter!(locations[[N+1], 1], locations[[N+1], 2], label="")
```

This is visibly suboptimal (for example, see the left quadrant) and the gap is 30%.

We would like to "tweak" the solution by "destroying" some of the edges, then "repairing" the solution, hopefully with lower cost!

We can destroy the routes by removing $k$ edges, as follows:


```julia
function split_routes(route_list::Vector, k::Int)
    num_edges = [length(r) - 1 for r in route_list]
    total_edges = sum(num_edges)
    cum_edges = cumsum(num_edges)
    splits = sort(shuffle(1:total_edges)[1:k])
    routes_to_split = [findfirst(x -> x >= s, cum_edges) for s in splits]
    places_to_split = [splits[i] - (routes_to_split[i] == 1 ? 0 : cum_edges[routes_to_split[i] - 1])
                       for i = eachindex(splits)]
    new_routes = []
    for (r, route) in enumerate(routes)
        if r in routes_to_split
            where_to_split = [places_to_split[i] for i = eachindex(places_to_split) if routes_to_split[i] == r]
            start = 1
            for j in where_to_split
                push!(new_routes, route[start:j])
                start = j + 1
            end
            push!(new_routes, route[start:end])
        else
            push!(new_routes, route)
        end
    end
    return new_routes
end
```

For example, we can remove 7 edges.


```julia
Random.seed!(1234)
fragments = split_routes(routes, 7)
```

The "fragments" look like this:


```julia
p = plot(size=(500,500), xlim=(0, L), ylim=(0, L))
p = scatter!(locations[1:N, 1], locations[1:N, 2], label="", legend=:topleft)
for route in fragments
    p = plot!(locations[route, 1], locations[route, 2], label="", color="black")
end
p = scatter!(locations[[N+1], 1], locations[[N+1], 2], label="")
```

We can then tweak our CVRP formulation to keep fragments intact and simply connect them optimally.


```julia
function repair_cvrp_mtz(distances::Matrix, Q::Int, fragments::Vector; verbose::Bool=true)
    n = size(distances, 1) - 1 # number of customers
    model = Model(() -> Gurobi.Optimizer(GRB_ENV))
    set_optimizer_attributes(model, "TimeLimit" => 5, "OutputFlag" => ifelse(verbose, 1, 0))
    @variable(model, x[1:(n+1), 1:(n+1)], Bin)
    # Minimize total distance, but don't count distance leaving the depot
    @objective(model, Min,
               sum(dist[i, j] * x[i, j] for i = 1:n, j = 1:(n+1)))
    @constraint(model, one_successor[i=1:n],
                sum(x[i, j] for j=1:(n+1)) == 1)
    @constraint(model, flow_conservation[i=1:(n+1)],
                sum(x[i, j] for j=1:(n+1)) == sum(x[j, i] for j=1:(n+1)))
    @constraint(model, no_self_edges[i=1:(n+1)],
                x[i, i] == 0)
    # variable designates the number of customers that have been visited before this one
    @variable(model, 0 <= u[i=1:n] <= Q)
    @constraint(model, [i=1:n, j=1:n; j != i], u[j] >= u[i] + 1 - (Q) * (1 - x[i, j]))
    
    # 
    # fix our route fragments - your code here
    
    
    optimize!(model)
    x_val = value.(x)
    routes = []
    while sum(x_val) > 1e-6
        current = findfirst(a -> a > 0.5, x_val[n+1, :])
        x_val[n+1, current] = 0.0
        route = []
        while true
            push!(route, current)
            current == n+1 && break
            next = findfirst(a -> a > 0.5, x_val[current, :])
            x_val[current, next] = 0.0
            current = next
        end
        push!(routes, route)
    end
    return routes, objective_value(model)
end 
```

Let's try to repair our fragments...


```julia
routes, cost = repair_cvrp_mtz(dist, 4, fragments)
println("Total cost: ", cost)
```


```julia
p = plot(size=(500,500), xlim=(0, L), ylim=(0, L))
p = scatter!(locations[1:N, 1], locations[1:N, 2], label="", legend=:topleft)
for route in routes
    p = plot!(locations[route, 1], locations[route, 2], label="", color="black")
end
p = scatter!(locations[[N+1], 1], locations[[N+1], 2], label="")
```

Looks like nothing changed, but we can do this iteratively!


```julia
Q = 4
routes, total_cost = cvrp_mtz(dist, Q, verbose=false)
println("Cost: ", total_cost)
Random.seed!(1234)
for t = 1:20
    fragments = split_routes(routes, 10)
    routes, total_cost = repair_cvrp_mtz(dist, Q, fragments, verbose=false)
    println("Cost: ", total_cost)
end
```

There is a tradeoff here: the more edges we remove, the larger the neighborhood we have to explore, and the more we will gain from an iteration, but the harder it becomes to optimize over the neighborhood.

850f237f550b3e9a3a9109e6620625f4d82f1db5
ipynb
Jupyter Notebook
Rec10_code/Recitation-10-Student.ipynb
leobix/mio2021
ec28e15b9f701745b3616ee9174d3573fbac0d14
MIT
Rec10_code/Recitation-10-Student.ipynb
leobix/mio2021
ec28e15b9f701745b3616ee9174d3573fbac0d14
MIT
Rec10_code/Recitation-10-Student.ipynb
leobix/mio2021
ec28e15b9f701745b3616ee9174d3573fbac0d14
MIT
Qwen/Qwen-72B
1. YES
2. YES