A trust-region filter-SQP method for mathematical programs with linear complementarity constraints

A trust-region filter-SQP method for mathematical programs with 
linear complementarity constraints (MPLCCs) is presented. The method 
is similar to that proposed by Liu, Perakis and Sun 
[Computational Optimization and Applications, 34, 5-33, 2006] 
but it solves the trust-region quadratic programming subproblems at each 
iteration and uses the filter technique to promote the global convergence. As a result, 
the method here is more robust since it admits the use of Lagrangian 
Hessian information and its performance is not affected by any penalty parameter. 
The preliminary numerical results on test problems generated by the QPECgen generator 
show that the presented method is effective.

1. Introduction. We consider to solve the following mathematical program with linear complementarity constraints (MPLCC): where f : R n+m → R is a twice continuously differentiable real-valued function, C ∈ R l×n , D ∈ R l×m , A ∈ R p×n , B ∈ R p×m , N ∈ R m×n , M ∈ R m×m are given matrices, and r ∈ R l , b ∈ R p , q ∈ R m are given vectors.
Mathematical programs with equilibrium constraints (MPEC) has extensive applications in practical areas such as traffic control, engineering design and economic modeling, e.g., see [21]. It is known to be a class of very difficult optimization problems. In recent years, there have been many works dealing with MPEC. Among these papers, some considered the existence and stationarity of its solution (see [8,24]), and some others proposed algorithms for MPEC, see, for example, [1,8,19,25].
As a special and simple case of MPEC, MPLCC has also attracted some attentions of optimization researchers. By reformulating the complementarity constraints into a system of semismooth equations, Fukushima et. al. [7] presented a sequential quadratic programming (SQP) algorithm for MPLCC. In Fukushima and Pang [8], the feasibility issues on the subproblems in SQP for MPLCC were studied. Fukushima and Tseng [9] proposed an active-set algorithm for MPLCC. Very recently, a primal-dual active set projected Newton method for MPLCC was studied by Chen and Goldfarb [2]. Some works for mathematical programs with complementarity constraints (MPCC), such as [14,17,22,24], may also apply to solve MPLCC.
Based on an inequality relaxation of the complementarity constraints, Liu et. al. [18] presented a robust SQP algorithm for MPLCC. The algorithm has some interesting features, such as, the feasibility of the problem is not required in advance, all subproblems are feasible, and convergence results do not rely on the MPEC-LICQ or the nondegeneracy conditions.
In order to avoid the practical difficulties associated with the setting of the penalty parameter, Fletcher and Leyffer [4] introduced the filter technique for globalizing methods for nonlinear programming, and showed numerically that a trustregion SQP algorithm with the technique performed very promising. Till now, many filter methods for solving nonlinear programming problems have been developed in literature [3,5,10,26,27,28]. Thus, it is a natural consideration to try to develop some effective filter methods for MPEC or MPCC. Actually, a filter method for MPEC has been developed by Leyffer and Munson [16] which solves two subproblems in each iteration, one is a linear program with equilibrium constraints (LPEC) and the other is an equality constrained quadratic program (EQP). The global convergence is obtained through the use of a three-dimensional filter. Long et. al. [20] present a projection-filter method for MPCC.
In this paper, we try to solve MPLCC with filter technique. A trust-region filter-SQP method for MPLCC is then proposed. The method is similar to that proposed by Liu et. al. [18] but it solves the trust-region quadratic programming subproblems at each iteration and uses the filter technique to promote the global convergence. Some numerical results are reported and compared with those in [18] and show that the new algorithm for MPLCC performs well especially for the test problems with non-quadratic objective functions.
The superscript k indicates that z k is generated in kth iteration and so forth. The subscript i means that z i is the ith component of z and C i the ith row of the matrix C. For simplicity, sometimes we use c i (z) and a i (z) to denote the value and the gradient of the ith constraint at z. The sets E, I and C are the index sets of equality constraints, inequality constraints and complementarity constraints, respectively. We refer N * to the neighborhood of z * . The paper is organized as follows. In Section 2, we describe the trust-region filter-SQP algorithm. The global convergence of the algorithm is proved in Section 3 and the preliminary numerical results are reported in the last section.
2. The new trust-region filter-SQP algorithm. Consider the inequality relaxation problem of (P): (N (τ )) where τ ≥ 0 is a relaxation parameter, e = (1 . . . 1) T . It is different from problem (P) in that the complementarity constraints of (P) are relaxed by inequalities (12). Some works, e.g., [15,18,24], have studied the above relaxation problem and have discussed the behavior of the sequence of stationary points of N (τ ) as τ → 0. In this paper, we attempt to solve the problem with filter technique.
The algorithm is iterative. Let z k = (x k , y k , w k ) be the current iteration point satisfying all linear and nonnegative constraints of the original problem and δ k be the current trust region radius. The trust-region quadratic programming subproblem in our algorithm denoted by QP (z k , τ k , δ k ) is as follows.
(QP (z k , τ k , δ k )) Adx + Bdy = 0, where g k = ∇f (x k , y k ), H k can be the Lagrangian Hessian at z k or its approximation. W k = diag(w k 1 , . . . , w k m ), Y k = diag(y k 1 , . . . , y k m ). We use the l ∞ -norm instead of the usual l 2 -norm here for easy computation. For simplicity, we set d = (dx, dy, dw).
It is noted that, if z k is a feasible point of problem (6)−(12), then problem (13)−(20) is feasible since d = 0 satisfies all the constraints of the problem. Then the solution d can be found. Filter technique is used in the algorithm to decide wether to accept z k + d as the new iteration point or not.
When z k is not feasible for problem (6)− (12), which may result in that the subproblem (13)−(20) is inconsistent. If this happens, we employ a restoration phase to find d such that z k + d is acceptable to the filter and to z k (see Definition 2.2), then z k + d will be taken as the new iteration point z k+1 .
Let F k denote the filter at the kth iteration, which means that F k is a subset of {(f i , h i )|i = 0, 1, . . . , k} in which no pair dominates any other one. Let F k = {j|j < k, (f j , h j ) ∈ F k } and π k = min j∈F k {h j }. The following definition is necessary to promote the global convergence of the filter method.
where β and γ are prescribed parameters such that 0 < γ < β < 1. In practice, β is generally selected to be close to 1 and γ close to 0. We also say that z is acceptable to z k .
for all j ∈ F k .
Since our method is a trust region method, we need to calculate the predicted reduction ∆q and the actual reduction ∆f , where and set η = ∆f /∆q, which is the ratio of the actual change to the predicted change in q(d). If d satisfies sufficient descent conditions: ∆f ≥ σ∆q > 0, where σ ≥ 0 is some constant, the resulting iteration is called an f-type iteration. If ∆q ≤ 0 or the subproblem is inconsistent, the purpose of the iteration is to reduce h and we refer the generated iteration as an h-type iteration. Now we introduce the restoration phase. If z k is infeasible,h k > 0, then it is possible that some constraints in (19) fail to be satisfied. In other words, QP (z k , τ k , δ k ) is inconsistent. In this situation, the algorithm need to enter the restoration phase to look for a d with the minimum violation on (19). The restoration techniques are motivated by [18]. We solve the problem A(z k , τ k , δ k ) to get an approximation solution with minimum violation for the relaxed constraints: where χ ∈ (0, 1) is a constant. By introducing some additional variables v ∈ R m + , the above problem can be equivalently transformed to the following linear program We should notice that B(z k , τ k , δ k ) is always consistent because z k satisfies the linear constraints and can be easily obtained by solving B(z k , τ k , δ k ).
In the restoration phase, the following modified quadratic programming problem M QP (z k , τ k , δ k ) is solved: The above subproblem M QP (z k , τ k , δ k ) is always feasible sinced k is a feasible solution. We denote the solution to M QP (z k , τ k , δ k ) byd k = (dx k ,dy k ,dw k ). It can be proved that, if ṽ k 1 −h k = 0, it must give rise to that z k +d k is acceptable to the filter and z k (See Lemma 3.8) as long as the trust region radius is small enough.
Now we are ready to present the trust-region filter-SQP algorithm as follows.

