Globally Optimal Facility Locations for Continuous-Space Facility Location Problems

The continuous-space singleand multi-facility location problem has attracted much attention in previous studies. This study focuses on determining the globally optimal facility locations for twoand higher-dimensional continuous-space facility location problems when the Manhattan distance is considered. Before we propose the exact method, we start with the continuous-space single-facility location problem and obtain the global minimizer for the problem using a statistical approach. Then, an exact method is developed to determine the globally optimal solution for the twoand higher-dimensional continuous-space facility location problem, which is different from the previous clustering algorithms. Based on the newly investigated properties of the minimizer, we extend it to multi-facility problems and transfer the continuous-space facility location problem to the discrete-space location problem. To illustrate the effectiveness and efficiency of the proposed method, several instances from a benchmark are provided to compare the performances of different methods, which illustrates the superiority of the proposed exact method in the decision-making of the continuous-space facility location problems.


Introduction
Over the last five decades, the facility location problem, also known as location analysis, has attracted much attention in mathematical science [1]. A large number of researchers have investigated both the formulations and the algorithms for diverse applications in the private and public sectors [2][3][4]. Concerning the private sector, organizations have to continuously search for new ways to reduce costs, improve customer satisfaction, and increase profitability due to global competition. Logistics operations, especially facility locations, such as industrial plants, bicycle-sharing stations, banks, distribution centers, warehouses, and fourth-generation/fifth-generation (4G/5G) base stations, have traditionally been an opportune field for cost-saving. Similarly, for the public sector, choosing appropriate locations for facilities, such as hospitals, ambulance stations, post stations, transport terminals, medical service centers, and relief centers, enables them to improve the service level and efficiency. Generally, a better option needs to be done to look for compromises on behalf of different goals.
As the travel cost and time can be analyzed by discrete or continuous aspects in the space, the facility location problem is commonly divided into two types, namely (i) the discrete-space facility location problem and (ii) the continuous-space facility location problem. Practically, for the discrete-space type, the locations of potential facilities can just be located among the given specific points, whereas in the continuous-space facility The facility location problem is an important issue in supporting physical distribution as it contributes significantly to the travel time or cost in logistics systems. Since the facility location problem can be divided into two types, namely discrete-space and continuous-space facility location problems, this study focuses on the continuous-space facility location problem.
For the continuous-space facility location problem, anywhere in the plane can be considered as the facility location [5]. Thus, there is an infinite number of potential locations for the facilities, thus making it quite difficult to solve [13]. Weber and Friedrich [14] have tried to locate a single industry or firm to minimize the transportation cost in the plane, known as the Weber problem, which has received much attention in the literature. Moreover, many studies have been developed for locating the facilities in the continuous space to minimize or maximize objective criteria from different perspectives [15][16][17][18][19][20]. However, due to an infinite number of possible locations for the facilities, it is a big challenge computationally, unless some appropriate approaches or algorithms can be developed.
To solve the above problem, many approaches have been developed previously [21,22], e.g., by Zhang [23]. The earliest approach is an iterative procedure that was offered by Weiszfeld [24], which was followed by some other variants [25][26][27]. The geometrybased approach is also used for solving continuous-space facility location problems [28]. Since the 1990s, the Voronoi diagram heuristic approach has been applied for solving the p-center problem in a continuous space [29,30]. Rather than searching for the optimal facility locations from an infinite number of possible locations, making a list of some good places that are considered as potential facility candidates helps solve the problem. Instead of large-scale computation to obtain the optimal solution, it is necessary to identify a finite set that contains the optimal solution, which is also known as the finite dominating set. As a consequence, search effort can be reduced dramatically so that large-scale computation is avoided. In this research area, the nodes of the network are identified as the finite dominating set for the p-median problem [31,32]. However, when the finite dominating set is very large, it is still computationally burdensome to find the optimal solution.
To overcome the above difficulty, more recently, the facility location problem is defined as a special clustering problem, where the sets of customers served by the same facility are considered as clusters [33]. The clustering-based algorithms, such as k-means clustering and fuzzy C-means (FCM) clustering algorithms, have been widely applied in previous studies. Although the clustering-based algorithms cannot guarantee the execution time and solution quality, they have proven useful in practice. Žalik [34] applied the FCM clustering algorithm to minimize the mean squared distance from each data point to its nearest center. Sheu [35,36] developed different versions of hybrid fuzzy clustering algorithms to group customer demands. Esnaf and Küçükdeniz [37] presented a fuzzy clustering-based hybrid method for a multi-facility location problem where the capacity of each facility was unlimited. Chen, Yeh [38] considered the Euclidean distance in searching for potential locations of temporary emergency medical centers using a clustering-based algorithm. Varghese and Gladston Raj [39] applied the k-means clustering algorithm so solve a multifacility location-allocation problem. Gao, Nayeem [40] also developed a clustering-based genetic algorithm for a multi-facility location problem, where the Manhattan distance was considered.
After clustering the customers or clients into several groups, the facility locations need to be identified in each of the groups. In this sense, the center-of-gravity (CG) method is one of the most widely used approaches in previous studies [41][42][43], which is also known as the weighted mean. In the previous studies using the CG methods, the goals are usually to minimize an objective function involving squared Euclidean distance, Euclidean distance, or Manhattan distance from the facilities to the demand points. For instance, Ohsawa [44] used the CG method so that the average squared Euclidean distance can be minimized from the facility to the demand points. Esnaf and Küçükdeniz [37] and Nadizadeh, Sahraeian [45] also applied the CG method to obtain the optimal facility location with the consideration of Euclidean distance. Onnela [46] also applied the CG method to minimize the total weighted Manhattan distance and squared Euclidean distance. However, the CG method can only be applied to minimize the total weighted squared Euclidean distance rather than the total weighted Manhattan distance. Accordingly, the CG method is invalid anymore to minimize the total weighted Manhattan distance.
In addition to the above discussions, we also provide a comprehensive review in Table 1, where relevant studies about this topic are summarized. The differences with previous studies are summarized in three aspects. Particularly, many efforts have been dedicated to various continuous-space multi-facility problems considering Euclidean distance. In this study, the Manhattan distance is considered and the optimal facility locations are determined. Besides, most of the previous studies applied a variety of clustering algorithms to obtain the facility locations, which can be considered as a generalized multi-Weber problem and referred to as an uncapacitated multi-facility location-allocation problem as stated by Copper [47]. Moreover, the problem can be interpreted as an enumeration of the Voronoi partitions of the customer set, which has been proven to be an NP-hard problem [48,49] and cannot guarantee the global optimum of the solutions. As a consequence, the importance of obtaining the globally optimal solution remains an urgent problem, posing the question: How do we find the set of facility locations in a continuous plane that is the global optimum so that the objective function can be minimized?
However, determining the global optimum of the facility locations in a continuous plane is quite difficult because of the infinite number of potential locations and continuous distance. Furthermore, different weights at demand points make the multi-facility determination challenging. Neither strategic approaches nor quantitative models to overcome the pitiful clustering algorithms in handling the continuous-space facility location problems have been investigated before. Therefore, we focus on determining the global optimum of the solutions for the continuous-space multi-facility location problem, which

