- About this Journal ·
- Abstracting and Indexing ·
- Aims and Scope ·
- Annual Issues ·
- Article Processing Charges ·
- Articles in Press ·
- Author Guidelines ·
- Bibliographic Information ·
- Citations to this Journal ·
- Contact Information ·
- Editorial Board ·
- Editorial Workflow ·
- Free eTOC Alerts ·
- Publication Ethics ·
- Reviewers Acknowledgment ·
- Submit a Manuscript ·
- Subscription Information ·
- Table of Contents

Abstract and Applied Analysis

Volume 2013 (2013), Article ID 624398, 15 pages

http://dx.doi.org/10.1155/2013/624398

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

Department of Mathematics, College of Science, Nanjing University of Aeronautics and Astronautics, Nanjing 210016, China

Received 1 August 2013; Accepted 27 August 2013

Academic Editor: Abdellah Bnouhachem

Copyright © 2013 Jian-Lin Jiang et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### 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 is called the unit ball of due to This way to define a gauge from a compact convex set was first introduced by Minkowski [29]. The gauge satisfies the following properties: (1) and ; (2) for any ; (3) for any and .It follows from and that any gauge is a convex function of . The distance measuring function can be derived from a gauge by When is symmetric around the origin, according to the definition of (1), we have Combining (3) and (4), it follows that which means that the distance measured by is symmetric. On the contrary, when is not symmetric around the origin, (4) does not hold any more, and thus, the distance measured by the gauge is asymmetric. Thus, when different compact convex sets are used as unit balls, different gauges (symmetric or asymmetric) can be generated and employed to measure distances in location problems, which depends on the requirements of practical applications.

Let denote the set of regional customers, and let be the location of the th facility. Each regional customer is simply assumed to be closed and convex. We denote the closest point in to the facility by Then, the closest distance between the facility and the regional customer , denoted by , can be represented by

When the gauge is used to measure distances, we have the following proposition for and which is similar to that in [16].

Proposition 1. *Let be the location of a facility; then, the closest point in (6) is well defined, and the closest distance 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 ; then, due to and in and the convexity of , it follows that , and thus, we have
Therefore, is convex with respect to , and the proof is complete.

Based on the notations introduced above, now the constrained multi-source location problem (abbreviated as CMLP) we consider in this paper takes the following formulation: where 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 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 , , 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: 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 , , and then . At the th iteration, let be the disjoint partition of in the sense that and (for ), and each () 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 , : 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 ( and ), if () is the nearest point for among all computed by (11), then and . If solved by (11) is the nearest facility for each regional customer in for any , then () are the desirable locations of facilities and stop. Otherwise, we set 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* , *do.**Step **1.* Set ( stores the number of reassignments);*Step **2.* Compute the facility of by solving CSLP (11), for ;*Step **3.* For do: for ; if and , then , ; . *Step **4.* If , then the iteration terminates with being the desirable locations for facilities and the customers in being served by ().

##### 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 contains customers, throughout this section, we ignore some superscripts and subscripts in (11) and consider the simplified model of (11) without confusion: According to Proposition 1, it follows that CSLP (11), or equivalently MCSLP (12), is convex problem of .

###### 3.2.1. LVI Reformulation of MCSLP

For the gauge in (12) which is defined by (1), there exists a dual gauge, , defined by Let be the unit ball of the dual gauge , which is also convex and compact and exactly the polar set of . The dual gauge of is again , that is, which can also be rewritten as For more details about gauge and dual gauge, as well as their unit balls, the readers can be referred to [27].

According to (15), MCSLP (12) is equivalent to the following min-max problem: where each is a vector in . Since is the closest point to in , we can introduce to replace . Hence, (16) is equivalent to Denote and let be the solution of (17); then, it follows that is the saddle point of the objective function ; that is, Thus, is the solution of the following linear variational inequality: A compact form of (20) is where , , Note that in (22) is a skew-symmetric matrix, then it is positive semidefinite, and thus, the linear variational inequality (21)-(22) is monotone.

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

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

*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 ; 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
which implies that
On the other hand, let be any vector in ; then, we have
We choose in (25) as the maximum point of the left term over ; then,
Combining (24) and (26), it follows that all terms in (24) are equal, and therefore,
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 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 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 where and are the -norm, respectively, and 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: A basic proposition of the projection mapping on a closed convex set is

