Nonlinear suboptimal attitude control law for near-space hypersonic vehicle: Eigenvalue analysis and weight matrices design

This paper considers a nonlinear suboptimal control problem for a near-space hypersonic vehicle’s (NSHV’s) attitude dynamics. The least-square and stable manifold methods ﬁrst solve an unconstrained approximately optimal control law corresponding to the nonlinear attitude model. Then, to further meet the dynamic performance requirement of the attitude control system, a novel strategy based on the Koopman operator, symplectic geometric theory, and the stable manifold theorem is proposed to approximate the eigenvalues of the closed-loop nonlinear unconstrained approximated optimal control system. The weight matrices in the optimal performance index, which directly determine the output responses of the nonlinear attitude dynamics, can be appropriately designed according to the eigenvalues. The ﬁnal control law considers the actuator constraints. The NSHV’s closed-loop attitude control system is proved to be locally exponentially stable, and the suboptimality of the control law is analyzed. Numerical simulation demonstrates the eﬀectiveness of the proposed scheme.


Introduction
Near-space hypersonic vehicles (NSHVs) operate in the atmospheric region from about 20 km to 100 km height at speeds more than 5 Mach. They have been intensively researched for more than half a century. Their enormous technological advantages in both military and civil fields have gained widespread attention [23]. In recent years, the enhanced research on the reusable launch systems for the NSHVs shows the potential economic benefit of the related cutting-edge technology, which makes NSHVs an ongoing research focus [14].
Among the research items on the NSHVs, attitude control is one of the most challenging issues because NSHVs' attitude control systems are complexly nonlinear. Their longitudinal and lateral dynamics are strongly coupled. Hence, it is difficult to use linear models to represent the attitude dynamics, which makes nonlinear control law design necessary [24].
Until now, several nonlinear control techniques have been applied to NSHVs' attitude models, such as backstepping control [4,27], and sliding mode control [25,26]. In these methods, virtual control law design as the intermediate process is requisite, making the design process and the stability analysis more complex. Besides, when initial attitude errors are relatively large, the magnitudes of the control commands may become too large to be practical.
Relative to the methods mentioned above, nonlinear optimal control schemes have the following two advantages. Firstly, they do not adopt the cascade control frame, and avoid the virtual control laws. Secondly, the optimal performance index is introduced, which pursues using less control energy to achieve satisfactory output tracking performance. These benefits are significant for the NSHVs' attitude control because, in hypersonic flight, the attitude dynamics must respond quickly and accurately with reasonable energy consumption.
Familiar nonlinear optimal control strategies include the adaptive dynamic programming (ADP) and the pseudospectral method. Their features are summarized as follows.
ADP is established on the Hamilton-Jacobi-Bellman (HJB) equation. It utilizes neural networks to approximate the solution of the HJB partial differential equation. The dynamic systems are often required to be globally Lipschitz continuous [11], and the NSHVs' longitudinal dynamic systems satisfy this condition. Nevertheless, when the NSHVs' longitudinal and lateral maneuvers are both considered, the attitude control systems become locally Lipschitz continuous, and ADP is no longer suitable in this case. In addition, adjusting the initial weight values of neural networks is crucial to solving problems smoothly, but the related theoretical guidance is deficient [2], which limits ADP's practicality.
The pseudospectral method has been practically used in astronautics [18], particularly in guidance law design. It is based on Pontryagin's maximum principle, and transforms optimal control problems into nonlinear programming problems (NLPs). In this scheme, control systems must be locally Lipschitz continuous, which is suitable for NSHVs' attitude dynamics. However, NLPs are very sensitive to the initial guesses of covectors [3]. Since the nonlinear dynamic responses of NSHVs' attitude are complex and fast, it is hard to estimate the initial guesses of the covectors appropriately. Therefore, the pseudospectral method cannot be directly applied to the NSHVs' attitude models.
Besides the above two common methods, another nonlinear optimal control scheme proposed in [19] is based on the stable manifold theory [17]. It aims to solve infinite horizon optimal control problems based on the Hamiltonian-Jacobi equation (HJE). In this framework, the canonical equation of the HJE is decomposed into a stable direction and an anti-stable direction, then the stable manifold theory can be used to approximate the optimal solution. However, the stable manifold theory assumptions are too strict for most nonlinear control systems, which limits its application.
The nonlinear optimal control techniques mentioned above only approximate optimal control laws, but cannot solve the eigenvalues of closedloop dynamics. Nevertheless, eigenvalue analysis is essential for studying control systems' time response characteristics [7]. It can also be used to design the weight matrices in the optimal performance indexes [9], which directly determine the output responses of dynamic systems. Hence, eigenvalue analysis is especially necessary for the optimal control law design of the NSHVs' attitude systems which require rapid and accurate dynamic responses.
In this paper, to deal with the control problems mentioned above, we propose a novel suboptimal control scheme for a NSHV's attitude system as follows.
At first, we solve the nonlinear unconstrained approximated optimal control problem corresponding to the NSHV's attitude model. In order to avoid the effects of different initial guess values on solving optimization problems, the proposed scheme is based on the stable manifold theorem [17]. The key step is to use the least square method to separate the complex nonlinear NSHV's attitude system into linear and nonlinear parts. Since the attitude control system is locally Lipschitz continuous, it is demonstrated that the separation is practical. Then the canonical equation of the HJE corresponding to the unconstrained optimal dynamics can be decomposed into a stable direction and an anti-stable direction. Meanwhile, the decomposed canonical equation satisfies the assumptions of the stable manifold theorem, such that the iteration process of the theorem can approximate the nonlinear unconstrained optimal solution.
After solving the approximated unconstrained optimal control law, the nonlinear unconstrained optimal dynamics is mapped into a highdimensional nonlinear optimal dynamics. Based on the Koopman operator method [1] and the symplectic geometric theory [6,15], it is proved that the two dynamics are equivalent. When the state vector of the high-dimensional dynamics is selected appropriately, the nonlinear unconstrained optimal control system can be further approximated by a linear system, whose eigenvalues reveal the dynamic response characteristics of the nonlinear system. The approximation process is based on the Weierstrass approximation theorem [5] and the stable manifold theorem [17]. Then, the weight matrices in the optimal performance index can be designed according to the eigenvalues to meet the performance requirement of the NSHV's nonlinear attitude control system.
Finally, the actual control law considers the actuator constraints, which are seen as the perturbation items in the stability analysis. It is proved that the NSHV's closed-loop nonlinear attitude control system with the actuator constraints is locally exponentially stable, and the suboptimality of the control law is analyzed.
The content in this paper is organized as follows. In Section 2, the NSHV's attitude model and the corresponding unconstrained optimal control problem derived from the HJE are introduced. The form of the saturated control law considering the actuator constraints is derived. In Section 3, the optimal control law for the unconstrained dynamics is approximated. In Section 4, the eigenvalues of the closed-loop nonlinear unconstrained optimal control system are approximated, based on which the weight matrices in the optimal performance index can be appropriately designed. Then the NSHV's actuator constraints are considered, the locally exponentially stability of the closed-loop control system with the actuator constraints is proved, and the suboptimality of the saturated control law is analyzed. The numerical simulation in Section 5 verifies the effectiveness of the proposed scheme. Section 6 summarizes conclusions.

