Abstract
The work presented in this paper is concerned with the identification of a piecewise linear aeroelastic system from inputoutput data. The main challenge with this problem is that the data are available only as a mixture of observations generated by a finite set of different interacting linear subsystems such that one does not know a prior which subsystem has generated which data, that is, the switching points of the freeplay nonlinearity. The linear part of the nonlinear aeroelastic system is represented by the orthonormal basis functions constructed by the physical poles of the linear part, and the nonlinear part is represented by a Hammerstein model. By a simple rearrangement of the data corresponding to the degreeoffreedom of freeplay and selecting a segment of the data, the identification of the physical poles could be reduced to a linear parametric problem. Afterwards, estimates of the unknown parameters of linear regression models are calculated by processing respective particles of inputoutput data. The iterative sequence of the switching points is constructed, and solved by a method synthesizing the noniterative and iterative algorithms. Then the parameters of the linear and nonlinear parts of the nonlinear system including the switching points are successfully obtained. A twodimensional airfoil with nonlinear structural freeplay in the pitch degreeoffreedom is presented to demonstrate the validity of the proposed identification algorithm.
1. Introduction
Aeroelasticity is one of the prime concerns for an aircraft or a flight vehicle. The interaction of the aerodynamics and the structural dynamics of an aircraft could result in several different aeroelastic behaviors. These phenomena often exhibit nonlinear behaviors such as limit cycle oscillations (LCO), chaotic motions, coexisting stable solutions, bifurcations and so on [14]. Nonlinear factors in aeroelastic systems is ubiquity such as the material nonlinearity, the geometric nonlinearity, the freeplay nonlinearity and so on in structure, and the shockboundary layer interaction, the flow separation, the unstable eddy current and so on in aerodynamics [5]. Freeplay nonlinearity can arise from worn hinges and loose attachments which are generally related with aircraft aging. The combination of freeplay nonlinearities with the aerodynamic and structural nonlinearities can vary the system’s response and change its instability from a supercritical to a subcritical one. These changes can result in catastrophic damages. To that end, many researchers have investigated numerically and experimentally the effects of freeplay nonlinearity on aeroelastic systems [610].
In order to study a nonlinear aeroelastic system with a freeplay nonlinearity, its accurate mathematics model needs to be constructed first. In engineering practice, freeplay such as gap in the bolt connection [9] may never be correctly measured, so the mathematic models of these systems have to be constructed by identification of experimental data. For example, the aeroelastic systems with freeplay can be identified by representing the freeplay as a polynomial [11] and a hyperbolic tangent function [12], respectively. It was pointed out that the locations of the switching points of the freeplay nonlinearity are of importance in the modeling process, which not only influence the location of bifurcation points, but also impact the dynamic behaviors [10]. Since the switching points in these identified models are eliminated artificially, their dynamic behaviors are different from the original systems in essence [11, 12]. Recently, a nonparameter method was developed to identify the switching points of freeplay [13]. But this method needs the bandwidth of the kernel function, and unfortunately the selection of the bandwidth depends strongly on the engineer experience and influences the precision of the switching points.
The objective of this work is to identify and construct the mathematic model of the nonlinear aeroelastic system with a freeplay nonlinearity. This is achieved and demonstrated by performing system identification of a twodimensional airfoil with a freeplay nonlinearity in the pitch degreeoffreedom. The system model makes use of pitch–plunge governing equations with quasisteady approximation of the aerodynamic loads.
2. The aeroelastic model
The model here is a twodimensional airfoil with nonlinear structural freeplay in the pitch degreeoffreedom. A schematic of the airfoil section with a trailingedge flap for control actuation [14] is shown in Fig. 1. This system, which has been studied extensively in the literature, can be modeled by a pair of coupled secondorder differential equations [15]:
where $h$ denotes the plunge displacement of the airfoil and $\alpha $ represents the pitch angle. Other variables include the nondimensional distance between the elastic axis and the center of mass ${x}_{\alpha}$, the mass of the wing $m$, the mass moment of inertia ${I}_{\alpha}$, semichord length $b$, structural damping coefficients in plunge and pitch ${c}_{h}$ and ${c}_{\alpha}$, and spring constants ${k}_{h}$ and ${k}_{\alpha}$. The lift ${L}_{ae}$ and moment ${M}_{ae}$ are determined by quasisteady aerodynamic theory. The parameters ${c}_{l\alpha}$ and ${c}_{l\alpha}$ are introduced to represent the lift and moment coefficients for angle of attack, ${c}_{l\beta}$ and ${c}_{m\beta}$ are introduced to represent the lift and moment coefficients for control surface position, $a$ is the nondimensional distance between the middle of the chord and the elastic axis, $\rho $ is density of air and $V$ is free stream velocity.
Fig. 1Airfoil model
a)
b)
With a structural freeplay gap, the pitchrestoring momentrotation relationships ${M}_{\alpha}\left(\alpha \right)$ are illustrated in Fig. 2 as:
where ${\delta}_{1}$ and ${\delta}_{2}$ are the switching points, ${M}_{0}$ is the restoring moment at the start of the freeplay and shown in Fig. 2.
Fig. 2Freeplay nonlinearity in the pitch direction model
Fig. 3Saturation function ω
By substituting Eq. (2) and Eq. (3) into Eq. (1), the transformed equations of motions in the state space form become:
where $\mathit{X}={\left[\begin{array}{cc}\begin{array}{cc}h& \alpha \end{array}& \begin{array}{cc}\dot{h}& \dot{\alpha}\end{array}\end{array}\right]}^{T}$, $T$ denotes the matrix or vector transpose,
$\omega =f\left(\alpha \right)={k}_{\alpha}\alpha {M}_{\alpha}\left(\alpha \right)$, $\mathit{M}=\left[\begin{array}{cc}m& m{x}_{\alpha}b\\ m{x}_{\alpha}b& {I}_{\alpha}\end{array}\right]$, $\mathit{K}=\left[\begin{array}{cc}{k}_{h}& \rho {V}^{2}b{c}_{l\alpha}\\ 0& {k}_{\alpha}\rho {V}^{2}{b}^{2}{c}_{m\alpha}\end{array}\right]$,
$\mathit{C}=\left[\begin{array}{cc}{c}_{h}+\rho Vb{c}_{l\alpha}& \rho V{b}^{2}{c}_{l\alpha}\left(0.5a\right)\\ \rho V{b}^{2}{c}_{m\alpha}& {c}_{\alpha}\rho V{b}^{3}{c}_{m\alpha}\left(0.5a\right)\end{array}\right]$, ${\mathit{F}}_{1}=\left[\begin{array}{c}\rho {V}^{2}b{c}_{l\beta}\\ \rho {V}^{2}{b}^{2}{c}_{m\beta}\end{array}\right]$, ${\mathit{F}}_{2}=\left[\begin{array}{c}0\\ 1\end{array}\right]$,
$\mathit{A}=\left[\begin{array}{cc}0& \mathit{I}\\ {\mathit{M}}^{1}\mathit{K}& {\mathit{M}}^{1}\mathit{C}\end{array}\right]$, ${\mathit{B}}_{1}=\left[\begin{array}{c}0\\ {\mathit{M}}^{1}{\mathit{F}}_{1}\end{array}\right]$, ${\mathit{B}}_{2}=\left[\begin{array}{c}0\\ {\mathit{M}}^{1}{\mathit{F}}_{2}\end{array}\right]$, $\mathit{I}$ denotes the identity matrix.
The saturation function $\omega $ is illustrated in Fig. 3 as:
Substituting Eq. (6) into Eq. (5) as:
where ${\mathit{D}}_{1}={\mathit{B}}_{2}{k}_{\alpha}{\delta}_{2}$, ${\mathit{D}}_{2}={\mathit{B}}_{2}{k}_{\alpha}$, ${\mathit{D}}_{3}={\mathit{B}}_{2}{k}_{\alpha}{\delta}_{1}$, $\mathit{H}={\mathit{B}}_{2}{M}_{0}$.
The airfoil model output is the pitch angle $\alpha $. Using a nonlinear feedback linear fractional transformation (LFT) [16], Eq. (5) represents as:
where $\mathit{P}$ is the nominal plant and:
where ${P}_{1k}\left(k=1,2\right)$ are the functions related to the input ${\left\{\beta ,\omega \right\}}^{T}$and output $\alpha $ signals. These transfer function matrices are built from the $\mathit{M}$, $\mathit{C}$, $\mathit{K}$, ${\mathit{F}}_{1}$ and ${\mathit{F}}_{2}$ matrices of the nominal aeroelastic model. ${F}_{l}\left[\bullet ,\bullet \right]$ denotes the low LFT. Eq. (8) is illustrated in Fig. 4 [16]. According to Eq. (5), ${P}_{11}$ and ${P}_{12}$ can be represented as:
where $\mathit{Q}=\left[\begin{array}{cc}\begin{array}{cc}0& 1\end{array}& \begin{array}{cc}0& 0\end{array}\end{array}\right]$ and $s$ is Laplace variable.
Fig. 4LFT with saturation function fα
According to LFT algebra, Eq. (8) can also be represented as:
where ${P}_{11}\beta $ is the linear part of the output $\alpha $, and ${P}_{12}f\left(\alpha \right)$ is the nonlinear part of the output $\alpha $ and a Hammerstein model, which comprises of a static nonlinearity followed by a linear time invariant system [16].
3. Identification of the poles of the linear subsystem
3.1. Model discretization
According to Ref. [17], ${P}_{11}$ and ${P}_{12}$ of Eq. (10) can be discretized as:
where ${t}_{s}$ is the sampling time and $q$ is the forward timeshift operator: $q{\alpha}_{k}={\alpha}_{k+1}$. ${\mathit{Q}}_{d}=\mathit{Q}$, ${\mathit{A}}_{d}={\mathrm{e}}^{\mathit{A}{t}_{s}}$, ${\mathit{B}}_{1d}={\mathit{A}}^{1}\left({\mathrm{e}}^{\mathit{A}{t}_{s}}\mathit{I}\right){\mathit{B}}_{1}$, ${\mathit{B}}_{2d}={\mathit{A}}^{1}\left({\mathrm{e}}^{\mathit{A}{t}_{s}}\mathit{I}\right){\mathit{B}}_{2}$.
${P}_{21}$ and ${P}_{22}$ of Eq. (12) can be represented in the rational polynomial form as:
where $A\left(q\right)=1+{a}_{1}{q}^{1}+\cdots +{a}_{n}{q}^{n}$, $B\left(q\right)={b}_{1}{q}^{1}+\cdots +{b}_{n}{q}^{n}$,
$C\left(q\right)={c}_{1}{q}^{1}+\cdots +{c}_{n}{q}^{n}$. $n$ is the order number of the characteristic polynomial of the matrix ${\mathit{A}}_{\mathit{d}}$.
According to Eq. (13), Eq. (11) can be discretized as (in Ref. [17]):
where ${C}_{1}\left(q\right)={c}_{1}+{c}_{2}{q}^{1}+\cdots +{c}_{n}{q}^{n+1}$, ${\alpha}_{k}$ and ${\beta}_{k}$ are the values of $\alpha \left(t\right)$ and $\beta \left(t\right)$ at time $k$, respectively.
According to Eq. (6) and Eq. (7), Eq. (14c) can be represented as:
where ${C}_{2}\left(q\right)={k}_{\alpha}{C}_{1}\left(q\right)$,${r}_{1}=\left({k}_{\alpha}{\delta}_{2}{M}_{0}\right)\left({c}_{1}+{c}_{2}+\cdots +{c}_{n}\right)$,
${r}_{3}=\left({k}_{\alpha}{\delta}_{1}{M}_{0}\right)\left({c}_{1}+{c}_{2}+\cdots +{c}_{n}\right)$. The parameters ${r}_{1}$ and ${r}_{3}$ are real constant numbers.
Compare the first or third linear subsystems, that is, Eq. (15a) or Eq. (15c), with the linear part ${P}_{21}$ of this nonlinear system, only one constant number is increased. But it doesn’t matter their poles which are the same as ones of the linear part ${P}_{21}$. However, the second linear subsystem, that is Eq. (15b), can be represented as ${A}_{0}\left(q\right){\alpha}_{k}=B\left(q\right){\beta}_{k}$ (${A}_{0}\left(q\right)=A\left(q\right){k}_{\alpha}C\left(q\right)$), whose poles are different to ones of the linear part ${P}_{21}$.
3.2. Identification of the poles
The airfoil model should be fully motivated. This means that the simulation processes are fully through three linear subsystems, that is, Eq. (15ac). Let us rearrange the outputs ${\left\{{\alpha}_{k}\right\}}_{k=1}^{N}$ of the system in an ascending order according to their values [18]. Defining $\stackrel{~}{\mathrm{N}}=\left\{1,2,\cdots ,N\right\}$ and a onetoone map $g:\stackrel{~}{N}\to \stackrel{~}{N}$ on it, that is, ${\stackrel{~}{\alpha}}_{i}={\alpha}_{g\left(i\right)}$ ($i=\mathrm{1,2},\cdots ,N$). According to Eq. (15ac), the rearranged outputs${\left\{{\stackrel{~}{\alpha}}_{i}\right\}}_{i=1}^{N}$ can be partitioned into three data sets: righthand side data set (${N}_{1}$ samples) with values higher than or equal to ${\delta}_{2}$, that is, ${\stackrel{~}{\alpha}}_{i}\ge {\delta}_{2}$; middle data set (${N}_{2}$ samples) with values higher than ${\delta}_{1}$ but lower than ${\delta}_{2}$, that is, ${\delta}_{1}<{\stackrel{~}{\alpha}}_{i}<{\delta}_{2}$; lefthand side data set (${N}_{3}$ samples) with values lower than or equal to ${\delta}_{1}$, that is, ${\stackrel{~}{\alpha}}_{i}\le {\delta}_{1}$. Here $N={N}_{1}+{N}_{2}+{N}_{3}$. Although the switching points of the freeplay nonlinearity are unknown, a constant number ${h}_{1}$ which is not only a little smaller than ${\stackrel{~}{\alpha}}_{N}$ but also such that ${h}_{1}>{\delta}_{2}$ can be selected. When ${\stackrel{~}{\alpha}}_{i}>{h}_{1}$, that is, ${\alpha}_{g\left(i\right)}>{h}_{1}$, the inputoutput relationship of the system in time $g\left(i\right)+1$ is a linear subsystem represented by Eq. (15a). Similarly, a constant number ${h}_{3}$ which is not only a little larger than ${\stackrel{~}{\alpha}}_{1}$ but also such that ${h}_{3}<{\delta}_{1}$ can also be selected. When ${\stackrel{~}{\alpha}}_{i}<{h}_{3}$, that is, ${\alpha}_{g\left(i\right)}<{h}_{3}$, the inputoutput relationship of the system in time$g\left(i\right)+1$is a linear subsystem represented by Eq. (15c). Selecting ${\stackrel{~}{\alpha}}_{i}$, $\forall i\in \left[{N}_{2}+{N}_{3}+{l}_{1},N\right]$ (${l}_{1}$ is a positive integer, and ${l}_{1}<{N}_{1}$) in the rearranged outputs ${\left\{{\stackrel{~}{\alpha}}_{i}\right\}}_{i=1}^{N}$, then the inputoutput relationship of the system in time $g\left(i\right)+1$is a linear subsystem represented by Eq. (15a). If the outputs ${\left\{{\alpha}_{k}\right\}}_{k=1}^{N}$ without measurement noise, the equation of this linear subsystem solved by least square estimation is:
where $\mathit{W}={\left[\begin{array}{cc}\begin{array}{cc}{\alpha}_{g\left({N}_{2}+{N}_{3}+{l}_{1}\right)+1}& {\alpha}_{g\left({N}_{2}+{N}_{3}+{l}_{1}+1\right)+1}\end{array}& \begin{array}{cc}\cdots & {\alpha}_{g\left(N\right)+1}\end{array}\end{array}\right]}^{T}$, $\widehat{\mathit{\mu}}$ is the estimation of the parameter:
${\mathit{U}}_{1}=\left[\begin{array}{cccc}\begin{array}{c}{\beta}_{g\left({N}_{2}+{N}_{3}+{l}_{1}\right)1}\\ {\beta}_{g\left({N}_{2}+{N}_{3}+{l}_{1}+1\right)1}\end{array}& \begin{array}{c}{\beta}_{g\left({N}_{2}+{N}_{3}+{l}_{1}\right)2}\\ {\beta}_{g\left({N}_{2}+{N}_{3}+{l}_{1}+1\right)2}\end{array}& \begin{array}{c}\cdots \\ \cdots \end{array}& \begin{array}{c}{\beta}_{g\left({N}_{2}+{N}_{3}+{l}_{1}\right)n}\\ {\beta}_{g\left({N}_{2}+{N}_{3}+{l}_{1}+1\right)n}\end{array}\\ \begin{array}{c}\vdots \\ {\beta}_{g\left(N\right)1}\end{array}& \begin{array}{c}\vdots \\ {\beta}_{g\left(N\right)2}\end{array}& \begin{array}{c}\ddots \\ \cdots \end{array}& \begin{array}{c}\vdots \\ {\beta}_{g\left(N\right)n}\end{array}\end{array}\right]$
${\mathit{U}}_{2}=\left[\begin{array}{ccccc}\begin{array}{c}{\alpha}_{g\left({N}_{2}+{N}_{3}+{l}_{1}\right)}\\ {\alpha}_{g\left({N}_{2}+{N}_{3}+{l}_{1}+1\right)}\end{array}& \begin{array}{c}{\alpha}_{g\left({N}_{2}+{N}_{3}+{l}_{1}\right)1}\\ {\alpha}_{g\left({N}_{2}+{N}_{3}+{l}_{1}+1\right)1}\end{array}& \begin{array}{c}\cdots \\ \cdots \end{array}& \begin{array}{c}{\alpha}_{g\left({N}_{2}+{N}_{3}+{l}_{1}\right)n+1}\\ {\alpha}_{g\left({N}_{2}+{N}_{3}+{l}_{1}+1\right)n+1}\end{array}& \begin{array}{c}1\\ 1\end{array}\\ \begin{array}{c}\vdots \\ {\alpha}_{g\left(N\right)}\end{array}& \begin{array}{c}\vdots \\ {\alpha}_{g\left(N\right)1}\end{array}& \begin{array}{c}\ddots \\ \cdots \end{array}& \begin{array}{c}\vdots \\ {\alpha}_{g\left(N\right)n+1}\end{array}& \begin{array}{c}\vdots \\ 1\end{array}\end{array}\right]$
$\mathit{U}=\left[\begin{array}{cc}{\mathit{U}}_{1}& {\mathit{U}}_{2}\end{array}\right].$
If the outputs ${\left\{{\alpha}_{k}\right\}}_{k=1}^{N}$ with measurement noise, the equation of this linear subsystem can be solved by the biaseliminated least square method [19]. The parameters of the first linear subsystem can be estimated, and its poles can also be solved. Similarly, selecting ${\stackrel{~}{\alpha}}_{i}\text{,}$$\forall i\in \left[1,{N}_{3}{l}_{2}\right]$ (${l}_{2}$ is a positive integer, and ${l}_{2}<{N}_{3}$) in the rearranged outputs ${\left\{{\stackrel{~}{\alpha}}_{i}\right\}}_{i=1}^{N}$, then the inputoutput relationship of the system in time $g\left(i\right)+1$ is a linear subsystem represented by Eq. (15c). The parameters of this subsystem can be solved like Eq. (16) and Eq. (17) or by the biaseliminated least square method [19] and are the same of ones of the first linear subsystem except the parameter ${r}_{1}\ne {r}_{3}$.
4. Identification of the switching points
Let us define ${\lambda}_{k}={P}_{22}f\left({\alpha}_{k}\right)$, which is the nonlinear part of the output ${\alpha}_{k}$ in Eq. (11) and also a discrete Hammerstein model [16], which comprises of a saturation function $f\left(\bullet \right)$ followed by a linear time invariant system ${P}_{22}\left(q\right)$ connected, shown in Fig. 5.
Define a switching function $h\left(u\right)$as follows:
The saturation function$f\left(\bullet \right)$is:
where ${g}_{1}\left({\alpha}_{k1},{\delta}_{1},{\delta}_{2}\right)={\alpha}_{k1}\left[h\left({\delta}_{2}{\alpha}_{k1}\right)h\left({\delta}_{1}{\alpha}_{k1}\right)\right]$;
${g}_{2}\left({\alpha}_{k1},{\delta}_{1},{\delta}_{2}\right)=h\left({\alpha}_{k1}{\delta}_{2}\right)$; ${g}_{3}\left({\alpha}_{k1},{\delta}_{1},{\delta}_{2}\right)=h\left({\delta}_{1}{\alpha}_{k1}\right)$; ${g}_{4}\left({\alpha}_{k1},{\delta}_{1},{\delta}_{2}\right)=1$;
${d}_{1}={k}_{\alpha}$; ${d}_{2}={k}_{\alpha}{\delta}_{2}$; ${d}_{3}={k}_{\alpha}{\delta}_{1}$; ${d}_{4}={M}_{0}$; the parameters ${d}_{1}$, ${d}_{2}$, ${d}_{3}$ and ${d}_{4}$ are unknown.
The linear system ${P}_{22}\left(q\right)$ in Eq. (14ac) can be represented by using the orthonormal basis functions ${\left\{{B}_{l}\left(q\right)\right\}}_{l=0}^{p1}$ [20] which are constructed by using the poles of linear part of the nonlinear system as follows:
where ${P}_{22}\left(q\right)\in \mathit{R}\left(\mathbb{Q}\right)$, $\mathbb{Q}$ is the unit circle [16] and the parameter ${e}_{l}\in R$ is unknown.
The inputoutput relationship of the system illustrated by Fig. 5 is:
where ${\lambda}_{k}$ and ${\eta}_{k}$ represent the system output and measurement noise at time $k$, respectively.
Fig. 5. Hammerstein model
The linear system ${P}_{21}\left(q\right)$ in Eq. (14)ac) represents by using the orthonormal basis functions ${\left\{{B}_{l}\left(q\right)\right\}}_{l=0}^{p1}$ [20] as follows:
When Eq. (21) and Eq. (22) are substituted into Eq. (14ac), the inputoutput relationship can be written as:
Let’s now define:
${B}_{0}\left(q\right){g}_{3}\left({\alpha}_{k1},{\delta}_{1},{\delta}_{2}\right),{B}_{0}\left(q\right){g}_{4}\left({\alpha}_{k1},{\delta}_{1},{\delta}_{2}\right),\cdots ,{B}_{p1}\left(q\right){g}_{1}\left({\alpha}_{k1},{\delta}_{1},{\delta}_{2}\right),$
${{B}_{p1}\left(q\right){g}_{2}\left({\alpha}_{k1},{\delta}_{1},{\delta}_{2}\right),\left.{B}_{p1}\left(q\right){g}_{3}\left({\alpha}_{k1},{\delta}_{1},{\delta}_{2}\right),{B}_{p1}\left(q\right){g}_{4}\left({\alpha}_{k1},{\delta}_{1},{\delta}_{2}\right)\right]}^{T}.$
When Eq. (24a) and Eq. (24b) are placed into Eq. (23), the later results in the regression vector:
Now, with the simulated data set ${\left\{{\beta}_{k},{\alpha}_{k}\right\}}_{k=1}^{N}$ of the nonlinear aeroelastic system with freeplay and defining the matrices ${\mathit{X}}_{N}\triangleq {\left[{\alpha}_{1},\cdots ,{\alpha}_{N}\right]}^{T}\text{,}$${\mathit{\Gamma}}_{N}\triangleq {\left[{\eta}_{1},\cdots ,{\eta}_{N}\right]}^{T}$ and ${\mathit{\Phi}}_{N}\triangleq \left[{\mathit{\varphi}}_{1},\cdots ,{\mathit{\varphi}}_{N}\right]$, we obtain:
First, the initial values of the switching points ${\delta}_{1}$ and ${\delta}_{2}$ are given and ${\mathit{\varphi}}_{k}$ is calculated. The estimated $\widehat{\mathit{\theta}}$ of the parameters $\mathit{\theta}$ are obtained by using the least square estimation or the biaseliminated least square method [19]. According to the definition of the parameter $\mathit{\theta}$, its first $p$ parameters are the coefficients of the linear part ${P}_{21}\left(q\right)$, and the else are the coefficients of the nonlinear part ${P}_{22}\left(q\right)f\left(\bullet \right)$. Since the parameters of the saturation function$f\left(\bullet \right)$ and the linear system ${P}_{22}\left(q\right)$ are coupling, the unique estimated parameters ${d}_{1}\text{,}$${d}_{2}\text{,}$${d}_{3}\text{,}$${d}_{4}$and ${b}_{l}$$\left(l=0,1,\cdots ,p1\right)$ can be obtained if the parameter $\mathit{d}={\left[\begin{array}{cc}\begin{array}{cc}{d}_{1}& {d}_{2}\end{array}& \begin{array}{cc}{d}_{3}& {d}_{4}\end{array}\end{array}\right]}^{T}$ are normalized, that is, ${\Vert \mathit{d}\Vert}_{2}=1$ [16]. Since the parameters $\mathit{d}$ are normalized, the identified parameters ${d}_{1}$, ${d}_{2}$, ${d}_{3}$ and ${d}_{4}$ lose the physical significance defining in the above, but the switching points ${\delta}_{1}$ and ${\delta}_{2}$ can be estimated by ${\widehat{\delta}}_{1}={\widehat{d}}_{3}/{\widehat{d}}_{1}$ and ${\widehat{\delta}}_{2}={\widehat{d}}_{2}/{\widehat{d}}_{1}$. Using the renewed switching points ${\widehat{\delta}}_{1}$ and ${\widehat{\delta}}_{2}$ repeat the above identified procedure until the convergence of the switching points ${\delta}_{1}$ and ${\delta}_{2}$. Then all parameters of this nonlinear model including the switching points can be obtained.
5. A simulated example
The parameters of a twodimensional airfoil with nonlinear structural freeplay in the pitch degreeoffreedom to be used in the numerical simulation are given in Table 1. The output of the simulated measured system is the pitch angle $\alpha \left(t\right)$, which is corrupted with additive Gaussian distributed white noise with a signaltonoise ratio of 20 dB, and the input is the flap deflection $\beta \left(t\right)$, which is a zeromean Gaussian distributed white noise with standard deviation 10. A key issue in the time marching integration of a piecewise linear system is accurately integrating to the “switching points” where the change in linear subdomains occurs. It is indicated that a standard time marching scheme, for example, the RungeKutta method with the uniform time step, may lead to inconsistent or even incorrect numerical results [21]. In this paper, the response of the system in the first 50 seconds is computed by using the Henon’s method [22, 23] whose time step is ${t}_{s}=\text{0.001 s}$.
It is easy to know the number of the degreeoffreedom for the discrete system. This sumulated model has two degreeoffreedoms. Assuming $A\left(q\right)$ and $B\left(q\right)$ are:
Table 1Parameters of a nonlinear airfoil model with freeplay
Parameters  Value  Parameters  Value 
$V$  6 m/s  $a$  0.6 
$b$  0.135 m  $\rho $  1.225 kg/m^{3} 
$m$  12.387 kg  ${x}_{\alpha}$  0.2466 
${I}_{\alpha}$  0.065 m^{2}kg  ${k}_{h}$  2844.4 N/m 
${c}_{\alpha}$  0.180 m^{2}kg/s  ${c}_{h}$  27.43 kg/s 
${c}_{l\alpha}$  6.28  ${c}_{m\alpha}$  0.628 
${c}_{l\beta}$  3.358  ${c}_{m\beta}$  0.635 
${k}_{\alpha}$  2.82 Nm/rad  ${\delta}_{1}$  0.05 
${\delta}_{2}$  0.25  ${M}_{0}$  0.282 Nm 
According to section 3.2, we may as well assume ${h}_{1}=$ 0.4 and ${h}_{2}=$ 0.1 since ${\stackrel{~}{\alpha}}_{N}=$ 0.7307 and ${\stackrel{~}{\alpha}}_{1}=$ 0.7454. Here selecting ${h}_{1}=$0.4 and ${h}_{2}=$ 0.1 make full of the inputoutput data in order to eliminate measurement noise. When ${h}_{1}=$ 0.4 is used, the parameter $\mathit{\mu}$ is estimated by the biaseliminated least square method [19], and given in Table 2; when ${h}_{2}=$ 0.1 is used, the parameter$\mathit{}\mathit{\mu}$ is estimated by the biaseliminated least square method [19], and also given in Table 2. The estimated parameters are the same in ${h}_{1}=$ 0.4 and ${h}_{2}=$ 0.1 cases except ${r}_{1}\ne {r}_{3}$. The parameter ${r}_{1}\ne {r}_{3}$ is caused by the asymmetric switching points of a freeplay nonlinearity.
Table 2Estimated parameters
${h}_{1}=0.4$  ${h}_{2}=0.1$  
${b}_{0}$  1.5059×10^{6}  1.5059×10^{6} 
${b}_{1}$  1.5265×10^{6}  1.5265×10^{6} 
${b}_{2}$  1.4923×10^{6}  1.4923×10^{6} 
${b}_{3}$  1.5107×10^{6}  1.5107×10^{6} 
${a}_{1}$  3.9931  3.9931 
${a}_{2}$  5.9798  5.9798 
${a}_{3}$  3.9801  3.9801 
${a}_{4}$  0.9935  0.9935 
${r}_{1}$  1.8882×10^{9}   
${r}_{3}$    6.2941×10^{10} 
Table 3Comparison between the identified modes and the true modes of the system
True modes  ${h}_{1}=0.4$  ${h}_{2}=0.1$  
${f}_{1}$, Hz  1.1660  1.1660  1.1660 
${f}_{2}$, Hz  2.6509  2.6509  2.6509 
${\zeta}_{1}$  0.2081  0.2081  0.2081 
${\zeta}_{2}$  0.1049  0.1049  0.1049 
According to the relationship of poles of the discrete system with the continue system $\gamma ={\mathrm{e}}^{{t}_{s}\kappa}$ ($\kappa $ and $\gamma $ denote a pole of the discrete system and the one of the corresponding continue system pole, respectively), the frequency ${f}_{i}$ and damping ${\zeta}_{i}$ ($i=1,2$) of these two physical modes of the linear part of the nonlinear airfoil are given in Table 3. In Table 3, the second column indicates the frequency and damping of the true modes. The third and fourth columns present the estimated frequency and damping of the true modes using the proposed method in section 3. Comparing the estimated modes and true modes in Table 3, it is obvious to see that the frequency ${f}_{i}$ and damping ${\zeta}_{i}$$\left(i=1,2\right)$ of these two physical modes are estimated consistently.
The poles of the linear part of the nonlinear model are the same with the ones of the first and third linear subsystem, that is Eq. (15a) and Eq. (15c). The poles and the orthonormal basis functions constructed by them [20] are:
The convergence of the switching points ${\delta}_{2}$ and ${\delta}_{1}$ are illustrated in Fig. 6, when the initial switching points ${\delta}_{2}$ and ${\delta}_{1}$ are (0.18, 0.04), (0.24, 0.02), (0.30, 0.03), (0.36, 0.04) and (0.40, 0.10). It is found from Fig. 6 that the switching points are converged to the true values only through the 6 steps iteration, and meanwhile it is stated that the convergence of the switching points is good when the different initial value of the switching points ${\delta}_{2}$ and ${\delta}_{1}$ are utilized. When the initial value of the switching point ${\delta}_{2}$ and ${\delta}_{1}$ are (0.40, 0.10), through the 20 steps iteration the estimated parameters of the system are: ${\widehat{\tau}}_{0}=\text{7.4606\xd7}{\text{10}}^{\text{3}}\text{,}$${\widehat{\tau}}_{1}=\text{4.9613\xd7}{\text{10}}^{\text{3}}\text{,}$${\widehat{\tau}}_{2}=\text{1.7787\xd7}{\text{10}}^{\text{2}}\text{,}$${\widehat{\tau}}_{3}=\text{5.2651\xd7}{\text{10}}^{\text{3}}$, ${\widehat{d}}_{1}=\text{0.9642}$, ${\widehat{d}}_{2}=\text{0.2424}$, ${\widehat{d}}_{3}=\text{0.04806}$, ${\widehat{d}}_{4}=\text{0.09645}$, ${\widehat{e}}_{0}=\text{1.6616\xd7}{\text{10}}^{\text{2}}$, ${\widehat{e}}_{1}=\text{1.5177\xd7}{\text{10}}^{\text{2}}\text{,}$${\widehat{e}}_{2}=\text{7.4796\xd7}{\text{10}}^{\text{2}}\text{,}$${\widehat{e}}_{3}=\text{1.6150\xd7}{\text{10}}^{\text{2}}$. The identified switching points are ${\widehat{\delta}}_{2}=\text{0.2514}$ and ${\widehat{\delta}}_{1}=\text{0.04985}$, and the true switching points are ${\delta}_{2}=\text{0.25}$ and ${\delta}_{1}=\text{0.05}$. The relative errors which are caused by the number computation are 0.5573 % and 0.3052 % and satisfied the engineering requirements.
Fig. 6Convergence of switching points δ2 and δ1: a) the initial values of switching point δ2 are 0.18, 0.24, 0.30, 0.36 and 0.4; the true value of switching point δ2 is 0.25, b) the initial values of switching point δ1 are 0.04, 0.02, 0.03, 0.04 and 0.10; the true value of switching point δ1 is 0.05
a)
b)
Fig. 7Output of nonlinear system (dotted line) and identified nonlinear model (solid line) with Gauss white noise input: a) with measurement noise, b) without measurement noise
a)
b)
Fig. 7 shows the output of the nonlinear system (dotted line) and the identified nonlinear model (solid line) with the Gauss white noise input. Because of the output of nonlinear system with the measurement noise in Fig. 7(a), there is a clear difference between them. In Fig. 7(b) they are almost coincided due to the outputs without the measurement noise. The validity of the identified nonlinear model is checked by several other input signals without shown here.
6. Conclusions
In this paper, a novel identification algorithm is proposed for the estimation of the nonlinear system with freeplay. In this algorithm, an orthonormal basis function expansion and a Hammerstein model are applied to represent the linear and nonlinear parts of the system, respectively; and an orthonormal basis function expansion and the constructional basis functions ${g}_{i}\left({\alpha}_{k1},{\delta}_{1},{\delta}_{2}\right)\left(i=1,2,\cdots ,4\right)$ are used to represent the linear and nonlinear functions of the Hammerstein model. An approach synthesized of the noniterative and iterative algorithms is implemented to estimate the parameters of the system including the switching points. Furthermore, a method is proposed to estimate physical poles of the system which could be essentially reduced to a linear parametric problem by a simple data rearrangement in an ascending order according to their values. Thus, the data which is concerned with the degreeoffreedom of the freeplay nonlinearity are partitioned into three data sets that correspond to distinct linear regression models. Later on the estimates of unknown parameters of linear regression models can be calculated by processing respective sets of the rearranged input and associated output observations. The consistent estimation of these physical poles plays a crucial role for the generation of the orthonormal basis functions and for the model approximation results.
References

