Jérémy Parriaux
email: jeremy.parriaux@esstin.uhp-nancy.fr

Gilles Millérioux
email: gilles.millerioux@esstin.uhp-nancy.fr

A constructive approach for the design of self-synchronizing dynamical systems: an application to communications

Keywords: switched systems, synchronization, communication, flatness, invertibility

This paper addresses the problem of self-synchronization of dynamical systems in a so-called master-slave configuration. The study is motivated by potential cryptographic applications. It is shown that the notion of flatness is central for guaranteeing self-synchronization and that the concept of transmissions zeros plays also an important role. We motivate the fact that switched linear systems have great interest in this context.

INTRODUCTION

Synchronization of dynamical systems is an important purpose in many fields like biology, mechanics, communications. Synchronization means coordinated behavior of different interconnected entities involved in an overall system. Many different definitions and related configurations, in terms of coupling, can be investigated. An exhaustive and interesting overview can be found in [START_REF] Blekhman | On self-synchronization and controlled synchronization[END_REF]. A special kind of synchronization is the selfsynchronization. By self-synchronization, it is meant a coordinated behavior which is achieved without any external control.

The configuration under consideration in this paper involves two parties: a so-called master system which forces, through a unidirectional coupling, a second system called slave. The configuration is borrowed from the field of communications and more specifically secure transmissions. In this context, cryptography plays a central role. It is the discipline which is mainly intended to protect information and to guarantee confidential exchanges through public channels. One of the cryptographic methods obeys the following principle. The transmitter, called the cipher, delivers a complex sequence (theoretically indistinguishable from a uniformly random sequence) used to conceal information. The information to be kept secret is, in some sense, "mixed" with the complex sequence so that the resulting sequence called cryptogram, which is conveyed to the receiver, cannot be understood by any unauthorized party. For proper information recovery, the receiver, called the decipher, must deliver the same complex sequence synchronized with the cipher. It is typically a master-slave configuration with unidirectional coupling.

The master is nothing but the cipher, the slave is nothing but the decipher. The coupling is achieved through the cryptograms. Some communication setups require that the synchronization must be guaranteed without any external control. In other words, self-synchronization must be achieved. Such a requirement can be motivated by the fact that, for instance, insertion of synchronization flags in the transmitted packets is forbidden for throughput purposes.

From the 90's, many "scrambling" methods resorting to synchronized dynamical systems have been proposed. A recent overview can be found in [START_REF] Banerjee | Chaos Synchronization and Cryptography for Secure Communications: Applications for Encryption[END_REF]. However, their well admitted poor efficiency regarding the security can be explained by the fact that they were disconnected from standard ciphering methods. Recently, in [START_REF] Millérioux | A connection between chaotic and conventional cryptography[END_REF] and [START_REF] Millérioux | Self-synchronizing stream ciphers and dynamical systems: state of the art and open issues[END_REF], the connection has been made and shows that using dynamical systems in a master-slave configuration makes sense from a cryptographic point of view under specific conditions. In particular, whenever the involved dynamical systems are flat, the communication scheme is structurally equivalent to a so-called self-synchronizing stream cipher. So far, the study reduced to analysis, no efficient constructive approach for the design was proposed. The aim of the present work is precisely to handle this problem. Discretetime switched linear systems are specifically addressed because they correspond to the so-called Maiorana McFarland construction which has proved to produce functions that have many interesting cryptographic properties (see [START_REF] Carlet | Boolean Models and Methods in Mathematics[END_REF]).

The outline of this paper is the following. In Section 2, strict necessary background on cryptography is provided. A special emphasis is put on the role of selfsynchronization in this context and a formal definition of finite-time self-synchronization is given. In Section 3, the design of admissible master-slave configurations, described by piecewise linear systems, achieving finite-time self-synchronization is detailed. A constructive approach for guaranteeing the self-synchronization is suggested. It is mainly based on the notion of nilpotent semigroups. A connection between the issue of guaranteeing selfsynchronization and the concept of flatness is brought out. Further considerations for the design are developed in Section 4 where it is shown that the concept of transmission zeros of a dynamical system play an important role. Finally, Section 5 is devoted to an example.

Notation 1 n stands for the identity matrix of dimension n, 0 stands for the zero matrix of appropriate dimension regarding the situation. We denote by {z} k2 k1 the sequence {z k1 , . . . , z k2 } when the initial and final times k 1 and k 2 are defined, otherwise the sequence is merely denoted {z}. A general encryption mechanism, also called cryptosystem or cipher, is illustrated in Fig. 1. We are given an alphabet A, that is, a finite set of basic elements named symbols. On the transmitter part, a plaintext (also called information or message) {u} ∈ U (U is called the message space) consisting of a string of symbols u k ∈ A is encrypted according to an encryption function e which depends on the key k e ∈ K (K is called the key space). The resulting ciphertext {y} ∈ C (C is called the ciphertext space), a string of symbols y k from an alphabet B usually identical to A, is conveyed through a public channel to the receiver. At the receiver side, the ciphertext y k is decrypted according to a decryption function d which depends on the key k d ∈ K. For a prescribed k e , the function e must be invertible. Cryptography distinguishes asymmetric and symmetric ciphers. Asymmetric cryptography is largely based upon computationally very demanding mathematical problems, for instance, integer factorization into primes. It is not discussed in this paper. In symmetric encryption, both keys are identical, that is k d = k e . That explains the terminology "symmetric". Next we describe a special class of symmetric ciphers, under consideration hereafter, called stream ciphers, for which synchronization is a central issue. It is shown that the problem can be tackled with efficiency through the control theory point of view.

CRYPTOGRAPHY AND SYNCHRONIZATION

For stream ciphers, the key k e and k d are replaced by a time-varying sequence called running key or keystream. They are denoted {x} (with samples x k ) at the transmitter part and { x} (with samples x k ) at the receiver part. As a result, stream ciphers require keystream generators. The keystreams {x} and { x} are produced by deterministic dynamical systems and must be synchronized. The secret key k e is some suitable selected parameters of the dynamical systems, the selection being based on security considerations. As mentioned in the introduction, some applications require that the synchronization is guaranteed without any external control that is, self-synchronization must be achieved. In such a case, the stream ciphers must have a special architecture and they are called Self-Synchronizing Stream Ciphers. An overview on this class of ciphers can be found in [START_REF] Millérioux | Self-synchronizing stream ciphers and dynamical systems: state of the art and open issues[END_REF], [START_REF] Daemen | The self-synchronizing stream cipher moustique[END_REF].

Self-synchronization and ciphering

Self-Synchronizing Stream Ciphers (written hereafter SSSC for short) admit at the transmitter and receiver ends the respective equations:

x k = g k e (y k-K , . . . , y k-1 ) y k = e(x k , u k ) (1) x k = g k e (y k-K , . . . , y k-1 ) u k = d( x k , y k ) (2)
g k e is the function that generates the keystreams {x} and { x}. It depends on K past values of y k .

The ciphertext y k is worked out through an encryption function e which must be invertible for any prescribed x k . The decryption is performed through a function d depending on the ciphertext y k and on the running key x k of the receiver. Such a function must obey the rule:

u k = d( x k , y k ) = u k if x k = x k
(3) According to (3), the synchronization of the keystreams {x} and { x} generated respectively at the transmitter and receiver sides is a condition for proper decryption. Since the function g k e is identical at the transmitter and receiver sides and shares the same arguments, namely the past ciphertexts y k-i (i = 1, . . . , K), it is clear that the generators synchronize automatically after a finite transient time of length K. This kind of self-synchronization is called finite time self-synchronization.

Actually, the model ( 1)-( 2) of an SSSC is a conceptual model, called canonical representation, that can correspond to different architectures. It turns out that resorting to dynamical systems instead of implementing directly the canonical function g k e is much more relevant for two major reasons (some detailed motivations can be found in [START_REF] Maurer | New approaches to the design of selfsynchronizing stream cipher. Advance in Cryptography[END_REF] and [START_REF] Daemen | Cipher and Hash function design, strategies based on linear and differential cryptanalysis[END_REF]). First, implementation of a complex function in a recursive way is, in general, much more computationally efficient than implementing the function itself. Secondly, the canonical representation (1)-(2) assumes that the synchronization delay K is bounded. This assumption limits the complexity of the ciphering which can be represented as a memoryless function. This requirement is not mandatory in practice, and it is acceptable that the synchronization delay is not a constant value but a random variable with a probability law that reaches one as time reaches infinity. In this case, self-synchronization is said to be statistical. Statistical To sum up, the problem can be formulated as follows. We are given a setup with two parts (Fig. 2). The first part consists of a dynamical system C, with input u k (playing the role of the plaintext), output y k (playing the role of the ciphertext) and state vector x k (playing the role of the keystream).

