Fabio Bagagiolo

Rosario Maggistro

Marta Zoppello

Keywords: Scallop theorem, Switching, Thermostat, Controllability

come

Introduction

The study of locomotion strategies in fluids is attracting increasing interest in recent literature, especially for its connection with the realization of artificial devices that can self-propel in fluids. Theories of swimming generally utilize either low Reynolds number approximation, or the assumption of inviscid ideal fluid dynamics (high Reynolds number). These two different regimes are also distinct in terms of the mechanism of locomotion [START_REF] Childress | Mechanics of swimming and Flying Cambridge[END_REF][START_REF] Lighthill | Mathematical Biofluiddynamics[END_REF].

In this paper we focus on swimmers immersed in these two kind of fluids which produce a linear dynamics. In particular we study the system describing the motion of a scallop for which it is well known [START_REF] Alouges | Optimal strokes for low Reynolds number swimmers: an example[END_REF][START_REF] Munnier | Generalized scallop theorem for linear swimmers[END_REF][START_REF] Purcell | Life at low Reynolds number[END_REF]] that the scallop theorem/paradox holds. This means that it is not capable to achieve any net motion performing cyclical shape changes, either in a viscous or in an inviscid fluid. Some authors tried to overcome this paradox changing the geometry of the swimmer, for example adding a degree of freedom, introducing the Purcell swimmer [START_REF] Purcell | Life at low Reynolds number[END_REF], or the three sphere swimmer [START_REF] Golestanian | Analytic results for the threesphere swimmer at low Reynolds number[END_REF]. Others, instead, supposed the scallop immersed in a non Newtonian fluid, in which the viscosity is not constant, ending up with a non reversible dynamics [START_REF] Cheng | Hydrodynamics of scallop locomotion: unsteady fluid forces on clapping shells[END_REF][START_REF] Qiu | Swimming by reciprocal motion at low Reynolds number[END_REF]. Inspired by this last approach, our aim is to propose some strategies which maintain the swimmer geometry and exploit instead a change in the dynamics. The idea is based on switching dynamics depending on the angular velocity of opening and closing of the scallop's valves. More precisely we analyze two cases: in the first one we suppose that if the modulus of the angular velocity is high, the fluid regime can be approximated by the ideal one, instead if this modulus is low the fluid can be considered as completely viscous. These assumptions are realistic since the Reynolds number changes depending on the characteristic velocity of the swimmer. In the second case we assume that the fluid reacts in a different way between the opening and closing of the valves: it facilitates the opening, so that it can be considered an ideal fluid, and resists the closing, like a viscous fluid. These last approximations model a fluid in which the viscosity changes with the sign of the angular velocity. More precisely we use two constant viscosities: one high (resp. one very small) if the angular velocity is negative (resp. positive). Moreover inspired by [START_REF] Qiu | Swimming by reciprocal motion at low Reynolds number[END_REF], where the scallop's opening and closing is actuated by an external magnetic field, in this last case we also introduce an hysteresis mechanism through a thermostat, see Fig 5 (see [START_REF] Visintin | Differential Models of Hysteresis Springer Hiderberg[END_REF] for mathematical models for hysteresis), to model a delay in the change of fluid's regime. In both cases we assume to be able to prescribe the angular velocity, using it as a control parameter and we prove that the system is controllable, i.e. the scallop is able to move both forward and backward using cyclical deformations. Furthermore we prove also that it is always possible to move between two fixed points, starting and ending with two prescribed angles. In the last part of the paper we show also some numerical examples to support our theoretical predictions.

The plan of the paper is the following. In Section 2 we present the swimmer model and derive its equation of motion both in the viscous and in the ideal approximation, proving the scallop theorem. Section 3 is devoted to the introduction of the switching strategies which lead to the controllability of the scallop system. Finally in Section 4 we present some numerical simulations showing different kind of controls that can be used.

The Scallop swimmer

In this section we are interested in analyzing the motion of an articulated rigid body immersed in a fluid that changes its configuration. In order to determine completely its state we need the position of its center of mass and its orientation. Their temporal evolution is obtained solving the Newton's equations coupled with the Navier-Stokes equations relative to the surrounding fluid. We will face this problem considering the body as immersed in two kinds of different fluids: one viscous at low Reynolds number in which we neglect the effects of inertia, and another one ideal inviscid and irrotational, in which we neglect the viscous forces in the Navier-Stokes equations. First of all we recall that in both cases a swimmer that tries to moves like a scallop, opening and closing periodically its valves, does not move at the end of a cycle. This situation is well known as scallop theorem (or paradox) [START_REF] Alouges | Optimal strokes for low Reynolds number swimmers: an example[END_REF][START_REF] Purcell | Life at low Reynolds number[END_REF]. In what follows we will consider a planar body composed by two rigid valves of elliptical shape, joined in order that they can be opened and closed. Moreover this body is constrained to move along one of the cartesian axes (the e x -axis) and is symmetric with respect to it. Finally we will neglect the interaction between the two valves. The configuration of the system is easily described by the position x of the juncture point along the e x -axis and by the angle θ that each valve forms with the axis ex x θ Fig. 1 The scallop configuration

The possible translation of the system is determined by the consecutive opening and closing of the valves. Our aim is to determine the net translation of the body, given the function of time describing the angular velocity θ.

Viscous fluid

Here we focus on the case in which the scallop is immersed in a viscous fluid. In this regime the viscous forces dominates the inertial ones that can be neglected, so the equations governing the dynamics of the fluid are the Stokes ones:

∆v -∇p = 0 together with the incompressibility condition div v = 0. Let us consider that the ellipses have major axis 2a and minor axis 2b with b << a, moreover let us suppose that θ ∈ (0, π 2 ) so that it remains acute. One of the main difficulties in computing explicitly the equation of motion is the complexity of the hydrodynamic forces exerted by the fluid on the swimmer as a reaction to its shape changes. Since in our assumptions the minor axis of the ellipse is very small with respect to the major one, i.e. b << a, we can consider the swimmer as one-dimensional, composed essentially by two links of length 2a (see Fig 1). In the case of slender swimmers, Resistive Force Theory (RFT) [START_REF] Gray | The propulsion of sea-urchin spermatozoa[END_REF] provides a simple and concise way to compute a local approximation of such forces, and it has been successfully used in several recent studies, see for example [START_REF] Becker | On selfpropulsion of micro-machines at low Reynolds number: Purcell's three-link swimmer[END_REF][START_REF] Friedrich | High-precision tracking of sperm swimming fine structure provides strong test of resistive force theory[END_REF]. From now on we use this approach as well, in order to obtain the forces acting on the swimmer, neglecting the interaction between the valves. Since the scallop is immersed in a viscous fluid the inertial forces are negligible with respect to the viscous ones, therefore the dynamics of the swimmer follows from Newton laws in which inertia vanishes:

F = 0 ( 1
)
where F is the total force exerted on the swimmer by the fluid. As already said we want to couple the fluid and the swimmer, using the local drag approximation of Resistive Force Theory. We denote by s the arc length coordinate on the i-th link (0 ≤ s ≤ 2a) measured from the juncture point and by v i (s) the velocity of the corresponding point. We also introduce the unit vectors e 1 = cos(θ) sin(θ) , e ⊥ 1 =

-sin(θ) cos(θ) , and e 2 = cos(θ) -sin(θ) ,

e ⊥ 2 =
-sin(θ) -cos(θ) in the directions parallel and perpendicular to each link and write the position of the point at arc length s as x i (s) = x 0 + se i where x is the coordinate of the joint between the two valves. By differentiation, we obtain,

v i (s) = ẋ 0 + s θi e ⊥ i . (2)
The density of the force f i acting on the i-th segment is assumed to depend linearly on the velocity. It is defined by

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

F = 2a 0 f 1 (s) ds + 2a 0 f 2 (s) ds = 0 (4)
Using ( 2) and ( 3) and since we are neglecting inertia we have

F x = -4aξ ẋ cos 2 (θ) -4aη ẋ sin 2 (θ) + 4a 2 η θ sin(θ) = 0 F y = 0 (5)
Observe that F y vanishes since the scallop is symmetric with respect to the e x axis. From ( 5) is now easy to determine the evolution of

x ẋ = V 1 (θ) θ = aη sin(θ) ξ cos 2 (θ) + η sin 2 (θ) θ (6)

Ideal Fluid

While in the previous subsection we faced the problem of the self-propulsion of the scallop immersed in a viscous fluid, here we focus on the case in which it is immersed in an ideal inviscid and irrotational fluid. Let us make the same assumptions on the parameters a and b that have been done in the previous section, moreover let us denote by Ω the region of the plane occupied by the swimmer in a reference configuration.

Assigning (x, θ) as functions of time let us call

f (x,θ) : Ω → R 2 ζ → f (x,θ) (ζ)
the function which maps each point of the swimmer ζ ∈ Ω in f (x,θ) (ζ) that is its position in the plane at time t. Supposing that θ can be assigned and that there are not other external forces, our aim is to find equations that describe the motion of x. To this end we call v the velocity of the fluid, its motion is given by the Euler equations for ideal fluids

v t + v • ∇v = -∇p (7)
with the incompressibility condition div v = 0. Moreover we impose a Neumann boundary condition, that is that the normal component of the velocity of the fluid has to be equal to the normal component of the velocity of the body.

v(f (x,θ) ) - ∂f (x,θ) ∂x ẋ + ∂f (x,θ) ∂θ θ , n (x,θ) = 0
where • denotes the scalar product, n (x,θ) is the external normal to the set f (x,θ) (Ω). To find the evolution of x we should solve the Lagrange equation

d dt ∂T b ∂ ẋ = ∂T b ∂x + F (8)
where T b is the kinetic energy of the body and F the external pressure force acting on the boundary of the swimmer. As already done in [START_REF] Bressan | Alberto Impulsive Control of Lagrangian Systems and locomotion in fluids[END_REF][START_REF] Mason | Propulsion and control of deformable bodies in a Ideal fluid[END_REF][START_REF] Munnier | Locomotion and control of a self-propelled shape-changing body in a fluid[END_REF] this force F can be reinterpreted as a kinetic term, precisely thanks to the fact that we are in an ideal fluid. Therefore the system body + fluid is geodetic with Lagrangian given by the sum of the kinetic energy of the body (T b ) and the one of the fluid (T f ):

