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

Abstract

This paper considers the locations of multiple facilities in the space 𝑅^𝑝, 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.

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 as 𝑙_𝑝-norms and block norms. References [3, 18] discuss the approximations about the weighted 𝑙_𝑝-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 space 𝑅^𝑝 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 𝑅^𝑝 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. Section 2 states the formulation of our problem, which is shown to be nonconvex and NP-hard. The spirit of the well-known location-allocation heuristic algorithm, which consists of a location phase and an allocation phase in each iteration, is also discussed in this section. In Section 3, the subproblems arising in location phase and allocation phase are solved. More specifically, for the subproblem arising in allocation phase, the regional customers are allocated to facilities according to the nearest reclassification (NCR) heuristic, whereas for the single-source subproblem arising in location phase, the relationship between the subproblem and a monotone linear variational inequality (LVI) is firstly built, and then, a projection-contraction (PC) method is adopted to find the optimal location of the facility. A new 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 space 𝑅^𝑝 for 𝑚 facilities, with the objective to minimize the sum of weighted closest distances between these facilities and 𝑛 regional customers. Note that the minisum single-source models with closest distance have been well justified in [6, 16] where the distances are particularly measured by 𝑙₁-norm and 𝑙₂-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 𝑙₁-norm used in [16] and the 𝑙₂-norm used in [6]) and the asymmetric distance are considered.

Let 𝐵 be a compact convex set in 𝑅^𝑝, and let the interior of 𝐵 contain the origin; then, the gauge of 𝐵 is defined by 𝛾(𝑥)𝑥=inf𝜆>0∣𝜆∈𝐵,∀𝑥∈𝑅^𝑝. (1)𝐵is called the unit ball of 𝛾(⋅) due to 𝐵=𝑥∈𝑅^𝑝∣𝛾(𝑥)≤1. (2) This way to define a gauge from a compact convex set was first introduced by Minkowski [29]. The gauge 𝛾(⋅) satisfies the following properties: (1)𝛾(𝑥)≥0 and 𝛾(𝑥)=0⇔𝑥=0; (2)𝛾(𝑡𝑥)=𝑡𝛾(𝑥) for any 𝑡≥0; (3)𝛾(𝑥+𝑦)≤𝛾(𝑥)+𝛾(𝑦) for any 𝑥 and 𝑦∈𝑅^𝑝.It follows from (2) and (3) that any gauge 𝛾(𝑥) is a convex function of 𝑥. The distance measuring function can be derived from a gauge by 𝐷𝑥,𝑦=𝛾𝑦−𝑥. (3) When 𝐵 is symmetric around the origin, according to the definition of (1), we have 𝛾(−𝑥)=𝛾(𝑥),∀𝑥∈𝑅^𝑝. (4) Combining (3) and (4), it follows that 𝐷𝑥,𝑦=𝐷𝑦,𝑥, (5) which means that the distance measured by 𝛾(⋅) is symmetric. On the contrary, when 𝐵 is not symmetric around the origin, (4) does not hold any more, and thus, the distance measured by the gauge 𝛾(⋅) is asymmetric. Thus, when different compact convex sets are used as unit balls, different gauges (symmetric or asymmetric) can be generated and employed to measure distances in location problems, which depends on the requirements of practical applications.

Let Λ={𝐴_𝑗⊂𝑅^𝑝∣𝑗=1,…,𝑛} denote the set of regional customers, and let 𝑥_𝑖(𝑖=1,…,𝑚) be the location of the 𝑖th facility. Each regional customer 𝐴_𝑗 is simply assumed to be closed and convex. We denote the closest point in 𝐴_𝑗 to the facility 𝑥_𝑖 by 𝑎_𝑗𝑥_𝑖∶=argmin𝐷𝑞,𝑥_𝑖∣𝑞∈𝐴_𝑗=argmin𝛾𝑥_𝑖−𝑞∣𝑞∈𝐴_𝑗. (6) Then, the closest distance between the facility 𝑥_𝑖 and the regional customer 𝐴_𝑗, denoted by 𝑑_𝑗(𝑥_𝑖), can be represented by 𝑑_𝑗𝑥_𝑖∶=min_{𝑞∈𝐴𝑗}𝐷𝑞,𝑥_𝑖=𝛾𝑥_𝑖−𝑎_𝑗𝑥_𝑖. (7)

When the gauge is used to measure distances, we have the following proposition for 𝑎_𝑗(𝑥) and 𝑑_𝑗(𝑥) which is similar to that in [16].

Proposition 1. Let 𝑥 be the location of a facility; then, the closest point 𝑎_𝑗(𝑥) in (6) is well defined, and the closest distance 𝑑_𝑗(𝑥) in (7) is a convex function of 𝑥.

Proof. Since 𝐴_𝑗 is a convex set and 𝛾(⋅) is a convex function, (6) is a convex problem, and thus, 𝑎_𝑗(𝑥) is well defined.
Now, we prove that 𝑑_𝑗(𝑥) is a convex function of 𝑥 as follows. Let 𝑥 and 𝑦 be two points in 𝑅^𝑝 and 𝜆∈[0,1]; then, due to 𝑎_𝑗(𝑥) and 𝑎_𝑗(𝑦) in 𝐴_𝑗 and the convexity of 𝐴_𝑗, it follows that 𝜆𝑎_𝑗(𝑥)+(1−𝜆)𝑎_𝑗(𝑦)∈𝐴_𝑗, and thus, we have 𝑑_𝑗𝜆𝑥+(1−𝜆)𝑦=min_{𝑎∈𝐴𝑗}𝐷𝑎,𝜆𝑥+(1−𝜆)𝑦≤𝐷𝜆𝑎_𝑗(𝑥)+(1−𝜆)𝑎_𝑗𝑦()𝑦(),𝜆𝑥+1−𝜆=𝛾𝜆𝑥+1−𝜆𝑦−𝜆𝑎_𝑗(𝑥)+(1−𝜆)𝑎_𝑗𝑦=𝛾𝜆𝑥−𝑎_𝑗(𝑥)+(1−𝜆)𝑦−𝑎_𝑗𝑦≤𝛾𝜆𝑥−𝑎_𝑗(𝑥)+𝛾(1−𝜆)𝑦−𝑎_𝑗𝑦=𝜆𝛾𝑥−𝑎_𝑗(𝑥)+(1−𝜆)𝛾𝑦−𝑎_𝑗𝑦=𝜆𝑑_𝑗(𝑥)+(1−𝜆)𝑑_𝑗𝑦. (8) Therefore, 𝑑_𝑗(𝑥) is convex with respect to 𝑥, and the proof is complete.

Based on the notations introduced above, now the constrained multi-source location problem (abbreviated as CMLP) we consider in this paper takes the following formulation: CMLP:min^𝑛_𝑋,𝑊^𝑚_𝑗=1𝑖=1𝑤_𝑖𝑗𝑑_𝑗𝑥_𝑖=^𝑛𝑚_𝑗=1𝑖=1𝑤_𝑖𝑗𝛾𝑥_𝑖−𝑎_𝑗𝑥_𝑖s.t.^𝑚_𝑖=1𝑤_𝑖𝑗=𝑠_𝑗𝑥,𝑗=1,2,…,𝑛,_𝑖∈Π_𝑖,𝑖=1,2,…,𝑚, (9) where 𝑠_𝑗≥0 is the given demand required by the 𝑗th customer, 𝑋=(𝑥^𝑇₁,…,𝑥^𝑇_𝑚)^𝑇 is the variable of the locations of facilities to be determined, 𝑊=(𝑤_𝑖𝑗)_𝑚×𝑛 is the undetermined variable of 𝑤_𝑖𝑗 which denotes the unknown allocation from the 𝑖th facility to the 𝑗th customer, and Π_𝑖 is the locational constraint for the 𝑖th facility which is assumed to be a convex and closed set in 𝑅^𝑝.

