```julia
using JuMP, Cbc, Plots     # The usual packages
using DelimitedFiles       # IO reading/writing files
using LinearAlgebra        # For convenience (includes function norm(x))
using Combinatorics        # To perform permutations
```

# MS-E2121 - Linear optimization
## Exercise session 8



### Demo exercise: Travelling salesman problem (TSP) - MTZ formulation
#### Examining the formulation

Show that the following formulation $P_{MTZ}$ is valid for the TSP defined on a directed graph $G = (N,A)$ with $N = \{1,\dots,n\}$ cities and arcs $A = \{(i,j) : i,j\in N, i\neq j\}$ between cities.
$$
P_{MTZ} = \left\{
   \begin{array}{ll}
   \displaystyle \sum_{j \in N \setminus \{i\}} x_{ij} = 1, & \forall i \in N \\
   \displaystyle \sum_{j \in N \setminus \{i\}} x_{ji} = 1, & \forall i \in N \\   
   \displaystyle u_{i} - u_{j} + (n-1) x_{ij} \leq n - 2, & \forall i,j \in N \setminus \{ 1 \} : i \neq j ~~(*)\\
    x_{ij} \in \{0,1\}, & \forall i,j \in N : i\neq j\\
  \end{array}
\right.
$$

where $x_{ij} = 1$ if city $j\in N$ is visited immediately after city $i \in N$, and $x_{ij} = 0$ otherwise. Constraints $(*)$ with the variables $u_i \in \mathbb{R}$ for all $i\in N$ are called *Miller-Tucker-Zemlin* (MTZ) subtour elimination constraints.


We want to show that
1. Constraints $(*)$ prevent subtours in any solution $x \in P_{MTZ}$.
2. Every TSP solution $x$ satisfies the constraints $(*)$.

To prove the first one, we first assume that a solution $x \in P_{MTZ}$ *has* a subtour with $k$ nodes and $k$ arcs between them, not going through node 1. For example, assume nodes 2, 3 and 5 form a subtour when $n=5$. Let's write the constraints $(*)$ corresponding to this subtour:

\begin{align}
    u_2 - u_3 + 4 \le 3 \\
    u_3 - u_5 + 4 \le 3 \\
    u_5 - u_2 + 4 \le 3
\end{align}

We observe that if $x_{ij} = 1$ such that $i,j \in N \setminus 1$, the constraint $(*)$ can be written as $u_i \le u_j - 1$, which for integer variables is the same as $u_i < u_j$. For a general result, we denote the nodes in the subtour by $\{i_1, ..., i_k\}$ and get $u_{i_1} < u_{i_k} < u_{i_1}$, which is a contradiction. This tells us that there can be no subtour ($k<n$) that doesn't contain node 1. A subtour containing node 1 would imply another subtour not containing node 1, so that is also forbidden by $(*)$. This proves the first part.



For the second part, we notice that the $u$-variables seem to imply an ordering for the nodes. Assume that all tours start from node 1 and $u_i, \ i \in N \setminus 1$ is the position of the node on the tour (the first node visited after the starting node 1 has $u$-value 2, the second one 3 and so on). For each arc $i \rightarrow j$, we have either $x_{ij}=0$ or $x_{ij}=1$.

If $x_{ij} = 0$: there is no arc from $i$ to $j$, and the constraint is $u_{i} - u_{j} \leq n - 2$, which holds, since we defined that $i$ and $j$ can't be 1, and also not greater than $n$. The upper bound of the difference between two $u$-values is thus $n-2$, and the constraint always holds if $x_{ij}=0$.

If $x_{ij} = 1$: there is an arc from $i$ to $j$, and the constraint is $u_{i} - u_{j} + n-1 \leq n - 2$ or $u_{i} - u_{j} \leq -1$, which looks familiar from before and actually holds when we do not have subtours since by our definition, $u_i - u_j = -1$. The key is that this time, we have no subtours that do not contain node 1, which is treated as a special case.

To conclude, we proved that this formulation does not allow subtours, and all valid solutions $x$ satisfy the constraints. Therefore, the formulation is valid.

#### MTZ and naive TSP implementation
We first write some helper functions, starting with one that computes the distances between coordinates:


```julia
## Function for getting the distances array
function get_dist(xycoord::Matrix{},n::Int)
    # Compute distance matrix (d[i,j]) from city coordinates
    
    dist = zeros(n,n)
    for i = 1:n
        for j = i:n
            d = norm(xycoord[i,:] - xycoord[j,:])
            dist[i,j] = d
            dist[j,i] = d
        end
    end
    return dist
end
```

A function to convert the adjacency matrix $x$ describing a tour to a vector representing the tour, starting from city 1.


```julia
# Get the optimal tour
# Input
#     x: solution matrix 
#     n: number of cities
# Returns
#     tour: ordering of cities in the optimal tour
function gettour(x::Matrix{Int}, n::Int)
    tour = zeros(Int,n+1)   # Initialize tour vector (n+1 as city 1 appears twice)
    tour[1] = 1             # Set city 1 as first one in the tour
    k = 2                   # Index of vector tour[k]
    i = 1                   # Index of current city 
    while k <= n + 1        # Find all n+1 tour nodes (city 1 is counted twice)
        for j = 1:n         
            if x[i,j] == 1  # Find next city j visited immediately after i
                tour[k] = j # Set city j as the k:th city in the tour
                k = k + 1   # Update index k of tour[] vector
                i = j       # Move to next city
                break       
            end
        end
    end 
    return tour             # Return the optimal tour 
end
```


```julia
## Defining the colors to be used
c_blue = palette(:auto)[1]                         # color :1
c_orange = palette(:auto)[2]                       # color :2
c_green = palette(:auto)[3];                       # color :3
```

The data is stored in csv-files, a common external data format. We use the DelimitedFiles package to read a file into a matrix ```data```.


```julia
# "data16a.csv" has 3 columns which are stored in (Nullable) Arrays
# data[1], data[2], data[3] after the function call
# data = CSV.read(...) below. The columns contain:
#
#     data[:,1]: all cities i in V
#     data[:,2]: x-coordinate of each city i in V
#     data[:,3]: y-coordinate of each city i in V
#

data  = readdlm("data16a.csv", ',')   
# data  = readdlm("data16b.csv", ',')   
n = 16                  # number of cities 
# println(data)         # Look at the data in compact form
V = data[2:n+1,1]       # All cities i in V
x = data[2:n+1,2]       # x-coordinates of cities i in V
y = data[2:n+1,3]       # y-coordinates of cities i in V
xycoord = [x y];        # n x 2 coordinate matrix
```


```julia
function tsp_naive(xycoord::Matrix{}, n::Int)
    
    # Create a model 
    m = Model(Cbc.Optimizer)  
    
    # Here the costs c are the distances between cities
    c = get_dist(xycoord,n)
    
    ## Variables

    # x[i,j] = 1 if we travel from city i to city j, 0 otherwise.
    @variable(m, x[1:n,1:n], Bin)
    
    ## Objective
    
    # Minimize length of tour
    @objective(m, Min, dot(c,x))

    ## Constraints
    
    # Ignore self arcs: set x[i,i] = 0  
    @constraint(m, sar[i = 1:n], x[i,i] == 0)

    # We must enter and leave every city exactly once             
    @constraint(m, ji[i = 1:n], sum(x[j,i] for j = 1:n if j != i) == 1)
    @constraint(m, ij[i = 1:n], sum(x[i,j] for j = 1:n if j != i) == 1)

    optimize!(m)
                                
    cost = objective_value(m)         # Optimal cost (length)
    sol_x = round.(Int, value.(x))    # Optimal solution vector
    
    return m, sol_x, cost
end;
```


```julia
## Solve the problem and evaluate time and memory with @time macro
(m_naive, x_naive, cost_naive) = @time tsp_naive(xycoord, n);

# Get the optimal tour
tour_naive = gettour(x_naive,n);
```


```julia
plt = scatter(xycoord[:,1],xycoord[:,2], 
      markercolor = c_blue,
      markerstrokewidth = 0,
      legend = false
)
for i in 1:length(tour_naive)-1
    annotate!(xycoord[tour_naive[i],1]+50, xycoord[tour_naive[i],2]+50, ("$(tour_naive[i])", 7, :left))
    plot!(([xycoord[tour_naive[i],1],xycoord[tour_naive[i+1],1]] , [xycoord[tour_naive[i],2],xycoord[tour_naive[i+1],2]]), c = c_blue, label = "")
end
plt
```


```julia
plt = scatter(xycoord[:,1],xycoord[:,2], 
      markercolor = c_blue,
      markerstrokewidth = 0,
      legend = false
)
for i in 1:n
    annotate!(xycoord[i,1]+50, xycoord[i,2]+50, ("$(i)", 7, :left))
    for j in 1:n
        if x_naive[i,j] == 1
            plot!(([xycoord[i,1],xycoord[j,1]] , [xycoord[i,2],xycoord[j,2]]), c = c_blue, label = "")
            break
        end
    end
end
plt
```


```julia
# Solve a directed, TSP instance (MTZ formulation)
# Input
#    xycoord: coordinates of city locations
#    n: number of cities
# Returns
#    tour: ordering of cities in the optimal tour
#    cost: cost (length) of the optimal tour
function tsp_mtz(xycoord::Matrix{}, n::Int)
    
    # Create a model 
    m = Model(Cbc.Optimizer)  
    
    # Here the costs c are the distances between cities
    c = get_dist(xycoord,n)
    
    ## Variables

    # x[i,j] = 1 if we travel from city i to city j, 0 otherwise.
    @variable(m, x[1:n,1:n], Bin)
    # Variables u for subtour elimination constraints
    @variable(m, u[2:n])   
    
    ## Objective
    
    # Minimize length of tour
    @objective(m, Min, dot(c,x))

    ## Constraints
    
    # Ignore self arcs: set x[i,i] = 0  
    @constraint(m, sar[i = 1:n], x[i,i] == 0)

    # We must enter and leave every city exactly once             
    @constraint(m, ji[i = 1:n], sum(x[j,i] for j = 1:n if j != i) == 1)
    @constraint(m, ij[i = 1:n], sum(x[i,j] for j = 1:n if j != i) == 1)
                                    
    # MTZ subtour elimination constraints
    @constraint(m, sub[i = 2:n, j = 2:n, i != j], u[i] - u[j] + (n-1)*x[i,j] <= (n-2))
                                
    optimize!(m)
                                
    cost = objective_value(m)         # Optimal cost (length)
    sol_x = round.(Int, value.(x))    # Optimal solution vector
    
    return m, sol_x, cost
end;
```


```julia
## Solve the problem and evaluate time and memory with @time macro
(m_mtz, x_mtz, cost_mtz) = @time tsp_mtz(xycoord, n);

# Get the optimal tour
tour_mtz = gettour(x_mtz,n);      
```


```julia
## Print the optimal tour and its cost
println("\nOptimal tour: $(tour_mtz)\n")
println("Optimal length: ", cost_mtz)
```


```julia
plt = scatter(xycoord[:,1],xycoord[:,2], 
      markercolor = c_blue,
      markerstrokewidth = 0,
      legend = false
)
for i in 1:length(tour_mtz)-1
    annotate!(xycoord[tour_mtz[i],1]+50, xycoord[tour_mtz[i],2]+50, ("$(tour_mtz[i])", 7, :left))
    plot!(([xycoord[tour_mtz[i],1],xycoord[tour_mtz[i+1],1]] , [xycoord[tour_mtz[i],2],xycoord[tour_mtz[i+1],2]]), c = c_blue, label = "")
end
plt
```

Now that we have solved the problem using the MTZ formulation, let's try solving the same problem starting with the naive implementation and using successive cutset or subtour elimination constraints:


```julia
## Initialisation
(m_naive, x_naive, cost_naive) = tsp_naive(xycoord, n);
subnodes = []
count = 0
stop = 0
lim = 100
methods = [:elimination,:cutset]
method = methods[2];
```


```julia
## Perform cuts to break subtours until we got an optimal
@time while stop == 0 && count < lim

    S = collect(permutations(subnodes,2))     # Possible connections present in the naive implementation
    NS = setdiff(V,subnodes)                  # Nodes that are still not included in the tour

    if method == :cutset
        ## Cutset constraints
        if length(S) > 0
            @constraint(m_naive,sum(m_naive[:x][subnodes[i],NS[j]] for i in 1:length(subnodes), j in 1:length(NS)) >= 1)
        end
    else
        ## Subtour elimination constraints
        if length(S) > 0
            @constraint(m_naive,sum(m_naive[:x][S[i][1],S[i][2]] for i in 1:length(S)) <= length(subnodes)-1)
        end
    end

    set_silent(m_naive)
    optimize!(m_naive)

    cost2 = objective_value(m_naive)          # Optimal cost (length)
    sol_x = round.(Int, value.(m_naive[:x]))     # Optimal solution vector

    tour2 = gettour(sol_x,n)            # Get the optimal tour
    
    if length(unique(tour2)) < n
        count = count + 1
        println("Method used: ",method)
        println("Subtours present in node 1, update subnodes vector to break the subtour: ", tour2')
        println("\nIteration $(count); not optimal.\n")
        subnodes = unique(tour2);
    else
        println("Method used: ",method)
        println("Optimal tour: ", tour2')
        println("\n Took $(count) iterations to find the optimal solution.")
        S = []
        stop = 1
    end;
end;
```

### Scheduling problem
TODO

7570b86b6a6769030f2b25506517f7e51ba35efc
ipynb
Jupyter Notebook
exercises/ex-08/exercise8_skeleton.ipynb
nnhjy/lp-julia
06c63cf06b6911bc4ab99a8c18b33bca2b61083c
MIT
exercises/ex-08/exercise8_skeleton.ipynb
nnhjy/lp-julia
06c63cf06b6911bc4ab99a8c18b33bca2b61083c
MIT
exercises/ex-08/exercise8_skeleton.ipynb
nnhjy/lp-julia
06c63cf06b6911bc4ab99a8c18b33bca2b61083c
MIT
Qwen/Qwen-72B
1. YES
2. YES