Cherifa Dad

Stephane Vialle

Mathieu Caujolle

Jean-Philippe Tavella

Michel Ianotto

Scaling of Distributed Multi-Simulations on Multi-Core Clusters

DACCOSIM is a multi-simulation environment for continuous time systems, relying on FMI standard, making easy the design of a multi-simulation graph, and specially developed for multi-core PC clusters, in order to achieve speedup and size up. However, the distribution of the simulation graph remains complex and is still the responsibility of the simulation developer. This paper introduces DACCOSIM parallel and distributed architecture, and our strategies to achieve efficient multi-simulation graph distribution on multi-core clusters. Some performance experiments on two clusters, running up to 81 simulation components (FMU) and using up to 16 multi-core computing nodes, are shown. Performances measured on our faster cluster exhibit a good scalability, but some limitations of current DACCOSIM implementation are discussed.

I. INTRODUCTION

Complex systems are characterized by the interconnection of numerous and heterogeneous cyber components (e.g. controllers) and physical components (e.g. power grids). For EDF (the major French utility company), the smart power grids will extensively rely on new control functions (i) to increase the grid efficiency, reliability, and safety, (ii) to enable better integration of new assets (e.g. distributed generation and alternative energy sources), (iii) to support market dynamics and manage new interactions between established and new energy players. Cyber-Physical Systems (CPS) design involves multiple teams working simultaneously on different aspects of the system. Especially, the development of a smarter electrical grid requires to reuse models and tools often based on separate areas of expertise.

This heterogeneity led EDF to investigate coupling standards such as the Functional Mock-up Interface (FMI) initiated by Daimler AG within the ITEA2 MODELISAR project and now maintained by the Modelica Association 1 . More precisely, EDF chose its FMI-CS (FMI for co-simulation) part because this operation mode allows to export models as active components called FMUs (Functional Mock-up Units), each FMU being a self-contained archive file including a model and a numerical solver. As an additional benefit, the model IP is readily protected.

For EDF, it is vital to develop agile modelling and simulation in order to design and validate distributed operating functions a long time before performing tests on experimental sites. The Distributed Architecture for Controlled CO-SIMulation (DACCOSIM software [START_REF] Galtier | FMI-Based Distributed Multi-Simulation with DAC-COSIM[END_REF]) developed by EDF and 1 https://www.fmi-standard.org/downloads CentraleSupelec is a part of the answer to this problematic. It is aimed at simulating large and complex multi-physics systems on multi-core PC clusters (the standard scalable computing architecture), despite some load unbalance and important communications inherent in our multi-simulations. This paper introduces the parallel and distributed architecture of DACCOSIM, and some scaling experiments on PC clusters (running up to 81 FMUs and using up to 16 multi-core cluster nodes). Finally we list some issues and future works to achieve better speedup and size up.

These researches are carried out by the RISEGrid2 institute, founded by EDF and CentraleSupelec with the ultimate goal to simulate the smart electric grids of the future.

II. CONTEXT AND RELATED WORKS

A FMI based multi-simulation is defined by a FMU graph, where all inputs and outputs need to be correctly initialized during a mainly iterative phase called co-initialization. EDF has achieved the co-initialization phase of DACCOSIM reusing its expertise in the Newton-Raphson algorithm (very popular for power flow calculations) in conjunction with recent works done at UC Berkeley on dependency cycles analysis [START_REF] Broman | Determinate composition of FMUs for co-simulation[END_REF] (which are now possible with the latest version 2.0 of the FMI standard). The global dependencies graph automatically deduced from the FMU graph is valuable to solve algebraic loops with a parallel Newton-Raphson algorithm staged at the initialization mode of the multi-simulation. Then, the FMU graph enters a loop of parallel multi-simulation time steps, concurrently running each FMU at each time step, without the need for additional power flow calculations between two consecutive steps. From this point of view, the current version 2.0 of the FMI standard seems sufficient for EDF as its current use cases do not report non convergence examples due to algebraic loops in the step mode.

So, our problematic is to distribute a FMU graph on a multi-core PC cluster. This looks like a task graph distribution problem, which is an important research field. For example, [START_REF] Daoud | A high performance algorithm for static task scheduling in heterogeneous distributed computing systems[END_REF] proposes a scheduling algorithm on heterogeneous distributed architectures for static tasks modeled by Directed Acyclic Graph (DAG), and [START_REF] Dummler | Mapping algorithms for multiprocessor tasks on multi-core clusters[END_REF] introduces mapping strategies of hierarchically structured multiprocessor tasks on multicore clusters, which is our target architecture. However, our DACCOSIM FMU graph is not hierarchically structured, and is a DAG with very few task dependencies. All FMU tasks can run in parallel when starting a new time step, and when all computations are finished all inter-FMU communications can start and occur in parallel [START_REF] Galtier | FMI-Based Distributed Multi-Simulation with DAC-COSIM[END_REF]. Then, when all communications are achieved, all FMUs can enter a new time step. There are no dependencies between our FMU task computations, and usual solutions of task graph distribution are not really adapted to our problem. The basic version of our FMU graph working is closer a classical Bulk Synchronous Parallel (BSP) model [START_REF] Valiant | A bridging model for parallel computation[END_REF], which is well known in parallel computing, but our FMU tasks are heterogeneous and load unbalanced. Moreover, a computationally big FMU can not be split into smaller ones as it would require to design new mathematical models. So, our problem is not a classical Simple Program Multiple Data (SPMD) scheme following a BSP model. Our FMU graph has heterogeneity and load unbalance typical of generic task graphs. It is also possible to consider our FMU graph like a kind of time stepped Multiple Program Multiple Data application running on a computing cluster and aiming to reach high performances [START_REF] Chang | Evaluating the Performance Limitations of MPMD Communication[END_REF]. Identification of an efficient distribution of our FMU graph thus requires to design a specific algorithmic solution.

Some others multi-simulation environments interconnect some continuous time based simulators, like EPOCHS [START_REF] Hopkinson | Epochs: a platform for agent-based electric power and communication simulation built from commercial off-the-shelf components[END_REF] and INSPIRE [START_REF] Georg | Inspire: Integrated co-simulation of power and ict systems for real-time evaluation[END_REF], or even some FMUs, like C2WT [START_REF]Model-Based Integration Platform for FMI Co-Simulation and Heterogeneous Simulations of Cyber-Physical Systems[END_REF], through a HLA logical event bus [START_REF] Kuhl | Creating Computer Simulation Systems: An Introduction to the High Level Architecture[END_REF]. This bus relies on a Run-Time Infrastructure (RTI), ensuring event routage and right synchronization between simulators. Many RTI implementations are distributed, and can activate concurrently simulators on different computing nodes. However, parallelism of the system depends on the abundance of safe events concentrated in a lookahead time interval [START_REF] Fujimoto | Parallel and Distributed Simulation Systems[END_REF], [START_REF] Wang | Exploiting event-level parallelism for parallel network simulation on multicore systems[END_REF], and our time stepped FMU graph contains more potential parallelism.

III. DACCOSIM ENVIRONMENT

A. DACCOSIM software suite DACCOSIM has been designed to achieve multisimulations of continuous time systems, discretized with time steps, and running solvers using constant or variable time steps (to maintain the requested accuracy with the minimal amount of computations, whatever the dynamic of the system). It consists in two complementary parts: a friendly Graphic User Interface (GUI) and a dedicated computation package. The GUI developed in Java facilitates the complex systems studies by designing the multi-simulation graph (Figure 1), i.e. the FMUs involved and the variables exchanged in-between, defining the resources used by the simulation (local machine or cluster), configuring the simulation case (duration, coinitialization method, time step control strategy. . . ) and implementing the graph into DACCOSIM master tasks managing the simulation.

The dedicated computation package controls all task execution issues relative to the multi-simulation: co-initialization, local or distributed computation steps, fixed or variable time step control strategies, detection of state events generated inside FMUs, inter-FMU communications, distributed and hierarchical decision process. . . The Java version of DACCOSIM relies on JavaFMI3 and is available for both Windows and Linux operating systems, whether 32-bit or 64-bit.

B. Parallel and distributed runtime architecture

