Fast Mesh Refinement in Pseudospectral Optimal Control

TL;DR

This paper introduces a novel Birkhoff interpolation-based mesh refinement method for pseudospectral optimal control that maintains low condition numbers and spectral accuracy even at very high polynomial orders, enabling efficient solutions to complex problems.

Contribution

It develops a Birkhoff interpolation approach to produce well-conditioned pseudospectral discretizations with condition numbers growing only as -20 words)

Findings

Condition number grows as -20 words)

Achieves spectral accuracy at high polynomial orders

Successfully solves a 1000th order orbit transfer problem

Abstract

Mesh refinement in pseudospectral (PS) optimal control is embarrassingly easy --- simply increase the order $N$ of the Lagrange interpolating polynomial and the mathematics of convergence automates the distribution of the grid points. Unfortunately, as $N$ increases, the condition number of the resulting linear algebra increases as $N^2$; hence, spectral efficiency and accuracy are lost in practice. In this paper, we advance Birkhoff interpolation concepts over an arbitrary grid to generate well-conditioned PS optimal control discretizations. We show that the condition number increases only as $\sqrt{N}$ in general, but is independent of $N$ for the special case of one of the boundary points being fixed. Hence, spectral accuracy and efficiency are maintained as $N$ increases. The effectiveness of the resulting fast mesh refinement strategy is demonstrated by using \underline{polynomials of over a thousandth order} to solve a low-thrust, long-duration orbit transfer problem.

Equations

$$t_i=\frac{1}{2}\left[(t_f+t_0)-(t_f-t_0)\cos \left(\frac{i\pi }{N}\right)\right],i=0,1,\dots ,N$$

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

$$\begin{array}{rlrlr}\mathbb{X}& ={\mathbb{R}}^{N_x}& \mathbb{U}& ={\mathbb{R}}^{N_u}& \\ \mathbf{x}& =(x_1,\dots ,x_{N_x})& \mathbf{u}& =(u_1,\dots ,u_{N_u})& \end{array}\}\text{(preamble)}$$

$$\displaystyle\begin{aligned} \overbrace{(G)}^{\text{\normalsize problem}}\left\{\begin{array}[]{ll}\text{Minimize }&\\[-12.50002pt] &\left.\begin{aligned} \makebox[76.3136pt][r]{$\displaystyle J[\boldsymbol{x}(\cdot),{\boldsymbol{u}}(\cdot),t_{0},t_{f}]:=$}&\\ \makebox[76.3136pt][r]{$\displaystyle E(\boldsymbol{x}_{0},\boldsymbol{x}_{f},t_{0},t_{f})$}&+\makebox[83.50421pt][l]{$\displaystyle\int_{t_{0}}^{t_{f}}F(\boldsymbol{x}(t),{\boldsymbol{u}}(t),t)\ dt$}\\ \end{aligned}\right\}\ \text{(cost)}\\[20.00003pt] \text{Subject to}&\\[-15.00002pt] &\left.\begin{aligned} \makebox[76.3136pt][r]{$\displaystyle\dot{\boldsymbol{x}}$}&=\makebox[83.50421pt][l]{$\displaystyle{\boldsymbol{f}}(\boldsymbol{x}(t),{\boldsymbol{u}}(t),t)$}\\ \end{aligned}\right\}\ \text{(dynamics)}\\[10.00002pt] &\left.\begin{aligned} \makebox[76.3136pt][r]{$\displaystyle{\boldsymbol{e}}^{L}$}&\leq\makebox[83.50421pt][l]{$\displaystyle{\boldsymbol{e}}(\boldsymbol{x}_{0},\boldsymbol{x}_{f},t_{0},t_{f})\leq{\boldsymbol{e}}^{U}$}\\ \end{aligned}\right\}\ \text{(events)}\\[10.00002pt] &\left.\begin{aligned} \makebox[76.3136pt][r]{$\displaystyle{\boldsymbol{h}}^{L}$}&\leq\makebox[83.50421pt][l]{$\displaystyle{\boldsymbol{h}}(\boldsymbol{x}(t),{\boldsymbol{u}}(t),t)\leq{\boldsymbol{h}}^{U}$}\\ \end{aligned}\right\}\ \text{(path)}\end{array}\right.\end{aligned}$$

$$x\in \mathbb{R},u\in \mathbb{R},\tau \in [-1,1]$$

