Forecasting bifurcations is a significant challenge and an important need in several applications. Most of the existing forecasting approaches focus on bifurcations in nonoscillating systems. However, subcritical and supercritical flutter (Hopf) bifurcations are very common in a variety of systems, especially fluid–structural systems. This paper presents a unique approach to forecast (nonlinear) flutter based on observations of the system only in the prebifurcation regime. The proposed method is based on exploiting the phenomenon of critical slowing down (CSD) in oscillating systems near certain bifurcations. Techniques are introduced to enhance the prediction accuracy for cases of low-frequency oscillations and large-dimensional dynamical systems. The method is applied to an aeroelastic system responding to gust loads. Numerical results are provided to demonstrate the performance of the method in predicting the postbifurcation regime accurately in both supercritical and subcritical cases.

## Introduction

Bifurcations occur in the dynamics of complex systems, from ecosystems to engineered systems. Such phenomena lead to various types of stability issues and can cause dramatic changes in the dynamics. Therefore, forecasting such bifurcations, i.e., predicting bifurcations and the bifurcation diagram with measurements only from the prebifurcation regime is a significant challenge and an important need. This is especially important for complex large-dimensional systems when an accurate model of the system is not easily available or when the system properties/parameters are unknown.

For complex nonlinear dynamical systems, it is usually not possible to have a model of the system accurate enough to predict bifurcations of the system without performing measurements in the postbifurcation regime. That is because existing model-based methods require extensive model calibration and system identification. However, even without knowing the system equations, tipping points in the dynamics can be forecasted using early warning indicators which can be extracted from system behavior in the prebifurcation regime [1].

Several forecasting methods using observing the recovery rate of the system from perturbations have been proposed [26] (see Ref. [7] for a review). However, these methods still have limitations and cannot be applied to more complex and large dimensional systems.

A common type of the bifurcations is Hopf. These bifurcations have been observed in a variety of systems such as aeroelastic systems [810], machine tools [11], automotive components [12], delayed predator–prey systems [13], nonlinear circuits [14], vibro-impact systems [15], etc. Although analyses of Hopf bifurcations using model-based methods have been done, forecasting of this type of bifurcations is much less studied [7]. Nonetheless, developing a forecasting method for this class of bifurcations is necessary especially because of the wide spread occurrence of such bifurcations in engineered systems.

There are several challenges in forecasting Hopf bifurcations. First, the system recovers to its equilibrium state with oscillations, and one cannot use all of the measured data for forecasting since adjacent points are at different phases. In this case, one can simply choose the local maxima of the oscillating recovery data, approximate the recovery rate, and use the method suggested by Lim and Epureanu [5,6]. However, that approach requires a large number of local peaks in the response during recovery and it is hard or impossible to use when there are only a few peaks in the measured recovery. Second, Hopf bifurcations may occur in large dimensional systems with several simultaneously active modes. Although many modes may be active in the recovery of a system from a perturbation, only one of these modes causes the bifurcation [16]. Accurate forecasting of the postbifurcation regime cannot be achieved without taking this into account. Moreover, most existing methods such as stochastic approaches for forecasting critical transitions focus exclusively on small perturbations and hence are able to forecast only low amplitudes. A method applicable for large perturbations is needed to forecast larger ranges of the bifurcation diagram [17].

In this paper, we present a novel approach to forecast Hopf bifurcations in large dimensional systems based on measurements of time series data in prebifurcation regime. The approach is mainly based on the phenomenon of CSD which accompanies such bifurcations. We use the technique presented by Lim and Epureanu [5,6] as the basis of our method, and extend it for bifurcation forecasting of more complex oscillating systems.

Numerical simulations are provided to demonstrate the accuracy of the proposed method. The method is applied to an aeroelastic system composed of an airfoil in pitch and plunge. Instead of using experimental data, we used surrogate date obtained as time series from simulation of a model of the nonlinear aeroelastic system. This model has two mechanical degrees-of-freedom (pitch and plunge). The aerodynamic part of the model leads to a more complex eight-dimensional nonlinear aeroelastic system.

To reflect more realistic situations where the system perturbations occur during operation, gust loads are applied as disturbances to the system while the system recoveries are measured. These perturbations are more accurate representations of realistic cases than simply choosing perturbed initial conditions in simulations. Both supercritical and subcritical bifurcations of the aeroelastic system are studied by changing the flow speed and structural nonlinearities in the system.

