Continuous Covering on Networks: Strong Mixed Integer Programming Formulations

Covering problems are well-studied in the domain of Operations Research, and, more specifically, in Location Science. When the location space is a network, the most frequent assumption is to consider the candidate facility locations, the points to be covered, or both, to be discrete sets. In this work, we study the set-covering location problem when both candidate locations and demand points are continuous sets on a network. This variant has received little attention, and the scarce existing approaches have focused on particular cases, such as tree networks and integer covering radius. Here we study the general problem and present a Mixed Integer Linear Programming formulation (MILP) for networks with edges' lengths no greater than the covering radius. The model does not lose generality, as any edge not satisfying this condition can be partitioned into subedges of appropriate lengths without changing the problem. We propose a preprocessing algorithm to reduce the size of the MILP, and devise tight big-$M$ constants and valid inequalities to strengthen our formulations. Moreover, a second MILP is proposed, which admits edges' lengths greater than the covering radius. As opposed to existing formulations of the problem (including the first MILP proposed herein), the number of variables and constraints of this second model does not depend on the lengths of the network's edges. This second model represents a scalable approach that particularly suits real-world networks, whose edges are usually greater than the covering radius. Our computational experiments show the strengths and limitations of our exact approach on both real-world and random networks. Our formulations are also tested against an existing exact method.


Introduction
Covering in Operations Research refers to the optimization problem of deciding the location of facilities to "cover" the points of the so-called demand set, which should fall within the radius coverage of at least one of the installed facilities. This classic problem finds applications in many different domains, including health care [1], surveillance of transport networks [2], computer networks security [3], crane location for construction [4], military evacuation systems [5], homeland defense [6], and urban air mobility [7].
Covering problems have taken many forms in the literature. A rough classification distinguishes between maximal covering location and set-covering location problems. The former aims at maximizing the covered demand with a fixed number of facilities (see e.g. [8]), while the latter seeks to minimize the number of Figure 1: Two cycle coverage points with respect to a cycle C of five nodes variants of the CSCP δ , which rely on different assumptions on the network and the covering radius. Here, we review and compare these FDS. Since we aim to propose a general exact algorithm, we show that it may not be viable to extend discretization methods to solve the CSCP δ for general networks and real radii.
First, Church and Meadows [18] studied the problem with demand at nodes, and identified the following points: The authors proved that F DS 1 := V ∪ N IP is an FDS for the network set-covering problem when the set of demand points is V and that of candidate locations is C(N ).
Secondly, Gurevich et al. [15] studied the continuous set-covering problem when the covering radius and the edge's lengths are natural numbers. They presented an FDS for the case of all edge lengths being one, which can be easily extended to the case of general edge lengths (see [16]), for some e ∈ E and v ∈ e; i = 0, . . . , 2 · l e .
Note that F DS 2 depends on the edge's length. Lastly, Fröhlich et al. [16] proposed a different FDS for the same version of the continuous set-covering with natural numbers. The authors defined the following set of cycle coverage points: CCP := p ∈ C(N ) : d(p, C) := min y∈C {d(p, y)} = δ − l C 2 mod δ, p / ∈ C, for a simple cycle C ⊆ C(N ) , where l C is the total length of the cycle C. Suppose that a cycle C is covered by a set of facilities. A cycle coverage point is the furthest point where a facility that contributes to cover C can be moved without compromising the coverage of the cycle (if the rest of the facilities remain unchanged).  Figure 1 also depicts two locations in C (marked with symbols 'x'), which correspond to two possible feasible locations for the remaining facility needed to cover C (note that the one at the bottom only yields a covering of the cycle if p 2 is located, while the other one together with either p 1 or p 2 can completely cover C).
The authors of [16] gave the following recursive definition of an FDS for the problem with natural numbers: S j+1 := S j ∪ {p ∈ C(N ) : d(p, y) = δ for some y ∈ bd(A(S j ))}; F DS 3 := S |J| , Figure 2: An instance of CSCP δ such that not all facilities in P * are at a distance δ from some node where bd(A(S j )) is the boundary of the area covered by S j , and J ⊆ E is the subset of edges to be covered.
As the author explained themselves, F DS 3 ⊆ F DS 2 . However, the cardinality of F DS 3 may be exponential in the input size, as the number of cycles in a network is in general exponential.
The example depicted in Figure 2 illustrates that none of F DS 1 and F DS 2 are FDS for the CSCP δ . A similar observation was already presented in [13] for a related problem. The figure shows eight nodes on a path, where all edges have equal length. If l e = 1 for all e ∈ E and δ = 1.2, P := {p 1 , p 2 , p 3 } is an optimal δ-cover.
It can be easily observed that there is no optimal solution in which p 2 is placed either at a node or at a distance δ from some of the eight nodes, which shows that F DS 1 is not a valid FDS. On the other hand, it is also easy to check that there is not a feasible solution with the three facilities located either at nodes or middle-points of edges, which proves that F DS 2 is also not a valid FDS. As opposed to F DS 2 , the assumption of the edge lengths and coverage radius being natural numbers is not fundamental in the definition of F DS 3 .
Conversely, F DS 3 is based on the idea of identifying those points at the "boundaries" of coverage areas, i.e., those delimiting the transition from covering/not covering a specific part of the network.
Consequently, F DS 3 could be extended to the general CSCP δ . However, such an extension potentially yields sets with many more candidates, due to the recursive construction of F DS 3

Contribution
As opposed to discretization methods, we directly tackle the CSCP δ for general networks and real radii.
Our main contribution is an exact integer programming approach for the CSCP δ , together with tailored algorithms and strategies to tackle it. Even if this problem has been known for decades, surprisingly, only a few partial results are known for some special cases and sub-classes of networks. To the best of our knowledge, only one MILP model [16] has been proposed so far which can address the general CSCP δ . Such a model can be applied to any network whose edges do not measure more than the covering radius. This condition does not restrict the applicability of the MILP in [16], as any network can be transformed into an equivalent one that satisfies it.
Here, we present an enhanced MILP formulation that relies on the same assumption as that in [16], but whose numbers of constraints and variables have smaller order of magnitude. In addition, preprocessing strategies to reduce the number of variables of the model are studied, and tailored algorithms are presented.
Approaches to strengthen this formulation are also presented, including big-M constants tightening and valid inequalities. The valid inequalities are constraints that reduce the feasible space without removing model solutions.
The introduction of a second MILP, which is scalable concerning the edge's lengths completes the main contributions of the paper. This second MILP is an adaptation of the first one we propose, with the difference that it does not require all edge lengths to be smaller than the covering radius. Finally, our computational experiments prove that the MILP model in [16] is not scalable. On the other hand, the preprocessing technique drastically reduces the size of the first model proposed herein.
Finally, we show in the experiments that the second model we propose is superior to both the model from [16] and our first model, in terms of the solution quality and solving time.
In the proposed setting, both the candidate facility locations and the demand points are continuous sets (in particular, they coincide with C(N )). The problem could be defined for a subset of demand edges, J ⊆ E, and/or a subset of candidate locations H ⊆ E. The theoretical results, model, and methods described in this paper apply to such cases, after straightforward adaptation.

