Aircraft longitudinal decoupling based on a singular perturbation approach

Abstract

Aircraft longitudinal dynamics is approximated by short-time mode and phugoid mode from experience. In this article, a rigorous mathematical method is provided based on the singular perturbation theory to deal with this decoupling problem. The longitudinal decoupling and singular perturbation are first introduced. The longitudinal equations are normalized and transformed into a canonical form to extract the perturbation coefficient. Thus, the entire dynamics model is partitioned into boundary-layer equations and slow equations according to singular perturbation theory, which presents a proof to the experience method. The simulation results show that the proposed decoupling approach is sufficiently excellent to approximate the underlying model both in time domain and frequency domain.

Introduction

Aircraft flight mechanics has always been an important issue for more than a century. The problem has a close relationship with flight controller design. Recently, unmanned aerial vehicles (UAVs) are widely spread and the autopilot study has been a hot topic. Thus, the classical flight mechanics draws attentions again. The nonlinear flight dynamic equations are generally linearized at the equilibrium or the operation point. Thus, the linearized equations are divided into lateral and longitudinal directional motion. For the longitudinal motion, it can be furthermore partitioned into short-time and phugoid modes according to the time scales.

There are many researchers who publish papers studying the longitudinal dynamics and the control problem. Argyle and Beard ¹ utilizes the total energy control system to regulate the airspeed dynamics. Liu and Chen ² investigated the aircraft dynamics under wind disturbances and developed a novel nonlinear disturbance observer for the particular UAV model with wind influences. Dlamini and Jones ³ investigated the robustness of the control system by analyzing how characteristic eigenvalues move as a result of damage and comparison to the non-fly-by-wire (FBW) aircraft is made. Alwi and Edwards ⁴ proposed a novel nonlinear fault-tolerant scheme for longitudinal control of an aircraft system, comprising an integral sliding mode control allocation scheme and a backstepping structure. Bronsvoort et al. ⁵ illustrated that inaccurate longitudinal aircraft intent and resulting excessive prediction errors are evident. Stojiljkovic et al. ⁶ presented an application of the root locus technique for the design of a feedback control system of an F-104A aircraft. Ishihara et al. ⁷ presented a method to determine the time-delay margin of a nonlinear adaptive outer-loop control system for a longitudinal aircraft model. Rosa et al. ⁸ described an application of a new fault detection and isolation technique based on set-valued observers to a linear parameter varying longitudinal aircraft dynamic model. Richardson et al. ⁹ designed a gain scheduled controllers for a mathematical second-order longitudinal aircraft model which represents an approximate BAE Hawk wind-tunnel model. Naldi and Marconi ¹⁰ solved numerically both minimum-time and minimum-energy optimal transition problems in order to compute reference maneuvers to be employed by the onboard flight control system to change the current flight condition. Gibbens and Medagoda ¹¹ introduced a method of model predictive control (MPC) using variable prediction time intervals to reduce the level of computation. Jan and Rudolf ¹² focused on the system behavior and the longitudinal dynamics of an aircraft, and an associated feedback controller to augment the damping pitch. Suresh et al. ¹³ proposed a novel neural network approach based on model-following direct adaptive control system design to improve damping. Nicolosi et al. ¹⁴ tested a prototype of this light aircraft in flight for a post-design performance optimization and for the assessment of flight qualities. Ataei and Wang ¹⁵ presented a robust method to design a nonlinear controller for longitudinal dynamics of a hypersonic aircraft model with parametric uncertainties. Oosterom and Babuška ¹⁶ provided the interpolation mechanism for the local flight control law parameters. Karmali and Shelhamer ¹⁷ described a practical and procedure for modeling and identifying the flight dynamics of UAVs. Adaptive control for the decoupling model of aircraft flight mechanics is also an interesting topic.^18–20 The singular perturbation method can also be found in the composite learning control of flexible-link manipulator.^21,22 B Xu et al.^23,24 studied the longitudinal dynamics of the hypersonic flight. These studies made use of the decoupled model, but did not give a reasonable interpretation of the approximation.

