Ordered p-median problems with neighborhoods

TL;DR

This paper introduces a new variant of the p-median problem where facility locations are uncertain within neighborhoods, proposing formulations and heuristics to optimize facility placement and assignment at minimal cost.

Contribution

It presents a novel p-median variant with neighborhood uncertainty, along with mixed-integer nonlinear formulations and heuristic methods for solution.

Findings

Developed multiple formulations for the problem.

Designed two heuristic algorithms for large instances.

Conducted extensive computational experiments demonstrating effectiveness.

Abstract

In this paper, we introduce a new variant of the $p$-median facility location problem in which it is assumed that the exact location of the potential facilities is unknown. Instead, each of the facilities must be located in a region around their initially assigned location (the neighborhood). In this problem, two main decisions have to be made simultaneously: the determination of the potential facilities that must be open to serve the demands of the customers and the location of the open facilities in their neighborhoods, at global minimum cost. We present several mixed integer non-linear programming formulations for a wide family of objective functions which are common in Location Analysis: ordered median functions. We also develop two math-heuristic approaches for solving the problem. We report the results of extensive computational experiments.

Tables

Table 1. Table 1. Integrality gaps for the formulations of DOMPN.

			Sc.1			Sc. 2			Sc. 3			Sc. 4
$n$	Pr.	$p$	${\mathrm{IG}}_{3I}$	${\mathrm{IG}}_{2I}$	${\mathrm{IG}}_{BEP}$	${\mathrm{IG}}_{3I}$	${\mathrm{IG}}_{2I}$	${\mathrm{IG}}_{BEP}$	${\mathrm{IG}}_{3I}$	${\mathrm{IG}}_{2I}$	${\mathrm{IG}}_{BEP}$	${\mathrm{IG}}_{3I}$	${\mathrm{IG}}_{2I}$	${\mathrm{IG}}_{BEP}$
10	C	2	3204.62%	442.38%	631.91%	3131.74%	605.50%	569.00%	998.10%	644.37%	383.48%	1796.27%	538.67%	371.76%
		3	1942.55%	483.36%	505.83%	2013.89%	561.03%	466.41%	582.44%	444.69%	296.41%	874.05%	418.28%	297.98%
		4	1151.43%	448.64%	408.55%	1289.11%	472.21%	387.08%	358.07%	323.78%	234.56%	458.06%	311.99%	233.18%
		5	773.88%	404.78%	340.68%	740.07%	404.48%	328.58%	258.86%	249.39%	195.38%	310.55%	247.76%	194.58%
	D	2	10294.86%	381.68%	629.23%	9961.19%	1031.13%	613.74%	3029.82%	1952.73%	528.34%	5391.17%	1604.94%	503.67%
		3	5303.63%	610.47%	510.75%	5410.29%	1315.60%	493.90%	1430.14%	1089.56%	399.89%	2246.47%	1062.61%	395.83%
		4	2780.29%	742.20%	430.02%	3047.22%	1074.29%	420.68%	752.36%	679.40%	314.18%	995.14%	671.09%	313.73%
		5	1615.81%	810.04%	369.91%	1497.97%	791.90%	355.09%	449.53%	432.78%	249.28%	555.35%	438.68%	250.08%
	K	2	12847.25%	521.31%	797.93%	12350.57%	1335.61%	767.53%	3705.93%	2390.52%	652.89%	6558.72%	1953.10%	606.63%
		3	7154.07%	860.28%	697.79%	7285.99%	1777.47%	651.89%	1880.91%	1433.23%	509.19%	2898.40%	1369.72%	475.52%
		4	4010.09%	1098.13%	609.63%	4333.66%	1524.65%	551.13%	1032.16%	931.75%	399.41%	1336.29%	900.58%	378.65%
		5	2439.10%	1219.91%	527.88%	2235.87%	1178.93%	463.39%	634.07%	610.26%	319.75%	759.07%	599.24%	297.79%
	M	2	17160.68%	371.06%	638.42%	16432.60%	1343.43%	630.91%	4947.70%	3187.44%	580.35%	8594.07%	2554.76%	551.56%
		3	8493.25%	719.43%	516.82%	8610.52%	1989.26%	510.23%	2217.45%	1688.09%	445.90%	3436.59%	1622.77%	433.57%
		4	4362.90%	1005.55%	448.08%	4716.87%	1659.11%	440.30%	1117.82%	1008.85%	355.60%	1440.78%	970.77%	345.31%
		5	2439.10%	1209.32%	393.95%	2235.87%	1178.93%	375.76%	634.07%	610.26%	284.65%	759.07%	599.24%	278.27%

Table 2. Table 2. Percentage of instances in which z 2 ​ I R ≥ z B ​ E ​ P R superscript subscript 𝑧 2 𝐼 𝑅 superscript subscript 𝑧 𝐵 𝐸 𝑃 𝑅 z_{2I}^{R}\geq z_{BEP}^{R} .

$n$	$p$	Sc. 1	Sc. 2	Sc. 3	Sc. 4
10	2	90.00%	32.50%	0.00%	0.00%
	3	50.00%	0.00%	0.00%	0.00%
	4	5.00%	0.00%	0.00%	0.00%
	5	0.00%	0.00%	0.00%	0.00%

Table 3. Table 3. Average results of comparing the formulations for the planar case for the n = 10 𝑛 10 n=10 instances.

			(${\mathrm{OMPN}}_{3I}$)			(${\mathrm{OMPN}}_{2I}$)			(${\mathrm{OMPN}}_{OT}$)			(${\mathrm{OMPN}}_{BEP}$)
Sc.	Pr.	$p$	time	#nodes	%gapR	time	#nodes	%gapR	time	#nodes	%gapR	time	#nodes	%gapR
1	C	2	217.19	118505.20	96.71%	26.73	51583.60	85.01%	0.39	79.60	32.97%	0.21	56.20	33.88%
		3	766.09	391417.00	93.85%	139.82	268504.80	84.06%	0.37	159.40	38.90%	0.26	103.20	36.97%
		4	3375.43	872945.80	88.82%	1088.98	1624482.00	78.99%	0.31	180.20	39.78%	0.19	54.80	36.34%
		5	3527.08	1109172.20	83.98%	2225.30	4064066.40	75.20%	0.29	217.80	45.23%	0.18	87.40	39.07%
	D	2	78.49	17176.80	98.90%	20.82	24835.40	86.82%	0.35	55.40	9.17%	0.26	34.00	4.49%
		3	142.18	15663.80	97.54%	21.56	30214.20	93.48%	0.35	6.60	10.01%	0.22	0.20	0.00%
		4	68.88	5336.60	95.00%	22.76	40081.00	90.40%	0.34	4.80	8.86%	0.15	0.20	0.00%
		5	79.83	10692.00	91.43%	147.91	304842.40	86.54%	0.43	45.00	15.07%	0.27	18.60	8.26%
	K	2	291.01	136621.00	99.14%	113.00	184032.00	90.47%	0.28	69.80	26.14%	0.34	88.40	23.22%
		3	2239.44	961277.00	98.20%	457.96	580536.00	95.25%	0.32	135.00	31.15%	0.30	103.40	26.99%
		4	3079.66	745095.80	96.55%	1538.17	1946764.20	93.38%	0.36	169.20	37.92%	0.31	123.40	32.72%
		5	613.05	246909.80	94.09%	4081.56	3501105.00	90.72%	0.38	259.40	42.68%	0.34	183.20	36.29%
	M	2	45.98	9877.60	99.32%	7.32	7234.40	86.93%	0.39	26.60	8.19%	0.26	5.20	2.21%
		3	35.95	6140.80	98.43%	10.44	6924.00	95.83%	0.23	6.80	2.43%	0.16	1.60	0.71%
		4	32.07	6413.60	96.76%	3.45	4077.80	93.77%	0.24	26.60	7.82%	0.15	0.20	0.00%
		5	24.75	4168.00	94.09%	4.42	4489.20	90.72%	0.27	9.20	5.37%	0.14	0.40	0.00%
2	C	2	185.79	243222.20	97.06%	35.65	160015.80	87.11%	0.18	100.40	33.21%	0.12	62.20	33.54%
		3	3368.98	2751134.40	94.71%	157.57	564670.20	82.55%	0.18	203.40	38.36%	0.15	165.00	31.62%
		4	3056.51	152278.20	90.92%	557.25	1636689.80	77.63%	0.18	255.40	38.03%	0.12	114.20	28.83%
		5	312.58	319355.20	82.92%	334.03	1082684.60	73.49%	0.18	234.00	39.00%	0.12	140.20	26.26%
	D	2	40.02	20172.60	99.00%	8.64	23279.00	95.43%	0.13	38.00	13.76%	0.14	31.80	8.55%
		3	76.92	30308.20	97.87%	11.71	34052.20	92.74%	0.13	22.00	10.90%	0.09	12.40	7.57%
		4	93.42	15935.20	95.85%	10.59	39267.20	89.51%	0.14	11.80	7.97%	0.13	0.60	4.02%
		5	66.53	15954.20	90.74%	14.53	70825.60	85.35%	0.16	40.00	14.52%	0.13	0.60	2.12%
	K	2	257.25	239132.80	99.21%	94.28	289128.40	96.41%	0.14	101.80	27.71%	0.16	105.80	25.70%
		3	1804.51	1410405.00	98.44%	320.74	836508.40	94.71%	0.13	116.00	32.21%	0.21	177.80	29.38%
		4	2053.02	1469281.00	97.08%	724.58	1705800.00	92.64%	0.17	177.00	37.55%	0.19	185.00	34.11%
		5	833.50	774058.20	93.53%	3164.10	6317876.00	89.75%	0.16	190.80	38.59%	0.18	245.80	34.51%
	M	2	37.72	18920.00	99.38%	5.01	9075.20	97.17%	0.16	42.60	11.33%	0.15	28.60	7.21%
		3	27.20	11529.40	98.63%	4.66	7610.40	95.35%	0.17	4.80	4.28%	0.09	1.00	1.78%
		4	24.24	11298.60	97.26%	4.86	12168.00	93.08%	0.19	59.00	4.15%	0.12	1.00	0.23%
		5	15.07	5040.80	93.53%	2.46	5522.40	89.75%	0.20	46.40	8.47%	0.12	5.60	3.31%
3	C	2	173.15	247260.00	88.93%	10.57	31894.00	83.40%	0.14	52.40	38.77%	0.10	59.80	29.83%
		3	1319.33	1565352.60	79.43%	87.54	246224.40	74.00%	0.13	97.00	39.18%	0.10	56.60	24.94%
		4	1143.93	996110.00	65.65%	128.58	429914.60	62.69%	0.14	117.60	34.17%	0.09	47.40	18.47%
		5	1035.39	1506863.00	53.49%	190.32	751796.20	52.24%	0.17	206.40	32.27%	0.14	150.00	17.85%
	D	2	48.22	36372.20	96.06%	4.12	9171.00	94.09%	0.13	62.20	33.80%	0.09	42.60	15.04%
		3	36.70	20929.40	90.86%	3.89	10671.20	88.43%	0.12	60.20	37.01%	0.12	18.80	8.69%
		4	340.46	506451.80	81.77%	4.29	11882.80	80.20%	0.10	74.20	40.74%	0.13	3.80	3.01%
		5	83.84	172707.40	69.92%	5.66	23412.20	69.08%	0.13	88.20	35.34%	0.16	22.20	4.15%
	K	2	558.29	466004.40	96.85%	131.97	270282.20	95.27%	0.12	49.80	41.50%	0.12	42.40	27.59%
		3	653.97	510418.80	93.17%	285.01	572239.00	91.36%	0.14	101.80	48.30%	0.15	91.60	29.63%
		4	309.64	307686.80	86.81%	1573.06	2854817.00	85.69%	0.16	163.60	54.60%	0.14	105.60	31.47%
		5	248.32	247774.40	77.51%	3732.59	1571394.80	76.88%	0.17	224.80	55.83%	0.16	147.20	28.04%
	M	2	50.20	37663.80	97.54%	3.95	6822.60	96.30%	0.16	80.40	37.33%	0.11	56.00	16.12%
		3	35.68	21078.00	93.93%	3.25	4576.40	92.31%	0.11	59.80	41.43%	0.15	17.40	9.80%
		4	22.94	16293.00	87.33%	1.94	3479.40	86.24%	0.14	58.80	48.34%	0.14	3.00	5.70%
		5	167.96	230957.40	77.51%	1.50	1669.40	76.88%	0.14	93.00	43.84%	0.14	11.60	3.94%
