J Waeytens
email: julien.waeytens@ifsttar.fr

B Rosić

P.-E Charbonnel

E Merliot

D Siegert

X Chapeleau

R Vidal

V Le Corvec

L.-M Cottineau

Model updating techniques for damage detection in concrete beam using optical fiber strain measurement device

Keywords: Structural Health Monitoring, damage detection techniques, model updating, optic fiber

The article aims at detecting and quantifying early structural damages using deterministic and probabilistic model updating techniques. To achieve this purpose, local information in a form of optical strain measurement is employed. The strategy consists in updating physical parameters associated to damages, such as Young's modulus, in order to minimize the gap between the numerical strain obtained from finite element solves and the strain sensor outputs. Generally, the damage estimation is an ill-posed inverse problem, and hence requires regularization. Herein, three model updating techniques are considered involving different type of regularization: classical Tikhonov regularization, Constitutive Relation Error based updating method and Bayesian approach.

An illustration of these three approaches is proposed for localizing and quantifying an early damage in a real 8 meter post-tensioned concrete beam. Numerical results show that all the methods properly localize the damaged area

Introduction

Structural Health Monitoring (SHM) has become of crucial importance in the civil engineering community for detecting and localizing damages, reducing the maintenance cost of structures and ensuring the safety of individuals. Three levels of damage assessment [START_REF] Rytter | Vibrational based inspection of civil engineering structures[END_REF] can be distinguished: damage detection (level 1), damage localization (level 2) and damage quantification (level 3). Levels 1 and 2 mostly rely on data-driven methods for computing damage indicators from measurements performed on vibrating structures. The damage indicators are generally designed to be sensitive to the variations of the modal signature of the structure under study and at least require a comparison with a reference state; different methods coexist and operate in time or frequency frameworks with deterministic or stochastic inputs. Basics and classical references on SHM can be found in [START_REF] Worden | The fundamental axioms of structural health monitoring[END_REF][START_REF] Brownjohn | Structural health monitoring of civil infrastructure[END_REF], whereas more recent works on vibration-based damage detection and localization are presented in [START_REF]New Trends in Vibration Based Structural Health Monitoring[END_REF][START_REF] Niezrecki | Structural Health Monitoring and Damage Detection[END_REF] and the references therein. Model updating in turn or model calibration techniques directly allows damage assessment up to level 3. The starting point of these methods is to consider a parametric model (generally derived by using a Finite Element framework) as sufficiently accurate to stand for the behavior of the real structure in an undamaged configuration; this reference model is then compared to the measurements performed on the actual damaged structure, and optimal parameters sets are sought as those ensure the best compatibility between model and measurements. The basic principle behind this consists in assuming that localized structural damage results in a local diminution of stiffness. Over the last few decades, uncertainty quantification or inverse problem resolution has gained substantially in interest and has led to a huge amount of approaches and strategies in various scientific domains. Among the most popular approaches for uncertainty quantification, one can non-exhaustively cite : the deterministic parametric residual minimization (deterministic approach, see e.g. [START_REF] Tarantola | Inverse problem theory and methods for model parameter estimation[END_REF] and the references therein, or [START_REF] Mottershead | Model updating in structural dynamics: a survey[END_REF][START_REF] Bonnet | Inverse problems in elasticity[END_REF] in linear mechanics), the probabilistic Bayesian approach (see e.g. [START_REF] Beck | Updating models and their uncertainties. i: Bayesian statistical framework[END_REF] or [START_REF] Litvinenko | Inverse problems and uncertainty quantification[END_REF]), and a non-probabilistic Fuzzy approach (see references in [START_REF] Zimmerman | Fuzzy set theory and its applications[END_REF][START_REF] Beer | Uncertain structural design based on nonlinear structural analysis[END_REF][START_REF] Degrauwe | Uncertainty propagation in structural analysis by fuzzy numbers[END_REF]). Bayesian approach was previously considered in a two-stage damage detection process consisting of Bayesian and immune genetic algorithms to determine the reduction of Young's moduli [START_REF] Guo | Structural damage identification based on bayesian theory and improved immune genetic algorithm[END_REF]. Similarly, in [START_REF] Papadimitriou | Bayesian methodology for structural damage identification and reliability assessmen[END_REF][START_REF] Jerez | A comparative study of sensor placement techniques for structural damage detection[END_REF] is considered the structural model selection and damage identification given measured data in Monte Carlo like way. However, most of these approaches utilize the expensive numerical algorithms which can be applied only for small dimensional problems. For further mention of existing methodologies on uncertainty quantification and damage assessment, the interested reader is referred to the rich review [START_REF] Simoen | Dealing with uncertainty in model updating for damage assessment: A review[END_REF]. This paper will focus on three of the existing methods: two of them from the deterministic parametric residual minimization approach and one from the Bayesian inference. In this article the new Bayesian technique following the ideas presented in [START_REF] Rosic | Parameter identification in a probabilistic setting[END_REF] is suggested to be used for high dimensional problems due to its small computational cost. The method can fully employ linear algebra apparatus and hence can re-use existing fast linear algebra computational libraries.

Post-tensioned concrete beams are widely used in civil engineering structures, in particular for independent span viaducts with prestessed beams. Due to this reason in the French research project named "SIPRIS", an 8 meter post-tensioned concrete beam is constructed and instrumented to simulate the damage scenario and repair actions in a laboratory surrounding. To study and to compare damage detection techniques, the post-tensioned concrete beam is instrumented by a large number of devices: accelerometers, inclinometers, displacement sensors, extensometers, and optic fibers placed along longitudinal reinforcement and groove in the concrete. Herein, we focus on the strain measurement provided by optic fibers. Combining model updating techniques and the local strain information, the aim of this paper is at detecting, localizing and quantifying early damages with the help of deterministic and probabilistic numerical procedures. For this purpose three model updating techniques are discussed and compared naming: classical Tikhonov regularization [START_REF] Tikhonov | Solutions to ill-posed problems[END_REF], Constitutive Relation Error (CRE) based updating method [START_REF] Charbonnel | A robust CRE-based approach for model updating using in situ measurements[END_REF][START_REF] Waeytens | Elastodynamics model updating for the monitoring of reinforced concrete beam: methodology and numerical implementation[END_REF][START_REF] Chamoin | Goal-oriented updating of mechanical models using the adjoint framework[END_REF][START_REF] Bonnet | Three-dimensional transient elastodynamic inversion using an error in constitutive relation functional[END_REF] and Bayesian framework [START_REF] Litvinenko | Inverse problems and uncertainty quantification[END_REF][START_REF] Rosic | Parameter identification in a probabilistic setting[END_REF]. They are analyzed with respect to the measure of the information gain obtained after the updating procedure. In this paper both deterministic and probabilistic approaches are assuming that the real value of the damage parameter is deterministic and give the estimate of its expected value given data. In the latter case the Bayesian posterior distribution represents the confidence bound on the estimate and may encounter uncertainties such as modelling and measurement errors. However, these uncertainties are not quantified per se as for this task further generalizations of the problem have to be provided, and these are not the subject of this paper.

The paper is organized as follows: in Section 2, the considered model updating techniques are summarized and in Section 3 the damage detection results are compared on a real 8 meter post-tensioned concrete beam.

Outlines of model updating techniques for damage detection

A classical damage theory [START_REF] Mazars | A description of micro-and macroscale damage of concrete structures[END_REF] takes into consideration structural damages via loss of Young's modulus. By assuming the mechanical behavior of the structure to be the static linear elasticity, we further employ the model updating strategies in order to determine a cartography of the Young's modulus, and thus to localize and to quantify its variations due to damage processes. Herein, the potential damaged area is discretized using a coarse mesh "H" and the Young's modulus is considered piecewise constant on the elements. The corresponding vector of Young's modulus to be updated is denoted by E. On the other side, a more refined mesh "h" is used to solve the mechanical problem by the Finite Element Method

K(E)U = F (1)
in which the nodal displacement U is equal to prescribed displacement U d on the kinematic boundary. In the finite element setting the rigidity matrix K (symmetric and positive definite) is parameterized by a vector E, whereas the vector U (resp. F) is associated to the nodal values of the displacement (resp.

the static loading). To this equation, one may add the mathematical prediction of the sensor outputs sim = BU, in which B is a strain-displacement matrix extracting the numerical strain at the sensor location and U is the solution of Eq. ( 1).

