The Vehicle Routing Problem with Service Level Constraints

TL;DR

This paper introduces a new vehicle routing problem with service level constraints, proposing algorithms that effectively balance costs and service levels, and demonstrating high-quality solutions and new optimal results through computational experiments.

Contribution

It presents a novel formulation and algorithms for the vehicle routing problem with service level constraints, extending existing models and providing new optimal solutions.

Findings

High-quality heuristic solutions for complex problem instances

Matching all optimal solutions for small and medium benchmarks

Improving algorithms for related vehicle routing problems

Abstract

We consider a vehicle routing problem which seeks to minimize cost subject to service level constraints on several groups of deliveries. This problem captures some essential challenges faced by a logistics provider which operates transportation services for a limited number of partners and should respect contractual obligations on service levels. The problem also generalizes several important classes of vehicle routing problems with profits. To solve it, we propose a compact mathematical formulation, a branch-and-price algorithm, and a hybrid genetic algorithm with population management, which relies on problem-tailored solution representation, crossover and local search operators, as well as an adaptive penalization mechanism establishing a good balance between service levels and costs. Our computational experiments show that the proposed heuristic returns very high-quality solutions for this difficult problem, matches all optimal solutions found for small and medium-scale benchmark instances, and improves upon existing algorithms for two important special cases: the vehicle routing problem with private fleet and common carrier, and the capacitated profitable tour problem. The branch-and-price algorithm also produces new optimal solutions for all three problems.

Tables

Table 1. Table 1: Synthesis of state-of-the-art heuristics and exact approaches for single- and multi-vehicle routing problems with profits.

	Min $\alpha \times$cost $-\beta \times$service level		Max service level s.t. cost constraints		Min cost s.t. service level constraints
	Ib) Profitable tour problem (PTP)		IIb) Orienteering problem (OP)		IIIb) Prize-collecting TSP (PCTSP)
	No dedicated heuristics, but:		Heuristic:		Heuristic:
	Reductions to elementary shortest path with		GRASP	[14]	Tabu search using GENIUS neighborhood	[49]
	possible negative cycles, or asymmetric TSP		Guided local search	[59]	Clustering search	[15]
T	Various approximation algorithms	[3]	Ant colony optimization and VNS	[53]	Lagrangian heuristic	[20]
S	Extensions into PCTSP with profits		Tabu search	[28]	Tabu Search (steel industry)	[45]
P	considered in the objective		and others…
	Exact:		Exact:		Exact:
	Branch-and-cut (+ capacity constraint)	[35]	Branch-and-cut	[24]	Branch-and-cut	[11]
	Bounding method via Lagrangian relax.	[20]	Branch-and-cut + Lagrangian relaxation	[36]	Bounding method via Lagrangian relax.	[20]
					Formulation + Additive bounds	[25]
	Ia) VRPPFCC and CPTP		IIa) Team-orienteering problem (TOP)		IIIa) Prize collecting VRP (PCVRP)
	Heuristic VRPPFCC:		Heuristic:		Limited number of works on this variant
	Hybrid GA with large neighborhoods	[66]	Pareto mimic algorithm	[38]
	Adaptive variable neighborhood search	[32]	Hybrid GA with large neighborhoods	[66]	Heuristic:
	Adaptive variable neighborhood search	[57]	Multi-start simulated annealing	[44]	Adaptive VNS for a related problem
V	Tabu with possible ejection chains	[18, 51]	Augmented large neighborhood search	[40]	with a minimum delivery level	[56]
R	and others…		Particle swarm optimization	[19]	Guided local search (steel industry)	[68]
P			Path relinking	[54]	Iterated local search (steel industry)	[58]
	Heuristic CPTP:		Guided local search	[59]	Two-stage self-adaptive VNS	[43]
	VNS and tabu search	[2]	Ant colony optimization	[37]
			and others…
	Exact:		Exact:		Exact:
	Branch-and-price	[1, 2]	Branch-and-price	[1, 2, 13]	None to our knowledge
			Branch-and-cut-and-price	[39]
			Branch-and-cut	[22]

Table 2. Table 2: Results of the exact methods for the VRP-SL instances

				Compact Formulation						Branch-and-Price
Data	$𝐧$	$𝐦$	${𝐧}_{\text{min}}$	LB0	T0(s)	LB	Tree	Gap	T(s)	LB0	T0(s)	LB	Tree	Gap	T(s)
A1	31	5	19.4	643.1	0.50	701.6	756k	0.78	3423.93	685.9	0.13	707.4	319	0.00	33.82
A2	32	6	17.6	610.0	0.37	652.7	650k	0.77	2288.52	643.8	0.12	656.1	2k	0.30	1443.06
A3	35	5	21.8	627.1	0.40	666.2	407k	0.00	1806.92	652.4	0.15	666.2	1k	0.00	224.77
A4	36	6	18.0	757.7	0.56	802.4	386k	0.28	2482.49	792.6	0.13	805.0	103	0.00	5.93
A5	38	5	18.8	629.3	0.55	667.7	584k	0.25	3212.89	649.7	0.21	667.9	3k	0.22	1861.53
A6	43	7	25.4	739.9	0.88	772.7	414k	1.94	4477.12	772.9	0.25	788.6	4k	0.11	2890.31
A7	44	7	28.2	956.9	1.66	989.3	560k	3.25	7200.00	1006.1	0.31	1020.8	8k	0.21	4118.13
A8	47	7	26.6	874.5	1.40	905.9	670k	4.47	7200.00	926.9	0.21	944.0	7k	0.44	3394.61
A9	53	7	34.2	993.4	1.94	1015.1	439k	5.30	7200.00	1045.0	0.67	1066.0	8k	0.64	4905.96
A10	59	9	37.4	1136.6	3.18	1163.1	390k	5.67	7200.00	1208.3	0.61	1226.3	14k	0.60	5981.17
A11	61	8	33.0	1033.1	3.62	1057.3	249k	5.90	7200.00	1094.7	0.74	1112.1	7k	1.04	4953.56
A12	62	9	35.2	1055.1	2.54	1080.0	350k	4.58	7200.00	1113.6	0.60	1128.3	18k	0.35	6236.87
A13	64	9	41.6	1030.2	2.17	1052.7	529k	4.40	7200.00	1074.9	0.56	1090.4	11k	1.00	5764.22
A14	79	10	43.2	1432.9	6.00	1453.1	177k	4.70	7200.00	1494.4	1.39	1508.6	7k	1.09	7200.00
B1	30	5	21.6	526.0	0.38	576.3	1786k	1.42	5794.55	550.3	0.15	569.7	8k	2.37	5762.20
B2	34	5	16.8	695.4	0.26	736.0	1128k	1.51	3476.90	709.6	0.31	733.9	6k	1.91	5760.99
B3	38	6	18.0	434.4	0.52	460.6	1673k	7.24	7200.00	466.5	0.60	477.9	10k	3.84	7200.00
B4	42	6	22.0	618.2	0.83	650.7	766k	2.29	6352.01	652.0	0.46	663.2	3k	0.44	2705.18
B5	44	5	22.8	548.8	0.82	587.0	1040k	6.13	7200.00	606.2	0.77	614.8	7k	1.68	7200.00
B6	49	7	30.0	599.1	1.17	630.1	1049k	5.74	7200.00	646.5	0.91	657.2	7k	1.69	5140.14
B7	50	7	35.4	863.7	0.58	908.5	1259k	4.87	7200.00	916.5	0.91	933.3	6k	2.23	5950.49
B8	55	7	37.2	563.0	1.29	574.3	1117k	9.76	7200.00	605.9	1.45	613.1	7k	3.65	7200.00
B9	56	9	35.2	1291.6	1.34	1316.1	724k	4.32	7200.00	1332.6	0.63	1347.1	14k	2.04	5761.09
B10	63	9	34.0	726.2	1.97	761.4	731k	6.22	7200.00	781.9	1.07	793.5	8k	2.27	7200.00
B11	66	10	42.0	878.6	1.85	907.8	787k	4.54	7200.00	923.4	1.21	936.1	8k	1.57	7200.00
B12	77	10	43.0	1003.3	3.77	1020.3	450k	6.49	7200.00	1060.6	1.91	1069.1	9k	2.02	7200.00
p03	100	3	53.8	542.1	3.48	554.4	176k	0.11	2734.21	545.0	92.27	548.4	209	1.21	7200.00
p06	50	2	30.0	370.4	0.44	387.0	1k	0.00	5.26	374.1	11.61	381.0	418	1.54	7200.00
p07	75	2.8	42.8	506.9	1.41	519.6	69k	0.00	301.26	511.4	28.14	515.0	519	0.89	6489.66
p08	100	3	55.2	531.3	3.33	544.3	109k	0.09	1784.18	534.4	44.19	537.6	232	1.34	7200.00
p09	150	4.6	88.2	669.4	17.11	684.5	2170k	1.58	7200.00	681.9	306.57	684.0	79	1.64	7200.00
p10	199	6	114.4	773.8	61.80	783.4	16k	3.31	7200.00	798.2	699.11	799.6	45	1.32	7200.00
p13	120	3	75.0	529.0	10.03	538.9	114k	15.77	7200.00	587.3	837.90	588.9	27	7.96	7200.00
p14	100	4	55.2	491.8	4.83	504.9	227k	14.13	7200.00	560.3	90.53	568.5	575	3.31	7200.00
p15	100	3.8	65.6	493.4	4.72	506.9	212k	13.92	7200.00	565.4	108.80	572.1	557	2.85	7200.00
p16	199	6.2	108.2	769.7	66.16	782.2	204k	2.69	7200.00	792.7	571.70	794.0	37	1.22	7200.00

