Evolutionary Parametric Identification of Dynamic Systems

This book presented techniques and experimental results which have been pursued for the purpose of evolutionary robotics. Evolutionary robotics is a new method for the automatic creation of autonomous robots. When executing tasks by autonomous robots, we can make the robot learn what to do so as to complete the task from interactions with its environment, but not manually pre-program for all situations. Many researchers have been studying the techniques for evolutionary robotics by using Evolutionary Computation (EC), such as Genetic Algorithms (GA) or Genetic Programming (GP). Their goal is to clarify the applicability of the evolutionary approach to the real-robot learning, especially, in view of the adaptive robot behavior as well as the robustness to noisy and dynamic environments. For this purpose, authors in this book explain a variety of real robots in different fields. For instance, in a multi-robot system, several robots simultaneously work to achieve a common goal via interaction; their behaviors can only emerge as a result of evolution and interaction. How to learn such behaviors is a central issue of Distributed Artificial Intelligence (DAI), which has recently attracted much attention. This book addresses the issue in the context of a multi-robot system, in which multiple robots are evolved using EC to solve a cooperative task. Since directly using EC to generate a program of complex behaviors is often very difficult, a number of extensions to basic EC are proposed in this book so as to solve these control problems of the robot


Parametric identification fundamentals
Consider a linear, time-invariant and casual dynamic system, with a single input and a single output, described by the following equation in the z-domain (Oppenheim & Schafer, 1989), (1) where X(z) and Y(z) denote the z-transforms of input and output respectively, and H(z) is a rational transfer function, with respect to the variable z, which describes the input-output dynamics.It should be noted that the selection of representing the true system in the zdomain is justified from the fact that data are always acquired in discrete time units.Due to one-to-one relationship between the z-transform and it's Laplace counterpart, it is easy to obtain a corresponding description in continuous time.The identification problem pertains to the estimation of a finitely parameterized transfer function model of a given structure, similar to that of H(z), by means of the available data set and taking under consideration the presence of noisy measurements.The estimated model must have similar properties to that of the true one, it should be able to simulate the dynamic system and, additionally, to predict future values of the output.Among a large number of ready-made models (known also as black-box models), ARMAX is widespread and has performed well in many engineering applications (Petsounis & Fassois, 2001).

The ARMAX model structure
A SISO ARMAX (na,nb,nc,nk) model has the following mathematical representation (2) where u t and y t , represent the sampled excitation and noise corrupted response signals, for time t = 1, ...,N respectively and e t is a white noise sequence with and , where and are Kronecker's delta and white noise variance respectively.N is the number of available data, q denotes the backshift operator, so that y t• q k =y t-k , and A(q), B(q), C(q) are polynomials with respect to q, having the following form (3) (4) (5) The term q -nk in ( 4) is optional and represents the delay from input to output.
In literature, the full notation for this specific model is ARMAX (na, nb, nc, nk), and it is totally described by the order of the polynomials mentioned above, the numerical values of their coefficients, the delay nk, as well as the white noise variance In Eq. ( 2) it is obvious that ARMAX consists of two transfer functions, one between input and output (6) which models the dynamics of the system, and one between noise and output (7) which models the presence of noise in the output.For a successful representation of a dynamic system, by means of ARMAX models, the stability of the above two transfer functions is required.This can be achieved by letting the roots of A(q) polynomial lie outside the unit circle with zero origin, in the complex plane (Ljung, 1999, Oppenheim & Schafer, 1989).In fact, there is an additional condition that must hold and that is the invertibility of the noise transfer function H(q) (Ljung, 1999, Box et al.,1994, Soderstrom & Stoica, 1989).For this reason, C(q) polynomial must satisfy the same requirement as A(q).

Formulation of PEM
For a given data set over the time t, it is possible to compute the output y t of an ARMAX model, at time t +1.This fact yields, for every time instant, to the formulation of one step ahead prediction-errors sequence, between the actual system's response and the one computed by the model ( 8) where p=[a i b i c i ] is the parameter vector to be estimated, for given orders na, nb, nc and delay nk, y t+1 the measured output, the model's output and the prediction error (also called model residual).The argument (1/p) denotes conditional probability (Box et al., 1994) and the hat indicates estimator/estimate.The evaluation of residuals is implemented through a scalar-valued function (see Introduction), which in general has the following form (9) Obviously, the parameter p which minimizes V N is selected as the most suitable (10) Unfortunately, V N cannot be minimized analytically due to the non-linear relationship between the model residuals êt (1/p) and the parameter vector p.This can be noted, by writing (2) in a slightly different form (11) The solution then has to be found by iterative, numerical techniques and this is the reason for the implementation of optimization algorithms within the PEM framework.