4	C	2	203.77	245152.40	94.31%	49.60	127865.40	82.70%	0.14	139.00	34.49%	0.12	121.20	32.57%
		3	1201.06	1114765.20	86.88%	529.54	728381.60	76.31%	0.17	249.80	33.46%	0.10	82.40	26.93%
		4	323.95	498328.20	73.29%	353.45	676140.80	64.98%	0.14	126.20	18.48%	0.08	63.40	8.65%
		5	1542.56	110782.40	62.29%	54.46	151331.20	56.15%	0.13	217.40	23.09%	0.07	51.40	6.05%
	D	2	423.35	210911.80	97.95%	13.32	29156.00	93.49%	0.19	166.80	32.28%	0.17	174.80	29.19%
		3	396.21	211169.80	94.31%	24.84	50853.00	89.21%	0.21	241.80	32.73%	0.20	243.20	28.73%
		4	213.08	121462.00	86.61%	18.59	45899.00	81.77%	0.20	277.40	26.03%	0.16	244.20	21.95%
		5	65.49	49704.20	76.53%	20.74	54514.80	72.06%	0.17	160.00	18.81%	0.14	168.20	14.17%
	K	2	700.32	457647.40	98.34%	349.44	456972.60	94.73%	0.16	158.60	40.93%	0.18	136.40	38.28%
		3	2865.25	1667538.20	95.59%	643.49	870765.40	91.64%	0.23	296.80	40.93%	0.24	248.40	37.75%
		4	2825.25	736735.60	90.02%	3361.02	4465105.80	86.42%	0.26	404.80	38.85%	0.22	304.20	35.83%
		5	1753.15	257832.20	82.12%	5823.20	2160439.20	78.71%	0.23	372.80	32.61%	0.20	350.80	29.26%
	M	2	576.80	281412.80	98.67%	6.96	12096.60	95.78%	0.17	150.40	38.42%	0.16	133.60	33.23%
		3	699.11	296994.00	96.20%	4.90	7174.60	92.80%	0.22	313.40	41.03%	0.17	272.60	35.29%
		4	170.92	109314.00	90.54%	6.66	12752.60	87.12%	0.19	306.20	36.70%	0.17	280.60	30.04%
		5	138.54	102782.60	82.12%	5.74	13638.80	78.71%	0.16	209.00	29.09%	0.20	433.00	22.64%

Table 4. Table 4. Average results of comparing the formulations for the 3-dimensional instances with n = 10 𝑛 10 n=10 .

			(${\mathrm{OMPN}}_{3I}$)			(${\mathrm{OMPN}}_{2I}$)			(${\mathrm{OMPN}}_{OT}$)			(${\mathrm{OMPN}}_{BEP}$)
Sc.	Pr.	$p$	time	#nodes	%gapR	time	#nodes	%gapR	time	#nodes	%gapR	time	#nodes	%gapR
1	C	2	473.84	368360.8	97.22%	35.74	52967	74.37%	0.24	82.2	24.03%	0.15	79.4	25.17%
		3	5763.65	845075.8	95.59%	1675.21	294038.4	79.75%	0.24	157.2	33.26%	0.15	126.6	35.43%
		4	5430.28	799358.6	92.80%	4792.30	494225	80.44%	0.26	216.8	38.54%	0.15	109	39.57%
		5	3842.67	611415.4	89.13%	4801.62	2054248.2	78.25%	0.23	203.2	41.80%	0.14	147	41.37%
	D	2	61.91	18930	99.19%	88.55	43946	68.97%	0.13	6.6	2.99%	0.15	2.6	0.91%
		3	371.31	47028.6	98.49%	199.28	68335	80.33%	0.18	29.6	6.27%	0.14	31.8	3.96%
		4	132.89	41835	97.18%	142.35	109411.8	87.92%	0.19	34.6	8.42%	0.13	16.6	5.93%
		5	128.02	28957	95.14%	192.70	260491.8	90.06%	0.24	79.8	8.15%	0.18	29.2	4.73%
	K	2	511.77	225792.8	99.33%	181.79	222769.2	76.49%	0.22	138.2	19.71%	0.22	146.8	17.53%
		3	4807.08	995386.6	98.87%	4377.40	1019791	86.45%	0.23	176	31.58%	0.26	235	29.14%
		4	5913.30	116629	98.05%	4861.38	2118727.8	92.08%	0.27	261.6	39.21%	0.27	226	36.34%
		5	4950.40	324797	96.85%	6245.90	2442139.8	93.58%	0.21	224.8	45.21%	0.23	308.8	41.78%
	M	2	57.47	18665.2	99.52%	14.19	13224	67.92%	0.21	32.4	3.85%	0.15	4.8	0.99%
		3	98.60	24726.8	99.08%	75.02	14259.6	80.26%	0.15	19.2	2.37%	0.11	36.2	1.60%
		4	42.06	18444	98.23%	9.60	11981	89.83%	0.17	14.2	2.58%	0.12	3	1.15%
		5	27.41	15402.4	96.85%	4.99	9265.4	93.43%	0.15	3.2	2.72%	0.14	3.8	0.22%
2	C	2	698.80	582674.6	97.05%	52.27	98390.6	84.39%	0.15	50.6	23.94%	0.13	67.6	25.04%
		3	6236.30	1034617.6	95.60%	1624.72	1452726	84.25%	0.18	123.2	32.22%	0.12	56.6	33.39%
		4	5647.40	696189.8	93.46%	3311.32	3658955.8	82.04%	0.20	176.4	36.36%	0.17	130.2	35.83%
		5	5778.93	20217.4	88.74%	5513.03	2033249.2	79.11%	0.20	179.8	38.00%	0.11	56.4	35.64%
	D	2	26.97	13871	99.14%	10.64	21167.8	91.26%	0.13	19	7.30%	0.16	10.6	5.64%
		3	26.09	16012.4	98.48%	7.30	14069.6	94.35%	0.12	28.4	7.39%	0.14	39	5.16%
		4	45.70	10844	97.38%	18.50	30920.4	92.60%	0.13	55.6	8.75%	0.13	38.4	7.08%
		5	23.17	4170.6	94.76%	11.40	40594.2	90.06%	0.16	88	9.53%	0.17	13.8	4.38%
	K	2	1785.21	167030.2	99.28%	255.95	299773.6	93.07%	0.13	41.2	20.37%	0.16	40.8	18.82%
		3	3734.82	368152.2	98.86%	2184.15	792684.6	95.79%	0.16	86.6	32.39%	0.16	78.2	30.66%
		4	2497.71	540841	98.16%	4728.84	1716593.4	94.82%	0.17	137.4	38.34%	0.17	121.8	36.34%
		5	3099.23	160480.8	96.57%	6860.23	3165215	93.49%	0.19	188.2	43.53%	0.18	163.8	41.23%
	M	2	22.84	12370.6	99.49%	5.46	7523.8	92.41%	0.14	13.4	6.93%	0.16	8.4	2.50%
		3	25.44	11016.8	99.06%	4.84	7306.8	96.49%	0.10	19.8	6.12%	0.16	17.2	3.92%
		4	16.44	7086.4	98.33%	2.74	4104.8	95.30%	0.14	19.2	6.04%	0.13	18.2	3.70%
		5	12.70	2887.4	96.57%	2.24	2401.8	93.49%	0.15	18.4	1.56%	0.10	3	0.71%
3	C	2	420.03	251532.6	90.68%	44.81	67897.8	84.83%	0.21	75	28.92%	0.12	63	22.81%
		3	4166.57	629476.2	85.05%	2144.77	627453.2	79.66%	0.24	161	36.97%	0.18	172.4	29.87%
		4	3496.18	360878.8	76.21%	2907.10	1019207.2	73.11%	0.30	273.6	37.35%	0.18	183	27.55%
		5	4390.44	566911.6	66.78%	5770.53	1238536.8	64.98%	0.28	307.2	37.85%	0.16	181.4	26.25%
	D	2	88.31	43064.4	97.12%	15.64	27226.4	95.31%	0.14	43.2	24.16%	0.15	37	10.38%
		3	257.27	72613.6	94.29%	42.11	46924.8	92.25%	0.25	191.6	27.63%	0.17	95.8	12.43%
		4	524.51	104272.8	89.46%	310.57	133618.4	88.08%	0.25	229.8	31.88%	0.21	150.6	11.61%
		5	284.34	89650	82.26%	161.15	228266	81.29%	0.30	324.4	34.84%	0.22	236.4	10.35%
	K	2	2833.08	429198.2	97.59%	273.07	318542.4	96.08%	0.18	70.4	31.86%	0.15	50.6	21.16%
		3	5614.46	389241.2	95.60%	2917.04	1364077.6	94.03%	0.28	219.6	43.42%	0.25	200.6	31.08%
		4	5881.77	60182	92.30%	6243.04	1470989.4	91.30%	0.37	342.4	50.21%	0.30	347.2	36.06%
		5	5833.66	40320.2	87.79%	7200.04	1	87.13%	0.33	357.4	55.10%	0.30	383.4	39.35%
	M	2	116.51	50768.6	98.27%	6.31	9523.2	97.18%	0.14	55.8	25.78%	0.14	28	10.32%
		3	174.75	52089.6	96.39%	50.50	19088.4	95.10%	0.23	196.4	27.91%	0.17	71	8.60%
		4	151.05	53351.2	93.05%	11.74	13976.8	92.15%	0.23	232.8	32.40%	0.14	94	8.53%
		5	74.87	32522.8	87.79%	9.07	12475	87.13%	0.26	281.4	37.69%	0.14	158	8.97%
4	C	2	535.61	315885.4	94.76%	64.32	98609.2	83.19%	0.17	101.2	26.22%	0.16	90.6	21.68%
		3	3061.28	171458.6	89.60%	1463.17	1489850.8	78.88%	0.19	132.8	25.60%	0.19	146.6	20.10%
		4	3104.27	151466.8	81.17%	1812.87	700664.4	73.08%	0.21	181.4	26.17%	0.26	213.8	19.90%
		5	3058.78	128605.2	71.34%	2007.26	496927.6	64.75%	0.18	175.4	22.72%	0.12	110	15.01%
	D	2	688.65	156903.4	98.35%	67.95	62430.4	94.61%	0.22	179.2	21.47%	0.26	221.8	19.05%
		3	1375.45	231147	96.31%	47.77	45659.8	92.40%	0.28	278	21.69%	0.25	258.8	18.76%
		4	683.20	149878.4	91.86%	206.62	98633	88.25%	0.31	354.6	21.45%	0.32	355	18.26%
		5	267.16	84209.6	85.28%	187.50	221145.2	81.69%	0.26	334.4	18.47%	0.25	338.4	15.05%
	K	2	4411.18	838744.6	98.63%	772.68	641139.2	95.51%	0.27	215	28.64%	0.26	205.2	26.63%
		3	6083.69	348711	97.14%	5170.44	649666	94.11%	0.27	248.2	36.85%	0.41	363.8	34.72%
		4	6244.01	434725.6	93.94%	7200.13	1	91.25%	0.39	384.2	39.64%	0.35	410.2	37.34%
		5	6069.28	142560.8	89.59%	7200.06	1	86.93%	0.39	493.8	39.47%	0.37	484.2	37.04%
	M	2	1010.07	173718.6	99.00%	19.32	23381.8	96.71%	0.21	183	26.69%	0.22	172	23.34%
		3	2168.18	133674.4	97.75%	33.67	26616	95.12%	0.30	296	24.86%	0.26	237.2	21.16%
		4	373.24	98523.8	94.47%	17.35	28016	92.02%	0.35	441.6	23.64%	0.29	346	19.04%
		5	1567.87	43192.4	89.60%	21.58	33498.6	86.93%	0.28	323.8	21.80%	0.24	332.6	16.50%

