On linear problems with complementarity constraints

A mathematical program with complementarity constraints is an optimization problem with equality/inequality constraints in which a complementarity type constraint is considered in addition. This complementarity condition modifies the feasible region so as to remove many of those properties that are usually important to obtain the standard optimality conditions, e.g., convexity and constraint qualifications. In this paper, in the linear case, we introduce a decomposition method of the given problem in a sequence of parameterized problems, that aim to force complementarity. Once we obtain a feasible solution, by means of duality results, we are able to eliminate a set of parameterized problems which are not worthwhile to be considered. Furthermore, we provide some bounds for the optimal value of the objective function and we present an application of the proposed technique in a non trivial example and some numerical experiments.


Introduction
In the field of equilibrium models, mathematical programs with complementarity constraints (MPCC) form a class of important, but extremely difficult, problems. MPCC's constitute a subclass of the well-known mathematical programs with equilibrium constraints (MPEC), widely studied in recent years. Their relevance comes from many applications like urban traffic control, economy, problems arising from the electrical sector or from structural engineering [4][5][6]26]. Their difficulty is due to the presence of the complementarity constraints, because the feasible region may not enjoy some fundamental properties: it may be not convex, even not connected and such that many of the standard constraint qualifications are violated at any feasible point. This last lack implies that the usual Karush-Kuhn-Tucker (KKT) conditions may not be fulfilled at an optimal solution, even in the linear case [15], denoted with LPCC. Specific constraint qualifications (CQs) for the MPCC were introduced to try to overcome this difficulty [8,9,15,22,23,25,30] and different notions of stationarity were also proposed [23]. Using these concepts, the behaviour of general nonlinear optimization algorithms for the MPCC was studied. For example, interior point methods [18,21], penalty approaches [12,20,28,29], relaxation methods [9,13,16,24,25], smoothing methods [1,7,27]. Also, many specific methods which combine the previous approaches for the MPCC were developed [11,14,19]; such methods have good convergence properties, but they are affected by the request of exact computation of KKT points [17].
The method that we propose is different from most of the aforementioned ones; the main difference is that it avoids the use of CQs and KKT conditions. In this paper we consider the linear case and a reformulation of the problem by means of a family of parameterized linear problems whose minimization leads to an optimal solution of LPCC. Exploiting the classic tools of the duality theory, we propose an iterative method which explores the set of parameters, excluding at each step a subset of them, by means of a suitable cut; indeed, the optimal values of the linear problems associated with such a subset are proved to be greater than or equal to the optimal value related to the current parameter. A similar approach for LPCC can be found in [10,31] where different kinds of cuts are considered. We define an algorithm which can be implemented in an interactive way taking advantage of some devices that speed up the solving procedure, owing to the decomposition of the given problem in a sequence of parameterized problems.
The paper is organized as follows. In Sect. 2, we introduce the problem and describe the decomposition method of LPCC in a family of parameterized problems. In Sect. 3, we exploit the classic tools of duality theory to obtain a necessary condition for improving the current solution and, moreover, we establish a sufficient optimality condition. In Sect. 4, we obtain upper and lower bounds of the optimum value of the problem that define a sequence of nested intervals by which the iterative procedure is stopped when their length is below a fixed tolerance. Finally, in Sect. 5, the iterative method is described and illustrated by means of an example and some numerical experiments involving problems with different dimensions are presented. The concluding Sect. 6 is devoted to some final remarks and the description of further developments.

A decomposition approach
We consider the following constrained minimization problem, whose objective function is linear and having a linear complementarity constraint besides affine ones: where A, B ∈ R m×n , c, d ∈ R n , b ∈ R m and ·, · denotes the scalar product in R n .

Remark 1
We observe that we can consider in the previous formulation equality constraints, provided that we replace them by a couple of inequality constraints. Moreover, in our setting, the feasible region includes the classic complementarity system [3] of the form z ≥ 0, q + Mz ≥ 0, z, q + Mz = 0. In fact, by setting y = q + Mz, the previous system is equivalent to Mz − y ≥ −q, −Mz + y ≥ q, z ≥ 0, y ≥ 0, z, y = 0.
We will assume that the feasible set K is nonempty and a global minimum point of P exists; call it (x 0 , y 0 ). Let us introduce the following penalized form for the gradient of the objective function of P: where α = (α 1 , . . . , α n ) ∈ Δ := {0, 1} n . For the sake of simplicity, we have used for the penalized form of the constants c j and d j (that are functions of α j ) the same name of the constants. For our purposes, we will assume that ρ j and σ j , j = 1, . . . , n, are large enough positive constants; the meaning of this assumption will be clear inside the proof of Theorem 1.
The given problem P can be associated with a family {P(α)} α∈Δ of subproblems Assumption 1 Let us suppose that, for large enough ρ j and σ j , j = 1, . . . , n, the objective function of P(α) is bounded from below on R, ∀α ∈ Δ.
By Assumption 1, it follows that ∀α ∈ Δ there exists a minimum point, say (x(α), y(α)), of P(α). A condition that guarantees Assumption 1 is that c, x + d, y is bounded from below on R, which obviously implies that the objective function of P is bounded from below on K , which in turn yields that a global minimum point of P exists.
We have the following result.