This work proposed a rigorous theoretical support to the aircraft longitudinal decoupling. The combination of the singular perturbation theory and the aircraft flight mechanics presents a novel perception of the classical decoupling issue. The simulation depicts that the proposed decoupling method can obtain an approximation of high quality to the original model.

Aircraft longitudinal dynamics

The aircraft flight mechanics is generally partitioned into lateral motion and longitudinal motion. The longitudinal motion denotes the dynamics on the vertical plane as shown in Figure 1.

Figure 1.

Illustration of an aircraft longitudinal motion.

Figure 1 depicts the longitudinal motion and the corresponding state variable. Pitch rate and angle of attack are the state variables of the short-period mode. Flight path angle is the difference of pitch angle and angle of attack. Flight path angle and airspeed are the state variables of the phugoid mode. Thrust and elevators are the control inputs responsible for pitch rotation and velocity regulation, respectively. The longitudinal equations tend to be transformed into linear time-invariant equations as follows

$$\stackrel{·}{\mathbf{X}}={\mathbf{A}}_0\mathbf{X}+{\mathbf{B}}_0\mathbf{U}$$

(1)

which has many specific expressions for different state variable selections. ²⁵ We choose the state vector X^T = [V θ q α]^T and input vector U^T = [δE δT]^T, and hence, the corresponding matrices in the state equations are

$${\mathbf{A}}_0=\left[\begin{array}{cccc}-D_V\hfill & -D_{\alpha }-g\hfill & 0\hfill & -g\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 1\hfill & 0\hfill \\ \multicolumn{1}{c}{M_V}& 0\hfill & 0\hfill & M_{\alpha }\hfill \\ \multicolumn{1}{c}{\frac{-L_V}{V_0}}& 0\hfill & 1\hfill & \frac{-L_{\alpha }}{V_0}\hfill \end{array}\right]\begin{array}{c}{\mathbf{B}}_0=\left[\begin{array}{cc}0\hfill & T_{\delta T}\hfill \\ \multicolumn{1}{c}{0}& 0\hfill \\ \multicolumn{1}{c}{M_{\delta E}}& 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill \end{array}\right]\hfill \end{array}$$

(2)

The elements are simplified for discussion under the assumptions of steady flow and level flight. ²⁶ The aerodynamic derivatives and notations are defined according to Stengel’s book. ²⁵

Phugoid and short-period decoupling

The state equations are generally decomposed into two blocks as follows

$$\left[\begin{array}{c}{\stackrel{·}{\mathbf{X}}}_1\hfill \\ \multicolumn{1}{c}{{\stackrel{·}{\mathbf{X}}}_2}\end{array}\right]=\left[\begin{array}{cc}{\mathbf{A}}_{11}\hfill & {\mathbf{A}}_{12}\hfill \\ \multicolumn{1}{c}{{\mathbf{A}}_{21}}& {\mathbf{A}}_{22}\hfill \end{array}\right]\left[\begin{array}{c}{\mathbf{X}}_1\hfill \\ \multicolumn{1}{c}{{\mathbf{X}}_2}\end{array}\right]+\left[\begin{array}{c}{\mathbf{B}}_1\hfill \\ \multicolumn{1}{c}{{\mathbf{B}}_2}\end{array}\right]\mathbf{U}$$

(3)

where the state vector is regrouped into ${\mathbf{X}}_1^{\mathrm{T}}$ = [V γ] and ${\mathbf{X}}_2^{\mathrm{T}}$ = [q α]. In classical flight mechanics, X₁ represents the the phugoid mode and X₂ the short-period mode. These two modes are decoupled in terms of the response time scale. The short-period mode describes the dynamics of the pitch rate and the angle of attack, which is often regulated as an inner loop of the longitudinal controllers. Therefore, the phugoid mode is responsible for altitude and airspeed regulation.

