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

Abstract

The aggregate constraint homotopy method uses a single smoothing constraint instead of 𝑚-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.

1. Introduction

In this paper, we consider the following nonlinear programming problem: (𝑥),(𝑥)min𝑓s.t.𝑔≤0, (1) where 𝑥∈𝑅^𝑛 is the variable, 𝑔(𝑥)=(𝑔₁(𝑥),…,𝑔_𝑚(𝑥))^𝑇, 𝑓∶𝑅^𝑛→𝑅, and 𝑔_𝑖(𝑥)∶𝑅^𝑛→𝑅, 𝑖=1,…,𝑚, are three times continuously differentiable, and 𝑚 is very large, but 𝑛 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.

Let 𝑔_max(𝑥)=max_{1≤𝑖≤𝑚}{𝑔_𝑖(𝑥)}; then (1) is equivalent to min𝑓(𝑥),s.t.𝑔_max(𝑥)≤0, (2) which has only one, but nonsmooth, constraint. In [17], the following aggregate function was introduced, which is a smooth approximation of 𝑔_max(𝑥) with a smoothing parameter 𝑡>0 and induced from the max-entropy theory: 𝑔(𝑥,𝑡)=𝑡ln^𝑚_𝑖=1𝑔exp^𝑖(𝑥)𝑡. (3) 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 ∇_𝑥𝑔(𝑥,𝑡)=^𝑚_𝑖=1𝑐_𝑖(𝑥,𝑡)∇𝑔_𝑖(𝑥),𝜕𝑔(𝑥,𝑡)=1𝜕𝑡𝑡𝑔(𝑥,𝑡)−^𝑚_𝑖=1𝑐_𝑖(𝑥,𝑡)𝑔_𝑖(𝑥),∇²_𝑥𝑥𝑔()=𝑥,𝑡_𝑚𝑖=1𝑐_𝑖(𝑥,𝑡)∇²𝑔_𝑖(𝑥)+1𝑡^𝑚_𝑖=1𝑐_𝑖(𝑥,𝑡)∇𝑔_𝑖(𝑥)∇𝑔_𝑖(𝑥)^𝑇−1𝑡^𝑚_𝑖=1𝑐_𝑖(𝑥,𝑡)∇𝑔_𝑖(𝑥)^𝑚_𝑖=1𝑐_𝑖(𝑥,𝑡)∇𝑔_𝑖(𝑥)^𝑇,∇²_𝑥𝑡𝑔()=1𝑥,𝑡𝑡^2𝑚_𝑖=1𝑐_𝑖(𝑥,𝑡)∇𝑔_𝑖(𝑥)^𝑚_𝑖=1𝑐_𝑖(𝑥,𝑡)𝑔_𝑖(𝑥)−^𝑚_𝑖=1𝑐_𝑖(𝑥,𝑡)𝑔_𝑖(𝑥)∇𝑔_𝑖(𝑥), (4) where 𝑐_𝑖()=𝑔𝑥,𝑡exp_𝑖(𝑥)/𝑡∑^𝑚_𝑖=1exp𝑔_𝑖(𝑥)∈(),/𝑡0,1^𝑚_𝑖=1𝑐_𝑖(𝑥,𝑡)=1, (5) are complicated combinations of gradients and Hessians of all constraint functions and, hence, are expensive to calculate when 𝑚 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 Section 2. 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 by 𝐼.

Unless otherwise specified, ‖⋅‖ denotes the Euclidean norm. Ω={𝑥∈𝑅^𝑛∣𝑔(𝑥)≤0} is the feasible set of (1), whereas Ω⁰={𝑥∈𝑅^𝑛∣𝑔(𝑥)<0} is the interior of Ω and 𝜕Ω=Ω⧵Ω⁰ is the boundary of Ω. The symbols 𝑅^𝑚₊ and 𝑅^𝑚₊₊ denote the nonnegative and positive quadrants of 𝑅^𝑚, respectively. The active index set is denoted by 𝐼(𝑥)={𝑖∈{1,…,𝑚}∣𝑔_𝑖(𝑥)=0} at 𝑥∈Ω. For a function 𝐹∶𝑅^𝑛→𝑅^𝑚, ∇𝐹(𝑥) is the 𝑛×𝑚 matrix whose (𝑖,𝑗)th element is 𝜕𝐹_𝑗(𝑥)/𝜕𝑥_𝑖; 𝐹⁻¹(𝑌)={𝑥∣𝐹(𝑥)∈𝑌} is the set-valued inverse for the set 𝑌⊆𝑅^𝑚.

Definition 1 (see [10]). The set Ω satisfies the normal cone condition if and only if the normal cone of Ω at any 𝑥∈𝜕Ω does not meet Ω⁰; that is, 𝑥+_{(𝑥)𝑖∈𝐼}𝑦_𝑖∇𝑔_𝑖(𝑥)∣𝑦_𝑖≥0∩Ω⁰=0/. (6)

Definition 2. If and only if there exists 𝜆^∗∈𝑅^𝑚₊ such that (𝑥^∗,𝜆^∗) satisfies ∇𝑓𝑥^∗+∇𝑔𝑥^∗𝜆^∗𝜆=0,_∗𝑖𝑔_𝑖𝑥^∗=0,𝜆^∗_𝑖≥0,𝑔_𝑖𝑥^∗≤0,𝑖=1,…,𝑚, (7) then 𝑥^∗ is called a KKT point of (1) and 𝜆^∗ is the corresponding Lagrangian multiplier.

Definition 3. Let 𝑈⊂𝑅^𝑛 be an open set and let 𝐹∶𝑈→𝑅^𝑚 be a differentiable map. We say 𝑦∈𝑅^𝑚 is a regular value of 𝐹 if and only if (𝑥)Range𝜕𝐹𝜕𝑥=𝑅^𝑚,∀𝑥∈𝐹⁻¹𝑦. (8)

Lemma 4 (parameterized Sard theorem [20]). Let 𝑈⊂𝑅^𝑛 and 𝑉⊂𝑅^𝑘 be two open sets and let 𝐹∶𝑈×𝑉→𝑅^𝑚 be a 𝐶^𝑟 differentiable map with 𝑟>max{0,𝑛−𝑚}. If 0∈𝑅^𝑚 is a regular value of 𝐹, then, for almost all 𝑎∈𝑉, 0 is a regular value of 𝐹_𝑎=𝐹(⋅,𝑎).

2. The Flattened Aggregate Constraint Function

In this section, we suppose that the following assumptions hold.

Assumption 5. Ω⁰ is nonempty and Ω is bounded.

