A variational inequality approach for constrained multifacility Weber problem under gauge

The classical multifacility Weber problem (MFWP) is one of the most important models in facility location. This paper considers more general and practical case of MFWP called constrained multifacility Weber problem (CMFWP), in which the gauge is used to measure distances and locational constraints are imposed to facilities. In particular, we develop a variational inequality approach for solving it. The CMFWP is reformulated into a linear variational inequality, whose special structures lead to new projection-type methods. Global convergence of the projection-type methods is proved under mild assumptions. Some preliminary numerical results are reported which verify the effectiveness of proposed methods.


1.
Introduction. The issue of the location of new facilities in modern logistics system is a strategic sense. The past decades have witnessed a fruitful development in the field of location theory as well as its wide applications in diverse fields. For a comprehensive review, see e.g. [6,18].
The classical multifacility Weber problem (MFWP) is to find the locations of m new facilities in order to minimize the sum of transportation costs between n customers and m facilities and transportation costs between m facilities. More specifically, the mathematical model of MFWP is the following:

MFWP: min
where 1) a j ∈ R 2 is the location of the j − th customer, j = 1, 2, · · · , n; 2) x i ∈ R 2 is the location of the i − th facility to be determined, i = 1, 2, · · · , m; 3) w ij > 0 is the weight associated with the transportation cost between the i − th facility and the j − th customer; 4) v il > 0 is the weight associated with the transportation cost between the i−th facility and the l − th facility; 5) · is the Euclidean distance.
MFWP is one of the most important models in the area of facility location. Due to its 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. [5,6,18,20,25].
In the literature, the distances of facility location are usually measured by norms such as l p -norms and block norms. [16,17] discussed the approximations about the weighted l p -norms based on statistical evidence. [30,31] 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, [33] has highlighted the fact that in numerous real situations the symmetry of the distance is violated. 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. In spite of the pioneering work of [33], the asymmetric distance has started to attract the researchers' interest until the recent several decades, and some progress has been made in both its theoretical and computational aspects recently, see e.g. [2,3,7,15,22,24].
As a matter of fact, not every point on the plane is possible to locate the new facilities due to some actual situations. For example, we can't build a hospital or medical center where there is a mountain or a lake. In addition, the location of a facility in health care system can't be very close to heavily polluting and noisy factories. Therefore, the locations of new facilities are usually restricted to be in some specific regions such as convex subsets of the plane after considering some relevant practical factors .
In this paper we consider the more general case of MFWP with the aim of enhancing its practical applicability. Firstly, the gauge rather than the widely-used symmetric norm is used to measure distance between points. By using gauge as the distance measuring function, both symmetric and asymmetric distances can be discussed in our model. Secondly, we introduce locational constraints on new facilities. By choosing appropriate locational constraints for facilities, e.g. the plane R 2 or the subsets of R 2 , both the unconstrained and constrained problems can be taken into consideration in our model.
The rest of this paper is organized as follows. The formulation of general problem is given in Section 2. In Section 3, it is solved by a variational inequality approach.
Specifically, the relationship between the problem and a linear variational inequality (LVI) is built and new projection-type methods are proposed. Then the convergence properties are proved under some mild conditions. Preliminary numerical results are reported in Section 4 to verify the efficiency of proposed projection-type methods. Finally, some conclusions are drawn in the last section.
2. Model formulation. Let B be a compact convex set in R 2 and the interior of B contain the origin, then the gauge of 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 [21]. 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 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 for measuring distances and the locational constraints on facilities, the constrained multifacility Weber problem (abbreviated as CMFWP) is formulated as follows: where 1) a j , x i , w ij and v il are defined as the same as those in (1); 2) X i denotes the locational constraints for the i − th facility and we assume X i , i = 1, · · · , m, are closed convex subsets of R 2 ; 3) γ(·) is the gauge generated by a given compact convex set B ⊂ R 2 . Remark 1. More explanations are required for the CMFWP (7). Firstly, the locational constraint X i can also be chosen as R 2 . If all X i , i = 1, · · · , m, are R 2 , the CMFWP (7) is a unconstrained location problem; otherwise, (7) is a constrained problem. Thus, both constrained and unconstrained cases can be considered in the CMFWP model. Secondly, the CMFWP (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 multifacility Weber problem (1).
Though the gauge is used and locational constraints are considered, fortunately we have the following proposition for CMFWP which reveals its convexity.
Proof. Note that all constraints X j , j = 1, · · · , m, are closed convex sets and γ(·) is a convex function, it follows that the CMFWP is a convex problem.
While CMFWP (7) is convex, the nondifferentiability occurs, e.g., when a new facility coincides with another new facility or one customer, which is usually called singularity. This unexpected characteristic of singularity is also shared by MFWP. As a classical and important facility location problem, MFWP has been well studied and many approaches have been suggested for solving it. The first attractive contribution to this aspect, denoted by Miehle algorithm (MA), is due to [20], where an extension of Weiszfeld procedure [32] was proposed for MFWP. Ostresh [23] proved later that the Miehle algorithm is a descent method. Rosen and Xue [26], however, constructed counterexamples showing that Miehle algorithm may stop or converge to one nonoptimal point. To deal with the undesired singularity, another well-known method, hyperboloid approximation procedure (HAP), which is perhaps the most standard method for MFWP, was proposed. The idea of HAP is to introduce a small positive number ε to the actual objective function to get a smooth approximating problem which is then solved by Miehle algorithm, see e.g. [9,27]. The global convergence of HAP can be found in Rosen and Xue [27]. 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. Consequently, in the implementation of HAP a large value of ε can be used initially and then the obtained solution is used as the initial iterate for HAP using a smaller value of ε. The process continues by successively reducing the value of ε until no significant improvement is possible. HAP with successively reducing the value of ε is called SHAP throughout this paper. For more details about HAP please refer to [4,25].
In this paper, a variational inequality (VI) approach, different from the existing approaches for MFWP, is developed for solving CMFWP. Specifically, the relationship of CMFWP and a linear variational inequality (LVI) is first built to show their equivalence, and then numerical methods are proposed to solve the related LVI and thus solve CMFWP. The advantages of the variational inequality approach for CMFWP include the following aspects. Firstly, the variational inequality approach provides a new research framework for CMFWP since the topic of variational inequality has been a mature research area (see, e.g., [10,11]) and a large number of efficient methods (e.g. the projection-type methods [12,28] and the proximal point methods [1,19]) have been proposed for solving VI. Secondly, the related LVI problem usually has a special structure and the variational inequality approach can take advantages of such structure to develop new efficient VI-based methods under the algorithmic framework of VI. Thirdly, the global convergence of new methods can be proved under the analytic framework of VI. Fourthly, the variational inequality approach is able to provide both the primal and dual solutions of CMFWP simultaneously. Fifthly and the most importantly, the variational inequality approach can avoid the undesirable singularity of nondifferential objective function and get the optimal solution of CMFWP.
In the spirit of the framework of variational inequality approach, new projectiontype methods are proposed in this paper to solve the related LVI. The iterative scheme of our projection-type methods is simple and easy to implement in computational aspects. The global convergence of projection-type methods is proved. Preliminary numerical experiments compare the projection-type methods with the well-known methods in facility location for MFWP (1) and numerical methods in variational inequality topic applicable to CMFWP (7), which verify the effectiveness of these projection-type methods.
3. The variational inequality approach. This section aims at solving CMFWP (7) by a variational inequality approach.
3.1. LVI reformulation of CMFWP. For the gauge γ(·) defined by (2), there 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 and 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 [24]. According to (10), CMFWP (7) is equivalent to the following min-max problem: , · · · , m and l = 1, · · · , n (l = i).
We first 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 matrix while others are all zero. 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 and introducing the matrix A and vector b where each zero block is a 1 × 2 zero matrix and there are totally m(m − 1) zero blocks in b, we can rewrite (11) as a compact form: Let (x * T , z * T ) T ∈ X × B d be a solution of (15), then it follows that thus (x * T , z * T ) T is a solution of the following LVI: A compact form of the above LVI is : where Therefore, the CMFWP (7) is reformulated into the LVI (18)- (19). Note that M in (19) is a skew-symmetric matrix, then it is positive semi-definite, and thus the LVI (18)-(19) is monotone.

Remark 2.
It is worth pointing out that the equivalence between CMFWP (7) and LVI (18)- (19) can also be obtained by the duality theory. In fact, the variable z ij (z il ) and the set B d ij (B d il ) in (14) are respectively the dual vector and dual ball which satisfy . Remark 3. The norms, especially l 1 , l 2 , l ∞ -norm and block norms, are frequently used to measure distances in the literature, see, e.g., [16,30,31]. It is noted that the gauge used in this paper is an extension of norms. When the gauge γ(·) is chosen as l 1 , l 2 , l ∞ -norm, the dual gauge γ d (·) will be l ∞ , l 2 , l 1 -norm, respectively; and when γ(·) is a block norm, its dual is also a block norm.

3.2.
New projection-type methods. Among numerous effective numerical algorithms for solving VI, especially LVI, is the projection-type method. The attractive characteristics of the projection-type method, e.g. simpleness and effectiveness, have motivated further development on VI, see e.g. [12,13,29]. In this paper we aim to propose new projection-type methods for solving the LVI (18)- (19).
For a given vector v ∈ R n and a nonempty closed convex set Π, the projection of v onto Π is It is well known (see e.g. [8]) that for any β > 0, u * is a solution of the LVI (18)- (19) if and only if Then the residual function, which measures the magnitude of u away from the solution set, is defined as We propose a new projection-type method for solving the LVI (18)-(19) as follows.
Algorithm 1. The new projection-type method 1 for LVI (18)-(19) Given an initial iterate u 0 = x 0 z 0 . For k = 0, 1, 2 · · · , do: Step Step 2. Generate the new iterate u k+1 = x k+1 z k+1 : Step 3. If x k+1 − x k ≥ ε, k = k + 1 and go to Step 1; otherwise, stop. W is an invertible matrix, r, s and α are positive constants. They satisfy some mild conditions so as to guarantee the convergence, which will be discussed later. The ε is positive and used as the stopping criterion parameter.
Remark 4. Existing algorithms available for solving the LVI (18)- (19) include those presented in [12,13], etc. For example, applying the well-known projectioncontraction method in [13] for the LVI (18)- (19) results in the following procedure to generate new iterate: where α k = e(u k , 1) 2 / e(u k , 1) + M T e(u k , 1) 2 and γ ∈ (0, 2). Instead of applying directly existing algorithms, we present the simple method (22)- (24). At (23), the latest available informationz k is used to compute x k+1 ; while in the computation of x k+1 in (26), z k+1 is not used. This idea inheriting from the Gauss-Seidel method for solving system of nonlinear equations contributes greatly to the improvement and effectiveness of Algorithm 1, which will be verified by the numerical results later.
Next we will analyze the convergence of Algorithm 1. Assume that the solution set of the LVI (18)-(19) is nonempty and u * = x * z * is a solution of the LVI (18)- (19).
Step 1 of Algorithm 1 as Equivalently, Let the matrix Q and the mapping F (u) defined as Theorem 3.1. Let W be an invertible matrix and r, s and α be positive constants such that H := QW −1 0 and G := Q T + Q − αW T HW 0, then the sequences {u k } and {ũ k } generated by Algorithm 1 satisfies Proof. Setting u = u * in (29), then we have It follows that Since Combining (32) with (33), we get Setting u =ũ k in (18), then it follows from (34) that and therefore the following inequality is true By a simple manipulation, we have where the inequality is according to (36). Thus, this theorem follows from (37) directly and the proof is complete.
In fact, the condition of Theorem 3.1 is easy to be satisfied. The following proposition provides one possibility for the choice of matrixes and parameters.
Multiplying this matrix with an invertible matrix and its transpose on both sides, we obtain sI .
sI share the same property of semidefinite.
Remark 5. For a special case of α = 1, the condition of rs > A T A (2−α) 2 becomes rs > A T A to satisfy Theorem 3.1. 3) The sequence {u k } converges to a solution of the LVI (18)- (19).
Proof. The first assertion follows from (30) directly. From the inequality (30) we get and thus we have lim Consequently, lim and the second assertion is thus proved.
According to the first and second assertion we can conclude that the sequence {ũ k } is also bounded. Letũ ∞ be a cluster point of {ũ k } and the subsequence {ũ kj } converges toũ ∞ . Then it follows from (41) that and consequently, This means thatũ ∞ is a solution of the LVI (18)- (19).
Let u * =ũ ∞ in (30), then we can get which means that u k −ũ ∞ is monotone decreasing. Combining it with (44), the sequence {u k } converges to the solution pointũ ∞ . In Step 1 of Algorithm 1, we firstly computez k and then updatex k using the information of obtainedz k . Alternatively, we can computex k firstly and then updatẽ z k . Then another projection-type method Algorithm 2 is proposed as follows.
W is an invertible matrix, r, s and α are positive constants. The ε is positive and used as the stopping criterion parameter.
Noted thatũ k = x k z k in Step 1 of Algorithm 2 is equivalent to Let Q and F (u) defined as which is the same as (29) except for different Q. Similar to Theorem 3.1 and Theorem 3.2, we can prove the convergence of Algorithm 2 under similar mild conditions. Specifically, when W = I − A T r 0 I , 0 < α < 2 and rs > A T A (2−α) 2 (for a special case, α = 1 and rs > A T A ), the condition is satisfied and Algorithm 2 is convergent to a solution of LVI (18)- (19). 4. Numerical results. This section reports some numerical results to verify the computational efficiency of the variational inequality approach. We will implement the proposed Algorithm 1 and 2 to solve a number of CMFWPs including unconstrained problems and constrained problems and also including problems with symmetric distances and asymmetric distances. To demonstrate the efficiency of Algorithm 1 and 2, in the first subsection we compare Algorithm 1 and 2 with the well-known methods in facility location area for MFWP and in the second subsection we compare them with the well-known methods in variational inequalities. For simplicity, we choose α = 1, r = 2 t A T A and s = 1 2 t A T A with t = 1.01 in the following implementations of Algorithm 1 and 2. All the programming codes were written by MATLAB 7.8 and were run on a Dell laptop (Ceron processor 1.86GHz). 4.1. Algorithm 1 and 2 for MFWP. Miehle algorithm (MA) and hyperboloid approximation procedure (HAP) are perhaps the most standard methods for solving MFWP. This paper provides a variational inequality framework to solve CMFWP and under the new framework Algorithm 1 and 2 are then proposed. Note that MFWP, as a classical model in facility location, is a special case of CMFWP. Hence, the variational inequality framework and Algorithm 1 and 2 are applicable to solving MFWP. In the following, we compare Algorithm 1 and 2 with MA and HAP for solving MFWP and show that even for such a special problem the new framework and the proposed Algorithm 1 and 2 are efficient.
First we compare Algorithm 1 and 2 with MA. The following two examples were constructed in [26]. x 0 1 = (1, 0), x 0 2 = (2, 0). It has been shown in [26] 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 : E1 with initial iterate x 0 1 and x 0 2 randomly chosen as different customers. E2 : E2 with initial iterate x 0 1 and x 0 2 randomly chosen as different 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 Algorithm 1 and 2 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 Algorithm 1 and 2, we construct the following randomly generated problems. E3 : E1 (E2 ) but all customers a j and the initial iterate x 0 1 , x 0 2 are randomly generated in [−250, 250] 2 , w ij and v il are randomly generated in [1,100]. Note that MA is a Jacobi-like method because the new locations of different facilities change at the same time at each iteration. That is to say, x k+1 i is a function of x k 1 , · · · , x k i−1 , x k i , · · · , x k m . An alternative, however, may potentially accelerate the convergence of MA by using the new information as it becomes available, i.e., x k+1 i is obtained from x k+1 1 , · · · , x k+1 i−1 , x k i , · · · , x k m . This idea inherits from the Gauss-Seidel method and thus such alternative is called as Gauss-Seidel-like Miehle algorithm (GMA).
We solve the five examples by MA, GMA, Algorithm 1 and Algorithm 2 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 four methods with the same initial iterate x 0 1 and x 0 2 . The initial z 0 ij and z 0 il in Algorithm 1 and 2 are randomly chosen in {ξ ∈ R 2 | ξ ≤ w ij } and {ξ ∈ R 2 | ξ ≤ v il }, respectively. For E1, E2, E1 and E2 , each example is solved by Algorithm 1 and 2 one hundred times with randomly generated z 0 ij and z 0 il , and the average numerical results are reported. For E3, we test it for one hundred scenes with each scene solved by four 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 Algorithm 1 and 2 need more iterations for E1, they can get their optimal solution and provide the optimal objective functional value. For E2, MA needs the most number of iterations among four methods. With the adoption of the Gauss-Seidel idea, GMA needs the least number of iterations. When we compare the points obtained by four methods, Algorithm 1 and 2 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. Algorithm 1 and 2, however, are able to deal with the singularity and in such case Algorithm 1 and 2 can still get the optimal solution, which is consistent with the analytic results.
For the randomly generated E3 , it is found that the average computing time and number of iterations of GMA are the least. Nevertheless, the solution obtained by GMA is the worst among four methods. MA needs much more computing time and number of iterations than the other methods, but the point obtained by it is still not the optimal one of E3 . Compared to MA and GMA, Algorithm 1 and 2 can get the optimal solution of E3 within reasonable computing time and number of iterations.
Next we will compare Algorithm 1 and 2 with HAP. An alternative of HAP using the new available information was proposed in [27]. This alternative is called Gauss-Seidel-like HAP (GHAP). The parameter ε is chosen as 50, 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 50 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 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, Algorithm 1 and Algorithm 2 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 Algorithm 1 and other methods. Actually the value of "Dobj." is the functional value of other methods minus that of Algorithm 1. Therefore, if this value is positive, the point obtained by other method is worse than that by Algorithm 1; otherwise, it is better. It is clear that the larger the value of "Dobj.", the worse the point obtained by the method. As well known for HAP, the larger the value of ε, the faster the convergence of HAP but the worse the accuracy of approximation to MFWP, which can be also illustrated in Column 3-6 of Table 2. The same conclusion can be easily drawn for GHAP. Table 2 also reveals that GHAP can accelerate the convergence of HAP with the adoption of Gauss-Seidel idea. The difference between the values of "Dobj." of HAP and GHAP is quite small for the same ε, which indicates that the points obtained by two methods are almost the same. However, the workload of GHAP is nearly half of HAP for the same ε.
Note that the value of ε in SHAP is successively reduced from 50 to 0.05. It implies that the points obtained by SHAP are highly possible to be identical to those obtained by HAP under ε=0.05. This can be verified by the almost same values of "Dobj." in Column 6 and Column 11. However, we find that the workload of SHAP is much less than that of HAP. Thus, SHAP can improve the efficiency of HAP greatly. A similar conclusion is also true for GHAP and SGHAP, indicated by Column 10 and Column 12 of Table 2.
According to the last two columns in Table 2, it can be seen that the performances of Algorithm 1 and 2 are quite similar. This founding is consistent with our expectation, since the difference between Algorithm 1 and Algorithm 2 is just different computation order ofx k andz k in Step 1.
When we compare Algorithm 1 and 2 with four HAP-type methods, it is easy to observe that Algorithm 1 and 2 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 Algorithm 1 and 2 are the optimal solution of MFWP. However, the points obtained by HAP-type methods are not the optimal solution, which is indicated by the positive values of "Dobj.". Secondly, we compare six methods from the computing workload. It is clear that the workloads of Algorithm 1 and 2 are comparable with HAP-type methods especially for large-scale MFWP. When m is small, the workloads of Algorithm 1 and 2 are greater than HAP-type methods except HAP under ε=0.05. However, with the value of m increasing, the increment speed of workload of HAP-type methods is much faster than Algorithm 1 and 2. When m is equal to 10, the workloads of HAP-type methods have exceeded Algorithm 1 and 2 except HAP under ε=50 and GHAP under ε=50 (note that the approximation accuracy of HAP under ε=50 and GHAP under ε=50 is low). 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 Algorithm 1 and 2 with some well-known methods in facility location area, their efficiency is verified. 4.2. Algorithm 1 and 2 for CMFWP. In this subsection we will apply Algorithm 1 and 2 to solve CMFWP (7). Note that Algorithm 1 and 2 are proposed for CMFWP under the variational inequality framework by solving LVI (18)- (19). There exist other methods for LVI, and thus they are also applicable to solving CMFWP. To verify the efficiency of Algorithm 1 and 2, we compare Algorithm 1 and 2 with some well-known methods for LVI to solve CMFWP. The first method to be compared is the projection and contraction (PC) method presented in [13] and the second one is the customized proximal point algorithm (PPA) proposed in [14]. In the customized PPA, when the metric proximal parameter G is chosen as it results in the following iterative scheme (denoted as PPA1) . Alternatively, we can also choose another metric proximal parameter G It results in the following scheme (denoted as PPA2) The convergence of PPA1 and PPA2 under the condition of rs > A T A and γ ∈ (0, 2) can be found in [14].
The parameters of different methods are given as follows.
• PC method: γ = 1.9; • PPA1 and PPA2: r = p · t A T A , s = 2 p · t * A T A , t = 1.01, p = 2; γ = 1.2. We first test CMFWP under Euclidean distances with n = 50, 100, 200, 500, 1000 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 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 CMFWPs for 100 scenes with each scene solved by five methods using the same randomly generated initial iterates. Table  3 reports the average computing time in seconds (CPU) and number of iterations (Iter.) of five methods to solve the CMFWPs under Euclidean distances (symmetric distances). According to Table 3, five methods can solve the randomly generated CMFWPs under Euclidean distance and with locational constraints efficiently. It is shown that the proposed Algorithm 1 and 2 are preferable to the PPA1 and PPA2 for solving the CMFWPs of all scales. This can be illustrated roughly as follows. For PPA1 and PPA2, while the iterative scheme of their first variable (z k for PPA1 andx k for PPA2) is simple to implement, the iteration of their second variable is complex, which will result in the unbalance of the execution of two variables. However, during the implementation of Algorithm 1 and 2, the iterative scheme of both variables,x k andz k , of Algorithm 1 and 2 is not only quite simple but also balanced and thus its performance is improved. In addition, it is also found that Algorithm 1 and 2 are much preferable to PC method due to using the latest available informationz k to computex k , which powerfully verifies the conclusion in Remark 4 in computational sense.
Since Algorithm 1 and 2 can solve the CMFWP whose distances are measured by asymmetric gauge, we test some randomly generated CMFWPs under gauge with different scales. The asymmetric gauge γ(·) is generated with the unit ball set as To verify the efficiency of Algorithm 1 and 2 for CMFWP under gauge, we also compare Algorithm 1 and 2 with the other three methods simultaneously. The performance of five methods is reported in Table 4. It is easy to observe that the proposed Algorithm 1 and 2 are efficient for CMFWP under gauge and they outperform PPA1, PPA2 and PC method for the randomly generated asymmetric problem. In comparison of Table 3 and Table 4, however, we found that for five methods due to the asymmetric property of gauge the computing time and number of iterations for the CMFWP under asymmetric gauge are much greater than those for CMFWP under Euclidean distances.

5.
Conclusion. This paper investigates the constrained multifacility Weber problem. The locational constraints on facilities are considered and the gauge is employed to measure the distances in this problem. In this paper, a variational inequality approach is developed for solving CMFWP. With the consideration of special structures of reformulated linear variational inequality, we propose new projectiontype methods. The convergence of proposed methods is proved under mild assumption. Preliminary numerical results justify the preferable merits of the proposed algorithms.
How to solve other facility location problems, such as multi-source location problem, by the approach of variational inequality deserves further and more extensive investigation in our future research.