• Author / Uploaded
• sumatrablackcoffee453

Citation preview

ARTICLE IN PRESS

Mechanical Systems and Signal Processing Mechanical Systems and Signal Processing 18 (2004) 431–455 www.elsevier.com/locate/jnlabr/ymssp

System identification using associated linear equations J.A. Vazquez Feijoo*, K. Worden, R. Stanway Dynamics Research Group, Department of Mechanical Engineering, University of Sheffield, Mappin Street, Sheffield S1 3 JD, UK Received 24 July 2002; received in revised form 14 May 2003; accepted 19 May 2003

Abstract In this paper, analytical models of non-linear systems are obtained by identifying the frequency response functions (FRFs) of their associated linear equations (ALEs). This allows the use of several methods of identification in the frequency domain usually applicable to linear systems. Among other advantages, the cumbersome multidimensional Fourier Transformation required in higher-order frequency response functions (HFRFs) analysis is eliminated. Two theoretical systems are used here as examples, an electrostrictive actuator and a Duffing oscillator. The concept of the non-linear gain constant arises as a simple means of identification. r 2003 Elsevier Ltd. All rights reserved. Keywords: Non-linear systems; System identification; Volterra series; Higher-order frequency response function; Associated linear equations

1. Introduction The higher-order frequency response functions (HFRFs) describe in the frequency domain, the characteristics of time-invariant continuous non-linear systems. Although a very convenient analytical tool, their practical application is limited to low levels of excitation because higherorder functions require multidimensional Fourier transformation. Obtaining HFRFs of orders higher than the third is time-consuming and the results are often difficult to interpret. In previous work, these systems were modelled by a sequence of continuous linear models named associated linear equations (ALEs) [1]. When the appropriate input signal is used, the output of each ALE is the corresponding order Volterra operator. The relationship between the Volterra operators and the HFRFs are well established. Because all the ALEs are linear models, *Corresponding author. E-mail address: mep97jav@sheffield.ac.uk (J.A. Vazquez Feijoo). 0888-3270/03/$ - see front matter r 2003 Elsevier Ltd. All rights reserved. doi:10.1016/S0888-3270(03)00078-5

ARTICLE IN PRESS J.A. Vazquez Feijoo et al. / Mechanical Systems and Signal Processing 18 (2004) 431–455

432

they each have a standard linear response function that is named here the associated frequency response function (AFRF) to emphasise that they correspond to a specific order ALE. To avoid the technical difficulties in obtaining the HFRFs, it is proposed here to find the AFRF of each ALE. From the AFRF of the nth-order ALE, the corresponding HFRF of the original system is obtained straightforwardly. Also, the associated AFRF is easier to interpret and provides the same information as the HFRF [1]. For any order, by the appropriate spectra from the system, it is possible to identify the AFRF and then the nth-order ALE. Thereafter, when all the Volterra operators are added, an analytical model is obtained. In this paper, two single-input single-output (SISO) theoretical systems are used as examples, a Hammerstein cascade which models an electrostrictive actuator and a Duffing oscillator. Once the spectra are obtained, the system parameters are extracted by a simple curve-fitting method. From the simulation results on both examples, it is emphasised that the difference between the firstorder ALE and all other superior orders is just the amplitude of the output signal. Therefore, the non-linear problem identification is reduced to finding a single scalar named here the non-linear gain constant Kn ; for each order. This parameter can be obtained by only a one-dimensional Fourier transformation independent of the ALE order.

2. Theory For a given non-linear system, using Gaussian white noise as the input and evaluating the appropriate spectra, an estimation of the first two higher-order frequency response functions (HFRFs) can be obtained [2], Syx ðoÞ ; ð1Þ H1 ðoÞDL1 ðoÞ ¼ Sxx ðoÞ H2 ðo1 ; o2 ÞDL2 ðo1 ; o2 Þ ¼

Sy0 xx ðo1 ; o2 Þ ; 2!Sxx ðo1 ÞSxx ðo2 Þ

ð2Þ

where Sxy is the cross power spectrum between the variables yðtÞ and xðtÞ; where xðtÞ represents the system input and yðtÞ represents the system output. The prime on y indicates that the mean of the signal has been removed. In theory, any HFRF can be obtained by the same procedure. In practice, the limitation is on the order of the multidimensional Fourier transformation that can be undertaken. In [1] it was established that if a non-linear system possesses a convergent Volterra representation,1 it can be split into a series of linear subsystems acting in parallel (each one producing a specific order Volterra operator). Those systems are modelled by ALEs. The nonlinear response, which for continuous non-linearities can be represented by polynomial terms, is restricted in the ALE models to act as the excitation in the linear subsystems. Because the subsystems are linear, by the appropriate selection of the input (the appropriate non-linear terms), Eq. (1) can be applied to any higher frequency order ALE and then by only a unidimensional Fourier transformation obtain a linear spectrum named here as the AFRF. For 1

Here for convergence it is understood that a finite number of operators exactly represents the system or at least a truncated series mimics the system up to a predefined accuracy.

ARTICLE IN PRESS J.A. Vazquez Feijoo et al. / Mechanical Systems and Signal Processing 18 (2004) 431–455

433

an SISO linear structural system, the parameters to evaluate can be considered to be reduced to three, the natural frequency o; the damping ratio z or the lost factor (in case of hysteretic damping). The methods presented here can be carried out in two steps, the parameter estimation of the linear part and a second step that is the identification of the non-linear terms that have to be included in the model. The trial input is a product of the non-linear terms to be tested; if a spectrum signal is obtained then it is assumed that this term is part of the system excitation of a certain frequency order. The non-linearity can appear either in the input (i.e. Hammerstein models) or in the output (the Duffing oscillator is a classic example). An example of each case is considered here. In both models viscous damping is assumed, therefore the Kelvin–Voigt [3] model describes the general structure linear system and therefore becomes the skeleton of the ALEs y. 1 þ 2zon y’ 1 þ o2n y1 ¼ AuðtÞ;

