A sequential quadratic programming algorithm that combines merit function and filter ideas

Abstract

A sequential quadratic programming algorithm for solving nonlinear programming problems is presented. The new feature of the algorithm is related to the definition of the merit function. Instead of using one penalty parameter per iteration and increasing it as the algorithm progresses, we suggest that a new point is to be accepted if it stays sufficiently below the piecewise linear function defined by some previous iterates on the (f,||C||2²)-space. Therefore, the penalty parameter is allowed to decrease between successive iterations. Besides, one need not to decide how to update the penalty parameter. This approach resembles the filter method introduced by Fletcher and Leyffer [Math. Program., 91 (2001), pp. 239-269], but it is less tolerant since a merit function is still used. Numerical comparison with standard methods shows that this strategy is promising.

sequential quadratic programming; merit functions; filter methods

A sequential quadratic programming algorithm that combines merit function and filter ideas^*

Francisco A.M. Gomes^*

Department of Applied Mathematics, IMECC, Universidade Estadual de Campinas 13081-970, Campinas, SP, Brazil. E-mail: [e-mail]

ABSTRACT

A sequential quadratic programming algorithm for solving nonlinear programming problems is presented. The new feature of the algorithm is related to the definition of the merit function. Instead of using one penalty parameter per iteration and increasing it as the algorithm progresses, we suggest that a new point is to be accepted if it stays sufficiently below the piecewise linear function defined by some previous iterates on the (f,||C||₂²)-space. Therefore, the penalty parameter is allowed to decrease between successive iterations. Besides, one need not to decide how to update the penalty parameter. This approach resembles the filter method introduced by Fletcher and Leyffer [Math. Program., 91 (2001), pp. 239-269], but it is less tolerant since a merit function is still used. Numerical comparison with standard methods shows that this strategy is promising.

Mathematical subject classification: 65K05, 90C55, 90C30, 90C26.

Key words: sequential quadratic programming, merit functions, filter methods.

1 Introduction

In this paper we are concerned with the problem

where f : ⁿ® is a C² nonlinear function, C : ⁿ® ^m represents a set of C² nonlinear constraints and we suppose that

-¥ < l_i < u_i < ¥, for i = 1,...,n.

Naturally, some of the components of x in (1) may be slack variables generated when converting inequality constraints to this form.

Algorithms based on the sequential quadratic programming (SQP) approach are among the most effective methods for solving (1). Some interesting algorithms of this class are given, for example, in [3, 4, 10, 21]. A complete coverage of such methods can be found in [6, 20].

Since SQP algorithms do not require the iterates to be feasible, they have to concern with two conflicting objectives at each iteration: the reduction of the infeasibility and the reduction of function f. Both objectives must be taken in account when deciding if the new iterate is to be accepted or rejected. To make this choice, most algorithms combine optimality and feasibility into one single merit function. Filter methods, on the other hand, usually require only one of these two goals to be satisfied, avoiding the necessity of defining a weighted sum of them. Both approaches have advantages and drawbacks. In this paper, these approaches are mixed with the objective of reducing their disadvantages while keeping their best features. To motivate the new algorithm, both methods are briefly introduced below.

1.1 A prototypical SQP algorithm

The algorithm presented in this paper will inherit the structure of the merit function based SQP method proposed in [10], hereafter called the GMM algorithm.

This method divides each iteration into two components, a normal and a tangential step. The first of these components is used to reduce the infeasibility, while the aim of the second is to improve the objective function. A trust region approach is used as the globalization strategy. All of the iterates are required to satisfy the bound constraints l < x < u, so the merit function chosen was the augmented Lagrangian, written here (and in [10]) in an unusual way as

In (2), q Î [0,1] is a ''penalty parameter'' used as a weight to balance the Lagrangian function (for the equality constrained subproblem), defined as

(x, l) = f(x) + C(x)^Tl,

with a measure of the infeasibility, given by

At iteration k, a new point x₊ = x_k + s is accepted if the ratio between the actual and the predicted reduction of the merit function (when moving from x_k to x₊) is greater than a positive constant.

The actual reduction of the augmented Lagrangian at the candidate point x₊ is defined as

A_red(x_k, s, q) = (x_k, q) - (x_k + s, q).

The predicted reduction of the merit function depends on the strategy used to approximately solve (1). The GMM algorithm approximates (1) by the quadratic programming problem

where H is a symmetric n×n matrix and A(x) = (ÑC₁(x), ¼, ÑC_m(x))^T is the Jacobian of the constraints.

In this case, denoting

as the approximation of j(x), the predicted reduction of the augmented Lagrangian merit function is given by

where

is the predicted reduction of the infeasibility and

is the predicted reduction of the Lagrangian.

The penalty parameter q plays a crucial role in the acceptance of the step. For a augmented Lagrangian of the form (2), (q - 1)/q can be viewed as the slope of the line that defines the forbidden region in the (j, )-plane, that is, the semi-space that contains all of the points that are not acceptable at the current iteration. This is illustrated in Figure 1, where the forbidden region is highlighted for different values of q.

Merit functions have been criticized for many reasons. First, it is not so easy to choose an initial value for q, since (x, l) and j(x) usually have very different meanings and units. Besides, it is necessary to decrease q as the algorithm progresses to force it to find a feasible solution. If the initial penalty parameter used is near to 1 and q is decreased slowly, the algorithm may take too many iterations to reach a feasible point. On the other hand, starting from a small q or decreasing this factor too quickly may force iterates to stay almost feasible, shortening the steps even when we are far from the optimal solution.

As shown in [10], the adoption of a non-monotone strategy for the reduction of q is very effective to avoid this premature step shortening. However, it also allows the algorithm to cycle between small and large penalty parameters, inducing some zigzagging in many cases. This undesired behavior can be controlled if an efficient choice of the non-monotonicity parameter is adopted, but this choice is also very problem dependent, so the criticism of the method still applies.

1.2 The filter method

To overcome some of the difficulties inherent to merit functions, Fletcher and Leyffer [8] introduced the idea of using a filter. Instead of combining infeasibility and optimality in one single function, the filter method borrows the idea of nondominance from multi-criteria optimization theory. Thus, a trial point is refused only if it is dominated by some other point, generated in a previous iteration.

This approach was promptly followed by many authors, mainly in conjunction with SLP (sequential linear programming), SQP and interior-point type methods (see, for instance, [1, 5, 6, 7, 9, 11, 12, 15, 16, 17, 22, 23, 24, 25]).

The SQP-filter algorithm presented in [6] illustrates how this nondominance criterion works. In this method, the objective function f is used to measure optimality, while infeasibility is measured by

Function f(x) differs from j(x) in the norm used. Besides, the violation of the bound constraints is also taken in account here, since the filter method do not require the iterates to satisfy these constraints.

For each approximate solution x, we can plot a point (f(x), f(x)) in the (f, f)-plane, just as we did with (j(x), (x)) in Figure 1. A third equivalent plane will be introduced in next section, as the algorithm presented in this paper works with (j, f) pairs.

Let be a set of previously generated pairs in the form (f_j, f_j). An iterate x_k is accepted by the filter whenever it satisfies

where g Î (0,1) is a constant.

One advantage of this type of Pareto dominance criterion is that it does no require the redefinition of a penalty parameter at each iteration. However, to avoid the acceptance of iterates that are too close to the points in , it was necessary to add to the last inequality of (6) a term that depends on the infeasibility, giving a merit function flavor to this specific filter.

The main disadvantage of the SQP-filter based on (6) is that this acceptance criterion is too tolerant. In fact, requiring only the infeasibility or the optimality to be improved makes possible the acceptance of points that are only marginally less infeasible but have a large increase in f(x) over the current iterate, or vice-versa.

