Julia Chartre

Lukas Rummelhard

Amaury Nègre
email: amaury.negre@gipsa-lab.grenoble-inp.fr

Jean-Alix David

Jérôme Lussereau

Christian Laugier

Situation Awareness for Intelligent Mobility in Dynamic Environments IRT Nanoelec Perfect Platform

Keywords: Perception, Bayesian Filtering, Sensor Fusion, Dynamic Occupancy Grid, Risk Assessment, Intelligent Vehicle, V2X, GPU Computing

In the field of mobile robotics and intelligent vehicles, a key requirement for basic functionalities, such as collision awareness systems, is the ability to perceive and model a relevant representation of the environment where objects are moving. This crucial task is a challenging step in terms of accuracy, complexity and uncertainty management. The objective of this paper is to describe the multi-sensor Bayesian perception approach which has been developed by Inria 1 in the scope IRT Nanoelec 2 , implemented and tested on the IRT experimental platform including an equipped Renault Zoe vehicle. These systems can be adapted to off road environments such as agriculture fields or construction area where workers lives are jeopardized by heavy moving machineries. Embedded multi-sensor Bayesian perception algorithms have been developed on several experimental hardware in cooperation with CEA 3 . For the purpose of computing accurate representation of the dynamic environment, these systems have been implemented on our intelligent vehicle, as well as on portable connected perception units. Shared perception have been implemented using Vehicle to Vehicle (V2V) and Vehicle to Infrastructure (V2I) communication protocols. Experimental results are presented and discussed in the paper.

INTRODUCTION

In the field of mobile robotics and intelligent vehicles, a key requirement for basic functionalities, such as collision awareness systems, is the ability to perceive and model a relevant representation of the environment where objects are moving. This crucial task is a challenging step in terms of accuracy, complexity and uncertainty management. Perception systems can provide accurate information of complex environments to drivers as well as fully autonomous vehicles. In off road environments and especially on construction areas and agriculture fields, collision risks between moving vehicles and workers are important. Safety can be greatly improved by using Advanced Driver Assistance Systems (ADAS) and/or Inria and CEA are working together on situation awareness systems using the knowledge of Inria on perception algorithms for Advanced Driver Assistance Systems solutions [START_REF] Rummelhard | Probabilistic Grid-based Collision Risk Prediction for Driving Application[END_REF], Rummelhard, 2014, 2015[START_REF] Lussereau | Integration of ADAS algorithm in a Vehicle Prototype in IEEE International Workshop on Advanced Robotics and its Social Impacts[END_REF] and the expertise of CEA on software / hardware integration [START_REF] Rakotovao | Real-Time Power-Efficient Integration of Multi-Sensor Occupancy Grid on Many-Core[END_REF][START_REF] Rakotovao | Integration of Multi-Sensor Fusion of Occupancy Grids into Automotive ECUs[END_REF]. IRT Nanoelec provides funding for technology transfer of research activities conducted by Inria as well as a secured platform that replicates realistic roads and intersections that are used to safely conduct experiments involving vehicles, pedestrians and other lightweight transportation systems.

In this paper we present a framework for the perception of the environment capable of merging data from different types of sensors and different sources of information. This data redundancy increases the reliability and robustness of the perception system compared to mono sensor solutions. In this framework, data are shared amongst distributed systems through the ITS-G5 geonetworking protocol and fused using Bayesian filtering algorithms. First, in section 2, we assess risks involving moving vehicles and workers in farming exploitation and construction areas. Then, in section 3, we present existing approaches and the principles of our perception framework, followed by its implementation on an experimental set up and the associated results.

WORKERS SAFETY

Farming Exploitations

Risks in agriculture fields have been studied and well summarized in Groupama "Attitude Prevention" booklet4 . The report mentions three main categories of collision risks:

-Collision risks with private cars when farm tractors are backing off from the exploitation into regular traffic. In this case, risks come from the poor visibility of the tractor driver and the speed difference between the two vehicles.

-Collision with a pedestrian on exploitations caused by low visibility from field dirt, machinery large blind spots and lack of attention of the driver due to task repetitiveness.

-Collision with another private vehicle of the exploitation caused by the same reasons mentioned in the second case.

The various level of customization and control over vehicles and worker safety equipment suggests different approaches. Cost and power requirements must also be adapted, whether it targets a human or a machine.

Construction Site

Risks in construction site are even more critical and very well detailed in INRS (French research and security Institute) report "risk prevention involving vehicles and machines driven on construction site"5 . Risks are similar to the ones in farming exploitations but have a much higher probability to occur due to vehicle density, the large number of workers, and how often they cross vehicle paths. INRS and similar institutes are recommanding to use driving assistance systems on construction machines, which is why solutions already exist on the market. Existing solutions use ultrasonic sensors, cameras, lidars or radio detection. These solutions are based on a single sensor, each kind having its own drawbacks:

-Ultrasonic solutions have a very short range (a few meters). -Camera images can often be hard to understand due to heavy dust and dirt on site. -Radio solutions imply that workers always wear badges, which might not be possible. -Lidar-based solutions generate too many false alarms due to dust and heavy rain.

These mono sensor solutions have very limited sensor fusion and filtering capabilities compared to multi-sensor systems that we describe in this paper. 3 SENSOR FUSION AND FILTERING

Introduction

In the field of perception of dynamic environments, the most classical approach has been to Detect and Track Moving Objects (DATMO), which leads to complex multiple targets object tracking literature [START_REF] Petrovskaya | Aware-ness of road scene participants for autonomous driving[END_REF], (Formtmann, 1980), [START_REF] Khan | An mcmc-based particle filter for tracking multiple interacting targets[END_REF]. Another common approach is the field of occupancy grids [START_REF] Elfes | Using occupancy grids for mobile robot perception and navigation[END_REF], (Movarec, 1988), which works on spatial occupancy without higher level segmentation. This approach presents significant advantages. The model is by design spatially dense, and properly represents information about free space, which is an important data in mobile robotics and intelligent vehicles. Furthermore the delicate data segmentation and recognition step required in object-based representation can be avoided.

The field of occupancy grid based interpretation of the environment is a developed study domain, overlapping various applications, such as intelligent environment management, automatic autonomous navigation or extended vehicle perception. The aim is to produce a compact, regularly subdivided, probabilistic estimation of the spatial occupancy, without requiring the concept of objects. These approaches have been rarely used in dynamic enrironemnt because it requires to enrich each cells with the estimation of velocities which has been computationnaly too heavy for real time applications [START_REF] Coué | Bayesian occupancy filtering for multitarget tracking: an automotive application[END_REF][START_REF] Thrun | Probabilistic Robotics[END_REF][START_REF] Gindele | Bayesian occupancy grid filter for dynamic environments using prior map knowledge[END_REF][START_REF] Danescu | Modeling and tracking the driving environment with a particle-based occupancy grid[END_REF].

The Inria CHROMA6 and eMotion7 research teams have been working on vehicle perception for more than 15 years. One of the main application domain is embedded perception and decision-making for driver assistance systems or autonomous cars. These research studies have led to the development of a perception framework based on Bayesian approaches, which is briefly described in the sequel.

Occupancy Grids

Our perception algorithms are based on a generic occupancy grid framework, initially developed within the Hybrid Sampling Occupancy Filter approach (HSBOF) (Negre, 2014) and further extended in the Conditional Monte Carlo Dense Occupancy Tracker (CMCDOT) (Rummelhard, 2015). This approach infers dynamics of the scene through a hybrid representation of the environment consisting of static and dynamic occupancy, empty spaces and unknown areas. This differentiation enables the use of state specific models (classic occupancy grids for motionless components and sets of moving particles for dynamic occupancy) as well as proper confidence estimation and management of data-less areas. The approach leads to a compact model that dramatically improves the accuracy of the results and the global efficiency in comparison to previous models. This grid-based approach can be computationally very heavy but highly parallelizable, resulting in real time performances on GPUs as shown in paragraph 4.5.2.

Collision Risk Assessment

Most methods for risk estimation detect and track dynamics objects in the scene. The risk is estimated through a Time to Collision (TTC) approach by projecting objects trajectories into the future. Our approach is based on the grid-based CMCDOT framework (Rummelhard, 2015). The idea is to compute estimations of the position in the near future of every static and dynamic part of the grid, instead of reasoning on objects themselves, as well as the trajectory of the vehicle. These estimations are iteratively computed over short time periods, until a potential collision is detected, to which is then associated a TTC. In each cell, the associated TTC are cumulated over different time periods (1, 2, 3 seconds) to estimate a cell-specific risk of collision. This way we generate risk grids, and global risk aggregates, on which we base our reasoning.

This strategy avoids solving the difficult problem of multi-objects detection and tracking while integrating the totality of the information available. It provides a probabilistic estimation of the risk associated to each part of the scene (Rummelhard, 2014).

Dynamic Object Clustering

So far has been described a method how to track spatial occupancy in the scene without object segmentation, but detection and tracking of moving objects (DATMO) is often required for higher level processing, like pedestrian detection. The CMCDOT introduces an object clustering method with limited additional computation. It is based on a particle filter resampling mechanism that associates a unique ID to all of the dynamic particles of an object.

These objects provide a higher level of representation, with much more condensed information that can be easily communicated between various connected systems: vehicles, embedded perception units or lightweight connected devices. This is a first step towards distributed perception.

EXPERIMENTAL SET UP

The described perception framework has been tested and improved throughout the years on several experimental platforms. In this contribution we present the most recent ones that were designed and tested within the scope of IRT Nanoelec.

For our experiments, we chose to use a single type of sensor : lidars, even though our perception framework can fuse heterogeneous kinds of sensors. Lidars are currently very expensive and fragile because of their internal moving mirrors but these issues will certainly be solved in the near future by the release of solid state lidars from companies like Quanergy8 or Leddartech9 . This new generation of lidars will not have any moving parts, reducing drastically their cost and making them more suitable in vibrating environment like farming or construction machines.

Intelligent Vehicle

The first platform that was designed for our experimentations is a Renault Zoe equipped with sensors and processing power. Such a small city car was chosen because of its low cost and ease of driving compared to heavy machines. The purpose of that intelligent vehicle is to run experimentations on the perception of the environment as well as V2V10 and V2I 11 experiments. The controls of the vehicle will be automated in the near future for active driving assistance system experiments. It has been equipped with a set of sensors: lidars, IMU, GPS and cameras. Our perception algorithms are based on lidars only, cameras have been integrated as well for object recognition and image-based processing algorithms. We chose IBEO Lux lidars which have 4 planar laser beams. Three are placed at the front and one at the rear of the car. They are synchronized using the IBEO fusion box. Fused lidar data are then filtered on a Nvidia Titan X GPU in a PC, in the trunk of the car. An IMU is used to get the displacement of the car. Such a solution is applicable when machines can either be industralized or customized with sensors and processing capabilities. It could be suitable for construction machines or tractors that could be equipped with perception capabilities like the Zoe.

