A Terrand-Jeanne

V Andrieu

V Dos

Santos Martins

C.-Z Xu

Adding integral action for open-loop exponentially stable semigroups and application to boundary control of PDE systems

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

Adding integral action for open-loop exponentially stable semigroups and application to boundary control of PDE systems

The use of integral action to achieve output regulation and cancel constant disturbances for infinite dimensional systems has been initiated by S. Pojohlainen in [START_REF] Pohjolainen | Robust multivariable pi-controller for infinite dimensional systems[END_REF]. It has been extended in a series of papers by the same author (see [START_REF] Pohjolainen | Robust controller for systems with exponentially stable strongly continuous semigroups[END_REF] for instance) and some other (see [START_REF] Xu | A robust pi-controller for infinitedimensional systems[END_REF]) always considering bounded control operator and following a spectral approach (see also [START_REF] Paunonen | Internal model theory for distributed parameter systems[END_REF]).

In the last two decades, Lyapunov approaches have allowed to consider a large class of boundary control problems (see for instance [START_REF] Bastin | Stability and boundary stabilization of 1-d hyperbolic systems[END_REF]). In this work our aim is to follow a Lyapunov approach to solve an output regulation problem. The results are separated into two parts.

In a first part, abstract Cauchy problems are considered. It is shown how a Lyapunov functional can be constructed for a linear system in closed loop with an integral controller when some bounds are assumed on the control operator and for an admissible measurement operator. This gives an alternative proof to the results of S. Pojohlainen in [START_REF] Pohjolainen | Robust multivariable pi-controller for infinite dimensional systems[END_REF] (and [START_REF] Xu | A robust pi-controller for infinitedimensional systems[END_REF]). It allows also to give explicit value to the integral gain that solves the output regulation problem.

In a second part, following the same Lyapunov functional design procedure, we consider a boundary regulation problem for a class of hyperbolic PDE systems. This result generalizes many others which have been obtained so far in the regulation of PDE hyperbolic systems (see for instance [START_REF] V Dos | Boundary control with integral action for hyperbolic systems of conservation laws: Stability and experiments[END_REF], [START_REF] Xu | Multivariable boundary pi control and regulation of a fluid flow system[END_REF], [START_REF] Bastin | Stability of linear density-flow hyperbolic systems under pi boundary control[END_REF], [START_REF] Trinh | Design of integral controllers for nonlinear systems governed by scalar hyperbolic partial differential equations[END_REF], [START_REF] Trinh | Stability and output regulation for a cascaded network of 2x2 hyperbolic systems with pi control[END_REF], [START_REF] Terrand-Jeanne | Regulation of the downside angular velocity of a drilling string with a p-i controller[END_REF]).

The paper is organized as follows. Section II is devoted to the regulation of the measured output for stable abstract Cauchy problems. It is given a general procedure, for an exponentially stable semigroup in open-loop, to construct a Lyapunov functional for the closed loop system obtained with an integral controller. Inspired by this procedure, the case All authors are with LAGEP-CNRS, Université Claude Bernard Lyon1, Université de Lyon, Domaine Universitaire de la Doua, 43 bd du 11 Novembre 1918, 69622 Villeurbanne Cedex, France. of boundary regulation is considered for a general class of hyperbolic PDE systems in Section III. The proof of the theorem obtained in the context of hyperbolic systems is given in Section IV.

This paper is an extended version of the paper presented in [START_REF] Terrand-Jeanne | Lyapunov functionals for output regulation of exponentially stable semigroups via integral action and application to a hyperbolic systems[END_REF]. Compare to this preliminary version all proofs are given and moreover more general classes of hyperbolic systems are considered.

Notation: subscripts t, s, tt, . . . denote the first or second derivative w.r.t. the variable t or s. For an integer n, I dn is the identity matrix in R n×n . Given an operator A over a Hilbert space, A * denotes the adjoint operator. D n is the set of diagonal matrix in R n×n .

II. GENERAL ABSTRACT CAUCHY PROBLEMS A. Problem description

Let X be a Hilbert space with scalar product , X and A : D(A) ⊂ X → X be the infinitesimal generator of a C 0semigroup denoted t → e At . Let B and C be linear operators, B from R m to X and C from D(C) ⊆ X to R m .

In this section, we consider the controlled Cauchy problem with output Σ(A, B, C) in Kalman form, as follows

ϕ t = Aϕ + Bu + w , y = Cϕ, (1)
where w ∈ X is an unknown constant vector and u : R + → R m is the controlled input. We consider the following exponential stability property for the operator A.

Assumption 1 (Exponential Stability):

The operator A generates a C 0 -semigroup which is exponentially stable. In other words, there exist ν and k both positive constants such that,

∀ϕ 0 ∈ X and t ∈ R + e At ϕ 0 X k exp(-νt) ϕ 0 X . (2)
We are interested in the regulation problem. More precisely we are concerned with the problem of regulation of the output y via the integral control

u = k i K i z , z t = y -y ref , (3)
where

y ref ∈ R m is a prescribed reference, z ∈ R m , K i ∈
R m×m is a full rank matrix and k i a positive real number. The control law being dynamical, the state space has been extended. Considering the system Σ(A, B, C) in closed loop, with integral control law given in (3) the state space is now X e = X × R m which is a Hilbert space with inner product

ϕ ea , ϕ eb Xe = ϕ a , ϕ b X + z a z b , where ϕ ea = ϕ a z a and ϕ eb = ϕ b z b . The associated norm is denoted • Xe . Let A e : (D(A) ∩ D(C)) × R m
→ X e be the extended operator defined as

A e = A BK i k i C 0 . (4)
The regulation problem to solve can be rephrased as the following.

Regulation problem: We wish to find a positive real number k i and a full rank matrix

K i such that ∀(w, y ref ) ∈ X × R m :
1) The system (1)-( 3) is well-posed. In other words, for all ϕ e0 = (ϕ 0 , z 0 ) ∈ X e there exists a unique (weak)

solution denoted ϕ e (t) = ϕ(t) z(t) ∈ C 0 (R + , X e ) defined
∀t 0 and initial condition ϕ e (0) = ϕ e0 . 2) There exists an equilibrium point denoted

ϕ e∞ = ϕ ∞ z ∞
∈ X e , depending on w and y ref , which is exponentially stable for the system (1)-( 3). In other words, there exist positive real numbers ν e and k e such that for all t 0 ϕ e (t) -ϕ e∞ Xe k e exp -νet ϕ e0 -ϕ e∞ Xe .

3) The output y is regulated toward the reference y ref .

More precisely,

∀ϕ e0 , lim t→+∞ |Cϕ(t) -y ref | = 0. (5)
We know with the work of S. Pohjolainen in [START_REF] Pohjolainen | Robust controller for systems with exponentially stable strongly continuous semigroups[END_REF] that when A generates an exponentially stable analytic semi-group, when B is bounded and when C is A-bounded, with a rank condition, the regulation may be achieved. This result has been extended to more general exponentially stable semi-groups in [START_REF] Xu | A robust pi-controller for infinitedimensional systems[END_REF].

Theorem 1 ( [START_REF] Xu | A robust pi-controller for infinitedimensional systems[END_REF]): Assume that X is separable and that A satisfies Assumption 1. Assume moreover that :

1) the operator B is bounded; 2) the operator C is A-admissible (see [START_REF] Tucsnak | Observation and control for operator semigroups[END_REF]), i.e.

• it is A-bounded :

|Cϕ| c( ϕ X + Aϕ X ) , ∀ ϕ ∈ D(A), (6)
for some positive real number c; • there exist T > 0 and c T > 0 such that

T 0 |Ce At ϕ| 2 dt c 2 T ϕ 2 X , ∀ ϕ ∈ D(A);
3) the rank condition holds. In other words operators A, B and C satisfy

rank{CA -1 B} = m; (7)
then there exists a positive real number k * i and a m×m matrix K i , such that for all 0 < k i < k * i the operator A e given in (4) is the generator of an exponentially stable C 0 -semigroup in the extended state space X e . More precisely, the system (1) in closed loop with (3) is well-posed and the equilibrium is exponentially stable. Moreover, for all w and y ref , equation ( 5) holds (i.e the regulation is achieved).

On another hand, if one wants to address nonlinear abstract Cauchy problems or unbounded operators, we may need to follow a Lyapunov approach. For instance in the context of boundary control, a Lyapunov functional approach has allowed to tackle feedback stabilization of a large class of PDEs (see for instance [START_REF] Bastin | Stability and boundary stabilization of 1-d hyperbolic systems[END_REF] or [START_REF] Coron | Control and nonlinearity[END_REF]).

It is well known (see for instance [START_REF] Jacob | Linear port-Hamiltonian systems on infinite-dimensional spaces[END_REF]Theorem 8.1.3]) that exponential stability of the operator A is equivalent to existence of a bounded positive and self adjoint operator P in L(X) such that Aϕ, Pϕ X + Pϕ, Aϕ X -µ ϕ X , ∀ ϕ ∈ D(A), [START_REF] V Dos | Boundary control with integral action for hyperbolic systems of conservation laws: Stability and experiments[END_REF] where µ is a positive real number. We assume that this Lyapunov operator P is given. The first question, we intend to solve is the following: Knowing the Lyapunov operator P, is it possible to construct a Lyapunov operator P e associated to the extended operator A e ?

To answer this question, we first give a construction based on a well-known technique in the nonlinear finite dimensional control community named the forwarding (see for instance [START_REF] Mazenc | Adding integrations, saturated controls, and stabilization for feedforward systems[END_REF], [START_REF] Sepulchre | Integrator forwarding: a new recursive nonlinear robust design[END_REF] or more recently [START_REF] Benachour | Forwarding design with prescribed local behavior[END_REF], or [START_REF] Astolfi | Integral action in output feedback for multi-input multi-output nonlinear systems[END_REF]).

B. A Lyapunov approach for regulation

Inspired by the forwarding techniques, the following result can be obtained.

Theorem 2 (Forwarding Lyapunov functional): Assume that all assumptions of Theorem 1 are satisfied and let P in L(X) be a positive self adjoint operator such that (8) holds. Then there exist a bounded operator M : X → R m and positive real numbers p and k * i , such that for all 0 < k i < k * i , there exists µ e > 0 such that the operator

P e = P + pM * M -pM * -pM p I d (9)
is positive and satisfies ∀ ϕ e = (ϕ, z) ∈ D(A) × R m A e ϕ e , P e ϕ e Xe + P e ϕ e , A e ϕ e Xe -µ e ( ϕ 2 X + |z| 2 ).

Proof: The operator A satisfying Assumption 1, 0 is in its resolvent set and consequently A -1 : X → D(A) is well defined and bounded. Let M : X → R m be defined by M = CA -1 which is well defined due to the fact that

D(A) ⊆ D(C) since C is A-bounded. Moreover, with (6) ∀ϕ ∈ X |Mϕ| = |CA -1 ϕ| c A -1 ϕ X + ϕ X c ϕ X ,
where c is a positive real number. Hence, M is a bounded linear operator. Moreover, M satisfies the following equation