The results show that the method accurately predicts the bifurcation point and also the postbifurcation regime in both supercritical and subcritical cases despite the fact that it uses only prebifurcation regime data and it does not use a model of the system.

## Forecasting Method

Consider a codimension one nonlinear system with the parameter $μ$ and amplitude of $r$. The change rate of system amplitude can be generally written as
$r˙=f(μ,r)$
(1)
Using a Taylor series around the bifurcation point ($μ=μc$), the recovery rate can be expressed as
$r˙=r (p(r)+α1(r) (μ−μc)+α2(r) (μ−μc)2+HOT)$
(2)

where $p(r)$, $α1(r)$, and $α2(r)$ are polynomial functions independent of parameter $μ$. This relation means neither that the dynamics of the system have been linearized in state space nor that the dynamics have small amplitudes. The linearization is only in the parameter space about the critical value $μc$ of the parameter. Thus, we are not restricted to small perturbations.

For parameter values close to the bifurcation, in this paper we assume that the system has second-order polynomial dependence on the parameter, and higher-order terms can be neglected. However, depending on the application, one may use higher-order terms in the approximation especially when the forecasting is used at parameters far from the tipping point ($μc$). Now consider that a perturbation is applied to the system. The recovery in Eq. (2) can be written as
$λ(μ,r)=r˙r=p(r)+α1(r) (μ−μc)+α2(r) (μ−μc)2$
(3)
Here, $λ(μ,r)$ is the rate of system recovery from perturbations expressed as
$λ(μ,r)=d (lnr)dt$
(4)

The polynomial used in Eq. (3) determines the minimum number of measurements required for forecasting the bifurcation diagram. Since in this study we use a second-order polynomial, forecasting is possible using at least three measurements in the prebifurcation regime (i.e., recoveries at three values of the bifurcation parameter). That being said, generally the more measurements are available in the prebifurcation regime, the more accurate the forecasting is.

Using measurements of the amplitude of the dynamics, one can employ the following finite-difference approximation to estimate $λ(μ,r)$ (see also Fig. 1) as:

$λ(μ,r)=d (lnr)dt ≅ln r+−ln r− 2Δt$
(5)

where $Δt$ is the time between samples, $r+$ is the value measured at time $t+Δt$, and $r−$ is the value measured at time $t−Δt$. The recovery rate plays the most important role in the forecasting method. The CSD phenomenon at the bifurcation point influences the dynamics near that point also. Specifically, the recovery rate decreases when the system approaches that point. At a fixed system amplitude $r$, the farther the system is from its corresponding postbifurcation regime, the weaker the CSD and the faster the recovery rate of the system is at that amplitude. This is used for predicting the postbifurcation regime (namely the bifurcation diagram).

Following Lim and Epureanu [5,6], suppose that we collect time series during recoveries at several different parameter values $μ1,μ2,...,μn$. At a given parameter value $μk$ (with $k=1…n$), we can choose a value $r=r̃$ and compute $λ(μk,r̃)$ using the measurements. Therefore, n equations are obtained from Eq. (3), one for each $k$, in the form
$λ(μk,r̃)=p(r̃)+α1(r̃)(μk−μc)+α2(r̃)(μk−μc)2$
(6)

The values of $λ(μk,r)$ in Eq. (6) are obtained from measurements using Eq. (5). Fitting a second-order polynomial to these values in the $μ−λ(μ,r̃)$ space, the unknown constants $α1(r̃)$, $α2(r̃)$, and $μc$ can be found. Individual values of $λ(μk,r)$ and the resulting fitted curve are conceptually shown in Fig. 2. The intersection at $μ=μ̃$ of the fitted curve with the $μ$ axis corresponds to a zero recovery rate. That is the most important point. This $μ̃$ value is the forecasted parameter for the corresponding $r̃$ value. This means that a fixed point exists in the postbifurcation regime at $r̃$ when the parameter value is $μ̃$. This procedure can be repeated for different values of $r̃$ and hence, the bifurcation diagram is predicted. This is conceptually shown in Fig. 3. Note that in this procedure the measurements are collected only in the prebifurcation regime, and the bifurcation diagram is predicted without any measurements in the postbifurcation domain.

