Benoît Perthame
email: benoit.perthame@upmc.fr

Nicolas Vauchelet
email: nicolas.vauchelet@upmc.fr

Incompressible limit of mechanical model of tumor growth with viscosity

Keywords: Tumor growth, Hele-Shaw equation, Free boundary problems, Porous media, Viscoelastic media Mathematical Classification numbers: 35J60, 35K57, 74J30, 92C10

Various models of tumor growth are available in the litterature. A first class describes the evolution of the cell number density when considered as a continuous visco-elastic material with growth. A second class, describes the tumor as a set and rules for the free boundary are given related to the classical Hele-Shaw model of fluid dynamics.

Following the lines of previous papers where the material is described by a purely elastic material, or when active cell motion is included, we make the link between the two levels of description considering the 'stiff pressure law' limit. Even though viscosity is a regularizing effect, new mathematical difficulties arise in the visco-elastic case because estimates on the pressure field are weaker and do not imply immediately compactness. For instance, traveling wave solutions and numerical simulations show that the pressure may be discontinous in space which is not the case for the elastic case.

The cell model with visco-elastic flow

We consider a mechanical model of tumor growth considered as a visco-elastic media. We denote the number density of tumor cells by n(x, t), the pressure by p(x, t) and we assume a Brinkman flow that means the macroscopic velocity field is given by ∇W for a potential W closely related to the pressure. With these assumptions, the model for tumor growth writes

∂ t n k -div(n k ∇W k ) = n k G p k , (1)
-ν∆W k + W k = p k (x, t) := Π k (n k ), (2)
where we choose the pressure law given by:

Π k (n) = k k -1 n k-1 , k > 2. (3)
Following [START_REF] Byrne | Individual-based and continuum models of growing cell populations: a comparison[END_REF][START_REF] Ranft | Fluidization of tissues by cell division and apoptosis[END_REF], we assume that growth is directly related to the pressure through a function G(•) that satisfies

G ∈ C 1 (R), G ′ (•) ≤ -α < 0, G(P M ) = 0 for some P M > 0. (4)
The pressure P M is usually called the homeostatic pressure . We complete equation ( 1), [START_REF] Allaire | Homogenization of the Navier-Stokes equations and derivation of Brinkman's law[END_REF] with a family of initial data n 0 k satisfying (for some constant C independent of k)

0 ≤ n 0 k , Π k (n 0 k ) ≤ p M , n 0 k L 1 (R d ) ≤ C. (5)
The viscosity coefficient, ν > 0, is supposed to be constant; when viscosity is neglected, that means equation [START_REF] Allaire | Homogenization of the Navier-Stokes equations and derivation of Brinkman's law[END_REF] with ν = 0, we recover Darcy's law for which an important literature is available, see [START_REF] Chatelain | Emergence of microstructural patterns in skin cancer: a phase separation analysis in a binary mixture[END_REF][START_REF] Sherratt | A new mathematical model for avascular tumour growth[END_REF][START_REF] Friedman | A hierarchy of cancer models and their mathematical challenges[END_REF][START_REF] Roose | Mathematical models of avascular tumor growth[END_REF][START_REF] Cristini | Nonlinear simulations of solid tumor growth using a mixture model : invasion and branching[END_REF][START_REF] Lowengrub | Nonlinear modelling of cancer: bridging the gap between cells and tumours[END_REF][START_REF] Perthame | The Hele-Shaw asymptotics for mechanical models of tumor growth[END_REF][START_REF] Perthame | Derivation of a Hele-Shaw type system from a cell model with active motion[END_REF] and the references therein. In that case only friction with the cell surrounding (extra-cellular matrix) is considered. Viscosity is a way to represent friction between cells themselves, considered as a Newtonian fluid and Brinkman's law has been derived rigorously for inhomogenous materials [START_REF] Allaire | Homogenization of the Navier-Stokes equations and derivation of Brinkman's law[END_REF]. Viscoelastic models for tumor growth, based on Stokes' or Brinkman's law have also been used in the context of tumor growth is [START_REF] Zheng | Nonlinear simulation of tumor necrosis, neovascularization and tissue invasion via an adaptive finite-element/level-set method[END_REF][START_REF] Bresch | A viscolelastic model for avascular tumor growth[END_REF][START_REF] Ranft | Fluidization of tissues by cell division and apoptosis[END_REF] with a major difference, namely the pressure does not follow a law-of-state (4) but follows from the tissue incompressibility. However, Stokes' or Brinkman's law are also used considering the tissue as 'compressible' [START_REF] Bittig | Dynamics of anisotropic tissue growth[END_REF][START_REF] Basan | Homeostatic competition drives tumor growth and metastasis nucleation[END_REF]. To use Laplacian in [START_REF] Allaire | Homogenization of the Navier-Stokes equations and derivation of Brinkman's law[END_REF], rather than Stokes viscosity terms, is to simplify the presentation and presentation of the mathematical ideas. Indeed, this is not central for our aim here, which is to explain the derivation of such 'incompressible' models from the 'compressible' equations. Note that the theory of mixtures allows for a general formalism containing both Darcy's law and Brinkman's law [START_REF] Byrne | A two-phase model of solid tumor growth[END_REF][START_REF] Ambrosi | On the closure of mass balance models for tumor growth[END_REF][START_REF] Preziosi | Multiphase modelling of tumour growth and extracellular matrix interaction: mathematical tools and applications[END_REF].

Our interest is in the 'stiff pressure law' limit of this model towards a free boundary model which generalizes the classical Hele-Shaw equation. That is the limit k → ∞ and we first explain formally what can be expected. The limit uses strongly the equation satisfied by the pressure. Multiplying equation (1) by Π ′ k (n) and using the chain rule, we deduce

∂ t p k -n k Π ′ k (n k )∆W k -∇p k • ∇W k = n k Π ′ k (n k )G(p k ).
From our choice for the law of state (3), we deduce that

nΠ ′ k (n) = kn k-1 = (k -1)Π k (n).
Injecting this expression into the above equation, we deduce that

   ∂ t p k -∇p k • ∇W k = k-1 ν p k Q k , Q k = W k -p k + νG(p k ) := W k -H -1 (p k ) (6)
where we have defined the function H, coming with some properties, as

H := (I -νG) -1 , p m := H(0) > 0; H is increasing, H ′ (•) < 1. (7)
Indeed, G is non-increasing and thus (I -νG) is invertible on [0, P M ] onto [-νG(0), P M ]. Furthermore, notice that (I -νG) ′ > 1.

Back to the limit k → +∞, at least when p k converges strongly, from (3), we first find the relation

p ∞ (1 -n ∞ ) = 0. ( 8
)
Letting k → +∞ and asuming we can pass into the limit in all terms, we formally deduce

p ∞ ∆W ∞ + G(p ∞ ) = 0.
Therefore, at the limit we can distinguish between two different regions. The first region is defined by the set

Ω(t) := {p ∞ (•, t) > 0} (9)
on which we have the system :

n ∞ = 1, (10)
-ν∆W ∞ + W ∞ = p ∞ , (11)
∆W ∞ + G(p ∞ ) = 0. ( 12
)
Thus the latter system reduces to :

n ∞ = 1, p ∞ = H(W ∞ ), x ∈ Ω(t), -ν∆W ∞ + W ∞ -H(W ∞ ) = 0.
On the second region, R d \ Ω(t), the limiting system writes

p ∞ = 0, ∂ t n ∞ -div(n ∞ ∇W ∞ ) = n ∞ G(0), -ν∆W ∞ + W ∞ = 0.
To establish rigorously this limit, we need some additional assumption on the initial data. Namely, we need that the family n 0 k is 'well-prepared'. By this, we mean that, for some open set Ω 0 ,

Π k (n 0 k ) -→ k→∞ p 0 ∞ = H(W ∞ ) a.e. in Ω 0 , n 0 k = 0 in R d \Ω 0 . ( 13
)
Note that, with the notation in [START_REF] Bittig | Dynamics of anisotropic tissue growth[END_REF], this assumption implies that Q 0 k ≡ 0 and n 0 k -→ k→+∞ 1 Ω 0 . For this purpose, the latter assumption can be slightly relaxed to n 0 k ≪ e -A/k for all A > 0 in R d \Ω 0 . With our present proof, we need to avoid the existence of a domain where n 0 k remains strictly between 0 and 1, a case which we leave open at this stage.

Our goal is to prove the Theorem 1.1 Under assumptions (4), ( 5) and (13), consider a solution of the system (1)- [START_REF] Ambrosi | On the closure of mass balance models for tumor growth[END_REF]. After extraction of subsequences, both the density n k and the pressure p k converge strongly in L 1 loc (0, T ) × R d , for all T > 0, as k → +∞ towards respectively n ∞ and p ∞ belonging to L 1 ∩ L ∞ (0, T ) × R d ; up to a subsequence, W k converges strongly in L 1 (0, T ), W 1,q loc (R d ) , for all q ≥ 1, towards W ∞ . Moreover, these functions satisfy