Table 3. Table 3: Performance of the HGS on the VRP-SL instance sets

Inst	n	m	Wor-10	Avg-10	Best-10	Gap${}_{\text{BKS}}$	Gap${}_{\text{LB}}$	T(s)	BKS	BKLB
A1	31	5	707.4	707.40	707.4	0.00	0.00	9.06	707.4	707.4
A2	32	6	658.2	658.20	658.2	0.00	0.32	8.72	658.2	656.1
A3	35	5	666.2	666.20	666.2	0.00	0.00	9.70	666.2	666.2
A4	36	6	805.0	805.00	805.0	0.00	0.00	9.51	805.0	805.0
A5	38	5	669.4	669.40	669.4	0.00	0.23	10.18	669.4	667.9
A6	43	6	789.6	789.60	789.6	0.00	0.13	12.78	789.6	788.6
A7	44	7	1023.0	1023.00	1023.0	0.00	0.22	13.90	1023.0	1020.8
A8	47	7	948.2	948.20	948.2	0.00	0.44	12.99	948.2	944.0
A9	53	7	1073.2	1073.20	1073.2	0.00	0.68	18.49	1073.2	1066.0
A10	59	9	1234.8	1234.50	1234.4	0.01	0.67	21.21	1234.4	1226.3
A11	61	8	1124.8	1124.10	1124.0	0.01	1.08	20.69	1124.0	1112.1
A12	62	10	1132.6	1132.44	1132.4	0.00	0.36	19.54	1132.4	1128.3
A13	64	9	1101.8	1101.70	1101.6	0.01	1.04	22.64	1101.6	1090.4
A14	79	10	1526.2	1525.30	1525.2	0.01	1.11	32.05	1525.2	1508.6
B1	30	5	584.8	584.80	584.8	0.00	0.57	9.53	584.8	581.5
B2	34	5	748.2	748.20	748.2	0.00	1.10	8.34	748.2	740.0
B3	38	5	497.2	497.20	497.2	0.00	3.98	10.24	497.2	478.2
B4	42	6	666.2	666.20	666.2	0.00	0.46	12.00	666.2	663.2
B5	44	5	625.0	625.00	625.0	0.00	1.66	11.71	625.0	614.8
B6	49	7	668.8	668.80	668.8	0.00	1.77	14.25	668.8	657.2
B7	50	7	956.2	956.20	956.2	0.00	2.38	15.13	956.2	934.0
B8	55	7	637.4	637.40	637.4	0.00	3.97	18.26	637.4	613.1
B9	56	9	1378.2	1376.72	1376.2	0.04	2.03	22.17	1376.2	1349.4
B10	63	9	812.0	812.00	812.0	0.00	2.33	20.20	812.0	793.5
B11	66	10	951.8	951.56	951.4	0.02	1.65	24.15	951.4	936.1
B12	77	10	1092.6	1091.94	1091.4	0.05	2.13	28.92	1091.4	1069.1
p03	100	–	555.0	555.00	555.0	0.00	0.10	22.67	555.0	554.4
p06	50	–	387.0	387.00	387.0	0.00	0.00	10.74	387.0	387.0
p07	75	–	519.6	519.60	519.6	0.00	0.00	17.51	519.6	519.6
p08	100	–	544.8	544.80	544.8	0.00	0.09	22.50	544.8	544.3
p09	150	–	695.6	695.60	695.6	0.00	1.50	46.12	695.6	685.3
p10	199	–	810.8	810.44	810.2	0.03	1.36	86.30	810.2	799.6
p13	120	–	640.8	640.20	639.8	0.06	8.70	44.05	639.8	588.9
p14	100	–	588.0	588.00	588.0	0.00	3.43	27.04	588.0	568.5
p15	100	–	588.8	588.80	588.8	0.00	2.93	26.78	588.8	572.1
p16	199	–	805.0	804.24	803.8	0.05	1.29	92.85	803.8	794.0
All						0.01	1.38	22.58

Table 4. Table 4: Performance of the HGS on VRPPFCC benchmark instances, in comparison with the current state-of-the-art algorithms

			UGHS – [66]			AVNS – [32]			This paper – HGS				This paper – B&P
Instance	$n$	$m$	Avg-10	Best-10	T(s)	Avg-10	Best-10	T(s)	Avg-10	Best-10	T(s)		LB0	LB	T(s)	BKS
CE-01	50	4	1119.66	1119.47	38.5	1119.47	1119.47	47.4	1119.47	1119.47	13.01		1112.17	1119.47	154.30	1119.47
CE-02	75	9	1815.63	1814.52	59.6	1817.06	1814.52	95.8	1814.52	1814.52	24.00		1795.13	1807.79	7200.00	1814.52
CE-03	100	6	1922.88	1919.05	476.8	1925.08	1919.05	236.7	1923.92	1919.05	45.16		1896.81	1904.30	7200.00	1919.05
CE-04	150	9	2509.82	2505.39	934.7	2518.14	2509.81	642.7	2510.35	2505.39	98.95		2478.50	2483.50	7200.00	2505.39
CE-05	199	13	3095.58	3081.59	1289.3	3101.40	3090.49	1437.2	3094.78	3084.40	213.69		3051.13	3053.92	7200.00	3081.59
CE-06	50	4	1207.47	1207.47	38.5	1207.47	1207.47	44.2	1207.47	1207.47	12.42		1199.81	1207.47	196.41	1207.47
CE-07	75	9	2012.33	2006.52	73.2	2009.58	2004.53	98.6	2006.32	2004.53	25.71		1984.80	1997.90	7200.00	2004.53
CE-08	100	6	2057.57	2052.05	500.3	2059.56	2052.05	239.6	2058.18	2052.05	48.30		2029.89	2037.66	7200.00	2052.05
CE-09	150	10	2428.19	2425.32	1241.0	2426.35	2420.71	953.9	2422.88	2420.57	123.43		2389.86	2394.60	7200.00	2419.84
CE-10	199	13	3387.12	3381.67	1229.9	3388.22	3373.84	1704.1	3384.03	3374.02	238.95		3340.29	3343.29	7200.00	3373.84
CE-11	120	6	2331.13	2330.94	1202.9	2330.94	2330.94	399.0	2330.95	2330.94	71.94		2317.21	2321.26	7200.00	2330.94
CE-12	100	8	1953.13	1952.86	150.8	1952.86	1952.86	183.6	1953.27	1952.86	34.28		1939.32	1946.63	7200.00	1952.86
CE-13	120	6	2859.07	2858.83	1184.9	2858.92	2858.83	437.8	2858.94	2858.83	65.11		2843.90	2847.50	7200.00	2858.83
CE-14	100	7	2213.02	2213.02	189.8	2213.78	2213.02	185.2	2213.41	2213.02	33.11		2193.70	2202.21	7200.00	2213.02
G-01	240	7	14151.51	14131.18	2405.9	14163.43	14129.48	2433.9	14135.93	14112.39	210.40		14024.30	14038.00	7200.00	14111.95
G-02	320	8	19190.77	19166.58	2409.9	19254.23	19140.69	6630.0	19134.06	19125.35	302.27		18974.40	18978.70	7200.00	19140.69
G-03	400	8	24588.29	24409.02	2418.1	24566.00	24406.67	11062.8	24372.81	24323.64	431.70		—	—	7200.00	24368.29
G-04	480	8	34517.47	34362.8	2421.1	34425.00	34231.56	15875.4	34116.56	34076.38	401.61		—	—	7200.00	34231.56
G-05	200	4	14296.07	14223.63	2408.0	14261.06	14229.50	1591.3	14268.87	14223.63	71.98		—	—	7200.00	14223.63
G-06	280	5	21488.29	21396.60	2411.4	21440.38	21357.16	4337.2	21376.19	21357.16	125.28		—	—	7200.00	21357.16
G-07	360	7	23463.05	23373.38	2414.9	23440.51	23263.22	7585.3	23285.33	23233.97	354.48		—	—	7200.00	23263.22
G-08	440	8	29918.06	29823.18	2415.6	29864.19	29657.38	12316.5	29603.54	29556.13	415.88		—	—	7200.00	29657.38
G-09	255	11	1332.63	1328.65	2323.4	1325.79	1320.29	3852	1324.06	1320.86	335.90		1308.14	1308.59	7200.00	1319.72
G-10	323	13	1603.82	1597.61	2342.3	1592.14	1588.05	6922.9	1587.42	1584.73	647.96		1567.66	1568.01	7200.00	1583.50
G-11	399	14	2192.68	2182.01	2405.2	2172.45	2163.50	12303.3	2166.52	2160.43	861.43		2137.91	2138.06	7200.00	2159.78
G-12	483	15	2529.84	2522.64	2407.2	2493.94	2483.06	19555.0	2486.56	2481.76	956.42		2455.51	2455.64	7200.00	2479.62
G-13	252	21	2261.50	2258.02	1412.7	2264.92	2261.66	2260.2	2266.41	2263.48	325.29		2232.70	2232.87	7200.00	2258.02
G-14	320	23	2687.50	2683.73	1935.7	2689.99	2684.66	4838.5	2696.11	2689.67	682.99		2654.94	2655.23	7200.00	2683.73
G-15	396	26	3152.00	3145.11	2301.4	3156.84	3150.67	8198.0	3153.76	3144.51	858.02		3107.12	3107.42	7200.00	3145.11
G-16	480	29	3632.04	3620.71	2450.1	3633.76	3624.56	12741.7	3635.54	3624.91	954.65		3578.85	3579.04	7200.00	3620.71
G-17	240	18	1671.72	1666.31	1805.3	1672.81	1666.31	1878	1667.81	1666.31	160.84		1666.31	1666.31	19.69	1666.31
G-18	300	22	2733.12	2730.55	2035.0	2734.86	2731.28	3916.6	2733.49	2730.97	459.92		2717.24	2717.95	7200.00	2730.55
G-19	360	26	3504.26	3497.20	1989.5	3499.55	3494.28	5225.4	3495.68	3491.24	839.89		3477.87	3478.79	7200.00	3494.27
G-20	420	31	4319.37	4312.45	2523	4314.48	4307.63	8081.6	4307.61	4300.77	1115.09		4277.33	4278.20	7200.00	4306.85
Avg. Gap(%)			0.445	0.210		0.345	0.058		0.141	-0.012			—	—
Avg. T(s)					1583.70			4656.22							6575.60
CPU			Xe 3.07GHz			Intel i5 2.6GHz			Intel i7 3.4GHz				Intel i7 3.4GHz