There are two methods to decouple the dynamics. ²⁵ One is truncated, which neglects all the coupling sub-matrices in A as the following form

$$\left[\begin{array}{c}{\stackrel{·}{\mathbf{X}}}_1\hfill \\ \multicolumn{1}{c}{{\stackrel{·}{\mathbf{X}}}_2}\end{array}\right]=\left[\begin{array}{cc}{\mathbf{A}}_{11}\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& {\mathbf{A}}_{22}\hfill \end{array}\right]\left[\begin{array}{c}{\mathbf{X}}_1\hfill \\ \multicolumn{1}{c}{{\mathbf{X}}_2}\end{array}\right]+\left[\begin{array}{cc}0\hfill & {\mathbf{B}}_{12}\hfill \\ \multicolumn{1}{c}{{\mathbf{B}}_{21}}& 0\hfill \end{array}\right]\left[\begin{array}{c}\delta E\hfill \\ \multicolumn{1}{c}{\delta T}\end{array}\right]$$

(4)

The other method is residualization, which assigns the time derivatives to zero for the short-period mode as follows

$$\left[\begin{array}{c}{\stackrel{·}{\mathbf{X}}}_1\hfill \\ \multicolumn{1}{c}{0}\end{array}\right]=\left[\begin{array}{cc}{\mathbf{A}}_{11}\hfill & {\mathbf{A}}_{12}\hfill \\ \multicolumn{1}{c}{{\mathbf{A}}_{21}}& {\mathbf{A}}_{22}\hfill \end{array}\right]\left[\begin{array}{c}{\mathbf{X}}_1\hfill \\ \multicolumn{1}{c}{{\mathbf{X}}_2}\end{array}\right]+\left[\begin{array}{cc}{\mathbf{B}}_{11}\hfill & {\mathbf{B}}_{12}\hfill \\ \multicolumn{1}{c}{{\mathbf{B}}_{21}}& 0\hfill \end{array}\right]\left[\begin{array}{c}\delta E\hfill \\ \multicolumn{1}{c}{\delta T}\end{array}\right]$$

(5)

Thus, the phugoid mode equations can be rewritten as

$$\begin{array}{c}\multicolumn{1}{c}{{\stackrel{·}{\mathbf{X}}}_1=({\mathbf{A}}_{11}-{\mathbf{A}}_{12}{\mathbf{A}}_{22}^{-1}{\mathbf{A}}_{21}){\mathbf{X}}_1}\\ \multicolumn{1}{c}{+({\mathbf{B}}_{11}-{\mathbf{A}}_{12}{\mathbf{A}}_{22}^{-1}{\mathbf{B}}_{21})\delta E+{\mathbf{B}}_{12}\delta T}\end{array}$$

(6)

where A₂₂ is nonsingular for an ordinary aircraft, and the short-period equation is simply partitioned as

$${\stackrel{·}{\mathbf{X}}}_2={\mathbf{A}}_{22}{\mathbf{X}}_2+{\mathbf{B}}_{21}\delta E$$

(7)

The two methods have the same short-period approximation and little difference in terms of the phugoid mode eigenvalues. Both of these simplifications as above, however, lack rigorous theoretical proof.

The decoupled structure is beneficial for controller design. A cascade controller is widely used in classical flight control as shown in Figure 2.

Figure 2.

An ordinary structure of the longitudinal dynamics.

In Figure 2, the short-period mode is controlled by the elevators, and the relevant controller is called the inner-loop controller. Hence, the phugoid mode can be regulated by the inner loop and the thrust. The entire closed system consists of two cascade stages. The components covered by gray square indicate the inner system.

Example of longitudinal decoupling

The decoupling is generally quite accurate, even though the coupling modes are neglected. The following example will show the validity of this approach.

Assume that a transport aircraft F-16 fighter cruises at the airspeed of 250 ft/s at low altitude. The longitudinal aerodynamic parameters in this situation are listed in Table 1. The data are taken from Lavretsky and Wise’s ²⁷ book.