ð3Þ

where ui ðtÞ is the trial input in the ALE (for the first-order ALE, u1 ðtÞ is the real system excitation). The equivalent expression for this model in the frequency domain is A : ð4Þ H1 ðoÞ ¼ 2 o þ i2zon o þ o2n With respect to the Hammerstein model, its general polynomial form is n X ai xi : y. þ 2zon y’ þ o2n y ¼

ð5Þ

i¼1

In [4], the ALEs for the structural Hammerstein model have to be y. i þ 2zon y’ i þ o2n yi ¼ ai xi :

ð6Þ

When identifying a system, the trial terms for testing if there is a non-linearity in the input of ith order is ui ðtÞ ¼ xi : In particular, the second-order ALE for the case that there is second-order nonlinearity coming only from the input is2 y. 2 þ 2zon y’ 2 þ o2n y2 ¼ a2 x2 : As the relationship between u2 ðtÞ and y2 ðtÞ is linear, according to Eq. (4), the FRF is a2 H12 ðO2 Þ ¼ ; 2 O2 þ i2zon O2 þ o2n

ð7Þ

ð8Þ

where O2 ¼ o1 þ o2 ; as it is valid for linear systems, Eq. (1) can be used to obtain H12 ðO2 Þ: Sy u ðO2 Þ Sy2 x2 ðO2 Þ H12 ðO2 ÞEL12 ðO2 Þ ¼ 2 2 ¼ ð9Þ Sx2 u2 ðO2 Þ Sx2 x2 ðO2 Þ the modal analysis extract from this signal the parameters value.3 The second-order HFRF for the Hammerstein model can be obtained by Harmonic probing as in [5] when applied to Note that x/yi is not a linear system. However, ui /yi is. It is known [1] that the natural frequency on and the damping ratio z for each ALE are of the same value. However, in the rest of this work it is supposed that each order possesses its own values of these parameters. This is with the intention of obtaining several estimations of these two parameters. 2 3

ARTICLE IN PRESS 434

J.A. Vazquez Feijoo et al. / Mechanical Systems and Signal Processing 18 (2004) 431–455

Eq. (7) yielding H2 ðo1 ; o2 Þ ¼

a2 : ðo1 þ o2 Þ þ i2zon ðo1 þ o2 Þ þ o2n2 2

ð10Þ

Comparison with Eq. (9) gives H2 ðo1 ; o2 Þ ¼ H12 ðO2 Þ and only a first-order Fourier transform is needed to obtain the second-order HFRF or in fact any other order observe that in general for the Hammerstein cascade, Sy xn ðOn Þ Hn ðo1 ; o2 ; y; on Þ ¼ H1n ðOn ÞE n ð11Þ Sxn xn ðOn Þ with On ¼ o1 þ o2 þ ? þ on : When the non-linearity belongs to the output something similar is possible. Consider the Duffing oscillator, k2 k3 1 ð12Þ y. þ 2zon y’ þ o2n y þ y2 þ y3 ¼ x: m m m From [1], the second-order ALE is found to be k1 y. 2 þ 2zon y’ 2 þ o2n y2 ¼  y21 : ð13Þ m If u2 ¼ y21 ; the third-order AFRF can be obtained from Eq. (1) as Sy3 y2 ðO2 Þ 1 : ð14Þ H12 ðO2 ÞE Sy3 y2 ðO2 Þ 1 1

Harmonic probing applied to Eq. (13) gives k2 =m H12 ðO2 Þ ¼ 2 O2 þ i2zon O2 þ o2n

ð15Þ

and to Eq. (12) gives H2 ðo1 ; o2 Þ ¼

ðk2 =mÞH1 ðo1 ÞH1 ðo2 Þ : ðo1 þ o2 Þ2 þ i2zon ðo1 þ o2 Þ þ o2n

From harmonic probing on the first ALE one has 1=m H1 ðoÞ ¼ : 2 o þ i2zon o þ o2n

ð16Þ

ð17Þ

If O2 ¼ o1 þ o2 ; then, H2 ðo1 ; o2 Þ ¼ H12 ðO2 ÞH1 ðo1 ÞH1 ðo2 Þ:

ð18Þ

From Eq. (12), it is clear that the non-linear behaviour comes from the non-linear terms containing the k2 and k3 ; respectively. It is then expected that the excitation of any higher-order ALE have to be any term of the polynomial expansion of yðtÞ when expressed by its Volterra operators, i.e., k2 k2 yðtÞ2 ¼ ð y1 ðtÞ2 þ 2y1 ðtÞy2 ðtÞ þ 2y1 ðtÞy3 ðtÞ þ ? þ y2 ðtÞ2 m m þ 2y2 ðtÞy3 ðtÞ þ ?Þ ð19Þ

ARTICLE IN PRESS J.A. Vazquez Feijoo et al. / Mechanical Systems and Signal Processing 18 (2004) 431–455

435

and k3 k3 yðtÞ3 ¼ ð y1 ðtÞ3 þ 3y1 ðtÞ2 y2 ðtÞ þ 3y1 ðtÞy2 ðtÞ2 þ ? þ y2 ðtÞ3 m m þ 3y2 ðtÞ2 y3 ðtÞ þ ?Þ:

ð20Þ

Many methods have been developed for extracting the system parameters from the spectral analysis of linear structural systems. Any of them, when appropriate for the particular linear spectrum, can be used to identify each AFRF. As the examples provided here are SISO, a simple circle-fitting method is used [6]. In general, for a linear Kelvin–Voigt SISO structural system the mobility is a1 ðoÞ ¼

iAo : o2 þ 2zon oi þ o2n

ð21Þ