$$\displaystyle(\textsf{$P$})\left\{\begin{array}[]{lrl}\emph{Minimize }&J[x(\cdot),u(\cdot)]=&E(x(-1),x(1))\\ \emph{Subject to}&\dot{x}(\tau)=&f(x(\tau),u(\tau))\\ &e(x(-1),x(1))=&0\end{array}\right.$$

$$x(\tau )=x^{\infty }(\tau ):=\underset{N\to \infty }{\lim }{x}^N(\tau ):=\underset{N\to \infty }{\lim }\sum _{m=0}^N{a}_m{P}_m(\tau )$$

$$x^N(\tau ):=\sum _{j=0}^N{x}^N({\tau }_j)L_j(\tau )$$

$$L_j({\tau }_i)={\delta }_{ij}$$

$${\pi }^N:=[{\tau }_0,{\tau }_1,\dots ,{\tau }_N]$$

$$-1\le {\tau }_0<{\tau }_1<\cdots <{\tau }_{N-1}<{\tau }_N\le 1$$

$$x^N(\tau ):=\sum _{j=0}^N{x}_j{L}_j(\tau )$$

$${\dot{x}}^N({\tau }_i)=\sum _{j=0}^N{x}_j{\dot{L}}_j({\tau }_i)$$

$${\boldsymbol{D}}:=\Big{(}\dot{L}_{j}(\tau_{i})\Big{)}_{0\leq i,j\leq N}$$

$$X:=[x_0,x_1,\dots ,x_N{]}^T$$

$$\dot{X}=\mathbf{D}X$$

$$u^N(\tau )=\sum _{j=0}^N{u}_j{C}_j(\tau )$$

$$C_j({\tau }_i)={\delta }_{ij}$$

$$f(x^{N}(\tau),u^{N}(\tau))\big{|}_{\pi^{N}}=[f(x_{0},u_{0}),f(x_{1},u_{1}),\ldots,f(x_{N},u_{N})]$$

$$f(X,U):=[f(x_0,u_0),f(x_1,u_1),\dots ,f(x_N,u_N){]}^T$$

$$U:=[u_0,u_1,\dots ,u_N{]}^T$$

$$X\in {\mathbb{R}}^{N+1},U\in {\mathbb{R}}^{N+1}$$

$$\displaystyle\left(\textsf{$P^{N}$}\right)\left\{\begin{array}[]{lrl}\emph{Minimize }&J^{N}[X,U]:=&E(x_{0},x_{N})\\ \emph{Subject to}&{\boldsymbol{D}}X=&f(X,U)\\ &e(x_{0},x_{N})=&0\end{array}\right.$$

$$\text{if }\frac{d}{dt}(z)$$

$$\text{if }\frac{d^2}{d{t}^2}(z)$$

$$\mathbf{M}\mathbf{D}X=\mathbf{M}f(X,U)$$

$${\tilde{x}}^N(\tau ):=x({\tau }_0)B_0(\tau )+\sum _{j=1}^{N-1}\ddot{x}({\tau }_j)B_j(\tau )+x({\tau }_N)B_N(\tau )$$

$${\tilde{x}}^N({\tau }_0)=x({\tau }_0),{\tilde{x}}^N({\tau }_N)=x({\tau }_N),{\ddot{\tilde{x}}}^N({\tau }_j)=\ddot{x}({\tau }_j),j=1,\dots ,N-1$$

$$B_0({\tau }_0)$$

$$B_0({\tau }_N)$$

$${\ddot{B}}_0({\tau }_i)$$

$${\overline{x}}^N(\tau ):=x({\tau }_0)B_0(\tau )+\sum _{i=1}^{N-1}\dot{x}({\tau }_i)B_i(\tau )+x({\tau }_N)B_N(\tau )$$

$${\overline{x}}^N({\tau }_0)=x({\tau }_0),{\overline{x}}^N({\tau }_N)=x({\tau }_N),{\dot{\overline{x}}}^N({\tau }_i)=\dot{x}({\tau }_i)$$

Full text

1 Introduction

In principle, mesh refinement in pseudospectral (PS) optimal control[1] is embarrassingly easy — simply increase the order $N$ of the Lagrange interpolating polynomial. For instance, in a Chebyshev PS method[2, 3, 4], the Chebyshev-Gauss-Lobatto (CGL) mesh points are given by[2],

[TABLE]

where, $[t_{0},t_{f}]$ is the time interval. Grid points generated by (1) for $N=5,10$ and $20$ over a canonical time-interval of $[0,1]$ are shown in Fig. 1.

The production of these mesh points is as fast as it takes to compute the cosine function. Furthermore, the mesh points are at optimal locations for a Chebyshev PS method because their placement is dictated by well-established rules of approximation theory[5]. The price for this speed and simplicity is an increase in the condition number of the underlying linear algebra that governs the iterations[1, 2, 6]. In fact, as shown in Fig. 2,

the condition number associated with the PS discretization of the dynamical equations,

[TABLE]

increases as $N^{2}$. The plots in this figure are the condition numbers of the non-singular, inner, square-portion of the various differentiation matrices[1, 2] associated with the indicated mesh selections. Note that the Legendre-Gauss-Radau mesh has the worst performance of all the indicated PS grid selections. An intuitive approach to manage the growth in the condition number is to break up the time interval into several segments and apply a PS method over each segment[7, 8, 9, 10]. When the segments are non-overlapping and information across the segments is passed across a single point, it results in a PS knotting method[8, 10]; see Fig. 3.

The special nodes where information across segments is transmitted are called PS knots. In the spectral algorithm for PS methods[10, 11], the number and location of the PS knots are determined adaptively. The adaptation process is based on an analysis of the solution generated by the prior iterate. Among other things, the analysis involves an estimation of the error and a determination of the frequency content of the system trajectory. Based on this analysis, either $N$ per segment is increased or new knots are placed near high frequency regimes[10, 11]. The entire process is repeated for the next iterate or until all stopping criteria are met[1, 8, 10, 11].

Although the PS knotting method manages the growth in the condition number by simply limiting the order of the approximating polynomial in a given segment, the price for this approach is quite severe. First, there is an increase in the computation time to determine the number and location of the knots; second, there is a significant reduction in the rate of convergence[12, 13, 14, 15] that adds to the increase in computational time and effort. The convergence rate varies inversely to the number of segments — the more the segments, the poorer the convergence rate. If the order of the approximating polynomial is limited to low orders (e.g., $N=4$), a PS method becomes theoretically equivalent to a Runge-Kutta method[16, 17, 18]. In fact, the “exponential” rate of convergence of an unknotted PS method is chiefly due to the high orders (i.e., $N\gg 10$) of the approximating polynomials[2, 5]. In many aerospace applications[19, 20, 21, 22], including flight implementations[1, 23, 24, 25, 26, 27, 28], very accurate solutions[29, 30] are generated for $N\leq 100$. In certain emerging areas of applications[31, 32, 33], there have been new requirements to generate solutions for $N\gg 100$. Because the unknotted PS method generates high condition numbers, various types of preconditioning are often used to alleviate the problem[34, 35]. In recent years, there have been major advancements in well-conditioned PS methods[36, 37, 38] for ordinary differential equations. **These methods offer a remarkable drop in the condition number from $\mathcal{O}(N^{2})$ to $\mathcal{O}(1)$. **In this paper, we adapt these advancements for PS optimal control using Birkhoff interpolants[36]. The dramatic drop in the condition number allows us to generate solutions for thousands of grid points (i.e., using approximating polynomials for state trajectories in excess of a thousand degrees) thereby generating a very fast mesh refinement strategy. Numerical results for a sample low thrust trajectory optimization problem are presented to illustrate the advancements.

The citation for the journal version of this paper is: J. Guid., Contr. & Dyn., 42/4, 2019, pp. 711–722.

2 A Brief Review of the State-of-the-Art in PS Optimal Control Methods

In order to provide a proper context for the issues in PS optimal control methods, we first briefly review the state of the art. To this end, consider a generic optimal control problem given by[39, 40, 41],

[TABLE]

where, the symbols denote the following:

• •

$\mathbb{X}$ and $\mathbb{U}$ are $N_{x}$- and $N_{u}$-dimensional real-valued state and control spaces respectively. We assume $N_{x}\in\mathbb{N}^{+}$ and $N_{u}\in\mathbb{N}^{+}$.

• •

$J$ is the scalar cost function. The arguments of $J$ are the optimization variables.

• •

The optimization variables are:

• –

${\boldsymbol{x}(\cdot)}$: the $N_{x}$-dimensional state trajectory,

• –

${{\boldsymbol{u}}(\cdot)}$: the $N_{u}$-dimensional control trajectory,

• –

$t_{0}$: the initial clock time, and

• –

$t_{f}$: the final clock time.

• •

$E$ is the scalar endpoint cost function. The arguments of $E$ are the endpoints. In the classical literature, $E$ is known as the “Mayer” cost function.

• •

The endpoints are the initial state $\boldsymbol{x}_{0}\equiv\boldsymbol{x}(t_{0})$, the final state $\boldsymbol{x}_{f}\equiv\boldsymbol{x}(t_{f})$, the initial time $t_{0}$ and the final time $t_{f}$.

• •

$F$ is the scalar running cost function. The arguments of $F$ are the instantaneous value of the state variable $\boldsymbol{x}(t)$, the instantaneous value of the control variable ${\boldsymbol{u}}(t)$ and time $t$.

• •

${\boldsymbol{f}}$ is the $N_{x}$-dimensional “dynamics function,” or more appropriately the right-hand-side of the dynamics equation. The arguments of ${\boldsymbol{f}}$ are exactly the same as the arguments of $F$.

• •

${\boldsymbol{e}}$ is the $N_{e}$-dimensional endpoint constraint function. The arguments of ${\boldsymbol{e}}$ are exactly the same as that of $E$.

• •

${\boldsymbol{e}}^{L}$ and ${\boldsymbol{e}}^{U}$ are the $N_{e}$-dimensional lower and upper bounds on the values of ${\boldsymbol{e}}$.

• •

${\boldsymbol{h}}$ is the $N_{h}$-dimensional path constraint function. The arguments of ${\boldsymbol{h}}$ are exactly the same as that of $F$.

• •

${\boldsymbol{h}}^{L}$ and ${\boldsymbol{h}}^{U}$ are the $N_{h}$-dimensional lower and upper bounds on the values of ${\boldsymbol{h}}$.

As shown in [42], many key aspects of PS optimal control can be understood by considering a distilled optimal control problem given by,

[TABLE]

The PS results generated for Problem $P$ easily generalize to Problem $G$, albeit with some tedious bookkeeping. In order to avoid the distractions of bookkeeping, we limit many of our subsequent discussions to Problem $P$ while noting that all of the results – with appropriate reindexing – apply to Problem $G$ as well. See also [1] for a more thorough review and survey of the key results.

In PS optimal control theory, the state trajectory $x(\cdot)$ is expressed – in principle – as a generalized Fourier expansion,

[TABLE]

where, $a_{m}\in\mathbb{R}$ is an unknown “spectral” coefficient of an $m^{th}$-order polynomial $P_{m}$. Note that (7) is an exact representation and not an approximation. The justification for (7) comes from the classic Stone-Weierstrass theorem[43] which guarantees a polynomial representation for $x(\cdot)$ under mild assumptions. For most aerospace applications, the “best” choices for $P_{m}$ are the Legendre and Chebyshev polynomials; hence, PS optimal control techniques resulting from these “big two” polynomials are known as the Legendre and Chebyshev PS methods respectively. A key principle in a generic PS method is that the spectral coefficients $a_{m}$ are computed indirectly by transforming (7) to an equivalent “pseudospectral” form[39, 44, 45],

[TABLE]

where, $\tau_{j},j=0,1,2,\ldots,N$ are discrete points in time, known as nodes, and $L_{j}$ are Lagrange interpolating polynomials. That is, $L_{j}(\tau),j=0,1,2,\ldots,N$ are such that they satisfy the Kronecker relationship,

[TABLE]

Equations (7) and (8) are known as the modal and nodal representations of $x(\cdot)$ respectively[39].

Let,

[TABLE]

be an arbitrary grid (or “mesh”) of points (see Fig. 4) such that

[TABLE]

Note that (11) allows the extreme points of the grid, $\tau_{0}$ and $\tau_{N}$, to be at or near the endpoints. A nodal representation of the state trajectory over $\pi^{N}$ is given by,

[TABLE]

Applying (9) to (12), it follows that $x_{i}=x^{N}(\tau_{i}),\ i=0,\ldots,N$. Differentiating (12) and evaluating $x^{N}(t)$ over the grid we get,

[TABLE]

Define

[TABLE]

as the PS differentiation matrix, and let,

[TABLE]

Then, using (14) and (15), (13) can be written as a matrix-vector product given by,

[TABLE]

Note that in sharp contrast to Runge-Kutta methods, PS methods rely almost entirely on an accurate discretization of $\dot{x}$ rather than $f(x,u)$, the right-hand-side of the differential equation, (2).

A nodal representation of the control trajectory over $\pi^{N}$ is given by,

[TABLE]

where, $C_{j}(\tau)$ is some interpolating function.

Remark 1

Note that the control trajectory $\tau\mapsto u^{N}(\tau)$ is not necessarily a polynomial.111Unfortunately, it is widely and erroneously reported in the optimal control literature that $C_{j}$ is a polynomial. That is, $C_{j}(\tau)$ is any function that satisfies,

[TABLE]

The non-polynomial form of $t\mapsto u^{N}(t)$ is exploited in [12, 13, 14, 15] for a strong proof of convergence; that is, a convergence proof that does not rely on unnecessary assumptions[46] that are typically not satisfied in many practical optimal control problems. Note as well that the Bellman pseudospectral method[47, 48] also capitalizes on the non-polynomial form of $t\mapsto u^{N}(t)$.

Substituting (12) and (17) in $f(x,u)$, and evaluating $f$ over $\pi^{N}$ we get,

[TABLE]

It is convenient to reuse $f$ as an overloaded function (see [39], page 8) defined by,

[TABLE]

where,

[TABLE]

Setting $\dot{X}=f(X,U)$ and using (16), we define Problem $P^{N}$ according to,

[TABLE]

Comparing (6) and (22), it is clear that Problem $P^{N}$ is strikingly similar to Problem $P$. In fact, one may view the production of Problem $P^{N}$ as simple as the process of replacing the continuous operator $d/dt$ in Problem $P$ by the matrix ${\boldsymbol{D}}$.

Remark 2

In the derivation of Problem $P^{N}$, no assumption was made on the structure of the grid $\pi^{N}$; hence, (22) represents a discretization of Problem $P$ for all PS methods[45].

Remark 3

Once $\pi^{N}$ is chosen for a given Problem $P$, the theoretical performance of a PS method depends only on the selection of the grid. The worst performance (i.e., rapid divergence) is characterized by a uniform distribution of nodes while the best performance (i.e., rapid convergence) requires node clustering at the end points[45].

Remark 3 is illustrated in Fig. 5.

This figure was generated using the data and results from [45]. The parameter $\beta$ (see Fig. 5 a)) denotes a measure of clustering: from optimal clustering (i.e., when $\beta=0$) to a uniform node distribution (i.e., when $\beta=1$). It is apparent from Fig. 5 b) that optimal clustering is not required for convergence but is indeed critical for the rate of convergence. Hence, it is desirable to maintain node clustering for a fast mesh refinement.

