Nomenclature
- AFS
-
active flutter suppression
- AFW
-
active flexible wing
- AAW
-
active aeroelastic wing
- ADRC
-
active disturbance rejection control
- BACT
-
benchmark active control technology
- MIMO
-
multiple input multiple output
-
$m$
-
total mass of the wing
-
${I_\alpha }$
-
rotational inertia of the wing to the elastic axis
-
${s_{h\alpha }}$
-
mass static moment of the wing to the elastic axis
-
${s_{h\beta }}$
-
mass static moment of the control surface to its rotation axis
-
${I_\beta }$
-
moment of inertia of the control surface to its rotation axis
-
$h$
-
sinking and floating displacement
-
$\alpha $
-
pitching angle
-
$\beta $
-
control surface deflection angle
-
${d_h}$
-
sinking and floating damping coefficient
-
${d_\alpha }$
-
pitching damping coefficient
-
${k_h}$
-
sinking and floating stiffness coefficient
-
${k_\alpha }$
-
pitching stiffness coefficient
-
${s_p}$
-
wing span
-
$L$
-
aerodynamic lift
-
${T_\alpha }$
-
aerodynamic moment to the elastic axis
-
$\delta $
-
control surface deflection instruction
-
$k$
-
proportional coefficient
-
$\omega $
-
natural frequency
-
$\varsigma $
-
damping ratio
-
$\boldsymbol{\varPhi} $
-
pseudo partial derivative
-
$\eta ,\varepsilon $
-
step length
-
$\lambda ,\mu $
-
regularisation coefficient
-
${T_s}$
-
sampling period
-
$R$
-
weighting coefficient
-
${\boldsymbol{Q}}$
-
weighting matrix
1.0 Introduction
Active flutter suppression (AFS) is one of the main active control technologies applied in aircraft design at present. It belongs to the comprehensive category of aeroservoelasticity and is one of the research hotspots in the field of aeroservoelasticity. It focuses on increasing the critical flutter velocity without significantly increasing the structural weight of the aircraft [Reference Qian, Huang, Hu and Zhao1, Reference Huang, Hu and Zhao2]. Active flexible wing (AFW) was proposed in the United States in the 1980s, and since then, active flutter suppression has been continuously invested in research and development. In this century, this concept has been developed into active aeroelastic wing (AAW) technology, which has promoted the aeroservoelasticity comprehensive technology to a higher theoretical and technical level. It is a new concept of modern aircraft design.
In recent years, important progress has been made in the research of AFS technology. Zeng et al. [Reference Zeng, Kukreja and Moulin3] developed a feedback control framework based on a wind tunnel experimental model to suppress flutter using robust control laws. Schmidt [Reference Schmidt4] aimed at Lockheed Martin’s BFF06 whole aircraft model, using the method named “identically located force and acceleration” to construct a double loop system, and increase the model damping, and achieved good results in flutter suppression. Huang et al. [Reference Huang, Zhao and Hu5] conducted a wind tunnel experimental study on closed-loop flutter modal identification and AFS of a 3D small aspect ratio wing model. The linear feedback controller designed through pole assignment optimisation effectively suppressed flutter, and explored the on-line closed-loop flutter identification of the wing model. Theis et al. [Reference Theis, Pfifer and Seiler6] proposed a systematic robust control design method for active flutter suppression. It extends the standard four-block mixed sensitivity formulation by a means that targets specific dynamic modes and adds damping. This enables a control design to augment the damping of critical flutter modes with minimal impact on the rigid-body autopilots. Vepa et al. [Reference Vepa and Kwon7] used the linearised aeroelastic model based on the doublet lattice method to initially design the control law family for the active flutter suppression system. Using these preliminary control laws and the nonlinear transonic small disturbance theory, the approximate optimal control law was selected in the transonic domain, and the flutter of a typical wing was successfully suppressed. Yang et al. applied active disturbance rejection control (ADRC) theory to the design of aeroelastic control law and studied the transonic AFS problem of benchmark active control technology (BACT) wing [Reference Yang, Huang, Zhao and Hu8] and 3D elastic wing [Reference Yang, Huang, Zhao and Hu9] considering parameter uncertainty and measurement noise interference. Although AFS technology has been widely used in aircraft design, it still faces the following difficulties.
First, the existing flutter suppression controller structure is complex: The complex structure of the mathematical model determines the complex structure of the controller, and the complex high-order nonlinear system model inevitably leads to the complex high-order nonlinear controller, and the problems of simplification and reduction of the controller and robustness become insurmountable design problems.
Second, there is no systematic design method for multiple input multiple output (MIMO) systems: For flutter suppression, it is impossible to find a general mode to design the controller in the existing control theories and methods. For different aeroservoelastic systems, it can only be analysed on a case-by-case basis, which undoubtedly requires designers to have strong expert experience and knowledge, and cannot be quickly designed in a process, increasing the design cycle.
Third, it is still necessary to rely on the model design control law: Many problems in aviation engineering, such as high-fidelity unsteady aerodynamic model, nonlinear time-variant aircraft model and various uncertain disturbances, need a lot of expert experience to establish the system mechanism model and global mathematical model to solve, and the cost is huge. In addition, for complex systems, due to the complexity of the system itself and various interferences, it is impossible to establish a global mathematical model of the system, even if the local model is not very accurate. Therefore, although the model-based control theory has rich results, it is weak in solving practical complex aeronautical engineering problems.
To solve the above problems, a model-free adaptive control method based on ridge regression is proposed in this paper. The active flutter suppression controller designed by using this method can achieve the purpose of flutter suppression only by using real-time measurement of the online input and output data of the closed-loop controlled system for continuous correction.
2.0 Modeling of aeroelastic system
Figure 1 shows a three-degree-of-freedom binary wing dynamic model with control surface. The equation of motion of the binary wing can be written as
\begin{equation}\left[ {\begin{array}{c@{\quad}c}m & {{s_{h\alpha }}}\\{{s_{h\alpha }}}&{{I_\alpha }}\end{array}} \right]\left\{ {\begin{array}{c}\displaystyle{\ddot h}\\\displaystyle{\ddot \alpha }\end{array}} \right\} + \left[ {\begin{array}{c@{\quad}c}{{d_h}} & 0\\0& {{d_\alpha }}\end{array}} \right]\left\{ {\begin{array}{c}{\dot h}\\{\dot \alpha }\end{array}} \right\} + \left[ {\begin{array}{c@{\quad}c}{{k_h}} &0\\0& {{k_\alpha }}\end{array}} \right] = - \left\{ {\begin{array}{c}{{s_{h\beta }}}\\{{s_{\alpha \beta }}}\end{array}} \right\}\ddot \beta + \left\{ {\begin{array}{{c}}{ - L}\\{{T_\alpha }}\end{array}} \right\}\end{equation}
Three-degree-of-freedom binary wing dynamic model.
where
$m$
is the total mass of the wing;
${I_\alpha }$
is the rotational inertia to the elastic axis;
${s_{h\alpha }}$
is the mass static moment of the wing to the elastic axis;
${s_{h\beta }}$
is the mass static moment of the control surface to its rotation axis;
${s_{\alpha \beta }} = \left( {\bar c - \bar a} \right)b{s_{h\beta }} + {I_\beta }$
, among them,
${I_\beta }$
is the moment of inertia of the control surface to its rotation axis;
$h$
is sinking and floating displacement;
$\alpha $
is pitching displacement;
$\beta $
is the control surface deflection angle;
${d_h}$
is sinking and floating damping coefficient;
${d_\alpha }$
is pitching damping coefficient;
${k_h}$
is sinking and floating stiffness coefficient;
${k_\alpha }$
is pitching stiffness coefficient; the aerodynamic forces acting on the wing with span of
${s_p}$
are calculated using Theodorson theory, among them,
$L$
is aerodynamic lift,
${T_\alpha }$
is the aerodynamic moment to the elastic axis, the specific algorithm is shown in Appendix.
The actuating system reflects the relationship between the control surface deflection instruction and the actual deflection angle of the control surface, which is described by the following second-order system
\begin{equation}\ddot \beta + 2\varsigma \omega \dot \beta + {\omega ^2}\beta \,{=}\,k{\omega ^2}\delta \end{equation}
where
$\delta $
represents control surface deflection instruction;
$k$
is the proportional coefficient;
$\omega $
represents the natural frequency of the system;
$\varsigma $
represents the damping ratio. According to Equations (1) and (2), the motion equation of the wing including the actuating system can be written as
\begin{equation}{\boldsymbol{M}}_{\boldsymbol{s}}\boldsymbol{\ddot{\boldsymbol{q}}}_{\boldsymbol{s}} + {\boldsymbol{D}}_{\boldsymbol{s}}\boldsymbol{\dot{\boldsymbol{q}}}_{\boldsymbol{s}} + {\boldsymbol{K}}_{\boldsymbol{s}}{\boldsymbol{q}}_{\boldsymbol{s}} = {\boldsymbol{F}}_{\boldsymbol{ae}} + {\boldsymbol{G}}_{\boldsymbol{s}}\delta \end{equation}
where
${\boldsymbol{M}}_{\boldsymbol{s}} = \left[ {\begin{array}{c@{\quad}c@{\quad}c}m& {{s_{h\alpha }}}& {{s_{h\beta }}}\\{{s_{h\alpha }}}& {{I_\alpha }}& {{s_{\alpha \beta }}}\\0& 0& 1\end{array}} \right]$
,
${\boldsymbol{D}}_{\boldsymbol{s}} = \left[ {\begin{array}{c@{\quad}c@{\quad}c}{{d_h}} &0 &0\\0 &{{d_\alpha }} &0\\0 &0 &{{d_\beta }}\end{array}} \right]$
,
${\boldsymbol{K}}_{\boldsymbol{s}} = \left[ {\begin{array}{c@{\quad}c@{\quad}c}{{k_h}} &0 &0\\0 &{{k_\alpha }} &0\\0 &0 &{{k_\beta }}\end{array}} \right]$
,
${\boldsymbol{F}}_{\boldsymbol{ae}} = \left[ {\begin{array}{c}{ - L}\\{{T_\alpha }}\\0\end{array}} \right]$
,
${\boldsymbol{G}}_{\boldsymbol{s}} = \left[ {\begin{array}{c}0\\0\\{{g_0}}\end{array}} \right]$
,
${\boldsymbol{q}}_{\boldsymbol{s}} = {\left[ {\begin{array}{c}h\,\, \alpha\,\, \beta \end{array}} \right]^{\text{T}}}$
,
${d_\beta } = 2\varsigma \omega $
,
${k_\beta } = {\omega ^2}$
,
${g_0} = k{\omega ^2}$
.
In order to design the optimal controller for flutter suppression, it is necessary to express the open-loop aeroelastic model in time-domain state space. At this time, in addition to the state variables related to the structure, there are also the constructed aerodynamic state variables. The state space realisation of the aeroelastic equation can be written as
\begin{equation}{\dot{\boldsymbol{X}}} = {\boldsymbol{AX}} + {\boldsymbol{B}}\delta \end{equation}
where
$\boldsymbol{A} = \left[ \begin{array}{c@{\quad}c@{\quad}c}{{\textbf{0}_{3 \times 3}}} &{{\boldsymbol{I}_{3 \times 3}}} &{{\textbf{0}_{3 \times 2}}}\\{ - {\boldsymbol{M}^{ - 1}}\boldsymbol{K}} &{ - {\boldsymbol{M}^{ - 1}}\boldsymbol{D}} &{ - {{\textit{V}}^2}{\boldsymbol{M}^{ - 1}}{\boldsymbol{E}_c}}\\{{\textit{V}}{\boldsymbol{K}_\boldsymbol{a}}} &{\boldsymbol{D}_\boldsymbol{a}} &{{\textit{V}}{\boldsymbol{A}_\boldsymbol{a}}}\end{array} \right]$
,
$\boldsymbol{B} = \left[ {\begin{array}{c}{{\textbf{0}_{3 \times 1}}}\\{{{\boldsymbol{M}}^{ - 1}}{\boldsymbol{G}}_{\boldsymbol{s}}}\\{{\textbf{0}_{2 \times 1}}}\end{array}} \right]$
,
$\boldsymbol{X} = {\left[ {\begin{array}{c}h \,\,\alpha\,\, \beta \,\, {\dot h}\,\, {\dot \alpha }\,\, {\dot \beta}\,\, {{x_{a1}}}\,\, {{x_{a2}}}\end{array}} \right]^{\textrm{T}}}$
, please refer to Appendix for the meaning of specific parameters and the calculation process. For the model-based optimal controller design, measurable output of the controlled object should be used to provide feedback information, while the aerodynamic state variable is unmeasurable, so the outputs selected are sinking and floating displacement of the wing, pitching angle of the wing and deflection angle of the control surface. The state equation and output equation of the controlled system can be written as
\begin{equation}\left\{ \begin{array}{l}{\dot{\boldsymbol{X}}} = {\boldsymbol{AX}} + {\boldsymbol{B}}\delta \\{\boldsymbol{Y}} = {\boldsymbol{CX}}\end{array} \right.\end{equation}
where
${\boldsymbol{C}} = \left[ {\begin{array}{*{20}{c}}{{{\boldsymbol{I}}_{3 \times 3}}}\quad {{\textbf{0}_{3 \times 5}}}\end{array}} \right]$
.
3.0 Design of adaptive filter for structural characteristics of discrete system
For MIMO discrete time system, it can be written in the form of Equation (6)
\begin{equation}\left\{ {\begin{array}{c}{{\boldsymbol{x}}\left( {k + 1} \right) = {\boldsymbol{g}}\left( {{\boldsymbol{x}}\!\left( k \right),{\boldsymbol{u}}\!\left( k \right)} \right)}\\{{\boldsymbol{y}}\left( {k + 1} \right) = {\boldsymbol{h}}\left( {{\boldsymbol{x}}\!\left( k \right),{\boldsymbol{u}}\!\left( k \right)} \right)}\end{array}} \right.\end{equation}
where
${\boldsymbol{g}}$
,
${\boldsymbol{h}}$
represent the unknown state equation and output equation, respectively,
${\boldsymbol{y}}\left( {k + 1} \right) = {\left[ {\begin{array}{*{20}{c}}{{y_1}\left( {k + 1} \right)} \quad \cdots\quad {{y_q}\left( {k + 1} \right)}\end{array}} \right]^{\text{T}}}$
is the system output,
${\boldsymbol{u}}\!\left( k \right) = {\left[ {\begin{array}{c}{{u_1}\!\left( k \right)} \cdots {{u_p}\!\left( k \right)}\end{array}} \right]^{\textrm{T}}}$
is the system input,
$p$
and
$q$
are the input and output dimensions of the system, respectively. We do not know the internal structure of the system, and can only obtain the online input data
${\boldsymbol{u}}$
and output data
${\boldsymbol{y}}$
of the system through sensor measurement. Equation (6) indicates that the output of the system at the next moment is not only related to the input of the current moment, but also related to the state of the current moment. However, we cannot measure the system state variable. Our idea is to carry out linear regression at each working point of the system according to the existing input and output data and then estimate some structural characteristics of the system. Therefore, we set up a linear regression model of input and output variation at each working point, as shown in Equation (7). Some works have shown that Equations (6) and (7) are equivalent under very broad conditions [Reference Hou and Bu10, Reference Bu, Hou, Yu and Fu11, Reference Hou and Zhu12].
\begin{equation}\Delta {y_n}\left({k + 1} \right)\,{=}\,\boldsymbol{\varPhi}_n^{\textrm{T}}\!\left( k \right)\Delta{\boldsymbol{u}}\!\left( k \right)\end{equation}
where
$\Delta {y_n}\left( {k + 1} \right) = {y_n}\left( {k + 1} \right) - {y_n}\!\left( k \right)$
,
$\Delta{\boldsymbol{u}}\!\left( k \right) = {\boldsymbol{u}}\!\left( k \right) - {\boldsymbol{u}}\left( {k - 1} \right)$
,
$1 \le n \le q$
,
${\boldsymbol{\varPhi} _n}\!\left( k \right) = {\left[ {\begin{array}{*{20}{c}}{{\phi _{n1}}} \quad \cdots \quad {{\phi _{np}}}\end{array}} \right]^{\textrm{T}}}$
represents the structural characteristics of the system corresponding to the nth output
${y_n}\left( {k + 1} \right)$
. Since it is obtained by linear regression using online input and output data at each operating point, we can call it pseudo partial derivative of the system. Observing Equation (7), we can see that the pseudo partial derivative
${\boldsymbol{\varPhi} _n}\!\left( k \right)$
depends on the system input variation
$\Delta{\boldsymbol{u}}\!\left( k \right)$
at the current time
$k$
and the nth system output variation
$\Delta {y_n}\left( {k + 1} \right)$
at the next time
$k + 1$
, both of which can be measured online, and
${\boldsymbol{u}}\!\left( k \right)$
needs to be designed by us. Next we need to construct a recursive format for the pseudo partial derivative of the system.
The estimation error is defined as
${e_n}\!\left( k \right) = \Delta {y_n}\left( {k + 1} \right) - \boldsymbol{\varPhi} _n^{\textrm{T}}\!\left( k \right)\Delta{\boldsymbol{u}}\!\left( k \right)$
, and the estimation loss is measured by Equation (8).
\begin{equation}f\left( {{\boldsymbol{\varPhi} _n}\!\left( k \right)} \right) = {\left\| {{e_n}\!\left( k \right)} \right\|^2} + \mu {\left\| {{\boldsymbol{\varPhi} _n}\!\left( k \right)} \right\|^2}\end{equation}
We need to use Newton method to minimise Equation (8), that is
\begin{equation}\mathop{\min }\limits_{\boldsymbol{\varPhi}} {\left\| {{e_n}\!\left( k \right)} \right\|^2} + \mu {\left\| {{\boldsymbol{\varPhi} _n}\!\left( k \right)} \right\|^2}\end{equation}
And then, we can obtain Equation (10).
\begin{equation}f\left( {{{\hat{\boldsymbol{\varPhi}} }_n}\left( {k + 1} \right)} \right)\,{=}\,f\left( {{{\hat{\boldsymbol{\varPhi}} }_n}\!\left( k \right) + \Delta } \right) = f\left( {{{\hat{\boldsymbol{\varPhi}} }_n}\!\left( k \right)} \right) + J_{f\left( \boldsymbol{\varPhi} \right)}^{\textrm{T}}\left( {{{\hat{\boldsymbol{\varPhi}} }_n}\!\left( k \right)} \right) \times \Delta + \frac{1}{2}{\Delta ^{\textrm{T}}}{H_{f\left( \boldsymbol{\varPhi} \right)}}\left( {{{\hat{\boldsymbol{\varPhi}} }_n}\!\left( k \right)} \right)\Delta \end{equation}
where
$\Delta = \varepsilon H_{f\left( \boldsymbol{\varPhi} \right)}^{ - 1}\left( {{{\hat{\boldsymbol{\varPhi}} }_n}\!\left( k \right)} \right){J_{f\left( \boldsymbol{\varPhi} \right)}}\left( {{{\hat{\boldsymbol{\varPhi}} }_n}\!\left( k \right)} \right)$
,
${J_{f\left( \boldsymbol{\varPhi} \right)}}$
is the Jacobi matrix about
$\boldsymbol{\varPhi} $
,
${H_{f\left( \boldsymbol{\varPhi} \right)}}$
is the Hessian matrix about
$\boldsymbol{\varPhi} $
. We can further obtain Equation (11).
\begin{equation}{\hat{\boldsymbol{\varPhi}} _n}\left( {k + 1} \right)\,{=}\,{\hat{\boldsymbol{\varPhi}} _n}\!\left( k \right) + \varepsilon H_{f\left( \boldsymbol{\varPhi} \right)}^{ - 1}\left( {{{\hat{\boldsymbol{\varPhi}} }_n}\!\left( k \right)} \right){J_{f\left( \boldsymbol{\varPhi} \right)}}\left( {{{\hat{\boldsymbol{\varPhi}} }_n}\!\left( k \right)} \right)\end{equation}
where
${J_{f\left( \boldsymbol{\varPhi} \right)}}\left( {{{\hat{\boldsymbol{\varPhi}} }_n}\!\left( k \right)} \right) = - 2\Delta{\boldsymbol{u}}\!\left( k \right){e_n}\!\left( k \right)$
,
${H_{f\left( \boldsymbol{\varPhi} \right)}}\left( {{{\hat{\boldsymbol{\varPhi}} }_n}\!\left( k \right)} \right) = 2\Delta{\boldsymbol{u}}\!\left( k \right)\Delta {\boldsymbol{u}^{\textrm{T}}}\!\left( k \right) + 2\mu \boldsymbol{I}$
,
$\varepsilon $
is the search step length of Newton method.
Observing Equation (8), the form is not the traditional mean square error function as shown in Equation (12). Instead, the term
$\mu {\left\| {{\boldsymbol{\varPhi} _n}\!\left( k \right)} \right\|^2}$
is added to Equation (12).
\begin{equation}f\left( {{\boldsymbol{\varPhi} _n}\!\left( k \right)} \right) = {\left\| {{e_n}\!\left( k \right)} \right\|^2}\end{equation}
This is because if Equation (12) is used as the loss function, the Hessian matrix in Equation (10) will change into
${H_{f\left( \boldsymbol{\varPhi} \right)}}\left( {{{\hat{\boldsymbol{\varPhi}} }_n}\!\left( k \right)} \right) = 2\Delta{\boldsymbol{u}}\!\left( k \right)\Delta {\boldsymbol{u}^{\textrm{T}}}\!\left( k \right)$
. At this point, every row and column of the Hessian matrix are linearly correlated, and its rank is 1. The inverse operation of the matrix will result in singularity. The loss function in the form of Equation (8) is actually a ridge regression technique, by adding regularisation term
$\mu {\left\| {{\boldsymbol{\varPhi} _n}\!\left( k \right)} \right\|^2}$
to the objective function, so that the inverse of Hessian matrix can be established.
By comparing the Hessian matrix forms obtained by the two loss functions, it can be seen that through ridge regression, the diagonal of the Hessian matrix obtained by Equation (12) is numerically loaded with a disturbance, which not only makes the Hessian matrix become full rank and be inversely solved, but also makes the numerical stability of the inverse operation good due to the existence of penalty of regularisation term. By substituting Hessian matrix obtained by Equation (8) into Equation (11), an adaptive filter can be obtained, as shown in Equation (13).
\begin{equation}{\hat{\boldsymbol{\varPhi}} _n}\left( {k + 1} \right) = {\hat{\boldsymbol{\varPhi}} _n}\!\left( k \right) + \eta {\left( {\Delta{\boldsymbol{u}}\!\left( k \right)\Delta {\boldsymbol{u}^{\textrm{T}}}\!\left( k \right) + \mu \boldsymbol{I}} \right)^{ - 1}}{e_n}\!\left( k \right)\Delta {\boldsymbol{u}^{\textrm{T}}}\!\left( k \right)\end{equation}
where
$\eta = - 4\varepsilon $
is the filter iteration step length. And then, we need to use the formula
\begin{equation}{\left( \boldsymbol{A} + \boldsymbol{BCD} \right)^{ - 1}} = {\boldsymbol{A}^{ - 1}} - {\boldsymbol{A}^{ - 1}}\boldsymbol{B}{\left( {\boldsymbol{D}{\boldsymbol{A}^{ - 1}}\boldsymbol{B} + {\boldsymbol{C}^{ - 1}}} \right)^{ - 1}}\boldsymbol{D}{\boldsymbol{A}^{ - 1}}\end{equation}
to deal with the inverse operation. Making
$\boldsymbol{A} = \mu \boldsymbol{I}$
,
$\boldsymbol{B} = \Delta{\boldsymbol{u}}\!\left( k \right)$
,
$\boldsymbol{C} = 1$
,
$\boldsymbol{D} = \Delta {\boldsymbol{u}^{\textrm{T}}}\!\left( k \right)$
, we can export more concise filtering device, as shown Equation (15).
\begin{equation}{\hat{\boldsymbol{\varPhi}} _n}\left( {k + 1} \right) = {\hat{\boldsymbol{\varPhi}} _n}\!\left( k \right) + \eta \frac{{\Delta {\boldsymbol{u}^{\textrm{T}}}\!\left( k \right)}}{{\mu + {{\left\| {\Delta{\boldsymbol{u}}\!\left( k \right)} \right\|}^2}}}{e_n}\!\left( k \right)\end{equation}
4.0 Design of data driven control law
Continuing the ridge regression thought in the previous section, the recursion of the control law is also the solution of one optimisation problem. We design the loss function of the optimisation problem to be the function of controlling signal variation, as shown in Equation (16). In order to achieve the purpose of tracking target signals, the loss function must be minimised.
\begin{equation}l\left( {\Delta{\boldsymbol{u}}\!\left( k \right)} \right) = {\left\| {y_n^ * \left( {k + 1} \right) - {y_n}\left( {k + 1} \right)} \right\|^2} + \lambda {\left\| {\Delta{\boldsymbol{u}}\!\left( k \right)} \right\|^2}\end{equation}
where
$y_n^ * $
is the tracking value of the nth output,
${\boldsymbol{y}^ * }\left( {k + 1} \right) = {\left[ {\begin{array}{*{20}{c}}{y_1^ * \left( {k + 1} \right)} \quad \cdots \quad {y_q^ * \left( {k + 1} \right)}\end{array}} \right]^{\textrm{T}}}$
. Defining control error corresponding to the nth output as
$e_n^ * \left( {k + 1} \right) = y_n^ * \left( {k + 1} \right) - {y_n}\left( {k + 1} \right)$
, we can obtain Equation (17) by substituting Equation (7) into it.
\begin{equation}e_n^ * \left( {k + 1} \right) = y_n^ * \left( {k + 1} \right) - {y_n}\!\left( k \right) - \hat{\boldsymbol{\varPhi}} _n^{\textrm{T}}\!\left( k \right)\Delta{\boldsymbol{u}}\!\left( k \right)\end{equation}
The loss function can be further written as
\begin{align}l\left( {\Delta{\boldsymbol{u}}\!\left( k \right)} \right) &= {\left\| {y_n^ * \left( {k + 1} \right) - {y_n}\!\left( k \right)} \right\|^2} - 2\left( {y_n^ * \left( {k + 1} \right) - {y_n}\!\left( k \right)} \right)\hat{\boldsymbol{\varPhi}} _n^{\textrm{T}}\!\left( k \right)\Delta{\boldsymbol{u}}\!\left( k \right)\nonumber\\ &\quad + \Delta {\boldsymbol{u}^{\textrm{T}}}\!\left( k \right){{\hat{\boldsymbol{\varPhi}} }_n}\!\left( k \right)\hat{\boldsymbol{\varPhi}} _n^{\textrm{T}}\!\left( k \right)\Delta{\boldsymbol{u}}\!\left( k \right) + \lambda {\left\| {\Delta{\boldsymbol{u}}\!\left( k \right)} \right\|^2}\end{align}
where
$\lambda {\left\| {\Delta{\boldsymbol{u}}\!\left( k \right)} \right\|^2}$
is still a regularisation term, to avoid the happening of singular and numerical instability. The gradient of the loss function to the control signal is the Jacobi matrix
\begin{equation}{J_{l\left( {\Delta{\boldsymbol{u}}} \right)}}\left( {\Delta{\boldsymbol{u}}\!\left( k \right)} \right) = - 2\left( {y_n^ * \left( {k + 1} \right) - {y_n}\!\left( k \right)} \right){\hat{\boldsymbol{\varPhi}} _n}\!\left( k \right) + 2{\hat{\boldsymbol{\varPhi}} _n}\!\left( k \right)\hat{\boldsymbol{\varPhi}} _n^{\textrm{T}}\!\left( k \right) + 2\lambda \Delta{\boldsymbol{u}}\!\left( k \right)\end{equation}
The optimal control signal gradient should be 0, and then we can obtain Equation (20).
\begin{equation}\Delta{\boldsymbol{u}}\!\left( k \right)\,{=}\,\left( {y_n^ * \left( {k + 1} \right) - {y_n}\!\left( k \right)} \right){\left( {{{\hat{\boldsymbol{\varPhi}} }_n}\!\left( k \right)\hat{\boldsymbol{\varPhi}} _n^{\textrm{T}}\!\left( k \right) + \lambda{\boldsymbol{I}}} \right)^{ - 1}}{\hat{\boldsymbol{\varPhi}} _n}\!\left( k \right)\end{equation}
Data-driven control system.
Binary wing model and actuator parameters
Otherwise,
${\boldsymbol{u}}\!\left( k \right) = {\boldsymbol{u}}\left( {k - 1} \right) + \Delta{\boldsymbol{u}}\!\left( k \right)$
. Using Formula (14) to deal with the inversion operation and making
$\boldsymbol{A} = \lambda{\boldsymbol{I}}$
,
$\boldsymbol{B} = {\hat{\boldsymbol{\varPhi}} _n}\!\left( k \right)$
,
$\boldsymbol{C} = 1$
,
$\boldsymbol{D} = \hat{\boldsymbol{\varPhi}} _n^{\textrm{T}}\!\left( k \right)$
, we can export more simple control signal recursive format as Equation (21).
\begin{equation}{\boldsymbol{u}}\!\left( k \right) = {\boldsymbol{u}}\left( {k - 1} \right) + \frac{{{{\hat{\boldsymbol{\varPhi}} }_n}\!\left( k \right)}}{{\lambda + \hat{\boldsymbol{\varPhi}} _n^{\textrm{T}}\!\left( k \right){{\hat{\boldsymbol{\varPhi}} }_n}\!\left( k \right)}}\left( {y_n^ * \left( {k + 1} \right) - {y_n}\!\left( k \right)} \right)\end{equation}
Optimal controller parameters
Root locus analysis of the uncontrolled system.
Data-driven controller parameters
System response (
$V = 20\,\,{\rm{m/s}}$
).
System state phase locus (
$V = 20\,\,{\rm{m/s}}$
).
Evaluation curve of system pseudo partial derivative (
$V = 20\,\,{\rm{m/s}}$
).
Observing Equation (21), we can see that the control signal
${\boldsymbol{u}}\!\left( k \right)$
at the current time
$k$
depends on the system pseudo partial derivative estimation
${\hat{\boldsymbol{\varPhi}} _n}\!\left( k \right)$
and system output
${y_n}\!\left( k \right)$
at the current time
$k$
, the target tracking value
$y_n^ * \left( {k + 1} \right)$
at the next time
$k + 1$
, and the control signal
${\boldsymbol{u}}\left( {k - 1} \right)$
at the previous time
$k - 1$
. Among them,
${y_n}\!\left( k \right)$
can be measured online,
$y_n^ * \left( {k + 1} \right)$
is known, and
${\hat{\boldsymbol{\varPhi}} _n}\!\left( k \right)$
can be estimated by the filter, which is a known quantity. It is necessary to point out that Equations (13) and (20) are the more essential forms, while Equations (15) and (21) are just more suitable for programming computation.
For convenience, the filter iteration step length corresponding to every output is equal, and the same to regularisation coefficient. Theoretically,
$q$
calculation results of recursive scheme should be exactly the same, but there will always be some disturbance in the measurement process. Therefore, the control signal can be averaged once, as shown in Equation (22).
\begin{equation}{\boldsymbol{u}}\!\left( k \right) = {\boldsymbol{u}}\left( {k - 1} \right) + \frac{1}{q}\sum\limits_{n = 1}^q {\frac{{{{\hat{\boldsymbol{\varPhi}} }_n}\!\left( k \right)}}{{\lambda + \hat{\boldsymbol{\varPhi}} _n^{\textrm{T}}\!\left( k \right){{\hat{\boldsymbol{\varPhi}} }_n}\!\left( k \right)}}} \left( {y_n^ * \left( {k + 1} \right) - {y_n}\!\left( k \right)} \right)\end{equation}
The model-free control scheme of MIMO discrete time system can be obtained by combining filter algorithm and control algorithm, as shown in Equation (23).
\begin{equation}\left\{ \begin{array}{l}{{\hat{\boldsymbol{\varPhi}} }_n}\!\left( k \right) = {{\hat{\boldsymbol{\varPhi}} }_n}\left( {k - 1} \right) + \eta \frac{{{e_n}\left( {k - 1} \right)\Delta {\boldsymbol{u}^{\textrm{T}}}\left( {k - 1} \right)}}{{\mu + {{\left\| {\Delta{\boldsymbol{u}}\left( {k - 1} \right)} \right\|}^2}}}\\{\boldsymbol{u}}\!\left( k \right) = {\boldsymbol{u}}\left( {k - 1} \right) + \frac{1}{q}\sum\limits_{n = 1}^q {\frac{{{{\hat{\boldsymbol{\varPhi}} }_n}\!\left( k \right)}}{{\lambda + \hat{\boldsymbol{\varPhi}} _n^{\textrm{T}}\!\left( k \right){{\hat{\boldsymbol{\varPhi}} }_n}\!\left( k \right)}}} \left( {y_n^ * \left( {k + 1} \right) - {y_n}\!\left( k \right)} \right)\end{array} \right.\end{equation}
where
$\Delta{\boldsymbol{u}}\left( { - 1} \right)\,{=}\,\boldsymbol{0}$
. Figure 2 shows the operation principle of the whole control system, in which
$\hat{\boldsymbol{\varPhi}} \!\left( k \right) = {\left[ {\begin{array}{*{20}{c}}{{{\hat{\boldsymbol{\varPhi}} }_1}\!\left( k \right)} \quad \cdots \quad {{{\hat{\boldsymbol{\varPhi}} }_q}\!\left( k \right)}\end{array}} \right]^{\textrm{T}}}$
.
5.0 Numerical simulation and discussion
For an aeroelastic system, when the flow velocity exceeds the critical flutter velocity, the system will have divergent motion, which leads to the failure of the structure. Therefore, in order to prevent divergent motion, stiffness of flutter system is often designed as hard characteristic in wind tunnel test. In this simulation example the limit spring on the pitching direction is used to limit the size of the pitch angle (
$ - {28^ \circ } \le \alpha \le {28^ \circ }$
), and the control surface is saturated with deflection (
$ - {30^ \circ } \le \beta \le {30^ \circ }$
). The micro-vibration of the wing near the equilibrium position is linear vibration, and large vibration is generated after the flutter. When the pitch angle exceeds the given threshold, the limiting spring will limit the pitching angle to increase further. Therefore, the wing exhibits limit cycle motion with finite amplitude when flutter occurs. The numerical simulation parameters are shown in Table 1.
The critical flutter velocity of the wing model can be determined by investigating the change of the eigenvalue of system matrix
$\boldsymbol{A}$
of the open-loop system in Equation (4) with the flow velocity, as shown in Fig. 3. Flutter critical velocity is
${V_f} = 17.5\,{\rm{m} \mathord{\left/ {\vphantom {m s}} \right.} \rm{s}}$
. The motion of the aeroelastic system under this wind speed is constant amplitude oscillation. At the critical flutter velocity, the motion of the system does not diverge due to the existence of the limiting spring, but behaves as a limit cycle.
The optimal controller is designed using LQR method in the research work of Darabseh et al. [Reference Darabseh, Tarabulsi and Mourad13]. The optimal control law is to seek a constant gain output feedback control law
\begin{equation}\delta = - {\boldsymbol{K}_{con}}\boldsymbol{Y} = - {\boldsymbol{K}_{con}}\boldsymbol{CX}\end{equation}
to minimise the following performance index.
\begin{equation}J = \frac{1}{2}\int_0^\infty {\left( {{\boldsymbol{X}^{\textrm{T}}}\boldsymbol{QX} + R{\delta ^2}} \right)dt} \end{equation}
where
${\boldsymbol{Q}}$
is the positive definite weighting matrix of the state variable, and
$R \gt 0$
is the weighting coefficient of the control variable. This objective function takes into account the requirements of both system response and control. According to the linear quadratic optimal control theory, if the performance index
$J$
is to be minimised, the control signal should be
\begin{equation}\delta = - {R^{ - 1}}{\boldsymbol{B}}\boldsymbol{PX}\end{equation}
where
${\boldsymbol{K}_{con}}\boldsymbol{C} = {R^{ - 1}}{\boldsymbol{B}}\boldsymbol{P}$
,
$\boldsymbol{P}$
is a symmetric matrix. The matrix satisfies the Riccati equation
\begin{equation}\boldsymbol{PA} + {\boldsymbol{A}^{\textrm{T}}}\boldsymbol{P} - \boldsymbol{PB}{R^{ - 1}}{\boldsymbol{B}^{\textrm{T}}}\boldsymbol{P} + \boldsymbol{Q} = \textbf{0}\end{equation}
where
${\boldsymbol{Q}}$
is the given weighted matrix,
$R$
is the given weighting coefficient.
For discrete system signal sampling, according to Nyquist-Shannon sampling theorem, sampling period is set as
${T_s} = 0.001\,{\rm{s}}$
. The related parameter settings of optimal controller and data-driven controller are shown in Tables 2 and 3, respectively.
Take different velocities to explore the ability of two controllers. Under the speed of
$V = 20\,{\rm{m} \mathord{\left/ {\vphantom {m s}} \right.} \rm{s}}$
, control was implemented at
$3.5\,{\rm{s}}$
respectively. The time history of the system motion before and after the control of the two controllers is shown in Fig. 4.
It can be seen from Fig. 4 that both controllers have good flutter suppression effect when the flow velocity just exceeds the critical velocity, and the response can be attenuated to 0 by only a few deflections of the control surface, but the deflection number of the optimal controller is more than that of the data-driven controller. The saturation of the control surface has a certain influence on the effect of the optimal controller, and the control surface reaches the upper limit of saturation in the initial stage. This is because the optimal controller has no adaptive ability, and the control strategy tends to be ‘drastic’ in the initial stage, resulting in the control surface easily reaching the upper limit of saturation. In contrast, data-driven controllers are much more ‘restrained’. The phase locus of the system before and after control is shown in Fig. 5. The phase locus shows that although both controllers can make the system state return to the equilibrium point quickly, the regression path of the optimal controller is too sharp, while that of the data-driven controller is relatively smooth. Figure 6 shows the evaluation curve of the filter on the system pseudo partial derivative. It can be seen that the estimation value can quickly converge to the real value after the controller starts to work.
Under the speed of
$V = 24\,{\rm{m} \mathord{\left/ {\vphantom {m s}} \right.} \rm{s}}$
, control was also implemented at
$3.5\,\rm{s}$
respectively. The time history of the system motion before and after the control of the two controllers is shown in Fig. 7.
System response (
$V = 24\,{\rm{m/s}}$
).
System state phase locus (
$V = 24\,{\rm{m/s}}$
).
Evaluation curve of system pseudo partial derivative (
$V = 24\,{\rm{m/s}}$
).
As shown in Fig. 7, the flutter of the wing becomes more intense with the increase of flow velocity. The flutter suppression effect of the optimal controller decreases significantly, and it takes a long time to attenuate the response, while the data-driven controller is obviously better. Due to its powerful adaptive ability, the data-driven controller can adjust the parameters in time according to the size of the system output, the amplitude of the control signal can change rapidly, the filter can master the system characteristics throughout the process, and the control cost can be accurately adjusted to avoid excessive control inertia. On the contrary, the optimal controller has no adaptive ability, the control signal has large inertia and the amplitude adjustment is slow, resulting in that the additional aerodynamic force generated by every deflection of the control surface is either too large or too small and the optimal controller needs to repeatedly adjust for a long time, which is the defect of traditional model-based optimal control. We can regard the reciprocal of the deflection number of the control surface as the capability margin of the controller in a certain sense. With the increase of the deflection number, the capability margin of the controller keeps approaching 0, which means that the controller is getting closer and closer to failure. Obviously, the data-driven controller still has a large margin space.
The design of optimal controller depends on the modeling of aeroelastic system, which is complicated for designers, while the design of data-driven controller avoids the modeling process. So from the practical point of view, the data-driven controller is much cheaper, which has a systematic design process, and it only needs to set a few initial parameters. In addition, for flutter, designers always want the response to decay quickly, otherwise prolonged oscillations increase the risk of fatigue in the wing structure. So the data-driven controller is also better in terms of the time it takes for the response to decay.
The phase locus of the system before and after control is shown in Fig. 8. At this time, the phase locus shows that when the wind speed increases, the equilibrium point regression paths of the two controllers become longer, while the optimal controller obviously increases more, and the sharp problem of its regression path still exists. Under the action of data-driven controller, the system state will find a relatively short path for equilibrium point regression. However, under the action of the optimal controller, the system state will be attracted by the limit cycle and can only return to the equilibrium point after continuous circling around the periphery. Figure 9 shows the evaluation curve of the filter on the system pseudo partial derivative under the speed of
$V = 24\,{\rm{m} \mathord{\left/ {\vphantom {m s}} \right.} \rm{s}}$
. It can be seen that the estimation value can still quickly converge to the real value after the controller starts to work, and the structural characteristics of the system at this time have changed due to the change of the flow velocity, so the pseudo partial derivative estimated by the filter is also different from that under the speed of
$V = 20\,{\rm{m} \mathord{\left/ {\vphantom {m s}} \right.} \rm{s}}$
.
6.0 Conclusion
In this paper, an adaptive control method based on ridge regression is proposed and applied to active flutter suppression, and the following conclusions are obtained:
-
(1) The controller designed by this method only uses the real-time measured online input and output data of the closed-loop controlled system, which can realise the system correction without identifying the system structure. Compared with the traditional feedback control, which requires model-based design, the design based on ridge regression is model-free.
-
(2) The controller structure designed by this method is relatively simple, easy to modify and program. The calculation amount is acceptable. It is very suitable for the digital control system widely used in the aviation industry.
-
(3) The focus of this method is to establish a regression model at each working point of the system. In essence, this is a kind of dynamic linearisation. Different from the traditional linearisation in a small range near the equilibrium point, this kind of dynamic linearisation is accurate.
-
(4) The cores of this method are inner product operation and regularisation. Inner product operation establishes projective information between data, and regularisation eliminates singularity through diagonal loading. It is worth pointing out that the method proposed in this paper is the foundation, and appropriate operations for these two cores can be used to extend the method.
Competing of interests
The authors have no competing interests to declare that are relevant to the content of this article.
Appendix
For subsonic incompressible flows, according to Theodorson theory, the aerodynamic lift force acting on a binary wing with span of
${s_p}$
and the aerodynamic torque against the elastic axis can be written as
\begin{align} L & \displaystyle = \pi {\rho _a}{b^2}{s_p}\left[ {\ddot h + V\dot \alpha - b\bar a\ddot \alpha - \left( {{V \mathord{\left/ {\vphantom {V \pi }} \right.} \pi }} \right){T_4}\dot \beta - \left( {{b \mathord{\left/ {\vphantom {b \pi }} \right.} \pi }} \right){T_1}\ddot \beta } \right]\nonumber\\ &\quad+ 2\pi {\rho _a}Vb{s_p}C\!\left( k \right)\left[ {V\alpha + \dot h + b\left( {0.5 - \bar a} \right)\dot \alpha + \left( {{V \mathord{\left/ {\vphantom {V \pi }} \right.} \pi }} \right){T_{10}}\beta + \left( {{b \mathord{\left/ {\vphantom {b {2\pi }}} \right.} {2\pi }}} \right){T_{11}}\dot \beta } \right]\\[-12pt] \nonumber\end{align}
\begin{align} {T_\alpha } &\displaystyle = \pi {\rho _a}{b^2}{s_p}\left\{ \begin{array}{l}\displaystyle b\bar a\ddot h - Vb\left( {0.5 - \bar a} \right)\dot \alpha - {b^2}\left( {{1 \mathord{\left/ {\vphantom {1 8}} \right.} 8} + {{\bar a}^2}} \right)\ddot \alpha - \left( {{{{V^2}} \mathord{\left/ {\vphantom {{{V^2}} \pi }} \right.} \pi }} \right)\left( {{T_4} + {T_{10}}} \right)\beta \\ \displaystyle + \left( {{{Vb} \mathord{\left/ {\vphantom {{Vb} \pi }} \right.} \pi }} \right)\left[ { - {T_1} + {T_8} + \left( {\bar c - \bar a} \right){T_4} - 0.5{T_{11}}} \right]\dot \beta + \left( {{{{b^2}} \mathord{\left/ {\vphantom {{{b^2}} \pi }} \right.} \pi }} \right)\left[ {{T_7} + \left( {\bar c - \bar a} \right){T_1}} \right]\ddot \beta \end{array} \right\}\nonumber\\ &\quad + 2\pi {\rho _a}V{b^2}{s_p}\left( {\bar a + 0.5} \right)C\!\left( k \right)\left\{ {V\alpha + \dot h + b\left( {0.5 - \bar a} \right)\dot \alpha + \left( {{V \mathord{\left/ {\vphantom {V \pi }} \right.} \pi }} \right){T_{10}}\beta + \left( {{b \mathord{\left/ {\vphantom {b {2\pi }}} \right.} {2\pi }}} \right){T_{11}}\dot \beta } \right\}\end{align}
where
${T_i}\left( {i = 1,4,7,8,10,11} \right)$
is the Theodorson constant dependent on the position of the elastic axis and the position of the hinge line of the control surface, the specific algorithm is as follows; Theodorson function
$C\!\left( k \right)$
is a complex function.
\begin{equation}{T_1} = - \frac{{2 + {{\bar c}^2}}}{3}\sqrt {1 - {{\bar c}^2}} + \bar c\arccos \bar c\end{equation}
\begin{equation}{T_4} = \bar c\sqrt {1 - {{\bar c}^2}} - \arccos \bar c\end{equation}
\begin{equation}{T_7} = \bar c\frac{{7 + 2{{\bar c}^2}}}{8}\sqrt {1 - {{\bar c}^2}} - \left( {\frac{1}{8} + {{\bar c}^2}} \right)\arccos \bar c\end{equation}
\begin{equation}{T_8} = - \frac{1}{3}\left( {1 + 2{{\bar c}^2}} \right)\sqrt {1 - {{\bar c}^2}} + \bar c\arccos \bar c\end{equation}
\begin{equation}{T_{10}} = \sqrt {1 - {{\bar c}^2}} + \arccos \bar c\end{equation}
\begin{equation}{T_{11}} = \left( {2 - \bar c} \right)\sqrt {1 - {{\bar c}^2}} + \left( {1 - 2\bar c} \right)\arccos \bar c\end{equation}
Jones form of Theodorson function is conducive to the time-domain state-space modeling of aeroelastic systems, which can be written as
\begin{equation}C\left( s \right) = 1 - \frac{{0.165s}}{{s + 0.0455\frac{V}{b}}} - \frac{{0.335s}}{{s + 0.3\frac{V}{b}}} = \frac{1}{2} + \frac{{0.0075\frac{V}{b}}}{{s + 0.0455\frac{V}{b}}} + \frac{{0.10055\frac{V}{b}}}{{s + 0.3\frac{V}{b}}}\,{=}\,\frac{1}{2} + \frac{{{z_1}\frac{V}{b}}}{{s + {p_1}\frac{V}{b}}} + \frac{{{z_2}\frac{V}{b}}}{{s + {p_2}\frac{V}{b}}}\end{equation}
where
${z_1} = 0.0075$
,
${z_2} = 0.10055$
,
${p_1} = 0.0455$
,
${p_2} = 0.3$
, and
$s$
is Laplace variable. If
$C\left( s \right)$
is regarded as the transfer function of a certain system, the corresponding state-space equation of the system can be written as
\begin{equation}\left\{ \begin{array}{l}\left[ {\begin{array}{*{20}{c}}{{{\dot x}_{a1}}}\\{{{\dot x}_{a2}}}\end{array}} \right] = \frac{V}{b}\left[ {\begin{array}{*{20}{c}}{ - {p_1}}\quad 0\\ 0\quad { - {p_2}}\end{array}} \right]\left[ {\begin{array}{*{20}{c}}{{x_{a1}}}\\{{x_{a2}}}\end{array}} \right] + \left[ {\begin{array}{*{20}{c}}1\\1\end{array}} \right]r = V{{\boldsymbol{A}}_a}{{\boldsymbol{x}}_a} + V{{\boldsymbol{K}}_a}{\boldsymbol{q}}_{\boldsymbol{s}} + {{\boldsymbol{D}}_a}{{\dot {\boldsymbol{q}}}_s}\\[8pt] {y_a} = \frac{V}{b}\left[ {\begin{array}{*{20}{c}}{{z_1}}\quad {{z_2}}\end{array}} \right]\left[ {\begin{array}{*{20}{c}}{{x_{a1}}}\\{{x_{a2}}}\end{array}} \right] + \frac{1}{2}r\end{array} \right.\end{equation}
\begin{equation}r = V\alpha + \dot h + b\left( {0.5 - \bar a} \right)\dot \alpha + \left( {{V \mathord{\left/ {\vphantom {V \pi }} \right.} \pi }} \right){T_{10}}\beta + \left( {{b \mathord{\left/ {\vphantom {b {2\pi }}} \right.} {2\pi }}} \right){T_{11}}\dot \beta = {{\boldsymbol{S}}_1}\boldsymbol{\dot{\boldsymbol{q}}}_{\boldsymbol{s}} + V{{\boldsymbol{S}}_2}{\boldsymbol{q}}_{\boldsymbol{s}}\end{equation}
where
$r$
is the system input,
${{\boldsymbol{S}}_1} = \left[ {\begin{array}{*{20}{c}}1 \quad{b\left( {0.5 - \bar a} \right)}\quad {\left( {{b \mathord{\left/ {\vphantom {b {2\pi }}} \right.} {2\pi }}} \right){T_{11}}}\end{array}} \right]$
,
${{\boldsymbol{S}}_2} = \left[ {\begin{array}{*{20}{c}}0 \quad 1 \quad{{{{T_{10}}} \mathord{\left/ {\vphantom {{{T_{10}}} \pi }} \right.} \pi }}\end{array}} \right]$
;
${x_{a1}}$
and
${x_{a2}}$
are the aerodynamic state variables,
${{\boldsymbol{x}}_a} = {\left[ {\begin{array}{*{20}{c}}{{x_{a1}}} \quad {{x_{a2}}}\end{array}} \right]^{\textrm{T}}}$
;
${{\boldsymbol{A}}_a} = \frac{1}{b}\left[ {\begin{array}{*{20}{c}}{ - {p_1}}\quad 0\\0\quad { - {p_2}}\end{array}} \right]$
,
${{\boldsymbol{K}}_a} = \left[ {\begin{array}{*{20}{c}}1\\1\end{array}} \right]{{\boldsymbol{S}}_2}$
,
${{\boldsymbol{D}}_a} = \left[ {\begin{array}{*{20}{c}}1\\1\end{array}} \right]{{\boldsymbol{S}}_1}$
. The calculation formula of aerodynamic force can be rewritten as
\begin{equation} - L = - \pi {\rho _a}{b^2}{s_p}\left[ {V\dot \alpha - \left( {{V \mathord{\left/ {\vphantom {V \pi }} \right.} \pi }} \right){T_4}\dot \beta - \left( {{b \mathord{\left/ {\vphantom {b \pi }} \right.} \pi }} \right){T_1}\ddot \beta } \right] - \pi {\rho _a}{b^2}{s_p}\left[ {\ddot h - b\bar a\ddot \alpha } \right] - 2\pi {\rho _a}Vb{s_p}{y_a}\end{equation}
\begin{align}{T_\alpha } &= \pi {\rho _a}{b^2}{s_p}\left\{ \begin{array}{l}\displaystyle - Vb\left( {0.5 - \bar a} \right)\dot \alpha - \left( {{{{V^2}} \mathord{\left/ {\vphantom {{{V^2}} \pi }} \right.} \pi }} \right)\left( {{T_4} + {T_{10}}} \right)\beta \\[4pt] \displaystyle+ \left( {{{Vb} \mathord{\left/ {\vphantom {{Vb} \pi }} \right.} \pi }} \right)\left[ { - {T_1} + {T_8} + \left( {\bar c - \bar a} \right){T_4} - 0.5{T_{11}}} \right]\dot \beta + \left( {{{{b^2}} \mathord{\left/ {\vphantom {{{b^2}} \pi }} \right.} \pi }} \right)\left[ {{T_7} + \left( {\bar c - \bar a} \right){T_1}} \right]\ddot \beta \end{array} \right\}\nonumber\\ &\quad+\pi {\rho _a}V{b^3}{s_p}\left[ {\bar a\ddot h - b\left( {{1 \mathord{\left/ {\vphantom {1 8}} \right.} 8} + {{\bar a}^2}} \right)\ddot \alpha } \right] + 2\pi {\rho _a}V{b^2}{s_p}\left( {\bar a + 0.5} \right){y_a}\end{align}
The aerodynamic force vector
${\boldsymbol{F}}_{\boldsymbol{ae}}$
is written as the sum of the non-cyclic part and the cyclic part, i.e
\begin{equation}{\boldsymbol{F}}_{\boldsymbol{ae}} = {\boldsymbol{F}_{nc}} + {\boldsymbol{F}_c}\end{equation}
The non-cyclic part
${\boldsymbol{F}_{nc}}$
can be written as
\begin{equation}{\boldsymbol{F}_{nc}} = - {\boldsymbol{M}_{nc}}\boldsymbol{\ddot{\boldsymbol{q}}}_{\boldsymbol{s}} - V{\boldsymbol{D}_{nc}}\boldsymbol{\dot{\boldsymbol{q}}}_{\boldsymbol{s}} - {V^2}{\boldsymbol{K}_{nc}}{\boldsymbol{q}}_{\boldsymbol{s}}\end{equation}
where
$${\boldsymbol{M}_{nc}} = {\rho _a}{b^2}{s_p}\left[ {\begin{array}{c@{\quad}c@{\quad}c}\displaystyle\pi &{ - \pi b\bar a} &{ - b{T_1}}\\[5pt] \displaystyle{ - \pi b\bar a} &{\pi {b^2}\left( {{1 \mathord{\left/ {\vphantom {1 8}} \right.} 8} + {{\bar a}^2}} \right)} &{ - {b^2}\left( {{T_7} + \left( {\bar c - \bar a} \right)T_1} \right)}\\[5pt] 0 \displaystyle&0 &0\end{array}}\right] $$
$${\boldsymbol{D}_{nc}} = {\rho _a}{b^2}{s_p}\left[ {\begin{array}{c@{\quad}c@{\quad}c}0 &\pi &{ - {T_4}}\\[4pt] 0 &{\pi b\left( {0.5 - \bar a} \right)} &{ - b\left( { - {T_1} + {T_8} + \left( {\bar c - \bar a} \right){T_4} - 0.5{T_{11}}} \right)}\\[4pt] 0 &0 &0\end{array}} \right]$$
$${\boldsymbol{K}_{nc}} = {\rho _a}{b^2}{s_p}\left[ {\begin{array}{c@{\quad}c@{\quad}c}0 &0 &0\\0 &0 &{{T_4} + {T_{10}}}\\0 &0 &0\end{array}} \right]$$
The cyclic part
${\boldsymbol{F}_c}$
can be written as
\begin{equation}{\boldsymbol{F}_c} = - {\rho _a}Vb{s_p}{\boldsymbol{R}_c}\left( {\frac{V}{b}\left[ {\begin{array}{*{20}{c}}{{z_1}}\,\,{{z_2}}\end{array}} \right]\left[ {\begin{array}{*{20}{c}}{{x_{a1}}}\\{{x_{a2}}}\end{array}} \right] + \frac{1}{2}\left( {\boldsymbol{S}_1{{\dot{\boldsymbol{q}}}_s} + V{\boldsymbol{S}_2}{\boldsymbol{q}}_{\boldsymbol{s}}} \right)} \right)\end{equation}
where
${\boldsymbol{R}_c} = {\left[ {\begin{array}{*{20}{c}}{2\pi } \,\,{ - 2\pi b\left( {0.5 + \bar a} \right)}\,\, 0\end{array}} \right]^{\textrm{T}}}$
. The cyclic part
${\boldsymbol{F}_c}$
can be further written as
\begin{equation}{\boldsymbol{F}_c} = - V{\boldsymbol{D}_c}\boldsymbol{\dot{\boldsymbol{q}}}_{\boldsymbol{s}} - {V^2}{\boldsymbol{K}_c}{\boldsymbol{q}}_{\boldsymbol{s}} - {V^2}{\boldsymbol{E}_c}{\boldsymbol{x}_a}\end{equation}
where
${\boldsymbol{D}_c} = \frac{1}{2}{\rho _a}b{s_p}{\boldsymbol{R}_c}{\boldsymbol{S}_1}$
,
${\boldsymbol{K}_c} = \frac{1}{2}{\rho _a}b{s_p}{\boldsymbol{R}_c}{\boldsymbol{S}_2}$
,
${\boldsymbol{E}_c} = {\rho _a}{s_p}{\boldsymbol{R}_c}\left[ {\begin{array}{*{20}{c}}{{z_1}}\,\, {{z_2}}\end{array}} \right]$
. In summary, the original wing motion equation can be written as
\begin{equation}\left( {{\boldsymbol{M}}_{\boldsymbol{s}} + {\boldsymbol{M}_{nc}}} \right)\boldsymbol{\ddot{\boldsymbol{q}}}_{\boldsymbol{s}} + \left( {{\boldsymbol{D}}_{\boldsymbol{s}} + V{\boldsymbol{D}_{nc}} + V{\boldsymbol{D}_c}} \right)\boldsymbol{\dot{\boldsymbol{q}}}_{\boldsymbol{s}} + \left( {{\boldsymbol{K}}_{\boldsymbol{s}} + {V^2}{\boldsymbol{K}_{nc}} + {V^2}{\boldsymbol{K}_c}} \right){\boldsymbol{q}}_{\boldsymbol{s}} = - {V^2}{\boldsymbol{E}_c}{\boldsymbol{x}_a} + {\boldsymbol{G}}_{\boldsymbol{s}}\delta \end{equation}
The state space expression of the aeroelastic equation of a binary wing can be obtained by combining the wing motion equation with the aerodynamic state equation.
\begin{equation}{\dot{\boldsymbol{X}}} = {\boldsymbol{AX}} + {\boldsymbol{B}}\delta \end{equation}
where
$\boldsymbol{A} = \left[ {\begin{array}{c@{\quad}c@{\quad}c}{{\boldsymbol{0}_{3 \times 3}}} &{{\boldsymbol{I}_{3 \times 3}}} &{{\boldsymbol{0}_{3 \times 2}}}\\[4pt] { - {\boldsymbol{M}^{ - 1}}\boldsymbol{K}} &{ - {\boldsymbol{M}^{ - 1}}\boldsymbol{D}} &{ - {V^2}{\boldsymbol{M}^{ - 1}}{\boldsymbol{E}_c}}\\[4pt] {V{\boldsymbol{K}_\boldsymbol{a}}} &{\boldsymbol{D}_\boldsymbol{a}} &{V{\boldsymbol{A}_\boldsymbol{a}}}\end{array}} \right]$
,
$\boldsymbol{B} = \left[ {\begin{array}{*{20}{c}}{{\boldsymbol{0}_{3 \times 1}}}\\{{\boldsymbol{M}^{ - 1}}{\boldsymbol{G}}_{\boldsymbol{s}}}\\{{\boldsymbol{0}_{2 \times 1}}}\end{array}} \right]$
,
$\boldsymbol{X} = {\left[ {\begin{array}{*{20}{c}} \displaystyle h \,\,\alpha \,\, \beta\,\, {\dot h}\,\, {\dot \alpha }\,\, {\dot \beta }\,\, {{x_{a1}}}\,\, {{x_{a2}}}\end{array}} \right]^{\text{T}}}$
,
${\boldsymbol{M}} = {\boldsymbol{M}}_{\boldsymbol{s}} + {{\boldsymbol{M}}_{nc}}$
,
$\boldsymbol{D} = {\boldsymbol{D}}_{\boldsymbol{s}} + V{\boldsymbol{D}_{nc}} + V{\boldsymbol{D}_c}$
,
$\boldsymbol{K} = {\boldsymbol{K}}_{\boldsymbol{s}} + {V^2}{\boldsymbol{K}_{nc}} + {V^2}{\boldsymbol{K}_c}$
.