T tot = T b + T f
The kinetic energy of the body is the sum of the kinetic energy of the two ellipses, that reads

T b = m ẋ2 + a 2 θ2 -2a ẋ θ sin θ + I θ2 (9)
Since we are dealing with an ideal fluid and thus inertial forces dominates over the viscous ones, in order to derive the kinetic energy of the fluid we will make use of the concept of added mass. In fluid mechanics, added mass or virtual mass is the inertia added to a system because an accelerating or decelerating body must move (or deflect) some volume of surrounding fluid as it moves through it. Added mass is a common issue because the object and surrounding fluid cannot occupy the same physical space simultaneously [START_REF] Bessel | Investigations on the length of the seconds pendulum Königlichen Akademie der Wissenschaften[END_REF]. For simplicity this can be modeled as some volume of fluid moving with the object, though in reality "all" the fluid will be accelerated, to various degrees. Therefore the kinetic energy of the fluid will be given by the sum of the kinetic energy of the added masses of the two ellipses.

T f = 1 2 v T 1 M 1 add v 1 + 1 2 v T 2 M 2 add v 2 (10)
where M i add are the added mass matrices relative to each ellipse which are diagonal, and v i the velocities of their centre of mass, expressed in the frame solidal to each ellipse with axes parallel and perpendicular to the major axis. Finally we can compute the total kinetic energy of the coupled system body+ fluid that is

T tot = m ẋ2 + a 2 θ2 -2a ẋ θ sin θ + I θ2 + + m 11 ẋ2 cos 2 θ + m 22 ẋ2 sin 2 θ + a 2 θ2 -2a ẋ θ sin θ + + m 3 3 θ2 (11)
Following a procedure introduced by Alberto Bressan in [START_REF] Bressan | Alberto Impulsive Control of Lagrangian Systems and locomotion in fluids[END_REF], in order to end up with a control system we perform a partial legendre transformation on the kinetic energy defining

p = ∂T tot ∂ ẋ = 2 ẋ m + m 11 cos 2 θ + m 22 sin 2 θ -2a θ sin θ(m + m 22 ) from which we derive ẋ = p + 2a θ sin θ(m + m 22 ) 2(m + m 11 cos 2 θ + m 22 sin 2 θ) (12)
There is a wide spread literature regarding the computation of added masses of planar contours moving in an ideal unlimited fluid. We will use in the rest of the paper the added mass coefficients for the ellipse computed in [START_REF] Newman | Marine hydrodynamics[END_REF]: the added mass in the direction of the major axis is m 11 = ρπb 2 , the one along the minor axis is m 22 = ρπa 2 . Notice now that writing the Hamilton equation relative to p, and recalling (11

) ṗ = ∂T tot ∂x = 0
thus, if we start with p(0) = 0, p remains null for all times and the evolution of

x becomes ẋ = V 2 (θ) θ = a sin θ(m + ρπa 2 ) m + ρπb 2 cos 2 θ + ρπa 2 sin 2 θ θ (13)
Theorem 1 (Scallop Theorem) Consider a swimmer dynamics of the type

ẋ = V (θ) θ ( 14
)
Then for every T -periodic deformation (i.e. stroke) one has

∆x = T 0 ẋ(t) dt = 0 ( 15
)
that is, the final total translation is null Proof Define the primitive of V by

F (θ) = θ 0 V (σ) dσ (16)
Then using ( 14)

∆x = T 0 V (θ(t)) θ(t) dt = T 0 d dt F (θ(t)) dt = F (θ(T )) -F (θ(0)) = 0
by the periodicity of t → θ(t).

Note that the dynamics (6) and (13) are of the type (14), therefore the scallop theorem is valid either in the viscous and in the ideal case.

Controllability

In this section we will give two different strategies to overcome the scallop theorem, both based on a switching mechanism. In particular we produce some partial and global controllability results for this switching systems.

Partial controllability in x

We have previously seen that if our scallop is immersed either in an ideal fluid or in a viscous one, if it experiences periodical shape changes it is not able to move after one cycle. Here we would like to find a way to overcome this problem. The main idea is to be able to change the dynamics during one periodical stroke and see if in this way we obtain a net motion and in particular some controllability. In order to do this we have to introduce the Reynolds number, a number which characterizes the fluid regime. It arises from the adimesionalization of the Navier-Stokes equations and it is defined by

Re = V Lρ η = V L ν ( 17
)
where V is the characteristic velocity of the body immersed in the fluid, L its characteristic length, ρ the density of the fluid, η its viscosity and ν = η ρ is the kinematic viscosity. The Reynolds number quantifies the relative importance of inertial versus viscous effects.

η = η(| θ|)