Location Estimation Based on the Manhattan Distance
In this section, we first illustrate the Manhattan distance and present the objective function based on the Manhattan distance. Then, we apply a statistical approach to obtain the global minimizer for the continuous-space single-facility location problem, where several properties of the minimizer can be obtained. With these properties of the minimizer, we can drive optimization procedures in searching for the globally optimal facility locations.

Objective Function Based on the Manhattan Distance
The Manhattan distance is usually applied in the urban area of interest due to the street configuration [40,56]. For instance, the determination of the bicycle-sharing stations or the ambulance stations needs to consider the Manhattan distance. For more details on the Manhattan distance, please refer to Gao and Cui [7]. Based on the Manhattan distance, the general objective function (minimization of the total weighted distance) can be formulated with the following notations.

Parameters:
• n, Number of demand points, (i = 1, 2, . . . , n); Generally, the objective function needs to be minimized so that the total weighted distance-related travel time or cost is optimized. Before the objectives are presented, the weights need to be normalized using the following equation.
where ω i is the normalized weight of demand point i. Then, the objective function for the two-dimensional single-facility location problem is given by where |x i − u| + |y i − v| is the distance measured along the axis at right angles. Please note that the normalized weight influences the objective function value, but it would not affect its final solution. Similarly, the same situation happens in the other objective functions. The objective function for the three-dimensional single-facility location problem is given by where |x i − u| + |y i − v| + |z i − w| is the distance measured along the axis at right angles, which is usually applied in the mansion building of interest and three-dimensional facility cases.