Theorem 1 If Assumption 1 holds, then
Proof Suppose that (x 0 , y 0 ) is a minimum point of P, so that f 0 = c, x 0 + d, y 0 .
Recalling that Δ is a finite set, let i = 1 and α i ∈ Δ := {α 1 , α 2 , . . . , α 2 n }. From the definition of c(α) and d(α), we have Let ver t R be the set of vertices of R.
(which happens independently on ρ and σ ), then x,ȳ = 0 and (x,ȳ) is a feasible solution of P. Therefore, and there exists j such that ρ i jx j > 0 (or σ i jȳ j > 0). Now, if and this means f 0 ≤ f (x,ȳ; α i ). Analogously, in the case σ i jȳ j > 0. Noticing that ver t R is a finite set, it is possible to find a couple (ρ i , σ i ) such that From Assumption 1 (by further increasing some components of ρ i and σ i , if necessary) a minimum point of P(α) exists and is attained in the set ver t R. From (3), it follows f 0 ≤ f ↓ (α i ). Now, consider α i+1 ; if f 0 ≤ f (x, y; α i+1 ), ∀(x, y) ∈ ver t R, with ρ = ρ i and σ = σ i , then set ρ i+1 = ρ i and σ i+1 = σ i ; in such a way, we obtain f 0 ≤ f ↓ (α i+1 ). Otherwise, choose ρ i+1 and σ i+1 by increasing the previous vectors ρ = ρ i and σ = σ i , in order to get f 0 ≤ f ↓ (α i+1 ). Set i = i + 1 and repeat this procedure for all i = 2, .., 2 n − 1. It turns out that Moreover, starting again from (x 0 , y 0 ) minimum point of P, let us define the following vector α 0 = α(x 0 , y 0 ): The above vector α 0 belongs to Δ and f 0 = f (x 0 , y 0 ; α 0 ). Since (x 0 , y 0 ) ∈ R is a feasible solution of the problem P(α 0 ), then Inequalities (4) and (5) imply (1) and this concludes the proof.
Based on Theorem 1, we can propose a decomposition approach for solving problem P. In fact, Theorem 1 establishes that the minimum f 0 of problem P can be achieved by first determining f ↓ (α) ∀α ∈ Δ, and secondly by minimizing f ↓ (α) with respect to α ∈ Δ; in other words, Eq. (1) decomposes the problem P in a sequence of subproblems P(α).
In view of the definition of α, ρ j and σ j , j = 1, . . . , n, and of problem P(α), by setting α j = 1 or α j = 0 one would expect x j = 0 or y j = 0, respectively, in an optimal solution of P(α); however, in general, an optimal solution of P(α) will not necessarily comply with such an expectation for a given α ∈ Δ. Anyway, if such an expectation does not occur for a given α ∈ Δ, then f ↓ (α) > f 0 , provided that we choose ρ and σ large enough. Consequently, we should work with the subsetΔ of Δ whose elements fulfill the following definition. Definition 1 LetΔ be the set of all α ∈ Δ such that the system is possible, when one sets x j = 0 for α j = 1 and y j = 0 for α j = 0, j = 1, . . . , n.
Clearly, if a minimum point of P exists, then the setΔ is nonempty.
Unfortunately, the cardinality of Δ, and also ofΔ, is in general so large that the above outlined decomposition is not, by itself, of use. We aim at solving P through the penalized problems P(α)'s, by running as less as possible on α ∈ Δ: at step k, having solved P(α k ), we try to determine α k+1 , such that An initial effort in this direction is described in the first part of next section.