The numerical solution sim compared to the sensor outputs mes gives the so-called model error assumed to be of much higher magnitude than the measurement error. In order to neutralize the gap between the numerical solution and the sensor outputs, we seek the optimal value of the vector E. This further coincides with an early damage detection given the local information on strain measurements. The process of estimating E given noisy measurement data is generally ill-posed [START_REF] Tarantola | Inverse problem theory and methods for model parameter estimation[END_REF], and hence requires regularization. This can be achieved in many different ways, three of which are considered in this paper. In the following sub-sections, the three approaches for regularizing the ill-posed inverse problem are discussed and presented.

Tikhonov-based technique

Damage localization and quantification can be formulated as a constrained minimization problem. The following approach has been proposed in [START_REF] Waeytens | Elastodynamics model updating for the monitoring of reinforced concrete beam: methodology and numerical implementation[END_REF]. We seek the vector E of Young modulus which minimizes the functional

J T (E) = 1 2 (BU -mes ) T (BU -mes ) + α T 2 s T (E -E ud ) T (E -E ud ) (2)
under the equilibrium constraint (i.e. the nodal displacements have to verify the mechanical problem given in Eq. ( 1)).

In Eq. ( 2) the first term of the functional J T measures the gap between the numerical solution and the sensor outputs, whereas the second term corresponds to a Tikhonov regularization [START_REF] Tikhonov | Solutions to ill-posed problems[END_REF]. The role of the latter one is to enforce the Young's modulus E to be in the vicinity of the undamaged Young modulus E ud , and thus to avoid the case of having several unphysical minima. For a physical homogeneity of the cost function terms and customization of their weights, the parameter s T and the normalized ponderation coefficient α T are introduced, respectively. The numerical value of the last one is determined according to the discrepancy principle due to Morozov [START_REF] Morozov | On the solution of functional equations by the method of regularization[END_REF]. In practice, α T is chosen such that the maximal value of the data misfit at the end of the model updating process is close to the measurement error and strictly superior. Commonly, to solve a constrained minimization problem in Eq. ( 2), one introduces the Lagrange multiplier Ψ and studies the Lagrangian function

L T (E, Ψ) = J T (E) -Ψ T (K(E)U -F) (3)
Writing the stationarity of the Lagrangian with respect to Ψ and U leads respectively to the discrete equilibrium equation ( 1) and to the so-called mechanical adjoint problem

KΨ = B T (BU -mes ) (4)
Note that the adjoint problem is quite similar to the direct problem in Eq.

(1). In both problems, the rigidity matrix and the boundary associated to the kinematic conditions are the same. Nevertheless, the specificity of the adjoint problem lies in the following facts: the kinematic conditions vanish, the loading is located at the sensor position and its amplitude is given by the data misfit.

The third stationary condition is obtained by differentiating the Lagrangian with respect to the Young's modulus, which further gives us the gradient of the

functional J T ∇J T =            ∂J ∂E 1 = -Ψ T ∂K ∂E 1 U + α T s T (E 1 -E ud ) ... ∂J ∂E nE = -Ψ T ∂K ∂E nE U + α T s T (E nE -E ud )            (5)
Its form simplifies in the linear elastic isotropic case for which some of the appearing partial derivatives read

Ψ T ∂K ∂E i U = ∂λ ∂E i Ωi div(ψ)div(u)dΩ + ∂µ ∂E i Ωi T r[ (ψ) (u)]dΩ with ∂λ ∂E i = ν (1 + ν)(1 -2ν) , ∂µ ∂E i = 1 2(1 + ν) (6)
In order to compute the functional gradient at a lower computational cost than the straightforward use of the finite-difference schemes the adjoint solution Ψ is employed.

Finally, to localize and to quantify the damage, we seek the vector E of Young's moduli minimizing the functional J T defined in Eq. ( 2). The minimization problem is solved in an iterative way by means of the steepest descent method.

At each iteration, the steps are the following:

1. Solve the direct problem (1) considering the vector E old of Young's moduli from the previous iteration, and get U;

2. Solve the adjoint problem (4) and get Ψ;

3. Compute the descent direction, i.e. the functional gradient ∇J T using Eq. ( 5) and ( 6), and determine the descent step β by a linesearch approach; 4. Update the vector E of Young's moduli:

E new = E old -β∇J T .

Constitutive Relation Error technique

The concept of the constitutive relation error was initially introduced for the estimation of the discretization error in the finite element simulations [START_REF] Ladeveze | Error estimate procedure in the finite-element method[END_REF].

Then, it has been extended to model updating purposes [START_REF] Bonnet | Inverse problems in elasticity[END_REF][START_REF] Charbonnel | A robust CRE-based approach for model updating using in situ measurements[END_REF][START_REF] Chamoin | Goal-oriented updating of mechanical models using the adjoint framework[END_REF][START_REF] Ladevèze | A localization method of stiffness errors for the adjustment of FE models[END_REF][START_REF] Ladeveze | Application of a posteriori error estimation for structural model updating[END_REF]. This framework is particularly adapted for the determination of the constitutive relation parameters, such as Young's moduli. In this technique, we distinguish a nodal displacement U satisfying the kinematic boundary conditions [START_REF] Mottershead | Model updating in structural dynamics: a survey[END_REF] and a nodal displacement V verifying the equilibrium equations in a finite element sense (8)

CU = U d (7) K(E)V = F (8)
where CU represents the extraction of the nodal displacement on the kinematic boundary with prescribed displacement U d .

The constitutive relation parameters, i.e. the vector E of Young's moduli, are determined by acquiring the energy gap between the kinematic field U and the static field V together with the misfit between the measurements to be minimized. The corresponding functional is called the "modified constitutive relation error", and is defined by

J CRE (E) = α CRE 2 s CRE (BU -mes ) T (BU -mes ) + 1 2 (U -V) T K(E)(U -V) (9)
The first term in Eq. ( 9) corresponds to the data misfit, as in the Tikhonovbased method, whereas the second term is associated to the constitutive relation error. The parameter s CRE ensures the physical homogeneity of both terms in the functional [START_REF] Beck | Updating models and their uncertainties. i: Bayesian statistical framework[END_REF]. As in Section 2.1, the normalized ponderation coefficient α CRE is determined using the discrepancy principle.

115

To solve the minimization problem, we introduce the Lagrangian given by:

L CRE (E, Ψ) = J CRE (E) -Ψ T (K(E)V -F) (10)
Finding its stationary points leads to the following system:

K(U -V) + α CRE s CRE B T (BU -mes ) = 0 K(U -V) + KΨ = 0 KV = F (11)
Hence the kinematically admissible field U, the statically admissible field V and the adjoint field Ψ are obtained solving:

(K + α CRE s CRE B T B)U = F + α CRE s CRE B T mes KV = F Ψ = V -U (12)
Let us note that, for a given vector E of Young's moduli, if the kinematically admissible field U and the statically admissible field V are equals, then the data misfit vanishes, i.e. BU -mes = 0.

Contrary to the previous method, herein the model parameters to be updated are chosen according to a constitutive relation error criteria instead of the functional gradient. First, we evaluate the contribution denoted ε i CRE of each subdomain Ω i (associated to H-size mesh), whose Young's modulus corresponds to E i , to the constitutive relation error term, involved in Eq. ( 9).

ε i CRE = 1 2 (U i -V i ) T K i (E i )(U i -V i ) ( 13
)
where U i , V i , and K i are the restrictions of U, V and K on the subdomain

Ω i .
Then, we deduce the highest constitutive relation error

ε max CRE = max i ε i CRE . The model parameter E i associated to the subdomain Ω i is updated if ε i CRE > 0.8ε max CRE .
The selection of the most erroneous parameters associated to the highest CRE contribution is called localization step and in practice, for each iteration, only one or two constitutive parameters are corrected. The regularization of the ill-posed inverse problem is indeed permitted by the alternating localization/correction step sequence that is adopted here.

