Facility Location with Tree Topology and Radial Distance Constraints

Let Gd (V, Ed) be an input disk graph with a set of facility nodes V and a set of edges Ed connecting facilities in V. In this paper, we minimize the total connection cost distances between a set of customers and a subset of facility nodes S⊆V and among facilities in S, subject to the condition that nodes in S simultaneously form a spanning tree and an independent set according to graphs Gd and Gd, respectively, where Gd is the complement of Gd. Four compact polynomial formulations are proposed based on classical and set covering p-Median formulations. However, the tree to be formed with S is modelled with Miller–Tucker–Zemlin (MTZ) and path orienteering constraints. Example domains where the proposed models can be applied include complex wireless and wired network communications, warehouse facility location, electrical power systems, water supply networks, and transportation networks, to name a few. e proposed models are further strengthened with clique valid inequalities which can be obtained in polynomial time for disk graphs. Finally, we propose Kruskal-based heuristics and metaheuristics based on guided local search and simulated annealing strategies. Our numerical results indicate that only theMTZ constrained models allow obtaining optimal solutions for instances with up to 200 nodes and 1000 users. In particular, tight lower bounds are obtained with all linear relaxations, e.g., less than 6% for most of the instances compared to the optimal solutions. In general, theMTZ constrainedmodels outperform path orienteering ones. However, the proposed heuristics and metaheuristics allow obtaining near-optimal solutions in signicantly short CPU time and tight feasible solutions for large instances of the problem.


Introduction
Facility location has been a major research topic within last decades [1][2][3][4][5][6]. In general, a facility location problem is mainly concerned with the optimal placement of facilities in a prede ned geographical area in such a way that customers (users) can connect to facilities at minimum (distance) costs. Example application domains include wireless and wired network communications, warehouse facility location, electrical power systems, water supply networks, and transportation networks, to name a few [1]. Naturally, the optimal placement is also conditioned to several factors which directly depend on the application itself, e.g., avoiding placing hazardous materials near housing, fair distribution of hospitals, public schools, military bases, radar installations, branch banks, shopping centers, waste disposal facilities, police and/or re departments in a city to meet user demands and avoiding interference between base stations in a wireless communication network. Notice that all of these examples require facilities to be placed separately, which naturally imposes the condition of a radial distance constraint. In principle, any of them can be represented by means of a unit disk graph, and then an independent set structure appears as a natural way to model these constraints. In particular, any disk graph representation need not be a connected graph.
Let G d (V, E d ) and G c (V, E c ) be, respectively, the input disk and complete graphs composed of a set of facility nodes V drawn in the Euclidean plane with edge sets E d and E c , and let S be a subset of V. In this paper, we consider the problem of assigning a set K of users to facility nodes in S while simultaneously forming a tree backbone T � (S, E c (T)) with nodes in S where E c (T) denotes the set of edges in E c spanned by T such that |E c (T)| � |S| − 1 and with the additional condition that no two adjacent nodes in V related with G d can belong to S. e latter represents the radial distance requirement that we handle with independent set constraints [7]. Notice that finding the optimal subset S leads to a hard combinatorial optimization problem. For this purpose, we assume that all candidate sites for facility placement are represented by means of a unit disk graph with some predefined distance facility radial coverage. If two facility nodes are located near a certain distance, then at most, one of them can be selected for the tree backbone. Finally, users can connect to the resulting subset S. Notice that forming a tree backbone with facilities is an efficient strategy to ensure connectivity in a network [5,[8][9][10][11][12]. A unit disk graph is the resulting graph obtained by intersecting several unit disks in the Euclidean plane where each vertex corresponds to a center of each disk and with edges connecting them whenever the corresponding vertices lie within a constant (unit) distance of each other.
However, an independent or stable set of an arbitrary graph G � G(V, E) with set of nodes V and set of edges E is a subset of vertices S ⊆ V, no two of which are adjacent.
is means that, for every two vertices in S, there is no edge connecting them. Notice that a maximal independent set is a subset S ⊆ V in which by adding any other vertex forces set S to loose the independent set property. However, a maximum independent set is a maximal independent set with largest cardinality among all maximal independent sets. Consequently, every maximum independent set is maximal although the converse does not necessarily hold. Finding the maximum independent set is an NP-hard optimization problem, and the number of maximal independent sets in a general arbitrary graph with n nodes is of the order of O(3 n/3 ) [13]. Furthermore, notice that finding a maximal independent set in an arbitrary graph G is equivalent to finding a maximal clique in its complement graph G. e reader must not be confused with the fact that if it is possible to list all maximal cliques of G in polynomial time, then listing all maximal cliques is also possible to achieve in its complement graph in polynomial time. is is not true in general. In particular, for unit disk graphs, all maximal cliques can be listed in polynomial time [14,15]. However, this is not possible to achieve with its complement graph which is no longer a unit disk graph. In fact, finding the maximum independent set in a unit disk graph has proved to be an NP-complete problem even if a disk representation of the graph is given [16][17][18].
In general, listing all maximal cliques in polynomial time can be achieved for sparse graphs [15,19]. e latter can be performed, for example, by using the Bron-Kerbosch algorithm [14,19].
We propose four compact polynomial formulations based on classical and set covering p-Median models in order to minimize the total connection cost distance between facilities and between customers and facilities simultaneously [10,20,21]. However, the tree to be formed with S is modelled with Miller-Tucker-Zemlin (MTZ) and path orienteering constraints, respectively [22,23]. Notice that path orienteering constraints have been recently proposed in [23] to avoid cycles in trees. Subsequently, we further strengthen the MILP models by adding clique valid inequalities which we obtain in polynomial time using the Bron-Kerbosch algorithm [14,15]. Notice that the fact that we can list all maximal cliques in unit disk graphs allows us to strengthen our proposed models in a straightforward manner, leading to tight gaps when compared to the optimal solutions. Finally, we propose Kruskal-based heuristics and metaheuristics adapted from a greedy approach proposed for the maximum independent set problem [24,25]. e greedy algorithm mainly consists of selecting the vertex of minimum (or maximum) degree of G and removing it together with its neighbors from the graph, and it iterates on this process on the remaining graph until no vertex remains. Notice that the output of this greedy algorithm always results in a maximal independent set. In particular, the metaheuristics are constructed based on guided local search and simulated annealing strategies [26][27][28][29]. As far as we know, placing facilities at a certain radial distance using independent sets on unit disk graphs and with tree backbone representations have never been considered in the literature before. Similarly, the solution methods are new to the proposed models. e paper is organized as follows: In Section 2, we review some related works. en, in Section 3, we introduce the mathematical formulations together with their strengthened versions and present two examples of feasible solutions for the problem. In Section 4, we present our proposed heuristics and metaheuristics. Subsequently, in Section 5, we conduct substantial numerical experiments in order to compare all the proposed models and algorithms. Finally, in Section 6, we give the main conclusions of the paper and provide some insights into future research.

Related Work
Facility location dates back to the well-known Fermat point problem which consists of finding a point such that the total distance to three other points in an Euclidean plane is minimized. is problem was originally proposed by Pierre de Fermat and later solved by Evangelista Torricelli [4,6]. Today, its extended version is known as the Geometric Median Problem and it is considered as a continuous facility location problem [1,2]. Discrete facility location has also been considered as a major research topic since several decades ago [3,5], and in particular, the distance requirement between facilities or between facilities and demand points has been recognized as a key decision factor. For instance, in [5], the authors illustrate and highlight the distance requirement with several reallife applications.
eir examples include distribution systems where transshipment among stocking points are 2 Complexity of significant size, locating emergency facilities such as ambulance stations and fire fighting units with respect to population centers, gasoline stations, fast food restaurants, and spacing between wireless radio base stations, construction of nuclear power plants or missile silos nearby, and dump sites for chemical waste or collected refuse, to name a few. For a more general and deep review related with facility location topic, the reader is referred to [1][2][3][4][5][6][30][31][32] and references therein. In particular, Farahani et al. [3] present a recent survey related with distance covering location problems together with a comprehensive review of models, solution approaches, and applications where it also highlights the importance of the distance factor required when designing facility networks.
Hereafter, we briefly discuss some works where the authors propose models which are closer to ours. Perhaps, a first model that is worth mentioning is the covering tour problem [33].
is problem consists of determining a minimum length Hamiltonian cycle which has to be formed with a subset of nodes from a larger set of nodes V such that every vertex of an another set W, that must be covered, is within a specified distance from the cycle. is problem is similar in nature to our problem in the sense that it requires to form a backbone network topology with a subset of facilities. In order to solve this problem, the authors proposed an exact branch and cut algorithm and a heuristic. Similarly, covering location problems related with paths and trees are proposed in [34], where the facilities are assumed to be small in size in proportion to their location regions which is the case of subway metro stations and highways for example and where the aim is to find a path through the facility network of minimum length such that the population coverage is maximized. In [34], the authors develop Lagrangian solution approaches to solve their models. e maximal covering location problem (MCLP) is also worth mentioning where the objective is to maximize the amount of demand covered within an acceptable service distance by locating a fixed number of new facilities [35]. e problem is solved with linear programming relaxations and branch and bound techniques as well. Variants of MCLP are also presented and discussed in [3]. e indirect covering tree problem proposed in [36] represents another approach which is similar in nature to our facility location problem. ere, it is assumed that a given backbone spanning tree network is given a priori, and then, two optimization problems are derived, namely, the minimum cost covering subtree and the maximal indirect covering subtree problems. In the former, the aim is to find the minimum cost collection of arcs that form a subtree and satisfy covering constraints for nodes of the network, whereas in the latter, the subtree maximizes the demand within a given distance of nodes of the subtree. Reduction techniques that were initially proposed to solve the location set covering problem are extended by the authors to solve these new problems. More variants of these problems are also explained and discussed in [3,36].
In [37], the authors study two continuous facility location problems. e first one consists of placing a fixed number of unit disks in an area in order to maximize the total weight of the points inside the union of the disks, whilst the second one deals with placing a single disk with center in a given constant region in order to minimize the total weight of the points inside the disk. e authors propose approximation algorithms to obtain solutions for these problems. Similarly, the authors in [12] deal with a continuous facility location problem while including backbone connectivity together with their related costs. In particular, the objective here is to minimize the weighted sum of three costs, namely, the fixed costs required to install facilities, the backbone network costs required to connect facilities, and the transportation costs incurred from providing services from the facilities to the service region. e authors analyze the behaviour of their proposed model and derive two asymptotically optimal configurations of facilities. Finally, they give a fast constant factor approximation algorithm which allows us to find the placement of a set of facilities in a convex polygon that minimizes the total sum of the costs. We note that the authors also highlight the importance of considering connectivity among facilities in this reference. Finally, in the domain of telecommunications, some related works, although not so similar, can be consulted for instance in references [30][31][32].

Description of the Problem and Mathematical Formulations
In this section, we first describe the facility location problem we are dealing with. For this purpose, we present and explain two feasible solutions of the problem for different radial coverage values. Subsequently, we present our MILP formulations.

Description of the Problem.
Consider the undirected graphs G d � (V, E d ) and G c � (V, E c ) as defined in Section 1. Notice that E c contains the edges obtained with all pairs of nodes v 1 and v 2 in V, (v 1 ≠ v 2 ) as G c is complete. However, E d is composed of the edges obtained with all pairs of nodes v 1 and v 2 in V, (v 1 ≠ v 2 ) which are located at a radial distance lower than or equal to d. e facility location problem we are dealing with consists of finding a subset of facilities S ⊆ V such that S is an independent set of G d and a tree backbone of G c and where each user in the set K is assigned to its nearest facility in S, thus minimizing the total distance costs between facilities and users and among facilities themselves. Notice that S does not necessarily need to be a maximal independent set. In Figure 1, we plot two examples of input disk graphs with |V| � 30 nodes within an area of one square km using radial coverage values of 0.2 and 0.5 km, respectively. More precisely, on the left, we plot each input disk graph G d and the independent set obtained from G d . Nodes in S are colored green. However, on the right, we Complexity 3 plot the spanning tree backbone formed with nodes in S. Users are omitted for the sake of clarity. In particular, notice that, for an input graph G d with low density, the cardinality of S should be higher. On the opposite, if the density of G d is high, the cardinality of S should be smaller. By density, we refer to the ratio obtained by dividing the number of edges of E d over the total number of edges of E c .

MILP Formulations.
In order to write a first compact polynomial formulation for the facility location problem described in Section 3.1, we first define set A c � (i, j), (j, i)/ i, j ∈ E c as the set of directed arcs obtained from edges in E c and H � (V, A c ) as the digraph obtained from G c . Next, we define the binary variables: us, a MILP model can be verified by means of the following proposition.  Figure 1: Input graphs and feasible solutions obtained for an instance with |V| � 30 nodes and |K| � 50 users using radial coverage values of 0.2 and 0.5 km: (a) independent set using 0.2 km; (b) minimum spanning tree backbone using 0.2 km; (c) independent set using 0.5 km; (d) minimum spanning tree backbone using 0.5 km. 4 Complexity facility location problem while satisfying the following conditions: (i) All users in K are assigned to facilities in S.
(ii) Subset S ⊆ V forms an independent set (not necessarily maximal) from G d . i∈V x kj , y j , Proof. To prove (i) and (ii), we focus on explaining how the set of constraints in P 1 together with its objective function imply each of the assertions. For this purpose, we define the input matrices C � (C kj ) and D � (D ij ), ∀k ∈ K, i, j ∈ V (i ≠ j), to be the distances between user k and facility j and between facilities i and j, respectively. In particular, we require matrix D � (D ij ) to be symmetric. us, the total distance is minimized in (2). It is easy to see from the objective function (2) that the first double sum will be less expensive when the cardinality of S, which can be computed by |S| � j∈V y j , is larger. However, the second term will be cheaper when the cardinality of S is smaller. Consequently, we seek for a tradeoff between these two terms leading to an optimal number of facilities in S which can be obtained by solving P 1 . Notice that S is composed of the nodes j ∈ V for which the binary variable y j � 1. Constraints (3) ensure that each user is connected to a unique facility, whilst constraints (4) ensure that a user is connected to a facility node that belongs to S. Next, constraints (5) are independent set constraints and ensure that no two adjacent facilities in V can belong to S according to E d . Furthermore, note that S is not conditioned to be a maximal independent set.
To prove condition (iii), we show that any feasible solution (x, y, z, u) of P 1 is an acyclic connected arborescence of H � (V, A c ). First, notice that the constraint (6) imposes the condition that any tree formed with nodes in S contains exactly |S| − 1 arcs. However, the constraints (10) ensure that the number of arcs entering node j ∈ V is at most one if and only if y j � 1. Similarly, the constraints (11) ensure that the number of incoming and outgoing arcs in node j ∈ V can be at most |V| − 1 if and only if y j � 1. Otherwise, if y j � 0, then the binary variables z ij � z ji � 0, ∀i, j ∈ V (i ≠ j). Finally, the constraints (7)-(9) together with the constraints (6) and (10) avoid cycles. To see how these constraints avoid cycles, let us assume that the subset en, the constraints (7) imply that Notice that, in general, the contradiction holds for any subset of nodes of S which implies that the digraph induced by z � (z ij ) cannot contain directed cycles. Furthermore, notice that constraints (7) do not consider the particular case when an arc is in the opposite direction in the sequence (v 1 , v 2 ), (v 2 , v 3 ), . . . , (v j , v 1 ) for some j ∈ S. is particular case is avoided by means of the constraints (6) and (10), which finally ensure that the digraph is connected with |S| − 1 directed arcs while avoiding the fact that a particular node has more than one incoming arc. Consequently, by dropping the directions of each arc, we obtain an undirected tree topology. Recall that the cost distance matrix D � (D ij ) ∀i, j ∈ V, (i ≠ j), is symmetric. Finally, the constraints (12) limit the domain of the decision variables, thus concluding the proof. □ Notice that P 1 is constructed based on a classical formulation of the p-Median problem and a directed Miller-Tucker-Zemlin characterization of the spanning tree polytope [10,21,22].
Regarding the complexity of the facility location problem associated with P 1 , we outline the following theorem.

Theorem 1.
e facility location problem associated with P 1 is NP-hard.
Proof. To see this, let us consider the following weighted version of the objective function (2) as where 0 ≤ α ≤ 2. In particular, when α � 0, any spanning tree will be feasible for P 1 at zero cost for the objective function of the problem. Consequently, in this case, P 1 can be equivalently written as the following optimization problem: x kj ≤ y j , ∀k ∈ K, j ∈ V, where the second quadratic sum in (14) is due to an equivalent form of writing independent set constraints [38]. For this purpose, parameter M should be a nonnegative real value such that M > 2 k∈K j∈V C kj . In general, notice that any independent set S implies i,j { }∈Ed y i y j � 0. Consequently, it is easy to see that P α 1 is a quadratic version of the classical uncapacitated facility location problem which is NP-hard [39].
Subsequently, when 0 ≤ α < 2, we aim to solve P α 1 with the implicit requirement that the number of facilities must be between 1 ≤ j∈V y j ≤ β(G d ), where β(G d ) denotes the cardinality of the maximum independent set of G d . Consequently, in this case, the problem is equivalent to finding a minimum spanning tree with variable number of nodes in an arbitrary graph which is also known to be an NP-hard optimization problem by reduction from the Steiner tree problem [40]. Finally, when α � 2, the objective function (13) equals zero and the solution can be found trivially in linear time. In this case, the cardinality of the optimal solution set S * will be equal to one and all users will be assigned to this unique facility node. Obviously, the tree obtained in this case does not contain any edge, and hence, its objective value equals zero as well.
□ Notice that, in particular, the objective function of P 1 is obtained with α � 1 which is set by default. e reason behind this is that we minimize the total cost distances with fair degree of importance for both users and facilities. However, below in the numerical result section, we further consider the more general case for values of α ∈ [0; 2] which reflects different cost structures.
An alternative formulation of P 1 can be obtained by removing constraints (7) and (9) together with the decision variables u j for each j ∈ V and by introducing the nonnegative precedence variables p ij for all i, j ∈ V (i ≠ j). Precedence variables are used to indicate if there exists a direct path from i to j. If such a path exists, then p ij � 1; otherwise, p ij � 0. Notice that precedence variables have recently been introduced in order to deal with subtour elimination constraints in combinatorial optimization problems defined under a tree topology structure [23]. us, an equivalent MILP model can be written as x kj ≤ y j , ∀k ∈ K, j ∈ V, i∈V x kj , y j , z ij ∈ 0, 1 { }, ∀k ∈ K, i, j ∈ V, Constraints (16) state that if there is a path from i to j, for all i, j ∈ V (i ≠ j) that belongs to the solution, then a direct path going from i to j must exist. Similarly, constraints (17) ensure that there should be at most one path either from i to j or from j to i for all i, j ∈ V, (i ≠ j). Finally, constraints (18) and (19) ensure that if there is a direct path from i to j and the arc path j to k belongs to the solution, then both path variables P ij and P ik must be equal.
Notice that both P 1 and P 2 are formulated based on a classical formulation of the p-Median problem while using MTZ and path orienteering constraints, respectively. Next, we present two formulations which are based on a set cover formulation of the p-Median problem [20,41] together with MTZ and path orienteering constraints as well. We denote these two models hereafter by P 3 and P 4 , respectively. In order to write a first model under this approach, we first construct for each by sorting in the ascending order 6 Complexity the different entries of each k-th row of the original cost matrix C � (C kj ) in (2) in such a way that C o k1 � 0 and C o k|V| � max C k,j , j ∈ V . Since the first entry in each row of these vectors equals zero, we need to introduce an artificial facility node labeled as "1" hereafter. Also, we need to expand the input graphs G d and G c with this additional node. us, we expand sets V, E d , and E c and denote them by , by adding to D an additional row and a column in position one with zero entries. Consequently, a first set covering-based formulation using MTZ constraints can be written as In P 3 , notice that the binary variable x kj for each k ∈ K, j ∈ V o , acts as a cumulative variable where x kj � 1 if and only if the allocation cost distance of user k is at least C o kj and x kj � 0 otherwise. For this purpose, the set covering constraints (24) ensure that variable x kj � 1 if there is no open facility at less than distance C o kj . e positive coefficients of the first double sum in the objective function (23) ensure that variable x kj � 0 otherwise. Next, the constraint (25) ensures that node "1" is excluded from the solution of the problem. Finally, all remaining constraints have the same effect as for P 1 . Analogously as for P 2 , we propose the following MILP formulation using path orienteering constraints: Ultimately, we derive strengthened versions for each of the proposed models and denote them hereafter by P s 1 , P s 2 , P s 3 , and P s 4 , respectively. In particular, from the MTZ constrained models, we remove constraints (5), (7), and (11) and replace them by the following sets of constraints: where the constraints (27) avoid cliques in the solution set S. For this purpose, we define the set Q containing all maximal cliques of G d . Recall that, in practice, all cliques can be obtained in polynomial time as G d is a unit disk graph [14,15]. Consequently, each row in the binary matrix MC(q, ·) corresponds to a maximal clique q ∈ Q. Constraints (28) are valid cuts obtained from constraints (7) and adapted for the spanning tree polytope [10,22]. Finally, constraints (29) and (30) are valid cuts for the constraints (11) [11].

Algorithms
In this section, we propose Kruskal-based heuristics and metaheuristics which are adapted from a greedy approach initially proposed for the maximum independent set problem [25]. In particular, the metaheuristics are constructed based on guided local search and simulated annealing greedy strategies [26][27][28][29].

Heuristic Approach.
ree variants of the same heuristic are proposed. e only difference between them relies on the greedy decision of adding new nodes to S at each iteration.
us, we only present the first one referred to as GMIN in Algorithm 1 and explain the other ones based on these steps. e heuristic requires the input matrices C � (C kj ), en, for each vertex v ∈ V, the following steps are performed. e set S is initialized as an empty set, and W is assigned set V. e set W is used as an auxiliary set. Next, v is added to S, and it is removed from W together with its neighbors. Subsequently, we obtain the minimum spanning tree with nodes in S using the Kruskal algorithm and assign each user k ∈ K to its nearest facility in S. Finally, we compute the objective function value of P 1 and save the current solution if its objective value is less than the best found so far. Initially, the objective value is set to infinity. en, in Step 2, we enter into a while loop and iterate until G d � (W, E d (W)) is empty, where E d (W) denotes the set of edges in E d induced by remaining nodes in W. Inside the while loop, we proceed similarly by selecting a new vertex v ∈ W with minimum or maximum degree to be added to S or by interchanging between minimum and maximum degrees.
is is the greedy decision which differentiates remaining variants of Algorithm 1. Hereafter, we denote the three heuristics by GMIN, GMAX, and IC, respectively. Finally, Algorithm 1 returns the best solution found and its objective function value.

Guided Local Search Approach.
e first metaheuristic we propose is based on a guided local search strategy and requires the same input parameters of Algorithm 1 [29]. Its pseudocode is presented in Algorithm 2, and it is explained as follows. In Step 1, first we set W � V and S � ∅, and iterate while G d � (W, E d (W)) is not empty. e following steps are performed inside the loop. First, we remove the node with minimum degree v from W and add it to S. en, we remove all neighbors of v from W and find the minimum spanning tree with current nodes in S using the Kruskal algorithm. Finally, we assign each user k ∈ K to its nearest facility in S, compute the objective function value of P 1 , and save S * � S if the objective value is less than the best found so far. Similarly, in Step 2, we enter into a second, while loop and iterate until the CPU time is greater than or equal to maxCPU which is an input parameter. Inside this loop, we penalize each value in vector Deg(S * ) corresponding to each node in S * by adding a random number uniformly distributed in (0, ρ). Initially, Deg(v) contains the degree of each node v ∈ V. Notice that this step acts as a guided local search strategy, and its main purpose is to avoid repeated nodes in the new solutions generated by the algorithm. en, we reset S � ∅ and W � V and construct a new independent set S from G d � (W, E d (W)) with cardinality less than or equal to R. Observe that, in this case, S is constructed without evaluating the objective function each time a new vertex is added to S. Once set S is constructed, we proceed to find the minimum spanning tree formed with nodes in S using the Kruskal algorithm, assign each user k ∈ K to its nearest facility in S, and compute the objective function value of P 1 . If this value is less than the best found so far, we update S * � S and reset the current CPU time to zero. Finally, we randomly generate a value c ∈ − 1, 0, 1 { } in order to decrease, maintain, or increase, respectively, the cardinality of the next independent set to be generated by the algorithm starting from |S * |. Notice that Algorithm 2 updates degree values in Deg(·) according to new and better solutions obtained following a simple local search rule in contrast to classical guided local search metaheuristics which use augmented objective functions in order to perform the local search [29].

Simulated Annealing
Approach. Next, we further consider another metaheuristic approach based on a simulated annealing (SA) greedy strategy. SA is a classical metaheuristic which was proved to be highly efficient when finding near-optimal solutions for hard combinatorial optimization problems [26][27][28]. SA randomly searches for ascent moves in order to escape from local minima and to find global optimal solutions. It is basically inspired on the annealing process of a material in metallurgy and consists of heating and cooling steps which must be controlled in such a way that certain equilibrium is reached in terms of temperature. In particular, the cooling step is analog to a slow decrease in the probability of accepting worse solutions in order to allow a major degree of diversity while performing the search process. e underlying idea in the SA algorithm is simple and can be described as follows. First, SA requires an initial feasible solution to start. Next, it randomly generates a neighbor solution. If this neighbor solution is better than the best found so far, it is accepted as a new best solution. Otherwise, SA allows us to move to a worse solution with a certain probability. In fact, this is the key ingredient that allows us to escape from local minima. In general, the probability of moving to a new solution is determined by Boltzmann distribution while varying the temperature which is gradually decreased. Initially, it is high enough to allow a high degree of diversity since probabilities are close to one. However, small temperature values allow probabilities to go down to zero which is equivalent to a pure local search method. e SA procedure we propose is depicted in Algorithm 3, and it is explained as follows. It requires the same parameters used by Algorithms 1 and 2. Initially, the 8 Complexity temperature T 0 is set to a large positive value T i . Inside the while loop of Step 2, we set S � ∅ and W � V, and for each v ∈ S * , which is obtained in step one, we pick randomly a number from [0, 1] and compare it with the input parameter η which is arbitrarily chosen from the interval [0, 1] as well.
If the random number is less than η, then we include node v ∈ S * in the new independent set S and remove from W node v and its neighbors. Otherwise, we skip node v and continue. is allows us to keep an η percentage of nodes in S * . Next, we start a while loop in order to add new nodes in Result: a feasible solution for P 1 .
Step 1:; Find the minimum spanning tree formed with nodes in S using the Kruskal algorithm; Assign each user k ∈ K to its nearest facility in S and compute the objective function value of P 1 ; Save the current solution found if its objective function value is less than the best found so far; Step 2:; Find the minimum spanning tree formed with nodes in S using the Kruskal algorithm; Assign each user k ∈ K to its nearest facility in S and compute the objective function value of P 1 ; Save the current solution found if its objective function value is less than the best found so far; Return the best solution found and its objective function value; ALGORITHM 1: Proposed heuristics.
Result: a feasible solution for P 1 .
Step 1:; Find the minimum spanning tree formed with nodes in S using the Kruskal algorithm; Assign each user k ∈ K to its nearest facility in S and compute the objective function value of P 1 ; Save the current solution found in S * if its objective function value is less than the best found so far; R � |S * |; Step 2:; while (cpuTime ≤ maxCPU) do Update degree vector Deg(·) according to the indices of each element in S * as Deg(S * ) � Deg(S * ) + (0, ρ); Find the minimum spanning tree formed with nodes in S using the Kruskal algorithm; Assign each user k ∈ K to its nearest facility in S and compute the objective function value of P 1 ; Save the current solution found in S * if its objective function value is less than the best found so far. If a better solution is found, set R � |S * | + c; Return best solution found and its objective function value; ALGORITHM 2: Guided local search-based metaheuristics.
Added nodes are chosen randomly from remaining nodes in W. Subsequently, we find the minimum spanning tree with nodes in S, assign each user k ∈ K to its nearest facility in S, and compute the objective function value of P 1 . If a better solution is found, we save this new solution in S * * and set S * � S * * and the current CPU time to zero in order to extend the duration of the algorithm. Otherwise, we compute the variable Δ as the difference between the current objective function value minus the previous one, generate a random number r from [0, 1], and compute the probability p � exp(− Δ/T 0 ). If r is less than p, we accept the new solution and save it in S * ; else, we set S * � S * * . Notice that the value of Δ will be positive when the objective value of the current solution is greater than the previous one, i.e., when the new solution is worse. In particular, when the initial temperature T 0 is large enough and the difference in Δ is low, the probability p will be close to one. On the opposite, when temperature values decrease and Δ is positive, the probability p will be close to zero and the algorithm will behave as a pure random local search algorithm. e next lines inside the while loop of Step 2 ensures that T 0 decreases at a δ rate. en, we check if T 0 < ϵ where ϵ is a small positive value. If so, we reset T 0 � T i and allow the algorithm to interchange between exploration and exploitation more frequently during the process. Finally, we randomly generate a value c ∈ − 1, 0, 1 { } to decrease, maintain, or increase, respectively, the size of the next independent set to be generated by the algorithm. Notice that the value of T 0 is of critical importance to handle efficiently the algorithm. If T 0 decreases rapidly, then a premature convergence to a local minimum may occur. On the opposite, if T 0 decreases slowly, then the algorithm will converge slowly as well. In order to guarantee that the algorithm converges to global optimal solutions with probability one, T 0 should be decreased at the logarithmic rate [42]. However, this can lead to poor convergence rates, so in practical settings, one is usually forced to decrease T 0 (t) at iteration t by T 0 (t) � δ * T 0 (t − 1) where 0.85 ≤ δ ≤ 0.96 [43].

Numerical Experiments
In this section, we perform substantial numerical experiments in order to compare all the proposed models and algorithms. For this purpose, we randomly generate connected Euclidean disk graph instances. More precisely, nodes are uniformly distributed in a plane within an area of one square kilometer, and each pair of nodes is connected if the distance of the edge connecting them is less than or equal to the predefined radial coverage. However, the complete version is obtained by connecting all pairs of nodes. Consequently, each entry in the input matrices C � (C kj ) and D � (D ij ) for all k ∈ K, i, j ∈ V, represents the distance between user k and facility j and distance between facilities i and j, (i ≠ j), respectively. Matrix D is symmetric. We generate two sets of instances with dimensions of |K| � 500, 600, 700, 800, 1000  [44]. In order to solve the strengthened versions of each MILP model, we find all cliques for each disk graph using the Bron-Kerbosch algorithm [14,19].
We implement a Matlab program using CPLEX 12.8 to solve all the proposed models P 1 , P s 1 , LP 1 , LP s 1 , P 2 , P s 2 , LP 2 , LP s 2 , P 3 , P s 3 , LP 3 , LP s 3 , P 4 , P s 4 , LP 4 , and LP s 4 [44]. e proposed algorithms GMIN, GMAX, IC, GL, and SA are also implemented in Matlab, where GL refers to Algorithm 2 and SA refers to Algorithm 3, respectively. e parameters in Algorithms 2 and 3 are calibrated to maxCPU � 100 s, η � 0.5, δ � 0.85, ϵ � 10 − 4 , T i � 1000, and c ∈ − 1, 0, 1 { }. e numerical experiments have been carried out on an Intel(R) 64 bits core (TM) with 3.40 GHz and 8 gigabytes of RAM. Tables 1-4, we present numerical results for the linear models. In these tables, we solve the first two sets of instances which use both radial coverage values. us, each row in each of these tables corresponds to a particular instance. Notice that, in Tables 2-4, we do not report numerical results for some of the instances because the linear models cannot be solved with CPLEX due to shortage of memory [44]. e column information in these four tables is exactly the same with the exception of Table 1 which includes four additional columns with instance dimensions. e information of each column is described as follows. In column 1, we present the instance number. en, only in Table 1, columns 2-5 report the number of users, the number of facility nodes, the density of each input graph, and the number of optimal facilities obtained with P 1 . e density of each graph is computed by dividing the number of edges of the input disk graph over the total number of edges. Next, in columns 6-10 of Table 1, which correspond to columns 2-6 in Tables 2-4, we present the optimal solution or best solution found with CPLEX in two hours for the MILP models, number of branch and bound nodes, CPU time in seconds, optimal solutions for the LP relaxations, and CPU time in seconds, respectively. Similarly, in columns 11-15 of Table 1 which correspond to  columns 7-11 in Tables 2-4, respectively, we present the same column information for the strengthened versions of the MILP models. Finally, the last two columns in Tables 2-4 report gaps that we compute by ((Opt − Opt LP )/Opt) * 100 where Opt refers to the optimal solution of each MILP model and Opt LP represents the optimal solution obtained with its linear relaxation.

Numerical Results for the Linear Models. In
From Table 1, we observe that the density of the input disk graphs increases with the radial coverage. On the opposite, we see that the number of optimal facilities obtained with P 1 is larger for the sparse graphs, which is obvious since sparse graphs contain more independent sets than dense graphs. In particular, we observe that, for the instances #16-#20, CPLEX cannot solve the MILP models to optimality within two hours of CPU time using a radial coverage of 0.2 km, which is also the case for the instances #17-#20 using a radial coverage of 0.5 km. All remaining instances are solved to optimality. Next, we observe that sparse graphs are harder to solve to optimality which is reflected in terms of branch and bound nodes for example. In particular, these values and the CPU times are lower for the strengthened models. Regarding the optimal solutions obtained with the LP relaxations, we see that these bounds are significantly tighter for the strengthened models as well.
is can be confirmed by the gap columns which show that the linear relaxations are less than 6% for most of the instances when compared to the optimal solutions. However, the CPU times for the strengthened linear relaxations are considerably higher.
In Table 2, first we observe that the optimal solutions obtained with P 2 and P s 2 are exactly the same as those obtained with P 1 and P s 1 in Table 1. In general, we observe that path orienteering constraints decrease CPLEX performance significantly. Consequently, we could only solve instances #1-#10 in this case. is observation is also valid for the LP relaxations which require more CPU times than for solving LP 1 and LP s 1 , respectively. Regarding the number of branch and bound nodes, we observe a slightly less number of nodes required to solve the MILP models in Table 2. Finally, we observe that the gaps are exactly the same in both Tables 1 and 2.
In Tables 3 and 4, we report numerical results obtained with P 3 , P s 3 , P 4 , and P s 4 , respectively. As it can be observed, in these tables, we cannot solve more than 15 and 10 instances with CPLEX in each table, respectively. From the obtained results, we observe that the optimal solutions are the same as in Tables 1 and 2. Similarly, we obtain the same gap values for the LP bounds when compared to the optimal solutions. In general, we observe that the strengthened models allow us to obtain optimal solutions more rapidly and verify that path orienteering constraints deteriorate CPLEX performance. Data: input matrices C � (C kj ) and D � (D ij ) and input graph G d � (V, E d ).
Result: a feasible solution for P 1 .
Step 1:; Set S � ∅ and W � V; while loop as in Step 1 of Algorithm 2; R � |S * | and T 0 � T i ; Step 2:; Find the minimum spanning tree formed with nodes in S using the Kruskal algorithm; Assign each user k ∈ K to its nearest facility in S and compute the objective function value of P 1 ;

if (A better solution is found) then
Save it in S * * and set cpuTime � 0 and S * � S * * ; else Compute Δ as the difference between the current objective value minus the previous one; Draw a random number r from [0, 1]; Accept the new solution and save it in S * ; else R � |S * | + c; Return best solution found and its objective function value; ALGORITHM 3: Simulated annealing-based metaheuristic.

Complexity 11
Finally, from the numerical results presented in Tables 1-4, we conclude that model P s 1 outperforms the set covering formulations. In general, we observe that MTZ constraints seem to work better than path orienteering ones.
In Tables 5 and 6, for each input disk graph instance in Table 1, we present the total number of maximal cliques, maximum clique number, total number of maximal independent sets, and maximum independent set number, respectively. We compute this information by using the Bron-Kerbosch algorithm [14,15]. Recall that the total number of cliques can be obtained in polynomial time for disk graphs. Consequently, the maximum clique number can be obtained straightforwardly by selecting the maximal clique with largest cardinality. However, the total number of maximal independent sets can be obtained with the Bron-Kerbosch algorithm listing all maximal cliques from its complement graph Notice that G d is no longer a disk graph, and therefore, we cannot expect to find all cliques in polynomial time in this case. In fact, the total number of maximal independent sets in an arbitrary graph with n nodes is of the order of O(3 n/3 ) [13]. Moreover, finding the maximum independent set is an NP-hard optimization problem. Although it is a difficult task to list all maximal independent sets, we still run the Bron-Kerbosch algorithm for two hours in order to give some insights with respect to the difficulty of doing this for graph instances with different radial coverage values. As it can be observed from Table 5, while using a radial coverage of 0.2 km, we can find all maximal cliques. On the opposite, we cannot list all maximal independent sets for any of the   instances in two hours. Consequently, in this case, we obtained the maximum independent set number by solving the following optimization problem:

Complexity
Similarly, from Table 6, while using a radial coverage of 0.5 km, we can find all maximal cliques and all maximal independent sets for half of the instances (e.g., instances #1-#10). is clearly shows that the sparser the input disk graphs are, the harder it is to list its maximal independent sets. We also see from Tables 5 and 6 that the maximum independent set numbers are slightly larger than the optimal number of facilities obtained with P 1 in Table 1. Notice that the maximum independent set number is an upper bound for the optimal number of facilities since no larger maximal set can be obtained from G d . is observation can be understood by simply looking at the objective function of the proposed P 1 model which seeks for a tradeoff between the two terms. In particular, when the number of facilities is close to the maximum independent set number, it means that the first double sum in the objective function of P 1 dominates the second term, which also means that the cost of assigning the total number of users to the facilities is significantly larger than the cost of the spanning tree required to connect facilities. Notice that the trade off between the two terms can vary with parameter α ∈ [0, 2] used in eorem 1. Numerical results while varying parameter α are discussed and reported below. Tables 7-9, we report numerical results for the proposed heuristics and metaheuristics for all the instances reported in Table 1 and for the larger sets of instances mentioned above. More precisely, numerical results for the instances in Table 1 are  reported in Tables 7 and 8, while Table 9 reports numerical results for the larger sets. In Table 7, for the sake of clarity, we repeat some information of Table 1. us, in columns 1-5, we present the instance number, number of users to be assigned to the facilities, number of nodes of the input graph, and optimal solution or minimum  - 16  14 Complexity Table  6: Maximal cliques and maximal independent sets obtained with the Bron-Kerbosch algorithm for the instances #1-#20 in Table 1 using a radial coverage value of 0.5 km. Instance number #     18 Complexity solution found with P 1 or P s 1 together with its minimum CPU time in seconds, respectively. en, in columns 6-9, 10-13 and 14-17, we report number of facilities, upper bounds, CPU times in seconds, and gaps for each heuristic GMIN, GMAX, and IC, respectively. e gaps are computed by ((UB − Opt)/Opt) * 100.

Numerical Results for the Algorithms. In
From Table 7, first we observe that the number of facilities is nearly the same as the optimal one for GMIN heuristics. In general, GMIN allows us to obtain a higher number of facilities than both IC and GMAX, whereas IC obtains more facilities than GMAX. We note that this can be explained as a consequence of the greedy decision that each heuristic performs at each iteration. Notice that GMIN removes less nodes from the graph than IC, which in turn removes less nodes than GMAX. Next, we further observe that GMIN obtains tighter bounds than the other heuristics although at a higher CPU time effort. is is also consequence of removing less nodes at each iteration since more nodes remain in subsequent iterations. Subsequently, by looking at the gap columns, we see that GMIN dominates the other ones with near-optimal gaps. In particular, the gaps are tighter for the instances that use a higher radial coverage value. Finally, we observe that the solutions obtained by GMIN require significantly less CPU time than those obtained with P 1 or P s 1 . Notice the fact that the first term in the objective function of P 1 dominating the second one implies that more facilities should be open, which should not always be the case if the cost to connect facilities under the spanning tree requirement is higher than the cost of assigning users to facilities. In order to deal with other situations where the cost of each term in the objective function of P 1 is different, later we further present numerical while varying the parameter α ∈ [0, 2] using the weighted objective function presented in eorem 1.
In Table 8, for the sake of clarity, we repeat columns 1-5 from Table 7, whereas in columns 6-9 and 10-13, we present numerical results for GL and SA, respectively. More precisely, in these columns, we report number of facilities, upper bounds, and CPU times in seconds and gaps that we compute by ((UB − Opt)/(Opt)) * 100.
From Table 8, we observe that the number of facilities is slightly larger for SA than for GL for most of the instances while using a radial coverage value of 0.2 km. However, for the coverage value of 0.5 km, most of these numbers remain nearly the same. Next, we see that the upper bounds obtained with SA are tighter than those obtained with GL when compared to the optimal solutions using both radial coverage values. In particular, this improvement is higher for 0.2 km which corresponds to the subset composed of the harder instances. e gap columns confirm these improvements. Finally, we observe that the CPU times obtained with GL and SA are significantly smaller than those required by CPLEX to solve the MILP models. In particular, we obtain negative gaps for the instances #18 and #20 using 0.2 km and for the instances #19 and #20 using 0.5 km which evidences the difficulty of finding optimal solutions for instances with higher dimensions.
In Table 9, we report similar results as those presented in Table 8, but for the larger set of instances. In columns 1-4, we report the instance number, number of users, number of facility nodes, and density of each input graph, respectively. en, in columns 5-7, we report number of facility nodes, upper bounds, and CPU times in seconds for GMIN heuristics. Similarly, in columns 8-11 and in columns 12-15, we report number of facilities, upper bounds, CPU times in seconds, and gaps for GL and SA, respectively. Notice that, in this case, we cannot solve the MILP models with CPLEX, and then we compute the gaps using as a reference the upper bounds obtained with GMIN. More precisely, we compute the gaps as ((Ub(GMIN) − Ub(GL))/Ub(GL)) * 100 and ((Ub(GMIN) − Ub(SA))/Ub(SA)) * 100 in columns 11 and 15, respectively.
From Table 9, we mainly observe that the number of facilities obtained with the three algorithms remain nearly the same. Next, we observe that the upper bounds obtained with SA are tighter than those obtained with GL, which in turn are tighter than those obtained with GMIN for most of the instances when using a radial coverage of 0.2 km. On the opposite, GL obtains tighter gaps than SA and GMIN when using a radial coverage of 0.5 km. In particular, a few negative gaps are reported when GMIN allows us to obtain better solutions compared to a metaheuristic. is is the case for the instances #6, #9, and #12, and for the instances #9 and #11 using 0.2 and 0.5 km, respectively. Finally, we observe that higher CPU times are required for GL compared to SA for dense instances. However, the opposite occurs for the sparse ones. Notice that the CPU times obtained with GMIN are deterministic which is not the case for the metaheuristics. In order to give more insights with respect to the upper bounds and CPU times obtained, later, we report some average results as well.
In Tables 10 and 11, we further report total number of maximal cliques, maximum clique cardinalities, number of maximal independent sets, and maximum independent set numbers for the large instances presented in Table 9 for the radial coverage values of 0.2 and 0.5 km, respectively.
Similarly as for Tables 5 and 6, the total number of maximal cliques and maximal independent sets are obtained with the Bron-Kerbosch algorithm when it is possible. Otherwise, we omit presenting this information. In particular, when it is not possible to obtain all maximal cliques of G d , its maximum clique number is obtained by solving the optimization problem (MIS) presented above using the set E d which corresponds to the set of edges of the complement graph of G d .
From Tables 10 and 11, we observe similar trends as in Tables 5 and 6, respectively. We observe that it is not possible to list all maximal cliques with the Bron-Kerbosch algorithm in two hours for some of the instances although it is known they can be obtained in polynomial time [14,15]. Finally, we see that the maximum independent sets are slightly larger than the number of facilities found with the proposed algorithms.
In order to give more insights with respect to the behaviours of P s 1 , LP s 1 , GMIN, GMAX, IC, GL, and SA algorithms while varying α using the objective function of eorem 1, in Figures 2-7, we plot upper bounds, optimal solutions, and number of facilities and CPU times in seconds Complexity Table 10: Maximal cliques and maximal independent sets obtained with the Bron-Kerbosch algorithm for the instances #1-#15 in Table 9 using a radial coverage value of 0.2 km. 1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  # Maximal C.  2184 2295 2130 2273 2073 4966 4798 4923 4743 4863 10974 12250 12843 13131 12086  Maximum C.  30  31  30  32  29  42  37  46  41  39  52  50  56  54  48  # Maximal I.S.  ---------------Maximum I.S.  28  28  28  28  27  29  28  28  29  29  30  30  29  30  30 -, not solved in 2 hours. Table 11: Maximal cliques and maximal independent sets obtained with the Bron-Kerbosch algorithm for the instances #1-#15 in Table 9 using a radial coverage value of 0.5 km.     Table 1 while varying α ∈ [0; 2] and using a radial coverage value of 0.2 km.   Table 1 while varying α ∈ [0; 2] and using a radial coverage value of 0.5 km.   Table 1 while varying α ∈ [1.7; 2] and using a radial coverage value of 0.2 km.    Table 1 while varying α ∈ [1.7; 2] and using a radial coverage value of 0.5 km.     Table 9 Table 9 while varying α ∈ [0; 2] and using a radial coverage value of 0.5 km. 24 Complexity for the instances #11 and #10 in Tables 1 and 9, respectively, while varying α ∈ [0; 2] and using radial coverage values of 0.2 and 0.5 km. In particular, in Figures 4 and 5, we restrict the values of α in the interval [1.7; 2] in order to decrease the number of optimal facilities while searching for a tradeoff between the two terms in the objective function. From Figures 2 and 3, we observe that the upper bounds and optimal solutions decrease with α. In particular, in Figure 2, the LP bounds remain near optimal for all values of α. Furthermore, we see that the upper bounds decrease in the following order: Ub(GMAX) > Ub(IC) > Ub(GMIN) > Ub(GL) > Ub(SA). Regarding the number of facilities (or medians), we observe that the order of the methods is as follows: Med(GMIN) > Med(SA) > Med(GL) > Med (IC) > Med(GMAX), where the optimal number of facilities is closer to Med(SA). However, the order of the methods with respect to the CPU times is as follows: CPU(P s 1 ) > CPU(GL) > CPU(SA). For the heuristics, the CPU times remain nearly the same.

Instance number
From Figure 3, we observe similar trends although, in this case, the curves are closer to each other with the exception in the CPU times which are significantly higher for    From Figure 4, we observe the following ordering for the upper bounds: Ub(GMIN) > Ub(IC) > Ub(GMAX) > Ub(GL) > Ub(SA). As it can be observed, in this case, GMAX allows us to obtain better solutions than GMIN. is can be explained by the fact that less number of facilities are required in the optimal solution of the problem. Regarding the number of medians, the ordering is where the optimal number of facilities is closer to Med(SA) or Med(GMAX). en, we observe that the CPU times required to solve P s 1 increase significantly for values of α ∈ [1.7; 2]. Notice that, in particular, we cannot solve P s 1 to optimality in two hours for α � 1.8, 1.9 { }, which is not the case for the instance #11 in Table 1. In general, the CPU times required by SA are higher than those required by GL. Finally, we observe that the lower bounds obtained with LP s 1 are not tight when compared to the optimal solution of the problem which explains somehow the difficulty in solving P s 1 . Similarly, from Figure 5, we see that GMAX outperforms GMIN and IC in terms of upper bounds. However, SA and GL outperform the proposed heuristics. Regarding the number of medians obtained, we observe that SA and GL obtain less number of facilities than the heuristics. Finally, the CPU times obtained show similar orders of magnitude for all the algorithms. In Figures 6 and 7, we report analog numerical results as in Figures 2 and 3 for the instance #10 presented in Table 9 using radial coverage values of 0.2 and 0.5 km, respectively. In these figures, we do not report optimal solutions for P s 1 as it is not possible to obtain these results for the instances with higher dimensions due to CPLEX shortage of memory.
From Figure 6, we observe that the upper bounds obtained with GMIN, GL, and SA are nearly the same and better than the other ones. Regarding the number of facilities, we observe that GMIN obtains the largest values, whilst GMAX obtains the lowest ones. e number of facilities of each remaining algorithm is between these values. Finally, we observe that the CPU times obtained with SA and GL algorithms are the largest and smallest ones, respectively. Consequently, the CPU times for the remaining algorithms are between these values as well. From Figure 7, we observe similar trends for the upper bounds and nearly constant values for the number of facilities obtained with each algorithm. In particular, GMIN obtains the largest values. Subsequently, we observe that the CPU times are different for most of the algorithms and in particular higher for the metaheuristics in some cases.
In order to give more insights with respect to the behaviours of GMIN, GL, and SA algorithms, we further present average upper bounds and CPU times in seconds for the instance #11 in Table 1 and for the instance #10 in Table 9 while using radial coverage values of 0.2 and 0.5 km. ese numerical results are presented in Figures 8-11 for each of the instances, respectively. In particular, the averages are obtained for 50 runs using each algorithm.
From Figures 8 and 9, we observe that the upper bounds obtained with SA are lower than GMIN and GL. However, the upper bounds obtained with GL are lower than those obtained with GMIN. Regarding the CPU times, the average values can be ordered in the opposite direction for the three algorithms being larger for SA and lower for GMIN. In particular, in Figure 9, the averages for SA and GL are closer to each other.
Next, in Figure 10, we observe that the upper bounds obtained with SA are the lowest ones. However, GMIN obtains the worst solutions. Regarding the CPU times obtained, we see that, in average, GMIN and SA require similar values. However,GL obtains solutions faster than the remaining algorithms. Subsequently, in Figure 11, we observe that GL obtains better solutions compared to GMIN and SA although at a higher CPU time. Finally, we observe that GL allows us to obtain better solutions than SA for higher radial coverage values, i.e., when dealing with dense input disk graphs. However, on the opposite, SA allows us to obtain better solutions for sparse graphs which resulted in more difficult instances to be solved by the proposed MILP models.

Conclusions
In this paper, we consider the problem of assigning a set K of users to a subset S of facilities chosen from a larger set V while simultaneously forming a tree backbone with S and with the additional condition that no two adjacent nodes in V of an input disk graph G d can belong to S. e latter condition is handled by means of independent set constraints. We model this problem on unit disk graphs and propose mixed integer linear programming models in order to minimize the total connection cost distance between facilities and between customers and facilities simultaneously. Four compact polynomial formulations are proposed based on classical and set covering p-Median formulations, whilst the tree backbone formed with S is modelled with Miller-Tucker-Zemlin and path orienteering constraints, respectively. e MILP models are further strengthened with clique valid inequalities which can be obtained in polynomial time for unit disk graphs. Finally, we propose Kruskal-based heuristics and metaheuristics which are adapted from a greedy approach initially proposed for the maximum independent set problem. In particular, the metaheuristics are constructed based on guided local search and simulated annealing strategies. Our main observations and conclusions on this paper can be outlined as follows: (1) Our numerical results indicated that only the Miller-Tucker-Zemlin constrained models allow us to obtain optimal solutions for instances with up to 26 Complexity 200 nodes and 1000 users. Notice that finding feasible solutions with certificate of optimality is an important achievement. Indeed, it allows us to compare the solutions obtained with the proposed algorithms against optimal solutions and also to compare possibly new algorithmic approaches as part of future research. (2) We obtained the same near-optimal solutions with all linear relaxations and gaps lower than 6% with the strengthened models for most of the instances compared to the optimal solutions. (3) We compared Miller-Tucker-Zemlin and path orienteering constrained models using novel formulations which are constructed based on classical and set covering p-Median models. We observed, from our numerical results, that Miller-Tucker-Zemlin constrained models outperform path orienteering ones.
(4) We also proposed a heuristic framework with three variants referred to as GMIN, IC, and GMAX and noticed that GMIN outperforms the other ones for most tested instances. Furthermore, we noticed that, in some cases, when varying the weights of the objective function of the proposed MILP models, GMAX outperforms the other ones. In particular, this occurs when the second term of the objective function of P s 1 is greater than the first one, i.e., when the cost incurred to construct the backbone tree is greater than the cost of assigning users to facilities.  Table 9 while using a radial coverage value of 0.5 km.
Finally, we verified that the heuristic approach allows one to obtain near-optimal solutions in short CPU time. (5) Finally, we proposed two metaheuristics based on guided local search and simulated annealing greedy strategies (GL and SA, respectively) that outperformed the proposed heuristics for most of the instances. In particular, we noticed that SA obtains better results than GL on sparse disk graphs. However, the opposite occurs for dense graphs. Finally, both GL and SA allowed to obtain nearoptimal solutions in significantly short CPU time and tight feasible solutions for large instances of the problem.
As future research, we plan to propose new modelling and algorithmic approaches while taking into account other network structures such as dominating and maximum leaf spanning trees. Ultimately, new modelling approaches should also be considered while including mathematical majorization concepts related with distances in trees and location theory [45], in order to compare with the models proposed in this paper.

Data Availability
All data were generated randomly as it is clearly explained in the manuscript. Consequently, no particular data were used to support this study.

Conflicts of Interest
e authors declare that they have no conflicts of interest.