- Email: [email protected]
Application of Differential Evolution in system identification of a servo-hydraulic system with a flexible load
Mechatronics 18 (2008) 513–528
Contents lists available at ScienceDirect
Mechatronics journal homepage: www.elsevier.com/locate/mechatronics
Application of Differential Evolution in system identification of a servo-hydraulic system with a flexible load Hassan Yousefi *, Heikki Handroos, Azita Soleymani Institute of Mechatronics and Virtual Engineering, Department of Mechanical Engineering, Lappeenranta University of Technology, P.O. Box 20, FIN-53851, Lappeenranta, Finland
a r t i c l e
i n f o
Article history: Received 18 December 2006 Accepted 15 March 2008
Keywords: Differential evolution System identification Servo-hydraulic system Flexible load
a b s t r a c t Electro Hydraulic Servo Systems (EHSS) are commonly used in industry. The systems are non-linear in nature and their dynamic equations have several unknown parameters. System identification is a prerequisite to analysis of a dynamic system. One of the most promising novel evolutionary algorithms for solving global optimization problems is the Differential Evolution (DE) algorithm. The DE algorithm is proposed for handling non-linear constraint functions with boundary limits of variables to find the best values for the unknown parameters of a servo-hydraulic system with a flexible load. The DE algorithm guarantees fast speed convergence and accurate solutions regardless of the initial conditions of parameters. Several tests are carried out to validate the accuracy of selected unknown parameters. The results indicate good agreements between the related simulated and measured states. Ó 2008 Elsevier Ltd. All rights reserved.
1. Introduction Electro Hydraulic Servo Systems (EHSS) are used in industry in a wide number of applications. Their ubiquity is due to their favourable size to power ratio, precise and accurate control, and the ability to apply very large forces. The dynamics of hydraulic systems are, However highly non-linear; the system may be subjected to non-smooth non-linearities due to control input saturation, directional change of valve opening, friction, and valve overlap. Aside from the non-linear nature of the hydraulic dynamics, there are many considerable model uncertainties, such as external disturbances and leakages that cannot be modelled exactly. Furthermore, the non-linear functions that describe the uncertainties may not be known. Problems involving global optimization over continuous spaces are ubiquitous throughout the scientific community. In general, the task is to optimize certain properties of a system by choosing pertinent system parameters. For convenience, the system parameters are usually represented as a vector. The standard approach to an optimization problem begins by designing an objective function that can model the problem’s objectives while incorporating any constraints. This study will only focus on optimization methods that use an objective function. In most cases, the objective function defines the optimization problem as a minimization task. For this reason, the following investigation is restricted to minimization problems. In such problems, the objective function is more accurately called a ‘‘cost” function. * Corresponding author. Tel.: +358 5 6212482; fax: +358 5 6212499. E-mail address: Yousefi@lut.fi (H. Yousefi). 0957-4158/$ - see front matter Ó 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.mechatronics.2008.03.005
When the cost function is non-linear and non-differentiable, central to every direct search method is an algorithm that generates variations in the parameter vectors. Once a variation is generated, a decision must then be made whether or not to accept the newly derived parameters. Most stand and direct search methods use the so-called greedy criterion to make this decision. Under the greedy criterion, a new parameter vector is accepted if, and only if, it reduces the value of the cost function. One of the most promising novel evolutionary algorithms for solving global optimization problems with continuous parameters is the Differential Evolution (DE) algorithm. DE was first introduced by Storn [1] and Schwefel [2]. The extensive application areas of DE algorithms are testimony to the simplicity and robustness that have fostered their widespread acceptance and rapid growth in the research community. Initially DE was mostly applied to scientific applications involving curve fitting, for example, fitting a non-linear function to photoemissions data [3]. DE enthusiasts then hybridized it with Neural Networks and Fuzzy Logic [4] to enhance and extend its performance. DE was subsequently applied to problems involving multiple criteria as a spreadsheet solver application [5]. New areas of interest emerged, such as heat transfer [6], and constraint satisfaction problems [7]. The popularity of DE continued to grow in areas such as electrical power distribution [8], and magnetics [9]. Followed by the extension of DE into the areas of environmental science [10], and linear system models [11]. Recently DE has been applied to problems involving medical science [12], multiple criteria [13], and noisy optimization problems [14]. DE as an easy and efficient evolutionary algorithm for model optimization was introduced [15]. A new field of DE application is efficient parameter
514
H. Yousefi et al. / Mechatronics 18 (2008) 513–528
estimation in bio-filter modelling [16] and recently DE was implemented to control and synchronization of chaotic systems [17]. The main aims of this study are utilizing the DE algorithm in system identification of a non-linear servo-hydraulic system with a flexible load and using several bench tests to demonstrate the reliability and robustness of the proposed method. The paper is organized as follows. Section 2 presents the dynamic equations of the system under study. More details of Differential Evolution is described in Section 3. Experimental and simulation results are discussed in Section 4 and finally conclusions are drawn in Section 5. 2. System model The servo-hydraulic system with a flexible load is comprised of a servo-valve, a hydraulic cylinder and two masses that are connected by a parallel combination of a spring and damper. The piston and mass load are connected together using a rigid shaft. A schematic diagram of the system is given in Fig. 1. 2.1. System model To be able to represent servo-valve approximation dynamics through a wider frequency range, a first or second order transfer function is used. For low frequencies (up to about 50 Hz), a first order approximation may be sufficient [18]. The relation between spool position XS and the input voltage U can be written as, Gv ðsÞ ¼
XS s1 ¼ : U s þ s2
ð1Þ
The two unknown parameters are s1 and s2. s1 is the valve gain and s1 2 is a time constant. Applying the second Newton law for each mass and considering the friction of the hydraulic cylinder results in, € P ¼ bðX_ P X_ L Þ kðX P X L Þ þ P 1 A1 P2 A2 F f ; m1 X € L ¼ bðX_ L X_ P Þ kðX L X P Þ; m2 X
ð2Þ ð3Þ
where unknown parameter b is the viscous friction coefficient and k is the spring constant, XP and XL are the positions of the mass load (piston load) and flexible load, m1 and m2 are the masses of the piston load and flexible load, and P1 and P2 are the pressures of hydraulic cylinder, respectively. Friction (Ff) in the hydraulic cylinder is taken into account as an external disturbance and the dynamics of the servo-valve is neglected. Contact between two surfaces occurs at surface asperities (microscopic roughness) [19]. Due to the tight sealing, hydraulic cylinders feature a strong dry friction effect. The behaviour of this friction force is rather complex [20,21]. Friction is usually modelled as a discontinuous static mapping between the velocity and the friction force that depends on the velocity’s
sign. It is often restricted to the Coulomb and viscous friction components. However, there are several important properties observed in systems with friction, which cannot be explained by static models only. This is basically due to the fact that friction does not have an instantaneous response to a change in velocity, i.e., it possesses internal dynamics. Examples of these complex properties are: (1) stick-slip motion, characterized by large friction at rest and low velocities, and small friction during rapid motion; (2) pre-sliding displacement, which shows that the friction behaves like a spring when the applied force is less than the static friction break-away force; (3) friction lag, which means that there is a hysteresis that characterizes the relationship between the friction and velocity. All these static and dynamic characteristics of the friction are captured by the analytic model of friction dynamics proposed in [22], which is called the LuGre model and defined by dz jX_ P j z; ð4Þ ¼ X_ P dt gðX_ P Þ X_ 2 1 P gðX_ P Þ ¼ ; ð5Þ F c þ ðF s F c Þe vs r0 dz ð6Þ þ kv X_ P ; F f ¼ r0 z þ r1 dt where X_ P (m/s) is the piston velocity, and Ff (N) is the friction force described by the linear combination of z, dz/dt and the viscous friction. Eq. (6) represents the dynamics of the friction. The parameter z in friction equations is an internal state, not the system’s state, so the parameter is only used to calculate the friction of a hydraulic cylinder. The function gðX_ P Þ describes part of the ‘‘steady-state” characteristics of the model for constant velocity motions: vs is the Stribeck velocity, Fs is the static friction, Fc is Coloumb friction, and kv is the viscous friction. Thus, the complete friction model is characterized by four static parameters and two dynamic parameters, the stiffness coefficient (r0) and damping coefficient (r1). The friction parameters are difficult to estimate since they appear non-linearly in the model [22]. To find the proper dynamic equation for the friction in a hydraulic cylinder six unknown parameters (vs, Fs, Fc, kv, r0, r1) must be identified. Applying the flow continuity equations for the servo-valve in volume between the orifices and their outlets yields the pressures rates [22],
Table 1 Known parameters m1 = 210 kg A1 = 8.04 104 m2 m01 = 2.13 104 m3 k = 42,950 N/m Pt = 0.9 MPa
Fig. 1. Schematic diagram of system.
m2 = 80 kg A2 = 4.24 104 m2 m02 = 1.07 104 m3 L=1m Ps = 14 MPa
515
H. Yousefi et al. / Mechatronics 18 (2008) 513–528
Fig. 2. Schematic diagram of Differential Evolution algorithm.
Fig. 3. Schematic procedure of system identification.
8 < P_ 1 ¼ be1 ðQ 1 A1 X_ P þ Q Li Q Le1 Þ; V1 : P_ 2 ¼ be2 ðQ 2 þ A2 X_ P Q Li Q Le2 Þ; V2
ð7Þ
where QLi is internal leakage flow. QLe1 and QLe2 are the external leakage flows of the valve ports, A1 and A2 are the areas of the piston, V1 and V2 are the volume between each side of the piston and the servo-valve, and parameters be1 and be2 are the bulk modulus of the hydraulic fluids in each side of the piston, respectively. As provided the flow is laminar, the governing equations of internal and
Table 2 Parameters for DE Symbol
Parameter
Value
D NP F CR
Number of parameters Number of population Scale factor Crossover control constant
16 90 0.8 0.7
516
H. Yousefi et al. / Mechatronics 18 (2008) 513–528
Fig. 4. Cost of system against the number of iterations.
Table 3 The values of unknown parameters s1 = 470
s2 = 0.490 1/s
b = 418.335 Ns/m CS = 2.36 105 m3/(sVPa1/2) Li = 1724 1013 a2 = 124.368 r0 = 97525.459 m/s Kv = 376.613 Ns/m Fs = 7485.084 N
Le1 = 2.401 1012 Le2 = 1.029 1013 a1 = 0.1972 a3 = 26.814 r1 = 385.167 Ns/m Fc = 247.804 N Vs = 0.026318 m/s
external leakage flows (i.e., leakage from one chamber to the other) are presented, 8 > < Q Li ¼ Li ðP2 P 1 Þ; ð8Þ Q Le1 ¼ Le1 ðP 1 pt Þ; > : Q Le2 ¼ Le2 ðP 2 pt Þ:
The three unknown parameters are Li, Le1 and Le2. Where Li is the laminar internal leakage flow coefficient, and Le1 and Le2 are laminar external leakage coefficients. There are lots of empirical formulae based on direct measurements for the effective bulk modulus, including the effects of entrained air and mechanical compliance. A commonly used equation for effective bulk modulus, be in a hydraulic cylinder is [23], P ð9Þ þ a3 ; be ¼ a1 Emax log a2 Pmax where the variable P is the pressure in each side of the valve and the constant parameters are Emax = 1.8e9 (Pa) and Pmax = 2.8 (MPa). The three unknown parameters, a1, a2 and a3, are constant and they must be identified for any experimental setup. The volumes between the valve and each side of the piston are calculated as,
Fig. 5. An estimated parameter (damping ratio) against the number of iterations.
H. Yousefi et al. / Mechatronics 18 (2008) 513–528
V 1 ¼ A1 X P þ V 01 ; V 2 ¼ A2 ðL X P Þ þ V 02 ;
ð10Þ
where V01 and V02 are the pipeline volumes at ports one and two, respectively. The non-linear equations of valve flows are usually described by the orifice equations with a linear relationship between the valve spool position, XS and the flow area (critical centre) as follows [23]: ( pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi C S X S ps P1 U P 0; ð11Þ Q1 ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi C S X S P 1 pt U < 0; ( pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi C S X S P 2 pt U P 0; ð12Þ Q2 ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi C S X S ps P2 U < 0; where ps and pt are the supply and tank pressures, respectively, and the unknown parameter CS is a parameter including discharge coefficient and fluid density. The above equations of flow are non-linear. Table 1 shows the known parameters of the system, the other parameters, 16 in total, are unknown.
517
The non-linearity of the constitutive equations as well as the sensitivity of the system’s parameters to the sign of the voltage fed to the valve make the system identification of these kinds of systems very complicated [24]. The Differential Evolution algorithm is used to estimate the unknown parameters of the system.
3. Differential Evolution A DE algorithm is a stochastic parallel direct search optimization method that is fast and reasonably robust. Differential Evolution is capable of handling non-differentiable, non-linear, and multimodal objective functions [25]. In a population of potential solutions within an n-dimensional search space, a fixed number of vectors are randomly initialized (initialization), and then evolved over time to explore the search space to locate the minima of the objective function.
Fig. 6. Mass load position, experimental (solid line) and simulation (dash line) under a unit step excitation.
Fig. 7. Flexible load position, experimental (solid line) and simulation (dash line) under a unit step excitation.
518
H. Yousefi et al. / Mechatronics 18 (2008) 513–528
New vectors are generated at each iteration, called a generation by the combination of vectors randomly chosen from the current population (mutation). The produced vectors are then mixed with a predetermined target vector. This operation is called recombination and produces the trial vector. The trial vector is accepted for the next generation if, and only if, it yields a reduction in the value of the cost function. This last operator is referred to as a selection [25]. There are four essential versions for DE that differs only in the manner in which new solutions are generated. In the study, a classic version (DE/rand/1/bin) is used. In this shorthand notation, the first term after ‘‘DE” specifies how the base vector is chosen and rand means that the base vectors are randomly chosen. The number that follows indicates how many vector differences are used. A DE version that uses uniform crossover is appended with the additional term ‘‘bin” for ‘‘binominal” distribution [26].
The DE optimization process is conducted by means of the following operations in the rest of Section 3. 3.1. Initialization In order to establish a starting point for the optimization process, an initial population must be created. Typically, each decision parameter in every vector of the initial population is assigned a randomly chosen value from within its feasible bounds: ðloÞ
ðhiÞ
yj;i;G¼0 ¼ yj;i þ rand½0; 1 ðyj ðhiÞ yj
ðloÞ
yj Þ;
ðloÞ yj
ð13Þ
where and are the upper and lower bound to the jth decision parameter, respectively. After the initial population has been created, it evolves through the operations of mutation, crossover and selection.
Fig. 8. Mass load acceleration, experimental (solid line) and simulation (dash line) under a unit step excitation.
Fig. 9. First pressure, experimental (solid line) and simulation (dash line) under a unit step excitation.
H. Yousefi et al. / Mechatronics 18 (2008) 513–528
519
3.2. Mutation
3.3. Crossover
Mutation generally refers to an operation that adds a uniform random variable to one or more vector parameters. The DE relies upon the population itself to supply increments of the appropriate magnitude and orientation. At every generation G, each vector in the population has to serve once as a target vector. For each target vector yj,i,G, a mutant vector mj,i is generated according to
In order to increase the diversity of the generated vectors, a crossover operation is introduced. Finally, the trial vector uj,i,G+1 is created by mixing the parameter of the parent vector yj,i,G and the mutant vector mj,i in the following form: ( mj;i If ðrandj ½0; 1Þ CR _ j ¼ jr Þ; ð15Þ uj;i;Gþ1 ¼ yj;i;G otherwise;
mj;i ¼ yj;r3 ;G þ F ðyj;r1 ;G yj;r2 ;G Þ;
ð14Þ
where yj;r1 ;G ; yj;r2 ;G and yj;r3 ;G are randomly chosen vectors from the set {1, . . . ,Np), Np is the population size, F is a user-defined constant (also known as the mutation scaling factor), which is typically chosen from the range (0, 1].
where randj[0, 1) generates a new uniformly distributed random number for each value of the j within the range of [0, 1). CR is known as a crossover rate constant and is a user-defined parameter within the range [0, 1]. The index jr is randomly chosen from the set {1, . . ., D}, which is used to ensure that uj,i,G+1 gets at least one
Fig. 10. Second pressure, experimental (solid line) and simulation (dash line) under a unit step excitation.
Fig. 11. Mass load position, experimental (solid line) and simulation (dash line) under a negative unit step excitation.
520
H. Yousefi et al. / Mechatronics 18 (2008) 513–528
parameter from mj,i. The index D shows the number of unknown parameters.
the next generation are as good as or better than the individuals of the current population.
3.4. Selection
3.5. Cost function
To decide whether or not the trail vector should become a member of the next generation, the trial vector uj,i,G+1 is compared to the target vector yj,i,G using a cost function value criterion. A vector with a minimum value of cost function is allowed to advance to the next generation. That is, ( ui;Gþ1 if F Co ðui;Gþ1 Þ 6 F Co ðyi;G Þ; ð16Þ yi;Gþ1 ¼ otherwise; yi;G
The cost function in Eq. (17) is defined as the summation of the square of the normalized error for each measured state of the system at each point of the test. In the study, a unit step input (+1 V) was impelled to the system for 4.113 (s). The sampling time was 1 ms, so the total number of data points is 4113. The error between the measured and simulation points can be calculated as follow: 4 113 X 5 D X X e2j þ Pðyi ; GÞ; ð17Þ F Co ðyi ; GÞ ¼
where FCo indicates the cost function and yi,G+1 is the next generation of vector yi,G. Using this selection procedure, all individuals of
where P(yi, G) is the penalty function. The penalty function and its effect are described below. Here, ej indicates the normalized error
i¼1
j¼1
i¼1
Fig. 12. Flexible load position, experimental (solid line) and simulation (dash line) under a negative unit step excitation.
Fig. 13. Mass load acceleration, experimental (solid line) and simulation (dash line) under a negative unit step excitation.
H. Yousefi et al. / Mechatronics 18 (2008) 513–528
between the measured and simulated data for each measured state. 3.6. Penalty function which miniTo obtain the unknown parameters of vector, p mizes the cost function, the parameters are subjected to the following constraints: 6p max ; min 6 p p
ð18Þ
min and p max are the minimum and maximum range limits where p . Direct search optimization methfor the unknown parameter of p ods do not usually provide a mechanism to restrict the parameters in the range defined by inequality Eq. (18); neither does the DE. However, an unconstrained optimization method can be trans-
521
formed to a constrained one using the concept of a penalty function. This function determines a penalty to be added to the value of the cost function whenever any parameter exceeds the range limits. The penalty function selected in this case is given by the following equation: 8 2 for ai aiðminÞ ; > < 20ððaiðminÞ ai Þ=aiðminÞ Þ for aiðminÞ 6 ai 6 aiðmaxÞ ; Pðai ; GÞ ¼ 0 > : 20ððaiðmaxÞ ai Þ=aiðmaxÞ Þ2 for ai aiðmaxÞ : ð19Þ In the penalty equation for each element of a parent or child vector (ai), the maximum and minimum are ai(min) and ai(max), respectively. The penalty function adds a large value to the cost function only for each out of range element of child vector.
Fig. 14. First pressure, experimental (solid line) and simulation (dash line) under a negative unit step excitation.
Fig. 15. Second pressure, experimental (solid line) and simulation (dash line) under a negative unit step excitation.
522
H. Yousefi et al. / Mechatronics 18 (2008) 513–528
3.7. Normalization The input and outputs of the system have different units with large differences in their values, so all the input and output values are normalized by pre-processing so that they fall in the interval [1, 1] using the following equation: qn ¼ 2
q minðqÞ 1; maxðqÞ minðqÞ
other hand, constrains can eliminate local minima that might otherwise trap vectors, thereby reducing the chance that DE will prematurely converge [26]. The next section describes the application of DE for system identification of a servo-hydraulic system with a flexible load. 4. Experimental results
ð20Þ
where qn is the normalized form of the vector q. The normalization process gives the same weight and chance for all the input and outputs of the system to be involved in system identification. One of the most important issues when working with DE is the questions of constraints. Constraints typically make optimization harder for DE because they can create forbidden regions on the objective function landscape that restrict the free movement of the vectors. Depending on how they are handled, constraints can often divide the search space into several disjoint islands. On the
The servo-hydraulic valve under study was a Bosch Rexroth servo solenoid valve with On-Board Electronics (OBE) and a nominal flow rate of 0.00067 m3/s. The dimensions of the cylinder were 22 32 1000(r1 r2 h) millimetres. An accumulator with an effective volume of 0.015 m3 was used. Structural damping was neglected because it was not important in the model. A total 7 states of the system were measurable but due to a lack € P ; P1 and P2 were diof sensors the following 5 states; X P ; X L ; X rectly measured by position, accelerometer and pressure sensors, respectively. The sampling frequency was 1000 Hz. The data
Fig. 16. Mass load position, experimental (solid line) and simulation (dash line) under a step excitation (+8 V).
Fig. 17. Flexible load position, experimental (solid line) and simulation (dash line) under a step excitation (+8 V).
H. Yousefi et al. / Mechatronics 18 (2008) 513–528
acquisition system was a commercially available Digital Signal Processor (DSP). The input voltage was fed to the servo-hydraulic valve using a SD113 I/O card. The control program was the C/C++ language program which was generated by RTW [27]. First order lowpass filters with a cutoff frequency of 50 HZ were used to remove sensor noise. In this section, the DE algorithm is first used to identify the unknown parameters of the system and then several tests are utilized to validate the accuracy of the estimated parameters. A positive unit step response in extension movement is used to identify the unknown parameters, and then the unknown parameters of the system are accurately determined using Differential Evolution. Several tests (negative unit step, different magnitude step, ramp, and noisy sine inputs) are used to determine whether the estimated parameters are correctly chosen. Finally, the effects of friction in extension and retraction are presented (see Fig. 2).
523
4.1. Positive unit step The positive unit step response of the system in extension movement is reported in the study. Fig. 3 shows the schematic procedure of system identification to tune the unknown parameters of the system. The input and five states of the system for off-line system identification were measured. The mathematical model of the system based on the parameters was found Using MATLAB/SIMULINK environment. The Differential Evolution code was written in C/C++ language and used an S-Function to find the best values of unknown parameters by changing the parameters of the mathematical model of the system to minimize the cost function. An S-function is a computer language description of a Simulink block. S-functions can be written in MATLAB, C, C++, Ada, or Fortran. C, C++, Ada, and Fortran S-functions are compiled as Mex-files using the Mex utility and they are dynamically linked into MATLAB when needed [28].
Fig. 18. Ramp input excitation.
Fig. 19. Mass load position, experimental (solid line) and simulation (dash line) under ramp excitation.
524
H. Yousefi et al. / Mechatronics 18 (2008) 513–528
The important factors in the DE algorithm are population size, NP, scaling factor, F, and crossover, CR. The iterations number is very closely related to these factors. If a small population size is chosen then the saturation problem will happen and the DE cannot find the global minimum of the defined cost function. The minimum and maximum of population size is 2 and 100 times the dimensionality of the problem, respectively, and a recommended size is 5–10 times that of the dimensionality of the problem [29], There were 16 unknown parameters in the study, so the selected size of population was 5 16 = 90. In general, F and CR affect the convergence speed and robustness of the search process. Their optimal values depend both on objective function characteristics and the population size and thus, the selection of optimal parameter values is an application dependent task [1,7]. After determining the size of the population, the two essential parameters, F and CR were monitored in the MATLAB/SINMULINK environment. Using an S-Function block in the SIMULINK environ-
ment gives the opportunity to monitor all of its input and output blocks. Here, the two parameters, F and CR are considered as constant input blocks of the DE S-function. When running SIMULINK, the values of two parameters depends on the speed of convergence can be changed. The initial parameters are shown in Table 2 [26]. Fig. 4 shows the amount of cost against the number of iterations. At earlier iterations the amount of cost is large, but the algorithm finds new child generations with lower cost. Fig. 4 shows the minimum cost for the initial population is 35.6 and it decreases to 5.6 when the number of iteration is 1484. The minimum cost approaches 1.28 at 29586 iterations and it remains constant, so 1.28 is the global minimum for the proposed cost function. The number of successful generations was 2073 times during 46020 iterations. Table 3 shows the best values for the unknown parameters. Fig. 5 indicates the variations of an estimated parameter (damping ratio) against the number of iterations. The damping ration ap-
Fig. 20. Flexible load position, experimental (solid line) and simulation (dash line) under ramp excitation.
Fig. 21. Noisy Sine excitation input.
H. Yousefi et al. / Mechatronics 18 (2008) 513–528
proaches the final value (b = 418.335 Ns/m) after the same number of iterations for minimizing the cost function and it shows the robustness of the estimated parameter. The DE ran for several times with different inputs. The speed of convergence with a constant population size for the same or different inputs is only related to the following items:
525
ages that cannot be modelled exactly; and the non-linear functions that describe them may not be known. Second, the small size of the system accumulator causes the bypass effects on the results of pressures and accelerations. The small vibrations in the measured parameters in Figs. 8–10 are related to that effect. 4.2. Negative unit step
(i) amount of scaling factor, (ii) amount of crossover, (iii) uniformity of the generated random number. Figs. 6–10 are plots of system states. It is obvious that the simulation and experimental results are close to each other and DE is a strong optimizer algorithm for finding the best parameters which constrain non-linear systems. The following reasons explain the small deviations in the simulation results compared to experimental values. First, EHSSs have large model uncertainties such as external disturbances and leak-
In this section, the input voltage was the negative unit step (1 V), so the system moves in the retraction direction. The values given in Table 3 were used for the unknown parameters. Figs. 11– 15 show good agreements between experimental and simulation results in retraction movement. It is obvious that almost all the simulation states of the system under study match the experimental results. To verify the validity of the estimated parameters, several tests with different inputs were carried out. Some of the results are discussed below.
Fig. 22. Mass load position, experimental (solid line) and simulation (dash line) under Noisy Sine excitation.
Fig. 23. Flexible load position, experimental (solid line) and simulation (dash line) under Noisy Sine excitation.
526
H. Yousefi et al. / Mechatronics 18 (2008) 513–528
4.3. Different inputs Figs. 16 and 17 show the position response of mass load and flexible load to the step input with different amplitude (8 V). The simulation predicts the states of the system very well, and its performance can be considered satisfactory. Figs. 18–20 show the ramp input and the position responses of the system. Figs. 19 and 20 show that the simulated positions of mass load and flexible load have good agreement with the experimental tests. Fig. 21 shows the plot of a noisy sine input. The position responses of the system are provided in Figs. 22 and 23. There are small deviations between the experimental and simulated positions due to the effect of the noisy input. As Figs. 22 and 23 indicate, the estimated parameters were chosen correctly. Another test for a close-loop system was done to validate the results of the system identification. To make a different input, a P-controller with position feedback of mass load was designed and performed in the simulation and experimental system. The
results are depicted in Figs. 24, 25. Fig. 24 shows the closed loop input fed to the system and Fig. 25 indicates the system response (mass load position) to the input. 4.4. Friction Friction is an important aspect of many control systems. To find the effect of friction, a constant voltage was applied to the servovalve and steady-state was reached after a certain time interval. The data obtained from the simulation is visualised in Fig. 26, where friction is plotted as a function of velocity. It can be seen that the friction displays two distinct zones both in extension and retraction. A sharp increase in friction occurs in areas close to the origin. The friction magnitude was found to be greater in retraction than in extension. The explanation relies on the particular configuration of the seals. As the cylinder is of a double action type, seals on the piston must be symmetrical in order to work in both directions. At the
Fig. 24. Input excitation using the P-controller.
Fig. 25. Mass load position, experimental (solid line) and simulation (dash line) under P-controller excitation.
H. Yousefi et al. / Mechatronics 18 (2008) 513–528
527
Fig. 26. Simulation results of friction in the hydraulic cylinder.
piston rod gland however, the wipers only work during retraction, when they have to remove any dirt from the exposed part of the rod. Naturally, friction is greater during the retraction movements [30]. The result shows that the friction effects have an important rule in servo-hydraulic systems. It can lead to tracking errors and undesired stick-slip motion, so it must be considered in system identification. 5. Conclusions In this study, the application of a DE algorithm to solve the minimal representation problem in system identification of a servohydraulic system with a flexible load has been studied. The results show that the DE algorithm accurately identified the time delay, structure and parameters of the system with a fast convergence rate. Several tests were performed to validate the accuracy and robustness of estimated parameters. The results indicate good agreements between simulated and measured states. The advantage of the DE algorithm over conventional system identification methods are its simplicity, and speed of calculation and convergence, even for identifying a number of unknown parameters while simultaneously minimizing several functions. Friction is an important aspect of many control systems. In the study the effect of friction has been studied. The results show that friction in the extension and retraction movements under the same velocities is different due to the fact that the hydraulic cylinder under study was asymmetric. References [1] Storn R, Price KV. Differential evolution – a simple and efficient adaptive scheme for global optimization over continuous spaces. Technical report, TR95-012, ICSI; 1995. [2] Schwefel HP. Evolution and optimum seeking. New York: John Wiley; 1995. [3] Cafolla AA. A new stochastic optimization strategy for quantitative analysis of core level photoemission data. Surf Sci 1998;561–565:402–4. [4] Schmitz GPJ, Aldrich C. Neuro-fuzzy modeling of chemical process systems with ellipsoidal radial basis function neural networks and genetic algorithms. Comput Chem Eng 1998;22:1001–4. [5] Bergey PK. An agent enhanced intelligent spreadsheet solver for multi-criteria decision making. In: Proceedings of AIS AMCIS 1999: Americas conference on information systems, Milwaukee, USA; 1999. p. 966–8.
[6] Babu BV, Sastry KKN. Estimation of heat transfer parameters in a trickle-bed reactor using differential evolution and orthogonal collocation. Comput Chem Eng 1999;23(3):327–39. [7] Storn R. System design by constraint adaptation and differential evolution. IEEE Trans Evolut Comput 1999;3(1):22–34. [8] Chang TT, Chang HC. An efficient approach for reducing harmonic voltage distortion in distribution systems with active power line conditioners. IEEE Trans Power Delivery 2000;15(3):990–5. [9] Stumberger G, Pahner U, Hameyer K. Optimization of radial active magnetic bearings using the finite element technique and the differential evolution algorithm. IEEE Trans Magnet 2000;36(4):1009–13. [10] Booty WG, Lam CLD, Wong IWS, Siconolfi P. Design and implementation of an environmental decision support system. Environ Model Softw 2000;16(5):453–8. [11] Cheng SL, Hwang C. Optimal approximation of linear systems by a differential evolution algorithm. IEEE Trans Syst Man Cybernet Part A 2001;31(6): 698–707. [12] Abbass HA. An evolutionary artificial neural networks approach for breast cancer diagnosis. Artif Intell Med 2002;25:265–81. [13] Babu BV, Jehan MML. Differential evolution for multi-objective optimization. In: Proceedings of 2003 conference on evolutionary computation (CEC-2003), Canberra, Australia; 2003. p. 2696–703. [14] Krink T, Filipic B, Fogel GB, Thomsen R. Noisy optimization problems. A particular challenge for differential evolution. In: Proceedings of congress on evolutionary computation, vol. 1; 2004. p. 332–9. [15] Mayer DG, Kinghorn BP, Archer AA. Differential evolution –an easy and efficient evolutionary algorithm for model optimisation? Agric Syst 2005;83:315–28. [16] Babu BV, Angira R. Modified differential evolution (MDE) for optimization of non-linear chemical processes? Comput Chem Eng 2006;30:989–1002. [17] Liu B, Wang L, Jin YH, Huang DX, Tang F. Control and synchronization of chaotic systems by differential evolution algorithm. Chaos Solitons Fract 2007;34(2):412–9. [18] Niksefat N, Sepehri N. Designing robust force control of hydraulic actuators despite system and environmental uncertainties? IEEE Control Syst Mag 2001:66–77. [19] Owen WS, Croft EA. The reduction of stick-slip friction in hydraulic actuators. IEEE/ASME Trans Mechatron 2003;8(3):362–71. [20] Lischinsky P, de Wit CC, Morel G. Friction compensation for an industrial hydraulic robot. IEEE Control Syst 1999:25–32. [21] Olsson CH, Åström KJ, de Wit CC, Lischinsky P. Friction models and friction compensation. Eur J Control 1998;4:176–95. [22] de Wit CC, Olsson CH, Åström KJ, Lischinsky P. A new model for control systems with friction. IEEE Trans Automatic Control 1995;40(3): 419–25. [23] Jelali M, Kroll A. Hydraulic servo-systems; modeling, identification and control. London: Springer; 2004. [24] Viersma TJ. An analysis, synthesis and design of hydraulic servosystems and pipelines. Amsterdam: Elsevier; 1980. [25] Weisstein EW, Vassilis P. Differential evolution. From MathWorld, a Wolfram web resource; 2006.. [26] Price K, Storn RM, Lampinen JA. Differential evolution: a practical approach to global optimization. Berlin: Springer; 2005.
528
H. Yousefi et al. / Mechatronics 18 (2008) 513–528
[27] User’s guide of real time workshop. The Mathworks Inc.; 2005. [28] SIMULINK dynamic system simulation for MATLAB. The Mathworks Inc.; 1998. [29] Price K, Storn R. Differential evolution – a simple evolution strategy for fast optimization. Dr. Dobb’s J 1997;78:18–24.
[30] Bonchis A, Ha QP, Corke P, Rye DC. Model-based friction compensation in hydraulic servo systems. In: Proceedings of the australian conference on robotics and automation, Brisbane; 1999. p. 184–9.
Contents lists available at ScienceDirect
Mechatronics journal homepage: www.elsevier.com/locate/mechatronics
Application of Differential Evolution in system identification of a servo-hydraulic system with a flexible load Hassan Yousefi *, Heikki Handroos, Azita Soleymani Institute of Mechatronics and Virtual Engineering, Department of Mechanical Engineering, Lappeenranta University of Technology, P.O. Box 20, FIN-53851, Lappeenranta, Finland
a r t i c l e
i n f o
Article history: Received 18 December 2006 Accepted 15 March 2008
Keywords: Differential evolution System identification Servo-hydraulic system Flexible load
a b s t r a c t Electro Hydraulic Servo Systems (EHSS) are commonly used in industry. The systems are non-linear in nature and their dynamic equations have several unknown parameters. System identification is a prerequisite to analysis of a dynamic system. One of the most promising novel evolutionary algorithms for solving global optimization problems is the Differential Evolution (DE) algorithm. The DE algorithm is proposed for handling non-linear constraint functions with boundary limits of variables to find the best values for the unknown parameters of a servo-hydraulic system with a flexible load. The DE algorithm guarantees fast speed convergence and accurate solutions regardless of the initial conditions of parameters. Several tests are carried out to validate the accuracy of selected unknown parameters. The results indicate good agreements between the related simulated and measured states. Ó 2008 Elsevier Ltd. All rights reserved.
1. Introduction Electro Hydraulic Servo Systems (EHSS) are used in industry in a wide number of applications. Their ubiquity is due to their favourable size to power ratio, precise and accurate control, and the ability to apply very large forces. The dynamics of hydraulic systems are, However highly non-linear; the system may be subjected to non-smooth non-linearities due to control input saturation, directional change of valve opening, friction, and valve overlap. Aside from the non-linear nature of the hydraulic dynamics, there are many considerable model uncertainties, such as external disturbances and leakages that cannot be modelled exactly. Furthermore, the non-linear functions that describe the uncertainties may not be known. Problems involving global optimization over continuous spaces are ubiquitous throughout the scientific community. In general, the task is to optimize certain properties of a system by choosing pertinent system parameters. For convenience, the system parameters are usually represented as a vector. The standard approach to an optimization problem begins by designing an objective function that can model the problem’s objectives while incorporating any constraints. This study will only focus on optimization methods that use an objective function. In most cases, the objective function defines the optimization problem as a minimization task. For this reason, the following investigation is restricted to minimization problems. In such problems, the objective function is more accurately called a ‘‘cost” function. * Corresponding author. Tel.: +358 5 6212482; fax: +358 5 6212499. E-mail address: Yousefi@lut.fi (H. Yousefi). 0957-4158/$ - see front matter Ó 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.mechatronics.2008.03.005
When the cost function is non-linear and non-differentiable, central to every direct search method is an algorithm that generates variations in the parameter vectors. Once a variation is generated, a decision must then be made whether or not to accept the newly derived parameters. Most stand and direct search methods use the so-called greedy criterion to make this decision. Under the greedy criterion, a new parameter vector is accepted if, and only if, it reduces the value of the cost function. One of the most promising novel evolutionary algorithms for solving global optimization problems with continuous parameters is the Differential Evolution (DE) algorithm. DE was first introduced by Storn [1] and Schwefel [2]. The extensive application areas of DE algorithms are testimony to the simplicity and robustness that have fostered their widespread acceptance and rapid growth in the research community. Initially DE was mostly applied to scientific applications involving curve fitting, for example, fitting a non-linear function to photoemissions data [3]. DE enthusiasts then hybridized it with Neural Networks and Fuzzy Logic [4] to enhance and extend its performance. DE was subsequently applied to problems involving multiple criteria as a spreadsheet solver application [5]. New areas of interest emerged, such as heat transfer [6], and constraint satisfaction problems [7]. The popularity of DE continued to grow in areas such as electrical power distribution [8], and magnetics [9]. Followed by the extension of DE into the areas of environmental science [10], and linear system models [11]. Recently DE has been applied to problems involving medical science [12], multiple criteria [13], and noisy optimization problems [14]. DE as an easy and efficient evolutionary algorithm for model optimization was introduced [15]. A new field of DE application is efficient parameter
514
H. Yousefi et al. / Mechatronics 18 (2008) 513–528
estimation in bio-filter modelling [16] and recently DE was implemented to control and synchronization of chaotic systems [17]. The main aims of this study are utilizing the DE algorithm in system identification of a non-linear servo-hydraulic system with a flexible load and using several bench tests to demonstrate the reliability and robustness of the proposed method. The paper is organized as follows. Section 2 presents the dynamic equations of the system under study. More details of Differential Evolution is described in Section 3. Experimental and simulation results are discussed in Section 4 and finally conclusions are drawn in Section 5. 2. System model The servo-hydraulic system with a flexible load is comprised of a servo-valve, a hydraulic cylinder and two masses that are connected by a parallel combination of a spring and damper. The piston and mass load are connected together using a rigid shaft. A schematic diagram of the system is given in Fig. 1. 2.1. System model To be able to represent servo-valve approximation dynamics through a wider frequency range, a first or second order transfer function is used. For low frequencies (up to about 50 Hz), a first order approximation may be sufficient [18]. The relation between spool position XS and the input voltage U can be written as, Gv ðsÞ ¼
XS s1 ¼ : U s þ s2
ð1Þ
The two unknown parameters are s1 and s2. s1 is the valve gain and s1 2 is a time constant. Applying the second Newton law for each mass and considering the friction of the hydraulic cylinder results in, € P ¼ bðX_ P X_ L Þ kðX P X L Þ þ P 1 A1 P2 A2 F f ; m1 X € L ¼ bðX_ L X_ P Þ kðX L X P Þ; m2 X
ð2Þ ð3Þ
where unknown parameter b is the viscous friction coefficient and k is the spring constant, XP and XL are the positions of the mass load (piston load) and flexible load, m1 and m2 are the masses of the piston load and flexible load, and P1 and P2 are the pressures of hydraulic cylinder, respectively. Friction (Ff) in the hydraulic cylinder is taken into account as an external disturbance and the dynamics of the servo-valve is neglected. Contact between two surfaces occurs at surface asperities (microscopic roughness) [19]. Due to the tight sealing, hydraulic cylinders feature a strong dry friction effect. The behaviour of this friction force is rather complex [20,21]. Friction is usually modelled as a discontinuous static mapping between the velocity and the friction force that depends on the velocity’s
sign. It is often restricted to the Coulomb and viscous friction components. However, there are several important properties observed in systems with friction, which cannot be explained by static models only. This is basically due to the fact that friction does not have an instantaneous response to a change in velocity, i.e., it possesses internal dynamics. Examples of these complex properties are: (1) stick-slip motion, characterized by large friction at rest and low velocities, and small friction during rapid motion; (2) pre-sliding displacement, which shows that the friction behaves like a spring when the applied force is less than the static friction break-away force; (3) friction lag, which means that there is a hysteresis that characterizes the relationship between the friction and velocity. All these static and dynamic characteristics of the friction are captured by the analytic model of friction dynamics proposed in [22], which is called the LuGre model and defined by dz jX_ P j z; ð4Þ ¼ X_ P dt gðX_ P Þ X_ 2 1 P gðX_ P Þ ¼ ; ð5Þ F c þ ðF s F c Þe vs r0 dz ð6Þ þ kv X_ P ; F f ¼ r0 z þ r1 dt where X_ P (m/s) is the piston velocity, and Ff (N) is the friction force described by the linear combination of z, dz/dt and the viscous friction. Eq. (6) represents the dynamics of the friction. The parameter z in friction equations is an internal state, not the system’s state, so the parameter is only used to calculate the friction of a hydraulic cylinder. The function gðX_ P Þ describes part of the ‘‘steady-state” characteristics of the model for constant velocity motions: vs is the Stribeck velocity, Fs is the static friction, Fc is Coloumb friction, and kv is the viscous friction. Thus, the complete friction model is characterized by four static parameters and two dynamic parameters, the stiffness coefficient (r0) and damping coefficient (r1). The friction parameters are difficult to estimate since they appear non-linearly in the model [22]. To find the proper dynamic equation for the friction in a hydraulic cylinder six unknown parameters (vs, Fs, Fc, kv, r0, r1) must be identified. Applying the flow continuity equations for the servo-valve in volume between the orifices and their outlets yields the pressures rates [22],
Table 1 Known parameters m1 = 210 kg A1 = 8.04 104 m2 m01 = 2.13 104 m3 k = 42,950 N/m Pt = 0.9 MPa
Fig. 1. Schematic diagram of system.
m2 = 80 kg A2 = 4.24 104 m2 m02 = 1.07 104 m3 L=1m Ps = 14 MPa
515
H. Yousefi et al. / Mechatronics 18 (2008) 513–528
Fig. 2. Schematic diagram of Differential Evolution algorithm.
Fig. 3. Schematic procedure of system identification.
8 < P_ 1 ¼ be1 ðQ 1 A1 X_ P þ Q Li Q Le1 Þ; V1 : P_ 2 ¼ be2 ðQ 2 þ A2 X_ P Q Li Q Le2 Þ; V2
ð7Þ
where QLi is internal leakage flow. QLe1 and QLe2 are the external leakage flows of the valve ports, A1 and A2 are the areas of the piston, V1 and V2 are the volume between each side of the piston and the servo-valve, and parameters be1 and be2 are the bulk modulus of the hydraulic fluids in each side of the piston, respectively. As provided the flow is laminar, the governing equations of internal and
Table 2 Parameters for DE Symbol
Parameter
Value
D NP F CR
Number of parameters Number of population Scale factor Crossover control constant
16 90 0.8 0.7
516
H. Yousefi et al. / Mechatronics 18 (2008) 513–528
Fig. 4. Cost of system against the number of iterations.
Table 3 The values of unknown parameters s1 = 470
s2 = 0.490 1/s
b = 418.335 Ns/m CS = 2.36 105 m3/(sVPa1/2) Li = 1724 1013 a2 = 124.368 r0 = 97525.459 m/s Kv = 376.613 Ns/m Fs = 7485.084 N
Le1 = 2.401 1012 Le2 = 1.029 1013 a1 = 0.1972 a3 = 26.814 r1 = 385.167 Ns/m Fc = 247.804 N Vs = 0.026318 m/s
external leakage flows (i.e., leakage from one chamber to the other) are presented, 8 > < Q Li ¼ Li ðP2 P 1 Þ; ð8Þ Q Le1 ¼ Le1 ðP 1 pt Þ; > : Q Le2 ¼ Le2 ðP 2 pt Þ:
The three unknown parameters are Li, Le1 and Le2. Where Li is the laminar internal leakage flow coefficient, and Le1 and Le2 are laminar external leakage coefficients. There are lots of empirical formulae based on direct measurements for the effective bulk modulus, including the effects of entrained air and mechanical compliance. A commonly used equation for effective bulk modulus, be in a hydraulic cylinder is [23], P ð9Þ þ a3 ; be ¼ a1 Emax log a2 Pmax where the variable P is the pressure in each side of the valve and the constant parameters are Emax = 1.8e9 (Pa) and Pmax = 2.8 (MPa). The three unknown parameters, a1, a2 and a3, are constant and they must be identified for any experimental setup. The volumes between the valve and each side of the piston are calculated as,
Fig. 5. An estimated parameter (damping ratio) against the number of iterations.
H. Yousefi et al. / Mechatronics 18 (2008) 513–528
V 1 ¼ A1 X P þ V 01 ; V 2 ¼ A2 ðL X P Þ þ V 02 ;
ð10Þ
where V01 and V02 are the pipeline volumes at ports one and two, respectively. The non-linear equations of valve flows are usually described by the orifice equations with a linear relationship between the valve spool position, XS and the flow area (critical centre) as follows [23]: ( pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi C S X S ps P1 U P 0; ð11Þ Q1 ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi C S X S P 1 pt U < 0; ( pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi C S X S P 2 pt U P 0; ð12Þ Q2 ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi C S X S ps P2 U < 0; where ps and pt are the supply and tank pressures, respectively, and the unknown parameter CS is a parameter including discharge coefficient and fluid density. The above equations of flow are non-linear. Table 1 shows the known parameters of the system, the other parameters, 16 in total, are unknown.
517
The non-linearity of the constitutive equations as well as the sensitivity of the system’s parameters to the sign of the voltage fed to the valve make the system identification of these kinds of systems very complicated [24]. The Differential Evolution algorithm is used to estimate the unknown parameters of the system.
3. Differential Evolution A DE algorithm is a stochastic parallel direct search optimization method that is fast and reasonably robust. Differential Evolution is capable of handling non-differentiable, non-linear, and multimodal objective functions [25]. In a population of potential solutions within an n-dimensional search space, a fixed number of vectors are randomly initialized (initialization), and then evolved over time to explore the search space to locate the minima of the objective function.
Fig. 6. Mass load position, experimental (solid line) and simulation (dash line) under a unit step excitation.
Fig. 7. Flexible load position, experimental (solid line) and simulation (dash line) under a unit step excitation.
518
H. Yousefi et al. / Mechatronics 18 (2008) 513–528
New vectors are generated at each iteration, called a generation by the combination of vectors randomly chosen from the current population (mutation). The produced vectors are then mixed with a predetermined target vector. This operation is called recombination and produces the trial vector. The trial vector is accepted for the next generation if, and only if, it yields a reduction in the value of the cost function. This last operator is referred to as a selection [25]. There are four essential versions for DE that differs only in the manner in which new solutions are generated. In the study, a classic version (DE/rand/1/bin) is used. In this shorthand notation, the first term after ‘‘DE” specifies how the base vector is chosen and rand means that the base vectors are randomly chosen. The number that follows indicates how many vector differences are used. A DE version that uses uniform crossover is appended with the additional term ‘‘bin” for ‘‘binominal” distribution [26].
The DE optimization process is conducted by means of the following operations in the rest of Section 3. 3.1. Initialization In order to establish a starting point for the optimization process, an initial population must be created. Typically, each decision parameter in every vector of the initial population is assigned a randomly chosen value from within its feasible bounds: ðloÞ
ðhiÞ
yj;i;G¼0 ¼ yj;i þ rand½0; 1 ðyj ðhiÞ yj
ðloÞ
yj Þ;
ðloÞ yj
ð13Þ
where and are the upper and lower bound to the jth decision parameter, respectively. After the initial population has been created, it evolves through the operations of mutation, crossover and selection.
Fig. 8. Mass load acceleration, experimental (solid line) and simulation (dash line) under a unit step excitation.
Fig. 9. First pressure, experimental (solid line) and simulation (dash line) under a unit step excitation.
H. Yousefi et al. / Mechatronics 18 (2008) 513–528
519
3.2. Mutation
3.3. Crossover
Mutation generally refers to an operation that adds a uniform random variable to one or more vector parameters. The DE relies upon the population itself to supply increments of the appropriate magnitude and orientation. At every generation G, each vector in the population has to serve once as a target vector. For each target vector yj,i,G, a mutant vector mj,i is generated according to
In order to increase the diversity of the generated vectors, a crossover operation is introduced. Finally, the trial vector uj,i,G+1 is created by mixing the parameter of the parent vector yj,i,G and the mutant vector mj,i in the following form: ( mj;i If ðrandj ½0; 1Þ CR _ j ¼ jr Þ; ð15Þ uj;i;Gþ1 ¼ yj;i;G otherwise;
mj;i ¼ yj;r3 ;G þ F ðyj;r1 ;G yj;r2 ;G Þ;
ð14Þ
where yj;r1 ;G ; yj;r2 ;G and yj;r3 ;G are randomly chosen vectors from the set {1, . . . ,Np), Np is the population size, F is a user-defined constant (also known as the mutation scaling factor), which is typically chosen from the range (0, 1].
where randj[0, 1) generates a new uniformly distributed random number for each value of the j within the range of [0, 1). CR is known as a crossover rate constant and is a user-defined parameter within the range [0, 1]. The index jr is randomly chosen from the set {1, . . ., D}, which is used to ensure that uj,i,G+1 gets at least one
Fig. 10. Second pressure, experimental (solid line) and simulation (dash line) under a unit step excitation.
Fig. 11. Mass load position, experimental (solid line) and simulation (dash line) under a negative unit step excitation.
520
H. Yousefi et al. / Mechatronics 18 (2008) 513–528
parameter from mj,i. The index D shows the number of unknown parameters.
the next generation are as good as or better than the individuals of the current population.
3.4. Selection
3.5. Cost function
To decide whether or not the trail vector should become a member of the next generation, the trial vector uj,i,G+1 is compared to the target vector yj,i,G using a cost function value criterion. A vector with a minimum value of cost function is allowed to advance to the next generation. That is, ( ui;Gþ1 if F Co ðui;Gþ1 Þ 6 F Co ðyi;G Þ; ð16Þ yi;Gþ1 ¼ otherwise; yi;G
The cost function in Eq. (17) is defined as the summation of the square of the normalized error for each measured state of the system at each point of the test. In the study, a unit step input (+1 V) was impelled to the system for 4.113 (s). The sampling time was 1 ms, so the total number of data points is 4113. The error between the measured and simulation points can be calculated as follow: 4 113 X 5 D X X e2j þ Pðyi ; GÞ; ð17Þ F Co ðyi ; GÞ ¼
where FCo indicates the cost function and yi,G+1 is the next generation of vector yi,G. Using this selection procedure, all individuals of
where P(yi, G) is the penalty function. The penalty function and its effect are described below. Here, ej indicates the normalized error
i¼1
j¼1
i¼1
Fig. 12. Flexible load position, experimental (solid line) and simulation (dash line) under a negative unit step excitation.
Fig. 13. Mass load acceleration, experimental (solid line) and simulation (dash line) under a negative unit step excitation.
H. Yousefi et al. / Mechatronics 18 (2008) 513–528
between the measured and simulated data for each measured state. 3.6. Penalty function which miniTo obtain the unknown parameters of vector, p mizes the cost function, the parameters are subjected to the following constraints: 6p max ; min 6 p p
ð18Þ
min and p max are the minimum and maximum range limits where p . Direct search optimization methfor the unknown parameter of p ods do not usually provide a mechanism to restrict the parameters in the range defined by inequality Eq. (18); neither does the DE. However, an unconstrained optimization method can be trans-
521
formed to a constrained one using the concept of a penalty function. This function determines a penalty to be added to the value of the cost function whenever any parameter exceeds the range limits. The penalty function selected in this case is given by the following equation: 8 2 for ai aiðminÞ ; > < 20ððaiðminÞ ai Þ=aiðminÞ Þ for aiðminÞ 6 ai 6 aiðmaxÞ ; Pðai ; GÞ ¼ 0 > : 20ððaiðmaxÞ ai Þ=aiðmaxÞ Þ2 for ai aiðmaxÞ : ð19Þ In the penalty equation for each element of a parent or child vector (ai), the maximum and minimum are ai(min) and ai(max), respectively. The penalty function adds a large value to the cost function only for each out of range element of child vector.
Fig. 14. First pressure, experimental (solid line) and simulation (dash line) under a negative unit step excitation.
Fig. 15. Second pressure, experimental (solid line) and simulation (dash line) under a negative unit step excitation.
522
H. Yousefi et al. / Mechatronics 18 (2008) 513–528
3.7. Normalization The input and outputs of the system have different units with large differences in their values, so all the input and output values are normalized by pre-processing so that they fall in the interval [1, 1] using the following equation: qn ¼ 2
q minðqÞ 1; maxðqÞ minðqÞ
other hand, constrains can eliminate local minima that might otherwise trap vectors, thereby reducing the chance that DE will prematurely converge [26]. The next section describes the application of DE for system identification of a servo-hydraulic system with a flexible load. 4. Experimental results
ð20Þ
where qn is the normalized form of the vector q. The normalization process gives the same weight and chance for all the input and outputs of the system to be involved in system identification. One of the most important issues when working with DE is the questions of constraints. Constraints typically make optimization harder for DE because they can create forbidden regions on the objective function landscape that restrict the free movement of the vectors. Depending on how they are handled, constraints can often divide the search space into several disjoint islands. On the
The servo-hydraulic valve under study was a Bosch Rexroth servo solenoid valve with On-Board Electronics (OBE) and a nominal flow rate of 0.00067 m3/s. The dimensions of the cylinder were 22 32 1000(r1 r2 h) millimetres. An accumulator with an effective volume of 0.015 m3 was used. Structural damping was neglected because it was not important in the model. A total 7 states of the system were measurable but due to a lack € P ; P1 and P2 were diof sensors the following 5 states; X P ; X L ; X rectly measured by position, accelerometer and pressure sensors, respectively. The sampling frequency was 1000 Hz. The data
Fig. 16. Mass load position, experimental (solid line) and simulation (dash line) under a step excitation (+8 V).
Fig. 17. Flexible load position, experimental (solid line) and simulation (dash line) under a step excitation (+8 V).
H. Yousefi et al. / Mechatronics 18 (2008) 513–528
acquisition system was a commercially available Digital Signal Processor (DSP). The input voltage was fed to the servo-hydraulic valve using a SD113 I/O card. The control program was the C/C++ language program which was generated by RTW [27]. First order lowpass filters with a cutoff frequency of 50 HZ were used to remove sensor noise. In this section, the DE algorithm is first used to identify the unknown parameters of the system and then several tests are utilized to validate the accuracy of the estimated parameters. A positive unit step response in extension movement is used to identify the unknown parameters, and then the unknown parameters of the system are accurately determined using Differential Evolution. Several tests (negative unit step, different magnitude step, ramp, and noisy sine inputs) are used to determine whether the estimated parameters are correctly chosen. Finally, the effects of friction in extension and retraction are presented (see Fig. 2).
523
4.1. Positive unit step The positive unit step response of the system in extension movement is reported in the study. Fig. 3 shows the schematic procedure of system identification to tune the unknown parameters of the system. The input and five states of the system for off-line system identification were measured. The mathematical model of the system based on the parameters was found Using MATLAB/SIMULINK environment. The Differential Evolution code was written in C/C++ language and used an S-Function to find the best values of unknown parameters by changing the parameters of the mathematical model of the system to minimize the cost function. An S-function is a computer language description of a Simulink block. S-functions can be written in MATLAB, C, C++, Ada, or Fortran. C, C++, Ada, and Fortran S-functions are compiled as Mex-files using the Mex utility and they are dynamically linked into MATLAB when needed [28].
Fig. 18. Ramp input excitation.
Fig. 19. Mass load position, experimental (solid line) and simulation (dash line) under ramp excitation.
524
H. Yousefi et al. / Mechatronics 18 (2008) 513–528
The important factors in the DE algorithm are population size, NP, scaling factor, F, and crossover, CR. The iterations number is very closely related to these factors. If a small population size is chosen then the saturation problem will happen and the DE cannot find the global minimum of the defined cost function. The minimum and maximum of population size is 2 and 100 times the dimensionality of the problem, respectively, and a recommended size is 5–10 times that of the dimensionality of the problem [29], There were 16 unknown parameters in the study, so the selected size of population was 5 16 = 90. In general, F and CR affect the convergence speed and robustness of the search process. Their optimal values depend both on objective function characteristics and the population size and thus, the selection of optimal parameter values is an application dependent task [1,7]. After determining the size of the population, the two essential parameters, F and CR were monitored in the MATLAB/SINMULINK environment. Using an S-Function block in the SIMULINK environ-
ment gives the opportunity to monitor all of its input and output blocks. Here, the two parameters, F and CR are considered as constant input blocks of the DE S-function. When running SIMULINK, the values of two parameters depends on the speed of convergence can be changed. The initial parameters are shown in Table 2 [26]. Fig. 4 shows the amount of cost against the number of iterations. At earlier iterations the amount of cost is large, but the algorithm finds new child generations with lower cost. Fig. 4 shows the minimum cost for the initial population is 35.6 and it decreases to 5.6 when the number of iteration is 1484. The minimum cost approaches 1.28 at 29586 iterations and it remains constant, so 1.28 is the global minimum for the proposed cost function. The number of successful generations was 2073 times during 46020 iterations. Table 3 shows the best values for the unknown parameters. Fig. 5 indicates the variations of an estimated parameter (damping ratio) against the number of iterations. The damping ration ap-
Fig. 20. Flexible load position, experimental (solid line) and simulation (dash line) under ramp excitation.
Fig. 21. Noisy Sine excitation input.
H. Yousefi et al. / Mechatronics 18 (2008) 513–528
proaches the final value (b = 418.335 Ns/m) after the same number of iterations for minimizing the cost function and it shows the robustness of the estimated parameter. The DE ran for several times with different inputs. The speed of convergence with a constant population size for the same or different inputs is only related to the following items:
525
ages that cannot be modelled exactly; and the non-linear functions that describe them may not be known. Second, the small size of the system accumulator causes the bypass effects on the results of pressures and accelerations. The small vibrations in the measured parameters in Figs. 8–10 are related to that effect. 4.2. Negative unit step
(i) amount of scaling factor, (ii) amount of crossover, (iii) uniformity of the generated random number. Figs. 6–10 are plots of system states. It is obvious that the simulation and experimental results are close to each other and DE is a strong optimizer algorithm for finding the best parameters which constrain non-linear systems. The following reasons explain the small deviations in the simulation results compared to experimental values. First, EHSSs have large model uncertainties such as external disturbances and leak-
In this section, the input voltage was the negative unit step (1 V), so the system moves in the retraction direction. The values given in Table 3 were used for the unknown parameters. Figs. 11– 15 show good agreements between experimental and simulation results in retraction movement. It is obvious that almost all the simulation states of the system under study match the experimental results. To verify the validity of the estimated parameters, several tests with different inputs were carried out. Some of the results are discussed below.
Fig. 22. Mass load position, experimental (solid line) and simulation (dash line) under Noisy Sine excitation.
Fig. 23. Flexible load position, experimental (solid line) and simulation (dash line) under Noisy Sine excitation.
526
H. Yousefi et al. / Mechatronics 18 (2008) 513–528
4.3. Different inputs Figs. 16 and 17 show the position response of mass load and flexible load to the step input with different amplitude (8 V). The simulation predicts the states of the system very well, and its performance can be considered satisfactory. Figs. 18–20 show the ramp input and the position responses of the system. Figs. 19 and 20 show that the simulated positions of mass load and flexible load have good agreement with the experimental tests. Fig. 21 shows the plot of a noisy sine input. The position responses of the system are provided in Figs. 22 and 23. There are small deviations between the experimental and simulated positions due to the effect of the noisy input. As Figs. 22 and 23 indicate, the estimated parameters were chosen correctly. Another test for a close-loop system was done to validate the results of the system identification. To make a different input, a P-controller with position feedback of mass load was designed and performed in the simulation and experimental system. The
results are depicted in Figs. 24, 25. Fig. 24 shows the closed loop input fed to the system and Fig. 25 indicates the system response (mass load position) to the input. 4.4. Friction Friction is an important aspect of many control systems. To find the effect of friction, a constant voltage was applied to the servovalve and steady-state was reached after a certain time interval. The data obtained from the simulation is visualised in Fig. 26, where friction is plotted as a function of velocity. It can be seen that the friction displays two distinct zones both in extension and retraction. A sharp increase in friction occurs in areas close to the origin. The friction magnitude was found to be greater in retraction than in extension. The explanation relies on the particular configuration of the seals. As the cylinder is of a double action type, seals on the piston must be symmetrical in order to work in both directions. At the
Fig. 24. Input excitation using the P-controller.
Fig. 25. Mass load position, experimental (solid line) and simulation (dash line) under P-controller excitation.
H. Yousefi et al. / Mechatronics 18 (2008) 513–528
527
Fig. 26. Simulation results of friction in the hydraulic cylinder.
piston rod gland however, the wipers only work during retraction, when they have to remove any dirt from the exposed part of the rod. Naturally, friction is greater during the retraction movements [30]. The result shows that the friction effects have an important rule in servo-hydraulic systems. It can lead to tracking errors and undesired stick-slip motion, so it must be considered in system identification. 5. Conclusions In this study, the application of a DE algorithm to solve the minimal representation problem in system identification of a servohydraulic system with a flexible load has been studied. The results show that the DE algorithm accurately identified the time delay, structure and parameters of the system with a fast convergence rate. Several tests were performed to validate the accuracy and robustness of estimated parameters. The results indicate good agreements between simulated and measured states. The advantage of the DE algorithm over conventional system identification methods are its simplicity, and speed of calculation and convergence, even for identifying a number of unknown parameters while simultaneously minimizing several functions. Friction is an important aspect of many control systems. In the study the effect of friction has been studied. The results show that friction in the extension and retraction movements under the same velocities is different due to the fact that the hydraulic cylinder under study was asymmetric. References [1] Storn R, Price KV. Differential evolution – a simple and efficient adaptive scheme for global optimization over continuous spaces. Technical report, TR95-012, ICSI; 1995. [2] Schwefel HP. Evolution and optimum seeking. New York: John Wiley; 1995. [3] Cafolla AA. A new stochastic optimization strategy for quantitative analysis of core level photoemission data. Surf Sci 1998;561–565:402–4. [4] Schmitz GPJ, Aldrich C. Neuro-fuzzy modeling of chemical process systems with ellipsoidal radial basis function neural networks and genetic algorithms. Comput Chem Eng 1998;22:1001–4. [5] Bergey PK. An agent enhanced intelligent spreadsheet solver for multi-criteria decision making. In: Proceedings of AIS AMCIS 1999: Americas conference on information systems, Milwaukee, USA; 1999. p. 966–8.
[6] Babu BV, Sastry KKN. Estimation of heat transfer parameters in a trickle-bed reactor using differential evolution and orthogonal collocation. Comput Chem Eng 1999;23(3):327–39. [7] Storn R. System design by constraint adaptation and differential evolution. IEEE Trans Evolut Comput 1999;3(1):22–34. [8] Chang TT, Chang HC. An efficient approach for reducing harmonic voltage distortion in distribution systems with active power line conditioners. IEEE Trans Power Delivery 2000;15(3):990–5. [9] Stumberger G, Pahner U, Hameyer K. Optimization of radial active magnetic bearings using the finite element technique and the differential evolution algorithm. IEEE Trans Magnet 2000;36(4):1009–13. [10] Booty WG, Lam CLD, Wong IWS, Siconolfi P. Design and implementation of an environmental decision support system. Environ Model Softw 2000;16(5):453–8. [11] Cheng SL, Hwang C. Optimal approximation of linear systems by a differential evolution algorithm. IEEE Trans Syst Man Cybernet Part A 2001;31(6): 698–707. [12] Abbass HA. An evolutionary artificial neural networks approach for breast cancer diagnosis. Artif Intell Med 2002;25:265–81. [13] Babu BV, Jehan MML. Differential evolution for multi-objective optimization. In: Proceedings of 2003 conference on evolutionary computation (CEC-2003), Canberra, Australia; 2003. p. 2696–703. [14] Krink T, Filipic B, Fogel GB, Thomsen R. Noisy optimization problems. A particular challenge for differential evolution. In: Proceedings of congress on evolutionary computation, vol. 1; 2004. p. 332–9. [15] Mayer DG, Kinghorn BP, Archer AA. Differential evolution –an easy and efficient evolutionary algorithm for model optimisation? Agric Syst 2005;83:315–28. [16] Babu BV, Angira R. Modified differential evolution (MDE) for optimization of non-linear chemical processes? Comput Chem Eng 2006;30:989–1002. [17] Liu B, Wang L, Jin YH, Huang DX, Tang F. Control and synchronization of chaotic systems by differential evolution algorithm. Chaos Solitons Fract 2007;34(2):412–9. [18] Niksefat N, Sepehri N. Designing robust force control of hydraulic actuators despite system and environmental uncertainties? IEEE Control Syst Mag 2001:66–77. [19] Owen WS, Croft EA. The reduction of stick-slip friction in hydraulic actuators. IEEE/ASME Trans Mechatron 2003;8(3):362–71. [20] Lischinsky P, de Wit CC, Morel G. Friction compensation for an industrial hydraulic robot. IEEE Control Syst 1999:25–32. [21] Olsson CH, Åström KJ, de Wit CC, Lischinsky P. Friction models and friction compensation. Eur J Control 1998;4:176–95. [22] de Wit CC, Olsson CH, Åström KJ, Lischinsky P. A new model for control systems with friction. IEEE Trans Automatic Control 1995;40(3): 419–25. [23] Jelali M, Kroll A. Hydraulic servo-systems; modeling, identification and control. London: Springer; 2004. [24] Viersma TJ. An analysis, synthesis and design of hydraulic servosystems and pipelines. Amsterdam: Elsevier; 1980. [25] Weisstein EW, Vassilis P. Differential evolution. From MathWorld, a Wolfram web resource; 2006.
528
H. Yousefi et al. / Mechatronics 18 (2008) 513–528
[27] User’s guide of real time workshop. The Mathworks Inc.; 2005. [28] SIMULINK dynamic system simulation for MATLAB. The Mathworks Inc.; 1998. [29] Price K, Storn R. Differential evolution – a simple evolution strategy for fast optimization. Dr. Dobb’s J 1997;78:18–24.
[30] Bonchis A, Ha QP, Corke P, Rye DC. Model-based friction compensation in hydraulic servo systems. In: Proceedings of the australian conference on robotics and automation, Brisbane; 1999. p. 184–9.