*Santa Claus needs to visit 15 chimneys as quickly as possible to deliver presents to the children who have been good in 2020 and stuck to social distancing rules.*

*Can you help him find the shortest route crossing each home exactly once and return to the first home to then go back to the North Pole from there?*

*Santa has asked you to obtain the best possible solution to this problem.*

Since Santa has asked for the best solution, we need to use an exact solution technique. We can therefore proceed to implement a mathematical model that implements the Miller, Tucker and Zemlin subtour elimination constraints (SECs).

We proceed to implement the model using PuLP. For the purposes of this notebook, we obtain the locations of the chimneys using our trusty blob maker.

In [1]:

```
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import seaborn as sns; sns.set()
from sklearn.datasets import make_blobs
%matplotlib inline
plot_size = 15
plot_width = 16
plot_height = 8
params = {'legend.fontsize': 'large',
'figure.figsize': (plot_width,plot_height),
'axes.labelsize': plot_size,
'axes.titlesize': plot_size,
'xtick.labelsize': plot_size*0.75,
'ytick.labelsize': plot_size*0.75,
'axes.titlepad': 25}
plt.rcParams.update(params)
```

In [2]:

```
num_cities = 15
```

In [3]:

```
center_box = (100, 200)
cities_coord, _ = make_blobs(n_samples=num_cities,
centers=2,
cluster_std=20,
center_box=center_box,
random_state = 2)
```

In [4]:

```
plt.scatter(cities_coord[:, 0], cities_coord[:, 1], s=plot_size*2, cmap='viridis');
```

In [5]:

```
from pulp import *
prob = LpProblem('prob', LpMinimize)
```

We are back again in using `PuLP`

.

In [6]:

```
from scipy.spatial import distance
dist_matrix = distance.cdist(cities_coord, cities_coord, 'euclidean')
dist_matrix.shape
```

Out[6]:

(15, 15)

In [7]:

```
N = [('N'+str(i+1)) for i in range(num_cities)]
```

In [8]:

```
c = {n_f:{n_t: round(dist_matrix[i][j],2) for j,n_t in enumerate(N)}
for i,n_f in enumerate(N)}
```

In [9]:

```
x = LpVariable.dicts('x',(N,N), 0,1,LpBinary)
```

In [10]:

```
prob += lpSum([lpSum([x[i][j]*c[i][j] for j in N if i!=j])
for i in N])
```

In [11]:

```
for i in N:
prob += lpSum([x[i][j] for j in N if i!=j]) == 1
for j in N:
prob += lpSum([x[i][j] for i in N if i!=j]) == 1
```

In [12]:

```
%time status = prob.solve()
print(f'\nSTATUS\n{LpStatus[status]}\n')
```

In [13]:

```
def plot_solution(_N, _x, _cities_coord):
vect_orig = []
vect_dest = []
u = v = np.zeros((len(vect_orig),len(vect_orig)))
for j,n_t in enumerate(_N):
for i,n_f in enumerate(_N):
if i!=j:
if _x[n_f][n_t].varValue > 0:
vect_orig.append(_cities_coord[i])
vect_dest.append(_cities_coord[j])
vect_orig = np.array(vect_orig)
vect_dest = np.array(vect_dest)
plt.scatter(_cities_coord[:, 0],
_cities_coord[:, 1],
s=plot_size*2,
cmap='viridis',
zorder = 10000);
for i in range(len(vect_orig)):
plt.arrow(vect_orig[i][0],
vect_orig[i][1],
vect_dest[i][0]-vect_orig[i][0],
vect_dest[i][1]-vect_orig[i][1],
color='black')
plot_solution(N, x, cities_coord)
```

In [14]:

```
for j,n_t in enumerate(N):
for i,n_f in enumerate(N):
if i!=j and x[n_f][n_t].varValue > 0:
print(f'{n_f} -> {n_t}')
```

In [15]:

```
from pulp import *
prob = LpProblem('prob', LpMinimize)
N = [('N'+str(i+1)) for i in range(num_cities)]
x = LpVariable.dicts('x',(N,N), 0,1,LpBinary)
c = {n_f:{n_t: dist_matrix[i][j] for j,n_t in enumerate(N)}
for i,n_f in enumerate(N)}
prob += lpSum([lpSum([x[i][j]*c[i][j] for j in N if i!=j])
for i in N])
for i in N:
prob += lpSum([x[i][j] for j in N if i!=j]) == 1
for j in N:
prob += lpSum([x[i][j] for i in N if i!=j]) == 1
```

In [16]:

```
for i in N:
for j in N:
if i!=j:
prob += x[i][j] + x[j][i] <= 1
```

In [17]:

```
%time status = prob.solve()
print(f'\nSTATUS\n{LpStatus[status]}\n')
```

In [18]:

```
plot_solution(N, x, cities_coord)
```

In [19]:

```
prob = LpProblem('prob', LpMinimize)
N = [('N'+str(i+1)) for i in range(num_cities)]
c = {n_f:{n_t: round(dist_matrix[i][j],2) for j,n_t in enumerate(N)} for i,n_f in enumerate(N)}
x = LpVariable.dicts('x',(N,N), 0,1,LpBinary)
prob += lpSum([lpSum([x[i][j]*c[i][j] for j in N if i!=j]) for i in N])
for i in N:
prob += lpSum([x[i][j] for j in N if i!=j]) == 1
for j in N:
prob += lpSum([x[i][j] for i in N if i!=j]) == 1
```

In [20]:

```
#we need to keep track of the order in the tour to eliminate the possibility of subtours
u = LpVariable.dicts('u', N, 0, num_cities-1, LpInteger)
```

In [21]:

```
U=num_cities
for i in N:
for j in N:
if i != j and (i != 'N1' and j!= 'N1'):
prob += u[i] - u[j] <= U*(1-x[i][j]) - 1
```

In [22]:

```
%time status = prob.solve()
print(f'\nSTATUS\n{LpStatus[status]}\n')
```

In [23]:

```
plot_solution(N, x, cities_coord)
```

In [24]:

```
N_unvisited = N.copy()
current = 'N1'
tour=[]
tour.append(N_unvisited.pop(N_unvisited.index(current)))
while len(N_unvisited) > 0:
for k in N_unvisited:
if x[current][k].varValue ==1:
tour.append( N_unvisited.pop(N_unvisited.index(k)))
current=k
break
tour.append('N1')
print('Optimal tour!')
print(' -> '.join(tour))
```