Efficient and More Accurate Representation of Solution Trajectories in Numerical Optimal Control

TL;DR

This paper introduces a novel trajectory representation method for optimal control that minimizes residual errors across the entire trajectory, resulting in higher accuracy and computational efficiency compared to traditional interpolation schemes.

Contribution

It proposes a new residual minimization approach that improves solution accuracy and reduces mesh resolution requirements in numerical optimal control.

Findings

Higher accuracy in trajectory solutions with fewer collocation points

Reduced computational complexity due to smaller problem dimensions

Enhanced robustness of solutions in optimal control applications

Abstract

We show via examples that, when solving optimal control problems, representing the optimal state and input trajectory directly using interpolation schemes may not be the best choice. Due to the lack of considerations for solution trajectories in-between collocation points, large errors may occur, posing risks if this solution is to be applied. A novel solution representation method is proposed, capable of yielding a solution of much higher accuracy for the same discretization mesh. This is achieved by minimizing the integral of the residual error for the overall trajectory, instead of forcing the errors to be zero only at collocation points. In this way, the requirement for mesh resolution can be significantly reduced, leaving the problem dimensions relatively small. This particular formulation also avoids some of the drawbacks found in the earlier work of integrated residual minimization, leading to more efficient computations.

Tables

Table 1. TABLE I: Typical numerical schemes

Method	Numerical Integration Scheme
(Order)
Euler (1)	$X_2^{\left(k\right)}=X_1^{\left(k\right)}+h_k{F}_1^{\left(k\right)}$
Trapezoidal (2)	$X_2^{\left(k\right)}=X_1^{\left(k\right)}+\frac{h_k}{2}\left(F_1^{\left(k\right)}+F_2^{\left(k\right)}\right)$
Hermite	$X_2^{\left(k\right)}=\frac{1}{2}\left(X_2^{\left(k\right)}+X_1^{\left(k\right)}\right)+\frac{h_k}{8}\left(F_1^{\left(k\right)}-F_3^{\left(k\right)}\right)$
Simpson (3)	$X_3^{\left(k\right)}=X_1^{\left(k\right)}+\frac{h_k}{6}\left(F_1^{\left(k\right)}+4F_2^{\left(k\right)}+F_3^{\left(k\right)}\right)$
LGR ($N^{\left(k\right)}$)	${ℐ}^{\left(k\right)}={\left[{𝒜}_{2:N+1}^{\left(k\right)}\right]}^{-1}$
	$X_{2:N+1}^{\left(k\right)}=X_1+\frac{\Delta t}{2}{ℐ}^{\left(k\right)}{F}_{1:N}^{\left(k\right)}$

Table 2. TABLE II: Continuity of the reconstructed solution (a.m.: at most; p.w.: piecewise)

Method	Dynamics ($\dot{\tilde{x}}$)	States ($\tilde{x}$)	Inputs ($\tilde{u}$)
Euler	p.w. constant	a.m.p.w. linear	same as dynamics
Trape.	a.m.p.w. linear	a.m.p.w. quad.
H-S	a.m.p.w. quad.	a.m.p.w. cubic
LGR	a.m.  order $N^{\left(k\right)}$	a.m. order $N^{\left(k\right)}$+1

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} J_c:=\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)=$$

$$F_i^{(k)}:=f(X_i^{(k)},U_i^{(k)},{\tau }_i^{(k)},t_0,t_f,p).$$

$$\tilde{x}^{(k)}(\mathcal{Z},t)=X_{1}^{(k)}+F_{1}^{(k)}(t-t_{1}^{(k)})\\ +\frac{1}{2}\bigg{(}-3F_{1}^{(k)}+4F_{2}^{(k)}-F_{3}^{(k)}\bigg{)}\frac{(t-t_{1}^{(k)})^{2}}{h_{k}}\\ +\frac{2}{3}\bigg{(}F_{1}^{(k)}-2F_{2}^{(k)}+F_{3}^{(k)}\bigg{)}\frac{(t-t_{1}^{(k)})^{3}}{h_{k}^{2}},$$

$$\dot{\tilde{x}}^{(k)}(\mathcal{Z},t)=F_{1}^{(k)}+\bigg{(}-3F_{1}^{(k)}+4F_{2}^{(k)}-F_{3}^{(k)}\bigg{)}\frac{t-t_{1}^{(k)}}{h_{k}}\\ +\bigg{(}2F_{1}^{(k)}-4F_{2}^{(k)}+2F_{3}^{(k)}\bigg{)}\bigg{(}\frac{t-t_{1}^{(k)}}{h_{k}}\bigg{)}^{2},$$