Remark 4

If $\pi^{N}$ is a convergent grid, many of the properties of the continuous-time operator $d/dt$ are faithfully represented by ${\boldsymbol{D}}$. For instance, let $z$ be a function that is not identically equal to zero; then,

[TABLE]

From (23a), it follows that ${\boldsymbol{D}}$ must be singular to accurately represent $d/dt$. Note also that ${\boldsymbol{D}}$ must be square so that higher-order derivatives can be properly represented analogous to a repetitive application of $d/dt$ as implied in (23b).

Remark 5

Node clustering is necessary for polynomial basis functions. It is possible to reduce node clustering through the use of non-polynomial basis functions[2, 49, 50]; however, the condition number of the resulting linear algebra still grows with an increase in $N$.

Problem $P^{N}$ is solved by the spectral algorithm[1, 11]. A key procedure for mesh refinement in the spectral algorithm for PS methods can be described as follows:

Select an initial sequence of integers $N_{0}<N_{1}<\cdots<N_{m}$ and generate a sequence of problems $P^{N_{k}},k=0,1,\ldots,m$. 2. 2.

Solve Problem $N_{0}$ using the relaxed elastic programming concept described in [11]. Set $i=1$. 3. 3.

Solve Problem $N_{i}$ using the results from Problem $N_{i-1}$ by a joint application of PS convergence theory for optimal control[13, 14, 15] and the covector mapping principle[39]. 4. 4.

