When Is Heterogeneous Distance-Decay Facility Location Tractable? A Structural Classification, Exact Methods, and a Real-World Study

Abstract: We study continuous planar facility location in which a demand point's captured value decays with distance, with the per-point decay scale varying across points. This heterogeneity is ubiquitous yet underexploited, and one nearest-facility objective unifies decay, clustering, and median goals, containing k-means, the Weber/p-median problem, and maximum covering as special cases. We make four contributions. (i) A tractability classification: the discrete objective is always monotone submodular, so the (1-1/e) greedy guarantee holds regardless of decay shape or heterogeneity, and the continuous cooperative objective is concave if and only if the decay is concave in distance; the clip max(0,d) in common coverage specifications is what destroys concavity, and the classification is tight. (ii) An exact discrete method: the candidate-discretized maximum-cover MIP has an empirically tight LP relaxation (~0% gap) and is solved by branch-and-bound in seconds for n <= 500. (iii) A force-as-gradient / large-neighborhood-search heuristic, within 0.5% of the discrete optimum, that outperforms the (1-1/e) greedy, Cooper-style alternating location-allocation, particle swarm optimization, and weighted k-means (30/30 per-instance wins at K=30, p<10^-9) and is competitive with bespoke solvers on k-means, Weber/p-median, and shape-demand instances. (iv) A real-world study: on 592,667 urban-delivery orders, ignoring the calibrated decay variation loses up to 9.7% of captured demand and relocates facilities by up to 37% of the map; a retail dataset calibrates the decay as exponential with scale R ~ 1.4 km.

• The paper establishes that discrete formulations are universally monotone submodular with a $(1-1/e)$ greedy guarantee, while continuous cooperative formulations are globally tractable exactly when all decay functions are concave in distance.
• The paper introduces FBM-LNS, combining force-based Lloyd ascent with large-neighborhood destroy-and-repair, which outperforms tested baselines by 1.5–3.5% and comes within 0.5% of fine-grid exact solutions.
• The paper shows that modeling heterogeneous decay scales is practically important: using a homogeneous mean scale reduced captured demand by up to 9.7% in a Shenzhen case study involving 592,667 delivery orders.

Problem and motivation

This paper studies continuous planar facility location in which each demand point's captured value decays with distance to its nearest facility, with a per-point decay scale $R_i$ that varies across points. The authors formulate the heterogeneous distance-decay maximum capture problem (PMC-HDD): place $p$ facilities in the plane to maximize

$$Z(\bm{X}) = \sum_{i} w_i \max_j \phi_i(\|p_i - X_j\|),$$

where $\phi_i$ is a non-increasing decay (linear, exponential, quadratic, or step) parameterized by scale $R_i$. The model deliberately unifies three canonical problems—maximum covering, $k$-means, and the Weber/$p$-median problem—as special cases of one nearest-facility objective distinguished only by the decay function. The motivating applications (emergency communication pre-positioning, ocean-monitoring platforms, wireless service planning) share two features: coverage is gradual rather than threshold-based, and distance sensitivity is heterogeneous, with dense urban cores demanding proximity while dispersed peripheries tolerate it.

The paper identifies five sources of difficulty: NP-hardness under step decay, non-convexity of the objective for $p \ge 2$, coupled continuous–combinatorial structure, an infinite candidate set with no known finite dominating set, and per-point heterogeneity. Each difficulty is addressed by a specific structural result.

Structural theory and the tractability classification

The theoretical core is a chain of results. The nearest-distance reduction shows that because $\phi_i$ is non-increasing, the max over facilities equals the decay evaluated at the nearest-facility distance, so allocation is always the Euclidean Voronoi partition regardless of decay form or heterogeneity. Convex-hull containment follows by metric projection: some optimal solution places every facility in $\operatorname{conv}\{p_i\}$.

The central classification theorem gives a sharp dichotomy:

• Discrete: after candidate discretization, both the cooperative ($p$0) and non-cooperative ($p$1) objectives are monotone submodular for any non-increasing decay, so the $p$2 greedy guarantee holds universally, improvable to $p$3 under bounded curvature.
• Continuous cooperative: $p$4 is concave in $p$5 if and only if every $p$6 is concave in the distance; the condition is tight, since a single non-concave $p$7 already yields non-concavity at $p$8, $p$9. Under concavity, every local maximum is global.
• Continuous non-cooperative: $$Z(\bm{X}) = \sum_{i} w_i \max_j \phi_i(\|p_i - X_j\|),$$0 is generically non-concave for $$Z(\bm{X}) = \sum_{i} w_i \max_j \phi_i(\|p_i - X_j\|),$$1 regardless of $$Z(\bm{X}) = \sum_{i} w_i \max_j \phi_i(\|p_i - X_j\|),$$2.