The mobility Nyquist plot (Fig. 1) is the familiar circle from which the system parameters can be extracted. The natural frequency on is the frequency diametrically opposite the origin. The damping ratio can be obtained from z¼

ðoa  ob Þ ; 2on

ð22Þ

where oa and ob are any two frequencies that subtend the same angle y from the natural frequency on to the left and to the right (Fig. 1). Finally, from the diameter D of the mobility circle, the modal constant is obtained as ð23Þ

A ¼ 2zon D:

Fig. 1. Mobility circle from the Nyquist plot.

ARTICLE IN PRESS 436

J.A. Vazquez Feijoo et al. / Mechanical Systems and Signal Processing 18 (2004) 431–455

3. Identification of a Duffing oscillator by ALEs 3.1. Description of the Duffing oscillator The Duffing oscillator is modelled as a single-degree-of-freedom (SDOF) soft spring–mass– damper system; the non-linear spring force is obtained as Fk ¼ 2:5  107 d þ 25  1010 d3 :

ð24Þ

The mass is 10 kg and the viscous damping coefficient is 1000 N=ðm=sÞ: The equation of motion is identical to Eq. (12) with k2 ¼ 0: When the data are substituted in this equation, 10yðtÞ . þ 1000yðtÞ ’ þ 2:5  107 yðtÞ  2:5  1010 y3 ¼ xðtÞ:

ð25Þ

From Eq. (25) it is possible to extract directly the natural frequency and the damping ratio: o2n ¼ 2:5  107 N=kg m-on ¼ 1581:1 rad=s;

ð26Þ

2zon ¼ 100 N s=kg m-z ¼ 0:0316:

ð27Þ

As there is no even power of yðtÞ in Eq. (25), all even Volterra operators vanish. And only the terms included in Eq. (20) are going to be the excitation to the higher-order ALEs. From Eq. (12), the lower non-zero order ALEs can be obtained as in [1]: y. 1 ðtÞ þ 100y’ 1 ðtÞ þ 2:5  106 y1 ðtÞ ¼ 0:1xðtÞ;

ð28Þ

y. 3 ðtÞ þ 100y’ 3 ðtÞ þ 2:5  106 y3 ðtÞ ¼ 25  109 y31 ðtÞ;

ð29Þ

y. 5 ðtÞ þ 100y’ 5 ðtÞ þ 2:5  106 y5 ðtÞ ¼ 7:5  109 y21 y3 ðtÞ;

ð30Þ

y. 7 ðtÞ þ 100y’ 7 ðtÞ þ 2:5  106 y7 ðtÞ ¼ 7:5  109 ð y21 y5 ðtÞ þ y23 y1 ðtÞÞ:

ð31Þ

In analogy with the first-order Kelvin–Voigt model (Eq. (3)) in which A is the input amplitude, the input amplitude of nth-order ALE ai is according to the following equations: a1 ¼ 0:1 kg1 ; a3 ¼ 2:5  109 rad=s2 m2 ; a5 ¼ 7:5  109 rad=s2 m2 ; a7 ¼ 7:5  109 rad=s2 m2 : 3.2. Identification process The ALEs are going to be found in ascendant frequency order. The first-order response is obtained as follows; the excitation level is adjusted until a linear response is obtained (see Section 6). For the present example a band-limited white noise is used (though a random signal can provide the appropriate spectrum). Eq. (1) gives the first-order AFRF in receptance form. In order to obtain the mobility the receptance is multiplied by io: Fig. 2 shows the circle obtained. The spectra is obtained from 100 averages of groups of 1024 points; a Hanning window is used.

ARTICLE IN PRESS J.A. Vazquez Feijoo et al. / Mechanical Systems and Signal Processing 18 (2004) 431–455

437

Fig. 2. Circle fitted to the mobility Nyquist plot of the first-order Duffing oscillator response.

The parameters obtained from circle fitting procedure are a1 ¼ 0:0956 kg1 ;

z ¼ 0:0307

and on ¼ 1565 rad=s:

According to Eq. (4), the first-order AFRF in receptance form is4 then H1 ðoÞ ¼

o2

0:0956 : þ 96:14oi þ 2:49331  106

ð32Þ

Applying the inverse Fourier transform gives the first-order ALE: y. 1 ðtÞ þ 96:1914y’ 1 ðtÞ þ 2:4933  105 y1 ðtÞ ¼ 0:0956xðtÞ:

ð33Þ

The next step is to identify the second order. At this point of the identification process, the polynomial non-linear terms that have to be added to the model are unknown. Therefore, it is necessary to test all the possible terms that can produce a second-order response. If the nonlinearity comes from the input, the unique possible non-linear term is u2 ¼ x2

ð34Þ

and if the non-linearity comes from the output, u2 ¼ y21 : 4

ð35Þ

Note that the corresponding H11 ðoÞ is actually the well-known linear FRF; so for this first order there is no need to repeat sub-indexes, and it should be no source of confusion. Denote this function as H1 ðoÞ:

ARTICLE IN PRESS 438

J.A. Vazquez Feijoo et al. / Mechanical Systems and Signal Processing 18 (2004) 431–455

Trying the input from Eq. (35) Eq. (1) is for this case Sy2 ð y2 Þ ðO2 Þ 1 : H12 ðO2 ÞEL1 ðO2 Þ ¼ Sð y2 Þð y2 Þ ðO2 Þ 1

ð36Þ

1

The results obtained for the second-order AFRF are shown in Fig. 3. The real and imaginary parts are observed to have zero mean. This behaviour suggests correctly that there is no secondorder excitation from the output ð y2 ðy21 ðtÞÞ ¼ 0Þ: The mobility Nyquist plot corresponding to Fig. 4 shows that there is no circular structure. A similar result is obtained when the input from Eq. (34) is tested. Therefore there is no second-order response from the system (as this is what is expected when observing Eq. (25)). Considering the third-order response from Eqs. (19) and (20), the input into the third-order ALE is ð37Þ u3 ðtÞ ¼ x3 ðtÞ or u3 ðtÞ ¼ y1 ðtÞ3 þ 2y1 ðtÞy2 ðtÞ:

ð38Þ

As long as there is no second-order component, only the first term contributes to the third-order response of the system; this leads to Syð y3 Þ ðO3 Þ 1 H13 ðO3 ÞEL1 ðO3 Þ ¼ : ð39Þ Sð y3 Þð y3 Þ ðO3 Þ 1

1

Fig. 5 shows a clear linear third-order response and the mobility Nyquist plot (Fig. 6) exhibits a true circular shape. Using the same procedure as in the first-order AFRF, the parameters extracted for the third-order ALE are a3 ¼ 2:659  109 rad=s2 m2 ; z ¼ 0:0ð26Þ8; on ¼ 15068:2 rad=s:

Fig. 3. From Eq. (43), spectra obtained for the second-order response of the Duffing oscillator.

ARTICLE IN PRESS J.A. Vazquez Feijoo et al. / Mechanical Systems and Signal Processing 18 (2004) 431–455

Fig. 4. Second-order mobility Nyquist plot of the Duffing oscillator.

Fig. 5. AFRF of the third-order ALE.

439

ARTICLE IN PRESS 440

J.A. Vazquez Feijoo et al. / Mechanical Systems and Signal Processing 18 (2004) 431–455

Fig. 6. Mobility circle of the AFRF of the third-order ALE of the Duffing oscillator.

The corresponding third-order AFRF is 2:659  109 H13 ðO3 Þ ¼ O23 þ 110iO3 þ 2:5033  106

ð40Þ

and the ALE y. 3 ðtÞ þ 110y’ 3 ðtÞ þ 2:5033  106 y3 ðtÞ ¼ 2:659  109 y1 ðtÞ3 : Harmonic probing on Eq. (25) gives the third-order HFRF as 2:659  109 H1 ðo1 ÞH1 ðo2 ÞH1 ðo3 Þ H3 ðo1 ; o2 ; o3 Þ ¼ : ðo1 þ o2 þ o3 Þ2 þ 110iðo1 þ o2 þ o3 Þ þ 2:5033  106

ð41Þ

ð42Þ

The next ALE for which a true response has been found is of fifth order. The unique response comes from (Eq. (20)) as ð43Þ u5 ðtÞ ¼ y1 ðtÞ2 y3 ðtÞ: The AFRF of the fifth ALE is Sy5 ð y2 y3 Þ ðO5 Þ 1 H15 ðO5 Þ ¼ : Sð y2 y3 Þð y2 y3 Þ ðO5 Þ 1

ð44Þ

1

The corresponding extracted parameters are a5 ¼ 7:4725  109 rad=s2 m2 ; z ¼ 0:0363;

on ¼ 1558:3 rad=s:

So the fifth-order ALE and its AFRF are 7:4865  109 H15 ðO5 Þ ¼ ; 2 O5 þ 114O5 i þ 2:47121  106 y. 5 ðtÞ þ 114y’ 5 ðtÞ þ 2:47121  106 y5 ðtÞ ¼ 7:4865  109 y1 ðtÞ2 y3 ðtÞ:

ð45Þ ð46Þ

ARTICLE IN PRESS J.A. Vazquez Feijoo et al. / Mechanical Systems and Signal Processing 18 (2004) 431–455

441

From Eq. (25), using harmonic probing, the fifth-order HFRF is H5 ðo1 ; o2 ; o3 ; o4 ; o5 Þ ¼ H15 ðo1 þ o2 þ o3 þ o4 þ o5 ÞH1 ðo1 ÞH1 ðo2 ÞH3 ðo3 ; o4 ; o5 Þ:

ð47Þ

Now, in order to measure how well the truncated Volterra series represents the system the mean square error (MSE) is used [8] MSE ¼

N 100 X ðni  n# i Þ2 : Ns2n i¼1

ð48Þ

Here N is the number of data points, s2n the variance of n that is the ‘‘true’’ value and n# is the ‘‘approximate value’’ that for the present case is provided by the truncated Volterra series. The MSE considers the variances and error of the data. The fifth-order Volterra representation gives an error of just 0.2953% for an excitation level up to 0:35 MN:

4. Identification of an electrostrictive actuator by ALEs 4.1. The electrostrictive-system parameters Following Appendix A, the electrostrictive actuator is represented by the equation y. þ 9128:7y’ þ 2:0833  103 y ¼ 2:4306  104 x2  7:027  1022 x6 :

ð49Þ

As in the square law version this system exhibits no first-order component in the output. To obtain a first-order Volterra operator, a DC component is going to be added at the output. This component is obtained by adding a constant signal a in the input [4]; the actual signal uðtÞ into the electrostrictive system is uðtÞ ¼ xðtÞ þ a;

ð50Þ

where a is an appropriate constant that for this case is a ¼ 8 kV: y. þ 9128:7y’ þ 2:0833  103 y ¼ 2:4306  104 ðx þ aÞ2  7:027  1022 ðx þ aÞ6 :

ð51Þ

Developing both parentheses in Eq. (51), one obtains y. þ 9128:7y’ þ 2:0833  103 y ¼ 7:027  1022 x6 þ 6ð2:4306  104 Þax5 þ 15½7:027  1022 a2 x4 þ 20ð7:027  1022 Þx3 a3 þ ð2:4306  104 þ 15½7:027  1022 a4 Þx2 þ ½2ð2:4306  104 Þa þ 6ð7:027  1022 Þa5 x þ ð2:4306  104 a2 þ ð7:027  1022 Þa6 Þ:

ð52Þ