Even though, the SQP-filter method can give us some good hints on how to improve the algorithms based on merit functions.

The first hint is that the same merit function that is reliable for points in the (j, f)-plane that are near to (j(x_k), f(x_k)) may be not so useful when the step is large, so the trial point is far from the current iterate. As illustrated in Figure 1, for values of q near to 1, the acceptance criterion based on a merit function cuts off a significant portion of the feasible region, including, in many cases, the optimal solution of the problem.

The second good idea behind the SQP-filter method is that a restoration step should be used sometimes. The objective of a restoration is to obtain a point that is less infeasible than the current one and is also acceptable for the filter. In [6, sec. 15.5], a restoration step is computed when the trust region quadratic subproblem is incompatible (i.e. has an empty feasible set). In our method, this strategy will be used whenever staying away from feasibility seems not worthwhile. In other words, infeasible iterates are welcome only if they generate large reductions of the objective function. If the decrease in f is small and the current point is very infeasible, it is better to move off and find a more feasible point.

The last lesson we can take from the SQP-filter method is that feasible points could never be refused. In [6, sec. 15.5], feasible points are never added to , so a feasible iterate will always satisfy the first inequality of (6). This assures that the optimal solution will always be accepted by the algorithm and a restoration will always succeed.

1.3 Motivation and structure of the paper

The objective of this paper is to present an algorithm that takes advantages from both the merit function and the filter ideas. The algorithm introduced here is a merit function SQP method, in the sense that it still combines feasibility and optimality in one single function and it still uses penalty parameters. However, not one but several penalty parameters are defined per iteration, each one related to a portion of the (f, j)-space. These parameters are also automatically computed from some (f, j)-pairs collected at previous iterations, so no update scheme need to be defined for q.

This paper is organized as follows. In the next section, we present the piecewise linear function we use to accept or reject points. Section 3 introduces the proposed algorithm. In section 4, we prove that the algorithm is well defined. Sections 5 and 6 contain the main convergence results. Finally, in section 7 some conclusions are presented, along with lines for future work.

Through the paper, we will omit some (or even all) of the arguments of a function, if this does not lead to confusion. Therefore, sometimes Q(H, x, s) will be expressed as Q(s), for example, if there is no ambiguity on H and x.

2 A merit function that uses several penalty parameters per iteration

As we have seen, a merit function deals with two different concepts: the infeasibility and the optimality of the current point.

In this paper, we will introduce a new merit function that uses information collected at previous iterations to accept new points. Since we compare points generated at different stages of the algorithm, this function cannot be based on the augmented Lagrangian, like in [10], as it would depend on the Lagrange multiplier estimates used and, obviously, these estimates change from one iteration to another. Therefore, we decided to adopt the so called smooth

₂ merit function, defined as:

The actual reduction of the

₂ merit function will be given by

where

Similarly, the predicted reduction of the merit function will be defined asin (3), replacing (5) by

where

Generally, for a trial point to be accepted, it is necessary that the actual reduction of the merit function satisfies

A_red(x, s, q) > hP_red(H, x, s, q),

where h Î (0, 1) is a given parameter.

However, this scheme based on a linear merit function is usually unreliable for trial points that are far from the current iterate. Therefore, we suggest the use of a piecewise linear function to accept or reject new points, which correspond to use several merit functions per iteration.

In order to define the new merit function, let F be a set of p points (j_i, f_i) in the (j, f)-plane. Suppose that these pairs are ordered so that j₁ < j₂ < ¼ < j_p. Suppose also that each point (j_i, f_i) in F is below the line segment joining (j_i-1, f_i-1) and (j_i+1, f_i+1), for i = 2, ¼, p - 1. Thus the piecewise linear function that passes through all of the points in F is convex.

For each point (j_i, f_i) in F, define another point (, ) by moving a little towards the southwest (a precise definition of and is given in (10) and (11) below). Let be the set of points (, ). The convex piecewise linear function that connects the points in is defined by

where g_s is a small positive constant, such as 10^-4.

This new function, illustrated in Figure 2, is formed by p + 1 line segments that can be viewed as merit functions in the form (7). The i-th of these functions is defined by the penalty parameter

and a particular choice of h that will be defined below.

The region between the piecewise linear function that passes through the points in F and the function defined by acts as a margin that prevents the acceptance of iterates that are not ''sufficiently'' better than the pairs in F.

At each iteration k, _k is generated defining, for each point (j_i, f_i) Î F_k, another point (, ) such that,

and

for some 0 < g_f < g_c < 1. Reasonable values for these constants are g_f = 10^-4 and g_c = 10^-3.

Our algorithm starts with F₀ = Æ. At the beginning of an iteration, say k, we define the temporary set _k as

As it will become clear in the next section, depending on the behavior of the algorithm, the pair (f(x_k), j(x_k)) may be permanently added to F_k+1 at the end of the iteration. In this case, each pair (f(x_i), j(x_i)) that is not below the line segment that joins (f(x_i-1), j(x_i-1)) to (f(x_i+1), j(x_i+1)) is removed from F_k+1, to keep the (F, j) function convex.

A new iterate x₊ = x_k + s_c is rejected if f(x_k + s_c) is above the piecewise-linear function (_k, j(x_k + s_c)) or if we predict a good reduction for the merit function, but the real reduction is deceiving (in the sense that (14) occurs).

To express the first of these conditions in the SQP jargon, we say that x₊ is not accepted if

where

The parameter is used in (12) to define the region between F_k and _k. However, a formula for cannot be written explicitly since we do not know in advance which of the terms in (10) and (11) will be used to define and _i. Fortunately, this formula will only be used in the proof of Lemma 15 and, in this case, a simple expression is known.

Some agreement between the model function and the objective function of (1) is an usual second requirement for accepting the step in filter methods. In the SQP-filter method presented in [6], for example, this extra condition is used whenever is sufficiently greater than a small fraction of the infeasibility at x_k. In this case, s_c is only accepted if /> g_g is satisfied. Most filter algorithms, such as those presented in [5, 7, 9, 16, 17, 22, 23] include similar tests.

The combination of infeasibility and optimality into a single merit function allow us to adopt a less stringent condition, rejecting x₊ only if

where

k > 0, g_mÎ (0, 1) and g_g Î (0, 1).

In words, (14) states that when we predict a good reduction for the optimality part of the merit function, the step is only accepted if there exists a penalty parameter Î [g_m, 1] such that A_red/P_red > g_g.

This condition seems to be somewhat inelegant, since one should expect that all of the iterates that do not belong to the forbidden region shown in Figure 2 are to be accepted. However, if we choose a small g_m, say 10^-4, (14) becomes less stringent than (12) in most cases, so it is seldom used. In fact, we generally have ₀ (º 0) < g_m < ₁, so this condition only applies when x_k is the leftmost point in . And even in this case, the forbidden region is only slightly enlarged.

Finishing this section, Figure 3 illustrates how the new merit function reduces the region of acceptable points in comparison with the filter. The comparison was made for problem 7 from the Hock and Schittkowski test collection [18]. In the figure, the current iterate, x_k = (-0.25, 2.27) is shown as a black dot. The points in the (j, f)-space used to define the set were (0.01, -1.5), (2.28, -2.21) and (20, -2.7). The forbidden region for the filter method, as defined in (6), is given in dark gray. The points that are acceptable for the filter but not for the new merit function are shown in light gray. The white region contains the points that are acceptable by both criteria. The contour lines of f(x) = log(1 + ) - x₂ are concave up, while concave down contour lines are related to the constraint C(x) = (1 + )² + - 4. The white concave down narrow region in the bottom of the graph corresponds to points that are almost feasible, i.e., points that satisfy j(x) < 0.01.