C x k+1 = f (x k , u k ) y k = h(x k , u k ) (4)
The output y k ensures a unidirectional coupling with the second part, the dynamical system D with state vector x k . It acts as an input for D.

D x k+1 = f ( x k , y k ) u k = h( x k , y k ) (5)
u k is the output of D. In a cryptographic context, it acts as the recovered information and must be equal to u k whenever x k = x k .

Definition 1. (Finite time self-synchronization). The unidirectional coupled system C-D is finite time self-synchronizing if, for all admissible input sequences, ∃K ∈ ∀x 0 , x 0 , ∀k ≥ K, x k = x k (6)

A delay r ∈ N can be allowed. If so, Equation (6) turns into ∃K ∈ N, ∀x 0 , x 0 , ∀k ≥ K, x k = x k+r (7)

Finally, the issue to be investigated is the following. How to design a master-slave setup C-D so that

• self-synchronization (6) (possibly ( 7)) can be guaranteed

• proper input recovery u k = u k is ensured whenever self-synchronization is achieved
It is the purpose of the next sections. Actually, as motivated in the introduction, we concentrate on the special class of switched linear systems.

FINITE TIME SELF-SYNCHRONIZATION AND SWITCHED SYSTEMS

The equations of the set-up read at the transmitter part

C x k+1 = A σ(k) x k + B σ(k) u k y k = C σ(k) x k + D σ(k) u k (8)
and at the receiver part

D x k+1 = A σ(k) x k + B σ(k) y k u k = C σ(k) x k + D σ(k) y k (9) with u k , u k ∈ F, y k ∈ F and x k , x k ∈ F n
where F is a field. In digital transmissions, F is a finite field of cardinality p q with p a prime and q a positive integer. When q = 1, all the operations, namely, addition, subtraction, multiplication and inversion are still defined like in the field of real numbers except that the results are computed modulo p.

The switching function σ is defined as

σ : k ∈ N → j = σ(k) ∈ {1, . . . , J} = J
At a given time k, the index j corresponds to the mode of the system given by the switching function σ. J is the number of modes. All the matrices, namely

A σ(k) ∈ F n×n , B σ(k) ∈ F n×1 , C σ(k) ∈ F 1×n and D σ(k) ∈ F belong to the respective finite sets (A j ) 1≤j≤J , (B j ) 1≤j≤J , (C j ) 1≤j≤J
and (D j ) 1≤j≤J . The switching function must depend on the output y k . The motivation of such a dependence lies in that the switching rule is the same for both systems C and D and must be self-synchronizing. Thus, it must depend on shared variables and so on the output y k or a finite sequence of delayed outputs. It is worth pointing out that the writing σ(k) is somehow abusive since the dynamical system is not a time-varying system. We denote by {v} the sequence of modes {v} = {σ(0), σ(1), . . .}. If the sequence has a finite length K, it is an element of the set denoted J K .

In the following, we derive conditions for guaranteeing selfsynchronization of the master-slave set-up C-D (8)-( 9) and propose constructive approaches for achieving finitetime self-synchronization.

General conditions

Theorem 1. The set-up (8)-( 9) is finite-time self-synchronizing whenever the three following conditions are fulfilled:

• ∀j ∈ J , D j = 0 (10)

• ∃K ∈ N, ∀x 0 , x 0 , ∀{v} ∈ J K , K-1 i=0 A vi = 0 (11) • Given the matrices (A j , B j , C j , D j ) of D, the system C reads      x k+1 = A σ(k) -B σ(k) (D σ(k) ) -1 C σ(k) x k +B σ(k) (D σ(k) ) -1 u k y k = -(D σ(k) ) -1 C σ(k) x k + (D σ(k) ) -1 u k ( 12
) Proof 1. Since D j = 0 for any j ∈ J (Condition (10)), the input u k can be derived from the output equation of C and reads

u k = D σ(k) y k + C σ(k) x k (13) Thus, u k -u k = C σ(k) x k + D σ(k) y k -D σ(k) y k -C σ(k) x k = C σ(k) ( x k -x k )
Let the reconstruction error be k = x k -x k . Then, from ( 9) and ( 12)

k+1 = A σ(k) x k + B σ(k) y k -(A σ(k) -B σ(k) (D σ(k) ) -1 C σ(k) )x k -B σ(k) (D σ(k) ) -1 u k = A σ(k) k -B σ(k) (y k -(D σ(k) ) -1 u k ) -B σ(k) (D σ(k) ) -1 u k + B σ(k) y k = A σ(k) k (14)
After iterating ( 14) K times and taking into account (11), k = 0 or equivalently x k = x k for any k ≥ K. Hence, according to Definition 1, the set-up ( 8)-( 9) is finite-time self-synchronizing. That completes the proof.

No constraint is imposed on B j and C j . Condition (11) means that regardless of the order of multiplication of the matrices A j , and so for any mode sequences, the product is zero after a finite number K of iterations. K is the delay of synchronization.

Remark 1. The condition D j = 0 for any j ∈ J means that the relative degree of the systems C and D is zero.

Remark 2. The system ( 8) is a right inverse for the system (9). Indeed, for any identical initial conditions x 0 = x 0 and for any identical mode sequence {v}, the system ( 8) drives ( 9) such that ∀k ≥ 0, u k = u k .

Theorem 1 does not provide a constructive solution for the selection of appropriate matrices A j which must fulfill the constraint ( 11). The purpose of the next paragraph is to obtain an equivalent constructive condition. It is based on the notion of nilpotent semigroups.

Constructive approach

Let us first recall two definitions:

Definition 2. (Semigroup). A semigroup S is a set together with an associative internal law. It is said to be finite if S has a finite number of elements.

Definition 3. (Nilpotent semigroup). A semigroup S is said to be nilpotent if it is such that any product of a finite number t ∈ N * of its elements (possibly the same element) is always 0. The smallest integer t is called the class of nilpotency of S.