More explanations are required for our model (9). First, the locational constraint Π_𝑖 in (9) can also be chosen as 𝑅^𝑝, and if all Π_𝑖(𝑖=1,…,𝑚) are 𝑅^𝑝, the CMLP (9) is a unconstrained problem, and thus, both the constrained and unconstrained problems are considered in our model. Second, mark that the minisum models analyzed in [6, 16] are two special cases of CMLP (9), where 𝑚=1, Π_𝑖=𝑅^𝑝, and 𝛾(𝑥) is particularly 𝑙₁-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: CMLP^′:min_𝑋∈Π𝐶(𝑋)=^𝑛_𝑗=1𝑠_𝑗min_{1≤𝑖≤𝑚}𝑑_𝑗𝑥_𝑖=^𝑛_𝑗=1𝑠_𝑗min_{1≤𝑖≤𝑚}𝛾𝑥_𝑖−𝑎_𝑗𝑥_𝑖, (10) where Π is the cartesian product of locational constraints; that is, Π=Π₁×⋯×Π_𝑚.

When all 𝐴_𝑗 are points not regions and all Π_𝑖 are 𝑅^𝑝, CMLP (9) reduces to the well-known multi-source Weber problem (MWP) which has wide applications in operations research, marketing, urban planning, and so forth; see, for example, [3, 30, 31]. Recall the fact that the multi-source Weber problem is nonconvex [32] and NP-hard [33], and therefore, heuristics algorithms are extremely popular and highly appreciated for overcoming the difficulty caused by its nonconvexity and NP-hardness; see, for example, [3, 30, 34–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,…,𝑚}, 𝒩={1,2,…,𝑛}, and then Λ={𝐴_𝑗∶𝑗∈𝒩}. At the (𝑘−1)th iteration, let {Λ₁^𝑘−1,Λ₂^𝑘−1,…,Λ_𝑚^𝑘−1} be the disjoint partition of Λ in the sense that ∪^𝑚_𝑖=1Λ_𝑖^𝑘−1=Λ and Λ_𝑖^𝑘−1∩Λ_𝑗^𝑘−1=∅ (for 𝑖≠𝑗), and each Λ_𝑖^𝑘−1 (𝑖=1,…,𝑚) in the partition is called one cluster. Then, at the 𝑘th iteration, the location phase finds the candidates of locations of facilities (denoted by 𝑥^𝑘₁,𝑥^𝑘₂,…,𝑥^𝑘_𝑚) by solving the following 𝑚 constrained single-source location problems (CSLP) with the closest distance under gauge for each cluster Λ_𝑖^𝑘−1,  𝑖=1,…,𝑚: 𝑥CSLP:_𝑘𝑖=argmin_{𝑥∈Π𝑖}𝐶^𝑘_𝑖(𝑥)∶=_{𝑗∈𝒩∶𝐴𝑗∈Λ𝑖𝑘−1}𝑠_𝑗𝛾𝑥−𝑎_𝑗(𝑥)=_{𝑗∈𝒩∶𝐴𝑗∈Λ𝑖𝑘−1}𝑠_𝑗min_{𝑞𝑗∈𝐴𝑗}𝛾𝑥−𝑞_𝑗. (11) After the location phase, the allocation phase then revises the current partition of Λ to generate a new disjoint partition of Λ={Λ^𝑘₁,Λ^𝑘₂,…,Λ^𝑘_𝑚} by the following nearest center reclassification (NCR) heuristic (see [38]): for some customer 𝐴_𝑗∈Λ_ℎ^𝑘−1 (𝑗∈{1,2,…,𝑛} and ℎ∈{1,2,…,𝑚}), if 𝑥^𝑘_𝑙 (𝑙≠ℎ) is the nearest point for 𝐴_𝑗 among all 𝑥^𝑘_𝑖 computed by (11), then Λ^𝑘_ℎ=Λ_ℎ^𝑘−1⧵{𝐴_𝑗} and Λ^𝑘_𝑙=Λ_𝑙^𝑘−1∪{𝐴_𝑗}. If 𝑥^𝑘_𝑖 solved by (11) is the nearest facility for each regional customer in Λ_𝑖^𝑘−1 for any 𝑖∈{1,2,…,𝑚}, then 𝑥^𝑘_𝑖 (𝑖=1,2,…,𝑚) are the desirable locations of facilities and stop. Otherwise, we set 𝑘=𝑘+1 and repeat the iterations.

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

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 Λ⁰={Λ⁰₁,Λ⁰₂,…,Λ⁰_𝑚}.
For  𝑘=1,2,…, do.
Step 1. Set 𝑡=0 (𝑡 stores the number of reassignments);
Step 2. Compute the facility 𝑥^𝑘_𝑖 of Λ_𝑖^𝑘−1 by solving CSLP (11), for 𝑖=1,2,…,𝑚;
Step 3. For 𝑗=1,2,…,𝑛 do: 𝑑_𝑖𝑗∶=𝛾(𝑥^𝑘_𝑖−𝑎_𝑗(𝑥^𝑘_𝑖)) for 𝑖=1,…,𝑚;  if 𝐴_𝑗∈Λ_ℎ^𝑘−1 and 𝑑_𝑙𝑗=min_{𝑖=1,…,𝑚;𝑖≠ℎ}{𝑑_𝑖𝑗}<𝑑_ℎ𝑗,  then Λ^𝑘_ℎ=Λ_ℎ^𝑘−1⧵{𝐴_𝑗}, Λ^𝑘_𝑙=Λ_𝑙^𝑘−1∪{𝐴_𝑗};  𝑡=𝑡+1.
Step 4. If 𝑡=0, then the iteration terminates with {𝑥^𝑘₁,…,𝑥^𝑘_𝑚} being the desirable locations for facilities and the customers in Λ_𝑖^𝑘−1 being served by 𝑥^𝑘_𝑖 (𝑖=1,…,𝑚).

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 Π_𝑖=𝑅^𝑝 and the gauge 𝛾(⋅) are the particular 𝑙₁-norm in [16] and 𝑙₂-norm in [6]. For the model under 𝑙₁-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 𝑙₂-norm in [6], an efficient Weiszfeld-type method is proposed, and the convergence of this method is analyzed. Similar to Weiszfeld procedure [2], one problem of the proposed method is that the singular case, that is, the current iterate happens to be within some location of regional customers, may occur during its implementation. Due to the use of the gradient of objective function in the iteration, this method will terminate unexpectedly once the singular case occurs. In order to tackle the undesirable singular case and make the Weiszfeld-type method computational effective, the authors suggest to ignore the gradient of ‖𝑥−𝑎_𝑗(𝑥)‖ if 𝑥∈𝐴_𝑗 and then add an extra descent check and a boundary check to the iteration. As pointed out by Theorem 1 in [6], however, the sequence generated by the proposed Weiszfeld-type method is possible to be convergent to a nonoptimal point which is on the boundary of the regional customer.