Problem description
The NSHV's attitude dynamics is modeled by [20] where α is the angle of attack, β is the sideslip angle, µ is the bank angle, p, q, and r are the roll, pitch, and yaw angular rates, respectively, M is the mass of the NSHV, S is the reference area, b and c are the reference lengths, I x , I y , I z , and I xz are the moments of inertia, V A is the airspeed,q is the dynamic pressure, g is the acceleration of gravity, C L,α , C Y,β , C l,β , C l,p , C l,r , C l,δa , C l,δr , C m,a , C m,q , C mδe , C n,β , C n,p , C n,r , C nδr , and C nδa are the aerodynamic coefficients which are the functions of Mach and α, δ r is the rudder, and δ e and δ a are the equivalent elevator and aileron, respectively, generated by the elevons δ L and δ R as follows, The elevons δ L and δ R and the rudder δ r are the actual control surfaces of the NSHV's attitude system.
To improve the performance of the NSHV's attitude control system, the actuator dynamics of δ L , δ R , and δ r are considered, based on which the magnitudes and rates of the control surfaces can be included in the performance index of the optimal control problem. The actuator dynamics are modeled by where δ Lv , δ Rv , and δ rv are the rates of the control surfaces, u ι (ι = 1, 2, 3) is the designed control surface command [28], the actual control surfaces δ L , δ R , and δ r are the responses of the actuator dynamics to u 1 , u 2 , and u 3 , respectively, O 3 is the 3-dimensional zero square matrix, I 3 is the 3dimensional identity matrix, A 1 = −ω 2 n I 3 , A 2 = −2ζω n I 3 , and B = ω 2 n I 3 , in which ω n is the natural frequency, and ζ is the damping ratio.
Besides, since u ι (ι = 1, 2, 3) is the control surface command, and the parameters of the actuator system (3) guarantees that the actual control surfaces δ L , δ R , and δ r can rapidly track reasonable u 1 , u 2 , and u 3 , respectively, the control surface command u ι (ι = 1, 2, 3) and the actual control surfaces δ L , δ R , and δ r satisfy the same constraints as follows where b lι < 0 and b uι > 0 (ι = 1, 2, 3) are the lower and upper bounds of the actual control surfaces δ L , δ R , and δ r , respectively.
In this paper, we design a suboptimal attitude tracking control law for (1)-(4) when the reference signals α c , β c , and µ c are step functions with respect to the time t. In addition, we aim to select the weight matrices in the optimal performance index according to eigenvalue analysis. Before solving these problems, (1)-(3) need to be transformed into an error dynamics as follows.
Corresponding to the actuator constraints (4), the error system (8) satisfies the constraints where e δ L , e δ R , and e δr are the components of e 3 , δ Lc , δ Rc , and δ rc are the components of x 3c in (6b), and b lι and b uι (ι = 1, 2, 3) are defined in (4).
To solve u e in (8), the following unconstrained optimal control problem is first considered. The unconstrained dynamics is defined aṡ where f and g are the same as in (8), u * e is the control law which does not consider the constraints (9), and e * is the state vector of the unconstrained system. The optimal performance index for (10a) is where Q is the semi-positive definite matrix considering the tracking error, the magnitudes of the angular rates, and the magnitudes and rates of the control surfaces, and R is the positive definite weight matrix considering the designed control surface commands. At each time t, (10b) is with regard to the current known x 1c . Since f (e * ) in (10a) satisfies f (0) = 0, the Hamilton Jacobi Equation (HJE) [12] corresponding to (10a)-(10b) is such that the optimal feedback control law u * e (t) for the unconstrained dynamics (10a) is Because the unconstrained system (10a) is complexly nonlinear, ρ(e * (t)) in (10) cannot be solved analytically and needs to be approximated. Denote the approximated solution of ρ(e * (t)) bŷ ρ(e * (t)), then the approximated optimal control law for the unconstrained optimal control problem (10) isû * Replace e * (t) in (11) with e(t), then the control law u e (t) for the original constrainted system (8)- (9) can be designed as u e (e(t)) = sat b (û * e (e(t))) (12a) where sat b (·) is the saturation function, which means the ιth (ι = 1, 2, 3) component of u e (e(t)) denoted by u eι (e(t)) is is the ιth component of u c which is solved in (6d), and b uι and b lι are the bounds defined in (4).
Then consider the actuator constraint (4), the NSHV's attitude dynamics (5) can be represented by where u c and u e are as in (6d) and (12), respectively, and the ιth (ι = 1, 2, 3) component of in which x 3ι and x 4ι are the ιth components of x 3 and x 4 , respectively, and b uι and b lι (ι = 1, 2, 3) are defined in (4).

