2.9 Numerical methods for nonlinear programming problems

Nowadays, strong and efficient mathematical programming techniques are available for solving a great variety of nonlinear problems, which are based on solid theoretical results and extensive numerical studies. Approximated functions, derivatives and optimal solutions can be employed together with optimization algorithms to reduce the computational time. The aim of this section is not to describe state-of-the-art algorithms in nonlinear programming, but to explain, in a simple way, a number of modern algorithms for solving nonlinear problems. These techniques are typically iterative in the sense that, given an initial point $x{}^{}0$ , a sequence of points, $\left\{x{}^{}k\right\}$ , is obtained by repeated application of an algorithmic rule. The objective is to make this sequence converge to a point $\overline{x}$ in a pre-specified solution set. Typically, the solution set is specified in terms of optimality conditions, such as those presented in 2.5 through 2.8.

We start by recalling a number of concepts in 2.9.1. Then, we discuss the principles of Newton-like algorithms for nonlinear systems in 2.9.2, and use these concepts for the solution of unconstrained optimization problems in 2.9.3. Finally, algorithms for solving general nonlinear problems are presented in 2.9.7, with emphasizes on sequential unconstrained minimization (SUM) and sequential quadratic programming (SQP) techniques.

2.9.1 Preliminaries

Two essential questions must be addressed concerning iterative algorithms. The first question, which is qualitative in nature, is whether a given algorithm in some sense yields, at least in the limit, a solution to the original problem; the second question, the more quantitative one, is concerned with how fast the algorithm converges to a solution. We elaborate on these concepts in this subsection.

The convergence of an algorithm is said to be asymptotic when the solution is not achieved after a finite number of iterations; except for particular problems such as linear and quadratic programming, this is generally the case in nonlinear programming. That is, a very desirable property of an optimization algorithm is global convergence:

Most modern mathematical programming algorithms are globally convergent. Locally convergent algorithms are not useful in practice because the neighborhood of convergence is not known in advance and can be arbitrarily small.

Next, what distinguishes optimization algorithms with the global convergence property is the order of convergence:

Since they involve the limit when $k\to \infty$ , $p$ and $\beta$ are a measure of the asymptotic rate of convergence, i.e., as the iterates gets close to the solution; yet, a sequence with a good order of convergence may be very ”slow” far from the solution. Clearly, the convergence is faster when $p$ is larger and $\beta$ is smaller. Near the solution, if the convergence rate is linear, then the error is multiplied by $\beta$ at each iteration. The error reduction is squared for quadratic convergence, i.e., each iteration roughly doubles the number of significant digits. The methods that will be studied hereafter have convergence rates varying between linear and quadratic.

It should also be noted that for most algorithms, the user must set initial values for certain parameters, such as the starting point and the initial step size, as well as parameters for terminating the algorithm. Optimization procedures are often quite sensitive to these parameters, and may produce different results, or even stop prematurely, depending on their values. Therefore, it is crucial for the user to understand the principles of the algorithms used, so that he or she can select adequate values for the parameters and diagnose the reasons of a premature termination (failure).

2.9.2 Newton-like Algorithms for nonlinear Systems

The fundamental approach to most iterative schemes was suggested over 300 years ago by Newton. In fact, Newton’s method is the basis for nearly all the algorithms that are described herein.

Suppose one wants to find the value of the variable $x\in {ℝ}^{n_x}$ such that

where $\varphi :{ℝ}^{n_x}\to {ℝ}^{n_x}$ is continuously differentiable. Let us assume that one such solution exists, and denote it by $x{}^{}\star$ . Let also $x$ be a guess for the solution. The basic idea of Newton’s method is to approximate the nonlinear function $\varphi$ by the first two terms in its Taylor series expansion about the current point $x$ . This yields a linear approximation for the vector function $\varphi$ at the new point $\overline{x}$ ,

Using this linear approximation, and provided that the Jacobian matrix $\frac{\partial \varphi }{\partial x}$ is nonsingular, a new estimate for the solution $x{}^{}\star$ can be computed by solving the previous equation for $\varphi \left(x\right)=0$

