\section{Sequential Quadratic Programming}\label{S:SequentialQuadraticProgramming}	
\setcounter{progress_ck_counter}{0}

\noi \emph{Sequential Quadratic Programming}\index{Sequential Quadratic Programming Technique (SQPT)}
is a powerful, iterative technique for solving nonlinear programming problems, well-suited for those
whose constraints are nonlinear.

\subsection{Method Derivation for Equality-type Constraints}
\label{SS:SQPEqualityDerivation}\index{Sequential Quadratic Programming Technique (SQPT)!equality-type constraints}

\noi An important connection exists between Newton's method and sequential quadratic programming\index{Newton's Method!role in sequential quadratic programming}. To best understand
the nature of this connection, we first develop
the technique in the case of equality-type
constraints. Thus, we focus on the problem

\begin{align}\label{E:SQPEqualityConstraints}
\mbox{ minimize } & f(\vecx )  \\
	\mbox{ subject to} \notag\\
   h_1(\vecx) &= 0 \notag\\
	h_2(\vecx) &= 0  \notag\\
	 \vdots & \notag \\
     h_p(\vecx) &= 0, \notag
	\end{align}
where $f$ and each $h_j$, $1 \leq j \leq p$, is a twice-differentiable function having domain
$S \subseteq \R^n$. In compact form, (\ref{E:SQPEqualityConstraints}) becomes

\begin{align}\label{E:SQPEqualityConstraintsVectorForm}
 \mbox{ minimize }  & f(\vecx )  \\
	\mbox{ subject to} & \notag\\
\vech(\vecx) &= \zerovec, \notag
	\end{align}	
where $\vech: S \to \R^p$ is the vector-valued function defined by
$\ds \vech(\vecx)= \begin{bmatrix} h_1(\vecx) \\ h_2(\vecx) \\ \vdots \\ h_p(\vecx) \end{bmatrix}$.

\noi
In a nutshell, to solve (\ref{E:SQPEqualityConstraintsVectorForm}), we will apply Newton's
Method to the associated Lagrangian function. At each iteration, we solve a quadratic programming
problem. The matrix associated with the objective function is the Hessian of the Lagrangian. Constraints are
formulated using the gradient of $L$, the Jacobian of $\vech$, and $\vech$ itself.

\newpar
The Lagrangian function, $L$, associated with (\ref{E:SQPEqualityConstraintsVectorForm}) is given by
\begin{equation}\label{E:SQPEqualityConstraintsLagrangian}
L(\vecx, \muvec)=f(\vecx)+\muvec^t \vech(\vecx),
\end{equation}
where $\vecx \in S$ and $\muvec \in \R^p$.

\newpar
Recall results (\ref{E:LagrangianCP1}) and (\ref{E:LagrangianCP1}) from
Section \ref{S:SaddlePoints}. These demonstrated that if
$\vecx_0$ is a KKT point of NLP (\ref{E:SQPEqualityConstraintsVectorForm}) with corresponding multiplier vector, $\muvec_0$, then the gradient of $L$ with respect
to both vector variables, $\vecx$ and $\muvec$, vanishes at $(\vecx_0,\muvec_0)$.
In other words, $(\vecx_0,\muvec_0)$ is a critical point of $L$.

\newpar
The nature of this critical point provides information useful for solving
(\ref{E:SQPEqualityConstraintsVectorForm}). If this critical point is a minimum
of $L$ on $S \times \R^p$, then for all feasible $\vecx$ in $S$,

\begin{align}\label{E:SQPLagrangianMinimum}
 f(\vecx_0)&= f(\vecx_0)+\muvec_0^t \vech(\vecx_0) \\
 &= L(\vecx_0,\muvec_0) \notag \\
 &\leq L(\vecx,\muvec_0) \notag \\
 &=f(\vecx)+\muvec_0^t \vech(\vecx) \notag \\
 &=f(\vecx).\notag
\end{align}	


\noi Hence, $\vecx_0$ is a solution of (\ref{E:SQPEqualityConstraintsVectorForm}). For many examples, such as the \emph{ConPro Manufacturing Company} NLP,
$(\vecx_0,\muvec_0)$ is not a minimum of $L$, but instead corresponds to a saddle point,
as defined in Section \ref{S:SaddlePoints}. Fortunately, Theorem
\ref{Th:SaddlePointsMainTheorem} tells us that this is also a sufficient condition
for $\vecx_0$ to be a solution of (\ref{E:SQPEqualityConstraintsVectorForm}).

\newpar
One means of estimating the critical point, $(\vecx_0,\muvec_0)$, is to
apply Newton's Method to $L$. Suppose
$\ds \vecw_1=\begin{bmatrix} \vecx_1 \\ \muvec_1 \end{bmatrix}$ is our initial
value, where $\vecx_1$ belongs to $S$ and
and $\muvec_1$ belongs to $\R^p$. Denote the Newton direction
of $L$\index{Lagrangian function!Newton direction} at $\vecw_1$ by
\begin{equation}\label{E:SQPNewtonDirection}
\Delta \vecw = \begin{bmatrix} \Delta \vecx \\ \Delta \muvec \end{bmatrix}.
\end{equation}

\newpar
Determining $\Delta \vecw$ requires evaluating at $\vecw_1$ the gradient and Hessian of $L$
in {\underbar{both}} vector variables, $\vecx$ and $\muvec$. Let $\nabla L$ and
$H_L$ denote these two quantities, where rows and columns are labeled so that differentiation takes place with respect
to the components of $\vecx$ first. Then straightforward calculations establish

\begin{equation}\label{E:SQPLagrangianGradient}
\nabla L(\vecx_1,\muvec_1) =
\begin{bmatrix} \nabla f(\vecx_1) + \jacobian_{\vech}(\vecx_1)^t\muvec_1 \\
\vech(\vecx_1)
\end{bmatrix},
\end{equation}
and