Table 5. Table 5: Performance of HGS on the CPTP benchmark instances

		UGHS – [66]			ILS – [66]			This paper – HGS
Instance	$n$	Avg-10	Best-10	T(s)	Avg-10	Best-10	T(s)	Avg-10	Best-10	T(s)		BKS
p03	100	254.07	254.07	215.82	253.99	254.07	191.63	254.07	254.07	12.32		254.07
p06	50	129.13	129.13	30.78	129.09	129.11	23.35	129.13	129.13	7.18		129.13
p07	75	192.51	192.56	103.94	192.37	192.56	85.80	192.52	192.56	8.36		192.56
p08	100	254.07	254.07	216.28	253.93	254.07	191.49	253.90	254.07	12.17		254.07
p09	150	319.61	319.72	295.34	319.23	319.72	289.75	319.67	319.72	18.44		319.72
p10	199	387.87	388.41	303.55	386.67	387.73	309.23	388.79	388.85	19.90		388.41
p13	120	180.20	180.39	235.98	178.33	180.32	255.97	180.08	180.32	15.15		180.39
p14	100	246.24	246.24	116.26	246.23	246.24	111.83	246.24	246.24	10.78		246.24
p15	150	327.99	328.36	295.02	326.91	327.81	289.85	328.28	328.37	15.24		328.36
p16	199	393.09	393.75	303.88	392.21	393.32	309.74	394.03	394.04	21.16		393.75
Avg. Gap(%)		0.029	0.000		0.172	0.014		0.012	-0.002
Avg. T(s)				211.68			205.86			14.07
CPU		Xe 3.07GHz			Xe 3.07GHz			Intel i7 3.4GHz

Table 6. Table 6: Performance of the proposed branch-and-price algorithm on the CPTP benchmark instances

		B&P1 – [1]			B&P2 – [1]			This paper – B&P
Instance	$n$	UB	LB	T(s)	UB	LB	T(s)	UB0	UB	LB	T(s)	BKUB
p03	100	256.90	254.07	1135.31	256.82	147.01	1135.38	258.17	255.36	254.07	1116.00	256.82
p06	50	129.75	129.13	625.08	129.64	118.09	602.38	131.77	129.39	129.13	287.89	129.64
p07	75	193.18	192.56	622.38	193.11	163.49	613.92	194.21	192.82	192.56	285.06	193.11
p08	100	256.89	254.07	1134.31	256.82	147.01	1135.00	258.17	255.36	254.07	1116.05	256.82
p09	150	324.25	316.33	1163.69	324.20	160.90	1163.85	323.66	320.83	319.72	1123.65	324.20
p10	199	392.03	377.01	1117.77	391.66	176.43	1117.77	392.07	390.43	388.85	596.37	391.66
p13	120	191.74	167.78	1916.08	191.70	116.06	1916.77	186.72	183.93	180.32	1216.18	191.70
p14	100	255.47	237.02	1110.31	255.41	116.24	1110.31	248.67	247.20	246.24	556.22	255.40
p15	150	331.98	324.19	893.54	331.93	190.18	893.77	332.36	330.67	328.37	603.52	331.93
p16	199	398.57	379.03	1118.46	398.27	178.90	1118.38	398.58	396.47	394.04	941.83	398.27
Average		273.07	263.12	1083.69	272.96	151.43	1080.75	272.44	270.25	268.74	784.28	272.96
CPU		Xe 2.26GHz			Xe 2.26GHz			Intel i7 3.4GHz

Table 7. Table 7: Sensitivity Analysis on the components of the HGS

	VRP-SL  Set S1			VRP-SL   Set S2			VRPPFCC			CPTP
	Best-10	Avg-10	T(s)	Best-10	Avg-10	T(s)	Best-10	Avg-10	T(s)	Best-10	Avg-10	T(s)
Standard	0.00	0.01	16.01	0.00	0.01	39.66	0.06	0.21	340.00	0.00	0.02	14.07
EOX	0.02	0.05	19.98	0.03	0.10	66.61	1.69	2.04	225.45	0.02	0.06	20.39
No SL	0.00	0.01	16.70	0.00	0.02	41.68	0.07	0.20	360.61	0.02	0.04	15.22
No INF	0.05	0.12	15.05	0.07	0.14	40.84	0.15	0.38	344.67	0.02	0.05	12.70
No DIV	0.01	0.05	13.62	0.01	0.08	29.04	0.12	0.37	155.33	0.02	0.17	11.68
No Repair	0.00	0.02	12.51	0.01	0.03	34.85	0.05	0.20	337.18	0.00	0.02	12.40
Pop $\downarrow$	0.00	0.01	12.99	0.00	0.02	31.16	0.05	0.23	343.94	0.00	0.03	11.24
Pop $\uparrow$	0.00	0.01	26.48	0.00	0.02	60.53	0.12	0.26	545.01	0.00	0.01	25.14
Feas $\uparrow$	0.00	0.01	15.12	0.00	0.02	40.01	0.10	0.26	348.12	0.01	0.03	14.06
Feas $\uparrow \uparrow$	0.01	0.04	14.93	0.01	0.05	40.61	0.16	0.37	380.05	0.01	0.05	14.23

Equations

$$\text{Minimize}\sum _{(v_i,v_j)\in \overline{\mathcal{E}}}{d}_{ij}{x}_{ij}$$

$$\text{Subject to}\sum _{v_i\in {\mathcal{V}}_k}{s}_i{y}_i$$

$$\sum _{v_i\in \overline{\mathcal{V}},i<k}{x}_{ik}+\sum _{v_j\in \overline{\mathcal{V}},j>k}{x}_{kj}$$

$$\sum _{v_j\in \overline{\mathcal{V}}}(f_{ji}-f_{ij})$$

$$\sum _{v_j\in {\mathcal{V}}^{\prime }}{f}_{0j}$$

$$\sum _{v_j\in {\mathcal{V}}^{\prime }}{f}_{n+1j}$$

$$f_{ij}+f_{ji}$$

$$x_{ij}$$

$$y_i$$

$$f_{ij}$$

$$z$$

$$\sum _{v_i\in {\mathcal{V}}_k}{s}_i{y}_i$$

$$f_{ij}\ge q_j{x}_{ij}\text{ and }$$

$$Z_{k}=\min_{y_{i}\in\{0,1\}}\left\{\sum_{v_{i}\in{\cal V}_{k}}q_{i}y_{i}\hskip 5.69046pt\bigg{|}\hskip 5.69046pt\sum_{v_{i}\in{\cal V}_{k}}s_{i}y_{i}\geq\alpha_{k}\sum_{v_{i}\in{\cal V}_{k}}s_{i}\right\}$$

$$\sum _{v_i\in V_k,v_j\in \mathcal{V}\setminus {\mathcal{V}}_k}{x}_{ij}$$

$$\text{Minimize}\sum _{r\in \Omega }{c}_r{\lambda }_r$$

$$\text{Subject to}\sum _{v_i\in {\mathcal{V}}_k}{s}_i{y}_i$$

$$\sum _{r\in \Omega }{a}_i^r{\lambda }_r$$

$$\sum _{r\in \Omega }{\lambda }_r$$

$${\lambda }_r$$

$$y_i$$

$${\overline{c}}_r$$

$$\mathcal{L}(P^{\prime })=\left(v_j,\overline{c}(P)+{\overline{c}}_{ij},q(P)+q_{v_j},\Pi (P)\cap N{G}_j\cup \{v_j\}\right)$$

$$\{v(P_1)=v(P_2)\}\wedge \{\overline{c}(P_1)\le \overline{c}(P_2)\}\wedge \{q(P_1)\le q(P_2)\}\wedge \{\Pi (P_1)\subseteq \Pi (P_2)\}.$$

$$\gamma$$

$${\beta }_i$$

$$\sum _{r\in \Omega }{b}_{ij}^r{\lambda }_r={\overline{x}}_{ij}$$