This procedure can be easily used in the cases that the system recovers from perturbation without oscillations (such as shown in Fig. 1). However, Hopf bifurcations involve oscillations of the system during its recovery from perturbations (schematically shown in Fig. 4). Hence, Eq. (5) cannot be used for all of the data collected in measurements. That is because data points may be adjacent to each other but they do not necessarily have the same phase. The solution used in the past by Lim and Epureanu [5,6] for this problem is to fix the phase by selecting the local maxima of the recovery data and using only those measured values in the same procedure (as for the nonoscillating case). This is a good approximation in the cases where the system oscillates with high frequencies [5,6]. However, in the case of low frequency oscillations, there may not be enough samples available to have a good approximation.

To overcome this limitation, we enhance the forecasting method to obtain recovery rates by nonlinear optimization. Distinct from Lim and Epureanu [5,6], at each parameter value, the recovery rate is assumed to be represented by a polynomial of order $p$ in $r$. One obtains
$λ(μ,r)=λ0(μ)+λ1(μ) r(μ,t)+...+λp(μ)rp(μ,t)$
(7)
where in general the coefficients $λi$ ($i=0...p$) are dependent on the bifurcation parameter. At a fixed parameter value ($μ$), these coefficients are constant over time. Using the definition of $λ(μ,r)=r˙/r$ from Eq. (3), one may rewrite Eq. (7) as
$r˙=r(λ0+λ1r+...+λprp)$
(8)

For a set of given values of the coefficients $λi$ ($i=0...p$), Eq. (8) can be integrated in time to obtain $r(t)$. When the coefficients $λi$ are accurate, the values obtained for $r(t)$ from the integration match the values $rm(t)$ obtained from measurements. In particular, this match has to hold for the local peak values in the measurements. Hence, the coefficients $λi$ can be obtained at each $μ$ value using a nonlinear least squares optimization which minimizes the difference between $rm(t)$ and $r(t)$ for all $t$ values which correspond to local peak responses at that $μ$ value.

Using this method, a more precise approximation of recovery rate $λ$ can be obtained when the time history has low oscillation frequency and only a few peaks. The order $p$ of the polynomial has to be chosen via a convergence test. In the current work, we used a value of four which was sufficient to capture both supercritical and subcritical bifurcations.

### Forecasting Hopf Bifurcations in Large Dimensional Systems.

In the Sec. 2, the forecasting method was developed for a single degree-of-freedom system oscillating during its recovery from perturbations. However, in large dimensional systems, there can be several modes active in the measured system recovery. In the great majority of Hopf bifurcations, including flutter, only one pair of conjugate eigenvalues is involved in the bifurcation (as shown in Fig. 5). Hence, the effects of the other modes can be neglected as they do not exhibit CSD. This is exploited in the proposed forecasting method to enhance the accuracy of bifurcation prediction.

In the aeroelastic systems, the inertial manifold near flutter is generally a two-dimensional nonlinear manifold. This manifold is tangent to the (two-dimensional) center space of the system. Hence, the center space is a good approximation of the inertial manifold near the bifurcation. Thus, we project the dynamics of the system on the center space. In this two dimensional space, we can capture and study the CSD behavior for forecasting when the system is close to the bifurcation. A basis in the center space is obtained as linear modes by using a measurement-based modal decomposition technique. This technique allows us to identify and separate the effects of the mode which exhibits CSD from the other modes. Hence, the dynamics of the system is projected onto the space spanned by the eigenvectors of the bifurcating mode while unwanted modes are filtered out. That leads to more accurate data for the mode which does have CSD, and that enhances the precision of forecasting the postbifurcation regime.

There are several modal identification approaches which can be used to extract mode shapes of a system from measured responses of the system to perturbations. Since we apply the forecasting method close to the bifurcation, the mode corresponding to the pair of eigenvalues with the smallest real part is the one which experiences CSD and is most closely involved in the bifurcation. Here, we use an eigensystem realization algorithm (ERA) [1821]. This ERA is designed for linear modal analysis and provides an approximation for the spatial correlations in the dynamics close to the bifurcation point. This approximation has acceptable accuracy for the aeroelastic system we consider and for the range of perturbations which occur in this system.

ERA is a time domain multi-input multi-output (MIMO) algorithm to identify modal parameters of the system based on the Hankel matrix and its singular value decomposition. The fundamentals of this method can be found in the literature [1821]. Here, only the major steps of the algorithm are presented as related to the aeroelastic system of interest. Specifically, ERA starts with the state-space representation of the system as
$x˙(t)=Ax(t)+Bu(t)y(t)=Cx(t)+Du(t)$
(9)

