    # Reachability as a Unifying Framework for Computing Helicopter Safe Operating Conditions and Autonomous Emergency Landing

We present a numeric method to compute the safe operating flight conditions for a helicopter such that we can ensure a safe landing in the event of a partial or total engine failure. The unsafe operating region is the complement of the backwards reachable tube, which can be found as the sub-zero level set of the viscosity solution of a Hamilton-Jacobi (HJ) equation. Traditionally, numerical methods used to solve the HJ equation rely on a discrete grid of the solution space and exhibit exponential scaling with dimension, which is problematic for the high-fidelity dynamics models required for accurate helicopter modeling. We avoid the use of spatial grids by formulating a trajectory optimization problem, where the solution at each initial condition can be computed in a computationally efficient manner. The proposed method is shown to compute an autonomous landing trajectory from any operating condition, even in non-cruise flight conditions.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

This paper presents a method to efficiently compute unsafe operating flight conditions for which no control sequence exists to safely land the helicopter in the event of partial or total engine failure. A chart called a height-velocity diagram or H-V diagram needs to be produced for each airframe, and a human pilot can utilize the diagram to avoid operating in unsafe conditions. If operating in the safe region, there exists a control sequence to initiate an autorotation, whereby the helicopter enters a glide slope so as to maintain rotor inertia and then flare to slow down near the ground and land safely (FAA, 2019, Chapter 11). H-V diagrams are ultimately determined through flight testing, which is inherently dangerous since the test objective is to define the unsafe boundaries of flight operations. Therefore, there is a need to accurately compute the safe and unsafe regions directly from the helicopter dynamics, and construct a H-V diagram prior to flight testing. The capability to produce H-V diagrams is desired for aircraft that are in the early design stages and have not yet flown, aircraft that are already flying and preparing for an H-V flight test event, as well as refining the H-V diagram of operational aircraft. The capability is required for aircraft that do not have a representative, high-fidelity flight simulation, as well as for more mature aircraft for which a high-fidelity simulation is available.

We consider the safe region as any initial flight condition where there exists a control sequence that can steer the system to a safe landing condition, i.e. minimal vertical and horizontal velocity at the ground level, rotor near level, etc. Determining the set of states of a dynamical system that can be driven into a particular final condition is commonly referred to as reachability analysis, and reachable sets can be determined from the sub-zero level set of the viscosity solution to a Hamilton–Jacobi (HJ) partial differential equation (PDE) (

Mitchell et al. (2005)).

Traditionally, these HJ PDEs are solved numerically by constructing a dense discrete grid of the solution space (Mitchell (2008)), and are supported by mature theory (Osher and Fedkiw (2006)). Despite this, HJ reachability analysis has suffered one critical draw-back: Computing the elements of a spatial grid scales poorly with dimension, and therefore have limited applicability for vehicle problems where the dimension of the state space is greater than four. The work of Harno and Kim (2018) attempted to compute helicopter safe operating regions using the grid-based method of Mitchell (2008), but was restricted to a two-dimensional state space, and therefore is not representative of actual helicopter motion and is not applicable for safety critical use.

Recent research based on generalizations of the Hopf formula (Hopf (1965)) have avoided a spatial grid of the state space and instead form a grid over time, where numerical solutions are obtained by constructing a trajectory optimization problem. These techniques were used in Darbon and Osher (2016) to compute solutions to the HJ equation for systems with state and time-independent Hamiltonians of the form , and Kirchner et al. (2018a) expanded the classes of systems to general linear systems in high-dimensions. Additionally, these methods were successfully applied to vehicle control problems such as collaborative pursuit-evasion (Kirchner et al. (2018b)) and coordination of heterogeneous groups of vehicles (Kirchner et al. (2020)).

The methods described above based on the Hopf formula can be seen as an optimization problem based on the costate trajectory of the system. We propose in this work to formulate an optimization of the state trajectory directly, which allows consideration of the non-linear dynamics encountered in helicopter motion. There is no need to directly formulate the optimal Hamiltonian of the system, since the Hamiltonian frequently resulting from helicopter dynamics does not have a standard form.

