Interpolation-Based Modeling Methodology for Efficient Aeroelastic Control of a Folding Wing

Abstract

The aeroelastic model of a folding wing varies with different configurations, so it actually represents a parameter-varying system. Firstly, a new approach based on interpolation of local models is proposed to generate the linear parameter-varying model of a folding wing. This model is capable of predicting the aeroelastic responses during the slow morphing process and is suitable for subsequent control synthesis. The underlying inconsistencies among local linear time-invariant (LTI) models are solved through the modal matching of structural modes and the special treatment of the rational functions in aerodynamic models. Once the local LTI models are represented in a coherent state-space form, the aeroservoelastic (ASE) model at any operating point can be immediately generated by the matrix interpolation technique. Next, based on the present ASE model, the design of a parameterized controller for suppressing the gust-induced vibration is studied. The receptance method is applied to derive fixed point controllers, and the effective independence method is adopted and modified for optimal sensor placement in variable configurations, which can avoid solving ill-conditioned feedback gains. Numerical simulation demonstrates the effectiveness of the proposed interpolation-based modeling approach, and the parameterized controller exhibits a good gust mitigation effect within a wide parameter-varying range. This paper provides an effective and practical solution for modeling and control of the parameterized aeroelastic system.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

Morphing aircraft has the ability to significantly change the shape or structure during flight so that a single aircraft can adapt to various mission scenarios [1]. Over the past decades, morphing aircraft has received widespread attention in the aviation industry and aircraft manufacturers. A variety of morphing concepts have emerged, including the changing airfoil camber, wingspan and twist morphing, and variable sweep angle [2, 3]. In 2003, the Defense Advanced Research Projects Agency (DARPA), the Air Force Research Laboratory (AFRL), and the National Aeronautics and Space Administration (NASA) jointly launched the Morphing Aircraft Structures Program. Three contractors of this program, Lockheed Martin, NextGen Aeronautics, and Raytheon Missile Systems, respectively, proposed the concepts of a folding wing, flexible skin morphing wing, and telescoping wing [1]. Among them, the folding wing proposed by Lockheed Martin is an innovative morphing design that allows the structure to switch smoothly between the unfolded and folded configurations so that the optimal configuration can be adopted according to different mission requirements: the unfolded configuration is for efficient cruising, and the folded configuration is for high-speed diving [4].

Substantial changes in the structure and aerodynamic shape of the folding wings introduce specific aeroelastic behaviors that do not occur on traditional fixed wings [5]. Previous studies have provided a variety of modeling and analysis methods to investigate these aeroelastic behaviors. Lee and Weisshaar [6] used ZAERO software to generate a linear aeroelastic model of a folding wing and investigated the hinge stiffness effect on flutter dynamic pressure. Later on, Lee and Chen [7] considered freeplay nonlinearity at the wing-fold hinge in this model and performed nonlinear flutter analysis to predict the limit cycle oscillation. Wang et al. [8] presented a general aeroelastic modeling method that uses the strip theory unsteady aerodynamic model and simplified structural model. This method can perform flutter analysis on folding wings with any number of wing segments. Zhao and Hu [9] proposed a parameterized aeroelastic modeling method, which enables rapid flutter analysis of a folding wing with different configurations. Tang and Dowell [10] introduced the component modal analysis method to efficiently and accurately derive the folding wing model. In order to explain the limit cycle oscillation measured in the experiment, Attar et al. [11] extended this model to include the geometric nonlinear effect. The amplitude and dominant response frequency of the limit cycle oscillation obtained by the computational analysis were in good agreement with the experimental results. The aforementioned parametric studies on folding wings were only applicable to the case of slow morphing, and the dynamic response calculation during the rapid morphing process needs to further consider the time-varying effects on the structural dynamic characteristics. Zhao and Hu [12] combined the Craig-Bampton substructure synthesis technique with the flexible multibody dynamic approach to investigate the transient responses of a folding wing during rapid morphing processes. Hu et al. [13] proposed an integrated model by incorporating the Kriging agent model of the unsteady aerodynamic force in the time domain with the flexible multibody structural model and calculated the aeroelastic responses of a folding wing during quasi-steady morphing processes. Later on, they used this technique to study the nonlinear aeroelastic response characteristics of a folding wing with cubic stiffness [14]. The shortcoming of the Kriging model they used is that the rational function expressions of the unsteady aerodynamic force are not unique, which brings risks to the use of agent-based interpolation. Recently, Verstraete et al. [15] built a simulating system of a folding wing, which used the unsteady vortex lattice method and the finite element method (FEM) to carry out the nonlinear aeroelastic analyses in multiple flight configurations.

It can be seen that a large body of work on folding wings focuses on flutter and dynamic response predictions. The involved modeling methods can perform aeroelastic analyses under fixed or time-varying parameters and can account for miscellaneous nonlinear effects. Due to the large variation of the folding wing configuration, the time-varying morphing process, and the nonlinearity in the structure or aerodynamics, most of the time-varying nonlinear modeling techniques are quite complicated and generally far away from the control-oriented modeling. Some studies on the morphing aircraft control trend to use fixed configuration controllers as a compromise [16–18]. When the wing shape changes, the controller needs to switch online between different parameters to ensure the stability and performance of the closed-loop system. Another promising solution is to describe the nonlinear aeroelastic system as a linear parameter-varying (LPV) model that approximately captures the complex behavior during the morphing process. The LPV model simplifies the nonlinear dynamic equations of the folding wing, especially the controller can be designed in the linear system theory framework [19, 20]. Theoretically, linearizing a nonlinear model around equilibrium points in the parameter space can directly generate an LPV model [21] or a set of linear time-invariant (LTI) models for interpolation [22]. In practice, however, the original nonlinear model may be completely opaque or overly complex, making this method difficult to implement [23]. Another approach is to use the global identification or local modeling technique [24]. The former allows to generate the LPV model in a single step, but it requires the control inputs and operating conditions to be constantly changing in one experiment; such experimental conditions do not always exist in practice. The latter is based on a set of LTI models that are estimated under different fixed parameters, and the LPV model is obtained by interpolating these local LTI models. Since interpolation can be efficiently implemented in real time, the local modeling technique is currently the most convenient and effective method to establish the LPV model of the aeroelastic system [25].

When using the local modeling technique, direct interpolation of the system matrices is usually infeasible. This is because the state-space representation of a certain system is not unique, which means that the state-space matrices may be expressed in an incoherent form at different parameter points. The common solution is to convert the state-space models into a coherent form through the state coordinate transformation. Two possible canonical forms can be chosen: the companion form [23] and the modal form [26]. However, the companion transformation requires the controllability of the system inputs, and this form is known to be poorly conditioned for large-scale systems [27]. In another solution, using the modal form requires pairing the decomposed modal matrices of the local models. Most existing methods cannot handle systems with more than one parameter and require the additional assumption that the number of complex and real poles does not change over the parameter range [28, 29]. There are other solutions for coherent representation of the local state-space models, including the balanced realization [30], the least-squares approximation [31], and the SMILE technique proposed by de Caigny et al. [32]. These methods more or less have problems such as relying on experience, being difficult to implement, or having harsh conditions.

A suitable aeroservoelastic open-loop model is of great significance for studying the aeroelastic control of folding wings. In order to solve the LPV modeling problem, this paper proposes a practical local modeling technique for a typical folding wing. This approach does not require the difficult coordinate transformation of state-space LTI models (usually necessary for interpolation) but deals with the structural finite element model and the doublet lattice-based aerodynamic model in a targeted manner. Firstly, a general algorithm for structure modal matching is presented, which converts the modal matching problem into a standard linear sum assignment problem (LSAP). LSAPs are solved immediately by the Hungarian algorithm so that the local structural modes are aligned to continuously change with the folding angle. Then, the rational function approximation (RFA) results of the generalized aerodynamic force (GAF) matrices are scaled to transform the local aerodynamic models into a coherent form suitable for interpolation. The above steps eliminate all possible inconsistencies in local models and ensure that the local system matrices are continuously dependent on the scheduling parameters. In this way, whether the scheduling parameters are folding angle or flight parameters (e.g., flow speed), the state-space representation naturally has a coherent form.

As the second task of this paper, the closed-loop analysis of active aeroelastic control for the folding wing is also studied through the present aeroservoelastic model. In order to suppress the gust-induced vibration at different configurations, a parameterized controller with a multi-input multioutput (MIMO) static output feedback structure is designed by using a receptance-based method. The receptance method was originally proposed by Ram and Mottershead [33, 34]. It provides a straightforward way for vibration control of linear systems through partial pole placement. Some studies have successfully used this method to suppress the unstable aeroelastic responses due to the flutter instability [35, 36], but its control effect on external disturbances has not been evaluated. The advantage of the receptance method is that the controller can be achieved only by transfer functions from the available sensors and actuators. Therefore, the tedious tasks of model order reduction and state observer design in modern control can be avoided. Moreover, the controller in the form of static output feedback has a simple structure, so it is easy to extend to the parameterized controller through interpolation of local controllers. In this work, an additional step of optimal sensor placement is introduced to find a proper sensor layout suitable for the variable configuration of the folding wing. For this purpose, we employed the effective independence method and modified it so that it can be used in the parameter-varying system. The sensor layout is optimized by iteratively constructing the independence distribution vector and eliminating the insignificant locations. As a result, the proper sensor layout avoids solving the ill-conditioned equations of the receptance method within the parameter range. Numerical examples demonstrate that the proposed modeling and control methods are effective and reliable for the parameterized aeroelastic system in variable configurations.

2. Description of the Parameterized Aeroelastic Modeling

The schematic diagram of a folding wing geometry is shown in Figure 1(a). The folding wing structure consists of three components: the fuselage, the inboard wing, and the outboard wing. The inboard wing and the outboard wing each have a trailing-edge control surface (see Figure 1(b)). These three substructures are connected by rotating hinges driven by the servomechanism. The folding angle $\theta$ between the inboard and the fuselage can be varied from 0 deg (fully unfolded configuration) to 120 deg (fully folded configuration). During the folding process, the outboard wing remains parallel to the $x$–$y$ plane, as shown in Figure 1(c).

[figures omitted; refer to PDF]

