Error Dynamics Based Dual Heuristic Dynamic Programming for Self-Learning Flight Control

Abstract

A data-driven nonlinear control approach, called error dynamics-based dual heuristic dynamic programming (ED-DHP), is proposed for air vehicle attitude control. To solve the optimal tracking control problem, the augmented system is defined by the derived error dynamics and reference trajectory so that the actor neural network can learn the feedforward and feedback control terms at the same time. During the online self-learning process, the actor neural network learns the control policy by minimizing the augmented system’s value function. The input dynamics identified by the recursive least square (RLS) and output of the critic neural network are used to update the actor neural network. In addition, the total uncertainty term of the error dynamics is also identified by RLS, which can compensate for the uncertainty caused by inaccurate modeling, parameter perturbation, and so on. The outputs of ED-DHP include the rough trim surface, feedforward and feedback terms from the actor neural network, and the compensation. Based on this control scheme, the complete knowledge of system dynamics and the reference trajectory dynamics are not needed, and offline learning is unnecessary. To verify the self-learning ability of ED-DHP, two numerical experiments are carried out based on the established morphing air vehicle model. One is sinusoidal signal tracking at a fixed operating point, and the other is guidance command tracking with a morphing process at variable operating points. The simulation results demonstrate the good performance of ED-DHP for online self-learning attitude control and validate the robustness of the proposed scheme

1. Introduction

The flight control system plays an essential role in guaranteeing the flight quality of air vehicles. However, the nonlinearity and model uncertainty of the vehicle dynamics become more significant because of increasing flight mission complexity [1,2]. Besides, it is difficult to model complex vehicle systems, and great effort is required in model identification. For these reasons, the demand for adaptive nonlinear attitude control methods with modeling inaccuracies is rising.

Reinforcement learning (RL) is a machine learning method. It possesses a self-learning capability by interacting with the environment. In recent years, RL has been applied to solve a large number of control problems in many fields, such as robotics [3], nuclear fusion [4], intelligent transportation [5], and smart grid [6]. RL is consistently used as a parameter adjuster [7], uncertainty compensator [8], and end-to-end controller [9,10] in the flight control system as well.

Adaptive/approximate dynamic programming (ADP) is one of the RL methods, bridging the gap between traditional optimal and adaptive control [11]. Different frameworks with ADP have been proposed to solve adaptive flight control problems. In some frameworks, linear/nonlinear control methods such as PID [12], sliding mode control (SMC) [13,14], nonlinear dynamic inverse (NDI) [15], and backstepping (BS) [16] are combined with ADP to improve the control quality. Moreover, ADP provides supplementary adjustments to expedite control performance. However, the linear/nonlinear controllers are dominant parts of the control system, so the system performance can be bad or even unstable without these controllers. Moreover, they are always designed with known dynamic models and may not work when high nonlinearity and large model uncertainty exist. Some other frameworks where the ADP-based controllers are dominant have been proposed. For example, Han et al. [17] proposed an online policy iteration (PI) method for adaptive control during the vehicle reentry phase. But the convergence condition of PI is that the initial control policy must make the system stable, so the neural networks in this method should be trained offline. In addition, an ADP method called adaptive critic designs (ACDs) [18,19] have been extensively studied. However, extra structures, such as the global system model, are also vital to approximate the system dynamics in ACDs. For this reason, an offline learning phase is additionally needed in system model identification.

To reduce the dependency on prior knowledge of the system model and offline phase, one framework called incremental model-based dual heuristic dynamic programming (iDHP) was proposed by Zhou et al. [20,21]. The critic neural network (NN) in iDHP approximates the costate function, which is a state dimension vector gradient and can bring more information than the value function. In iDHP, the incremental model is identified by the RLS method, and the incremental model is provided for the neural network updating. It should be noted that the RLS works with high-frequency sample data and relatively slow-varying system assumption. Given the relatively high training failure rate in the practical application of iDHP, Li et al. [22] proposed a switching control system consisting of an iDHP-based controller and a PID controller, which was introduced as a safety control policy. Furthermore, to improve the iDHP’s adaptability to a wide range of dynamic variations, Liu et al. [23] proposed an adaptive learning rate-based iDHP method where a cost function was designed to evaluate the controller stability based on the wavelet analysis. Lee et al. [24] presented two iDHP controller designs with and without the cascaded actor structure to realize the air vehicle altitude control. Heyer et al. [25] optimized the iDHP with an additional network structure to enhance the learning stability and applied the method to the high-dimensional aerospace models.

To further optimize the identification ability of RLS in iDHP and make full use of the iDHP’s advantage in learning efficiency, a novel approach called error dynamics-based dual heuristic dynamic programming is proposed for online nonlinear flight control. One merit of ED-DHP is that complete knowledge of the vehicle dynamics and reference trajectory dynamics is unnecessary. Unlike iDHP, RLS is introduced to identify error dynamics containing the input distribution matrix and the total disturbance term. What’s more, the state of ED-DHP is augmented by tracking error and reference trajectory, so the actor NN can learn the feedforward and feedback control terms. The contributions of this paper are summarized below:

• An ED-DHP controller for the fast loop is designed to control the morphing air vehicle attitude system with different work points without offline learning;
• The error dynamics of the fast loop are derived. Based on these, the control effectiveness matrix and the total disturbance term are identifiable by RLS. This enhances learning efficiency effectively. The ADP system is augmented by tracking error and reference trajectory, so the feedforward and feedback control terms can be learned by the actor NN, on the basis that the complete knowledge of system dynamics and the reference trajectory dynamics are unknown;
• The outputs of ED-DHP include the rough trim surface, feedforward and feedback terms, and compensation. The robustness of the method can be enhanced.

The remainder of this paper is organized as follows: Section 2 describes the attitude tracking problem of the air vehicle; Section 3 presents the control structure of the system, including design and analysis of the ED-DHP framework; In Section 4, some simulation results are given to demonstrate the effectiveness of ED-DHP; Section 5 states the conclusion of this article.