In this section, a variational inequality approach is proposed to solve the general CSLP (11), where the locational constraints are imposed to the facility and the gauge is used as distance measuring function. Note that the study of variational inequality has received much attention due to its various applications arising in engineering, operations research, economics, transportation, and so forth; see, for example, [39–45]. Specifically, the CSLP (11) considered in this paper is first reformulated into an equivalent linear variational inequality (LVI), and then, a projection-contraction (PC) method is adopted to solve the LVI. Consequently, a sequence will be generated by the variational inequality approach, which is shown to be convergent to the optimal location of the facility 𝑥^𝑘_𝑖 of (11) even in the singular case. In addition, the closest points to the facility and the dual vectors with respect to the gauge can also be obtained from the generated sequence.

For convenience and succinctness, with the assumption that Λ_𝑖^𝑘−1 contains 𝑑 customers, throughout this section, we ignore some superscripts and subscripts in (11) and consider the simplified model of (11) without confusion: MCSLP:𝑥=argmin_𝑥∈Π𝐶(𝑥)∶=^𝑑_𝑗=1𝑠_𝑗𝛾𝑥−𝑎_𝑗(𝑥). (12) According to Proposition 1, it follows that CSLP (11), or equivalently MCSLP (12), is convex problem of 𝑥.

3.2.1. LVI Reformulation of MCSLP

For the gauge 𝛾(⋅) in (12) which is defined by (1), there exists a dual gauge, 𝛾^𝑑(⋅), defined by 𝛾^𝑑(𝑧)=max𝑧^𝑇𝑥∣𝛾(𝑥)≤1. (13) Let 𝐵^𝑑 be the unit ball of the dual gauge 𝛾^𝑑(⋅), which is also convex and compact and exactly the polar set of 𝐵. The dual gauge of 𝛾^𝑑(⋅) is again 𝛾(⋅), that is, 𝛾(𝑥)=max𝑧^𝑇𝑥∣𝛾^𝑑(𝑧)≤1, (14) which can also be rewritten as 𝛾(𝑥)=max_{𝑧∈𝐵𝑑}𝑧^𝑇𝑥. (15) 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: min_𝑥∈Πmax_{𝑧𝑗∈𝐵𝑑𝑠𝑗}^𝑑_𝑗=1𝑧^𝑇_𝑗𝑥−𝑎_𝑗(𝑥), (16) where each 𝑧_𝑗 is a vector in 𝐵^𝑑_𝑠𝑗={𝜉∈𝑅^𝑝∣𝛾^𝑑(𝜉)≤𝑠_𝑗}. Since 𝑎_𝑗(𝑥) is the closest point to 𝑥 in 𝐴_𝑗, we can introduce 𝑦_𝑗∈𝐴_𝑗 to replace 𝑎_𝑗(𝑥). Hence, (16) is equivalent to min_{𝑥∈Π,𝑦𝑗∈𝐴𝑗}max_{𝑧𝑗∈𝐵𝑑𝑠𝑗}^𝑑_𝑗=1𝑧^𝑇_𝑗𝑥−𝑦_𝑗. (17) Denote 𝑦=𝑦^𝑇₁,…,𝑦^𝑇_𝑑^𝑇,𝑧=𝑧^𝑇₁,…,𝑧^𝑇_𝑑^𝑇,𝐴=𝐴₁×⋯×𝐴_𝑑,ℬ^𝑑=𝐵^𝑑_𝑠1×⋯×𝐵^𝑑_𝑠𝑑, (18) and let (𝑥^∗,𝑦^∗,𝑧^∗)∈Π×𝐴×ℬ^𝑑 be the solution of (17); then, it follows that (𝑥^∗,𝑦^∗,𝑧^∗) is the saddle point of the objective function ∑^𝑑_𝑗=1𝑧^𝑇_𝑗(𝑥−𝑦_𝑗); that is, ^𝑑_𝑗=1𝑧^𝑇_𝑗𝑥^∗−𝑦^∗_𝑗≤^𝑑_𝑗=1𝑧_𝑗^∗𝑇𝑥^∗−𝑦^∗_𝑗≤^𝑑_𝑗=1𝑧_𝑗^∗𝑇𝑥−𝑦_𝑗,∀𝑥,𝑦,𝑧∈Π×𝐴×ℬ^𝑑. (19) Thus, (𝑥^∗,𝑦^∗,𝑧^∗) is the solution of the following linear variational inequality: 𝑥^∗∈Π,𝑦^∗∈𝐴,𝑧^∗∈ℬ^𝑑,𝑥−𝑥^∗𝑇𝑑_𝑗=1𝑧^∗_𝑗𝑦≥0,∀𝑥∈Π,_𝑗−𝑦^∗_𝑗^𝑇−𝑧^∗_𝑗≥0,∀𝑦_𝑗∈𝐴_𝑗,𝑧_𝑗−𝑧^∗_𝑗^𝑇−𝑥^∗−𝑦^∗_𝑗≥0,∀𝑧_𝑗∈𝐵^𝑑_𝑠𝑗. (20) A compact form of (20) is LVIΩ,𝑀,𝑞:𝑢^∗∈Ω,𝑢−𝑢^∗𝑇𝑀𝑢^∗+𝑞≥0,∀𝑢∈Ω, (21) where 𝑢=(𝑥^𝑇,𝑦^𝑇,𝑧^𝑇)^𝑇, Ω=Π×𝐴×ℬ^𝑑, 𝑀=0𝑁−𝑁^𝑇0,𝐼𝑁=²⋯⋯𝐼₂−𝐼₂0⋯00⋱⋱⋮⋮⋱⋱00⋯0−𝐼₂,𝑞=0. (22) Note that 𝑀 in (22) is a skew-symmetric matrix, then it is positive semidefinite, and thus, the linear variational inequality (21)-(22) is monotone.

Based on the deduction above, we know that if (𝑥^∗,𝑦^∗,𝑧^∗) is the solution of (17), that is, 𝑥^∗ is the solution of the MCSLP (12), then (𝑥^∗,𝑦^∗,𝑧^∗) will be the solution of the LVI (21)-(22). Further, we can prove that the MCSLP (12) and the LVI (21)-(22) are equivalent in the following theorem.

