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

Authors: Jian-Lin Jiang; Saeed Assani; Kun Cheng; Xiao-Xing Zhu
Journal: Abstract and Applied Analysis
 (2013)

---

## Abstract

This paper considers the locations of multiple facilities in the spaceRp, with the aim of minimizing the sum of weighted distances between facilities and regional customers, where the proximity between a facility and a regional customer is evaluated by the closest distance. Due to the fact that facilities are usually allowed to be sited in certain restricted areas, some locational constraints are imposed to the facilities of our problem. In addition, since the symmetry of distances is sometimes violated in practical situations, the gauge is employed in this paper instead of the frequently used norms for measuring both the symmetric and asymmetric distances. In the spirit of the Cooper algorithm (Cooper, 1964), a new location-allocation heuristic algorithm is proposed to solve this problem. In the location phase, the single-source subproblem with regional demands is reformulated into an equivalent linear variational inequality (LVI), and then, a projection-contraction (PC) method is adopted to find the optimal locations of facilities, whereas in the allocation phase, the regional customers are allocated to facilities according to the nearest center reclassification (NCR). The convergence of the proposed algorithm is proved under mild assumptions. Some preliminary numerical results are reported to show the effectiveness of the new algorithm.


---

## Body

## 1. 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–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–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 of importance in facility location; see, for example, [6, 16, 17]. To realize precisely the desired best-off nature, we need to consider the closest distance; that is, the proximity from a regional customer to a facility is evaluated by the closest distance to this facility. Thus, in this paper, we consider the closest distance to measure such a proximity. Note that the three kinds of distances have been well justified in [14, 15] and the reader can be referred to them for more details.Focusing on the real-life applications with vast eyes, the regional customers and the closest distances are highly essential to be considered. An example is about the communication network where several central servers are required to be sited. The demand points on the network are partially connected forming groups, each containing a large number of demand points. The points in the same group are connected to one another, and thus, each group becomes a regional customer in the plane. The servers need to be connected to the closest points in each regional customer, and the rest of the regional customer will be connected to the servers through that connection. Hence, the regional customers and the closest distances need to be used in the locations of these central servers. For more applications, we refer to, for example, [6, 15, 16].As a matter of fact, in the literature of facility location, the distance is usually measured by norms such aslp-norms and block norms. References [3, 18] discuss the approximations about the weighted lp-norms based on statistical evidence. References [19, 20] investigate the use of block norms to obtain good approximation of actual distances. As it is known, the symmetry property of the norm assures that the distance from one point to another is always equal to the distance back. Nevertheless, as one of the first in-depth studies of mathematical location problems, [21] highlights the fact that in numerous real situations the symmetry of the distance is violated, for example, transportation in rush-hour traffic, flight in the presence of wind, navigation in the presence of currents, transportation on an inclined terrain, and so forth. For about two decades after the work of [21], however, no further research on this topic seems to be published, and only in the recent twenty years have the asymmetric distance problems started to attract the researchers’ interest, for example, [22–28].In this paper, we are interested in the locations of multiple facilities in the spaceRp with the aim of minimizing the sum of weighted distance between these facilities and regional customers, where the distance between a facility and a regional customer is evaluated by the closest distance. In addition, we formulate this problem in a quite general way with the aim of enhancing its practical applicability. First, we recognize that, usually, not all the points in the space Rp are possible locations; that is, new facilities are often allowed to be sited within the confines of the restricted areas. Therefore, we introduce locational constraints so that both the unconstrained and constrained problems are taken into consideration in this paper. Second, since the distance in many practical situations is not necessarily symmetric, the gauge is used to measure the distance instead of the widely used norms. With the more general distance measuring function used, both the symmetric and asymmetric distances can be considered in our problem.The rest of this paper is organized as follows. Section2 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 location-allocation 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.


## 2. Model Description

This paper focuses on finding locations in the spaceRp for m facilities, with the objective to minimize the sum of weighted closest distances between these facilities and n 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 l1-norm and l2-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 l1-norm used in [16] and the l2-norm used in [6]) and the asymmetric distance are considered.LetB be a compact convex set in Rp, and let the interior of B contain the origin; then, the gauge of B is defined by
(1)γ(x)=inf{λ>0∣xλ∈B},∀x∈Rp.Bis called the unit ball of γ(·) due to
(2)B={x∈Rp∣γ(x)≤1}.
This way to define a gauge from a compact convex set was first introduced by Minkowski [29]. The gauge γ(·) satisfies the following properties: (1)
γ(x)≥0 and γ(x)=0⇔x=0;(2)
γ(tx)=tγ(x) for any t≥0;(3)
γ(x+y)≤γ(x)+γ(y) for any x and y∈Rp.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
(3)D(x,y)=γ(y-x).
When B is symmetric around the origin, according to the definition of (1), we have
(4)γ(-x)=γ(x),∀x∈Rp.
Combining (3) and (4), it follows that
(5)D(x,y)=D(y,x),
which means that the distance measured 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⊂Rp∣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
(6)aj(xi):=argmin{D(q,xi)∣q∈Aj}=argmin{γ(xi-q)∣q∈Aj}.
Then, the closest distance between the facility xi and the regional customer Aj, denoted by dj(xi), can be represented by
(7)dj(xi):=minq∈AjD(q,xi)=γ(xi-aj(xi)).When the gauge is used to measure distances, we have the following proposition foraj(x) and dj(x) which is similar to that in [16].Proposition 1.
Letx 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.
SinceAj is a convex set and γ(·) is a convex function, (6) is a convex problem, and thus, aj(x) is well defined.
Now, we prove thatdj(x) is a convex function of x as follows. Let x and y be two points in Rp and λ∈[0,1]; then, due to aj(x) and aj(y) in Aj and the convexity of Aj, it follows that λaj(x)+(1-λ)aj(y)∈Aj, and thus, we have
(8)dj(λx+(1-λ)y)=mina∈AjD(a,λx+(1-λ)y)≤D(λaj(x)+(1-λ)aj(y),λx+(1-λ)y)=γ(λx+(1-λ)y-(λaj(x)+(1-λ)aj(y)))=γ(λ(x-aj(x))+(1-λ)(y-aj(y)))≤γ(λ(x-aj(x)))+γ((1-λ)(y-aj(y)))=λγ(x-aj(x))+(1-λ)γ(y-aj(y))=λdj(x)+(1-λ)dj(y).
Therefore, dj(x) is convex with respect to x, 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:(9)CMLP:minX,W∑j=1n‍∑i=1m‍wijdj(xi)=∑j=1n‍∑i=1m‍wijγ(xi-aj(xi))s.t.∑i=1m‍wij=sj,j=1,2,…,n,xi∈Πi,i=1,2,…,m,
where sj≥0 is the given demand required by the jth customer, X=(x1T,…,xmT)T is the variable of the locations of facilities to be determined, W=(wij)m×n is the undetermined variable of wij which denotes the unknown allocation from the ith facility to the jth customer, and Πi is the locational constraint for the ith facility which is assumed to be a convex and closed set in Rp.More explanations are required for our model (9). First, the locational constraint Πi in (9) can also be chosen as Rp, and if all Πi(i=1,…,m) are Rp, 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 m=1, Πi=Rp, and γ(x) is particularly l1-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:(10)CMLP′:minX∈ΠC(X)=∑j=1n‍sjmin1≤i≤mdj(xi)=∑j=1n‍sjmin1≤i≤mγ(xi-aj(xi)),
where Π is the cartesian product of locational constraints; that is, Π=Π1×⋯×Πm.When allAj are points not regions and all Πi are Rp, 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–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 ℳ={1,2,…,m}, 𝒩={1,2,…,n}, and then Λ={Aj:j∈𝒩}. At the (k-1)th iteration, let {Λ1k-1,Λ2k-1,…,Λmk-1} be the disjoint partition of Λ in the sense that ∪i=1mΛik-1=Λ and Λik-1∩Λjk-1=∅ (for i≠j), and each Λik-1 (i=1,…,m) in the partition is called one cluster. Then, at the kth iteration, the location phase finds the candidates of locations of facilities (denoted by x1k,x2k,…,xmk) by solving the following m constrained single-source location problems (CSLP) with the closest distance under gauge for each cluster Λik-1,  i=1,…,m:
(11)CSLP:xik=argminx∈Πi{=∑{j∈N:Aj∈Λik-1}‍sjminqj∈Ajγ(x-qj)Cik(x)=∑{j∈N:Aj∈Λik-1}‍sjminqj∈Ajγ(x-qj):=∑{j∈𝒩:Aj∈Λik-1}‍sjγ(x-aj(x))=∑{j∈𝒩:Aj∈Λik-1}‍sjminqj∈Ajγ(x-qj)}.
After the location phase, the allocation phase then revises the current partition of Λ to generate a new disjoint partition of Λ={Λ1k,Λ2k,…,Λmk} by the following nearest center reclassification (NCR) heuristic (see [38]): for some customer Aj∈Λhk-1 (j∈{1,2,…,n} and h∈{1,2,…,m}), if xlk (l≠h) is the nearest point for Aj among all xik computed by (11), then Λhk=Λhk-1∖{Aj} and Λlk=Λlk-1∪{Aj}. If xik solved by (11) is the nearest facility for each regional customer in Λik-1 for any i∈{1,2,…,m}, then xik (i=1,2,…,m) are the desirable locations of facilities and stop. Otherwise, we set k=k+1 and repeat the iterations.


