TSP Essay

A fun study of some heuristics for the Travelling Salesman Problem.

---

I've always found the Travelling salesman problem particularly interesting. After studying a few papers on the subject (in this repository too), mainly the superb Clarke & Wright's Savings Algorithm from 1997, I felt like trying a few techniques.

To see the original ipynb format (probably better): https://github.com/rsalmei/tsp-essay/blob/master/tsp_essay.ipynb

To run this yourself, you need:

• a recent python (I'm using 3.6.5)
• a virtualenv
• pip install -r requirements.txt
• jupyter notebook tsp_essay.ipynb

Enjoy :)

Setup environment

import pandas as pd
import numpy as np

%matplotlib inline
%config InlineBackend.figure_format = 'retina'
import matplotlib.pyplot as plt
from pylab import rcParams

/Users/rogerio/.pyenv/versions/3.6.5/lib/python3.6/importlib/_bootstrap.py:219: RuntimeWarning: numpy.dtype size changed, may indicate binary incompatibility. Expected 96, got 88
  return f(*args, **kwds)
/Users/rogerio/.pyenv/versions/3.6.5/lib/python3.6/importlib/_bootstrap.py:219: RuntimeWarning: numpy.dtype size changed, may indicate binary incompatibility. Expected 96, got 88
  return f(*args, **kwds)

Customize graphics style

print(plt.style.available)

['seaborn-dark', 'seaborn-darkgrid', 'seaborn-ticks', 'fivethirtyeight', 'seaborn-whitegrid', 'classic', '_classic_test', 'fast', 'seaborn-talk', 'seaborn-dark-palette', 'seaborn-bright', 'seaborn-pastel', 'grayscale', 'seaborn-notebook', 'ggplot', 'seaborn-colorblind', 'seaborn-muted', 'seaborn', 'Solarize_Light2', 'seaborn-paper', 'bmh', 'tableau-colorblind10', 'seaborn-white', 'dark_background', 'seaborn-poster', 'seaborn-deep']

plt.style.use('classic')
rcParams['figure.figsize'] = 16, 10

Initialize data

np.random.seed(18)

size, cmin, cmax = 150, -100, 100
data = pd.DataFrame(dict(x=[0], y=[0])).append(
    pd.DataFrame((np.random.random_sample(2*size)*(cmax-cmin)+cmin).reshape(-1,2), columns=['x', 'y']), ignore_index=True)
data

<style scoped> .dataframe tbody tr th:only-of-type { vertical-align: middle; }

.dataframe tbody tr th {
    vertical-align: top;
}

.dataframe thead th {
    text-align: right;
}

</style>

	x	y
0	0.000000	0.000000
1	30.074848	1.090675
2	75.720294	-63.631955
3	70.446614	50.027257
4	33.220333	97.579090
5	-48.606315	-94.338815
6	27.143823	69.462477
7	47.234925	-95.838578
8	-77.679374	-40.455252
9	37.394038	72.325211
10	-60.273128	31.437806
11	39.931126	-29.521506
12	57.992110	62.809588
13	-60.515743	89.891059
14	-4.679923	33.320070
15	-57.737593	61.505265
16	-36.062776	-41.212805
17	-39.616364	65.577880
18	10.506738	-61.664519
19	42.851348	27.974071
20	6.640766	14.505583
21	-97.533493	39.047596
22	6.931842	-84.173057
23	-97.272500	-44.429392
24	99.477056	-54.879872
25	-78.996929	-96.739844
26	-31.674510	61.639612
27	-56.536075	86.920994
28	97.699575	-97.796661
29	56.980551	37.072724
...	...	...
121	89.841848	-55.198411
122	15.997366	-84.731802
123	57.065474	-1.286114
124	8.397404	45.994062
125	-95.971359	91.328450
126	33.324722	-82.136313
127	34.783261	52.531436
128	45.139227	-60.250218
129	35.987597	83.219655
130	67.314264	13.003851
131	-36.886862	-35.034144
132	42.068799	65.334579
133	-69.769883	-43.590844
134	33.143474	-27.783041
135	-41.797934	75.064175
136	5.324715	37.861594
137	68.512768	-50.739106
138	-64.144841	-50.464696
139	-26.845766	79.781956
140	0.068680	88.957681
141	18.895939	-56.449630
142	-26.104390	40.905810
143	-78.601095	8.400365
144	28.941010	83.297786
145	-61.078038	33.677658
146	44.080113	-41.024718
147	-61.088577	2.871590
148	15.453097	21.393731
149	26.546590	68.767646
150	-63.012893	11.913804

151 rows × 2 columns

Plot problem

plt.scatter(data.x, data.y, c='r', marker='o', s=40)
plt.scatter(data.iloc[0].x, data.iloc[0].y, marker='*', c='r', s=400)
plt.axis([-100,100,-100,100])
plt.grid(True)

Compute distance matrix

distance_func = lambda a,b: np.sqrt((a.x-b.x)**2 + (a.y-b.y)**2)
def compute_distances(dfunc=distance_func):
    for i in range(size+1):
        current=data.iloc[i]
        data[i]=dfunc(current, data)
        #data[str(i)][i]=np.nan

%time compute_distances()
data

CPU times: user 346 ms, sys: 5.93 ms, total: 352 ms
Wall time: 357 ms

<style scoped> .dataframe tbody tr th:only-of-type { vertical-align: middle; }

.dataframe tbody tr th {
    vertical-align: top;
}

.dataframe thead th {
    text-align: right;
}

</style>

	x	y	0	1	2	3	4	5	6	7	...	141	142	143	144	145	146	147	148	149	150
