Sylvain Koos

Jean-Baptiste Mouret

Stéphane Doncieux

Automatic system identification based on coevolution of models and tests

In evolutionary robotics, controllers are often designed in simulation, then transferred onto the real system. Nevertheless, when no accurate model is available, controller transfer from simulation to reality means potential performance loss. It is the reality gap problem. Unmanned aerial vehicles are typical systems where it may arise. Their locomotion dynamics may be hard to model because of a limited knowledge about the underlying physics. Moreover, a batch identification approach is difficult to use due to costly and time consuming experiments.

An automatic identification method is then needed that builds a relevant local model of the system concerning a target issue. This paper deals with such an approach that is based on coevolution of models and tests. It aims at improving both modeling and control of a given system with a limited number of manipulations carried out on it. Experiments conducted with a simulated quadrotor helicopter show promising initial results about test learning and control improvement.

• to obtain a good simulator of the system by local identification of its dynamics; • to find a relevant controller for a given task/problem. Two populations are brought into play: models and tests. Tests are made of a controller, an initial state and a target

I. INTRODUCTION

In evolutionary robotics ( [START_REF] Meyer | Evolutionary Robotics: A Survey of Applications and Problems[END_REF], [START_REF] Nolfi | Evolutionary robotics[END_REF]), evolutionary algorithms are used to design controllers or robots' morphologies: a selective pressure based on the description of the expected behavior leads the search to a solution or a set of solutions. Such an approach relies on the evaluation of the different solutions on a set of performance objectives. But an experimental setup means in particular that the robot can easily be set to a given state. The system can also show risky behaviors when exploring the solution space. To avoid these problems and, moreover, to speed up the evaluation of the solutions, this evaluation step is often performed within a simulator instead of the real system (see [START_REF] Nolfi | How to evolve autonomous robots: Different approaches in evolutionary robotics[END_REF] for discussion).

Once the algorithm has converged, it may be transferred onto the real system. But, because of the discrepancies between the simulator and the real system, this step often means performance loss or even total loss of the previously selected behavior. An accurate simulation model of the system is therefore needed for a good transfer. This transfer problem, also called reality gap [START_REF] Jakobi | Noise and the Reality Gap: The Use of Simulation in Evolutionary Robotics[END_REF], becomes a critical issue for systems that can't be easily accurately modeled like birdsized unmanned aerial vehicles (UAVs) and small helicopters. Because of their sizes and their speeds, these flying robots bring into play little-known dynamics with low/medium Reynolds numbers ( [START_REF] De Margerie | Artificial evolution of the morphology and kinematics in a flapping-wing mini-UAV[END_REF]).

In order to get a good simulation model of the system, one can resort to learning its dynamics with a system identification method (see e.g. [START_REF] Abbeel | Learning vehicular dynamics, with application to modeling helicopters[END_REF]). Such an approach consists Sylvain Koos, Jean-Baptiste Mouret and Stéphane Doncieux are with the Université Pierre et Marie Curis (UPMC) -Paris 06, Institut des Systèmes Intelligents et de Robotique (ISIR), CNRS UMR 7222, F-75005, Paris, France; e-mail: {koos,mouret,doncieux}@isir.fr. This work is part of the ROBUR project.

in performing some experiments on the system to generate data from which a model is learned. Once an accurate model is found, it can be used as a simulator of the system, to design an optimal controller for instance. Nevertheless, a batch identification approach assumes a convenient mean to generate data on the system. For physical systems like UAVs, this assumption requires costly and time consuming experimental setup (pilot, wind tunnel,...). We want then to reduce the number of experiments performed on the system.

To this end, the active learning approach (see [START_REF] Mitrovic | Adaptive optimal control for redundantly actuated arms[END_REF], [START_REF] Robbel | Active Learning in Motor Control[END_REF], [START_REF] Thrun | Exploration in active learning[END_REF]) provides an alternative to batch identification: the data generation becomes part of the identification process in order to perform only meaningful experiments concerning the target issue (trajectory tracking, forward locomotion, ...). One promising way in active learning is coevolution, a method that relies on evolutionary concepts. Bongard et al. introduces in [START_REF] Bongard | Nonlinear System Identification Using Coevolution of Models and Tests[END_REF] an estimation-exploration algorithm that implements an identification method based on coevolution of models and tests. The models are evaluated from previously generated test data and the test (a controller here) that discriminates at most between these simulators is transferred on the system to generate new meaningful identification data. This process is iterated until a good simulator is found. Such a method allows to reduce the number of tests required to explore efficiently the state space.

In the same vein, de Jong et al. addresses in [START_REF] Jong | Ideal Evaluation from Coevolution[END_REF] theorical issues about a coevolution process with multiobjective evaluation called Pareto-Coevolution. It brings into play learners and tests. A learner is evaluated by interaction terms, reflecting each its quality on a given test. In such a context, the existence of a minimal set of tests that leads to optimal learners is proved.

In Bongard's approach, the models are only evaluated from a global test dataset by comparing real data generated from start with predicted ones. But to obtain simulators that are relevant for the control issue, it seems more intuitive to evaluate them directly by some interactions with test controllers, then in a Pareto-Coevolution context.

Our approach is an implementation of the estimationexploration algorithm [START_REF] Bongard | Nonlinear System Identification Using Coevolution of Models and Tests[END_REF] based on Pareto-Coevolution concepts [START_REF] Jong | Ideal Evaluation from Coevolution[END_REF] to improve in a joint process both modeling and control of a robotic system with few experiments carried out on it. The problem includes then two main goals: state. If a test is selected the corresponding trajectory of the real system is generated. Models are learned from trajectories' data. As in Pareto-Coevolution, each model is evaluated by some interactions, also called performance objectives, one by test trajectory. The control task is then directly taken into account for model evaluation. Tests are selected according to some objectives computed with the set of the best models like discrimination power as in [START_REF] Bongard | Nonlinear System Identification Using Coevolution of Models and Tests[END_REF] and, optionally, control quality concerning the defined issue.

An application to the control of a simulated quadrotor helicopter based on [START_REF] Adigbli | Nonlinear attitude and position control of a micro quadrotor using sliding mode and backstepping techniques[END_REF] has then been performed in order to evaluate our approach. A physical simulator of quadrotor dynamics plays the role of the "real" system. We want then to find good models of this simulator according to the selected tests, which have to be relevant on the specified control task.

After a review of some related works in the next section, section III deals with the coevolution process and the different elements which are brought into play: the models and the tests. The implementation of the algorithm is detailed in section IV, in particular the evaluation objectives for each kind of elements. Next, section V details the application to the quadrotor's control. Section VI shows the main results obtained with this application that are discussed in section VII. The paper ends with some conclusions in section VIII.

II. PREVIOUS WORK

System identification is a highly multidisciplinary issue as it's resorted to when modeling any partially-known or unknown system. Whatever the system is, the approach is similar: data are generated and a model is learned. Nevertheless, the data generation can raise some problems with systems like UAVs, for which each experiment means costly experimental setup.

In [START_REF] Abbeel | Learning vehicular dynamics, with application to modeling helicopters[END_REF] and [START_REF] Abbeel | An Application of Reinforcement Learning to Aerobatic Helicopter Flight[END_REF], a Markov decision process model of an helicopter is learned with data collected from a trajectory initially performed by a pilot. The data are then generated before the identification task. A controller is next designed from the obtained model used as a simulator before it is transferred onto the UAV. The approach shows good results, but such a batch approach is not suitable for more open behaviors because it requires to collect all experimental data about the given problem beforehand.

In order to limit the number of experiments, the modeling process has to be led to meaningful parts of the state space of the system concerning the addressed issue. To this end, the data generation can be included in the identification process with an active learning approach [START_REF] Mitrovic | Adaptive optimal control for redundantly actuated arms[END_REF], [START_REF] Robbel | Active Learning in Motor Control[END_REF], [START_REF] Thrun | Exploration in active learning[END_REF]: from a given target (trajectory, behavior,...) the system explores effectively its environment to generate pertinent identification data. The data generation becomes then an integral part of the identification process. Among active learning methods, coevolution is a promising approach based on evolutionary concepts. Coevolution designates the joint evolution of two (or more) populations affected by their interactions.

Bongard et al. introduce in [START_REF] Bongard | Nonlinear System Identification Using Coevolution of Models and Tests[END_REF] an estimation-exploration algorithm that implements a system identification method based on coevolution of models and tests: the models try to model at most the system from previously generated data while the tests try to create new pedagogical data by discriminating at most between the models. For their robotic applications, a test means a controller designed from the learned model for a specified task. The method has been successfully applied to forward locomotion for a simulated quadrupedal robot. Models are sets of parameters and control is provided by neural networks.

The estimation-exploration algorithm has been later implemented on a real four-legged machine [START_REF] Bongard | Resilient Machines Through Continuous Self-Modeling[END_REF]: the robot infers its own structure by performing actions that discriminate at most between the built models. The models are still sets of parameters. From these models the machine can generate forward locomotion. Besides, if a leg is removed, the coevolution process makes models adapt to this new structure and new locomotion's behaviors can then be generated again. In contrary to the previous approach, the tests are only actions and the controller search is performed after the identification task.

Such a method leads to two main advantages. First, as an experiment on the real system is only conducted when a pertinent test is selected, the number of experiments is drastically reduced in comparison with a batch identification approach. The modeling task can secondly be led to meaningful parts of the state space by the evaluation process of the tests, which affects the data generation. We can then obtain a local model good enough to design relevant controllers for a given control task. We can also explicitly improve the models towards this specific task. The evaluation process can bring into play interactions with test controllers to quantify the model's proximity with the physical system concerning the control task. This is the context of Pareto-Coevolution, introduced by de Jong et al. in [START_REF] Jong | Ideal Evaluation from Coevolution[END_REF]. It consists in finding a relevant minimal set of tests to improve the quality of a learner population concerning some underlying objectives of a given problem. As these objectives are often unknown, the learners are evaluated by a set of interactions with tests. Each interaction evaluates the performance of one learner on a given test: each test gives rise to an objective in a multiobjective meaning. In this Pareto-Coevolution context, it is then proved that a minimal set of tests that distinguishes at most between the models' interaction values is an ideal evaluation set and leads the learners to a maximization of the underlying objectives.

There is nevertheless a difference concerning the implementation of the coevolution process between Pareto-Coevolution and estimation-exploration's approach. Coevolution is a priori simultaneous in Pareto-Coevolution context, that is, both populations are evaluated in the same time. In [START_REF] Bongard | Nonlinear System Identification Using Coevolution of Models and Tests[END_REF], this simultaneity aspect is replaced by a more sequential coevolution in two main phases. Our approach, detailed in the next section, is based on this sequential process that we call double evolutionary loop. before detailing the modeling issue. After the description of these elements, our approach is introduced. An outline of the whole method is shown on Figure 1.

III. APPROACH The coevolution process relies on two populations: model and test populations. We first describe what is meant by test

Notations: at time t, the state of the system is called η t and ω t designates a vector of motor commands.

A. Test description

A candidate test Υ is made up of a controller, an initial state, a target state. A controller C is a function which returns motor commands ω t from the current state η t and the target state of the system η d :

• ω t = C(η t , η d )
A controller C can be evaluated on the real system as follows: with a given initial state η 1 and a given target state η d , the command returned by C is applied on the system at each time step and the successive states of the system are recorded. The corresponding trajectory T is then summed up by a dataset:

• T = {(ω t , η t , η t+1 ), t ∈ [0, n]},
where n is the length (i.e. the number of steps) of the control sequence If a candidate test is selected, the controller is transferred onto the real system, the corresponding trajectory is generated and the data are linked to the selected test structure. A whole test is thus generated by each experiment performed on the real system.

For simplicity, we assume in this paper that a candidate test can be summed up by N parameters: proportional-integralderivative (PID) controller's gains, neural network's weights, initial state, target state, etc. Three evolutionary operators can then be defined on this vector representation:

• a mutation operator modifies the value of one parameter; • a crossover operator builds a new candidate test by mixing two parameter vectors of the current population; • a reproduction operator copies a candidate test of the current population to the next one without change.

B. Model description

A model M is a function that predicts ηt+1 , the next state of the system, from its current state η t and the applied commands ω t :

• ηt+1 = M(ω t , η t )
The models are evaluated on the selected test trajectories. The model structure doesn't take any prior information on the system into account. Moreover we assume that the model structure allows to define three operators: mutation, crossover and reproduction. Nevertheless, as the definition of these operators is highly dependent on the model structure, we'll only detail them in the application description (see V-A).

C. A double evolutionary loop

In order to initialize the algorithm, we assume that an initial test Υ 0 is available, that is a controller C 0 , an initial state η 0 1 , a target state η 0 d and the corresponding trajectory T 0 of the real system. We define then a set of selected tests Ω. At the initialization, Ω = {Υ 0 }.

The models are learned from test data and evaluated by comparing their prediction to each selected test trajectory. The test management is a little bit more complicated. The candidate tests are ranked thanks to a set of objectives concerning the models. As in [START_REF] Bongard | Nonlinear System Identification Using Coevolution of Models and Tests[END_REF], candidate tests that discriminate at most between the models are selected. But test evaluation can include other objectives too like control performance in order to lead the search to more efficient controllers and then the identification to more meaningful state space areas concerning the control issue. The influence of such a quality objective will be discussed later.

The coevolution process, or double evolutionary loop, is implemented by the sequence below and iterated as long as no good pair (model, controller) is found:

1) model loop: the best models are looked for from data of the selected tests' set Ω (step 1 in Figure 1), next saved (step 2); 2) test loop: candidate tests are evolved from the previously found best models (step 3); 3) data generation: one among the more relevant candidate tests is transferred onto the real system (step 4) and added with the generated data to Ω (step 5). Both loops are evolutionary algorithms and bring into play the three operators presented above. There is a question left: considering one model and a set of selected tests, how can this model pertinently be evaluated? Likewise, how can we select the best candidate test to transfer onto the system? The next section details the evaluation and selection processes for both populations.

IV. EVALUATION PROCESS

As models and tests are evaluated on a set of objectives, each evolutionary loop boils down to a multiobjective search of optima in a given objective space. In this section, some concepts and tools concerning multiobjective problems are introduced first in order to describe next the evaluation objectives of each population.

A. Multiobjective concepts & NSGA-II ranking method

Let i be an individual evaluated by the objectives obj i j to minimize, j ∈ [1, N ]. The dominance relation between the individuals is a partial order relation defined as follows:

• i 1 dominates i 2 ⇐⇒ ∀j ∈ [1, N ], obj i1 j ≤ obj i2 j ∃j ∈ [1, N ], obj i1 j < obj i2 j
For a given problem, the set of the non-dominated individuals, i.e. the best solutions, is called Pareto front, noted P. On a given population, a non-dominated set S can be defined that contains its non-dominated individuals.

An ideal point z * [START_REF] Deb | Multi-Objective Optimization Using Evolutionary Algorithms[END_REF] can be defined on a given set of objectives. Its coordinates z * j are the minimal values for each objective on the Pareto front P:

• z * j = min i (obj i j ), j ∈ [1, N ],
i ∈ P As the non-dominated set S of a population is an approximation of P, an approximate ideal point ẑ * can be defined too:

• ẑ * j = min i (obj i j ), j ∈ [1, N ], i ∈ S From these definitions, a multiobjective evolutionary algorithm addresses two issues:

• to make S converge to the Pareto front P;

• to maintain a good spread of solutions in S. Among the recent multiobjective evolutionary algorithms, NSGA-II shows good performance on both issues with an reasonable computational cost (see [START_REF] Khare | Performance Scaling of Multiobjective Evolutionary Algorithms[END_REF] for discussion). The algorithm we used in model and test loops is based on the NSGA-II ranking method 1 .

The NSGA-II ranking method [START_REF] Deb | A fast and elitist multiobjective genetic algorithm: NSGA-II[END_REF] is based on the dominance relation between individuals. For a given population, the individuals of the non-dominated set S get 1 as rank value. If these individuals are temporarily removed from the population, we can define a new non-dominated set on it. The corresponding individuals get 2 as rank value and so forth.

Such a ranking method doesn't explicitly maintain diversity among the non-dominated sets. Deb et al. suggest the use of a crowding-distance term which favours individuals in lesser crowded regions of the objective space [START_REF] Deb | A fast and elitist multiobjective genetic algorithm: NSGA-II[END_REF]. Moreover, individuals located on the edges of the non-dominated set are given an infinite crowding distance and thus systematically selected.

Both evolutionary loops are designed as follows: three steps are iterated for a certain number of generations.

1) construction of the offspring of the current population based on elitism and ranking selection with the three operators: mutation, crossover, reproduction (for details, see parts III-A and III-B); 2) evaluation of the new individuals by the corresponding set of objectives; 3) sorting of the new population according to the NSGA-II ranking method and addition of the non-dominated individuals to the non-dominated set S.2 At the end of each loop, the approximate ideal point ẑ * is constructed from the found non-dominated set. For each non-dominated set, we define then a best trade-off for our specific issue by choosing the nearest individual to ẑ * in S according to an Euclidean distance.

From this shared loop structure, the two populations differ by their evaluation objectives, that are described in the next parts.

B. Models' evaluation objectives

As described in part III-B, the models are learned from the selected test trajectories. An objective can then be defined that reflects the prediction performance of a given model against the data of a trajectory (see Figure 2). We call it open-loop objective and its value is the Mean Squared Error (MSE) of the model prediction on the whole trajectory. Nevertheless, the simulation model and the system have to show similarities that concern the control issue too. We introduce therefore a closed-loop objective, computed in a closed-loop process as schematized in Figure 3. With several elements of a given test -a controller C, an initial state η 1 , a target state η d -and a given model M, the two calculations below are iterated for t ∈ [1, n]:

1) ω t = C(η t , η d ); 2) ηt+1 = M(ω t , η t ).
The trajectory T ′ obtained in this way can be compared with the true trajectory T of the test: the MSE between the predicted states ηt+1 of T and η t+1 of T ′ is the value of the closed-loop objective.

These two objectives seem relatively close, but some preliminary results have shown that they are not redundant. As the closed-loop objective is computed from a temporal sequence of system states, it accumulates the model's errors while the open-loop one doesn't. Let N be the number of tests at the beginning of the loop, a model is then evaluated by 2N objectives.

C. Tests' evaluation objectives

The test loop takes place like the model one. It boils down to the choice of a controller, an initial state and a target state to generate a new trajectory of the real system, i.e. new pertinent data. To that end, we want to select a new candidate test that discriminates at most between the individuals of the model non-dominated set S M [START_REF] Bongard | Nonlinear System Identification Using Coevolution of Models and Tests[END_REF].

Let M be a model and Υ a candidate test. The corresponding simulated trajectory can be built according to the previously described closed-loop process. Given the control task, each trajectory can be summed up by the distance δ Υ M between the final state reached by the system η n and the target state η d of the test:

• δ Υ M = η n -η d
The first objective we introduce is the discrimination power of a test [START_REF] Bongard | Nonlinear System Identification Using Coevolution of Models and Tests[END_REF]. A test Υ is even more discriminating as the variance of the δ Υ M values is large on the set P M . The objective value to minimize is then the inverse of this variance.

But, with this objective only, candidate tests whose trajectories are near instabilities of the system are preferred [START_REF] Bongard | Action-selection and crossover strategies for selfmodeling machines[END_REF]. The neighborhoods of the instabilities are indeed parts of the state space where two models can easily be discriminated, even if they are very close. Moreover, it is difficult to model the system's behavior near an instability. To avoid such modeling problems, we assume that the relevant behaviors we look for stay away from instabilities. A second objective is then needed to take into account this assumption.

Consider a model of the system, the best trade-off M * for instance. Consider the candidate test generated from the triplet Υ = (C, η 1 , η d ). A corresponding trajectory T can be built with M * . Consider yet the candidate test generated from Υ ′ = (C ′ , η 1 , η d ), where C ′ is a slightly mutated version of C. The corresponding trajectory T' built with M * is expected to be similar to T. If it's not, there is an instability nearby.

The second objective to minimize is then the variance of the distances δ Υ ′ M * values computed each with a mutant C ′ of the evaluated controller C. This value is even smaller as the corresponding trajectory stays away from the instabilites of the system. A third objective can be added: the control quality of the test. As we want to obtain a pertinent controller for the system, the use of such an objective is indeed attractive. The term δ Υ M defined above can be used to quantify this objective to minimize. Its usefulness will be discussed later.

By choice, the controller population is randomly generated at each beginning of the test loop. The addition of a selected test between the iterations i and i+1 can indeed modify a lot the computation of the discrimination power once this new test is learned. Nevertheless, as the multiobjective algorithm used here doesn't resort to archives, we could also use the population of the previous test loop as initial population.

Notations: the test non-dominated set is called S Υ and the corresponding best trade-off

Υ * = (C * , η * 1 , η * d ).
V. APPLICATION TO THE QUADROTOR'S CONTROL UAVs' study is currently an important robotic field because of their various applications. Moreover, UAVs are systems that can hardly be accurately modeled and experiments on it means costly facilities like a wind tunnel. Batch identification methods are then difficult to use.

We have envisaged an application to quadrotor control from Adigbli et al. [START_REF] Adigbli | Nonlinear attitude and position control of a micro quadrotor using sliding mode and backstepping techniques[END_REF], which designed controllers tested with a physical model of the UAV. To validate our approach in a relatively simple context, the "real" system we want to model and control is this physical model. We have then two main points to check:

1) the best models we obtained are good simulators of the "real" system in a neighborhood of the test trajectories; 2) the controllers of the selected tests are relevant concerning the given control task. We use Linear Genetic Programming as learning method and detail it in the first part. The section deals then with the physical model and the used controller.

A. A learning method: LGP

A system identification approach goes always with a learning method. In evolutionary context, it is often resorted to Genetic Algorithms for parametric identification ( [START_REF] Kristinsson | System identification and control using genetic algorithms[END_REF], [START_REF] Tan | System identification and linearisation using genetic algorithms with simulated annealing[END_REF]) and (often tree-based) Genetic Programming (GP) for symbolic/non-parametric identification ( [START_REF] Bongard | Automated reverse engineering of nonlinear dynamical systems[END_REF], [START_REF] Gray | Nonlinear model structure identification using genetic programming[END_REF]).

Here, we are only interested in non-parametric identification. We have selected a page-based Linear Genetic Programming (LGP) [START_REF] Heywood | Page-based linear genetic programming[END_REF], that is a linear version of GP. Previous comparisons between tree-based GP and LGP on various problems (see [START_REF] Brameier | Linear Genetic Programming[END_REF]) show indeed very promising results with LGP and LGP's implementation is simple. A model is then a sequence of instructions as pictured in Figure 4, i.e. elementary calculations (addition, subtraction, multiplication, division) involving the inputs and constants. In such a program, the same symbol can represent the corresponding input or a temporary vector. Each model contains 5 pages and each page is made up of 10 instructions.

LGP brings into play three operators:

• a mutation operator modifies single instructions of a model; • a page-based crossover swaps single pages between two models [START_REF] Heywood | Dynamic page based crossover in linear genetic programming[END_REF]; • a reproduction operator copies a model of the current population to the next one without change. As in tree-based GP methods, an individual contains introns, instructions that don't affect its prediction whatever the inputs. In LGP, the intron sequences are yet easier to detect and to remove temporarily. It allows then to accelerate significantly the evaluation step with a low additional computational cost.

B. Physical model of the quadrotor -Backstepping controller

A quadrotor helicopter is an UAV that is propelled by four rotors as schematized in Figure 5. In this part, we only address the issue of modeling and controlling the quadrotor's attitude 4 , that is the angles Φ, Θ and Ψ. The state vector η t is thus : As shown in [START_REF] Adigbli | Nonlinear attitude and position control of a micro quadrotor using sliding mode and backstepping techniques[END_REF], let I x (resp. I y or I z ) and τ Φ (resp. τ Θ or τ Ψ ) be respectively the moment of inertia and the torque applied, both about the axis x (resp. y or z), let J R be the torque of inertia of each rotor, the attitude's dynamics can be modeled by the equations below. In the application, this model is the "real" system we want to simulate and control.

• η t = (Φ t , Θ t , Ψ t , Φt , Θt , Ψt )
Φ = τΦ Ix -JR π Ix Θ + Θ Ψ ( Iy-Iz Ix ) Θ = τΘ Iy -JR π Iy Φ + Φ Ψ ( Iz -Ix Iy ) Ψ = τΨ Iz + Φ Θ ( Ix-Iy Iz )
From the applied torques τ Φ , τ Θ , τ Ψ and the current state of the system η t , we can compute its next state by applying the fourth-order Runge-Kutta method to the model and calculating first the new rotational speeds Φt+1 , Θt+1 , Ψt+1 , then the new angles Φ t+1 , Θ t+1 , Ψ t+1 .

The used controller is derived from the backstepping controller described in [START_REF] Adigbli | Nonlinear attitude and position control of a micro quadrotor using sliding mode and backstepping techniques[END_REF]. It consists in one proportionalintegral-derivative controller (PID) controller per angle. Each controller can therefore be summed up by 9 parameters: 3 gains per angle. The torque vector τ ctrl obtained with the PID controllers is next translated in a command vector ω with a conversion chart. The control sequence lasts 2 seconds and commands are sent to the motors each 20 ms.

As described in III-A, a test contains an initial state and a target state too. In order not to increase unnecessarily the number of parameters for this application, we address the issue of the attitude stabilization from the initial state η = 0 to the target state η d = (-30 • , 30 • , -70 • , 0, 0, 0). Thus, only the 9 controller parameters are needed to define one test.

VI. RESULTS

In this section, we want first to show some important points on test learning. In the first part, the tests are evaluated with two objectives only. Some results on the addition of a control quality objective are next presented in a second part. The parameter values used for the performed experiments are shown in Table I (values in brackets are only used in the second part).

A. Global results on the learning of test trajectories

The first points to verify concerns the learning method LGP: does it allow to learn a good model of the physical simulation on a test trajectory? To answer this question, we record the closed-loop objective values obtained with the test Υ 0 and the best trade-off models selected at each iteration of the double loop process in 5 runs 5 . The corresponding means and standard deviations are shown in Table II. The test Υ 0 is significantly learned after the first and the second iterations. Moreover, the models contain in average a higher percentage of intron sequences to 50%. These results validate the choice of LGP as learning method and the size of the models.

Nevertheless the later iterations are accompanied with a performance loss. In fact, only one run performs badly from the 3 rd iteration, but 4 runs do from the 4 th and 5 th ones. The last one doesn't show any significant loss. This perfomance loss is mainly due to the growing objectives number. After some iterations, the model evaluation space becomes indeed too high. The non-dominated set can't converge efficiently to the Pareto front in the model loop with the specified number of iterations. Thus, we often find models that are on the edges of the model Pareto front, then good on some test trajectories and bad on the others.

It is difficult to calculate relevant stastitics on several runs for the other selected tests as they are often different between two runs. The Figure 6 shows information on the selected tests for the more successful run among the 5 ones. Each line corresponds to the closed-loop objective values on a given selected test computed with the successive best tradeoff models. According to the Figure 6, in this run, not only the first test but also the next ones are learned on a permanent basis. A little performance loss appears at the 6 th iteration with the addition of the test Υ 5 . Nevertheless, it shows that the algorithm can find relatively good trade-off models in spite of the growing objective number.

B. Influence of a control quality objective

The second main point that we want to achieve consists in finding relevant controllers for the "real" system. To that end, we compare two versions of the algorithm concerning the test evaluation:

1) with two objectives; 2) with a control quality objective more.

To compare these two versions, we have conducted 25 runs of the test loop for each version on a same model nondominated set, that is after the same run of the model loop.

The control quality value have then been recorded on the "real" trajectory of the test selected at the end of the test loop. Some statistics are shown in Table III. The tests selected in the version 2 show control quality values significantly lower than in the version 1. This is mainly due to the fact that in version 1, selected tests can be relevant about the discrimination power and the instability proximity issues, but totally irrelevant concerning the control task. The addition of a control quality objective is then necessary to lead the search to meaningful control trajectories.

Nevertheless, in both versions, a controller significantly better than C 0 is rarely found. There are two main explanations. On the one hand, due to the evaluation space that grows with the number of iterations in the model loop, the algorithm has been run with a low number of iterations [START_REF] Mitrovic | Adaptive optimal control for redundantly actuated arms[END_REF]. On the other hand, the initial tests have been designed from not so bad controllers that perform relatively well in simulation, then tests with control quality of the same order are often found, but better ones aren't.

VII. DISCUSSION

As tests and models are evaluated by several factors, an intuitive and meaningful way to address model and test selections was to translate our problem into a multiobjective one: each learning/control goal becomes a full-fledged objective. But, an iterative addition of tests raises some questions in such a context.

The first point concerns multiobjective algorithms. Roughly, the actual algorithms are indeed relevant with a number of objectives less than 5, a limit already reached at the 3 th iteration of the evolutionary double loop, that is after the addition of 2 tests only. An important issue is then to reduce the number of objectives. To that end, one can try to merge the open-loop objective on all trajectories because this objective can be defined on any dataset and not only on a trajectory. Nevertheless, it would not solve the whole problem. Another solution consists in limiting the number of tests. We'll come back to this issues later on.

The second point involves the number of generations used for the model loop. The higher the number of tests is, the higher the model evaluation space is and the more it takes time to explore efficiently the model Pareto front from the previous one. On the one hand, too few generations don't allow the convergence of the non-dominated set to the Pareto front. Thus, no good trade-off is found. On the other hand, a too high number of generations leads to overfitting on the first trajectories and the model non-dominated set can't be efficiently extended in the added dimensions. These issues are all the more crucial that the NSGA-II ranking method favours individuals that are on the edges of the nondominated set thanks to an infinite crowding distance term.

To address both problems, one can resort to a number of generations that grows with iterations, but there is no easy way to design such a function. A more pertinent approach introduces a diversity term that encourages models to spread over the Pareto front by adding a new evaluation objective [START_REF] Mouret | Overcoming the boostrap problem in evolutionary robotics using behavioral diversity[END_REF]. Although it adds one more objective, in contrary to the crowding-distance this method allows to favour individuals that would be dominated else, but which are yet pertinent as good trade-offs between the different objectives.

Another approach consists in using a vector of tests with a constant size. At each double loop iteration, some new tests are added to this vector by replacing the oldest ones. It will then fix the evaluation space dimension of the models. Moreover, as old data are progressively "forgotten", such a method would allow the models to adapt to changes of the system dynamics like system damages, collisions, addition/removal of system parts, etc. The process could become less dependent on the initial test Υ 0 too.

The test evaluation can also be improved. On the one hand, staying away from instabilities is convenient, but it penalizes perhaps some relevant behaviors depending on the system. On the other hand, objectives that evaluate if a trajectory is risky or not can be meaningful when transferring onto the real system.

A last unsolved question deals with the first test Υ 0 . In the application, we design indeed not so bad first controllers C 0 from [START_REF] Adigbli | Nonlinear attitude and position control of a micro quadrotor using sliding mode and backstepping techniques[END_REF]. It would be interesting to use worse controllers, or even incomplete initial trajectories in order to study the robustness of the coevolution process in more realistic situations. This issue will be examined in future work.

VIII. CONCLUSIONS AND FUTURE PROSPECTS

This paper addressed the reality gap problem in the case of controller transfer, an important issue in robotics. An identification method of nonlinear systems based on coevolution of models and tests has been implemented. The goal consists finally in automatically modeling and controlling systems (or parts of a system) without the need of an accurate prior model and/or a lot of experiments.

The application to the control of the simulated quadrotor helicopter shows promising results on the modeling and the control issue. Nevertheless, in spite of its intuitive advantages, the iterative addition of tests asks some problems in a multiobjective context. An alternative with a fixed number of tests could be considered. Moreover, previous works on neural networks controllers [START_REF] Mouret | Incremental evolution of target-following neuro-controllers for flapping-wing animats[END_REF] might be used to address the control issue in a more generic way.

Fig. 1 .

 1 Fig.1. Outline of the whole algorithm -A test is made of a controller, an initial state and a target state -Besides, a selected test also contains the corresponding trajectory of the real system.

Fig. 2 .

 2 Fig. 2. Outline of the open-loop process.

Fig. 3 .

 3 Fig. 3. Outline of the closed-loop process.

Fig. 4 .

 4 Fig. 4. Example of a LGP model -(X1, X2) in input, X3 as outputmarked instructions are introns.

Fig. 5 .

 5 Fig. 5. Outline of the quadrotor helicopter.

Fig. 6 .

 6 Fig. 6. Graph of the closed-loop objective values on the selected tests computed with the best trade-off model selected at each iteration of the evolutionary double loop.

TABLE I MAIN

 I PARAMETER VALUES.

	global	model loop		test loop
	process
	number of	number of	population	number of	population	number of
	iterations	generations	size	generations	size	objectives
	7	15000	250	500 (750)	100	2 (3)

TABLE II CLOSED

 II -LOOP OBJECTIVE VALUES OBTAINED ON Υ 0 WITH THE BEST TRADE-OFF MODEL M * SELECTED AT EACH ITERATION.

	it.	0	1	2	3	4	5	6	7
	mean	0.204	0.024	0.029	0.129	0.341	0.357	0.394	0.389
	st.dev.	0.175	0.030	0.030	0.167	0.250	0.327	0.265	0.403

TABLE III STATISTICS

 III ON CONTROL QUALITY WITH EACH VERSION.

		number of	control quality	control quality
		experiments	mean	standard deviation
	version 1	25	0.72	3.11
	version 2	25	0.19	0.29

The selection method is not the same as in[START_REF] Deb | A fast and elitist multiobjective genetic algorithm: NSGA-II[END_REF], but as our approach is a priori not dependent on this implementation choice, it shouldn't have a strong influence on the results

An individual that becomes dominated by another one is, of course, removed from S.

A non-dominated individual remains non-dominated when a new objective is added.

For any details concerning the physical model and the backstepping controller, see[START_REF] Adigbli | Nonlinear attitude and position control of a micro quadrotor using sliding mode and backstepping techniques[END_REF].

The best trade-off model at the iteration 0 is the best model of the initial random population.


Eric Kergosien
email: eric.kergosien@univ-pau.fr

Marie-Noelle Bessagnet
email: marie-noelle.bessagnet@univ-pau.fr

Mauro Gaio
email: mauro.gaio@univ-pau.fr

Semantic cartography: towards helping experts in their indexation task

Keywords: Knowledge Engineering, Thesaurus, Knowledge Representation

Our approach aims at helping experts in their indexation work using the relationship between concepts included in descriptive notices defined by using an external semantic structure (taxonomy, thesaurus, etc).

Introduction

Thanks to the experts' indexation efforts, documents have rich descriptions, notably descriptive notices, described on a well-known semantic structure-base (taxonomy, thesaurus...). We propose to the experts an approach which would make validating their indexation work possible. This approach has two steps: (i) describe information thanks to expert knowledge which is automatically extracted from notices; (ii) give the possibility to navigate within the collection via the identified knowledge representing the indexation work. In a first part (&2), we expose our objectives. Thus, we can develop our approach (&3) to design a specific semantic structure within the PIV1 project. We describe our experimentations for exploration in a corpus indexed by librarians using RAMEAU 2 . We make propositions (&4) to help the validation of the indexation work within the corpus based on the enriched thesaurus.

Objectives

In the aim to propose to experts a validation tool for the use of controlled vocabularies (thesaurus, etc) which they apply to fit their documents analyze, we develop in our approach two preliminary steps. The first step allows extracting and structuring knowledge of our corpus made up of descriptive notices and the controlled vocabulary used to describe these notices. In this step, we connect to research work such as [START_REF] Schatz | Interactive term suggestion for users of digital libraries: Using subject thesauri and co-occurrence lists for information retrieval[END_REF] who attaches to mix terms from a thesaurus and terms from other sources in order to better point information retrieval within management system. One of our objectives is the design of a process to transform a classic controlled vocabulary in a knowledge base [START_REF] Elghoul | Méthodologie de conception d'un SIAD pour la gestion documentaire, aide à l'indexation, aide à la construction du thésaurus, aide à la recherché et aide[END_REF]. The second step proposes a map's representation of this structure, a concepts map [START_REF] Card | Readings in information visualization: using vision to think[END_REF]. In the end, we use techniques of semantic cartography [START_REF] Kohonen | Self-Organizing Maps[END_REF] to tackle in a synthetic way the complete indexation work of a given collection.

The domain-specific thesaurus PIV and a visualization tool

The aim is to represent information thanks to expert knowledge automatically extracted from notices. Indeed, we suppose that we can propose a description of expert knowledge extracting a list of terms (root-words) used in notices by experts to describe documents and the controlled vocabulary used to choose these root-words (RAMEAU for us). According to the AFNOR3 , a thesaurus is a documentary language based on a hierarchical structuring for one or more knowledge domains; notions are represented by terms from one or more natural language and relationships between notions by conventional signs. The first step automatically identifies and extracts « root-words » of the controlled vocabulary (named "subject authorities" in RAMEAU) used to describe the content of the document within XML descriptive notices. These root-words are selected by librarians within RAMEAU and used in notices via predefined tags (figure 1). <DEE>Stations climatiques --Eaux minérales --Barèges (Hautes-Pyrénées) --18e s. </DEE> Fig. 1 -Extract of descriptive notice 1 Each tag DEE contains several root-words separated by element « --». For each term found in a notice, we attach within the XML Topics Map structure [START_REF] Pepper | XML Topic Maps (XTM) 1.0 Specification[END_REF] a link to the document. Authorities represent the conceptual level and document the physical level. This list of root-words is a first step towards the definition of a thesaurus depicting part of librarians' knowledge on the collection; it remains to identify set of relations between these terms in order to develop the sub-thesaurus PIV. By exploiting the thesaurus structure (RAMEAU for the example), we automatically improve the above vocabulary with: « generic », « used for » and « related » (for thesaurus) terms; relations which are linked. The above process gives us an improved organization of knowledge i.e. Domain-specific Thesaurus. Then, we propose a cartographic representation allowing experts to tackle in a synthetic way the complete indexation work of a given collection, realized by different librarians. Up until now this had been very difficult to represent due to the large number of documents and associated descriptive terms.We propose an interface which proposes a global view (i) and a local view (ii). The global view allows librarians to apprehend the expert knowledge structured in the thesaurus in an integral way; the local view, based-on concepts map representation which highlights a subset of the semantic structure [START_REF] Ziti | Hypermedia Exploration with Interactive Dynamic Maps[END_REF]. We explicitly represent the existing relations between the terms of the thesaurus (conceptual level) and the documents of the collection related to these terms (physical level). The combination of both views facilitates navigation and thus re-reading of the indexation work completed on the collection.

Suggestions for assistance in the validation of indexation work

The first feedback from experts of the media library is encouraging, showing the possibility to propose a synthetic display of their indexation work.

Following the automatic thesaurus creation, an automatic control process begin in order to make sure that the terms used in the descriptive notices indexing the collection, are well labelled as root-words in the given controlled vocabulary. This process identifies 4 types of errors due to the usage of the controlled vocabulary: (i) Documents without descriptive notices; (ii) Incorrect information on authorities; (iii) Management of non-selected terms; (iv) Management of homonymy. Concerning the cases linked to name errors and missing words, certain solutions are proposed, by choosing a correct term which is listed. It is possible to index the documents by using terms that are not part of the semantic structure used. If we take the example of a photo of la Place Royale à Pau, the expert can choice the term Place Royale (Pau), a specific term which the committee in charge of updating RAMEAU did not choose to include in the thesaurus. Concerning errors linked to homonymy, we limit correction assistance by suggesting, a list of terms containing the term, which reports the error.

Conclusion

In this article, we have presented our research work on the structuring and adaptation of a controlled vocabulary as an indexing tool. We worked on a collection which includes 750 files and associated descriptive notices. We validated our first experiments with experts of the domain, i.e. librarians. This feedback encourages us to predict new tools permitting to validate the content of descriptive notices. The defined thesaurus representing this indexation work offers a first step for expert users to navigate in the collection through their own representation. Our structure, including a verification phase of the used terms allows in a second step to offer tools for correcting any errors, thereby facilitating the validation of descriptive notices. Our current work focuses on an approach that aims to define an ontology based on a thesaurus. We hope to integrate new knowledge characterizing the territory (adding "localized named entities" and links between concepts). The aim is to offer an interface which permits all types of user, wanting to discover a territory described by documents, to navigate through the collection thanks to domain ontology.

6

PIV : Pyrénnées Itinéraires Virtuels

RAMEAU : http://rameau.bnf.fr/ is a french thesaurus defined by the Bibliothèque Nationale de France (BNF) for or the majority of the French libraries

AFNOR is a French institution for standardization. Documentation : règles d'établissement des thésaurus monolingues. NF Z47-100, 1981.


N Eckert
email: nicolas.eckert@cemagref.fr

