Efficient Implementation of Rate Constraints for Nonlinear Optimal Control

TL;DR

This paper introduces a novel on-mesh method for implementing rate constraints in nonlinear optimal control, improving solution properties and computational efficiency over traditional approaches.

Contribution

The paper presents a general on-mesh implementation technique for rate constraints applicable to all collocation methods, enhancing solution quality and computational speed.

Findings

Achieves over 30% computational speedup.

Prevents singular control arcs with the new method.

Applicable to all collocation methods.

Abstract

We propose a general approach to directly implement rate constraints on the discretization mesh for all collocation methods, for both state and input variables. Unlike conventional approaches that may lead to singular control arcs, the solution of this on-mesh implementation has better properties. Moreover, computational speedups of more than 30% can be achieved by exploiting the properties of the resulting linear constraint equations.

Tables

Table 1. TABLE I : Numerical differentiation schemes

Method	No. of Data	Numerical Differentiation
	Points ($r$)
Trapezoidal	2	$Z_i^{\prime }=\frac{Z_{i+1}-Z_i}{\Delta t\Delta {\tau }_i}$
($h$)	(equal spaced)
		$Z_i^{\prime }=\frac{-3Z_i+4Z_{i+1/2}-Z_{i+1}}{\Delta t\Delta {\tau }_i}$
Hermite	3	$Z_{i+1/2}^{\prime }=\frac{Z_{i+1}-Z_i}{\Delta t\Delta {\tau }_i}$
Simpson ($h$)	(equal spaced)
		$Z_{i+1}^{\prime }=\frac{Z_i-4Z_{i+1/2}+3Z_{i+1}}{\Delta t\Delta {\tau }_i}$
LGR	N+1	$Z_{1:N+1}^{\prime }=\frac{2}{\Delta t}{𝒜}_{LGR}{Z}_{1:N+1}$
($p$/$hp$)	(LGR bases)

Table 2. TABLE II : Contribution to the NLP dimensions with an N 𝑁 N -point mesh for different rate constraint implementations on input variables

	add-state	on-mesh
2 point,	$N$ linear	$2N$ linear
collocated	inequality constraints	inequality constraints
(Trapezoidal)	(defect constraints)	(pre-computed)
3 point,	$2N-1$ linear	$4N-4$ linear
collocated	inequality constraints	inequality constraints
(H-S)	(defect constraints)	(pre-computed)
$p$ point,	$p\left(N-1\right)$ linear	$2p\left(N-1\right)$ linear
collocated	inequality constraints	inequality constraints
(LGR)	(defect constraints)	(pre-computed)

Table 3. TABLE III : Contribution to the NLP dimensions with an N 𝑁 N -point mesh for different rate constraint implementations on state variables

	add-path constraint	on-mesh
2 point,	$2N$ nonlinear	$2N$ linear
collocated	inequality constraints	inequality constraints
(Trapezoidal)	(path constraints)	(pre-computed)
3 point,	$4N-2$ nonlinear	$6N-6$ linear
collocated	inequality constraints	inequality constraints
(H-S)	(path constraints)	(pre-computed)
$p$ point,	$2p\left(N-1\right)$ nonlinear	$2p\left(N-1\right)$ linear
collocated	inequality constraints	inequality constraints
(LGR)	(path constraints)	(pre-computed)

Table 4. TABLE IV: Mesh refinement criteria

	$d$	$h$	$v$	$\gamma$	$\alpha$	Path Constraint
	[ft]	[ft]	[ft/s]	[deg]	[deg]	[m]
${\eta }_{tol}$	1	0.5	0.1	0.5	0.5	-
${\epsilon }_{g_{tol}}$	1	0.5	0.1	0.5	0.5	$1\times {10}^{-5}$

Equations

$$\underset{x,u,p,t_0,t_f}{\min }\Phi (x(t_0),t_0,x(t_f),t_f,p)+{\int }_{t_0}^{t_f}L(x(t),u(t),t,p)dt$$

