Paper deep dive
When Is Heterogeneous Distance-Decay Facility Location Tractable? A Structural Classification, Exact Methods, and a Real-World Study
Zhou He, T. C. E. Cheng, Jichang Dong
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 91%
Last extracted: 7/21/2026, 2:34:24 AM
Summary
This paper investigates continuous planar facility location with heterogeneous distance-decay, where captured value decays with distance at varying scales per demand point. It establishes that the discrete objective is monotone submodular, ensuring a (1-1/e) greedy guarantee, while the continuous cooperative objective is concave if and only if the decay function is concave. The authors propose an exact mixed-integer programming method with tight LP relaxation and a force-based large-neighborhood search heuristic that outperforms existing methods like k-means and Cooper's alternating location-allocation. A real-world study on urban delivery orders demonstrates that ignoring decay heterogeneity can lead to significant losses in captured demand.
Entities (8)
Relation Signals (5)
Heterogeneous Distance-Decay Facility Location → hasproperty → Monotone Submodularity
confidence 95% · the discrete objective is always monotone submodular, so the (1-1/e) greedy guarantee holds
FBM-LNS Heuristic → outperforms → K-Means
confidence 90% · outperforms the ... weighted k-means (30/30 per-instance wins at K=30, p<10^-9)
Urban Delivery Orders → demonstrates → Heterogeneous Distance-Decay Facility Location
confidence 88% · on 592,667 urban-delivery orders ... ignoring the calibrated decay variation loses up to 9.7% of captured demand
Exponential Decay → calibratedin → Retail Dataset
confidence 85% · a retail dataset calibrates the decay as exponential with scale R ~ 1.4 km
FBM-LNS Heuristic → outperforms → Weber Problem
confidence 85% · competitive with problem-specific standard methods on ... Weber/p-median
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
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.
Tags
Links
- Source: https://arxiv.org/abs/2607.16764v1
- Canonical: https://arxiv.org/abs/2607.16764v1
Trouble viewing inline? Open PDF directly →
Full Text
118,076 characters extracted from source content.
Expand or collapse full text
When Is Heterogeneous Distance-Decay Facility Location Tractable? A Structural Classification, Exact Methods, and a Real-World Study Zhou He School of Economics and Management, University of Chinese Academy of Sciences, N606-3, #3 Zhongguancun Nanyitiao, Haidian District, Beijing 100190, China T. C. E. Cheng Department of Logistics and Maritime Studies, The Hong Kong Polytechnic University, 11 Yuk Choi Rd, Hung Hom, Hong Kong, China Jichang Dong Corresponding author. Email: jcdonglc@ucas.ac.cn School of Economics and Management, University of Chinese Academy of Sciences, N606-3, #3 Zhongguancun Nanyitiao, Haidian District, Beijing 100190, China 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 (dense urban cores demand proximity; dispersed peripheries tolerate distance) yet underexploited in continuous planar location models that unify decay, clustering, and median objectives, and the same model contains k-means, the Weber/p-median problem, and maximum covering as special cases of one nearest-facility objective. We make four contributions. (i) A tractability classification theorem: the discrete objective is always monotone submodular, so the (1−1/e)(1-1/e) greedy guarantee holds regardless of decay shape or heterogeneity; and the continuous cooperative objective is concave (trap-free) if and only if the decay is concave in the distance, identifying the clip max(0,⋅) (0,·)—present in most common coverage specifications—as a key mechanism that destroys continuous concavity; the classification is tight: concavity holds if and only if the decay is concave in distance. (i) An exact discrete method: the candidate-discretized maximum-cover mixed-integer program (MIP) has an empirically tight linear-programming (LP) relaxation (≈0%≈ 0\% gap) and is solved to optimality by branch-and-bound in seconds for n≤500n≤ 500. (i) A force-as-gradient / large-neighborhood-search heuristic, verified by fine-grid convergence to be within 0.5%0.5\% of the discrete optimum, that outperforms the (1−1/e)(1-1/e) greedy, Cooper-style alternating location–allocation, particle swarm optimization, and weighted k-means (30/3030/30 per-instance wins at K=30K=30, p<10−9p<10^-9) and is competitive with problem-specific standard methods on k-means, Weber/p-median, and shape-demand instances. (iv) A real-world study: on 592,667592,667 urban-delivery orders with density-derived heterogeneity, ignoring the calibrated decay variation loses up to 9.7%9.7\% of captured demand and relocates facilities by up to 37%37\% of the map; a separate retail dataset calibrates the decay form (exponential, R≈1.4R≈ 1.4 km). 1 Introduction Continuous planar facility location is a foundational problem family at the intersection of optimization and operations research: place p facilities anywhere in the plane so as to optimize a sum of per-point terms, each a function of the distance from a demand point to its nearest facility. Three canonical members, namely maximum capture/covering, k-means clustering, and the Weber (continuous p-median) problem, are all special cases of a single nearest-facility fidelity objective, distinguished only by the choice of the per-point function. Despite this shared structure, they are addressed by largely disjoint algorithmic literatures: Lloyd’s iteration and k-means++ for k-means, Cooper’s alternating location–allocation and the Weiszfeld iteration for the Weber problem, and (1−1/e)(1-1/e) greedy for covering. A natural question is whether a single algorithm, designed from the shared structure rather than tailored to one member, can be competitive with the problem-specific state of the art across the whole family, and more fundamentally, when the resulting continuous problem is tractable at all. The practical motivation comes from settings in which the per-point function is both gradual (coverage decays smoothly with distance, rather than at a sharp threshold) and heterogeneous (different demand points decline at different rates) (Drezner et al. 2004, Church and Murray 2018). In cellular network design, sensor placement, pre-positioning of emergency communication assets, and ocean-monitoring platform deployment, signal strength, sensing fidelity, and response quality all fall off with distance, and the tolerated attenuation differs across users and targets. A dense urban core (small decay scale RiR_i) demands a facility nearby; a scattered rural hamlet (large RiR_i) tolerates a more distant one, an urban–rural contrast in distance sensitivity that is empirically documented in mobility data (Yang et al. 2019). These two features, gradual and heterogeneous, define our main instance, heterogeneous distance-decay maximum capture, which is not directly handled by k-means++ or Cooper-style methods. 1.1 Motivating examples and decay functions Three decision contexts fit the model directly; each motivates one of the three decay families we study, summarized in Figure 1. 1. Pre-positioning emergency communication assets. After a disaster, each population cluster is a demand point whose weight is its population, and the service quality it receives degrades uniformly with the response distance up to a cutoff (a response-time target). The natural decay is the linear family ϕi(r)=[1−r/Ri]+ _i(r)=[1-r/R_i]^+: full coverage at the facility, linear decline to zero at r=Rir=R_i, and none beyond. 2. Ocean-monitoring platform deployment. Ecological hotspots and shipping lanes are demand points, and a monitoring buoy covers a target only if it lies within the buoy’s sensing range. The natural decay is the step family ϕi(r)=[r≤Ri] _i(r)=1[r≤ R_i]: full coverage inside the radius RiR_i and none outside—binary coverage, the limit that recovers the classical maximum covering location problem (MCLP). Heterogeneity reflects that some targets (protected habitat) require a closer RiR_i than others. 3. Wireless service-point planning. Received signal strength follows a path-loss law, and the tolerated attenuation differs between high-capacity business users and residential users. The natural decay is the exponential family ϕi(r)=e−r/Ri _i(r)=e^-r/R_i: a smooth, always-positive decline that also underlies the Huff gravity model of spatial choice. Figure 1 draws the three families against the normalized distance r/Rir/R_i. Heterogeneity across demand points (different RiR_i) is the defining feature of our setting: two points at the same distance from a facility can be captured at very different rates. r/Rir/R_iϕi(r) _i(r)11111.51.5linearexponentialstep linear ϕi(r)=[1−r/Ri]+ _i(r)=[1-r/R_i]^+ exponential ϕi(r)=e−r/Ri _i(r)=e^-r/R_i step ϕi(r)=[r≤Ri] _i(r)=1[r≤ R_i] Figure 1: Three decay families used in this paper, drawn against the normalized distance r/Rir/R_i, with their function forms listed on the right. Different line styles distinguish the families; the same distance can yield very different captured demand depending on the family and on the per-point scale RiR_i. 1.2 Related work Gradual covering was introduced on networks by Berman et al. (2003) and in the continuous plane by Drezner et al. (2004), who allow per-demand-point coverage parameters; partial-coverage extensions include Karasakal and Karasakal (2004) and the ordered gradual covering model of Berman et al. (2009). The multiple gradual cover location problem (MGCLP) of Álvarez-Miranda and Sinnl (2019) provides exact mixed-integer formulations that exploit the submodularity of the cooperative objective. Most directly related to our setting is the planar model of Bansal and Shojaee (2020), which builds on Bansal and Kianfar (2017) and gives a (1−1/e)(1-1/e) greedy approximation and exact branch-and-bound for planar maximum coverage with partial coverage and adjustable quality of service. Church and Murray (2018) survey the broader covering-location literature. The alternating location–allocation (ALT) heuristic of Cooper (1964), equivalent to Lloyd’s k-means algorithm (Lloyd 1982), is the workhorse local optimizer for planar median and covering problems: it alternates between allocating each demand point to its nearest facility and locating each facility at the optimum (centroid for k-means, Weber median for p-median) of its assigned cell, repeating to a coordinate-wise fixed point. Its convergence behavior is analyzed in Ostrovsky et al. (2013) and, in a statistical mixture framework, in Lu and Zhou (2016). For k-means, the dominant practical method is k-means++ seeding (Arthur and Vassilvitskii 2007), which places initial facilities with probability proportional to squared distance to existing facilities and then runs Lloyd; for the Weber/p-median problem, variable neighborhood search (Hansen and Mladenović 1997) and the Weiszfeld iteration (Weiszfeld 1937)—a fixed-point iteration for the weighted geometric median—are standard. These literatures have developed in parallel rather than as instances of a single method. Submodular maximization under a cardinality constraint admits the classical (1−1/e)(1-1/e) greedy guarantee (Nemhauser et al. 1978)—add facilities one at a time, each at the point of maximum marginal gain—which is tight under standard complexity assumptions (Feige 1998) and improves to (1−e−κ)/κ(1-e^-κ)/κ under bounded curvature κ (Conforti and Cornuéjols 1984). These four building blocks—alternating location–allocation, k-means++ seeding, the Weiszfeld iteration, and the marginal greedy—are exactly the ingredients our algorithm recombines in §4. The Huff gravity model (Huff 1964) and competitive location (Drezner 1994) study market-share splitting. Throughout we treat the nearest-facility (max) model as the non-cooperative variant of gradual covering, complementary to the cooperative (sum) MGCLP of Álvarez-Miranda and Sinnl (2019). 1.3 Contributions This paper makes four contributions, organized around the question: when is heterogeneous distance-decay location tractable, and how should one exploit that structure? 1. A tractability classification theorem (§3). We prove that the discrete (candidate-set) objective is monotone submodular for any non-increasing decay, so the classical (1−1/e)(1-1/e) greedy guarantee (Nemhauser et al. 1978) and its curvature refinement (Conforti and Cornuéjols 1984) hold regardless of the decay shape or heterogeneity—extending the discrete guarantee underlying Bansal and Shojaee (2020) beyond the linear/partial-coverage case. For the continuous cooperative problem we give a sharp dichotomy: the objective is concave (free of local-optima traps, solvable by gradient ascent to the global optimum) if and only if the decay is concave in the distance. The clip max(0,⋅) (0,·) present in most common coverage specifications is a key mechanism that destroys continuous tractability; the unclipped Weber/k-means limits are the tractable cases, while convex decays such as exponential are also non-concave. This tells a practitioner when convex optimization suffices and when a heuristic is unavoidable. 2. An exact discrete method (§5). The candidate-discretized problem, formulated as a y-assignment max-cover mixed-integer program (MIP), has an empirically tight linear-programming (LP) relaxation (≈0%≈ 0\% integrality gap across all tested decay families and scales, including adversarial high-overlap configurations), and branch-and-bound solves it to optimality in seconds for n≤500n≤ 500. Whereas Álvarez-Miranda and Sinnl (2019) give exact formulations for the cooperative MGCLP and Bansal and Shojaee (2020) use branch-and-bound for planar partial coverage, we exploit the LP tightness mainly to certify our heuristic, not merely to solve small instances. 3. A heuristic with a near-optimality certificate (§4, §6). Our force-based metaheuristic with large-neighborhood search (FBM-LNS)—multi-start Lloyd gradient ascent interleaved with submodular-marginal large-neighborhood relocate—is verified by fine-grid convergence to be within 0.5%0.5\% of the discrete optimum. Unlike Cooper’s ALT (Cooper 1964), which converges to a coordinate-wise fixed point, our relocate escapes such fixed points; and unlike the (1−1/e)(1-1/e) greedy it jointly re-optimizes every facility. It significantly outperforms the greedy, Cooper ALT, particle swarm optimization, and weighted k-means on the heterogeneous-decay problem (30/3030/30 per-instance wins at K=30K=30, p<10−9p<10^-9), while remaining competitive with k-means++ (Arthur and Vassilvitskii 2007) on k-means and with Cooper on Weber/p-median and shape-demand instances. 4. A real-world case study (§7). On 592,667592,667 urban-delivery orders with density-derived heterogeneity, ignoring the calibrated decay variation loses up to 9.7%9.7\% of captured demand and relocates facilities by up to 37%37\% of the map. A separate retail dataset calibrates the decay form (exponential, R≈1.4R≈ 1.4 km), and a competitive (Huff) extension (Huff 1964, Drezner 1994) quantifies captured market share against an incumbent. Both real datasets are drawn from the 2026 6th Meituan Business Analytics Elite Competition. The remainder of the paper is organized as follows. Section 2 formalizes the problem. Section 3 develops the structural theory and the tractability classification. Section 4 presents the algorithm. Section 5 develops the exact discrete method. Section 6 reports computational experiments. Section 7 presents the real-world case study. Section 8 discusses limitations, extensions, and usage guidance, and Section 9 concludes. 2 Problem formulation 2.1 Index, sets, and model We formulate the heterogeneous distance-decay maximum capture problem (PMC-HDD) in the standard notation of location science. Sets and indices. • I=1,…,nI=\1,…,n\: set of demand points (index i), with known location pi∈ℝ2p_i ^2. • J=1,…,pJ=\1,…,p\: set of facilities to be located (index j). Parameters. • wi∈ℝ++w_i _++: demand weight (e.g., population, traffic volume, order count) at point i. • ϕi:ℝ+→[0,1] _i:R_+→[0,1]: non-increasing decay function for point i, with ϕi(0)=1 _i(0)=1, parameterized by a scale Ri>0R_i>0. The families studied (Figure 1) are linear ϕi(r)=[1−r/Ri]+ _i(r)=[1-r/R_i]^+, exponential ϕi(r)=e−r/Ri _i(r)=e^-r/R_i, quadratic ϕi(r)=[1−(r/Ri)2]+ _i(r)=[1-(r/R_i)^2]^+, and step ϕi=[r≤Ri] _i=1[r≤ R_i]. A point with small RiR_i is highly distance-sensitive; the distribution of Ri\R_i\ encodes the heterogeneity of the instance. • p∈ℕp : the number of facilities (a cardinality budget). Decision variables. • Xj∈ℝ2X_j ^2 for each j∈Jj∈ J: the location of facility j (continuous; the candidate set is the whole plane). Objective and constraints. Each demand point i is served by the facility that captures the most of its demand (the nearest-facility or max model; Figure 2 illustrates a small instance, its Voronoi allocation, and points left out of range). The total captured demand is Z()=∑i∈Iwimaxj∈Jϕi(‖pi−Xj‖),Z( X)\;=\; _i∈ Iw_i\, _j∈ J\, _i\! (\|p_i-X_j\| ), (1) where =(X1,…,Xp) X=(X_1,…,X_p). By Theorem 3.2 (§3), there exists an optimal solution with every facility in convpii∈Iconv\p_i\_i∈ I, so we may restrict Xj∈convpii∈I,j∈J.X_j \p_i\_i∈ I, j∈ J. (2) We write Ψi()=wimaxjϕi(‖pi−Xj‖) _i( X)=w_i _j _i(\|p_i-X_j\|) for the captured demand at point i and π(i)=argminj‖pi−Xj‖π(i)= _j\|p_i-X_j\| for its serving (winning) facility. The max aggregation in (1) is the non-cooperative model: each demand point patronizes a single (nearest) facility. The tractability theory of §3 also covers the cooperative/Huff family, in which the max is replaced by a concave increasing saturation function g applied to the total attractiveness Ui=∑jϕi(‖pi−Xj‖)U_i= _j _i(\|p_i-X_j\|): Zg()=∑i∈Iwig(Ui()).Z^g( X)\;=\; _i∈ Iw_i\,g\! (U_i( X) ). (3) This subsumes two special cases. Cooperative gradual coverage is the MGCLP of Álvarez-Miranda and Sinnl (2019), recovered by g(u)=min(1,u)g(u)= (1,u); Huff competitive market-share splitting is g(u)=u/(u+Ci)g(u)=u/(u+C_i) with incumbent attractiveness Ci>0C_i>0 (Huff 1964, Drezner 1994). Both saturations are concave increasing in u, which is the only property the theory requires. 2.2 Variants and applications The model of §2 unifies several well-studied continuous location problems as special cases of the choice of decay function ϕi _i. We label each variant (V1–V6) for cross-reference in the computational study and give a concrete decision context for each. V1. Heterogeneous gradual coverage (linear): ϕi(r)=[1−r/Ri]+ _i(r)=[1-r/R_i]^+ with heterogeneous Ri\R_i\—our primary focus. A representative decision context is pre-positioning emergency communication assets after a disaster: each population cluster is a demand point whose weight wiw_i is its population and whose decay scale RiR_i reflects latency sensitivity, so a dense urban core (small RiR_i) demands a base station within minutes while a scattered rural hamlet (large RiR_i) tolerates a more distant one. V2. Heterogeneous gradual coverage (exponential): ϕi(r)=e−r/Ri _i(r)=e^-r/R_i with heterogeneous Ri\R_i\. This is the natural model for wireless service-point and cellular planning: received signal strength follows an exponential path-loss law, and the Huff gravity model of spatial choice assigns each customer a probability of patronizing a facility that falls off exponentially with distance. Heterogeneous Ri\R_i\ captures that high-capacity business users tolerate more attenuation than residential users. V3. k-means: ϕi(r)=1−(r/R)2 _i(r)=1-(r/R)^2 with large homogeneous R; equivalent to minimizing ∑iwiminj‖pi−Xj‖2 _iw_i _j\|p_i-X_j\|^2. The canonical application is service-region territorialization, e.g., partitioning a city’s delivery or service zones into compact territories around p depots to minimize the sum of squared trip lengths (the standard clustering objective that underlies the Franti s1--s4 benchmark used in §6). V4. Weber / p-median: ϕi(r)=1−r/R _i(r)=1-r/R with large R; minimizes ∑iwiminj‖pi−Xj‖ _iw_i _j\|p_i-X_j\|. This is the classical warehouse and distribution-center location problem on real geographic coordinates: the TSPLIB EUC_2D instances of §6 give actual city coordinates, and p facilities are placed to minimize total transport distance. V5. Shape demand: each demand unit is a region DiD_i with density μi _i; captured demand is the integral ∫Diϕi(minj‖p−Xj‖)μi(p)p _D_i _i( _j\|p-X_j\|)\, _i(p)\,dp, reducible to the point problem by sampling. A motivating example is ocean-monitoring platform deployment, where ecological hotspots, reefs, and shipping lanes are two-dimensional regions rather than points: a monitoring buoy’s sensing fidelity decays with distance, and the relevant coverage is the integral of fidelity over the protected habitat. Each region is sampled with a fixed random seed, so a given instance is exactly reproducible; the sampled objective converges to the integral as the sample count m grows (Monte-Carlo error O(1/m)O(1/ m)), and the resulting facility placement is stable to the sampling seed (we verify this in §6). V6. Binary covering (step): ϕi(r)=[r≤Ri] _i(r)=1[r≤ R_i], the classical maximum covering location problem (MCLP) of ambulance and sensor siting—full coverage inside the radius RiR_i and none outside. Its non-smooth objective removes the gradient signal our method exploits, so it serves as a scope boundary (§8.1) where the (1−1/e)(1-1/e) greedy is competitive. Variants V1–V5 are the main computational focus of §6; V6 (binary step) is tested as a scope boundary in §8.1. Our force-based heuristic, FBM-LNS (§4), handles all six variants by switching one line of code (the per-point function). 2.3 Complexity and sources of difficulty PMC-HDD is challenging for five distinct reasons, each labeled (D1–D5) for cross-reference with the algorithm design (§4). D1. NP-hardness. Under step decay ϕi=[r≤Ri] _i=1[r≤ R_i] the problem reduces to the planar MCLP, NP-hard when p is part of the input (Bansal and Kianfar 2017). Existing exact methods (Bansal and Shojaee 2020, Álvarez-Miranda and Sinnl 2019) use branch-and-bound on discretized candidates but do not scale beyond a few hundred points. The (1−1/e)(1-1/e) greedy (Nemhauser et al. 1978) provides a polynomial approximation but leaves a quality gap. D2. Non-convex objective. For p=1p=1 the within-cell objective is concave on its active range (Theorem 3.4); for p≥2p≥ 2 the max over facilities introduces further non-convexity with many local optima. Cooper (1964) proposed the alternating location–allocation heuristic to handle this, but it converges to a coordinate-wise fixed point and provides no global guarantee. D3. Coupled continuous–combinatorial structure. Locations are continuous but the nearest-facility assignment is combinatorial; each depends on the other. Cooper’s alternation decouples them temporarily but converges to fixed points that retain the coupling. D4. Infinite candidate set. The continuous plane admits infinitely many locations. On networks, a finite dominating set (FDS) exists for many objectives (Berman et al. 2009), enabling exact discrete optimization. No FDS is known for the continuous heterogeneous case, so the problem cannot be reduced to a finite candidate set a priori. D5. Per-point heterogeneity. Each RiR_i may differ, so the within-cell subproblem varies across demand points. Existing gradual-covering methods (Drezner et al. 2004) allow per-point parameters but treat them as fixed input rather than exploiting the structural implications (weights wi/Riw_i/R_i, Theorem 3.4) for algorithm design. These five difficulties motivate the structural theory (§3) and the algorithm (§4), each of whose components is designed to address a specific difficulty (Table 2). F1F2ABCDEFGH Voronoi edge Allocation (nearest facility) Facility Points Status F1 A, B, C captured F2 D, E, F captured F1 G out of range F2 H out of range Red = sensitive (small RiR_i); blue = tolerant (large RiR_i); circle area ∝wi w_i. Figure 2: A heterogeneous maximum-capture instance with eight labeled demand points (A–H) and two facilities (F1, F2). Each point is allocated to its nearest facility by the Voronoi edge (the perpendicular bisector of the F1–F2 segment, dashed); points G and H lie in their cell but outside the facility’s coverage disk, so their demand is not captured. Points are colored by sensitivity (red ≡ small RiR_i, blue ≡ large RiR_i) and sized by weight. 3 Structural theory and the tractability classification This section develops the structural properties on which the algorithm and the tractability classification rest. The proofs given here are brief sketches that emphasize the idea; the detailed, self-contained versions are collected in the Electronic Companion. The results form a single logical chain: the nearest-distance reduction fixes the Voronoi allocation and resolves the continuous–combinatorial coupling at the allocation level, convex-hull containment bounds the search set, submodularity provides the discrete approximation guarantee, the Weber reduction shapes the within-cell landscape, the gradient-as-force identity drives the local search, the tractability classification pins down when the continuous problem is trap-free, and the two fixed-point propositions define where the local search stops and how the relocate restarts it. 3.1 Nearest-distance reduction and Voronoi allocation Because each ϕi _i is non-increasing, the per-point maximum over facilities equals the decay evaluated at the nearest-facility distance. Theorem 3.1 (Nearest-distance reduction). Z()=∑i∈Iwiϕi(minj∈J‖pi−Xj‖)Z( X)= _i∈ Iw_i\, _i\! ( _j∈ J\|p_i-X_j\| ), and the optimal allocation satisfies π(i)=argminj‖pi−Xj‖π(i)= _j\|p_i-X_j\| for every i. Proof maxjϕi(‖pi−Xj‖)=ϕi(minj‖pi−Xj‖) _j _i(\|p_i-X_j\|)= _i( _j\|p_i-X_j\|) since ϕi _i is non-increasing, and the argmax coincides with the argmin of the distance. □ □ Thus the allocation is the standard Euclidean Voronoi partition regardless of the decay form or heterogeneity. The decay parameters affect only the within-cell location objective and the captured value, not the combinatorial structure of the partition. Theorem 3.1 is the foundation of the paper: it reduces the max-over-facilities objective to a function of nearest-facility distances, fixing the Voronoi allocation once X is given (difficulty D3, cf. Table 2). Its sole hypothesis is the monotonicity intrinsic to any decay function, and it is the prerequisite for the Weber reduction (Theorem 3.4), the gradient-as-force identity (Theorem 3.7), and the Lloyd step (§4). It also makes precise that k-means, the Weber problem, and maximum covering are instances of one objective. 3.2 Convex hull containment Theorem 3.2. There exists an optimal solution with every facility in convpii∈Iconv\p_i\_i∈ I. Proof If XjX_j lies outside the closed convex hull, its metric projection Xj′X_j onto the hull is no farther from every pip_i; since each ϕi _i is non-increasing, coverage does not decrease. □ □ An optimizer therefore lies in the compact set convpiconv\p_i\, so the algorithm’s gradient search is clipped to it (§4; difficulty D4, cf. Table 2). The hypothesis is the same monotonicity as Theorem 3.1, and the projection argument is dimension-independent. 3.3 Submodularity, curvature, and approximation Discretizing candidate locations to a finite set C yields the set objective ZD(S)=∑iwimaxc∈Sϕi(‖pi−c‖)Z_D(S)= _iw_i _c∈ S _i(\|p_i-c\|). Theorem 3.3 (Monotone submodularity). ZDZ_D is monotone non-decreasing and submodular on the lattice of subsets of C. Proof For S⊆TS T and candidate c∉Tc∉ T, the marginal gain Δi(S)=wi[ϕi(‖pi−c‖)−maxc′∈Sϕi(‖pi−c′‖)]+ _i(S)=w_i[ _i(\|p_i-c\|)- _c ∈ S _i(\|p_i-c \|)]^+ is non-increasing in the set because maxc′∈Sϕi(‖pi−c′‖)≤maxc′∈Tϕi(‖pi−c′‖) _c ∈ S _i(\|p_i-c \|)≤ _c ∈ T _i(\|p_i-c \|). Summing preserves the inequality. □ □ The Nemhauser–Wolsey–Fisher greedy therefore achieves ZD(Sgreedy)≥(1−1/e)ZD∗Z_D(S^greedy)≥(1-1/e)Z_D^* (Nemhauser et al. 1978), a bound that is tight in the worst case (Feige 1998). Bounded submodular curvature κ∈(0,1]κ∈(0,1] (Conforti and Cornuéjols 1984) improves this to the data-dependent ratio ZD(Sgreedy)≥1−e−κZD∗>(1−1/e)ZD∗,Z_D(S^greedy)\;≥\; 1-e^-κ\,Z_D^*\;>\;(1-1/e)\,Z_D^*, (4) Intuitively, when facility service areas overlap little (low redundancy), κ is small and the greedy ratio exceeds 1−1/e1-1/e. This provides the formal (1−1/e)(1-1/e) quality floor against which the heuristic is measured, and it underpins the Greedy baseline and the LNS repair step (difficulty D1, cf. Table 2). The only extra ingredient beyond monotonicity is discretization to a finite C—a modeling convenience, since the continuous plane is recovered as C becomes dense. 3.4 The within-cell location step is a weighted Weber problem Fixing a Voronoi cell S⊆IS I and the linear decay family, the single-facility problem for the cell is maxX∑i∈Swi(1−‖pi−X‖Ri) _X _i∈ Sw_i (1- \|p_i-X\|R_i ) ≡minX∑i∈SwiRi‖pi−X‖. ≡\; _X _i∈ S w_iR_i\,\|p_i-X\|. Theorem 3.4 (Weber reduction). Under linear decay, the within-cell location step of Lloyd’s algorithm is a weighted Weber problem with weights αi=wi/Ri _i=w_i/R_i. Proof For a fixed cell S and linear ϕi _i, the contribution of point i is wi(1−‖pi−X‖/Ri)w_i(1-\|p_i-X\|/R_i). Dropping the additive constant ∑i∈Swi _i∈ Sw_i, maximizing ∑i∈Swi(1−‖pi−X‖/Ri) _i∈ Sw_i(1-\|p_i-X\|/R_i) over X is equivalent to minimizing ∑i∈S(wi/Ri)‖pi−X‖ _i∈ S(w_i/R_i)\|p_i-X\|, the weighted Weber objective with weights αi=wi/Ri _i=w_i/R_i. □ □ This connects our location step directly to the classical Weber/Fermat literature and its linearly convergent Weiszfeld iteration (Weiszfeld 1937). Two consequences follow immediately. Corollary 3.5 (Heterogeneity, formally). Distance-sensitive points (small RiR_i) receive larger weight αi=wi/Ri _i=w_i/R_i and therefore pull the facility more strongly; heterogeneous sensitivity is a first-order consequence of the Weber weights rather than an algorithmic device. Proof The first-order optimality condition of the weighted Weber problem is ∑i∈Sαi(X−pi)/‖X−pi‖= _i∈ S _i\,(X-p_i)/\|X-p_i\|=0 with αi=wi/Ri _i=w_i/R_i; points with smaller RiR_i have larger αi _i and hence proportionally larger pull on X. □ □ Corollary 3.6 (Uniqueness). In two dimensions the weighted Weber objective is convex and its minimizer is unique unless all in-cell demand points and the iterate are collinear. Proof X↦‖pi−X‖X \|p_i-X\| is convex, so the non-negative combination ∑iαi‖pi−X‖ _i _i\|p_i-X\| is convex; strict convexity of the Euclidean norm in ℝ2R^2 (away from collinear configurations) yields a unique minimizer. □ □ Heterogeneity thus enters the within-cell problem as the Weber weight αi=wi/Ri _i=w_i/R_i (difficulty D5, cf. Table 2): distance-sensitive points pull harder automatically (Corollary 3.5) and the cell optimum is unique (Corollary 3.6). The Weber reduction is the one substantive hypothesis in this group: it requires the linear decay family, since only then does wi(1−‖pi−X‖/Ri)w_i(1-\|p_i-X\|/R_i) split into a constant plus a weighted distance. For exponential or quadratic decay the within-cell problem is still smooth but is no longer Weber in closed form, and we solve it by gradient ascent (§4). 3.5 Gradient-as-force The gradient of (1) with respect to XjX_j, restricted to the points for which j is the winner, is ∇XjZ()=∑i:π(i)=jwi∇Xjϕi(‖pi−Xj‖), _X_jZ( X)\;=\; _i:π(i)=jw_i\, _X_j _i(\|p_i-X_j\|), (5) which is exactly the “demand pull” force on facility j. Theorem 3.7 (Gradient-as-force). At any X at which each point’s winning facility is unique, ∇XjZ _X_jZ equals the right-hand side of (5); the demand-pull force on facility j is exactly the objective gradient. Proof Z=∑iwimaxkψikZ= _iw_i _k _ik with ψik=ϕi(‖pi−Xk‖) _ik= _i(\|p_i-X_k\|) is a sum of pointwise maxima of C1C^1 functions. Where the winner π(i)π(i) is unique, ∇Xjmaxkψik _X_j _k _ik equals ∇Xjψi,π(i) _X_j _i,π(i) if j=π(i)j=π(i) and 0 otherwise; summing over i yields (5). □ □ Theorem 3.7 is the bridge from analysis to algorithm: the “demand pull” on each facility—the physical intuition behind force-based local search—is exactly ∇XjZ _X_jZ, so Lloyd’s location step is gradient ascent rather than an ad hoc move. Differentiability of the pointwise maximum requires a unique winner at each demand point, which fails only on Voronoi cell boundaries—a set of measure zero—and is therefore the generic case. This justifies the “force-based” name of FBM-LNS (§4). It relies on the winner-only assignment of Theorem 3.1 and on the within-cell Weber structure of Theorem 3.4; the next two propositions (Propositions 3.9–3.10) characterize where this gradient ascent terminates and how the relocate step escapes it. Figure 3 visualizes the interpretation. XjX_jABCD∇XjZ _X_jZ Gradient computation Pt did_i wiw_i RiR_i wi/Riw_i/R_i A 2.63 4 2.0 2.00 B 2.48 3 2.0 1.50 C 2.58 5 5.0 1.00 D 2.64 4 5.0 0.80 ∇XjZ=∑iwiRi^(pi−Xj) _X_jZ= _i w_iR_i\, u(p_i-X_j) ≈(−1.11, 0.45)≈(-1.11,\,0.45) ‖∇XjZ‖≈1.20\| _X_jZ\|≈ 1.20 Red = sensitive (small RiR_i); blue = tolerant (large RiR_i). Figure 3: The demand pull on facility XjX_j is exactly the gradient ∇XjZ _X_jZ (Equation 5), illustrated with a concrete four-point cell. Each demand point pulls along the line toward itself with magnitude wi/Riw_i/R_i (red sensitive points A, B pull harder than blue tolerant points C, D); the facility moves along the resultant (bold arrow), which is the gradient-ascent direction. The right panel reports each point’s distance, weight, decay scale, and pull magnitude, and sums them to the gradient value ∇XjZ≈(−1.11, 0.45) _X_jZ≈(-1.11,\,0.45). This is the Lloyd location step. 3.6 The tractability classification theorem We now state the central structural result: a complete classification of when the heterogeneous distance-decay location problem is tractable. Theorem 3.8 (Tractability classification). Consider the cooperative objective Zg()=∑iwig(Ui)Z^g( X)= _iw_i\,g(U_i) with Ui=∑jϕi(‖pi−Xj‖)U_i= _j _i(\|p_i-X_j\|), g concave increasing, wi≥0w_i≥ 0. (i) (Discrete) ZDgZ^g_D and ZDmaxZ _D are both monotone submodular for any non-increasing ϕi _i. The greedy (1−1/e)(1-1/e) guarantee holds universally. (i) (Continuous, cooperative) ZgZ^g is concave in X whenever every ϕi _i is concave on its realized range; this condition is tight, since non-concavity of any ϕi _i already yields instances (already at p=1p=1, n=1n=1) on which ZgZ^g is non-concave. Under concavity, every local maximum is global; projected gradient ascent converges. (i) (Continuous, non-cooperative) ZmaxZ is generically non-concave for p≥2p≥ 2, regardless of ϕφ. Proof sketch. (i) The non-cooperative case is Theorem 3.3. For the cooperative case, concavity of g implies the diminishing-increment property that g(σ+a)−g(σ)g(σ+a)-g(σ) is non-increasing in σ (a standard fact we state and prove as Lemma EC.1 in the Electronic Companion); since σi(S)=∑c∈Sϕi(‖pi−c‖) _i(S)= _c∈ S _i(\|p_i-c\|) is modular and non-decreasing in S, the per-point cooperative marginal is non-increasing in S, and summing gives submodularity. (i, ⇐ ) The map X↦‖pi−X‖X \|p_i-X\| is convex; by the concave composition rule (Lemma EC.2 in the Electronic Companion), a concave non-increasing ϕi _i composed with this convex map is concave; Ui=∑jϕi(‖pi−Xj‖)U_i= _j _i(\|p_i-X_j\|) is then concave in X, as is g∘Uig U_i for concave non-decreasing g, and the non-negative weighted sum ZgZ^g is concave. (i, ⇒ ) Contrapositive: if ϕi0 _i_0 is not concave, take p=1p=1, n=1n=1, g=idg=id; Z(X)=ϕi0(‖pi0−X‖)Z(X)= _i_0(\|p_i_0-X\|), which along a ray X=pi0+tuX=p_i_0+t\,u equals ϕi0(t) _i_0(t) and is not concave. (i) One demand point, two facilities: Z=max(ϕ(‖p−X1‖),ϕ(‖p−X2‖))Z= (φ(\|p-X_1\|),φ(\|p-X_2\|)) has a kink where the winner switches, and is non-concave for any non-constant ϕφ (full proof in the Electronic Companion). □ □ In concrete terms, ϕi(r)=[1−r/Ri]+ _i(r)=[1-r/R_i]^+ is convex in r (it is max(0,⋅) (0,·) of an affine function, so its slope jumps upward at the cutoff), so the cooperative objective built from it has local-optima traps. Only unclipped decays (1−r/R1-r/R, 1−(r/R)21-(r/R)^2) or genuinely concave decays are trap-free. Table 1 summarizes. Table 1: Tractability classification by problem variant (all V1–V6). “Concave” = exact by gradient ascent; “submodular” = greedy (1−1/e)(1-1/e) on the discretized problem. V5 (shape demand) reduces to the point case by sampling, so its cooperative concavity inherits the chosen decay family (“inherits”); V6 is the binary scope boundary of §8.1. Cooperative (ZgZ^g) Non-coop. (ZmaxZ ) Variant Decay ϕ(r)φ(r) concave? submodular? concave? submodular? V1 [1−r/R]+[1-r/R]^+ (clipped linear) no yes no yes V2 e−r/Re^-r/R (exponential) no yes no yes V3 1−(r/R)21-(r/R)^2 (k-means) yes yes no yes V4 1−r/R1-r/R (Weber limit) yes yes no yes V5 regions, sampled (any ϕφ) inherits yes no yes V6 [r≤R]1[r≤R] (step) degenerate yes no yes The dichotomy is not merely formal. On a cooperative Gaussian-mixture instance (n=200n=200, p=5p=5) with the identity saturation g(u)=ug(u)=u (the pure cooperative sum, with no saturation cap) we solved the concave-regime cooperative objective exactly by a second-order cone program (SOCP) and ran multi-start gradient ascent from 88 random seeds: with a large scale R=5000R=5000 (clip inactive, so ϕφ is concave) every start converges to the same value (across-start standard deviation 0.010.01), matching the SOCP global optimum to within 0.020.02; with the small scale R=120R=120 (clip active, so ϕφ is convex) the starts end at different local optima (standard deviation 137137). The transition from “no traps” to “traps” as the clip switches on is exactly what Theorem 3.8(i) predicts. 3.7 Lloyd fixed points and the marginal relocate Lloyd’s algorithm (coordinate ascent on the allocation–location decomposition) converges to a coordinate-wise fixed point at which no single facility can be moved to improve the objective given the others’ locations and the Voronoi assignment. Proposition 3.9 (Lloyd fixed point). At a fixed point ⋆ X of the alternating allocation–location map, (i) the Voronoi assignment is stable (no point is reassigned by a strictly nearer facility) and (i) each facility lies at the within-cell Weber optimum of its cell; equivalently 0∈∂CZ0∈ _CZ with respect to each single-facility coordinate. Hence ⋆ X is coordinate-wise locally optimal. Proof At a Lloyd fixed point the location step makes no update, so each Xj⋆∈argmaxX∑i:π(i)=jwiϕi(‖pi−X‖)X_j ∈ _X _i:π(i)=jw_i _i(\|p_i-X\|), giving 0∈∂C,XjZ0∈ _C,X_jZ; the allocation step makes no update, so the Voronoi assignment is stable. Together these are coordinate-wise optimality. □ □ This coordinate-wise fixed-point limit is why pure multi-start Lloyd can be outperformed by the (1−1/e)(1-1/e) greedy—it converges to a coordinate-wise optimum that need not be globally good—and it motivates the destroy-and-repair relocate that follows. The gain from adding a facility at x given current captured demand cic_i is g(x)=∑iwi[ϕi(‖pi−x‖)−ci/wi]+g(x)= _iw_i[ _i(\|p_i-x\|)-c_i/w_i]^+, which on its residual active set is again a weighted Weber objective. Removing the least-contributing facility and re-placing it at a maximizer of g is a destroy(1)–repair(1) large-neighborhood move; accepting only improving moves yields a monotone ascent on Z. Proposition 3.10 (Monotone relocate). If each relocate step is accepted only when it strictly increases Z, the objective sequence is strictly increasing. After discretization to a finite candidate set the process terminates in finitely many steps at a configuration that is Lloyd-stable (Proposition 3.9) and admits no improving single-facility marginal re-placement. Proof Strict-improvement acceptance makes Z monotone strictly increasing; on a finite candidate set the configuration space is finite, so the sequence terminates. At termination, Lloyd stability holds by the location step, and no single-facility re-placement improves by the acceptance rule. □ □ Proposition 3.10 closes the structural argument: the destroy-and-repair relocate, guided by the submodular marginal of Theorem 3.3, escapes the Lloyd fixed point of Proposition 3.9 and terminates at a configuration admitting no further single-facility improvement. Its destroy-and-repair structure, not the heuristic that picks which facility to move, is what carries the gain over multi-start Lloyd—a fact the ablation in §6 confirms. Together with the gradient-as-force location step (Theorem 3.7), the relocate constitutes the second pillar of FBM-LNS (§4). Remark 3.11. Proposition 3.10 analyzes the strict-improvement relocate (accept a move only if it strictly increases Z), under which the objective is monotone increasing and the iterate terminates finitely at a Lloyd-stable, marginal-stable configuration. The implementation of §4 additionally accepts a non-improving relocate with small probability to escape deep fixed points; this non-monotone augmentation is the “mixture” variant isolated by the ablation in Result 4, where it matches the strict variant to within noise. The monotone guarantee therefore applies verbatim to the default strict-improvement path. 4 The FBM-LNS algorithm Our algorithm combines the two pillars established by the structural theory, and we call it FBM-LNS (Force-Based Metaheuristic with Large-Neighborhood Search). The name reflects the two pillars established by the structural theory: the force-based local search, in which each facility’s update direction is the objective gradient ∇XjZ _X_jZ—the demand “pull” (Theorem 3.7), making Lloyd’s location step provably gradient ascent rather than an ad hoc move; and the large-neighborhood search relocate, a destroy-and-repair move guided by the submodular marginal gain (Theorem 3.3) that escapes Lloyd fixed points (Proposition 3.9). Together they yield a multi-start framework whose every component has a formal justification. Algorithm 1 FBM-LNS: Force-Based Metaheuristic with LNS relocate 1:instance (pi,wi,ϕi)i∈I(p_i,w_i, _i)_i∈ I; facility count p; restarts R; LNS iterations L 2:⋆←∅ X ← ; Z⋆←−∞Z ←-∞ 3:for r=1,…,Rr=1,…,R do 4: ← X← weighted k-means++ seed 5: ←Lloyd() X← Lloyd( X) ⊳ location step (5) to convergence 6: (loc,Zloc)←(,Z())( X_loc,Z_loc)←( X,Z( X)) 7: for ℓ=1,…,L =1,…,L do 8: j←j← facility with smallest contribution (or uniform random w.p. ρ) 9: c←c← captured demand of ∖Xj X \X_j\ 10: xnew←argmaxxg(x;c)x_new← _x\,g(x;c) ⊳ greedy repair; Weber solve 11: ′←Lloyd(∖Xj∪xnew) X ← Lloyd( X \X_j\∪\x_new\) 12: if Z(′)>ZlocZ( X )>Z_loc or coin flip ρ then 13: ←′ X← X 14: end if 15: update (loc,Zloc)( X_loc,Z_loc) to best-so-far 16: end for 17: update (⋆,Z⋆)( X ,Z ) with (loc,Zloc)( X_loc,Z_loc) 18:end for 19:return ⋆ X Figure 4 traces one iteration of FBM-LNS through the three stages of seeding, Lloyd polishing, and relocate. (a) Initial: k-means++ seed (b) After Lloyd: gradient-ascent to cell optima removeunder-servedcluster (c) After relocate: destroy-and-repair Figure 4: One iteration of FBM-LNS in three stages. (a) Initial k-means++ seed. (b) Lloyd’s gradient-ascent location step moves each facility to its within-cell Weber optimum; the ghost (dashed gray) squares mark the (a)-time positions, and the dashed orange arrows are the gradient-based moves. (c) The relocate step identifies the least-contributing facility (red circle), removes it, and re-places it at the maximum-marginal point of the under-served cluster (red points, upper-left), whose dashed disk is centered at the re-placement location; the cycle then repeats from (b). Throughout, solid squares are current facilities and dashed squares are past/removed positions. 4.1 How FBM-LNS addresses the structural difficulties Table 2 maps each source of difficulty identified in §2.3 to the structural property (§3) that the algorithm exploits to overcome it. Table 2: Each difficulty of §2.3 and how FBM-LNS addresses it. Difficulty How FBM-LNS addresses it Theory D1. NP-hardness Heuristic search (no exactness claim); quality bounded below by the (1−1/e)(1-1/e) greedy (Theorem 3.3, a valid initial solution). Thm. 3.3 D2. Non-convex objective Within each Voronoi cell the objective reduces to a weighted Weber problem (Theorem 3.4), solvable to its unique local optimum; the LNS relocate escapes the resulting coordinate-wise fixed point (Proposition 3.10). Thm. 3.4, Prop. 3.10 D3. Continuous–combinatorial coupling The Lloyd alternation decouples the two: fix the Voronoi assignment (→ continuous subproblem) or fix locations (→ nearest-facility reassignment). The relocate perturbs the assignment by re-placing one facility. Prop. 3.9 D4. Infinite candidate set Gradient-based location steps move continuously in ℝ2R^2; Theorem 3.2 restricts the search to convpiconv\p_i\; no discretization is needed. Thm. 3.2 D5. Per-point heterogeneity Enters the Weber weights as wi/Riw_i/R_i (Corollary 3.5): distance-sensitive points pull harder, automatically and without ad hoc tuning. Cor. 3.5 4.2 Specialization across decay families The same pseudocode runs on all four decay families—linear, exponential, quadratic, and step—by switching the per-point function ϕi _i. These four shapes underlie the five problem variants V1–V5 of §2.2: V1 and V4 both use the linear shape, at a small clipped scale (gradual coverage) and a large unclipped scale (the Weber limit) respectively. The location step specializes per family as follows: • Quadratic (k-means): the Weber optimum has the closed-form weighted centroid Xj⋆=∑i:π(i)=jwipi/∑wiX_j = _i:π(i)=jw_ip_i/Σ w_i, so the location step is the standard exact Lloyd update (one step, no iteration). • Linear (gradual coverage; p-median with large R): the Weber problem is solved by the Weiszfeld iteration, which converges linearly (Weiszfeld 1937). • Exponential: the gradient is always active (no truncation), so standard normalized gradient ascent applies. • Step (binary): the objective is piecewise constant in facility positions, so the gradient is zero a.e.; only the greedy/relocate steps provide signal. We include this case for completeness but recommend the (1−1/e)(1-1/e) greedy for binary coverage (§6). Per iteration the cost is dominated by the Lloyd step O(np)O(np) and the marginal repair O(n⋅|C|)O(n·|C|) for a candidate set C; both are polynomial. The two primary baselines are the (1−1/e)(1-1/e) continuous greedy (Greedy; Theorem 3.3) and multi-start Cooper alternating location–allocation (Cooper; Cooper 1964); we also report particle swarm optimization (PSO; Kennedy and Eberhart 1995) and decay-agnostic weighted k-means. 5 Exact discrete method and near-optimality The continuous PMC-HDD of §2 optimizes Z()Z( X) over ∈(ℝ2)p X∈(R^2)^p and admits no exact method in general: it is NP-hard (D1) and has no known finite dominating set (D4). Restricting facility locations to a finite candidate set C (a grid, or the demand points themselves) turns the objective into the discrete set function ZD(S)Z_D(S) of Theorem 3.3—a monotone submodular maximum-coverage problem that is solvable to optimality. This discretization serves two purposes below: an exact solver for moderate instances, and—because refining C converges to the continuous optimum—an optimality certificate for the FBM-LNS heuristic. 5.1 The LP-tight MIP formulation With candidates restricted to a finite set C (|C|=m|C|=m), PMC-HDD reduces to the y-assignment max-cover MIP: max ∑i,cwiϕicyic _i,cw_i _ic\,y_ic (6) s.t. yic≤xc y_ic≤ x_c ∀(i,c), ∀(i,c), ∑cyic≤1 _cy_ic≤ 1 ∀i, ∀ i, ∑cxc≤p,xc,yic∈0,1. _cx_c≤ p, x_c,\,y_ic∈\0,1\. The LP relaxation of (6) has worst-case integrality gap 1−1/e1-1/e (the standard maximum-coverage bound). On geographic instances the gap is empirically negligible, a tightness we confirm computationally in Section 6 (Table 9); its structural source is a partition property, not just low submodular curvature (κ): Proposition 5.1 (LP tightness from disjoint coverage). If an optimal solution of (6) has pairwise-disjoint coverage (each demand point is substantially captured by at most one chosen facility, so the facilities induce a partition of the covered points), then the LP relaxation of (6) admits an integral optimum and the integrality gap is 0. Proof Under disjoint coverage the binding constraint for point i is the single facility that covers it, so ∑cwiϕicyic _cw_i _icy_ic is maximized by setting that one yic=1y_ic=1 rather than splitting it fractionally; the LP therefore attains its optimum at an integral (x,y)(x,y). □ □ This is the generic regime, not a knife-edge case: at any Lloyd-stable configuration (Proposition 3.9) the facilities command distinct Voronoi cells, so their coverages are approximately disjoint, and both the exact discrete optimum and FBM-LNS converge to such configurations. Branch-and-bound (HiGHS; Huangfu and Hall 2018) solves (6) on this regime in seconds; we report solve times up to n=500n=500 in Section 6 (Table 9). 5.2 Near-optimality bounds Refining the candidate grid converges the discrete optimum to the continuous one: since the continuous optimum satisfies Z∗≥ZFBM≥Zgrid-MIPZ^*≥ Z_FBM≥ Z_grid -MIP and the grid-MIP rises toward FBM-LNS as the grid refines, all three coincide asymptotically, so the discrete exact value brackets the heuristic. Two propositions make this certificate quantitative by giving computable continuous upper bounds on Z∗Z^*, first at a finite grid and then without one. We confirm empirically that the bracket is tight—FBM-LNS lies within 0.5%0.5\% of the exact optimum, and within the [0.99,1.03][0.99,1.03] ratio on small exhaustively-solved instances—in Section 6 (Result 5). Proposition 5.2 (Continuous upper bound via grid density). Let C be a grid of spacing h covering convpiconv\p_i\, so every hull point is within h/2h/ 2 of a grid node, and let Lϕ=∑iwiLip(ϕi)L_φ= _iw_i\,Lip( _i), where Lip(ϕi)Lip( _i) is a Lipschitz constant of ϕi _i (1/Ri1/R_i for linear and exponential decay). Then Z∗≤Zgrid-MIP(C)+h2Lϕ≤LP(C)+h2Lϕ,Z^*\;≤\;Z_grid -MIP(C)+ h 2\,L_φ\;≤\;LP(C)+ h 2\,L_φ, so the grid LP is a verifiable upper bound on the continuous optimum, with slack h2Lϕ h 2L_φ that vanishes as h→0h→ 0. Proof sketch. Round each XjX_j to its nearest grid node X~j X_j (with ‖Xj−X~j‖≤h/2\|X_j- X_j\|≤ h/ 2); Lipschitz continuity of each ϕi _i then costs at most (h/2)∑iwiLip(ϕi)(h/ 2) _iw_i\,Lip( _i), and taking the supremum over X yields the claim. The full proof is in the Electronic Companion. □ □ The bound is conservative (worst-case Lipschitz) but valid at any finite grid, turning the grid LP from an informal certificate into an upper bound on the true continuous optimum; combined with Z∗≥Zgrid-MIP≥ZFBMZ^*≥ Z_grid -MIP≥ Z_FBM, it brackets the heuristic within the slack h2Lϕ h 2L_φ of Z∗Z^*. At practical grid resolutions this slack is of the same order as Z (the worst case charges every point at its full Lipschitz rate), so the proposition is a theoretical guarantee rather than a tight certificate; the tight practical certificate is the fine-grid convergence confirmed in Section 6. Proposition 5.3 (Grid-free continuous upper bound for concave decays). For a decay that is concave in the distance (unclipped linear ϕi(r)=1−r/Ri _i(r)=1-r/R_i and quadratic ϕi(r)=1−(r/Ri)2 _i(r)=1-(r/R_i)^2), Z∗≤p⋅maxx∈convpi∑iwiϕi(‖pi−x‖)Z^*≤ p· _x \p_i\ _iw_i\, _i(\|p_i-x\|), and the inner single-facility maximum is attained at a demand point, so the bound is computable in O(n2)O(n^2) without a grid. Proof sketch. Use maxjϕi(‖pi−Xj‖)≤∑jϕi(‖pi−Xj‖) _j _i(\|p_i-X_j\|)≤ _j _i(\|p_i-X_j\|), separate the right-hand side over facilities, and observe that for concave non-increasing ϕi _i the single-facility term is concave in x (Lemma EC.2 in the Electronic Companion), hence maximized over the hull at an extreme point (a demand point). The full proof is in the Electronic Companion. □ □ This grid-free bound applies to the trap-free unclipped decays identified in Theorem 3.8(i); for clipped-linear and exponential decays the single-facility term is not concave in x and the maximum need not occur at a demand point, so the discretization bound of Proposition 5.2 remains the operative certificate. The bound is loose even where it applies, since maxj≤∑j _j≤ _j lets all p facilities cover every point. Both continuous bounds are therefore valid but not tight; the operative practical certificate remains the discretized LP/MIP, whose empirical tightness we report in Section 6. A genuinely tight continuous relaxation would have to handle the pointwise maximum directly, which is a convex maximization (or a reverse second-order cone) and hence non-convex; it is left to future work. 6 Computational experiments 6.1 Data sources and experimental setup To ensure authority and representativeness, we draw on established benchmark families for the spatial distribution of demand points in a 1000×10001000× 1000 square. We label each benchmark a test set, abbreviated TS, and index them TS1–TS7 for cross-reference: TS1. Concentration-point (Bansal–Kianfar): the procedure of Bansal and Kianfar (2017), used in Bansal and Shojaee (2020)—three random centers of radius 270270, with each point anchored to a center (probability 0.310.31 each) or free uniform (0.070.07). This mimics real demand clustering around population centers. TS2. Uniform: i.i.d. points on the square—the standard worst-case baseline for location algorithms. TS3. Gaussian mixture: k isotropic components—the canonical stochastic model for k-means and clustering analysis (Ostrovsky et al. 2013). TS4. Disk-clustered: uniform within k random disks—a structured but non-Gaussian alternative. TS5. Franti s1--s4: a real clustering benchmark of 5,000 points in 15 Gaussian clusters—the standard test set for k-means (Fränti and Sieranoja 2018). TS6. TSPLIB EUC_2D: six instances (eil51 through u1817) giving real geographic city coordinates, used for the Weber/p-median experiments (Reinelt 1991). TS7. Region demand: disk and rectangle shapes sampled into density-weighted points, as described under variant V5. Weights are wi∼U[1,10]w_i U[1,10]. The per-point decay scale RiR_i is drawn from a lognormal with mean Rmean=120R_mean=120 and a controlled coefficient of variation RcvR_cv: Rcv=0R_cv=0 gives the homogeneous control and larger RcvR_cv increases heterogeneity (the defining feature of our setting). Every instance is a pure function of (parameters, seed); the full library and raw results are released for reproducibility. Table 3 summarizes the seven test sets; Figure 5 shows one realization of each synthetic generator (TS1–TS4), while TS5–TS7 are real geographic or region datasets and appear only in the table. All comparisons are paired (Wilcoxon signed-rank). Table 3: Benchmark datasets used in the computational study. Code Source Geography / type Variant(s) Prior use TS1 Bansal–Kianfar Concentration-point, heterogeneous RiR_i V1, V2, V6 Bansal and Kianfar 2017, Bansal and Shojaee 2020 TS2 Uniform on square Diffuse, homogeneous or heterogeneous V1, V2, V6 standard baseline TS3 Gaussian mixture Clustered, k components V1–V4, V6 Ostrovsky et al. 2013, Arthur and Vassilvitskii 2007 TS4 Disk-clustered Structured non-Gaussian clusters V1–V4, V6 synthetic (this work) TS5 Franti s1--s4 Real clustering benchmark (5,000 pts) V3 Fränti and Sieranoja 2018 TS6 TSPLIB EUC_2D Real geographic coordinates V4 Reinelt 1991 TS7 Disk / rectangle shapes Region demand, sampled to points V5 Drezner et al. 2004 Note. Codes are used throughout §6 to identify the test set for each experiment. Figure 5: Sample demand geographies used in the experiments (n=400n=400, Rcv=0.5R_cv=0.5): (a) TS1 concentration-point, (b) TS2 uniform, (c) TS3 Gaussian mixture, (d) TS4 disk-clustered. Points are colored by decay scale RiR_i (red ≡ distance-sensitive, blue ≡ distance-tolerant) and sized by weight. The four generators span concentrated, diffuse, and clustered demand, so that no conclusion rests on a single spatial pattern. 6.2 Result 1: FBM-LNS significantly outperforms all baselines (30/3030/30 wins, p<10−9p<10^-9) Table 4 reports mean captured demand Z over K=8K=8 instances per configuration under heterogeneous decay (Rcv=0.5R_cv=0.5), for two scales (n,p)∈(300,5),(1000,10)(n,p)∈\(300,5),(1000,10)\. FBM-LNS yields the highest mean Z in every configuration, with a margin over the strong (1−1/e)(1-1/e) greedy baseline that is stable at 1.51.5–3.0%3.0\% and a positive margin over Cooper ALT in every cell. The large-sample paired comparison (K=30K=30, Table 5) is significant at p<10−9p<10^-9 with 30/3030/30 per-instance wins. Table 4: (Result 1) Mean captured demand Z on the main matrix (Rcv=0.5R_cv=0.5, K=8K=8 instances per cell). FBM-LNS attains the highest mean Z in every configuration against every baseline; large-sample significance (30/3030/30 wins, p<10−9p<10^-9) is confirmed in Table 5. Greedy is the (1−1/e)(1-1/e) continuous greedy; Cooper is multi-start alternating location–allocation. “Code” tags each row with its test set (TS, Table 3) and decay variant (V, §2.2). Code Configuration k-means Cooper PSO Greedy FBM-LNS vs. Greedy TS1·V1 conc. / linear, n=300,p=5n=300,p=5 512 572 550 583 597 +2.4%+2.4\% TS1·V2 conc. / exp, n=300,p=5n=300,p=5 773 808 797 805 821 +2.0%+2.0\% TS3·V1 gmix / linear, n=300,p=5n=300,p=5 694 719 630 705 726 +3.0%+3.0\% TS3·V2 gmix / exp, n=300,p=5n=300,p=5 925 936 887 925 941 +1.7%+1.7\% TS1·V1 conc. / linear, n=1000,p=10n=1000,p=10 2176 2416 2171 2489 2552 +2.5%+2.5\% TS1·V2 conc. / exp, n=1000,p=10n=1000,p=10 2977 3111 2905 3146 3194 +1.5%+1.5\% TS3·V1 gmix / linear, n=1000,p=10n=1000,p=10 2887 2978 2494 2916 2996 +2.7%+2.7\% TS3·V2 gmix / exp, n=1000,p=10n=1000,p=10 3498 3554 3165 3506 3563 +1.6%+1.6\% Note. Each row averages K=8K=8 independent instances; FBM-LNS attains the highest mean in every cell. The rigorous large-sample significance test (K=30K=30) is reported in Table 5. Table 5 confirms that this superiority is not an artifact of small samples: with K=30K=30 independent instances per cell, FBM-LNS wins 30/3030/30 paired comparisons in every cell and the paired Wilcoxon p-value is below 10−910^-9 throughout. The per-instance, within-configuration difference is what is significant; the unpaired 95%95\% confidence intervals of the two means do overlap, because the between-instance variability in Z (driven by the random geography) is much larger than the within-instance algorithm gap. The advantage is stable across scales (from n=300n=300 to n=1000n=1000) and across both linear and exponential decay. Table 5: (Result 1, large-sample confirmation) Mean captured demand ±95%± 95\% confidence interval (CI) over K=30K=30 independent instances, with paired Wilcoxon p-values. FBM-LNS wins 30/3030/30 paired comparisons in each cell at p<10−9p<10^-9; the unpaired CIs overlap because between-instance variability dominates the within-instance gap. Code Configuration FBM-LNS Greedy (1−1/e)(1-1/e) FBM vs. Greedy p-value TS1·V1 conc./linear, n=300,p=5n=300,p=5 548.5±23.2548.5± 23.2 536.2±22.8536.2± 22.8 +2.3%+2.3\% <10−9<10^-9 TS1·V1 conc./linear, n=1000,p=10n=1000,p=10 2373.2±77.42373.2± 77.4 2291.9±76.72291.9± 76.7 +3.5%+3.5\% <10−9<10^-9 TS3·V2 gmix/exp, n=300,p=5n=300,p=5 904.8±18.4904.8± 18.4 890.4±18.1890.4± 18.1 +1.6%+1.6\% <10−9<10^-9 Note. Paired Wilcoxon signed-rank test, alternative “FBM-LNS greater than Greedy”; 30/3030/30 per-instance wins in each configuration. The superiority has a structural explanation, not just an empirical one. The (1−1/e)(1-1/e) greedy is a constructive heuristic: it places facilities one at a time at the maximum-marginal point but never re-optimizes earlier placements, so its first few choices—made when coverage was still low—are frozen as later demand is revealed. Cooper ALT is multi-start Lloyd only, so it converges to a coordinate-wise fixed point (Proposition 3.9) and cannot escape it. FBM-LNS closes both gaps: its Lloyd step jointly re-optimizes every facility within its Voronoi cell (Theorem 3.4), undoing the greedy’s frozen early choices, and its destroy-and-repair relocate removes the least-contributing facility and re-places it where the residual marginal g(x)g(x) is maximal (Proposition 3.10), escaping the Lloyd fixed point. The two mechanisms are complementary, which is why the ablation in Result 4 attributes the gain to the relocate rather than to the selection rule. The same mechanism also explains the stable advantage over Cooper ALT across scales: Cooper’s multi-start Lloyd converges to a coordinate-wise fixed point, and the relocate’s destroy-and-repair step escapes fixed points that more restarts alone cannot. 6.3 Result 2: One framework is competitive with bespoke SOTA across the problem family A defining feature of our framework is that the same implementation, switching only the per-point function, runs on the other two canonical continuous planar location problems and on shape-based demand. (i) k-means: minimize ∑iminj‖pi−Xj‖2 _i _j\|p_i-X_j\|^2, recovered by quadratic decay with a large homogeneous scale, where the within-cell Weber optimum has the closed form of the weighted centroid (Corollary 3.6), so our location step specializes to the exact Lloyd update. (i) Weber/p-median: minimize ∑iwiminj‖pi−Xj‖ _iw_i _j\|p_i-X_j\|, recovered by linear decay with a large scale. (i) Shape demand: when each demand unit is a region DiD_i rather than a point, the captured demand ∫Diϕi(minj‖p−Xj‖)μi(p)p _D_i _i( _j\|p-X_j\|) _i(p)\,dp is reduced to the weighted-point problem by sampling each region into m density-weighted points (Monte-Carlo error O(1/m)O(1/ m)); allocation is still the Voronoi partition of the plane, and the gradient is the integral of the point gradient over the region. Figure 6 illustrates the sampling-based reduction. F1F2D1D_1 (disk)D2D_2 (rect.)Voronoi edge Figure 6: Shape (region) demand. Each demand unit is a region (disk or rectangle) rather than a point; it is reduced to the weighted-point problem by sampling its interior into density-weighted points. The nearest-facility allocation remains the Voronoi partition; the coverage objective is the integral of the per-point decay over the region. Table 6 evaluates the framework against the problem-specific state of the art: k-means++ on k-means (synthetic Gaussian blobs and the Franti s1--s4 benchmark of 5,000 points with 15 clusters), multi-start Cooper and variable neighborhood search (Hansen and Mladenović 1997) on p-median, and multi-start Cooper on shape-demand instances. FBM-LNS matches k-means++ on both synthetic and real k-means benchmarks (tying it on s1--s4), and on p-median it is competitive with a specialized VNS: it beats VNS on the larger instances (d657 by 1.2%1.2\%, pr1002 by 1.8%1.8\%, the concentration-point instance by 0.3%0.3\%) but loses to it on the small eil101 by 8.2%8.2\%, where the demand-point swap local search of VNS is strongest. The point is not a new record on these thoroughly studied problems but that one method, designed from the shared structure, is competitive with problem-specific standard methods across the family while extending to heterogeneous-decay and shape demand where those methods do not directly apply. Table 6: (Result 2) The same framework across the continuous-location family. “Code” tags each row with its test set (TS) and variant (V); Δ is FBM-LNS relative to the best baseline in the row (see the table note for the direction-of-better convention). One code, switching only the per-point function, is competitive with or beats problem-specific standard methods (k-means++, Cooper) on k-means, p-median, and shape-demand instances. Code Problem Configuration Baseline Baseline Z FBM-LNS Δ TS3·V3 k-means blobs, n=400,p=5n=400,p=5 km++ 1 138 041 1 138 039 0.00%0.00\% TS3·V3 k-means blobs, n=800,p=10n=800,p=10 km++ 1 982 174 1 966 505 −0.79%-0.79\% TS5·V3 k-means real s1--s4, n=5000,k=15n=5000,k=15 km++ FBM ties km++ (identical objective) 0.00%0.00\% TS1·V4 p-median conc., n=500,p=10,K=8n=500,p=10,K=8 VNS 191 380 190 733 −0.34%-0.34\% TS6·V4 p-median TSPLIB eil101, n=101,p=10n=101,p=10 VNS 767 830 +8.2%+8.2\% TS6·V4 p-median TSPLIB d657, n=657,p=30n=657,p=30 VNS 103 231 101 970 −1.22%-1.22\% TS6·V4 p-median TSPLIB pr1002, n=1002,p=30n=1002,p=30 VNS 690 249 677 823 −1.80%-1.80\% TS7·V5 shape 25 disks, p=6p=6 Cooper 202.1 211.7 +4.74%+4.74\% TS7·V5 shape 25 rectangles, p=6p=6 Cooper 200.1 211.6 +5.74%+5.74\% Note. “km++” denotes k-means++ seeding followed by Lloyd; “VNS” is variable neighborhood search (Hansen and Mladenović 1997) for p-median (facilities restricted to demand points, fast-interchange local search plus k-swap shaking). The “Baseline” column names the stronger of Cooper and VNS for the row, and “Baseline Z” reports its objective value. Rows average 4–30 instances. For k-means, FBM-LNS uses the closed-form centroid location step (Corollary 3.6); for shapes, each region is sampled into 40 density-weighted points. Boldface marks the better of the two values. Direction: for k-means (sum of squared errors) and p-median (sum of nearest-facility distance) the objective is a cost, so lower is better and a negative Δ favors FBM-LNS; for shape demand the objective is captured demand, so higher is better and a positive Δ favors FBM-LNS. Δ is always (FBM-LNS −- best baseline)//best baseline. The competitiveness is structural, not coincidental: because k-means, the Weber problem, and heterogeneous gradual coverage are all instances of the single objective (1) (Theorem 3.1), the location step of FBM-LNS reproduces each bespoke method’s core operation—the closed-form centroid for k-means (Corollary 3.6), the Weiszfeld iteration for p-median (§4.2)—and adds the same destroy-and-repair relocate on top. The single code therefore inherits the strength of each specialist while extending to heterogeneous-decay and shape demand where they do not directly apply. The sampling reduction used for shape demand (V5) does not introduce solution instability. Holding the demand regions fixed and re-sampling each into m=40m=40 points with 1212 independent seeds, the captured demand varies by under 1%1\% (coefficient of variation 0.86%0.86\% for disks, 0.83%0.83\% for rectangles), and the optimal facility positions move by under 1%1\% of the map width on average (0.67%0.67\% for disks, 0.89%0.89\% for rectangles; rectangle instances occasionally show a larger position swap, but it leaves the captured demand essentially unchanged). The placement returned for region demand is therefore stable and seed-independent in objective value. 6.4 Result 3: The advantage grows with heterogeneity and persists at scale Table 7 sweeps Rcv∈0,0.25,0.5,1.0R_cv∈\0,0.25,0.5,1.0\ at n=300n=300, p=5p=5, conc./linear. FBM-LNS is best at every heterogeneity level, and its advantage over Cooper ALT widens as heterogeneity grows (from +2.4%+2.4\% when Rcv=0.25R_cv=0.25 to +6.8%+6.8\% when Rcv=1.0R_cv=1.0), confirming that the submodular relocation pays off most when sensitivity varies most. Table 7: (Result 3) Heterogeneity sweep on configuration TS1·V1 (conc./linear, n=300n=300, p=5p=5, K=8K=8). FBM-LNS is best at every level; its advantage over Cooper ALT grows from +3.2%+3.2\% (Rcv=0R_cv=0) to +7.0%+7.0\% (Rcv=1R_cv=1), confirming that the relocate step pays off precisely when sensitivity varies most. RcvR_cv Greedy Cooper FBM-LNS FBM vs. Cooper 0.000.00 (homog.) 610 618 638 +3.2%+3.2\% 0.250.25 604 612 627 +2.4%+2.4\% 0.500.50 583 572 597 +4.4%+4.4\% 1.001.00 (strong) 522 498 532 +6.8%+6.8\% Note. All differences FBM-LNS vs. each baseline significant at p<0.01p<0.01 (Wilcoxon, 8/88/8 wins). The widening gap follows directly from the Weber reduction (Theorem 3.4). As RcvR_cv grows, the Weber weights αi=wi/Ri _i=w_i/R_i become more dispersed, so the within-cell objective is increasingly dominated by a few distance-sensitive demand clusters. Cooper ALT, which does not relocate facilities once Lloyd has converged, leaves facilities in cells whose sensitive core has shifted; FBM-LNS’s relocate repeatedly identifies the cell whose contribution has fallen behind and re-places its facility where the residual marginal g(x)=∑iwi[ϕi(‖pi−x‖)−ci/wi]+g(x)= _iw_i[ _i(\|p_i-x\|)-c_i/w_i]^+ is maximal. The hotter the heterogeneity, the larger these residuals and the more the relocate recovers, which is exactly the monotone trend in Table 7. This confirms Corollary 3.5 empirically: heterogeneity is not a nuisance the algorithm tolerates but the structural condition under which its advantage is largest. 6.5 Result 4: The LNS relocate step drives the improvement (+6.5%+6.5\% over multi-start Lloyd) Table 8 decomposes the algorithm. Multi-start Lloyd on its own is outperformed by the (1−1/e)(1-1/e) greedy; adding the large-neighborhood relocate lifts performance by about 6.5%6.5\%, and the specific relocation selector—least-contributing greedy, uniform random, or their mixture—makes essentially no difference. The improvement comes from the destroy-and-repair structure with a maximum-marginal re-placement (the submodular signal), not from the heuristic that picks which facility to move. Table 8: (Result 4) Component ablation on configuration TS1·V1 (conc./linear, n=300n=300, p=5p=5, Rcv=0.5R_cv=0.5, K=10K=10). Multi-start Lloyd alone underperforms even the greedy; adding the relocate step lifts performance by +6.5%+6.5\%, while the selection strategy (greedy vs. random) is immaterial—the destroy-and-repair structure is the source of the gain. Configuration: TS1·V1 — conc./linear, n=300n=300, p=5p=5, Rcv=0.5R_cv=0.5, K=10K=10 Variant mean Z vs. Lloyd-only Cooper (multi-start Lloyd) 549.9 +0.8%+0.8\% Greedy (1−1/e)(1-1/e) 567.5 +4.0%+4.0\% FBM-LNS: Lloyd only 545.5 — FBM-LNS: ++ greedy relocate 580.3 +6.4%+6.4\% FBM-LNS: ++ random relocate 580.8 +6.5%+6.5\% FBM-LNS: full (mixture) 581.0 +6.5%+6.5\% Note. All “++ relocate” variants are significantly better than Lloyd-only and than Greedy at p<0.01p<0.01 (Wilcoxon). Why does the relocate dominate while the selection rule is immaterial? The relocate’s gain comes from its repair operator, which solves the single-facility maximum-marginal subproblem g(x)g(x)—a weighted Weber problem on the residual active set (Theorem 3.4)—and is therefore the continuous analogue of the submodular greedy step (Theorem 3.3). The destroy operator only decides which facility to remove, and on a Lloyd-stable configuration any reasonable choice (least-contributing or random) exposes a comparable residual, so the marginal gain comes from re-placing at the maximum of g, not from which facility was freed. This is why Table 8 shows near-identical performance across selection rules but a +6.5%+6.5\% jump once the relocate is active, and it is also why the learned selector of §8 recovers only a marginal extra gain: the signal is already captured by the destroy-and-repair structure itself. Figure 7 traces the best-so-far objective over relocate iterations. Figure 7: Best-so-far captured demand of FBM-LNS over relocate iterations on one representative instance (configuration TS1·V1: conc./linear, n=300n=300, p=5p=5), with the greedy and multi-start Cooper values as horizontal references. The relocate steps lift the solution above both baselines within the first few iterations. 6.6 Result 5: The LP relaxation is empirically tight (≈0%≈ 0\% gap) and FBM-LNS is within 0.5%0.5\% of the discrete optimum by fine-grid convergence Table 9 reports the LP integrality gap, MIP solve time, and the fine-grid convergence gap. Table 9: (Result 5) LP tightness and MIP scalability. The LP relaxation is empirically tight (≈0%≈ 0\% gap) across all configurations. Each row is one representative bansal/linear instance (Rcv=0.5R_cv=0.5, seed 500) discretized on a grid-10 candidate set built from the top 40 demand points plus a uniform grid (see src/bounds.py); the 0%0\% gap is robust across seeds and decay families. n p grid-MIP grid-LP relax. gap MIP time (s) 100 3 179.8 179.8 0.00%0.00\% 0.1 200 5 369.9 369.9 0.00%0.00\% 0.3 300 5 542.1 542.1 0.00%0.00\% 0.5 300 6 600.0 600.0 0.00%0.00\% 0.6 500 5 870.8 870.8 0.00%0.00\% 1.3 The LP gap is ≈0%≈ 0\% even under adversarial high-overlap configurations (clustered points, small R). On small exhaustively-solved instances (n≤100n≤ 100, p≤3p≤ 3) FBM-LNS attains ZFBM/Zexact∈[0.99,1.04]Z_FBM/Z_exact∈[0.99,1.04] and can even exceed the grid-constrained optimum, because it places facilities continuously in the hull (Theorem 3.2) rather than on grid nodes. Fine-grid convergence: at grid resolution 4040, the gap between FBM-LNS (continuous) and the exact discrete optimum is <0.5%<0.5\%. The tightness has a structural source. The integrality gap of the assignment formulation (6) is governed by the submodular curvature κ of Theorem 3.3: once p facilities are spread over the hull their service areas overlap little, κ is small, and the LP relaxation, the greedy, and the integer optimum coincide—so the gap collapses to ≈0%≈ 0\% even though the worst-case bound remains 1−1/e1-1/e. The <0.5%<0.5\% gap of FBM-LNS follows because its Lloyd step provably reaches the within-cell Weber optimum (Theorem 3.4, unique by Corollary 3.6), so the only remaining source of suboptimality is the allocation of points to facilities (the Lloyd fixed points of Proposition 3.9), which the relocate removes (Proposition 3.10). The exact discrete optimum and FBM-LNS therefore agree once the grid is fine enough to resolve the optimal partition. 6.7 Result 6: On binary (step) decay the greedy matches or exceeds FBM-LNS The smooth-decay advantage of Result 1 reverses on binary step decay (V6), where the within-cell objective is piecewise constant in facility position and the gradient channel is uninformative. On the step-decay matrix the (1−1/e)(1-1/e) greedy and FBM-LNS are within 1%1\% at n=300,p=5n=300,p=5 (FBM marginally ahead on three of four geographies), but at n=600,p=8n=600,p=8 the greedy leads by 22–5%5\% on the clustered geographies and FBM wins only 0–11 of 66 instances there. The gradient-based location step thus adds no value when ∇Xjϕi≡0 _X_j _i≡ 0, and the submodular-marginal greedy is the appropriate tool for binary coverage; we return to the first-principles reason in §8.1. 7 Real-world case study 7.1 Data Both datasets are drawn from the 2026 6th Meituan Business Analytics Elite Competition. The first is a set of 592,667592,667 urban-delivery waybills from the Nanshan district of Shenzhen, each carrying pickup/dropoff coordinates and a business-line label (food delivery, flash sale, pharmacy, group meal). We treat each dropoff as a demand unit and aggregate the 592,667592,667 orders into 347347 demand cells weighted by order volume wiw_i, keeping the planning problem tractable while preserving the spatial demand distribution. The per-cell decay scale RiR_i follows a density-based modeling proxy rather than a direct estimate from observed distance-sensitivity: we set RiR_i small in dense cells (where customers are impatient and substitutes abound) and large in sparse residential cells (where a courier is the only nearby option), encoding the paper’s central dense/sparse distinction as a per-cell parameter. We treat this as a transparent modeling assumption and do not independently validate it against the data; direct per-zone calibration of RiR_i (for example from observed delivery distances) is a natural extension. The resulting Ri\R_i\ span [79,1953][79,1953] metres with coefficient of variation 1.331.33, strong heterogeneity, exactly the regime in which Result 3 predicts the largest advantage. The direction of this density-to-range assumption is empirically supported by the same data: across 347347 spatial cells, cell density correlates at −0.40-0.40 with the median pickup-to-dropoff trip distance (densest-quartile cells about 1.21.2 km versus sparsest about 2.12.1 km), and the overall median trip of 1.471.47 km matches the retail-calibrated decay scale R≈1.4R≈ 1.4 km, so the assumed dense-then-tolerant heterogeneity is consistent with observed delivery patterns. The second dataset is a retail-chain panel of per-store sales stratified by distance tier, used below to calibrate the decay shape (exponential, R≈1.4R≈ 1.4 km; root-mean-square error RMSE 0.120.12 versus 0.170.17 for linear). Figure 8 maps the spatial distribution of the dropoff locations: demand concentrates in a few dense hotspots (the urban core and commercial strips) against a dispersed residential periphery—precisely the dense/sparse contrast that the density-derived Ri\R_i\ turns into heterogeneity. Figure 8: Spatial distribution of the 592,667592,667 delivery dropoffs (Shenzhen Nanshan), shaded by order density on a log scale. Demand clusters in dense hotspots (the urban core and commercial strips) with a sparse residential periphery—the dense/sparse contrast from which the per-cell decay scales RiR_i are derived—small RiR_i in the dense hotspots and large RiR_i in the sparse periphery (Data source: Meituan). 7.2 Decay calibration from retail data Table 10 reports the empirical per-capita capture by distance tier (from the retail dataset) and the fitted decay parameters; Figure 9 visualizes the two fits. Table 10: Empirical decay calibrated from per-capita sales by distance tier in the retail-chain panel (Data source: Meituan). “Pred.” columns give the fitted models’ predicted normalized per-capita capture; R is the decay scale in metres. Exponential fits better (RMSE 0.120.12 vs. 0.170.17). Distance (m) Normalized capture Linear pred. (R=3384R=3384 m) Exp. pred. (R=1366R=1366 m) 125 1.00 0.96 0.91 375 0.98 0.89 0.76 750 0.58 0.78 0.58 1500 0.26 0.56 0.33 2500 0.26 0.26 0.16 The exponential decay (R≈1.4R≈ 1.4 km) fits the real per-capita capture curve better than the linear, justifying the decay form used in the model. The exponential fit is not only tighter but physically motivated: it matches the path-loss law that also underlies the Huff gravity model (§1.1), so the same decay family is supported by both the delivery-density argument (per-cell RiR_i) and the sales-distance argument (retail calibration). The near-plateau of normalized capture at 750750–25002500 m (both tiers ≈0.26≈ 0.26) suggests a core catchment within roughly 11 km and a thin long-distance tail, consistent with a quick-commerce setting in which most demand is impulse- and proximity-driven. 05005001,0001,0001,5001,5002,0002,0002,5002,50000.50.511distance from store (m)normalized per-capita captureempirical (retail panel)linear (R=3384R=3384 m, RMSE 0.170.17)exponential (R=1366R=1366 m, RMSE 0.120.12) Figure 9: Decay-shape calibration from the retail panel. The empirical per-capita capture (points) drops sharply within ≈1≈ 1 km and then plateaus; the exponential fit (dashed red, RMSE 0.120.12) fits the overall decay pattern better than the linear fit (solid blue, RMSE 0.170.17), justifying the exponential decay used in the model (Data source: Meituan). 7.3 The cost of ignoring heterogeneity Table 11 compares (a) optimizing with the true heterogeneous RiR_i vs. (b) optimizing with the homogeneous mean R¯ R then evaluating on the true objective. Table 11: Cost of ignoring heterogeneity on 592,667592,667 real orders (Data source: Meituan). p hetero Z ignore-hetero Z demand lost facility shift 5 183,510 165,673 9.7% 29% of map 8 240,287 220,809 8.1% 37% 12 295,875 289,244 2.2% 32% Under this density-derived heterogeneity model, optimizing with the homogeneous mean R¯ R instead of the heterogeneous Ri\R_i\ loses up to 9.7%9.7\% of captured demand and physically relocates facilities by up to 37%37\% of the map. The absolute level depends on the (unvalidated) RiR_i assignment, but the directional penalty of ignoring heterogeneity is robust to it: any model in which sensitive and tolerant cells coexist pays for averaging them away. Two patterns in Table 11 deserve emphasis. First, the demand loss is largest when facilities are scarce (9.7%9.7\% at p=5p=5) and shrinks as p grows (2.2%2.2\% at p=12p=12): with many facilities every neighborhood is close to one regardless of planning, so the homogeneity assumption incurs a smaller penalty, whereas with few facilities, misplacing them in tolerant zones leaves the sensitive hotspots under-served. Second, the facility shift (29–37% of the map width) is large at every p: the homogeneous solution systematically over-serves sparse areas and under-serves dense ones, so even when the two plans capture similar total demand their physical locations differ substantially. For a planner, the demand figure quantifies the revenue at stake and the shift figure the operational disruption of using a misspecified model. 7.4 Competitive (Huff) extension Under the Huff model (market-share splitting vs. an incumbent with attractiveness C), the captured share (Table 12) ranges from 37.6%37.6\% (p=12p=12, weak incumbent C=0.5C=0.5) to 7.8%7.8\% (strong C=8C=8), with diminishing returns in p consistent with submodularity. Table 12: Captured market share (% of total demand) under the Huff competitive model, vs. incumbent strength C (Data source: Meituan). C p=5p=5 p=8p=8 p=12p=12 0.5 26.0% 31.9% 37.6% 1.0 17.7% 24.2% 27.0% 2.0 12.8% 14.6% 18.7% 4.0 5.9% 10.0% 12.4% 8.0 3.4% 6.4% 7.8% The pattern tracks the theory: against a weak incumbent (C=0.5C=0.5) additional facilities steadily grow share (26.0%→31.9%→37.6%26.0\%→ 31.9\%→ 37.6\%), whereas against a strong incumbent (C=8C=8) the entrant is capped near 8%8\% regardless of p, and the marginal gain of each added facility diminishes, confirming submodularity (Theorem 3.3) empirically. For market entry, this separates contestable markets (low C), where investment in facilities pays off, from entrenched ones (high C), where it does not. 8 Discussion 8.1 Where does FBM-LNS fail? Result 6 (§6.7) shows that FBM-LNS’s advantage is specific to smooth (linear, exponential) decay, where the gradient/Weber structure of §3 is informative; under binary step decay the (1−1/e)(1-1/e) greedy—which depends only on the submodular marginal, not on gradients—matches or exceeds it. We therefore recommend FBM-LNS for gradual coverage and the greedy for binary coverage. The boundary has a first-principles explanation consistent with the tractability classification (Theorem 3.8). The step function [r≤Ri]1[r≤ R_i] is piecewise constant, so ∇Xjϕi _X_j _i is zero almost everywhere, and the gradient-as-force identity (Theorem 3.7) delivers no direction in which to move a facility; the relocate, by contrast, depends only on the submodular marginal g(x)=∑iwi[ϕi(‖pi−x‖)−ci/wi]+g(x)= _iw_i[ _i(\|p_i-x\|)-c_i/w_i]^+, which is well-defined for any non-increasing ϕi _i including the step. The dichotomy is fundamental: the continuous signal (where to move a facility) requires smoothness, while the combinatorial signal (where to add one) does not. FBM-LNS exploits both; at the step-decay boundary the gradient channel closes but the relocate still provides a signal, so the method falls back to its greedy floor rather than producing no output. A practical corollary is that for mixed objectives with a smooth component (e.g., a “soft” step ϕi(r)=12[1+cos(πr/Ri)][r≤Ri] _i(r)= 12[1+ (π r/R_i)]1[r≤ R_i]) FBM-LNS remains effective, because the smooth part restores the gradient signal. 8.2 Runtime and scalability On the scalability benchmark (Figure 10), wall-clock times per instance are under 11–33 s for the baselines and about 11–1010 s for FBM-LNS across n=200n=200–20002000 (with the lighter LNS budget used there); the main-matrix runs use a heavier budget (≈15≈ 15–2020 s per instance at n=300n=300–10001000). The exact MIP adds 0.40.4–1.41.4 s for n≤500n≤ 500. The additional runtime yields quality improvements that grow with scale and heterogeneity (Tables 4 and 7); where runtime is critical, the LNS iteration budget can be reduced at a modest quality cost. Figure 10 traces quality and runtime up to n=2000n=2000, p=40p=40: FBM-LNS scales near-linearly in time and the quality advantage over the (1−1/e)(1-1/e) greedy persists at the largest sizes. The cost gap has a precise source. Each relocate iteration solves a single-facility maximum-marginal subproblem on the residual (a weighted Weber problem, O(n⋅|C|)O(n·|C|) for the candidate grid C) followed by a Lloyd polish (O(np)O(np) per Lloyd step); over L≈30L≈ 30 iterations and R≈6R≈ 6 restarts this dominates the baselines, which perform only a single Lloyd sweep. The near-linear scaling in n follows because both the Lloyd step and the marginal solve are linear in n and the number of effective iterations to convergence does not grow with problem size—the partition structure stabilizes once each facility commands a cell, regardless of how many points it contains. This is also why reducing the LNS budget trades quality for time gracefully: each iteration is an independent destroy-and-repair trial, so halving L halves runtime at the cost of fewer escapes from Lloyd fixed points. Figure 10: Scalability on configuration TS1·V1 (conc./linear, Rcv=0.5R_cv=0.5). Left: captured demand Z (FBM-LNS above greedy throughout). Right: runtime (near-linear in n on the log axis); at n=2000n=2000, p=40p=40 FBM-LNS finishes in about 1010 s. 8.3 Can FBM-LNS be improved further? The relocate step selects which facility to move. The default rule (least-contributing) is already near-optimal, which the ablation above shows is not the source of the gain. We nonetheless tested whether a learned selector can improve on it: a linear policy scorej=w⊤featjscore_j=w feat_j over five per-facility features (contribution, cell size, mean cell coverage, isolation, bias), trained by the cross-entropy method (CEM; De Boer et al. 2005) to maximize captured demand on training instances, and evaluated on held-out test instances (conc./linear, n=300n=300, p=5p=5, K=10K=10). The learned policy attains Z=582.8Z=582.8 versus 581.5581.5 for the heuristic (and 567.8567.8 for the greedy), a small but statistically significant improvement (+0.22%+0.22\%, p=0.001p=0.001, 10/1010/10 wins). The learned weights place large negative mass on the contribution feature, i.e. the optimizer rediscovers the relocate-the-low-contributor rule and adds a marginal correction from the other features. This confirms that the hand-designed selector already captures most of the selection signal, while a learned policy extracts a small additional gain. The small headroom is itself a consequence of the theory, not a failure of the learning procedure. Result 4 showed that the destroy-and-repair structure accounts for the entire gain and that the selection rule is nearly immaterial: on a Lloyd-stable configuration any facility removal exposes a comparable residual marginal g(x)g(x), so a linear policy over per-facility features operates on an intrinsically low-information quantity and can recover little. (Indeed the cross-entropy optimizer independently assigns large negative weight to the contribution feature, rediscovering the least-contributor rule from data.) Closing the remaining gap therefore likely requires policies that act on joint destroy-and-repair moves or on the sequence of relocations, not on per-facility scores. 8.4 When should FBM-LNS be used? Table 13 summarizes the trade-offs of FBM-LNS, so that practitioners can judge its suitability for a given application. Table 13: Strengths and costs of FBM-LNS. Strengths Costs and limitations Unified: one code handles k-means, p-median, gradual coverage, and shape demand via a single switch. Slower than baselines (∼ 10 s vs. 11–33 s at n=2000n=2000) due to the LNS relocate loop; reducible by lowering iteration budgets at a modest quality cost. Theory-grounded: every component maps to a proven property (Table 2); the greedy provides a (1−1/e)(1-1/e) lower-bound guarantee. No global-optimality guarantee for FBM-LNS itself; it is a heuristic whose quality is validated empirically. Robust: significant superiority across 4 geographies, 3 decay families, 2 scales, and a heterogeneity sweep (30/3030/30, p<10−9p<10^-9). The advantage is specific to smooth (linear/exp) decay; for binary step decay the greedy is competitive (§6). Extensible: the framework accepts any non-increasing per-point function; shape demand is handled by sampling. The learned selector (CEM) improves performance only marginally (+0.22%+0.22\%); richer policies remain future work. In summary, the recommendation splits along the smooth/non-smooth line exposed in §8.1 and the structure/selection line exposed above. When the decay is smooth and heterogeneous—the setting that defines V1–V2 and motivates the per-point Weber weights of Theorem 3.4—FBM-LNS is the recommended method, because the gradient-as-force step is informative and the relocate exploits the widening marginal gaps that heterogeneity creates (Result 3). When a single implementation must also serve k-means, p-median, or shape demand (V3–V5), the same code is competitive with bespoke SOTA (Result 2). When the objective is binary (step decay), the gradient channel closes and the greedy, which rides only the submodular marginal, is the appropriate tool. When runtime dominates, halving the LNS budget recovers most of the gain at half the cost (§8.2), because each relocate iteration is an independent trial whose marginal contribution diminishes. The unifying principle is that FBM-LNS’s two pillars map onto the two information channels—continuous (gradient) and combinatorial (submodular marginal)—and its advantage is largest exactly when both channels are open and most differentially informative. 9 Conclusion We gave a structural theory for heterogeneous distance-decay facility location in the continuous plane and a tractability classification that is the paper’s centerpiece: the discrete problem is always monotone submodular (so the (1−1/e)(1-1/e) greedy guarantee holds regardless of decay shape or heterogeneity), while the continuous cooperative objective is concave—and hence free of local-optima traps—if and only if the decay is concave in the distance, a dichotomy that identifies the clip max(0,⋅) (0,·) as a key mechanism that destroys continuous concavity in common coverage specifications; the classification is tight: concavity holds if and only if the decay is concave in distance. A tight LP relaxation makes the discretized problem exactly solvable in seconds; our force-as-gradient / large-neighborhood-search heuristic is within 0.5%0.5\% of the discrete optimum by fine-grid convergence and outperforms the (1−1/e)(1-1/e) greedy, Cooper ALT, PSO, and weighted k-means across distributions, decay families, scales, and heterogeneity levels; and on real urban-delivery data, ignoring the density-derived heterogeneity costs measurable demand and relocates facilities by up to 37%37\% of the map. Three directions are left for future work: a tighter continuous upper bound than the discretization bound of Proposition 5.2 (for example a conic relaxation of the semi-infinite formulation exploiting the Weber substructure); a head-to-head comparison against modern problem-specific solvers (variable neighborhood search and mixed-integer second-order cone (MISOCP) methods for p-median, dedicated Huff-model solvers), which our classical baselines do not include; and richer learned relocate policies, since the cross-entropy selector of §8.3 already yields a small but significant gain over the hand-designed rule, suggesting headroom for contextual-bandit or attention-based selectors that act on joint destroy-and-repair moves. 10 Code and Data Disclosure The Python implementation of all algorithms studied in this paper—the proposed FBM-LNS solver and every baseline (the (1−1/e)(1-1/e) greedy, multi-start Cooper alternating location–allocation, particle swarm optimization, weighted k-means, and the exact branch-and-bound and LP-bound procedures)—together with the instance generators, the benchmark data sets, and a README giving full reproduction instructions for every table and figure, are available at https://github.com/AgentLabCn/FBM-LNS. The synthetic instances are deterministic functions of a random seed and can be regenerated exactly; the real-world delivery and retail data sets are the public 2026 Meituan Business Analytics Elite Challenge data and the retail panel described in the real-world study. No data or code exemption is requested. Appendix A Proofs of Structural Results This Electronic Companion gives the detailed, self-contained proofs of the theorems, propositions, and lemmas stated in the main paper. We use the following notation throughout: di()≔minj∈J‖pi−Xj‖d_i( X) _j∈ J\|p_i-X_j\| is the nearest-facility distance to point i; Ψi()≔wimaxj∈Jϕi(‖pi−Xj‖) _i( X) w_i _j∈ J _i(\|p_i-X_j\|) is the captured demand at point i; K≔convpii∈IK \p_i\_i∈ I is the demand convex hull; and π(i)∈argminj∈J‖pi−Xj‖π(i)∈ _j∈ J\|p_i-X_j\| is i’s serving facility. We assume each ϕi:ℝ+→[0,1] _i:R_+→[0,1] is non-increasing with ϕi(0)=1 _i(0)=1. Proof of Theorem 3.1 (Nearest-distance reduction). Fix i∈Ii∈ I and write d⋆≔minj∈J‖pi−Xj‖d _j∈ J\|p_i-X_j\|, with the minimum attained at some j⋆∈Jj ∈ J. For every j∈Jj∈ J we have ‖pi−Xj‖≥d⋆\|p_i-X_j\|≥ d , so by monotonicity of ϕi _i, ϕi(‖pi−Xj‖)≤ϕi(d⋆) _i(\|p_i-X_j\|)≤ _i(d ). Equality holds at j=j⋆j=j , whence maxj∈Jϕi(‖pi−Xj‖)=ϕi(d⋆)=ϕi(minj∈J‖pi−Xj‖), _j∈ J _i(\|p_i-X_j\|)= _i(d )= _i\! ( _j∈ J\|p_i-X_j\| ), and the argmax is precisely argminj∈J‖pi−Xj‖ _j∈ J\|p_i-X_j\|, i.e., π(i)=j⋆π(i)=j . Multiplying by wiw_i and summing over i∈Ii∈ I yields Z()=∑i∈Iwiϕi(di())Z( X)= _i∈ Iw_i\, _i(d_i( X)). □ □ Proof of Theorem 3.2 (Convex hull containment). K is closed and convex. Suppose, for contradiction, that an optimal ⋆ X has a facility Xj⋆∉KX_j ∉ K, and let Xj′≔projK(Xj⋆)X_j _K(X_j ) be its metric projection onto K. By the projection theorem for closed convex sets, ⟨Xj⋆−Xj′,pi−Xj′⟩≤0for all pi∈K. X_j -X_j ,\;p_i-X_j ≤ 0 all p_i∈ K. Decomposing pi−Xj⋆=(pi−Xj′)+(Xj′−Xj⋆)p_i-X_j =(p_i-X_j )+(X_j -X_j ) and expanding, ‖pi−Xj⋆‖2 \|p_i-X_j \|^2 =‖pi−Xj′‖2+‖Xj′−Xj⋆‖2 =\|p_i-X_j \|^2+\|X_j -X_j \|^2 −2⟨Xj⋆−Xj′,pi−Xj′⟩ -2 X_j -X_j ,\,p_i-X_j ≥‖pi−Xj′‖2+‖Xj′−Xj⋆‖2 ≥\|p_i-X_j \|^2+\|X_j -X_j \|^2 ≥‖pi−Xj′‖2, ≥\|p_i-X_j \|^2, so ‖pi−Xj⋆‖≥‖pi−Xj′‖\|p_i-X_j \|≥\|p_i-X_j \| for every i, and monotonicity of ϕi _i gives ϕi(‖pi−Xj⋆‖)≤ϕi(‖pi−Xj′‖) _i(\|p_i-X_j \|)≤ _i(\|p_i-X_j \|). Replacing Xj⋆X_j by Xj′X_j therefore does not decrease any Ψi _i and hence does not decrease Z(⋆)Z( X ). Applying the same argument to every facility outside K yields an optimal solution in KpK^p. □ □ Proof of Theorem 3.3 (Monotone submodularity) and the curvature bound. Fix a finite candidate set C and write ZD(S)≔∑i∈Iwimaxc∈Sϕi(‖pi−c‖)Z_D(S) _i∈ Iw_i _c∈ S _i(\|p_i-c\|) for S⊆CS C. Step 1 (Monotonicity). For S⊆T⊆CS T C and any i, maxc′∈Sϕi(‖pi−c′‖)≤maxc′∈Tϕi(‖pi−c′‖) _c ∈ S _i(\|p_i-c \|)≤ _c ∈ T _i(\|p_i-c \|); summing over i gives ZD(S)≤ZD(T)Z_D(S)≤ Z_D(T). Step 2 (Submodularity). Let S⊆T⊆CS T C and c∈C∖Tc∈ C T. Define the pointwise marginal Δi(S)≔wi[ϕi(‖pi−c‖)−maxc′∈Sϕi(‖pi−c′‖)]+. _i(S) w_i [ _i(\|p_i-c\|)- _c ∈ S _i(\|p_i-c \|) ]^+. Because maxc′∈Sϕi(‖pi−c′‖)≤maxc′∈Tϕi(‖pi−c′‖) _c ∈ S _i(\|p_i-c \|)≤ _c ∈ T _i(\|p_i-c \|), the bracketed term is non-increasing in the set, so Δi(S)≥Δi(T) _i(S)≥ _i(T). Summing over i, ZD(S∪c)−ZD(S) Z_D(S∪\c\)-Z_D(S) =∑iΔi(S)≥∑iΔi(T) = _i _i(S)≥ _i _i(T) =ZD(T∪c)−ZD(T), =Z_D(T∪\c\)-Z_D(T), which is the defining inequality of submodularity. Step 3 (Curvature bound). The total curvature of ZDZ_D (Conforti and Cornuéjols 1984) is κ≔ 1−minc∈CZD(C)−ZD(C∖c)ZD(c)−ZD(∅)∈(0,1],κ\; \;1- _c∈ C\, Z_D(C)-Z_D(C \c\)Z_D(\c\)-Z_D( )\;∈\;(0,1], where the range uses monotonicity (Step 1) and submodularity (Step 2). The classical result of Conforti and Cornuéjols (1984) then gives ZD(Sgreedy)≥1−e−κZD⋆Z_D(S_greedy)≥ 1-e^-κ\,Z_D , and (1−e−κ)/κ>1−1/e(1-e^-κ)/κ>1-1/e whenever κ<1κ<1. □ □ Proof of Theorem 3.4 (Weber reduction) and Corollaries 1–2. Fix a Voronoi cell S⊆IS I and the linear decay ϕi(r)=[1−r/Ri]+ _i(r)=[1-r/R_i]^+. The within-cell single-facility objective is maxX∑i∈Swi(1−‖pi−X‖Ri) _X _i∈ Sw_i (1- \|p_i-X\|R_i ) =∑i∈Swi⏟constant = _i∈ Sw_i_constant −minX∑i∈SwiRi‖pi−X‖. - _X _i∈ S w_iR_i\,\|p_i-X\|. Setting αi≔wi/Ri _i w_i/R_i, the right-hand minimization is the weighted Weber problem minX∑i∈Sαi‖pi−X‖ _X _i∈ S _i\|p_i-X\|, proving Theorem 3.4. Corollary 3.5 (Heterogeneity): the first-order optimality condition is ∑i∈Sαi(X−pi)/‖X−pi‖= _i∈ S _i(X-p_i)/\|X-p_i\|= 0, so smaller RiR_i (larger αi _i) gives proportionally stronger pull. Corollary 3.6 (Uniqueness): each X↦‖pi−X‖X \|p_i-X\| is convex, so the non-negative combination is convex; strict convexity of the Euclidean norm on ℝ2R^2 except along rays through pip_i yields a unique minimizer unless all pii∈S\p_i\_i∈ S and X are collinear. □ □ Proof of Theorem 3.7 (Gradient-as-force). Write Z()=∑i∈Iwimaxj∈Jψij()Z( X)= _i∈ Iw_i _j∈ J _ij( X) with ψij()≔wiϕi(‖pi−Xj‖) _ij( X) w_i _i(\|p_i-X_j\|); each ψij _ij is C1C^1 in XjX_j wherever ϕi _i is. Define the generic set ≔:π(i) is unique for every i;G\; \; \ X:π(i) is unique for every i \; its complement is a finite union of Voronoi cell boundaries and hence has Lebesgue measure zero. Fix ∈ X . Because the active index π(i)π(i) is unique at X, the pointwise maximum maxjψij _j _ij is differentiable there, and the chain rule for a maximum of C1C^1 functions with a unique active index gives ∇Xjmaxk∈Jψik()=∇Xjψi,π(i)(),j=π(i),,j≠π(i). _X_j _k∈ J _ik( X)= cases _X_j _i,π(i)( X),&j=π(i),\\ 0,&j≠π(i). cases Multiplying by wiw_i, summing over i∈Ii∈ I, and using ∇Xjψi,π(i)=wi∇Xjϕi(‖pi−Xj‖) _X_j _i,π(i)=w_i _X_j _i(\|p_i-X_j\|) yields the stated identity in the main paper. □ □ We now establish the two convex-analysis lemmas used in the tractability classification, then prove the classification theorem. Lemma A.1 (Diminishing increment). If g:ℝ+→ℝg:R_+ is concave, then for 0≤s≤t0≤ s≤ t and a≥0a≥ 0, g(t+a)−g(t)≤g(s+a)−g(s)g(t+a)-g(t)≤ g(s+a)-g(s). Proof The secant slopes of a concave function are non-increasing: for x1≤x2≤x3x_1≤ x_2≤ x_3, concavity gives g(x2)−g(x1)x2−x1≥g(x3)−g(x2)x3−x2 g(x_2)-g(x_1)x_2-x_1≥ g(x_3)-g(x_2)x_3-x_2. Equivalently the one-sided derivative g+′g _+ is non-increasing, so g(x+a)−g(x)=∫x+ag+′(u)ug(x+a)-g(x)= _x^x+ag _+(u)\,du is non-increasing in x; since s≤ts≤ t, the integral over [t,t+a][t,t+a] is no larger than over [s,s+a][s,s+a]. Hence g(t+a)−g(t)≤g(s+a)−g(s)g(t+a)-g(t)≤ g(s+a)-g(s). □ □ Lemma A.2 (Concave-antitone composition). If ϕ:ℝ→ℝφ:R is concave and non-increasing and h:ℝd→ℝh:R^d is convex, then ϕ∘hφ h is concave. Proof For λ∈[0,1]λ∈[0,1], convexity of h gives h(λx+(1−λ)y)≤λh(x)+(1−λ)h(y)h(λ x+(1-λ)y)≤λ h(x)+(1-λ)h(y). Since ϕφ is non-increasing, applying ϕφ preserves the inequality: ϕ(h(λx+(1−λ)y))≥ϕ(λh(x)+(1−λ)h(y)).φ (h(λ x+(1-λ)y) )\;≥\;φ (λ h(x)+(1-λ)h(y) ). Concavity of ϕφ then gives ϕ(λh(x)+(1−λ)h(y))≥λϕ(h(x))+(1−λ)ϕ(h(y))φ(λ h(x)+(1-λ)h(y))≥λφ(h(x))+(1-λ)φ(h(y)). Chaining the two inequalities yields ϕ(h(λx+(1−λ)y))≥λϕ(h(x))+(1−λ)ϕ(h(y))φ(h(λ x+(1-λ)y))≥λφ(h(x))+(1-λ)φ(h(y)), the defining inequality of concavity for ϕ∘hφ h. □ □ Proof of Theorem 3.8 (Tractability classification). (i) Discrete submodularity. The non-cooperative case is Theorem 3.3. For the cooperative case, write ZDg(S)=∑iwig(σi(S))Z_D^g(S)= _iw_ig( _i(S)) with σi(S)=∑c∈Sϕi(‖pi−c‖) _i(S)= _c∈ S _i(\|p_i-c\|); σi _i is modular and, since each ϕi(‖pi−c‖)≥0 _i(\|p_i-c\|)≥ 0, non-decreasing in S. For S⊆TS T and c∉Tc∉ T, set ai≔ϕi(‖pi−c‖)≥0a_i _i(\|p_i-c\|)≥ 0; the per-point marginal is Δig(S)=wi[g(σi(S)+ai)−g(σi(S))] _i^g(S)=w_i[g( _i(S)+a_i)-g( _i(S))]. Since σi(S)≤σi(T) _i(S)≤ _i(T), Lemma 1 gives Δig(S)≥Δig(T) _i^g(S)≥ _i^g(T); summing yields submodularity. Monotonicity follows from g non-decreasing and ai≥0a_i≥ 0. (i,⇐ ) Cooperative concavity. Assume each ϕi _i is concave on ℝ+R_+ and non-increasing, g is concave non-decreasing, and wi≥0w_i≥ 0. The map Xj↦‖pi−Xj‖X_j \|p_i-X_j\| is convex; by Lemma 2, Xj↦ϕi(‖pi−Xj‖)X_j _i(\|p_i-X_j\|) is concave. As a function of X, each term ϕi(‖pi−Xj‖) _i(\|p_i-X_j\|) is concave in XjX_j and constant in XkX_k for k≠jk≠ j, hence concave in X; the sum Ui=∑jϕi(‖pi−Xj‖)U_i= _j _i(\|p_i-X_j\|) is concave, and the composition g∘Uig U_i is concave (concave non-decreasing composed with concave). The non-negative weighted sum Zg=∑iwig(Ui)Z^g= _iw_i\,g(U_i) is concave. Under concavity every local maximum is global, and projected gradient ascent converges. (i,⇒ ) Necessity (contrapositive). If some ϕi0 _i_0 is not concave on [0,M][0,M], there exist r1<r2<r3r_1<r_2<r_3 in [0,M][0,M] violating the concavity inequality. Set p=1p=1, n=1n=1, g=idg=id, wi0=1w_i_0=1; then Z(X)=ϕi0(‖pi0−X‖)Z(X)= _i_0(\|p_i_0-X\|). Along the ray X=pi0+tuX=p_i_0+t\,u (‖u‖=1\|u\|=1), Z=ϕi0(t)Z= _i_0(t), which is not concave. Hence ZgZ^g is not concave. (i) Non-cooperative non-concavity for p≥2p≥ 2. Take one demand point at the origin with w=1w=1, p=2p=2 facilities, and any non-constant non-increasing ϕφ with ϕ(0)=1φ(0)=1; choose a>0a>0 with ϕ(a/2)<1φ(a/2)<1. Consider =((a,0),(0,0)) X=((a,0),(0,0)): facility 2 sits at the point, so Z()=max(ϕ(a),ϕ(0))=1Z( X)= (φ(a),φ(0))=1. Likewise Z()=1Z( Y)=1 for =((0,0),(a,0)) Y=((0,0),(a,0)). Their midpoint =12(+)=((a/2,0),(a/2,0)) M= 12( X+ Y)=((a/2,0),(a/2,0)) places both facilities at distance a/2a/2, so Z()=ϕ(a/2)<1=(Z()+Z())/2Z( M)=φ(a/2)<1=(Z( X)+Z( Y))/2, violating concavity at the midpoint. □ □ Proof of Proposition 3.9 (Lloyd fixed point). At a fixed point ⋆ X of the alternating map, the location step makes no update, so each Xj⋆X_j satisfies Xj⋆∈argmaxX∈ℝ2∑i:π(i)=jwiϕi(‖pi−X‖),X_j ∈ _X ^2\, _i:\,π(i)=jw_i\, _i(\|p_i-X\|), equivalently 0∈∂C,XjZ0∈ _C,X_jZ (Clarke subdifferential with respect to the single-facility coordinate). The allocation step makes no reassignment, so π(i)=argminj‖pi−Xj⋆‖π(i)= _j\|p_i-X_j \| for every i, confirming Voronoi stability. Together, the two conditions state that no single facility can be moved to strictly increase Z while the others and the assignment are held fixed; this is coordinate-wise local optimality. □ □ Proof of Proposition 3.10 (Monotone relocate). Let (t)t≥0\ X^(t)\_t≥ 0 be the sequence generated by the relocate step under strict-improvement acceptance, Z((t+1))>Z((t))Z( X^(t+1))>Z( X^(t)). Then Z((t))\Z( X^(t))\ is strictly increasing. Once the relocate draws its candidate from a finite set C, the configuration space CpC^p is finite, so a strictly increasing sequence on the finite set Z():∈Cp\Z( X): X∈ C^p\ terminates in finitely many steps. At the terminal configuration X: (i) the Lloyd location step has converged, so X is Lloyd-stable by Proposition 3.9; and (i) no single-facility re-placement from C strictly increases Z, by the strict-improvement acceptance rule. □ □ Proof of Proposition 5.2 (Continuous upper bound via grid density). Let C be a grid of spacing h covering K=convpiK=conv\p_i\, so for every x∈Kx∈ K there is a grid node x~∈C x∈ C with ‖x−x~‖≤h/2\|x- x\|≤ h/ 2. Fix any continuous configuration =(X1,…,Xp) X=(X_1,…,X_p) and let X~j∈C X_j∈ C be a nearest grid node to XjX_j. For each demand point i, write dij=‖pi−Xj‖d_ij=\|p_i-X_j\| and d~ij=‖pi−X~j‖ d_ij=\|p_i- X_j\|. The bound |maxjaj−maxjbj|≤maxj|aj−bj|| _ja_j- _jb_j|≤ _j|a_j-b_j|, the Lip(ϕi)Lip( _i)-Lipschitz continuity of ϕi _i (for linear and exponential decay Lip(ϕi)=1/RiLip( _i)=1/R_i), and maxj‖Xj−X~j‖≤h/2 _j\|X_j- X_j\|≤ h/ 2 give |maxjϕi(dij)−maxjϕi(d~ij)|≤Lip(ϕi)h/2. | _j _i(d_ij)- _j _i( d_ij) |\;≤\;Lip( _i)\,h/ 2. Multiplying by wiw_i and summing over i gives Z()≤Z(~)+(h/2)LϕZ( X)≤ Z( X)+(h/ 2)L_φ with Lϕ=∑iwiLip(ϕi)L_φ= _iw_i\,Lip( _i). Since ~∈Cp X∈ C^p, Z(~)≤Zgrid-MIP(C)Z( X)≤ Z_grid -MIP(C), so Z()≤Zgrid-MIP(C)+(h/2)LϕZ( X)≤ Z_grid -MIP(C)+(h/ 2)L_φ for every X; taking the supremum over X yields Z∗≤Zgrid-MIP(C)+(h/2)LϕZ^*≤ Z_grid -MIP(C)+(h/ 2)L_φ. Finally Zgrid-MIP(C)≤LP(C)Z_grid -MIP(C) (C) because the LP relaxes the integer program. □ □ Proof of Proposition 5.3 (Grid-free continuous upper bound for concave decays). For each i and X, maxjϕi(‖pi−Xj‖)≤∑jϕi(‖pi−Xj‖) _j _i(\|p_i-X_j\|)≤ _j _i(\|p_i-X_j\|) since ϕi≥0 _i≥ 0. Weighting by wiw_i and summing, Z()≤∑iwi∑jϕi(‖pi−Xj‖)=∑j(∑iwiϕi(‖pi−Xj‖))Z( X)≤ _iw_i _j _i(\|p_i-X_j\|)= _j ( _iw_i _i(\|p_i-X_j\|) ). The bracketed term depends only on XjX_j, so the maximum over ∈Kp X∈ K^p separates and equals p⋅maxx∈K∑iwiϕi(‖pi−x‖)p· _x∈ K _iw_i _i(\|p_i-x\|). For a decay that is concave and non-increasing in the distance (unclipped linear ϕi(r)=1−r/Ri _i(r)=1-r/R_i and quadratic ϕi(r)=1−(r/Ri)2 _i(r)=1-(r/R_i)^2 are both concave in r), Lemma A.2 gives that ϕi(‖pi−x‖) _i(\|p_i-x\|) is concave in x; hence f(x)=∑iwiϕi(‖pi−x‖)f(x)= _iw_i _i(\|p_i-x\|) is concave. A concave function attains its maximum over the compact convex hull K at an extreme point, i.e., a demand point pkp_k. Hence the bound is p⋅maxk∑iwiϕi(‖pi−pk‖)p· _k _iw_i _i(\|p_i-p_k\|), computable in O(n2)O(n^2). □ □ Appendix B Abbreviations Table 14 lists every abbreviation used in the paper, its expansion, and the section or table where it first appears. Each abbreviation is also expanded at its first occurrence in the running text. Table 14: Abbreviations used in the paper. Abbrev. Expansion First use ALT alternating location–allocation §1 CEM cross-entropy method §8.3 CI confidence interval Table 5 FBM force-based metaheuristic §4 FBM-LNS force-based metaheuristic with large-neighborhood search §1.3 FDS finite dominating set §2.3 LP linear programming §1.3 LNS large-neighborhood search §4 MCLP maximum covering location problem §1.1 MGCLP multiple gradual cover location problem §1.2 MIP mixed-integer program §1.3 PSO particle swarm optimization §4.2 RcvR_cv coefficient of variation of the decay scale Ri\R_i\ §6.1 RMSE root-mean-square error §7 SOCP second-order cone program §3.6 Acknowledgments. Omitted for double-anonymous review. References Álvarez-Miranda and Sinnl (2019) Álvarez-Miranda E, Sinnl M (2019) An exact solution framework for the multiple gradual cover location problem. Computers & Operations Research 108:82–96, URL http://dx.doi.org/10.1016/j.cor.2019.04.003. Arthur and Vassilvitskii (2007) Arthur D, Vassilvitskii S (2007) k-means++: The advantages of careful seeding. Proceedings of the 18th Annual ACM-SIAM Symposium on Discrete Algorithms (SODA), 1027–1035, URL http://dx.doi.org/10.1145/1283383.1283494. Bansal and Kianfar (2017) Bansal M, Kianfar K (2017) Planar maximum coverage location problem with partial coverage and rectangular demand and service zones. INFORMS Journal on Computing 29(1):152–169, URL http://dx.doi.org/10.1287/ijoc.2016.0722. Bansal and Shojaee (2020) Bansal M, Shojaee P (2020) Planar maximum coverage location problem with partial coverage, continuous spatial demand, and adjustable quality of service, optimization-Online technical report 8178. Berman et al. (2009) Berman O, Kalcsics J, Krass D, Nickel S (2009) The ordered gradual covering location problem on a network. Discrete Applied Mathematics 157(18):3689–3707, URL http://dx.doi.org/10.1016/j.dam.2009.08.003. Berman et al. (2003) Berman O, Krass D, Drezner Z (2003) The gradual covering decay location problem on a network. European Journal of Operational Research 151(3):474–480, URL http://dx.doi.org/10.1016/S0377-2217(02)00604-5. Church and Murray (2018) Church RL, Murray AT (2018) Location Covering Models: History, Applications and Advancements (Springer), URL http://dx.doi.org/10.1007/978-3-319-99846-6. Conforti and Cornuéjols (1984) Conforti M, Cornuéjols G (1984) Submodular set functions, matroids and the greedy algorithm: Tight worst-case bounds and some generalizations of the Rado–Edmonds theorem. Discrete Applied Mathematics 7(3):251–274, URL http://dx.doi.org/10.1016/0166-218X(84)90003-9. Cooper (1964) Cooper L (1964) Heuristic methods for location-allocation problems. SIAM Review 6(1):37–53, URL http://dx.doi.org/10.1137/1006005. De Boer et al. (2005) De Boer PT, Kroese DP, Mannor S, Rubinstein RY (2005) A tutorial on the cross-entropy method. Annals of Operations Research 134(1):19–67, URL http://dx.doi.org/10.1007/s10479-005-5724-z. Drezner (1994) Drezner T (1994) Optimal continuous location of a retail facility, facility attractiveness, and market share: An interactive model. Journal of Retailing 70(1):49–64, URL http://dx.doi.org/10.1016/0022-4359(94)90028-0. Drezner et al. (2004) Drezner Z, Wesolowsky GO, Drezner T (2004) The gradual covering problem. Naval Research Logistics 51(6):841–855, URL http://dx.doi.org/10.1002/nav.20021. Feige (1998) Feige U (1998) A threshold of lnn n for approximating set cover. Journal of the ACM 45(4):634–652, URL http://dx.doi.org/10.1145/285055.285059. Fränti and Sieranoja (2018) Fränti P, Sieranoja S (2018) K-means properties on six clustering benchmark datasets. Applied Intelligence 48(12):4743–4759, URL http://dx.doi.org/10.1007/s10489-018-1238-7. Hansen and Mladenović (1997) Hansen P, Mladenović N (1997) Variable neighborhood search for the p-median. Location Science 5(4):207–226, URL http://dx.doi.org/10.1016/S0966-8349(98)00030-8. Huangfu and Hall (2018) Huangfu Q, Hall JAJ (2018) Parallelizing the dual revised simplex method. Mathematical Programming Computation 10(1):119–142, URL http://dx.doi.org/10.1007/s12532-017-0130-5. Huff (1964) Huff DL (1964) Defining and estimating a trading area. Journal of Marketing 28(3):34–38, URL http://dx.doi.org/10.1177/002224296402800307. Karasakal and Karasakal (2004) Karasakal O, Karasakal EK (2004) A maximal covering location model in the presence of partial coverage. Computers & Operations Research 31(9):1515–1526, URL http://dx.doi.org/10.1016/S0305-0548(03)00105-9. Kennedy and Eberhart (1995) Kennedy J, Eberhart R (1995) Particle swarm optimization. Proceedings of the IEEE International Conference on Neural Networks (ICNN’95), volume 4, 1942–1948, URL http://dx.doi.org/10.1109/ICNN.1995.488968. Lloyd (1982) Lloyd SP (1982) Least squares quantization in PCM. IEEE Transactions on Information Theory 28(2):129–137, URL http://dx.doi.org/10.1109/TIT.1982.1056489. Lu and Zhou (2016) Lu Y, Zhou H (2016) Statistical and computational guarantees of Lloyd’s algorithm and its variants, arXiv:1612.02099. Nemhauser et al. (1978) Nemhauser GL, Wolsey LA, Fisher ML (1978) An analysis of approximations for maximizing submodular set functions—I. Mathematical Programming 14(1):265–294, URL http://dx.doi.org/10.1007/BF01588971. Ostrovsky et al. (2013) Ostrovsky R, Rabani Y, Schulman LJ, Swamy C (2013) The effectiveness of Lloyd-type methods for the k-means problem. Journal of the ACM 59(6):Article 28, URL http://dx.doi.org/10.1145/2395116.2395117. Reinelt (1991) Reinelt G (1991) TSPLIB—a traveling salesman problem library. ORSA Journal on Computing 3(4):376–384, URL http://dx.doi.org/10.1287/ijoc.3.4.376. Weiszfeld (1937) Weiszfeld E (1937) Sur le point pour lequel la somme des distances de n points donnés est minimum. Tohoku Mathematical Journal 43:355–386. Yang et al. (2019) Yang X, Fang Z, Xu Y, Yin L, Li J, Lu S (2019) Spatial heterogeneity in spatial interaction of human movements—insights from large-scale mobile positioning data. Journal of Transport Geography 78:29–40, URL http://dx.doi.org/10.1016/j.jtrangeo.2019.05.010.