Marie Cottrell

Madalina Olteanu

Fabrice Rossi

Joseph Rynkiewicz

Nathalie Villa-Vialaneix

Neural Networks for Complex Data

come

Introduction

In many real world applications of machine learning and related techniques, the raw data are not anymore in a standard and simple tabular format in which each object is described by a common and fixed set of numerical attributes. This standard vector model, while useful and efficient, has some obvious limitations: it is limited to numerical attributes, it cannot handle objects with non uniform descriptions (e.g., situations in which some objects have a richer description than others), relations between objects (e.g., persons involved in a social network), etc.

In addition, it is quite common for real world applications to have some dynamic aspect in the sense that the data under study are the results of a temporal process. Then, the traditional hypothesis of statistical independence between observations does not hold anymore: new hypothesis and theoretical analysis are needed to justify the mathematical soundness of the machine learning methods in this context.

Artificial neural networks provide some of the most efficient techniques for machine learning and data mining [START_REF] Osowski | Mlp and svm networks -a comparative study[END_REF]. As other solutions, they were mainly developed to handle vector data and analyzed theoretically in the context of statistically independent observations. However, the last decade has seen numerous efforts to overcome those two limitations [START_REF] Hammer | Neural methods for non-standard data[END_REF]. We survey in this article some of the resulting solutions. We will focus our attention on the two major artificial neural network models: the Multi-Layer Perceptron (MLP) and the Self-Organizing Map (SOM).

Multi-Layer Perceptrons

The Multi-Layer Perceptron (MLP) is one the most well known artificial neural network model (see e.g., [START_REF] Bishop | Neural Networks for Pattern Recognition[END_REF]). On a statistical point of view, MLP can be considered as a parametric family of regression functions. Technically, if the data set consists in vector observations in R p , that is if each object is described by a vector x = (x 1 , • • • , x p ) T , the output of a one hidden layer perceptron with k hidden neurons is given by

F θ (x) = β + k ∑ i=1 a i ψ w T i x + b i , (1)
where the w i are vectors of R p , and the β , a i and b i are real numbers (θ denote the vector of all parameters obtained by concatenating the w i , a i and b i ). In this equation, ψ is a bounded transfer function which introduces some non linearity in F θ . Given a set of training examples, that is N pairs (X i ,Y i ), the learning process consists in minimizing over θ a distance between Y i the target value and F θ (X i ) the predicted value. Given an error criterion (such as the mean squared error), an optimal value for θ is determined by any optimization algorithm (such as quasi Newton methods see e.g. [START_REF] Boyd | Convex Optimization[END_REF]), leveraging the well know backpropagation algorithm [START_REF] Werbos | Beyond regession: New tools for prediction and analysis in the behavior sciences[END_REF] which enables a fast computation of the derivatives of F with respect to θ . The use of a one hidden layer perceptron model is motivated by approximation results such as [START_REF] Hornik | Multilayer feedforward networks are universal approximators[END_REF] and by learnability results such as [START_REF] White | Connectionist nonparametric regression: mutilayer feedforward networks can learn arbitrary mappings[END_REF] (in statistical community, learning is called estimation and learnability consistency).

Model selection issues for MLP

It is well known since the seminal paper of [START_REF] Lapedes | Nonlinear signal processing using neural networks: Prediction and system modeling[END_REF] that MLP are an efficient solution for modeling time series whenever the linear model proves to be inadequate. The simplest approach consists in building a non linear auto-regressive model: given a real valued time series (Y t ) t∈N , one builds training pairs

Z t = (U t ,Y t ), where U t is a vector in R p defined by U t = (Y t-1 , • • • ,Y t-p )).
Then a MLP is used to learn the mapping between the U t (the past of the time series in a time window of length p) and Y t (the current value of the time series), as in any regression problem.

In order to avoid overlearning and/or large computation time, the question of selecting the correct number of neurons or, more generally, the question of model selection arises immediately. Standard methods used by the neural-networks community are based on pruning: one trains a possibly too large MLP and then removes useless neurons and/or connection weights. Heuristic solutions include Optimal Brain Damage [32] and Optimal Brain Surgeon [START_REF] Hassibi | Optimal brain surgeon and general network pruning[END_REF], but a statistically founded method, SSM (Statistical Stepwise Method), was introduced by [START_REF] Cottrell | Neural modeling for time series: a statistical stepwise method for weight elimination[END_REF]. The method relies on the minimization of the Bayesian Information Criterion (BIC). Shortly after, [START_REF] Yao | On least square estimation for stable nonlinear ar processes[END_REF] and [START_REF] Rynkiewicz | Modèles de réseaux de neurones pour l'analyse des séries temporelles ou la régression[END_REF] proved the consistency (almost surely) of BIC in the case of MLPs with one hidden layer. These results, established for time series, allow to generalize the consistency results in [START_REF] White | Connectionist nonparametric regression: mutilayer feedforward networks can learn arbitrary mappings[END_REF] for the iid case.

The convergence properties of BIC may be generalized even further. A first extension is given in [START_REF] Rynkiewicz | Asymptotic law of likelihood ratio for multilayer perceptron models[END_REF]. The noise is supposed to be Gaussian and the transfer function ψ is supposed to be bounded and three times derivable. Then [START_REF] Rynkiewicz | Asymptotic law of likelihood ratio for multilayer perceptron models[END_REF] shows that under some mild hypothesis, the maximum of the likelihood-ratio test statistic (LRTS) converges toward the maximum of the square of a Gaussian process indexed by a class of limit score functions. The theorem establishes the tightness of the likelihood-ratio test statistic and, in particular, the consistency of penalized likelihood criteria such as BIC. Some practical applications of such methods can be found in [START_REF] Mangeas | Neural model selection: how to determine the fittest criterion[END_REF]. The hypothesis on the noise was relaxed in [START_REF] Rynkiewicz | General bound of overfitting for mlp regression models[END_REF]. The noise is no longer supposed to be Gaussian, but only to admit exponential moments. Under this more general assumption, BIC criterion is still consistent (in probability).

On the basis of the theoretical results above, a practical procedure for MLP identification is proposed. For a one hidden layer perceptron with k hidden units, we first introduce

T n (k) = min θ (E n (θ ) + a n (k, θ )) ,
where E n (θ ) is the mean squared error of the MLP for parameter θ and a n (k, θ ) is a penalty term. Then we proceed as follows:

1. Determination of the right number of hidden units.

(a) begin with one hidden unit, compute T n (1), (b) add one hidden unit if T n (k + 1) ≤ T n (k), (c) if T n (k + 1) > T n (k) then stop and keep k hidden units for the model. 2. Prune the weights of the MLP using classical techniques like SSM [START_REF] Cottrell | Neural modeling for time series: a statistical stepwise method for weight elimination[END_REF].

Note that the choice of the penalty term a n (k, θ ) is very important. On simulated data, good results have been reported for

a n (k, θ ) from a n (k, θ ) = E n (θ ) log(n) n to a n (k, θ ) = E n (θ ) √ (n) n
(see [START_REF] Rynkiewicz | Modèles de réseaux de neurones pour l'analyse des séries temporelles ou la régression[END_REF], [START_REF] Rynkiewicz | Consistent estimation of the architecture of multilayer perceptrons[END_REF]). Let us also mention that the tightness of the LRTS and, in particular, the consistency of the BIC criterion were recently established for more complex neural-networks models such as mixtures of MLPs [START_REF] Olteanu | Estimating the number of components in a mixture of multilayer perceptrons[END_REF] and mixtures of experts [START_REF] Olteanu | Asymptotic properties of mixtureof-experts models[END_REF].

Modeling and forecasting nonstationary time series

As mentioned in the previous section, MLP are a useful tool for modeling time series. However, most of the results cited above are available for iid data or for stationary time series. In order to deal with highly nonlinear or nonstationary time series, a hybrid model involving hidden Markov models (HMM) and multilayer perceptrons (MLP hereafter) was proposed in [START_REF] Rynkiewicz | Hybrid hmm/mlp models for time series prediction[END_REF]. Let us consider (X t ) t∈N a homogeneous Markov chain valued in a finite state-space E = {e 1 , • • • , e N } and (Y t ) t∈N the observed time series. The hybrid HMM/MLP model can be written as follows:

Y t+1 = F X t+1 (Y t , • • • ,Y t-p+1 ) + σ X t+1 ε t+1 , (2)
where

F X t+1 ∈ {F e 1 , • • • , F e N }
is a regression function of order p. In this case, F e i is the i-th MLP of the model, parameterized by the weight vector

w i . σ X t+1 ∈ {σ e 1 , • • • , σ e N } is a strictly positive number and (ε t ) t∈N is a iid sequence of standard Gaussian variables.
The estimating procedure as well as the statistical properties of the parameter estimates were established in [START_REF] Rynkiewicz | Estimation of hybrid hmm/mlp models[END_REF]. The proposed model was successfully applied in modeling difficult data sets such as ozone peaks [START_REF] Dutot | A 24-h forecast of ozone peaks and exceedance levels using neural classifiers and weather predictions[END_REF] or financial shocks [START_REF] Maillet | Nonlinear analysis of shocks when financial markets are subject to changes in regime[END_REF].

Functional data

The original MLP model is limited to vector data for an obvious reason: each neuron computes its output as a non linear transformation ψ applied to a (shifted) inner product w T x + b (see equation ( 1)). However, as first pointed out in [START_REF] Sandberg | Notes on weighted norms and network approximation of functionals[END_REF], this general formula applies to any data space on which linear forms can be defined: give a data space X and a set of linear functions W from X to R, one can define a general neuron with the help of w ∈ W , as calculating ψ(w(x) + b).

This generalization is particularly suitable for functional data, that is for data in which each object is described by one or several functions [START_REF] Ramsay | Functional Data Analysis[END_REF]. This type of data is quite common for instance in multiple time series setting (where each object under study evolves through time and is described by the temporal evolutions of its characteristics) or in spectrometry. A functional neuron [START_REF] Rossi | Functional multi-layer perceptron: a nonlinear tool for functional data anlysis[END_REF] can then be defined as calculating ψ (b + f w dµ), where f is the observed function and w is a parameter function. Results in [START_REF] Rossi | Functional multi-layer perceptron: a nonlinear tool for functional data anlysis[END_REF] show that MLP based on this type of neurons share many of the interesting properties of classical MLP, from the universal approximation to statistical consistency (see also [START_REF] Rossi | Theoretical properties of projection based multilayer perceptrons with functional inputs[END_REF] for an alternative functional neuron with similar properties). In addition, the parameter functions w can be represented by standard numerical MLP, leading to a hierarchical solution in which a top level MLP for functional data is obtained by using a numerical MLP in each of its functional neurons. Experimental results in [START_REF] Rossi | Functional multi-layer perceptron: a nonlinear tool for functional data anlysis[END_REF][START_REF] Rossi | Functional data analysis with multi layer perceptrons[END_REF] show the practical relevance of this technique.

Self-Organizing Maps

As the MLP, Kohonen's Self-Organizing Map (SOM) is one of the most well known artificial neural network model [START_REF] Kohonen | Self-Organizing Maps[END_REF]. The SOM is a clustering and visualization model in which a set of vector observations in R p is mapped to set of M neurons organized in a low dimensional prior structure, mainly a two dimensional grid or a one dimensional string. Each neuron c is associated to a codebook vector p c in R p (p c is also called a prototype). As in all prototype based clustering methods, each p c represents the data points that have been assigned to the corresponding neuron, in the sense that p c is close to those points (according to the Euclidean distance in R p ). The distinctive feature of the SOM is that each prototype p c is also somewhat representative of data points assigned to other neurons, based on the geometry of the prior structure: if neurons c and d are neighbours in the prior structure, then p c will be close to data points assigned to neuron d (and vice versa). On the contrary, if c and d are far away from each other in the prior structure, the data points assigned to one neuron will not influence the prototype of the other neuron. This has some very important consequences in terms of visualization capabilities, as illustrated in [START_REF] Vesanto | Som-based data visualization methods[END_REF] for instance.

The original SOM algorithm has been designed for vector data, but numerous adaptations to more complex data have been proposed. We survey here three specific extensions, respectively to time series, functional data and categorical data. Another important extension not covered here is proposed in [START_REF] Hammer | A general framework for unsupervised processing of structured data[END_REF] which is built upon processing of multiple time series with recursive versions of the SOM. The authors show that trees and graphs can be clustered by those versions of the SOM, using a temporal coding of the structure. Recent advances in this line of research include e.g. [START_REF] Hagenbuchner | Graph selforganizing maps for cyclic and unbounded graphs[END_REF]. Other specific adaptation include the symbol strings SOM described in [START_REF] Somervuo | Online algorithm for the self-organizing map of symbol strings[END_REF].

Time series with metadata

While the SOM is a clustering algorithm, it has been used frequently in supervised context as a component of a complex model. We described briefly here one such model as an example of complex time series processing with the SOM. Let us consider a time series with two time scales, i.e., that can be written down with two subscripts. The date is denoted by ( j, h) where j represents the slow time scale and corresponds for instance to the day (or month or year) while h = 1, . . . , H corresponds to the observed values (e.g. the hours or half-hours of the day, the days of the month, the months of the year, etc.). Then the time series is denoted (c j ) j≥0 = (c j,1 , . . . , c j,H ) j≥0 . We assume in addition that the slow time scale is associated with metadata. For instance, if each j corresponds to a day in a year and one knows the day of the week, the month, etc. Metadata are supposed to be available prior a prediction.

The original time series c j,h takes value in R, but the dual time scale leads naturally to a vector valued time series representation, that is to the c j ∈ R H . In this point of view, given the past of the vector valued time series, one has to predict a future vector value, that is a complete vector of H values. This could be seen as a long term forecasting problem for which a usual solution would be to iterate one-step ahead forecasts. However, this leads generally to unsatisfactory solutions either because of a squashing behaviour (convergence of the forecasting to the mean value of series) or to a chaotic behaviour (for nonlinear methods).

An alternative solution is explored in [START_REF] Cottrell | Forecasting of Curves Using Kohonen Classification[END_REF]. It consists in forecasting separately, on the one hand, the mean and variance of the time series on next slow time scale step (that is, on the next j), and on the other hand, the profile of the fast time scale. The prediction of the mean and of the variance is done by any classical technique. For the profile, a SOM is used as follows. The vector values of the time series, i.e., the (c j ) j≥0 , are centred and normalized with respect to the fast time scale, that is are transformed into profiles defined by

q j = 1 σ j (c j,1 -µ j , . . . , c j,H -µ j ) , (3)
where µ j = 1 H ∑ H h=1 c j,h and σ 2 j = 1 H ∑ H h=1 (c j,hµ j ) 2 are respectively the mean and the variance of c j . The profiles are clustered with a SOM leading to some prototype profiles p c . Each prototype is associated to the metadata of the profiles that has been assigned to the corresponding neuron.

Then a vector value is predicted as follows: the mean µ and variance σ are obtained by a standard forecasting model for the slow time scale. Then the metadata of the vector to predict is matched against the metadata associated to neurons: assume for instance, that metadata are days of the week, and the we try to predict a Sunday. Then one collects all the neurons to which Sunday profiles have been assigned. Finally, a weighted average of the matching prototypes is computed and rescaled according to µ and σ . As shown in [START_REF] Cottrell | Forecasting of Curves Using Kohonen Classification[END_REF] this technique enables both some stable and meaningful full day predictions, while integrating non numerical metadata.

Functional data

The dual time scale approach described in the previous section has become a standard way of dealing with time series in a functional way, as shown in e.g. [START_REF] Besse | Autoregressive forecasting of some functional climatic variations[END_REF]. But as pointed out in Section 2.3, functional data arise naturally in other contexts such as spectrometry. Then, the SOM has been naturally adapted to functional data in other contexts than time series. In those contexts, in addition to the normalization technique described above that produces profiles, one can use functional transformation such as derivative calculations in order to drive the clustering process by the shapes of the functions rather than mainly by their average values [START_REF] Rossi | Clustering functional data with the SOM algorithm[END_REF].

Another adaptation consists in integrating the SOM with optimal segmentation techniques that represent functions or time series with simple models, such as piecewise constant functions for instance. The main idea it to a apply a SOM to functional data using any functional distance (from the L 2 norm to more advanced Sobolev norms [START_REF] Villmann | Sobolev metrics for learning of functional data -mathematical and theoretical aspects[END_REF]) with an additional constraint that prototypes must be simple, e.g., piecewise constant. This leads to interesting visualization capabilities in which the complexity of the display is automatically globally adjusted [START_REF] Hébrail | Exploratory analysis of functional data via clustering and optimal segmentation[END_REF].

Categorical data

In surveys, it is quite standard that the collected answers are categorical variables with a finite number of possible values. In this case, a specific adaptation of the SOM algorithm can be defined, in the same way that Multiple Correspondence Analysis is related to Principal Component Analysis. More precisely, useful encoding methods for categorical data are the Burt Table (BT), which is the full contingency table between all pairs of categories of the variables, or the Complete Disjunctive Table (CDT), that contains the answers of each individual coded as 0/1 against dummy variables that correspond to all the categories of all variables. Then, a Multiple Correspondence Analysis of the BT or of the CDT is nothing else than a Principal Component Analysis on BT or CDT, previously transformed to take into account a specific distance between the rows and a weighting of the individuals [START_REF] Roux | Geometric Data Analysis: From Correspondence Analysis to Structured Data Analysis[END_REF]. The SOM can be adapted to categorical data using this approach, as described in [START_REF] Cottrell | How to use the Kohonen algorithm to simultaneously analyse individuals in a survey[END_REF] and [START_REF] Cottrell | SOM-based algorithms for qualitative variables[END_REF]. The same transformation on BT or CDT is achieved and a SOM using the rows of the transformed tables can thus be trained. This training provides an organized clustering of all the possible values of the categorical variables on a prior structure such as a two dimensional grid. Moreover, if a simultaneous representation of the individuals and of the values is needed, two coupled SOM can be trained and superimposed. The aforementioned articles present various real-world use cases from socio-economic field.

Kernel and dissimilarity SOM

The extensions of artificial neural networks model described in the previous sections are ad hoc in the sense that they are constructed using specific features of the data at hand. This is a strength but also a limitation as they are not universal: given a new data type, one has to design a new adaptation of the general technique. In the present section, we present more general versions of the SOM that are based on a dissimilarity or a kernel on the input data. Assuming the existence of such a measure is far weaker than assuming the data are in a vector format. For instance, it is simple to define a dissimilarity/similarity between the vertices of a graph, a data structure that is very frequent in real world problems [START_REF] Newman | The structure and function of complex networks[END_REF], while representing directly those vertices as vectors is generally difficult.

Dissimilarity SOM

Let us assume that the data under study belong to a set X on which a dissimilarity d is defined: d is a function from X × X to R + that maps a pair of objects x and y to a non negative real number which measures how different x and y are. Hypothesis on d are minimal: it has to by symmetric (d(x, y) = d(y, x)) and such that d(x, x) = 0.

As pointed out above, dissimilarities are readily available on sets of non vector data. A classical example is the string edit distance [START_REF] Levenshtein | Binary codes capable of correcting deletions, insertions and reversals[END_REF] which defines a distance1 on symbol strings. More general edit distances can be defined, such as for instance the graph edit distance which measure distances between graphs [START_REF] Bunke | Graph-based tools for data mining and machine learning[END_REF].

As the hypothesis on X are minimal, one cannot assume anymore that vector calculation are possible in this set. Then, the learning rules of the SOM do not apply as they are based on linear combination of the prototypes with the data points. To circumvent this difficulty, [START_REF] Kohohen | Self-Organizing maps of symbol strings[END_REF] suggest to chose the values of the prototypes p c in the set of observations (X i ) i . This leads to a batch version of the SOM which proceeds as follows. After a random initialization of the prototypes, each observation is assigned to the neuron with the closest propotype (according to the dissimilarity measure) and the prototypes are then updated. For each neuron, the updated p c is chosen among the observations as the minimizer the following distortion

∑ i Γ (N(X i ), c)d(X i , p) (4)
where N(X i ) is X i 's neuron and Γ is a decreasing function of the distance between neurons in the prior structure. This modification of the SOM algorithm is known as the median SOM and is closely related to the earlier median version of the standard k-means algorithm [START_REF] Kaufman | Clustering by means of medoids[END_REF].

In the case where (X i ) i is a small sample, the constraint to chose the prototypes in the data can be seen as too strong. Then, [START_REF] El Golli | Une adaptation des cartes auto-organisatrices pour des données décrites par un tableau de dissimilarités[END_REF] suggests to associate several prototypes (a given number q) to each neuron. A neuron is represented by a subset of size q from (X i ) and the different steps of the SOM algorithm are modified accordingly. A fast implementation is described in [START_REF] Rossi | Fast algorithm and implementation of dissimilarity self-organizing maps[END_REF].

A successful application of the dissimilarity SOM on real world data concerns school-to-work transitions. In [START_REF] Massoni | Career-path analysis using optimal matching and self-organizing maps[END_REF], we were interested in identifying career-path typologies, which is a challenging topic for the economists working on the labor market. The data was issued from the "Generation'98" survey by the CEREQ. The data sample contained information about 16040 young people having graduated in 1998 and monitored during 94 months after having left school. The labor-market statuses had nine categories, from permanent contracts to unemployed and including military service, inactivity or higher education.

The dissimilarity matrix was computed using optimal matching distances [START_REF] Abbott | Sequence analysis and optimal matching methods in sociology[END_REF], which are currently the main stream in economy and sociology. The most striking opposition appeared between the career-paths leading to stableemployment situations and the "chaotic" ones. The stable positions were mainly situated in the west region of the map. However, the north and south regions were quite different: in the north-west region, the access to a permanent contract (red) was achieved after a fixed-term contract (orange), while the south-west classes were only subject to transitions through military service (purple) or education (pink). The stability of the career paths was getting worse as we moved to the east of the map. In the north-east region, the initial fixed-term contract was getting longer until becoming precarious, while the south-east region was characterized by the excluding trajectories: unemployment (light blue) and inactivity (dark blue).

Two other extensions of the SOM to dissimilarity data have been proposed; they both avoid the use of constrained prototypes. The oldest one is based on deterministic annealing [START_REF] Graepel | Self-organizing maps: generalizations and new optimization techniques[END_REF] while a more recent one uses the so-called relational Fig. 1 Career-path visualization with the dissimilarity SOM [START_REF] Massoni | Career-path analysis using optimal matching and self-organizing maps[END_REF]: colors correspond to the nine different categories approach that relies on pseudo-Euclidean spaces [START_REF] Hammer | Topographic mapping of large dissimilarity data sets[END_REF][START_REF] Hammer | Topographic processing of relational data[END_REF]. Both approaches lead to better results for datasets where the ratio between the number of observations and the number of neurons is small.

Kernel SOM

An alternative approach to dissimilarities is to rely on kernels. Kernels can be seen as a generalization of the notion of similarity. More precisely, a kernel on a set X is a symmetric function K from X × X to R that satisfies a positivity property:

∀N ∈ N * , ∀(x i ) 1≤i≤N ∈ X N , ∀(α i ) 1≤i≤N ∈ R N , N ∑ i, j=1 α i α j K(x i , x j ) ≥ 0.
For such a kernel, there is a Hilbert space H (called the feature space of the kernel) and a mapping φ from X , such that the inner product in H corresponds to the kernel via the mapping, that is [START_REF] Aronszajn | Theory of reproducing kernels[END_REF] :

φ (x), φ (x ) H = K(x, x ). (5)
Then K can be interpreted as a similarity on X (values close to zero correspond to unrelated objects) and defines indirectly a distance between objects in X as follows:

d K (x, x ) = φ (x) -φ (x ) H = K(x, x) + K(x , x ) -2K(x, x ). (6)
As shown in e.g. [START_REF] Shawe-Taylor | Kernel methods for pattern analysis[END_REF], kernels are a very convenient way to extend standard machine learning methods to arbitrary spaces. Indeed, the feature space H comes with the same elementary operations as R p : linear combination, inner product, norm and distance. Then, one has just to work in the feature space as if it were the original data space. The only difficulty comes from the fact that φ and H are not explicit in general, mainly because H is an infinite dimensional functional space. Then one has to rely on equation ( 5) to implement a machine learning algorithm in H completely indirectly using only K. This is the so called Kernel trick.

In the case of the batch version of the SOM, this is quite simple [START_REF] Boulet | Batch kernel SOM and related laplacian methods for social network analysis[END_REF]. Indeed, assignments of data points to neurons are based on the Euclidean distance in the classical numerical case: this translates directly into the distance in the feature space, which is calculated solely using the kernel (see equation ( 6)). Prototypes update is performed as weighted averages of all data points: weights are computed with the Γ function introduced in equation ( 4) as a proxy for the prior structure. It can be shown that those weights, which are computed using the assignments only, are sufficient to define the prototypes and that they can be plugged into the distance calculation, without needing an explicit calculation of φ . Variants of this scheme, especially stochastic ones, have been studied in [START_REF] Andras | Kernel-Kohonen networks[END_REF][START_REF] Donald | The kernel self organising map[END_REF]. It should also be noted that the relational approach mentioned in the previous section [START_REF] Hammer | Topographic mapping of large dissimilarity data sets[END_REF][START_REF] Hammer | Topographic processing of relational data[END_REF] can be seen a relaxed kernel SOM, that is an application of a similar algorithm in situations where the function K is not positive.

While kernels are very convenient, the positivity conditions might seem very strong at first. It is indeed much stronger than the conditions imposed to a dissimilarity, for instance. Nevertheless, numerous kernels have been defined on complex data [START_REF] Gärtner | Kernel for Structured Data[END_REF], ranging from kernels on strings based on substrings [START_REF] Lodhi | Text classification using string kernels[END_REF] to kernel between the vertices of a graph such as the heat kernel [START_REF] Kondor | Diffusion kernels on graphs and other discrete structures[END_REF][START_REF] Smola | Kernels and regularization on graphs[END_REF] (see [START_REF] Boulet | Batch kernel SOM and related laplacian methods for social network analysis[END_REF] for a SOM based application of this kernel to a medieval data set of notarial acts). Two graphs can also be compared via a kernel based on random walks [17] or on subtrees comparisons [START_REF] Ramon | Expressivity versus efficiency of graph kernels[END_REF].

Conclusion

Present days data are becoming more and more complex, according to several criteria: structure (from simple vector data to relational data mixing a network structure with categorical and numerical descriptions), time evolution (from a fixed snapshot of the data to ever changing dynamical data) and volume (from small datasets with a handful of variables and one thousand of objects to terabytes and more datasets). Adapting artificial neural networks to those new data is a continuous challenge which can be solved only by mixing different strategies as outlined in this paper: adding complexity to the models enable to tackle non standard behavior (such as non-stationarity), theoretical guarantees limit the risk of overfitting, new models can be tailor made for some specific data structures such as graph or functions, while generic kernel/dissimilarity models can handle almost any type of data. The ability to combine all those strategies demonstrates once again the flexibility of the artificial neural network paradigm.

A distance is a dissimilarity that satisfies in addition the strong hypothesis of the triangle inequality: d(x, y) ≤ d(x, z) + d(z, y).


C M Verrelli
email: verrelli@ing.uniroma2.it

G Damm

Robust Transient Stabilization Problem for a Synchronous Generator in a Power Network

Keywords: Robust nonlinear control, adaptive control, L 2 and L ∞ disturbance attenuation, large-scale systems, power systems control

The robust transient stabilization problem (with stability proof) of a synchronous generator in an uncertain power network with transfer conductances is rigorously formulated and solved. The generator angular speed and electrical power are required to be kept close, when mechanical and electrical perturbations occur, to the synchronous speed and mechanical input power, respectively, while the generator terminal voltage is to be regulated, when perturbations are removed, to its pre-fault reference constant value. A robust adaptive nonlinear feedback control algorithm is designed on the basis of a third order model of the synchronous machine: only two system parameters (synchronous machine damping and inertia constants) along with upper and lower bounds on the remaining uncertain ones are supposed to be known. The conditions to be satisfied by the remote network dynamics for guaranteeing L2 and L∞ robustness and asymptotic relative speed and voltage regulation to zero are weaker than those required by the single machine-infinite bus approximation: dynamic interactions between the local deviations of the generator states from the corresponding equilibrium values and the remote generators states are allowed.

Introduction

Power networks are among the most complex large-scale, interconnected nonlinear systems (see [START_REF] Bergen | Power systems analysis[END_REF], [START_REF] Lu | Nonlinear control systems and power system dynamics[END_REF] for the development of dynamical models of increasing complexity). They have continuously increased in size, power and number of components, but with a large acceleration in the last two decades. A rather difficult control problem, which attracted in the past the interest of the control community and has recently become a major concern as the dramatic blackouts in North America, Italy and Germany confirm, is represented by the power systems transient stabilization: dynamic nonlinear systems interconnected through a large-scale network are required to be controlled, on the basis of partial measurements of the local state variables only, in spite of time varying perturbations of almost all system parameters. In particular, the transient stabilization problem consists in the design of suitable excitation feedback controls such that: i) the angular speed and electrical power of each generator are kept close to the synchronous speed and mechanical input power, respectively, when mechanical and electrical perturbations, such as load shedding, generation tripping or short circuits, occur; ii) the terminal voltage of each generator is regulated to its pre-fault reference constant value when perturbations are removed. Decentralized linear controllers were first designed on the basis of linear approximations around operating conditions. This strategy was motivated by the fact that power systems were relatively small since production and consumption areas were not far. However, this situation has largely changed in the last two decades since the small power systems have been interconnected in order to increase performance and security: large amounts of electric power fluxes are transmitted from one region to another to cope with timevarying consumption and production levels. As a consequence, the overall system complexity has largely increased and the behaviour of such scattered plants has become unpredictable with disregarded nonlinear network phenomena arising: linear controllers may be no longer able to guarantee the power grid stability (see [START_REF] Guo | Global transient stability and voltage regulation for power systems[END_REF]) and to handle the severe disturbances and contingencies typically occurring in power networks. In the recent years, several nonlinear algorithms have been proposed for power systems control. For a particular power systems structure M consisting of a group of generators tied together by a strong network of transmission lines and linked to a single generator g m by a comparatively weak set of tie lines, the well-known single machine-infinite bus approximation, which models all the remaining network as a fixed voltage source and an impedance, is advantageous in the nonlinear control design: transient stabilization of the generator g m can be achieved, with a stability proof, even without requiring the knowledge of critical parameters ( [START_REF] Guo | Global transient stability and voltage regulation for power systems[END_REF], [START_REF] Damm | Adaptive nonlinear excitation control of synchronous generators[END_REF], [START_REF] Damm | Adaptive nonlinear output feedback for transient stabilization and voltage regulation of power generators with unknown parameters[END_REF], [START_REF] Damm | Transient stabilization and voltage regulation of a synchronous generator[END_REF], [START_REF] Marino | Adaptive nonlinear excitation control of synchronous generators with unknown mechanical power[END_REF], [8], [START_REF] Shen | Adaptive nonlinear excitation control with L 2 disturbance attenuation for power systems[END_REF], [START_REF] Tan | Augmentation of transient stability using a superconducting coil and adaptive nonlinear control[END_REF], [START_REF] Wang | Robust nonlinear coordinated control of power systems[END_REF], [START_REF] Wang | Transient stabilization of power systems with an adaptive control law[END_REF]). In particular, L 2 and L ∞ robustness and asymptotic relative speed regulation to zero are guaranteed in [8] despite uncertainties in all system parameters. In this paper, we rigorously formulate and solve the transient stabilization problem (with stability proof) for the particular power systems structure M above considered. The single machine-infinite bus model, used in [START_REF] Guo | Global transient stability and voltage regulation for power systems[END_REF], [START_REF] Damm | Adaptive nonlinear excitation control of synchronous generators[END_REF], [START_REF] Damm | Adaptive nonlinear output feedback for transient stabilization and voltage regulation of power generators with unknown parameters[END_REF], [START_REF] Damm | Transient stabilization and voltage regulation of a synchronous generator[END_REF], [START_REF] Marino | Adaptive nonlinear excitation control of synchronous generators with unknown mechanical power[END_REF], [8], [START_REF] Shen | Adaptive nonlinear excitation control with L 2 disturbance attenuation for power systems[END_REF], [START_REF] Tan | Augmentation of transient stability using a superconducting coil and adaptive nonlinear control[END_REF], [START_REF] Wang | Robust nonlinear coordinated control of power systems[END_REF], [START_REF] Wang | Transient stabilization of power systems with an adaptive control law[END_REF] -which neglects the transient behaviour of the other generators as well as the interconnections between them, does not capture the typical multivariable nature with nonlinear complex coupling of power systems (see for instance the inter-area oscillations) and does not take into account the effect of the generator g m dynamics on the remote nework dynamics -is not used since interactions between the dynamics of the generator g m and the remote network machines are allowed: the transient behaviour of the remote generators states is allowed to depend on the local deviations of the generator states from the corresponding equilibrium values. Following the theoretical developments in [START_REF] Marino | Robust adaptive state-feedback tracking control for nonlinear systems[END_REF], [START_REF] Marino | Nonlinear output feedback tracking with almost disturbance decoupling[END_REF], [START_REF] Shen | Robust H ∞ control of uncertain nonlinear systems via state feedback[END_REF] (even though they do not apply to the model considered in this paper in which an uncertain function multiplies the control input), a robust adaptive nonlinear feedback control is designed for the generator g m which does not assume the knowledge of the overall system parameters excepting for the machine damping and inertia constants. On the basis of upper and lower bounds on the uncertain model parameters, L 2 and L ∞ robustness and transient stabilization are guaranteed under a set of assumptions on the network dynamics which are weaker than those required by the single machine-infinite bus approximation. The presented approach in conjunction with the special structure of power systems considered here allows us to rigorously formulate and solve the robust transient stabilization problem (with L 2 and L ∞ robustness) of a synchronous generator in a power network in the presence of transfer conductances and uncertainties in almost all system parameters under possible dynamic interactions between the local deviations of the generator states from the corresponding equilibrium values and the remote generators states: this constitutes the main original contribution of the paper with respect to the results obtained in the literature for more general structures of power systems (see for instance [START_REF] Shen | Adaptive nonlinear excitation control with L 2 disturbance attenuation for power systems[END_REF], [START_REF] Guo | Nonlinear decentralized control of large-scale power systems[END_REF], [START_REF] Lu | Nonlinear decentralized disturbance attenuation excitation control via new recursive design for multi-machine power systems[END_REF], [START_REF] Machowski | Decentralized stabilityenhancing control of synchronous generator[END_REF], [START_REF] Ortega | Transient stabilization of multimachine power systems with nontrivial transfer conductances[END_REF], [START_REF] Wang | Dissipative Hamiltonian realization and energy-based L 2 -disturbance attenuation control of multimachine power systems[END_REF], [START_REF] Wang | Adaptive H ∞ excitation control of multimachine power systems via the Hamiltonian function method[END_REF], [START_REF] Wang | Robust decentralized nonlinear controller design for multimachine power systems[END_REF], [START_REF] Xi | Nonlinear decentralized controller design for multimachine power systems using Hamiltonian function method[END_REF], [START_REF] Xi | Nonlinear decentralized saturated controller design for power systems[END_REF], [START_REF] Zecevic | Robust decentralized exciter control with linear feedback[END_REF]). Simulation results with reference to a general 3-machine, 9-buses power network show that the proposed robust nonlinear excitation control prevents each network machine from going out of step in the presence of electrical parameter perturbations and unmodelled dynamics and improves the performance with respect to the controller in [8] based on the single machine-infinite bus approximation.

System dynamic model

A power system consisting of n generators interconnected through a transmission network is described by the 3n-order nonlinear models in [START_REF] Sauer | Power system dynamics and stability[END_REF],

1 ≤ i ≤ n, δi = ω i ωi = - D i 2H i ω i + ω 0 2H i P mi - ω 0 2H i P ei (1)
Ė′ qi = k ci T ′ d0i u f i - E ′ qi T ′ d0i - (x di -x ′ di ) T ′ d0i I di
in which the first two equations represent the i-th generator mechanical dynamics involving the power angle δ i (rad), the relative angular speed ω i (rad/s), the active electrical power P ei (p.u.), the mechanical input power P mi (p.u.), the synchronous speed ω 0 (rad/s), the damping constant D i (p.u.) (which is an often neglected non-negative parameter) and the inertia constant H i (s), while the third equation constitutes the i-th generator electrical dynamics involving the transient EMF E ′ qi (p.u.) in the quadrature axis, the input u f i (p.u.) to the thyristor amplifier, the gain k ci of the excitation amplifier, the direct axis transient open circuit time constant T ′ d0i (s), the direct axis reactance x di (p.u.) and the direct axis transient reactance x ′ di (p.u.). The i-th generator electrical equations are

P ei = E ′2 qi G ii + E ′ qi n j=1,j =i E ′ qj G ij cos (δ ij ) + E ′ qj B ij sin (δ ij ) = v di I di + v qi I qi Q ei = -E ′2 qi B ii + E ′ qi n j=1,j =i E ′ qj G ij sin (δ ij ) -E ′ qj B ij cos (δ ij ) -x ′ di I 2 di + I 2 qi = v qi I di -v di I qi I di = -E ′ qi B ii + n j=1,j =i E ′ qj G ij sin (δ ij ) -E ′ qj B ij cos (δ ij ) (2)
I qi = E ′ qi G ii + n j=1,j =i E ′ qj G ij cos (δ ij ) + E ′ qj B ij sin (δ ij ) δ ij = δ i -δ j V ti = x ′ di I qi 2 + E ′ qi -x ′ di I di 2 = v 2 di + v 2
qi in which (for the i-th generator): Q ei (p.u.) is the reactive electrical power, I di (p.u.) is the direct axis current, I qi (p.u.) is the quadrature axis current, V ti (p.u.) is the terminal voltage with (d, q)-components (v di , v qi ), G ij (p.u.) and B ij (p.u.) are the i-th row and the j-th column element of nodal conductance and susceptance matrices, respectively, at the internal nodes after eliminating all physical buses, which depend on x ′ di (p.u.), on the transformer reactance x T i (p.u.), on the loads and on the transmission line reactance x ij (p.u.) between the i-th generator and the j-th generator. The nodal conductance and susceptance matrices represent the full power network of transmission lines and loads connecting the power generators: the network reduction of the load busses makes unsatisfactory the assumption of negligible transfer conductances (see [START_REF] Ortega | Transient stabilization of multimachine power systems with nontrivial transfer conductances[END_REF]), which is typically used in power systems control design. The pre-fault equilibrium, 1 ≤ i ≤ n,

δ i = δ is ω i = 0 E ′ qi = E ′ qi0 with P ei = P mi
guarantees, when all model parameters take their unperturbed values, that the voltages V ti , 1 ≤ i ≤ n, are equal to their corresponding reference constant values V * ti , 1 ≤ i ≤ n. In this paper, we study a particular structure of power systems, i.e., a group G M of n -1 generators, tied together by a strong network of transmission lines, which is linked to a generator g m (in the following referred to as r-th generator) by a comparatively weak set of tie lines. First note that the electrical power P ei affinely appears in the relative angular speed dynamics in (1): it represents an actually measured variable (thus available for feedback as a measured output) which is to be kept close to the mechanical input power during perturbations (thus to be controlled as a controlled output). Therefore it is advantageous to explicitly compute its dynamics in order to obtain a triangular structure for system (1) which is suitable for control design. Let us restrict the analysis to the r-th generator and let us compute [G rj and B rj , 1 ≤ j ≤ n, are assumed to be constant]

Ṗer = Ė′ qr I qr + E ′ qr G rr + E ′ qr n j=1,j =r Ė′ qj G rj cos (δ rj ) + Ė′ qj B rj sin (δ rj ) -E ′ qr n j=1,j =r E ′ qj G rj sin (δ rj ) -E ′ qj B rj cos (δ rj ) ω r -ω j .
We assume that there exist bounded connected open sets D δ and D p such that E ′ qr (0) > 0, δ r (t) ∈ D δ and P er (t) ∈ D p imply I qr (t) > c Ir > 0, which describes the whole practical operating region of the generator (see [START_REF] Shen | Adaptive nonlinear excitation control with L 2 disturbance attenuation for power systems[END_REF]). The r-th generator third order dynamic model can be thus written as

δr = ω r ωr = - D r 2H r ω r + ω 0 2H r P mr - ω 0 2H r P er (3)
u f r is the control input, θ 1r = 1 T ′ d0r , θ 2r = (x dr -x ′ dr ) T ′ d0r , θ 3r = Grr(x dr -x ′ dr ) T ′ d0r , θ 4r = Grr T ′ d0r , θ 5r = kcr T ′ d0r , θ 6r = Grrkcr T ′ d0r and the term R r = n j=1,j =r Ė′ qj G rj cos (δ rj ) + B rj sin (δ rj ) + n j=1,j =r ω j E ′ qj G rj sin (δ rj ) -E ′ qj B rj cos (δ rj )
represents the effect of the network remote dynamics (the group G M of generators) on the r-th generator, where Ẽ′

qj = E ′ qj -E ′ qj0 , 1 ≤ j ≤ n, j = r.
The expression of the function R r (t) will be crucial in the control design since it will allow us to rigorously formulate and solve the robust transient stabilization problem (with L 2 and L ∞ robustness) of a synchronous generator in a power network characterized by the special structure in exam.

In the practice, the exact values of the model parameters are hard to obtain, and in particular P mr , G rr , B rr are lumped parameters which account for unmodelled dynamics such as turbine and load dynamics. Those parameters may undergo sudden on-line variations due to mechanical and electrical perturbations and faults. In the following, we will suppose the parameters ω 0 , D r , H r to be known1 and will assume that: the uncertain piecewise continuous parameters T ′ d0r (t), x dr (t), x ′ dr (t), k cr (t) and the uncertain constant parameter G rr are within the corresponding known positive bounds (T ′ d0rm , T ′ d0rM ), (x drm , x drM ), (x ′ drm , x ′ drM ), (k crm , k crM ), (G rrm , G rrM ); the uncertain constant parameter B rr is within the corresponding known bounds (B rrm , B rrM ); the mechanical input power P mr (t) ∈ D p is a class C 1 function satisfying: P mrM ≥ P mr (t) ≥ P mrm and | Ṗmr (t)| ≤ ṖMr , with P mrm , P mrM , ṖMr known positive reals. Physical considerations concerning transmission lines, loads and mechanical turbines make the above assumptions reasonable in the transient stabilization problem of power networks: upper and lower bounds on the equivalent conductances and susceptances can be computed, while the mechanical input power is usually the output of a slower first order dynamic system (controlled by a higher level control). On the other hand, the knowledge of bounds on the uncertain model parameters somehow restricts the on-line variations which can be tolerated by the controller and this is in line with the local result which will be obtained in the following.

Control problem formulation

A suitable control problem formulation is introduced in this section, which will allow us to solve the transient stabilization problem for the synchronous generator g m by quantitatively characterizing the robustness with respect to both permanent and vanishing model parameter perturbations. Let the pre-fault constant value for the power angle δ r satisfying δ rs ∈ D δ and let θ irm , θ irM be the known positive bounds on the uncertain parameters θ ir (t), 1 ≤ i ≤ 6, with Pmr (t) a suitable estimate of the uncertain mechanical power P mr . We will introduce in the following certain assumptions which are to be satisfied by the remaining part of the network in order to guarantee desired transient and asymptotic closed loop properties: they will impose conditions on the performance achieved by the controllers of the remote network generators (and in particular restrictions on the disturbance caused by remote network dynamics in response to the local deviations of the generator states from the corresponding equilibrium values) and will represent less restrictive conditions than those required by the single machine-infinite bus approximation in which no interaction between the local r-th generator deviations from (δ rδ rs , ω r , P er -P mr ) = 0 and the remote network dynamics is taken into account. They will allow dynamic (and not only static as in [START_REF] Guo | Nonlinear decentralized control of large-scale power systems[END_REF] and [START_REF] Wang | Robust decentralized nonlinear controller design for multimachine power systems[END_REF]) interactions between the local deviations of the generator states from the corresponding equilibrium values and the remote generators states: assumptions imposing only static interactions between the local deviations of the generator states from the corresponding equilibrium values and the remote generators states may not be able to comply with the typical instability phenomena in which the behaviour of generators in the network becomes oscillatory with increasing amplitudes. Let Σ r (y r1 , y r2 , y r3 , R r ) be the system consisting of the j-th generators (j = r, 1 ≤ j ≤ n) with inputs the y r (t) vector components y r1 (t), y r2 (t), y r3 (t) (i.e. the local r-th generator deviations from (δ rδ rs , ω r , P er -P mr ) = 0) and output the function R r (t). The subsequent assumption ii) will be satisfied if the controllers of the remote network generators make the system Σ r (y r1 , y r2 , y r3 , R r ) Input to Output Stable with uniform gain functions φµr (•), φνr (•), φρr (•), while the subsequent assumption iv) will be satisfied if the controllers of the remote network generators force the system Σ r (y r1 , y r2 , y r3 , R r ) to have uniform finite L 2 gains (recall that functions in L 2 represent signals having finite energy over the infinite time interval (t 0 , +∞) and therefore the L 2 gain can be interpreted as the ratio between the energies of the input and the output). Finally, the subsequent condition vi) will impose that, in the case of vanishing perturbations, the pre-fault state values for the remaining part of the network are asymptotically recovered provided that asymptotic regulation (S3) (for the r-th generator) is achieved. The transient stabilization problem addressed in this paper is rigorously formulated as follows.

Definition 1 (Transient Stabilizing Control): Assume that (t ≥ 0): i) for each j = r, 1 ≤ j ≤ n, the j-th generator power angle δ j (t), the relative angular speed ω j (t) and the quadrature axis transient EMF E ′ qj (t) are piecewise differentiable functions of time t and boundedness of the r-th generator variables δ r (t), ω r (t), P er (t) implies boundedness of E ′ qj (t);

ii) there exist µ r , ν r , ρ r (unknown) non-negative reals, φµr (•), φνr (•), φρr (•) known K ∞ functions2 and g r (t) (unknown) bounded non-negative real-valued function of time t such that the following inequality: Define [ε jr , ε xr , ε Br (1 ≤ j ≤ 6) are positive reals]:

|R r (t)| ≤ sup
y r (t) = δ r (t) -δ rs , ω r (t), P er (t) -P mr (t) T ξ r (t) = y r (t) T , P mr (t) -Pmr (t) T w dr (t) = Ṗmr (t), θ 1rM -θ 1rm + ε 1r , θ 2rM -θ 2rm + ε 2r , θ 3rM -θ 3rm + ε 3r , θ 4rM -θ 4rm + ε 4r , max{θ 5rM -θ 5rm + ε 5r , θ 6rM -θ 6rm + ε 6r }, x ′ drM -x ′ drm + ε xr , B rrM -B rrm + ε Br , sup 0≤τ ≤t {g r (τ )}, µ r , ν r , ρ r T .
A bounded piecewise continuous real-valued control law u f r (•) is called a transient stabilizing control for the r-th generator if it guarantees the closed loop system to satisfy the following:

(S1) L ∞ disturbance attenuation property, i.e.

ω 2 r (t) + (P er (t) -P mr (t)) 2 ≤ h 1r (ξ r (0))e -crt + 1 k r γ 1r ( w dr (•) ∞ )
holds for all t ≥ 0, where h 1r (ξ r (0)) ≥ 0, c r > 0 and γ 1r (r) is a class K ∞ function;

(S2) L 2 disturbance attenuation property, i.e. According to the previous definition, the effect of model parameters variations from their unperturbed values is quantified by properties (S1)-(S2), which give bounds on L ∞ and L 2 regulation errors. While the formulation introduced in Definition 1 leads to a closed loop robustness characterization, the adaptive formulation introduced in the following definitions allows us to obtain asymptotic relative speed regulation to zero in the case of definitely constant permanent parameter perturbations and asymptotic voltage regulation in the case of vanishing parameter perturbations.

Definition 2 (Transient Adaptive Stabilizing Control): Assume that, in addition to assumptions i)-ii), there exist non-negative reals t 0 , M Rr , γ δr , γ ωr and γ pr such that:

iii) the uncertain machine parameters P mr (t), T ′ d0r (t), x dr (t), x ′ dr (t), k cr (t) are constant for all t ≥ t 0 ; iv) the following inequality: Definition 3 (Adaptive Regulating Control): Assume that, in addition to assumptions i)-ii), iv) v) there exist a non-negative real t0 such that all the model parameters take their unperturbed values (i.e. all the perturbations are removed) for all t ≥ t0 ; vi) under the asymptotic regulation property (S3), for each j = r,

lim t→+∞ t t0 R 2 r (τ )dτ ≤ M Rr + γ δr lim t→+∞ t t0 [δ r (τ ) -δ rs ] 2 dτ + γ ωr lim t→+∞ t t0 [ω r (τ )]
lim t→+∞ [δ j (t) -δ js ] = 0 lim t→+∞ [ω j (t)] = 0 lim t→+∞ E ′ qj (t) -E ′ qj0 = 0.
A transient stabilizing control u f r (•) is called an adaptive regulating control for the r-th generator if it guarantees, the additional property:

(S4) asymptotic voltage regulation, i.e.

lim t→+∞ [V tr (t) -V * tr ] = 0.

Nonlinear robust design and stability analysis

By virtue of techniques similar to those used in [START_REF] Damm | Transient stabilization and voltage regulation of a synchronous generator[END_REF], we design, in this section, a transient adaptive stabilizing control3 for the r-th generator according to Definitions 1 and 2. The subsequent control design is motivated by the choice of the quadratic function V r and W r which will be introduced in the following. Define the power angle regulation and relative angular speed tracking errors δr = δ rδ rs ωr = ω rω * r with δ rs being the power angle constant reference value and ω * r being the relative angular speed time-varying reference signal

ω * r = - 5 4 k δr δr . (4)
We design the active electrical power time-varying reference signal P * er as

P * er = 2H r ω 0 5 4 k ωr ωr + 5 4 k δr ω r - D r 2H r ω * r + δr + 1 k ωpr ωr + Pmr (5)
which relies on the estimate Pmr of the uncertain mechanical power P mr satisfying the estimation law4

Pmr = φ r + 2H r ω 0 5 4 k per + k r 4 + 1 k r + ω 2 0 16H 2 r k ωpr ω r (6) φr
= 5 4 k per + k r 4 + 1 k r + ω 2 0 16H 2 r k ωpr -φ r + D r ω 0 ω r + P er - 2H r ω 0 5 4 k per + k r 4 + 1 k r + ω 2 0 16H
) > 0] Proj[ζ, ẑr ] =              ζ if z rm ≤ ẑr ≤ z rM ζ if ẑr < z rm and ζ ≥ 0 ζ if ẑr > z rM and ζ ≤ 0 ζ 1 - zrm 2 -ẑ 2 r zrm 2 -(zrm-εzr) 2 if ẑr < z rm and ζ < 0 ζ 1 - ẑ2 r -zrM 2 (zrM +εzr) 2 -zrM 2
if ẑr > z rM and ζ > 0 whose properties are the following:

1. z rm -ε zr ≤ ẑr (t) ≤ z rM + ε zr , ∀t ≥ 0; 2. Proj[ζ, ẑr , •, •, •] is Lipschitz continuous; 3. |Proj[ζ, ẑr , z rm , z rM , ε zr ]| ≤ |ζ|;
0≤τ ≤t {g r (τ )} 2 + µ 2 r + ν 2 r + ρ 2 r
so that, according to (4)-( 5) and property 1. of the projection algorithm, we can establish that the robust adaptive nonlinear feedback control algorithm ( 4)-( 7) is a transient stabilizing control for the r-th generator. Suppose that, in addition to assumptions i)-ii), assumptions iii)-iv) hold and consider the quadratic function

W r = V r + 1 2 6 j=1 β jr θ2 jr + 1 2 β xr x′2 dr + 1 2 β Br B2
ir whose time derivative along the trajectories of the error system (8), for all t ≥ t 0 , is given by

Ẇr = - 5 4 k δr δ2 r - 5 4 k ωr + D r 2H r ω2 r + ω 0 2H r Pmr ωr - 1 k ωpr ω2 r - 5
then, by integrating [START_REF] Tan | Augmentation of transient stability using a superconducting coil and adaptive nonlinear control[END_REF] we obtain

lim t→+∞ t t0 k δr δ2 r (τ ) + k ωr ω2 r (τ ) + k pr P 2 er (τ ) + k per P 2 mr (τ ) dτ ≤ W r (t 0 ) + M Rr k Rr (12)
so that, according to (4)-( 5) and Barbalat's Lemma, we can establish that the robust adaptive nonlinear feedback control algorithm (4)-( 7), ( 9) is a transient adaptive stabilizing control for the r-th generator. Finally, from (2) it is straightforward to establish that, under the additional assumptions v) (implying iii)) and vi), asymptotic voltage regulation (S4) is achieved so that the robust adaptive nonlinear feedback control algorithm (4)-( 7), ( 9) is an adaptive regulating control for the r-th generator. The result holds for any initial condition (of the r-th generator) and positive control parameter k r maintaining, according to (S1), δ r (t) ∈ D δ and P er (t) ∈ D p (guaranteeing I qr (t) > c Ir > 0) for all t ≥ 0.

The main result of this paper can be summarized in the following theorem, which somehow 6 extends the recent theoretical contribution in [8].

Theorem: The robust adaptive nonlinear feedback control algorithm ( 4)-( 7), ( 9) is:

• a transient stabilizing control for the r-th generator;

6 Here the machine damping and inertia constants are assumed to be known.

• a transient adaptive stabilizing and an adaptive regulating control for the r-th generator when (see [START_REF] Wang | Robust nonlinear coordinated control of power systems[END_REF]) for any initial condition (of the r-th generator) and positive control parameter k r maintaining, according to (S1), δ r (t) ∈ D δ and P er (t) ∈ D p (guaranteeing I qr (t) > c Ir > 0 so that singularities in the controller are avoided) for all t ≥ 0.

k
The robust adaptive nonlinear feedback control algorithm ( 4)-( 7), ( 9) relies on: the available signals 7 δ r (t), ω r (t), (I dr (t), I qr (t)) (depending on a common reference for measuring the machine rotor angle), P er (t), Q er (t) (recall [START_REF] Lu | Nonlinear control systems and power system dynamics[END_REF] and see for instance [START_REF] Guo | Nonlinear decentralized control of large-scale power systems[END_REF], [START_REF] Lu | Nonlinear decentralized disturbance attenuation excitation control via new recursive design for multi-machine power systems[END_REF], [START_REF] Wang | Robust decentralized nonlinear controller design for multimachine power systems[END_REF], [START_REF] Xi | Nonlinear decentralized saturated controller design for power systems[END_REF]); the positive control parameters k δr , k ωr , k pr , k per , k ωpr , k r , k Rr , β jr , β xr , β Br , ε jr , ε xr , ε Br (1 ≤ j ≤ 6) whose role may be evaluated by examining both the closed loop error equations and the corresponding functions V r and W r with their time derivatives. The parameters k δr , k ωr , k pr , k per direct affect the dynamics of δr , ωr , Per and Pmr , respectively; the parameters k ωpr , k r , k Rr characterize the robustifying terms in ( 5), ( 6) and ( 7); the parameters β jr , β xr , β Br (1 ≤ j ≤ 6) are the adaptation gains for θjr , x′ dr , Brr (1 ≤ j ≤ 6) and the larger they are chosen the slower the adaptations for θjr , x′ dr , Brr (1 ≤ j ≤ 6) result.

Simulation results

In this section we illustrate the performance and the robustness of the feedback control algorithm ( 4)-( 7), [START_REF] Shen | Adaptive nonlinear excitation control with L 2 disturbance attenuation for power systems[END_REF] in the presence of unmodelled dynamics: the proposed control is applied to each generator of the popular Western System Coordinating Council (WSCC) 3-machine, 9-bus system reported in [START_REF] Sauer | Power system dynamics and stability[END_REF] and [START_REF] Anderson | Power system control and stability[END_REF] [D i = 0, 1 ≤ i ≤ 3] and described by the two-axis model

([26]) [1 ≤ i ≤ 3] δi = ω i ωi = - D i 2H i ω i + ω 0 2H i P mi - ω 0 2H i P ei Ė′ di = - E ′ di T ′ q0i + (x qi -x ′ di ) T ′ q0i I qi Ė′ qi = k ci T ′ d0i u f i - E ′ qi T ′ d0i - (x di -x ′ di ) T ′ d0i I di P ei = E ′ qi I qi + E ′ di I di = v di I di + v qi I qi Q ei = E ′ qi I di -E ′ di I qi -x ′ di I 2 di + I 2 qi = v qi I di -v di I qi I di = n j=1 E ′ qj G ij sin (δ ij ) -E ′ qj B ij cos (δ ij ) + E ′ dj G ij cos (δ ij ) + E ′ dj B ij sin (δ ij ) I qi = n j=1 E ′ qj G ij cos (δ ij ) + E ′ qj B ij sin (δ ij ) -E ′ dj G ij sin (δ ij ) + E ′ dj B ij cos (δ ij ) V ti = E ′ di + x ′ di I qi 2 + E ′ qi -x ′ di I di 2 = v 2 di + v 2
qi from which the model (1) has been derived by neglecting the dynamics of the fast damperwinding E ′ di and by using the simplification x qi = x ′ di with x qi (p.u.) being the quadrature axis reactance. The initial conditions for the state variables are computed by systematically solving the load-flow equations of the network and by computing the values of the algebraic variables. For 1 ≤ i ≤ 3, the functions φµi (•), φνi (•), φρi (•) are set equal to I(•) [0,+∞)

[I(•) is the identity function] i.e. φji (q) = q, j = µ, ν, ρ (q ∈ [0, +∞)), the control parameters are chosen as

k δi = k ωi = k pei = k ωpi = 1, k pi = 720, k i = 0.001, k Ri = 0.1, β 1i = β 2i = β 3i = β 4i = β 5i = β 6i = β xi = β Bi =
s(t) = f (t), if f (t) > s(t -h) s(t -h), otherwise , t ≥ h s(0) = f (0)
with h the integration step and t = lh (l ∈ N 0 ). In order to avoid division by zero 8 , for 1 ≤ i ≤ 3, η(P ei ) and η(I qi ) replace P ei and I qi in the control algorithm ( 4)-( 7), ( 9), respectively, where

η(ξ) = ξ, if ξ ≥ 0.05 0.05, otherwise.
The goal of the simulation is to verify the effects of a three-phase fault occurring near bus 7 at the end of line 5-7 at t = 0.001 s, which is cleared at 0.084 s by opening line 5-7 (circuit breakers reclose at t = 1 s). Figure 1-2 show a satisfactory performance of the proposed control even with the physical saturations (which are hit several times according to Fig. 2(b)): despite the considered severe perturbation, synchronous speeds are quickly restored while regulation of both electrical power and terminal voltage is guaranteed. For comparison, the same simulation is performed by applying (to each generator of the network) the nonlinear robust adaptive control (with control parameters k δ = 0.1, k ω = 3, k p = 120, k = 0.001, β i = 0.0000001, 1 ≤ i ≤ 7 and initial conditions 9 θi (0) = 0 i = 1, 2, 4, 5, θ6 (0) =

u f r0 sin(δr(0)) v2 (0)
, θ3 (0) = P mr , θ7 (0) = 1)

u f r = θ6 sin(δ r ) v2 - k 4 sin(δ r ) Per v2 2 P * er = θ7 v1 + k 4 ωr v2 1 8
Recall that the simulation is carried out by using the fourth order model and by applying the proposed controller to each generator of the network. 9 We denote by u f r0 the nominal input value.

v1 = -θ1 ω r + θ3 + k δ ω r + k 4 ωr ω 2 r + k 4 ωr + k ω ωr + δr θi = β i Proj φ 2i ωr + φ3i Per , θi , θi (0) ∈ [θ im , θ iM ], i = 1, 3, 4, 5 θ2 = β 2 Proj φ32 Per -ωr Per , θi , θ2 (0) ∈ [θ 2m , θ 2M ] θ6 = β 6 Proj -v 2 Per , θi , θ6 (0) ∈ [θ 6m , θ 6M ] θ7 = β 7 Proj v1 ωr , θi , θ7 (0) ∈ [θ 7m , θ 7M ] Proj[ζ, θi ] =                ζ if θ im ≤ θi ≤ θ iM ζ if θi < θ im and ζ ≥ 0 ζ if θi > θ iM and ζ ≤ 0 ζ 1 - θim 2 -θ2 i θim 2 -(θim-ρi) 2
if θi < θ im and ζ < 0 designed in [8] on the basis of a (third order) single machine-infinite bus model. As illustrated by the simulation results reported in Figs. 34, neglecting the transient behaviour of the other generators and the interconnections between them may be critical for the control design in [8]: in the presence of the considered perturbation, synchronous speeds are not quickly restored while electrical power and terminal voltage regulations are not satisfactorily obtained.

ζ 1 -

Conclusions

The transient stabilization problem (with stability proof) has been rigorously formulated in Section 3 and solved in Section 4 for a particular power systems structure M consisting of a group of generators tied together by a strong network of transmission lines and linked to a single generator g m by a comparatively weak set of tie lines. A robust adaptive nonlinear feedback control ( 4)-( 7), ( 9) has been designed for the generator g m which does not assume the knowledge of the overall system parameters excepting for the machine damping and inertia constants: an innovative design technique has been used since available techniques developed in [8], [START_REF] Marino | Robust adaptive state-feedback tracking control for nonlinear systems[END_REF], [START_REF] Marino | Nonlinear output feedback tracking with almost disturbance decoupling[END_REF], [START_REF] Shen | Robust H ∞ control of uncertain nonlinear systems via state feedback[END_REF] do not apply to the model (3). The proposed controller guarantees the L 2 and L ∞ disturbance attenuation and asymptotic regulation properties (S1)-(S4) under assumptions i)-vi) on the network dynamics generalizing those required by the single machine-infinite bus approximation (which does not capture the typical multivariable nature with nonlinear complex coupling of power systems and does not take into account the effect of the generator g m dynamics on the remote nework dynamics) and allowing dynamic interactions between the local deviations of the generator states from the corresponding equilibrium values and the remote generators states (in order to comply with the typical instability phenomena in which the behaviour of generators in the network becomes oscillatory with increasing amplitudes).

In our view the result of this paper constitutes a first step towards rigorous mathematical formulations and solutions to the multimachine transient stabilization problem.

  0≤τ ≤t {g r (τ )} + µ r φµr max 0≤τ ≤t {|δ r (τ )δ rs |} + ν r φνr max 0≤τ ≤t {|ω r (τ )|} +ρ r φρr max 0≤τ ≤t {|P er (τ ) -P mr (τ )|} holds uniformly in u f r (•).

T 0 ω 2 rγ

 2 (τ ) + (P er (τ ) -P mr (τ )) 2 dτ ≤ h 2r (ξ r (02r ( w dr (τ ) )dτ holds for any given T > 0, where h 2r (ξ r (0)) ≥ 0 and γ 2r (r) is a class K ∞ function.

4 .

 4 (z r -ẑr )Proj[ζ, ẑr , z rm , z rM , ε zr ] ≥ (z r -ẑr )ζ [which hold provided that the uncertain constant z r and the initial condition ẑr (0) belong to the compact set [z rm , z rM ]]. Define the estimation errors (1 ≤ j ≤ 6) θjr = θ jr -θjr x′ dr = x ′ dr -x′ dr Brr = B rr -Brr so that, on the basis of the r-th generator dynamics (3), we obtain the error dynamics δr = -5 4 k δr δr + ωr ωr = -5 4

θ2 i -θiM 2 (φ 21 =-

 221 θiM +ρi) 2 -θiM 2 if θi > θ iM and ζ > 0 θ im (θ imρ i ) ≥ 0, θ iM (θ iM + ρ i ) ≥ 0 v2 = -k p Per + θ2 ωr -k 4 ω2 r Per -φ30 --ω r , φ 23 = 1, φ 24 = φ 25 = 0 φ30 = P er ω r cotg(δ r )k δ θ7 v1 + θ1 ω r -θ3 θ7 + k 2 ωr v1φ31 = ω r φs , φ32 = P er φs , φ33 = -φs , φ34 = -P er , φ35 = ω r sin 2 (δ r ) rδ rs , ωr = ω r + k δ δr , Per = P er -P * er

Figure 1 :

 1 Figure 1: Proposed control algorithm [Generator 1 (solid), Generator 2 (dot), Generator 3 (dash)]: a) Power angle regulation errors δ iδ is (1 ≤ i ≤ 3); b) Relative angular speeds ω i (1 ≤ i ≤ 3); c) Active electrical powers P ei (1 ≤ i ≤ 3).

Figure 2 :

 2 Figure 2: Proposed control algorithm [Generator 1 (solid), Generator 2 (dot), Generator 3 (dash)]: a) Terminal voltages V ti (1 ≤ i ≤ 3); b) Control signals u f i (1 ≤ i ≤ 3).

  Ṗer = -θ 1r P erθ 2r I dr I qrθ 3r I dr P er

			I qr	-θ 4r	P 2 er I 2 qr	+	θ 5r I 2 qr + θ 6r P er I qr	u f r
	-Q er + B rr	P 2 er I 2 qr	+ x ′ dr (I 2 dr + I 2 qr ) ω r +	P er I qr	R

r in which δ r , ω r , P er are the state variables,

  We design the robust control law u f r with the robustifying terms v r , π r as rely on the estimates θjr , x′ dr , Brr of the uncertain parameters θ jr , x ′ dr , B rr whose estimation laws (1 ≤ j ≤ 6) [functions ξ θj (t), ξ x (t), ξ b (t) are yet to be chosen] θ1r = Proj ξ θj (t), θjr , θ jrm , θ jrM , ε jr , θ jrm ≤ θjr (0) ≤ θ jrM

	u f r =	I qr qr + θ6r P er θ5r I 2		v r -	k r I qr 2I 4 qr + 2P 2 er v 2 r Per qr + θ6r P er 4 θ5r I 2 2 θ 5rm I 2 qr + θ 6rm P er
	v r = -	5 4	k pr Per +	ω 0 2H r	ωr +	5D r k δr 4ω 0	ω r +	2H r ω 0	ω r + θ1r P er + θ2r I dr I qr + θ3r	I dr P er I qr
	+ θ4r	P 2 er I 2 qr		+ Brr	P 2 er ω r I 2 qr	+ x′ dr I 2 dr + I 2 qr ω r + Q er ω r -	5k δr 4ω 0	D r ω r -ω 0 ( Pmr -P er )
	-	2H r ω 0	5 4	k ωr +	1 k ωpr		5 4	k ωr +	D r 2H r	ωr + δr +	ω 0 2H r	Per +	1 k ωpr	ωr + π r	(7)
	π r = -	k Rr 4	P 2 er I 2 qr		Per -	k r 4	Per P 2 er + I 2 dr I 2 qr +	I 2 dr P 2 er I 2 qr	+	P 4 er qr I 4	1 + ω 2 r + ω 2 r I 2 dr + I 2 qr	2
	-	k r 4		5 4	k ωr +	5 4	k δr +	5 4	k per +	k r 4	+	1 k r	+	ω 2 0 16H 2 r	k ωpr +	1 k ωpr	2 Per
	-	k r 4	P 2 er I 2 qr	Per	φ2 µr max 0≤τ ≤t { δr (τ )} + φ2 νr max 0≤τ ≤t {ω r (τ )} +	5 4	k δr max 0≤τ ≤t { δr (τ )}
	+ φ2 ρr max 0≤τ ≤t { Per (τ )} + P mrM -P mrm +	ṖMr k per	+	2H r ω 0	5 4	k ωr +	1 k ωpr	+	5 4	k δr max 0≤τ ≤t {ω r (τ )}
	+		2H r ω 0	1 +	25 16	k 2 δr +	5D r k δr 4ω 0	max 0≤τ ≤t { δr (τ )} + 1
	which 5 ẋ′ dr		= Proj ξ x (t), x′ dr , x ′ drm , x ′ drM , ε xr , x ′ drm ≤ x′
																	ωr = ̟r +	ω 0 2Hr	Pmr
	in which
																̟r = -	Dr 2Hr	ωr -	ω 0 2Hr	Per
	is available for feedback.

2 r k ωpr ω r P mrm ≤ Pmr (0) ≤ P mrM and introduce the active electrical power tracking error Per = P er -P * er . dr (0) ≤ x ′ drM Ḃrr = Proj ξ b (t), Brr , B rrm , B rrM , ε Br , B rrm ≤ Brr (0) ≤ B rrM will be designed by using the projection algorithm Proj[ζ, ẑr , z rm , z rM , ε zr ] in [8] [(z rmε zr

,

  θ1r P er Per + β 1r θ1r -θ2r I dr I qr Per + β 2r θ2r -θ3r I dr P er Per Brr , B rrm , B rrM , ε Br , B rrm ≤ Brr (0) ≤ B rrM then, by using property 4. of the projection algorithm and by completing the squares, we obtain

	ẋ′ dr	= Proj -	I 2 dr + I 2 qr ω r Per β xr	, x′ dr , x ′ drm , x ′ drM , ε xr , x ′ drm ≤ x′ dr (0) ≤ x ′ drM
	Ḃrr = Proj -	P 2 er ω r Per β Br I 2 qr
			Ẇr ≤ -	5 4	k δr	δ2 r -	5 4	k ωr	ω2 r -	5 4	k pr	P 2 er -	5 4	k per	P 2 mr +	R r (t) 2 k Rr	.	(10)
	4 Since δr , ωr , Per , Pmr are bounded, according to (8) and assumption i), δr , ωr , Ṗ er , Ṗ mr are k pr P 2 er + π r Per bounded so that δr (t), ωr (t), Per (t), Pmr (t) are uniformly continuous for all t ≥ t 0 . On the
	-other hand, by virtue of assumptions iii)-iv), if k r 4 2I 4 qr + 2P 2 er θ 5r I 2 qr + θ 6r P er v 2 r θ5r I 2 qr + θ6r P er 2 θ 5rm I 2 qr + θ 6rm P er P 2 er + k r 4 + 1 k r + ω 2 0 16H 2 r k ωpr + 1 k ωpr Pmr Per -k Rr > max 4γ δr k δr , 4γ ωr k ωr , with	+ k per + P er I qr R r Per -k r 4 + 1 k r 16γ pr 5 4 k pr 16γ pr , k per	5 4 +	k ωr + ω 2 0 16H 2 r	5 4 k ωpr k δr + P 2 5 4 mr k per
		γδr = γ δr + + Ṗmr Pmr -I qr 25 8 γ ωr k 2 δr + 4 2H r ω 0 25 16 k 2 δr + 1 + 5 4 2 D r k δr ω 0 γ pr -x ′ dr I 2 dr + I 2 qr ω r Per + β xr ẋ′ dr -θ4r P 2 er Per I 2 qr + β 4r θ4r -Brr P 2 er ω r Per qr I 2 γωr = 2γ ωr + 4 2H r ω 0 5 4 k ωr + 5 4 k δr + 2 1 k ωpr γ pr	+ β 3r θ3r + β Br Ḃrr
		+ θ5r	 	I 2 qr v r Per θ5r I 2 qr + θ6r P er	-β 5r θ5r	  + θ6r	 	P er v r Per θ5r I 2 qr + θ6r P er	-β 6r θ6r	  .
																	that
	θ1r = Proj -	P er Per β 1r	, θ1r , θ 1rm , θ 1rM , ε 1r , θ 1rm ≤ θ1r (0) ≤ θ 1rM
	θ2r = Proj -	I dr I qr Per β 2r	, θ2r , θ 2rm , θ 2rM , ε 2r , θ 2rm ≤ θ2r (0) ≤ θ 2rM
	θ3r = Proj -	I dr P er Per I qr β 3r		, θ3r , θ 3rm , θ 3rM , ε 3r , θ 3rm ≤ θ3r (0) ≤ θ 3rM
	θ4r = Proj -	P 2 er Per I 2 qr β 4r	, θ4r , θ 4rm , θ 4rM , ε 4r , θ 4rm ≤ θ4r (0) ≤ θ 4rM	(9)
	θ5r = Proj	I 2 qr v r Per θ5r I 2 qr + θ6r P er β 5r	, θ5r , θ 5rm , θ 5rM , ε 5r , θ 5rm ≤ θ5r (0) ≤ θ 5rM
	θ6r = Proj	P er v r Per θ5r I 2 qr + θ6r P er β 6r	, θ6r , θ 6rm , θ 6rM , ε 6r , θ 6rm ≤ θ6r (0) ≤ θ 6rM

If we choose the yet to be defined functions ξ θj (t), ξ x (t), ξ b (t) (1 ≤ j ≤ 6) so

  348000 while the initial conditions for the parameter estimates are set equal to the corresponding unperturbed values. The time function s(t) = max 0≤τ ≤t {f (τ )} is implemented as

The parameters Dr, Hr typically represent mechanical characteristics of the physical machine.

The requirement that the K∞ functions φµr(•), φνr(•), φρr(•) are known, which makes them to play the role of design functions, somehow restricts the remote network disturbance effects tolerated by the resulting controller, in line with the local result we are going to obtain.

Note that the control techniques developed in[8],[START_REF] Marino | Robust adaptive state-feedback tracking control for nonlinear systems[END_REF],[START_REF] Marino | Nonlinear output feedback tracking with almost disturbance decoupling[END_REF],[START_REF] Shen | Robust H ∞ control of uncertain nonlinear systems via state feedback[END_REF] do not apply to the model (3) due to the presence of the uncertain term (θ 5r I 2 qr + θ 6r Per) multiplying the control input u f r .

Note that the dynamics of the relative angular speed ωr can be rewritten in the more compact form

Note that infinite memory is formally required by the controller[START_REF] Marino | Adaptive nonlinear excitation control of synchronous generators with unknown mechanical power[END_REF] to compute the time function max 0≤τ ≤t {f (τ )} with f (•) : R + 0 → R piecewise differentiable (see Section 5 for a finite memory implementation of the controller).

Methods for measuring the power angle δr can be found in[START_REF] Barrera-Cardiel | Microcontroller-based power-angle instrument for a power systems laboratory[END_REF],[START_REF] Chen | A new approach to real time measurement of power angles of generators at different locations for stability control[END_REF],[START_REF] De Mello | Measurement of synchronous machine rotor angle from analysis of zero sequence harmonic components of machine terminal voltage[END_REF].

Acknowledgment

The research was supported by the Italian Ministry of Education, University and Research and through a European Community Marie Curie Fellowship (in the framework of the CTS), c.n. HPMT-CT-2001-00278. The first author is indebted to Prof. R. Marino and P. Tomei for helpful discussions concerning this paper.