and go to
Step 1; Step , update H k to H k+1 , let k = k + 1 and go to Step 1, else δ k = 0.5δ k and go to Step 4.
As the algorithm enters into the restoration phase,d k andd k are generated. If h k − ṽ k 1 = 0, then update τ k and go to the normal phase. Otherwise, δ k would be decreased till we find ad k such that z k +d k is acceptable to F k and z k .
In this algorithm, the point z k is included in the filter at the end of the iteration if and only if that iteration is an h-type iteration. Hence,h j > 0, j ∈ F k and π k > 0. Ifh j = 0, an f-type iteration would be derived. We can see that, ifh j = 0, and z j is not a KKT point then ∆q > 0 holds. Thus, ifh j = 0, the resulting iteration is an f-type iteration and z j is not added into the filter.
3. Global convergence. For proving the global convergence of the algorithm, we need the following standard assumptions:

Assumption
(1) All the trial points z in the algorithm lie in a non-empty close and bounded set Z.
(2) There exists an M > 0 such that H k 2 ≤ M for all k.
For any feasible point z = (x, y, w) and any relaxation parameter τ , define the index sets and D * = D(z * ). The MPEC-MFCQ plays an important role in proving the global convergence of algorithms for MPEC.
Proposition 1. Let z * solve the MPLCC (P), and assume that (P) satisfies MPEC-MFCQ at z * . Then we have the following results: (1) z * is a feasible point of (P).
(2) For every partition (D y , D w ) of D * , the system is inconsistent, that is, its set of solutions is empty.
Proof. It can be proved similar to Proposition 4.2 in [16].
Lemmas 3.2−3.3 and Corollary 1 have been proved in [6] (see Lemmas 1−2 and the Corollary in [6]). These results are important and still hold for our algorithm for MPLCC. We state them again here for reader's convenience.
Proof. We can obtain the results from Lemma 3 in [6] due to Proof. It is a direct consequence from Lemma 3.4 and h(z k + d) ≤ βπ k .
Lemma 3.6. Let standard assumptions hold. If z * is a feasible point of (P) but not an optimal point, the MPEC-MFCQ holds at z * , then there exist a neighborhood N * of z * and a constant such that, for any z ∈ N * ∩ Z, QP (z, τ, δ) has a feasible solution d which is an f-type step and satisfies for all δ ≤ κ, where κ is a constant related to z, and τ > w • y ∞ .
Proof. Since z * is feasible point of (P) and is not optimal, at which MPEC-MFCQ holds, by Definition 3.1 and Proposition 1, there exist d with d = 1 and a partition (D y , D w ) of D * such that (63)−(69) hold. By the continuity, there exists a neighborhood N * of z * and > 0 such that where g is evaluated for any z ∈ N * ∩G and a i is a constant vector because (56)−(58) is linear. Part A: Let d α = αδd, where α ∈ [0, 1]. Thus, d α ∞ ≤ δ, that is, d α satisfies the trust region constraint of QP (z, τ, δ). For i ∈ E, since z ∈ N * ∩ G, c i (z) = 0. c i (z) + d T α a i = c i (z) + αδd T a i = 0, so d α satisfies the equality constraints (15)−(16) of QP (z, τ, δ).