Outline
The rest of the paper is organized as follows. Section 2 presents useful notation and the theoretical development upon which our model is built. Then, our first MILP model is introduced in Section 3, while strategies to strengthen this model are described in the next section. The network processing algorithms that complement our MILP are detailed in Section 5. A second MILP model, which we call reduced formulation and is a modification of the first MILP, is presented in Section 6. Finally, Section 7 describes our computational experiments, and reports and analyzes the obtained results. Section 8 closes the paper with some conclusions.

Covering characterization
This section presents several notation, definitions, assumptions, observations, and results related to the CSCP δ . On the one hand, Proposition 2.1 gives a characterization of the δ-covers of a network, which is based on the individual coverage of each edge of the network. Then, this result is refined to obtain a second necessary and sufficient covering condition in Proposition 2.2. It distinguishes between two alternative possibilities for covering each edge, namely complete or partial, and will be useful for our MILP formulation and methods.
The rest of the section is oriented to characterize the so-called partial and complete covers. The idea of these sets is to delimit the areas of the network where a facility, if placed, would completely cover a given edge, and those where a facility would reach the edge (but maybe not completely covering it).
We first introduce some related notation, definitions and assumptions. We assume that V is totally ordered by the binary relation . Every edge e ∈ E has a unique representation, and v a v b . From now on, we take e = (v a , v b ) indifferently as a continuum in C(N ) or as an edge ending at v a , v b . We extend the edges' length function to l : C(N ) → R + as a length measure on the continuum of points. For two points p, p ∈ C(N ), we denote by Π(p, p ) ⊆ 2 C(N ) and Π * (p, p ) ⊆ Π(p, p ) the set of paths and shortest paths, respectively, connecting p and p . Any path π ∈ Π(p, p ) is indifferently treated as a continuum in C(N ), then l π := l(π) is the length of π. The distance between p and p , d(p, p ), is the length of a shortest path connecting them: d(p, p ) := min{l π : π ∈ Π(p, p )} = l π * for any π * ∈ Π * (p, p ).
In particular, if p and p belong to the same edge, we denote by l(p, p ) the length of the unique path in that edge connecting them. We work under the following assumption: If Assumption 2.1 did not hold, we could consider a set I ⊂ 2 C(N ) , that would contain, for each e = (v a , v b ) ∈ E, the following continuum sets of points (segments): • If δ ≥ l e , e ∈ I; • If δ < l e , let n := le δ + 1. We We consider N = (V , E = I), where V contains the endpoints of I. The new network N satisfies that δ ≥ l e for all e ∈ E , and it is isomorphic to N with respect to the length function. Indeed, since N is obtained by subdividing edges in N , C(N ) = C(N ) and a set of points δ-covers N if and only it δ-covers N . Therefore, Assumption 2.1 always holds after the network N is transformed into N (via a preprocessing step). Such transformation yields a network with more nodes and edges, which has a direct impact on the size of optimization models. In Section 6, we present a model that avoids this effect.

Observations
In the following, we give several observations of optimal δ-covers, which guide our quantitative analysis of covering conditions and the resulting MILP model of Section 3.
Observation 2.1. For an edge e = (v a , v b ) ∈ E and a point p ∈ C(N ), one of the following cases holds: 1. p cannot δ-cover any point in e; 2. p can δ-cover the whole e; 3. p can δ-cover a continuous portion of e containing either v a or v b ; 4. p can δ-cover two continuous portions of e, which do not intersect, each contains either v a or v b .
Observation 2.2. (a similar statement was proven in [16]) There exists an optimal δ-cover that satisfies: i) Each edge e ∈ E has at most two facilities (due to Assumption 2.1); ii) If there are two facilities in the edge e, we can assume without loss of generality that they are located at end nodes v a , v b (this follows from i) and Assumption 2.1).
As a consequence, the set of candidate facilities of a δ-cover is in one-to-one correspondence to the edges and nodes of the network.
With the above observations, we can already give a high-level description of the covering characterization behind our model. Namely, if we fix a set of facilities on some network edges, there would be some edges completely covered regardless of the exact facility locations within their edges. Some other edges would be partially covered from the left-end node and/or from the right-end node, and how much depends on the actual facility locations. So we have variables that specify the exact facility locations on each edge. Finally, we stipulate that the cover from the left and the cover from the right better exceed the edge length. One difficulty is that the corresponding covering function is not linear However, we show that such a function is a piece-wise linear function, which can be modeled by a MILP. The aim of the remainder of this section is to give a mathematical specification of the above characterization.

Covering conditions
We give a sufficient and necessary condition that the network is δ-covered by installed facilities.
δ-covered by P if and only if either there exists p ∈ P ∩ e or Moreover, the set P is a δ-cover of N if and only if for each e ∈ E, either there exists p ∈ P ∩ e or (2) is satisfied.
Proof. If there exists p ∈ P ∩ e, then e is δ-covered by p due to Assumption 2.1. Otherwise, for each i ∈ {a, b}, (2) can be rewritten as follows: Note that max{δ − l π * i , 0} represents the maximum length that can be δ-covered by P (specifically, from p * i ) after passing through v i . Since the path(s) that δ-cover e must contain v a and/or v b , the edge is covered if and only if these "maximum lengths" for v a and v b add up to more than l e . and v b respectively. In this case, the min functions inside (2) are attained respectively at p and p . In the middle, Figure 3b shows a single location p that covers e through two different paths, which traverse v a and v b respectively. These paths form a cycle that contains p and e. In this case, the two min operations inside (2) are attained at the same point, p. Finally, Figure 3c illustrates the case in which a single facility located at p covers e through one of its end nodes, v a . Here, one of the max operators in (2) is equal to zero (p is further from v b than δ).