MAϕ = Cϕ , ∀ϕ ∈ D(A) . (11)
Let K i = (CA -1 B)

+ 2p(z -Mϕ) (Cϕ -MAϕ) + k i ϕ, PBK i z X + k i PBK i z, ϕ X -2p(z -Mϕ) MBK i k i z.
Employing equation [START_REF] Paunonen | Internal model theory for distributed parameter systems[END_REF] and MBK i = I dm , the former inequality becomes A e ϕ e , P e ϕ e Xe + P e ϕ e , A e ϕ e Xe =

Aϕ, Pϕ X + Pϕ, Aϕ X + k i ϕ, PBK i z X + k i PBK i z, ϕ X -2p(z -Mϕ) k i z. ( 13
)
Let PBK i 2 X = α which is well defined due to the boundedness assumption on B. Given a, b positive constants, the following inequalities hold

ϕ, PBK i z X 1 2a ϕ 2 X + aα 2 |z| 2 , (14)
z Mϕ

1 2b ϕ 2 + b M 2 2 |z| 2 , (15)
it yields given (8) that A e ϕ e , P e ϕ e Xe + P e ϕ e , A e ϕ e Xe -µ

+ k i a + pk i b ϕ 2 X + k i p(-2 + b M 2 ) + aα |z| 2 . ( 16
)
We pick b sufficiently small such that

-2 + b M 2 < 0. (17)
In a second step, we select a sufficiently small and p sufficiently large such that

p(-2 + b M 2 ) + aα < 0. (18)
Finally, picking k * i sufficiently small such that

-µ + k * i a + pk * i b < 0 ( 19
)
the result is obtained with

µ e = min µ -ki a -pki b , p(2 -b M 2 ) -aα . 2

C. Discussion on the result

A direct interest of the Lyapunov approach given in Theorem 2, is that it allows to give an explicit value for k * i which appears in Theorem 1. We may compute the largest value of k * i following this route. First of all, from ( 19)

k * i = sup
a,b,p,such that (17)-( 18)

µp p a + p 2 b , (20)
= µ sup a,b,p,such that (17)-( 18)

ab pa + b (21)
On another hand, taking the the value of a and b given by ( 17) and ( 18), one can rewritten them with 0 < β < 1 and

0 < θ < 1 as b = 2 M 2 β, a = 2(1 -β)θ p α . (22)
Then

ab pa + b = 2(1 -β)βθp (1 -β)θ M 2 p 2 + αβ (23)
which right expression is a function of p taking its maximum value when

p = αβ (1 -β)θ M 2 then k * i = sup
a,b,p,such that (17)-( 18)

µp p a + p 2 b , (24)
= µ sup 0<θ<1, 0<β<1 β(1 -β)θ √ α M (25)
It is reached for β = 1 2 and θ = 1 and this yields

k * i = µ 2 M √ α = µ 2 CA -1 PB(CA -1 B) -1 X .
Of course, this optimal value depends on the considered Lyapunov operator P solution of [START_REF] V Dos | Boundary control with integral action for hyperbolic systems of conservation laws: Stability and experiments[END_REF]. Note that a possible solution to this equation with µ = 1 is given for all (ϕ 1 , ϕ 2 ) in X 2 by (see [START_REF] Jacob | Linear port-Hamiltonian systems on infinite-dimensional spaces[END_REF])

ϕ 1 , Pϕ 2 X = lim t→+∞ t 0 e As ϕ 1 , e As ϕ 2 X ds.
Due to (2) it is well defined and positive. Note also that we have

P X k 2 2ν .
This implies, the following corollary. Corollary 1 (Explicit integral gain): Given a system Σ(A, B, C) satisfying the assumptions of the Theorem 1, points 1), 2), and 3) of Theorem 1 hold with

K i = (CA -1 B) -1 and k * i = ν CA -1 k 2 B(CA -1 B) -1 . ( 26
)
An interesting question would now to know in which aspect this value may be optimal.

D. Illustration on a parabolic systems

Consider the problem of heating a bar of length L = 10 with both endpoints at temperature zero. We control the heat flow in and out around the points s = 2, 5, and 7 and measure the temperature at points 3, 6, and 8. The problem is to find an integral controller such that the measurements at s = 3, 6, and 8 are regulated to (for instance) 1, 3, and 2, respectively. Thus the control system is governed by the following PDE

φ t (s, t) = φ ss (s, t) + 1 [ 3 2 , 5 2 ] (s)u 1 (t) + 1 [ 9 2 , 11 2 ] (s)u 2 (t) + 1 [ 13 2 , 15 2 ] (s)u 3 (t), (s, t) ∈ (0, 10) × (0, ∞) (27) where φ : [0, +∞) × [0, 10] → R with boundary conditions φ(0, t) = φ(10, t) = 0 φ(s, 0) = φ 0 (s), (28)
where

1 [a,b] : [0, 10] → R denotes the characteristic function on the interval [a, b], i.e., 1 [a,b] (s) = 1 ∀ s ∈ [a, b], 0 ∀ s ∈ [a, b].
The output and the reference are given as

y(t) =   φ(t, 3) φ(t, 6) φ(t, 8)   , y ref =   1 3 2   .
Let the state space be the Hilbert space X = L 2 ((0, 10), R) with usual inner product, and let the input space and the output space be equal to R 3 . Clearly, from (28), we get the semigroup generator A : D(A) → X , the input operator B : R 3 → X and the output operator C : D(A) → R 3 as follows:

D(A) = {ϕ ∈ H 2 (0, 10) | ϕ(0) = ϕ(10) = 0}, and Aϕ = ϕ ss ∀ ϕ ∈ D(A), Bu = 1 [ 3 2 , 5 2 ] u 1 + 1 [ 9 2 , 11 2 ] u 2 + 1 [ 13 2 , 15 2 ] u 3 , and
Cϕ =   ϕ(3) ϕ(6) ϕ(8)   .
Moreover, note that with Sobolev embedding, an integration by part and by completing the square, we have for all ϕ in D(A)

sup s∈(0,10) |ϕ(s)| ≤ c 10 0 ϕ(s) 2 ds + c 10 0 ϕ s (s) 2 ds ≤ c ϕ X + c 10 0 |ϕ(s)ϕ ss (s)|ds ≤ 3 2 c ϕ X + 1 2 c ϕ ss X .
Hence C is A-bounded. Moreover, by direct computation we find that which gives CA -1 ≤ 6.2466. We have For the open-loop system, consider the Lyapunov operator P = I d . Then the growth rate may be taken as µ = π 2 50 . It is easy to see that K i = 4.2433, and B ≤ √ 3. Putting together the numerical values into the formula (26) allows to estimate the tuning parameter

CA -1 B = -1 10   14
K i =   -1.
k * i = ω 2 BK i CA -1 ≈ 2.1498 * 10 -3 .
With Corollary 1, the integral controller (3) with 0 < k i < 2.1498 * 10 -3 stabilizes exponentially the equilibrium along solutions of the closed-loop system and drives asymptotically the measured temperatures to the reference values for any initial condition.

III. CASE OF BOUNDARY REGULATION FOR HYPERBOLLIC PDES

In the following section we adapt this framework to hyperbolic PDE systems with boundary control.

A. System description

To illustrate the former abstract theory, we consider the case of hyperbolic partial differential equations as studied in [START_REF] Coron | Dissipative boundary conditions for one-dimensional nonlinear hyperbolic systems[END_REF]. More precisely, the system is given by a one dimensional n×n hyperbolic system The output to be regulated to a prescribed value denoted by y ref , is given as a disturbed linear combination of the boundary conditions. Namely, the outputs to regulate are in R m given as

φ t (s, t) + Λ 0 (s)φ s (s, t) + Λ 1 (s)φ(s, t) = 0 s ∈ (0, 1), t ∈ [0, +∞), (29) where φ : [0, +∞) × [0, 1] → R n Λ 0 (s) = diag{λ 1 (s), . . . , λ n (s)} λ i (s) > 0 ∀i ∈ {1, . . . , } λ i (s) < 0 ∀i ∈ { + 1, . . . , n}, where the maps Λ 0 is in C 1 ([0, 1]; D n ) and Λ 1 is in C 1 ([0, 1]; R n×n ) with the initial condition φ(0, s) = φ 0 (s) for s in [0, 1] where φ 0 : [0, 1] → R n and with the boundary conditions φ + (t, 0) φ -(t, 1) = K φ + (t, 1) φ -(t, 0) + Bu(t) + w b (30) = K 11 K 12 K 21 K 22 φ + (t, 1) φ -(t, 0) + B 1 B 2 u(t) + w b (31) where φ = φ + φ - with φ + in R , φ -in R n-
y(t) = L 1 φ + (t, 0) φ -(t, 1) + L 2 φ + (t, 1) φ -(t, 0) + w y , (32)
where L 1 and L 2 are two matrices in R m×n and w y is an unknown disturbance in R m . We wish to find a positive real number k i and a full rank matrix K i such that

u(t) = k i K i z(t) , z t (t) = y(t) -y ref , z(0) = z 0 (33)
where z(t) takes value in R m and z 0 ∈ R m solves the regulation problem ∀y ref ∈ R m . The state space denoted by X e of the system (29)-(30) in closed loop with the control law (33) is the Hilbert space defined as:

X e = (L 2 (0, 1), R n ) × R m ,
equipped with the norm defined for ϕ e = (φ, z) in X e as:

v Xe = φ L 2 ((0,1),R m ) + |z|.
We introduce also a smoother state space defined as:

X e1 = (H 1 (0, 1), R n ) × R m .

B. Output regulation result

In this section, we give a set of sufficient conditions allowing to solve the regulation problem as described in the introduction. Our approach follows what we have done in the former section. Following [2, Proposition 5.1, p161] we consider the following assumption.

Assumption 2 (Input-to-State Exponential Stability): There exist a C 1 function P : [0, 1] → D n , a real numbers µ > 0, P , P and a positive definite matrix S in R n×n such that (P (s)Λ 0 (s)) s -P (s)Λ 1 (s) -Λ 1 (s)P (s) -µP (s), (34) P I dn P (s) P I dn , ∀s ∈ [0, 1],

and

-K + P (1)Λ 0 (1)K + + K -P (0)Λ 0 (0)K --S. (36)
where

K + = I d 0 K 21 K 22 , K -= K 11 K 12 0 I dn- (37)
As it will be seen in the following section, this assumption is a sufficient condition for exponential stability of the equilibrium of the open loop system. It can be found in [START_REF] Bastin | Stability and boundary stabilization of 1-d hyperbolic systems[END_REF] in the case in which S may be semi-definite positive. The positive definiteness of S is fundamental to get an input-tostate stability (ISS) property of the open loop system with respect to the disturbances on the boundary. More general results are given in [START_REF] Prieur | Iss-lyapunov functions for timevarying hyperbolic systems of balance laws[END_REF]. The second assumption is related to the rank condition. Let Φ : [0, 1] → R n×n be the matrix function solution to the system Φ s (s) = Λ 0 (s) -1 Λ 1 (s)Φ(s), Φ(0) = I dn .