0	0.000000	0.000000	0.000000	30.094619	98.906970	86.402846	103.078947	106.124389	74.577630	106.846484	...	59.528289	48.525504	79.048708	88.182216	69.747483	60.216973	61.156032	26.391096	73.713707	64.129271
1	30.074848	1.090675	30.094619	0.000000	79.199277	63.440275	96.539672	123.683115	68.434599	98.436519	...	58.616164	68.857475	108.921496	82.214930	96.802687	44.383034	91.180819	25.020186	67.768880	93.714821
2	75.720294	-63.631955	98.906970	79.199277	0.000000	113.781493	166.719068	128.062552	141.682044	42.996311	...	57.276463	145.932898	170.304863	154.196790	167.877766	38.886864	152.116365	104.218532	141.236354	157.968538
3	70.446614	50.027257	86.402846	63.440275	113.781493	0.000000	60.390171	187.123389	47.464297	147.701132	...	118.299618	96.980912	154.751471	53.194391	132.536951	94.792692	139.732470	62.001337	47.732738	138.795084
4	33.220333	97.579090	103.078947	96.539672	166.719068	60.390171	0.000000	208.633848	28.765741	193.924735	...	154.693357	82.044400	143.027539	14.908663	113.910385	139.028595	133.655082	78.229684	29.574281	128.838562
5	-48.606315	-94.338815	106.124389	123.683115	128.062552	187.123389	208.633848	0.000000	180.468687	95.852974	...	77.408944	137.103775	107.028154	193.825565	128.622553	106.925989	98.008518	132.278609	179.587519	107.224850
6	27.143823	69.462477	74.577630	68.434599	141.682044	47.464297	28.765741	180.468687	0.000000	166.517540	...	126.181957	60.422310	122.108842	13.951546	95.203204	111.777718	110.540955	49.469965	0.916230	106.958325
7	47.234925	-95.838578	106.846484	98.436519	42.996311	147.701132	193.924735	95.852974	166.517540	0.000000	...	48.524090	155.169851	163.402757	180.068054	168.837654	54.904573	146.552646	121.463982	165.901224	154.159518
8	-77.679374	-40.455252	87.582604	115.486087	155.140639	173.575323	177.065594	61.226474	151.887491	136.641483	...	97.890812	96.330687	48.864310	163.348464	75.969024	121.760819	46.394717	111.798720	150.972491	54.384039
9	37.394038	72.325211	81.420208	71.609564	141.255984	39.870685	25.596449	187.544551	10.642469	168.451485	...	130.096654	70.846519	132.443410	13.851032	105.784607	113.546950	120.509879	55.456475	11.415928	117.179733
10	-60.273128	31.437806	67.979303	95.308474	165.929113	132.034914	114.523783	126.316558	95.328899	166.605100	...	118.287542	35.456251	29.438717	103.192151	2.380087	127.044938	28.577853	76.389427	94.504923	19.715298
11	39.931126	-29.521506	49.658978	32.159786	49.440745	85.200944	127.277634	109.727672	99.806534	66.718059	...	34.170205	96.543753	124.450615	113.353321	119.151104	12.228572	106.086252	56.493675	99.196287	110.970068
12	57.992110	62.809588	85.487596	67.739189	127.678315	17.846642	42.691441	189.891666	31.557530	159.012444	...	125.504070	86.902225	147.030837	35.549018	122.582092	104.762143	133.314565	59.370370	32.004987	131.272965
13	-60.515743	89.891059	108.363082	126.854891	205.254899	136.895075	94.050826	184.614411	90.008480	214.722389	...	166.498686	59.863977	83.473428	89.699398	56.216213	167.568594	87.021354	102.289533	89.588215	78.017229
14	-4.679923	33.320070	33.647120	47.398608	125.951935	76.961852	74.603291	135.004884	48.156250	139.201679	...	92.813901	22.727764	78.008534	60.234037	56.399249	88.908315	64.101842	23.400342	47.240087	62.136653
15	-57.737593	61.505265	84.359512	106.587745	182.948985	128.697068	97.850218	156.111362	85.253574	189.146277	...	140.662914	37.749133	57.056254	89.376138	28.027384	144.496514	58.729353	83.461441	84.596487	49.871252
16	-36.062776	-41.212805	54.763301	78.509679	114.009092	140.246209	155.123626	54.586750	127.452314	99.611656	...	57.031756	82.720230	65.352699	140.457750	78.957867	80.143110	50.692451	81.076899	126.552884	59.571347
17	-39.616364	65.577880	76.615369	94.949802	173.198517	111.156111	79.556659	160.169188	66.873109	183.298716	...	135.330717	28.129780	69.203161	70.810371	38.447726	135.533075	66.280727	70.603715	66.239800	58.542554
18	10.506738	-61.664519	62.553213	65.735267	65.243228	126.758990	160.855321	67.542303	132.178220	50.167978	...	9.877943	108.908434	113.354726	146.129710	119.224623	39.410315	96.388788	83.205405	131.414712	104.013969
19	42.851348	27.974071	51.174082	29.765012	97.324363	35.324803	70.268165	152.725068	44.362306	123.890225	...	87.756612	70.157849	123.019616	57.045691	104.085773	69.009729	106.928212	28.177386	43.931320	107.075533
20	6.640766	14.505583	15.953423	27.002147	104.295043	73.027225	87.222021	122.062865	58.656932	117.574315	...	72.005774	42.062065	85.460217	72.316444	70.380430	66.972524	68.721275	11.184979	57.798040	69.701862
21	-97.533493	39.047596	105.059492	133.133829	201.395046	168.338555	143.256758	142.076751	128.333542	197.869084	...	150.583972	71.453269	36.023443	133.992086	36.848832	162.683708	51.351098	114.357458	127.589765	43.908023
22	6.931842	-84.173057	84.458001	88.348756	71.789888	148.471716	183.643479	56.460868	154.959354	41.957394	...	30.194834	129.368139	126.038572	168.910884	136.066655	56.936572	110.469669	105.910143	154.193375	118.848436
23	-97.272500	-44.429392	106.938815	135.238395	174.055293	192.488336	192.859504	69.709024	168.673588	153.379595	...	116.788665	111.117040	56.032174	179.566372	86.085715	141.393610	59.553834	130.536369	167.763931	65.941462
24	99.477056	-54.879872	113.611114	89.159232	25.317636	108.849770	166.233837	153.250426	143.851022	66.384153	...	80.596406	157.941751	188.987345	155.139934	183.358596	57.103298	170.635735	113.479903	143.553342	175.682608
25	-78.996929	-96.739844	124.896406	146.517791	158.219947	209.461105	224.393766	30.485314	197.203121	126.235072	...	105.859884	147.458288	105.140954	209.914619	131.642742	135.100458	101.208433	151.249294	196.296113	109.823061
26	-31.674510	61.639612	69.301633	86.482119	165.004877	102.779232	74.182119	156.894728	59.336275	176.142212	...	128.461820	21.468972	70.968457	64.368608	40.576327	127.588116	65.718093	61.973721	58.655820	58.777123
27	-56.536075	86.920994	103.689859	121.935621	200.394455	132.233699	90.386990	181.433181	85.481724	210.165367	...	162.003470	55.167786	81.561966	85.553841	53.436712	162.768923	84.172605	97.346099	85.042773	75.286305
28	97.699575	-97.796661	138.236731	119.799035	40.624081	150.315118	205.740751	146.346747	181.531632	50.502624	...	88.992079	185.918792	205.814807	193.708387	206.145160	78.090334	188.010037	144.813096	181.125415	194.589012
29	56.980551	37.072724	67.979188	44.929107	102.433443	18.685684	65.004371	168.575143	44.037785	133.268118	...	100.979544	83.173313	138.580254	54.064519	118.107395	79.155743	122.922888	44.388740	43.940802	122.602601
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
121	89.841848	-55.198411	105.443930	82.100886	16.448190	106.998208	162.932381	143.874477	139.539885	58.881007	...	70.956942	150.597316	180.049519	151.294774	175.144991	47.906471	161.716165	106.770982	139.190062	166.938975
122	15.997366	-84.731802	86.228731	86.969380	63.340601	145.343471	183.122614	65.314090	154.596635	33.153364	...	28.430318	132.504217	132.749649	168.527388	141.284883	51.951419	116.690174	106.126929	153.861518	124.831864
123	57.065474	-1.286114	57.079966	27.095074	65.076925	53.029397	101.700144	140.802458	76.815809	95.062127	...	67.081494	93.259771	136.011933	89.137095	123.208583	41.806413	118.227181	47.391617	76.412903	120.801706
124	8.397404	45.994062	46.754360	49.862068	128.647717	62.180150	57.246772	151.468612	30.036557	147.053904	...	102.980237	34.874979	94.773546	42.586471	70.558705	94.050645	81.779271	25.592168	29.120939	79.125814
125	-95.971359	91.328450	132.481650	155.017749	231.280675	171.466411	129.342815	191.613623	125.041868	235.668275	...	187.170662	86.161684	84.727760	125.170250	67.388112	192.696071	95.086406	131.553313	124.577837	85.982245
126	33.324722	-82.136313	88.639219	83.290415	46.257927	137.277982	179.715433	82.834751	151.724741	19.525516	...	29.461763	136.642539	143.959295	165.492170	149.414715	42.495197	127.044144	105.061245	151.056109	134.634173
127	34.783261	52.531436	63.003389	51.655794	123.165636	35.751163	45.074759	168.892546	18.574746	148.891587	...	110.133009	61.987590	121.669896	31.316133	97.697766	94.016941	107.969947	36.649855	18.205969	105.895608
128	45.139227	-60.250218	75.283720	63.163602	30.767480	113.144091	158.278711	99.750986	130.955022	35.650012	...	26.517063	123.726290	141.508196	144.459031	141.790525	19.254651	123.566617	86.873475	130.350663	130.017410
129	35.987597	83.219655	90.667624	82.341544	152.131793	47.845157	14.623649	196.680299	16.354579	179.411128	...	140.711171	75.139046	136.852089	7.047020	108.977736	124.507642	126.014265	65.146838	17.262478	122.006644
130	67.314264	13.003851	68.558809	39.098566	77.095448	37.155675	91.188634	157.987432	69.290986	110.679059	...	84.664752	97.496483	145.987959	80.085853	130.046105	58.812516	128.801988	52.535424	69.076799	130.331715
131	-36.886862	-35.034144	50.872702	76.084644	116.181782	136.952243	150.004295	60.451548	122.553957	103.796215	...	59.752355	76.701618	60.221539	135.409580	72.845897	81.188287	44.972966	76.964773	121.649555	53.727826
132	42.068799	65.334579	77.707085	65.353914	133.284620	32.243053	33.436564	183.623444	15.485298	161.255931	...	123.969250	72.417874	133.426864	22.248948	107.895462	106.378312	120.594648	51.373083	15.897325	117.881047
133	-69.769883	-43.590844	82.267845	109.386510	146.864011	168.597197	174.745343	54.984118	148.907085	128.140356	...	89.593395	95.112353	52.735908	160.762448	77.755833	113.878913	47.266508	107.172530	147.990855	55.914424
134	33.143474	-27.783041	43.247974	29.036321	55.659054	86.290015	125.362154	105.416788	97.430420	69.499101	...	32.011960	90.710902	117.456747	111.160293	112.494937	17.174169	99.092814	52.261883	96.775793	104.028296
135	-41.797934	75.064175	85.916806	103.139593	181.788752	115.002981	78.324082	169.539751	69.168959	192.703401	...	144.843457	37.590972	76.148121	71.216505	45.657050	144.401073	74.725486	78.474180	68.633958	66.618645
136	5.324715	37.861594	38.234185	44.324593	123.517116	66.248509	65.911644	142.777814	38.401684	140.114958	...	95.282659	31.576193	88.946645	51.207196	66.534434	87.892156	75.066808	19.333251	37.490693	73.097992
137	68.512768	-50.739106	85.255242	64.527513	14.770714	100.784918	152.459316	124.971254	127.121243	49.866912	...	49.944368	131.723944	158.555876	139.756261	154.660813	26.293040	140.251971	89.545938	126.661060	145.685918
138	-64.144841	-50.464696	81.616458	107.402542	140.483568	167.968724	177.191816	46.544432	150.718768	120.267375	...	83.256175	98.972952	60.614178	162.964343	84.198224	108.635878	53.423778	107.235560	149.804162	62.388769
139	-26.845766	79.781956	84.177525	97.119896	176.316042	101.740598	62.647221	175.475253	54.966966	190.605668	...	143.705771	38.883215	88.169982	55.897455	57.423469	140.088304	84.188922	72.099782	54.516591	76.903493
140	0.068680	88.957681	88.957708	92.849238	170.313722	80.427803	34.254355	189.649310	33.363548	190.720508	...	146.621116	54.717565	112.598467	29.421860	82.430589	137.231302	105.598414	69.293345	33.297406	99.574313
141	18.895939	-56.449630	59.528289	58.616164	57.276463	118.299618	154.693357	77.408944	126.181957	48.524090	...	0.000000	107.252558	117.094805	140.107971	120.493838	29.532534	99.581774	77.919459	125.450782	106.689342
142	-26.104390	40.905810	48.525504	68.857475	145.932898	96.980912	82.044400	137.103775	60.422310	155.169851	...	107.252558	0.000000	61.745510	69.477160	35.712774	107.881768	51.676834	45.910194	59.568512	46.933719
143	-78.601095	8.400365	79.048708	108.921496	170.304863	154.751471	143.027539	107.028154	122.108842	163.402757	...	117.094805	61.745510	0.000000	131.053150	30.757098	132.263062	18.364522	94.947451	121.244564	15.979246
144	28.941010	83.297786	88.182216	82.214930	154.196790	53.194391	14.908663	193.825565	13.951546	180.068054	...	140.107971	69.477160	131.053150	0.000000	102.789037	125.240877	120.721578	63.356419	14.726106	116.409592
145	-61.078038	33.677658	69.747483	96.802687	167.877766	132.536951	113.910385	128.622553	95.203204	168.837654	...	120.493838	35.712774	30.757098	102.789037	0.000000	128.991014	30.806070	77.510705	94.389526	21.849691
146	44.080113	-41.024718	60.216973	44.383034	38.886864	94.792692	139.028595	106.925989	111.777718	54.904573	...	29.532534	107.881768	132.263062	125.240877	128.991014	0.000000	113.962007	68.670000	111.183576	119.462961
147	-61.088577	2.871590	61.156032	91.180819	152.116365	139.732470	133.655082	98.008518	110.540955	146.552646	...	99.581774	51.676834	18.364522	120.721578	30.806070	113.962007	0.000000	78.750857	109.645851	9.244708
148	15.453097	21.393731	26.391096	25.020186	104.218532	62.001337	78.229684	132.278609	49.469965	121.463982	...	77.919459	45.910194	94.947451	63.356419	77.510705	68.670000	78.750857	0.000000	48.655456	79.036577
149	26.546590	68.767646	73.713707	67.768880	141.236354	47.732738	29.574281	179.587519	0.916230	165.901224	...	125.450782	59.568512	121.244564	14.726106	94.389526	111.183576	109.645851	48.655456	0.000000	106.081385
150	-63.012893	11.913804	64.129271	93.714821	157.968538	138.795084	128.838562	107.224850	106.958325	154.159518	...	106.689342	46.933719	15.979246	116.409592	21.849691	119.462961	9.244708	79.036577	106.081385	0.000000