$$\displaystyle\dot{x}(t)=f(x(t),u(t),t,p),\$$

$$\displaystyle c(x(t),u(t),t,p)\leq 0,\$$

$$\displaystyle\phi(x(t_{0}),t_{0},x(t_{f}),t_{f},p)=0,\$$

$$x^{(k)}(\tau )\approx {\bar{x}}^{(k)}(\tau ):=\sum _{j=1}^{N^{(k)}}{\mathcal{X}}_j^{(k)}{\mathcal{B}}_j^{(k)}(\tau ),$$

$$\begin{gathered} \underset{X,U,p,t_0,t_f}{\min }\Phi (X_1^{(1)},t_0,X_f^{(K)},t_f,p) \\ +\sum _{k=1}^K\sum _{i=1}^{N^{(k)}}{w}_i^{(k)}L(X_i^{(k)},U_i^{(k)},{\tau }_i^{(k)},t_0,t_f,p) \end{gathered}$$

$$\sum _{j=1}^{N^{(k)}}{\mathcal{A}}_{ij}^{(k)}{X}_j^{(k)}+{\mathcal{D}}_i^{(k)}f(X_i^{(k)},U_i^{(k)},{\tau }_i^{(k)},t_0,t_f,p)=$$

$$c(X_i^{(k)},U_i^{(k)},{\tau }_i^{(k)},t_0,t_f,p)\le$$

$$\varphi (X_1^{(1)},t_0,X_f^{(K)},t_f,p)=$$

$${\dot{u}}_L\le \frac{du}{dt}(t)$$

$${\dot{x}}_L\le \frac{dx}{dt}(t)$$

$$\dot{u}(t)=\nu (t)\text{ with }{\dot{u}}_L\le \nu (t)\le {\dot{u}}_U.$$

$${\dot{x}}_L\le f(x(t),u(t),t,p)\le {\dot{x}}_U.$$

$$\underset{x_1,u,\nu }{\min }{\int }_0^{t_f}{g}_1(x_1(t))+g_2(x_1(t),u_1(t))dt$$

$$\dot{x_1}(t)=g_3(x_1(t),u_1(t))$$

$${\dot{u}}_1(t)={\nu }_1(t)$$

$${\dot{u}}_{1_L}\le {\nu }_1(t)\le {\dot{u}}_{1_U}$$

$$\begin{array}{rl}& H(x_1(t),u_1(t),{\lambda }_{x_1}(t),{\lambda }_{u_1}(t),{\nu }_1(t)):=g_1(x_1(t))\\ & +g_2(x_1(t),u_1(t))+{\lambda }_{x_1}(t)g_3(x_1(t),u_1(t))+{\lambda }_{u_1}(t){\nu }_1(t).\end{array}$$

$${\nu }_1^{\ast }(t)\in \arg \underset{v}{\min }H(x_1^{\ast }(t),u_1^{\ast }(t),{\lambda }_{x_1}^{\ast }(t),{\lambda }_{u_1}^{\ast }(t),v).$$