3 An SQP algorithm

In the general framework of a trust region sequential quadratic programming algorithm, a step s_c is obtained approximating problem (1), in a neighborhood of an iterate x_k, by a quadratic programming (QP) problem.

In our case, this QP problem has the form

where Q(H, x, s) is defined by (8), x_k is supposed to belong to

W = {x Î n | l < x < u}

and H_k is an approximation of the Hessian of the Lagrangian at x_k. It should be noticed that H_k does not need to be positive definite, so this problem must be carefully handled. The infinity norm was chosen here so the bound constraints and the trust region constraint of (16) can be grouped.

We will use the term j-stationary to say that a point satisfies the first order optimality conditions of

Unfortunately, if x_k is not j-stationary, the constraints of (16) may be inconsistent, so this problem may not have a solution. Some algorithms, such as the SQP-filter method presented in [6], include a restoration step, called to find a new point that makes (16) compatible. Another common practice to overcome this difficulty is to directly divide the step s_c into two components. The first of these components, called normal step, or simply s_n, is obtained as the solution of the feasibility problem

where b_dÎ (0, 1] is a given constant. If M(x_k, s_n) = 0, then x_k can be substituted by x_k + s_n in (16) to make this problem feasible, so it can be solved by any QP algorithm. Otherwise, the second component of s_c, called the tangential step, or s_t, is computed so Q is reduced but the predicted reduction of the infeasibility obtained so far is retained. In other words, s_c is the solution of the (now consistent) problem

Usually, b_d is set to some value around 0.8 so the trust region is enlarged from (18) to (19). This is done to prevent s_n from being the the only solution of (19) when this point is in the border of the trust region of (18).

To insure a sufficient decrease of M, a Cauchy point, , is computed. This Cauchy point is based on a decent direction for j(x) given by P_w(x_k - Ñj(x_k)), the orthogonal projection of x_k - Ñj(x_k) on W. The solution of (18) is required to keep at least ninety percent of the reduction obtained by .

A similar procedure is adopted for (19). In this case, , the Cauchy point, is obtained from a descent direction for f(x) on the tangent space, given by P_x(-ÑQ(s_n)), the orthogonal projection of -ÑQ(s_n) on the set

Again, the decrease on Q obtained by the solution of (19) must not be less than a fixed percentage of the reduction supplied by the Cauchy point.

Besides using this two-step scheme, the algorithm presented here also performs a restoration whenever the step becomes too small and the current point is very infeasible. Although this feasibility reinforcement seems to be unnecessary, the restoration plays a very important role in accelerating the method and keeping the trust region radius sufficiently large.

This role is better explained by an example. Going back to Figure 2, let's suppose that the current iterate x_k is represented by () and that the optimal point has a function value greater than . In this case, unless the trust region radius is sufficiently large and our quadratic model is really good (so the algorithm can jump over a large portion of the forbidden region), it will be necessary to perform a considerable number of short-step iterations to traverse from the southeast to the northwest part of the figure. To shorten the path between the current point and the desired solution, it is necessary to focus only on reducing the infeasibility, and this is exactly what the restoration step does.

One may notice that the use of a standard filter is not useful for circumventing this difficulty. In fact, the problem is aggravated in the presence of (almost) right angles in the frontier of the forbidden region. To see why this happens, let's suppose now that the restoration is used only to ensure that problem (16) is compatible and that x_k is in the vicinity of an almost horizontal segment of the filter envelope. In this case, to escape from this region, it may be necessary to severely reduce the trust region radius before (16) becomes inconsistent, so the restoration is called.

The main steps of the algorithm are given below. We start from k = 0 and take F₀ = Æ as the initial set of points used to define the piecewise linear function (F). An initial point x₀Î W, an initial trust-region radius D₀> D_min and an initial symmetric matrix H₀ need also to be given.

Algorithm 1. A new SQP algorithm

1. WHILE the stopping criteria are not satisfied

1.1.

k F_k {(f(x_k), j(x_k))};

1.2. IF ||C(x_k)|| = 0 (x_k is feasible),

1.2.1. s_n 0;

1.3. ELSE

1.3.1. Compute d_n (a descent direction for j(x)):

d_nP_w(x_k - g_nÑ_j(x_k)) - x_k;

1.3.2. Determine (the decrease step for j(x)), the solution of

1.3.3. Compute s_n (the normal step) such that

l<x_k + s_n < u,

||s_n||_¥ < b_dD_k, and

M(x_k, 0) - M(x_k, s_n) > b_m[M(x_k, 0) - M(x_k, )];

1.4. Compute d_t (a descent direction for f(x) on the tangent space):

d_tP_x(-g_tÑQ(s_n));

1.5. Determine (the decrease step for f(x)), the solution of

1.6. Compute a trial step s_c such that

A(x_k)s_c = A(x_k)s_n,

l<x_k + s_c < u,

||s_c||_¥ < D_k, and

Q(s_n) - Q(s_c) = b_q[Q(s_n) - Q()];

1.7. IF (f(x_k + s_c) > (_k, j(x_k + s_c))) OR

((H_k, x_k, s_c) > kj(x_k) AND < g_m),

1.7.1. D_k a_R min{D_k, ||s_c||_¥}; (reduce D)

1.8. ELSE

1.8.1. r_k

(x_k, s_c)/(H_k, x_k, s_c);

1.8.2. IF (H_k, x_k, s_c) < kj(x_k) OR r_k < g_r,

1.8.2.1. F_k+1

_k; (include (f(x_k), j(x_k)) in F)

1.8.3. ELSE F_k+1

F_k;

1.8.4. Accept the trial point:

Determine H_k+1;

kk + 1;

1.9. IF D_k < D_rest AND j(x_k) > _h,

1.9.1. Compute a restoration step s_r so that

(j(x_k + s_r) < _h AND f(x_k + s_r) < (_k, j(x_k + s_r)))

OR x_k + s_r is j-stationary but infeasible;

1.9.2. F_k+1

_k; (include (f(x_k), j(x_k)) in F)

1.9.3. Accept the new point:

x_k+1

x_k + s_r;

D_k+1 max{b_rD_rest, D_min};

Determine H_k+1;

kk + 1;

The constants used here must satisfy 0 < b_d< 1, 0 < b_m < 1, 0 < b_q < 1, k > 0, 0 < g_r < g_g < h < 1, g_n > 0, g_t > 0, 0 < D_min < D_rest, 0 < a_R < 1, a_A> 1, _h > 0 and b_r > 0. Parameters k, g_n, g_t, D_min, D_rest, _h and b_r are problem dependent and must be chosen according with some measure of problem data. Reasonable values for the remaining parameters might be b_d = 0.8, b_m = 0.9, b_q = 0.9, g_r = 0.01, g_g = 0.05, h = 0.9, a_R = 0.25 and a_A = 2.5. The constant h should not be confused with the parameter defined in (12).

Algorithm 1 seems to be somewhat inefficient, since it requires the solution of several quadratic programming problems per iteration. However, as it will become clear at section 7, only steps 1.3.3 and 1.6 usually require the solution of a quadratic problem. Besides, we do not need to solve these problems exactly, so the algorithm is competitive with modern nonlinear programming codes.