Proposition 1. In order for (11) to be fulfilled, the set of dynamical matrices {A j , j ∈ J must generate a nilpotent semigroup. The delay of synchronization K equals the class of nilpotency of this semigroup.

A theorem, useful for the construction of semigroups with a given class of nilpotency is stated in the book [START_REF] Radjavi | Simultaneous Triangularization[END_REF] and recalled below.

Theorem 2. (Levitsky's theorem). Any semigroup of nilpotent matrices can be triangularized.

In other words, there is a common basis in which all the matrices of the semigroup are upper triangular with zeros on the diagonal.

Remark 3. The product of t nilpotent matrices which commute pairwise is 0 but the product of t nilpotent matrices is not, in general, nilpotent. Indeed, we observe that 0 1 0 0 0 0 1 0 = 1 0 0 0 . Theorem 2 provides a generalization of this special case, should each matrix be nilpotent is only a necessary condition.

Hence, based on Levitsky's theorem, the construction of the family (A j ) 1≤j≤J which fulfills (11) follows three successive steps

• choose an invertible matrix T ∈ F n×n • choose a set of J upper triangular matrices Ā j with zero on the diagonal • for all j ∈ J , compute

A j = T -1 Ā j T
The matrix T may possibly be the identity matrix.

Remark 4. Because of Levitzky's theorem, the consideration of a semigroup of n-dimensional matrices is equivalent to the consideration of the corresponding set of upper triangular matrices. And yet, for triangular matrices, it is clear that the nilpotency class is at most n. As a result, the delay of synchronization K is upper bounded by n.

Connection with flatness

Flatness is an important concept in control theory. It was introduced in [START_REF] Fliess | Flatness and defect of non-linear systems: introductory theory and examples[END_REF] and a deep insight can be found in the book Sira-Ramirez and Agrawal [2004]. In this section, we show that the constructive approach proposed for designing a finite-time self-synchronizing master-slave system amounts to designing a flat system C with flat output y k .

Definition 4. (Flat dynamical system). A system with input u k and state vector x k is said to be flat if there is a set of independent variables y k , referred to as flat output, such that all the system variables can be expressed as a function of the flat output and a finite number of its backward and/or forward iterates. In particular, there exist two functions F and G such that

x k = F(y k+t1 , . . . , y k+t2 ) u k = G(y k+t 1 , . . . , y k+t 2 )
where t 1 , t 2 , t 1 , t 2 ∈ Z.

Proposition 2. The system C resulting from the conditions (10)-( 11)-( 12) is flat with flat output y k Proof 2. The state of the switched system (9) can be written, at time k + K

x k+K = K-1 i=0 A σ(k+K-1-i) x k + K-1 i=0   K-1 j=i+1 A σ(k+K-j)   B σ(k+i) y k+i
Therefore, if (11) holds, any state at time k ≥ 0 reads:

x k+K = K-1 i=0   K-1 j=i+1 A σ(k+K-j)   B σ(k+i) y k+i (15)
And yet, k = 0 or equivalently x k = x k for any k ≥ K. Hence, after a shift of K, one obtains

x k = x k = K-1 i=0   K-1 j=i+1 A σ(k-j)   B σ(k-K+i) y k-K+i (16)
which gives the function F.

On the other hand, since D j = 0 for any j ∈ J , the input u k reads like (13). Substituting the expression ( 16) of x k into (13) gives the function G. That completes the proof.

Relation ( 16), and in the general case, the function F, gives explicitly the function g k e of ( 1) and ( 2). As a result, and as pointed out in Section 2.2, there is an equivalence between the recursive part of both equations ( 4) and ( 5) and the function F. The equivalence applies under flatness conditions. In the special case of switched linear systems, ( 8) and ( 9) can be equivalently rewritten into the respective canonical forms (1) and ( 2)

       x k = K-1 i=0   K-1 j=i+1 A σ(k-j)   B σ(k-K+i) y k-K+i y k = C σ(k) x k + D σ(k) u k (17)        x k = K-1 i=0   K-1 j=i+1 A σ(k-j)   B σ(k-K+i) y k-K+i u k = C σ(k) x k + D σ(k) y k (18)
It is worth pointing out that, from a computational point of view, the recursive form ( 8)-( 9) is more relevant than ( 17)-( 18)

TRANSMISSION ZEROS AND SURJECTIVITY

For cryptographic purposes (basically a consideration regarding the entropy of sequences), it is relevant that the maps x k → A j x k j ∈ J are surjective. In other words, we want to guarantee that ∀j ∈ J , rank(A j ) = n (19) The problem lies in that, according to Theorem 1, the matrices (A j , B j , C j , D j ) of the system C are not designed directly but are derived from (A j , B j , C j , D j ) of D. Hence, we must find out a condition on the matrices (A j , B j , C j , D j ) so that ( 19) is ensured. It turns out that the notion of transmission zeros are relevant to this end.

A definition of transmission zeros can be found for example in [START_REF] Schrader | Research on system zeros: a survey[END_REF]. It is recalled below and particularized for a SISO system. Definition 5. Let us consider a SISO linear system with state space realization (A, B, C, D). The transmission zeros are the complex numbers {s i } which satisfy rank

A -s i 1 n B C D < n + 1 ( 20
)
where it is recalled that 1 n stands for the identity matrix of dimension n.

Before proceeding further, we must introduce some notation. Consider the matrix T and the corresponding matrices A j = T -1 Ā j T derived from Ā j (j ∈ J ) as explained in Section 3.2 devoted to the constructive approach.

Let us write Ā j as

Ā j =        0 a 1 j 0 a 2 j A * j . . . . . . 0 0 a n-1 j • • • 0        (21)
where A * j denotes the coefficients above the n -1 diagonal entries a m j (m = 1, . . . , n -1) located above the zero diagonal. Let have

T B j = [b 1 j • • • b n j ] T (22)
b m j stands for the m th component of the column vector T B j .

C j T -1 = [c 1 j • • • c n j ] (23)
c m j stands for the m th component of the row vector C j T -1 .

Proposition 3. The surjectivity of each map x k → A j x k (j ∈ {1, . . . , J}) of C is guaranteed whenever

c 1 j b n j n-1 m=1 a m j = 0 (24)
Proof 3. According to Remark 2, C is a right inverse for D. Furthermore, let us recall that (see Remark 1) the relative degree of C and D is zero. We conclude that each realization (A j , B j , C j , D j ) (j ∈ J ) of D has n transmission zeros s i and the s i 's are nothing but the n eigenvalues λ i of A j of C. They are the roots of

Ψ j (s) = det R = 0 with R = A j -s1 n B j C j D j (25)
R is often called the Rosenbrock's matrix. Ψ j (s) is a polynomial, its constant monomial is Ψ j (0) and corresponds to the product n i=1 of the roots of Ψ j (s) and so corresponds to the product n i=1 λ i of the eigenvalues of A j of C. Hence, surjectivity of x k → A j x k j ∈ J is guaranteed whenever Ψ j (0) = 0.

The following equalities apply

Ψ j (0) = det A j B j C j D j = det T -1 Ā j T B j C j D j = det T -1 0 0 1 Ā j T B j C j T -1 D j T 0 0 1 = det Ā j T B j C j T -1 D j (26)
Consider a partitioned matrix with four sub-blocks E, F , G, H of compatible dimensions. We recall a result concerning its determinant.

det E F G H = det(H) • det(E -F H -1 G)
Taking into account the special structure (21) of Ā j , ( 22) and ( 23), it turns out that basic manipulations yield

Ψ j (0) = det Ā j T B j C j T -1 D j = c 1 j b n j n-1 m=1 a m j (27)
That completes the proof.

EXAMPLE

This section gives an example that illustrates the construction of a finite-time self-synchronizing setup. We propose to design a finite time self-synchronizing system of dimension n = 3 and with J = 3 modes. We consider matrices defined over the finite field F = Z/7Z. This means that the only coefficients allowed for the matrices are elements in the set {0, . . . , 6} and that the operations of additions and multiplications are performed modulo 7.

The design starts with the setting of the matrices A j , B j , C j , D j which must fulfill the three conditions of Theorem 1, the condition (11) being replaced by the constructive approach provided in Section 3.2. We add the condition (24) on surjectivity. First, for simplicity, we choose D j = 1 for any j ∈ {1, 2, 3}.

Secondly, we choose a set of three 3-dimensional matrices Ā j in the form of strict upper triangular matrices and with non zero entries located above the diagonal in order to guarantee the surjectivity. After the setting is completed, a sequence {u} is applied to C. As expected, the self-synchronization is achieved after a finite transient time, so does the recovery of the sequence of inputs (see Figure 3). The transient time before selfsynchronization is of length K = 3 since the class of nilpotency t of the set (A 1 , A 2 , A 3 ) equals 3.

  Fig. 1. General encryption mechanism

Fig. 2 .

 2 Fig. 2. Dynamical system-based cryptosystemsself-synchronization is more general than the finite-time one. Its interest lies in a broader choice of candidate dynamical systems. The resulting flexibility is important in view of matching additional constraints, besides the self-synchronization, in particular regarding the security of the communication setup. Statistical self-synchronization is not addressed in the present paper, it is detailed for the Boolean case in[START_REF] Parriaux | Synchronization of boolean dynamical systems: a spectral characterization[END_REF].

  of basis A j = T -1 Ā j T , we get that choose arbitrary matrices B j and C j except the fact that the first entry c 1 j of C j T -1 and last entry b n j of T B j are not zero to fulfill the surjectivity condition (24). , C 2 = (1 5 6) , C 3 = (4 6 6)D 1 = D 2 = D 3 = 1

Fig. 3 .

 3 Fig. 3. Time evolution of {u} and { u} of the setup C-D

This work was supported in part by the Institut


Laurent Arcese
email: laurent.arcese@bourges.univ-orleans.fr

Matthieu Fruchard
email: matthieu.fruchard@bourges.univ-orleans.fr

Antoine Ferreira
email: antoine.ferreira@ensi-bourges.fr

Endovascular Magnetically-Guided Robots: Navigation Modeling and Optimization

Keywords: Endovascular navigation, magnetic steering, nonlinear modeling, optimal trajectory, nonlinear controller and observer

This paper deals with the benefits of using a nonlinear model-based approach for controlling magnetically guided therapeutic microrobots in the cardiovascular system. Such robots used for minimally invasive interventions consist of a polymer binded aggregate of nanosized ferromagnetic particles functionalized by drug-conjugated micelles. The proposed modeling addresses wall effects (blood velocity in minor and major vessels' bifurcations, pulsatile blood flow and vessel walls, and effect of robot-to-vessel diameter ratio), wall interactions (contact, van der Waals, electrostatic and steric forces), non-Newtonian behaviour of blood and different driving designs as well. Despite nonlinear and thorough, the resulting model can both be exploited to improve the targeting ability and be controlled in closed-loop using nonlinear control theory tools. In particular, we infer from the model an optimization of both the designs and the reference trajectory to minimize the control efforts. Efficiency and robustness to noise and model parameter's uncertainties are then illustrated through simulations results for a bead pulled robot of radius 250µm in a small artery.

I. INTRODUCTION

Emerging real-time high spatial resolution medical imagers, using either computed tomography [START_REF] Achenbach | Contrast-enhanced coronary artery visualization by dual-source computed tomographyinitial experience[END_REF], optical coherence tomography or Magnetic Resonance Imaging (MRI) devices [START_REF] Uecker | Real-time mri at resolution of 20 ms[END_REF], [START_REF] Zhang | Real-time magnetic resonance imaging of temporomandibular joint dynamics[END_REF], are promising for the development of automated medical interventions. In particular, it opens a wide field in minimally invasive diagnostic, therapy or in-vivo biosensing. Deflected ferromagnetic catheters [START_REF] Mcneil | Functional design features and initial performance characteristics of a magnetic-implant guidance system for stereotactic neurosurgery[END_REF] or unthetered capsule [START_REF] Carpi | Magnetic manoeuvring of endoscopic capsules by means of a robotic navigation system[END_REF] are already used for medical interventions at the millimeter scale. Recent developments in nanotechnology have allowed the exploration of robots capable of propulsion at the micro and nano scale. One of the primary challenges are the locomotion and power at these scales. Since the smallest adequate motors are some hundreds micrometers long, embedded actuation does not seem to be a suitable choice for artificial robots swimming in the blood flow. Furthermore, traditional onboard power supplies, i.e batteries, cannot be scaled at these dimensions.

Two main approaches have been proposed to solve this issue using harvesting energy from the biological environment and an externally energy field. The first approach employs chemical locomotion to harvest energy from chemicals in the environment such as H 2 O 2 [START_REF] Sitti | Miniature devices: Voyage of the microrobots[END_REF]. Magnetotactic bacteria actuated thanks to embedded ferromagnetic material has also been demonstrated [START_REF] Martel | Flagellated magnetotactic bacteria as controlled mri-trackable propulsion and steering systems for medical nanorobots operating in the human microvasculature[END_REF]. To date the most systematically approach uses wireless power transmission through magnetic fields [START_REF] Mcneil | Functional design features and initial performance characteristics of a magnetic-implant guidance system for stereotactic neurosurgery[END_REF], [START_REF] Honda | Micro swimming mechanisms propelled by external magnetic fields[END_REF]. Different designs have been proposed in the literature for such microrobots with distinguishing physical features [START_REF] Nelson | Microrobots for minimally invasive medecine[END_REF]. A first propulsion scheme is biologically inspired using either Helical Tail (HT) [START_REF] Edd | Biomimetic propulsion for a swimming surgical microrobot[END_REF], [START_REF] Zhang | Artificial bacterial flagella for micromanipulation[END_REF], beating flagella [START_REF] Dreyfus | Microscopic artificial swimmers[END_REF], or Elastic Flagella (EF) [START_REF] Wiggins | Trapping and wiggling: Elastohydrodynamics of driven microfilaments[END_REF], [START_REF] Evans | Propulsion by passive filaments and active flagella near boundaries[END_REF]; a second one is referred as Bead Pulling (BP) [START_REF] Choi | Twodimensional locomotion of a microrobot with a novel stationary electromagnetic actuation system[END_REF], [START_REF] Floyd | Two-dimensional contact and noncontact micromanipulation in liquid using untethered mobile magnetic microrobot[END_REF]. Actually, the most advanced design relies on bead pulling [START_REF] Mathieu | Method of propulsion of a ferromagnetic core in the cardiovascular system through magnetic gradients generated by an mri system[END_REF] since in-vivo experiments were conducted in the carotid artery of a living swine [START_REF] Martel | Automatic navigation of an untethered device in the artery of a living animal using a conventional clinical magnetic resonance imaging system[END_REF].

All these contributions emphasize the difficulties of controlling magnetic microrobots in fluids through magnetic fields actuation. Achieving motion control of magnetic microrobot necessitates nonlinear dynamics modeling and analysis [START_REF] Abbott | Modeling magnetic torque and force for controlled manipulation of soft-magnetic bodies[END_REF], [START_REF] Li | Dynamics modeling and analysis of a swimming microrobot for controlled drug delivery[END_REF], [START_REF] Pawashe | Modeling and experimental characterization of an untethered magnetic micro-robot[END_REF]. Some models are very precise and thorough but rely on finite or discrete elements methods [START_REF] Vartholomeos | Simulation platform for selfassembly structures in mri-based nanorobotic drug delivery systems[END_REF]: while these studies are useful to predict the microrobot's trajectory, they are awkward for parametric analysis and their high computational costs are not compatible with on-line control purposes and experiments. At the opposite, most of the literature dedicated to the control of magnetically-guided microrobots neglects wall effects and interactions. These studies mainly address the control of linearized models. Thus, stability and robustness issues, pointed out by the authors in [START_REF] Tamaz | Real-time mri-based control of a ferromagnetic core for endovascular navigation[END_REF], [START_REF] Choi | Position stabilization of microrobot using pressure signal in pulsating flow of blood vessel[END_REF], could have been expected.

On top of dealing with the nonlinearities of the model, control approaches have to solve some actuators sizing problems [START_REF] Vartholomeos | Simulation platform for selfassembly structures in mri-based nanorobotic drug delivery systems[END_REF], [START_REF] Jun | Nanoscaling laws of magnetic nanoparticles and their applicabilities in biomedical sciences[END_REF] due to both limitations on actuators and the increasing prevalence of blood flow drag over magnetic force as the robot navigates in smaller vessels [START_REF] Abbott | How should microrobots swim?[END_REF]. In [START_REF] Mathieu | Method of propulsion of a ferromagnetic core in the cardiovascular system through magnetic gradients generated by an mri system[END_REF] and [START_REF] Vartholomeos | Simulation platform for selfassembly structures in mri-based nanorobotic drug delivery systems[END_REF], the authors show that an optimal robot-to-vessel diameter ratio should be respected so as to minimize the drag force prevalence. Hence, designing the microrobot as a degradable binded aggregate of magnetic nanoparticles can be a technological way to navigate in smaller vessels.

The present paper addresses the benefits of using a nonlinear model-based approach for designing and controlling a magnetically-guided functionalized microrobot in the blood flow, especially at vessel bifurcations (see Figure 1). The contribution is threefold. First, in Section II, we give a precise forces balance, which includes wall effects (bifurcations, pulsatile blood flow), wall interactions (contact, van der Waals, electrostatic, and steric forces induced by the robot's functionalization) and non-Newtonian behaviour of blood. In particular, we address the modeling of bifurcation effects (minor and major ones) using analytical approximations of Computational Fluid Dynamics (CFD) simulations. We also address the design of nonlinear observer and controller that guarantee the stability of the tracking along any navigation trajectory. Section III-A is then devoted to the exploitation of the proposed nonlinear modeling to optimize both the design and the reference trajectory that the microrobot should track in the macro and micro-vasculature. In particular, we show that the optimal size for different swimming methods (BP, EF and HT designs) should not be as big as expected in [START_REF] Mathieu | Method of propulsion of a ferromagnetic core in the cardiovascular system through magnetic gradients generated by an mri system[END_REF] due to the influence of wall interactions. The navigation can also be optimized so as to minimize the control efforts. These two results suggest that embolization hazards can be reduced while improving the targeting of medical interventions. In Section III-C, simulation results illustrate the performances and robustness to noise measurements and parametric errors of the proposed approach in a small artery, with concerns about the limitations of a clinical MRI gradient coils. Finally, we discuss some open issues in Section III-D and conclude in Section IV.

II. METHODOLOGY A. Analytical model

Let consider a spherical ferromagnetic microrobot, of radius r, immersed in blood vessels. The microrobot is governed by several forces including a magnetic motive force F m , the blood hydrodynamic drag force F d , the microrobot's apparent weight W a , the contact force F c , the van der Waals and the electrostatic interaction forces, F vdw and F el , between the microrobot and the vessel walls and finally the steric interaction force F s between the robot's polymer chains and the endothelial glycocalyx. The Knudsen number is large enough to neglect the effect of Brownian motion. The translational motion of a robot of mass m is given by:

m d v dt = F m + F d + W a + F c + F vdw + F el + F s ( 1
)
where v is the translational velocity of the robot. 1) Motive force: Three main schemes of magnetically actuated microrobots have been developped (see Figure 2), namely bead pulled, elastic flagellated, and helical tailed robots. The motive force F m for each class of robot is detailed hereinafter.

• Bead pulling (BP) approach relies on the use of magnetic gradient coils for propulsion. The magnetic gradient induces a magnetic motive force F m on the robot (see Figure 2(a)):

F m = τ m V ( M .∇) B (2)
where V and V m are the robot's and the ferromagnetic volumes, τ m = Vm V , and M is the bead's magnetization, ∇ B is the magnetic gradient. Within an MRI device, the orientation of the robot does not change due to the magnetic torque Γ m which tends to align the robot's magnetization M along the MRI uniform magnetic field B 0 = 3T , which also ensures a saturation magnetization M sat .

• Biologically inspired flagellated robots (see Figure 2(b))

have been designed through various concepts, see [START_REF] Nelson | Microrobots for minimally invasive medecine[END_REF] for a survey. In magnetically-driven Elastic Flagellated propulsion (EF), the wave motion of the elastic tail is induced by the magnetic torque applied on the attached magnetic bead:

Γ m = τ m V M × B (3)
With an oscillating magnetic input B of norm B 0 at frequency ω, the post-transient resulting motive torque has an amplitude γ 0 = τ m V M B 0 . For small-amplitude motion (under assumption that γ0L kκ ≪ 1), the elastohydrodynamics are described by an hyperdiffusion equation and the resulting motive force is given by [START_REF] Evans | Propulsion by passive filaments and active flagella near boundaries[END_REF]: . The latter is a decreasing function of S p , valid only for S p ≥ 3.4 [START_REF] Evans | Propulsion by passive filaments and active flagella near boundaries[END_REF]. It should be noticed that more effective concepts have been proposed, e.g. planar wave flagellated robot [START_REF] Behkam | Design methodology for biomimetic propulsion of miniature swimming robots[END_REF], where the increased efficiency requires a distributed actuation.

F m = ωζ ⊥ γ 2 0 L 4 4k 2 κ Z(S p ) < B > < B > (4)
• The rotating magnetic torque (3) can also be exploited to generate a corkscrew effect for a bead attached to an helical tail (HT) illustrated by Figure 2(c). The post-transient motive force is given by [START_REF] Abbott | How should microrobots swim?[END_REF]:

F m = 2πnr[r cos ϕ(ζ ⊥ -ζ )ω Γm Γm - ζ ⊥ sin 2 ϕ+ζ cos 2 ϕ sin ϕ ( v-v f )] (5
) where ϕ and n are the pitch angle and the number of helical turns, v and v f respectively denote the robot's and the fluid's velocity. The viscous drag coefficients for flows respectively normal and parallel to the tail's axis are given by:

ζ ⊥ = 4πη ln 0.36πr κ sin ϕ + 0.5 > ζ = 2πη ln 0.36πr κ sin ϕ (6)
In [START_REF] Carpi | Magnetic manoeuvring of endoscopic capsules by means of a robotic navigation system[END_REF], the first term is a pure driving force resulting from the corkscrew effect, and the second term is a slow down force induced by the drag exerting on the tail. 2) Hydrodynamic drag force: Hydrodynamics effects are classically modeled using the Navier-Stokes equation, whose analytical solution is lacking. A numerical approach implies time-consuming computations, that are not compatible with real-time approaches. Assuming that the fluid is incompressible, the drag force exerted on a sphere is approximated by:

F d = - 1 2 ρ f ( v -v f ) 2 ∞ AC d∞ ( v -v f ) ∞ ( v -v f ) ∞ (7)
where ρ f denotes the fluid density, A is the robot's frontal area, and C d is the R e -dependent drag coefficient; for a sphere in a fluid of infinite extent in laminar flow, we have:

C d∞ = 24 Re∞ + 6 1 + √ Re∞ + 0.4, Re∞ = 2rρ f v -v f ∞ η (8)
• Non-Newtonian fluid. In the case of blood, the fluid's viscosity η depends both on the vessel's diameter d (in 10 -6 m) and on the hematocrit rate h d (see Figure 3) according to the three following empirical relations [START_REF] Pries | Biophysical aspects of blood flow in the microvasculature[END_REF]:

η = η pl d 2 (d -1.1) 2 1 + (η 0.45 -1)d 2 (d -1.1) 2 (1 -h d ) c -1 (1 -0.45) c -1 (9)
with parameters η pl and η 0.45 denoting respectively the plasma viscosity and the relative apparent blood viscosity for a fixed discharge hematocrit of 0.45, given by: η 0.45 = 6 e -0.085d + 3.2 -2.44 e -0.06d 0.645

The shape of the viscosity dependance on hematocrit is:

c = 10 11 d 12 -(0.8 + e -0.075d ) d 12 d 12 + 10

Capillaries Arterioles Small arteries Arteries

Fig. 3. Blood viscosity exhibits non-Newtonian behaviour in smaller vessels.

• Wall effects. Equation ( 7) neglects wall effects, which has to be corrected considering the impact of the vessel walls on the robot's velocity for endovascular applications. This wall effect is commonly expressed as a ratio between the terminal relative velocity of the robot ( v -v f ) and its velocity ( v -v f ) ∞ in an infinite extent of fluid [START_REF] Kehlenbeck | Empirical relationships for the terminal settling velocity of spheres in cylindrical columns[END_REF]:

v -v f ( v -v f ) ∞ = 1 -λ α0 1 + λ λ0 α0 (10)
with diameter ratio λ = 2r/D and D denoting the vessel diameter (in m). Parameters α 0 and λ 0 are functions of Reynolds number, but are commonly set to 1.5 and 0.29, respectively. Thus, equation ( 7) becomes:

F d = - ρ f 2   1 + λ λ 0 α 0 1 -λ α 0 ( v -v f )   2 AC d ( v -v f ) ( v -v f ) (11)
Boundary layer at the wall-fluid interface also results in a parabolic profile of blood flow (see Figure 4). The pulsative blood's velocity is modeled by an affine combination of a time-varying periodic flow with a spatial parabolic profile p p (δ). So as to simplify the analytical expression, with no loss of generality, we only consider the first terms in the time-varying Fourier series of the physiological pulse. In the case of an artery, such an approximation leads to: where n is the unit vector normal to the wall and pointing outside the vessel. Let M 0 , N and M denote respectively the center of gravity of the microrobot and its projections along n onto the vessel wall at equilibrium and onto the microrobot's outer surface. The signed distance δ between the wall and the robot is thus defined by δ = --→ N M . n. It should be noted that setting a f > 1 in ( 12) results in cyclic occurrence of negative values for v f (t), what properly models the diastolic backflow occurring in arteries. Pulsatile flow caused by cardiac pumping also induces a periodic 10% deformation of the vessel diameter D(t) (see Figure 4), synchronized with the blood velocity v f (t):

v f (t) = V f (1 + a f sin (w f t + φ f ))p p (δ) (12)
θ = 0 • θ = 20 • θ = 40 • θ = 60 • θ = 80 • θ = 90 •
D(t) = D (1 + 0.1 sin (w f t + φ f )) (13)
• Bifurcation effects. The velocity profile strongly depends on the geometry of the vessel. Minor bifurcations lead to slight perturbations in the velocity profile and amplitude (see Figure 5): the usual parabolic profile of the blood flow supplies a good approximation. On the contrary, major bifurcation entails major changes (Figure 6) and the parabolic profile in ( 12) is no more valid. On Figure 6(a), before the apex of the bifurcation, the increase of the vessel's cross-sectional area results in skewed and flattened profiles. In the vicinity of the apex region, the velocity at the parent branch centerline promptly cancels out (from θ = 25 • to 30 • ) and the maximum velocity is reached close to the upper wall. After the apex, the maximum velocity progressively shifts to the inner wall (see Figure 6(c)). The flow then gradually returns to a steady-state parabolic profile in the daughter branch. Since these numerical simulations, consistent with the dedicated literature, are not compatible with real-time control constraints, it is necessary to develop analytical models of the velocity profiles that fits CFD results. Even though one can find a polynomial approximation for each angle value θ, there are two severe obstructions to such an approximation: the problem is ill-conditioned and it does not allow for a smooth change of velocity profiles as θ grows. We therefore propose a rational approximation of the blood velocity profile in the form p r (δ) = aδ 2 + bδ + c + dδ+e (δ-p) 2 +q 2 and solve a well-conditioned nonlinear inverse problem for some θ values. Figure 6 shows that the interpolation results in smooth profiles that fit quite well the CFD results. Note that adding another rational fraction to p r (δ) results in a better approximation, at the expense of a more complicated inverse problem. Since such in-vivo measurements are awkward, mathematical models and numerical simulations will become the main investigative tools.

3) Apparent weight: Apparent weight, i.e. combined action of weight and buoyancy, is also acting on the spherical robot:

W a = V (ρ -ρ f ) g (14)
where ρ = τ m ρ m + (1 -τ m )ρ poly with ρ m and ρ poly the magnetic material's and polymer's densities. 4) Contact force: Since the robot moves in low velocity, the duration of the impact occurs within a few milliseconds. Hence, the contact force can be expressed by a modified Hertzian contact law [START_REF] Choi | Low-velocity impact analysis of composite laminates using linearized contact law[END_REF] under the assumption of no friction between the robot and the wall during the impact:

F c = -k|δ| 3/2 H(δ) n : loading (15)
with H the Heaviside step function. Stiffness k is given by : where E c and σ c (respectively E w and σ w ) are the Young's modulus and Poisson's ratio of the microrobot (respectively of the vessel wall).

k = 4 √ r 3 1-σ 2 c Ec + 1-σ 2
The contact force during unloading is different from the one during loading and it is given by:

F c = -F δm δ -δ 0 δ m -δ 0 p H(δ) n : unloading (16)
where F δm and δ m are the maximum contact force and the corresponding maximum deformation of the wall. δ 0 is the permanent wall's deformation from loading/unloading cycle and exponent q belongs to [1.5; 2.5].

5) Van der Waals and electrostatic forces:

The microrobot and the wall also interact each other through van der Waals and electrostatic forces. When the microrobot is close to the wall, van der Waals force is dominant. As the robot moves away from the wall, the van der Waals force rapidly decreases and the electrostatic force prevails.

The van der Waals potential between the spherical robot and the wall is given by [START_REF] Iimura | Simulation of entrainment of agglomerates from plate surfaces by shear flows[END_REF]:

V vdw = - A h 6 1 h + 1 2 + h + ln h 2 + h n ( 17
)
where A h is the Hamaker constant and h = |δ|/r is the normalized distance from the robot to the wall. Then, the attractive van der Waals interaction force is obtained by differenciating [START_REF] Mathieu | Method of propulsion of a ferromagnetic core in the cardiovascular system through magnetic gradients generated by an mri system[END_REF]:

F vdw = -(∇ V vdw )H(-δ) (18)
The electrostatic force between the robot and the wall considered as an uncharged surface is given by [START_REF] Hays | Electrostatic adhesion of non-uniformly charged dielectric spheres[END_REF]:

F el = q 2 4πεε 0 (r + |δ|) 2 H(-δ) n (19)
with q the robot's charge, and ε 0 , ε the vacuum and the relative blood's permittivities. [START_REF]Role of Electrostatics in Adhesion[END_REF] gives a maximum allowable charge q(µC) = 4πr 2 × 30(100r) -0.3 for a sphere of radius r. 6) Steric force: Therapeutic microrobots are generally functionalized via polymer chain grafting, e.g. doxorubicin linked via a polymeric micelle to the robot's outer wall. This polymer brush, characterized by an equilibrium length L 0 , affects the interaction of the microrobot with the endothelial glycocalyx of the vessels walls. For 0.2L 0 < |δ| < 0.9L 0 , the repulsive interaction potential per unit area between the functionalized robot and the endothelium is given by [START_REF] Butt | Steric forces measured with the atomic force microscope at various temperatures[END_REF]:

V s, p/p (z) = - 25 π k b T L 0 τ 3 2 e -2π(z/L0) n ( 20
)
where k b is the Boltzmann constant, T is the absolute temperature, τ is the grafting density of the polymer, z the distance between the two flat surfaces, and index p/p stands for a planeto-plane interaction. The repulsive steric force per unit area between two flat surfaces is then deduced from (20):