where $x$ is an $N$ -dimensional vector of state variables, $y$ is an $n$-dimensional vector of measured state variables, $u$ is an $m$-dimensional excitation vector, and $A$, $B$, $C$, and $D$ are state space matrices (constant over time).

Consider that the system is exposed to an initial gust which provides a perturbation/excitation of the type usually occurring during the system operation. The solution of Eq. (9) for the system recovering from this perturbation with $u(t)=0$ from the initial conditions $x0$ can be expressed as
$y(t)=C Φ(t) x0$
(10)
where $Φ(t)$ is the state transition matrix from state $x0$ at time $t=0$ to state $x(t)=Φ(t) x0$ at time $t$. Further consider that the response of the system measured in time and is sampled with a $Δt$ time increment. Equation (10) at time $t=kΔt$ (for $k=0,1,2,...$) can be written as
$yk(t)=C Φk(Δt) x0$
(11)
Repeating the measurements for a sequence of $l$ distinct nonzero initial conditions $X0=[x0,1,x0,2, ... ,x0,l]$ and letting $S=Φ(Δt)$, Eq. (11) can be rewritten as
$Yk(t)=C SkX0$
(12)

where $Yk$ is a matrix with dimension $n×l$.

The objective of the ERA is to find $S$, $X0$, and $C$ such that Eq. (9) is satisfied as closely as possible by using a sequence of experimentally measured $Yk$ (for $k=0,1,2,...$). The measurements are organized in generalized Hankel matrices for each $k$ as follows:
$H(k)=[Y(k)Y(k+1).Y(k+v−1)Y(k+1)Y(k+2).Y(k+v)⋮⋮⋱⋮Y(k+r−1)Y(k+r).Y(k+r+v−2)]nr×lv$
(13)

where $r$ and $v$ are chosen depending on the particular application.

The singular value decomposition of $H(0)$ can be written as $H(0)=P¯Z¯J¯T$, with the singular values in $Z¯$ ordered in decreasing order. Only the first $N$ of them are nonzero. These singular values are grouped in a truncated version of $Z¯$ denoted by $Z$. Matrices $J¯$ and $P¯$ are truncated accordingly to obtain $J$ and $P$ by keeping only their first $N$ columns. It can be shown that one solution for the matrices $S$, $X0$, and $C$ can be expressed as
$S=Z−1/2PTH(1) J Z−1/2 X0=Z1/2JTEl C=EnTP Z1/2$
(14)

where $EnT=[In 0n 0n ... 0n]$ and $ElT=[Il 0l 0l ... 0l]$, while $Ii$ and $Oi$ are identity and zero matrices of order $i$, respectively.

Equation (14) is not a unique solution. For example, for any nonsingular matrix $T$, $S̃=T−1S T$, $X̃0=T−1X0$, and $C̃=CT$ are also a solution. Letting $T$ be the eigenvector matrix of $S$, and transforming the computed matrices to modal coordinates, it can be shown that $C̃$ contains the mode shapes of the system [21]. Furthermore, natural frequencies and damping ratios of continuous system can be obtained as follows:
$ηi=σi±iωd=1Δt ln(s̃i)$
(15)

where $s̃i$ is the ith diagonal element of matrix $S̃$ (which is a diagonal matrix).

## Aeroelastic Model

The proposed forecasting method is applied to an aeroelastic system. The system is modeled as a 2D airfoil oscillating in pitch and plunge as shown in Fig. 6. The free stream velocity is the bifurcation parameter in this demonstration of the forecasting approach. The system is considered to be exposed to gusts which create perturbations. The recovery of the system from these perturbations is measured and used for forecasting.

An actual model of this system is not needed for forecasting. However, a model can provide simulated or surrogate measurement data. Closely following Ref. [22], the equations of motion for this system can be written as:
$ξ¨+xaα¨+2ζξω¯U*ξ˙+(ω¯U*)2 G(ξ)=−1πμCL(τ)−Pg(τ),$