Table 5. Table 5. Comparison of ( OMPN O ​ T subscript OMPN 𝑂 𝑇 {\rm OMPN}_{OT} ) and ( OMPN B ​ E ​ P subscript OMPN 𝐵 𝐸 𝑃 {\rm OMPN}_{BEP} ) for instances with n = 20 , 30 𝑛 20 30 n=20,30 .

Sc.	$n$	$p$	TimeBEP	StDevBEP	TimeOT	StDevOT	%NonSolvedBEP	%NonSolvedOT	%DevBest
3	20	2	9.73	2.71	11.36	3.58	0%	0%	0%
		5	253.35	18.65	449.32	28.15	0%	2.56%	0.01%
		10	46.97	13.85	77.52	17.47	0%	0%	0%
	30	2	59.64	6.48	148.23	14.70	0%	0%	0%
		5	2931.44	36.11	3099.25	32.91	75%	77.5%	1.63%
		10	2861.03	37.34	3070.86	34.40	77.5%	80%	3.75%
4	20	2	26.45	4.40	30.03	4.95	0%	0%	0%
		5	1865.88	41.15	1874.01	41.33	40%	40%	0.30%
		10	9.51	2.44	22.13	6.17	0%	0%	0%
3	30	2	735.58	36.69	849.51	37.10	15%	15%	0.17%
		5	2742.49	39.04	2836.95	36.91	75%	75%	1.15%
		10	2745.59	39.12	2789.27	38.12	75%	75%	3.28%

Table 6. Table 6. Results of Math-Heuristic Approaches and ( OMPN B ​ E ​ P subscript OMPN 𝐵 𝐸 𝑃 {\rm OMPN}_{BEP} ) in the US dataset.

Sc.	$p$	Pr.	H0	H1	H2	BEP	t0	t1	t2	tBEP	%VarFixed	G1	G2	G0
SC1	2	M	395.3482	394.8909	395.3482	394.891	2.5	31.52	13.28	32.44	89.15%	0.00%	0.12%	0.12%
		C	19.105	18.0278	18.0278	18.0278	0.08	0.28	4.2	8.18	90.41%	0.00%	0.00%	5.64%
		K	272.8855	270.8245	270.7348	270.7348	0.84	12.88	33.24	22.92	89.84%	0.03%	0.00%	0.79%
		D	207.2299	207.2299	207.2299	207.2298	2.99	18.14	12.41	22.12	89.88%	0.00%	0.00%	0.00%
	5	M	222.3974	221.607	221.5985	222.6594	2.68	56.03	160.26	$>3600$	81.54%	-0.47%	-0.48%	-0.12%
		C	18.619	18.2246	16.2491	16.2491	0.17	0.35	9.57	37.59	84.90%	10.84%	0.00%	12.73%
		K	167.6711	163.564	160.1906	160.1138	3.98	24.07	239.04	$>3600$	83.05%	2.11%	0.05%	4.51%
		D	120.9466	120.0572	119.8601	119.7391	2.66	67.81	116.23	$>3600$	83.02%	0.27%	0.10%	1.00%
	10	M	146.2992	144.6164	144.2134	141.8635	8.4	145.96	245.51	$>3600$	85.33%	1.90%	1.63%	3.03%
		C	27.0357	22.9692	19.9944	19.9944	0.19	3.06	15.9	55.57	84.83%	12.95%	0.00%	26.04%
		K	118.7461	118.0626	117.7095	125.0503	6.34	70.23	1760.77	$>3600$	85.80%	-5.92%	-6.24%	-5.31%
		D	86.6926	85.2866	85.7954	85.5946	10.84	94.08	357.67	$>3600$	86.32%	-0.36%	0.23%	1.27%
SC2	2	M	399.6586	395.1866	398.2659	395.1789	2.58	96.14	27.78	462.06	77.53%	0.00%	0.78%	1.12%
		C	23.3894	21.8305	22.1818	21.7935	0.11	0.45	6.72	10.48	85.33%	0.17%	1.75%	6.82%
		K	275.9869	274.8279	274.6485	274.2601	1.4	25.67	32.68	365.95	78.30%	0.21%	0.14%	0.63%
		D	211.7095	209.7885	211.8238	209.7885	2.59	138.32	56.4	262.4	75.94%	0.00%	0.96%	0.91%
	5	M	228.6722	221.514	226.7708	223.3972	637.69	$>3600$	689.41	$>3600$	67.96%	-0.85%	1.49%	2.31%
		C	29.0877	23.0926	22.0907	21.5931	0.55	1.68	15.67	46.04	71.11%	6.49%	2.25%	25.77%
		K	171.2099	168.7613	167.6214	166.7242	21.55	163.7	1025.78	$>3600$	64.07%	1.21%	0.54%	2.62%
		D	129.7708	125.1129	125.4487	127.2981	266.3	600.94	641.76	$>3600$	69.92%	-1.75%	-1.47%	1.91%
	10	M	153.013	153.013	155.285	151.5515	3600.1	$>3600$	3208.72	$>3600$	76.33%	0.96%	2.40%	0.96%
		C	44.5115	35.0981	29.8835	29.842	0.66	16.66	46.47	53.42	74.23%	14.98%	0.14%	32.96%
		K	132.6658	132.6658	131.9971	142.0924	$>3600$	$>3600$	$>3600$	$>3600$	76.33%	-7.11%	-7.65%	-7.11%
		D	100.474	100.474	97.6098	108.2469	$>3600$	$>3600$	3215.14	$>3600$	76.33%	-7.74%	-10.90%	-7.74%
SC3	2	M	404.466	394.9373	398.7388	395.7311	2.94	108.55	49.53	$>3600$	62.96%	-0.20%	0.75%	2.16%
		C	28.0977	24.1051	24.1924	24.1051	0.11	0.66	6.03	10.05	74.69%	0.00%	0.36%	14.21%
		K	280.5643	278.0229	278.9015	278.0228	1.42	32.54	41.81	1455.24	52.65%	0.00%	0.32%	0.91%
		D	216.5167	210.8067	214.7534	210.6248	2.97	74.29	83.65	1488.5	72.21%	0.09%	1.92%	2.72%
	5	M	232.5595	223.4184	235.2989	240.8416	915.52	$>3600$	702.09	$>3600$	60.38%	-7.80%	-2.36%	-3.56%
		C	34.7835	28.4267	27.4501	26.5732	0.23	3.57	13.95	36.93	64.23%	6.52%	3.19%	23.60%
		K	183.8955	171.2615	178.0249	173.4423	128.12	1912.6	1310.85	$>3600$	54.76%	-1.27%	2.57%	5.68%
		D	134.1824	133.0089	132.521	138.3565	1069.6	$>3600$	1206.69	$>3600$	57.96%	-4.02%	-4.40%	-3.11%
	10	M	167.5751	167.5751	171.1057	177.7685	$>3600$	$>3600$	3208.25	$>3600$	70.41%	-6.08%	-3.89%	-6.08%
		C	50.7627	47.0082	39.3866	38.5431	1.1	4.02	25.13	32.2	70.00%	18.01%	2.14%	24.07%
		K	144.9262	144.9262	146.6831	161.303	$>3600$	$>3600$	2336.69	$>3600$	70.41%	-11.30%	-9.97%	-11.30%
		D	110.2827	110.2827	110.5119	116.5976	$>3600$	$>3600$	3214.45	$>3600$	70.41%	-5.73%	-5.51%	-5.73%

Equations

$$\sum _{i=1}^n{\lambda }_i{D}_{(i)},$$

$$\mathrm{DOMP}$$

$$\mathrm{OMPN}$$

