ADMM-type methods for generalized multi-facility Weber problem

The well-known multi-facility Weber problem (MFWP) is one of fundamental models in facility location. With the aim of enhancing the practical applicability of MFWP, this paper considers a generalized multi-facility Weber problem (GMFWP), where the gauge is used to measure distances and the locational constraints are imposed to new facilities. This paper focuses on developing efficient numerical methods based on alternating direction method of multipliers (ADMM) to solve GMFWP. Specifically, GMFWP is equivalently reformulated into a minmax problem with special structure and then some ADMM-type methods are proposed for its primal problem. Global convergence of proposed methods for GMFWP is established under mild assumptions. Preliminary numerical results are reported to verify the effectiveness of proposed methods.


1.
Introduction. Let a 1 , · · · , a n be the locations of n customers in the plane. The well-known multi-facility Weber problem (MFWP) is to find the locations of m new facilities x 1 , · · · , x n such that the sum of transportation costs between n customers and m facilities and the costs between m facilities is minimized. More specifically, the mathematical model of MFWP is as follows: where w ij > 0 is the weight associated with the unit transportation cost between x i and a j , v il > 0 is the weight associated with the unit cost between x i and x l , and · is the Euclidean distance.
As is known, MFWP is one of the fundamental models in the field of facility location. Due to the wide applications in operations research, marketing, urban planing etc., MFWP has attracted wide attentions and many authors are devoted to both promoting its theoretical development and presenting effective numerical algorithms, see, e.g. [7,8,24,29]. The first attractive contribution to this aspect, denoted by Miehle algorithm (MA), is due to [24], where an extension of Weiszfeld procedure [36] was proposed for MFWP. Ostresh [26] proved later that the Miehle algorithm is a descent method. Rosen and Xue [31], however, constructed counterexamples showing that Miehle algorithm may stop or converge to one nonoptimal point.
How to improve the performance of MA and make it computationally preferable thus become the main challenges in the study of MFWP. The famous hyperboloid approximation procedure (HAP) is proposed in this effort, which is perhaps the most standard method for MFWP at present. The idea of HAP is to introduce a small positive scalar ε to the actual objective function to get a smooth approximating problem which is then solved by Miehle algorithm, see e.g. [9,30]. The global convergence of HAP can be found in [30]. Recall that for HAP the larger the value of ε, the faster the convergence to optimal value of approximating function but the worse the accuracy of approximation to the actual nondifferentiable function. So HAP can be improved by successively reducing the value of ε, which is called SHAP. For more details about HAP please refer to [6,29].
Note that the distances in MFWP (1) are Euclidean distances, i.e., measured by l 2 -norm. As in MFWP (1), the distances of facility location are usually measured by norms such as l p -norms and block norms in the literature. [23] discussed the approximations about the weighted l p -norms based on statistical evidence, and [35] investigated the use of block norms to obtain good approximation of actual distances. The symmetry property of the norm assures that the distance from one point to another is always equal to the distance back. Nevertheless, [38] has highlighted the fact that the symmetry of the distance is violated in numerous real situations. The examples include transportation in rush-hour traffic, flight in the presence of wind, navigation in the presence of currents, transportation on an inclined terrain, etc. Ever since the pioneering work of [38], the asymmetric distance has attracted the researchers' interest, and some progress has been made in its theoretical and computational aspects recently, see e.g. [2,3,20,28].
It is also clear that the locations of new facilities in MFWP can be chosen anywhere in the plane. However, it is not always this case in practice; that is to say, not every point in the plane is suitable to locate the new facilities due to some actual situations. For example, a fire station can not be built within an area where nobody lives; a hospital or medical center can not be very close to heavily polluting and noisy factories. Therefore, after we consider some relevant practical factors the locations of new facilities are usually restricted in some specific regions such as convex subsets of the plane.
In this paper we consider more general case of MFWP to enhance its practical applicability. Firstly, the gauge rather than the widely-used symmetric norm is used to measure distance between points. With the usage of gauge, both symmetric and asymmetric distances are allowed to be considered in our model. Secondly, the locational constraints are imposed for new facilities. By choosing appropriate locational constraints for facilities, both constrained and unconstrained location problems can be investigated by our model.
The rest of this paper is organized as follows. The general MFWP is formulated in Section 2. In Section 3, it is reformulated into a minmax problem and then the corresponding primal problem and dual problem are given. Then in Section 4, some ADMM-type methods are proposed for solving the primal problem. Preliminary numerical results are reported in Section 5 to verify the efficiency of proposed ADMM-type methods. Finally, some conclusions are drawn in the last section.
2. Model formulation. Let B be a compact convex set in R 2 whose interior contains the origin, then the gauge based on B is defined by B is called the unit ball of γ(·) due to This definition of a gauge from a compact convex set was first introduced by Minkowski [25]. The gauge γ(·) satisfies the following properties: for any x and y ∈ R 2 .
It follows from 2) and 3) that any gauge γ(x) is a convex function of x. The distance measuring function can be derived from a gauge by When B is symmetric around the origin, according to the definition of (2) we have It follows that D(x, y) = D(y, x), which means the distance measured by γ(·) is symmetric. On the contrary, when B is not symmetric around the origin, (5) does not hold any more and then the distance measured by γ(·) is asymmetric. Thus, with different compact convex sets used as unit balls, various gauges (symmetric or asymmetric) can be employed to measure distances in location problems, depending on the requirements of practical applications.
With considering the gauge to measure distances and the locational constraints on facilities, the generalized multi-facility Weber problem (abbreviated as GMFWP) is formulated as: where γ(·) is the gauge generated by a given compact convex set B ⊂ R 2 and X i denotes the locational constraint for the i−th facility. For simplification, all X i , i = 1, · · · , m, are assumed to be closed convex subsets of R 2 throughout this paper. Note that X i can also be chosen as R 2 . If all X i are R 2 , the GMFWP becomes an unconstrained location problem.
Remark 1. The GMFWP (7) includes some well-known location problems as its special cases: if m = 1, X = R 2 and γ(·) is the Euclidean distance, (7) is the classical Weber problem; if m = 1 and γ(·) is the Euclidean distance, (7) is a constrained Weber problem; if all X i = R 2 and γ(·) is the Euclidean distance, (7) is the classical multi-facility Weber problem (1).
Though the gauge is used and locational constraints are considered, fortunately we have the following proposition for GMFWP which reveals its convexity.
Proof. Note that all constraints X j , j = 1, · · · , m, are closed convex sets and the function γ(·) is also convex, then it follows immediately that the GMFWP is convex.
While GMFWP (7) is convex, it is obvious that the nondifferentiability will occur similarly as MFWP. This means that both GMFWP and MFWP share the unexpected characteristic of singularity.
3. The primal and dual reformulations of GMFWP. For the gauge γ(·) defined by (2), there always exists a dual gauge γ d (·) defined by Let B d be the unit ball of the dual gauge γ d (·), which is also convex and compact. B d is exactly the polar set of B. It is well known that the dual gauge of γ d (·) is again γ(·), i.e., which can be rewritten as For more details about the gauge and dual gauge, as well as their properties, the readers can be referred to [28]. According to (10), GMFWP (7) is equivalent to the following min-max problem: where each z ij is a vector in B d ij = {ξ ∈ R 2 | γ d (ξ) ≤ w ij } for i = 1, · · · , m and j = 1, · · · , n; each z il is a vector in B d il = {ξ ∈ R 2 | γ d (ξ) ≤ v il } for i = 1, · · · , m and l = 1, · · · , n (l = i).
We define the following matrix . .
which consists of mn × m blocks. Each block of A 1 is either 2 × 2 identity matrix (I 2 ) or 2 × 2 zero matrix (blank block), and each block column includes n identity matrices while the others are all zero matrices. We also define m matrices N t , t = 1, · · · , m, as Similar to A 1 , each block of N t is also either 2 × 2 identity matrix or 2 × 2 zero matrix. Based on these N t , t = 1, · · · , m, we define the following matrix By introducing the following B d , z, X and x , · · · , z T 1 n · · · , z T m1 , · · · , z T m n , z T 12 , · · · , z T 1 m , · · · , z T m1 , · · · , z T m m−1 ) T , and introducing the matrix A and vector b T , · · · , a n T , · · · , a 1 T , · · · , a n T , 0, · · · , 0) T , where each zero block is a 1 × 2 zero matrix and there are totally m(m − 1) zero blocks in b, (11) can be written in a compact form as follows: We further rewrite (15) as where δ X and δ B d are the indicator functions of X and B d , respectively.
where δ * X is the conjugate of δ X . According to (17) and (18), it is clear that the min-max problem (16) is a primal-dual reformulation of the following primal problem: (P) min or of the following dual problem: It thus follows that we can solve the primal problem (P) or the dual problem (D) instead of solving (16) directly. When we compare the scales of (P) and (D), it is found that the number of variables of (P) is much smaller than that of (D) (the number of variables of (P) is m and the number of variables of (D) is m(m+n−1)). This observation encourages us to solve the GMFWP (7) by solving the primal problem (P) in an efficient way.
Remark 2. The l p -norm (especially l 1 , l 2 , l ∞ -norm) and block norm are frequently used to measure distances in the literature, see, e.g., [23,35]. It is noted that the gauge used in this paper is an extension of norm. When the gauge γ(·) is chosen as l 1 , l 2 , l ∞ -norm, the dual γ d (·) will be l ∞ , l 2 , l 1 -norm, respectively; and when γ(·) is a block norm, its dual is also a block norm.
Remark 3. The paper [19] considered similar generalized MFWP. However, in [19] a variational inequality (VI) approach was developed and new projection-type methods under VI were proposed for solving GMFWP. This paper reformulates GMFWP as a primal problem of the equivalent min-max problem and will propose ADMM-type methods to solve it. The numerical experiments to be reported later verify the computational preference of ADMM-type methods to projection-type methods.
Remark 4. A recent paper [18] also applied ADMM technique to solve facility location problem. In fact, these two papers solve two different facility location problems. The paper [18] considered the nonconvex constrained multi-source Weber problem, where ADMM was employed to solve the subproblem of the proposed location-allocation algorithm. Instead, this paper considers a generalized MFWP, which is convex itself. Thus we need to modify ADMM in a different way here.
4. ADMM-type methods for the primal problem (19). This section aims to develop efficient numerical methods for the primal problem (P) by exploiting its favourable structure. With introducing an auxiliary variable y, the primal problem (P) can be rewritten as Let L β (x, y, λ) be the augmented Lagrangian function of (21), which has the form of where λ is the Lagrange multiplier (dual variable) and β > 0 is a penalty parameter for the linear constraint.
Revisiting the formulation of (21), we can easily see that its objective function is separated into two individual functions without coupled variables. In other words, (P) is formulated as a linearly constrained two-block separable convex minimization problem. A large number of separable convex problems arises from various practical areas such as image processing and statistical learning, and thus this type of problems has been well studied in the past several decades. It is evident that the favorable separable structure is very beneficial for developing efficient numerical methods. To exploit the separable structure, operator splitting methods have been well discussed in the literature. The alternating direction method of multipliers (ADMM) proposed originally in [11,13] is one of the most popular operator splitting methods, which is closely related to the well-known Douglas-Rachford splitting method in partial differential equation. For the comprehensive review of ADMM and its novel applications in various areas, the reader can be referred to [1,5,14,15,21,32,37] and references cited therein.
This paper concentrates on presenting some efficient ADMM-type numerical methods to solve the two-block separable problem (21) and thus solve the GM-FWP (7). Applying the ADMM to (21), it will result in the following ADMM-type iterative scheme.
In order to accelerate the convergence of ADMM, the penalty parameter β can also vary self-adaptively according to certain strategies, see, e.g., [5]. For the simplification of our discussion, the parameter β is set to be a constant positive number throughout this paper. Next we will discuss how to solve the subproblems of ADMM-type method (23).
I) The variable x k+1 . According to (22) and (23a), x k+1 can be obtained by solving The x-subproblem can not be easily solved because of general structure of the matrix A. We know that the efficiency of ADMM-type methods heavily depends on how fast the subproblems are solved. To overcome this difficulty, we apply the linearization technique for x-subproblem. We linearize the quadratic term β (24) and add a proximal term r 2 x − x k 2 to approximate the objective function as follows The linearized x-subproblem can be also seen as a special application of PPA because of the additional proximal term r 2 x − x k 2 . It is obvious that (25) is quite simple enough to have the following closed-form solution where P X (v) is the projection of v onto the set X, i.e., Thus, no inner iteration is required for the x-subproblem by using the linearization strategy.
II) The variable y k+1 . Let p = Ax k+1 − b − 1 β λ k , then it is easy to see that y k+1 is the solution to the following problem Using Moreau's decomposition, we can characterize the closed-form solution to y-subproblem (27) by the following proposition.
Proposition 2. For β > 0 and the vector p defined in (27), the solution to the y-subproblem is given as Proof. Let ϕ : R 2m(m+n−1) → R be a function and Prox dϕ be the proximity operator of function ϕ, i.e., Thus, according to (27) the solution to y-subproblem is It follows from Moreau's decomposition that where the last equation is from the proposition of conjugate that δ * * = δ. For the second term Prox βδ (βp) of (30), according to the definition of proximity operator it satisfies Substituting (31) into (30), we obtain Combining it with (29), this proposition is derived directly.
Based on the previous discussion on the iterative scheme and the corresponding subproblems, now we are ready to present our first ADMM-type method for (19), denoted as ADMM-a, as follows. Algorithm 1. ADMM-a for (19) Given the parameter β > 0, r > 0 and initial iterate w 0 = (x 0 , y 0 , λ 0 ). For k = 0, 1, 2 · · · , do: As we have explained, the proposed ADMM-a is in fact a linearized ADMM method for (19), in which the x-subproblem is linearized and thus its closed-form solution can be derived easily. The linearized ADMM has been well studied in the past several years, see, e.g., [15,22,39,34]. Our proposed ADMM-a is a special scenario of that of the generic linearized ADMM and thus its convergence can be found in the literature. Let {w k } be the sequence generated by ADMM-a with w k = (x k , y k , λ k ), and we state the convergence result for (19) as follows with omitting the proof.
To ensure the convergence of ADMM-a, it follows from Theorem 4.1 that the parameter β and r should satisfy rI − βA T A 0. Note that where λ max (A T A) is the largest eigenvalue of A T A. It means that r should not be too small. On the other hand, however, r can not be too large; otherwise, the linearized x-subproblem (25) will be far from the original x-subproblem (23a), which may result in a slow convergence of ADMM-a. Then, it leads to a natural concern about how to choose the value of the parameter r. In the following, a self-adaptive adjustment of the value of r is presented which can alleviate the restriction of (35) for r to some extent. Based on this self-adaptive manner, we propose the second ADMM-type method (denoted as ADMM-b) for (19). It consists of two steps: the prediction step and the correction step, where the prediction step produces a tentative pointw k and the correction step generates the next iterate w k+1 . The iterative framework of ADMM-b is elaborated in the following. (19) Given the parameter β > 0 and initial iterate w 0 = (x 0 , y 0 , λ 0 ). For given w k = (x k , y k , λ k ) and a small constant r > 0, w k+1 = (x k+1 , y k+1 , λ k+1 ) is generated via the following scheme.