Covering delimitation and simplification
For an optimization or search problem, delimitation refers to the reduction of the candidate space. As we have seen in Section 1.1, FDS is studied in related works as a way for reduction of CSCP δ (under some assumption), and hence it is a kind of delimitation. Instead of FDS, we consider a different delimitation that can be used for the general CSCP δ , and which allows us to obtain a reduced MILP formulation.
The characterization in Proposition 2.1 is based on the individual covering of every edge in the network.
When considering possible locations to cover a fixed edge, we can restrict ourselves to its surroundings within the radius δ. Delimiting those parts of the network that could "contribute" to covering a particular edge or node reduces the search space. We then introduce three kinds of delimitation: potential covers, complete covers and partial covers. We will represent the covering condition under such delimitation. In effect, the covering condition has a simplified form compared to its general form in Proposition 2.1, and an adequate preprocessing procedure can reduce and strengthen our MILP model. We find that, for every node, there may exist potential covers, i.e., a set of edges and nodes where, if a facility is located, it can possibly δ-cover this node. The potential covers in the following definition delimit the edges and nodes of C(N ) that can contribute to covering a particular node of the network.
Definition 2.1. For each v ∈ V , the potential covers of v are the candidate facility locations to cover v: Clearly, v is not reachable within the radius δ for any facility installed outside F(v). Regarding the covering of edges, an edge incident to v could be covered by some of the facilities in F(v).
We find that, for every edge e, there may exist complete covers, i.e., a set of edges and nodes where, if a facility is located, it can always δ-cover the edge e, regardless of the exact facility location. Once there is a facility within a complete cover, the whole edge e is guaranteed to be covered by this facility. The following definition serves to delimit the network to such complete covers.
If a facility is placed at F c (e) (either at a node in V c (e) or at a point on an edge belonging to E c (e)), we can immediately conclude that e is δ-covered. Note that any facility placed at e can completely cover e if and only if any facility placed at e can completely cover e . That is, E c is symmetric over E. On the other hand, it is obvious that e ∈ E c (e), and v a , v b ∈ V c (e), for all e = (v a , v b ) ∈ E.
Given a node v, we can characterize the complete covers of the incident edges to the node v. This helps us refine the potential covers of this node. The following definition identifies those candidate facility locations in the potential covers of the node v that cannot completely cover any incident edges to v. Definition 2.3. We define the following sets for each v ∈ V : We call these sets the partial covers of E(v), where E(v) := {e ∈ E : v ∈ e} is the set of incident edges to v.
The set F p (v) contains those candidate locations that can contribute to partially covering some of the edges in E(v). Note that, if a facility is placed at F(v) \ F p (v), then this facility completely covers E(v).
Definitions 2.1, 2.2, and 2.3 provide us with a refined covering condition. Indeed, the following proposition is a consequence of Proposition 2.1 and the aforementioned definitions, and will be used to characterize coverings in the MILP formulation presented in Section 3.
Moreover, For the second statement of the proposition, we have it suffices to see that we can take p ∈ P ∩ (E p (v i ) \ V ) in the second inner min operator of the right-hand We prove that the end nodes of e can be excluded from the second inner min operator.
, and we can exclude it from the second inner min operator (this node is already considered by the first inner min operator). We consider now the other end node of e . If d(v i , v ī ) ≤ δ then, similarly, v ī can be excluded from the second inner min.
Otherwise, we know that the outer min is not attained at v ī , as Remark 2.1. The covering conditions described both in Propositions 2.1 and 2.2 would be also applicable if only a subset of E, J ⊆ E, is to be covered. Indeed, these covering conditions are based on the individual coverage of the edges, so it would be sufficient to apply them just to the edges in J. As a consequence, our methods, including the MILP formulation and algorithms presented in the next sections, apply to this more general version of the CSCP δ .
In order to exploit the newly defined potential, complete, and partial covers in our formulation, from a practical viewpoint, we need to have some characterizations that can operate in a computer. Definition 2.1, which introduces potential covers, satisfies this requirement. Indeed, it just depends on distances between pairs of nodes, which we can easily calculate. Conversely, Definition 2.2 presents complete covers with a condition that must be satisfied by "the infinitely many points of an edge", which is not directly computable.
Finally, the elements in the partial covers defined by Definition 2.3 can be easily calculated once both potential and complete covers are known.

Characterization of complete covers
In the following, we focus on characterizing the complete covers, which will be useful for our MILP formulation and tailored algorithms, (see forthcoming Sections 3 and 5). To begin with, we note that Proposition 2.1 already gives us a characterization of the nodes in V c (e). Indeed, it is easy to observe that, for a given e ∈ E, v ∈ V c (e) if and only if P := {v} δ-covers e. We then focus on the sets E c (e). On the one hand, it is clear that, for every edge Moreover, if we define E c (v) ⊆ E(v) as follows, we have a tighter set containing E c (e).
Definition 2.4. The edges that can completely cover a node v ∈ V are: It is clear that the following observation holds.
With this observation, we can limit the search of E c (e) in the set E c (v a ) ∩ E c (v b ), as we will show in Section 5, the latter set is easy to compute.
and constant: which satisfies the following equation: The function d v (q) represents the distance between v and a point p ∈ e such that q = l(v a , p ), where l(v a , p ) measures the length of the continuum (v a , p ) ⊆ e . The inner terms in the minimization that defines d v (q) coincide for q = Q ve . Informally, Q ve is the "bottleneck" coordinate on e (see original definition in [19]), for which the distance to v is the same if we go through it follows that 0 ≤ Q ve ≤ l e . Note that d v (q), r v (q), and Q ve depend also on the edge e .
Proof. e ∈ E c (e) if and only e is δ-covered by any point p ∈ e . Take P = {p } in Proposition 2.1, e is δ-covered by P, if and only if holds for all p ∈ e . Let q = l(v a , p ), q ∈ [0, l e ], be the measure of the sub-edge (v a , p ) ⊆ e . Then, the We present the following lemma without proof. The lemma is a direct consequence of Definition 2.5.
, and it decreases for q ∈ [Q ve , l e ]. Moreover, d v (q) admits the following piece-wise linear representation: Note that with Observation 2.3, to find E c (e), we can check . Then, we consider v ∈ V , and we want to characterize r v for all points on edges e ∈ E such that e ∈ E c (v). We present the following result without proof.
, then r v (q) admits the following piece-wise linear representation: Finally, we have a tractable version of Lemma 2.1.
Proof. From Lemma 2.1, e ∈ E c (e) if and only if min q∈[0,l e ] r(q) ≥ l e . Due to Lemma 2.2, the minimum argument must be some of the breakpoints in the piece-wise linear description of r v (q). Then, it suffices to ≥ l e and r vi (l e ) ≥ l e always holds for all i ∈ {a, b}. It suffices thus to check the following condition Then, the result follows.

