Bifurcation Forecasting for Large Dimensional Oscillatory Systems: Forecasting Flutter Using Gust Responses

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 [2–6] (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 [8–10], 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.

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.

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:

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).

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.

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.

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) [18–21]. 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.

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).

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

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

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.

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

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.

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].

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

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.

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

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.

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.

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 [2–4].

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.

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