Numerical Verification

In this section it is our goal to first demonstrate that problem (P) can be solved to good accuracy using a finite difference method. Next, as was done in [7,8] we will, through an eigenvalue computation, verify that the computed solution satisfies the SSC and is thus a local minimizer. Then, in order to check the properties of the specific example shown analytically above, we will compute an additional eigenvalue. We start by presenting the finite-dimensional analogue of (P) and then outline how the algebraic problems are solved. We define the following discretization of problem (P).

$$\begin{eqnarray*} &&\mbox{minimize}\quad f_h(y_h,u_h)=\\ &&\frac{dxdt}2\sum^{m}... ...igg(\sum^{m}_{i=0}\gamma_{i}(a_y(t_i)y_{n,i}+a_u(t_i)u_i) \bigg) \end{eqnarray*}$$

     subject to(P$_h$)

$$\begin{eqnarray*} &&\quad\frac{y_{j,i+1}-y_{j,i}}{dt}=\frac12(y_{j-1,i}-2y_{j,i}... ...xdt\sum^{m}_{i=0}\sum^{n}_{j=0}\beta_{j}\gamma_{i}y_{j,i} \le 0. \end{eqnarray*}$$

Here $x_j=jdx$, $dx=\pi /n$, $t_i=idt$, $dt=T/m,\beta_0 =\beta_n=\frac12, \beta_j =1$ otherwise; analogously for $\gamma$.

The discrete control problem (P$_h$) is essentially a generic nonlinear optimization problem of the form

$$\min F^h(z)\quad\mbox{subject to}\quad G^h(z)=0,\quad H^h(z)\le0$$ (4.1)

where $z$ comprises the discretized control and state variables. $G^h(z)$ symbolizes the state equation and boundary conditions while $H^h(z)$ stands for both pointwise control bounds and the integral state constraint, the only constraints of inequality type prescribed above. We state the well-known SSC for (4.1), assuming $z\in\mathbf{R}^{N_h}$, $G^h:\mathbf{R}^{N_h}\to\mathbf{R}^{M_h}$, $M_h<N_h$. Let $z^*$ be an admissible point satisfying the first-order necessary optimality conditions with associated Lagrange multipliers $\mu^*$ and $\lambda^*$. Let further

$$N(z^*)=(\nabla G^h(z^*),\nabla H_a(z^*))$$

be a column-regular ${N_h}\times(M_h+P_h)$ matrix where $M_h+P_h<N_h$ and $\nabla H_a(z^*)$ denotes the gradients of the $P_h$ active inequality constraints with positive Lagrange multipliers. For (4.2) we have $N_h=(m+1)(n+2)$ and $M_h=(m+1)(n+1)$ resulting in $m+1$ degrees of freedom which are further reduced by one through the active integral nonnegativity constraint on $y$ and by any active bounds on $u$. Let finally $N=QR$ be a QR decomposition and $Q=(Q_1,Q_2)$ a splitting into the first $M_h+P_h$ and the remaining columns. The point $z^*$ is a strict local minimizer if a $\gamma>0$ exists such that, see, for example, [13]

$$\lambda_{\min}(L_2(z^*))=\gamma>0.$$ (4.2)

Here $L_2(z^*)$ is the projected Hessian of the Lagrangian

$$L_2(z^*)=Q^T_2(\nabla^2F^h(z^*)-\mu^{*T}\nabla^2G^h(z^*))Q_2.$$

No Hessian of $H^h$ appears on the right due to its linearity. To clarify the relationship between the way we proceed here and which is standard in optimization and the analysis of the previous section we add the following explanations. In order to verify the SSC in the discrete case we have to determine the smallest eigenvalue of the Hessian on the tangent space of the active constraints. We do this by explicitly computing the orthogonal projection matrix $Q_2$ onto the tangent space and forming $L_2(z^*)$. Due to the verified regularity of the computed solution or the nondegeneracy of the active constraints the tangent space is equal to the nullspace of the active constraint gradients and thus the smallest eigenvalue of $L_2(z^*)$ corresponds to the minimal value of its quadratic form on the space of all $(y,u)$ satisfying the linearized equation as well as having vanishing $u$ components corresponding to indices $i$ for which the solution is at the bound which coincidentally also is zero. These components do include the ones corresponding to the interval $[0,1/4]$. Next, we will detail how condition (4.2) will be checked.

As was already done in [7,8] the control problems are written in the form of AMPL [6] scripts. This way, a number of nonlinear optimization codes can be utilized for their solution. It had been an observation in our previous work that from the then available codes only LOQO [14] was able to solve all the problems effectively and for sufficiently fine discretizations. This has changed. As recent comparisons [9] have shown, the trust region interior point method KNITRO [2] which became available only recently, may outperform LOQO on such problems. It was used for the computations reported below. The following procedure is independent of the solver used.

After computing a solution an AMPL stub (or $*.nl$) file is written as well as a file with the computed Lagrange multipliers. This allows to check the SSC (4.2) with the help of a Fortran, alternatively, a C or Matlab, program. This program reads the files and verifies first the necessary first-order optimality conditions, the column regularity of $N(z^*)$ and the strict complementarity. For this, it utilizes routines provided by AMPL which permit evaluation of the objective and constraint gradients. Next, the the QR decomposition of $N(z^*)$ is computed by one of the methods exploiting sparsity. We have utilized the algorithm described in [12]. AMPL also provides a routine to multiply the Hessian of the Lagrangian times a vector. This is called with the columns of $Q_2$ and thus $L_2(z^*)$ can be formed. Its eigenvalues are computed with LAPACK routine DSYEV and the smallest eigenvalue $\gamma=\gamma_h$ is determined. The use of this eigenvalue routine is possible since the order of the matrices corresponding to the "free" control variables is moderate. In case of distributed control problems when this number may be on the order of the state variables, a sparse solver, preferably just for finding the minimal eigenvalue, will have to be used.

With the procedure described above the SSC for problem (P$_h$) can be checked for constant or variable $\alpha_0$. In the nonconstant case and for $\alpha_0$ below the bound given above an additional eigenvalue problem is solved. Let $Q$ in the previous section be split into $Q=(Q_1,Q_2)$ where now $Q_1$ corresponds to the equality constraints only and thus has $M_h$ columns. Then, in analogy to (4.2) we define $L_2(z^*)$ and call its smallest (leftmost on the real line) eigenvalue $\delta_h$. We will have to obtain a negative $\delta_h$ for sufficiently negative $\alpha_0$. As described above, this $Q_2$ projects onto the nullspace of the equality constraints only and thus $L_2(z^*)$ is the projection of the Hessian onto the larger subspace of pairs $(y,u)$ satisfying the linearized equation only and for which $u$ may be nonzero everywhere on $[0,T]$.

Problem (4.2) was solved as described above for two discretizations which were chosen to be about equidistant in both coordinates. In Table 1 are the errors of both state and control listed in two different norms.

Table 1: Solution errors for problem 4.2

$1/dx$	$1/dt$	$\quad \Vert u-\bar{u}\Vert _\infty$	$\Vert u-\bar{u}\Vert _2$	$\Vert y-\bar{y}\Vert _\infty$	$\Vert y-\bar{y}\Vert _2$
127	41	3.739e-2	4.590e-5	6.096e-3	2.730e-7
192	61	1.331e-2	1.152e-5	1.770e-3	1.850e-7

In Table 2 we list the eigenvalues for both discretizations and various values of $\alpha_0$. As can be seen the sign change for $\delta_h$ occurs in both cases between $-10.5$ and $-11$ while the estimate (3.5) above yielded a bound of $-10.845$. The computed state is shown in the same figure.

Table 2: Minimal eigenvalues for problem 4.2

$1/dt$	$\alpha_0$	$\gamma_h$	$\delta_h$
41	-11	5.804e-5	-3.846e-5
	-10.5	5.804e-5	4.320e-5
	-10	5.804e-5	4.878e-5
61	-30	3.589e-5	-1.784e-3
	-15	3.589e-5	-3.869e-4
	-11	3.589e-5	-2.121e-5
	-10.5	3.589e-5	2.337e-5
	-10	3.589e-5	3.279e-5
	1	3.589e-5	3.279e-5

Figure 1: Example 4.2, Optimal control and state