3.1 Gap metric

Let $K$ be a linear operator in Hilbert space $H$ [Reference French33]:

(5) \begin{align} K \,:\,H \to H. \end{align}

$H \times H$ is the set of all pairs $u,v$ ; $u$ and $v$ are in $H$ . $H \times H$ is also a Hilbert space with the inner product derived from $H$ . The domain of $K$ , $D\left( K \right)$ , is the set of all $u$ in $H$ with the property that $Ku$ is in $H$ . The graph $G\left( K \right)$ of $K$ is the set of all pairs $u,Ku$ with $u$ in $D\left( k \right)$ and is defined as $G\left( K \right) = \{ \left( {u,Ku} \right) \in H \times H|u \in D\left( K \right)\} $ .

Because $K$ is linear, $G\left( k \right)$ is a subspace of the product Hilbert space $H \times H$ . $K$ is considered closed if its graph $\left( K \right)$ is a closed subspace of $H \times H$ .

(6) \begin{align} \delta \left( {{K_1},{K_2}} \right) = \delta \left( {G\left( {{K_1}} \right),G\left( {{K_2}} \right)} \right) = {\rm{max}}\left( {{\delta _{12}}\left( {G\left( {{K_1}} \right),G\left( {{K_2}} \right)} \right),{\delta _{21}}\left( {G\left( {{K_2}} \right),G\left( {{K_1}} \right)} \right)} \right)\!, \end{align}

where ${\delta _{12}}\left( {G\left( {{K_1}} \right),G\left( {{K_2}} \right)} \right) = \mathop {{\rm{sup}}}\limits_{u \in D\left( {{K_1}} \right),u \ne 0} \mathop {{\rm{inf}}}\limits_{v \in D\left( {{K_2}} \right)} \frac{{u - {v^2} + {K_1}u - {K_2}{v^2}}}{{\sqrt {{u^2} + {K_1}{u^2}} }}.$

${\delta _{21}}\left( {G\left( {{K_2}} \right),G\left( {{K_1}} \right)} \right)$ is defined similarly.

Example 1: ${K_1} = \left[ {{\rm{cos}}\theta, - {\rm{sin}}\theta ;\,{\rm{sin}}\theta, {\rm{cos}}\theta } \right]$ is a rotation by $\theta $ radians. ${K_2} = \left[ {\lambda, 0;\,0,\lambda } \right]$ is a scaling by a factor of $\lambda $ . $u$ and $v$ are the unit vectors on the unit circle ${\mathbb{R}^2}$ . Figure 2 shows a gap $\delta \left( {{K_1},{K_2}} \right)$ .

Figure 2. Geometric representation of operators ${K_1}$ and ${K_2}$ with gap.

We consider linear systems as mappings in Hilbert space. For a multiinput multioutput linear system, which is represented by state matrices $A$ and $B$ and input matrices $C$ and $D$ , the transfer function is expressed as $P\left( s \right) = C{\left( {sI - A} \right)^{ - 1}}B + D$ , where $I$ is the identity matrix. Notably, $P\left( s \right)$ is a linear operation in Hilbert space.

(7) \begin{align} {\delta _{gap}}\left( {{P_1},{P_2}} \right) = {\rm{max}}\left( {{\delta _{12}}\left( {{P_1},{P_2}} \right),{\delta _{21}}\left( {{P_2},{P_1}} \right)} \right)\!, \end{align}

where

\begin{align*} {\delta _{12}}\left( {{P_1},{P_2}} \right) = \delta \left( {G\left( {{P_1}} \right),G\left( {{P_2}} \right)} \right) = {{\rm{\Pi }}_{G\left( {{P_1}} \right)}},{{\rm{\Pi }}_{G\left( {{P_2}} \right)}} = \mathop {{\rm{inf}}}\limits_{Q \in {H_\infty }} \left[ {\begin{array}{*{20}{l}}{{M_1}}\\[3pt] {{N_1}}\end{array}} \right] - \left[ {\begin{array}{*{20}{l}}{{M_2}}\\[3pt] {{N_2}}\end{array}} \right]{Q_\infty }, \end{align*}

and $\left( {{M_1},{N_1}} \right)$ , $\left( {{M_2},{N_2}} \right)$ represent the left prime decomposition of the linear transfer functions ${P_1}$ and ${P_2}$ , respectively. $Q$ is a Hilbert matrix.

The properties of the gap metric include the following (refer to Ref. (Reference French33,Reference Tao, Li and Li36) for more details):

1. (1) The gap metric is a measure of the distance between two linear systems.

2. (2) The gap metric lies between 0 and 1, $0 \le {\delta _{gap}}\left( {{P_1},{P_2}} \right) \le 1$ . The smaller the gap between the two systems, the closer their dynamic characteristics are.

3.2 Gap metric-based NPS algorithm for LPV modeling

The GM was used to analyse the generalised stability of parametric matrix families, which requires state space models. Therefore, it is crucial to develop an efficient LPV to cover nonlinear dynamics.