Evaluate the (decay of the) spectral coefficients for convergence and stopping criteria[10, 11]. Exit if all stopping criteria are met; else, set $i=i+1$, and go to Step 3.

Many of the advancements in the spectral algorithm are implemented in DIDO©[30, 39], a state-of-the-art MATLAB® optimal control toolbox for solving optimal control problems. As noted earlier, the condition number of the linear matrix equations involved in the spectral algorithm grows as $N_{k}^{2},k=0,1,\ldots,m$; hence, there is a loss in spectral accuracy and convergence rate on those problems that require a finer mesh. Among many, one of the current remedies for this problem is to use a pre-conditioner, ${\boldsymbol{M}}$, such that the discretized differential equation in Problem $P^{N}$ can be written as,

[TABLE]

The matrix ${\boldsymbol{M}}$ must be chosen so that the condition number of the matrix equations resulting from the product ${\boldsymbol{M}}{\boldsymbol{D}}$ is lower than when ${\boldsymbol{M}}$ is the identity matrix. If ${\boldsymbol{M}}$ is not a diagonal matrix, then the sparsity of the right-hand-side of (24) decreases leading to an increase in the computational burden of the resulting linear algebra. Yet another remedy for managing the growth in the condition number is to use PS knots[8]. In fact, the full spectral algorithm incorporates an automatic determination of the number and location of knots based on a detection of rapid changes in the control function [10, 11]. As noted in Section I, this approach is far from satisfactory. Although pre-conditioning and PS knots can be used jointly to mitigate the side-effects of both concepts, the resulting algorithm involves too many tuning parameters and inconsistent performance results. Clearly, there has been a lot of room for improvement in PS optimal control methods, particularly for those problems that require fast and accurate solutions over a fine mesh. As will be apparent shortly, the use of Birkhoff interpolating polynomials[51, 52] offers a new PS approach where solutions over thousands of points can be generated in a fast and stable manner.

3 Birkhoff Interpolation Over an Arbitrary Grid

