A Heuristic Algorithm for Constrained Multi-Source Location Problem with Closest Distance under Gauge : The Variational Inequality Approach

and Applied Analysis 3 When B is symmetric around the origin, according to the definition of (1), we have γ (−x) = γ (x) , ∀x ∈ R p . (4) Combining (3) and (4), it follows that D(x, y) = D (y, x) , (5) whichmeans that the distancemeasured by γ(⋅) is symmetric. On the contrary, when B is not symmetric around the origin, (4) does not hold any more, and thus, the distance measured by the gauge γ(⋅) is asymmetric. Thus, when different compact convex sets are used as unit balls, different gauges (symmetric or asymmetric) can be generated and employed to measure distances in location problems, which depends on the requirements of practical applications. LetΛ = {Aj ⊂ R p | j = 1, . . . , n}denote the set of regional customers, and let xi(i = 1, . . . , m) be the location of the ith facility. Each regional customer Aj is simply assumed to be closed and convex. We denote the closest point in Aj to the facility xi by aj (xi) := argmin {D (q, xi) | q ∈ Aj} = argmin {γ (xi − q) | q ∈ Aj} . (6) Then, the closest distance between the facility xi and the regional customer Aj, denoted by dj(xi), can be represented by dj (xi) := min q∈Aj D(q, xi) = γ (xi − aj (xi)) . (7) When the gauge is used tomeasure distances, we have the following proposition for aj(x) and dj(x) which is similar to that in [16]. Proposition 1. Let x be the location of a facility; then, the closest point aj(x) in (6) is well defined, and the closest distance dj(x) in (7) is a convex function of x. Proof. Since Aj is a convex set and γ(⋅) is a convex function, (6) is a convex problem, and thus, aj(x) is well defined. Now, we prove that dj(x) is a convex function of x as follows. Let x and y be two points in R and λ ∈ [0, 1]; then, due to aj(x) and aj(y) inAj and the convexity ofAj, it follows that λaj(x) + (1 − λ)aj(y) ∈ Aj, and thus, we have


Introduction
Due to the large number of practical applications in various fields such as marketing, urban planning, supply chain management, and transportation, the continuous facility location problem has aroused the attention of many researchers ever since the pioneering work [1,2].For a comprehensive review on this topic, see, for example, [3,4].More specifically, the continuous facility location problem in a space is to seek the optimal locations for facilities serving customers (also called demand points), with certain objectives such as minimizing the sum of distances between facilities and customers.
The vast majority of literature treats the locations of facilities and customers as points in a space.Hence, the distances between facilities and customers are just the point-to-point distance without any ambiguity.In practical applications, however, regional demand arises frequently in such scenarios as uncertain demand, mobile demand, or cumbersome discrete situation whose number of demand points is extremely large.For such scenarios, many authors (e.g., [5][6][7][8][9][10][11]) promote that the regional customer, that is, the locations of customers are geometrically connected regions rather than points, should be considered.Therefore, in this paper, we consider the case of regional customer.
The question to be emphasized here, however, is how to measure the proximity from a regional customer to a facility.In the literature, different kinds of distances have been used to measure the required proximity.For example, the average distance evaluated by the proximity between the facility and some mean point in the interior of a regional customer (e.g., [10][11][12][13]) and the farthest distance measured by the proximity between the facility and the farthest point on the boundary of a regional customer (e.g., [8,9,14,15]).Definitely, the regional customers are treated by the average fashion when the average distance is considered, while the farthest distance realizes the worst-off nature in the sense that the regional customers are represented by their respective farthest points.In some real-life applications, however, the best-off nature is The rest of this paper is organized as follows.Section 2 states the formulation of our problem, which is shown to be nonconvex and NP-hard.The spirit of the well-known location-allocation heuristic algorithm, which consists of a location phase and an allocation phase in each iteration, is also discussed in this section.In Section 3, the subproblems arising in location phase and allocation phase are solved.More specifically, for the subproblem arising in allocation phase, the regional customers are allocated to facilities according to the nearest reclassification (NCR) heuristic, whereas for the single-source subproblem arising in location phase, the relationship between the subproblem and a monotone linear variational inequality (LVI) is firstly built, and then, a projection-contraction (PC) method is adopted to find the optimal location of the facility.A new locationallocation heuristic algorithm is proposed in Section 4 for solving our targeted problem, and its convergence is proved in Section 5. Preliminary numerical results are reported in Section 6 to verify the efficiency of the proposed algorithm.Finally, some conclusions are drawn in Section 7.