Letting $d=-(\frac{\partial \varphi }{\partial x}|x{)}^{-1}\varphi \left(x\right)$ , we get the update $\overline{x}=x+d$ .

Of course, $\varphi$ being a nonlinear function, one cannot expect that $\varphi \left(\overline{x}\right)=0$ , but there is much hope that $\overline{x}$ be a better estimate for the root $x{}^{}\star$ than the original guess $x$ . In other words, we might expect that

If the new point is an improvement, then it makes sense to repeat the process, thereby defining a sequence of points $\{x{}^{}0,x{}^{}1,\dots \}$ . An algorithm implementing Newton’s method is as follows:

It can be shown that the rate of convergence for Newton’s method is quadratic (see Definition 2.24). Loosely speaking, it implies that each successive estimate of the solution doubles the number significant digits, which is a very desirable property for an algorithm to possess.

But can anything go wrong with Newton’s method?

Lack of Global Convergence First and foremost, if the initial guess is not sufficiently close to the solution, i.e., within the region of convergence, Newton’s method may diverge. Said differently, Newton’s method as presented above does not have the global convergence property (see Definition 2.23 and Example 2.20 hereafter). This is because $d{}^{}k:=\nabla \varphi {(x{}^{}}^k\left)-\top \varphi \right(x{}^{}k)$ may not be a valid descent direction far from the solution, and even if , a unit step size might not give a descent in $\varphi$ . Globalization strategies, which aim at correcting this latter deficiency, will be presented in 2.9.4 in the context of unconstrained optimization.

Singular Jacobian Matrix A second difficulty occurs when the Jacobian matrix $\nabla \varphi \left(x{}^{}k\right)$ becomes singular during the iteration process, since the correction defined by the Newton’s algorithm is not defined in this case. Note that the assumption in Theorem 2.21 guarantees that $\nabla \varphi {\left(x\right)}^{\top }$ cannot be singular. But when the Jacobian matrix is singular at the solution point $x{}^{}\star$ , then Newton’s method looses its quadratic convergence property.

Computational Efficiency Finally, at each iteration, Newton’s method requires (i) that the Jacobian matrix $\nabla \varphi \left(x{}^{}k\right)$ be computed, which may be difficult and/or costly especially when the expression of $\varphi \left(x\right)$ is complicated, and (ii) that a linear system be solved. The analytic Jacobian can be replaced by a finite-difference approximation, yet this is costly as $n_x$ additional evaluations of $\varphi$ are required at each iteration. With the objective of reducing the computational effort, quasi-Newton methods generate an approximation of the Jacobian matrix, based on the information gathered from the iteration progress. To avoid solving a linear system for the search direction, variants of quasi-Newton methods also exist that generate an approximation of the inverse of the Jacobian matrix. Such methods will be described in 2.9.5 in the context of unconstrained optimization.

2.9.3 Unconstrained Optimization

We now turn to a description of basic techniques used for iteratively solving unconstrained problems of the form:

Many unconstrained optimization algorithms work along the same lines. Starting from an initial point, a direction of movement is determined according to a fixed rule, and then a move is made in that direction so that the objective function value is reduced; at the new point, a new direction is determined and the process is repeated. The main difference between these algorithms rest with the rule by which successive directions of movement are selected. A distinction is usually made between those algorithms which determine the search direction without using gradient information (gradient-free methods), and those using gradient (and higher-order derivatives) information (gradient-based methods). Here, we shall focus our attention on the latter class of methods, and more specifically on Newton-like algorithms.

Throughout this subsection, we shall assume that the objective function $f$ is twice continuously differentiable. By theorem 2.3, a necessary condition for $x{}^{}\star$ to be a local minimum of $f$ is $\nabla f\left(x{}^{}\star \right)=0$ . Hence, the idea is to devise an iterative scheme that finds a point satisfying the foregoing condition. Following the techniques discussed earlier in 2.9.2, this can be done by using a Newton-like algorithm, with $\varphi$ corresponding to the gradient $\nabla f$ of $f$ , and $\nabla {\varphi }^{}\top$ to its Hessian matrix ${\nabla }^{2}f$ .