Assessing the impact of climate change on snow avalanche activity in France over the last 60 winters using hierarchical Bayesian spatio-temporal change point models

Keywords: Snow Avalanches, Climate Change, Hierarchical Bayesian Modelling, Change Point Models, Avalanche occurrences, Runout Altitudes

Assessing the impact of climate change on snow avalanche activity in France over the last 60 winters using hierarchical Bayesian spatio-temporal change point models

Nicolas Eckert

Avalanche runout altitudes are studied in the whole French territory. A change in runout altitude regime has occurred in France around 1976. Between 1946and 1976, a decrease of 55 m has affected the mean runout altitude, but the probability of a high magnitude event has remained constant. After the change point, the mean runout altitude has regained its initial state, whereas the probability of a high magnitude avalanche has been divided by two. A retreat of avalanche is therefore engaged in France since nearly 30 years. This especially concerns high magnitude events, since the return period associated with an avalanche reaching its minimal altitude on a mean path has increased from 20 to 40 years over the last 30 years.

Avalanche occurrences and runout altitudes are therefore differently influenced by changes in constraining climatic factors, so that the interest of a joint temporal modelling of the two phenomena would be limited. One possible explanation is that dry snow avalanches are progressively replaced by wet snow avalanches because of climate worming, thus keeping constant the number of events, but reducing their magnitude by modifying snow rheology. To confirm this statement, further research is needed to compare and explicitly correlate the obtained annual effects with climatic data such as precipitation and temperature series. This will allow improving our knowledge on climate change in the alpine space and its consequence on avalanche hazard.

INTRODUCTION

Climate change is likely to modify the frequency and magnitude of both ordinary and extreme snow avalanches. Papers dealing with this question generally analyze the evolution of the snow cover rather than the changes in avalanche activity. Snow cover under a given climatic scenario has though been modeled, with the main result being that an increase of the proportion of wet snow avalanches compared to the dry snow ones seems realistic [START_REF] Martin | Impact of climate change on avalanche hazard[END_REF]. Analyses of real avalanche data generally focus on the last few decades. [START_REF] Schneebeli | Destructive snow avalanches and climate change in the Swiss Alps[END_REF] have for instance seen no change in the number of catastrophic avalanches around Davos, Switzerland, during the 20th century. [START_REF] Keylock | The North Atlantic Oscillation and snow avalanching in Iceland[END_REF] has pointed out correlations between avalanche activity in Iceland and the North Atlantic Oscillation.

From the statistical point of view, extracting a climatic signal from an avalanche series is fairly complex. First, avalanche activity is characterized by several variables (occurrence, runout distance, snow volume, etc.) and one or several of them has to be chosen. Second, the different effects that explain the observed data have to be separated in order to extract the common annual component that is related to the large-scale climatic fluctuations. Third, the isolated latent climatic component has to be modeled to extract the structured part of the interannual fluctuations.

Step 2 can be seen as an analysis of variance (ANOVA), with one of the considered factors being the annual component. Depending on the variable studied, a transformation may be necessary so as to relate the observation to a Gaussian and generally Markovian random field. The third step is a classical time series analysis that can put to use a large variety of models, depending of the temporal patterns investigated, for instance memory effects, monotonic trends, or change points for the mean, the variance and/or the extreme values of the series considered. It must be emphasized that steps two and three cannot be separated because the uncertainty related to the estimation of the temporal term must be taken into account for the estimation of the time trend. Steps two and three should therefore be implemented in a hierarchical framework, so as to isolate the temporal component and explicitly model it within the same statistical treatment [START_REF] Clark | Computational Statistics: Hierarchical Bayes and MCMC methods in the environmental[END_REF]. The Bayesian approach is appropriate, providing simulationbased MCMC algorithms for inference [START_REF] Brooks | Markov Chain Monte Carlo Method and its application[END_REF].

In two recent papers, a hierarchical Bayesian framework has been used to study the fluctuations in France of avalanche occurrences [START_REF] Eckert | A spatio temporal modelling framework for assessing the fluctuations of avalanche occurrences under a changing climate: application to 60 years of data in the northern French Alps[END_REF] and runout altitudes [START_REF] Eckert | Assessing the impact of climate change on snow avalanche runout altitudes in France over the last 61 years using a hierarchical change point model[END_REF]. These two variables quantify frequency and magnitude of damageable events and are therefore crucial for hazard assessment. For both of them, the best time series model has been selected among a relatively large class using the DIC criterion. In this paper, the results provided in each case by the best model are summarized. Then, the annual and structured effects are analyzed, looking for possible correlations between avalanche occurrences and runout distances fluctuations over the last 60 years.

FRENCH AVALANCHE DATA

The "Enquête Permanente sur les Avalanches -EPA" is a chronicle describing the avalanche events on about 3,800 determined paths in the French Alps and Pyrenees. Its aim is not to register all the events on all the French avalanche paths, but to be as exhaustive as possible on the paths considered. Avalanche counts are registered, along with several qualitative and quantitative information. The field data are collected by forest rangers and are stored by the Cemagref in a database. In this paper, all runout altitudes acquired since 1946 are considered, whereas avalanche occurrences are only modeled in the northern French Alps, though the most active region at the country scale. Since World War II, observation errors are reasonable in the database, which justifies the choice of a time period 60 obs T = years for both data sets. Note that here a year is a winter, the year 2000 corresponding, for example, to the winter of 2000/2001. Aberrant values were detected using simple statistical tests. This decreases the sample size, but provides more reliable estimates. For runout altitudes and avalanche occurrences, respectively 19,960 and 21,682 events were kept into the analysis. The high number of runout altitudes discarded is explained by multiple error sources (imperfect localization of certain paths in the past, changes in maps available for observation, etc.).

FLUCTATIONS OF AVALANCHE OCCURENCES IN THE NORTHERN FRENCH ALPS

Hierarchical shifting level model

For avalanche occurrences, a spatio-temporal approach is conduced at the scale of the township. The number of avalanches in the township j during the winter t is noted jt a . M is the total number of townships in the region studied and obs T the length of the observation period. A non homogenous Poisson observation model Eckert, Assessing the impact of climate change on snow avalanches using Bayesian change point models is postulated, which is classical for rare discrete events (Figure 1). The parameters Eckert, Assessing the impact of climate change on snow avalanches using Bayesian change point models

Application

Box plots give practical and concise representations of the marginal posterior distribution of each annual component t g (Figure 2, top). The lowest value is obtained for winter 1963/1964, for which only five avalanches were registered. On the contrary, the maximal value is obtained for winter 1994/1995, which was a harsh winter during which the maximal annual number of avalanches was recorded in the EPA database.

Winter 1998/1999, the last catastrophic winter in the European Alps, also shows a very high value for t g , though not the highest. This is explained by the fact that here all events are considered, and not only the most catastrophic ones. The ratio