Let us recall that if v(t, x) is a solution of the Navier Stokes equations, the function u(t, x) = v(ct, x), c > 0 is still a solution of the Navier Stokes equations but with a different viscosity. Now assume that the absolute value of the speed θ is very high, this means that rescaling the time of the solution of the Navier Stokes equations, we end up with a viscosity η that is very small and therefore the Reynolds number is large. In this case the inertial forces dominates over the viscous ones, so we can consider the scallop immersed in an ideal fluid and thus use the dynamics [START_REF] Munnier | Generalized scallop theorem for linear swimmers[END_REF]. Then we suppose that at a certain point of the cycle the absolute value of the angular velocity is very small. In this case we have a solution of the Navier Stokes equations with a very high viscosity η. Thus we can suppose that the scallop is immersed in a Stokes fluid, since the viscous effects dominates the inertial ones and use the dynamics [START_REF] Cheng | Hydrodynamics of scallop locomotion: unsteady fluid forces on clapping shells[END_REF]. This situation is well represented by a switching system in which the change of the dynamics is determined by the modulus of the angular velocity θ: if it is big (i.e | θ| > M with M > 0) we use the ideal approximation and the corresponding dynamics; if it is small (i. Our aim is to prove that using this kind of switching we are able to have a net displacement, both forward or backward, using periodic continuous functions θ

According to what said before we can prescribe the angular velocity θ and thus use it as a control function u. Therefore we write the system as a control system that is

           ẋ(t) = V w(t) (θ(t))u(t), θ(t) = u(t) w(t) = h[u](t) x(0) = x 0 , θ(0) = θ 0 w(0) = w 0
where u is continuous and

h[u] = 2 if |u| > M 1 if |u| < M
Fig. 2 The rule of the classical switching Moreover let us call F i the primitives of the functions V i , for i = 1, 2. They are :

F 1 = aη arctan h( η-ξ η cos θ) η(η -ξ) , F 2 = -a m + a 2 ρπ arctan h( √ (a 2 -b 2 )ρπ cos θ √ m+a 2 ρπ ) ρπ(a 2 -b 2 )
Theorem 2 With the previous switching scheme we are able to overcome the Scallop paradox, thus to move both forward and backward. More precisely there are r > 0 small enough (see remark 1), a final time T > 0 and a continuous T -periodic control function u(t), which make the system move between two fixed points along the x axis, x 0 and x f ∈]x 0 -r, x 0 + r[, in the time T .

Proof First case: u(0) > M In this case we start with the ideal approximation (i.e w 0 = 2)

V w(t) (θ(t)) =                V 2 (θ(t)) 0 < t < t 1 , V 1 (θ(t)) t 1 < t < t 2 V 2 (θ(t)) t 2 < t < t 3 V 1 (θ(t)) t 3 < t < t 4 V 2 (θ(t)) t 4 < t < T (18)
with

t 1 := inf {T > t > 0 | u(t) = M } and t 2 := inf {T > t > t 1 | u(t) = -M } and t 3 := inf {T > t > t 2 | u(t) = -M } and t 4 := inf {T > t > t 3 | u(t) = M }
assuming that inf (∅) = +∞. The net motion is then calculated as

∆x =(F 2 -F 1 )(θ(t 1 )) + (F 2 -F 1 )(θ(t 3 )) -(F 2 -F 1 )(θ(t 2 )) -(F 2 -F 1 )(θ(t 4 )). ( 19
)
taking into account that θ(0) = θ(T ) and that (F 2 -F 1 )(θ(t i )) does not appear in the equation if t i = +∞. We want to prove that we are able to move choosing a suitable periodic evolution for our control function θ = u. Let us call the unknowns θ i := θ(t i ), for i = 1 . . . 4.

First of all we show that ∆x as function of (θ

1 , θ 2 , θ 3 , θ 4 ) is surjective in ]0, π 2 [×]0, π 2 [×]0, π 2 [×]0, π 2
[. We are going to prove that

∇(∆x) =     -(V 2 -V 1 )(θ 1 ) (V 2 -V 1 )(θ 2 ) -(V 2 -V 1 )(θ 3 ) (V 2 -V 1 )(θ 4 )     = 0 in (θ 1 , θ 2 , θ 3 , θ 4 ) ∈]0, π 2 [×]0, π 2 [×]0, π 2 [×]0, π 2 [
so that ( 19) is a submersion and surjective as required.

Recall that the function (F 2 -F 1 )(•) is always increasing indeed

(F 2 -F 1 )(θ) ∂θ = - aη ξ cos 2 θ + η sin 2 θ + ma + ρπa 2 m + ρπb 2 cos 2 θ + ρπa 2 sin 2 θ sin θ = sin θ cos 2 θ ma(η -ξ) + ρπ(ξa 2 -ηb 2 ) (m + ρπb 2 cos 2 θ + ρπa 2 sin 2 θ)(ξ cos 2 θ + η sin 2 θ) > 0 for θ ∈]0, π 2 [ and b << a (20)
From this immediatly follows that ∇(∆x) = 0. The surjectivity ensures us that for any fixed ∆x in a neighborhood of zero we are always able to find a (θ 1 , θ 2 , θ 3 , θ 4 ) which realize the desired displacement. Moreover, thanks to the symmetry properties of the function defining the displacement, also each of the 4uplets (θ 1 , θ 4 , θ 3 , θ 2 ), (θ 3 , θ 4 , θ 1 , θ 2 ) and (θ 3 , θ 2 , θ 1 , θ 4 ) realizes the same displacement. Supposing ∆x > 0 and recalling that the function (F 2 -F 1 )(•) is increasing, then the angles (θ 1 , θ 2 , θ 3 , θ 4 ) will have a suitable order that can or not be coherent with the switching rule and the periodicity of θ. If their sorting is appropriate we will choose a control θ = u such that θ(t i ) = θ i . Otherwise at least one of the 4-uplets above will be right. Thus defining (θ 1 , θ 2 , θ 3 , θ 4 ) this latter uple, we take a control u such that θ(t i ) = θ i . This choice of the control will lead us to obtain the desired positive displacement.

For example suppose that the uplet (θ 1 , θ 2 , θ 3 , θ 4 ) which realizes the desired positive displacement, satisfy

θ 3 > θ 4 > θ 1 > θ 2 . Indeed (F 2 -F 1 )(θ 1 ) -(F 2 -F 1 )(θ 2 ) > 0 =⇒ ∆x > 0 (F 2 -F 1 )(θ 4 ) -(F 2 -F 1 )(θ 3 ) < 0
To respect the switching scheme in the time interval (t 2 , t 3 ) the function θ should decrease and thus θ 2 > θ 3 . The latter is not satisfied by (θ 1 , θ 2 , θ 3 , θ 4 ), but taking (θ 1 , θ 2 , θ 3 , θ 4 ) = (θ 3 , θ 4 , θ 1 , θ 2 ), we have the same ∆x and the switching scheme is now respected. Analogous arguments can be used if ∆x < 0. Second case: -M < u(0) < M In this case we start form the viscous approximation (i.e w 0 = 1). Using arguments similar to the ones used before to compute ∆x and to prove its surjectivity, and redefining accordingly the times t i for i = 1 . . . 4. we have that

∆x =(F 2 -F 1 )(θ(t 2 )) + (F 2 -F 1 )(θ(t 4 )) -(F 2 -F 1 )(θ(t 1 )) -(F 2 -F 1 )(θ(t 3 )). ( 21
)
going on as before, exploiting the surjetivity and the symmetry of the last function, we are able to find a control u that realizes the desired displacement.

Third case: u(0) < -M This case is analogous to the first one. In conclusion we have proved that wherever we start on the switching diagram we are able to achieve a net displacement either positive or negative and then we have the controllability.

Remark 1 Note that the value of r in the the last theorem is the maximal value that the function

|∆x(θ 1 , θ 2 , θ 3 , θ 4 )| can assume in ]0, π 2 [×]0, π 2 [×]0, π 2 [×]0, π 2
[. Thus the constant r is independent from x and θ. To cover distances |∆x| ≥ r we should divide the spatial interval in N subintervals of length less than r, each one realized by a u of period T N . Repeating N times this control u we are able to reach the desired displacement.

η = η(sign( θ))

While in the previous subsection we supposed that the change in the fluid regime was linked to the magnitude of the modulus of the angular velocity, here we would like to link the two fluids approximations to the sign of θ. This model can be used to describe a different response of the fluid to the opening and closure of the scallop' valves. For instance the situation in which the fluid has a pseudoelastic nature that assists the valve opening but resist the valve closing [START_REF] Cheng | Hydrodynamics of scallop locomotion: unsteady fluid forces on clapping shells[END_REF]. Thus, according to our assumption the viscosity of the fluid changes between the opening and the closing of the valves, switching from one constant value to another one. This can be represented by a switching scheme as in Fig. 4. If the valves are opening ( θ > 0) we suppose that the fluid is not opposing resistance, as if the scallop is immersed in an ideal fluid; instead when the valves are closing ( θ < 0) the fluid is opposing a big resistance, and we can consider the scallop immersed in a viscous fluid.

The system can be written as a control system, in which the control function u(t) is the angular velocity θ:

           ẋ(t) = V w(t) (θ(t))u(t), θ(t) = u(t) w(t) = h[u](t) x(0) = x 0 , θ(0) = θ 0 w(0) = w 0
where the control u is continuous and now Proof Let us suppose to start with the ideal approximation, so that we are opening the valves u(0) > 0 and w 0 = 2,

h[u] = 2 if u > 0 1 if u < 0
V w(t) (θ(t)) =      V 2 (θ(t)) 0 < t < t 1 , V 1 (θ(t)) t 1 < t < t 2 , V 2 (θ(t)) t 2 < t < T. (22)
with

t 1 := inf {T > t > 0 | u(t) = 0} and
t 2 := inf {T > t > t 1 | u(t) = 0}
with inf (∅) = +∞. The net motion can be computed as

∆x = F 2 (θ(t 1 ))+F 1 (θ(t 2 ))-F 1 (θ(t 1 ))-F 2 (θ(t 2 )). ( 23
)
recalling as before that θ(0) = θ(T ). We want to prove that we are able to move choosing a suitable periodic evolution for our control function θ = u. Let us call θ 1 := θ(t 1 ) and θ 2 := θ(t 2 ), first of all we show that ∆x as function of ( 23) is a submersion and surjective as required.

θ 1 , θ 2 ) is surjective in ]0, π 2 [×]0, π 2 [. Like before we prove that ∇(∆x) = (V 2 -V 1 )(θ 1 ) (V 1 -V 2 )(θ 2 ) = 0 in (θ 1 , θ 2 ) ∈]0, π 2 [×]0, π 2 [ hence (
Notice that

∆x = (F 2 -F 1 )(θ 1 ) -(F 2 -F 1 )(θ 2 )
If we chose a control such that θ 1 > θ 2 then ∆x will be positive, while if θ 1 < θ 2 then ∆x will be negative. But since we need to respect the switching rule the last case could not be achieved because after t 1 θ = u < 0 and thus we are closing the valves therefore θ(t 2 ) = θ 2 will be necessarily less than θ(t 1 ) = θ 1 .

The case where u(0) < 0 is analogous to the previous one.

In conclusion we have proved that for every choice of w 0 we are able to achieve a net displacement but only forward.

Thermostatic case

In this section we introduce a variant of the previous switching in order to be able to move both forward and backward and therefore have a result of partial controllability in x. Our approach is to link the variation of u = θ ∈ R by a delayed thermostat, an operator with memory, introduced rigorously in [START_REF] Visintin | Differential Models of Hysteresis Springer Hiderberg[END_REF], consisting of two different thresholds for passing separately from one edge to the another one and vice-versa. This idea was inspired by [START_REF] Qiu | Swimming by reciprocal motion at low Reynolds number[END_REF] in which the Scallop opening and closing is actuated by an external magnetic field, and thus a delay mechanism is reasonable. We suppose that the dynamics V depends on the angle θ ∈]0, π 2 [, and also depends on a discrete variable w ∈ {1, 2}, whose evolution is governed by a delayed thermostatic rule, subject to the evolution of the control u. In Fig. 5 the behavior of such a rule is explained, correspondingly to the choice of a fixed threshold parameter ε > 0. The output w ∈ {1, 2} may jump from 2 to 1 only when the input u is equal to -ε, and must jump when U coming from the right (i.e. from values larger than or equal to -ε), possibly goes below the threshold -ε; it may jump from 1 to 2 only when u is equal to ε, and must jump when it comes from the left (i.e. from values smaller than or equal to ε) possibly goes above the threshold ε. In all other situations it remains locally constant in time. In particular, when u > ε then w is equal to 2, and when u < -ε then w is equal to 1.

Fig. 5 The thermostatic approximation

The controlled evolution is then given by

           ẋ(t) = V w(t) (θ(t))u(t), θ(t) = u(t) w(t) = h ε [u](t) x(0) = x 0 , θ(0) = θ 0 w(0) = w 0 (24)
where h ε [•] represents the thermostatic delayed relationship between the input u and the output w. Note that the initial value w 0 ∈ {1, 2} must be coherent with the thermostatic relation: w 0 = 2 (resp. w 0 = 1) whenever θ0 > ε (resp. θ0 < -ε).

We start now to analyse the value of the displacement ∆x depending of the value of u proving the following result:

Theorem 4 Let x f ∈]x 0 -r, x 0 + r[ with r > 0 small enough. Then, there always exits a time T > 0 and a continuous T -periodic control function θ = u (hence a periodic θ) such that one can move from x 0 to x f in time T when the delayed thermostat is taken into account. In other words the system (24) is partially controllable in x.

Proof First case -ε < u(0) < ε and w 0 = 1 then we have

V w(t) (θ(t)) = V 1 (θ(t)) 0 < t < t 1 V 2 (θ(t)) t 1 < t < T. ( 25
)
where t 1 is the first time for which u goes through ε, i.e.

t 1 := inf {T > t > 0 | u(t) = ε}
and T is the final time. The displacement is then

∆x = F 1 (θ(t 1 )) -F 1 (θ(0)) + F 2 (θ(0)) -F 2 (θ(t 1 )). ( 26
)
recallin as before that θ(0) = θ(T ). We call θ(t 1 ) = θ 1 and we want to prove that we able to obtain ∆x = c, ∀ |c| < r using a suitable periodic control function. In order to do this we show that ∆x(θ 1 ) is surjective in a neighborhood of zero. First of all we compute the derivative and show that it is different from 0 and negative.

∂∆x ∂θ 1 = V 1 (θ 1 ) -V 2 (θ 1 ) = aη sin θ 1 ξ cos 2 θ 1 + η sin 2 θ 1 - (ma + ρπa 2 ) sin θ 1 m + ρπb 2 cos 2 θ 1 + ρπa 2 sin 2 θ 1 = sin θ 1 cos 2 θ 1 -ma(η -ξ) -ρπ(ξa 2 -ηb 2 ) (m + ρπb 2 cos 2 θ 1 + ρπa 2 sin 2 θ 1 )(ξ cos 2 θ 1 + η sin 2 θ 1 ) = 0 for θ 1 ∈]0, π 2 [
Notice also that since in our assumptions b is negligible with respect to a, i.e b << a, we have that ηb 2 << ξa 2 and thus the derivative is always negative and consequently the ∆x is decreesing. We are interested in θ 1 ∈]0, π 2 [. Since the derivative of the function defining the displacement is different from 0 in ]0, π 2 [, (26) is locally invertible. Thus, since the inverse image of 0 is θ 0 then the inverse image of a neighborhood of 0 is a neighborhood of θ 0 . Finally, recalling that ∆x is decreasing, we can conclude that (26) can be positive or negative i.e. if we chose a control such that θ 1 < θ 0 the displacement will be positive instead if θ 1 > θ 0 it will be negative. In both cases the switching rule is respected thanks to the presence of the thermostat.

V w(t) (θ(t)) = V 2 (θ(t)) 0 < t < t 1 V 1 (θ(t)) t 1 < t < T. ( 27
)
where t 1 is the first time for which u goes through -ε

t 1 := inf {T > t > 0 | u(t) = -ε}
and T the final time. The displacement is

∆x = F 2 (θ(t 1 )) -F 2 (θ(0)) + F 1 (θ(T )) -F 1 (θ(t 1 )). ( 28
)
Calling again θ(t 1 ) := θ 1 also in the case we verify the surjectivity showing that the derivative of the displacement is different from zero.

Hence (28) is locally invertible and the inverse image of a neighborhood of 0 is a neighborhood of θ 0 . We can conclude as in the previous case that (28) can be either positive or negative choosing a suitable control.

Third case u(0) > ε and w 0 = 2,

V w(t) (θ(t)) =      V 2 (θ(t)) 0 < t < t 1 , V 1 (θ(t)) t 1 < t < t 2 , V 2 (θ(t)) t 2 < t < T. ( 29
)
with

t 1 := inf {T > t > 0 | u(t) = -ε} and t 2 := inf {T > t > t 1 | u(t) = ε}
The net motion is

∆x = F 2 (θ(t 1 ))+F 1 (θ(t 2 ))-F 1 (θ(t 1 ))-F 2 (θ(t 2 )). (30)
recalling that θ(0) = θ(T ). Also in this case we want to prove that we are able to move both forward or backward. Therefore we show that ∆x is surjective in ]0, π 2 [×]0, π 2 [ as in the non hysteretic case. We compute the gradient and show that it is never null

∇(∆x) = (V 2 -V 1 )(θ 1 ) (V 1 -V 2 )(θ 2 ) = 0 in (θ 1 , θ 2 ) ∈]0, π 2 [×]0, π 2 [
hence (30) is a submersion and surjective as required.

Notice that

∆x = (F 2 -F 1)(θ 1 ) -(F 2 -F 1 )(θ 2 )
and recall that the function (F 2 -F 1)(•) is always increasing. Hence, if we use a control such that θ 1 > θ 2 the ∆x will be positive, while if θ 1 < θ 2 then ∆x will be negative. Also in this case both the alternatives can be achieved respecting the switching rule. Therefore we are able to obtain the desired displacement.

Fourth case

The case where u(0) < -ε is analogous to the previous one.

In conclusion we have proved that for every choice of w 0 we are always able to find a periodic and continuous control θ = u that allows us to obtain the desired displacement. The system (24) is then partially controllable in x.

Remark 2

The introduction of the thermostat is essential because allows us to achieve displacements of every sign and thus the controllability result in x. This fact is strictly linked to the presence of the thresholds, indeed we are allowed to move between them without changing dynamics and therefore obtain values θ 1 < θ 2 either θ 1 > θ 2 , and thus move both forward and backward.

Remark 3 Note that the maximal value of r in the last theorem is |∆x(

π 2 )| if -ε < u(0) < ε, and |∆x( π 2 , 0)| if -ε < u(0) or u(0) > ε.
Thus it is always independent from x and θ. To cover distances |∆x| ≥ r we should divide the spatial interval in N subintervals of length less than r, each one realized by a u of period T N . Repeating N times this control u we are able to reach the desired displacement.

Global controllability result

In this subsection we are interested in studying whether it is feasible for the system of the scallop to move between two fixed configurations ((x 0 , θ 0 ) and (x f , θ f )). This part add something to the previous one, since we are prescribing both the initial and final positions and angles. The following holds: Theorem 5 Let A and B be two fixed positions along the x-axis and θ 0 , θ f two fixed angles. Then, we are always able to find a suitable control function u(t) such that the scallop system moves between A and B passing from θ(0) = θ 0 to θ(T ) = θ f , where T is a suitable big enough final time. Moreover such function u(t) respects the switching rules modeling the dependence of the viscosity η from | θ| (Fig. 2) and from sign( θ) with the thermostat, (Fig. 5). In other words the system (24) is controllable.

Proof Let u(t) the periodic function that makes the system move between A and B with final angle θ 0 during a time t . We have proved the existence of such a function with both switching rules, in the previous subsection. Now whatever w(t ) we open or close the valves respecting the switching rule in Fig. 2 or Fig. 5 respectively until we reach the desired angle θ f . We call t the time in which we have θ f and C the point in which we are arrived. Now starting from C with w(t ) we move to B using another periodic u(t) (hence θ(t) periodic), whose existence is ensured from Theorem 4 In this section we will show, through numerical simulations, that our theoretical pretictions on the controllability of the Scallop along x are good. Moreover we will also describe how it is possible to obtain the same results removing the continuity hypothesis on θ. In what follows the pictures are all relative to the controllability result which follows the thermostatic switching scheme (see Fig 5) that is the most interesting one. Similar results can be obtained analogously using the other switching described in Fig. 2.

Let us suppose to start with w(0) = 2 which means θ(0) > ε, the following pictures show a possible choice of the control θ to obtain a displacement ∆x = 1 cm, using the following parameters: a = 2 cm, b = 0.1 cm, η = 2 N sm -2 , ξ = 1 N sm -2 m = 1 g and ρ = 1 gcm -3 . More precisely in these simulations we decided to use a periodic polynomial control θ(t) that can be uniquely determined imposing the following constraints. Fig. 8 The lpolynomial control θ(t), the resulting periodic angle θ(t) and the corresponding x displacement in function of time.

It is easy to see that (since we want a positive displacement θ 1 > θ 2 ) θ respects the thermostatic switching rule and that after a time T = 7 s we have gained the desired displacement of 1 cm.

Starting from these simulations we want to build a piecewise constant control, instead of a continuous one, to obtain the same displacement. We note that in the case of delayed thermostat a discontinuous input is in general not allowed due to the presence of memory. The main difficulty of using a discontinuous control is to chose the switching times. Having in mind the previous simulations we can take the switching times of the continuous control and build a piecewise constant control which satisfies the constraints (31). Referring to the simulations in Fig. 8 Fig. 9 The piecewise constant control θ(t), the resulting angle θ and the corresponding x displacement in function of time.

These simulations actually prove that the displacement does not depend on the whole control trajectory but only on the values that the angle θ and its derivative θ assume in the switching times.

Conclusions

In this paper we analyze the system of a scallop proposing some strategies to overcome the famous scallop theorem. The main idea is to introduce a switching in the dynamics related to the variation of the angular velocity. This is done in two different ways, fact that helps to brake the reversibility of the equation of motion producing a net displacement. Original tool is also the introduction of the thermostat to model a delay in the change of fluid regime, and we show that it is crucial to gain both forward and backward motion. Moreover numerical simulations suggest also a way to use the switching schemes without necessarily using a continuous control input. Namely they show that it is possible to obtain the same displacement using a piecewise constant control (see Fig. 9) and choosing the switching times according to the ones used in the continuous case (see Fig. 8).

  e | θ| < M with M > 0) we use instead the viscous approximation and the relative dynamics.The switching rule in Fig 2 should also consider what happens when | θ| = M . However in the sequel we are going to exhibit a function θ which stays in M or -M for only a set of times of null measure.

Fig. 3

 3 Fig. 3 This figure shows a possible choice of θ(t) which realizes a positive displacement and respects the switching scheeme

Fig. 4

 4 Fig. 4 The rule of the classical switching

Fig. 6 A

 6 Fig. 6 A possible choice of θ(t) starting from 0 < u(0) < ε which realizes a positive displacement

Fig. 7

 7 Fig. 7 This figure represents one of the cases considered in the proof of Theorem 5

  θ(0) = θ0 θ(t 1 ) = -ε θ(t 2 ) = ε θ(T ) = θ0 θ(0) = θ 0 θ(t 1 ) = θ 1 θ(t 2 ) = θ 2 θ(T ) = θ 0 (31)where θ 1 and θ 2 are determined by the numerical inversion of the function ∆x (26) and we chose t 1 = 2 s, t 2 = 6 s and T = 7 s.

Acknowledgements The work has been developed within the OptHySYS project of the University of Trento that is gratefully acknowledged. Moreover we thank also Gruppo Nazionale Analisi Matematica Probabilitá e Applicazioni (GNAMPA) for partial financial support.

Compliance with Ethical Standards: Funding: this study was funded by University of Trento and Gruppo Nazionale Analisi Matematica Probabilitá e Applicazioni (GNAMPA).


Guillermo Gallardo

William Wells

Rachid Deriche

Demian Wassermann

Efficient Population-Representative Whole-Cortex Parcellation Based on Tractography

Introduction

The human brain is arranged in areas based on criteria such as cytoarchitecture or extrinsic connectivity. Current hypotheses attribute specialized functions to several areas of this patchwork. Hence, parcellating the cortex into such areas and characterizing their interaction is key to understanding brain function. Diffusion MRI enables the exploration of physical connections through axonal bundles, namely extrinsic connectivity. Current theories hold that brain function is determined by extrinsic connectivity. However, obtaining a population-representative parcellation based on extrinsic connectivity remains challenging [START_REF] Jbabdi | Long-range connectomics[END_REF]. Particularly, whole-cortex parcellation methods [START_REF] Moreno-Dominguez | A hierarchical method for whole-brain connectivity-based parcellation[END_REF][START_REF] Parisot | Information Processing in Medical Imaging[END_REF] are computationally expensive and need tuning of several parameters. Our main contribution is an efficient technique to create single-subject and population-representative parcellations based on tractography. Our method creates a dendrogram using only one parameter: the minimum size of each parcel. Then, by choosing cutting criteria, we can explore different parcellation granularities without recomputing the dendrogram. Experiments show that our extrinsic based parcellations are consistent within subjects with anatomical [START_REF] Desikan | An automated labeling system for subdividing the human cerebral cortex on MRI scans into gyral based regions of interest[END_REF] and functional [START_REF] Barch | Function in the human connectome: task-fMRI and individual differences in behavior[END_REF] parcellations existent in the literature.

Methods

We use 27 males aged 31-35 from the group S500 of the Human Connectome Project (HCP) processed with the minimum pipeline [START_REF] Glasser | The minimal preprocessing pipelines for the Human Connectome Project[END_REF].

For each subject we create tractograms using 15000 particles [START_REF] Behrens | Characterization and propagation of uncertainty in diffusionweighted MR imaging[END_REF]. We avoid superficial cortico-cortical fibers [START_REF] Reveley | Superficial white matter fiber systems impede detection of long-range cortical connections in diffusion MR tractography[END_REF] by shrinking each surface 3mm into the white matter. Once the tractogram of a seed is created, the voxel's value is the odds of being connected to the seed. To compare and average these odds values we transform them into a vectorial space using the logit function [START_REF] Pohl | Using the logarithm of odds to define a vector space on probabilistic atlases[END_REF]. This transform allows us to treat tractograms efficiently in a sound mathematical framework.

The proposed transformation allows us to cluster tractograms efficiently. We use as clustering algorithm a modified Agglomerative Hierarchical Clustering (AHC). Since we work in a vectorial space, we use the euclidean distance and the centroid as similarity and linkage functions, improving performance. Finally, to get local coherence [START_REF] Schmahmann | Fiber Pathways of the Brain[END_REF], we add a parameter to constraint the minimum size of the resulting clusters: clusters smaller than this size are only merged with neighbors, i.e. physically close clusters in the cortex. We create each subject's parcellation clustering all of his tractograms. Also, taking advantage of working in a vectorial space and seed correspondence across subjects, we create a population parcellation clustering the population-averaged tractograms of each seed (Figure 1).

Results

Figure 2 shows single-subject parcellations using a minimum cluster size of 1mm 2 . Figure 3 (a) shows 9 regions from an anatomical atlas [START_REF] Desikan | An automated labeling system for subdividing the human cerebral cortex on MRI scans into gyral based regions of interest[END_REF] projected over the subject 1 parcellation in Figure 2 and how many of our parcels were contained by those projections. Figure 3 (b) shows the projection of some of subject 1 parcels over functional responses to particular stimuli [START_REF] Barch | Function in the human connectome: task-fMRI and individual differences in behavior[END_REF]. Each color encodes the response to a different stimulus thresholded with a Z-score > 5. To compare them we calculate the Dice coefficient between our parcels and the activations. Several projected parcels have high coefficients with at least one functional response. Figure 4

Conclusions

Our pure extrinsic based parcellation has good agreement with anatomical and functional parcellations from the literature. Particularly, the motor and sensory cortex appear to be found. Having the tractograms in a vectorial space allowed us to work efficiently with them and to create a population-representative parcellation.

  (a) shows the same as figure3(a) for the population clustering of the 27 subjects. Figure4(a) shows the same as figure3(a) for the population clustering of our 27 subjects. Figure4(b) compares our population clustering with average functional responses in the Unrelated100 population from the HCP.

Figure 1 :

 1 Figure 1: Clustering Process: the tractograms are transformed into a vectorial space where we can manipulate and average them.

Figure 2 :

 2 Figure2: Results for the single-subject case. Some regions from[START_REF] Desikan | An automated labeling system for subdividing the human cerebral cortex on MRI scans into gyral based regions of interest[END_REF] appear to be found.

Figure 3 :

 3 Figure 3: Comparison between our automatic generated parcels; an anatomical atlas[START_REF] Desikan | An automated labeling system for subdividing the human cerebral cortex on MRI scans into gyral based regions of interest[END_REF]) and a functional study over the subject[START_REF] Barch | Function in the human connectome: task-fMRI and individual differences in behavior[END_REF]).