Approximated optimal control law for unconstrained dynamics
In this section, we approximate the solution of the unconstrained optimal control problem (10) corresponding to the NSHV's attitude model. Our contribution is to use the least-square method to transform the unconstrained dynamics (10a) and (10d) into linear and nonlinear parts, such that the stable manifold method [17,19] can be applied to approximate ρ(e * (t)) in (10d).

Separating unconstrained dynamics into linear and nonlinear parts
In order to apply the stable manifold method [17,19] to approximate the solution of the unconstrained optimal control problem (10), the unconstrained dynamics (10a) and (10d) needs to be firstly separated into linear and nonlinear parts as follows.
where A ∈ R 12×12 , A and ρ(e * (t)) need to be solved, and N (e * ) is given by When the system (10) satisfies the following assumption, the least squared method can be used to realize the separation.
Assumption 1 For the unconstrained optimal control problem (10), when the initial error e * (t 0χ ) ∈ E 0 , where E 0 is a compact set, the approximated optimal control lawû * e in (11) exists, and makes e * (t) ∈ E for all t ∈ [t 0χ , t f χ ) (16) where E is a compact set, and E 0 ⊂ E.
Based on Assumption 1, the nonlinear vector field in the unconstrained dynamics (10a) can be separated into linear and nonlinear parts by the following steps.
Step 1. Discretize (17) to obtain the discrete vector field where ∆T is the small enough sampling time, and f d can be obtained by the Runge-Kutta method. When e * (k) ∈ E, according to the Weierstrass approximation theorem [5], f d (e * (k), ∆T ) can be uniformly approximated by polynomials, which means that for given ∆T , e * (f ) (k + 1) can be represented by where A d ∈ R 12×12 needs to be solved, and N d (e * (k)) is Step 2. To solve A d and N d (e * (k)) in (18b) over the compact set E, collect the data set as follows where e * l ∈ E (l = 1, 2, · · · , L), e * + l is generated by , and L is large enough such that E can sufficiently cover E. According to (18), the data set (19) can be represented by Step 3. Define the optimal performance index According to the least squares method, where [·] T is the transpose of matrix, and [·] ‡ is the pseudoinverse. Since the row vectors of E is linearly independent according to the nonlinear discrete vector field (18), A d in (21c) is the unique solution of (21a)-(21b) [13]. As a result, for any e * (k) = e * l ∈ E (l = 1, 2, · · · , L), e * (f ) (k + 1) in (18) satisfies where N d (e * l ) is given by in which b N is determined by (21). Then consider the case that e * (k) ∈ E and e * (k) / ∈ E. In this case, the bound of the norm of N d (e * (k)) in (18b) is analyzed as follows. Let Since f d is locally Lipschitz continuous over the domain E, where e * (f ) (k + 1) is as in (18a), l 1 is the Lipschitz constant. Besides, sinceê * ∈ E, according to (22) In addition, the inequality holds over the domain E. As a result, according to (18b) and (24), the bound of the norm of N d (e * ) for any e * (k) ∈ E and e * (k) / ∈ E satisfies Step 4. Apply the bilinear transformation [16] to (18b) to obtain the continuous systeṁ where A ∈ R 12×12 is solved by the bilinear transformation [16], and N (e * ) is as in (15b), which has a proper bound over the compact set E. Substitute the vector field (26) and the control law (10d) into the unconstrained dynamics (10a), then the unconstrained control system (10a) and (10d) can be represented by (15).