At each iteration, a new iterate $x{}^{}k+1$ is obtained such that the linear approximation to the gradient at that point is zero,

The linear approximation yields the Newton search direction:

As discussed in 2.9.2, if it converges, Newton’s method exhibits a quadratic rate of convergence when the Hessian matrix ${\nabla }^{2}f\left(x{}^{}\star \right)$ is nonsingular at the solution point. However, since the Newton iteration is based on finding a zero of the gradient vector, there is no guarantee that the step will move towards a local minimum, rather than a saddle point or even a maximum. To preclude this, the Newton steps should be taken downhill, i.e., the following descent condition should be satisfied at each iteration,

Interestingly enough, with the Newton direction iteration, the descent condition becomes

That is, a sufficient condition to obtain a descent direction at $x{}^{}k$ is that the Hessian matrix ${\nabla }^{2}f\left(x{}^{}k\right)$ be positive definite. Moreover, if ${\nabla }^{2}f\left(x{}^{}\star \right)$ is positive definite at a local minimizer $x{}^{}\star$ of $f$ , then the Newton iteration converges to $x{}^{}\star$ when started sufficiently close to $x{}^{}\star$ . (Recall that, by Theorem 2.8, positive definiteness of ${\nabla }^{2}f\left(x{}^{}\star \right)$ is a sufficient condition for a local minimum of $f$ to be a strict local minimum.)

We now discuss two important improvements to Newton’s method, which are directly related to the issues discussed in Subsection 2.9.2, namely (i) the lack of global convergence, and (ii) computational efficiency.

2.9.4 Globalization Strategies

Up to this point, the development has focused on the application of Newton’s method. However, even in the simplest one-dimensional applications, Newton’s method has deficiencies (see, e.g., Example 2.20). Methods for correcting global convergence deficiencies are referred to as globalization strategies. It should be stressed than an efficient globalization strategy should only alter the iterates when a problem is encountered, but it should not impede the ultimate behavior of the method, i.e., the quadratic convergence of a Newton’s method should be retained.

In unconstrained optimization, one can detect problems in a very simple fashion, by monitoring whether the next iterate $x{}^{}k+1$ satisfies a descent condition with respect to the actual iterate $x{}^{}k$ , e.g., $f\left(x{}^{}k+1\right)<f\left(x{}^{}k\right)$ . Then, either one of two globalization strategies can be used to correct the Newton step. The first strategy, known as line search method, is to alter the magnitude of the step; the second one, known as trust region method, is to modify both the step magnitude and direction. We shall only concentrate on the former class of globalization strategies subsequently.

A line search method proceeds by replacing the full Newton step $x{}^{}k+1=x{}^{}k+\alpha d{}^{}k$ with

where the step-length $\alpha \ge 0$ is chosen such that the objective function is reduced,