Figure 4 :

 4 Figure 4: Comparison between a population-representative parcellation of 27 subjects; an anatomical atlas (Desikan 2006) and a functional study over the Unrelated100 HCP population.


Valérie Pascual
email: valerie.pascual@inria.fr

Laurent Hascoët
email: laurent.hascoet@inria.fr

Mixed-language Automatic Differentiation

As AD usage is spreading to larger and more sophisticated applications, problems arise for codes that use several programming languages. Many AD tools have been designed with one application language in mind. Only a few use an internal representation that promotes language-independence, at least conceptually. When faced with the problem of building (with AD) the derivative code of a mixed-language application, end-users may consider using several AD tools, one per language. However, this leads to several problems:

• Different AD tools may implemented very different AD models such as overloading-based versus source-transformation based or association-by-address versus association-by-name. These models are often not compatible.

• When selecting the source-transformation model (for efficiency of the differentiated code), performance of the differentiated code strongly depends on the quality of data-flow analysis, which must be global on the code. A global analysis with separate AD tools would require inter-tool communication at the level of data-flow analysis, which does not exist at present.

 provides such an internal representation, accessible at present from C or Fortran sources.

Other AD tools provide a language-independent internal representation. OpenAD provides such a representation based on the XAIF formalism. However, this gives birth to two separate tools, OpenAD/F for Fortran, and ADIC2 for C. Still, it seems that their is no deep reason to prevent OpenAD application to mixed-language codes. We are lacking information about common architecture between TAF and TAC++ that would allow such mixed-language AD.

 was the first AD tool to support algorithmic differentiation in tangent mode of mixed-language components, specifically C++ and Fortran. As Rapsodia uses operator overloading, it performs no global analysis of the code. To our knowledge the extension of mixed-language differentiation with Rapsodia to adjoint mode is not yet provided.

Language standards and interoperability

Language standards often say very little about interoperability with other languages, leaving much freedom to compilers. Still, usage has progressively brought de facto standards, in particular between C and Fortran. The Fortran 2003 standard has specified its interaction with C in more detail. As far as AD is concerned, AD tools should in any case not commit to any specific interoperability strategy, and in particular to parameter-passing behaviors. Those might change with new languages and versions of languages. Instead, an AD tool must be able to handle a set of behaviors (hopefully small) from which one can describe all reasonable ways of parameter-passing.

The main issue raised by analysis and transformation of mixed-language codes is parameter-passing. Other issues are related to matching elements across languages, mainly types and procedures. These matching issues seem less complex than parameter-passing, but they should be discussed first.

Interoperability between types (and between variables of these types) across languages relies on identical memory representations built by compilers. Obviously interoperable types must match in the sense that they have the same structure, number of fields, and these fields must recursively be of interoperable types. Compilers often grant a natural interoperability between structured types. However Fortran 2003 provides the bind attribute to tell at compile time that a Fortran type has a C equivalent. It is essential to identify interoperable types in both languages, in particular because it may help distinguish candidate interoperable procedures according to the types of their arguments.

Similarly, interoperability between procedures relies on the compiler identifying the called procedure from another language, usually comparing procedure names, arguments number, and interoperable arguments types. Interoperable procedures in Fortran 2003 are declared with an explicit interface. The bind attribute allows one to associate the Fortran procedure name to the C procedure name. With Fortran 77 and Fortran 90, different conventions exist to associate both names: either by adding an underscore character at the end of the Fortran name, either with the same name, or with the name in uppercase, depending on the compilers.

The parameter-passing strategy is already a property of each given, single language. Mixed-language calls may be at the interface between two different parameter-passing strategies, which adds extra complexity. Still all resulting parameter-passing behaviors should be eventually expressed using a small number of elementary tactics.

Inside a given language, parameter-passing may use one of a few classical strategies: call by value, call by reference, call by value-result, call by sharing, call by name . . . . Call by value is the most common strategy. In call by value, the argument expression is evaluated, and the resulting value is copied to the corresponding variable in the function. If the function assigns values to its parameters, only its local copy is assigned and the argument passed into the function call is unchanged when the function returns. In call by reference, a procedure receives a reference to a variable and can modify the variable used as argument. Call by value-result, also named call by copy-restore, is a special case of call by reference. It differs from call by reference when two arguments alias one another. Under call by reference, writing to one will affect the other. Call by value-result gives the function distinct copies, but result in the callers environment depends on which of the aliased arguments is copied back first. Call by sharing is a terminology used by languages such as Python, Java and other object oriented languages. It is analogous to call by value where the passed value is either the argument when it is of a primitive type or its address when is is an object. Whereas in call by value, the arguments are computed before calling the function, in call by name, each occurrence of the formal argument is replaced with the actual argument in the style of macro-expansion.

Consider now the classical mix of Fortran and C. Parameter-passing mechanisms differ in the two languages. In Fortran, call by value-result and call by reference are used. Fortran 2003 introduces the VALUE attribute to specify call by value. In C, call by value is the only parameter-passing mechanism. All parameters are passed by value, except for arrays, which are translated into the address of the first element. In C, one simulates a call by reference by passing a pointer to this parameter. In mixed language calls, the caller and the called procedure must agree on how parameters are passed, e.g. passing a pointer from C to a Fortran procedure if the parameter has no VALUE attribute.

We believe that every mixed-language parameter-passing strategy can boil down to a few simple behaviors at the time of entry into and return from the called procedure. At call time, we define what we call the passed argument, which may be the same as the actual argument, or the memory pointed to by the actual argument, or conversely the address of the actual argument, depending on the mixed-language strategy that must be captured. Then the internal memory corresponding to this passed argument is copied into the internal memory corresponding to the called procedure's formal argument. At return time, the internal memory corresponding to the formal argument may be either copied back to the actual parameter, or not copied in which case it will vanish when the called procedure is popped from the call stack. When there is a back copy, it follows the link from the passed argument back to the actual argument: if they are the same, the copy is written into the actual argument, and if the passed argument is the destination of the actual, then the copy is written at the address designated by the actual argument.