The restoration step is called whenever the trust region radius becomes too small compared to j(x_k). In this case, we need to find a point that is less infeasible than x_k and that is also acceptable for the piecewise linear merit function. To ensure that such a point can always be obtained, it would be necessary to use an algorithm for global minimization. Of course, this alternative is unaffordable, so the restoration is computed in practice by an algorithm for solving the box constrained problem (17). The only drawback of this approach is that we cannot guarantee that a stationary but infeasible point would not be reached. Therefore, we say that the main algorithm fails if this happens.

If x_k is feasible, then the condition < kj(x_k) is never satisfied, since is always greater or equal to zero. Besides, the condition < g_r is also never satisfied when x_k is feasible and f(x_k + s_c) < (_k, j(x_k + s_c)). Therefore, all of the points in F_k are infeasible, although _k may contain a feasible point. This result is very important for two reasons. First, it prevents the optimal solution of problem 1 from being refused by the algorithm. Moreover, it also assures that the algorithm is well defined, as stated in the next section.

4 The algorithm is well defined

An iteration of algorithm 1 ends only when a new point x_k + s is below the piecewise linear function (_k, j(x_k + s)), besides satisfying some other conditions stated at steps 1.7 or 1.9.1. While such a point is not found, the trust region radius is reduced and the iteration is repeated. It is not obvious that an acceptable point will be obtained, as we may generate a sequence of points that are always rejected by the algorithm. In this section, we prove that the algorithm is well defined, i.e. a new iterate x_k+1 can always be obtained unless the algorithm stops by finding a j-stationary but infeasible point or a feasible but not regular point.¹

In the following lemma, we consider the case where x_k is infeasible.

Lemma 1. If x_k is not j-stationary, then after a finite number of repetitions of steps 1.1 to 1.9, a new iterate x_k+1 is obtained by the algorithm.

Proof. At each iteration k, if f(x_k + s_c) < (_k, j(x_k + s_c)) and one of the conditions

(for some

> g_m) is satisfied, then x_k + s_c is accepted and we move to iteration k + 1. Otherwise, D_k is reduced and after some unfruitful steps, D_k < D_rest and j(x_k) > _h, so a restoration is called.

Suppose that a j-stationary but infeasible point is never reached (otherwise the algorithm fails). As the restoration generates a sequence of steps {s_j} converging to feasibility, and since _k does not include feasible points (because x_k is infeasible and no feasible point is included in F_k), there must exist an iterate x_k + s_r that satisfies j(x_k + s_r) < min{1, _h}, so we can proceed to the next iteration.

Now, in order to prove that the algorithm is also well defined when x_k is feasible, we need to make the following assumptions.

A1. f(x) and C_i(x) are twice-continuously differentiable functions of x.

A2. The sequence of Hessian approximations {H_k} is bounded.

As a consequence of A1 and A2, the difference between the actual and the predicted reduction of the merit function is proportional to D², so the step is accepted for a sufficiently small trust region radius, as stated in the following lemma.

Lemma 2. Suppose that A1 and A2 hold and that x_k is feasible and regular for problem 1 but the KKT conditions do not hold. Then, after a finite number of trust region reductions, the algorithm finds a new point x_k + s_c that satisfies f(x_k+s_c) < (_k, j(x_k+s_c)) and A_red(x_k, s_c, ) > g_g P_red(H_k, x_k, s_c, ) for some

> g_m.

Proof. Since x_k is feasible, s_n = 0. Supposing that x_k is regular and non-stationary, there must exist a vector d_t ¹ 0 satisfying

Let us define, for all D > 0,

p(D) = t(D)d_t,

where

Clearly, x + d_t Î W, so we have that ||t(D)d_t||_¥ = D whenever D < ||d_t||_¥. Define, in this case,

Since Q() < Q(p(D)), by elementary properties of one-dimensional quadratics, there exists D₁Î (0,||d_t||_¥] such that, for all D Î (0, D₁),

Moreover, since x_k is feasible and As_n = 0, we have that M(x_k, 0) = M(x_k, s_c) = 0, so

Once x_k is feasible, (j(x_k), f(x_k)) is the first pair in _k. Thus, there exists ₂Î (0, D₁] such that, for D < ₂, we need to consider only the portion of (_k, j) defined on the interval [0, j₂]. This linear function may be rewritten so the condition

f(x_k + s_c) < (_k, j(x_k + s_c))

is equivalent to

where,

and

₁ > 0 is given by (9).

Now, by A1, A2 and the definition of P_red, we have

So, using (23) and (24) we deduce that

Thus, for D < min{ (1 - h₁) b_qc

₁/|₁|, ₂} = ₃, the inequality (22) necessarily takes place.

Now, using the fact that = 1 for x_k feasible and replacing ₁ by 1 in (25), we can conclude that, for

the condition A_red(x_k, s_c(D), ) > g_gP_red(H_k, x_k, s_c(D), ) is also satisfied and the step is accepted.

5 The algorithm converges to a feasible point

As mentioned in the last section, our algorithm can stop if a j-stationary but infeasible point is found. Naturally, this unexpected behavior of the algorithm makes somewhat pretentious the title of this section.

Formally, what we will prove is that, supposing that a j-stationary but infeasible point is never reached and that the restoration always succeeds, an infinite sequence of iterates converges to feasibility.

In the proofs of the lemmas presented here, we will suppose that A1 and the following assumption are satisfied.

A3. The sequence of iterates {x_k} lies within a closed and bounded domain W₀.

This requirement is easy to fulfill, as we usually can define finite lower and upper limits for the variables. Besides, as mentioned in [6, p. 730], assumptions A1 and A3 together ensure that, for all k,

for some constants f^min, f^max and j^max > 0. Our analysis will be based on the fact that the rectangle A₀ = [0, j^max] × [f^min, f^max] is covered by a finite number of rectangles with area greater than a small constant. Therefore, each time we expand the forbidden region (see fig. 2) by adding to it a small rectangle, we drive the iterates towards feasibility.

Let us start investigating what happens to j(x) when an infinite sequence of iterates is added to F, following the skeleton of Lemma 15.5.2 of [6].

Lemma 3. Suppose that A1 and A3 hold and that {k_i} is any infinite subsequence at which the iterate

is added to F. Then

Proof. Let us suppose, for the purpose of obtaining a contradiction, that there exists an infinite subsequence {k_j} Í {k_i} for which

where > 0.

At iteration k_j, the (j, f)-pair associate with is included in F at position m, which means that j_m-1< (º j_m) < j_m+1 and f_m-1> (º f_m) > f_m+1. Thus, as long as the pair (, ) remains in F, no other (j, f)-pair is accepted within the rectangle

Notice that, by (10), (11) and (27), the area of this rectangle is, at least,

Assume now that (, ) is excluded from F by another pair (, ), included in F at an iteration k_l > k_j. This case is illustrated in Fig. 4. Notice that (, ) cannot fall in regions I and V since, in this case, (, ) will not be excluded from F. It can be easily verified that the worst case occurs when (, ) lies on p₁(j) or p₂(j).

Suppose (, ) lies on p₂(j), as depicted in Fig. 4. In this case, the rectangle r_m will be entirely above ₂, the line that connects (, ) to (_m+1, _m+1). Since ₂ will be included in the new piecewise linear function (), no point within r_m can ever be reached by a new iterate.

The same idea can be applied in the case (, ) lies on p₁(j). Therefore, once (, ) is included in F, r_m will always be above (). Since the area of this rectangle is at least and the set A₀ is completely covered by at most Surf(A₀)/ of such rectangles, it is impossible for an infinite subsequence of {k_i} to satisfy (27), and the conclusion follows.

Finally, we are going to consider the case where no point is added to F_k for k sufficiently large.