MILP formulation
We present the variables of the MILP formulation first. For each v ∈ V , there is one candidate facility (fixed location); and for each e ∈ E, there is another one (location within the interior of e, e \ {v a , v b }). Then, the finite set F := E ∪ V will be used to index the candidate facilities. In our MILP formulation, there are two decisions associated with each f ∈ F. One is to decide if a facility is installed at f . The second is only necessary for those facilities installed at the interior of edges, and consists in determining their locations within the corresponding edges.
To represent the first of the above decisions, we define the following binary variables, which we call the placement variables: We identify the set of installed facilities with F 1 = {f ∈ F : y f = 1}. To represent the second of the above decisions, we define the following continuous variables, which we name the coordinate variables: We use v and e to denote nodes and edges where facilities are installed, and use v and e to denote nodes and edges to be covered, respectively. We refer to v a , v b as the end nodes of e = (v a , v b ) ∈ E; given We use the necessary and sufficient condition of δ-covering in Proposition 2.2, and the second result in this proposition regarding the distance function. Other than the placement and coordinate variables, some additional variables are used, which we present next.
We sometimes refer to r v as the "residual cover" at node v, since it represents the maximum remaining length that can be covered after reaching v from "a sufficiently close" facility. Remembering Definition 2.5, these variables must satisfy: In an optimal solution, it can be r v = r v (q e ) for an edge e ∈ E as stated above. In this case, r v is the maximum remaining length that can be covered after reaching v from the closest facility. However, we do not impose this equality in our formulation, since it is enough for guaranteeing the coverage of e = (v a , v b ) ∈ E that the sum of the residuals r va + r v b exceeds l e (see Proposition 2.2). That is, if l e is already exceeded by r va + r v b for some r vi < r vi (q e ), then the coverage condition of Proposition 2.2 will hold.
We denote by M * the sufficiently large big-M constant associated with index * , the value of which will be determined later. Our formulation of the CSCP δ reads as follows: Constraints (6b) and (6c) model the logic or constraint w e = ∨ f ∈Fc(e) y f . Constraints (6d) and (6e) enforce the x v is the product of the w e variables such that e ∈ E(v). Constraints (6f) prevent two facilities in a solution from being installed respectively at the interior of an edge and one of their end nodes. Constraints (6g) bound the coordinate variables with the corresponding edge length, and set them to zero if no facility is located at its interior. The covering condition in Proposition 2.2 is enforced by (6h). If w e = 1, then the condition is satisfied (e is covered by F c (e)). Otherwise, the inequality (5) of the proposition has to be satisfied. The rest of the constraints of the model (6i)-(6n), together with variables r, x, q, and z, aim at modeling (5). To begin with, (6i) impose that, for each v ∈ V , one of the following statements holds: i) All incident edges to v, e ∈ E(v), are completely covered by facilities placed at their complete covers, iii) A sufficiently close facility to v is installed at e ∈ E p (v) and v is reached through v i of e (z ve i = 1), that is, If the case i) above holds, then the covering condition in (6h) is satisfied for all e ∈ E(v), regardless of the value of the residual cover variables. Otherwise, suppose that x v = 0 and w e = 0 for some e ∈ E(v). In this case, the corresponding constraint (6h) is "active", that is, the inequality (5) of Proposition 2.2 has to be satisfied for e. Since x v = 0, constraints (6i) impose that there is a facility among those installed at F p (v) that is sufficiently close one to v. This facility is the one bounding the residual variables r v (see constraints (6m)-(6n)), which represent the terms in the left-hand side of (5). Constraints (6j) (resp. (6k)) ensure that z vv (resp. z ve i ) can be one only if facility is installed at v (resp. e ). Due to (6i), for every fixed node v ∈ V , at most one of the constraints in (6l)-(6n) will be active. If the covering conditions (6h) are "inactive" and r v is not needed to guarantee the coverage of any e ∈ E(v). presents refined values of these big-Ms.
We observe that the number of variables and constraints in (6) can be reduced. Namely, for each v ∈ V always holds for all p ∈ e . For such nodes and candidate facilities, we do not need both variables, z ve a and z ve b , and corresponding constraints in (6n) (we know beforehand that one of these constraints would never be active if a facility is located at e ). Therefore, EI p (v) would only contain one of the pairs (e , a) or (e , b).

Comparative insights with respect to an existing MILP
To the best of our knowledge, the only existing MILP for the CSCP δ was proposed in [16]. The authors used a similar observation to ours with respect to optimal δ-covers. They noted that every edge contains at most one facility. Indeed, in their setting, if two facilities are located at both end-nodes of an edge e = (v a , v b ) one of them is considered to be "hosted" by an adjacent edge, e ∈ E(v a ) ∪ E(v b ) (by optimality, neither v a nor v b is a leaf). Their location variables are indexed then by E.
However, this approach has issues with symmetry, which leads to redundant solutions. Indeed, many different solutions to the MILP model represent the same facility locations to the CSCP δ , since there exist many combinations of the edges "hosting" the facilities that are located at nodes. Let us consider an example: let v ∈ V be a node, e ∈ E(v) be an incident edge, and p ∈ e be a point on e. Even though the distance d(p, v) is equal to for a very small > 0, the point p is still located at e. However, when d(p, v) = 0, the point p is located at the node v and, consequently, at every edge in E(v). When a facility is located at a node, the discontinuity there leads to the question: which node or edge do we choose to represent this facility? In our MILP model, we have a specific node facility variable (i.e., y v ) which prevents edge facility variables (i.e., y e , e ∈ E(v)) from "hosting" facilities at v, thanks to constraint (6f).
On the other hand, we consider model size in terms of the number of variables. A second main difference between MILP (6) and the MILP in [16] is that the latter uses binary variables to identify the two edges containing the facilities that cover a given edge.  actually correspond to the same characterization of δ-cover, but are modeled in a slightly different way.
Namely, the authors of [16] defined the "residual covers" for each edge (where a facility might be placed) and node of the network. The interested readers might consult [16] and the MILP therein, which we do not reproduce here for the sake of concision.
Nonetheless, Table 1 shows a comparative summary of the two formulations, based on the number of variables and constraints. This summary considers an upper bound on the size of MILP (6). That is, we take and e ∈ E-however, this would never be the case, as the partial and complete covers are complementary. On the other hand, Table 1 considers the MILP in [16] with J = E (the set of edges to be covered).