Obviously, the equations of motion of the folding wing depend on the parameter $\theta$. Conventional nonparameterized FEM is only applicable for fixed structural configurations. When the folding angle changes, the finite element model has to be reestablished, which leads to a huge amount of repetitive modeling work. To solve this problem, Zhao and Hu [9] proposed an efficient substructure synthesis method to construct finite element models under different folding angles. It only requires modeling each substructure once. This paper adopts this method to establish the parameterized finite element model under arbitrary folding angles. To this end, the folding wing is divided into three substructures: the fuselage, the inboard wing, and the outboard wing. Finite element models of the substructures are established in their respective local coordinate systems, as shown in Figure 2. The four-node CQUAD4 elements in NASTRAN are used to discretize each substructure. For an arbitrary folding angle $\theta$, a coordinate transformation is used to express each substructural model in the global coordinate system. Substructure synthesis is performed according to the compatibility conditions for the forces and displacements at the interfaces of the substructures. In this way, the finite element model of the entire structure at an arbitrary folding angle $\theta$ is quickly obtained. The synthesized structural model of the folding wing is expressed exactly as a regular FEM equation as follows: $$\begin{array}{}\text{(1)}& {\mathbf{M}}_{\text{s}}{\ddot{\mathbf{x}}}_{\text{s}}+{\mathbf{C}}_{\text{s}}{\stackrel{.}{\mathbf{x}}}_{\text{s}}+{\mathbf{K}}_{\text{s}}{\mathbf{x}}_{\text{s}}={\mathbf{f}}_{\text{s}}+{\mathbf{f}}_{\text{g}},\end{array}$$where ${\mathbf{x}}_{\text{s}}$ is the displacement vector in global coordinates. ${\mathbf{M}}_{\text{s}}$, ${\mathbf{C}}_{\text{s}}$, and ${\mathbf{K}}_{\text{s}}$ are the mass matrix, damping matrix, and stiffness matrix, respectively, all of which depend on the parameter $\theta$. ${\mathbf{f}}_{\text{s}}$ and ${\mathbf{f}}_{\text{g}}$ are the unsteady aerodynamic forces induced by structural motion and gust disturbance, respectively.

[figures omitted; refer to PDF]

The unsteady aerodynamic model is established by using the doublet lattice method (DLM) [37]. In general, the aerodynamic boxes are independent of the finite element meshes, but in this paper, the aerodynamic boxes of the folding wing coincide with the finite element meshes. The aerodynamic influence coefficient (AIC) matrix generated by DLM is used to calculate the pressure coefficients distributed on the lifting surface under simple harmonic motion, and then the equivalent aerodynamic forces acting on the structure are derived from the force and displacement transfer relationships between structural nodes and aerodynamic grid points. Based on this, the unsteady aerodynamic forces in Equation (1) have the following form [38]: $$\begin{array}{}\text{(2)}& {\mathbf{f}}_{\text{s}}=q_{\infty }{\mathbf{G}}_{\text{ks}}^{\text{T}}{\mathbf{S}}_{\text{kj}}{\mathbf{A}}_{\text{aic}}^{-1}k,M_{\infty }{\mathbf{D}}_{\text{jk}}k{\mathbf{G}}_{\text{ks}}{\mathbf{x}}_{\text{s}},\text{(3)}& {\mathbf{f}}_{\text{g}}=q_{\infty }{\mathbf{G}}_{\text{ks}}^{\text{T}}{\mathbf{S}}_{\text{kj}}{\mathbf{A}}_{\text{aic}}^{-1}k,M_{\infty }{\mathbf{\Phi }}_{\text{jg}}k\frac{w_{\text{g}}}{U_{\infty }},\end{array}$$where ${\mathbf{A}}_{\text{aic}}$ is the AIC matrix, which is a function of the reduced frequency $k=\omega b_{\text{R}}/U_{\infty }$ and the Mach number $M_{\infty }$. $\omega$ is the angular frequency of harmonic vibration, $b_{\text{R}}$ is the reference half chord length, and $U_{\infty }$ is the flow speed. $q_{\infty }=0.5{\rho }_{\infty }{U}_{\infty }^2$ is the dynamic pressure, and ${\rho }_{\infty }$ is the air density. ${\mathbf{S}}_{\text{kj}}$ is the integration matrix used to concentrate the distributed pressures to the aerodynamic grid points. ${\mathbf{D}}_{\text{jk}}$ is the substantial differentiation matrix derived from the airflow boundary conditions. ${\mathbf{G}}_{\text{ks}}$ is the spline matrix obtained by the infinite plate spline method [39], which is used for mesh deformation interpolation and equivalent aerodynamic force calculation. ${\mathbf{\Phi }}_{\text{jg}}$ is the gust mode vector, including the time delay and amplitude ratio of the harmonic gust excitation on each aerodynamic box [40]. $w_{\text{g}}$ is the gust velocity disturbance at the gust reference point, and the gust reference point is selected at the leading-edge of the fuselage, i.e., $x=0$. It should be noted that in Equations (2) and (3), except that ${\mathbf{D}}_{\text{jk}}$ and ${\mathbf{S}}_{\text{kj}}$ are not related to the folding angle $\theta$, the rest of the matrices depend on $\theta$.

The structural motion can be expressed as the superposition of each modal motion. Taking the control surface deflection modes into consideration, the following modal coordinate transformation is introduced: $$\begin{array}{}\text{(4)}& {\mathbf{x}}_{\text{s}}={\mathbf{\Phi }}_{\text{h}}{\mathbf{\Phi }}_{\text{c}}\begin{array}{c}\mathbf{\xi }\\ \mathbf{\delta }\end{array},\end{array}$$where ${\mathbf{\Phi }}_{\text{h}}$ is the structure modal matrix and $\mathbf{\xi }$ is the vector of corresponding modal displacements. ${\mathbf{\Phi }}_{\text{c}}$ is the matrix of control surface modes, which contains rigid displacements due to the unit deflection of each control surface. $\mathbf{\delta }$ is the vector of actual deflections of the control surfaces.

Combining Equations (1)–(4), the aeroelastic equation of the folding wing in modal coordinates is obtained as follows: $$\begin{array}{}\text{(5)}& {\mathbf{M}}_{\text{hh}}\ddot{\mathbf{\xi }}+{\mathbf{C}}_{\text{hh}}\stackrel{.}{\mathbf{\xi }}+{\mathbf{K}}_{\text{hh}}\mathbf{\xi }+{\mathbf{M}}_{\text{hc}}\ddot{\mathbf{\delta }}=q_{\infty }{\mathbf{Q}}_{\text{hh}}\mathbf{\xi }+q_{\infty }{\mathbf{Q}}_{\text{hc}}\mathbf{\delta }+q_{\infty }{\mathbf{Q}}_{\text{hg}}\frac{w_{\text{g}}}{U_{\infty }},\end{array}$$where ${\mathbf{M}}_{\text{hh}}$, ${\mathbf{C}}_{\text{hh}}$, and ${\mathbf{K}}_{\text{hh}}$ are the modal mass, modal damping, and modal stiffness matrices in the diagonal form, respectively. ${\mathbf{M}}_{\text{hc}}$ is the mass coupling matrix, which represents the inertial effect induced by the control surfaces. ${\mathbf{Q}}_{\text{hh}}$, ${\mathbf{Q}}_{\text{hc}}$, and ${\mathbf{Q}}_{\text{hg}}$ are the frequency-domain GAF matrices due to structural modes, control surface modes, and gust modes, respectively. Matrices in Equation (5) are expressed as follows: $$\begin{array}{c}\text{(6)}& {\mathbf{M}}_{\text{hh}}={\mathbf{\Phi }}_{\text{h}}^{\text{T}}{\mathbf{M}}_{\text{s}}{\mathbf{\Phi }}_{\text{h}},{\mathbf{C}}_{\text{hh}}={\mathbf{\Phi }}_{\text{h}}^{\text{T}}{\mathbf{C}}_{\text{s}}{\mathbf{\Phi }}_{\text{h}},\\ {\mathbf{K}}_{\text{hh}}={\mathbf{\Phi }}_{\text{h}}^{\text{T}}{\mathbf{K}}_{\text{s}}{\mathbf{\Phi }}_{\text{h}},\\ {\mathbf{M}}_{\text{hc}}={\mathbf{\Phi }}_{\text{h}}^{\text{T}}{\mathbf{M}}_{\text{s}}{\mathbf{\Phi }}_{\text{c}},\\ {\mathbf{Q}}_{\text{hh}}={\mathbf{\Phi }}_{\text{h}}^{\text{T}}{\mathbf{G}}_{\text{ks}}^{\text{T}}{\mathbf{S}}_{\text{kj}}{\mathbf{A}}_{\text{aic}}^{-1}k,M_{\infty }{\mathbf{D}}_{\text{jk}}k{\mathbf{G}}_{\text{ks}}{\mathbf{\Phi }}_{\text{h}},\\ {\mathbf{Q}}_{\text{hc}}={\mathbf{\Phi }}_{\text{h}}^{\text{T}}{\mathbf{G}}_{\text{ks}}^{\text{T}}{\mathbf{S}}_{\text{kj}}{\mathbf{A}}_{\text{aic}}^{-1}k,M_{\infty }{\mathbf{D}}_{\text{jk}}k{\mathbf{G}}_{\text{ks}}{\mathbf{\Phi }}_{\text{c}},\\ {\mathbf{Q}}_{\text{hg}}={\mathbf{\Phi }}_{\text{h}}^{\text{T}}{\mathbf{G}}_{\text{ks}}^{\text{T}}{\mathbf{S}}_{\text{kj}}{\mathbf{A}}_{\text{aic}}^{-1}k,M_{\infty }{\mathbf{\Phi }}_{\text{jg}}k.\end{array}$$

So far, the parameterized aeroelastic model of the folding wing in modal coordinates has been initially established. All the structural and aerodynamic matrices in Equation (5) depend on parameter $\theta$. For any fixed $\theta$, system analysis and response calculation are easy to carry out. Note that the aerodynamic forces in the current model are expressed in the frequency domain. If time-domain aeroelastic analyses are desired, the GAF matrices can be transformed into the time domain through the RFA technique [41]. However, the current model is not suitable for simulation under the time-varying folding angle. It is extremely inefficient to perform online substructure synthesis, structure modal decomposition, and aerodynamic matrix calculation at each time step. Therefore, a modeling method for quickly obtaining an aeroelastic model at any point in the parameter space is desired.

3. Interpolation-Based Modeling Methodology

For the folding wing system, both the folding angle and the aerodynamic parameters (e.g., flow speed) will undergo variations. The response analysis and control synthesis are expected to be carried out under the LPV framework. Using a set of local models to compute an interpolating LPV model is a practical and efficient modeling strategy [25]. This method begins with the discretization of a given parameter space, which generates a set of grid points called operating points; then, the LTI models (local models) are prebuilt for fixed operating conditions at each operating point; finally, the LPV model is computed by interpolating these local LTI models. Although the local LTI models are obtained for fixed parameters and do not incorporate the time-varying effects, in case the parameter variations are slow relative to the dynamic characteristics of the system, the dynamic parameter-dependent part in the system can be ignored [32, 42]. Therefore, in this paper, while applying the local modeling technique, we assume that the folding process and flight conditions undergo slow and smooth variations.