151 rows × 153 columns

Helper functions

estimate_cost = lambda route: sum(data.iloc[i][j] for i,j in zip(route, route[1:]))
get_coords = lambda route: list(zip(*[(data.iloc[i].x,data.iloc[i].y) for i in route]))

Scenario: center depot, one van, nearest neighbor heuristic

Apply a greedy NN algorithm

def nearest_neighbor(unserved, stop_condition=lambda route, target: False):
    current=0 #depot
    result_path=[]
    while True:
        result_path.append(current)
        unserved.remove(current)
        if not unserved:
            break

        current=data.iloc[unserved,current+2].idxmin()
        if stop_condition(result_path, int(current)):
            if len(result_path)>1:
                break

    result_path.append(0)
    return result_path

%time route = nearest_neighbor(list(range(size+1)))
print('cost={}\nroute={}'.format(estimate_cost(route),route))

CPU times: user 116 ms, sys: 8.29 ms, total: 125 ms
Wall time: 122 ms
cost=2623.828598616268
route=[0, 53, 60, 20, 66, 148, 55, 51, 1, 120, 104, 123, 74, 93, 130, 42, 29, 19, 91, 107, 81, 11, 134, 117, 57, 30, 34, 67, 75, 46, 39, 131, 16, 114, 72, 96, 38, 22, 122, 92, 36, 126, 113, 7, 58, 111, 28, 108, 63, 2, 137, 50, 128, 69, 141, 18, 146, 37, 54, 95, 24, 121, 35, 83, 3, 49, 12, 73, 105, 129, 90, 144, 9, 86, 132, 31, 70, 127, 149, 6, 71, 41, 140, 76, 102, 32, 100, 97, 139, 88, 78, 26, 17, 135, 27, 13, 110, 89, 59, 119, 40, 15, 56, 145, 10, 106, 82, 94, 115, 87, 98, 44, 64, 52, 14, 33, 136, 124, 65, 142, 150, 147, 143, 109, 112, 23, 47, 101, 45, 138, 99, 133, 8, 48, 62, 80, 103, 43, 68, 61, 125, 21, 25, 118, 5, 79, 116, 77, 85, 84, 4, 0]

Plot result path

coords=get_coords(route)
plt.plot(coords[0], coords[1], 'b-', linewidth=2)
plt.scatter(data.x, data.y, c='r', s=30)
plt.scatter(data.iloc[0].x, data.iloc[0].y, marker='*', c='r', s=400)
plt.axis([-100,100,-100,100])
plt.grid(True)

Let's improve it with a 2-opt heuristic

def twoOptSwap(route,i,j):
    t=route[i:j+1];
    t.reverse();
    route[i:j+1]=t

def twoOpt(route):
    cost=estimate_cost(route)

    def inner():
        nonlocal cost
        for i in range(1,len(route)-2):
            nodei=data.iloc[route[i]]
            for j in range(i+1,len(route)-1):
                nodej=data.iloc[route[j]]
                save=nodei.iloc[route[i-1]+2]+nodej.iloc[route[j+1]+2] - (nodej.iloc[route[i-1]+2]+nodei.iloc[route[j+1]+2])
                if save>0:
                    twoOptSwap(route,i,j)

                    cost-=save
                    print('exchanging {}-{},{}-{} with {}-{},{}-{} => save={} => cost: {}'.format(
                            route[i-1],route[i],route[j],route[j+1],route[i-1],route[j],route[i],route[j+1],save,cost))
                    return False
        return True

    while True:
        if inner():
            break

%time twoOpt(route)

exchanging 53-64,60-52 with 53-60,64-52 => save=2.8803537373675177 => cost: 2620.948244878901
...
exchanging 1-120,81-107 with 1-81,120-107 => save=26.496627104257605 => cost: 2100.7350743873462
CPU times: user 3min 11s, sys: 1.65 s, total: 3min 13s
Wall time: 3min 21s

Plot improved path

coords=get_coords(route)
plt.plot(coords[0], coords[1], 'b-', linewidth=2)
plt.scatter(data.x, data.y, c='r', s=30)
plt.scatter(data.iloc[0].x, data.iloc[0].y, marker='*', c='r', s=400)
plt.axis([-100,100,-100,100])
plt.grid(True)

Incredible! A clean path, without any crossing parts!

Scenario: center depot, one van, 2-opt heuristic

Just for fun, let's pass a random route to 2-opt heuristic

route = [x for x in range(size+1)]+[0]
print('cost={}\nroute={}'.format(estimate_cost(route),route))

cost=16332.039263088578
route=[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 145, 146, 147, 148, 149, 150, 0]

Plot result path

coords=get_coords(route)
plt.plot(coords[0], coords[1], 'b-', linewidth=2)
plt.scatter(data.x, data.y, c='r', s=30)
plt.scatter(data.iloc[0].x, data.iloc[0].y, marker='*', c='r', s=400)
plt.axis([-100,100,-100,100])
plt.grid(True)

Let's see what 2-opt heuristic can do

%time twoOpt(route)

exchanging 0-4,1-5 with 0-1,4-5 => save=11.966404491447804 => cost: 16320.07285859713
...
exchanging 47-23,112-8 with 47-112,23-8 => save=11.556726776627453 => cost: 2181.5435075388827
CPU times: user 19min 11s, sys: 8.18 s, total: 19min 19s
Wall time: 19min 44s

Plot improved path

coords=get_coords(route)
plt.plot(coords[0], coords[1], 'b-', linewidth=2)
plt.scatter(data.x, data.y, c='r', s=30)
plt.scatter(data.iloc[0].x, data.iloc[0].y, marker='*', c='r', s=400)
plt.axis([-100,100,-100,100])
plt.grid(True)

WOW, very cool, a clean path too, almost as good as the previous one.

But extremely slower... It took ~19 minutes, while NN + 2-opt took only ~3 minutes.

Scenario: center depot, several vans, no clustering

Let's apply a sequential NN algorithm, but including a max cost constraint