[ ] 2 2 2 2 g temp v g h r VAR u σ σ σ σ = + + +
quantifies the contribution of the interannual variability to the total variability of avalanche occurrences. Its estimate is 0.17. A little less than 20% of the variability of avalanche occurrences is therefore explained by the interannual fluctuations, which is relatively small, but not insignificant. This low value stems principally from the large spatial variability of avalanche occurrences in the northern French Alps.

FLUCTUATIONS OF AVALANCHE RUNOUT ALTITUDES IN FRANCE

Change point model

For avalanche runout altitudes ijt z , a simpler modelling approach is employed. A normalization by the altitude of the valley floor j min z is used to compare data from one path to another. This defines a Runout Altitude Index (RAI), with ( ) N j the number of avalanches recorded on the path [ ]

1, j M ∈
. M is here the total number of paths under survey (Figure 4). Interannual fluctuations can be summarized by the mean annual index

t t i j t m E RAI   = 
 and the annual probability of high-magnitude avalanches ( )

1 t i j t p P RAI = = .
Since a significant number of avalanches reach the valley floor, the annual distribution of the RAI is modeled by a mixture of two distributions. All slopes except 1 p b are significantly non-zero at the 95% credibility level, and both linear trends are nearly continuous around τ (Figure 3, middle and bottom). The two linear trends explain 53-82% of the variability of mean-and high-magnitude events, depending on the period considered. These high values indicate that the underlying linear trends significantly quantify the variations of mean-and high-magnitude avalanche runout altitudes in France during the studied period before and after the change point.

Obtaining the mean runout altitude fluctuations implies inverting the RAI. This highlights that, in contrast to the mean annual probability of reaching the valley floor, the mean annual runout altitude computed with our model depends on the path considered. With the mean altitude of the valley floor, 1223.3 m asl, the mean behavior at the state level is obtained. Under this assumption, the mean runout altitude has decreased, from a little more than 1400 m asl in 1946 to nearly 1350 m asl in 1977. Since the change point, it has increased again, up to nearly 1400 m asl in 2005. The mean avalanche runout altitude is therefore not different for now than it was 60 years ago in France. However, a clear increase seems to be engaged for nearly 30 years, which is consistent with the general context of climate warming. Moreover, the 1960s and 1970s in France were a short cold period differing from the general context of glacial retreat since the end of the Little Ice Age, which is also coherent with the decrease of mean runout altitudes at that time. These results give good confidence in the ability of the mean RIA to be a meaningful climatic indicator. For the probability of reaching the valley floor, the decrease is strong and continuous since the change point. If the return period corresponding to the valley floor is computes, it appears that it has increased from 20 years in 1980 to nearly 40 years in 2005. The same minimal altitude from a French path is now reached two times less often than 25 years ago. This is important in terms of hazard mitigation, since it suggests that French mountain valleys are now globally less exposed to avalanche hazard than they were earlier. For instance, this indicates that the usual assumption of an underlying stationary process while computing reference scenarios from the available data may lead to over-pessimistic decisions.

Eckert, Assessing the impact of climate change on snow avalanches using Bayesian change point models

DISCUSSION AND CONCLUSION

In the previous sections, a local and small scale approach and a more global approach have been proposed for capturing annual patterns in snow avalanche occurrences and runout altitudes. The two analyzed data sets are not fully coherent (northern French Alps versus the entire French territory), but certain similarities exist in terms of modelling. Totally different structured trends have been obtained. Empirical correlations between estimates are therefore very low, both for annual and structured trends. The only exception is that the mean RAI strongly depends on the annual probability of reaching the valley floor; so that the correlation coefficient ρ equals 0.62 and 0.7 between respectively the annual and structured terms (Figure 5). This shows that avalanche occurrences and runout altitudes are influenced differently by changes in constraining climatic factors. Moreover, this indicates that the interest of a joint temporal modelling of avalanche occurrences and runout altitudes is limited, since there is little information shared by the two phenomena. One possible explanation in accordance with [START_REF] Martin | Impact of climate change on avalanche hazard[END_REF]'s work is that dry snow avalanches are progressively replaced by wet snow avalanches because of climate worming, thus keeping constant the number of events, but reducing their magnitude by modifying snow rheology. To confirm this statement, further research is needed to compare and explicitly correlate the obtained results with climatic data such temperature series. This will allow improving our knowledge on climate change in the alpine space and its consequence on avalanche hazard. Conversely, it will then also be possible to use snow avalanches as proxy indicators that point out further signals of climate change.

Finally, possible unstationarities are for the moment not taken into account in the avalanche engineering practices. Reference scenarios for hazard management are for instance always computed under the assumption of a stationary process. These results show that this should be reconsidered in the near future given the strong modification of runout altitude regime that has occurred during the last 25 years.

Figure 1 .

 1 Figure 1. Model for avalanche occurrences. The circled nodes represent stochastic latent variables and the diamonds represent the overall parameters, while the rectangles indicate observed values.

Figure 2 .

 2 Figure 2. Box plots of the annual components in avalanche occurrences, mean runout altitude and probability of reaching the valley floor in France. No monotonic trend such as a constant increase or decrease over the period studied is visible on Figure 2 (top). This indicates that there is for the moment no irreversible change in avalanche activity in the northern French Alps. On the other hand a combination of different time structures exists. First, memory effects can be suspected, with many sequences of consecutive years showing close values, for example between 1949 and 1954. Second, cyclic variations seem to occur, with a succession of four cycles of around 15 years, with maxima around 1951, 1965, 1980 and 1995. Finally, abrupt variations also exist, with very low values directly following very high values, for example,between 1962 and 1963 or between 1987 and 1988. The interannual trend t trendg provided by the shifting level model consists therefore in a complex combination of pseudo-periodic cycles and abrupt changes(Figure 3, top). The ratio

Figure 3 .

 3 Figure 3. Box plots of the structured annual components separated from random noise in avalanche occurrences, mean runout altitude and probability of reaching the valley floor in France.

Figure 4 .

 4 Figure 4. Change point model for mean and extreme snow avalanche runout altitudes 4.2. Application

Figure 5 .

 5 Figure 5. Scatter plots and empirical correlations between point estimates obtained with the two models

  's are modeled as nonexchangeable latent variables with a model simple, but one flexible enough to capture a monotonic trend and various types of changes in mean and variance. τ is the year of a possible change point separating two periods of different runout regimes. Before and after the change point, both the

								1 RAI is a discrete random number taking the value 1 if the avalanche i ijt
	reaches its possible minimal altitude, with an annual probability t p .	RAI	ijt	2	is a continuous random
	number whose annual distribution is possibly skewed and can therefore be modeled by a Beta distribution
	with an annual parameter pair (	) α β . The annual mean t , t t m is easily obtained from the model's
	parameters. As a consequence, the triplet (	) p m β fully characterizes the RAI annual distribution. To , , t t t
	capture the systematic variations of mean-and high-magnitude avalanche runout distances, the t m 's and the
	t m 's and the t p 's are broken down into random noises and structured linear trends. The random noises,
	with variances (	2	1	,	2	2	2 1 , ,	2	2	)

t p