Model interpolation requires that all local models have a coherent form; that is, the system matrix should change continuously over the considered parameter range. Note that Equation (1) is expressed in the physical coordinate system. The matrices in Equation (1) are continuously dependent on the folding angle $\theta$. However, this is not the case in the modal coordinate system of Equation (5). It can be inferred that only the structure modal matrix ${\mathbf{\Phi }}_{\text{h}}$ may cause inconsistency. Since the mode crossing is a common phenomenon in the parameter-varying system, the structural modes with the same order cannot be directly interpolated because they often do not change continuously (do not belong to the same mode branch). Besides, the mode shapes may exhibit abrupt sign changes, which will also make interpolation impossible. Another inconsistency occurs in the subsequent RFA procedure for transforming the frequency-domain aerodynamic forces into the time domain. The commonly used minimum-state (MS) method generates nonunique coefficient matrices that cannot be directly interpolated. This section presents a structure modal matching algorithm and a coherent RFA representation method, which eliminates the underlying inconsistencies in local models and paves the way for model interpolation.

3.1. Modal Matching and Alignment

To perform the model interpolation, a set of local aeroelastic models of the folding wing are established under the $n_{\text{t}}$ folding angles ${\theta }_1,\cdots ,{\theta }_{n_{\text{t}}}$, and the following structural eigenvalue problems are solved: $$\begin{array}{}\text{(7)}& {\mathbf{K}}_{\text{s},l}{\mathbf{\Phi }}_{\text{h},l}={\mathbf{M}}_{\text{s},l}{\mathbf{\Phi }}_{\text{h},l}{\mathbf{\Omega }}_{\text{h},l}^2,\hspace{1em}l=1,\cdots ,n_{\text{t}},\end{array}$$where $$\begin{array}{c}\text{(8)}& {\mathbf{\Omega }}_{\text{h},l}=\mathrm{diag}\begin{array}{ccc}{\omega }_{1,l}& \cdots & {\omega }_{n_{\text{h}},l}\end{array},{\mathbf{\Phi }}_{\text{h},l}=\begin{array}{ccc}{\mathbf{\phi }}_{1,l}& \cdots & {\mathbf{\phi }}_{n_{\text{h}},l}\end{array}.\end{array}$$

${\mathbf{\Omega }}_{\text{h},l}$ is a diagonal matrix whose elements are angular natural frequencies of the structure. The column vectors in matrix ${\mathbf{\Phi }}_{\text{h},l}$ are structural mode shapes normalized for unit-generalized masses. $n_{\text{h}}$ is the number of the retained structural modes.

Typical eigenvalue solvers generally sort structural modes by natural frequencies in ascending order. When the mode crossing occurs, modes in matrix ${\mathbf{\Phi }}_{\text{h},l}$ should be rearranged before interpolation to ensure continuous change between parameter points. Considering that the structural modes in the same branch have similar vibration patterns, we can compare the natural frequencies and mode shapes under the adjacent parameters and rearrange the orders to match these dynamic properties. For the two adjacent folding angles ${\theta }_{l-1}$ and ${\theta }_l$, define the following distance metric to measure the dynamic similarity between the two modes: $$\begin{array}{}\text{(9)}& d_{i,j}={\omega }_{i,l-1}-{\omega }_{j,l}\cdot 1-\sqrt{{\text{MAC}}_{i,j}},\hspace{1em}i,j=1,\cdots ,n_{\text{h}},\end{array}$$where the modal assurance criterion (MAC) is calculated as $$\begin{array}{}\text{(10)}& {\text{MAC}}_{i,j}=\frac{{{\mathbf{\phi }}_{i,l-1}^{\text{T}}{\mathbf{\phi }}_{j,l}}^2}{{\mathbf{\phi }}_{i,l-1}^{\text{T}}{\mathbf{\phi }}_{i,l-1}\cdot {\mathbf{\phi }}_{j,l}^{\text{T}}{\mathbf{\phi }}_{j,l}}.\end{array}$$

The distance metric uses the linear distance of two natural frequencies and is weighted by the MAC. MAC takes value in the interval $0,1$; a larger MAC value indicates a higher degree of linear correlation between the associated mode shapes. Assuming that the modal data at point $l-1$ has been sorted in a proper order, the goal of modal matching is to compare the modal data at point $l$ with point $l-1$ and give a correct pairing by minimizing the total distance metric. From the perspective of the graph theory, the structural modes to be investigated can be naturally divided into the following two sets of vertices: $$\begin{array}{c}\text{(11)}& V_{l-1}=\begin{array}{ccc}{\omega }_{1,l-1}& \cdots & {\omega }_{n_{\text{h}},l-1}\end{array},V_l=\begin{array}{ccc}{\omega }_{1,l}& \cdots & {\omega }_{n_{\text{h}},l}\end{array}.\end{array}$$

Define the edge set $E_l=V_{l-1}\times V_l$; that is, any two vertices in $V_{l-1}$ and $V_l$ are connected by a unique edge, thus forming a complete bipartite graph $G_l=V_{l-1},V_l;E_l$. The distance metric $d_{i,j}$ given in Equation (9) specifies the cost of each edge. The modal matching problem is reformulated as finding a perfect matching in the graph $G_l$, which minimizes the sum of the costs. This problem is a typical LSAP and can be solved by the Hungarian algorithm, which is a mature solution in the field of multiobject tracking [43]. The Hungarian algorithm was originally presented by Kuhn [44]. An improved version based on the shortest augmenting path technique [45] is currently more popular because it achieves an $O{n}_{\text{h}}^3$ worst-case time complexity. Therefore, the modal matching problem can be efficiently and completely solved.

Figure 3 briefly illustrates the modal matching and alignment procedure. According to the modal matching results, the precalculated natural frequencies and modal vectors are aligned (reordered) successively for $l=2,\cdots ,n_{\text{t}}$. In order to ensure numerical continuity, each aligned mode shape vector should be adjusted to $c_{i,l}\cdot {\mathbf{\phi }}_{i,l}$, where the coefficient $c_{i,l}$ is calculated by a signum function as $$\begin{array}{}\text{(12)}& c_{i,l}=\mathrm{sign}{\mathbf{\phi }}_{i,l-1}^{\text{T}}{\mathbf{\phi }}_{\text{i},l},\hspace{1em}i=1,\cdots ,n_{\text{h}}.\end{array}$$

[figure omitted; refer to PDF]

After modal matching and alignment, the set of aeroelastic models established by Equation (5) have a coherent form under all folding angles. Next, a method to construct coherent time-domain models is given.

3.2. Coherent RFA Representation

The GAF matrices in Equation (5) are expressed in the frequency domain. To obtain the time-domain aeroelastic model, the frequency-domain aerodynamic forces need to be converted to the Laplace domain through the RFA technique. The RFA technique uses the tabulated GAF data at several discrete reduced frequencies to fit a specific rational function of the Laplace variable $s$. The most commonly used MS method gives a general expression of RFA as [41] $$\begin{array}{}\text{(13)}& \mathbf{Q}p={\mathbf{A}}_0+{\mathbf{A}}_{1}p+{\mathbf{A}}_2{p}^2+\mathbf{D}{p\mathbf{I}-\mathbf{R}}^{-1}\mathbf{E}p,\end{array}$$where $p={sb}_{\text{R}}/U_{\infty }$ is the nondimensional Laplace variable and $\mathbf{Q}$ is the combination of the following three GAF terms: $$\begin{array}{}\text{(14)}& \mathbf{Q}p=\begin{array}{ccc}{\mathbf{Q}}_{\text{hh}}p& {\mathbf{Q}}_{\text{hc}}p& {\mathbf{Q}}_{\text{hg}}p\end{array}.\end{array}$$

Matrix $\mathbf{R}$ in Equation (13) is a diagonal matrix containing manually specified negative numbers called the aerodynamic lag roots, and the other coefficient matrices are solved by the nonlinear least-squares iterative procedure. Although the fitting algorithm (irrespective of the way it is implemented) gives an explicit RFA expression, the $\mathbf{D}$ and $\mathbf{E}$ matrices in the lag term are not unique. Specifically, for any nonsingular matrix $\mathbf{T}$, the following equality always holds: $$\begin{array}{}\text{(15)}& \mathbf{D}{p\mathbf{I}-\mathbf{R}}^{-1}\mathbf{E}=\mathbf{D}\mathbf{T}{p\mathbf{I}-{\mathbf{T}}^{-1}\mathbf{R}\mathbf{T}}^{-1}{\mathbf{T}}^{-1}\mathbf{E}.\end{array}$$

When building the parameterized system model, RFAs are performed independently at each parameter point. The nonuniqueness of the RFA expression brings inconsistency between local aerodynamic models, which will lead to incorrect interpolation results. In fact, the right side of Equation (15) contains all possible expressions. It can be further proved that ${\mathbf{T}}^{-1}\mathbf{R}\mathbf{T}=\mathbf{R}$ if and only if $\mathbf{T}$ is a diagonal matrix. Therefore, considering the prescribed matrix $\mathbf{R}$, the uncertainty of the RFA expression is limited to the case where $\mathbf{T}$ is a diagonal matrix, and the nonuniqueness is manifested in the arbitrary scaling of the column vectors of $\mathbf{D}$ and the row vectors of $\mathbf{E}$. Assuming that the matrix $\mathbf{R}$ is consistent or changes continuously with the parameters, then a unique scaling matrix $\mathbf{T}$ can be constructed according to the RFA results under adjacent parameters so that the unique and coherent RFA expression can be obtained.

Consider the case where the folding angle is the variable parameter. For the two adjacent folding angles ${\theta }_{l-1}$ and ${\theta }_l$, take the matrix ${\mathbf{D}}_{l-1}$ at point $l-1$ as the standard reference, and the diagonal matrix ${\mathbf{T}}_l$ at point $l$ is constructed as follows: $$\begin{array}{}\text{(16)}& {\mathbf{T}}_l=\mathrm{diag}\begin{array}{ccc}{t}_{1,l}& \cdots & t_{n_{\text{a}},l}\end{array},\end{array}$$where $n_{\text{a}}$ is the number of aerodynamic lag roots and the diagonal element $t_{i,l}$ is calculated as follows: $$\begin{array}{}\text{(17)}& t_{i,l}=\mathrm{sign}{\mathbf{d}}_{i,l-1}^{\text{T}}{\mathbf{d}}_{i,l}\cdot \frac{{{\mathbf{d}}_{i,l-1}}_2}{{{\mathbf{d}}_{i,l}}_2},\hspace{1em}i=1,\cdots ,n_{\text{a}},\end{array}$$where ${\mathbf{d}}_{i,l-1}$ and ${\mathbf{d}}_{i,l}$ are the $i$-th column of matrices ${\mathbf{D}}_{l-1}$ and ${\mathbf{D}}_l$, respectively. After obtaining ${\mathbf{T}}_l$, adjust the matrices ${\mathbf{D}}_l$ and ${\mathbf{E}}_l$ according to Equation (15). The corresponding column vectors of the scaled matrix ${\mathbf{D}}_l$ and the reference matrix ${\mathbf{D}}_{l-1}$ thus have the same 2-norm and consistent direction.