Lemma 4. Suppose that assumptions A1 and A3 hold. Suppose also that, for all k > k₀, x_k is never included in F_k. Then,

Proof. Since x_k is not included in F_k, no restorations are made and both conditions stated at step 1.8.2 of algorithm 1 are never satisfied for k > k₀. Therefore, we have

for all k > k₀, which means that the objective function always decrease between infeasible iterations. Since A1 and A3 imply f^min< f(x_k) < f^max, we must have

Then, (28) follows from (29) and (30).

6 The algorithm finds a critical point

Finally, we are able to prove the convergence of the algorithm to a stationary point for (1). In order to do that, we will need to make one additional assumption on the choice of the normal step s_n.

A4. There exist k_n, k_N > 0 such that, if {} is a subsequence of iterates that converges to a feasible point, the choice of s_n at step 1.3.3 of algorithm 1 satisfies

||s_n(x_k, D_k)|| < k_n||C(x_k)||₂

for k_i > k_N.

This requirement is also easy to fulfill since it only applies when the infeasibility is small and, in this case, it is reasonable to suppose that the normal step will also be small.

In the following lemma, derived from Lemma 6.1 of [10], we show that in the neighborhood of a feasible, regular and non-stationary point, the directional derivative of the quadratic model (8) along d_t is bounded away from zero.

Lemma 5. Suppose that A2 and A4 hold and that {} is an infinite subsequence that converges to the feasible and regular point x^*Î W, which is not stationary for (1). Then there exists k₁, c₁ > 0 such that

for all x Î { |k > k₁}. Moreover, ||d_t(H, x, D)|| is bounded and bounded away from zero for all x Î { | k > k₁}.

Proof. For all x Î {}, we have that

d_t(H,x,D) = P_x(-g_tÑQ(s_n(x, D))) = P_x(-g_t [Hs_n(x, D) + Ñf(x)]).

By the contractive property of the orthogonal projections,

||P_x(-g_t[Hs_n(x, D) + Ñf(x)]) - P_x(-g_tÑf(x))||₂< g_t ||H||₂||s_n(x, D)||₂.

So, by A2 and A4, we have that

for k > k_N. By the continuity of Ñf(x) and the fact that {} converges, we deduce that

Notice that P_x(-g_tÑf()) is the solution of

Now, define (-g_t Ñf(x^*)) as the solution of

Since x^* is regular but is not a stationary point for (1), it follows that z = 0 is not a solution for (34). So, (-g_t Ñf(x^*)) ¹ 0. Moreover, since z = 0 is feasible for (34), we have that

which implies that Ñf(x^*)^T(-g_tÑf(x^*)) < 0.

Using the fact that P_x(-g_tÑf(x)) is a continuous function of x and s_n for all regular x (see [10]), we can define ₂, ₃, ₄ > 0 and ₂Î such that, for all x Î { | k > ₂}, we have

Now, from (32), (33) and (35), the continuity of C(x) and the feasibility of x^*, there exists ₃> max{₂, k_N} such that, whenever x Î { | k > ₃},

Therefore, ||d_t(H, x, D)|| is bounded and bounded away from zero for all x Î { | k > ₃}.

Finally, since d_t Î (A(x)), assumptions A2 and A4 hold, and ||d_t|| is bounded, we have that, for all x Î { | k > ₃},

where g₃ > 0. Then, (31) follows defining c₁ = 4/4 and choosing k₁ > ₃ such that ||C(x)|| < 4/(4g₃).

Using Lemma 5, we state in the next lemma that, in the neighborhood of a feasible, regular and non-stationary point, the decrease of the quadratic model (8) is proportional to the trust region radius D.

Lemma 6. Suppose that A2 and A4 hold and that {} is an infinite subsequence that converges to the feasible and regular point x^*Î W, which is not stationary for (1). Then there exists c₂, k₂ > 0 and D₁Î (0, D_min) such that

for all x Î { | k > k₂}.

Proof. See Lemma 6.2 of [10].

Now, we are able to present a crucial lemma, derived from Lemma 6.3 of [10], that relates to the trust region radius in the neighborhood of a feasible point.

Lemma 7. Suppose that A1, A2 and A4 hold and that {} is an infinite subsequence that converges to the feasible and regular point x^*Î W, which is not stationary for (1). Then there exists , c₃, k₃ > 0 and D₁Î (0, D_min) such that, for k_i> k₃, if

we have that

Proof. By Lemma 6, assumptions A1 and A4 and the convergence of {}, we have that

Q(0) - Q(s_c) >Q(s_n) - Q(s_c) - |Q(0) - Q(s_n)| >c₂ min {D, D₁} - g₄ ||C(x)||

for all x Î { | k > k₂}, where c₂, k₂ and D₁ are defined as in Lemma 6 and g₄ > 0. Therefore, (38) follows if we choose c₃ < c₂ and k₃> k₂ such that < a² (c₂ - c₃)²/(), where a = min{1, D₁/ D}.

We next examine what happens if D is bounded away from zero and an infinite subsequence of points is added to F.

Lemma 8. Suppose that A1, A2, A3 and A4 hold and that {} is an infinite subsequence at which

is added to F. Suppose furthermore that the limit points of this sequence are feasible and regular, that the restoration always terminates successfully and that

> D₂, where D₂ is a positive scalar. Then there exists a limit point of this sequence that is a feasible, regular and stationary point for (1).

Proof. From assumption A3, we know that {} has a convergent subsequence, say {}. Let us suppose that the limit point of this subsequence is not stationary for (1).

From Lemma 3 we know that there exists

₅Î such that, for k_i > ₅,

Thus, a restoration is never called for k_i > ₅. So, the hypothesis that is added to implies that A_red(, s_c, ) > g_gP_red(, s_c, ) for some > g_m, and that one of the inequalities stated at step 1.8.2 of the algorithm is satisfied at iteration k_i.

Suppose, for the purpose of obtaining a contradiction, that {} converges to a point that is not stationary for (1). So, from Lemma 3 and (38), there exists ₆> ₅ such that j() < and

for all k_i > ₆.

Using Lemma 3 again, we can deduce that there exists

₇> ₆ such that j() < (c₃/k) min{D₁, D₂} and the condition < kj(x_k) is never satisfied for k_i > ₇.

Therefore, r_k < g_r must hold. To show that this is not possible, let us write the inequality A_red(, s_c, ) > g_gP_red(, s_c, ), as

Using the hypothesis that r_k < g_r and the fact that (, s_c) > 0, we have

Then, taking k₄ > k₃ (defined in Lemma 7), we deduce from (38) that, for k_i > k₄,

(1 - ) (j() - j(+ s_c)) > (g_g - g_r) c₃ min {D₁, D₂}.

But, since, g_g > g_r and lim_{i ® ¥} = 0, we must have

which contradicts the fact that

> g_m. Therefore, {} must converge to a stationary point for (1).

Supposing again that D is bounded away from zero, we will now complete our analysis investigating what happens when no iterates are added to F for k sufficiently large.

Lemma 9. Suppose that A1, A2, A3 and A4 hold, that x_k is always accepted but F_k remains unchanged for k > k₅ and that D_k> D₃, for some positive D₃. Suppose also that the limit points of the infinite sequence {x_k} are feasible and regular. Then there exists a limit point of {x_k} that is a feasible, regular and stationary point of (1).

Proof. Assumption A3 implies that there exists a convergent subsequence {}. If the limit point of this subsequence is not stationary for (1), then from Lemma 7, we have

for all k_i > max{k₅, k₃}. Moreover, since is always accepted and F_k is not changed, we deduce, from step 1.8.2 of the algorithm, that r_k> g_r. Therefore,