Now let us look at the inactive inequality constraints i /
∈ A * I . By the continuity of c i and the boundedness of a i , there exist constantsc 1 andā such that c i (z) ≤c 1 and a T i d ≤ā, (18) could be handled by the same arguments as above. Next we consider the relaxed constraint (19) of QP (z, τ, δ). Notice that d α = (αδdx, αδdy, αδdw). For i ∈ D y , because d = 1 and α ∈ [0, 1],

Applying Lemma 3.3 and Lemma 3.4, the minimum value of Φ(α) occurs at
if δ ≤ σ/2γm. To sum up, for δ ≤ κ, δd is a feasible point of QP (z, τ, δ) and is an f-type step, where Lemma 3.7. Let the standard assumptions hold. Then the inner iteration loop (step1 → step 2 → step1) of QP (z k , τ k , δ k ) terminates in a finite number of iterations.
Proof. If d = 0 solves QP (z k , τ k , δ k ), then z k is a KKT point of N (τ k ), and the algorithm is continued with τ k updated. Otherwise, if the inner iteration loop does not terminate in a finite number of iterations, by the algorithm, δ k → 0. Then two cases need to be considered, depending on wetherh k = 0 orh k > 0.
(1)h k = 0. Obviously, d = 0 is a feasible point of QP (z k , τ k , δ k ) and z k is a feasible point of N (τ k ). Hence, if z k is not a KKT point of N (τ k ), by the similar arguments as that in the proof of Lemma 3.6, there have a κ and a feasible solution d of QP (z k , τ k , δ) such that d is an f-type step and acceptable to filter for δ, where δ ≤ min{κ, βπ k /m}. Becauseh k = 0, τ k > w k • y k ∞ and κ > 0, a suitable δ could be located and the inner iteration loop terminates finitely.
(2)h k > 0, there exists an i, i ∈ C such that w k i y k i > τ k . Thus, Then inner iteration loop terminates and the algorithm starts the restoration phase because the ith constraint does not hold.
Lemma 3.8. In the restoration phase, if ṽ k 1 −h k = 0, thend k generated by M QP (z k , τ k , δ k ) is always an acceptable step when δ k is small enough.
Proof. The algorithm enters the restoration phase because of the violation on (19), thush k > 0. By the constraints of B(z k , τ k , δ k ), for α ∈ (0, 1) We can observe that if dk < u : Thus, h(z k +d k ) ≤ βh k holds for a proper β ∈ (0, 1). In this way, for any δ k < min{u, βπ k /m}, z k +d k can be accepted by the filter and z k from Lemma 3.5.
Theorem 3.9. The algorithm terminates with one of the following outcomes: (1) A KKT point of N (τ k ) is found, where τ k is less than the tolerance termination; (2) Let z * be any accumulation point of the sequence {z k }. If the MPEC-MFCQ holds at z * and the algorithm enters the restoration phase in a finite number of iterations, then z * is a strong stationary point of (P); (3) If MPEC-MFCQ does not hold at z * and the algorithm enters the restoration phase in a finite number of iterations, then z * is a singular stationary point; (4) The algorithm enters the restoration phase again and again, and there exists an accumulation point z * that is an infeasible stationary point.
Proof. (1) The first outcome corresponds to the normal termination of the algorithm.
(2) By Lemma 3.7 and 3.8, the inner loop in each iteration is always terminated in a finite number of iterations. All iterates lie in Z, and Z is bounded, so the sequence has one or more accumulation points.
Case A: The sequence {z k } contains an infinite number of h-type iterations and we consider this subsequence. For simplicity, this subsequence is still denoted as {z k }. Because in an h-type iteration (f k , h k ) is always added into the filter, so it follows by Corollary 1 that h k → 0. Only h-type iterations can reset π k , so there exists a subsequence {z kj } of {z k }, on which π kj+1 = h kj < π kj is set. Because Z is bounded, there exists an accumulation point z * and a subsequence {z kj i , i ∈ I}(I is the index set) of {z kj } on which z kj i → z * , h kj i → 0 and π kj i+1 = h kj i < π kj i . So z * is a feasible point. For sufficient large i, i ∈ I, z kj i ∈ N * ∩ Z as defined in Lemma 3.6. If z * is not a KKT point and the MPEC-MFCQ holds, by Lemma 3.5 and 3.6 there exists a QP (z kj i , τ, δ) which can generate an f-type step d and d is acceptable to the filter for This contradicts the fact that the filter includes only h-type step. Case B: The main sequence {z k , k ∈ K} ( K is the index set) contains only a finite number of h-type iterations. There exists an index k 0 such that all iterations are f-type iterations and π k = π k0 which is a constant for all k ≥ k 0 . It follows that (f k+1 , h k+1 ) is always acceptable to (f k , h k ) and ∆f ≥ σ∆q ≥ 0 so that f k is strictly monotonically decreasing for k > k 0 . So h k → 0 and any accumulation point is feasible from Lemma 3.2. Because f (x) is bounded on Z it follows that k ∆f k is convergent and z k → z * . For sufficient large k > k 0 , z k ∈ N * ∩ Z as defined in Lemma 3.6. If z * is not a KKT point of (P) and the MPEC-MFCQ holds, by Lemma 3.5 and 3.6, there exists a QP (z k , τ, δ) which can generate an f-type step d and d is acceptable to the filter for From Lemma 3.6 which contradicts the fact that k ∆f k is convergent. Hence, by the definition 4.1 in [18], z * is a strong stationary point of (P).
(3) From (2) we can notice that the limit point is a feasible point of (P) even if the MPEC-MFCQ does not hold. So there is a limit point that is a singular stationary point [18].
(4) If the algorithm enters the restoration phase again and again, there exists a subsequence {(f k , h k )} ∈F k such that ṽ k 1 −h k = 0. From (82) and (85) Since ∞ i=0 τ k = τ 0 1−σ < +∞ and (82) holds, it can be seen that h k andh k are convergent. For the boundness of Z, z k converge to z * . As k → ∞, τ k → 0, we can obtain (y * ) T w * > 0. Hence this limit point is an infeasible stationary point [18]. 4. Numerical results. We implemented the algorithm in MATLAB version 6.5 on a LENOVO personal computer with a WINDOWSxp system. During the implementation, we generated the same initial point as [18] and used the same techniques as the subproblem (141)-(146) in [18] for coping with the termination errors in solving problem B(z k , τ k , δ k ). For comparison to [9,18], our test problems also included mathematical programs with quadratic and non-quadratic objective functions. The programs with quadratic objective functions were generated by a generator named QPECgen presented in [12], and the non-quadratic objective functions were derived by adding a cubic function 1 3 to the corresponding quadratic objective functions. The parameters were chosen as the following: τ 0 = max{1, (y 0 ) T w 0 /m}, χ = 0.9, ξ = 0.1,δ = 10 3 , δ = 10 −3 , σ = 0.3, β = 0.99, γ = 0.25, = 5 × 10 −7 and u = 1000. We report our numerical results on the QPECgen problems in Table 1, Table  2 and Table 3. The columns in iL and f * L are respectively the number of iterations and the optimal value presented in [18], and i, f * are those obtained by the trust-region filter SQP algorithm in this paper. The columns in f gen and f 0 are respectively the values of objective function at (x gen , y gen ) and the initial point. s d and m M are "second deg" and "mono M " for short, they are the same as those in [18]. Table 1 shows that: 1. Liu et. al. [18] solves 9 test problems successfully, while the trust-region filter SQP algorithm does 11.
2. The new algorithm takes less iterations. 3. In order to further observe the performance of the algorithm in solving TP12, we set x 0 = x gen , y 0 = y gen and w 0 = N x 0 + M y 0 − q as done in [18]. From the last line of Table 1, we can see that a much smaller optimal value than f gen has been obtained. Thus, for quadratic test problems, the new algorithm is competitive to that of [18].
In Table 2, we report the residues on these test problems, where N orm1 = (x 0 , y 0 ) − (x gen , y gen ) ∞ , N orm2 = (x * , y * ) − (x gen , y gen ∞ , N orm3 = ((Cx + Dy − c) + , Ax + By − b, N x + M y − w − q, y + , w + ) , N orm4 = l 2 norm of residues of KKT conditions of N (τ * ), Obviously, Norm4 measures the optimality of the solution to N (τ * ) for τ * ≤ , Norm3 and Norm5 are respectively measures of the maximal violation on the linear constraints and the complementarity constraints, Norm5 examines the feasibility of the derived solution to the original problem. We report in Table 3 our numerical results on the MPLCCs with non-quadratic objective (NP1-12). The initial points were the same as problems TP1-12.
We can note that: 1. The new algorithm takes less iterations for most of the test problems. 2. For each problem of NP1-12, the new algorithm derived a much smaller optimal value than that of [18].
3. By comparing Table 1 with Table 3, we can see that the new algorithm performs much better for non-quadratic objective than for quadratic objective.
We can observe the process of the algorithm by taking TP2 as an example and noticing that where Knorm and Krest record the iteration number in the normal phase and the iteration number in the restoration phase, respectively. For example, the three 10s in Knorm show that the trust region radius was decreased three times in the 10th iteration before a step was accepted. The three 1s in Krest mean that the restoration phase finded an acceptable step after the trust region radius was decreased three times in the first iteration. It seems that the restoration phase makes sense in the algorithm and it is not started when the algorithm is approaching the limit point. Similar results can also be observed in other testing problems.