In contrast to a Lagrange interpolant which is based on interpolating the values of a function, a Birkhoff interpolant is based on interpolating the values of a heterogenous mix of various orders of derivatives of a function[51]. That is, a Birkhoff interpolant generalizes Lagrange and Hermite interpolants[52]. For second-order boundary value problems (BVPs) with Dirichlet boundary conditions, a Birkhoff interpolation polynomial $\widetilde{x}^{N}(\tau)$ is given by[36],

[TABLE]

where, $[\tau_{0},\tau_{1},\ldots,\tau_{N}]:=\pi^{N}$ is an arbitrary grid as before, and $B_{k},k=0,\ldots,N$ are polynomials of order $N$ or less such that, if they exist, must satisfy the interpolation requirements,

[TABLE]

Equation (26) generates the following conditions for the Birkhoff polynomials,

[TABLE]

That is, for a second-order BVP, a Birkhoff interpolating polynomial may be based on a mix of the second derivatives of a function and its values at the boundary points. Because an optimal control problem also generates a BVP (through an application of the necessary conditions), it seems reasonable to simply rewrite (25) for first-order systems so that state-space representation may be used. This simple adaptation generates a proposal for a Birkhoff interpolant given by,

[TABLE]

where, $B_{k},k=0,\ldots,N$ are polynomials (possibly different from those in (25)) that must satisfy,

[TABLE]

The conditions for (29) that are analogous to (27) are given by,

[TABLE]

The satisfaction of the conditions given by (30) has some well-known problems related to the existence and uniqueness of a solution[51, 52]; hence, we do not pursue this further. If the conditions on $B_{0}$ or $B_{N}$ are removed from (30), then it can be shown that a unique Birkhoff interpolant exists[53]. In doing the latter, a Birkhoff interpolant is given by,

[TABLE]

where, $B_{k},k=0,\ldots,N$ must now only satisfy,

[TABLE]

Equation (32) generates fewer conditions than (30); these are given by,

[TABLE]

In the same manner, an equivalent Birkhoff interpolant is given by,

[TABLE]

where, $B_{k},k=0,\ldots,N$ are yet another set of polynomials that must satisfy,

[TABLE]

The conditions for $B_{k},k=0,\ldots,N$ in (34) are given by,

[TABLE]

It is easy to show that a Birkhoff interpolant exists for both cases given by (31) and (34).

Remark 6

The Birkhoff interpolants given by (31) and (34) are over an arbitrary grid, $\pi^{N}$. Endpoint node clustering is necessary for convergence (see Fig. 5). Consequently, Birkhoff interpolants for Jacobi-Gauss-Lobatto, Jacobi-Gauss-Radau, and Jacobi-Gauss meshes follow from (31) and (34) by an appropriate choice of $\pi^{N}$.

It will be apparent shortly that it is convenient to define two subsets of the grid $\pi^{N}$ given by,

[TABLE]

Thus, $\pi^{N}$ may be represented either by $\pi^{N}=[\tau_{0},\pi^{N}_{a}]$ or $\pi^{N}=[\pi^{N}_{b},\tau_{N}]$. Furthermore, we define the matrices,

[TABLE]

where, ${\boldsymbol{B}}_{a}$ and ${\boldsymbol{B}}_{b}$ satisfy (33) and (36) respectively.

Theorem 1** (Wang-Samson-Zhao)**

Let $\omega\in\left\{a,b\right\}$. Then,

[TABLE]

where, ${\boldsymbol{I}}_{N}$ is an $N\times N$ identity matrix.

Proof 3.2**.**

The proof of this result for $\omega=a$ is implied in [36] by Theorems 3.2 and 4.1. For the purposes of completeness and clarity, we prove this theorem for an arbitrary grid by the same procedures used in [36]. Furthermore, because the matrix ${\boldsymbol{B}}_{b}$ is not presented in [36], we first prove this result for $\omega=b$. The proof for $\omega=a$ follows by a trivial extension.

For any $N^{th}$-order polynomial $\phi$, we can write,

[TABLE]

Setting $\phi(\tau)=B_{j}(\tau)$ in (39) and evaluating the resulting equation at $\tau_{i}$, we get,

[TABLE]

We can rewrite (40) for Case (b) evaluation as,

[TABLE]

Using (36) to evaluate (41) we get,

[TABLE]

From (38), it is clear that (42) is the same as ${\boldsymbol{D}}_{\omega}{\boldsymbol{B}}_{\omega}={\boldsymbol{I}}_{N}$ for $\omega=b$. It is also apparent that the proof for $\omega=a$ follows by a similar procedure.

4 A Family of Well-Conditioned PS Optimal Control Methods

In computational optimal control, it is apparent that it is important to generate a reasonably accurate control trajectory $[t_{0},t_{f}]\ni t\mapsto{\boldsymbol{u}}$. Furthermore, it is more important to generate an accurate value of the initial control than its final value because any error in ${\boldsymbol{u}}(t_{0})$ will be amplified in terms of propagating the state trajectory through the dynamics $\dot{\boldsymbol{x}}={\boldsymbol{f}}(\boldsymbol{x},{\boldsymbol{u}}(t))$; see [30, 39] for details. This implies that the constraint at the initial point $\dot{\boldsymbol{x}}(t_{0})={\boldsymbol{f}}(\boldsymbol{x}(t_{0}),{\boldsymbol{u}}(t_{0}))$ must be well-represented. In [36], initial value problems (IVPs) are solved using the matrix ${\boldsymbol{B}}_{a}$ and the differential equation is not imposed at the initial point. Because this approach is not satisfactory for optimal control applications, we need to modify existing results on Birkhoff interpolation for optimal control applications.

4.1 Birkhoff PS Optimal Control Method Over an Arbitrary Grid Based on (31)

In applying the results of the previous section to Problem $P$, we set,

[TABLE]

where, $x_{0}$ and $v_{j},j=1,\ldots,N$ are the unknown optimization variables over an arbitrary grid $\pi^{N}$, and $B_{k},k=0,1,\ldots,N$ satisfy (33). Differentiating both sides of (43) we get,