Fig. 1 illustrates these behaviors for a few example multi-language calls, and also for pure Fortran calls distinguishing the scalar case from the array case and for a C call using pointers. For each situation, we give the choice of the passed argument and the choice about back copy that implement the desired behavior. Anticipating on the next section, Fig. 1 also shows for each situation, the Translator object used by Tapenade to specify these choices to the data-flow analyses. We will now investigate how to extend our AD tool to analyze mixed-language codes, and how differentiation must be adapted.

Extension of Tapenade algorithms for interoperability

Tapenade was originally designed to support different imperative languages. The motivation was to share the model of the tool and its implementation between these languages (at present, Fortran and C). We believe that this architecture also lets us deal with mixed-language source for a minimal implementation effort, affecting only a few components of the tool.

Tapenade represents a code as a call graph whose nodes represent procedures and arrows represent calls. Each call graph node contains a flow graph, in which nodes are blocks of elementary instructions (in particular calls), and arrows represent control jumps. At present, Tapenade called on C source or on Fortran source builds an internal representation of the same nature, using the same components for procedures, instructions, variables, types . . . We exploit this unicity of representation to deal with mixed-language codes, in particular mixing Fortran and C.

Let's first take a look at matching of types and procedures. Two interoperable entities are represented with the same internal representation in the AD tool. This representation distinguishes each component of structured types, and distinguishes pointer variables from their pointee destination variables. Representation of arrays on the other hand do not distinguish array elements, and therefore we need not worry about their possibly different memory layout in Fortran and C. Sometimes it is not possible to preserve this nice structural matching, for instance for the complex Fortran type which should match a 2-fields structure in C. Then the correspondence must be enforced by implementation, and possibly incurs some loss of information.

The question of procedure matching amounts to finding, for a given procedure call, the node of the call graph that will be effectively called. This depends on the mixed-language conventions on procedure names and on type matching. Procedure name conventions may vary, and Tapenade offer parameterization to define one convention, using command-line arguments or directives at the call site. It also interprets attributes of Fortran 2003. For instance, these parameters let us specify that inside a Fortran code all calls to a procedure named FOO will connect to a C procedure named foo .

Fortran calls Fortran, scalars

Fortran calls Fortran, arrays Parameter passing comes into play during data-flow analysis. All data-flow analyses (e.g. in-out, activity, liveness) inherit from a base analysis class that provides primitives to transfer data-flow information between a caller procedure and a callee. Most adaptions to mixed-language code must be done in this base class. This information transfer is driven by an object we call a Translator, which describes how actual arguments are matched with formal arguments. Variables, and in particular those in arguments, may have a finer structure. In addition to the classical scalar and arrays that existed for Fortran 77, languages have introduced variables of structured type and pointers. Since the data-flow properties that we analyze may be different for each component of these structured objects, we distinguish all of their elementary components. For example, a formal argument of type mystruct *arg where mystruct is a record: struct mystruct { int numElems; float *elems; } is represented by 4 elementary formal arguments, for *(arg->elems), arg->elems, arg->numElems, and the toplevel arg respectively. To each elementary formal argument of the called procedure, the Translator associates the corresponding elementary actual argument at the call site, which is either an expression or an elementary variable known to the calling procedure. In addition, the Translator associates to each elementary formal argument a boolean that specifies whether it must be copied upon return. For each example situation in Fig. 1, the Translator that implements the desired behavior is shown below the textual description of the behavior, as a set of arrows from formal elementary argument to actual elementary argument and back copy boolean. The rule of thumb is that the Translator associates the root formal argument with the passed argument.