$$\begin{gathered} {\tilde{u}}^{(k)}(\mathcal{Z},t)=\frac{2}{h_k^2}(t-\frac{1}{2}{t}_1^{(k)}-\frac{1}{2}{t}_3^{(k)})(t-t_3^{(k)})U_1^{(k)} \\ -\frac{4}{h_k^2}(t-t_1^{(k)})(t-t_3^{(k)})U_2^{(k)} \\ +\frac{2}{h_k^2}(t-t_1^{(k)})(t-\frac{1}{2}{t}_1^{(k)}-\frac{1}{2}{t}_3^{(k)})U_3^{(k)}. \end{gathered}$$

$${\epsilon }_r(t):=\dot{\tilde{x}}(\mathcal{Z},t)-f(\tilde{x}(\mathcal{Z},t),\tilde{u}(\mathcal{Z},t),t,p).$$

$${\eta }_j:={\int }_{t_j}^{t_{j+1}}\Vert {\epsilon }_r(s){\Vert }_{2}ds,$$

$${\sigma }_{j,q}:={\int }_{t_j}^{t_{j+1}}|{\epsilon }_{r_q}(s)|ds,\text{ for }q=1,\dots ,n,$$

$${\epsilon }_{c_{\zeta }}(t):=\{\begin{array}{ll}0& \text{if }{c}_{\zeta }(\tilde{z}(t))\le 0\\ c_{\zeta }(\tilde{z}(t))& \text{if }{c}_{\zeta }(\tilde{z}(t))>0\end{array},\text{ for }\zeta =1,\dots ,n_g.$$

$$\underset{\hat{x},\hat{u},t,p}{\min }{\int }_{t_0}^{t_f}r(\hat{x}(t),\hat{u}(t),t,p)dt$$

$$r(\hat{x}(t),\hat{u}(t),t,p):=\Vert \dot{\hat{x}}(t)-f(\hat{x}(t),\hat{u}(t),t,p){\Vert }_2^2.$$

$$F_2^{(k)}=$$

$$F_3^{(k)}=$$

$$\underset{X,U,p,t_0,t_f}{\min }\sum _{k=1}^{K}R(X^{(k)},U^{(k)},{\tau }^{(k)},{\tau }_q^{(k)},t_0,t_f,p)$$

$$\begin{array}{rl}\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)& \\ +\Phi (X_1^{(1)},t_0,X_f^{(K)},t_f,p)\le & J_c\end{array}$$

$$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)=$$

$$\begin{gathered} R(X^{(k)},U^{(k)},{\tau }^{(k)},{\tau }_q^{(k)},t_0,t_f,p):= \\ \sum _{\iota =1}^{N_q^{(k)}}{w}_{\iota }^{(k)}\Vert \dot{\tilde{x}}(\mathcal{Z},t_{q_{\iota }}^{(k)})-f(\tilde{x}(\mathcal{Z},t_{q_{\iota }}^{(k)}),\tilde{u}(\mathcal{Z},t_{q_{\iota }}^{(k)}),t_{q_{\iota }}^{(k)},p){\Vert }_2^2 \end{gathered}$$

$$\underset{x,u,t_f}{\min }{t}_f+0.01{\int }_0^{t_f}{u}_1(t{)}^2+u_2(t{)}^{2}dt.$$

Full text

©2019 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

Solving optimal control problems (OCPs) are central in the field of trajectory optimization and real-time optimization-based control. In practice, most OCPs need to be solved with numerical schemes. In many situations, indirect methods can be difficult to implement, since they require analytic expressions of optimality conditions, which can be hard to derive. Direct methods have consequently become the standard for solving practical optimal control problems [1].

Direct methods often require transcription of the infinite-dimensional OCP into a nonlinear programming (NLP) problem of a finite dimension, via the introduction of a discretization mesh. If only the control trajectories are discretized as decision variables, the method is called sequential, of which direct single shooting is an example. In simultaneous approaches, such as direct multiple shooting and direct collocation, the state trajectories are also discretized and included as decision variables. A comprehensive comparison of different solution strategies are available in [2, 3], demonstrating the advantages of simultaneous methods in solving large-scale problems using sparse NLP formulations, compared to sequential methods.