def sequential_nn(max_cost):
    unserved=list(range(size+1))
    result=[]
    while True:
        result_path = nearest_neighbor(unserved, lambda route, target: estimate_cost(route+[target,0])>=max_cost)
        result.append(result_path)
        print('van {}: {}'.format(len(result), estimate_cost(result_path)))

        if not unserved:
            break
        unserved.append(0)
    return result

%time result_paths = sequential_nn(400)

van 1: 389.0742611200774
van 2: 398.3533734081084
van 3: 397.4195886051539
van 4: 390.9917686566429
van 5: 371.6762623788105
van 6: 320.29683749922486
van 7: 364.8394134778564
van 8: 320.99556479051967
van 9: 210.1189842618746
van 10: 276.59421133431306
CPU times: user 1.48 s, sys: 24.4 ms, total: 1.51 s
Wall time: 1.54 s

Plot result paths

for path in result_paths:
    coords=get_coords(path)
    plt.plot(coords[0], coords[1], linewidth=2)
plt.scatter(data.x, data.y, c='r', s=30)
plt.scatter(data.iloc[0].x, data.iloc[0].y, marker='*', c='r', s=400)
plt.axis([-100,100,-100,100])
plt.grid(True)

Let's improve them with a 2-opt heuristic

def improve_seq_nn():
    for i,path in enumerate(result_paths):
        print('Starting path {}/{}'.format(i+1, len(result_paths)))
        twoOpt(path)

%time improve_seq_nn()

Starting path 1/10
exchanging 55-19,51-91 with 55-51,19-91 => save=2.6873147586143844 => cost: 386.386946361463
...
exchanging 67-39,75-0 with 67-75,39-0 => save=4.275672113835263 => cost: 368.50842741951476
Starting path 2/10
exchanging 17-97,78-0 with 17-78,97-0 => save=0.09782392307953103 => cost: 398.25554948502884
...
exchanging 88-78,26-0 with 88-26,78-0 => save=1.6805548383893267 => cost: 394.249851279636
Starting path 3/10
Starting path 4/10
exchanging 124-127,65-0 with 124-65,127-0 => save=0.5258127779938491 => cost: 390.465955878649
exchanging 14-65,33-0 with 14-33,65-0 => save=3.343148155406432 => cost: 387.1228077232426
Starting path 5/10
exchanging 0-18,141-69 with 0-141,18-69 => save=6.7520080457503155 => cost: 364.9242543330602
exchanging 50-121,146-0 with 50-146,121-0 => save=41.908681513869396 => cost: 323.0155728191908
Starting path 6/10
exchanging 47-45,101-138 with 47-101,45-138 => save=1.4275147235708587 => cost: 318.869322775654
Starting path 7/10
exchanging 3-35,83-77 with 3-83,35-77 => save=3.633104522812438 => cost: 361.20630895504394
...
exchanging 90-129,144-0 with 90-144,129-0 => save=5.489882024385494 => cost: 347.690672350374
Starting path 8/10
Starting path 9/10
Starting path 10/10
CPU times: user 3.25 s, sys: 37.7 ms, total: 3.29 s
Wall time: 3.31 s

Plot improved paths

for path in result_paths:
    coords=get_coords(path)
    plt.plot(coords[0], coords[1], linewidth=2)
plt.scatter(data.x, data.y, c='r', s=30)
plt.scatter(data.iloc[0].x, data.iloc[0].y, marker='*', c='r', s=400)
plt.axis([-100,100,-100,100])
plt.grid(True)

Great! Several clean paths, but with many superpositions, let's improve them once more.

Scenario: center depot, several vans, KMeans clustering

Compute delivery clusters

from sklearn.cluster import KMeans