We denote Φ(s) = Φ 11 (s) Φ 12 (s) Φ 21 (s) Φ 22 (s) and

Φ + (1) = Φ 11 (1) Φ 12 (1) 0 I dn- , Φ -(1) = I d 0 Φ 21 (1) Φ 22 (1)
Assumption 3 (Rank condition 1): The matrix in R n×n Φ -(1) -KΦ + (1) is full rank and so is the matrix T defined as

T 1 = (L 1 Φ -(1) + L 2 Φ + (1)) (Φ -(1) -KΦ + (1)) -1 B. (38)
Another rank condition has to be introduced. This one is used when solving the forwarding equation. Let Ψ : [0, 1] → R n×n be the matrix function solution to the system

Ψ s (s) = Ψ(s) (Λ 1 (s) -Λ 0s (s)) Λ 0 (s) -1 , Ψ(0) = I dn . ( 39
)
Assumption 4 (Rank condition 2): The matrix in R n×n

Ψ(1)Λ 0 (1)K + -Λ 0 (0)K - ( 40
)
is full rank and so is the matrix

T 2 = -L 1 B + M Λ 0 (0) B 1 0 -Ψ(1)Λ 0 (1) 0 B 2
where

M = (L 1 K + L 2 ) (Λ 0 (0)K --Ψ(1)Λ 0 (1)K + ) -1 . ( 41
)
With these assumptions, the following result may be obtained. Theorem 3 (Regulation for hyperbolic PDE systems): Assume that Assumptions 2, 3 and 4 are satisfied then with

K i = T -1 2 there exists k * i > 0 such that for all 0 < k i < k * i
the output regulation is obtained. More precisely, for all (w b , w y , y ref ) in R p × R m × R m , the following holds. 1) For all (φ 0 , z 0 ) in X e (resp. X 1e ) which satisfies the boundary conditions (30) (resp. the C 1 compatibility condition), there exists a unique weak solution to (29)-( 30)-(33) that we denote v and which belongs to C 0 ([0, +∞); X e ) (Respectively, strong solution in:

C 0 ([0, +∞); X e1 ) ∩ C 1 ([0, +∞); X e ) ). (42)
2) There exists an equilibrium state denoted v ∞ in X e which is globally exponentially stable in X e for system (29)-( 30)-(33). More precisely, we have for all t 0:

v(t) -v ∞ Xe k exp(-νt) v 0 -v ∞ Xe . (43)
3) Moreover, if v 0 satisfies the C 1 -compatibility condition and is in X 1e , the regulation is achieved, i.e.

lim t→+∞ |y(t) -y ref | = 0. ( 44
)
The next section is devoted to the proof of this result.

C. About this result

The first assumption needed in Theorem 3 is Assumption 2. When considering only integral control laws, there is no hope to obtain the result without assuming exponential stability of the open loop system. Assumption 2 is slightly more restrictive than exponential stability since it requires an ISS property with respect to the input u. In the case in which this assumption is not satisfied for a given hyperbolic system, a possibility is to modify the boundary condition via a static output feedback (or proportional feedback) following the route of [START_REF] Bastin | Stability and boundary stabilization of 1-d hyperbolic systems[END_REF] in order to satisfy this assumptions.

One interest of our approach is that, part of the exponential stability of the closed loop system, only Assumptions 3 and 4 which are rank conditions involving the boundary conditions have to be satisfied. In the case in which the two above mentioned assumptions are not satisfied, we may obtain these properties by adding a proportional feedback and consequently changing the value of K in T 1 and T 2 to obtain these rank conditions. These Assumptions 3 and 4 are version of Point 3) in Theorem 1.

In the particular case in which Λ 0 is constant and Λ 1 = 0, the matrix function Φ(s) and Ψ(s) are simply equal to identity for all s in [0, 1]. In that case, it yields

T 1 = (L 1 + L 2 )(I dn -K) -1 B, (45)
and,

T 2 = -L 1 B + (L 1 K + L 2 )(K --K + ) -1 B 1 -B 2 = -L 1 B + (L 1 K + L 2 ) K 11 -I d K 12 -K 21 I dn--K 22 -1 × I d 0 0 -I d B = -L 1 B + (L 1 K + L 2 ) × I d 0 0 -I d K 11 -I d K 12 -K 21 I dn--K 22 -1 B = -L 1 B -(L 1 K + L 2 )(I dn -K) -1 B = [-L 1 (I dn -K) -(L 1 K + L 2 )] (I dn -K) -1 B = -(L 1 + L 2 )(I dn -K) -1 B. (46)
Hence when Λ 0 is constant and Λ 1 = 0, Assumption 3 and Assumption 4 are equivalent. Also, an interesting aspect of this Lyapunov approach is that explicit values of the supremum value of the gain k * i may be given. For instance, as in [START_REF] Trinh | Design of integral controllers for nonlinear systems governed by scalar hyperbolic partial differential equations[END_REF] consider the very particular case of a transport equation. In this case the system is simply

φ t (s, t) + φ s (s, t) = 0, s ∈ (0, 1), t ∈ [0, +∞) φ(t, 0) = u(t) + w b y(t) = φ(1, t) + w y We can apply Theorem 3 with n = 1, Λ 0 (s) = -1, Λ 1 (s) = 0, K = 0, B = 1, L 1 = 0, L 2 = 1. This yields Ψ(s) = 1, Ψ(s) = 1, T 1 = 1, T 2 = -1.
Hence, Assumptions 3 and 4 are satisfied. Assumption 2 is satisfied for all µ > 0 with P (s) = e -µs , S = 1, P = 1, P = e -µ . In that case, employing theorem 3, it yields that there exists k * i > 0 such that for all 0 < k i < k * i with u(t) = -k i z, ż = y, the output regulation is obtained and so the output converges asymptotically to zero. Following the proof of Theorem 3, equation (85) gives

k * i = µe -µ .
This bound is better then the one obtained in [START_REF] Trinh | Design of integral controllers for nonlinear systems governed by scalar hyperbolic partial differential equations[END_REF] for the linear transport equation (its maximal value is obtained for µ = 1 and is 1 √ e . Note however that similar to the bound of [START_REF] Trinh | Design of integral controllers for nonlinear systems governed by scalar hyperbolic partial differential equations[END_REF], the result obtained with our novel Lyapunov functional is far from the value we get following a frequency approach ( π 2 in this case). Recently in [START_REF] Coron | PI controllers for 1-D nonlinear transport equation[END_REF], the Lyapunov functional obtained in [START_REF] Trinh | Design of integral controllers for nonlinear systems governed by scalar hyperbolic partial differential equations[END_REF] has been modified to reach this optimal value of the integral gain. A natural question for future research topic is to know if it is possible to modify the Lyapunov functional obtained in Theorem 3 following the methods of [START_REF] Coron | PI controllers for 1-D nonlinear transport equation[END_REF] to remove the conservatism. D. Illustration in a 2 × 2 hyperbolic system Theorem 3 generalizes many available results on output regulation via integral action for hyperbolic PDEs available in the literature. For instance, the case of 2 × 2 linear hyperbolic systems has been considered in [START_REF] Trinh | Boundary pi controllers for a star-shaped network of 2x2 systems governed by hyperbolic partial differential equations (long version)[END_REF], [START_REF] V Dos | Boundary control with integral action for hyperbolic systems of conservation laws: Stability and experiments[END_REF], (see also [START_REF] Bastin | Stability and boundary stabilization of 1-d hyperbolic systems[END_REF]Section 2.2.4]). The case of cascade of such systems is also considered in [START_REF] Trinh | Stability and output regulation for a cascaded network of 2x2 hyperbolic systems with pi control[END_REF]. Note also that in [START_REF] Terrand-Jeanne | Regulation of inhomogeneous drilling model with a p-i controller[END_REF], this procedure is applied on a Drilling model which is composed of a hyperbolic PDE coupled with a linear ordinary differential equation.

In order to compare the way we improve existing results, the same example as in [START_REF] V Dos | Boundary control with integral action for hyperbolic systems of conservation laws: Stability and experiments[END_REF] is considered. In this context, the linearized de Saint-Venant equations can be written in the form of ( 29)-(30). After normalization, one gets :

Λ 0 (s) = c 0 0 -d and Λ 1 (s) = 0 2×2 , ∀s (47)
where c > 0 and d > 0 and

K = 0 k 0 k 1 0 and B = b 0 0 0 b 1 , (48)
with b 0 = 0 and b 1 = 0. For the system (29)-(30) with these parameters, it is shown in [START_REF] V Dos | Boundary control with integral action for hyperbolic systems of conservation laws: Stability and experiments[END_REF] that the output of dimension m = 2 defined in (32) with

L 1 = c c+d 0 0 -1 c+d and L 2 = 0 d c+d 1 c+d 0 ( 49
)
can be regulated with an integral control law provided

|k 0 k 1 | < 1 , |k 0 | < 1 , |k 1 | < c d . ( 50
)
On another hand, employing ([7]- [START_REF] V Dos | Boundary control with integral action for hyperbolic systems of conservation laws: Stability and experiments[END_REF]), Assumptions 2 is satisfied assuming that |k 0 k 1 | < 1. Moreover, with equations ( 45) and ( 46), it yields,

T 1 = -T 2 = 1 c + d c d 1 -1 1 -k 0 -k 1 1 -1 b 0 0 0 b 1 .
This matrix is well defined and full rank if |k 0 k 1 | < 1 and consequently Assumptions 3 and 4 are always satisfied. Hence, employing Theorem 3, both outputs defined in (49) can be regulated with an integral control law with the only assumption that |k 0 k 1 | < 1.

Then

K i = T -1 2 = ϑ b 1 (1 -k 0 ) b 1 (d + ck 0 ) b 0 (1 -k 1 ) b 0 (c + dk 1 ) (51) with ϑ = -(c + d) 2 (1 -k 0 k 1 ) 2 b -1 0 b -1 1 [(1 -k 0 )(c + dk 1 ) + (1 -k 1 )(d + ck 0 )] (52)
and

k * i = √ µP |M |Ψ c T -1 2 (53)
µ is given in [START_REF] V Dos | Boundary control with integral action for hyperbolic systems of conservation laws: Stability and experiments[END_REF], and P the lower bound of the Lyapunov can be deduced easily from the expression of the Lyapunov function involved. M has been defined above, with T 2 . As Ψ is the identity matrix, Ψ is 1. Remark that in [START_REF] V Dos | Boundary control with integral action for hyperbolic systems of conservation laws: Stability and experiments[END_REF], K i is diagonal and here is full matrix. Note that some other choices of K i are possible as long as

T 2 K i + K i T 2 > 0.
To conclude, a work is needed to transpose this approach to the global de Saint-Venant equations this is the aim of another paper.

IV. PROOF OF THEOREM 3