Model Description
This paper focuses on finding locations in the space   for  facilities, with the objective to minimize the sum of weighted closest distances between these facilities and  regional customers.Note that the minisum single-source models with closest distance have been well justified in [6,16] where the distances are particularly measured by  1 -norm and  2 -norm, respectively.In our multi-source model, however, the gauge is used to measure the distances between facilities and regional customers, and thus, both the symmetric distance (including the  1 -norm used in [16] and the  2 -norm used in [6]) and the asymmetric distance are considered.
Let  be a compact convex set in   , and let the interior of  contain the origin; then, the gauge of  is defined by is called the unit ball of (⋅) due to This way to define a gauge from a compact convex set was first introduced by Minkowski [29].The gauge (⋅) satisfies the following properties: (1) () ≥ 0 and () = 0 ⇔  = 0; (2) () = () for any  ≥ 0; (3) ( + ) ≤ () + () for any  and  ∈   .
It follows from (2) and (3) that any gauge () is a convex function of .The distance measuring function can be derived from a gauge by When  is symmetric around the origin, according to the definition of (1), we have Combining ( 3) and ( 4), it follows that which means that the distance measured by (⋅) is symmetric.
On the contrary, when  is not symmetric around the origin, (4) does not hold any more, and thus, the distance measured by the gauge (⋅) is asymmetric.Thus, when different compact convex sets are used as unit balls, different gauges (symmetric or asymmetric) can be generated and employed to measure distances in location problems, which depends on the requirements of practical applications.Let Λ = {  ⊂   |  = 1, . . ., } denote the set of regional customers, and let   ( = 1, . . ., ) be the location of the th facility.Each regional customer   is simply assumed to be closed and convex.We denote the closest point in   to the facility   by Then, the closest distance between the facility   and the regional customer   , denoted by   (  ), can be represented by When the gauge is used to measure distances, we have the following proposition for   () and   () which is similar to that in [16].