Theorem 3. The MCSLP (12) and the LVI (21)-(22) are equivalent in the sense that they have the same solution of 𝑥∈Π.

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 (𝑥^∗,𝑦^∗,𝑧^∗) is a solution of (17) if and only if (𝑥^∗,𝑦^∗,𝑧^∗) is a solution of LVI (21)-(22).
Let (𝑥^∗,𝑦^∗,𝑧^∗) be the solution of (17); then, according to the deduction above, we know that (𝑥^∗,𝑦^∗,𝑧^∗) is the solution of LVI (21)-(22).
On the other hand, let (𝑥^∗,𝑦^∗,𝑧^∗) be the solution of LVI (21)-(22) and 𝜙(𝑥,𝑦,𝑧)=∑^𝑑_𝑗=1𝑧^𝑇_𝑗(𝑥−𝑦_𝑗); then, the inequality (19) is true, which means that (𝑥^∗,𝑦^∗,𝑧^∗) is the saddle point of 𝜙(𝑥,𝑦,𝑧).
Note that (𝑥^∗,𝑦^∗,𝑧^∗) is the saddle point of 𝜙(𝑥,𝑦,𝑧) if and only if (𝑥^∗,𝑦^∗,𝑧^∗)∈Ω and max_{𝑧∈ℬ𝑑}𝜙𝑥^∗,𝑦^∗,𝑧=𝜙𝑥^∗,𝑦^∗,𝑧^∗=min_{𝑥∈Π,𝑦∈𝑌}𝜙𝑥,𝑦,𝑧^∗, (23) which implies that min_{𝑥∈Π,𝑦∈𝑌}max_{𝑧∈ℬ𝑑}𝜙𝑥,𝑦,𝑧≤max_{𝑧∈ℬ𝑑}𝜙𝑥^∗,𝑦^∗𝑥,𝑧=𝜙^∗,𝑦^∗,𝑧^∗=min_{𝑥∈Π,𝑦∈𝑌}𝜙𝑥,𝑦,𝑧^∗≤max_{𝑧∈ℬ𝑑}min_{𝑥∈Π,𝑦∈𝑌}𝜙𝑥,𝑦,𝑧. (24) On the other hand, let 𝑧^′ be any vector in ℬ^𝑑; then, we have min_{𝑥∈Π,𝑦∈𝑌}𝜙𝑥,𝑦,𝑧^′≤min_{𝑥∈Π,𝑦∈𝑌}max_{𝑧∈ℬ𝑑}𝜙𝑥,𝑦,𝑧. (25) We choose 𝑧^′ in (25) as the maximum point of the left term over 𝑧^′∈ℬ^𝑑; then, max_{𝑧∈ℬ𝑑}min_{𝑥∈Π,𝑦∈𝑌}𝜙𝑥,𝑦,𝑧≤min_{𝑥∈Π,𝑦∈𝑌}max_{𝑧∈ℬ𝑑}𝜙𝑥,𝑦,𝑧. (26) Combining (24) and (26), it follows that all terms in (24) are equal, and therefore, 𝜙𝑥^∗,𝑦^∗,𝑧^∗=min_{𝑥∈Π,𝑦∈𝑌}max_{𝑧∈ℬ𝑑}𝜙𝑥,𝑦,𝑧, (27) which implies that (𝑥^∗,𝑦^∗,𝑧^∗) 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 𝑧_𝑗 and the set 𝐵^𝑑_𝑠𝑗(𝑗=1,…,𝑑) in (16) are, respectively, the dual vector and dual ball in the space 𝑅^𝑝 which satisfy 𝑧_𝑗∈𝐵^𝑑_𝑠𝑗.

The norms especially 𝑙₁, 𝑙₂, and 𝑙_∞ are frequently used to measure distances in the literature; see, for example, [18, 19]. It should be noted that the gauge used in this paper is an extension of norms which include 𝑙₁, 𝑙₂, and 𝑙_∞. When the gauge 𝛾(⋅) is chosen as the 𝑙₁, 𝑙₂, and 𝑙_∞-norm, the dual gauge 𝛾^𝑑(⋅) will be the 𝑙_∞, 𝑙₂, and 𝑙₁-norm, respectively. Let 𝐵_𝑠𝑗,2=𝜉∈𝑅^𝑝∣𝜉₂≤𝑠_𝑗,ℬ₂=𝐵_𝑠1,2×⋯×𝐵_𝑠𝑑,2, (28) where ‖⋅‖₂ is the Euclidean norm; then, as a particular case of our single-source location problem (11), the problem under 𝑙₂-norm analyzed in [6] can be reformulated into the LVI (21)-(22) in which ℬ^𝑑 is equal to ℬ₂ and Π=𝑅^𝑝.

Further let 𝐵_𝑠𝑗,1=𝜉∈𝑅^𝑝∣𝜉₁≤𝑠_𝑗,𝐵_𝑠𝑗,∞=𝜉∈𝑅^𝑝∣𝜉_∞≤𝑠_𝑗, (29) where ‖⋅‖₁ and ‖⋅‖_∞ are the 𝑙₁,𝑙_∞-norm, respectively, and ℬ₁=𝐵_𝑠1,1×⋯×𝐵_𝑠𝑑,1,ℬ_∞=𝐵_𝑠1,∞×⋯×𝐵_𝑠𝑑,∞. (30) Then, the CSLP (11) under 𝑙₁-norm (the minisum model discussed in [16]) and the CSLP under 𝑙_∞-norm are equivalent to the LVI (21)-(22), where ℬ^𝑑 is, respectively, equal to ℬ_∞ and ℬ₁ and the locational constraints Π are both 𝑅^𝑝.

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].

Let 𝑊 be a nonempty closed convex set of 𝑅^𝒬. For a given 𝑣∈𝑅^𝒬, the projection of 𝑣 onto 𝑊 denoted by 𝑃_𝑊(𝑣) is the unique solution of the following problem: 𝑃_𝑊(𝑣)=argmin‖𝑢−𝑣‖₂∣𝑢∈𝑊. (31) A basic proposition of the projection mapping on a closed convex set is 𝑣−𝑃_𝑊(𝑣)^𝑇𝑢−𝑃_𝑊(𝑣)≤0,∀𝑣∈𝑅^𝒬,∀𝑢∈𝑊. (32)

It is well known (see, e.g., [50]) that for any 𝛽>0, 𝑢^∗ is the solution of LVI(Ω,𝑀,𝑞) if and only if 𝑒𝑢^∗,𝛽∶=𝑢^∗−𝑃_Ω𝑢^∗−𝛽𝑀𝑢^∗+𝑞=0. (33) In the literature of variational inequalities, 𝑒(𝑢,𝛽) is usually called the error bound of LVI, and it quantitatively measures how much 𝑢 fails to be the solution of LVI(Ω,𝑀,𝑞). Therefore, 𝑒(𝑢,𝛽) can serve as the stopping criterion for solving LVI(Ω,𝑀,𝑞) iteratively.

Let 𝑒(𝑢)=𝑒(𝑢,1),𝑔(𝑢)=𝑀^𝑇𝑒(𝑢)+,𝜑(𝑢)(𝑢)𝑀𝑢+𝑞=𝑒^𝑇𝑀𝑢+𝑞, (34) and Ω^∗ be the set of solutions of LVI(Ω,𝑀,𝑞); then, for the positive semidefinite (not necessarily symmetric) matrix 𝑀, the following theorem can be obtained.

Theorem 5 (Lemma 1 and Theorem 2 in [48]). Let 𝑢∈Ω, 𝑢^∗∈Ω^∗, 𝑔(𝑢), and 𝜑(𝑢) be defined as (34). Then, it holds that 𝑢−𝑢^∗𝑇𝑔(𝑢)≥𝜑(𝑢)≥‖𝑒(𝑢)‖²₂. (35)

For 𝑢∈Ω⧵Ω^∗, it follows from Theorem 5 that −𝑔(𝑢) is a descent direction of the unknown function ‖𝑢−𝑢^∗‖²₂. 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, 𝛼₁, 𝛼₂ (𝛼₁>𝛼₂) and 𝑢⁰∈Ω. Set 𝑘=0.
Step 1. Calculate 𝑒(𝑢^𝑘). If ‖𝑒(𝑢^𝑘)‖<𝜀, stop.
Step 2. Calculate 𝑔(𝑢^𝑘) and set 𝛼(𝑢^𝑘) as 𝛼𝑢^𝑘=𝑒𝑢^𝑘2𝑒𝑢^𝑘+𝑀^𝑇𝑒𝑢^𝑘2. (36)
Step 3. Calculate 𝑢^𝑘+1 as 𝑢^𝑘+1=𝑃_Ω𝑢^𝑘−𝜌𝛼𝑢^𝑘𝑔𝑤^𝑘,𝜌∈(0,2). (37)
Step 4. Adjust 𝛽 as follows 3𝛽=2𝛽𝛼𝑢^𝑘≥𝛼₁,23𝛽𝛼𝑢^𝑘≤𝛼₂,𝛽otherwise, (38) and set 𝑀=𝛽𝑀,𝑞=𝛽𝑞. (39) Let 𝑘=𝑘+1 and go to Step 1.