Since the NLP solver will only return samples of the solution at a finite number of points in the discretization mesh, additional steps and care must be taken when representing the continuous-time results for time instances other than these sampled points. The common practice today is to directly use interpolation schemes, which are selected in accordance with the type of discretization mesh [4, 5].

Collocation methods with direct interpolation have a major drawback: the residual error of the ordinary differential equations (ODEs) is forced to zero only at collocation points. In general, no guarantees on accuracy and constraint satisfaction can be derived for system trajectories in-between collocation points. Posterior analysis is needed to identify intervals where errors are high and mesh refinement procedures are used to modify the discretization mesh. The problem has to be solved iteratively until all the errors are within tolerances.

We present a solution representation method that can significantly improve the solution accuracy for the same discretization mesh, so that results of higher quality are obtainable with relatively coarse meshes. This is achieved by minimizing the ODE residual error integrated over the solution trajectory. Section II provides the background information for solving optimal control problems numerically with the direct collocation method. Section III introduces the fundamental concept of the residual minimization method and motivates the development of the proposed scheme, which is presented in Section IV. The benefits of the method are demonstrated in Section V with two example problems, followed by concluding remarks in Section VI.

II Numerical Optimal Control

Optimization-based control often requires the solution of OCPs 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}}$).

II-A Direct collocation methods

Direct collocation methods can be categorized into fixed-order $h$ methods (e.g. Euler, Trapezoidal, and Hermite-Simpson (H-S) as in [4]), and variable higher-order $p$/$hp$ methods (e.g. Legendre-Gauss-Radau (LGR) as in [6]). Here, we aim to provide a high level overview. With a mesh of size $N:=\sum_{k=1}^{K}N^{(k)}$, the states can be approximated as

[TABLE]