To summarize, the determination of the vector E of Young's moduli minimizing the functional J CRE , defined in Eq. ( 9), is achieved iteratively. At each iteration, the steps are the following:

1. Consider the vector E old of Young's moduli from the previous iteration; 6. Deduce the new vector E new of Young's moduli.

Bayesian technique

The previously described Tikhonov and Constitutive Relation Error regularization procedures mathematically promote the idea of loss functions which measure the error in predicting the observation data in a purely deterministic way. However, the estimation of Young's moduli E, i.e. localization of damage, can be considered in another setting which allows the prior expert beliefs on unknown parameters to enter the model. Prior beliefs describe the quantity of matter before gaining any evidence, and are commonly described by a probability density function p(E). This function plays the regularization role for an inverse problem considered in a Bayesian setting

π(E|z) = constL(z|E)p(E) (14)
In this manner the process of damage localization becomes well posed and relies on the computation of the probability density function π(E|z) that best fits the data z := mes . However, one may note that the process of updating the prior belief p(E) to a more realistic description π(E|z), called posterior, comes Following the derivation in [START_REF] Litvinenko | Inverse problems and uncertainty quantification[END_REF][START_REF] Rosic | Parameter identification in a probabilistic setting[END_REF] one may rewrite the Bayes rule in terms of quadratic minimization

E a = min ϕ∈Φ ||E -w(ϕ, E)|| 2 (15)
in which the space Φ spans the set of all measurable functions w(ϕ, E) := ϕ(z(E)). As the complete space Φ cannot be represented by computations, the previous minimization is further relaxed to

E a = min ϕ∈Φ n ||E -w(ϕ, E)|| 2 (16)
in which the closed subspace Φ n ∈ Φ is spanned by polynomials of up to n-th degree

ϕ(z) = H 0 + H 1 z + ... + H k z (⊗k) + ... + H n z (⊗n) (17)
Here, the polynomial coefficients H k are symmetric and k linear tensors, whereas z (⊗k) is symmetric tensor product of z's taken k times with itself.

Satisfying the stationary condition of Eq. ( 16) one obtains the linear system of equations

        1 < z > ... < z ⊗n > < z > < z ⊗2 > ... < z ⊗(n+1) > ... ... ... ... < z ⊗n > < z ⊗(n+1) > ... < z ⊗2n >                 H 0 H 1 ... H n         =         < E > < E ⊗ z > ... < E ⊗ z ⊗n >         (18)
in which < • > denotes the mathematical expectation of the quantity of consideration:

< q >= Ω q(ω)P (dω) (19)
i.e. the high dimensional integral over the probability space Ω with the probability measure P .

In a special case when n = 1 the system in Eq. ( 18) reduces to the so-called linear Bayesian formula (see [START_REF] Litvinenko | Inverse problems and uncertainty quantification[END_REF])

E a = E f + H 1 (z -ẑ) (20)
in which H 1 plays the role of the Kalman gain, and E a and E f are the posterior and prior random variables, respectively.

Employing the random variables instead of probability density functions, the mathematical formulation given in Eq. ( 18) offers many more possibilities for the numerical computation than the Bayes rule in Eq. ( 14). For example, if the random variables E are represented by a set of samples (realizations) then Eq. ( 20) is equivalent to the well-known ensemble Kalman filter. Otherwise, if the functional approximation approach is used as described in [START_REF] Rosic | Parameter identification in a probabilistic setting[END_REF] then Eq. ( 20)

reduces to the polynomial chaos based update. The same holds for the higher order update following from Eq. ( 17) and [START_REF] Rosic | Parameter identification in a probabilistic setting[END_REF].

In this paper the random variables are represented by a set of samples (realizations) E(ω 1 ), ..., E(ω N ) such that the system matrix and right hand side in Eq. ( 18) are easily evaluated. In this manner the computation of the posterior via Eq. ( 17) becomes rather simple. This kind of approach is similar to the ensemble Kalman filter, and hence favorable for the damage detection where the computation of the predicted measurement is difficult and expensive.

To summarize, the determination of the vector E of Young's moduli in a probabilistic manner is done in several steps:

1. Choose the prior distributions for the vector of Young's moduli;

2. Represent random variables E(ω) by a set of samples E(ω 1 ), ..., E(ω N );

3. Predict the measurement data ẑ(ω 1 ), ..., ẑ(ω N ) given E(ω 1 ), ..., E(ω N );

4. Solve the system in Eq. ( 18); 5. Make the estimation of posterior by using Eq. ( 17) , i.e.:

E a = E f + H 1 (z -ẑ) + H 2 (z ⊗2 -ẑ⊗2 ) + h.o.t (21)
Note that in the previous equation the measurement error ε 1 as well as possible modelling errors ε 2 are taken into account. Namely, the predicted measurement ẑ is modelled as ẑ = z f + ε in which z f is the prediction of the structural response given uncertain Young's moduli E f , whereas ε is independent uncertain term encountering the data-model discrepancies in a following

manner ε = ε 1 + ε 2 .
The distribution of ε is taken a priori (e.g. normally distributed), but later can be updated in a similar manner as shown previously.

The main goal of this paper is not to quantify the existing discrepancies (i.e.

their distribution), but to estimate the deterministic value of the Young's modulus and its confidence bounds given structural model. In order to identify aleratoric and epistemic uncertainties, one has to pose different kind of question than the one posed here. For more information the reader is addressed to e.g. [START_REF] Schennach | Nonparametric identification and semiparametric estimation of classical measurement error models without side information[END_REF][START_REF] Goulet | Estimation of modelling errors in structural system identification[END_REF] 3 In our experiment, an optic fiber (represented in blue in Figure 1) was put along a longitudinal reinforcement at a distance of 30mm from the bottom part of the beam. To measure the strain in the concrete beam induced by the static loading, two measurements of Rayleigh signal were necessary. The first measurement (reference) is signal obtained on the unloaded configuration, whereas the second measurement is a signal coming from the configuration loaded by 10kN . Note that the room temperature is kept constant during signal Rayleigh measurements, which further allows the temperature effect to be neglected. In Figure 2, the brown curve shows the strain measurement obtained by the fiber optic which is in accordance with a three points bending test as expected. However, due to a crack a small peak appears at 0.7m from the middle of the beam.

Note that this peak can be detected by the fiber optic because the Rayleigh sensing technique has a high spatial resolution. So, this technique of measurement permits to localize accurately the location of cracks. However, the quantification of cracks (length and opening) is less obvious. Indeed, the strain transfer from the tested material to the core of the fiber depends on the structure of the optical cable, and its protective and adhesive layers. Each fiber optic cable has its own mechanical transfer function. For accurate crack opening measurements, it is essential to know the mechanical transfer function. However, its determination is not easy and requires to combine experimental results and numerical modeling. Another approach is to use model updating techniques for damage detection as those presented above (Section 2).

Damage detection results

The goal is to localize and to quantify a potential damage on the 8 meter post-tensioned concrete beam using strain sensor outputs and model updating techniques. Before the updating process, the undamaged Young's modulus E ud = 25GP a is considered in the whole structure. To restrict the number of model parameters to be updated, the potential damaged area, i.e. lower part of the beam, is decomposed into 21 subdomains along the longitudinal axis. We aim at updating the Young's modulus in each subdomain. The results obtained from the three different model updating techniques presented in Section 2 are compared. Instead of considering the primary optic fiber outputs (about 1cm resolution) as input in all the model updating techniques, we use a 40cm averaged of the optic fiber outputs (see black diamond in Figure 2). Hence, we get 35 averaged strain sensor outputs. The measurement error is estimated at 2µm/m.

In Figure 2, we observe that the simulated strain is closed to the strain sensor outputs. The highest data misfit is located at 0.7m from the middle of the beam (x = 0.7 in Figure 2), which may correspond to a damaged area. In Tikhonovbased model updating, we take s T = 1/(E 2 ud ) and we find that the normalized ponderation coefficient α T = 2 × 10 -10 allows the highest data misfit at the end of the model updating process to be closed to the measurement error and strictly superior. The simulated strain after Tikhonov updating is represented in Figure 2 (blue circle) and the updated values of the Young modulus is given in Figure 3. The highest reduction of Young modulus, estimated at 8%, is located in the subdomain 13. It exactly corresponds to the subdomain with a crack which was brought out by the optic fiber outputs (see Section 3.2).

Concerning the constitutive relation error technique, to verify the discrepancy principle we take α CRE = 4 × 10 -7 and s CRE = 1/(LHe). In Figure 4, the contribution of the constitutive relation error in each subdomain is represented at the first and the second iteration of the model updating technique. At iteration 0, the highest error ε max CRE corresponds to the error in the subdomain 13. The red line is associated to the 0.8 ε max CRE threshold. This threshold being exceeded only by the error in the subdomain 13, only the Young's modulus in subdomain 13 is updated. The most import decrease of the functional J CRE is achieved at iteration 0 when updating the Young's modulus E 13 . The functional J CRE goes from 1.02 × 10 -5 to 8.53 × 10 -6 . After iteration 1, the model updating process should be stopped. Indeed, after iteration 1 the modifications of the Young's moduli and the decreases of the functional are not significant, i.e. few MPa variation for E and about 10 -9 variation for J CRE . In Figure 5, the cartography of the updated Young modulus is given. As observed in Figure 3, contrary to the Tikhonov strategy the constitutive relation error technique only updates the Young's modulus in subdomain 13. Lastly, we note that both deterministic methods give the same estimation of the damage level in this subdomain, i.e. 8% reduction of the Young's modulus.

Finally, the identification of the local damage by approximations of Bayesian inference is undergone by assuming the normalized Young's moduli to be independent and uncorrelated Gaussian random variables with mean equal to 0.9 and standard deviation 0.1, see Figure 6. This choice is taken against logical reason saying that the normal distribution is not a natural choice for the description of the positive defined material properties. However, one may argue that the probability of negative value occurring for small standard deviations compared to the mean value is relatively small. The posterior is estimated using quadratic and linear Bayesian formula and ensemble of 100 samples (i.e. 100 forward simulations). The posterior results are similar to those obtained by the previous two procedures. Namely, the damaged area is appearing to be at the place of the subdomain 13. In this section the Young's modulus obtains the lowest value and predicts the damage of around 9%, see Figure 7. The information coming from the measurement data are reducing uncertainty in the prior, and hence the posterior is characterized by much smaller standard deviation. Comparing the modes of the numerically simulated strain given linear and quadratic posteriors to the measured value of the strain, one may observe that the linear estimate is better in most of subdomains besides in the 13th.

There the quadratic estimate is better, see Figure 9. This also can be supported by graphs of 5 sigma regions in Figure 10 which are drawing the same conclusion. Namely, the quadratic probability region is including all the values of the measured strain, whereas the linear does not. The reason for this is more than obvious: linear update is the optimal estimate when the relation between the measurement and the parameter being identified is linear. However, note that the quadratic estimate is also not the optimal but closer. To obtain the best estimate one would have to compute higher degree updates. This is the subject of the further study which has also to take into account more realistic models for the prior.

Conclusions and prospects

Three model updating techniques, involving different kind of regularization, have been compared; their respective ability to localize and to quantify early damages has been illustrated using the results of an optical fiber sensor instrumenting a 8 meters post-tensioned concrete beam. Even though the procedures have used different kinds of estimation, they successed to localize the damage area, as well as to give similar estimation of its value. In the present work, Young modulus was considered as a damage indicator. In practice, the proposed techniques can be employed in a long term monitoring strategy. Each month, the model updating process can be run to actualize the cartography of Young's moduli in the structure. Hence, the site manager can follow the time evolution of the damage indicator. In case of abnormal variation, an alarm could be generated and a visual inspection could be recommended.

2 .

 2 Determine the kinematically admissible field U solving the first equation in (12); 3. Determine the statically admissible field V solving the second equation in (12); 4. Evaluate the constitutive relation error ε i CRE from Eq. (13) in each potential damaged subdomain Ω i and get the highest constitutive relation error ε max CRE ; 5. Update only the model parameters E i verifying ε i CRE > 0.8ε max CRE . If only one model parameter is selected, a standard linesearch is employed to minimize the functional J CRE . Otherwise, the steepest descent method is applied.

  on the expense of solving a probabilistic problem. Namely, p(E) is corrected by the information gain L(z|E) describing of how likely the measurement data are given the model. The last one is essentially the measure of distance zẑ, e.g. for normally distributed likelihood one has L(z|E) = N (z -ẑ, C ), which implicitly includes the prediction of the measurement data ẑ given E.

3 . 1 .

 31 . Comparison of damage detection results on a 8 meter post-tensioned concrete beam Description of the static test on the 8 meter post-tensioned concrete beam Let us consider a 8 meter post-tensioned concrete beam represented in Figure 1 which is simply supported on two elastomer supports. The length, the height and the depth are L = 8m, H = 0.6m and e = 0.3m. A static loading of10kN is applied at the middle of the beam. Initially, only the two longitudinal cables are activated to prestress the concrete beam. After damaging strongly the structure by means of an overload (higher than 10kN ), the curved cable was used to close the cracks and to repair the beam. In the considered experiments, we focus on the detection of early damage which may not be detected by visual inspection.

Figure 1 :

 1 Figure 1: Post-tensioned concrete beam of 8 meter length by means of 3 steel cables (in red) -Optic fiber (in blue) on a longitudinal armature

Figure 2 :

 2 Figure 2: Strain sensor outputs from optic fiber (in brown) and simulated strain before and after Tikhonov-based model updating process

Figure 3 :

 3 Figure 3: Normalized Young modulus E/E ud in each subdomain after deterministic model updating strategies -Comparison of Tikhonov and Constitutive Relation Error results.

Figure 4 :

 4 Figure 4: Constitutive Relation Error in each subdomain at iteration 0 (at left) and at iteration 1 (at right) -At iteration 0, Young modulus in subdomain 13 has to be updated.

Figure 5 :

 5 Figure 5: Cartography of Young modulus after Constitutive Relation Error-based model updating process -the damaged area is in red and it corresponds to a 8% damage.

Figure 6 :

 6 Figure 6: 95% interval describing prior distribution

Figure 7 :

 7 Figure 7: Probability densitiy functions of linear and quadratic Bayesian estimate of Young's modulus in subdomain 13 compared to the prior assumption.

Figure 8 :

 8 Figure 8: Comparison of posterior modes obtained after linear and quadratic Bayesian update to the observed value of the strain

Figure 9 :

 9 Figure 9: 95% probability interval describing linear (left) and quadratic (right) posteriors

Figure 10 :

 10 Figure 10: Five strandard deviations region of estimated strain simulated using: linear posterior (left) and quadratic posterior (right)

Acknowledgments

This research was supported by Conseil Régional d'Ile de France. This work was conducted as part of the project "SIPRIS", which involves several partners: Advitam (Vinci), ASF (Highways of the South of France), CEDRAT, ESIEE-Paris, NECS, SYROKKO and IFSTTAR.


Rachid Riah

Mirko Fiacchini

Mazen Alamir

Domain of attraction estimation of cancer chemotherapy model affected by state proportional uncertainty

Keywords: Cancer, Chemotherapy, Uncertainty, Invariance, Set theory

This paper proposes an analytical method for the characterisation of the set of all initial tumor and normal states of a cancer dynamical model for which there exists a successful drug administration therapy. The dynamical model considered simulates the behaviour of the tumor and the normal cells populations growth in presence of a chemotherapy drug. This model is assumed to be affected by a state proportional uncertainty. The successful drug administration therapy is attained by respecting the tumor eradication and avoiding the lower admissible bound of the normal cells population which is viewed as an indicator of the patients health. The problem is posed in terms of domain of attraction computation and is solved by applying set theory and invariance tools.

I. INTRODUCTION

On cancer research, the mathematical modeling is viewed as the pillar that have a big impact on the success of the cancer treatment. Some of the models are based on the evolution of different cells, see [START_REF] Afenya | Acute leukemia and chemotherapy: a modeling viewpoint[END_REF], [START_REF] De Pillis | Mixed immunotherapy and chemotherapy of tumors: modeling, applications and biological interpretations[END_REF], [START_REF] Hahnfeldt | Tumor development under angiogenic signaling a dynamical theory of tumor growth, treatment response, and postvascular dormancy[END_REF], [START_REF] Matveev | Application of optimal control theory to analysis of cancer chemotherapy regimens[END_REF]. Such models may incorporate the effect of external drugs on the evolution of the cancer, as chemotherapy and immunotherapy drug or their combination. In the recent years, mathematical models of cancer have been widely used also for controlbased tumoral therapy design, applying optimal control [START_REF] De Pillis | Chemotherapy for tumors: An analysis of the dynamics and a study of quadratic and linear optimal controls[END_REF], [START_REF] Ledzewicz | Optimal control for combination therapy in cancer[END_REF], or feedback control, see [START_REF] Alamir | On probabilistic certification of combined cancer therapies using strongly uncertain models[END_REF], [START_REF] Chareyron | Mixed immunotherapy and chemotherapy of tumors: Feedback design and model updating schemes[END_REF].

In this paper, the cancer chemotherapy model to be studied is inspired from the work [START_REF] Afenya | Acute leukemia and chemotherapy: a modeling viewpoint[END_REF], [START_REF] Matveev | Application of optimal control theory to analysis of cancer chemotherapy regimens[END_REF]. It represents the evolution of normal and tumor cells populations under chemotherapy drug therapy. In [START_REF] Afenya | Acute leukemia and chemotherapy: a modeling viewpoint[END_REF], this model is used to represent the evolution of acute myeloblastic leukemia (AML) and in [START_REF] Matveev | Application of optimal control theory to analysis of cancer chemotherapy regimens[END_REF] it has been used to simulate the cancer evolution. The normal cells represent a part of the innate immune system of the human body and can be interpreted as the aggregation of NK cells, CD8+ cells and a circulating lymphocytes (or white blood cells) as defined in [START_REF] Chareyron | Mixed immunotherapy and chemotherapy of tumors: Feedback design and model updating schemes[END_REF].

The uncertainties affecting the cancer models, either parametric [START_REF] Alamir | Robust feedback design for combined therapy of cancer[END_REF], or additive are crucial for the analysis and control design in order to find the successful drug administration therapy. Therefore, the main objective of the paper, whose preliminary version is [START_REF] Riah | Invariance-based analysis of cancer chemotherapy[END_REF], consists in the development of a new method for chemotherapy design considering that the model is affected by state proportional uncertainties. The state proportional uncertainties are the uncertainties which affect proportionally the dynamic of each state of the model. This method is based on the properties of robust invariant sets for Convex Difference Inclusions (CDI) systems.

In this paper, our contribution is substantially based on set theory and invariance for control. Stability and robustness will be addressed using the properties of invariant sets. Many results regarding invariance and related topics have been provided in the literature, see [START_REF] Bertsekas | Infinite-time reachability of state-space regions by using feedback control[END_REF], [START_REF] Gilbert | Linear systems with state and control constraints: The theory and application of maximal output admissible sets[END_REF], [START_REF] Kolmanovsky | Theory and computation of disturbance invariant sets for discrete-time linear systems[END_REF], [START_REF] Fiacchini | On the computation of convex robust control invariant sets for nonlinear systems[END_REF] and [START_REF] Blanchini | Set-Theoretic Methods in Control[END_REF]. Invariance is related to reachability [START_REF] Matthias | Reachability analysis of nonlinear systems with uncertain parameters using conservative linearization[END_REF], and viability. Typically, the viability kernel is much harder to compute [START_REF] Shahab | Computing the viability kernel using maximal reachable sets[END_REF].

The new method provides the domain of attraction of the cancer chemotherapy model affected by a state proportional uncertainty. This domain contains all tumor and normal states for which a set of admissible administration chemotherapy profile exists, that drives theses states to the safe region. The safe region is defined to be the set for which the number of tumor cells is small enough and the normal cells number higher than the minimal admitted level. A minimal level is imposed as a health condition.

The paper is organized as follows: First the mathematical model of cancer is recalled in Section II. Section III presents the problem statement addressed in the paper. The properties and numerical methods related to controlled robust invariant sets and convex difference inclusions systems are given in Section IV. The application to the cancer problem of the methods developed in Section IV is proposed in Section V. Section VI illustrates the results given in Section V. Section VII ends the paper providing conclusions and giving hints for future investigations.

Notation:

Given n ∈ N, define N n = {m ∈ N : 1 ≤ m ≤ n}. Given A ∈ R n×m , A i with i ∈ N n denotes its i-th row, A i, j
with i ∈ N n and j ∈ N m a value that corresponds to the i-th row and the j-th column.

II. MATHEMATICAL CANCER MODEL

We consider the mathematical model that has been developed in [START_REF] Matveev | Application of optimal control theory to analysis of cancer chemotherapy regimens[END_REF], [START_REF] Afenya | Acute leukemia and chemotherapy: a modeling viewpoint[END_REF] and [START_REF] Afenya | Some perspectives on modeling leukemia[END_REF] in order to describe the dynamical behaviour of the normal and tumor cells populations under a chemotherapy drug. This model involves the following cells populations

• L: tumor cells population,

• N: normal cells population.

The assumptions considered on the model are given in the following, see [START_REF] Matveev | Application of optimal control theory to analysis of cancer chemotherapy regimens[END_REF], [START_REF] Afenya | Acute leukemia and chemotherapy: a modeling viewpoint[END_REF] and [START_REF] Afenya | Some perspectives on modeling leukemia[END_REF] for more details.

1) The growth dynamics of both tumor and normal cells populations are the same anywhere in the body; 2) The tumor and normal cells populations follow a process of Gompertizian growth;

3) The tumor cells exhibit a negative and inhibiting effect on the growth and development of the normal cells; 4) The chemotherapy drug kills both tumor and normal cells and affects the normal more than the tumor cells; 5) The drug spreads instantaneously within the body. The resulting model, based on the above assumptions, is defined as follows

L = αL ln θ L L -Π 1 (u)L, L(0) = L 0 , Ṅ = β N ln θ N N -Ξ(L)N -Π 2 (u)N, N(0) = N 0 , (1)
where the description of each term is given in Table I. Following the reasonings in [START_REF] Afenya | Acute leukemia and chemotherapy: a modeling viewpoint[END_REF], [START_REF] Matveev | Application of optimal control theory to analysis of cancer chemotherapy regimens[END_REF] and [START_REF] Afenya | Some perspectives on modeling leukemia[END_REF], the assumptions below can be done to simplify the model.

• The function Ξ(L) is taken to be linear as in [START_REF] Afenya | Acute leukemia and chemotherapy: a modeling viewpoint[END_REF] (i.e.

Ξ(L) = γL, where γ is a given parameter); • We consider that the loss function Π (u) is linear, taking into account the saturation phenomena, i.e, Π 1 (u) = ku, k is the fraction of tumor cells that are killed due to the chemotherapy drugs; • We assume that the loss function Π 2 (u) is linear, i.e, Π 1 (u) = lu, l is the fraction of the normal cells that are destroyed by the effects of the chemotherapy drugs (it is assumed that l < k); • We assume also that the dynamics of the cells populations are affected by state proportional uncertainties. Physically, the state proportional uncertainties can be viewed as a disturbance that grows proportionally with the number of the cells populations. It is in fact reasonable to suppose that these uncertainties should depend on the number of the cells populations, since its effect on the dynamic of the model must be suitable for small and high number of the cells populations.

The resulting model becomes

L = αL ln θ L L -kuL + w 1 L, L(0) = L 0 , Ṅ = β N ln θ N N -γLN -luN + w 2 N, N(0) = N 0 , (2)
where (w 1 w 2 ) T ∈ W ⊆ R 2 is an unknown uncertainty belonging to a bounded set W . The bounded set W can be taken as a box, a polyhedral set, an ellipsoidal set or their intersection.

As proposed in [START_REF] Matveev | Application of optimal control theory to analysis of cancer chemotherapy regimens[END_REF], we introduce the change of variables x 1 = ln θ L L and x 2 = ln θ N N , and we get the equivalent system