A nonlinear model in the LPV form can be represented using three methods: Jacobian linearisation, state transformations and function substitution [Reference Cavanini, Ippoliti and Camacho37]. Among them, Jacobian linearisation is widely used. However, it relies on a set of linearised models corresponding to several equilibrium points, which must cover the entire flight envelope to ensure good fidelity with the original nonlinear model. Therefore, to avoid choosing similar equilibrium points and reduce unnecessary computational complexity, the gap metric is introduced to measure the similarity among linear systems based on Jacobian linearisation, which is called the NPS algorithm. The NPS algorithm employs gap metric analysis to identify optimal nominal points within the sub-envelopes, ensuring accurate LPV modeling of nonlinear systems by minimising computational complexity and enhancing model fidelity.

First, the flight envelope is divided into subenvelopes according to the gap metric value. Then, a point in each subenvelope with the smallest difference in model characteristics with other state points is selected. This point is called the nominal point. Finally, the LPV model is obtained by fitting the linear model at nominal points using the following steps:

Suppose the large flight envelope is denoted as $\phi $ , and its range is $\left[ {{M_{{\rm{min}}}},{M_{{\rm{max}}}}} \right] \times \left[ {{h_{{\rm{min}}}},{h_{{\rm{max}}}}} \right]$ .

Step 1. In the flight envelope, the average grid points are divided at certain intervals of ${\Delta}M$ and ${\Delta}h$ , and a total of ${N_{Ma}} \times {N_h}$ state points are obtained.