Approximated optimal control law based on stable manifold method for unconstrained dynamics
The stable manifold method [17,19] can be applied to (15) to approximate ρ(e * (t)). The detailed process is summarized in the following steps, in which Steps 1-5 are the same as in [19].
Step 1. According to (15), the HJE (10c)-(10d) can be represented by Step 2. The canonical equation of the HJE (27) is and N x (e * ) and N ρ (e * , ρ) are the nonlinear terms Step 3. Solve the Riccati equation where P is a symmetric matrix, and is the unique solution of (29a).
Step 4. Solve the Lyapunov equation where V is a symmetric matrix, and is the solution of (29b).
Step 5. Let where I is the identity matrix. Define the transformation then the dynamics of e ′ and ρ ′ are The transformation (30) decouples the linear part of (28a) into stable and anti-stable directions, which correspond to U and −U T in (30b), respectively. Then the system (30b) satisfies the following assumptions [17].
Based on the above assumptions, the following iterative integration in the stable manifold theorem [17] can approximate ρ ′ in (30b). Let Then, when k = 0, 1, 2, · · · , the solution of e ′ k+1 (t) and ρ ′ k+1 (t) are given by where ρ ′ k (∞) = 0. In practice, ∞ in the integral interval of (33b) can be replaced with a large enough positive real number t f .
4 Eigenvalue analysis, weight matrices design, and system performance analysis In this section, we aim to design an eigenvalue analysis scheme for the closed-loop nonlinear unconstrained optimal dynamics (10a) and (10d) corresponding to the NSHV's attitude model, based on which we can design the weight matrices in the optimal performance index (10b) to meet the NSHV's dynamic response performance requirement. Then the final NSHV's attitude control law considers the actuator saturation, and the corresponding stability and suboptimality are analyzed.
The eigenvalue analysis and weight matrices design procedure is summarized as the following steps.
Step 1. Based on the Koopman operator method [1], the unconstrained optimal control system (10) is mapped into a high-dimensional optimal control system. According to the symplectic geometry [6,15], it is proved that the two systems are equivalent.
Step 2. Based on the Weierstrass approximation theorem [5] and the stable manifold theorem, the high-dimensional optimal control system can be approximated by a linear system, whose eigenvalues reveal the dynamic response characteristic of the original system (10).
Step 3. Based on the approximated eigenvalues obtained in Step 2, the weight matrices Q and R in (10b) can be appropriately designed.