Table 1.

F-16 longitudinal aerodynamics derivatives.

Derivative	Value	Derivative	Value
DV	0.038	MV	0.0001
Mα	–0.759	Lα	158.0
LV	0.25	MδE	–0.011
Dα	–51.158	TδT	10.1

Figure 3 shows the exact eigenvalues and the approximated modes on the same complex plane. It depicts that the decoupling method is able to approximate the exact dynamics.

Figure 3.

Comparison of the exact modes and the approximated modes.

Singular perturbation issue

Normalized equations

Define a normalized time as

$$t_r\stackrel{\Delta }{=}{D}_{V}t$$

(8)

and scale the velocity with V₀

$$V_r\stackrel{\Delta }{=}V/V_0$$

(9)

Therefore, the normalization equations can be derived from equation (1) as

$$\frac{d{\mathbf{X}}_r}{d{t}_r}={\mathbf{A}}_{r0}{\mathbf{X}}_r+{\mathbf{B}}_{r0}{\mathbf{U}}_r$$

(10)

where

$${\mathbf{A}}_{r0}=\left[\begin{array}{cccc}-1\hfill & -\frac{g}{D_V{V}_0}\hfill & 0\hfill & -\frac{D_{\alpha }+g}{D_V{V}_0}\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 1\hfill & 0\hfill \\ \multicolumn{1}{c}{\frac{M_V{V}_0}{D_V^2}}& 0\hfill & 0\hfill & \frac{M_{\alpha }}{D_V^2}\hfill \\ \multicolumn{1}{c}{-\frac{L_V}{D_V}}& 0\hfill & 1\hfill & -\frac{L_{\alpha }}{D_V{V}_0}\hfill \end{array}\right]$$

(11)

$${\mathbf{B}}_{r0}^{\mathrm{T}}=\left[\begin{array}{cccc}0\hfill & 0\hfill & \frac{M_{\delta E}}{D_V^2}\hfill & 0\hfill \\ \multicolumn{1}{c}{\frac{T_{\delta T}}{D_V{V}_0}}& 0\hfill & 0\hfill & 0\hfill \end{array}\right]$$

(12)

and

$${\mathbf{X}}^{\mathrm{T}}\stackrel{\Delta }{=}{\left[\begin{array}{cccc}{V}_r\hfill & \theta \hfill & q_r\hfill & \alpha \hfill \end{array}\right]}^{\mathrm{T}}$$

(13)

$${\mathbf{U}}_r^{\mathrm{T}}\stackrel{\Delta }{=}{\left[\begin{array}{cc}\delta E\hfill & \delta T_r\hfill \end{array}\right]}^{\mathrm{T}}$$

(14)

where q_r = q/D_V. In the equations above, physical variables are scaled and the longitudinal equations become dimensionless.

Standard singular perturbation problem

The description of the singular perturbation theory is from Khalil and Grizzle’s ²⁸ book. The standard singular perturbation problem can be expressed in the following form

$$\begin{array}{c}\stackrel{·}{x}=f(t,x,z,\epsilon )\hfill \\ \epsilon \stackrel{·}{z}=g(t,x,z,\epsilon )\hfill \end{array}$$

(15)

where ε denotes a small perturbation constant. When ε approaches 0, the second differential equation becomes an algebraic one as follows

$$0=g(t,x,z,0)$$

(16)

If the solution of equation (16) can be obtained in the following form

$$z=h(t,x)$$

(17)

Equation (15) can be reduced to

$$\stackrel{·}{x}=f(t,x,h(t,x),0)$$

(18)

which is called the slow model, while the second equation in equation (15) is named the fast model.

To clarify such an approximation due to the time scale, subtract h(t, x) from z as

$$y\stackrel{\Delta }{=}z-h(t,x)$$

(19)

The underlying system of equation (15) is transformed into