(8) \begin{align} \left\{ {\begin{array}{*{20}{l}}{{N_{Ma}} = \left( {{M_{{\rm{max}}}} - {M_{{\rm{min}}}}} \right)/{\Delta}M + 1}\\[3pt] {{N_h} = \left( {{h_{{\rm{max}}}} - {h_{{\rm{min}}}}} \right)/{\Delta}h + 1} \end{array}} \right.\!. \end{align}

Step 2. The corresponding LTI set is obtained using the Taylor series expansion method at the grid points, which can be denoted as ${P_{ij}},i = 1,2, \ldots {N_{Ma}}$ , $j = 1,2, \ldots {N_h}$ .

Step 3. Calculate the gap measure between adjacent points using the expression ${\delta _{gap}}\left( {{P_{i,j}},{P_{i,j + 1}}} \right)$ or ${\delta _{gap}}\left( {{P_{i,j}},{P_{i + 1,j}}} \right)$ .

Step 4. Divide the envolpe into subenvelopes.

The altitude interval is divided only at ${M_{{\rm{min}}}}$ to avoid divergence between the Mach and altitude directions. When the cumulative value of the gap metric reaches a threshold $\tau $ , expressed as

(9) \begin{align} \mathop \sum \nolimits_{i = {i_0}}^k {\delta _{gap}}\left( {{P_{i,j}},{P_{i + 1,j}}} \right) \ge \tau, j = 1\left( {M = {M_{{\rm{min}}}}} \right)\!,\end{align}

the corresponding interval numbers can be considered the altitude boundaries of the subenvelope to divide the entire envelope, denoted by the set ${h_B}$ . Here, ${h_B} = \left\{ {{h_{B1}},{h_{B2}}, \ldots, {h_{Bn}}} \right\}$ , where each ${h_{Bi}}$ is a boundary value in the altitude direction. ${i_0}$ is the starting index for the summation, and $k$ is the number of intervals required for the cumulative gap metric to reach the threshold $\tau $ . After reaching the threshold $\tau $ , the next accumulation starts at ${i_0} = k + 1$ .

The speed interval is then divided at the altitude boundaries in the set ${h_B}$ . Similarly, if

(10) \begin{align} \mathop \sum \nolimits_{j = {j_0}}^k {\delta _{gap}}\left( {{P_{i,j}},{P_{i,j + 1}}} \right) \ge \tau, i \in \left\{ {i\left| {i = 1,i{\Delta}h \in {h_B}} \right.} \right\}\!,\end{align}

the corresponding interval numbers can be considered the Mach boundaries. Here, ${j_0}$ is the starting index of the summation. After reaching the threshold $\tau $ , the next accumulation starts at $k + 1$ . For each altitude boundary, the Mach values are divided into multiple intervals. They can be considered as the boundaries under the corresponding altitudes, denoted by the sets ${M_{B1}}$ ,…, ${M_{Bn}}$ . The subintervals are then numbered in ascending order of altitude and Mach. Thus, $\phi $ can be divided into $L$ subenvelopes based on these boundaries, denoted as ${\phi _l},l = 1,2, \ldots, L$ .

This process ensures a systematic division of the envelope $\phi $ into subenvelopes with multiple Mach intervals under each altitude boundary, thereby providing detailed and flexible segmentation, as shown in Fig. 3.

Figure 3. Division of the envelope into subenvelopes.

Step 5. Select nominal points.

The points within the subenvelope ${\phi _l}$ that is least different from the other point linear models, called the nominal point ${E_l}$ , are calculated. Notably, the nominal points must satisfy Equation (11)

(11) \begin{align} {\bar \delta _{gap}}\left( {{E_l}} \right) = \mathop {{\rm{min}}}\limits_{a = 1,2,...,{n_M} \times {n_h}} \left\{ {{1 \over {{n_M} \times {n_h}}}\left[ {\mathop \sum \limits_{b = 1,b \ne a}^{{n_M} \times {n_h}} {\delta _{gap}}\left( {{E_a},{E_b}} \right)} \right]} \right\}\!, \end{align}

where ${n_M}$ and ${n_h}$ are the numbers of state points in ${\phi _l}$ along the velocity and altitude dimensions, respectively, and $a$ and $b$ index all state points in the ${\phi _l}$ . ${E_a}$ is set to the nominal point being considered, and ${E_b}$ represents all other state points within the subenvelope, except the nominal point.

This ensures that the selected nominal point ${E_l}$ best represents the subenvelope ${\phi _l}$ , resulting in an accurate LPV model.

Step 6. Obtain the LPV model of the aircraft by fitting the LTI system corresponding to the nominal points.

Figure 4. Variation of gap metrics with Mach numbers at various altitudes.

Figure 5. Variation of gap metrics with altitudes at various Mach numbers.

In this study, the NPS algorithm was used to derive the LPV model. This is achieved by linearising a set of LTI systems covering the region of interest, resulting in an LPV model that locally approximates the original nonlinear system around specific equilibrium points.

3.3 LPV model

The large flight envelope considered herein is $h = \left[ {24,34} \right]$ km, ${\rm{\;\;}}M = \left[ {5,10} \right]$ , which corresponds to a classic hypersonic flight envelope [Reference Gao, Chen, Sun, Huang, Wang and Chen38]. Using the algorithm proposed in section 3.2, we analysed the trend of the gap metrics under constant speed or altitude conditions. Figures 4 and 5 show the variation in the gap metrics with altitude and speed, respectively. As shown in Fig. 4, increasing the Mach number generally reduces the gap metric, indicating that the difference in model characteristics due to speed decreases with increasing speed, and this difference is more obvious at low altitudes. At different speeds, the trend of the clearance measures with altitude was similar, as shown in Fig. 5. With increasing altitude, the differences in the model characteristics due to the change in altitude increase, especially at lower Mach numbers. This indicates that the impact of altitude on the model characteristics is more pronounced at lower speeds.

Here, we consider ${\Delta}M = 0.5$ , ${\Delta}h = 1$ , and $\tau = 0.37$ ; thus, $N_{Ma} = N_{h} = 11$ , which corresponds to 121 grid points. According to gap metric theory, the flight envelope is divided into 16 subenvelopes, and the nominal points of each subenvelope are calculated, as shown in Fig. 6. The blue and red points in the figure indicate the grid and nominal points, respectively. Taking the second sub-envelope as an example, the gap metrics between each state point and other state points are calculated, and the mean value is recorded (Fig. 7). Here, point ${E_2}$ ( $h = 25{\rm{km}}$ , $M = 6.5$ ) has the smallest gap metric compared to the other state points, indicating that the model characteristics at point ${E_2}$ have the least variation from the rest of the points within the subenvelope. Therefore, the point with the minimal gap metric value ${E_2}$ can be selected as the nominal point.

Figure 6. Flight envelope division and schematic of the nominal points.

Figure 7. Gap metric mean value of the second subenvelope.

To simplify the range of values, the time-varying parameters $M$ , $h$ are normalised as follows:

(12) \begin{align} \left\{ {\begin{array}{*{20}{l}}{{\delta _M} = \left( {2M - \left( {{M_{{\rm{max}}}} + {M_{{\rm{min}}}}} \right)} \right)/\left( {{M_{{\rm{max}}}} - {M_{{\rm{min}}}}} \right)}\\[3pt] {{\delta _h} = \left( {2h - \left( {{h_{{\rm{max}}}} + {h_{{\rm{min}}}}} \right)} \right)/\left( {{h_{{\rm{max}}}} - {h_{{\rm{min}}}}} \right)}\end{array}} \right.\!, \end{align}

where ${\delta _M} \in \left[ { - 1,1} \right]$ , ${\delta _h} \in \left[ { - 1,1} \right]$ .

Several LTI systems can be obtained by Jacobi linearisation at 16 equilibrium points in the flight envelope, and the LPV model is established by fitting them with ${\delta _M}$ and ${\delta _h}$ as the scheduling parameters:

(13) \begin{align} \left\{ {\begin{array}{*{20}{l}}{\dot X = A\left( {{\delta _M},{\delta _h}} \right)X + B\left( {{\delta _M},{\delta _h}} \right)U}\\[3pt] {Y = CX}\end{array}} \right.\!, \end{align}

with $X = {\left[ {v,\gamma, h,\alpha, q} \right]^T}$ , $U = {\left[ {\sigma, {\delta _e}} \right]^T}$ ,

\begin{align*} A\left( {v,h} \right) = \left[ \begin{array}{l@{\quad}l@{\quad}l@{\quad}l@{\quad}l} {{a_{11}}} {} & { - \mu /{r^2}} {} & {{a_{13}}} & {}{{a_{14}}} & {}0\\[3pt] {{a_{21}}} & {}0 {}& {{a_{23}}} & {}{{a_{24}}} {} & 0\\[3pt] 0 & {}v & {}0 & {}0 & {}0\\[3pt] { - {a_{21}}} & {}0 & {}{{a_{43}}} & {}{{a_{44}}} & {}1\\[3pt] {{a_{51}}} & {}0 & {}{{a_{53}}} & {}{{a_{54}}} & {}{{a_{55}}}\end{array} \right], B\left( {v,h} \right) = \left[ {\begin{array}{l@{\quad}l}{{b_{11}}} & {}0\\[3pt] {{b_{21}}} & {}0\\[3pt] 0 & {}0\\[3pt] {{b_{41}}} & {}0\\[3pt] 0 & {}{{b_{52}}}\end{array}} \right], C = \left[ {\begin{array}{l@{\quad}l@{\quad}l@{\quad}l@{\quad}l}1 & {}0 & {}0 & {}0 {}& 0\\[3pt] 0 & {}1 & {}0 & {}0 & {}0\\[3pt] 0 & {}0 & {}1 & {}0 {}& 0\\[3pt] 0 & {}0 & {}0 & {}1 & {}0\\[3pt] 0 & {}0 & {}0 & {}0 & {}1\end{array}} \right]\!, \end{align*}

where $\sigma $ is the fuel equivalent ratio, ${\delta _e}$ is the elevator angle, and ${a_{11}}$ , ${a_{13}}$ , ${b_{11}}$ , etc., are functions of ${\delta _M}$ and ${\delta _h}$ (see Appendix).

4.1 Design of active disturbance rejection controller

Figure 8 shows the gain-scheduling framework, which comprises three key components: LADRC, LPV modeling based on gap metrics and gain scheduling using the GM. LADRC is central to vehicle attitude control, which is the core control method in this study. LPV models, which vary with vehicle altitude and speed, provide the foundation for calculating a stable region through the GM. The GM then dictates how control gains should be scheduled according to changes in altitude and speed, ensuring that the ADRC can meet the performance requirements of flight control.

Figure 8. Framework of the gain-scheduling control.

The primary objective of this study is to control the attitude (angle-of-attack) of the aircraft. Therefore, an LPV model is extracted for the short-period mode, with only the elevator deflection angle ${\delta _e}$ as input and the angle-of-attack $\alpha $ and pitch rate $q$ as state variables. Considering the control of the angle-of-attack, we designed an active disturbance rejection controller. According to the nonlinear model, the differential equation that describes the dynamics of $\alpha $ can be expressed as follows:

(14) \begin{align} \ddot \alpha & = \dot q - \ddot \gamma = M/{I_y} - \ddot \gamma \nonumber \\[3pt] & = 0.5\rho {v^2}sc\left( {{C_{M\alpha }} + {C_{M{\delta _e}}} + {C_{Mq}}} \right)/{I_y} - \ddot \gamma \nonumber \\[3pt] & = 0.5\rho {v^2}sc \times 0.0292{\delta _e}/{I_y} + 0.5\rho {v^2}sc\left( {{C_{M\alpha }} + {C_{Mq}}} \right)/{I_y} - \ddot \gamma \end{align}

Let ${b_0} = 0.5\rho {v^2}sc \times 0.0292/{I_y}$ , ${x_1} = \alpha $ , ${x_2} = \dot \alpha $ . According to the total disturbance concept of active disturbance rejection control, let $f = 0.5\rho {v^2}sc\left( {{C_{M\alpha }} + {C_{Mq}}} \right)/{I_y} - \ddot \gamma $ be the total disturbance, defined as the new state ${x_3}$ .

Equation (14) can then be extended to the state-space form as follows:

(15) \begin{align} \left\{ {\begin{array}{*{20}{l}}{{{\dot x}_1} = {x_2}}\\[3pt] {{{\dot x}_2} = {x_3} + {b_0}{\delta _e}}\\[3pt] {{{\dot x}_3} = \dot f} \end{array}} \right.\!\!. \end{align}

Based on the idea of online disturbance compensation, the following extended state observer is designed to estimate the total disturbance in real time:

(16) \begin{align} \left\{ {\begin{array}{*{20}{l}}{{{\dot z}_1} = {z_2} + {l_1}\left( {{x_1} - {z_1}} \right)}\\[3pt] {{{\dot z}_2} = {z_3} + {l_2}\left( {{x_1} - {z_1}} \right) + {b_0}{\delta _e}}\\[3pt] {{{\dot z}_3} = {l_3}\left( {{x_1} - {z_1}} \right)}\end{array}} \right.\!\!, \end{align}

where ${z_i}\left( {i = 1,2,3} \right)$ is the estimate of ${x_i}\left( {i = 1,2,3} \right)$ , and $L = {\left[ {{l_1},{l_2},{l_3}} \right]^\intercal}$ is the observer gain. To simplify the tuning process, the observer gains can be parameterised as $L = {\left[ {{l_1},{l_2},{l_3}} \right]^\intercal} = {\left[ {3{\omega _o},3{\omega _o}^2,{\omega _o}^3} \right]^\intercal}$ [Reference Gao17], where the observer bandwidth ${\omega _o}$ is the only tunable parameter.

The linear feedback control law is adopted as follows:

(17) \begin{align} {\delta _e} = \frac{{{k_p}\left( {{\alpha _{ref}} - {z_1}} \right) - {k_d}{z_2} - {z_3}}}{{{b_0}}}, \end{align}

Thus, the closed-loop equation can be expressed as follows:

(18) \begin{align} \left\{ {\begin{array}{*{20}{l}}{{{\dot X}_{cl}} = {A_{cl}}\left( {{\delta _M},{\delta _h}} \right){X_{cl}} + {B_{cl}}\left( {{\delta _M},{\delta _h}} \right){\alpha _{ref}}}\\[3pt] {{Y_{cl}} = {C_{cl}}{X_{cl}}} \end{array}} \right. \end{align}

with ${X_c} = {\left[ {\alpha, q,{z_1},{z_2},{z_3}} \right]^{\rm{T}}}$ ,

\begin{align*} {A_{cl}}\left( {{\delta _M},{\delta _h}} \right) = \left[ {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c@{\quad}c}{{a_{44}}} & {}1 & {}0 & {}0 & {}0\\[3pt] {{a_{54}}} & {}{{a_{55}}} & {}{ - \frac{{{b_{52}}{k_p}}}{{{b_0}}}} & {}{ - \frac{{{b_{52}}{k_d}}}{{{b_0}}}} & {}{ - \frac{{{b_{52}}}}{{{b_0}}}}\\[3pt] {3{\omega _o}} & {}0 & {}{ - 3{\omega _o}} & {}1 & {} 0\\[3pt] {3{\omega _o}^2} & {}0 & {}{ - 3{\omega _o}^2 - {k_p}} & {}{ - {k_d}} & {}0\\[3pt] {{\omega _o}^3} & {}0 & {}{ - {\omega _o}^3} & {}0 & {}0\end{array}} \right], {B_{cl}}\left( {{\delta _M},{\delta _h}} \right) = \left[ {\begin{array}{*{20}{l}}0\\[3pt] {\frac{{{b_{52}}{k_p}}}{{{b_0}}}}\\[3pt] 0\\[3pt] {\begin{array}{*{20}{l}}{{k_p}}\\[3pt] 0\\[3pt] {}\end{array}}\end{array}} \right]\!, \end{align*}

\begin{align*} {C_{cl}}\left( {{\delta _M},{\delta _h}} \right) = \left[ {\begin{array}{l@{\quad}l@{\quad}l@{\quad}l@{\quad}l}1 & {}0 & {}0 & {}0 & {}0\\[3pt] 0 & {}1 & {}0 & {}0 & {}0\\[3pt] 0 & {}0 & {}1 & {}0 & {}0\\[3pt] 0 & {}0 & {}0 & {}1 & {}0\\[3pt] 0 & {}0 & {}0 & {}0 & {}1\end{array}} \right]\!. \end{align*}

Next, we designed a suitable controller parameter gain-scheduling algorithm to tune ${k_p}$ and ${k_d}$ for different operating conditions.

4.2 Guardian map theory and stability analysis

Guardian maps are scalar maps that project a set of real matrices of order $N$ onto a region of a complex plane. Let $\nu $ map ${{\rm{R}}^{n \times n}}$ onto the entire complex plane $\mathbb{C}$ , and ${\rm{\Omega }}$ is a known region in the complex plane. The stability sets of interest are as follows [Reference Saydy, Tits and Abed27]:

(19) \begin{align} S\left( {\rm{\Omega }} \right) = \left\{ {A \in {{\rm{R}}^{n \times n}}\,:\,\lambda \left( A \right) \subset {\rm{\Omega }}} \right\}\!. \end{align}

(20) \begin{align} \nu \left( A \right) = 0{\rm{\;}}if\;and\;only\;if\;A \in \partial S\left( {\rm{\Omega }} \right)\!, \end{align}

where $\mathop S\limits^\_ \left( {\rm{\Omega }} \right)$ denotes the closure of the set $S\left( {\rm{\Omega }} \right)$ , and $\partial S\left( {\rm{\Omega }} \right)$ is the boundary of $S\left( {\rm{\Omega }} \right)$ .

Figure 9 shows the GMs of some regions.

1. (1) Hurwitz stability ( ${\rm{\Omega }} = {\mathbb{C}_ - }$ ): A GM is expressed as follows:

(21) \begin{align} \nu _H\left( A \right) = {\rm{det}}\left( {A \otimes I} \right){\rm{det}}\left( A \right)\!. \end{align}

where $ \otimes $ denotes the bialternate product [Reference Saydy, Tits and Abed27] (see Appendix).

1. (2) Stability margin: The open $\sigma $ -shifted left half-plane region ( ${\rm{R}}{{\rm{e}}_{}}\left( z \right) \lt \sigma $ ) is guarded by the following equation:

(22) \begin{align} {\nu _\sigma }\left( A \right) = {\rm{det}}\left( {A \otimes I - \sigma I \otimes I} \right){\rm{det}}\left( {A - \sigma I} \right)\!. \end{align}

1. (3) The region enclosed by the circle centred at the origin of the complex plane with radius ${\omega _n}$ . A GM that corresponds to this conic damping ratio region is expressed as follows:

(23) \begin{align} {\nu _{{\omega _n}}}\left( A \right) = {\rm{det}}\left( {A \otimes A - {\omega _n}^2I \otimes I} \right) \cdot {\rm{det}}\left( {{A^2} - {\omega _n}^2I} \right)\!. \end{align}

1. (4) A conic surface with an interior angle of 2 $\gamma $ , and a region with a damping ratio greater than $\xi = {\rm{cos}}\gamma $ . The corresponding GMs are expressed as follows:

(24) \begin{align} {\nu _\xi }\left( A \right) = {\rm{det}}\left( {{A^2} \otimes I + \left( {1 - 2{\xi ^2}} \right)A \otimes A} \right){\rm{det}}A. \end{align}

Figure 9. Guardian maps (GMs) for some stable regions.

The GMs for the stabilisation region of the flight control system considering the stable margin, natural frequency and damp ratio are as follows:

(25) \begin{align} \nu \left( A \right) = {\nu _\sigma }\left( A \right) \cdot {\nu _\xi }\left( A \right) \cdot {\nu _{{\omega _n}}}\left( A \right)\!. \end{align}

Based on these, a GM related to a specific region can be obtained. To analyse the largest variation range of the parameter that guarantees the stability of the parameterised matrices, we introduce Lemma 1 [Reference Saydy, Tits and Abed27].

(26) \begin{align} \left\{ {\begin{array}{*{20}{l}}{{r^ - } = {\rm{sup\;}}\left\{ {r \lt {r_0}:{\nu _{\rm{\Omega }}}\left[ {A\left( r \right)} \right] = 0} \right\}}\\[3pt] {{r^ + } = {\rm{inf\;}}\left\{ {r \gt {r_0}:{\nu _{\rm{\Omega }}}\left[ {A\left( r \right)} \right] = 0} \right\}}\end{array}} \right. \end{align}

${r^ - }$ is the largest real root of ${\nu _{\rm{\Omega }}}\left[ {A\left( r \right)} \right] = 0$ on the interval $\left[ {{r_{{\rm{min}}}},{r_0}} \right)$ (if it does not exist, then ${r^ - } = {r_{{\rm{min}}}}$ ), and ${r^ + }$ is the smallest real root of ${\nu _{\rm{\Omega }}}\left[ {A\left( r \right)} \right] = 0$ on the interval $\left[ {{r_0},{r_{{\rm{max}}}}} \right)$ (if it does not exist, then ${r^ + } = {r_{{\rm{max}}}}$ ).

When $r$ is scalar, given a stable region ${\rm{\Omega }}$ , the range of values of $A\left( r \right)$ with respect to the region stable $r$ can be obtained from Lemma 1 using GMs $\nu $ .

We consider the stability of two-parameter polynomial families of real matrices (or polynomials) relative to a domain ${\rm{\Omega }}$ for which $S\left( {\rm{\Omega }} \right)$ is endowed with a polynomic GM ${\nu _{\rm{\Omega }}}$

(27) \begin{align} A\left( {{r_1},{r_2}} \right) = \mathop \sum \limits_{i = 0}^{i = k} \mathop \sum \limits_{j = 0}^{j = l} r_1^ir_2^j{A_{ij}}, \end{align}

where ${r_1}$ and ${r_2}$ are the real parameters. To analyse the necessary and sufficient conditions for the stability of $A\left( {{r_1},{r_2}} \right)$ over a specific rectangle-shaped domain $\left( {{r_1},{r_2}} \right) \in \left[ {{\alpha _1},{\beta _1}} \right] \times \left[ {{\alpha _2},{\beta _2}} \right]$ Lemma 2 [Reference Saussié, Saydy, Akhrif and Berard39] is introduced. The aim is to find the largest stable open rectangle with a fixed side. First, let ${r^0} = \left( {r_1^0,r_2^0} \right)$ be a nominal parameter vector such that $A\left( {r_1^0,r_2^0} \right)$ is ${\rm{\Omega }}$ stable. By freezing one of the parameters ( ${r_2} = r_2^0$ ), the largest (open) interval containing $r_1^0$ , ${I_{{\rm{max}}}}\left( {r_2^0} \right) = \left( {r_1^ -, r_1^ + } \right)$ can be computed such that $\forall {r_1} \in {I_{{\rm{max}}}}\left( {{r_1},r_2^0} \right) \in S\left( {\rm{\Omega }} \right)$ (Lemma 1). Let ${\nu _1}\left( {{r_2}} \right) = {p_{\left[ {{r_2}} \right]}}\left( {r_{1s}^ - } \right){p_{\left[ {{r_2}} \right]}}\left( {r_{1s}^ + } \right)$ , wherein ${p_{\left[ {{r_2}} \right]}}\left( {{r_1}} \right) = {\nu _{\rm{\Omega }}}\left[ {A\left( {{r_1},{r_2}} \right)} \right]$ represents a polynomial with ${r_1}$ as the variable for each ${r_2}$ .

4.3 Controller parameter gain-scheduling algorithm based on GM

Parameter tuning based on GM is an offline tuning method that adjusts the control parameters by determining whether the closed-loop poles meet the performance requirements within the target region. The target region adopted herein is defined as Equation (25) in section 4.2, encompassing the constraints on the closed-loop poles residing in the left half of the complex plane.

Gain scheduling under all operating conditions is a 2D scheduling problem that considers both altitude and Mach number. The basic concept is as follows: The traversal is performed over the envelope of the Mach number $\left[ {5,10} \right] \times \left[ {24,34} \right]$ km. A controller is first designed at the left boundary point of the domain ( $M = 5$ , $h = 24{\rm{km}}$ , i.e. ${\delta _M} = - 1$ , ${\delta _h} = - 1$ ). Then, ${\delta _h}$ is fixed as -1, and Lemma 1 is employed to calculate the variable range for which the controller can stabilise the system. Next, new controller parameters are searched at the right boundary point of the stable range, and the corresponding stable variable range is calculated. This process is iterated until the stable variable range covers the entire Mach interval. Finally, Lemma 2 is employed to determine the maximum stable altitude region. The process is iterated at the upper bound of the altitude interval to generate a set of controller gains that ensure stability across the full operational envelope.

The procedure is as follows and is illustrated in detail in Algorithm 1.

Algorithm 1 Mach-altitude adaptive gain-scheduling algorithm based on GM

Step 1: Initialisation

1. (1) Define the target stabilisation region ${{\rm{\Omega }}_t}$ and the corresponding GM ${\nu _{{{\rm{\Omega }}_t}}}$ .

2. (2) Specify the ranges of Mach numbers $\left( {{\delta _{M,{\rm{min}}}},{\delta _{M,{\rm{max}}}}} \right)$ , and altitudes $\left( {{\delta _{h,{\rm{min}}}},{\delta _{h,{\rm{max}}}}} \right)$ .

3. (3) Let $i = 1$ , $\delta _h^i = {\delta _{h,{\rm{min}}}}$ .

Step 2: Controller parameter initialisation

Fix $\delta _h^i$ . Let $j = 1$ and ${\delta _{M0}} = {\delta _{M,{\rm{min}}}}$ . Find an initial controller parameter ${K_j} = \left[ {{k_p},{k_d}} \right]$ that ensures that the characteristic roots of $A\left( {{\delta _{M0}},{\delta _{h0}},{K_j}} \right)$ lie within ${{\rm{\Omega }}_t}$ .

Step 3: Robustness analysis

1. (1) Calculate ${\nu _{{{\rm{\Omega }}_t}}}\left[ {A\left( {{\delta _M},\delta _h^i,{K_j}} \right)} \right] = 0$ and use Lemma 1 to calculate the maximum interval $\left( {{\delta _{M,}}_j^ -, {\delta _{M,}}_j^ + } \right)$ for which ${\nu _{{{\rm{\Omega }}_t}}}\left[ {A\left( {{\delta _M},\delta _h^i,{K_j}} \right)} \right]$ is stable with respect to ${{\rm{\Omega }}_t}$ .

2. (2) If $\delta _{M,j}^ + \geqslant {\delta _{M,{\rm{max}}}}$ , go to Step 6; otherwise, go to the next step with ${\delta _{M,j}} = \delta _{M,j}^ + $ .

Step 4: Find the central controller parameter ${K_{j + 1}}$

The central value of the stable region is obtained iteratively using a 2D search algorithm to determine the central controller parameters ${K_{j + 1}}$ , which provide the maximum margin. The maximum search times and the termination accuracy are specified as ${N_s}$ and $\varepsilon $ , respectively.

1. (1) Let ${\delta _{M,j}} = \delta _{M,j}^ + $ , $m = 1$ and ${K^m} = {K_j}$ .

2. (2) While $m \lt {N_s}$ , fix ${k_d}$ and use Lemma 1 to find the largest stability interval for ${k_p}$ as $k_{}^{m + 1} = \left( {k_p^{m - } + k_p^{m + }} \right)/2$ .

3. (3) Fix ${k_p}$ and seek the largest stability interval for ${k_d}$ as $k_d^{m + 1} = \left( {k_d^{m - } + k_d^{m + }} \right)/2$ .

4. (4) When ${K^{m + 1}} - {K^m} \le \varepsilon \left( {1 + {K^m}} \right)$ , go to Step 5; otherwise, let $m = m + 1$ and go to Step 2.

5. (5) ${K_{j + 1}} = {K^{m + 1}}$ . Recalculate the GM ${\nu _{{{\rm{\Omega }}_t}}}\left[ {A\left( {{\delta _M},{\delta _h},{K_{j + 1}}} \right)} \right] = 0$ and seek the next stability interval $\left( {{\delta _{M,}}_{j + 1}^ -, {\delta _{M,}}_{j + 1}^ + } \right)$ containing ${\delta _{M,j}}$ .

6. (6) If $\delta _{M,j + 1}^ + \geqslant {\delta _{M,{\rm{max}}}}$ , proceed to Step 5; otherwise, repeat Step 4 with $j = j + 1$ .

Step 5: Determine the range of ${\delta _h}$ $\left( {\mathop \delta \limits_{ - h}^i, \mathop \delta \limits_h^{ - i} } \right)$ when ${K_j}$ stablises $A\left( {{\delta _M},\delta _h^i,{K_j}} \right)$ with respect to ${{\bf{\Omega }}_t}$

1. (1) Use Lemma 2 to find the largest stability domain $\left( {{\delta _{M,}}_j^ -, {\delta _{M,}}_j^ + } \right) \times \left( {\underline \delta _{h,j}^i,\bar \delta _{h,j}^i} \right)$ containing.

2. (2) Take their largest common intersection as $\left( {\mathop \delta \limits_{ - h}^i, \mathop \delta \limits_h^{ - i} } \right)$ .

Step 6: Terminating indicator

If $\bar \delta _h^i \gt {\delta _{h,{\rm{max}}}}$ , the end. Otherwise let $\delta _h^{i + 1} = \bar \delta _h^i$ , $i = i + 1$ , and return to Step 2.

This process generates a set of controller parameters $\left\{ {{K_1},{K_2}, \ldots {K_j}} \right\}$ and the corresponding range sets $\left\{ {\left[ {{\delta _{M,{\rm{min}}}},\delta _{M,1}^ + } \right),...,\left( {\delta _{M,j}^ -, \delta _{M,j}^ + } \right),...,\left( {\delta _{M,j}^ -, {\delta _{ma,{\rm{max}}}}} \right]} \right\} \times \left( {\underline \delta _h^i,\bar \delta _h^i} \right)$ that meet the stability requirements. An analytical expression for the controller parameters $K$ as a function of the variable ( ${\delta _M}$ , ${\delta _h}$ ) can be derived by curve fitting.

Remark: Assuming the length of each side of a rectangle is constrained to be less than the Lipschitz constant $L$ , it is guaranteed that smaller rectangles can adequately cover the entire envelope.

5.2 Robustness analysis

To account for uncertainties in the vehicle model, random variations of up to 20% were introduced to each aerodynamic parameter, and Monte Carlo simulations were performed 1,000 times to evaluate the robustness of the control system. Monte Carlo simulations were performed 1,000 times to evaluate the robustness of the control system. Figure 20 shows the angle-of-attack $\alpha $ and elevator deflection ${\delta _e}$ curves of the proposed algorithm. The results demonstrate that within the specified uncertainty range, the angle-of-attack consistently follows the commanded signal, and the controller exhibits robust and reliable performance in maintaining control effectiveness.

Figure 20. Change curve of $\alpha $ and ${\delta _e}$ with 20% aerodynamic parameters bias.

To evaluate the aircraft’s response to wind gusts, we simulated the climb phase under gust disturbance. The Horizontal Wind Model (HWM) [Reference Drob, Emmert, Meriwether, Makela, Doornbos, Conde and Klenzing40,Reference Drob, Emmert, Crowley, Picone, Shepherd, Skinner and Vincent41], which is widely recognised and adopted by the U.S. Department of Defense [Reference Madden42], was employed to mathematically represent continuous gusts. This model characterises the horizontal wind field in the stratosphere based on changes in altitude and latitude. The tangential wind velocity is expressed as follows:

(30) \begin{align} {v_{{\rm{wind}}}}\left( {h,\varphi } \right) = {v_0}\left( \varphi \right) + a\left( \varphi \right)\left( {h - {h_0}} \right) + b\left( \varphi \right){\left( {h - {h_0}} \right)^2}, \end{align}

where ${v_{{\rm{wind}}}}$ is the tangential wind speed at altitude $h$ and latitude $\varphi $ , ${h_0}$ is the reference altitude (set to 24 km herein), and $\varphi $ was set to 45°. The baseline wind speed, ${v_0}\left( \varphi \right)$ , is defined as $10 + 20 \cdot \left| {{\rm{sin}}\left( {\varphi \cdot \pi /180} \right)} \right|$ . The coefficients $a\left( \varphi \right)$ and $b\left( \varphi \right)$ , which represent the linear and quadratic altitude variations, are defined as $a\left( \varphi \right) = 2 + 0.1 \cdot \left| {{\rm{sin}}\left( {\varphi \cdot \pi /180} \right)} \right|$ and $b\left( \varphi \right) = 0.05 \cdot {\rm{cos}}\left( {\varphi \cdot \pi /180} \right)$ , respectively.

Figure 21 shows the tangential wind speed profile. The wind disturbance was applied between 100 and 200s to evaluate the system’s robustness under realistic atmospheric conditions. Figure 22 shows the angle-of-attack $\alpha $ response. The wind caused an initial sharp deviation in the angle-of-attack $\alpha $ response, but it gradually stabilised with a minor error. As shown in Table 5, the ITAE and IAE values demonstrate the robustness and accuracy of the proposed controller compared with the GM-PID controller, particularly under wind disturbance. Figure 23 compares Mach-altitude flight envelopes with and without wind disturbance. Overall, the system maintained stability and consistency under disturbances, demonstrating its robustness and reliability in overcoming wind disturbances.

Figure 21. Wind gust velocity profile.

Figure 22. Angle-of-attack ( $\alpha $ ) response under wind gust disturbance.

Table 5. Performance index comparison under wind disturbance

Figure 23. Effects of wind disturbance on flight envelopes.