Minimum Distance Approach
In this subsection, we provide the three-dimensional single-facility location methods for the Manhattan distance. With the objective function of Ψ 3 in (3), let ( u, v, w) be the optimal facility location. Then, we have The minimizer of Ψ 3 in (3), denoted by ( u, v, w), can be obtained by the following estimating equations that need to be solved for u, v, and w according to Section 1.3 of Hettmansperger and McKean [57].
Then, the optimal values u, v, and w can be calculated separately. It is easily seen that they are the weighted medians of the x-axis, y-axis, and z-axis observations, which will be detailed later. Note that the weighted median was first suggested by Edgeworth [58] and since then it has been widely used in many applications [59]. As an illustration, we briefly introduce the conventional median first and then the weighted median. As the values u, v, and w can be obtained separately, we only consider the weighted median for the x-axis observations. Then, the weighted medians for the y-axis and z-axis observations are easily obtained using the same method. Definition 1. Given a set of observations x 1 ,x 2 , . . . , x n , the empirical cumulative distribution function F n is defined as where I(A) represents the indicator function defined as Definition 2. Csorgo [60] defined the left sample quantile function (inverse cumulative distribution function) as below.
Using this, the conventional median F n −1 (1/2) is obtained as: where Definition 3. Wasserman [61]defined the right sample quantile function (inverse cumulative distribution function), which is given by Note that the definition by Wasserman [61] is slightly different to that of Csorgo [61]. It is easily seen that Thus, the sample quantile by Csorgo [60] is called the left quantile, while the sample quantile by Wasserman [61] is the right quantile. Rychlik [62] and Hosseini [63] showed that Based on the quantile function by Wasserman [61], we have the corresponding median F n −1 (1/2), which is obtained as: where It is obvious that the difference between (13) and (15) is the case when F n (1/2) = 0.5 with even n. It is obvious that (13) takes the left bound value and (15) takes the right bound value. Moreover, both of them are the minimizers (medians) to the total Manhattan distance of the x-axis observations.
The above definitions on the empirical distribution and the sample quantile do not consider the weights of the observations. Thus, the definitions with weights are given below.

Definition 4.
Given a set of observations x 1 ,x 2 , . . . , x n with corresponding positive weights ω 1 ,ω 2 , . . . , ω n such that ∑ n i =1 ω i = 1, we have the empirical cumulative distribution function G n (x) with weights, which is defined as Note that the above G n includes the conventional empirical cumulative distribution function F n in (8), as a special case when ω i = 1/n. Similar to the definition of the sample quantile function in Csorgo [60], we define the sample quantile function with weights. Definition 5. Given a set of observations x 1 ,x 2 , . . . , x n with corresponding positive weights ω 1 ,ω 2 , . . . , ω n such that ∑ n i =1 ω i = 1, the sample left and right quantiles with weights are given by Next, our goal is to obtain the minimizer of the Ψ 3 in (3). As the weighted medians for the x-axis, y-axis, and z-axis observations can be calculated separately, we only focus on the weighted median of x-axis observations, which are given by where ω (j) is the weight for x (j) . According to the equations in (14) and (20), we can deduce the following equation. Then, we have It is obvious that the above x (k) minimizes the weighted Manhattan distance. However, it is not a unique minimizer when there is a tied value at ∑ k j=1 ω (j) = 1/2. Thus, we have two cases and both of them are the minimizers to the objective function: j=1 ω (j) < 1/2 and 1/2 = ∑ k j=1 ω (j) Because the minimizer of Ψ 3 , denoted by u, v, and w can be calculated separately. Based on the weighted median obtained in this subsection, the globally optimal location of the facility in a three-dimensional space is given by In addition, the globally optimal location of the facility in a two-dimensional plane is easily obtained by (22) and (23). For the higher-dimensional (d = 4, 5, . . . ,) globally optimal facility location problem with Manhattan distance, the main task is to calculate each of the weighted medians for all the axes. Then, the optimal facility location in a higher-dimensional hyperspace can be obtained.