$$\begin{array}{c}\multicolumn{1}{c}{\stackrel{·}{x}=f(t,x,y+h(t,x),\epsilon )}\\ \multicolumn{1}{c}{\epsilon \stackrel{·}{y}=g(t,x,y+h(t,x),\epsilon )-\epsilon \frac{\partial h}{\partial t}}\\ \multicolumn{1}{c}{-\epsilon \frac{\partial h}{\partial x}f(t,x,y+h(t,x),\epsilon )}\end{array}$$

Define a fast time variable τ = t/ε, and the second equation in equation (20) will be rewritten as

$$\begin{array}{c}\multicolumn{1}{c}{\frac{\mathrm{dy}}{d\tau }=g(t_0+\epsilon \tau ,x,y+h(t_0+\epsilon \tau ,x),\epsilon )-\epsilon \frac{\partial h}{\partial t}}\\ \multicolumn{1}{c}{-\epsilon \frac{\partial h}{\partial x}f(t_0+\epsilon \tau ,x,y+h(t_0+\epsilon \tau ,x),\epsilon )}\end{array}$$

where t₀ is the initial time, and given ε = 0, we can obtain the boundary-layer equation

$$\frac{\mathrm{dy}}{d\tau }=g(t_0,x_0,y+h(t_0,x_0),0)$$

(22)

which indicates that if ε is quite small, the system will converge quickly to the manifold z = h(t, x) because of the short time scale.

In the viewpoint of geometry, the system dynamics can be approximated by the fast manifold and the slow manifold. Equation (22) is an expression of the fast manifold, and the characteristics of the entire system are dominated by the first equation of equation (20) after it arrives at the slow manifold.

Longitudinal equation transformation

Consider the following second-order subsystem of the normalized longitudinal dynamics

$$\frac{d}{d{t}_r}\left[\begin{array}{c}{q}_r\hfill \\ \multicolumn{1}{c}{\alpha }\end{array}\right]=\left[\begin{array}{cc}0\hfill & -a_2\hfill \\ \multicolumn{1}{c}{1}& -a_1\hfill \end{array}\right]\left[\begin{array}{c}{q}_r\hfill \\ \multicolumn{1}{c}{\alpha }\end{array}\right]+\left[\begin{array}{cc}{a}_3\hfill & a_5\hfill \\ \multicolumn{1}{c}{-a_4}& 0\hfill \end{array}\right]\left[\begin{array}{c}{V}_r\hfill \\ \multicolumn{1}{c}{\delta E}\end{array}\right]$$

(23)

where

$$\begin{array}{c}{a}_1\stackrel{\Delta }{=}\frac{L_{\alpha }}{D_V{V}_0},a_2\stackrel{\Delta }{=}-\frac{M_{\alpha }}{D_V^2},a_3\stackrel{\Delta }{=}\frac{M_V{V}_0}{D_V^2}\hfill \\ a_4\stackrel{\Delta }{=}\frac{L_V}{D_V},a_5\stackrel{\Delta }{=}\frac{M_{\delta E}}{D_V^2}\hfill \end{array}$$

(24)

The corresponding characteristic matrix is

$$s\mathbf{I}-\mathbf{A}=\left[\begin{array}{cc}s\hfill & a_2\hfill \\ \multicolumn{1}{c}{-1}& s+a_1\hfill \end{array}\right]$$

(25)

where A denotes the underlying system matrix. Hence, the first-order and second-order determinant factors are

$$D_1\left(s\right)=1,D_2\left(s\right)=s^2+a_{1}s+a_2$$

(26)

Assuming a₂ no less than $a_1^2/4$ , which means that the system is underdamped, there is another matrix

$$\tilde{\mathbf{A}}=\left[\begin{array}{cc}-\frac{a_1}{2}\hfill & \sqrt{a_2-\frac{a_1^2}{4}}\hfill \\ \multicolumn{1}{c}{-\sqrt{a_2-\frac{a_1^2}{4}}}& -\frac{a_1}{2}\hfill \end{array}\right]$$