Strengthening
In this section, we analyze modifications of the MILP (6) that can yield a tighter linear relaxation of this formulation. Namely, we tighten our big-M constraints (6l)-(6n) by devising small constants M v , M vv , M ve i , δ vv , and δ ve i . We also present several families of valid inequalities. Valid inequalities define conditions that have to be satisfied by any feasible solution, and yield tighter linear programming relaxations, see e.g. [20].

Constants tightening
From the MILP formulation, it is easy to yield the following observation. For v ∈ V , it suffices for a facility f ∈ F p (v) to contribute to the residual cover r v at most U v := max e∈E(v) l e . Indeed, r v aims at ensuring that the inequality (5) of Proposition 2.2 is satisfied for all e ∈ E(v). We define Since U v is a valid upper bound for the residual cover variable r v , the big-Ms in the constraints (6m) and (6n) should guarantee that Taking the minimums of the above big-Ms, we can now tighten the big-M constants of the MILP (6) as (a) A facility is located at v (b) Solution with two facilities Figure 4: Illustration of valid inequalities (9) follows: where the last equations in the definition of M vv and M ve i follow from the definition of δ vv and δ ve i , respectively.

Valid inequalities "Leafs" inequalities
If a node v ∈ V has degree one, we can assume without loss of generality that no facility is located at v nor at its incident edge. Indeed, an equivalent δ-cover could be built by just moving such a facility to the unique neighbor of v in N . More than valid inequalities, the following are valid variable elimination: "Adjacent edges" inequalities Consider a node v ∈ V of degree two. If there is a facility at v, then no facility is placed at the edges incident to v (we recall the model constraints (6f)). Otherwise, we can assume that at most one facility is placed at these edges in an optimal solution, which can be enforced by the following valid inequalities: Figure 4 illustrates the above inequalities. Figure 4a shows the case in which a facility is located at v.
Otherwise, if two facilities are placed at e and e respectively, we can build an equivalent solution by moving one of these facilities to the end node of the corresponding edge that is not v, as depicted in Figure 4b. We recall that the last statement holds due to our Assumption 2.1.

"Neighborhood" inequalities
Let us now consider a node v ∈ V , and suppose that there are several facilities placed at different edges in E(v) in a feasible solution. Take e * ∈ E(v) containing a facility f * such that d(f * , v) = min{d(f, v) : (a) A facility is located at v (b) Solution with several facilities Figure 5: Illustration of valid inequalities (10) f is installed at e ∈ E(v)}. The following proposition gives an equivalent feasible solution where the facilities at the edges e ∈ E(v) such that e = e * are moved to the nodes. • y e = 0 for all e ∈ E(v) such that e = e * , andŷ e = 1; • y f =ŷ f otherwise. Otherwise, suppose the facility at e was partially covering e through node v. In this case, the facility at e * covers at least the same part of e (it is closer to v). Since this facility remains unchanged in the new solution, we can guarantee that e is still covered.
As a consequence of Proposition 4.1, the following inequalities are valid: Figure 5 illustrates the valid inequalities (10). In particular, Figure 5b illustrates the equivalent solution given in Proposition 4.1. It is easy to observe that these new inequalities are a generalization of inequalities (9).

Network processing
The network processing algorithm analyzes the network N to compute the parameters and sets needed to construct the MILP model (6), which we recall next: 1. V c (e), E c (e) for all edges e ∈ E; The above data is computed by Algorithms 1, 2 and 3. Algorithms 1 and 2 contain auxiliary functions, which are called within the main Algorithm 3. Algorithm 1 computes the sets E(v) and V(v) (which are not directly used in the MILP but necessary to obtain E p (v) and V p (v)), and the distances Algorithm 1 also computes the sets E c (v), which will serve as intermediate sets to finally obtain E c (e) in Algorithm 3. The main task in Algorithm 3 is to compute the sets V c (e) and E c (e). To that aim, this algorithm calls both Algorithm 1 and the procedure "mutual" described in Algorithm 2. Once V c (e) and E c (e) are known, the computation of V p (v) and E p (v) in Algorithm 3 easily follows by definition.
In the following, we present Algorithm 1, which defines the function "nodeCover ( Algorithm 2 describes the procedure "mutual", which determines, given e = (v a , v b ) ∈ E and a candidate edge for the complete cover e ∈ E c (v a ) ∩ E c (v b ), whether e ∈ E c (e). This algorithm is based on Proposition 2.3 in Section 2.
Network processing Algorithm 3 computes all the sets that are needed by the MILP formulation. The algorithm starts with empty sets E c (e), V c (e), E p (v), V p (v), for e ∈ E and v ∈ V . In Line 3, the algorithm loops through all nodes v ∈ V and computes the function "nodeCover(N, δ, v)", storing its output. Then, the algorithm calculates the sets V c (e) for e ∈ E, by applying the symmetric relation between these sets and the sets E c (v) from "nodeCover(N, δ, v)". After that, in Line 10, the algorithm loops through all edges e ∈ E c (e )) by calling the procedure "mutual". Since e ∈ E c (e) if and only if e ∈ E c (e ), the loop only runs over pairs such that e < e (we assume an total order on the elements of E). The loop starting in line 18, iterates on each node v ∈ V and looks for v ∈ V(v) such that there exists an e ∈ E(v) but e / ∈ E c (v ). The nodes v found are added to V p (v). Finally, the loop in line 28 also iterates on v ∈ V , and looks for e = (v a , v b ) ∈ E(v) such that there exists e ∈ E(v) but e / ∈ E c (e ). Each edge found is added to E p (v), and, right after that, the set EI p (v) may be updated after checking the dominance rule described at the end of Section 3. We have the following complexity result for Algorithm 3. The time complexity of the network processing Algorithm 3 is O(|E| 2 + |V ||E|(D + log |V |)).
Proof. We first analyze the time complexity of the procedure nodeCover described in Algorithm 1. The main while loop is a modification of the Dijkstra algorithm, and it can be implemented with time complexity O((|E| + |V |) log |V |), see [21]. Therefore, the overall time complexity of nodeCover is also The network processing Algorithm 3 has four outer loops, and next we analyze the complexity of each outer loop. The first outer loop runs the nodeCover algorithm over the nodes, so its complexity is Algorithm 1: single node δ-cover algorithm: nodeCover 2 Output: Boolean value indicating whether e ∈ E c (e).