$${\nu }_1^{\ast }(t)=\{\begin{array}{ll}{\dot{u}}_{1_U}& \text{if }{\lambda }_{u_1}^{\ast }(t)<0\\ ?& \text{if }{\lambda }_{u_1}^{\ast }(t)=0\\ {\dot{u}}_{1_L}& \text{if }{\lambda }_{u_1}^{\ast }(t)>0\end{array}$$

$${\dot{z}}_L-Z_i^{\prime }$$

$$Z_i^{\prime }-{\dot{z}}_U$$

$$\underset{x_1,x_2,u}{\min }{\int }_0^5{x}_1^2(t)+x_2^2(t)dt$$

$$\dot{x_1}(t)=x_2(t),\dot{x_2}(t)=u(t)\in [-1,1]\forall t\in [0,5].$$

$$\dot{\alpha }(t)=\nu (t).$$

Full text

©2020 IEEE. Personal use of this material is permitted. Permission from IEEE must be obtained for all other uses, in any current or future media, including reprinting/republishing this material for advertising or promotional purposes, creating new collective works, for resale or redistribution to servers or lists, or reuse of any copyrighted component of this work in other works.

I Introduction

Optimization-based control strategies, such as model predictive control (MPC), can be seen in an increasing number of applications. For many engineering problems, constraints may need to be imposed on the rate of changes for the state and/or input variables, to account for physical actuation limitations (e.g. the maximum rotation rate of flight control surfaces on aircraft) or to fulfill certain ride comfort requirements (e.g. the maximum longitudinal and lateral accelerations experienced by passengers).

In optimal control, the underlying optimization problem can often be formulated and implemented in a number of different ways. Under a linear framework, many implementations are computationally comparable, thus straightforward approaches are often used. For example, rate constraints on input variables are generally implemented through additional dynamic equations [1, 2], and rate constraints on state variables are commonly addressed with additional path constraints [3]. However, under a nonlinear framework, this way of implementing input rate constraints are known to result in numerical difficulties and introducing fluctuations and ringing phenomena in the solution due to singular control [4]. To improve the solution quality, additional regularization terms may be added to the optimal control problem (OCP) formulation [5]; however, this practice often leads to computational challenges by needing to solve the problem repetitively with appropriate weightings.

In our previous work, we proposed implementing rate constraints directly on the discretization mesh with linear constraints [6]. In this paper, we present a more in-depth analysis of this method. We demonstrate that the proposed method will not introduce singular arcs to the problem, resulting in solutions of higher accuracy than the conventional approach. The computational comparisons with the conventional implementation as well as the regularization approach are also notably expanded with further insights.

Sections II–III aim at providing a brief introduction on solving OCPs with direct collocation methods. Following this, different approaches for implementing rate constraints in the OCP are introduced and analysed in Section IV. This is followed by two classical examples of different complexity in Section V, where the pros and cons of each implementation are demonstrated. In Section VI, we provide concluding remarks and some guidelines for implementation.

II Optimal Control Problem

Generally speaking, optimization-based control requires the solution of optimal control problems with the objective functional expressed in the general Bolza form:

[TABLE]

with $x:\mathbb{R}\rightarrow\mathbb{R}^{n}$ is the state trajectory of the system, $u:\mathbb{R}\rightarrow\mathbb{R}^{m}$ is the control input trajectory, $p\in\mathbb{R}^{s}$ are static parameters, $t_{0}\in\mathbb{R}$ and $t_{f}\in\mathbb{R}$ are the initial and terminal time. $\Phi$ is the Mayer cost functional ($\Phi$: $\mathbb{R}^{n}\times\mathbb{R}\times\mathbb{R}^{n}\times\mathbb{R}\times\mathbb{R}^{s}\to\mathbb{R}$), $L$ is the Lagrange cost functional ($L:\mathbb{R}^{n}\times\mathbb{R}^{m}\times\mathbb{R}\times\mathbb{R}^{s}\to\mathbb{R}$), $f$ is the dynamic constraint ($f:\mathbb{R}^{n}\times\mathbb{R}^{m}\times\mathbb{R}\times\mathbb{R}^{s}\to\mathbb{R}^{n}$), $c$ is the path constraint ($c:\mathbb{R}^{n}\times\mathbb{R}^{m}\times\mathbb{R}\times\mathbb{R}^{s}\to\mathbb{R}^{n_{g}}$) and $\phi$ is the boundary condition ($\phi:\mathbb{R}^{n}\times\mathbb{R}\times\mathbb{R}^{n}\times\mathbb{R}\times\mathbb{R}^{s}\to\mathbb{R}^{n_{q}}$).

III Direct collocation method

Most optimal control problems need to be solved with numerical discretization schemes in practice. With a direct method, the OCP is first discretized through a transcription process, after which the resulting nonlinear programming (NLP) problem is numerically solved. Thanks to its simplicity in implementation, direct methods have become the de facto standard for solving practical optimal control problems [7].

One approach in direct methods is to solve the dynamics equations, the path constraints and the boundary conditions altogether on a discretization mesh. This is often referred to as direct collocation methods. Moreover, it can be further categorized into fixed-order $h$ methods (e.g. Euler, Trapezoidal, Hermite-Simpson (H-S) and the Runge-Kutta (RK) family) [4], and variable higher-order $p$/$hp$ methods [8, 9]. Here, we aim to provide a high level overview, which is valid for both $h$ and $p$/$hp$ methods.

With a mesh of size $N=\sum_{1}^{K}N^{(k)}$, the states can be approximated as

[TABLE]

with mesh interval $k$ $\in$ $\{1$, $\ldots$, $K\}$, $N^{(k)}$ denoting the number of collocation points for interval $k$, and $\mathcal{B}_{j}^{(k)}(\cdot)$ are basis functions. For classical $h$ methods, $\tau\in\mathbb{R}^{N}$ takes on values in the interval $[0,1]$ representing $[t_{0},t_{f}]$, and $\mathcal{B}_{j}^{(k)}(\cdot)$ are chosen to be elementary B-splines of various orders. For $p$/$hp$ methods, $\mathcal{B}_{j}^{(k)}(\cdot)$ are Lagrange interpolating polynomials over the normalized time interval $\tau$ $\in$ $[-1,1]$. We use $X_{j}^{(k)}$ and $U_{j}^{(k)}$ to represent the approximated states and inputs at collocation points, e.g.  $X_{j}^{(k)}=\bar{x}^{(k)}(\tau_{j}^{(k)})\in\mathbb{R}^{n}$, where $\tau_{j}^{(k)}$ is the $j^{\text{th}}$ collocation point in mesh interval $k$.

Consequently, the optimal control problem (1) can be approximated by

[TABLE]

where $w_{j}^{(k)}$ are the quadrature weights for the respective discretization method chosen, $\mathcal{A}$ is the numerical differentiation matrix with $\mathcal{A}_{ij}$ the element $(i,j)$ of the matrix, and $\mathcal{D}$ a constant matrix. The discretized problem can then be solved with off-the-shelf NLP solvers, such as interior point solver IPOPT [10].

The NLP solver outputs a discretized solution $\mathcal{Z}\coloneqq(X,U,p,\tau,t_{0},t_{f})$ as sampled data points, however it does not necessarily indicate that the solution is a discrete-time control sequence. In-between the sampled points, interpolating splines can be used to construct an approximation of the continuous-time optimal trajectory $t\mapsto\tilde{\mathcal{Z}}(t)\coloneqq(\tilde{x}(t),\tilde{u}(t),t,p)$ in accordance to the discretization method employed. The quality of the interpolated solution needs to be assured through error analysis. If necessary, modifications must be made accordingly to the discretization mesh, until the solutions obtained with the new mesh fulfills all predefined error tolerance levels (e.g. the absolute local error $\eta_{tol}$ and the absolute local constraint violation $\varepsilon_{c_{tol}}$). This process is known as mesh refinement.

IV Implementations of Rate Constraints

In many problems, constraints of the form

[TABLE]

may need to be implemented to restrict the rate of change for the state and/or input variables.

IV-A Conventional Implementation

For input variables, a common approach is to introduce $u$ as an additional state variable, and $\nu$ as the new input with a simple bound through the dynamic equation

[TABLE]

For rate constraints on the state variable $x$, additional path constraints are needed:

[TABLE]

For simplicity, we refer to (4) as the add-state implementation, and (5) as the add-path constraint implementation.

Unfortunately these conventional implementations exhibit many shortcomings. These are mainly:

The number of state variables and constraint equations are increased, resulting in a larger NLP. In addition, the index of the DAE (differential-algebraic equations) of the transcribed problem may also increase, leading to a problem that is often more difficult to solve numerically. 2. 2.

When (4) is used, singular arcs may occur and affect the solution quality. This can occur if the original control input $u$ appears nonlinearly in the Lagrange cost or other system dynamics and the new control $\nu$ appears linearly instead.

Detailed justification for the first point can be found in [4]. Here, we focus the discussions on the issue of singular arcs, which also falls with a much larger concern regarding the quality of the solution.

Generally speaking, a singular arc is an interval in the OCP solution where the optimality conditions yield no information about the optimal control function. Precise mathematical definitions and an analysis of singular arcs can be found in [11] and [12]. For ease of demonstration, consider the following OCP, which is simplified but still sufficiently general, with $x_{1}\in\mathbb{R}$ the state variable, and $u_{1}\in\mathbb{R}$ the control input. As per the conventional approach, the rate constraint on the original control input is implemented with a new control input $\nu_{1}\in\mathbb{R}$:

[TABLE]

Here we follow the same hypotheses as in [13]: $g_{1}$, $g_{2}$ and $g_{3}$ are continuous, continuously differentiable for all $u_{1}\in\mathcal{U}$, and Lipschitz in $x_{1}$. Also, the admissible control set $\mathcal{U}$ is assumed to be a bounded set in some Euclidean space.

Proposition 1**.**

If the OCP (6) has a linear objective and dynamics with respect to the original control input $u_{1}$, i.e.  if $g_{2}$ and $g_{3}$ are both functions that only have strictly linear input (i.e. in the form of $g_{\iota}(x_{1}(t),u_{1}(t))=\tilde{g}_{\iota}(x_{1}(t))+c_{st}u_{1}(t)$, with $c_{st}$ a constant), the resulting optimal control $\nu_{1}^{\ast}$ will not contain a singular arc. However, if $u_{1}$ appears nonlinearly in the objective and/or dynamics, i.e. if $g_{2}$ and/or $g_{3}$ are arbitrary nonlinear functions, there exists problems where singular arcs will occur for some intervals of the solution.

Proof.

First, we formulate the Hamiltonian of the system, with $\lambda(t)$ the costate of the dynamics

[TABLE]

From Pontryagin’s minimum principle, we know that if the state and costate are optimal, the optimal control $\nu_{1}^{\ast}$ minimizes the Hamiltonian, i.e.

[TABLE]

Substituting the Hamiltonian (7) in (8) yields the optimal control

[TABLE]

We first show that this implementation yields a singular arc free solutions with linear OCPs. From the first order necessary conditions for optimality we have $\dot{\lambda}_{u_{1}}(t)=-\frac{\partial H}{\partial u_{1}}$. Thus, if both $g_{2}$ and $g_{3}$ only contain input terms in the form of $c_{st}u_{1}(t)$, $\dot{\lambda}_{u_{1}}(t)$ will then be constant and $\lambda_{u_{1}}(t)$ will be a linear straight line. Relating this to (9) we can see that the optimal control will exhibit bang-bang behaviour with at most one switch depending on the crossing of $\lambda_{u_{1}}(t)$ with the $x$ axis. Therefore, the solution is free of singular arcs if the objective and dynamics with respect to the original control input $u_{1}$ are all linear. This is the reason why in linear optimal control problems, the conventional implementation can be used without issues.

To show that this implementation will suffer from singular arc problems when nonlinear OCPs are considered, we assume that $g_{2}$ and/or $g_{3}$ are now arbitrary nonlinear functions. Thus, $\lambda_{u_{1}}(t)$ can be a function of any shape and the optimal control will not be uniquely defined on intervals where $\lambda_{u_{1}}(t)=0$, a.k.a. the singular arc. The problem in Section V-A is an example where such an issue arises. ∎

For a direct collocation method to yield the correct solution for singular control problems, one might have to use a multi-phase formulation and additionally impose the singular arc condition specifically on the phases with singular control. For example, if one takes the same approach (as in the proof) for the example problem in Section V-A, the condition for singular control to occur is when $\lambda_{x_{2}}=0$. Repetitively taking time derivatives of this equation would yield the singular arc condition $u(t)=x_{1}(t)$. We again note that this approach requires analytical derivations and would become increasingly challenging for complex real-world problems.

An ad-hoc method sometimes used in practice for dealing with the singular arc is to augment the original objective with an additional regularization term (e.g. in [5]), often in the form of $\rho||\nu||_{\mathcal{L}_{1}}$ or $\rho||\nu||_{\mathcal{L}_{2}}^{2}$. With relatively large values of the penalty weight $\rho$, the fluctuations on the singular arc can be suppressed, but at the cost of obtaining sub-optimal trajectories. To get closer to the optimal from this point, the problem may need to be repetitively solved with the penalty weight gradually reduced.

IV-B On-mesh Implementation

To mitigate the above-mentioned shortcomings, a method is proposed to directly impose algebraic rate constraints for input variables on the discretization grid. Based on previous work [4], we generalize this on-mesh approach for all collocation methods ($h$, $p$ and $hp$ type), as well as for state variables.

Since the treatment for state variables $x$ and input variables $u$ are similar, for simplicity we will use $z$ to represent the variable on which the rate constraints are imposed. If $Z_{i}$ represents the discretized version of $z$ at time instance $i$, then the numerical differentiation of $z$ at that grid point ($Z_{i}^{\prime}$) can be calculated using $s$-point finite difference approximations, with $s$ the number of data points in the interval (including endpoints). See Table I for the formulations of some of the most commonly used discretization methods, with $\Delta\tau_{i}=\tau_{i+1}-\tau_{i}$, $\Delta t=t_{f}-t_{0}$, and $\mathcal{A}_{LGR}$ is the Legendre-Gauss-Radau (LGR) differentiation matrix. Details on the determination of the numerical differentiation equations are available in [14].

Note that for $p$/$hp$ methods, the numerical differentiation for all grid points on the polynomial ($i=1,\ldots,N^{(k)}$) are obtained altogether. It is also worth mentioning that if Legendre-Gauss-Radau (LGR) collocation is used, the end-point value for the control ($U_{N+1}^{(K)}$) may need to be approximated.

It is then straightforward to implement the rate constraints as linear constraints

[TABLE]

for all possible values of $i$. This approach will be referred to as the on-mesh implementation.

The on-mesh implementation of rate constraints has several benefits in comparison to the conventional add-state and add-path constraint approaches. Firstly, we compare the solution quality in terms of singular arcs. A challenge arises here, since the singular arc problem is commonly analyzed with the original OCP (1), but (10) is a discretized formulation that does not have the continuous form. Therefore, we do not yet have a mathematically rigorous proof that the on-mesh implementation will be singular-arc free using Pontryagin’s minimum principle. However, one would observe that, without introducing a dynamic constraint in the form of $\dot{u}(t)=\nu(t)$, the singular control situation as described in Proposition 1 will not occur, at least not in the same way. In addition, in our computational experience, we had not encountered a single case where the on-mesh implementation causes a singular-arc free problem to become singular. In contrast, the on-mesh rate constraint method was able to convert many well-known singular control problems to be singular arc free ones, with one example demonstrated in Section V-A.

Another major advantage in comparison to the conventional implementation is regarding the computational cost. For systems with nonlinear dynamics, rate constraints on state variables implemented through (5) will be nonlinear path constraints relating different state variables at the same time instance. Thus their Jacobian and Hessian contributions can make the solution of the OCP computationally demanding. In contrast, on-mesh implementation of rate constraints with (10) are linear constraints, with no contribution to the Hessian.

In addition, note that the rate constraints (10) only depend on the numerical differentiation schemes. Thus, once a discretization scheme for the OCP has been chosen, and the corresponding discretization mesh has been determined, the Jacobian contributions of the rate constraint equations can be pre-computed during the transcription process. Therefore, although the NLP dimension increases more rapidly with the on-mesh implementation as shown in Table II and III, the computational complexity for obtaining the derivative information with respect to the rate constraint equations is actually lower than the conventional approach. Altogether, the computational advantages can be rather significant, as demonstrated with the example problem.

A remark is appropriate when comparing the on-mesh implementation in Table III to Table II: For Hermite-Simpson discretization, specifically, the increase in the size of the NLP for implementation on input variables is less than that on state variables. This is because, when the control $u$ is discretized as a quadratic function of time, the rate of change w.r.t. time ($\dot{u}$) is linear, thus extreme values only occur at the end-points of each interval ($U_{i}$ and $U_{i+1}$). In this special case only, the rate constraints relating to the middle points ($U^{\prime}_{i+1/2}$) can be neglected.

V Example Problems

The problem of singular arcs is often demonstrated with toy problems in the literature (e.g. the first example), as they are much more illustrative and free from influence of other factors. However, this common practice often results in it being neglected by engineers working on complex problems. To show that it really matters, a second real-world example is also presented here to demonstrate the acclaimed benefits of the on-mesh rate constraint implementation in terms of solution quality and computational efficiency.

V-A Second order singular regulator

First, we consider a simple regulator problem originally presented in [15]. This can be considered as regulation control of a double integrator system, with a constraint on the acceleration.

[TABLE]

In this original OCP formulation, the optimal control is in the form of bang-singular. Figure 1 shows that a numerical implementation of this OCP using direct collocation would yield fluctuations and ringing phenomena for solutions at collocation points, as well as for the trajectories in-between.

However, upon noticing that the second differential equation is equivalent to the add-state implementation of a rate constraint $-1\leq\dot{x_{2}}(t)\leq 1$, we can remove that differential equation from the OCP and use the on-mesh rate constraint method instead. As illustrated in the figure, this approach successfully yields a stable and accurate solution in correspondence to the reference (analytical) optimal input trajectory. In other words, the proposed on-mesh rate constraint implementation has successfully converted the classical second order singular regulator problem into a singular-arc-free formulation.

V-B Aircraft go-around in the Presence of Windshear

Based on the developments by [16, 17, 18], a problem is presented in [4] where the aircraft needs to stay as high above the ground as possible after encountering a severe windshear during landing.

Details about the system dynamics, variable simple bounds, boundary conditions, aerodynamic modelling, as well as static parameter values, are available in the references above. The problem has state variables being the horizontal distance $d$, the altitude $h$, the true airspeed $V$, the flight path angle $\gamma$, and the input variable being the angle of attack $\alpha$ . We also emphasise that the problem requires the implementation of a rate constraint $\lvert\dot{\alpha}(t)\rvert\leq 3$ deg/s.

To avoid discontinuities and to assist convergence, a static optimization parameter $h_{min}$ is introduced to represent the minimum altitude. The objective can then be expressed as $\Phi(x(t_{0}),t_{0},x(t_{f}),t_{f},p):=-h_{min}$ together with a new path constraint $h(t)\geq h_{min}$.

Figure 2 illustrates the solution to the problem using Hermite-Simpson discretization.

All figures presented in this paper are the outcome of a mesh refinement scheme that minimizes errors to the tolerance as specified in Table IV.

V-B1 Implementation of Rate Constraints for Input Variables

Constraint $\lvert\dot{\alpha}(t)\rvert\leq 3$ applies directly on the rate of change for the control input $\alpha$. Using the conventional approach, $\alpha$ can be treated as an additional state variable, and $\nu$ introduced as the new control input with the dynamics

[TABLE]

Thus the rate constraints for $\alpha$ can be implemented as simple bounds on $\nu$: $-3\leq\nu(t)\leq 3$ [deg/s].

As mentioned earlier, due to the fact that the original control input $\alpha$ appears nonlinearly in the system, whereas the new input $\nu$ appears linearly, singular arc behaviour can occur, which is shown in Figure 3(a), with large fluctuations in the solution. In contrast, when the rate constraints are directly implemented on the discretization mesh instead (Figure 3(b)), the optimal control input trajectory can be obtained with little ambiguity.

Comparing the solutions from the two implementations, it is interesting to observe that, although the integrated values (i.e. angle of attack) along the singular arc solution at the collocation points are generally the same, the interpolated trajectory from the add-state method is actually distorted by the fluctuations of its rate values.

With the LGR orthogonal collocation method, improvements are relatively minor. Because the end-point value for the control input is only approximated (not a collocation point), the errors have distortion effects on all previous points of the polynomial (Figure 4(b)). On the other hand, thanks to this extra level of continuity imposed by higher order polynomials, the problem of singular arc behaviour is far less pronounced with the conventional add-states implementation (Figure 4(a)), when compared to $h$ type discretization methods.

The results regarding computation times presented in Figure 5 and 7 were all obtained on an Intel Core i7-4770 computer, running 64-bit Windows 10 with Matlab 2017a. The OCPs were transcribed into NLP problems using the optimal control software ICLOCS2 [19] and solved with the NLP solver IPOPT compiled with the sparse linear solver MA57 [20]. The computation times are the averages of 10 independent runs, all starting with a very rough initial guess obtained using linear interpolation of initial and terminal conditions.

From Figure 5 it can be seen that, for the computation time per iteration, the on-mesh implementation saw a slight advantage in comparison to the conventional approach. This is because the on-mesh implementation explicitly exploits the fact that the linear rate constraints have no contribution to the Hessian, and the contributions to the Jacobian are constants and can be pre-computed. The scale of the benefit also grows with the size of the mesh, from about 5% for a coarse mesh to around 10% for the dense mesh.

Figure 6 presents the computation performance of the problem when regularized with an additional $\rho||\nu||_{\mathcal{L}_{2}}^{2}$ term. It can be seen that a relative large penalty weight is required to suppress the singular arc fluctuations, but with a larger $\rho$ the result diverges quickly from the optimal. Also note that for a single solve with regularization, the norm of angle of attack rate ($||\dot{\alpha}^{\ast}-\dot{\alpha}||_{\mathcal{L}_{2}}$) never reaches the accuracy level obtained by the on-mesh implementation with the same discretization mesh. Thus, to obtain a good solution, $\rho$ needs to be gradually reduced, making the process complicated and computationally inefficient — it is also difficult to guarantee solution quality.

V-B2 Implementation of Rate Constraints for State Variables

We will additionally impose a rate constraint for the velocity state as $-5\leq\dot{V}(t)\leq 5$ [ft/s2]. With this new formulation, the minimum altitude achievable is slightly lower.

From Figure 7, it is obvious that the two methods are not computationally comparable. Due to the reasons explained in the end of Section IV-B, although the increase in NLP dimension is higher for the on-mesh implementation, the resulting (larger) NLP problems with linear rate constraints are actually much easier to solve. Consequently, regardless of the discretization method, the computation time per iteration recorded for the on-mesh implementations are all significantly (more than 30%) lower than the conventional method.

VI Conclusions

Through both the mathematical analysis and a computation example, we demonstrated that mathematically equivalent formulations for rate constraints on state and input variables may not have the same solution quality and computational complexity in numerical implementations. For all collocation methods tested, and for both state and input variables, the proposed approach to directly implement rate constraints on the discretization mesh appears to be an attractive alternative for nonlinear optimization based control problems.

In contrast to conventional approaches, the proposed method is numerically verified to not introduce singular control arcs, a phenomena which may lead to severe distortions and fluctuations in the optimal control trajectories. Additionally, this on-mesh implementation replaces the additional dynamic equations and nonlinear path constraints in conventional implementations with linear rate equations. Thus, there is no contribution to the Hessian and the contribution to the Jacobian can be pre-computed, enabling faster iterations. Based on observations, the scale of reduction in computation time seems to grow quite quickly with the increase in mesh size (number of collocation points), making the method especially suitable for the solution of large scale problems. A possible downside of the proposed scheme is that the method cannot be directly implemented in most of the existing OCP packages through their current interfaces. However, we believe that this implementation can be easily added by the authors of the software, and ensure that all computational benefits are fully realized.