There have been various proposals to compute the H-V diagram using trajectory optimization (see e.g. Johnson (1977); Lee et al. (1988); Carlson et al. (2006); Yomchinda (2013); Bibik and Narkiewicz (2012)), but they suffer several drawbacks. These methods minimize cost functionals with weighted cost terms and a terminal state constraint, which do not guarantee necessary ground landing conditions are met. Therefore, the solution of these optimization problems do not, in general, give the reachable set. Since the functional is parameterized by the terminal state at precisely terminal time, , time is an unknown free parameter that must be solved for. When combined with the terminal constraint, it is not guaranteed that a candidate trajectory is feasible for a particular time and numeric convergence issues can be observed.

We propose to construct an implicit surface representation of the desired safe landing condition and, in doing so, create a trajectory optimization that gives the backwards reachable tube (Bansal et al. (2017)). We consider the case where the time-to-land need not be known a priori by considering that the intersection of the helicopter with the ground could occur at any time on the interval . We utilize Lagrange polynomials to create a collocation scheme that converges rapidly with the number of time samples and can represent the required half-infinite time intervals. From this we compute the optimal landing sequence to autonomously achieve a safe landing from any initial condition that is within the backwards reachable tube.

## 2 Reachability and Safe Operating Regions

We consider helicopter dynamics

 ddsx(s)=f(x(s),u(s)), (1)

for where is the system state and is the control input, constrained to lie in the admissible control set . We denote by the state trajectory that evolves in time with measurable control input according to starting from initial state at . The trajectory is a solution of in that it satisfies almost everywhere:

 ddsγ(s;x,u(⋅)) (2) γ(0;x,u(⋅)) =x.

We denote by as the goal set that represents the set of acceptable safe landing conditions. We seek to determine if a control sequence exists that drives the system into at exactly time . The set of all initial states where there exists a control to drive the system to the set at exactly time is called the backward reachable set (BRS) of the system and is defined as

 ΘS(t)={x:∃u(⋅)∈U,γ(t;x,u(⋅))∈Θ0}. (3)

### 2.1 The Connection of the BRS to the HJ Equation

We represent the set of safe landing conditions, , implicitly with the function such that

 Θ0={x∈Rn|J(x)≤0} (4)

and use it to construct a cost functional for the system trajectory , given terminal time as

 K(t,x,u(⋅))=∫t0IU(u(s))ds+J(γ(t;x,u(⋅))), (5)

where the running cost function

is the characteristic function for the set

and is defined by

The value function is defined as the minimum cost, , among all admissible controls for a given state as

 φ(x,t)=inf% u(⋅)∈UK(t,x,u(⋅)). (6)