This system possesses a finite Volterra representation of sixth order [4]. As just up to the fourth order is going to be identified here, the corresponding ALEs are y0 ¼ 2:11  102 ;

ð53Þ

y. 1 þ 9128:7y’ 1 þ 2:0833  103 y1 ¼ 3:75x;

ð54Þ

ARTICLE IN PRESS 442

J.A. Vazquez Feijoo et al. / Mechanical Systems and Signal Processing 18 (2004) 431–455

y. 2 þ 9128:7y’ 1 þ 2:0833  103 y1 ¼ 1:9989  104 x2 ;

ð55Þ

y. 3 þ 9128:7y’ 3 þ 2:0833  103 y3 ¼ 1:9989  1043 x3 ;

ð56Þ

y. 4 þ 9128:7y’ 4 þ 2:0833  103 y4 ¼ 6:7459  1013 x4 :

ð57Þ

From Eqs. (53)–(57), the amplitudes to identify are a0 ¼ 2:11  102 m=s2 ;

a1 ¼ 3:75 m=s2 V;

a2 ¼ 1:9989  104 m=s2 V2 ;

a3 ¼ 1:9989  104 m=s2 V3 ;

a4 ¼ 6:7459  1013 m=s2 V4 : 4.2. Identifying the electrostrictive system For this system, once more the input to the system is a band-limited zero-mean white noise. As in the case of the Duffing oscillator, an amplitude at which the system responds linearly is selected (Section 6). From here the next step is to apply Eq. (1) to the data. Observe that although the input is zero-mean, the output possesses a DC component. This mean component is the response to the a which is part of the system. This response is taken into account by a constant y0 that is the zero-order Volterra operator. To evaluate y0 ; the system output is measured when there is no excitation into the system ðxðtÞ ¼ 0Þ obtaining y0 ¼ 7:3632  106 m=s2 :

ð58Þ

The next step is to excite the system by the appropriate amount, to estimate the first-order operator from the output data (Refer to Section 6). Once the linearity is assured, Eq. (1) is applied. The parameters extracted from the Nyquist plot (Fig. 7) are a1 ¼ 3:3944 m=s2 V;

z ¼ 0:1171;

on ¼ 44638 rad=s:

The electrostrictive actuator first-order ALE and its corresponding AFRF when estimated according to Eq. (3) are y. 1 ðtÞ þ 1:0579  104 y’ 1 ðtÞ þ 1:9926  109 y1 ðtÞ ¼ 3:3944xðtÞ; H1 ðoÞ ¼

o2

3:3944 : þ 10579io þ 1:9926  109

ð59Þ ð60Þ

The second-order AFRF is once more tested for u2 ðtÞ ¼ x2 ðtÞ and u2 ðtÞ ¼ y21 ðtÞ: It is found that only the former produces a second-order response as shown in Fig. 8. The parameters extracted are a2 ¼ 1:6394  109 m=s2 V2 ;

z ¼ 0:1278;

on ¼ 43769 rad=s:

The second-order ALE and AFRF according to Eq. (3) are y. 2 þ 11346:4y’ 2 þ 1:928  109 y2 ¼ 1:7794  104 x2 ;

ð61Þ

ARTICLE IN PRESS J.A. Vazquez Feijoo et al. / Mechanical Systems and Signal Processing 18 (2004) 431–455

443

Fig. 7. Circle fitting of the first-order harmonic response of the electrostrictive actuator.

H12 ðO2 Þ ¼

1:7794  104 : O22 þ 11346iO2 þ 1:928  109

ð62Þ

According to Eq. (7), the corresponding HFRF is H2 ðo1 ; o2 Þ ¼ H12 ðO2 Þ ¼

1:7794  109 : ðo1 þ o2 Þ2 þ 11346iðo1 þ o2 Þ þ 1:928  109

ð63Þ

For the third-order ALE, the unique term found that produces a response of this order is u2 ðtÞ ¼ xðtÞ3 ; and the results are a3 ¼ 7:23  109 m=s2 V3 ;

z ¼ 0:11;

on ¼ 43725 rad=s;

so, y. 3 þ 9653:4y’ 3 þ 1:9254  109 y3 ¼ 7:23  109 x3

ð64Þ

and the third-order HFRF is H3 ðo1 ; o2 ; o3 Þ ¼ H13 ðO3 Þ ¼

7:23  109 : ð65Þ ðo1 þ o2 þ o3 Þ2 þ 9653iðo1 þ o2 þ o3 Þ þ 1:9254  109

The fourth-order parameters, ALE and its corresponding HFRF are A4 ¼ 6:6768  1013 m=s2 V3 ;

z ¼ 0:1168;

on ¼ 44533 rad=s;

ARTICLE IN PRESS 444

J.A. Vazquez Feijoo et al. / Mechanical Systems and Signal Processing 18 (2004) 431–455

Fig. 8. Circle fitting of the second-order harmonic response of the electrostrictive actuator.

y. 4 þ 10528y’ 4 þ 1:9832  109 y4 ¼ 6:6908  1013 x4 ;

ð66Þ

H4 ðo1 ; o2 ; o3 ; o4 Þ ¼ H14 ðO4 Þ ¼

6:6908  1013 : ð67Þ ðo1 þ o2 þ o3 þ o4 Þ2 þ 10528iðo1 þ o2 þ o3 þ o4 Þ þ 1:9832  109

The electrostrictive actuator from which the data were modelled suffers saturation at 15 kV: The theoretical model presented here does not consider saturation. This is in order to be able to assess up to which input amplitude a certain order Volterra system provides an accurate representation. Here, the fourth-order Volterra representation proved to be accurate up to a 1% mean square error for an input amplitude of 20 kV: This voltage is higher than the saturation, therefore no higher-order Volterra representation may be needed.

