This is a simplified version of problem faced by logistics company.

Problem: Given there are # dispatchers and # of locations for delivery, optimize the paths for each dispatcher.

Solution:

Linear Programming. Computing adjacency matrices for all dispatcher in a single run is very expensive and slow. One could adopt divide and conquer strategy for scalability. For example, during the first iteration, set number of dispatcher to 1 to get overall shortest path, then cut it into 2 and rerun the optimizer for each of the parts. More constraints can be introduced to set the minimum and maximum work to be carried out by each of the dispatcher.

```
from pulp import *
# X to be solved. X is the adjacency matrices. #num_dispatcher X adjacency matrix
adj_mats = LpVariable.dicts("adj_mats", (range_num_dispatcher,
range_num_loc,
range_num_loc),0,1,LpInteger)
prob = LpProblem("Logistics Problem", LpMinimize)
for i in range_num_loc:
for k in range_num_dispatcher:
# Ensuring no self-recursion.
prob += choices[k][i][i] == 0, ""
# Ensuring one route will only be taken up by one dispatcher.
for j in range_num_loc[int(i):]:
prob += lpSum([adj_mats[k][i][j]
for k in range_num_dispatcher if i != j] +
[adj_mats[k][j][i]
for k in range_num_dispatcher if i != j]) <= 1, ""
if int(i) > 0:
# Ensuring every location is covered.
prob += lpSum([adj_mats[k][i][j] for j in range_num_loc
for k in range_num_dispatcher if i != j]) == 1, ""
prob += lpSum([adj_mats[k][j][i] for j in range_num_loc
for k in range_num_dispatcher if i != j]) == 1, ""
for k in range_num_dispatcher:
# Ensuring every dispatcher will go through the distribution centre.
prob += lpSum([adj_mats[k]["0"][j] for j in range_num_loc[1:]]) == 1, ""
prob += lpSum([adj_mats[k][i]["0"] for i in range_num_loc[1:]]) == 1, ""
for i in range_num_loc[1:]:
# Ensuring exit of location if entered.
prob += lpSum([adj_mats[k][i][j] for j in range_num_loc] + [adj_mats[k][j][i] for j in range_num_loc]) == lpSum([adj_mats[k][i][j] for j in range_num_loc]) * 2 , ""
# Objective function - Minimize travelling.
prob += lpSum([adj_mats[k][i][j] * euc_mat[int(i), int(j)] for k in range_num_dispatcher for i in range_num_loc for j in range_num_loc])
prob.solve()
```