Algorithm 2. ADMM-b for
Step 1. Prediction step. Step Step 1.2. If the following inequality is not satisfied, then set r = 2r and go to Step 1.1.
Remark 6. It is obvious that the prediction step of ADMM-b is as the same as that of ADMM-a. In other words, ADMM-a provides a tentative point for ADMMb in the prediction step, and then based on such a point the new iterate w k+1 is generated in the correction step. The workload of correction step can be almost negligible since the term A T A(x k −x k ) has been already computed when we check the inequality (37) in Step 1.2.
Remark 7. The parameters β and r are required to satisfy the condition (37) in ADMM-b. Thus, for given β, instead of using a fixed value the parameter r can be adjusted in a self-adaptive manner based on (37). When we compare the two conditions (35) and (37), it is easy to see that (35) implies (37). Thus, the restriction on β and r of ADMM-b is more relaxed than ADMM-a with the adoption of selfadaptive strategy.
ADMM-b is an application of linearization technique to our problem (21), and the convergence of ADMM-b can be obtained from [15]. Analogous to ADMM-a, we omit the convergence proof of ADMM-b in detail and just give its convergence result.
Recall that the two ADMM-type methods proposed above are linearized ADMM. The linearization technique is quite important for ADMM when one or more subproblems are hard to solve. With the adoption of linearization technique, the difficult subproblems may be simple enough to have closed-form solutions or can be easily solved up to a high precision, as ADMM-a and ADMM-b have shown. On the other hand, however, the following two points should be concerned for the linearized ADMM. The first one is that the linearized subproblem are not the original ADMM subproblem, which will affect the performance of linearized ADMM. The second point is how to choose an appropriate value for r in x-subproblem (25) so as to enhance the convergence. While the inequalities (35) and (37) ensure the convergence of ADMM-a and ADMM-b, respectively, we really don't know which value of r can result in fast convergence of two methods.
The mentioned two points thus inspire us to consider whether an ADMM-type method without the adoption of linearization technique can be developed for solving (21). Thanks the favorable structure of (21), we can present our third ADMMtype method which does not adopt linearization technique for subproblems in the following.
First, similar to the variable y in (21) we introduce another auxiliary variable z to (21), then its formulation can be rewritten as It is clear that the formulation (39) is a three-block separable convex minimization problem whose variables are composed of x, y and z.
The augmented Lagrangian function of (39) is written as According to the augmented Lagrangian function (40), we can apply the ADMM for three-block separable problem to (39) which results in the following iteration scheme.
In the following we will discuss how to solve the subproblems of (41) in detail.