## 3. 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 tom clusters by the nearest center reclassification (NCR) heuristic, and the location phase will find the optimal location for each cluster by solving m CSLPs (11).


### 3.1. 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.Algorithm 2 (the implementation of NCR).
Given an initial partitionΛ0={Λ10,Λ20,…,Λm0}.
For  k=1,2,…, do.
Step
1. Set t=0 (t stores the number of reassignments);
Step
2. Compute the facility xik of Λik-1 by solving CSLP (11), for i=1,2,…,m;
Step
3. For j=1,2,…,n do:  dij:=γ(xik-aj(xik))
fori=1,…,m;ifAj∈Λhk-1 and dlj=mini=1,…,m;i≠h{dij}<dhj,thenΛhk=Λhk-1∖{Aj}, Λlk=Λlk-1∪{Aj};t=t+1.
Step
4. If t=0, then the iteration terminates with {x1k,…,xmk} being the desirable locations for facilities and the customers in Λik-1 being served by xik (i=1,…,m).


### 3.2. The Variational Inequality Approach for CSLP (11)

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 Πi=Rp and the gauge γ(·) are the particular l1-norm in [16] and l2-norm in [6]. For the model under l1-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 l2-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 ∥x-aj(x)∥ if x∈Aj 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–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 xik 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Λik-1 contains d customers, throughout this section, we ignore some superscripts and subscripts in (11) and consider the simplified model of (11) without confusion:
(12)MCSLP:x=argminx∈Π{C(x):=∑j=1d‍sjγ(x-aj(x))}.
According to Proposition 1, it follows that CSLP (11), or equivalently MCSLP (12), is convex problem of x.


#### 3.2.1. LVI Reformulation of MCSLP

For the gaugeγ(·) in (12) which is defined by (1), there exists a dual gauge, γd(·), defined by
(13)γd(z)=max{zTx∣γ(x)≤1}.
Let Bd be the unit ball of the dual gauge γd(·), which is also convex and compact and exactly the polar set of B. The dual gauge of γd(·) is again γ(·), that is,
(14)γ(x)=max{zTx∣γd(z)≤1},
which can also be rewritten as
(15)γ(x)=maxz∈BdzTx.
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:
(16)minx∈Πmaxzj∈Bsjd∑j=1d‍zjT(x-aj(x)),
where each zj is a vector in Bsjd={ξ∈Rp∣γd(ξ)≤sj}. Since aj(x) is the closest point to x in Aj, we can introduce yj∈Aj to replace aj(x). Hence, (16) is equivalent to
(17)minx∈Π,yj∈Ajmaxzj∈Bsjd∑j=1d‍zjT(x-yj).
Denote
(18)y=(y1T,…,ydT)T,z=(z1T,…,zdT)T,A=A1×⋯×Ad,ℬd=Bs1d×⋯×Bsdd,
and let (x*,y*,z*)∈Π×A×ℬd be the solution of (17); then, it follows that (x*,y*,z*) is the saddle point of the objective function ∑j=1d‍zjT(x-yj); that is,
(19)∑j=1d‍zjT(x*-yj*)≤∑j=1d‍zj*T(x*-yj*)≤∑j=1d‍zj*T(x-yj),∀(x,y,z)∈Π×A×ℬd.
Thus, (x*,y*,z*) is the solution of the following linear variational inequality:
(20)x*∈Π,y*∈A,z*∈ℬd,(x-x*)T(∑j=1d‍zj*)≥0,∀x∈Π,(yj-yj*)T(-zj*)≥0,∀yj∈Aj,(zj-zj*)T(-(x*-yj*))≥0,∀zj∈Bsjd.
A compact form of (20) is
(21)LVI(Ω,M,q):u*∈Ω,(u-u*)T(Mu*+q)≥0,∀u∈Ω,
where u=(xT,yT,zT)T, Ω=Π×A×ℬd,
(22)M=(0N-NT0),N=(I2⋯⋯I2-I20⋯00⋱⋱⋮⋮⋱⋱00⋯0-I2),q=0.
Note that M 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(x*,y*,z*) is the solution of (17), that is, x* is the solution of the MCSLP (12), then (x*,y*,z*) 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 x∈Π.Proof.
Since the MCSLP (12) is equivalent to (17), we need to prove that (17) and the LVI (21)-(22) are equivalent. In the following, we will prove that (x*,y*,z*) is a solution of (17) if and only if (x*,y*,z*) is a solution of LVI (21)-(22).
Let(x*,y*,z*) be the solution of (17); then, according to the deduction above, we know that (x*,y*,z*) is the solution of LVI (21)-(22).
On the other hand, let(x*,y*,z*) be the solution of LVI (21)-(22) and ϕ(x,y,z)=∑j=1d‍zjT(x-yj); then, the inequality (19) is true, which means that (x*,y*,z*) is the saddle point of ϕ(x,y,z).
Note that(x*,y*,z*) is the saddle point of ϕ(x,y,z) if and only if (x*,y*,z*)∈Ω and
(23)maxz∈ℬdϕ(x*,y*,z)=ϕ(x*,y*,z*)=minx∈Π,y∈Yϕ(x,y,z*),
which implies that
(24)minx∈Π,y∈Ymaxz∈ℬdϕ(x,y,z)≤maxz∈ℬdϕ(x*,y*,z)=ϕ(x*,y*,z*)=minx∈Π,y∈Yϕ(x,y,z*)≤maxz∈ℬdminx∈Π,y∈Yϕ(x,y,z).
On the other hand, let z′ be any vector in ℬd; then, we have
(25)minx∈Π,y∈Yϕ(x,y,z′)≤minx∈Π,y∈Ymaxz∈ℬdϕ(x,y,z).
We choose z′ in (25) as the maximum point of the left term over z′∈ℬd; then,
(26)maxz∈ℬdminx∈Π,y∈Yϕ(x,y,z)≤minx∈Π,y∈Ymaxz∈ℬdϕ(x,y,z).
Combining (24) and (26), it follows that all terms in (24) are equal, and therefore,
(27)ϕ(x*,y*,z*)=minx∈Π,y∈Ymaxz∈ℬdϕ(x,y,z),
which implies that (x*,y*,z*) is the solution of (17).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 zj and the set Bsjd(j=1,…,d) in (16) are, respectively, the dual vector and dual ball in the space Rp which satisfy zj∈Bsjd.The norms especiallyl1, l2, and l∞ 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 l1, l2, and l∞. When the gauge γ(·) is chosen as the l1, l2, and l∞-norm, the dual gauge γd(·) will be the l∞, l2, and l1-norm, respectively. Let
(28)Bsj,2={ξ∈Rp∣∥ξ∥2≤sj},ℬ2=Bs1,2×⋯×Bsd,2,
where ∥·∥2 is the Euclidean norm; then, as a particular case of our single-source location problem (11), the problem under l2-norm analyzed in [6] can be reformulated into the LVI (21)-(22) in which ℬd is equal to ℬ2 and Π=Rp.Further let(29)Bsj,1={ξ∈Rp∣∥ξ∥1≤sj},Bsj,∞={ξ∈Rp∣∥ξ∥∞≤sj},
where ∥·∥1 and ∥·∥∞ are the l1,l∞-norm, respectively, and
(30)ℬ1=Bs1,1×⋯×Bsd,1,ℬ∞=Bs1,∞×⋯×Bsd,∞.
Then, the CSLP (11) under l1-norm (the minisum model discussed in [16]) and the CSLP under l∞-norm are equivalent to the LVI (21)-(22), where ℬd is, respectively, equal to ℬ∞ and ℬ1 and the locational constraints Π are both Rp.


#### 3.2.2. A Projection-Contraction Method for LVI (21)-(22)