Dowell E. H., Tang D. Nonlinear aeroelasticity and unsteady aerodynamics. AIAA Journal, Vol. 40, Issue 9, 2002, p. 16971707.

Gilliat H. C., Strganac T. W., Kurdila A. J. An investigation of internal resonance in aeroelastic systems. Nonlinear Dynamics, Vol. 31, Issue 1, 2003, p. 122.

Raghothama A., Narayanan S. Nonlinear dynamics of a twodimensional airfoil by incremental harmonic balance method. Journal of Sound and Vibration, Vol. 226, Issue 3, 1999, p. 493517.

Liu L., Wong Y. S., Lee B. H. K. Application of the centre manifold theory in nonlinear aeroelasticity. Journal of Sound and Vibration, Vol. 234, Issue 4, 2000, p. 641659.

Dowell E. H., Edwards J. W., Strganac T. Nonlinear aeroelasticity. Journal of Aircraft, Vol. 40, Issue 5, 2003, p. 857874.

Tang D., Dowell E. H. Aeroelastic airfoil with freeplay at angle of attack with gust excitation. AIAA Journal, Vol. 48, Issue 2, 2010, p. 427442.

Gordon J. T., Meyer E. E., Minogue R. L. Nonlinear stability analysis of control surface flutter with freeplay effects. Journal of Aircraft, Vol. 45, Issue 6, 2008, p. 19041916.