When the folding angle and Mach number are variable parameters, the RFA data should be generated on the two-dimensional parameter grid points. In this case, the processing method for coherent RFA representation should be as follows: first, fix the first Mach number and successively adjust the coefficient matrices for the sequential folding angles; then, fix each folding angle and successively adjust the coefficient matrices for the sequential Mach numbers in the same way.

After obtaining a set of coherent local RFAs, the state-space representation of the LPV system can be generated by combining Equation (13) with the aeroelastic equation given in Equation (5), as shown below: $$\begin{array}{}\text{(18)}& \stackrel{.}{\mathbf{x}}=\mathbf{A}\mathbf{\rho }\mathbf{x}+\mathbf{B}\mathbf{\rho }\mathbf{u}+{\mathbf{B}}_{\text{w}}\mathbf{\rho }{\tilde{\mathbf{w}}}_{\text{g}},\text{(19)}& \mathbf{y}=\mathbf{C}\mathbf{\rho }\mathbf{x}+{\mathbf{D}}_{\text{w}}\mathbf{\rho }{\tilde{\mathbf{w}}}_{\text{g}},\end{array}$$where $\mathbf{\rho }$ is the vector of scheduling parameters including the folding angle $\theta$. The state vector $\mathbf{x}$ contains modal displacement vector $\mathbf{\xi }$ and its time derivative $\stackrel{.}{\mathbf{\xi }}$, as well as the augmenting aerodynamic states introduced by the RFA. The motion of the control surfaces is generally determined by the actuator model, and the above state-space equation also integrates the actuator equation of the state-space form [46] so that the control input vector $\mathbf{u}$ represents the control surface deflection commands, while the actual deflections in $\mathbf{\delta }$ are incorporated into $\mathbf{x}$ as the augmenting states. The gust input vector ${\tilde{\mathbf{w}}}_{\text{g}}$ contains gust velocity $w_{\text{g}}$ and its time derivative ${\stackrel{.}{w}}_{\text{g}}$. There are various standard methods for generating gust disturbances [47]. Two types of gust models, the 1-cos discrete gust and the Dryden continuous gust, are used in the following time-domain simulations. The system output vector $\mathbf{y}$ may include responses at the FEM nodes (monitoring points) and/or internal loads (bending moment, etc.). For more details about constructing the state-space matrices, see Ref. [48].

Taking the fixed altitude flight as an example, the scheduling parameters in $\mathbf{\rho }$ thus consist only of the folding angle $\theta$ and the flow speed $U_{\infty }$. Depending on the local modeling approach, the LPV model shown in Equations (18) and (19) is established by local LTI models on the two-dimensional parameter grid. The methods to deal with structural modes and RFA coefficients proposed in this section successfully solved the possible inconsistencies in the local models. Whether it is the change of the folding angle or the flow speed, the local models are guaranteed to be numerically continuous. In this way, by interpolating the state-space matrices, the model at any point in the parameter space can be quickly obtained, and the LPV simulation and control synthesis can also be performed efficiently.

4. Receptance-Based Active Aeroelastic Control

In this section, the receptance method [33, 34] is applied to design a parameterized controller for the folding wing, which is expected to reduce the structural vibration and additional loads induced by gust disturbances. This method uses the receptance transfer function extracted from the analysis model or the identified model to design the controller and achieves the active vibration control through partial pole placement. The receptance-based controller is theoretically solved under fixed parameters, and the control strategy under variable parameters can be realized by interpolation of the local controllers. Compared with the standard LPV control design such as gain-scheduled control [49, 50], interpolation of the fixed point controllers is easy to implement, but the system performance cannot be guaranteed when parameters change rapidly. However, under the assumption of slow parameter variation, the interpolation approach can avoid introducing the conservativeness involved in the LPV approaches and thus obtain a better control effect than the gain-scheduled control [51].

In order to avoid blindly selecting the sensor locations, an additional step of optimal sensor placement is introduced in the control design strategy. As will be seen below, this step also avoids solving the ill-conditioned equations of the receptance method.

4.1. Receptance-Based Control

The receptance-based controller is more suitable to be derived using the frequency-domain equation. Converting Equation (5) into the Laplace domain, we obtain $$\begin{array}{}\text{(20)}& s^2{\mathbf{M}}_{\text{hh}}+s{\mathbf{C}}_{\text{hh}}+{\mathbf{K}}_{\text{hh}}-q_{\infty }{\mathbf{Q}}_{\text{hh}}s\mathbf{\xi }s=-s^2{\mathbf{M}}_{\text{hc}}+q_{\infty }{\mathbf{Q}}_{\text{hc}}s\mathbf{\delta }s+\mathbf{d}s,\end{array}$$where $\mathbf{d}s\in {ℂ}^{n_{\text{h}}}$ denotes the gust disturbances. Generally, the dynamic model of the actuators driving the control surfaces is represented by a specific transfer function as $$\begin{array}{}\text{(21)}& \mathbf{\delta }s={\mathbf{H}}_{\text{c}}s\mathbf{u}s,\end{array}$$where $\mathbf{u}s\in {ℂ}^{n_{\text{c}}}$ denotes the actuator commands and $\mathbf{\delta }s\in {ℂ}^{n_{\text{c}}}$ stands for the actual deflections of the control surfaces. The output equation of the system is $$\begin{array}{}\text{(22)}& \mathbf{y}s={\mathbf{\Phi }}_{\text{s}}\mathbf{\xi }s,\end{array}$$where $\mathbf{y}s\in {ℂ}^{n_{\text{s}}}$ is the output vector and ${\mathbf{\Phi }}_{\text{s}}\in {ℝ}^{n_{\text{s}}\times n_{\text{h}}}$ is the sensor modal matrix. Let ${\mathbf{G}}_{\text{d}}$ and ${\mathbf{G}}_{\text{v}}\in {ℝ}^{n_{\text{c}}\times n_{\text{s}}}$ be the proportional displacement and velocity output feedback gains to be solved, respectively. The feedback control law can be written as $$\begin{array}{}\text{(23)}& \mathbf{u}s=s{\mathbf{G}}_{\text{v}}+{\mathbf{G}}_{\text{d}}\mathbf{y}s.\end{array}$$

Set the control inputs and external disturbances to be zero and solve Equation (20) to obtain the eigenvalues ${\lambda }_i$ and eigenvectors ${\mathbf{w}}_i$ of the open-loop system, expressed as $$\begin{array}{}\text{(24)}& {\lambda }_i^2{\mathbf{M}}_{\text{hh}}+{\lambda }_i{\mathbf{C}}_{\text{hh}}+{\mathbf{K}}_{\text{hh}}-q_{\infty }{\mathbf{Q}}_{\text{hh}}{\lambda }_i{\mathbf{w}}_i=0,\hspace{1em}i=1,\cdots ,2n_{\text{h}}.\end{array}$$

Due to the nonlinearity of the ${\mathbf{Q}}_{\text{hh}}$ term, Equation (24) does not express a typical quadratic eigenvalue problem. The number of eigenvalues is generally not equal to twice the number of the structural modes. Consider that Equation (24) is also the equation used in linear flutter analysis. The classical $p$-$k$ method only finds $n_{\text{h}}$ solutions (ignoring the complex conjugate) of the equation under each dynamic pressure, which is enough to reveal the governing modes of the aeroelastic system [52]. This paper uses the RFA method instead of the flutter analysis method to calculate the open-loop eigenvalues (more accurate unsteady aerodynamic forces), but only the $2n_{\text{h}}$ eigenvalues closely related to the structural motion are retained.

The target for partial pole placement is to assign a part of closed-loop eigenvalues to the desired values ${\mu }_i$, $i=1,\cdots m$, through the feedback gains, where the assigned eigenvalues are distinct from all open-loop eigenvalues; meanwhile, the rest of the closed-loop eigenvalues are unchanged, i.e., ${\mu }_i={\lambda }_i$, $i=m+1,\cdots ,2n_{\text{h}}$. Substituting Equations (21)–(23) into Equation (20), the assigned eigenvalues ${\mu }_i$ and eigenvectors ${\mathbf{v}}_i$ of the closed-loop system should satisfy the following equation: $$\begin{array}{}\text{(25)}& {\mathbf{v}}_i=\mathbf{H}{\mu }_i\mathbf{P}{\mu }_i{\mu }_i{\mathbf{G}}_{\text{v}}+{\mathbf{G}}_{\text{d}}{\mathbf{\Phi }}_{\text{s}}{\mathbf{v}}_i,\hspace{1em}i=1,\cdots ,m,\end{array}$$where $\mathbf{H}{\mu }_i$ and $\mathbf{P}{\mu }_i$ are the values of the receptance transfer function $\mathbf{H}s$ and the control channel transfer function $\mathbf{P}s$, respectively, defined as $$\begin{array}{c}\text{(26)}& \mathbf{H}s={s^2{\mathbf{M}}_{\text{hh}}+s{\mathbf{C}}_{\text{hh}}+{\mathbf{K}}_{\text{hh}}-q_{\infty }{\mathbf{Q}}_{\text{hh}}s}^{-1},\mathbf{P}s=-s^2{\mathbf{M}}_{\text{hc}}+q_{\infty }{\mathbf{Q}}_{\text{hc}}s{\mathbf{H}}_{\text{c}}s.\end{array}$$

The following in this subsection describes how to solve the feedback gains. Introduce the weighting parameter vector ${\mathbf{\alpha }}_i\in {ℂ}^{n_{\text{c}}}$ as follows: $$\begin{array}{}\text{(27)}& {\mathbf{\alpha }}_i={\mu }_i{\mathbf{G}}_{\text{v}}+{\mathbf{G}}_{\text{d}}{\mathbf{\Phi }}_{\text{s}}{\mathbf{v}}_i,\hspace{1em}i=1,\cdots ,m.\end{array}$$

The value in vector ${\mathbf{\alpha }}_i$ indicates the degree of participation of each control input for controlling the $i$-th mode. Equation (25) can therefore be written as $$\begin{array}{}\text{(28)}& {\mathbf{v}}_i=\mathbf{H}{\mu }_i\mathbf{P}{\mu }_i{\mathbf{\alpha }}_i,\hspace{1em}i=1,\cdots ,m.\end{array}$$