A reduced formulation for networks with long edges
The MILP (6) assumes l e ≤ δ for all e ∈ E, however, in many real-world networks, some edge lengths are greater than the covering radius. To reuse the previous results, one approach is to transform a network with long edges into another network with edge lengths at most δ. The transformation is by subdividing edges of the original network into smaller pieces, so the optimal cover does not change. This transformation enables us to apply MILP (6) or the MILP in [16] on the transformed network. We note that, in [16], there is a recursive transformation. However, the recursive transformation is not constructive, so we cannot employ it in practice.
Another trivial transformation is by subdividing edges as suggested in Section 2. However, this strategy is a trivial heuristic, and it increases the number of edges and nodes of the transformed network, and thus the number of variables and constraints of the MILP model by a nonlinear factor.
In this section, we present an alternative approach to tackle networks with edge lengths greater than the covering radius. Instead of transforming the network, this approach directly treats long edges in the formulation by using specific sets of constraints and variables. We highlight that the approach is also applicable to "long paths". That is, if there is a path in the network with all their middle-nodes having degree two, we can represent it by a single edge of length equal to the total length of the path. Indeed, the CSCP δ does not change after this transformation.
The main idea of the reduced formulation is to assume a predefined covering of those edges that are long enough. Such covering consists in placing facilities every 2δ distance units on the long edge. Let us consider e ∈ E such that l e > 2δ. An edge satisfying this condition is called a long edge. We denote bŷ l e := l e − 2δ l e /(2δ) the length of the last piece of e after dividing it into pieces of measure 2δ. We calll e the tail of e. The following proposition guarantees the correctness of the reduced formulation.
Proposition 6.1. Let N be an undirected network, e = (v a , v b ) ∈ E be a long edge, and P be a feasible δ-cover of N . Define P with p ∈ P for all p ∈ P \e. Let p e ∈ P ∩e be such that l(v a , p e ) = min p ∈P ∩e l(v a , p ), and let q e := l(v a , p e ) (here q e represents a length, although it will also be a variable of the reduced MILP that we introduce afterwards). Note that q e ∈ [0, 2δ] (otherwise, P would not be a δ-cover). The set P can be completed in such a way that it δ-covers N and |P| ≤ |P |, as follows: -p e ∈ P; -p ∈ P for all p ∈ e, p > p e , such that l(p e , p) = 2δ · k, for some k ∈ N; -If ∃p ∈ P ∩ e such that d(p e , p ) > 2δ l(p e , v b )/(2δ) then v b ∈ P.
Proof. It is easy to observe that |P| ≤ |P |. First, since facilities are placed every 2δ distance on e, the original covering P cannot contain fewer facilities than P. Since the rest of the facilities are just taken from P , |P| ≤ |P |. On the other hand, P has to δ-cover N for the same reason. That is, the facilities of P that were not in P cover at least as much as the ones originally in P . The last part of the proposition easily follows from construction, and is illustrated by Figure 6.
In the following, we present our reduced formulation, which is an adaptation of MILP (6). We treat long edges specifically, improving the scalability of our approach. Edges e ∈ E such that δ < l e ≤ 2δ are subdivided into two sub edges of length smaller than δ. Therefore, we assume that, for every e ∈ E, either l e ≤ δ or l e > 2δ. In the former case, all the constraints and variables of the model remain unchanged. In the latter, we introduce new variables and constraints to the model, while dropping some of the constraints originally in (6).
The objective function also needs adaptation. We introduce all these modifications next.
Let e = (v a , v b ) ∈ E be a long edge and, for any feasible solution, let q e be as in Proposition 6.1. That is, when e is a long edge, we use the former variable q e of MILP (6) to represent the position of the left-most facility on e with respect to v a . The placement variables y va , y v b , and y e will be used as well, with slightly different meanings to those in (6), as we will explain later on. We introduce an indicator variable u e ∈ {0, 1} to distinguish between two possible ranges in the domain of q e . If 0 ≤ q e ≤l e , then u e = 0; otherwisê l e ≤ q e ≤ 2δ and u e = 1. This can be modeled with the following constraints: Note that the coefficients on the last term in the objective already account for the facilities installed at v a for each e = (v a , v b ) ∈ L, while they do not do so for v b . This will condition the values of the placement variables in an optimal solution, namely, y va = 0 for all e = (v a , v b ) ∈ L. The facilities installed at these nodes will be tracked by the variables q e , namely, if q e = 0 then a facility would be installed at v a . To complete the modeling of the CSCP δ , we need to ensure that the covering of the long edge fits into the covering of the rest of the network. Namely, some parts of the network might be covered by facilities placed on e, and part of e (namely its tail or the portion between v a and p e of Proposition 6.1) could be covered by facilities placed outside e.
We focus first on the case 0 ≤ q e ≤l e . We need to ensure that both the segment (v a , p e ) and the tail of e are covered. For the tail, we know that there is a facility at a distancel e − q e from v b . This facility covers a length δ of the remaining fragment on e on its right-hand side, which has length equal tol e − q e . The rest of δ = 1 and le = 3.5. Therefore, to cover e, we need 3 facilities (|P ∩ e| = 3). As le/(2δ) + 2 = 3, the upper bound can be tight. Similarly, we can give an example with a tight lower bound, by simply setting le = 2 in the previous network.
such fragment should be covered, which can be imposed by the following constraint: To ensure the covering of the segment (v a , p e ), we have: Let us consider now the casel e < q e ≤ 2δ. There is a facility installed at a distance of 2δ − q e +l e from v b .
Then, to ensure that the tail of e is covered, we need to cover the fragment between this facility and v b . Since the facility already covers a length δ on this fragment, the following constraint enforces the covering of the tail: To ensure the covering of the segment (v a , p e ), we have the same equation as before: In summary, the following constraints ensure that the edge e is fully covered: Constraint (16) gathers (13) and (14) where, if e ∈ L, the term τ ve l i (q e ) in (6n) is replaced by d(v, v i ) + 1 i=a q e + 1 i=b (2δu e +l e − q e ). We enforce (17) because e always contains a facility if e ∈ L. On the other hand, we need to include constraints (18) to guarantee that the variables r va and r v b can take positive values. Indeed, if w e = 1 for e ∈ L, it may happen that x va = 1 or x v b = 1, which will will imply, respectively, r va = 0 or r v b = 0 due to (6l). We compute the complete and partial cover sets in the same way as for the original MILP model. Note that no edge or node can completely cover e if e ∈ L.
The following theorem is on the scalability of the reduced MILP above. Proof. The number of constraints and variables of the reduced model does not grow with the edge lengths, except for a constant factor of 2 for those edges e ∈ E such that δ < l e ≤ 2δ.