Dimitriadis G. Bifurcation analysis of aircraft with structural nonlinearity and freeplay using numerical continuation. Journal of Aircraft, Vol. 45, Issue 3, 2008, p. 893905.

Tang D., Attar P., Dowell E. H. Flutter/limit cycle oscillation analysis and experiment for wingstore model. AIAA Journal, Vol. 44, Issue 7, 2006, p. 16621675.

Liu L., Wong Y. S. Nonlinear aeroelastic analysis using the point transformation method, Part 1: Freeplay model. Journal of Sound and Vibration, Vol. 253, Issue 2, 2002, p. 447469.

Abdelkefi A., Vasconcellos R., Marques F. D., Hajj M. et al. Modeling and identification of freeplay nonlinearity. Journal of Sound and Vibration, Vol. 331, Issue 8, 2012, p. 18981907.

Kukreja S. L., Brenner M. Nonlinear blackbox modeling of aeroelastic systems using structure detection: application to F/A18 data. Journal of Guidance, Control, and Dynamics, Vol. 30, Issue 2, 2007, p. 557564.

Popescu C. A., Wong Y. S., Lee B. H. K. An expert system for predicting nonlinear aeroelastic behavior of an airfoil. Journal of Sound and Vibration, Vol. 319, Issue 5, 2009, p. 13121329.