f s,p/p (z) = -(∇ V s,p/p ) = -50k b T τ 3 2 e -2π(z/L0) n (21)
Equation ( 21) is integrated using the Derjaguin approximation to infer the overall interaction between a sphere at distance |δ| from a flat surface modeling the vessel wall:

F s (δ) = F s,s/p (δ) = 2πr ∞ |δ| f s,p/p (z)dz = -50rL 0 k b T τ 3 2 e -2π(|δ|/L0) H(-δ) n (22)

B. Nonlinear Controller and Observer Design

Now, let x 1 , x 2 , (x 3 , x 4 ) denote respectively robot's position and velocity along ı axis (respectively along k axis). Assuming that positions x 1 and x 3 can be measured thanks to a realtime medical imager, let y denotes the state's measure. Using expressions of forces given by ( 2), ( 11), ( 14), ( 15), ( 16), ( 18), [START_REF] Abbott | Modeling magnetic torque and force for controlled manipulation of soft-magnetic bodies[END_REF], [START_REF] Vartholomeos | Simulation platform for selfassembly structures in mri-based nanorobotic drug delivery systems[END_REF] and adequate projection of local frame along the geometry of the vessel, equation ( 1) can be written in the canonical form:

(S)    Ẋ = Z Ż = F (X, Z) + U Y = X (23)
where state vectors are given by X = (x 1 , x 3 ) and Z = (x 2 , x 4 ), the input vector U = M ρ (∇B x , ∇B z ) is directly related to magnetic gradients, and F is the function on R 2 that stems from the projection of F tot .

Denoting X ref , Ẋref , Ẍref respectively a reference trajectory, velocity and acceleration, the backstepping control law [START_REF] Kanellakopoulos | A toolkit for nonlinear feedback design[END_REF] for system (S):

U = Ẍref + (k 1 + k 2 ) Ẋref + (1 + k 1 k 2 )X ref -(1 + k 1 k 2 )X -(k 1 + k 2 )Z -F (X, Z) (24)
where k 1 and k 2 are the strictly positive backstepping gains, ensures the Lyapunov stabilization of the microrobot along the desired optimal trajectory (see [START_REF] Arcese | Dynamic behavior investigation for trajectory control of a microrobot in blood vessels[END_REF] for more details). The backstepping control law [START_REF] Choi | Position stabilization of microrobot using pressure signal in pulsating flow of blood vessel[END_REF] requires the knowledge of the robot's velocity in addition to its position, measured by the imager. Since the velocity is not measured, we build an high gain observer to reconstruct it:

(S ′ ) Ẋ = Ẑ + LG x ( X -Y ) Ż = F ( X, Ẑ) + U + L 2 G z ( X -Y ) ( 25
)
where L is the high gain, and matrices G x , G y are conveniently chosen [START_REF] Gauthier | A simple observer for nonlinear systems. application to bioreactors[END_REF].

III. RESULTS

The goal of the present section is to derive guidelines for optimizing the design, trajectory and control of magnetically propelled microrobots towards any bifurcation branch.

A. Model-based optimization

Limitations on the hardwares used to generate the magnetic fields and gradients induce strong restrictions on the available motive force. For BP design, the optimization is based on a maximum magnetic field gradient ||∇ B|| = 80mT.m -1 available on recent clinical MRI devices, with a uniform magnetic field of amplitude B 0 = 3T . For EF and HT designs, we consider a maximum synchronized frequency ω 2π = 30Hz and a rotating magnetic field of amplitude B 0 = 2mT [START_REF] Zhang | Artificial bacterial flagella for micromanipulation[END_REF]. To acquire results that apply for a broad range of blood vessels, III. We perform a model-based optimization to determine both the optimal capabilities of each robot's design and the less demanding reference trajectory so as not to reach actuators' saturations. The total force F tot an endovascular robot is subjected to is the sum of all the forces acting on the robot apart from the motive force, i.e.

F tot = F d + F el + F vdw + F s + W a .
Whatever the selected design (namely BP, EF or HT robot) the same external forces apply: the only difference stands in the source of their motive magnetic force F m (cf Section II-A). It is thus necessary to optimize both the design and the trajectory in order to maximize the ratio Fm Ftot defined by Martel in [START_REF] Mathieu | Method of propulsion of a ferromagnetic core in the cardiovascular system through magnetic gradients generated by an mri system[END_REF]. Let ∆ denotes the distance between the microrobot and the vessel centerline. For each vessel, a family of non-dimensional 3D curves is drawn as a function of λ = 2r/D (for the design optimization) and ∆ (for the optimal trajectory), as shown in Figures 7(a) to 7(l). Physiological parameters' ranges, measured or estimated on patients, are given in Table II. The 3D curves may therefore differ slightly from Figures 7(a) to 7(l) depending on the parameters' values, yet the general conclusions remain valid.