Optimality conditions
Suppose that we have solved, at step k, the problem P(α k ); we try to determine α k+1 , such that (6) holds. To this aim, we introduce the dual problem of P(α) given by: Proof Observe thatλ is a feasible solution of P * (α k+1 ) so that the maximum of such a problem is greater than or equal to that of P * (α k ): By the strong duality theorem, the minimum of P(α k+1 ) is greater than or equal to that of P(α k ).
By the previous theorem, we infer that a necessary condition for (6) to hold is that In other words, having an optimal solution of P(α k ) for some α k ∈ Δ, and, thus, a feasible solution to the given problem P, a necessary condition for obtaining a problem P(α k+1 ) with a minimum f ↓ (α k+1 ) < f ↓ (α k ) (namely, a feasible solution of P "better" than the current one), is that α k+1 belongs to the set Δ\Δ k . Theorem 2 provides a criterion for eliminating the subset Δ k from the subsequent iterations; observe that Δ k cannot be empty as it contains at least α k .
We need to deepen the analysis of (7); to this aim, let us reconsider the inequalities in the definition of Δ k :λ Taking into account the definition of c(α) and d(α), the two inequalities are equivalent to: where A j and B j are the j-th column of A and B, respectively.

Remark 2
If the left-hand side of the j-th inequality of (8a) is non positive, then the inequality is fulfilled for any choice of α j , i.e. both for α j = 0 and α j = 1. Otherwise, let us suppose that it is positive; in such a case, since ρ j > 0 is large enough, we can assume -without any loss of generality -that Hence, the j-th inequality of (8a) is fulfilled if α j = 1, but not if α j = 0. Obviously, we can make analogous remarks on the j-th inequality of (8b) and we obtain that if Let us introduce the following sets of indexes where, both in (9) and in (10), the second equality is true because ρ j > 0 and σ j > 0, j = 1, . . . , n. Now we establish some results expressed in terms of the sets I x (λ) and I y (λ). In the next results we assume that α k ∈ Δ andλ is a maximum point of P * (α k ).
Proof Inequalities (8) can be equivalently written as: They are fulfilled by α k becauseλ is an optimal solution of P * (α k ). Ab absurdo, if I x (λ) ∩ I y (λ) = ∅, then there exists at least one index j ∈ I x (λ) ∩ I y (λ) such that the j-th inequalities of (11) are simultaneously satisfied neither by α k j = 0 nor by α k j = 1 and this is a contradiction.

Theorem 3 (sufficient optimality condition)
If then an optimal solution (x k ,ȳ k ) of the problem P(α k ) is an optimal solution of problem P. Proof Then all α ∈ Δ satisfy the inequalities (11); i.e., Δ = Δ k which implies Δ\Δ k = ∅. Hence, by Theorem 2, there is no way to improve the current solution (x k ,ȳ k ), so that it is an optimal solution of P.
The following result establishes a condition equivalent to the necessary condition (7). Notice that we are interested only in the case where I x (λ) ∪ I y (λ) = ∅. Indeed, if I x (λ) ∪ I y (λ) = ∅, the optimality of the current solution is proved by Theorem 3.
Proof Let us consider the inequalities in (13) in the equivalent form (8). Therefore (13) holds if at least one of the inequalities of (8) is not fulfilled. Observe that, if i / ∈ I x (λ) or j / ∈ I y (λ), then the i-th inequality in (8a) or the j-th in (8b) is fulfilled whatever α i or α j may be, respectively. Then the system (8) is equivalent to the following The above system is possible iff which coincides with (14).
Define the relaxed problem of P obtained by dropping the complementarity constraints, i.e.
and denote byf the optimal value of R P (possibly −∞). The next result deepens the meaning of the sufficient optimality condition given by Theorem 3.

Proof
The assumption I x (λ) ∪ I y (λ) = ∅ is equivalent to say thatλA ≤ c ,λB ≤ d. Therefore, sinceλ ≥ 0, we have thatλ is a feasible solution for the dual of R P: if f = −∞, we achieve a contradiction. Suppose thatf > −∞; then, by strong duality f ↓ (α) = λ , b and, by weak duality, On the other hand, f 0 ≥f and this completes the proof.