When ${\mathbf{\alpha }}_i$ is given in advance, the closed-loop eigenvector ${\mathbf{v}}_i$ is determined by Equation (28), and the pole placement problem is transformed into finding feedback gains satisfying Equation (27). In addition, the unchanged closed-loop eigenvalues also satisfy the following condition: $$\begin{array}{}\text{(29)}& \mathbf{P}{\lambda }_i{\lambda }_i{\mathbf{G}}_{\text{v}}+{\mathbf{G}}_{\text{d}}{\mathbf{\Phi }}_{\text{s}}{\mathbf{w}}_i=0,\hspace{1em}i=m+1,\cdots ,2n_{\text{h}}.\end{array}$$

Since $\mathbf{P}s$ can take various forms, Equation (29) is modified to the following sufficient condition: $$\begin{array}{}\text{(30)}& {\lambda }_i{\mathbf{G}}_{\text{v}}+{\mathbf{G}}_{\text{d}}{\mathbf{\Phi }}_{\text{s}}{\mathbf{w}}_i=0,\hspace{1em}i=m+1,\cdots ,2n_{\text{h}}.\end{array}$$

Combining Equations (27) and (30), the following linear equation is deduced to solve the feedback gain matrix: $$\begin{array}{}\text{(31)}& \mathbf{G}\mathbf{Y}=\mathbf{X},\end{array}$$where $$\begin{array}{c}\text{(32)}& \mathbf{G}=\begin{array}{cc}{\mathbf{G}}_{\text{v}}& {\mathbf{G}}_{\text{d}}\end{array},\mathbf{X}=\begin{array}{cccccc}{\mathbf{\alpha }}_1& \cdots & {\mathbf{\alpha }}_m& 0& \cdots & 0\end{array},\\ \mathbf{Y}=\begin{array}{cccccc}{\mu }_1{\mathbf{\Phi }}_{\text{s}}{\mathbf{v}}_1& \cdots & {\mu }_m{\mathbf{\Phi }}_{\text{s}}{\mathbf{v}}_m& {\lambda }_{m+1}{\mathbf{\Phi }}_{\text{s}}{\mathbf{w}}_{m+1}& \cdots & {\lambda }_{2n_{\text{h}}}{\mathbf{\Phi }}_{\text{s}}{\mathbf{w}}_{2n_{\text{h}}}\\ {\mathbf{\Phi }}_{\text{s}}{\mathbf{v}}_1& \cdots & {\mathbf{\Phi }}_{\text{s}}{\mathbf{v}}_m& {\mathbf{\Phi }}_{\text{s}}{\mathbf{w}}_{m+1}& \cdots & {\mathbf{\Phi }}_{\text{s}}{\mathbf{w}}_{2n_{\text{h}}}\end{array}.\end{array}$$

In case the number of sensors $n_{\text{s}}$ is greater than or equal to the number of structural modes $n_{\text{h}}$, the proper selection of sensor location will make the sensor modal matrix ${\mathbf{\Phi }}_{\text{s}}$ have a full column rank. Thus, the coefficient matrix $\mathbf{Y}$ in Equation (31) is also a full column rank, and the solution of the equation must exist. That is to say, in theory, the weighting parameters ${\mathbf{\alpha }}_i$ can be arbitrarily chosen (appear in conjugate pairs), and the real feedback gains that meet the pole placement conditions can always be obtained by Equation (31).

According to the discussion in Ref. [53], the minimum control effort can be achieved by choosing the closed-loop eigenvectors ${\mathbf{v}}_i$ equal to the open-loop eigenvectors ${\mathbf{w}}_i$. When the number of control surfaces $n_{\text{c}}$ is less than the number of structural modes $n_{\text{h}}$, Equation (28) shows that it is unlikely to find an ${\mathbf{\alpha }}_i$ that exactly meets this condition. Therefore, in this study, the following formula is used to approximately obtain the ${\mathbf{\alpha }}_i$ so that the closed-loop eigenvectors ${\mathbf{v}}_i$ and the open-loop eigenvectors ${\mathbf{w}}_i$ are as close as possible. $$\begin{array}{}\text{(33)}& {\mathbf{\alpha }}_i={{\mathbf{L}}_i^{\text{T}}{\mathbf{L}}_i}^{-1}{\mathbf{L}}_i^{\text{T}}{\mathbf{w}}_i,\hspace{1em}i=1,\cdots ,m,\end{array}$$where ${\mathbf{L}}_i$ is the value of the transfer function $\mathbf{H}s\mathbf{P}s$ at point ${\mu }_i$.

4.2. Optimal Sensor Placement

The number and placement of sensors are critical to determining the existence of control gains. Based on the locations of FEM nodes where sensors are placed, the sensor modal matrix ${\mathbf{\Phi }}_{\text{s}}$ in Equation (22) is formed from the structure modal matrix ${\mathbf{\Phi }}_{\text{h}}$. Normally, it is easy to select a satisfactory ${\mathbf{\Phi }}_{\text{s}}$ to meet the control requirements through trial and error. However, for the folding wing model, the ${\mathbf{\Phi }}_{\text{s}}$ corresponding to a certain sensor layout varies with the folding angle. A sensor layout suitable for one fixed folding angle may not meet the control requirements at other folding angles.

In order to solve the sensor placement problem of the morphing structure, based on the effective independence method proposed by Kammer [54], this paper generalizes the independence distribution vector to a function form. Then, the problem of optimal sensor placement is settled by optimizing the modified distribution vector in an iterative manner. The optimization strategy of the effective independence method is to quantitatively evaluate all possible sensor locations and iteratively eliminate insignificant locations to obtain the final sensor layout. First, define the Fisher information matrix ${\mathbf{Q}}_{\text{s}}\in {ℝ}^{n_{\text{h}}\times n_{\text{h}}}$ as a performance index of the sensor distribution as $$\begin{array}{}\text{(34)}& {\mathbf{Q}}_{\text{s}}={\mathbf{\Phi }}_{\text{s}}^{\text{T}}{\mathbf{\Phi }}_{\text{s}}.\end{array}$$

Only if the information matrix has a relatively small condition number, the problem of observing target structural modes through the sensor outputs is well-conditioned, and Equation (31) does not appear to be ill-conditioned, and the numerical solution exists. Moreover, a larger value of ${\mathbf{Q}}_{\text{s}}$ means that the signal energy output by sensors is larger, which is beneficial for improving the sensor noise immunity and implementing active control. It is commonly suggested to use the trace norm or determinant as an overall index of ${\mathbf{Q}}_{\text{s}}$. In order to evaluate the contribution of each sensor individually, the independence distribution vector ${\mathbf{e}}_{\text{D}}\in {ℝ}^{n_{\text{s}}}$ is defined as the diagonal elements of the projection matrix formulated by ${\mathbf{\Phi }}_{\text{s}}$, as shown below: $$\begin{array}{}\text{(35)}& {\mathbf{e}}_{\text{D}}=\mathrm{diag}{\mathbf{\Phi }}_{\text{s}}{{\mathbf{\Phi }}_{\text{s}}^{\text{T}}{\mathbf{\Phi }}_{\text{s}}}^{-1}{\mathbf{\Phi }}_{\text{s}}^{\text{T}}.\end{array}$$

Kammer [54] pointed out that each element in ${\mathbf{e}}_{\text{D}}$ represents the fractional contribution of the corresponding sensor location to the linear independence of modes in ${\mathbf{\Phi }}_{\text{s}}$ so that removing sensor locations with smaller independence contributions can maintain the determinant of the information matrix. Poston and Tolson [55] later demonstrated that the $i$-th element $e_{\text{D}i}$ of the vector ${\mathbf{e}}_{\text{D}}$ has the following explicit expression: $$\begin{array}{}\text{(36)}& e_{\text{D}i}=\frac{\det {\mathbf{Q}}_{\text{s}}-\det {\mathbf{Q}}_{\text{s}i}}{\det {\mathbf{Q}}_{\text{s}}},\hspace{1em}i=1,\cdots ,n_{\text{s}},\end{array}$$where ${\mathbf{Q}}_{\text{s}i}$ is the information matrix constructed after removing the $i$-th row of ${\mathbf{\Phi }}_{\text{s}}$. From the above formula, it can be seen that $0\le e_{\text{D}i}\le 1$, and $e_{\text{D}i}$ equal to 0 indicates that the determinant of the information matrix remains unchanged after removing the corresponding sensor location, and $e_{\text{D}i}$ equal to 1 indicates that the information matrix is singular after removing the corresponding sensor location. Generally, removing the sensor location with the smallest independence distribution value will maximize the determinant of the information matrix while reducing the number of sensors.

In the folding wing model, ${\mathbf{\Phi }}_{\text{s}}$ is a function of the folding angle $\theta$. In order to apply the effective independence method to the morphing structure, the independence distribution vector defined in Equation (35) should also generalize to the function ${\mathbf{e}}_{\text{D}}\theta$. This paper adopts a harmonic mean fashion to construct the mean value ${\overline{\mathbf{e}}}_{\text{D}}$ of the function ${\mathbf{e}}_{\text{D}}\theta$ as an integrated performance index for the sensor placement, where the $i$-th element ${\overline{e}}_{\text{D}i}$ is calculated as $$\begin{array}{}\text{(37)}& {\overline{e}}_{\text{D}i}=1-{\theta }_{\max }-{\theta }_{\min }{{\int }_{{\theta }_{\min }}^{{\theta }_{\max }}\frac{1}{1-e_{\text{D}i}\theta }d\theta }^{-1},\hspace{1em}i=1,\cdots ,n_{\text{s}},\end{array}$$where ${\theta }_{\min }$ and ${\theta }_{\max }$ are the lower and upper bounds of the folding angle range, respectively. Compared with the arithmetic mean, the above formula is more biased toward larger values in the function curve. This is because the independence distribution close to 1 means that the corresponding sensor location is of great significance; the above index ensures that as long as the function $e_{\text{D}i}\theta$ has a value close to 1, then its mean value ${\overline{e}}_{\text{D}i}$ is also close to 1. In practice, the structural models are established at $n_{\text{t}}$ discrete folding angles, and the definite integral in Equation (37) is approximated by a finite sum as $$\begin{array}{}\text{(38)}& {\overline{e}}_{\text{D}i}=1-n_{\text{t}}{\sum _{l=1}^{n_{\text{t}}}\frac{1}{1-e_{\text{D}i,l}}}^{-1},\hspace{1em}i=1,\cdots ,n_{\text{s}},\end{array}$$where $e_{\text{D}i,l}$ represents the $i$-th element of the independence distribution vector at the $l$-th folding angle. To ensure numerical accuracy, it is recommended to increase the number of modal matrices ${\mathbf{\Phi }}_{\text{s}}$ by interpolating the coherent local models.

