Benchmarking the MTZ Model

Santa Claus tried to use the MTZ model in order to plan a route for all the children in London. But soon after he clicked the run button, his server farm caught fire!

Should we perhaps check how well the MTZ model performs for different problem sizes?

In this notebook we will see how the solution time requirements of the TSP/MTZ formulation vary as we begin to increase the size of the problem.

To do this, we will use a for-loop that will measure how long it takes to run the model for a range of problem sizes.

As the solution times might vary for different problem instances of the same size, we will be generating 5 different instances, by varying the value of our blob generator's random_state parameter.

As we will be solving the model many times, we are going to define a function that encapsulates the entire process of instance generation, model formulation and execution.

The resulting time_solver functions returns the run-time, in seconds.

import time
import gc
from sklearn.datasets import make_blobs
from scipy.spatial import distance
from pulp import *

def time_solver(num_cities, rand_seed):

    # We are saving the time before we start the operation.
    start = time.time()

    # As before, we are generating our cities. Every time we 
    # are using a different combination of city numbers and
    # random seeds
    cities_coord, _ = make_blobs(n_samples=num_cities, 
                           centers=2, 
                           cluster_std=20, 
                           center_box=center_box, 
                           random_state = rand_seed)
    
    # Again, calculating our distance matrix...
    dist_matrix = distance.cdist(cities_coord, cities_coord, 'euclidean')

    # Formulating our model
    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
    u = LpVariable.dicts('u', N, 0, num_cities-1, LpInteger)
    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

    # Solving the problem...
    status = prob.solve()

    # Storing the time that we obtained a solution....
    end = time.time()
    
    # Asking the Garbage Collector to free up memory 
    # (don't worry about this!)
    gc.collect()
    
    # And ultimately returning the 
    return end - start

The following code implements our nested for-loop that actually calls the solver. As you can see if you run it, the code returns the run-time for every single instance of the same size, and then an average value for all 5.

center_box = (100, 200) 

print(f'Starting Iterations\n')

for n_cities in range(5,40,2):
    print("")
    sum_=0
    for iter_ in range(5):
        elapsed = time_solver(n_cities,iter_)
        print(f'#{iter_+1}: Nodes={n_cities} --> {elapsed}')
        sum_+= elapsed
    avg=sum_/5
    print(f'#AVG for {n_cities}: {round(avg,2)}')

Starting Iterations


#1: Nodes=5 --> 0.04161190986633301
#2: Nodes=5 --> 0.02799510955810547
#3: Nodes=5 --> 0.034867048263549805
#4: Nodes=5 --> 0.04032421112060547
#5: Nodes=5 --> 0.03560805320739746
#AVG for 5: 0.04

#1: Nodes=7 --> 0.13392996788024902
#2: Nodes=7 --> 0.11475205421447754
#3: Nodes=7 --> 0.22600769996643066
#4: Nodes=7 --> 0.3792400360107422
#5: Nodes=7 --> 0.1381242275238037
#AVG for 7: 0.2

#1: Nodes=9 --> 0.5028700828552246
#2: Nodes=9 --> 0.16100406646728516
#3: Nodes=9 --> 0.3553199768066406
#4: Nodes=9 --> 0.3422539234161377
#5: Nodes=9 --> 0.2715470790863037
#AVG for 9: 0.33

#1: Nodes=11 --> 0.785477876663208
#2: Nodes=11 --> 0.38356900215148926
#3: Nodes=11 --> 1.4543190002441406
#4: Nodes=11 --> 0.6511647701263428
#5: Nodes=11 --> 0.5775840282440186
#AVG for 11: 0.77

#1: Nodes=13 --> 0.45705103874206543
#2: Nodes=13 --> 0.7822351455688477
#3: Nodes=13 --> 5.205451011657715
#4: Nodes=13 --> 2.501317262649536
#5: Nodes=13 --> 1.008293867111206
#AVG for 13: 1.99

#1: Nodes=15 --> 1.3093912601470947
#2: Nodes=15 --> 2.3266632556915283
#3: Nodes=15 --> 2.6898627281188965
#4: Nodes=15 --> 19.195463180541992
#5: Nodes=15 --> 11.021780014038086
#AVG for 15: 7.31

#1: Nodes=17 --> 7.777148962020874
#2: Nodes=17 --> 7.428498983383179
#3: Nodes=17 --> 3.4429140090942383
#4: Nodes=17 --> 5.9074249267578125
#5: Nodes=17 --> 169.41908311843872
#AVG for 17: 38.8

#1: Nodes=19 --> 4.574414253234863
#2: Nodes=19 --> 5.637892246246338
#3: Nodes=19 --> 41.972410917282104
#4: Nodes=19 --> 8.959167242050171
#5: Nodes=19 --> 30.69031310081482
#AVG for 19: 18.37

