Alexander V Terekhov
email: avterekhov@gmail.com

Jean-Baptiste Mouret
email: mouret@isir.umpc.fr

Christophe Grand
email: christophe.grand@isir.upmc.fr

Stochastic Multi-objective Optimization for Aggressive Maneuver Trajectory Planning on Loose Surface

Keywords: autonomous mobile robots, trajectory planning, stochastic optimization

In this work we regard the problem of trajectory planning for aggressive maneuver of a wheeled mobile robot on loose surface. Our approach is inspired by previously reported analysis of professional rally racers actions during sharp turn. Using numerical simulations we obtain a set of solutions representing all possible compromises between speed and accuracy of maneuver. We choose a particular solution and extend it by producing a continuous mapping from desired trajectory turn angles to parameters of the control inputs.

INTRODUCTION

Autonomous field robotics is an attractive and challenging area of research. The problems range from obstacles avoidance, path planning and trajectory tracking under space limitations, to enhancement of the obstacle climbing capabilities on highly uneven terrain. The most famous application, of such researches is the planetary exploration. In most cases, the considered velocities are low and involve only small dynamical perturbations on the system. Thus, most of the developed approaches are based on kinematic modeling. Since recently research projects started to analyze problems induced by the locomotion of ground vehicle at higher velocity (about 10 m/s). One of the famous milestone among past researches was the Great DARPA Challenge that took place in 2005 (see [START_REF] Buehler | The 2005 DARPA Grand Challenge -The Great Robot Race[END_REF]). This kind of challenge requires to find a way to deal with the problems, like locomotion on natural terrain, the control of highly dynamic behavior, the presence of large wheel-ground slippage, problems of perception at high rate frequency, the planning of feasible reference trajectory, etc.

In this paper we focus on the problem of path planning on slipping terrain that allows for aggressive maneuvers at high velocity. The general approach of our work is to control a wheeled mobile robot at high velocity, performing steering maneuvers by taking into account the wheelground slippage. We expect that admitting significant slippage angle during vehicle displacement would allow to noticeably increase the velocity of the maneuver.

In a past work, we have proposed methods to stabilize the control of a mobile robot performing a path tracking task on a slipping terrain (see [START_REF] Lucet | Stabilization algorithm for a high speed car-like robot achieving steering maneuver[END_REF]

, however

The work was funded by DGA, grant REI 2008.34.0018 see also [START_REF] Peters | Mobile robot path tracking of aggressive maneuvers on sloped terrain[END_REF]). In this method the robot was driven by standard kinematic controller, but the final control inputs were corrected by an algorithm that stabilizes the yaw dynamics. The method has been extended to skid-steering robots using a sliding mode controller (see Lucet et al. (2008a)). This method allows to deal with the slippage when the stabilized trajectory is produced by a classical controller. An evident limitation of that work is that it does not allow to take advantage of slippage effect in order to produce aggressive maneuvers like over-steer. In this work, we propose a methodology that deals with the problem of feasible dynamic trajectory planning, i.e., that is consistent with the dynamics of the vehicle.

Recently few attempts were made to approach the problem of highly maneuver steering by getting inspiration from the performance of professional rally drivers (see Velenis et al. (2007a,b)). The driver's actions were recorded during a sharp steering maneuver with 90 degrees turning angle. The authors proposed a reasonable approximation of the driver's actions: the steering, the break and the throttle were shown to be nearly piecewise linear functions of time. Using a computer model the authors have shown that such approximation produces vehicle motion similar to the one exhibited by rally racers.

In order to use the solutions proposed by Velenis et al. (2007a) as a reference trajectory for the vehicle it is necessary to have a set of possible maneuvers including solutions for different trajectory turn angle, road conditions, etc. In fact, one needs a mapping from the set of task conditions to the set of corresponding maneuvers. It is highly desirable that such mapping is sufficiently smooth. Indeed, it may happen, for example, that in course of a maneuver the initial estimation of the road turn angle is discovered to be incorrect and thus another maneuver should have been used. In case of continuous mapping the discrepancy between the performed and the necessary maneuver is not so dramatic and thus, a correction can be made. In addition, the continuous mapping could be used for interpolation of solutions absent in the initial set. The ability of skilled human rally drivers to perform maneuver in different conditions make us think that such a continuous mapping may be possible.

In the current study we aim towards adapting the solution proposed in Velenis et al. (2007a) for a mobile robot with two independent steering wheels and four driving wheels. We want to extend the obtained solution for a range of the angles of the road turn and thus to obtain a continuous mapping from the turn angles to the parameters of the control inputs. For this purpose, we first obtain a mathematical model of the robot, then find a solution for a given road turn angle, e.g 90 degrees, and finally extend the solution continuously for other turn angles. To obtain the solutions we use algorithms of stochastic optimization.

MATHEMATICAL MODEL

For our purpose we need a computationally inexpensive model of a mobile robot, which at the same time captures the main features of its behavior. Currently we are only interested in the robot's movement when all four wheels stay in contact with the ground. For the sake of simplicity we regard only planar motion of the robot (see Fig. 1).

The position of the robot can be defined by the location of its center of mass x, y and the heading angle ϕ between the x axis and the longitude axis of the robot. The motion of the robot satisfies the equations:

ẋ = V m cos ϕ -V l sin ϕ ẏ = V m sin ϕ + V l cos ϕ φ = ω J ω = T M Vm = F m + M ωV l M Vl = F l -M ωV m (1)
where ω is the angular velocity of the robot's trunk, V m and V l are the projections of the linear velocity of the robot's center of mass on the medial (longitude) and lateral axes respectively, F m , F l and T are the total medial and lateral forces and the total torsion torque, defined as the following:

F m = F m1 cos α 1 -F l1 sin α 1 + F m2 cos α 2 - F l2 sin α 2 + F m3 + F m4 F l = F l1 cos α 1 + F m1 sin α 1 + F l2 cos α 2 + F m2 sin α 2 + F l3 + F l4 T = L(-F m1 cos α 1 + F l1 sin α 1 + F m2 cos α 2 - F l2 sin α 2 -F m3 + F m4 )+ d(F l1 cos α 1 + F m1 sin α 1 + F l2 cos α 2 + F m2 sin α 2 -F l3 -F l4 )
F mi , F mi are lateral and medial projections of the tangential forces of wheel-road interaction for each wheel reference frame (see Fig. 1), α 1 , α 2 are the steering angles of the left and right wheels, respectively.

For the forces F mi and F li we use the brush model, which is relatively simple computationally, but at the same time captures main features of the wheel-road interaction.

The details of the brush model can be found in [START_REF] Pacejka | Tire and Vehicle Dynamics[END_REF]. Roughly, the medial and lateral tangential forces are defined as nonlinear function of the lateral and medial projections of the slippage velocity, V sli and V smi , that is the velocity of the contact point of the wheel (in case of no sliding this velocity is zero):

F mi = µF ni f (V smi /V i ) F li = µF ni f (V sli /V i )
where V i is absolute value of the velocity of the center of i-th wheel, F ni is the normal force at the i-th wheel contact point and µ is the coefficient of Coulomb friction. The function f depends on the tangential stiffness of the tire c p .

To determine the normal forces we used the method described in [START_REF] Velenis | Modeling aggressive maneuvers on loose surfaces: The cases of trail-braking and pendulum-turn[END_REF]. The resultant normal forces for i-th wheel is given by the equations:

F ni = M 4dL (dLg + hL i Vl + hd i Vm ) (2)
where h is the height of the center of mass of the robot.

The equation ( 2) describes weight redistribution caused by the acceleration of the center of mass of the robot in case when the pitch and roll angles of the robot are close to zero. However, to compute the medial and lateral acceleration one must provide the total medial and lateral tangential forces, which, according to the brush model, depend on the normal forces themselves. In order to solve this problem we made an assumption that the weight redistribution (2) does not happen instantly but with a characteristic time τ , which roughly correspond to the characteristic time of the suspension system of the robot. We appended the dynamic equations of the robot (1) with the following:

τ ȧl = F l /M -a l τ ȧm = F m /M -a m
and used the values a l , a m instead of VL , Vm in the equation (2).

For sake of simplicity, we ignored the dynamics of the wheels and assumed that the velocity controller in each wheel tracks the desired wheel rotation velocity perfectly. We make the same assumption regarding the functioning of the steering system. Thus, 2 steering angles and 4 wheel rotation velocities represent the control inputs to the model. In order to reduce the dimension of the control inputs we regarded a steering angle α of the front wheels and the linear velocities of the front and rear wheels, V F , V R , correspondingly. The steering angles and velocities of the wheels were computed using Ackermann rule:

Rω 1 = V F 1 - d 2L sin α 2 + sin 2 α 1/2 Rω 3 = V R 1 - d 2L sin α α 1 = arctan sin α 1 -d 2L sin α (3)
where ω 1 , ω 3 are the left front and left rear wheel angular velocities; α 1 is the steering angle of the left front wheel; R is the wheel radius. The formulas for right wheels and steering angle are the same as (3), with the exception that every "plus" is substituted with "minus".

The parameters of the model used in the simulations are presented in Tab. 1.

90 DEGREES TURN

We first regard the case of 90 degrees turn. We assume that the desired path is composed of three parts: two linear parts orthogonal to each other and one circular part connecting them (see Fig. 4). The robot is located 30 m from the point of turn, directed towards that point and has initial velocity of 10 m/s. The simulation is run for 10 seconds, during which the robot is controlled in a feedforward manner. The goal is to perform a maneuver, minimizing the deviation from the desired path and maximizing the average movement velocity. These two objectives are clearly conflicting. To deal with the conflicting objectives, one usually fixes extremal value of one of them, for example, the maximum deviation, at a reasonable value and then solves the optimal problem with single objective and the additional inequality constraint. A significant drawback of this approach is that the choice of the extrema value is to great extent arbitrary, that provokes questions like how good the choice is and how the result of optimization depends on it. Another way to approach this problem is to use the framework of the multiobjective optimization.

Recent research in the stochastic optimization proposed numerous algorithms to simultaneously optimize several objectives (see [START_REF] Deb | Multi-objectives optimization using evolutionnary algorithms[END_REF]); most of them rely on the concept of Pareto dominance, defined as follows:

A solution p * is said to dominate another solution p, if both conditions 1 and 2 are true:

(1) the solution p * is not worse than p with respect to all objectives;

(2) the solution p * is strictly better than p with respect to at least one objective.

The non-dominated set of the entire feasible search space is the globally Pareto-optimal set (Pareto front). It represents the set of optimal trade-offs, that is solutions that cannot be improved with respect to one objective without decreasing their score with respect to another one.

Pareto-based multi-objective optimization algorithms aim at finding the best approximation of the Pareto front, both in terms of distance to the Pareto front and of uniformity of its sampling. The whole set of optimal trade-offs is found using only one execution of such algorithm and the choice of the final solution is left to the researcher/engineer. Compared to classic approaches to multi-objective optimization, which most of the time aggregate objectives (e.g. with a weighted sum) then rely on a single-objective optimization algorithm, multi-objective optimization algorithm don't require to tune the relative importance of each objective. Moreover, depending of the shape of the Pareto front, some optimal trade-offs cannot be found using the most classic aggregation methods (see [START_REF] Deb | Multi-objectives optimization using evolutionnary algorithms[END_REF]).

Typical algorithms of multiobjective optimization first generate a set of N random points, called a population.

Then they enter a loop of four steps until a convergence criteria is met (in this work, a fixed number of iterations is performed):

(1) Sort population with respect to dominance such that non-dominated candidate solutions are ranked 1, those which are only dominated by non-dominated ones ranked 2, etc. Candidate solutions that are attributed the same rank are then sorted with regard to a diversity measure (in objective space) to favor solutions in the less crowded parts. (2) Keep only the best N solutions (during the first iteration, this step is useless). (3) Use the sorted population to generate new candidate solutions by perturbing the kept ones (e.g. by adding a Gaussian noise). (4) Merge the newly generated candidate solutions and the previous population; this gives the new population.

In this work, we employed NSGA-2 (see [START_REF] Deb | A Fast Elitist Non-Dominated Sorting Genetic Algorithm for Multi-Objective Optimization: NSGA-II[END_REF]), one of the most efficient stochastic multi-objective optimization algorithm (see [START_REF] Zitzler | Comparison of Multiobjective Evolutionary Algorithms: Empirical Results[END_REF]); it follows the previously introduced steps. The computations were performed using the Sferes v2 framework described in [START_REF] Mouret | Sferes v2: Evolvin' in the multi-core world[END_REF].

Solving the mentioned optimization problem in the space of all possible control functions is barely possible. So, we used piecewise linear functions of time, similar to those proposed in Velenis et al. (2007a). We parametrized the time change of the steering angle α, the velocity of the front wheels V F and the velocity of the rear wheels V R with the duration of each linear part and the amplitude of its change. Schematic representation for the steering angle is given in Fig. 2. The velocities were allowed to change withing the range from 1 to 10 m/s, the steering angle varied from -60 to 60 degrees. The time intervals were forced to be not longer then 5 s and not shorter then 0.4 s.

In addition, they were checked to not exceed 10 s interval The circles denote the selected solutions. on the whole. For the velocities V F and V B the initial and final values were set to 10 m/s, for the steering angle they were set to 0 degrees. We had 6 parameters for each control input (see Fig. 2) and 18 parameters on the whole.

The estimated Pareto front is provided in Fig. 3. It can be clearly seen that the objectives of speed and accuracy of the maneuver are indeed conflicting and thus each point on the Pareto front represents a compromise between them. The tail of the front (the right side, e.g. the points with high values of maximum deviation) corresponds to the solutions, in which almost no decrease of the average speed can be observed while the steering is performed exclusively by the control of the steering angle of the front wheels. Because of the high deviation from the reference trajectory these solutions are unacceptable for the robot control. In turn, the solutions with small value of maximum deviation are associated with significant deceleration and, as it can be seen from Fig. 3 (the plots, corresponding to 90 degrees turn), require significant change of all three control inputs. In general, we are interested in the most accurate performance possible. However, after some instant (left from the circle in Fig. 3) the average velocity of the maneuver drops dramatically without significant gain in accuracy. Thus we have selected that solution as a reasonable compromise at the point when further increase of performance would require significant drop of velocity. The performance of the selected solution is illustrated in Fig. 4. It must be noted that the produced solution is far non-trivial as it is associated with large slippage angle (the angle between the longitude axis of the robot and the instant linear velocity). It must be noted that the obtained trajectory is qualitatively similar to the one, obtained for humancontrolled mobile robot, reported in [START_REF] Huang | An experimental platform for motion estimation and maneuver characterization in high speed off-road driving[END_REF]. More detailed discussion of the obtained solution is given below.

CONTINUATION

At the second stage we want to obtain a continuous mapping from the desired turn angles to the parameters of the control inputs. We would like to note that the latter can be hardly achieved simply by applying to each turn angle the procedure, described in the previous paragraph. Indeed, since the described problem is expected to have numerous local extrema, such approach, applied for two close values of the turn angle, would lead to rather distant solutions. In order to avoid such behavior the solutions must be forced to stay close to each other, for example, by constraining the distance between the solutions in the space of parameters. Thus, here again we have conflicting objectives: (i) to minimize the distance between the close solutions and at the same time to maximize (ii) average velocity of the maneuver and (iii) its accuracy. This problem can be regarded in the framework of multiobjective optimization. In this case the Pareto front should be a surface in the three-dimensional space of the objective functions. However, solving the problem with 3 objective functions requires an order more computation then with 2. Thus, we decided not to optimize the average velocity of maneuver assuming that minimization of the distance between the solutions would in addition force the average velocities to be close to each other. From the Fig. 7A one can see that this assumption is reasonable enough.

To find the mapping between the turning angles and parameters of the control inputs, we start with solution selected for 90 degrees turn (Fig. 3). Then we compute the Pareto front for 95 degrees turn. The difference between the parameters is defined as following:

ρ(p 1 , p 2 ) = 18 i=1 (p 1i -p 2i ) 2
where vectors of parameters p 1 and p 2 were normalized in such a way that p ji ∈ [0, 1].

An example of the Pareto front for 95 degrees is provided in Fig. 5. It can be seen that the Pareto front has a flat tail (to the right from the circle) associated with nearly linear relationship between the difference from the previous solution and the accuracy of maneuver. The tail proceed to the right until it reaches the solution selected for the 90 degrees turn. The tail is almost flat, meaning that until some point significant gain in accuracy is achieved without deviating a lot from the previous solution. In turn, after that point dramatic growth of difference happens without bringing considerable increase in accuracy. Thus, it makes sense to select the corresponding solution as the best compromise between the two objectives.

The selected solution was used to compute the Pareto front for the 100 degrees turn and so on until 135 degrees. Similar procedure was performed in the direction of decrease of the turning angle until 50 degrees. Examples of Pareto fronts obtained this way for different turn angles are given in Fig. 5. It can be clearly seen that the described relationship between the objective functions holds for another turn angles as well.

Fig. 6 illustrates control inputs selected for different angles of turn. One can see that all three control functions are qualitatively the same for different angles of turn. In general, the maneuver starts with braking with both front and rear wheels, so that the linear velocity of the robot drops down (Fig. 7A). The rear wheels velocity reaches its minimum within about 0.5 seconds. After that a 2 seconds "plato" starts. During this time the rear wheels velocity increases just slightly. It must be noted that profile of the rear wheels velocity is quiet the same for all regarded angles. The front wheels velocity has nearly triangular form: first it decreases until some minimum and then increases with nearly the same magnitude of acceleration. The velocity profile is essentially the same for different turn angles, however it is shifted forward in time for sharper turns.

The steering angle has more complicated structure, which is significantly effected by the angle of turn. For all turn angles the steering starts with rather steep turn of the wheels, followed by the part with relatively slow change of the steering angle, after which it returns to zero rather steeply. For smaller angles of turn the initial steep part is shorter and has smaller amplitude then for larger. It ends approximately at the moment when the front wheels velocity reaches its minimum. The part of slow change is 1-1.5 seconds long and is associated with increase of the steering angle for the turn angles below 90 degrees and decrease otherwise. The start and the end of the whole steering process nearly coincide with those of the plato in the rear wheels velocity. On general, the controller commands appears to depend on the turn angle in rather continuous manner.

The examples of basic kinematic characteristics of the robot motion during maneuver are given in Fig. 7. One can see that the linear velocity of the robot has stereotypical profile, composed of constant deceleration, transients and acceleration. The first and the last part of the profile are nearly the same for all steering angles, while the transients differ in such manner that higher angle of turn is associated with lower average velocity. One can see that the decrease of front wheels angular velocity leads to almost unnoticeable change in the linear velocity of the robot. It happens because the tangential forces of the front and rear wheels nearly compensate each other. For easier understanding one might imagine a steady vehicles, in which the front and the rear wheels rotated in opposite directions with the same velocity.

The angular velocity is significantly more variable. For the smaller angles it seems that the robot slightly overturns and compensate it by negative value of the angular velocity at the end of maneuver. This compensation arises exclusively because of the properties of the wheel-road interaction, since no compensatory steering can be observed on Fig. 6. The slippage angle illustrates significant oversteering in the beginning of maneuver and slight understeering at the end.

CONCLUSION

In this paper we developed the approach to aggressive maneuver planning proposed in Velenis et al. (2007a,b) and inspired by the behavior of rally drivers. We slightly modified it in order to adapt to mobile robot application and gave an estimation of the Pareto front of the speedaccuracy relationship of the maneuver. The Pareto front estimate allowed us to propose a criteria for choosing the "best compromise" between those two objectives. After that we proposed a way to extend the solution for a range of different angles of the road turn. We show that the control inputs have rather stereotypical profiles, which depend continuously on the turn angle. The further development of this work can be performed in several direction. First of all, it seems reasonable to extend the current approach for different road conditions, i.e. different values of Coulomb friction coefficient. Next, it can be combined with existing algorithms of trajectory stabilization, like those proposed in [START_REF] Lucet | Stabilization algorithm for a high speed car-like robot achieving steering maneuver[END_REF]. And, finally, we would be interested in introducing feedback elements into the obtained control As the final goal of the current research we consider implementing the developed control algorithms on a mobile robot that is being currently designed at our institute.

Fig. 1 .

 1 Fig. 1. Schematic representation of the robot.

Fig. 3 .

 3 Fig. 2. Schematic illustration of parametrization used for the steering angle control input.

Fig. 4 .

 4 Fig. 4. The desired and actual trajectory for the selected solution. Black circles denote robot position every 250 ms; the straight lines show the direction of the robot.

Fig. 5 .

 5 Fig.5. Examples of a Pareto front for different angles. The circles denote the selected solutions. on the whole. For the velocities V F and V B the initial and final values were set to 10 m/s, for the steering angle they were set to 0 degrees. We had 6 parameters for each control input (see Fig.2) and 18 parameters on the whole.

Fig. 6 .

 6 Fig. 6. Control parameters in solutions selected for different angles of the turn: A velocity of the front left wheel, B velocity of the rear left wheel, C steering angle of the front left wheel. inputs, which are purely time-dependent at the moment.As the final goal of the current research we consider implementing the developed control algorithms on a mobile robot that is being currently designed at our institute.

Fig. 7 .

 7 Fig. 7. Characteristics of robot maneuver in solutions selected for different angles of the turn: A linear velocity, B angular velocity, C slippage angle.

Table 1 .

 1 Parameters of the model.

	M 40 kg
	J	3 kg•m 2
	L 0.50 m
	d	0.25 m
	h	0.10 m
	µ	0.6
	cp	10 5 N/m 2
	τ	0.05 s


Paul Tonelli

Jean-Baptiste Mouret

Stéphane Doncieux

Influence of Promoter Length on Network Convergence in GRN-based Evolutionary algorithms

Genetic Regulation Networks (GRNs) are a model of the mechanisms by which a cell regulates the expression of its different genes depending on its state and the surrounding environment. These mechanisms are thought to greatly improve the capacity of the evolutionary process through the regulation loop they create. Some Evolutionary Algorithms have been designed to offer improved performance by taking advantage of the GRN mechanisms. A recent hypothesis suggests a correlation between the length of promoters for a gene and the complexity of its activation behavior in a given genome. This hypothesis is used to identify the links in in-vivo GRNs in a recent paper and is also interesting for evolutionary algorithms. In this work, we first confirm the correlation between the length of a promoter (binding site) and the complexity of the interactions involved on a simplified model. We then show that an operator modifying the length of the promoter during evolution is useful to converge on complex specific network topologies. We used the Analog Genetic Encoding (AGE) model in order to test our hypothesis.

Introduction

Evolutionary Algorithms (EA) are nowadays able to offer improved solutions for many problems and sometimes outperform engineering methods. Though, we are unable to obtain solutions as complex as what in-vivo evolution has produced. Some of the evolutionary mechanisms behind biological evolution are still unexplained. It is believed EA could benefit from understanding these mechanisms [START_REF] Banzhaf | GuidelinesFrom artificial evolution to computa-tional evolution: a research agenda[END_REF]. There are many of these mechanisms, as for example splicing, through which a single gene can encode multiple proteins [START_REF] Tischer | The human gene for vascular endothelial growth factor. Multiple protein forms are encoded through alternative exon splicing[END_REF][START_REF] Modrek | A genomic view of alternative splicing[END_REF]. Genetic Regulation Networks (GRNs) are a model of other of theses mechanisms by which a cell regulates the expression of its different genes depending on its state [START_REF] Lee | Others: Transcriptional regulatory networks in Saccharomyces cerevisiae[END_REF]. A strong hypothesis is that the complexity of the interactions obtained thanks to the GRNs is, at least, partly responsible for the diversity increase of many living organisms [START_REF] Levine | Transcription regulation and animal diversity[END_REF], as well as improvements in the genome evolvability and robustness [START_REF] Wagner | Robustness and evolvability in living systems[END_REF]. Though, the complexity of these mechanisms is still not fully understood [START_REF] Wittkopp | Evolution of cis-regulatory sequence and function in Diptera[END_REF].

We believe that understanding which in-vivo characteristics improve the performance of biological evolution is a key to designing efficient EA. In order to do so, it is necessary to identify which elements of these mechanisms are necessary and must be reproduced to improve performance. In this work, we investigated if one of these mechanisms is relevant: a mutation operator adding or removing one base in string based evolutionary algorithms. This operator is usually considered minor, as in-vivo RNA is read three bases at a time (codon) during protein synthesis, and adding or deleting a base in DNA shifts the sequence, creating a protein totally different from the original. Though, in the case of non-coding DNA, this operator may be more useful. In the following work, we asked ourselves two questions. First, does the add/remove operator provide a sufficient mechanism to obtain a correlation between the length of cis-regulatory sequences as stated in [START_REF] Kristiansson | Evolutionary forces act on promoter length: identification of enriched cis-regulatory elements[END_REF]. And more important, does this operator improve the performance of the algorithm in any way ?

2 Related Works

Gene Regulation Networks

GRNs rely on multiple mechanisms. One of them is the possibility for a protein, called a transcription factor to bind itself to a sequence of DNA located before or after a given gene. When a relevant protein binds itself to the site, the expression of the gene will either be enhanced (enhancer) or blocked (inhibitor). Most genes present such cis-regulatory sequences, which contain multiple binding sites for various proteins. These mechanisms can first be seen as a mean to create "programs", enhancing the cell capacity to order the synthesis of proteins or reactions to adapt to specific stimuli [START_REF] Alon | An introduction to systems biology: design principles of biological circuits[END_REF].

GRNs are also believed to speed up evolution. A single mutation in a cisregulatory region can have varying impact on an organism without modifying the gene itself, for example by removing or creating a new interaction between the regulated gene cis-regulatory sequence and another transcription factor. Duplication of a transcription factor gene or binding site also creates new interactions in a genome [START_REF] Teichmann | Gene regulatory network growth by duplication[END_REF]. Therefore interactions provided by the GRNs provide another level impacted by evolution.

Existing Methods

The objective to find ways to harness the properties of these GRNs to improve the performance of evolutionary algorithms is stated in [START_REF] Banzhaf | GuidelinesFrom artificial evolution to computa-tional evolution: a research agenda[END_REF]. More precisely, the goal is to understand which GRN properties improve the evolvability of living organisms. Some properties have already been highlighted in several articles [START_REF] Izquierdo | The Evolution of Evolvability in Gene Transcription Networks[END_REF][START_REF] Tanay | Conservation and evolvability in regulatory networks: The evolution of ribosomal regulation in yeast[END_REF][START_REF] Chen | The evolution of gene regulation by transcription factors and microRNAs[END_REF]. Algorithms have been created to take advantage of them as Artificial Ontogeny [START_REF] Bongard | Evolving modular genetic regulatory networks[END_REF] or lately PBGA [START_REF] Bellas | Using promoters and functional introns in genetic algorithms for neuroevolutionary learning in non-stationary problems[END_REF]. We tried to find a model closely related to the biological mechanisms while avoiding the overhead of more complex biology based models like HeRoN [START_REF] Gonçalves | A Computational Model of Gene Regulatory Networks and its Topological Properties[END_REF].

Research of relevant properties of GRNs

The first step to take advantage of the GRN capabilities is to understand their properties and the implications for evolvability. Examples of these properties are found in [START_REF] Izquierdo | The Evolution of Evolvability in Gene Transcription Networks[END_REF] which tries to understand how varying goals coupled to specific evolution mechanisms can change the evolvability of a genome to speed up convergence on specific problems. It highlights the fact that the specificity of transcription factors to multiple binding sites is, in itself, a way to encode evolutionary information. The mechanisms behind GRNs are quite complex and multiple parameters are still unknown. Here, we restrict our study to algorithms using both a string based encoding while keeping simple enough matching mechanisms. AGE is the only existing example we could find to fulfill these conditions. A work similar to this one was done by C. Mattiussi [START_REF] Mattiussi | Evolutionary synthesis of analog networks[END_REF], who studied the impact of a mutation operator allowing the duplication of sequences in a genome to obtain convergence of the algorithm. Here, we do a similar work with the possibility to incrementally modify the length of the cis-regulatory sequences.

AGE

Analog Genetic Encoding is an indirect encoding mechanism which was designed to use some mechanisms of the GRNs [START_REF] Mattiussi | Evolutionary synthesis of analog networks[END_REF] to generate networks by evolving a string based EA. It has recently shown impressive results for the reverse engineering of existing in-vivo GRNs [START_REF] Marbach | Replaying the evolutionary tape: Biomimetic reverse engineering of gene networks[END_REF]. AGE features complex generation mechanisms, many of which are not relevant to our problem, therefore, we chose to focus only on part of these mechanisms which will be described here. Our goal, as in [START_REF] Izquierdo | The Evolution of Evolvability in Gene Transcription Networks[END_REF] is to assess the capacity of a set of nodes to converge to a specific network topology. In AGE, a network is composed by a set of genes. Each gene is equivalent to a node in a network (see figure 1). Each gene is composed by output and input sequences, each of which is located in the gene and located between a start and an end sequence (which could be compared to start and stop codons). The links in the network are defined by comparing the input and output sequences of different nodes. Similar sequences will be considered as a strong link between two nodes while two totally different sequences will be considered as an absence of link. As expressed in [START_REF] Mattiussi | Evolutionary synthesis of analog networks[END_REF], this is quite similar to the process by which micro RNA can repress the expression of other genes by blocking their DNA or RNA expression [START_REF] Ruvkun | Molecular biology: glimpses of a tiny RNA world[END_REF]. For this model, the strength of the link is computed using a local alignment score [START_REF] Smith | Identification of common molecular subsequences[END_REF].

Experiments

For our experiment, we considered a fixed set of nodes (three to five, depending on the experiment). Each node is composed of two sequences of nucleotides. One input sequence which models the promoter sequence, and one output sequence, loosely modelling the transcription factor / microRNA. Our algorithm mutates theses two sequences by using three possible operators. The first (second) one is the addition (deletion) of a nucleotide at some point in the sequence. Another point is that we do not extend our algorithm to differentiating inhibitory and enhancing links.

The fitness of a network is obtained by comparing the network topology to a reference network. An optimal link strength (1.0) is defined by two sequence containing a common subsequence of three bases. The absence of link is defined as two sequences with a longest common subsequence between two sequences of one base. The intermediary is considered as a link with a strength of 0.5. We made two different experiments related to the length of the sequence of the binding sites. The first one is used to test the correlation between the length of a sequence and the number of links this sequence has in the network. The second experiment tests the effect of switching on and off the add/remove operator on the convergence speed for a "complex" network. All the runs were done using the simplified version of AGE described previously where the only additional mutation operators possible are the exchange of a base in the sequence for another and the addition / deletion operators. In order to improve the performance, the size of the alphabet was set at 7 (this gives the best results for our networks). There is no order relationship between the bases in the alphabet (the replacing mutation operator switches randomly from one base to the other). The selection algorithm used for all the experiments was NSGA 2 [START_REF] Deb | A fast and elitist multiobjective genetic algorithm: NSGA-II[END_REF], a commonly used tournament based multi objective selection algorithm as further experiments required multiple objectives. Each run was repeated at least 10 times. Figure 2 sums up all the parameters used for the experiments.

Convergence of Sequence Length

The first experiment was done in two steps. In the first step, we tried to evolve two simple networks of 3 nodes with homogeneous (2 output links and 2 input links per node) or heterogeneous links (2 outputs for each node but 1 to 3 inputs per node) as shown in figure 3 and we analysed the length of all the sequences of the first individual to reach the optimal fitness in each run. The results of these experiments show a convergence to a size of 7 for each sequence in the homogeneous network. In the experiment trying to converge on a 3 node heterogeneous network (center box of figure fig:networkstype),we obtained a correlation between the length of the sequence and the number of links connected to the node. As the differences were not significant, we made a similar experiment with a network containing 5 nodes and heterogeneous connections (nodes had 3 outputs each and respectively 1, 2, 3, 4 and 5 inputs). All the results (mean length of the sequences and standard deviation) are shown in figure 4.

For the results on the 5 nodes network, we have a significant difference (using Wilcoxon T-test) as the probability of the two sequences being from the same data is less than 1% between all the nodes having different numbers of inputs apart from between sequences 2 and 3 where this probability is 7.5%. These results confirm our first hypothesis, which is that the length of the sequences illustrates the complexity of the interactions the node is involved in. The add/remove mutation operator is therefore sufficient to obtain these results. This is also a confirmation of the results stated in [START_REF] Kristiansson | Evolutionary forces act on promoter length: identification of enriched cis-regulatory elements[END_REF].

How the operator affects performance of the EA

The first experiment showed that the length of the cis-regulatory sequence depends on the complexity of the interactions the node is involved in. Our next step was to test if a network correctly initialized could converge without length modification during the run or if the modification of the length improves performance. In evident cases (simple 3 nodes networks), the add/remove operator doesn't have a significant impact on the results unless initial parameters are initialized to abnormal values. Therefore, we experimented on a network which the algorithm has difficulties to converge to. In order to place ourselves in the worst possible scenario for our hypothesis (necessity of the add remove operator), we took an homogeneous symmetric target network where each node has 4 inputs and 4 outputs and tried to compare the performance. To do so, we first randomly initialized a 5 node network and made it converge to the target network with the addition / deletion operator enabled. It converged in all of the runs and showed a mean length of 8 bases per sequence and a standard deviation around 1. We then ran the same experiments with length modification operators disabled and a fixed initial length of 8 corresponding to the "optimal" length and compared it to a similar network without the length modification enabled.

The results were that the runs without sequence length modification were significantly faster than the runs with the operator enabled. However, this is a case where we specified the optimal length before running the algorithm, which is an unusual situation. Therefore, we made several additional runs of both algorithms by changing only the initial sequence length to compare their performance. The results are given in figure 5 and show that, if the initial length is not optimal, the add/remove operator is a good way to avoid seeing the algorithm get stuck because of insufficient initial complexity. It also helps the convergence rate in non optimal bootstrap cases as can be seen in figure 6. Therefore, we believe that, as the optimal situation of both an homogeneous network and predefined optimal sequence length is unusual, it is usually a good idea to enable the sequence length modification. The alternative being to define the optimal length by another mean before running the algorithm.

Conclusion and Discussion

In this work we have first confirmed that there is a correlation between cisregulatory sequence length and the complexity of interactions the following gene is involved in, but also that the dynamic modification of sequence length is a useful operator (sometimes to allow convergence of the EA, or to facilitate convergence in complex situations). We have also illustrated that a simple evolution mechanism is able to take advantage of these operators, at least in certain cases.

The results shown in the first experiment converge around a length of 7 for both the 3 nodes network while the 5 nodes network converge on a length of 8. Therefore, it could be argued that the optimal length is usually around 8 and that, as this length is sufficient, the add/remove operator can be disabled to improve performance. This is partly true for these sandbox networks (in experiment 2, the best performance is achieved for a bootstrap length of 9). The goal of AGE and other EA is eventually to solve complex problems, with potentially many more nodes and links. In these situations, using a default a length of 8 is raises a risk that the network might not offer enough complexity to converge, as was the case with the 5 bases long runs of the third experiment and therefore be unable to converge.

Fig. 1 .

 1 Fig.1. Left: Two sample nodes in our simplified implementation of AGE; Right: the sequence comparison algorithm for definition of the link between the output sequence of the bottom node and the input sequence of the top node. Only contiguous bases are considered, the longest subsequence for this comparison is therefore of 3 bases even if both also have the F in common.

20 Fig. 2 .

 202 Fig. 2. Summary of the parameters.

Fig. 3 .

 3 Fig. 3. Left: homogeneous network of three vertex (the top circles show the outputs, the bottom ones the inputs); Center: heterogeneous (different number of inputs and same numbrer of outputs) network with 3 nodes; Right: heterogeneous network with 5 nodes.

Fig. 4 .

 4 Fig.[START_REF] Lee | Others: Transcriptional regulatory networks in Saccharomyces cerevisiae[END_REF]. Left: sequence length for an homogeneous network containing 3 nodes; Center: sequence length for an heterogeneous network containing 3 nodes; Right: sequence length for an heterogeneous network containing 5 nodes.

Fig. 5 .Fig. 6 .

 56 Fig. 5. Performance of the EA on a complex network for different bootstrap lengths. bootstrap length 5 6 7 8 9 10 11 12 13 add/remove enabled 1 1 1 1 1 1 1 1 0.8 add/remove disabled 0 1 1 1 1 1 1 1 0.4


Christophe Grand
email: grand@isir.fr

Pascal Martinelli

Jean-Baptiste Mouret

Stéphane Doncieux

J.-B Mouret

Flapping-Wing Mechanism for a bird-Sized UAVS: Design, Modeling and Control

come   L'archive ouverte pluridisciplinaire

Introduction

The capabilities of Mini-UAVs have drastically increased thanks to recent advances in terms of energy storage, effector power and electronic miniaturization, but they still remain far below the maneuverability and energy efficiency exhibited by birds and bats. An European kestrel for instance can stay above a point, fly forward at varying speeds, glide or soar to save energy, while demonstrating maneuverability capacities that far exceed those of the most efficient acrobatic aircraft. We have chosen to study flapping-flight for bird-sized UAVs, i.e., with a wing-span ranging from 0.5 to 1 meter and a weight of 500 g. This article describes current research efforts targeted at designing a flappingwing platform within the ROBUR project ([Doncieux et al., 2006]). Whereas current flapping-flight artifacts show a maneuverability similar to that of fixed-wing engines, our overall objective is to design an aircraft whose capabilities will more closely resemble that of a kestrel, or at least of a pigeon. To this end, the wing kinematics have to be carefully controlled, for instance through the implementation of a neural network-based closed-loop control. As a consequence, the mechanical instantiation of our artificial bird must be able to move the wings on a wide range of periodical and non periodical trajectories, while remaining as energyefficient as possible. This makes the corresponding wing-beat mechanism different from many mechanisms found in the literature [START_REF] Pornsin-Sirirak | Microbat: A Palm-Sized Electrically Powered Ornithopter[END_REF], Vest and Katz, 1999, Raney and Slominski, 2004, Hunt et al., 2005]), because they generate periodical or quasi-sinusoidal movements only. The design of such an improved device is especially challenging since, to properly size its actuators and me-chanical parts, we must know the order of magnitude of torques, angle ranges and dimensions, which strongly depend on the wing-beat kinematics and the morphology of the artificial bird. The optimal kinematics, as well as the required degrees of freedom and the required power to fly at a given speed, being still open scientific questions, we want to address them while designing such mechanism.

This article deals with the design, the modeling and the control of the flappingwing mechanism developed within the ROBUR project. In a first part, we briefly describe the results obtained on the optimal design of bird morphology and kinematics. Based on these results the wing beating mechanism that we have designed is then presented. In a second and a third part, the system modeling and its control are described.

Mechanical design

In a first step, the morphology and kinematics the most adapted to several flying speeds have been sought using an evolutionary algorithm. The results of this optimization give us the first insights of the typical mechanical parameters required for a basic horizontal flight of a 0.5 kg UAV at different speed ([de Margerie et al., 2007]). Optimization lead to a minimum energy consumption of 20-50 W/kg for a 10-12 m/s speed, with wing-beat frequencies between 3 to 5 Hz. Table 1 describes the obtained angle ranges for the different DOFs, for Paretooptimal kinematics. This means that the UAV should at least allow the DOFs to move in these ranges for a basic flight along a straight line.

Speed (m/s) Dihedral Shoulder incidence Wrist incidence 6-8 15-50 0-30 10-50 10-12 25-45 0-15 8-15 16-20 30-65 0-5 1-10

Table 1 Angular ranges (in degrees) for the three DOFs (wing folding is not represented) for Pareto-optimal individuals.

Capitalizing on these results, we designed a wing-beating mechanism that allows a wide range of dihedral and twist variations, with a high energetic and mechanical efficiency. To reach such capabilities, an innovative mechanism was developed as shown in figure 2. In this mechanism, the wings motion is produced by four position-controlled brushless motors (each motor is about 30W power and 100g weight). These four motors are associated two by two, constituting a parallel mecanism that uses two connected rod-crank mechanisms. Thus, the wings can be moved to follow an arbitrary trajectory and the power required to execute the quasi-sinusoidal movements is especially low. Each pair of motors is used to control: (1) the dihedral (DI) motion and (2) the shoulder incidence (SINC) motion.

Kinematical modelling

The wing-beating system is composed of two identical mechanisms. This elementary part is a parallel system actuated by two motors. The first one, located in the front plane, is directely used to control the dihedral motion of the wings and the second one, located on the rear plane, controls the shoulder incidence motion through two conical gears. Each of these elementary systems is made symmetrical by using a drive-belt component (see figure 2). Thus, in this first prototype, the flapping motion is identical for left and right wings (both for dihedral and twist motion).

The elementary parallel mechanism is composed of five rods connected through six revolute joints. The figure 3 shows the kinematic schema of this parallel mechanism. Joint J 1 and J 2 are coupled by the drive-belt, joint J 3 and J 4 are actuated by two motors, and the last two joints (J 5 , J 6 ) remain free. The mobility index of the mechanism is given by the classical Grubler equation :

m = n ∑ i=1 f i -3(n -b)
where b is the number of moving bodies, n the number of joints, f i the number of degrees of freedom of each joint i. In this case, the mobility index is m = 3 (5 bodies, 6 revolute joints), but the overall system mobility is reduced to m = 2 when the symmetry condition from the drive-belt is considered. As these two degrees of freedom are actuated by two motors, the mechanism motion is completely constrained.

J' 3 J' 4 J' 5 J' 6 J 1 J 2 J 3 J 4 J 5 J 6
Wing Wing

ϑ ϑ

Motor 2 Motor 1 Fig. 3 Kinematic schema of the parallel mechanism For analysis purpose, the system can be simplified by the kinematically equivalent system depicted in figure 4(b). Because of the symmetrical motion on joint J 1 and J 2 (angle θ ), the motion of points A 3 and A 4 is constrained to an horizontal axis (A i denotes the center of joint J i ). Thus, we can consider only the 3 rods mechanism and express the kinematic relation between λ and the input angles α 1 and α 2 . The solution is obtained by solving the kinematical closure-form equations :

b cos α 1 + L cos γ -b cos α 2 = λ b sin α 1 + L sin γ -b sin α 2 = 0 (1)
Eliminating γ from these equations gives the following expression:

L 2 = λ + b(cos α 1 -cos α 2 ) 2 + b 2 (sin α 1 -sin α 2 ) 2 (2)
Then, λ can be determined from this second order equation. When the solution existence condition is satisfied (L > √ 2b), we obtain the following unique solution:

λ = L 2 -b 2 (sin α 1 -sin α 2 ) 2 + b(cos α 1 -cos α 2 ) (3)
And considering the schema on the figure 4(a), the wing flapping angle θ is given as a function of λ : The reader should notice that the parameters α 1 and α 2 are not directly the control inputs. But, if we consider u = [u 1 u 2 ] t the input vector composed of the two motor angles (corresponding to the joint angles of J 3 and J 4 ), there can be computed as function of the parameters α 1 , α 2 and the kinematical configuration of the mechanism characterized by the angle θ :

θ = sin -1 L -λ 2a (4) a b L L λ ϑ ϑ A 3 A 4 A 1 A 2 A 5 A 6 α 1 α 2 u 1 u 2 a b γ (a)
u 1 = (α 1 -π 2 ) -θ u 2 = (α 2 -π 2 ) + θ
However for the mathematical description of the kinematical model, it is more efficient to consider the input parameters (α 1 ,α 2 ) instead of the motor angles (u 1 , u 2 ).

The system presents two control inputs α 1 and α 2 for one state parameter θ . So, the mechanism is over-actuated and we need to determine the compatible angles. Thus, let us introduce a new set of input variables α and ϕ, respectively the mean input angle and the half-phase angle :

α = 1 2 (α 2 + α 1 ) ϕ = 1 2 (α 2 -α 1 )
and

α 1 = α -ϕ α 2 = α + ϕ
The equation (3) becomes:

λ = L 2 -4b 2 cos 2 α sin 2 ϕ + 2b sin α sin ϕ (5)
And the relations between the motor angular positions (u 1 and u 2 ) and the new variables become:

u 1 = (α -π 2 ) -(θ + ϕ) u 2 = (α -π 2 ) + (θ + ϕ) and α = 1 2 (u 2 + u 1 ) + π 2 ϕ + θ = 1 2 (u 2 -u 1 )
4 Motion control of the flapping angle

In this part, we demonstrate that the previous choice of variables α and ϕ leads to a simplified control model of the wings-flapping motion. It allows to define quasisinusoidal motion given by the desired flapping frequency and amplitude.

The equation ( 5) can be differentiated with respect to α:

dλ dα = 2b cos α sin ϕ   1 + 2b sin α sin ϕ L 2 4b 2 -cos 2 α sin 2 ϕ   (6)
By considering the schema depicted in figure 3, one can conclude that the wing angle θ is at its extremums when λ is also at its extremums. In order to find the extremum positions α min and α max , the relation dλ dα = 0 must be solved. This leads to the following solutions: cos α = 0 sin ϕ = ± L 2 b Thus, for any values of ϕ there are only two extremum positions at each periode that are solutions of cos α = 0. When introducing this result in the equations ( 5) and (4), we can find the maximum values for λ and consequently for θ :

λ max = L ± 2 b sin ϕ a sin θ max = b sin ϕ (7)
The figure 5 shows some trajectories θ obtained for different phase angles ϕ. The motion is quasi-sinusoidal, its amplitude depends on ϕ and its frequency can be modulated thanks to the velocity term α. So, if the desired quasi-sinusoidal trajectory is specified by its frequency f θ and its amplitude θ max , the velocity α and the phase are determined as follows: In order to extent this simple quasi-sinusoidal control to an arbitrary flapping motion, we need to investigate the differential kinematic model. The flapping velocity is given by θ that can be expressed as a fonction of the input velocities ( α, φ):

α = 2π f θ ϕ = sin -1 a b sin θ max (8)
θ = dθ dλ dλ dα α + dλ dϕ φ = J α α + J ϕ φ (9)
The corresponding jacobian terms are determined from equations ( 4) and ( 6):

                   J α = dθ dλ dλ dα = - b cos α sin ϕ a 1 -L-λ 2a 2   1 + 2b sin α sin ϕ L 2 4b 2 -cos 2 α sin 2 ϕ   J ϕ = dθ dλ dλ dϕ = - b sin α cos ϕ a 1 -L-λ 2a 2   1 - 2b cos α cot α sin ϕ L 2 4b 2 -cos 2 α sin 2 ϕ   (10)
Then, the control principle is to compute the instantanous velocities terms ( α, φ) as a fonction of the desired flapping velocity θ c . This desired velocity becomes the control input used to follow a desired trajectory θ (t) with a classical feedforward controller:

θ c = K 1 (θ -θ m ) + θ (11)
where θ m is the mesured flapping angle, K 1 is a positive gain and θ is the feedforward velocity computed from the trajectory θ (t). Then, the control law becomes:

   φ = K 2 (ϕ d -ϕ m ) α = 1 J α θ c -J ϕ φ (12)
where ϕ d is the phase angle determined from the maximum absolute values of the flapping angle trajectory θ (t) on a given time horizon.

Conclusion

These results constitute the first stones of a long work towards a fully maneuverable flapping-wing UAV. Kinematical and morphological data from an evolutionary optimization process have been used to properly dimension a flapping-mechanism able to move the wing dihedral and incidence to follow arbitrary kinematics. The kinematic model of this innovative parallel mechanism has been detailed in this paper.

Based on this model, a simple control law for quasi-sinusoidal motions has been developed. Its extension to a velocity model based controller that is able to follow various cyclic trajectories has been proposed.

  Figure (1.a) shows the wing panels and their DOFs (dihedral (DI), sweep (SW), shoulder incidence (SINC) and wrist incidence (WINC)). In figure (1.b) possible morphologies corresponding to boundary values of wing area (0.1-0.4 m 2 ) and wing aspect ratio (4.5-10) are illustrated.

Fig. 1

 1 Fig.1The morphology of a simulated UAV.

  Fig. 2 Wing-beating mechanism.

Fig. 4

 4 Fig. 4 Detailed kinematic schema

Fig. 5

 5 Fig. 5 Exemple of θ trajectory obtained with kinematic based control


Jean-Baptiste Mouret
email: jean.baptiste.mouret@lip6.fr

Stéphane Doncieux
email: stephane.doncieux@lip6.fr

Jean-Arcady Meyer
email: jean-arcady.meyer@lip6.fr

Incremental Evolution of Target-Following Neuro-controllers for Flapping-Wing Animats

Using an incremental multi-objective evolutionary algorithm and the ModNet encoding, we generated working neuro-controllers for target-following behavior in a simulated flapping-wing animat. To this end, we evolved tail controllers that were combined with two closed-loop wing-beat controllers previously generated, and able to secure straight flight at constant altitude and speed. The corresponding results demonstrate that a wing-beat strategy that consists in continuously adapting the twist of the external wing panel leads to better manoeuvring capabilities than another strategy that adapts the beating amplitude. Such differences suggest that further improvements in flying control should better rely on some sort of automatic incremental evolution procedure than on any hand-designed decomposition of the problem.

Introduction

Birds continuously demonstrate capabilities which would be of great interest for most Unmanned Aerial Vehicles (UAVs). They are highly agile, able to take off without any runway, to settle on a branch, and to exploit thermals like a sailplane. Consequently, bioinspired flapping-wing platforms could represent useful trade-offs able to benefit from both the manoeuvrability of helicopters and the energy efficiency of standard airplanes.

However, taking inspiration from birds to design a flapping-wing UAV requires a deep understanding of complex aerodynamic principles that nature learned to exploit in about 150 millions years of evolution. In particular, the numerous degrees of freedom of such UAV must be carefully synchronized to produce adequate thrust and lift forces. Additionally, the corresponding rhythmic movements must be continuously adapted to unpredictable changes in the surrounding air mass. To tackle such issues, that are currently not solved by traditional engineering approaches, the ROBUR project of the AnimatLab [START_REF] Doncieux | The ROBUR project: towards an autonomous flapping-wing animat[END_REF] aims at designing the morphology and control of a flapping-wing animat through artificial evolution.

In a previous work [START_REF] Mouret | Evolution of closed-loop neurocontrollers for flapping-wing animats[END_REF], we used a multi-objective evolutionary algorithm to generate wing-beat neuro-controllers able to secure a straight and horizontal flight at constant speed, even in cases where the flying animat was artificially slowed down. Two efficient control strategies emerged, but with no indication about which one should be used in more challenging conditions.

The goal of the research effort described in this article was to extend this work to the control of target-following behavior. To this end, we used an incremental approach that already proved to be efficient at designing wheeled [START_REF] Urzelai | Incremental robot shaping[END_REF][START_REF] Chavas | Incremental evolution of neural controllers for robust obstacle-avoidance in Khepera[END_REF], legged [START_REF] Kodjabachian | Evolution and development of neural networks controlling locomotion, gradient-following, and obstacle-avoidance in artificial insects[END_REF] and flying [START_REF] Barlow | Incremental evolution of autonomous controllers for unmanned aerial vehicles using multi-objective genetic programming[END_REF] robots. We thus capitalized on the previously evolved wing-beat controllers, and let evolution combine their effect with that of newly generated tail controllers that would force the animat to orient itself towards a targeted direction. We also assumed that the corresponding results would help better assessing the relative advantages of the two wing-beat strategies just mentioned.

Fig. 1. The simulated bird is modelled using cones, cylinders and rigid panels. A wings internal panel can be moved along the dihedral and the twist axis, while its external panel can be moved along the twist and the sweep axis.

The corresponding experiments called upon a realistic aerodynamic simulator that computes lift and drag forces whatever the local airflow direction. The underlying aerodynamic model has been validated in a wind-tunnel for a fixed-wing UAV.

The simulated bird was made of cones and cylinders which made up its body, and of rigid panels that composed its wings and tail (figure 1). The total wingspan was 124 cm. Additional relevant details are to be found in [START_REF] Mouret | Evolution of closed-loop neurocontrollers for flapping-wing animats[END_REF].

This article starts with a summary of the previous results we obtained with wingbeat controllers. The next section describes the evolution of tail neuro-controllers for target-following. Sample trajectories and typical neural networks are then exposed, and the corresponding results are discussed.

Wing-beat controllers

Contrary to airplanes which use their wings to sustain themselves, and a propeller to create thrust, birds use their wings to create both an upward and a thrust forces. A wing-beat is made of two distinct phases, the down-stroke and the up-stroke. During the down-stroke, the wings are fully extended and powered downward. The twist is tilted down during this phase, particularly towards the tip. As a consequence, the lift force, created by the pressure difference around the airfoil, is oriented forward and upward. During the up-stroke, the wing is partially folded, to reduce the drag. Additional information on these bird-flight kinematics can be found in [START_REF] Norberg | Vertebrate Flight[END_REF][START_REF] Hedrick | Estimates of circulation and gait change based on a three-dimensional kinematic analysis of flight in cockatiels (nymphicus hollandicus) and ringed turtle doves (streptopelia risoria)[END_REF].

We previously evolved wing-beat neuro-controllers able to exploit such forces and to maintain a flying animat at a constant speed and altitude despite external disturbances. The design of these controllers drew inspiration from the work of biologists on Central Pattern Generators (CPGs) [START_REF] Delcomyn | Neural basis for rhythmic behaviour in animals[END_REF][START_REF] Cohen | The nature of the coupling between segmental oscillators of the lamprey spinal generator for locomotion : A mathematicaal model[END_REF][START_REF] Matsuoka | Mechanisms of frequency and pattern control in the neural rhythm generators[END_REF] and called upon both non-linear oscillators [START_REF] Cohen | The nature of the coupling between segmental oscillators of the lamprey spinal generator for locomotion : A mathematicaal model[END_REF][START_REF] Matsuoka | Mechanisms of frequency and pattern control in the neural rhythm generators[END_REF] and standard McCulloch and Pitt's neurons. These controllers could use a speed sensor as input, and four actuators as outputs: the dihedral and the twist of a wing's internal panel, the twist and the sweep of its external panel. The symmetry between wing-beats was forced.

Evolutionary algorithms are the only means that allow to optimize both the topology and the parameters of this kind of neural-networks. The problem to be solved requiring multiple trade-offs -from the energy consumption minimization to the maximization of accelerations -the multi-objective evolutionary algorithm MOGA [START_REF] Fonseca | Genetic algorithms for multiobjective optimization : formulation, discussion and generalization[END_REF] was used, together with ModNet [START_REF] Doncieux | Evolving neural networks for the control of a lenticular blimp[END_REF], a modular encoding scheme adapted to the task of evolving neural networks. The corresponding fitness function depended upon six objectives [START_REF] Mouret | Evolution of closed-loop neurocontrollers for flapping-wing animats[END_REF], which were evaluated in two stages.

Two different classes of optimal strategies emerged that both relied on the same kinematic principles. According to the first strategy, the twist of a wing's external panel is increased when an acceleration is required, hence increasing the thrust component of the lift force. The wing being folded during the up-stroke, the twist of this panel is not changed. The analysis of the corresponding controllers showed that they implement a simple proportional control of the external twist as a function of the difference between the current speed and the targeted speed.

According to the second wing-beat strategy, all degrees of freedom of a wing exhibit a sinusoidal movement. When an acceleration is required, the amplitude of the oscillations is increased and, consequently, the magnitude of the upward and forward components of the lift force are increased.

Videos of some animats exhibiting these strategies can be downloaded from our website: http://animatlab.lip6.fr/RoburEvolvingEn. To the best of our knowledge, this is the first time that closed-loops controllers are obtained for flapping-wing flight.

Target following

The detailed kinematics used by birds to change the direction of their flight largely remain an open question, the answer to which probably varies across bird species. It has been suggested that the tail might be used for such use, in a manner similar to the use of elevators in an aircraft [START_REF] Thomas | On the tails of birds[END_REF]. However, some birds succeed to fly without their tail, and mostly rely on wing movements for that. For instance, it has been shown that pigeons use down-stroke velocity asymmetries and rapid alternating wing movements to turn [START_REF] Warrick | Kinematic, aerodynamic and anatomical mechanisms in the slow, maneuvering flight of pigeons[END_REF][START_REF] Warrick | Bird maneuvering flight: Blurred bodies, clear heads1[END_REF].

Likewise, although the control of standard airplanes and UAVs is a widely studied topic, described in many textbooks like [START_REF] Mclean | Automatic Flight Control Systems[END_REF] for instance, it remains to be proved that classical control methods can be directly used to control an artificial flapping-wing bird because the dynamics of such an engine are complex to model, and because wing-beats generate a lot of parasite movements, especially when orienting.

Additionally, it turns out that radio-controlled ornithopters built by hobbyists exploit their tail to execute simple manoeuvres, as many radio-controlled airplanes do.

For all these reasons, and because the control of a tail can easily be decoupled from the wing-beat control, in a first approach towards implementing useful flying capacities, we chose to evolve tail controllers that could be combined with the wing-beat controllers previously evolved (section 2). More precisely, the objective of these additional controllers was to orient the artificial bird towards a target point, for instance a GPS way-point or a visual landmark, while keeping its altitude constant. To this end, we used exactly the same methodology, calling upon MOGA and ModNet softwares, that the one evocated above and that led to efficient flapping-wing controllers.

To the best of our knowledge, only two papers previously dealt with the generation of kinematics for a flapping-wing artificial bird able to turn [START_REF] Wu | Realistic modeling of bird flight animations[END_REF][START_REF] Shim | Evolving flying creatures with path following behaviors[END_REF]. But the goal of these research efforts was to create visually convincing movements and not to design closed-loop controllers. Consequently, the optimization of the kinematics parameters characterizing each trajectory were computationally too greedy to be tested as competitive approaches to the one that has been chosen here.

Sensors and actuators

In the absence of any a priori knowledge concerning the sensors that birds use to control their flying manoeuvres, we allowed evolution to incorporate, or not, four sensors in the neural controllers it would generate: an altitude sensor -assessing the difference between the current and targeted altitude -a direction sensor -reorienting the animats relative direction to the target -a roll sensor and a pitch sensor.

Furthermore, in the absence of precise specifications of real sensors to be embedded on a real platform, in this preliminary stage we chose to use ideal sensors that would prevent evolution from exploiting specific characteristics of specific devices. Besides the wing characteristics evocated above, an ideal servo-mechanism was used to move each of the two panels that constituted the animats V-shaped tail. To provide an intuitive control of the bird, these two actuators were mixed together in a way similar to that used in radio-controlled sailplanes. Thus, two virtual actuators are provided, one to control the pitch angle and one to act on the yaw and roll angles. Figure 2 displays some typical reactions of the simulated bird to various tail configurations.

Fitness

A multi-objective evolutionary algorithm was used to seek controllers able to secure a constant-altitude flight while pointing towards a given target. Figure 3 shows an overview of the corresponding control loop.

Let us denote by N the total length of the evaluation. The bird started its flight at 15m.s -1 and the position of its center of gravity g(n) was measured at each time-step n. The altitude objective was written as:

O alt = - 1 N n=N ∑ n=0 |g z (n) -Z| (1)
where g z (n) denotes the altitude of the artificial bird at time-step n and Z the desired altitude. Let first define the vector v(n) which links the current position g(n) to the target position T:

v(n) = T -g(n) (2)
. We then defined the second objective, which rewarded how close the simulated bird was to the desired direction:

O tar = - 1 N n=N ∑ n=0 |atan2(v y (n), v x (n)) -θ(n)| (3)
The atan2(y, x) function calculates the arc tangent of y x except that the signs of both arguments are used the determine the quadrant of the result.

These two objectives had to be maximized, with an optimal value of 0. They might be evaluated for k different targets. In this case, the results of successive evaluations were summed together.

O alt = 1 k i=k ∑ i=0 O alt,i (4)
O tar = 1 k i=k ∑ i=0 O tar,i (5)

Experimental setup

We used the objectives just defined in conjunction with the MOGA algorithm and the ModNet encoding. Population's size was 350 with 60 % of elites. Six different targets were used, towards which the animat was expected to fly. The total evaluation length was 18000 time-steps, simulating 54 seconds of flight (9 seconds for each target). Thanks to the use of a so-called model-module pool [START_REF] Doncieux | Evolving neural networks for the control of a lenticular blimp[END_REF][START_REF] Mouret | Evolution of closed-loop neurocontrollers for flapping-wing animats[END_REF], the neural controllers that were evolved could incorporate and modify across successive generations any number of three different modules:

-a "derivative" module, which computed an approximation of the derivative of its input signal; -an "integral" module, which computed an approximation of the integral of its input signal; -a generic module, made of standard McCulloch and Pitt's neurons, with an evolvable structure.

Using two previously evolved controllers, each exhibiting a different wing-beat strategy (section 2), three evolutionary runs of 500 generations were performed to evolve tail controllers likely to complement them. About 12 hours on 20 Pentium at 2Ghz were required for each run.

Results

First wing-beat controller

Starting with a wing-beat controller that adapted the twist of the external panel to maintain a targeted flying speed, the individuals that constituted the Pareto front of the last generation segregated in two populations: those that obtained good results on the altitude objective, and those that were efficient with respect to the target objective. The behavior of two randomly-selected individuals representative of each of these populations is shown on Figure 4.

While both individuals succeeded to execute light turns, especially during the evaluation period (9 s), they didn't deal successfully with larger changes of directions beyond Trajectories of another individual (Top) with good results on the target objective (Bottom).The targets that were used here were different from those that served to evaluate the individuals. For flights not exceeding the corresponding evaluation period (9 sec), the corresponding behaviors were satisfactory. Beyond this period, some trajectories led to stalling, particularly when the targeted orientation angle was high.

the evaluation period. They used a different approach to handle the latter case. The first individual, which had the best fitness for the altitude objective, stopped orienting when the required change in direction was too large. As a consequence, it did not loose altitude, but at the price of not aligning correctly with the target. The second individual, which had a better fitness on the target objective, often started orienting in the right direction, but ended stalling along a spiral trajectory. However, because such events occurred after the evaluation period, this individual was not much penalized with respect to the altitude objective.

We performed a Multiperturbation Shapley value Analysis (MSA) [START_REF] Keinan | Fair attribution of functional contribution in artificial and biological networks[END_REF] to understand the inner workings of the first individual's tail controller (figure 5). The most useful neurons were those numbered 0, 1, 2, 3, 4 and 5, while the other neurons didn't seem to contribute a lot to the animat's orienting behavior. In particular, neurons 1, 3 and 4 had a large contribution to the first objective, while neurons 0 and 2 mostly contributed to the second objective only. This result indicates that the two objectives were decoupled by the evolutionary process, a conclusion that is confirmed by the close ob- servation of the organization of this controller (Figure 5), which appears as almost split in two separate networks. The first one, on the left side of the figure, controls the flight direction using the target's direction sensor. The second one, on the right side, controls the elevator. These two independent controllers are linked by a connexion with a very low weight, which explains the null Shapley value of neuron 11.

By assigning a null or negative contribution to neurons 5, 12 and 11, the MSA also indicates that the roll sensor is not useful to this controller. Such could be also the case with the pitch sensor, as it will be shown later for an other controller, but the MSA is not conclusive on this point because one cannot decide if the utility of neuron 3 must be attributed to information brought by the altitude sensor, by the pitch sensor, or by both sensors. Be that as it may, we believe that the role of the roll sensor would be much greater if the animat had to fly in an unstable air mass, for instance to secure a constant roll angle despite external perturbations. This hypothesis could be tested by adding perturbations during the evaluation procedure, as we did in [START_REF] Mouret | Evolution of closed-loop neurocontrollers for flapping-wing animats[END_REF].

Second wing-beat controller

Figure 6 shows the behavior of two individuals populating the Pareto front of the last generation, starting with a wing-beat controller that adapted the amplitude of the oscillations to maintain a targeted flying speed. These individuals were randomly-selected among those that respectively obtained good results on the altitude objective, and good results on the target objective.

The first individual exhibits a behavior similar to the one displayed on figure 4 because it stopped orienting when the required changes of direction were large. However, it stalled after having flied for about 200m. Moreover, the error between the animat's direction and the targeted one was larger on figure 6 (a) than on figure 4 (a).

The second individual got the best fitness evaluation according to the angle objective. Surprisingly, although the animat succeeded to perform the largest turns (labeled Both individuals performed correctly during the evaluation period (9s). 9 and 10), it couldnt avoid climbing up. This was due to the fact that the wing-beat controller reached a saturated state, according to which the wings were steadily flapped with a maximum strength, thus producing a maximum thrust. By flying faster, the bird generated more lift and, accordingly, went up. This saturated state was not reached for every target, as demonstrated by the stalling trajectory 7.

The analysis of the corresponding neural controller (figure 7) indicates that it is also split in two sub-networks, one controlling the direction and the other controlling the altitude. For this controller, both the roll and the pitch sensors seem useless. Again, these two sensors would probably be required for flights in an unstable air mass.

Objective space

Figure 8 displays the two objective scores attained by each animat generated during the three evolutionary runs, for the two wing-beat strategies.

It thus turns out that individuals exploiting the first strategy are distributed in most of the objective space, with a higher concentration close to the Pareto front. The shape of this front expresses the necessity of a trade-off, since no individual gets an optimal score on both objectives. The evolutionary process generated a lot of such trade-off solutions situated near the optimum, i.e., at coordinates 0,0. The exploration of the objective space by individuals exploiting the second strategy was quite different. First, the best scores thus attained are substantially lower than in the previous case, especially on the target objective. Second, the explored solutions cover a much limited range of possible scores.

The horizontal line with a target objective value of -0.6 corresponds to individuals that did not turn at all. These individuals quickly disappeared when we used the first wing-beat controller, but they were still present after 500 generations when the second controller was used.

All these results suggest that optimizing a target-following behavior is more difficult with the second controller than with the first.

Discussion

The above results demonstrate that, once neural controllers for straight-line flight have been evolved, it is possible to capitalize on the corresponding networks to evolve additional tail controllers that exhibit minimal target-following capacities. However, as mentioned above, better results would probably be obtained by relaxing the symmetry constraint that was imposed here on wing movements, and by exploiting the manoeuvrability capacities thus offered. Such approach might entail the joint evolution of wingbeat and tail control, a task that seems highly challenging according to current technology. To raise its chances of success, the recourse to some sort of automatic incremental methodology seems mandatory. Indeed, it has been shown here that one cannot rely on fundamental principles or empirical knowledge to decide which, among two available wing-beat controllers, would be better suited to pave the way to additional flying manoeuvrability. Nevertheless, their aptitudes for this endeavour ultimately turned out to be quite different.

Two main reasons motive the use of incremental evolutionary approaches to evolutionary robotics:

-The bootstrap problem. In many real-life situations, an intermediate action is required -e.g. pushing a button to switch-on a light -before receiving any reward -e.g. going to the light to get food. By decomposing the problem into sub-tasks, one may guide the evolutionary process towards satisfying solutions that would otherwise be hard or impossible to discover. -The search space problem. By imposing intermediary stages, one reduces the size of the search space, hence speeding-up the evolutionary process.

Some promising work has been carried out to automatize such a procedure using cooperative co-evolution [START_REF] Haith | Coevolution for problem simplification[END_REF][START_REF] Edwin | Ideal evaluation from coevolution[END_REF]. We plan to assess the applicability of similar approaches to flying behavior in the near future.

Conclusion

Although evolving flying robots seems a greater challenge than evolving crawling, walking, or swimming artefacts, we have shown that a suitable evolutionary algorithm, combined with an efficient coding and a two-stage approach, made it possible to generate close-loop controllers for a target-following flying behavior. However our results suggest that future improvements should better rely on some sort of automatic incremental evolution procedure than on any hand-designed decomposition of the considered problem.

Fig. 2 .

 2 Fig. 2. A v-tail is used to control the artificial bird. (a) When the the panels are in the neutral position, the animat move along a straight line. (b) When the left panel is raised and the right one lowered, the animat turns right. (c) Symetrically, when the left panel is lowered and the right one rised, the animat turns left. (d) When both panels are rised, the animat goes up. (e) Symetricallly, when both panels are lowered, the animat goes down.

Fig. 3 .

 3 Fig. 3. Overview of the control loop. The wing-beat controller has been evolved in a previous work.

Fig. 4 .

 4 Fig. 4. (a) Top: Examples of target-following trajectories for an individual using the first type of wing-beat controller. Bottom: This individual obtained good results on the altitude objective. (b)Trajectories of another individual (Top) with good results on the target objective (Bottom).The targets that were used here were different from those that served to evaluate the individuals. For flights not exceeding the corresponding evaluation period (9 sec), the corresponding behaviors were satisfactory. Beyond this period, some trajectories led to stalling, particularly when the targeted orientation angle was high.

Fig. 5 .

 5 Fig. 5. Left: the neural network that produced the trajectories shown in figure 4 (a). Right: the results of the corresponding MSA. The y-axis gives the numbers identifying the neurons. The x-axis gives the Shapley values. High values indicate important contributions to each of the two considered objectives (altitude and target). Grey levels characterize which neuron is important to the realization of which objective.

Fig. 6 .

 6 Fig. 6. (a) Top: Examples of target-following trajectories for an individual using the second type of wing-beat controller. Bottom: This individual obtained good results on the altitude objective. (b) Trajectories of another individual (Top) with good results on the target objective (Bottom). Both individuals performed correctly during the evaluation period (9s).

Fig. 7 .Fig. 8 .

 78 Fig. 7. Left: The neural network that produced the trajectories shown in figure 6 (a). Right: Its MSA results.

Acknowledgements This work was supported by a UPMC BQR grant.