∂ t n ∞ -div(n ∞ ∇W ∞ ) = n ∞ G(p ∞ ), n ∞ (t = 0) = n 0 ∞ = 1 {Ω 0 } , (14)
-ν∆W ∞ + W ∞ = p ∞ , (15)
p ∞ = H(W ∞ )1 {p∞>0} , p ∞ (1 -n ∞ ) = 0, ( 16
)
p ∞ p ∞ -W ∞ -νG(p ∞ ) = 0, a.e. ( 17
)
The first relation in ( 16) is equivalent to the statement [START_REF] Perthame | Existence of solutions of the hyperbolic Keller-Segel model[END_REF] and replaces the usual 'complementary relation' in Hele-Shaw flow, p ∞ (∆p ∞ + G(p ∞ )), see [START_REF] Perthame | The Hele-Shaw asymptotics for mechanical models of tumor growth[END_REF][START_REF] Perthame | Derivation of a Hele-Shaw type system from a cell model with active motion[END_REF][START_REF] Elliotta | A variational inequality approach to Hele-Shaw flow with a moving boundary[END_REF].

Because the function H(•) does not vanish, we conclude from the first relation in [START_REF] Lowengrub | Nonlinear modelling of cancer: bridging the gap between cells and tumours[END_REF], that p ∞ is discontinuous. This is a major difference with elastic materials (Darcy's law), then p ∞ is continous in space, and this is illustrated by traveling wave solutions we build in Section 3. The pressure jump is however related to the potential W ∞ , a difference with models including surface tension where the jump is related to the free boundary curvature, see [START_REF] Alikakos | Convergence of the Cahn-Hilliard Equation to the Hele-Shaw Model[END_REF][START_REF] Escher | Classical solutions for Hele-Shaw models with surface tension[END_REF] and the reference therein.

We first prove Theorem 1.1 in several steps. In a first step, we derive a priori estimates. Because they do not give compactness for the pressure, we analyze possible oscillations using a kinetic formulation. From properties of solutions of the corresponding kinetic equation, we conclude that strong compactness occurs. All these steps are in Section 2. The one dimensional traveling wave profiles are presented in Section 3 with numerical illustrations. The final Section is devoted to a conclusion and presentation of some perspective.

Proof of the Hele-Shaw limit

We divide the proof of our main result Theorem 1.1 in several steps. We begin with several bounds which are useful for the sequel. Then, in order to prove strong convergence of the pressure p k , we analyze possible oscillations using the kinetic formulation of (6) in the spirit of [START_REF] Perthame | Existence of solutions of the hyperbolic Keller-Segel model[END_REF].

Estimates

Lemma 2.1 (A priori estimates) Under previous assumptions, for all T > 0, the uniform bounds with respect to k hold

n k , p k and W k ∈ L ∞ (0, T ); L 1 ∩ L ∞ (R d ) , p k ≤ P M , W k ∈ L ∞ (0, T ); W 1,q (R d ) , for 1 ≤ q ≤ ∞, D 2 W k ∈ L ∞ (0, T ); L q (R d ) , for 1 < q < ∞, ∂ t W k ∈ L 1 (0, T ); L q (R d ) , for 1 ≤ q ≤ ∞, ∂ t ∇W k ∈ L 1 (0, T ); L q (R d ) , for 1 < q < d d-1 .
For some nonnegative constant C(T ), independent of k, we have

k T 0 R d p k p k -W k -νG(p k ) dxdt ≤ C(T ). ( 18
)
We can draw several consequences of this Lemma. First, after extracting subsequences, it is immediate that the following convergences hold as k → ∞:

n k → n ∞ ≤ 1, p k → p ∞ ≤ P M weakly -⋆ in L ∞ (0, +∞) × R d ,
and these limits belong to L ∞ (0, T ); L 1 (R d ) for all t > 0. Also, we have

W k → W ∞ < P M , ∇W k → ∇W ∞ locally in L q (0, T ) × R d , 1 ≤ q < ∞.
Passing to the limit in ( 2) and in the left hand side of (1), we get

-ν∆W ∞ + W ∞ = p ∞ . ( 19
)
The second consequence concerns the backward flow with velocity ∇W k defined as

d ds X (k) (x,t) (s) = -∇W k (X (k) (x,t) (s), s), X (k) (x,t) (t) = x, (20)
as well as the forward flow

d dt Y (k) (x) (t) = -∇W k (Y (k) (x) (t), t), Y (k) (x) (t = 0) = x. (21)
Even though, ∇W k is not uniformly Lipschitz continuous but slightly less, and according to DiPerna-Lions theory [START_REF] Diperna | Ordinary differential equations, transport theory and Sobolev spaces[END_REF], these flow are well defined a.e. and, after extraction of subsequences as in Lemma 2.1, it converges a.e. to the limiting flows defined by (41) for the backward flow and by (28) for the forward flow.

The third conclusion uses a combination of the above flow with equation ( 6). We have

p k (x, t) = 0 for x ∈ R d \Ω k (t), Ω k (t) = Y (k) (t)[Ω 0 ] (22) Proof. 1st step. A priori bounds in L 1 ∩ L ∞ . Clearly n k is nonnegative provided n k (t = 0) ≥ 0. Inte- grating, we deduce a bound for n k in L ∞ (0, T ); L 1 (R d ) , uniformly with respect to k.
By definition of p k in (3), we clearly have that Π ′ k (n k ) ≥ 0 when k > 1. We can apply the maximum principle of [START_REF] Tang | Composite waves for a cell population system modelling tumor invasion[END_REF]Lemma 2.1] to obtain the uniform bound

0 ≤ p k ≤ P M .
Therefore, still using relation (3), we have

n k = k-1 k p k 1/(k-1) and n k is uniformly bounded in L ∞ (0, +∞) × R d . Then, writing p k ≤ n k n k k-2 ∞ , we deduce an uniform bound of (p k ) k in L ∞ (0, T ); L 1 (R d ) . 2nd step. Representation of W k . Using elliptic regularity on (2), we conclude that for all t ∈ [0, T ], W k (t, •) is bounded in W 2,q (R d ).
Moreover, denoting by K the fondamental solution of -ν∆K + K = δ 0 , we have

W k = K ⋆ p k , K(x) = 1 4π ∞ 0 e -π |x| 2 4sν + s 4π ds s d/2 . ( 23
)
We recall that

   K ∈ L q (R d ), ∀1 ≤ q < d d-2 , (1 ≤ q ≤ +∞ for d = 1), ∇K ∈ L q (R d ), ∀1 ≤ q < d d-1 , and that K ≥ 0, R d K(x) dx = 1, which we use below.
Taking the convolution of ( 6), we deduce

∂ t W k = K ⋆ [∇p k • ∇W k + k -1 ν p k Q k ]. ( 24
)
3rd step. Bounds on Q k . Then, by definition of Q k and using ( 6), we compute

∂ t Q k -∇Q k • ∇W k + k -1 ν 1 -νG ′ (p k ) p k Q k = -|∇W k | 2 + K ⋆ [∇p k • ∇W k + k -1 ν p k Q k ].
Therefore, from a standard computation, we deduce

∂ t |Q k | -∇|Q k |.∇W k + k -1 ν 1 -νG ′ (p k ) p k |Q k | ≤ |∇W k | 2 + |K ⋆ [∇p k .∇W k ]| + k -1 ν K ⋆ [p k |Q k |].
We may integrate in x and t. Because p k and W k are uniformly bounded in L 1 ∩ L ∞ , and |G ′ | ≥ α from (4), we find

α(k -1) T 0 R d p k |Q k | dxdt ≤ R d |Q k (x, 0)|dx - R d |Q k (x, T )|dx + T 0 R d |∇W k | 2 dxdt + T 0 R d -|Q k | ∆W k + |K ⋆ [∇p k • ∇W k ]| dxdt.
The three first terms in the right hand side are all controlled uniformly and, to conclude the bound [START_REF] Perthame | Derivation of a Hele-Shaw type system from a cell model with active motion[END_REF], we have to estimate last two terms. Using (2), the first term is

- T 0 R d |Q k |∆W k dxdt = 1 ν T 0 R d |Q k | (p k -W k ) dxdt ≤ 1 ν T 0 R d |Q k | p k dxdt,
and this term is controlled, for k large enough, by the αk term in the left hand side. The second term is