We are also working on connected technologies, for vehicle to vehicle and vehicle to infrastructure use cases. All of our experimental platforms are equipped with communication capabilities. The communication protocol that was chosen is the ITS-G5 protocol, the current European standard for V2X communications (V2V and V2I communications). It is a geonetworking peer-to-peer solution that does not require any infrastructure. We chose the implementation made by ITRI (Industrial Technology Research Institute) 12 because of its SDK, 12 https://www.itri.org.tw Linux toolchain and it was the most open and customizable solution.

Object detected by the Zoe from our perception algorithms can be sent through the ITRI unit to other connected devices including other vehicles or infrastructures as explained in distributed perception use cases described in paragraph 4.4.

Connected Devices

Using this ITS-G5 geonetworking communication protocol, other lightweight connected devices has been implemented, like a connected traffic cone. It reads its position from a GPS and broadcasts it to neighboring connected vehicles. It can be placed next to dangerous areas, warning vehicles of a risk ahead of time, when danger may be hidden by other vehicles in front of the car or by a curve. This architecture can be used for any other connected devices like worker hats and small V2X boxes that could be mounted in any vehicles. These solutions would be especially useful when it is impossible to integrate sensors or processing units into existing vehicles. In order for such distributed systems to work, all of the connected vehicles and devices must be synchronized. In our case, this is done through the Network time protocol (NTP), the current time and PPS (1 pulse per second) signal are read from Garmin 18x LVC GPSs.

Perception Units

The third type of equipment that are developed are perception units. The idea is to design portable units, capable of perceiving the environment the same way our intelligent vehicle can, through a set of sensors and a processing unit running our Bayesian filters. Like all of our other connected platforms, they are able to communicate information through ITRI units. These units could be placed in the field or on infrastructures, and communicate detected objects to neighboring vehicles and connected devices. Intelligent vehicles capable of perceiving the environment could receive information from perception units and increase their field of perception, as described in paragraph 4.4.3.

Perception units are equipped with a GPS to localize the unit, a lidar to detect objects, an embedded GPU and a V2X communication device. We chose a Quanergy M8 lidar that generates point clouds of the whole environment for a reasonable price.

Table 2: Sensors integrated on a perception unit.

Improving Off Road Safety

Intelligent Machines

From these experimental platforms, modular solutions can be designed to solve the use cases described in paragraph 2. Off road machines could be equipped with sensors and processing power the same way our intelligent Zoe was. For example, on construction sites, as it is shown in Figure 11, moving machines could detect workers and other vehicles through their own sensors and embedded processing power.

Hardware

Perception Unit Likewise, tractors could be able to detect incoming cars when backing off in traffic as shown in Figure 12.

GPS

V2X Solutions

Another way of estimating collision risk on construction site would be to use connected devices that would signal one another locations as shown in Figure 13. Small low energy connected devices would be placed in worker hats and in vehicles, connected traffic cones would be used to signal dangerous areas and machines.

Distributed Perception

Sharing or distributing the perception of the environment is possible using the experimental platforms described in paragraph 4. It is possible for an intelligent vehicle to share detected objects with other vehicles and for perception units to warn connected vehicles of hidden risks.

For example, two connected vehicles on a construction site could communicate the location of workers to each other reducing the size of their respective blind spots. Figure 14 shows an example of a V2V communication where the vehicle on the right sends the location of worker B to the vehicle on the left, preventing it from backing into the worker. Perception units could be useful in many use cases, for example when tractors are backing off their exploitations. Farmers would place the perception unit on dangerous crossing, the unit would perceive incoming cars from its sensors and warn the driver of the tractor of incoming cars.

P a g e | 8

Another use case is the distribution of the perception between an intelligent vehicle and a perception unit. In dangerous crossing, some of the environement might be hidden by signes or bushes, the perception units could be placed at stategic locations to perceive hidden areas are send object locations to drivers. For example, in Figure 16, worker A is hidden behind the building and is walking towards the digger. The perception unit in orange sees the worker through its sensors and sends its location to the connected digger. Even though the digger is equipped with perception capabilities, the worker is not in its field of perception, preventing the digger from seeing him.

Figure 16: Shared perception between a perception unit and a digger. The field of perception of the unit is in orange, the one of the digger in blue. The worker is seen by the perception units that sends its location to the connected digger.

Our solutions can solve many more use cases, these are some examples of how the environment can be perceived depending on the specificities of each situation and the ability to equip vehicles with sensors.

CMCDOT Software Implementation

ROS

We use the Robot Operating System (ROS) framework 13 that lets us prototype solutions quickly, write reusable code and take advantage of the large community of academic users. A variety of ROS drivers for many kind of sensors are available, it provides many useful tools for code troubleshooting and for graphical visualization of data. We use ROS to encapsulate CMCDOT algorithms, manage and synchronize incoming data from various set of sensors (lidars, cameras, IMU, GPS).

13 www.ros.org

CUDA

CMCDOT algorithms have been highly parallelized in Cuda kernels and optimized to run on multiple Nvidia GPUs: Titan X and embedded Jetson TX1. Our algorithms run in real time, in our experiments they are more than able to fuse and merge lidar data outputs, which are at 20Hz (50 ms per frame) as shown in Table 3.

Table 3: Computational time of sensor data fusion in the grid and CMCDOT algorithms.

Close-to-Market Solutions

Within the IRT Nanoelec we are working closely with our partners at CEA in order to bring our solutions closer to the actual market requirements. CEA conducts efforts to implement our perception framework on low cost, low power architectures like embedded many core. They work on optimized integer implementations of the data fusion that better fits these hardware architectures [START_REF] Rakotovao | Real-Time Power-Efficient Integration of Multi-Sensor Occupancy Grid on Many-Core[END_REF], (Rakotovao, 2106).

EXPERIMENTATIONS

IRT Nanoelec provides a platform, the PTL 14 where experiments are conducted in a secured and dedicated environment that replicates actual roads and infrastructure. Experiments are conducted on two use cases:

-Risk assessment from a Renault Zoe -Distributed perception between the Renault Zoe and a perception unit

Risk Assessment a Renault Zoe

A crash test dummy was designed by Laurence Boissieux from the Inria Grenoble SED team, the dummy simulates a pedestrian crossing a street.

Figure 17 shows our collision risk avoidance experiments on the PTL site with the Renault Zoe and the dummy. We are able to merge and filter data from four IBEO lidars into a single occupancy grid as shown in Figure 18. White dots are lidar impacts merged into a single occupancy grid. Black areas correspond to empty spaces and the gray one to hidden areas, where it is impossible to know if there are objects or not. The red box represents an object with a high risk of collision. This risk was computed based on the Zoe trajectory, cell occupancy and velocity.

The collision risks are shown on the upped image of Figure 19 as red circles. As expected, the dummy was detected as a collision risk for the car. Red arrows indicates the direction of the target. The lower graph corresponds to the probability risks at 3 time horizons, 1 (red), 2 (yellow) and 3 (green) seconds. Figure 19: On the upper camera image, red circles correspond to collision risks and arrows the direction of moving targets. On the lower image, graphs correspond to collision risk probabilities for 3 time horizons: 1s (red), 2s (yellow) and 3s (green).

Distributed perception between a Renault Zoe and a Perception Unit

In these experiments, we want to demonstrate how the field of perception of a vehicle can be enlarged with perception units. The map at the lower right of Figure 20 shows a top down view of the PTL, our experimental site. The Renault Zoe has been parked next to the green building so that most of the environment on the left of the car was hidden. The field of perception of the Zoe is the area shown in blue, we placed a perception unit (red dot) at the lower left angle of the green building. Its field of perception is represented in red. A, B and C yellow crosses represent pedestrians moving in front of the vehicle. Pedestrian A is invisible to the driver because of the green building.

Figure 20 shows the front camera of the vehicle along with lidar impacts on the environment. The field of perception of combined lidars is greater than in the camera, which explains the lack of camera images on each side (black rectangles instead). Red dots correspond to impacts from the lidar on the left of the car, the green ones from the lidar centered at the front of the car and the blue ones from the lidar on the right. Pedestrian B and C have lidar impacts (on their feet) since they are in the field of perception of the vehicle. The perception unit is able to send location of pedestrians A and B through V2X communication protocol. These objects are represented by the yellow dashed boxes. Pedestrians A and B and their associated speed have been detected and sent to the car. This shows that we are able to enlarge the field of perception of drivers using stationary perception units located in the field.

CONCLUSIONS

In this paper we presented a detection and tracking framework for the perception of dynamic environments, based on Bayesian occupancy grids.

We briefly introduced the HSBOF and CMCDOT frameworks that can filter data temporally and spatially from heterogenous kinds of sensors. From the grid representation of the environment where cells states can be static, dynamic, free or unknow, we were able to estimate spatial occupancy, velocity and then cluster objects.

The experimental platforms, including an intelligent vehicle along with connected devices and perception units has been described. The experimental results on pedestrian detection with the Renault Zoe, show that we are able to robustly detect collision risks in real time using highly parallelized algorithms on GPUs. It has also been shown that collision detection algorithms can be run on embedded GPUs in connected perception units, and that perception between intelligent vehicles and perception units can be distributed using V2X communications.

A couple of frequent collision risks that occure on farming exploitations and construction sites has been described. Multiple solution for each of these off-road use cases has been explained using our perception algorithms and/or V2X communications depending on the type of equipment that can be integrated into existing vehicles and safety devices.

Future work includes the enhancement of object classification. Some work has already been completed using lidar data to compute the location of objects and camera to run deep learning techniques (neural network algorithm) to differentiate pedestrians from vehicles, bikes, trucks, etc. We also started implementing higher level object tracking algorithms to track multiple objects over longer period of time and anticipate their behavior and social interactions.

1

  French National Institute for Research on Computer Science and Control 2 Technological Research Institutes Nanoelec, http://www.irtnanoelec.fr/ 3 French Alternative Energies and Atomic Energy Commission Vehicle to Vehicle (V2V) or Vehicle to Infrastructure (V2I) communications.

Figure 2 :

 2 Figure 2: Collision risk in a busy construction areas.

Figure 1 :

 1 Figure 1: Collision risk when tractors are backing off into traffic

Figure 3 :

 3 Figure 3: Data representation in CMCDOT formulation. The environment is divided into cells, with which static, dynamic, empty and unknown variables are associated. The dynamic part of each cell is represented by a set of particles.This methodology is particularly suitable for heterogeneous sensor data fusion (camera, lidars, radars etc…) both spatially and temporally. The occupancy of each cell over time can be estimated from various sensors data whose specific uncertainty (noise, measurement errors) are taken into consideration. Filtered data of each cell are thus much more robust leading to a very reliable global occupancy of the environment reducing by far the number of false detection.

Figure 4 :

 4 Figure 4: Collision risk estimation over time for a specific cell. The cell position is predicted according to its velocity, along with the moving vehicle. This risk profile is computed for every cell and then used to integrate over time the global risk.

Figure 5 :

 5 Figure 5: Sensor location in the Zoe.

Figure 6 :

 6 Figure 6: System architecture of our intelligent vehicle.

Figure 7 :

 7 Figure 7: Connected traffic cone.

Figure 8 :

 8 Figure 8: Connected Traffic cone architecture.

Figure 10 :

 10 Figure 10: Perception Unit System Architecture.