\begin{equation}\label{E:SQPLagrangianHessian}
H_L(\vecx_1,\muvec_1) =
\begin{bmatrix}
H_{L,\vecx}(\vecx_1,\muvec_1)  & \jacobian_{\vech}(\vecx_1)^t \\
\jacobian_{\vech}(\vecx_1) & 0_{p \times p}
\end{bmatrix}.
\end{equation}

\noi Here we have used $H_{L,\vecx}$ to denote the Hessian of $L$ in $\vecx$ alone. In terms of $f$
and $\vech$,
\begin{equation*}
H_{L,\vecx}(\vecx_1,\muvec_1)=
H_f(\vecx_1) + \sum_{j=1}^p \mu_{1,j} H_{h_j}(\vecx_1),
\end{equation*}
where $\mu_{1,j}$ denotes the $j$th component of $\muvec_1$.

\newpar Using the Newton direction formula, Equation (\ref{E:unconstrained_numeric_Newton_direction}) from Section \ref{S:NewtonsMethod}, we obtain

\begin{equation}\label{E:SQPNewtonDirection}
-\begin{bmatrix}
H_{L,\vecx}(\vecx_1,\muvec_1)  & \jacobian_{\vech}(\vecx_1)^t \\
\jacobian_{\vech}(\vecx_1) & 0_{p \times p}
\end{bmatrix}
\begin{bmatrix} \Delta \vecx \\ \Delta \muvec \end{bmatrix}
=\begin{bmatrix} \nabla f(\vecx_1) +  \jacobian_{\vech}(\vecx_1)^t \muvec_1 \\
\vech(\vecx_1)
\end{bmatrix}.
\end{equation}

\newpar
Equation (\ref{E:SQPNewtonDirection}) bears a resemblance to Equation
(\ref{E:QPKKTPartionedMatrix}) from Section \ref{S:QP}, which arose in the
 context of quadratic programming. To see this connection more clearly,
 we will associate to (\ref{E:SQPEqualityConstraintsVectorForm})
 the following \emph{quadratic subproblem}\index{quadratic subproblem}:

\begin{align}\label{E:SQPNewtonDirectionQPFormulation}
 \mbox{ minimize } \mbox{  } & \frac{1}{2} (\Delta \vecx)^t H_{L,\vecx}(\vecx_1,\muvec_1) \Delta \vecx +(\nabla f(\vecx_1) + \jacobian_{\vech}(\vecx_1)^t\muvec_1)^t(\Delta \vecx)  \\
	\mbox{ subject to  } &
	\jacobian_{\vech}(\vecx_1)(\Delta \vecx) = -\vech(\vecx_1). \notag
	\end{align}	

\noi
For $\Delta \vecx$ to be a KKT point of (\ref{E:SQPNewtonDirectionQPFormulation}) with
associated multiplier vector, $\Delta \muvec$, the feasibility condition must be
satisfied, i.e., $\ds \jacobian_{\vech}(\vecx_1)(\Delta \vecx) = -\vech(\vecx_1)$, and the
gradient of the associated Lagrangian must vanish, or, in other words,
\begin{equation*}
H_{L,\vecx}(\vecx_1,\muvec_1) \Delta \vecx +
(\nabla f(\vecx_1) + \jacobian_{\vech}(\vecx_1)^t\muvec_1)
+\jacobian_{\vech}(\vecx_1)^t \Delta \muvec = \zerovec.
\end{equation*}
\noi
But these conditions combine to yield matrix Equation (\ref{E:SQPNewtonDirection}).

\newpar
Thus, the values of $\Delta \vecx$ and $\Delta \muvec$ can be obtained
by solving either the quadratic subproblem (\ref{E:SQPNewtonDirectionQPFormulation}), or
the matrix equation (\ref{E:SQPNewtonDirection}). Using these values, we obtain
the Newton direction of $L$ at $\ds \vecw_1=\begin{bmatrix} \vecx_1 \\ \muvec_1 \end{bmatrix}$.
Namely, $\ds \Delta \vecw = \begin{bmatrix} \Delta \vecx \\ \Delta \muvec \end{bmatrix}$.
Therefore, $\vecw_2=\vecw_1 + \Delta \vecw$ so that
$\vecx_2=\vecx_1+\Delta \vecx$ and $\muvec_2=\muvec_1+\Delta \muvec$. This completes
the first iteration of
the \emph{Sequential Quadratic Programming Technique}.
We now summarize this process.

\newpar
{\underbar{The Sequential Programming Technique (SQPT) for NLPs with}}
\newline \underbar{{Equality-type Constraints}}\index{Sequential Quadratic Programming Technique (SQPT)!summary of steps} \newline
To obtain an approximate solution of NLP (\ref{E:SQPEqualityConstraintsVectorForm}) under
the assumption that $f$ and each component of $\vech$ are twice-differentiable on $S$:

\be
\item Starting with $\vecx_1$ in $S$, along with $\muvec_1$ in $\R^p$, use
the solution of (\ref{E:SQPNewtonDirectionQPFormulation}), which is
also the solution of
(\ref{E:SQPNewtonDirection}), to compute $\Delta \vecx$ and $\Delta \muvec$.
Of course, successfully doing so
hinges on whether
$\ds \begin{bmatrix}
H_{L,\vecx}(\vecx_1,\muvec_1)  & \jacobian_{\vech}(\vecx_1)^t \\
\jacobian_{\vech}(\vecx_1) & 0_{p \times p}
\end{bmatrix}$ is invertible.
\item Set
$\vecx_2=\vecx_1+\Delta \vecx$ and $\muvec_2=\muvec_1+\Delta \muvec$.
\item Return to (1), replace $\vecx_1$ and $\muvec_1$
with $\vecx_2$ and $\muvec_2$, respectively, and repeat the process.
\item Terminate the process when the difference between successive
iterates, $\vecx_{k+1}-\vecx_k$ is smaller than some specified tolerance.
\ee