Assumption 6. For any 𝑥∈𝜕Ω, {∇𝑔_𝑖(𝑥)∣𝑖∈𝐼(𝑥)} are positive independent; that is, _{(𝑥)𝑖∈𝐼}𝑦_𝑖∇𝑔_𝑖(𝑥)=0,𝑦_𝑖≥0,𝑖∈𝐼(𝑥)⟹𝑦_𝑖=0,𝑖∈𝐼(𝑥). (9)

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: 𝑔_{(𝜃,𝑐1,𝑐2,𝛼)}(𝑥,𝑡)=𝜃𝑡ln∑^𝑚_𝑖=1𝜑𝑔_𝑖(𝑥)𝑔,𝑡exp^𝑖(𝑥)−𝜀(𝑡)𝜃𝑡+exp,𝜃𝑡 (10) where 𝜀(𝑡)=𝑐₁𝑡+𝑐₂, (11)𝑐₁>0, 𝑐₂>0, 0<𝜃≤1 are some adjusting parameters, 𝜑(𝑧,𝑡)∶𝑅×[0,1]→[0,1] is a 𝐶³ bivariate function satisfying 𝜑(𝑧,𝑡)=0 for 𝑧≤−𝛼𝜀(𝑡) with 𝛼>1, and 𝜑(𝑧,𝑡)=1 for 𝑧≥−𝜀(𝑡). For simplicity of discussion, we write 𝑔_{(𝜃,𝑐1,𝑐2,𝛼)}(𝑥,𝑡) as 𝑔(𝑥,𝑡). By its definition, 𝜑(𝑧,𝑡) can be expressed as 𝜑()=(𝑡)𝑧,𝑡(𝑡),0,for𝑧≤−𝛼𝜀with𝛼>1,1,for𝑧≥−𝜀𝜑(𝑧,𝑡),otherwise, (12) where, for any 𝑡∈[0,1], 𝜑(𝑧,𝑡) satisfies 𝜑((𝑡))𝜕−𝛼𝜀,𝑡=0,𝜑(−𝛼𝜀(𝑡),𝑡)𝜕𝜕𝑧=0,₂𝜑(−𝛼𝜀(𝑡),𝑡)𝜕𝑧²𝜕=0,³𝜑(−𝛼𝜀(𝑡),𝑡)𝜕𝑧³=0,𝜑((𝑡))𝜕−𝜀,𝑡=1,𝜑(−𝜀(𝑡),𝑡)𝜕𝜕𝑧=0,₂𝜑(−𝜀(𝑡),𝑡)𝜕𝑧²𝜕=0,³𝜑(−𝜀(𝑡),𝑡)𝜕𝑧³=0. (13) There exist many kinds of functions 𝜑(𝑧,𝑡), such as polynomial functions and spline functions. Since polynomial functions have simple mathematical expressions and many excellent properties, throughout this paper, 𝜑(𝑧,𝑡) is chosen as 𝜑()=((𝑡𝑧,𝑡−20𝑧+𝜀))⁷(𝛼𝜀(𝑡)−𝜀(𝑡))⁷+−70(𝑧+𝜀(𝑡))⁶(𝛼𝜀(𝑡)−𝜀(𝑡))⁶+−84(𝑧+𝜀(𝑡))⁵(𝛼𝜀(𝑡)−𝜀(𝑡))⁵+−35(𝑧+𝜀(𝑡))⁴(𝛼𝜀(𝑡)−𝜀(𝑡))⁴+1. (14)

Let 𝐼_𝜀(𝑥,𝑡)={𝑖∣𝑔_𝑖(𝑥)>−𝛼𝜀(𝑡)}; by the definition of 𝜑(𝑧,𝑡), 𝑔(𝑥,𝑡) in (10) can be rewritten as 𝑔(𝑥,𝑡)=𝜃𝑡ln_{𝑖∈𝐼𝜀()𝑥,𝑡}𝜑𝑔_𝑖(𝑥)𝑔,𝑡exp^𝑖(𝑥)−𝜀(𝑡)𝜃𝑡+exp.𝜃𝑡 (15) The gradient of 𝑔(𝑥,𝑡) with respect to 𝑥 is ∇_𝑥𝑔(𝑥,𝑡)=_{𝑖∈𝐼𝜀()𝑥,𝑡}𝜆_𝑖(𝑥,𝑡)+𝜉_𝑖(𝑥,𝑡)∇𝑔_𝑖(𝑥), (16) where 𝜆_𝑖()=𝜑𝑔𝑥,𝑡_𝑖(𝑥),𝑡exp𝑔_𝑖(𝑥)/𝜃𝑡∑_{𝑗∈𝐼𝜀()𝑥,𝑡}𝜑𝑔_𝑗(𝑥),𝑡exp𝑔_𝑗(𝑥)((𝑡)),/𝜃𝑡+exp−𝜀/𝜃𝑡𝜉_𝑖()=𝑔𝑥,𝑡𝜃𝑡𝜕𝜑_𝑖(𝑥),𝑡/𝜕𝑧exp𝑔_𝑖(𝑥)/𝜃𝑡∑_{𝑗∈𝐼𝜀()𝑥,𝑡}𝜑𝑔_𝑗(𝑥),𝑡exp𝑔_𝑗(𝑥)((𝑡))./𝜃𝑡+exp−𝜀/𝜃𝑡 (17)(18) Then the gradient and Hessian of 𝑔(𝑥,𝑡) only relate to a part of constraint functions; that is, 𝑖∈𝐼_𝜀(𝑥,𝑡).

Proposition 8. For any 𝜃>0, 𝑡∈(0,1], 𝑥∈Ω with 𝜃𝑡↓0, 𝑡→𝑡∈[0,1], 𝑥→𝑥∈𝜕Ω,(1)𝜆_𝑖(𝑥,𝑡)→0  for  𝑖∉𝐼(𝑥);(2)∑_{𝑖∈𝐼(𝑥)}𝜆_𝑖(𝑥,𝑡)→1;(3)𝜉(𝑥,𝑡)→0.