Proposition 1.
Let  be the location of a facility; then, the closest point   () in ( 6) is well defined, and the closest distance Proof.Since   is a convex set and (⋅) is a convex function, (6) is a convex problem, and thus,   () is well defined.Now, we prove that   () is a convex function of  as follows.Let  and  be two points in   and  ∈ [0, 1]; then, due to   () and   () in   and the convexity of   , it follows that   () + (1 − )  () ∈   , and thus, we have Therefore,   () is convex with respect to , and the proof is complete.
Based on the notations introduced above, now the constrained multi-source location problem (abbreviated as CMLP) we consider in this paper takes the following formulation: where   ≥ 0 is the given demand required by the th customer,  = (  1 , . . .,    )  is the variable of the locations of facilities to be determined,  = (  ) × is the undetermined variable of   which denotes the unknown allocation from the th facility to the th customer, and Π  is the locational constraint for the th facility which is assumed to be a convex and closed set in   .
More explanations are required for our model (9).First, the locational constraint Π  in (9) can also be chosen as   , and if all Π  ( = 1, . . ., ) are   , the CMLP ( 9) is a unconstrained problem, and thus, both the constrained and unconstrained problems are considered in our model.Second, mark that the minisum models analyzed in [6,16] are two special cases of CMLP (9), where  = 1, Π  =   , and () is particularly  1 -norm in [16] and Euclidean norm in [6].
It is noted that, with the presupposition that each facility is capable of providing sufficient services for the targeted customers, each customer is ultimately served only by the nearest facility in order to minimize the total sum of weighted distances.Therefore, the mathematical model CMLP also has the following form: where Π is the cartesian product of locational constraints; that is, When all   are points not regions and all Π  are   , CMLP (9) reduces to the well-known multi-source Weber problem (MWP) which has wide applications in operations research, marketing, urban planning, and so forth; see, for example, [3,30,31].Recall the fact that the multi-source Weber problem is nonconvex [32] and NP-hard [33], and therefore, heuristics algorithms are extremely popular and highly appreciated for overcoming the difficulty caused by its nonconvexity and NP-hardness; see, for example, [3,30,[34][35][36][37].In particular, the classical location-allocation heuristic, also called the Cooper algorithm, has received much attention ever since it was presented originally by Cooper in [34] for MWP, whose attractive characteristic is that each iteration consists of a location phase and an allocation phase.Now, as a more general problem of MWP, CMLP is also nonconvex and NP-hard.Hence, in this paper, we are interested in applying the location-allocation heuristic to solve the CMLP (9).Accordingly, some location subproblems and allocation subproblems occur.To clarify it, let M = {1, 2, . . ., }, N = {1, 2, . . ., }, and then Λ = {  :  ∈ N}.At the ( − 1)th iteration, let After the location phase, the allocation phase then revises the current partition of Λ to generate a new disjoint partition of Λ = {Λ  1 , Λ  2 , . . ., Λ   } by the following nearest center reclassification (NCR) heuristic (see [38]): for some customer   ∈ Λ −1 ℎ ( ∈ {1, 2, . . ., } and ℎ ∈ {1, 2, . . ., }), if is the nearest point for   among all    computed by (11), then (11) is the nearest facility for each regional customer in Λ −1  for any  ∈ {1, 2, . . ., }, then    ( = 1, 2, . . ., ) are the desirable locations of facilities and stop.Otherwise, we set  =  + 1 and repeat the iterations.

The Subproblems in Location and Allocation Phases
In this section, we will discuss the subproblems arising in the location phase and allocation phase.The allocation phase will partition the customers to  clusters by the nearest center reclassification (NCR) heuristic, and the location phase will find the optimal location for each cluster by solving  CSLPs (11).

Nearest Center Reclassification for Allocation of Customers.
The implementation of NCR heuristic to allocate regional customers can be executed by the following framework; see [38] for more details about this heuristic.

The Variational Inequality Approach for CSLP
According to the spirit of location-allocation heuristic algorithm, our central task for the CMLP is to solve CSLP (11) in location phase by an efficient means.Recall that the CSLP is a generalized problem of the minisum models discussed in [6,16], where Π  =   and the gauge (⋅) are the particular  1norm in [16] and  2 -norm in [6].For the model under  1 -norm in [16], by taking advantage of the piecewise linearity of the objective function, this model can be reduced to a standard minisum problem which can be easily solved by obtaining a median point for each coordinate separately.For the minisum model under  2 -norm in [6], an efficient Weiszfeld-type method is proposed, and the convergence of this method is analyzed.Similar to Weiszfeld procedure [2], one problem of the proposed method is that the singular case, that is, the current iterate happens to be within some location of regional customers, may occur during its implementation.Due to the use of the gradient of objective function in the iteration, this method will terminate unexpectedly once the singular case occurs.In order to tackle the undesirable singular case and make the Weiszfeld-type method computational effective, the authors suggest to ignore the gradient of ‖ −   ()‖ if  ∈   and then add an extra descent check and a boundary check to the iteration.As pointed out by Theorem 1 in [6], however, the sequence generated by the proposed Weiszfeld-type method is possible to be convergent to a nonoptimal point which is on the boundary of the regional customer.
In this section, a variational inequality approach is proposed to solve the general CSLP (11), where the locational constraints are imposed to the facility and the gauge is used as distance measuring function.Note that the study of variational inequality has received much attention due to its various applications arising in engineering, operations research, economics, transportation, and so forth; see, for example, [39][40][41][42][43][44][45].Specifically, the CSLP (11) considered in this paper is first reformulated into an equivalent linear variational inequality (LVI), and then, a projection-contraction (PC) method is adopted to solve the LVI.Consequently, a sequence will be generated by the variational inequality approach, which is shown to be convergent to the optimal location of the facility    of ( 11) even in the singular case.In addition, the closest points to the facility and the dual vectors with respect to the gauge can also be obtained from the generated sequence.
For convenience and succinctness, with the assumption that Λ −1  contains  customers, throughout this section, we ignore some superscripts and subscripts in (11) and consider the simplified model of ( 11) without confusion: According to Proposition 1, it follows that CSLP (11), or equivalently MCSLP (12), is convex problem of .

LVI Reformulation of MCSLP.
For the gauge (⋅) in ( 12) which is defined by ( 1), there exists a dual gauge,   (⋅), defined by Let   be the unit ball of the dual gauge   (⋅), which is also convex and compact and exactly the polar set of .The dual gauge of   (⋅) is again (⋅), that is, which can also be rewritten as For more details about gauge and dual gauge, as well as their unit balls, the readers can be referred to [27].
According to (15), MCSLP ( 12) is equivalent to the following min-max problem: where each   is a vector in Denote and let ( * ,  * ,  * ) ∈ Π××B  be the solution of ( 17); then, it follows that ( * ,  * ,  * ) is the saddle point of the objective function Thus, ( * ,  * ,  * ) is the solution of the following linear variational inequality: A compact form of ( 20) is where Note that  in ( 22) is a skew-symmetric matrix, then it is positive semidefinite, and thus, the linear variational inequality ( 21)-( 22) is monotone.
Based on the deduction above, we know that if ( * ,  * ,  * ) is the solution of (17), that is,  * is the solution of the MCSLP (12), then ( * ,  * ,  * ) will be the solution of the LVI ( 21)- (22).Further, we can prove that the MCSLP (12) and the LVI ( 21)- (22) are equivalent in the following theorem.Theorem 3. The MCSLP (12) and the LVI (21)-( 22) are equivalent in the sense that they have the same solution of  ∈ Π.
Remark 4. It is worth pointing out that the equivalence between the MCSLP (12) and the linear variational inequality ( 21)-( 22) can also be obtained by the duality theory and the variable   and the set     ( = 1, . . ., ) in ( 16) are, respectively, the dual vector and dual ball in the space   which satisfy   ∈     .
The norms especially  1 ,  2 , and  ∞ are frequently used to measure distances in the literature; see, for example, [18,19].It should be noted that the gauge used in this paper is an extension of norms which include  1 ,  2 , and  ∞ .When the gauge (⋅) is chosen as the  1 ,  2 , and  ∞ -norm, the dual gauge   (⋅) will be the  ∞ ,  2 , and  1 -norm, respectively.Let where ‖ ⋅ ‖ 2 is the Euclidean norm; then, as a particular case of our single-source location problem (11), the problem under  2 -norm analyzed in [6] can be reformulated into the LVI ( 21)- (22) in which B  is equal to B 2 and Π =   .Further let where ‖ ⋅ ‖ 1 and ‖ ⋅ ‖ ∞ are the  1 ,  ∞ -norm, respectively, and Then, the CSLP (11) under  1 -norm (the minisum model discussed in [16]) and the CSLP under  ∞ -norm are equivalent to the LVI ( 21)- (22), where B  is, respectively, equal to B ∞ and B 1 and the locational constraints Π are both   .

A Projection-Contraction Method for LVI (21)-(22).
Among numerous effective numerical algorithms for solving VI, especially LVI, one famous one is the projectioncontraction (PC) method which was originally proposed by Uzawa [46].The attractive characteristics of the PC method, for example, simpleness and effectiveness, have motivated further development on VI especially in computational aspects; see, for example, [39,[47][48][49].In this section, we will summarize some concepts and results about linear variational inequalities and then adopt the projection-contraction method in [48] for solving LVI ( 21)- (22).More details about the proposed PC method can be referred to [48].
Let  be a nonempty closed convex set of  Q .For a given V ∈  Q , the projection of V onto  denoted by   (V) is the unique solution of the following problem: A basic proposition of the projection mapping on a closed convex set is It is well known (see, e.g., [50]) that for any  > 0,  * is the solution of LVI(Ω, , ) if and only if In the literature of variational inequalities, (, ) is usually called the error bound of LVI, and it quantitatively measures how much  fails to be the solution of LVI(Ω, , ).Therefore, (, ) can serve as the stopping criterion for solving LVI(Ω, , ) iteratively.
Step 2. Calculate (  ) and set (  ) as Step 3. Calculate  +1 as Step 4. Adjust  as follows and set Let  =  + 1 and go to Step 1.
Remark 7. In Step 4 of Algorithm 6, the parameter  is selfadaptive during the iterations according to the value of (  ).

A Location-Allocation Heuristic
for CMLP (9) Recall the fact that for the well-known multi-source Weber problem (MWP), heuristics algorithms are extremely popular and frequently used for overcoming its nonconvexity and NP-hardness.In particular, the location-allocation heuristic algorithm has drawn much attention ever since its presentation by Cooper [34].Note that the targeted CMLP ( 9) is an extension of the MWP and it is harder than MWP, and thus, in this paper, we also focus on applying the locationallocation heuristic algorithm for solving the CMLP in the spirit of Cooper's work.
Our previous analysis indicates that each iteration of the location-allocation heuristic algorithm to be presented consists of an allocation phase and a location phase.The allocation task generates a new disjoint partition of all the regional customers according to the principle of NCR as in the Cooper algorithm, and the location phase identifies the optimal locations for the current partition of customers via implementing the variational inequality approach for solving  CSLPs.
Mark that the CMLP (9) differs from MWP mainly in that the customers are represented by regions rather than points.Consequently, the CSLPs involved in the location phase are constrained location problems with regional demand and closest distances under gauge.No doubt that the numerical implementation of the heuristic algorithm to be presented is expected to be more complicated than the location-allocation algorithms for MWP.Therefore, how to accelerate the convergence of the proposed heuristic deserves further consideration.To achieve this objective, we here consider a particular strategy for the initial partition of regional customers or the initial locations of facilities.In practical implementation, we suggest to choose the solution of the following constrained multi-source Weber problem (CMWP) as the initial locations of facilities for CMLP: CMWP: min where   's are geometric centers of the regional customers.Then, we apply the NCR to determine an initial partition of regional customers according to the solution of (42).For solving the constrained multi-source Weber problem (42), we employ the location-allocation heuristic algorithm in [35].As we will show by numerical experiments, this initialization strategy can accelerate the convergence of the proposed algorithm greatly.
In the spirit of Cooper's work, the new heuristic algorithm is ready to be presented for solving the targeted CMLP (9), and its iterative framework can be elaborated as follows.
Step 2 (allocation phase).Update the partition of regional customers Λ from Λ −1  to Λ   based on the spirit of NCR heuristic.
Step 3 If ‖  −  −1 ‖ < , the current locations and partition are heuristic locations of facilities and heuristic partition of customers.Otherwise, set  =  + 1 and go to Step 1.
Remark 10.At the ( + 1)th iteration, it is recommended to use    (and the corresponding  and ) in Step 1 as the initial iterate in the variational inequality approach for solving  +1  ( = 1, . . ., ), considering the fact that Λ   usually differs from Λ −1  slightly in practical implementation.
Remark 11.Compared to the main body of the proposed location-allocation heuristic algorithm, the workload of the initialization is relatively less.However, this initialization strategy can reduce its number of iterations and computing time, which will be verified by the numerical experiments to be reported in Section 6.2.Hence, the convergence of the proposed algorithm is accelerated greatly by this initialization strategy.

Convergence of the Proposed Heuristic Algorithm
In this section, we analyze the convergence of the proposed location-allocation heuristic (Algorithm 9).For simplification of our discussion, some notations are introduced as follows.Let  =  1 × ⋅ ⋅ ⋅ ×   , and recall Π = Π 1 × ⋅ ⋅ ⋅ × Π  .
During the implementation of Algorithm 9, we denote the map L : Π ×  → Π ×  as the location operation in Step 1 and the map A : Π× → Π ×  as the allocation operation in the Step 2. It follows that where   and  +1 are, respectively, the locations of facilities in the th and ( + 1)th iteration,   is the variable of the closest points to   in the partition of Λ  , C is the closest points to the new  +1 in the partition of Λ  , and  +1 is the closest points to the new  +1 in the new partition of Λ +1 .Then, the iterate scheme of the location-allocation heuristic is Let ( 0 ,  0 ) denote the iterative sequence generated by the location-allocation heuristic for CMLP with the initial iterate ( 0 ,  0 ).During the implementation of the proposed heuristic algorithm, we choose the initial iterate in location phase for solving CSLP as Remark 10 indicates.We first give the following proposition which reveals the monotonicity of the generated sequence ( 0 ,  0 ).Proposition 12. ( 0 ,  0 ) is strictly monotone in the sense that ( +1 ,  +1 ) < (  ,   ) if  +1 ̸ =   .
Proof.Since  +1 ̸ =   , there exists at least one  ∈ M such that  +1  ̸ =    .For such 's according to the following convex optimization problem we have Based on Remark 10, we know that if    is the solution of ( 46), then  +1  will be equal to    .Therefore,  +1  ̸ =    implies that    is not the solution of ( 46), and thus,  +1  ( +1 ) <    (  ).It follows that  ( +1 , C ) <  (  ,   ) . ( On the other hand, based on the principle of NCR in the allocation phase of Algorithm 9, we also have Theorem 16.Assume that  +1 ̸ =   for any  ∈ ; then, all the accumulation points of the sequence {  } belong to Ξ. Proof.Let  be an accumulation point of {  }.Due to   ∈ Π×  and the compactness of Π × , we know that (1)  ∈ Π ×  and (2) there exists a subsequence {  }  which is convergent to .According to Theorem 13, we know that lim According to Lemma 14, the map L is closed at  ∈ Π × ; then, it follows that  1 = L().Further, due to  ∉ Ξ, V 1 will be not equal to .Otherwise, we can choose  as the initial iterate of the location-allocation heuristic algorithm, then, the sequence generated by the algorithm will be constant.It follows from the first case (i.e.,  +1 =   ) that  will be a local solution, which contradicts with  ∉ Ξ.Therefore, we can obtain V 1 ̸ = .Together with the monotone proposition of L (48), we have the inequality  ( 1 ) <  () . (63) On the other hand, note that A( +1 , C ) = ( +1 ,  +1 ); then, by the monotonicity of A (49), it follows that and thus, by taking the limit for (64) and by the continuity of , we know ( 2 ) ≤ ( 1 ).Combining this with (63), we obtain However, note that  2 and  are two accumulation points of {  }, and according to the second assertion of Theorem 13, ( 2 ) = (), which will contradict with (65).Therefore, our assumption is wrong, and thus,  ∈ Ξ.
As a result, we have the following convergence theorem for the sequence generated by the proposed locationallocation algorithm.
Theorem 17.The sequence ( 0 ,  0 ) generated by the proposed location-allocation heuristic algorithm either converges to a point in Ξ or all accumulation points of ( 0 ,  0 ) belong to Ξ.

Numerical Results
This section reports some preliminary numerical results to verify the theoretical assertions proved in previous sections.Section 3.1, reports some numerical results of the proposed variational inequality approach for the CSLP (11) (or equivalently ( 12)) which includes (1) the results of the comparison between our approach and the Weiszfeld-type method by solving the example in [6] and some randomly generated unconstrained examples under Euclidean distances and (2) the results of our approach for solving some randomly generated constrained examples under a gauge.These numerical results demonstrate the efficiency of the proposed variational inequality approach for CSLP.In the second subsection, we apply the proposed location-allocation heuristic algorithm to solve some randomly generated examples of the CMLP (9).In particular, the effectiveness of the initialization strategy adopted in this heuristic for accelerating convergence will be justified.All the programming codes are written by Matlab 2012b and were run on an ASUS notebook (Intel Core2 Duo T6670 2.20 GHz).

Numerical Results of Variational Inequality Approach for CSLP.
When applying the variational inequality approach for solving CSLP and MCSLP (12), theoretically, the initial iteration  0 in Algorithm 6 can be chosen arbitrarily in Ω.In practical implementation, however, we choose  0 judiciously similar to the initialization strategy in the location-allocation heuristic: let   ( = 1, . . ., ) be the centers of regional customers, solve the following single-source Weber problem (SWP): by the projection-contraction method in [35], and then use its solution as the initial iterate for Algorithm 6.We call this the initialization strategy of variational inequality approach.In addition, throughout our experiments of VI approach, the  1 and  2 in Algorithm 6 are chosen as 1 and 0, respectively.We first solve the example given in [16] by the proposed variational inequality approach and the Weiszfeld-type method in [16].
In order to clarify the comparison of two methods, we choose the same stopping criterion as ‖ +1 −   ‖ ≤ 10 −4 (throughout this section, ‖ ⋅ ‖ is the  ∞ -norm).We test this example for 100 times with the same initial iterate for the two methods which is randomly generated in [0, 5] × [0, 3], and the numerical results including the location of new facility, the closest points to the facility, number of iterations, and computing time in units of second are reported in Table 1.According to Table 1, it follows that both methods can get the optimal location of the new facility.Though our variational inequality approach needs smaller iteration numbers and less computing time, both methods are efficient for solving Example 1 given in [6].
Recall that the sequence generated by the Weiszfeld-type method in [16] is possible to be convergent to a nonoptimal point on the boundary of the regional customer, as indicated in [16].The variational inequality approach, however, can obtain the optimal location of the new facility, which is guaranteed by the theoretical analysis in Section 3.2.To illustrate this, we compare two methods by solving the following particular example.
Example 19.Similar to Example 18, the number of customers  = 5, and all customers are unit squares whose edges are parallel to the axes.Let  1 be the distance between any two neighboring customers, and we set  1 equal to 1, 0.1, 0.01, and 0.001, respectively.The weights   ( = 1, . . ., 5) are randomly generated in the area of regional customers.
Note that in Example 19 the parameter  1 reflects how close the customers are away from its neighborhoods.We test this example for a large number of times with the stopping criterion ‖ +1 −  ‖ ≤ 10 −4 , and the average numerical results are reported in Table 2.In this table, each row reports the average results by testing Example 19 for one hundred times.The column of "No. of Iter.0 " gives the average iteration times of initialization in the variational inequality approach, and the column of "No. of Iter." reports the average iteration times of both methods.The two columns of "CPU" give the average computing time in units of second for variational inequality approach (including the computing time for initialization) and Weiszfeld-type method, respectively.The columns of "Obj." give the average objective functional value obtained by the two methods, and the column of "Impro.Percent" gives the improvement percentage in objective functional values of the VI approach to Weiszfeld-type method.Remark that the convergence of VI approach to the optimal location of new facility is guaranteed, and then the column of "Freq.Num." reports the frequency among one hundred times that the Weiszfeld-type method can get the same solution as VI approach; that is, it does not converge to the nonoptimal solution on the boundary of the customer.
According to Table 2, it follows that both the VI approach and the Weiszfeld-type method are efficient for solving this particular example, and both of them need a small number of iterations and little computing time.In comparison of two methods, we can find that VI approach needs more iteration times and computing time than Weiszfeld-type method.From the column of "Impro.Percent, " however, we can conclude that VI approach can obtain a better solution (in fact, the solution obtained by VI approach is the optimal location of new facility) than Weiszfeld-type method.In addition, according to the last column, we find that when  1 decreases, that is, the customers become closer and closer, the frequency that Weiszfeld-type method obtains the optimal solution gets smaller and smaller.When  1 is 0.001, which implies that the customers are quite close to the neighborhoods, this frequency is totally less than 20.In other words, for this particular example with  1 = 0.001, the sequence generated by Weiszfeld-type method has a great possibility (more than 80%) to be convergent to a nonoptimal solution which is on the boundary of the customer.On the contrary, when  1 is 1, which means that the customers are enough far away from one another, Weiszfeld-type method can obtain the optimal location of facility in most cases, and the frequency even exceeds 90.
Since our main effort in this paper is to solve the general location problem under gauge and locational constraint, it is necessary to apply the variational inequality approach to solve some CSLPs.In particular, we test a large number of randomly generated CSLPs with the number of customers  from 10 to 2000.In the experiments, all regional customers are assumed to be square units, and their edges are parallel to the coordinate axes.The geometric centers of all regional customers are randomly generated in [−100, 100] 2 ; the weights of the regional customer are all randomly chosen in (1,5); the locational constraint is ‖ − ‖ ≤ , where the center  is randomly generated in [−100, 100] 2 , and the radius  is randomly generated in (1,5); the stopping criterion of VI approach is chosen as and the initial iterate is randomly generated in [−100, 100] 2 ; the gauge (⋅) is generated with the unit ball set as For each , we test one hundred randomly generated CSLPs, and the average numerical results are reported in Table 3.
To illustrate the effect of the initialization strategy of VI approach, we also report the results of VI approach without the initialization strategy.The columns of "VI approach with Initial." and "VI approach without Initial., " respectively, report the average number of iterations and average computing time of variational inequality approach with and without the initialization strategy.According to Table 3, it is easy to conclude that the variational inequality approach is effective for solving CSLP under gauge considering the difficulty of this problem.In addition, the number of iterations of "VI approach with Initial." is less than that of "VI approach without Initial., " which shows that the initialization strategy can accelerate the convergence of the variational inequality approach.This strategy, however, does not necessarily reduce the computing time of VI approach, especially when  is small, for example,  = 10, 20, 50, and 100, which can be explained as follows.When the number of customers  is small, the variational inequality problem is small scale, and thus, it can be solved in a short time.In this case, the computational workload of initialization plays an important role in the total workload, and therefore, the computing time of "VI approach with Initial." is greater than that of "VI approach without Initial." due to the computational iterations for initialization.With  increasing, the scale of VI problem as well as the number of iterations becomes larger.Then, in the comparison of the workload of VI approach, the workload of initialization can almost be ignored, and therefore, the computing time of "VI approach with Initial." will be smaller than that of "VI approach without Initial." As a matter of fact, Table 3 reveals the computational necessity of the initialization strategy of variational inequality approach for large-scale CSLP; for example, the iteration number and computing time are reduced about 1/5 by the initialization strategy when  = 2000.

Numerical Results of Heuristic Algorithm for CMLP.
This subsection applies the proposed location-allocation heuristic algorithm (Algorithm 9) to solve a large number of CMLPs  (9) and also verifies the necessity of the initialization strategy in Algorithm 9.In the experiments, we again generate a large number of CMLPs with unit square customers and assume that the edges of these regions are parallel to the coordinate axes.The geometric centers of all the demand regions are randomly generated in [−100 100] 2 , and all the demands,   ( = 1, 2, . . ., ), are randomly generated in [ 1 10 ].The gauge (⋅) is also defined with the unit ball set as (68).We test the scenario with  = 100, 200, 500, and 1000 and  = 2, 4, 6, 8, and 10; the locational constraints are ‖ −   ‖ ≤   ( = 1, . . ., ), where the radius   is randomly generated in [1,10], and the center   is given in advance as follows: To show the significance of the initialization strategy, we compare the numerical performance of the locationallocation heuristic algorithm with initialization strategy (denoted by "Algorithm 9 with Initial.") and without this initialization strategy (denoted by "Algorithm 9 without Initial.").In the initialization step of Algorithm 9, the location-allocation algorithm in [35] is adopted to solve the corresponding CMWP (42), where a PC method is proposed to solve the subproblems in location phase, and the numbers of iterations of the algorithm in [35] (denoted by "Iter 0 ") and the average iteration numbers of PC method in one iteration of the algorithm (denoted by "Iter 0 -PC") are reported.The columns of "Iter." and "CPU, " respectively, report the number of iterations and computing time of Algorithm 9 with and without initialization strategy.Since the efficiency of Algorithm 9 is mainly determined by the number of iterations of the variational inequality approach, we also report the average number of iterations of Algorithm 6 in one iteration of Algorithm 9 (denoted by "Iter.-PC").For each given pair (, ), we test the CMLP for 100 times, and the computational performance is reported in Table 4.
It follows from Table 4 that the proposed locationallocation heuristic, with or without the initialization strategy, is capable of tackling the CMLP (9) efficiently, even for large-scale cases.Also, the necessity of the initialization strategy is evident.In fact, this strategy reduces both the number of iterations and the computing time by about 50%.
Another interesting fact obtained from Table 4 deserves further illustration.Recall that with the number of new facilities () increasing, the number of subproblems (CSLPs) in location phase increases too.According to Table 4, however, we find that for fixed number of customers (), with  increasing, the number of iterations for solving  CSLPs in one iteration of Algorithm 9 does not increase but almost decreases with , especially for large-scale CMLP.This can be illustrated roughly as follows.For fixed , when  increases, the average number of customers in each Λ   becomes smaller, which implies that the scale of the involved CSLP (11) in location phase is smaller.According to Table 3, it follows that we need smaller number of iterations for small-scale CSLP, and thus, the total number of iterations for solving  CSLPs decreases.Similarly, due to the same reason, the computing time for solving CMLP also decreases with  increasing, as reported in the column of "CPU" in Table 4.

Conclusion
In this paper, we are interested in the locations of multiple facilities in the space   with regional demands, where the closest distance is used to measure the proximities between facilities and customers.With locational constraints introduced for the locations of new facilities and with the gauge used as the distance measuring function, the problem considered in this paper has much more applications in practice.Due to its nonconvexity and NP-hardness, a new location-allocation heuristic algorithm is proposed to solve this problem, and its convergence is proved under mild assumptions.Some preliminary numerical experiments are reported to verify the computational efficiency of the proposed algorithm.

Table 2 :
Numerical results of VI approach and Weiszfeld-type method for Example 19.

Table 4 :
Numerical results of location-allocation heuristic for CMLP.