K ⋆ [∇p k • ∇W k ] = ∇K ⋆ [p k • ∇W k ] -K ⋆ [p k ∆W k ].
Using the uniform bounds on p k , we have that

p k • ∇K * p k is uniformly bounded, with respect to k, in L ∞ (0, T ); L q (R d ) , 1 ≤ q ≤ ∞, and thus, ∇K ⋆ [p • ∇W k ] is also uniformly bounded in L ∞ (0, T ); L q (R d ) , 1 ≤ q ≤ ∞. Finally, p k ∆W k is also uniformly bounded in L ∞ (0, T ); L q (R d ) , 1 ≤ q ≤ ∞.
This immediately concludes the proof of estimates [START_REF] Perthame | Derivation of a Hele-Shaw type system from a cell model with active motion[END_REF].

4th step. Estimate on ∂ t W k . Finally, using the above estimate and equation ( 24), we deduce that

∂ t W k is uniformly bounded with respect to k in L ∞ (0, T ); L q (R d ) , 1 ≤ q ≤ ∞.
For the estimate for ∂ t ∇W k , we can use again the above calculation and write

∂ t ∇W k = -D 2 K ⋆ [p k • ∇W k ] + ∇K ⋆ [p k ∆W k ] + k -1 ν ∇K * [p k Q k ].
Since D 2 K is a bounded operator in L 1 , we conclude the last bound in Lemma 2.1.

Which oscillations for the pressure?

We deduce from Lemma 2.1 that, up to a subsequence, the sequence (W k ) k converges strongly in L 1 ((0, T ), W 1,q loc ). However, we only get weak convergence for the pressure (p k ) k and the density (n k ) k . Here, we give an argument showing that the only obstruction to strong compactness, is oscillations of p k between the values p k ≈ 0 and p k ≈ H(W ∞ ). Lemma 2.2 Let T > 0 and let H be defined in [START_REF] Bresch | A viscolelastic model for avascular tumor growth[END_REF] with the assumptions (4). Consider real numbers β 1 > 0, β 2 > 0 small enough, and let p k be as in Lemma 2.1, then we have

meas β 1 ≤ p k (x, t) ≤ H(W ∞ (x, t)) -β 2 -→ k→+∞ 0, meas p k (x, t) ≥ H(W ∞ (x, t)) + β 2 -→ k→+∞ 0,
where meas denotes the Lebesgue measure.

Proof. Let 0 < β 1 < β 2 < p m , p m being defined in (7), we have for all k ∈ N T 0 R d 1 {β 1 ≤p k ≤H(W∞)-β 2 } dxdt ≤ T 0 R d p k β 1 1 {β 1 ≤p k ≤H(W∞)-β 2 } dxdt. (25)
From assumption (4), the function I -νG is increasing and by definition ( 7), (

I -νG)(H(W ∞ )) = W ∞ ≥ 0 (the nonnegativity is because W ∞ is a solution of (19)). Therefore, on the set {p k ≤ H(W ∞ ) -β 2 }, we have, for some ω(β 2 ) > 0, (I -νG)(p k ) ≤ (I -νG)(H(W ∞ ) -β 2 ) ≤ W ∞ -ω(β 2 ), W ∞ -(I -νG)(p k ) ≥ ω(β 2 ).
Thus we can estimate

T 0 R d p k 1 {β 1 ≤p k ≤H(W∞)-β 2 } dxdt ≤ 1 ω(β 2 ) {β 1 ≤p k ≤H(W∞)-β 2 } p k |(I -νG)(p k ) -W ∞ | dxdt ≤ 1 ω(β 2 ) T 0 R d p k |(I -νG)(p k ) -W k | dxdt + {β 1 ≤p k } p M |W ∞ -W k |dxdt .
Additionally, using estimate [START_REF] Perthame | Derivation of a Hele-Shaw type system from a cell model with active motion[END_REF], and the strong convergence of W k , we deduce that

lim k→+∞ T 0 R d p k 1 {β 1 ≤p k ≤H(W∞)-β 2 } dxdt = 0. ( 26
)
We notice, for future use, that in the same spirit we also have that

lim k→+∞ T 0 R d p k 1 {p k ≤H(W∞)-β 2 } dxdt = 0. ( 27
)
Thus estimates ( 25)-( 26) prove the first statement of Lemma 2.2.

The second statement can be proved in the same way.

Strong convergence of the pressure

However, we need strong convergence to recover the asymptotic limit, in particular the equation satisfied by p ∞ . A difficulty here is that we do not have estimates on the derivatives on p, unlike in [START_REF] Perthame | The Hele-Shaw asymptotics for mechanical models of tumor growth[END_REF][START_REF] Perthame | Derivation of a Hele-Shaw type system from a cell model with active motion[END_REF]. Then we develop another strategy based on estimate [START_REF] Perthame | Derivation of a Hele-Shaw type system from a cell model with active motion[END_REF] to obtain the following strong convergence result :

Lemma 2.3 (Strong convergence of p k ) Up to a subsequence, p k converges strongly locally in L 1 (0, T ) × R d towards p ∞ . Moreover, p ∞ = H(W ∞ )1 {p∞>0} a.e.
Furthermore, we have

Ω(t) = {p ∞ (•, t) = H(W ∞ (•, t))} = R d \{p ∞ (•, t) = 0}
is the image of Ω 0 by the limiting flow Y (x) (t), defined by

d dt Y (x) (t) = -∇W ∞ (Y (x) (t), t), Y (x) (t = 0) = x. (28)
Finally, we have for all T > 0,

k T 0 R d p k (x, t)|Q k (x, t)|dxdt -→ k→+∞ 0. ( 29
)
Proof. The strategy is to pass to the limit in the equation ( 6) for p k and to combine this information with the possible oscillations of p k as described by Lemma 2.2. For that, we need a representation of the weak limit of p k which we can obtain thanks to a kinetic representation.

1st step. Representation of nonlinear weak limits. Our first result is that there is a measurable function 0 ≤ f (x, t) ≤ 1 such that for all smooth function S : [0, ∞) → R, we have, up to a subsequence,

S(p k ) ⇀ k→+∞ S(0)(1 -f ) + S(H(W ∞ ))f, (30)
and

S(0)(1 -f ) + S(H(W ∞ ))f = ∞ 0 S ′ (ξ)χ(ξ) dξ + S(0), χ(x, ξ, t) = f (x, t)1 {0<ξ<H(W∞(x,t))} . (31)
Interpreted in terms of Young measures, this means that p k oscillates between the values 0 and H(W ∞ (x, t)) with the weights 1 -f (x, t) and f (x, t). Notice that for S(p) = p, we find

p ∞ = f H(W ∞ ). ( 32
)
To prove these results, we define

χ k (x, ξ, t) = 1 {0<ξ<p k (x,t)}
and we write

S(p k ) -S(0) = ∞ 0 S ′ (ξ)χ k (x, ξ, t)dξ. (33)
We can extract a subsequence, still denoted (p k ) k , such that

1 {0<ξ<p k } converges in L ∞ ((0, ∞) × R d ) -
weak⋆ towards a function χ(x, ξ, t) which satisfies 0 ≤ χ(x, ξ, t) ≤ 1. Then S(p k ) converges weakly to S(0)

+ ∞ 0 S ′ (ξ)χ(x, ξ, t)dξ.
We define,

f (x, t) = w-lim 1 {p k (x,t)≥pm/2}
where we recall that p m is defined in [START_REF] Bresch | A viscolelastic model for avascular tumor growth[END_REF]. Since H(W ∞ ) > p m , we may use Lemma 2.2 to conclude (30)-(31).

2nd step. Equation satisfied by χ k . We use the equation ( 6)

∂ t p k -∇p k • ∇W k = k -1 ν p k Q k , Q k = W k -p k + νG(p k ).
For any function S ∈ C 2 (R; R), multiplying it by S ′ (p k ) leads to

∂ t S(p k ) -∇S(p k ) • ∇W k = (k -1)p k Q k S ′ (p k ).
Denoting δ the Dirac mass, we can rewrite the later equation as

∂ t ∞ 0 S ′ (ξ)χ k dξ -∇ ∞ 0 S ′ (ξ)χ k dξ • ∇W k = ∞ 0 S ′ (ξ)µ k (x, ξ, t) dξ, (34)
µ k (x, ξ, t) := k -1 ν p k Q k δ {ξ=p k } = k -1 ν p k [W k -p k + νG(p k )]δ {ξ=p k } . (35)
Eliminating the test function S ′ (•), this is equivalent to write

∂ t χ k -∇χ k • ∇W k = µ k . (36)
However, this formula is not enough to pass to the limit k → ∞ and we need the divergence form,