Clearly, the optimal choice for $\alpha$ would be ${\alpha }^{\star }=\mathrm{\text{argmin}}\{l\left(\alpha \right)=f(x{}^{}k+\alpha d{}^{}k\alpha \ge 0\}$ . The resulting minimization problem can be solved by any one-dimensional exact minimization technique (e.g., Newton’s method itself). However, such techniques are costly in the sense that they often require many iterations to converge and, therefore, many function (or even gradient) evaluations.

In response to this, most modern algorithms implement so-called inexact line search criteria, which aim to find a step-length giving an ”acceptable” decrease in the objective function. Note that sacrificing accuracy, we might impair the convergence of the overall algorithm that iteratively employs such a line search. However, by adopting a line search that guarantees a sufficient degree of descent in the objective function, the convergence of the overall algorithm can still be established.

We now describe one popular definition of an acceptable step-length known as Armijo’s rule. Other popular approaches are the quadratic and cubic fit techniques, as well as Wolfe’s and Glodstein’s tests. Armijo’s rule is driven by two parameters $0<{\kappa }_1<1$ and ${\kappa }_2>1$ , which respectively manage the acceptable step-length from being too large or too small (Typical vales are ${\kappa }_1=0.2$ and ${\kappa }_2=2$ ). Define the line search function $l\left(\alpha \right):=f(x{}^{}k+\alpha d{}^{}k)$ , for $\alpha \ge 0$ , and consider the modified first-order approximation $\widehat{l}\left(\alpha \right)=l\left(0\right)+{\kappa }_1\alpha l^{\prime }\left(0\right)$ . A step-length $\overline{\alpha }\in (0,1)$ is deemed acceptable if the following conditions hold:

The first condition prevents the step-length $\overline{\alpha }$ from being too large, whereas the second condition prevents $\overline{\alpha }$ from being too small. The acceptable region defined by the Armijo’s rule is shown in Fig. 1.16. below.

2.9.5 Recursive Updates

Another limitation of Newton’s method when applied to unconstrained optimization problems is that the Hessian matrix of the objective function is needed at each iteration, then a linear system must be solved for obtaining the search direction. For many applications, this can be a costly computational burden. In response to this, quasi-Newton methods attempt to construct this information recursively. However, by so doing, the quadratic rate of convergence is lost.

The basic idea for many quasi-Newton methods is that two successive iterates $x{}^{}k,x{}^{}k+1$ together with the corresponding gradients $\nabla f\left(x{}^{}k\right),\nabla f\left(x{}^{}k+1\right)$ , yield curvature information by means of the first-order approximation relation

with $d{}^{}k=x{}^{}k+1-x{}^{}k$ . In particular, given $n_x$ linearly independent iteration increments ${\delta }^{}0,\dots ,{\delta }^{}{n}_x-1$ an approximation of the Hessian matrix can be obtained as

or for the inverse Hessian matrix as

where ${\gamma }^{}k=\nabla f\left(x{}^{}k+1\right)-\nabla f\left(x{}^{}k\right)$

Note that when the objective function is quadratic, the previous relations are exact. Many interesting quasi-Newton methods use similar ways, although more sophisticated, to construct an approximate Hessian matrix $B^k$ that progressively approaches the inverse Hessian. One of the most popular class of quasi-Newton methods (known as the Broyden family) proceeds as follows:

where $0\le \xi \le 1$ . It is easily seen that when supplemented with a line search strategy, $d{}^{}{k}^{\top }{\gamma }^{}k<0$ at each $k$ , and hence the Hessian matrix approximations are guaranteed to exist. Moreover, it can be shown that the successive approximates remain positive definite provided that $B^0$ is itself positive definite.

By setting $\xi =0$ , the previous equation yields the Davidon-Fletcher-Powell (DFP) method, which is historically the first quasi-Newton method, while setting $\xi =1$ gives the Broyden-Fletcher-Goldfard-Shanno (BFGS) method, for which there is substantial evidence that it is the best general purpose quasi-Newton method currently known.

2.9.6 Summary

A Newton-like algorithm including both a line search method (Armijo’s rule) and Hessian recursive update (DFP update) is as follows:

The standard unconstrained optimization algorithm in MATLAB (i.e. the function fminunc) is an implementation of quasi-Newton’s method, with DFP or BFGS update, and a line search strategy.

2.9.7 Constrained Nonlinear Optimization

In this subsection, we turn our attention to algorithms for iteratively solving constrained problems of the form:

Many modern deterministic algorithms for constrained NLP problems are based on the (rather natural) principle that, instead of solving a difficult problem directly, one had better solve a sequence of simpler, but related, subproblems, which converges to a solution of the original problem either in a finite number of steps or in the limit. Working along these lines, two classes of algorithms can be distinguished for solution of NLP problems with equality and/or inequality constraints. On the one hand, penalty function and interior-point methods consist of solving the problem as a sequence of unconstrained problems (or problems with simple constraints), so that algorithms for unconstrained optimization can be used. These methods, which do not rely on the KKT theory described earlier in 2.6 through 2.8, shall be briefly presented in 2.9.8 and 2.9.9. On the other hand, Newton-like methods solve NLP problems by attempting to find a point satisfying the necessary conditions of optimality (KKT conditions in general). Sequential quadratic programming (SQP), which shall be presented in 2.9.10, represents one such class of methods.

2.9.8 Penalty Function Methods

Methods using penalty functions transform a constrained problem into a single unconstrained problem or a sequence of unconstrained problems. This is done by placing the constraints into the objective function via a penalty parameter in a way that penalizes any violation of the constraints. To illustrate it, consider the NLP problem:

where $X\subset {ℝ}^{n_x}$ , $x\in {ℝ}^{n_x}$ and $f:X\to ℝ$ , $g:X\to {ℝ}^{n_g}$ and $h:X\to {ℝ}^{n_h}$ are defined on $X$ .

In general, a suitable penalty function $\alpha \left(x\right)$ for the minimization problem is defined by:

where $\varphi$ and $\psi$ are continuous functions satisfying the conditions:

Typically, $\varphi$ and $\psi$ are of the forms

with $p$ a positive integer (taking $p\ge 2$ provides continuously differentiable penalty functions). The function $f\left(x\right)+\mu \alpha \left(x\right)$ is referred to as the auxiliary function.

The conclusions of Example 2.22 that the solution of the penalty problem can be made arbitrarily close to the solution of the original problem, and the optimal auxiliary function value arbitrarily close to the optimal objective value, by choosing sufficiently large, is formalized in the following:

Note that the assumption that $X$ is compact is necessary, for it possible that the optimal objective values of the original and penalty problems are not equal otherwise. Yet, this assumption is not very restrictive in most practical cases as the variables usually lie between finite upper and lower bounds. Note also that no restriction is imposed on $f$ , $g$ and $h$ other than continuity. However, the application of an efficient solution procedure for the (unconstrained) auxiliary problems may impose additional restriction on these functions (see 2.9.3).

Under the conditions that (i) $f$ , $g$ , $h$ and $\varphi$ , $\psi$ are continuously differentiable, and (ii) $\overline{x}$ is a regular point the solution to the penalty problem can be used to recover the Lagrange multipliers associated with the constraints at optimality. In the particular case where $X={ℝ}^{n_x}$ , we get

The larger $\mu$ , the better the approximation of the Lagrange multipliers:

From a computational viewpoint, superlinear convergence rates might be achievable, in principle, by applying Newton’s method (or its variants such as quasi-Newton methods). Yet, one can expect ill-conditioning problems when $\mu$ is taken very large in the penalty problem. With a large $\mu$ , more emphasis is placed on feasibility, and most procedures for unconstrained optimization will move quickly towards a feasible point. Even though this point may be far from the optimum, both slow convergence and premature termination can occur due to very small step size and finite precision computations (round-off errors).

As a result of the above mentioned difficulties associated with large penalty parameters, most algorithms using penalty functions employ a sequence of increasing penalty parameters. With each new value of the penalty parameter, an optimization technique is employed, starting with the optimal solution corresponding to the previously chosen parameters value. Such an approach is often referred to as sequential unconstrained minimization (SUM) technique. This way, a sequence of infeasible points is typically generated, whose limit is an optimal solution to the original problem (hence the term exterior penalty function approach).

To conclude our discussion on the penalty function approach, we give an algorithm to solve a general NLP problem, where the penalty function used are of the form specified by $\varphi$ and $\psi$

2.9.9 Interior-Point Methods

Similar to penalty functions, barrier functions can also be used to transform a constrained problem into an unconstrained problem (or into a sequence of unconstrained problems). These functions act as a barrier and prevent the iterates from leaving the feasible region. If the optimal solution occurs at the boundary of the feasible domain, the procedure moves from the interior to the boundary of the domain, hence the name interior-point methods. To illustrate these methods, consider the NLP problem:

where $X$ is a subset of ${ℝ}^{n_x}$ , and $f:X\to ℝ$ , $g:X\to {ℝ}^{n_g}$ are continuous on ${ℝ}^{n_x}$ . Note that equality constraints, if any, should be accommodated within the set $X$ . (In the case of affine equality constraints, one can possibly eliminate them after solving for some variables in terms of the others, thereby reducing the dimension of the problem.) The reason why this treatment is necessary is because barrier function methods require the set $\{x\in {ℝ}^{n_x}:g(x)<0\}$ to be nonempty; this would obviously be not possible if the equality constraints $h\left(x\right)=0$ were accommodated within the set of inequalities as $h\left(x\right)\le 0$ and $h\left(x\right)\ge 0$ .

A barrier problem formulates as:

where $𝜃\left(\mu \right):=\inf \left\{f\right(x)+\mu b(x):g(x)<0,x\in X\}$ . Ideally, the barrier function $b$ should take value zero on the region $\{x:g(x)\le 0\}$ , and value $\infty$ on its boundary. This would guarantee that the iterates do not leave the domain $\{x:g(x)\le 0\}$ provided the minimization problem started at an interior point. However, this discontinuity poses serious difficulties for any computational procedure. Therefore, this ideal construction of $b$ is replaced by the more realistic requirement that $b$ be nonnegative and continuous over the region $\{x:g(x)\le 0\}$ and approach infinity as the boundary is approached from the interior:

where $\varphi$ is a continuous function over $\{z:z<0\}$ that satisfies the conditions

In particular, $\mu b$ approaches the ideal barrier function described above as $\mu$ approaches zero.

Typically barrier functions are

The following barrier function, known as Frisch’s logarithmic barrier function, is also widely used

Actually the Frisch’s logarithmic barrier $\varphi \left(z\right)=-\ln (-z)$ does not satisfy the nonnegativity requirement for $z<-1$ . However the requirement on $\varphi$ can be relaxed and it is sufficient that $\varphi$ be positive close to $z=0$ .

The function $f\left(x\right)+\mu b\left(x\right)$ is referred to as the auxiliary function.

Given $\mu >0$ , evaluating $𝜃\left(\mu \right)=\mathrm{\text{inf}}\left\{f\right(x)+\mu b(x):g(x)<0,x\in X\}$ seems no simpler than solving the original problem because of the constraint $g\left(x\right)<0$ . However, starting the optimization from a point in the region $S:=\{x:g(x)<0\}\cap X$ yields an optimal point in $S$ , even when the constraint $g\left(x\right)<0$ is ignored. This is because $b$ approaches infinity as the iterates approach the boundary of $\{x:g(x)\le 0\}$ from within $S$ , hence preventing them from leaving the set $S$ . This is formalized in the following:

Under the conditions that (i) $f$ , $g$ and $\varphi$ are continuously differentiable, and (ii) $x{}^{}\star$ is a regular point, the solution to the barrier function problem can be used to recover the Lagrange multipliers associated with the constraints at optimality. In the particular case where $X={ℝ}^{n_x}$ , we get:

The approximation of the Lagrange multipliers, gets better as $\mu$ gets closer to $0$ ,

The use of barrier functions for solving constrained NLP problems also faces several computational difficulties. First, the search must start with a point $x\in X$ such that $g\left(x\right)<0$ , and finding such a point may not be an easy task for some problems. Also, because of the structure of the barrier function $b$ , and for small values of the parameter $\mu$ , most search techniques may face serious ill-conditioning and difficulties with round-off errors while solving the problem to minimize $f\left(x\right)+\mu b\left(x\right)$ over $x\in X$ , especially as the boundary of the region $\{x:g(x)\le 0\}$ is approached. Accordingly, interior-point algorithms employ a sequence of decreasing penalty parameters $\left\{{\mu }^k\right\}\to 0$ as $k\to \infty$ ; with each new value ${\mu }^k$ , an optimal solution to the barrier problem is sought by starting from the previous optimal solution. As in the exterior penalty function approach, it is highly recommended to use suitable second-order Newton or quasi-Newton methods for solving the successive barrier problems.

We describe below a scheme using barrier functions for optimizing a NLP problem.

Note that although the constraint $g\left(x\right)<0$ may be ignored, it is considered in the problem formulation as most line search methods use discrete steps, and a step could lead to a point outside the feasible region (where the value of the barrier function is a large negative number), when close to the boundary. Therefore, the problem can effectively be treated as an unconstrained optimization problem only if an explicit check for feasibility is made.

In recent years, there has been much excitement because some variants of the interior-point algorithm can be shown to be polynomial in time for many classes of convex programs. Moreover, interior-point codes are now proving to be highly competitive with codes based on other algorithms, such as SQP algorithms presented subsequently. A number of free and commercial interior-point solvers are given in Tab. 1.1. below.

2.9.10 Sequential Quadratic Programming

Sequential quadratic programming (SQP) methods, also known as successive, or recursive, quadratic programming, employ Newton’s method (or quasi-Newton methods) to directly solve the KKT conditions for the original problem. As a result, the accompanying subproblem turns out to be the minimization of a quadratic approximation to the Lagrangian function subject to a linear approximation to the constraints. Hence, this type of process is also known as a projected Lagrangian, or the Newton-Lagrange, approach. By its nature, this method produces both primal and dual (Lagrange multiplier) solutions.

To present the concept of SQP, consider first the equality constrained nonlinear program:

where $x\in {ℝ}^{n_x}$ , and $f:{ℝ}^{n_x}\to ℝ$ , $h:{ℝ}^{n_x}\to {ℝ}^{n_h}$ are twice continuously differentiable. We shall also assume throughout that the equality constraints are linearly independent at a solution $x{}^{}\star$ . (The extension for including inequality constraints is considered subsequently.)

By Theorem 2.13, the first-order necessary conditions of optimality for an equality constrained NLP require a primal solution $x{}^{}\star \in {ℝ}^{n_x}$ and a Lagrange multiplier vector ${\lambda }^{}\star \in {ℝ}^{n_h}$ such that:

where $\mathcal{ℒ}(x,\lambda ):=f\left(x\right)+{\lambda }^{}\top h\left(x\right)$ . Now, consider a Newton-like method to solve the latter system of $n_x+n_h$ nonlinear equations. Given an iterate $(x{}^{}k,{\lambda }^{}k)$ , a new iterate $(x{}^{}k+1,{\lambda }^{}k+1)$ is obtained by solving the first-order approximation:

Substituting ${\nabla }_x\mathcal{ℒ}(x{}^{}k,{\lambda }^{}k)=\nabla f\left(x{}^{}k\right)+\nabla h\left(x{}^{}k\right){\lambda }^{}k$ , the first row of the system of equations is:

The Lagrange multiplier ${\lambda }^{}k$ can be simplified and thus:

The second row of the system of equations is independent of the Lagrange multipliers, it is:

Therefore the system can be simplified as:

Where we have defined the Newton step on $x$ as $d{}^{}k=x{}^{}k+1-x{}^{}k$ . The simplified system can be solved for $(d{}^{}k,{\lambda }^{}k)$ if a solution exists. Setting $x{}^{}k+1:=x{}^{}k+d{}^{}k$ , and incrementing $k$ by $1$ , we can then repeat the process until $d{}^{}k=0$ happens to solve the linear system. When this occurs, if at all, we shall have found a stationary point to the equality constrained NLP.

Interestingly enough, a quadratic programming (QP) minimization subproblem can be employed in lieu of the foregoing linear system to find any optimal solution for the linearized NLP problem. Given $(x{}^{}k,{\lambda }^{}k)$ we define the subproblem $\text{QP}(x{}^{}k,{\lambda }^{}k)$ as:

First, note that an optimum $d{}^{}k$ to $\text{QP}(x{}^{}k,{\lambda }^{}k)$ , if it exists, together with the set of Lagrange multipliers ${\lambda }^{}k+1$ associated with the linearized constraints satisfies the linearized KKT conditions. Second, the objective function of $\text{QP}(x{}^{}k,{\lambda }^{}k)$ represents not just a quadratic approximation for $f\left(x\right)$ but also incorporates an additional term $\frac{1}{2}d{}^{}{k}^{\top }\{{\sum }_{i=1}^{n_h}{\lambda }_i^k{\nabla }^2{h}_i\left(x{}^{}k\right)\}d{}^{}k$ to represent the curvature of the constraints.