3.1 Basic concept of projection

As described in Section 2.0, the flexible aircraft can be represented as a grid-based LPV model of the form given in Equation (1). If traditional LTI model order reduction methods are independently applied at each grid point, the local transformation matrices will be selected separately. This may result in inconsistencies in the preserved states or their order among different grid points [Reference Theis, Seiler and Werner23]. Here, the concept of oblique projection is introduced to avoid this issue [Reference Villemagne and Skelton38].

Definition 1. A matrix $\Pi \in {\mathbb{R}^{n \times n}}$ is a projection matrix (or idempotent matrix) of a real value vector space ${\mathbb{R}^n}$ if

(2) \begin{align} { \Pi ^2} = \Pi \end{align}

Denoting ${\rm{Im}}\left( \Pi \right) = {{\mathcal{S}}_1}$ and ${\rm{Ker}}\left( \Pi \right) = {{\mathcal{S}}_2}$ , $\Pi $ defines the projection onto ${{\mathcal{S}}_1}$ along ${{\mathcal{S}}_2}$ , where ${\rm{Im}}\left( \cdot \right)$ indicates the imagine space and ${\rm{Ker}}\left( \cdot \right)$ represents the null space. The real value vector space ${\mathbb{R}^n}$ can be decomposed as the direct sum of the imagine space and null space of the projection matrix.

(3) \begin{align} {\mathbb{R}^n} = {{\mathcal{S}}_1} \oplus {{\mathcal{S}}_2}\end{align}

When ${{\mathcal{S}}_2}$ = ${\mathcal{S}}_1^ \bot $ , the projection in Equation (2) is an orthogonal projection and can be expressed as

(4) \begin{align} {\Pi _{V,V}} = V{V^{\rm{T}}}\end{align}

where ${( \cdot )^ \bot }$ denotes orthogonal complement, $V \in {\mathbb{R}^{n \times r}}$ is an orthogonal matrix whose columns span ${{\mathcal{S}}_1}$ . For the general case where ${{\mathcal{S}}_2}$ may be distinct from ${\mathcal{S}}_1^ \bot $ , the projection is an oblique projection and can be expressed as

(5) \begin{align}{\Pi _{V,W}} = V{({W^{\rm{T}}}V)^{ - 1}}{W^{\rm{T}}}\end{align}

where $V \in {\mathbb{R}^{n \times r}}$ and $W \in {\mathbb{R}^{n \times r}}$ are two full-column rank matrices whose columns span respectively ${{\mathcal{S}}_1}$ and ${{\mathcal{S}}_2}$ .

Based on the above definition of projections, several important corollaries are given below and will be used in the projection-based LPV model order reduction.

3.2 Projection-based method

For a high order model as shown in Equation (1), based on the projection ${{\rm{\Pi }}_{V\left( \rho \right),W\left( \rho \right)}}$ , the model order reduction can be started from approximating $x$ as $V\left( \rho \right){x_r}$ .

(6) \begin{align}\frac{{d\left( {V\left( \rho \right){x_r}} \right)}}{{dt}} = A\left( \rho \right)V\left( \rho \right){x_r} + B\left( \rho \right)u + r\left( t \right)\end{align}

where $x \in {\mathbb{R}^{{n_f}}}$ and ${x_r} \in {\mathbb{R}^{{n_r}}}$ with ${n_r} \lt {n_f}$ . $V\left( \rho \right) \in {\mathbb{R}^{{n_f} \times {n_r}}}$ and the related subspace ${\mathcal{V}}$ spanned by it is called basis space. $r\left( t \right)$ is the residual introduced by the approximation. Constrain the residual to be orthogonal to a subspace ${\mathcal{W}}$ defined by a test basis $W\left( \rho \right) \in {\mathbb{R}^{{n_f} \times {n_r}}}$ as $W{(\rho )^{\rm{T}}}r\left( t \right) = 0$ . This leads to the governing equation of the Petrov-Galerkin projection model order reduction [Reference Benner, Grivet-Talocia, Quarteroni, Rozza, Schilders and Silveira39]. Considering the LTI model at a fixed grid point ${\rho _k}$ , and $V\left( {{\rho _k}} \right)$ and $W\left( {{\rho _k}} \right)$ are constant matrices with ${\rm{rank}}\left( {{W^{\rm{T}}}\left( {{\rho _k}} \right)V\left( {{\rho _k}} \right)} \right) = {n_r}$ , and there is

(7) \begin{align}{W^{\rm{T}}}\left( {{\rho _k}} \right)V\left( {{\rho _k}} \right){\dot x_r} = {W^{\rm{T}}}\left( {{\rho _k}} \right)\left( {A\left( {{\rho _k}} \right)V\left( {{\rho _k}} \right){x_r} + B\left( {{\rho _k}} \right)u} \right)\end{align}

The high-dimensional state $x$ can be projected to the reduced-order state by ${x_r} = {({W^{\rm{T}}}\left( {{\rho _k}} \right)V\left( {{\rho _k}} \right))^{ - 1}}{W^{\rm{T}}}\left( {{\rho _k}} \right)x$ as