∂ t S(p k ) -div S(p k )∇W k ] + S(p k ) W k -p k ν = (k -1)p k Q k S ′ (p k ) = ∞ 0 S ′ (ξ)µ k (dξ).
Therefore, using (33) and the fact that

S(p k )p k = ∞ 0 S(ξ) + ξS ′ (ξ) χ k dξ, we have ∞ 0 S ′ (ξ) ∂ t χ k -div[χ k • ∇W k ] + χ k W k -ξ ν dξ - ∞ 0 S(ξ) -S(0) ν χ k dξ = ∞ 0 S ′ (ξ)µ k (dξ). (37) Because χ k (ξ) = -∂ ∂ξ ∞ ξ χ k (x, η, t
)dη, and integrating by parts, we have

∞ 0 S(ξ) -S(0) ν χ k dξ = ∞ 0 S ′ (ξ) ν ∞ ξ χ k (x, η, t)dηdξ.
Therefore, (37) is equivalent to our final formulation

∂ t χ k -div[χ k ∇W k ] + χ k W k -ξ ν - 1 ν ∞ ξ χ k (x, η, t)dη = µ k . ( 38
)
One can simplify this relation and write

∂ t χ k -div[χ k ∇W k ] + χ k W k -ξ ν - (p k -ξ) + ν = µ k .
Finally, (37) is equivalent to

∂ t χ k -div[χ k ∇W k ] + χ k W k -p k ν = µ k .
In particular, integrating in ξ we recover the expected formula

∂ t p k -div[p k ∇W k ] + p k ν [W k -p k ] = µ k dξ.
3rd step. Equation satisfied by f . We may pass to the limit in (38). For all T > 0, the sequence

µ k is uniformly bounded in L 1 (R d × R × [0, T ]
) thanks to estimate [START_REF] Perthame | Derivation of a Hele-Shaw type system from a cell model with active motion[END_REF]. Thus we can extract a subsequence converging, in the weak sense of measures, towards a measure denoted µ in

M b (R d × R × [0, T ]). Because Q k (x, ξ, t) = W k -ξ + νG(ξ) is positive for ξ ≤ p m , we have µ(x, ξ, t) ≥ 0 for ξ ≤ p m .
Therefore passing to the limit k → +∞ into (38), in the sense of distributions,

∂ t χ -div[χ • ∇W ∞ ] + χ W ∞ -ξ ν - 1 ν ∞ ξ χ(x, η, t)dη = µ.
This last equation can also be written with (31)

∂ t χ -div[χ • ∇W ∞ ] + χ W ∞ -ξ ν -f (x, t) (H(W ∞ ) -ξ) + ν = µ,
and thus

∂ t χ -div[χ • ∇W ∞ ] + χ W ∞ -H(W ∞ ) ν = µ. ( 39
)
Using the assumption [START_REF] Elliotta | A variational inequality approach to Hele-Shaw flow with a moving boundary[END_REF], this equation is complemented with the initial condition

χ(x, ξ, t = 0) = 1 Ω 0 1 {0<ξ<H(W∞(x,t=0))} and f (x, t = 0) = 1 Ω 0 := f 0 (x).
It is useful to keep in mind the equivalent form of this equation,

∂ t χ -∇χ • ∇W ∞ + χ p ∞ -H(W ∞ ) ν = µ ≥ 0.
and thus, using (32),

∂ t χ -∇χ • ∇W ∞ + χ H(W ∞ ) f -1 ν = µ ≥ 0. ( 40
)
We can also integrate (39) and recover

∂ t p ∞ -div[p ∞ • ∇W ∞ ] + p ∞ ν [W ∞ -H(W ∞ )] = µdξ.
4th step. The set {g(x, t) = 1 and ξ < p m }. It is useful to consider the function

g(x, t) = f 0 X (x,t) (s = 0) ,
with the characteristics defined by

d ds X (x,t) (s) = -∇W ∞ (X (x,t) (s), s), X (x,t) (t) = x. (41)
This function g is the solution of the transport equation

∂ t g -∇g • ∇W ∞ = 0, g 0 = f 0 .
Using (40) and 0 ≤ f ≤ 1, we find

∂ t f -∇f • ∇W ∞ = µ(x, ξ, t) + χ H(W ∞ ) 1 -f ν ≥ 0. ( 42
)
From the comparison principle, we conclude that f (x, t) ≥ g(x, t) and we conclude that,

   f (x, t) = g(x, t) = 1, in the set {g(x, t) = 1}, µ(x, ξ, t) = 0 in the set {g(x, t) = 1 and ξ < p m }. ( 43
)
5th step. Strong convergence of p k . Another wording for step 4, is that

Ω(t) = Y (x) (t)[Ω 0 ] = {p ∞ (•, t) > 0}, with Y (x) (t) the limiting flow of Y (k) (x) (t
) defined in [START_REF] Preziosi | Multiphase modelling of tumour growth and extracellular matrix interaction: mathematical tools and applications[END_REF]. Indeed, from [START_REF] Ranft | Fluidization of tissues by cell division and apoptosis[END_REF] and the strong convergence of the flow, we infer that

p ∞ (•, t) = 0 in Y (x) (t)[R d \Ω 0 ].
Then we have f (x, t) = 1 Ω(t) = 1 {p∞(x,t)>0} . We recall that by definition, f = w-lim k→+∞ 1 {p k ≥pm/2} . We show that it implies the strong convergence locally in L 1 ((0, T )×R d ) of p k towards H(W ∞ )1 {p∞>0} . Let U be an open bounded subset of R d , we have

T 0 U |p k -H(W ∞ )1 {p∞>0} | dx ≤ I k + II k + III k , (44)
with

I k = T 0 U 1 {p k ≥pm/2} |p k -H(W ∞ )| dx, II k = T 0 U 1 {p k <pm/2} p k dx, III k = T 0 U H(W ∞ ) 1 {p k ≥pm/2} (1 -1 {p∞>0} ) + 1 {p k <pm/2} 1 {p∞>0} dx.
For the first term I k , we have that

I k ≤ T 0 U 1 {p k ≥pm/2} |p k -H(W k )| dx + T 0 U 1 {p k ≥pm/2} |H(W k ) -H(W ∞ )| dx ≤ 2 p m T 0 U p k |p k -H(W k )| dx + C T 0 U |W k -W ∞ | dx.
Using estimate [START_REF] Perthame | Derivation of a Hele-Shaw type system from a cell model with active motion[END_REF], we deduce that the first term of the right hand side goes to 0 as k → +∞. From the local strong convergence of W k towards W ∞ , the second term of the right hand side converges to 0 too. We conclude that lim k→+∞ I k = 0. Moreover, it has been proved in Lemma 2.2, see equation ( 27), that lim k→+∞ II k = 0. For the last term, we have, using the fact that W ∞ is bounded in L ∞ , that for some nonnegative constant C,

III k ≤ C U 1 {p k ≥pm/2} 1 -1 {p∞>0} + 1 -1 {p k ≥pm/2} 1 {p∞>0} dx
We have shown in the 4th step above that 1 {p k ≥pm/2} converges weakly towards 1 {p∞>0} . Then passing to the limit k → +∞ in the latter inequality, we deduce that lim k→+∞ III k = 0. We conclude from (44) that, for any open bounded subset U ,

T 0 U |p k -H(W ∞ )1 {p∞>0} | dx -→ k→+∞ 0.
By uniqueness of the weak limit, we deduce that p ∞ = H(W ∞ )1 {p∞>0} a.e.

6th step. Derivation of (29). From definition (35), this limit is now a consequence of

k T 0 R d p k (x, t)|Q k (x, t)| dxdt = T 0 R d (0,∞) |µ k (x, ξ, t)| dξdxdt.
But µ k vanishes for k → ∞ because from (42) we infer that µ = 0 both when f = 1 and f = 0. Therefore, we find (29).

Proof of Theorem 1.1

The proof of the Theorem 1.1 can now be easily deduced from Lemma 2.3. First, up to a subsequence, we have that p k converges a.e. towards p ∞ . On the one hand, recalling that the sequence (p k ) is uniformly bounded in L ∞ , we use the Lebesgue dominated convergence Theorem to show that, for any bounded open U ,

T 0 U p k |p k -W k -νG(p k )| dx -→ k→+∞ T 0 U p ∞ |p ∞ -W ∞ -νG(p ∞ )| dx.
On the other hand, we have from estimate (18) that