Equivalent unconstrained optimal control systems based on symplectic geometry
In this section, the Koopman operator [1] is used to map the unconstrained optimal control system (10) into a high-dimensional unconstrained optimal control system. Then the symplectic geometric theory is used to prove that the two systems are equivalent. According to the Koopman operator theory, the high-dimensional system is constructed as follows. Let Ψ (e * ) = e * T ψ 1 (e * ) ψ 2 (e * ) · · · ψN (e * ) T (35a) where ψ j (e * ) (j = 1, 2, · · · ,N ) is the monomial in e * , Ψ (e * ) is the observable vector, whose dimension is N =N + 6, and C is the output matrix. Define the high-dimensional dynamicṡ where f (e * ) and g are the same as in (10a), and u * he is the control input. The optimal performance index for (36a) is given by The HJE corresponding to (36a)-(36b) is where λ is the covector. The optimal control law u * he (t) for (36a)-(36b) is The following Theorem 1 is based on the symplectic geometry. It proves that the optimal control systems (36) and (10) are equivalent, which means that the unconstrained optimal control system (10) is equivalent to Theorem 1 Consider the optimal control systems (36) and (10), and suppose the transformation between the vector spaces (Ψ , λ) and (e * , ρ) is given by φ : then φ is symplectic, which means the hypothesis (37a) holds, such that and (36) is equivalent to (10).
Proof For any e * ∈ E, where E is the compact set defined in Assumption 1, Ψ (e * ) in (35a) satisfies Ψ ∈ E h , where E h is a compact set determined by E and the vector value function Ψ (e * ). Consider the diffeomorphism, h : which is a natural diffeomorphism [6], and satisfies where the calculation of ·, · is the same as the dot product of two vectors,ė * corresponds to (10), andΨ corresponds to (36). Then according to (39) and the definition 6.3.1 in [15], (38) defines a cotangent lift, such that the transformation φ in (37a) is a symplectic transformation. Hence, according to Proposition 2.4.2 [15], where · is the dot product, • denotes the function composition, X H is the Hamiltonian vector field corresponding to H in (10c), Dφ is the differentiation of φ with respect to Ψ and λ, and X H•φ is the Hamiltonian vector field corresponding to H • φ.
Besides, according to (36), (10) and (38), Substituting (40b) into (40a) yields Since the vector field X H is the first half elements of X H is the same asė * in (10). On the other hand, since the differential Dφ corresponding to the map φ in (37a) is