[TABLE]

Substituting (44) in the dynamic constraint $\dot{x}^{N}(\tau)=f(x^{N}(\tau),u^{N}(\tau))$, we get,

[TABLE]

Evaluating (45) over the grid $\pi^{N}=[\tau_{0},\pi^{N}_{a}]$ we get,

[TABLE]

where, $I_{a},V_{a}$ and $X_{a}$ are given by,

[TABLE]

The vector $X_{a}$ can be evaluated using (43); this generates the constraint,

[TABLE]

where,

[TABLE]

Collecting all the relevant equations, we define

[TABLE]

Remark 4.3**.**

The Birkhoff equality constraint $I_{a}\,V_{a}+x_{0}\dot{B}_{0}(\tau_{0})-f(x_{0},u_{0})=0$ in Problem $P^{N}_{a}$ imposes the differential equation $\dot{x}(\tau)=f(x(\tau),u(\tau))$ at the initial point $\tau_{0}$. This Birkhoff constraint is equivalent to imposing the same differential equation via the Lagrange condition $\sum_{j=0}^{N}x_{j}\dot{L}_{j}(\tau_{0})-f(x^{N}(\tau_{0}),u^{N}(\tau_{0}))=0$. This equivalency is proved as Proposition 4.4.

Proposition 4.4**.**

Let $x_{j}=x^{N}(\tau_{j})$. At $\tau=\tau_{0}$, the Lagrange and Birkhoff interpolants satisfy the condition,

[TABLE]

Proof 4.5**.**

The derivative $\dot{B}_{j}(\tau_{i})$ can be obtained by an application of the PS differentiation matrix; hence, we have,

[TABLE]

Setting $i=j=0$ in (58) we get,

[TABLE]

where, the last equality in (59) follows from the fact that $B_{0}(\tau_{0})=1$; see (33).

Similarly, setting $i=0$ in (58) we get,

[TABLE]

Evaluating the right-hand-side of (57) using (59) and (60) we have,

[TABLE]

Next, from (43), we have,

[TABLE]

Evaluating the left-hand-side of (57) using (62) we get,

[TABLE]

It is evident that the right-hand-sides of (61) and (63) are identical.

4.2 Birkhoff PS Optimal Control Method Over an Arbitrary Grid Based on (34)

Although many of the results are similar to those resulting from (43), we provide a brief development here to highlight the key differences between the two approaches. To this end, we set,

[TABLE]

Substituting the (64) in the dynamic constraint we get,

[TABLE]

Evaluating (65) over the grid $\pi^{N}=[\pi^{N}_{b},\tau_{N}]$ we get,

[TABLE]

where, $I_{b},V_{b}$ and $X_{b}$ are given by,

[TABLE]

Using (64) to evaluate $X_{b}$ we get,

[TABLE]

where,

[TABLE]

Thus, we now have a new discretization of Problem $P$ given by,

[TABLE]

Remark 4.6**.**

In Problem $P^{N}_{b}$, the differential constraint $\dot{x}(\tau)-f(x(\tau),u(\tau))=0$ at the initial point $\tau_{0}$ is imposed via the constraint $V_{b}-f(X_{b},U_{b})=0$. The Birkhoff equality constraint $I_{b}\,V_{b}+x_{N}\dot{B}_{N}(\tau_{N})-f(x_{N},u_{N})=0$ imposes the differential constraint at the final time. This Birkhoff condition is equivalent to imposing the same differential constraint via the Lagrange condition $\sum_{j=0}^{N}x_{j}\dot{L}_{j}(\tau_{N})-f(x^{N}(\tau_{N}),u^{N}(\tau_{N}))=0$. The proof of this equivalency is similar to Proposition 1; hence, it is omitted.

4.3 Several Pre-Conditioned Lagrange PS Optimal Control Methods

It is possible to develop several pre-conditioned Lagrange PS optimal control methods based on using Birkhoff matrices as pre-conditioners. We first consider the case where we partition the mesh $\pi^{N}$ as $[\tau_{0},\pi^{N}_{a}]$.

4.3.1 Left-Preconditioned Lagrange PS Method Based on $[\tau_{0},\pi^{N}_{a}]$ Partitioning

It is convenient to rewrite (12) as,

[TABLE]

Differentiating (77) and substituting the result in the dynamic constraint $\dot{x}=f(x,u)$, we get,

[TABLE]

Evaluating (78) over the grid $\pi^{N}=[\tau_{0},\pi^{N}_{a}]$ we get,

[TABLE]

where, $X_{a}$ and $U_{a}$ are defined in (48) and ${\boldsymbol{l}}_{0}$ is defined by,

[TABLE]

Note that (79) is merely a partitioned version of the discrete differential constraint stipulated in Problem $P^{N}$ defined in (22). Pre-multiplying both sides of (79b) by ${\boldsymbol{B}}_{a}$ and using Theorem 1 we get,

[TABLE]

As a result of (81), Problem $P^{N}$ transforms according to,

[TABLE]

4.3.2 Right-Preconditioned Lagrange PS Method Based on $[\tau_{0},\pi^{N}_{a}]$ Partitioning

Instead of pre-multiplying equations to achieve better conditioning, we can generate an alternative method by post-multiplying only the left-hand-side of (79b) by $B_{a}$. To clarify this point, we substitute (49) in (79b); this results in,

[TABLE]

Comparing (87) with (47), it is evident that we need to prove that $x_{0}{\boldsymbol{D}}_{a}{\boldsymbol{b}}_{0}=-x_{0}{\boldsymbol{l}}_{0}$. This result is proved as Proposition 4.7.

Proposition 4.7**.**

${\boldsymbol{D}}_{a}{\boldsymbol{b}}_{0}+{\boldsymbol{l}}_{0}={\bf 0}$**

Proof 4.8**.**

From (58), we have,

[TABLE]

Evaluating (88) for $j=0$ and $i=1,\ldots N$, we get,