Among numerous effective numerical algorithms for solving VI, especially LVI, one famous one is the projection-contraction (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–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].LetW be a nonempty closed convex set of R𝒬. For a given v∈R𝒬, the projection of v onto W denoted by PW(v) is the unique solution of the following problem:
(31)PW(v)=argmin{∥u-v∥2∣u∈W}.
A basic proposition of the projection mapping on a closed convex set is
(32)(v-PW(v))T(u-PW(v))≤0,∀v∈R𝒬,∀u∈W.It is well known (see, e.g., [50]) that for any β>0, u* is the solution of LVI(Ω,M,q) if and only if
(33)e(u*,β):=u*-PΩ[u*-β(Mu*+q)]=0.
In the literature of variational inequalities, e(u,β) is usually called the error bound of LVI, and it quantitatively measures how much u fails to be the solution of LVI(Ω,M,q). Therefore, e(u,β) can serve as the stopping criterion for solving LVI(Ω,M,q) iteratively.Let(34)e(u)=e(u,1),g(u)=MTe(u)+(Mu+q),φ(u)=e(u)T(Mu+q),
and Ω* be the set of solutions of LVI(Ω,M,q); then, for the positive semidefinite (not necessarily symmetric) matrix M, the following theorem can be obtained.Theorem 5 (Lemma 1 and Theorem 2 in [48]).
Letu∈Ω, u*∈Ω*, g(u), and φ(u) be defined as (34). Then, it holds that
(35)(u-u*)Tg(u)≥φ(u)≥∥e(u)∥22.Foru∈Ω∖Ω*, it follows from Theorem 5 that -g(u) is a descent direction of the unknown function ∥u-u*∥22. We state the projection-contraction method in [48] as follows which is used to solve the LVI (21)-(22).Algorithm 6 (the projection-contraction method for LVI (21)-(22)).
Step
0. Let ε>0, β=1, α1, α2 (α1>α2) and u0∈Ω. Set k=0.
Step
1. Calculate e(uk). If ∥e(uk)∥<ε, stop.
Step
2. Calculate g(uk) and set α(uk) as
(36)α(uk)=∥e(uk)∥2∥e(uk)+MTe(uk)∥2.
Step
3. Calculate uk+1 as
(37)uk+1=PΩ[uk-ρα(uk)g(wk)],ρ∈(0,2).
Step
4. Adjust β as follows
(38)β={32βα(uk)≥α1,23βα(uk)≤α2,βotherwise,
and set
(39)M=βM,q=βq.
Let k=k+1 and go to Step 1.Remark 7.
In Step 4 of Algorithm6, the parameter β is self-adaptive during the iterations according to the value of α(uk). Note that M is skew-symmetric, and thus, α(uk) can also be rewritten as
(40)α(uk)=∥e(uk)∥2(∥e(uk)∥2+∥MTe(uk)∥2),
which is shown to be in [0,1]. It follows from (39)-(40) that the two terms ∥e(uk)∥2 and ∥MTe(uk)∥2 in the denominator of α(uk) are balanced by the self-adaptive parameter β.Theorem 8 (Theorem 3 in [48]).
Letu* be a solution of LVI (21)-(22); then, the sequence {uk} generated by Algorithm 6 satisfies
(41)∥uk+1-u*∥2≤∥uk-u*∥2-ρ(2-ρ)∥I+MT∥2∥e(uk)∥2.As a result,{e(uk)} converges to zero, and thus, all accumulation points of {uk} are the solutions of LVI (21)-(22). However, it follows from (41) that ∥uk+1-u*∥ ≤ ∥uk-u*∥, which implies that {uk} has only one accumulation point. Thus, the sequence {uk} generated by Algorithm 6 will converge to the optimal solution of LVI (21)-(22).


## 3.1. 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.Algorithm 2 (the implementation of NCR).
Given an initial partitionΛ0={Λ10,Λ20,…,Λm0}.
For  k=1,2,…, do.
Step
1. Set t=0 (t stores the number of reassignments);
Step
2. Compute the facility xik of Λik-1 by solving CSLP (11), for i=1,2,…,m;
Step
3. For j=1,2,…,n do:  dij:=γ(xik-aj(xik))
fori=1,…,m;ifAj∈Λhk-1 and dlj=mini=1,…,m;i≠h{dij}<dhj,thenΛhk=Λhk-1∖{Aj}, Λlk=Λlk-1∪{Aj};t=t+1.
Step
4. If t=0, then the iteration terminates with {x1k,…,xmk} being the desirable locations for facilities and the customers in Λik-1 being served by xik (i=1,…,m).


## 3.2. The Variational Inequality Approach for CSLP (11)

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 Πi=Rp and the gauge γ(·) are the particular l1-norm in [16] and l2-norm in [6]. For the model under l1-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 l2-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 ∥x-aj(x)∥ if x∈Aj 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–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 xik 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Λik-1 contains d customers, throughout this section, we ignore some superscripts and subscripts in (11) and consider the simplified model of (11) without confusion:
(12)MCSLP:x=argminx∈Π{C(x):=∑j=1d‍sjγ(x-aj(x))}.
According to Proposition 1, it follows that CSLP (11), or equivalently MCSLP (12), is convex problem of x.


### 3.2.1. LVI Reformulation of MCSLP

For the gaugeγ(·) in (12) which is defined by (1), there exists a dual gauge, γd(·), defined by
(13)γd(z)=max{zTx∣γ(x)≤1}.
Let Bd be the unit ball of the dual gauge γd(·), which is also convex and compact and exactly the polar set of B. The dual gauge of γd(·) is again γ(·), that is,
(14)γ(x)=max{zTx∣γd(z)≤1},
which can also be rewritten as
(15)γ(x)=maxz∈BdzTx.
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:
(16)minx∈Πmaxzj∈Bsjd∑j=1d‍zjT(x-aj(x)),
where each zj is a vector in Bsjd={ξ∈Rp∣γd(ξ)≤sj}. Since aj(x) is the closest point to x in Aj, we can introduce yj∈Aj to replace aj(x). Hence, (16) is equivalent to
(17)minx∈Π,yj∈Ajmaxzj∈Bsjd∑j=1d‍zjT(x-yj).
Denote
(18)y=(y1T,…,ydT)T,z=(z1T,…,zdT)T,A=A1×⋯×Ad,ℬd=Bs1d×⋯×Bsdd,
and let (x*,y*,z*)∈Π×A×ℬd be the solution of (17); then, it follows that (x*,y*,z*) is the saddle point of the objective function ∑j=1d‍zjT(x-yj); that is,
(19)∑j=1d‍zjT(x*-yj*)≤∑j=1d‍zj*T(x*-yj*)≤∑j=1d‍zj*T(x-yj),∀(x,y,z)∈Π×A×ℬd.
Thus, (x*,y*,z*) is the solution of the following linear variational inequality:
(20)x*∈Π,y*∈A,z*∈ℬd,(x-x*)T(∑j=1d‍zj*)≥0,∀x∈Π,(yj-yj*)T(-zj*)≥0,∀yj∈Aj,(zj-zj*)T(-(x*-yj*))≥0,∀zj∈Bsjd.
A compact form of (20) is
(21)LVI(Ω,M,q):u*∈Ω,(u-u*)T(Mu*+q)≥0,∀u∈Ω,
where u=(xT,yT,zT)T, Ω=Π×A×ℬd,
(22)M=(0N-NT0),N=(I2⋯⋯I2-I20⋯00⋱⋱⋮⋮⋱⋱00⋯0-I2),q=0.
Note that M 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(x*,y*,z*) is the solution of (17), that is, x* is the solution of the MCSLP (12), then (x*,y*,z*) 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 x∈Π.Proof.
Since the MCSLP (12) is equivalent to (17), we need to prove that (17) and the LVI (21)-(22) are equivalent. In the following, we will prove that (x*,y*,z*) is a solution of (17) if and only if (x*,y*,z*) is a solution of LVI (21)-(22).
Let(x*,y*,z*) be the solution of (17); then, according to the deduction above, we know that (x*,y*,z*) is the solution of LVI (21)-(22).
On the other hand, let(x*,y*,z*) be the solution of LVI (21)-(22) and ϕ(x,y,z)=∑j=1d‍zjT(x-yj); then, the inequality (19) is true, which means that (x*,y*,z*) is the saddle point of ϕ(x,y,z).
Note that(x*,y*,z*) is the saddle point of ϕ(x,y,z) if and only if (x*,y*,z*)∈Ω and
(23)maxz∈ℬdϕ(x*,y*,z)=ϕ(x*,y*,z*)=minx∈Π,y∈Yϕ(x,y,z*),
which implies that
(24)minx∈Π,y∈Ymaxz∈ℬdϕ(x,y,z)≤maxz∈ℬdϕ(x*,y*,z)=ϕ(x*,y*,z*)=minx∈Π,y∈Yϕ(x,y,z*)≤maxz∈ℬdminx∈Π,y∈Yϕ(x,y,z).
On the other hand, let z′ be any vector in ℬd; then, we have
(25)minx∈Π,y∈Yϕ(x,y,z′)≤minx∈Π,y∈Ymaxz∈ℬdϕ(x,y,z).
We choose z′ in (25) as the maximum point of the left term over z′∈ℬd; then,
(26)maxz∈ℬdminx∈Π,y∈Yϕ(x,y,z)≤minx∈Π,y∈Ymaxz∈ℬdϕ(x,y,z).
Combining (24) and (26), it follows that all terms in (24) are equal, and therefore,
(27)ϕ(x*,y*,z*)=minx∈Π,y∈Ymaxz∈ℬdϕ(x,y,z),
which implies that (x*,y*,z*) is the solution of (17).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 zj and the set Bsjd(j=1,…,d) in (16) are, respectively, the dual vector and dual ball in the space Rp which satisfy zj∈Bsjd.The norms especiallyl1, l2, and l∞ 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 l1, l2, and l∞. When the gauge γ(·) is chosen as the l1, l2, and l∞-norm, the dual gauge γd(·) will be the l∞, l2, and l1-norm, respectively. Let
(28)Bsj,2={ξ∈Rp∣∥ξ∥2≤sj},ℬ2=Bs1,2×⋯×Bsd,2,
where ∥·∥2 is the Euclidean norm; then, as a particular case of our single-source location problem (11), the problem under l2-norm analyzed in [6] can be reformulated into the LVI (21)-(22) in which ℬd is equal to ℬ2 and Π=Rp.Further let(29)Bsj,1={ξ∈Rp∣∥ξ∥1≤sj},Bsj,∞={ξ∈Rp∣∥ξ∥∞≤sj},
where ∥·∥1 and ∥·∥∞ are the l1,l∞-norm, respectively, and
(30)ℬ1=Bs1,1×⋯×Bsd,1,ℬ∞=Bs1,∞×⋯×Bsd,∞.
Then, the CSLP (11) under l1-norm (the minisum model discussed in [16]) and the CSLP under l∞-norm are equivalent to the LVI (21)-(22), where ℬd is, respectively, equal to ℬ∞ and ℬ1 and the locational constraints Π are both Rp.


### 3.2.2. A Projection-Contraction Method for LVI (21)-(22)

Among numerous effective numerical algorithms for solving VI, especially LVI, one famous one is the projection-contraction (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–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].LetW be a nonempty closed convex set of R𝒬. For a given v∈R𝒬, the projection of v onto W denoted by PW(v) is the unique solution of the following problem:
(31)PW(v)=argmin{∥u-v∥2∣u∈W}.
A basic proposition of the projection mapping on a closed convex set is
(32)(v-PW(v))T(u-PW(v))≤0,∀v∈R𝒬,∀u∈W.It is well known (see, e.g., [50]) that for any β>0, u* is the solution of LVI(Ω,M,q) if and only if
(33)e(u*,β):=u*-PΩ[u*-β(Mu*+q)]=0.
In the literature of variational inequalities, e(u,β) is usually called the error bound of LVI, and it quantitatively measures how much u fails to be the solution of LVI(Ω,M,q). Therefore, e(u,β) can serve as the stopping criterion for solving LVI(Ω,M,q) iteratively.Let(34)e(u)=e(u,1),g(u)=MTe(u)+(Mu+q),φ(u)=e(u)T(Mu+q),
and Ω* be the set of solutions of LVI(Ω,M,q); then, for the positive semidefinite (not necessarily symmetric) matrix M, the following theorem can be obtained.Theorem 5 (Lemma 1 and Theorem 2 in [48]).
Letu∈Ω, u*∈Ω*, g(u), and φ(u) be defined as (34). Then, it holds that
(35)(u-u*)Tg(u)≥φ(u)≥∥e(u)∥22.Foru∈Ω∖Ω*, it follows from Theorem 5 that -g(u) is a descent direction of the unknown function ∥u-u*∥22. We state the projection-contraction method in [48] as follows which is used to solve the LVI (21)-(22).Algorithm 6 (the projection-contraction method for LVI (21)-(22)).
Step
0. Let ε>0, β=1, α1, α2 (α1>α2) and u0∈Ω. Set k=0.
Step
1. Calculate e(uk). If ∥e(uk)∥<ε, stop.
Step
2. Calculate g(uk) and set α(uk) as
(36)α(uk)=∥e(uk)∥2∥e(uk)+MTe(uk)∥2.
Step
3. Calculate uk+1 as
(37)uk+1=PΩ[uk-ρα(uk)g(wk)],ρ∈(0,2).
Step
4. Adjust β as follows
(38)β={32βα(uk)≥α1,23βα(uk)≤α2,βotherwise,
and set
(39)M=βM,q=βq.
Let k=k+1 and go to Step 1.Remark 7.
In Step 4 of Algorithm6, the parameter β is self-adaptive during the iterations according to the value of α(uk). Note that M is skew-symmetric, and thus, α(uk) can also be rewritten as
(40)α(uk)=∥e(uk)∥2(∥e(uk)∥2+∥MTe(uk)∥2),
which is shown to be in [0,1]. It follows from (39)-(40) that the two terms ∥e(uk)∥2 and ∥MTe(uk)∥2 in the denominator of α(uk) are balanced by the self-adaptive parameter β.Theorem 8 (Theorem 3 in [48]).
Letu* be a solution of LVI (21)-(22); then, the sequence {uk} generated by Algorithm 6 satisfies
(41)∥uk+1-u*∥2≤∥uk-u*∥2-ρ(2-ρ)∥I+MT∥2∥e(uk)∥2.As a result,{e(uk)} converges to zero, and thus, all accumulation points of {uk} are the solutions of LVI (21)-(22). However, it follows from (41) that ∥uk+1-u*∥ ≤ ∥uk-u*∥, which implies that {uk} has only one accumulation point. Thus, the sequence {uk} generated by Algorithm 6 will converge to the optimal solution of LVI (21)-(22).


## 3.2.1. LVI Reformulation of MCSLP

For the gaugeγ(·) in (12) which is defined by (1), there exists a dual gauge, γd(·), defined by
(13)γd(z)=max{zTx∣γ(x)≤1}.
Let Bd be the unit ball of the dual gauge γd(·), which is also convex and compact and exactly the polar set of B. The dual gauge of γd(·) is again γ(·), that is,
(14)γ(x)=max{zTx∣γd(z)≤1},
which can also be rewritten as
(15)γ(x)=maxz∈BdzTx.
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:
(16)minx∈Πmaxzj∈Bsjd∑j=1d‍zjT(x-aj(x)),
where each zj is a vector in Bsjd={ξ∈Rp∣γd(ξ)≤sj}. Since aj(x) is the closest point to x in Aj, we can introduce yj∈Aj to replace aj(x). Hence, (16) is equivalent to
(17)minx∈Π,yj∈Ajmaxzj∈Bsjd∑j=1d‍zjT(x-yj).
Denote
(18)y=(y1T,…,ydT)T,z=(z1T,…,zdT)T,A=A1×⋯×Ad,ℬd=Bs1d×⋯×Bsdd,
and let (x*,y*,z*)∈Π×A×ℬd be the solution of (17); then, it follows that (x*,y*,z*) is the saddle point of the objective function ∑j=1d‍zjT(x-yj); that is,
(19)∑j=1d‍zjT(x*-yj*)≤∑j=1d‍zj*T(x*-yj*)≤∑j=1d‍zj*T(x-yj),∀(x,y,z)∈Π×A×ℬd.
Thus, (x*,y*,z*) is the solution of the following linear variational inequality:
(20)x*∈Π,y*∈A,z*∈ℬd,(x-x*)T(∑j=1d‍zj*)≥0,∀x∈Π,(yj-yj*)T(-zj*)≥0,∀yj∈Aj,(zj-zj*)T(-(x*-yj*))≥0,∀zj∈Bsjd.
A compact form of (20) is
(21)LVI(Ω,M,q):u*∈Ω,(u-u*)T(Mu*+q)≥0,∀u∈Ω,
where u=(xT,yT,zT)T, Ω=Π×A×ℬd,
(22)M=(0N-NT0),N=(I2⋯⋯I2-I20⋯00⋱⋱⋮⋮⋱⋱00⋯0-I2),q=0.
Note that M 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(x*,y*,z*) is the solution of (17), that is, x* is the solution of the MCSLP (12), then (x*,y*,z*) 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 x∈Π.Proof.
Since the MCSLP (12) is equivalent to (17), we need to prove that (17) and the LVI (21)-(22) are equivalent. In the following, we will prove that (x*,y*,z*) is a solution of (17) if and only if (x*,y*,z*) is a solution of LVI (21)-(22).
Let(x*,y*,z*) be the solution of (17); then, according to the deduction above, we know that (x*,y*,z*) is the solution of LVI (21)-(22).
On the other hand, let(x*,y*,z*) be the solution of LVI (21)-(22) and ϕ(x,y,z)=∑j=1d‍zjT(x-yj); then, the inequality (19) is true, which means that (x*,y*,z*) is the saddle point of ϕ(x,y,z).
Note that(x*,y*,z*) is the saddle point of ϕ(x,y,z) if and only if (x*,y*,z*)∈Ω and
(23)maxz∈ℬdϕ(x*,y*,z)=ϕ(x*,y*,z*)=minx∈Π,y∈Yϕ(x,y,z*),
which implies that
(24)minx∈Π,y∈Ymaxz∈ℬdϕ(x,y,z)≤maxz∈ℬdϕ(x*,y*,z)=ϕ(x*,y*,z*)=minx∈Π,y∈Yϕ(x,y,z*)≤maxz∈ℬdminx∈Π,y∈Yϕ(x,y,z).
On the other hand, let z′ be any vector in ℬd; then, we have
(25)minx∈Π,y∈Yϕ(x,y,z′)≤minx∈Π,y∈Ymaxz∈ℬdϕ(x,y,z).
We choose z′ in (25) as the maximum point of the left term over z′∈ℬd; then,
(26)maxz∈ℬdminx∈Π,y∈Yϕ(x,y,z)≤minx∈Π,y∈Ymaxz∈ℬdϕ(x,y,z).
Combining (24) and (26), it follows that all terms in (24) are equal, and therefore,
(27)ϕ(x*,y*,z*)=minx∈Π,y∈Ymaxz∈ℬdϕ(x,y,z),
which implies that (x*,y*,z*) is the solution of (17).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 zj and the set Bsjd(j=1,…,d) in (16) are, respectively, the dual vector and dual ball in the space Rp which satisfy zj∈Bsjd.The norms especiallyl1, l2, and l∞ 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 l1, l2, and l∞. When the gauge γ(·) is chosen as the l1, l2, and l∞-norm, the dual gauge γd(·) will be the l∞, l2, and l1-norm, respectively. Let
(28)Bsj,2={ξ∈Rp∣∥ξ∥2≤sj},ℬ2=Bs1,2×⋯×Bsd,2,
where ∥·∥2 is the Euclidean norm; then, as a particular case of our single-source location problem (11), the problem under l2-norm analyzed in [6] can be reformulated into the LVI (21)-(22) in which ℬd is equal to ℬ2 and Π=Rp.Further let(29)Bsj,1={ξ∈Rp∣∥ξ∥1≤sj},Bsj,∞={ξ∈Rp∣∥ξ∥∞≤sj},
where ∥·∥1 and ∥·∥∞ are the l1,l∞-norm, respectively, and
(30)ℬ1=Bs1,1×⋯×Bsd,1,ℬ∞=Bs1,∞×⋯×Bsd,∞.
Then, the CSLP (11) under l1-norm (the minisum model discussed in [16]) and the CSLP under l∞-norm are equivalent to the LVI (21)-(22), where ℬd is, respectively, equal to ℬ∞ and ℬ1 and the locational constraints Π are both Rp.


## 3.2.2. A Projection-Contraction Method for LVI (21)-(22)

Among numerous effective numerical algorithms for solving VI, especially LVI, one famous one is the projection-contraction (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–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].LetW be a nonempty closed convex set of R𝒬. For a given v∈R𝒬, the projection of v onto W denoted by PW(v) is the unique solution of the following problem:
(31)PW(v)=argmin{∥u-v∥2∣u∈W}.
A basic proposition of the projection mapping on a closed convex set is
(32)(v-PW(v))T(u-PW(v))≤0,∀v∈R𝒬,∀u∈W.It is well known (see, e.g., [50]) that for any β>0, u* is the solution of LVI(Ω,M,q) if and only if
(33)e(u*,β):=u*-PΩ[u*-β(Mu*+q)]=0.
In the literature of variational inequalities, e(u,β) is usually called the error bound of LVI, and it quantitatively measures how much u fails to be the solution of LVI(Ω,M,q). Therefore, e(u,β) can serve as the stopping criterion for solving LVI(Ω,M,q) iteratively.Let(34)e(u)=e(u,1),g(u)=MTe(u)+(Mu+q),φ(u)=e(u)T(Mu+q),
and Ω* be the set of solutions of LVI(Ω,M,q); then, for the positive semidefinite (not necessarily symmetric) matrix M, the following theorem can be obtained.Theorem 5 (Lemma 1 and Theorem 2 in [48]).
Letu∈Ω, u*∈Ω*, g(u), and φ(u) be defined as (34). Then, it holds that
(35)(u-u*)Tg(u)≥φ(u)≥∥e(u)∥22.Foru∈Ω∖Ω*, it follows from Theorem 5 that -g(u) is a descent direction of the unknown function ∥u-u*∥22. We state the projection-contraction method in [48] as follows which is used to solve the LVI (21)-(22).Algorithm 6 (the projection-contraction method for LVI (21)-(22)).
Step
0. Let ε>0, β=1, α1, α2 (α1>α2) and u0∈Ω. Set k=0.
Step
1. Calculate e(uk). If ∥e(uk)∥<ε, stop.
Step
2. Calculate g(uk) and set α(uk) as
(36)α(uk)=∥e(uk)∥2∥e(uk)+MTe(uk)∥2.
Step
3. Calculate uk+1 as
(37)uk+1=PΩ[uk-ρα(uk)g(wk)],ρ∈(0,2).
Step
4. Adjust β as follows
(38)β={32βα(uk)≥α1,23βα(uk)≤α2,βotherwise,
and set
(39)M=βM,q=βq.
Let k=k+1 and go to Step 1.Remark 7.
In Step 4 of Algorithm6, the parameter β is self-adaptive during the iterations according to the value of α(uk). Note that M is skew-symmetric, and thus, α(uk) can also be rewritten as
(40)α(uk)=∥e(uk)∥2(∥e(uk)∥2+∥MTe(uk)∥2),
which is shown to be in [0,1]. It follows from (39)-(40) that the two terms ∥e(uk)∥2 and ∥MTe(uk)∥2 in the denominator of α(uk) are balanced by the self-adaptive parameter β.Theorem 8 (Theorem 3 in [48]).
Letu* be a solution of LVI (21)-(22); then, the sequence {uk} generated by Algorithm 6 satisfies
(41)∥uk+1-u*∥2≤∥uk-u*∥2-ρ(2-ρ)∥I+MT∥2∥e(uk)∥2.As a result,{e(uk)} converges to zero, and thus, all accumulation points of {uk} are the solutions of LVI (21)-(22). However, it follows from (41) that ∥uk+1-u*∥ ≤ ∥uk-u*∥, which implies that {uk} has only one accumulation point. Thus, the sequence {uk} generated by Algorithm 6 will converge to the optimal solution of LVI (21)-(22).


## 4. 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 location-allocation 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 solvingm 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:
(42)CMWP:min(x1,…,xm)∈Π1×⋯×ΠmC′(X)=∑j=1n‍sjmin1≤i≤mγ(xi-gj),
where gj’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.Algorithm 9 (a location-allocation heuristic algorithm for CMLP).
Initialization: Solve (42) by the location-allocation heuristic in [35] and use its heuristic solutions as the initial locations of facilities (x10,x20,…,xm0). Then, the initial partition of regional customers, which is denoted by Λ0={Λ10,Λ20,…,Λm0}, is generated by the spirit of NCR heuristic (Step 3 in Algorithm 2). Set k=1.
Step
1 (location phase). Solve the involved CSLP (11) and find the location of facility xik for Λik-1 by the variational inequality approach. Denote Xk=(x1kT,…,xmkT)T.
Step
2 (allocation phase). Update the partition of regional customers Λ from Λik-1 to Λik based on the spirit of NCR heuristic.
Step
3 If ∥Xk-Xk-1∥<ε, the current locations and partition are heuristic locations of facilities and heuristic partition of customers. Otherwise, set k=k+1 and go to Step 1.Remark 10.
At the(k+1)th iteration, it is recommended to use xik (and the corresponding y and z) in Step 1 as the initial iterate in the variational inequality approach for solving xik+1 (i=1,…,m), considering the fact that Λik usually differs from Λik-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 Section6.2. Hence, the convergence of the proposed algorithm is accelerated greatly by this initialization strategy.


## 5. Convergence of the Proposed Heuristic Algorithm

In this section, we analyze the convergence of the proposed location-allocation heuristic (Algorithm9). For simplification of our discussion, some notations are introduced as follows. Let A=A1×⋯×An, and recall Π=Π1×⋯×Πm. For any X=(x1T,…,xmT)T∈Π and C=(c1T,…,cnT)T∈A, we can define an ordered pair (X,C) and we can also define the function ω(X,C), in the current partition of customers as follows:
(43)ω(X,C)=∑i=1m∑j=1nwijγ(xi-cj)=∑i=1m∑{j∈𝒩:Aj∈Λi}sjγ(xi-cj),
which represents the objective functional value of CMLP (9) at (X,C).During the implementation of Algorithm9, we denote the map ℒ:Π×A→Π×A as the location operation in the Step 1 and the map 𝒜:Π×A→Π×A as the allocation operation in the Step 2. It follows that
(44)ℒ(Xk,Ck)=(Xk+1,C~k),𝒜(Xk+1,C~k)=(Xk+1,Ck+1),
where Xk and Xk+1 are, respectively, the locations of facilities in the kth and (k+1)th iteration, Ck is the variable of the closest points to Xk in the partition of Λk, C~k is the closest points to the new Xk+1 in the partition of Λk, and Ck+1 is the closest points to the new Xk+1 in the new partition of Λk+1. Then, the iterate scheme of the location-allocation heuristic is
(45)(Xk+1,Ck+1)=𝒜ℒ(Xk,Ck).LetS(X0,C0) denote the iterative sequence generated by the location-allocation heuristic for CMLP with the initial iterate (X0,C0). 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 S(X0,C0).Proposition 12.
S
(
X
0
,
C
0
) is strictly monotone in the sense that ω(Xk+1,Ck+1)<ω(Xk,Ck) if Xk+1≠Xk.Proof.
SinceXk+1≠Xk, there exists at least one i∈ℳ such that xik+1≠xik. For such i’s according to the following convex optimization problem
(46)xik+1=argminx∈Πi{Cik+1(x)=∑{j∈𝒩:Aj∈Λik}‍sjγ(x-aj(x))},
we have
(47)ω(Xk+1,C~k)≤ω(Xk,Ck).
Based on Remark 10, we know that if xik is the solution of (46), then xik+1 will be equal to xik. Therefore, xik+1≠xik implies that xik is not the solution of (46), and thus, Cik+1(xk+1)<Cik(xk). It follows that
(48)ω(Xk+1,C~k)<ω(Xk,Ck).
On the other hand, based on the principle of NCR in the allocation phase of Algorithm 9, we also have
(49)ω(Xk+1,Ck+1)≤ω(Xk+1,C~k).
By Combining (48) and (49), it follows that
(50)ω(Xk+1,Ck+1)<ω(Xk,Ck).
The proof is complete.Based on the monotonicity of the generated sequenceS(X0,C0), the following theorem can be proved.Theorem 13.
Let{Sk}:=S(X0,C0). Then, the generated sequence {Sk} satisfies that(1)
ω(Sk)→ω(S) for some S∈Π×A,(2)
all accumulation points of{Sk} have the same objective functional values.Proof.
After a finite number of iterations, ifXJ+1=XJ, then the iterates after SJ will be constant, and thus, {Sk} is convergent to (XJ,CJ)∈Π×A, and the two assertions are both true.
Below, we will discuss the case thatXk+1≠Xk for any k∈N. First, we prove the first assertion. Since Sk∈Π×A and Π×A is a compact space, it follows from the Bolzano-Weierstrass theorem that there exists a subsequence of {Sk} which, say {Sk}K, converges to an element S∈Π×A; that is,
(51)limK→∞{Sk}K⟶S,S∈Π×A.
Note that ω(X,C) is a continuous function according to (43); then,
(52)limK→∞ω({Sk}K)⟶ω(S),S∈Π×A.
Due to that {ω(Sk)} is a monotone sequence (Proposition 12) and has lower bound, then {ω(Sk)} is convergent. Thus, any subsequence of {ω(Sk)} will be convergent to the same value. Note that {ω({Sk}K)} is a subsequence of {ω(Sk)} and it is convergent to ω(S); then, it follows that
(53)limk→∞ω(Sk)⟶ω(S),S∈Π×A.
The second assertion can easily be proved. LetP be an accumulation point of {Sk}; then, there exists one subsequence {Sk}K′ which converges to P, and due to the continuity of ω(X,C), we have
(54)ω({Sk}K′)⟶ω(P).
The first assertion has shown that {ω(Sk)} is convergent, and note that {ω({Sk}K′)} is a subsequence of {ω(Sk)}; then, it follows that
(55)ω(P)=limK′→∞ω({Sk}K′)=limk→∞ω(Sk).
Thus, all accumulation points of {Sk} have the same objective functional values equal to limk→∞ω(Sk).Lemma 14.
ℒ
:
Π
×
A
→
Π
×
A which is defined in (44) is a closed map over Π×A.Proof.
Note that the CSLP (11) is a convex problem, then
(56)argminx∈Πi{∑{j∈𝒩:Aj∈Λik-1}‍sjminqj∈Ajγ(xi-qj)},i=1,2,…,m,
are continuous. Since Π×A is a compact space and also a Hausdorff space and every continuous map from a compact space to a Hausdorff space is closed, it follows that ℒ is closed over Π×A.Lemma 15.
Letυ0 be a given vector in Π×A and Δ:={υ∈Π×A∣ω(υ)≤ω(υ0)}. Then, Δ is a compact set.Proof.
It is known that every closed subset of a compact space is also compact, and therefore, it is enough to prove thatΔ is a closed set.
For any sequence{υk} with υk∈Δ, since Π×A is compact, according to Bolzano-Weierstrass, there exists a convergent subsequence {υk}K of {υk} such that
(57)limK→∞{υk}K⟶υ.
Due to the continuity of ω, it follows that
(58)limK→∞ω({υk}K)=ω(υ).
On the other hand,{υk}K∈Δ implies
(59)ω({υk}K)≤ω(υ0).
By combining (58), (59), and the continuity of ω, the following inequality is obtained:
(60)ω(υ)≤ω(υ0),
and accordingly, υ∈Δ. This means that Δ is closed, and the proof is complete.Now, we are ready to prove the convergence of the proposed location-allocation heuristic (Algorithm9). Let Ξ⊆Π×A be the nonempty local solution set of CMLP (9). Recall that in the location-allocation Cooper algorithm for MWP, if XJ+1=XJ occurs after a finite number of iterations, the iterates after XJ will be constant. Then, no further improvement is possible for MWP, and it follows from [32, 34] that the XJ is a local solution of MWP. Similarly, in the proposed location-allocation heuristic algorithm for CMLP, if XJ+1=XJ occurs, the iterates after SJ=(XJ,CJ) will also be constant. Then, exactly as in the location-allocation Cooper algorithm for MWP, no further improvement is possible for CMWP, and a local solution, namely, SJ, is obtained. Hence in this case, {Sk} is convergent to the SJ∈Ξ. However, it is not assured that XJ+1=XJ always occurs during the implementation of Algorithm 9, and therefore, we assume that Xk+1≠Xk for any k∈N and prove the convergence in this case.Theorem 16.
Assume thatXk+1≠Xk for any k∈N; then, all the accumulation points of the sequence {Sk} belong to Ξ.Proof.
LetS be an accumulation point of {Sk}. Due to Sk∈Π×A and the compactness of Π×A, we know that (1)  S∈Π×A and (2) there exists a subsequence {Sk}K which is convergent to S. According to Theorem 13, we know that
(61)limk→∞ω(Sk)=ω(S),S∈Π×A.
So, it is enough to prove S∈Ξ. We prove this by contradiction. Assume that S∉Ξ; that is, S is not a solution, and we consider the subsequence {Sk+1}K. Denote Δ={υ∈Π×A∣ω(υ)≤ω(S0)}. Due to Proposition 12, it follows that for all k∈N we have (Xk+1,C~k)∈Δ and (Xk+1,Ck+1)∈Δ. According to Lemma 15, Δ is compact, and thus, there exists K′⊂K such that
(62)limK′→∞(Xk+1,C~k)K′=υ1,limK′→∞(Xk+1,Ck+1)K′=υ2.
According to Lemma 14, the map ℒ is closed at S∈Π×A; then, it follows that υ1=ℒ(S). Further, due to S∉Ξ, v1 will be not equal to S. Otherwise, we can choose S 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., XJ+1=XJ) that S will be a local solution, which contradicts with S∉Ξ. Therefore, we can obtain v1≠S. Together with the monotone proposition of ℒ (48), we have the inequality
(63)ω(υ1)<ω(S).
On the other hand, note that𝒜(Xk+1,C~k)=(Xk+1,Ck+1); then, by the monotonicity of 𝒜 (49), it follows that
(64)ω(Xk+1,Ck+1)K′≤ω(Xk+1,C~k)K′,
and thus, by taking the limit for (64) and by the continuity of ω, we know ω(υ2)≤ω(υ1). Combining this with (63), we obtain
(65)ω(υ2)<ω(S).
However, note that υ2 and S are two accumulation points of {Sk}, and according to the second assertion of Theorem 13, ω(υ2)=ω(S), which will contradict with (65). Therefore, our assumption is wrong, and thus, S∈Ξ.As a result, we have the following convergence theorem for the sequence generated by the proposed location-allocation algorithm.Theorem 17.
The sequenceS(X0,C0) generated by the proposed location-allocation heuristic algorithm either converges to a point in Ξ or all accumulation points of S(X0,C0) belong to Ξ.


## 6. Numerical Results

This section reports some preliminary numerical results to verify the theoretical assertions proved in previous sections. Section3.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).


### 6.1. Numerical Results of Variational Inequality Approach for CSLP

When applying the variational inequality approach for solving CSLP and MCSLP (12), theoretically, the initial iteration u0 in Algorithm 6 can be chosen arbitrarily in Ω. In practical implementation, however, we choose u0 judiciously similar to the initialization strategy in the location-allocation heuristic: let gj(j=1,…,d) be the centers of regional customers, solve the following single-source Weber problem (SWP):
(66)x*=argminx∈Π{∑j=1d‍γ(x-gj)}
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].Example 18.
Here,d=5; that is, there are five regional customers, and all customers are unit squares whose sides are parallel to the axes. The geometric centers of the five customers are (0.5,0.5), (4.5,0.5), (0.5,2.5), (2.5,2.5), and (4.5,2.5), and sj=1, j=1,…,5.In order to clarify the comparison of two methods, we choose the same stopping criterion as∥xk+1-xk∥≤10-4 (throughout this section, ∥·∥ is the l∞-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.Table 1
Numerical results for Example1 given in [16].
Main results
VI approach
Weiszfeld-type method
x
(2.5000, 1.9484)
(2.5000, 1.9484)
a
1
(
x
)
(1.0000, 1.0000)
(1.0000, 1.0000)
a
2
(
x
)
(4.0000, 1.0000)
(4.0000, 1.0000)
a
3
(
x
)
(4.0000, 2.0000)
(4.0000, 2.0000)
a
4
(
x
)
(2.5000, 2.0000)
(2.5000, 2.0000)
a
5
(
x
)
(1.0000, 2.0000)
(1.0000, 2.0000)
Objective value
6.6027
6.6027
Number of iterations
32.81/23.77
237.60
CPU
0.0124
0.0610
“32.81/23.77” means 32.81 Iter. for initialization and 23.77 Iter. for VI approach.According to Table1, 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 Example18, the number of customers d=5, and all customers are unit squares whose edges are parallel to the axes. Let d1 be the distance between any two neighboring customers, and we set d1 equal to 1, 0.1, 0.01, and 0.001, respectively. The weights sj(j=1,…,5) are randomly generated in the area of regional customers.Note that in Example19 the parameter d1 reflects how close the customers are away from its neighborhoods. We test this example for a large number of times with the stopping criterion ∥xk+1-xk∥≤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.Table 2
Numerical results of VI approach and Weiszfeld-type method for Example19.
d
1
VI approach
Weiszfeld-type method
Impro.
Freq.
No. of Iter.0
No. of Iter.
CPU
Obj.
No. of Iter.
CPU
Obj.
percent
num.
1
25.11
59.96
0.015
15.4465
52.99
0.0117
15.5664
0.78
90
25.04
58.97
0.015
15.1803
27.28
0.0055
15.2926
0.74
92
25.35
60.89
0.015
15.3693
31.11
0.0070
15.4834
0.74
92
25.24
59.80
0.015
15.0924
35.45
0.0080
15.1961
0.69
90
23.61
55.13
0.014
14.6321
34.27
0.0069
14.7055
0.50
93
0.1
32.95
278.26
0.060
4.2157
28.59
0.0082
4.6791
10.99
63
32.11
270.81
0.059
4.0864
60.38
0.0135
4.4318
8.45
66
34.07
301.67
0.067
4.2248
36.63
0.0081
4.6001
8.88
63
34.16
278.24
0.065
4.1189
34.03
0.0067
4.4854
8.89
63
32.63
280.00
0.062
4.0842
38.20
0.0069
4.3619
6.80
64
0.01
33.11
427.34
0.089
3.0633
25.85
0.0063
3.8005
24.07
28
35.28
407.61
0.086
3.0615
39.92
0.0087
3.7669
23.04
25
39.26
382.82
0.080
3.1259
26.53
0.0064
3.6954
18.22
29
35.46
425.33
0.089
2.9975
35.06
0.0077
3.6500
21.77
35
32.96
437.24
0.091
3.0691
57.90
0.0130
3.6925
20.31
33
0.001
35.44
151.52
0.035
3.0477
18.88
0.0045
3.8310
25.70
15
36.87
140.58
0.032
3.0224
17.20
0.0048
3.7128
22.84
14
33.43
91.72
0.022
2.9797
17.53
0.0042
3.9617
32.96
17
33.79
115.79
0.027
2.9600
15.20
0.0044
3.7221
25.75
19
37.19
96.08
0.022
3.1936
16.47
0.0047
4.1362
29.52
18According to Table2, 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 d1 decreases, that is, the customers become closer and closer, the frequency that Weiszfeld-type method obtains the optimal solution gets smaller and smaller. When d1 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 d1=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 d1 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 customersd 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 ∥x-O∥≤r, where the center O is randomly generated in [-100,100]2, and the radius r is randomly generated in (1,5); the stopping criterion of VI approach is chosen as
(67)∥xk+1-xk∥≤10-4,∥e(uk+1)∥≤10-4,
and the initial iterate is randomly generated in [-100,100]2; the gauge γ(·) is generated with the unit ball set as
(68)9(x+23)2+12y2=16.
For each d, 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.Table 3
Numerical results of VI approach for CSLP (12).
d
VI approach with Initial.
VI approach without Initial.
No. of Iter.0
No. of Iter.
CPU
No. of Iter.
CPU
10
26.59
20.30
0.0108
25.26
0.0080
20
24.51
32.80
0.0239
38.20
0.0193
50
27.25
57.69
0.0847
64.89
0.0762
100
25.09
92.53
0.2664
106.84
0.2442
200
29.83
145.97
1.2196
170.66
1.3389
500
27.97
224.74
13.7783
252.99
15.9382
1000
26.71
259.03
56.8598
284.54
62.0384
2000
23.12
281.59
245.5680
346.52
295.2305According to Table3, 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 d is small, for example, d=10,20,50,and100, which can be explained as follows. When the number of customers d 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 d 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 d=2000.


### 6.2. Numerical Results of Heuristic Algorithm for CMLP

This subsection applies the proposed location-allocation heuristic algorithm (Algorithm9) 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 [-100100]2, and all the demands, sj(j=1,2,…,n), are randomly generated in [110]. The gauge γ(·) is also defined with the unit ball set as (68). We test the scenario with n=100,200,500,and1000 and m=2,4,6,8,and10; the locational constraints are ∥x-Oj∥≤rj(j=1,…,m), where the radius rj is randomly generated in [1,10], and the center Oj is given in advance as follows:
(69)m=2:O1=(-50,0)T,O2=(50,0)T;m=4:O1=(-50,-50)T,O2=(50,-50)T,vvvvvvviO3=(-50,50)T,O4=(50,50)T;m=6:O1=(-50,-50)T,O2=(0,-50)T,vvvvvvviO3=(50,-50)T,O4=(-50,50)T,vvvvvvviO5=(0,50)T,O6=(50,50)T;m=8:O1=(-50,-50)T,O2=(0,-50)T,vvvvvvviO3=(50,-50)T,O4=(-50,50)T,vvvvvvviO5=(0,50)T,O6=(50,50)T,vvvvvvviO7=(-25,0)T,O8=(25,0)T;m=10:O1=(-50,-50)T,O2=(0,-50)T,vvvvvvviO3=(50,-50)T,O4=(-50,50)T,vvvvvvviO5=(0,50)T,O6=(50,50)T,vvvvvvviO7=(-60,0)T,O8=(-20,0)T,vvvvvvviO9=(20,0)T,O10=(60,0)T.
The initial locations of facilities are randomly generated in [-100100]2, and the stopping criterion used in Algorithms 6 and 9 is chosen as
(70)∥xk+1-xk∥<10-4.To show the significance of the initialization strategy, we compare the numerical performance of the location-allocation heuristic algorithm with initialization strategy (denoted by “Algorithm9 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 “Iter0”) and the average iteration numbers of PC method in one iteration of the algorithm (denoted by “Iter0-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 (n,m), we test the CMLP for 100 times, and the computational performance is reported in Table 4.Table 4
Numerical results of location-allocation heuristic for CMLP.
n
m
Algorithm9 with Initial.
Algorithm9 without Initial.
Iter.0
Iter.0-PC
Iter.
Iter.-PC
CPU
Iter.
Iter.-PC
CPU
100
2
2.83
73.33
1.20
289.00
0.4868
3.04
176.70
0.7018
4
3.61
99.88
1.29
176.80
0.2652
2.65
218.47
0.4992
6
4.55
99.05
1.53
121.75
0.1710
3.50
130.96
0.2580
8
3.60
126.23
1.82
149.37
0.2250
3.02
170.32
0.2836
10
3.84
138.61
2.26
257.51
0.3212
4.65
240.48
0.5118
200
2
3.24
81.22
1.24
1620.80
5.4664
3.49
1192.33
11.4596
4
4.01
112.96
1.71
447.08
1.2932
4.37
433.10
2.5694
6
4.65
141.10
2.23
262.97
0.7958
3.88
305.75
1.4538
8
4.26
130.25
2.65
246.80
0.6268
4.83
287.12
1.1734
10
4.67
171.65
2.68
281.40
0.5866
4.46
346.18
1.0232
500
2
4.03
82.35
1.63
1998.20
64.5468
4.26
1815.90
139.2036
4
4.66
104.28
2.27
1304.48
12.3492
5.25
1445.82
26.5510
6
5.27
133.55
2.62
727.13
5.3322
7.60
875.26
18.2956
8
7.60
169.97
2.84
549.60
2.8050
5.24
429.83
3.8626
10
5.42
187.95
3.47
508.13
2.3742
5.02
507.70
3.3946
1000
2
4.49
56.69
1.91
718.17
93.2321
5.63
1587.31
365.6085
4
4.81
80.60
2.10
1376.13
46.2636
5.49
1576.85
113.8830
6
6.84
96.75
3.43
683.35
16.1398
7.65
1344.53
57.7762
8
7.20
161.50
3.28
563.69
9.3570
6.06
760.93
20.0896
10
6.63
171.89
4.05
615.07
7.5880
6.47
700.53
12.9448It follows from Table4 that the proposed location-allocation 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 Table4 deserves further illustration. Recall that with the number of new facilities (m) increasing, the number of subproblems (CSLPs) in location phase increases too. According to Table 4, however, we find that for fixed number of customers (n), with m increasing, the number of iterations for solving m CSLPs in one iteration of Algorithm 9 does not increase but almost decreases with m, especially for large-scale CMLP. This can be illustrated roughly as follows. For fixed n, when m increases, the average number of customers in each Λik 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 m CSLPs decreases. Similarly, due to the same reason, the computing time for solving CMLP also decreases with m increasing, as reported in the column of “CPU” in Table 4.


## 6.1. Numerical Results of Variational Inequality Approach for CSLP

When applying the variational inequality approach for solving CSLP and MCSLP (12), theoretically, the initial iteration u0 in Algorithm 6 can be chosen arbitrarily in Ω. In practical implementation, however, we choose u0 judiciously similar to the initialization strategy in the location-allocation heuristic: let gj(j=1,…,d) be the centers of regional customers, solve the following single-source Weber problem (SWP):
(66)x*=argminx∈Π{∑j=1d‍γ(x-gj)}
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].Example 18.
Here,d=5; that is, there are five regional customers, and all customers are unit squares whose sides are parallel to the axes. The geometric centers of the five customers are (0.5,0.5), (4.5,0.5), (0.5,2.5), (2.5,2.5), and (4.5,2.5), and sj=1, j=1,…,5.In order to clarify the comparison of two methods, we choose the same stopping criterion as∥xk+1-xk∥≤10-4 (throughout this section, ∥·∥ is the l∞-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.Table 1
Numerical results for Example1 given in [16].
Main results
VI approach
Weiszfeld-type method
x
(2.5000, 1.9484)
(2.5000, 1.9484)
a
1
(
x
)
(1.0000, 1.0000)
(1.0000, 1.0000)
a
2
(
x
)
(4.0000, 1.0000)
(4.0000, 1.0000)
a
3
(
x
)
(4.0000, 2.0000)
(4.0000, 2.0000)
a
4
(
x
)
(2.5000, 2.0000)
(2.5000, 2.0000)
a
5
(
x
)
(1.0000, 2.0000)
(1.0000, 2.0000)
Objective value
6.6027
6.6027
Number of iterations
32.81/23.77
237.60
CPU
0.0124
0.0610
“32.81/23.77” means 32.81 Iter. for initialization and 23.77 Iter. for VI approach.According to Table1, 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 Example18, the number of customers d=5, and all customers are unit squares whose edges are parallel to the axes. Let d1 be the distance between any two neighboring customers, and we set d1 equal to 1, 0.1, 0.01, and 0.001, respectively. The weights sj(j=1,…,5) are randomly generated in the area of regional customers.Note that in Example19 the parameter d1 reflects how close the customers are away from its neighborhoods. We test this example for a large number of times with the stopping criterion ∥xk+1-xk∥≤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.Table 2
Numerical results of VI approach and Weiszfeld-type method for Example19.
d
1
VI approach
Weiszfeld-type method
Impro.
Freq.
No. of Iter.0
No. of Iter.
CPU
Obj.
No. of Iter.
CPU
Obj.
percent
num.
1
25.11
59.96
0.015
15.4465
52.99
0.0117
15.5664
0.78
90
25.04
58.97
0.015
15.1803
27.28
0.0055
15.2926
0.74
92
25.35
60.89
0.015
15.3693
31.11
0.0070
15.4834
0.74
92
25.24
59.80
0.015
15.0924
35.45
0.0080
15.1961
0.69
90
23.61
55.13
0.014
14.6321
34.27
0.0069
14.7055
0.50
93
0.1
32.95
278.26
0.060
4.2157
28.59
0.0082
4.6791
10.99
63
32.11
270.81
0.059
4.0864
60.38
0.0135
4.4318
8.45
66
34.07
301.67
0.067
4.2248
36.63
0.0081
4.6001
8.88
63
34.16
278.24
0.065
4.1189
34.03
0.0067
4.4854
8.89
63
32.63
280.00
0.062
4.0842
38.20
0.0069
4.3619
6.80
64
0.01
33.11
427.34
0.089
3.0633
25.85
0.0063
3.8005
24.07
28
35.28
407.61
0.086
3.0615
39.92
0.0087
3.7669
23.04
25
39.26
382.82
0.080
3.1259
26.53
0.0064
3.6954
18.22
29
35.46
425.33
0.089
2.9975
35.06
0.0077
3.6500
21.77
35
32.96
437.24
0.091
3.0691
57.90
0.0130
3.6925
20.31
33
0.001
35.44
151.52
0.035
3.0477
18.88
0.0045
3.8310
25.70
15
36.87
140.58
0.032
3.0224
17.20
0.0048
3.7128
22.84
14
33.43
91.72
0.022
2.9797
17.53
0.0042
3.9617
32.96
17
33.79
115.79
0.027
2.9600
15.20
0.0044
3.7221
25.75
19
37.19
96.08
0.022
3.1936
16.47
0.0047
4.1362
29.52
18According to Table2, 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 d1 decreases, that is, the customers become closer and closer, the frequency that Weiszfeld-type method obtains the optimal solution gets smaller and smaller. When d1 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 d1=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 d1 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 customersd 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 ∥x-O∥≤r, where the center O is randomly generated in [-100,100]2, and the radius r is randomly generated in (1,5); the stopping criterion of VI approach is chosen as
(67)∥xk+1-xk∥≤10-4,∥e(uk+1)∥≤10-4,
and the initial iterate is randomly generated in [-100,100]2; the gauge γ(·) is generated with the unit ball set as
(68)9(x+23)2+12y2=16.
For each d, 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.Table 3
Numerical results of VI approach for CSLP (12).
d
VI approach with Initial.
VI approach without Initial.
No. of Iter.0
No. of Iter.
CPU
No. of Iter.
CPU
10
26.59
20.30
0.0108
25.26
0.0080
20
24.51
32.80
0.0239
38.20
0.0193
50
27.25
57.69
0.0847
64.89
0.0762
100
25.09
92.53
0.2664
106.84
0.2442
200
29.83
145.97
1.2196
170.66
1.3389
500
27.97
224.74
13.7783
252.99
15.9382
1000
26.71
259.03
56.8598
284.54
62.0384
2000
23.12
281.59
245.5680
346.52
295.2305According to Table3, 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 d is small, for example, d=10,20,50,and100, which can be explained as follows. When the number of customers d 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 d 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 d=2000.


## 6.2. Numerical Results of Heuristic Algorithm for CMLP

This subsection applies the proposed location-allocation heuristic algorithm (Algorithm9) 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 [-100100]2, and all the demands, sj(j=1,2,…,n), are randomly generated in [110]. The gauge γ(·) is also defined with the unit ball set as (68). We test the scenario with n=100,200,500,and1000 and m=2,4,6,8,and10; the locational constraints are ∥x-Oj∥≤rj(j=1,…,m), where the radius rj is randomly generated in [1,10], and the center Oj is given in advance as follows:
(69)m=2:O1=(-50,0)T,O2=(50,0)T;m=4:O1=(-50,-50)T,O2=(50,-50)T,vvvvvvviO3=(-50,50)T,O4=(50,50)T;m=6:O1=(-50,-50)T,O2=(0,-50)T,vvvvvvviO3=(50,-50)T,O4=(-50,50)T,vvvvvvviO5=(0,50)T,O6=(50,50)T;m=8:O1=(-50,-50)T,O2=(0,-50)T,vvvvvvviO3=(50,-50)T,O4=(-50,50)T,vvvvvvviO5=(0,50)T,O6=(50,50)T,vvvvvvviO7=(-25,0)T,O8=(25,0)T;m=10:O1=(-50,-50)T,O2=(0,-50)T,vvvvvvviO3=(50,-50)T,O4=(-50,50)T,vvvvvvviO5=(0,50)T,O6=(50,50)T,vvvvvvviO7=(-60,0)T,O8=(-20,0)T,vvvvvvviO9=(20,0)T,O10=(60,0)T.
The initial locations of facilities are randomly generated in [-100100]2, and the stopping criterion used in Algorithms 6 and 9 is chosen as
(70)∥xk+1-xk∥<10-4.To show the significance of the initialization strategy, we compare the numerical performance of the location-allocation heuristic algorithm with initialization strategy (denoted by “Algorithm9 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 “Iter0”) and the average iteration numbers of PC method in one iteration of the algorithm (denoted by “Iter0-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 (n,m), we test the CMLP for 100 times, and the computational performance is reported in Table 4.Table 4
Numerical results of location-allocation heuristic for CMLP.
n
m
Algorithm9 with Initial.
Algorithm9 without Initial.
Iter.0
Iter.0-PC
Iter.
Iter.-PC
CPU
Iter.
Iter.-PC
CPU
100
2
2.83
73.33
1.20
289.00
0.4868
3.04
176.70
0.7018
4
3.61
99.88
1.29
176.80
0.2652
2.65
218.47
0.4992
6
4.55
99.05
1.53
121.75
0.1710
3.50
130.96
0.2580
8
3.60
126.23
1.82
149.37
0.2250
3.02
170.32
0.2836
10
3.84
138.61
2.26
257.51
0.3212
4.65
240.48
0.5118
200
2
3.24
81.22
1.24
1620.80
5.4664
3.49
1192.33
11.4596
4
4.01
112.96
1.71
447.08
1.2932
4.37
433.10
2.5694
6
4.65
141.10
2.23
262.97
0.7958
3.88
305.75
1.4538
8
4.26
130.25
2.65
246.80
0.6268
4.83
287.12
1.1734
10
4.67
171.65
2.68
281.40
0.5866
4.46
346.18
1.0232
500
2
4.03
82.35
1.63
1998.20
64.5468
4.26
1815.90
139.2036
4
4.66
104.28
2.27
1304.48
12.3492
5.25
1445.82
26.5510
6
5.27
133.55
2.62
727.13
5.3322
7.60
875.26
18.2956
8
7.60
169.97
2.84
549.60
2.8050
5.24
429.83
3.8626
10
5.42
187.95
3.47
508.13
2.3742
5.02
507.70
3.3946
1000
2
4.49
56.69
1.91
718.17
93.2321
5.63
1587.31
365.6085
4
4.81
80.60
2.10
1376.13
46.2636
5.49
1576.85
113.8830
6
6.84
96.75
3.43
683.35
16.1398
7.65
1344.53
57.7762
8
7.20
161.50
3.28
563.69
9.3570
6.06
760.93
20.0896
10
6.63
171.89
4.05
615.07
7.5880
6.47
700.53
12.9448It follows from Table4 that the proposed location-allocation 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 Table4 deserves further illustration. Recall that with the number of new facilities (m) increasing, the number of subproblems (CSLPs) in location phase increases too. According to Table 4, however, we find that for fixed number of customers (n), with m increasing, the number of iterations for solving m CSLPs in one iteration of Algorithm 9 does not increase but almost decreases with m, especially for large-scale CMLP. This can be illustrated roughly as follows. For fixed n, when m increases, the average number of customers in each Λik 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 m CSLPs decreases. Similarly, due to the same reason, the computing time for solving CMLP also decreases with m increasing, as reported in the column of “CPU” in Table 4.


## 7. Conclusion

In this paper, we are interested in the locations of multiple facilities in the spaceRp 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.



---

*Source: 624398-2013-10-27.xml*