C Girerd
email: cgirerd@eng.ucsd.edu

A V Kudryavtsev

P Rougeot

P Renaud

K Rabenorosoa

B Tamadazte
email: brahim.tamadazte@femto-st.fr

SLAM-based Follow-the-Leader Deployment of Concentric Tube Robots

This paper deals with an original combination of navigation by simultaneous localization and mapping (SLAM) and follow-the-leader (FTL) deployment for the control of concentric tube robots (CTRs). The objective is to make possible the automatic navigation of such continuum robots within an anatomical structure, by simultaneously managing both collision avoidance and map updating of the environment, i.e., of the surrounding anatomical structures. Indeed, CTRs designed to achieve a FTL deployment to date were deployed on path identified during a planning task. In this work, this limitation is overcome by embedding, in the control scheme, an online trajectory estimation, therefore allowing an autonomous deployment of CTRs with an approximate FTL behavior. The proposed hybrid control law integrates a stability criterion to ensure at the same time CTR stability during deployment, required for a FTL deployment. The interest of the approach is demonstrated through simulation using realistic anatomical data and experiments, with a 3 degrees-of-freedom CTR.

I. INTRODUCTION

Navigating within the human body remains challenging when targeted organs are located in confined spaces or are located in a fragile environment. Several robot paradigms have been proposed in order to make such navigation possible: active catheters for vascular navigation and interventions [1], concentric tube robots (CTRs) [START_REF] Burgner-Kahrs | Continuum robots for medical applications: A survey[END_REF] including designs with embedded actuators [START_REF] Chikhaoui | Kinematics and performance analysis of a novel concentric tube robotic structure with embedded soft micro-actuation[END_REF], hyper-redundant structures such as the i 2 Snake [START_REF] Berthet-Rayne | The i2snake robotic platform for endoscopic surgery[END_REF]. CTRs interestingly offer a high compactness and the ability to deploy along predefined paths, including in a follow-the-leader (FTL) manner [START_REF] Gilbert | Concentric tube robots as steerable needles: Achieving follow-the-leader deployment[END_REF], [START_REF] Garriga-Casanovas | Complete follow-the-leader kinematics using concentric tube robots[END_REF] where the CTR body follows its tip.

Traditionally, navigation can be achieved with two approaches, either open-loop based on an offline trajectory planning [START_REF] Kuntz | Motion planning for a three-stage multilumen transoral lung access system[END_REF] using 3D models reconstructed from medical images [START_REF] Fichera | Through the eustachian tube and beyond: A new miniature robotic endoscope to see into the middle ear[END_REF], or using automatic navigation [START_REF] Dwyer | A continuum robot and control interface for surgical assist in fetoscopic interventions[END_REF]. Implementation of image-based deployment relies on the use of an external imaging system for an eye-to-hand configuration, or on an eye-in-hand one when a visual sensor is mounted at the tip of the robot. Intracorporeal navigation requires a high degree of accuracy to ensure safety. Image-guided control with eye-in-hand configuration is particularly well-suited for continuum robots [START_REF] Wu | Model-free image guidance for intelligent tubular robots with pre-clinical feasibility study: Towards minimally invasive trans-orifice surgery[END_REF]: it makes target recognition easier, allows a close visualization of the surgical site, and leads to less registration issues compared to eye-to-hand approaches.

Controlling a CTR with an eye-in-hand configuration was indeed recently demonstrated to provide a high degree of accuracy during positioning tasks with model-based [START_REF] Kudryavtsev | Eye-in-hand visual servoing of concentric tube robots[END_REF] or model-free [START_REF] Wu | Model-free image guidance for intelligent tubular robots with pre-clinical feasibility study: Towards minimally invasive trans-orifice surgery[END_REF] visual servoing schemes. However, the autonomous navigation of CTRs is still an open challenge [START_REF] Burgner-Kahrs | Continuum robots for medical applications: A survey[END_REF]. This seems however mandatory given i) the inter-patient variability which can make hazardous CTR control based on a priori knowledge which is not patient specific, ii) the possible anatomy variations between planning and inspection phases, and iii) the existence of physiological motions such as breathing-induced motions, leading to time-varying displacements and deformations of the patient anatomy. Considering the environment of CTRs in real-time during navigation is needed, and can be well illustrated by the task of the olfactory cleft inspection, where contacts with tissues are not admissible according to the clinical requirements [START_REF] Girerd | In vivo inspection of the olfactory epithelium: Feasibility of robotized optical biopsy[END_REF].

SLAM methods allow a robot to simultaneously compute the sensor pose and update a map of an unknown environment [START_REF] Durrant-Whyte | Simultaneous localization and mapping: part i[END_REF], [START_REF] Bailey | Simultaneous localization and mapping (slam): Part ii[END_REF]. SLAM has been greatly improved to obtain accurate and robust approaches offering satisfying localization and mapping. In this work, our contribution is to combine state-of-art SLAM with closed-loop control based on an eye-in-hand configuration to make safe autonomous approximate FTL deployment of CTRs in confined spaces. Safety is also considered by checking CTR stability within the controller scheme. To develop and assess this proposition, we consider a 3-tube CTR design, previously identified as relevant for navigation in the nasal cavity [START_REF] Girerd | In vivo inspection of the olfactory epithelium: Feasibility of robotized optical biopsy[END_REF]. The tube geometries of the robot were selected for FTL deployment using a dataset of nasal cavities. The control law based on visual SLAM allows to approximate the traditional FTL deployment of the CTR without any subject-specific planning task. Simultaneously, the SLAM algorithm provides a reconstruction of the anatomical environment of interest. Such algorithms have been used in combination with CTR [START_REF] Dwyer | A continuum robot and control interface for surgical assist in fetoscopic interventions[END_REF], but to the best of our knowledge, this is the first exploitation of visual SLAM for the deployment of CTRs.

The paper is organized as follows. Section II introduces required elements related to CTR kinematics as well as the stability criterion included in the developed controller. Section III focuses on the SLAM algorithm and on the formulation of the autonomous navigation controller of a CTR. Validation through simulation using realistic anatomical data thanks to a developed open-source simulator, and experiments, are described in Sections IV and V, respectively.

II. CONCENTRIC TUBE ROBOTS: BACKGROUND A. Kinematics

During a deployment in free space, as considered in this work, the kinematics of a CTR composed of n tubes is formulated as a boundary value problem governed by differential algebraic equations [START_REF] Webster | Design and kinematic modeling of constant curvature continuum robots: A review[END_REF]. Using the variables and notations reported in Table I, they are expressed as

k it ψ i = k ib k b n j=1 k jb κ i κ j sin(ψ i -ψ j ), (1)
with i ∈ {1,. . .,n} and k b = n i=0 k ib . The derivatives are taken with respect to s. The tube angles are known at their proximal ends, where they are attached to their actuators, while their torsion equals zero at their free ends (Fig. 1(a), 1(b)). This leads to the boundary conditions (2):

ψ i (0) = ψ i (β i ) -β i ψ i (0), ψ i (L i + β i ) = 0.
(2) Future 3D position of the end-effector of the CTR J inv T Jacobian matrix for the "Translating" task J inv R Jacobian matrix for the "Rotating" task J inv S Specific Jacobian matrix for the "Rotating" task λ T Gain for the "Translating" task λ R Gain for the "Rotating" task η Priority ratio between the "Translating" and "Rotating" tasks

The shape of the robot is determined by integration of Eq. (3), with p 1 the position and R 1 the orientation of the backbone Bishop frame along the robot, and ûB the Bishop frame curvature, with . denoting the conversion from R 3 to the cross product matrix.

p 1 = R 1 e 3 R 1 = R 1 ûB (3)
Eq (3) are associated to the boundary conditions of Eq (4),

p 1 (0) = 0 R 1 (0) = R z (ψ 1 (0)). (4)

B. Stability Conditions

CTRs can exhibit instabilities during deployment [START_REF] Hendrick | Designing snapfree concentric tube robots: A local bifurcation approach[END_REF], [START_REF] Peyron | A numerical framework for the stability and cardinality analysis of concentric tube robots: Introduction and application to the follow-the-leader deployment[END_REF]. Global stability criterion exist only for CTRs made of two constant-curvature tubes. Local stability criterion exist for CTRs made of any number of piecewise constant curvature tubes [START_REF] Hendrick | Designing snapfree concentric tube robots: A local bifurcation approach[END_REF]. The latter is derived from the linearization of the system of equations (1) around the equilibrium configuration to assess. The resulting system has the form

K t ψ(L 1 + β 1 ) = W 2 K t ψ(0), (5)
where

K t = diag(k 1t • • • k nt ).
The matrix W 2 depends on the tube curvatures, deployed and transmission lengths, and on their bending and torsional stiffnesses. A CTR is stable if det(W 2 ) > 0. As it is a local stability criterion, it must be verified at any time during deployment, with the corresponding deployed lengths.

C. FTL Deployment

Situations of FTL deployment were reported in [START_REF] Gilbert | Concentric tube robots as steerable needles: Achieving follow-the-leader deployment[END_REF], [START_REF] Garriga-Casanovas | Complete follow-the-leader kinematics using concentric tube robots[END_REF], [START_REF] Bergeles | Concentric tube robot design and optimization based on task and anatomical constraints[END_REF]. When the CTR is composed of tubes with constant curvatures, one requirement is to arrange them in a plane, so that they have aligned or opposite curvatures. This is the situation considered in this paper. Only the tube translations are thus used during the CTR deployment. To ensure that the CTR respects a FTL deployment, the condition det(W 2 ) > 0 must hold so that robot stays in its plane of deployment.

III. NAVIGATION USING VISUAL SLAM

Our approach is intended to use SLAM paradigm to ensure CTR deployment. The overall strategy for CTR navigation, summarized in Algorithm 1, is detailed hereafter. In the remaining of the paper, the CTR considered is made of an inner straight tube and two outer curved tubes, that conform to make a straight, curved, and straight links as illustrated in Fig. 2, compatible with a deployment in the nasal cavity which is the clinical motivation of this work [START_REF] Girerd | In vivo inspection of the olfactory epithelium: Feasibility of robotized optical biopsy[END_REF].

(a)

β 1 β 2 β 3 ψ 1 (β 1 ) ψ 2 (β 2 ) ψ 3 (β 3 ) L 3 + β 3 L 2 + β 2 L 1 + β 1 l 1 l 2 l 3 s = 0 (b) ψ 1 ψ 2 x B x 1 x 2 y B y 1 y 2 Fig. 1.
(a) 3-tube CTR straightened for ease of understanding, and (b) cross-sectional view of link 2 of the CTR with the orientation of the material base frames of tubes 1 and 2 relatively to the Bishop frame.

A. Selection of a SLAM Algorithm