[TABLE]

The left-hand-side of (89) is zero from (33); hence, the proposition follows from (80) and the fact that $B_{0}(\tau_{0})=1$ by definition.

As a result of Proposition 4.7, (87) simplifies to $V_{a}=f(X_{a},U_{a})$. Using Proposition 4.4 to replace (79a), Problem $P^{N}$ now transforms to a “right-preconditioned” method given by,

[TABLE]

It is apparent that (95) is identical to Problem $P^{N}_{a}$ defined by (56). In other words, the Birkhoff-PS method (for Case (a)) may be viewed as a right-pre-conditioned Lagrange-PS method.

Proposition 4.9**.**

If the variable $V_{a}$ is eliminated from Problem right-$P^{N}_{a}$, it reduces to Problem left-$P^{N}_{a}$.

Proof 4.10**.**

Comparing (86) with (95), it follows that this proposition is proved if we can show that:

• (1)

${\boldsymbol{B}}_{a}f(X_{a},U_{a})-x_{0}{\boldsymbol{B}}_{a}{\boldsymbol{l}}_{0}=x_{0}\,{\boldsymbol{b}}_{0}+{\boldsymbol{B}}_{a}V_{a}$, and

• (2)

the equation $\sum_{i=0}^{N}x_{i}\dot{L}_{i}(\tau_{0})=f(x_{0},u_{0})$ is equivalent to $I_{a}\,V_{a}=f(x_{0},u_{0})-x_{0}\dot{B}_{0}(\tau_{0})$

Hence, we prove this proposition in two parts.

Part (1):**

Substituting (47) in (49), we get,

[TABLE]

where, the last equality in (96) follows from Proposition (4.7) and Theorem 1.

Part (2):**

The proof of this part is a direct consequence of Proposition 4.4.

Remark 4.11**.**

It is clear that two additional left and right pre-conditioned Lagrange-PS methods can be similarly obtained by using the (Case (b)) partitioning $\pi^{N}=[\pi^{N}_{b},\tau_{N}]$. For the purposes of brevity, these cases are not derived.

Remark 4.12**.**

For Case (a), it is straightforward to show that $B_{0}(\tau)=1$ for all $\tau\in[\tau_{0},\tau_{N}]$. Consequently, ${\boldsymbol{b}}_{0}$ defined in (50) is a vector of ones and $\dot{B}_{0}(\tau_{0})=0$. Similarly, for Case (b), $B_{N}(\tau)=1$ for all $\tau\in[\tau_{0},\tau_{N}]$; hence, ${\boldsymbol{b}}_{N}$ defined in (70) is a vector of ones and $\dot{B}_{N}(\tau_{N})=0$.

Remark 4.13**.**

All PS methods proposed in this section are for an arbitrary grid. It is possible to generate a very large sub-class of each of these methods by choosing slightly different and well-known grid points. Popular examples of such grid points are Chebyshev-Gauss-Lobatto (CGL), Legendre-Gauss-Lobatto (LGL), Chebyshev-Gauss-Radau (CGR), Legendre-Gauss-Radau (LGR), Chebyshev-Gauss (CG) and Legendre-Gauss (LG). See [1, 42, 45] for further details on advantages and pitfalls associated with the choice of these mesh points as they apply to computational optimal control problems.

Because it is possible to generate a very large family of PS methods based on the results of the preceding subsections, we limit the scope of the present paper to illustrating only a small set of these methods. A more extensive discussion of other options and numerical results are discussed in [54]. Nonetheless, we note that a selection of the proper grid for optimal control problems is largely based on whether or not the horizon is finite or infinite[39].

5 A Numerical Demonstration of the Well-Conditioning of Birkhoff PS Discretizations

A quick examination of (95) shows that the linear equation given by (49) is the only equation in a Birkhoff PS method (for Case (a)) that is independent of the problem data functions $E,e$ and $f$. Consequently, if the original continuous-time problem is well-conditioned, the conditioning of the discretized problem is driven by the condition numbers of the linear equation given by (49) for Case (a), and by (69) for Case (b). Because all conclusions drawn with respect to the condition number of Case (a) are identical to that of Case (b), we focus the entirety of our discussions in this section to that of Case (a) only.

A very powerful statement on the condition number of the Birkhoff PS discretization can be drawn for the special case of a known value of the initial condition, $x_{0}$. In this case, we can re-write (49) as,

[TABLE]

where, the right-hand-side of (97) is a known quantity. Hence, the condition number of a Birkhoff PS discretization is driven by the condition number of the ($N\times 2N$)-matrix ${\boldsymbol{C}}^{\text{Birk}}_{N}$ defined by,

[TABLE]

where ${\boldsymbol{I}}_{N}$ is the $N\times N$ identity matrix. As evident from Fig. 6,

the condition number of ${\boldsymbol{C}}^{\text{Birk}}_{N}$ is remarkably flat with respect to $N$; i.e., of $\mathcal{O}(1)$. The condition number of a corresponding Lagrange PS discretization may be obtained by multiplying both sides of (97) by ${\boldsymbol{D}}_{a}$ and invoking Theorem 1. This results in

[TABLE]

with the right-hand-side of (99) considered to be a known quantity. Hence, the condition number of a Lagrange PS discretization is driven by the condition number of the ($N\times 2N$)-matrix,

[TABLE]

It is clear from Fig. 7

that the condition number of ${\boldsymbol{C}}^{\text{Lagr}}_{N}$ is $\mathcal{O}(N^{2})$ and consistent with the plot shown in Fig. 2.

In a generic optimal control problem, $x_{0}$ is an optimization variable. In this general case, the right-hand-side of (97) is not known; hence, we need to examine the condition number of the $N\times(2N+1)$-matrix,

[TABLE]

It is apparent from Fig. 8 that the condition number of ${\boldsymbol{A}}_{N}^{\text{Birk}}$ varies only as $\mathcal{O}\left(\sqrt{N}\right)$.