(27)

$s\mathbf{I}-\tilde{\mathbf{A}}$ has the same determinant factors with sI-A, which proves that

$$s\mathbf{I}-\mathbf{A}\approx s\mathbf{I}-\tilde{\mathbf{A}}$$

(28)

where the symbol ≈ means that two polynomial matrices are equivalent. This equivalence can result in that A which is similar to $\tilde{\mathbf{A}}$ . ²⁹ Furthermore, there exists a transformation matrix P such that

$${\mathbf{P}}^{-1}\mathbf{AP}=\tilde{\mathbf{A}}$$

(29)

Therefore, equation (23) can be transformed into

$$\frac{d}{d{t}_r}\left[\begin{array}{c}{x}_1\hfill \\ \multicolumn{1}{c}{x_2}\end{array}\right]=\tilde{\mathbf{A}}\left[\begin{array}{c}{x}_1\hfill \\ \multicolumn{1}{c}{x_2}\end{array}\right]+{\mathbf{P}}^{-1}\left[\begin{array}{cc}{a}_3\hfill & a_5\hfill \\ \multicolumn{1}{c}{-a_4}& 0\hfill \end{array}\right]\left[\begin{array}{c}{V}_r\hfill \\ \multicolumn{1}{c}{\delta E}\end{array}\right]$$

(30)

where [q_rα]^T = P[x₁x₂]^T. Matrix equation (29) can be expanded into the following homogeneous linear algebraic equations

$$\begin{array}{c}{\left[\begin{array}{cccc}0\hfill & 0\hfill & 0\hfill & 0\hfill \end{array}\right]}^{\mathrm{T}}=\hfill \\ \left[\begin{array}{cccc}-\frac{a_1}{2}\hfill & -1\hfill & \sqrt{a_2-\frac{a_1^2}{4}}\hfill & 0\hfill \\ \multicolumn{1}{c}{a_2}& \frac{a_1}{2}\hfill & 0\hfill & \sqrt{a_2-\frac{a_1^2}{4}}\hfill \\ \multicolumn{1}{c}{-\sqrt{a_2-\frac{a_1^2}{4}}}& 0\hfill & -\frac{a_1}{2}\hfill & -1\hfill \\ \multicolumn{1}{c}{0}& -\sqrt{a_2-\frac{a_1^2}{4}}\hfill & a_2\hfill & \frac{a_1}{2}\hfill \end{array}\right]\left[\begin{array}{c}{q}_{11}\hfill \\ \multicolumn{1}{c}{q_{12}}\\ \multicolumn{1}{c}{q_{21}}\\ \multicolumn{1}{c}{q_{22}}\end{array}\right]\hfill \end{array}$$

(31)

where the solutions are not unique. Pick one particular solution and rearrange the elements to express P^–1 as follows

$$\left[\begin{array}{c}{q}_{11}\hfill \\ \multicolumn{1}{c}{q_{12}}\\ \multicolumn{1}{c}{q_{21}}\\ \multicolumn{1}{c}{q_{22}}\end{array}\right]=\left[\begin{array}{c}-1/\sqrt{a_2-\frac{a_1^2}{4}}\hfill \\ \multicolumn{1}{c}{a_1/\left(2\sqrt{a_2-\frac{a_1^2}{4}}\right)}\\ \multicolumn{1}{c}{0}\\ \multicolumn{1}{c}{1}\end{array}\right],{\mathbf{P}}^{-1}=\left[\begin{array}{cc}{q}_{11}\hfill & q_{12}\hfill \\ \multicolumn{1}{c}{q_{21}}& q_{22}\hfill \end{array}\right]$$

(32)

Thus, equation (30) can be re-expressed as

$$\frac{d}{d{t}_r}\left[\begin{array}{c}{x}_1\hfill \\ \multicolumn{1}{c}{x_2}\end{array}\right]=\tilde{\mathbf{A}}\left[\begin{array}{c}{x}_1\hfill \\ \multicolumn{1}{c}{x_2}\end{array}\right]+\left[\begin{array}{cc}{q}_{11}{a}_3-q_{12}{a}_4\hfill & q_{11}{a}_5\hfill \\ \multicolumn{1}{c}{-a_4}& 0\hfill \end{array}\right]\left[\begin{array}{c}{V}_r\hfill \\ \multicolumn{1}{c}{\delta E}\end{array}\right]$$

(33)

To derivate the singular perturbation form, make the following polar transformation

$$x_1=r\cos \lambda ,x_2=r\sin \lambda$$

(34)

Combined with equations (33) and (34), the underlying system (30) can be transformed furthermore into

$$\begin{array}{c}\frac{\mathrm{dr}}{d{t}_r}=-\frac{a_1}{2}r+\cos \lambda ((q_{11}{a}_3-q_{12}{a}_4)V_r+q_{11}{a}_5\delta E)\hfill \\ -\frac{a_4{V}_r\sin \lambda }{a_1}\hfill \\ \frac{d\lambda }{d{t}_r}=-\sqrt{a_2-\frac{a_1^2}{4}}-\frac{a_4\cos \lambda +(q_{11}{a}_3-q_{12}{a}_4)\sin \lambda }{r}{V}_r\hfill \\ -\frac{q_{11}{a}_5\sin \lambda }{r}\delta E\hfill \end{array}$$

(35)

To analyze the open-loop response, assign δE to 0, because feedback will change the time constant and destroy the singular perturbation structure. Select the perturbation constant ε = 1/a₁, and equation (35) will be rewritten as the following singular perturbation model

$$\begin{array}{c}\epsilon \stackrel{·}{r}=-\frac{1}{2}r+V_{r}A\cos (\lambda +\varphi )\hfill \\ \epsilon \stackrel{·}{\lambda }=-\sqrt{\frac{a_2}{a_1^2}-\frac{1}{4}}-\frac{A\sin (\lambda +\varphi )}{r}{V}_r\hfill \end{array}$$

(36)

where

$$\begin{array}{c}A\stackrel{\Delta }{=}\sqrt{\frac{{(q_{11}{a}_3-q_{12}{a}_4)}^2+a_4^2}{a_1^2}}\hfill \\ \varphi \stackrel{\Delta }{=}\arctan \frac{a_4}{q_{11}{a}_3-q_{12}{a}_4}\hfill \end{array}$$

(37)

Assign the derivative to zero, and the following expression can be obtained according to equations (16) and (17)

$$h\left(V_r\right)=\left[\begin{array}{c}{r}_0\hfill \\ \multicolumn{1}{c}{{\lambda }_0}\end{array}\right]=\left[\begin{array}{c}A{a}_1{a}_2^{-1/2}{V}_r\hfill \\ \multicolumn{1}{c}{-\arctan \sqrt{\frac{4a_2}{a_1^2}-1}-\varphi }\end{array}\right]$$

(38)

and make another transformation as

$$\begin{array}{c}r\cos \lambda =y_r\cos y_{\lambda }+r_0\cos {\lambda }_0\hfill \\ r\sin \lambda =y_r\sin y_{\lambda }+r_0\sin {\lambda }_0\hfill \end{array}$$

(39)

which derives the boundary-layer equations as follows

$$\begin{array}{c}\frac{d{y}_r}{d\tau }=-\frac{1}{2}{y}_r\hfill \\ \frac{d{y}_{\lambda }}{d\tau }=-\sqrt{\frac{a_2}{a_1^2}-\frac{1}{4}}\hfill \end{array}$$

(40)

As boundary-layer model equation (40) is global asymptotically stable for r, it follows that

$$\begin{array}{c}{x}_1=r\cos \lambda \to r_0\cos {\lambda }_0\hfill \\ x_2=r\sin \lambda \to r_0\sin {\lambda }_0\hfill \end{array}$$

(41)

and thus the slow model of system equation (10) is

$$\begin{array}{c}\multicolumn{1}{c}{\frac{d}{d{t}_r}\left[\begin{array}{c}{V}_r\hfill \\ \multicolumn{1}{c}{\theta }\end{array}\right]=\left[\begin{array}{cc}-1\hfill & -\frac{g}{D_V{V}_0}\hfill \\ \multicolumn{1}{c}{0}& 0\hfill \end{array}\right]\left[\begin{array}{c}{V}_r\hfill \\ \multicolumn{1}{c}{\theta }\end{array}\right]+\left[\begin{array}{c}{T}_{\delta T}\hfill \\ \multicolumn{1}{c}{0}\end{array}\right]T}\\ \multicolumn{1}{c}{+\left[\begin{array}{cc}0\hfill & -\frac{g}{D_V{V}_0}\hfill \\ \multicolumn{1}{c}{1}& 0\hfill \end{array}\right]\mathbf{P}\left[\begin{array}{c}{r}_0\cos {\lambda }_0\hfill \\ \multicolumn{1}{c}{r_0\sin {\lambda }_0}\end{array}\right]}\end{array}$$

which is equivalent to the residualization method. It proves that equations (6) and (7) provide the same decoupling approach with singular perturbation.

Simulation

For the example in section “Example of longitudinal decoupling,” its relevant matrices in the normalized equations are calculated as

$$\begin{array}{c}{\mathbf{A}}_r=\left[\begin{array}{cccc}-1\hfill & -0.1287\hfill & 0\hfill & -1.9983\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 1\hfill & 0\hfill \\ \multicolumn{1}{c}{0.0001}& 0\hfill & 0\hfill & -525.6233\hfill \\ \multicolumn{1}{c}{-6.5789}& 0\hfill & 1\hfill & -16.6316\hfill \end{array}\right]\hfill \\ {\mathbf{B}}_r^{\mathrm{T}}=\left[\begin{array}{cccc}0\hfill & 0\hfill & -7.6177\hfill & 0\hfill \\ \multicolumn{1}{c}{1.0632}& 0\hfill & 0\hfill & 0\hfill \end{array}\right]\hfill \end{array}$$

(43)

where the modes are shown in Figure 4.

Figure 4.

Comparison of the exact modes and the approximated modes for the normalized equations.

We can see from Figure 4 that the approximated modes almost have the same eigenvalues with the exact ones. The approximation accuracy can also be checked in the time domain from Figures 5 and 6.

Figure 5.

Impulse responses of the exact dynamics and the phugoid approximation: (a) velocity plot and (b) pitch angle plot.

Figure 6.

Impulse responses of the exact dynamics and the short-time approximation: (a) pitch rate plot and (b) angle of attack plot.

The eigenvalue approximations are reflected in the time domain as shown in Figures 5 and 6. Figure 5 shows that the phugoid dynamics has a perfect approximation in the natural frequency although the damp ratio is a little smaller than that of the exact dynamics. Figure 6 shows that the short-time reduced equations provide a perfect approximation in terms of time response. Figure 5 shows that the phugoid reduced dynamics give an acceptable response except for the damping ratio.

Conclusion

The aircraft longitudinal dynamics are ordinary differential equations with four state variables. In the classical flight control theory, the equations are generally partitioned into two modes, phugoid mode and short-time mode, but there is no rigorous theoretical support to this decoupling. We present a decoupling approach based on the singular perturbation theory. The underlying equations are transformed into a canonical form. The perturbation parameter is abstracted from the canonical form and generates a perturbation problem. The short-time mode and phugoid mode is, therefore, naturally decoupled. Finally, the simulation results depict that the proposed decoupling method can provide an approximation of high quality to the exact dynamics. In the future, the adaptive control for the decoupling model of aircraft flight mechanics deserves further discussion.