#1: Nodes=21 --> 1.181546926498413
#2: Nodes=21 --> 14.072715997695923
#3: Nodes=21 --> 23.084228992462158
#4: Nodes=21 --> 143.51868081092834
#5: Nodes=21 --> 361.980357170105
#AVG for 21: 108.77

#1: Nodes=23 --> 333.8613290786743
#2: Nodes=23 --> 12.9081289768219
#3: Nodes=23 --> 28.929744005203247
#4: Nodes=23 --> 628.4628422260284
#5: Nodes=23 --> 42.0759961605072
#AVG for 23: 209.25

#1: Nodes=25 --> 2.7486460208892822
#2: Nodes=25 --> 6.412477254867554
#3: Nodes=25 --> 66.76630091667175
#4: Nodes=25 --> 1315.6705391407013
#5: Nodes=25 --> 139.22304797172546
#AVG for 25: 306.16

#1: Nodes=27 --> 91.64272999763489
#2: Nodes=27 --> 49.908385038375854
#3: Nodes=27 --> 2531.7875142097473
#4: Nodes=27 --> 42.82020926475525
#5: Nodes=27 --> 23.829869031906128
#AVG for 27: 548.0

#1: Nodes=29 --> 3.3185040950775146
#2: Nodes=29 --> 996.8490288257599
#3: Nodes=29 --> 11.784288883209229
#4: Nodes=29 --> 50.01945614814758
#5: Nodes=29 --> 97.68729782104492
#AVG for 29: 231.93

#1: Nodes=31 --> 52.46259784698486
#2: Nodes=31 --> 76.54280209541321
#3: Nodes=31 --> 153.19842791557312
#4: Nodes=31 --> 63.102813959121704
#5: Nodes=31 --> 463.2422709465027
#AVG for 31: 161.71

#1: Nodes=33 --> 170.7737557888031
#2: Nodes=33 --> 80.32510685920715
#3: Nodes=33 --> 55.51656484603882
#4: Nodes=33 --> 67.96415996551514

---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
<ipython-input-3-66c518aa5686> in <module>
      7     sum_=0
      8     for iter_ in range(5):
----> 9         elapsed = time_solver(n_cities,iter_)
     10         print(f'#{iter_+1}: Nodes={n_cities} --> {elapsed}')
     11         sum_+= elapsed

<ipython-input-2-ef1af55bc5ee> in time_solver(num_cities, rand_seed)
     34 
     35     # Solving the problem...
---> 36     status = prob.solve()
     37 
     38     # Storing the time that we obtained a solution....

~/opt/anaconda3/lib/python3.8/site-packages/pulp/pulp.py in solve(self, solver, **kwargs)
   1897         #time it
   1898         self.solutionTime = -clock()
-> 1899         status = solver.actualSolve(self, **kwargs)
   1900         self.solutionTime += clock()
   1901         self.restoreObjective(wasNone, dummyVar)

~/opt/anaconda3/lib/python3.8/site-packages/pulp/apis/coin_api.py in actualSolve(self, lp, **kwargs)
     99     def actualSolve(self, lp, **kwargs):
    100         """Solve a well formulated lp problem"""
--> 101         return self.solve_CBC(lp, **kwargs)
    102 
    103     def available(self):

~/opt/anaconda3/lib/python3.8/site-packages/pulp/apis/coin_api.py in solve_CBC(self, lp, use_mps)
    151         args.extend(cmds[1:].split())
    152         cbc = subprocess.Popen(args, stdout = pipe, stderr = pipe, stdin=devnull)
--> 153         if cbc.wait() != 0:
    154             raise PulpSolverError("Pulp: Error while trying to execute, use msg=True for more details" +  \
    155                                     self.path)

~/opt/anaconda3/lib/python3.8/subprocess.py in wait(self, timeout)
   1077             endtime = _time() + timeout
   1078         try:
-> 1079             return self._wait(timeout=timeout)
   1080         except KeyboardInterrupt:
   1081             # https://bugs.python.org/issue25942

~/opt/anaconda3/lib/python3.8/subprocess.py in _wait(self, timeout)
   1802                         if self.returncode is not None:
   1803                             break  # Another thread waited.
-> 1804                         (pid, sts) = self._try_wait(0)
   1805                         # Check the pid and loop as waitpid has been known to
   1806                         # return 0 even without WNOHANG in odd situations.

~/opt/anaconda3/lib/python3.8/subprocess.py in _try_wait(self, wait_flags)
   1760             """All callers to this function MUST hold self._waitpid_lock."""
   1761             try:
-> 1762                 (pid, sts) = os.waitpid(self.pid, wait_flags)
   1763             except ChildProcessError:
   1764                 # This happens if SIGCLD is set to be ignored or waiting

KeyboardInterrupt: