Helon Vicente

Hultmann Ayala

Didace Habineza
email: habinezadidace@gmail.com

Micky Rakotondrabe

Leandro Dos
email: leandro.coelho@pucpr.br

Santos Coelho

Nonlinear black-box system identification through coevolutionary algorithms and radial basis function artificial neural networks

Keywords:

The present work deals with the application of coevolutionary algorithms and artificial neural networks to perform input selection and related parameter estimation for nonlinear black-box models in system identification. In order to decouple the resolution of the input selection and parameter estimation, we propose a problem decomposition formulation and solve it by a coevolutionary algorithm strategy. The novel methodology is successfully applied to identify a magnetorheological damper, a continuous polymerization reactor and a piezoelectric robotic micromanipulator. The results show that the method provides valid models in terms of accuracy and statistical properties. The main advantage of the method is the joint input and parameter estimation, towards automating a tedious and error prone procedure with global optimization algorithms.

Introduction

The field of black-box nonlinear system identification [START_REF] Sjöberg | Nonlinear black-box modeling in system identification: a unified overview[END_REF] deals with the creation of data-driven mathematical abstractions of dynamic systems with little or no information about its intrinsic properties. The range of its applications spans various domains [START_REF] Valdivia | Blackbox modeling of three-phase voltage source inverters for system-level analysis[END_REF][START_REF] Valdivia | Black-box behavioral modeling and identification of dc-dc converters with input current control for fuel cell power conditioning[END_REF][START_REF] Kim | Nonlinear system identification of large-scale smart pavement systems[END_REF][START_REF] Da Costa Lopes | A control-oriented model of a PEM fuel cell stack based on NARX and NOE neural networks[END_REF][START_REF] Gil | Recursive subspace system identification for parametric fault detection in nonlinear systems[END_REF][START_REF] Tutunji | Parametric system identification using neural networks[END_REF][START_REF] Garrido | A heuristic approach to outputonly system identification under transient excitation[END_REF]. The model construction procedure in nonlinear black-box system identification is experimental and involves subjective decisions [START_REF] Nørgård | Neural Networks for Modelling and Control of Dynamic Systems: A Practitioner's Handbook[END_REF]. After the definition of the experiment protocol, the excitation signal should be devised so that the measured data contains information about the system dynamics. Once data is available, the model structure and related parameters are set so it is possible to evaluate the model quality and assess its adherence to validation metrics. If the model fails to meet the validation criteria, depending on the outputs in terms of model accuracy and residual properties the engineer should reiterate and construct a new model. Being so, the data-driven modeling activity involves subjective decisions and at times tedious and * Corresponding author.

error-prone activities. In the present contribution we try to address one of such issues, focusing on the automatic definition of the inputs of the model.

In this context, Computational Intelligence (CI) [START_REF] Engelbrecht | Computational Intelligence: An Introduction[END_REF] paradigms are a source of algorithms in order to alleviate from the engineer some procedures when creating the simulation artifacts. For a review on swarm and evolutionary computing applied to system identification, see [START_REF] Gotmare | Swarm and evolutionary computing algorithms for system identification and filter design: A comprehensive review[END_REF]. Parameter estimation of infinite impulse response filter models through system identification is made in [START_REF] Panda | IIR system identification using cat swarm optimization[END_REF] with a cat swarm algorithm. In [START_REF] Sa'ad | A robust structure identification method for evolving fuzzy system[END_REF] the authors introduce a structure identification procedure to evolve fuzzy models using incremental partitioning learning. The cuckoo search algorithm has been applied to optimize an adaptive Hammerstein model in [START_REF] Gotmare | Nonlinear system identification using a cuckoo search optimized adaptive hammerstein model[END_REF]. Solar radiation prediction, which is important in the context of renewable energy dispatch, is made in [START_REF] Wang | Solar radiation prediction using different techniques: model evaluation and comparison[END_REF] using different types of Artificial Neural Networks (NNs). In [START_REF] Alvernaz | Autoencoder-augmented neuroevolution for visual doom playing[END_REF] the Covariance Matrix Adaptation Evolution Strategy (CMA-ES) is implemented to perform neuroevolution of NN controllers for playing the Doom game in visual mode. Evolutionary Algorithms (EAs) are used to discover deep architectures for NNs in [START_REF] Real | Large-scale evolution of image classifiers[END_REF] and test it in image classification problems, where no human intervention is needed for describing the architecture. Pan evaporation prediction is carried out in [START_REF] Wang | Pan evaporation modeling using six different heuristic computing methods in different climates of china[END_REF], showing overall better results for CI methods such as NN and neuro-fuzzy, and in [START_REF] Wang | Evaporation modelling using different machine learning techniques[END_REF] with a genetic fuzzy hybrid method. In [START_REF] Barret Zoph | Neural architecture search with reinforcement learning[END_REF] the authors use reinforcement learning to generate recurrent NNs for the tasks of image classification and language modeling. Genetic programming has been applied to evolve deep NN structures also for the task of image classification in [START_REF] Suganuma | A genetic programming approach to designing convolutional neural network architectures[END_REF]. In [START_REF] Qiao | A self-organizing deep belief network for nonlinear system modeling[END_REF] a self-organizing strategy for constructing deep belief networks is proposed and tested for an artificial system, where the architecture is defined by the method iteratively. A variable structure fuzzy wavelet neural network controller has been proposed in [START_REF] Solgi | Variable structure fuzzy wavelet neural network controller for complex nonlinear systems[END_REF], which can also be applied to nonlinear adaptive modeling. In [START_REF] Qin | Improving the estimation of daily aerosol optical depth and aerosol radiative effect using an optimized artificial neural network[END_REF] the authors use a hybrid training procedure for NNs with Genetic Algorithm (GA) and backpropagation to predict aerosol optical depth.

Bringing CI methods to the context of system identification is however not straightforward. It is necessary to address how these methods adhere to the system identification procedure towards the automation of model generation, which would yield more efficient and cheaper model building activity. Indeed, by automating the tasks in model construction one needs to investigate application specific issues, such as how to avoid to obtain models overly complex which might not be applicable to design feedback control laws. In fact, there are many sources of inspiration which could be drawn in order to build more accurate models with system identification relying more in CI methods.

CoEAs [START_REF] Potter | Cooperative coevolution: An architecture for evolving coadapted subcomponents[END_REF] are based on the interaction of different populations and their successes in order to evolve the individuals for decision making [START_REF] Popovici | Coevolutionary principles[END_REF]. Individuals are compared not by their fitness as in traditional EAs, but by their outcomes in the interaction with other populations. Many problems have been tackled by such algorithms by decomposition as results appeared in the literature [START_REF] Wang | A differential evolution with cooperative coevolutionary selection operation for high-dimensional optimization[END_REF] and also for evolving NNs [START_REF] Gomez | Accelerated neural evolution through cooperatively coevolved synapses[END_REF]. CoEAs have been applied in the context of system identification or function approximation. Reference [START_REF] Delgado | Coevolutionary genetic fuzzy systems: a hierarchical collaborative approach[END_REF] implements a CoEA with GAs to evolve Takagi-Sugeno (TS) fuzzy models and further extended to soft sensor design in [START_REF] Delgado | A neuro-coevolutionary genetic fuzzy system to design soft sensors[END_REF] to deal with high-dimensional input space. In [START_REF] Parras-Gutierrez | Short, medium and long term forecasting of time series using the L-Co-R algorithm[END_REF][START_REF] Parras-Gutierrez | A Radial Basis Function Neural Network-Based Coevolutionary Algorithm for Short-Term to Long-Term Time Series Forecasting[END_REF] the authors implement two populations of a GAs to represent both the lag space and the related parameters of Radial Basis Functions Neural Networks (RBFNNs) for the purpose of time-series forecasting.

In the present paper, we propose the use of a CoEA to perform model-based input selection and related parameter estimation for the special case of system identification and RBFNN models. The coevolutionary approach was able to build accurate and valid models as shown later in the results. This approach is flexible and we may employ different or improved EAs in order to achieve better results. Different types of nonlinear mappings may be used as the metaheuristic paradigm for optimization is very flexible in facing different classes of problems and does not require the cost function to be neither differentiable or convex. The novelties of the present work are highlighted in the following. In the work of Wang and Li [START_REF] Wang | A coevolutionary differential evolution with harmony search for reliability-redundancy optimization[END_REF], the authors propose a CoEA framework for solving the reliability-redundancy problem by dividing the resolution of the problem with continuous and integer decision variables to different EAs. In this contribution, the authors use a Differential Evolution (DE) algorithm for the continuous and an improved Harmony Search (HS) algorithm for the integer variables. The problem of input and parameter estimation can be formulated similarly as in the work from Wang and Li, where in our case the continuous variables represent the NN parameters and the binary variables represent the presence/absence of a lagged input/output in the model. Herein we introduce an improved architecture for the model building problem in system identification based on the concept of CoEA. Specifically, we propose the use of the DE algorithm for optimizing the continuous variables and a recently introduced Adaptive Binary Harmony Search (ABHS) for choosing the inputs of the model. To the best of our knowledge the mentioned CoEA approach has not been applied so far for the purpose of model input selection and related parameter estimation in nonlinear black-box system identification problems. It is important to highlight that, by using a similar method as in [START_REF] Wang | A coevolutionary differential evolution with harmony search for reliability-redundancy optimization[END_REF] proposed to other domain, it is possible to solve concomitantly the problem of selecting the model inputs and estimating the related model parameters. This is important in the modeling activity as saves the engineer time spent on defining the model architecture and estimating the parameters in an iterative procedure. We test the novel methodology with three case studies: a MR damper, a continuous polymerization reactor, and a hysteretic piezoelectric actuator that the authors proposed in a previous work [START_REF] Ayala | Nonlinear black-box system identification through neural networks of a hysteretic piezoelectric robotic micromanipulator[END_REF]. The results are validated with the tests based on higher-order correlation functions of the residual and the input and the multiple correlation coefficient, showing competitive results when compared to results found in the literature. Moreover, from the results it can be seen that the hysteretic behavior of the piezoelectric micromanipulator is adequately captured, what motivates the application of the proposed method for building sensorless controllers which are important in this domain. It is also important to highlight that the framework proposed in this paper leverages the automation of the data-driven model building activity, what is desirable with the increasing complexity of systems and the pervasiveness of simulation in various industries for building more accurate models using less human resources. Fig. 1 summarizes the approach presented in this paper.

The remainder of the paper is organized as follows. We state the concepts used in the scope of nonlinear black-box system identification input selection and parameter estimation, together with the formulation of the problem in Section 2. The EAs stepby-step procedures are given in Section 3 while in Section 4 the proposed methodology is detailed. The case studies used as references for checking the accuracy of the method are given in Section 5. The results of the proposed methodology are detailed in Sections 6 and 7 concludes the document, giving general remarks about the study and also future research directions which will be pursued.

Nonlinear black-box system identification

This section introduces the mathematical formulation of the NN models used for system identification and their validation metrics. Moreover, we state the problem formulation in this context for input selection and parameter estimation which will be the focus of the CoEA methodology proposed for system identification.

In the present work we restrict our attention to Nonlinear AutoRegressive with eXogenous inputs (NARX) models which are defined as

y(t) = F [y(t -1), y(t -2), . . . , y(t -n y ), u(t -1), u(t -2), . . . , u(t -n u )] + ξ (t). (1)
where u(t) and y(t) are the input and output of the system at a given instant t, considered to be measured. The residual is given by ξ (t) = y(t) -ŷ(t), where ŷ(t) is the predicted output. Let F [•] be in general a nonlinear function mapping in R n φ → R, from the n φ model inputs to the predicted output. In this paper we define

F [•] as a RBFNN, thus we have ŷ(t) = F [r r r(t)] = M m=1 w m φ(r r r(t), c c c m , σ m ), (2)
where M ∈ N + is the number of neurons in the hidden layer and r r r(t) ∈ R nr is the input vector at a given instant t; c c c m ∈ R nr , σ m ∈ R + , and w m ∈ R are the RBFNN parameters set during training, respectively the center and the width of the mth hidden node and the output weights. We use the Gaussian activation function for φ(•) so that we have

φ(r r r(t), c c c m , σ m ) = exp - 1 2σ 2 m nr i=1 (r i (t) -c m,i ) 2 .
(

Now we are ready to state the following two problems which will be tackled in the present work.

Problem 1 (Input Selection). For systems like the ones described in Eq. ( 1), select the model inputs of lagged terms on u(t) and y(t) such that the model represents the physical system accordingly.

Problem 2 (Parameter Estimation). Estimate the related parameters of F [•] for the inputs chosen for the model Eq. ( 1), in a way that the model is valid according to some predefined measure.

Both problems are critical whenever building a model. It is well known that inadequate inputs lead unavoidably to poor model abstractions. Moreover, often this part is left dependent on the engineer which chooses in a tedious trial and error procedure a suitable set of inputs. Thus, automatic methods to select appropriate inputs which are able to give valid models are highly desirable. Also, the parameter estimation naturally impacts the overall model quality. The engineer builds a model for each new set of inputs of the models chosen, as this part is dependent on input selection. In the following section, we introduce a novel methodology for input selection and parameter estimation based on the concept of coevolution and problem decomposition, which aims at the resolution of Problems 1 and 2 concomitantly in the context of system identification.

Let us introduce some model validation metrics in the following. The multiple correlation coefficient (R 2 ) is calculated as [START_REF] Haber | Structure identification of nonlinear dynamic systems-a survey on input/output approaches[END_REF]

R 2 = 1 - N t=1 ξ (t) 2 N t=1 (y(t) -ȳ) 2 (4)
where the upper bar denotes the mean value of the sequence. The R 2 measures the adherence of the model to the measured data, such that R 2 = 1 means perfect data reconstruction and R 2 > 0.9 may be considered sufficient for many applications [START_REF] Schaible | Fuzzy logic models for ranking process effects[END_REF].

It can be calculated for predictions in (i) One-Step-Ahead (OSA) and (ii) Free-Run simulation (FR). While (i) uses the most recent measured data to perform predictions, in (ii) predictions are made recursively. In general FR-based metrics are more reliable when assessing the model predictive capability and show more easily if the model is valid for longer prediction horizons.

The statistical correlation can also be used for model validation. It indicates whether the estimation procedure was able to explain the data provided for the estimation phase. Usually OSA predictions are taken for calculating the correlations, as they are more often used as metrics for estimating parameters due to computational issues. According to the set of tests given in [START_REF] Billings | Properties of neural networks with applications to modelling non-linear dynamical systems[END_REF] for NNs, we have

             φ ξ ξ (τ ) = δ(τ ), φ uξ (τ ) = 0, ∀τ , φ ξ (ξ u) (τ ) = 0, τ ≥ 0, φ (u 2 ) ′ ξ (τ ) = 0, ∀τ , φ (u 2 ) ′ ξ 2 (τ ) = 0, ∀τ , (5)
where δ(•) is the Kronecker delta function, (u 2 )

′ (t) = (u(t)) 2 -ū2 , (ξ u) = ξ (t +1)u(t +1
) and φ ab is the normalized cross-correlation function between two sequences {a} and {b}, which is given by

φ ab (τ ) = N-τ t=1 [a(t) -ā] b(t + τ ) -b N t=1 [a(t) -ā] 2 N t=1 [b(t) -b] 2 1/2 . (6)

Evolutionary algorithms

This section briefly outlines the DE and HS algorithms, which compose the CoEA system identification approach. As will be shown, they are used as global search methods in order to jointly define the inputs and the parameters of the black-box models represented by Eq. ( 1). In the following we depict the step-bystep procedure of both algorithms. All initial solutions are initialized randomly within the search space using an uniform distribution, and the optimization problems are all set as minimization with no loss of generality.

Differential evolution

The notation used for the description of the DE algorithm is given in Table 1. The steps necessary for implementing computationally the algorithm are as follows.

• Step 1: Perform mutation;

v v v i (g) = x x x r 1 (g) + F x x x r 2 (g) -x x x r 3 (g) (7)
• Step 2: Perform binary crossover;

u i,j (g) = v i,j (g), if rand < CR x i,j (g), otherwise, ( 8
)
where rand is a random number distributed uniformly in the range [0, 1]. • Step 3: Evaluation of the newly generated solutions and selection;

x x x i (g + 1) = u u u i (g), if f (u u u i (g)) < f (x x x i (g)) x x x i (g), otherwise, (9)
• Step 4: Increment generation counter g = g + 1 and check termination criterion. If g < G, go to Step 1, otherwise terminate.

The DE algorithm just exposed is used in the present paper for evolving the NN parameters, as will be discussed in the next section.

Adaptive binary harmony search algorithm

The ABHS algorithm used in the scope of the present work was proposed in [START_REF] Wang | An improved adaptive binary harmony search algorithm[END_REF]. In this paper, the authors study many adaptive mechanisms to adjust the project parameters of the HS algorithm with binary encoding. In the present work it will be used to denote the evolve the set of candidate lags. The description of the algorithm is given below, while in Table 2 we state the variables used.

• Step 1: Adaptively update PAR and HMCR:

PAR = PAR 0 + PAR f -PAR 0 G • g, (10a) HMCR = HMCR 0 + HMCR f -HMCR 0 G • g. (10b)
•

Step 2: New harmony improvisation, considering Boolean decision variables -Harmony memory consideration;

1. Bit selection strategy; 2. Individual selection strategy;

h new j = h pj , if r 1 < HMCR round(r 2 ), otherwise (11)
CR Crossover probability, in [0, 1] x x x i (g) i-th individual vector v v v i (g) i-th mutant vector u u u i (g) i-th trial vector r 1 , r 2 ,
h new j = h tj , if r 1 < HMCR round(r 2 ), otherwise (12)
where t = round (r 3 • HMS).

-Pitch adjustment;

h new j = h bj , if r 1 < PAR h new j , otherwise (13)
where b is the index of the best harmony found so far by the algorithm.

•

Step 3: Update HM by replacing the worst stored harmony by h h h new , if it is better than any of the ones stored.

• Step 4: Update generation counter g = g + 1. Check termination criterion, go to Step 1 if g < G or terminate otherwise.

The ABHS algorithm is used to perform the input selection in the CoEA scheme given in the next section.

Proposed methodology

Now we shall focus to state the proposed methodology based on CoEA to perform input selection and related parameter estimation. We adopt the CoEA paradigm with two populations, namely H and X , in order to decompose the problem of system identification into lag selection and parameter estimation respectively. Being so, H is composed of a set of candidate lags and is evolved by a binary optimization algorithm while X represents the corresponding RBFNN model parameters (centers and widths) evolved by the continuous optimization algorithm. We shall denote hereafter both populations at the gth iteration respectively by H(g) and X (g).

The individuals in H and X coevolve during the procedure, and their individuals' fitnesses are calculated according to the ability of the individuals of both populations. It means that if one population performs inappropriately, then unavoidably the performance of the other population will also decrease. In this way, the individuals cooperate as both fitnesses depend on the composed solution formed by them. Each individual of the populations h h h ∈ {0, 1} D H in H and x x x ∈ R D X in X is evaluated according to Eq. ( 14), which represents the objective function to be minimized based on the Mean Squared Error (MSE)

f (w w w) = MSE = 1 N N t=1 [ξ (t)] 2 (14)
where N represents the amount of data used and w w w denotes the solution formed by the concatenation of both decision vectors

w w w h h h ⊤ x x x ⊤ ⊤ . (15)
We show in Fig. 2 how the individuals of each population are composed. From this illustration we can see that in H the absence or presence of a given lag is indicated by binary variables and in X the individuals encode the corresponding model parameters for the RBFNN, composed by the centers and widths (the output weights are defined by QR factorization [START_REF] Golub | Matrix Computations[END_REF]).

The dimensions of the decision variables x x x and h h h depend on the project parameters given by (i) n max y (maximum number of lags in y(t)); (ii) n max u (maximum number of lags in u(t)); and (iii) M (number of neurons of the RBFNN). Being so, the total dimensions of the individuals in H and X are given as defined in Eq. ( 16):

D H = n max y + n max u , (16a) D X = (n max y + n max u + 1) • M. (16b)
The objective function aims at the definition of the model parameters based on the MSE of the OSA prediction, which is a quite standard metric for building models. Here, we could in principle employ different types of models but we rather focus on the case of RBFNNs. They are proven to be global approximators [START_REF] Park | Universal approximation using radial-basis-function networks[END_REF] being very convenient from the architecture point of view, should the designer select solely the complexity of the RBFNN. In this approach we adopt the framework for coevolution as described in [START_REF] Wang | A coevolutionary differential evolution with harmony search for reliability-redundancy optimization[END_REF] to solve the optimization problem with the cost function in Eq. [START_REF] Gotmare | Nonlinear system identification using a cuckoo search optimized adaptive hammerstein model[END_REF]. In that work, the authors tackle the problem of reliability-redundancy optimization, which is formed by both integer and continuous variables. They adopted a population to handle the integer variables and another to cope with the continuous ones, which coevolve concomitantly in a cooperative fashion. They proposed a modified HS algorithm to deal with the integer variables and a standard DE to manipulate the continuous variables. After half of the overall evolution process, the integer optimization stops and the best candidate is held for the rest of the half of the algorithm, where the continuous variables are further refined.

The framework in [START_REF] Wang | A coevolutionary differential evolution with harmony search for reliability-redundancy optimization[END_REF] is not directly applicable to our problem as H is formed with binary individuals, not integers. Thus, we substituted the modified HS by an ABHS algorithm recently proposed [START_REF] Wang | An improved adaptive binary harmony search algorithm[END_REF], with individual selection strategy. In principle we could employ any other metaheuristic optimization algorithm capable of dealing with binary decision variables. In this case there are many possibilities. ABHS was one of the possible choices, but the best candidate so as to change as little as possible the original framework given by [START_REF] Wang | A coevolutionary differential evolution with harmony search for reliability-redundancy optimization[END_REF] and not increase the number of project parameters. Figs. 3 and4 show the overall procedure and the representation of the individuals for the coevolutionary system identification approach.

It is important to highlight that we face the problem of input selection and related parameter estimation for system identification as different subproblems which are solved separately in a cooperative fashion. We adopt one population to represent the lags and the other that encodes the parameters of a RBFNN. Being so, we are able to cope with the problem of selecting the input variables and the parameters of the model concomitantly. In this way the designer is not obliged to set the orders of the Fig. 4. Coevolutionary approach adopted for system identification. In this figure we summarize the individual representation for both populations and how they relate to the RBFNN model for black-box system identification.

model beforehand, what is very convenient specially when e.g. (i) little knowledge about the system is available, (ii) little time is available to spend on trial and error methods for input selection, (iii) there are many inputs and outputs on the system studied, what makes trial and error methods for input selection even more time-consuming.

Case studies

This section is devoted to describe the benchmarks used in the scope of the present work. We describe three dynamic systems given below which will be modeled using their input-output data through the CoEA approach for system identification.

(i) MR Damper: it is a device used for structural vibration control, as it may change its viscosity actively. The data has been used for identification in e.g. [START_REF] Wang | Identification of Hammerstein systems without explicit parameterisation of non-linearity[END_REF][START_REF] Ayala | Cascaded evolutionary algorithm for nonlinear system identification based on correlation functions and radial basis functions neural networks[END_REF]. The input is the velocity at the ends of the damper, while the output is the force which it reacts to excitation, that are at 200 Hz.

(ii) Continuous polymerization reactor: it describes freeradical polymerization of MMA (methyl methacrylate) with AIBN (azobisisobutyronitrile) as initiator and toluene as solvent and is supposed that the reaction is made in a jacketed continuous stirred tank reactor. The nonlinear state space model based in first-principles is given in [START_REF] Doyle | Nonlinear model-based control using second-order Volterra models[END_REF]:

ż1 (t) = 10[6 -z 1 (t)] -2.4568 z 2 (t), (17a) ż2 (t) = 80u(t) -10.1022z 2 (t), (17b) ż3 (t) = 0.0024121z 1 (t) z 2 (t) + 0.112191z 2 (t) -10z 3 (t), (17c) ż4 (t) = 245.978z 1 (t) z 2 (t) -10z 4 (t), (17d) y(t) = z 4 (t) z 3 (t) , ( 17e
)
where u(t) is the dimensionless volumetric flow rate of the initiator; z 1 (t) is a dimensionless variable representing monomer concentration; z 2 (t) represents the initiator concentration; and y(t) is the output of the system, representing the number-average molecular weight. We generated 50,000 data at 5 Hz, with u(t) ∼ U[0.007, 0.015] [START_REF] Rhodes | Determining the model order of nonlinear input/output systems[END_REF], by solving the nonlinear state equations with the Runge-Kutta method. With the aforementioned input, the range of the output is roughly between 2.6 × 10 4 and 3.4 × 10 4 [START_REF] Rhodes | Determining the model order of nonlinear input/output systems[END_REF]. As in [START_REF] Madár | Genetic programming for the identification of nonlinear input-output models[END_REF], we use 960 data for the purpose of model building and validation. This system has been used for model input selection in [START_REF] Rhodes | Determining the model order of nonlinear input/output systems[END_REF][START_REF] Madár | Genetic programming for the identification of nonlinear input-output models[END_REF]. Reference [START_REF] Rhodes | Determining the model order of nonlinear input/output systems[END_REF] used this simulated system to exemplify an input selection algorithm for nonlinear models based on the false nearest neighbors algorithm using only the input and output data. In [START_REF] Madár | Genetic programming for the identification of nonlinear input-output models[END_REF], the authors apply a novel algorithm to build linear-in-the-parameters models for system identification using genetic programming and orthogonal least squares.

(iii) Piezoelectric micromanipulator: the microactuator, called unimorph, is a free-clamped cantilever composed of two layers: the piezoelectric layer based on lead-zirconate-titanate material (PZT), and a nonpiezoelectric layer based on copper material, see Fig. 5a. The nonpiezoelectric layer is called passive layer. When a voltage u is applied to the piezoelectric layer, it expands or contracts. Due to the constraint between the two layers, this expansion/contraction yields a global deflection (displacement) y of the cantilever (Fig. 5b). Such principle is widely exploited to perform precise positioning [START_REF] Agnus | Robotic microassembly and micromanipulation at femto-st[END_REF][START_REF] Binnig | Single-tube three-dimensional scanner for scanning tunneling microscopy[END_REF][START_REF] Lescano | Kinematic analysis of a meso-scale parallel robot for laser phonomicrosurgery[END_REF][START_REF] Mitsuhashi | Relay using multilayer piezoelectric actuator, Japan[END_REF][START_REF] Devasia | A survey of control issues in nanopositioning[END_REF] because of the nanometric resolution and the high dynamics of the deformation. However, the relation between the voltage u and the output displacement y is known to be hysteretic which drastically affects the final precision. Moreover, this hysteresis was demonstrated to be dynamic, such that the piezoelectric actuator finally exhibits nonlinearities [START_REF] Rakotondrabe | Multivariable classical Prandtl-Ishlinskii hysteresis modeling and compensation and sensorless control of a nonlinear 2-dof piezoactuator[END_REF]. The experimental setup, presented in Fig. 5c is composed of (i) the piezoelectric actuator with sizes of 15 mm × 2 mm × 0.3 mm (length × width × total thickness). The thickness of the piezoelectric layer is 0.2 µm while that of the passive layer is 0.1 µm; (ii) an optical sensor (LK2420 from Keyence company) which is used to measure the deflection (displacement) of the actuator. Its resolution is tuned to be 40 nm and its bandwidth to 5 kHz, which are sufficient for the tests

Results

In the present section, we describe the results obtained by the application of the proposed methodology based on CoEAs for system identification. Three examples are used to test it, namely (i) a MR damper, (ii) a continuous polymerization reactor and (iii) a piezoelectric micromanipulator which were detailed in Section 5.

The common project parameters used for the simulations in the case studies mentioned above are stated in Table 3, as indicated by [START_REF] Geem | Optimal cost design of water distribution networks using harmony search[END_REF][START_REF] Simon | Evolutionary Optimization Algorithms: Biologically-Inspired and Population-Based Approaches to Computer Intelligence[END_REF]. Specifically, for each case study (i), (ii) and (iii) we used G = 200, G = 400 and G = 100 respectively and the quantity of neurons on each RBFNN as M = 10 for cases (i) and (ii) and M = 13 for (iii), according to the complexity of the system to be modeled. It is important to highlight that the number of generations and model complexity needed for a given system should be set in order to meet model validation requirements. In the present work we established the number of generations so as not to waste computational resources with stalled evaluations after some experimentation. For all case studies we normalized the data in the range [-1,1]. The search range was set for the c m,i in [-4,4] and σ m in [0.01,20]. The initial population was generated randomly according to a uniform distribution. Now we focus on the statistical analysis of the final outcome of the proposed methodology. As the search procedures are stochastic, it is necessary to evaluate the outcome of the procedure with different initial conditions. Thus, each case study reported in Section 5 was identified with 30 different initial conditions of the CoEA algorithm. Table 4 shows the statistics of the values obtained by all runs of the CoEA optimization. It is possible to see that the standard deviation is relatively low when considering many different initial conditions. This indicates that the proposed algorithm was able to converge to similar solutions in the objective space irrespective to the initial conditions, what shows the robustness of the stochastic search approach. Fig. 6 shows the progression of the objective function value at each generation considering all runs. The mean value at each generation is plotted, together with a shaded area representing the bounds set with the standard deviation. From this graph it is possible to confirm what was mentioned before for Table 4 that the solutions vary slightly, even though working on completely different initial conditions. We further notice that the algorithm was able to converge after the defined number of generations as expected.

Let us now focus on the analysis of the multiple correlation coefficients obtained by the models, as given by Eq. (4). To this end, we denote by these coefficients obtained in estimation and validation phases respectively by R 2 e and R 2 v . In Table 5 we depict the metrics in terms of R 2 for the best models found, after running the procedure with 30 different initial conditions for the algorithm. By best model found we mean that we selected the run with minimal value found for R 2 v in FR. We can see that in the case of the polymerization reactor and the piezoelectric manipulator, the values of the R 2 metric are close to unity for estimation and validation phases and in OSA and FR, what represents excellent modeling capability. On the other hand, the MR damper case study presented R 2 close to unity in OSA and R 2 > 0.9 in FR what should be considered satisfactory as defined by [START_REF] Schaible | Fuzzy logic models for ranking process effects[END_REF] and by our practical appreciation. The excellent accuracy reported in Table 5 is also confirmed in the plots of the measured versus predicted variables. The predictions in OSA and FR can been seen in Figs. 7, 8 and 9 for all case studies, which illustrate R 2 coefficients reported in Table 5 as previously discussed.

Hereafter the results of the solutions presented in Table 5 are detailed. Table 6 shows the set of lags selected by the proposed methodology. In the case of the polymerization reactor, the results reported are in line with previous works which also used the same case study. In [START_REF] Rhodes | Determining the model order of nonlinear input/output systems[END_REF] the authors apply the false nearest neighbors algorithm to define the order of nonlinear input/output systems, as part of the algorithm herein presented performs. In this work the authors suggest the use of y(t -1), u(t -1), u(t -2) for identification. In the results reported in [START_REF] Madár | Genetic programming for the identification of nonlinear input-output models[END_REF], the authors use genetic programming together with orthogonal least squares to discover higher-order polynomial model structures for nonlinear system identification. They have found different sets of lags for various structures, allowing n y , n u up to 4 in the overall search procedure. For the most accurate results in average reported, 60% of the trials returned orders as in [START_REF] Rhodes | Determining the model order of nonlinear input/output systems[END_REF]. Thus, we emphasize that the present methodology, which performs joint input selection and related parameter estimation, presents similar results which are compatible with the ones previously reported in [START_REF] Rhodes | Determining the model order of nonlinear input/output systems[END_REF][START_REF] Madár | Genetic programming for the identification of nonlinear input-output models[END_REF]. Specifically, the coevolutionary based approach selected the same lags for u(t) and one additional lag on y(t), out of 10 possibilities which were allowed in each u(t), y(t) if compared to [START_REF] Rhodes | Determining the model order of nonlinear input/output systems[END_REF] per the CoEA algorithm parameterization (see Table 3). With respect to the MR damper case study, the lags found in the present work are simpler, when compared to the ones found by another method for joint input and parameter estimation method previously reported by the authors in [START_REF] Ayala | Cascaded evolutionary algorithm for nonlinear system identification based on correlation functions and radial basis functions neural networks[END_REF], as here we have one less lagged variable in each input and output. It is important to highlight that the set of lags are not unique when considering real-world measured data and the computational burden is much higher in [START_REF] Ayala | Cascaded evolutionary algorithm for nonlinear system identification based on correlation functions and radial basis functions neural networks[END_REF] to obtain the models, as for this method to generate a single model we need a couple of hours and the current method instead needs a couple of seconds. This is due to the fact that the problem is decomposed by the CoEA formulation, as opposed to the approach in [START_REF] Ayala | Cascaded evolutionary algorithm for nonlinear system identification based on correlation functions and radial basis functions neural networks[END_REF].

Now we analyze, for the micromanipulator case study, the ability of the generated model to capture its hysteretic behavior. It is important to highlight that this is another validation requirement for this particular system as this characteristic is known beforehand and important for its real-world application. The path following control algorithms are typically constructed

Table 6

Best regressors found by the coevolutionary based methodology.

Case study

Regressors chosen MR damper

y(t -1) y(t -2) u(t -1) u(t -2) u(t -6) u(t -7)
Continuous polymerization reactor

y(t -1) y(t -2) u(t -1) u(t -2)
Piezoelectric actuator

y(t -1) y(t -2) y(t -3) y(t -4) y(t -8) u(t -1) u(t -2) u(t -6) u(t -8) u(t -10)
on the basis of hysteresis models such as Bouc-Wen, Preisach, Phase-Preisach combined, Prandtl-Ishlinskii, among other approaches [START_REF] Rakotondrabe | Multivariable classical Prandtl-Ishlinskii hysteresis modeling and compensation and sensorless control of a nonlinear 2-dof piezoactuator[END_REF]. The goal is to use the information contained in the model about the hysteresis of the system in order to adequately compensate it for improving the accuracy of the regulation. In Fig. 10 we can see the plot of the output versus the input for the piezoelectric manipulator case study. This is very similar to the one presented in [START_REF] Ayala | Nonlinear black-box system identification through neural networks of a hysteretic piezoelectric robotic micromanipulator[END_REF]. In this paper the authors have solved this same identification problem testing manually many architectures; in contrast, in the present work we have selected the inputs and parameters of the model concomitantly, what represents a substantial gain in computational effort. It can thus be seen that the hysteretic behavior has been adequately captured. This reinforces the great opportunity in the application of data-driven modeling tools in order to design the mathematical abstraction for systems with higher complexity as demonstrated in a recent work which investigated 2-DOF piezoelectric micromanipulators [START_REF] Ayala | Modeling of a 2-dof piezoelectric micromanipulator at high frequency rates through nonlinear black-box system identification[END_REF] with ad-hoc model structures for NN. Thus, machine learning-based feedback control such as [START_REF] De | Discrete time control based in neural networks for pendulums[END_REF] can be directly applied on the basis of the accurate models constructed with less intervention from the engineer, as we have demonstrated. Finally, Fig. 11 shows the plot of the statistical tests based on the autocorrelation of the residuals and the higher order cross-correlation between the residuals and the system's inputs, given in Eq. [START_REF] Da Costa Lopes | A control-oriented model of a PEM fuel cell stack based on NARX and NOE neural networks[END_REF]. In this case it can be seen that the models obtained are also valid in the statistical sense, since the residual analysis points that the dynamics present in the data has been adequately captured. This is important to demonstrate that the models that were built have acquired the dynamics present in the OSA residual information, which was the metric established for the CoEA approach, what confirms the results reported in terms of prediction quality.

Conclusion

In the present paper, we showed the application of a novel coevolutionary based methodology for black-box system identification with RBFNN models. We have shown that the present methodology has proven to be suitable for choosing the model orders and related parameters -saving thus time spent on choosing manually the orders of the system or resorting to decoupled methods for input selection and parameter estimation. We chose the ABHS and DE algorithms tuned with standard parameters to compose the overall model building scheme. The methodology was tested with three nonlinear case studies for system identification, giving accurate and valid results as the R 2 for OSA and FR showed as well as the statistical tests based on correlation functions of the residuals. A piezoelectric manipulator has been used in the scope of the present work, where accurate models have been built with real acquired data by using the present methodology, and the hysteretic behavior of the system is adequately captured.

In future works, other metaheuristic algorithms may be used in order to evaluate the performance of each combination for the task of coevolutionary system identification. To this end, we may make use of statistical tests in order to evaluate each combination of algorithms that compose the coevolutionary approach [START_REF] Derrac | A practical tutorial on the use of nonparametric statistical tests as a methodology for comparing evolutionary and swarm intelligence algorithms[END_REF] or test different approaches [START_REF] Puchalsky | Agribusiness time series forecasting using wavelet neural networks and metaheuristic optimization: An analysis of the soybean sack price and perishable products demand[END_REF]. Being so, we may establish what is the best symbiotic relation among algorithms in the CoEA approach. This will be important to devise new research directions on the improvement of each of the EAs that compose the CoEA framework. There is relevant work in the field of evolutionary computation in the creation of novel, more efficient and complex, metaheuristic optimization algorithms based on krill herd [START_REF] Gandomi | Krill herd: A new bio-inspired optimization algorithm[END_REF], beetles [START_REF] Kallioras | Pity beetle algorithm -a new metaheuristic inspired by the behavior of bark beetles[END_REF] and thermal exchange [START_REF] Kaveh | A novel meta-heuristic optimization algorithm: Thermal exchange optimization[END_REF] or yet DE improvements [START_REF] Tsili | A novel multiobjective lognormal-beta differential evolution approach for the transformer design optimization[END_REF] to mention a few. It is important that the field of data-driven modeling takes advantage of these new methodologies. Another further study that will deserve attention is the automatic definition of the structure of the nonlinear mapping F [•] and also newer approaches such as [START_REF] Zhang | Identification of continuous-time nonlinear systems: The nonlinear difference equation with moving average noise (NDEMA) framework[END_REF]. In the case of RBFNNs this regards to the automatic definition of the number of neurons. Some studies have been devoted to the generation of metrics to measure complexity in the scope of NNs applied to system identification [START_REF] Ugalde | Neural network design and model reduction approach for black box nonlinear system identification with reduced number of parameters[END_REF][START_REF] Ugalde | Computational cost improvement of neural network models in black box nonlinear system identification[END_REF], which may be applied in studies directed to the extension of the present methodology. A possibility in the specific case of RBFNNs is to adopt the two-steps methodology for training with a clustering technique that defines automatically the number of classes, as the X-means [START_REF] Pelleg | X-means: Extending K-means with efficient estimation of the number of clusters[END_REF] or hybrid approaches [START_REF] Subudhi | Nonlinear system identification using memetic differential evolution trained neural networks[END_REF]. Another idea which should be explored in the context of system identification is creation of ensembles, which is the use of many towards having more precise predictions [START_REF] Ribeiro | Enhanced ensemble structures using wavelet neural networks applied to short-term load forecasting[END_REF].

Declaration of competing

No author associated with this paper has disclosed any potential pertinent conflicts which may be perceived to have impending with this work. For full disclosure statements to https://doi.org/10.1016/j.asoc.2019.105990.

Fig. 1 .

 1 Fig. 1. Summary of the proposed methodology. The CoEA estimates the orders of the model and the model parameters in a cooperative fashion. The response of the model is compared to the one measured, and this information is used to evaluate the individuals in each population.

r 3 Uniformly

 3 random integers, mutually different with the target vector index, in the range[START_REF] Sjöberg | Nonlinear black-box modeling in system identification: a unified overview[END_REF] NP]

Fig. 2 .

 2 Fig. 2. Encoding of the individuals (a) h h h and (b)x x x for the proposed coevolutionary approach.

Fig. 3 .

 3 Fig. 3. Pseudocode for the coevolutionary system identification approach based on ABHS and DE.

Fig. 5 .

 5 Fig. 5. The experimental setup. (a) and (b): principle of functioning of the piezoelectric actuator. (c): the experimental setup diagram.

Table 3

 3 Parameters used in the CoEA for system identification.Parameter description ValueNumber of candidate lags on y(t) and u(t

Fig. 6 .

 6 Fig. 6. Evolution at each generation of the mean and standard deviation of the objective function for the (i) MR damper , (ii) polymerization reactor and (iii) piezoelectric actuator. In dotted lines we see the mean of the objective function value measured at each iteration of the algorithm, while in red we see area covered by the standard deviation around the mean value.

Fig. 7 .

 7 Fig. 7. Predictions and errors for the MR damper case study.

Fig. 8 .

 8 Fig. 8. Predictions and errors for the polymerization reactor case study.

Fig. 9 .

 9 Fig. 9. Predictions and errors for the piezoelectric actuator case study.

Fig. 10 .

 10 Fig.10. Plot of the u(t) versus the real output and the free simulation result, for the piezoelectric actuator case study. Note that the hysteretic behavior has been adequately captured.

Fig. 11 .

 11 Fig. 11. Tests based on the correlation of the residuals and inputs for the (i) MR damper case study, (ii) the continuous polymerization reactor and (iii) the piezoelectric actuator.

Table 1

 1 Description of the nomenclature for the variables used in the DE optimization algorithm.

	Symbol	Description
	NP	Number of individuals in a population
	D	Problem dimension
	g	Generation counter
	G	Total number of generations
	F	Scale factor, in [0, 1]

Table 2

 2 Description of the nomenclature for the variables used in the HS optimization algorithm.

	Symbol	Description
	HMCR 0 , HMCR f Initial and final harmony memory considering rate,
		respectively
	PAR 0 , PAR f	Initial and final pitch adjusting rate, respectively
	bw j	Pitch bandwidth for the jth dimension
	HM	Harmony matrix
	HMS	Number of solutions HM
	h h h i	i-th binary harmony stored in HM
	h h h

new

New solution generated according to ABHS strategies

r 1 , r 2 , r 3

Random numbers uniformly distributed in the range [0, 1]

Table 4

 4 Minimum (min.), maximum (max.), mean and standard deviation (std. dev.) of the results in terms of f (w w w) for the case studies analyzed.

	Case study	Min.	Max.	Mean	Std. dev.
	MR damper	10.5049	11.8069	11.2856	0.3273
	Pol. reactor (×10 4 )	0.9806	4.7628	2.4681	0.8897
	Piezoelectric actuator	0.012570	0.01273	0.01266	4.1034E-05

Table 5

 5 Multiple correlation coefficients for the case studies analyzed in OSA and FR, for estimation (est.) and validation (val.) phases.

	Case study	OSA		FR
		R 2	R 2	R 2	R 2
		(est.)	(val.)	(est.)	(val.)
	MR damper	0.9968	0.9971	0.9716	0.9591
	Polymerization reactor	0.9989	0.9993	0.9987	0.9992
	Piezoelectric actuator	0.99998	0.99998	0.99995	0.99995

Acknowledgments

authors would like to thank the financial support from CAPES (Coordination for the Improvement of Higher Education Personnel -Brazilian Government) grant no. 88887 . 163347/2018-00, the National Council of Scientific and Technologic Development of Brazil -CNPq (Grants number: 303906/2015-4-PQ, 404659/2016-0-Univ, 430395/2018-3-Univ, and 400119/2019-6-ERC) and FA (Fundação Araucária) PRONEX:042/2018 for their financial support.


Henok Moges

Thanos Manos
email: thanos.manos@u-cergy.fr

Charalampos Skokos

On the behavior of the Generalized Alignment Index (GALI) method for regular motion in multidimensional Hamiltonian systems

We investigate the behavior of the Generalized Alignment Index of order k (GALI k ) for regular orbits of multidimensional Hamiltonian systems. The GALI k is an efficient chaos indicator, which asymptotically attains positive values for regular motion when 2 ≤ k ≤ N , with N being the dimension of the torus on which the motion occurs. By considering several regular orbits in the neighborhood of two typical simple, stable periodic orbits of the Fermi-Pasta-Ulam-Tsingou (FPUT) β model for various values of the system's degrees of freedom, we show that the asymptotic GALI k values decrease when the index's order k increases and when the orbit's energy approaches the periodic orbit's destabilization energy where the stability island vanishes, while they increase when the considered regular orbit moves further away from the periodic one for a fixed energy. In addition, performing extensive numerical simulations we show that the index's behavior does not depend on the choice of the initial deviation vectors needed for its evaluation.

I. INTRODUCTION

Dynamical systems constitute a rather substantial research area with many different applications in several fields of physics. In particular, Hamiltonian systems of many degrees of freedom (dof) have been broadly used for the study and better understanding of energy transport and equipartition phenomena (see for example [START_REF] Mackay | Hamiltonian Dynamical Systems[END_REF][START_REF] Wiggins | Global Bifurcations and Chaos: Analytical Methods[END_REF][START_REF] Lichtenberg | Regular and Chaotic Dynamics[END_REF][START_REF] Simó | Hamiltonian Systems with Three or More Degrees of Freedom[END_REF][START_REF] Bountis | Complex Hamiltonian Dynamics[END_REF]). The latter phenomena are associated with the dynamical nature of the motion fostered in the respective phase space, namely the regular or chaotic evolution of the system's orbits [START_REF] Lichtenberg | Regular and Chaotic Dynamics[END_REF].

The rapid and efficient detection of the regular or chaotic nature of motion in many dof systems has been a very active research topic over the years, starting with the theoretical introduction of Lyapunov Exponents (LEs) [START_REF] Lyapunov | The general problem of the stability of motion (Taylor and Francis[END_REF][START_REF] Oseledec | [END_REF], and the development of efficient algorithms for their calculation [8,9] (see also [10] for a recent review), which in time expanded and evolved as a subfield of dynamical systems leading to a gallery of chaos detection methods (see e.g. [11][12][13] for comparisons among different methods and [START_REF]Chaos detection and predictability[END_REF] for a collection of recent reviews of several chaos indicators).

In this study, we employ the Generalized Alignment Index (GALI) [15], a chaos detection method which has been shown to be a very efficient and fast tool for such purposes, especially in high dimensional problems (see e.g. [16,17]). GALI's main advantages are (i) the quick distinction between regular and chaotic motion, (ii) the early detection of weak chaotic behavior, (iii) the determination of the dimensionality of the torus on which the quasi-periodic motion occurs, and (iv) the prediction of slow diffusion [16].

The GALI requires the evolution of several deviation vectors from a given orbit. The number k of the vectors used defines the order of the index i.e. GALI k . In the case of regular motion GALI k , for k not greater than the dimension N of the tangent space of the torus on which the motion takes place, eventually oscillates around a positive value, while for chaotic orbits it tends exponentially fast to zero [15]. In the case of stable periodic orbits of Hamiltonian systems, GALI k tends to zero following particular power laws, while for unstable periodic orbits it becomes zero exponentially fast [18]. A detailed discussion of the theory of GALIs, along with their applications to various dynamical systems can be found in [19].

In this paper, we seek out to investigate and understand the GALIs' prediction efficiency and performance when the initial condition (IC) of an orbit moves gradually from a regular motion region to a chaotic one in Hamiltonian systems, by studying in detail the variation of the asymptotic values of the indices. This transition can take place either by changing the total energy of the system or by moving the regular orbit's IC towards the edge of a stability island. In addition, we follow the time evolution of the deviation vectors needed for the computation of the GALIs and examine if different initial distributions of the coordinates of these vectors affect the final GALI value. We perform our investigations in phase space regions around some simple, stable periodic orbits of the Fermi-Pasta-Ulam-Tsingou (FPUT) model [START_REF] Fermi | Studies of the nonlinear problems. I[END_REF][START_REF] Ford | [END_REF], which describes a chain of harmonic oscillators coupled through nonlinear interactions. The dynamics of this system has been studied extensively in the last decades (see e.g. [START_REF]Focus Issue: The 'Fermi-Pasta-Ulam' Problem-The First Fifty Years[END_REF][START_REF]The Fermi-Pasta-Ulam Problem[END_REF] and references therein) and is considered nowadays a prototypical, multidimensional, nonlinear model.

The paper is organized as follows. In Sect. II we present the FPUT β system, discuss the numerical integration of its equations of motion and its variational equations, and briefly present the definition of the GALI method, its properties, as well as its computation. Then, in Sect. III we study in detail the asymptotic behavior of the GALIs for regular orbits in the phase space neighborhood of two types of stable periodic orbits of the FPUT β model considering systems with different dof. In Sect. IV we focus our attention on the properties of the deviation vectors used for the computation of the GALIs, while in Sect. V we summarize our results and discuss our findings.

II. MODEL AND NUMERICAL TECHNIQUES

In our study, we consider the one-dimensional (1D) FPUT β model [START_REF] Fermi | Studies of the nonlinear problems. I[END_REF][START_REF] Ford | [END_REF], which represents a lattice chain of N identical particles with quartic nearest-neighbor interactions. The system's Hamiltonian function is given by

H N = N i=1 1 2 p 2 i + N i=0 1 2 (x i+1 -x i ) 2 + β 4 (x i+1 -x i ) 4 , (1)
where x i is the displacement of the i-th particle from its equilibrium position and p i is the corresponding conjugate momentum. In our investigation we impose fixed boundary conditions to system (1), i.e. x 0 = x N +1 = p 0 = p N +1 = 0 and set β = 1. The time evolution of a phase space orbit with ICs at time t = t 0 z(t 0 ) = ( x(t 0 ), p(t 0 )), where x(t 0 ) = (x 1 (t 0 ), x 2 (t 0 ), . . . , x N (t 0 )) and p(t 0 ) = (p 1 (t 0 ), p 2 (t 0 ), . . . , p N (t 0 )), is governed by the Hamilton equations of motion ẋi = dxi dt = ∂H N ∂pi , ṗi = dpi dt = -∂H N ∂xi , i = 1, 2, . . . , N . In order to compute the GALI k chaos indicator for a particular orbit we have to follow the evolution of the orbit itself, along with a set of k, initially linearly independent, deviation vectors v i , i = 1, 2, . . . , k, corresponding to k different perturbations of the studied orbit. The time evolution of an initial deviation vector v(t 0 ) = δz(t 0 ) = ( δx(t 0 ), δp(t 0 )) = (δx 1 (t 0 ), . . . , δx N (t 0 ), δp 1 (t 0 ), . . . , δp N (t 0 )) from an orbit with ICs z(t 0 ) = ( x(t 0 ), p(t 0 )) is governed by the so-called variational equations (see e.g. [10])

˙ v(t) = J 2N • D 2 H (z(t)) • v(t 0 ), (2)
where J 2N = 0 N I N -I N 0 N , with I N and 0 N being respectively the identity and the zero N × N matrices, and

D 2 H (z(t)) is the 2N × 2N Hessian matrix with elements D 2 H (z(t)) i,j = ∂ 2 H ∂zi∂zj z(t)
, i, j = 1, 2, . . . , 2N .

Since Hamiltonian (1) can be split in two integrable parts H N ( x, p) = A( p) + B( x), with A( p) being the kinetic energy which is a function of only the momenta p i , and B( x) being the potential energy depending only on the coordinates x i , we implement an efficient fourth-order symplectic integration scheme called ABA864 [24][25][26] for integrating the system's equations of motion. Symplectic integrators are numerical schemes specifically designed to preserve the symplectic structure of Hamiltonian systems. One of their main advantages is that they keep the error in the computed value of the Hamiltonian (which is typically referred to as the system's 'energy') bounded for all times. In our numerical simulations, we typically integrate system (1) up to a final time t f = 10 8 . By adequately adjusting the used integration time step we always keep the absolute relative energy error below 10 -8 . By implementing the so-called 'tangent map' technique [27][28][29] we use the ABA864 scheme for integrating also the system's variational equations [START_REF] Wiggins | Global Bifurcations and Chaos: Analytical Methods[END_REF], following in this way the time evolution of the set of deviation vectors needed for the computation of the GALIs.

A. The GALI method

Let us now briefly discuss the GALI method and its properties considering our N dof autonomous Hamiltonian [START_REF] Mackay | Hamiltonian Dynamical Systems[END_REF]. By following an orbit z(t) = ( x(t), p(t)) and a set of k initially independent deviation vectors v i (t), i = 1, 2, . . . , k, with 2 ≤ k ≤ 2N , the Generalized Alignment Index of order k (GALI k ) represents at any time t the volume of the generalized parallelogram having as edges the k unit deviation vectors vi (t) = vi(t) vi(t) [15]. This volume is computed as the norm of the wedge product of these vectors

GALI k (t) = v1 (t) ∧ v2 (t) ∧ . . . ∧ vk (t) . (3)
We note that the number k of the considered deviation vectors should not exceed the dimension 2N of the system's phase space, because in that case the vectors will be by default linearly depended, and the corresponding volume (and consequently the value of GALI k ) will be zero. The GALI is a generalization of the Smaller Alignment Index (SALI) [30][31][32], with GALI 2 being equivalent to SALI [15]. Both indices have been proven to be very efficient chaos indicators and have been successfully implemented in studies of various dynamical systems (see e.g. [19] and references therein).

In an N dof Hamiltonian system regular motion typically occurs on N -dimensional (N D) tori. The system's dynamics lead any deviation vector of the regular orbit to eventually fall on the N D tangent space of the torus [15]. Thus, asymptotically the volume defined by k initially linearly independent deviation vectors, with N < k ≤ 2N , will vanish as these vectors will become linearly dependent, while, in general, for 2 ≤ k ≤ N the volume will not become zero as the vectors will remain linearly independent. The general behavior of the GALI k for regular orbits lying on an N D torus of an N dof Hamiltonian system is given by [15] GALI

k (t) ∝    constant if 2 ≤ k ≤ N t -2(k-N ) if N < k ≤ 2N. ( 4
)
We note that the power law decay in (4) changes in the case of motion on a low-dimensional torus [16] (i.e. a torus with dimensionality lower than N ), and if some of the initial deviation vectors are already on the torus' tangent space [15], but we will not discuss such cases here. Tori of regular motion exist around stable periodic orbits in conservative (or Hamiltonian) systems. For such orbits GALI k always decreases to zero following some specific power laws [18]

GALI k (t) ∝    t -(k-1) if 2 ≤ k ≤ 2N -1 t -2N if k = 2N. (5)
On the other hand, for chaotic orbits and for unstable periodic orbits GALI k decays exponentially fast to zero [15,18]

GALI k ∝ e -[(σ1-σ2)+(σ1-σ3)+•••+(σ1-σ k )]t , (6)
where

σ 1 ≥ σ 2 ≥ • • • ≥ σ k are
A = v1 v2 • • • vk =      v1,1 v2,1 • • • vk,1 v1,2 v2,2 • • • vk,2 . . . . . . . . . v1,2N v2,2N • • • vk,2N      (7)
having as columns the coordinates of the k unitary vectors vi (t) = vi(t) vi(t) = (v i,1 , vi,2 , . . . , vi,2N ) [16]. In particular

GALI k = k i=1 z i (8)
where z i , i = 1, . . . , k, are the so-called singular values of A, obtained through the SVD procedure. We compute the value of GALI k by following the time evolution of k initially independent random unit deviation vectors v1 (0), v2 (0), . . . , vk (0). In order to statistically analyze the behavior of GALIs we average the indices over several different choices of the set of initial deviation vectors. The random choice of the initial vectors leads to different GALI k (0) values. Thus, in order to fairly and adequately compare the behavior of the indices for different initial sets of vectors we normalize the GALIs evolution by registering the ratio GALI k (t)/GALI k (0), i.e. we measure the change of the volume defined by the k deviation vectors with respect to the initially defined volume. Another option is to start the evolution of the dynamics by considering a set of k orthonormal vectors so that GALI k (0) = 1. As we will see later on, both approaches lead to similar results, so in our study we will follow the former procedure, unless otherwise stated.

In order to illustrate the basic behavior of the GALIs for regular and chaotic orbits we present in Fig. 1(a 1) with N = 5 and total energy H5 = 5. All the other coordinates of both orbits are initially set to zero. The straight lines in (a) correspond to functions proportional to t -2 , t -4 , t -6 , t -8 and t -10 , as indicated in the panel's legend. Both axes of (a) and the vertical axis of (b) are logarithmic.

[chaotic] orbit of Hamiltonian (1) with N = 5. The regular orbit of Fig. 1(a) lies on a 5D torus and consequently GALI 2 , GALI 3 , GALI 4 and GALI 5 eventually oscillate around some constant non-zero value, which decreases with increasing k. One of our main aims in this work is to analyze in depth this behavior. Note that the saturation of the indices to their practically constant limiting values happens at latter times for higher k values. The GALIs with k ≥ 6 eventually tend to zero following asymptotic power laws, which are indicated by straight lines in Fig. 1(a). All these behaviors are in accordance to Eq. ( 4). On the other hand, as expected from (6) the GALIs decay exponentially fast to zero in the case of the chaotic orbit of Fig. 1(b). For the sake of completeness, in Fig. 1(b) we also demonstrate the behavior of GALIs for chaotic motion. Since we focus here on the behavior of GALIs for regular orbits, we do not investigate further the case of Fig. 1(b), by e.g. computing the corresponding LEs. Such investigations were systematically performed, for example, in [15,16].

III. THE BEHAVIOR OF THE GALI FOR REGULAR ORBITS

In general, regular motion occurs in the vicinity of stable periodic orbits in Hamiltonian systems. Thus, in order to study the behavior of the GALIs for regular orbits we first locate some stable periodic orbits of system [START_REF] Mackay | Hamiltonian Dynamical Systems[END_REF] and then compute the GALIs for orbits in the neighborhood of these orbits. The stability type of a periodic orbit is determined by the eigenvalues of the so-called monodromy matrix M(T ), which is obtained from the solution of the variational equations of the periodic orbit for one period T (for more details, explicit equations and explanations see e.g. [33] and references therein). The monodromy matrix is symplectic, and its columns are linearly independent solutions of the equations that govern the evolution of deviation vectors from the periodic orbit. If all the eigenvalues of M(T ) are on the unit circle in the complex plane the corresponding periodic orbit is stable, while otherwise it is unstable. We note that there exist several different types of instabilities [30], but we will not discuss this issue here.

In an N D autonomous Hamiltonian system, two eigenvalues are always equal to λ = 1 [30], which means that in practice the remaining 2(N -1) eigenvalues define the periodic orbit's stability. Thus, we can reduce our investigation to a 2(N -1)D subspace of the whole phase space through the well-known method of the Poincaré Surface of Section (PSS) (see e.g. [START_REF] Lichtenberg | Regular and Chaotic Dynamics[END_REF]), where the corresponding monodromy matrix has 2(N -1) eigenvalues, none of which is by default

λ = 1.
In what follows we investigate the behavior of the GALIs for regular orbits in the neighborhood of two simple periodic orbits (SPOs) of system (1), which we refer to as SPO1 and SPO2. The dynamics of these orbits was discussed in [34,35].

A. Regular motion in the neighborhood of SPO1

The first SPO we study is called SPO1 in [35] and it is obtained by considering an FPUT β lattice (1) with N being an odd integer, so that all particles at evennumbered positions are kept stationary at all times, while the odd-numbered particles are always displaced symmetrically to each other, i.e.

x2j (t) = 0, x2j-1 (t) = -x 2j+1 (t) = x(t), (9)
for j = 1, 2, . . . , N -1 2 . By inserting conditions (9) in the system's equations of motion we end up with a second order nonlinear differential equation for variable x(t) which describes the oscillations of all moving particles of SPO1. In particular, for j = 1, 3, 5, . . . , N we have

ẍj (t) = -2x(t) -2β x3 (t), (10)
while xj (t) = 0, j = 2, 4, 6, . . . , (N -1) for all times. The stability analysis of the SPO1 [34,35] showed that for small values of H N (1) the orbit is stable, but it becomes unstable when the energy increases beyond a certain threshold H c N . More specifically, in [34] it was shown that the energy density threshold H c N /N decreases with increasing N following an asymptotic law H c N /N ∝ N -1 (see for example Fig. 1 of [34]). In Fig. 2 we see the arrangement of the eigenvalues of the monodromy matrix of the SPO1 orbit of Hamiltonian (1) with N = 11 when H N /N = 0.1 [Fig. 2(a)] and H N /N = 0.2 [Fig. 2(b)]. In the first case the SPO1 is stable as all eigenvalues are on the unit circle, while in the latter the orbit is unstable as two eigenvalues are off the unit circle. We note that the transition from stability to instability for the SPO1 happens at H c N /N ≈ 0.1755.

To investigate how GALIs of different orders k behave for regular orbits we perturb the stable SPO1 with ICs xi (0), pi (0), i = 1, 2, . . . , N and energy density H N /N to obtain a nearby orbit with ICs X i (0), P i (0), ensuring that the new orbit is at the same energy density level. The phase space distance D between these two orbits in this case is given by

D = N i=1 (x i (0) -X i (0)) 2 + (p i (0) -P i (0)) 2 1/2 . ( 11
)
Having this regular orbit lying on a N D torus means that GALI k with 2 ≤ k ≤ N will asymptomatically approach some almost constant value [see Eq. ( 4)]. In order to statistically analyze these asymptotic GALI k values we follow the evolution of GALIs for n v = 10 different random sets of initial deviation vectors whose coordinates are drawn from a uniform distribution in the interval [-0.5, 0.5], and compute the evolution of the average (over the n v sets of deviation vectors) values GALI k (t)/GALI k (0) . The integration of the orbit and of the deviation vectors is performed until the GALI values show a clear saturation. Then we estimate the asymptotic GALI value by finding the mean value of GALI k (t)/GALI k (0) over the last n t = 20 recorded values, captured during, approximately, the one-fifth of the last decade of integration. We denote this quantity as GALI k and estimate its error through the standard deviation of the considered n t values.

The outcome of this analysis for a regular orbit with distance (11) D = 0.12 from the stable SPO1 and H N /N = 0.01 is seen in Fig. 3. In particular, Fig. 3(a) depicts the evolution of GALI k (t)/GALI k (0) for some selected orders k (black curves).The gray area around these curves represents one standard deviation. We observe that GALIs of higher orders converge to lower values and need more time to settle to these values. In Fig. 3(b) we see more clearly how these final asymptotic values, GALI k , decrease with increasing k.

Let us now investigate how these behaviors change when the studied regular orbit is taken further and further away from the stable SPO1, i.e. as we increase D. A representative case is presented in Fig. 4 where we plot the evolution of GALI 2 [Fig. 4(a)] and GALI 4 [Fig. 4(b)] for several orbits in the neighborhood of the stable SPO1 with H N /N = 0.01 and D 1 = 0.008, D 2 = 0.01, D 3 = 0.06, D 4 = 0.1, D 5 = 0.22 and D 6 = 0.4, when one set of random initial deviation vectors is used. First of all, we observe that the behavior of the GALIs for the regular orbit closest to the SPO1 (D 1 = 0.008) follows the same power laws as the periodic orbit itself, namely (5), i.e. GALI 2 ∝ t -1 and GALI 2 ∝ t -3 . As D increases, i.e. the IC of the regular orbit moves further away from the stable SPO1 inside the stability island surrounding the periodic orbit, GALIs start deviating from the power law decay observed for D = D 1 , and finally saturate to a positive value. This deviation starts earlier for larger D values and consequently the asymptotic values of GALIs increase as D grows. Eventually, for very large values of D, in the particular case discussed in Fig. 4 for D = D 6 = 0.4, the perturbed orbit becomes chaotic as it is located outside the stability island. In that case, the GALIs decrease exponentially fast to zero in accordance with [START_REF] Lyapunov | The general problem of the stability of motion (Taylor and Francis[END_REF].

In Fig. 4 we show the increase of only the asymptotic values of GALI 2 and GALI 4 with respect to D. This behavior is clearly seen for the whole spectrum of GALIs in Fig. 5 where we plot the GALI k values versus k for regular orbits around the SPO1 orbit with H N /N = 0.01 having D values which are relatively not too small, in order to avoid the power law decays observed in Fig. 4 for D = 0.008, but also relatively not too large in order for the studied orbit to be well inside the stability island.

So far we saw how the asymptotic GALI values change when we increase the distance D of the studied orbit from the stable SPO1 for constant energy H N (and obviously constant energy density H N /N ). In order to see how GALI k values change when we increase the orbit's energy density, we begin by estimating the size of the stability island around the SPO1 orbit with different H N /N values by finding the largest D value for which regular motion is observed. We denote this value D m . The outcome of this process is shown in Fig. 6. The D m value (and consequently the size of the stability island) decreases as H N /N grows, and vanishes when the SPO1 periodic orbit destabilizes and becomes unstable. Recall that in this case the destabilization of the SPO1 takes place for H c N /N ≈ 0.1755. Thus, the region below the curve of Fig. 6 corresponds to regular motion.

In Fig. 7 we present the values of GALI k as a function of the GALI's order k for regular orbits in the neighborhood of the stable SPO1 for increasing energy densities E = H N /N , i.e. moving the orbit's IC in the region below the curve of Fig. 6 towards the destabilization energy of SPO1. In particular, we consider regular orbits with D = 0.1 for E 1 = 0.01, E 2 = 0.04 and with D = 0.01 for E 3 = 0.12 and E 4 = 0.16. As the energy density increases the regular region in Fig. 6 is shrinking, thus the considered regular orbit is located closer and closer to the stable SPO1 orbit itself and the D value should also decrease in order to avoid transition to chaotic motion. Since the regular orbit approaches the stable SPO1 the GALI k values shown in Fig. 7 decrease for orbits close to the periodic one, in accordance to what we observed in Fig. 5.

The behavior we observe in the N = 11 case for the asymptotic GALI values of regular orbits located further and further away from the stable SPO1 (Fig. 5), as well as the one we see for regular orbits whose energies approach the destabilization energy of the SPO1 (Fig. 7), remain the same also for other values of N . As a testimony to that we present in Fig. 8 results for the dependence of the GALI k values on the index's order k for the SPO1 orbit of system [START_REF] Mackay | Hamiltonian Dynamical Systems[END_REF] with N = 21. Note that in this case the destabilization energy density of the SPO1 is H N /N ≈ 0.0675.

B. Regular motion in the neighborhood of SPO2

The findings of Sect. III A are not restricted to only one type of SPO but are rather general. To demonstrate this we perform a similar analysis for regular orbits in the neighborhood of what was called SPO2 in [35]. This periodic orbit appears in FPUT β systems (1) with N = 5 + 3m, m = 0, 1, 2, . . ., particles where every third particle remains always stationary and the two particles in between move in opposite directions Similarly to the case of SPO1 inserting conditions ( 12) and ( 13) in the equations of motion of system (1) leads to a single differential equation

x3j (t) = 0, j = 1, 2, 3, . . . , N -2 3 , (12)
ẍj (t) = -3x(t) -9β x3 (t), ( 14
)
for the moving particles, while xj (t) = 0 for the stationary particles with j = 3, 6, 9, . . . , N -2. In [35], where the stability of this orbit was studied in detail, it was found that the destabilization energy density threshold H c N /N decreases when the number of particles grows, as H c N /N ∝ N -2 (see Fig. 2(a) of [35]).

For N = 11 the stability region of SPO2 is shown in Fig. 9, which can be directly compared with the stability region of SPO1 in Fig. 6, as both cases have the same number of particles. The shapes of both stability regions look similar but the extend of the SPO2 region in the energy density axis is much smaller as the stable SPO2 becomes unstable for H c N /N ≈ 0.01395. In Fig. 10 we present results regarding the asymptotic GALI values, GALI k , for regular orbits inside the stability island around the SPO2 orbit for N = 11 [Figs. 10(a) and (b)] and for N = 20 [Figs. 10(c) and(d)]. We observe similar behaviors to the ones encountered in the neighborhood of SPO1. In particular, the GALI k values increase as we move away from the stable periodic orbit [Figs. 10(a) and(c)], similarly to what we observed in Figs. 5 and8(a), while they decrease as the energy of the regular orbit approaches the destabilization energy of SPO2 [Figs. 10(b) and(d)] as we saw in Figs. 7 and8(b) for the SPO1 case. We note that for N = 20 the destabilization threshold is in this case H c N /N ≈ 0.00425.

IV. STATISTICAL ANALYSIS OF DEVIATION VECTORS

In the case of chaotic orbits it is known that all, initially different, deviation vectors will eventually become aligned to the direction defined by the maximum LE [8][9][10]. This behavior leads to the eventual vanishing of the GALIs through the exponential decay (6) [15]. On the other hand, for regular orbits we know that deviation vectors eventually fall on the tangent space of the torus on which the motion occurs [15], but it is not clear how this behavior affect the actual GALI values. For this reason we investigate here in detail the behavior of deviation vectors for regular orbits. An important finding is that the evolution of the GALIs does not practically depend on the initial choice of deviation vectors needed for their computation. A particular example illustrating this property is shown in Fig. 11 where we present the evolution of GALIs for a regular orbit close to the stable SPO1 (D = 0.12) of Hamiltonian (1) for N = 11 and H N /N = 0.01. In Fig. 11(a) we show the evolution of GALI k (t)/GALI k (0) for some selected values of k when the coordinates of the unit initial deviation vectors are chosen from a uniform distribution in the interval [-0.5, 0.5] (black curves) or from a normal distribution with mean 0 and standard deviation 1 (light black curves). The two curves practically coincide. The presented results are averaged over n v = 10 random sets of initial deviation vectors and the gray area around the curves indicate one standard deviation, as was done for example in Fig. 3(a). In Fig. 11(b) we show a similar computation, but now the initial deviation vectors, whose initial coordinates were created in a similar way as in Fig. 11(a), were orthonormalized at the beginning of their evolution. This choice sets GALI k (0) = 1 and for this reason we report in Fig. 11(b) the average value GALI k (t) . Also for this case we see that the method of creation of the initial, orthonormal deviation vectors does not affect the evolution of the GALIs. Furthermore, by comparing Figs. 11( In Fig. 12 we plot the probability density distributions of the coordinates of the unit deviation vectors needed for the evaluation of GALI 11 for a regular orbit close to SPO1 of Hamiltonian [START_REF] Mackay | Hamiltonian Dynamical Systems[END_REF] with N = 11. In particular, we consider the orbit with H N /N = 0.01, D = 0.12 in Fig. 12(a) and H N /N = 0.002, D = 0.04 in Fig. 12(b). The distributions are created by the coordinates of the set of 11 deviation vectors obtained at 10 snapshots when the GALI 11 has reached its asymptotic value. Black curves correspond to initially unit deviation vectors, while gray curves to orthonormal ones. Both curves practically overlap.

The results of Figs. 11 and12 clearly indicate that the initial distribution of the deviation vector coordinates does not play a role in the evolution of GALIs. Thus, we will continue our investigations by considering only initial unit deviation vectors having their coordinates generated from a uniform distribution.

The distribution of the deviation vectors coordinates evolve in time as we see in Fig. 13 where we divide into five time intervals the evolution of GALI 11 for a regular orbit close to the stable SPO1 of Hamiltonian (1) for N = 11, with H N /N = 0.004 and D = 0.04 [Fig. 13(a)]. Each interval has the same length in logarithmic scale. The corresponding coordinate distributions are shown in Fig. 13(b), while in Fig. 13(c) we present similar distributions but for another, random set of initial deviation vectors. From Figs. 13(b) and (c) we see the evolution of the distributions from a rather spread 'triangular' shape to a more concentrated one, along with clear similarities between the distributions obtained by different initial sets of deviation vectors. Thus, in general, not only the final distribution obtained when GALIs have reached their asymptotic behavior, but also the time evolution of the coordinate distributions is practically independent of the initial choice of deviation vectors.

In Figs. 14 and 15 we investigate the dependence of the final deviation vector coordinate distribution on the regular orbit's distance D from the stable SPO in the case of SPO1 [Fig. 14(a)] and SPO2 [Fig. 15(a)]. We see that in both cases we have a transition from a rather pointy distribution with high concentrations in the distribution's middle and edges (small D values) to a more 'triangular' shape as we approach the boundaries of the stability island (large D values). In accordance to this behavior we see in Figs. 14(b) and 15(b) that when the orbit's energy density increases towards the SPO's destabilization energy the distributions move to a more pointy and concentrated shape, because the regular orbit is located closer to the stable SPO, as we have already pointed out in the discussion of Fig. 7.

V. SUMMARY AND DISCUSSION

In this work we investigated in detail the behavior of the GALI chaos indicator for regular motion of multidimensional Hamiltonian systems. Thus, our study completes in some sense previous works on the GALI method [15,16,18] where other aspects of the index were investigated, like for example its behavior for periodic orbits [18]. In particular, we considered several regular orbits in the vicinity of two basic SPOs of the FPUT β model (1), the so-called SPO1 and SPO2 orbits [35], for various numbers of the system's dof. We showed that the time evolution of the GALIs, along with the distribution of the coordinates of the deviation vectors needed for their computation, are not influenced by the initial choice of the deviation vectors (Figs. 11 and 12). In order to do so we considered various sets of initial deviation vectors. More specifically, we used unit vectors whose coordinates were drawn randomly from a uniform or a normal distribution. In that case, since the initial value GALI k (0) is different for each choice of vectors, we registered the evolution of GALI k (t)/GALI k (0) . An al-ternative option we considered was vectors, whose coordinates were again randomly generated from a uniform or a normal distribution, which were also orthonormalized at the beginning of our numerical simulations, setting in this way GALI k (0) = 1. For this choice we followed the evolution of GALI k (t) .

The value of GALI k eventually saturates to a positive value for regular motion on an N D torus for 2 ≤ k ≤ N . The asymptotic GALI k values depend on (i) the order k of the index, i.e. GALIs of higher order attain lower values, (ii) the phase space distance D (11) of the regular orbit from the nearby SPO, i.e. the GALI k values increase when D grows, and (iii) the orbit's energy H N (1), or equivalently the energy density H N /N , i.e. as we approach the destabilization energy of the SPO the GALI k values decrease. Furthermore, we showed that also the shape of the coordinate distributions of the deviation vectors depend on D and H N /N (Figs. 14 and 15), although these distributions are not influenced by the initial choice of the vectors.

8 GALI 9 GALI 10 tFIG. 1 .

 89101 FIG. 1. The time evolution of the normalized GALI indices GALI k (t)/GALI k (0), k = 2, 3, . . . , 10 for (a) a regular orbit with ICs x1 = -1.03003, x3 = -x5 = 1.04003, p1 = 0.29284, and (b) a chaotic orbit with ICs x1 = 1.00097, x3 = -x5 = 1.04003, p1 = 0.57060, of system (1) with N = 5 and total energy H5 = 5. All the other coordinates of both orbits are initially set to zero. The straight lines in (a) correspond to functions proportional to t -2 , t -4 , t -6 , t -8 and t -10 , as indicated in the panel's legend. Both axes of (a) and the vertical axis of (b) are logarithmic.

FIG. 2 .

 2 FIG. 2. The arrangement on the complex plane of the eigenvalues λi, i = 1, 2, . . . , 20, of the monodromy matrix of the SPO1 (9) of Hamiltonian (1) with N = 11 for (a) the stable SPO1 with HN /N = 0.1 and x(0) = -0.4112, and (b) the unstable SPO1 with HN /N = 0.2 and x(0) = -0.5626. The monodromy matrix is evaluated on the PSS (a) x1 = -0.411, p1 > 0, and (b) x1 = -0.560, p1 > 0. The critical energy density for which the SPO1 encounters its first transition from stability to instability is H c N /N ≈ 0.176.

  FIG. 3.(a) The time evolution of the average, over nv = 10 sets of random initial deviation vectors, GALI k (t)/GALI k (0) (black curves) for k =2, 3, 5, 7, 9 and 11, for a regular orbit in the vicinity of the stable SPO1 with HN /N = 0.01 and D = 0.12. The shaded gray area around each curve denotes 1 standard deviation. (b) Estimation of the asymptotic GALI values, GALI k , as a function of their order k for the orbit of panel (a). The error bars denote 1 standard deviation.

3 FIG. 4 . 4 FIG. 5 .

 3445 FIG. 4. The time evolution of (a) GALI2(t)/GALI2(0) and (b) GALI4(t)/GALI4(0) for orbits with HN /N = 0.01 and distances D1 = 0.008, D2 = 0.01, D3 = 0.06, D4 = 0.1, D5 = 0.22 and D6 = 0.4 from the stable SPO1, for 1 set of initial deviation vectors. The straight lines correspond to functions proportional to (a) t -1 , and (b) t -3 .

FIG. 6 .

 6 FIG. 6. Estimation of the extent in phase space of the stability island around SPO1: The maximum value of D, denoted by Dm, for which regular motion occurs, as a function of the energy density HN /N for the Hamiltonian (1) with N = 11.

4 FIG. 7 . 4 FIG. 8 .

 4748 FIG. 7. Similar to Fig. 5 but for regular orbits with different energy densities E = HN /N : E1 = 0.01, E2 = 0.04, E3 = 0.12 and E4 = 0.16. The phase space distance D (11) is D = 0.1 for E1 and E2, and D = 0.01 for E3 and E4.

FIG. 9 .

 9 FIG. 9. Similar to Fig. 6 but for the SPO2 of Hamiltonian (1) with N = 11.

4 FIG. 10 .

 410 FIG. 10. Estimation of the asymptotic GALI values, GALI k , as a function of their order k for regular orbits in the vicinity of the stable SPO2 of system (1) with N = 11 [(a) and (b)] and N = 20 [(c) and (d)]. The considered regular orbits have (a) E = HN /N = 0.001 and distances D (11) from the stable SPO2 D1 = 0.01, D2 = 0.04, D3 = 0.08 and D4 = 0.12, (c) E = HN /N = 0.002 with D1 = 0.01, D2 = 0.04, D3 = 0.08 and D4 = 0.12, (b) E1 = 0.001 and D = 0.12, E2 = 0.004 and D = 0.12, E3 = 0.008 and D = 0.04, and E4 = 0.012 and D = 0.01, (d) E1 = 0.0004 and D = 0.1, E2 = 0.0008 and D = 0.04, E3 = 0.002 and D = 0.04, and E4 = 0.004 and D = 0.002.

FIG. 11 .

 11 FIG. 11. The time evolution of the average over nv = 10 sets of initial deviation vectors of (a) GALI k (t)/GALI k (0) and (b) GALI k (t) for k = 2, 3, 5, 7, 9, 11, for a regular orbit in the neighborhood of the stable SPO1 (D = 0.12, HN /N = 0.01) of Hamiltonian (1) with N = 11. The coordinates of the initially (a) unit [(b) orthonormalized] deviation vectors were chosen from a uniform (black curves) and a normal distribution (light black curves). The 2 curves practically coincide in panels (a) and (b), while the gray area around them indicate 1 standard deviation. (c) Similar to Fig. 3(b) but for the 4 different cases of panels (a) and (b). The cases are denoted in the legend by the panels name, (a) or (b), and by 'u' or 'n' for uniform and normal distributions. The 4 curves practically overlap.

  a) and (b) we observe that also the choice of unit or orthonormal initial deviation vectors does not have any practical influence on the evolution of the indices. This becomes more evident from the results of Fig.11(c) where we plot the estimations of the asymptotic GALI values, GALI k , obtained from the results of Figs.11(a) and (b) [in the same way that Fig.3(b)was obtained from Fig.3(a)], as the four curves practically overlap.

FIG. 12 .

 12 FIG.12. Probability density distributions of the coordinates vi of 11 initially unit (black curves -'u') or orthonormalized (gray curves -'o') deviation vectors for a regular orbit of Hamiltonian[START_REF] Mackay | Hamiltonian Dynamical Systems[END_REF] with N = 11 close to the stable SPO1, with (a) HN /N = 0.01, D = 0.12, and (b) HN /N = 0.002, D = 0.04. The distributions are created from data obtained for 10 snapshots of the deviation evolution when the GALI11 values have practically reached their asymptotic behavior. In both panels the two curves practically overlap.

5 FIG. 13

 513 FIG. 13. (a) The time evolution of GALI11 for a regular orbit in the neighborhood of the stable SPO1 orbit of Hamiltonian (1) for N = 11 with HN /N = 0.004 and D = 0.04. The evolution is divided into 5 intervals, I1: 0 ≤ log 10 t < 1.5, I2: 1.5 ≤ log 10 t < 3, I3: 3 ≤ log 10 t < 4.5, I4: 4.5 ≤ log 10 t < 6 and I5: 6 ≤ log 10 t ≤ 7.5. The coordinate distributions of the 11 unit deviation vectors for the intervals I1, I2, I3, I4 and I5 are shown in (b) for the deviation vectors used in (a), while in (c) we see results similar to the ones of (b) but for another random set of initial deviation vectors.

4 FIG. 14 . 4 FIG. 15 .

 414415 FIG. 14. The final coordinate distributions of the 11 unit deviation vectors used for the computation of GALI11 for regular orbits in the neighborhood of the stable SPO1 of Hamiltonian (1) with N = 11. The considered regular orbits have (a) E = HN /N = 0.01 and distances D (11) from the stable SPO1 D1 = 0.01, D2 = 0.06, D3 = 0.12 and D4 = 0.22, and (b) E1 = 0.001 and D = 0.015, E2 = 0.01 and D = 0.015, E3 = 0.1 and D = 0.015, and E4 = 0.175 and D = 0.001.

ACKNOWLEDGMENTS

H. M. was partially funded by the University of Cape Town (UCT) International and Refugee Grant, as well as by the Woldia University and MoSHE. Ch. S. acknowledges support by the UCT's Research Committee (URC). We also thank the High Performance Computing facility of UCT (http://hpc.uct.ac.za) and the Center for High Performance Computing (https://www.chpc.ac. za) for providing the computational resources needed for this work.