The procedure for optimal sensor placement based on the effective independence method is as follows: (1) select an initial candidate set of sensor locations and construct independence distribution vectors at each folding angle; (2) construct the integrated distribution vector according to Equation (38); (3) delete the sensor location with the smallest value in the integrated distribution vector and then reconstruct it; and (4) repeat the above process, deleting one location each time until the required number of sensors is reached.

5. Numerical Results and Discussions

5.1. Validation of the Interpolation-Based Modeling

In this subsection, the effectiveness of the developed interpolation-based modeling method is validated, and the LPV baseline model for control synthesis and response prediction is established. As mentioned above, the procedure of interpolation-based modeling includes two steps: the structure modal matching and the treatment of the RFA coefficient matrices.

In the simulation, a set of 61 local aeroelastic models are generated according to fixed folding angles, ranging from 0 deg to 120 deg with 2 deg intervals. The initial aeroelastic model shown in Equation (5) includes the structural matrices and the GAF matrices, and the coherency of these matrices cannot be guaranteed. To verify the proposed modal matching method, orders of the first 6 structural modes at each parameter point are randomly assigned, as shown in Figure 4(a). The MAC value in Figure 4(b) reveals that the mode shapes between adjacent folding angles have extremely low linear correlations. Hence, all the system matrices in Equation (5) exhibit discontinuities. After implementing the modal matching algorithm described in Section 3.1, it can be seen from Figure 5 that the scrambled structural modes are successively aligned. Thus, the obtained coherent structure modal matrices are continuously dependent on the folding angle.

[figures omitted; refer to PDF]

[figures omitted; refer to PDF]

Figure 6 exhibits the evolution of the first ten mode branches. It is clear that the structural modes are significantly affected by the folding angle, and the proposed modal matching algorithm exactly tracks each mode branch. Complex mode crossing phenomena occur in the high-order modes. Even so, the algorithm still gives the correct matching result.

[figure omitted; refer to PDF]

In the following, retain the first 6 matched structural modes and continue to construct the state-space aeroelastic model of the folding wing. At the incompressible flow condition, the GAF matrices are computed at 16 reduced frequencies in the range of 0 to 1.5. Then, the MS method is applied to fit the GAF matrices to the RFA expression shown in Equation (13). The normal RFA is performed independently at each folding angle, and the aerodynamic lag roots required in the MS method are obtained by the following empirical formula: $$\begin{array}{}\text{(39)}& R_i=-1.7k_{\max }{\frac{i}{n_{\text{a}}+1}}^2,\hspace{1em}i=1,\cdots ,n_{\text{a}},\end{array}$$where the number of aerodynamic lag roots $n_{\text{a}}$ is set to 6 and the maximum reduced frequency $k_{\max }$ is 1.5. In the MS procedure, a constant initial matrix $\mathbf{D}$ and 300 iterations are used for all parameter points.

After the normal RFA procedure is completed, the coefficient matrices are adjusted to a coherent form according to the method described in Section 3.2. All results shown below use the matched structural modes, and the verification work focuses on the influence of the coherent and normal RFA representations on the modeling results. As shown in Figure 7, the directly computed GAF data (solid lines) vary smoothly with the folding angle. The normal RFA (open triangles) and the coherent RFA (filled triangles) are in good agreement with the directly computed GAF data. Although the normal RFA gives accurate fitting results at scattered parameter points, the coefficient matrices $\mathbf{D}$ and $\mathbf{E}$ obtained by normal RFA have dramatic jumps as the parameter varies (see Figure 8(a)), which makes the aeroelastic state-space interpolation impossible. By comparison, it can be concluded from Figure 8(b) that the proposed method for coherent RFA representation successfully adjusts the RFA coefficient matrices to a coherent form without changing the fitting results. The continuously varying system matrices pave the way for state-space model interpolation.

[figures omitted; refer to PDF]

[figures omitted; refer to PDF]

Figure 9 is the comparison of the frequency response functions of the interpolated model and the exact model. In the figure, the flight conditions are fixed at flow speed of 30 m/s and sea level, and the system model at $\theta =65$ deg is obtained by linear interpolation of the local models at $\theta =64$ deg and $\theta =66$ deg. The system inputs are the deflections of the outboard and inboard control surfaces as well as the vertical gust disturbance, and the outputs are the first two modal displacements of the structure (numbered according to the mode branches in Figure 6). It can be seen that the coherent system matrices construct an accurate interpolated model. However, for the incoherent system matrices constructed by unprocessed RFA coefficient matrices, the interpolated model here brings serious modeling errors, as shown by red dashed lines in Figure 9.

[figure omitted; refer to PDF]

Using the parameterized modeling approach proposed in this paper, an aeroelastic LPV model with two scheduling parameters is established as the baseline model for the subsequent simulation and control implementation. Here, the folding angle and the flow speed are taken as the scheduling parameters. In order to generate the local models, the two-dimensional parameter space is divided into regular grids, in which the folding angle parameter is still in the range of 0 deg to 120 deg with 2 deg intervals, and the flow speed is in the range of 10 m/s to 60 m/s with 1 m/s intervals. Figure 10 shows the frequency response curves for varying folding angles and flow speeds, in which the system input and output are the vertical gust disturbance and the first structure modal displacement, respectively.

[figures omitted; refer to PDF]

5.2. Implementation of the Receptance-Based Control

The controller is designed to work at subcritical conditions, so only parameters below the flutter speed need to be considered. The first step in control design is to obtain the eigenvalues and eigenvectors of the open-loop system. Since the GAF matrices have been converted to the RFA representation, solving the characteristic equation of Equation (24) is equivalent to solving the eigenvalues of the state-space matrix $\mathbf{A}$ in Equation (18). Similar to the structure modal data, the initial eigenvalues at discrete parameter points may not have a coherent order. Obviously, the structure modal matching algorithm in Section 3.1 is also applicable to the aeroelastic systems, and it only needs to replace the natural frequency $\omega$ and the mode shape $\mathbf{\phi }$ in Equation (9) with the eigenvalue $\lambda$ and the eigenvector $\mathbf{w}$ of the aeroelastic system. Through the modal matching procedure, the eigenvalues that vary smoothly with the parameter can be obtained. This is useful to compute a continuous parameterized controller through uniform pole placement.

According to the number of structural modes, 6 pairs of complex conjugate eigenvalues are calculated. Based on this, the flutter characteristics of the folding wing within the folding angle range are also obtained. As an example for $\theta =90$ deg, Figure 11 shows the root locus of the matched eigenvalues with respect to the flow speed. The flutter speed is identified by the critical point at which the real part of an eigenvalue trajectory changes from negative to positive. At $\theta =90$ deg, the 2nd mode firstly becomes unstable, and the corresponding flutter speed $U_{\text{f}}=45.1\text{m}/\text{s}$. Figure 12 is a contour map where the real part of the eigenvalues is equal to zero on the two-dimensional parameter plane. In the entire folding angle range, only the 2nd and 4th modes have eigenvalues with nonnegative real parts. Combining the critical contour lines of these two modes, the flutter boundary and the stable region in the parameter space can be obtained.

[figure omitted; refer to PDF]

The sensor layout should be determined before calculating the control gains. In the current example, the necessary condition for the existence of the receptance controller is that the number of sensors is greater than or equal to 6. Therefore, we intend to find 7 sensor locations that maximize the information matrix ${\mathbf{Q}}_{\text{s}}$ defined in Section 4.2. To apply the effective independence method, the initial candidate set of sensor locations is selected as the normal direction of the structural nodes on the inboard and outboard wings. The 7 optimal locations are iteratively retained from 388 initial locations. Figure 13 shows the final sensor layouts obtained by the standard method (fixed configuration) and the modified method (parameter-varying version). Figure 13(a) shows the optimal locations at fixed configurations for $\theta =0$ deg and $\theta =60$ deg, and Figure 13(b) shows the optimal locations obtained by integrating all folding angles.

[figures omitted; refer to PDF]

The information matrices related to the above three sensor layouts are compared and shown in Figures 14 and 15. As shown in Figure 14, the sensor layout obtained at fixed $\theta =0$ deg causes a large condition number of the information matrix at 21.3 deg; the sensor layout obtained at fixed $\theta =60$ deg also encounters the same problem at 104.3 deg. Either of these two sensor layouts will lead to ill-conditioned solutions of Equation (31), and the resulting control gains will produce high control effort or invalid control at certain folding angles. As a comparison, the modified method proposed in this paper provides an improved sensor layout, and the condition number of the information matrix has a relatively small value in the entire folding angle range. Figure 15 also shows the determinant of the information matrix varying with the folding angle. It can be seen that under fixed folding angles of $\theta =0$ deg and $\theta =60$ deg, the sensor layouts obtained by the standard method are optimal at the each given folding angle, but it cannot guarantee performances for other angles. The modified method does not produce the maximum determinant at most folding angles, but the obtained sensor layout achieves good comprehensive performance in the entire folding angle range.

[figure omitted; refer to PDF]

The next step in control design is to manually assign the eigenvalues of the closed-loop system. The system response induced by gust disturbances is dominated by low-frequency modes. Therefore, the first 4 eigenvalues of the open-loop system are intended to be assigned to specified values, while the 5th and 6th eigenvalues remain unchanged. The assigned closed-loop eigenvalues are specified by the known open-loop eigenvalues plus the real and imaginary part increments, so there are 8 increments that need to be given at each parameter point. In the two-dimensional parameter space, the vector $\Delta$ composed of the increments of the open-loop eigenvalues is uniformly given by the following polynomial: $$\begin{array}{}\text{(40)}& \Delta x,y=\sum _{i+j=0}^n{\mathbf{p}}_{ij}{x}^i{y}^j,\end{array}$$where $x$ and $y$ are the normalized parameters of $\theta$ and $U_{\infty }$, ${\mathbf{p}}_{ij}$ is the coefficient vector, and $n$ is the degree of the polynomial. The way of specifying the closed-loop eigenvalues by an explicit expression can make the controller change continuously with the parameters. Set $n=3$, and then all undetermined coefficients constitute 80 free variables, and these variables will determine the final control gains and closed-loop models. The influence of pole placement on disturbance rejection can be evaluated by calculating the $H_2$-norm or $H_{\infty }$ -norm of the closed-loop system. Unfortunately, analytically assigning poles to minimize the system norm is still an open problem, so this study uses the optimization method to obtain the coefficients in Equation (40). The optimization is implemented by using the Nelder-Mead simplex algorithm [56], which is a direct search method that does not use gradients. The objective function is set as the average value of the $H_2$-norm of all local closed-loop systems within the flutter boundary.