$${\overline{c}}_r=\sum _{(v_i,v_j)\in \mathcal{E}}{b}_{ij}^r\left(c_{ij}-\frac{{\beta }_i+{\beta }_j}{2}-{\rho }_{ij}\right)$$

$$\varphi (r)={\varphi }^D(r)-{\varphi }^P(r)+w^Q\max \{0,{\varphi }^Q(r)-Q\}.$$

$${\varphi }^{\text{textsc}cost}(S)$$

$$\text{with }{\varphi }^k(S)$$

$${\varphi }_P^{\text{textsc}div}(S)=\frac{{\sum }_{S_2\in N_P(S)}\delta (S,S_2)}{n^{\text{textsc}close}},$$

$$\delta (S_1,S_2)=1-\frac{|E(S_1)\cap E(S_2)|}{|E(S_1)\cup E(S_2)|}.$$

$$\phi^{\textsc{bias}}_{P}(S)=R(S,P,\phi^{\textsc{cost}})+\Big{(}1-\frac{n^{\textsc{elite}}}{|P|}\Big{)}R(S,P,\phi^{\textsc{div}}).$$

$$w^c=\{\begin{array}{ll}{w}^c\times 1.2& \text{ if }{\xi }^c\le {\xi }^{\text{textsc}min},\\ w^c\times 0.85& \text{ if }{\xi }^c\ge {\xi }^{\text{textsc}max},\\ w^c& \text{ otherwise.}\end{array}$$

$$m=\lceil \frac{{\sum }_{k=1}^K(Q_k^{\text{textsc}min}+Q_k^{\text{textsc}max})}{2Q}\rceil ,$$

Full text

The Vehicle Routing Problem with Service Level Constraints

Teobaldo Bulhões, Minh Hoàng Hà, Thibaut Vidal

The Vehicle Routing Problem with Service Level Constraints

**Teobaldo Bulhões ***

Instituto de Computação, Universidade Federal Fluminense, Niterói, Brazil

tbulhoes@ic.uff.br

Minh Hoàng Hà

FPT Technology Research Institute, FPT University, Hanoi, Viet Nam

hoanghm@fpt.edu.vn

Rafael Martinelli

Departamento de Engenharia Industrial, Pontifícia Universidade Católica do Rio de Janeiro, Brazil

martinelli@puc-rio.br

Thibaut Vidal

Departamento de Informática, Pontifícia Universidade Católica do Rio de Janeiro, Brazil

vidalt@inf.puc-rio.br

Abstract. We consider a vehicle routing problem which seeks to minimize cost subject to service level constraints on several groups of deliveries. This problem captures some essential challenges faced by a logistics provider which operates transportation services for a limited number of partners and should respect contractual obligations on service levels. The problem also generalizes several important classes of vehicle routing problems with profits. To solve it, we propose a compact mathematical formulation, a branch-and-price algorithm, and a hybrid genetic algorithm with population management, which relies on problem-tailored solution representation, crossover and local search operators, as well as an adaptive penalization mechanism establishing a good balance between service levels and costs. Our computational experiments show that the proposed heuristic returns very high-quality solutions for this difficult problem, matches all optimal solutions found for small and medium-scale benchmark instances, and improves upon existing algorithms for two important special cases: the vehicle routing problem with private fleet and common carrier, and the capacitated profitable tour problem. The branch-and-price algorithm also produces new optimal solutions for all three problems.

Keywords. Routing, collaborative logistics, service level constraints, integer programming, genetic algorithms.

1 Introduction

In a recent industrial application, the authors have been confronted with a complex variant of the vehicle routing problem (VRP) which received, until now, only limited attention in the academic literature. We take the viewpoint of a third-party logistics provider (3PL), which operates long-haul transportation services for a number of business partners. The company operates on a planning horizon and delivers products to various delivery locations as requested by the partners a few days in advance. A strict delivery deadline, in the form of a last possible delivery day, is set for each transportation request. The company has established agreements with each partner specifying, among others, a minimum level of on-time deliveries for its group of requests.

The efficient management of collaborative logistics has stimulated a rich set of studies over the years [34, 42, 55]. In our case, even in the presence of precise information on delivery requests, the company faces the following optimization challenge: because of limited available resources, i.e., fleet size and vehicle capacity, all services cannot be realistically fulfilled. It is necessary to select a subset of deliveries and determine cost-efficient vehicle routes in such a way that the overall activity is profitable and that the agreements with the partners are respected. A natural strategy of the company consists in attempting to fulfill the service levels on a rolling horizon with some safety margin, hence ensuring overall satisfaction of contractual clauses on a larger time period (e.g., one month of activity). As such, the firm seeks to balance quality of service and operational costs, subject to a minimum threshold on quality of service for some groups of requests. These group requirements create linking constraints between the service selection decisions, which need to be carefully considered during routing optimization.

We now introduce a simple and deterministic variant of the VRP which captures some essential decisions in this situation. The study of this simplified problem will help to identify key properties and methods. The VRP with service levels (VRP-SL) can be formulated as follows. Let ${\cal G}=({\cal V},{\cal E})$ be a complete undirected graph with $|{\cal V}|=n+1$ nodes. The node $v_{0}\in{\cal V}$ represents a depot, where a fleet of $m$ identical vehicles is based. Each other node $v_{i}$ for $i\in\{1,\ldots,n\}$ represents a customer, associated with a demand $q_{i}$, a profit $p_{i}$, and a service weight $s_{i}$ which represents its relative importance in the group service level constraint.

The set of customers is distributed into $K$ subsets: ${\cal V}-\{v_{0}\}=\bigcup_{k=1,\dots,K}{\cal V}_{k}$, such that ${\cal V}_{k}\cap{\cal V}_{k^{\prime}}=\varnothing$ for any $k\neq k^{\prime}$. Each subset represents the deliveries of one partner and is associated with a requested service level $\alpha_{k}$. Any edge $(i,j)\in{\cal E}$ represents a possible trip between a node $v_{i}\in{\cal V}$ and a node $v_{j}\in{\cal V}$ with a distance cost $d_{ij}$. The goal of the VRP-SL is to find up to $m$ vehicle routes starting and ending at the depot, such that

• –

each customer is serviced at most one time,

• –

the total demand quantity of any route does not exceed a vehicle capacity $Q$,

• –

the service level of each group $k$ is attained, i.e., the total service weight of the deliveries to this group reaches $\alpha_{k}\sum_{v_{i}\in{\cal V}_{k}}s_{i}$, and

• –

the sum of travel costs and lost profits is minimized.

This problem belongs to the wide class of vehicle routing problems with profits, which also includes the team orienteering problem (TOP), the profitable VRP (VRPP), the VRP with private fleet and common carrier (VRPPFCC) and the capacitated profitable tour problem (CPTP). Interestingly, as highlighted in Section 2, this problem fills a gap in the literature, since most known multi-vehicle problems with customer selections either aim to maximize service levels subject to distance constraints (TOP) or seek a weighted optimization of distance and service levels, through penalties for outsourcing or lost profits (VRPPFCC and CPTP). To this date, very few works on deterministic settings [58, 68, 43] have addressed multi-vehicle routing optimization subject to a service level (SL) constraint. Finally, the VRP-SL is finally a natural extension of the generalized VRP (GVRP, see [29, 5, 10, 31]), which also models a rich set of VRP applications.

To address the VRP-SL, we introduce a compact ILP formulation and a branch-and-price algorithm which can solve to optimality small- and medium-scale instances, as well as a hybrid population metaheuristic inspired by the Unified Hybrid Genetic Search (UHGS) framework of [62, 64]. Previous applications of UHGS have led to efficient algorithms for several VRPs, including the GVRP, TOP, VRPP and VRPPFCC, among others [64, 66]. However, UHGS relies heavily on the fact that problem objectives and constraints are either cumulated on or separately applied to each route, hence allowing to optimize customer-selection decisions independently via a dedicated route evaluation operator. This decomposition does not apply to the VRP-SL without a Lagrangian relaxation of the group constraints, which would make it impossible to find the optimal solution in some cases and impede the overall method performance. Because of these characteristics, we had to revise most of the operators of the method while keeping the general principles. We thus use a two-chromosome solution representation, which contains a service level chromosome and a giant-tour chromosome. We derive a new crossover, use dedicated local search moves as well as a penalization strategy to represent, optimize and inherit customer-selection decisions. The contributions of this article are the following:

We introduce the VRP-SL, a rich VRP connected with important applications in collaborative logistics, generalizing many classes of routing problems with profits. 2. 2.

We propose a compact mathematical formulation for the problem, a branch-and-price algorithm, as well as an efficient hybrid genetic search (HGS) which exploits problem-tailored selection representation, crossover, local searches, and penalty management operators to find a good balance between service levels and costs. 3. 3.

We conduct extensive computational experiments to evaluate the performance of the proposed methods on new benchmark instances for the VRP-SL as well as classical benchmark instances for the VRPPFCC and CPTP. The proposed HGS finds all known optimal solutions for the considered problems, outperforms previous methods for the VRPPFCC and CPTP, and generates solutions of consistent quality on large instances. Several solutions of the VRPPFCC and CPTP are also proven optimal for the first time by the branch-and-price algorithm.

The remaining parts of the article are organized as follows.

Section 2 reviews the related literature. Section 3 presents some mathematical formulations for the problem. Sections 4 and 5 describe the branch-and-price and the hybrid genetic algorithm, respectively. Section 6 reports the experimental analyses, and Section 7 concludes.

2 Related literature and subproblems

Vehicle routing problems with profits are the subject of an extensive literature, surveyed in [3, 23, 60, 63]. Generally, one seeks to jointly minimize cost and maximize customer’s service levels, two objectives which are conflicting when the profits do not render all deliveries profitable, or in the presence of additional vehicle constraints (e.g., distance or capacity limits). The related literature can thus be classified according to three fundamental solution techniques for this bi-objective problem:

• I)