Properties of the Minimizer
Next, we derive the properties of the globally optimal locations for the continuousspace multi-facility location problem. As obtained of the minimizer for Ψ 3 , the globally optimal single-facility location is the weighted median of the axis observations. In the first case, there is only one minimizer x (k) , whereas the second case contains an infinite number of minimizers between x (k) and x (k+1) . Suppose that we choose only one minimizer among x (k) , x (k+1) , (x (k) + x (k+1) )/2 for the second case, we can construct the candidate facility locations. To better observe the characteristic of the optimal single-facility location, we use four examples (i.e., Examples 1-4) to illustrate the property of the optimal facility location, and the datasets are provided in Table 2.
As shown in Figure 1, the black circles stand for the demand points and black hollow squares represent the optimal facility locations. In Example 1 (see Figure 1a), the optimal facility location is located at the second demand point due to its coordinate values and weights, which are given by As shown in Figure 1c, there is only one specific location for the facility due to the singular weighted median on each of the axis. The optimal facility location does not locate at the demand point, which is given by  There are four optimal facility locations (black hollow squares) in Example 2 (see Figure 1b) because of weight construction. Then, the weighted medians are given by It is easy to obtain the optimal facility locations through permutation and combination of the weighted medians on different axes, which are (3,3), (3,6), (3, 4.5), (5,3), (5,6), (5, 4.5), (4,3), (4,6), and (4, 4.5). It is obvious that two of them, i.e., (3,3) and (5,6), are located at the demand points, but the others are not. As they have the same objective function value, any one among them can be considered as the location for the facility.

Example 3
As shown in Figure 1c, there is only one specific location for the facility due to the singular weighted median on each of the axis. The optimal facility location does not locate at the demand point, which is given by

Example 4
In this case, the weighted median on the axes are given by Thus, there are three black hollow squares in Figure 1d, where any of them can be considered as the optimal location for the facility.
As presented in Figure 1 for Examples 1-4, we can conclude that the globally optimal facility location is based on the coordinate values of the observations on different axes. After constructing the set of mesh points that are consist of the coordinate combinations, it is obvious that the optimal single-facility location is an element from the set of mesh points. Then, we have the following formula Similarly, if we choose only one minimizer among x (k) and x (k+1) for the second case in (39), the optimal single-facility location is also an element from the other set of mesh points, which is given by Note that the number of candidate locations in (34) is smaller than that in (33) but the objective function obtained through (34) is the same as that through (33). In this sense, the number of candidate locations in (34) is suggested in this study.

Identification of the Candidate Locations
As presented in (33), the globally optimal single-facility location is an element from the set of mesh points. Then, we investigate the candidate locations for the multi-facility location problem. In the multi-facility location problem, each demand point is assigned to its closest facility. Suppose that we have obtained the optimal facility locations, then the set of demand points served by the same facility is considered as a group, which is independent of other groups. In this sense, the optimal facility location is still one of the mesh points within this group. After considering the optimal facility locations for all groups, these optimal facility locations are from the mesh points that are constructed based on the observation coordinates. Then, we have where m is the number of facilities, which is indexed by l.
As an illustration, we provide an example to show the candidate locations (i.e., mesh points) for the facilities. In the example, we consider three demand points (see Table 3). As shown in Figure 2, we construct the mesh points based on these three demand-point coordinates. Accordingly, 25 mesh points are represented by black hollow rhombuses that also contain three demand points. These 25 mesh points are considered as the candidate facility locations. It is easy to conclude that the number of candidate facilities E is related to the number of demand points n, where the relationship between them is given by   Based on the conclusion about the number of the candidate locations for the twodimensional multi-facility location problem, it is easy to deduce the number of candidate locations for the three-dimensional multi-facility location problem, which is given by If we consider the minimizer from ( ) and ( ) for the second case presented in (39), the number of mesh points can be reduced and the number of candidate facilities E for the two-and three-dimensional multi-facility location problems are given below Finally, we construct the candidate facility locations and transfer the continuousspace multi-facility location problem to the discrete-space multi-facility location problem. Next, our goal is to select the optimal locations from the mesh points, which is the global optimum solution to the problem. When the demand points have the same coordinate value on one axis, the number of candidate facilities E can be smaller. Suppose that there are p (p ≤ 2n − 1) different x-axis values and q (q ≤ 2n − 1) different y-axis values, the number of candidate facilities E can be smaller, which is given by Based on the conclusion about the number of the candidate locations for the twodimensional multi-facility location problem, it is easy to deduce the number of candidate locations for the three-dimensional multi-facility location problem, which is given by If we consider the minimizer from x (k) and x (k+1) for the second case presented in (39), the number of mesh points can be reduced and the number of candidate facilities E for the two-and three-dimensional multi-facility location problems are given below Finally, we construct the candidate facility locations and transfer the continuous-space multi-facility location problem to the discrete-space multi-facility location problem. Next, our goal is to select the optimal locations from the mesh points, which is the global optimum solution to the problem.