Figure 9 :

 9 Figure 9: Perception unit hardware components.

Figure 13 :

 13 Figure 13: Example of a construction site equipped with connected devices in worker hats, vehicles and traffic cones. The digger on the right is warned by the presence of

Figure 14 :

 14 Figure 14: Distributed perception between two vehicles able to perceive the environment in front of them. The digger on the right perceives worker B though its sensors and sends the location of this worker to the digger on the left to prevent him from backing into worker B.

Figure 11 :

 11 Figure 11: Workers detection from intelligent construction machines. The sensor field of perception is represented in blue. Worker B would be detected as a collision risk.

Figure 12 :

 12 Figure 12: Vehicle detection from an intelligent tractor. The sensor field of perception is represented in blue. The car would be detected as a collision risk by the processing unit of the tractor.

Figure 15 :

 15 Figure 15: Detection of a private vehicle from a perception unit in orange. The driver of the tractor receives a warning from the unit.

Figure 17 :

 17 Figure 17: Risk collision experiments with our dummy crossing a street. © Inria / Photo H. Raguet.

Figure 18 :

 18 Figure 18: Filtered occupancy grid constructed from lidar data.

Figure 20 :

 20 Figure 20: Pedestrian detection in a distributed perception environment. Pedestrian A and B are detected by the perception unit, their locations are sent to the car.

Table 1 :

 1 MCG 2016 -5 th International Conference on Machine Control & Guidance "Facing complex outdoor challenges by inter-disciplinary research" Vichy, France, October 5-6 th , 2016 Sensors types integrated in the Zoe

P a g e | 5

  MCG 2016 -5 th International Conference on Machine Control & Guidance "Facing complex outdoor challenges by inter-disciplinary research" Vichy, France, October 5-6 th , 2016

	Nvidia	Grid fusion	CMCDOT
	Boards	in ms	processing time in
			ms
	Titan X	0.3	7
	Jetson TX1	0.7
	14 Technology Liaison Platform

P a g e | 9

https://www.groupama.fr/assuranceagricole/exploitation/prevention-agricole.html

http://www.inrs.fr/media.html?refINRS=R%20434

https://team.inria.fr/chroma/en/

https://team.inria.fr/e-motion/en/

http://www.quanergy.com/

http://leddartech.com/

Vehicle to Vehicle

ACKNOLEDGEMENTS

This work is part of the IRT Nanoelec Program which benefits from French State aid under the "Programme Investissements d'Avenir" bearing reference ANR-10-AIRT-05.

The authors would like to thank Olivier Simonin for his insightful comments and discussions, as well as Nicolas Turro, Jean-Francois Cuniberto and Laurence Boissieux for their collaboration and continuous support on the experimental platforms.

P a g e | 12

Groupama, farming risks prevention, https://www.groupama.fr/assuranceagricole/exploitation/prevention-agricole.html INRS risk prevention involving by vehicles and machines driven on construction site http://www.inrs.fr/media.html?refINRS=R%20434 Robot Operating Sytem (ROS), http://www.ros.org/


François Alouges
email: francois.alouges@polytechnique.edu

Antonio Desimone
email: desimone@sissa.it

Laetitia Giraldi
email: laetitia.giraldi@inria.fr

Marta Zoppello
email: marta.zoppello@unitn.it

Purcell magneto-elastic swimmer controlled by an external magnetic field

Keywords: Application of nonlinear analysis and design, Tracking, Control in system biology

teaching and research institutions in France or abroad, or from public or private research centers.

Purcell magneto-elastic swimmer controlled by an external magnetic field F. Alouges * A. DeSimone * * L. Giraldi * * * M. Zoppello * * * *

INTRODUCTION

In the last decade, Micro-motility has become a subject of growing interest, both for the biological understanding of micro-organisms and technological applications. In the latter direction, the topic addresses several challenges as for instance the conception of artificial self-propelled and/or easily controllable microscopic robots. Such kind of devices could revolutionize the biomedical applications [START_REF] Peyer | Bioinspired magnetic swimming microrobots for biomedical applications[END_REF] as for instance it could be useful to minimize invasive microsurgical operations [START_REF] Nelson | Microrobots for minimally invasive medicine[END_REF].

One of the few possibilities recently studied in the literature, is to consider an artificial swimmer that possesses a magnetic flexible tail, and use an external magnetic field to act on this flagellum [START_REF] Dreyfus | Microscopic artificial swimmers[END_REF]; [START_REF] Gao | Magnetically powered flexible metal nanowire motors[END_REF][START_REF] Gao | Cargotowing fuel-free magnetic nanoswimmers for targeted drug delivery[END_REF]; [START_REF] Pak | Highspeed propulsion of flexible nanowire motors: Theory and experiments[END_REF]. This particular design is influenced by the locomotion of spermatozoa which achieves their propulsion by propagating travelling wave along their flagellum.

On the other hand, there exists now a quite wide literature that makes a connection between problem of swimming at the micro-scale, and the mathematical control theory. Starting from the pioneering work of Shapere and Wilczek [START_REF] Shapere | Efficiencies of selfpropulsion at low Reynolds number[END_REF], and Montgomery [START_REF] Montgomery | A tour of subriemannian geometries, theirs geodesics and applications[END_REF], the dynamics of self-propelled microscopic artificial swimmers has been considered for instance in Alouges et al. (2013[START_REF] Alouges | Optimal strokes for low Reynolds number swimmers : an example[END_REF]; Alouges and Giraldi (2013); [START_REF] Gérard-Varet | Rough wall effect on microswimmer[END_REF]) where the rate of shape changes of the swimmer is considered as a natural control.

The aim of the present paper is therefore to study thoroughly the mathematics of a magneto-elastic artificial swimmer, as the ones proposed in [START_REF] Alouges | Can magnetic multilayers propel microswimmers mimicking sperm cells?[END_REF], when controlled with an external oscillating magnetic field.

As in [START_REF] Gutman | Simple model of a planar undulating magnetic microswimmer[END_REF], we provide an estimate of the displacement of the swimmer with respect to small amplitude of the fields, generalizing the approach of E. Gutman and Y. Or, to the case of a 3-link swimmer.1

In our model the links are supposed to be uniformly magnetized and linked together with rotational springs. When an external magnetic field is applied its shape is expected to experience a deformation which hopefully leads to a global displacement. It was shown in [START_REF] Alouges | Can magnetic multilayers propel microswimmers mimicking sperm cells?[END_REF] that the dynamics of the a magneto-elastic swimmer is governed by a system of ODEs which is an affine control system with respect to the magnetic field. We adopt the same approach of [START_REF] Alouges | Can magnetic multilayers propel microswimmers mimicking sperm cells?[END_REF] to get the equations of motion. At the micro-scale, the flow is characterized by small Reynolds number. Thus, we assume that the surrounding fluid is governed by Stokes equations which implies that hydrodynamic forces and torques are linear with respect to the velocity distribution on the boundary (rates of deformation and displacement). In that case, the Resistive Force Theory (RFT) (see [START_REF] Gray | The propulsion of seaurchin spermatozoa[END_REF]) provides a simple and concise way to compute a local approximation of hydrodynamic forces involved in our system. The magnetic behavior of the segments is modeled by assuming that their magnetization is always parallel to the segment with fixed magnitude and stray fields, especially magnetic interactions between different segments, are neglected. Only the magnetic torque induced by the external magnetic field on each segment is considered.

After briefly recalling the equation of motion in Section 2, we provide, in Section 3 an estimate of the displacement of the swimmer for an oscillating magnetic field of small amplitude, which steers it along one direction. The last section 4 of the paper is devoted to discuss the main obstruction to get classical controllability results.

MODELING

The model of the magneto-elastic N -link swimmer was already introduced in [START_REF] Alouges | Can magnetic multilayers propel microswimmers mimicking sperm cells?[END_REF]. The two link one was studied in [START_REF] Gutman | Simple model of a planar undulating magnetic microswimmer[END_REF]; [START_REF] Giraldi | Local controllability of the two-link magneto-elastic micro-swimmer[END_REF].Here, we recall briefly the equation of motion focusing on the case N = 3.