Remark 7. In Step 4 of Algorithm 6, the parameter 𝛽 is self-adaptive during the iterations according to the value of 𝛼(𝑢^𝑘). Note that 𝑀 is skew-symmetric, and thus, 𝛼(𝑢^𝑘) can also be rewritten as 𝛼𝑢^𝑘=𝑒𝑢^𝑘2𝑒𝑢^𝑘2+𝑀^𝑇𝑒𝑢^𝑘2, (40) which is shown to be in [0,1]. It follows from (39)-(40) that the two terms ‖𝑒(𝑢^𝑘)‖² and ‖𝑀^𝑇𝑒(𝑢^𝑘)‖² in the denominator of 𝛼(𝑢^𝑘) are balanced by the self-adaptive parameter 𝛽.

Theorem 8 (Theorem 3 in [48]). Let 𝑢^∗ be a solution of LVI (21)-(22); then, the sequence {𝑢^𝑘} generated by Algorithm 6 satisfies 𝑢^𝑘+1−𝑢^∗2≤𝑢^𝑘−𝑢^∗2−𝜌2−𝜌𝐼+𝑀^𝑇2𝑒𝑢^𝑘2. (41)

As a result, {𝑒(𝑢^𝑘)} converges to zero, and thus, all accumulation points of {𝑢^𝑘} are the solutions of LVI (21)-(22). However, it follows from (41) that ‖𝑢^𝑘+1−𝑢^∗‖ ≤ ‖𝑢^𝑘−𝑢^∗‖, which implies that {𝑢^𝑘} has only one accumulation point. Thus, the sequence {𝑢^𝑘} 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 solving 𝑚 CSLPs.

Mark that the CMLP (9) differs from MWP mainly in that the customers are represented by regions rather than points. Consequently, the CSLPs involved in the location phase are constrained location problems with regional demand and closest distances under gauge. No doubt that the numerical implementation of the heuristic algorithm to be presented is expected to be more complicated than the location-allocation algorithms for MWP. Therefore, how to accelerate the convergence of the proposed heuristic deserves further consideration. To achieve this objective, we here consider a particular strategy for the initial partition of regional customers or the initial locations of facilities. In practical implementation, we suggest to choose the solution of the following constrained multi-source Weber problem (CMWP) as the initial locations of facilities for CMLP: CMWP:min_{(𝑥1,…,𝑥𝑚)∈Π1×⋯×Π𝑚}𝐶^′(𝑋)=^𝑛_𝑗=1𝑠_𝑗min_{1≤𝑖≤𝑚}𝛾𝑥_𝑖−𝑔_𝑗, (42) where 𝑔_𝑗’s are geometric centers of the regional customers. Then, we apply the NCR to determine an initial partition of regional customers according to the solution of (42). For solving the constrained multi-source Weber problem (42), we employ the location-allocation heuristic algorithm in [35]. As we will show by numerical experiments, this initialization strategy can accelerate the convergence of the proposed algorithm greatly.

In the spirit of Cooper’s work, the new heuristic algorithm is ready to be presented for solving the targeted CMLP (9), and its iterative framework can be elaborated as follows.

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 (𝑥⁰₁,𝑥⁰₂,…,𝑥⁰_𝑚). Then, the initial partition of regional customers, which is denoted by Λ⁰={Λ⁰₁,Λ⁰₂,…,Λ⁰_𝑚}, is generated by the spirit of NCR heuristic (Step 3 in Algorithm 2). Set 𝑘=1.
Step 1 (location phase). Solve the involved CSLP (11) and find the location of facility 𝑥^𝑘_𝑖 for Λ_𝑖^𝑘−1 by the variational inequality approach. Denote 𝑋^𝑘=(𝑥₁^𝑘𝑇,…,𝑥_𝑚^𝑘𝑇)^𝑇.
Step 2 (allocation phase). Update the partition of regional customers Λ from Λ_𝑖^𝑘−1 to Λ^𝑘_𝑖 based on the spirit of NCR heuristic.
Step 3 If ‖𝑋^𝑘−𝑋^𝑘−1‖<𝜀, the current locations and partition are heuristic locations of facilities and heuristic partition of customers. Otherwise, set 𝑘=𝑘+1 and go to Step 1.

Remark 10. At the (𝑘+1)th iteration, it is recommended to use 𝑥^𝑘_𝑖 (and the corresponding 𝑦 and 𝑧) in Step 1 as the initial iterate in the variational inequality approach for solving 𝑥_𝑖^𝑘+1 (𝑖=1,…,𝑚), considering the fact that Λ^𝑘_𝑖 usually differs from Λ_𝑖^𝑘−1 slightly in practical implementation.

Remark 11. Compared to the main body of the proposed location-allocation heuristic algorithm, the workload of the initialization is relatively less. However, this initialization strategy can reduce its number of iterations and computing time, which will be verified by the numerical experiments to be reported in Section 6.2. Hence, the convergence of the proposed algorithm is accelerated greatly by this initialization strategy.

5. Convergence of the Proposed Heuristic Algorithm

In this section, we analyze the convergence of the proposed location-allocation heuristic (Algorithm 9). For simplification of our discussion, some notations are introduced as follows. Let 𝐴=𝐴₁×⋯×𝐴_𝑛, and recall Π=Π₁×⋯×Π_𝑚. For any 𝑋=(𝑥^𝑇₁,…,𝑥^𝑇_𝑚)^𝑇∈Π and 𝐶=(𝑐^𝑇₁,…,𝑐^𝑇_𝑛)^𝑇∈𝐴, we can define an ordered pair (𝑋,𝐶) and we can also define the function 𝜔(𝑋,𝐶), in the current partition of customers as follows: 𝜔(𝑋,𝐶)=^𝑚𝑛_𝑖=1𝑗=1𝑤_𝑖𝑗𝛾𝑥_𝑖−𝑐_𝑗=^𝑚_{𝑖=1𝑗∈𝒩∶𝐴𝑗∈Λ𝑖}𝑠_𝑗𝛾𝑥_𝑖−𝑐_𝑗, (43) which represents the objective functional value of CMLP (9) at (𝑋,𝐶).

During the implementation of Algorithm 9, we denote the map ℒ∶Π×𝐴→Π×𝐴 as the location operation in the Step 1 and the map 𝒜∶Π×𝐴→Π×𝐴 as the allocation operation in the Step 2. It follows that ℒ𝑋^𝑘,𝐶^𝑘=𝑋^𝑘+1,𝐶^𝑘,𝒜𝑋^𝑘+1,𝐶^𝑘=𝑋^𝑘+1,𝐶^𝑘+1, (44) where 𝑋^𝑘 and 𝑋^𝑘+1 are, respectively, the locations of facilities in the 𝑘th and (𝑘+1)th iteration, 𝐶^𝑘 is the variable of the closest points to 𝑋^𝑘 in the partition of Λ^𝑘, 𝐶^𝑘 is the closest points to the new 𝑋^𝑘+1 in the partition of Λ^𝑘, and 𝐶^𝑘+1 is the closest points to the new 𝑋^𝑘+1 in the new partition of Λ^𝑘+1. Then, the iterate scheme of the location-allocation heuristic is 𝑋^𝑘+1,𝐶^𝑘+1=𝒜ℒ𝑋^𝑘,𝐶^𝑘. (45)

Let 𝑆(𝑋⁰,𝐶⁰) denote the iterative sequence generated by the location-allocation heuristic for CMLP with the initial iterate (𝑋⁰,𝐶⁰). 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 𝑆(𝑋⁰,𝐶⁰).

Proposition 12. 𝑆(𝑋⁰,𝐶⁰) is strictly monotone in the sense that 𝜔(𝑋^𝑘+1,𝐶^𝑘+1)<𝜔(𝑋^𝑘,𝐶^𝑘) if 𝑋^𝑘+1≠𝑋^𝑘.

Proof. Since 𝑋^𝑘+1≠𝑋^𝑘, there exists at least one 𝑖∈ℳ such that 𝑥_𝑖^𝑘+1≠𝑥^𝑘_𝑖. For such 𝑖’s according to the following convex optimization problem 𝑥_𝑖^𝑘+1=argmin_{𝑥∈Π𝑖}𝐶_𝑖^𝑘+1(𝑥)=_{𝑗∈𝒩∶𝐴𝑗∈Λ𝑘𝑖}𝑠_𝑗𝛾𝑥−𝑎_𝑗(𝑥), (46) we have 𝜔𝑋^𝑘+1,𝐶^𝑘≤𝜔𝑋^𝑘,𝐶^𝑘. (47) Based on Remark 10, we know that if 𝑥^𝑘_𝑖 is the solution of (46), then 𝑥_𝑖^𝑘+1 will be equal to 𝑥^𝑘_𝑖. Therefore, 𝑥_𝑖^𝑘+1≠𝑥^𝑘_𝑖 implies that 𝑥^𝑘_𝑖 is not the solution of (46), and thus, 𝐶_𝑖^𝑘+1(𝑥^𝑘+1)<𝐶^𝑘_𝑖(𝑥^𝑘). It follows that 𝜔𝑋^𝑘+1,𝐶^𝑘<𝜔𝑋^𝑘,𝐶^𝑘. (48) On the other hand, based on the principle of NCR in the allocation phase of Algorithm 9, we also have 𝜔𝑋^𝑘+1,𝐶^𝑘+1≤𝜔𝑋^𝑘+1,𝐶^𝑘. (49) By Combining (48) and (49), it follows that 𝜔𝑋^𝑘+1,𝐶^𝑘+1<𝜔𝑋^𝑘,𝐶^𝑘. (50) The proof is complete.

Based on the monotonicity of the generated sequence 𝑆(𝑋⁰,𝐶⁰), the following theorem can be proved.

Theorem 13. Let {𝑆_𝑘}∶=𝑆(𝑋⁰,𝐶⁰). Then, the generated sequence {𝑆_𝑘} satisfies that(1)𝜔(𝑆_𝑘)→𝜔(𝑆) for some 𝑆∈Π×𝐴,(2)all accumulation points of {𝑆_𝑘} have the same objective functional values.

Proof. After a finite number of iterations, if 𝑋^𝐽+1=𝑋^𝐽, then the iterates after 𝑆_𝐽 will be constant, and thus, {𝑆_𝑘} is convergent to (𝑋^𝐽,𝐶^𝐽)∈Π×𝐴, and the two assertions are both true.
Below, we will discuss the case that 𝑋^𝑘+1≠𝑋^𝑘 for any 𝑘∈𝑁. First, we prove the first assertion. Since 𝑆_𝑘∈Π×𝐴 and Π×𝐴 is a compact space, it follows from the Bolzano-Weierstrass theorem that there exists a subsequence of {𝑆_𝑘} which, say {𝑆_𝑘}_𝐾, converges to an element 𝑆∈Π×𝐴; that is, lim_𝐾→∞𝑆_𝑘𝐾⟶𝑆,𝑆∈Π×𝐴. (51) Note that 𝜔(𝑋,𝐶) is a continuous function according to (43); then, lim_𝐾→∞𝜔𝑆_𝑘𝐾⟶𝜔(𝑆),𝑆∈Π×𝐴. (52) Due to that {𝜔(𝑆_𝑘)} is a monotone sequence (Proposition 12) and has lower bound, then {𝜔(𝑆_𝑘)} is convergent. Thus, any subsequence of {𝜔(𝑆_𝑘)} will be convergent to the same value. Note that {𝜔({𝑆_𝑘}_𝐾)} is a subsequence of {𝜔(𝑆_𝑘)} and it is convergent to 𝜔(𝑆); then, it follows that lim_𝑘→∞𝜔𝑆_𝑘⟶𝜔(𝑆),𝑆∈Π×𝐴. (53)
The second assertion can easily be proved. Let 𝑃 be an accumulation point of {𝑆_𝑘}; then, there exists one subsequence {𝑆_𝑘}_𝐾′ which converges to 𝑃, and due to the continuity of 𝜔(𝑋,𝐶), we have 𝜔𝑆_𝑘𝐾′⟶𝜔(𝑃). (54) The first assertion has shown that {𝜔(𝑆_𝑘)} is convergent, and note that {𝜔({𝑆_𝑘}_𝐾′)} is a subsequence of {𝜔(𝑆_𝑘)}; then, it follows that 𝜔(𝑃)=lim_{𝐾′→∞}𝜔𝑆_𝑘𝐾′=lim_𝑘→∞𝜔𝑆_𝑘. (55) Thus, all accumulation points of {𝑆_𝑘} have the same objective functional values equal to lim_𝑘→∞𝜔(𝑆_𝑘).

Lemma 14. ℒ∶Π×𝐴→Π×𝐴 which is defined in (44) is a closed map over Π×𝐴.

Proof. Note that the CSLP (11) is a convex problem, then argmin_{𝑥∈Π𝑖𝑗∈𝒩∶𝐴𝑗∈Λ𝑖𝑘−1}𝑠_𝑗min_{𝑞𝑗∈𝐴𝑗}𝛾𝑥_𝑖−𝑞_𝑗,𝑖=1,2,…,𝑚, (56) are continuous. Since Π×𝐴 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 Π×𝐴.

Lemma 15. Let 𝜐₀ be a given vector in Π×𝐴 and Δ∶={𝜐∈Π×𝐴∣𝜔(𝜐)≤𝜔(𝜐₀)}. 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 {𝜐_𝑘} with 𝜐_𝑘∈Δ, since Π×𝐴 is compact, according to Bolzano-Weierstrass, there exists a convergent subsequence {𝜐_𝑘}_𝐾 of {𝜐_𝑘} such that lim_𝐾→∞𝜐_𝑘𝐾⟶𝜐. (57) Due to the continuity of 𝜔, it follows that lim_𝐾→∞𝜔𝜐_𝑘𝐾=𝜔(𝜐). (58)
On the other hand, {𝜐_𝑘}_𝐾∈Δ implies 𝜔𝜐_𝑘𝐾≤𝜔𝜐₀. (59) By combining (58), (59), and the continuity of 𝜔, the following inequality is obtained: 𝜔(𝜐)≤𝜔𝜐₀, (60) 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 (Algorithm 9). Let Ξ⊆Π×𝐴 be the nonempty local solution set of CMLP (9). Recall that in the location-allocation Cooper algorithm for MWP, if 𝑋^𝐽+1=𝑋^𝐽 occurs after a finite number of iterations, the iterates after 𝑋^𝐽 will be constant. Then, no further improvement is possible for MWP, and it follows from [32, 34] that the 𝑋^𝐽 is a local solution of MWP. Similarly, in the proposed location-allocation heuristic algorithm for CMLP, if 𝑋^𝐽+1=𝑋^𝐽 occurs, the iterates after 𝑆_𝐽=(𝑋^𝐽,𝐶^𝐽) 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, 𝑆_𝐽, is obtained. Hence in this case, {𝑆_𝑘} is convergent to the 𝑆_𝐽∈Ξ. However, it is not assured that 𝑋^𝐽+1=𝑋^𝐽 always occurs during the implementation of Algorithm 9, and therefore, we assume that 𝑋^𝑘+1≠𝑋^𝑘 for any 𝑘∈𝑁 and prove the convergence in this case.

Theorem 16. Assume that 𝑋^𝑘+1≠𝑋^𝑘 for any 𝑘∈𝑁; then, all the accumulation points of the sequence {𝑆_𝑘} belong to Ξ.

Proof. Let 𝑆 be an accumulation point of {𝑆_𝑘}. Due to 𝑆_𝑘∈Π×𝐴 and the compactness of Π×𝐴, we know that (1)  𝑆∈Π×𝐴 and (2) there exists a subsequence {𝑆_𝑘}_𝐾 which is convergent to 𝑆. According to Theorem 13, we know that lim_𝑘→∞𝜔𝑆_𝑘=𝜔(𝑆),𝑆∈Π×𝐴. (61) So, it is enough to prove 𝑆∈Ξ. We prove this by contradiction. Assume that 𝑆∉Ξ; that is, 𝑆 is not a solution, and we consider the subsequence {𝑆_𝑘+1}_𝐾. Denote Δ={𝜐∈Π×𝐴∣𝜔(𝜐)≤𝜔(𝑆₀)}. Due to Proposition 12, it follows that for all 𝑘∈𝑁 we have (𝑋^𝑘+1,𝐶^𝑘)∈Δ and (𝑋^𝑘+1,𝐶^𝑘+1)∈Δ. According to Lemma 15, Δ is compact, and thus, there exists 𝐾^′⊂𝐾 such that lim_{𝐾′→∞}𝑋^𝑘+1,𝐶^𝑘_𝐾′=𝜐₁,lim_{𝐾′→∞}𝑋^𝑘+1,𝐶^𝑘+1_𝐾′=𝜐₂. (62) According to Lemma 14, the map ℒ is closed at 𝑆∈Π×𝐴; then, it follows that 𝜐₁=ℒ(𝑆). Further, due to 𝑆∉Ξ, 𝑣₁ will be not equal to 𝑆. Otherwise, we can choose 𝑆 as the initial iterate of the location-allocation heuristic algorithm, then, the sequence generated by the algorithm will be constant. It follows from the first case (i.e., 𝑋^𝐽+1=𝑋^𝐽) that 𝑆 will be a local solution, which contradicts with 𝑆∉Ξ. Therefore, we can obtain 𝑣₁≠𝑆. Together with the monotone proposition of ℒ (48), we have the inequality 𝜔𝜐₁<𝜔(𝑆). (63)
On the other hand, note that 𝒜(𝑋^𝑘+1,𝐶^𝑘)=(𝑋^𝑘+1,𝐶^𝑘+1); then, by the monotonicity of 𝒜 (49), it follows that 𝜔𝑋^𝑘+1,𝐶^𝑘+1_𝐾′≤𝜔𝑋^𝑘+1,𝐶^𝑘_𝐾′, (64) and thus, by taking the limit for (64) and by the continuity of 𝜔, we know 𝜔(𝜐₂)≤𝜔(𝜐₁). Combining this with (63), we obtain 𝜔𝜐₂<𝜔(𝑆). (65) However, note that 𝜐₂ and 𝑆 are two accumulation points of {𝑆_𝑘}, and according to the second assertion of Theorem 13, 𝜔(𝜐₂)=𝜔(𝑆), which will contradict with (65). Therefore, our assumption is wrong, and thus, 𝑆∈Ξ.

As a result, we have the following convergence theorem for the sequence generated by the proposed location-allocation algorithm.

Theorem 17. The sequence 𝑆(𝑋⁰,𝐶⁰) generated by the proposed location-allocation heuristic algorithm either converges to a point in Ξ or all accumulation points of 𝑆(𝑋⁰,𝐶⁰) belong to Ξ.

6. Numerical Results

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

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 𝑢⁰ in Algorithm 6 can be chosen arbitrarily in Ω. In practical implementation, however, we choose 𝑢⁰ judiciously similar to the initialization strategy in the location-allocation heuristic: let 𝑔_𝑗(𝑗=1,…,𝑑) be the centers of regional customers, solve the following single-source Weber problem (SWP): 𝑥^∗=argmin_𝑥∈Π^𝑑_𝑗=1𝛾𝑥−𝑔_𝑗 (66) 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 𝛼₁ and 𝛼₂ 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, 𝑑=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 𝑠_𝑗=1, 𝑗=1,…,5.

In order to clarify the comparison of two methods, we choose the same stopping criterion as ‖𝑥^𝑘+1−𝑥^𝑘‖≤10⁻⁴ (throughout this section, ‖⋅‖ is the 𝑙_∞-norm). We test this example for 100 times with the same initial iterate for the two methods which is randomly generated in [0,5]×[0,3], and the numerical results including the location of new facility, the closest points to the facility, number of iterations, and computing time in units of second are reported in Table 1.

Table 1 

Numerical results for Example 1 given in [16].

According to Table 1, it follows that both methods can get the optimal location of the new facility. Though our variational inequality approach needs smaller iteration numbers and less computing time, both methods are efficient for solving Example 1 given in [6].

Recall that the sequence generated by the Weiszfeld-type method in [16] is possible to be convergent to a nonoptimal point on the boundary of the regional customer, as indicated in [16]. The variational inequality approach, however, can obtain the optimal location of the new facility, which is guaranteed by the theoretical analysis in Section 3.2. To illustrate this, we compare two methods by solving the following particular example.

Example 19. Similar to Example 18, the number of customers 𝑑=5, and all customers are unit squares whose edges are parallel to the axes. Let 𝑑₁ be the distance between any two neighboring customers, and we set 𝑑₁ equal to 1, 0.1, 0.01, and 0.001, respectively. The weights 𝑠_𝑗(𝑗=1,…,5) are randomly generated in the area of regional customers.

Note that in Example 19 the parameter 𝑑₁ reflects how close the customers are away from its neighborhoods. We test this example for a large number of times with the stopping criterion ‖𝑥^𝑘+1−𝑥^𝑘‖≤10⁻⁴, 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.₀” 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 Example 19.

According to Table 2, it follows that both the VI approach and the Weiszfeld-type method are efficient for solving this particular example, and both of them need a small number of iterations and little computing time. In comparison of two methods, we can find that VI approach needs more iteration times and computing time than Weiszfeld-type method. From the column of “Impro. Percent,” however, we can conclude that VI approach can obtain a better solution (in fact, the solution obtained by VI approach is the optimal location of new facility) than Weiszfeld-type method. In addition, according to the last column, we find that when 𝑑₁ decreases, that is, the customers become closer and closer, the frequency that Weiszfeld-type method obtains the optimal solution gets smaller and smaller. When 𝑑₁ 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 𝑑₁=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 𝑑₁ is 1, which means that the customers are enough far away from one another, Weiszfeld-type method can obtain the optimal location of facility in most cases, and the frequency even exceeds 90.

Since our main effort in this paper is to solve the general location problem under gauge and locational constraint, it is necessary to apply the variational inequality approach to solve some CSLPs. In particular, we test a large number of randomly generated CSLPs with the number of customers 𝑑 from 10 to 2000. In the experiments, all regional customers are assumed to be square units, and their edges are parallel to the coordinate axes. The geometric centers of all regional customers are randomly generated in [−100,100]²; the weights of the regional customer are all randomly chosen in (1,5); the locational constraint is ‖𝑥−𝑂‖≤𝑟, where the center 𝑂 is randomly generated in [−100,100]², and the radius 𝑟 is randomly generated in (1,5); the stopping criterion of VI approach is chosen as 𝑥^𝑘+1−𝑥^𝑘≤10⁻⁴,𝑒𝑢^𝑘+1≤10⁻⁴, (67) and the initial iterate is randomly generated in [−100,100]²; the gauge 𝛾(⋅) is generated with the unit ball set as 92𝑥+3²+12𝑦²=16. (68) For each 𝑑, we test one hundred randomly generated CSLPs, and the average numerical results are reported in Table 3. To illustrate the effect of the initialization strategy of VI approach, we also report the results of VI approach without the initialization strategy. The columns of “VI approach with Initial.” and “VI approach without Initial.,” respectively, report the average number of iterations and average computing time of variational inequality approach with and without the initialization strategy.

Table 3 

Numerical results of VI approach for CSLP (12).

According to Table 3, it is easy to conclude that the variational inequality approach is effective for solving CSLP under gauge considering the difficulty of this problem. In addition, the number of iterations of “VI approach with Initial.” is less than that of “VI approach without Initial.,” which shows that the initialization strategy can accelerate the convergence of the variational inequality approach. This strategy, however, does not necessarily reduce the computing time of VI approach, especially when 𝑑 is small, for example, 𝑑=10,20,50,and100, which can be explained as follows. When the number of customers 𝑑 is small, the variational inequality problem is small scale, and thus, it can be solved in a short time. In this case, the computational workload of initialization plays an important role in the total workload, and therefore, the computing time of “VI approach with Initial.” is greater than that of “VI approach without Initial.” due to the computational iterations for initialization. With 𝑑 increasing, the scale of VI problem as well as the number of iterations becomes larger. Then, in the comparison of the workload of VI approach, the workload of initialization can almost be ignored, and therefore, the computing time of “VI approach with Initial.” will be smaller than that of “VI approach without Initial.” As a matter of fact, Table 3 reveals the computational necessity of the initialization strategy of variational inequality approach for large-scale CSLP; for example, the iteration number and computing time are reduced about 1/5 by the initialization strategy when 𝑑=2000.

6.2. Numerical Results of Heuristic Algorithm for CMLP

This subsection applies the proposed location-allocation heuristic algorithm (Algorithm 9) to solve a large number of CMLPs (9) and also verifies the necessity of the initialization strategy in Algorithm 9. In the experiments, we again generate a large number of CMLPs with unit square customers and assume that the edges of these regions are parallel to the coordinate axes. The geometric centers of all the demand regions are randomly generated in −100100², and all the demands, 𝑠_𝑗(𝑗=1,2,…,𝑛), are randomly generated in 110. The gauge 𝛾(⋅) is also defined with the unit ball set as (68). We test the scenario with 𝑛=100,200,500,and1000 and 𝑚=2,4,6,8,and10; the locational constraints are ‖𝑥−𝑂_𝑗‖≤𝑟_𝑗(𝑗=1,…,𝑚), where the radius 𝑟_𝑗 is randomly generated in [1,10], and the center 𝑂_𝑗 is given in advance as follows: 𝑚=2∶𝑂₁=(−50,0)^𝑇,𝑂₂=(50,0)^𝑇;𝑚=4∶𝑂₁=(−50,−50)^𝑇,𝑂₂=(50,−50)^𝑇,𝑂₃=(−50,50)^𝑇,𝑂₄=(50,50)^𝑇;𝑚=6∶𝑂₁=(−50,−50)^𝑇,𝑂₂=(0,−50)^𝑇,𝑂₃=(50,−50)^𝑇,𝑂₄=(−50,50)^𝑇,𝑂₅=(0,50)^𝑇,𝑂₆=(50,50)^𝑇;𝑚=8∶𝑂₁=(−50,−50)^𝑇,𝑂₂=(0,−50)^𝑇,𝑂₃=(50,−50)^𝑇,𝑂₄=(−50,50)^𝑇,𝑂₅=(0,50)^𝑇,𝑂₆=(50,50)^𝑇,𝑂₇=(−25,0)^𝑇,𝑂₈=(25,0)^𝑇;𝑚=10∶𝑂₁=(−50,−50)^𝑇,𝑂₂=(0,−50)^𝑇,𝑂₃=(50,−50)^𝑇,𝑂₄=(−50,50)^𝑇,𝑂₅=(0,50)^𝑇,𝑂₆=(50,50)^𝑇,𝑂₇=(−60,0)^𝑇,𝑂₈=(−20,0)^𝑇,𝑂₉=(20,0)^𝑇,𝑂₁₀=(60,0)^𝑇. (69) The initial locations of facilities are randomly generated in −100100², and the stopping criterion used in Algorithms 6 and 9 is chosen as 𝑥^𝑘+1−𝑥^𝑘<10⁻⁴. (70)

To show the significance of the initialization strategy, we compare the numerical performance of the location-allocation heuristic algorithm with initialization strategy (denoted by “Algorithm 9 with Initial.”) and without this initialization strategy (denoted by “Algorithm 9 without Initial.”). In the initialization step of Algorithm 9, the location-allocation algorithm in [35] is adopted to solve the corresponding CMWP (42), where a PC method is proposed to solve the subproblems in location phase, and the numbers of iterations of the algorithm in [35] (denoted by “Iter₀”) and the average iteration numbers of PC method in one iteration of the algorithm (denoted by “Iter₀-PC”) are reported. The columns of “Iter.” and “CPU,” respectively, report the number of iterations and computing time of Algorithm 9 with and without initialization strategy. Since the efficiency of Algorithm 9 is mainly determined by the number of iterations of the variational inequality approach, we also report the average number of iterations of Algorithm 6 in one iteration of Algorithm 9 (denoted by “Iter.-PC”). For each given pair (𝑛,𝑚), we test the CMLP for 100 times, and the computational performance is reported in Table 4.

Table 4 

Numerical results of location-allocation heuristic for CMLP.

It follows from Table 4 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 Table 4 deserves further illustration. Recall that with the number of new facilities (𝑚) increasing, the number of subproblems (CSLPs) in location phase increases too. According to Table 4, however, we find that for fixed number of customers (𝑛), with 𝑚 increasing, the number of iterations for solving 𝑚 CSLPs in one iteration of Algorithm 9 does not increase but almost decreases with 𝑚, especially for large-scale CMLP. This can be illustrated roughly as follows. For fixed 𝑛, when 𝑚 increases, the average number of customers in each Λ^𝑘_𝑖 becomes smaller, which implies that the scale of the involved CSLP (11) in location phase is smaller. According to Table 3, it follows that we need smaller number of iterations for small-scale CSLP, and thus, the total number of iterations for solving 𝑚 CSLPs decreases. Similarly, due to the same reason, the computing time for solving CMLP also decreases with 𝑚 increasing, as reported in the column of “CPU” in Table 4.

7. Conclusion

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

Acknowledgments

Jian-lin Jiang is supported by NSFC no. 11101211, JSNSF no. BK2011719, the Fundamental Research Funds for the Central Universities no. NZ2012306, and the project sponsored by SRF for ROCS, SEM.