The dedicated computation package of DACCOSIM includes a parallel and distributed runtime architecture, designed to take maximal advantage of any cluster of multi-core nodes. A DACCOSIM simulation executes a series of time steps composed of three stages: the time step computations of the FMUs (independent computations), the communication of the FMU outputs to the connected inputs (many small communications), and the simulation control by the hierarchical control master (information gathering, next operations decision and order broadcasting). All these operations include potential parallelism exploited by DACCOSIM runtime.

As inside one computing step all FMUs can achieve their computation concurrently, DACCOSIM architecture distributes FMUs on cluster nodes, encapsulates FMUs of a same node in different threads and implements a hierarchical (and DACCOSIM uses FMUs in co-simulation mode: they embed their solvers and are implemented as dynamic libraries (enhanced with meta-data XML files). A FMU wrapper thread encapsulates each FMU, calls its computation function to achieve a time step progress, and sends each of its outputs to the connected inputs. Three others threads are associated to each FMU, as illustrated on Figure 3: one receipts the output values coming from other computing nodes (currently across ZMQ middleware), one achieves similar receipt but from FMUs located on the same node (not going through the middleware to run faster), and one stores simulation results on disk (asynchronously and per block). Received input values are stored into buffers and distributed to the real FMU inputs before the start of each FMU computation step, so that the inputs remain stable during the computations.

When using variable time steps, the control master gathers and analyses some time step results. It decides if they are valid and if the simulation can enter the next step with the same or a larger time step, or if the simulation has to roll back and continue with a smaller time step (to reach the required accuracy). Two orchestration modes are available in DACCOSIM (see Figure 4). The ordered mode executes the three phases of each time step in order: (1) FMU computations, (2) communication and control with the hierarchical master, and (3) inter-FMU communications if the master has validated the time step results. The overlapped mode overlaps the inter-FMU communications with the control master operations (phases 2 and 3). If the control master requires a rollback, the received input values are forgotten and the FMU states at the beginning of the time step are restored. But when the dynamic of the simulated system is limited (no turbulence), there are few rollbacks and the overlapped mode reduces the simulation time. When the dynamic is high, many rollbacks can appear and a lot of inter-FMU communication phases can be achieved needlessly. Currently the orchestration mode has to be set by the user when running the multi-simulation. In the future it could be auto-tuned dynamically depending on the frequency of the rollbacks. A deployment tool (dacrun) completes the DACCOSIM software suite, and makes easy the deployment of a multisimulation on a set of cluster nodes. It is compliant with the OAR4 environment, allowing to allocate nodes and run interactive or batch jobs on a PC cluster, but can be easily adapted to any similar cluster management environment. However, design of a multi-simulation graph distribution on different cluster nodes is the responsibility of the multi-simulation developer. An automatic distribution of the FMUs is under investigation, according to the lessons learnt from our experiments (see next sections).

IV. PARALLEL AND DISTRIBUTED EXPERIMENTS

This section introduces our scalable benchmarks and PC clusters, and our methodology to identify the best distribution patterns of elementary multi-simulation benchmarks. Then we detail our scaling approach, and our size up experiments on larger benchmarks.

A. Testbed introduction

Our test case is a simplified industrial case provided by EDF R&D, representing heat transfers in a set of n three-floor buildings, with two zones per floor separated by an indoor wall. The entire multi-simulation includes 1 + 10 × n FMUs, interconnected as illustrated in Figures 5 and6. The unique FMU BC is in charge of the thermal boundary conditions (e.g. actual temperatures recorded in a French suburb) which are supposed identical for all the floors in every building. A complementary FMU CS models the crawl space temperature for the lower floor of a building, and each building having a specific FMU CS. The FMUs ZNx (resp. ZSx) represent the Northern (resp. Southern) part of every floor, each being designed with Modelica differential equations modeling physical phenomena such as conduction, convection and solar radiation. The FMUs IWx detail the behavior of the wall between the floor zones depending on different insulating properties. All these FMUs are equation-based only, and modeled by

B. Speedup tracking strategy

Identifying the optimal distribution of a multi-simulation graph on a multi-core PC cluster encounters several locks. First, strong differences appear after measuring the computing load of each FMU. FMUs are heterogenous tasks leading to complex load balancing between computing nodes. Moreover, splitting time consuming FMUs into several smaller ones would require if possible to redesign their mathematical models and would imped to quickly reuse existing simulation components.

Second, installing only one FMU per multi-core node would be a waste, requiring many under-used nodes and increasing the inter-node communication time. But running less than k FMUs in parallel on a k-core node leads to a global computation time greater than the longer FMU one. Figure 7 shows the computation time of a mix of big and medium (half a big) load FMUs on one 4-cores/8-threads Nehalem node. Up to 4 FMUs, the computation time should match the horizontal thin lines, but it appears to follow the thick curves (we also observed this behavior on other multi-core nodes). Most of parallel codes are memory bound, excepted when achieving very regular memory accesses and implemented with tiling (subdivided and blocked in cache memories). When codes are not tiled, parallel computing time on each multi-core node becomes hard to predict, specially with heterogenous tasks. Computing a load balanced distribution of a FMU set on a multi-core cluster is thus not obvious.

Third, a multi-simulation graph achieves many communications between FMUs, as well as between each FMU and the hierarchical master. Our measures have shown there are numerous and small communications, sensitive to the internode network latency. Communication times can be significant compared to computation times (like in our cl3 benchmark), and inter-node communications remain longer than intra-node communications. So, the distribution strategy of our multisimulation graphs has to take into account both computations and communications issues.

To achieve our first investigations we adopted an experimental methodology: (1) measuring each FMU computation time (running our real multi-simulation with only one FMU per node to avoid FMUs disturb each others), (2) measuring computation times of different mix of FMUs on one multicore node, (3) deducing FMU distributions with best load balancing on different numbers of nodes, and (4) tuning these load balanced oriented distribution in order to decrease the number of inter-node communications when communications appear to be important. We applied this methodology to our cl3 benchmark with 1 building. Figure 8 shows the multisimulation times we measured on our Nehalem cluster function of the number of used nodes. Two distribution strategies (s1 focuses on load balancing while s2 takes into account both computations and communications) and two orchestration modes are considered (see section III). Finally, the best distribution we identified and experimented for our cl3 benchmark with 1 building uses 3 nodes and is illustrated on Figure 10 (FMUs with the same color will be located on the same node). Figure 9 shows the little tuning (sacrifice) of load balancing successfully achieved by strategy s2 on 3 nodes to increase the number of intra-node communications (decreasing the number of inter-node ones). For our cl5 benchmark with 1 building, including huge computations and negligible communications, the best distribution we identified on our Nehalem cluster uses 4 nodes and is illustrated on Figure 5. One of our objectives is now to automatize this investigation methodology, in order to quickly analyze new FMU graphs and define some efficient distribution patterns, like the optimal distributions of our cl3 and cl5 benchmarks with 1 building. Then, these patterns will be reused in the distribution of larger multi-simulation graphs.

C. Scaling and size up experiments

All experiments introduced in this section use the overlapped orchestration mode, and avoid to go through the middleware for internal node communications. This configuration appears the most efficient on our benchmarks and clusters. We measured multi-simulation time with result storage (approximately 850M B per run), however we stored these results on the computing nodes local disks, in order to avoid fluctuations due to network traffic on our LAN.

We conducted many experiments to measure size up of our multi-simulation benchmarks: we increased both the problem size (number of buildings simulated) and the number of computing ressources (number of cluster nodes), attempting to keep constant the execution time. When processing larger multi-simulations we replicate the best distribution pattern identified for one building (see figure 5 and6). So we attempt to obtain:

T (b buildings, b × n 1 nodes) = T (1, 1 × n 1 )
With n1: the optimal number of nodes to simulate 1 building

Our previous experiments has shown the optimal distribution for one cl5 building requires 4 nodes, so we simulated 2, 3 and 4 buildings on 2 × 4, 3 × 4 and 4 × 4 nodes. Figure 11 shows this approach has fully succeeded for the cl5 benchmark. We obtained approximately the same execution time: 53s on our Nehalem cluster with an Eth 1Gb/s interconnect and 45s on At the opposite, the cl3 benchmark includes significant amount of communications compared to the amount of computations. We replicated our 3 nodes best distribution pattern, and we simulated 2, 3, 4 and 5 buildings on 2 × 3, 3 × 3, 4 × 3 and 5 × 3 nodes. Figure 12 shows we obtained unperfect but interesting performances on our 10Gb/s cluster (bottom thick curve). Multi-simulation times increase by +31%: from 58s up to 76s when processing 1 building on 3 nodes up to 5 buildings on 15 nodes. On our 1Gb/s cluster the times increase by +72% (top thick curve). This experiment shows the impact of the cluster network bandwidth, but we also observed the FMUs exchange mainly small messages. So, a high bandwidth and a low latency seem both necessary to achieve good size up of DACCOSIM simulations.

Thin curves on figure 12 show the execution time when the multi-simulation results are not stored. No result storage means no asynchronous IO thread running in parallel with all others on each node, and no time spent writing on disks. On cl3 benchmark it leads to smaller execution times but size up benchmarks still exhibit regular increase of execution time. So, we can probably improve our result storage mechanism, but it is not the main limitation of our size up benchmark.

Finally, figures 13 and 14 exhibit multi-simulation times of cl3 benchmark for 4, 5 and 8 buildings, when increasing the number of computing nodes. We replicated our best distribution pattern (1 building on 3 nodes for cl3 benchmark), to deploy for example 4 buildings on 12 = 4 × 3 nodes. Then we grouped 2 buildings per 3-node block to deploy 4 buildings on 6 nodes, and finally we grouped the 4 buildings on 3 nodes. So, all our deployments are based on the previously identified optimal deployment of 1 building. This straightforward approach does not track a global optimized deployment, but allows to quickly configure large scale experiments.

Logarithmic scales on both axis lead to expect straight lines with -1 slope for ideal distributions, corresponding to theoretical times T (n nodes) = T (1 node)/n. Experimental time curves roughly look like straight lines, but the slope appears stronger on 10Gb/s cluster, showing again the impact of the interconnect performance on our distributed multisimulations. When simulating 4 or 5 buildings (middle and bottom curves), it was bearable to run and measure execution times on 1 multi-core node and to compute speedup. On 10Gb/s cluster we achieved significant speedup close to 7 on 12 nodes and to 9.7 on 15 nodes (see figure 14), compared to multi-threaded executions on one node. When simulating 8 buildings we did not waste time to run long computations on one node: large problems are not intended for mononode executions (and perhaps could not fit into one node memory). But we observe the 8-building curves are similar to the 4-building ones with better decrease. So, DACCOSIM can efficiently run larger problems on greater number of nodes and scale our multi-simulations (with a fast cluster network).

Executions on 6 nodes (for 4 and 8 buildings) should be 2 times longer than executions on 12 nodes. Computing ressources are 2 times less numerous, but communication achievement is also different. The BC FMU is connected to 6 local FMUs instead of 3 and achieves 2 times more intra-node communications, and the network switch is solicited on 6 ports instead of 12. Deep investigations would be required to analyse the origin of these abnormally good performances on 6 nodes, but they point out our straightforward deployment (replicating the 1-building best one) is not optimal. A global deployment of b buildings (i.e. 1+10×b FMUs) on n multi-core nodes could achieve better performances, but is not currently available.

V. CONCLUSION AND FUTURE WORK

Our FMI based multi-simulation environment (DAC-COSIM) is multithreaded and distributed, and allows the user to distribute FMU graphs on clusters of multi-core nodes. According to our distribution strategies, our experiments have achieved significant speedup and satisfying size up on an Ethernet 10Gb/s cluster of hexa-cores, running up to 81 FMUs and using up to 16 nodes to simulate heat transfer inside buildings. Complementary large benchmarks replicating an elementary distribution have shown DACCOSIM can scale (efficiently running larger problems on more computing nodes). In order to improve DACCOSIM performances, our next works will consist in: (1) exchanging more compact messages, (2) using high performance MPI communication library instead of 0MQ middleware, (3) exploiting low latency and high bandwidth Infiniband cluster network by relying on MPI.

Automatic distribution of the FMU graph is also under investigation. Numerous experiments have allowed to define a first distribution strategy, that we aim to automatize. DAC-COSIM should be able to point out the ideal number of nodes and the associated distribution of the FMU graph, and also to compute the optimal distribution when the number of nodes is constrained by the cluster availability.

Finally, the use case considered in this paper was based on multi-floor buildings with equation-based models only. In the future, EDF R&D intends to run hybrid multi-simulations involving cyber and physical models mixing piecewise constant signals and continuous time signals. As FMI-CS has been designed only for continuous time signals, some extensions to the standard FMI are required and some proposals are on the table (see [START_REF] Cremona | FIDE -An FMI Integrated Design Environment[END_REF] from UC Berkeley). EDF also leads a French working group on this topic with partners (CentraleSupélec, Dassault Systèmes and CEA List). From the point of view of the optimal distribution of a multi-simulation on clusters, the control FMUs will certainly be very light and ought to be handled by the deployment algorithm without disturbing the distribution of the heavily calculative FMUs. ACKNOWLEDGMENT Authors thank Region Lorraine for its constant support.

Fig. 1 .

 1 Fig. 1. DACCOSIM GUI: FMU graph definition and configuration

Fig. 2 .Fig. 3 .

 23 Fig. 2. Distributed DACCOSIM architecture, with a hierarchical master

Fig. 4 .

 4 Fig. 4. Two orchestration modes of a multi-simulation step

Fig. 5 .Fig. 6 .

 56 Fig. 5. Best distribution (4 nodes) of cl5 benchmark with 1 building

Fig. 7 .Fig. 8 .Fig. 9 .Fig. 10 .

 78910 Fig. 7. cl3 benchmark on one Nehalem node

Fig. 11 .

 11 Fig. 11. Size up experiments on cl5 benchmark

Fig. 13 .

 13 Fig. 13. cl3 benchmark time with 4, 5 and 8 buildings on 1Gb/s cluster

http://www.supelec.fr/342_p_36889/risegrid.html

SIANI, University of Las Palmas, Spain: JavaFMI (2016), https://bitbucket.org/siani/javafmi/wiki/Home

https://oar.imag.fr/


G Le Quilliec

H.-P Lieurade

M Bousseau

M Drissi-Habti

G Inglebert

P Macquet

L Jubin

Fatigue Behaviour of Welded Joints Treated by High Frequency Hammer Peening: Part I, Experimental Study

Keywords: high frequency hammer peening, hammering, fatigue testing, steel welded assemblies, experimental investigation, post-weld treatment

High frequency hammer peening is a recent improvement method which is probably one of the most effective for treating welded assemblies. A number of experimental results relating to this process are presented in this article. These results lead to better understand the mechanisms of the process, to outline the influence of the operating parameters and to confirm the role played by the initial quality of the welds. In the long run, the aim of this study is to propose an industrially applicable approach for the fatigue design of welded assemblies treated by high frequency hammer peening. The numerical study carried out in parallel is presented in the IIW document: XIII-2395-11.

Introduction

In comparison with the base and filler materials used, welded structures have relatively low fatigue strength. This can be explained by the number of cycles required for a micro-crack to initiate. The initiation stage is often predominant in the fatigue life of metallic materials. It becomes very weak in the specific case of welded assemblies. This is mainly caused by the usual presence of residual tensile stresses, high stress raisers and an often poor surface roughness at the weld toe. It is therefore essential to have a design that takes into account the specific features in the case of exposure to cyclic loading. In addition, the precise study of the loading that will be undergone by the welded assembly will enable better design of the welded connection and make it possible to justify whether or not post treatments should be applied. These treatments will make possible to save weight. This gain will be all the more significant since the stresses will be limited and ultimate tensile strength will be high.

One of the most recent (and probably one of the most efficient) improvement methods is high frequency hammer peening or HFH [START_REF] Le Quilliec | Traitement des assemblages soudés par martelage haute fréquence[END_REF]. With respect to conventional hammer peening, HFH stands out thanks to its improved simplicity and user comfort. Applying this treatment has been shown to be very competitive making the most of the advantages provided by certain materials in the field of welding and more particularly high strength steels. However, there is a certain number of technical difficulties that must be overcome to apply HFH whether for new or old structures. Using the various results available for this process is uneasy for the fatigue design of an HFH treated welded assembly for a given configuration and treatment parameters. Moreover, premature failures are often observed and the cause of these early failures is not clearly identified.

The work presented in this article falls within the framework of a project which should in the long run propose an industrially applicable approach with the aim of increasing the guaranteed fatigue life of new structures and to extend the one of existing structures. The experimental study of the HFH carried out on a cruciform welded assembly will be presented in the body of the text. The numerical study carried out in parallel is presented in the document IIW: XIII-2395-11.

2 Analysis of the as-welded state of the selected welded assembly

Configuration of the welded assembly

Cruciform non-carrying-fillet welded joints are one of the most common structural details of welded structures as well as the most susceptible to fatigue phenomenon. This assembly configuration was therefore selected for the HFH study in S355K2 grade rolled structural steel, a grade that is just as commonly used (as per standard NF EN 10025-2, 2005).

According to this standard, the yield strength of the base metal at ambient temperature in the rolling direction must be at least 355 MPa for an ultimate tensile strength ranging between 470 and 630 MPa. The thickness of the sheet used for the main plate and the attachments is 15 mm (Figure 1). The weld test coupons manufactured are each 1 m long and 250 mm wide with a total height at the attachments of 115 mm. The edges were not prepared in any special way prior to one-pass welding (cross passes) with a flux-cored type filler metal which has a similar grade (T42 4 MC 2H5) to that of the base metal. Finally, all the assemblies were performed with the automatic MAG process which makes it possible to ensure a good quality as well as proper repeatability of the operation on each run. The welding operation was carried out in flat position which results in run profiles that have a higher quality with regard to fatigue strength.

Main plate

Attachments

Geometrical measurements of the as-welded weld run

The detailed geometry of welded assemblies was analysed using a laser triangulation portable dimensional measuring device.

The overall geometry of a welded assembly is quickly characterized and inspected with this technology. Figure 2 illustrates the position of the system during measurement.

The statistical analysis carried out after measurements on all the welded assemblies is presented in table 1 (a detailed definition of the various parameters is given in figure 3). In accordance with the even appearance of these assemblies, the standard deviations obtained reveal a low dispersion of the values. The mean toe radius, given the position of the weld is greater than standard quality, of approximately 1.5 mm in MAG welding. The depth of the undercut measured, which is of approximately the same magnitude of resolution of the measuring device (0.05 mm), is insignificant.

Monotonic tensile tests

Standard NF EN 10025-2 [2005] only provides the ranges of the main mechanical characteristics that must be fulfilled for the selected S355K2 grade steel. Consequently, monotonic tensile tests were performed on the base metal used as well as on weld metal specimens. Welded samples has been obtained from a butt-welded assembly which has been especially carried out with the same base and filler materials as for the cruciform assembly. As opposed to the cruciform assemblies, the single-V butt-welded assemblies allow samples of traditional sized tensile specimens to be taken. However, the thermal cycle as well as the dilution phenomena between the items to be assembled and the filler product vary depending on the configuration. The mechanical behaviour obtained in the weld zone of the butt-welded assembly consequently differs from that of the cruciform assembly. However it is worth noting that this validation "practice" is common in the welding field.

According to the engineering tensile curves presented in figure 4, the Young's modulus measured is approximately 205 GPa for the base metal and 210 GPa for the weld metal. The offset yield strength at 0.2 % of the strain (R p0.2 ) and the ultimate tensile strength are respectively 386 and 541 MPa for the base metal whereas those of the weld metal are 407 and 512 MPa. The presence of a significant Lüders plateau should also be noted for the weld metal.

Low-cycle fatigue tests

In addition to the monotonic tensile tests, low-cycle fatigue tests were carried out on the base metal. The results obtained will be used during high cycle fatigue design of the welded assemblies. In addition, the cyclic tensile tests performed make it possible to characterize the type of hardening, mainly kinematic in this case as shown in figure 5.

All the cyclic properties estimated from the experimental tests are given in table 2. The values estimated by the Universal Material Law (UML), proposed by [START_REF] Bäumel | Materials Data for Cyclic Loading -Supplement I[END_REF], and the medians method, proposed by [START_REF] Meggiolaro | Statistical evaluation of strain-life fatigue crack initiation predictions[END_REF], are also given. These two methods are based on statistical studies of many experimental test results. The Manson-Coffin curves obtained are presented in figure 6.

σ f ε f b c K n [MPa] [%]
[MPa] According to experimental results 702.5 27.2 -7.82 • 10 -2 -5.

HFH treatment device

HFH treatment equipment

The Esonix-type UIT (Ultrasonic Impact Treatment) device was selected to perform high frequency hammer peening treatments (Figure 7). This choice was mainly motivated by the many results given in the literature relating to this hammer model. Comparisons may be made at the end of the treatment.

Power supply and excitation frequency

The tool is powered and controlled by an alternating current generator, Esonix PLC05 type. A water cooling device is used to limit transducer heating during the treatment. The generator has 5 power levels which vary the amplitude of the sonotrode vibrations. Whatever is the level of power selected, the excitation frequency remains equal to the resonance frequency of the transducer that is, 27 kHz for the model used.

Indenters

The set of indenters crosses a rubber strip inserted in the removable support in order to keep the indenters in their casing when the tool is inactive (Figure 8).

Indenters with 3 mm in diameter are recommended for the treatment of steel welded assemblies. The detail of their geometry was measured and presented in figure 9. The measurements carried out indicate that the rounded edge has a main radius of approximately 2.4 mm and a secondary blend radius of approximately 0.86 mm. On average, each indenter weights 1.37 g and is 25.6 mm long. The geometry measured will be reused during the numerical simulation of the treatment.

Contact force

The transducer, sonotrode and indenter assembly (hammer peening device) is suspended in the casing in order to isolate the operator from residual vibrations. The suspension spring is pre-stressed upon assembly so that there needs to be a certain pressure for the hammer peening system to be released from stop position during the treatment Depending on the position of the tool, its weight acts in addition to the pressure applied by the operator. The mean contact force of the indenters on the treated surface was measured for the following three directions:

• In vertical position, indenters directed downwards, the force varies between 60 and 70 N depending on whether the spring is close to the lower or upper stop;

• In horizontal position, it varies between 35 and 45 N, and the weight is no longer involved in the process;

• Finally, in vertical position, indenters directed upwards, the weight acts in opposition and the mean contact force only varies between 10 and 20 N.

Configuration of the treatment device

A specific assembly was adapted on a conventional milling machine in order to best control the treatment parameters (Figure 10). The tool is assembled vertically above the part to be treated which may be flanged in a sloped position on the moving table . A mobile support guides the hammer in vertical translation with a low stop. Depending on the vertical position of the table, this low stop helps to adjust and constantly maintain the contact force value. In addition to the direction of the tool or the part and the contact force, this assembly helps to simply control the pass speed by taking advantage of the automatic feed system of the milling machine table.

Standard parameters and treatment conditions

The high frequency hammer peening treatment tests were carried out in various configurations. Only one parameter varies at a time with the following standard values for the stationary parameters:

• Vertical position, indenters directed downwards;

• Contact force of the indenters of 70 N;

• Feed rate of 400 mm/min;

• One single pass;

• Maximum power level (5/5);

• 23°angle between the tool axis and the normal at the assembly main plate;

• Four 3 mm diameter indenters (or three 4.8 mm diameter indenters) aligned in the treatment direction for the tests on the base metal;

• Reduced to only two 3 mm diameter indenters (or one 4.8 mm diameter indenter) for the welded assemblies.

The number of indenters for the treatment of welded assemblies was decreased due to the difficulties experienced in following the profile of the weld toes. There is greater freedom of movement in the vertical direction obtained using only 2 indenters and no untreated zone was observed.

It is worth noting that the indenters are consumables. The 3 mm diameter indenters were renewed on a regular basis in order to prevent any failure during treatment. Otherwise, if they are used for an extended period this would cause either their strain for low feed rates stemming from a greater heating (causing a blue stain at the end of the indenters through surface oxidation), or fatigue failure in their housing for higher feed rates (refer to figure 11).

Finally, ear protection is recommended in order to attenuate the noise generated by the ultrasonic transducer and by the repeated impacts of the indenters between the sonotrode and the surface of the part to be treated. For information, a 96 dBA sound level was recorded at a distance of 50 cm from the edge of the tool during treatment. This level is higher than that expected based on the literature. 4 Analysing the effects of high frequency hammer peening

Aspect of the grooves treated

The general aspect of the grooves treated was measured for a set of base metal samples. Figure 12 shows the development of the groove for 3 different treatment power levels. Hammer peening striations appear whatever is the power used. For a greater hammering power, a border is formed at the edge of the hammered zone given by the greater depths of penetration reached (case c). This border expands at a lower treatment speed (case d).

Furthermore, in the case of 2 passes at 315 mm/min, damage is created in the hollow of the groove. The one pass treatment at 160 mm/min which corresponds to a theoretical coverage rate closer to case f does not feature this damage, however there is a heightened groove depth. These differences may stem from the fact that the heating has less time to be dissipated at lower treatment speed. Finally, the aspect of the groove obtained with the 4.8 mm diameter indenters is given for example. The kinetic energy transmitted to the indenters is then distributed over a larger contact surface which tends to explain the moderate depth of the groove obtained in this last case.

The operation was then repeated on the cruciform welded assemblies and presented in figure 13. The striations marked by hammer peening can be found on every grooves. The borders are less pronounced even at lower treatment speed. The hardness at the weld toe may explain the limited formation of borders and the moderate depth of the grooves obtained in comparison to the base metal which is softer. It can also be noted that at low power (case a), the initial weld toe is still present in places. Furthermore, a "flaking" type damage appears on occasions and more specifically for cases d and f . A lower coverage rate (case c et e) is therefore to be recommended in order to limit the flaking phenomenon on surface. Finally, it should be noted that the treatment with the 4.8 mm diameter indenters results in a significant improvement of the weld toe radius. However, the treatment depth obtained remains limited in this last case. An excessively closed toe angle may prevent this type of indenters from reaching the weld toe itself.

Profiles of the grooves and surface roughness

Measurements by mechanical contact

A first series measurements was carried out mechanical contact using a "Somicronic-Surfascan 35" profilometer. These measurements were performed on a length varying from 1.6 to 5 mm based on the accessibility of the probed zone at the rate of one measurement point every 2 µm. Various profiles measured in the transverse direction to the groove on a welded assembly hammered with the standard parameters are listed in figure 14. Other measurements were carried out in the longitudinal direction of hammer peening on this same sample and recorded in figure 15. The results obtained reveal that there is a certain uniformity for each of the two measurement directions:

• In the hammer peening direction, the shape defects observed correspond to striation separating each of the imprints left by the indenters;

• In the transverse direction, the general profile obtained is much more even. It is very similar to the geometry of the indenters sloping in the base metal direction corresponding to the direction of the tool upon treatment.

From a macroscopic standpoint, it can be noted that the groove radius in the zone with greater stress raiser corresponds more to the blend radius of the indenters than their main tip radius (Figure 9).

On a smaller scale, the roughness deduced from the previous data is given in table 3. Maximum roughness height R t remains generally less than 4 µm in the transverse direction to the groove and ranges between 10 and 20 µm in the collinear direction. These values associated with the blend radius obtained at the weld toe are especially interesting with respect to the fatigue strength of welded assemblies, given the irregularities usually present in this zone.

No Tab. 3: Transverse (1 to 4) and longitudinal (5 to 7) roughness measurements in the groove of a weld toe hammered with the standard parameters.

Optical measurements

• Methodology:

An optical measurement analysis was also carried out using a "Stil-Micromesure 2" station. The optical pen selected enables lateral precision less than 0.1 µm for a vertical precision a few nanometres. However, the measurement range is limited to 300 µm in depth and a maximum slope of the scanned surface of +/-28°. Given that the slope at each edge of the hammered groove is too steep to be detected by the selected optical pen, the measurement zone was limited to the centre of the groove over a width of 2 mm (width of the grove: 3 mm). The measurement spacing was set to 1 µm in both directions.

The three-dimensional surface obtained on a base metal sample (hammered with the standard parameters) is given in figure 16. Impressions 1 to 9 correspond to those of the last of the 4 indenters which contributed to generate the groove. • Estimation of the hammer peening frequency:

Whereas a constant feed rate was dictated at 400 mm/min for this sample, the impacts generated by the treatment device are visibly not quite periodic. The zone of index 0 separating the 3 rd and 4 th impacts corresponds moreover to a space that was not hammered by the last indenter.

All these impacts make it possible to estimate a mean hammer peening frequency of approximately 48 Hz for the analysed sample by considering:

f moy = N × v × n l (1)
With:

N the number of indenters n the number of impacts identified v the feed rate (in m/s) l the analysed length (in m)

In the case of the welded assembly treated with the standard parameters, an optical microscopic analysis helped to identify a set of 49 impacts over an analysed length of 3 cm (that is, a distance between impacts of approximately 60 µm). By applying the relation proposed ( 1), we obtained a mean hammer peening frequency of 21.8 Hz (for 2 indenters used). This result is similar to that estimated on the base metal from figure 16 if only the frequency per indenter is considered, that is 12 Hz on the base metal and 10.9 Hz on the welded assembly. These hammer peening frequencies are however markedly less than those given in the literature. For example, the frequency measured by [START_REF] Statnikov | Comparison of the efficiency of 27, 36 and 44 kHz UIT tools[END_REF] range between 200 and 350 Hz for 4 indenters, that is between 50 and 87.5 Hz per indenter.

Microstructure and surface defects caused

In addition to the peeling (or flaking) observed on surface for high coverage rates, figure 17 reveals the presence of material folds which appear more occasionally and essentially for higher treatment powers (levels 4 and 5). Figure 18 gives an example of the microstructure identified in the treated weld toe. As observed in the literature, the grains are moderately elongated in depth and more pronounced on surface. It should be noted that the depth of the groove caused by the tool for the standard parameters is approximately 0.2 mm.

Residual stresses

X-ray diffraction measurements were carried out using a chrome tube and a linear detector (making it possible to register the radiation of the chrome for the family of lattice planes (211) of phase α of the steel), on 3 different specimens:

• Base metal hammered along a line with standard parameters (Figure 19);

• Cruciform welded assemblies treated at the weld toe with standard parameters (Figure 20);

• Cruciform welded assemblies treated then submitted to fatigue loading before the residual stress measurements (Figure 21).

The collimator used leads to an irradiated volume of 2.5 (longitudinally) × 1.0 (transversally) × 0.006 (depth) mm 3 . The decision was made not to carry out any specific preparation of the groove surface for these measurements except for cleaning with hydrochloric acid. The analyses were carried out along the longitudinal and transverse directions of the groove.

The results, presented in figures 19, 20 and 21, are only qualitative given the non-flat geometry of the groove and its surface roughness (striations) which affects the accuracy of the results. However it can be observed that:

• Relatively high compressive stresses are present on the edge of the groove in both directions;

• The longitudinal stresses in the hollow of the groove are strongly compressed (between -300 and -600 MPa) whereas the transverse stresses are subject to low compression;

• Fatigue loading does not result in significant relaxation of the transverse residual stresses located in the hollow of the groove. 5 Fatigue tests

Configuration of the tests specimens

After treatment, the weld test coupons (Figure 1) were cut and machined in order to make the fatigue test specimens. The geometry of the specimens is specified on figure 22.

Figure 22: Dimensions of the fatigue specimens (in mm).

Fatigue test results

Reverse tensile tests (R = -1) were carried out on the as-welded cruciform welded assemblies or after hammer peening treatment. Figure 23 highlights the satisfactory fatigue strength of the as-welded assembly, in compliance with the previously mentioned welding quality, as well as the remarkable increase of fatigue life in the case of the treated assemblies, especially for long lives.

On the same figure, the base metal fatigue curve at 50 % of failure probability, estimated by extrapolation of the low cycle fatigue curve, was presented in red doted line. In theory, this curve represents the optimum fatigue strength which may be reached by the treated welded assemblies. As a matter of fact, if treatment gives rise to greater fatigue strength of the hammered zone than that of the base metal, the latter becomes the dimensioning element of the assembly.

Fracture surface

The stress ratio R = -1 results in a caulking of the crack surfaces which restricts analysis. However on figure 24, several radial lines can be seen around the hammered surface. These radial lines are due to multiple initiations and joining of micro cracks in the hollow of the groove caused by hammer peening. For the specimens broken outside of the hammered zone, a single initiation point was noted. This initiation point is located on the edge at the end of the calibrated zone (moderate stress raiser location).

Results of additional tests

An additional test campaign was carried out on the specimens treated either by manual process or the milling machine by using the same hammer peening device. The treatment was applied in two cases:

• Cruciform welded assemblies composed of DH36 steel (grade for marine application equivalent to S355);

• Butt-welded assemblies composed of S690QL steel.

DH36 cruciform welded assemblies

It can be noted on figure 25 that the improved fatigue strength of the treated assemblies (blue points) is not significant in comparison to the results obtained on the as-welded assemblies (black triangles). However, a significantly better gain for the manually treated assemblies can be noted with a tool direction and treatment speed as close as possible to that of treatments on the milling machine.

The slight improvement obtained on these cruciform assemblies can partly be explained by the quality of the weld and the relatively low connecting angle at weld toe.

mm

Figure 25: Fatigue results at on cruciform joints composed of DH36 steel for treatment on a milling machine (top) and manually (bottom) [Bousseau, 2010a].

S690QL butt-welded assemblies

In this case as well, a slight improvement is noted on the treated assemblies (red points and blue triangles on figure 26) compared with the results of the as-welded assembly (black circles). Here, the connecting angle at weld toe and the initial general quality of the welded assembly does not explain the low gain obtained. However, the analysis of the fracture surface highlights a crack initiation located at the bottom of hammer peening folds (Figures27 et 28). Figure 28 in particular reveals a partially propagated crack, initiating once again at a fold at the bottom of the groove.

The formation of a hammer peening fold can be explained through a false direction of the tool at the time of treatment of the butt-welded assemblies. As a matter of fact, in the case in question, the angle of incidence of the indenters was set at 6°with respect to the normal at the assembly sheet. In this direction, the indenters first impact the weld toe from the side of the weld zone and fold a part of the metal under the indenter itself, in the stress raiser zone, creating the folds observed.

Figure 1 :

 1 Figure 1: Selected configuration of the cruciform joints.

Figure 2 :

 2 Figure 2: Laser triangulation dimensional analysis (WISC process).

Figure 3 :

 3 Figure 3: Definitions of the various parameters characterizing the geometry of the weld run.

Figure 4 :

 4 Figure 4: Engineering tensile curves of the base metal and the weld metal.

Figure 6 :

 6 Figure 6: Experimental and estimated Manson-Coffin curves.

Figure 7 :

 7 Figure 7: High frequency hammer peening tool used.

Figure 8 :

 8 Figure 8: Detail of the indenter support.

Figure 9 :

 9 Figure 9: Measured geometry of the 3 mm diameter indenters.

Figure 10 :

 10 Figure 10: Semi-automatic treatment device.

Figure 11 :

 11 Figure 11: Fracture surface of a 3 mm indenter.

  Figure 12: Change of the hammered groove on the base metal for various treatment parameters.

Figure 13 :

 13 Figure 13: Change of the hammered groove on the welded assembly for the various treatment parameters.

Figure 14 :Figure 15 :

 1415 Figure 14: Profiles measured transversally in the groove (case of a welded assembly).

Figure 16 :

 16 Figure 16: 3D surface and depth level (along z), expressed in false colours of the central zone of a groove from a base metal hammered.

Figure 17 :

 17 Figure 17: Fold type defects on the base metal sample treated with standard parameters.

Figure 18 :

 18 Figure 18: General microstructure of the weld zone treated with standard parameters and magnification at the impacted surface on base metal side.

Figure 19 :Figure 20 :Figure 21 :

 192021 Figure 19: Transverse (left hand side) and longitudinal (right hand side) residual stresses with respect to the groove, measured on a base metal sample hammered with the standard parameters.

  Figure 23: Fatigue results of the welded joints at R = -1 and related estimations.

Figure 24 :

 24 Figure 24: Fracture surface of a welded assembly treated by HFH and broken in the hammered zone.

Figure 26 :Figure 27 :

 2627 Figure26: Fatigue results at R = 0.1 on butt-welded joints composed of S690QL for treatment on a milling machine[Bousseau, 2010b].

Figure 28 :

 28 Figure 28: Microstructure of a S690QL butt-welded assembly treated by HFH and detail of the crack initiation zone at the bottom of the groove.

  Cyclic tensile curve of the base metal (5 cycles at ∆ε t = 1.209 %).

	σ [MPa]
	400
	200
	0
	-200
	-400						ε [%]
	0	0.2	0.4	0.6	0.8	1.0	1.2
	10 -1 10 0 10 2 Figure 5: 10 -2 ∆ε/2 [%]	10 3		10 4		10 5	10 6	N
							51 • 10 -1 845.3 0.142
	UML		812	59	-8.7 • 10 -2	-5.8 • 10 -1	893.1	0.15
	Medians method		812	45	-9 • 10 -2

-5.9 • 10 -1 --Tab. 2: Estimated cyclic properties of the base metal. Trend curve ∆ε t /2 Trend curve ∆ε e /2 Trend curve ∆ε p /2 ∆ε t /2 estimated from the UML ∆ε t /2 estimated from the medians method

Conclusions

The aim of this study was to better understand the effect of HFH parameters on the fatigue life improvement provided of steel welded assemblies. The following points are worthy to be noted:

• With optimized hammer peening treatment conditions, the results obtained are well in line with those of the literature;

• Even though the welded assemblies have good behaviour in the as-welded state, a significant improvement in fatigue strength was noted after treatment;

• The residual compressive stresses measured in the groove created by HFH are relatively low in the loading direction;

• The improvement of the local geometry and moving the initiation point with respect to the initial weld toe are significant parameters;

• The reduced number of indenters (2 instead of 4 normally used) does not seem to have affected these results;

• The estimated hammer peening frequency is markedly lower than those indicated in literature;

• The results presented on several assembly batches seem to highlight the importance, on the one hand, of the initial weld toe geometry and, on the other hand, the conditions of use of the tool (direction, pressure, etc.) in relation to the risk of hammer peening folds appearing.

This study therefore reveals the importance of developing a consistent and robust approach with the aim of defining the optimum operating conditions of the process.


J Serres

T Raharijaona

E Vanhoutte

F Ruffier

Event-based visual guidance inspired by honeybees in a 3D tapered tunnel

In view of neuro-ethological findings on honeybees and our previously developed vision-based autopilot, in-silico experiments were performed in which a "simulated bee" was make to travel along a doubly tapering tunnel including for the first time event-based controllers. The "simulated bee" was equipped with:

• a minimalistic compound eye comprising 10 local motion sensors measuring the optic flow magnitude, • two optic flow regulators updating the control signals whenever specific optic flow criteria changed, • and three event-based controllers taking into account the error signals, each one in charge of its own translational dynamics. A MORSE/Blender based simulator-engine delivered what each of 20 "simulated photoreceptors" saw in the tunnel lined with high resolution natural 2D images. The "simulated bee" managed to travel safely along the doubly tapering tunnel without requiring any speed or distance measurements, using only a Gibsonian point of view, by:

• concomitantly adjusting the side thrust, vertical lift and forward thrust whenever a change was detected on the optic flow-based signal errors, • avoiding collisions with the surface of the doubly tapering tunnel and decreasing or increasing its speed, depending on the clutter rate perceived by motion sensors.

I. INTRODUCTION

Winged insects are able to navigate through unfamiliar environments without any need for sensory systems such as Global Positioning Systems (GPS) or Inertial Measurement Units (IMUs), or emissive proximity sensors such as ultrasonic or laser range finders, radar, or scanning LIght Detecting And Ranging (LIDAR) sensors, using only the Optic Flow (OF) ( [START_REF] Gibson | The perception of the visual world[END_REF]; [START_REF] Nakayama | Optical velocity patterns, velocitysensitive neurons, and space perception: a hypothesis[END_REF]; [START_REF] Koenderink | Facts on optic flow[END_REF]). Flying insects make use of the OF to jointly avoid lateral obstacles and control their speed ( [START_REF] Srinivasan | Honeybees as a model for the study of visually guided flight, navigation, and biologically inspired robotics[END_REF]; [START_REF] Roubieu | A biomimetic vision-based hovercraft accounts for bees complex behaviour in various corridors[END_REF]), and to cruise and land ( [START_REF] Franceschini | A bio-inspired flying robot sheds light on insect piloting abilities[END_REF]; [START_REF] Srinivasan | Honeybees as a model for the study of visually guided flight, navigation, and biologically inspired robotics[END_REF]; [START_REF] Expert | Flying over uneven moving terrain based on optic-flow cues without any need for reference frames or accelerometers[END_REF]). Electrophysiological studies on flies motion sensitive neurons have shown that OF detection is texture-and contrastindependent [START_REF] Franceschini | Directionally selective motion detection by insect neurons[END_REF], and can be explained by a "time-of-travel" scheme ( [START_REF] Lichtsteiner | A 128× 128 120 db 15 µs latency asynchronous temporal contrast vision sensor[END_REF]; [START_REF] Benosman | Asynchronous frameless event-based optical flow[END_REF]; [START_REF] Sabiron | Low-speed optic-flow sensor onboard an unmanned helicopter flying outside over fields[END_REF]; Fig. 6 in [START_REF] Floreano | Miniature curved artificial compound eyes[END_REF]), which differs from the Hassenstein-Reichardt model originally developed in 1956 [START_REF] Hassenstein | Systemtheoretische analyse der zeit-, reihenfolgen-und vorzeichenauswertung bei der bewegungsperzeption des rüsselkäfers chlorophanus[END_REF]. The output of "time-of-travel"-based OF sensors can be only updated when a change in the optical contrast with time is detected. An asynchronous signal is therefore generated in the output of a "time-of-travel"-based OF sensor, which can be conveniently used in event-based control systems [START_REF] Durand | Further results on event-based PID controller[END_REF]. Because of the slight delay in the camera read-out, event-based collision avoidance algorithms have been only developed so far under open loop conditions in slightly cluttered corridors, but they were expected to be ready for use in the near future in closed loop on-board robotic platforms ( [START_REF] Clady | Asynchronous visual event-based time-to-contact[END_REF]; [START_REF] Milde | Bioinspired event-driven collision avoidance algorithm based on optic flow[END_REF]). Event-based control systems have also been used for visual odometric purposes in the field of robotics [START_REF] Censi | Low-latency event-based visual odometry[END_REF].

In the present paper, we present for the first time an event-based version of the ALIS autopilot (ALIS stands for "AutopiLot using an Insect-based vision System"), which was originally designed with fixed-step digital controllers to perform tunnel-following tasks [START_REF] Portelli | Modelling honeybee visual guidance in a 3-D environment[END_REF]. The ALIS autopilot is based on the principle of dual OF regulators ( [START_REF] Portelli | Modelling honeybee visual guidance in a 3-D environment[END_REF], [START_REF] Roubieu | A biomimetic vision-based hovercraft accounts for bees complex behaviour in various corridors[END_REF], [START_REF] Expert | Flying over uneven moving terrain based on optic-flow cues without any need for reference frames or accelerometers[END_REF]) using ventral, dorsal, and lateral OFs. In section 2, the set-up of the simulation and the simulated agent mimicking a honeybee equipped with a minimalistic visual system including only 10 local OF measurements are described. In section 3, we describe the trajectories of the simulated agent fitted with the event-based ALIS autopilot in a doubly tapering tunnel (with a tapering angle of 7 • in both planes). The results of the simulations performed show that the event-based ALIS autopilot generates a small number of events making the "simulated bee" travel along the tapered tunnel, thus greatly decreasing the number of control signal computations required.

II. SIMULATION SET-UP

In-silico experiments were carried out on two standard computers running: (i) a first one equipped with the MATLAB/Simulink c software program at a sampling rate of 1kHz for OF sensing purposes and feeding the eventbased control systems, and (ii) a second one equipped with the MORSE software program (MORSE stands for "Modular OpenRobots Simulation Engine" [19], [START_REF] Echeverria | Modular openrobots simulation engine: Morse[END_REF]) to simulate the visual interactions with the virtual 3D environment. The two computers were linked together by a TCP/IP protocol in which the maximum delay was limited to 3ms. The virtual 3D environment was built with the Blender software program [START_REF]Blender release 2.76 under GNU GPLv2[END_REF]. The four surfaces of the tunnel were lined with high resolution photographs of natural panoramic scenes [START_REF] Brinkworth | Biomimetic motion detection[END_REF]. The 20 scalar photoreceptor signals were transmitted via TCP/IP to a second computer processing the 10 local optic flow computations, the 3 event-based controllers and the 3 translational dynamics of the "simulated bee" using the Matlab/Simulink software program. The Matlab/Simulink program then sent the 3D position (X, Y, Z) of the "simulated bee" back to the MORSE simulator. bees 3D trajectories starting at the initial coordinates x 0 = 0.1m; y 0 = 0.5m; z 0 = 0.5m (in black), x 0 = 0.1m; y 0 = 0.4m; z 0 = 0.5m (in blue) and x 0 = 0.1m; y 0 = 0.5m; z 0 = 0.3m (in red) plotted every 300ms, (C) Trajectories in the vertical plane (x,z), and in (D) in the horizontal plane (x,y), plotted every 300ms. (E) Forward speed Vx profiles.

A. Optic flow

The OF vector field perceived by a simulated agent depends in particular on the structure of the 3D environment ([1]; [START_REF] Nakayama | Optical velocity patterns, velocitysensitive neurons, and space perception: a hypothesis[END_REF]; [START_REF] Koenderink | Facts on optic flow[END_REF]). The OF can be defined by a vector field of the apparent motion of objects, surfaces, and edges in a visual scene generated by the relative motion between an agent and the scene. The OF field ω (1) is a combination between two components: a translational ω T and a rotational OF ω R [START_REF] Koenderink | Facts on optic flow[END_REF].

ω = ω T + ω R (1)
It has been established that hymenopterans stabilize their gaze by compensating for any body rotations [START_REF] Viollet | Feed-forward and visual feedback control of head roll orientation in wasps (polistes humilis, vespidae, hymenoptera)[END_REF]. Since the simulated bee's head is assumed to be stabilized along the tunnel axis, each OF sensor will receive a purely translational OF ( ω R = 0). The translational OF (expressed in rad/s) can be defined as follows:

ω T = - V -( V • d) • d D ( 2
)
where d is a unit vector describing the viewing direction, V is the translational velocity vector, and D is the distance from the object seen by an OF sensor.

B. Simulated flying agent

The behaviour of the "simulated bee" introduced into the in-silico experiments was based on findings which may explain how a flying honeybee controls its speed and avoids obstacles ( [START_REF] Franceschini | A bio-inspired flying robot sheds light on insect piloting abilities[END_REF]; [START_REF] Srinivasan | Honeybees as a model for the study of visually guided flight, navigation, and biologically inspired robotics[END_REF]; [START_REF] Roubieu | A biomimetic vision-based hovercraft accounts for bees complex behaviour in various corridors[END_REF]); [START_REF] Expert | Flying over uneven moving terrain based on optic-flow cues without any need for reference frames or accelerometers[END_REF]). Each translational axis of the bee's dynamic model was uncoupled, as occurs in a quadrotor. Each of the bee's translational dynamics can be defined by a first order transfer function (see [START_REF] Portelli | Modelling honeybee visual guidance in a 3-D environment[END_REF] for details) with a time constant of 0.22s [START_REF] Franceschini | A bio-inspired flying robot sheds light on insect piloting abilities[END_REF]. Pitch angle, roll angle, and wing stroke amplitude were bounded on the basis of data previously published on bees [START_REF] Portelli | Modelling honeybee visual guidance in a 3-D environment[END_REF]. In this indoor study, the simulated agent was not subjected to any wind. In addition, the agent's rotational dynamics and the various mechanical couplings have been not considered in this work-in-progress paper yet.