Determination of the Globally Optimal Facility Locations
After the construction of the candidate locations for the facilities, what we need to do is to select the optimal facility locations from the mesh points. We formulate the mathematical model for the multi-facility location problem, which is also called the location-allocation problem. Besides, the construction of the candidate locations for the facilities makes it possible to determine the globally optimal locations for the capacitated facilities. Then, we present the different models for different facility location problems.

Model for the Uncapacitated Multi-Facility Problem
With the number of demand points n, we have E candidate facility locations. Let I be the set of demand points, which is indexed by i with i ∈ I. Let L be the set of candidate facility locations, which is indexed by l with l ∈ L. Then, the distance from the candidate facility location l to the demand point i is known, which is denoted by D li . Given the number of facilities m > 1, our goal is to select m locations among E candidate facility locations and allocate the demand points to these facilities to minimize the total cost. To formulate this problem, we need to define two more decision variables, which are given by g l 1 if candidate facility l is selected 0 otherwise t il 1 if demand point i assigned to candidate facility l 0 otherwise Then, we have the following mathematical model 1 for the general multi-facility location-allocation problem: 1: The objective function in (41) minimizes the total weighted travel cost. Constraint in (42) means that we must locate exactly m facilities. Constraint (43) states that a demand point can only be serviced by one facility. Constraint (44) indicates that the demand point can only be assigned to the opened facility.

Model for the Uncapacitated Multi-Facility Problem with Fixed Cost
In the above mathematical model, the fixed cost of opening a facility is not considered. However, in practice, the fixed cost of opening a facility is inevitable. Thus, the fixed cost, denoted by F, needs to be considered in the multi-facility location-allocation problem. In this sense, the number of facilities is also a decision variable. Then, the mathematical model 2 is given by 2: The objective function in (45) minimizes the total cost including the weighted travel cost and fixed cost. The first constraint in (46) means that the number of facilities cannot be more than n. Constraints (47) and (48) are the same as the Constraints (43) and (44).

Model for the Capacitated Multi-Facility Problem with Fixed Cost
Generally, the clustering-based method is applied to solve the uncapacitated continuous-space multi-facility location problem. Even though some studies have applied the adjusted clustering-based algorithms to solve the capacitated multi-facility problem [30,31], it is still impossible to guarantee the global optimum of the solution. In this sense, obtaining the globally optimal locations for the facilities is warranted when the practical factor of the limited capacities is considered. Based on the constructed candidate locations for the facilities, the globally optimal locations for the capacitated facilities also can be selected from the mesh points. Let r i be the quantity of demand at demand point i and Cap l be the capacity of the facility l. Note that r i can be considered as the weight of demand point i. Then, the capacitated continuous-space multi-facility problem can be formulated as the following mathematical model 3.

3:
Min The objective function in (49) minimizes the total cost including the travel cost and fixed cost. The first constraint in (50) guarantees the maximum number of facilities. Constraints (51) and (52) are the same as the Constraints (43) and (44). Constraint (53) restricts the facility capacity.