I)
The variable x k+1 . It follows from (40) and (41a) that x k+1 is the solution of the following optimization problem Now this x-subproblem is quite simple since the coefficient of x in the quadratic term is the identity matrix, and thus it has the closed-form solution as follows II) The variable y k+1 . It is easy to see that y k+1 satisfies where q = Az k − b − 1 β λ k 1 . Analogy to Proposition 2, using Moreau's decomposition we can obtain the closed-form solution of y-subproblem as follows III) The variable z k+1 . We can get z k+1 via solving This is an unconstrained optimization problem with regard to z and its optimal solution can be obtained by solving g(z k+1 ) = 0, where g(·) is the gradient function of (44), i.e., Therefore, the solution to z-subproblem is given as where the vector u = A T λ k 1 − λ k 2 + βA T (y k+1 + b) + βx k+1 . IV) The Lagrange multipliers λ k+1 1 and λ k+1 2 . The two multipliers are updated according to (41d) and (41e), respectively, which are trivial computation. Now we are at the stage to present the following ADMM-type method, denoted as ADMM-c, to solve (19).

Algorithm 3. ADMM-c for (19)
Given the parameter β > 0 and initial iterate w 0 = (x 0 , y 0 , z 0 , λ 0 ). For k = 0, 1, 2 · · · , do: With introducing a new auxiliary variable z to the model (21), all subproblems of ADMM-c have closed-form solutions, as indicated by (45a)-(45e). Recall that the xsubproblem (24) is hard to solve due to the existence of the general matrix A. With the variable z introduced, however, the matrix A has moved from the x-subproblem (42) to the z-subproblem (44). While ADMM-c does not adopt the linearization technique as in ADMM-a and ADMM-b, both two subproblems (42) and (44) are now easy to solve. This idea contributes to the improvement and effectiveness of ADMM-c for GMFWP, which will be verified by the numerical results later.
One essential concern with regard to ADMM-c is its convergence. The proposed ADMM-c is the application of classical ADMM for three-block separable convex problem (39). A recent paper [4], however, has revealed that the convergence of classical ADMM for the generic three-block separable convex problem is not guaranteed, even though its numerical efficiency has been verified empirically by many satisfactory applications (e.g. [27,33]). Fortunately, due to the favorable structure of (39), the proposed ADMM-c has global convergence, i.e., for arbitrary initial iterate w 0 the iterate sequence {w k } generated by ADMM-c is convergent to the optimal solution to (39). To see this, let us have a closer look at the iterative scheme of x and y. First, we couple x with y to form the vector (x, y). Based on the augmented Lagrangian function (40), we apply the ADMM for two-block convex problem, then the iterative scheme turns to be as follows.
By dropping the terms irrelevant to x and y, the subproblem (46a) can be written as where q is defined as in (43). It means that for (x k+1 , y k+1 ) in (46a) the variables x k+1 and y k+1 are obtained by solving the x-subproblem and y-subproblem of ADMM-c independently. It follows that ADMM-c can be regarded as the classical ADMM for two-block separable convex problem, and thus its global convergence can be established as follows.
Proof. Let u = (x T , y T ) T and u * = (x * T , y * T ) T . It follows from Theorem 1 in [18] that we have where and (x * , y * , z * ) is the optimal solution to (39). According to (47) Therefore, the sequence {w k } converges to w * = (x * , y * , z * , λ * ). Since (x * , y * , z * ) is the optimal solution to (39), it is clear that x * is the optimal solution to (19). The proof is complete.
Recall that as the linearization technique is adopted, the parameters β and r of ADMM-a and ADMM-b should satisfy some particular conditions to ensure the global convergence, e.g., rI − βA T A 0 for ADMM-a. Compared with ADMM-a and ADMM-b, the global convergence condition of ADMM-c is more loose than two previous ADMM-type methods.
5. Numerical results. This section reports some numerical results to verify the computational efficiency of ADMM-type methods. We will implement the proposed ADMM-a, ADMM-b and ADMM-c to solve a number of GMFWPs including unconstrained problems and constrained problems and also including problems with symmetric distances and asymmetric distances. To demonstrate the efficiency of ADMM-type methods, in the first subsection we compare ADMM-a, ADMM-b and ADMM-c with the well-known methods in facility location for MFWP and in the second subsection we compare them with some methods in the variational inequality area. The parameters of ADMM-type methods are given as follows: β = 0.01 and r = β ·λ max (A T A)+1 for ADMM-a; β = 0.01 and r = β ·λ(A T A)+10 for ADMM-b, where λ(·) denotes the average of eigenvalues; β 1 = 0.05 (for the first constraint) and β 2 = 50 (for the second constraint) for ADMM-c. All the programming codes were written by MATLAB 9.0 and were run on a ASUS laptop (Intel Core i7-6700HQ 2.60GHz). 5.1. ADMM-type methods for MFWP. Miehle algorithm (MA) and hyperboloid approximation procedure (HAP) are the most important methods for solving MFWP. This paper proposes ADMM-type methods to solve GMFWP. Note that MFWP, as a classical model in facility location, is a special case of GMFWP.
Hence, ADMM-type methods are applicable to solving MFWP. In the following, we will compare ADMM-a, ADMM-b and ADMM-c with MA and HAP for solving MFWP and show that even for such a special problem the proposed ADMM-type methods are efficient.
First we compare ADMM-a, ADMM-b and ADMM-c with MA. The following two examples were constructed in [31].
It has been shown in [31] that the iterate sequences generated by MA stop at a nonoptimal point of E1 and converge to a nonoptimal point of E2. This means that the convergence of MA to the optimal solution is not guaranteed.
We further construct the following two examples.
E1 and E2 : E1 and E2 with the initial iterate x 0 1 and x 0 2 chosen as the locations of two customers.
As the location of customers are chosen as initial iterates in E1 and E2 , the undesirable singularity will occur. The two examples are generated to show the difference of the performance of ADMM-type methods with MA in case of singularity.
Note that the location of customers, the weights and the initial iterates are particularly chosen in E1, E2, E1 and E2 . For well justification of the efficiency of ADMM-type methods, we construct the following randomly generated problems.
Note that MA is a Jacobi-like method because the locations of new facilities change at the same time. An alternative method potentially accelerates the convergence of MA by applying new information of facility locations once it becomes available. This leads to the so-called Gauss-Seidel-like Miehle algorithm (GMA).
We solve the five examples by MA, GMA, ADMM-a, ADMM-b and ADMM-c and the numerical results are reported in Table 1, where "CPU" means the computing time, "Iter." gives the number of iterations and "Obj." reports the objective functional values. The stopping criterion is chosen as Each example is solved by five methods with the same initial iterate x 0 1 and x 0 2 . The initial y 0 , z 0 and λ 0 in ADMM-a, ADMM-b and ADMM-c are zero. For E1 , E2 and E3, we test each example for one hundred scenes with each scene solved by five methods and report the average results.
For E1, MA achieves a point, (1.50591 0.00000) and (1.50591 0.00000), in the first iteration and then stops at it from the second iteration. However, it is not the optimal solution of E1, which is (0,0) and (0,0). When GMA is applied to solve E1, the point obtained by GMA, (1.65749,0) and (1.65749,0), is also not the optimal one. In fact, this point is even worse than that from MA. While ADMM-a, ADMMb and ADMM-c need more iterations for E1, they can get their optimal solution and provide the optimal objective functional value. For E2, with the adoption of the Gauss-Seidel idea GMA needs much less number of iterations than MA. When we compare the points obtained by five methods, ADMM-a, ADMM-b and ADMM-c provide the optimal solution but the other two methods do not.
For E1 and E2 , the undesirable singularity occurs. Thus, MA and GMA can not proceed the iterations and then their numerical results are labeled as "/" in  Table 1. ADMM-a, ADMM-b and ADMM-c, however, are able to deal with the singularity and in such case ADMM-a, ADMM-b and ADMM-c can still get the optimal solution, which is consistent with the analytic results. For the randomly generated E3 , it is found that the solution obtained by GMA is the worst among five methods. MA and GMA need more than 15 times computing time than proposed ADMM-type methods, but the points are still not the optimal one of E3 . Compared to MA and GMA, ADMM-a, ADMM-b and ADMM-c can get the optimal solution of E3 .
Next we will compare ADMM-a, ADMM-b and ADMM-c with HAP. An alternative of HAP using the new available information was proposed in [30]. This alternative is called Gauss-Seidel-like HAP (GHAP). The parameter ε is chosen as 5, 0.5, 0.05, respectively, and we report the numerical results of HAP and GHAP under different values of ε. SHAP successively reduces the value of ε during its implementation to accelerate the convergence of HAP under small value of ε. This technique is also applicable to GHAP and GHAP with such technique is called SGHAP in this paper. For SHAP and SGHAP, the ε is reduced from 5 to 0.05 with its value multiplied by 0.1 each time. For convenience, we call HAP, GHAP, SHAP and SGHAP as HAP-type methods.
We test MFWP with n = 500 and m = 2, 4, 6, 8, 10. All customers are randomly generated in [−250, 250] 2 ; the weights between customers and facilities are randomly chosen in [1,5] and the weights between facilities are chosen as 100. The stopping criterion is set as x k+1 − x k < 10 −4 . The initial iterate is randomly generated in [−250, 250] 2 . For each pair (n, m), we solve 100 randomly generated MFWPs by HAP, GHAP, SHAP, SGHAP, ADMM-a, ADMM-b and ADMM-c using the same randomly generated customers and initial iterates. The numerical results are reported in Table 2, where "CPU" and "Iter." are as the same as those in Table 1 and "Dobj." reports the difference of objective functional values between ADMM-a and other methods. Actually the value of "Dobj." is the functional value of other methods minus that of ADMM-a. Therefore, if this value is positive, the point obtained by other method is worse than that by ADMM-a; otherwise, it is better. It is clear that the smaller the value of "Dobj.", the better the point obtained by the method. According to the last three columns in Table 2, it can be seen that the performances of ADMM-a, ADMM-b and ADMM-c are quite similar, and ADMM-c outperforms ADMM-a and ADMM-b a little.
When we compare ADMM-a, ADMM-b and ADMM-c with four HAP-type methods, it is easy to observe that ADMM-a, ADMM-b and ADMM-c are preferable to HAP-type methods especially for large-scale problems. Firstly, we compare them from the quality of obtained points. We have known that the points obtained by ADMM-a, ADMM-b and ADMM-c are almost the optimal solution of MFWP. However, the points obtained by HAP-type methods are clearly not the optimal solution, which is indicated by the positive values of "Dobj.". Secondly, we compare seven methods from the computing workload. It is clear that the workloads of ADMM-a, ADMM-b and ADMM-c are less than HAP-type methods especially for large-scale MFWP. With the value of m increasing, the increment speed of workload of HAP-type methods is much faster than ADMM-type methods. In addition, we also expect to decrease ε of HAP-type methods to improve the quality of obtained points. However, it is found during our experiments that though the quality of obtained points is improved, the computing time and number of iterations increase drastically with smaller ε.
Thus, with the comparison of ADMM-a, ADMM-b and ADMM-c to some wellknown methods in facility location area, their efficiency is verified. 5.2. ADMM-type methods for GMFWP. In this subsection we will apply ADMM-type methods to solve GMFWP (7). Recently in the paper [19], GMFWP is reformulated into a linear variational inequality (LVI). Note that there exist some methods for LVI, which are also applicable to solving GMFWP. To verify the efficiency of ADMM-a, ADMM-b and ADMM-c, we compare them with some methods for LVI to solve GMFWP. The first method to be compared is the projection and contraction (PC) method presented in [17] : where γ is a positive constant. The second one is the customized proximal point algorithm (PPA) proposed in [16]: where r, s and γ are positive constants. The convergence of PPA under the condition of rs > A T A and γ ∈ (0, 2) can be found in [16]. The paper [19] proposed a new projection-type method under the LVI framework for GMFWP, whose iterative scheme is as follows.
where W is an invertible matrix and r, s and γ are positive constants. Let the matrix Q defined as Q = rI 0 A sI .
The parameters of the above three methods are given as follows.
• PC method: γ = 1.9; • PPA: r = p · t A T A , s = 2 p · t * A T A , t = 1.01, p = 2, γ = 1.2; • Projection-type method: We first test GMFWP under Euclidean distances with n = 50, 100, 500, 1000, 2000 and m = 2, 4, 6, 8, 10. All customers are randomly generated in [−250, 250] 2 ; the weights between customers and facilities are randomly chosen in [1,3] and the weights between facilities are chosen in [1,5]. The locational constraints are set as , 250] 2 and r i is randomly chosen in [10,15], i = 1, · · · , m. The initial iterate is randomly generated in [−250, 250] 2 and the stopping criterion is given as For each given pair (n, m), we test the GMFWPs for 50 scenes with each scene solved by six methods using the same randomly generated initial iterates. Table 3 reports the average computing time in seconds (CPU) and the average number of iterations (Iter.) of six methods to solve the GMFWPs under Euclidean distances (symmetric distances). According to Table 3, six methods can solve the randomly generated GMFWPs under Euclidean distance and with locational constraints efficiently. It is shown that the proposed ADMM-type methods are preferable to the PC method, PPA and projection-type method for solving the GMFWPs of all scales. For large-scale GMFWP, the workload of ADMM-type methods is only about 10% that of other methods. In addition, it is also found that ADMM-c is preferable to ADMM-a and ADMM-b for large-scale problems.
Since ADMM-type methods can solve the GMFWP whose distances are measured by asymmetric gauge, we test some randomly generated GMFWPs under gauge with different scales. The asymmetric gauge γ(·) is generated with the unit ball set as (x + 2) 2 + 2y 2 ≤ 8.
To verify the efficiency of ADMM-type methods for GMFWP under gauge, we also compare ADMM-a, ADMM-b and ADMM-c with the other three methods simultaneously. The performance of six methods is reported in Table 4.
It is easy to observe that the proposed ADMM-a, ADMM-b and ADMM-c are efficient for GMFWP under gauge and they outperform the PC method, PPA and projection-type method for the randomly generated asymmetric problem.
When we compare the performance of three ADMM-type methods for largescale problems in Table 3 and Table 4, we found that ADMM-c is a little better than the other two methods in Table 3 but much better than them in Table 4. This implies that ADMM-c is robust to the measuring distance function (norm or gauge). Moreover, it follows from Table 3-4 that the computing time and number of iterations of three ADMM-type methods is not so sensitive to the size of GMFWP. All of these desirable properties contribute to making the proposed ADMM-type methods especially ADMM-c promising in practice. According to Table 3-4, we also found that due to the asymmetric property of gauge the computational load of six methods for GMFWP under gauge are much greater than those for GMFWP under Euclidean distances. 6. Conclusion. A generalized multi-facility Weber problem is investigated in this paper. This problem considers the locational constraints on facilities and employs the gauge to measure the distances. With the consideration of the primal problem of reformulated minmax problem, some ADMM-type methods are proposed for solving the generalized problem. The convergence of proposed methods is established under mild assumption. Preliminary numerical results justify the preferable merits of proposed ADMM-type methods.
How to solve other problems arising in facility location computationally by the ADMM approach deserves further and more extensive investigation in our future research.