(8) \begin{align}\begin{array}{*{20}{l}}{} {}{{{\dot x}_r} = {A_r}\left( {{\rho _k}} \right){x_r} + {B_r}\left( {{\rho _k}} \right)u}\\{} {}{{y_r} = {C_r}\left( {{\rho _k}} \right){x_r} + {D_r}\left( {{\rho _k}} \right)u}\end{array}\end{align}

where

(9) \begin{align}\begin{array}{*{20}{l}}{} {}{{A_r}\left( {{\rho _k}} \right) = {{({W^{\rm{T}}}\left( {{\rho _k}} \right)V\left( {{\rho _k}} \right))}^{ - 1}}{W^{\rm{T}}}\left( {{\rho _k}} \right)A\left( {{\rho _k}} \right)V\left( {{\rho _k}} \right) \in {\mathbb{R}^{{n_r} \times {n_r}}}}\\{} {}{{B_r}\left( {{\rho _k}} \right) = {{({W^{\rm{T}}}\left( {{\rho _k}} \right)V\left( {{\rho _k}} \right))}^{ - 1}}{W^{\rm{T}}}\left( {{\rho _k}} \right)B\left( {{\rho _k}} \right) \in {\mathbb{R}^{{n_r} \times {n_u}}}}\\{} {}{{C_r}\left( {{\rho _k}} \right) = C\left( {{\rho _k}} \right)V\left( {{\rho _k}} \right) \in {\mathbb{R}^{{n_y} \times {n_r}}}}\\{} {}{{D_r}\left( {{\rho _k}} \right) = D\left( {{\rho _k}} \right) \in {\mathbb{R}^{{n_y} \times {n_u}}}}\end{array}\end{align}

The original system can be reconstructed by $\tilde x = V\left( {{\rho _k}} \right){x_r}$ . Pre-multiplying the state equation in Equation (8) by $V\left( {{\rho _k}} \right)$ , the following equivalent system is obtained.

(10) \begin{align} \begin{array}{*{20}{l}} {\dot{\tilde{x}} = {\Pi _{V\left( {{\rho _k}} \right),W\left( {{\rho _k}} \right)}}\left( {A\left( {{\rho _k}} \right)\tilde{x} + B\left( {{\rho _k}} \right)u} \right)}\\{} {\tilde{y} = C\left( {{\rho _k}} \right)\tilde{x} + D\left( {{\rho _k}} \right)u}\end{array} \end{align}

Model order reduction techniques for LTI systems such as balanced truncation and proper orthogonal decomposition both fall under the category of Petrov-Galerkin projection-based methods. If ${\mathcal{W}} = {\mathcal{V}}$ , the related method is called a Galerkin projection method and modal truncation belongs to this category.

Remark 3. If $T = \left[ {{V_r},{V_t}} \right]$ and ${T^{ - 1}} = {[{W_r},{W_t}]^T}$ are constructed based on singular value decomposition of the product of controllability Gramian ${W_c}$ and observability Gramian ${W_o}$ of the system

(11) \begin{align}{W_{\rm{c}}}{W_{\rm{o}}}T = T{\Sigma ^2}\end{align}

where $\Sigma $ is a diagonal matrix and the entries on the diagonal are known as Hankel singular values. The transformation $x = T{x_r}$ provides a system with states hierarchically ranked according to their ability to capture the input-output characteristics [Reference Moore40]. In this case, select $V = {V_r}$ and $W = {W_r}$ , and then the related model order reduction method is balanced truncation. The details for constructing ${W_{\rm{c}}}$ and ${W_{\rm{o}}}$ are described in the literature [Reference Brunton and Kutz41].

Now considering the LPV model for the whole parameter space, the basis $V\left( \rho \right)$ and test basis $W\left( \rho \right)$ are parameter-dependent matrices. Still approximate $x$ as $V\left( \rho \right){x_r}$ and there is

(12) \begin{align}\frac{{d\left( {V\left( \rho \right){x_r}} \right)}}{{dt}} = V\left( \rho \right){\dot x_r} + \mathop \sum \limits_{i = 1}^{{n_\rho }} \frac{{\partial V\left( \rho \right)}}{{\partial {\rho _i}}}{\dot \rho _i}{x_r}\end{align}

where ${n_\rho }$ is the number of scheduling parameters. Then Equation (6) becomes

(13) \begin{align}V\left( \rho \right){\dot x_r} = \left( { - \mathop \sum \limits_{i = 1}^{{n_\rho }} \frac{{\partial V\left( \rho \right)}}{{\partial {\rho _i}}}{{\dot \rho }_i} + A\left( \rho \right)V\left( \rho \right)} \right){x_r} + B\left( \rho \right)u + r\left( t \right)\end{align}

Thus, the system matrix ${A_r}\left( {{\rho _k}} \right)$ in Equation (9) will include a term containing parameter change rates and can be expressed as

(14) \begin{align}{A_r}\left( \rho \right) = {({W^{\rm{T}}}\left( \rho \right)V\left( \rho \right))^{ - 1}}{W^{\rm{T}}}\left( \rho \right)\left( {A\left( \rho \right)V\left( \rho \right) - \mathop \sum \limits_{i = 1}^{{n_\rho }} \frac{{\partial V\left( \rho \right)}}{{\partial {\rho _i}}}{{\dot \rho }_i}} \right) \in {\mathbb{R}^{{n_r} \times {n_r}}}\end{align}