Approximated eigenvalues of unconstrained optimal control system
In this section, we aim to approximate the eigenvalues of the equivalent systems (36) and (10), and apply the eigenvalue analysis to reveal the time response characteristics of the nonlinear dynamics (10). We first separate the high-dimensional nonlinear dynamics (36) into linear and nonlinear parts, then based on the stable manifold method, the eigenvalues of (36) and (10) can be approximated. The first step is to use the least-square method in analogy to Section 3.1 to separate the vector field in (36a) into the linear and nonlinear parts as followṡ where A Ψ is obtained by the least-square method. The control gain matrix in (36a) is denoted by g Ψ , where I is the 12-dimensional identity matrix. Then according to A Ψ , g Ψ in (43), and the optimal performance index (36b), the Riccati equation for (36) is constructed as where R Ψ is given by in which ON ×3 is theN -by-3 zero matrix, and P Ψ is the unique solution of (44a). As a result, the dynamics of Ψ along the approximated optimal dynamics (10)-(11) can be represented bẏ The system (45) satisfies the following assumptions [17].
Consider k = 0 in (48b). According to the Weierstrass approximation theorem [5], the components of N Ψ (Ψ 0 ) along the solution Ψ 0 (t) in (48a) can be approximated by the polynomial functions of e * . Therefore, when the components of Ψ in (35a) are properly chosen, N Ψ (Ψ 0 ) can be represented by where M 0 ∈ R N ×N is the linear map, which can be solved by the least-squares method, and r 0 (s) is the approximation error. Then, substituting (49) into (48b) yields where the system matrix can be diagonalized into the diagonal matrix as follows, where Λ Ψ is the diagonal matrix whose components are the eigenvalues of F Ψ , and the block matrices satisfy the following equation.
where I is the identity matrix. Let then according to (51a) and (50), with the initial value The solution of (51b) is The initial value vector Υ 1 (t 0χ ) can be represented by where Λ r1 is the diagonal matrix, whose diagonal elements are the ratios between the components of Υ 1 (t 0χ ) and Υ 0 (t 0χ ). Substituting (51f) into (51d) yields Then according to (51a), (50c), and (51g), the solution Ψ 1 (t) can be represented by which means Ψ 1 (t) can be approximated by Ψ 0 (t), such that when k = 1, the terms N Ψ (Ψ 1 ) can also be approximated by Ψ 0 (t) as follows where M 1 ∈ R N ×N is the linear map, which can be solved by the least-squares method, and r 1 (s) is the residual error. The above analysis procedure (49)-(53) can be repeatedly applied to the iteration process in (48b) when k = 1, 2, · · · . The nonlinear term N Ψ (Ψ k (s)) can be represented by where M k ∈ R N ×N is the linear map, and r k (s) is the residual error. M k and r k (s) are solved in the iteration process. Furthermore, the solution of Ψ k (t) can be represented by where Λ rk is obtained in the above iteration analysis process, and r ik (t) is According to the Weierstrass approximation theorem [5], when the components of Ψ (e) in (35a) are properly chosen, r ik (t) in (55) can be small. Besides, according to the stable manifold theory, the approximated optimal solution in (45) can be represented by where r Ψ (t) is the approximated error, when k is large enough, r Ψ (t) is small. Then, let according to (55)-(57), (48a), and (37a), the unconstrained optimal dynamics e * in (10) can be represented by = C e Ψ 0 (e * (t)) + r e (t) (58a) where C e = CM Ψ is the output matrix, and r e (t) = C(T Ψ r ik (t) + r Ψ (t)) is the residual error, when the components of Ψ (e * ) in (35a) are properly chosen, and k is large enough, r e (t) is small. According to the linear system (48a) and (58a), the eigenvalues of the system matrix F Ψ in (48a) can approximate the eigenvalues of the unconstrained optimal dynamics (10), and reveal its time response characteristics. The whole eigenvalue approximation process is summarized in Figure 1. When the actuator saturation does not occur over the time interval [t 0χ , t fχ ), (58a) can be used to interpret the dynamic response of (14). The corresponding approximated optimal solution of (14) is given by where C Ψ = CM Ψ T Ψ is the output matrix. Then the κth (κ = 1, 2, · · · , 12) element of e(t) denoted by e κ can be approximated by where ν κ,j is the element in the κth row and jth column of the matrix C Ψ , γ j is the jth diagonal element of Λ Ψ , ̟ κ is the κth row vector of the matrix T −1 Ψ , ξ κ,j = ν κ,j {̟ κ Ψ(e(t 0χ ))}, and r eκ (t) is the κth element of r e (t) in (58a).
Then the eigenvalue analysis for the NSHV's closed-loop approximated optimal attitude dynamics (14) without actuator saturation is similarly to [8], which can be summarized as follows.
3. The exponential item exp(γ j (t − t χ0 )) determines the nature of the time response of (14).
4. The scalar ν κ,j represents the extent to which each state contributes to the jth mode of (14), which is determined by the initial value Ψ(e(t 0χ )) and the nonlinear term N Ψ (Ψ ) in (45a) according to the solving process for (58).