It is well known (see, e.g., [50]) that for any , is the solution of LVI() if and only if In the literature of variational inequalities, is usually called the error bound of LVI, and it quantitatively measures how much fails to be the solution of LVI(). Therefore, can serve as the stopping criterion for solving LVI() iteratively.

Let 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
*

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 , , , () and . Set .*Step **1.* Calculate . If , stop.*Step **2.* Calculate and set as
*Step **3.* Calculate as
*Step **4.* Adjust as follows
and set
Let 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
which is shown to be in . 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
*

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 ≤ , 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: 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 .*Step **1 (location phase).* Solve the involved CSLP (11) and find the location of facility for by the variational inequality approach. Denote .*Step **2 (allocation phase).* Update the partition of regional customers from to based on the spirit of NCR heuristic.*Step **3* If , the current locations and partition are heuristic locations of facilities and heuristic partition of customers. Otherwise, set and go to Step 1.

*Remark 10. *At the 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 (), considering the fact that usually differs from 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: 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 where and are, respectively, the locations of facilities in the th and th iteration, is the variable of the closest points to in the partition of , is the closest points to the new in the partition of , and is the closest points to the new in the new partition of . Then, the iterate scheme of the location-allocation heuristic is

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 if .*

*Proof. *Since , there exists at least one such that . For such ’s according to the following convex optimization problem
we have
Based on Remark 10, we know that if is the solution of (46), then will be equal to . Therefore, implies that is not the solution of (46), and thus, . It follows that
On the other hand, based on the principle of NCR in the allocation phase of Algorithm 9, we also have
By Combining (48) and (49), it follows that
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 , 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 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,
Note that is a continuous function according to (43); then,
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

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
The first assertion has shown that is convergent, and note that is a subsequence of ; then, it follows that
Thus, all accumulation points of have the same objective functional values equal to .

Lemma 14. * which is defined in (44) is a closed map over . *

*Proof. *Note that the CSLP (11) is a convex problem, then
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
Due to the continuity of , it follows that

On the other hand, implies
By combining (58), (59), and the continuity of , the following inequality is obtained:
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 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 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 always occurs during the implementation of Algorithm 9, and therefore, we assume that for any and prove the convergence in this case.

Theorem 16. *Assume that 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 and there exists a subsequence which is convergent to . According to Theorem 13, we know that
So, it is enough to prove . We prove this by contradiction. Assume that ; that is, is not a solution, and we consider the subsequence . Denote . Due to Proposition 12, it follows that for all we have and . According to Lemma 15, is compact, and thus, there exists such that
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., ) that will be a local solution, which contradicts with . Therefore, we can obtain . Together with the monotone proposition of (48), we have the inequality

On the other hand, note that ; then, by the monotonicity of (49), it follows that
and thus, by taking the limit for (64) and by the continuity of , we know . Combining this with (63), we obtain
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 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 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 be the centers of regional customers, solve the following single-source Weber problem (SWP): by the projection-contraction method in [35], and then use its solution as the initial iterate for Algorithm 6. We call this the initialization strategy of variational inequality approach. In addition, throughout our experiments of VI approach, the 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, ; 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 , , , , and , and , .

