# The Flattened Aggregate Constraint Homotopy Method for Nonlinear Programming Problems with Many Nonlinear Constraints

Authors: Zhengyong Zhou; Bo Yu
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

The aggregate constraint homotopy method uses a single smoothing constraint instead ofm-constraints to reduce the dimension of its homotopy map, and hence it is expected to be more efficient than the combined homotopy interior point method when the number of constraints is very large. However, the gradient and Hessian of the aggregate constraint function are complicated combinations of gradients and Hessians of all constraint functions, and hence they are expensive to calculate when the number of constraint functions is very large. In order to improve the performance of the aggregate constraint homotopy method for solving nonlinear programming problems, with few variables and many nonlinear constraints, a flattened aggregate constraint homotopy method, that can save much computation of gradients and Hessians of constraint functions, is presented. Under some similar conditions for other homotopy methods, existence and convergence of a smooth homotopy path are proven. A numerical procedure is given to implement the proposed homotopy method, preliminary computational results show its performance, and it is also competitive with the state-of-the-art solver KNITRO for solving large-scale nonlinear optimization.


---

## Body

## 1. Introduction

In this paper, we consider the following nonlinear programming problem:(1)min⁡f(x),s.t.g(x)≤0,
where x∈Rn is the variable, g(x)=(g1(x),…,gm(x))T, f:Rn→R, and gi(x):Rn→R, i=1,…,m, are three times continuously differentiable, and m is very large, but n is moderate. It has wide applications, and a typical situation is the discretized semi-infinite programming problem.From the mid-1980s, much attention has been paid to interior point methods for mathematical programming, and many results on theory, algorithms, and applications on linear programming, convex programming, complementarity problems, semidefinite programming, and linear cone programming were obtained (see monographs [1–6] and references therein). For nonlinear programming, the typical algorithms used were the Newton-type methods to the perturbed first-order necessary conditions combined with line search or trust region methods with a proper merit function (e.g., [7–9]). The general conditions of global convergence for these methods required that the feasible set be bounded and that the Jacobian matrix be uniformly nonsingular. Another typical class of globally convergent methods for nonlinear programming was probability-one homotopy methods (e.g., [10–12]), whose global convergence can be established under some weaker conditions than the ones for Newton-type methods. The excellent feature is that, unlike line search or trust region methods, they do not depend on the descent of a merit function and so are insensitive to the local minimum of the merit function, in which any search direction is not a descent direction of the merit function.In [10, 11], Feng et al. proposed a homotopy method for nonlinear programming (1), which was called the combined homotopy interior point (abbreviated by CHIP) method; its global convergence was proven under the normal cone condition (see below for its definition) for the feasible set as well as some common conditions. On the basis of the CHIP method, some modified CHIP methods were presented in [13, 14]; the global convergence was established under the quasinormal cone and pseudocone condition for the feasible set, respectively. In [12], Watson described some probability-one homotopy methods for the unconstrained and inequality constrained optimization, whose global convergence was established under some weaker assumptions. Recently, Yu and Shang proposed a constraint shifting combined homotopy method in [15, 16], in which not only the objective function but also the constraint functions were regularly deformed. The global convergence was proven under the condition that the initial feasible set, which approaches the feasible set of (1) as the homotopy parameter changes from 1 to 0, not necessarily the feasible set of (1), satisfies the normal cone condition.Letgmax⁡(x)=max⁡1≤i≤m⁡{gi(x)}; then (1) is equivalent to
(2)min⁡f(x),s.t.gmax⁡(x)≤0,
which has only one, but nonsmooth, constraint. In [17], the following aggregate function was introduced, which is a smooth approximation of gmax⁡(x) with a smoothing parameter t>0 and induced from the max-entropy theory:
(3)g^(x,t)=tln⁡⁡∑i=1mexp⁡⁡(gi(x)t).
It is also known as exponential penalty function (see [18]). By using it for all constraint functions of problem (1), an aggregate constraint homotopy (abbreviated by ACH) method was presented by Yu et al. in [19], whose global convergence was obtained under the condition that the feasible set satisfies the weak normal cone condition. Although the ACH method has only one smoothing constraint, the gradient and Hessian of the aggregate constraint function
(4)∇xg^(x,t)=∑i=1m‍ci(x,t)∇gi(x),∂g^(x,t)∂t=1t(g^(x,t)-∑i=1m‍ci(x,t)gi(x)),∇xx2g^(x,t)=∑i=1m‍ci(x,t)∇2gi(x)+1t∑i=1m‍ci(x,t)∇gi(x)(∇gi(x))T-1t∑i=1m‍ci(x,t)∇gi(x)∑i=1m‍ci(x,t)(∇gi(x))T,∇xt2g^(x,t)=1t2(∑i=1m‍ci(x,t)∇gi(x)∑i=1m‍ci(x,t)gi(x)-∑i=1m‍ci(x,t)gi(x)∇gi(x)),
where
(5)ci(x,t)=exp⁡⁡(gi(x)/t)∑i=1mexp⁡⁡(gi(x)/t)∈(0,1),∑i=1mci(x,t)=1,
are complicated combinations of gradients and Hessians of all constraint functions and, hence, are expensive to calculate when m is very large.Throughout this paper, we assume that the nonlinear programming problem (1) possesses a very large number of nonlinear constraints, but a small number of variables, and the objective and constraint functions are not sparse. For such a problem, the number of constraint functions can be so large that the computation of gradients and Hessians of all constraint functions is very expensive and cannot be stored in memory and, hence, the general numerical methods for solving nonlinear programming are not efficient. Although active set methods only need to calculate gradients and Hessians of a part of constraint functions, require lower storage, and have faster numerical solution, the working set is difficult to estimate without knowing the internal structure of the problem.In this paper, we present a new homotopy method called the flattened aggregate constraint homotopy (abbreviated by FACH) method for nonlinear programming (1) by using a new smoothing technique, in which only a part of constraint functions is aggregated. Under the normal cone condition for the feasible set and some other general assumptions, we prove that the FACH method can determine a smooth homotopy path from a given interior point of the feasible set to a KKT point of (1), and preliminary numerical results demonstrate its efficiency.The rest of this paper is organized as follows. We conclude this section with some notations, definitions, and a lemma. The flattened aggregate constraint function with some properties is given in Section2. The homotopy map, existence and convergence of a smooth homotopy path with proof are given in Section 3. A numerical procedure for tracking the smooth homotopy path, and numerical test results with some remarks are given in Section 4. Finally, we conclude the paper with some remarks in Section 5.When discussing scalars and scalar-valued functions, subscripts refer to iteration step so that superscripts can be used for exponentiation. In contrast, for vectors and vector-valued functions, subscripts are used to indicate components, whereas superscripts are used to indicate the iteration step. The identity matrix is represented byI.Unless otherwise specified,∥·∥ denotes the Euclidean norm. Ω={x∈Rn∣g(x)≤0} is the feasible set of (1), whereas Ω0={x∈Rn∣g(x)<0} is the interior of Ω and ∂Ω=Ω∖Ω0 is the boundary of Ω. The symbols R+m and R++m denote the nonnegative and positive quadrants of Rm, respectively. The active index set is denoted by I(x)={i∈{1,…,m}∣gi(x)=0} at x∈Ω. For a function F:Rn→Rm, ∇F(x) is the n×m matrix whose (i,j)th element is ∂Fj(x)/∂xi; F-1(Y)={x∣F(x)∈Y} is the set-valued inverse for the set Y⊆Rm.Definition 1 (see [10]).
The setΩ satisfies the normal cone condition if and only if the normal cone of Ω at any x∈∂Ω does not meet Ω0; that is,
(6){x+∑i∈I(x)yi∇gi(x)∣yi≥0}∩Ω0=∅.Definition 2.
If and only if there existsλ*∈R+m such that (x*,λ*) satisfies
(7)∇f(x*)+∇g(x*)λ*=0,λi*gi(x*)=0,λi*≥0,gi(x*)≤0,i=1,…,m,
then x* is called a KKT point of (1) and λ* is the corresponding Lagrangian multiplier.Definition 3.
LetU⊂Rn be an open set and let F:U→Rm be a differentiable map. We say y∈Rm is a regular value of F if and only if
(8)Range[∂F(x)∂x]=Rm,∀x∈F-1(y).Lemma 4 (parameterized Sard theorem [20]).
LetU⊂Rn and V⊂Rk be two open sets and let F:U×V→Rm be a Cr differentiable map with r>max⁡⁡{0,n-m}. If 0∈Rm is a regular value of F, then, for almost all a∈V, 0 is a regular value of Fa=F(·,a).


## 2. The Flattened Aggregate Constraint Function

In this section, we suppose that the following assumptions hold.Assumption 5.
Ω
0 is nonempty and Ω is bounded.Assumption 6.
For anyx∈∂Ω, {∇gi(x)∣i∈I(x)} are positive independent; that is,
(9)∑i∈I(x)yi∇gi(x)=0,yi≥0,i∈I(x)⟹yi=0,i∈I(x).Assumption 7.
Ω satisfies the normal cone condition (see Definition 1).The first assumption is the Slater’s condition and the boundedness of the feasible set, which are two basic conditions. The second assumption provides the regularity of constraints, which is weaker than the linear independence constraint qualification. The last assumption, the normal cone condition of the feasible set, is a generalization of the convex condition for the feasible set. Indeed, ifΩ is a convex set, then it satisfies the normal cone condition. Some simple nonconvex sets, satisfying the normal cone condition, are shown in [10].In this paper, we construct a flattened aggregate constraint function as follows:(10)g^(θ,c1,c2,α)(x,t)=θtln⁡⁡(∑i=1m‍φ(gi(x),t)exp⁡(gi(x)θt)+exp⁡⁡(-ɛ(t)θt)),
where
(11)ɛ(t)=c1t+c2,c1>0, c2>0, 0<θ≤1 are some adjusting parameters, φ(z,t):R×[0,1]→[0,1] is a C3 bivariate function satisfying φ(z,t)=0 for z≤-αɛ(t) with α>1, and φ(z,t)=1 for z≥-ɛ(t). For simplicity of discussion, we write g^(θ,c1,c2,α)(x,t) as g^(x,t). By its definition, φ(z,t) can be expressed as
(12)φ(z,t)={0,forz≤-αɛ(t)withα>1,1,forz≥-ɛ(t),φ¯(z,t),otherwise,
where, for any t∈[0,1], φ¯(z,t) satisfies
(13)φ¯(-αɛ(t),t)=0,∂φ¯(-αɛ(t),t)∂z=0,∂2φ¯(-αɛ(t),t)∂z2=0,∂3φ¯(-αɛ(t),t)∂z3=0,φ¯(-ɛ(t),t)=1,∂φ¯(-ɛ(t),t)∂z=0,∂2φ¯(-ɛ(t),t)∂z2=0,∂3φ¯(-ɛ(t),t)∂z3=0.
There exist many kinds of functions φ¯(z,t), such as polynomial functions and spline functions. Since polynomial functions have simple mathematical expressions and many excellent properties, throughout this paper, φ¯(z,t) is chosen as
(14)φ¯(z,t)=-20(z+ɛ(t))7(αɛ(t)-ɛ(t))7+-70(z+ɛ(t))6(αɛ(t)-ɛ(t))6+-84(z+ɛ(t))5(αɛ(t)-ɛ(t))5+-35(z+ɛ(t))4(αɛ(t)-ɛ(t))4+1.LetIɛ(x,t)={i∣gi(x)>-αɛ(t)}; by the definition of φ(z,t), g^(x,t) in (10) can be rewritten as
(15)g^(x,t)=θtln⁡⁡(∑i∈Iɛ(x,t)φ(gi(x),t)exp⁡⁡(gi(x)θt)hhhhhhh+exp⁡⁡(-ɛ(t)θt)∑i∈Iɛ(x,t)φ(gi(x),t)exp⁡⁡(gi(x)θt)).
The gradient of g^(x,t) with respect to x is
(16)∇xg^(x,t)=∑i∈Iɛ(x,t)(λi(x,t)+ξi(x,t))∇gi(x),
where
(17)λi(x,t)=φ(gi(x),t)exp⁡(gi(x)/θt)∑j∈Iɛ(x,t)φ(gj(x),t)exp⁡⁡(gj(x)/θt)+exp⁡⁡(-ɛ(t)/θt),(18)ξi(x,t)=θt∂φ(gi(x),t)/∂zexp⁡⁡(gi(x)/θt)∑j∈Iɛ(x,t)φ(gj(x),t)exp⁡⁡(gj(x)/θt)+exp⁡⁡(-ɛ(t)/θt)‍.
Then the gradient and Hessian of g^(x,t) only relate to a part of constraint functions; that is, i∈Iɛ(x,t).Proposition 8.
For anyθ>0, t∈(0,1], x∈Ω with θt↓0, t→t¯∈[0,1], x→x¯∈∂Ω,(1)
λi(x,t)→0 for  i∉I(x¯);(2)
∑i∈I(x¯)λi(x,t)→1;(3)
ξ(x,t)→0.Proof.
(1) Becausex¯∈∂Ω, there exists i*∈I(x¯) such that gi*(x¯)=0. By the continuity of φ(z,t) and g(x), we know that lim⁡x→x¯,t→t¯⁡φ(gi*(x),t)=φ(gi*(x¯),t¯)=1 and lim⁡x→x¯⁡(gi(x)-gi*(x))=gi(x¯)<0 for i∉I(x¯); hence
(19)0≤λi(x,t)≤exp⁡⁡(gi(x)/θt)φ(gi*(x),t)exp⁡⁡(gi*(x)/θt)=exp⁡⁡((gi(x)-gi*(x))/θt)φ(gi*(x),t)⟶0.
Hence, item 1 is satisfied for θt↓0, t→t¯, and x→x¯.
(2) Bylim⁡x→x¯,t→t¯⁡(-ɛ(t)-gi*(x))=-ɛ(t¯)<0, we have
(20)1≥∑i∈Iɛ(x,t)‍λi(x,t)=∑i∈Iɛ(x,t)φ(gi(x),t)exp⁡⁡(gi(x)/θt)∑j∈Iɛ(x,t)φ(gj(x),t)exp⁡⁡(gj(x)/θt)+exp⁡⁡(-ɛ(t)/θt)≥φ(gi*(x),t)exp⁡⁡(gi*(x)/θt)φ(gi*(x),t)exp⁡⁡(gi*(x)/θt)+exp⁡⁡(-ɛ(t)/θt)≥φ(gi*(x),t)1+exp⁡⁡((-ɛ(t)-gi*(x))/θt)⟶1.
Hence, ∑i∈Iɛ(x,t)λi(x,t)→1. Together with item 1, we have ∑i∈I(x¯)λi(x,t)→1 for θt↓0, t→t¯, and x→x¯.
(3) Ifgi(x)≥-ɛ(t), then ∂φ(gi(x),t)/∂z=0 by its definition, and hence
(21)ξi(x,t)=0;
else
(22)|ξi(x,t)|≤θt|∂φ(gi(x),t)/∂z|exp⁡⁡(gi(x)/θt)exp⁡⁡(-ɛ(t)/θt)≤θtM,
where M=max⁡(z,t)∈R×[0,1]⁡|∂φ(z,t)/∂z|. Then, we have ξ(x,t)→0 for θt↓0, t→t¯, and x→x¯.Proposition 9.
For any given0<θ≤1,(1)
-ɛ(t)≤g^(x,t)≤max⁡{gmax⁡(x),-ɛ(t)}+θtln⁡(mIɛ+1),(2)
gmax⁡(x)≤g^(x,t)≤gmax⁡(x)+θtln⁡(mIɛ+1) when  gmax⁡(x)≥-ɛ(t),(3)
g^(x,t)=-ɛ(t) when  gmax⁡(x)≤-αɛ(t) with  α>1,where mIɛ denotes the cardinality of Iɛ(x,t).Proof.
(1) Byφ(z,t)≥0, we have
(23)g^(x,t)=θtln⁡⁡(∑i∈Iɛ(x,t)φ(gi(x),t)exp⁡⁡(gi(x)θt)hhhhhhh+exp⁡⁡(-ɛ(t)θt)∑i∈Iɛ(x,t)φ(gi(x),t)exp⁡⁡(gi(x)θt))≥θtln⁡⁡(exp⁡(-ɛ(t)θt))=-ɛ(t).
Then the left inequality can be obtained. By φ(z,t)≤1, we have
(24)g^(x,t)=θtln⁡⁡(∑i∈Iɛ(x,t)φ(gi(x),t)exp⁡(gi(x)θt)hhhhhhhh+exp⁡(-ɛ(t)θt)∑i∈Iɛ(x,t)φ(gi(x),t)exp⁡(gi(x)θt))≤θtln⁡⁡(∑i∈Iɛ(x,t)exp⁡⁡(gi(x)θt)+exp⁡⁡(-ɛ(t)θt)).
Then the right inequality can be obtained by Proposition 2.3 in [19].
(2) For anyx,t with gmax⁡(x)≥-ɛ(t), let gi*(x)=gmax⁡(x); then φ(gi*(x),t)=1 by its definition, and hence
(25)g^(x,t)=θtln⁡⁡(∑i∈Iɛ(x,t)φ(gi(x),t)exp⁡⁡(gi(x)θt)hhhhhhhh+exp⁡(-ɛ(t)θt)∑i∈Iɛ(x,t)φ(gi(x),t)exp⁡⁡(gi(x)θt))≥θtln⁡⁡(φ(gi*(x),t)exp⁡⁡(gi*(x)θt))=θtln⁡⁡(exp⁡(gi*(x)θt))=gi*(x)=gmax⁡(x);
then the left inequality is true. The right inequality can be obtained by item 1.
(3) It is trivial by the definition ofg^(x,t).Proposition 10.
LetΩθ(t)={x∣g^(x,t)≤0} and Ωθ(t)0={x∣g^(x,t)<0}; one has(1)
Ωθ(t)⊂Ω for  t∈(0,1];(2)
for any bounded and closed set  Ω^⊂Ω0, there exists a  θ∈(0,1], such that  Ω^⊂Ωθ(1)0;(3)
Ωθ(t)→Ω as  t↓0.Proof.
(1) Ifx∉Ω, which is equivalent to gmax⁡(x)>0, by Proposition 9(2),
(26)0<gmax⁡(x)≤g^(x,t),
which means x∉Ωθ(t); then we have Ωθ(t)⊂Ω.
(2) By the continuity ofg(x) and the fact that Ω^ is a bounded closed set, there exists a point x*∈Ω^ such that gmax⁡(x) reaches its maximum in Ω^, and gmax⁡(x*)<0. For any θ∈(0,1] satisfying θ≤-max⁡{gmax⁡(x*),-ɛ(1)}/(2ln⁡(m+1)), by Proposition 9(1), for any x∈Ω^,
(27)g^(x,1)≤max⁡⁡{gmax⁡(x),-ɛ(1)}+θln⁡(m+1)≤max⁡⁡{gmax⁡(x*),-ɛ(1)}+θln⁡⁡(m+1)=12max⁡⁡{gmax⁡(x*),-ɛ(1)}<0,
which means Ω^⊂Ωθ(1)0.
(3) For anyx∈Ω, that is, gmax⁡(x)≤0, by the right inequality of Proposition 9(1),
(28)lim⁡t↓0g^(x,t)≤lim⁡t↓0max⁡⁡{gmax⁡(x),-ɛ(t)}≤max⁡⁡{gmax⁡(x),-c2}≤0,
where the second inequality comes from -ɛ(t)=-c1t-c2≤-c2<0 for t∈[0,1], which means Ω⊂lim⁡t↓0⁡Ωθ(t). Then, together with item 1,
(29)lim⁡t↓0Ωθ(t)=Ω.Proposition 11.
There exists aθ∈(0,1] such that, for any t∈(0,1], x∈∂Ωθ(t); then ∇xg^(x,t)≠0.Proof.
If not, for anyθ∈(0,1], there exist corresponding tθ∈(0,1] and xθ∈∂Ωθ(tθ) such that ∇xg^(xθ,tθ)=0, which means that there must exist three sequences {θk}k=1∞⊂(0,1], {tk}k=1∞⊂(0,1], and {xk}k=1∞ with xk∈∂Ωθk(tk), such that θk↓0, tk→t¯, xk→x¯∈∂Ω, and
(30)∇xg^(xk,tk)=∑i∈Iɛ(xk,tk)(λi(xk,tk)+ξi(xk,tk))∇gi(xk)⟶0
as k→∞. Then, by Proposition 8, we know that λ(xk,tk)→λ¯, ξ(xk,tk)→0 as k→∞, and
(31)∑i∈lim⁡k→∞⁡Iɛ(xk,tk)λ¯i=∑i∈I(x¯)λ¯i=1,
where the first equality comes from that λ¯i=0 for i∉I(x¯). Therefore, we have ∑i∈I(x¯)λ¯i∇gi(x¯)=0 by taking limits on (30); this is a contradiction to Assumption 6.Proposition 12.
For any bounded closed setΩ^⊂Ω0, there exists a θ∈(0,1] such that
(32){x+λ∇xg^(x,1)∣λ>0}∩Ω^=∅,
for any x∈∂Ωθ(1).Proof.
If not, there exist a bounded closed setΩ^ and four sequences {θk}k=1∞⊂(0,1], {xk}k=1∞ with xk∈∂Ωθk(1), {x^k}k=1∞⊂Ω^, and {λk}k=1∞>0 such that θk↓0, xk→x¯∈∂Ω, x^k→x^∈Ω^, λi(xk,1)→λi*≥0, ξi(xk,1)→0 as k→∞, and
(33)x^k=xk+λk∑i∈Iɛ(xk,1)(λi(xk,1)+ξi(xk,1))∇gi(xk).
By Proposition 8, we have
(34)∑i∈lim⁡k→∞⁡Iɛ(xk,1)λi*=∑i∈I(x¯)λi*=1,
where the first equality comes from that λi*=0 for i∉I(x¯). Then, we have
(35)∑i∈Iɛ(xk,1)(λi(xk,1)+ξi(xk,1))∇gi(xk)⟶∑i∈I(x¯)λi*∇gi(x¯).
By Assumption 6 and using (34),
(36)∑i∈I(x¯)λi*∇gi(x¯)≠0.
Hence {λk}k=1∞ is a bounded sequence, and there must exist a subsequence converging to λ¯>0; then we have by (33)
(37)x^=x¯+λ¯∑i∈I(x¯)λi*∇gi(x¯),
which contradicts Assumption 7.


## 3. The Flattened Aggregate Constraint Homotopy Method

In [19], the following aggregate constraint homotopy was introduced:
(38)H(x,λ,t)=((1-t)(∇f(x)+λ∇xg^θ(x,t))+t(x-x0)λg^θ(x,t)-tλ0g^θ(x0,1))=0,
where g^θ(x,t)=θtln⁡⁡∑i=1mexp⁡⁡(gi(x)/θt) with θ∈(0,1], λ is the Lagrangian multiplier of the aggregate constraint function g^θ(x,t), and (x0,λ0) is the starting point in Ω^×R++. Under the weak normal cone condition, which is similar to the normal cone condition, it was proved that the ACH determines a smooth interior path from a given interior point to a KKT point. Then, the predictor-corrector procedure can be applied to trace the homotopy path from (x0,λ0,1) to (x*,λ*,0), in which (x*,λ*) is a KKT point of (1).


### 3.1. The Flattened Aggregate Constraint Homotopy

Using the flattened aggregate constraint functiong^(x,t) in (10), we construct the following flattened aggregate constraint homotopy:
(39)H(x,λ,t)=((1-t)(∇f(x)+λ∇xg^(x,t))+t(x-x0)λg^(x,t)-tλ0g^(x0,1))=0,
where λ is the Lagrangian multiplier of the flattened aggregate constraint and (x0,λ0) is the starting point and can be randomly chosen from Ω^×R++.We give the main theorem on the existence and convergence of a smooth path from(x0,λ0,1) to (x*,λ*,0), in which (x*,λ*) is a solution of the KKT system of (1), and hence the global convergence of our proposal, namely, the FACH method, can be proven.Theorem 13.
Suppose that Assumptions5–7 hold and θ∈(0,1] satisfies Propositions 10–12 for the bounded closed set Ω^⊂Ωθ(1)0⊂Ω0; then for almost all starting points w0=(x0,λ0)∈Ω^×R++, the zero-points set H-1(0)={(x,λ,t)∣H(x,λ,t)=0} defines a smooth path Γw0, which starts at (x0,λ0,1) and approaches the hyperplane t=0. Furthermore, let (x*,λ*,0) be any limit point of Γw0, and y*=λ*λ(x*,0)  (λ(x*,0) is a limit of λ(x,t) as x→x*,t↓0), then x* is a KKT point of (1) and y* is the corresponding Lagrangian multiplier.Proof.
ConsiderH as a map of the variable (w0,x,λ,t), for any (w0,t)∈Ω^×R++×(0,1], which means x0∈Ωθ(1)0; hence g^(x0,1)<0, because
(40)∂H(w0,x,λ,t)∂w0=(-tI0-tλ0(∇xg^(x0,1))T-tg^(x0,1))
is nonsingular and the Jacobi matrix of H(w0,x,λ,t) is of full row rank. Using the parameterized Sard theorem, Lemma 4, for almost all w0∈Ω^×R++, 0 is a regular value of H(x,λ,t). From
(41)H(x,λ,1)=(x-x0λg^(x,1)-λ0g^(x0,1)),
we know that (x0,λ0) is the unique and simple solution of H(x,λ,1)=0. Then, by the implicit function theorem, there exists a smooth curve Γw0⊂Ωθ(t)0×R++×(0,1] starting from (x0,λ0,1) and being transversal to the hyperplane t=1.
SinceΓw0 can be extended in Ωθ(t)0×R++×(0,1) until it converges to the boundary of Ωθ(t)0×R++×(0,1), there must exist an extreme point (x¯,λ¯,t¯)∈∂(Ωθ(t¯)×R+×[0,1]). Let (x¯,λ¯,t¯) be any extreme point of Γw0 other than (x0,λ0,1); then only the following five cases are possible:(1)
(x¯,λ¯,t¯)∈Ω×R+×{0}, λ¯<+∞;(2)
(x¯,λ¯,t¯)∈Ωθ(1)×R+×{1}, λ¯<+∞;(3)
(x¯,λ¯,t¯)∈∂Ωθ(t¯)×R++×(0,1), λ¯<+∞;(4)
(x¯,λ¯,t¯)∈Ωθ(t¯)×{0}×(0,1);(5)
(x¯,λ¯,t¯)∈Ωθ(t¯)×{+∞}×[0,1].
Case (2) means thatH(x,λ,1)=0 has another solution (x¯,λ¯) except (x0,λ0), or (x0,λ0) is a double solution of H(x,λ,1)=0, which contradicts the fact that (x0,λ0) is the unique and simple solution of H(x,λ,1)=0. Case (3) means g^(x¯,t¯)=0, 0<λ¯<∞, and 0<t¯<1; hence
(42)λ¯g^(x¯,t¯)-t¯λ0g^(x0,1)=-t¯λ0g^(x0,1)<0,
which contradicts the last equation of (39). Case (4) means g^(x¯,t¯)≤0, λ¯=0, and 0<t¯<1; hence
(43)λ¯g^(x¯,t¯)-t¯λ0g^(x0,1)=-t¯λ0g^(x0,1)<0,
which contradicts the last equation of (39). Then, Cases (2), (3), and (4) are impossible.
BecauseH(x,λ,1)=0 has only one and also unique solution (x0,λ0), Case (2) is impossible. By the continuity and the last equation of (39), we know that Cases (3) and (4) are impossible.
If Case (5) holds, there must exist a sequence{(xk,λk,tk)}k=1∞ on Γw0 such that xk→x¯, λk→∞, and tk→t¯. By the last equation of (39), we have
(44)g^(xk,tk)=tkλ0g^(x0,1)λk⟶0
as k→∞, which means x¯∈∂Ωθ(t¯). The following two cases may be possible.(1)
Ift¯=1, let α¯ be any accumulation point of {(1-tk)λk}, which is known to exist from (39) and lim⁡k→∞⁡∇xg^(xk,tk)≠0; then we have
(45)x0=x¯+α¯∇xg^(x¯,1).
If α¯=0, then x¯=x0 contradicts x0∈Ωθ(1)0; else contradicts Proposition 12.(2)
Ift¯<1, by (39), we have
(46)(1-t¯)lim⁡k→∞λk∇xg^(xk,tk)=-(1-t¯)∇f(x¯)-t¯(x¯-x0),
the right-hand is finite; however, by lim⁡k→∞⁡∇xg^(xk,tk)≠0, we know that the left-hand is infinite; this is a contradiction.
As a conclusion, Case (1) is the only possible case. This implies thatΓw0 must approach the hyperplane t=0. Because {λ(xk,tk)}k=1∞(defined in (17)), {xk}k=1∞, and {λk}k=1∞ are bounded sequences, we know that {(xk,λk,λ(xk,tk))}k=1∞ has at least one accumulation point as k→∞. Let (x*,λ*,λ(x*,0)) be any accumulation point; yi*=λ*λi(x*,0). By (39) and the fact that ξ(x,t)→0 as t↓0, we have
(47)∇f(x*)+∑i=1myi*∇gi(x*)=0.
By the fact thatlim⁡t↓0⁡Ωθ(t)=Ω, we have x*∈Ω. If x*∈Ω0, we know that
(48)lim⁡x→x*,t↓0g^(x,t)≤lim⁡x→x*,t↓0(max⁡⁡{gmax⁡(x),-ɛ(t)}+θtln⁡⁡(mIɛ+1))=max⁡{gmax⁡(x*),ɛ(0)}=max⁡{gmax⁡(x*),-αc2}<0,
where the first inequality comes from Proposition 9(1), the second equality comes from the continuity, and the third equality comes from the definition of ɛ(t) in (11); hence λ*=0 by the last equation of (39); else x*∈∂Ω, lim⁡x→x*,t↓0⁡g^(x,t)=gmax⁡(x*)=0 by Proposition 9(2) and λi(x*,0)=0 for i∉I(x*) by Proposition 8. Thus we have
(49)yi*gi(x*)=0,i=1,…,m.
Summing up,(x*,y*) is a solution of the KKT system of (1), which means that x* is a KKT point of (1) and y* is the corresponding Lagrangian multiplier.


### 3.2. The Modified Flattened Aggregate Constraint Homotopy

From Proposition8(3), we know that ξ(x,t)→0 as t↓0 for any x∈Ω and, hence, we can use the following modified flattened aggregate constraint homotopy (MFACH) instead of the FACH:
(50)H(x,λ,t)=((1-t)(∇f(x)+λh(x,t))+t(x-x0)λg^(x,t)-tλ0g^(x0,1))=0,
where
(51)h(x,t)=∑i∈Iɛ(x,t)λi(x,t)∇gi(x).Remarks. (i) Since ξ(x,t) is dropped in (51), the expressions of the homotopy map (50) and its Jacobian and hence the corresponding code become simpler. Moreover, the computation of H in (50) and its Jacobi matrix DH is a little cheaper than that for the FACH method.(ii) To guarantee thatH(x,λ,t) in (39) be a C2 map, φ(z,t) must be C3; hence, if  φ¯(z,t) is chosen as a polynomial, the total degree should be seven; in contrast, for the homotopy map in (50), φ(z,t) is only needed to be C2. Then φ¯(z,t) should satisfy
(52)φ¯(-αɛ(t),t)=0,∂φ¯(-αɛ(t),t)∂z=0,∂2φ¯(-αɛ(t),t)∂z2=0,φ¯(-ɛ(t),t)=1,∂φ¯(-ɛ(t),t)∂z=0,∂2φ¯(-ɛ(t),t)∂z2=0
and φ¯(z,t) can be chosen as
(53)φ¯(z,t)=6(z+ɛ(t))5(αɛ(t)-ɛ(t))5+15(z+ɛ(t))4(αɛ(t)-ɛ(t))4+10(z+ɛ(t))3(αɛ(t)-ɛ(t))3+1.(iii) The existence and convergence of the homotopy path defined by the MFACH can be proven in a similar way with that of the FACH.


## 3.1. The Flattened Aggregate Constraint Homotopy

Using the flattened aggregate constraint functiong^(x,t) in (10), we construct the following flattened aggregate constraint homotopy:
(39)H(x,λ,t)=((1-t)(∇f(x)+λ∇xg^(x,t))+t(x-x0)λg^(x,t)-tλ0g^(x0,1))=0,
where λ is the Lagrangian multiplier of the flattened aggregate constraint and (x0,λ0) is the starting point and can be randomly chosen from Ω^×R++.We give the main theorem on the existence and convergence of a smooth path from(x0,λ0,1) to (x*,λ*,0), in which (x*,λ*) is a solution of the KKT system of (1), and hence the global convergence of our proposal, namely, the FACH method, can be proven.Theorem 13.
Suppose that Assumptions5–7 hold and θ∈(0,1] satisfies Propositions 10–12 for the bounded closed set Ω^⊂Ωθ(1)0⊂Ω0; then for almost all starting points w0=(x0,λ0)∈Ω^×R++, the zero-points set H-1(0)={(x,λ,t)∣H(x,λ,t)=0} defines a smooth path Γw0, which starts at (x0,λ0,1) and approaches the hyperplane t=0. Furthermore, let (x*,λ*,0) be any limit point of Γw0, and y*=λ*λ(x*,0)  (λ(x*,0) is a limit of λ(x,t) as x→x*,t↓0), then x* is a KKT point of (1) and y* is the corresponding Lagrangian multiplier.Proof.
ConsiderH as a map of the variable (w0,x,λ,t), for any (w0,t)∈Ω^×R++×(0,1], which means x0∈Ωθ(1)0; hence g^(x0,1)<0, because
(40)∂H(w0,x,λ,t)∂w0=(-tI0-tλ0(∇xg^(x0,1))T-tg^(x0,1))
is nonsingular and the Jacobi matrix of H(w0,x,λ,t) is of full row rank. Using the parameterized Sard theorem, Lemma 4, for almost all w0∈Ω^×R++, 0 is a regular value of H(x,λ,t). From
(41)H(x,λ,1)=(x-x0λg^(x,1)-λ0g^(x0,1)),
we know that (x0,λ0) is the unique and simple solution of H(x,λ,1)=0. Then, by the implicit function theorem, there exists a smooth curve Γw0⊂Ωθ(t)0×R++×(0,1] starting from (x0,λ0,1) and being transversal to the hyperplane t=1.
SinceΓw0 can be extended in Ωθ(t)0×R++×(0,1) until it converges to the boundary of Ωθ(t)0×R++×(0,1), there must exist an extreme point (x¯,λ¯,t¯)∈∂(Ωθ(t¯)×R+×[0,1]). Let (x¯,λ¯,t¯) be any extreme point of Γw0 other than (x0,λ0,1); then only the following five cases are possible:(1)
(x¯,λ¯,t¯)∈Ω×R+×{0}, λ¯<+∞;(2)
(x¯,λ¯,t¯)∈Ωθ(1)×R+×{1}, λ¯<+∞;(3)
(x¯,λ¯,t¯)∈∂Ωθ(t¯)×R++×(0,1), λ¯<+∞;(4)
(x¯,λ¯,t¯)∈Ωθ(t¯)×{0}×(0,1);(5)
(x¯,λ¯,t¯)∈Ωθ(t¯)×{+∞}×[0,1].
Case (2) means thatH(x,λ,1)=0 has another solution (x¯,λ¯) except (x0,λ0), or (x0,λ0) is a double solution of H(x,λ,1)=0, which contradicts the fact that (x0,λ0) is the unique and simple solution of H(x,λ,1)=0. Case (3) means g^(x¯,t¯)=0, 0<λ¯<∞, and 0<t¯<1; hence
(42)λ¯g^(x¯,t¯)-t¯λ0g^(x0,1)=-t¯λ0g^(x0,1)<0,
which contradicts the last equation of (39). Case (4) means g^(x¯,t¯)≤0, λ¯=0, and 0<t¯<1; hence
(43)λ¯g^(x¯,t¯)-t¯λ0g^(x0,1)=-t¯λ0g^(x0,1)<0,
which contradicts the last equation of (39). Then, Cases (2), (3), and (4) are impossible.
BecauseH(x,λ,1)=0 has only one and also unique solution (x0,λ0), Case (2) is impossible. By the continuity and the last equation of (39), we know that Cases (3) and (4) are impossible.
If Case (5) holds, there must exist a sequence{(xk,λk,tk)}k=1∞ on Γw0 such that xk→x¯, λk→∞, and tk→t¯. By the last equation of (39), we have
(44)g^(xk,tk)=tkλ0g^(x0,1)λk⟶0
as k→∞, which means x¯∈∂Ωθ(t¯). The following two cases may be possible.(1)
Ift¯=1, let α¯ be any accumulation point of {(1-tk)λk}, which is known to exist from (39) and lim⁡k→∞⁡∇xg^(xk,tk)≠0; then we have
(45)x0=x¯+α¯∇xg^(x¯,1).
If α¯=0, then x¯=x0 contradicts x0∈Ωθ(1)0; else contradicts Proposition 12.(2)
Ift¯<1, by (39), we have
(46)(1-t¯)lim⁡k→∞λk∇xg^(xk,tk)=-(1-t¯)∇f(x¯)-t¯(x¯-x0),
the right-hand is finite; however, by lim⁡k→∞⁡∇xg^(xk,tk)≠0, we know that the left-hand is infinite; this is a contradiction.
As a conclusion, Case (1) is the only possible case. This implies thatΓw0 must approach the hyperplane t=0. Because {λ(xk,tk)}k=1∞(defined in (17)), {xk}k=1∞, and {λk}k=1∞ are bounded sequences, we know that {(xk,λk,λ(xk,tk))}k=1∞ has at least one accumulation point as k→∞. Let (x*,λ*,λ(x*,0)) be any accumulation point; yi*=λ*λi(x*,0). By (39) and the fact that ξ(x,t)→0 as t↓0, we have
(47)∇f(x*)+∑i=1myi*∇gi(x*)=0.
By the fact thatlim⁡t↓0⁡Ωθ(t)=Ω, we have x*∈Ω. If x*∈Ω0, we know that
(48)lim⁡x→x*,t↓0g^(x,t)≤lim⁡x→x*,t↓0(max⁡⁡{gmax⁡(x),-ɛ(t)}+θtln⁡⁡(mIɛ+1))=max⁡{gmax⁡(x*),ɛ(0)}=max⁡{gmax⁡(x*),-αc2}<0,
where the first inequality comes from Proposition 9(1), the second equality comes from the continuity, and the third equality comes from the definition of ɛ(t) in (11); hence λ*=0 by the last equation of (39); else x*∈∂Ω, lim⁡x→x*,t↓0⁡g^(x,t)=gmax⁡(x*)=0 by Proposition 9(2) and λi(x*,0)=0 for i∉I(x*) by Proposition 8. Thus we have
(49)yi*gi(x*)=0,i=1,…,m.
Summing up,(x*,y*) is a solution of the KKT system of (1), which means that x* is a KKT point of (1) and y* is the corresponding Lagrangian multiplier.


## 3.2. The Modified Flattened Aggregate Constraint Homotopy

From Proposition8(3), we know that ξ(x,t)→0 as t↓0 for any x∈Ω and, hence, we can use the following modified flattened aggregate constraint homotopy (MFACH) instead of the FACH:
(50)H(x,λ,t)=((1-t)(∇f(x)+λh(x,t))+t(x-x0)λg^(x,t)-tλ0g^(x0,1))=0,
where
(51)h(x,t)=∑i∈Iɛ(x,t)λi(x,t)∇gi(x).Remarks. (i) Since ξ(x,t) is dropped in (51), the expressions of the homotopy map (50) and its Jacobian and hence the corresponding code become simpler. Moreover, the computation of H in (50) and its Jacobi matrix DH is a little cheaper than that for the FACH method.(ii) To guarantee thatH(x,λ,t) in (39) be a C2 map, φ(z,t) must be C3; hence, if  φ¯(z,t) is chosen as a polynomial, the total degree should be seven; in contrast, for the homotopy map in (50), φ(z,t) is only needed to be C2. Then φ¯(z,t) should satisfy
(52)φ¯(-αɛ(t),t)=0,∂φ¯(-αɛ(t),t)∂z=0,∂2φ¯(-αɛ(t),t)∂z2=0,φ¯(-ɛ(t),t)=1,∂φ¯(-ɛ(t),t)∂z=0,∂2φ¯(-ɛ(t),t)∂z2=0
and φ¯(z,t) can be chosen as
(53)φ¯(z,t)=6(z+ɛ(t))5(αɛ(t)-ɛ(t))5+15(z+ɛ(t))4(αɛ(t)-ɛ(t))4+10(z+ɛ(t))3(αɛ(t)-ɛ(t))3+1.(iii) The existence and convergence of the homotopy path defined by the MFACH can be proven in a similar way with that of the FACH.


## 4. The FACH-S-N Procedure and Numerical Results

### 4.1. The FACH-S-N Procedure

In this section, we give a numerical procedure, FACH-S-N procedure, to trace the flattened aggregate constraint homotopy path by secant predictor and Newton corrector steps. It consists of three main steps: the predictor step, the corrector step, and the end game.The predictor step is an approximate step along the homotopy path: it uses a predictor directiondk and a steplength hk to get a predictor point. The first predictor direction uses the tangent direction, and others use the secant direction. The steplength is determined by several parameters. It is set to no more than 1, which can ensure that the predictor point is close to the homotopy path and hence stays in the convergence domain of Newton’s method in the corrector step. If the angle of the predictor direction and the previous one βk is greater than π/4, the steplength will be decreased to avoid using the opposite direction as the predictor direction. If the corrector criteria are satisfied with no more than four Newton iterations for three times in succession, the steplength will be increased or kept invariable. Otherwise, the steplength will be decreased.Once a predictor point(w(k,0),t(k,0)) is calculated, one or more Newton iterations are used to bring the predictor point back to the homotopy path in the corrector step. The corrector points (w(k,i),t(k,i)), i=1,…,5, are calculated by (w(k,i+1),t(k,i+1))=(w(k,i),t(k,i))+d(k,i+1), i=0,…,4, where the step d(k,i+1) is the solution of an augmented system defined by the homotopy equation and the direction perpendicular to the predictor direction. The corrector step terminates when d(k,i) and H(w(k,i),t(k,i)) satisfy the tolerance criterions.At each predictor step and corrector step, the feasibility of the predictor point(wk,tk) and the corrector point (w(k,i),t(k,i)) needs to be checked. If tk<tc in the predictor step or tk<0 in the corrector step, a damping step is used to get a new point (w(k,0),0). Then, if w(k,0) is feasible, the end game, a more efficient strategy than predictor-corrector steps when the homotopy parameter is close to 0, is invoked. Starting with w(k,0), Newton's method is used to solve
(54)Ftc(x,λ)=(∇f(x)+λh(x,tc)λg^(x,tc))=0,
where tc is a small positive constant. For other situations, the steplength will be decreased to make new predictor-corrector steps (see Algorithm 1).Algorithm 1:The FACH-S-N procedure.
Input  Give θ∈(0,1], c1>0, c2>0, α>1, tend and tc; starting point w0∈Ω^×R++; initial steplength h1, steplength
contraction factor Bmin⁡, steplength expansion factor Bmax⁡; tracking tolerances Htol and Hfinal for correction. Set
t0=1, k=1, IT=0, Ngood=2, β1=0.
Step  1. The predictor step.
Step  1.1. Compute the predictor direction.
If k=1, compute d by solving
(H′(x0,λ0,1)(d0)T)d=-d0,
where d0=(0,…,0,-1)∈Rn+2, set the predictor direction d1=d/∥d∥;
else, compute the predictor direction dk, the angle βk between dk and dk-1:
dk=(wk-1,tk-1)-(wk-2,tk-2)∥(wk-1,tk-1)-(wk-2,tk-2)∥,βk=arccos((dk)Tdk-1).
Step  1.2. Adjust the steplength.
Adjust the steplength hk as follows:
if the predictor point is infeasible, set hk=Bmin⁡(1)hk, Ngood=0;
else if the corrector step fails or βk>π/4, set hk=Bmin⁡(1)hk-1, Ngood=0;
else if i=5, set hk=Bmin⁡(2)hk-1, Ngood=0;
else if i=4, set hk=hk-1, Ngood=Ngood+1;
else if i=3, set Ngood=Ngood+1,
if Ngood>2, set hk=min⁡{1,Bmax⁡(2)hk-1}, else, set hk=hk-1;
else, set Ngood=Ngood+1,
if Ngood>2, set hk=min⁡{1,Bmax⁡(1)hk-1}, else, set hk=hk-1.
If hk<10-10, stop the algorithm with an error flag.
Step  1.3. Compute the predictor point and check its feasibility.
Compute the predictor point
(w(k,0),t(k,0))=(wk-1,tk-1)+hkdk.
If w(k,0)∉Ωθ(t(k,0))0×R++ or t(k,0)≥1,
the predictor point is infeasible, goto Step  1.2;
else if t(k,0)≤tend,
adjust the steplength hk, compute the point (w(k,0),0), if w(k,0)∈Ω×R+, goto Step  3, else, the predictor
point is infeasible, goto Step  1.2;
else, goto Step  2.
Step  2. The corrector step.
Set i=0.
Repeat
If i=5, set IT=IT+5, wk=wk-1, tk=tk-1, dk+1=dk, k=k+1, the corrector step fails, goto Step  1.2.
Compute the Newton step d(k,i+1) by solving
(H′(w(k,i),t(k,i))(dk)T)d(k,i+1)=(-H(w(k,i),t(k,i))0),
the corrector point
(w(k,i+1),t(k,i+1))=(w(k,i),t(k,i))+d(k,i+1);
set i=i+1.
If w(k,i)∉Ωθ(t(k,i))0×R++ or t(k,i)≥1,
set IT=IT+i, wk=wk-1, tk=tk-1, dk+1=dk, k=k+1, the corrector step fails, goto Step  1.2;
else if t(k,i)<0,
compute the point (w(k,0),0) by a damping Newton step, set IT=IT+i, if w(k,0)∈Ω×R+, goto Step  3,
else, set wk=wk-1, tk=tk-1, dk+1=dk, k=k+1, the corrector step fails, goto Step  1.2.
Until ∥H(w(k,i),t(k,i))∥inf⁡≤Htol and ∥d(k,i)∥≤Htol.
If tk<tc, return with x*=xk, stop the algorithm; else, set Htol=min⁡{Htol,tk}, k=k+1, goto Step  1.
Step  3. The end game.
Set l=0.
Repeat
If l=5 or ∥d(k,l)∥>∥d(k,l-1)∥ for l>1, or w(k,l)∉Ω×R+, set tend=0.3*tend, wk=wk-1, tk=tk-1,
dk+1=dk, IT=IT+l, k=k+1, the corrector step fails, goto Step  1.2.
Compute the Newton step d(k,l+1) by solving Ftc′(w(k,l))d(k,l+1)=-Ftc(w(k,l)), the corrector point
w(k,l+1)=w(k,l)+d(k,l+1);
set l=l+1.
Until ∥Ftc(w(k,l))∥inf⁡≤Hfinal and ∥d(k,l)∥≤Hfinal.
Set IT=IT+l, return with x*=x(k,l), stop the algorithm.


### 4.2. The Numerical Experiment

Although there exist so many test problems, such as the CUTEr test set [21], we cannot find a large collection of test problems with moderate variables and very many complicated nonlinear constraints. In this paper, six test problems are chosen to test the algorithm. Problem 4.1 is chosen from the CUTEr test set and it is used to illustrate a special situation; others are derived from the discretized semi-infinite programming problems. We also give two artificial test problems 4.2 and 4.3 and use three problems 4.4–4.6 in [22]. The numbers of variables in problem 4.2 and the number of constrains in problems 4.2–4.6 can be arbitrary. For each test problem, the gradients and Hessians of the objective and constraint functions are evaluated analytically.The FACH-S-N procedure was implemented in MATLAB. To illustrate its efficiency, we also implemented the ACH method in [19] using a similar procedure. In addition, we downloaded the state-of-the-art solver KNITRO, which provides three algorithms for solving large-scale nonlinear programming problems, and we used the interior-point direct method with default parameters to compare with the FACH and MFACH methods. For any iterate point xk, if gmax⁡(xk)<10-6, it was treated as a feasible point. The test results were obtained by running MATLAB R2008a on a desktop with Windows XP Professional operation system, Intel(R) Core(TM) i5-750 2.66 GHz processor, and 8 GB of memory. The default parameters were chosen as follows.(i)
Parameters for the flattened aggregate constraint function:θ=0.01, c1=0.05, c2=0.5*10-5, and α=2.(ii)
Parameters in the end game section:tc=10-6 and tend=0.1.(iii)
Step size parameters:h0=0.1, Bmin⁡=[0.5,0.75], and Bmax⁡=[3,1.5].(iv)
Tracking tolerances:Htol=10-5 and Hfinal=10-12.(v)
Initial Lagrangian multipliers:1 for ACH, FACH, and MFACH methods.For each problem with different parameters, we list the value of objective functionf(x) and the maximal function of constraints gmax⁡(x) at x*, the number of iterations IT, the number of evaluations of gradients of individual constraint functions in the FACH and MFACH methods Ns (in contrast, Ns is m*IT in the ACH method), and the running time in seconds Time. For problems that were not solved by the conservative setting, we also give the reason for failure. The notation “fail1” indicates that the steplength in predictor step is smaller than 10-10; it is generally due to poor conditioned Jacobian matrix. The notation “fail2” means out of memory. The notation “fail3” means no result in 5000 iterations or 3600 seconds.Example 14 (see [21]).
Consider(55)f(x)=sin⁡(x1-1+1.5π)+∑2≤i≤1000100sin⁡(-xi+1.5π+xi-12),gi(x)={x1-πi=1;xi-12-xi-πi=2,…,1000;-x1-πi=1001;-xi-12+xi-πi=1002,…,2000,x0=(1,…,1)∈R1000.Remark 15.
In this example, the starting point happens to be an unconstrained minimizer of the objective function and we found that thex-components of iterative points generated by the FACH and MFACH methods (not other homotopy methods) remain invariant. This is not an occasional phenomenon. In fact, if (x0,0)∈Ω0×R+m is a solution of the KKT system, that is, x0 is a stationary point of the objective function f(x), when the parameters α, c1, and c2 in FACH and MFACH methods satisfy gmax⁡(x0)≤-α(c1+c2), the x-components of points on the homotopy path Γw0 will remain invariant, which can be derived from the fact that (x0,tλ0ɛ(1)/ɛ(t),t) is the solution of (39) and (50) for t∈(0,1]. Moreover, since
(56)H′(x0,λ0,1)=(I000g^(x0,1)*),
the x-component of d1 is 0. By a similar discussion, the x-component of d(1,i) is 0, and hence x1=x(1,i)=x0. In this analogy, we know that xk=x(k,i)=x0 for any k and i in algorithm FACH-S-N.Example 16.
Consider(57)f(x)=x32+x42,gi,j(x)=(ti-x1)2x32+(tj′-x2)2x42-1,ti=i(m-1),i=0,…,m-1,tj′=j(m-1),j=0,…,m-1,x0=(0,0,100,100)∈R4.Example 17.
Consider(58)f(x)=1n∑1≤k≤n(xk-1)2,gi(x)=∏1≤k≤ncos⁡⁡(tixk)+ti∑1≤k≤nxk3,ti=0.5+πi(m-1),i=0,…,m-1,x0=(-2,…,-2)∈Rn.Example 18 (see [22]).
Consider(59)f(x)=x12+x22+x32,gi(x)=x1+x2exp⁡⁡(x3ti)+exp⁡⁡(2ti)-2sin⁡(4ti),ti=i(m-1),i=0,…,m-1,x0=(-200,-200,200)∈R3.Example 19 (see [22]).
Consider(60)f(x)=(x1-2x2+5x22-x23-13)2+(x1-14x2+x22+x23-29)2,gi(x)=x12+2x2ti2+exp⁡⁡(x1+x2)-exp⁡⁡(ti),ti=i(m-1),i=0,…,m-1,x0=(0,-45)∈R2.Example 20 (see [22]).
Consider(61)f(x)=x123+x12+x22,gi(x)=(1-x12ti2)2-x1ti2-x22+x2,ti=i(m-1),i=0,…,m-1,x0=(-1,100)∈R2.
To explain the numerical efficiency, we make the following remarks by numerical results in Tables 1, 2, 3, 4, 5 and 6.(i)
Ifgmax⁡(x)≤-αɛ(t) for any (x,λ,t)∈Γw0, the FACH and MFACH methods do not need to calculate the gradient and Hessian of any constraint functions.(ii)
For problems whose gradients and Hessians of constraint functions are expensive to evaluate, the performance of the FACH and MFACH methods is much better than the ACH method based on similar numerical tracing procedures.(iii)
Compared to the interior-point direct algorithm of the state-of-the-art solver KNITRO, the test results show that the FACH and MFACH methods perform worse whenm is small but much better when m is large for most problems. In addition, we can see that their time cost increases more slowly than KNITRO solver as m increases.(iv)
The functionφ(z,t) is important for the flattened aggregate constraint function. Theoretically, the parameters c1, c2, and α in the FACH and MFACH methods can be chosen freely. However, they do matter in the practical efficiency of the FACH and MFACH methods and should be suitably chosen. If these parameters are too large, then too many gradients and Hessians of individual constraint functions need to be evaluated and, hence, cause low efficiency. On the other hand, if they are too small, the Hessian of the flattened aggregated constraint function (10) may become ill-conditioned. In our numerical tests, we fixed c1=0.05, c2=0.5*10-5, and α=2. In addition, the function φ¯(z,t) can be defined in many ways, and preliminary numerical experiments show that the algorithms with different functions φ¯(z,t) have similar efficiencies.(v)
In algorithm FACH-S-N, we gave only a simple implementation of the ACH, FACH, and MFACH methods. To improve implementation of the FACH and MFACH methods, a lot of work needs to be done on all processes of numerical path tracing, say, schemes of predictor and corrector, steplength updating, linear system solving, and end game. Other practical strategies in the literature for large-scale nonlinear programming problems (e.g., [9, 23–26]) are also very important for improving the efficiency.Table 1
Test results for Example14.
Method
f
(
x
*
)
g
max
⁡
(
x
*
)
IT
N
s
Time
ACH
−99901.0000
−2.1416
35
—
85.5496
KNITRO
−99901.0000
−2.1416
0
—
2.8733
FACH
−99901.0000
−2.1416
4
0
8.3912
MFACH
−99901.0000
−2.1416
4
0
8.4150Table 2
Test results for Example16.
m
Method
f
(
x
*
)
g
max
⁡
(
x
*
)
IT
N
s
Time
1
0
2
ACH
1.0000
−1.3863e −08
394
—
0.2423
KNITRO
1.0000
−7.1233e −10
18
—
0.2288
FACH
1.0000
−1.3863e −08
541
721
0.3270
MFACH
1.0000
−1.3863e −08
567
823
0.2692
1
0
4
ACH
1.0000
−1.3863e −08
397
—
1.9730
KNITRO
1.0000
−8.0005e −10
20
—
0.7227
FACH
1.0000
−1.3863e −08
541
721
0.4359
MFACH
1.0000
−1.3863e −08
567
823
0.4676
1
0
6
ACH
1.0000
−1.3863e −08
404
—
249.3395
KNITRO
1.0000
−7.8830e −10
13
—
2170.4301
FACH
1.0000
−1.3863e −08
541
805
10.9486
MFACH
1.0000
−1.3863e −08
567
907
11.2815Table 3
(a) Test results for Example17 with n=100. (b) Test results for Example 17 with n=500. (c) Test results for Example 17 with n=1000. (d) Test results for Example 17 with n=2000.
(a)
m
Method
f
(
x
*
)
g
max
⁡
(
x
*
)
IT
N
s
Time
1
0
2
ACH
1.4918
−1.5282e −03
243
—
36.9898
KNITRO
0.2603
−1.9708e −06
32
—
4.1262
FACH
1.4915
−6.6613e −16
156
96
0.4067
MFACH
1.4915
6.1062e −15
144
85
0.3186
1
0
3
ACH
1.4918
−1.5282e −03
243
—
325.6083
KNITRO
0.2603
−2.7042e −05
19
—
25.1720
FACH
1.4915
−6.6613e −16
156
96
0.6852
MFACH
1.4915
6.1062e −15
144
85
0.5342
1
0
4
ACH
1.4918
−1.5282e −03
243
—
3218.5422
KNITRO
0.2603
−2.5434e −04
18
—
271.1651
FACH
1.4915
−6.6614e −16
156
96
5.1941
MFACH
1.4915
6.1062e −15
144
86
4.0120
(b)
m
Method
f
(
x
*
)
g
max
⁡
(
x
*
)
IT
N
s
Time
1
0
2
ACH
1.2524
−1.0693e −02
405
—
1371.0198
KNITRO
0.1417
−2.0927e −05
31
—
95.1782
FACH
1.2510
1.2490e −14
159
45
3.9951
MFACH
1.2510
1.2490e −14
139
25
2.9676
1
0
3
ACH
—
—
—
—
fail3
KNITRO
—
—
—
—
fail3
FACH
1.2510
1.2490e −14
159
45
6.5173
MFACH
1.2510
1.2490e −14
139
25
4.6898
1
0
4
ACH
—
—
—
—
fail2
KNITRO
—
—
—
—
fail2
FACH
1.2510
1.2490e −14
159
45
31.7107
MFACH
1.2510
1.2490e −14
139
25
25.1842
(c)
m
Method
f
(
x
*
)
g
max
⁡
(
x
*
)
IT
N
s
Time
1
0
2
ACH
—
—
—
—
fail3
KNITRO
0.1099
−5.6905e −05
27
—
359.4507
FACH
1.1880
2.3370e −14
197
46
20.8919
MFACH
1.1880
2.4258e −14
197
46
21.2641
1
0
3
ACH
—
—
—
—
fail2
KNITRO
—
—
—
—
fail2
FACH
1.1880
2.3370e −14
197
46
28.8674
MFACH
1.1880
2.4258e −14
197
46
27.5140
1
0
4
ACH
—
—
—
—
fail2
KNITRO
—
—
—
—
fail2
FACH
1.1880
2.3370e −14
197
46
101.3199
MFACH
1.1880
2.4258e −14
197
46
99.1354
(d)
m
Method
f
(
x
*
)
g
max
⁡
(
x
*
)
IT
N
s
Time
1
0
2
ACH
—
—
—
—
fail3
KNITRO
—
—
—
—
fail3
FACH
1.1406
−6.2839e −14
260
47
126.6290
MFACH
1.1406
4.0967e −14
258
46
125.7973
1
0
3
ACH
—
—
—
—
fail2
KNITRO
—
—
—
—
fail2
FACH
1.1406
−6.2839e −14
260
47
146.6079
MFACH
1.1406
4.0967e −14
258
46
145.7293
1
0
4
ACH
—
—
—
—
fail2
KNITRO
—
—
—
—
fail2
FACH
1.1406
−6.2839e −14
260
47
368.5331
MFACH
1.1406
4.0967e −14
258
46
369.3701Table 4
Test results for Example18.
m
Method
f
(
x
*
)
g
max
⁡
(
x
*
)
IT
N
s
Time
1
0
2
ACH
5.3347
−4.4409e −16
1605
—
0.7546
KNITRO
1.2000e +05
−3.9900e +02
2
—
0.0141
FACH
5.3347
1.3323e −15
815
288
0.2977
MFACH
5.3347
1.3323e −15
820
293
0.3004
1
0
4
ACH
5.3347
−4.4409e −16
1605
—
7.0581
KNITRO
1.2000e +05
−3.9900e +02
2
—
0.3163
FACH
5.3347
1.3323e −15
815
288
1.3127
MFACH
5.3347
1.3323e −15
820
294
1.3185
1
0
6
ACH
5.3347
1.3323e −15
1606
—
575.4663
KNITRO
1.2000e +05
−3.9900e +02
2
—
3243.2201
FACH
5.3347
−9.2480e −08
960
6444
42.7664
MFACH
5.3347
−1.1278e −07
931
7094
41.7933Table 5
Test results for Example19.
m
Method
f
(
x
*
)
g
max
⁡
(
x
*
)
IT
N
s
Time
1
0
2
ACH
97.1589
0.0000
213
—
0.0964
KNITRO
97.1592
−7.5266 − 05
68
—
0.5162
FACH
97.1589
0.0000
170
191
0.0795
MFACH
97.1589
0.0000
160
183
0.0608
1
0
4
ACH
97.1589
0.0000
227
—
0.5310
KNITRO
102.2095
−6.7465 − 01
263
—
15.4868
FACH
97.1589
0.0000
249
10926
0.2846
MFACH
97.1589
0.0000
251
11321
0.2783
1
0
6
ACH
97.1589
0.0000
229
—
50.8763
KNITRO
188.9496
−8.3841e − 01
17
—
2206.4660
FACH
97.1589
0.0000
334
427520
9.5724
MFACH
97.1589
0.0000
253
486506
6.8910Table 6
Test results for Example20.
m
Method
f
(
x
*
)
g
max
⁡
(
x
*
)
IT
N
s
Time
1
0
2
ACH
2.4305
0.0000
342
—
0.1749
KNITRO
2.4305
−1.3822e −10
21
—
0.0291
FACH
2.4305
2.2204e −16
348
387
0.1368
MFACH
2.4305
2.2204e −16
341
374
0.1276
1
0
4
ACH
2.4305
−6.6632e −09
344
—
0.9804
KNITRO
2.4305
−1.6180e −10
38
—
1.0329
FACH
2.4305
−6.6632e −09
296
14710
0.2610
MFACH
2.4305
−6.6632e −09
358
23559
0.2802
1
0
6
ACH
2.4305
−4.9783e −08
344
—
94.3943
KNITRO
2.6182
−1.2609e −04
61
—
2372.2093
FACH
—
—
—
—
fail1
MFACH
2.4305
−4.9021e −08
516
2231308
16.3756


## 4.1. The FACH-S-N Procedure

In this section, we give a numerical procedure, FACH-S-N procedure, to trace the flattened aggregate constraint homotopy path by secant predictor and Newton corrector steps. It consists of three main steps: the predictor step, the corrector step, and the end game.The predictor step is an approximate step along the homotopy path: it uses a predictor directiondk and a steplength hk to get a predictor point. The first predictor direction uses the tangent direction, and others use the secant direction. The steplength is determined by several parameters. It is set to no more than 1, which can ensure that the predictor point is close to the homotopy path and hence stays in the convergence domain of Newton’s method in the corrector step. If the angle of the predictor direction and the previous one βk is greater than π/4, the steplength will be decreased to avoid using the opposite direction as the predictor direction. If the corrector criteria are satisfied with no more than four Newton iterations for three times in succession, the steplength will be increased or kept invariable. Otherwise, the steplength will be decreased.Once a predictor point(w(k,0),t(k,0)) is calculated, one or more Newton iterations are used to bring the predictor point back to the homotopy path in the corrector step. The corrector points (w(k,i),t(k,i)), i=1,…,5, are calculated by (w(k,i+1),t(k,i+1))=(w(k,i),t(k,i))+d(k,i+1), i=0,…,4, where the step d(k,i+1) is the solution of an augmented system defined by the homotopy equation and the direction perpendicular to the predictor direction. The corrector step terminates when d(k,i) and H(w(k,i),t(k,i)) satisfy the tolerance criterions.At each predictor step and corrector step, the feasibility of the predictor point(wk,tk) and the corrector point (w(k,i),t(k,i)) needs to be checked. If tk<tc in the predictor step or tk<0 in the corrector step, a damping step is used to get a new point (w(k,0),0). Then, if w(k,0) is feasible, the end game, a more efficient strategy than predictor-corrector steps when the homotopy parameter is close to 0, is invoked. Starting with w(k,0), Newton's method is used to solve
(54)Ftc(x,λ)=(∇f(x)+λh(x,tc)λg^(x,tc))=0,
where tc is a small positive constant. For other situations, the steplength will be decreased to make new predictor-corrector steps (see Algorithm 1).Algorithm 1:The FACH-S-N procedure.
Input  Give θ∈(0,1], c1>0, c2>0, α>1, tend and tc; starting point w0∈Ω^×R++; initial steplength h1, steplength
contraction factor Bmin⁡, steplength expansion factor Bmax⁡; tracking tolerances Htol and Hfinal for correction. Set
t0=1, k=1, IT=0, Ngood=2, β1=0.
Step  1. The predictor step.
Step  1.1. Compute the predictor direction.
If k=1, compute d by solving
(H′(x0,λ0,1)(d0)T)d=-d0,
where d0=(0,…,0,-1)∈Rn+2, set the predictor direction d1=d/∥d∥;
else, compute the predictor direction dk, the angle βk between dk and dk-1:
dk=(wk-1,tk-1)-(wk-2,tk-2)∥(wk-1,tk-1)-(wk-2,tk-2)∥,βk=arccos((dk)Tdk-1).
Step  1.2. Adjust the steplength.
Adjust the steplength hk as follows:
if the predictor point is infeasible, set hk=Bmin⁡(1)hk, Ngood=0;
else if the corrector step fails or βk>π/4, set hk=Bmin⁡(1)hk-1, Ngood=0;
else if i=5, set hk=Bmin⁡(2)hk-1, Ngood=0;
else if i=4, set hk=hk-1, Ngood=Ngood+1;
else if i=3, set Ngood=Ngood+1,
if Ngood>2, set hk=min⁡{1,Bmax⁡(2)hk-1}, else, set hk=hk-1;
else, set Ngood=Ngood+1,
if Ngood>2, set hk=min⁡{1,Bmax⁡(1)hk-1}, else, set hk=hk-1.
If hk<10-10, stop the algorithm with an error flag.
Step  1.3. Compute the predictor point and check its feasibility.
Compute the predictor point
(w(k,0),t(k,0))=(wk-1,tk-1)+hkdk.
If w(k,0)∉Ωθ(t(k,0))0×R++ or t(k,0)≥1,
the predictor point is infeasible, goto Step  1.2;
else if t(k,0)≤tend,
adjust the steplength hk, compute the point (w(k,0),0), if w(k,0)∈Ω×R+, goto Step  3, else, the predictor
point is infeasible, goto Step  1.2;
else, goto Step  2.
Step  2. The corrector step.
Set i=0.
Repeat
If i=5, set IT=IT+5, wk=wk-1, tk=tk-1, dk+1=dk, k=k+1, the corrector step fails, goto Step  1.2.
Compute the Newton step d(k,i+1) by solving
(H′(w(k,i),t(k,i))(dk)T)d(k,i+1)=(-H(w(k,i),t(k,i))0),
the corrector point
(w(k,i+1),t(k,i+1))=(w(k,i),t(k,i))+d(k,i+1);
set i=i+1.
If w(k,i)∉Ωθ(t(k,i))0×R++ or t(k,i)≥1,
set IT=IT+i, wk=wk-1, tk=tk-1, dk+1=dk, k=k+1, the corrector step fails, goto Step  1.2;
else if t(k,i)<0,
compute the point (w(k,0),0) by a damping Newton step, set IT=IT+i, if w(k,0)∈Ω×R+, goto Step  3,
else, set wk=wk-1, tk=tk-1, dk+1=dk, k=k+1, the corrector step fails, goto Step  1.2.
Until ∥H(w(k,i),t(k,i))∥inf⁡≤Htol and ∥d(k,i)∥≤Htol.
If tk<tc, return with x*=xk, stop the algorithm; else, set Htol=min⁡{Htol,tk}, k=k+1, goto Step  1.
Step  3. The end game.
Set l=0.
Repeat
If l=5 or ∥d(k,l)∥>∥d(k,l-1)∥ for l>1, or w(k,l)∉Ω×R+, set tend=0.3*tend, wk=wk-1, tk=tk-1,
dk+1=dk, IT=IT+l, k=k+1, the corrector step fails, goto Step  1.2.
Compute the Newton step d(k,l+1) by solving Ftc′(w(k,l))d(k,l+1)=-Ftc(w(k,l)), the corrector point
w(k,l+1)=w(k,l)+d(k,l+1);
set l=l+1.
Until ∥Ftc(w(k,l))∥inf⁡≤Hfinal and ∥d(k,l)∥≤Hfinal.
Set IT=IT+l, return with x*=x(k,l), stop the algorithm.


## 4.2. The Numerical Experiment

Although there exist so many test problems, such as the CUTEr test set [21], we cannot find a large collection of test problems with moderate variables and very many complicated nonlinear constraints. In this paper, six test problems are chosen to test the algorithm. Problem 4.1 is chosen from the CUTEr test set and it is used to illustrate a special situation; others are derived from the discretized semi-infinite programming problems. We also give two artificial test problems 4.2 and 4.3 and use three problems 4.4–4.6 in [22]. The numbers of variables in problem 4.2 and the number of constrains in problems 4.2–4.6 can be arbitrary. For each test problem, the gradients and Hessians of the objective and constraint functions are evaluated analytically.The FACH-S-N procedure was implemented in MATLAB. To illustrate its efficiency, we also implemented the ACH method in [19] using a similar procedure. In addition, we downloaded the state-of-the-art solver KNITRO, which provides three algorithms for solving large-scale nonlinear programming problems, and we used the interior-point direct method with default parameters to compare with the FACH and MFACH methods. For any iterate point xk, if gmax⁡(xk)<10-6, it was treated as a feasible point. The test results were obtained by running MATLAB R2008a on a desktop with Windows XP Professional operation system, Intel(R) Core(TM) i5-750 2.66 GHz processor, and 8 GB of memory. The default parameters were chosen as follows.(i)
Parameters for the flattened aggregate constraint function:θ=0.01, c1=0.05, c2=0.5*10-5, and α=2.(ii)
Parameters in the end game section:tc=10-6 and tend=0.1.(iii)
Step size parameters:h0=0.1, Bmin⁡=[0.5,0.75], and Bmax⁡=[3,1.5].(iv)
Tracking tolerances:Htol=10-5 and Hfinal=10-12.(v)
Initial Lagrangian multipliers:1 for ACH, FACH, and MFACH methods.For each problem with different parameters, we list the value of objective functionf(x) and the maximal function of constraints gmax⁡(x) at x*, the number of iterations IT, the number of evaluations of gradients of individual constraint functions in the FACH and MFACH methods Ns (in contrast, Ns is m*IT in the ACH method), and the running time in seconds Time. For problems that were not solved by the conservative setting, we also give the reason for failure. The notation “fail1” indicates that the steplength in predictor step is smaller than 10-10; it is generally due to poor conditioned Jacobian matrix. The notation “fail2” means out of memory. The notation “fail3” means no result in 5000 iterations or 3600 seconds.Example 14 (see [21]).
Consider(55)f(x)=sin⁡(x1-1+1.5π)+∑2≤i≤1000100sin⁡(-xi+1.5π+xi-12),gi(x)={x1-πi=1;xi-12-xi-πi=2,…,1000;-x1-πi=1001;-xi-12+xi-πi=1002,…,2000,x0=(1,…,1)∈R1000.Remark 15.
In this example, the starting point happens to be an unconstrained minimizer of the objective function and we found that thex-components of iterative points generated by the FACH and MFACH methods (not other homotopy methods) remain invariant. This is not an occasional phenomenon. In fact, if (x0,0)∈Ω0×R+m is a solution of the KKT system, that is, x0 is a stationary point of the objective function f(x), when the parameters α, c1, and c2 in FACH and MFACH methods satisfy gmax⁡(x0)≤-α(c1+c2), the x-components of points on the homotopy path Γw0 will remain invariant, which can be derived from the fact that (x0,tλ0ɛ(1)/ɛ(t),t) is the solution of (39) and (50) for t∈(0,1]. Moreover, since
(56)H′(x0,λ0,1)=(I000g^(x0,1)*),
the x-component of d1 is 0. By a similar discussion, the x-component of d(1,i) is 0, and hence x1=x(1,i)=x0. In this analogy, we know that xk=x(k,i)=x0 for any k and i in algorithm FACH-S-N.Example 16.
Consider(57)f(x)=x32+x42,gi,j(x)=(ti-x1)2x32+(tj′-x2)2x42-1,ti=i(m-1),i=0,…,m-1,tj′=j(m-1),j=0,…,m-1,x0=(0,0,100,100)∈R4.Example 17.
Consider(58)f(x)=1n∑1≤k≤n(xk-1)2,gi(x)=∏1≤k≤ncos⁡⁡(tixk)+ti∑1≤k≤nxk3,ti=0.5+πi(m-1),i=0,…,m-1,x0=(-2,…,-2)∈Rn.Example 18 (see [22]).
Consider(59)f(x)=x12+x22+x32,gi(x)=x1+x2exp⁡⁡(x3ti)+exp⁡⁡(2ti)-2sin⁡(4ti),ti=i(m-1),i=0,…,m-1,x0=(-200,-200,200)∈R3.Example 19 (see [22]).
Consider(60)f(x)=(x1-2x2+5x22-x23-13)2+(x1-14x2+x22+x23-29)2,gi(x)=x12+2x2ti2+exp⁡⁡(x1+x2)-exp⁡⁡(ti),ti=i(m-1),i=0,…,m-1,x0=(0,-45)∈R2.Example 20 (see [22]).
Consider(61)f(x)=x123+x12+x22,gi(x)=(1-x12ti2)2-x1ti2-x22+x2,ti=i(m-1),i=0,…,m-1,x0=(-1,100)∈R2.
To explain the numerical efficiency, we make the following remarks by numerical results in Tables 1, 2, 3, 4, 5 and 6.(i)
Ifgmax⁡(x)≤-αɛ(t) for any (x,λ,t)∈Γw0, the FACH and MFACH methods do not need to calculate the gradient and Hessian of any constraint functions.(ii)
For problems whose gradients and Hessians of constraint functions are expensive to evaluate, the performance of the FACH and MFACH methods is much better than the ACH method based on similar numerical tracing procedures.(iii)
Compared to the interior-point direct algorithm of the state-of-the-art solver KNITRO, the test results show that the FACH and MFACH methods perform worse whenm is small but much better when m is large for most problems. In addition, we can see that their time cost increases more slowly than KNITRO solver as m increases.(iv)
The functionφ(z,t) is important for the flattened aggregate constraint function. Theoretically, the parameters c1, c2, and α in the FACH and MFACH methods can be chosen freely. However, they do matter in the practical efficiency of the FACH and MFACH methods and should be suitably chosen. If these parameters are too large, then too many gradients and Hessians of individual constraint functions need to be evaluated and, hence, cause low efficiency. On the other hand, if they are too small, the Hessian of the flattened aggregated constraint function (10) may become ill-conditioned. In our numerical tests, we fixed c1=0.05, c2=0.5*10-5, and α=2. In addition, the function φ¯(z,t) can be defined in many ways, and preliminary numerical experiments show that the algorithms with different functions φ¯(z,t) have similar efficiencies.(v)
In algorithm FACH-S-N, we gave only a simple implementation of the ACH, FACH, and MFACH methods. To improve implementation of the FACH and MFACH methods, a lot of work needs to be done on all processes of numerical path tracing, say, schemes of predictor and corrector, steplength updating, linear system solving, and end game. Other practical strategies in the literature for large-scale nonlinear programming problems (e.g., [9, 23–26]) are also very important for improving the efficiency.Table 1
Test results for Example14.
Method
f
(
x
*
)
g
max
⁡
(
x
*
)
IT
N
s
Time
ACH
−99901.0000
−2.1416
35
—
85.5496
KNITRO
−99901.0000
−2.1416
0
—
2.8733
FACH
−99901.0000
−2.1416
4
0
8.3912
MFACH
−99901.0000
−2.1416
4
0
8.4150Table 2
Test results for Example16.
m
Method
f
(
x
*
)
g
max
⁡
(
x
*
)
IT
N
s
Time
1
0
2
ACH
1.0000
−1.3863e −08
394
—
0.2423
KNITRO
1.0000
−7.1233e −10
18
—
0.2288
FACH
1.0000
−1.3863e −08
541
721
0.3270
MFACH
1.0000
−1.3863e −08
567
823
0.2692
1
0
4
ACH
1.0000
−1.3863e −08
397
—
1.9730
KNITRO
1.0000
−8.0005e −10
20
—
0.7227
FACH
1.0000
−1.3863e −08
541
721
0.4359
MFACH
1.0000
−1.3863e −08
567
823
0.4676
1
0
6
ACH
1.0000
−1.3863e −08
404
—
249.3395
KNITRO
1.0000
−7.8830e −10
13
—
2170.4301
FACH
1.0000
−1.3863e −08
541
805
10.9486
MFACH
1.0000
−1.3863e −08
567
907
11.2815Table 3
(a) Test results for Example17 with n=100. (b) Test results for Example 17 with n=500. (c) Test results for Example 17 with n=1000. (d) Test results for Example 17 with n=2000.
(a)
m
Method
f
(
x
*
)
g
max
⁡
(
x
*
)
IT
N
s
Time
1
0
2
ACH
1.4918
−1.5282e −03
243
—
36.9898
KNITRO
0.2603
−1.9708e −06
32
—
4.1262
FACH
1.4915
−6.6613e −16
156
96
0.4067
MFACH
1.4915
6.1062e −15
144
85
0.3186
1
0
3
ACH
1.4918
−1.5282e −03
243
—
325.6083
KNITRO
0.2603
−2.7042e −05
19
—
25.1720
FACH
1.4915
−6.6613e −16
156
96
0.6852
MFACH
1.4915
6.1062e −15
144
85
0.5342
1
0
4
ACH
1.4918
−1.5282e −03
243
—
3218.5422
KNITRO
0.2603
−2.5434e −04
18
—
271.1651
FACH
1.4915
−6.6614e −16
156
96
5.1941
MFACH
1.4915
6.1062e −15
144
86
4.0120
(b)
m
Method
f
(
x
*
)
g
max
⁡
(
x
*
)
IT
N
s
Time
1
0
2
ACH
1.2524
−1.0693e −02
405
—
1371.0198
KNITRO
0.1417
−2.0927e −05
31
—
95.1782
FACH
1.2510
1.2490e −14
159
45
3.9951
MFACH
1.2510
1.2490e −14
139
25
2.9676
1
0
3
ACH
—
—
—
—
fail3
KNITRO
—
—
—
—
fail3
FACH
1.2510
1.2490e −14
159
45
6.5173
MFACH
1.2510
1.2490e −14
139
25
4.6898
1
0
4
ACH
—
—
—
—
fail2
KNITRO
—
—
—
—
fail2
FACH
1.2510
1.2490e −14
159
45
31.7107
MFACH
1.2510
1.2490e −14
139
25
25.1842
(c)
m
Method
f
(
x
*
)
g
max
⁡
(
x
*
)
IT
N
s
Time
1
0
2
ACH
—
—
—
—
fail3
KNITRO
0.1099
−5.6905e −05
27
—
359.4507
FACH
1.1880
2.3370e −14
197
46
20.8919
MFACH
1.1880
2.4258e −14
197
46
21.2641
1
0
3
ACH
—
—
—
—
fail2
KNITRO
—
—
—
—
fail2
FACH
1.1880
2.3370e −14
197
46
28.8674
MFACH
1.1880
2.4258e −14
197
46
27.5140
1
0
4
ACH
—
—
—
—
fail2
KNITRO
—
—
—
—
fail2
FACH
1.1880
2.3370e −14
197
46
101.3199
MFACH
1.1880
2.4258e −14
197
46
99.1354
(d)
m
Method
f
(
x
*
)
g
max
⁡
(
x
*
)
IT
N
s
Time
1
0
2
ACH
—
—
—
—
fail3
KNITRO
—
—
—
—
fail3
FACH
1.1406
−6.2839e −14
260
47
126.6290
MFACH
1.1406
4.0967e −14
258
46
125.7973
1
0
3
ACH
—
—
—
—
fail2
KNITRO
—
—
—
—
fail2
FACH
1.1406
−6.2839e −14
260
47
146.6079
MFACH
1.1406
4.0967e −14
258
46
145.7293
1
0
4
ACH
—
—
—
—
fail2
KNITRO
—
—
—
—
fail2
FACH
1.1406
−6.2839e −14
260
47
368.5331
MFACH
1.1406
4.0967e −14
258
46
369.3701Table 4
Test results for Example18.
m
Method
f
(
x
*
)
g
max
⁡
(
x
*
)
IT
N
s
Time
1
0
2
ACH
5.3347
−4.4409e −16
1605
—
0.7546
KNITRO
1.2000e +05
−3.9900e +02
2
—
0.0141
FACH
5.3347
1.3323e −15
815
288
0.2977
MFACH
5.3347
1.3323e −15
820
293
0.3004
1
0
4
ACH
5.3347
−4.4409e −16
1605
—
7.0581
KNITRO
1.2000e +05
−3.9900e +02
2
—
0.3163
FACH
5.3347
1.3323e −15
815
288
1.3127
MFACH
5.3347
1.3323e −15
820
294
1.3185
1
0
6
ACH
5.3347
1.3323e −15
1606
—
575.4663
KNITRO
1.2000e +05
−3.9900e +02
2
—
3243.2201
FACH
5.3347
−9.2480e −08
960
6444
42.7664
MFACH
5.3347
−1.1278e −07
931
7094
41.7933Table 5
Test results for Example19.
m
Method
f
(
x
*
)
g
max
⁡
(
x
*
)
IT
N
s
Time
1
0
2
ACH
97.1589
0.0000
213
—
0.0964
KNITRO
97.1592
−7.5266 − 05
68
—
0.5162
FACH
97.1589
0.0000
170
191
0.0795
MFACH
97.1589
0.0000
160
183
0.0608
1
0
4
ACH
97.1589
0.0000
227
—
0.5310
KNITRO
102.2095
−6.7465 − 01
263
—
15.4868
FACH
97.1589
0.0000
249
10926
0.2846
MFACH
97.1589
0.0000
251
11321
0.2783
1
0
6
ACH
97.1589
0.0000
229
—
50.8763
KNITRO
188.9496
−8.3841e − 01
17
—
2206.4660
FACH
97.1589
0.0000
334
427520
9.5724
MFACH
97.1589
0.0000
253
486506
6.8910Table 6
Test results for Example20.
m
Method
f
(
x
*
)
g
max
⁡
(
x
*
)
IT
N
s
Time
1
0
2
ACH
2.4305
0.0000
342
—
0.1749
KNITRO
2.4305
−1.3822e −10
21
—
0.0291
FACH
2.4305
2.2204e −16
348
387
0.1368
MFACH
2.4305
2.2204e −16
341
374
0.1276
1
0
4
ACH
2.4305
−6.6632e −09
344
—
0.9804
KNITRO
2.4305
−1.6180e −10
38
—
1.0329
FACH
2.4305
−6.6632e −09
296
14710
0.2610
MFACH
2.4305
−6.6632e −09
358
23559
0.2802
1
0
6
ACH
2.4305
−4.9783e −08
344
—
94.3943
KNITRO
2.6182
−1.2609e −04
61
—
2372.2093
FACH
—
—
—
—
fail1
MFACH
2.4305
−4.9021e −08
516
2231308
16.3756


## 5. Conclusions

By introducing a flattened aggregate constraint function, a flattened aggregate constraint homotopy method is proposed for nonlinear programming problems with few variables and many nonlinear constraints, and its global convergence is proven. By greatly reducing the computation of gradients and Hessians of constraint functions in each iteration, the proposed method is very competitive for nonlinear programming problems with a large number of complicated constraint functions.



---

*Source: 430932-2014-05-29.xml*