Proof. (1) Because 𝑥∈𝜕Ω, there exists 𝑖^∗∈𝐼(𝑥) such that 𝑔_𝑖∗(𝑥)=0. By the continuity of 𝜑(𝑧,𝑡) and 𝑔(𝑥), we know that lim_{𝑥→𝑥,𝑡→𝑡}𝜑(𝑔_𝑖∗(𝑥),𝑡)=𝜑(𝑔_𝑖∗(𝑥),𝑡)=1 and lim_𝑥→𝑥(𝑔_𝑖(𝑥)−𝑔_𝑖∗(𝑥))=𝑔_𝑖(𝑥)<0 for 𝑖∉𝐼(𝑥); hence 0≤𝜆_𝑖()≤𝑔𝑥,𝑡exp_𝑖(𝑥)/𝜃𝑡𝜑𝑔_𝑖∗(𝑥),𝑡exp𝑔_𝑖∗(𝑥)=𝑔/𝜃𝑡exp_𝑖(𝑥)−𝑔_𝑖∗(𝑥)/𝜃𝑡𝜑𝑔_𝑖∗(𝑥),𝑡⟶0. (19) Hence, item 1 is satisfied for 𝜃𝑡↓0, 𝑡→𝑡, and 𝑥→𝑥.
(2) By lim_{𝑥→𝑥,𝑡→𝑡}(−𝜀(𝑡)−𝑔_𝑖∗(𝑥))=−𝜀(𝑡)<0, we have 1≥_{𝑖∈𝐼𝜀()𝑥,𝑡}𝜆_𝑖()=∑𝑥,𝑡_{𝑖∈𝐼𝜀()𝑥,𝑡}𝜑𝑔_𝑖(𝑥),𝑡exp𝑔_𝑖(𝑥)/𝜃𝑡∑_{𝑗∈𝐼𝜀()𝑥,𝑡}𝜑𝑔_𝑗(𝑥),𝑡exp𝑔_𝑗(𝑥)((𝑡))≥𝜑𝑔/𝜃𝑡+exp−𝜀/𝜃𝑡_𝑖∗(𝑥),𝑡exp𝑔_𝑖∗(𝑥)/𝜃𝑡𝜑𝑔_𝑖∗(𝑥),𝑡exp𝑔_𝑖∗(𝑥)((𝑡))≥𝜑𝑔/𝜃𝑡+exp−𝜀/𝜃𝑡_𝑖∗(𝑥),𝑡1+exp−𝜀(𝑡)−𝑔_𝑖∗(𝑥)/𝜃𝑡⟶1. (20) Hence, ∑_{𝑖∈𝐼𝜀(𝑥,𝑡)}𝜆_𝑖(𝑥,𝑡)→1. Together with item 1, we have ∑_{𝑖∈𝐼(𝑥)}𝜆_𝑖(𝑥,𝑡)→1 for 𝜃𝑡↓0, 𝑡→𝑡, and 𝑥→𝑥.
(3) If 𝑔_𝑖(𝑥)≥−𝜀(𝑡), then 𝜕𝜑(𝑔_𝑖(𝑥),𝑡)/𝜕𝑧=0 by its definition, and hence 𝜉_𝑖(𝑥,𝑡)=0; (21) else 𝜉_𝑖(≤𝑔𝑥,𝑡)𝜃𝑡𝜕𝜑_𝑖(𝑥),𝑡/𝜕𝑧exp𝑔_𝑖(𝑥)/𝜃𝑡exp(−𝜀(𝑡)/𝜃𝑡)≤𝜃𝑡𝑀, (22) where 𝑀=max_{(𝑧,𝑡)∈𝑅×[0,1]}|𝜕𝜑(𝑧,𝑡)/𝜕𝑧|. Then, we have 𝜉(𝑥,𝑡)→0 for 𝜃𝑡↓0, 𝑡→𝑡, and 𝑥→𝑥.

Proposition 9. For any given 0<𝜃≤1,(1) −𝜀(𝑡)≤𝑔(𝑥,𝑡)≤max{𝑔_max(𝑥),−𝜀(𝑡)}+𝜃𝑡ln(𝑚_𝐼𝜀+1),(2) 𝑔_max(𝑥)≤𝑔(𝑥,𝑡)≤𝑔_max(𝑥)+𝜃𝑡ln(𝑚_𝐼𝜀+1)  when  𝑔_max(𝑥)≥−𝜀(𝑡),(3) 𝑔(𝑥,𝑡)=−𝜀(𝑡)  when  𝑔_max(𝑥)≤−𝛼𝜀(𝑡)  with  𝛼>1,where 𝑚_𝐼𝜀 denotes the cardinality of 𝐼_𝜀(𝑥,𝑡).

Proof. (1) By 𝜑(𝑧,𝑡)≥0, we have 𝑔(𝑥,𝑡)=𝜃𝑡ln_{𝑖∈𝐼𝜀()𝑥,𝑡}𝜑𝑔_𝑖(𝑥)𝑔,𝑡exp^𝑖(𝑥)−𝜀(𝑡)𝜃𝑡+exp−𝜀(𝑡)𝜃𝑡≥𝜃𝑡lnexp(𝑡).𝜃𝑡=−𝜀 (23) Then the left inequality can be obtained. By 𝜑(𝑧,𝑡)≤1, we have 𝑔(𝑥,𝑡)=𝜃𝑡ln_{𝑖∈𝐼𝜀()𝑥,𝑡}𝜑𝑔_𝑖(𝑥)𝑔,𝑡exp^𝑖(𝑥)−𝜀(𝑡)𝜃𝑡+exp𝜃𝑡≤𝜃𝑡ln_{𝑖∈𝐼𝜀()𝑥,𝑡}𝑔exp^𝑖(𝑥)−𝜀(𝑡)𝜃𝑡+exp.𝜃𝑡 (24) Then the right inequality can be obtained by Proposition 2.3 in [19].
(2) For any 𝑥,𝑡 with 𝑔_max(𝑥)≥−𝜀(𝑡), let 𝑔_𝑖∗(𝑥)=𝑔_max(𝑥); then 𝜑(𝑔_𝑖∗(𝑥),𝑡)=1 by its definition, and hence 𝑔(𝑥,𝑡)=𝜃𝑡ln_{𝑖∈𝐼𝜀()𝑥,𝑡}𝜑𝑔_𝑖(𝑥)𝑔,𝑡exp^𝑖(𝑥)−𝜀(𝑡)𝜃𝑡+exp𝜑𝑔𝜃𝑡≥𝜃𝑡ln_𝑖∗(𝑥)𝑔,𝑡exp^𝑖∗(𝑥)𝑔𝜃𝑡=𝜃𝑡lnexp^𝑖∗(𝑥)𝜃𝑡=𝑔_𝑖∗(𝑥)=𝑔_max(𝑥); (25) then the left inequality is true. The right inequality can be obtained by item 1.
(3) It is trivial by the definition of 𝑔(𝑥,𝑡).

Proposition 10. Let Ω_𝜃(𝑡)={𝑥∣𝑔(𝑥,𝑡)≤0} and Ω_𝜃(𝑡)⁰={𝑥∣𝑔(𝑥,𝑡)<0}; one has(1)Ω_𝜃(𝑡)⊂Ω  for  𝑡∈(0,1];(2)for any bounded and closed set  Ω⊂Ω⁰, there exists a  𝜃∈(0,1], such that  Ω⊂Ω_𝜃(1)⁰;(3)Ω_𝜃(𝑡)→Ω  as  𝑡↓0.

Proof. (1) If 𝑥∉Ω, which is equivalent to 𝑔_max(𝑥)>0, by Proposition 9(2), 0<𝑔_max(𝑥)≤𝑔(𝑥,𝑡), (26) which means 𝑥∉Ω_𝜃(𝑡); then we have Ω_𝜃(𝑡)⊂Ω.
(2) By the continuity of 𝑔(𝑥) and the fact that Ω is a bounded closed set, there exists a point 𝑥^∗∈Ω such that 𝑔_max(𝑥) reaches its maximum in Ω, and 𝑔_max(𝑥^∗)<0. For any 𝜃∈(0,1] satisfying 𝜃≤−max{𝑔_max(𝑥^∗),−𝜀(1)}/(2ln(𝑚+1)), by Proposition 9(1), for any 𝑥∈Ω, 𝑔(𝑥,1)≤max𝑔_max(𝑥)(1()𝑔,−𝜀)+𝜃ln𝑚+1≤max_max𝑥^∗(1()=1,−𝜀)+𝜃ln𝑚+12max𝑔_max𝑥^∗,−𝜀(1)<0, (27) which means Ω⊂Ω_𝜃(1)⁰.
(3) For any 𝑥∈Ω, that is, 𝑔_max(𝑥)≤0, by the right inequality of Proposition 9(1), lim_𝑡↓0𝑔(𝑥,𝑡)≤lim_𝑡↓0max𝑔_max(𝑥)(𝑡𝑔,−𝜀)≤max_max(𝑥),−𝑐₂≤0, (28) where the second inequality comes from −𝜀(𝑡)=−𝑐₁𝑡−𝑐₂≤−𝑐₂<0 for 𝑡∈[0,1], which means Ω⊂lim_𝑡↓0Ω_𝜃(𝑡). Then, together with item 1, lim_𝑡↓0Ω_𝜃(𝑡)=Ω. (29)