work_per_van=20 #how many deliveries each van will start, ie granularity; a better approach would be capacity.
model = KMeans(n_clusters=1+size//work_per_van, init='k-means++', random_state=0)
%time model.fit(data[['x', 'y']])

CPU times: user 75.7 ms, sys: 4.05 ms, total: 79.7 ms
Wall time: 92.8 ms





KMeans(algorithm='auto', copy_x=True, init='k-means++', max_iter=300,
    n_clusters=8, n_init=10, n_jobs=1, precompute_distances='auto',
    random_state=0, tol=0.0001, verbose=0)

Plot clusters

plt.scatter(data.x, data.y, marker='o', c=model.labels_, s=200)
plt.scatter(model.cluster_centers_[:,0], model.cluster_centers_[:,1], marker='+', c='r', s=500)
plt.scatter(data.iloc[0].x, data.iloc[0].y, marker='*', c='r', s=400)
plt.axis([-100,100,-100,100])
plt.grid(True)

Let's apply a greedy NN algorithm sequentially, on each of the clusters

def greedy_nn_clusters():
    result=[]
    for i in range(model.n_clusters):
        unserved=data[model.labels_ == i].index.tolist()
        if 0 not in unserved:
            unserved.append(0)
        result_path = nearest_neighbor(unserved)
        result.append(result_path)
        print('van {}: {}'.format(i+1, estimate_cost(result_path)))
    return result

%time result_paths = greedy_nn_clusters()

van 1: 409.3080162039939
van 2: 266.4745473994821
van 3: 413.3606064604762
van 4: 446.80922078430706
van 5: 457.66477785115126
van 6: 455.80304195030294
van 7: 459.76894476196577
van 8: 376.5008949566861
CPU times: user 164 ms, sys: 10.1 ms, total: 174 ms
Wall time: 174 ms

Plot result paths

plt.scatter(data.x, data.y, marker='o', c=model.labels_, s=200)
plt.scatter(model.cluster_centers_[:,0], model.cluster_centers_[:,1], marker='+', c='r', s=500)
for path in result_paths:
    coords=get_coords(path)
    plt.plot(coords[0], coords[1], 'k-', linewidth=2)
plt.scatter(data.iloc[0].x, data.iloc[0].y, marker='*', c='r', s=400)
plt.axis([-100,100,-100,100])
plt.grid(True)

Let's improve them with a 2-opt heuristic

def improve_result_paths():
    for i, path in enumerate(result_paths):
        print('Starting path {}/{}'.format(i+1, len(result_paths)))
        twoOpt(path)

%time improve_result_paths()

Starting path 1/8
exchanging 48-112,99-138 with 48-99,112-138 => save=11.5162692633929 => cost: 397.791746940601
exchanging 45-118,25-5 with 45-25,118-5 => save=13.984640981159586 => cost: 383.8071059594414
Starting path 2/8
exchanging 148-52,136-65 with 148-136,52-65 => save=1.6057804518926915 => cost: 264.8687669475894
exchanging 148-124,52-142 with 148-52,124-142 => save=6.904413311483367 => cost: 257.964353636106
Starting path 3/8
exchanging 63-24,108-111 with 63-108,24-111 => save=11.657922849361384 => cost: 401.7026836111148
...
exchanging 113-36,126-0 with 113-126,36-0 => save=7.772149456496223 => cost: 378.1902685773841
Starting path 4/8
exchanging 147-109,150-0 with 147-150,109-0 => save=5.601665728075147 => cost: 441.2075550562319
...
exchanging 62-109,143-150 with 62-143,109-150 => save=0.5147852919035643 => cost: 394.93412566866925
Starting path 5/8
exchanging 9-149,6-144 with 9-6,149-144 => save=0.0011015802705607314 => cost: 457.66367627088067
...
exchanging 105-129,144-90 with 105-144,129-90 => save=3.918490824545831 => cost: 407.2051777421843
Starting path 6/8
exchanging 30-69,34-116 with 30-34,69-116 => save=1.1688372549307786 => cost: 454.63420469537215
...
exchanging 22-116,79-38 with 22-79,116-38 => save=5.039989423215701 => cost: 380.39750649818404
Starting path 7/8
exchanging 0-55,1-107 with 0-1,55-107 => save=14.962441964752642 => cost: 444.8065027972131
...
exchanging 0-107,55-42 with 0-55,107-42 => save=6.044195398723886 => cost: 403.197665386089
Starting path 8/8
exchanging 0-71,26-59 with 0-26,71-59 => save=17.738760993960298 => cost: 358.76213396272584
...
exchanging 17-78,26-88 with 17-26,78-88 => save=0.17155006280330554 => cost: 319.53217682612376
CPU times: user 2.81 s, sys: 64.8 ms, total: 2.87 s
Wall time: 3.03 s

Plot improved paths

plt.scatter(data.x, data.y, marker='o', c=model.labels_, s=200)
plt.scatter(model.cluster_centers_[:,0], model.cluster_centers_[:,1], marker='+', c='r', s=500)
for path in result_paths:
    coords=get_coords(path)
    plt.plot(coords[0], coords[1], 'k-', linewidth=2)
plt.scatter(data.iloc[0].x, data.iloc[0].y, marker='*', c='r', s=400)
plt.axis([-100,100,-100,100])
plt.grid(True)

Cool clean paths, with clearly defined cluster subsections!

Let's change the depot to the farthest point

Find farthest point and set it as new depot

maxid=data[0].idxmax()
print('farthest point={}; distance={}; coords={}'.format(maxid, data.iloc[maxid,2], tuple(data.iloc[maxid][['x','y']].tolist())))
data.drop(data.columns.difference(['x', 'y']), axis=1, inplace=True)
data.iloc[0]=data.iloc[maxid]
data.drop(maxid, inplace=True)
data.reset_index(drop=True, inplace=True)
size-=1
data

farthest point=84; distance=138.29710566715653; coords=(99.07656719022526, 96.48794364952255)

<style scoped> .dataframe tbody tr th:only-of-type { vertical-align: middle; }

.dataframe tbody tr th {
    vertical-align: top;
}

.dataframe thead th {
    text-align: right;
}

</style>

	x	y
0	99.076567	96.487944
1	30.074848	1.090675
2	75.720294	-63.631955
3	70.446614	50.027257
4	33.220333	97.579090
5	-48.606315	-94.338815
6	27.143823	69.462477
7	47.234925	-95.838578
8	-77.679374	-40.455252
9	37.394038	72.325211
10	-60.273128	31.437806
11	39.931126	-29.521506
12	57.992110	62.809588
13	-60.515743	89.891059
14	-4.679923	33.320070
15	-57.737593	61.505265
16	-36.062776	-41.212805
17	-39.616364	65.577880
18	10.506738	-61.664519
19	42.851348	27.974071
20	6.640766	14.505583
21	-97.533493	39.047596
22	6.931842	-84.173057
23	-97.272500	-44.429392
24	99.477056	-54.879872
25	-78.996929	-96.739844
26	-31.674510	61.639612
27	-56.536075	86.920994
28	97.699575	-97.796661
29	56.980551	37.072724
...	...	...
120	89.841848	-55.198411
121	15.997366	-84.731802
122	57.065474	-1.286114
123	8.397404	45.994062
124	-95.971359	91.328450
125	33.324722	-82.136313
126	34.783261	52.531436
127	45.139227	-60.250218
128	35.987597	83.219655
129	67.314264	13.003851
130	-36.886862	-35.034144
131	42.068799	65.334579
132	-69.769883	-43.590844
133	33.143474	-27.783041
134	-41.797934	75.064175
135	5.324715	37.861594
136	68.512768	-50.739106
137	-64.144841	-50.464696
138	-26.845766	79.781956
139	0.068680	88.957681
140	18.895939	-56.449630
141	-26.104390	40.905810
142	-78.601095	8.400365
143	28.941010	83.297786
144	-61.078038	33.677658
145	44.080113	-41.024718
146	-61.088577	2.871590
147	15.453097	21.393731
148	26.546590	68.767646
149	-63.012893	11.913804

150 rows × 2 columns

Plot relocated depot

plt.scatter(data.x, data.y, c='r', marker='o', s=40)
plt.scatter(data.iloc[0].x, data.iloc[0].y, marker='*', c='r', s=400)
plt.axis([-100,100,-100,100])
plt.grid(True)

Recompute distance matrix

%time compute_distances()
data

CPU times: user 365 ms, sys: 5.95 ms, total: 371 ms
Wall time: 374 ms

<style scoped> .dataframe tbody tr th:only-of-type { vertical-align: middle; }

.dataframe tbody tr th {
    vertical-align: top;
}

.dataframe thead th {
    text-align: right;
}

</style>

	x	y	0	1	2	3	4	5	6	7	...	140	141	142	143	144	145	146	147	148	149
0	99.076567	96.487944	0.000000	117.736469	161.814392	54.573525	65.865272	241.298748	76.842016	199.190980	...	172.681309	136.965856	198.314834	71.365094	172.030897	148.102471	185.517910	112.392284	77.646716	182.827181
1	30.074848	1.090675	117.736469	0.000000	79.199277	63.440275	96.539672	123.683115	68.434599	98.436519	...	58.616164	68.857475	108.921496	82.214930	96.802687	44.383034	91.180819	25.020186	67.768880	93.714821
2	75.720294	-63.631955	161.814392	79.199277	0.000000	113.781493	166.719068	128.062552	141.682044	42.996311	...	57.276463	145.932898	170.304863	154.196790	167.877766	38.886864	152.116365	104.218532	141.236354	157.968538
3	70.446614	50.027257	54.573525	63.440275	113.781493	0.000000	60.390171	187.123389	47.464297	147.701132	...	118.299618	96.980912	154.751471	53.194391	132.536951	94.792692	139.732470	62.001337	47.732738	138.795084
4	33.220333	97.579090	65.865272	96.539672	166.719068	60.390171	0.000000	208.633848	28.765741	193.924735	...	154.693357	82.044400	143.027539	14.908663	113.910385	139.028595	133.655082	78.229684	29.574281	128.838562
5	-48.606315	-94.338815	241.298748	123.683115	128.062552	187.123389	208.633848	0.000000	180.468687	95.852974	...	77.408944	137.103775	107.028154	193.825565	128.622553	106.925989	98.008518	132.278609	179.587519	107.224850
6	27.143823	69.462477	76.842016	68.434599	141.682044	47.464297	28.765741	180.468687	0.000000	166.517540	...	126.181957	60.422310	122.108842	13.951546	95.203204	111.777718	110.540955	49.469965	0.916230	106.958325
7	47.234925	-95.838578	199.190980	98.436519	42.996311	147.701132	193.924735	95.852974	166.517540	0.000000	...	48.524090	155.169851	163.402757	180.068054	168.837654	54.904573	146.552646	121.463982	165.901224	154.159518
8	-77.679374	-40.455252	223.598080	115.486087	155.140639	173.575323	177.065594	61.226474	151.887491	136.641483	...	97.890812	96.330687	48.864310	163.348464	75.969024	121.760819	46.394717	111.798720	150.972491	54.384039
9	37.394038	72.325211	66.246298	71.609564	141.255984	39.870685	25.596449	187.544551	10.642469	168.451485	...	130.096654	70.846519	132.443410	13.851032	105.784607	113.546950	120.509879	55.456475	11.415928	117.179733
10	-60.273128	31.437806	172.115792	95.308474	165.929113	132.034914	114.523783	126.316558	95.328899	166.605100	...	118.287542	35.456251	29.438717	103.192151	2.380087	127.044938	28.577853	76.389427	94.504923	19.715298
11	39.931126	-29.521506	139.199730	32.159786	49.440745	85.200944	127.277634	109.727672	99.806534	66.718059	...	34.170205	96.543753	124.450615	113.353321	119.151104	12.228572	106.086252	56.493675	99.196287	110.970068
12	57.992110	62.809588	53.124047	67.739189	127.678315	17.846642	42.691441	189.891666	31.557530	159.012444	...	125.504070	86.902225	147.030837	35.549018	122.582092	104.762143	133.314565	59.370370	32.004987	131.272965
13	-60.515743	89.891059	159.728596	126.854891	205.254899	136.895075	94.050826	184.614411	90.008480	214.722389	...	166.498686	59.863977	83.473428	89.699398	56.216213	167.568594	87.021354	102.289533	89.588215	78.017229
14	-4.679923	33.320070	121.472588	47.398608	125.951935	76.961852	74.603291	135.004884	48.156250	139.201679	...	92.813901	22.727764	78.008534	60.234037	56.399249	88.908315	64.101842	23.400342	47.240087	62.136653
15	-57.737593	61.505265	160.668816	106.587745	182.948985	128.697068	97.850218	156.111362	85.253574	189.146277	...	140.662914	37.749133	57.056254	89.376138	28.027384	144.496514	58.729353	83.461441	84.596487	49.871252
16	-36.062776	-41.212805	192.935581	78.509679	114.009092	140.246209	155.123626	54.586750	127.452314	99.611656	...	57.031756	82.720230	65.352699	140.457750	78.957867	80.143110	50.692451	81.076899	126.552884	59.571347
17	-39.616364	65.577880	142.095606	94.949802	173.198517	111.156111	79.556659	160.169188	66.873109	183.298716	...	135.330717	28.129780	69.203161	70.810371	38.447726	135.533075	66.280727	70.603715	66.239800	58.542554
18	10.506738	-61.664519	181.264492	65.735267	65.243228	126.758990	160.855321	67.542303	132.178220	50.167978	...	9.877943	108.908434	113.354726	146.129710	119.224623	39.410315	96.388788	83.205405	131.414712	104.013969
19	42.851348	27.974071	88.630841	29.765012	97.324363	35.324803	70.268165	152.725068	44.362306	123.890225	...	87.756612	70.157849	123.019616	57.045691	104.085773	69.009729	106.928212	28.177386	43.931320	107.075533
20	6.640766	14.505583	123.553571	27.002147	104.295043	73.027225	87.222021	122.062865	58.656932	117.574315	...	72.005774	42.062065	85.460217	72.316444	70.380430	66.972524	68.721275	11.184979	57.798040	69.701862
21	-97.533493	39.047596	204.828975	133.133829	201.395046	168.338555	143.256758	142.076751	128.333542	197.869084	...	150.583972	71.453269	36.023443	133.992086	36.848832	162.683708	51.351098	114.357458	127.589765	43.908023
22	6.931842	-84.173057	202.802977	88.348756	71.789888	148.471716	183.643479	56.460868	154.959354	41.957394	...	30.194834	129.368139	126.038572	168.910884	136.066655	56.936572	110.469669	105.910143	154.193375	118.848436
23	-97.272500	-44.429392	241.682957	135.238395	174.055293	192.488336	192.859504	69.709024	168.673588	153.379595	...	116.788665	111.117040	56.032174	179.566372	86.085715	141.393610	59.553834	130.536369	167.763931	65.941462
24	99.477056	-54.879872	151.368346	89.159232	25.317636	108.849770	166.233837	153.250426	143.851022	66.384153	...	80.596406	157.941751	188.987345	155.139934	183.358596	57.103298	170.635735	113.479903	143.553342	175.682608
25	-78.996929	-96.739844	262.768240	146.517791	158.219947	209.461105	224.393766	30.485314	197.203121	126.235072	...	105.859884	147.458288	105.140954	209.914619	131.642742	135.100458	101.208433	151.249294	196.296113	109.823061
26	-31.674510	61.639612	135.315373	86.482119	165.004877	102.779232	74.182119	156.894728	59.336275	176.142212	...	128.461820	21.468972	70.968457	64.368608	40.576327	127.588116	65.718093	61.973721	58.655820	58.777123
27	-56.536075	86.920994	155.906450	121.935621	200.394455	132.233699	90.386990	181.433181	85.481724	210.165367	...	162.003470	55.167786	81.561966	85.553841	53.436712	162.768923	84.172605	97.346099	85.042773	75.286305
28	97.699575	-97.796661	194.289484	119.799035	40.624081	150.315118	205.740751	146.346747	181.531632	50.502624	...	88.992079	185.918792	205.814807	193.708387	206.145160	78.090334	188.010037	144.813096	181.125415	194.589012
29	56.980551	37.072724	72.816502	44.929107	102.433443	18.685684	65.004371	168.575143	44.037785	133.268118	...	100.979544	83.173313	138.580254	54.064519	118.107395	79.155743	122.922888	44.388740	43.940802	122.602601
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
120	89.841848	-55.198411	151.967201	82.100886	16.448190	106.998208	162.932381	143.874477	139.539885	58.881007	...	70.956942	150.597316	180.049519	151.294774	175.144991	47.906471	161.716165	106.770982	139.190062	166.938975
121	15.997366	-84.731802	199.355838	86.969380	63.340601	145.343471	183.122614	65.314090	154.596635	33.153364	...	28.430318	132.504217	132.749649	168.527388	141.284883	51.951419	116.690174	106.126929	153.861518	124.831864
122	57.065474	-1.286114	106.417566	27.095074	65.076925	53.029397	101.700144	140.802458	76.815809	95.062127	...	67.081494	93.259771	136.011933	89.137095	123.208583	41.806413	118.227181	47.391617	76.412903	120.801706
123	8.397404	45.994062	103.789897	49.862068	128.647717	62.180150	57.246772	151.468612	30.036557	147.053904	...	102.980237	34.874979	94.773546	42.586471	70.558705	94.050645	81.779271	25.592168	29.120939	79.125814
124	-95.971359	91.328450	195.116154	155.017749	231.280675	171.466411	129.342815	191.613623	125.041868	235.668275	...	187.170662	86.161684	84.727760	125.170250	67.388112	192.696071	95.086406	131.553313	124.577837	85.982245
125	33.324722	-82.136313	190.341615	83.290415	46.257927	137.277982	179.715433	82.834751	151.724741	19.525516	...	29.461763	136.642539	143.959295	165.492170	149.414715	42.495197	127.044144	105.061245	151.056109	134.634173
126	34.783261	52.531436	77.883271	51.655794	123.165636	35.751163	45.074759	168.892546	18.574746	148.891587	...	110.133009	61.987590	121.669896	31.316133	97.697766	94.016941	107.969947	36.649855	18.205969	105.895608
127	45.139227	-60.250218	165.759126	63.163602	30.767480	113.144091	158.278711	99.750986	130.955022	35.650012	...	26.517063	123.726290	141.508196	144.459031	141.790525	19.254651	123.566617	86.873475	130.350663	130.017410
128	35.987597	83.219655	64.469107	82.341544	152.131793	47.845157	14.623649	196.680299	16.354579	179.411128	...	140.711171	75.139046	136.852089	7.047020	108.977736	124.507642	126.014265	65.146838	17.262478	122.006644
129	67.314264	13.003851	89.322100	39.098566	77.095448	37.155675	91.188634	157.987432	69.290986	110.679059	...	84.664752	97.496483	145.987959	80.085853	130.046105	58.812516	128.801988	52.535424	69.076799	130.331715
130	-36.886862	-35.034144	189.166893	76.084644	116.181782	136.952243	150.004295	60.451548	122.553957	103.796215	...	59.752355	76.701618	60.221539	135.409580	72.845897	81.188287	44.972966	76.964773	121.649555	53.727826
131	42.068799	65.334579	64.964742	65.353914	133.284620	32.243053	33.436564	183.623444	15.485298	161.255931	...	123.969250	72.417874	133.426864	22.248948	107.895462	106.378312	120.594648	51.373083	15.897325	117.881047
132	-69.769883	-43.590844	219.388219	109.386510	146.864011	168.597197	174.745343	54.984118	148.907085	128.140356	...	89.593395	95.112353	52.735908	160.762448	77.755833	113.878913	47.266508	107.172530	147.990855	55.914424
133	33.143474	-27.783041	140.678536	29.036321	55.659054	86.290015	125.362154	105.416788	97.430420	69.499101	...	32.011960	90.710902	117.456747	111.160293	112.494937	17.174169	99.092814	52.261883	96.775793	104.028296
134	-41.797934	75.064175	142.494221	103.139593	181.788752	115.002981	78.324082	169.539751	69.168959	192.703401	...	144.843457	37.590972	76.148121	71.216505	45.657050	144.401073	74.725486	78.474180	68.633958	66.618645
135	5.324715	37.861594	110.573318	44.324593	123.517116	66.248509	65.911644	142.777814	38.401684	140.114958	...	95.282659	31.576193	88.946645	51.207196	66.534434	87.892156	75.066808	19.333251	37.490693	73.097992
136	68.512768	-50.739106	150.366053	64.527513	14.770714	100.784918	152.459316	124.971254	127.121243	49.866912	...	49.944368	131.723944	158.555876	139.756261	154.660813	26.293040	140.251971	89.545938	126.661060	145.685918
137	-64.144841	-50.464696	219.627654	107.402542	140.483568	167.968724	177.191816	46.544432	150.718768	120.267375	...	83.256175	98.972952	60.614178	162.964343	84.198224	108.635878	53.423778	107.235560	149.804162	62.388769
138	-26.845766	79.781956	127.025683	97.119896	176.316042	101.740598	62.647221	175.475253	54.966966	190.605668	...	143.705771	38.883215	88.169982	55.897455	57.423469	140.088304	84.188922	72.099782	54.516591	76.903493
139	0.068680	88.957681	99.293840	92.849238	170.313722	80.427803	34.254355	189.649310	33.363548	190.720508	...	146.621116	54.717565	112.598467	29.421860	82.430589	137.231302	105.598414	69.293345	33.297406	99.574313
140	18.895939	-56.449630	172.681309	58.616164	57.276463	118.299618	154.693357	77.408944	126.181957	48.524090	...	0.000000	107.252558	117.094805	140.107971	120.493838	29.532534	99.581774	77.919459	125.450782	106.689342
141	-26.104390	40.905810	136.965856	68.857475	145.932898	96.980912	82.044400	137.103775	60.422310	155.169851	...	107.252558	0.000000	61.745510	69.477160	35.712774	107.881768	51.676834	45.910194	59.568512	46.933719
142	-78.601095	8.400365	198.314834	108.921496	170.304863	154.751471	143.027539	107.028154	122.108842	163.402757	...	117.094805	61.745510	0.000000	131.053150	30.757098	132.263062	18.364522	94.947451	121.244564	15.979246
143	28.941010	83.297786	71.365094	82.214930	154.196790	53.194391	14.908663	193.825565	13.951546	180.068054	...	140.107971	69.477160	131.053150	0.000000	102.789037	125.240877	120.721578	63.356419	14.726106	116.409592
144	-61.078038	33.677658	172.030897	96.802687	167.877766	132.536951	113.910385	128.622553	95.203204	168.837654	...	120.493838	35.712774	30.757098	102.789037	0.000000	128.991014	30.806070	77.510705	94.389526	21.849691
145	44.080113	-41.024718	148.102471	44.383034	38.886864	94.792692	139.028595	106.925989	111.777718	54.904573	...	29.532534	107.881768	132.263062	125.240877	128.991014	0.000000	113.962007	68.670000	111.183576	119.462961
146	-61.088577	2.871590	185.517910	91.180819	152.116365	139.732470	133.655082	98.008518	110.540955	146.552646	...	99.581774	51.676834	18.364522	120.721578	30.806070	113.962007	0.000000	78.750857	109.645851	9.244708
147	15.453097	21.393731	112.392284	25.020186	104.218532	62.001337	78.229684	132.278609	49.469965	121.463982	...	77.919459	45.910194	94.947451	63.356419	77.510705	68.670000	78.750857	0.000000	48.655456	79.036577
148	26.546590	68.767646	77.646716	67.768880	141.236354	47.732738	29.574281	179.587519	0.916230	165.901224	...	125.450782	59.568512	121.244564	14.726106	94.389526	111.183576	109.645851	48.655456	0.000000	106.081385
149	-63.012893	11.913804	182.827181	93.714821	157.968538	138.795084	128.838562	107.224850	106.958325	154.159518	...	106.689342	46.933719	15.979246	116.409592	21.849691	119.462961	9.244708	79.036577	106.081385	0.000000

150 rows × 152 columns

Scenario: farthest depot, one van, nearest neighbor heuristic

Apply a greedy NN algorithm

%time route = nearest_neighbor(list(range(size+1)))
print('cost={}\nroute={}'.format(estimate_cost(route),route))

CPU times: user 111 ms, sys: 13.4 ms, total: 124 ms
Wall time: 118 ms
cost=2598.9001149206997
route=[0, 84, 77, 83, 3, 49, 12, 73, 104, 128, 89, 143, 9, 85, 131, 31, 70, 126, 148, 6, 71, 41, 139, 76, 101, 32, 99, 96, 138, 87, 78, 26, 17, 134, 27, 13, 109, 88, 59, 118, 40, 15, 56, 144, 10, 105, 82, 93, 114, 86, 97, 44, 64, 53, 60, 20, 66, 147, 55, 51, 1, 119, 103, 122, 74, 92, 129, 42, 29, 19, 90, 106, 81, 11, 133, 116, 57, 30, 34, 67, 75, 46, 39, 130, 16, 113, 72, 95, 38, 22, 121, 91, 36, 125, 112, 7, 58, 110, 28, 107, 63, 2, 136, 50, 127, 69, 140, 18, 145, 37, 54, 94, 24, 120, 35, 123, 135, 33, 14, 52, 65, 141, 149, 146, 142, 108, 111, 23, 47, 100, 45, 137, 98, 132, 8, 48, 62, 80, 102, 43, 68, 61, 124, 21, 25, 117, 5, 79, 115, 4, 0]

Plot result path

coords=get_coords(route)
plt.plot(coords[0], coords[1], 'b-', linewidth=2)
plt.scatter(data.x, data.y, c='r', s=30)
plt.scatter(data.iloc[0].x, data.iloc[0].y, marker='*', c='r', s=400)
plt.axis([-100,100,-100,100])
plt.grid(True)

Let's improve it with a 2-opt heuristic

%time twoOpt(route)

exchanging 84-4,77-0 with 84-77,4-0 => save=15.967192548480725 => cost: 2582.932922372219
...
exchanging 3-49,12-83 with 3-12,49-83 => save=1.0654890882007315 => cost: 2111.9918767877866
CPU times: user 3min 36s, sys: 1.67 s, total: 3min 38s
Wall time: 3min 44s

Plot improved path

coords=get_coords(route)
plt.plot(coords[0], coords[1], 'b-', linewidth=2)
plt.scatter(data.x, data.y, c='r', s=30)
plt.scatter(data.iloc[0].x, data.iloc[0].y, marker='*', c='r', s=400)
plt.axis([-100,100,-100,100])
plt.grid(True)

A great optimized route too.

Scenario: farthest depot, one van, 2-opt heuristic

Just for fun, let's pass a random route to 2-opt heuristic

route = [x for x in range(size+1)]+[0]
print('cost={}\nroute={}'.format(estimate_cost(route),route))

cost=16502.514122712284
route=[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 145, 146, 147, 148, 149, 0]

Plot result path

coords=get_coords(route)
plt.plot(coords[0], coords[1], 'b-', linewidth=2)
plt.scatter(data.x, data.y, c='r', s=30)
plt.scatter(data.iloc[0].x, data.iloc[0].y, marker='*', c='r', s=400)
plt.axis([-100,100,-100,100])
plt.grid(True)

Let's see what 2-opt heuristic can do

%time twoOpt(route)

exchanging 0-2,1-3 with 0-1,2-3 => save=6.263294946633721 => cost: 16496.25082776565
...
exchanging 31-70,126-131 with 31-126,70-131 => save=1.7915970025161734 => cost: 2147.6785119636334
CPU times: user 30min 1s, sys: 27.7 s, total: 30min 29s
Wall time: 33min 24s

Plot improved path

coords=get_coords(route)
plt.plot(coords[0], coords[1], 'b-', linewidth=2)
plt.scatter(data.x, data.y, c='r', s=30)
plt.scatter(data.iloc[0].x, data.iloc[0].y, marker='*', c='r', s=400)
plt.axis([-100,100,-100,100])
plt.grid(True)

Yes, it's beautiful what this heuristic can do.

Only again extremely slower. It's best suited to improve a path, not to build one.

Scenario: farthest depot, several vans, no clustering

Let's apply a sequential NN algorithm, but including a max cost constraint

%time result_paths = sequential_nn(400)

van 1: 385.64822798111754
van 2: 396.65132151762236
van 3: 356.93203827425646
van 4: 396.0161561312648
van 5: 387.061496278307
van 6: 395.6147269472501
van 7: 369.199156243418
van 8: 394.6734767469922
van 9: 397.84112119751114
van 10: 398.459895479128
van 11: 398.88707002367426
van 12: 393.1275886242274
van 13: 399.0252974546198
van 14: 369.5354695202354
van 15: 399.04095929693733
van 16: 398.1246462968952
van 17: 390.2323088432243
van 18: 391.6922432861831
van 19: 396.62966713208937
van 20: 397.37088575736055
van 21: 398.69882895151926
van 22: 398.71167524672734
van 23: 401.9265755043099
van 24: 405.5614801848182
van 25: 405.6059533951813
van 26: 406.6405588064364
van 27: 409.6579507334966
van 28: 416.22815658152126
van 29: 426.0013563494297
van 30: 427.5654816187484
van 31: 428.5781354316992
van 32: 438.7764370279704
van 33: 439.2553078706148
van 34: 445.7400163235422
van 35: 447.1961602342642
van 36: 461.09599384348525
van 37: 464.4374453992292
van 38: 471.7901899743616
van 39: 482.2725996546565
van 40: 482.59749534967307
van 41: 483.3659132930797
van 42: 494.26350457825447
van 43: 525.5364809069648
van 44: 551.3280289580272
CPU times: user 1.22 s, sys: 48.8 ms, total: 1.27 s
Wall time: 1.43 s

Plot result paths

for path in result_paths:
    coords=get_coords(path)
    plt.plot(coords[0], coords[1], linewidth=2)
plt.scatter(data.x, data.y, c='r', s=30)
plt.scatter(data.iloc[0].x, data.iloc[0].y, marker='*', c='r', s=400)
plt.axis([-100,100,-100,100])
plt.grid(True)

Max cost 400 doesn't quite work here, as the distances are much larger

%time result_paths = sequential_nn(500)

van 1: 493.1613509811791
van 2: 498.37402590912467
van 3: 499.57415652250484
van 4: 490.0564767221547
van 5: 497.67366460125766
van 6: 491.0474410501304
van 7: 476.1857519522164
van 8: 487.87440569372217
van 9: 462.17325815341474
van 10: 486.2976818275437
van 11: 498.3091236408933
van 12: 466.4397450832517
van 13: 464.4374453992292
van 14: 494.26350457825447
van 15: 525.5364809069648
van 16: 551.3280289580272
CPU times: user 573 ms, sys: 20.1 ms, total: 593 ms
Wall time: 588 ms

Plot result paths

for path in result_paths:
    coords=get_coords(path)
    plt.plot(coords[0], coords[1], linewidth=2)
plt.scatter(data.x, data.y, c='r', s=30)
plt.scatter(data.iloc[0].x, data.iloc[0].y, marker='*', c='r', s=400)
plt.axis([-100,100,-100,100])
plt.grid(True)

Let's improve them with a 2-opt heuristic

%time improve_result_paths()

Starting path 1/16
exchanging 84-134,77-0 with 84-77,134-0 => save=17.58671573136749 => cost: 475.57463524981165
...
exchanging 131-9,85-73 with 131-85,9-73 => save=0.22871091978264957 => cost: 415.7579280346732
Starting path 2/16
exchanging 114-10,93-0 with 114-93,10-0 => save=4.063910634805211 => cost: 494.31011527431946
exchanging 114-105,10-82 with 114-10,105-82 => save=2.900319415888454 => cost: 491.409795858431
Starting path 3/16
exchanging 4-13,123-0 with 4-123,13-0 => save=19.13464498164811 => cost: 480.43951154085676
exchanging 65-14,52-33 with 65-52,14-33 => save=1.587339787656024 => cost: 478.8521717532007
Starting path 4/16
exchanging 19-16,55-0 with 19-55,16-0 => save=0.7537655368855951 => cost: 489.3027111852691
...
exchanging 46-130,39-75 with 46-39,130-75 => save=1.7704562419552907 => cost: 455.4046302803684
Starting path 5/16
exchanging 37-112,94-0 with 37-94,112-0 => save=11.94699370814547 => cost: 485.7266708931122
exchanging 58-110,28-107 with 58-28,110-107 => save=9.304516959883209 => cost: 476.422153933229
Starting path 6/16
exchanging 62-142,144-0 with 62-144,142-0 => save=23.25445433373693 => cost: 467.79298671639344
Starting path 7/16
Starting path 8/16
exchanging 44-100,48-0 with 44-48,100-0 => save=1.0746379367560621 => cost: 486.7997677569661
...
exchanging 44-98,8-48 with 44-8,98-48 => save=6.17952602155286 => cost: 476.7280985588525
Starting path 9/16
Starting path 10/16
exchanging 113-5,72-0 with 113-72,5-0 => save=1.4836706701319144 => cost: 484.8140111574118
Starting path 11/16
Starting path 12/16
Starting path 13/16
Starting path 14/16
Starting path 15/16
Starting path 16/16
CPU times: user 2.16 s, sys: 40 ms, total: 2.2 s
Wall time: 2.24 s

Plot improved paths

for path in result_paths:
    coords=get_coords(path)
    plt.plot(coords[0], coords[1], linewidth=2)
plt.scatter(data.x, data.y, c='r', s=30)
plt.scatter(data.iloc[0].x, data.iloc[0].y, marker='*', c='r', s=400)
plt.axis([-100,100,-100,100])
plt.grid(True)

Nice! Several clean paths, but some vans had to do only one distant delivery.

As we see, it is better to increase the max cost again.

In a real system, there should be a way to estimate it based on a full day of work (generally 8 hours). Let's try again.

%time result_paths = sequential_nn(800)

van 1: 792.9220639985251
van 2: 793.9102149907774
van 3: 749.3577273914822
van 4: 769.1224381882283
CPU times: user 1.03 s, sys: 18.5 ms, total: 1.05 s
Wall time: 1.05 s

Plot result paths

for path in result_paths:
    coords=get_coords(path)
    plt.plot(coords[0], coords[1], linewidth=2)
plt.scatter(data.x, data.y, c='r', s=30)
plt.scatter(data.iloc[0].x, data.iloc[0].y, marker='*', c='r', s=400)
plt.axis([-100,100,-100,100])
plt.grid(True)

Let's improve them with a 2-opt heuristic

%time improve_result_paths()

Starting path 1/4
exchanging 84-51,77-0 with 84-77,51-0 => save=11.952005586800027 => cost: 780.9700584117251
...
exchanging 66-51,147-126 with 66-147,51-126 => save=6.004072981449653 => cost: 732.2204189744828
Starting path 2/4
exchanging 35-69,29-0 with 35-29,69-0 => save=15.964056951390774 => cost: 777.9461580393865
...
exchanging 35-129,90-42 with 35-90,129-42 => save=14.532045586175727 => cost: 750.7675349695343
Starting path 3/4
exchanging 4-124,123-0 with 4-123,124-0 => save=19.23021437016152 => cost: 730.1275130213206
...
exchanging 65-14,52-33 with 65-52,14-33 => save=1.587339787656024 => cost: 701.2274950765787
Starting path 4/4
exchanging 19-21,54-0 with 19-54,21-0 => save=12.646642051225456 => cost: 756.4757961370028
CPU times: user 3.71 s, sys: 44.8 ms, total: 3.75 s
Wall time: 3.81 s

Plot improved paths

for path in result_paths:
    coords=get_coords(path)
    plt.plot(coords[0], coords[1], linewidth=2)
plt.scatter(data.x, data.y, c='r', s=30)
plt.scatter(data.iloc[0].x, data.iloc[0].y, marker='*', c='r', s=400)
plt.axis([-100,100,-100,100])
plt.grid(True)

Nice, clean paths. But with several overlaps.

In a real system, it should be analysed not only the cost (distance), but also the number of points in the route. Each one introduces a delay, as the driver has to park, walk, wait the customer, etc.

Scenario: farthest depot, several vans, KMeans clustering

Compute delivery clusters

from sklearn.cluster import KMeans

work_per_van=20 #how many deliveries each van will start, ie granularity; a better approach would be capacity.
model = KMeans(n_clusters=1+size//work_per_van, init='k-means++', random_state=0)
%time model.fit(data[['x', 'y']])

CPU times: user 91.4 ms, sys: 10.6 ms, total: 102 ms
Wall time: 111 ms





KMeans(algorithm='auto', copy_x=True, init='k-means++', max_iter=300,
    n_clusters=8, n_init=10, n_jobs=1, precompute_distances='auto',
    random_state=0, tol=0.0001, verbose=0)

Plot clusters

plt.scatter(data.x, data.y, marker='o', c=model.labels_, s=200)
plt.scatter(model.cluster_centers_[:,0], model.cluster_centers_[:,1], marker='+', c='r', s=500)
plt.scatter(data.iloc[0].x, data.iloc[0].y, marker='*', c='r', s=400)
plt.axis([-100,100,-100,100])
plt.grid(True)

Let's apply a greedy NN algorithm sequentially, on each of the clusters

%time result_paths = greedy_nn_clusters()    

van 1: 390.25298877225174
van 2: 687.7516428319581
van 3: 382.7248611037869
van 4: 505.82568809746317
van 5: 478.6115152097322
van 6: 644.9250908928291
van 7: 595.3986361541317
van 8: 665.6387979350595
CPU times: user 154 ms, sys: 16.6 ms, total: 170 ms
Wall time: 165 ms

Plot result paths

plt.scatter(data.x, data.y, marker='o', c=model.labels_, s=200)
plt.scatter(model.cluster_centers_[:,0], model.cluster_centers_[:,1], marker='+', c='r', s=500)
for path in result_paths:
    coords=get_coords(path)
    plt.plot(coords[0], coords[1], 'k-', linewidth=2)
plt.scatter(data.iloc[0].x, data.iloc[0].y, marker='*', c='r', s=400)
plt.axis([-100,100,-100,100])
plt.grid(True)

Let's improve them with a 2-opt heuristic

%time improve_result_paths()

Starting path 1/8
exchanging 77-4,83-29 with 77-83,4-29 => save=0.8194748732687884 => cost: 389.43351389898294
...
exchanging 6-143,128-89 with 6-128,143-89 => save=5.407506824441722 => cost: 324.93008877322006
Starting path 2/8
exchanging 30-5,34-0 with 30-34,5-0 => save=4.481243727960788 => cost: 683.2703991039973
...
exchanging 115-95,121-18 with 115-121,95-18 => save=8.322195150021315 => cost: 594.5702251203929
Starting path 3/8
exchanging 32-59,99-41 with 32-99,59-41 => save=4.114432829335271 => cost: 378.61042827445164
exchanging 17-78,26-87 with 17-26,78-87 => save=0.17155006280330554 => cost: 378.43887821164833
Starting path 4/8
exchanging 42-19,129-106 with 42-129,19-106 => save=3.8728268587745447 => cost: 501.9528612386886
...
exchanging 37-94,54-0 with 37-54,94-0 => save=1.3747236340181104 => cost: 495.2855372216575
Starting path 5/8
exchanging 0-147,123-65 with 0-123,147-65 => save=11.990854640529989 => cost: 466.6206605692022
...
exchanging 135-65,123-0 with 135-123,65-0 => save=3.919229108305899 => cost: 435.93197839715947
Starting path 6/8
exchanging 56-21,15-0 with 56-15,21-0 => save=6.923848236857282 => cost: 638.0012426559717
...
exchanging 102-43,118-40 with 102-118,43-40 => save=15.652843812142969 => cost: 602.6373464143627
Starting path 7/8
exchanging 63-24,107-110 with 63-107,24-110 => save=11.657922849361384 => cost: 583.7407133047703
...
exchanging 112-36,91-125 with 112-91,36-125 => save=4.722248082275083 => cost: 559.3163570654888
Starting path 8/8
exchanging 48-111,98-137 with 48-98,111-137 => save=11.5162692633929 => cost: 654.1225286716666
...
exchanging 8-132,137-98 with 8-137,132-98 => save=7.784370993699362 => cost: 615.4951078634247
CPU times: user 2.17 s, sys: 49.8 ms, total: 2.22 s
Wall time: 2.39 s

Plot improved paths

plt.scatter(data.x, data.y, marker='o', c=model.labels_, s=200)
plt.scatter(model.cluster_centers_[:,0], model.cluster_centers_[:,1], marker='+', c='r', s=500)
for path in result_paths:
    coords=get_coords(path)
    plt.plot(coords[0], coords[1], 'k-', linewidth=2)
plt.scatter(data.iloc[0].x, data.iloc[0].y, marker='*', c='r', s=400)
plt.axis([-100,100,-100,100])
plt.grid(True)

Very cool subsections! Aproximately the same number of deliveries in each route!

But some routes are way longer than others. Some vans maybe will have to do an extra hour or two.

---

That's it!

Thank you for getting this far, hope you liked! :)