Two level Differential Evolution algorithms for ARMA parameters estimation

Two Level Differential Evolution Algorithms for ARMA Parameters Estimation Momoh-Jimoh E. Salami*, Ismaila B. Tijani**, Ayodele .I. Abdullateef* Musa A. Aibinu*** ***Federal University of Technology, Minna, Nigeria maibinu@gmail.com *IMSRU, Faculty of Engineering, IIUM, Malaysia **Electronics Eng. Dept., ADMC, HCT, Abu Dhabi, UAE momoh@iium.edu.my, tijani.ismail@hct.ac.ae, isqeel@hotmail.co.uk Abstract— The problem of determining simultaneously the model order and coefficient of an Autoregressive Moving Average (ARMA) model is examined in this paper. An Evolutionary Algorithm (EA) comprising two-level Differential Evolution (DE) optimization scheme is proposed. The first level searches for the appropriate model order while the second level computes the optimal/sub-optimal corresponding parameters. The performance of the algorithm is evaluated using both simulated ARMA models and practical rotary motion system. The results of both examples show the effectiveness of the proposed algorithm over a well known conventional technique. Keywords—ARMA, Differential optimization, system identification. I. Evolution, Evolutionary INTRODUCTION Parametric modeling approach based on Autoregressive Moving Average (ARMA) technique has been a successful technique for signal prediction, filtering, system modeling and identification in almost all fields of study such as biomedical signal processing; image processing ; building and built environment industry [1], nuclear plant; communication [2] etc. It has been applied to determine an unknown system by the knowledge of the input and output data, or to predict the future values based on past output values, or for filtering purpose, as well as to find the frequency content (spectral estimation) or response of a system One of the problems associated with the use of this parametric modeling approach is the inability to obtain an appropriate method of model order estimation. Since the optimal model order is not known a priori, the traditional approaches have always been to evaluate various model orders based on final prediction error (FPE), Akaike information criterion (AIC), Minimum description length (MDL) and Hannan and Quinn (HQ) criterion [3-5]. FPE criterion method has been shown to only work favorably for simulated AR model, when used for real life data, the approach tends only to favor low order selection [3]. Kashyap, [5] has shown that AIC technique is statistically inconsistent, and the probability of error in choosing the correct order does not tend towards zero as the data length increases. Furthermore, both FPE and AIC tend to poorly estimate the order if the SNR is very high. MDL and HQ were proposed so as to counteract the over fitting c 978-1-4673-6322-8/13/$31.00 2013 IEEE tendency of AIC. These two methods have higher penalty factor for high model orders, thus favoring the selection of lower orders but the methods fail to work satisfactory for short data length [3]. Recently, the use of evolutionary algorithm to address either the parameters estimation or model order determination problem has been reported [6-7]. It is however noted that, most of the studies have been devoted to the use of the Genetic algorithm (GA) for parameters estimation [8]. In a recent study by Zaer et al., [8], a combination of GA and an iterative reduction process is proposed. Though, this approach has been shown to provide accurate estimation of the model order and parameters, the procedure is based on subjective decision making in the reduction of the model order pairs (p, q). In this study, a two-level evolutionary algorithm based on Differential Evolution (DE) is proposed. Due to reported advantages of DE over GA [9] among which are simplicity and its use of real parameter’s values, it is conjectured that, by using a DE to search for appropriate model order, and another DE to estimate optimal/sub-optimal model parameters, an effective and less mathematical intensive algorithm could be developed for ARMA parametric estimation. DE algorithm originally proposed by Storn and Price [10] is one of the recent and efficient evolutionary based optimization techniques. This is a population based algorithm like genetic algorithms using the similar operators; crossover, mutation and selection. Among the main advantages of DE are: finding the true global minimum of a multi modal search space regardless of the initial parameter values, fast convergence, and using a few control parameters. It has been reported that DE gives better and comparative performance with genetic algorithm on real-life problems [9, 11]. The rest of the paper is organized as follows. A brief description of ARMA model and problem formulation is given in Section 2. The overview of DE algorithm is presented in Section 3, which is followed by the description of the proposed two level DE algorithms in Section 4. Section 5 and 6 provide application examples, results and discussion, respectively. The paper is concluded in Section 7. 1337 II. ARMA MODEL DEVELOPMENT: PROBLEM FORMULATION ϕ ig, j=0 = L j + ( H j − L j ) The ARMA model generally involves representation of inputoutput relationship of a system by a difference equation of the form [12] p q k =1 k =0 y ( n ) = − ¦a k y ( n − k ) + ¦b k x ( n − k ) (1) where a k and b k are the model coefficients, p and q are real-valued model order for the AR and MA parts respectively. Taking the z-transform of both sides of this equation gives: −1 −q Y (z ) b o + b 1 z + ... + b q z (2) = X (z ) 1 + a 1 z −1 + ... + a p z − p where Y(z) is the system zeros of size q (i.e. order of Y(z)), and X(z) gives the system poles of size p (i.e. order of X(z)). Usually, the output y(n) is observed with an additive white Gaussian noise, w (n ) , such that: y w (n ) = y (n ) + w (n ) (3) rand(0,1) (5) Mutation is the first step in the generation of new vector known as child/trial vector, χ i . Based on common “DE/rand/1/bin” strategy, it involves random selection of three vectors with indexes r 1 , r 2 , r 3 such that r 1 ≠ r 2 ≠ r 3 . The mutant vector υ i is given as: υ i = ϕ r + F (ϕ r − ϕ r ) 3 2 (6) 1 where F is the scaling factor which is a positive real number in the range of (0,1) that controls the rate at which the population evolves. The mutation process is complemented with recombination process known as crossover process. This ensures that each parameter of the differential mutant vector is accepted into the trial vector with some probability, CR, known as crossover constant. That is, CR ∈ [0,1] is a probability value that controls the fraction of parameter values that are copied from the mutant into the trial vector as follows: For a given input-output data pair of a system with unknown order (i.e. p and q are unknown a priori), the problem of ARMA model development involves simultaneous determination of the order (p and q), and estimation of the corresponding model parameters: a k and bk . In order to of (7) is to ensure that the trial vector does not duplicate ϕ i , j , address this duo design problem, a two-level evolutionary algorithm based on DE is proposed in this study as detailed in the next Section. and to ensure that at least one parameter is altered. Comparison of the objective values of both parent vector, ϕ i , III. DE ALGORITHM OVERVIEW Given a linear time invariant state-space model of the system Similar to other evolutionary optimization process, the goal of DE is to search for a set of decision variables ϕ j , ϕj j = 1 , 2 , ..., J , ϕ ∈ ℜJ which minimizes the objective function f k : min f ( ϕ ) subject to bounds on the decision variables: L j ≤ ϕ j ≤ H j j = 1 ... J and constraints g c ( ϕ ) ≤ βc where J is dimension of the decision variables, L (4) j and H if ( rand ( 0 ,1 ) ≤ CR or j = Rand otherwise The randomly chosen index, Rand j j (7) in the selection condition and trial vector, χ i , is carried out under selection process. If the trial vector has an equal or lower objective function value than that of its target vector (parent vector), it replaces the target vector in the new generation; otherwise, the target vector retains its place in the population. The creation of the new generation vector, g + 1 , from the old generation, g , is given as ­χ g ϕ ig +1 = ® ig ¯ϕ i if f ( χ ig ) ≤ f ( ϕ ig ) (8) otherwise Once the new population is generated, the process of mutation, crossover (recombination) and selection is repeated until the pre-specified termination criterion is met. j are lower and upper bounds on the decision variables, while g c ( ϕ ) is constrain to be met by the solutions. The space spanned by the decision variables is called search space, ℜ J , while the space spanned by the objective values is known as solutions space, ℜ s . Like other EAs, DE solves this problem by using the basic concepts of initialization, mutation, and selection as detailed in [13]. Initialization involves generation of initial population for the DE process. Giving a lower and upper bound on the decision variables, the initialization of the jth parameter of ith vector for initial generation (g=0), is expressed as: 1338 ­υ i , j χi , j = ® ¯ϕ i , j IV. PROPOSED TWO LEVEL ALGORITHM FOR ARMA MODEL DEVELOPMENT The DE process described in the previous section is adopted to develop a two-level evolutionary algorithm for simultaneous determination and estimation of an ARMA model order and coefficients respectively. Figure 1 shows the block diagram of the proposed algorithm. The DE Level 2 follows the standard DE algorithm as described in Section 3 to estimate 2013 IEEE 8th Conference on Industrial Electronics and Applications (ICIEA) the parameters of a given ARMA model order by the DE level allowable maximum ( p max / q max ) and minimum ( p min / q min ) 1. model order, respectively. ii. Apart from the bounds on the decision variables, additional constraint is added to ensure the resulting model is proper, that is, q <= p . iii. A floating point transformation suggested by [14] is adopted for the mutation process. The integers values of the decision variable is transformed to floating point equivalent by dividing the selected mutation candidates with specified upper bound ( p max / q max ). Then, the mutation and crossover processes given in equations (6) and (7) are applied to generate child vector. The final decision vector is obtained by . re-transformation Figure 1 Block diagram of proposed two level DE algorithm of the resulting floating values to corresponding integer values by multiplication and truncation process. The goal of the DE level 1 is to search for set of decision ϕj iv. It is observed that, the conventional selection process given , (the parameters of the ARMA model) that in equation (8) does not discriminate situation when two minimizes the objective function given as the mean square models give rise to the same objective values. This may variables, violate the general principle of parsimony in statistics [15] error between the measured, y w and estimated output, ŷ which favors simplicity of model. To address this problem, the min 1 N f ( ϕ ) := ¦ (y w (n ) − ŷ (n )) N n =0 (9) where N is the length of given input-output data set. The optimal parameter set and the corresponding objective value are passed to the DE level 1 for the model order selection process is carried out using both the objective value and model order as follows: Modified selection process: a trial vector dominates (replaces) its parent vector if its objective values is better (less) than that of the parent vector, and in the event that both trial and parent determination. Unlike general optimization problem followed by the DE level 2, the level 1 DE process is a typical combinatorial problem involving finite discrete parameter space, and due to the use of real parameter values, DE has been shown to be effective for vectors’ objective values are equal, the model order of both candidates are examined. Then, the trial vector dominates if its model order is lesser than that of the parent vector. This is expressed as: this type of problem [14]. Hence, the general DE algorithm reported in Section 3 is modified to give the DE level 1 ϕi g +1 algorithm appropriate for this type of task, and yet meet the ­ χ ig ° = ® χ ig °ϕ g ¯ i if f (χ ig ) = f (ϕ ig elseif f (χ otherwise g i ) =f ) (ϕ ) & Γ(χ ) < Γ(ϕ ) g i g i g i (10) peculiarity of the ARMA model order constraints. The major where Γ(∗) is an additional function that computes the model modification is as follows: order . i. The initialization of population is achieved using random By implication, the algorithm tends to give preference to low integer generation function “randint” instead of “rand()” in order model without compromising the model accuracy. This equation (5). The upper, H and lower, L, bounds represent the 2013 IEEE 8th Conference on Industrial Electronics and Applications (ICIEA) 1339 is particularly useful for many practical applications where The DC motor driven rotary motion system represents a reduction of model complexity is highly desired. typical practical problem in system identification. It consists V. of servo motor driven by an amplifier and position encoder APPLICATION EXAMPLES The performance of the proposed algorithm is evaluated on two known parametric models (ARMA (3, 2) and ARMA (5, 4)) similarly used in [8], and a practical DC motor driven rotary motion system shown in Figure 2. ARMA (3, 2) given in equation (11) is an example of low order system with p = 3 and q = 2 , while ARMA (5,4) in equation (12) represents high order system with p = 5 and q = 4 . That is ARMA (3, 2) attached to the shaft as the feedback sensor. Although, an approximate lump mathematical model based on the physical laws of motion and Kirchhoff’s laws could be derived for the system as reported in [16], the actual model order and parameters of the complete system components are not always known a priori. The system is excited with PRBS (Figure 2) signal to obtain input-output data pairs required for system identification process. Given this as a system with unknown satisfies the difference equation. order and coefficients, the proposed algorithm is applied to y(n) + 0.3 y(n − 1) + 0.25 y(n − 2) + 0.5 y(n − 3) (11) = x(n) + 0.9 y(n − 1) + 0.6 y(n − 2) identify both the model order and coefficients of the system. The necessary pre-specified optimization parameters for both while ARMA (5, 4) is described by levels of the algorithm are: population size, 10*D (where D is y (n ) + 0.28 y (n − 1) + 0.619 y (n − 2) + 0.24 y (n − 3) + .. + 0.1 y ( n − 4) + 0.32 y (n − 5) = x(n ) + 0.27 y (n − 1) + . + 0.37 y ( n − 2) − 0.13 y ( n − 3) + 0.27 y (n − 4) (12) A Pseudo Random Binary Signal (PRBS) shown in Figure 3, is used to excite the dynamics of the model to obtain two sets of input-output data pair requires for the identification. One of the data set is corrupted with a white Gaussian noise to size decision variables for each level), mutation constant, F=0.75; crossover constant CR=0.25; and generation size 10 and 100 for level 1 and level 2, respectively, L=1, H=20 for the level 1, and L=-1,H=1 for the level 2. For comparative study, Minimum description length (MDL) together with prediction error model (PEM) by Ljung [17] is applied to estimate the two ARMA models in equation (11) and (12), and to estimate the plant model. The results and discussion are presented in Section VI. examine the effect of output noise. VI. RESULTS AND DISCUSSION The estimated ARMA models for the two simulated models: ARMA (3,2) and ARMA(5,4) with noise-free data are given in equations (13) and (14), respectively. Similarly, the estimated models with snr = 10 are given in equation (15) and (16) for ARMA (3, 2) and ARMA (5, 4), respectively. (1 + 0 . 3 z + 0 . 2502 z + 0 . 5002 z )y (n ) = (1 + 0 . 8999 z + 0 . 6002 z )x (n ) (1 + 0 . 313 z + 0 . 232 z + 0 . 505 z )y (n ) = (1 + 0 . 923 z + 0 . 572 z )x (n ) −1 Figure 2 DC-Motor driven rotary motion systems [16]. 1.5 Amplitude (units) −2 −1 0.5 -1 −1 5 10 15 20 25 Time (s) 30 35 40 45 50 −3 −3 −2 −2 −1 −4 −3 −2 −5 −3 −3 Figure 3: PRBS input excitation signal 1340 (14) −2 −2 −1 -0.5 (13) (1 + 0 . 22 z + 0 . 203 z + 0 . 237 z + 0 . 084 + 0 . 364 z )y (n ) 15) = (1 + 0 . 1853 z + 0 . 4033 z − 0 . 1317 z + 0 . 27 z )x (n ) (1 + 0 . 25 z + 0 . 233 z + 0 . 244 z + 0 . 104 + 0 . 36 z )y (n ) 16) = (1 + 0 . 24 z + 0 . 369 z − 0 . 19 z + 0 . 297 z )x (n ) −1 0 0 −3 −2 −1 1 -1.5 −2 −1 2013 IEEE 8th Conference on Industrial Electronics and Applications (ICIEA) −4 −4 −5 −4 As indicated in equations (13) to (16), the proposed algorithm has been able to estimate correctly the model order, and associated parameters. Average parameter error which is computed as the average of error between the actual and estimated parameters, and overall fit is used to compare the overall performance of the proposed algorithm and MDLPEM on the two given ARMA models. The fit is the percentage of the output variation that is reproduced by the model and is defined mathematically as: (p=3,q=4). The responses of the model obtained using the proposed algorithm and that from the MDL-PEM technique for the plant model development are shown in Figure 4, and the overview of the performance comparison is presented in Table II. As shown in Figure 4 and presented in Table II, up to 9% improvement in model fit is achieved with the proposed algorithm over the use of MDL-PEM approach. Measured and simulated model output 40 30 ( )= y − mean ( y ) ∗ 100 (17) where y is the measured output (plant response), and ŷ is the simulated model output. 20 position (degree) Fit 1 − y − ŷ PEM: Fit>> 69.714% Measured DE based Algorithm Fit>>77.96% 10 0 -10 -20 The comparative results of the proposed DE-based algorithm and MDL-PEM for the two ARMA models are given in Table I. The MDL approach fails to estimate correctly the model -30 -40 25 30 35 40 45 50 Time (s) Figure 4: Measures and simulated plant responses with MDL-PEM and DEbased algorithm order of the ARMA (5, 4) despite several runs using the inbuilt MATLAB identification toolbox. As shown in Table I, TABLE II. PERFORMANCE COMPARISON FOR PLANT MODEL better average parameter errors and fit are achieved by the proposed algorithm compared to the MDL-PEM. TABLE I. DE Algorithm MDL+PEM model order 3,2 3,4 fit (%) 77.76 69.7 PERFORMANCE COMPARISON FOR SIMULATED MODELS VII. CONCLUSION Similarly, the estimated plant model produced by the proposed algorithm is given as: (1 − 0 . 9218 z − 0 . 366 z + 0 . 286 z )y (n ) = (0 . 570 + 0 . 995 z + 0 . 799 z )x (n ) −1 −2 −3 −1 −2 (18) while the resulting model using MDL-PEM for the plant model identification is given as: (1 − 0 . 3426 = (1 . 845 z z −1 −1 − 0 . 345 z + 1 . 952 z −2 −2 + 0 . 30966 z + 1 . 253 z −3 −3 )y (n ) + 0 . 3362 z −4 )x (n ) (19) It is observed from equations (18) and (19) that, the proposed algorithm was able to produce lower order model (p=3,q=2) compared to the model produced by the MDL-PEM A two-level DE based optimization algorithm has been proposed in this study for ARMA parameters estimation. The algorithm provides simultaneous model order and associated coefficients estimation of the model from a given input-output data set. Performance of the algorithm has been evaluated on both simulated data and real plant model. The results and comparative study have shown the effectiveness of this proposed design approach in overcoming the challenges associated with ARMA technique. This is expected to facilitate the applications of this technique in various emerging engineering and scientific problems. Future study will be directed towards further applications of the algorithm to practical emerging problems. In addition, its extension to handle multi-input multi output systems will constitute part of the future study. Acknowledgment This research was supported by the Fundamental Research Grant Scheme (FRGS), Research Management Center, IIUM, Malaysia. FRGS12-065-0214. 2013 IEEE 8th Conference on Industrial Electronics and Applications (ICIEA) 1341 REFERENCES [1] [2] [3] [4] [5] [6] [7] [8] [9] Aibinu, A.M., Shae, A. A., Salami, M. J. E., Salami, A. F., Bamgbopa I. A., & Lawal, W. A. (2008b). Development of a new method of crack modeling and prediction algorithm. Proc. 3rd International conference on mechatronics (ICOM-08), Kuala Lumpur, Malaysia, pp 434-438. Fattah, S. A. & Zhu, W. P. (2008). An algorithm for the identication of autoregressive moving average systems from noisy observations. Proc. of Canadian Conf. on Electrical and Computer Engineering , pp.18151818. Marple, S. L. (1987). Digital spectral analysis with applications. Prentice Hall, Inc. Manolakis, D. G, Ingle, V. K. & Kogon, S. M. (2000). Statistical and adaptive signal processing. Mc. Graw Hill. Kashyap R. L. & Chellappa, R. (1981). Stochastic models for closed boundary analysis. Representation and Reconstruction. In IEEE Trans. on Information Theory, Vol. IT-27, No. 5 pp. 627-637. S. Rolf, J. Sprave, and W. Urfer, Model identification and parameter estimation of ARMA models by means of evolutionary algorithms, Comput. Intell. Financial Eng. 23 (1997), pp.237–243. M. Chiogna, C. Gaetan, and G. Masarotto, Automatic identification of seasonal transfer func-tion models by means of iterative stepwise and genetic algorithms, J. Time Series Anal. 29 (2007), pp. 37–50. Za'er S. Abo-Hammour, Othman M.K. Alsmadi, Adnan M. Al-Smadi, Maha I. Zaqout and Mohammad S. Saraireh (2012): ARMA model order and parameter estimation using genetic algorithms, Mathematical and Computer Modelling of Dynamical Systems: Methods, Tools and Applications in Engineering and Related Sciences, 18:2, 201-221 Nurhan Karaboga and Bahadir Cetinkaya, (2004). “Performance Comparison of Genetic and Differential Evolution Algorithms for 1342 [10] [11] [12] [13] [14] [15] [16] [17] [18] Digital FIR Filter Design T. Yakhno (Ed.): ADVIS 2004, LNCS 3261, pp. 482–488, 2004.Springer-Verlag Berlin Heidelberg 2004 Nurhan Karaboga and Bahadir Cetinkaya, (2004). “Performance Comparison of Genetic and Differential Evolution Algorithms for Digital FIR Filter Design T. Yakhno (Ed.): ADVIS 2004, LNCS 3261, pp. 482–488, 2004.Springer-Verlag Berlin Heidelberg 2004 Storn, R. and Price, K., (1997), “Differential Evolution – A Simple and Efficient Heuristic for Global Optimization over Continuous Spaces” Journal of Global Optimization 11: 341–359, 1997. 341 1997 Kluwer Academic Publishers. Printed in the Netherlands Vitaliy Feoktistov,(2006). “Differential Evolution:In Search Of Solutions”. 2006 Springer Science and Business Media, LLC. Proakis, J. G. & Manolakis, D. G. (2007). Digital signal processing: principles, algorithms and applications. ( 4th edn). Pearson Prentice Hall. Tijani I.B.. Flight control system with MODE based H-infinity for small scale autonomous helicopter. PhD thesis submitted to Mechatronics engineering department, IIUM,Malaysia, October 2012. Kenneth V. Price, Rainer M. Storn and Jouni A. Lampinen, 2005. Differential Evolution: A Practical Approach to Global Optimization. Springer-Verlag Berlin Heidelberg 2005, Printed in Germany. Gauch Jr., Huch G. (1993): "Prediction, Parsimony and Noise," American Scientist 81: 468-478. I.B. Tijani and Rini Akmeliawati, “Support vector regression based friction modeling and compensation in motion control system”. Eng. Appl. of Artif. Intell., vol. 25, issue 5, Elsevier, August, 2012 (ISI) Ljung L. Issue in system identification. IEEE Control systems, vol. 11, no. 1, pp. 25-29, 1991. 2013 IEEE 8th Conference on Industrial Electronics and Applications (ICIEA)