- Email: [email protected]
IDENTIFIABILITY ANALYSIS FOR NONLINEAR FEEDBACK LIMIT-CYCLE SYSTEMS USING HARMONIC BALANCE METHODS
14th IFAC Symposium on System Identification, Newcastle, Australia, 2006
IDENTIFIABILITY ANALYSIS FOR NONLINEAR FEEDBACK LIMIT-CYCLE SYSTEMS USING HARMONIC BALANCE METHODS Sangho Ko, Robert R. Bitmead
Department of Mechanical and Aerospace Engineering, University of California, San Diego, 9500 Gilman Drive, La Jolla, California, 92093-0411, USA. {sko,rbitmead}@ucsd.edu
Abstract: This paper proposes a method for analyzing identifiability of inherently closed-loop systems with nonlinear feedback by determining the order of persistent excitation of the input to the linear part of the dynamics. The method is based on the harmonic balance and the double Fourier series which is used for approximately expressing the output of the nonlinearity in the feedback loop. It is shown that a multiplicative periodic reference signal enhances the identifiability of the system by the additional sinusoids with interweaving frequencies to the original unforced c case. Copyright 2006 IFAC Keywords: Identifiability, Harmonic balance, Fourier analysis, Limit cycles, Closed-loop identification
1. INTRODUCTION
lean operating combustor typically shows high frequency pressure oscillations, mainly due to the nonlinear coupling between heat release and acoustics, which is referred to as thermoacoustic instability.
This study aims to develop a hands-on method of identifiability analysis for limit-cycle systems having a nonlinear feedback element as shown in Figure 1 where a periodic reference signal r(t) = 1 + cos ωr t is multiplied to the output of the nonlinearity Ψ(·). In this paper, we study if the input to the linear transfer function G(s) is persistently exciting to estimate the unknown parameters in G(s) with = 0 (r(t) ≡ 1, the unforced case), and see the effect of the periodic reference signal cos ωr t on the identifiability in terms of the order of persistent excitation.
Fig. 1. Multiplicative reference signal model
One example of such systems having a nonlinear feedback element, which has drawn much attention from aerospace academia and industries, is the combustor dynamics of gas turbine engines operating at low fuel-to-air ratios for reduced NOx formation (Annaswamy and Ghoniem, 2002). This
To cope with the various adverse effects on engine performance (e.g., reduced structural durability and engine blowout) arising from these very high frequency pressure oscillations, there have been extensive efforts to understand and model the underlying dynamics to actively control this
r(t)
u(t)
G(s)
y(t)
Ψ(·)
1341
phenomenon (Annaswamy and Ghoniem, 2002). The most recent and a promising identification method for the combustor dynamics is the socalled Grey-Box approach employed by Peracchio and Proscia (1999), Murray et al. (1998), Dunstan (2003), and McManus et al. (2004), considering its very complex dynamics which cannot be captured solely from physical laws. The model used in Dunstan (2003) and Murray et al. (1998) for this approach is shown in Fig. 2 where the part within the dotted line corresponds to G(s) of Fig. 1. The identification task for this model is now to estimate the parameters N, M, ζi , ωi , τ , and the nonlinearity Ψ(·) by using the measurement of the acoustic pressure p(t) (and/or the heat-release rate q(t)).
N s2+2ζ1ω1s+ω12
+ +
M s2+2ζ2ω2s+ω22 d dt
the identifiability of the system in nonlinear limitcycle operation by using the Harmonic Balance method (Mees, 1981). For the case of nonlinear limit-cycle operation, the questions we ask in this study are (Q1) Can we answer whether the input u(t) is persistently exciting for identifying G(s) for the unforced case r(t) = 1, by using only the pressure measurement {p(t)}? (Q2) What is the effect of the multiplicative periodic forcing (r(t) = 1 + cos ωr t) on identifiability for G(s) ?
p(t)
2. HARMONIC BALANCING AND DOUBLE FOURIER SERIES
d dt
When the unforced system in Fig. 1 operates in a limit cycle, the signals in the closed-loop must necessarily be periodic and, in many cases, the frequencies contained in those signals may not be harmonically related. Then, for such general cases, the assumed input to G(s) ∞ ∞ αi sin(ωi t + φi ), with |αi |2 < ∞. u(t) =
e−sτ q(t)
Ψ(·)
i=1
i=0
(1) must be equal to the output of the nonlinearity given by u(t) = Ψ(y(t)) (2) where ∞ αi G(jωi ) sin ωi t + φi + ∠G(jωi ) . y(t) =
Fig. 2. Combustion instability model For the unforced case (r(t) = 1), according to the experiment data in Dunstan (2003) and McManus et al. (2004), the pressure p(t) shows strong tonal signal consisting of a couple of harmonics. This situation may pose an information-lacking identification problem from the viewpoint of system identifiability in that the input u(t) to G(s) may not be persistently exciting. In general, for identifying Nu unknowns, we need a persistently exciting input of order Nu .
i=1
(3) The unknown coefficients αi s and other system parameters can be determined by using (1) and (2). This is called the harmonic balance method. Within authors’ knowledge, a computationally tractable method for expressing the output of the nonlinearity Ψ(y(t)) when y(t) is given by (3) does not exist. However, when y(t) is composed of two additive sinusoids, which corresponds to the current combustion instability problem, we can apply the double Fourier series (Gelb and Vander Velde, 1968) to express the output. That is, for a sum of two sinusoidal inputs
In order to see the effect of fuel modulation on identifiability of the combustor model in lean conditions, Dunstan and Bitmead (2003) and Dunstan (2003) used a multiplicative reference signal model of Fig. 1 for the case of linear regime of operation by replacing the saturation nonlinearity Ψ(·) by a static gain. For a periodic reference signal r(t) = 1 + cos ωr t (with a small ), they used a lifting technique (Khargonekar et al., 1985), changing the linear periodic system into a linear time-invariant system, and concluded that the presence of a multiplicative periodic reference signal can engender identifiability by introducing sidebands into the original baseband system with r(t) = 1.
y(t) = A sin ψ1 + B sin ψ2 ,
(4)
the nonlinearity output can be expressed by u(t) = Ψ(A sin ψ1 + B sin ψ2 ) ∞ ∞ Pmn sin(mψ1 + nψ2 )+ = m=0 n = −∞(m = 0) n = 0(m = 0)
Qmn cos(mψ1 + nψ2 ) .
As a continuing effort to the previous study of Dunstan and Bitmead (2003), this paper studies
(5)
1342
The Fourier coefficients Pmn and Qmn are given by π 1 Ψ y(t) sin(mψ1 +nψ2 )dψ1 dψ2 Pmn = 2 2π −π π 1 Qmn = 2 Ψ y(t) cos(mψ1 +nψ2 )dψ1 dψ2 2π −π (6)
Remark 1. Since we consider a system in a nonlinear limit-cycle, if the input range to the nonlinearity is [−ym , ym ] (from the measurement), then for a fixed (but unknown) d1 , the unknown parameter d2 of the nonlinearity N (·) in Fig 3 should be approximately in the range d2 ∈ { d2 | N (−ym ) ≥ 0.5d1 }, because outside of this range the system would be almost in a linear operation.
For the case of static single-valued non-linearities, it can be shown that
Remark 2. Since we are interested in only the frequency content in u(t) = N (y) + b and the parameter d1 scales only the magnitude of u(t), we can fix d1 arbitrarily.
• Qmn = 0 for (m + n) odd • Pmn = 0 for (m + n) even • In addition, if the non-linear function Ψ(·) is odd, all Qmn = 0.
3.1 Unforced Case: r(t) = 1 3. APPLICATION TO COMBUSTION INSTABILITY
From the assumption (A1), we further approximate the output y(t) of G(s) as a sum of two sinusoidal signals
In order to apply the harmonic balance method to the current combustion instability problem and to answer to the questions raised at the end of Section 1, we assume the followings:
y(t)
2
αi G(jωi ) sin ωi t + φi + ∠G(jωi ) ,
i=1
(8) when the input to G(s) is (1). By applying the double Fourier series, we obtain u(t) = Ψ y(t) ∞ ∞ Pmn sin(mψ1 + nψ2 ). (9) =
(A1) The linear transfer function G(s) has two dominant harmonics (as shown in Fig. 2) such that • |G(jω)| 0, for ω ∈ / [ω1 ±1 ]∪[ω2 ±2 ] • ∠G(jω) −π, through [ωi ±i ], i = 1, 2. i is a very small number compared to ωi .
m=0
(A2) The nonlinear function Ψ(·) is of the form Ψ(·) = N (·) + b, where N (·) is a nonlinear odd function and b is a constant. Furthermore, N (·) is of a negative saturation type (Fig. 3) which can be represented by 2d1 · tan−1 (y/d2 ) π iN a2i+1 y 2i+1
where, for iN = 3, the Fourier coefficients Pmn can be expressed as Pmn =
N (y) = −
(7)
J(i, m, n)=
(10)
2i+1 2(i+1)−i m
(i) ¯ δ(m−im )δ(n−i n )Cim in .
im =0 in =0 im +in =odd
where d1 , d2 > 0 are unknown and iN < ∞ is a positive integer.
(11) By comparing (9) with (1), we can express the assumed amplitude αi as a function of Pmn for each ωi . In (11), δ(x) represents the Kronecker delta (e.g., δ(x) = 1 for x = 0, δ(x) = 0 otherwise), and ¯ − in ) δ(n − in ) + (−1)in · δ(n + in ). δ(n
N(y) N(−ym) = 0.9 d1 N(−ym) = 0.7 d1 N(−ym) = 0.5 d1 N(−ym) = 0.3 d1 ym
d2 increases
iN 1 a2i+1 J(i, m, n) 2π 2 i=0
with
i=0
d1
n = −∞ n+m = odd
(i)
Cim in are dependent on the coefficients Ai αi |G(jωi )| for i = 1, 2 as shown in Appendix.
y −ym
For (11), the following formula was used π sin pψ1 cos qψ2 sin(mψ1 + nψ2 )dψ1 dψ2
N(y) = − d1(2/π)tan−1(y/d2)
−π
−d1
= π 2 [δ(n − q) + δ(q + n)] × [δ(m − p) − δ(m + p)] δ c (p)δ c (m) Fig. 3. Nonlinearity model
with δ c (x) 1 − δ(x).
1343
Because of a finite number of the Kronecker delta functions in J(i, m, n) in (11) (in turn, in Pmn ), the input u(t) given by (9) actually consists of only a finite number of sinusoidal signals. Practically, when counting the number of different frequencies contained in u(t), any frequency components higher than an interesting frequency range (usually set by the sampling frequency) and the sinusoidal components having small magnitude should be ignored.
p(t) which can be fitted into a sum of two sinusoid using the measured data, and suppose that the fitted model is p(t) = f1 · sin ω1 t + f2 · cos ω2 t
(15)
with f1 = −0.6062 × 10−3 , f2 = −0.1584 × 10−3 , ω1 = 170(Hz), and ω2 = 2 · ω1 = 340(Hz). Then, the input to the nonlinearity is y(t) = p(t ˙ − τ)=
2
Ai sin ωi t + φi + ∠G(jωi )
i=1
(16) 3.2 The Effect of Multiplicative Reference Signal
where • A1 = f1 · ω1 = 0.6595, A2 = f2 · ω2 = 0.3447 • φ1 = −τ ω1 + 0.5π, φ2 = 2φ1 • ∠G(jω1 ) = ∠G(jω2 ) = −π.
Now, we see the effect of the periodic forcing. Using r(t) = 1 + cos(ωr t + θr ) with nonzero changes the input u(t) to um (t): um (t) = r(t) · u(t) = u(t) + u (t)
The following parameters of the multiplicative reference signal are chosen
(12)
• = 0.2, ωr = 80 (Hz)
where r(t) = 1 + cos(ωr t + θr ) u (t) = u(t) cos(ωr t + θr ).
(13)
to generate side-band peaks between the resonant frequencies by a subharmonic forcing.
With ψi ωi t + φi + ∠G(jωi ) for i = 1, 2 and ψ3 ωr t + θr , we have ∞ ∞ Pmn sin(mψ1 +nψ2 −ψ3 ) u (t)= 2 m=0
By substituting these parameters into (9), (10), and (11), we obtain Fig. 4 and 5 which show the magnitude of the sinusoids contained in u(t) for an arbitrarily fixed d1 = 1 and the different values of d2 ( d2 = 0.1584 for N (−ym ) = 0.9d1 in Fig. 4 and d2 = 1.0 for N (−ym ) = 0.5d1 in Fig. 5).
n=−∞ n+m=odd
+ sin(mψ1 +nψ2 +ψ3 ) .
N(ym) = −0.9 d1
1.0
(14) Compared to the input u(t) given in (9), the new input um (t) contains the additional sinusoidal signals with interweaving frequencies 1 mω1 + nω2 ± ωr to the original frequencies mω1 + nω2 . Therefore, it can be said that the multiplicative periodic forcing enhances the identifiability of the system. This will be demonstrated by using a numerical example in the next section.
without periodic ref signal peridodic ref signal effect
Amplitude
0.8
0.6
0.4
0.2
Selecting the optimal values of and ωr for enhancing the identifiability should be compromised with the following limitations (Dunstan, 2003):
0
0
500
1000
1500
Frequency (Hz)
Fig. 4. Spectral contents of u(t): N (−ym ) = 0.9d1
• needs to be small for the small flow variation due to the fuel modulation
N(ym) = −0.5 d1
0.4
without periodic ref signal periodic ref signal effect
• ωr less than the first harmonic ω1 is preferred (subharmonic forcing) due to the physically limited bandwidth of pulsating valves for fuel modulation
Amplitude
0.3
0.2
4. NUMERICAL EXAMPLE
0.1
For the combustion instability model in Fig. 2, the only available measurement data is the pressure
0
0
500
1000
1500
Frequency (Hz)
Fig. 5. Spectral contents of u(t): N (−ym ) = 0.5d1
1 This result is similar to the side-band spectral peaks at both ω1 + ωr and ω1 − ωr of the experimental data in McManus et al. (2004).
From these, we observe the following facts:
1344
• The high frequency components ωi = i · ω1 for i ≥ 3 are caused by the nonlinearity N (·), and the smaller d2 , the bigger the nonlinearity effect, as shown Fig. 3. Therefore, in a sense, the nonlinearity enhances the identifiability by adding the high frequency sinusoids.
REFERENCES Annaswamy, Anuradha M. and Ahmed F. Ghoniem (2002). Active control of combustion instability: theory and practice. IEEE Control Systems Magazine 22(6), 37–54. Dunstan, Wayne J. (2003). System identification of nonlinear resonant systems. PhD thesis. University of California, San Diego. La Jolla, CA, USA. Dunstan, Wayne J. and Robert R. Bitmead (2003). Identification of resonant systems using periodic multiplicative reference signals. In: 13th IFAC Symposium on System Identification. Rotterdam, Netherlands. pp. 671–676. Gelb, Arthur and Wallace E. Vander Velde (1968). Multiple-input describing functions and nonlinear system design. McGraw-Hill, Inc. New York, NY, USA. Khargonekar, Pramod P., Kameshwar Poolla and Allen Tannenbaum (1985). Robust control of linear time-invariant plants using periodic compensation. IEEE Transactions on Automatic Conrol 30(11), 1088–1096. McManus, Keith, Fei Han, Wayne Dunstan, Corneliu Barbu and Minesh Shah (2004). Modeling and control of combustion dynamics in industrial gas turbines. In: Proceedings of ASME Turbo Expo 2004. Vienna, Austria. pp. 567–575. Mees, A. I. (1981). Dynamics of feedback systems. John-Wiley & Sons. New York, USA. Mees, Alistair I. and Leon O. Chua (1979). The Hopf bifurcation theorem and its applications to nonlinear oscillations in circuits and systems. IEEE Transactions on Circuits and Systems CAS-26(4), 235–254. Moiola, Jorge L., H´ector G. Chiacchiarini and Alfredo C. Desages (1996). Bifurcations and Hopf degeneracies in nonlinear feedback systems with time delay. International Journal of Bifurcation and Chaos 6(4), 661–672. Murray, R. M., C. A. Jacobson, R. Casa, A. I. Khibnik, C. R. Johnson Jr., R. Bitmead, A. A. Peracchio and W. M. Proscia (1998). System identification for limit cycling systems: a case study for combustion instabilities. In: Proceeding of the American Control Conference. Pittsburgh, PA, USA. pp. 2004– 2008. Peracchio, A. A. and W. M. Proscia (1999). Nonlinear heat-release/acoustic model for thermoacoustic instability in lean premixed combustors. Journal of Engineering for Gas Turbines and Power 121(3), 415–421. Thothadri, M., R. A. Casas, F. C. Moon, R. D’andrea and C. R. Johnson Jr. (2003). Nonlinear system identification of mulidegree-of-freedom systems. Nonlinear Dynamics 32(3), 307–322.
• As d2 varies, the frequency components in u(t) do not change noticeably, except for the high frequency range. Therefore, the order of persistent excitation is weakly dependent on the parameter d2 for the low frequency range. • If an interesting frequency band is the range less than 400 Hz and we ignore some components having small amplitude, then the order of persistent excitation for both Fig. 4 and Fig. 5 is approximately 5 (one constant and two sinusoids), with r(t) = 1. Therefore, from Fig. 2, it can be said that for the 5 unknowns (ω1 , ζ1 = ζ2 , M , N , and τ ), the input is persistently exciting. • The additional periodic reference signal increases the order of persistent excitation by 4, with the additional 2 interweaving sinusoids. In this case, the linear model G(s) with two different damping ratios ζ1 = ζ2 can be identified.
5. CONCLUDING REMARKS In this paper, a method for analyzing identifiability was introduced for nonlinear feedback limitcycle systems by using the harmonic balance. With the double Fourier series, the approximate expression of the input to the linear part of the dynamics was derived as a sum of sinusoids with different frequencies and magnitudes. From this, we can determine the order of persistent excitation by counting the number of different frequencies of the sinusoids contained in the input to the linear part of the dynamics. Also, the enhanced identifiability through a multiplicative periodic reference signal was confirmed by a numerical example. The following tasks are remained as future works: 1. Developing an identification method for identifying the combustion dynamics modelled by a nonlinear feedback limit-cycle system, by using only the pressure measurement 2. Application of the existing identification methods in Van Pelt and Bernstein (2001) and Thothadri et al. (2003) for nonlinear feedback systems 3. Connection with the Hopf bifurcation (Mees and Chua (1979) and Moiola et al. (1996))
1345
9A4 A3 35 A72 +3A21 A52 + 1 2 +A61 A2 16 4 2
7 2 5 105A 105A41 A32 21A A (3) 2 1 2 2 − − C03 =2π − 64 32 32
7 2 5 21A 7A A (3) 2 1 2 + C05 =2π 2 64 32 7
A (3) 2 C07 =2π − 2 64
Van Pelt, Tobin H. and Dennis S. Bernstein (2001). Non-linear system identification using Hammerstein and non-linear feedback models with piecewise linear static maps. International Journal of Control 74(18), 1807–1823.
(3)
C01 =2π 2
APPENDIX
9A3 A4 35 A71 +3A51 A22 + 1 2 +A1 A62 16 4 2
5 2 6 105A 105A31 A42 105A A A 1 2 (3) 1 2 2 − − C12 =π − 16 32 8
6 3 4 105A 21A A A 1 (3) 2 1 2 + C14 =π 2 16 32
7A1 A62 (3) 2 C16 =π − 32
105A61 A2 105A41 A32 −105A21 A52 (3) − − C21 =π 2 16 32 8
2 5 4 3 35A1 A2 105A1 A2 (3) + C23 =π 2 32 8
2 5 21A A (3) 1 2 C25 =π 2 − 32
(0)
(3)
C01 =2π 2 A2
C10 =π 2
(0)
C10 =π 2 A1
3A32 3A21 A2 + 2 4
3 A (1) C03 =2π 2 − 2 4
3 3A1 A22 3A1 (1) + C10 =π 2 4 2
2 3A1 A2 (1) C12 =π 2 − 2
2 3A A2 (1) C21 =π 2 − 1 2
3 A (1) 1 2 C30 =π − 4 (1) C01 =2π 2
30A21 A32 15A41 A2 5A52 + + 8 8 8
5 2 3 10A 5A A (2) 1 2 C03 =2π 2 − 2 − 16 8 5
A2 (2) 2 C05 =2π 16
5 30A31 A22 15A1 A42 5A1 (2) 2 + + C10 =π 8 8 8
3 2 4 5A1 A2 30A1 A2 (2) − C12 =π 2 − 8 2
4 5A1 A2 (2) C14 =π 2 8
2 3 5A41 A2 30A (2) 1 A2 − C21 =π 2 − 8 2
2 3 10A A (2) 1 2 C23 =π 2 8
10A31 A22 5A51 (2) 2 − C30 =π − 16 8
3 2 10A1 A2 (2) C32 =π 2 8 4
5A1 A2 (2) C41 =π 2 8 5
A (2) 1 C50 =π 2 16 (2)
C01 =2π 2
105A51 A22 105A31 A42 −21A71 − − 64 32 32
5 2 3 4 35A 105A A A (3) 1 2 1 2 + C32 =π 2 32 8
35A31 A42 (3) 2 C34 =π − 32
6 105A41 A32 21A (3) 1 A2 2 + C41 =π 16 32
4 3 35A1 A2 (3) C43 =π 2 − 32
7 21A51 A22 7A1 (3) + C50 =π 2 64 32
5 2 21A1 A2 (3) C52 =π 2 − 32
6 7A A2 (3) C61 =π 2 − 1 32
7 A (3) 1 2 C70 =π − 64 (3)
C30 =π 2
1346
IDENTIFIABILITY ANALYSIS FOR NONLINEAR FEEDBACK LIMIT-CYCLE SYSTEMS USING HARMONIC BALANCE METHODS Sangho Ko, Robert R. Bitmead
Department of Mechanical and Aerospace Engineering, University of California, San Diego, 9500 Gilman Drive, La Jolla, California, 92093-0411, USA. {sko,rbitmead}@ucsd.edu
Abstract: This paper proposes a method for analyzing identifiability of inherently closed-loop systems with nonlinear feedback by determining the order of persistent excitation of the input to the linear part of the dynamics. The method is based on the harmonic balance and the double Fourier series which is used for approximately expressing the output of the nonlinearity in the feedback loop. It is shown that a multiplicative periodic reference signal enhances the identifiability of the system by the additional sinusoids with interweaving frequencies to the original unforced c case. Copyright 2006 IFAC Keywords: Identifiability, Harmonic balance, Fourier analysis, Limit cycles, Closed-loop identification
1. INTRODUCTION
lean operating combustor typically shows high frequency pressure oscillations, mainly due to the nonlinear coupling between heat release and acoustics, which is referred to as thermoacoustic instability.
This study aims to develop a hands-on method of identifiability analysis for limit-cycle systems having a nonlinear feedback element as shown in Figure 1 where a periodic reference signal r(t) = 1 + cos ωr t is multiplied to the output of the nonlinearity Ψ(·). In this paper, we study if the input to the linear transfer function G(s) is persistently exciting to estimate the unknown parameters in G(s) with = 0 (r(t) ≡ 1, the unforced case), and see the effect of the periodic reference signal cos ωr t on the identifiability in terms of the order of persistent excitation.
Fig. 1. Multiplicative reference signal model
One example of such systems having a nonlinear feedback element, which has drawn much attention from aerospace academia and industries, is the combustor dynamics of gas turbine engines operating at low fuel-to-air ratios for reduced NOx formation (Annaswamy and Ghoniem, 2002). This
To cope with the various adverse effects on engine performance (e.g., reduced structural durability and engine blowout) arising from these very high frequency pressure oscillations, there have been extensive efforts to understand and model the underlying dynamics to actively control this
r(t)
u(t)
G(s)
y(t)
Ψ(·)
1341
phenomenon (Annaswamy and Ghoniem, 2002). The most recent and a promising identification method for the combustor dynamics is the socalled Grey-Box approach employed by Peracchio and Proscia (1999), Murray et al. (1998), Dunstan (2003), and McManus et al. (2004), considering its very complex dynamics which cannot be captured solely from physical laws. The model used in Dunstan (2003) and Murray et al. (1998) for this approach is shown in Fig. 2 where the part within the dotted line corresponds to G(s) of Fig. 1. The identification task for this model is now to estimate the parameters N, M, ζi , ωi , τ , and the nonlinearity Ψ(·) by using the measurement of the acoustic pressure p(t) (and/or the heat-release rate q(t)).
N s2+2ζ1ω1s+ω12
+ +
M s2+2ζ2ω2s+ω22 d dt
the identifiability of the system in nonlinear limitcycle operation by using the Harmonic Balance method (Mees, 1981). For the case of nonlinear limit-cycle operation, the questions we ask in this study are (Q1) Can we answer whether the input u(t) is persistently exciting for identifying G(s) for the unforced case r(t) = 1, by using only the pressure measurement {p(t)}? (Q2) What is the effect of the multiplicative periodic forcing (r(t) = 1 + cos ωr t) on identifiability for G(s) ?
p(t)
2. HARMONIC BALANCING AND DOUBLE FOURIER SERIES
d dt
When the unforced system in Fig. 1 operates in a limit cycle, the signals in the closed-loop must necessarily be periodic and, in many cases, the frequencies contained in those signals may not be harmonically related. Then, for such general cases, the assumed input to G(s) ∞ ∞ αi sin(ωi t + φi ), with |αi |2 < ∞. u(t) =
e−sτ q(t)
Ψ(·)
i=1
i=0
(1) must be equal to the output of the nonlinearity given by u(t) = Ψ(y(t)) (2) where ∞ αi G(jωi ) sin ωi t + φi + ∠G(jωi ) . y(t) =
Fig. 2. Combustion instability model For the unforced case (r(t) = 1), according to the experiment data in Dunstan (2003) and McManus et al. (2004), the pressure p(t) shows strong tonal signal consisting of a couple of harmonics. This situation may pose an information-lacking identification problem from the viewpoint of system identifiability in that the input u(t) to G(s) may not be persistently exciting. In general, for identifying Nu unknowns, we need a persistently exciting input of order Nu .
i=1
(3) The unknown coefficients αi s and other system parameters can be determined by using (1) and (2). This is called the harmonic balance method. Within authors’ knowledge, a computationally tractable method for expressing the output of the nonlinearity Ψ(y(t)) when y(t) is given by (3) does not exist. However, when y(t) is composed of two additive sinusoids, which corresponds to the current combustion instability problem, we can apply the double Fourier series (Gelb and Vander Velde, 1968) to express the output. That is, for a sum of two sinusoidal inputs
In order to see the effect of fuel modulation on identifiability of the combustor model in lean conditions, Dunstan and Bitmead (2003) and Dunstan (2003) used a multiplicative reference signal model of Fig. 1 for the case of linear regime of operation by replacing the saturation nonlinearity Ψ(·) by a static gain. For a periodic reference signal r(t) = 1 + cos ωr t (with a small ), they used a lifting technique (Khargonekar et al., 1985), changing the linear periodic system into a linear time-invariant system, and concluded that the presence of a multiplicative periodic reference signal can engender identifiability by introducing sidebands into the original baseband system with r(t) = 1.
y(t) = A sin ψ1 + B sin ψ2 ,
(4)
the nonlinearity output can be expressed by u(t) = Ψ(A sin ψ1 + B sin ψ2 ) ∞ ∞ Pmn sin(mψ1 + nψ2 )+ = m=0 n = −∞(m = 0) n = 0(m = 0)
Qmn cos(mψ1 + nψ2 ) .
As a continuing effort to the previous study of Dunstan and Bitmead (2003), this paper studies
(5)
1342
The Fourier coefficients Pmn and Qmn are given by π 1 Ψ y(t) sin(mψ1 +nψ2 )dψ1 dψ2 Pmn = 2 2π −π π 1 Qmn = 2 Ψ y(t) cos(mψ1 +nψ2 )dψ1 dψ2 2π −π (6)
Remark 1. Since we consider a system in a nonlinear limit-cycle, if the input range to the nonlinearity is [−ym , ym ] (from the measurement), then for a fixed (but unknown) d1 , the unknown parameter d2 of the nonlinearity N (·) in Fig 3 should be approximately in the range d2 ∈ { d2 | N (−ym ) ≥ 0.5d1 }, because outside of this range the system would be almost in a linear operation.
For the case of static single-valued non-linearities, it can be shown that
Remark 2. Since we are interested in only the frequency content in u(t) = N (y) + b and the parameter d1 scales only the magnitude of u(t), we can fix d1 arbitrarily.
• Qmn = 0 for (m + n) odd • Pmn = 0 for (m + n) even • In addition, if the non-linear function Ψ(·) is odd, all Qmn = 0.
3.1 Unforced Case: r(t) = 1 3. APPLICATION TO COMBUSTION INSTABILITY
From the assumption (A1), we further approximate the output y(t) of G(s) as a sum of two sinusoidal signals
In order to apply the harmonic balance method to the current combustion instability problem and to answer to the questions raised at the end of Section 1, we assume the followings:
y(t)
2
αi G(jωi ) sin ωi t + φi + ∠G(jωi ) ,
i=1
(8) when the input to G(s) is (1). By applying the double Fourier series, we obtain u(t) = Ψ y(t) ∞ ∞ Pmn sin(mψ1 + nψ2 ). (9) =
(A1) The linear transfer function G(s) has two dominant harmonics (as shown in Fig. 2) such that • |G(jω)| 0, for ω ∈ / [ω1 ±1 ]∪[ω2 ±2 ] • ∠G(jω) −π, through [ωi ±i ], i = 1, 2. i is a very small number compared to ωi .
m=0
(A2) The nonlinear function Ψ(·) is of the form Ψ(·) = N (·) + b, where N (·) is a nonlinear odd function and b is a constant. Furthermore, N (·) is of a negative saturation type (Fig. 3) which can be represented by 2d1 · tan−1 (y/d2 ) π iN a2i+1 y 2i+1
where, for iN = 3, the Fourier coefficients Pmn can be expressed as Pmn =
N (y) = −
(7)
J(i, m, n)=
(10)
2i+1 2(i+1)−i m
(i) ¯ δ(m−im )δ(n−i n )Cim in .
im =0 in =0 im +in =odd
where d1 , d2 > 0 are unknown and iN < ∞ is a positive integer.
(11) By comparing (9) with (1), we can express the assumed amplitude αi as a function of Pmn for each ωi . In (11), δ(x) represents the Kronecker delta (e.g., δ(x) = 1 for x = 0, δ(x) = 0 otherwise), and ¯ − in ) δ(n − in ) + (−1)in · δ(n + in ). δ(n
N(y) N(−ym) = 0.9 d1 N(−ym) = 0.7 d1 N(−ym) = 0.5 d1 N(−ym) = 0.3 d1 ym
d2 increases
iN 1 a2i+1 J(i, m, n) 2π 2 i=0
with
i=0
d1
n = −∞ n+m = odd
(i)
Cim in are dependent on the coefficients Ai αi |G(jωi )| for i = 1, 2 as shown in Appendix.
y −ym
For (11), the following formula was used π sin pψ1 cos qψ2 sin(mψ1 + nψ2 )dψ1 dψ2
N(y) = − d1(2/π)tan−1(y/d2)
−π
−d1
= π 2 [δ(n − q) + δ(q + n)] × [δ(m − p) − δ(m + p)] δ c (p)δ c (m) Fig. 3. Nonlinearity model
with δ c (x) 1 − δ(x).
1343
Because of a finite number of the Kronecker delta functions in J(i, m, n) in (11) (in turn, in Pmn ), the input u(t) given by (9) actually consists of only a finite number of sinusoidal signals. Practically, when counting the number of different frequencies contained in u(t), any frequency components higher than an interesting frequency range (usually set by the sampling frequency) and the sinusoidal components having small magnitude should be ignored.
p(t) which can be fitted into a sum of two sinusoid using the measured data, and suppose that the fitted model is p(t) = f1 · sin ω1 t + f2 · cos ω2 t
(15)
with f1 = −0.6062 × 10−3 , f2 = −0.1584 × 10−3 , ω1 = 170(Hz), and ω2 = 2 · ω1 = 340(Hz). Then, the input to the nonlinearity is y(t) = p(t ˙ − τ)=
2
Ai sin ωi t + φi + ∠G(jωi )
i=1
(16) 3.2 The Effect of Multiplicative Reference Signal
where • A1 = f1 · ω1 = 0.6595, A2 = f2 · ω2 = 0.3447 • φ1 = −τ ω1 + 0.5π, φ2 = 2φ1 • ∠G(jω1 ) = ∠G(jω2 ) = −π.
Now, we see the effect of the periodic forcing. Using r(t) = 1 + cos(ωr t + θr ) with nonzero changes the input u(t) to um (t): um (t) = r(t) · u(t) = u(t) + u (t)
The following parameters of the multiplicative reference signal are chosen
(12)
• = 0.2, ωr = 80 (Hz)
where r(t) = 1 + cos(ωr t + θr ) u (t) = u(t) cos(ωr t + θr ).
(13)
to generate side-band peaks between the resonant frequencies by a subharmonic forcing.
With ψi ωi t + φi + ∠G(jωi ) for i = 1, 2 and ψ3 ωr t + θr , we have ∞ ∞ Pmn sin(mψ1 +nψ2 −ψ3 ) u (t)= 2 m=0
By substituting these parameters into (9), (10), and (11), we obtain Fig. 4 and 5 which show the magnitude of the sinusoids contained in u(t) for an arbitrarily fixed d1 = 1 and the different values of d2 ( d2 = 0.1584 for N (−ym ) = 0.9d1 in Fig. 4 and d2 = 1.0 for N (−ym ) = 0.5d1 in Fig. 5).
n=−∞ n+m=odd
+ sin(mψ1 +nψ2 +ψ3 ) .
N(ym) = −0.9 d1
1.0
(14) Compared to the input u(t) given in (9), the new input um (t) contains the additional sinusoidal signals with interweaving frequencies 1 mω1 + nω2 ± ωr to the original frequencies mω1 + nω2 . Therefore, it can be said that the multiplicative periodic forcing enhances the identifiability of the system. This will be demonstrated by using a numerical example in the next section.
without periodic ref signal peridodic ref signal effect
Amplitude
0.8
0.6
0.4
0.2
Selecting the optimal values of and ωr for enhancing the identifiability should be compromised with the following limitations (Dunstan, 2003):
0
0
500
1000
1500
Frequency (Hz)
Fig. 4. Spectral contents of u(t): N (−ym ) = 0.9d1
• needs to be small for the small flow variation due to the fuel modulation
N(ym) = −0.5 d1
0.4
without periodic ref signal periodic ref signal effect
• ωr less than the first harmonic ω1 is preferred (subharmonic forcing) due to the physically limited bandwidth of pulsating valves for fuel modulation
Amplitude
0.3
0.2
4. NUMERICAL EXAMPLE
0.1
For the combustion instability model in Fig. 2, the only available measurement data is the pressure
0
0
500
1000
1500
Frequency (Hz)
Fig. 5. Spectral contents of u(t): N (−ym ) = 0.5d1
1 This result is similar to the side-band spectral peaks at both ω1 + ωr and ω1 − ωr of the experimental data in McManus et al. (2004).
From these, we observe the following facts:
1344
• The high frequency components ωi = i · ω1 for i ≥ 3 are caused by the nonlinearity N (·), and the smaller d2 , the bigger the nonlinearity effect, as shown Fig. 3. Therefore, in a sense, the nonlinearity enhances the identifiability by adding the high frequency sinusoids.
REFERENCES Annaswamy, Anuradha M. and Ahmed F. Ghoniem (2002). Active control of combustion instability: theory and practice. IEEE Control Systems Magazine 22(6), 37–54. Dunstan, Wayne J. (2003). System identification of nonlinear resonant systems. PhD thesis. University of California, San Diego. La Jolla, CA, USA. Dunstan, Wayne J. and Robert R. Bitmead (2003). Identification of resonant systems using periodic multiplicative reference signals. In: 13th IFAC Symposium on System Identification. Rotterdam, Netherlands. pp. 671–676. Gelb, Arthur and Wallace E. Vander Velde (1968). Multiple-input describing functions and nonlinear system design. McGraw-Hill, Inc. New York, NY, USA. Khargonekar, Pramod P., Kameshwar Poolla and Allen Tannenbaum (1985). Robust control of linear time-invariant plants using periodic compensation. IEEE Transactions on Automatic Conrol 30(11), 1088–1096. McManus, Keith, Fei Han, Wayne Dunstan, Corneliu Barbu and Minesh Shah (2004). Modeling and control of combustion dynamics in industrial gas turbines. In: Proceedings of ASME Turbo Expo 2004. Vienna, Austria. pp. 567–575. Mees, A. I. (1981). Dynamics of feedback systems. John-Wiley & Sons. New York, USA. Mees, Alistair I. and Leon O. Chua (1979). The Hopf bifurcation theorem and its applications to nonlinear oscillations in circuits and systems. IEEE Transactions on Circuits and Systems CAS-26(4), 235–254. Moiola, Jorge L., H´ector G. Chiacchiarini and Alfredo C. Desages (1996). Bifurcations and Hopf degeneracies in nonlinear feedback systems with time delay. International Journal of Bifurcation and Chaos 6(4), 661–672. Murray, R. M., C. A. Jacobson, R. Casa, A. I. Khibnik, C. R. Johnson Jr., R. Bitmead, A. A. Peracchio and W. M. Proscia (1998). System identification for limit cycling systems: a case study for combustion instabilities. In: Proceeding of the American Control Conference. Pittsburgh, PA, USA. pp. 2004– 2008. Peracchio, A. A. and W. M. Proscia (1999). Nonlinear heat-release/acoustic model for thermoacoustic instability in lean premixed combustors. Journal of Engineering for Gas Turbines and Power 121(3), 415–421. Thothadri, M., R. A. Casas, F. C. Moon, R. D’andrea and C. R. Johnson Jr. (2003). Nonlinear system identification of mulidegree-of-freedom systems. Nonlinear Dynamics 32(3), 307–322.
• As d2 varies, the frequency components in u(t) do not change noticeably, except for the high frequency range. Therefore, the order of persistent excitation is weakly dependent on the parameter d2 for the low frequency range. • If an interesting frequency band is the range less than 400 Hz and we ignore some components having small amplitude, then the order of persistent excitation for both Fig. 4 and Fig. 5 is approximately 5 (one constant and two sinusoids), with r(t) = 1. Therefore, from Fig. 2, it can be said that for the 5 unknowns (ω1 , ζ1 = ζ2 , M , N , and τ ), the input is persistently exciting. • The additional periodic reference signal increases the order of persistent excitation by 4, with the additional 2 interweaving sinusoids. In this case, the linear model G(s) with two different damping ratios ζ1 = ζ2 can be identified.
5. CONCLUDING REMARKS In this paper, a method for analyzing identifiability was introduced for nonlinear feedback limitcycle systems by using the harmonic balance. With the double Fourier series, the approximate expression of the input to the linear part of the dynamics was derived as a sum of sinusoids with different frequencies and magnitudes. From this, we can determine the order of persistent excitation by counting the number of different frequencies of the sinusoids contained in the input to the linear part of the dynamics. Also, the enhanced identifiability through a multiplicative periodic reference signal was confirmed by a numerical example. The following tasks are remained as future works: 1. Developing an identification method for identifying the combustion dynamics modelled by a nonlinear feedback limit-cycle system, by using only the pressure measurement 2. Application of the existing identification methods in Van Pelt and Bernstein (2001) and Thothadri et al. (2003) for nonlinear feedback systems 3. Connection with the Hopf bifurcation (Mees and Chua (1979) and Moiola et al. (1996))
1345
9A4 A3 35 A72 +3A21 A52 + 1 2 +A61 A2 16 4 2
7 2 5 105A 105A41 A32 21A A (3) 2 1 2 2 − − C03 =2π − 64 32 32
7 2 5 21A 7A A (3) 2 1 2 + C05 =2π 2 64 32 7
A (3) 2 C07 =2π − 2 64
Van Pelt, Tobin H. and Dennis S. Bernstein (2001). Non-linear system identification using Hammerstein and non-linear feedback models with piecewise linear static maps. International Journal of Control 74(18), 1807–1823.
(3)
C01 =2π 2
APPENDIX
9A3 A4 35 A71 +3A51 A22 + 1 2 +A1 A62 16 4 2
5 2 6 105A 105A31 A42 105A A A 1 2 (3) 1 2 2 − − C12 =π − 16 32 8
6 3 4 105A 21A A A 1 (3) 2 1 2 + C14 =π 2 16 32
7A1 A62 (3) 2 C16 =π − 32
105A61 A2 105A41 A32 −105A21 A52 (3) − − C21 =π 2 16 32 8
2 5 4 3 35A1 A2 105A1 A2 (3) + C23 =π 2 32 8
2 5 21A A (3) 1 2 C25 =π 2 − 32
(0)
(3)
C01 =2π 2 A2
C10 =π 2
(0)
C10 =π 2 A1
3A32 3A21 A2 + 2 4
3 A (1) C03 =2π 2 − 2 4
3 3A1 A22 3A1 (1) + C10 =π 2 4 2
2 3A1 A2 (1) C12 =π 2 − 2
2 3A A2 (1) C21 =π 2 − 1 2
3 A (1) 1 2 C30 =π − 4 (1) C01 =2π 2
30A21 A32 15A41 A2 5A52 + + 8 8 8
5 2 3 10A 5A A (2) 1 2 C03 =2π 2 − 2 − 16 8 5
A2 (2) 2 C05 =2π 16
5 30A31 A22 15A1 A42 5A1 (2) 2 + + C10 =π 8 8 8
3 2 4 5A1 A2 30A1 A2 (2) − C12 =π 2 − 8 2
4 5A1 A2 (2) C14 =π 2 8
2 3 5A41 A2 30A (2) 1 A2 − C21 =π 2 − 8 2
2 3 10A A (2) 1 2 C23 =π 2 8
10A31 A22 5A51 (2) 2 − C30 =π − 16 8
3 2 10A1 A2 (2) C32 =π 2 8 4
5A1 A2 (2) C41 =π 2 8 5
A (2) 1 C50 =π 2 16 (2)
C01 =2π 2
105A51 A22 105A31 A42 −21A71 − − 64 32 32
5 2 3 4 35A 105A A A (3) 1 2 1 2 + C32 =π 2 32 8
35A31 A42 (3) 2 C34 =π − 32
6 105A41 A32 21A (3) 1 A2 2 + C41 =π 16 32
4 3 35A1 A2 (3) C43 =π 2 − 32
7 21A51 A22 7A1 (3) + C50 =π 2 64 32
5 2 21A1 A2 (3) C52 =π 2 − 32
6 7A A2 (3) C61 =π 2 − 1 32
7 A (3) 1 2 C70 =π − 64 (3)
C30 =π 2
1346