5. The non-linear gain constant From the results of the two cases presented here, and considering that variations in the parameter estimates appear as a consequence of the curve-fitting process, all the ALEs show the same natural frequency on and the same damping ratio z as the theory predicts [1]. Table 1

ARTICLE IN PRESS J.A. Vazquez Feijoo et al. / Mechanical Systems and Signal Processing 18 (2004) 431–455

445

Table 1 Average of the system parameters obtained from the ALEs Duffing oscillator

on z

Electrostrictive actuator

Average

e%

True value

Average

e%

True value

1577.8 0.03293

0.2 4.2

1581.1 0.0316

44239.75 0.119

3.1 19

45643.54 0.1

contains the average value of the ALE parameters for both cases and the mean square error compared to the actual parameters. As on and z are the same, independent of the ALE order (e.a. observe Eqs. (16) and (17)), H1n ðoÞ an ¼ ¼ Kn ; ð68Þ H1 ðoÞ a1 where Kn defined here, is the non-linear gain constant, which can be used for identification as follows. For a given system, two different inputs xðtÞ and x0 ðtÞ are tried. From the first excitation xðtÞ; the nth-order AFRF is according to Eq. (1): Sy u ðOn Þ ; ð69Þ H1n ðOn Þ ¼ n n Sun un ðOn Þ where un is the excitation to the nth-order ALE produced when the system input is xðtÞ: At this stage of the identification process, the first-order ALE is already known, then this first-order response can be simulated as an independent linear system. The second excitation x0 ðtÞ goes into this system and is made numerically equal to x0 ðtÞ ¼ un : The corresponding first-order AFRF is obtained as Sy0 un ðOn Þ : ð70Þ H10 ðOn Þ ¼ 1 Sun un ðOn Þ As the system parameters do not depend on the input and just on the system characteristics, when Eqs. (69) and (70) are substituted into Eq. (68), the non-linear gain constant Kn can be obtained as Kn ¼

Syn un ðOn Þ=Sun un ðOn Þ ; Sy01 un ðOn Þ=Sun un ðOn Þ

ð71Þ

Kn ¼

Syn un ðOn Þ : Sy01 un ðOn Þ

ð72Þ

so,

Once more all the possible terms of second order, third, fourth etc., have to be tested. From here on only the terms that produce an nth-order output component are going to be mentioned. For the second order K2 ¼

Sy2 x2 ðOn Þ : Sy01 x2 ðOn Þ

ð73Þ

ARTICLE IN PRESS 446

J.A. Vazquez Feijoo et al. / Mechanical Systems and Signal Processing 18 (2004) 431–455

The first step is to estimate y1 ðtÞ: The second step is to make x0 ðtÞ ¼ u2 ¼ x2 ðtÞ and analytically compute y01 ½u2 ðtÞ : Then apply Eq. (72). The inputs for the third and fourth order are, respectively, x0 ðtÞ ¼ xðtÞ3

ð74Þ

x0 ðtÞ ¼ xðtÞ4 :

ð75Þ

and

The results from the simulation are plotted in Figs. 9 and 10. The non-linear gain constant Kn is obtained from a mean value calculated around the resonance zone. The same figures include two extra lines that show the 95% confidence that the true value falls into the interval enclosed by these lines. This interval is calculated as if the spectral distribution were Gaussian. If kn is the spectra mean, then, 1:96spKn  kn p1:96s; where s is the standard deviation of the spectrum. According to Eq. (68), the amplitude ratio an can then be obtained as a n ¼ Kn a 1 :

ð76Þ

Table 2 shows the Kn obtained and its corresponding an : Observe that the true value resides within the 95% confidence. The non-linear constants of the Duffing oscillator are obtained by the same procedure.

Fig. 9. Results obtained from Eq. (76) for the second-order non-linear gain constant K2 of the electrostrictive actuator.

ARTICLE IN PRESS J.A. Vazquez Feijoo et al. / Mechanical Systems and Signal Processing 18 (2004) 431–455

447

Fig. 10. Results obtained from Eq. (76) for the third-order non-linear gain constant K3 of the electrostrictive actuator.

Table 2 Non-linear gain constants of the electrostrictive actuator and the amplitude ratio obtained from them

Kn calculated Kn real an calculated

Second-order

Third-order

Fourth-order

5:2075  105 5:3304  105 1:7915  104

2:1779  109 1:91  109 8:6  109

2:0708  1013 1:799  1013 7:7655  1013

Figs. 11 and 12 show a sample of the results for the Duffing oscillator and Table 3 shows the amplitude ratio an obtained from Kn :

6. Obtaining the appropriate amplitude for the Nth-order spectral response This section deals with the criteria to select the proper amplitude for isolating the higher-order responses. In both methods used here, the process of identification goes from the lower to the higher-order response. In order to obtain the first Volterra operator, an amplitude should be found that guarantees that the system response is dominantly linear. There are several methodologies to determine if the system is linear or not, i.e. the correlation function [7,9] or the Hilbert transformation [6]. The former has some restrictions that prevent it from being used

ARTICLE IN PRESS 448

J.A. Vazquez Feijoo et al. / Mechanical Systems and Signal Processing 18 (2004) 431–455

Fig. 11. Results obtained from Eq. (76) for the third-order non-linear gain constant K3 of the Duffing oscillator.

Fig. 12. Results obtained from Eq. (76) for the fifth-order non-linear gain constant K5 of the Duffing oscillator.

ARTICLE IN PRESS J.A. Vazquez Feijoo et al. / Mechanical Systems and Signal Processing 18 (2004) 431–455

449