Visual SLAM algorithms are designed to take advantage of the rich information available in the image data. The related literature is varied. Among the recent developments in this field, one can mention the DSO-SLAM approach [START_REF] Engel | Direct sparse odometry[END_REF]. It can be considered as an intermediate method between purely sparse methods such as ORB-SLAM [START_REF] Mur-Artal | Orb-slam: a versatile and accurate monocular slam system[END_REF] and dense approaches such as DTAM [START_REF] Lovegrove | Dtam: Dense tracking and mapping in real-time[END_REF]. In this work, we opt for the DSO-SLAM algorithm for three main reasons. First, purely sparse methods are very sensitive to the illumination changes and motion blur. Second, dense approaches generally require a GPU implementation and associated powerful hardware. Third, as reported in [START_REF] Yang | Challenges in monocular visual odometry: Photometric calibration, motion bias, and rolling shutter effect[END_REF], DSO-SLAM is more accurate than both ORB-SLAM [START_REF] Mur-Artal | Orb-slam: a versatile and accurate monocular slam system[END_REF] and LSD-SLAM [START_REF] Engel | Lsd-slam: Large-scale direct monocular slam[END_REF] in camera pose estimation, thanks to the combination of a direct probabilistic model which minimizes a photometric error, with a consistent joint optimization of all model parameters.

B. Task 1: Translating

Each tube is only actuated in translation, so the joint vector is

q = L 1 + β 1 L 2 + β 2 L 3 + β 3 .
While moving in the plane, the end-effector has 3 degrees-of-freedom (DOF) (Fig. 2): translations along x EE and z EE axes, and one rotation θ about the y EE axis. For the CTR deployment, the primary task is to move along the optical axis of the camera that coincides with z EE , i.e., the axis of the end-effector. The last link composed by the inner tube is straight, so for this task referred to with the index T, only the inner tube is actuated with velocity v, and the inverse kinematic model is:

q T = 0 1 0 0 0 0 0 0 0 0 v 0 = J inv T 0 v 0 . ( 6
)

C. Task 2: Rotating

Obviously, CTR tip translation is not sufficient. Due to the curved nature of intracorporeal passageways, rotations around its tip are needed to change the deployment direction, and avoid obstacles faced by the robot. The CTR tip corresponds to the camera center in the eye-in-hand configuration. For this task, designated by the index R, the end-effector pose is formulated as a function of link lengths. Thereby, using the geometric relations reported in Fig. 2, we obtain:

l 3 cos(θ) r sin(θ) l 1 l 3 , γ 3 l 2 , γ 2 l 1 , γ 1 θ r r cos(θ) l 3 sin(θ) z f x f z EE z c x E E y c θ = l 2 γ 2 r = 1/γ 2 p E E
x

E E z E E θ = r -r cos(θ) + l 3 sin(θ) l 1 + r sin(θ) + l 3 cos(θ) l 2 /r , (7)
where l i is the length of the i-th link of the CTR, and r, θ are the curvature radius and the arc angle of the second link. Thus, the analytical Jacobian J inv R is calculated as follows:

J inv R = ∂x E E ∂l 1 ∂x E E ∂l 2 ∂x E E ∂l 3 ∂z E E ∂l 1 ∂z E E ∂l 2 ∂z E E ∂l 3 ∂θ ∂l 1 ∂θ ∂l 2 ∂θ ∂l 3 -1 = -1 tan(γ 2 l 2 ) 1 l 3 sin(γ 2 l 2 ) 0 0 1 γ 2 1 sin(γ 2 l 2 ) 0 -1 γ 2 + l 3 tan(γ 2 l 2 )
.

(8) The Jacobian J inv R is expressed as a function of link lengths, and not in the joint variables. As introduced in [START_REF] Webster | Design and kinematic modeling of constant curvature continuum robots: A review[END_REF], those are related by a so-called specific Jacobian matrix, which takes here the expression:

q = 0 0 1 0 1 1 1 1 1 l 1 l 2 l 3 = J inv S l 1 l 2 l 3 . (9)
The "Rotating" task, q R is then expressed as:

q R = J inv S J inv R 0 0 ω , ( 10
)
where ω is the angular velocity. All the tube translations are required to perform this task.

D. Combining the Tasks

The activation of the rotating task and its direction depend on the environment. The robot will indeed start rotating when its tip gets closer to an obstacle. The switching mechanism between the two tasks is developed hereafter in five main steps. We first introduce the following notations:

• H = [h 1 h 2 • • • h N ]
is the set of 3D points given by the SLAM algorithm, that corresponds to the estimation of the surface of the environment;

• X = [x 1 x 2 • • • x M ] is a projection of points of H that are
close to the future tip position;

• p E E ∈ R 3 is the current 3D position of the end-effector; • p * E E ∈ R 3
is the future 3D position of the end-effector.

Step 1: Determination of plane π: As a first step, a plane π is built perpendicular to the optical axis of the camera and located at a distance d π from the end-effector (Fig. 3). The normalized vector z c corresponds to the camera optical axis. The point p * E E can then be expressed as:

p * E E = p E E + d π z c , (11)
and the plane π is defined through the normal z c and the point p * E E :

π : p ∈ R 3 | z c • (p -p * E E ) = 0 . ( 12
)
Step 2: Determination of X: All the 3D points h ∈ H that are located at a distance lower than d * π to π are projected on it. A set of 2D points X is obtained, as visible in Fig. 3. We define a constraint for the future position of the end-effector p * E E as follows: there should be no physical contact between the robot and the environment. Thus, the distance between X (and thus H ) and p * E E should be superior to d min X . This leads to the following expression for X:

X : {h ∈ H | abs(z c • (h -p * E E )) < d * π and (h -p * E E ) ∧ z c > d min X }. (13)
X is thus an approximation of the object contour at the future desired position p * E E of the end-effector.

Step 3: Determination of the repulsive force F: To ensure that the robot tip remains inside the object contour, a repulsive force F is generated [START_REF] Khatib | Real-time obstacle avoidance for manipulators and mobile robots[END_REF], [START_REF] Haddadin | Reactive motion generation for robots in dynamic environments[END_REF]. It is similar to the one treated by Coulomb's law: each point of X exerts a force on the point p * E E that is inversely proportional to the square of the distance between them and p * E E , so that they repel the tip of the robot (Eq. ( 14)). This is illustrated in Fig. 4. As only planar motions of the CTR are considered, the projection F of the repulsive force onto y c is used (Fig. 3).

F =          1 Card(X) Card(X) j=1 y c x j -p * E E x j -p * E E 3 if Card(X) > 0 0 otherwise (14)
Step 4: Filtering of F: The point cloud reconstructed with the SLAM algorithm can be sensitive to image noise or other external disturbances. The computation of the repulsive force is thus alike. Consequently, we filter it using a moving median filter. The force at the t-th moment of time is then given by Eq. [START_REF] Webster | Design and kinematic modeling of constant curvature continuum robots: A review[END_REF] where δ is the filtering coefficient: Step 5: Combination of tasks: Depending on the repulsive force F, the implication of q T for the "Translating" task and q R for the "Rotating" task differs during the navigation. Thus, let us introduce a parameter η = f (F), so that

F t = δF t-1 + (1 -δ)F. ( 15
)
p E E p * E E X camera CTR 3D point cloud H π x c y c z c d π d * π d X
q = λ T η q T Translating + λ R (1 -η) q R Rotating , (16)
where λ T and λ R are the task gains. The parameter η lies in the interval [0.5,1]. When it equals 0.5, the same priority is assigned to the "Translating" and the "Rotating" tasks. In the extreme case where it equals 1, only the "Translating" task is activated. The chosen mathematical function enabling such possibility is a sigmoid. In this work, we use the formulation given in [START_REF] Peyron | A numerical framework for the stability and cardinality analysis of concentric tube robots: Introduction and application to the follow-the-leader deployment[END_REF]:

η = 1.5 - 1 1 + e -a |F | , (17)
with a a scalar that controls the slope of the function. Thus, the final expression for the designed control law is:

q = λ T η J inv T 0 v 0 + λ R (1 -η) J inv S J inv R 0 0 sign(F) , (18)
where ω = sign(F) is used for the angular velocity.

E. Stability Assessment During Deployment

Stability conditions during deployment are ensured by integrating the criterion presented in section II-B. An analytical expression of det(W 2 ) is available for the CTR geometry under consideration [START_REF] Girerd | Towards optical biopsy of olfactory cells using concentric tube robots with follow-the-leader deployment[END_REF]:

det(W 2 ) = cosh(l 1 α) -α sinh(l 1 α) β 2 k 3t + β 3 k 2t k 2t + k 3t , (19)
with

α = -κ 2 κ 3 k 3b k 2b (k 2t + k 3t ) (k 1b + k 2b + k 3b )k 2t k 3t . (20)
Equation ( 19) is evaluated during deployment. The CTR navigation is stopped if det(W 2 ) < , with > 0. det(W 2 ) varies in a smooth way as linear actuation values are smoothly varied [START_REF] Hendrick | Designing snapfree concentric tube robots: A local bifurcation approach[END_REF]. The margin therefore ensures that the deployment can be stopped before the robot is unstable, i.e., det(W 2 ) < 0.

Algorithm 1 Control algorithm for the CTR navigation. compute J inv R (Eq. ( 8))

8:

compute joint velocities q (Eq. ( 18))

9:

evaluate det(W 2 ) using q

10: if (det(W 2 ) > ) then 11:
send q to the robot define plane π in front of camera (Eq. ( 12))

19:

point cloud X ← H close to π (Eq. ( 13))

20:

compute F from X and p * E E (Eq. ( 14))

21:

return filtered F (Eq. ( 15))

22: end

IV. VALIDATION BY SIMULATION A. Developed Simulator

A realistic, generic and open-source simulator is developed specifically for validations1. It implements the CTR whose parameters were described previously, as well as a face and a nasal cavity given the proposed application. The 3D models are loaded from STL files, built from CT scans of subjects to reflect the true shape of the nasal corridor. Finally, at the tip of the robot, a virtual camera is simulated with the following intrinsic parameters: focal length (α x , α y ) = (856.145, 856.145) pixels, principal point (u 0 , u 0 ) = (320, 240) pixels, which provides images in PNG format at the resolution of 640 × 480 pixels. A screenshot of the developed simulator is depicted in Fig 5.

B. Simulation Results

In order to evaluate the feasibility of the proposed approach that combines SLAM and FTL deployment, the developed simulator was used. It was configured to reproduce experimental conditions, i.e., the number, length and curvature of the tubes, the camera parameters, and the environment of navigation. This investigation by simulation allows to assess the performance of both methods and materials without, for instance, taking camera limitations into account. The parameters selected for the simulation are II. The automatic navigation starts with the tubes fully retracted in the actuation unit. They are then deployed in the nasal cavity through the nostril as illustrated in Fig. 5. The tip of the robot, with the camera placed on it, is located under the nostril of the virtual subject. The controller then deploys the robot in accordance with Algorithm 1, with an initialization step that consists in moving the tip of the robot forward to initialize the SLAM algorithm. The results presented in the following are illustrated in the provided video. Fig. 7 shows that the robot successfully goes through the nasal cavity, while simultaneously estimating the geometry of the environment. The red curve represents the path taken by the camera during navigation relatively to the point cloud. The latter is the complete set of key-points detected by the SLAM algorithm on successive images. An example of keypoints detected on a single frame is visible at the bottom-left illustration of Fig. 7, with the bottom-right image showing the same frame without the key-points. The values of η and F over time are visible in Fig. 6a. The deployed lengths of the tubes are also represented in Fig. 6b. As visible in these figures, only the "Translating" task is active during deployment step 1, at the entry of the nasal cavity. The tasks "Translating" and "Rotating" are then active during deployment step 2 to avoid the curved part of the nasal cavity, and then only the "Translating" task is active during deployment step 3 in the olfactory cleft. As shown in Fig. 6c, the robot remains stable during deployment, which allows the simulation to be fully conducted. Finally, Fig. 8 shows the superposition of all intermediate configurations adopted Superposition of all intermediate CTR shapes adopted during deployment for the simulation, with measurements of the outer envelope. An approximate FTL deployment is observed.

