Fast Mesh Refinement in Pseudospectral Optimal Control

Fast Mesh Refinement in Pseudospectral Optimal Control

Introduction

In principle, mesh refinement in pseudospectral (PS) optimal controlPSReview-ARC-2012 is embarrassingly easy --- simply increase the order $N$ of the Lagrange interpolating polynomial. For instance, in a Chebyshev PS methodboyd,fahroo:cheb-jgcd,cheb-costate, the Chebyshev-Gauss-Lobatto (CGL) mesh points are given byboyd,

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

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 theoryatap. The price for this speed and simplicity is an increase in the condition number of the underlying linear algebra that governs the iterationsPSReview-ARC-2012, boyd, shen-book. In fact, as shown in Fig. (ref),

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

increases as $N^2$. The plots in this figure are the condition numbers of the non-singular, inner, square-portion of the various differentiation matricesPSReview-ARC-2012,boyd 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 segmentfahrooRoss:patching,knots,acc:hybrid,auto-knots. When the segments are non-overlapping and information across the segments is passed across a single point, it results in a PS knotting methodknots,auto-knots; see Fig. (ref).

The special nodes where information across segments is transmitted are called PS knots. In the spectral algorithm for PS methodsauto-knots,spec-alg, 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 regimesauto-knots,spec-alg. The entire process is repeated for the next iterate or until all stopping criteria are metPSReview-ARC-2012,knots,auto-knots,spec-alg.

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 convergenceTAC:linearizable,kang-ijrc,Kang_2008_convergence,kang-rate 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 methodconway-GL,hnw-ode,conway:survey. In fact, the “exponential” rate of convergence of an unknotted PS method is chiefly due to the high orders (i.e., $N\gg10$) of the approximating polynomialsboyd,atap. In many aerospace applicationsKarp-JWST,Kepler-micro-slew,RobStevens-Tether, CDC-Workshop, including flight implementationsPSReview-ARC-2012,zpm:IEEE,bhatt:opm,TRACE-IEEE-Spectrum,TRACE-JGCD-2014,Minelli-AAS,zpm:NASA-report, very accurate solutionsTEI-JGCD-2011,scaling are generated for $N \le 100$. In certain emerging areas of applicationsross-challenge,nolcos16,napa16, 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 problemhesthavan,elbarbary. In recent years, there have been major advancements in well-conditioned PS methodswang,olver,du 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 interpolantswang. 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.

\vskip 0.4in \fbox{

}

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 byross-book, longuski,clarke-2013book, \newlength{\leftside} \newlength{\rightside}

where, the symbols denote the following:

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

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 PSReview-ARC-2012 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,

where, $a_m \in \mathbb R$ is an unknown “spectral” coefficient of an $m^{th}$-order polynomial $P_m$. Note that (ref) is an exact representation and not an approximation. The justification for (ref) comes from the classic Stone-Weierstrass theoremTao:epsilon 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 (ref) to an equivalent “pseudospectral” formross-book, arb-grid-aas, arb-grid,

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,

Equations (ref) and (ref) are known as the modal and nodal representations of $x(\cdot)$ respectivelyross-book.

Let,

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

Note that (ref) 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,

Applying (ref) to (ref), it follows that $x_i = x^N(\tau_i), \ i = 0, \ldots, N$. Differentiating (ref) and evaluating $x^N(t)$ over the grid we get,

Define

as the PS differentiation matrix, and let,

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

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, (ref).

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

where, $C_j(\tau)$ is some interpolating function.

Substituting (ref) and (ref) in $f(x,u)$, and evaluating $f$ over $\pi^N$ we get, $$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)] $$ It is convenient to reuse $f$ as an overloaded function (see ross-book, page 8) defined by, $$ f(X,U) := [f(x_0, u_0), f(x_1, u_1), \ldots, f(x_N, u_N)]^T$$ where,

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

Comparing (ref) and (ref), 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 (ref) is illustrated in Fig. (ref).

