- Email: [email protected]
Identification of a Hammerstein model for an aerospace electrohydraulic servovalve
19th IFAC Symposium on Automatic Control in Aerospace September 2-6, 2013. Würzburg, Germany
Identification of a Hammerstein model for an aerospace electrohydraulic servovalve Alexandro G. Brito ∗ Waldemar C. Leite Filho ∗ Elder M. Hemerly ∗∗ ∗
Institute of Aeronautics and Space, S˜ ao Jos´e dos Campos, SP, CEP:12220-904, Brazil (e-mails: [email protected] / [email protected]). ∗∗ Technological Institute of Aeronautics, S˜ ao Jos´e dos Campos, SP, CEP:12220-900, Brazil (e-mail: [email protected]) Abstract: This paper discusses the black-box identification of a real aerospace electrohydraulic actuator. A Hammerstein model structure is adopted due to its simplicity and capacity to represent a wide set of nonlinear phenomena. The identification strategy begins with a frequencydomain analysis so that a polynomial NARX model structure can be defined. Then, an extended least squares procedure is applied to obtain the coefficient values. Keywords: NARX/NARMAX models; FRF’s; Multisine excitation; Aerospace actuators. 1. INTRODUCTION
2. NARMAX MODELS
The electrohydraulic servovalves are widely applied to a variety of applications and the knowledge of its dynamics is essential to the control system design. Several works cover the aspects about these devices, including their internal moving parts, the motion equations and the most important nonlinearities [Merritt (1967); Jelali and Kroll (2003)]. Servovalves are especially useful as the main actuator element of aerospace systems; particularly in launchers, they operate in the thrust vector control of solid propelled vehicle or in the gimbaled chamber for liquid propelled rockets. Since their influence over the aerospace vehicle control loops, a proper modeling of the actuator linear and nonlinear behaviors is of primary importance to the launcher performance.
The NARMAX (Nonlinear AutoRegressive with Moving Average and eXogeneous inputs) models proposed by Leontaritis & Billings (1985) is a very general nonlinear representation, which includes the Volterra and Wiener series, the Bilinear, output-affine, rational and other model classes [Chen and Billings (1989)]. To define such NARMAX models, consider u(k) and y(k) the input and output signals of a general nonlinear system, respectively, each belonging to the real Banach spaces U and Y. Assume Uk = {(u(1), . . . , u(k)) : u(i) ∈ U, i = 1, . . . , k}, and the function response that describes the input-output dynamical behavior fx0 : Uk → Y, beginning from an initial state x0 . If the following conditions:
In this paper, a full discussion about the electrohydraulic actuator black-box modeling is performed. As an example, the servovalve used in the attitude control of the VLS (the Brazilian launch vehicle) is considered. The methodology is based on an initial frequency-domain analysis so that the model structure can be inferred. Then, a discrete Hammerstein formulation based on the NARMAX models is fitted. Finally, a model validation step is performed. The paper is organized as follows: the section 2 discusses the definitions associated to the NARMAX model structure selection and identification. In section 3, a whole discussion about the mathematical modeling of the electrohydraulic actuator is presented. Section 4 elaborates all the concepts used for the model identification, including the structure selection, the numerical procedure and the validation. Finally, some conclusions are presented in Section 5.
i) fx0 is a finitely realizable function response, and; ii) a linearized model exists when the system is operated in the neighborhood of an equilibrium point, hold, then there exists a nonlinear function F (·) so that y(k) = F (y(k − 1), . . . , y(k − ny ), u(k − 1), . . . , u(k − nu )). (1) Despite of the wide range of structures that the nonlinear functional F (·) can assume, only the input-output polynomial NARMAX models are discussed in this paper, since that the aim is to obtain a linear-in-the-parameter representation. A polynomial NARMAX model can be constructed by adopting the following format for the function F (·) y(k) =
y ,nu d X m nX X
m=0 p=0 n1 ,nm p Y
y(k − ni )
×
⋆ We would like to thank SIA/FINEP and FUNDEP by the financial support.
978-3-902823-46-5/2013 © IFAC
459
cp,m−p (n1 , · · · , nm ) ×
i=1
m Y
u(k − ni )
(2)
i=p+1
where d is the polynomial degree and 10.3182/20130902-5-DE-2040.00119
2013 IFAC ACA September 2-6, 2013. Würzburg, Germany
ny ,nu
X
≡
n1 ,nm
ny ny X X
n1 =1 n2 =1
···
nu X
.
(3)
nm =1
It is easy to see that the model (2) is linear in the parameters cp,m−p , thus a Least Squares identification procedure can be used. Initially, the polynomial is written as y(k) = Ψ(k − 1, :)θ + ξ(k) nθ X θi ψi (k − 1) + ξ(k) =
(4)
are functions of valve displacement xv which is related to the input voltage. The mechanical actuator mounting can introduce nonlinearities in this displacement, mainly deadzones. The load flow as a function of valve position and load pressure is given by the nonlinear relation s 1 xv QL = Cd wxv (PS − PL ) (6) ρ |xv | where PS is the supply pressure and w ≡ ∂A/∂xv ([w] = m2 /m) is the area gradient of the valve.
i=1
where the regressor matrix Ψ is constructed with all the monomial of (2). By using N measurement samples, (4) becomes a matrix equation. Due to the measurement noise, the parameters θi cannot be calculated through the ordinary least squares procedure, otherwise these parameters can be biased. Instead, other alternatives such as the Instrumental Variables or the Extended Least Squares should be considered [Ljung (1987)]. All the following identification results are obtained by adopting this last strategy. 3. ELECTROHYDRAULIC ACTUATORS The EHA consists of a powerful component in control systems that blends the versatility of electrical components with hydraulic actuators performance at high power levels. In this way, several dynamic effects – including the nonlinear one – are present and one should pay attention to some details in the modeling procedure. Fig. 1 presents a simplified block diagram for this system. The servovalve stage converts the input voltage (Vi ), the supply pressure (PS ) and the load pressure (PL ) in a load flow (QL ) to be applied to the hydraulic cylinder. This cylinder is in charge of producing the load force (fL ) necessary to provide the load stroke displacement (xp ). A thorough modeling discussion for the electrohydraulic actuator is present in Merritt (1967). In this paper, only the most important nonlinear phenomena are presented.
Fig. 1.
The fluid flowing through uncompensated valve orifices causes forces with a direction such that it tends to close the valve port. The magnitude of this force is given by F1 = 2Cd Cv A0 (P1 − P2 )cos(θ) (7) ◦ where, typically θ ≈ 69 is the jet angle of vena contracta, Cd ≈ 0.61 is the discharge coefficient, Cv ≈ 0.98 is the empirical factor called velocity coefficient and ∆P = P1 − P2 the pressure drop. By using Eq. 7 and the numerical values above, one can yield the usual form of the steadystate flow force equation F1 = 0.43w∆P xv = Kf xv . (8) Merritt (1967) comments that on larger single stage EHA, this force can exceed 20lbs. (9kgf ). The reduction of this steady-force flow force can be obtained by using twostage configuration or geometric compensating techniques. These compensations can lead to a nonlinear function between the flow force and the stroke displacement. 3.2 The cylinder dynamics The continuity equation combined with the equation of state (ρ = ρ(P, T )) leads to the expression for the load flow in the cylinder chamber dxp Vp dPL QL = Ap + (9) dt 4β dt where β is the Bulk modulus (for mineral oils and for common values for pressure and temperature, β is typically 1400 to 1600 MPa) and Ap and Vp are, respectively, the section area and the chamber volume of the cylinder. The term Vp /(4β) is equivalent to a linear hydraulic capacitance, known as CH . Since the servovalve is generally attached to the cylinder, the high pressure lines are sufficiently short and additional hydraulic capacitances due to these lines can be neglected. Finally, the load force (fL ) is given by fL = Ap PL .
Electrohydraulic actuator block diagram
3.1 The servovalve dynamics
(10)
4. EXPERIMENTAL RESULTS
The turbulent flow through an orifice occurs at high Reynolds number and is modeled by applying the Bernoulli’s equation. This analysis yields the well-known volumetric flow rate r 2 Q = Cd A0 (P1 − P2 ) (5) ρ where Cd is the discharge coefficient, A0 the orifice area, ρ the mass density of fluid and (P1 − P2 ) the pressure drop. Considering a typical four-way spool valve where the orifice areas depend on the valve geometry, their four areas 460
This paper focuses on the identification of an electrohydraulic actuator used in the attitude control of VLS, the Brazilian Launcher. The objective is to explain some nonlinear behaviors that occurred in previous flights, such as unexpected high amplitude limit-cycles. Initially, a previous analysis in the frequency domain is performed before the identification step. This procedure gives useful information about the system nonlinear behavior to the user, so that he can choose a proper Hammerstein
2013 IFAC ACA September 2-6, 2013. Würzburg, Germany
structure. In this way, the further identification process has a higher probability of success. The nonlinearity detection approach adopt in this paper is based on the frequency domain. Other technique to analyze the presence of nonlinear behavior in the frequency domain is to apply the coherence concept. However, this tool only indicates the presence of nonlinearities, but it is not able to inform which is the involved type of nonlinearities. 4.1 Nonlinearity analysis in the frequency domain For better results for the nonlinearity analysis, the user should choose a proper excitation signal so that the main nonlinear behavior can be detected. A class of excitation commonly used for this task is the odd-odd multisine signal, which has been used in several practical problems involving the identification of nonlinear systems [Pintelon and Schoukens (2001)]. Let u(k) be a discrete periodic input with fundamental frequency ω0 and a sampling frequency ωs = Nu ω0 , where Nu is the data length. Denote I = {l | l = 1, 3, 5, · · · , N < Nu /2} the set of all the odd numbers up to N and F ⊂ I a subset with some of these odd numbers previously chosen by the user. Then, the oddodd multisine signal is defined as X 2πk u(k) = Al sin lω0 + φl , φl ∼ U (−π, π). (11) ωs
lines ωe and at some even non-excited frequency lines ωp ; iii) system containing only odd nonlinearities: then, there exists useful output power at the excited frequency lines ωe and at some odd non-excited frequency lines ωi ; iv) system contains both odd and even nonlinearities: then, there exists useful output power at the excited frequency lines ωe , at some odd non-excited frequency lines ωi and at some even non-excited frequency lines ωp . Then, the nonlinearity detection is performed as follows. Firstly, an odd-odd multisinusoidal input is applied and both the input and the output signals are acquired. The next step is to calculate and analyze the Fourier transform of these signals. If the system contains useful power spectrum only at the excited lines from the input, it is linear. If the system also exhibits useful power spectra at the odd non-excited frequency lines, then it contains some odd nonlinearity such as cubic terms, saturation, etc. Finally, if it has useful power spectrum at the even non-excited frequency lines, the system contains some even nonlinearity such as quadratic terms, absolute-value function, etc. For the electrohydraulic actuator discussed in this paper, an odd-odd multisine with the characteristics presented in Table. 1 was applied. Table 1.
l∈F
where U (−π, π) denotes a uniformly distributed random number between −π and π. A sketch of this signal is presented in Fig. 2. It is important to note that, in spite of
⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆
Odd-odd multisine signal to be used in the actuator nonlinear analysis
sampling frequency fs = 40Hz fundamental frequency f0 = 0.005Hz start excitation frequency f1 = 0.005Hz stop excitation frequency fN = 10Hz 1000 odd frequency lines 800 effectively excited freq. lines (randomly chosen) input RMS value uRM S = 2
By applying this signal to the VLS actuator, one can obtain the output power spectra presented in Fig. 3. Some conclusions can be drawn from it. Firstly, the linear Output Power Spectrum
Fig. 2.
−20
Representation of the Odd-Odd Multisine Signal.
excit. even odd
the phase is a random number for each harmonic, y(k) is a deterministic signal since this random phase is maintained constant along the signal realization.
Anmplitude (dB)
−40 −60 −80 −100
i) pure linear system: then, there exists useful output power only at the excited frequency lines ωe ; ii) system containing only even nonlinearities: then, there exists useful output power at excited frequency 461
−120 −3 10
−2
10
−1
0
10
10
1
10
2
10
Coherence 1 0.8 0.6 g(f)
The nonlinearity detection is based on the behavior of the nonlinear systems when a multisinusoidal input is applied. This topic is richly discussed by Vanhoenacker et al. (2001) and it is only summarized herein. Consider an input-output system and odd-odd multisine signal input as described by (2). Define the set M = {m | m = N, N + 2, · · · , 2N − 1} and denote the excited frequency lines as ωe = {lω0 | l ∈ F}, the odd non-excited frequency lines as ωi = {lω0 | l ∈ (F ∪ M), F ∪ F = I} and the even nonexcited frequency lines as ωp = {lω0 | l = 0, 2, · · · , 2N }. Then, the system output satisfies one of the conditions:
0.4 0.2 0 −3 10
−2
10
−1
0
10
10
1
10
2
10
Freq (Hz)
Fig. 3.
Nonlinearity detection. Top: the output power spectrum. Bottom: The coherence for the output spectrum
contributions, which are represented as the blue crosses,
2013 IFAC ACA September 2-6, 2013. Würzburg, Germany
are in accordance with the expected one for this type of servovalve. However, it is easy to see that several nonexcited frequency lines present useful power amplitude, what indicates the presence of nonlinearities. This is especially true for the odd non-excited frequency lines – marked as red circles. Then, one can infer that the actuator has an odd nonlinear nature, according to the discussed above. In the bottom plot of Fig. 3, the coherence function for the output Fourier transform is presented. A well-known result is that the coherence function is unitary for the entire spectrum if the system is linear. Notice that this is false for the present case, mainly for the odd-nonexcited frequency lines where this value is considerably dropped. This confirms the nonlinear behavior, however is too difficult to infer the type of nonlinearities by a simple inspection of the coherence plot. 4.2 Time-domain model identification Since the conclusions drew from the last section, it is expected that the odd behavior dominates the nonlinear contributions. This information is very useful to delimitate the class of models to be considered before the application of any identification method to the data.
y = Ψθ + e
y = [ y(1) y(2) · · · y(L) ]T , e = [ e(1) e(2) · · · e(L) ]T , θ = [ a 1 a 2 b 1 b 2 b 3 ν 1 ν 2 ν 3 ν 4 ]T , and each column (regressor ) of Ψ is formed by the respective monomial of (14) for k = 1, · · · , L. The model estimation through the traditional least squares, although possible, will provide biased coefficients due to the output measurement noise. The alternative is to use the extended least squares where some regressors related to the noise e(k) are also included in Ψ so that the model coefficient identification is not affected by the correlation with the noise. In this way, the biasness is avoided. By applying this approach to the actuator data, one can obtain the coefficient values shown in Table 2. The extended least squares procedure considered 10 additional noise regressors in Ψ in an iterative procedure repeated until the coefficient convergence. Table 2.
⋆ Least Squares identification
Model Simulation
Amplitude
adopted for the both parts. Due to the odd nature of the electrohydraulic actuator (Fig. 3), an odd fifth degree polynomial is adopted: f (u(k)) = u(k) + γ3 u3 (k) + γ5 u5 (k). (12) The linear part was chosen as a third order discrete transfer function of format z(a1 z + a2 ) . (13) H(z) = 3 z − b1 z 2 − b2 z − b3 The difference equation related to the Hammerstein model in Fig. 4 is given by y(k) =a1 u(k − 1) + ν1 u3 (k − 1) + ν2 u5 (k − 1) a2 u(k − 2) + ν3 u3 (k − 2) + ν4 u5 (k − 2) b1 y(k − 1) + b2 y(k − 2) + b3 y(k − 3), (14) where the coefficients ν are formed by the product between the terms a and γ of the respective monomial.
462
est
0
95
96
97
98 Time (s)
99
100
101
0.3 Cξ,y
0.2
Cξ2,y2
0.1 0 −0.1
0
2
4
6
8
10
12
14
1.5 Cξ,ξ
1
Cξ2,ξ2
0.5 0 −0.5
Eq. (14) is linear in the parameters, hence a least squares identification procedure can be applied. By considering L data points and an output measurement noise e(k), (14) can be written as
y y
0.5
−0.5
Correlation
A hammerstein model.
−1.8096 × 10−3 1.8155 × 10−5 3.3126 × 10−4 −6.0943 × 10−6
ν1 ν2 ν3 ν4
The identified model was initially simulated with the same data set to study the estimation error behavior. The results are summarized in Fig. 5. In the top plot, the model output response is compared to the real data. Notice that both data present similar behavior. The other two plots represent, respectively, the correlation between the output and the identification residues (middle), and the residue cross-correlation (bottom). One can note a low residue cross-correlation, but still higher than the confidence limits. The same can be seen for the correlation between the model output and the identification residues.
Correlation
Fig. 4.
The model coefficients for the VLS electrohydraulic actuator.
9.3872 × 10−2 −7.2915 × 10−3 6.2472 × 10−1 1.3682 × 10−1 −1.0478 × 10−1
a1 a2 b1 b2 b3
⋆ Model structure selection A Hammerstein model is characterized by a static nonlinearity followed by a discrete linear part, as showed in Fig. 4. For this paper, a polynomial representation is
(15)
where
0
2
4
6
8
10
12
14
Fig. 5. Model response and residue analysis (95 % confidence limits). Due to the results above, one could conclude that the model is not adequate because the correlation levels are
2013 IFAC ACA September 2-6, 2013. Würzburg, Germany
Model Validation
1.5
Amplitude
est
0 −0.5 97
97.5
98
98.5
99 Time (s)
99.5
100
100.5
101
Correlation
0.3 Cξ,y
0.2
Cξ2,y2
0.1 0 −0.1
0
2
4
6
8
10
12
14
1.5 Cξ,ξ
1
Cξ2,ξ2
0.5 0 −0.5
Signal analysis
y y
0.5
Correlation
higher than the confidence bounds. A well-known identification result is that if the residue correlation function is significantly out of confidence bounds, this means that unmodeled dynamics still persist on the data. However, these conclusions can be only drawn if the input signal has a significant random aspect. To gain more information about the behavior in Fig. 5, an input correlation analysis was performed. The results are presented in Fig. 6. In the top plot, the input cross-correlation shows that the input is not totally random due to the points significantly higher than the confidence bounds. Similar conclusions can be obtained from the bottom plot, which describes the correlation between the data input and output. This nonrandom nature of the input signal can partially explain the output model behavior.
0
2
4
6
8
10
12
14
C
u,u
Correlation
Fig. 7.
Cu2,u2
1
parameter bias is avoided. The identification procedure is based on an initial frequency-domain analysis step. The objective is to gain some useful information about the system behavior so that one can define a proper structure for the Hammerstein model.
0.5
0
−0.5
0
2
4
6
8
10
12
14
The results showed that the identified model is able to represent the general nonlinear behavior. However, if more detailed nonlinear information is necessary, the user should consider another more complete model formulation. In spite of this, the proposed model can be used in several practical applications, including the attitude control loop design and analysis.
0.4 C
u,y
Correlation
0.3
Cu2,y2
0.2 0.1 0 −0.1
0
Fig. 6.
2
4
6
8
10
12
Representation of the Irregular Odd Multisine Signal.
14
REFERENCES
Input correlation analysis (95 % confidence limits).
4.3 Validation The next step is the model validation where another data set was used to evaluate the model adequacy. This is presented in Fig. 7. It is easy to see that the model presented a similar result when compared to Fig. 5. However, the model response is still quite different from the data set what indicates that a richer model can be necessary in some contexts. In spite of the model can still present unmodeled dynamics, one can note that its behavior is useful for many applications. Not only a raw nonlinear analysis can be performed by using this model, but it can be used in simulated control loops so that proper controllers can be designed. For instance, this model can be very useful for the launch vehicle attitude control. A common practice is to adopt a simplified linear model, however this can just provide approximated results. For a better comprehension of the system, the proposed model can be more adequate. 5. CONCLUSION This paper presented the model identification of an aerospace electrohydraulic servovalve. A Hammerstein model, where a static nonlinearity is followed by a linear system, is used. A polynomial NARMAX representation is considered and the model coefficients are evaluated through an extended least squares procedure, so that the 463
Chen, S. and Billings, S. (1989). Representations of nonlinear systems: the NARMAX model. Intern. Journ. Control, 49(3), 1013–1032. Jelali, M. and Kroll, A. (2003). Hydraulic servo-systems. Springer. Ljung, L. (1987). System identification - theory for the user. Prentice Hall. Merritt, H.E. (1967). Hydraulic control systems. John Willey & Sons. Pintelon, R. and Schoukens, J. (2001). System identification - a frequency domain approach. IEEE Press. Vanhoenacker, K., Dobrowiecki, T., and Schoukens, J. (2001). Design of multisine excitations to characterize the nonlinear distortions during FRF-measurements. IEEE Trans. Inst. and Measurement, 50(5), 1097–1102.
Identification of a Hammerstein model for an aerospace electrohydraulic servovalve Alexandro G. Brito ∗ Waldemar C. Leite Filho ∗ Elder M. Hemerly ∗∗ ∗
Institute of Aeronautics and Space, S˜ ao Jos´e dos Campos, SP, CEP:12220-904, Brazil (e-mails: [email protected] / [email protected]). ∗∗ Technological Institute of Aeronautics, S˜ ao Jos´e dos Campos, SP, CEP:12220-900, Brazil (e-mail: [email protected]) Abstract: This paper discusses the black-box identification of a real aerospace electrohydraulic actuator. A Hammerstein model structure is adopted due to its simplicity and capacity to represent a wide set of nonlinear phenomena. The identification strategy begins with a frequencydomain analysis so that a polynomial NARX model structure can be defined. Then, an extended least squares procedure is applied to obtain the coefficient values. Keywords: NARX/NARMAX models; FRF’s; Multisine excitation; Aerospace actuators. 1. INTRODUCTION
2. NARMAX MODELS
The electrohydraulic servovalves are widely applied to a variety of applications and the knowledge of its dynamics is essential to the control system design. Several works cover the aspects about these devices, including their internal moving parts, the motion equations and the most important nonlinearities [Merritt (1967); Jelali and Kroll (2003)]. Servovalves are especially useful as the main actuator element of aerospace systems; particularly in launchers, they operate in the thrust vector control of solid propelled vehicle or in the gimbaled chamber for liquid propelled rockets. Since their influence over the aerospace vehicle control loops, a proper modeling of the actuator linear and nonlinear behaviors is of primary importance to the launcher performance.
The NARMAX (Nonlinear AutoRegressive with Moving Average and eXogeneous inputs) models proposed by Leontaritis & Billings (1985) is a very general nonlinear representation, which includes the Volterra and Wiener series, the Bilinear, output-affine, rational and other model classes [Chen and Billings (1989)]. To define such NARMAX models, consider u(k) and y(k) the input and output signals of a general nonlinear system, respectively, each belonging to the real Banach spaces U and Y. Assume Uk = {(u(1), . . . , u(k)) : u(i) ∈ U, i = 1, . . . , k}, and the function response that describes the input-output dynamical behavior fx0 : Uk → Y, beginning from an initial state x0 . If the following conditions:
In this paper, a full discussion about the electrohydraulic actuator black-box modeling is performed. As an example, the servovalve used in the attitude control of the VLS (the Brazilian launch vehicle) is considered. The methodology is based on an initial frequency-domain analysis so that the model structure can be inferred. Then, a discrete Hammerstein formulation based on the NARMAX models is fitted. Finally, a model validation step is performed. The paper is organized as follows: the section 2 discusses the definitions associated to the NARMAX model structure selection and identification. In section 3, a whole discussion about the mathematical modeling of the electrohydraulic actuator is presented. Section 4 elaborates all the concepts used for the model identification, including the structure selection, the numerical procedure and the validation. Finally, some conclusions are presented in Section 5.
i) fx0 is a finitely realizable function response, and; ii) a linearized model exists when the system is operated in the neighborhood of an equilibrium point, hold, then there exists a nonlinear function F (·) so that y(k) = F (y(k − 1), . . . , y(k − ny ), u(k − 1), . . . , u(k − nu )). (1) Despite of the wide range of structures that the nonlinear functional F (·) can assume, only the input-output polynomial NARMAX models are discussed in this paper, since that the aim is to obtain a linear-in-the-parameter representation. A polynomial NARMAX model can be constructed by adopting the following format for the function F (·) y(k) =
y ,nu d X m nX X
m=0 p=0 n1 ,nm p Y
y(k − ni )
×
⋆ We would like to thank SIA/FINEP and FUNDEP by the financial support.
978-3-902823-46-5/2013 © IFAC
459
cp,m−p (n1 , · · · , nm ) ×
i=1
m Y
u(k − ni )
(2)
i=p+1
where d is the polynomial degree and 10.3182/20130902-5-DE-2040.00119
2013 IFAC ACA September 2-6, 2013. Würzburg, Germany
ny ,nu
X
≡
n1 ,nm
ny ny X X
n1 =1 n2 =1
···
nu X
.
(3)
nm =1
It is easy to see that the model (2) is linear in the parameters cp,m−p , thus a Least Squares identification procedure can be used. Initially, the polynomial is written as y(k) = Ψ(k − 1, :)θ + ξ(k) nθ X θi ψi (k − 1) + ξ(k) =
(4)
are functions of valve displacement xv which is related to the input voltage. The mechanical actuator mounting can introduce nonlinearities in this displacement, mainly deadzones. The load flow as a function of valve position and load pressure is given by the nonlinear relation s 1 xv QL = Cd wxv (PS − PL ) (6) ρ |xv | where PS is the supply pressure and w ≡ ∂A/∂xv ([w] = m2 /m) is the area gradient of the valve.
i=1
where the regressor matrix Ψ is constructed with all the monomial of (2). By using N measurement samples, (4) becomes a matrix equation. Due to the measurement noise, the parameters θi cannot be calculated through the ordinary least squares procedure, otherwise these parameters can be biased. Instead, other alternatives such as the Instrumental Variables or the Extended Least Squares should be considered [Ljung (1987)]. All the following identification results are obtained by adopting this last strategy. 3. ELECTROHYDRAULIC ACTUATORS The EHA consists of a powerful component in control systems that blends the versatility of electrical components with hydraulic actuators performance at high power levels. In this way, several dynamic effects – including the nonlinear one – are present and one should pay attention to some details in the modeling procedure. Fig. 1 presents a simplified block diagram for this system. The servovalve stage converts the input voltage (Vi ), the supply pressure (PS ) and the load pressure (PL ) in a load flow (QL ) to be applied to the hydraulic cylinder. This cylinder is in charge of producing the load force (fL ) necessary to provide the load stroke displacement (xp ). A thorough modeling discussion for the electrohydraulic actuator is present in Merritt (1967). In this paper, only the most important nonlinear phenomena are presented.
Fig. 1.
The fluid flowing through uncompensated valve orifices causes forces with a direction such that it tends to close the valve port. The magnitude of this force is given by F1 = 2Cd Cv A0 (P1 − P2 )cos(θ) (7) ◦ where, typically θ ≈ 69 is the jet angle of vena contracta, Cd ≈ 0.61 is the discharge coefficient, Cv ≈ 0.98 is the empirical factor called velocity coefficient and ∆P = P1 − P2 the pressure drop. By using Eq. 7 and the numerical values above, one can yield the usual form of the steadystate flow force equation F1 = 0.43w∆P xv = Kf xv . (8) Merritt (1967) comments that on larger single stage EHA, this force can exceed 20lbs. (9kgf ). The reduction of this steady-force flow force can be obtained by using twostage configuration or geometric compensating techniques. These compensations can lead to a nonlinear function between the flow force and the stroke displacement. 3.2 The cylinder dynamics The continuity equation combined with the equation of state (ρ = ρ(P, T )) leads to the expression for the load flow in the cylinder chamber dxp Vp dPL QL = Ap + (9) dt 4β dt where β is the Bulk modulus (for mineral oils and for common values for pressure and temperature, β is typically 1400 to 1600 MPa) and Ap and Vp are, respectively, the section area and the chamber volume of the cylinder. The term Vp /(4β) is equivalent to a linear hydraulic capacitance, known as CH . Since the servovalve is generally attached to the cylinder, the high pressure lines are sufficiently short and additional hydraulic capacitances due to these lines can be neglected. Finally, the load force (fL ) is given by fL = Ap PL .
Electrohydraulic actuator block diagram
3.1 The servovalve dynamics
(10)
4. EXPERIMENTAL RESULTS
The turbulent flow through an orifice occurs at high Reynolds number and is modeled by applying the Bernoulli’s equation. This analysis yields the well-known volumetric flow rate r 2 Q = Cd A0 (P1 − P2 ) (5) ρ where Cd is the discharge coefficient, A0 the orifice area, ρ the mass density of fluid and (P1 − P2 ) the pressure drop. Considering a typical four-way spool valve where the orifice areas depend on the valve geometry, their four areas 460
This paper focuses on the identification of an electrohydraulic actuator used in the attitude control of VLS, the Brazilian Launcher. The objective is to explain some nonlinear behaviors that occurred in previous flights, such as unexpected high amplitude limit-cycles. Initially, a previous analysis in the frequency domain is performed before the identification step. This procedure gives useful information about the system nonlinear behavior to the user, so that he can choose a proper Hammerstein
2013 IFAC ACA September 2-6, 2013. Würzburg, Germany
structure. In this way, the further identification process has a higher probability of success. The nonlinearity detection approach adopt in this paper is based on the frequency domain. Other technique to analyze the presence of nonlinear behavior in the frequency domain is to apply the coherence concept. However, this tool only indicates the presence of nonlinearities, but it is not able to inform which is the involved type of nonlinearities. 4.1 Nonlinearity analysis in the frequency domain For better results for the nonlinearity analysis, the user should choose a proper excitation signal so that the main nonlinear behavior can be detected. A class of excitation commonly used for this task is the odd-odd multisine signal, which has been used in several practical problems involving the identification of nonlinear systems [Pintelon and Schoukens (2001)]. Let u(k) be a discrete periodic input with fundamental frequency ω0 and a sampling frequency ωs = Nu ω0 , where Nu is the data length. Denote I = {l | l = 1, 3, 5, · · · , N < Nu /2} the set of all the odd numbers up to N and F ⊂ I a subset with some of these odd numbers previously chosen by the user. Then, the oddodd multisine signal is defined as X 2πk u(k) = Al sin lω0 + φl , φl ∼ U (−π, π). (11) ωs
lines ωe and at some even non-excited frequency lines ωp ; iii) system containing only odd nonlinearities: then, there exists useful output power at the excited frequency lines ωe and at some odd non-excited frequency lines ωi ; iv) system contains both odd and even nonlinearities: then, there exists useful output power at the excited frequency lines ωe , at some odd non-excited frequency lines ωi and at some even non-excited frequency lines ωp . Then, the nonlinearity detection is performed as follows. Firstly, an odd-odd multisinusoidal input is applied and both the input and the output signals are acquired. The next step is to calculate and analyze the Fourier transform of these signals. If the system contains useful power spectrum only at the excited lines from the input, it is linear. If the system also exhibits useful power spectra at the odd non-excited frequency lines, then it contains some odd nonlinearity such as cubic terms, saturation, etc. Finally, if it has useful power spectrum at the even non-excited frequency lines, the system contains some even nonlinearity such as quadratic terms, absolute-value function, etc. For the electrohydraulic actuator discussed in this paper, an odd-odd multisine with the characteristics presented in Table. 1 was applied. Table 1.
l∈F
where U (−π, π) denotes a uniformly distributed random number between −π and π. A sketch of this signal is presented in Fig. 2. It is important to note that, in spite of
⋆ ⋆ ⋆ ⋆ ⋆ ⋆ ⋆
Odd-odd multisine signal to be used in the actuator nonlinear analysis
sampling frequency fs = 40Hz fundamental frequency f0 = 0.005Hz start excitation frequency f1 = 0.005Hz stop excitation frequency fN = 10Hz 1000 odd frequency lines 800 effectively excited freq. lines (randomly chosen) input RMS value uRM S = 2
By applying this signal to the VLS actuator, one can obtain the output power spectra presented in Fig. 3. Some conclusions can be drawn from it. Firstly, the linear Output Power Spectrum
Fig. 2.
−20
Representation of the Odd-Odd Multisine Signal.
excit. even odd
the phase is a random number for each harmonic, y(k) is a deterministic signal since this random phase is maintained constant along the signal realization.
Anmplitude (dB)
−40 −60 −80 −100
i) pure linear system: then, there exists useful output power only at the excited frequency lines ωe ; ii) system containing only even nonlinearities: then, there exists useful output power at excited frequency 461
−120 −3 10
−2
10
−1
0
10
10
1
10
2
10
Coherence 1 0.8 0.6 g(f)
The nonlinearity detection is based on the behavior of the nonlinear systems when a multisinusoidal input is applied. This topic is richly discussed by Vanhoenacker et al. (2001) and it is only summarized herein. Consider an input-output system and odd-odd multisine signal input as described by (2). Define the set M = {m | m = N, N + 2, · · · , 2N − 1} and denote the excited frequency lines as ωe = {lω0 | l ∈ F}, the odd non-excited frequency lines as ωi = {lω0 | l ∈ (F ∪ M), F ∪ F = I} and the even nonexcited frequency lines as ωp = {lω0 | l = 0, 2, · · · , 2N }. Then, the system output satisfies one of the conditions:
0.4 0.2 0 −3 10
−2
10
−1
0
10
10
1
10
2
10
Freq (Hz)
Fig. 3.
Nonlinearity detection. Top: the output power spectrum. Bottom: The coherence for the output spectrum
contributions, which are represented as the blue crosses,
2013 IFAC ACA September 2-6, 2013. Würzburg, Germany
are in accordance with the expected one for this type of servovalve. However, it is easy to see that several nonexcited frequency lines present useful power amplitude, what indicates the presence of nonlinearities. This is especially true for the odd non-excited frequency lines – marked as red circles. Then, one can infer that the actuator has an odd nonlinear nature, according to the discussed above. In the bottom plot of Fig. 3, the coherence function for the output Fourier transform is presented. A well-known result is that the coherence function is unitary for the entire spectrum if the system is linear. Notice that this is false for the present case, mainly for the odd-nonexcited frequency lines where this value is considerably dropped. This confirms the nonlinear behavior, however is too difficult to infer the type of nonlinearities by a simple inspection of the coherence plot. 4.2 Time-domain model identification Since the conclusions drew from the last section, it is expected that the odd behavior dominates the nonlinear contributions. This information is very useful to delimitate the class of models to be considered before the application of any identification method to the data.
y = Ψθ + e
y = [ y(1) y(2) · · · y(L) ]T , e = [ e(1) e(2) · · · e(L) ]T , θ = [ a 1 a 2 b 1 b 2 b 3 ν 1 ν 2 ν 3 ν 4 ]T , and each column (regressor ) of Ψ is formed by the respective monomial of (14) for k = 1, · · · , L. The model estimation through the traditional least squares, although possible, will provide biased coefficients due to the output measurement noise. The alternative is to use the extended least squares where some regressors related to the noise e(k) are also included in Ψ so that the model coefficient identification is not affected by the correlation with the noise. In this way, the biasness is avoided. By applying this approach to the actuator data, one can obtain the coefficient values shown in Table 2. The extended least squares procedure considered 10 additional noise regressors in Ψ in an iterative procedure repeated until the coefficient convergence. Table 2.
⋆ Least Squares identification
Model Simulation
Amplitude
adopted for the both parts. Due to the odd nature of the electrohydraulic actuator (Fig. 3), an odd fifth degree polynomial is adopted: f (u(k)) = u(k) + γ3 u3 (k) + γ5 u5 (k). (12) The linear part was chosen as a third order discrete transfer function of format z(a1 z + a2 ) . (13) H(z) = 3 z − b1 z 2 − b2 z − b3 The difference equation related to the Hammerstein model in Fig. 4 is given by y(k) =a1 u(k − 1) + ν1 u3 (k − 1) + ν2 u5 (k − 1) a2 u(k − 2) + ν3 u3 (k − 2) + ν4 u5 (k − 2) b1 y(k − 1) + b2 y(k − 2) + b3 y(k − 3), (14) where the coefficients ν are formed by the product between the terms a and γ of the respective monomial.
462
est
0
95
96
97
98 Time (s)
99
100
101
0.3 Cξ,y
0.2
Cξ2,y2
0.1 0 −0.1
0
2
4
6
8
10
12
14
1.5 Cξ,ξ
1
Cξ2,ξ2
0.5 0 −0.5
Eq. (14) is linear in the parameters, hence a least squares identification procedure can be applied. By considering L data points and an output measurement noise e(k), (14) can be written as
y y
0.5
−0.5
Correlation
A hammerstein model.
−1.8096 × 10−3 1.8155 × 10−5 3.3126 × 10−4 −6.0943 × 10−6
ν1 ν2 ν3 ν4
The identified model was initially simulated with the same data set to study the estimation error behavior. The results are summarized in Fig. 5. In the top plot, the model output response is compared to the real data. Notice that both data present similar behavior. The other two plots represent, respectively, the correlation between the output and the identification residues (middle), and the residue cross-correlation (bottom). One can note a low residue cross-correlation, but still higher than the confidence limits. The same can be seen for the correlation between the model output and the identification residues.
Correlation
Fig. 4.
The model coefficients for the VLS electrohydraulic actuator.
9.3872 × 10−2 −7.2915 × 10−3 6.2472 × 10−1 1.3682 × 10−1 −1.0478 × 10−1
a1 a2 b1 b2 b3
⋆ Model structure selection A Hammerstein model is characterized by a static nonlinearity followed by a discrete linear part, as showed in Fig. 4. For this paper, a polynomial representation is
(15)
where
0
2
4
6
8
10
12
14
Fig. 5. Model response and residue analysis (95 % confidence limits). Due to the results above, one could conclude that the model is not adequate because the correlation levels are
2013 IFAC ACA September 2-6, 2013. Würzburg, Germany
Model Validation
1.5
Amplitude
est
0 −0.5 97
97.5
98
98.5
99 Time (s)
99.5
100
100.5
101
Correlation
0.3 Cξ,y
0.2
Cξ2,y2
0.1 0 −0.1
0
2
4
6
8
10
12
14
1.5 Cξ,ξ
1
Cξ2,ξ2
0.5 0 −0.5
Signal analysis
y y
0.5
Correlation
higher than the confidence bounds. A well-known identification result is that if the residue correlation function is significantly out of confidence bounds, this means that unmodeled dynamics still persist on the data. However, these conclusions can be only drawn if the input signal has a significant random aspect. To gain more information about the behavior in Fig. 5, an input correlation analysis was performed. The results are presented in Fig. 6. In the top plot, the input cross-correlation shows that the input is not totally random due to the points significantly higher than the confidence bounds. Similar conclusions can be obtained from the bottom plot, which describes the correlation between the data input and output. This nonrandom nature of the input signal can partially explain the output model behavior.
0
2
4
6
8
10
12
14
C
u,u
Correlation
Fig. 7.
Cu2,u2
1
parameter bias is avoided. The identification procedure is based on an initial frequency-domain analysis step. The objective is to gain some useful information about the system behavior so that one can define a proper structure for the Hammerstein model.
0.5
0
−0.5
0
2
4
6
8
10
12
14
The results showed that the identified model is able to represent the general nonlinear behavior. However, if more detailed nonlinear information is necessary, the user should consider another more complete model formulation. In spite of this, the proposed model can be used in several practical applications, including the attitude control loop design and analysis.
0.4 C
u,y
Correlation
0.3
Cu2,y2
0.2 0.1 0 −0.1
0
Fig. 6.
2
4
6
8
10
12
Representation of the Irregular Odd Multisine Signal.
14
REFERENCES
Input correlation analysis (95 % confidence limits).
4.3 Validation The next step is the model validation where another data set was used to evaluate the model adequacy. This is presented in Fig. 7. It is easy to see that the model presented a similar result when compared to Fig. 5. However, the model response is still quite different from the data set what indicates that a richer model can be necessary in some contexts. In spite of the model can still present unmodeled dynamics, one can note that its behavior is useful for many applications. Not only a raw nonlinear analysis can be performed by using this model, but it can be used in simulated control loops so that proper controllers can be designed. For instance, this model can be very useful for the launch vehicle attitude control. A common practice is to adopt a simplified linear model, however this can just provide approximated results. For a better comprehension of the system, the proposed model can be more adequate. 5. CONCLUSION This paper presented the model identification of an aerospace electrohydraulic servovalve. A Hammerstein model, where a static nonlinearity is followed by a linear system, is used. A polynomial NARMAX representation is considered and the model coefficients are evaluated through an extended least squares procedure, so that the 463
Chen, S. and Billings, S. (1989). Representations of nonlinear systems: the NARMAX model. Intern. Journ. Control, 49(3), 1013–1032. Jelali, M. and Kroll, A. (2003). Hydraulic servo-systems. Springer. Ljung, L. (1987). System identification - theory for the user. Prentice Hall. Merritt, H.E. (1967). Hydraulic control systems. John Willey & Sons. Pintelon, R. and Schoukens, J. (2001). System identification - a frequency domain approach. IEEE Press. Vanhoenacker, K., Dobrowiecki, T., and Schoukens, J. (2001). Design of multisine excitations to characterize the nonlinear distortions during FRF-measurements. IEEE Trans. Inst. and Measurement, 50(5), 1097–1102.