λ T = 1, λ R = 1, v = 0.
by the CTR, used to compute a discrete FTL deployment error [START_REF] Neumann | Considerations for follow-theleader motion of extensible tendon-driven continuum robots[END_REF]. It illustrates the desired approximate FTL behavior, with an outer envelope slightly larger than the diameter of the largest tube. The validation using simulation on realistic anatomical data is thus satisfactory, and illustrates the interest and performance of the proposed approach.

V. VALIDATION BY EXPERIMENT A. Experimental Setup

The performance of the proposed navigation approach is also assessed experimentally. The CTR actuation unit is composed of 3 motorized linear stages (stepper motors, LIMES 90-55-HSM, OWIS (Staufen, Germany)) for the tube translation. The main control software runs on a computer with GNU/Linux (Ubuntu 16.04). Another computer running Microsoft Windows 7 is used for communication with the actuation unit for compatibility reasons with the controllers. UDP socket is used for communication between the two computers. A low-cost miniature camera with a rolling shutter (Misumi MD-B802L-55) is attached at the extremity of tube 1. The impact of its weight and the radial offset induced between the tube centerline and the sensor of the camera are neglected. The camera provides 640×480 pixels images at 25 images per second. Note that the camera parameters are very similar to those used on the virtual camera for the developed simulator, except for the shutter technology.

The characteristics of the CTR tubes are visible in Table III. They allow to have a straight, curved, and straight CTR shape after assembly. The environment of navigation is an origami tunnel with non-flat walls and a checkerboard-like image on the inner surface. Fig. 12 represents the origami tunnel with the robot during the navigation task. We opted for the artificial tunnel with chessboard markers because of the performance of the miniature camera. The image quality is indeed low for standard components of one or two millimeters in size [START_REF] Yang | Challenges in monocular visual odometry: Photometric calibration, motion bias, and rolling shutter effect[END_REF]. With such experiment, we can here validate the navigation approach in presence of image noise.

B. Experimental Results

The parameters used for the experimental validation are λ T = 1, λ R = 1, v = 0.5 mm.s -1 , d π = 30 mm, d * π = 5 mm, d min X = 10 mm, δ = 0.95, a = 10 and = 0.1. The CTR tip axis is initially roughly aligned with the tunnel. Due to limitations on the tube strokes, the experiment starts with the CTR partially deployed, so that its distal end is located in the tunnel. The experimental results presented in the following are illustrated in the provided video. The evolution of η and F over time are visible in Fig. 9a. The initially deployed lengths and their variation over time are visible in Fig. 9b. Three deployment steps can be identified on these figures, with the "Translating" task active during deployment steps 1, 2 and 3, and the "Rotating" task only active during deployment step 2. During deployment, the CTR remains stable, as visible in Fig. 9c, allowing the experiment to be fully conducted. The detection of key-points by the camera is successful, as visible in the bottom-left images of Fig. 10. The color of the key-points illustrates that their relative depth is correctly estimated, with the red points being close to the camera, and the blue ones being far from it. The top image of Fig. 10 shows the 3D point cloud faced by the camera during navigation, which is a part of the tunnel wall. Compared to the navigation in the simulated environment, the navigation is here assessed with inaccuracies in the reconstructed 3D point cloud of the environment. This can be observed with the scattering of the tunnel mapping in Fig. 10. Finally, Fig. 11 represents the superposition of all intermediate robot shapes during deployment. Larger displacements of the robot body compared to the simulation are visible. This is due to Fig. 12. CTR with the RGB camera attached at its distal end successfully navigating through the origami tunnel.

the initial deployed lengths of the tubes, as our actuation system allowed a limited stroke for them. This lead to an angle between the CTR tip and the tunnel, that forced a reconfiguration of the CTR to avoid a collision at the tip. Such reconfiguration, that induces a CTR body motion, allows the tip of the robot to be tangent to the local curvature of the tunnel, which is the expected behavior during automatic navigation. Even in such sub-optimal conditions, the controller still manages to avoid contacts at the tip of the robot, with no contact observed either along the CTR body. The navigation is thus still successful with the robot crossing the tunnel, as visible in Fig. 12, so the experiment outlines the interest and robustness of the proposed approach.

VI. DISCUSSION AND CONCLUSION

In this paper, an automatic intracorporeal navigation paradigm for CTR deployment was presented. The method combines an adapted monocular-based visual SLAM combined with a CTR that can follow the leader. It was shown that the proposed method is able to manage simultaneously the automatic navigation in an unknown environment, the collision avoidance, as well as the 3D mapping of the robot's environment, i.e., the anatomical structure. The proposed methods and materials were tested both in simulation, using a specifically-developed simulator using realistic anatomical data, and in an experimental scenario using a CTR with a camera mounted in an eye-in-hand configuration.

In the presented work, the CTR of interest has a trivial inverse kinematic model. Future developments will be focused on the extension to other CTRs that present redundancies, that can or cannot follow the leader, and more generally to the navigation of any continuum robot. Their entire shape could then be optimized based on the reconstructed environment acquired during deployment, to avoid collisions.

Further work will also be undertaken to improve the time-computation of the proposed method, the image quality (i.e., using higher quality miniature cameras) as well as the performance of the intracorporeal navigation method in more realistic conditions, i.e., in human cadavers. It is also foreseen to extend this method with other types of medical imaging probes such as intravascular ultrasound (IVUS) or optical coherence topography (OCT). R

Fig. 2 .

 2 Fig.2. Geometry of the CTR under consideration. l i and γ i are the length and curvature of i-th link of the CTR.

Fig. 3 .

 3 Fig. 3. The points of H that are closer than d * π to the plane π are projected on it. The resulting set of 2D points X is used for motion planning.

FFig. 4 .

 4 Fig.[START_REF] Berthet-Rayne | The i2snake robotic platform for endoscopic surgery[END_REF]. Computation of the repulsive force used as a switch between the "Translating" and "Rotating" tasks. In (a) , the detected contour of the environment along with the circle of radius d mi n X are represented. In (b) , some points of the contour move inside the circle, creating a repulsive force visible in (c).

Fig. 5 .

 5 Fig. 5. Snapshot of the developed simulator. The camera image is displayed at the bottom-left in the snapshot. 1avkudr.github.io/visa-model-ctr

  25 mm.s -1 , d π = 3 mm, d * π = 0.4 mm, d min X = 2 mm, δ = 0.95, a = 10 and = 0.1. The characteristics of the tubes of the CTR are visible in Table

Fig. 6 .

 6 Fig. 6. (a) Normalized repulsive force F (filtered) and η, (b) inserted tube lengths, and (c) stability of the CTR over time, with the deployment steps identified and labeled, in the case of the simulation with deployment in the nasal cavity.

1 2 3 Fig. 7 .

 137 Fig. 7.Intranasal navigation in the simulated environment, with the deployment steps labeled. Top: reconstructed shape of nasal cavity, with the red lines showing the camera trajectory. Bottom: images coming from the virtual camera, with and without the detected features.

d

  t i p = 0.65 mm d ma x = 2.43 mm d b a s e = 1.60 mm

Fig. 8 .

 8 Fig. 8.Superposition of all intermediate CTR shapes adopted during deployment for the simulation, with measurements of the outer envelope. An approximate FTL deployment is observed.

Fig. 9 .

 9 Fig. 9. (a) Normalized repulsive force F (filtered) and η, (b) tube lengths inserted, and (c) stability of the CTR over time, with the deployment steps identified and labeled, in the case of the experiment with deployment in the origami tunnel.

Fig. 10 .Fig. 11 .

 1011 Fig. 10. Navigation in the origami tunnel. Top: reconstructed shape, with the red lines showing the camera trajectory. Bottom: images coming from the virtual camera, with and without the detected features.

This work has been supported by the EIPHI Graduate School (ANR-17-EURE-0002) and ANR NEMRO (ANR-14-CE17-0013).


Automatic Tip-Steering of Concentric Tube Robots in the Trachea Based on Visual SLAM Cédric Girerd 1,2 * , Andrey V. Kudryavtsev 1 * , Patrick Rougeot 1 , Pierre Renaud 3 , Kanty Rabenorosoa 1 and Brahim Tamadazte 1,4 Abstract-Concentric tube robots are continuum robots particularly well-suited for minimally-invasive surgery. Their small diameter, natural free channel and snake-like shape unlock the potential for applications that require a high degree of dexterity in confined spaces. The synthesis and planning of such robots are generally based on the patient's anatomy and the type of surgery to be performed. However, their practical deployment can be challenging. Open-loop approaches are prone to errors due to modeling accuracy or physiological motions of the patient during the procedure, and human-in-the-loop approaches require specific skills from the surgeon. In this paper, we propose an automatic tip-steering approach for concentric tube robots based on visual SLAM (simultaneous localization and mapping). This allows a closed-loop control scheme to compensate for open-loop errors, and provides to the surgeon the opportunity to focus on the medical procedure, simplifying the management of the deployment. A method is proposed in the case of a robot geometry that achieves approximate follow-the-leader deployment. It is validated by performing an automatic tip-steering task in the trachea using a realistic simulator for concentric tube robots.

Index Terms-Medical robotics, concentric tube robots, tip steering, visual SLAM.

I. INTRODUCTION

T ASK automation is an interesting area in medical robotics. In particular, the automation of intracorporeal deployment tasks in a time-varying and patient-specific environment is an exciting, yet open research topic. In this paper, we focus on patient-specific automatic deployment of concentric tube robots (CTRs). CTRs are a particular type of continuum robots [START_REF] Dupont | Design and control of concentric-tube robots[END_REF] made of precurved, nested tubes that interact in bending and torsion, defining the robot backbone shape. The latter can be changed by relative rotations and translations of the tubes. The shape variety of a CTR is limited by the initial geometry of its tubes, whose synthesis is based on surgical task requirements and patient's anatomy. While CTR tube design can take patient inter-variability into account, the shape of the robot still has to be adjusted by modification of its joint variables when used on a specific patient. This is usually performed offline during a planning task that requires 3-D preoperative images (e.g., computerized tomography scan) of the patient [START_REF] Leibrandt | Concentric tube robots: Rapid, stable path-planning and guidance for surgical use[END_REF]. A first approach based on point clouds has been proposed [START_REF] Kuntz | Planning high-quality motions for concentric tube robots in point clouds via parallel sampling and optimization[END_REF], which can adapt to changes of geometry of the anatomy during the planning, and has the potential to be used online. To lower the duration and cost of the procedure, we presented an approach that allowed online tip-steering of a CTR using monocular visual SLAM (Simultaneous Localization And Mapping) [START_REF] Girerd | SLAM-based followthe-leader deployment of concentric tube robots[END_REF]. The focus was put on the tip of the robot, for which obstacles were avoided, while the rest of the body had an approximate follow-the-leader (FTL) deployment, thanks to tube geometries that allow such behavior. In this paper, we now investigate the use of this approach in our opensource simulator (see https://avkudr.github.io/visa-model-ctr) with tip-steering of a CTR in a realistic model of trachea, for applications such as vocal cords surgery [START_REF] Renevier | Endoscopic laser surgery: Design, modeling, and control[END_REF] or laryngeal surgery [START_REF] Friedrich | Teleoperated tubular continuum robots for transoral surgery-feasibility in a porcine larynx model[END_REF], that motivate this work.

II. MATERIALS AND METHODS

A. Robot geometry

The proposed method first requires tube geometries that strictly respect FTL design rules [START_REF] Gilbert | Concentric tube robots as steerable needles: Achieving follow-the-leader deployment[END_REF], which are then deployed in an approximate FTL manner, using a tube actuation sequence computed online during the tip-steering task. Given the planar nature of the buccal cavity that leads to the trachea, the case of a planar robot made out of constant-curvature tubes arranged in a plane, with either aligned or opposite curvatures, is considered. This constrains the tube angles, that are not actuated, as opposed to their translation, that allow to change the orientation of the CTR tip in its deployment plane. One possible access path to the trachea is composed of successive straight, curved and straight sections. A CTR made out of 3 tubes is thus selected, that form a robot that has a straight, curved, and straight link. It is visible in Fig. 1, along with relevant parameters.

Let q be the set of coordinates defining the deployed lengths of the CTR tubes with respect to the base frame. In the example of Fig. 1, with tubes 1, 2 and 3 represented in blue, orange and violet, with tube 1 being the innermost tube, q = l 1 + l 2 + l 3 l 1 + l 2 l 1 , with l i the length of the CTR links, with link 1 being the most proximal link. The tube deployment sequence is obtained by a control algorithm that relies on a 3-D reconstruction of the environment, using a visual SLAM approach. For this purpose, DSO-SLAM [START_REF] Engel | Direct sparse odometry[END_REF] is selected for its ability to compute both the localization of the camera and the sparse mapping of the surrounding environment directly using photometric information instead of geometric ones (e.g., key-points, descriptors). This makes it a suitable and robust approach in the case low-textured environments as is the case for medical applications. A control law that uses the information provided by the DSO-SLAM algorithm is detailed in the following for tip-steering purposes.

B. Tip-steering task

Automatic tip steering is achieved by combining joint velocities computed from two distinct tasks: the "Translation" (T) and "Orientation" (O) tasks. The functions that represent these tasks, q T and q O , respectively, are associated to the fixed gains λ T and λ O , respectively, as well as the control gains η and (1 -η), respectively, that evolve during deployment, and depend on the presence of obstacles close to the CTR tip. The complete control law is given by Eq. (1). Each task is defined in the following.

q =   l 1 + l 2 + l 3 l 1 + l 2 l 1   = λ T η qT Translation + λ O (1 -η) qO Orientation (1)
1) Translation task: The primary task of the robot is to progress inside the anatomical area, by moving its tip forward. This corresponds to a translation motion of the inner straight tube of the CTR, using Eq. ( 2), with v the deployment speed of the inner tube.

qT =   v 0 0   (2)
2) Orientation task: When facing an obstacle, the "Orientation" task is activated, during which the CTR rotates about its tip at a velocity ω in order to change its tip insertion angle. This task is formulated in Eq. ( 3):

qO = J inv S J inv O   0 0 ω   , (3)
which involves the Jacobian matrix J inv O , given by Eq. ( 4), that represents the relationship between the CTR tip pose, i.e., translation and orientation of the tip in the plane, and the deployed link lengths. The matrix J inv O is obtained as follows:

J inv O =    ∂x EE ∂l 1 ∂x EE ∂l 2 ∂x EE ∂l 3 ∂z EE ∂l 1 ∂z EE ∂l 2 ∂z EE ∂l 3 ∂θ ∂l 1 ∂θ ∂l 2 ∂θ ∂l 3    -1 , (4)
and the inverse of the so-called Specific Jacobian matrix J inv S , given by Eq. ( 5), expresses the relationship between the link lengths and the tube lengths.

J inv S =   1 1 1 1 1 0 1 0 0   (5)
3) Combining the tasks: The balance between the "Translation" and "Orientation" tasks depends on the presence of obstacles close to the CTR tip. Only smooth deviations of the environment walls are considered, with no frontal obstacles faced at the tip. It is based on a repulsive force F , computed as follows. First, a point P * tip = P tip + d π z tip is created, as well as a plane π, normal to z tip , that contains P * tip (see Fig. 2). All the 3D reconstructed points given by the SLAM algorithm that are at a distance lower than d * π to π, and closer than a distance d min X to the (P tip , P * tip ) axis of the camera, are projected on the plane π. A force F k is computed during the control iteration k, which is inversely proportional to the square of their distances to P * tip , so that they repel the tip of the robot. As the CTR has a planar shape and obstacles are avoided on its deployment plane, only the component of F k along x tip is considered. The 3D point cloud given by the SLAM algorithm can be subject to image noise, so F k is filtered over time using the expression given in Eq. ( 6), with δ a weighting filtering coefficient.

F = δF k-1 + (1 -δ)F k (6)
The obtained filtered repulsive force F is finally used to compute η (see Eq. ( 7)), required to compute the gains η and (1 -η) of the "Translation" and "Orientation" tasks, respectively (see Eq. ( 1)). In Eq. ( 7), the scalar a controls the slope of η. As η lies in the interval [0.5, 1], the "Translation" task, weighted with the gain η, is always activated during the CTR control, while the "Orientation" task, weighted with the gain (1 -η), is only activated in the presence of obstacles.

η = 1.5 - 1 1 + e -a|F | (7)
Finally, we choose an angular velocity ω = sign(F ) for the "Orientation" task, allowing the CTR tip to be oriented in a direction that depends on the location of the obstacle faced. Furthermore, since CTRs can be subject to instabilities, a local stability criterion is implemented, with the condition det(W 2 ) > 0 that must hold during deployment, with W 2 which depends on the tube design parameters and deployed lengths [START_REF] Gilbert | Elastic stability of concentric tube robots: A stability measure and design test[END_REF]. The stability of the CTR is first assessed offline for sampled deployed tube lengths that enable the area of interest to be reached. It is then evaluated during the deployment task as a safety measure, as it is a critical condition for the CTR body to remain in its plane of deployment. The deployment is stopped if the CTR approaches instabilities, i.e. if det(W 2 ) < , with > 0 a stability margin.

III. RESULTS

The proposed control method is assessed using an inhouse simulator [START_REF] Kudryavtsev | Towards a generic simulator for continuum robot control[END_REF]. It contains an anatomical model of the trachea, the CTR, and a virtual camera attached to its tip. The characteristics of the tubes used for the simulation are reported in Table I. They allow to form a straight (0.00 mm -1 ), curved (0.02 mm -1 ) and straight (0.00 mm -1 ) links, with a curved value that is well-adapted to the considered anatomical area. As the robot tubes do not undergo any torsion during deployment, a torsionally-rigid kinematic CTR model [START_REF] Dupont | Design and control of concentric-tube robots[END_REF] Fig. 3 represents the 3-D reconstruction of the trachea by the SLAM algorithm. It illustrates the obtained reconstructed shape at the end of the deployment (see the provided video for the reconstructed shape of the trachea during the entire deployment). As visible in this figure, the robot successfully deploys inside the trachea, while reconstructing it in 3D. The bottom-left and right images depict the camera view in the represented pose, with and without the detected features. In a complementary way, Fig. 4 is a planar representation of the reconstructed trachea, that corresponds to a projection of the detected keypoints on the CTR deployment plane. It is represented at different stages of the deployment, with the red lines representing the camera path, the blue points representing the earliest keypoints detected by the camera, and the yellow ones representing the latest ones.

Fig. 5a represents the values of η and F during deployment, and Fig. 5b represents the deployed lengths of the tubes over time. Three different steps can be observed during deployment, that are labeled from 1 to 3 in these figures. During the first one, F is positive, and all the tubes are moving forward. Indeed, the robot first faces the external wall of the throat at the top in the camera frame, which forces its tip to rotate in one direction and increase the length of its curved link. During the second step, F becomes negative, and tubes 2 and 3 are moving backwards. This corresponds to the CTR tip facing the wall of the trachea at the bottom in the camera frame, which forces the tip to rotate in the opposite direction, in order to decrease the length of its curved link. This phenomenon of the CTR tip rotating in two directions in steps 1 and 2 can also be observed in Fig. 4 and in the top image of Fig. 3, with the red line representing the camera path that has an "S"-like shape. It is also visible in the provided video. Then, during the last deployment step, F equals zero, and only the "Translation" task is active, with the inner tube moving forward, while tubes 2 and 3 remain still. The deployment inside the considered anatomical area was fully conducted and was successful, with no instability encountered by the robot, as visible in Fig. 5c. Finally, the impact of robot pose errors on the deployment is evaluated. The tracking accuracy of commercial systems was assessed in [START_REF] Ameler | A comparative evaluation of steamvr tracking and the optitrack system for medical device tracking[END_REF], with reported errors in the order of magnitude of 1 mm and 1 degree in the pose estimation of the device. We thus generated Gaussian distributions of errors centered on 0, with 6σ = 1 mm and 1 degree for the translation and rotation pose errors, respectively, and truncated these distributions at ±3σ. We then performed 30 independent deployment simulations, injecting pose errors on the device frame R device , defined as the translation of the base frame R base (see Fig. 1) by -100 mm along the z-axis. Despite the pose errors added to the initial unoptimized robot pose, all 30 deployments inside the trachea were successful, which tends to confirm the overall interest of the proposed approach.

IV. CONCLUSION

In this paper, our visual SLAM-based tip-steering method was presented and implemented in a practical case, for deployment in the trachea. It allows an approximate FTL deployment of a CTR, using a control algorithm that combines a tip translation and orientation task. The tasks are weighted based on the observations of an eye-in-hand camera, located at the tip of the robot. The performance of our approach was evaluated in our in-house simulator, with successful deployments, even in the presence of pose errors. Future work will be focused on experimental evaluations in a realistic anatomical model, collision assessment along the entire robot body, extension of the proposed method to CTRs with tubes that don't strictly respect FTL design rules, and application to other types of continuum robots.

Fig. 1 .

 1 Fig.1. Geometry of the CTR under consideration. l i and γ i are the length and curvature of i-th CTR link.

Fig. 2 .

 2 Fig. 2. The points of H that are closer than d * π to the plane π are projected on it. The resulting set of 2D points X is used for motion planning.

  was implemented in the simulator. The parameters selected for the simulation are λ T = 1, λ O = 1, v = 0.25 mm.s -1 , d π = 3 mm, d * π = 0.4 mm, d min X = 2 mm, δ = 0.95, a = 10 and = 0.1. The simulation starts with the tubes partially deployed, with the CTR tip at the entry of the throat, as this region does not present any deployment challenge.