Weighted sum – minimizing a weighted difference of costs and service levels;

• II)

Constraints on cost/distance – maximizing service levels subject to distance constraints (independently for each route);

• III)

Constraints on service levels – minimizing distance subject to service level constraints.

Table 1 presents an overview of these main classes of methods, for the single and multi-vehicle routing problems with profits (variants of the TSP and VRP, respectively). Objectives I and II are the subjects of a vast literature, which includes a very large variety of exact methods and metaheuristics. The TOP, in particular, seeks the maximization of service levels subject to distance constraints, and has been studied in dozens of articles in the past decade (see [1, 2, 13, 19, 22, 37, 38, 40, 44, 54, 59, 66], among others). In contrast, a larger methodological gap remains for objective III, which aims to minimize travel cost subject to service level constraints. The multi-vehicle case of this objective has been mostly studied in the context of hot strip mill scheduling for steel production [68, 69]. More recently, [43] investigated a special case with one group of customers (all of them) and one service level constraint, usually referred to as the prize-collecting VRP (PCVRP). The authors proposed a self-adaptive VNS and reported computational experiments on new instances. These instances, however, are now unavailable. The VRP-SL generalizes this problem by introducing multiple groups and associating profits to deliveries. By doing so, the problem remains concise and simple to define, but generalizes many VRP classes:

• –

The PCTSP and PCVRP correspond to VRP-SL instances with a single group.

• –

The CPTP, PTP and VRPPFCC all correspond to instances with 0% service levels.

• –

The CVRP corresponds to instances with 100% service levels.

• –

The generalized VRP (GVRP) can be reduced to the VRP-SL by imposing a small service level for each group, hence forcing at least one delivery.

• –

Any instance of the periodic vehicle routing problem (PVRP) such that, for each customer $i$, any combination of $f_{i}$ different days in a set $D_{i}$ of available days is feasible, can be transformed into a VRP-SL instance with $\sum_{i}|D_{i}|$ customers. This is done by duplicating each customer $i$ into as many vertices as possible visit days, defining a group for the resulting visits with a service level $\alpha_{i}=f_{i}/|D_{i}|$, and setting a large cost $M$ for the edges that link two vertices from different days. This reduction follows the same principles as the reduction from PVRP instances with frequency $f_{i}=1$ into GVRP instances, discussed in [5], but is more general as it allows to deal with frequency values greater than 1.

Finally, since the name prize-collecting VRP is not consistently used in the literature and far from self-explanatory, we opted for the name VRP with service levels (VRP-SL) for the proposed problem with several groups, which is clearer and more distinctive.

3 Mathematical Formulations

Compact formulation. We first propose a mixed integer linear formulation of the VRP-SL. Our mathematical model is based on the two-commodity flow formulation of [7], which has already been extended to several closely related VRP variants in [30, 31]. The graph ${\cal G}$ is first extended into $\overline{{\cal G}}=(\overline{{\cal V}},\overline{{\cal E}})$ by adding a new vertex $v_{n+1}$, representing a copy of the depot. We thus define $\overline{{\cal V}}={\cal V}\cup\{v_{n+1}\}$, ${\cal V}^{\prime}=\overline{{\cal V}}\setminus\{v_{0},v_{n+1}\}$, $\overline{{\cal E}}={\cal E}\cup\{(v_{i},v_{n+1}),v_{i}\in{\cal V}^{\prime}\}$, and $d_{i,n+1}=d_{0,i}$ for all $v_{i}\in{\cal V}^{\prime}$.

For each edge $(v_{i},v_{j})\in\overline{{\cal E}}$, we define a binary variable $x_{ij}$, set to $1$ if and only if a vehicle travels on this edge, as well as two flow variables $f_{ij}$ and $f_{ji}$. When a vehicle travels from $v_{i}$ to $v_{j}$, the flow $f_{ij}$ represents the current load in the vehicle, and the flow $f_{ji}$ represents the residual capacity of the vehicle ($f_{ji}=Q-f_{ij}$). Finally, $y_{i}$ is a binary variable which is set to $1$ if and only if $v_{i}\in{\cal V}\setminus\{v_{0}\}$ is serviced. The VRP-SL can be formulated as:

[TABLE]

The objective of Equation (1) aims to minimize the sum of transportation costs and lost profits.

With this objective, the optimal solution value is always non-negative (no negative terms). Moreover, the optimal CVRP solution is a feasible solution of this model, with a cost greater or equal to the VRP-SL optimum. Constraints (2) impose the service levels for each group. Constraints (3) ensure that each vertex of ${\cal V}\setminus\{v_{0}\}$ is visited at most once. These constraints also connect the path variables ($x_{i}$) to the customer selection variables ($y_{i}$) in order to evaluate the profits. Constraints (4)–(7) define a feasible two-commodity flow between the source $v_{0}$ and the sink $v_{n+1}$. Specifically, Constraints (4) state that the inflow minus the outflow at each vertex $v_{i}\in{\cal V}^{\prime}$ is equal to $2q_{i}$ if $v_{i}$ is used, and [math] otherwise. The outflow at the source vertex $v_{0}$, computed in Constraint (5), is set to the total demand of the vertices that are serviced in the solution, and the outflow at the sink $v_{n+1}$, calculated in Constraint (6), corresponds to the total capacity of the vehicle fleet. Finally, Constraints (7) establish the link between the flows $f_{ij}$ and $f_{ji}$, and Constraint (11) sets a bound on the number of vehicles.

The linear relaxation of the VRP-SL can be strengthened via some simple valid inequalities. If $s_{i}$ is integer for all $v_{i}\in{\cal V}\setminus\{v_{0}\}$, then the service level constraints (2) can be transformed into:

[TABLE]

Moreover, the following flow inequalities from [7, Equation 64] are used:

[TABLE]

The first inequality explicitly forces $f_{ij}$ to contain the delivery quantity $q_{j}$ for the next customer, and the second inequality forces the residual capacity $f_{ji}$ to be greater than the delivery quantity $q_{i}$ of the previous customer.

Finally, for each group $k\in\{1,\dots,K\}$, we define $Z_{k}$ as the minimum load quantity which allows to satisfy the service level constraint:

[TABLE]

These values can be evaluated in pseudo-polynomial time. Then, the following capacity cut is valid for each group $k$:

[TABLE]

Note that we also investigated in preliminary experiments a similar formulation based on one-commodity flows [27] for directed graphs, hence associating two variables $x_{ij}$ and $x_{ji}$ for each edge ($v_{i},v_{j}$). Although both formulations should produce similar bounds [41], the suggested two-commodity flow formulation of the VRP-SL led to overall better results in our context, possibly due to the smaller number of binary variables.

Set Partitioning Formulation. We now present a set-partitioning based formulation of the VRP-SL, which corresponds to a Dantzig-Wolfe decomposition of the model of Equations (1)–(11). Similar formulations have been successfully used in the VRP literature in the past years to obtain stronger linear relaxations [6, 17, 26, 47]. The drawback of this formulation comes from its exponential number of variables, which must be tackled using a column generation algorithm.

Let $\Omega$ be the set of all feasible routes for the problem. A route $r\in\Omega$ is a closed walk starting and ending at the depot, visiting a set of customers only once and respecting the vehicle’s capacity (this definition will be extended afterwards to allow visiting a customer more than once). Then for each route, we define a binary variable $\lambda_{r}$ indicating whether the route $r$ is used in the solution. The resulting formulation is as follows:

[TABLE]

The objective function (16) minimizes the total cost of the active routes plus the lost profits. Constraints (17) are the same as Constraints (2). In Constraints (18), $a^{r}_{i}$ (boolean) represents the number of times route $r$ visits customer $i$. Then each constraint forces one of the routes which visits customer $i$ to be active if variable $y_{i}$ is set to $1$. Constraints (19) limit the number of active routes to the number of available vehicles.

4 Branch-and-price algorithm

The set partitioning formulation has an exponential number of variables. To circumvent this issue, we use a column generation algorithm which starts with no variables and solves a pricing sub-problem to generate new ones. For the formulation presented in the last section, the pricing sub-problem is an Elementary Shortest Path Problem with Resource Constraints (ESPPRC). Ideally, one would want to price only elementary routes, but since the ESPPRC is known to be strongly NP-hard [21], we solve a pseudo-polynomial relaxation, the Shortest Path Problem with Resource Constraints (SPPRC) [16], allowing the routes to visit a given customer more than once. For both problems, the objective is to find a route with negative reduced cost. Given the dual variables $\beta_{i}$ and $\gamma$ associated to Constraints (18) and (19), a constant $b^{r}_{ij}$ representing the number of times route $r$ traverses edge $(v_{i},v_{j})$ and setting $\beta_{0}=\gamma$, the original reduced cost of a route and its reformulation as a function of the edges are stated in Equations (22) and (23):

[TABLE]