within mesh interval $k$ $\in$ $\{1,\ldots,K\}$, where $N^{(k)}$ is 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 OCP (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 or structure-exploiting NLP solvers.

II-B Representing the results

The NLP solver generates a discretized solution $\mathcal{Z}\coloneqq(X,U,p,t_{0},t_{f})$ as sampled data points. Interpolating splines may be used to construct an approximation of the continuous-time solution $\tilde{z}(t)\coloneqq(\tilde{x}(\mathcal{Z},t),\tilde{u}(\mathcal{Z},t),t,p)$, with $\tilde{x}(\mathcal{Z},\cdot)$, $\tilde{u}(\mathcal{Z},\cdot)$ the approximated state and input trajectories.

II-B1 Representation via direct interpolation

Conventionally, the interpolation of the solution corresponds to the discretization scheme used in the transcription process. Thus, we must analyze how the state approximation (2) enters the optimal control problem formulation (3a).

It is not difficult to discover that the only dependency on the basis function in (3a) appears in the defect constraint (3b), through the first term representing the numerical differentiation of the approximated function $\bar{x}(\cdot)$.

For most commonly-used numerical schemes, the numerical differentiation formulation has an equivalent integration form. Both forms are presented in Table I, where $h_{k}:=\Delta t(\tau_{N}^{(k)}-\tau_{1}^{(k)})$, $\Delta t:=t_{f}-t_{0}$, and

[TABLE]

For each numerical scheme, direct interpolation of the OCP solution is possible using splines with the type and order in accordance with Table II.

For example, with Hermite-Simpson transcription, the reconstructed state trajectory inside mesh interval $k$ using cubic splines will be

[TABLE]

the dynamics trajectory with quadratic splines will be

[TABLE]

and the control trajectory with quadratic splines will have the expression

[TABLE]

for all $t\in[t_{1}^{(k)},t_{3}^{(k)}]$. The whole trajectory $\tilde{x}(\mathcal{Z},\cdot)$, $\dot{\tilde{x}}(\mathcal{Z},\cdot)$ and $\tilde{u}(\mathcal{Z},\cdot)$ can then be expressed as piecewise polynomials.

For $p/hp$ methods, Lagrange interpolating polynomials are often used as basis functions during the transcription process. An alternative version, namely barycentric Lagrange interpolation, is often used instead for solution interpolation, due to its improved numerical stability.

II-B2 Evaluation of errors

The quality of the interpolated solution needs to be assured through error analysis, assessing the level of accuracy and constraint satisfaction. Firstly, any valid trajectory $\tilde{z}(\cdot)$ must satisfy the system dynamics (1b) with a good level of accuracy. Therefore, one measure for the error due to discretization and interpolation is through the calculation of the ODE residual $\varepsilon_{r}(t)\in\mathbb{R}^{n}$ defined as

[TABLE]

For the discretized problem, the error in the state variables over each interval in-between collocation points can then be estimated with the integral

[TABLE]

as a single metric for a multi-variable problem, or

[TABLE]

for each dynamics equation separately. $\eta\in\mathbb{R}^{N}$ or $\sigma\in\mathbb{R}^{n\times N}$ are typically referred to as the absolute local error [4]. The operator $\|\cdot\|_{2}$ is the vector 2-norm. The integral can be practically estimated by high order quadrature.

In addition, numerical discretization inevitably leads to possible constraint violations of the trajectories in-between the collocation points. For path and box constraints that are expressed semi-explicitly as (1c), the absolute local constraint violation $\varepsilon_{c_{\zeta}}(t)\in\mathbb{R}^{n_{g}}$ may be straight-forwardly estimated by

[TABLE]

Once the distributions of errors are calculated, appropriate modifications can be made to the discretization mesh, to iteratively resolve the problem until the obtained solution fulfills all predefined error tolerances ($\eta_{tol}$ and $\varepsilon_{c_{tol}}$). This process is called mesh refinement (MR). Common approaches for mesh refinement include adding intervals and/or changing the polynomial order. The NLP formulated based on the new mesh is warm started using the previous solution from the coarser mesh. This can often lead to significantly faster convergence, compared to a fine uniform mesh without MR, thus reducing the overall computation time.

II-C Problems associated with direct reconstruction

Practical experience has shown that trajectory interpolation in accordance with the discretization scheme is not the best choice. In many cases large discretization errors and constraint violations occur inside the intervals in-between collocation points. Furthermore, if the optimal control trajectory is discontinuous, direct interpolation using polynomials can often result in a Gibbs-like phenomenon, inducing non-physical oscillations in the solution. A typical example for this to happen would be in problems with bang-bang control, if the switch happens in-between mesh intervals.

These issues are fundamentally rooted in the direct collocation formulation. Firstly, states, dynamics and controls can rarely all be approximated accurately by polynomials. Even in the simple case where $f(x(t),u(t))=\dot{x}(t)=ax(t)+u(t)$ and $u(t)=1$ are both polynomials (thus can be represented exactly by polynomials), the corresponding state trajectory $x(t)=x(0)e^{at}+\int_{0}^{t}e^{a(t-s)}u(s)\ ds$ is clearly not a polynomial and approximation errors should be expected.

It is then important to note that driving the defect constraint (3b) to zero (or machine precision) at collocation points does not imply that the polynomial functions used for the state and input approximations in the NLP will satisfy the dynamic equations and constraints in-between collocation points. In fact, the opposite can and often does occur.

It is well-known in the field of curve fitting that if a function cannot be exactly represented by a polynomial, forcing the polynomial to exactly go though some sampled data points generally results in larger errors in comparison to fitting using least squares criteria. The same analogy can be applied here: forcing the defect constraints to be zero at collocation points will generally result in larger overall defect errors for the whole trajectory, in comparison to a method that minimizes the integral of the defect errors in a least squares manner. This observation motivated the development of the integrated residual minimization method.

III Method of integrated residual minimization

Integrated residual minimization is motivated by the recently-proposed method in [7], which is a generalization of the least-squares approach for solving differential equations to solving dynamic optimization problems. The idea is that instead of forcing the ODE residuals (8) to be zero at collocation points with (3b), the method tries to minimize the square of the 2-norm of the ODE residuals for the represented solution polynomials integrated along the whole trajectory, i.e.

[TABLE]

with

[TABLE]

As presented in [7], the expressions for state and input functions $\hat{x}$ and $\hat{u}$ can be polynomials of any standard types, with polynomial coefficients $P_{j,q}$ as decision variables.

This concept shares some resemblance with the method of direct error enforcement in [8], where the residual errors are computed with a different formulation and only at a single non-collocation reference point for each mesh interval.

Remark 1

The choice of representation in [7] with $P_{j,q}$ as decision variables increases the computational complexity of the problem in comparison to direct collocation:

• •

One extra decision variable is required for every state and input variable in every mesh segment.

• •

Simple bounds need to be implemented as general inequality constraints.

• •

Additional computations are required to obtain the initial guesses of the decision variables from an estimation of the solution trajectory.

• •

The magnitudes of decision variables may span a wide numerical range. This is detrimental in terms of ensuring consistent numerical accuracy in computations.

• •

If finite differences are used for obtaining the derivative information, the calculations can be less accurate.

• •

Proper scaling of decision variables can be difficult.

• •

State continuity in-between mesh segments might need to be enforced with additional equality constraints.

Therefore, we need to develop a method that avoids the above-listed drawbacks and, to a great extent, retains the computational efficiency of direct collocation.

IV The proposed scheme

Based on the above observations, we propose a method to generate solution trajectories that can be orders of magnitude more accurate than direct interpolation for the ODE defect error, without increasing the size of the discretization mesh. The method retains the same decision variables as in (3), namely $\mathcal{Z}\coloneqq(X,U,p,t_{0},t_{f})$, and uses the interpolation polynomial formula $\tilde{x}(\mathcal{Z},\cdot)$, $\dot{\tilde{x}}(\mathcal{Z},\cdot)$ and $\tilde{u}(\mathcal{Z},\cdot)$ to directly represent $\hat{x}(\cdot)$, $\dot{\hat{x}}(\cdot)$ and $\hat{u}(\cdot)$ in (9) and (10).

For example, consider Hermite-Simpson discretization. The input trajectory inside mesh interval $k$ can be represented by the polynomial as in (7), based on the values of the decision variables $U_{1}^{(k)}$, $U_{2}^{(k)}$ and $U_{3}^{(k)}$. However, for the state trajectory, one challenge arises. For solutions to (3), continuity of state variables are automatically fulfilled when using (5) as the interpolation equation; however, this is not generally the case for arbitrary solutions that violate the defect constraint (3b).

To avoid imposing additional path constraints for state continuity, we make use of the original Hermite-Simpson numerical integration scheme (in Table I), and obtain

[TABLE]

As a check, substituting $t=t_{2}^{(k)}=t_{1}^{(k)}+h_{k}/2$ and (11) into (5) will result in $X_{2}^{(k)}$, and substituting $t=t_{3}^{(k)}=t_{1}^{(k)}+h_{k}$ and (12) into (5) will result in $X_{3}^{(k)}$. Thus, with $F_{1}^{(k)}$, $F_{2}^{(k)}$ and $F_{3}^{(k)}$ calculated based on (4), (11) and (12), respectively, the interpolation formula (5) guarantees state trajectory continuity without imposing additional constraints.

Thus, an optimization problem for representing the OCP solution can be formulated as

[TABLE]

with (13b) the constraint for the objective of the original OCP (3) and $J_{c}\in\mathbb{R}$ the value of the cost obtained from direct collocation. $R$ is the residual cost: for certain problems, this can be calculated precisely with analytical expressions; for most practical problems, quadrature rules of sufficiently high order can be used, i.e.

[TABLE]

with $t_{q_{\iota}}^{(k)}:=\frac{t_{f}^{(k)}-t_{0}^{(k)}}{2}\tau_{q_{\iota}}^{(k)}+\frac{t_{f}^{(k)}+t_{0}^{(k)}}{2}$ and $\tau_{q}^{(k)}\in\mathbb{R}^{N_{q}^{(k)}}$ the quadrature mesh for approximating the integral inside a mesh interval, where $w_{\iota}^{(k)}$ are the corresponding quadrature weights. Typically a Gaussian quadrature of order $N_{q}^{(k)}\geq 4N^{(k)}+1$ is required for a good accuracy [7].

When comparing (13) to the original OCP formulation of direct collocation (3), it is straightforward to note that, if $\mathcal{Z}$ is a solution to (3), then $\mathcal{Z}$ will also be a feasible point for (13). If no better feasible solution is obtainable, (13) is at least guaranteed to have one solution, namely $\mathcal{Z}$.

In terms of computational complexity, the proposed formulation (13) avoids the shortcomings as listed in Remark 1. In addition, unlike the penalty-barrier finite element (PBF) method proposed in [7], which requires tailored solvers for a good performance, (13) can be efficiently solved with the same off-the-shelf sparse NLP solvers as direct collocation. The transcription and the majority of the computational components can be shared between the two, and warm starting techniques can be exploited to accelerate the computations.

V Example Problems

Here, we present two example problems to demonstrate the main advantages of the proposed scheme. Both OCPs are transcribed using the optimal control software ICLOCS2 [9], and numerically solved to a tolerance of $10^{-9}$ with NLP solver IPOPT [10] (version 3.12.9). With extremely coarse meshes, the emphasis of the comparison will not be on yielding solutions that look similar to the optimal trajectory. Instead, the goal is to obtain sub-optimal solutions that, when applied, can result in low discrepancies between the represented solution and the implementation outcome.

V-A Two-Link Robot Arm

The two-link robot arm problem presented here was adapted from [11, Ex. 2, Sect. 12.4.2]. Consider a system consisting of two identical beams with the same property (mass: $m=1$ kg, length: $l=1$ m, and moment of inertia), connected at two actuated joints. The objective is to reposition a payload of mass $M=1$ kg in minimum time, with the addition of a regularization term:

[TABLE]

The system has angular rates $\omega_{\phi}$, $\omega_{\psi}$, and angles $\phi$, $\chi=\phi-\psi$ as state variables, and nondimensionalized torque $u_{1}$ and $u_{2}$ as inputs. Furthermore, the variable simple bounds and boundary conditions are imposed in accordance to the reference, except that $\chi(t_{f})=0.5$ rad, and $\phi(t_{f})=0.522$ rad.

Figures 1 and 2 illustrate the solutions to the two-link robot arm problem problem generated with the two different solution representation methods.

Presented alongside are the outcomes from the actual implementation of the resultant input trajectory on the same dynamic model, solved with a non-stiff variable-order ODE solver (Matlab ode113) with a time step 100 times smaller than the discretization grid of the optimization problem. Observe that:

• •

Despite a very small tolerance and successful termination of the NLP solver, the collocation solution and interpolation of the solution exhibit large errors, leading to significant deviations to the state trajectories when the inputs are directly applied. In contrast, only minor discrepancies can be observed for the solutions represented using integrated residual minimization, on the same coarse grid with relatively low-order discretization.

• •

Although the constraints are implemented in the exact same way, the proposed method, to a greater extent, alleviates the issues of constraint violations inside the mesh intervals. This is because these constraint violations are often related to the large ODE defect errors in-between collocation points, which are directly dealt with by the residual minimization scheme.

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

Based on previous developments [12], 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. See the illustration in Figure 3.

A simplified windshear model is used with wind speed contributions represented by a horizontal and a vertical component. Other details about the aerodynamic modelling, parameter values, simple bounds and boundary conditions are the same as in [4]. A static parameter $h_{min}$ is introduced to represent the minimum altitude. The objective is therefore to minimize $-h_{min}$ together with path constraint $h(t)\geq h_{min}$.

The angle of attack $\alpha$ is the actual control input to the physical system (aircraft); however, in order to implement a constraint on its rate of change, $\nu$ is introduced as angle of attack rate and serves as the control input with $\dot{\alpha}(t)=\nu(t)$. This implementation is known to exhibit singular arc behaviour [4, 13, 14], leading to fluctuations and ringing phenomena in the solutions.

The solutions to this problem are collectively shown in Figures 4 and 5.

In addition to the advantages identified in the previous example, solution representation via integrated residual minimization have clear benefits in suppressing fluctuations. This ringing phenomenon is frequently observed in direct collocation solutions of singular control problems, as well as the directly interpolated solution trajectories.

VI Conclusions

Although interpolation of collocation solutions can sometimes yield good results, it is very difficult to guarantee accuracy without posterior error assessments and mesh design iterations. As shown by the examples, despite successfully solving the NLP to negligibly small tolerance, the validity of the solution may still be questionable with large discrepancies. The flaws are rooted in collocation schemes, where the ODE defect errors are forced to zero at collocation points, regardless of the errors inside the intervals.

The proposed solution representation method of integrated residual minimization fundamentally addresses this shortcoming by instead minimizing the integrated ODE residual error along the whole trajectory. As a result, solutions of higher accuracy are obtainable with the same discretization mesh, allowing the mesh to be relatively coarse. This benefit is clearly demonstrated with the example problems: despite being highly nonlinear, moderately complex and solved on a coarse low-order mesh, only minor differences are observed between the represented solution and actual implementation.

Since solving OCPs with numerical methods are essentially multi-objective optimization problems, one will inevitably face the trade-off between minimizing the objective (for optimality) and minimizing the discretization errors (for accuracy). By utilizing the objective value from a collocation solution, our proposed approach is yet to offer complete flexibility for the user in managing this trade-off process. Further development of this concept into a standalone method could potentially offer more benefits.