Fig. 3 .Fig. 4 .Fig. 5 .

 345 Fig. 3. Top: reconstructed 3-D shape of the environment given by the DSO-SLAM algorithm, with the red lines showing the camera trajectory. Bottom: frame obtained by the virtual camera during deployment, with and without the detected keypoints.

ACKNOWLEDGMENT

This work was supported by the Equipex ROBOTEX project (ANR-10-EQPX-44-01), the EIPHI Graduate School (ANR-17-EURE-0002) and ANR NEMRO (ANR-14-CE17-0013).


J U Alvarez-Muñoz

J Escareno

N Marchand
email: nicolas.marchand@inpg.fr

J F Guerrero-Castellanos

T Raharijaona
email: thibaut.raharijaona@univ-amu.fr

M Rakotondrabe

Quaternion Modeling and Observer-based Torque Compensation of an Aerial Manipulator

Keywords: Aerial Manipulation, Euler-Lagrange and Quaternion Modeling, Stochastic vs Deterministic State Estimation, Nonlinear Control

This paper proposes a simple solution regarding the stabilization of a quadcopter unmanned aerial vehicle endowed with a manipulator arm. The manipulator robot is attached below the rotors plane and this one induces torques producing stability issues. The present study deals with the stabilization of the full system (quadcopter and arm) by means of a set of nonlinear control techniques. First, a mathematical model is proposed for the system. Then, an attitude control, consisting on a bounded quaternion-based feedback allows the quadcopter attitude stabilization while compensating adverse torques from manipulator's motion. A simpleto-implement strategy is proposed to estimate the actual torque for compensation purposes. Then, the formulation of a nonlinear control, which drives the aerial vehicle to a desired position is presented. Both controls consist on saturation functions. Simulation results validate the proposed control strategy and compare the results with different manipulator torque estimations.

INTRODUCTION

Aerial manipulation has been an active area of research in recent times, mainly because the active tasking of Unmanned Aerial Vehicles (UAVs) increases the employability of these vehicles for various applications. For active task one would consider manipulation, grasping, transportation, etc. Unlike fixed wings UAVs, that are incapable of driving their velocity to zero, VTOL (Vertical Take-Off and Landing) vehicles with four, six, or eight rotary wings (afterwards called multirotors) are ideally suited for aerial manipulation or grasping. However, there are many challenges in aerial grasping for these vehicles. The biggest challenge arises from their limited payload. While multiple robots can carry payloads with cables or grippers [START_REF] Mellinger | Cooperative grasping and transport using multiple quadrotors[END_REF], their end effectors and grippers have to be light weight themselves and capable of grasping complex shapes. Secondly, the dynamics of the robot is significantly altered by the addition of payloads. Indeed this is also an This work was also supported by the CODE-TRACK project (ANR-17-CE05-0014-01).

attraction in assembly because aerial robots can use this to sense disturbance forces and moments, like in [START_REF] Mohammadi | Cooperative aerial tele-manipulation with haptic feedback[END_REF]. Besides, for payload transport, it is necessary that the robots are able to estimate the inertia of the payload and adapt to it to improve tracking performance. Numerous approaches have been proposed to deal with such a problem. In Orsag et al. (2013a), a Lyapunov based model Reference Adaptive Control is used to stabilize a quadrotor with a multi degree of freedom (DOF) manipulator. [START_REF] Jimenez-Cano | Control of an aerial robot with multilink arm for assembly tasks[END_REF] presents a Newton-Euler approach to model and control a quadrotor through a Variable Parameter Integral Backstepping (VPIB). [START_REF] Kim | Aerial manipulation using a quadrotor with a two dof robotic arm[END_REF] presents aerial manipulation using a quadrotor and a 2-DOF robot arm. The dynamic model of the whole system is obtained by the Euler-Lagrange formulation. Then, an adaptive sliding mode controller is designed. The effectiveness of the proposed method is experimentally showed by picking up and delivering an object. In [START_REF] Thomas | Toward image based visual servoing for aerial grasping and perching[END_REF] the problem is solved through an autonomous avian-inspired grasping method. Finally, in [START_REF] Yuksel | Differential flatness and control of protocentric aerial manipulators with any number of arms and mixed rigid-/elasticjoints[END_REF] a new class of aerial manipulator is presented, which is called protocentric. It consists on a PVTOL equipped with parallel manipulator arms attached to the Center Of Mass (CoM) of the aerial vehicle. A control law has been proposed for the case of rigid joints and validated through simulations.

The present paper proposes a bounded quaternion-based attitude controller tolerant to external disturbances via the observer-based compensation of the manipulator's torque. For this, a mathematical model considering the couplings between the aerial and manipulation subsystems is presented. Moreover, unlike the cited works, our design for attitude stabilization is based on quaternion formalism, avoiding singularities. Concerning the observer scheme, two approaches are presented to evaluate their effectiveness for the estimation of the manipulator's torque: a Luenberger observer and a linear Kalman filter.

With quaternion parametrization one proposes a constructive control law for the attitude and position stabilization. First, the design of a smooth almost globally asymptotically control law for attitude stabilization which takes into account the arm motion effects is introduced. Then, a globally asymptotically nonlinear controller for the translational dynamics is proposed. In general the control law is based in the usage of nested and sum of saturation functions. Simulation results validate the proposed strategy.

The paper is structured as follows. In section 2, the attitude model of the quadcopter with the manipulator arm is presented. Section 3 presents two strategies to estimate the angular position of each link in the robot manipulator. Then, the attitude and position control laws are formulated in section 4. Section 5 presents the simulation scenario and results. Finally, section 6 .

MATHEMATICAL MODELING

Unit quaternion and attitude kinematics

Consider two orthogonal right-handed coordinate frames: the body coordinate frame, B(x b , y b , z b ), located at the center of mass of a rigid body and the inertial coordinate frame, N (x n , y n , z n ), located at some point in the space (for instance, the earth NED frame). The rotation of the body frame B with respect to the fixed frame N is represented by the attitude matrix R ∈ SO(3) = {R ∈ R 3×3 : R T R = I, det R = 1}. The cross product between two vectors ξ, ∈ R 3 is represented by a matrix multiplication [ξ × ] = ξ× , where [ξ × ] is the well known skew-symmetric matrix.

The n-dimensional unit sphere embedded in R n+1 is denoted as S n = {x ∈ R n+1 : x T x = 1}. Members of SO(3) are often parameterized in terms of a rotation β ∈ R about a fixed axis e v ∈ S 2 by the map U :

R×S 2 → SO(3) defined as U(β, e v ) := I 3 + sin(β)[e × v ] + (1 -cos(β))[e × v ] 2
(1) Hence, a unit quaternion, q ∈ S 3 , is defined as

q := cos β 2 e v sin β 2 = q 0 q v ∈ S 3 (2)
where q v = (q 1 q 2 q 3 ) T ∈ R 3 and q 0 ∈ R are known as the vector and scalar parts of the quaternion respectively. The quaternion q represents an element of SO(3) through the map R : S 3 → SO(3) defined as

R := I 3 + 2q 0 [q × v ] + 2[q × v ] 2
(3) Remark 2.1. R = R(q) = R(-q) for each q ∈ S 3 , i.e. even quaternions q and -q represent the same physical attitude.

Denoting by ω = (ω 1 ω 2 ω 3 ) T the angular velocity vector of the body coordinate frame, B relative to the inertial coordinate frame N expressed in B, the kinematics equation is given by q0 qv = 1 2

-q T v I 3 q 0 + [q × v ] ω = 1 2 Ξ(q)ω (4)
The attitude error is used to quantify mismatch between two attitudes. If q defines the current attitude quaternion and q d the desired quaternion, i.e. the desired orientation, then the error quaternion that represents the attitude error between the current orientation and the desired one is given by q e := q -1 d ⊗ q = (q e0 q T ev ) T

(5) where q -1 is the complementary rotation of the quaternion q which is given by q -1 := (q 0 -q T v ) T and ⊗ denotes the quaternion multiplication, Shuster (1993).

Model of a quadcopter carrying a manipulator arm

The attitude dynamics and kinematics for the quadcopter have been reported in many works e.g. [START_REF] Guerrero-Castellanos | Bounded attitude stabilization: Real-time application on four-rotor mini-helicopter[END_REF]; [START_REF] Alaimo | Mathematical modeling and control of a hexacopter[END_REF]. In these works it is considered that the quadcopter mass distribution is symmetric. However, the mass distribution of a quadcopter with a manipulator is no longer symmetrical and varies with the movement of the arm. Consider a quadcopter with a manipulator arm with n links attached to its lower part. If the dynamics of the arm is neglected, the attitude kinematics and dynamics is given by q0 qv = 1 2 Ξ(q)ω ( 6)

J ω = -ω × Jω + Γ T (7) where J ∈ R 3×3 is the inertia matrix of the rigid body expressed in the body frame B and Γ T ∈ R 3 is the vector of applied torques. Γ T depends on the (control) couples generated by the actuators, the aerodynamic couples such as gyroscopic couples, the gravity gradient or, as in the case of the present work, the couple generated by the movement of a robot manipulator placed under the body. Here, only the control couples, gyroscopic couples and the couple generated by the manipulator is considered in the control design. Consequently, Γ T = Γ + Γ arm + Γ G (8) where Γ and Γ G are the control torques and the vector of gyroscopic couples, [START_REF] Castillo | Real-time stabilization and tracking of a four rotor mini rotorcraft[END_REF], respectively. On the other hand, the vector Γ arm = (τ 1 , τ 2 , τ 3 ) T is the torque exerted by the total propulsive force being applied at the quadcopter geometric center which is displaced from the center of mass.

The total torque exerted by the manipulator arm can be found from its dynamic model by:

M (θ ai ) θai + C(θ ai , θai ) θai + N (θ ai , θai ) = τ i (9)
where θ ai , θai and θai are the angular position, velocity and acceleration vectors, respectively on each link, M ∈ R 3×3 is the manipulator inertia matrix, C is a matrix containing the Coriolis and centrifugal force terms, N includes gravity terms acting at the joints and τ i is the vector of the actuator torques.

To derive the dynamics of the system, the Lagrange equation is used

d dt ∂L ∂ θai - ∂L ∂θ ai = τ i (10) L = K -U (11
) where K and U are the total kinetic and potential energy of the system and τ i represents the generalized force.

The kinetic energy can be computed from:

K = 1 2 m i v T i v i + 1 2 I i θ 2 ai ( 12
)
where m i , v i and I i are the mass, the velocity vector and the inertia moments for every link in the manipulator, respectively.

Finally, the gravitational forces on the manipulator are computed from:

N (θ ai , θai ) = ∂U ∂θ ai ( 13
)
where U : R n → R is the potential energy of the manipulator. For the two link-manipulator under consideration here, the potential energy is given by U (θ ai ) = U (θ a1 ) + U (θ a2 ) (14) Since each link has its own frame, the computed torque has a rotation matrix according to the axis in which the links are moving. Besides, the quadcopter rotation matrix is equally considered. Then, the computation of the manipulator's torque considering the dynamics is given by: Γ arm = 1 + R 1 τ 2 ) (15) where R was given by (3) and R 1 is a rotation matrix on the y axis.

In this case, let consider the scheme in Fig. 1, which shows a 2-dof arm manipulator, where l c1 and l c2 are the distances from the respective joint axes to the center of mass of each link, l 1 and l 2 are the length of the links, and θ ai measures the angular displacement in the x axis. Then, since servomotors are used as actuators for the manipulator arm, these can be considered as first order systems. For this, a parameter identification is performed. In general, the found system has the form:

θ aim = Ku(t) -a θaim a ( 16
)
where θ aim and θaim are the angular position and velocity of the servo, a is a time constant linked to the time response of the servo, K is the gain of the system and u(t) is the velocity input generated by the torque.

ESTIMATION STRATEGY AND TORQUE RECONSTRUCTION

This section presents two observers schemes to estimate joints angular displacements: a Luenberger state observer (LSO) and an augmented-state Kalman filter (ASKF).

The estimation of joint angles in the manipulator arm is computed combining the data coming from the first order model of the manipulator actuators and the end effector position data tracked by a positioning system, since the first order model does not fully describe the behavior of the arm manipulator (non-modeled dynamics, actuators malfunction, etc.). Then, in order to know the manipulator links angular positions with respect to the base body (quadcopter), the inverse kinematics of the manipulator arm is used.

Manipulator links Luenberger observer design

The objective of the observer is to estimate the angular positions and velocities on the arm manipulator, using the expression in (16).

A positioning system is used to obtain the end-effector position information and then, the inverse kinematics is applied to know the angle on each link. With this in consideration, the expression that describes a link angle is given by: θ aV = θ a + µ V (17) where θ aV is the estimated angle computed with the positioning system and the inverse kinematics, θ a is the real angle and µ V is a noise of minimal value.

On the other hand, the observer allows the computation of the angular velocity. The expressions that represent the observer are given by: θδ = a θδ + Ku δ + L(θ aV -θ) (18) θ = θV (19) where θδ is the estimated angle on a link in the manipulator, a and K are parameters of the first order system previously presented and L is a positive tuning parameter.

Augmented-state Kalman Filter

An augmented-state discrete linear Kalman filter (AS-LKF) is designed regarding the estimation of the disturbances arising during the arm's planar joint displacement. The disturbed nominal model ( 16) is extended to include the velocity couplings between the aerial robot and the robotic arm θaim = -aθ aim + b i u i (t) + ρ(t) i (20) where ρ(t) corresponds to external disturbances (e.g. wind gusts), and whose dynamics is considered as ρ ≈ 0 (21) Thus, it is possible to rewrite (20) and ( 21) into state-space representation θaim ρi = -a i ρ i 0 0

θ aim ρ i + k i 0 u(t) (22)
The Linear Kalman Filter (LKF) is derived from a continuous system ẋ(t) = Ax(t) + Bu(t) E α(t)β(t) T = 0 (26) Including the disturbance as state the transition matrix is augmented and the state-space equations are written as

+ P ρ + Mα(t) y(t) = Cx(t) + β(t) (23)
Ẋe (t) = AX e (t) + B(u i ) + Mα i (27)
Y e (t) = CX e (t) + β i (28) In the previous regard, let us modify the extended model of the servomotor ( 22) θaim ρi = -a i ρ i 0 0

θ aim ρ i + k i 0 u(t) + 1 0 0 1 α θ α ρ (
29) Finally, the classical linear Kalman filter is applied to (29).

Torque Reconstruction

Based on the information provided by the LSO and ASKF presented before, it is possible to compute the manipulator torque from the equation ( 15) and use this term as Γ arm into the attitude control law (32).

ATTITUDE AND POSITION CONTROL

Attitude Control: Problem statement

The objective is to design a control law which drives the quadcopter carrying the manipulator arm to attitude stabilization. Let q d denotes the constant quadcopter stabilization orientation, then the control objective is described by the following asymptotic conditions q → [±1 0 0 0] T , ω → 0 as t → ∞ (30) Furthermore, it is known that actuator saturation reduces the benefits of the feedback. Then, besides the asymptotic stability, the control law also takes into account the physical constraints of the system in order to apply only feasible control signals to the actuators.

Attitude control with manipulator arm

In this subsection, a control law that stabilizes the system described by ( 6) and ( 7) is proposed. The goal is to design a control torque that is bounded. Theorem 4.2. Consider a rigid body rotational dynamics described by ( 6) and ( 7) with the following bounded 32) with i ∈ {1, 2, 3} and where σ Mi1 and σ Mi2 are saturation functions. Assuming δ i < M i2 -M i1 and M i1 ≥ 3λ i ρ i . λ i and ρ i are positive parameters. Then the inputs (32) asymptotically stabilize the rigid body to the origin (1 0 T 0 T ) T (i.e. q 0 = 1, q v = 0 and ω = 0) with a domain of attraction equal to S 3 × R 3 \ (-1 0 T 0 T ) T .

control inputs Γ = (Γ 1 Γ 2 Γ 3 ) T such that Γ i = -σ Mi2 (Γ armi + σ Mi1 (λ i [ω i + ρ i q i ])) (

Position control: Problem statement

The objective is to design a control law which stabilizes the quadcopter to a desired position. Once the control law has stabilized the attitude of the system, the position control law should stabilize the quadcopter to a desired position, lim t→∞ (p, v) = (p d , 0).

Position stabilization strategy

The schematic representation of a quadcopter carrying a manipulator arm can be seen in Fig. 3, where the inertial reference frame N (x n , y n , z n ), the body reference frame B(x b , y b , z b ), the force u (thrust) and the weight vector mg are depicted. The dynamics of the whole system is obtained with the Newton-Euler formalism and the kinematics is represented using the quaternions formalism, given by Σ T :

     ṗ = v m T v = -m T g + R 0 0 u (33) Σ O : q = 1 2 Ξ(q)ω J ω = -ω × Jω + Γ T (34)
where p and v are linear position and velocity vectors, m T is the total mass of the system, g is the gravity, R is the rotation matrix, given in (3). Note that the rotation matrix R can also be given in function of the Euler angles, given in [START_REF] Castillo | Real-time stabilization and tracking of a four rotor mini rotorcraft[END_REF]. Now, assume that using the control law (32) one can stabilize the yaw dynamics, that is ψ = 0, then after a sufficiently long time, system (33) becomes:

ṗx ṗy ṗz = v x v y v z , ( 35
) vx vy vz =   -u m T sin θ u m T sin φ cos θ u m T cos φ cos θ -g   (36)
With an appropriate choice of the target configuration, it will be possible to transform ( 35)-( 36) into three independent linear triple integrators. For this, take

φ d := arctan r 2 r 3 + g , θ d := arcsin -r 1 r 2 1 + r 2 2 + (r 3 + g) 2 (37)
where r 1 , r 2 and r 3 will be defined after. Then, choose as positive thrust the input control

u = m T r 2 1 + r 2 2 + (r 3 + g) 2 (38)
By taking ( 35) and ( 36), it follows:

Σ x : ṗ1 = p 2 ṗ2 = p 3 ṗ3 = r 1 Σ y : ṗ4 = p 5 ṗ5 = p 6 ṗ6 = r 2 Σ z : ṗ7 = p 8 ṗ8 = p 9 ṗ9 = r 3 (39)
Note that u will be always positive, and u ≥ mg, in order to compensate the system's weight.

Since the chains of integrators given in (39) have the same form, a control law can be proposed as in [START_REF] Cruz-José | Estabilización global de mini naves aéreas tipo vtol[END_REF], and can be established by the next theorem: Theorem 4.3. Consider the quadcopter translational dynamics expressed in (35)(36). Then, the thrust input u, with control laws r 1 , r 2 , r 3 as in (40), where σ M1 (•) is defined in (31) with M 1 = 1 and ς i are given by (41), a (1,2,3) , b (1,2,3) , c (1,2,3) > 0 tuning parameters.

r 1 := -ς 1 {a 3 σ M 1 [ 1 ς 1 (a 2 p 1 + p 2 + p 3 )] + a 2 σ M 1 [ 1 ς 1 (a 1 p 2 + p 3 )] + a 1 σ M 1 [ 1 ς 1 (p 3 )]}, r 2 := -ς 2 {b 3 σ M 1 [ 1 ς 1 (b 2 p 4 + p 5 + p 6 )] + b 2 σ M 1 [ 1 ς 2 (b 1 p 5 + p 6 )] + b 1 σ M 1 [ 1 ς 2 (p 6 )]}, r 3 := -ς 3 {c 3 σ M 1 [ 1 ς 1 (c 2 p 7 + p 8 + p 9 )] + c 2 σ M 1 [ 1 ς 3 (c 1 p 8 + p 9 )] + c 1 σ M 1 [ 1 ς 3 (p 9 )]} (40) ς 1 = r1 /(a 1 + a 2 + a 3 ), ς 2 = r2 /(b 1 + b 2 + b 3 ), ς 3 = r3 /(c 1 + c 2 + c 3 ) (41)
The proof of this Theorem is not presented here, but it can be derived from the seminal work of [START_REF] Marchand | Global stabilization of multiple integrators with bounded controls[END_REF], [START_REF] Teel | Global stabilization and restricted tracking for multiple integrators with bounded controls[END_REF] and [START_REF] Johnson | Nested saturation with guaranteed real poles[END_REF]. Remark 4.4. In the above Theorem, the stabilization goal is the origin. In the case where the asymptotic condition is different from the origin, the variables p 2 , p 5 , p 8 should be replaced in the control law ( 40) by e 1 = p 2 -p d x , e 2 = p 5 -p d y , e 3 = p 8 -p d z , respectively. In this case p d

x , p d y , p d z represent the desired position in the space.

SIMULATION RESULTS

This section presents numerical results to show the effectiveness of the proposed control-estimation strategy of a quadrotor having a 2-DOF robotic manipulator. The simulation model features the parameters depicted by the Table 1.

Simulation scenario

A set of simulations were carried out in order to compare the performance of the proposed control-estimation strategies, scenario 1 and scenario 2, respectively. For the attitude control law (32) and the position control law given in Fig. 4 shows the simulation results using the two estimation strategies: Luenberger observer and Kalman filter.

From top to bottom it shows the real and the estimated manipulator torques, the torque errors, and the angular and linear positions of the aerial robot during the simulation.

A general comparison can be made from the obtained results. First, the Luenberger observer or the Kalman filter, estimate the torque generated by the manipulator arm, then this estimation is used by the attitude control law as a feed-forward term in order to compensate the perturbation. From the results, the Kalman filter allows a better torque estimation. Even more, the Kalman filter is able to estimate the external constant disturbance. Consequently,

CONCLUSIONS

In this paper, a novel model for a quadcopter carrying a manipulator arm was proposed. Besides, a control law was designed to asymptotically stabilize the attitude and position of the system. Moreover, this work has presented two strategies for aiding the solution through the design of a Luenberger observer or a Kalman filter to develop a feed-forward term which allows the estimation of the moments and torques exerted by the manipulator. Simulation results show the effectiveness of the proposed estimationcontrol strategy face to the disturbances coming from the manipulator. With the use of two strategies, a comparison was performed, showing that the usage of the Kalman filter allows a better stabilization of the system. As a future work, experimental torque and mass estimations for picking up and delivery an object will be pursued.

Figure 1 .

 1 Figure 1. Manipulator arm with two degrees of freedom.

Figure 2 .

 2 Figure 2. Observer-based manipulator's torque reconstruction that considers that the pair AC verifies the controllability property. The signals α (unmodeled dynamics and parametric uncertainties) and β (sensors noise) stand for a white Gaussian random process with zero-mean (E [α(t)] = 0) and E [β(t)] = 0)) with constant power spectral density (PSD) W (t) and V (t) defining respectively, the process covariance matrix, Q = E α(t)α(t + τ ) T = W ∆(τ ) (24) , and the sensor covariance matrix,R = E β(t)β(t + τ ) T = V ∆(τ )(25) It is also assumed that both stochastic processes are not correlated, i.e.E α(t)β(t) T = 0 (26) Including the disturbance as state the transition matrix is augmented and the state-space equations are written as

Figure 3 .

 3 Figure 3. Schematic configuration of a quadrotor carrying a manipulator arm.

  Definition 4.1. Given a positive constant M , a continuous, nondecreasing function σ M : R → R is defined by (1)σ M = s if |s| < M ; (2)σ M = sign(s)M elsewhere. (31) Note that the components of Γ armi are always bounded, i.e. | Γ armi |< δ i . Then, one has the following result.

  40), where max|Γ armi | = 0.65N m and δ i = 0.1, the next parameters values are proposed: M 11,21,31 = 1.5, M 12,22 = 1, M 32 = 0.8, λ 1,2 = 0.035, λ 3 = 0.033, ρ 1,2 = 10.5 and ρ 3 = 11. For the control (40),a 1 = b 1 = 5.2, c 1 = 0.45, a 2 = b 2 = 2.3, c 2 = 2, a 3 = b 3 = 0.1, c 3 = 0.05 and r1,2,3 = 6. The simulations consists on two parts. First, the links of the manipulator arm are initialized to θ ai = (0 • 0 • ) T and the quadrotor is driven to the position p d = (0 0 1) T , then, at time 7s the links change to θ ai = (45 • 45 • ) T , at 15s the manipulator links perform a signal tracking given by θ ai = ((30 sin(2t) + 40) • (4 sin(1.5t) + 5) • ) T , at 22s the angles change to θ ai = (0 • 90 • ) T . Finally, at time 26s a constant disturbance is exerted to the manipulator and the simulation finishes.