1) Design optimization: BP, EF and HT designs can be optimized with respect to the respective sets of parameters {r}, {r, L, κ, k κ , S p }, and {r, n, θ, κ}, so as to maximize the motive force F m . To compare equitably helical and elastic tail designs, we set a maximal ratio between the tail's length L and the bead's radius r corresponding to n = 10 helical turns for an helical tail with an optimal pitch angle ϕ = 35 • . Without violating the slender body and the small-amplitude assumptions, the elastic flagellum's radius κ and stiffness k κ are then optimized so as to maximize the motive force (for S p = 3.4). As the bead's radius r strongly affects both the motive and the total force, this optimization will result from the evolution of the ratio Fm Ftot as a function of the diameter ratioλ = 2r D depicted on Figures 7(a) to 7(l). Figures 7(a) to 7(l) depict the evolution of the ratio Fm Ftot , for BP, EF and HT designs (from left to right), as they navigate in a broad range of blood vessels: from artery (Figures 7(a The drag force F d scales in ∼ r 2 in large vessels and in ∼ r in smaller ones, where the drag [START_REF] Zhang | Artificial bacterial flagella for micromanipulation[END_REF] tends to the Stokes flow expression. Yet the motive force F m scales in r 3 for BP design (see [START_REF] Uecker | Real-time mri at resolution of 20 ms[END_REF]) and in r 2 for EF and HT designs (cf (4) and ( 5)). The smaller the robot's radius is, the more the drag decreases and it is even worst for the motive force. So a high radius, and in turn a high λ, improves the performances. However, drag is also affected by the partial occlusion of the vessel, through λ, which penalizes too high λ values. It results in a vessel-dependent optimal range for λ ranging from λ ∈ [0.15; 0.55] in an artery to λ ∈ [0.3; 0.6] in a small artery for BP design (see Figures 7(a) and7(d)). Previous studies [START_REF] Mathieu | Method of propulsion of a ferromagnetic core in the cardiovascular system through magnetic gradients generated by an mri system[END_REF], [START_REF] Vartholomeos | Simulation platform for selfassembly structures in mri-based nanorobotic drug delivery systems[END_REF] have predicted respective optimal ratios λ = 0.42 and λ = 0.25 in an artery: we show that taking advantage of the interaction forces results in a lower optimal value λ = 0.15, what reduces embolization hazards. In smaller vessels, we conversely predict almost the same profiles and optimal ratios than the ones evaluated in [START_REF] Vartholomeos | Simulation platform for selfassembly structures in mri-based nanorobotic drug delivery systems[END_REF].

EF design does not provide enough thrust in any vessel configuration ( Fm Ftot < 1), even if its efficiency increases as the robot enters smaller vessels (additional studies show that EF robot can resist only a maximum v f = 15µm.s -1 in a capillary). The same effect occurs for HT design, whose motive force also grows in r 2 , although Figures 7(i) and 7(l) illustrate that the thrust is high enough to drive such a robot in arterioles and capillaries, as emphasized in [START_REF] Abbott | Robotics in the small: Part 1 microrobotics[END_REF]. In macrovasculature, Figures 7(a) and 7(d) show that BP design is the best fitted to navigate in vessels where the drag is important, provided that the robot's trajectory is defined alongside the upper wall.

2) Reference trajectory: As the robot approaches the vessel's wall (as |∆| increases), the drag force decreases because of the parabolic profile [START_REF] Dreyfus | Microscopic artificial swimmers[END_REF] of the blood flow, while interaction forces F el , F vdw and F s increase and can either counterbalance (for the upper wall) or amplify (for the lower wall) the apparent weight W a contribution. Hence, the symmetry is broken in the ∆-axis, as shown in Figure 7. Next to the upper wall, attractive interaction forces equipoise the robot's apparent weight, which minimizes the control efforts. This effect is clearly visible in arteries (Figures 7(a) to 7(f)): since blood velocity is high, drag is dominant, and it is worthwhile navigating close to the wall. In arterioles and capillaries (Figures 7(g) to 7(l)), the problem is different: blood velocity is smaller and so the drag (though higher than available motive force for BP and ET designs). Interaction forces are consequently dominant, which prevents the robot from moving in the vicinity of the wall so to benefit from a reduced drag.

According to the forces balance presented in Section II-A and Figures 7(a) to 7(l), we can define freely the number of anchor points, e.g. four points in the present simulation. These points correspond to an optimal ratio Fm Ftot , as explained in the previous section, and are distributed along the path from the start point (injection seed) to the end point (e.g. tumor target). This path has to be close to the upper vessel wall to minimize the control efforts in arteries, whereas it should be centered in smaller vessels. To build an optimal reference trajectory the controller can track, these anchor points have to be joined in a smooth (C 2 ) curve. This interpolation is provided by B-splines for flexibility concerns: a change in one anchor point's position will result only in a local modification of the path, and the reference velocity along this path can be rescaled easily. This can be of particular interest if Fm Ftot < 1. For instance, figures 7(g) and 7(j) show that actual MRI devices cannot steer a BP robot in the smallest vessels because the drag force prevails. However, this limitation can be circumvented by computing a reference trajectory so as to synchronize the reference velocity with the blood flow.

B. Force balance and sensitivity study

The force balance plotted on Figure 8(a) for a small artery shows that the steric and van der Waals interactions can be omitted if considering the control of a microcarrier in large vessels, since they are respectively around 3 and 6 orders of magnitude under other forces. It is nonetheless necessary to take it into account in small aterioles and capillaries where they impact severely the force balance. Since the proposed model involves numerous uncertain parameters, a sensitivity study is also performed on Figure 8(b). The force balance is highly sensitive -linearly-to magnetization M , andnonlinearly-to blood permittivity ε, charge q, blood velocity V f . We have already showed in [START_REF] Arcese | Adaptive backstepping and mems force sensor for an mri-guided microrobot in the vasculature[END_REF] that such uncertainties on linear parameters like M or q 2 ε can be successfully addressed using adaptive backstepping.

C. Simulation results

The following simulation is performed by taking into account some limitations of a clinical MRI system. The magnetic field is strong enough (3T ) to assume that the particles reach their saturation magnetization. In order not to exceed the capacity of MRI systems, the applied control law is now corrected as

U a (t) = U(t) k(t) with k(t) = max 1, U(t)
Usat and a 100 mT m -1 s -1 slew rate.

The studied artery's radius is around 1.5mm and exhibits some minor and one major bifurcation. The most effective design is thus BP design (Figure 7(d)), and the optimal reference trajectory is set near the outer wall, as stated in Section III-A2. The radius of the robot is hence set to 250µm so as to respect the optimal ratio λ = 0.15, as stated in Section III-A1. Performances and stability of the controller with respect to noise measurement and parameters uncertainties are now illustrated by a closed-loop simulation, whose nominal parameters are given in Table III.

MRI and computed tomography devices currently have a spatial resolution around some hundreds of micrometers. Besides, since the microrobot is composed of a ferromagnetic nanoparticles aggregate, it induces artifacts on the imager of the MRI device. Hence we suppose that the position measurement is with an accuracy to within 100µm (i.e 10% of the measured signal), modeled by a white gaussian noise. Moreover, we assume that a 5% uncertainty affects the blood velocity, what models the modeling errors generated by i) the presence of minor bifurcations, ii) errors related to the approximation of CFD velocity profiles by a rational profile, and iii) the difficulty to estimate correctly the blood velocity.

Figure 9 (a) illustrates that the proposed controller-observer is stable. After the transient phase and despite the 100µm noise amplitude, the estimation error (Figure 9 (b)) is smaller than 10µm (i.e. 10% of nominal value). During the convergence of the observer, within t = 0s and t = 0.2s, the tracking error (Figure 9 (b)) increases. After this transient phase, one can notice that the observer reconstructs correctly the state. At t = 0.5s, important aberrant values are injected in the imaging measurements to test the robustness of the controller-observer

2 .N -1 .m -2 ] Blood velocity V f 0.025 [m.s -1 ] parameters a f , w f , φ f 1.15, 2π, π 2
Initial condition on X, Z X 0 , Z 0 (0, 0.001) T , (0, 0) T Initial condition on X, Ẑ X0 , Ẑ0 (0, 0) T , (0, 0) pair to bounded disturbances. It causes the microrobot to collide with the wall and the control inputs reach the actuator saturation (Figure 10 (a)), but the controller quickly rejects the disturbance to track the reference trajectory. After t = 0.6s, despite noise measurement and modeling error, the tracking is efficient and the error does not exceed 150µm. This is an excellent result since the model is very sensitive to errors on blood's velocity.

D. Discussion

The dynamic behaviour of the robot in blood vessels is highly affected by various physiological parameters (see Table II and Figure 8), e.g. the blood relative permittivity ǫ, viscosity η, density ρ f , etc. These parameters may vary with respect to many other factors such as temperature, patient age and health, etc... The selected control approach must be sufficiently robust to these parameters uncertainties to track correctly the pre-planned reference trajectory. Even though off-line identification of some biological parameters is a solution, invivo estimation remains a difficult issue. That is why we have developed an adaptive control [START_REF] Arcese | Adaptive backstepping and mems force sensor for an mri-guided microrobot in the vasculature[END_REF] to on-line estimate some physiological parameters, which greatly improves the robustness to modeling errors on linear parameters. The estimation of non-linearly varying parameters requires further study and remains an outstanding problem.

Poor performances of EF design should be compared to the maximum velocity that other similar designs can reach within the same constraints: 4µm.s -1 for a paramagnetic beating flagellated robot [START_REF] Dreyfus | Microscopic artificial swimmers[END_REF], 30µm.s -1 for E. Coli, and over 200µm.s -1 for magnetotactic bacteria [START_REF] Martel | Flagellated magnetotactic bacteria as controlled mri-trackable propulsion and steering systems for medical nanorobots operating in the human microvasculature[END_REF]. In the last two cases, the improved efficiency is a result of the distributed actuation in natural flagella. The shape anisotropy and flagellum design of both EF and HT designs make it difficult to i) manoeuver in vessels and anticipate wall interactions with the tail, and to ii) consider a progressive downscaling so as to benefit from an optimal ratio while navigating from macro to microvasculature. Since such designs are well suited for explorating small vessels, it is tempting to embed such 5µm devices inside a bigger agglomerate. The latter could be bead pulled in arteries, then break down to release the flagellated microrobots in the arterioles. However such a design requires both strong magnetic gradients and oscillating magnetic field: in any case, a clinical MRI device cannot fulfill these requirements.

Spatial resolution of current medical imagers is an important limitation, especially within real-time constraints [START_REF] Achenbach | Contrast-enhanced coronary artery visualization by dual-source computed tomographyinitial experience[END_REF], [START_REF] Zhang | Real-time magnetic resonance imaging of temporomandibular joint dynamics[END_REF]. However, future imagers like 7T MRI or 256 slice computed tomography devices are promising for high spatial and temporal resolution imaging. This will be a key point for implementing targeted minimally invasive therapies using microrobots that require to navigate into smaller vessels. Using an MRI device to provide both imaging and thrust also requires multiplexing approach, what induces both a lower imaging resolution and a degraded tracking during imaging period [START_REF] Tamaz | Real-time mri-based control of a ferromagnetic core for endovascular navigation[END_REF]. Imaging post-processing -notably time required to perform accurate localization of the microrobot using its magnetic signature despite a size below the MRI resolution [START_REF] Shapiro | Mri detection of single particles for cellular imaging[END_REF]-and delays can also severely degrade the proposed approach. A way to overcome these downsides is either to upgrade the MRI devices with additionnal gradient coils, either to use an additionnal device (stereotactic device from Stereotaxis company www.stereotaxis.com, or Octomag device from Aeon Scientific www.aeon-scientific.com) inside the imager for propulsion purpose. Slowing down the reference trajectory may also reduce discrepancies induced by imaging delays. However, addressing properly stability concerns induced by delays requires a Lyapunov-Razumikhin/Krasovskii based backstepping.

Controlling such microrobots in the vasculature raises safety concerns. The two main limitations are inadequate embolization and magnetic hazards. Designing the microrobot as a degradable polymer aggregate is an answer to the first issue. Using time-varying magnetic fields can induce neurostimulation: at the present time, there are no worldwide regulations for time-varying magnetic gradients. Our works consequently deal with standard slew rates available on clinical MRI devices.

We approximate the blood velocity given by CFD results by analytical functions for quite simple vascular geometry. A more generic and robust approach relies on the design of an observer of the blood velocity, i.e. on an on-line estimation of the blood flow from the microrobot's position and heartbeat measurements, which is still an open problem.

IV. CONCLUSION

This paper proposes a physics-based modeling framework to optimize the navigation of magnetically guided therapeutic microrobots in blood vessels. The robot consists of a polymer binded aggregate of nanosized ferromagnetic particles grafted with micelles to drugs. The model of the robot takes into account wall effects, wall interactions and non-Newtonian behaviour of blood. The benefits of using a model based approach are illustrated through the definition of both an optimal size of the microrobot and an optimal trajectory so as to minimize the control efforts. The proposed optimization framework makes it possible to generalize the navigation process to complex vasculature networks. Taking into account the geometry of the vessels, the blood flow velocity at bifurcations and blood properties, we can optimize targeting with an appropriate backstepping control law coupled with a high gain observer. The performances and robustness of the approach are validated by simulation results including noise, disturbances, and model parameters uncertainties.

Fig. 1 .

 1 Fig. 1. Complexity of endovascular navigation for a magnetic microrobot from macro to micro-vasculature.

  where L, ζ ⊥ = 4πη ln L 2κ +0.193 and κ ≪ r are respectively the length, the normal viscous drag coefficient, and the radius of the flagellum. k κ = πEκκ 4 4 denotes the bending stiffness of the tail, < • > denotes the time average operator over one period, and Z is a scaling function dependent on the Sperm number S p = L 4 ωζ ⊥ kκ 1/4

Fig. 2 .

 2 Fig. 2. Different microrobots designs. (a) Bead pulling: control input is the magnetic gradient ∇ B. (b) Elastic flagella: control inputs are the time-average magnetic field direction and the frequency. (c) Helical tail: control inputs are the magnetic torque direction and the frequency.

11

 11

Fig. 4 .

 4 Fig.4. Scheme of a blood vessel with minor bifurcations.

  (a) Blood velocity profiles. (b) Velocity fields.

Fig. 5 .

 5 Fig. 5. CFD results of a parent vessel with minor bifurcations. (a) Blood velocity profiles: CFD results (dotted lines) and parabolic approximation (solid line). (b) Blood velocity fields: CFD results using COMSOL software.

  Blood's velocity (m/s) (Steady state profile in parent branch) Distance to the upper wall (m) 10 -3 θ = 5 • θ = 10 • θ = 15 • θ = 20 • θ = 25 • θ = 30 • θ = 0 • (a) Blood velocity profiles before the apex. (b) Velocity profiles in a major bifurcation. Distance to the upper wall (m) Blood's velocity (m/s) Steady state profile in daughter branch 10 -3 θ = 35 • θ = 40 • θ = 45 • (c) Blood velocity profiles after the apex.

Fig. 6 .

 6 Fig. 6. CFD results of a parent vessel with a major bifurcation. (a) Blood velocity profiles before the apex: CFD results (dotted lines) and rational approximation (solid lines). (b) Blood velocity profiles: CFD results using FLUENT software. (c) Blood velocity profiles after the apex: CFD results (dotted lines) and rational approximation (solid lines).

) to 7 (

 7 c)), small artery (Figures 7(d) to 7(f)) and arteriole (Figures 7(g) to 7(i)), to capillary (Figures 7(j) to 7(l)).

Fig. 7 .

 7 Fig. 7. Evolution of the ratio Fm/Ftot as a function of λ and ∆. From left to right: bead pulled, elastic flagellated and helical tailed robots. From top to bottom, the robots navigate from artery through small artery and arterioles to capillary.

  Force balance as a function to the distance to the vessel centerline.

  Sensitivity to some uncertain parameters.

Fig. 8 .

 8 Fig. 8. (a) Force balance in a small artery. (b) Sensitivity study in a small artery.

Fig. 9 .Fig. 10 .

 910 Fig. 9. (a) XZ trajectory: reference trajectory (black dotted line), robot's trajectory (blue solid line) and estimated trajectory (red dashed line). (b) Error: tracking error (blue solid line) estimation error (red dashed line)

This work was supported by European Union's 7th Framework Program and its research area ICT-2007.3.6 Micro/nanosystems under the project NANOMA (Nano-Actuactors and Nano-Sensors for Medical Applications). L. Arcese and M.