The present agent was equipped with a set of 10 OF sensors: each of these sensors consisted of just two photoreceptors driving an Local Motion Sensor (LMS) based on a "time-of-travel" scheme [START_REF] Franceschini | Directionally selective motion detection by insect neurons[END_REF]. The visual axes of the two adjacent photoreceptors were assumed to be separated by the angle ∆φ, and each photoreceptor's angular sensitivity was assumed to be a Gaussoid function with an angular width at half height of ∆ρ = ∆φ = 4.2 • . 6 OF sensors were located in the horizontal plane oriented at azimuthal angles ϕ of ±45 • , ±90 • , and ±135 • , and the other 4 were oriented in the vertical plane at elevation angles θ of +45 • , ±90 • and -135 • .

III. EVENT-BASED DISCRETE TIME CONTROLLERS

In this section, we describe how we tuned the eventbased control approach introduced in [START_REF] Durand | Further results on event-based PID controller[END_REF] for OF regulation purposes. Three control input signals were computed, corresponding to the three uncoupled translational degrees of freedom on the surge (x), sway (y) and heave axes (z). We consider in this work as an "event", an eventrelated condition based on the error signals. Lead controllers are introduced into the heave and sway control systems to improve their stability. The natural integrator that relates the simulated agent's position to its speed makes both the heave and sway control systems reach a zero steady state error. However, a PI controller is introduced into the surge control system because OF is proportional to speed, but inversely proportional to distance [START_REF] Nakayama | Optical velocity patterns, velocitysensitive neurons, and space perception: a hypothesis[END_REF]. Consequently, we have to add an integral action to cancel the steady state error.