\noi
Under ideal circumstances, the sequence of iterates,
\begin{equation*}
\left\{ (\vecx_1,\muvec_1), (\vecx_2,\muvec_2),\ldots \right\}
\end{equation*}

\noi converges to $\ds (\vecx_0,\muvec_0)$, where $\vecx_0$ is a solution
of (\ref{E:SQPEqualityConstraintsVectorForm}) having associated
multiplier vector, $\muvec_0$.

\newpar
It is important to recognize that the SQPT, for the case of equality constraints.
generates a sequence of iterates identical to that obtained by applying Newton's Method to the associated Lagrangian function.

\newpar
To demonstrate this new technique, we consider the NLP

\begin{align}\label{E:QPBasic Example}
 \mbox{minimize} f(x_1,x_2 )&= \sin(x_1-x_2) \\
	\mbox{ subject to} & \notag\\
	x_1^2+x_2^2 &=1 \notag
	\end{align}	

\noi The solution of (\ref{E:QPBasic Example}) can of course be obtained by solving
the constraint for one decision variable, substituting the result into the objective, and minimizing
the resulting function of a single variable. Doing so leads to a solution of
$\ds \vecx_0=\begin{bmatrix} -1/\sqrt{2} \\ 1/\sqrt{2} \end{bmatrix}$ with corresponding objective value
$f(\vecx_0)=-\sin(\sqrt{2})$.

\newpar
While NLP (\ref{E:QPBasic Example}) is a very simple example, it is also well suited for demonstrating our new technique, due to the nature of the objective function and the nonlinear constraint.

\newpar
We use as our initial value $\ds \vecx_1=\begin{bmatrix} 1 \\ 2 \end{bmatrix} $
and $\muvec_1=1$. (For the sake of consistency with our previous derivations, we
express $\muvec_1$ and $\vech$ using vector notation even though both are scalar-valued for this particular problem.) The quantities needed to compute $\Delta \vecx$ and $\Delta \muvec$
consist of the following:

\be
\item $\ds \vech(\vecx)=x_1^2+x_2^2-1$
\item $\ds \nabla f(\vecx)=\begin{bmatrix} \cos(x_1-x_2) \\ -\cos(x_1-x_2) \end{bmatrix}$
\item $\ds \jacobian_{\vech}(\vecx) = \begin{bmatrix} 2x_1 & 2x_2 \end{bmatrix}$
\item $\ds L(\vecx,\muvec)= \sin(x_1-x_2) + \muvec (x_1^2+x_2^2-1)$
\item $\ds H_{L,\vecx}(\vecx,\muvec)=
\begin{bmatrix} -\sin(x_1-x_2) +2\muvec & \sin(x_1-x_2) \\
\sin(x_1-x_2) & -\sin(x_1-x_2) +2\muvec \end{bmatrix}$
\ee

\noi
Evaluation of these quantities at $\ds \vecx_1=\begin{bmatrix} 1 \\ 2 \end{bmatrix}$
and $\muvec_1=1$ and substitution of the results into (\ref{E:SQPNewtonDirection})
yields the following matrix equation:

\begin{equation}\label{E:SQPBasicExampleMatrixEquation}
-\begin{bmatrix} \sin(1)+2 & -\sin(1) & 2 \\
-\sin(1) & \sin(1)+2 & 4 \\
2 & 4 & 0 \end{bmatrix}
\begin{bmatrix} \Delta \vecx \\ \Delta \muvec \end{bmatrix}
= \begin{bmatrix} \cos(1)+2 \\ -\cos(1)+4 \\ 4 \end{bmatrix}.
\end{equation}

\noi Solving (\ref{E:SQPBasicExampleMatrixEquation}), we obtain
$\ds \begin{bmatrix} \Delta \vecx \\ \Delta \muvec \end{bmatrix}
\approx \begin{bmatrix} -.7 \\ -.650 \\ -.550 \end{bmatrix}$, from
which it follows that
\begin{align*}
\vecx_2 &\approx \vecx_1+\Delta \vecx \\
&= \begin{bmatrix} .3 \\ 1.35 \end{bmatrix}
\end{align*}
and $\muvec_2=\muvec_1+\Delta \muvec \approx .45$.


\newpar
Table \ref{Ta:SQPIterations} provides values of $\vecx_k$ and $\muvec_k$ for the first
seven iterations of the SQPT applied
to NLP (\ref{E:QPBasic Example}) using $\ds \vecx_1=\begin{bmatrix} 1 \\ 2 \end{bmatrix} $
and $\muvec_1=1$.

\begin{table}[h]
\caption{Results of the Sequential Quadratic Programming Technique applied to NLP (\ref{E:QPBasic Example})}
\label{Ta:SQPIterations}
\begin{tabular}{|c|c|c|}
\hline
$k$ & $\vecx_k^t$ &      $\muvec_k$ \\ \hline
1 &   $[1,2]$  &     $1$ \\ \hline
2 & $[.3,1.35]$ & .45\\  \hline
3 & $[-.161,1.115]$ & 1.90\\  \hline
4 & $[-.474,.949]$ & .095\\  \hline
5 & $[-.678, .780]$ & .075\\  \hline
6 & $[-.730,.691]$ & .104\\  \hline
7 & $[-.706,.708]$ & .110\\  \hline
8 & $[-.707,.707]$ & .110\\  \hline
\end{tabular}
\end{table}


\noi The values of $\vecx_k$ clearly converge to the solution of (\ref{E:QPBasic Example}), namely
$\ds \vecx_0=\begin{bmatrix} -1/\sqrt{2} \\ 1/\sqrt{2} \end{bmatrix}$. But to what do
the values of $\muvec_k$ converge? The answer of course is the corresponding Lagrange multiplier.