for all cases i.e. the signal should have even order components. Hilbert transformation requires the previous estimation of a Fourier transformation. An objective way of assessing linearity is to verify the proportional property of the linear systems. The steps taken here are to excite the system at two levels, X0 ; and 2X0 ; then obtain the MSE between the two responses. The input is a band-limited white noise. The first-order amplitude is straightforward to obtain in both cases. The AFRF of the ALEs is more difficult to isolate as the other order increases. On the one hand, if the input amplitude into the system is too low, the higher-order response can be indistinguishable from the noise and on the other, the magnitude of the harmonics that the non-linear systems produces scales with a power of the corresponding order of the input amplitude [1]. The relative magnitude between consecutive non-linear operators is more significant as the input amplitude increases. If the system is free of noise, small amplitudes should produce well-separated non-linear-order components, but with noise, the minimum amplitude is limited. The best results are found when

Table 3 Non-linear gain constants of the Duffing oscillator and the amplitude ratio obtained from them

Kn calculated Kn true an calculated

Third-order

Fifth-order

Seventh-order

2:51  1010 25  1010 2:39  109

7:6843  1010 7:5  1010 7:3446  109

7:8131  1010 7:5  1010 7:4693  109

Fig. 13. AFRF models obtained for H13 ðO2 Þ from the frequency spectra analysis obtained from circle fitting, the harmonic constant K3 of the Duffing oscillator.

ARTICLE IN PRESS 450

J.A. Vazquez Feijoo et al. / Mechanical Systems and Signal Processing 18 (2004) 431–455

Fig. 14. AFRF models obtained for H15 ðO3 Þ from the frequency spectra analysis obtained from circle fitting, the harmonic constant K5 of the Duffing oscillator.

Fig. 15. AFRF models obtained for H12 ðO2 Þ from the frequency spectra analysis obtained from circle fitting, the harmonic constant K2 of the electrostrictive actuator.

ARTICLE IN PRESS J.A. Vazquez Feijoo et al. / Mechanical Systems and Signal Processing 18 (2004) 431–455

451

Fig. 16. AFRF models obtained for H13 ðO3 Þ from the frequency spectra analysis obtained from circle fitting, the harmonic constant K2 of the electrostrictive actuator.

the MSE between 2yðxðtÞÞ and yð2xðtÞÞ is of around 5%; this corresponds to a root-mean-square (RMS) of 1070 kN for the Duffing oscillator and of 4 kV for the electrostrictive actuator. The data obtained are the sum of all the Volterra operators (Eq. (39)). So from these data it is necessary to estimate the individual values. The estimation is carried out in ascending order according to the following formulas. First order: y1 ðtÞDyðtÞ  y0 ;

ð77Þ

Second-order ALE: y2 ðtÞDyðtÞ  y0  y1 ðtÞ;

ð78Þ

Third order: y3 ðtÞ ¼ yðtÞ  y0  y1 ðtÞ  y2 ðtÞ;

ð79Þ

etc. For the Duffing oscillator neglect the term containing y0 as this system does not possess any DC term. The evaluation of both methods can be assessed observing Figs. 13–16 where the estimated AFRF with model fits are shown. Also, the mean square error of each model with respect to the spectra is shown in Table 4.

ARTICLE IN PRESS 452

J.A. Vazquez Feijoo et al. / Mechanical Systems and Signal Processing 18 (2004) 431–455

Table 4 Mean square error obtained from the models of the Duffing oscillator and the electrostrictive actuator Duffing oscillator Model Circle fitting Non-linear gain constant Kn

H13 0.0315 0.0351

Electrostrictive actuator H15 0.3897 0.3771

H12 0.3897 0.3771

H13 0.9918 0.6256

H14 1.5678 1.5709

7. Discussion and conclusions If a non-linear system possesses convergent Volterra series representation, it can be identified using the decomposition in ALEs. For this objective, two different approaches have been used here; to find each AFRF and the non-linear gain constant Kn directly (the ratio between the nth-order AFRF and the first-order one). Each method possesses its own strengths and weaknesses. Because all the higher-order ALEs have the same damping ratio and natural frequency, the direct AFRF identification procedure offers a means to evaluate both parameters several times. In the examples presented here, four higher-order ALEs were identified, so that for each system there were four different evaluations of the damping ratio z and the natural frequency on from different data groups (each group is produced for a specific order signal). As the unique different parameter between the ALEs is the modal constant, the non-linear gain constant offers a second evaluation of the amplitude ratio an : For the Duffing oscillator, both methods seem to produce accurate results, very close to each other. The electrostrictive actuator results seem to be good for the first two orders, the third- and fourth-order ALEs are a little harder to obtain. The accuracy in the case of the Duffing oscillator is observed to be higher, this could be because of the nature of the systems. The Duffing oscillator possesses a ‘‘natural’’ first-order response and the same can be said for higher-order responses. The electrostrictive actuator is forced to produce a first-order output by adding the DC. This characteristic makes the response of the high-order ALEs a little closer in magnitude to each other, leaving a smaller amplitude interval to isolate a particular-order response. To be able to apply these methods, it is necessary to have the possibility to test the system in several levels of excitation until the proper level for testing is obtained (Section 6). The method’s capability of obtaining higher frequency order responses lay strongly in an accurate first-order modelling. This is because the trial inputs for the non-linear part are conformed by products of lower order signal produced by simulation. As the order becomes higher, the error tends to find a clean spectrum. Appendix A. The Hammerstein model of an electrostrictive actuator The theoretical electrostrictive actuator is based on the characteristics of the Lead-LanthanumZirconium-Titanate (PLZT) electostrictive material. When the grain size is of the order of 4:5 mm;

ARTICLE IN PRESS J.A. Vazquez Feijoo et al. / Mechanical Systems and Signal Processing 18 (2004) 431–455

453

its strain-electric field e–E curve in the d33 direction is shown in Fig. 17 that also contains the average strain for a given electric field E: Its saturation point is at a field of 22 kV=cm: It is usual to model the electrostrictive effect using the quadratic law as eðtÞ ¼ d33 Ef2 ðtÞ;

ðA:1Þ