Lower and upper bounds for the minimum
The sufficient optimality condition (12) given in Theorem 3 is a very restrictive condition. Indeed, it directly implies that the minimum value of the relaxed problem R P is equal to the one of P, as proved by Proposition 2.
Therefore, in this section, we propose an alternative iterative approach that leads, not only to a different sufficient optimality condition, but mainly to the possibility to evaluate the difference between the current value of the objective function of P and its minimum value, i.e., f ↓ (α k ) − f 0 . We will achieve this purpose by defining a finite sequence of lower and upper bounds of the minimum of P.
We restrict our attention to the case where it is possible to find vectors of upper bounds, say X = (X 1 , . . . , X n ) and Y = (Y 1 , . . . , Y n ), for x and y respectively, in such a way that the optimal value of problem P does not change. Therefore, alternatively to {P(α)} α∈Δ , we can associate with the given problem P the following family {Q(α)} α∈Δ of subproblems: We assume that Q(α) admits a solution, ∀α ∈ Δ. Denote by S(α) the feasible set of Q(α). The dual of Q(α), say Q * (α), is given by: An optimal basic solution of Q * (α) can be immediately derived from an optimal basic solution of P * (α), as proved by the following proposition.
Letλ be an optimal basic solution of P * (α) and let Γ := (A 1 , B 1 , −I L ) with det Γ = 0 be a basis matrix of order m related toλ, i.e., such that where B 2 ) and, accordingly to the previous partitions, ; moreover, I L is the submatrix of the m−identity matrix related to the columns j, such thatλ j = 0. Note that, sinceλ is an optimal solution of P * (α) and we have supposed ρ j and σ j large enough, then the constraints (16), being fulfilled as equalities, do not contain ρ j and σ j . Therefore, (16) can be written asλ Therefore, in correspondence of the equalities (17), according to (15) we haveμ 1 = 0,ν 1 = 0, whereμ := (μ 1 ,μ 2 ),ν := (ν 1 ,ν 2 ). We now note that (λ,μ,ν) = (λ,μ 1 ,μ 2 ,ν 1 ,ν 2 ) fulfills the following equalities in the constraints of Q * (α): Note that H 1 is a block triangular matrix (where the first block is given by Γ ) with det H 1 = ± det(A 1 , B 1 , −I L ) = ± det Γ = 0, which yields det H = 0, and this shows that (λ,μ,ν) is a feasible basic solution of Q * (α). Next we show that (λ,μ,ν) is also an optimal solution of Q * (α). We first note that ∀α ∈ Δ we have that in other words, λ , b is an upper bound of the feasible values of the objective function of Q * (α). Now, we will prove that the objective function of Q * (α) assumes at (λ,μ,ν) exactly the value λ , b and hence (λ,μ,ν) is an optimal solution of Q * (α). To this aim, we rewrite (15) as follows: When α j = 0, ν j does not affect the objective function of Q * (α); as we have already observed, if α j = 0, the j-th inequality ofλA ≤ c(α) in P * (α) becomesλA j ≤ c j and by the second implication of (18a) we obtainμ j = 0. Similar results hold for α j = 1: μ j does not affect the objective function of Q * (α) andν j = 0. In conclusion, in the objective function of Q * (α), we have n j=1μ j X j (1−α j )+ n j=1ν j Y j α j = 0 and this completes the proof.
Let us observe that the feasible set of Q * (α) does not depend on α; call this set S * , ∀α ∈ Δ, and notice that S * = ∅ by Proposition 3. Denote by V := ver t S * , the set of all vertices of S * , or, equivalently, of all basic solutions of Q * (α). Then, we have the following result.

Theorem 5
The minimum f 0 of problem P equals the minimum of the problem: Proof The following relations are readily seen to hold: The last equality and the introduction of the scalar variable f prove the thesis of the theorem.

Remark 3
Observe that ifV is any subset of V , then by (20) we have: Therefore, the minimum in (21) is a lower bound t of f 0 . At every α met in the solution process, an optimal basic solution of Q(α) and, hence, an optimal basic solution of Q * (α) is available. Let V k be the set of the basic solutions of Q * (α k ) considered until step k. Thus, V k which is initially empty, gains a new element. Every time this happens, problem (21) may be solved to find a new lower bound on f 0 , say it t k . Moreover, the minimum f ↓ (α k ) of problem P(α k ) for every α k ∈ Δ, is obviously an upper bound on the optimal value f 0 ; let T k be the minimum of all the previously found upper bounds. Obviously, at k-th step, the equality T k = t k is a sufficient condition for optimality; moreover, |T k − t k | is un upper bound of the difference between the current value of the objective function of P and its minimum. We can decide to stop the iterative procedure if such a difference is small enough, in a sense that can be specified case by case according to the meaning of the given complementarity problem.

The iterative method
The analysis developed in the previous sections allows us to define an iterative method for the minimization of the problem P.
Step (2) (computation of a cut on the set Δ).
Determine an optimal solution (λ,μ,ν) of the problem Q * (α k ) and let t the minimum in (21) obtained by adding the (basic) solution (λ,μ,ν) to the setV . If t > t k−1 , then t k = t; otherwise t k = t k−1 . If T k = t k then (x(α k ), y(α k )) is an optimal solution of P; hence → STOP. Otherwise, set k = k + 1 and go to Step 1.