From (39) and (40) we conclude that f() - f(+ s_c) > g_rc₃ min{D₁, D₃} for all k_i sufficiently large, which contradicts the compactness assumption A3. Thus, the limit point of {} must be stationary.

In the last part of this section, we will discuss the behavior of the algorithm when D ® 0. We will start showing that the predicted reduction of the quadratic model is sufficiently large when D is small.

Lemma 10. Suppose that A2 and A4 hold and that {} is an infinite subsequence that converges to the feasible and regular point x^*Î W, which is not stationary for (1). Suppose also that j_k satisfies (37) and that

for k_i> k₆, where c₃, and D₁ are defined as in Lemma 7. Then > kj().

Proof. Suppose, for the purpose of obtaining a contradiction, that < kj(x_k) for some k_i > k₆. Then, from (38), we have

c₃ min {D, D₁} << kj() < k D²,

which is impossible because of (41). Thus > kj() must hold.

The purpose of the next four Lemmas is to prove that there exists a sufficiently small trust region radius so the step is always accepted and D is not reduced further at step 1.7.1 of algorithm 1.

The first lemma shows the relation between the predicted reduction of the infeasibility and D.

Lemma 11. Suppose that assumption A1 holds and that x_k is not j-stationary. Then there exists D₆, c₄ > 0 such that

if D_kÎ (0, D₆).

Proof. The proof of this Lemma is based on the same arguments used to obtain (21). However, here we deal with the reduction of the infeasibility, instead of the reduction of the objective function, so some modifications need to be done.

Firstly, we must notice that d_n ¹ 0, since we suppose that x_k is not j-stationary. Thus, we can redefine (20) as

t(D) = max {t > 0 | [x_k, x_k + td_n] Î W and ||td_n|| < b_d D}.

Now, using the fact that for D sufficiently small we have ||t(D)d_n|| = b_d D and defining

there must exist D₆Î (0, ||d_n||] such that

for all D Î (0, D₆). Therefore, for the normal step s_n computed at step 1.3.3 of algorithm 1, we have

M(0) - M(s_n) > b_db_m cD_k.

But, since A(x_k)s_c = A(x_k)s_n, we deduce from (4) that

and the desired inequality follows.

In order to prove that x_k + s_c will be accepted, we need to consider how and are computed. Let us begin using the previous lemma to show that, for a small D, , defined in (10), will depend on the predicted reduction of the infeasibility.

Lemma 12. Suppose that A1 holds and that x_k is not j-stationary. Then there exists D₇ > 0 such that

g_c(x_k, s_c) > g_fj(x_k),

if j(x_k) < _h and D_k Î (0, D₇).

Proof. Lemma 11 ensures that

Defining D₇ = min{ g_c c₄ / (g_f

_h), D₆}, where D₆ is given in Lemma 11, we have that

for all D Î (0, D₇), so the desired result follows.

Using Lemma 11 again, we can also show that , defined in (11), will depend on if D is sufficiently small.

Lemma 13. Suppose that A1, A2 and A4 hold, that {} is an infinite subsequence that converges to the feasible and regular point x^*Î W, which is not stationary for (1), and that

_k is given by (10). Then there exists D₈ > 0 such that

g_f(H_k, x_k, s_c) > (j(x_k) - _k),

if j(x_k) < min{_h, } and D_kÎ (0, D₈), where

is defined as in Lemma 7.

Proof. From Lemma 7 we deduce that, if D_k Î (0, D₁], then

g_f (H_k, x_k, s_c) > g_f c₃D_k.

Now, defining D₈ = min{g_f c₃ / _h, D₁}, we have

g_f (H_k, x_k, s_c) > _hD₈D > _h

> j(x_k) > (j(x_k) - _k)

and the desired conclusion follows.

Let us prove now that, for a infeasible x, the ratio between the actual and the predicted reduction of the merit function is sufficiently large if D is small, not matter how the penalty parameter q is chosen.

Lemma 14. Suppose that A1, A2 and A4 hold, that q Î [0, 1] and g Î (0, 1) are given parameters and that {} is an infinite subsequence that converges to the feasible and regular point x^*Î W, which is not stationary for (1). Then there exists D₉ > 0 such that, for k_i sufficiently large,

A_red(x_ki, s_c, q) > gP_red(, , s_c, q)

for all Î (0, D₉), if

is infeasible.

Proof. Since lim_{i ® ¥} j() = 0 and, at the beginning of iteration k, the trust region radius satisfies > D_min, there must exist k₇> k₃ (defined in Lemma 7) such that, for k_i > k₇, the condition j() < is satisfied, so (38) holds. Besides, since is infeasible, (42) also holds if we take D₉< D₆ (defined in Lemma 11). Therefore, from the definition of P_red, we have that

P_red(, , s_c, q) > qc₃ min {, D₁} + (1 - q) c₄.

where D₁, c₃ and c₄ are defined as in Lemmas 7 and Lemma 11.

For any q Î [0, 1], the above inequality implies that

P_red(, , s_c, q) > min {, D₁} min{c₃, c₄} > 0.

But, from A1 and A4, we also have that

|A_red() - P_red()| < c₅.

for some c₅ > 0. From the last two inequalities, we deduce that, for < D₁,

Therefore, defining D₉ = min{(1 - g) min{c₃, c₄}/c₅, D₁, D₆}, we obtain the required result.

In our last lemma, we will use the previous results to prove that, if D ® 0, there is no infinite subsequence that converges to a point that is not stationary for (1).

Lemma 15. Suppose that A1, A2, A3 and A4 hold. Suppose also that the limit points of the infinite sequence {x_k} are feasible and regular and that lim_{k® ¥} D_k = 0. Then there exists a limit point of {x_k} that is a stationary point of (1).

Proof. Assumption A3 implies that there exists a convergent subsequence {}. Let us suppose, for the purpose of obtaining a contradiction, that the limit point of this subsequence is not stationary for (1).

Now, we need to consider separately two mutually exclusive situations. First, let us suppose that is feasible. In this case, Lemma 2 assures that, for D < ₄ (defined in (26)), the step is accepted and the trust region radius need not to be reduced further.