This figure was generated using the data and results from arb-grid. The parameter $\beta$ (see Fig. (ref) 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. (ref) 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.

From (ref), 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 (ref).

Problem $P^N$ is solved by the spectral algorithmPSReview-ARC-2012,spec-alg. A key procedure for mesh refinement in the spectral algorithm for PS methods can be described as follows:

Many of the advancements in the spectral algorithm are implemented in DIDO$^\copyright$scaling,ross-book, a state-of-the-art MATLAB$^\circledR$ 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,

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 (ref) 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 knotsknots. 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 auto-knots,spec-alg. 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 polynomialslorentz,schoenberg offers a new PS approach where solutions over thousands of points can be generated in a fast and stable manner.

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 functionlorentz. That is, a Birkhoff interpolant generalizes Lagrange and Hermite interpolantsschoenberg. For second-order boundary value problems (BVPs) with Dirichlet boundary conditions, a Birkhoff interpolation polynomial $\widetilde{x}^N(\tau)$ is given bywang,

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,

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

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 (ref) for first-order systems so that state-space representation may be used. This simple adaptation generates a proposal for a Birkhoff interpolant given by,

where, $B_k, k = 0, \ldots, N$ are polynomials (possibly different from those in (ref)) that must satisfy,

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

The satisfaction of the conditions given by (ref) has some well-known problems related to the existence and uniqueness of a solutionlorentz,schoenberg; hence, we do not pursue this further. If the conditions on $B_0$ or $B_N$ are removed from (ref), then it can be shown that a unique Birkhoff interpolant existsfinden. In doing the latter, a Birkhoff interpolant is given by,

where, $B_k, k = 0, \ldots, N$ must now only satisfy,

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

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

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

The conditions for $B_k, k = 0, \ldots, N$ in (ref) are given by,

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

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

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,

where, ${\boldsymbol B}_a$ and ${\boldsymbol B}_b$ satisfy (ref) and (ref) respectively.

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 scaling,ross-book 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 wang, 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.

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

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

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 (ref). Differentiating both sides of (ref) we get,

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

Evaluating (ref) over the grid $\pi^N = [\tau_0, \pi^N_a]$ we get,

where, $I_a, V_a$ and $X_a$ are given by,

The vector $X_a$ can be evaluated using (ref); this generates the constraint,

where,

Collecting all the relevant equations, we define

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

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

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

Evaluating (ref) over the grid $\pi^N = [\pi^N_b, \tau_N]$ we get,

where, $I_b, V_b$ and $X_b$ are given by,

Using (ref) to evaluate $X_b$ we get,

where,

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

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]$.

Left-Preconditioned Lagrange PS Method Based on $[\tau_0, \pi^N_a]$ Partitioning

It is convenient to rewrite (ref) as,

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

Evaluating (ref) over the grid $\pi^N = [\tau_0, \pi^N_a]$ we get,

where, $X_a$ and $U_a$ are defined in (ref) and ${\boldsymbol l}_0$ is defined by,

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

As a result of (ref), Problem $P^N$ transforms according to,

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 (ref) by $B_a$. To clarify this point, we substitute (ref) in (ref); this results in,

Comparing (ref) with (ref), 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 (ref).

As a result of Proposition (ref), (ref) simplifies to $V_a = f(X_a, U_a)$. Using Proposition (ref) to replace (ref), Problem $P^N$ now transforms to a “right-preconditioned” method given by,

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

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 koeppen. 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 infiniteross-book.

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

A quick examination of (ref) shows that the linear equation given by (ref) 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 (ref) for Case (a), and by (ref) 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 (ref) as,

where, the right-hand-side of (ref) 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,

where ${\boldsymbol I}_N$ is the $N\times N$ identity matrix. As evident from Fig. (ref),

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 (ref) by ${\boldsymbol D}_a$ and invoking Theorem 1. This results in

with the right-hand-side of (ref) 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,

It is clear from Fig. (ref)

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

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

It is apparent from Fig. (ref) 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 (ref)) have similar behavior; see, for example, Fig. (ref). 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. PSReview-ARC-2012, Section 4.4.4 in Rossross-book, and Ref. advances). Radau-based mesh points --- first introduced in Radau-GNC05 and Radau-JGCD 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$) pointPSReview-ARC-2012,Radau-JGCD. Consequently, a Legendre-Radau grid is mathematically justified only for some limited choices of boundary conditions in a finite-horizon optimal control problemPSReview-ARC-2012,advances. Using a Radau mesh for finite-horizon optimal control applications with arbitrary boundary conditions produces well-known convergence issues and control-chatter problemsPSReview-ARC-2012,ross-book,advances. Pure Gauss points are even more limiting with regards to their suitability for solving generic optimal control problemsPSReview-ARC-2012,advances.

An Illustrative Example

There are many emerging problems in aerospace trajectory optimization that require a fine meshross-challenge,nolcos16,napa16,conway:book-chapter. A detailed discussion of these problems is beyond the scope of this paper; hence, we use a proxy low-thrust orbital maneuvering problembellman-low-t, conway:book-chapter 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. (ref)) given by,

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 hurdleross-challenge,ross-book. 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. (ref), 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 (ref)), 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. (ref). To illustrate these features, we solve Problem $O_\text{Xfer}$ using the following canonical units:

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

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

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

The time to accomplish the maneuver shown in Fig. (ref) 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. (ref) is given by a polynomial of over a thousand degree! The “discrete” steering angle at 1024 CGL points is shown in Fig. (ref).

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. (ref). 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. (ref); 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. (ref) as a matter of completeness.

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 \le 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 \big($\sim\mathcal{O}\left(N^2\right) $ \big) 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.