2. Mathematical Model of Air Vehicle

The kinematical and dynamical equations of the air vehicle are shown below.

$$\dot{\Omega }=g_s\omega +\frac{1}{mv}{F}_s$$

(1)

$$\dot{\omega }=-{\left(I^{\sigma (t)}\right)}^{-1}\mu I^{\sigma \left(t\right)}\omega +{\left(I^{\sigma (t)}\right)}^{-1}M$$

(2)

where $\Omega ={[\alpha ,\beta ,{\gamma }_V]}^{\mathrm{T}}$ and $\omega ={[{\omega }_x,\hspace{0.33em}{\omega }_y,\hspace{0.33em}{\omega }_z]}^{\mathrm{T}}$ are the system state vectors. $\alpha$, $\beta$, and ${\gamma }_V$ denote the attack angle, sideslip angle, and bank angle separately; $\omega$ is the angular rate vector with roll rate ${\omega }_x$, yaw rate ${\omega }_y$, and pitch rate ${\omega }_z$. $m$ and $v$ are the mass and the velocity of the vehicle.

For the slow loop subsystem (1), the concrete expressions of $g_s$ and $F_s$ are as follows:

$$\hspace{0.33em}{g}_s=\left[\begin{array}{ccc}-\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{\alpha }\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{\beta }& \mathrm{s}\mathrm{i}\mathrm{n}\mathrm{\alpha }\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{\beta }& 1\\ \mathrm{s}\mathrm{i}\mathrm{n}\mathrm{\alpha }& \mathrm{c}\mathrm{o}\mathrm{s}\mathrm{\alpha }& 0\\ \mathrm{c}\mathrm{o}\mathrm{s}\mathrm{\alpha }\mathrm{s}\mathrm{e}\mathrm{c}\mathrm{\beta }& -\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{\alpha }\mathrm{s}\mathrm{e}\mathrm{c}\mathrm{\beta }& 0\end{array}\right], \hspace{0.33em}\hspace{0.33em}{F}_s=\left[\begin{array}{c}-F_y/\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{\beta }\\ F_z\\ (F_y\mathrm{c}\mathrm{o}\mathrm{s}{\gamma }_V-F_z\mathrm{s}\mathrm{i}\mathrm{n}{\gamma }_V)\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{\theta }+F_y\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{\beta }\end{array}\right]$$

For the fast loop subsystem (2), the concrete expressions of $I^{\sigma \left(t\right)}$, $\mu$, and $M$ are as follows:

$$I^{\sigma \left(t\right)}=\left[\begin{array}{ccc}{I}_x^{\sigma (t)}& -I_{xy}^{\sigma (t)}& -I_{xz}^{\sigma (t)}\\ -I_{xy}^{\sigma (t)}& I_y^{\sigma (t)}& -I_{yz}^{\sigma (t)}\\ -I_{zx}^{\sigma (t)}& -I_{zy}^{\sigma (t)}& I_z^{\sigma (t)}\end{array}\right], \mu =\left[\begin{array}{ccc}0& -{\omega }_z^{}& {\omega }_y^{}\\ {\omega }_z^{}& 0& -{\omega }_x^{}\\ -{\omega }_y^{}& {\omega }_x^{}& 0\end{array}\right],$$

$$M=\left[\begin{array}{c}{M}_x\\ M_y\\ M_z\end{array}\right]=\left[\begin{array}{c}{C}_{mx}^{\sigma (t)}\left(Ma,\alpha ,\beta ,\delta \right)qSl\\ C_{my}^{\sigma (t)}\left(Ma,\alpha ,\beta ,\delta \right)qSl\\ C_{mz}^{\sigma (t)}\left(Ma,\alpha ,\beta ,\delta \right)qSl\end{array}\right]=\left[\begin{array}{c}\left(C_{mx\_basic}^{\sigma (t)}\left(Ma,\alpha ,\beta \right)+C_{mx\_\delta }^{\sigma (t)}\left(Ma,\alpha ,\beta ,\delta \right)\right)qSl\\ \left(C_{my\_basic}^{\sigma (t)}\left(Ma,\alpha ,\beta \right)+C_{my\_\delta }^{\sigma (t)}\left(Ma,\alpha ,\beta ,\delta \right)\right)qSl\\ \left(C_{mz\_basic}^{\sigma (t)}\left(Ma,\alpha ,\beta \right)+C_{mz\_\delta }^{\sigma (t)}\left(Ma,\alpha ,\beta ,\delta \right)\right)qSl\end{array}\right]$$

where $\sigma (t)$ describes different morphing configurations of the air vehicle; $\theta$ is the trajectory inclination angle; $F_y$ is the normal force and $F_z$ is the lateral force; $I^{\sigma \left(t\right)}$ is the moment of inertia matrix; $q=\rho v^2/2$ denotes the dynamic pressure; $S$ and $l$ are the reference area and reference aerodynamic chord, respectively. $C_{m\ast \_basic}^{\sigma (t)}\left(Ma,\alpha ,\beta \right)$ is the aerodynamic moment coefficient, and $C_{m\ast \_\delta }^{\sigma (t)}\left(Ma,\alpha ,\beta ,\delta \right)$ is the control moment coefficient induced by the control surface; $Ma$ is the Mach number and $\delta$ is the control surface vector.

Remark. 

the research object is a morphing air vehicle, and the moment of inertia ${\overline{I}}^{\sigma \left(t\right)}$ and the aerodynamic coefficient $C_{m\ast }^{\sigma (t)}$ change with different configurations, which can be derived by linear interpolation.

However, uncertainties always present for the nominal parameters during flight. Therefore, to ensure the robustness of the air vehicle’s flight control, the uncertain parameters in systems (1) and (2) cannot be neglected, described as:

$$\left\{\begin{array}{l}I={\overline{I}}^{\sigma (t)}(1+\mathrm{\Delta }I)\\ m=\overline{m}(1+\mathrm{\Delta }m)\\ \mathrm{\rho }=\overline{\mathrm{\rho }}(1+\mathrm{\Delta }\rho )\\ C_{m\ast }^{\sigma (t)}={\overline{C}}_{m\ast }^{\sigma (t)}(1+\mathrm{\Delta }{C}_m)\end{array}\right.$$

(3)

where the expressions $\left(\overline{•}\right)$ are the nominal value of the key parameters, and the parametric uncertainties are assumed to be bounded. Uncertain parameters can affect the system modeling and control effect. This uncertainty model will be applied in chapter 4 to verify the effectiveness and robustness of ED-DHP.

3. Online Tracking Controller Based on ED-DHP

3.1. Control Structure of the System

Based on the multiple-timescale features, the dynamics can be divided into the slow loop subsystem and the fast loop subsystem. So the controller for each system can be designed separately. As illustrated in Figure 1, the slow loop subsystem control law is designed by dynamic inversion with the tracking differentiator (TD) in this article to get the virtual control input of the slow loop subsystem. The ED-DHP is applied to angular rate control of the pitch, yaw, and roll rate, augmented with the slow outer loop. Rate control exhibits the lowest learning complexity because the angular rates and control surfaces are under direct dynamic relation. The control surfaces consist of the rough trim value, the output of actor NN, and the uncertainty compensator proposed in Section 3.2.

The virtual control vector of the slow loop is designed by dynamic inversion:

$${\omega }_{ref}=g_s{}^{-1}\left(\left[\begin{array}{c}{\dot{\alpha }}_{ref}+k_{\alpha }({\alpha }_{ref}-\alpha )\\ {\dot{\beta }}_{ref}+k_{\beta }({\beta }_{ref}-\beta )\\ {\dot{\gamma }}_{Vref}+k_{{\gamma }_V}({\gamma }_{Vref}-{\gamma }_V)\end{array}\right]-f_s{C}_b^a\left(\left[\begin{array}{c}{a}_x\\ a_y\\ a_z\end{array}\right]+C_I^b\left[\begin{array}{c}{g}_x\\ g_y\\ g_z\end{array}\right]\right)\right)$$

(4)

where ${\omega }_{ref}={[{\omega }_{x,ref}\hspace{0.33em}\hspace{0.33em}{\omega }_{y,ref}\hspace{0.33em}\hspace{0.33em}{\omega }_{z,ref}]}^{\mathrm{T}}$ is the reference trajectory vector of the fast loop subsystem. ${[a_x\hspace{0.33em}\hspace{0.33em}{a}_y\hspace{0.33em}\hspace{0.33em}{a}_z]}^{\mathrm{T}}$ is the acceleration vector, which can be measured by an accelerometer. ${[g_x\hspace{0.33em}\hspace{0.33em}{g}_y\hspace{0.33em}\hspace{0.33em}{g}_z]}^{\mathrm{T}}$ is the gravity vector. $C_b^a$ and $C_I^b$ are the coordinate transformation matrix calculated by attitude angles with pitch angle $\phi$, yaw angle $\psi$, and roll angle $\gamma$, as well as attack angle $\alpha$ and sideslip angle $\beta$.

$$C_I^b=\left[\begin{array}{ccc}\cos \psi \cos \phi & \cos \psi \sin \phi & -\sin \psi \\ -\cos \gamma \sin \phi +\sin \gamma \sin \psi \cos \phi & \cos \gamma \cos \phi +\sin \gamma \sin \psi \sin \phi & \sin \gamma \cos \psi \\ \sin \gamma \sin \phi +\cos \gamma \sin \psi \cos \phi & -\sin \gamma \cos \phi +\cos \gamma \sin \psi \sin \phi & \cos \gamma \cos \psi \end{array}\right]$$

$$C_b^a=\left[\begin{array}{ccc}\cos \alpha \cos \beta & \sin \beta & \sin \alpha \cos \beta \\ -\cos \alpha \sin \beta & \cos \beta & -\sin \alpha \sin \beta \\ -\sin \alpha & 0& \cos \alpha \end{array}\right]$$

And the matrix function $f_s$ can be written as:

$$f_s=\frac{1}{V}\left[\begin{array}{ccc}0& -\sec \beta & 0\\ 0& 0& 1\\ 0& \tan \beta +\cos {\gamma }_V\tan \theta & -\sin {\gamma }_V\tan \theta \end{array}\right]$$

As we can see from the virtual control vector (4), the derivative of guidance commands is necessary to improve the dynamic performance of the system. To alleviate the “items of explosion” problem, the tracking differentiator was introduced [26], which is described as:

$$\left\{\begin{array}{l}{\dot{x}}_1=x_2\\ {\dot{x}}_2=-r^2\mathrm{sign}\left(x_1-\zeta \right){\left|x_1-\zeta \right|}^{\tau }-r{x}_2\end{array}\right.$$

(5)

where $\zeta$ is the guidance command, $x_2$ is the differential value. $r$ and $\tau$ are positive parameters to be designed.

3.2. The Fast Loop Error Dynamics of Attitude System and Model Identification

The dynamical equations of the fast loop can be rewritten as:

$$\dot{x}=f_1+f_2+gu+\mathrm{\Delta }f$$

(6)

The symbols in Equation (6) are given as:

$$x=\omega , f_1=-{\left({\overline{I}}^{\sigma (t)}\right)}^{-1}R{\overline{I}}^{\sigma (t)}\omega ,$$

$$f_2={\left({\overline{I}}^{\sigma (t)}\right)}^{-1}\left[\begin{array}{c}{\overline{C}}_{mx\_basic}^{\sigma (t)}\left(Ma,\alpha ,\beta \right)qSl\\ {\overline{C}}_{my\_basic}^{\sigma (t)}\left(Ma,\alpha ,\beta \right)qSl\\ {\overline{C}}_{mz\_basic}^{\sigma (t)}\left(Ma,\alpha ,\beta \right)qSl\end{array}\right], gu={\left({\overline{I}}^{\sigma (t)}\right)}^{-1}\left[\begin{array}{c}{\overline{C}}_{mx\_\delta }^{\sigma (t)}\left(Ma,\alpha ,\beta ,\delta \right)qSl\\ {\overline{C}}_{my\_\delta }^{\sigma (t)}\left(Ma,\alpha ,\beta ,\delta \right)qSl\\ {\overline{C}}_{mz\_\delta }^{\sigma (t)}\left(Ma,\alpha ,\beta ,\delta \right)qSl\end{array}\right]$$

The disturbance term $\mathrm{\Delta }f$ expressed parameter uncertainty and modeling error, which causes the inaccurate trim value and external disturbance.

The trim surface is set as $u_p\approx -g^{-1}{f}_2$, which can be roughly obtained by the newton iteration method offline, and the standard control input $u\in {ℝ}^m$ is given by:

$$u=u_p+u_a$$

(7)

where $u_a$ is the output vector of actor NN in ED-DHP.

Substitute Equation (7) into Equation (6) and subtract the reference trajectory dynamics ${\dot{x}}_{ref}$, we can obtain:

$$\begin{array}{l}\dot{e}=\dot{x}-{\dot{x}}_{ref}=f_1+f_2+gu+\mathrm{\Delta }f-{\dot{x}}_{ref}\\ =f_1+f_2+g(u_p+u_a)+\mathrm{\Delta }f-{\dot{x}}_{ref}\\ =f_1(x)+g{u}_a+\mathrm{\Delta }f-{\dot{x}}_{ref}\end{array}$$

(8)

And the problem can be written in a discrete form by Euler’s method:

$$\left\{\begin{array}{l}e(t+1)\approx e(t)-Z(t)+\mathrm{\Delta }t{f}_1(x(t))+\mathrm{\Delta }tg(t)u_a(t)+\mathrm{\Delta }t\mathrm{\Delta }f\\ x_{ref}(t+1)=\varphi (x_{ref}(t))\\ Z(t)=x_{ref}(t+1)-x_{ref}(t)\end{array}\right.$$

(9)

where $f_1$ is known. $g$, $\mathrm{\Delta }f$ and reference trajectory function $\varphi$ are unrevealed, however.

The first equation in Equation (9) can be rewritten as follows:

$$Y(t)=\varsigma {(t)}^{\mathrm{T}}H(t)$$

(10)

where

$$Y(t)=e(t)-e(t-1)+Z(t-1)-\mathrm{\Delta }t{f}_1(x(t-1)),$$

$$\varsigma (t)=\left[\begin{array}{c}\mathrm{\Delta }tg(t-1)\\ \mathrm{\Delta }t\mathrm{\Delta }f\end{array}\right]=\left[\begin{array}{c}{\varsigma }_1(t)\\ {\varsigma }_2(t)\end{array}\right]\in {ℝ}^{(m+1)\times n}, H(t)=\left[\begin{array}{c}{u}_a(t-1)\\ 1\end{array}\right]\in {ℝ}^{(m+1)\times 1}.$$

$Y(t)$ can be directly calculated from the measured vehicle attitude. $H(t)$ is the information matrix, which is only related to the previous output vector of the actor NN. $\varsigma (t)$ is the parameter matrix, and it can be identified in real-time by the RLS method:

$$\epsilon (t)=Y(t)-\widehat{Y}(t)=\varsigma {(t)}^{\mathrm{T}}H(t)-\widehat{\varsigma }{(t)}^{T}H(t)$$

(11)

$$\widehat{\varsigma }(t)=\widehat{\varsigma }(t-1)+\frac{Cov(t-1)H(t)\epsilon {(t)}^{\mathrm{T}}}{{\gamma }_{\mathrm{RLS}}+H{(t)}^{\mathrm{T}}Cov(t-1)H(t)}$$

(12)

$$Cov(t)=\frac{1}{{\gamma }_{\mathrm{RLS}}}\left(Cov(t-1)-\frac{Cov(t-1)H(t)H{(t)}^{\mathrm{T}}Cov(t-1)}{{\gamma }_{\mathrm{RLS}}+H{(t)}^{\mathrm{T}}Cov(t-1)H(t)}\right)$$

(13)

where $\epsilon (t)\in {ℝ}^{n\times 1}$ is the prediction error and ${\gamma }_{\mathrm{RLS}}$ is the forgetting factor. $Cov(t)\in {ℝ}^{(m+1)\times (m+1)}$ is the estimation covariance matrix.

3.3. ED-DHP Framework

For the system dynamics (9), the perfect control input contains the feedforward input $u_d^{\ast }$ and the feedback input $u_e^{\ast }$ if the system dynamics is known [27].

$$\left\{\begin{array}{l}{u}_d^{\ast }={(\mathrm{\Delta }tg)}^{-1}(\varphi (x_{ref}(t))-\mathrm{\Delta }t{f}_1(x_{ref})-\mathrm{\Delta }t\mathrm{\Delta }f)\\ u_e^{\ast }=-\frac{\gamma }{2}{R}^{-1}{(\mathrm{\Delta }tg(t))}^{\mathrm{T}}\frac{\partial V(e(t+1))}{\partial e(t+1)}\\ u_a^{\ast }=u_d^{\ast }+u_e^{\ast }\end{array}\right.$$

(14)

The symbols in Equation (14) are given as follows. Function $J$ is known as the utility and is selected based on the quadratic energy function in this paper. $V$ is defined as a performance measure function or value function in RL. It represents the sum of expected future reward from the current state. Matrices $Q\in {ℝ}^{n\times n}$ and $R\in {ℝ}^{m\times m}$ are positive definite. $\gamma$ is a discount factor.

$$J(e(t),u_e(t))=\frac{1}{2}\left[e{(t)}^{\mathrm{T}}Qe(t)+u_e{(t)}^{\mathrm{T}}R{u}_e(t)\right]$$

(15)

$$V(t)={\displaystyle \sum _{i=t}^{+\infty }{\gamma }^{(i-t)}J(e(i),u_e(i))}$$

(16)

For learning the control strategy of the attitude system without some prior knowledge, the augmented system is constructed in terms of $e(t)$ and $x_{ref}(t)$ [28,29] as:

$$X(t)=\left[\begin{array}{c}e(t)\\ x_{ref}(t)\end{array}\right]\in {ℝ}^{2n}$$

(17)

$$\left[\begin{array}{c}e(t+1)\\ x_{ref}(t+1)\end{array}\right]=\left[\begin{array}{c}e(t)-Z(t)+\mathrm{\Delta }t{f}_1(x(t))\\ \varphi (x_{ref}(t))\end{array}\right]+\left[\begin{array}{c}\mathrm{\Delta }tg(t)\\ 0\end{array}\right]u_a(t)\equiv F(X(t))+G(t)u_a(t)$$

(18)

So we could have the following controller for the fast loop subsystem as follows: (Here, the second part $u_a$ is the output vector of actor NN, and the third $u_b=-{\widehat{\varsigma }}_1{}^{-1}{\widehat{\varsigma }}_2$ is designed to compensate for the total uncertainty.)

$$\begin{array}{l}{u}^{*}=u_p+u_d^{\ast }+u_e^{\ast }-{(\mathrm{\Delta }tg)}^{-1}(\mathrm{\Delta }t\mathrm{\Delta }f)\\ \approx u_p+\underset{u_a:\mathrm{ED}\text{-}\mathrm{DHP}}{\underbrace{u_d+u_e}}+\underset{u_b:\mathrm{RLS}}{\underbrace{-{\widehat{\varsigma }}_1{}^{-1}{\widehat{\varsigma }}_2}}\\ =u_p+u_a+u_b\end{array}$$

(19)

For the augmented system (18), the cost function $J$ and value function $V$ can be rewritten as:

$$J_X(X(t),u_a(t))=\frac{1}{2}\left[X{(t)}^T{Q}_{X}X(t)+u_a{(t)}^{T}R{u}_a(t)\right]$$

(20)

$$\begin{array}{l}{V}_X(t)={\displaystyle \sum _{i=t}^{+\infty }{\gamma }^{(i-t)}{J}_X(X(t),u_a(t))}\\ =J_X(X(t),u_a(t))+\gamma {\displaystyle \sum _{i=t+1}^{+\infty }{\gamma }^{(i-t-1)}{J}_X(X(i),u_a(i))}\\ =J_X(X(t),u_a(t))+\gamma V_X(t+1)\end{array}$$

(21)

where

$$Q_X=\left[\begin{array}{cc}Q& 0\\ 0& 0\end{array}\right]\in {ℝ}^{2n\times 2n}$$

(22)

Figure 2 shows the schematic of ED-DHP. The critic NN approximates the derivatives of the value function $V(X(t))$ concerning the augmented state $X(t)$, defined as $\lambda (X(t))$.

As can be observed from this figure, the previous derivatives of the value function $\widehat{\lambda }(X(t-1))$ and the augmented state $X(t-1)$ are required to be stored and used in updating the NNs, without predicted ones. The control policy estimated by the actor NN is updated by minimizing the value function. A two-layer neural network with one hidden layer containing 10 neurons is used as the critic NN and actor NN in this article, shown in Figure 3. The hidden-layer activation function is tanh. For the critic NN, $w_{\mathrm{critic}1,ij}$, and $w_{\mathrm{critic}2,jk}$ represent the weights from the ith input node to the jth hidden layer node and the jth hidden node to the kth output node, respectively. For the actor NN, $w_{\mathrm{actor}1,ij}$, and $w_{\mathrm{actor}2,jk}$ represent the weights from the ith input node to the jth hidden layer node and the jth hidden node to the kth output node. The intermediate variables for the jth hidden layer node and the outputs of the neural networks can be written as:

$$p_{c,j}={\displaystyle \sum _{i=1}^n{e}_i(t)w_{\mathrm{critic}1,ij}(t)}+{\displaystyle \sum _{i=n+1}^{2n}{x}_{ref,i-n}(t)w_{\mathrm{critic}1,ij}(t)},\hspace{0.33em}\hspace{0.33em}j=1,\dots ,10$$

(23)

$$q_{c,j}=\tanh (p_{c,j})=\frac{e^{p_{c,j}}-e^{-p_{c,j}}}{e^{p_{c,j}}+e^{-p_{c,j}}},\hspace{0.33em}\hspace{0.33em}j=1,\dots ,10$$

(24)

$${\widehat{\lambda }}_k={\displaystyle \sum _{j=1}^{10}{q}_{c,j}{w}_{\mathrm{critic}2,jk}(t)},\hspace{0.33em}\hspace{0.33em}k=1,\dots ,2n$$

(25)

$$p_{a,j}={\displaystyle \sum _{i=1}^n{e}_i(t)w_{\mathrm{actor}1,ij}(t)}+{\displaystyle \sum _{i=n+1}^{2n}{x}_{ref,i-n}(t)w_{\mathrm{actor}1,ij}(t)},\hspace{0.33em}\hspace{0.33em}j=1,\dots ,10$$

(26)

$$q_{a,j}=\tanh (p_{a,j})=\frac{e^{p_{a,j}}-e^{-p_{a,j}}}{e^{p_{a,j}}+e^{-p_{a,j}}},\hspace{0.33em}\hspace{0.33em}j=1,\dots ,10$$

(27)

$$u_{a,k}={\displaystyle \sum _{j=1}^{10}{q}_{aj}{w}_{\mathrm{actor}2,jk}(t)},\hspace{0.33em}\hspace{0.33em}k=1,\dots ,m$$

(28)

It is desired to update the critic NN by minimizing the square Bellman error:

$$E_C(t)=\frac{1}{2}{e}_C^{\mathrm{T}}(t)e_C(t)$$

(29)

where

$$\begin{array}{l}{e}_C(t)=\frac{\partial \left[V_X(t-1)-J_X(X(t-1),u_a(t-1))-\gamma V_X(t)\right]}{\partial X(t-1)}=\lambda (X(t-1))-\frac{\partial J_X(X(t-1),u_a(t-1))}{\partial X(t-1)}-\gamma \frac{\partial V_X(t)}{\partial X(t-1)}\\ =\lambda (X(t-1))-\left(X{(t-1)}^{\mathrm{T}}{Q}_X+u_a{(t-1)}^{\mathrm{T}}R\frac{\partial u_a(t-1)}{\partial X(t-1)}\right)-\gamma \lambda (X(t))\left(\frac{\partial X(t)}{\partial X(t-1)}+\frac{\partial X(t)}{\partial u_a(t-1)}\frac{\partial u_a(t-1)}{\partial X(t-1)}\right)\end{array}$$

(30)

In Equation (30), $\lambda (X(t-1))$ and $\lambda (X(t))$ can be approximated by critic NN. $\partial u_a(t-1)/\partial X(t-1)$ is calculated by actor NN. The transition matrices $\partial X(t)/\partial X(t-1)$ and $\partial X(t)/\partial u_a(t-1)$ can be obtained from Equation (9). Because the partial derivative $f_1$ with respect to the angle rate of other channels is close to 0, we can obtain the transition matrices as follows: (The input dynamics are identified by RLS.)

$$\frac{\partial X(t)}{\partial X(t-1)}\approx {\left[\begin{array}{cc}\mathrm{I}& 0\\ 0& \mathrm{I}\end{array}\right]}_{2n\times 2n}$$

(31)

$$\frac{\partial X(t)}{\partial u_a(t-1)}\approx {\left[\begin{array}{c}\mathrm{\Delta }tg(t-1)\\ 0\end{array}\right]}_{2n\times m}\approx {\left[\begin{array}{c}{\widehat{\varsigma }}_1(t-1)\\ 0\end{array}\right]}_{2n\times m}$$

(32)

The weights of critic NN are updated by the gradient descent method with a learning rate ${\eta }_C$.

$$\begin{array}{l}{w}_{\mathrm{critic}2}(t+1)=w_{\mathrm{critic}2}(t)-{\eta }_C{e}_C{(t)}^{\mathrm{T}}\frac{\partial \widehat{\lambda }(X(t-1))}{\partial w_{\mathrm{critic}2}(t)}\\ \approx w_{\mathrm{critic}2}(t)-{\eta }_C{\left(\widehat{\lambda }(X(t-1))-\left(X{(t-1)}^{\mathrm{T}}{Q}_X+u_a{(t-1)}^{\mathrm{T}}R\frac{\partial u_a(t-1)}{\partial X(t-1)}\right)-\gamma \widehat{\lambda }(X(t))\left(\left[\begin{array}{cc}\mathrm{I}& 0\\ 0& \mathrm{I}\end{array}\right]+\left[\begin{array}{c}{\widehat{\varsigma }}_1(t-1)\\ 0\end{array}\right]\frac{\partial u_a(t-1)}{\partial X(t-1)}\right)\right)}^{\mathrm{T}}\frac{\partial \widehat{\lambda }(X(t-1))}{\partial w_{\mathrm{critic}2}(t)}\end{array}$$

(33)

$$\begin{array}{l}{w}_{\mathrm{critic}1}(t+1)=w_{\mathrm{critic}1}(t)-{\eta }_C{e}_C{(t)}^{\mathrm{T}}\frac{\partial \widehat{\lambda }(X(t-1))}{\partial w_{\mathrm{critic}1}(t)}\approx w_{\mathrm{critic}1}(t)-\\ {\eta }_C{\left(\widehat{\lambda }(X(t-1))-\left(X{(t-1)}^{\mathrm{T}}{Q}_X+u_a{(t-1)}^{\mathrm{T}}R\frac{\partial u_a(t-1)}{\partial X(t-1)}\right)-\gamma \widehat{\lambda }(X(t))\left(\left[\begin{array}{cc}\mathrm{I}& 0\\ 0& \mathrm{I}\end{array}\right]+\left[\begin{array}{c}{\widehat{\varsigma }}_1(t-1)\\ 0\end{array}\right]\frac{\partial u_a(t-1)}{\partial X(t-1)}\right)\right)}^{\mathrm{T}}\frac{\partial \widehat{\lambda }(X(t-1))}{\partial q_c}\frac{\partial q_c}{\partial p_c}\frac{\partial p_c}{\partial w_{\mathrm{critic}1}(t)}\end{array}$$

(34)

To minimize the value function in Equation (21), we can update the weights of actor NN by the gradient descent method, similarly. The learning framework of ED-DHP is shown in Algorithm 1.

$$\begin{array}{l}{w}_{\mathrm{actor}2}(t+1)=w_{\mathrm{actor}2}(t)-{\eta }_A\frac{\partial V_X(t-1)}{\partial u_a(t-1)}\frac{\partial u_a(t-1)}{\partial w_{\mathit{actor}2}(t)}\approx w_{\mathrm{actor}2}(t)-{\eta }_{\mathrm{A}}\left[\frac{\partial J(t-1)}{\partial u_a(t-1)}+\gamma \widehat{\lambda }(X(t))\frac{\partial X(t)}{\partial u_a(t-1)}\right]\frac{\partial u_a(t-1)}{\partial w_{\mathrm{actor}2}(t)}\\ \approx w_{\mathrm{actor}2}(t)-{\eta }_A\left[u_a{(t-1)}^{T}R+\gamma \widehat{\lambda }(X(t)){\left[\begin{array}{c}{\widehat{\varsigma }}_1(t-1)\\ 0\end{array}\right]}_{2n\times m}\right]\frac{\partial u_a(t-1)}{\partial w_{\mathrm{actor}2}(t)}\end{array}$$

(35)

$$\begin{array}{l}{w}_{\mathrm{actor}1}(t+1)=w_{\mathrm{actor}1}(t)-{\eta }_A\frac{\partial V_X(t-1)}{\partial u_a(t-1)}\frac{\partial u_a(t-1)}{\partial w_{\mathit{actor}1}(t)}\approx w_{\mathrm{actor}1}(t)-{\eta }_{\mathrm{A}}\left[\frac{\partial J(t-1)}{\partial u_a(t-1)}+\gamma \widehat{\lambda }(X(t))\frac{\partial X(t)}{\partial u_a(t-1)}\right]\frac{\partial u_a(t-1)}{\partial w_{\mathrm{actor}1}(t)}\\ \approx w_{\mathrm{actor}1}(t)-{\eta }_{\mathrm{A}}\left[u_a{(t-1)}^{\mathrm{T}}R+\gamma \widehat{\lambda }(X(t)){\left[\begin{array}{c}{\widehat{\varsigma }}_1(t-1)\\ 0\end{array}\right]}_{2n\times m}\right]\frac{\partial u_a(t-1)}{q_a}\frac{\partial q_a}{\partial p_a}\frac{\partial p_a}{\partial w_{\mathrm{actor}1}(t)}\end{array}$$

(36)

Algorithm 1 Learning Framework

Require:

Simulation parameters: sampling time $\mathrm{\Delta }t$ and the total simulation time $T$;

The controller parameters for the slow loop subsystem;

The forgetting factor of RLS;

The positive definite matrices;

The discount factor and the learning rates of two neural networks.

Initialize:

The initial states and inputs of the air vehicle;

The parameter matrix and the estimation covariance matrix of RLS;

Initialize the parameters of two neural networks randomly.

Online Learning:

1: for t = 0 to $(T/\mathrm{\Delta }t)-1$ do;

2:  Calculate reference angular rate vector ${\omega }_{ref}$ of the fast loop subsystem with (4);

3:  Obtain the angular rate vector $\omega$ through the sensors;

4:  Update the RLS model with (11)–(13);

5:  Update the parameters of the critic NN with (33) and (34);

6:  Update the parameters of the actor NN with (35) and (36);

7:  Calculate the output vector of the actor NN $u_a$ with (26)–(28);

8:  Calculate the control input vector of the air vehicle with (19);

9:  Get the next system states by taking action $u$;

10: end for.

4. Numerical Simulation and Discussion

In this section, the ED-DHP’s capability of online self-learning without complete knowledge of system dynamics and reference trajectory dynamics is demonstrated. Two numerical simulations are conducted to verify the design of the control scheme, sinusoidal signal tracking at a fixed operating point and command tracking with a morphing process at variable operating points.

4.1. Example 1

The scheme firstly controls the angle of attack $\alpha$, sideslip angle $\beta$, and bank angle ${\gamma }_V$ to track the sinusoidal reference signals under the condition of a fixed operating point. The initial conditions are $\alpha (0)=2^{\circ }$, $\beta (0)={\gamma }_V(0)=0^{\circ }$, and ${\omega }_x(0)={\omega }_y(0)={\omega }_z(0)=0^{\circ }/\mathrm{s}$. The primary control surfaces are ${\delta }_z(0)=-{11}^{\circ }$ and ${\delta }_x(0)={\delta }_y(0)=0^{\circ }$, which are close to the sketchy trim values. The reference commands are designed as:

$$\left\{\begin{array}{l}{\alpha }_{ref}=2+2\times \sin (0.4\pi t)\\ {\beta }_{ref}=1\times \sin (0.4\pi t)\\ {\gamma }_V{}_{ref}=1\times \sin (0.4\pi t)\end{array}\right.$$

(37)

The controller parameters for the slow loop subsystem are set as $K_{\alpha }=K_{\beta }=K_{{\gamma }_V}=3.5$. The positive parameters of TD are designed as $r=120$ and $\tau =0.9$. The forgetting factor of RLS is ${\gamma }_{\mathrm{RLS}}=1.0$, to ensure consistency of the identification. It was found that if ${\gamma }_{\mathrm{RLS}}$ is set to a value less than 1, the identification speed is improved with the reduction of stability of ED-DHP. So the parameters identified onlineshould have the correct sign and, to some extent, the correct relative magnitude in the self-learning process. The initial parameter matrix of RLS is set as follows, (which only conforms to the polarity of the vehicle):

$$\varsigma (0)=\left[\begin{array}{c}\begin{array}{ccc}-0.3& 0.0& 0.0\\ 0.0& -0.3& 0.0\\ 0.0& 0.0& -0.3\end{array}\\ \begin{array}{ccc}0.05& 0.05& 0.05\end{array}\end{array}\right]\in {ℝ}^{(m+1)\times n}{|}_{m=n=3}$$

The initial estimation covariance matrix is set as:

$$Cov(0)=\left[\begin{array}{c}\begin{array}{cccc}100& 0& 0& 0\end{array}\\ \begin{array}{cccc}0& 100& 0& 0\end{array}\\ \begin{array}{cccc}0& 0& 100& 0\end{array}\\ \begin{array}{cccc}0& 0& 0& 100\end{array}\end{array}\right]\in {ℝ}^{(m+1)\times (m+1)}{|}_{m=n=3}$$

The ED-DHP is implemented with the sampling time $\mathrm{\Delta }t=0.01s$, and the integration period of RK4 is 0.001s. The weights of the critic NN and actor NN are randomly initialized in $[-0.5,0.5]$. The positive definite matrices of ED-DHP are $Q=\mathrm{diag}\{800,800,800\}$ and $R=\mathrm{diag}\{0.1,0.2,0.2\}$. The learning rate of critic NN is ${\eta }_C=0.004$ and of the actor NN is ${\eta }_{\mathrm{A}}=0.002$, respectively. The discount factor $\gamma$ is set as 0.99.

As is shown in Figure 4, Figure 5, Figure 6 and Figure 7, both $\alpha$, $\beta$, and ${\gamma }_V$ can follow the reference signals after less than 2 s. From Figure 4, we see that during the learning process, the max dynamic tracking errors of $\alpha$, $\beta$, and ${\gamma }_V$ are 0.6 degrees, 1 degree, and 0.8 degrees, respectively.

After the parameters of the neural networks converge to a steady state, the tracking errors of the three angles are maintained within 0.05 degrees. The same observation is made from the parameters of two NNs, as illustrated in Figure 7. The ED-DHP shows good stability even though the learning rate does not change throughout the tracking process. However, it is worth noting that the angular rate loop has a large amplitude oscillation, especially the yaw rate loop and the roll rate loop. The phenomenon is mainly aroused by the cross-channel coupling and characteristics of the air vehicle’s roll channel.

4.2. Example 2

To further verify the effectiveness and robustness of the proposed control scheme, the guidance command tracking with morphing process at variable operating points is considered, as well as the model uncertainty.

The initial tracking errors of the three channels are 1 degree, and the control parameters of the fast loop and slow loop are all the same as in Section 4.1. Three cases with different parameters perturbation are carried into the numerical simulation, shown in

Table 1. The parameter uncertainty for 3 cases.

Case 0	Case 1	Case 2
$\mathrm{\Delta }I=0\%$	$\mathrm{\Delta }I=-5\%$	$\mathrm{\Delta }I=5\%$
$\mathrm{\Delta }m=0\%$	$\mathrm{\Delta }m=-5\%$	$\mathrm{\Delta }m=5\%$
$\mathrm{\Delta }\rho =0\%$	$\mathrm{\Delta }\rho =10\%$	$\mathrm{\Delta }\rho =-10\%$
$\mathrm{\Delta }{C}_m=0\%$	$\mathrm{\Delta }{C}_m=30\%$	$\mathrm{\Delta }{C}_m=-30\%$

. The switch of the air vehicle’s configuration occurs between 20 s~23 s. To ensure the stability of initial learning, TD starts to act on the system at 5 s.

As is shown in Figure 8, the whole duration of the simulation is 30 s, and the guidance commands are not concrete expressions. The controller makes the vehicle track the reference value of the attack angle of 2 degrees and maintains a sideslip angle and bank angle of around 0 degrees. After that, the vehicle is lifted at an attack angle of 8 degrees, morphing and maintained at the current angle of attack. Additionally, in this simulation condition, the flight speed of the vehicle is fast, and the dynamic pressure changes dramatically within the simulation time.

$\alpha$, $\beta$, and ${\gamma }_V$ can follow the reference signals after about 2 s. In Figure 9, the max tracking errors of $\alpha$, $\beta$, and ${\gamma }_V$ are both within 0.5 degrees before 2 s. The tracking errors of the attack angle are shown in

Table 2. Tracking errors of attack angle at different time periods.

Time	Case 0	Case 1	Case 2
2 s~10 s	$<{0.005}^{\circ }$	$<{0.01}^{\circ }$	$<{0.0025}^{\circ }$
10 s~20 s	$<{0.19}^{\circ }$	$<{0.22}^{\circ }$	$<{0.21}^{\circ }$
20 s~25 s	$<{0.08}^{\circ }$	$<{0.1}^{\circ }$	$<{0.06}^{\circ }$
25 s~30 s	$<{0.02}^{\circ }$	$<{0.005}^{\circ }$	$<{0.04}^{\circ }$

.

Figure 10 shows that the drastic changes in control surfaces are mainly concentrated in the first two seconds of the learning phase. In the lifting process, the symbol change of the pitch channel’s equivalent control surface has occurred, which means the desperate transformation of the operating points. This can be furthermore reflected in the sharp change of critic NN parameters in Figure 11, while the parameters of actor NN are only slightly adjusted.

The parameters of critic NN are obviously affected during the morphing process. However, the actor NN is adjusted slightly. It shows that the updating of NNs is mainly related to the tracking error. When actor NN owns strong control ability, that is, the tracking errors are very small during the morphing, the network parameters can adapt to the process only by fine-tuning.

5. Conclusions

This article presents the ED-DHP method for online nonlinear flight attitude control. ED-DHP defines the augmented system state, which contains tracking error and reference trajectory, to solve the tracking problem. This can make the actor neural network learn the feedforward and feedback control terms simultaneously without requiring complete knowledge of the vehicle dynamics and the reference trajectory dynamics. An RLS method is used to identify not only the input dynamics but also the total uncertainty based on the error dynamics in real-time.

This proposed method is applied in two numerical experiments based on a morphing air vehicle model: sinusoidal signal tracking at a fixed operating point and guidance command tracking with morphing process at variable operating points. The results indicate the robustness and adaptiveness of ED-DHP. They also reveal that ED-DHP can learn the multichannel attitude control policy online without offline training.

In the meantime, there are some questions worth pondering:

• Augmented state increases the dimension of the system state space, which may affect the convergence speed of the algorithm without offline training;
• The systems with input constraints can be further considered by designing the new cost function to make this method better apply to real vehicle systems;
• ED-DHP is suitable for the control of fast loop subsystems. It is necessary to optimize the error dynamic and RLS method to expand the application scenarios of ED-DHP.