REAL X SUBROUTINE GEE(V) . . . REAL V CALL GEE(X) -Passed argument is X -Upon return, value of V is copied back into X Translator: V -> X (Back copy) REAL Y( 1 0 0 ) SUBROUTINE GEE(B) . . . REAL B( 2
The "Fortran calls C" situation deserves further comment: since the C formal argument is a pointer to a float, there are in fact two elementary formal arguments, one for a and one for *a. The same happens for "C calls C". As a is associated with the passed argument &X, *a is naturally associated with *(&X) in other words with X. It is in fact thought this second elementary argument that the actual value of (or information regarding) X is propagated to the callee. Consequently, even if no back copy is done upon return into &X itself, every write into *a in foo is automatically reflected into X.

The way each data-flow analysis, which computes a given data-flow value, uses the specification from the Translator can be sketched as follows: immediately before the call, and for each elementary formal argument (left column of Translator), we retrieve the corresponding elementary actual argument (right column), and we retrieve the current data-flow value for it. The initial data-flow value of the elementary formal argument is set to this retrieved value. Analysis can then run on the called procedure. Upon return from the call, the top-level elementary arguments that bear the "No back copy" retain the data-flow value they had before the call. For all other elementary arguments the data-flow value of the elementary formal argument upon return is back copied into the data-flow value of the elementary actual argument. This description applies to forward data-flow analyses. Adaption to backward analyses require only minor technical changes. The main impact of dealing with mixed-language code on actual differentiation is related to arguments being copied to the formal arguments and not being copied back. For adjoint differentiation, the consequence is that the copied formal argument must not be identified with the actual argument of the call. In particular, the adjoint formal argument must be a different variable from the adjoint actual argument, and the adjoint code must follow the pattern illustrated in Fig. 2: variable u is the original formal argument, which is not copied back because of its "value" attribute. The differentiated actual argument is (conceptually) ub0, whereas the differentiated formal argument is ub. Only at the end of BAR B, the differentiated formal argument ub is added into the differentiated actual argument ub0. Consequently, ub0 must be passed a reference xb, even if the corresponding non-differentiated argument was passed a dereferenced *x. In other words the adjoint of a passed by value argument must sometimes be passed by reference.

When differentiating a procedure that calls a function F that returns a value, the last argument of the differentiated procedure contains the result of F. This argument and its corresponding passed parameter must be declared and used according to the parameter-passing mechanism, e.g. by passing the address of the argument when calling a Fortran function from C, instead of just passing the argument for a C function.

Conclusion and further work

Tapenade is now able to differentiate codes that mix C and Fortran, with calls in either direction. Apart from short non-regression tests, we are validating this extension on the "CalculiX" finite element library. The architecture of Tapenade allowed us to implement this functionality quite easily, with modifications in only a few component of the tool. We still need to study interoperability of global variables, e.g. using the bind statement of Fortran 2003. As the implementation is still at an early stage, we need to experiment on a few representative examples.

-

  Passed argument is *y -Upon return, value of V is copied back into *y Translator: V -> *y (Back copy) f l o a t y ; SUBROUTINE BAR(V) . . . REAL, VALUE : : V bar ( y ) ; -Passed argument is y -Upon return, no copy takes place into y Translator: V -> y (No back copy) Fortran calls C C calls C REAL X void f o o ( f l o a t * a ) . . . CALL FOO(X) -Passed argument is address of X -Upon return, no copy takes place into &X Translator: a -> &X (No back copy) *a -> X (Back copy) f l o a t * y ; void bar ( f l o a t * a ) . . . bar ( y ) ; -Passed argument is y -Upon return, no copy takes place into y Translator: a -> y (No back copy) *a -> *y (Back copy)

-Figure 1 :

 1 Figure 1: Mixed-language calls and the Translator that implements their behaviors

  void bar ( f l o a t a , f l o a t * b ) ; void f o o ( f l o a t * x , f l o a t * y ) { bar ( * x , y ) ; } void f o o b ( f l o a t * x , f l o a t * xb , f l o a t * y , f l o a t * yb ) { b a r b ( * x , xb , y , yb ) ; SUBROUTINE BAR( u , v ) BIND(C) IMPLICIT NONE REAL, VALUE : : u REAL : : v u = 2 * u v = u * u END SUBROUTINE SUBROUTINE BAR B( u , ub0 , v , vb ) BIND( c ) IMPLICIT NONE REAL, VALUE : : u REAL : : ub , ub0 REAL : : v REAL : : vb u = 2 * u ub = 2 * u * vb * vb = 0 . 0 ub = 2 * ub ub0 = ub0 + ub END SUBROUTINE BAR B

Figure 2 :

 2 Figure 2: Mixed-language differentiation with Tapenade, C calling Fortran case


Mohamad Taki Asghar
email: mohamadtaki.asghar@arcelormittal.com

Marc Jungers
email: marc.jungers@univlorraine.fr

Irinel-Constantin Morarescu
email: constantin.morarescu@univlorraine.fr

Ahmed Khelassi
email: ahmed.khelassi@arcelormittal.com

Julien Francken
email: julien.francken@arcelormittal.com

Mohamad Taki Asghar
email: mohamadtaki.asghar@univlorraine.fr

Irinel-Constantin Morȃrescu

Tandem cold rolling mill modeling for multi-variable control synthesis

Keywords: Tandem Cold Rolling, Continuous Tandem, Mathematical Model, State Space Representation, differential-algebraic system

come   L'archive ouverte pluridisciplinaire

INTRODUCTION

Metal rolling is a metallurgic process in which the thickness of a metallic strip is reduced. This process is important since the manufacturing of a wide types of products requires various thin and smooth metal sheets. Tandem Cold Mill (TCM) [START_REF] Roberts | Cold Rolling of Steel. Manufacturing Engineering and Materials Processing[END_REF] consists of consecutive stands where each one is usually a two high mill stand, which means that the stand contains a pair of work rolls. Writing down the mathematical formulation of all the subparts of TCM, we end up with a system having a large dimension consisting of non-linear differential equations and non-linear algebraic constraints. Moreover, the transportation of the strip with variable speed from one stand to the consecutive one implies the presence of state-dependent delays.

Due to the evolution of the technology, the specifications required for the steel sheets continuously change asking for thinner sheets with an increased performance assessment (mainly hardness, but also flatness...) (Jelali, 2007). In order to accomplish this objective one needs a grasp understanding of the process allowing the use of advanced control techniques, that can take into account the forthcoming introduction of possible new actuators. Moreover ArcelorMittal produces annually around 90 millions tonnes of crude steel. Any improvement of the global production cost by control theory, that is without changing the device, will induce significantly positive fallouts and will be at the heart of international competitiveness.

Due to its wide applications, the rolling process attracted the interest of many researchers and different mathematical models This work has been supported by a grant of the French Agency for Research and Technology (ANRT).

have been proposed [START_REF] Orowan | The calculation of roll pressure in hot and cold flat rolling[END_REF]Bland andFord, 1948, 1952;Bryant, 1973;[START_REF] Roberts | Cold Rolling of Steel. Manufacturing Engineering and Materials Processing[END_REF]. As detailed in the next section the modeling of complex systems occurring in the rolling process received also a lot of attention (see for instance [START_REF] Pietrzyk | Thermal-mechanical modeling of the flat rolling process[END_REF] for the hardness variation model or [START_REF] Sims | Calculation of roll force and torque in hot rolling mills[END_REF] for the force-mill characteristics).

Even if the control theory presents for industrial applications important evolutions during the last decades, as reported in the literature [START_REF] O'dwyer | Pi and pid controller tuning rules for time delay process: a summary[END_REF][START_REF] Silva | PID Controllers for Time Delay Systems[END_REF] more than 95% of the control-loop in the industry are controlled by (a set of) single input single output (SISO) Proportional-Integral-Derivative (PID) controllers. That can be explained by their ease of implementation and tuning. This observation holds true for the rolling process. Nevertheless, the evolution of the required steel sheet specifications asks for performances beyond the capacities of such a controller. It is thus necessary to adopt new control strategies aiming at: reaching the required specifications of the steel sheets; taking into account the complex interactions inside the system; reducing the energy cost of the rolling process and finally being ready for the further introduction of new actuators such as flexible lubrication [START_REF] Laugier | Flexible lubrification concept, the future of cold rolling lubrification[END_REF][START_REF] Alves | Flexible lubrification for controlling friction in cold rolling, crucial to be successful for the AHSS challenge[END_REF].We are entering together this way by developing new multivariable controllers for TCM.

Even if several techniques are dedicated to complex non-linear systems and have been applied on TCM (see [START_REF] Pittner | Tandem Cold Metal Rolling Mill Control[END_REF]), a larger class of flexible tools is related to multivariable linear systems [START_REF] Skogestad | Multivariable Feedback Control: Analysis and Design[END_REF]. In this paper, we propose to provide a rigorous linearization methodology taking into account the interstands transport delays to obtain an appropriate model for the design of multivari-able controllers. The proposed methodology being generic can be adapted to TCM with new additional actuators.

Our main motivation is given by the closed-loop performance improvements in presence of variable friction due to the flexible lubrication added in the process. The main contributions of this paper are thus the following:

(1) we establish and linearize a Multi Input Multi Output (MIMO) model that takes in consideration a variable friction for the overall TCM system; (2) we validate the linearized model on appropriate variation intervals of various quantities occurring into dynamics;

(3) we propose a preliminary multi-variable controller in absence of a priori fixed performance objectives.

The remaining part of the paper is organized as follows. Section 2 provides a detailed presentation of the subsystems consisting the dynamics of each stand of the TCM. The mathematical description of the overall TCM model as well as its linearization is presented in Section 3. This section also provides a numerical validation of the linearized system on appropriate variation intervals of the physical quantities occurring in the process. A preliminary multi-variable control scheme is designed and illustrated in Section 3. The paper finishes with some concluding remarks and perspectives.

TANDEM COLD MILL MODELING

The cold rolling tandem mill is typically made up of five stands, associated in cascade, traversed by the product steel strip which passes through their independent driven work rolls.

The objective of such process is to reduce strip thickness and to obtain a desired product having certain thickness, rigidity and flatness. The structure of these stands is similar, thus the description of only one of them will be given in this section. The whole model concatenates the models of each part of the stand and the actuators, among them, several are phenomenological approximations and then not unique. First of all, let us present a block diagram representing the different functions and components of the chosen model representing the behavior of one mill stand in Figure 1. In the following, each block is briefly presented.

Stand Motor: Most of the literature (see (Geddes, 1998;Pittner andSimaan, 2010, 2011;[START_REF] Alves | Flexible lubrification for controlling friction in cold rolling, crucial to be successful for the AHSS challenge[END_REF][START_REF] Alves | Modeling, simulation and identification for control of tandem cold metal rolling[END_REF]) on the modeling of TCM assumes that the stand motor is a first order system. In this work, the motor is modeled as a second order to reflect the presence of an internal pre-controller.

Stand Actuator: This block is a high order non-linear system describing typically two main parts of the actuator: the servo valve and the hydraulic cylinder. An identified second order model was used to approximate and describe mathematically this actuator. Yield Stress -SMATCH Model: Yield stress (flow stress, deformation resistance) model reveals the hardness variation of material during rolling process. The reduction of strip thickness after each stand increases the yield stress quantity.

The relation describing the yield stress with the strain and the strain rate is complex and non-linear. It may be algebraic or dynamic. Several approximations more or less sophisticated are available in the literature. Among them, several models are obtained via implemented statistical methods interpolating the flow stress in function of the strain.

In this paper, SMATCH law (Ferreira et al., 2008;Dbouk et al., 2013Dbouk et al., , 2014) ) is considered in the process modeling and could be described as follows:

k(ε p ) = λ 0 + (λ 1 + λ 2 ε p ) 1 -λ 3 e -λ4εp , (1)
where λ 0 , λ 1 , λ 2 , λ 3 and λ 4 are constants determined statistically for each specific steel, k(•) is the equivalent yield stress and ε p is the equivalent plastic strain. Mass Flow Conservation: This law is simply conserving material volume (mass flow) between the input and output of each mill stand. Assuming that the strip width does not change, the product of strip thickness and speed is conserved during the rolling process. Forward Slip: The rolled strip undergoes backward and forward slips respectively at the entry and exit of each roll bite. The backward slip is neglected while the forward slip at the exit is taken into consideration. The forward slip is the difference between the strip exit speed and work roll speed divided by the work roll speed. Sims' Gauge Model: Based on the experimental curve of force-mill stretch characteristics exhibiting a non-linear part for low forces and a linear one for high forces, an affine approximation is given to estimate the strip exit thickness in function of the work roll gap variation and the total rolling force [START_REF] Roberts | Cold Rolling of Steel. Manufacturing Engineering and Materials Processing[END_REF][START_REF] Sims | Calculation of roll force and torque in hot rolling mills[END_REF]. Notice that the slope depends on the mill modulus. Inter-stand Tension Stress: The passage of strip from one stand to another, accompanied by the variation of strip speeds at the entry and exit of each stand leads to the variation of the tension stress between each two consecutive stands. The tension stress is given by Hook's law applied to the strip. The elasticity of the strip depends on Young's Modulus for steel. Note that the non-linearity of both entry and exit speed models leads to strongly non-linear model of tension stress. Rolling Mathematical Model -(Roll-Gap): Through the various research in the domain of cold rolling during the past few decades, several models where given and developed for expressing the required rolling force and torque in certain analytical models. Many physical quantities (tension stress, flow stress, strip speed, thickness,• • • ) at the entry and the exit of each stand intervene in this model as seen in Figure 2.

Several studies tried to reduce the complexity of Orowan model [START_REF] Orowan | The calculation of roll pressure in hot and cold flat rolling[END_REF][START_REF] Pietrzyk | Thermal-mechanical modeling of the flat rolling process[END_REF] by introducing some assumptions for cold rolling process. As a result, more tractable models such as Bland andFord (Bland andFord, 1948, 1952;Bland and Sims, 1953), Bryant and Osborn (Bryant, 1973), W. Roberts [START_REF] Roberts | Cold Rolling of Steel. Manufacturing Engineering and Materials Processing[END_REF], A. Tselikov [START_REF] Pietrzyk | Thermal-mechanical modeling of the flat rolling process[END_REF] were proposed. to regulate is the exit thickness at the exit of each stand, because it is very sensitive and a source of perturbation for the whole system. To better regulate this variable we add an integral of the exit thickness. That will allow the multivariable controller to have an integral action, similar to the integral term of the PID and to vanish the static position error.

TCM MULTIVARIABLE MODEL

Non-linear Form

As it has been shown, TCM consists of several physical systems describing all the process. Some of them are strongly non-linear especially for inter-stand tension stress and Roll-Gap models.

Therefore by collecting all these dynamics, we can derive a state space non-linear multi-variable representation describing the whole TCM dynamics. This representation encapsulates all internal interactions and couplings between sub-systems. The obtained global model will depend on differential and algebraic states, input and output vectors as follows.

d dt x(t) = F (x(t), u(t), z(t), µ(t), d(t), h en (t)) , 0 = G (x(t), u(t), z(t), µ(t), d(t), h en (t)) , y(t) = L (x(t), u(t), z(t), µ(t), d(t), h en (t)) , (2)
where the control vector u(t) ∈ R 10 , including motor and actuator inputs. The differential state vector x(t) ∈ R 29 , which contains dynamics of inter-stand tension stress, motor and actuator dynamics, exit thickness integral actions. The algebraic state vector z(t) ∈ R 15 includes the rolling forces, rolling torques and forward slips. The exogenous vector d(t) ∈ R 3 includes the entry annealed thickness, and TCM entry and exit tension stresses. The output vector y(t) ∈ R 50 , which includes all needed physical quantities in the TCM process such as tension forces, rolling forces, work roll speeds... etc. Entry thickness vector h en (t) ∈ R 4 is the vector representing the entry thicknesses for all stands except first one, which is an exogenous input. Friction vector µ(t) ∈ R 5 , which includes the frictions between strip and work rolls in all stands.

The non-linear function F in the first equation represents the evolution of the system differential dynamics, while the nonlinear function G in the second one represents the algebraic constraints of the model. The last function L represents the system outputs.

Modeling the TCM process aims at approaching as closely as possible the behavior of the real tandem. The class of nonlinear multi-variable models has a strong capacity to represent high order and complex industrial processes. Nevertheless, dedicated non-linear control strategy may imply very sophisticated tools, like State-Dependent Riccati Equations (SDRE) [START_REF] Pittner | Tandem Cold Metal Rolling Mill Control[END_REF], sliding modes or flatness. Here our approach is to linearize the TCM non-linear model by taking into account the presence of state-dependent delays to permit the use of the wide range of linear control tools. This method will allow to reach a compromise between good approximation and availability of tools. Linearization will permit to consider one of the linear multi-variable control strategies dedicated to such types of systems containing state and time dependent delays. Furthermore these strategies are robust and can be adjusted to certain industrial specifications booklet.

Linearized Form

In this paper, the linearization will be carried out in the neighborhood of a nominal operating point which is the steady state of the non-linear differential-algebraic system given by Equation (2). In other words, the nominal operating point is a solution of the following set of algebraic equations. 0 = F (x eq (t), u eq (t), z eq (t), µ eq (t), d eq (t), h en,eq (t)) , 0 = G (x eq (t), u eq (t), z eq (t), µ eq (t), d eq (t), h en,eq (t)) , y eq (t) = L (x eq (t), u eq (t), z eq (t), µ eq (t), d eq (t), h en,eq (t)) .

(3)

Solving the above equations may have many solutions. The user in rolling process precises the nominal operating values of some physical quantities (inter-stand tension stress between stands, exit thickness in each stand, forward slips and others) for each rolled product and thus this limits the number of solutions. The solution depending on the rolling scenario is defined for certain product.

During rolling process and due to several internal and external perturbations, the system deviates from its rolling operating point. In the vicinity of the steady state, the linearized form of the system governs the evolution of these deviations, which should be rejected by the controller. The linearization of the set of Equations ( 2) consists in obtaining the expansion of Taylor series of first order at the operating point.

The validation of the linearization answers if the size of the vicinity where the difference between the non-linear model and the linearized one is enough small for the application.

The linearization phase of the non-linear Roll-Gap model is numerically validated by varying quasi-statically one input of the Roll-Gap algorithm and fixing the others one. It is then possible to compute the partial derivative by the rate of differences.

Figure 3 shows the non-linear variation of the non-linear specific rolling force with respect to the exit thickness (blue curve) around its nominal value for an arbitrary stand i. The straight line is the linear approximation of the non-linear model. As it is clear, for a large variation of the exit thickness (up to ±8% of the nominal value), the linear approximation (black line) remains acceptable with a small error less than 3% at the extremities. This means that the Roll-Gap linear model can be considered as enough correct. To go deeper into the validation, the same method is applying by varying quasi-statically two inputs of the Roll-Gap algorithm and fixing the others ones. Figure 4 shows the variation of the specific rolling force with respect to both entries, exit thickness (up to ±8%) and friction (up to ±20%). The non-linear surface (blue) shows slight curvilinear aspect at the extremities.

Friction varies in a wide range, since in cold rolling there is no determined variation of friction between work rolls and the strip. The plane representing the linear surface (green) shows a very good approximation of the non-linear one (multi-colored).

In order to better evaluate and validate this approximation, it is important to introduce the gap between both surfaces. Figure 5 shows that the difference is less than ±4% with high friction variation (up to ±20%). The error is small comparing to the inputs wide variations.

As it has been precised, the linearized model will be operated around its steady state operating point, and any variation or deviation will take place in its neighborhood. So it is better and easier to express any quantity X which have an equilibrium value X eq by its variation value denoted ∆X = X -X eq . This variation is due to the external perturbation and variation in the exogenous components, the model will deviate relative to its equilibrium point. See Equations ( 5)-( 7). The entry thickness of the each stand h en,i is the delayed exit thickness of the previous stand. ∆h en,i+1 (t) = ∆h ex,i (t -τ i,i+1 ),

where h ex (•) is the exit thickness vector containing all exit thicknesses h ex,i (•). τ i,i+1 represents the time delay between stand i and i + 1. The exit thickness is delayed due to the transportation of the steel over the inter-stand distance and then enters the next stand. Thus a part of the vector y will be delayed and looped as an additional input of the TCM system as shown in Figure 6. (7)

1978), the exit thickness is function of work roll position (differential dynamic in vector x) and total rolling force (algebraic dynamic in vector z). So, the entry Thickness vector ∆h en can be expressed by certain combination of delayed state differential and algebraic vectors.

∆h en (t) = 4 i=1 Λ x,i x(t -τ i,i+1 ) + Λ z,i z(t -τ i,i+1 ), (8)
where Λ x,i and Λ z,i are constant matrices of appropriate dimensions. Replacing Equation ( 8) in ( 5)-( 7) and regrouping all matrices in Equation ( 5)-( 7), the obtained form is given by Equations ( 9)-( 11), where τ 0,1 = 0, Ãi , Bz,i , Bd,i , Bµ,i , Ci , Dd,i , Dµ,i are matrices of appropriate dimensions. This form introduces the state dependent delay effect on all physical quantities. Note that in other studies (Pittner andSimaan, 2010, 2011;[START_REF] Alves | Flexible lubrification for controlling friction in cold rolling, crucial to be successful for the AHSS challenge[END_REF][START_REF] Alves | Modeling, simulation and identification for control of tandem cold metal rolling[END_REF] delays are usually approximated by transfer functions following the padé approximation or multiple first order lags, by assuming that the delay is constant to allow a frequency domain representation. Here such an assumption is not consider to emphasize the fact that the delays are state-dependent.

Simplistic Multi-variable Control

In order to highlight the efficiency of the model and also to outline the perspectives of our research, a simulation is presented with a typical variation of the strip entry thickness at the first stand. To have a sense, it is nevertheless necessary to associate a control law with the obtained model. A simplistic Linear Quadratic Regulator (LQR) (Anderson and Moore, 1971) is chosen for illustration purposes and to show the feasibility of our method. The LQR strategy leads to a state-feedback controller solution of an optimization problem with a quadratic cost function (specifying the requirements) under the linear constraint characterizing the system. It is well known by its well robustness and performance. The class of controllers and their tunings will be tackled during the following steps of our research. In that conditions, the following simplifying assumptions are imposed:

• the full state of the system is accessible;

• the weights of quadratic cost function are fixed to very simple values; • the linear system is controllable; in fact, this is verified numerically on the linear system; • only for the LQR synthesis, the delays are imposed to be trivial.

The obtained LQR is applied on the linearized system with (non trivial) state-dependent delays. The exogenous signals vary near their nominal values, the annealed entry thickness (strip entry thickness at the first stand) varies about 8%. The TCM entry traction stress varies around 2% to 3%, while the TCM exit traction stress varies about 1%. The friction in each stand varies as a noise around a certain value. Figure (7) shows the normalized variation of the exit thickness of each stand, with respect to the steady exit thickness. The percentage error at the exit of the first four stands is less than 2%, while at the tandem exit is less than 1%. The delay effect is clear in the strip propagation within the presented figures. This controller synthesis did not take into account the presence of inter-stand delays, though it shows a good feasibility when tested on the delayed system. The control of TCM is going to be the main objective in the future work. It should consider all system complexities and different aspects.

All delays are variable, same as friction in stands are not constant. Thus the control method has to be robust enough and should follow certain specifications. Different MIMO control methods are available such as robust control (Geddes and [START_REF] Postlethwaite | Gauge control incontrol in tandem cold rolling mills: a multivariable case study using H ∞ optimization[END_REF][START_REF] Postlethwaite | Gauge control incontrol in tandem cold rolling mills: a multivariable case study using H ∞ optimization[END_REF][START_REF] Pittner | Tandem Cold Metal Rolling Mill Control[END_REF].

CONCLUSION

The present work has proposed a global mathematical modeling of the TCM process. The elaborated MIMO model is strongly non-linear with state dependent delays. Moreover this modeling takes into account stands' friction variation in the process during rolling. It is linearized around an operating point depending on the rolled product. The linearization was validated for sufficient variation around the steady state. The linearized model is completed by the presence of state-dependent delays. Future works will focus on the design of MIMO control strategies that consider robustness, performance industrial specifications and also variable frictions.

Fig. 1 .

 1 Fig. 1. Mill Stand Block Diagram

Fig. 2

 2 Fig. 2. Rolling Bite ArcelorMittal Maizières with the help of academic partners has developed the so called Roll-Gap Model for cold rolling of steel which computes numerically the rolling force, torque and forward slip for a given rolling situation (Dbouk et al., 2013). Moreover it expresses high complexity and nonlinearity. It is based on Slab method considering sever elastic deformation of work rolls. As a contribution in (Hitchcock, 1935), the flattened work roll radius was assumed to be elliptical. In Roll-Gap model, influence function technique is used to calculate the radial deformation of work roll radius due to rectangular compressive stresses. Inside the rolling zone, Roll-Gap model decomposes the rolled strip and work roll cylinder into finite elements and study the group of stresses exerted on each element. This permits to calculate precisely the elastic deformation of the work rolls and the plastic deformation of the strip. This model includes convergence algorithm with relaxation to validate the obtained values. Roll-Gap has been validated via large number of simulations and practical experiments. It has been simulated and compared to basic rolling algorithms. Inter-stand Delay: The transportation of the rolled strip with variable speed from one stand to the consecutive one induces the transport delays, that are state-dependent and vary in time. Exit Thickness Integral Action: One of the main variablesto regulate is the exit thickness at the exit of each stand, because it is very sensitive and a source of perturbation for the whole system. To better regulate this variable we add an integral of the exit thickness. That will allow the multivariable controller to have an integral action, similar to the integral term of the PID and to vanish the static position error.

  Fig. 3. Comparison between Specific Force Roll Gap and its linearized model in function of exit thickness.

Fig. 4 .

 4 Fig. 4. Comparison between specific force and its linearized model in function of exit thickness and friction.

Fig. 5 .

 5 Fig. 5. Error between non-linear and linearized specific force surfaces.

Fig. 6 .

 6 Fig. 6. TCM Linearized System looped by delayed Exit Thickness

Fig. 7 .

 7 Fig. 7. TCM stands' exit thicknesses ∆h ex,i (t) with a simplistic LQR

  A∆x(t) +B z ∆z(t) + B µ ∆µ(t) +B d ∆d(t) + B h ∆h en (t) +B u ∆u(t), (5) 0 = H x ∆x(t)+H z ∆z(t) + H µ ∆µ(t)+H d ∆d(t) + H h ∆h en (t), (6) ∆y(t) = C∆x(t) +D z ∆z(t) + D µ ∆µ(t)+D d ∆d(t) + D h ∆h en (t).

	d dt	∆x(t) =
		Roberts,