Other system matrices ${B_r}\left( \rho \right)$ , ${C_r}\left( \rho \right)$ and ${D_r}\left( \rho \right)$ have similar expressions as those given in Equation (9), just replacing ${\rho _k}$ with $\rho $ . It is noted when extending the projection-based LTI reduction methods to grid-based LPV models, two key issues need to be addressed: firstly, how to obtain a basis $V$ that allows for continuous interpolation, and secondly, how to minimise the impact of parameter variation rates on the error caused by model order reduction. These two issues will be discussed in the next subsection.

3.3 LPV model order reduction algorithm

In practice, although the basis $V\left( {{\rho _k}} \right)$ obtained from LTI model reduction at discrete grid points may not be directly applicable for continuous interpolation, a new basis $V\left( \rho \right)$ suitable for continuous interpolation can be reconstructed using the oblique projection operator. Additionally, the error introduced by parameter variation can be mitigated by incorporating constraints during the reconstruction process.

Take the modal decomposition as an example. According to Remark 2, the eigenvector matrix and its inverse can be expressed in component form as

(15) \begin{align}\begin{array}{*{20}{l}}{} {}{T\left( \rho \right) = \left[ {{v_1}\left( \rho \right),{v_2}\left( \rho \right), \cdots, {v_{{n_f}}}\left( \rho \right)} \right]}\\{} {}{{T^{ - 1}}\left( \rho \right) = {{\left[ {{w_1}\left( \rho \right),{w_2}\left( \rho \right), \cdots, {w_{{n_f}}}\left( \rho \right)} \right]}^{\rm{T}}}}\end{array}\end{align}

According to Corollary 1, the expression ${v_i}\left( \rho \right)w_i^{\rm{T}}\left( \rho \right)$ in the modal decomposition is actually a continuous oblique projection function, denoted as ${\Pi _i}\left( \rho \right)$ . Therefore,

(16) \begin{align}A\left( \rho \right) = \mathop \sum \limits_{i = 1}^{{n_f}} {\lambda _i}\left( \rho \right){v_i}\left( \rho \right)w_i^{\rm{T}}\left( \rho \right) = \mathop \sum \limits_{i = 1}^{{n_f}} {\lambda _i}\left( \rho \right){\Pi _i}\left( \rho \right)\end{align}

For any two component vectors in Equation (15) at grid point ${\rho _k}$ , considering the sum of the oblique projection matrices formed by each of them, we have

(17) \begin{align}{\Pi _i}\left( {{\rho _k}} \right) + {\Pi _j}\left( {{\rho _k}} \right) = {v_i}\left( {{\rho _k}} \right)w_i^{\rm{T}}\left( {{\rho _k}} \right) + {v_j}\left( {{\rho _k}} \right)w_j^{\rm{T}}\left( {{\rho _k}} \right) = \left[ {{v_i}\left( {{\rho _k}} \right),{v_j}\left( {{\rho _k}} \right)} \right]X \cdot {X^{ - 1}}{\left[ {{w_i}\left( {{\rho _k}} \right),{w_j}\left( {{\rho _k}} \right)} \right]^{\rm{T}}}\end{align}

where matrix $X$ represents an arbitrary invertible linear transformation of vectors ${v_i}\left( {{\rho _k}} \right)$ and ${v_j}\left( {{\rho _k}} \right)$ . According to Corollary 2, ${\Pi _i}\left( {{\rho _k}} \right) + {\Pi _j}\left( {{\rho _k}} \right)$ is a uniquely defined oblique projection matrix and its subspaces will not change with the transformation $X$ . This property can be generalised to the sum of oblique projection matrices constructed from any number of component vectors. Therefore, the local oblique projection matrix formed by $V\left( {{\rho _k}} \right)$ can be expressed as follows.

(18) \begin{align}\Pi \left( {{\rho _k}} \right) = \mathop \sum \limits_{i = 1}^{{n_f}} {\Pi _i}\left( {{\rho _k}} \right) = \mathop \sum \limits_{i = 1}^{{n_r}} {\Pi _i}\left( {{\rho _k}} \right) + \mathop \sum \limits_{i = {n_r} + 1}^{{n_f}} {\Pi _i}\left( {{\rho _k}} \right) = {\Pi _r}\left( {{\rho _k}} \right) + {\Pi _t}\left( {{\rho _k}} \right)\end{align}

As long as the partition of $\Pi \left( {{\rho _k}} \right)$ into ${\Pi _r}\left( {{\rho _k}} \right)$ and ${\Pi _t}\left( {{\rho _k}} \right)$ remains consistent at each grid point, a continuous oblique projection matrix function ${\Pi _r}\left( \rho \right)$ can be derived by interpolating ${\Pi _r}\left( {{\rho _k}} \right)$ at all grid points.