$$x_{ij}=\left\{\begin{array}[]{cl}1&\mbox{if client $i$ is allocated to facility $j$ ($i\neq j$) or if facility $j$ is open ($i=j$),}\\ 0&\mbox{otherwise}\end{array}\right.$$

$$\displaystyle\mathcal{X}_{R}=\Big{\{}x\in[0,1]^{n\times n}:\displaystyle\sum_{j=1}^{n}x_{ij}=1,\forall i=1,\ldots,n,\displaystyle\sum_{j=1}^{n}x_{jj}=p,x_{ij}\leq x_{jj},\forall i,j=1,\ldots,n\Big{\}}$$

$$\displaystyle\mathcal{D}=\Big{\{}(d,\bar{a})\in\mathbb{R}^{n\times n}_{+}\times\overline{\mathcal{N}}:d_{ij}\geq\|a_{i}-\bar{a}_{j}\|,i,j=1,\ldots,n,i\neq j\Big{\}},$$

$$\min$$

$$z_i=\sum _{i=1}^n{d}_{ij}{x}_{ij},i,j=1,\dots ,n,$$

$$x\in \mathcal{X},(d,\bar{a})\in \mathcal{D}.$$

$$Z\ge e^t(X-Y),\forall e\in \mathrm{Ext}(P^{\ast }),$$

$$c(\mathcal{B})=\sum _{l=1}^n{\lambda }_l{D}_{(l)}+\sum _{b\in B}f(b)\ge \sum _{l=1}^m{\lambda }_l{D}_{(l)}\ge \sum _{l=1}^m{\lambda }_l{D}_i=D_i\left(\sum _{l=1}^m{\lambda }_l\right)>UB$$

$$\begin{gathered} \Vert A_{i}x-b_i{\Vert }_2\le c_i^{t}x+d_i,\forall i=1,\dots ,M, \\ x\in {\mathbb{R}}^N, \end{gathered}$$

$$\mathcal{D}=\{(d,\bar{a})\in {\mathbb{R}}_{+}^{n\times n}\times {\mathbb{R}}^{n\times d}:d_{ij}\ge \Vert a_i-{\bar{a}}_j{\Vert }_{\nu },r_j\ge \Vert a_j-{\bar{a}}_j{\Vert }_{\tau },i,j=1,\dots ,n\}$$

$$w_{ij}^{k}=\left\{\begin{array}[]{cl}1&\mbox{if customer $i$ is allocated to facility $j$ and its distance, $\|a_{i}-\bar{a}_{j}\|$,}\\ &\mbox{is sorted in the $k$th position.}\\ 0&\mbox{otherwise}.\end{array}\right.$$

$${\mathrm{OMPN}}_{3I}$$

$$\sum _{j,k=1}^n{w}_{ij}^k=1,\forall i=1,\dots ,n,$$

$$\sum _{i,j=1}^n{w}_{ij}^k=1,\forall k=1,\dots ,n,$$

$$\sum _{k=1}^n{w}_{ij}^k\le x_{jj},\forall i,j=1,\dots ,n,$$

$$\sum _{j=1}^n{x}_{jj}=p,$$

$$\sum _{i,j=1}^n{d}_{ij}{w}_{ij}^{k-1}\ge \sum _{i,j=1}^n{d}_{ij}{w}_{ij}^k,\forall k=2,\dots ,n,$$

$$w_{ij}^k\in \{0,1\},\forall i,j,k,=1,\dots ,n,$$

$$x_{jj}\in \{0,1\},\forall j=1,\dots ,n.$$

$$(d,\bar{a})\in \mathcal{D}.$$

$${\theta }_{ij}^k=d_{ij}{w}_{ij}^k,i,j,k=1,\dots ,n.$$

$$\sum _{i,j,k=1}^n{\lambda }_k{\theta }_{ij}^k+\sum _{j=1}^n{f}_j{x}_{jj}.$$

$${\theta }_{ij}^k\ge$$

$$\sum _{i,j=1}^n{\theta }_{ij}^{k-1}\ge$$

$$w_{ij}^k=0,\forall i,j,k=1,\dots ,n,\text{mbox}suchthat\underset{j\ne i}{\min }{\hat{D}}_{ij}>\frac{UB}{\sum _{l=k}^n{\lambda }_l}.$$

$$\sum _{k=1}^n{w}_{jj}^k=x_{jj},\forall j=1,\dots ,n,$$

$$s_{ik}=\left\{\begin{array}[]{cl}1&\mbox{if the distance supported by the $i$th customer is sorted in the $k$th position.}\\ 0&\mbox{otherwise}.\end{array}\right.$$

$${\xi }_k=\sum _{i=1}^n{s}_{ik}\sum _{j=1}^n{d}_{ij}{x}_{ij},k=1,\dots ,n.$$

$${\mathrm{OMPN}}_{2I}$$

$${\xi }_k\ge {\xi }_{k+1},\forall k=1,\dots ,n-1,$$

$$\sum _{k=1}^n{\xi }_k=\sum _{i,j=1}^n{d}_{ij}{x}_{ij},$$

$${\xi }_k\ge d_{ij}{x}_{ij}-{\hat{D}}_{ij}\left(1-s_{ik}\right),\forall i,k=1,\dots ,n,$$

$$\sum _{i=1}^n{s}_{ik}=1,\forall k=1,\dots ,n$$

$$\sum _{k=1}^n{s}_{ik}=1,\forall i=1,\dots ,n,$$

$${\xi }_i\ge 0,\forall i,k=1,\dots ,n,$$

Full text

1. Introduction

Facility location concerns with the optimal placement of one or several new facilities/plants to satisfy the customers’ demands. In discrete facility location problems, the positions of both the customers and the potential new facilities are part of the input as well as the travel costs between them. On the other hand, in continuous facility location problems, although the (geographical) coordinates (in a given $d$-dimensional space) of the customers are provided, the information about the potential location of the facilities is unknown, in the sense that the facilities can be located at any place of the given space. Both the discrete and the continuous versions of facility location problems have been widely studied in the literature (see the monographs [14, 31] and the references therein). Several versions of these facility location problems have been analyzed, by considering different objective functions [43], by fixing either the number of facilities to be located (as in the $p$-median or $p$-center problems) [18] or maximum capacities for the facilities (capacitated facility location)[23, 36], or assuming uncertainty in the demands of the customers (see [2, 9, 10] for a recent review), amongst many others.

In this paper, we propose a unified framework for facility location problems in which the underlying problem is a discrete facility location problem. However, because of locational imprecision or unaccuracy, the new facilities are allowed to be located not only in the exact location of the potential facilities, but in certain regions around each of them, the neighborhoods. In case the initial placements of the potential facilities are exact enough, that is, their neighborhoods are singletons (with a single element which coincides with the initial placement of the potential facilities), the problem becomes the discrete location version of the problem. On the other hand, if the neighborhoods are large enough, the problem turns into the continuous location version of the problem, allowing the facilities to be located in the entire space. Otherwise, different shapes and sizes for the neighborhoods allow one to model how imprecise the locational information provided is. The goal is, apart from the discrete location decision of the problem (placement of facilities among the given set and allocations customers-plants), to find the optimal location of the open facilities in the neighborhoods. The main difference between this problem and its underlying discrete facility location problem, is that in the latest, the travel distances between facilities and customers are assumed to be known, while in the neighborhood version of the problem, as in the continuous case, those distances depend on the place where the facility is located in the neighborhood. Hence, in this problem, the matrix of travel costs is not provided, but a distance measure to compute the travel costs between customers and facilities is given. This problem, as far as we know, has not been fully investigated in Location Analysis, although some attempts have been presented in [8] and [19] where sensitivity analyses were performed by allowing the customers to move around disc-shaped neighborhoods on the plane. Also, this problem can be seen as a constrained version of the classical multifacility location problem, which have been only partially studied in the literature (see [6]). This framework will be called Facility Location with Neighborhoods, a terminology borrowed from the neighborhood versions of the Minimum Spanning Tree problem [4, 13] and Traveling Salesman problem [11, 15].

The importance of analyzing this family of problems comes from its wide range of applications. It is well known that discrete facility location problems are useful in many real-world applications (see [25, 30], amongst many others). However, in many situations, as for instance in the design of telecommunication networks, where a set of servers must be located to supply connection to a set of customers, the exact location of a server may not be exactly provided. In contrast, a region where the decision maker wishes to locate each of the facilities (a corridor, a room, or any other bounded space) can be easily given. In such a case, a robust worst-case decision would not reflect reality, since the decision maker does not known the location of the facility because a lack of certainty but because it allows locational flexibility to the decision. An optimal design may be obtained if the new facilities are allowed to be located in adequately chosen neighborhoods.

In this paper, we provide suitable mathematical programming formulations for the neighborhood versions of a widely studied family of objective functions in facility location problems: ordered median (OM) functions. In these problems, $p$ facilities are to be located by minimizing a flexible objective function that allows one to model different classical location problems. For instance, OM problems allow modeling location problems in which the customers support the median ($p$-median) or the maximum ($p$-center) travel costs, among many other robust alternatives. OM problems were introduced in Location Analysis by Puerto and Fernández [37] and several papers have analyzed this family of objective functions in facility location: discrete problems [21, 24, 35], continuous problems [3, 5], network/tree location problems [20, 42, 45], hub location problems [39], stochastic facility location problems [47], multiobjecive location[17], etc (see [41] for a recent overview on the recent developments on ordered median location problems). In particular, we analyze the neighborhood version of OM location problems for the so-called monotone case. We study the still general case in which the neighborhoods are second-order cone representable regions. These sets allow one to model as particular cases polyhedral neighborhoods or $\ell_{\tau}$-norm balls. The distance measure to represent travel costs between customers and facilities are assumed to be $\ell_{\nu}$-norm based distances. Within this framework we present four different Mixed Integer Second Order Cone Optimization (MISOCO) models.

The current limitations of the on-the-shelf solvers to solve mixed integer nonlinear problems, and the difficulty of solving even the underlying problem (the classical $p$-median problem is NP-hard), makes the resolution of the problem under study a hard challenge. For that reason, we also develop two math-heuristic algorithms based on different location-allocation schemes, which are able to solve larger problems.

Our paper is organized in five sections. In Section 2 we introduce the problem and some general properties are stated. Section 3 is devoted to provide four different mixed integer non linear programming formulations of the problem. At the end of the section, we run some computational experiments in order to compare the four formulations. In Section 4 the two math-heuristic approaches are described, and the results of some computational experiments are reported. Finally, some conclusions are presented in Section 6.

2. DOMP with Neighborhoods

In this section we introduce the Ordered Median Problem with Neighborhoods (OMPN) in which the underlying discrete facility location problem is the Discrete Ordered $p$-Median Problem (DOMP).

For the sake of presentation, we first describe the DOMP problem. The input data for the problem is:

• •

$\mathcal{A}=\{a_{1},\ldots,a_{n}\}\subseteq\mathbb{R}^{d}$: set of coordinates of the customers. We assume, as usual in the location literature, that the coordinates of potential facilities coincides with $\mathcal{A}$.

• •

D=\Big{(}d(a_{i},a_{j})\Big{)}_{i,j=1}^{n}\in\mathbb{R}^{n\times n}: Travel cost matrix between facilities.

• •

$\lambda_{1},\ldots,\lambda_{n}\geq 0$: Ordered median function weights.

The goal of DOMP is to select, from the elements of $\mathcal{A}$, a subset of $p$ facilities, $\mathcal{B}\subset\mathcal{A}$ with $|\mathcal{B}|=p$, that minimizes the ordered median objective function:

[TABLE]

where $D_{i}=\min_{b\in\mathcal{B}}d(a_{i},b)$ (the smallest travel cost to supply customer $i$ from the open facilities), and $D_{(i)}$ represent the $i$-th largest element in the set $\{D_{1},\ldots,D_{n}\}$, i.e. $D_{(i)}\in\{D_{1},\ldots,D_{n}\}$ with $D_{(1)}\geq\cdots\geq D_{(n)}$.

The DOMP can be stated as the following optimization problem:

[TABLE]

We will assume that the $\lambda$-weights verify $\lambda_{1}\geq\cdots\geq\lambda_{n}\geq 0$, dealing with the so-called convex ordered median problem. Most of the main well-known objective functions in Locational Analysis are part of this family, as for instance:

• •

Median ($\lambda=(1,\ldots,1)$): $\sum_{i=1}^{n}D_{i}$.

• •

Center ($\lambda=(1,0,\ldots,0)$: $\max_{i=1,\ldots,n}\;D_{i}$.

• •

$K$–Centrum $\lambda=(1,\stackrel{{\scriptstyle K}}{{\ldots}},1,0,\ldots,0)$: $\sum_{i=1}^{K}D_{(i)}$.

• •

Cent-Dianα ($\lambda=(1,1-\alpha,\cdots,1-\alpha)$): $\alpha\max_{i=1,\ldots,n}D_{i}+(1-\alpha)\sum_{1\leq i\leq n}D_{i}$, for $0\leq\alpha\leq 1$.

Ordered median functions are continuous and symmetric (in the sense that they are invariant under permutations). Furthermore, if $\lambda_{1}\geq\ldots\geq\lambda_{n}\geq 0$, ordered median functions are convex, fact that will be exploited throughout this paper. The interested reader is referred to [41] for a complete description of the properties of ordered median functions.

A few formulations and exact solution approaches for DOMP have been developed since the problem was introduced. In particular Boland et. al [7] formulated the problem as a (non convex) quadratic problem with quadratic constraints. A suitable three index (pure) binary programming reformulation with $O(n^{3}$) variables and $O(n^{2})$ linear constraints was provided by linearizing the bilinear terms. A second formulation, reducing to two the indices of the variables in the formulation, was also presented in the same paper by using a different linearization strategy, and that allows reducing the number of binary variables to $O(n^{2})$. Puerto [36], Marin et. al [28], Marin et. al [27] and Labbé et. al [24] provided alternative formulations for the problem with two and three indices, that need, in a preprocessing phase, sorting the elements in the matrix $D$ (and removing duplicates). All the above mentioned formulations are valid for general ordered median problems. Concerning the convex case, Ogryzack and Tamir [34] presented a different formulation which exploits the monotonicity of the $\lambda$-weights by applying a $k$-sum representation of the ordered median function (see also the recent paper [40] for further details on the powerful of this representation in a wide variety of optimization problems). Finally, in Blanco et. al [5] the authors derived a formulation that also avoid using the binary variables for sorting the involved distances. Also, a few heuristic approaches are available in the literature for the DOMP problem (see [12, 38, 44]).

Observe also that in the DOMP, once the travel costs matrix is provided, the locational coordinates of the customers are not needed, and then, the problem does not depend on the dimension of the space where the customers live.

For the OMPN framework, instead of providing a travel cost matrix between customers, we consider a travel distance measure $d:\mathbb{R}^{d}\times\mathbb{R}^{d}\rightarrow\mathbb{R}_{+}$ induced by a norm $\|\cdot\|$, i.e., $d(a,b)=\|a-b\|$, for $a,b\in\mathcal{A}$.

Also, each potential facility, $a\in\mathcal{A}$, is associated to a convex set, $\mathcal{N}(a)\subset\mathbb{R}^{d}$, with $a\in\mathcal{N}(a)$, its neighbourhood. We denote by $\overline{\mathcal{N}}=\displaystyle\prod_{a\in\mathcal{A}}\mathcal{N}(a)$, the space of neighborhoods. We also consider in this case, set-up costs for opening facilities (which may be neighborhood-dependent), which we denote by $f(a)$ for each $a\in\mathcal{A}$.

The goals of the Ordered Median Problem with Neighborhoods are:

• •

to find the indices of the $p$ facilities to open: $\mathcal{B}=\{b_{1},\ldots,b_{p}\}$,with $b_{j}\in\mathcal{A}$ for $j=1,\ldots,p$,

• •

to locate the facilities into their neighbourhoods: $\bar{b}_{1},\ldots,\bar{b}_{p}$ with $\bar{b}_{j}\in\mathcal{N}(b_{j})$, $j=1,\ldots,p$, and

• •

to allocate customers to their closest open facilities $\bar{b}_{1},\ldots,\bar{b}_{p}$,

by minimizing an ordered median function of the travel distances plus set-up costs.

Observe that the optimization problem to solve for the OMPN is similar to (${\rm DOMP}$):

[TABLE]

but now, $D_{i}=\min_{b\in\mathcal{B}}d(a_{i},\bar{b})$, i.e. the travel distance from a customer to its closest facility depends on the position of the facilities in their neighborhoods. So both the discrete location (open facilities and allocation scheme) and the continuous location decisions (coordinates of the new facilities) are involved in the problem.

We use the classical notation for the variables in $p$-median problem:

[TABLE]

for $i,j=1,\ldots,n$.

Note that, using the above family of variables, the set of open facilities and assignments between customers and $p$ facilitites can be represented by the set $\mathcal{X}=\mathcal{X}_{R}\cap\{0,1\}^{n\times n}$, where

[TABLE]

is the so-called $p$-median polytope.

Observe also that, the above settings easily extend to the case in which the possible connections between demand points and facilities is induced by a graph.

On the other hand, the set of distances, will be represented by the following set:

[TABLE]

where $d_{ij}$ (when one tries to minimize some aggregating function of the travel-costs) represents the distance between the customer located at $a_{i}$ and the facility located at $\bar{a}_{j}$, for all $i,j=1,\ldots,n$.

Note that the set $\mathcal{D}$ can be easily adapted to the case in which each customer uses a different travel distance measure (norm), and the structure of $\mathcal{D}$ remains the same.

With the above notation, the general OMPN can be compactly formulated as:

[TABLE]

where $f_{j}$ denotes the set-up cost of the facility initially located at $a_{j}$, $j=1,\ldots,n$ and $z_{i}$ represents the minimum distance between customer located at $a_{i}$ and the open facilities.

Observe that (2.1)–(2.3) is a mixed integer non linear programming problem (MINLP), whose continuous relaxation is not convex nor concave due to the bilinear constraint (2.2) and probably to the constraints in $\mathcal{D}$. In case the neighborhoods are convex, the set $\mathcal{D}$ is also convex (because of the convexity of the norm). Hence, if the discrete location variables $x$ were known, the problem (also because the convexity of the ordered median function) becomes a continuous convex problem. On the other hand, if the distances were known, the problem becomes a DOMP, so several formulations can be applied to solve the problem. In the OMPN, both $\mathcal{X}$ and $\mathcal{D}$ are part of the final decision. Thus, both the difficulties of handling the DOMP problem and the continuous problem are inherited to the OMPN. In particular, since the $p$-median problem (or the $p$-center problem) is known to be NP-hard [22] which is a particular case of OMPN, the OMPN is also NP-hard.

The simplest OMPN problem, apart from the DOMP case (where the neighborhoods can be seen as singletons), is obtained when the set $\mathcal{D}$ is a polyhedron (and then, defined by a set of linear inequalities). Since the geometry of $\mathcal{D}$ depends on the distance measure induced by $\|\cdot\|$ and the shapes of the neighborhoods, $\mathcal{D}$ will be a polyhedron when these two features can be linearly represented. The norms which are polyhedrally-representable are called block (or polyhedral) norms (see [31]) which are characterized by the fact that their unit balls are polytopes, i.e., $P=\{z\in\mathbb{R}^{d}:\|z\|\leq 1\}$ is a bounded polyhedron. On the other hand, the neighborhoods, because they are assumed to be compact and convex sets, their polyhedral representability is assured if and only if they are also polytopes. In those cases, both the set $\mathcal{D}$ and $\mathcal{X}$ are identified with sets of linear inequalities (and integrality constraints in $\mathcal{X}$). Furthermore, as can be checked in [33] or [41], the ordered median function can be also modeled by using a set of linear inequalities and equations and by adding a set of $O(n^{2})$ binary variables to our model. The above observations are summarized in the following result.

Theorem 2.1**.**

Let $\lambda_{1}\geq\cdots\geq\lambda_{n}\geq 0$, $\|\cdot\|$ a block norm and $\mathcal{N}(a)$ a polyhedron, for each $a\in\mathcal{A}$. Then OMPN can be formulated as a mixed-integer linear programming problem.

Proof.

The proof follows noting that constraints in the form $Z\geq\|X-Y\|$, as those that appear in the description of $\mathcal{D}$, can reformulated as:

[TABLE]

where ${\rm Ext}(P^{*})$ the set of extreme points of $P^{*}=\{v\in\mathbb{R}^{d}:v^{t}b_{g}\leq 1,g=1,\ldots,|Ext(P)|\}$, the unit ball of the dual norm of $P$ (see [31, 46] ). ∎

The following example illustrates the new framework under study.

Example 2.2*.*

Let us consider a set of customers/potential facilities with coordinates in the plane $\mathcal{A}=\{(0,5)$, $(1,1)$, $(1,6)$, $(1,4)$, $(5,3)$, $(10,4)$, $(6.5,0)$, $(8,6)\}$, and travel distances measured with the Euclidean norm. The solutions for the $2$-median, the $2$-center and the $2$-$4$-center ($K$-center with $K=4$) are drawn in Figure 1 (there, stars represent open facilities, customers are identified with dots and lines are the allocation patterns).

Note that, as expected, the solutions for the three DOMP problems highly depend on the $\lambda$-weights, being the optimal set of open facilities different for the three cases.

Let us now consider, for each demand point, a neighbourhood defined as the Euclidean disk with radii $r\in\{1,0.6,1,0.6,2.4,2.4,0.8,1.6\}$ (see Figure 2).

The new facilities, now, are not restricted to be exactly located in the given coordinates but in a disk around them. We consider the radius of its neighborhood (disk) as a mesure of the set-up cost for each facility, that is $f(a_{i})=r_{i}$. The solutions of the neighbourhood version of the $2$-median, $2$-center and $2$-$4$-center problems are shown in Figure 3.

In what follows, we derive some structural properties of DOMP that are inherited to the OMPN.

For each $i,j=1,\ldots,n$, we denote by $\widehat{D}_{ij}=\max\{\|a_{i}-\bar{a}_{j}\|:\bar{a}_{j}\in\mathcal{N}_{j}\},$ and $\widehat{d}_{ij}=\max\{\|a_{i}-\bar{a}_{j}\|:\bar{a}_{j}\in\mathcal{N}_{j}\}$, upper and lower bounds for the distances between the $i$th customer and the $j$th potential facility, respectively.

Properties 2.3**.**

The following properties are satisfied:

1. (1)

There exists an optimal solution of (${\rm OMPN}$) in which the $p$ smaller travel distances equal [math]. 2. (2)

Let $\mathcal{B}\subseteq\mathcal{A}$ a set of $p$ facilities such that its ordered cost $C(\mathcal{B})\leq UB$ and such that $\min_{j\neq i}\widehat{d}_{ij}>\dfrac{UB}{\displaystyle\sum_{i=1}^{m}\lambda_{i}}$ for some $m=2,\ldots,n$, then the $i$-th client is sorted at most in position $m$ in the whole sorted set of optimal distances.

Proof.

1. The result follows from the observation that if the facility $a\in\mathcal{A}$ is an open facility, the travel costs between $a$ and $a$ are zero.

2. Assume that the $i$th customer is sorted in position $r\geq m$ in the sorting sequence of distances, i.e., $D_{(1)}\geq\ldots\geq D_{(m)}\geq D_{(r)}=D_{i}$. Then, we have that:

[TABLE]

which contradicts the hyphotesis. ∎

3. MINLP Formulations for the OMPN

In this section, we describe different mathematical programming formulations for solving general OMPN. In particular, we extend the formulations presented in [5], [7] and [34] to our problem. As mentioned above, the main difference between the DOMP and the OMPN problem is that in the OMPN the distances are not part of the input, but part of the decision. Hence, the formulations for the DOMP based on preprocessing the travel distances matrix (as those proposed in [24],[28], [27] or [36]) cannot be applied to our framework.

Observe that, in OMPN, an adequate representation of $\mathcal{D}$ is crucial for the development of efficient solution approaches for the problem. We assume that the neighborhoods belong to a family of convex sets that allows us to represent most of the convex shapes which are useful in practice, and that can be efficiently handled by commercial optimization solvers: second order cone (SOC)-representable sets [26]. SOC-representable sets are convex sets defined by second-order cone constraints in the form:

[TABLE]

where $A_{1},\ldots,A_{i}\in\mathbb{R}^{M_{i}\times N}$, $b_{i}\in\mathbb{R}^{M_{i}}$, $c_{i}\in\mathbb{R}^{N}$, $d_{i}\in\mathbb{R}$, for $i=1,\ldots,M$, and $\|\cdot\|_{2}$ is the Euclidean norm. Most of the state-of-the-art solvers are capable to efficiently solve optimization problems involving SOC constraints by means of quadratic constraints with positive definite matrices, second order cone constraints (in the form $x^{t}x\leq y^{2}$, for $y\geq 0$) or rotated second order cone constraints ($x^{t}x\leq yz$ with $y,z\geq 0$). SOC constraints allow one to represent, not only Euclidean balls, but any $\ell_{\tau}$-norm ball (see [5] for further details on the explicit representation of $\ell_{\tau}$-norm based distance constraints as a set of SOC constraints for any $\tau\in\mathbb{Q}$ with $\tau\geq 1$). Clearly, any polyhedron is SOC-representable (setting $A$ and $b$ equal to zero) so any intersection of $\ell_{\tau}$-norm balls and polyhedra is suitable to be represented as a set of second order cone constraints. Hence, both our neighborhoods and the distances involved in our problem will be defined as SOC-constraints, being then $\mathcal{D}$ a SOC-representable set.

For the sake of simplicity, and without loss of generality, we assume that the neighborhood of each $a\in\mathcal{A}$ is a $\ell_{\tau}$-norm ball, i.e. $\mathcal{N}(a)=\{z\in\mathbb{R}^{d}:\|z-a\|_{\tau}\leq r_{a}\}$, for some $r_{a}\in\mathbb{R}_{+}$ and $\tau\in\mathbb{Q}$ with $\tau\geq 1$.

Also, we consider that the travel distances are induced by a $\ell_{\nu}$-norm with $\nu\in\mathbb{Q}$ and $\nu\geq 1$. With these settings, we explicitly describe $\mathcal{D}$ as follows:

[TABLE]

where $r_{j}$ denotes the radius of the neighborhood $\mathcal{N}(a_{j})$, i.e., $r_{j}=r_{a_{j}}$.

The following result, whose proof is straightforward from [5, Theorem 2], allows us to efficiently represent the set $\mathcal{D}$ when the involved norms are $\ell_{\tau}$-based norms.

Properties 3.1**.**

Let $\tau=\frac{r_{\tau}}{s_{\tau}}\geq 1$ and $\nu=\frac{r_{\nu}}{s_{\nu}}\geq 1$ with $r_{\tau},s_{\tau},r_{\nu},s_{\nu}\in\mathbb{Z}_{+}$ and $\gcd(r_{\tau},s_{\tau})=\gcd(r_{\nu},s_{\nu})=1$. Then, $\mathcal{D}$ is representable as a set of $(n^{2}+n)(2d+1)$ linear inequalities and $nd(n\log\;r_{\nu}+\log\;r_{\tau})$ second order cone contraints.

3.1. The three index formulation

The first formulation is based on the one proposed in [7], which uses, apart from the $x_{jj}$-variables described above, the following set of sorting/allocation binary variables for the DOMP:

[TABLE]

This formulation reads as follows:

[TABLE]

The objective function assigns to each sorted distance its adequate weight $\lambda$. (3.1) (resp. (3.2)) ensures that each demand point (resp. each position) is assigned to a unique facility and a unique position (resp. demand point). (3.3) assures that allocation is not allowed unless the plant is open, and (3.4) restrict the problem to open exactly $p$ facilities. Constraints (3.5) allows us a correct definition of the $w$-variables in which the sorting of the distances is imposed (the $(k-1)$th is at least as larger as the $k$-th distance).

Although the above formulation is valid for OMPN, both the objective function and the set of constraints (3.5) are quadratic and non-convex. We introduce a new set of variables to account for the non linear terms in the above formulation:

[TABLE]

Using the $\theta$-variables, the objective function can be reformulated as:

[TABLE]

The correct definition of the new variables and satisfaction of constraints (3.5) is assured by the following sets of linear constraints:

[TABLE]

where the first set of constraints comes from the McCormick linear reformulation [29] of the bilinear terms defining the $\theta$-variables, and the second is the reformulation of (3.5) with the new variables.

The formulation above, has $O(n^{3})$ variables and $O(n^{2})$ constraints. Properties 2.3 allow us to strengthen the formulation (${\rm OMPN}_{3I}$). In particular, if $UB$ is a known upper bound for the optimal value of OMPN, then:

[TABLE]

Also, because of the relationship between the $w$ and $x$ variables we get that

[TABLE]

are valid equations for (${\rm OMPN}_{3I}$).

3.2. The $2$-index formulation

The second formulation, also based on the one presented in [7], considers an alternative representation of the sorting variables. It uses two different sets of variables. The first one allows us to sort the distances of supplying each of the customers:

[TABLE]

while the second represents the sorting (non decreasing) sequence of distances:

[TABLE]

This representation allows us to simplify the formulation to the following with $O(n^{2})$ variables and $O(n^{2})$ constraints.

[TABLE]

[TABLE]

The correct definition of the $\xi$-variables is assured by constraints (3.6)–(3.8), while contraints (3.9) and (3.10) allows us the adequate modeling of the $s$-variables.

As in (${\rm OMPN}_{3I}$), to avoid nonconvex terms in the formulation, the bilinear terms $d_{ij}x_{ij}$ can be linearized by introducing a new variable $\theta_{ij}=d_{ij}x_{ij}$ and replacing (3.7) and (3.8) by:

[TABLE]

3.3. The $K$-sum formulation

Ogryczak and Tamir presented in [34] some linear programming formulations for the problem of minimizing the sum of the $K$ largest linear functions (which is a particular case of ordered median function). In the same paper, the approach is extended to the minimization of convex ordered median functions by means of a telescopic sum of $K$-sum functions. In the next formulation, we apply this idea to formulate the OMPN. For the sake of readability, we first formulate the $K$-center problem.

Let $\lambda=(1,\stackrel{{\scriptstyle K)}}{{\ldots}},1,0,\stackrel{{\scriptstyle n-K)}}{{\ldots}},0)$. The ordered median function associated to this particular choice of $\lambda$ is known as the $K$-center function. With our notation, provided the set of distances $D_{1},\ldots,D_{n}$, the $K$-center problem consists of minimizing $\Theta_{K}(D)=\sum_{i=1}^{K}D_{(i)}$. Such an objective function, is proved in [34] to be equivalent to the following expression

[TABLE]

where $z_{+}=\max\{0,z\}$ for $z\in\mathbb{R}$, and the optimal value $t^{*}$ into the above expression coincides with $D_{(K)}$ (the $K$-th largest distance). Hence, to minimize $\Theta_{K}(D)$ one may proceed by solving:

[TABLE]

where the variable $z_{i}$ is identified with $(D_{i}-t)_{+}$ in the above formulation, for $i=1,\ldots,n$. Thus, incorporating the whole information to define the distances in our location problem, the $K$-center location problem with neighborhoods can be formulated as:

[TABLE]

The above formulation can be extended to general convex ordered median functions. Observe that if $\lambda_{1}\geq\cdots\geq\lambda_{n}\geq\lambda_{n+1}:=0$ one may represent the ordered median function of the distances $D_{1},\ldots,D_{n}$ by using a telescopic sum:

[TABLE]

where $\Delta_{k}=\lambda_{k}-\lambda_{k+1}\geq 0$ for $k=1,\ldots,n-1$.

Thus, the convex ordered objective functions can be equivalently rewritten as a weighted sum of $K$-sums, being then suitable to be represented as in the $K$-center problem. With such an observation and introducing new $t$-variables (one for each of the $K$-sums involved) and $z$-variables, in this case with two indices to account not only for the customer ($i$) but also for the $K$-sum representation, one obtain the following valid formulation for the OMPN:

[TABLE]

Observe that this formulation has also $O(n^{2})$ variables and $O(n^{2})$ constraints, but, as will be shown in the computational results, it has a better performance than (${\rm OMPN}_{2I}$) since it uses, intrinsically, the especial structure of the convex ordered median objective.

3.4. The BEP formulation

Finally, we present a formulation, based on the one provided in [5] and that, as the one in the previous subsection, is only valid for the convex case. The idea under the formulation comes from the observation that, because $\lambda_{1}\geq\cdots\geq\lambda_{n}\geq 0$, the evaluation of the ordered median function on a set of distances, is reached when choosing, among all possible permutations of the indices, $\mathcal{P}_{n}$, the one that maximizes the weighted sum, that is:

[TABLE]

Then, if the permutations of $\{1,\ldots,n\}$ are represented by using the set of binary variables

[TABLE]

verifying that $\sum_{i=1}^{n}p_{ik}=1$ (for all $k=1,\ldots,n$) and $\sum_{k=1}^{n}p_{ik}=1$ (for all $i=1,\ldots,n$).

Then, using these variables, the ordered median sum of a given set of values $D_{1},\ldots,D_{n}$ is equivalent to:

[TABLE]

The optimization problem above is an assignment problem, then, the total unimodularity of the constraints matrix assures that its optimal value coincides with the one of its dual problem which reads:

[TABLE]

Merging the above representation of the ordered median function into the location problem, the OMPN is reformulated as:

[TABLE]

3.5. Comparison of Formulations

In this section we provide some theoretical results that allow us to compare the tightening of the convex relaxations of each of the provided formulations. Let us denote by $z^{R}_{3I}$, $z^{R}_{2I}$, , $z^{R}_{OT}$ and $z^{R}_{BEP}$ optimal values of the continuous relaxations of formulations (${\rm OMPN}_{3I}$), (${\rm OMPN}_{2I}$), (${\rm OMPN}_{OT}$) and (${\rm OMPN}_{BEP}$), respectively.

Proposition 3.2**.**

The following relations are verified:

[TABLE]

Proof.

Denote by $F_{3I}$, $F_{2I}$, $F_{OT}$ and $F_{BEP}$ the feasible regions of (${\rm OMPN}_{3I}$), (${\rm OMPN}_{2I}$), (${\rm OMPN}_{OT}$) and (${\rm OMPN}_{BEP}$) obtained when relaxing the integrality conditions of the models.

• •

Let us consider the mapping $\pi:\mathbb{R}^{n}_{+}\times\mathbb{R}^{n\times n}_{+}\times[0,1]^{n\times n}\times\mathcal{X}_{R}\times\mathcal{D}\rightarrow\mathbb{R}^{n^{3}}_{+}\times[0,1]^{3}\times[0,1]^{n}\times\mathcal{D}$ defined as:

[TABLE]

First, let us check that $\pi(F_{2I})\subseteq F_{3I}$, which would prove the first inequality. Let $(\theta,\xi,s,x,(d,\bar{a}))\in F_{2I}$, and define $(\bar{\theta},\bar{x},(d,\bar{a}))=\pi(\theta,\xi,s,x,(d,\bar{a}))$, i.e.:

[TABLE]

By construction, the constraints (3.1)-(3.4) are verified:

• –

$\displaystyle\sum_{j,k=1}^{n}\bar{w}_{ij}^{k}=\displaystyle\sum_{j,k=1}^{n}s_{ik}x_{ij}=\displaystyle\sum_{j=1}^{n}x_{ij}\displaystyle\sum_{k=1}^{n}s_{ik}\stackrel{{\scriptstyle\eqref{domp2:5}}}{{=}}\displaystyle\sum_{j=1}^{n}x_{ij}\stackrel{{\scriptstyle x\in\mathcal{X}_{R}}}{{=}}1$.

• –

$\displaystyle\sum_{i,j=1}^{n}\bar{w}_{ij}^{k}=\displaystyle\sum_{i,j=1}^{n}s_{ik}x_{ij}=\displaystyle\sum_{i=1}^{n}s_{ik}\displaystyle\sum_{j=1}^{n}x_{ij}\stackrel{{\scriptstyle x\in\mathcal{X}_{R}}}{{=}}\displaystyle\sum_{i=1}^{n}s_{ik}\stackrel{{\scriptstyle\eqref{domp2:4}}}{{=}}1$.

• –

$\displaystyle\sum_{k=1}^{n}\bar{w}_{ij}^{k}=\displaystyle\sum_{k=1}^{n}s_{ik}x_{ij}=x_{ij}\displaystyle\sum_{k=1}^{n}s_{ik}\stackrel{{\scriptstyle\eqref{domp2:5}}}{{=}}x_{ij}\stackrel{{\scriptstyle x\in\mathcal{X}_{R}}}{{=}}x_{jj}$.

• –

$\displaystyle\sum_{j=1}^{n}\bar{x}_{jj}=\displaystyle\sum_{j=1}^{n}x_{jj}=\stackrel{{\scriptstyle x\in\mathcal{X}_{R}}}{{=}}p$.

• –

$\bar{\theta}_{ij}^{k}=\xi_{k}s_{ik}x_{ij}\stackrel{{\scriptstyle\eqref{domp2:7},\eqref{domp2:8}}}{{\geq}}(d_{ij}-\widehat{D}_{ij}(2-s_{ik}-x_{ij}))\,s_{ik}x_{ij}=d_{ij}-\widehat{D}_{ij}(1-s_{ik}x_{ij})+(\widehat{D}_{ij}-d_{ij})(1-s_{ik}x_{ij})+\widehat{D}_{ij}s_{ik}x_{ij}(s_{ik}+x_{ij})\stackrel{{\scriptstyle\bar{w}_{ij}^{k}=s_{ik}x_{ij},s_{ik},x_{ij}\leq 1,d_{ij}\leq\widehat{D}_{ij}}}{{\geq}}d_{ij}-\widehat{D}_{ij}(1-\bar{w}_{ij}^{k})$.

• –

$\displaystyle\sum_{i,j=1}^{n}\bar{\theta}_{ij}^{k}=\displaystyle\sum_{i,j=1}^{n}\xi_{k}s_{ik}x_{ij}=\xi_{k}\displaystyle\sum_{i=1}^{n}s_{ik}\displaystyle\sum_{j=1}^{n}x_{ij}\stackrel{{\scriptstyle x\in\mathcal{X}_{R}}}{{=}}\xi_{k}\displaystyle\sum_{i=1}^{n}s_{ik}\stackrel{{\scriptstyle\eqref{domp2:4}}}{{=}}\xi_{k}\stackrel{{\scriptstyle\eqref{domp2:1}}}{{\geq}}\xi_{k+1}=\displaystyle\sum_{i,j=1}^{n}\bar{\theta}_{ij}^{k+1}$.

Then, $\pi(\theta,\xi,s,x,(d,\bar{a}))\in F_{3I}$, so $\pi(F_{2I})\subset F_{3I}$, i.e., any solution of the convex relaxation of (${\rm OMPN}_{2I}$) induces a solution of the convex relaxation of (${\rm OMPN}_{3I}$). Furthermore, the objective values for $(\theta,\xi,s,x,(d,\bar{a}))$ in (${\rm OMPN}_{2I}$) and $\pi(\theta,\xi,s,x,(d,\bar{a}))$ in (${\rm OMPN}_{3I}$) coincides:

[TABLE]

Thus, $z_{2I}^{R}\geq z_{3I}^{R}$.

• •

Let $(u,v,D,x,(d,\bar{a}))\in\mathbb{R}^{n}\times\mathbb{R}^{n}\times\mathbb{R}^{n}_{+}\times\mathcal{X}_{R}\times\mathcal{D}$ be the optimal solution of the continuous relaxation of (${\rm OMPN}_{BEP}$). Let $p_{ik}$ the optimal dual variables associated to constraint (3.21). By optimality conditions they must verify:

[TABLE]

Let us construct the following vector in $\mathbb{R}^{n}_{+}\times\mathbb{R}^{n\times n}_{+}\times[0,1]^{n\times n}\times\mathcal{X}_{R}\times\mathcal{D}$:

[TABLE]

By construction, $\bar{s}$ clearly verifies (3.9) and (3.10). Also, note from the construction of the BEP formulation that for given $D_{1},\ldots,D_{n}$, the problem

[TABLE]

is equivalent to

[TABLE]

which is an assignment problem related to the best sorting on the variables based on their costs given by $D_{1},\ldots,D_{n}$. Because the monotonicity and nonnegativity of the $\lambda$-weights, this is equivalent to compute the ordered median sum $\displaystyle\sum_{i=1}^{n}\lambda_{i}D_{(i)}=\displaystyle\sum_{i,k=1}^{n}\lambda_{k}D_{i}p_{ik}$ (where $p$ are the corresponding solution to the problem above indicating if $p_{ik}=1$ that element $i$ is sorted in the $k$th position). Hence, $\bar{\xi}_{k}=\displaystyle\sum_{i=1}^{n}p_{ik}D_{i}\geq\displaystyle\sum_{i=1}^{n}p_{ik+1}D_{i}=\bar{\xi}_{k+1}$ (constraint (3.6)). The proof of the verification of remainder constraints are straightforward. Also, the reader can easily check that the objective values of both solutions coincide. Thus, $z^{R}_{BEP}\geq z^{R}_{3I}$

• •

Let $(u,v,D,x,(d,\bar{a}))\in\mathbb{R}^{n}\times\mathbb{R}^{n}\times\mathbb{R}^{n}_{+}\times\mathcal{X}_{R}\times\mathcal{D}$ be the optimal solution of the continuous relaxation of (${\rm OMPN}_{BEP}$). Let us construct a feasible solution for the continuous relaxation of (${\rm OMPN}_{OT}$). Let $(\bar{t},\bar{z})\in\mathbb{R}^{n}\times\mathbb{R}^{n\times n}_{+}$ the solution to the problem

[TABLE]

By the construction in Subsection 3.3, the vector $(\bar{t},\bar{z},D,x,(d,\bar{a}))$ is a feasible solution to the continuous relaxation of (${\rm OMPN}_{OT}$) with same objective value than $(u,v,D,x,(d,\bar{a}))$ in the continuous relaxation of (${\rm OMPN}_{BEP}$), being then $z_{OT}^{R}\leq z_{BEP}^{R}$. Observe that the oposite direction can be derived with a similar reasoning. Thus, $z_{OT}^{R}=z^{R}_{BEP}$.

∎

In the above result, the relationship between $z^{R}_{2I}$ and $z^{R}_{BEP}$ (or $z^{R}_{OT}$) is not stated. One may think that the continuous relaxation of (${\rm OMPN}_{BEP}$) is tightener than (${\rm OMPN}_{2I}$), because the first exploits the monotonicity of the $\lambda$-weights. However, that is not always true as illustrated in the following example.

Example 3.3*.*

Let us consider five points in $\mathbb{R}^{2}$, $\mathcal{A}=\{(2,92),$ $(33,70),$ $(65,50),$ $(73,69),$ $(40,63)\}$ and neighborhoods defined as Euclidean disks with radii $\{2,1,0.05,5,1\}$. If the distance measure is the Euclidean norm and we choose as ordered median function the one with $\lambda=(1,1,1,1,1)$ and $p=2$, we get that:

[TABLE]

where $z_{OPT}$ is the optimal value of the OMPN.

3.6. Computational Comparison of Relaxations

We have run a series of experiments to study the computational performance of the formulations (${\rm OMPN}_{3I}$), (${\rm OMPN}_{2I}$), (${\rm OMPN}_{OT}$) and (${\rm OMPN}_{BEP}$) and also to know the computational limits of the OMPN. We have randomly generated instances of $n$ demand points in $[0,100]^{2}$ and $[0,100]^{3}$ with $n$ ranging in $\{5,6,7,8,9,10,20,30\}$. Five random instances were generated for each number of points. We solved OMPN problems $p$ (number of new facilities to be located) ranging in $\{2,3,5\}$ (provided that $p<n$). Euclidean distances were considered to measure the distances between points. We considered neighborhoods defined as discs (for the planar instances) or 3-dimensional Euclidean balls (for the $3$-dimensional instances) with randomly generated radii. The sizes of the neighborhoods were generated under four scenarios:

**Scenario 1.: **

Radii generated in $[0,5]$.

**Scenario 2.: **

Radii generated in $[5,10]$.

**Scenario 3.: **

Radii generated in $[10,15]$.

**Scenario 4.: **

Radii generated in $[15,20]$.

In Figure 4 we show, for one of our $8$-points instances, the neighborhoods for each the four scenarios. Note that Scenario 1 slightly differs from the DOMP while Scenario 4 is closer to the continuous problem. As will be observed from our experiments, the computational difficulty of solving problems with larger radii is higher than the one of those with small radii.

The set-up cost of each facility was assumed to be the radius of its neighborhood. It can be interpreted as the cost of covering the neighborhood (larger as $r$ increases).

The four formulations were coded in C, and solved using Gurobi 7.01 in a Mac OSX El Capitan with an Intel Core i7 processor at 3.3 GHz and 16GB of RAM. A time limit of 1 hour was set in all the experiments.

Also, four different convex ordered median problems were solved for each of the instances:

$p$**-Median (M):: **

$\lambda=(1,\ldots,1)$.

$p$**-Center (C):: **

$\lambda=(1,0,\ldots,0)$.

$p$**-KCenter (K):: **

$\lambda=(\overbrace{1,\ldots,1}^{\lfloor\frac{n}{2}\rfloor},0,\ldots,0)$.

$p$**-Cent-Dian0.5 (D):: **

$\lambda=(1,0.5,\ldots,0.5)$.

For the small instances ($n\leq 10$) we report here only the results obtained for $n=10$ (the interested reader may check the complete result of the experiments for $n=5,\ldots,10$ in bit.ly/resultsDOMPN)

First, we run the continuous relaxation of the four formulations to analyze their integrality gaps. The average results are shown in Table 1. We report the integrality gaps ($IG=\dfrac{z^{*}}{z^{R}}$) for each of the formulations ((${\rm OMPN}_{3I}$), (${\rm OMPN}_{2I}$) and (${\rm OMPN}_{BEP}$) (obviating (${\rm OMPN}_{OT}$) whose integrality gap coincides with the one of (${\rm OMPN}_{BEP}$) by Theorem 3.2). The table summarizes the average integrality gaps for the each of the scenarios. As remarked above, there is no order relation between $z^{R}_{2I}$ and $z^{R}_{BEP}$. We have bolfaced those values in which the average integrality gaps of (${\rm OMPN}_{2I}$) is smaller than the one for (${\rm OMPN}_{BEP}$). Note that it only happens for a few combinations of $n$, $p$, problem types and scenarios. In particular, it mostly occurs for small values of $p$ and only for Scenario 1. In Table 2 we show the percentage of instances (of all of those with fixed $n$ and $p$) for which $z_{2I}^{R}\geq z_{BEP}^{R}$. In the total instances, this percentage is $9.97\%$, while among the instances of Scenario 1 is $26.15\%$.

Tables 3 and 4 show the results for the planar and 3-dimensional problems, respectively. For each of the formulations, we provide the average CPU times (time), the number of nodes explored in the branch-and-bound tree (#nodes) and the deviation with respect to the solution obtained at the end of the root node of the search tree (%gapR).

As can be observed from the results, formulations (${\rm OMPN}_{OT}$) and (${\rm OMPN}_{BEP}$) are much less time consuming than (${\rm OMPN}_{3I}$) and (${\rm OMPN}_{2I}$) in all the cases. Also, the solutions obtained after exploring the root node of (${\rm OMPN}_{OT}$) and (${\rm OMPN}_{BEP}$) are tighten than the rest. Consequently, the number of explored nodes to find the optimal solution or to certify optimality is higher in the two first formulations. Observe that the results are as expected since the first two formulations do not exploit the convexity of monotone ordered median problems. Observe that the sorting constraints in the first two formulations involve binary variables while in the two last formulations no need of new binary variables are needed for this task.

Since (${\rm OMPN}_{OT}$) and (${\rm OMPN}_{BEP}$) seems to have a similar computational behavior for the small-size instances, we have performed a series of experiments for medium-size instances to compare these two formulations. The results are shown in Table 5, where now $n$ ranges in $\{20,30\}$ and $p$ in $\{2,5,10\}$. As can be observed, the performance (in terms of CPU time) of both formulation is similar, but (${\rm OMPN}_{BEP}$) seems to need, in average, less CPU time to solve the problems in most of the problems, and the standard deviations (StDev) of the consuming times for (${\rm OMPN}_{BEP}$) are smaller than those for (${\rm OMPN}_{OT}$). Furthermore, we were able to solve all the instances before the 1 hour time limit, but $2.56\%$ of them by using (${\rm OMPN}_{BEP}$), while (${\rm OMPN}_{OT}$) was not able to solve $11.34\%$ of the the problems. Moreover, in all the instances, (${\rm OMPN}_{BEP}$) obtained better upper bounds for the optimal value of the problems in all the instances (the deviation of the best upper bounds obtained with the OT formulation with respect to the best solution obtained with the BEP formulation is shown in column %DevBest).

4. Math-heuristics for the OMPN

In this section we describe two mathematical programming location-allocation based heuristic approaches for solving the OMPN for larger instances. Some heuristics have been proposed for solving ordered $p$-median problems (see [12]). However, most of them are based on the use of “fast” procedures to compute the overall cost of opening/closing certain sets of facilities. Note that when the travel cost matrix is provided and a set of open facilities is obtained, one can easily evaluate, for each customer, its cheapest facility (or its second cheapest facility), and once all of them are computed, evaluate an ordered median function can be also efficiently performed. In the OMPN, even if the open facilities are known, the allocation costs depend on the final location of the facilities (which depends also on the customers allocated to each facility). Hence, the developed heuristics for the DOMP are no longer valid for the OMPN problem. We propose two alternative local search math-heuristic which allows us to solve larger instances of the problem at smaller computational costs than the exact approaches, at the price of not warrantying the optimality of the solution.

The two heuristics procedures that we propose are part of the well-known familiy of location-allocation procedures. These type of schemes based on performing changes over a set of $p$ facilities candidates for being opened, trying to improve the incumbent objective value. For the sake of that we need to compute the overall cost of opening a given set of facilities $J\subseteq\{1,\ldots,n\}$ with $|J|=p$. Observe that if a set of open facilities is known, we have to compute the minimum (ordered weighted) cost of allocating the customers to those facilities. As mentioned above the computation of such a cost will involve the computation of the allocation customer-open facility and also the position of the open facilities inside their neighborhoods. Although different formulations can be used for such a task, we will use the one based on formulation (${\rm OMPN}_{BEP}$).

Proposition 4.1**.**

Let $J\subset N=\{1,\ldots,n\}$ with $|J|=p$. Then, the cost of using $J$ as open facilities of the OMPN problem can be computed by solving the following mixed integer non linear programming problem:

[TABLE]

Proof.

The proof easily follows noting that (${\rm ALLOC}(J)$) is nothing but the simplification of (${\rm OMPN}_{BEP}$) when the values of $x_{jj}$ are known and fixed to $x_{jj}=1$ if $j\in J$ and $x_{jj}=0$, otherwise. ∎

For each $J$, (${\rm ALLOC}(J)$) can be reformulated as a mixed integer second order cone constraint problem with $(n-p)p$ binary variables (instead of the $n^{2}$ in (${\rm OMPN}_{BEP}$). Furthermore, a variable fixing strategy can be applied to (${\rm ALLOC}(J)$) in order to reduce the number of binary variables of the problem.

Properties 4.2**.**

Let $J\subseteq N$ with $|J|=p$, $i\in N\backslash J$, $j\in J$ and $x^{*}\in\mathcal{X}$ optimal allocation solutions of (${\rm ALLOC}(J)$).

1. (1)

If $\exists k\in J\backslash\{j\}$ such that $\widehat{D}_{ik}<\widehat{d}_{ij}$, then $x^{*}_{ij}=0$. 2. (2)

If $\min_{k\neq j}\widehat{d}_{ik}>\widehat{D}_{ij}$, then $x^{*}_{ij}=1$. 3. (3)

If $\{j^{\prime}\in J:\widehat{D}_{ik}\geq\widehat{d}_{ij^{\prime}},\forall k\neq j^{\prime}\}=\{j\}$, then $x^{*}_{ij}=1$.

Proof.

1. (1)

Let us assume that $x^{*}_{ij}=1$. Then, $d_{ij}^{*}=\|a_{i}-\bar{a}^{*}_{j}\|\geq\widehat{d}_{ij}$. By hypothesis there exists $k\in J$ ($k\neq j$) with $\widehat{D}_{ik}<\widehat{d}_{ij}$. Hence, $d_{ij}^{*}>\widehat{D}_{ik}\geq d_{ik}=\|a_{i}-\bar{a}_{k}^{*}\|$, so $d_{ij}^{*}\neq\min_{j^{\prime}\in J}\|a_{i}-\bar{a}_{k}\|$ contradicting the optimality of the solution. 2. (2)

If $x^{*}_{ij}=0$, then, there exists $k\in J$, $k\neq j$ such that $d_{ij}^{*}\geq d_{ik}^{*}>\widehat{D}_{ij}\geq d_{ij}^{*}$. Thus, $x^{*}_{ij}=1$. 3. (3)

If applying 1, all the facilities except $j$ must verify $x_{ij^{\prime}}^{*}=0$, then the unique choice for allocating $i$ is $j$.

∎

As we will show in our computational experiments, the above strategies for fixing variables allows us to fix an average of $80\%$ of the binary variables in the test problems.

Using the above-described formulation, we implemented two different heuristic algorithms. Both algorithms will move through different feasible solutions in order to improve an initial feasible solution. This initial solution is constructed by either solving the standard DOMP problem with a set of weights based on the distances between centers of the neighborhoods (a convex combination of $\widehat{D}_{ij}$ and $\widehat{d}_{ij}$), or solving the OMPN problem for simpler neighborhoods (as polyhedral neighborhoods) and polyhedral distances (which may require less computational effort than general $\ell_{\tau}$-norm based metrics or neighborhoods). Hence, we consider that an initial solution $x^{0}\in\mathcal{X}$ is known.

4.1. Math-heuristic Algorithm 1

Given a feasible solution $\bar{x}\in\mathcal{X}$, the first algorithm searches, for each facility $j_{0}$ in $J$, the best replacement by a facility in $N\backslash J$. Two different options are possible here. First, to construct the new set of open facilities $J^{\prime}=J\cup\{i\}\backslash\{j_{0}\}$ for each $i\in N\backslash J$, solve (${\rm ALLOC}(J)$) for such a $J^{\prime}$ and keep the best possible change for $j_{0}$. The second option is to solve a single mixed integer non linear programming problem which decides (through the binary variable $\xi_{i}$, whether the non-opened facility $i$ is interchanged by $j_{0}$ to obtain the best improvement:

[TABLE]

Note that constraints (4.1) are the linearization of the bilinear terms as in the previus formulations, but obviating the facility that wants to be replaced ($j_{0}$). For the candidates to replace $j_{0}$, constraints (4.2) assures that in case $j$ is chosen for the replacement, and a customer $i$ is allocated to $j$, then the travel cost for $i$ is $d_{ij}$, otherwise, the constraint is redundant. With respect to the variables $\xi$ that model the selection of the facility to be swapped with $j_{0}$, (4.3) ensures that unchosen facilities cannot serve any customer. (4.4) assures that a single choice for $j_{0}$ is possible. Although (${\rm BestRepl}(j)$) is similar to (${\rm OMPN}_{BEP}$), the number of binary variables in the problemiss $(n-p)n$ instead of $n^{2}$.

In our experiments, we have checked that solving (${\rm BestRepl}(j)$) required more CPU time than solving the $n-p$ problems in the form (${\rm ALLOC}(J)$), although for problems in which $n-p\ll n$ (i.e., when $p$ is large), the compact formulation may consume less CPU time than loading and solving $n-p$ problems in the shape of (${\rm ALLOC}(J)$).

In what follows we describe our math-heuristic procedure, whose pseudocode is shown in Algorithm 1. Given an initial set of $p$ open facilities, it iterates by interchanging open facilities with other potential facilities trying to improve the best upper bound. At each iteration an open facility is selected to be replaced and the best replacement is chosen. After checking all the open facilities, if an improvement is found when compared to the best upper bound, the latest and the set of open facilities are updated. The procedure repeats the same scheme until a termination criterion is fulfilled. In our case, two stopping criteria are considered: maximum number of iterations and maximum number of iterations without improvement in the solution. In order to reduce the computation times required for solving (${\rm ALLOC}(J)$) or (${\rm BestRepl}(j)$), we consider a randomized version of the algorithm in which instead of finding best replacements for all the open facilities, a random one is selected at that phase of the approach.

A crucial point of local search heuristics is the quality of an initial feasible solution (in the $x$-variables). We compute the solution of the DOMP problem but using at costs between facilities $i$ and $j$ a convex combination of the lower and upper bounds $\widehat{d}_{ij}$ and $\widehat{D}_{ij}$ which provide good results in practice.

4.2. Math-Heuristic Algorithm 2

The second heuristic is based on alternating the location and allocation decisions. Initially, a DOMP is solved by fixing $\bar{a}=a$, and precomputing the distances between the facilities. Once a solution is obtained, the optimal open facilities are kept and given as input to (${\rm ALLOC}(J)$). Then, the variables $\bar{a}$ are updated with the obtained solution and the process is repeated until stabilization. In order to escape from local optima, the scheme is applied again but forbidding the use of one of the facilities opened in the first stage. The process iterates until no improvements are found.

The pseudocode for this approach is shown in Algorithm 2.

5. Experiments

In order to test the performance of the math-heuristic approaches, we have run some experiments over the real dataset instance of 2-dimensional coordinates (normalized longitude and latitude) of geographical centers of 49 states of the Unites States (we exclude Alaska, Hawaii and those outside Northamerica). We considered as neighborhoods Euclidean disks with radii based on the areas of each state. For each state (indexed by $j$), the area (in km2), $A_{j}$, was obtained and we construct the radius $r^{0}_{j}=\sqrt{\frac{A_{j}}{\pi}}$. The coordinates and the discs built applying this strategy are drawn in Figure 5. Then, three different scenarios were considered:

**S1: **

: $r_{j}=r^{0}_{j}$, for $j=1,\ldots 49$.

**S2: **

: $r_{j}=2\times r^{0}_{j}$, for $j=1,\ldots 49$.

**S3: **

: $r_{j}=3\times r^{0}_{j}$, for $j=1,\ldots 49$.

The interested reader may download the datasets at http://bit.ly/datasetUS.

We implemented in Gurobi under the C API the two math-heuristic approaches and we compare the running times and the best values obtained with these procedures and those obtained with the exact (${\rm OMPN}_{BEP}$) formulation (with a time limit of 1 hour). We solved the $p$-Median (M), $p$-Center (C), $p$-$25$-center (K) and $p$-Centdian (D) with $p\in\{2,5,10\}$. The results are reported in Table 6. In such a table, the first column indicates the scenario ($1$, $2$ or $3$), the second column (Pr.) shows the problem type and the third column indicates the number of facilities to be open, $p$. The values of the solutions obtained by using the different aproaches as well as their CPU running times (in seconds) are reported:

• •

Initial solution obtained by solving the nominal DOMP problem and solving (${\rm ALLOC}(J)$) for the obtained open facilities: H0 and t0.

• •

Best solution obtained by the math-heuristic approach 1: H1 and t1.

• •

Best solution obtained by the math-heuristic approach 2: H2 and t2.

• •

Best solution obtained by exact formulation (${\rm OMPN}_{BEP}$) within the time limit: BEP and tBEP.

We also report in the 12th column (%VarFixed) the average number of binary variables fixed in the first heuristic, and the pertentage deviations of the obtained solutions with respect to the best solution found with the exact formulation within the time limit: G1, G2 and G0 for the first heuristic, the second heuristic and the initial solution, respectively.

One can observe from the results that, the CPU times needed to run the math-heuristic approaches are much smaller than those needed to solve the OMPN problem with the MINLP formulation. In those cases in which all the approaches were able to solve the problem before the time limit of one hour, the highest deviation with respect to the optimal solutions was $15\%$ for the first heuristic and $3.2\%$ for the second one. In those cases in which the exact approach was not able to certify optimality in one hour, the math-heuristic approaches found a better solution for the problem. In the first heuristic, we apply the fixing variables strategy each time (${\rm ALLOC}(J)$) is solved. The average percentage of binary variables that are fixed with this strategy, is at least $84\%$ for scenario SC1, $75\%$ for SC2 and $52\%$ for SC3. Observe also that the initial solution based on fixing the open facilities to the solution of the DOMP problem and then compute the location on the neighborhoods and the allocation of the customers according to these positions, is in some case far of being a close-to-optimal choice, with percentage deviations of $33\%$ in some cases.

Note that the two math-heuristic approaches are still very time consuming. One may not forget that both proposed approaches are based on solving mixed integer non linear programming problems which are known to be NP-hard. The advantage of the two approaches is that they provided good quality solutions at the first iterations, which are competitive with the exact solutions (in terms of gap).

In Figure 6 we show the best solutions obtained for the test problem for $p=5$ under the center objective function for scenario SC1. The initial solution for this problem is drawn in Figure 7 (the solutions for $p=2,5,10$ can be found in bit.ly/resultsDOMPN). The reader can observe that small modification of the coordinates of the potential facilities (through neighborhoods) may produce different location-allocation solutions.

6. Conclusions

A unified version of the classical $p$-median problem is introduced in this paper which includes as particular cases the discrete and the continuous $p$-median and $p$-center problems. The problem considers that each facility can be located not only in the exact given position but in a neighborhood around it. Also, ordered median objective functions are modeled for the problem. Several mathematical programming formulations are proposed based on formulations for the discrete ordered median problem obtained from different sources. Two location-allocation approaches for solving the problem are presented. Although the optimization problems needed to solve are still NP-hard, the reduced dimensionality of them allows us to provide good quality solution in more reasonable times.

Several extensions are posible within this new framework. The first is the development of decomposition approaches for solving the OMPN. Lagrangean decomposition (relaxing the ordering constraints) combined with Benders decomposition (to separate the discrete and the continuous decisions) may produce exact solutions in better CPU times. On the other hand, although we analyze the ordered $p$-median problem, the results in this paper can be extended to other discrete location problems. For instance, capacitated [36] or multiperiod [1, 32] location problems can be embedded into the neighborhoods framework. Other interesting related problem which is left for further research is the consideration of location-routing problems with neighborhoods. That problem would involve not only the discrete facility location problem with neighborhoods but also the TSP with neighborhoods, then, the combination of the methods proposed in this paper with those provided in [16] may be applicable to the problem. Also, the case in which the neighborhood of each facility is the union of convex sets would be an interesting next step within this framework. In particular, it would model the case in which two facilities may belong to the same neighborhood. The extended MINLP formulations for such a problem will become disjunctive MINLP for which some techniques are available in the literature. Another approach that would extend the version introduced in this paper is the one in which $k_{j}$ facilities are allowed to be located at the $j$-th neighborhood to allocate the demand points. In such a case, a nested multifacility $p$-median problem is considered for which more sophisticated strategies should be developed to solve even small-size instances.

Acknowledgements

The author was partially supported by project MTM2016-74983-C2-1-R (MINECO, Spain), the research group SEJ-534 (Junta de Andalucía) and the research project PP2016-PIP06 (Universidad de Granada).