Computational results
In this section, we present the computational experiments testing the existing and proposed formulations and strengthening techniques for CSCP δ and its discrete variant (facilities on nodes).

Experiment Setup
We describe the setup of the experiments including the benchmarks, development environment, implementation of algorithms and solution statistics. The computational results and source code are publicly released on our project website: https://github.com/lidingxu/cflg/, where we provide a bash file to reproduce the experiments in Linux systems. Those benchmarks that we generated for this study, or that were publicly available already, are also available at the repository. Benchmarks. We use three different benchmarking sets: two come from the literature, and the other has been generated synthetically. For every instance, we set the coverage radius δ equal to the average of the edge lengths. We describe these benchmarks next.
Kgroup. It consists of 23 prize-collecting Steiner tree problem instances from [22], and the benchmark includes the graphs and edge lengths of these instances. These random geometric instances are designed to have a local structure somewhat similar to street maps. Nodes correspond to random points in the unit square. The number of nodes ranges from 22 to 241. There is an edge between two nodes if their distance is no more than a prescribed threshold which depends on the number of nodes, and the length of an edge is the Euclidean distance between the two points. It is divided into two sets, Kgroup A and Kgroup B. The first one consists of 12 small instances with up to 45 nodes, and the second one consists of 11 large instances with up to 241 nodes.
City. It consists of real data of 9 street networks for some German cities, and it was first used in [13]. The number of nodes ranges from 132 to 771. The length of each edge is the length of the underlying street segment.
Random. It consists of 24 random network instances generated via Erdős-Rényi binomial method with the package "Networkx" (see [23] Problem preprocessing. Networks of instances are modified in a problem preprocessing step to be amenable to MILP models.
Given an original network of each instance, in the first preprocessing step, we delete any degree-two node and concatenate its adjacent edges to a new edge, as long as the deletion does not yield a self-loop. Such a node can be treated as an interior point of the new edge. We refer to the preprocessed network without any such degree-two node as the degree-two-free network.
Even after the first preprocessing step, the degree-two-free network may not correspond to the actual  feasible solution of CSCP δ . We name this discrete restriction by the discrete facility SCP δ . This model solves the discrete facility SCP δ , which solely sets y e = 0 for all e ∈ E in SF model.
The above models are summarized in Table 2. Both EF and F0 consider that any two points in the network can possibly cover each other, and do not utilize the complete and partial cover information. They have been already compared in Section 3.1, and hence F0 should have fewer variables and constraints than EF. We are interested in the dual gaps obtained after the models are solved within the time limit for these models.
Performance metrics and statistical tests. We describe the performance metrics and the ways to compute their statistics. These statistics will be used to evaluate the model performance.
Let v be a dual lower bound and v be a primal upper bound obtained after solving some of the models described above, the relative dual gap is defined as: A smaller relative dual gap indicates better primal and dual behavior of the model.
Let n sd be the number of nodes of the subdivided network of that instance, note that a trivial primal solution is the set of the nodes of the subdivided network (for which edge length is at most δ). Therefore, to normalize the primal solution value, we define the relative primal bound v r := v n sd .
If v r < 1, then the model finds a solution better than the trivial one.
In order to evaluate model performance, we compute shifted geometric means (SGMs) of performance metrics, which provides a measure for relative differences. This avoids statistics being dominated by outliers with large absolute values as is the case for the arithmetic mean. The SGM also avoids an over representation of results with small absolute values. The SGM of values v 1 , ..., v M ≥ 0 with shift s ≥ 0 is defined as We say an instance is affected by a model, if solving this model finds a feasible solution; the instance is solved by this model, if solving this model finds an optimal solution. If an instance is unaffected, usually the model is too large to be read into the MILP solver.
We record the following performance metrics of each instance for each model, and compute the benchmarkwise SGMs: 1. t: the total running time in CPU seconds, with a shifted value set to 1 second; 2. σ: the relative dual gap, with a shifted value set to 0.01; 3. v r : the relative primal bound, with a shifted value set to 0.01.
For an unaffected instance, we set by default t = 1800, σ = 1 and v r = 1. Note that the time does not include the preprocessing time, since we find that the preprocessing is usually at most 0.5 seconds.
We will discuss the computational results, which are divided into two parts. In the first part, we compare the five models EF, F0, F, SF, and RF. We evaluate the performance metrics of these models. The second part compares RF and SFD, quantifying the facilities that are saved by allowing continuous location. In the following, we will analyze the aggregated results. We refer to Table 5 in the appendix for the detailed instance-wise results, where decimal representations of gaps and bounds are converted to percentages.

Comparative Analysis of Continuous Models
We compare five continuous models for the CSCP δ , namely EF, F0, F, SF, and RF. For each benchmark, radius and model, we record a triple of integers S/A/T: S denotes the number of solved instances, A denotes the number of affected instances, and T denotes the number of total instances in the benchmark. Moreover, we also report the average SGMs of the dual gaps, solving times and relative primal bounds among all instances in the benchmark. Table 3 summarizes these results.
First, we notice that EF cannot affect any instance in any benchmark; RF, SF and RF can affect all instances, i.e., solutions are provided by these models; RF is the model that solves the most number of instances (11), and SF is the second best one (10).
Secondly, we compare EF and F0. F0 is obviously superior to EF. With F0, 39 among 56 instances of small radius (resp. 42 among 56 instances of large radius) can be read by the CPLEX solver, while the instances modeled by EF are too large to read. Therefore, better solutions than trivial solutions are found by F0: on average, for instance of small radius (resp. large radius), F0 finds solutions that use 25.2% (resp. 74.1%) fewer facilities than the trivial solution.
Then, we compare F0, F and SF. With the delimitation of complete and partial covering sets, F and SF can affect all instances (especially those in Kgroup B, of which F0 could just read one). With the strengthening technique, SF has only marginal improvement in the relative primal bound, and solving time, while F is even slightly better than SF in the dual gap. We observe, in our experiment, that adding valid inequalities might slow down the internal solving process of CPLEX.
Finally, we compare SF and RF. RF outperforms SF in all performance metrics. Moreover, RF is the best one among those models affecting all instances. Indeed, for many instances, their degree-two-free networks may contain long edges, and RF avoids introducing too many variables and constraints for modeling their coverage.
Kgroup B is the hardest benchmark. The best model RF still has an average dual gap of 59.1% and 154.2% relative primal bound for instances of small radius, and this means that RF cannot produce better solutions than the trivial one.
We find that for all the models (except for EF), the average dual gaps and solving times of instances of large radius are smaller than those of instances of small radius. This shows that the large radius has a positive effect on the model performance, and an instance of the small radius may be more difficult than the same instance of large radius. This is because, with larger radius, the network after processing is smaller (see Table 5).
In Figure 7, we show scatter plots of the relative dual gaps and the relative primal bounds of affected instances between different settings. For every plot, there is a line in which the points have equal (X,Y)-values.
If points fall below the line, then the Y-axis model performs better for the corresponding instances. Note that when comparing F0 and F, the plots do not consider the unaffected instances of F0 which are affected or solved by F. Moreover, F0 even closes more duality gaps than F, but F can find better primal solutions. These plots give an overview of all affected instances and support the above analysis.
To summarize, we have shown that the two proposed techniques-that to delimit the coverage areas from a given point in Section 2, and that to cover long edges in Section 6,-can reduce the model size drastically.
Among the five models tested, RF features the best overall performance, which is achieved by directly modeling covers on long edges. On the other hand, delimiting the covering sets to the potential, complete and partial covers also reduces the model size, which allows F to read all the tested instances.

Comparative Analysis of Continuous and Discrete Models
In CSCP δ , the facilities are located either at nodes or edges, while in the discrete variant considered in this section facilities can only be located at nodes. Our objective is to evaluate the number of facilities that can be saved by allowing continuous location. Since the discrete model studied here, SFD, is a discrete restriction of CSCP δ , every optimal solution is a feasible solution of CSCP δ . We solve SFD for the discrete facility SCP δ and compare the results with the best model for CSCP δ , RF.
In addition to the previous performance statistics, we also record for each instance, a new relative primal bound for the continuous model defined as: where v d is the best solution found by SFD. If v r < 1, then the continuous model (in this case, RF) finds a solution better than the one found by the discrete model.      In Figure 8, we also show scatter plots of the relative dual gaps (σ) and the relative primal bounds (v r ) for those instances affected by both SFD and RF. These plots complement the averaged results of Table 4 by giving information on all affected instances, and support the above analysis. By allowing location at edges, the continuous model can reduce the number of installed facilities. However, it becomes more challenging to solve the problem. The results suggest that calling SFD and passing its solution as a warm-start to RF can make sense as a two-step optimization approach.

Conclusions
In this work, we use an integer programming approach to solve CSCP δ , propose various MILP formulations for this problem and test these formulations against an existing MILP formulation on several benchmarks from the literature.
The existing works mainly consider discretization methods and FDS. Discretization methods are indeed preprocessing procedures that restrict CSCP δ to an equivalent set-covering problem with continuous demands and candidate facilities in FDS. But FDS is only computable when CSCP δ satisfies some assumption, so it is not practical to employ FDS for general CSCP δ . On the other hand, to delimit the search space of MILP models, we use alternative preprocessing procedures. We explore the delimitation that relaxes the concept of FDS, which also restricts the candidate space of facilities (from the full network to a still continuous sub network). We learn the following ideas to tackle similar problems (possibly with more complex constraints): i).
Integer programming methods are more viable and flexible for general graphs, as a partial delimitation of the problems is useful for strengthening the models (via separation of valid inequalities, tightening big-M constants, variable fixing), while the discretization methods require a full delimitation and a complete characterization of FDS. ii). Our delimitation is applicable as long as the facility location is continuous.
Specifically, we devise and implement four models for CSCP δ : F0, F and SF, which belong to the same family of models, and RF, which is the reduced one. These models mainly differ by preprocessing procedures applied. We find that MILP solvers cannot read or build the MILP model from [16] for any instances in our test bed, and this is due to the model having a large number of constraints and variables. So we mainly compare our four models. We find the MILP size is the main barrier to scalability. The delimitation of those parts of the network that can be covered from a specific location has been revealed as a very effective technique to reduce the model size. In addition, avoiding breaking long edges in the reduced model also results in better scalability. In conclusion, RF is the best model: it can find good solutions with a small dual gap.
Meanwhile, the model SF is easily casted into SFD for the discrete restriction of CSCP δ . We find that allowing continuous facilities decreases the number of installed facilities but increases the solving time significantly. We note that SFD finds an optimal solution for the discrete facility SCP δ quickly, which is a primal solution for CSCP δ . Therefore, SFD can be called as a fast MILP-based primal heuristic for CSCP δ .
As for future studies, devising efficient heuristics to be integrated into MILP solvers can be useful to improve the primal performance of the proposed models. For instance, different relaxations of CSCP δ can be worth exploring, such as that where demand only happens at nodes (i.e. only nodes are to be covered).
Every solution of CSCP δ would be a solution of such relaxation, and hence the optimal value of the latter is a valid dual lower bound of the optimal value of CSCP δ . If solving this combinatorial relaxation is efficient and provides a stronger dual lower bound than the LP relaxation of CSCP δ , we can utilize this result and integrate the combinatorial dual bound into the MILP solver, which leads to a combinatorial branch-and-bound algorithm.
Moreover, we can investigate the potential use of FDS to further delimit the search space and integrate FDS in the preprocessing procedure.