The partition of oblique projection matrices can be achieved by modal matching. This involves grouping the same modes at different grid points and further dividing them into two categories based on whether they are retained or truncated. Distinguishing the same modes (and the related projection matrices) at different grid points is realised through the metric proposed in Ref. (Reference Liu, Gao, Li and Lu25).

(19) \begin{align}{\left[ {{M_k}} \right]_{i,j}} = \frac{{\left| {{\Pi _i}\left( {{\rho _k}} \right){\Pi _j}\left( {{\rho _{k + 1}}} \right){\Pi _i}\left( {{\rho _k}} \right)} \right|}}{{\left| {{\Pi _i}\left( {{\rho _k}} \right)} \right|}}\end{align}

Denoting ${\Pi _j}\left( {{\rho _{k + 1}}} \right) = \left( {{v_j}\left( {{\rho _k}} \right) + {\rm{\Delta }}v} \right){({w_j}\left( {{\rho _k}} \right) + {\rm{\Delta }}w)^{\rm{T}}}$ , the metric can be rewritten as

(20) \begin{align}{\left[ {{M_k}} \right]_{i,j}} = \left( {\delta _{i,j}^2 - {\delta _{i,j}}{\rm{\Delta }}{w^{\rm{T}}}{\rm{\Delta }}v + w_i^{\rm{T}}\left( {{\rho _k}} \right){\rm{\Delta }}v{\rm{\Delta }}{w^{\rm{T}}}{v_i}\left( {{\rho _k}} \right)} \right)\end{align}

where ${\delta _{i,j}} = w_i^{\rm{T}}\left( {{\rho _k}} \right){v_j}\left( {{\rho _k}} \right)$ is the Kronecker delta that is defined based on the orthogonality condition of the eigenvectors. If the grid points of the LPV model are sufficiently dense, the last two terms in the above equation are both negligible higher-order infinitesimals, and the metric ${[{M_k}]_{i,j}}$ should approach the value of ${\delta _{i,j}}$ very closely. Thus, the following conditions can be used to determine whether the oblique projection matrix of two adjacent grid points should be matched to the same group.