$xarα2ξ¨+α¨+2ζαζξU*α˙+1U*2 M(α)=2πμ rα2CM(τ)+Qg(τ),$
(16)
where $α$ is the pitch angle, $h$ is the plunge displacement, $ξ=h/b$ is the nondimensional plunge displacement, $τ=Ut/b$ is the nondimensional time, and the overdot indicates differentiation with respect to $τ$. $U$ is the free stream velocity, and $U*=U/(bωα)$ is the nondimensional free stream velocity. $ω¯=ωξ/ωα$ is the ratio of the natural frequencies of uncoupled plunge and pitch modes, namely $ωξ$ and $ωa$. The damping ratios for the pitch and plunge are $ζα$ and $ζξ$, and $rα$ is the radius of gyration about the elastic axis. $M(α)$ and $G(ξ)$ are the structural elastic torque and elastic force created by the pitch and plunge motions. They are considered decoupled and nonlinear to account for structural nonlinearities as follows:
$M(α)=m1α+m3α3+m5α5G(ξ)=g1ξ+g3ξ3+g5ξ5$
(17)

where $mi$ and $gi$ ($i=1$, $3$, $5$) are coefficients characterizing the airfoil elasticity in pitch and plunge, respectively.

The effects of the aerodynamics are modeled by $C(τ)L$ and $C(τ)M$ which are the coefficients of lift and moment. These coefficients are given for incompressible flow by Fung [23].

The effects of the gust are modeled by $P(τ)g$ and $Q(τ)g$ which are the lift force and pitch moment due to a gust profile, and are given by the following expressions [24]:
$P(τ)g=2μ∫0τψ˙(τ−σ)wG(σ)dσQ(τ)g=2μ rα2(12+ah)∫0τψ˙(τ−σ)wG(σ)dσ$
(18)

where $wG(τ)$ is the vertical gust velocity distribution, and $ψ$ is the Kussner function [24].

The aerodynamic states are contributing many degrees-of-freedom to the system, and that is reflected in the existence of the integral terms in the model for the aerodynamic forces (such as the terms in Eq. (18) involving the Kussner function). For enhanced computational efficiency, new variables (augmented states) and some simplifications can be introduced [22]. These transform Eq. (15) into a set of eight first order ordinary differential equations as follows:
$x˙=f(x,U*,τ)$
(19)
where $x=[x1,x2,...,x8]T$ is the vector of state space variables, namely $x1=α$, $x2=α˙$, $x3=ξ$, $x4=ξ˙$, $x5=w1$, $x6=w2$, $x7=w3$, and $x8=w4$. The augmented variables $wi$ are defined as [22]
$w1=∫0τe−ε1(τ−σ)α(σ)dσ , w2=∫0τe−ε2(τ−σ)α(σ)dσw3=∫0τe−ε1(τ−σ)ξ(σ)dσ , w4=∫0τe−ε2(τ−σ)ξ(σ)dσ$
(20)

where $ε1=0.0455$ and $ε2=0.3$ are parameters obtained to approximate the properties of the integral terms in the aerodynamics (such as the ones involving the Kussner function).

Equation (19) is used to simulate the dynamic of the system and obtain recovery data from perturbations (gusts) in the prebifurcation regime. These data are used as surrogate measurements to demonstrate the forecasting method. Note that the actual model in Eq. (19) is not needed for the forecasting, but only for generating surrogate measurement data.

## Results and Discussion

In this section, the forecasting method is demonstrated in a numerical example. The parameters used for the aeroelastic system are $rα=0.7$, $ω¯=0.2$, $ah=−0.5$, $xa=0.25$, $μ=100$, $ζα=ζξ=0$. A wind gust of (1 − cos) type (see Fig. 7) is used as perturbation to the system [24]. The vertical gust velocity distribution is zero except for $0≤τ≤2τG$ where it is given by

$wG(τ)=w02(1−cos πττG)$
(21)

where $w0$ is the gust intensity, and $τG$ is the gust gradient expressed in half-chord unit length.

After the airfoil passes the gust, i.e., after $τ=2τG$, the system starts to recover from the perturbation caused by the gust. This response of the system can be obtained by solving Eq. (19) in time, and the results can be used as surrogate data to demonstrate the forecasting method. Also, depending on the structural nonlinearities, the system may face supercritical or subcritical Hopf bifurcations. The forecasting method is applied for each case in Secs. 4.1 and 4.2.

### Supercritical Bifurcations.

For supercritical bifurcations, the coefficients of the elasticity of the airfoil in Eq. (17) are selected as $m1=1$, $m3=1.5$, $m5=0$, $g1=2$, $g3=0$, and $g5=0$. Based on Eq. (3), we discuss results obtained using the lowest number of measurements (i.e., three measurements). Three parameter values ($U*=7.55$, $7.50$, and $7.45$) are chosen in prebifurcation regime, and time histories of the dynamics of the aeroelastic system in response to gust perturbations are obtained by solving the equations of motion numerically. We confirm that the forecasted bifurcation diagrams have good accuracy. Note that increasing the number of measurements does not change these results.

The applied gust parameters are chosen as $w0=0.1$ and $τG=5$ for all three flow speeds. Note that in the proposed method, the perturbations do not have to be all the same and they do not have to be dependent on the bifurcation parameter. For the aeroelastic system that means that the intensity of the gusts which provide the perturbation does not have to be the same and the gusts do not have to be related to the far field flow speed. In addition, perturbations do not have to have the same source (e.g., gust). They could be caused by any perturbation (e.g., maneuvering). However, each perturbation creates an amplitude for the resulting (perturbed) dynamics. The largest forecasted amplitude in the forecasted bifurcation diagram cannot be larger than the amplitude created by these perturbations. For the aeroelastic system, we used only gust perturbations. Also, since the gust intensity can be any value, we used the same value.

To extract the modal parameters of the system, the ERA method is used. We assume that the measured states in the aeroelastic system are pitch, plunge, and their velocities (i.e., $α$, $α˙$, $ξ$, and $ξ˙$). Therefore, the vector $y$ in Eq. (9) is $y(τ)=[α(τ),α˙(τ), ξ(τ), ξ˙(τ)]T$. The parameters $r$ and $v$ are chosen to have a value of 500 in this study, and $Δτ=0.5$.

The modal parameters are only computed at the largest flow speed value and used for all measurements (at all flow speeds). That is consistent with the assumption made in the forecasting method that the inertial manifold varies slowly with the bifurcation parameter. Therefore, we only need to find the center space at the nearest parameter value to the bifurcation point and use that for the entire forecasting procedure.

The forecasting method is applied to obtained time series data of pitch and plunge displacements. Figure 8 shows the recovery rate $λ$ versus pitch amplitude obtained for a fixed flow speed value. Figure 9 shows the plot of $U*−λ$ for the pitch displacement. In this figure, the intersection between each curve and the horizontal axis is the forecasted flow speed $U*$ value in the postbifurcation regime for its corresponding amplitude.

Figure 10 shows the actual and the predicted bifurcation diagrams of pitch and plunge displacements, demonstrating that the method predicts the postbifurcation regime accurately.

To understand the importance of applying the modal decomposition and projection to the center space, the forecasted bifurcation diagram of the pitch displacement when these are not applied is shown in Fig. 11. The prediction is not accurate especially at large amplitudes which are affected more by the modes that are not involved in the bifurcation. The smaller amplitudes have better accuracy since they correspond to later times in the recovery, when the components of the response which do not have CSD are lower (they had time to decay). This effect can be observed in Fig. 12, where the system response of pitch and plunge displacements are compared before and after modal decomposition. The difference between the signals used for forecasting at small amplitudes is negligible. However, at large amplitudes, the difference is larger and affects the forecasting accuracy. Note in Fig. 12(b) that the plunge amplitudes greater than 0.1 cannot be predicted using original data because these values are highly affected by the modes without CSD. These results show the importance of separating modes for Hopf bifurcation predictions.

### Subcritical Bifurcations.

Prediction of subcritical bifurcations has an even greater importance than supercritical ones because it can cause sudden and dramatic change in the LCO amplitude. To demonstrate the forecasting method in such a case, the coefficients of Eq. (17) are chosen as $m1=1$, $m3=−1.5$, $m5=50$, $g1=2$, $g3=0$, and $g5=0$. Then, gust perturbations are applied to the system and data is collected in recovery periods at flow speeds $U*=7.55$, $7.50$, and $7.45$ in the prebifurcation regime. These surrogate measurement data is used to forecast the bifurcation diagram for pitch and plunge displacements as shown in Fig. 13. In this figure, the very largest amplitudes in the bifurcation diagram could not be predicted accurately because these very large values of amplitudes are not present in the recovery time series. In fact, gust perturbations of low intensity are not strong enough to perturb the system into such large amplitudes. Nonetheless, the important fact is that the forecasting method can predict the bifurcation diagram sufficiently to indicate that the system is approaching a subcritical bifurcation.

As the gust intensity is increased, larger perturbations are created and the forecasting method predicts larger ranges of the bifurcation diagram. As an example, the forecasted post-bifurcation regime for pitch and plunge displacements using larger perturbations are shown in Fig. 14.

In summary, the results show that the proposed method is able to forecast both supercritical and subcritical bifurcation diagrams for both pitch and plunge displacements with a good accuracy. When large amplitude perturbations are not available, the method predicts accurately the linear bifurcation point as well as the type of bifurcation (supercritical or subcritical). This is an important advantage of the proposed method compared to other techniques such as the ones based on variance [24].

## Conclusions

A new approach of forecasting supercritical and subcritical flutter (Hopf bifurcations) was presented. The proposed method was developed based on the CSD phenomenon present in oscillating systems near certain types of bifurcations. The recovery rate of the system was extracted using an optimization approach which improves the forecasting accuracy especially when the system has transient recoveries with low frequency oscillations. Moreover, it was shown that applying modal decomposition to the collected transient recovery and keeping only the effects of the mode involved in the bifurcation can increase the forecasting accuracy. In addition, the modal decomposition has the added benefit of alleviating the effects of measurement noise (a feature of most ERA-like methods). The forecasting method relies on local maxima in the measurements, and the accuracy of these measurements is enhanced by the ERA. Future studies are ongoing to analyze and characterize the robustness of the method to measurement and process noise. For noisy measurements, the forecasting accuracy is expected to decrease. When measurement data are contaminated by noise, the number of measurements may be increased to enhance accuracy. The number of flow speeds and the number of measured recoveries at each flow speed can be increased to alleviate the effects of noise and enhance accuracy. For the case of noisy measurements, more than just three flow speeds may need to be measured, depending on the level of noise and the level of desired accuracy.

The method was applied to a nonlinear aeroelastic model in response to gust perturbations. The results show that the method accurately predicts the bifurcation point and also the postbifurcation regime (i.e., the bifurcation diagram) in both supercritical and subcritical cases despite the fact that it uses only prebifurcation regime data.

The technique presented can be used in bifurcation forecasting of large dimensional systems and improving forecasting accuracy especially for larger amplitudes of the bifurcation diagram. This was demonstrated by using a model for the aeroelastic system with eight state variables.

Since bifurcations can cause dramatic changes in system dynamics, this type of forecasting which does not require exploring the postbifurcation regime provides great advantages to a variety of applications especially when safety and maximum system performance is needed.

## Acknowledgment

The authors acknowledge the National Science Foundation (NSF) for the generous support of this work (Grant No. 1334908).

## References

References
1.
Veraart
,
A. J.
,
Faassen
,
E. J.
,
Dakos
,
V.
,
van Nes
,
E. H.
,
Lürling
,
M.
, and
Scheffer
,
M.
,
2012
, “
Recovery Rates Reflect Distance to a Tipping Point in a Living System
,”
Nature
,
481
(
7381
), pp.
357
359
.
2.
Jeffries
,
C.
, and
Wiesenfeld
,
K.
,
1985
, “
Observation of Noisy Precursors of Dynamical Instabilities
,”
Phys. Rev. A
,
31
(
2
), pp.
1077
1084
.
3.
Guttal
,
V.
, and
Jayaprakash
,
C.
,
2008
, “
Changing Skewness: An Early Warning Signal of Regime Shifts in Ecosystems
,”
Ecol. Lett.
,
11
(
5
), pp.
450
460
.
4.
Horsthemke
,
W.
,
1984
,
Noise Induced Transitions
,
Springer
,
Berlin
.
5.
Lim
,
J.
, and
Epureanu
,
B. I.
,
2011
, “
Forecasting a Class of Bifurcations: Theory and Experiment
,”
Phys. Rev. E
83
(
1
), p.
016203
.
6.
Lim
,
J.
, and
Epureanu
,
B. I.
,
2012
, “
Forecasting Bifurcation Morphing: Application to Cantilever-Based Sensing
,”
Nonlinear Dyn.
,
67
(
3
), pp.
2291
2298
.
7.
Scheffer
,
M.
,
Bascompte
,
J.
,
Brock
,
W. A.
,
Brovkin
,
V.
,
Carpenter
,
S. R.
,
Dakos
,
V.
,
Held
,
H.
,
Van Nes
,
E. H.
,
Rietkerk
,
M.
, and
Sugihara
,
G.
,
2009
, “
Early-Warning Signals for Critical Transitions
,”
Nature
,
461
(
7260
), pp.
53
59
.
8.
Lee
,
B. H. K.
,
Jiang
,
L. Y.
, and
Wong
,
Y. S.
,
1999
, “
Flutter of an Airfoil With a Cubic Restoring Force
,”
J. Fluids Struct.
,
13
(
1
), pp.
75
101
.
9.
Lee
,
B. H. K.
,
Price
,
S. J.
, and
Wong
,
Y. S.
,
1999
, “
Nonlinear Aeroelastic Analysis of Airfoils: Bifurcation and Chaos
,”
Progr. Aerosp. Sci.
,
35
(
3
), pp.
205
334
.
10.
Liu
,
L.
,
Wong
,
Y. S.
, and
Lee
,
B. H. K.
,
2000
, “
Application of the Center Manifold Theory in Non-Linear Aeroelasticity
,”
J. Sound Vib.
,
234
(
4
), pp.
641
659
.
11.
Kalmár-Nagy
,
T.
,
Stépán
,
G.
, and
Moon
,
F. C.
,
2001
, “
Subcritical Hopf Bifurcation in the Delay Equation Model for Machine Tool Vibrations
,”
Nonlinear Dyn.
,
26
(
2
), pp.
121
142
.
12.
Jafri
,
F. A.
,
Shukla
,
A.
, and
Thompson
,
D. F.
,
2007
, “
A Numerical Bifurcation Study of Friction Effects in a Slip-Controlled Torque Converter Clutch
,”
Nonlinear Dyn.
,
50
(
3
), pp.
627
638
.
13.
Song
,
Y.
, and
Wei
,
J.
,
2005
, “
Local Hopf Bifurcation and Global Periodic Solutions in a Delayed Predator–Prey System
,”
J. Math. Anal. Appl.
,
301
(
1
), pp.
1
21
.
14.
Mees
,
A.
, and
Chua
,
L. O.
,
1979
, “
The Hopf Bifurcation Theorem and Its Applications to Nonlinear Oscillations in Circuits and Systems
,”
IEEE Trans. Circuits Syst.
,
26
(
4
), pp.
235
254
.
15.
Luo
,
G. W.
, and
Xie
,
J. H.
,
1998
, “
Hopf Bifurcation of a Two-Degree-Of-Freedom Vibro-Impact System
,”
J. Sound Vib.
,
213
(
3
), pp.
391
408
.
16.
,
A.
, and
Epureanu
,
B. I.
,
2015
, “
Forecasting Subcritical and Supercritical Flutter Using Gust Responses
,”
ASME
Paper No. IMECE2015-53105.
17.
D'Souza
,
K.
,
Epureanu
,
B. I.
, and
Pascual
,
M.
,
2015
, “
Forecasting Bifurcations From Large Perturbation Recoveries in Feedback Ecosystems
,”
PloS One
,
10
(
9
), p.
e0137779
.
18.
Hannan
,
E. J.
, and
Quinn
,
B. G.
,
1979
, “
The Determination of the Order of an Autoregression
,”
J. R. Stat. Soc. Ser. B (Methodol.)
,
41
(
2
), pp.
190
195
.
19.
Juang
,
J. N.
,
1987
, “
Mathematical Correlation of Modal-Parameter-Identification Methods Via System-Realization Theory
,”
NASA TM
, Hampton, VA, NASA Technical Memorandum 87720.
20.
Juang
,
J. N.
, and
Pappa
,
R. S.
,
1985
, “
An Eigensystem Realization Algorithm for Modal Parameter Identification and Model Reduction
,”
J. Guid., Control, Dyn.
,
8
(
5
), pp.
620
627
.
21.
Pappa
,
R. S.
,
1994
, “
Eigensystem Realization Algorithm User's Guide for VAX/VMS Computers
,” NASA TM, Hampton, VA, Report No. 109066, p.
1994
.
22.
Lee
,
B. H. K.
,
Gong
,
L.
, and
Wong
,
Y. S.
,
1997
, “
Analysis and Computation of Nonlinear Dynamic Response of a Two-Degree-Of-Freedom System and Its Application in Aeroelasticity
,”
J. Fluids Struct.
,
11
(
3
), pp.
225
246
.
23.
Fung
,
Y. C.
,
2002
,
An Introduction to the Theory of Aeroelasticity
,
Dover Publications
, New York.
24.
Dessi
,
D.
, and
Mastroddi
,
F.
,
2008
, “
A Nonlinear Analysis of Stability and Gust Response of Aeroelastic Systems
,”
J. Fluids Struct.
,
24
(
3
), pp.
436
445
.