where d33 is the electrostrictive constant when the electric field is applied in the same direction as the controlled deformation. Considering the mid-points curve in Fig. 17, it is found that d33 ¼ 0:13 cm2 =kV2 : However, no second-order law is able to model the electrostrictive response up to saturation. It is found that a sixth order is the lower-order approximation that produces a good approximation: 0 Ef2 þ fEf6 e ¼ d33

ðA:2Þ

by interpolation of the mid-points e ¼ 0:92  106 Ef6 þ 0:14Ef2 :

ðA:3Þ

The voltage x needed to obtain a specific electric field E is x ¼ EL; where L is the length of the electroceramic in the 3-direction. The electrostrictive actuator system model is shown in Fig. 18 under an external excitation force F ðtÞ: The structural stiffness is represented by k and the electrostrictive stiffness is k0 : The viscous dashpot and the stiffness of other structural elements are represented by well-known mathematical models. The electrostrictive effect does not produce stress in the ceramic, but simply changes the length of the ceramic element. The stress is produced because of the stiffnessstrain that is obtained from the difference of the actual length of the ceramic at the time t and the

Fig. 17. Strain–electric field response of the PLZT: the quadratic law partially models the electrostrictive effect.

ARTICLE IN PRESS 454

J.A. Vazquez Feijoo et al. / Mechanical Systems and Signal Processing 18 (2004) 431–455

Fig. 18. 2A. The Wiener model.

position yðtÞ of the mass m at the same time. Then, F 0 ¼ k0 ðy  DðtÞÞ;

ðA:4Þ

where DðtÞ is the increment in length of the ceramic produced by the electrostrictive effect: DðtÞ ¼ eðtÞL:

ðA:5Þ

Here eðtÞ is the strain obtained from Eq. (A.1), and L is the ceramic length in the deformation direction. Adding all the forces in the direction of y; one obtains Fk þ Fct þ F 0 ¼ Fi ;

ðA:6Þ

where Fi is the inertia force, Fk is the equivalent structural stiffness of all other elements apart from the electroceramic, Fc is the viscous damping force and F 0 is the force exerted by the ceramic stress–strain. When expressions for all the forces are substituted in Eq. (A.6) one obtains my. þ cy’ þ ðkt þ k0 Þy ¼ k0 Di :

ðA:7Þ

The equivalent stiffness of an axially loaded member is EA k0 ¼ ; ðA:8Þ L where E is the Young’s modulus. Substituting Eqs. (A.1), (A.5) and (A.8) into (A.7) we obtain EA my. þ cy’ þ ðkt þ k0 Þy ¼ ðA:9Þ d33 Ef2 Li : L Finally, using (A.4) and the total system stiffness k; EA my. þ cy’ þ ky ¼ 2 d33 xðtÞ2 : ðA:10Þ L Dividing by the mass, an equation that agrees with Eq. (3) is obtained: EA ðA:11Þ y. þ 2zon y’ þ o2n y ¼ 2 d33 xðtÞ2 : Lm For a more accurate model for the electrostrictive effect, the sixth-degree polynomial in Eq. (A.2) is used. In the same way as for the square law, the corresponding Hammerstein model is

ARTICLE IN PRESS J.A. Vazquez Feijoo et al. / Mechanical Systems and Signal Processing 18 (2004) 431–455

455

obtained as 0 k0 d33 k0 f ðA:12Þ x2 þ 5 x6 : Lm Lm 0 The data used in this equation are, L ¼ 12 mm; A ¼ 0:3  4 mm; m ¼ 4:8 g; E ¼ 100 GPa; d33 ¼ 2 6 2 4 2 2 10 6 41 6 cm =kV ¼ 9:2  10 m =V6 : 0:14  10 cm =kV ¼ 1:4e  15 m =V and f ¼ 0:78  10 In the particular case that the mass m is suspended by the electrostrictive plate, Fkt is zero and then, EA k ¼ k0 ¼ ¼ 107 N=m: L Then, rffiffiffiffi k ¼ 45643:54 rad=s: on ¼ m The final version of the electrostrictive model is therefore,

y. þ 2zon y’ þ o2n y ¼

y. þ 9128:7y’ þ 2:0833  109 y ¼ 2:4306  104 x2  7:027  1022 x6 :

ðA:13Þ

References [1] J.A. Vazquez-Feijoo, K. Worden, R. Stanway, On perturbation analysis and the Volterra series, Proceedings of the Seventh International Conference of Structural Dynamics; Recent Advances, University of Southampton, 2000, pp. 415–425. [2] S.J. Gifford, Volterra series analysis of non-linear structures, Ph.D. Thesis, University of Edingorgh, 1975. [3] C.W. Bert, Material damping: an introductory review of mathematical models, measures and experimental techniques, Journal of Sound and Vibration 29 (2) (1959) 129–153. [4] J.A. Vazquez-Feijoo, K. Worden, R. Stanway, Linearization of electrostrictive-type systems, Proceedings from the Conference on Healt Monitoring and System Identification, Madrid, Spain, 2000, pp. 115–126. [5] I.J. Leontaritis, S.A. Billings, Input–output parametric models for non-linear systems; Part I: deterministic nonlinear systems, International Journal of Control 41 (1971) 303–328. [6] E.J. Ewins, Modal Testing: Theory and Practice, Research Studies Press Ltd., John Willey, NY, 1970. [7] K. Worden, G.R. Tomlinson, Nonlinearity in Structural Dynamics, Detection, Identification and Modelling, Institute of Physics Publishing, Bristol, Philadelphia, 2001. [8] J.S. Bendat, Nonlinear Systems Techniques and Applications, Wiley, Interscience, New York, 1998. [9] S.A. Billings, K.M. Tsang, Spectral analysis of block-structured non-linear system, Mechanics and Systems Signal Processing 4 (1978) 117–130.