ẋ1 = -αx 1 + ku -w 1 , x 1 (0) = x 10 , ẋ2 = -β x 2 + γθ L e -x 1 + lu -w 2 , x 2 (0) = x 20 . (3)
The discrete-time system modelling the cancer evolution is obtained by sampling the continuous system (3), with numerical integration schemes using Euler's method and sampling time of T s = 1 day. It results in

x + 1 = (1 -T s α)x 1 + T s ku -T s w 1 = f 1 (x, u) -T s w 1 , x 1 (0) = x 10 , x + 2 = (1 -T s β )x 2 + T s γθ L e -x 1 + T s lu -T s w 2 = f 2 (x, u) -T s w 2 , x 2 (0) = x 20 . (4)
The model (4) involves two state variables x 1 , x 2 , one manipulated variable u and two uncertainties variables w 1 and w 2 . Notice that in this space, the uncertainty appears as additive terms. The consistent values for the parameters involved in this model are inferred from [START_REF] Afenya | Acute leukemia and chemotherapy: a modeling viewpoint[END_REF] and [START_REF] Matveev | Application of optimal control theory to analysis of cancer chemotherapy regimens[END_REF], then modified to reproduce the cells evolutions given in [START_REF] Afenya | Acute leukemia and chemotherapy: a modeling viewpoint[END_REF] and summarized in Table II.

III. PROBLEM STATEMENT

Consider the nonlinear model ( 2) that describes the dynamics of the cells populations in presence of chemotherapy drug and state proportional uncertainties.

The system is subject to the following constraints on system states, control inputs and uncertainty variables.

(L N) T ∈ Y ⊆ R 2 , (5)
u ∈ U ⊆ R = {u ∈ R m : u min ≤ u ≤ u max }, (6)
w ∈ W ⊆ R 2 . (7)
In cancer therapy drug administration, the main objective is to reduce the tumor cells population L, while maintaining the number of the normal cells population N above a prescribed level, see [START_REF] Chareyron | Mixed immunotherapy and chemotherapy of tumors: Feedback design and model updating schemes[END_REF], [START_REF] De Pillis | Mixed immunotherapy and chemotherapy of tumors: modeling, applications and biological interpretations[END_REF]. The initial values of (L N) T belong to Y . The number of the normal cells population N is used in [START_REF] Matveev | Application of optimal control theory to analysis of cancer chemotherapy regimens[END_REF] as a measure of the patient health.

In this paper we would like to determine all the initial tumor and normal states for which there exists an appropriate drug injection profile u satisfying the constraint (6) for all uncertainties realizations w belonging to W . This therapy must lead to a substantial regression of the tumor size, i.e. L becomes small enough, while avoiding that the health measure reaches dangerous values for the patient, i.e. N ≥ N min , with N min denoting the minimal admitted value of the normal cells.

In the space of x 1 and x 2 , the state constraint (5) becomes

(x 1 x 2 ) T ∈ X ⊆ R n , ( 8
)
where X is determined from the set Y applying the change of variable defined above.

In this space, we would like to determine the set of the initial points (x 1 , x 2 ) ∈ X, such that the dynamics of the system (4) can be driven to the safe region for all param. value param. value uncertainties realizations w belonging to W , by applying an admissible chemotherapy drug profile u satisfying the constraint [START_REF] Blanchini | Set-Theoretic Methods in Control[END_REF]. The safe region is the set in the state space x 1 , x 2 defined by x 1 high enough and x 2 ≤ x 2max , where x 2max = ln θ N N min denotes the maximal admitted value of x 2 . These initial tumor and normal states represent the domain of attraction of the cancer system (4). Its computation is substantially based on the properties of positively invariant sets and convex inclusions, that will be introduced in the next section.

α 3.96 × 10 -4 day -1 β 3.33 × 10 -2 day -1 θ N 1.4 × 10 8 cells θ L 3 × 10 8 cells k 8 × 10 -2 l 15 × 10 -3 γ 10 -9 (cells.day) -1

IV. CONVEX DIFFERENCE INCLUSIONS SYSTEMS AND DOMAIN OF ATTRACTION

In this section, the method that provides the domain of attraction of a nonlinear systems is stated. This method is based on the approximation of the nonlinear system with a Convex Difference Inclusions (CDI) system and also on the properties of positively invariant sets. In practice, this method computes the maximal contractive invariant convex set as an estimation of the domain of attraction of a nonlinear system.

The convex difference inclusions systems are defined by a particular class of set valued maps as a dynamic function. The set valued maps that can be used to define this class of systems is such that, given a point in the state space, its image through the map is a convex and compact set.

Let us consider the uncertain difference inclusions system

x + ∈ F W (x, u) = F (x, u) +W (9)
where

x ∈ X ⊆ R n is the state, x + is the successor, u ∈ U ⊆ R m is the control input and F (•, •) is a set valued map on R n , that is F (x, u) represents a function which relates a set to every point (x u) T ∈ R n × R m , i.e. F (x, u) ⊆ R n for all x ∈ R n .
The sets X, W and U represent the set of state constraints, the set of additive uncertainties and the set of admissible control inputs respectively. The set W is an additive uncertainty bounding set satisfying the following assumption. Assumption 1: The set W ⊂ R n is compact, convex and 0 ∈ W .

The interest of an uncertain difference inclusions systems lies in the fact that they can be used to approximate the uncertain nonlinear systems x + = f (x, u) + w, i.e. f (x, u) + w ∈ F W (x, u) for all x ∈ X, all w ∈ W and for all u ∈ U.

Since a robust invariant set for the CDI system is robust invariant also for the nonlinear system then the approximation may be useful in order to determine the domain of attraction of the nonlinear system. Although some conservatism is in general introduced, it is arbitrarily small for the case under analysis.

Notice that the important tool to deal with convex closed sets and convexity is the support function. For some properties of support functions readers can see [START_REF] Rockafellar | Convex Analysis[END_REF].

Definition 1: Given a set Ω ⊆ R n , the support function of

Ω evaluated at η ∈ R n is φ Ω (η) = sup x∈Ω η T x.
The support function permits to transform the set-inclusion condition in terms of linear inequalities, see [START_REF] Rockafellar | Convex Analysis[END_REF] for instance. From the definition of the support function, we get this property.

Property 1: Given two closed, convex sets Ω ⊆ R n and Γ ⊆ R n , then x ∈ Ω if and only if η T x ≤ φ Ω (η) for all η ∈ R n , and Γ ⊆ Ω if and only if φ Γ (η) ≤ φ Ω (η), for all η ∈ R n .

In order to define the uncertain difference inclusions systems, let us introduce this assumption.

Assumption 2: Assume that the set valued map F determining the system dynamics ( 9) is such that F (x, u) is compact and convex for all (x u) ∈ R n × R m and for every η ∈ R n , the function F(x, u, η) : R n × R m × R n -→ R defined as

F(x, u, η) = φ F (x,u) (η) = sup z∈F (x,u) η T z, (10)
is convex with respect to (x u) on X ×U. Thus, in practice, F(x, u, η) is the support function of the set F (x, u) evaluated at η ∈ R n and then

F (x, u) = {z ∈ R n : η T z ≤ F(x, u, η), ∀η ∈ R n }, ( 11
)
and it is convex in (x u) T . Since the set-valued map F W (x, u) determining the system (9) is the sum of two convex sets, then it is convex also with respect to (x u) T ∈ R n × R m . Thus, the set valued map determining the system ( 9) is such that

F W (x, u) = {z ∈ R n : η T z ≤ F(x, u, η) + φ W (η), ∀η ∈ R n }. ( 12
)
The dynamical systems [START_REF] De Pillis | Mixed immunotherapy and chemotherapy of tumors: modeling, applications and biological interpretations[END_REF] for which Assumption 2 holds are tightly related to the uncertain CDI systems defined in [START_REF] Fiacchini | Convex difference inclusions for systems analysis and design[END_REF], [START_REF] Fiacchini | Invariant sets computation for convex difference inclusions systems[END_REF]. Now, we recall the standard definitions of the controlled robust invariant set for the generic nonlinear system, adapted here for the case of potentially unbounded sets that do not necessarily contain the origin and set valued maps.

Definition 2 ( [START_REF] Blanchini | Set-Theoretic Methods in Control[END_REF]): The closed convex set Ω ⊆ R n is a controlled robust invariant set for the uncertain nonlinear system x + = f W (x, u) = f (x, u) + w if for every x ∈ Ω, there exists u ∈ U such that f W (x, u) ∈ Ω for all w ∈ W ; it is controlled robust invariant for the CDI system (9) if for every x ∈ Ω, there exists u ∈ U such that F W (x, u) ⊆ Ω.

Whatever the additive uncertainties w ∈ W affecting the trajectory started in a controlled robust invariant set Ω, there exists an admissible control input, such that this trajectory remains inside Ω.

For the system (9), the robust controlled invariance of Ω ⊆ X can be tested using the controlled robust one-step backward operator, defined below.

Definition 3: Consider the closed convex set Ω and assume that Assumption 2 holds for the set valued map F (•, •) determining the dynamic system [START_REF] De Pillis | Mixed immunotherapy and chemotherapy of tumors: modeling, applications and biological interpretations[END_REF]. The controlled robust one-step operator is defined as follows

Q W (Ω,U, X) = {x ∈ X : ∃u ∈ U : F W (x, u) ⊆ Ω}. (13)
Hence, the controlled robust one-step operator associates to every set Ω the set of points that are mapped inside Ω through F W (•, •) applying an appropriate control input u ∈ U. Thus, by definition, a closed convex set Ω ⊆ X is robust invariant if and only if Ω ⊆ Q W (Ω,U, X) and the controlled robust one-step operator can be used to check controlled robust invariance and to compute increasing sequences of nested controlled robust invariant sets as for linear or CDI systems. Indeed, Algorithm 1, standard for generating a sequence of convex sets that converges to the domain of attraction, see [START_REF] Blanchini | Set-Theoretic Methods in Control[END_REF], can be applied also in this context.

Algorithm 1 Sequence of convex sets converging to the maximal controlled robust invariant set for the system (9)

Input: Initial convex closed controlled robust invariant set Ω 0 ⊆ X. 1: for k ∈ N N do 2: Compute Ω k+1 = Q W (Ω k ,U, X) ∩ X 3: if Ω k+1 = Ω k then Ω max = Ω k return 4:
end if 5: end for Notice that Algorithm 1 can be used starting from the set of state constraints Ω 0 = X or from an initial controlled robust invariant set Ω 0 . Thus, in the latter case, Ω k are controlled robust invariant set and converges to the domain of attraction for the system (9), and hence for all nonlinear systems approximated by [START_REF] De Pillis | Mixed immunotherapy and chemotherapy of tumors: modeling, applications and biological interpretations[END_REF]. Finally, being our objective to have a sequence of Ω k that are polyhedral, we introduce the following results, proved in [START_REF] Fiacchini | Convex difference inclusions for systems analysis and design[END_REF], which is functional for this purpose.

Proposition 1: Consider the set-valued map F W (x, u) determining the system dynamics [START_REF] De Pillis | Mixed immunotherapy and chemotherapy of tumors: modeling, applications and biological interpretations[END_REF]. Given the state constraint set X, the control input set U and the polytope Ω = {x ∈ R n : Hx ≤ h}, with H ∈ R n h ×n , the controlled robust onestep operator is

Q W (Ω,U, X) = {x ∈ X : ∃u ∈ U : F(x, u, H i T ) ≤ h i -φ W (H T i ), ∀i ∈ N n h }. (14)
Proof: From Definition 2, we have that x ∈ Q W (Ω,U, X) if and only if ∃u ∈ U such that F W (x, u) ⊆ Ω, which is the condition given to determine the controlled robust one-step operator. Thus, using the Property 1 and Eq.( 12) the latter condition is equivalent to x ∈ Q W (Ω,U, X) if and only if ∃u ∈ U such that φ F (x, u)(η) + φ W (η) ≤ φ Ω (η) for all η ∈ R n . For polyhedral case, this condition is equivalent also to x ∈ Q W (Ω,U, X) if and only if ∃u ∈ U such that F(x, u, H i T ) ≤ h iφ W (H T i ), ∀i ∈ N n h which is proved in [START_REF] Fiacchini | Convex difference inclusions for systems analysis and design[END_REF]. Finally, the controlled robust one-step operator can be defined as given by Eq.( 14).

Notice that, as proved in [START_REF] Fiacchini | Convex difference inclusions for systems analysis and design[END_REF] for general convex closed sets, for every polytopic Ω, the set Q W (Ω,U, X) is closed and convex. Moreover, if F(•, •, η) are piecewise affine functions of (x u) T , then Q W (Ω,U, X) is a polyhedron, that is the intersection of a finite number of halfspaces, for every polytopic (or polyhedral) Ω.

Finally, once the controlled robust one-step operator is given, Algorithm 1 can be used to compute the domain of attraction of the nonlinear systems x + = f (x, u) + w, with w ∈ W .

V. CONTROL ROBUST ONE-STEP OPERATOR FOR CANCER

CHEMOTHERAPY

In this section, the controlled robust one-step operator for the uncertain cancer chemotherapy model ( 4) is stated. Notice that the model ( 4) is in the form x + = f (x, u) + w. Then given the set Ω = {x ∈ R 2 : Hx ≤ h}, with H ∈ R n h ×2 , the controlled robust one-step operator of this model can be computed using Eq.( 14). Eq.( 14) depends on the bounding functions F(x, u, H i ) for all H i T ∈ R 2 , i ∈ N n h , that must first be determined.

First, we will compute the bounding functions F(x, u, η) for all η ∈ R 2 . This function has to be convex with respect to (x u) T , and such that the associated set-valued maps F (x, u) as [START_REF] Fiacchini | On the computation of convex robust control invariant sets for nonlinear systems[END_REF], satisfies f (x, u) ∈ F (x, u), with f (x, u) in ( 4), for every (x u) T ∈ X ×U. Thus, according to Property 1, the functions F(x, u, η) must satisfy the following inequality η T f (x, u) ≤ F(x, u, η), for every η ∈ R 2 . Moreover, since the nonlinearity in (4) involves only f 2 (x, u) then the bounding functions of η T f (x, u) are related only to η 2 , since the dynamics of x 1 is linear. A possible choice is

• if η 2 ≥ 0: F(x, u, η) = η T f (x, u), since F(x, u, η) is already convex in this case; • if η 2 < 0: F(x, u, η) = η 1 (1 -T s α)x 1 + T s ku + η 2 (1 -T s β )x 2 +T s γθ L (ax 1 + b) + T s lu , ( 15
)
with a and b such that (ax 1 + b) ≤ e -x 1 , obtained for instance as the tangent to the graph of e -x 1 at one point; for all η ∈ R 2 .

We are interested in convex piecewise affine bounding functions for computational purposes. For this it is sufficient to replace e -x 1 with a convex piecewise affine upper bound, which is easily obtainable, for the case of η 2 ≥ 0. In fact, it is sufficient to choose a set of q ∈ N parameters c i ∈ R,

d i ∈ R, with i ∈ N q , such that, e -x 1 ≤ max i∈N q {c i x 1 + d i }, for all
x 1 ∈ R. To obtain those parameters, it is sufficient to define z i with i ∈ N q+1 such that z i < z i+1 and then

c i = e -z i+1 -e -z i z i+1 -z i , d i = z i+1 e -z i -z i e -z i+1 z i+1 -z i
are the parameters of an affine function such that e -x 1 ≤ c i x 1 + d i for every x 1 ∈ [z i , z i+1 ] and e -z i = c i z i + d i for all i ∈ N q , see Figure 1.

Remark 1: Notice that, in order to increase the precision, the number of parameters q can be taken as big as desired. Hence, arbitrary precision can be attained.

Concerning the case of η 2 < 0, no modification is required with respect to [START_REF] Kolmanovsky | Theory and computation of disturbance invariant sets for discrete-time linear systems[END_REF], since the function ( 15) is already piecewise affine function. Then, we obtain