According to the above steps, the control gains are calculated one by one at each design parameter point. Based on the interpolation of the local controllers, a parameterized gust alleviation controller for the folding wing is finally established. Figure 16 compares the $H_2$-norm distribution of the open-loop and closed-loop systems in the parameter space under the vertical gust disturbance. The figure uses different colors to indicate the system norm. The $H_2$-norm of the open-loop system at the flutter boundary approaches infinity, so only the gust alleviation results within the flutter boundary are investigated. In general, increasing the folding angle can passively reduce the influence of the vertical gust disturbance. After the active control system works, the response energy of the closed-loop system is significantly reduced, and the controller shows a powerful control capability within a wide parameter-varying range.

[figures omitted; refer to PDF]

Figure 17 shows the frequency responses of the open-loop and closed-loop systems at different folding angles. In each subfigure, the folding angle is fixed, and a series of curves represent the frequency response amplitude at different flow speeds from 10 m/s to the flutter speed $U_{\text{f}}$. The open-loop and closed-loop frequency responses are represented by gray and colored curves, respectively. It can be seen that changing the first 4 eigenvalues of the system can effectively reduce the low-frequency responses. The first-order bending vibration of the wing is greatly reduced by the parameterized controller.

[figure omitted; refer to PDF]

In order to further verify the controller performance, the time-domain dynamic responses of the folding wing to gust excitations are calculated. Figures 18 and 19 show the dynamic responses of the folding wing to the 1-cos discrete gust under two sets of simulation parameters. The frequency of the 1-cos gust is set to 10 Hz, and the maximum gust velocity is 1 m/s. The simulation results show that the displacement and wing-root bending moment of the folding wing are significantly reduced due to the driving of the inboard and outboard control surfaces. In the flight condition of $\theta =0$ deg and $U_{\infty }=15\text{m}/\text{s}$, the maximum value of the wing-tip displacement at node 641 is reduced by 62.00%, while the maximum value of the wing-root bending moment is reduced by 43.88%. In the condition of $\theta =90$ deg and $U_{\infty }=30\text{m}/\text{s}$, the maximum values of the displacement and the bending moment are reduced by 79.99% and 60.34%, respectively. Table 1 lists the detailed data of the open-loop and closed-loop responses at four different folding angles.

[figures omitted; refer to PDF]

[figures omitted; refer to PDF]

Table 1

Maximum values of the wing-tip displacement ($D_n$), wing-root bending moment ($M_x$), and outboard flap deflection ($\delta$) in discrete gust responses.

$\theta$ (deg)	Variable	Maximal value of the variable
		$U_{\infty }=15\text{m}/\text{s}$			$U_{\infty }=30\text{m}/\text{s}$
		Open-loop	Closed-loop	Reduction of the maximal value	Open-loop	Closed-loop	Reduction of the maximal value
0	$D_n$ (mm)	35.21	13.38	62.00%	77.57	35.08	54.77%
	$M_x$ (N∙mm)	340.86	191.28	43.88%	731.16	618.80	15.37%
	$\delta$ (deg)	—	14.62	—	—	12.18	—
50	$D_n$ (mm)	27.00	7.99	70.40%	Inf.	12.21	—
	$M_x$ (N∙mm)	292.28	117.68	59.74%	Inf.	294.88	—
	$\delta$ (deg)	—	10.98	—	—	9.21	—
90	$D_n$ (mm)	16.57	3.11	81.21%	28.24	5.65	79.99%
	$M_x$ (N∙mm)	205.04	48.23	76.48%	338.76	134.35	60.34%
	$\delta$ (deg)	—	11.31	—	—	6.86	—
120	$D_n$ (mm)	11.82	2.88	75.65%	19.26	4.79	75.15%
	$M_x$ (N∙mm)	125.75	14.93	88.12%	199.03	19.75	90.08%
	$\delta$ (deg)	—	8.73	—	—	4.66	—

For the continuous gust control, Table 2 shows an overview of the folding wing responses to the Dryden gust, as well as the statistical comparisons of the open-loop and closed-loop responses at different folding angles. The scale of the Dryden gust is set to 5 m, and the root mean square (RMS) value of the gust velocity is 0.5 m/s. The RMS values of the wing-tip displacement, the wing-root bending moment, and the outboard control surface deflection are listed in Table 2.

Table 2

RMS values of the wing-tip displacement ($D_n$), wing-root bending moment ($M_x$), and outboard flap deflection ($\delta$) in continuous gust responses.

$\theta$ (deg)	Variable	RMS value of the variable
		$U_{\infty }=15\text{m}/\text{s}$			$U_{\infty }=30\text{m}/\text{s}$
		Open-loop	Closed-loop	Reduction of the RMS value	Open-loop	Closed-loop	Reduction of the RMS value
0	$D_n$ (mm)	22.47	5.90	73.75%	67.73	16.67	75.39%
	$M_x$ (N∙mm)	272.53	111.24	59.18%	634.32	303.91	52.09%
	$\delta$ (deg)	—	5.82	—	—	4.77	—
50	$D_n$ (mm)	17.01	3.95	76.80%	Inf.	7.14	—
	$M_x$ (N∙mm)	217.84	67.33	69.09%	Inf.	160.25	—
	$\delta$ (deg)	—	4.18	—	—	3.94	—
90	$D_n$ (mm)	9.60	1.50	84.33%	15.43	2.44	84.18%
	$M_x$ (N∙mm)	125.33	36.05	71.24%	209.29	114.84	45.13%
	$\delta$ (deg)	—	4.80	—	—	4.02	—
120	$D_n$ (mm)	5.72	1.23	78.53%	9.29	3.01	67.56%
	$M_x$ (N∙mm)	58.12	6.20	89.33%	102.92	10.04	90.25%
	$\delta$ (deg)	—	3.34	—	—	2.00	—

Through the fixed parameter simulations, it is clearly seen that the controller is valid at all parameter points and has good alleviation effects in both the discrete and continuous gusts. Next, the time-varying system simulation and control for the folding wing system are performed. As shown in Figure 20, during the simulation time of 60 seconds, both the folding angle and the flow speed change slowly and uniformly, where the folding angle changes from 0 deg to 120 deg and the flow speed changes from 10 m/s to 30 m/s. The time-varying system simulation shows that the structural vibration and additional loads are alleviated in each time period. In conclusion, the above results verify that the proposed parameterized controller is valid for a wide range of parameters. Under the assumption of slow parameter variation, extending the receptance method to the parameter-varying system is able to design a reliable and effective gust alleviation controller.

[figures omitted; refer to PDF]

6. Conclusions

To efficiently investigate the aeroelasticity and control of a folding wing, this paper presents an interpolation-based modeling strategy for parameterized aeroelastic systems. The key steps involved are the structure modal matching and the manipulation of the RFA coefficient matrices. The main advantage of the proposed method is that the coherent local models are obtained before constructing the aeroelastic state-space matrices, which makes the coherent state-space representation much easier. Based on the developed modeling strategy, the LPV model of the folding wing system is constructed within the local modeling framework. Numerical examples demonstrate that the interpolation of incoherent local models brings serious modeling errors, while the proposed coherent representation method gives accurate modeling results. Next, the aeroelastic control for the folding wing under various flight conditions and structural configurations is studied. For this purpose, the original receptance-based control method for fixed configuration is extended to the parameter-varying system, and a modified version of the effective independence method is derived to select an optimal sensor layout suitable for all folding angles. The simulation results show that the designed parameterized controller for gust alleviation achieves satisfactory closed-loop performance in the given parameter space.

Acknowledgments

This work was supported by the National Natural Science Foundation of China under Grant No. 11472128.

References

[1] T. A. Weisshaar, "Morphing aircraft systems: historical perspectives and future challenges," Journal of Aircraft, vol. 50 no. 2, pp. 337-353, DOI: 10.2514/1.C031456, 2013.

[2] S. Barbarino, O. Bilgen, R. M. Ajaj, M. I. Friswell, D. J. Inman, "A review of morphing aircraft," Journal of Intelligent Material Systems and Structures, vol. 22 no. 9, pp. 823-877, DOI: 10.1177/1045389X11414084, 2011.

[3] R. M. Ajaj, C. S. Beaverstock, M. I. Friswell, "Morphing aircraft: the need for a new design philosophy," Aerospace Science and Technology, vol. 49, pp. 154-166, DOI: 10.1016/j.ast.2015.11.039, 2016.

[4] M. H. Love, P. S. Zink, R. L. Stroud, D. R. Bye, C. Chase, "Impact of actuation concepts on morphing aircraft structures," Proceedings of the 45th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics and Materials Conference, .

[5] I. Wang, "Aeroelastic and flight dynamics analysis of folding wing systems, [Ph.D. thesis]," 2013.

[6] D. H. Lee, T. A. Weisshaar, "Aeroelastic studies on a folding wing configuration," Proceedings of the 46th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics and Materials Conference, .

[7] D. H. Lee, P. C. Chen, "Nonlinear aeroelastic studies on a folding wing configuration with free-play hinge nonlinearity," Proceedings of the 47th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics and Materials Conference, .

[8] I. Wang, S. C. Gibbs, E. H. Dowell, "Aeroelastic model of multisegmented folding wings: theory and experiment," Journal of Aircraft, vol. 49 no. 3, pp. 911-921, DOI: 10.2514/1.C031589, 2012.

[9] Y. Zhao, H. Hu, "Parameterized aeroelastic modeling and flutter analysis for a folding wing," Journal of Sound and Vibration, vol. 331 no. 2, pp. 308-324, DOI: 10.1016/j.jsv.2011.08.028, 2012.

[10] D. Tang, E. H. Dowell, "Theoretical and experimental aeroelastic study for folding wing structures," Journal of Aircraft, vol. 45 no. 4, pp. 1136-1147, DOI: 10.2514/1.32754, 2008.

[11] P. J. Attar, D. Tang, E. H. Dowell, "Nonlinear aeroelastic study for folding wing structures," AIAA Journal, vol. 48 no. 10, pp. 2187-2195, DOI: 10.2514/1.44868, 2010.

[12] Y. Zhao, H. Hu, "Prediction of transient responses of a folding wing during the morphing process," Aerospace Science and Technology, vol. 24 no. 1, pp. 89-94, DOI: 10.1016/j.ast.2011.09.001, 2013.

[13] W. Hu, Z. Yang, Y. Gu, "Aeroelastic study for folding wing during the morphing process," Journal of Sound and Vibration, vol. 365, pp. 216-229, DOI: 10.1016/j.jsv.2015.11.043, 2016.