The proof of Theorem 3 is divided into three steps. In a first part, it is shown that with Assumption 3, it can be shown that the closed loop system (29)-( 30)-(33) admits a steady state. In a second step, it is established that the desired regulation is obtained provided the steady state is exponentially stable. Finally, the construction of an appropriate Lyapunov functional is performed to show the exponential stability of the equilibrium.

A. Stabilization implies regulation

In this first subsection, we explicitly give the equilibrium state of the system (29)-( 30)-(33). We show also that if we assume that k i and K i are selected such that this equilibrium point is exponentially stable along the closed loop, then the regulation is achieved.

1) Definition of the equilibrium: The first step of the study is to exhibit equilibrium denoted φ ∞ , z ∞ of the disturbed hyperbolic PDE in closed loop with the boundary integral control (i.e. system (29)-( 33)).

We have the following proposition. Proposition 1: Assumption 3 is a necessary and sufficient condition for the existence of an equilibrium of the system (29)-( 30)-(33). Moreover, if Assumption 3 holds then point 1) of Theorem 3 holds.

Proof: First of all, equilibria are such that

φ ∞s (s) = -Λ 0 (s) -1 Λ 1 (s)φ ∞ (s), for all s in [0, 1]. Hence, φ ∞ (s) = Φ(s)φ ∞ (0). (54)
Hence,

φ ∞+ (0) φ ∞-(1) = Φ -(1)φ ∞ (0), φ ∞+ (1) φ ∞-(0) = Φ + (1)φ ∞ (0)
Moreover, with z t = 0, we have

L 1 φ ∞+ (0) φ ∞-(1) + L 2 φ ∞+ (1) φ ∞-(0) = y ref -w y ,
Hence,

(L 1 Φ -(1) + L 2 Φ + (1)) φ ∞ (0) = y ref -w y , (55)
On another side, boundary conditions (30) gives

(Φ -(1) -KΦ + (1)) φ ∞ (0) = Bk i K i z ∞ + w b (56)
For all w y and y ref both in R m , w b in R p , by Assumption 3 and since the matrix K i is full rank the former equation and (55) admit a unique solution (z ∞ , φ ∞ (0)) given as

z ∞ = K -1 i k i T -1 1 (y ref -w y ) - K -1 i k i T -1 1 (L 1 Φ -(1) + L 2 Φ + (1)) × (Φ -(1) -KΦ + (1)) -1 w b (57)
and,

φ ∞ (0) = (Φ -(1) -KΦ + (1)) -1 k i BK i z ∞ + (Φ -(1) -KΦ + (1)) -1 w b . (58)
Finally, in that case, we can introduce φ(s, t) = φ(s, t)φ ∞ (s) and z(t) = z(t) -z ∞ . It can be checked that φ, z satisfies the following system: φt (s, t) + Λ 0 (s) φs (s, t) + Λ 1 (s) φ(s, t) = 0 , s ∈ (0, 1),

zt = L 1 φ+ (t, 0) φ-(t, 1) + L 2 φ+ (t, 1) φ-(t, 0) (59)
with the boundary conditions

φ+ (t, 0) φ-(t, 1) = K φ+ (t, 1) φ-(t, 0) + Bu(t), (61)
u(t) = k i K i z(t). ( 62
)
As it is shown in [START_REF] Bastin | Stability and boundary stabilization of 1-d hyperbolic systems[END_REF], for each initial condition ṽ0 = ( φ0 , z0 ) in X e which satisfies the boundary conditions (30), there exists a unique weak solution that we denoted ṽ and which belongs to C 0 ([0, +∞); X e ). Moreover, if the initial condition ṽ0 satisfies also the C 1 -compatibility condition (see [START_REF] Bastin | Stability and boundary stabilization of 1-d hyperbolic systems[END_REF] for more details) and lies in X e1 then the solution lies in the set defined in (42).

2) Sufficient conditions for Regulation: In the following, we show that the regulation problem can be rephrased as a stabilization of the equilibrium state introduced previously.

Proposition 2: Assume Assumption 3 holds and that there exist a functional V e : X e → R + , and positive real numbers µ e and L e such that:

v ∞ -v 2 Xe L e V e (v) L e v ∞ -v 2 Xe . ( 63
)
Assume moreover that for all v 0 in X e and all t 0 in R + such that the solution v of system (29)-( 30)-(33) initiated from v 0 is C 1 at t = t 0 , we have:

Ve (t) -µ e V e (t), (64)
where with a slight abuse of notation V e (t) = V e (v(t)). Then points 1), 2) and 3) of Theorem 3 hold.

Proof: Point 1) is directly obtained from Proposition 1. The proof of point 2) is by now standard. Let v 0 be in X e1 and satisfies the C 0 and C 1 -compatibility conditions. It yields that v is C 1 for all t. Consequently, (64) is satisfied for all t 0. With Grönwall's lemma, this implies that: V e (v(t)) e -µet V e (v 0 ). Hence with (63), this implies that (43) holds with k = L e and ν = µe 2 for initial conditions in X e1 . X e1 being dense in X e , the result holds also with initial condition in X e and point 2) is satisfied.

On another hand, we have

y(t) -y ref = L 1 φ + (t, 0) φ -(t, 1) + L 2 φ + (t, 1) φ -(t, 0) + w y -y ref , (65)
= L 1 φ+ (t, 0) φ-(t, 1) + L 2 φ+ (t, 1) φ-(t, 0) , (66)
with φ(t, x) = φ(t, x) -φ ∞ . To show that equation ( 44) holds, we need to show that the right hand side of the former equation tends to zero. This may be obtained provided the initial condition is in X 1 . Indeed, let v 0 be in X 1 and satisfies C 1 -compatibility conditions. With (42), we know that v t ∈ C([0, ∞); X e ). Moreover, v t satisfies the dynamics system (29)-( 30)-(33) with w b = 0, w y = 0, y ref =0 (simply differentiate with time these equations). Hence, v t (t) Xe converges exponentially toward 0 and in particular φt (t, •) (L 2 (0,1),R n ) ke -νt v t (0) .

On another hand, employing (29), it yields:

φs (t, •) L 2 ((0,1),R n ) = Λ -1 0 ( φt (t, •) + Λ 1 (•)φ(t, •)) L 2 ((0,1),R n ) . Hence, φs (t, •) L 2 ((0,1),R n ) c φt (t, •) L 2 ((0,1),R n ) + φ(t, •)) L 2 ((0,1),R n ) . ( 67
)
where c is a positive constant. Consequently φs (t, •) L 2 ((0,1),R n ) converges also to zero and so is φ(t, •) H 1 ((0,1),R n ) . With Sobolev embedding

sup x∈[0,1] | φ(t, x)| C φ(t, •) H 1 ((0,1),R n ) ,
where C is a positive real number. It implies that:

lim t→+∞ | φ(t, 1)| + | φ(t, 0)| = 0.
Consequently, with (65), it yields that (44) holds and point 3) is satisfied.

2

With this proposition in hand, to prove the Theorem 3, it is sufficient to construct a Lyapunov functional V e which satisfies (63)-( 64) along C 1 -solutions of ( 29)-( 30)-(33) or equivalently along C 1 -solutions of (59)-(61). This is considered in the next section following the route of Section II-B.

B. Lyapunov functional construction 1) Open loop ISS: Inspired by the Lyapunov functional construction introduced in [START_REF] Coron | Dissipative boundary conditions for one-dimensional nonlinear hyperbolic systems[END_REF] (see also [START_REF] Bastin | Stability and boundary stabilization of 1-d hyperbolic systems[END_REF]), we know that typical Lyapunov functionals allowing to exhibit stability property for this type of hyperbolic PDE are given as functional V : L 2 ((0, 1), R n ) → R + defined as

V (ϕ) = 1 0 ϕ(s) P (s)ϕ(s)ds , (68)
where P : [0, 1] → D n is a C 1 function. Typically in [START_REF] Coron | Dissipative boundary conditions for one-dimensional nonlinear hyperbolic systems[END_REF], these functions are taken as exponential.

With a slight abuse of notation, we write V (t) = V ( φ(•, t)) and we denote by V (t) the time derivative of the Lyapunov functional along solutions which are C 1 in time. In our context, with Assumption 2, it yields the following proposition.

Proposition 3: If Assumption 2 holds, there exists a positive real number c such that for every solution φ of (29)-(30) initiated from ( φ0 , z0 ) in X e which satisfies (61)

V (t) -µV (t) + c|u(t)| 2 . (69)
Proof: First of all, with (29), With an integration by part, this implies

V (t) = - 1 0 2φ(t, s) P (s)Λ 0 (s)φ s (t, s)ds
V (t) = 1 0 φ(t, s) [(P (s)Λ 0 (s)) s ] φ(t, s)ds - 1 0 φ(t, s) (P (s)Λ 1 (s) + Λ 1 (s) P (s))φ(t, s)ds -φ(t, 1) P (1)Λ 0 (1)φ(t, 1) + φ(t, 0) P (0)Λ 0 (0)φ(t, 0).
With (34), it gives

V (t) -µV (t) -φ(t, 1) P (1)Λ 0 (1)φ(t, 1)
+ φ(t, 0) P (0)Λ 0 (0)φ(t, 0).

With the boundary condition (31) and (36), this implies

V (t) -µV (t) -φ + (1) φ -(0) S φ + (1) φ -(0) + 2 φ + (1) φ -(0) Qu(t) + u(t) Ru(t), (70)
where,

R = -0 B 2 (P (1)Λ 0 (1) + Λ 0 (1)P (1)) 0 B 2 + B 1 0 (P (0)Λ 0 (0) + Λ 0 (0)P (0)) 0 B 1 , (71)
and,

Q = -K + (P (1)Λ 0 (1) + Λ 0 (1)P (1)) 0 B 2 + K -(P (0)Λ 0 (0) + Λ 0 (0)P (0)) B 1 0 .
Since S is positive definite, selecting c sufficiently large, it yields

-S Q Q R -c I dm 0 .
Consequently, (70) implies that (69) holds.

2) Forwarding approach to deal with the integral part: Following the route of Section II-B, a Lyapunov functional is designed from V adding some terms to take into account the state of the integral controller. Let the operator M : L 1 ((0, 1); R n ) → R m be given as

Mϕ = 1 0 M Ψ(s)ϕ(s)ds ( 72
)
where Ψ is the matrix function defined in (39), and M is a matrix in R m×n defined in (41). Following the Lyapunov functional construction in Theorem 2, we consider the candidate Lyapunov functional V e : L 2 ((0, 1); R n ) × R m given as

V e (ϕ, z) = V (ϕ) + p(z -Mϕ) (z -Mϕ). (73)
In the following theorem, it is shown that by selecting properly K i , k i and p, this function is indeed a Lyapunov functional for the closed loop system. Again, with a slight abuse of notation, we write V e (t) = V e ( φ(•, t), z(t)) and we denote by Ve (t) the time derivative of the Lyapunov functional along solutions which are C 1 in time. Proposition 4: Assume that Assumptions 2 and 3 hold. Then there exists a matrix K i in R m×m and k * i > 0 such that for all 0 < k i < k * i , there exist positive real numbers L e and µ e such that for all (ϕ, z) in X

e 1 L e ϕ 2 X + |z| 2 V e (ϕ, z) L e ϕ 2 X + |z| 2 , ( 74
)
and along C 1 solution of the system (59)-( 61)-( 62)

Ve (t) -µ e V e (t) , ∀t ∈ R + . (75)
Proof: With (35), it yields for all ϕ in L 2 ((0, 1); R n ), P ϕ 2

pk i < µP |M | 2 Ψ 2 , k 2 i < pk i c T -1 2 ,
this yields the existence of µ e such that equation (75) holds. This is obtained for all k i < k * i when

k * i = √ µP |M |Ψ c T -1 2 , (85)
and p < µP

k i |M | 2 Ψ 2 .
2 With this proposition, the proof of Theorem 3 is completed.

V. CONCLUSION

In the last three decades, the regulation problem has been studied for different classes of distributed parameter systems. Most of existing results follow a semigroup approach and the perturbation theory for linear operator. In this paper we have shown that is was also possible to construct Lyapunov functionals to address the regulation problem in the case in which is used an integral action. This framework allows to explicitly give an integral gain. Moreover, it is no more necessary to impose boundedness of control or measurement operators to guarantee the regulation. This is applied to PDE hyperbolic systems and this allows to generalize many available results in this field.

A

  Terrand-Jeanne, Vincent Andrieu, Valérie dos Santos Martins, C.-Z Xu To cite this version: A Terrand-Jeanne, Vincent Andrieu, Valérie dos Santos Martins, C.-Z Xu. Adding integral action for open-loop exponentially stable semigroups and application to boundary control of PDE systems. IEEE Transactions on Automatic Control, 2020, 10.1109/TAC.2019.2957349. hal-01971584 I. INTRODUCTION

0

  It is easy to see that the above matrix is regular. Consequently all Assumptions of Theorem 1 hold. With Corollary 1, it is possible to compute explicitly the integral controller gain. By direct computation we have for all ϕ in X CA -1 ϕ = (s -10)ϕ(s)ds + 3 0 (3 -s)ϕ(s)ds

3 5 10 0 6 0 5 100 8 0

 10658 (s -10)ϕ(s)ds + (6 -s)ϕ(s)ds4 (s -10)ϕ(s)ds + (8 -s)ϕ(s)ds     ,

  and where w b in R p is an unknown disturbance, u(t) is a control input taking values in R m and K, B are matrices of appropriate dimensions.

- 1 0

 1 φ(t, s) P (s)Λ 1 (s) + Λ 1 (s) P (s) φ(t, s)ds

L 2 ( 1 0M 1 0M 1 0M 1 0M

 21111 (0,1);R n ) V (ϕ) P ϕ 2 L 2 ((0,1);R n ) .(76)Let Ψ > 0 be such that|Ψ(s)| Ψ , ∀s ∈ [0, 1].Note that for all ϕ in L 2 ((0, 1); R n , by Cauchy-Schwartz inequality,|Mϕ| |M | Ψ ϕ L 2 ((0,1);R n ) .(77)Hence, for each p > 0 equation (74) holds.Note that along C 1 solution to system (59), we haveM φt (t, •) = M(-Λ 0 (•) φs (t, •) -Λ 1 (•) φ(t, •)) = Ψ(s)(-Λ 0 (s) φs (t, s) -Λ 1 (s) φ(t, s))dsWith an integration by part this impliesM φt (t, •) = (Ψ(s)Λ 0 (s)) s φ(t, s)ds -Ψ(s)Λ 1 (s) φ(t, s))ds -M Ψ(1)Λ 0 (1) φ(t, 1) -Λ 0 (0) φ(t, 0) .This gives,M φt (t, •) = (Ψ s (s)Λ 0 (s) + Ψ(s)(Λ 0s (s) -Λ 1 (s))) φ(t, s)ds -M Ψ(1)Λ 0 (1) φ(t, 1) -Λ 0 (0) φ(t, 0) .With the definition of Ψ, it yieldsM φt (t, •) = -M Ψ(1)Λ 0 (1) φ(t, 1) -Λ 0 (0) φ(t, 0) .With the boundary condition (61), it yieldsM φt (t, •) = -M (Ψ(1)Λ 0 (1)K + -Λ 0 (0)K -) φ+ (t, the definition of M , it implies M φt (t, •) = (L 1 K +

  A e ϕ e , P e ϕ e Xe + P e ϕ e , A e ϕ e Xe = Aϕ, Pϕ X + Pϕ, Aϕ X

	Moreover, we have
	-1 which exists due to the third assump-
	tion of Theorem 1. Note that,
	ϕ e , P e ϕ e Xe = ϕ, Pϕ X + p(z -Mϕ) (z -Mϕ), (12)
	hence P e is positive. This candidate Lyapunov functional is
	similar to the one given in [4, Equation (34)]. It is selected
	following a forwarding approach.

  Bu(t).Let now, K i = T -1 2 . Hence, this gives with u = k i K i z,Ve (t) -µV (t) + ck 2 i |K i z(t)| 2 -2p|z(t)| 2 k i + 2pk i (Mφ(•, t)) z(t), (80)With (77), and completing the square it yields for all ϕ in L 2 ((0, 1); R n ) and z in R m ,2(Mϕ) z |Mϕ| 2 + |z| 2 , Ve (t) -µV (t) + pk i |M | 2 Ψ 2 φ(•, t) 2 L 2 ((0,1);R n ) + ck 2 i |K i | 2 -pk i |z(t)| 2 . (83) |K i | 2 -pk i |z(t)| 2 . (84)

						(81)
			|M | 2 Ψ	2 ϕ 2 L 2 ((0,1);R n ) + |z| 2 .	(82)
	Merging the last two inequality yields,
	With (76), this yields
	Ve (t)	-µ + pk i	|M | 2 Ψ P	2	V (t)
	+ ck 2 i Note that if
						L 2 ) φ+ (t, 1) φ-(t, 0)
						+ M Λ 0 (0)	B 1 0	-Ψ(1)Λ 0 (1)	0 B 2	u(t)
						On another hand,
						z t (t) = (L 1 K + L 2 ) φ+ (t, 1) φ-(t, 0) + L 1 This gives,
						M φt (t, •) = z t (t) -L 1 Bu(t)
						+ M Λ 0 (0)	B 1 0	-Ψ(1)Λ 0 (1)	0 B 2	u(t).
						Hence, it yields

M φt (t, •) = z t (t) + T 2 u(t).

(78)

We recognize here equation (

11

) when u = 0. This gives with (

69

)

Ve (t) -µV (t) + c|u(t)| 2 -2p(z(t) -Mφ(•, t)) T 2 u(t). (79)

Alexandre Terrand-Jeanne graduated in electrical engineering from ENS Cachan, France, in 2013. After one year in the robotic laboratory "Centro E.Piaggio" in Pisa, Italy, he is currently a doctoral student in LAGEP, university of Lyon 1. His PhD topic concerns the stability analysis and control laws design for systems involving hyperbolic partial differential equations coupled with nonlinear ordinary differential equations. This work is under the supervision of V. Dos Santos Martins, V. Andrieu and M. Tayakout-Fayolle.


Martine Ben Amar
email: benamar@lps.ens.fr

Nonlinear visco-elasticity of soft tissues under cyclic deformations

Keywords: Finite elasticity and visco-elasticity; cyclic deformations, hysteresis, living fibrotic tissues

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

Nonlinear visco-elasticity of soft tissues under cyclic deformations

Introduction

As the description of fibrous tissues evolves and makes progress [START_REF] Gasser | Hyperelastic modelling of arterial layers with distributed collagen fibre orientation[END_REF][START_REF] Ogden | Biomechanical Modelling at the molecular, cellular and tissue levels CISM[END_REF], experimental set-ups have also been improved allowing the automatic determination of stress-strain results which can be compared to modelling of complex hyper-elastic tissues. Some theoretical models originate from a microscopic analysis of the structural properties of the sample under tension and the changes which may occur when a tensional state is applied. The case of elongation of tendons is an example [START_REF] Ciarletta | Amar A finite dissipative theory of temporary interfibrillar bridges in the extra-cellular matrix of ligaments and tendons[END_REF][START_REF] Shearer | Sequential straightening 280 and loading viscoelasticity: A model of the time-dependent behaviour of tendon fascicles incorporating fibril recruitment[END_REF]; the micro-structures have been analyzed in details, the link existing between the main components: collagen fibers, extra-cellular matrix and proteo-glycans, has been incorporated in the formalism with the aim to deduce a stress-strain compared to experimental data [START_REF] Pioletti | On the independence of time and strain effects in the stress relaxation of ligaments and tendons[END_REF][START_REF] Natali | Anisotropic elasto-damage constitutive model for the biomechanical analysis of tendons[END_REF]. Elaborated models can be established where the involved parameters rely on microscopic quantities as done for composite materials [START_REF] Muliana | A concurrent micromechanical model for predicting nonlinear viscoelastic responses of composites reinforced with solid spherical particles[END_REF] . Hysteresis and plasticity behaviors can also be recovered. Such studies require a profound knowledge of the biological structures and of the evolving interactions between them under stress conditions. However, they remain specific to the system under study and need to be validated by stress-strain experimental curves, when they are available. Using different approaches, biophysics groups have elaborated a characterization of the material properties going from individual cells [START_REF] Legoff | A Global Pattern of Mechanical Stress Polarizes Cell Divisions and Cell Shape in the Growing Drosophila Wing Disc[END_REF] to tissues. After investigated in detail isolated cell behavior during mitosis for example or migration of connected cells in an epithelium on a solid substrate [START_REF] La Porta | Osmotic stress affects functional properties of human melanoma cell lines[END_REF], they show that 1 collections of cells, interacting with each other and/or with the extra-cellular matrix, may exhibit emergent properties, which do not necessarily mirror the properties of their constitutive el- ements. An integrated approach is then necessary at the scale of the tissue leading to continuous visco-elastic modelling where fibers are treated using the formalism of nematic liquid crystal [START_REF] Kruse | Asters, Vortices, and Rotating Spirals in Active Gels of Polar Filaments[END_REF]11]. Nonetheless, not all biological tissues exhibit a well defined organization. Let us mention the brain and fat tissues, [START_REF] Mihai | A comparison of hyperelastic constitutive models applicable to brain and fat tissues[END_REF][START_REF] Mihai | A family of hyperelastic models for human brain tissue[END_REF][START_REF] Jiang | Measuring the linear and nonlinear elastic properties of brain tissue with shear waves and inverse analysis[END_REF], or fibrotic tissues [START_REF] Amar | Towards a unified approach in the modeling of fibrosis: A review with research perspectives[END_REF] which play an important role in organogenesis but also in the development of pathologies: wound-healing [START_REF] Wu | Amar Growth and remodelling for profound circular wounds in skin[END_REF][START_REF] Bowden | A morphoelastic model for dermal wound closure[END_REF], capsules around implants [START_REF] Ben Amar | Morpho-elasticity of inflammatory fibrosis: the case of capsular contracture[END_REF], sarcoma (the generic word for cancer of tissues) and various tumors of the brain. Even if these tissues present an order at the 40 microscopic scale, like bundles of fibers with an averaged orientation [START_REF] Gasser | Hyperelastic modelling of arterial layers with distributed collagen fibre orientation[END_REF], phenomenological constitutive laws with a small number of parameters can be derived. Often, they intervene in the complex geometry of organs and a simplified constitutive law takes all its value to understand biological processes such as growth [START_REF] Amar | Anisotropic growth shapes intestinal tissues during embryo-genesis[END_REF][START_REF] Balbi | Morpho-elasticity of intestinal vill[END_REF] in embryo-genesis or in pathologies. Parameters of constitutive laws may be experimentally deduced from compression/tension-shear tests under a sequence of multiple cycles [START_REF] Budday | Mechanical characterization of human brain tissue[END_REF]. These experiments exhibit hysteretic cycles which suggest that visco-elasticity plays an important role. 50 However, if macroscopic constitutive laws exist for tissues like the Mooney-Rivlin [START_REF] Mooney | A theory of large elastic deformation[END_REF], Ogden [START_REF] Ogden | Large deformation isotropic elasticity on the correlation of theory and experiment for incompressible rubberlike solids[END_REF][START_REF] Ogden | Non-linear elastic deformations Dover Publications and Ellis Horwood[END_REF], Gent [START_REF] Gent | A new constitutive relation for rubber[END_REF], or the GOH [START_REF] Gasser | Hyperelastic modelling of arterial layers with distributed collagen fibre orientation[END_REF] models, macroscopic viscoelasticity [START_REF] Winemann | Nonlinear viscoelastic solids-a review Math[END_REF] do not offer the same diversity. Indeed, a brief history of nonlinear dissipation modeling can be found in [START_REF] Destrade | Proper formulation of viscous dissipation for nonlinear waves in solids The[END_REF], going back to the beginning of the last century and a presentation of the so-called Landau model, commonly used in acoustics [START_REF] Ye | Simulation of non-linear transient elastography: finite element model for the propagation of shear waves in homogeneous soft tissues[END_REF]. In the spirit of Landau, the theory of visco-elasticity rests on the definition of a dissipation energy rate, formulated in the same way as the elastic energy density, then introducing two new "Lamé" parameters for the viscoelastic part. But, as shown by Destrade and coworkers, to deduce correctly the viscous nominal stress or the viscous Cauchy stress is not so trivial. Fundamental properties such as objectivity (frame-indifference) or symmetry of the Cauchy stress are not always verified in the acoustic literature where the Landau's model is commonly accepted. Applying these symmetry arguments, they derive the correct viscous tensor for the Landau's model. Having in mind the fact that hyper-elasticity and visco-elasticity must represent the same intrinsic tissue properties, being two sides of the same coin, it is proposed here a treatment to systematically derive the viscous tensor for an arbitrary representation of the hyper-elastic tissue properties. A variational method based on the extremum of the dissipation rate also called the Rayleighian [START_REF] Landau | of A Course of Theoretical Physics[END_REF][START_REF] Onsager | s variational principle in soft matter[END_REF] in polymer physics is used. For each biomechanical constitutive law we show that we can construct easily its visco-elastic counter-part, and for each nominal or Cauchy stress, also its visco-elastic counterpart. As an example, a test on classical models is done for fibrotic mammary tissues at very low stretch, less than 1.2.