In order to clarify the comparison of two methods, we choose the same stopping criterion as (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 , and the numerical results including the location of new facility, the closest points to the facility, number of iterations, and computing time in units of second are reported in Table 1.

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

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

*Example 19. *Similar to Example 18, the number of customers , 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 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 , 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 ” gives the average iteration times of initialization in the variational inequality approach, and the column of “No. of Iter.” reports the average iteration times of both methods. The two columns of “CPU” give the average computing time in units of second for variational inequality approach (including the computing time for initialization) and Weiszfeld-type method, respectively. The columns of “Obj.” give the average objective functional value obtained by the two methods, and the column of “Impro. Percent” gives the improvement percentage in objective functional values of the VI approach to Weiszfeld-type method. Remark that the convergence of VI approach to the optimal location of new facility is guaranteed, and then the column of “Freq. Num.” reports the frequency among one hundred times that the Weiszfeld-type method can get the same solution as VI approach; that is, it does not converge to the nonoptimal solution on the boundary of the customer.

According to Table 2, it follows that both the VI approach and the Weiszfeld-type method are efficient for solving this particular example, and both of them need a small number of iterations and little computing time. In comparison of two methods, we can find that VI approach needs more iteration times and computing time than Weiszfeld-type method. From the column of “Impro. Percent,” however, we can conclude that VI approach can obtain a better solution (in fact, the solution obtained by VI approach is the optimal location of new facility) than Weiszfeld-type method. In addition, according to the last column, we find that when 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 , 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 ; the weights of the regional customer are all randomly chosen in ; the locational constraint is , where the center is randomly generated in , and the radius is randomly generated in ; the stopping criterion of VI approach is chosen as and the initial iterate is randomly generated in ; the gauge is generated with the unit ball set as For each , we test one hundred randomly generated CSLPs, and the average numerical results are reported in Table 3. To illustrate the effect of the initialization strategy of VI approach, we also report the results of VI approach without the initialization strategy. The columns of “VI approach with Initial.” and “VI approach without Initial.,” respectively, report the average number of iterations and average computing time of variational inequality approach with and without the initialization strategy.

According to Table 3, it is easy to conclude that the variational inequality approach is effective for solving CSLP under gauge considering the difficulty of this problem. In addition, the number of iterations of “VI approach with Initial.” is less than that of “VI approach without Initial.,” which shows that the initialization strategy can accelerate the convergence of the variational inequality approach. This strategy, however, does not necessarily reduce the computing time of VI approach, especially when is small, for example, , 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 by the initialization strategy when .

##### 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 , and all the demands, , are randomly generated in . The gauge is also defined with the unit ball set as (68). We test the scenario with and ; the locational constraints are , where the radius is randomly generated in , and the center is given in advance as follows: The initial locations of facilities are randomly generated in , and the stopping criterion used in Algorithms 6 and 9 is chosen as

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 “”) and the average iteration numbers of PC method in one iteration of the algorithm (denoted by “-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 times, and the computational performance is reported in Table 4.

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 .

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.

#### References

- A. Weber,
*UBer Den Standort Der Industrien, 1. Teil: Reine Theorie Des Standortes*, Mohr Siebeck, Tübingen, Germany, 1909. - E. Weiszfeld, “Sur le point pour lequel la somme des distances de n points donnes est minimum,”
*Tohoku Mathematical Journal*, vol. 43, pp. 355–386, 1937. View at Google Scholar - R. F. Love, J. G. Morris, and G. O. Wesolowsky,
*Facilities Location: Models and Methods,*, vol. 7, North-Holland Publishing, Amsterdam, Netherlands, 1988. View at MathSciNet - F. Plastria, “Continuous location problems: research, results and questions,” in
*Facility Location: A Survey of Applications and Methods*, Z. Drezner, Ed., pp. 225–260, Springer, New York, NY, USA, 1995. View at Google Scholar - C. D. Bennett and A. Mirakhor, “Optimal facility location with respect to several regions,”
*Journal of Regional Science*, vol. 14, no. 1, pp. 131–136, 1974. View at Google Scholar - J. Brimberg and G. O. Wesolowsky, “Minisum location with closest Euclidean distances,”
*Annals of Operations Research*, vol. 111, pp. 151–165, 2002. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - Z. Drezner and G. O. Weslowsky, “Optimal location of a facility relative to area demands,”
*Naval Research Logistics Quarterly*, vol. 27, no. 2, pp. 199–206, 1980. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - J. Jiang and X. Yuan, “A Barzilai-Borwein-based heuristic algorithm for locating multiple facilities with regional demand,”
*Computational Optimization and Applications*, vol. 51, no. 3, pp. 1275–1295, 2012. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - A. Suzuki and Z. Drezner, “The
*p*-center location problem in an area,”*Location Science*, vol. 4, no. 1-2, pp. 69–82, 1996. View at Google Scholar · View at Scopus - G. O. Wesolowsky and R. F. Love, “Location of facilities with rectangular distances among point and area destinations,”
*Naval Research Logistics Quarterly*, vol. 18, pp. 83–90, 1971. View at Google Scholar - E. Carrizosa, M. Muñoz-Márquez, and J. Puerto, “The weber problem with regional demand,”
*European Journal of Operational Research*, vol. 104, no. 2, pp. 358–365, 1998. View at Google Scholar · View at Scopus - R. E. Stone, “Some average distance results,”
*Transportation Science*, vol. 25, no. 1, pp. 83–91, 1991. View at Publisher · View at Google Scholar · View at MathSciNet - J. Puerto and A. M. Rodríguez-Chía, “On the structure of the solution set for the single facility location problem with average distances,”
*Mathematical Programming*, vol. 128, no. 1-2, pp. 373–401, 2011. View at Publisher · View at Google Scholar · View at MathSciNet - Z. Drezner and G. O. Wesolowsky, “Location models with groups of demand points,”
*INFOR Journal*, vol. 38, no. 4, pp. 359–372, 2000. View at Google Scholar · View at Scopus - O. Berman, Z. Drezner, and G. O. Wesolowsky, “Location of facilities on a network with groups of demand points,”
*IIE Transactions*, vol. 33, no. 8, pp. 637–648, 2001. View at Publisher · View at Google Scholar · View at Scopus - J. Brimberg and G. O. Wesolowsky, “Note: facility location with closest rectangular distances,”
*Naval Research Logistics*, vol. 47, no. 1, pp. 77–84, 2000. View at Google Scholar · View at MathSciNet - S. Nickel, J. Puerto, and A. M. Rodriguez-Chia, “An approach to location models involving sets as existing facilities,”
*Mathematics of Operations Research*, vol. 28, no. 4, pp. 693–715, 2003. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - R. F. Love and J. G. Morris, “Mathematical models of road travel distances,”
*Management Science*, vol. 25, no. 2, pp. 130–139, 1979. View at Google Scholar · View at Scopus - J. E. Ward and R. E. Wendell, “A new norm for measuring distance which yields linear location problems,”
*Operations Research*, vol. 28, pp. 836–844, 1980. View at Google Scholar - J. E. Ward and R. E. Wendell, “Using block norms for location modeling,”
*Operations Research*, vol. 33, no. 5, pp. 1074–1090, 1985. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - C. Witzgall, “Optimal location of a central facility, mathematical models and concepts,”
*Report*8388, National Bureau of Standards, Washington, DC, USA, 1964. View at Google Scholar - C. Michelot and O. Lefebvre, “A primal-dual algorithm for the Fermat-Weber problem involving mixed gauges,”
*Mathematical Programming*, vol. 39, no. 3, pp. 319–335, 1987. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - R. Durier, “On Pareto optima, the Fermat-Weber problem, and polyhedral gauges,”
*Mathematical Programming*, vol. 47, no. 1, pp. 65–79, 1990. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - E. Carrizosa, E. Conde, M. Munoz-Márquez, and J. Puerto, “Simpson points in planar problems with locational constraints. The polyhedral-gauge case,”
*Mathematics of Operations Research*, vol. 22, no. 2, pp. 291–300, 1997. View at Google Scholar · View at Scopus - S. Nickel, “Restricted center problems under polyhedral gauges,”
*European Journal of Operational Research*, vol. 104, no. 2, pp. 343–357, 1998. View at Google Scholar · View at Scopus - M. Cera, J. A. Mesa, F. A. Ortega, and F. Plastria, “Locating a central hunter on the plane,”
*Journal of Optimization Theory and Applications*, vol. 136, no. 2, pp. 155–166, 2008. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - F. Plastria, “Asymmetric distances, semidirected networks and majority in Fermat-Weber problems,”
*Annals of Operations Research*, vol. 167, pp. 121–155, 2009. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - I. Norman Katz and S. R. Vogl, “A Weiszfeld algorithm for the solution of an asymmetric extension of the generalized Fermat location problem,”
*Computers & Mathematics with Applications*, vol. 59, no. 1, pp. 399–410, 2010. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - H. Minkowski,
*Theorie der Konvexen Körper*, Gesammelte Abhandlungen, Teubner, Berlin, 1911. - Z. Drezner,
*Facility Location: A Survey of Applications and Methods*, Springer, New York, NY, USA, 1995. View at MathSciNet - A. Ghosh and G. Rushton,
*Spatial Analysis and Location-Allocation Models*, Van Nostrand Reinhold, New York, NY, USA, 1987. - L. Cooper, “Solutions of generalized location equilibrium models,”
*Journal of Regional Science*, vol. 7, pp. 1–18, 1967. View at Google Scholar - N. Megiddo and K. J. Supowit, “On the complexity of some common geometric location problems,”
*SIAM Journal on Computing*, vol. 13, no. 1, pp. 182–196, 1984. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - L. Cooper, “Heuristic methods for location-allocation problems,”
*SIAM Review*, vol. 6, pp. 37–53, 1964. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - J.-L. Jiang and X.-M. Yuan, “A heuristic algorithm for constrained multi-source Weber problem: the variational inequality approach,”
*European Journal of Operational Research*, vol. 187, no. 2, pp. 357–370, 2008. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - J.-l. Jiang, K. Cheng, C.-C. Wang, and L.-p. Wang, “Accelerating the convergence in the single-source and multi-source Weber problems,”
*Applied Mathematics and Computation*, vol. 218, no. 12, pp. 6814–6824, 2012. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - Y. Levin and A. Ben-Israel, “A heuristic method for large-scale multi-facility location problems,”
*Computers & Operations Research*, vol. 31, no. 2, pp. 257–272, 2004. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - K. Fukunaga,
*Introduction to Statistical Pattern Recognition*, Academic Press, Boston, Mass, USA, 1990. View at MathSciNet - F. Facchinei and J. S. Pang,
*Finite-Dimensional Variational Inequalities and Complementarity Problems*, Springer, New York, NY, USA, 2004. - M. C. Ferris and C. Kanzow, “Engineering and economic applications of complementarity problems,”
*SIAM Review*, vol. 39, no. 4, pp. 669–713, 1997. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - K. Addi, B. Brogliato, and D. Goeleven, “A qualitative mathematical analysis of a class of linear variational inequalities via semi-complementarity problems: applications in electronics,”
*Mathematical Programming*, vol. 126, no. 1, pp. 31–67, 2011. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - A. Bnouhachem, H. Benazza, and M. Khalfaoui, “An inexact alternating direction method for solving a class of structured variational inequalities,”
*Applied Mathematics and Computation*, vol. 219, no. 14, pp. 7837–7846, 2013. View at Publisher · View at Google Scholar · View at MathSciNet - A. Barbagallo and P. Mauro, “Time-dependent variational inequality for an oligopolistic market equilibrium problem with production and demand excesses,”
*Abstract and Applied Analysis*, vol. 2012, Article ID 651975, 35 pages, 2012. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - J. Gwinner, “Three-field modelling of nonlinear nonsmooth boundary value problems and stability of differential mixed variational inequalities,”
*Abstract and Applied Analysis*, vol. 2013, Article ID 108043, 10 pages, 2013. View at Publisher · View at Google Scholar · View at MathSciNet - Y.-B. Zhao and J.-Y. Yuan, “An alternative theorem for generalized variational inequalities and solvability of nonlinear quasi-${P}_{x}^{M}2a;$-complementarity problems,”
*Applied Mathematics and Computation*, vol. 109, no. 2-3, pp. 167–182, 2000. View at Publisher · View at Google Scholar · View at MathSciNet - H. Uzawa, “Iterative methods for concave programming,” in
*Studies in Linear and Nonlinear Programming*, K. J. Arrow, L. Hurwicz, and H. Uzawa, Eds., pp. 154–165, Stanford University Press, Stanford, Calif, USA, 1958. View at Google Scholar - B. S. He, “A new method for a class of linear variational inequalities,”
*Mathematical Programming*, vol. 66, no. 2, pp. 137–144, 1994. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - B. S. He, “A modified projection and contraction method for a class of linear complementarity problems,”
*Journal of Computational Mathematics*, vol. 14, no. 1, pp. 54–63, 1996. View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - N. Xiu, C. Wang, and J. Zhang, “Convergence properties of projection and contraction methods for variational inequality problems,”
*Applied Mathematics and Optimization*, vol. 43, no. 2, pp. 147–168, 2001. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - B. C. Eaves, “On the basic theorem of complementarity,”
*Mathematical Programming*, vol. 1, no. 1, pp. 68–75, 1971. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet