Abstract
Let be an input disk graph with a set of facility nodes V and a set of edges 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 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 and , respectively, where is the complement of . 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. The 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 the MTZ 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, the MTZ constrained models outperform path orienteering ones. However, the proposed heuristics and metaheuristics allow obtaining near-optimal solutions in significantly short CPU time and tight feasible solutions for large instances of the problem.
1. Introduction
Facility location has been a major research topic within last decades [1–6]. In general, a facility location problem is mainly concerned with the optimal placement of facilities in a predefined 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 fire 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 and be, respectively, the input disk and complete graphs composed of a set of facility nodes V drawn in the Euclidean plane with edge sets and , 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 with nodes in S where denotes the set of edges in spanned by T such that and with the additional condition that no two adjacent nodes in V related with can belong to S. The 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–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 with set of nodes V and set of edges E is a subset of vertices , no two of which are adjacent. This means that, for every two vertices in S, there is no edge connecting them. Notice that a maximal independent set is a subset 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 [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 . The 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. This 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–18]. In general, listing all maximal cliques in polynomial time can be achieved for sparse graphs [15, 19]. The 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]. The 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–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.
The paper is organized as follows: In Section 2, we review some related works. Then, 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.
2. 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. This 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 real-life applications. Their examples include distribution systems where transshipment among stocking points are 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–6, 30–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]. This 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. This 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. The 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]. The problem is solved with linear programming relaxations and branch and bound techniques as well. Variants of MCLP are also presented and discussed in [3].
The indirect covering tree problem proposed in [36] represents another approach which is similar in nature to our facility location problem. There, 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. The 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. The 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. The 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–32].
3. 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.
3.1. Description of the Problem
Consider the undirected graphs and as defined in Section 1. Notice that contains the edges obtained with all pairs of nodes in V, as is complete. However, is composed of the edges obtained with all pairs of nodes in V, which are located at a radial distance lower than or equal to d. The facility location problem we are dealing with consists of finding a subset of facilities such that S is an independent set of and a tree backbone of 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 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 and the independent set obtained from . Nodes in S are colored green. However, on the right, we 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 with low density, the cardinality of S should be higher. On the opposite, if the density of is high, the cardinality of S should be smaller. By density, we refer to the ratio obtained by dividing the number of edges of over the total number of edges of .
(a)
(b)
(c)
(d)
3.2. 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 as the set of directed arcs obtained from edges in and as the digraph obtained from . Next, we define the binary variables:
Thus, a MILP model can be verified by means of the following proposition.
Proposition 1. The linear model allows us to obtain a feasible solution with minimum cost for the facility location problem while satisfying the following conditions:(i)All users in K are assigned to facilities in S.(ii)Subset forms an independent set (not necessarily maximal) from .(iii)Subset forms a tree backbone from where is the set of edges spanned by T where :
Proof. To prove (i) and (ii), we focus on explaining how the set of constraints in together with its objective function imply each of the assertions. For this purpose, we define the input matrices and , , , to be the distances between user k and facility j and between facilities i and j, respectively. In particular, we require matrix to be symmetric. Thus, 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 , 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 . Notice that S is composed of the nodes for which the binary variable . 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 . Furthermore, note that S is not conditioned to be a maximal independent set.
To prove condition (iii), we show that any feasible solution of is an acyclic connected arborescence of . First, notice that the constraint (6) imposes the condition that any tree formed with nodes in S contains exactly arcs. However, the constraints (10) ensure that the number of arcs entering node is at most one if and only if . Similarly, the constraints (11) ensure that the number of incoming and outgoing arcs in node can be at most if and only if . Otherwise, if , then the binary variables , . 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 is solution to . Then, the constraints (7) imply that since . In turn, this implies that which cannot hold. Notice that, in general, the contradiction holds for any subset of nodes of S which implies that the digraph induced by 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 for some . This particular case is avoided by means of the constraints (6) and (10), which finally ensure that the digraph is connected with 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 , is symmetric. Finally, the constraints (12) limit the domain of the decision variables, thus concluding the proof.
Notice that 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 , we outline the following theorem.
Theorem 1. The facility location problem associated with is NP-hard.
Proof. To see this, let us consider the following weighted version of the objective function (2) aswhere . In particular, when , any spanning tree will be feasible for at zero cost for the objective function of the problem. Consequently, in this case, can be equivalently written as the following optimization problem:where the second quadratic sum in (14) is due to an equivalent form of writing independent set constraints [38]. For this purpose, parameter should be a nonnegative real value such that . In general, notice that any independent set S implies . Consequently, it is easy to see that is a quadratic version of the classical uncapacitated facility location problem which is NP-hard [39].
Subsequently, when , we aim to solve with the implicit requirement that the number of facilities must be between , where denotes the cardinality of the maximum independent set of . 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 , 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 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 is obtained with which is set by default. The 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 which reflects different cost structures.
An alternative formulation of can be obtained by removing constraints (7) and (9) together with the decision variables for each and by introducing the nonnegative precedence variables for all . Precedence variables are used to indicate if there exists a direct path from i to j. If such a path exists, then ; otherwise, . 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]. Thus, an equivalent MILP model can be written asConstraints (16) state that if there is a path from i to j, for all 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 , . 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 and must be equal.
Notice that both and 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 and , respectively. In order to write a first model under this approach, we first construct for each , the vectors by sorting in the ascending order the different entries of each k-th row of the original cost matrix in (2) in such a way that and . 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 and with this additional node. Thus, we expand sets V, , and and denote them by , , and , respectively. Since the artificial node has to be excluded from the solution of the problem, the extra edges added to and are irrelevant for the solution of the problem. Similarly, we construct a new matrix , , 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 asIn , notice that the binary variable for each , , acts as a cumulative variable where if and only if the allocation cost distance of user k is at least and otherwise. For this purpose, the set covering constraints (24) ensure that variable if there is no open facility at less than distance . The positive coefficients of the first double sum in the objective function (23) ensure that variable 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 . Analogously as for , 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 , , , and , 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 . Recall that, in practice, all cliques can be obtained in polynomial time as is a unit disk graph [14, 15]. Consequently, each row in the binary matrix corresponds to a maximal clique . 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].
From the path orienteering-based models, we only remove constraints (5) and (11), and replace them with (27) and (29) and (30), respectively. Hereafter, we denote the linear programming (LP) relaxations of , , , , , , , and by , , , , , , , and , respectively.
4. 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–29].
4.1. Heuristic Approach
Three variants of the same heuristic are proposed. The only difference between them relies on the greedy decision of adding new nodes to S at each iteration. Thus, we only present the first one referred to as in Algorithm 1 and explain the other ones based on these steps. The heuristic requires the input matrices , , , , and the input disk graph . Then, for each vertex , the following steps are performed. The set S is initialized as an empty set, and W is assigned set V. The set W is used as an auxiliary set. Next, 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 to its nearest facility in S. Finally, we compute the objective function value of 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. Then, in Step 2, we enter into a while loop and iterate until is empty, where denotes the set of edges in induced by remaining nodes in W. Inside the while loop, we proceed similarly by selecting a new vertex with minimum or maximum degree to be added to S or by interchanging between minimum and maximum degrees. This is the greedy decision which differentiates remaining variants of Algorithm 1. Hereafter, we denote the three heuristics by , , and , respectively. Finally, Algorithm 1 returns the best solution found and its objective function value.
|
4.2. Guided Local Search Approach
The 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 and , and iterate while is not empty. The following steps are performed inside the loop. First, we remove the node with minimum degree from W and add it to S. Then, we remove all neighbors of from W and find the minimum spanning tree with current nodes in S using the Kruskal algorithm. Finally, we assign each user to its nearest facility in S, compute the objective function value of , and save 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 which is an input parameter. Inside this loop, we penalize each value in vector corresponding to each node in by adding a random number uniformly distributed in . Initially, contains the degree of each node . 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. Then, we reset and and construct a new independent set S from 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 to its nearest facility in S, and compute the objective function value of . If this value is less than the best found so far, we update and reset the current CPU time to zero. Finally, we randomly generate a value in order to decrease, maintain, or increase, respectively, the cardinality of the next independent set to be generated by the algorithm starting from . Notice that Algorithm 2 updates degree values in 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].
|
4.3. 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–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. The 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. The 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 temperature is set to a large positive value . Inside the while loop of Step 2, we set and , and for each , which is obtained in step one, we pick randomly a number from and compare it with the input parameter η which is arbitrarily chosen from the interval as well. If the random number is less than η, then we include node in the new independent set S and remove from W node and its neighbors. Otherwise, we skip node and continue. This allows us to keep an η percentage of nodes in . Next, we start a while loop in order to add new nodes in the solution until is empty or . Added nodes are chosen randomly from remaining nodes in W. Subsequently, we find the minimum spanning tree with nodes in S, assign each user to its nearest facility in S, and compute the objective function value of . If a better solution is found, we save this new solution in and set 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 , and compute the probability . If r is less than p, we accept the new solution and save it in ; else, we set . 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 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. The next lines inside the while loop of Step 2 ensures that decreases at a δ rate. Then, we check if where ϵ is a small positive value. If so, we reset and allow the algorithm to interchange between exploration and exploitation more frequently during the process. Finally, we randomly generate a value to decrease, maintain, or increase, respectively, the size of the next independent set to be generated by the algorithm. Notice that the value of is of critical importance to handle efficiently the algorithm. If decreases rapidly, then a premature convergence to a local minimum may occur. On the opposite, if 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, 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 at iteration t by where [43].
|
5. 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 and for all , , represents the distance between user k and facility j and distance between facilities i and , respectively. Matrix D is symmetric. We generate two sets of instances with dimensions of users and facility nodes for radial coverage values of 0.2 and 0.5 km. Each set is composed of 20 instances. Then, we further generate two larger sets of instances with dimensions of users and nodes for both radial coverage values of 0.2 and 0.5 km as well. The larger sets are composed of 15 instances each one and are only solved with the proposed algorithms since the linear models cannot be solved with CPLEX due to shortage of memory [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 , , , , , , , , , , , , , , , and [44]. The proposed algorithms , , , , and are also implemented in Matlab, where refers to Algorithm 2 and refers to Algorithm 3, respectively. The parameters in Algorithms 2 and 3 are calibrated to , , , , , and . The numerical experiments have been carried out on an Intel(R) 64 bits core (TM) with 3.40 GHz and 8 gigabytes of RAM.
5.1. Numerical Results for the Linear Models
In 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. Thus, 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]. The column information in these four tables is exactly the same with the exception of Table 1 which includes four additional columns with instance dimensions. The information of each column is described as follows. In column 1, we present the instance number. Then, 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 . The 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 where refers to the optimal solution of each MILP model and represents the optimal solution obtained with its linear relaxation.
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 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. This 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 and are exactly the same as those obtained with and 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. This observation is also valid for the LP relaxations which require more CPU times than for solving and , 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 , , , and , 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. Finally, from the numerical results presented in Tables 1–4, we conclude that model 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 , where . Notice that 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 [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:
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). This 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 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 . This observation can be understood by simply looking at the objective function of the proposed 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 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 used in Theorem 1. Numerical results while varying parameter α are discussed and reported below.
5.2. Numerical Results for the Algorithms
In 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.
Thus, 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 solution found with or together with its minimum CPU time in seconds, respectively. Then, 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 , , and , respectively. The gaps are computed by .
From Table 7, first we observe that the number of facilities is nearly the same as the optimal one for heuristics. In general, allows us to obtain a higher number of facilities than both and , whereas obtains more facilities than . We note that this can be explained as a consequence of the greedy decision that each heuristic performs at each iteration. Notice that removes less nodes from the graph than , which in turn removes less nodes than . Next, we further observe that obtains tighter bounds than the other heuristics although at a higher CPU time effort. This 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 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 require significantly less CPU time than those obtained with or . Notice the fact that the first term in the objective function of 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 is different, later we further present numerical while varying the parameter using the weighted objective function presented in Theorem 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 and , respectively. More precisely, in these columns, we report number of facilities, upper bounds, and CPU times in seconds and gaps that we compute by .
From Table 8, we observe that the number of facilities is slightly larger for than for 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 are tighter than those obtained with 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. The gap columns confirm these improvements. Finally, we observe that the CPU times obtained with and 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. Then, in columns 5–7, we report number of facility nodes, upper bounds, and CPU times in seconds for 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 and , 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 . More precisely, we compute the gaps as and 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 are tighter than those obtained with , which in turn are tighter than those obtained with for most of the instances when using a radial coverage of 0.2 km. On the opposite, obtains tighter gaps than and when using a radial coverage of 0.5 km. In particular, a few negative gaps are reported when allows us to obtain better solutions compared to a metaheuristic. This 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 compared to for dense instances. However, the opposite occurs for the sparse ones. Notice that the CPU times obtained with 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 , its maximum clique number is obtained by solving the optimization problem presented above using the set which corresponds to the set of edges of the complement graph of .
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 , , , , , , and algorithms while varying α using the objective function of Theorem 1, in Figures 2–7, we plot upper bounds, optimal solutions, and number of facilities and CPU times in seconds for the instances #11 and #10 in Tables 1 and 9, respectively, while varying 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 in order to decrease the number of optimal facilities while searching for a tradeoff between the two terms in the objective function.
(a)
(b)
(c)
(a)
(b)
(c)
(a)
(b)
(c)
(a)
(b)
(c)
(a)
(b)
(c)
(a)
(b)
(c)
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: UbUbUbUbUb. Regarding the number of facilities (or medians), we observe that the order of the methods is as follows: MedMedMed Med Med, where the optimal number of facilities is closer to Med. However, the order of the methods with respect to the CPU times is as follows: CPUCPUCPU. For the heuristics, the CPU times remain nearly the same.
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 . In particular, the number of facilities is exactly the same for all the methods. Finally, from Figures 2 and 3, we observe that the number of facilities decreases for values of . Consequently, in Figures 4 and 5, we plot for the same instances numerical results as in Figures 2 and 3, respectively, but for values of .
From Figure 4, we observe the following ordering for the upper bounds: UbUbUb UbUb. As it can be observed, in this case, allows us to obtain better solutions than . This 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 MedMed, MedMedMed, where the optimal number of facilities is closer to Med or Med. Then, we observe that the CPU times required to solve increase significantly for values of . Notice that, in particular, we cannot solve to optimality in two hours for , which is not the case for the instance #11 in Table 1. In general, the CPU times required by are higher than those required by . Finally, we observe that the lower bounds obtained with are not tight when compared to the optimal solution of the problem which explains somehow the difficulty in solving . Similarly, from Figure 5, we see that outperforms and in terms of upper bounds. However, and outperform the proposed heuristics. Regarding the number of medians obtained, we observe that and 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 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 , , and are nearly the same and better than the other ones. Regarding the number of facilities, we observe that obtains the largest values, whilst obtains the lowest ones. The number of facilities of each remaining algorithm is between these values. Finally, we observe that the CPU times obtained with and 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, 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 , , and 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. These 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.
(a)
(b)
(a)
(b)
(a)
(b)
(a)
(b)
From Figures 8 and 9, we observe that the upper bounds obtained with are lower than and . However, the upper bounds obtained with are lower than those obtained with . Regarding the CPU times, the average values can be ordered in the opposite direction for the three algorithms being larger for and lower for . In particular, in Figure 9, the averages for and are closer to each other.
Next, in Figure 10, we observe that the upper bounds obtained with are the lowest ones. However, obtains the worst solutions. Regarding the CPU times obtained, we see that, in average, and require similar values. However, obtains solutions faster than the remaining algorithms. Subsequently, in Figure 11, we observe that obtains better solutions compared to and although at a higher CPU time. Finally, we observe that allows us to obtain better solutions than for higher radial coverage values, i.e., when dealing with dense input disk graphs. However, on the opposite, allows us to obtain better solutions for sparse graphs which resulted in more difficult instances to be solved by the proposed MILP models.
6. 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 can belong to S. The 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. The 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 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 , , and and noticed that 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, outperforms the other ones. In particular, this occurs when the second term of the objective function of 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. 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 ( and , respectively) that outperformed the proposed heuristics for most of the instances. In particular, we noticed that obtains better results than on sparse disk graphs. However, the opposite occurs for dense graphs. Finally, both and allowed to obtain near-optimal 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
The authors declare that they have no conflicts of interest.
Acknowledgments
The authors acknowledge financial support from projects FONDECYT (nos. 11180107 and 3190147).