A. Event-based PI controller on the surge axis

A discrete time PI controller was designed. The modelling of the proportional part was quite straightforward, and the backward difference approximation method was used to model the integral part. The resulting code is: The event-based controller is tuned with the values K p = 20.10 -3 , K i = 6.10 -2 and T s = 1ms.

B. Event-based PD controllers on the sway and the heave axes

Two discrete time PD controllers were also designed to deal with the sway (y) and heave axes (z). Modelling the proportional part was quite straightforward, and the backward difference approximation method was used to model the derivative part. The following code was thus obtained for each controller: The event-based controllers on the sway axis were tuned with the values Kp = 20.10 -4 , T d = 0.1s, and on the heave axis, they were tuned with the values Kp = 20.10 -3 , T d = 0.1s. The input sign was defined on the basis of Fig. 2 adapted from [START_REF] Portelli | Modelling honeybee visual guidance in a 3-D environment[END_REF]. The parameter elim was chosen in order to significantly reduce the control updates, while making it possible to perform reference tracking of the OF

Side thrust

Direction of the largest OF Fig. 2. The control direction selector automatically selects the tunnel surface to be followed (wall, ground or roof) by multiplying the control signal (the output from the Positioning controller) by a direction factor u(3) that depends on the direction of the largest OF signal. Note that the sway and heave dynamics can be driven alternately, depending on which OF (side or vertical) is the largest at any given time. Adapted from [START_REF] Portelli | Modelling honeybee visual guidance in a 3-D environment[END_REF].

IV. RESULTS

In Fig. 1, the simulated environment was a doubly tapering tunnel 6m long, 1m wide, 1m high with a slope angle of 7 • . Fig. 1A gives a perspective view of the tunnel. The simulated bee entered the tunnel at zero speed, with the initial coordinates x 0 = 0.1m and various couples/pairs of y 0 and z 0 (Fig. 1B). Within the first 800ms, the simulated bee was controlled in the open loop mode to enable all the sensors measure the OFs. From the time 800ms onwards, the OF was regulated in the closed loop mode via the event-based controllers. Fig. 1C shows three trajectories in the vertical plane (x, z) and Fig. 1D shows these trajectories in the horizontal plane (x, y), plotted every 300ms. The simulated bee can be seen to have gradually increased both its ground clearance (Fig. 1C) and its right and left clearance (Fig. 1D), while the forward speed (Fig. 1E) increased automatically up to 1.5m/s. Fig. 3A gives the control input signal and the control updates on the surge axis with x 0 = 0.1m, y 0 = 0.5m and z 0 = 0.3m. The control input signals are also plotted on the sway and heave axes in Fig. 3A and3B. It can be seen that the number of control updates is greatly reduced in comparison with fixed-step digital controllers. Fig. 4A plots the OF with respect to the set-point of 540 • /s driving the surge dynamics. In Fig. 4B, the OF is plotted with respect the set-point of 315 • /s driving the sway and heave dynamics. Both systems of OF regulation gave satisfactory results.