(21) \begin{align}\left\{ {\begin{array}{*{20}{l}}{{\rm{if}}{{\left[ {{M_{}}} \right]}_{i,j}} - 1 \lt tol,{\rm{\;\;}}{\Pi _i}\left( {{\rho _k}} \right){\rm{and}}{\Pi _j}\left( {{\rho _{k + 1}}} \right){\rm{match}},}\\{{\rm{else\;if}}{{\left[ {{M_k}} \right]}_{i,j}} \lt tol,{\rm{\;\;}}{\Pi _i}\left( {{\rho _k}} \right){\rm{and}}{\Pi _j}\left( {{\rho _{k + 1}}} \right){\rm{donotmatch}},}\\{{\rm{else}},{\rm{\;\;thedistancebetween}}{\rho _k}{\rm{and}}{\rho _{k + 1}}{\rm{istoolarge}}.}\end{array}} \right.\end{align}

where $tol$ is the tolerance set by considering the higher order terms in Equation (19). After completing the matching process, the oblique projection matrices can be further classified based on the requirements for retaining or truncating specific modes, and ${\Pi _r}\left( {{\rho _k}} \right)$ can be constructed by referring to Equation (18).

Then, ${\Pi _r}\left( {{\rho _k}} \right)$ can be used to formulate a continuous transformation to achieve LPV model order reduction by multiplying with a continuous matrix function ${V_0}\left( \rho \right)$ .

(22) \begin{align}\begin{array}{*{20}{l}}{{V_r}\left( \rho \right)} {}{ = {\Pi _r}\left( \rho \right){V_0}\left( \rho \right)}\\{{W_r}\left( \rho \right)} {}{ = \Pi _r^{\rm{T}}\left( \rho \right){V_r}\left( \rho \right){{\left( {V_r^{\rm{T}}\left( \rho \right){V_r}\left( \rho \right)} \right)}^{ - 1}}}\end{array}\end{align}

Thus, the reduced-order LPV model with state consistency can be obtained by substituting ${V_r}\left( {{\rho _k}} \right)$ and ${W_r}\left( {{\rho _k}} \right)$ into $V$ and $W$ respectively in Equation (9) at each grid point ${\rho _k}$ . The selection of the matrix function ${V_0}\left( \rho \right)$ should ensure that the matrix function ${\Pi _r}\left( \rho \right){V_0}\left( \rho \right)$ has the same column rank as ${\Pi _r}\left( \rho \right)$ , so that the image space spanned by ${V_r}\left( \rho \right)$ is equal to the one spanned by ${V_0}\left( \rho \right)$ . In general, a constant basis can be used for generating ${V_0}$ .

(23) \begin{align}\begin{array}{*{20}{l}}{} {}{U{\rm{\Sigma }}{N^{\rm{T}}} = {\rm{svd}}\left( {{\Pi _r}\left( {{\rho _1}} \right),{\Pi _r}\left( {{\rho _2}} \right), \cdots, {\Pi _r}\left( {{\rho _{{n_k}}}} \right)} \right)}\\{} {}{{V_0}\left( {{\rho _k}} \right) \equiv {U_r} = U\left( {:\,,1:\,{n_r}} \right)}\end{array}\end{align}

where ${\rm{svd}}\left( \cdot \right)$ denotes the singular value decomposition. This procedure can be generalised to multi-dimensional LPV parameters spaces. For systems with slowly varying parameters, the parameter variation rate is negligible and the related term in Equation (14) may not be taken into account.

In this paper, the scheduling parameter $\rho $ is the airspeed, and thus the parameter space is one-dimensional, allowing for further derivation of a parameter-dependent matrix function ${V_0}$ . Assuming ${V_0}$ is already the desired ${V_r}$ , substituting it into Equation (22) and differentiating the equation yields

(24) \begin{align}\frac{{\partial {\Pi _r}\left( \rho \right)}}{{\partial \rho }}{V_r}\left( \rho \right) = \left( {I - {\Pi _r}\left( \rho \right)} \right)\frac{{\partial {V_r}\left( \rho \right)}}{{\partial \rho }}\end{align}

Since ${\Pi _r}\left( \rho \right)$ is already obtained, the above equation is a differential equation with respect to ${V_r}$ . The solution to this equation is not unique as ${\rm{rank}}\left( {I - {\Pi _r}\left( \rho \right)} \right) \lt {n_f}$ . Thus, the following constraint is introduced.

(25) \begin{align}{\Pi _r}\left( \rho \right)\frac{{\partial {V_r}\left( \rho \right)}}{{\partial \rho }} = 0\end{align}

Pre-multiplying $W_r^{\rm{T}}$ to both sides of the above equation yields

(26) \begin{align}\begin{array}{*{20}{l}}{W_r^{\rm{T}}\left( \rho \right){\Pi _r}\left( \rho \right)\frac{{\partial {V_r}\left( \rho \right)}}{{\partial \rho }}} {}{ = W_r^{\rm{T}}\left( \rho \right){V_r}\left( \rho \right)W_r^{\rm{T}}\left( \rho \right)\frac{{\partial {V_r}\left( \rho \right)}}{{\partial \rho }} = W_r^{\rm{T}}\left( \rho \right)\frac{{\partial {V_r}\left( \rho \right)}}{{\partial \rho }} = 0}\end{array}\end{align}

In addition to making the equation solvable, this constraint also eliminates terms related to $\dot \rho $ in Equation (14). Substituting Equation (25) into Equation (24) yields

(27) \begin{align}\frac{{\partial {\Pi _r}\left( \rho \right)}}{{\partial \rho }}{V_r}\left( \rho \right) = \frac{{\partial {V_r}\left( \rho \right)}}{{\partial \rho }}\end{align}

Using finite differences to approximate differentiation at point $\rho = {\rho _k}$ , we have

(28) \begin{align}\begin{array}{*{20}{l}}{{{\left. {\frac{{\partial {\Pi _r}\left( \rho \right)}}{{\partial \rho }}} \right|}_{{\rho _k}}} \approx \frac{{{\Pi _r}\left( {{\rho _{k + 1}}} \right) - {\Pi _r}\left( {{\rho _k}} \right)}}{{{\rho _{k + 1}} - {\rho _k}}}}\\{{{\left. {\frac{{\partial {V_r}\left( \rho \right)}}{{\partial \rho }}} \right|}_{{\rho _k}}} \approx \frac{{{V_r}\left( {{\rho _{k + 1}}} \right) - {V_r}\left( {{\rho _k}} \right)}}{{{\rho _{k + 1}} - {\rho _k}}}}\end{array}\end{align}

Substituting above equations into Equation (27) and solving it provide the construction method for the continuous transformation matrix as follows.

(29) \begin{align}{V_r}\left( {{\rho _{k + 1}}} \right) = {\Pi _r}\left( {{\rho _{k + 1}}} \right){V_r}\left( {{\rho _k}} \right)\end{align}

At point $k = 1$ , ${V_r}\left( {{\rho _1}} \right)$ can be derived based on Equation (23). The transformations ${V_r}\left( {{\rho _k}} \right)$ and ${W_r}\left( {{\rho _k}} \right)$ at all grid points can be derived by combining Equations (22) and (29), and substituting ${\Pi _r}\left( {{\rho _k}} \right)$ into this combination. Then the reduced-order model can be generated according to Equation (9). In summary, the projection-based model order reduction method can be outlined as Algorithm 1.

Algorithm 1 Oblique projection-based LPV model order reduction

4.1 Main process

Assuming that high-frequency modes, which are stable and exceed the bandwidth of the aircraft, have a negligible impact on aircraft dynamics, these modes are truncated using the modal truncation method to derive an intermediate model with 46 states from the 680th-order high-fidelity model. Then, by focusing on preserving the input-output characteristics, the model is further reduced using balanced truncation, ultimately resulting in a model with only 19 states. The final retained order is determined with reference to the $\nu $ -gap criterion, keeping it close to zero within the frequency range of [0,100] rad/s. This guarantees that the controller designed for reduced-order model can also achieve similar performance when applied to the high-fidelity model. The $\nu $ -gap between systems ${G_1}$ and ${G_2}$ is defined as follows [Reference Vinnicombe42].

(30) \begin{align}{\delta _\nu }\left( {{G_1},{G_2}} \right)\left( {j\omega } \right) = \bar \sigma \left( {{{\left( {I + {G_2}G_2^{\rm{*}}} \right)}^{ - 1/2}}\left( {{G_2} - {G_1}} \right){{\left( {I + G_1^{\rm{*}}{G_1}} \right)}^{ - 1/2}}} \right)\left( {j\omega } \right)\end{align}

where $\bar \sigma \left( \cdot \right)$ denotes the maximum singular value, ${( \cdot )^{\rm{*}}}$ represents conjugate transpose.

The $\nu $ -gaps between the high-fidelity model and the reduced-order models (46th-order and 19th-order respectively) for inputs and outputs of wingtip control surfaces and IMUs are shown in Fig. 4. In addition, the $\nu $ -gap between the 56th-order bottom-up model and the high-fidelity model is also provided for comparison. Comparing the results of the 46th-order model with those of the 56th-order bottom-up model indicates that integrating the aerodynamic-elastic model first and then applying the projection-based LPV model order reduction method allows for retaining fewer modes without increasing the $\nu $ -gap. When the model is reduced to the 19th-order, the $\nu $ -gap slightly increases, but near the critical flutter frequencies, its values are still far less than 1. This is because the dynamics of the aircraft in this frequency range are primarily dominated by the flutter modes, and the balanced truncation effectively preserves this characteristic of the system. From the comparison of pole migrations between the 19th-order and high-fidelity models in Fig. 5, it can be seen that, although the states of the model obtained through balanced truncation no longer have clear physical meaning, the pole trajectories associated with the two flutter modes remain highly consistent with those of the corresponding states in the original model.

Figure 4. ${\rm{\nu }}$ -gap between the high-fidelity model and reduced-order models.

Figure 5. Pole migrations of the high-fidelity model and the 19th-order model.

4.2 Comparison of frequency domain responses

The comparison of the frequency domain responses between the high-fidelity model and reduced-order models at 54 m/s are shown in Figs 6 and 7. In the concerned frequency range, there is no significant difference between the 46th-order model and the 56th-order bottom-up model, both of which accurately reflect the input-output response of the high-fidelity model. The plots of the 19th-order model exhibit varying differences from the high-fidelity model across different input-output pairs and frequency ranges. These differences arise because the contributions of energy associated with different input-output channels to the system vary across different frequency bands, while balanced truncation tends to retain the components that significantly contribute to the system’s energy. For instance, because the longitudinal position of the main wing’s control surfaces are very close to the centre of gravity, their contributions to the pitch attitude are minimal. Consequently, the 19th-order model does not preserve the characteristics from WL4 to the pitch rate, leading to a notable difference in the corresponding Bode plot compared to the high-fidelity model, as shown in Fig 7. However, the low-frequency responses of the rigid-body attitude, as well as the wingtip vibration responses near the flutter frequencies, are effectively captured by the 19th-order model.

Figure 6. Bode plots from the elevator and aileron to rigid-body motions.

Figure 7. Bode plots from the WL4 control surface to rigid-body and flexible motions.

4.3 Comparison of time domain responses

To further verify the fidelity of the reduced-order models, open-loop time-domain simulations are also employed. Figure 8 shows the responses of both the high-fidelity and reduced-order models to the WL4 control surface doublet command at 48 m/s. At this speed, the aircraft exhibits static stability. The disturbance from the wingtip control surface induces large-amplitude flexible vibrations and small-amplitude rigid-body attitude motions, both of which gradually decay after the disturbance ends. Apart from the inability to accurately capture the ultra-high-frequency response at the wingtip caused by sudden signal changes, the reduced-order models are generally able to reflect the aircraft’s dynamic response well. The bottom-up model shows smaller errors in pitch motion and speed variations, while the 46th-order and 19th-order models are more aligned with the high-fidelity model in terms of roll motion and wingtip vibrations.

Figure 8. Responses of different models to the doublet deflection of WL4 at 48 m/s.

Figure 9. Responses of different models to the step throttle increment at 51 m/s.

Figure 9 illustrates the responses of different models to a step throttle increment at an initial speed of 51 m/s. The increase in throttle leads to a significant rise in the aircraft’s speed, exceeding the first critical flutter speed of 52 m/s and approaching the second critical flutter speed of 55 m/s. During this process, the wingtip vibrations of the aircraft gradually increase and tend to diverge. All reduced-order models still adequately reflect the aircraft’s dynamic response. However, the fidelity of the 19th-order and 46th-order models is obvious higher than that of the bottom-up model. Overall, the model order reduction approach described in this paper achieves a lower-order reduction compared to the bottom-up model while preserving the key dynamic characteristics of the high-fidelity model. The resulting 19th-order model can effectively represent the dynamic behaviour of the high-fidelity model and will be used for subsequent controller design.

5.1 Dynamic output feedback control design

Figure 10 shows the structure of the closed-loop weighted control system used for the design of the LPV dynamic output feedback controller. The reference $r = \left[ {{r_{{V_a}}};{r_\theta };{r_\phi }} \right]$ contains target values of flight velocity, pitch angle and roll angle. The measurement output of the aircraft model is $y = \left[ {{V_a},\theta, q,\phi, p,\beta, {q_{l6}},{q_{r6}}} \right]$ . In addition to the signals that need to be tracked, it also includes pitch rate $q$ , roll rate $p$ , sideslip angle $\beta $ , as well as pitch rates ${q_{l6}}$ and ${q_{r6}}$ from the IMU L6 and R6.

Figure 10. Weighted interconnection of the closed-loop control system.

Selecting induced ${L_2}$ norm as the performance index, for the generalised plant depicted in Fig. 11, the problem of LPV dynamic output feedback controller can be defined as finding a parameter-dependent controller in the state-space representation, which can stabilise the closed-loop system and minimise induced ${L_2}$ gain from exogenous signals $w$ to controlled outputs $z$ . According to the result in literature [Reference Apkarian and Adams43], it can be formulated as a convex optimisation problem with linear matrix inequality constraints.

Table 1. Design of weighting or constraint coefficients

Figure 11. General framework of LPV control system.

In order to obtain the desired dynamic performance for reference tracking, an ideal model block ${W_{ideal}}$ is introduced in Fig. 10 to generate the ideal response for the reference $r$ .

(31) \begin{align}\begin{array}{*{20}{l}}{{W_{ideal}} = {\rm{diag}}\left( {{W_{ideal,{V_a}}},{W_{ideal,\theta }},{W_{ideal,\phi }}} \right)}\\{{W_{ideal,{V_a}}} = \frac{{0.81}}{{{s^2} + 1.44s + 0.81}}}\\{{W_{ideal,\theta }} = \frac{{2.25}}{{{s^2} + 2.4s + 2.25}}}\\{{W_{ideal,\phi }} = \frac{{2.25}}{{{s^2} + 2.4s + 2.25}}}\end{array}\end{align}

where diag( $ \cdot $ ) denotes a diagonal block composed of diagonal elements specified inside the parentheses. The ideal models corresponding to each reference signal are defined as standard second-order systems. The weighting function ${W_p}$ is used to limit the steady-state error between the ideal and actual responses and can be expressed as

(32) \begin{align}{W_p} = {\rm{diag}}\left( {{W_{p,{V_a}}},{W_{p,\theta }},{W_{p,\phi }}} \right)\end{align}

where ${W_{p,{V_a}}}$ , ${W_{p,\theta }}$ and ${W_{p,\phi }}$ are designed based on the desired error responses of speed, pitch angle and roll angle, respectively. To avoid introducing excessive dynamic states, they are set as constants, as shown in Table 1. The weighting matrix ${W_u}$ is used to limit the amplitude of controller outputs.

(33) \begin{align}{W_u} = {\rm{diag}}\left( {{W_{u,{\delta _{cs}}}},{W_{u,{\delta _{cs}}}},{W_{u,{\delta _{cs}}}},{W_{u,{\delta _t}}}} \right)\end{align}

Its diagonal components correspond sequentially to the ailerons, elevators, rudder and throttle, with specific values provided in Table 1. In addition, ${W_n}$ is used to account for the influence of measurement noise and is defined as

(34) \begin{align}{W_n} = {\rm{diag}}\left( {0.1,0.01,0.01,0.01,0.01,0.01,0.01,0.01} \right)\end{align}

where the diagonal elements represent the amplitude of the noise, corresponding to the measurement signal $y$ . By interconnecting the aircraft dynamic model and weighting functions in Fig. 10, the generalised plant $P\left( \rho \right)$ can be obtained. It can be represented by the following state-space model.

(35) \begin{align} \begin{array}{*{20}{l}}{\dot{\bar{x}} = \bar A\left( \rho \right)\bar x + {{\bar B}_u}\left( \rho \right)u + {{\bar B}_w}\left( \rho \right)w}\\{z = {{\bar C}_z}\left( \rho \right)\bar x + {{\bar D}_{zu}}\left( \rho \right)u + {{\bar D}_{zw}}\left( \rho \right)w}\\{\bar y = {{\bar C}_y}\left( \rho \right)\bar x + {{\bar D}_{yu}}\left( \rho \right)u + {{\bar D}_{yw}}\left( \rho \right)w}\end{array} \end{align}

where $u$ represents the control input, $\bar y = \left[ {y;r} \right]$ represents the generalised measurement output, $z = \left[ {{e_u};{e_p}} \right]$ represents the regulated output, and $w = \left[ {r;n} \right]$ represents the external disturbance (including reference signals). Then, based on literature [Reference Apkarian and Adams43], the linear matrix inequalities can be formulated to solve for the parameter-dependent controller.

5.2 Results and discussions

By using the methods and parameter settings described in the previous subsection, a 25th-order controller is obtained based on the reduced-order model, with corresponding performance index $\gamma $ of 1.20. Note that Algorithm 1 can also be used for order reduction of the LPV controller, which is further reduced to 21st-order. The derived controller has a fastest pole of 57.3 Hz, allowing it to be implemented on FLEXOP’s onboard hardware (with a sampling frequency of 200 Hz).

5.2.1 Frequency domain analysis

The pole trajectories of the open-loop high-fidelity model and the closed-loop system formed by connecting this model with the 21st-order controller are shown in Fig. 12. It can be seen that, within the concerned speed and frequency ranges, under the influence of the designed controller, the flutter modes that originally become unstable with increasing velocity are stabilised, and other flexible as well as rigid modes also remain stable.

Figure 12. Comparison of pole migrations of open-loop and closed-loop models.

The open-loop and closed-loop frequency domain responses of the full-order model are shown in Fig. 13. The input is selected as WL4, and the outputs are chosen to be the z-direction acceleration and pitch rate from the IMU-L6 located at the left wingtip. The response magnitude of the closed-loop system shows a slight increase only in a small range around 1.5 rad/s (the design frequency of the ideal models). While in other frequency bands, especially near the flutter frequencies, the values are obvious lower than those of the open-loop system. This indicates that the disturbance from the wingtip control surface to the wingtip vibration will be significantly reduced, as the controller effectively increases the damping of the flutter modes.

Figure 13. Comparison of Bode plots between open-loop and closed-loop models.

To evaluate the robustness of the controller, the worst-case gain and phase margins of the closed-loop system at different flight velocities are presented in Fig.14. First, consider the margins of the individual loops. Due to the inherently large damping of the rigid-body motion modes and their control by multiple control surfaces, the $p$ and $q$ loops clearly exhibit strong robustness, with minimal impact from flight velocity. The loop-at-a-time margins for the channels of outermost control surfaces and sensors are close in their value and gradually decrease as the flight velocity increases. Under the effect of the LPV controller, the rate of decline in loop-at-a-time margins with flight speed is much smoother compared to the damping decrease of the flutter modes of the aircraft. Even at a velocity close to 70 m/s, these loops still maintain a gain margin of approximately 5 dB and a phase margin of 20 degrees. Multi-loop margins address situations where multiple loops may have uncertainties or variations either individually or simultaneously, providing a more comprehensive robustness evaluation. In Fig.14(b), the inputs include the throttle and all control surface channels, while the outputs consist of all the sensor signals. Due to the mutual interactions between different loops, the margins associated with the multi-loop inputs and outputs are reduced compared to those of the individual loops. However, as a result of the integrated rigid-flexible coupling controller design, the closed-loop system is able to maintain a certain level of robustness within the designed velocity range.

Figure 14. Single-loop and multi-loop margins.

5.2.2 Time domain analysis

To verify the attitude tracking performance of the controller, closed-loop simulations are conducted using the full-order high-fidelity LPV model. To better reflect engineering practice, Gaussian white noise has been added to the sensor measurement signals. Three velocity points, 48 m/s, 54 m/s and 65 m/s, are selected to demonstrate the capability of the designed controller in scenarios where the original model is stable, a single flutter mode is unstable and two flutter modes are unstable. The step responses of pitch angle tracking and roll angle tracking are given in Figs 15 and 16. With the parameter-dependent controller, while maintaining stability, the aircraft achieves great attitude tracking in accordance with the ideal response, with very small tracking errors and performance that is almost unaffected by changes in velocity. Due to the damping effect of the controller, the symmetric and anti-symmetric responses of both wingtips during pitch and roll tracking are primarily caused by rigid-body motion. The sudden change in commands only induces small wingtip vibrations, which quickly decay to near zero.

Figure 15. Responses of pitch angle tracking under different velocities.

Figure 16. Responses of roll angle tracking under different velocities.

To validate the flutter suppression capability of the designed controller with continuous changes in scheduling parameter, velocity tracking simulations are performed within the design parameter range. Under steady-level flight at 45 m/s, a staircase command that increases by 5 m/s every 10 seconds is given. The aircraft’s velocity, pitch rate, z-direction accelerations of the IMU-L6 and IMU-R6, and the control surface responses are illustrated in Fig. 17. The aircraft velocity closely tracks the ideal response, while the wingtip vibrations are maintained within a relatively small range. At each step of the command, very small-amplitude high-frequency oscillations can be observed, which are quickly dissipated. Compared to the open-loop simulations shown in Fig. 9, the aircraft is able to smoothly pass through the critical flutter speeds under the damping effect of the controller. Due to the low-order model’s omission of the high-frequency characteristics of the high-fidelity model, although the predicted stability range of the original two flutter modes is extended to over 70 m/s, the aircraft experiences extremely high-frequency oscillations that tend to diverge when the airspeed exceeds 68 m/s. However, compared to the critical flutter speeds of 61 m/s and 63 m/s reported in the literatures [Reference Patartics, Lipták, Luspay, Seiler, Takarics and Vanek16, Reference Takarics and Vanek17], the flutter-free flight envelope of the aircraft has experienced a noticeable improvement. Although retaining more high-frequency modes during the model order reduction process can further mitigate the effects of high-frequency distortion, but it also increases the controller’s frequency, making it impractical for engineering implementation. The controller designed in this paper achieves a good balance between performance and practical applicability.

Figure 17. Responses during velocity tracking.