The condition numbers for other choices of grid selections (see Remark 4.13) have similar behavior; see, for example, Fig. 2. In this paper we largely focus on Lobatto-based grid selections because they are the correct choice for solving generic finite-horizon optimal control problems (see for example, Ref. [1], Section 4.4.4 in Ross[39], and Ref. [42]). Radau-based mesh points — first introduced in [55] and [56] for optimal control applications — are more appropriate for infinite-horizon optimal control problems because the non-unit weight function associated with the Legendre-Radau-grid vanishes at the $+1$ (or $-1$) point[1, 56]. Consequently, a Legendre-Radau grid is mathematically justified only for some limited choices of boundary conditions in a finite-horizon optimal control problem[1, 42]. Using a Radau mesh for finite-horizon optimal control applications with arbitrary boundary conditions produces well-known convergence issues and control-chatter problems[1, 39, 42]. Pure Gauss points are even more limiting with regards to their suitability for solving generic optimal control problems[1, 42].

6 An Illustrative Example

There are many emerging problems in aerospace trajectory optimization that require a fine mesh[31, 32, 33, 57]. A detailed discussion of these problems is beyond the scope of this paper; hence, we use a proxy low-thrust orbital maneuvering problem[47, 57] as a prototypical challenge to illustrate the numerics of the new PS methods developed in the preceding sections. That is, we choose this problem for illustrative purposes only and not necessarily to address the myriad of issues that are associated with the specifics of low-thrust orbital maneuvering. To this end, consider a minimum-time circle-to-circle continuous thrust orbit transfer problem (see Fig. 9) given by,

[TABLE]

where, $r$ is the radial position of the spacecraft, $\theta$ is the true anomaly, $v_{r}$ and $v_{t}$ are the radial and transverse components of the spacecraft velocity, $\alpha$ is the thrust steering angle, $\mu$ is the gravitational parameter and $A$ is the magnitude of the constant thrust acceleration. What makes Problem $O_{\text{Xfer}}$ challenging is that if $A$ is low and $\left(r^{f}/r^{0}\right)$ is high, then the number of orbital revolutions to attain the final orbit is high; hence, from an elementary application of the Nyquist-Shannon sampling theorem, it is apparent that the number of nodes must be high. In principle, solving a problem with a large number of nodes is not a significant technological hurdle[31, 39]. As an elementary illustration of the technology requirements, we note that a million variables only requires 8 MB of computer memory. What makes a high-node problem hard to solve is science and not technology; that is, the problem is the high condition number of the differentiation matrix associated with a PS method. For example, from Fig. 2, the condition number of a PS differentiation matrix for $N=1000$ is about $10^{6}$; i.e., $\mathcal{O}(N^{2})$. Because it is equivalent to a right-preconditioned Lagrange PS method (see (95)), a Birkhoff PS method mitigates the problem of high condition numbers. For example, for $N=1000$, we expect the condition number to be about $30$; see Fig. 8. To illustrate these features, we solve Problem $O_{\text{Xfer}}$ using the following canonical units:

[TABLE]

The boundary conditions for a representative LEO-to-GEO orbit transfer problem (in canonical units) is given by,

[TABLE]

A low-thrust solution ($A=5\times 10^{-4}$) resulting from an application of Birkhoff PS method is shown in Fig. 10.

This solution was obtained using a CGL mesh and the right preconditioning method defined by (95). Results for other cases and problems including solutions from the vast variety of other possible Birkhoff PS methods outlined in Remark 4.13 are discussed in [54]. In order to limit the scope of this section, we only discuss this particular case.

The time to accomplish the maneuver shown in Fig. 10 was approximately 1187 time units; this translates to approximately 13 days of continuous low-thrust maneuvering. The number of CGL points was $N=1024$. In other words, $r(t)$ shown in Fig. 10 is given by a polynomial of over a thousand degree! The “discrete” steering angle at 1024 CGL points is shown in Fig. 11.

When this steering angle is interpolated through the CGL mesh, we get the continuous function $t\mapsto\alpha$ indicated by the solid line in Fig. 11. Using the interpolated function $t\mapsto\alpha$, the dynamical equations can be easily integrated using ode45 (from MATLAB) to generate a propagated trajectory. This propagated trajectory is included in Fig. 10; however, this plot is not noticeable in this figure because the differences are too small to be visually apparent.

Consequently, we provide an error plot in Fig. 12 as a matter of completeness.

7 Conclusions

State trajectories in aerospace dynamical systems are absolutely continuous; hence, by an elementary application of the Stone-Weierstrass theorem, they can be expressed to any arbitrary precision in terms of some finite-order polynomial. In using the big two polynomials, namely Legendre and Chebyshev, pseudospectral (PS) optimal control theory has found numerous applications that include include flight implementations in aerospace engineering. In many practical PS optimal control applications, very accurate solutions can be generated using polynomials of order $N\leq 100$, thanks to the spectral (i.e., near exponential) rate of convergence. Emerging applications have sought to generate solutions in excess of $N>1000$. The challenge with such high orders is the high condition number ($\sim\mathcal{O}\left(N^{2}\right)$ ) of the associated differentiation matrix. This is largely a science problem rather than a technological one. That is, a computer’s processing capability is not a limiting factor because 8MB of memory for a million variables is no longer a technological hurdle. New advancements in the science associated with Birkhoff interpolation provide a dramatic drop in the condition number: from $\mathcal{O}\left(N^{2}\right)$ to \mathcal{O}\big{(}\sqrt{N}\big{)} or even \mathcal{O}\big{(}1\big{)} in some cases. In advancing these new ideas, the problem with condition numbers is now relegated to specific applications and not to PS optimal control techniques. Nonetheless, a significant amount of new research on Birkhoff PS methods remains to be done.

Acknowledgment

We take this opportunity to thank the anonymous reviewers for providing many detailed comments. Their comments greatly improved the quality of this paper.