\begin{progress_ck}
Verify that $\ds \vecx_0=\begin{bmatrix} -1/\sqrt{2} \\ 1/\sqrt{2} \end{bmatrix}$ is a KKT point
of (\ref{E:QPBasic Example}) with corresponding multiplier,
$\ds \muvec_0=\frac{\cos(\sqrt{2})}{\sqrt{2}} \approx .110$. Then use the \mf{NewtonsMethod} procedure from Section \ref{SS:NewtonsMethodMapleImplementation} to verify that the values
in Table \ref{Ta:SQPIterations} coincide with those obtained by applying Newton's
Method to the Lagrangian, $L$, using an initial
value
$\ds \begin{bmatrix} \vecx_1 \\ \muvec_1 \end{bmatrix}=\begin{bmatrix} 1 \\ 2 \\ 1 \end{bmatrix}$.
Note that the procedure
must be modified slightly to reflect the fact the Lagrangian is a function
of three variables, as opposed to two. The three consist of the two components of $\vecx$, along with $\muvec$.
\end{progress_ck}

\subsection{The Convergence Issue}\index{Sequential Quadratic Programming Technique (SQPT)!convergence issue}

\noi The SQPT is based upon Newton's Method, which, as
we discovered in Section \ref{SS:NewtonsMethodConvergenceIssues}, exhibits
a quadratic rate of convergence\index{Newton's Method!quadratic rate of convergence} to the desired
minimum, under appropriate conditions. Table \ref{Ta:SQPIterations} exhibits this fast rate of convergence. At the same time,
a numeric technique such as Newton's Method can generate a sequence of iterates converging to a critical point that is not a minimum.

\newpar
It should not come as a surprise then
that a sequence of SQPT iterates,
$\ds \left\{ (\vecx_1,\muvec_1), (\vecx_2,\muvec_2),\ldots \right\}$, can
have the same undesirable property. We can demonstrate this phenomenon by
applying the SQPT to
NLP (\ref{E:QPBasic Example}) using a different initial value.
If $\ds \vecx_1=\begin{bmatrix} -1 \\ -2 \end{bmatrix} $
and $\muvec_1=-1$, then the resulting sequences of iterates do not converge
to $\ds \vecx_0=\begin{bmatrix} -1/\sqrt{2} \\ 1/\sqrt{2} \end{bmatrix}$
and $\ds \muvec_0=\frac{\cos(\sqrt{2})}{\sqrt{2}}$, respectively. Instead, the sequence
$\left\{\vecx_k\right\}$ converges to the second KKT point of NLP (\ref{E:QPBasic Example}),
$\ds \vecx_0=\begin{bmatrix} 1/\sqrt{2} \\ -1/\sqrt{2} \end{bmatrix}$, and
the sequence $\left\{\muvec_k\right\}$ converges to the corresponding multiplier vector,
$\ds \muvec_0=-\frac{\cos(\sqrt{2})}{\sqrt{2}} \approx -.110$. Note that while $\vecx_0$
in this case is a KKT point, it is not a solution of (\ref{E:QPBasic Example}). Mere comparison
of objective values shows this to be the case.

\subsection{Inequality-Type Constraints}\index{Sequential Quadratic Programming Technique (SQPT)!inequality-type constraints}

\noi In the case where only equality-type constraints are present, solving
NLP (\ref{E:SQPEqualityConstraints}) using the SQPT reduces to
solving a sequence of quadratic programming problems of the form (\ref{E:SQPNewtonDirectionQPFormulation}), each of which has only equality-type constraints.

\newpar
If an NLP has inequality-type constraints as
well, a modification of the Sequential Programming Technique requires
us to solve at each iteration a quadratic programming subproblem where inequality-type
constraints are also present. Instead of formally deriving the exact quadratic programming
problem as we did in Section \ref{SS:SQPEqualityDerivation}, we instead state the end
result.

\newpar
Suppose we add to NLP (\ref{E:SQPEqualityConstraintsVectorForm}) inequality-type constraints, expressing it as
\begin{align}\label{E:SQPInequalityConstraintsVectorForm}
 \mbox{ minimize }  & f(\vecx )  \\
	\mbox{ subject to} & \notag\\
	\vecg(\vecx) &\leq \zerovec \notag\\
\vech(\vecx) &= \zerovec. \notag
	\end{align}	

\newpar
Here, $\vecg: S \to \R^m$ is the vector-valued function each of whose components are
twice-differentiable. The Lagrangian, $L$, is given by
\begin{equation*}
L(\vecx,\lambdavec,\muvec)=f(\vecx)+\lambdavec^t \vecg(\vecx)+\muvec^t \vech(\vecx),
\end{equation*}
whose Hessian in $\vecx$ we once again denote by $H_{L,\vecx}$.

\newpar
Choose $\vecx_1$ in $S$, $\lambdavec_1$ in $\R^m$ with $\lambdavec_1 \geq \zerovec$, and $\muvec_1$ in $\R^p$. Then set $\vecx_2=\vecx_1+\Delta \vecx$, where $\Delta \vecx$ is the solution the following:

 \begin{align}\label{E:SQPInequalityConstraintsQPFormulation}
 \mbox{ minimize } \mbox{  } & \frac{1}{2} (\Delta \vecx)^t H_{L,\vecx}(\vecx_1,\muvec_1) \Delta \vecx +\nabla f(\vecx_1)^t \Delta \vecx  \\
	\mbox{ subject to} & \notag\\
    &\jacobian_{\vecg}(\vecx_1)(\Delta \vecx) \leq -\vecg(\vecx_1) \notag \\
	&\jacobian_{\vech}(\vecx_1)(\Delta \vecx) = -\vech(\vecx_1). \notag
	\end{align}

\noi
The multiplier vectors $\lambdavec_2$ and $\muvec_2$ are assigned the values associated
with the solution of (\ref{E:SQPInequalityConstraintsQPFormulation}). Note that this process differs from the situation when only equality-type constraints were present. There, the multiplier vector was labeled
$\Delta \muvec$, which we then added to $\muvec_1$ to obtain $\muvec_2$.


\newpar
For example, consider the 3-variable NLP,

\begin{align}\label{E:QPBasic Example2}
 \mbox{minimize} f(x_1,x_2,x_3 )&= x_1^2+2x_1^2x_2x_3+x_2^2-x_3^2 \\
	\mbox{ subject to} & \notag\\
    x_1+2x_2-x_3 &\leq 2 \notag \\
   x_1^2 \phantom{x2,x2,x2} &\geq 1 \notag \\
	x_1^2+x_2^2+x_3^2 &=4. \notag
	\end{align}	


\noi We will use as initial values,
$\ds \vecx_1=\begin{bmatrix} 1 \\ 2 \\ 3 \end{bmatrix}$,
$\ds \lambdavec_1=\begin{bmatrix} 1 \\ 1 \end{bmatrix}$,
and $\muvec_1=-1$.

\newpar
Straightforward calculations lead to the following expressions:
\be
\item $\ds -\vecg(\vecx)=-\begin{bmatrix} x_1 + 2x_2 -x_3-2 \\ -x_1^2+1 \end{bmatrix}$
\item $\ds -\vech(\vecx)=-\left(x_1^2+x_2^2+x_3^2-4 \right)$
\item $\ds \nabla f(\vecx)=\begin{bmatrix} 2x_1+4x_1x_2x_3
\\ 2x_1^2x_3+2x_2 \\ 2x_1^2x_2-2x_3 \end{bmatrix}$
\item $\ds \jacobian_{\vecg}(\vecx) = \begin{bmatrix} 1 & 2 & -1 \\
0 & -2x_2 & 0 \end{bmatrix}$
\item $\ds \jacobian_{\vech}(\vecx) = \begin{bmatrix} 2x_1 & 2x_2 & 2x_3\end{bmatrix}$
\item \begin{align*}
L(\vecx,\lambdavec,\muvec) &=f(\vecx)+\lambdavec^t \vecg(\vecx)+\muvec^t \vech(\vecx) \\
&=x_1^2+2x_1^2x_2x_3+x_2^2-x_3^2 +
\lambda_1\left(x_1+2x_2-x_3 - 2\right) \\
+& \lambda_2\left( -x_2^2+1 \right) +
\mu_1 \left(x_1^2+x_2^2+x_3^2-4\right)
\end{align*}
\item $\ds H_{L,\vecx}(\vecx,\lambdavec,\muvec)=
\begin{bmatrix}
2+4x_2x_3+2\mu_1 & 4x_1x_3 & 4x_1x_2 \\
4x_1x_3 & 2-2\lambda_2+2\mu1 & 2x_1^2 \\
4x_1x_2 & 2x_1^2 & -2+2\mu_1
\end{bmatrix}$.
\ee

\newpar
Evaluating these quantities at $\vecx_1$, $\lambdavec_1$, $\muvec_1$, we obtain
\be
\item $\ds -\vecg(\vecx_1)=\begin{bmatrix} 0 \\ 3 \end{bmatrix}$
\item $\ds -\vech(\vecx_1)=-10$
\item $\ds \nabla f(\vecx_1)=\begin{bmatrix} 26 \\ 10 \\ -2 \end{bmatrix}$
\item $\ds \jacobian_{\vecg}(\vecx_1) = \begin{bmatrix} 1 & 2 & -1 \\
0 & -4 & 0 \end{bmatrix}$
\item $\ds \jacobian_{\vech}(\vecx_1) = \begin{bmatrix} 2 & 4 & 6 \end{bmatrix}$
\item $\ds H_{L,\vecx}(\vecx_1,\lambdavec_1,\muvec_1)=
\begin{bmatrix}
24 & 12 & 8 \\ 12 & -2 & 2 \\ 8 & 2 & -4
\end{bmatrix}$.
\ee

\newpar
Now we substitute these results into (\ref{E:SQPInequalityConstraintsQPFormulation})
and solve the resulting quadratic programming problem. To find the KKT point and corresponding multipliers, we could certainly construct the associated Lagrangian and use methods from
Section \ref{SS:SolvingConvexNLPs}. In an effort to streamline the solution process,
however, we will first determine the KKT point
using the matrix form of Maple's \mf{QPSolve} command. If we enter
the preceding quantities in Maple so that items 1, 2, 3, 4, 5, and 6 correspond
to \mf{d, b, p, C, A,} and \mf{Q}, respectively, then \mf{QPSolve([p,Q], [C,d,A,b])} returns a value
of $\ds \Delta \vecx \approx \begin{bmatrix} -.4000 \\ -.7500 \\ -1.0335 \end{bmatrix}$.

\newpar
While the \mf{QPSolve} command does not return the multiplier values, we can obtain them
by using $\Delta \vecx$ along with the given functions. Substituting
$\Delta \vecx$ into each side of the inequality
$\ds \jacobian_{\vecg}(\vecx_1)(\Delta \vecx) \leq -\vecg(\vecx_1)$ establishes
that the constraint corresponding to the second row of $\ds \jacobian_{\vecg}(\vecx_1)$
is binding. Thus $\lambda_{2,1}$, the first component of $\lambdavec_2$, equals zero by
complementary slackness.

\newpar
In the spirit of Section \ref{SS:QPInequalityConstraints},
let us define $\widetilde{C}$ to be the submatrix of $\ds \jacobian_{\vecg}(\vecx_1)$ formed
using the row corresponding to the second constraint, which is binding. In this example,
$\ds \widetilde{C}=\begin{bmatrix} 0 & -4 & 0 \end{bmatrix}$. Moreover, let
$\ds \widetilde{\vecd}$ be the subvector of
$\ds \vecg(\vecx_1)=\begin{bmatrix} 0 \\ 3 \end{bmatrix}$
corresponding to the second constraint, i.e.,  $\ds \widetilde{\vecd}=-3$.

\newpar
Now consider the matrix equation
\begin{equation}\label{QPBasiceExampleBorderedHessian}
\begin{bmatrix}
H_{L,\vecx}(\vecx_1,\lambdavec_1,\muvec_1) & \jacobian_{\vech}(\vecx_1)^t & \widetilde{C}^t \\
\jacobian_{\vech}(\vecx_1) & 0 & 0 \\
\widetilde{C} & 0 & 0
\end{bmatrix}
\begin{bmatrix} \Delta \vecx \\ \muvec_2 \\ \lambda_{2,2}
\end{bmatrix} =
\begin{bmatrix} -\nabla f(\vecx_1) \\ -\vech(\vecx_1) \\ -\widetilde{\vecd}
\end{bmatrix}.
\end{equation}

\noi
This equation must be satisfied by the triple $\ds \Delta \vecx , \muvec_2 , \lambda_{2,2}$
in order for (\ref{E:SQPInequalityConstraintsQPFormulation})
to have a solution $\Delta \vecx$ that yields a non-binding first
constraint, i.e., a constraint for which $\ds \lambda_{2,1}=0$.

\newpar
Since $\ds \lambda_{2,1}=0$, we use only the unknown value, $\lambda_{2,2}$
in the vector of unknowns. In this example, substitution of all known quantities
into Equation (\ref{QPBasiceExampleBorderedHessian}) yields

\begin{equation*}
\begin{bmatrix}
24 & 12 & 8 & 2 & 0 \\
12 & -2 & 2 & 4 & -4 \\
8 & 2 & -4 & 6 & 0  \\
2 & 4 & 6 & 0 & 0 \\
0 & -4 & 0 & 0 & 0
\end{bmatrix}
\begin{bmatrix} \Delta \vecx \\ \muvec_2 \\ \lambda_{2,2}
\end{bmatrix} =
\begin{bmatrix} -26 \\ -10 \\ 2 \\ -10 \\ 3
\end{bmatrix},
\end{equation*}
from which it follows that
$\ds \Delta \vecx \approx \begin{bmatrix} -.4 \\ -.75 \\ -1.0355
\end{bmatrix}$,
$\muvec_2 \approx .4268$, and $\lambda_{2,2} \approx 1.587$.

\newpar
Thus, one iteration of the technique is complete with the results
\newline $\ds \vecx_2=\vecx_1+\Delta \vecx \approx \begin{bmatrix} .6 \\ 1.25 \\ 1.9665
\end{bmatrix}$,
$\ds \lambdavec_2 \approx \begin{bmatrix} 0 \\ 1.587
\end{bmatrix}$, and $\muvec_2 \approx .4268$.

\newpar
At this point, we repeat the process, returning to
(\ref{E:SQPInequalityConstraintsQPFormulation}), replacing $\vecx_1$, $\lambdavec_1$,
and $\muvec_1$, with $\vecx_1$, $\lambdavec_1$,
and $\muvec_1$, respectively, and so on. The results after only the third iteration are
given by $\ds \vecx_4 \approx \begin{bmatrix} .003 \\ 1.00 \\ 1.735
\end{bmatrix}$,
$\ds \lambdavec_4 = \approx \begin{bmatrix} 0 \\ 1.995
\end{bmatrix}$, and $\muvec_4 \approx 1.011$. The exact values, which
can be determined by finding the KKT point and corresponding multipliers for
the original NLP (\ref{E:QPBasic Example2}), are given
by $\ds \vecx_0= \begin{bmatrix} 0 \\ 1 \\ \sqrt{3} \end{bmatrix}$,
$\ds \lambdavec_0 = \begin{bmatrix} 0 \\ 2
\end{bmatrix}$, and $\muvec_0 =1$. Thus, a mere three iterations in this case
produces quite accurate results.

\subsection{A Maple Implementation of the
Sequential \mbox{Quadratic} \mbox{Programming} Technique}\index{Sequential Quadratic Programming Technique (SQPT)!Maple implementation}

\noi The worksheet \textbf{Sequential Quadratic Programming Technique.mw}
provides a systematic means by
which to execute the SQPT. Here we use it
to complete the first iteration of NLP (\ref{E:QPBasic Example2})

\begin{maplelist}
\item \mf{with(VectorCalculus):with(LinearAlgebra):with(Optimization):}
\item \mf{f:=(x1,x2,x3)-> x1\^{}2+2x1\^{}2*x2*x3+x2\^{}2-x3\^{}2: \newline
\# Enter objective function.}
\item \mf{Delf:=unapply(Gradient(f(x1,x2,x3),[x1,x2,x3]),[x1,x2,x3]): \newline
\# Gradient of f.}
\item \mf{h:=(x1,x2,x3)-> <x1\^{}2+x2\^{}2+x3\^{}2-4>: \newline
\# Vector-valued form of equality constraint function.}
\item \mf{Jh:=unapply(Jacobian(h(x1,x2,x3),[x1,x2,x3]),[x1,x2,x3]): \newline
\# Jacobian function of h.}
\item \mf{g:=(x1,x2,x3)-><x1+2*x2-x3-2,-x\^{}2+1>: \newline
\# Vector-valued form of inequality constraint functions.}
\item \mf{Jg:=unapply(Jacobian(g(x1,x2,x3),[x1,x2,x3]),[x1,x2,x3]): \newline
\# Jacobian function of g.}
\item \mf{lambda:=<lambda1,lambda2>: \newline
\# Create vector of multipliers.}
\item \mf{L:=unapply(f(x1,x2,x3)+Transpose(lambda).g(x1,x2,x3),\newline
   +<mu>.h(x1,x2,x3), [x1,x2,x3,mu,lambda1,lambda2]):\newline
 \# Create Lagrangian Function.}
\item \mf{HLx:=unapply(Hessian(L(x1,x2,x3,mu,lambda1,lambda2),
\newline [x1,x2,x3]),[x1,x2,x3,mu,lambda1,lambda2]): 
\newline \# Hessian in x1,x2,x3 of the Lagrangian.}
\item \mf{X1:=1,2,3:mu1:=-1:Lambda1:=1,1: 
\newline \# Initial choices for variables and multipliers. Note 
\newline Lambda1 is case sensitive since lambda1 has already 
\newline been defined as a variable in the Lagrangian.}
\item \mf{Qsol:=QPSolve([Delf(X1),HLx(X1,mu1,Lambda1)],
\newline [Jg(X1),-g(X1),Jh(X1),-h(X1)]);\newline
\# Solve first quadratic subproblem in iteration process.}
\mo{Qsol:=\left[-8.1547, \begin{bmatrix} -.3994 \\ -.7500 \\ -1.0335 \end{bmatrix} \right]}
\item \mf{Jg(X1).<Qsol[2][1],Qsol[2][2],Qsol[2][3]>+g(X1); \newline
\# Determine which inequality constraints in solution of 
\newline subproblem are binding.}
\mo{(-.8659)\vece_x+(0)\vece_y}
\item \mf{$\sim$C:=SubMatrix(Jg(X1),[2],[1,2,3]); \newline
\# The second constraint is binding, so we form a submatrix
\newline of Jg(X1) using only the second row. Furthermore, we
\newline know the first component of Lambda2 is zero by
\newline complementary slackness.}
\mo{\sim C:=\begin{bmatrix} 0 & -4 & 0 \end{bmatrix}}
\item \mf{$\sim$d:=SubVector(g(X1),[2]); \newline
\# Use entry of Jg(X1) corresponding to binding constraint.}
\item \mf{B:=<<<HLx(X1,mu1,Lambda1)|Transpose(Jh(X1))|Transpose($\sim$C)>>, \newline
<Jh(X1)|ZeroMatrix(1,2)>,<$\sim$C|ZeroMatrix(1,2)>>;\newline
\# Create coefficient matrix for finding changes in x1, x2, 
\newline and x3, along with new multiplier values.}
\mo{B:=\begin{bmatrix}
24 & 12 & 8 & 2 & 0 \\
12 & -2 & 2 & 4 & -4 \\
8 & 2 & -4 & 6 & 0  \\
2 & 4 & 6 & 0 & 0 \\
0 & -4 & 0 & 0 & 0
\end{bmatrix}}
\item \mf{w:=evalf(MatrixInverse(B).<-Delf(X1),-h(X1),-$\sim$ d>);}
\mo{w:=\begin{bmatrix} -.3994 \\ -.7500 \\ -1.0335 \\ .4268 \\ 1.5869 \end{bmatrix}}
\item \mf{X2:=X1[1]+w[1,1],X1[2]+w[2,1],X1[3]+w[3,1];\newline
\# Create X2 using each entry of X1, added to the 
\newline corresponding entry in w.}
\mo{X2:=.6006,1.2500,1.9665}
\item \mf{mu2:=w[4,1];}
\mo{\mu2:=.4268}
\item \mf{Lambda2:=0,w[5,1];\newline
\# Lambda2 has its first entry equal to zero since the first 
\newline constraint in the quadratic subproblem was not binding at 
\newline the solution. We use the fifth entry of w to form the 
\newline second component of Lambda2. Now return to Step 1 and 
\newline repeat the process again starting at X2,mu2,Lambda2.}
\mo{\Lambda2=0,1.5869}
\end{maplelist}

\newpar
Subsequent iterations are handled in a similar manner.

\subsection{An Improved Version of the SQPT}\index{Sequential Quadratic Programming Technique (SQPT)!improved version with merit function (MSQPT)}

\noi
Numerous refinements of the SQPT exist, each intended to improve upon the original method. Perhaps the simplest of these stems from the fact that at each iteration of the SQPT,
the value of $\vecx_{k+1}$ is obtained by adding the entire quantity $\Delta \vecx$
to $\vecx_k$. While $\Delta \vecx$ is the solution of the quadratic
subproblem, the value $\vecx_{k+1}$ may be no less ``infeasible'' than was $\vecx_k$ if it
was infeasible to start with.
Consequently, instead of adding the entire quantity $\Delta \vecx$
to $\vecx_k$ in order to obtain $\vecx_{k+1}$, we should consider adding only a fraction
of its value.

\newpar
The simplest means for carrying out this process is to utilize a decision-making
rule, whereby at each iteration, we choose the fraction of $\Delta \vecx$ that
best decreases the objective value while simultaneously minimizing the ``infeasibility error.''
This task is accomplished by introducing a new function,

\begin{equation}\label{E:SQPMeritFunction}
M(\vecx)=f(\vecx)+\rho P(\vecx),
\end{equation}
where $P$ is a \emph{penalty function}\index{penalty function} that measures the extent of constraint violation and where $\rho$ is a fixed, large positive real number. The function, $M$, is referred to as
a \emph{merit function}\index{merit function}. Corresponding to NLP (\ref{E:SQPInequalityConstraintsVectorForm}), various different choices exist for
$P$, with one of the most common being

\begin{equation}\label{SQPPenaltyFunction}
P(\vecx)= \sum_{i=1}^m \max(g_i(\vecx),0)^2 +  \sum_{j=1}^p h_j(\vecx)^2.
\end{equation}

\noi
Here, \emph{max} denotes the maximum function,
\begin{equation}\label{SQPMaximumFunction}
\max(x,0)=\begin{cases}
		x,	&\text{if $x>0$};\\
		0,	&\text{otherwise.}
	\end{cases}
\end{equation}

\noi
Note that if $\vecx$ is feasible for (\ref{E:SQPInequalityConstraintsVectorForm}),
then $P(\vecx)=0$, and that the value of $P$ increases with the overall amount of constraint violation.

\newpar
With this choice of $P$, the merit function becomes
\begin{equation}\label{E:SQPMeritFunction}
M(\vecx)=f(\vecx)+\rho \left(\sum_{i=1}^m \max(g_i(\vecx),0)^2 +  \sum_{j=1}^p h_j(\vecx)^2\right).
\end{equation}

\newpar
Our decision-making rule is to scale $\Delta \vecx$ at each iteration by an amount, $t_0$, where
$t_0$ is the smallest positive local minimum of $\phi(t)=M(\vecx_k+t\Delta \vecx)$. Then
we set $\vecx_{k+1}=\vecx_k + t_0 \Delta \vecx$. However, we continue to determine $\lambdavec_{k+1}$
and $\muvec_{k+1}$ as we did for the original technique. The new method based upon an appropriate
scaling of $\Delta \vecx$, is known as the
\emph{Merit Function Sequential Programming Technique}, or MSQPT.

\newpar
To illustrate how this modification is an improvement upon the original technique, we revisit NLP
(\ref{E:QPBasic Example2}). Recall that
$\ds \vecx_1 =\begin{bmatrix} 1 \\ 2 \\ 3
\end{bmatrix}$
and that
at the end of the first iteration,
$\ds \Delta \vecx \approx \begin{bmatrix} -.4 \\ -.75 \\ -1.034
\end{bmatrix}$. Using $\rho=10$, the function $\phi$ becomes

\begin{align}\label{E:MSQPTBasicExampleMeritFunction}
\phi(t)&= M(\vecx_1+t\Delta \vecx) \\
&=f(1-.4t,2-.75t,3-1.034t) \notag \\
&+10\left(\sum_{i=1}^2 \max(g_i(1-.4t,2-.75t,3-1.034t),0)^2\right. \notag \\
&+\left. h(1-.4t,2-.75t,3-1.034t)^2 \right). \notag
\end{align}

\noi
Due to the piecewise nature of $P$, the formula for $\phi$ is quite complicated. However,
its graph, shown in Figure (\ref{Fi:SQPMeritFunction}), clearly indicates the desired
positive local minimum.

\begin{figure}[h]
\centering
\psfrag{taxislabel}[Bc][Bc]{$t$}
\includegraphics[width=.8\linewidth]{images/f09x05x01.eps}
\caption{Graph of $\phi(t)= M(\vecx_1+t\Delta \vecx)$. }
\label{Fi:SQPMeritFunction}
\end{figure}


\newpar
A single-variable implementation of Newton's Method demonstrates
that the minimum occurs at $t_0 \approx 1.3062$. Thus,
\begin{align}\label{E:MSQPTBasicExampleEndingValue}
\vecx_2 &= \vecx_1 + t_0 \Delta \vecx \\
&\approx \begin{bmatrix} 1 \\ 2 \\ 3 \end{bmatrix} + 1.3062 \begin{bmatrix} -.4 \\ -.75 \\ -1.034
\end{bmatrix} \notag \\
&\approx \begin{bmatrix} .478 \\ 1.020 \\ 1.649 \end{bmatrix}. \notag
\end{align}

\newpar
Comparing this value to that obtained using the SQPT,
given by
\newline $\ds \vecx_2=\vecx_1+\Delta \vecx \approx \begin{bmatrix} .6 \\ 1.25 \\ 1.966
\end{bmatrix}$, we see that the MSQPT yields a first iterate significantly closer to
the NLP's solution of $\ds \vecx_0= \begin{bmatrix} 0 \\ 1 \\ \sqrt{3} \end{bmatrix}$.

\newpar
The Maple worksheet, \textbf{Sequential Quadratic Programming Technique.mw},
is easily modified to incorporate use of the merit function. After the objective
and constraint functions, \mf{f}, \mf{g}, and \mf{h} have been defined, the penalty and merit functions are constructed using the following input:

\begin{maplelist}
\item \mf{P:=unapply(add(piecewise(g(x1,x2,x3)>=0,g(x1,x2,x3)[i]\^{}2,
\newline i=1..2) +h(x1,x2,x3)\^{}2,[x1,x2,x3])):}
\item \mf{M:=unapply(f(x1,x2,x3)+10*P(x1,x2,x3),[x1,x2,x3]):}
\end{maplelist}

\noi To account for the scaling of $\Delta \vecx$, we modify the worksheet after the computation
of \mf{w} as follows:

\begin{maplelist}
\item \mf{w:=evalf(MatrixInverse(B).<-Delf(X1),-h(X1),-($\sim$ d)>);}
\mo{w:=\begin{bmatrix} -.3994 \\ -.7500 \\ -1.0335 \\ .4268 \\ 1.5869 \end{bmatrix}}
\item \mf{phi:=t->M(X1[1]+t*w[1,1],X1[2]+t*w[2,1],X1[3]+t*w[3,1]):}
\item \mf{t0:=NewtonsMethod(phi(t),1,5,.01); \newline
\# NewtonsMethod procedure modified for a function of one 
\newline variable, using an initial value 1, maximum number of
\newline five iterations, and tolerance of .01.}
\mo{t0:=1.3062}
\item \mf{X2:=X1[1]+t0*w[1,1],X1[2]+t0*w[2,1],X1[3]+to*w[3,1];}
\mo{\begin{matrix} .478 & 1.020 & 1.649 \end{matrix}}
\end{maplelist}

\newpar
Subsequent iterations are handled in a similar manner.

\section*{Exercises Section \protect\ref{S:SequentialQuadraticProgramming}}

\addcontentsline{toc}{subsection}{Exercises}

\input{Exercises/Exercises_SequentialQuadraticProgramming}	
\endinput