When relaxing the ESPPRC into the SPPRC, the bounds obtained by column generation deteriorate in most cases. Some alternative relaxations have thus been proposed to find a better balance between efficiency and solution quality [16, 33, 46]. We use the ng-route relaxation [8], which has been successfully applied to multiple VRP variants in the last years. For each vertex $i$, a set $NG_{i}\subseteq{\cal V}^{\prime}$ is defined to represent the “memory” of $i$. These $NG_{i}$ sets usually contain a subset of vertices closest from $i$. The pricing sub-problem is solved by a forward dynamic programing algorithm, which maintains a memory of past visits to prohibit some vertices, but uses the $NG_{i}$ sets to reduce this memory size and thus the size of the state space. During the algorithm, each path $P$ has an associated label $\mathcal{L}(P)$ containing the last customer visited $v(P)$, the total reduced cost $\overline{c}(P)$, the current load $q(P)$ and the customers which have been visited and remembered $\Pi(P)$. When extending a path from customer $v(P)$ to customer $v_{j}$, the extension is only allowed if $j\notin\Pi(P)$, and the label for the new path $P^{\prime}$ can be obtained as:

[TABLE]

Lastly, we use a dominance rule to fathom labels which cannot lead to an optimal solution. A label $\mathcal{L}(P_{1})$ dominates a label $\mathcal{L}(P_{2})$ if

[TABLE]