T 0 U p k |p k -W k -νG(p k )| dx -→ k→+∞ 0.
We deduce that p ∞ p ∞ -W ∞ -νG(p ∞ ) = 0 a.e. that is [START_REF] Perthame | Existence of solutions of the hyperbolic Keller-Segel model[END_REF]. We may apply the strong convergence for transport equations, as in [START_REF] Ben Belgacem | Compactness for nonlinear transport equations[END_REF][START_REF] Diperna | Ordinary differential equations, transport theory and Sobolev spaces[END_REF], to conclude that, since the term G(p k ) converges strongly, n k , which solves the transport equation (1), itself converges strongly. Note in particular that, from assumption (13), we have n 0 k -→ k→+∞ 1 Ω 0 . Passing to the limit in the equation (1), we recover the limit equation for n ∞ [START_REF] Escher | Classical solutions for Hele-Shaw models with surface tension[END_REF]. Finally, passing to the limit in the relation, 16) is then a direct consequence of Lemma 2.3.

n k p k = k k -1 1-1 k-1 p k/(k-1) k , we deduce that (1 -n ∞ )p ∞ = 0. The relation (

One dimensional traveling waves

In order to examplify Theorem 1.1 and to give a simple case, with a solution that can be build analytically, we look for a one dimensional traveling wave solution to the Hele-Shaw limit.

Because, traveling waves are defined up to a translation, we may set, in the moving frame, Ω(t) = R + . Then, the system rewrites p = 0, -σn ′ -(nW ′ ) ′ = nG(0), -νW ′′ + W = 0, for x > 0, (45)

n = 1, -νW ′′ + W -H(W ) = 0, p = H(W ), for x < 0. ( 46
)
Moreover, the jump condition at the interface x = 0 implies -σ[n] -[nW ′ ] = 0, which leads to the traveling velocity σ = -W ′ (0).

We denote W 0 := W (0). For x > 0, we have

W (x) = W 0 e -x/ √ ν , (47)
from which we deduce that σ = W 0 √ ν .

Then we can rewrite the first equation in (45) as

-n ′ (x) W 0 √ ν - W 0 √ ν e -x/ √ ν = n(x) G(0) + W 0 ν e -x/ √ ν .
Taking the limit x → 0 leads to n(0) = 0. Moreover, since n ′ ≤ 0, we deduce that n = 0 on (0, +∞). For x < 0, we solve the second order ODE for W with boundary condition W (0) = W 0 and W ′ (0) = -W 0 / √ ν. As an example, we choose for the growth term the function G(p) = P M -p, and thus

H(W ) = W + νP M 1 + ν . (48)
Then equation (46) for W rewrites :

-(ν + 1)W ′′ + W = P M .
The only bounded solution on (-∞, 0) such that W (0) = W 0 is given by

W (x) = P M + (W 0 -P M )e x/ √ ν+1 .
Moreover, the continuity of the derivative implies, from (47), that W ′ (0) = -W 0 / √ ν. We deduce the value for W 0 :

W 0 = √ ν √ ν + √ ν + 1 P M .
Then we conclude that for x < 0,

W (x) = P M 1 - 1 1 + ν/(ν + 1) e x/ √ ν+1 .
The pressure is then given by :

p(x) = P M 1 - 1 ν + 1 + ν(ν + 1) e x/ √ ν+1 .
and the traveling velocity

σ = P M √ ν + √ ν + 1 .
We notice that the pressure is nonnegative and has a jump at the interface x = 0. The height of the jump is given by P M 1 -1/(ν + 1 + ν(ν + 1)) . We observe moreover that σ is a decreasing function of ν. Letting ν → 0, we recover the result for the Hele-Shaw model for purely elastic tumors [START_REF] Perthame | Traveling Wave Solution of the Hele-Shaw Model of Tumor Growth with Nutrient[END_REF][START_REF] Tang | Composite waves for a cell population system modelling tumor invasion[END_REF].

Figure 1: Plot of the density n (dashed line), pressure p (line with dot), W (continuous line). Left : for ν = 1 and at final time T = 25 s. We notice a jump for the density from 0 to 1 at the front and a jump of the pressure. Right : for ν = 0 and at final time T = 12.5 s. In this case, we have p = W and there is no jump on the pressure; moreover, the velocity of the front is faster. This observation is compatible with the interpretation that viscosity acts as a friction.

Numerical simulations. Finally, we present numerical simulations of the system (1)-(3) in one dimension. We use a discretization thanks to a cartesian grid of a bounded domain [-L, L] of the real line. Equation ( 1) is discretized by a finite volume upwind scheme. Equation ( 2) is discretized thanks to finite difference scheme. Since we focus on the case where k is large, we use k = 100 in the numerical computation. For the initial data, we choose n 0 = 1 [-0.2,0.2] . The growth function G is chosen as in (48) with P M = 1.

In Figure 1, we display the shape of the density n, the pressure p and W obtained by the numerical simulation. The figure on the left displays the result with a viscosity coefficient ν = 1. For the comparison, we plot on the right of Figure 1, the shape in the case without viscosity (ν = 0). Comparing both figures, we observe that in the case ν = 1, we have a jump of the pressure at the interface of the solid tumor, whereas in the case ν = 0, the pressure is continuous at the interface.

We display in Figure 2 the first steps of the formation of the propagating front with the initial data n 0 = 1 [-0.2,0.2] . For this simulation we take ν = 1 and k = 100. The dynamics is represented thanks to the plot at 4 successives times of the density n, pressure p and W . After a transitory regime during which the pressure increases until reaching its maximal value P M = 1, the shape of the traveling waves is obtained and the front of the tumor invades the healthy tissue.

Conclusion

A geometric model, also called incompressible, has been derived from a cell density model (also called compressible) when the pressure law is stiff. Because the viscosity is considered here, the limiting problem is a free boundary problem for the set Ω(t) of non-zero pressure. The limiting system for the pressure consists in an algebraic relation between the pressure p ∞ and the limiting potential W ∞ (17), coupled with an elliptic equation for the potential W ∞ set in the whole space [START_REF] Friedman | A hierarchy of cancer models and their mathematical challenges[END_REF]. This is a major difference with the case where viscosity is neglected, the so-called Hele-Shaw system [START_REF] Friedman | A hierarchy of cancer models and their mathematical challenges[END_REF][START_REF] Lowengrub | Nonlinear modelling of cancer: bridging the gap between cells and tumours[END_REF]; then, the pressure is given by an elliptic equation for the pressure in the moving domain Ω(t). A paradox is that the effect of keeping viscosity generates a jump of the pressure at the interface of the region defining the tumor, unlike in [START_REF] Perthame | The Hele-Shaw asymptotics for mechanical models of tumor growth[END_REF][START_REF] Perthame | Derivation of a Hele-Shaw type system from a cell model with active motion[END_REF] where the Hele-Shaw problem is complemented with Dirichlet boundary conditions and therefore the pressure is continuous. This point is also observed in the numerical simulations in Section 3. The velocity of the propagating front of the tumor is given by the equation satisfied by the density [START_REF] Escher | Classical solutions for Hele-Shaw models with surface tension[END_REF]. Because the pressure is discontinuous, it has weaker regularity that in the inviscid case treated in [START_REF] Perthame | The Hele-Shaw asymptotics for mechanical models of tumor growth[END_REF][START_REF] Perthame | Derivation of a Hele-Shaw type system from a cell model with active motion[END_REF] and we need to develop a new strategy of proof to derive the incompressible limit. Our approach is based on a kinetic formulation of the equation satisfied by the pressure. This work also opens several additional questions. First, the case of general initial data is not treated here because we assume that n 0 vanishes outside Ω 0 . Then, it would be interesting to consider the case with active motion as in [START_REF] Perthame | Derivation of a Hele-Shaw type system from a cell model with active motion[END_REF]. In such a case, equation ( 1) is replaced by a parabolic equation. Then the structure of the problem is different but the limiting system should be the same, except the equation for the density which implies then a faster propagation of the region Ω(t). Finally, it is formally clear from ( 10)-( 12) that letting ν → 0, we recover the Hele-Shaw system. However, a rigorous proof of this fact requires compactness of the sequence which we is not directly available with the method developed here.

Figure 2 :

 2 First steps of the formation of the propagating front with k = 100 and ν = 1. The density n (dashed line), the pressure p (line with dot) and the potential W (continuous line) are represented at 4 successives times : a) t = 0.1 s, b) t = 1.25 s, c) t = 3.75 s and d) t = 12.5 s.

Acknowledgment. This work has been supported by the french "ANR blanche" project Kibord: ANR-13-BS01-0004.


Fengwei Chen

Hugues Garnier

Marion Gilson

Parameter and Time-delay Identification of Continuous-time Models from Non-uniformly Sampled Data

This paper considers the problem of continuous-time model identification with arbitrary time-delay from irregularly sampled data. The proposed method estimates the plant and the time-delay in a separable way, when estimating one of them, the other is assumed to be fixed. More precisely, the plant is estimated by the iterative instrumental variable SRIVC method while the time-delay is estimated by the Gauss-Newton method. Because of the nonlinear relationship between the loss function and the time-delay, a low-pass filter is employed to extend the global convergence region for the time-delay estimation. Numerical examples are presented to illustrate the properties of the proposed method.

I. INTRODUCTION

Irregularly sampled data occurs quite frequently in many application areas. For example, in astronomy, particularly in stellar physics, the obtained data always has unknown periodicity. In even-based sampling (e.g. Lebesgue sampler), the sampler is triggered by the arrival of a certain event process. When one wants to model a system from irregularly sampled data, the use a continuous-time (CT) model is preferable, since the CT model parameters are independent of the sampling period (see e.g. [START_REF]Identification of Continuous-time Models from Sampled Data[END_REF], [START_REF] Chen | Identification of continuous-time transfer function models from nonuniformly sampled data in presence of colored noise[END_REF]). Many processes have intrinsic delay, as it is the case in thermal and chemical processes, therefore it is a common practice to assume a time-delay in the model. The problem of time-delay estimation based on a discrete-time (DT) model has been studied a lot. Usually the time-delay is assumed to be an integer number of the sampling period, thus the delayed output of the system can be obtained by adjusting the phase-shift. However, this flexibility is lost in the situation of non-uniform sampling because of the time-varying sampling intervals. To handle irregularly sampled data, the algorithm must be capable of estimating arbitrary time-delay. Several methods have been proposed for CT model parameter and time-delay estimation. A linear filter-based method was introduced in [START_REF] Ahmed | Parameter and delay estimation of continuous-time models using a linear filter[END_REF], where the time-delay along with the transfer function parameters are estimated in an iterative way through simple linear regression. Iterative global nonlinear least-squares and instrumental variable methods were suggested in [START_REF] Yang | Identification of continuous-time systems with multiple unknown time delays by global nonlinear least-squares and instrumental variable methods[END_REF], where the plant parameters and time-delay are estimated in a separable way. More recently, a wavelet correlation method was introduced in [START_REF] Ni | Time delay estimation for MIMO dynamical systems -With time -frequency domain analysis[END_REF] to estimate the time-delay for MIMO dynamical systems, by calculating and Fengwei Chen, Hugues Garnier and Marion Gilson are with the Centre de Recherche en Automatique de Nancy, Université de Lorraine, CNRS, 2 rue Jean Lamour, 54519 Vandoeuvre-lès-Nancy, France. Email: fengwei.chen, hugues.garnier, marion.gilson@univ-lorraine.fr handling the cross-correlation between the CWT coefficients of the input and the output data. In this paper, the time-delay and the plant parameters are also estimated in a separable way, but a low-pass filter is introduced to improve the convergence performance for the time-delay estimation (see [START_REF] Ferretti | On the identifiability of the time delay with least-squares methods[END_REF]). More precisely, the plant and the time-delay are estimated by the iterative instrumental variable (IV) and the Gauss-Newton methods, respectively. For the IV method, the simplified refined instrumental variable method for CT systems (SRIVC) is adopted. It was first developed by Young and Jakeman in the open-loop situation [START_REF] Young | Refined instrumental variable methods of time-series analysis: Part III, extensions[END_REF] and recent developments can be found in [START_REF] Garnier | An optimal IV technique for identifying continuous-time transfer function model of multiple input systems[END_REF], [START_REF] Young | Refined instrumental variable identification of continuous-time hybrid Box-Jenkins models[END_REF]. This paper is organized in the following way. We first define the parameter estimation problem in Section II. Subsequently, the SRIVC method for plant parameter identification is recalled in Section III. Then the time-delay estimation problem is discussed in Section IV. The final algorithm for time-delay system identification is detailed in Section V. Finally, in Section VI, numerical examples are presented to illustrate the properties of the proposed method.

II. PROBLEM STATEMENT

Consider the following CT output error (COE) model

   x(t) = G(p, θ o )u(t -τ o ) = B(p, θ o ) A(p, θ o ) u(t -τ o ) y(t k ) = x(t k ) + e(t k ) (1)
where u(t), x(t) and τ o are the excitation signal, the noisefree response and the time-delay, respectively. p denotes the differentiation operator, i.e. px = dx/dt. B(p, θ o ) and A(p, θ o ) are polynomials assumed to be coprime of the following form

B(p, θ o ) = b o 0 p m + b o 1 p m-1 + • • • + b o m A(p, θ o ) = a o 0 p n + a o 1 p n-1 + • • • + 1 (n > m)
The unknown parameters are stacked columnwise in

θ o = [a o 0 • • • a o n-1 b o 0 • • • b o m ] T (2)
The input-output signals of the system are observed at irregular time-instant t k , for k = 1, 2, • • • , N , the sampling interval is denoted as

h k = t k+1 -t k (3)
In most practical measurement situations, it is reasonable to consider that x(t k ) is corrupted by a DT measurement noise e(t k ), which is assumed to be white here. The identification problem can be stated as: assume the orders n and m are known, the identification objective is to estimate the unknown parameters and time-delay from irregularly sampled data Z N = {u(t k ); y(t k )} N k=1 .

III. INSTRUMENTAL VARIABLE METHOD FOR PLANT PARAMETER ESTIMATION

Amongst the different implementation of the optimal IV approach, the SRIVC method is known to be one of the most effective and is used here to identify the CT transfer function parameters. In this section, the time-delay is assumed to be known (τ = τ o ). The SRIVC method makes use of a low-pass filter to generate the (filtered) input/output timederivatives of the differential equation model. The filter F (p) takes the following form

F (p, θj ) = 1 A(p, θj ) (4)
Then, the COE model can be equivalently reformulated in the following regression form

y F (t k , θj ) = φ T F (t k , θj , τ o ) θj + ε(t k , θj , τ o ) (5)
with