A practically important corollary is that the clip $$Z(\bm{X}) = \sum_{i} w_i \max_j \phi_i(\|p_i - X_j\|),$$3 present in most common coverage specifications destroys continuous concavity: clipped-linear decay $$Z(\bm{X}) = \sum_{i} w_i \max_j \phi_i(\|p_i - X_j\|),$$4 is convex in $$Z(\bm{X}) = \sum_{i} w_i \max_j \phi_i(\|p_i - X_j\|),$$5 (its slope jumps upward at the cutoff), so only unclipped decays such as $$Z(\bm{X}) = \sum_{i} w_i \max_j \phi_i(\|p_i - X_j\|),$$6 or $$Z(\bm{X}) = \sum_{i} w_i \max_j \phi_i(\|p_i - X_j\|),$$7 are trap-free. This is demonstrated numerically on a Gaussian-mixture instance where multi-start gradient ascent matches an SOCP global optimum when the clip is inactive but scatters across local optima (across-start standard deviation 137 versus 0.01) once the clip activates.

Two further results shape the algorithm. The Weber reduction shows that under linear decay, the within-cell location step is a weighted Weber problem with weights $$Z(\bm{X}) = \sum_{i} w_i \max_j \phi_i(\|p_i - X_j\|),$$8, so heterogeneity enters formally as stronger pull from distance-sensitive points, with a unique cell optimum away from collinear configurations. The gradient-as-force identity establishes that the demand-pull force on each facility equals exactly $$Z(\bm{X}) = \sum_{i} w_i \max_j \phi_i(\|p_i - X_j\|),$$9 wherever winners are unique—a measure-zero exception on Voronoi boundaries—making Lloyd's location step provably gradient ascent. Finally, Lloyd fixed points are characterized as coordinate-wise optima, and a monotone destroy-and-repair relocate that accepts only strict improvements terminates finitely at configurations admitting no improving single-facility re-placement.

The FBM-LNS algorithm

FBM-LNS (force-based metaheuristic with large-neighborhood search) combines the two pillars: multi-start weighted $\phi_i$0-means++ seeding followed by Lloyd gradient ascent, interleaved with a large-neighborhood relocate that removes the least-contributing facility and re-places it at the maximizer of the residual marginal gain $\phi_i$1—itself a weighted Weber problem on the residual active set. A small probability of accepting non-improving relocations escapes deep fixed points; the strict-improvement variant retains the monotone guarantee. The same code handles all decay families by switching one line, specializing to the closed-form centroid for $\phi_i$2-means and the Weiszfeld iteration for linear decay.

Exact discrete method and near-optimality certificates

Restricting facilities to a finite candidate set yields a $\phi_i$3-assignment max-cover MIP whose LP relaxation has worst-case gap $\phi_i$4 but is empirically tight. The paper attributes this tightness structurally: at Lloyd-stable configurations, facilities command distinct Voronoi cells with approximately disjoint coverage, so submodular curvature is low and LP, greedy, and integer optimum coincide. Branch-and-bound solves instances up to $\phi_i$5 in seconds, with reported integrality gaps of $\phi_i$6 across all tested configurations.

Two propositions give computable continuous upper bounds: a grid-density bound via Lipschitz rounding slack $\phi_i$7, valid for any finite grid but conservative, and a grid-free bound for concave decays computable in $\phi_i$8 but loose (it allows all $\phi_i$9 facilities to cover every point). Both bounds are conceded to be valid but not tight; the operative practical certificate is fine-grid convergence, which places FBM-LNS within $R_i$0 of the exact discrete optimum. A genuinely tight continuous relaxation would require handling the pointwise maximum directly—a convex maximization—and is left open.

Computational results

Experiments span seven test sets (concentration-point, uniform, Gaussian mixture, disk-clustered, Franti s1–s4, TSPLIB EUC_2D, region shapes), four decay families, and controlled heterogeneity via the coefficient of variation $R_i$1 of the decay scales. Six findings stand out:

1. Superiority on heterogeneous decay. FBM-LNS attains the highest mean captured demand in every configuration against $R_i$2-means++, Cooper ALT, PSO, and the $R_i$3 greedy, with margins of 1.5–3.5% over greedy. At $R_i$4 paired instances it wins 30/30 comparisons at $R_i$5; the paper notes honestly that unpaired confidence intervals overlap because between-instance variability dominates the within-instance algorithmic gap.
2. Cross-family competitiveness. One implementation ties $R_i$6-means++ on the Franti s1–s4 benchmark, beats specialized variable neighborhood search on larger $R_i$7-median instances (d657 by 1.2%, pr1002 by 1.8%) but loses on eil101 by 8.2%, and beats Cooper by roughly 4.7–5.7% on shape-demand instances reduced by sampling (seed-stable to within ~1%).
3. Heterogeneity amplifies the advantage. FBM-LNS's margin over Cooper grows from +2.4% at $R_i$8 to +6.8% at $R_i$9, consistent with the Weber-weight interpretation: dispersed $k$0 create larger residual marginals for the relocate to recover.
4. The relocate drives the gain. Multi-start Lloyd alone underperforms even the greedy; adding the relocate lifts performance by about 6.5%, while the selection rule (least-contributing, random, or mixture) is immaterial—the destroy-and-repair structure with maximum-marginal repair carries the improvement.
5. LP tightness and near-optimality. Integrality gaps are ≈0% even under adversarial high-overlap clustering, and FBM-LNS lies within 0.5% of the grid-MIP optimum at resolution 40, sometimes exceeding it because it optimizes continuously within the hull.
6. Scope boundary at binary decay. On step decay the advantage reverses: the greedy leads by 2–5% on clustered geographies at $k$1, $k$2. The explanation is first-principles—the step function's gradient vanishes almost everywhere, closing the continuous channel while the submodular-marginal channel remains open.

Real-world case study

Using data from the 2026 Meituan Business Analytics Elite Competition, the study aggregates 592,667 Shenzhen Nanshan delivery orders into 347 demand cells with density-derived decay scales spanning [79, 1953] m (coefficient of variation 1.33). A separate retail panel calibrates the decay form: exponential with $k$3 km fits better than linear (RMSE 0.12 vs. 0.17), consistent with the observed median trip length of 1.47 km.

Under this heterogeneity model, optimizing with the homogeneous mean scale instead of the heterogeneous $k$4 loses up to 9.7% of captured demand at $k$5 (shrinking to 2.2% at $k$6) and relocates facilities by 29–37% of the map width at every $k$7. The paper is explicit that the density-to-range assignment of $k$8 is a transparent modeling assumption not independently validated against the data, though its direction is supported by a −0.40 correlation between cell density and median trip distance. A Huff competitive extension shows captured share ranging from 37.6% against a weak incumbent to 7.8% against a strong one, with diminishing returns in $k$9 confirming submodularity empirically.

Limitations and open questions

The paper concedes several limitations at the points where they bear on results. FBM-LNS itself carries no global-optimality guarantee; its quality rests on empirical validation and the fine-grid certificate. Its advantage is specific to smooth decay, and the recommended tool switches to the greedy under binary coverage. The continuous upper bounds are loose, and a tight continuous relaxation remains open. The real-world $p$0 values derive from an unvalidated density proxy. Baseline comparisons use classical methods; modern problem-specific solvers such as MISOCP-based $p$1-median methods and dedicated Huff solvers are not included. A learned cross-entropy relocate selector improves performance only marginally (+0.22%, $p$2), which the authors attribute to the intrinsically low information content of per-facility scores given the ablation result—suggesting that progress requires policies acting on joint destroy-and-repair moves rather than per-facility features.

Conclusion

The paper delivers a tight structural classification for heterogeneous distance-decay facility location: universal discrete submodularity with the $p$3 guarantee, and continuous cooperative concavity if and only if the decay is concave in distance, with the clip mechanism identified as the practical destroyer of tractability. These results ground an algorithm whose every component maps to a proven property, an exact discretized method with empirically zero integrality gap, and quantified evidence—from synthetic benchmarks and 592,667 real orders—that exploiting decay heterogeneity matters materially, costing up to 9.7% of demand when ignored. The framework's scope is precisely delimited along the smooth/non-smooth boundary, making its applicability conditions explicit for practitioners.