Numerical Examples
In this section, several illustrative instances from a benchmark composed of 27 instances (http://neo.lcc.uma.es/vrp/, accessed on 1 May 2020) referred to Augerat, Belenguer [64] are provided to evaluate the performance of the proposed method. Each instance provides customers' (demand point) geographical coordinates and quantities of demand. The proposed exact approach is used to determine the globally optimal facility locations under different numbers of facilities for four instances (see Table 4) from the benchmark (i.e., Instances 24-27). Given the parameter values, the proposed models (i.e., 1, 2, 3) are implemented in the IBM ILOG CPLEX Optimization Studio (Version: 12.6). All the experiments are run on a computer with an Intel(R) Core(TM) i7-7700 CPU@3.6 GHz and 8 GB memory under the Windows 10 Pro system. In what follows, the numerical results are provided to validate the effectiveness of the proposed method. Then, the proposed method is compared with the previous clustering-based methods, where the comparison results are also provided to illustrate the advantages of the proposed method.

Numerical Results
First of all, the model 1 is applied and the facility locations are presented in Figures 3-6 for Instances 24-27. As shown in Figures 3-6, different numbers of the facilities (i.e., 3, 4, 5, and 6) are tested for each of the instances, where the black hollow triangle stands for the demand point and the solid square is represented as the facility. It is obvious that the facility locations are well scattered in the planning area. Since the difference between the models 1 and 2 is subject to whether the fixed cost is considered, the same results but with different objective function values are obtained based on the models 1 and 2.

Numerical Results
First of all, the model ℳ1 is applied and the facility locations are presented in Figures 3-6 for Instances 24-27. As shown in Figures 3-6, different numbers of the facilities (i.e., 3, 4, 5, and 6) are tested for each of the instances, where the black hollow triangle stands for the demand point and the solid square is represented as the facility. It is obvious that the facility locations are well scattered in the planning area. Since the difference between the models ℳ1 and ℳ2 is subject to whether the fixed cost is considered, the same results but with different objective function values are obtained based on the models ℳ1 and ℳ2.  Similarly, the globally optimal facility locations can also be obtained when the facility capacities are considered in model 3. As shown in Figure 7, the facility locations are determined for Instance 24 when different numbers of facilities are considered, which verifies the effectiveness of the proposed approach in obtaining the globally optimal facility locations, even though the facility capacities are involved. Note that similar results can be obtained when Instances 25-27 are considered. It is also obvious that the globally optimal facility locations obtained from model 3 are different from those in models 1 and 2, which indicates that the facility capacities have a great influence on the multi-facility location determination problem.   Next, the influence of different numbers of facilities and fixed costs on the objective function values are tested for models 1, 2, and 3. Regarding model 1, the objective function value is only affected by the number of facilities. In this sense, the relationship between the number of facilities and the objective function value is presented in Figure 8 given different shipping costs. It is obvious that the increase in the number of facilities always decreases the total cost, no matter how the shipping cost changes.  When the fixed cost is considered at facilities, in addition to the number of facilities, the fixed cost also has a great influence on the objective function value. To further provide a visual illustration of the influence of the parameters m and F on the total cost, a crossvalidation experiment is carried out when the number of facilities m goes from 2 to 18 and fixed cost F goes from 60 to 360. Then, the results under different shipping costs are presented in Figures 9-13. Specifically, as shown in Figure 9, the growth in the fixed cost from 60 to 360 increases the total cost. However, when the number of facilities m goes from 2 to 18, the total cost decreases first, since more facilities lead to lower transportation costs, and then the total cost increases because too many facilities contribute to marginally lower transportation costs, but result in higher fixed costs. Similar results are obtained in Figures 10-13. As a consequence, an appropriate number of facilities needs to be determined so that the globally minimum total cost can be obtained.
To investigate the influence of the capacitated facilities on the optimal facility locations and objective function values, the proposed three models 1, 2, and 3 are compared in terms of Instance 24 and the corresponding comparison results are provided in Table 5. Regarding the models 1 and 2, the same facility locations are obtained because the fixed cost does not influence the determination of the facility locations. When the facility capacities are considered, the globally optimal facility locations in model 3 can still be obtained, which is different from that in model 2. Besides, a higher total cost is obtained in model 3 compared with model 2, which verifies that the facility capacities affect the determination of facility locations. Similarly, the globally optimal facility locations can also be obtained when the facility capacities are considered in model ℳ3. As shown in Figure 7, the facility locations are determined for Instance 24 when different numbers of facilities are considered, which verifies the effectiveness of the proposed approach in obtaining the globally optimal facility locations, even though the facility capacities are involved. Note that similar results can be obtained when Instances 25-27 are considered. It is also obvious that the globally optimal facility locations obtained from model ℳ3 are different from those in models ℳ1 and ℳ2, which indicates that the facility capacities have a great influence on the multi-facility location determination problem.  jective function value is only affected by the number of facilities. In this sense, the relationship between the number of facilities and the objective function value is presented in Figure 8 given different shipping costs. It is obvious that the increase in the number of facilities always decreases the total cost, no matter how the shipping cost changes. When the fixed cost is considered at facilities, in addition to the number of facilities, the fixed cost also has a great influence on the objective function value. To further provide a visual illustration of the influence of the parameters and F on the total cost, a crossvalidation experiment is carried out when the number of facilities goes from 2 to 18 and fixed cost F goes from 60 to 360. Then, the results under different shipping costs are presented in Figures 9-13. Specifically, as shown in Figure 9, the growth in the fixed cost from 60 to 360 increases the total cost. However, when the number of facilities goes from 2 to 18, the total cost decreases first, since more facilities lead to lower transportation costs, and then the total cost increases because too many facilities contribute to marginally lower transportation costs, but result in higher fixed costs. Similar results are obtained in Figures 10-13. As a consequence, an appropriate number of facilities needs to be determined so that the globally minimum total cost can be obtained.        To investigate the influence of the capacitated facilities on the optimal facility locations and objective function values, the proposed three models ℳ1, ℳ2, and ℳ3 are compared in terms of Instance 24 and the corresponding comparison results are provided in Table 5. Regarding the models ℳ1 and ℳ2, the same facility locations are obtained because the fixed cost does not influence the determination of the facility locations. When the facility capacities are considered, the globally optimal facility locations in model ℳ3 can still be obtained, which is different from that in model ℳ2. Besides, a higher total cost is obtained in model ℳ3 compared with model ℳ2, which verifies that the facility capacities affect the determination of facility locations.

Comparison with Previous Methods
When the potential facility candidates are not given, the general strategy to determine the set of facility locations is the application of clustering-based algorithms. Among them, the K-means clustering algorithm is one of the most popular methods in handling continuous multi-facility location problems. However, the minimizer of the total weighted Manhattan distance is the weighted median rather than the weighted mean in each cluster. In this sense, the weighted median is used to update the cluster center in each of the iterations in the clustering-based algorithms [18]. In addition, Gao, Nayeem [40] also developed a clustering-based genetic algorithm for the continuous multi-facility location problem considering the Manhattan distance.
To illustrate the outperformance of the proposed approach, different methods are compared for Instances 24-27 in handling the continuous multi-facility location problems. Given different numbers of facilities, the comparison results including objective function values and the facility locations are presented in Figure 13 and Tables 6-9, respectively. Since different randomly generated initial solutions result in different clustering results, we repeat each simulation five times and pick up the best one among them. As shown in Figure 13, the performance of minimizing the total transportation cost using different meth-ods is presented. The proposed method outperforms other clustering-based algorithms. Besides, given any number of facilities from 2 to 12, the method proposed by Gao [65] demonstrates better performance in reducing the total transportation cost compared with the k-means method and clustering-based genetic algorithm, since the weighted median is used to update the cluster center. Besides, from Tables 6-9, the facility locations obtained using the method proposed by Gao [65] are closer to the globally optimal solution compared with the other two existing methods. In this sense, the method proposed by Gao [65] is more appropriate to determine the optimal solution when the problem size is quite large. As a consequence, regarding the Manhattan distance, there is definitely strong evidence that the proposed exact method can determine the globally optimal facility locations for the continuous multi-facility location problems. Moreover, the application of the weighted median in the clustering-based methods is also reasonable when the Manhattan distance is considered.
it to high-dimension multi-facility problems and transferred the continuous-space location problem to the discrete-space location problem. As a consequence, three mathematical models were proposed to formulate different types of continuous-space multi-facility location problem and the corresponding globally optimal solutions are obtained, which is different from the previous clustering-based methods. To illustrate the effectiveness and efficiency of the proposed exact method, several illustrative instances from a benchmark were used to compare the performances of different methods. The comparison results show the superiority of the proposed exact method in the decision-making of the continuousspace multi-facility location problems.
Despite the above contributions and insights, this work still has several limitations. This study only investigated a continuous-space multi-facility location problem based on the Manhattan distance. Moreover, it is difficult to determine the global optimum when the problem size is quite large. Accordingly, some meaningful problems of interest can be explored in more depth in future. It is important to find the candidate locations of the facilities for the Euclidean distance. Besides, we also need to develop a heuristic algorithm to determine the facility locations based on all of the candidate locations when the problem size is large. Furthermore, the maximum distance (i.e., blast radius) from the facility to the demand point will be considered in the continuous-space multi-facility location problem.
Author Contributions: Conceptualization, X.G. and C.P.; software, X.G. and X.C.; writing-original draft preparation, X.G. and E.X.; writing-review, editing, and revising, X.G., C.P., G.H. and D.Z.; funding acquisition, C.P. All authors have read and agreed to the published version of the manuscript.