Weight matrices design and system performance analysis
In this section, we first design Q and R in (10b) according to the eigenvalues analysis for (58c). Then we consider the NSHV's attitude control system with actuator saturation, and analyze the corresponding closed-loop system stability and the suboptimality of the saturated control law.
The unconstrained optimal dynamics (see(10a) and (10d)) The equivalent high-dimensional optimal dynamics (see (36e)) Separating the high-dimensional optimal dynamics into the linear and nonlinear parts (see (45)) Approximating the eigenvalues of the unconstrained optimal dynamics (see(58a)) The design process for Q and R is similar to the linear-quadratic-regulator (LQR) method [9], which is summarized as follows.
Step 1. Let Q andR be diagonal matrices [9] as follows where Q considers the time responses of e α , e β , and e µ , e α max , e β max , and e µ max are the expected maximum values of e α , e β , and e µ , respectively, q 1 , q 2 , and q 3 are the weight values, which satisfies q 1 + q 2 + q 3 = 1, u e1 max , u e2 max , and u e3 max are the expected maximum values of u e1 , u e2 , and u e3 , respectively, r 1 , r 2 , and r 3 are the weight values, which satisfies r 1 + r 2 + r 3 = 1.
Step 2. Let where c r is a constant which needs to be designed [9].
Step 3. Set c r at different values, solve F Ψ according to (44)-(45), and calculate the eigenvalues of F Ψ . Then according to the eigenvalue analysis, an appropriate value of c r can be selected to satisfy the dynamic response requirement of the NSHV's attitude dynamics.
Finally, consider the NSHV's attitude error dynamics (14a) with the actuator constraints (14b). According to the selected Q and R, the control law with the actuator constraints is designed as in (12).
The following theorem proves the closed-loop stability of the NSHV's attitude error control system.

Theorem 2
The NSHV's attitude error control system (14) under the control law (12) is locally exponentially stable.
Proof Consider the unconstrained approximated optimal dynamics (10), since the approximated optimal control law of (10) has been accurately approximated by (11), the optimal performance index (10b) can be fit into the following cost-to-go function [12,21] where V (e * (t)) is the polynomial function of e * (t), which is a Lyapunov function, and can be used to prove the locally exponential stability of the unconstrained optimal control system (10) [21]. Then replace e * (t) in the unconstrained optimal dynamics (10) and the Lyapunov function in (60) with e(t). Compare the constrained dynamics (14) and the unconstrained dynamics (10), let where O 1×3 is the 3-dimensional zero vector, and ∆ e3 is given by ∆ e3 = he(e 4 ) − e 4 (61b) and compare the constrained control law (12) and the unconstrained control law (11), let ∆u = sat b (û * e (e(t))) −û * e (t) (61c) then the NSHV's attitude error dynamics (14) under the control law (12) can be represented bẏ e = f (e) + gû * e + ∆e + g∆u (61d) which can be seen as the perturbed system relative to the unconstrained optimal dynamics (10). Therefore, the locally exponential stability of the closed-loop system (14) under the control law (12) can be proved according to the perturbed system theory [10].
The suboptimality of the saturated control law (12) can be analyzed as follows [22]. Case 1. If the control law in (12) satisfies u e (e(t)) =û * e (e(t)) over the time interval [t 0χ , t fχ ), then u e (e(t)) is the approximated optimal control law in [t 0χ , t fχ ).
Case 2. If the last actuator saturation begins at t s ∈ [t 0χ , t fχ ), then u e (e(t)) is approximately optimal in the interval [t s , t fχ ).