F(x, u, η)=        η 1 (1 -T s α)x 1 + T s ku + η 2 (1 -T s β )x 2 +T s γθ L max i∈N q {c i x 1 + d i } + T s lu if η 2 ≥ 0, F(x, u, η) as in (15) if η 2 < 0. ( 16
)
Figure 1 shows the convex upper bound and the concave lower bound of e -x 1 implicitly used to determine F(x, u, η). Notice that F(x, u, η) is a convex upper bound of η T f (x, u), for all (x u) T ∈ X × U and for all η ∈ R 2 . Once a convex upper bound of η T f (x, u) is determined, the controlled robust one-step operator will be given using the Proposition 1.

Remark 2: By applying Algorithm 1, we found that the sequence of polytopes Ω k = {x ∈ R n : H k x ≤ h k } generated by the controlled robust one step operator never have a facet determined by H k i with H k i,2 < 0 except the case related to the trivial constraints x 2 ≤ 0. Since the latter constraints can be neglected in the computation, then the lower bound of e -x 1 is never used. This would mean that the sequence of polytopes obtained and the domain of attraction are affected only by the mismatches between e -x 1 and the piecewise function max i∈N q {c i x 1 + d i }, mismatch that can be done arbitrarily small, as notice in Remark 1. Hence, the desired precision can be achieved by employing sufficiently close piecewise approximations of e -x 1 . Consequently, the conservatism introduced by approximating the uncertain nonlinear cancer system (4) by the uncertain CDI system (9) can be reduced increasing the precision as indicated in Remark 1.

According to Remark 2, the piecewise affine bound function of η T f (x, u) to be employed is the one given for η 2 0. Therefore, the related controlled robust one-step operator defined in Proposition 1becomes

Q W (Ω,U, X) = x ∈ X | ∃u ∈ U : ∀i ∈ N n h , ∀ j ∈ N q , H i,1 (1 -T s α) + H i,2 T s δ c j x 1 + H i,2 1 -T s β x 2 ≤ h i -H i,2 T s δ d j -H i,1 T s k + H i,2 T s l u -φ W (H T i ) . (17)
Then Q W (Ω,U, X) maps polytopes in polytopes and the Algorithm 1 generates a sequence of polytopes that converge to the maximal controlled robust invariant polytopic set for the system (4).

VI. SIMULATIONS

Let us assume that the chemotherapy drug profile is constrained to take value from U such that u min = 0 and u max = 1. We also consider that the minimal admitted level of the normal cells is N min = 10 7 Cells, which implies,using the change of variable that x 2max = 2.30. Since the exact values of the model parameters are not available in the literature, the value of N min is also taken approximately as the others parameters described and noted in Section II.

Let us also assume that the polyhedral set of the state constraints in the domain of L and N is Y = {y = (L N) T ∈ R 2 : L ≥ 0; N min ≤ N}, which implies, using the change of variables that X = {x ∈ R 2 :

x 1 ≥ 0; 0 ≤ x 2 ≤ x 2max }.
The state proportional uncertainty w is assumed to belong to the set W = {w ∈ R 2 : w 2 0.06; 0 w 1 ; w 2 0}. This set corresponds to the worst case that can affect the dynamics of the cells populations. In the model (2), the tumor cells population L without chemotherapy drug diverges. The worst case for the uncertainties increases the divergence of the tumor cells population L to infinity, i.e. w 1 is taken positive. The normal cells population needs to grow and the worst case of the uncertainties prevent this growth, i.e. w 2 is taken negative.

Once the set of constraint states X and the initial invariant set Ω 0 are selected, Algorithm 1 is fully automatic and no parameters need to be selected. Now, applying Algorithm 1 with the controlled robust onestep operator given by Eq.( 17), starting from an initial set Ω 0 , we get a sequence of convex sets that converges to the maximal controlled robust invariant set for the cancer system (4).

In Figure 2, the sequence of nested controlled robust invariant sets that converge to the maximal controlled robust invariant set is depicted in black line whereas the maximal controlled robust invariant set is depicted in blue line. This sequence of polytopes are computed starting from the set Ω 0 = {x ∈ X : x 1 ≥ a} with a big enough, for which the invariance condition Ω 0 ⊆ Q W (Ω 0 , X,U) holds. Consequently, the controlled robust one-step operator can be used to check controlled robust invariance.

Notice that Ω k+1 is defined as the set that contains all the states for which there exists a set of admissible chemotherapy drug such that these states will be mapped in Ω k in one step regardless of state proportional uncertainties. In the same figure, the maximal controlled invariant set for the cancer system (4) is depicted in red line.

Figure 3 shows the maximal controlled robust invariant set in blue line and the maximal controlled invariant set in red line. Notice that the maximal controlled robust invariant set is calculated using Algorithm 1 starting from the set of state constraints X. The sequence of polytopes is depicted in black line.

By applying the change of variable defined above to the sets given in Figure 2, one gets the tumor-normal cells domains of attraction of the cancer system (2), that are depicted in Figure 4. In blue line, the domain of attraction of the uncertain cancer system (2) is showed. This contains all the tumor and normal states that can be driven to the set of healthy states considering that the dynamical cancer chemotherapy model is affected by state proportional uncertainties. In the same figure, the tumor-normal cells domain of attraction of the cancer system (2) without the state proportional uncertainties is given by the red line.

Remark from Figure 4 the drastic reduction of the domain of attraction caused by the state proportional uncertainties. Therefore, in order to achieve a successful chemotherapy treatment, the identification of the real values of the model parametersis crucial.

These domains of attraction represent, with arbitrary precision, the exact set of all the initial tumor and normal states for which there exists a specified set of administration therapy profile. This therapy leads to the eradication of the tumor cells without violating the constraint on the patient's health. Therefore, given the initial state of a patient, one could infer on the existence of a successful chemotherapy treatment. Furthermore, the amount of drugs that should be delivered and the therapy length can be calculated.

VII. CONCLUSION

In this paper, we developed a method to determine the domain of attraction for cancer chemotherapy model. This domain contains all the initial tumor and normal states for which there exists a set of admissible drug administration therapy. This method could be used to characterize which patients can be successfully cured. This domain is obtained in the presence of state proportional uncertainties affecting the model and saturation constraints which limit the amount of the drug injection. In future work, alternative modelling frameworks could be employed.

1 Fig. 1 :

 11 Fig. 1: Bounds of the exponential function.

Fig. 2 :

 2 Fig. 2: Comparison between the maximal controlled invariant polyhedral set (red line) and the maximal controlled robust invariant polyhedral set (blue line) provided with polyhedral sequences considering state proportional uncertainty w ∈ W = {w ∈ R 2 : w 2 0.06; 0 w 1 ; w 2 0} (black line), computed starting from an initial set Ω 0 . The admitted level of x 2 is equivalent to the admitted level of the normal cells, which is 10 -1 θ N cells. The chemotherapy drugs are constrained to take value between u max = 1 and u min = 0.

Fig. 3 :

 3 Fig.3: Comparison between the maximal controlled invariant polyhedral set (red line) and the maximal controlled robust invariant polyhedral set (blue line) provided with polyhedral sequences considering state proportional uncertainty w ∈ W = {w ∈ R 2 : w 2 0.06; 0 w 1 ; w 2 0} (black line), computed starting from the state constraints set X. The admitted level of x 2 is equivalent to the admitted level of the normal cells, which is 10 -1 θ N cells. The chemotherapy drugs are constrained to take value between u max = 1 and u min = 0.

Fig. 4 :

 4 Fig. 4: Comparison between the domain of attraction of the deterministic cancer model (red line) and the domain of attraction of the cancer model (blue line) affected by state proportional uncertainty w ∈ W = {w ∈ R 2 : w 2 0.06; 0 w 1 ; w 2 0}. The admitted level of the normal cells is 10 -1 θ N cells. The chemotherapy drugs are constrained to take value between u max = 1 and u min = 0.

TABLE I :

 I Model terms description[START_REF] Afenya | Acute leukemia and chemotherapy: a modeling viewpoint[END_REF].

GIPSA-lab, Grenoble Campus, 11 rue des Mathématiques, BP 46, 38402 Saint Martin d'Hères Cedex, France. {rachid.riah, mirko.fiacchini, mazen.alamir}@gipsa-lab.fr

This work was supported by the French national project CATS-INSERM.