Proposition 11. There exists a 𝜃∈(0,1] such that, for any 𝑡∈(0,1], 𝑥∈𝜕Ω_𝜃(𝑡); then ∇_𝑥𝑔(𝑥,𝑡)≠0.

Proof. If not, for any 𝜃∈(0,1], there exist corresponding 𝑡_𝜃∈(0,1] and 𝑥^𝜃∈𝜕Ω_𝜃(𝑡_𝜃) such that ∇_𝑥𝑔(𝑥^𝜃,𝑡_𝜃)=0, which means that there must exist three sequences {𝜃_𝑘}^∞_𝑘=1⊂(0,1], {𝑡_𝑘}^∞_𝑘=1⊂(0,1], and {𝑥^𝑘}^∞_𝑘=1 with 𝑥^𝑘∈𝜕Ω_𝜃𝑘(𝑡_𝑘), such that 𝜃_𝑘↓0, 𝑡_𝑘→𝑡, 𝑥^𝑘→𝑥∈𝜕Ω, and ∇_𝑥𝑔𝑥^𝑘,𝑡_𝑘=_{𝑖∈𝐼𝜀𝑥𝑘,𝑡𝑘}𝜆_𝑖𝑥^𝑘,𝑡_𝑘+𝜉_𝑖𝑥^𝑘,𝑡_𝑘∇𝑔_𝑖𝑥^𝑘⟶0 (30) as 𝑘→∞. Then, by Proposition 8, we know that 𝜆(𝑥^𝑘,𝑡_𝑘)→𝜆, 𝜉(𝑥^𝑘,𝑡_𝑘)→0 as 𝑘→∞, and _{𝑖∈lim𝑘→∞𝐼𝜀𝑥𝑘,𝑡𝑘}𝜆_𝑖=_{(𝑖∈𝐼𝑥)}𝜆_𝑖=1, (31) where the first equality comes from that 𝜆_𝑖=0 for 𝑖∉𝐼(𝑥). Therefore, we have ∑_{𝑖∈𝐼(𝑥)}𝜆_𝑖∇𝑔_𝑖(𝑥)=0 by taking limits on (30); this is a contradiction to Assumption 6.

Proposition 12. For any bounded closed set Ω⊂Ω⁰, there exists a 𝜃∈(0,1] such that 𝑥+𝜆∇_𝑥𝑔(𝑥,1)∣𝜆>0∩Ω=0/, (32) for any 𝑥∈𝜕Ω_𝜃(1).

Proof. If not, there exist a bounded closed set Ω and four sequences {𝜃_𝑘}^∞_𝑘=1⊂(0,1], {𝑥^𝑘}^∞_𝑘=1 with 𝑥^𝑘∈𝜕Ω_𝜃𝑘(1), {𝑥^𝑘}^∞_𝑘=1⊂Ω, and {𝜆^𝑘}^∞_𝑘=1>0 such that 𝜃_𝑘↓0, 𝑥^𝑘→𝑥∈𝜕Ω, 𝑥^𝑘→𝑥∈Ω, 𝜆_𝑖(𝑥^𝑘,1)→𝜆^∗_𝑖≥0, 𝜉_𝑖(𝑥^𝑘,1)→0 as 𝑘→∞, and 𝑥^𝑘=𝑥^𝑘+𝜆^𝑘_{𝑖∈𝐼𝜀𝑥𝑘,1}𝜆_𝑖𝑥^𝑘,1+𝜉_𝑖𝑥^𝑘,1∇𝑔_𝑖𝑥^𝑘. (33) By Proposition 8, we have _{𝑖∈lim𝑘→∞𝐼𝜀𝑥𝑘,1}𝜆^∗_𝑖=_{(𝑖∈𝐼𝑥)}𝜆^∗_𝑖=1, (34) where the first equality comes from that 𝜆^∗_𝑖=0 for 𝑖∉𝐼(𝑥). Then, we have _{𝑖∈𝐼𝜀𝑥𝑘,1}𝜆_𝑖𝑥^𝑘,1+𝜉_𝑖𝑥^𝑘,1∇𝑔_𝑖𝑥^𝑘⟶_{(𝑖∈𝐼𝑥)}𝜆^∗_𝑖∇𝑔_𝑖(𝑥). (35) By Assumption 6 and using (34), _{(𝑖∈𝐼𝑥)}𝜆^∗_𝑖∇𝑔_𝑖(𝑥)≠0. (36) Hence {𝜆^𝑘}^∞_𝑘=1 is a bounded sequence, and there must exist a subsequence converging to 𝜆>0; then we have by (33) 𝑥=𝑥+𝜆_{(𝑖∈𝐼𝑥)}𝜆^∗_𝑖∇𝑔_𝑖(𝑥), (37) which contradicts Assumption 7.

3. The Flattened Aggregate Constraint Homotopy Method

In [19], the following aggregate constraint homotopy was introduced: 𝐻()=()(𝑥)𝑥,𝜆,𝑡1−𝑡∇𝑓+𝜆∇_𝑥𝑔_𝜃(𝑥,𝑡)+𝑡𝑥−𝑥⁰𝜆𝑔_𝜃(𝑥,𝑡)−𝑡𝜆⁰𝑔_𝜃𝑥⁰,1=0, (38) where 𝑔_𝜃(𝑥,𝑡)=𝜃𝑡ln∑^𝑚_𝑖=1exp(𝑔_𝑖(𝑥)/𝜃𝑡) with 𝜃∈(0,1], 𝜆 is the Lagrangian multiplier of the aggregate constraint function 𝑔_𝜃(𝑥,𝑡), and (𝑥⁰,𝜆⁰) is the starting point in Ω×𝑅₊₊. 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 (𝑥⁰,𝜆⁰,1) to (𝑥^∗,𝜆^∗,0), in which (𝑥^∗,𝜆^∗) is a KKT point of (1).

3.1. The Flattened Aggregate Constraint Homotopy