The value function in satisfies the dynamic programming principle (Evans (2010)) and also satisfies the following initial value Hamilton–Jacobi (HJ) equation with being the viscosity solution of

 {∂φ∂s(x,s)+H(x,∇xφ(x,s))=0,φ(x,0)=J(x), (7)

for , where the Hamiltonian is defined by

 H(x,p)=supu∈U⟨−f(x,u),p⟩, (8)

where we denote by the costate variable. [Mitchell (2007a)] The zero sub level sets of the viscosity solution is an implicit surface representation of the finite time backwards reachable set, defined in . We note here the important distinction that the BRS in defines the set of initial states that can be driven into the set at precisely time . It is possible for the system to be driven into the set at an earlier time and then later exit the set . This motivated the characterization of the backwards reachable tube in Mitchell et al. (2005).

### 2.2 The Backwards Reachable Tube

We are instead interested if the system can be driven into at any time on the range . This set of initial states, coined the backwards reachable tube (BRT) in Bansal et al. (2017), is defined by

 ΘT(t)={x:∃u(⋅)∈U,∃s∈[0,t],γ(s;x,u(⋅))∈Θ0}. (9)

The seminal work of Mitchell et al. (2005) showed that the backwards reachable tube is found from the zero sub level set of , which is the viscosity solution to the following modified Hamilton–Jacobi equation

 {∂φ∂s(x,s)+max(0,H(x,∇xφ(x,s)))=0,φ(x,0)=J(x), (10)

with the same as given above in . The intuition is that the value function cannot decrease111In Mitchell et al. (2005), a minimum operator is used since time is defined in that work on the interval ., which does not allow the level sets to contract. It can be seen that any trajectory that enters the set is not allowed to escape by “freezing” it in time. The reachable set and reachable tube are connected through the following relation (Mitchell, 2007b, Proposition 1)

 ΘT(t)=⋃s∈[0,t]ΘS(s). (11)

The boolean set operation in implies a corresponding property (Osher and Fedkiw (2006)) of the implicit surface representations

 φT(x,t)=mins∈[0,t]φS(x,s). (12)

The backwards reachable tube, , characterizes the controllably safe regions of flight. Therefore, the set of unsafe states is found from

 ΘC:=limt→∞Rn∖ΘT(t).

## 3 Computing the Reachable Tube Through Optimization

Under a mild set of regularity conditions (Subbotin, 1995, Ch. 7, p. 63), is the unique viscosity solution of (Subbotin, 1995, Th. 8.1, p. 70). The uniqueness of the solution, , implies that the viscosity solution is equivalent to the value function and can be found by minimizing with constraints given by . The HJ of the BRT in can be constructed by augmenting the dynamics of with

 ddsx(s)=g(x(s),u(s),a(s)):=a(s)f(x(s),u(s)), (13)

where is measurable scalar function (Mitchell et al. (2005)). We denote by the state trajectory of the augmented system that satisfies

 =g(ζ(s;x,u(⋅),a(⋅)),u(s),a(s)), ζ(0;x,u(⋅),a(⋅)) =x. (14)

Note that since is constrained to , a value of makes the augmented dynamics of equivalent to the original dynamics of and if the dynamics stop completely. We denote by

 σ(s):=∫s0a(ℓ)dℓ,

the pseudotime variable, and denote by as a quasi-inverse of in the sense that

 σ†(σ(s))=s.

The formal definition of this quasi-inverse is given in (Mitchell et al., 2005, Lemma 6), and it was shown that the trajectories of have the following relations

 γ(σ(s);x,u(σ†(⋅)))=ζ(s;x,u(⋅),a(⋅)),

for any (Mitchell et al., 2005, Lemma 4), and consequently visits only a subset states of the (Mitchell et al., 2005, Corollary 5). Therefore, we can solve the modified optimization problem to find the reachable tube:

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩minu(⋅),a(⋅)J(ζ(t;x,u(⋅),a(⋅)))Subject toddsζ(s)=g(ζ(s),u(s),a(s))ζ(0)=xu(⋅)∈Ua(⋅)∈[0,1]˙a(⋅)≤0, (15)

where the last line of is an optional regularization term since the optimal is not unique.

### 3.1 Time-Domain Transformation

We use collocation to approximate the trajectory and first perform a time domain transformation which maps into the computation interval with the invertible transform defined as

 s=Γ−1(τ)=clog(21−τ), (16)

where the choice of is a scaled version of that proposed in Garg et al. (2011). The dynamics are similarly transformed with

 ddτx(τ) =ddτΓ(τ)g(x(τ),u(τ),a(τ)),

and with the transform defined in becomes

 ddτx(τ)=c(1−τ)g(x(τ),u(τ),a(τ)).

Hereafter, we assume that the trajectory, , is a solution to the transformed dynamics satisfying

 ddτζ(τ;x,u(⋅),a(⋅)) =c(1−τ) ×g(ζ(τ;x,u(⋅),a(⋅)),u(τ),a(τ)) ζ(−1;x,u(⋅),a(⋅)) =x.

### 3.2 Finite Trajectory Approximation

We propose a polynomial trajectory approximation, which is commonly referred to as pseudospectral optimal control, and was introduced in Elnagar et al. (1995) and later refined in Ross and Karpenko (2012) and Garg et al. (2010). We denote by as the approximation to with Lagrange polynomials as

 ζ(τ)≈ζN(τ)=N∑j=0xjLj(τ), (17)

defined by a set of points, , sampled on the time grid

 πN={τj:j=0,…,N},

where each . The Lagrange basis functions are given as

 Lj(τ):=N∏k=0k≠jτ−τkτj−τk,

and it follows that

 ddτζN(τi)=N∑j=0xj˙Lj(τi)=N∑j=0Dijxj,

where we denote by as the Gauss differentiation matrix constructed with each element, with row and column , given by

 Dij:=˙Lj(τi). (18)

We select as Legendre-Gauss-Radau (LGR) points with and the corresponding differentiation matrix, , and quadrature weights, , are found from (Shen et al., 2011, Section 3.3). Note that since LGR points do not include a point at , we avoid a singularity in at . We denote by as the first column of corresponding to the boundary condition and denote by as the remaining columns representing the interior nodes such that

 D=[D0DI]. (19)

We denote by

as the concatenated vector of all collocation points for

given as

 X:=(x1,⋯,xN)⊤, (20)

and likewise denote as the concatenated vector of all the collocated control input points

 U:=(u0,⋯,uN)⊤,

where each is the control input at each time . And similarly, for the augmented inputs, , we have

 A:=(a0,⋯,aN)⊤.

Recall that . We denote by as the concatenated vector of the equations of motion,

 G(X,U,A):=⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣a0c(1−τ0)f(x,u0)a1c(1−τ1)f(x1,u1)⋮aNc(1−τN)f(xN,uN)⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦,

evaluated at each point in . It follows the differentiation matrix of the concatenated of state is

 DI=DI⊗In,

where denotes the Kronecker product and is the identity matrix. Likewise, we have

 D0=D0⊗In.

The dynamic equality constraint is now written as

 DsX+D0x=G(X,U,A).

Recall that is the vector of LGR quadrature weights. It was proposed in Garg et al. (2010)

to estimate the terminal state with

 xN+1=x+(w⊤⊗In)G(X,U,A).

With

as the given initial condition, we construct the non-linear programming problem (NLP):

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩minX,U,AJ(xN+1)Subject toDsX+D0x=G(X,U,A)xN+1=x+(w⊤⊗In)G(X,U,A)uj∈Uaj∈[0,1]aj−aj−1≤0, (21)

which has a numerical solution that approximates when is sufficiently large. The time of minimum cost in can be found from the following quadrature:

 t∗=w⊤⎡⎢ ⎢ ⎢ ⎢⎣a0c(1−τ0)⋮aNc(1−τN)⎤⎥ ⎥ ⎥ ⎥⎦. (22)

## 4 Helicopter Dynamics

We define the state vector as with the vertical velocity, ; horizontal velocity, ; rotor angular speed, ; height above ground level, ; horizontal displacement, ; and engine power, . The control inputs are given as = with thrust coefficient, , effectively the collective control, and , the aircraft pitch angle. Following Lee (1985), the dynamics are as follows:

 f:⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩˙y=g−ρm(πR2)(ΩR)2CTcosθ+12ρmfehVfy˙v=ρm(πR2)(ΩR)2CTsinθ−12ρmfezVfv˙Ω=1IRΩ(P−1ηρ(πR2)(ΩR)3Cp)˙h=−y˙z=v˙P=−1κP (23)

where denotes the engine response time constant, the rotor power efficiency factor, the mass of the helicopter, the velocity magnitude of the fuselage, the rotor solidity, the rotor radius, the rotor inertia, the density of the air, and the acceleration due to gravity. and is the flat plat drag area in the vertical and horizontal directions, respectively. denotes the power coefficient defined as

 Cp=18σRcd+CTλ,

where is the drag coefficient of the rotor airfoil and

 λ=vsinθ−ycosθ+νΩR

is the inflow ratio. The advance velocity, , is given by

 UT=vcosθ+ysinθ.

The variable is induced velocity and we use the inflow model

 ν=ΩR√12CT,

which is the ideal inflow in hover. More sophisticated inflow models exist but are outside the scope of this work, and the reader is encouraged to read Johnson (1977) and Chen and Hindson (1986) for more details.

## 5 Results

We ensure safe landing when the conditions and are met as the helicopter is sufficiently close to the ground; in this case . The parameters and are found from the structural specifications of the airframe. We select the function that satisfies the terminal conditions with

 J(x)=max(|y|−ymax,|v|−vmax,|h|−1),

where, for this example, we chose and . The remaining initial states are chosen as the trim conditions such that , and initial engine power was set , simulating an instant engine failure and providing a safety factor as it gives the largest unsafe set. The input is bounded by a blade stall condition with (Lee et al. (1988)), where . The rest of the coefficients for the model in are from (Yomchinda, 2013, Table 2-1, Table A-1). The number of interior time samples, , is fixed at and the constant in is set to . To take advantage of the inherent sparsity of the formulated optimization, we use the NLP code IPOPT (Wächter and Biegler (2006)), where the constraint Jacobian was computed using automatic differentiation by the methods of Andersson et al. (2019).

The value function for the system was computed using the method of Section 3, and an H-V diagram is produced by finding the unsafe regions, which is all areas where the value function is greater than zero. Figure 1 shows the H-V diagram for the dynamics given in , where dark regions on the diagram represent initial height above ground level (AGL) and forward speed such that a safe landing in the event of an engine failure is impossible. Note that a disjoint lobe of the unsafe set appears on the bottom right of Figure 1 and is known as the high-speed unsafe region. Figure 1: Shown in green is the set of unsafe H-V flight conditions as computed from the method of Section 3.

### 5.1 Autonomous Autorotation

The arguments that minimize , denoted as and , are the optimal landing trajectory and control sequence, respectively, provided that the helicopter is operating in the reachable tube. As an example we consider a helicopter in a trimmed cruise at with a forward speed of . Figure 2 shows the state evolution as computed at the time of engine failure. The computed trajectory results in a safe landing in as found from . We note the trajectory of the rotor speed, , in Figure (c)c, where an increase in rotor speed is observed due to translational kinetic energy being transfered into rotational kinetic energy to successfully decelerate, or “flare”, at the terminal phase of the landing procedure.

## 6 Conclusion and Future Work

Presented is a numeric method for computing the backwards reachable tube for general non-linear systems. We are motivated by locating the unsafe operating states of a helicopter so they can be avoided by the pilot. By remaining outside of the unsafe set, it ensures it is possible to initiate an autorotation sequence after engine failure and guide the helicopter to a safe landing. The method can additionally be used to compute control inputs in order to complete an autonomous landing via autorotation. Future work includes augmenting with human pilot dynamics as in McRuer and Krendel (1974) and include pilot reaction delay (Kirchner (2019)). Additionally, we will investigate integrating failure detection to initiate the autonomous autorotation landing procedure.

## References

• FAA (2019) (2019). Rotorcraft Flying Handbook, FAA Manual H-8083-21B. Federal Aviation Administration.
• Andersson et al. (2019) Andersson, J.A.E., Gillis, J., Horn, G., Rawlings, J.B., and Diehl, M. (2019). CasADi – A software framework for nonlinear optimization and optimal control. Mathematical Programming Computation, 11(1), 1–36.
• Bansal et al. (2017) Bansal, S., Chen, M., Herbert, S., and Tomlin, C.J. (2017). Hamilton-Jacobi reachability: A brief overview and recent advances. In IEEE 56th Annual Conference on Decision and Control (CDC), 2242–2253. IEEE.
• Bibik and Narkiewicz (2012) Bibik, P. and Narkiewicz, J. (2012). Helicopter optimal control after power failure using comprehensive dynamic model. Journal of Guidance, Control, and Dynamics, 35(4), 1354–1362.
• Carlson et al. (2006) Carlson, E.B., Xue, S., Keane, J., and Kevin, B. (2006). H-1 upgrades height-velocity diagram development through flight test and trajectory optimization. In Annual Forum Proceedings-American Helicopter Society, volume 62, 729. American Helicopter Society.
• Chen and Hindson (1986) Chen, R.T. and Hindson, W.S. (1986). Influence of dynamic inflow on the helicopter vertical response. Technical Report 88327, National Aeronautics and Space Administration.
• Darbon and Osher (2016) Darbon, J. and Osher, S. (2016).

Algorithms for overcoming the curse of dimensionality for certain Hamilton-Jacobi equations arising in control theory and elsewhere.

Research in the Mathematical Sciences, 3(1), 19.
• Elnagar et al. (1995) Elnagar, G., Kazemi, M.A., and Razzaghi, M. (1995). The pseudospectral Legendre method for discretizing optimal control problems. IEEE Transactions on Automatic Control, 40(10), 1793–1796.
• Evans (2010) Evans, L.C. (2010). Partial Differential Equations. American Mathematical Society, Providence, R.I.
• Garg et al. (2011) Garg, D., Hager, W.W., and Rao, A.V. (2011). Pseudospectral methods for solving infinite-horizon optimal control problems. Automatica, 47(4), 829–837.
• Garg et al. (2010) Garg, D., Patterson, M., Hager, W.W., Rao, A.V., Benson, D.A., and Huntington, G.T. (2010). A unified framework for the numerical solution of optimal control problems using pseudospectral methods. Automatica, 46(11), 1843–1851.
• Harno and Kim (2018) Harno, H.G. and Kim, Y. (2018). Safe flight envelope estimation for rotorcraft: A reachability approach. In 18th International Conference on Control, Automation and Systems (ICCAS), 998–1002. IEEE.
• Hopf (1965) Hopf, E. (1965). Generalized solutions of non-linear equations of first order. Journal of Mathematics and Mechanics, 14, 951–973.
• Johnson (1977) Johnson, W. (1977). Helicopter optimal descent and landing after power loss. Technical Report NASA-TM-X-73244, National Aeronautics and Space Administration.
• Kirchner (2019) Kirchner, M.R. (2019). A level set approach to online sensing and trajectory optimization with time delays. In IFAC-PapersOnline, volume 52, 301–306.
• Kirchner et al. (2018a) Kirchner, M.R., Hewer, G., Darbon, J., and Osher, S. (2018a). A primal-dual method for optimal control and trajectory generation in high-dimensional systems. In IEEE Conference on Control Technology and Applications, 1575–1582.
• Kirchner et al. (2018b) Kirchner, M.R., Mar, R., Hewer, G., Darbon, J., Osher, S., and Chow, Y.T. (2018b). Time-optimal collaborative guidance using the generalized Hopf formula. IEEE Control Systems Letters, 2(2), 201–206.
• Kirchner et al. (2020) Kirchner, M.R., Debord, M.J., and Hespanha, J.P. (2020). A Hamilton-Jacobi formulation for optimal coordination of heterogeneous multiple vehicle systems. arXiv preprint arXiv:2003.05792.
• Lee (1985) Lee, A.Y.N. (1985). Optimal Landing of a Helicopter in Autorotation. Ph.D. thesis, Department of Aeronautics and Astronautics. Stanford University.
• Lee et al. (1988) Lee, A.Y., Bryson, A.E., and Hindson, W.S. (1988). Optimal landing of a helicopter in autorotation. Journal of Guidance, Control, and Dynamics, 11(1), 7–12.
• McRuer and Krendel (1974) McRuer, D.T. and Krendel, E.S. (1974). Mathematical models of human pilot behavior. AGARDDograph 188, Advisory Group on Aerospace Research and Development.
• Mitchell et al. (2005) Mitchell, I., Bayen, A.M., and Tomlin, C.J. (2005). A time-dependent Hamilton-Jacobi formulation of reachable sets for continuous dynamic games. IEEE Transactions on Automatic Control, 50(7), 947–957.
• Mitchell (2007a) Mitchell, I.M. (2007a). A toolbox of level set methods. Technical Report TR-2007-11, UBC Department of Computer Science.
• Mitchell (2008) Mitchell, I.M. (2008). The flexible, extensible and efficient toolbox of level set methods. Journal of Scientific Computing, 35(2), 300–329.
• Mitchell (2007b) Mitchell, I.M. (2007b). Comparing forward and backward reachability as tools for safety analysis. In International Workshop on Hybrid Systems: Computation and Control, 428–443. Springer.
• Osher and Fedkiw (2006) Osher, S. and Fedkiw, R. (2006). Level Set Methods and Dynamic Implicit Surfaces, volume 153. Springer Science & Business Media.
• Ross and Karpenko (2012) Ross, I.M. and Karpenko, M. (2012). A review of pseudospectral optimal control: From theory to flight. Annual Reviews in Control, 36(2), 182–197.
• Shen et al. (2011) Shen, J., Tang, T., and Wang, L.L. (2011). Spectral Methods: Algorithms, Analysis and Applications, volume 41. Springer Science & Business Media.
• Subbotin (1995) Subbotin, A.I. (1995). Generalized Solutions of First Order PDEs: The Dynamical Optimization Perspective. Birkhäuser.
• Wächter and Biegler (2006) Wächter, A. and Biegler, L.T. (2006). On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming. Mathematical Programming, 106(1), 25–57.
• Yomchinda (2013) Yomchinda, T. (2013). Real-time Path Planning and Autonomous Control for Helicopter Autorotation. Ph.D. thesis, Department of Aerospace Engineering, Pennsylvania State University.