φ T F (t k , θj , τ o ) = -y (n) F (t k , θj ) • • • -y (1)
F (t k , θj ), u (m) F (t k -τ o , θj ) • • • u F (t k -τ o , θj ) (6)
where ε(t k , θj , τ o ) is the equation error, (•) F denotes the filtering operation, i.e. (•) F = F (p, θj )(•). The IV vector ψ T F (t k , θj , τ o is chosen as

ψ T F (t k , θj , τ o ) = -x (n) F (t k , θj , τ o ) • • • -x (1)
F (t k , θj , τ o ) u (m) F (t k -τ o , θj ) • • • u F (t k -τ o , θj ) (7)
where x(t k , θj , τ o ) is the noise-free response of the following auxiliary model

x(t k , θj , τ o ) = B(p, θj ) A(p, θj ) u(t k -τ o ) (8)
From N measured input-output data, the SRIVC parameter estimates are given by

θj+1 = Ψ F ( θj , τ o )Φ T F ( θj , τ o ) -1 Ψ F ( θj , τ o )y F ( θj ) (9) with y F ( θj ) = y F (t s , θj ) • • • y F (t N , θj ) T Φ F ( θj , τ o ) = φ F (t s , θj , τ o ) • • • φ(t N , θj , τ o ) Ψ F ( θj , τ o ) = ψ F (t s , θj , τ o ) • • • ψ F (t N , θj , τ o )
where t s ≥ τ o . As it can be noted from ( 9), the SRIVC method uses an adaptive procedure where the parameter estimates are iteratively refined. At each iteration, an auxiliary model is used to generate the instrumental variables and filter based on the parameters obtained at the previous iteration. When performing the CT filtering operation, the intersample behavior is required to reconstruct a 'CT input signal' from the irregularly sampled data. In many situations, as it is the case in computer controlled systems, the intersample behavior of the input u(t k ) is always known, such as piece-wise constant (zero-order hold) or piece-wise linear (first-orderhold). However, the intersample behavior of the sampled output y(t k ) is missing. It is then a common practice to assume y(t) is first-order hold. The 4th-order Runge-Kutta (RK4) method can then be used to compute the filtered timederivatives for the output signal at the required irregular timeinstants.

IV. EXPANSION OF GLOBAL CONVERGENCE REGION IN TIME-DELAY ESTIMATION

In this section, the plant-independent time-delay estimation is investigated, so θ = θ o is assumed in the sequel. τ is assumed to be a 'pure delay' and appears as an explicit parameter in the CT model. With these assumptions, the estimated error can be formulated as

ǫ(t, θ o , τ ) = y(t) -G (p, θ o ) u(t -τ ) (10)
The time-delay can be estimated by minimizing a certain criterion of ǫ(t, θ o , τ ). Here the energy of the estimated error is considered. The estimate of τ can be obtained by the following fitting

τ = arg min τ J(θ o , τ ) = arg min τ 1 2 ∞ -∞ ǫ 2 (t, θ o , τ )dt (11)
Numerical methods are typical ways to solve [START_REF] Young | Refined instrumental variable identification of continuous-time hybrid Box-Jenkins models[END_REF], such as the Gauss-Newton or the Levenberg-Marquardt algorithm.

Since the cost function J(θ o , τ ) is multi-modal with respect to τ , the numerical method suffers from local minima and the final estimate τ is highly dependent on the initial value τ 0 .

To increase the chance of converging to the global minimum, the use of a low-pass filter was suggested in [START_REF] Ferretti | On the identifiability of the time delay with least-squares methods[END_REF]. This result is used in this paper to improve the convergence performance of the proposed method. Let L(p) be a CT low-pass filter with the cutoff frequency ω LF c . The estimate τ can be obtained by minimizing the energy of the filtered estimated error

τ = arg min τ J(θ o , τ ) (12)
J(θ o , τ ) = 1 2 ∞ -∞ [L(p)ǫ(t, θ o , τ )] 2 dt (13)
From the Parseval theorem, J(θ o , τ ) also has the following frequency-domain interpretation (see e.g. [START_REF] Ljung | System Identification -Theory for the User[END_REF], [START_REF] Ferretti | On the identifiability of the time delay with least-squares methods[END_REF])

J(θ o , τ ) = 1 4π +∞ -∞ G(iω, θ o )e -iωτ o -G(iω, θ o )e -iωτ 2 × Φ u (ω) + Φ e (ω) |L(iω)| 2 dω ( 14
)
where Φ u (ω) and Φ e (ω) are the power spectral densities (PSD) of u(t) and e(t), respectively. Let the time-delay error denoted as

δτ = τ -τ o (15)
Let us now define V (θ o , τ ) as the noise-free version of ( 14)

V (θ o , δτ ) = 1 4π +∞ -∞ G(iω, θ o )e -iωτ o 2 1 -e -iωδτ 2 × Φ u (ω) |L(iω)| 2 dω (16)
From ( 16), it can be observed that V (θ o , δτ ) depends on δτ and the following transfer function

T (iω) = L(iω)G(iω, θ o )e -iωτ o (17)
Since T (iω) requires the knowledge of the true model G(iω, θ o ), which is not known in practical situations, the following theorem provides an alternative to compute V (θ o , δτ ) from filtered data. Theorem 1 (see [START_REF] Ferretti | On the identifiability of the time delay with least-squares methods[END_REF]): Assume the PSD of u(t) is denoted by Φ u (ω) and let z(t) be the noise-free output of the following filter

T (iω) = L(iω)G(iω, θ o )e -iωτ o
Then an estimate of V (θ o , δτ ) can be obtained by

V (θ o , δτ ) = R z (0) -R z (δτ ) (18)
where R z is the autocorrelation function of z(t).

Proof: V (θ o , δτ ) can be given as

V (θ o ,δτ ) = 1 4π +∞ -∞ |T (iω)| 2 1 -e -iωδτ 2 Φ u (ω)dω = 1 4π +∞ -∞ Φ z (ω) 2 -e -iωδτ -e iωδτ dω (19)
where

Φ z (ω) = |T (iω)| 2 Φ u (ω)
is the PSD of z(t). By using the well-known Wiener-Khinchin theorem, (19) can be reformulated as

V (θ o , δτ ) =R z (0) - 1 2 R z (-δτ ) - 1 2 R z (δτ ) = R z (0) -R z (δτ ) (20)
The second equation in (20) uses the fact that R z (δτ ) is an even function.

To illustrate the low-pass filtering effects on the global convergence, let us consider the following example.

Example 1: Consider the following system

x(t) = b o 0 a o 0 p + 1 u(t -τ o ) (21)
a) The system (with a o 0 = 0.5, b o 0 = 0.5 and τ o = 5) is observed for N = 20, 000 time-instants, the sampling period h k is uniformly distributed in [0.01, 0.09]s. b) To approximate Φ u (ω) = 1, u(t k ) is chosen to be a zero-mean random sequence with the variance of

E u 2 (t k ) = 1/h k c) Three different low-pass filters are considered 1) Filter 1: L 1 (iω) = 1 (No filtering) 2) Filter 2: L 2 (iω) = 1 |ω| ≤ ω LF c 0 |ω| > ω LF c
3) Filter 3:

L 3 (iω) = G -1 (iω, θ o )e iωτ o |ω| ≤ ω LF c 0 |ω| > ω LF c
Actually the cut-off frequency ω LF c of L(iω) can take any positive value, given the fact the purpose of the low-pass filtering is to extend the global convergence region. It is then advised to choose ω LF c to be lower than the system bandwidth ω b . In this example ω LF c = 0.25rad/s while ω b = 0.5rad/s.

Note that when filter L 3 (iω) is chosen, T (iω) becomes an ideal low-pass filter. V (θ o , δτ ) has then the following timedomain representation (see [START_REF] Ferretti | On the identifiability of the time delay with least-squares methods[END_REF])

V (θ o , δτ ) = 1 4π ω LF c -ω LF c 1 -e -iωδτ 2 dω = ω LF c π - sin(ω LF c δτ ) πδτ (22)
Let δτ max denote the maximum delay uncertainty (or admissible region uncertainty), for which |δτ | ≤ δτ max , such that V (θ o , δτ ) is unimodal. Equation ( 22) can then be transformed to make a cardinal sine function appear

r = sin(t)/t (23)
A cardinal sine function has many extrema. An approximation of the t-coordinate (t ≥ 0) of the i-th extremum is

t i = (i + 1/2)π - 1 (i + 1/2)π (24)
Because t 0 = 0 is the global minimum, t 1 can then be used to compute the maximum delay uncertainty δτ max . By equating δτ L 1 max =1. 9 δτ L 2 max =8. 6

t 1 to ω LF c δτ L3 max , we have δτ L3 max ≈ t 1 /ω LF c ≈ 4.5/ω LF c (25)
δτ L 3 max =8. 5 δτ V ( θ o , δτ) /Rz(0) FilterL1 FilterL2 FilterL3
Fig. 1. Normalized cost function V (θ o , δτ )/Rz(0) for the 3 considered filters.

The normalized cost function V (θ o , δτ )/R z (0) versus δτ is plotted in Figure 1. It can be observed that without any low-pass filtering, δτ L1 max ≈ 1.9. When a low-pass filtering is applied, the admissible region is extended to δτ L2 max ≈ 8.6 and δτ L3 max ≈ 8.5. The global convergence region is clearly enlarged by using a low-pass filter.

From Figure 1, it can also be noticed that δτ L2 max ≈ δτ L3 max ≈ 4.5/ω LF c . This means that L 2 (iω) and L 3 (iω) have almost the same effect. This phenomenon can be explained as follows: when ω LF c is much smaller than the system bandwidth, if ω < ω LF c , G(iω, θ o ) ≈ constant, its corresponding T (iω) can also be regarded as an ideal low pass filter. In the remainder of the paper, the ideal low-pass filter is used

L(iω) = 1 |ω| ≤ ω LF c 0 |ω| > ω LF c ( 26
)
Remark. 1: If the data are regularly sampled, the filtering operation can be carried out in the frequency domain by using the Fourier transform and its inverse, by adding a rectangular window to the frequency data. In the present case of irregularly sampled data, the non-uniform fast Fourier transform (NUFFT) and its inverse can be used (see e.g. [START_REF] Ware | Fast approximate Fourier transforms for irregularly spaced data[END_REF]).

V. SEPARABLE METHOD FOR TIME-DELAY SYSTEM IDENTIFICATION

From the previous statements, θ is estimated by the SRIVC method while τ is estimated by a numerical method. So it is natural to use a separable way to estimate θ and τ , when estimating each of them, the other is assumed to be fixed. The regression model once the additional filter L(p) is applied, takes the following form

ȳF (ρ) = ΦT F (ρ, τ )θ + ǭ(ρ) (27)
where ρ = θ T , τ T , and x means x is filtered by L(p). The cost function to be minimized is given as

JN (ρ) = 1 2(N -s + 1) ǭT (ρ)ǭ(ρ) ( 28
)
where s is chosen to guarantee that t s ≥ τ . The parameters can be estimated by the following separable method (see [START_REF] Ngia | Separable nonlinear least-squares methods for efficient off-line and on-line modeling of systems using Kautz and Laguerre filters[END_REF]). Compared with the original separable non-linear leastsquares method, we propose to estimate the linear part by the SRIVC method instead of the simple least-squares method, and perform a low-pass filtering operation in the first few iterations of the algorithm to widen the convergence region for the initial time-delay search. Theorem 2: Let

θ(τ ) = Ψ F (ρ) ΦT F (ρ) -1 Ψ F (ρ) ȳF (θ) (29)
Then τ can be obtained as

