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 = 30
plot_width = 19
plot_height = 17
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_customer = 40
num_facility = 5
```

In [3]:

```
center_box = (100, 200)
cust_coord, _ = make_blobs(n_samples=num_customer,
centers=2,
cluster_std=20,
center_box=center_box,
random_state = 2)
facil_coord, _ = make_blobs(n_samples=num_facility,
centers=1,
cluster_std=30,
center_box=center_box,
random_state = 50)
```

In [4]:

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

In [5]:

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

In [6]:

```
from scipy.spatial import distance
dist_matrix = distance.cdist(cust_coord, facil_coord, 'euclidean')
```

In [7]:

```
S = [('S'+str(i+1)) for i in range(num_facility)]
D = [('C'+str(i+1)) for i in range(num_customer)]
```

In [8]:

```
n = 2
```

We are using the above parameter variable to set the number of facilities to be opened.

In [9]:

```
d = {customer:{facility:dist_matrix[i][j] for j,facility in enumerate(S)} for i,customer in enumerate(D)}
```

In [10]:

```
x = LpVariable.dicts('x', (D,S), lowBound = 0, upBound = 1, cat = LpInteger)
y = LpVariable.dicts('y', S, lowBound = 0, upBound = 1, cat = LpInteger)
```

In [11]:

```
prob += lpSum([lpSum([x[i][j]*d[i][j] for j in S]) for i in D])
```

In [12]:

```
for i in D:
prob += lpSum([x[i][j] for j in S]) == 1
prob += lpSum([y[j] for j in S]) == n
for i in D:
for j in S:
prob += x[i][j] <= y[j]
```

Note that in the 2nd constrained we use n to define the number of facilities to be opened.

In [13]:

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

In [14]:

```
chosen_facil = []
for j,facility in enumerate(S):
if y[facility].varValue == 1:
chosen_facil.append(j)
print(f'We will be opening facility {facility}')
chosen_facil
```

We will be opening facility S2 We will be opening facility S4

Out[14]:

[1, 3]

```
from pulp import *
from scipy.spatial import distance
# Calculate distance matrix
dist_matrix = distance.cdist(cust_coord, facil_coord, 'euclidean')
# Initialise problem
prob = LpProblem('prob', LpMinimize)
# Define sets of customers and facilities
n_customers = len(cust_coord)
n_facil = len(facil_coord)
S = [('S'+str(i+1)) for i in range(num_facility)]
D = [('C'+str(i+1)) for i in range(num_customer)]
# Prepare distance matrix parameters
d = {customer:{facility:dist_matrix[i][j] for j,facility in enumerate(S)} for i,customer in enumerate(D)}
# Define decision variables
x = LpVariable.dicts('x', (D,S), lowBound = 0, upBound = 1, cat = LpInteger)
y = LpVariable.dicts('y', S, lowBound = 0, upBound = 1, cat = LpInteger)
# Objective function
prob += lpSum([lpSum([x[i][j]*d[i][j] for j in S]) for i in D])
# Constraint 1: Each customer will be served by 1 facility
prob += lpSum([y[j] for j in S]) == 1
# Constraint 2: Only 1 facility will be established
prob += lpSum([y[j] for j in S]) == 1
# Constraint 3: Only established facilities can serve customers
for i in D:
for j in S:
prob += x[i][j] <= y[j]
# Solve the model
status = prob.solve()
```

In [15]:

```
plt.scatter(cust_coord[:, 0],
cust_coord[:, 1],
s=plot_size*2,
cmap='viridis',
zorder = 10000);
plt.scatter(facil_coord[:, 0],
facil_coord[:, 1],
s=plot_size*5,
cmap='viridis');
plt.scatter(facil_coord[chosen_facil, 0],
facil_coord[chosen_facil, 1],
s=plot_size*10,
cmap='viridis',
zorder=99999)
```

Out[15]:

<matplotlib.collections.PathCollection at 0x13e46fd90>

In [16]:

```
vect_orig = []
vect_dest = []
u = v = np.zeros((len(vect_orig),len(vect_orig)))
for j,facility in enumerate(S):
for i,customer in enumerate(D):
if x[customer][facility].varValue == 1:
vect_orig.append(facil_coord[j])
vect_dest.append(cust_coord[i])
vect_orig = np.array(vect_orig)
vect_dest = np.array(vect_dest)
```

In the above code we are iterating through our sets of facilities and customers.

Wherever we find a pair that is connected, we will append the origin and destination coordinates of a link between them to `vect_orig`

and `vect_dest`

, respectively.

In [17]:

```
plt.scatter(cust_coord[:, 0],
cust_coord[:, 1],
s=plot_size*2,
cmap='viridis',
zorder = 10000);
plt.scatter(facil_coord[:, 0],
facil_coord[:, 1],
s=plot_size*5,
cmap='viridis');
plt.scatter(facil_coord[chosen_facil, 0],
facil_coord[chosen_facil, 1],
s=plot_size*10,
cmap='viridis',
zorder=99999)
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')
```

We have now replotted the graph, but with the addition of arrows between facilities and the customers that they serve.

The eagled-eyed among you might have notices that we are now also using the `zorder`

parameter, which allows us to define the order with which the elements of the graph are served.