[14] W. Hu, Z. Yang, Y. Gu, X. Wang, "The nonlinear aeroelastic characteristics of a folding wing with cubic stiffness," Journal of Sound and Vibration, vol. 400, pp. 22-39, DOI: 10.1016/j.jsv.2017.04.002, 2017.

[15] M. L. Verstraete, B. A. Roccia, D. T. Mook, S. Preidikman, "A co-simulation methodology to simulate the nonlinear aeroelastic behavior of a folding-wing concept in different flight configurations," Nonlinear Dynamics, vol. 98 no. 2, pp. 907-927, DOI: 10.1007/s11071-019-05234-9, 2019.

[16] T. M. Seigler, D. A. Neal, J. S. Bae, D. J. Inman, "Modeling and flight control of large-scale morphing aircraft," Journal of Aircraft, vol. 44 no. 4, pp. 1077-1087, DOI: 10.2514/1.21439, 2007.

[17] B. Yan, Y. Li, P. Dai, S. Liu, "Aerodynamic analysis, dynamic modeling, and control of a morphing aircraft," Journal of Aerospace Engineering, vol. 32 no. 5, article 04019058,DOI: 10.1061/(ASCE)AS.1943-5525.0001047, 2019.

[18] R. Huang, Z. Yang, X. Yao, Y. Zhao, H. Hu, "Parameterized modeling methodology for efficient aeroservoelastic analysis of a morphing wing," AIAA Journal, vol. 57 no. 12, pp. 5543-5552, DOI: 10.2514/1.J058211, 2019.

[19] D. H. Baldelli, D. H. Lee, R. S. S. Pena, B. Cannon, "Modeling and control of an aeroelastic morphing vehicle," Journal of Guidance, Control, and Dynamics, vol. 31 no. 6, pp. 1687-1699, DOI: 10.2514/1.35445, 2008.

[20] T. Yue, L. Wang, J. Ai, "Gain self-scheduled _H_ ∞ control for morphing aircraft in the wing transition process based on an LPV model," Chinese Journal of Aeronautics, vol. 26 no. 4, pp. 909-917, DOI: 10.1016/j.cja.2013.06.004, 2013.

[21] A. Marcos, G. J. Balas, "Development of linear-parameter-varying models for aircraft," Journal of Guidance, Control, and Dynamics, vol. 27 no. 2, pp. 218-228, DOI: 10.2514/1.9165, 2004.

[22] H. Pfifer, S. Hecker, "Generation of optimal linear parametric models for LFT-based robust stability analysis and control design," IEEE Transactions on Control Systems Technology, vol. 19 no. 1, pp. 118-131, DOI: 10.1109/TCST.2010.2076329, 2011.

[23] G. Ferreres, "Computation of a flexible aircraft LPV/LFT model using interpolation," IEEE Transactions on Control Systems Technology, vol. 19 no. 1, pp. 132-139, DOI: 10.1109/TCST.2010.2078510, 2011.

[24] M. Lovera, M. Bergamasco, F. Casella, "LPV modelling and identification: an overview," Robust Control and Linear Parameter Varying Approaches,DOI: 10.1007/978-3-642-36110-4_1, 2013.

[25] C. Poussot-Vassal, C. Roos, "Generation of a reduced-order LPV/LFT model from a set of large-scale MIMO LTI flexible aircraft models," Control Engineering Practice, vol. 20 no. 9, pp. 919-930, DOI: 10.1016/j.conengprac.2012.06.001, 2012.

[26] C. Roos, "Generation of flexible aircraft LFT models for robustness analysis," IFAC Proceedings Volumes, vol. 42 no. 6, pp. 349-354, DOI: 10.3182/20090616-3-IL-2002.00060, 2009.

[27] C. Poussot-Vassal, F. Demourant, "Dynamical medium (large)-scale model reduction and interpolation with application to aircraft systems," The French Aerospace Lab, vol. 4, 2012.

[28] J. Theis, B. Takarics, H. Pfifer, G. Balas, H. Werner, "Modal matching for LPV model reduction of aeroservoelastic vehicles," Proceedings of the AIAA Atmospheric Flight Mechanics Conference, .

[29] A. K. Al-Jiboory, G. Zhu, S. S. M. Swei, W. Su, N. T. Nguyen, "LPV modeling of a flexible wing aircraft using modal alignment and adaptive gridding methods," Aerospace Science and Technology, vol. 66, pp. 92-102, DOI: 10.1016/j.ast.2017.03.009, 2017.

[30] M. Lovera, G. Mercere, "Identification for gain-scheduling: a balanced subspace approach," Proceedings of the 2007 American Control Conference, .

[31] F. Ferranti, Y. Rolain, "A local identification method for linear parameter-varying systems based on interpolation of state-space matrices and least-squares approximation," Mechanical Systems and Signal Processing, vol. 82, pp. 478-489, DOI: 10.1016/j.ymssp.2016.05.037, 2017.

[32] J. de Caigny, R. Pintelon, J. F. Camino, J. Swevers, "Interpolated modeling of LPV systems," IEEE Transactions on Control Systems Technology, vol. 22 no. 6, pp. 2232-2246, DOI: 10.1109/TCST.2014.2300510, 2014.

[33] Y. M. Ram, J. E. Mottershead, "Receptance method in active vibration control," AIAA Journal, vol. 45 no. 3, pp. 562-567, DOI: 10.2514/1.24349, 2007.

[34] Y. M. Ram, J. E. Mottershead, "Multiple-input active vibration control by partial pole placement using the method of receptances," Mechanical Systems and Signal Processing, vol. 40 no. 2, pp. 727-735, DOI: 10.1016/j.ymssp.2013.06.008, 2013.

[35] K. V. Singh, L. A. McDonough, R. Kolonay, J. E. Cooper, "Receptance-based active aeroelastic control using multiple control surfaces," Journal of Aircraft, vol. 51 no. 1, pp. 335-342, DOI: 10.2514/1.C032183, 2014.

[36] H. Liu, X. Gao, X. Wang, "Parametric active aeroelastic control of a morphing wing using the receptance method," Journal of Fluids and Structures, vol. 98, article 103098,DOI: 10.1016/j.jfluidstructs.2020.103098, 2020.

[37] E. Albano, W. P. Rodden, "A doublet-lattice method for calculating lift distributions on oscillating surfaces in subsonic flows," AIAA Journal, vol. 7 no. 2, pp. 279-285, DOI: 10.2514/3.5086, 1969.

[38] Y. Zhao, R. Huang, Advanced Aeroelasticity and Control, 2015.

[39] R. L. Harder, R. N. Desmarais, "Interpolation using surface splines," Journal of Aircraft, vol. 9 no. 2, pp. 189-191, DOI: 10.2514/3.44330, 1972.

[40] M. Karpel, B. Moulin, P. C. Chen, "Dynamic response of aeroservoelastic systems to gust excitation," Journal of Aircraft, vol. 42 no. 5, pp. 1264-1272, DOI: 10.2514/1.6678, 2005.

[41] M. Karpel, "Design for active flutter suppression and gust alleviation using state-space aeroelastic modeling," Journal of Aircraft, vol. 19 no. 3, pp. 221-227, DOI: 10.2514/3.57379, 1982.

[42] Z. Alkhoury, M. Petreczky, G. Mercere, "Comparing global input-output behavior of frozen-equivalent LPV state-space models," IFAC-PapersOnLine, vol. 50 no. 1, pp. 9766-9771, DOI: 10.1016/j.ifacol.2017.08.2182, 2017.

[43] R. Burkard, M. Dell'Amico, S. Martello, Assignment Problems,DOI: 10.1137/1.9780898717754, 2009.

[44] H. W. Kuhn, "The Hungarian method for the assignment problem," Naval Research Logistics, vol. 2 no. 1–2, pp. 83-97, DOI: 10.1002/nav.3800020109, 1955.

[45] R. Jonker, A. Volgenant, "A shortest augmenting path algorithm for dense and sparse linear assignment problems," Computing, vol. 38 no. 4, pp. 325-340, DOI: 10.1007/BF02278710, 1987.

[46] E. Livne, "Alternative approximations for integrated control/structure aeroservoelastic synthesis," AIAA Journal, vol. 31 no. 6, pp. 1100-1108, DOI: 10.2514/3.49052, 1993.

[47] F. M. Hoblit, Gust Loads on Aircraft: Concepts and Applications,DOI: 10.2514/4.861888, 1988.

[48] Y. Zhao, C. Yue, H. Hu, "Gust load alleviation on a large transport airplane," Journal of Aircraft, vol. 53 no. 6, pp. 1932-1946, DOI: 10.2514/1.C033713, 2016.

[49] P. Apkarian, P. Gahinet, "A convex characterization of gain-scheduled H ∞ controllers," IEEE Transactions on Automatic Control, vol. 40 no. 5, pp. 853-864, DOI: 10.1109/9.384219, 1995.

[50] P. Apkarian, P. Gahinet, G. Becker, "Self-scheduled H ∞ control of linear parameter-varying systems: a design example," Automatica, vol. 31 no. 9, pp. 1251-1261, DOI: 10.1016/0005-1098(95)00038-X, 1995.

[51] D. J. Stilwell, W. J. Rugh, "Stability preserving interpolation methods for the synthesis of gain scheduled controllers," Automatica, vol. 36 no. 5, pp. 665-671, DOI: 10.1016/S0005-1098(99)00193-4, 2000.

[52] H. J. Hassig, "An approximate true damping solution of the flutter equation by determinant iteration," Journal of Aircraft, vol. 8 no. 11, pp. 885-889, DOI: 10.2514/3.44311, 1971.

[53] B. Mokrani, F. Palazzo, J. E. Mottershead, S. Fichera, "Multiple-input multiple-output experimental aeroelastic control using a receptance-based method," AIAA Journal, vol. 57 no. 7, pp. 3066-3077, DOI: 10.2514/1.J057855, 2019.

[54] D. C. Kammer, "Sensor placement for on-orbit modal identification and correlation of large space structures," Journal of Guidance, Control, and Dynamics, vol. 14 no. 2, pp. 251-259, DOI: 10.2514/3.20635, 1991.

[55] W. L. Poston, R. H. Tolson, "Maximizing the determinant of the information matrix with the effective independence method," Journal of Guidance, Control, and Dynamics, vol. 15 no. 6, pp. 1513-1514, DOI: 10.2514/3.11419, 1992.

[56] J. C. Lagarias, J. A. Reeds, M. H. Wright, P. E. Wright, "Convergence properties of the Nelder-Mead simplex method in low dimensions," SIAM Journal on Optimization, vol. 9 no. 1, pp. 112-147, DOI: 10.1137/S1052623496303470, 1998.