V. CONCLUSIONS

In this paper, we present results of work-in-progress on an event-based autopilot mimicking the direct Optic Flow (OF) feedback control loops which may occur in honeybees. Likewise, honeybees tiny brains contain smart OF sensors, called Velocity-Tuned motion-sensitive descending neurons (or VT neurons): these are spiking neurons which are sensitive to the OF [START_REF] Ibbotson | Evidence for velocity-tuned motion-sensitive descending neurons in the honeybee[END_REF]. Flies are also equipped with a sensory system based on Lobula Plate Tangential Cells (LPTC), which are also spiking neurons coding the local motion detection signals triggered inside the insects brain [START_REF] Krapp | Estimation of self-motion by optic flow processing in single visual interneurons[END_REF]: each spike can be taken to constitute an event occurring in the insects neuronal system. Our proposed event-based control system is currently based on error signals, but in a near future it will be likely to become a major tool which can be used to understand more clearly how insects process the event-based visual information they receive in order to react appropriately: this information ranges from optical contrasts and local motion detection to self-induced wing movements generating the insects 3D movements in any tunnel configuration.

The results of the present in-silico experiments show that the OF-based regulatory strategy consisting in updating the control signals solely when visual contrasts are detected provides a suitable basis for traveling along an unknown tapered tunnel. It is now proposed to perform further studies involving a more sophisticated OF-based autopilot traveling along more complex tunnels, as well as to implement the event-based ALIS autopilot in the near future on-board a 380-g X4 MaG quadrotor [START_REF] Manecy | RT-MaG: An open-source SIMULINK toolbox for linux-based real-time robotic applications[END_REF].