Remark 4
The following observations are worth noting: (a) In the solution of problem P(α), we need to choose the values of the vectors ρ = (ρ 1 , ρ 2 , . . . , ρ n ) and σ = (σ 1 , σ 2 , . . . , σ n ) to penalize the costs c(α) and d(α). To obtain such a penalization, the components of ρ and σ must be of order of magnitude greater than that of the elements of c and d. In the subsequent Example 1, any components of ρ and σ will be set equal to 10 10 . (b) In the formulation of problem Q(α) and hence of its dual Q * (α), we need to choose the values of the upper bounds X = (X 1 , . . . , X n ) and Y = (Y 1 , . . . , Y n ) of which we have assumed the existence. The selection of such vectors is a crucial aspect of the method. Even if it is valid only for the particular case, a suggestion for this choice is given in Example 1 at Step 3 of Iteration 1. (c) If, at a certain iteration k, we can prove that f ↓ (α) ≥ f ↓ (α k ) ∀α ∈ Δ k+1 , then f ↓ (α k ) is the optimal value of P and the current solution (x(α k ), y(α k )) is an optimal solution. Suppose for example that |I x (λ)| + |I y (λ)| = 1; from the inequality (14) it follows that there is a unique index j such that any α to be considered in the sequel has α j = 0 if |I x (λ)| = 1, and α j = 1 if |I y (λ)| = 1. If, by solving R P with the additional condition x j = 0 when α j = 1 or y j = 0 when α j = 0, we obtain a minimum greater than or equal to f ↓ (α k ), the current solution is an optimal one. Otherwise, such a minimum is a lower bound of f 0 and it will replace the current lower bound, computed at Step 3, if it is better. (d) Recall that in the proposed decomposition method we should work with the subset Δ of Δ, introduced in Definition 1. If a vector α / ∈Δ, the optimal value of the corresponding P(α) is of the same magnitude of ρ j 's and σ j 's. In this case, we skip Step 1 and we generate a new α. We refer to this case as a null step.
Example 1 Let us apply the iterative method to the following problem P: For the solution of some of the steps, the numerical software MATLAB has been used.
Go to Step 3.
By using MATLAB we have developed some numerical tests. The decomposition method requires to process the set of binary vectors α ∈ Δ; therefore, a crucial aspect is the choice of α k ∈ Δ k in Step 1. Since Δ k is given by a system of inequalities, each one defined by (14) (see Theorem 4), to choose a new α k ∈ Δ k we solve such a system by using the MATLAB routine INTLINPROG, where the objective function is a constant or is arbitrary.
We have tested randomly generated problems of different dimensions (m is the number of linear constraints and 2n is the number of total variables). In order to generate the test problems we used the MATLAB function RAND for the random construction of matrices A, B and vectors c, d, b. With the aim of avoiding very dense matrices we converted to zero a fixed proportion of elements. The results of the tests with n = 50, 100, 200 and m = n/2 are shown in the following Tables 1, 2 and 3 which report in thef column the optimal value of the relaxed problem R P, in the f val column the optimal value of the complementarity problem P, in the niter column the number of iterations the method took to terminate (a maximum number of 350 iterations has been set), in the null column the number of iterations that do not generate a cut (null steps) and in the time column the total computational time. Tables 1, 2 and 3 may be compared with Tables 6.1  when the set Δ k turns out to be empty; for two problems (Problem n.1 in Table 2 and Problem n.8 in Table 3) the maximum number of iterations is reached. We observe that the number of iterations requested by our method is in general comparable with the number of LPs considered in Tables 6.1-6.2-6.3 of [10] and with the number of nodes reported in Tables 5-6 of [31]. We point out that for problems where m ≈ n the convergence is not always reached. In Table 4 we report the results of some experiments with m = 40 and n = 50. In half of the cases we have convergence, while in two tests (5)(6)(7)(8)(9) of the other half, where the maximum number of iterations is reached, we obtain a good approximation of the optimal value. For the remaining three cases (2-4-7) the algorithm is not successful; we remark that in such cases there is a significant number of null steps.

Conclusions
We have introduced a decomposition method for a linear problem with complementarity constraints in a sequence of parameterized problems. By means of suitable cuts we have proposed an iterative method that leads to an optimal solution of the given problem or to an approximation of it providing an estimate of the error. For problems of different dimensions we have implemented some numerical experiments which show that in most of the cases the method converges linearly w.r.t. the dimension n of the problem. A completely unified MATLAB code fully implementing the whole algorithm is still in progress. This is a possible further development, together with some more numerical testing.