Numerical simulation
In this section, the numerical simulation verifies the effectiveness of the proposed schemes in this paper, including the eigenvalue analysis, weight matrices design, and the NSHV's attitude response under the actuator saturation.
The parameters mentioned above are in the unit of degree. In the simulation, the unit is transformed from the degree into the rad.
According to the data in Tables 3(a), when c r = 0.7, the error response of attack of angle e α (t) = α(t) − α c can be accurately approximated byê αI as followŝ e αI (t) = According to (63a)-(63b), for c r = 0.7 and c r = 0.2, the time responses of e α (t) can be approximated by the 5th order system and the 4th order system, respectively. Besides, according to the decay rates of components in (63), the components which determine the overshooting of e α (t) can be selected as follows. When c r = 0.7, the overshooting of e α (t) is determined byê o1 (t) as followŝ e o1 (t) = which consists of one decaying exponential damped sinusoidal response. Figure 3 demonstrates that the components in (64a) and (64b) can accurately approximate the overshooting of the time response of e α (t) when c r = 0.7 and c r = 0.2, respectively.
According to the data in Tables 3(a), when c r = 0.7, the error response of roll angle e µ (t) = µ(t)−µ c can be accurately approximated byê µI (t) as followŝ e µI (t) = According to (65), the time response e µ (t) can be approximated by second-order systems when c r = 0.7 and c r = 0.2.
Based on (63) and (65), the time response performance of e α (t) and e µ (t) can be calculated, and are listed in Table 4. The data in Table 4 illustrate that, when c r = 0.2, both e α (t) and e µ (t) have more rapid response speed than c r = 0.7, and when c r = 0.2, the overshoot of e α (t) is smaller than c r = 0.7. In conclusion, to improve the output response performance of the NSHV's attitude control system, c r = 0.2 is more proper than c r = 0.7.
In Figures 4(a), the overshoot of α(t) is 6.2649%, the rising time of α(t) is 3.3892s, and the settling time is 10.0386s. In Figures 4(b), the magnitude of β(t) is less than 0.02 • . In Figures 4(c), the rising time of µ(t) is 2.8740s, the settling time is 4.8464s, and there is no overshoot. Figure 4(d) illustrates the time responses of the angular rates. The magnitudes of angular rates p, q and r are less than 1.9619 • /s, 1.3695 • /s 0.3071 • /s, and they converge to the components of x 2c in (62a), respectively. Figure 4(e) shows the time responses of the control surfaces. Their peaks occurs at the initial stage, and converge to the components of x 3c in (62b), respectively. Figure 5 demonstrates the Hamiltonian value H in the HJE (10c). When the actuator saturation does not occur, the closer H is to zero, the better the suboptimality. In Figure 5, the value of H decays to zero quickly after the saturation of δ L , which means the suboptimality is good.
The above numerical simulation demonstrates that the suboptimal control law (13a), (6d) and (12) can effectively make the NSHV's attitude control model (1)-(4) track a step reference signal. Besides, the eigenvalue analysis method proposed in this paper can correctly reveals the time response characteristics of the closed-loop nonlinear NSHV's attitude dynamics when the actuator saturation does not occur, and help to design the weight matrices in (10b) for the NSHV's suboptimal attitude control problem. Sideslip angle β(t) c Bank angle µ(t) d Angular rates p(t), q(t) and r(t) e Control surfaces δ L (t), δ R (t) and δ r (t)