Optimization algorithms
In this section, the hybrid optimization algorithm is presented.The new method is a combination of a stochastic and a deterministic algorithm.The stochastic component belongs to the Evolutionary Algorithms (EA's) and the deterministic one to the quasi-Newton methods for optimization.

Evolution strategies
In general, EA's are methods that simulate natural evolution for the task of global optimization (Baeck, 1996).They originate in the theory of biological evolution described by Charles Darwin.In the last forty years, research has developed EA's so that nowadays they can be clearly formulated with very specific terms.Under the generic term Evolutionary Algorithms lay three categories of optimization methods.These methods are Evolution Strategies (ES), Evolutionary Programming (EP) and Genetic Algorithms (GA) and share many common features but also approximate natural evolution from different points of view.The main features of ES are the use of floating-point representation for the population and the involvement of both recombination and mutation operators in the search procedure.Additionally, a very important aspect is the deterministic nature of the selection operator.The more advanced and powerful variations are the multi-membered versions, the so-called ( + )-ES and ( , )-ES which present self-adaptation of the strategy parameters.

The quasi-Newton BFGS optimization method
Among the numerous deterministic optimization techniques, quasi-Newton methods are combining accuracy and reliability in a high level (Nocedal & Wright, 1999).They are derived from the Newton's method, which uses a quadratic approximation model of the objective function, but they require significantly less computations of the objective function during each iteration step, since they use special formulas in order to compute the Hessian matrix.The decrease of the convergence rate is negligible.The most popular quasi-Newton method is the BFGS method.This name is based on its discoverers Broyden, Fletcher, Goldfarb and Shanno (Fletcher, 1987).

Description of the hybrid algorithm
The optimization procedure presented in this paper focuses in interconnecting the advantages presented by EA's and mathematical programming techniques, and aims at combining high convergence rate with increased reliability in the search for the global optimum in real parameter optimization problems.The proposed algorithm is based on the distribution of the local and the global search for the optimum.The method consists of a super-positioned stochastic global search and an independent deterministic procedure, which is activated under conditions in specific members of the involved population.Thus, while every member of the population contributes in the global search, the local search is realized from single individuals.Similar algorithmic structures have been presented in several fully stochastic techniques that simulate biological procedures of insect societies.Such societies are distributed systems that, in spite of the simplicity of their individuals, present a highly structured social organization.As a result, such systems can accomplish complex tasks that in most cases far exceed the individual's capabilities.The corresponding algorithms use a population of individuals, which search for the optimum with simple means.The synthesis, though, of the distributed information enables the overall procedure to solve difficult optimization problems.Such algorithms were initially designed to solve combinatorial problems (Dorigo et al., 2000), but were soon extended to optimization problems with continuous parameters (Monamarche et al., 2000, Rjesh et al., 2001).A similar optimization technique presenting a hybrid structure has been already discussed in (Kanarachos , 2002), and it's based on a mechanism that realizes cooperation between the (1,1)-ES and the Steepest Descent method.The proposed methodology is based on a mechanism that aims at the cooperation between the ( + )-ES and the BFGS method.The conventional ES (Baeck, 1996, Schwefel, 1995), is based on three operators that take on the recombination, mutation and selection tasks.In order to maintain an adequate stochastic character of the new algorithm, the recombination and selection operators are retained with out alterations.The improvement is based on the substitution of the stochastic mutation operator by the BFGS method.The new deterministic mutation operator acts only on the non-privileged individuals in order to prevent loss of information from the corresponding search space regions, while three other alternatives were tested.In these, the deterministic mutation operator is activated by: • every individual of the involved population, • a number of privileged individuals, and • a number of randomly selected individuals.The above alternatives led to three types of problematic behavior.Specifically, the first alternative increased the computational cost of the algorithm without the desirable effect.
The second alternative led to premature convergence of the algorithm to local optima of the objective function, while the third generated unstable behavior that led to statistically low performance.

Efficiency of the hybrid algorithm
The efficiency of the hybrid algorithm is compared to that of the (15 +100)-ES, the (30, 0.001, 5, 100)GA, as well as the (60, 10, 100)meta-EP method, for the Fletcher & Powell test function, with twenty parameters.Progress of all algorithms is measured by base ten logarithm of the final objective function value The results are presented in Table 1.

Description of the proposed method
The proposed method for the parameter estimation of ARMAX(na,nb,nc,nk) models consists of two stages.In the first stage, Linear Least Squares are used to estimate an ARX(na,nb,nk) model of the form ( 14) based upon the observation that the nonlinear relationship between the model residuals and the parameter vector would be overcome if C(q) polynomial was monic (see ( 11)).
Considering the same loss function, as in ( 9), by expressing the ARX model in ( 14) as ( 15) with ( 16) being the regression vector and p*=[a i b i ] the parameter vector, the minimizing argument for the model in ( 14) can be found analytically, by setting the gradient of V N equal to zero, which yields The parameter vector p, as computed from ( 17), can be used as a good starting point for the values of A(q) and B(q) coefficients.In other estimation methods, like the two Stage Least Squares or the Multi-Stage Least Squares (Petsounis & Fassois, 2001) the ARX model is estimated with sufficiently high orders.In the presented method this is not necessary, since the resulted from (17) values are constantly optimized within the hybrid algorithm.Additionally, this stage cannot be viewed as initial ''guess'', since the information that is used does not deal with the ARMAX model in question.
It is rather an adaptation of the hybrid optimization algorithm to become problem specific.
In the second stage the ARMAX(na, nb, nc, nk) model is estimated by means of the hybrid algorithm.The parameter vector now becomes (18) www.intechopen.com Frontiers in Evolutionary Robotics 282 with c i denoting the additional parameters, due to the presence of C(q) polynomial.The values c i are randomly chosen from the normal distribution.The hybrid algorithm is presented below, counter j :=0; initialize (0) evaluate (0) while T{ (j)}≠ true do recombine: '(j)= (j)) evaluate: (j) mutate: (j)=m( (j)) evaluate: (j) select: (j+1)=s( (j)) j:=j+1 end while where j is the iteration counter, (j) the current population, T the termination criterion, r the recombination operator, m the mutation operator (provided by the BFGS) and s the selection operator (see Sec. 3).The evaluation of parameter vector at each iteration is realized via the calculation of objective function.For the successful realization of the hybrid algorithm, two issues must be further examined: the choice of the predictor, which modulates the residual sequence and the choice of the objective function, from which this sequence is evaluated at each iteration.An additional topic is the choice of the ''best'' model among a number of estimated ones.This topic is covered by statistical tests for order selection.

Choice of predictor
It is obvious that in every iteration of the hybrid algorithm the parameter vector (j) is evaluated, in order to examine its quality.Clearly, this vector formulates a corresponding ARMAX model with the ability to predict the output .For parametric models there is a large number of predictor algorithms, whose functionality depends mostly on the kind of the selected model, as well as the occasional scope of prediction.For the ARMAX case, a well-suited, one step-ahead predictor is stated in (Ljung, 1999) and has the following form: In Eq. ( 19) the predictor can be viewed as a sum of filters, acting upon the data set and producing model's output at time t+1.Both of these filters have the same denominator dynamics, determined by C(q), and they are required to be stable, in order to predict stable outputs.This is achieved by letting the roots of C(q) have magnitude greater than one, requirement which coincides with the invertibility property of H(q) transfer function (see Sec. 2).

Choice of objective function
The choice of an appropriate objective function performs vital role in any optimization problem.In most of the cases the selection is problem-oriented and, despite its importance, this topic is very often undiscussed.However, in optimization theory stands as the starting point for any numerical algorithm, deterministic or stochastic, and is in fact the tool for transmitting any given information of the test case into the algorithm, in a way that allows functionality.For the ARMAX parameter estimation problem, the objective function that has been designed lies in the field of quadratic criterion functions, but takes a slightly different form, which enforces the adaptivity of the hybrid optimization algorithm to the measured data.The objective function (of ) is computed from the following pair of equations ( 20) where ( 21) Equation ( 21) can be considered as the ratio of the absolute error integral to the absolute response integral.When fit reaches the value 100, the predicted time-series is identical with the measured one.In this case, of results to zero, which is the global minimum point.Nevertheless, it must be noted that for a specific parameter vector, the global minimum value of the corresponding objective function is not always equal to zero, since the selected ARMAX structure may be unable to describe the dynamics of the true system.The proposed method, as already mentioned, guarantees the stability of the estimated ARMAX model, by penalizing the objective function when at least one root of A(q) or C(q) polynomials lies within the unit circle of the complex plane.Thus, the resulted models satisfy the required conditions stated in Sec. 2.

Model order selection
The selection of a specific model among a number of estimated ones, is a matter of crucial importance.The model which shall be selected for the description of the true system's dynamics, must have as small over-determination as possible.There is a large number of statistical tests that determine model order selection, but the most common are the Residual Sum of Squares (RSS) and the Bayesian Information Criterion (BIC).The RSS criterion is computed by a normalized version of (9), that is, (22) where ||.|| denotes the Euclidian norm.The RSS criterion generally leads to overdetermination of model order, as it usually decreases for increasing orders.The BIC criterion overcomes this fact by penalizing models with relatively high model order (23) Clearly, both of the methods indicate the ''best'' model, as the one minimizing ( 22) and ( 23) respectively.

Implementation of the method
In this Section the proposed methodology is implemented to the identification of a testing apparatus described below, by means of SISO ARMAX models.The process of parametric modelling consists of three fundamental stages: In the first stage, the delay nk shall be determined, in the second stage ARMAX(na,nb,nc,nk) models will be estimated, using the method described in Sec. 4, while in the third, the corresponding (selected) ARMAX model will be further examined and validated.

The testing apparatus
The testing apparatus is presented in Fig. 2. A motion control card, through a motor drive unit, which controls a brushless servomotor, guides a flexible robotic arm.A piezoelectric sensor is mounted on the arm's free end and acquires its transversal acceleration, by means of a DAQ device.The transfer function, considered for estimation in this study, is that relating the velocity of the servomotor with the acceleration of the arm.The velocity signal selected to be a stationary, zero-mean white noise sequence.The sampling frequency was 100 Hz and the number of recorded data was N = 5000, for both input and output signals.

Post-treatment of data
The sampled acceleration signal found to be strongly non-stationary, with no fixed variance and mean.Thus, stationarity transformations made before the ARMAX estimation phase.Firstly, the BOX-COX (Box et al., 1994) transformation was used with BC = 1.1, for the stabilization of variance, and afterwards difference transformations with d= 2 for the stabilization of mean value were implemented.The resulted acceleration signal, as well as the velocity one, was zero-mean subtracted.The final input-output data set is presented in Fig. 3.For the estimation of acceleration's spectral density, Thompson's multi-paper method [23] has been implemented, with time-bandwidth product n T = 4, and number of Fast Fourier Transforms N FFT = 2 13 .The estimated spectral density is presented in Fig. 4. Clearly, there are three spectral peaks appearing in the graph, corresponding to natural frequencies of the system at about 2, 5 and 36 Hz, while an extra area of natural frequency can be considered at about 17 Hz.An additional inference that can be extracted, is the high-pass trend of the system.For the determination of delay, ARMAX(k, k, k, nk) models were estimated, with k = 6, 7, 8, 9 and nk = 0,1,2,3.The resulted models were evaluated using the RSS and BIC criterions.Figure 5 presents the values of the two criterions, with respect to model order, for the various delays.The models with delay nk = 3 presented better performance and showed smaller deviations between them.Thus the delay of the system was set to nk = 3.This is an expected value, due to the flexibility of the robotic arm, which corresponds to delayed acceleration responses in it's free end.

Estimation of ARMAX(na, nb, nc, 3) models
The selection of an appropriate ARMAX model, capable of describing input-output dynamics, and flexible enough to manage the presence of noise, realizes through a threephase procedure:

•
In the first phase, ARMAX(k, k, k, 3) models where estimated, for k = 6 : 14.The low bound for k is justified from the fact that the three peaks in the estimated spectrum (see Fig. 4), correspond to three pairs of complex, conjugate roots of A(q) characteristic polynomial.The upper bound was chosen in order to avoid over-determination.The resulted models qualified via BIC and RSS criteria and the selected model was ARMAX (7,7,7,3).

•
The second phase dealt with the determination of C(q) polynomial.Thus, ARMAX (7, 7, k, 3) models were estimated, for k = 2 : 16.Again, BIC and RSS criteria qualified ARMAX (7,7,7,3) as the best model, and also the one with the lowest variance of residuals' sequence.

•
In the third stage ARMAX(7, k, 7, 3) models were estimated for the selection of B(q) polynomial.Using the same criteria, the ARMAX(7, 6, 7, 3) model was finally selected for the description of the system.The implementation of the proposed methodology in the above procedure presented satisfying performance, as the hybrid optimization algorithm presented quick convergence rate despite model order, the resulted models were stable and invertible and overdetermination was avoided.

Validation of ARMAX(7,6,7,3)
For an additional examination of the ARMAX (7,6,7, 3) model, some common tests of it's properties have been implemented.Firstly, the sampled autocorrelation function of model residuals was computed for 1200 lags and it is presented in Fig. 6.It is clear that, except few lags, the residuals are uncorrelated (within the 95% confidence interval) and can be considered white.In Fig. 7, the frequency response of the transfer function G(q) (see ( 6)) is presented.The high-pass performance of the system is obvious and coincides with the same result that was extracted from the estimated spectral density in Fig. 4. Figure 8 displays a simulation of the system, using a fresh data set that was not used in the estimation tasks (the validation set).The dash line represents model's simulated acceleration, while the continuous line the measured one.
Figure 8.One second simulation of the system Finally, in Table 2 the natural frequencies in Hz and the corresponding percentage damping of the model are presented.While the three displayed frequencies were detected successfully, the selected model was unable to detect the low frequency of 2 Hz, probably due to it's high-pass performance.Yet, this specific frequency was unable to detected, even from higher order estimated models.

Conclusion
In this paper a new method for the estimation of SISO ARMAX models was presented.The proposed methodology lies in the context of Evolutionary system identification.It consists of a hybrid optimization algorithm, which interconnects the advantages of its deterministic and stochastic components, providing superior performance in PEM, as well as a two-stage estimation procedure, which yields only stable models.The method's main characteristics can be summarized as follows: • improvement of PEM is implemented through the use of a hybrid optimization algorithm, • initial ''guess'' is not necessary for good performance, • convergence in local minima is avoided, • computational complexity is sufficiently decreased, compared to similar methods for Evolutionary system identification.Furthermore, the method has competitive convergence rate to conventional gradient-based techniques, • stability is guaranteed in the resulted models.The unstable ones are penalized through the objective function, • it is successive, even in the presence of noise-corrupted measurement.The encouraging results suggest further research in the field of Evolutionary system identification.Specifically, efforts to design more flexible constraints are taking place, while the implementation of the method to Multiple Input-Multiple Output structures is also a topic of current research.Furthermore, the extraction of system's valid modal characteristics (natural frequencies, damping ratios), by means of the proposed methodology, is an additive problem of crucial importance.Evolutionary system identification is an growing scientific domain and presents an ongoing impact in the modelling of dynamic systems.Yet, many issues have to be taken under consideration, while the knowledge of classical system identification techniques and, additionally, signal processing and statistics methods, is necessary.Besides, system identification is a problem-specific modelling methodology, and any possible knowledge of the true system's performance is always useful.
Figure 1 presents the topology of the Fletcher & Powell test function for n = 2.The maximum number of objective function evaluations is 2•10 5 .In order to obtain statistically significant data, a sufficiently large number of independent tests must be performed.Thus, the results of N = 100 runs for each algorithm were collected.The expectation is estimated by the average:

Figure 2 .
Figure 2. The testing apparatus

Figure 3 .
Figure 3.The input-output data set

Figure 4 .
Figure 4. Acceleration's estimated spectrum For the subsequent tasks, and after excluding the first 500 points to avoid transient effects, the data set was divided into two subsets: the estimation set, used for the determination of

Figure 5 .
Figure 5. Selection of system's delay

Table 1 .
Results on the Fletcher & Powell function for n = 20

Table 2 .
Natural frequencies and corresponding damping