O’Neil T., Strganac T. W. Aeroelastic response of a rigid wing supported by nonlinear springs. Journal of Aircraft, Vol. 35, Issue 4, 1998, p. 616622.

Prazenica R. J., Reisenthel P. H., Kurdila A. J., Brenner M. J. Volterra kernel extrapolation for modeling nonlinear aeroelastic systems at novel flight conditions. Journal of Aircraft, Vol. 44, Issue 1, 2007, p. 149162.

Baldelli D. H., Lind R., Brenner M. Nonlinear aeroelastic/aeroservoelastic modeling by blockoriented identification. Journal of Guidance, Control, and Dynamics, Vol. 28, Issue 5, 2005, p. 10561064.

Toth R., Heuberger P. S. C., Vandenhof P. M. J. Discretization of linear parametervarying statespace representations. IET Control Theory Appl., Vol. 4, Issue 10, 2010, p. 20822096.

Pupeikis R. On the identification of Hammerstein systems having saturationlike functions with positive slopes. Informatica, Vol. 17, Issue 1, 2006, p. 5568.

Zheng W. X. Transfer function estimation from noisy input and output data. International Journal of Adaptive Control and Signal Processing, Vol. 12, Issue 4, 1998, p. 365380.

Ninness B., Gustafsson F. A unifying construction of orthonormal bases for system identification. IEEE Transactions on Automatic Control, Vol. 42, Issue 4, 1997, p. 515521.

Liu L., Dowell E. H. Hamonic balance approach for an airfoil with a freeplay control surface. AIAA Journal, Vol. 43, Issue 4, 2005, p. 802815.

Henon M. On the numerical computation of Poincare maps. Physical D, Vol. 5, Issue 23, 1982, p. 412414.

Conner M. D., Virgin L. N., Dowell E. H. Accurate numerical integration of state space models for aeroelastic systems with freeplay. AIAA Journal, Vol. 34, Issue 11, 1996, p. 22022205.
About this article
This research is supported by the National Natural Science Foundation of China (Grant No. 10872089 and 11102085) and the Specialized Research Fund for the Doctoral Program of Higher Education of China (Grant No. 201132180002).