Fig. 1 .

 1 Fig. 1. (A)Perspective view of the doubly tapering tunnel, (B) Simulated bees 3D trajectories starting at the initial coordinates x 0 = 0.1m; y 0 = 0.5m; z 0 = 0.5m (in black), x 0 = 0.1m; y 0 = 0.4m; z 0 = 0.5m (in blue) and x 0 = 0.1m; y 0 = 0.5m; z 0 = 0.3m (in red) plotted every 300ms, (C) Trajectories in the vertical plane (x,z), and in (D) in the horizontal plane (x,y), plotted every 300ms. (E) Forward speed Vx profiles.

  \% Inputs ysp = u(1);\% Optic flow SetPoint = 540 [deg/s] y = u(2);\% Optic flow Measurement e = ysp -y; \% Update control signal hact = hact + Ts;\% Ts=1e-3s if (abs(e-e_old)>elim) \% elim = 1 \% Calculate control signal up = Kp * e; ui = ui + Ki * hact * e end u = up + ui; \% Update e_old=e; y_old=y; hact=0;

  \% Inputs ysp = u(1); \% Optic flow SetPoint = 315 [deg/s] y = u(2); \% Optic flow Measurement sign = u(3); e = ysp -y; \% Update control signal hact = hact + Ts; \% Ts=1e-3s if (abs(e-e_old)>elim) \% elim = 1 \% N=100 \% Calculate control signal up= sign * Kp * e; ud= Td/(N * hact + Td) * ud -Kp * Td * N/(N * hact + Td) * (y -y_old); end u = up + sign * ud; \% Update e_old=e; y_old=y; hact=0;

Fig. 3 .

 3 Fig. 3. Simulated control signals and control updates during a 3D trajectory starting at the initial coordinates x 0 = 0.1m; y 0 = 0.5m; z 0 = 0.3m

Fig. 4 .

 4 Fig. 4. Simulated systems of OF regulation in the case of a 3D trajectory starting at the initial coordinates x 0 = 0.1m; y 0 = 0.5m; z 0 = 0.3m (A) Forward feedback signal determined by taking the largest sum of two diametrically opposed OF sensors (horizontal OF sensors, vertical OF sensors) (B) Positioning feedback signal determined by taking the largest output from the OF sensors (right OF sensors; left OF sensors, ventral OF sensors, dorsal OF sensors).

ACKNOWLEDGMENT

We thank R. Brinkworth and D. O'Carroll (Adelaide Uni.) for kindly making their panoramic images available to us, and J. Blanc for improving the English manuscript.

This work was supported partly by CNRS (Life Science;