We consider a magneto-elastic 3-link Purcell swimmer moving in a plane. This two-dimensional setting is suitable for the study of slender, essentially one-dimensional swimmers exploring planar trajectories as explained in Alouges et al. (2013[START_REF] Alouges | Can magnetic multilayers propel microswimmers mimicking sperm cells?[END_REF]. The swimmer consists of 3 rigid segments, each of length L with articulated joints at their ends (see Fig. 1), moving in the horizontal 2d-plane of the lab-frame. Because of the symmetric geometry of the swimmer, we use slightly different notation and variables than in [START_REF] Alouges | Can magnetic multilayers propel microswimmers mimicking sperm cells?[END_REF]. Indeed we call x = (x, y) the coordinates of the central point of the second segment, θ the angle that it forms with the x-axis, α 2 the relative angle between the first and second segments and finally α 3 the relative between the third and the second segments (see Fig 1). Therefore the position and the orientation of the swimmer are characterized by the triplet (x, y, θ), while its shape is given by (α 2 , α 3 ).

As in [START_REF] Alouges | Can magnetic multilayers propel microswimmers mimicking sperm cells?[END_REF] the three segments are uniformly magnetized and linked together with torsional springs, with elastic constant K, that tend to align the segments one with another. Those produce torques when the segments are not fully aligned. In what follows we assume that H(t) :=

H x (t) H y (t) 0
is horizontal in such a way that the motion holds in the horizontal plane.

• • • • • A 1 A 2 A 3 A 4 x α 2 α 3 θ e x e y

H(t)

Fig. 1. The magneto-elastic Purcell swimmer of shape (α 2 , α 3 ) at the position x in the plane subject to an external magnetic fields H.

Equations of motion

As it was noticed in [START_REF] Alouges | Can magnetic multilayers propel microswimmers mimicking sperm cells?[END_REF] the equations which govern the dynamics of the swimmer form a system of ODEs, which is affine with respect to the magnetic field H(t). The hydrodynamic forces acting on the i-th link, are approximated using RTF, with parallel (resp. perpendicular) drag coefficients ξ i (resp. η i ). In our particular case this system describes the evolution of the position and the shape variables and thus consists of five equations. Those are obtained by writing the balance of forces for the whole swimmer and the balance of torques for the subsystems consisting of the three, two and one rightmost segments.

We call those subsystems S 1 , S 2 , S 3 respectively (S 1 is therefore the whole system). The motion of the swimmer holds in the horizontal plane since only horizontal forces and vertical torques apply. The final system reads:

           F h = 0 , e z • T A1 h + T A1 e + T A1 m = 0 , e z • T A2 h + T A2 e + T A2 m = 0 , e z • T A3 h + T A3 e + T A3 m = 0 .
(1)

Here, F h denotes the total hydrodynamic force acting on the swimmer, T Ai h (resp. T Ai e and T Ai m ) is the hydrodynamic (resp. elastic and magnetic) torque with respect to A i (see Fig. 1), acting on the subsystem S i .

Following the construction made in [START_REF] Alouges | Can magnetic multilayers propel microswimmers mimicking sperm cells?[END_REF] system (1) becomes

M h (θ, α 2 , α 3 )    ẋ θ α2 α3    = -K      0 0 0 α 2 α 3      -M x m (θ, α 2 , α 3 ) H x (t) -M y m (θ, α 2 , α 3 ) H y (t) , (2)
with M h 5 × 5 matrix, M x m and M y m vectors in R 5 all depending on (θ, α 2 , α 3 ). All these matrices can be computed explicitly following the approach given in [START_REF] Alouges | Can magnetic multilayers propel microswimmers mimicking sperm cells?[END_REF]. Finally system (2) can be rewritten as

   ẋ θ α2 α3    =f 0 (θ, α 2 , α 3 ) + f x (θ, α 2 , α 3 )H x (t) + f y (θ, α 2 , α 3 )H y (t) (3)
where,

f 0 = -M -1 h K      0 0 0 α 2 α 3      f j = -M -1 h M j m , j = x, y. (4)
Notice however that the dynamics of (θ, α 2 , α 3 ) is independent of x and can be decoupled. Indeed by block decomposing the matrix M h as

M h = A h B h B T h C h (5)
(A h , B h and C h being respectively 2 × 2, 2 × 3 and 3 × 3 matrices), and considering the first two rows of the system (2), we can solve for ẋ as

ẋ = -A -1 h (θ, α 2 , α 3 ) B h (θ, α 2 , α 3 ) θ α2 α3 = G 1 (θ, α 2 , α 3 ) θ + G 2 (θ, α 2 , α 3 ) α2 + G 3 (θ, α 2 , α 3 ) α3 . (6)
Moreover, the subsystem associated with the shape and the orientation of the swimmer becomes

(-B T h A -1 h B h + C h )   θ α2 α3   = -K 0 α 2 α 3 -Mx m H x (t) -My m H y (t) , (7)
that we rewrite inverting the left hand side matrix as

  θ α2 α3   = g 0 (θ, α 2 , α 3 )+ + g x (θ, α 2 , α 3 )H x (t) + g y (θ, α 2 , α N )H y (t) . (8)
The whole dynamical system (3) is an affine control system with drift where the controls are the two components of the magnetic field. The explicit expression of the dynamics are formally computed by using a symbolic computation software as Mathematica.

STEERING ALONG ONE DIRECTION WITH SMALL SINUSOIDAL MAGNETIC FIELDS

As it was noticed by experiments in [START_REF] Poper | On the dynamics of magnetically driven elastic filaments[END_REF], in the rest we show, that a swimmer with an initial symmetric shape, remains symmetric. Moreover by using an asymptotic expansion we provide an estimate of the swimmer displacement for a prescribed small sinusoidal magnetic field.

A symmetry obstruction

In this part, we assume that the drag coefficients of the 3 links are identical and that they have all a uniform magnetization,

ξ i = ξ and η i = η , M i = M ∀i ∈ {1, 2, 3}
. These assumptions are suitable to describe the behavior of a magnetic filament. We consider the set of symmetric shapes S := (α 2 , α 3 ) ∈ [0, 2π] 2 , α 2 = α 3 , (see Fig. 2). As it was experimentally observed in [START_REF] Poper | On the dynamics of magnetically driven elastic filaments[END_REF], the following Proposition 1 proves that for a swimmer which possesses an initial shape belonging to S, its shape remains symmetric, regardless of the applied magnetic field. Proposition 1. Let T > 0, if at the initial time α 2 (0) = α 3 (0) then for any magnetic field t → (H x (t), H y (t)), applied to the system, the shape of the swimmer remains symmetric i.e.,

α 2 (t) = α 3 (t) , ∀t ∈ [0, T ] .
Proof. This result is based on a symmetry argument. The dynamics must be invariant by the rotation R of angle π about the center O of the lab-frame of the whole system (swimmer and magnetic field). Notice that this rotation changes x to -x, leaves θ invariant, interchanges α 2 and α 3 , and reverses the magnetic field H and the magnetization along the swimmer. Therefore, if the function t → (x(t), θ(t), α 2 (t), α 3 (t)) is solution of the system (3) magnetized in one direction (say +), for an external magnetic field H(t), then the trajectory (-x(t), θ(t), α 3 (t), α 2 (t)) is the solution corresponding to the magnetic field -H(t) of the system magnetized in the opposite direction (say -).

We summarize this by saying that

(x(t), θ(t), α 2 (t), α 3 (t), H(t), +) R -→ (-x(t), θ(t), α 3 (t), α 2 (t), -H(t), -) (9)
where the last component corresponds to the direction of the magnetization along the swimmer.

Similarly, we consider a second transformation T which reverses only the magnetization and the external magnetic field. We remark that the equations of motion (3) depend only on the product M e k ∧ H(t), if the function t → (x(t), θ(t), α 2 (t), α 3 (t)) is solution of (3) with a prescribed magnetic field and magnetization, then it remains a solution with opposite magnetic field and magnetization, so that

(x(t), θ(t), α 2 (t), α 3 (t), H(t), +) T -→ (x(t), θ(t), α 2 (t), α 3 (t), -H(t), -) . (10)
The geometric transformations R and T are sketched in Fig. 2. Checking more formally those symmetry properties of the systems can of course be done on (2) but is left to the reader. Now, composing R and T we have, using the notation above (x(t), θ(t), α 2 (t), α 3 (t), H(t), +)

T •R ---→ (-x(t), θ(t), α 2 (t), α 3 (t), H(t), +) , (11)
which means that by uniqueness of the solution of (3), a swimmer starting at position x(0) = 0 with a symmetric shape (α 2 (0) = α 3 (0)), under whatever driving magnetic field H(t) verifies x(t) = -x(t) , and α 2 (t) = α 3 (t) . It hence experiences no displacement and stays symmetric.

Small oscillating magnetic field

In this subsection we focus on a swimmer satisfying

η 2 = η 3 = η ξ 2 = ξ 3 = ξ with η 1 = η ξ 1 = ξ (12)
These latter assumptions allow to overcome the previous symmetry obstruction. They are suitable in the case of a swimmer with an head.

Starting from a swimmer with a horizontal shape (θ = α 2 = α 3 = 0), we can use the horizontal component of the external magnetic field as a "stabilizer" whereas the oscillating vertical component produces the shape deformation and the motion. In order to understand further what happens when such a field is applied we make the following perturbation analysis. We assume that (H x (t), H y (t)) = (1, sin(ωt)) , (13) and compute the asymptotic expansion of the swimmer displacement with respect to small after a period 2π w . Linearizing the system of equations ( 8) for small angles (θ, α 2 , α 3 ), to first order in , i.e.,

θ α 2 α 3 =   θ α2 α3   + o( ) we get that the triplet   θ α2 α3   satisfies the equation,   θ α2 α3   = A   θ α2 α3   + b sin(ωt) (14)
with A = ∇ (g 0 + g x ) (0, 0, 0) , b = g y (0, 0, 0) . (15) Here, A is the 3 × 3 matrix which depends on the drag coefficients, (η 1 , η) and (ξ 1 , ξ), on the magnetization M and on the elastic constant K. We find that its explicit expression is given by

A = δ a 1 1 a 1 2 a 1 3 a 2 1 a 2 2 a 2 3 a 3 1 a 3 2 a 3 3 where a 1 1 = M η 1 (5η + η 1 ) a 1 2 = (19K + 9M )ηη 1 + 2Kη 2 1 a 1 3 = 2(8K + 3M )ηη 1 + (5K + 3M )η 2 1 a 2 1 = -M 4η 2 + 13η 1 η + η 2 1 a 2 2 = -4(K + M )η 2 -(42K + 23M )η 1 η -2Kη 2 1 a 2 3 = -(28K + 9M )ηη 1 -(5K + 3M ) -η 2 1 a 3 1 = -6M (2ηη 1 + η 2 1 ) a 3 2 = -4(7K + 3M )ηη 1 -5Kη 2 1 a 3 3 = -16(2K + M )ηη 1 -(16K + 11M )η 2 1 and δ = 6 L 3 ηη1(8η+7η1) .
The solution of the system ( 14) is thus given by

θ(t) α2 (t) α3 (t) = 1 2i A + (ω) exp (iωt)b -A -(ω) exp (-iωt)b + c(ω) exp (At) (16)
where A ± (ω) := (-A ± iωI) -1 and c(ω

) = A -(ω)b - A + (ω)b.
The first part of the solution corresponds to a periodic solution, while the last is an exponentially decaying perturbation, as we shall see now. Indeed, by applying Routh-Hurwitz criterion on the characteristic polynomial of A, we prove that the real part of its eigenvalues are all negative. This provides the stability of the asymptotic periodic solution

Let us recall the Routh-Hurwitz criterion.

Lemma 1. (Routh-Hurwitz criterion, [START_REF] Gantmacher | Applications of the Theory of Matrices[END_REF]). Let P (s) = a 3 s 3 + a 2 s 2 + a 1 s + a 0 = 0 a polynomial of third degree. If the two following conditions are satisfied

(1) All the coefficients have the same sign, (2) a 2 a 1 > a 3 a 0 , then the real part of the roots of P are strictly negative.

Here the characteristic polynomial of A reads

p A (λ) = det(A -λId) = a 3 λ 3 + a 2 λ 2 + a 1 λ + a 0 , where
a 0 = detA = - 432M 3K 2 + 4KM + M 2 (2η + η 1 ) d < 0 , a 1 = - 36 d (M 2 10η 2 + 28ηη 1 + η 2 1 + K 2 16η 2 + 64ηη 1 + η 2 1 + KM 31η 2 + 98ηη 1 + 3η 2 1 ) < 0 , a 2 = trA = -12 d M 2η 2 + 17ηη 1 + 5η 2 1 + K 2η 2 + 37ηη 1 + 9η 2 1 < 0 , a 3 = -1 < 0 d = L 9 η 2 8ηη 1 + 7η 2 1 so that condition 1. of Lemma 1 is satisfied. Moreover, a 3 a 0 -a 2 a 1 = -1 L 9 η 3 η 2 1 (8η + 7η 1 ) 2 432 2(K + M )(16K 3 + 31KM + 10M 2 )η 4 + + 5 144K 3 + 339K 2 M + 217KM 2 + 42M 3 η 3 η 1 + 2514K 3 + 5015K 2 M + 2867KM 2 + 506M 3 η 2 η 2 1 + + 613K 3 + 1309K 2 M + 802KM 2 + 150M 3 ηη 3 1 + (9K + 5M ) K 2 + 3KM + M 2 η 4 1 < 0 . (17)
Therefore, we deduce that the steady state of the equation ( 14) is stable, i.e., exp (At) → 0 , as t → ∞ ,

and the solution of ( 16) exponentially converges to the

periodic solution   θ∞ (t) α∞ 2 (t) α∞ 3 (t)   = 1 2i A + (ω) exp (iωt)b -A -(ω) exp (-iωt)b , (19)
and in particular θ ∼ θ∞ oscillates around 0 indicating that the swimmer stays nearly horizontal, stabilized by the horizontal component of the magnetic field. Similarly, the fact that the shape variables (α 2 , α 3 ) are periodic (and small) indicates that the swimmer stays nearly straight.

In order to go further, and compute the (asymptotic) net displacement of the swimmer after one period of the oscillating external field, we linearize as well the equation ( 6) to first order in (θ, α 2 , α 3 ) near (0, 0, 0).

Noting, ẋ = G x (θ, α 2 , α 3 )•   θ α2 α3   and ẏ = G y (θ, α 2 , α 3 )•   θ α2 α3   where G x (resp. G y ) is the 1 × 3 matrix composed of (G i • e x ) i=1,•••3 (resp. (G i • e y ) i=1,•••3 ), we obtain ∆x = 2π ω 0 G x (0, 0, 0) + t θ, α2 , α3 ∇G x (0, 0, 0) θ∞ α∞ 2 α∞ 3 dt +o( 2 ) . (20) Since, t →   θ∞ (t) α∞ 2 (t) α∞ 3 (t)   is periodic, the latter equality reads ∆x = 2 2π ω 0 t θ∞ , α∞ 2 , α∞ 3 ∇G x (0, 0, 0)   θ∞ α∞ 2 α∞ 3   dt +o( 2 ) , (21)
and a straight forward computation leads to express

∇G x (0, 0, 0) = L 2 (2η + η 1 )     2(η -η 1 ) -η 1 η - (6ηη 1 -4ηξ + η 1 ξ 1 (2ξ + ξ 1 ) - η 1 (2η + ξ 1 ) (2ξ + ξ 1 ) - η(η 1 -ξ 1 ) (2ξ + ξ 1 ) (2η 2 + 4ηη 1 -3η 1 ξ) (2ξ + ξ 1 ) η 1 (η -ξ) (2ξ + ξ 1 ) η(η + η 1 + ξ) (2ξ + ξ 1 )     .
Similarly, the same formula holds for ∆y by substituting G x for G y and in this case, ∇G y (0, 0, 0) = 0 0 0 0 0 0 0 0 0 , thus, ∆y = o( 2 ). It follows that the leading term, with respect to small angles, of the trajectory of the swimmer along the y-axis is negligible after one period of the oscillating fields compare to the one along the x-axis.

From now on, we focus on the x-displacement of the swimmer, ∆x and we prove that the leading term of order 2 does not vanish after one period of the oscillating fields.

Plugging ( 19) into (21) and noting that the two terms vanish because of periodicity,

2π ω 0 t A + (ω) exp (iωt)b ∇G x (0, 0, 0) A + (ω) exp (iωt)b = 0 , 2π ω 0 t A -(ω) exp (-iωt)b ∇G x (0, 0, 0) A -(ω) exp (-iωt)b = 0 we obtain ∆x = 2 2π ω 0 ω 4 t b t A + (ω) ∇G x (0, 0, 0) -t ∇G x (0, 0, 0) A -(ω)b . (22)
Notice already that since A + (0) = A -(0) = -A -1 , ∆x tends to 0 when ω tends to 0. A very low frequency produces no net motion (at order 2 ), even after one period.

Moreover, the 3×3 matrix (∇G x (0, 0, 0 -t ∇G x (0, 0, 0)) is skew-symmetric and not null. Therefore, 0 is an eigenvalue of multiplicity 1. Let us denote by u its associated eigenvector. A direct computation, still using Mathematica, leads to

u =   η 1 ξ + ηξ 1 -2ηη 1 2η 2 + 4η 1 η -2ξη -ξ 1 η -3η 1 ξ 6ηη 1 -2ξη 1 -4ηξ 1   .
Thus, to ensure that ( 22) is not null, it is sufficient to prove that the three of vectors {u, A -(ω)b, A + (ω)b} are independent. But, for large frequencies ω, we can expand the matrix A ± (ω) as

A ± (ω) = ± I iω - A ω 2 + o( 1 ω 2 ) . (23)
and

det(u, A -(ω)b, A + (ω)b) = det(u, b ω , Ab ω 2 ) + o( 1 w 4 ) = -L 9 ω 3 η 3 η 2 1 (8η + 7η 1 ) 2 (3ηη 1 -2ηξ 1 -η 1 ξ) -1 × ηη 1 (113K + 38M ) + 216M 2 (2η + η 1 ) η 2 ξ 1 4η 2 (5K + 2M ) + η 2 1 (29K + 8M ) + ηη 1 (η 1 -η) K 4η 2 + 37ηη 1 + 13η 2 1 + 6η 1 M (2η + η 1 ) -η 1 ξ 8η 3 (2K + M ) + 2η 2 η 1 (40K + 13M ) + ηη 2 1 (53K + 14M ) + η 3 1 (13K + 6M ) .
This determinant does not vanish identically and we then obtain that by prescribing an oscillating magnetic field as ( 13), the magneto-elastic Purcell swimmer moves along the x-axis. Notice that here the assumption on the drag coefficients ( 12) is crucial.

Moreover, ∆x tends to 0 as ω → 0 and ω → ∞. This suggests the existence of an optimal frequency to drive the swimmer as was already observed in [START_REF] Alouges | Can magnetic multilayers propel microswimmers mimicking sperm cells?[END_REF] (see Fig. 5 and7).

DISCUSSION

This section underlines the challenge that we have to face in order to control this magnetic micro-swimmer.

The previous result indicates that with a small sinusoidal magnetic field we are able to control the direction of the swimmer's displacement, but it does not imply neither global or local controllability properties.

The latter local property is classically obtained by verifying the Kalman condition at an equilibrium point.

Thus let us consider the system (3), around (X e , u e ) := ((x, 0, 0, 0), (0, 0)), which is an equilibrium point. At such a point, the swimmer is aligned with the horizontal field and thus, the torque due to the horizontal field vanishes leading to f x (X e ) = 0. Moreover, the system being invariant under translations, f 0 does not depend on (x, y), from which we deduce that the matrix range R of the matrix ∇f 0 (X e ) is at most of dimension 3. The Lie bracket ad f0 (f x )(X e ) = [f 0 , f x ](X e ) vanishes since f 0 (X e ) = 0 and f x (X e ) = 0. By induction, for k > 1 ad k f0 (f x )(X e ) = [f 0 , ad k-1 f0 (f x )](X e ) = 0 . As far as ad k f0 (f y )(X e ) is concerned, we have, still by induction, for all k > 1 ad k f0 (f y )(X e ) = [f 0 , ad k-1 f0 (f y )](X e ) ∈ R . Therefore Span ad k f0 (f x )(X e ), ad k f0 (f y )(X e ), k ≥ 0 = Span f y (X e ), ad k f0 (f y )(X e ), k ≥ 0 ⊂ Span (f y (X e )) + R is at most of dimension 4. Thus the Kalman condition is not satisfied. It turns out that for non horizontal straight swimmers, i.e. X e = (x, θ, 0, 0) with θ = 0, the same situation occurs due to the fact that f x (X e ) = -tan(θ)f y (X e ). ( 24)

Moreover let us notice that [f x , f y ](X e ) = f 3 y,Xe f y (X e )

where f 3 y,Xe is the third component of the vector f y (X e ). Therefore the Lie algebra Lie{f 0 , f x , f y }| Xe span a vector space of dimension 4. By changing the reference frame, a similar argument holds for all equilibrium points such as {(x, θ, 0, 0), θ ∈ [0, 2π]}. This means that also the classical LARC condition is not satisfied and then the Sussmann condition does not hold (see [START_REF] Coron | Control and Nonlinearity[END_REF]).

Here we have underlined the fact that the the magnetoelastic Purcell swimmer model is singular at the straight position which makes hard to get a local controllability result. Of course if instead the swimmer starts at a non straight position, thanks to the boundedness of the drift, the magnetic field can compensate it and drive the swimmer.

CONCLUSION

In this paper, by prescribing a particular oscillating field, we make an asymptotic expansion of the displacement of the swimmer, proving that this particular field allows to steer the swimmer along one direction. Moreover we highlight the difficulties to get controllability result by showing that the classical conditions are not satisfied. It indicates that sophisticated techniques (see for instance [START_REF] Giraldi | Local controllability of the two-link magneto-elastic micro-swimmer[END_REF]) are required to obtain such controllability result.

Fig. 2 .

 2 Fig. 2. two R and T of the system. Remark 3.1. Similar argument holds for a swimmer composed by an odd number of links.

The 3-link swimmer was already proposed in[START_REF] Purcell | Life at low Reynolds number[END_REF]. It is one of the more famous systems, which was further considered in the literature[START_REF] Becker | On self-propulsion of micro-machines at low Reynolds number: Purcell's three-link swimmer[END_REF];[START_REF] Giraldi | Controllability and optimal strokes for N-link microswimmer[END_REF];[START_REF] Passov | Dynamics of Purcell's threelink microswimmer with a passive elastic tail[END_REF];[START_REF] Tam | Optimal strokes patterns for Purcell's three link swimmer[END_REF] for instance.


Marta Zoppello
email: marta.zoppello@unitn.it

Antonio Desimone
email: desimone@sissa.it

François Alouges
email: francois.alouges@polytechnique.edu

Laetitia Giraldi
email: laetitia.giraldi@inria.fr

Pierre Martinon
email: pierre.martinon@inria.fr

Franc ¸ois Alouges
email: fran-cois.alouges@polytechnique.edu

Optimal control of slender microswimmers

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Mathematical setting of the problem

In this section we describe the kinematics of the so called N-link swimmer, inspired by the Purcell's 3-link swimmer. A discrete representation od the swimmer's curvature is provided by he angles between successive inks. These angles are considered as freely prescribed shape parameters. We then write the balance of total viscous force and torque, i.e. the equations of motion, solving for the time evolution of position and orientation of the swimmer in response to a prescribed history of (concentrated) curvatures along the swimmer's body.

Kinematics of the N-link swimmer

Here we are interested in essentially one-dimensional swimmers moving in a plane. This setting is suitable for the study of slender, one-dimensional swimmers exploring planar trajectories. The general case is a bit more involved because of the nonadditivity of three-dimensional rotations, see e.g. [START_REF] Alouges | Optimally swimming stokesian robots[END_REF], but it can be handled with similar techniques.

Our swimmer is composed of N rigid links with joints at their ends (see Fig. 1.1), moving in the plane (e x , e y ) (2d lab-frame). We set e z := e x × e y . The i-th link is the segment with end points x i and x i+1 . It has length L i > 0 and makes an angle θ i with the vector e x . The lenght of the sticks is chosen such that the size of the swimmer is of order of µm. We define by x i := (x i , y i ) (i = 1, • • • , N) the coordinates of the first end of each link. Notice that, for i ∈ {2 • • • N}, x i is a function of x 1 , θ k and L k , with k ∈ {1 • • • i -1}:

x i := x 1 + i-1 ∑ k=1 L k cos(θ k ) sin(θ k ) . (1.1)
The swimmer is described by two kind of variables:

• the state variables which denote the position and the orientation of one selected link, labeled as the i * -th one; • the shape variables which describe the relative angles between successive links.

For each link with i > i * , this is the angle relative to the preceding one, denoted by α i = θ iθ i-1 , for i * < i ≤ N. For i < i * this is the the angle relative to the following one, denoted by α i = θ i+1θ i , for 1 ≤ i < i * .

For example, if the triplet (x 1 , θ 1 ) is the state of the swimmer then the vector

(α 2 = θ 2 -θ 1 , • • • , α N = θ N -θ N-1
) describes the shape of the swimmer. We will use these coordinates in the rest of the chapter.

ex ey

Equations of motion

The dynamics of the swimmer is governed by a system of three ODEs. This system represents the Newton laws, in which inertia is neglected, namely

F = 0 , e z • T x 1 = 0 , (1.2)
where F is the total force that the fluid exerts on the swimmer and T x 1 is the corresponding total torque computed with respect to the point

x 1 .
In what follows we use the local drag approximation of Resistive Force Theory (RTF), to couple the fluid and the swimmer. According to this approximation the hydrodinamic forces are linear in the velocities of each point. More precisely denoting by s the arc length coordinate on the i-th link (0 ≤ s ≤ L i ), by v i (s) the velocity of the corresponding point, and calling e i = cos(θ i ) sin(θ i ) and e ⊥ i = sin(θ i ) cos(θ i ) the unit vectors in the directions parallel and perpendicular to the i-th link respectively, we can write x i (s) = x i + se i . Differentiating, we obtain,

v i (s) = ẋi + s θi e ⊥ i . (1.3)
According to RTF, the density of the force f i acting on the i-th segment depends linearly on the velocity and can be written as

f i (s) := -ξ (v i (s) • e i ) e i -η v i (s) • e ⊥ i e ⊥ i , (1.4)
where ξ and η are the drag coefficients in the directions of e i and e ⊥ i respectively, measured in N s m -2 . We thus obtain

             F = N ∑ i=1 L i 0 f i (s) ds , e z • T x 1 = e z • N ∑ i=1 L i 0 (x i (s) -x 1 ) × f i (s) ds .
(1.5)

Using (1.3) and (1.4) into (1.5), the total force and torque can be expressed as

F = - N ∑ i=1 L i ξ (ẋ i • e i ) e i + L i η(ẋ i • e ⊥ i ) + L 2 i 2 η θi e ⊥ i , (1.6)
and

e z • T x 1 = - N ∑ i=1 L 2 i 2 η ẋi • e ⊥ i + L 3 i 3 η θi + (x i -x 1 ) × L i ξ (ẋ i • e i ) e i + L i η(ẋ i • e ⊥ i ) + L 2 i 2 η θi e ⊥ i • e z . (1.7)
Moreover, the differentiation of (1.1) gives

ẋi = ẋ1 + i-1 ∑ k=1 L k θk e ⊥ k , (1.8)
which is linear in ẋ1 and ( θk ) 1≤k≤N . This implies that also (1.6) and (1.7) are linear in ẋ1 and θi for i ∈ [1 • • • N], and therefore system (1.2) reads

F e z • T x 1 = M (θ 1 , • • • , θ N )        ẋ1 θ1 θ2 . . . θN        =   0 0 0   . (1.9)
We point out that for all i ∈ {2, • • • , N}, α i = θ iθ i-1 , equations (1.6) and (1.7) can be written using the relative angles

(α i ) i=2,•••,N instead of the variables (θ i ) 2≤i≤N recalling that θ i = θ i-1 + α i , i = 2, • • • , N .
(1.10)

To this end, we introduce the matrix C defined by

C =                1 0 • • • • • • • • • • • • 0 0 1 . . . .
0 0 0 • • • 0 -1 1                (1.11)
and obtain

C        ẋ1 θ1 θ2 . . . θN        =        ẋ1 θ1 α2 . . . αN       
.

(1.12) Thus, by setting

N(θ 1 , α 2 , • • • , α N ) := M (θ 1 , θ 2 (θ 1 , α 2 , • • • , α N ), • • • , θ N (θ 1 , α 2 , • • • , α N )) C -1 , (1.13)
system (1.9) can be rewritten in the equivalent form

N (θ 1 , α 2 , • • • , α N )        ẋ1 θ1 α2 . . . αN        =   0 0 0   .
(1.14)

We observe that we can decompose the 3

× (N + 2) matrix N (θ 1 , α 2 , • • • , α N ) in blocks into a 3 × 3 sub-matrix A (θ 1 , α 2 , • • • , α N ) and a 3 × (N -1) sub-matrix B (θ 1 , α 2 , • • • , α N ), according to N = (A | B) . ( 1

.15)

The matrix A is known as the "grand-resistance-matrix" of a rigid system evolving at frozen shape, i.e., with αi ≡ 0, i = 2, . . . , N, see [START_REF] Happel | Low Reynolds number hydrodynamics with special applications to particulate media[END_REF]. it can be easily verified that it is symmetric and negative definite [START_REF] Happel | Low Reynolds number hydrodynamics with special applications to particulate media[END_REF] and hence invertible. Therefore the equations of motion of the swimmer turn out to be affine system without drift. Indeed, solving (1.14) for (ẋ 1 , θ1 ) leads to

ẋ1 θ1 = -A -1 (θ 1 , α 2 , • • • , α N ) B (θ 1 , α 2 , • • • , α N )    α2 . . . αN   
which can be rewritten in the form

ẋ1 θ1 = N ∑ i=2 gi (θ 1 , α 2 , • • • , α N ) αi , (1.16)
where the N -1 vector fields {g i } N i=2 , are the columns of the 3 × (N -1) matrix -A -1 B.

The equation above encode the link between the displacement (both translation and rotation) of the swimmer and its deformation. More precisely, prescribing the shape functions t → (α 2 , • • • , α N )(t), the motion of the swimmer is obtained by solving the system (1.16). In what follows we call stroke a time-periodic shape change, i.e., the functions t → α i (t), i = 2, • • • , N are all periodic, with the same period.

In order to solve (1.16) numerically, we need to compute the vector fields gi explicitly. To this end, we notice that the total force F and the total torqueT x 1 depend linearly on (ẋ i ) 1≤i≤N and ( θi ) 1≤i≤N and that these quantities depend in turn linearly on (ẋ 1 , θ1 , • • • , θN ) in view of (1.8). Therefore, we can rewrite (1.6) and (1.7) as

F = P 1             ẋ1 . . . ẋN -- θ1 . . . θN             = P 1 Q      ẋ1 θ1 . . . θN      , e z • T x 1 = P 2             ẋ1 . . . ẋN -- θ1 . . . θN             = P 2 Q      ẋ1 θ1 . . . θN      , (1.17)
where

P 1 := -m 1 • • • -m N | η 2 L 2 1 e ⊥ 1 • • • η 2 L 2 N e ⊥ N with m i := L i (ξ e i ⊗ e i + ηe i ⊥ ⊗ e i ⊥ ) for i = 1 • • • N, P 2 := • • • -(L 2 i ηe ⊥ i + (x i -x 1 ) × m i ) T • • • | • • • ηL 2 i ( L i 3 + (x i -x 1 )×e ⊥ i •e z 2 ) • • • ,
and, finally,

Q =               1 0 0 0 • • • 0 1 L 1 e ⊥ 1 0 0 • • • 0 1 L 1 e ⊥ 1 L 2 e ⊥ 2 0 • • • 0 . . . . . . . . . . . . • • • 0 1 L 1 e ⊥ 1 L 2 e ⊥ 2 • • • L N-1 e ⊥ N-1 0 0 . . . Id 0               .
We thus have that the matrix M in (1.13) is

M = P 1 Q P 2 q
and we can compute N = C -1 M, where C -1 is explicitly given as

C -1 =                1 0 • • • • • • • • • • • • 0 0 1 . . . . . . . . . . . . . . . 0 0 1 . . . . . . . . . 0 0 1 . . . . . . . . . . . . . . . 1 . . . . . . 0 . . . . . . . . . . . . . . . . . . 0 0 0 1 • • • 1 1 1                . (1.18)
Matrices A and B are obtained from the columns of N as in (1.15) and, finally, the vectors gi are simply the columns of -A -1 B. Finally the dynamics of the swimmer is expressed as

       α2 . . . αN ẋ1 θ1        = N-1 ∑ i=1 b i gi (θ 1 , α 2 , • • • , α N ) αi+1 . (1.19)
where b i is the i-th vector of the canonical basis of R N-1 .

Applications of the N-link swimmer

The N-link swimmer model is very useful and can be used as a discrete approximation of a swimmer's flexible tail whose shape is controlled by curvature. We show in this section, how curvature control can be implemented in our model in a concrete case reproducing the motion of a sperm cell analyzed in [START_REF] Friedrich | High-precision tracking of sperm swimming fine structure provides strong test of resistive force theory[END_REF].

Curvature approximation

Here, we describe the method to approximate the curvature of a beating tail with the discrete N link swimmer model. Let L > 0 be the total length of the flexible tail and let r(s,t) be the position, in the body frame of the swimmer (see Figure 1.2), at time t > 0 of the point of arc-length coordinate s ∈ [0, L] along the tail. We also define the angle between the tangent vector to the tail at the point r(s,t) and the x-axis in the lab-frame as Ψ (s,t). It is wll known that the derivative of Ψ (s,t) with respect to s is the local curvature of the curve. We discretize the swimmer's tail into N equal parts of lenght L i = L/N, and define the angles (θ i ) 1≤i≤N by averaging the function Ψ (s,t) on the interval [iL/N, (i + 1)L/N]

θ i (t) = N L iL N (i-1)L N Ψ (s,t) ds , i = 1 . . . N . (1.20)
Finally, differentiating (1.20) with respect to time we get the angular velocities θi ,

i = 1, • • • , N, θi (t) = N L iL N (i-1)L N ∂Ψ (s,t) ∂t ds , i = 1 . . . N . (1.21) -2 0 2 4 6 8
x 10

N-link approximation of sperm cell swimmer

We now focus on reproducing with our model the motion of a sperm cell and compare to the one reported in [START_REF] Friedrich | High-precision tracking of sperm swimming fine structure provides strong test of resistive force theory[END_REF]. To preform this comparison, we have to to take into account the presence of the head of the sperm cell. In order to do this, we modify the first segment of the N-link swimmer so that it has its own translational and rotational viscous drag. Indeed, we denote by x 1 the position of the central point of the head and θ 1 the angle that the direction of first segment (attached to the head) (e 1 ) makes with the horizontal axis. The movement of the head generates a viscous force and torque that are given by

F head = -ξ head (ẋ 1 • e 1 )e 1 -η head (ẋ 1 • e ⊥ 1 )e ⊥ 1 , (1.22)
and

T head • e z = -ζ head θ1 . (1.23)
We also fix that the length of the head is L head = 10 µm and we assume again that L is the length of the tail which is attached to one of the extremities of the head segment. As suggested in [START_REF] Riedel-Kruse | How molecular motors shape the flagellar beat[END_REF], the wave profile along the tail of the sperm cell swimmer was obtained from experimental data, keeping only the two first Fourier modes, and we use the method described before in section 1.2.1 to approximate the tail's motion. More precisely, we describe the shape of the wave shown in Fig. 1.3 by

r(s,t) = L head 2 e 1 (t) + s 0 cos(Ψ (u,t))e 1 (t) + sin(Ψ (u,t)e ⊥ 1 (t)du . (1.24)
where

Ψ (s,t) = K 0 s + 2A 0 s cos(ωt - 2πs λ ) .
(1.25)

In the previous equations, K 0 represents the mean flagellar curvature while ω, λ and A 0 are the frequency, the wave-length and the amplitude of the wave respectively. Following [START_REF] Friedrich | High-precision tracking of sperm swimming fine structure provides strong test of resistive force theory[END_REF], in the next numerical simulations we use the following values for the wave parameters:

A 0 = 15.2•10 3 rad m -1 , K 0 = 19.1•10 3 rad m -1 , ω = 200 rad , s -1 and λ = 71.6 • 10 -6 m.
Except the first segment, we discretize the rest of the tail with N -1 segments of extremities (x i , x i+1 ) for i = 2, • • • , N. We use the method described in section 1.2.1, to approximate the beating wave and obtain the shapes shown in Figure 1.4 for one period (0 ≤ t ≤ 2π ω ).

-2 0 2 4 6 8

x 10 With the above notation, the equations of motion become

         F = F head + N ∑ i=1 L i 0 f i (s) ds , T x 1 = T head + N ∑ i=1 L i 0 f i (s) × (x i (s) -x 1 ) ds .
(1.26)

where

L i = L/N is the length of each segment (x i , x i+1 ) for i = 2, • • • , N, while the first segment, also of size L 1 = L/N is given by (x 1 + L head 2 e 1 , x 2
). Since the two previous formulas (1.26) are linear in θ1 and ẋ1 , we end up with the same compact expression of the equations of motion (1.16). More precisely, the matrix P 1 and P 2 defined in system (1.17) are replaced by 6 below, summarize our results that are in perfect agreement with those of [START_REF] Friedrich | High-precision tracking of sperm swimming fine structure provides strong test of resistive force theory[END_REF] (see Figure 3 for the trajectory and Figure 4 for the various speeds).

P head 1 := -ξ head e 1 ⊗ e 1 + η head e ⊥ 1 ⊗ e ⊥ 1 -m 1 • • • -m N | η 2 L 2 1 e ⊥ 1 • • • η 2 L 2 N e ⊥ N and P head 2 := -p 1 • • • -p N | -ζ head + q 1 q 2 • • • q N , with m i := L i (ξ e i ⊗ e i + ηe ⊥ i ⊗ e ⊥ i ) for i = 1 • • • N, and
p i := (L 2 i ηe ⊥ i + (x i -x 1 ) × m i ) T , q i := ηL 2 i ( L i 3 + (x i -x 1 )×e ⊥ i •e z 2 ), for i = 1 • • • N.

Controllability

This section is devoted to the controllability of the N-link swimmer, which is its ability to move between two fixed configurations prescribing (controlling) its shape parameters. More precisely we prove that there exist control shape functions which allow the swimmer to move everywhere in the plane.

Theorem 1.1. Consider the N-link swimmer described in section 1.1 evolving in the space R 2 . Then for almost every lengths of the sticks (L i ) i=1,•••,N and for any initial configuration

(x i 1 , θ i 1 , α i 2 , • • • , α i N ) ∈ R 2 × [0, 2π] N , any final configu- ration (x f 1 , θ f 1 , α f 2 , • • • , α f N )
and any final time T > 0, there exists a shape function

(α 2 , • • • , α N ) ∈ W 1,∞ ([0, T ]), satisfying (α 2 , • • • , α N )(0) = (α i 2 , • • • , α i N ) and (α 2 , • • • , α N )(T ) = (α f 2 , • • • , α f N )
and such that if the self-propelled swimmer starts in position (x i 1 , θ i 1 ) with the shape

(α i 2 , • • • , α i N ) at time t = 0, it ends at position (x f 1 , θ f 1 ) and shape (α f 2 , • • • , α f N ) at time t = T by changing its shape along (α 2 , • • • , α N )(t).
Proof. The proof of the theorem is divided into three steps. First of all, we show the analyticity of the dynamics vector fields. Then, we prove the controllability of the 3-link swimmer (Purcell swimmer), exploiting the Chow theorem and the Orbit theorem. Finally, we generalize the result to the case of N links. We start by recalling some classical results used in the proof.

Classical results in geometric control

Theorem 1.2. (Chow (see [START_REF] Coron | Control and Nonlinearity[END_REF])) Let m,n ∈ N and let (g i ) i=1,n be C ∞ vector fields on R n . Consider the control system, of state trajectory q,

q = m ∑ i=1 u i g i (q), (1.27
)

with input function u = (u i ) i=1,m ∈ L ∞ ([0, +∞[, B R n (0, δ )) for some δ > 0.
Let O an open and connected set of R n and assume that Lie q (g 1 , ...g m ) = R n q ∈ O Then the system (1.27) is controllable, i.e., for every q 0 , q 1 in O and for every T > 0 exists u ∈ L ∞ ((0, T ), B R n (0, δ )) such that q(0) = q 0 and q(T ) = q 1 and q(t) ∈ O for every t ∈ [0, T ].

If the vector fields are analytic, we can apply the Orbit Theorem to extend the dimension property of the Lie algebra defined by the dynamics vector fields on the whole orbit.

Theorem 1.3. (Orbit (see [START_REF] Jurdjevic | Geometric control theory[END_REF]) Let M be an analytic manifold, and F a family of analytic vector fields on M . Then a) each orbit of F is an analytic submanifold of M , and b) if N is an orbit of F , then the tangent space of N at x is given by Lie x (F ). In particular the dimension of Lie x (F ) is constant as x varies on N.

In our case, the manifold in which the state and the shape of the swimmer evolve is defined by

M := [0, 2π] N-1 × R 2 × [0, 2π].
The vector fields of the dynamics are denoted by

g i (θ 1 , α 2 , • • • , α N ) := b i gi (θ 1 , α 2 , • • • , α N )
.

We say that the Lie algebra of the family of vector fields

{g i } i=1,•••,N-1 is fully generated at the point q = (α 2 , • • • , α N , x 1 , y 1 , θ 1 ) ∈ M if the tangent space of the manifold, T q M , is equal to the Lie algebra Lie((g i )) i=1,•••,N-1 )(q).

Regularity

We first prove that the vector fields (g i ) are analytic on M . From (1.6) and (1.7), the entries of the matrices A and B are analytic functions on [0, 2π] N . Since the coefficients of A -1 are obtained by multiplication and division of those of A, and because det(A) = 0 (A is symmetric and negative defined), the entries of A -1 remain analytic functions on [0, 2π] N . Thus, the (g i ) i=1,•••,N := A -1 B are analytic on [0, 2π] N .

Controllability of the Purcell Swimmer (N=3)

Setting N = 3 in (1.19) the dynamics becomes       α2 α3 ẋ1 ẏ1 θ1       = g 1 (θ 1 , α 2 , α 3 ) α2 + g 2 (θ 1 , α 2 , α 3 ) α3 .
(1.28)

To prove the controllability of this system we want use Theorem 1.2. Therefore we compute the Lie algebra of the vector fields g 1 and g 2 for any θ 1 ∈ [0, 2π] at (α 2 , α 3 ) = (0, 0), for a swimmer whose sticks have the length L 1 = L 3 = L and L 2 = 2L where L > 0. First we have

g 1 (θ 1 , 0, 0) =           1 0 9L sin(θ 1 ) 64
-9L cos(θ 1 ) 64 27 32           g 2 (θ 1 , 0, 0) =           0 1 -7L sin(θ 1 ) 64 7L cos(θ 1 ) 64 -5 32          
Then, the iterated Lie brackets are equals to

[g 1 , g 2 ](θ 1 , 0, 0) = 0, 0, 7L(η-ξ ) cos(θ 1 ) 128ξ , 7L(η-ξ ) sin(θ 1 ) 128ξ , 0 T , [g 1 , [g 1 , g 2 ]](θ 1 , 0, 0) =              0 0 - L(126η 2 +31ξ η-76ξ 2 )sin(θ1) 4096ηξ L(126η 2 +31ξ η-76ξ 2 )cos(θ1) 4096ηξ - 3(9η 2 -4ξ η+4ξ 2 ) 2048ηξ              , [g 2 , [g 1 , g 2 ]](θ 1 , 0, 0) =              0 0 L(36η 2 -103ξ η+148ξ 2 )sin(θ1) 4096ηξ - L(36η 2 -103ξ η+148ξ 2 )cos(θ1) 4096ηξ 3(9η 2 -4ξ η+4ξ 2 ) 2048ηξ              .
To see if they are linearly independent we compute the determinant of the matrix whose columns are the 5 previous vector fields that is equal to

21L 2 (η -ξ ) 2 (45η + 112ξ ) 9η 2 -4ηξ + 4ξ 2 536870912η 2 ξ 3 . (1.29)
Since the drag coefficients ξ and η are positive, this determinant is null only when ξ = η. This would imply an isotropic drag, as we would have if we use spheres instead of sticks. Thus in our case the Lie algebra of the vector fields g 1 and g 2 is fully generated at the point (θ 1 , 0, 0), for any

θ 1 ∈ [0, 2π]. Notice that any point (α 2 , α 3 , x 1 , θ 1 ) ∈ [0, 2π] 2 × R 2 × [0, 2π]
belongs to the orbit of the point (0, 0, x 1 , θ 1 ). Since the vector fields are analytic, the Orbit Theorem 1.3 guarantees that the Lie algebra of g 1 and g 2 is fully generated everywhere in the manifold

M = [0, 2π] 2 × R 2 × [0, 2π].
To conclude, by Chow Theorem (1.2) we get the controllability of the Purcell swimmer.

Controllability of the N-link swimmer

The third step is to generalize the previous controllability result to the N-link swimmer, whose dynamics is described by (1.19). It is easy to see that the vector fields g i generate the tangent space of the manifolds [0

, 2π] N-1 , Span(g 1 , • • • , g N-1 ) = R N-1 .
(1.30)

We can obtain the vector fields g 1 and g 2 starting from the Purcell's one defined in (1.28) as follows: we add N -2 rows of zeroes, take sticks of null length L i = 0 for 4 ≤ i ≤ N -1, while keeping the three sticks L 1 = L 3 = L and L 2 = 2L. In this case, for any

(x 1 , θ 1 ) ∈ R 2 ×[0, 2π], subsection 1.3.3 shows that g 1 (θ 1 , 0, • • • 0), g 2 (θ 1 , 0, • • • 0) and their iterated Lie brackets [g 1 , g 2 ](θ 1 , 0, • • • 0), [g 1 , [g 1 , g 2 ]](θ 1 , 0, • • • 0), and [g 2 , [g 1 , g 2 ]](θ 1 , 0, • • • 0) are linearly independent.
Therefore, the Lie algebra of the family

(g i ) i=1,•••,N-1 at the point (θ 1 , 0, • • • , 0) is equal to the tangent space T (0,•••,0,x 1 ,θ 1 ) M .
Then, by analyticity of the vector fields g i , the Orbit Theorem ensures that the Lie algebra is fully generated everywhere for a swimmer whose length of sticks verify L 1 = L 3 = L, L 2 = 2L and L i≥4 = 0. We call D (0,•••,0) , the function that maps (L 1 , • • • , L N ) to the determinant of the vectors g 1 , • • • , g N-1 and their iterated Lie brakets at the point (0, • • • , 0). Since the vector fields g i depend analytically on the sticks length L i , we get the analyticity of the function D (0,•••,0) . Thus for any L > 0, the value of D (0,•••,0) at the point (L, 2L, L, 0 • • • 0) is not null. By analyticity, it remains non null almost everywhere in R N . Therefore, we obtain that the Lie algebra has full rank for almost every swimmer. Finally, Chow Theorem gives the controllability in the Theorem 1.1.

Minimum time optimal control problem for the N-link swimmer

We present in subsection 1.4.1 the minimum time optimal control problem for the N-link swimmer, which is well defined from the controllability result proven in section 1.3. Then in subsection 1.4.2 we present the numerical method used to solve this problem.

Minimum Time Problem

For any time t > 0, we use the following notation: the state of the swimmer is

z(t) := (α 2 , • • • , α N , x 1 , θ 1 )(t), the control functions are u(t) := ( α2 , • • • , αN )(t)
and the dynamics is f(z(t), u(t)) = ∑ N-1 i=1 g i (z(t)) αi+1 (t). We now assume that the swimmer starts at the initial configuration z i , and we fix a final state z f . Our aim is to find a swimming strategy that minimizes the time to swim between the initial and the final configuration, i.e., (OCP)

       inf t f , ż(t) = f(z(t), u(t)) , ∀t ∈ [0,t f ] , u(t) ∈ U := [-1, 1] N , ∀t ∈ [0,t f ] , z(0) = z i , z(t f ) = z f .
By applying Filippov-Cesary Theorem ( [START_REF] Trelat | Contrôle optimal : théorie and applications[END_REF]) which ensures the existence of a solution of the minimum time problem for controllable systems, there exists a minimal time such that the constraints are satisfied i.e., the infimum can be written as a minimum.

Numerical Optimization

In order to solve this optimal control problem, we use a direct approach. This approach transforms the infinite dimensional optimal control problem (OCP) into a finite dimensional optimization problem (NLP). This is done with a discretization procedure on the dynamics equation summarized below:

t ∈ [0,t f ] → {t 0 = 0, . . . ,t N = t f } z(•), u(•) → X = {z 0 , . . . , z N , u 0 , . . . , u N-1 ,t f } Criterion → min t f Dynamics → (ex : Euler) z i+i = z i + h f (z i , u i ) Controls → -1 ≤ u i ≤ 1 I/F Cond. → Φ(z 0 , z N ) = 0
We therefore obtain a nonlinear programming problem on the discretized state and control variables

(NLP) min F(z) = t f LB ≤ C(z) ≤ UB
All tests were run using the software BOCOP ( [START_REF] Bonnans | Bocop -A collection of examples[END_REF]). The discretized nonlinear optimization problem is solved by the well-known solver IPOPT [START_REF] Wächter | On the implementation of a primal-dual interior point filter line search algorithm for large-scale nonlinear programming[END_REF] with MUMPS [START_REF] Amestoy | A fully asynchronous multifrontal solver using distributed dynamic scheduling[END_REF], while the derivatives are computed by sparse automatic differentiation with ADOL-C [START_REF] Walther | Getting started with Adol-c[END_REF] and COLPACK [START_REF] Gebremedhin | Exploiting sparsity in jacobian computation via coloring and automatic differentiation: a case study in a simulated moving bed process[END_REF]. In the numerical experiments, we used a Midpoint (implicit 2nd order) discretization with 1000 time steps. Execution times on a Xeon 3.2GHz CPU were a few minutes.

Numerical simulations for the Purcell's swimmer (N= 3)

We present in this section the numerical simulations regarding the Purcell swimmer (3 links). Without making any assumptions on the structure of the optimal trajectory, we obtain an optimal solution with periodic strokes. Comparing this stroke to the one of Purcell ([10], [START_REF] Becker | On self-propulsion of micro-machines at low Reynolds number: Purcell's three-link swimmer[END_REF]), we observe that it gives a better displacement speed.

In the rest of the chapter, in order to match the notations used in [START_REF] Becker | On self-propulsion of micro-machines at low Reynolds number: Purcell's three-link swimmer[END_REF], we we the following coordinates (see Fig 1 This reformulation gives the new dynamics

    β1 β3 ẋ2 θ2     = M (θ 2 , β 1 )     α2 α3 ẋ1 θ1     , M (θ 2 , β 1 ) =       -1 0 0 0 0 0 1 0 0 0 sin(θ 2 ) + cos(β 1 ) 0 1 0 -sin(θ 2 ) -cos(β 1 ) -cos(θ 2 ) 0 0 1 cos(θ 2 ) -1 0 0 0 1       .
As a result, the dynamics (1.19) reads in this case as

    β1 β3 ẋ2 θ2     = f1 (θ 2 , β 2 , β 3 ) β1 + f2 (θ 2 , β 2 , β 3 ) β3 (1.31) where for i = 1, 2 fi (θ 2 , β 1 , β 3 ) = M (θ 2 , β 1 , ) gi (θ 1 , α 2 , α 3 ) . (1.32)
Observe that since the new state variables are the image of the former ones through a one-to-one mapping, the controllability result in subsection 1.3.3 holds also for (1.31).

The classical Purcell stroke

We recall the stroke presented by Purcell in [START_REF] Purcell | Life at low Reynolds number[END_REF] in order to compare it to the optimal strategy given by our numerical results. Let us denote by ∆ θ the angular excursion, which means that β 1 and

β 3 belong to [-∆ θ 2 , ∆ θ 2 ]
. The Purcell stroke is defined by this periodic of deformation over [0, T ]:

(β 1 (t), β 3 (t)) =        ( 4∆ θ T t -∆ θ 2 , ∆ θ 2 ) 0 ≤ t ≤ T 4 ( ∆ θ 2 , -4∆ θ T t + 3∆ θ 2 ) T 4 ≤ t ≤ T 2 (-4∆ θ T t + 5∆ θ 2 , -∆ θ 2 ) T 2 ≤ t ≤ 3T 4 (-∆ θ 2 , 4∆ θ T t -7∆ θ 2 ) 3T 4 ≤ t ≤ T .
In what follows, we call the "classical" Purcell stroke the one corresponding to ∆ θ = π 3 , with T = 4∆ θ chosen in order to satisfy the constraints on the controls of (OCP), i.e., u i (t) := βi (t) ∈ [-1, 1].

Comparison of the optimal stroke and Purcell stroke

We set the initial position (x 2 , θ 2 ) = (0, 0, 0) and the final position (x 2 , θ 2 ) = (-0.25, 0, 0). Moreover we constrain the angles β 1 and β 3 in [-π 6 , π 6 ] for all time. Solving the minimum time problem numerically with the direct method gives us a solution that is actually periodic, as shown on Fig. 1.8. We notice that the xdisplacement is not monotone: during each stroke, the swimmer alternately moves forward, closer to the target, and goes partially backward. Now we isolate only one stroke from this solution, and compare it with the Purcell stroke. We show on Fig. 1.9 the angles functions β 1 and β 3 , as well as the phase portrait. Note that the time required to complete our candidate for an optimal stroke is shorter than for the Purcell one (roughly 2.5 versus 4.1). We illustrate on Fig. 1.10 the shape changes in the plane for the Purcell and optimal stroke. Finally we make a comparison between the two x-displacement, Fig. 1.11 shows the x-displacement of the swimmer with the classical Purcell stroke (dashed) and the optimal stroke (solid). Both trajectories were recomputed in Matlab using the same ODE solver, and the results for the Purcell stroke match the ones in [START_REF] Becker | On self-propulsion of micro-machines at low Reynolds number: Purcell's three-link swimmer[END_REF]. The final time t f = 15.3252 is the one given by the numeric simulation to reach x 2 = (-0.25 , 0). We see that using Purcell strokes, the swimmer only reaches (≈ -0.18, 0), which confirms that our optimal stroke allows a greate xdisplacement.

More precisely, each optimal stroke gives a x-displacement close to the Purcell stroke, however the cycle of deformation is performed in less time. Therefore, for a given time frame, more optimal strokes can be performed, leading to an overall greater displacement. In Fig. 1.11, almost 3.5 Purcell strokes are performed, while 6 optimal strokes are completed within the same time. Remark: The initial shape of the swimmer is not identical for both strategies, however the increasing gap between the two curves clearly shows that the optimal stroke is faster.

We also observe that the optimal stroke consistently gives a swimming speed better by 20% than the Purcell stroke.

Conclusions

In this chapter we study the N-link swimmer, and use the Resistive Force Theory to derive its dynamics. In this context, we prove that for N greater than 3 and for almost any N-uplet of sticks lengths, the swimmer is globally controllable in the whole plane. Then, we focus on finding a swimming strategy that leads the N-link swimmer from an fixed initial position to a given final position, in minimum time. As a consequence of the controllability result, we show that there exists a shape change function which allows to reach the final state in a minimal time. We formulate this optimal control problem and solve it with a direct approach (BOCOP) for the case N = 3 (Purcell swimmer). Without any assumption on the structure of the trajectory, we obtain a periodic solution, from which we identify an optimal stroke. Comparing this optimal stroke with the Purcell one confirms that it is better, actually giving a greater displacement speed.

Fig. 1 . 1

 11 Fig. 1.1 Coordinates of the N-link swimmer

Fig. 1 . 2

 12 Fig. 1.2 The discrete approximation by the N-link swimmer (blue curve), N = 15 of a continuous tail (red curve).

Fig. 1 . 3

 13 Fig. 1.3 The prescribed continuous wave (red curve) and its discrete approximation by the N-link swimmer (blue curve), N = 15.

Fig. 1 . 4

 14 Fig. 1.4 Flagellar beating during one period. The red curve represents the tail as described by formula 1.24 while the blue links describe the tail according to our discrete approximation.

  Moreover we use the following values for the drag coefficients • for the head, ξ head = 40.3 • 10 3 pN s m -1 , η head = 46.1 • 10 3 pN s m -1 , and ζ head = 0.84 • 10 -6 pN s m • for the links rcomposing the tail, ξ = 0.38 • 10 9 pN s m -2 , η ξ = 1.89. The graphs in Figures 1.5 and 1.

- 1 )Fig. 1 . 5

 115 Fig. 1.5 Above translational speed of the swimmer head in the tangent and perpendicular directions, and below rotational speed.

Fig. 1 . 6

 16 Fig. 1.6 Trajectory of the head of the sperm-cell during one period.

  .

Fig. 1 . 7

 17 Fig. 1.7 Purcell's 3-link swimmer.

Fig. 1 . 8

 18 Fig. 1.8 Angles and x-displacement for a whole periodic trajectory.

Fig. 1 . 9

 19 Fig. 1.9 Angles and phase portrait -Purcell stroke and optimal stroke.

Fig. 1 .

 1 Fig. 1.10 Purcell stroke (above) and optimal stroke (below).

Fig. 1 .

 1 Fig. 1.11 x displacement for one Purcell and one optimal stroke.