On the other hand, if is not j-stationary, Lemmas 12 and 13 assure that, for < min{D₇, D₈} = D₁₀, we have > g_fj_i /g_c and > (j_i - / g_f, so the parameter used in (12) is given by

where i is defined in such a manner that _i-1< j(x₊) < . Thus, 0 < < 1 and f( + s_c) > (, j( + s_c)) is equivalent to (12).

Now, using Lemmas 10 and 14, we can deduce that, for < D₉, the step is always accepted. Consequently, > a_R min{₄, D₉, D₁₀}, which contradicts the hypothesis that lim_{k ® ¥} D_k = 0, so we conclude that the limit point of the subsequence {} is a stationary point of (1).

Finally, let us state a theorem that puts together all of the results presented so far.

Theorem 1. Suppose that A1, A2, A3 and A4 hold and that {x_k} is an infinite sequence generated by Algorithm 1. Then either the restoration converges to a j-stationary but infeasible point of (1), or lim_{k ® ¥} j(x_k) = 0. Moreover, if the restoration always succeeds and all of the limit points of {x_k} are regular, there exists a limit point x^* that is a regular and stationary point for (1). In particular, if all of the j-stationary points are feasible and regular, then there exists a subsequence of {x_k} that converges to a feasible, regular and stationary point of (1).

Proof. This result is a direct consequence of Lemmas 3, 4, 8, 9 and 15.

7 Numerical experience

The description of some steps of Algorithm 1 was intentionally left vague to suggest the reader that several alternative implementations are available. In this section, we describe one of such possible implementations and present the numerical results obtained by applying the algorithm to some problems from the CUTEr collection [13].

7.1 Algorithmic details

The computational effort of algorithm 1 may be decomposed into three main parts. The first is related to the reduction of the infeasibility and includes steps 1.3.1 to 1.3.3. The aim of the second part, that comprises steps 1.4 to 1.6, is to improve the objective function. Finally, a restoration is called, at step 1.9.1, if the infeasibility needs to be drastically reduced. Each one of these parts is briefly described below.

Taking a closer look at steps 1.3.1 and 1.3.2, one may notice that vectors d_n and can be easily determined, since the first involves computing a projection on a box and the second requires only the solution of a one-dimensional problem.

The normal step, s_n, can be obtained by solving the bound-constrained least squares problem (18), replacing the word reduce by minimize. In our experiments, the Quacon package [2] was used for this purpose. The computation of s_n is declared successful when both conditions

M(x_k ,0) - M(x_k, s_n) > b_m[M(x_k, 0) - M(x_k, )]

and

||g_p(s_n)||₂< 0.001||g_p()||₂

are satisfied, where b_m = 0.9 and g_p is the projected gradient of the quadratic function minimized. Otherwise, the limit of max{1000, 6n} iterations is reached.

Vector s_c is computed by applying the MINOS 5.4 [19] solver to the quadratic problem (19), stopping the algorithm when both

and condition (36) are satisfied, where, again, g_p is the projected gradient and l is the vector of Lagrange multipliers. A limit of max{1000, 3n} iterations also applies.

Vector d_t is obtained, at step 1.4, using the MINOS package to approximately solve the problem

where g_t = 0.001. Since d_t is much more difficult to obtain than d_n, we decided to compute this vector only if s_c fails to satisfy conditions (36) and (44). In this case, after determining d_n and , s_c is computed again.

To determine the restored point s_r, successive normal steps were computed until the conditions stated at step 1.9.1 were satisfied. A more sophisticated restoration procedure could be devised, but this very simple scheme seemed to be satisfactory.

In all of the experiments, exact Hessians were computed, using the Lagrange multipliers supplied by MINOS as the approximate dual variables for problem (1).

Algorithm 1 terminates when ||C(x_k)||_¥< 10^-6 and one of the conditions

1. ||Ñ(x_k, l_k)||_¥ < _g max{1, ||Ñ(x_typ, l_typ)||_¥},

2. ||s_c||_¥< _x max{1, ||x_k||_¥}, or

3.

< _s and < _o (for three successive iterations),

is satisfied, where

_g = 10^-6, _x = 10^-18, _s = _o = 10^-8, and x_typ and l_typ are typical primal and dual steps, determined at the first iterations of the algorithm.

The algorithm also stops if, after a restoration step, Ñj(x_k + s) < 10^-10. In this case, it is more than likely that a j-stationary but infeasible point was found.

The remaining parameters used are D₀ = 10⁵, D_min = 10^-5, g_n = 1.0, g_t = 10^-8, g_r = 10^-2, g_g = 5.10^-2, k = 10^-4, h = 0.9, a_A = 2.5, a_R = 0.25, b_d = 0.8, b_q = 0.9, b_r = 16.0, g_s = 10^-4, g_f = 10^-6, g_c = 10^-3, g_m = 10^-4, D_rest = 10^-1 and _h = 10².

7.2 Algorithm performance

To analyze the behavior of the algorithm just described, a set of 82 small to medium-size problems was extracted from the CUTEr collection [13]. The selected problems are presented in Tables 1, 2 and 3. The number of variables of the original CUTEr problem is given by n, while s is the number of slack variables required to convert the problem to the form shown in (1). The number of constraints is given by m.

The algorithm, hereafter called GMM99, was implemented in FORTRAN 77 and the executable program was generated using the ifort 8.1 compiler, under the Fedora 3 Linux operating system. To evaluate its performance, the algorithm was compared to Lancelot (release B), the well known nonlinear programming package distributed along with the Galahad library [14]. The tests were performed on a Dell Optiplex GX280 computer, using an Intel Pentium 4, 540 processor, with a clock speed of 3.2GHz, 1MB of cache memory, a 800MHz front side bus and the Intel 915G chipset. The Lancelot default parameters were used, except for the maximum number of iterations that was increased to 10000. Exact first and second derivatives were computed by both methods.

The new algorithm failed to obtain a feasible point for problems LUBRIF, ROTDISC and YORKNET. Besides, it also could not find an objective function as good as Lancelot for problems BRAINPC5, CATENARY, DIXCHLNG, DTOC2, EG3, ELATTAR, NGONE and OPTCTRL3. These 11 occurrences were classified as GMM99 failures.

Lancelot attained the limit of 10000 iterations without finding a feasiblepoint for problems CRESC132, DITTERT, LAUNCH, NET3, OPTMASS andSMMPSF. For problems DRUGDIS, GAUSSELM, HALDMADS, HVYCRASH, READING3 and SINROSNB, the new algorithm found a better objective function. Finally, for problems BRIDGEND and ROTDISC Lancelot exceeded the limit of 3000 seconds of running time. These 14 cases were treated as Lancelot failures.

Table 4 shows a summary of the results obtained comparing the new algorithm to Lancelot. According to this table, GMM99 outperformed Lancelot in 47 problems (out of 82), while Lancelot achieved the best performance in 32 cases (including GMM99 failures). As mentioned above, both algorithms failed to solve problem ROTDISC.

Table 5 shows the number of times GMM99 and Lancelot obtained the best solution for each type of objective function and constraint. In this table, the problems are also divided according to their size.

For this specific selection of problems, GMM99 performed better for problems with a linear objective function, while Lancelot was the best option for solving problems with a quadratic function and nonlinear constraints. These results are probably related to the fact that MINOS was used to compute the tangential step in GMM99. In general, GMM99 had a slightly better performance.

The new algorithm was also compared with GMM, the standard trust-region SQP method presented in [10]. The results are shown in Tables 6 and 7.

The GMM algorithm was not able to find a feasible solution for problemsCATENARY, CRESC132, LAKES, LUBRIF, ROTDISC and STEENBRC. Moreover, the objective function at the solution was worse than the value obtained by GMM99 for problems AUG2D, COSHFUN, EIGMINB, HAGER1, HAGER2, HIMMELBI, NGONE and SSNLBEAM. For the OPTCTRL3 problem, the algorithm found a solution too far from the one supplied by Lancelot. These problems were classified as GMM failures.

The 9 GMM99 failures cited in Table 6 include problems LUBRIF, ROTDISC and YORKNET, for which a feasible point was not found, problem OPTCTRL3, because Lancelot obtained a much better solution, and also problems BRAINPC5, ELATTAR, GAUSSELM, HVYCRASH and SINROSNB, forwhich GMM has found a lower objective function.

Naturally, the definition of failure used here is arguable, since the fact that an algorithm has found a better solution does not imply that the other method failed to find a local optimum. Besides, problems CATENARY, DIXCHLNG, DTOC2, EG3 and NGONE, classified as GMM99 failures when comparing this algorithm to Lancelot, are not considered as failures in Table 6, since GMM was not able to find a better solution in all of these cases. Even though, we decided to use the term failure to stress that an algorithm that finds a lower solution performs better than the rival, no matter the time spent to find the solution.

Finally, Table 7 shows how the good performance of GMM and GMM99 is distributed among different types of problems.

Tables 6 and 7 clearly show the superiority of GMM99 over GMM. The new algorithm is more robust and attained the best performance for all of the problem classes. Taking into consideration that GMM also uses Minos to compute the horizontal step and Quacon to compute the vertical step, so the core of the algorithm is the same of GMM99, these results suggest that the use of the new acceptance criterion, combined with a restoration procedure, can improve the performance of the SQP algorithm.

For some problems, no points were included in F, so the behavior of GMM and GMM99 were very similar. Considering the whole set of 82 problems, the mean number of points in F at the last iteration of GMM99 was 5. However, for 7 problems, the cardinality of the F set was greater or equal to 20.

The restoration was called at least once for 24 problems. The purpose of the restoration is to avoid the excessive reduction of the step size when x_k is very infeasible. This strategy has shown to be effective for 18 problems, while the performance was not significantly altered for two problems and slightly worsened for 4 problems.

8 Conclusions

In this paper, we depict the general framework of a new SQP algorithm that combines ideas from both merit functions and the filter introduced by Fletcher and Leyffer in [8]. The use of several penalty parameters defined automaticallyby the previous iterates avoids the premature reduction of q as well as the zigzagging that can occur when a non-monotone strategy is used to update this parameter. The new method is also less tolerant than the filter method, since we do not accept points that marginally reduce the infeasibility or the objective function.

The algorithm presented here is unique in the sense that it combines

• The use of a piecewise-linear function to accept or reject iterates. This function can be interpreted as a merit function that selects the penalty parameter according to the position of the predicted point in the (j, f)-plane, but also as a convex piecewise filter.

• The effective use of the points in F to define the forbidden region in the (j, f)-plane. In addition to dominated points, filter-SQP methods cut off the entire northwestern part of the allowed region when the predicted reduction of the objective function is positive. The GMM99 algorithm reject points if they belong to the forbidden region depicted in Figure 2 or if condition (14) is satisfied. However, this second condition is so mild that it is seldom used in practice.

• The use of the restoration to prevent the step shortening inherent to methods that use points in the (j, f)-plane to define the forbidden region.

The results shown in the previous section suggest that the algorithm, as implemented, is superior to the method presented in [10] and slightly better than Lancelot for small to medium-size problems. As a future work, we intend to test other strategies to compute the restoration, as well as other solvers to compute the normal and tangential steps, in order to deal with larger problems. Perhaps, the library used to solve the quadratic subproblems can be selected dynamically, after knowing the size of the problem and the type of the objective function and constraints.

Acknowledgments. The author would like to thank Prof. José Mario Martínez for his insightful comments and Tadao Ando Jr. for the revision of the text.

Received: 05/III/07. Accepted: 12/IV/07.

#716/07.

• [1] C. Audet and J.E. Dennis, A pattern search filter method for nonlinear programming without derivatives. SIAM J. Optimization, 14 (2004), 980-1010.
• [2] R.H. Bielschowsky, A. Friedlander, F.A.M. Gomes, J.M. Martínez and M. Raydan, An adaptative algorithm for bound constrained quadratic minimization. Investigación Operativa, 7 (1997), 67-102.
• [3] R.H. Byrd, J.C. Gilbert and J. Nocedal, A trust region method based on interior point techniques for nonlinear programming. Math. Programming, 89 (2000), 149-185.
• [4] R.H. Byrd, M.E. Hribar and J. Nocedal, An interior point algorithm for large-scale nonlinear programming. SIAM J. Optimization, 9 (2000), 877-900.
• [5] C.M. Chin and R. Fletcher, On the global convergence of an SLP-filter algorithm that takes EQP steps. Math. Programming, 96 (2003), 161-177.
• [6] A. Conn, N.I.M. Gould and P.L. Toint, Trust Region Methods. SIAM, Philadelphia (2000).
• [7] R. Fletcher, N.I.M. Gould, S. Leyffer, P.L. Toint and A. Wächter, Global convergence of a trust-region SQP-filter algorithm for general nonlinear programming. SIAM J. Optimization, 13 (2002), 635-659.
• [8] R. Fletcher and S. Leyffer, Nonlinear programming without a penalty function. Math. Programming, 91 (2002), 239-269.
• [9] R. Fletcher, S. Leyffer and P.L. Toint, On the global convergence of a filter-SQP algorithm. SIAM J. Optimization, 13 (2002), 44-59.
• [10] F.A.M. Gomes, M.C. Maciel and J.M. Martínez, Nonlinear programming algorithms using trust regions and augmented Lagrangians with nonmonotone penalty parameters. Math. Programming, 84 (1999), 161-200.
• [11] C. Gonzaga, E. Karas and M. Vanti, A globally convergent filter method for nonlinear programming. SIAM J. Optimization, 14 (2003), 646-669.
• [12] N.I.M. Gould, S. Leyffer and P.L. Toint, A multidimensional filter algorithm for nonlinear equations and nonlinear least-squares. SIAM J. Optimization, 15 (2005), 17-38.
• [13] N.I.M. Gould, D. Orban and P.L. Toint, CUTEr (and SifDec), a constrained and unconstrained testing environment, revisited. ACM Trans. Math. Software, 29 (2003), 373-394.
• [14] N.I.M. Gould, D. Orban and P.L. Toint, GALAHAD, a library of thread-safe Fortran 90 packages for large-scale nonlinear optimization. ACM Trans. Math. Software, 29 (2003), 353-372.
• [15] N.I.M. Gould, C. Sainvitu and P.L. Toint, A filter-trust-region method for unconstrained optimization. SIAM J. Optimization, 16 (2006), 341-357.
• [16] N.I.M. Gould and P.L. Toint, Global convergence of a hybrid trust-region SQP-filter algorithm for general nonlinear programming. in System Modeling and Optimization XX, E. Sachs and R. Tichatschke (eds.), Kluwer, Dordrecht (2003), pp. 23-54.
• [17] N.I.M. Gould and P.L. Toint, Global convergence of a non-monotone trust-region filter algorithm for non-linear programming. Tech. report RAL-TR-2003-003, CLRC Rutheford Appleton Laboratory, Oxfordshire, UK (2003).
• [18] W. Hock and K. Schittkowski, Test examples for nonlinear programming codes. Lecture Notes in Economics and Mathematical Systems, 187 (1981), Springer, Heidelberg.
• [19] R.B. Murtagh and M.A. Saunders, Large-scale linearly constrained optimization. Math. Programming, 14 (1978), 41-72.
• [20] J. Nocedal and S.J. Wright, Numerical Optimization. Springer, New York (1999).
• [21] T.D. Plantenga, A trust-region method for nonlinear programming based on primal interior point techniques. SIAM J. Sci. Computing, 20 (1999), 282-305.
• [22] M. Ulbrich, S. Ulbrich and L.N. Vicente, A globally convergent primal-dual interior point filter method for nonconvex nonlinear programming. Math. Programming, 100 (2004), 379-410.
• [23] S. Ulbrich, On the superlinear local convergence of a filter-SQP method. Math. Programming, 100 (2004), 217-245.
• [24] A. Wächter and L.T. Biegler, Line search filter methods for nonlinear programming: motivation and global convergence. SIAM J. Optimization, 16 (2005), 1-31.
• [25] A. Wächter and L.T. Biegler, On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming. Math. Programming, 106 (2006), 25-57.

*

This work was supported by FAPESP grant 2004/05891-1.

1

We say that a point is regular if the linear independence constraint qualification (LICQ) holds.

Publication Dates

• Publication in this collection
  22 Nov 2007
• Date of issue
  2007

History

• Accepted
  12 Apr 2007
• Received
  05 Mar 2007