Using the flattened aggregate constraint function 𝑔(𝑥,𝑡) in (10), we construct the following flattened aggregate constraint homotopy: 𝐻()=()(𝑥)𝑥,𝜆,𝑡1−𝑡∇𝑓+𝜆∇_𝑥𝑔(𝑥,𝑡)+𝑡𝑥−𝑥⁰𝜆𝑔(𝑥,𝑡)−𝑡𝜆⁰𝑔𝑥⁰,1=0, (39) where 𝜆 is the Lagrangian multiplier of the flattened aggregate constraint and (𝑥⁰,𝜆⁰) is the starting point and can be randomly chosen from Ω×𝑅₊₊.

We give the main theorem on the existence and convergence of a smooth path from (𝑥⁰,𝜆⁰,1) to (𝑥^∗,𝜆^∗,0), in which (𝑥^∗,𝜆^∗) 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 Assumptions 5–7 hold and 𝜃∈(0,1] satisfies Propositions 10–12 for the bounded closed set Ω⊂Ω_𝜃(1)⁰⊂Ω⁰; then for almost all starting points 𝑤⁰=(𝑥⁰,𝜆⁰)∈Ω×𝑅₊₊, the zero-points set 𝐻⁻¹(0)={(𝑥,𝜆,𝑡)∣𝐻(𝑥,𝜆,𝑡)=0} defines a smooth path Γ_𝑤0, which starts at (𝑥⁰,𝜆⁰,1) and approaches the hyperplane 𝑡=0. Furthermore, let (𝑥^∗,𝜆^∗,0) be any limit point of Γ_𝑤0, and 𝑦^∗=𝜆^∗𝜆(𝑥^∗,0)  (𝜆(𝑥^∗,0) is a limit of 𝜆(𝑥,𝑡) as 𝑥→𝑥^∗,𝑡↓0), then 𝑥^∗ is a KKT point of (1) and 𝑦^∗ is the corresponding Lagrangian multiplier.

Proof. Consider 𝐻 as a map of the variable (𝑤⁰,𝑥,𝜆,𝑡), for any (𝑤⁰,𝑡)∈Ω×𝑅₊₊×(0,1], which means 𝑥⁰∈Ω_𝜃(1)⁰; hence 𝑔(𝑥⁰,1)<0, because 𝜕𝐻𝑤⁰,𝑥,𝜆,𝑡𝜕𝑤⁰=−𝑡𝐼0−𝑡𝜆⁰∇_𝑥𝑔𝑥⁰,1^𝑇−𝑡𝑔𝑥⁰,1 (40) is nonsingular and the Jacobi matrix of 𝐻(𝑤⁰,𝑥,𝜆,𝑡) is of full row rank. Using the parameterized Sard theorem, Lemma 4, for almost all 𝑤⁰∈Ω×𝑅₊₊, 0 is a regular value of 𝐻(𝑥,𝜆,𝑡). From 𝐻(𝑥,𝜆,1)=𝑥−𝑥⁰𝜆𝑔(𝑥,1)−𝜆⁰𝑔𝑥⁰,,1 (41) we know that (𝑥⁰,𝜆⁰) is the unique and simple solution of 𝐻(𝑥,𝜆,1)=0. Then, by the implicit function theorem, there exists a smooth curve Γ_𝑤0⊂Ω_𝜃(𝑡)⁰×𝑅₊₊×(0,1] starting from (𝑥⁰,𝜆⁰,1) and being transversal to the hyperplane 𝑡=1.
Since Γ_𝑤0 can be extended in Ω_𝜃(𝑡)⁰×𝑅₊₊×(0,1) until it converges to the boundary of Ω_𝜃(𝑡)⁰×𝑅₊₊×(0,1), there must exist an extreme point (𝑥,𝜆,𝑡)∈𝜕(Ω_𝜃(𝑡)×𝑅₊×[0,1]). Let (𝑥,𝜆,𝑡) be any extreme point of Γ_𝑤0 other than (𝑥⁰,𝜆⁰,1); then only the following five cases are possible:(1)(𝑥,𝜆,𝑡)∈Ω×𝑅₊×{0}, 𝜆<+∞;(2)(𝑥,𝜆,𝑡)∈Ω_𝜃(1)×𝑅₊×{1}, 𝜆<+∞;(3)(𝑥,𝜆,𝑡)∈𝜕Ω_𝜃(𝑡)×𝑅₊₊×(0,1), 𝜆<+∞;(4)(𝑥,𝜆,𝑡)∈Ω_𝜃(𝑡)×{0}×(0,1);(5)(𝑥,𝜆,𝑡)∈Ω_𝜃(𝑡)×{+∞}×[0,1].
Case (2) means that 𝐻(𝑥,𝜆,1)=0 has another solution (𝑥,𝜆) except (𝑥⁰,𝜆⁰), or (𝑥⁰,𝜆⁰) is a double solution of 𝐻(𝑥,𝜆,1)=0, which contradicts the fact that (𝑥⁰,𝜆⁰) is the unique and simple solution of 𝐻(𝑥,𝜆,1)=0. Case (3) means 𝑔(𝑥,𝑡)=0, 0<𝜆<∞, and 0<𝑡<1; hence 𝜆𝑔𝑥,𝑡−𝑡𝜆⁰𝑔𝑥⁰,1=−𝑡𝜆⁰𝑔𝑥⁰,1<0, (42) which contradicts the last equation of (39). Case (4) means 𝑔(𝑥,𝑡)≤0, 𝜆=0, and 0<𝑡<1; hence 𝜆𝑔𝑥,𝑡−𝑡𝜆⁰𝑔𝑥⁰,1=−𝑡𝜆⁰𝑔𝑥⁰,1<0, (43) which contradicts the last equation of (39). Then, Cases (2), (3), and (4) are impossible.
Because 𝐻(𝑥,𝜆,1)=0 has only one and also unique solution (𝑥⁰,𝜆⁰), 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 {(𝑥^𝑘,𝜆^𝑘,𝑡_𝑘)}^∞_𝑘=1 on Γ_𝑤0 such that 𝑥^𝑘→𝑥, 𝜆^𝑘→∞, and 𝑡_𝑘→𝑡. By the last equation of (39), we have 𝑔𝑥^𝑘,𝑡_𝑘=𝑡_𝑘𝜆⁰𝑔𝑥⁰,1𝜆^𝑘⟶0 (44) as 𝑘→∞, which means 𝑥∈𝜕Ω_𝜃(𝑡). The following two cases may be possible.(1)If 𝑡=1, let 𝛼 be any accumulation point of {(1−𝑡_𝑘)𝜆^𝑘}, which is known to exist from (39) and lim_𝑘→∞∇_𝑥𝑔(𝑥^𝑘,𝑡_𝑘)≠0; then we have 𝑥⁰=𝑥+𝛼∇_𝑥𝑔(𝑥,1). (45) If 𝛼=0, then 𝑥=𝑥⁰ contradicts 𝑥⁰∈Ω_𝜃(1)⁰; else contradicts Proposition 12.(2)If 𝑡<1, by (39), we have 1−𝑡lim_𝑘→∞𝜆^𝑘∇_𝑥𝑔𝑥^𝑘,𝑡_𝑘=−1−𝑡∇𝑓(𝑥)−𝑡𝑥−𝑥⁰, (46) the right-hand is finite; however, by lim_𝑘→∞∇_𝑥𝑔(𝑥^𝑘,𝑡_𝑘)≠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 Γ_𝑤0 must approach the hyperplane 𝑡=0. Because {𝜆(𝑥^𝑘,𝑡_𝑘)}^∞_𝑘=1(defined in (17)), {𝑥^𝑘}^∞_𝑘=1, and {𝜆^𝑘}^∞_𝑘=1 are bounded sequences, we know that {(𝑥^𝑘,𝜆^𝑘,𝜆(𝑥^𝑘,𝑡_𝑘))}^∞_𝑘=1 has at least one accumulation point as 𝑘→∞. Let (𝑥^∗,𝜆^∗,𝜆(𝑥^∗,0)) be any accumulation point; 𝑦^∗_𝑖=𝜆^∗𝜆_𝑖(𝑥^∗,0). By (39) and the fact that 𝜉(𝑥,𝑡)→0 as 𝑡↓0, we have ∇𝑓𝑥^∗+^𝑚_𝑖=1𝑦^∗_𝑖∇𝑔_𝑖𝑥^∗=0. (47)
By the fact that lim_𝑡↓0Ω_𝜃(𝑡)=Ω, we have 𝑥^∗∈Ω. If 𝑥^∗∈Ω⁰, we know that lim_{𝑥→𝑥∗,𝑡↓0}𝑔(𝑥,𝑡)≤lim_{𝑥→𝑥∗,𝑡↓0}max𝑔_max(𝑥),−𝜀(𝑡)+𝜃𝑡ln𝑚_𝐼𝜀𝑔+1=max_max𝑥^∗(0𝑔,𝜀)=max_max𝑥^∗,−𝛼𝑐₂<0, (48) 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 𝜀(𝑡) in (11); hence 𝜆^∗=0 by the last equation of (39); else 𝑥^∗∈𝜕Ω, lim_{𝑥→𝑥∗,𝑡↓0}𝑔(𝑥,𝑡)=𝑔_max(𝑥^∗)=0 by Proposition 9(2) and 𝜆_𝑖(𝑥^∗,0)=0 for 𝑖∉𝐼(𝑥^∗) by Proposition 8. Thus we have 𝑦^∗_𝑖𝑔_𝑖𝑥^∗=0,𝑖=1,…,𝑚. (49)
Summing up, (𝑥^∗,𝑦^∗) is a solution of the KKT system of (1), which means that 𝑥^∗ is a KKT point of (1) and 𝑦^∗ is the corresponding Lagrangian multiplier.

3.2. The Modified Flattened Aggregate Constraint Homotopy

From Proposition 8(3), we know that 𝜉(𝑥,𝑡)→0 as 𝑡↓0 for any 𝑥∈Ω and, hence, we can use the following modified flattened aggregate constraint homotopy (MFACH) instead of the FACH: 𝐻()=()(𝑥)()𝑥,𝜆,𝑡1−𝑡∇𝑓+𝜆ℎ𝑥,𝑡+𝑡𝑥−𝑥⁰𝜆𝑔(𝑥,𝑡)−𝑡𝜆⁰𝑔𝑥⁰,1=0, (50) where ℎ(𝑥,𝑡)=_{𝑖∈𝐼𝜀()𝑥,𝑡}𝜆_𝑖(𝑥,𝑡)∇𝑔_𝑖(𝑥). (51)

Remarks. (i)  Since 𝜉(𝑥,𝑡) 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 that 𝐻(𝑥,𝜆,𝑡) in (39) be a 𝐶² map, 𝜑(𝑧,𝑡) must be 𝐶³; hence, if  𝜑(𝑧,𝑡) is chosen as a polynomial, the total degree should be seven; in contrast, for the homotopy map in (50), 𝜑(𝑧,𝑡) is only needed to be 𝐶². Then 𝜑(𝑧,𝑡) should satisfy 𝜑((𝑡))𝜕−𝛼𝜀,𝑡=0,𝜑(−𝛼𝜀(𝑡),𝑡)𝜕𝜕𝑧=0,₂𝜑(−𝛼𝜀(𝑡),𝑡)𝜕𝑧²=0,𝜑((𝑡))𝜕−𝜀,𝑡=1,𝜑(−𝜀(𝑡),𝑡)𝜕𝜕𝑧=0,₂𝜑(−𝜀(𝑡),𝑡)𝜕𝑧²=0 (52) and 𝜑(𝑧,𝑡) can be chosen as 𝜑()=6((𝑡𝑧,𝑡𝑧+𝜀))⁵(𝛼𝜀(𝑡)−𝜀(𝑡))⁵+15(𝑧+𝜀(𝑡))⁴(𝛼𝜀(𝑡)−𝜀(𝑡))⁴+10(𝑧+𝜀(𝑡))³(𝛼𝜀(𝑡)−𝜀(𝑡))³+1. (53)

(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 direction 𝑑^𝑘 and a steplength ℎ_𝑘 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 𝛽^𝑘 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 (𝑤^(𝑘,0),𝑡_(𝑘,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 (𝑤^(𝑘,𝑖),𝑡_(𝑘,𝑖)), 𝑖=1,…,5, are calculated by (𝑤^{(𝑘,𝑖+1)},𝑡_{(𝑘,𝑖+1)})=(𝑤^(𝑘,𝑖),𝑡_(𝑘,𝑖))+𝑑^{(𝑘,𝑖+1)}, 𝑖=0,…,4, where the step 𝑑^{(𝑘,𝑖+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 𝑑^(𝑘,𝑖) and 𝐻(𝑤^(𝑘,𝑖),𝑡_(𝑘,𝑖)) satisfy the tolerance criterions.

At each predictor step and corrector step, the feasibility of the predictor point (𝑤^𝑘,𝑡_𝑘) and the corrector point (𝑤^(𝑘,𝑖),𝑡_(𝑘,𝑖)) needs to be checked. If 𝑡_𝑘<𝑡_𝑐 in the predictor step or 𝑡_𝑘<0 in the corrector step, a damping step is used to get a new point (𝑤^(𝑘,0),0). Then, if 𝑤^(𝑘,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 𝑤^(𝑘,0), Newton's method is used to solve 𝐹_𝑡𝑐()=(𝑥)𝑥,𝜆∇𝑓+𝜆ℎ𝑥,𝑡_𝑐𝜆𝑔𝑥,𝑡_𝑐=0, (54) where 𝑡_𝑐 is a small positive constant. For other situations, the steplength will be decreased to make new predictor-corrector steps (see Algorithm 1).

Input  Give 𝜃∈(0,1], 𝑐1>0, 𝑐2>0, 𝛼>1, 𝑡end and 𝑡𝑐; starting point 𝑤0∈Ω×𝑅++; initial steplength ℎ1, steplength    contraction factor 𝐵min, steplength expansion factor 𝐵max; tracking tolerances 𝐻tol and 𝐻final for correction. Set    𝑡0=1, 𝑘=1, IT=0, 𝑁good=2, 𝛽1=0. Step  1. The predictor step.   Step  1.1. Compute the predictor direction.      If 𝑘=1, compute 𝑑 by solving                    𝐻′𝑥0,𝜆0𝑑,10𝑇𝑑=−𝑑0,      where 𝑑0=(0,…,0,−1)∈𝑅𝑛+2, set the predictor direction 𝑑1=𝑑/‖𝑑‖;      else, compute the predictor direction 𝑑𝑘, the angle 𝛽𝑘 between 𝑑𝑘 and 𝑑𝑘−1:               𝑑𝑘=(𝑤𝑘−1,𝑡𝑘−1)−(𝑤𝑘−2,𝑡𝑘−2)‖(𝑤𝑘−1,𝑡𝑘−1)−(𝑤𝑘−2,𝑡𝑘−2)‖,𝛽𝑘=arccos𝑑𝑘𝑇𝑑𝑘−1.   Step  1.2. Adjust the steplength.        Adjust the steplength ℎ𝑘 as follows:        if the predictor point is infeasible, set ℎ𝑘=𝐵min(1)ℎ𝑘, 𝑁good=0;        else if the corrector step fails or 𝛽𝑘>𝜋/4, set ℎ𝑘=𝐵min(1)ℎ𝑘−1, 𝑁good=0;        else if 𝑖=5, set ℎ𝑘=𝐵min(2)ℎ𝑘−1, 𝑁good=0;        else if 𝑖=4, set ℎ𝑘=ℎ𝑘−1, 𝑁good=𝑁good+1;        else if 𝑖=3, set 𝑁good=𝑁good+1,          if 𝑁good>2, set ℎ𝑘=min{1,𝐵max(2)ℎ𝑘−1}, else, set ℎ𝑘=ℎ𝑘−1;        else, set 𝑁good=𝑁good+1,          if 𝑁good>2, set ℎ𝑘=min{1,𝐵max(1)ℎ𝑘−1}, else, set ℎ𝑘=ℎ𝑘−1.       If ℎ𝑘<10−10, stop the algorithm with an error flag.   Step  1.3. Compute the predictor point and check its feasibility.         Compute the predictor point                   𝑤(𝑘,0),𝑡(𝑘,0)=𝑤𝑘−1,𝑡𝑘−1+ℎ𝑘𝑑𝑘.      If 𝑤(𝑘,0)∉Ω𝜃(𝑡(𝑘,0))0×𝑅++ or 𝑡(𝑘,0)≥1,        the predictor point is infeasible, goto Step  1.2;      else if 𝑡(𝑘,0)≤𝑡end,        adjust the steplength ℎ𝑘, compute the point (𝑤(𝑘,0),0), if 𝑤(𝑘,0)∈Ω×𝑅+, goto Step  3, else, the predictor        point is infeasible, goto Step  1.2;      else, goto Step  2. Step  2. The corrector step.    Set 𝑖=0.    Repeat      If 𝑖=5, set IT=IT+5, 𝑤𝑘=𝑤𝑘−1, 𝑡𝑘=𝑡𝑘−1, 𝑑𝑘+1=𝑑𝑘, 𝑘=𝑘+1, the corrector step fails, goto Step  1.2.      Compute the Newton step 𝑑(𝑘,𝑖+1) by solving                𝐻′𝑤(𝑘,𝑖),𝑡(𝑘,𝑖)𝑑𝑘𝑇𝑑(𝑘,𝑖+1)=−𝐻𝑤(𝑘,𝑖),𝑡(𝑘,𝑖)0,      the corrector point                𝑤(𝑘,𝑖+1),𝑡(𝑘,𝑖+1)=𝑤(𝑘,𝑖),𝑡(𝑘,𝑖)+𝑑(𝑘,𝑖+1);      set 𝑖=𝑖+1.      If 𝑤(𝑘,𝑖)∉Ω𝜃(𝑡(𝑘,𝑖))0×𝑅++ or 𝑡(𝑘,𝑖)≥1,        set IT=IT+𝑖, 𝑤𝑘=𝑤𝑘−1, 𝑡𝑘=𝑡𝑘−1, 𝑑𝑘+1=𝑑𝑘, 𝑘=𝑘+1, the corrector step fails, goto Step  1.2;      else if 𝑡(𝑘,𝑖)<0,        compute the point (𝑤(𝑘,0),0) by a damping Newton step, set IT=IT+𝑖, if 𝑤(𝑘,0)∈Ω×𝑅+, goto Step  3,        else, set 𝑤𝑘=𝑤𝑘−1, 𝑡𝑘=𝑡𝑘−1, 𝑑𝑘+1=𝑑𝑘, 𝑘=𝑘+1, the corrector step fails, goto Step  1.2.   Until ‖𝐻(𝑤(𝑘,𝑖),𝑡(𝑘,𝑖))‖inf≤𝐻tol and ‖𝑑(𝑘,𝑖)‖≤𝐻tol.   If 𝑡𝑘<𝑡𝑐, return with 𝑥∗=𝑥𝑘, stop the algorithm; else, set 𝐻tol=min{𝐻tol,𝑡𝑘}, 𝑘=𝑘+1, goto Step  1. Step  3. The end game.    Set 𝑙=0.    Repeat      If 𝑙=5 or ‖𝑑(𝑘,𝑙)‖>‖𝑑(𝑘,𝑙−1)‖ for 𝑙>1, or 𝑤(𝑘,𝑙)∉Ω×𝑅+, set 𝑡end=0.3∗𝑡end, 𝑤𝑘=𝑤𝑘−1, 𝑡𝑘=𝑡𝑘−1,      𝑑𝑘+1=𝑑𝑘, IT=IT+𝑙, 𝑘=𝑘+1, the corrector step fails, goto Step  1.2.      Compute the Newton step 𝑑(𝑘,𝑙+1) by solving 𝐹′𝑡𝑐(𝑤(𝑘,𝑙))𝑑(𝑘,𝑙+1)=−𝐹𝑡𝑐(𝑤(𝑘,𝑙)), the corrector point                    𝑤(𝑘,𝑙+1)=𝑤(𝑘,𝑙)+𝑑(𝑘,𝑙+1);      set 𝑙=𝑙+1.                Until ‖𝐹𝑡𝑐(𝑤(𝑘,𝑙))‖inf≤𝐻final and ‖𝑑(𝑘,𝑙)‖≤𝐻final.                Set IT=IT+𝑙, return with 𝑥∗=𝑥(𝑘,𝑙), stop the algorithm.

Algorithm 1 

The FACH-S-N procedure.

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 𝑥^𝑘, if 𝑔_max(𝑥^𝑘)<10⁻⁶, 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, 𝑐₁=0.05, 𝑐₂=0.5∗10⁻⁵, and 𝛼=2.(ii)Parameters in the end game section: 𝑡_𝑐=10⁻⁶ and 𝑡_end=0.1.(iii)Step size parameters: ℎ₀=0.1, 𝐵_min=[0.5,0.75], and 𝐵_max=[3,1.5].(iv)Tracking tolerances: 𝐻_tol=10⁻⁵ and 𝐻_final=10⁻¹².(v)Initial Lagrangian multipliers: 1 for ACH, FACH, and MFACH methods.

For each problem with different parameters, we list the value of objective function 𝑓(𝑥) and the maximal function of constraints 𝑔_max(𝑥) at 𝑥^∗, the number of iterations 𝐼𝑇, the number of evaluations of gradients of individual constraint functions in the FACH and MFACH methods 𝑁_𝑠 (in contrast, 𝑁_𝑠 is 𝑚∗𝐼𝑇 in the ACH method), and the running time in seconds 𝑇𝑖𝑚𝑒. For problems that were not solved by the conservative setting, we also give the reason for failure. The notation “fail¹” indicates that the steplength in predictor step is smaller than 10⁻¹⁰; it is generally due to poor conditioned Jacobian matrix. The notation “fail²” means out of memory. The notation “fail³” means no result in 5000 iterations or 3600 seconds.

Example 14 (see [21]). Consider 𝑓(𝑥)=sin𝑥₁+−1+1.5𝜋_{2≤𝑖≤1000}100sin−𝑥_𝑖+1.5𝜋+𝑥²_𝑖−1,𝑔_𝑖(𝑥)=𝑥₁𝑥−𝜋𝑖=1;_2𝑖−1−𝑥_𝑖−𝜋𝑖=2,…,1000;−𝑥₁−𝜋𝑖=1001;−𝑥²_𝑖−1+𝑥_𝑖𝑥−𝜋𝑖=1002,…,2000,₀=(1,…,1)∈𝑅¹⁰⁰⁰. (55)

Remark 15. In this example, the starting point happens to be an unconstrained minimizer of the objective function and we found that the 𝑥-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 (𝑥⁰,0)∈Ω⁰×𝑅^𝑚₊ is a solution of the KKT system, that is, 𝑥⁰ is a stationary point of the objective function 𝑓(𝑥), when the parameters 𝛼, 𝑐₁, and 𝑐₂ in FACH and MFACH methods satisfy 𝑔_max(𝑥⁰)≤−𝛼(𝑐₁+𝑐₂), the 𝑥-components of points on the homotopy path Γ_𝑤0 will remain invariant, which can be derived from the fact that (𝑥⁰,𝑡𝜆⁰𝜀(1)/𝜀(𝑡),𝑡) is the solution of (39) and (50) for 𝑡∈(0,1]. Moreover, since 𝐻^′𝑥⁰,𝜆⁰=0𝑔𝑥,1𝐼00₀∗,,1 (56) the 𝑥-component of 𝑑¹ is 0. By a similar discussion, the 𝑥-component of 𝑑^(1,𝑖) is 0, and hence 𝑥¹=𝑥^(1,𝑖)=𝑥⁰. In this analogy, we know that 𝑥^𝑘=𝑥^(𝑘,𝑖)=𝑥⁰ for any 𝑘 and 𝑖 in algorithm FACH-S-N.

Example 16. Consider 𝑓(𝑥)=𝑥²₃+𝑥²₄,𝑔_𝑖,𝑗(𝑥)=𝑡_𝑖−𝑥₁²𝑥²₃+𝑡^′_𝑗−𝑥₂²𝑥²₄𝑡−1,_𝑖=𝑖√√𝑚−1,𝑖=0,…,𝑡𝑚−1,_′𝑗=𝑗√√𝑚−1,𝑗=0,…,𝑥𝑚−1,₀=(0,0,100,100)∈𝑅⁴. (57)

Example 17. Consider 𝑓(𝑥)=1𝑛_{1≤𝑘≤𝑛}𝑥_𝑘−1²,𝑔_𝑖(𝑥)=_{1≤𝑘≤𝑛}cos𝑡_𝑖𝑥_𝑘+𝑡_{𝑖1≤𝑘≤𝑛}𝑥³_𝑘,𝑡_𝑖=0.5+𝜋𝑖()𝑥𝑚−1,𝑖=0,…,𝑚−1,₀=(−2,…,−2)∈𝑅^𝑛. (58)

Example 18 (see [22]). Consider 𝑓(𝑥)=𝑥²₁+𝑥²₂+𝑥²₃,𝑔_𝑖(𝑥)=𝑥₁+𝑥₂exp𝑥₃𝑡_𝑖+exp2𝑡_𝑖−2sin4𝑡_𝑖,𝑡_𝑖=𝑖()𝑥𝑚−1,𝑖=0,…,𝑚−1,₀=(−200,−200,200)∈𝑅³. (59)

Example 19 (see [22]). Consider 𝑓(𝑥)=𝑥₁−2𝑥₂+5𝑥²₂−𝑥³₂−13²+𝑥₁−14𝑥₂+𝑥²₂+𝑥³₂−29²,𝑔_𝑖(𝑥)=𝑥²₁+2𝑥₂𝑡²_𝑖+exp𝑥₁+𝑥₂−exp𝑡_𝑖,𝑡_𝑖=𝑖()𝑥𝑚−1,𝑖=0,…,𝑚−1,₀=(0,−45)∈𝑅². (60)

Example 20 (see [22]). Consider 𝑓(𝑥)=𝑥²₁3+𝑥₁2+𝑥²₂,𝑔_𝑖(𝑥)=1−𝑥²₁𝑡²_𝑖²−𝑥₁𝑡²_𝑖−𝑥²₂+𝑥₂,𝑡_𝑖=𝑖()𝑥𝑚−1,𝑖=0,…,𝑚−1,₀=(−1,100)∈𝑅². (61) To explain the numerical efficiency, we make the following remarks by numerical results in Tables 1, 2, 3, 4, 5 and 6.(i)If 𝑔_max(𝑥)≤−𝛼𝜀(𝑡) for any (𝑥,𝜆,𝑡)∈Γ_𝑤0, 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 when 𝑚 is small but much better when 𝑚 is large for most problems. In addition, we can see that their time cost increases more slowly than KNITRO solver as 𝑚 increases.(iv)The function 𝜑(𝑧,𝑡) is important for the flattened aggregate constraint function. Theoretically, the parameters 𝑐₁, 𝑐₂, 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 𝑐₁=0.05, 𝑐₂=0.5∗10⁻⁵, and 𝛼=2. In addition, the function 𝜑(𝑧,𝑡) can be defined in many ways, and preliminary numerical experiments show that the algorithms with different functions 𝜑(𝑧,𝑡) 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 Example 14.

Table 2 

Test results for Example 16.

Table 3 

(a) Test results for Example 17 with 𝑛=100. (b) Test results for Example 17 with 𝑛=500. (c) Test results for Example 17 with 𝑛=1000. (d) Test results for Example 17 with 𝑛=2000.

Table 4 

Test results for Example 18.

Table 5 

Test results for Example 19.

Table 6 

Test results for Example 20.

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.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgment

The research was supported by the National Natural Science Foundation of China (11171051, 91230103).