arg min ρ JN (ρ) = arg min τ JN (τ ) (30) τ = arg min τ JN (τ ) (31)
where

JN (τ ) = JN (ρ) θ= θ(τ )
The proposed Gauss-Newton method is then given as

τ j+1 =τ j -µ j ∇ 2 JN (τ j ) -1 ∇ JN τ j (32) θj+1 = Ψ F ( θj , τ j+1 ) ΦT F ( θj , τ j+1 ) -1 × Ψ F ( θj , τ j+1 ) ȳF ( θj ) (33)
where

∇ JN (τ ) = 1 N -s + 1 ǭT τ ǭ(ρ) ∇ 2 JN (τ ) ≈ 1 N -s + 1 ǭT τ ǭτ -ǭT τ ǭθ ǭT θ ǭθ -1 ǭT θ ǭτ ǫ τ = ∂ǫ(ρ) ∂τ ρ= ρj = pG(p, θj )u(τ j ) ǫ θ = ∂ǫ(ρ) ∂θ ρ= ρj = -Ψ T F ( ρj )

A. The Proposed SRIVC-based Algorithm

The complete algorithm (tdsrivc) for time-delay system identification can be summarized as 1) Initialization: Set boundaries ∆τ min , ∆τ max , τ min , τ max , the cut-off frequencies ω SVF c , ω LF c , and ∆ JN (τ ) min . Initialize F (p) by the following state-variable-filter (SVF)

F (p) = 1 (p + ω SVF c
) n Based on the initial value τ 1 and F (p), use srivc (see e.g. [START_REF] Young | Refined instrumental variable identification of continuous-time hybrid Box-Jenkins models[END_REF] for more details about srivc method and the choice of ω SVF c ) algorithm to compute θ1 . 2) Initialization of the time-delay search for j=1:M a) Compute ∆τ j

∆τ j = -∇ 2 JN τ j -1 ∇ JN τ j (34) b) Perform the following i) Compute τ j+1 = τ j + ∆τ j . If τ j+1 / ∈ [τ min , τ max ], let ∆τ j = ∆τ j /2 and repeat this step. ii) If |∆τ j | > ∆τ max , let |∆τ j | = ∆τ max .
If |∆τ j | < ∆τ min , go to step (c). iii) Estimate θ using srivc Update auxiliary model and prefilter

G(p, θj ) = B p, θj / Â p, θj F (p, θj ) = 1/ Â p, θj
Compute the filtered signals in iv) Compute JN τ j+1 . If JN τ j+1 ≥ JN (τ j ), let ∆τ j = ∆τ j /2 and go to step iii. c) Check the stop condition, if JN (τ j ) -JN τ j+1 ≥ ∆ JN (τ ) min , go to step (a), else break. end 3) Refined plant parameter and delay estimation for j=M :convergence Repeat step 2) without applying the low-pass filter L(iω). end Note that in the above algorithm, tdsrivc runs in the end without low-pass filtering to obtain a refined estimate. The reason for using a two-step scheme is that the lowpass filtering increases the convergence region at the price of slightly decreasing the accuracy of the final estimate, so a refinement step is suggested. In addition, the maximum number of iterations M for the initial time-delay search can be typically chosen between 5 to 10.

Ψ F θj , τ j+1 , ΦF θj , τ j+1 , ȳF θj Then estimate θ θj+1 = Ψ F θj , τ j+1 ΦT F θj , τ j+1 -1 × Ψ F θj , τ j+1 ȳF θj

B. The Possibility of Handling Arbitrary Time-delay

As discussed in Section III, when performing CT filtering, the actual input of the filter is a 'CT signal', which is reconstructed from the irregularly sampled data and the provided inter-sample behavior. With this reconstructed CT signal, it is possible to predict the output of a filter at any desired time-instant with the specified time-delay. In this way, the proposed method can therefore handle arbitrary time-delay identification.

VI. NUMERICAL EXAMPLES In this section, all the systems to be identified are assumed to be excited by a pseudo random binary sequence (PRBS), which is generated from a 9-stage shift register and the clock is set to 0.5s. The time-varying sampling period h k is assumed to be uniformly distributed in the following interval

h k ∼ U [0.01, 0.09]s
The stopping rule for the iterative algorithm is that the relative change in JN (τ ) or JN (τ ) is lower than 10 -4 ; the maximum iteration number of stage 2) is M = 10; the timedelay boundaries are τ min = 0 and τ max = 10.

A. Noise-free Case

To illustrate the positive effect of the low-pass filtering strategy on the convergence performance, let us first consider the noise-free situation. The second-order system to be studied is given as

y(t) = 2 0.25p 2 + 0.7p + 1 u(t -5) (35)
Different initializations for the cut-off frequency of the initial filter ω SVF c used in the SRIVC method and τ 0 are tested one by one

ω SVF c = {0.2π, 0.4π, 0.6π, • • • , 3π} rad/s τ 0 = {0, 1, 2, • • • , 10}s
for two cut-off frequency of the low-pass filter (ω LF c = 0.4π, ω LF c = 0.2π). For each run, the number of observations is N = 10, 000. The criterion of global convergence is fit > 99, where fit is the normalized root-mean-square error defined as follows fit = 100 1 -y(t) -x(t) 2 y(t) -E {y(t)} 2 where x(t k ) is the noise-free response of the estimated model, y(t k ) is the sampled output. Two identification schemes has been tested:

A. Identification without low-pass filtering. B. Identification with low-pass filtering. The estimated results are presented in Fig. 2 for the two chosen values of ω LF c . From this figure we see that the results for scheme A (no filtering) are very poor. Convergence is reached when τ 0 lies in the interval [START_REF] Garnier | An optimal IV technique for identifying continuous-time transfer function model of multiple input systems[END_REF][START_REF] Ljung | System Identification -Theory for the User[END_REF]. Meanwhile, a SVF filter with lower cut-off frequency slightly increases the probability of global convergence. The results for Scheme B are overall quite satisfying. It can be noted that when ω LF c decreases, the global minimum is more easy to access. One question remains here to know how ω L should be chosen. It has concluded that the global minimum is easier to access for smaller values of ω LF c . Then, one may ask, can ω LF c be arbitrarily small? The answer is no, the reason can be given as: smaller values for ω LF c mean that more information of the sampled data is removed, which can deteriorate the plant parameter estimation. A rule of thumb can be to choose the cut-off frequency of the ideal low-pass filter as 1 10 to 1 2 of the system bandwidth.

B. Noisy-output Case

In this subsection, it is assumed that the measured output is corrupted by a discrete-time white noise, the signal to noise ratio (SNR) is 15dB. The SNR is defined as SNR = 10log P x P e

where P x and P e represent the average power of the noisefree signal x(t k ) and the disturbance e(t k ). The number of observed data is N = 2000. The cut-off frequency of the initial SVF filter is fixed to ω SVF c = 2π. Due to the fact that the time-delay of the system to be studied is lower than 10, any choice of ω LF c leads to δτ max ≥ 5 is reasonable, here we choose ω LF c = 0.1π. The initial guess of the time-delay is assumed to be uniformly distributed in the following interval τ 0 ∼ U [0, 10]s In this noisy-output case, the criterion for global convergence fit satisfies the following inequality fit ≤ 100 1 -√ P e √ P x + P e = 100 1 -1 + 10 SNR/10 -1/2 ≈ 82.5

The criterion for 'success' is fit > 81. The estimated results are given in Table I. This table shows that when no filter is applied (scheme A), the global convergence ratio is less than 60%, while the success ratio surges to 84% ∼ 100% when a filtering option (scheme B) is adopted. The good performance does not come for free, the cost is the slight increase of the iteration number to converge.

VII. CONCLUSION

In this paper, identification of continuous-time system with arbitrary time-delay from irregularly sampled data has been considered. The proposed algorithm estimate the plant parameters and the time-delay in a separable way, more precisely, the plant is estimated by the SRIVC method and the time-delay is estimated by the Gauss-Newton method, respectively. An ideal low-pass filter is introduced to extend the global convergence region. Numerical examples have been used to show that, under mild initializations, the global convergence ratio can be increased to 80%∼100% after the use of a low-pass filter.

Fig. 2 .

 2 Fig. 2. Estimated results with and without low-pass filtering. Ring -fit>99. Dot -fit≤99.

TABLE I ESTIMATED

 I MODELS FROM IRREGULARLY SAMPLED DATA. N iter -THE NUMBER OF ITERATIONS, P -THE RATIO OF GLOBAL CONVERGENCE. THE MEAN VALUE, THE STANDARD DEVIATION AND N iter ARE COMPUTED FROM THE SUCCESSFUL ESTIMATED MODELS.