Figure 4 .

 4 Figure 4. General behavior of the system during the simulation. the computed control signals allow a better orientation stabilization and then, the non-desired displacement along the 3d space is reduced.

Table 1 .

 1

	System	Description	Value	Units
		Mass (m)	590	g
		Distance (d)	25	cm
	Quadcopter	Inertial moment x (J φ )	0.0039 Kg • m 2
		Inertial moment y (J θ )	0.0039 Kg • m 2
		Inertial moment z (J ψ ) 0.0073 Kg • m 2
		Mass ma	160	g
	Manipulator	Length links l 1 and l 2	12	cm


Pauline Butaud

Morvan Ouisse
email: morvan.ouisse@femto-st.fr

Gaël Chevallier

Emmanuel Foltête

Design of thermally adaptive composite structures for damping and stiffness control

Keywords: damping, viscoelastic, composite, thermally adaptive

Viscoelastic materials are widely used to control vibrations. However, their mechanical properties are known to be frequency and temperature-dependent. Thus, in a narrow frequency bandwidth, there is an optimal temperature that corresponds to a maximum loss factor and it is tricky to get a high damping level over a wide frequency range. Furthermore, an optimal temperature for a maximum structural damping leads to a poor static stiffness because the peak of the loss factor is obtained during the glass transition when the storage modulus is decreasing. Additionally, in industrial applications, the requirements might change according to the system life-cycle. For instance, the stabilization functions that are used for optronics applications require high stiffness for positioning steps, and high damping for filtering functions. To achieve this goal, engineers usually use several viscoelastic materials with functionally graded damping properties. This allows obtaining a high loss factor over a wide frequency range. This solution is however not adaptive. In order to be able to adjust the properties in real time, we suggest in this paper to use a single material which properties are functionally graded thanks to a non-homogeneous temperature field over the structure. A composite structure has been numerically designed integrating a viscoelastic core and a heat control device. The optimal temperature field has been obtained based on the static and dynamic elastic strain energy densities that reflects the compromise between structural damping over a wide frequency band and high static rigidity.

INTRODUCTION

Vibration and noise reduction in dynamic systems can be achieved by a number of means including damping, which can be split in two major categories, namely passive damping and active damping. Their efficiency is qualitatively represented in Fig. 1.(a-b). Passive damping does not require external energy. The dissipation energy of the structure of interest is increased by adding friction devices, or electromechanics damping, or with viscoelastic materials (polymers) which have the natural capacity to dissipate the vibration energy under deformation. This kind of treatments lead to reliable, low cost and robust vibration control [START_REF] Rao | Recent applications of viscoelastic damping for noise control in automobiles and commercial airplanes[END_REF][START_REF] Trindade | Experimental analysis of active-passive vibration control using viscoelastic materials and extension and shear piezoelectric actuators[END_REF] . However, the efficiency of the passive damping devices is quite low (Fig. 1.(a)) and can be dependent of the temperature and the frequency. On the other hand, active control, using for example piezoelectric materials [START_REF] Sunar | Recent advances in sensing and control of flexible structures via piezoelectric materials technology[END_REF] or adaptive dry friction, [START_REF] Bouaziz | Vibration reduction of an assembly by control of the tightening load[END_REF] has many advantages such as high precision, performance controllable systems and a quite significant efficiency (Fig. 1.(b)). However active controllers are very sensitive to variations and uncertainties of system parameters [START_REF] Trindade | Piezoelectric active vibration control of damped sandwich beams[END_REF][START_REF] Gaul | Friction control for vibration suppression[END_REF] even if distributed strategies can enhance performances while increasing robustness [START_REF] Collet | Adaptive metacomposites for vibroacoustic control applications[END_REF][START_REF] Tateo | Design variables for optimizing adaptive metacomposite made of shunted piezoelectric patches distribution[END_REF] . In this work, we want to demonstrate that the tuning of the damping properties of a sandwich composite structure can be achieved through the control of the temperature field of the viscoelastic core and can be adaptated to its frequency environement (Fig. 1.(c)). A Shape Memory Polymer (SMP), which belongs to the class of smart materials, is chosen as viscoelastic core for its high damping power. These materials have found growing interest [START_REF] Hager | Shape memory polymers: past, present and future developments[END_REF] throughout the last decades. The SMPs have the ability of changing their shape in response to an external stimulus, most typically thermal activation. It turns out that efficient Shape Memory effects are associated to fast transitions between stable glassy and rubbery states with large elasticity gap, inducing high loss factor values at glass transition [START_REF] Tsai | Shape memory effects of poly (ethylene terephthalate-coethylene succinate) random copolymers[END_REF][START_REF] Barwood | The effect of organoclay addition on the properties of an acrylate based, thermally activated shape memory polymer[END_REF] . In SMP-related open literature, various experimental results provide loss factors values at glass transition. The typical values which are reported vary between 0.5 12 and more than 2.5, [START_REF] Ortega | Structureproperty relationships in photopolymerizable polymer networks: Effect of composition on the crosslinked structure and resulting thermomechanical properties of a (meth) acrylate-based system[END_REF] most of them being between 1 and 2 [START_REF] Tandon | Durability assessment of styrene-and epoxy-based shapememory polymer resins[END_REF][START_REF] Xie | Facile tailoring of thermal transition temperatures of epoxy shape memory polymers[END_REF][START_REF] Ellson | Tunable thiol-epoxy shape memory polymer foams[END_REF][START_REF] Butaud | Investigations on the frequency and temperature effects on mechanical properties of a shape memory polymer (veriflex)[END_REF] . This intrinsic property has been investigated for the design of composite structures for damping applications in a previous work of Butaud et al. [START_REF] Butaud | Sandwich structures with tunable damping properties: On the use of shape memory polymer as viscoelastic core[END_REF] : a sandwich structure composed of aluminum skins and SMP core was designed and the ability of the SMP to reach very high loss factor values at glass transition was used to control the wave propagation. Results and perspectives are recalled in section 2 of this paper, whose objective consists in pushing the boundaries of this previous study where an awesome damping was observed at the cost of a decrease of the structure stiffness. Instead of applied an homogeneous temperature control of the sandwich, we propose here to control the SMP core with a non-homogeneous temperature field in order to find a compromise between stiffness and damping. The organization of this paper is as follows. The damping properties of the Shape Memory Polymer core sandwich are described in Section 2. Section 3 shows the tuning strategy of the SMP damping properties by the non-homogeneous temperature field in order to obtain a stiffness/damping compromise. The results of this approach are presented in Section 4. Section 5 gives some final conclusions and remarks.

DAMPING PROPERTIES OF SMP CORE SANDWICH

Elaboration procedures, mechanical characterization, and modelling of the viscoelastic behavior of the SMP of interest, can be found in Butaud et al. [START_REF] Butaud | Experimental investigations on viscoelastic properties of a shape memory polymer[END_REF] The tBA/PEGDMA, studied by Srivastava et al. [START_REF] Srivastava | Thermally actuated shape-memory polymers: Experiments, theory, and numerical simulations[END_REF] is chosen for this study. It is elaborated at the FEMTO-ST Department of Applied Mechanics following the procedure described in Yakacki et al. [START_REF] Yakacki | Unconstrained recovery characterization of shape-memory polymer networks for cardiovascular applications[END_REF] The SMP has been synthesized by manually mixing 95 wt% of tBA with 5 wt% of the crosslinking agent PEGDMA, and 0.5 wt% of the total weight of DMPA photoinitiator. The dynamic mechanical properties of the SMP have been experimentally measured, a viscoelastic model allowing the description of the rheological properties of the tBA/PEGDMA has been proposed and implemented in a FE software for computing the vibration response of a sandwich structure embedding a SMP layer 22 (see Fig. 2). This composite sandwich has been experimentally and numerically studied. The sandwich was in free-free conditions, a broadband random excitation was applied between 100 and 10000 Hz. The details of the experimental tests and the damping control of the structure can be found in Butaud et al. [START_REF] Butaud | Sandwich structures with tunable damping properties: On the use of shape memory polymer as viscoelastic core[END_REF] and the main informations are recalled here.

The Fig. 3 illustrates the damping capacities of the SMP, through the frequency response functions of the sandwich composite at four representative temperatures. At 20 • C, the SMP is in the glassy state, it has very large storage modulus (E = 2100 M P a) but very low damping (tan(δ) = 0.01). The sandwich structure is then rigid but not damped and the FRF exhibits resonance peaks. At 50 • C, the SMP is close to the glass transition. The storage modulus is around 1000 MPa, and the loss factor has an average value of 0.2 over the frequency range. Thus the structure heated at 50 • C is comparable to a classical sandwich composite exhibiting damping effect through viscoelastic behavior of the core. At 80 • C, the tBA/PEGDMA is in the transition state where the material changes from a glassy state to a rubbery state, the material modulus decreases while the loss factor increases up to 2.5 which leads to the cancelling of the sandwich structure resonances. Almost no dynamic effects are visible in the whole frequency range: the curve is completely smoothed, the SMP damping is really efficient. In the rubbery state, at 130 • C, the storage modulus has a low value. The loss factor varies with frequency from 0.2 (at 100 Hz) to 2 (at 10 kHz), hence the sandwich FRF is not as damped as at 80 • C, and some resonances are visible at low frequencies. Fig. 3 shows the impressive damping properties of the tBA/PEGDMA which are furthermore controllable by tuning the temperature of the sandwich core using the rainbow calibration curve. [START_REF] Butaud | Sandwich structures with tunable damping properties: On the use of shape memory polymer as viscoelastic core[END_REF] A good correlation is observed between the experimental results and the simulation results [START_REF] Butaud | Sandwich structures with tunable damping properties: On the use of shape memory polymer as viscoelastic core[END_REF] hence validating of the viscoelastic model chosen for the FE simulation. The drawback of this damping temperature control is the impact on the static stiffness of the structure. Indeed at 20 • C the SMP is rigid with a storage modulus around 2000 MPa but at 80 • C the storage modulus is really low around 1 MPa. The next section presents, through a numerical study, the proposed strategy to preserve this impressive damping while preserving the static stiffness.

TEMPERATURE FIELD OPTIMIZATION

Simulating the temperature control of the SMP core in a sandwich structure for vibration damping remains costly since the non-homogeneous temperature field requires the use of a fine 3D FE model. In order to demonstrate the efficiency of the proposed strategy, we focus on a "slice" of the sandwich, corresponding to a 2D composite beam with plane strain assumption. The computational cost is hence drastically reduced, and the concept will be easily adapted to 3D models. The beam is clamped at x = 0 and a distributed harmonic force f = 1 N/m 2 is applied at x = 200 mm (Fig. 4). The Fig. 5 shows the spatial mean displacement of the beam, at 20 • C and at 77 • C (the temperature for an optimal damping in the frequency range of interested). At 20 • C, the sandwich beam is rigid but not damped and the dynamic response exhibits four resonance peaks at 72, 445, 1215 and 2298 Hz. At 77 • C, no dynamic effect is visible, the curve is completely smoothed. However, at low frequencies, the displacement amplitude is 10 times larger at 77 • C than at 20 • C. The previous works of Butaud et al. [START_REF] Butaud | A smp core sandwich composite with tunable damping properties[END_REF] have highlighted the interest to apply a non-homogeneous temperature field on the SMP core, indeed thanks to an arbitrary non-homogeneous temperature field, a first compromise was found between static stiffness and damping. The work presented here tries to find an optimal non-homogeneous temperature field.

Numerical model definition

The structure is meshed with 2D quadratic elements with one element in the thickness of the aluminum skins and three elements in the SMP core thickness. The Finite Element model has 1,000 elements and 8,822 degrees of freedom. The materials are assumed to be isotropic. The aluminum material properties for the skins and the SMP core material properties are taken from Butaud et al. [START_REF] Butaud | Sandwich structures with tunable damping properties: On the use of shape memory polymer as viscoelastic core[END_REF] The structure is studied in the [1-3,500] Hz frequency range of interest.

From elastic strain energy density to temperature field

The idea proposed in this paper is to define the temperature field on the basis of the elastic strain energy density of the static and dynamic structure states. The elastic strain energy density is defined by

W = 1 2 (σ : ε), (1)
with σ the stress tensor and ε the strain tensor.

In the following, W is estimated only on the SMP core, a mean is done on the nodes of each x coordinate. W is evaluated when the SMP is in the glassy state at 20 • C. The quantities of interest are normalized between 0 and 1 for an easier comparison. W will be noted W S for the static state and W D for the dynamic state.

Static stiffness

With the objective to ensure a minimum static stiffness to the structure, the static response of the structure has to be considered when searching for a compromise. The static elastic strain energy density of the structure is computed with the FE model. On the basis of this static behavior, the objective consists in defining a temperature field such that the stiffness remains unchanged. This can be achieved only by keeping the zones in which W S is large at low temperature. On zones where the static strain energy density is low, a higher temperature may be used to increase the damping. It has been chosen to define only two temperatures (ie "cold" at 20 • C and "hot" at 77 • C) with an arbitrary threshold τ = 0.3 corresponding to W S = 1/τ . Hence all areas where W S > 1/τ will be kept at 20 • C while others will be heated at 77 • C. This strategy is illustrated in Fig. 6, using the inverse of W S as quantity of interest (hence all areas with 1/W S > 1/τ will be heated).

The choice has been done, up to now, to establish a discrete temperature field in order to avoid intermediate temperatures as 50 • C since, as seen in section 2, these intermediate temperatures are neither interesting in terms of damping nor in terms of stiffness. The non-homogeneous temperature field found here is so 20 • C for 0 x 0.13 m and 77 • C for 0.13 x 0.2 m. Obviously, in the practical implementation of the temperature field, no discontinuity will be observed: depending on the thermal properties of the core and skins of the composite, a smoothing of the thermal field will be observed. This point is currently under investigation.

Dynamic damping

The temperature defined in the previous section is based on the static behavior of the structure. In this section, we propose a similar strategy to define some temperature fields adapted to the dynamic behavior of the structure. The responses of the first four modes are assumed to be representative of the dynamic state. The elastic strain energy densities are computed for the four first modes of the structure. On the basis of these dynamic behaviors, the objective consists in defining a temperature field such that the damping is increased. This can be achieved only by applying high temperature on the zones in which W D is large. On zones where the dynamic strain energy density is low, a smaller temperature may be used to increase the stiffness. As previous, it has been chosen to define only two temperatures (ie "cold" at 20 • C and "hot" at 77 • C) with an arbitrary threshold τ = 0.3 corresponding to W D = τ . Hence all areas where W D < τ will be kept at 20 • C while others will be heated at 77 • C. This strategy is illustrated in Fig. 7, using W D as quantity of interest (hence all areas with W D > τ will be heated).

Temperature fields for a compromise between static stiffness and dynamic damping

We propose now a strategy for combining the temperature fields defined in the two previous sections, in order to achieve a compromise between static stiffness and dynamic damping. The strategy consists in maintaining the 20 • C temperature where for most cases the temperature is low. Besides, if for the majority, dynamic modes and static included, the temperature is at 77 • C, then the SMP core is heated. Obviously, the temperature field depends on the threshold τ (until then at 0.3). A high value of τ corresponds to higher weight on the static rigidity in the decision, while low values of τ will lead to more damping with lower static rigidity. The discrete temperature fields estimated for τ =0.2, 0.22, 0.25 and 0.3 are shown in Fig. 8. For τ = 0.2, high damping is expected, a large part of the structure has to be heated at 77

RESULTS: STIFFNESS/DAMPING COMPROMISE

Direct harmonic analyses are performed for the four temperature fields determined previously on the [1-3,500] Hz frequency range. The Fig. 9 presents the results for the two extreme temperatures 20 • C and 77 • C and for the four discrete temperature fields associated to each selected threshold. The stiffness/damping compromise obtained is evident: according to the selected threshold, the beam is more or less rigid and more or less damped.

The 0.2 threshold seems quite interesting, indeed the dynamic response is still very damped while the static response is largely improved compare to the 77 • C response. The user can make a choice concerning which is the most important for the considered application, stiffness or damping, and adapt the temperature field accordingly. Moreover it is possible to adjust the temperature fields by adding importance criteria on static state or on one mode or another, to control the frequency response.

CONCLUSIONS

This study highlights a new approach to control the stiffness and the damping thanks to a temperature field in a viscoelastic sandwich core. The methodology which is proposed in this paper is applicable to any viscoelastic material exhibiting frequency -and temperature -dependent high damping properties. The tuning of the stiffness/damping compromise is performed through the elastic strain energy density value. By adjusting the temperature field, the dynamic and static responses are updated and an optimal solution can be found according to the expectations. The structure proposed is adaptive, providing a way to change the composite properties in real time, by playing on the temperature, in order to adapt the structure to its frequency environment. This device opens up news possibilities on the damping and stiffness adaptive control.

Figure 1 .

 1 Figure 1. (a) Passive damping, (b) active damping, (c) thermal adaptive structure proposed in the paper.

Figure 2 .

 2 Figure 2. Sandwich structure used in numerical simulation and in experimental tests.

Figure 3 .

 3 Figure 3. Mean square velocity with respect to z direction, computed on the composite sandwich surface at 20 • C, 50 • C, 80 • C and 130 • C. 18

Figure 4 .

 4 Figure 4. Beam sandwich structure used in numerical simulation.

Figure 5 .

 5 Figure 5. Spatial mean displacement with respect to z direction, computed on the sandwich at 20 • C and 77 • C.

Figure 6 .

 6 Figure 6. Inverse of the static elastic strain energy density and the associated temperature field with a selected threshold to 0.3.

Figure 7 .

 7 Figure 7. Elastic strain energy density for the fourth first modes, and the associated temperature field with a selected threshold to 0.3.

Figure 8 .

 8 Figure 8. Temperature fields in the beam for different thresholds on the elastic strain energy densities.

2 Figure 9 .

 29 Figure 9. Spatial mean displacement with respect to z direction, estimated on the beam at 20 • C, 77 • C and for different thresholds on the elastic strain energy densities.

  • C. When τ increases, the number of cold zones is growing and the static stiffness is increasing. These temperature fields are used in the next section to check the efficiency of the proposed strategy.

		x 10 -3												x 10 -3
	3												3
	1 2	77 20	77		20 77 20			77		20	77	20	1 2	77 20	77	20 77 20 77 20 7720	77	20	77	20	77 20 77 20
	0												0
		0	0.02	0.04	0.06	0.08	0.1	0.12	0.14	0.16	0.18	0.2		0	0.02	0.04	0.06	0.08	0.1	0.12	0.14	0.16	0.18	0.2
		-3												-3
		x 10												x 10
	3												3
	2												2
		20	77	20 77 20 77 20 7720	77	20 77 20 77	20				20	77	20	77		20		77		20
	1												1
	0												0
		0	0.02	0.04	0.06	0.08	0.1	0.12	0.14	0.16	0.18	0.2		0	0.02	0.04	0.06	0.08	0.1	0.12	0.14	0.16	0.18	0.2

ACKNOWLEDGMENTS

This work has been performed in cooperation with the Labex ACTION program (ANR-11-LABX-0001-01). The authors would like to acknowledge the "Fond Interministériel Unique" that funds the project CLIMA, and the "Bourse L'Oral-UNESCO Pour les Femmes et la Science" that encourages women to make sciences.