Other techniques can be used to further improve the overall column generation efficiency. We use two approaches. The first one is a simple heuristic pricing. This algorithm stores only the label with best reduced cost for each customer $v_{i}$ and load $q$ during the dynamic programming. The second approach is dual stabilization [50]. We use a parameter $\alpha\in[0,1[$, as shown in Equations (25)–(26), to avoid a large variation in the values of the dual variables between two iterations of the column generation. Our algorithm starts with $\alpha=0.9$ and reduces the value of $\alpha$ by $0.1$ each time the pricing sub-problem returns an invalid route.

[TABLE]

Finally, to improve the bounds obtained by the column generation algorithm, we embedded it into a branch-and-bound (B&B) procedure. At each node of the branch-and-bound tree, the column generation algorithm is called to obtain the solution of the linear relaxation, taking into account possible fixed variables due to branching. This algorithm framework is usually called branch-and-price (B&P). A key difference with classic B&B algorithms relates to the fact that it is not possible, in our context, to branch on the $\lambda_{r}$ variables, since a fixing of $\lambda_{r}=0$ would result in repricing the variable. To overcome this difficulty, the B&P branches on the original $x_{ij}$ variables as well as the $y_{i}$ variables, choosing at each node on the tree the most fractional variable to branch, but always giving priority to $y_{i}$ variables. Fixing a $y_{i}$ variable is straightforward. The $x_{ij}$ variables, however, are not explicitly present in the formulation, and thus we add Equation (27) to the set partitioning formulation for each fixed edge. In this equation, $b^{r}_{ij}$ represents the number of times route $r$ traverses edge $(v_{i},v_{j})$, and $\overline{x}_{ij}$ is the value which should be fixed for $x_{ij}$.

[TABLE]

The presence of Equation (27) introduces a new dual variable $\rho_{ij}$ which must be considered when computing the reduced costs, thus changing Equation (23) into Equation (28).

[TABLE]

This combination of techniques leads to an efficient exact method, whose performance will be analyzed in Section 6.

5 Population-based metaheuristic

The VRP-SL is known to be NP-hard as it generalizes the CVRP, and even sophisticated exact approaches can only solve small- and medium-scale problem instances within a reasonable CPU time. To fill this methodological gap and to solve the larger instances which arise in practice, we introduce a dedicated hybrid genetic search with advanced diversity control.

5.1 General structure of the method

The method, illustrated in Algorithm 1, uses the same resolution strategy as the unified hybrid genetic search (UHGS) of [62, 64]. Starting from an initial population, it iteratively selects two parents to generate an offspring individual via a crossover operator. This offspring is improved by means of a local search procedure and inserted in the population. This sequence of operations is performed until the termination of the method, once $\mathit{It}_{\text{NI}}$ successive iterations without improvement have been performed.

The method exploits penalized infeasible solutions in the population, which is divided into two sub-populations of feasible and infeasible individuals. Whenever one sub-population reaches a maximum size, a number of individuals are eliminated to retain the best solutions. A repair procedure is also applied on infeasible solutions to restore feasibility and generate additional feasible individuals. Finally, the approach uses an adaptive diversity management, which has been shown to be particularly successful when solving VRPs [62]. Parents and survivors selections are driven by two criteria, cost and diversity contribution, rather than solely on cost as in traditional GAs, and additional diversification phases are implemented after every $\mathit{It}_{\textsc{div}}$ iterations without improvement to provide new solution characteristics to the population.

Finally, the proposed algorithm also significantly differs from UHGS in the definition of its basic building blocks: solution representation, crossover, local search moves, distance measure between individuals, and penalties allowed. These components are described below.

5.2 Solution Representation and Evaluation

In the proposed HGS, each individual in the population is represented by two chromosomes: a service level chromosome, which gives the current service level of each group $k\in\{1,\dots,K\}$, and the giant-tour chromosome, which provides a permutation of visits for the serviced customers, without occurrences of the depot. This solution representation is illustrated in Figure 1. It is incomplete in the sense that some additional information, the locations of the visits to the depot, is needed to perform cost evaluations. Still, this information can be quickly recovered by means of a dynamic-programming algorithm, called Split, which optimally subdivides the giant-tour chromosome into separate routes. In the proposed HGS, the Split algorithm strictly respects the sequence of visits, i.e., it cannot select or exclude customer visits and modify the service levels. This is the classical context of application of the algorithm of [9] and [52], which reduces the splitting problem into the search of a shortest path in an acyclic graph in which each arc represents a possible route, i.e., a sequence of consecutive visits in the giant tour, connected to the depot. The reader is referred to [52] and [61] for a detailed description of the Split algorithm as well as an efficient $\mathcal{O}(n)$ implementation.

We now describe the cost function used for route and solution evaluations. As observed in our computational experiments and in [65], the exploration of penalized infeasible solutions during the search has a significantly positive impact on solution quality. Thus, the cost $\phi(r)$ of a route $r$ involves the distance, the total profit associated to the customers which are visited, and a possible excess of load in a route multiplied by a penalty factor $w^{Q}$. Let $\phi^{D}(r)=\sum_{i=1}^{|r|-1}d_{r(i),r(i+1)}$, $\phi^{Q}(r)=\sum_{i=1}^{|r|}q_{r(i)}$ and $\phi^{P}(r)=\sum_{i=1}^{|r|}p_{r(i)}$ be, respectively, the total distance, load and profit collected in route $r$, then the route cost is defined as:

[TABLE]

We also explore infeasible solutions with respect to service level constraints. These linking constraints involve the whole solution, and thus the penalty is defined at the level of the solution evaluation rather than the route cost. The penalized cost $\phi^{\textsc{cost}}(S)$ of a solution $S$, described as a set of routes $r\in S$, can thus be evaluated as:

[TABLE]

In this equation, $\phi_{P}=\sum_{i=1}^{n}p_{i}$ is the total profit of all customers (a constant), $\phi^{k}(S)$ is the weight ratio of group $k$ in solution $S$, and $w^{S}_{k}$ for $k\in\{1,\ldots,K\}$ are the penalty coefficients associated to service level violations. These penalty coefficients are automatically adjusted by the method during the search, as explained in Section 5.4.

5.3 Generation of New Individuals

Each new solution is generated by a successive application of the Selection, Crossover, and Education operators, followed by a possible Repair.

Selection and Crossover

To generate a new solution, the algorithm first selects two parents $P_{1}$ and $P_{2}$ in the population via a binary tournament based on the biased fitness measure described in Section 5.4. A new offspring solution $C$ is then obtained by crossover of $P_{1}$ and $P_{2}$. For this purpose, we propose an adapted order crossover (AOX), which extends the well-known order crossover with the ability to transmit customer selection and visit sequence decision from both parents. This crossover operator is illustrated in Figure 2.

In a first step, the crossover inherits the service level information from both parents. This is done by crossing the service level chromosomes of both parents using an extended intermediate recombination [48], i.e., a target weight ratio $\alpha^{T}_{k}(C)$ is randomly chosen between $\alpha_{k}(P_{1})$ and $\alpha_{k}(P_{2})$ for each group $k$, where $\alpha_{k}(P)$ is the weight ratio of group $k$ in individual $P$.

Then, in a second step, the giant-tour chromosome of the child $C$ is initialized with the longest size among both parent and inherits, as in the order crossover (OX), a fragment of $P_{1}$.

In the third and final step, the giant-tour chromosome of $C$ is completed by sweeping circularly the deliveries of $P_{2}$ and inheriting them, starting one index after the end of the fragment from $P_{1}$. Each insertion of a visit $i$ of a group $k$ is done under the condition that $i$ does not already exist in $C$, and that the target service level $\alpha^{T}_{k}(C)$ has not yet been reached. To complete the representation, the service level chromosome of the child is finally derived from the giant-tour chromosome.

Education

The goal of the crossover operator was to generate new solutions which inherit common characteristics from both parents while introducing a significant level of randomness. As such, the crossover operator is not the main force which drives solution improvement, this role being assumed by a subsequent local search-based education procedure.

The local search (LS) is applied on the complete solution representation, including the visits to the depot. Therefore, the Split algorithm has to be run beforehand. The LS uses the same classical vehicle routing neighborhoods as in [62]: 2-Opt, 2-Opt*, Swap, Relocate as well as generalized Swap and Relocate involving two consecutive nodes, and limited to moves between close services. These classical neighborhoods only involve services which are present in the current solution. To also optimize the decision subset related to customer selections in the LS, we include three additional neighborhoods:

• –

Remove: If $u$ is a visited customer, then remove $u$ from the solution.

• –

Add: If $u$ is a visited customer and $v$ is a non-visited customer, add $v$ after $u$.

• –

Replace: If $u$ is a visited customer and $v$ is a non-visited customer, replace $u$ by $v$.

All neighborhoods are explored in a random order with a first improvement move acceptance policy. The LS stops when no improving move can be found in the entire neighborhood, and the resulting solution is converted back into a giant-tour and service level individual representation, which is inserted in the population.

Repair

Finally, it is possible for a solution to remain infeasible after education. When this situation happens, a Repair operator is applied with $50\%$ chance. As in [62], this operator simply consists in running the local search with higher penalty values, first by a factor of 10, and then 100, with the aim of converging towards a feasible solution.

5.4 Population Management

The population is formed of two sub-populations, designed to host feasible and infeasible individuals, respectively. The algorithm starts by generating $4\times\mu$ random initial solutions. Each sub-population is then managed to contain at least $\mu$ and at most $\mu+\lambda$ individuals. Whenever a population reaches its maximum size, $\lambda$ individuals are eliminated to produce the next generation. This is done by iteratively selecting out either a clone solution, with a distance of [math] to another solution, or the worst solution according to a biased fitness function $\phi^{\textsc{bias}}(S)$ when there are no more clones. The biased fitness measure evaluates every solution $S$ based on its cost $\phi^{\textsc{cost}}$ (Equation 30) and its contribution to the sub-population diversity, defined as:

[TABLE]

where $N_{P}(S)$ is the set of the $n^{\textsc{close}}$ individuals in the sub-population $P$ closest to $S$ with respect to a distance measure $\delta(S_{1},S_{2})$. Our distance measure counts the percentage of common edges between two solutions. Let $E(S_{1})$ and $E(S_{2})$ be the set of edges used in the solutions $S_{1}$ and $S_{2}$, then the distance is expressed as:

[TABLE]

Finally, let $n^{\textsc{elite}}$ be a parameter controlling how elitist $\phi^{\textsc{bias}}(S)$ is, and let $R(S,P,f)$ be an application which returns the rank of an individual $S$ in the sub-population $P$ relatively to a measure $f$. Then, the biased fitness of $S$ in $P$ is evaluated as in UHGS [64]:

[TABLE]

Each distance measure can be evaluated in $\mathcal{O}(n)$ time. As such, when a new individual enters the population, its distance from every other solution can be computed in $\mathcal{O}(n|P|)$ time, and the computational complexity of updating the biased fitness measures is $\mathcal{O}(n|P|+|P|\log|P|)$.

Adaptive Penalty Coefficients

The exploration of infeasible solutions contributes positively to the search if the ratio of feasible solutions is adequately controlled. To that extent, the penalty parameters are adapted to achieve a ratio of feasible solutions within a predefined interval $[\xi^{\textsc{min}},\xi^{\textsc{max}}]$.

We rely on $1+K$ penalty parameters for the VRP-SL: one for the capacity constraints, and $K$ penalties for the service level constraints, one for each group. Let $\xi^{c}$ be the ratio of feasible individuals with respect to a constraint $c$, measured among the last 100 individuals generated by local search. In order to drive the search towards feasible solutions, every 100 iterations we update the penalty coefficient of constraint $c$ using the following rule:

[TABLE]

Initially, all penalty parameters are set to $10$. We aim to obtain around 25% feasible solutions after LS. Since all constraints need to be respected to obtain a globally feasible solution, we used $\xi^{\textsc{min}}=0.15^{\frac{1}{1+K}}$ and $\xi^{\textsc{max}}=0.35^{\frac{1}{1+K}}$ to achieve this goal.

Diversification procedure

Finally, after each consecutive $It_{\textsc{div}}$ iterations without improvement of the best solution, we apply the same diversification procedure as in [62] to only keep the best $\mu/3$ individuals in each subpopulation and reintroduce new random initial solutions. This procedure complements the biased fitness function so as to avoid a premature convergence of the method due to the strong intensification of the LS.

6 Experimental Analyses

This section aims to

1. introduce a set of instances for the VRP-SL derived from classical vehicle routing instances,
2. evaluate the performance of the proposed methods on the VRP-SL instances,
3. evaluate the performance of the metaheuristic and the branch-and-price on classical problems generalized by the VRP-SL, namely the VRPPFCC and CPTP, and finally
4. examine the impact of some key parameters and design choices. All algorithms were coded in C/C++. The exact algorithms were run on a single thread of a 3.07-GHz Intel Xeon CPU, and the metaheuristic was run on a single thread of a 3.4-GHz Intel Core i7 CPU. We used CPLEX 12.7 for the resolution of the compact formulation and for the linear programs in the B&P algorithm.

Benchmark instances for the VRP-SL. To compensate for the unavailability of benchmark instances for the VRP-SL, we derived two sets of instances from classical CVRP and prize-collecting VRP test sets:

• –

The first set (S1) has been generated from a subset of 26 instances from [4], and includes between 31 and 80 vertices. The profit of each client has been set to $h\times q_{i}$, where $h$ is a random variable uniformly generated in the interval $[0.75,2.25]$.

• –

The second set (S2) has been derived from the 10 capacitated profitable tour instances of [1], with 51 to 200 vertices. The capacity of the vehicles has been set to 500, and the number of vehicles has been set to

[TABLE]

where $Q_{k}^{\textsc{min}}$ is the subset of services of group $k$ with smallest delivery quantity which allows to satisfy the service level requirements, and $Q_{k}^{\textsc{max}}$ is the sum of demands of all customers of this group. The original profits were multiplied by a factor of $0.5$ to obtain a good balance between profits and distance.

For each instance, we considered five configurations for the assignment of visits to groups: {1, 2R, 2C, 5R, 5C}. The first number corresponds to the number of groups, and the second letter, when applicable, corresponds to their distribution: random (R) or clustered (C). In all cases, the service levels for each group have been randomly selected in $\{0.45,0.55,0.65,0.75,0.85,0.95,1\}$, and the weight of each customer, reflecting its importance for the service level constraint, coincides with the demand quantity ($s_{i}$ = $q_{i}$). Overall, this leads to $26\times 5=130$ instances of set S1, and $10\times 5=50$ instances of set S2. For each instance, the distance values between customers are rounded to the nearest integer. For the sake of brevity, the results presented in the paper are aggregated per group of five instances. All instances and detailed results are available in the electronic companion of this paper, also available at https://w1.cirrelt.ca/~vidalt/en/VRP-resources.html.

6.1 Exact solutions and lower bounds

The two exact methods have been tested on each VRP-SL instance using a single core with a time limit of two hours. Since the VRP-SL instances were designed to be challenging for both exact and heuristic approaches, only a subset of the problems could be solved to optimality. To speed up the resolution, we used the best integer solution found by the metaheuristic as a warm start for both exact methods, hence limiting the size of the branch-and-bound tree.

Table 2 presents average results for each group of instances. The first columns list the characteristics of each group, followed by the minimal number of visits $n_{\textsc{min}}$ required to satisfy service level constraints, the value LB0 of the lower bound and the processing time T0 at the root node, the value LB of the best overall lower bound, the number of nodes in the search tree, the final integrality gap and the total CPU time. For each instance, the best lower bound is highlighted in boldface. We do not indicate the best upper bound found by each method, since no improvement was found over the initial value obtained by the metaheuristic.

As illustrated by the experiments, the branch-and-price algorithm outperformed the compact formulation-based method on most instances. More precisely, it performed better on all instances of the set S1 with the exception of B1 and B2. For the set S2, the compact formulation performed better for the first five groups of instances and worse for the remaining ones. As expected, the branch-and-price algorithm is very effective when the solution includes short routes, while the compact formulation performs better on instances with few vehicles. The total number of instances solved to optimality for configurations {1, 2R, 2C, 5R, 5C} were {10, 9, 9, 9, 8} for the compact formulation and {15, 9, 9, 8, 9} for the branch-and-price algorithm. These results show that the exact methods’ performances are rather insensitive to the distribution of nodes in $V$, as well as the number of groups $K$. However, the complexity generally increases with $n$, $m$ and $n_{\textsc{min}}$, as visualized in Figure 3.

6.2 Performance of the hybrid genetic algorithm

In this subsection, we report the results obtained with the proposed HGS on the VRP-SL instances. For each instance, the algorithm was run ten times with different seeds, using the same parameter setting and termination criterion as in [62] for the CVRP, that is $(n^{\textsc{elite}},\mu,\lambda)=(10,25,40)$, $It_{\textsc{max}}=2\times 10^{4}$ and $It_{\textsc{div}}=0.4\times It_{\textsc{max}}$. Moreover, to accelerate the convergence, the parameters governing the population size have been halved when dealing with problem instances containing 200 or more services. Table 3 reports, for each group of instances, the worst, average and best solution quality of HGS over 10 runs (Wor-10, Avg-10 and Best-10), the percentage gap between the average solutions and the best ones found (Gap${}_{\textsc{BKS}}$), the percentage gap between the average solutions and the best lower bounds found by the exact methods (Gap${}_{\textsc{LB}}$), the average CPU time in seconds, the best known solutions (BKS) and lower bounds (BKLB).

In the absence of results from previously published heuristics, we consider three main indicators of method performance: its ability to reach known optimal solutions, the stability of the solution quality over several runs, illustrated by the gap between the average solution quality and the BKS, and the deviation from the lower bounds produced by the mathematical programming algorithms. As observed in our experiments, HGS finds all of the 70 known optimal solutions on all test runs. The method also returns solutions of consistent high quality: it found for 25/36 groups of instances the value of the best known solution on all runs of all instances. The percentage gaps between the average and best known solutions are close to zero ($0.01\%$ overall), and the instances with fewer groups appear to be generally easier to solve. The gaps to the lower bounds are also small ($1.38\%$ in average), and thus the solutions of HGS are guaranteed to be close to the optima. This gap is more likely to be due to the quality of the lower bounds, since over all 2h-runs of the exact approaches not a single best solution of the metaheuristic was improved. Finally, the average CPU time never exceeds 93 seconds, the worst case being observed on problem p16 with 199 services.

6.3 Comparative analyses on key subproblems

As discussed in Section 2, the VRP-SL generalizes several important problem classes. Two such problems in particular, the VRPPFCC and CPTP, have been the focus of a wide literature, opening the way to some comparative performance analyses.

Experiments on the VRPPFCC

We rely on the two sets of instances from [12]. Set CE includes up to 199 customers, while Set G includes larger instances with up to 483 customers. For these instances, the convention is to compute all distances with double precision and report the final solution with two digits. Our heuristic method is compared to the two best current metaheuristics in the literature: the UGHS proposed by [66], which uses an exhaustive solution representation (with all visits) with a route-evaluation operator in charge of customer selections, and the recent AVNS of [32]. To the best of our knowledge there are no known exact results for this problem in the literature. The results of our methods are displayed in Table 4. The column BKS reports the best known solutions found in the literature before this paper. Avg-10 and Best-10 represent the average and best solutions found by the metaheuristic over ten runs. For each instance, the best heuristic method in highlighted in boldface. For the branch-and-price results, columns LB0 and LB follow the same convention as Table 2.

In these experiments, the proposed metaheuristic appears to outperform previous methods in terms of solution quality, with an average gap of 0.141%, in comparison to 0.445% for UHGS with the exhaustive solution representation, and 0.345% for AVNS. During these tests, eight new best known solutions (BKS) have been found, as underlined in the table. Finally, the average CPU time is markedly faster than previous approaches which were run on processors of a similar generation, with 340 seconds on average, compared to 1584 and 4656 seconds for the other methods. The proposed branch-and-price algorithm is the first in the literature to report optimal VRPPFCC solutions for three instances: CE-01, CE-06 and G-17. On the other hand, on six instances, the algorithm was not able to solve the root node within a time limit of two hours. This is due to their size (up to 480 customers), which can be considered very large for the current exact methods. Finally, the average integrality gap calculated over the tractable instances was 0.667%, confirming the good performance of the approach.

Experiments on the CPTP

For this problem, we compare the proposed HGS with the previous two best methods in the literature, the UHGS with exhaustive solution representation and the multi-start ILS of [66]. We also compare the proposed branch-and-price algorithm with the exact approach of [1], which generated, to this date, the best bounds for the CPTP. In that work, the authors proposed a branch-and-price algorithm using a q-route relaxation, and reported detailed results with and without a primal heuristic. We rely on the ten test cases of [2], each case being used to produce $12$ instances with a different fleet size and vehicle capacity, for a total of $120$ instances. Tables 5 and 6 present a summary of our computational experiments on these instances, using the same conventions as previously. Each line in the table corresponds to an average measure over a group of $12$ instances.

These instances are generally smaller, with $50$ to $199$ service locations, and thus the performance differences between state-of-the-art heuristics are less marked. Still, we observe that the proposed HGS obtains average solutions of similar or better quality (0.012% gap compared to 0.029% and 0.172% gap) in a fraction of the CPU time of previous algorithms ($14$ seconds in average, compared to $211$ or $205$ seconds). Three previous BKS were improved, leading to an average gap of $-0.002\%$ for the best solution quality of 10 runs. For the branch-and-price algorithm, we set the time limit to 3600 seconds in order to make a fair comparison with [1]. The proposed B&P found similar or better solutions for all instances tested. It improved the bounds for 37 instances, with an average improvement of 0.629%, and proved optimality for 103 instances, including ten new optimality certificates.

6.4 Sensitivity analyses

Finally, this section reports additional sensitivity analyses on the impact of key components of the proposed HGS. For this purpose, we compare the results of the standard method described in Section 5 against several alternative configurations obtained by modifying one operator, design choice, or group of parameters:

• EOX –

An edge-recombination crossover (EOX) is used. As described in [67], this crossover maintains, for each vertex $i$, an adjacency list of non-visited vertices adjacent to $i$ in at least one parent. After a random choice for the first vertex, EOX iteratively inserts the adjacent vertex with the shortest adjacency list. Ties are broken randomly, and a random vertex is chosen whenever the adjacency list is empty. As in our adaptation of OX, target service levels are inherited for the groups from the parents, and any vertex belonging to a group for which the target service level has already been attained is eliminated from the adjacency lists.

• No SL –

Service levels are not used to filter service insertions in the crossover.

• No INF –

All penalty coefficients are set to a large value to avoid infeasibility.

• No DIV –

Individual diversity contributions are not counted in the biased fitness.

• No Rep –

The Repair operator is not applied.

• Pop $\downarrow$ –

Smaller population: $(\mu^{elite},\mu,\lambda)=(4,12,20)$.

• Pop $\uparrow$ –

Larger population: $(\mu^{elite},\mu,\lambda)=(16,50,80)$.

• Feas $\uparrow$ –

50% feasible solutions as a target: $(\xi^{\textsc{min}},\xi^{\textsc{max}})=(0.4^{\frac{1}{1+K}},0.6^{\frac{1}{1+K}})$.

• Feas $\uparrow\uparrow$ –

75% feasible solutions as a target: $(\xi^{\textsc{min}},\xi^{\textsc{max}})=(0.65^{\frac{1}{1+K}},0.85^{\frac{1}{1+K}})$.

Each of these algorithm configurations was run 10 times on every benchmark instance for the VRP-SL, VRPPFCC and CPTP. Table 7 presents the average percentage gap, best percentage gap and time of each method for each set of instances.

From these experiments, it appears that the proposed method is a sort of “local optimum” in terms of design choice and parameter settings, in the sense that any change of its main operators and parameters impacts negatively the method performance. Still, some design choices have a much larger impact than others. In particular, using the EOX crossover operator strongly deteriorates the method performance for VRPPFCC instances, while speeding-up the resolution. Such a speed-up may be a symptom of premature convergence due, in this case, to the crossover. Both diversity management and infeasible-solution management contribute significantly to the performance of the method (configurations No INF and No DIV). This confirms the earlier observations of [62, 65]. Still, although the management of infeasible solutions is critical, deactivating the repair operator or changing the target level of feasible solutions has little impact (configurations No Rep, Feas $\uparrow$ and Feas $\uparrow\uparrow$). The control of the service levels in the crossover has a beneficial effect on performance for the CPTP (configuration No SL). Finally, HGS is relatively insensible to reasonable changes of population size (configurations Pop $\downarrow$ and Pop $\uparrow$).

7 Conclusions

In this article, we have introduced the VRP-SL, an important VRP variant arising in collaborative logistics operations, which aims to take into account the requirements of various partners via service level constraints on groups of deliveries. To establish a basis for further research, we introduced a first set of benchmark instances, a compact mathematical formulation, a branch-and-price algorithm and a first effective hybrid genetic search. The service level constraints tend to make the selection of services more complex, and thus new problem-tailored search operators, solution representation, crossover, LS moves, and penalty management strategies were introduced in HGS. Thanks to these elements, the proposed heuristic finds all known optimal solutions for the VRP-SL, and outperforms previous algorithms for two important special cases, the VRPPFCC and the CPTP, which have been intensively studied in past literature. Finally, the proposed branch-and-price algorithm was able to produce tight bounds for all problems and new optimality certificates for 63 instances, outperforming the existing exact approaches.

The research perspectives are numerous. First, the new algorithms can still be improved via the addition of new families of cuts, better relaxations, new neighborhoods and other heuristic strategies. Moreover, the VRP-SL is only a simplification of an intricate real-life application, and the assumptions about the time constraints, the dynamics and stochasticity of the problem were voluntarily simplified to allow for reproducibility. Guaranteeing, in a stochastic and on-line context, the satisfaction of contractual service levels for prize-collecting problems is an important challenge, as the violation of such obligations can lead to large penalties or lost contracts. To circumvent this risk while mitigating the cost of robustness, it is possible to search for robust solutions on a larger planning horizon, and consider alternative transportation modes as a recourse. These aspects, and the interactions between them, will be considered in future works.

8 Acknowledgments

This research is partially supported by CNPq and CAPES in Brazil, and the National Foundation for Science and Technology Development (NAFOSTED) in Vietnam.