A cascade of hysteric loops are obtained when cycles of extensions/compression are applied. Strong nonlinearities appear that are analyzed with this formalism.The objective is to recover the stress-strain curves for the uniaxial stretch of a cuboid for low values of the stretch: λ ∼ 1 as shown in Fig. [START_REF] Gasser | Hyperelastic modelling of arterial layers with distributed collagen fibre orientation[END_REF].

Constitutive laws for hyper-elastic tissues

This section is just a reminder of 4 different models of nonlinear elasticity, which will serve as a basis for the derivation of a nonlinear dissipation formalism. We focus on the Mooney-Rivlin [START_REF] Mooney | A theory of large elastic deformation[END_REF], Ogden [START_REF] Ogden | Large deformation isotropic elasticity on the correlation of theory and experiment for incompressible rubberlike solids[END_REF][START_REF] Ogden | Non-linear elastic deformations Dover Publications and Ellis Horwood[END_REF], Landau [START_REF] Landau | A Course of Theoretical Physics[END_REF] and GOH [START_REF] Gasser | Hyperelastic modelling of arterial layers with distributed collagen fibre orientation[END_REF][START_REF] Ogden | Biomechanical Modelling at the molecular, cellular and tissue levels CISM[END_REF] models of hyper-elastic tissues. Their standard representation varies, sometimes introducing principal stretches or principal invariants of the Cauchy-Green tensors. .

The Mooney-Rivlin model

We begin with the Mooney-Rivlin model, which is restricted to two invariants I 1 and I 2 , for incompressible materials since

I 3 =1. W MR = µ 2 {(I 1 -1) + ρ(I 2 -1)} (1)
with

I 1 =Tr[C] and I 2 = 1/2(Tr[C] 2 -Tr[C 2 ]
), C representing the right Cauchy-Green strain tensor: C = F T F if F is the deformation gradient. Tr[A] means the trace of the tensor A. The shear modulus µ is positive, ρ can be negative, but only to some extent. Keeping in mind the standard tensorial relationships (see [START_REF] Holzapfel | Nonlinear Solid Mechanics a continuum approach for engineering[END_REF],p 41-42).:

∂Tr[C] ∂C = I; ∂Tr[C α ] ∂C = αC T(α-1) = αC α-1 and ∂Det[C] ∂C = C -1
(2) The Cauchy stress is then:

σ e = 2C ∂W MR ∂C -p = µ (1 + 3ρI 1 )C -ρC 2 -pI (3)
where p is the Lagrange parameter which ensures the incompressibility property.

The Ogden model

The Ogden model, introduced for rubber elasticity, turns out to be more and more employed in recent works on living tissues such as [START_REF] Mihai | A comparison of hyperelastic constitutive models applicable to brain and fat tissues[END_REF][START_REF] Ben Amar | Morpho-elasticity of inflammatory fibrosis: the case of capsular contracture[END_REF][START_REF] Budday | Mechanical characterization of human brain tissue[END_REF]. Its classical formulation is given in term 100 of the principal stretches λ i for each axis of coordinate X i but different variants can be found in the literature:

W Og = µ p 2α p λ 2α p 1 + λ 2α p 2 + λ 2α p 3 -3 + µ p 2α p ρ λ -2α p 1 + λ -2α p 2 + λ -2α p 3 -3 (4)
which can replaced by:

W Og = µ p 2α p Tr[C α p ] -3 + ρ 2 Tr[C α p ] 2 -Tr[C 2α p ] -6 (5)
when α p is integer. Only the two-term Ogden model is presented but in principle a series expansion is possible varying α p and µ p . According to Eq.( 3), we can deduce the nominal stress:

σ e = µ p {(C α p + ρ(C α p Tr[C α p ] -C 2α p )} -pI (6)
Let us give the first component as:

σ e,1 = µ p λ 2α p 1 (1 + ρ(λ 2α p 2 + λ 2α p 3 )) -p (7)

The non-linear Landau model

The nonlinear Landau model is just a nonlinear version of the Hookean model where the standard deformation tensor is replaced by the Green-Lagrange strain tensor, E = 1 2 (F T F-I) = 1 2 (C -I). It seems to have been proposed by Brillouin [START_REF]Brillouin Sur les tensions de radiation (On tensions of radiation)[END_REF] and turns out to be attributed to Landau, nowadays. For this model, the strain energy density is commonly given in term of the Green-Lagrange strain tensor, E = 1 2 (F T F -I) = 1 2 (C -I) as an expansion of third order. It is justified by the fact that it is directly related to the linear Hookean energy. Restricted to the first leading contributions, it reads:

W La = λ 0 2 (Tr[E]) 2 + µ 0 Tr[E 2 ] (8)
where λ 0 and µ 0 are equivalent to the Lamé coefficients in linear elasticity (see [START_REF] Landau | A Course of Theoretical Physics[END_REF]). This expansion is commonly called the Landau model in acoustic [START_REF] Ye | Simulation of non-linear transient elastography: finite element model for the propagation of shear waves in homogeneous soft tissues[END_REF]. It is different from the Mooney-Rivlin expansion for hyper-elasticity since such expansion involves I 1 , I 2 1 and I 2 . Indeed after elementary algebra, it reads:

W La = c 1 4 (I 1 -3) + c 2 2 (I 2 -3) + c 3 8 (I 2 1 -9) (9)
with c 1 = -(2µ 0 + 3λ 0 ), c 2 = -µ 0 and c 3 = (2µ 0 + λ 0 ), while the Mooney-Rivlin model does not includeI 2 1 . In addition, incompressibility in linear elasticity means λ 0 infinite, indicating that this model may appear inadequate for living tissues. It is not difficult to deduce the stress tensor using Eq.( 3):

σ e = λ 0 2 (I 1 -3)C + µ 0 (C 2 -C) -pI (10)
These models treat all directions equivalently, they are isotropic so a priori do not involve fibers. The possible anisotropy of the fiber network which may play a role for fibrotic tissues is discarded. Let us consider now the GOH model.

GOH model

The Neo-Hookean model is recovered for ρ = 0 in the Mooney-Rivlin definition (see Eq.( 1)), for α p = 2 in the oneterm Ogden model and for c 2 = c 3 = 0 in the Landau model. The Mooney-Rivlen model corresponds to α p = 1 for the Ogden model. The discussion for the GOH model is postponed to the next paragraph.

Nonlinear viscoelasticity

According to Destrade et al [START_REF] Destrade | Proper formulation of viscous dissipation for nonlinear waves in solids The[END_REF], some inconsistencies appear in the literature to establish the viscous dissipation tensor with the correct symmetry and the objectivity property although some issues can be found in [START_REF] De Pascalis | On nonlinear viscoelastic deformations: a reappraisal of Fung's quasi-linear viscoelastic model[END_REF]. Special attention is put on this question, focusing on the Landau's model. Based on fundamental principles of nonlinear mechanics, they derive a viscous stress tensor as:

S La,vd = λ0 (Tr( Ė)F + 2µ 0 F Ė (11)
where λ0 and μ0 are chosen by analogy with Eq.( 8) but having different values. Here, we suggest another way to recover this result, perhaps more simple and which can be generalized to arbitrary hyper-elasticity modeling with its extension to viscoelasticity. Considering first the Landau model,we define an energy dissipation rate per unit time, in the same spirit as Eq.( 8):

R La = λ0 2 (Tr[ Ė]) 2 + μ0 Tr[ Ė2 ]) (12)
The full energy variation rate also called the Rayleighian [START_REF] Onsager | s variational principle in soft matter[END_REF][START_REF] Landau | of A Course of Theoretical Physics[END_REF] of the system is then:

R = V dV dW E dt -p J + R Ė (13)
= V dV ( ∂W E ∂E -p ∂J ∂E ) Ė + R Ė (14)
Its extremum, with respect to an arbitrary elastic strain rate requires to consider :

δR = V dV ∂W E ∂E -p ∂J ∂E + ∂R Ė ∂ Ė :δ Ė (15)
The first two terms in Eq.( 15) is nothing else that F -1 S e where S e is the nominal stress. By analogy, we can define the viscous stress by imposing that the last term in Eq.( 15) is equal to F -1 S vd so that the total stress is S = S e + S vd . The integrand becomes F -1 S : δ Ė = Tr[F -1 S δ Ė] since the tensor F -1 S is selfadjoint. We will show now that, by imposing the stationarity of the Rayleighian with respect to any variation of the local velocity, we recover the classical law of the nonlinear elasticity, exactly in the same way as the stationary energy principle. Using the cyclic property of the Trace, it reads:

δR = V dVTr[F -1 S (δ ḞT F + F T δ Ḟ) = V dVTr[S T δ Ḟ] = V dV S : Grad ∂ r ∂t ( 16
)
By standard integration by part for the last integral, the classical equilibrium equation is recovered Div S = 0 but for a dissipative system and the border condition:

F(t) • N = AS • N
, where A is the transverse area of the sample, perpendicular to the normal N. Defining Σ = F(t)/A, the Kelvin-Voigt model is recovered for nonlinear hyper-elastic tissues. In addition, in the case of the Landau's model, the Destrade et al. result is found without difficulty (Eq.( 11)) since:

F -1 S vd = ∂R Ė ∂ Ė = λ0 Tr[ Ė]I + 2 μ0 Ė (17)
With this variational method, based on the extremum of the dissipation rate, no assumption is required. The viscous stress S vd simply results from the first principles of mechanics of dissipative systems with a proper choice of dissipation energy. Concerning hyper-elastic model given in term the Cauchy-Green tensor C a similar dissipation energy rate R Ċ is expected to be given in term of Ċ, then the variational method again applies with this formulation and we get:

σ vd = 2C ∂R Ċ ∂ Ċ or S vd = 2F ∂R Ċ ∂ Ċ (18)
Illustration of this formulation is given in the next section with application to the cuboid stretching.

Viscous dissipative stress tensor and stress distribution

120 for a cuboid under uniaxial tension

Coming back to the Ogden model for the hyper-elastic energy, Eq. ( 4), the dissipation rate can be formulated as follows:

R Og = μp 2α q Tr[ Ċα q ] + ρ 2 Tr[ Ċα q ] 2 -Tr[ Ċ2α q ] ( 19
)
The nominal stresses, deduced from the energy density W Og (Eq.( 5) reads:

S Og = µ p F C α p -1 (1 + ρ(Tr[C α p ] -C α p ) -pF -T (20)
and the viscous stress tensor becomes :

S Og,vd = μp F Ċα q -1 1 + ρ Tr[ Ċα q ] -Ċα q (21)
We consider the simplest case of a cuboid in the vertical direction (called one), with a ramp linear in time λ1 (being the velocity of elongation divided by the initial length of the sample). We assume incompressibility and both boundaries 2 and 3 are free of stress. The stretch λ 2 and λ 3 are equal and the incompressibility imposes the following relations:

λ 2 = λ 3 = 1 √ λ 1 and λ2 = λ3 = - 1 2 λ 2 λ1 λ 1 = - 1 2 λ1 λ 3/2 1 ( 22
)
Let us define each component of the nominal stress by separating the elastic and viscous part such: S Og,i = S ip/λ i , it reads:

S 1 = µ p λ 2α p -1 1 1 + 2ρλ -α p 1 ( 23
)
S 2 = µ p λ -α p +1/2 1 1 + ρ(λ 2α p 1 + λ -α p 1 ) (24)
For the viscous stress tensor, it reads: Cancellation of the total stress in the second and third components, free at the boundaries, gives the value of p and we obtain finally the stress component of interest as:

S 1,vd = μp 2 α q -1 λ α q 1 λα q -1 1      1 + 2 ρ       -λ1 λ 2 1       α q      (25)
S 2,vd = μp λ -1/2 1       -λ1 λ 2 1       α q -1      1 + ρ 2λ 1 λ1 α q +       -λ1 λ 2 1       α q      (26
S Og,1 = S 1 + S 1,vd - 1 λ 3/2 1 (S 2 + S 2,vd ) (27)
So the total stress is decomposed into an elastic and a viscouselastic part which may be the analogous of a non-linear version of the Kelvin-Voigt model. Let us give the result for the Mooney-Rivlin model [START_REF] Mooney | A theory of large elastic deformation[END_REF]:

S 1 = µ p       (λ 1 - 1 λ 2 1 ) + ρ(1 - 1 λ 3 1 )       (28)
while for the viscous stress, we take α q = 2 since α q = 1 has no physical meaning. This emphasizes that, for the dissipation, the degree of nonlinearities is not the same as for the elasticity.

In this case, we derive

S v1 = μp λ1 (λ -4 1 + 2λ 2 1 + ρλ 2 1 λ -8 1 (1 + 8λ 6 1 ) (29)
For the Landau dissipation model we get:

S v1 = λ1 λ 4 1 λ0 (λ 3 1 -1) 2 + μ0 (1 + 2λ 6 1 ) (30)
All these models are isotropic which may be sufficient for low stretch values. However, it may be important to consider also fibrous anisotropic tissues which play an important role in pathologies but also in embryogenesis [START_REF] Clément | Viscoelastic Dissipation Stabilizes Cell Shape Changes during Tissue Morphogenesis[END_REF] as demonstrated recently by laser nano-indentation [START_REF] Vuong-Brender | The interplay of stiffness and force anisotropies drive embryo elongation[END_REF] .

The GOH model

Nonlinearity and visco-elasticity are expected to be stronger for fibrotic tissues. It is why we also present the results for the GOH model [START_REF] Gasser | Hyperelastic modelling of arterial layers with distributed collagen fibre orientation[END_REF] and its possible visco-elastic counterpart. This model induces strong nonlinearities and has been proved to exhibit the correct behavior in a previous work on the same kind of fibrotic tissues [START_REF] Ben Amar | Morpho-elasticity of inflammatory fibrosis: the case of capsular contracture[END_REF]. For multiple families of fibers, it allows to introduce an averaged orientation and a possible dispersion for each family:

W GOH = µ 2 (Tr[C] -3) + Σ j k j,1 2k j,2 e k j,2 [C:H-1] 2 -1 (31)
where k i, j 's are positive dimensionless stiffness parameters. k j,1 is related to the stiffness of the fibers in the small strain regime while k j,2 is related to large strain stiffening behavior of the fibers [START_REF] Destrade | Proper formulation of viscous dissipation for nonlinear waves in solids The[END_REF]. Considering no specific orientation for the fibers, the model simplifies becoming:

W GOH = µ 2 Tr[C] -3 + k 1 2k 2 e k 2 (Tr[C]-3)-1) (32)
We can imagine that the visco-elasticity dissipation rate obeys the same rules giving:

R GOH,vd = μ 8 Tr[ Ċ2 ] + k1 k2 (e k2 (Tr[˙C 2 ]-Ċ2 ) (33)
It is then easy to calculate the viscous stress, for example in the direction 1

S GOH,vd = μλ 2 1 λ1 1 + k1 e 4 k2( ( λ1 λ 1 ) 2 +( λ2 λ 2 ) 2 ) ( 34
)
with equivalent expression for the other components. Since our apparatus imposes ramps of extension, linear in time, the force to apply is given by the initial horizontal surface 130 of the sample times the value of the nominal stress S 1 which is plotted in Fig. [START_REF] Ogden | Biomechanical Modelling at the molecular, cellular and tissue levels CISM[END_REF][START_REF] Ciarletta | Amar A finite dissipative theory of temporary interfibrillar bridges in the extra-cellular matrix of ligaments and tendons[END_REF] for different values of the elastic model and different values for the ramp velocities. Independently of the model used for the dissipative contribution, there is a jump in tension as soon as the loading ramp is applied. This jump at zero static stretch is not observed in our data, hidden by the calibration. If this model may represent conveniently the nonlinear counterpart of the Voigt viscoelasticity model, for the Maxwell model it does not seem so obvious. So we will accept without any more demonstration the linear Maxwell model [START_REF] Belytschko | Nonlinear Finite Elements for Continua and Structures[END_REF] joined 140 to the nonlinear Voigt model. One important question concerns the existence of pre-stretch or pre-stress for these tissues. The pre-stress accumulated by the pathology seems to have disappeared but pre-stretch is possible. In case of pre-stretch one has to revisit our determination of S 1 by a modification of λ 1 by λ 1 /λ p , λ 2 by λ 2 / λ p and λ 3 by λ 3 / λ p .

The nonlinear Maxwell-Kelvin-Voigt Model

A simple glance to our Fig. [START_REF] Gasser | Hyperelastic modelling of arterial layers with distributed collagen fibre orientation[END_REF] indicates that it will be impossible to recover the data data with only the Kelvin-Voigt approach whatever the model. The experiment is realized with a constant ramp so λ1 is a constant. This has been checked by analyzing the data automatically recorded by a camera, the applied load is effectively applied with a constant ramp and the absolute value is the same in the ascending and descending parts. [START_REF] Ye | Simulation of non-linear transient elastography: finite element model for the propagation of shear waves in homogeneous soft tissues[END_REF][START_REF] Landau | of A Course of Theoretical Physics[END_REF]. For clarity, we only show the results for the first and last ascending ramps obtained with a velocity of 0.01 mm/s and the first spindle of curves, with a maximal extension of 0.1 mm. The extension is given in mm and the length of the sample is 1cm, corresponding to a stretch of λ 1 = 1.01. [START_REF] Ye | Simulation of non-linear transient elastography: finite element model for the propagation of shear waves in homogeneous soft tissues[END_REF][START_REF] Landau | of A Course of Theoretical Physics[END_REF]. For clarity, we only show the results for the first and last ascending ramps obtained with a velocity of 0.01 mm/s and for the last spindle of curves, with a maximal extension of 0.3 mm, corresponding to a stretch of λ 1 = 1.03. The extension is given in mm and the length of the sample is 1cm. Analytical results are given by red curves for the first round trip, in green for the last one. The extension is given in mm and the length of the sample is 1cm. Exact values of the parameters can be found in Section 7

As a consequence, the viscous tensor contribution will consist in adding or subtracting a constant since the length extension remains small. So there is no hope that the models presented here may reproduce the data. In addition, there is a difficulty to extract the true value at time t = 0 due to calibration of the apparatus. Whatever the value found at λ = 1, it is obvious that we do not have the degree of freedom necessary to make a cycle. This is true when λ 1 increases but also when it decreases. The Maxwell model which introduces a viscous dissipation damping put in series with our Voigt nonlinear approach, may help to recover the cycles observed in Fig.

η dS dt + S = S e1 + S v1 (1)
Since the left-hand side is directly deduced from linear elasticity, the choice of the Cauchy stress may be preferred although the force acting on the sample is directly proportional to the 150 nominal stress. In [START_REF] Belytschko | Nonlinear Finite Elements for Continua and Structures[END_REF] and references included, the choice is made of the nominal stress and this is the choice we finally adopt in the last section. The ramp being imposed, fixing λ1 , this equation which as one degree of freedom allows to fix one value for the stress S for the ascending ramp and for the decreasing ramp, then a cycle can be recovered giving hope to recover the data. Only a good fitting with a limited number of parameters can confirm the strategy. Notice that if we can estimate the parameters which intervene in the elastic and viscoelastic tensors, it is not the case for η. In particular we do

Data versus Modelling

In this section, preliminary results on cyclic deformations are presented for small extension of fibrotic and surgical tissues. All the spectra: uniaxial stress/extension realized at very small velocities 0.01mm/s exhibit the same appearance as in Fig( 1) both for fibrotic and artificial surgical tissues. Our samples consist in surgical fibrotic mammary tissues [START_REF] Ben Amar | Morpho-elasticity of inflammatory fibrosis: the case of capsular contracture[END_REF] which vary from 170 one patient to another, depending on the grade of the pathology and the human variability. An extensive study will be the subject of future work. What is known is that these tissues contain a very disorganized collagenous structure. However, the sample treated here presents generic features which is not shared by elastomers for example, with the same dimensions. The details on the apparatus and experimental procedure can be found in [38]. The main difficulty of such tensional test are drainage effects consisting of exchange of water between the tissue and the bath which will spoil the measurements. So the ramp velocity 180 is chosen very small to avoid these exchanges difficult to evaluate quantitatively. As shown in Fig. [START_REF] Gasser | Hyperelastic modelling of arterial layers with distributed collagen fibre orientation[END_REF], nonlinearities of the stress/strain answer are very strong, even for a very low stretch value, less than 1.01. This is an amazing result but it has been observed on 4 samples. Indeed, if it strongly suggests dissipation, standard models as shown in Fig. [START_REF] Gasser | Hyperelastic modelling of arterial layers with distributed collagen fibre orientation[END_REF] show more easily an increase of the stress in the negative ramp compared to the positive ramp. One may think of the so-called Mullins effect which has been introduced in the biomechanical literature to explain hysteresis. It is based on an internal degree of freedom which varies with the applied stretch. An illustration for rubber-like materials has been given by Tommasi et al. ([39,[START_REF] Tommasi | A micromechanics-based model for the Mullins effect[END_REF]). These authors consider the microscopic structure of polymeric gels which can breakdown under stretching. This obviously modifies the elastic energy of the samples which contains less fibers due to the breaking of the chains. In [START_REF] Ben Amar | Morpho-elasticity of inflammatory fibrosis: the case of capsular contracture[END_REF], such a scenario is also investigated with eventually a reorientation of the fiber network in the direction of the stretching. However all these changes have little chance to occur for the rather low stretch value, shown in Fig. [START_REF] Gasser | Hyperelastic modelling of arterial layers with distributed collagen fibre orientation[END_REF] and for these models, the stress always 200 cancel for zero stretch. Although, there is an uncertainty on the real zero in our experiment, these models cannot explain the systematic shift towards low stress values. Only dissipation can explain the gap observed at low vanishing stretching since all models give a value proportional to λ1 . Unfortunately, the calibration of the zero appears to be uncertain and so we use the value at zero stretch as a parameter to include in the fitting process. The final results are shown in Figs. [START_REF] Pioletti | On the independence of time and strain effects in the stress relaxation of ligaments and tendons[END_REF][START_REF] Shearer | Sequential straightening 280 and loading viscoelasticity: A model of the time-dependent behaviour of tendon fascicles incorporating fibril recruitment[END_REF][START_REF] Natali | Anisotropic elasto-damage constitutive model for the biomechanical analysis of tendons[END_REF][START_REF] Muliana | A concurrent micromechanical model for predicting nonlinear viscoelastic responses of composites reinforced with solid spherical particles[END_REF][START_REF] Legoff | A Global Pattern of Mechanical Stress Polarizes Cell Divisions and Cell Shape in the Growing Drosophila Wing Disc[END_REF] allowing the following conclusions:

• At extremely stretch values, hysteretic cycles on fibrotic • The best models for the dissipation energy rates is given by 220 Eqs. [START_REF] Ye | Simulation of non-linear transient elastography: finite element model for the propagation of shear waves in homogeneous soft tissues[END_REF][START_REF] Landau | of A Course of Theoretical Physics[END_REF], so the dissipation based on the Neo-Hookean formulation. In particular, the Landau model, according to Eq.( 30) (with λ0 0) gives relatively unstable results.

• A pre-stretch is necessary, larger than than one during a descending ramp, after an ascending ramp, smaller than one after a descending ramp. Although remaining very close to one, the quality of the fitting depends tremendously on its value.

• The coefficient η of Eq.( 35) is different for ascending/descending ramps but remains stable during the nine examples shown in fig [START_REF] Legoff | A Global Pattern of Mechanical Stress Polarizes Cell Divisions and Cell Shape in the Growing Drosophila Wing Disc[END_REF].

Finally, here are the results for the sample we have selected, beginning by the intrisic values of the material: the numerical parameters µ p (Eq.( 28) and μp (Eq.( 29) that we can estimate are: µ p = 10 5 Pas and μp = 0.1710 8 Pas.s [START_REF] Vuong-Brender | The interplay of stiffness and force anisotropies drive embryo elongation[END_REF] while for the viscous parameter according to Eq.( 35), we get two different values, one for the ascending part, the other for the descending part as for both pre-stretch values:

η A = 0.44s; η D = 0.08s; λ p,A = 0.955 and λ p,D = 1.006 [START_REF] Belytschko | Nonlinear Finite Elements for Continua and Structures[END_REF] So, in this example, elasticity and visco-elasticity are in competition but nonlinearities for the stresses (right-hand side of Eq.( 35) can be discarded. The stress/strain answer is essentially sensitive to the η coefficient and the Maxwell viscosity. The pre-strain is probably due to a weak effect of orientation/ disorientation when the sample is stretched or compressed. However, it is a parameter which dramatically influences the agreement between the model and the experiment. Here we give an averaged value since it very slightly changes from one set of curves 240 from another.

Conclusion

This approach is based on traditional models of biomechanics. Other modelings can be found and complexity can be added via the introduction of internal variables, maybe useful at larger stretch for fibrotic tissues. The theory of nematic elastomers have pointed out that the fibers play a crucial role in the sample response to unidirectional extension. It is complicated by the fact that, for incompressible materials, the extension in one direction induces automatically compression in one or two other 250 directions with a different answer according to fiber orientations. In addition, if the fibers have a given preferred direction in average, this one may change during loading, becoming an internal variable. It is highly appreciated if measurements give values which are typical of the material and not of the measure itself, but in this simple set-up, measurements affect the internal structure even at rather low stretch value. A more systematic study of these tissues (either fibrotic or artificially reconstructed for plastic surgery),at larger extension values, will be the subject of future work.

Acknowledgments

MBA would like to thank the Isaac Newton Institute for Mathematical Sciences, Cambridge, for support and hospitality during the programme "Growth, Form and Self-Organization" where work on this paper was undertaken. MBA acknowledges the partial support from the Simons Foundation, the EPSRC grant no EP/K032208/1 and from ANR under the contract MecaTiss (ANR-17-CE30-0007).

MBA thanks Khawla Sleiman and Fréderic Lechenault for their participation to the experiment and Dr Michael Atlan for providing surgical tissues.

30

 30

Figure 2 :

 2 Figure 2: Nominal stress S 1 as a function of the vertical stretch in the Ogden model given by Eqs.(4,10). Being a scaling parameter, µ p is chosen to 1. Notice the change of curvature in the stress answer for rather large exponents due to nonlinearities but in this example, the extension is ten times the values reached in Fig.(1)

)Figure 3 :

 3 Figure 3: Nominal stress S 1 taking into account both the elastic and viscoelastic stresses according to Eq.[START_REF] Destrade | Proper formulation of viscous dissipation for nonlinear waves in solids The[END_REF][START_REF] Onsager | s variational principle in soft matter[END_REF]. The linear ramp is λ1 = 1, the viscous coefficients are: λ0 = μ0 = 0.5.

Figure 4 :Figure 5 :Figure 6 :

 456 Figure4: The 9 stress/extension curves obtained experimentally for one specific sample of mammary fibrotic tissue. The measurement corresponds to the nominal stress. For clarity, we only show the results for ascending ramp obtained with a velocity of 0.01 mm/s. The extension is given in mm and the length of the sample is 1cm.

Figure 7 :

 7 Figure 7: Superposition of the nominal stress obtained experimentally and theoretically in the Neo-Hookean model, see Eqs.[START_REF] Ye | Simulation of non-linear transient elastography: finite element model for the propagation of shear waves in homogeneous soft tissues[END_REF][START_REF] Landau | of A Course of Theoretical Physics[END_REF]. For clarity, we only show the results for the first and last ascending ramps obtained with a velocity of 0.01 mm/s and for the last spindle of curves, with a maximal extension of 0.3 mm, corresponding to a stretch of λ 1 = 1.03. The extension is given in mm and the length of the sample is 1cm.

Figure 8 :

 8 Figure 8: The stress/extension curves obtained in the first and third rounds of the last set (see Fig.(1). Experimental data are represented by black lines.Analytical results are given by red curves for the first round trip, in green for the last one. The extension is given in mm and the length of the sample is 1cm. Exact values of the parameters can be found in Section 7

160

  not know if η depends on the whole apparatus or on the tissue mechanical properties. η has the dimension of a time. Let us now analyze the data shown in Fig.[START_REF] Gasser | Hyperelastic modelling of arterial layers with distributed collagen fibre orientation[END_REF].

  210tissues are explained by a competition between the elastic energy and the dissipation. • In this range of forcing values, the Neo-Hookean model, possibly improved by the Mooney-Rivlin,gives the most stable results between different cycles and different extensions (see the agreement between data and the Neo-Hookean models in Fig.(6) for ascending ramps, in Fig.(7 ) for descending ramps and for two cycles of the last set of curves,Fig.(8).

  260

���� ���� ���� ���� ���� ���� ���� -��� -��� ��� ��� ��� ��� �� � � �� � ��

	Figure 1: A typical set of measurements, force versus extension, concerning a
	surgical fibrotic tissue. Several cycles are made with the same maximal exten-
	sion but with 3 different values for the final extension. There is a systematic
	shift towards lower values as the number of cycles increases.