Transport Analytics Training Series - Last Revision: October 2022
Transport Analytics Training Series - Last Revision: October 2022
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: