Output-only parameter identification of a colored-noise-driven Van der Pol oscillator -- Thermoacoustic instabilities as an example

TL;DR

This paper investigates output-only parameter identification of nonlinear oscillators driven by colored noise, demonstrating that accurate parameter estimation is possible with proper filtering, using Van der Pol oscillators and thermoacoustic systems as examples.

Contribution

It introduces a method for accurate parameter identification of nonlinear oscillators under colored noise forcing, challenging the white noise assumption.

Findings

Parameter identification is accurate with band-pass filtering.

Colored noise forcing impacts are mitigated by filtering.

Method applies to thermoacoustic instabilities in combustion chambers.

Abstract

The problem of output-only parameter identification for nonlinear oscillators forced by colored noise is considered. In this context, it is often assumed that the forcing noise is white, since its actual spectral content is unknown. The impact of this white noise forcing assumption upon parameter identification is quantitatively analyzed. First, a Van der Pol oscillator forced by an Ornstein-Uhlenbeck process is considered. Second, the practical case of thermoacoustic limit cycles in combustion chambers with turbulence-induced forcing is investigated. It is shown that in both cases, the system parameters are accurately identified if time signals are appropriately band-pass filtered around the oscillator eigenfrequency.

Equations

$$\ddot{x}+{\omega }_0^{2}x=[2\nu -\kappa x^2]\dot{x}+\xi (t),$$

$$\dot{\xi }(t)=-\frac{1}{{\tau }_{\xi }}\xi (t)+\frac{\sqrt{\gamma }}{{\tau }_{\xi }}\zeta (t),$$

$$H(s)=\frac{\hat{\xi }(s)}{\hat{\zeta }(s)}=\frac{\sqrt{\gamma }}{1+{\tau }_{\xi }s},$$

$$S_{\xi \xi }(\omega )=|H{|}^2{S}_{\zeta \zeta }=\frac{\Gamma }{2\pi }\frac{\gamma }{1+{\omega }^2{\tau }_{\xi }^2},$$

$${\Gamma }_{\text{e}}=2\pi S_{\xi \xi }({\omega }_0)=\Gamma \frac{\gamma }{1+{\omega }_0^2{\tau }_{\xi }^2},$$

$$\gamma =\frac{{\tau }_{\xi }({\omega }_2-{\omega }_1)}{\text{atan}({\omega }_2{\tau }_{\xi })-\text{atan}({\omega }_1{\tau }_{\xi })}.$$

$$k_{\xi \xi }(t)=\Gamma \frac{\gamma }{2{\tau }_{\xi }}{e}^{-\frac{t}{{\tau }_{\xi }}},{\tau }_{\xi }=\frac{1}{k_{\xi \xi }(0)}{\int }_0^{\infty }|k_{\xi \xi }(t)|dt.$$

$$\dot{A}=A\left(\nu -\frac{\kappa }{8}{A}^2\right)+\frac{{\Gamma }_{\text{e}}}{4{\omega }_0^{2}A}+\mu (t),\text{with}⟨\mu {\mu }_{\tau }⟩=\frac{\delta (\tau )\pi S_{\xi \xi }({\omega }_0)}{{\omega }_0^2}=\frac{{\Gamma }_{\text{e}}}{2{\omega }_0^2}\delta (\tau ).$$

$$\dot{A}=A\left(\nu -\frac{\kappa }{8}{A}^2\right)+\frac{\Gamma }{4{\omega }_0^{2}A}+\mu (t),\text{with}⟨\mu {\mu }_{\tau }⟩=\frac{\delta (\tau )\pi S_{\zeta \zeta }({\omega }_0)}{{\omega }_0^2}=\frac{\Gamma }{2{\omega }_0^2}\delta (\tau ).$$

$$P_{\mathrm{ou}}(A)={\mathcal{N}}_{\mathrm{ou}}A\exp \left[\frac{4{\omega }_0^2}{{\Gamma }_{\text{e}}}\left(\frac{\nu A^2}{2}-\frac{\kappa A^4}{32}\right)\right],$$

$$P_{\text{w}}(A)={\mathcal{N}}_{\text{w}}A\exp \left[\frac{4{\omega }_0^2}{\Gamma }\left(\frac{\nu A^2}{2}-\frac{\kappa A^4}{32}\right)\right],$$

$$H=\sqrt{1-B},$$

$${\nabla }^2\hat{p}(s,x)-{\left(\frac{s}{c}\right)}^2\hat{p}(s,x)=-s\frac{(\gamma -1)}{c^2}\hat{q}(s,x)\text{in the domain},$$

$$\frac{\hat{p}(s,x)}{\hat{\mathbf{u}}(s,x)\cdot \mathbf{n}}=Z(s,x)\text{on boundaries},$$

$$\hat{p}(s,x)=s\frac{(\gamma -1)}{4\pi r{c}^2}{\int }_{V_f}\hat{q}(s,y)e^{\frac{s}{c}|x-y|}{d}^{3}y,$$

$$p(x,t)=\sum _{i=1}^{\infty }{\eta }_i(t){\psi }_i(x),$$

$${\ddot{\eta }}_j+{\omega }_j^2{\eta }_j=g_j({\eta }_j,{\dot{\eta }}_j)+{\xi }_j,$$

$$\ddot{\eta }+{\omega }_0^2\eta =[2\nu -\kappa {\eta }^2]\dot{\eta }+\xi (t)$$

$$\hat{\xi }(s)=H(s)\hat{\zeta }(s)=\frac{\sqrt{\gamma }{\tau }^{2}s}{(1+\tau s{)}^2}\hat{\zeta }(s),$$

$$S_{\xi \xi }(\omega )=\frac{\Gamma }{2\pi }\frac{\gamma {\omega }^2{\tau }^4}{(1+{\omega }^2{\tau }^2{)}^2},$$

$$f_{\text{max}}=\frac{1}{2\pi \tau }.$$

$${\Gamma }_{\text{e}}=2\pi S_{\xi \xi }({\omega }_0)=\Gamma \frac{\gamma {\omega }_0^2{\tau }^4}{(1+{\omega }_0^2{\tau }^2{)}^2}.$$

$$\gamma =\frac{2({\omega }_2-{\omega }_1)}{\tau }{\left(\text{atan}({\omega }_2\tau )-\text{atan}({\omega }_1\tau )-\frac{{\omega }_2\tau }{1+{\omega }_2^2{\tau }^2}+\frac{{\omega }_1\tau }{1+{\omega }_1^2{\tau }^2}\right)}^{-1}.$$

$$k_{\xi \xi }(t)=\Gamma \frac{\gamma {\tau }^2}{4\sqrt{2\pi }}(\tau -t)e^{-\frac{t}{\tau }},$$

$${\tau }_{\xi }=\frac{1}{k_{\xi \xi }(0)}{\int }_0^{\infty }|k_{\xi \xi }(t)|dt=\frac{2\tau }{e},$$

$$\dot{A}=A\left(\nu -\frac{\kappa }{8}{A}^2\right)+\frac{{\Gamma }_{\text{e}}}{4{\omega }_0^{2}A}+\mu (t),⟨\mu {\mu }_{\tau }⟩=\frac{\pi S_{\xi \xi }({\omega }_0)\delta (\tau )}{{\omega }_0^2}=\frac{{\Gamma }_{\text{e}}}{2{\omega }_0^2}\delta (\tau ),$$

$$P_{\text{c}}(A)={\mathcal{N}}_{\text{c}}A\exp \left[\frac{4{\omega }_0^2}{{\Gamma }_{\text{e}}}\left(\frac{\nu A^2}{2}-\frac{\kappa A^4}{32}\right)\right],$$

$$\{\begin{array}{l}{\ddot{\eta }}_1+{\alpha }_1{\dot{\eta }}_1+{\omega }_1^2{\eta }_1=[{\beta }_1-{\kappa }_1{\eta }_1^2]\dot{{\eta }_1}+[{\beta }_2-{\kappa }_2{\eta }_2^2]\dot{{\eta }_2}+\xi \\ {\ddot{\eta }}_2+{\alpha }_2{\dot{\eta }}_2+{\omega }_2^2{\eta }_2=[{\beta }_1-{\kappa }_1{\eta }_1^2]\dot{{\eta }_1}+[{\beta }_2-{\kappa }_2{\eta }_2^2]\dot{{\eta }_2}+\xi \\ p=p_1+p_2={\psi }_1{\eta }_1+{\psi }_2{\eta }_2.\end{array}$$

Full text

1 Introduction

System identification (SI) is a long-standing problem that has fostered much research effort Ljung (1999); Pintelon and Schoukens (2012). A wide variety of SI methods have been developed in different frameworks (control theory, machine learning, information theory), and tailored to the specific situation at hand. In each case, the following questions, among others, must be considered to choose the adequate SI method: is it possible to apply a forcing to the system of interest and observe its response (input-output SI), or is it only possible to measure a given observable (output-only or “blind” SI)? Is a model of the system already available, with parameters to be identified (parameter identification), or has the model itself to be uncovered (model identification)? Does the system exhibit nonlinear and/or transient behavior or can linear time invariance (LTI) be assumed? Is the output corrupted by measurement noise? Is the system itself subject to dynamic noise, i.e. external stochastic forcing? Is the system chaotic?

Classical input-output SI techniques generally employ a state-space representation and estimate the parameters of a (postulated or physically derived) model by minimizing the error between the predicted and measured values of some output-based quantity, using e.g. maximum likelihood (ML), prediction error method (PEM), or least-squares (LS) Hamilton (1994); Shumway and Stoffer (2011); Pillonetto et al. (2014); Sovardi et al. (2016). Popular model classes include auto-regressive / moving-average (AR, MA, ARMAX) models Liu et al. (2010), finite impulse response (FIR) models Polifke (2014), output error models (OEM) Ding et al. (2010) and Volterra series King et al. (2016). If physical insight is lacking, SI can take care of selecting an adequate model among several candidates, although a careful trade-off between accuracy and simplicity is needed; this kind of Occam’s razor principle is typically applied with probabilistic (Bayesian) approaches Beck (2010) or sparsity-promoting algorithms Chen et al. (2014). Methods based on machine learning use kernels Pillonetto and De Nicolao (2010) rather than postulating a model in the first place.

Output-only SI methods have to rely on partial information, either because the system cannot be arbitrarily driven, or because the input cannot be measured. Standard tools include Kalman filters Kalman (1960), synchronization methods Yu and Parlitz (2008), modal identification Nagarajaiah and Basu (2009) and reduced-order modeling Rowley and Dawson (2017). Empirical dynamic modeling Ye et al. (2015) allows for model-free output-only SI. As for input-output SI, sparse identification techniques are available for output-only SI Brunton et al. (2016).

Rather than identifying a model or its parameters, some techniques allow the determination of a number of characteristics of a system: distinguish between its chaotic and stochastic nature Zunino et al. (2012), unveil time delays Zunino et al. (2010) or discover hidden patterns Crutchfield (2012) based on information theory (e.g. entropy and complexity); detect causality with convergent cross mapping Sugihara et al. (2012); analyze periodicity and intermittency in noisy signals using recurrence quantification analysis Eckmann et al. (1987); Zbilut et al. (1998); Suresha et al. (2016).

The presence of measurement noise and dynamic noise often complicates the task of SI, deteriorating both its accuracy when identifying parameters and its ability to select plausible models, even though state-space representations can explicitly account for noise. See Reynders and De Roeck (2008); Zhang (2011); Kwasniok (2012) for some efforts towards better noisy SI. However, one can take advantage of the very presence of dynamic noise to extract information and perform output-only SI: inherent stochastic forcing drives the system away from its deterministic equilibrium trajectory and make it visit states that would not been visited otherwise. As proposed in Friedrich et al. (2011), these enriched statistics can then be processed to reconstruct the coefficients of the system’s Langevin equation or corresponding Fokker-Planck equation Risken (1984) and identify the governing parameters.

In the present study this approach is adopted for output-only model-based SI for stochastically driven nonlinear oscillators: the parameters of a given analytical model are identified from the output signal of the system forced by a non-measurable random input. Of course, in the case of linear harmonic oscillators, the system parameters (linear damping rate and resonance frequency) can be readily obtained, e.g. by estimating peak frequency and corresponding quality factor, which is not possible when nonlinearities are active. In this context, accurate and robust output-only parameter identification requires:

1. (i)

an adequate model of the system, 2. (ii)

a model for the driving noise, 3. (iii)

an appropriate data pre-processing.

These aspects are pictured in fig. 1, where a summary of the present work is sketched. In regards to point (i), the selected model for this work is a Van der Pol oscillator (henceforth “VDP”), which is a canonical model used in many different disciplines such as electronics (since the pioneering work Van der Pol (1920)), biology and medicine Van Der Pol and Van der Mark (1928); Jewett and Kronauer (1998); Lucero and Schoentgen (2013), neurology FitzHugh (1961); Nagumo et al. (1962), optics Wirkus and Rand (2002); Barland et al. (2003), seismology Cartwright et al. (1999), sociology and economics Glass and Mackey (1988) or thermoacoustics dynamics in turbulent combustors Noiray and Denisov (2016), the latter being the application discussed in more detail in the second part of the paper. The stochastic differential equation of a Van der Pol oscillator driven by additive noise reads:

[TABLE]

where $x$ represents the state of the system, $f_{0}=\omega_{0}/2\pi$ the natural oscillation frequency, $\nu$ the linear growth rate, $\kappa$ the saturation constant and $\xi(t)$ the additive driving noise.

Concerning point (ii), the simplest model for $\xi$ in eq. 1 is the white noise because it greatly simplifies the analytical derivations. However, a real stochastic forcing is always “colored”, i.e. it always features a non-zero autocorrelation time and a non-constant spectral distribution. One can find a wide collection of studies where the color of the noise plays a fundamental role in the system dynamics, in topics such as economics, biology and mechanical configurations Perelló and Masoliver (2002); Qing-Lin et al. (2015); Sapsis and Athanassoulis (2008), as well as in the specific case of oscillators Masoliver and Porra (1993); Xu et al. (2011); Spanos (1978). In the field of thermoacoustics, one can for instance refer to Tony et al. (2015) or Waugh et al. (2011), the latter investigating the effect of different types of noise on limit-cycle triggering. This suggests that it is essential to take the noise color into account in system identification.

In the first part of the present work, the widely used Ornstein-Uhlenbeck process is used as the driving source of the Van der Pol oscillator. Afterwards another type of noise is introduced for the specific case of thermoacoustic instabilities in turbulent combustors. In both cases, the associated system dynamics and statistics are scrutinised and the effect of noise color on parameter identification is addressed. The need of properly filtering the output data to reliably identify the parameters – item (iii) in the aforementioned list – is then discussed.

2 Van der Pol oscillator driven by Ornstein-Uhlenbeck noise

2.1 Effect of colored noise on oscillations statistics

In this section, the noise that drives the Van der Pol oscillator is generated by an Ornstein-Uhlenbeck (OU) process. It is widely used in various contexts to account for finite correlation time effects of a stochastic forcing. One therefore considers that $\xi$ in eq. 1 satisfies the following Langevin equation:

[TABLE]

where $\zeta$ is a unit-variance Gaussian white noise of intensity $\Gamma$, $\gamma$ is a constant coefficient, which will be used later in the paper to adjust the power of the noise $\xi$, and $\tau_{\xi}$ denotes its characteristic time constant. In the frequency domain, the OU process $\widehat{\xi}$ results from filtering $\widehat{\zeta}$ with the following transfer function

[TABLE]

where $s=i\omega$ is the Laplace variable. The power spectrum of $\xi$ is given by

[TABLE]

It is useful to define the quantity

[TABLE]

which is the power spectral density of $\xi$ at the oscillator eigenfrequency, referred to as “effective OU noise intensity” in the remainder of the paper.

Considering that the target of this study is to quantitatively compare white and colored noise forcing on the oscillator, it is necessary to set a criterion regarding the input power. It is convenient to adjust the intensity of $\xi$ by using the coefficient $\gamma$ such that the powers provided by $\xi$ and by a white noise of intensity $\Gamma$ in a band $[\omega_{1};\omega_{2}]$ are equal, i.e. $\int_{\omega_{1}}^{\omega_{2}}S_{\xi\xi}d\omega=\int_{\omega_{1}}^{\omega_{2}}{\Gamma}/{2\pi}d\omega$, which yields:

[TABLE]

A sensible choice is to define this “iso-power band” around the oscillator resonance frequency $\omega_{0}$: $[\omega_{1};\omega_{2}]=[\omega_{0}-\Delta\Omega;\omega_{0}+\Delta\Omega]$. Henceforth, $\Delta\Omega$ can vary between 0 (band degenerating in the single angular frequency $\omega_{0}$) and $\omega_{0}$ (band $[0;2\omega_{0}]$). The frequency range $\Delta\Omega$ will be referred to as “iso-power semi-bandwidth”. One can see in fig. 2 how this parameter affects the forcing noise power spectrum.

The parameter $\tau_{\xi}$ is a direct measure of how much “colored” the noise is: the shorter $\tau_{\xi}$, the closer to a white noise $\xi$ is. As $\tau_{\xi}$ goes to zero, the cut-off frequency goes to infinity, leading to a constant power spectrum, i.e. a white noise source. This is illustrated in fig. 3 (red spectra), together with the fact that the power spectral density of the oscillator response (blue spectra) is accordingly affected. Note that in the limit $\tau_{\xi}\rightarrow 0$, one gets $\gamma\rightarrow 1$ and $S_{\xi\xi}(\omega)\rightarrow{\Gamma}/{2\pi}=S_{\zeta\zeta}(\omega)$. The characteristic time $\tau_{\xi}$ is the noise correlation time, obtained via the autocorrelation function of $\xi$:

[TABLE]

Such OU process is now considered as being the driving force of the Van der Pol oscillator given by eq. 1. It is convenient to investigate the system in terms of its slowly-varying amplitude and phase dynamics with $x(t)\approx A(t)\cos{[\omega t+\varphi(t)]}=A(t)\cos{\phi(t)}$. This coordinate change is legitimate provided that $\nu\ll\omega_{0}$. Performing deterministic and stochastic averaging Stratonovich (1967) yields the following stochastic differential equation for the amplitude $A$:

[TABLE]

It is important to underline that the averaging procedure is valid only if $\tau_{\xi}\ll\tau_{A}$, where the amplitude correlation time $\tau_{A}$ is related to the system growth rate by $\tau_{A}\simeq\pi/|\nu|$ (see Lax (1967); Lax et al. (2006); Noiray (2016)). One can refer to fig. 4 where the important time scales of the considered system are presented. It is also interesting to compare eq. 8 to its white-noise-driven oscillator counterpart

[TABLE]

The two equations only differs by the fact that $\Gamma_{\text{e}}$ substitutes $\Gamma$. In the limit $\tau_{\xi}\rightarrow 0$, $\Gamma_{\text{e}}\rightarrow\Gamma$, and eq. 8 tends to eq. 9.

Considering the Fokker-Planck equation associated with eq. 8, one can derive the stationary probability distribution (PDF) for the amplitude of the VDP oscillator driven by an OU noise:

[TABLE]

and for the white-noise driven VDP oscillator:

[TABLE]

where $\mathcal{N}_{\mathrm{ou}}$ and $\mathcal{N}_{\mathrm{w}}$ are two normalization constants such that $\int_{0}^{\infty}P(A)dA=1$.

Apart from the normalization constants, $P_{\mathrm{ou}}$ and $P_{\mathrm{w}}$ differ by the factor $\Gamma/\Gamma_{\text{e}}={(1+\omega_{0}^{2}\tau_{\xi}^{2})}/{\gamma}$ in the exponential, which is depicted in fig. 5.

Figure 6 compares the amplitude PDFs of the oscillator driven by white noise (shaded area) and by the colored noise (solid lines) for the same system parameters $\nu$, $\kappa$, $\omega_{0}$ and $\Gamma$. Two different iso-power semi-bandwidths $\Delta\Omega$ (columns), which were already considered in fig. 2, as well as two different values of the linear growth rates $\nu$ (rows) are considered. In the case of a wide iso-power band, one can observe that $P_{\mathrm{ou}}$ significantly deviates from $P_{\mathrm{w}}$ when $\tau_{\xi}$ increases. One can note that for large enough $\tau_{\xi}$ and for $\Delta\Omega<\omega_{0}$, $P_{\mathrm{ou}}$ tends to a limit case distribution111It can be proven that $\lim_{{\tau_{\xi}}\to\infty}(1+\omega_{0}^{2}\tau_{\xi}^{2})/\gamma=\omega_{0}^{2}/(\omega_{0}^{2}-\Delta\Omega^{2})$, so except for the case $\Delta\Omega=\omega_{0}$, this limit is finite and the $P_{\mathrm{ou}}$ asymptotically tends to a limit PDF. Remember that $\tau_{\xi}\ll\tau_{A}$ must anyway hold to have a valid derivation of the equations.. On the other hand, no significant difference among the PDFs can be noticed when $\Delta\Omega/\omega_{0}=0.1$.

To obtain a quantitative measure of the difference between the two PDFs, one can make use of the Hellinger distance:

[TABLE]

where $B=\int_{-\infty}^{+\infty}\sqrt{p(x)q(x)}\mathrm{d}x$ is the Bhattacharyya coefficient. The Hellinger distance $H$ is a statistic quantity that measures the difference between two PDFs of the same random variable $p(x)$ and $q(x)$, and ranges from 0 when $p(x)=q(x)$, to 1 when they do not overlap. In the following, $H$ is computed to compare $P_{\text{w}}$ and $P_{\text{ou}}$ in a systematic way for different points $(\Delta\Omega,\tau_{\xi},\nu)$ of the space of iso-power semi-bandwidth, correlation time and growth rate. The results are presented as colormaps in fig. 7.

The linear growth rate $\nu$ has a minor effect: all the maps in fig. 7.a are similar, but $H$ is slightly higher when $\nu<0$, due to the shift of the amplitude of maximum probability $A_{\mathrm{m}}$ observed in this case (see again fig. 6). Focusing on the other two parameters in fig. 7.b, $H$ is large in the upper-right corner of the map, i.e. for high values of $\Delta\Omega$ and $\tau_{\xi}$. The influence of $\tau_{\xi}$ is intrinsically related to the noise color: as discussed earlier, the shorter $\tau_{\xi}$, the closer is $\xi$ to a white noise. That is why the region of match between $P_{\mathrm{c}}$ and $P_{\mathrm{w}}$ is wider for short correlation times. In case of a long $\tau_{\xi}$, the bandwidth $\Delta\Omega$ has a strong influence, leading for large values to a significant difference between $P_{\mathrm{c}}$ and $P_{\mathrm{w}}$. A large $\Delta\Omega$ means that the equality of power between white noise and colored noise is set in a wide band around the oscillator eigenfrequency. If $\tau_{\xi}$ is long enough to let the oscillator frequency $f_{0}$ fall in the decaying part of $S_{\xi\xi}$, the power spectral density of the two forcing noise is sensibly different around $f_{0}$ (see again fig. 2), and the response of the system significantly changes.

2.2 Parameter identification and white-noise assumption

In this section, the influence of the finite correlation time $\tau_{\xi}$ of the driving OU noise upon parameter identification strategies is investigated. The problem is the following: the noise driving the oscillator is never white in practice. Therefore, the use of a white noise driven oscillator model as a base for parameter identification can be brought into question.

One alternative would be to adopt a model featuring a noise source with finite correlation time as exemplified in the previous section with the OU process. However, this would not make any difference if the adopted SI method relies on the statistics of the signal. In fact, looking at eq. 10, one can see that the analytical expression for $P_{\text{ou}}$ produces self-similar probability distributions. In other words, different combinations of $\Gamma$ and $\tau_{\xi}$, lead to the same output amplitude statistics. However, if one is only interested in identifying the linear growth rate $\nu$ and the saturation coefficient $\kappa$ one should presumably be able to use a white noise driven VDP model as a basis for the SI. Still, it has to be verified if the presence of non-zero autocorrelation time $\tau_{\xi}$ can affect the identification process: even though the amplitude PDFs of the two models are the same, the output time traces and spectra are different, especially for some combinations of parameters.

To verify the possibility of achieving a robust parameter identification of the linear growth rate $\nu$ and the saturation coefficient $\kappa$ using a white noise approximation, the following test is performed. A Van der Pol oscillator (see eq. 1), having the true parameters $\nu=\nu_{\text{t}}$ and $\kappa=\kappa_{\text{t}}$ and forced with an OU noise of intensity $\Gamma_{\text{t}}$ and correlation time $\tau_{\xi}$, is simulated in Simulink®, and then the slowly-varying envelope $A(t)$ and phase $\varphi(t)$ of the output signal $x(t)$ are extracted.

A parameter identification using the white noise driven model is then attempted, making use of the approaches 3 and 4 proposed in Noiray and Schuermans (2013a). They consist in finding the optimum parameters $\nu$, $\kappa$ and $\Gamma$ giving the best fit of $P(A)$ and $P(A\dot{\varphi})$ for method 3, and of the drift and diffusion coefficients of the Fokker-Planck equation for method 4. However, the identified parameters significantly differ from the actual values: $\nu_{\text{id}}=2.1\nu_{\text{t}}$, $\kappa_{\text{id}}=1.9\kappa_{\text{t}}$, $\Gamma_{\text{id}}=1.7\Gamma_{\text{t}}$ with approach 3 and $\nu_{\text{id}}=1.5\nu_{\text{t}}$, $\kappa_{\text{id}}=1.6\kappa_{\text{t}}$, $\Gamma_{\text{id}}=1.5\Gamma_{\text{t}}$ using the approach 4.

As will become apparent, the parameter identification failed because of the lack of pre-processing of the data. In fact it is wrong to assume that the measured output spectrum $S_{xx}(\omega)$ can be generated by an equivalent white noise source: the actual driving noise spectral power distribution $S_{\xi\xi}(\omega)$ leaves some peculiar signature in $S_{xx}(\omega)$. However, it is indeed possible to reproduce over a limited band around the oscillator frequency the actual output of the colored noise driven VDP with a white noise forcing, because $S_{\xi\xi}$ is a smooth function of frequency. This is exemplified in fig. 8.b, where one can see the spectrum of a colored noise driven VDP (thick grey line), overlaying the one of a VDP driven by a white noise of intensity $\Gamma_{\text{e}}$ (thin black line).

The next attempt is, therefore, to bandpass filter the signal obtained from the simulation in the band $f_{0}\pm\Delta f$, using a $\Delta f$ progressively narrower222Note that $\Delta f$ is not related to $\Delta\Omega$: the first is the filter semi-width adopted to pre-process the data for parameter identification, the second is a semi-bandwidth arbitrarily chosen to define the driving noise intensity.. The obtained identification of $\nu$ is presented in fig. 8.a as a function of $\Delta f$. If $\Delta f=f_{0}$, the identified parameters values are close to the ones obtained using the unfiltered data. Decreasing $\Delta f$, the identified growth rate $\nu_{\text{id}}$ converges to the actual one $\nu_{\text{t}}$ for $\Delta f/f_{0}=0.3$. The same trend is found for the saturation constant $\kappa$. This indicates that it is necessary to filter the data around the frequency of interest in order to perform a reliable model-based output-only parameter identification.

One might be tempted to reduce further the filter bandwidth, in order to decrease even more the driving noise modeling inaccuracy. However, one can see that below $\Delta f=0.05$ the estimated $\nu$ again deviates from the actual one. This fact is explained through the other two panels of fig. 8. In the panel b, the spectrum of the signal generated by the simulation of the OU noise driven VDP oscillator is presented, together with two different filter widths. The corresponding filtered time traces of the oscillation amplitude, used as data for the parameter identification, are plotted in panel c, superimposed to the unfiltered oscillator signal (grey). One can observe that if a too narrow band is considered, the signal is altered and substantially deviates from the original: the amplitude time trace follows the general trend, but does not capture anymore the high frequency content. This affects the statistics and dynamics of the data and, therefore, the outcome of the parameter identification. Hence, one must refrain from filtering too much the signal, to preserve the core information of the original signal.

In the next step, the parameter identification is performed for different colored noise parameters, to ensure that an adequate filtering is the means of achieving a reliable identification. In fig. 9 the results of this test are presented. Each panel includes the identification result (method 4 in Noiray and Schuermans (2013a) is adopted) of 100 different simulations of the system, each corresponding to a different combination of noise parameters $\Delta\Omega$ and $\tau_{\xi}$. The identification inaccuracy is given in terms of relative error $\varepsilon=|\nu_{\text{t}}-\nu_{\text{id}}|/\nu_{\text{t}}$. In the left column, the identification results when using the raw data are presented. The iso-power semi-bandwidth $\Delta\Omega$ does not noticeably affect this error, as it just changes the value of $\Gamma_{\text{e}}$ to be identified. The noise correlation time has a dramatic impact on the identification error for long correlation times. However, the error vanishes if $\tau_{\xi}$ is very short, as in this case the driving noise gets closer a white one. In the right column of fig. 9, the same signals are bandpass filtered in a band $[f_{0}(1\pm 0.5)]$ before the parameter identification is run. One can notice how the identification is considerably enhanced, leading to very accurate results regardless of the parameters of the noise source. This result consolidates the confidence on the output-only parameter identification methods, as even without knowing the noise parameters $\Gamma$, $\Delta\Omega$ and $\tau_{\xi}$ it is possible to obtain the correct oscillator parameters $\nu$ and $\kappa$ by just applying an adequate filter to the output signal.

Summing up, it can be stated that, for a OU noise driven VDP oscillator, the parameter identification based on a white noise approximation will accurately estimate the linear growth rate $\nu$ and the saturation constant $\kappa$ if the signal is filtered before the analysis. The filtering bandwidth has to be:

• •

narrow enough, to have a satisfactory approximation of the real noise with a white one over the considered band,

• •

not too narrow, to preserve the amplitude dynamics of the signal.

A sensible strategy is to use a progressively narrower filter for the data pre-processing, and repeat the identification process until the obtained parameters reach a plateau.

In the next part of this work, the study will be carried out using a different type of noise source, which is also closer to the actual stochastic forcing characteristic of thermoacoustic systems.

3 Thermoacoustic instabilities: Modeling

3.1 Practical context

In gas turbine, aeronautics and aerospace applications, the race for more efficient, less polluting, more fuel- and operation-flexible systems is ongoing, towed by customers needs and environmental regulations Lieuwen (2012). The thermoacoustic instabilities taking place in the combustion chambers of these engines constitute a major difficulty to overcome Poinsot (2016), because their resulting high amplitude acoustic levels induce high cycle fatigue of the combustor components and reduce their lifetime. The mechanisms ruling the constructive interaction between flames and acoustic modes are complex and the occurrence of these instabilities at a given engine operating point is hard to predict.

Therefore, the development of reliable predictive methods is of primary importance. Currently, brute force Large Eddy Simulations cannot be routinely used in a combustor design optimisation context due to their prohibitive computational costs. Therefore, a significant portion of the research efforts concentrate on the development of Helmholtz solvers and low-order thermoacoustic network models that are combined with experiments or computationally-cheaper numerical simulations Schuermans et al. (2010); Han et al. (2015); Silva et al. (2017); Nicoud et al. (2007); Bourgouin et al. (2015); Oberleithner et al. (2015); Campa and Camporeale (2014); Schmid et al. (2013); Ghirardo et al. (2015).

In parallel, it is also important to establish robust system identification methods in order to validate the aforementioned linear-stability prediction tools. It has been recently shown that thermoacoustic linear growth rates can be extracted from limit cycle dynamic pressure data recorded in real systems Noiray and Schuermans (2013a, b); Noiray (2016); Noiray and Denisov (2016), and compared to the ones obtained using predictive thermoacoustic methods. Such network model validation is performed in Bothien et al. (2015).

In the context of the present work, this section deals with output-only parameter identification methods applied to thermoacoustic systems, where the measurable output is the acoustic pressure at one location in the combustion chamber while the unknown input is the stochastic forcing resulting from the intense turbulence in the combustor. This last contribution is often modelled as an additive forcing, and assumed to be a white noise. In reality, this noise is not delta-correlated as explained in section 3.2. Therefore, in section 3.3 a more accurate model of the actual noise is introduced and the equations for the Van der Pol forced by this specific noise source are derived. The impact of the selected model on the effectiveness of the parameter identification is afterward scrutinised in section 4.1.

Regarding the system modeling, a single thermoacoustic mode description is often adopted in order to keep the number of system parameters to be identified to a minimum. This allows the use of a single oscillator as a model of the system. However the raw data, i.e. the acoustic pressure at a given location in the combustion chamber, result from the superposition of the contributions from all the combustor eigenmodes. Consequently, it cannot be directly treated and requires pre-processing to isolate the information corresponding to the single eigenmode considered for parameter identification. This can be done by bandpass filtering the data Noiray and Denisov (2016) or by performing a modal projection if several simultaneous records at different locations in the chamber are available Noiray and Schuermans (2013b). These data manipulations can, however, change the outcome of output-only parameter identification methods, because the signal and its statistics can be sensibly altered. This problem is considered in section 4.2.

3.2 Colored random excitation

In thermoacoustics, the acoustic pressure satisfies the Helmholtz equation with heat release rate source in the volume of the domain and the impedance conditions on boundaries:

[TABLE]

[TABLE]

where $\widehat{p}$ and $\widehat{u}$ are the acoustic pressure and velocity fluctuations, $s$ the Laplace variable, $x$ the position, $c$ the local speed of sound, $\gamma$ the specific heat ratio, $\widehat{q}$ the heat release rate fluctuation, $\mathbf{n}$ the outward normal to the boundary and $Z$ the acoustic impedance. This equation stands if the Mach number is low. If the flame is placed in an open environment, waves generated by the reaction zone are radiated away without reflections. In reference Hirsch et al. (2007), the radiated sound field in this situation is modelled as function of the turbulence-induced heat release rate fluctuation and compared to experimental data. The formal solution of eq. 13 for a fluctuating heat release rate source in an open environment is:

[TABLE]

where $x$ is the observer position in the far field, $r\approx|x|$ is the distance of the observer from the flame. This equation is valid when the flame brush, which extends over the volume $V_{f}$, is compact with respect to the considered acoustic wavelength. An example of the far-field acoustic power spectral density $S_{pp}$ in such configuration, i.e. the so-called combustion noise Strahle (1978), is given in fig. 10, together with integrated heat release oscillation power spectrum $S_{qq}$. In this situation, the heat release rate fluctuations $\widehat{q}$ generating the sound field are only due to the non-coherent turbulent component $\widehat{q}_{\text{n}}$ (see fig. 11.a).

The combustion noise spectrum $S_{pp}$ features a maximum at frequency $f_{\mathrm{max}}$ and a bandpass signature, in contrast with the low-pass character of $S_{qq}$, having $f_{\mathrm{max}}$ as cut-off frequency. The two spectra are related to each other by eq. 15, which is the topic of e.g. Rajaram and Lieuwen (2009); Ihme et al. (2009).

All the authors, from the fundamental theoretical work by Clavin and Siggia Clavin and Siggia (1991) to the systematic study by Rajaram and Lieuwen Rajaram and Lieuwen (2009), agree on the shape of the combustion noise spectrum $S_{pp}$. In Rajaram and Lieuwen (2009) it has been shown that the normalized combustion noise power spectra of different burners operating under different conditions collapse on top of each other, indicating a general scaling law (see fig. 12.a). The combustion noise spectrum features a maximum, and varies like $S_{pp}(f)\propto f^{2}$ on the left side, and like $f^{-r}$, with $2<r<3.4$, on the right side. The peak frequency of the combustion noise spectrum can be estimated making use of experimental relations such as the one proposed in Shivashankara et al. (1975), involving dimensions, flow properties and chemical quantities. Alternatively, it has been observed in Rajaram and Lieuwen (2009) that the Strouhal number $St=f_{\mathrm{max}}{L_{\mathrm{F}}}/{U_{\mathrm{avg}}}$ is almost in any case close to 1, where $L_{\mathrm{F}}$ is the flame length and $U_{\text{avg}}$ the average velocity of the reactants mixture. Hence $f_{\mathrm{max}}\approx{U_{\mathrm{avg}}}/{L_{\mathrm{F}}}$, which is shown in the inset of fig. 12.a.

As exemplified in fig. 12, the acoustic signature dramatically changes when the flame is placed within a combustion chamber. In fig. 12b, a single mode dominates the spectrum, with a sharp peak of frequency $f_{0}$, surrounded by several side peaks, which correspond to the other thermoacoustic eigenmodes. One can conveniently express the acoustic pressure at a given location $x$ as

[TABLE]

where $\eta_{i}$ denotes the amplitude of the $i^{th}$ mode and $\psi_{i}$ the spatial distribution of the corresponding natural acoustic mode of the chamber.

This spatial projection leads to a set of coupled stochastic nonlinear differential equations for the modes $\boldsymbol{\eta}(t)=\left[\eta_{1}(t),\cdots,\eta_{j}(t),\cdots\right]^{T}$. However, it is often possible to describe the dynamics of a single mode $j$ by neglecting the influence of other modes Culick (2006). In this case, the mode amplitude $\eta_{j}$ is given by the nonlinear stochastic oscillator

[TABLE]

where $\omega_{j}=2\pi f_{j}$ is the angular frequency of the $j^{\text{th}}$ natural acoustic mode, $\xi_{j}(t)$ is the additive stochastic forcing coming from turbulence-induced processes. The term $g_{j}$ is a non-linear function which includes, amongst others, the effects of acoustic damping mechanisms and coherent heat release rate fluctuations (this last contribution is coherent in the sense that it depends on the acoustic field). One can see in fig. 11.b a diagram depicting the coherent feedback $\widehat{q}_{\text{a}}$ from a flame located in a combustion chamber. At the same time the flame is also influenced by the turbulent flow. The resulting heat release fluctuation is the aforementioned $\widehat{q}_{\text{n}}$, which was the only source in case of an open flame. The turbulence-induced flow perturbations exhibit a much smaller spatial correlation than the acoustic ones, which are correlated over the entire combustor. These quantities $\widehat{q}_{\text{a}}(\omega)$ and $\widehat{q}_{\text{n}}(\omega)$ can be measured in dedicated test rigs equipped with loudspeakers and microphones, as explained in e.g. Paschereit et al. (2002), and can be used afterwards in network models providing predictions of the system stability.

3.3 Colored noise driven Van der Pol oscillator

In the following, it is assumed that the non-linear function $g_{j}$ in eq. 17 results from a linear acoustic damping and a non-linear flame feedback: $g({\eta},\dot{\eta})=\dot{q}_{\mathrm{a}}-\alpha\dot{\eta}$, where $\alpha$ is the damping constant, and subscripts are omitted from now on. The flame response is expanded up to the third order in acoustic amplitude, which is often sufficient to characterise supercritical thermoacoustic bifurcations Lieuwen (2003); Boujo et al. (2016): ${q}_{\mathrm{a}}=\beta\eta-{\kappa}\eta^{3}/3$. This assumption yields the already presented Van der Pol oscillator equation:

[TABLE]

where $\nu=(\beta-\alpha)/2$ is the linear growth rate.

Considering $\xi$ as a white noise, i.e. a delta correlated forcing, simplifies the modeling approach and has been used in most of the studies dealing with stochastically forced thermoacoustic limit cycles (e.g., again Lieuwen (2003); Boujo et al. (2016)).

In the remainder of the paper the random forcing $\xi$ is assumed to result from the non-coherent heat release rate fluctuations $q_{\text{n}}$ only. As a result, $S_{\xi\xi}$ follows the same power law as the combustion noise and is therefore proportional to $S_{pp}$ Ihme et al. (2009); Liu and Echekki (2015).

In order to keep the problem tractable, $\xi$ is defined by

[TABLE]

where $\zeta$ is a unit-variance Gaussian white noise of intensity $\Gamma$, $\gamma$ is a constant used to adjust the power of the process $\xi$ and $\tau$ is its characteristic time constant. The resulting power spectrum is given by $|H|^{2}S_{\zeta\zeta}$:

[TABLE]

that features a maximum at

[TABLE]

One can again define an “effective colored noise intensity”:

[TABLE]

This model is a close approximation of actual experimental data, as shown in fig. 10 (-$\Circle$-). This model is also close to others provided in literature, like in Liu and Echekki (2015), but, thanks to its simplicity, it allows for the analytical derivation that follows.

As done for the OU case, the colored noise power is equated to the one of a white noise of intensity $\Gamma$ in the band $[\omega_{1};\omega_{2}]=\omega_{0}\pm\Delta\Omega$, which yields:

[TABLE]

One can see in fig. 13 how the parameter $\Delta\Omega$ affects the forcing noise power spectrum.

The characteristic time $\tau$ is related to the noise correlation time $\tau_{\xi}$, that can be obtained via the autocorrelation function of $\xi$,

[TABLE]

[TABLE]

where $e=\exp{(1)}\simeq 2.718$.

The value of $\tau_{\xi}$ is related to the “color” of the noise: it determines where the maximum of the noise spectrum $f_{\text{max}}$ is located compared to the oscillator eigenfrequency $f_{0}$, affecting, as presented in fig. 14, the response of the VDP. Focusing in a band around $f_{0}$, one can see how the oscillator is forced either by a source having the power increasing with frequency (“blue” noise), or almost constant (close to a white noise), or decreasing (“pink” noise). The resulting output $p$ is, accordingly, substantially different.

The VDP equation is again recast in amplitude-phase coordinates. In this case, this substitution is legitimate as, in most of the practical cases, the thermoacoustic systems satisfy the condition $\nu\ll\omega_{0}$. This means that the right hand side of eq. 17 is much smaller than the left one and then $\eta(t)\approx A(t)\cos{[\omega t+\varphi(t)]}=A(t)\cos{\phi(t)}$. Adopting the colored noise model (20) for $\xi$, deterministic and stochastic averaging yields the stochastic differential equation:

[TABLE]

with $\Gamma_{\text{e}}$ given by eq. 22. Again, the averaging method is valid if the correlation times are such that $\tau_{\xi}\ll\tau_{A}$ Stratonovich (1967). This is generally verified for practical cases. The amplitude correlation time is related to the growth rate by $\tau_{A}\simeq\pi/|\nu|$ Noiray (2016). Taking $\nu=10$ rad/s, for instance, $\tau_{A}=314$ ms, while the noise correlation time $\tau_{\xi}=2\tau/e\approx 1/e\pi f_{\mathrm{max}}$ is generally smaller than 1 ms ($f_{\mathrm{max}}\geq 50$ Hz, see fig. 12.a).

The stationary probability distribution for the amplitude of the bandpass noise driven VDP oscillator is then:

[TABLE]

where $\mathcal{N}_{\mathrm{c}}$ is the normalization constant to have $\int_{0}^{\infty}P_{\text{c}}(A)dA=1$.

Like for the OU case, eqs. 26 and 27 have the same structure as their white noise driven system counterparts eqs. 9 and 11, with the effective colored noise intensity $\Gamma_{\text{e}}$ (eq. 22) replacing the white noise intensity $\Gamma$. Therefore, $P_{\mathrm{c}}$ and $P_{\mathrm{w}}$ differ only for the factor $\Gamma/\Gamma_{\text{e}}$ in the exponential. In fig. 15 one can see a map of this factor, as function of the iso-power bandwidth $\Delta\Omega$ and of the source noise correlation time $\tau_{\xi}$.

Comparing this map with the one for the OU noise (fig. 5), one can notice how $\Gamma$ and $\Gamma_{\text{e}}$ might differ whatever the correlation time. This is due to the fact that this type of noise does not converge to a white one for short $\tau_{\xi}$. Another difference is that this coefficient can be lower than one.

In line with this map, $P_{\mathrm{c}}$ and $P_{\mathrm{w}}$ can show significant differences, as depicted in fig. 16. To compare quantitatively $P_{\text{c}}$ and $P_{\text{w}}$, the Hellinger distance is plotted in fig. 17. As for the OU noise, for small $\Delta\Omega$, $H$ tends to 0. However in this case, for large $\Delta\Omega$, $H$ is large whatever the correlation time of the noise source.

4 Thermoacoustic instabilities: parameter identification

In this section, the white noise approximation is assessed in the context of parameter identification. As discussed before, the dynamics of a thermoacoustic mode can be seen as a SISO system. Although the output, represented by pressure oscillation, is easily accessible via experimental measures, the input, resulting from turbulence, is not known. Therefore, this system necessitates output-only parameter identification methods.

4.1 Assessment of the white noise approximation

Following the same procedure as in section 2.2, 100 test cases with fixed oscillator parameter $\nu=\nu_{\text{t}}$ and $\kappa=\kappa_{\text{t}}$, but different noise parameters $\tau_{\xi}$ and $\Delta\Omega$, are analysed to ensure that the identification methods relying on the white noise assumption are not biased by the actual noise spectrum and autocorrelation. The relative error $\varepsilon=|\nu_{\text{t}}-\nu_{\text{id}}|/\nu_{\text{t}}$ on the estimated oscillator linear growth rate $\nu_{\text{id}}$ is presented in fig. 18.

Like for the OU noise case, the identification might fail if the unfiltered data are used (left panel). It is interesting to notice that, compared to the OU case, the error is generally less severe. This is due to the spectral distribution of the bandpass noise, rapidly decaying in power at high and low frequencies. Another peculiar aspect is the distribution of the errors in the map. While for the OU noise low $\tau_{\xi}$ means a quasi-white noise forcing and, therefore, a small identification error, here short $\tau_{\xi}$ corresponds to a blue noise forcing.

As in the OU case, filtering the data prior to parameter identification improves the identification results (right panel of fig. 18). This is, again, due to a more accurate approximation of the real forcing spectrum with a white one in the considered frequency range.

The bandwidth of the filter adopted in the pre-processing of the data has to be chosen with care in order not to discard essential amplitude dynamics. In addition, practical acoustic spectra often feature neighboring peaks around the main one, due to the coexistence of several thermoacoustic modes in the combustor. This is a further constraint when one analyses experimental data and performs single-mode output-only parameter identification. These two aspects are covered in the following.

4.2 Effect of data preprocessing on parameter identification

A typical combustor acoustic pressure spectrum features several peaks (fig. 12.b). The different modes acting in the domain are mutually coupled, each one influencing the response of the others. However, if the neighboring peaks are not too close, one can analyse one mode at a time, isolating its dynamic from that of the other modes. This is easily done by bandpass filtering the data and simplifies the system identification, since neither the parameters of neighboring modes, nor the coupling coefficients have to be taken into account.

Figure 19 shows a typical situation and the effects of a different filter bandwidth. A wider portion of this spectrum has already been shown in fig. 12.b. This experimental spectrum features a strong peak, corresponding to the dominant mode eigenfrequency, surrounded by two others small peaks. In order to identify the mode parameters accurately, removing the other modes effect, the signal is filtered around the main peak, i.e. in the band $[f_{0}-\Delta f;f_{0}+\Delta f]$. The maximum bandwidth is the one that discards neighboring peaks while keeping the main peak and its tails ($\Delta f/f_{0}=0.20$ in this case). One could also choose narrower bands ($\Delta f/f_{0}=0.10$ or $\Delta f/f_{0}=0.025$ in this example), obtaining different resulting time signals. Looking at the central panels of fig. 19, one can see that in the first case (green), the dynamics on time scales comparable to the amplitude correlation time $\tau_{A}=\pi/|\nu|$ is preserved: compared to the widest filter (blue), only high-frequency amplitude oscillations are lost. This means that the essential dynamics are unaffected. In the second case (red), the general trend is followed, but too much information has been lost to reliably identify the parameters $\nu$ and $\kappa$.

In the following, a “toy model” of two coupled oscillators driven by colored noise is used to illustrate this issue:

[TABLE]

The total output $p$, which is the sum of the outputs of the two oscillators $\eta_{1}$ and $\eta_{2}$ weighted by $\psi_{1}$ and $\psi_{2}$, features a spectrum, plotted in fig. 20, that is similar to the experimental pressure spectrum shown in fig. 12. This figure also highlights what is hidden behind a single-mode approximation.

Note the difference between the spectra of $p_{2}$ without coupling (theoretical, thick blue) and with coupling (numerical, thin blue), especially for $f=f_{1}$. This difference appears because the oscillators are coupled and the linearly unstable oscillator #1, characterized by a limit-cycle at $f_{1}=\omega_{1}/2\pi$, is forcing oscillator #2, having eigenfrequency $f_{2}=1.3f_{1}=\omega_{2}/2\pi$. At the same time, the linearly stable mode (oscillator #2) contributes to $S_{pp}$ around the eigenfrequency $f_{1}$ of the unstable mode (oscillator #1). Therefore, the response of the system at $f=f_{1}$ is not due to the oscillator #1 only. However, if the two peaks are distant enough and one is stronger than the other, these mutual contributions are negligible, compared to the direct output of the oscillator #1 at its natural frequency (more than 20 dB of difference in this example). Restricting the discussion to this case, one can adopt the aforementioned single-mode approximation, and attempt a parameter identification on one mode at a time.

To test the sensibility of the identification results to the filter bandwidth, the output signal is filtered with different bandwidth around the first eigenfrequency $f_{1}$. The aim is to extract the linear growth rate of the first unstable oscillator, which has the true value $\nu_{\text{t}}\approx(\beta_{1}-\alpha_{1})/2$. For this purpose, three different methods of Noiray and Schuermans (2013a) are used. The results are presented in fig. 21. One can observe that, whatever the adopted identification method, for too narrow filter bandwidth, the identified growth rate is far from the true one, whereas it converges to $\nu_{\text{t}}$ for large enough windows. On the other hand, when the filter is too wide, the effect of the neighboring mode starts to bias the identification. Therefore, when one analyses experimental data around a frequency of interest, there exist, for the filter bandwidth: i) a lower limit, given by the need not to alter the amplitude statistics, ii) an upper limit, given by the distance from the neighboring peaks. These constraints have to be satisfied in parallel with the one regarding the validity of the white noise approximation (section 4.1). However, in most of the practical cases, neighboring peaks are close and the maximum filter bandwidth to satisfy condition ii) is narrow enough that the effect of noise color can be safely neglected.

On the other hand, it is clear that any identification attempt on a mode that is both highly unstable and very close to another mode will fail because the filter to adopt to isolate one mode dynamics would be so narrow that condition i) is not fulfilled. In this situation a two-mode model would be required for parameter identification. As already suggested, it is advisable to iterate the parameter identification varying the applied filter bandwidth: one can be confident on the result if a plateau is observed.

5 Conclusion

In this work, the effects of the color of a stochastic excitation driving a Van der Pol Oscillator has been investigated. First, an Ornstein-Uhlenbeck process has been considered as the driving source. Then, a noise model, mimicking the stochastic forcing exerted by turbulence in thermoacoustic systems, has been used. It has been shown that in both cases the envelope statistics is the same as the one obtained with a white noise forcing, provided that an equivalent effective noise intensity is considered. Then, the approximation of a colored noise by a white one has been assessed in the context of data analysis and parameter identification. The main conclusion is that one can reliably identify the linear growth rate and saturation constant by band-pass filtering the data around the oscillator eigenfrequency.

This result is valid regardless of the parameters values and nature of the forcing noise. This fact consolidates the output-only parameter identification methods proposed in Noiray and Schuermans (2013a), because in real cases it might be impossible to determine the spectral distribution of the forcing noise.

Acknowledgement

This research is supported by the Swiss National Science Foundation under Grant 160579.