Rachid Riah
email: rachid.riah@gipsa-lab.fr

Mirko Fiacchini
email: mirko.fiacchini@gipsa-lab.fr

Mazen Alamir
email: mazen.alamir@gipsa-lab.fr

Iterative method for estimating the robust domains of attraction of non-linear systems: Application to cancer chemotherapy model with parametric uncertainties

Keywords: Robust domains of attraction, parametric uncertainties, invariance, parameters-dependent CDI system, cancer chemotherapy model

In this paper, we present an iterative procedure method for estimating the robust domains of attraction of non-linear systems. This method is based on the approximation of the uncertain non-linear system with a parameters-dependent Convex Difference Inclusions (CDI) system and the classical iterative methods for linear systems, which are introduced in this paper. A robust one-step operator computing a sequence of convex sets is derived, and the polyhedral case is discussed. An algorithm summarizing the iterative procedure based on the robust one-step operator is given, which is the theoretical contribution of this paper. This method is applied to cancer chemotherapy model considering parametric uncertainties and it is shown that drastic reduction of the robust domain of attraction of the cancer chemotherapy model has happened and this is caused by the parametric uncertainties. It is also proved that the chemotherapy aggressive is not the effective treatment for all the patients.

Introduction

Approximating the region of attraction of non-linear systems is an important task in model analysis and controlled design/evaluation, and several works have been devoted to this issue [START_REF] Chesi | Domain of attraction: analysis and control via SOS programming[END_REF]. These regions of attraction have the property of invariance. The importance of invariant sets in control and analysis of dynamical systems is due to the implicit stability and robustness properties of these regions of the state space. Many results regarding invariance and related topics have been provided in the literature: see, for instance, the notable pioneering contribution [START_REF] Bertsekas | Infinite-time reachability of state-space regions by using feedback control[END_REF], the works [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], concerning the maximal invariant set, and [START_REF] Raković | Invariant approximations of the minimal robust positively invariant set[END_REF] regarding the minimal one. The problem of obtaining invariant sets for discrete-time non-linear systems is dealt with using ellipsoids in [START_REF] Magni | A stabilizing modelbased predictive control algorithm for nonlinear systems[END_REF] and polytopes in [START_REF] Alamo | Convex invariant sets for discrete-time Lur'e systems[END_REF][START_REF] Fiacchini | On the computation of convex robust control invariant sets for nonlinear systems[END_REF]. For a recent monograph on the subject, reader can consult [START_REF] Blanchini | Set-Theoretic Methods in Control[END_REF]. The mathematical theory that may be used to address these issues is the viability theory, see [START_REF] Aubin | Viability theory[END_REF] for more details.

According to our knowledges, there are few contributions for estimating the region of attraction of non-linear systems using set-theoretic methods. An important result in this field is the theoretical and computational methods for a wide class of nonlinear systems and the corresponding iterative procedures to approximate them given in [START_REF] Fiacchini | Convex difference inclusions for systems analysis and design[END_REF][START_REF] Fiacchini | On the computation of convex robust control invariant sets for nonlinear systems[END_REF]. In [START_REF] Fiacchini | Invariant sets computation for convex difference inclusions systems[END_REF], invariant sets computation for Convex Difference Inclusions (CDI) systems are investigated. In [START_REF] Aubin | Differential inclusions: set-valued maps and viability theory[END_REF], CDI systems and set-valued maps are discussed for systems control purposes. In [START_REF] Riah | Invariance-based analysis of cancer chemotherapy[END_REF], the estimation of the region of attraction of non-linear system is based on the approximation of the non-linear system with CDI system as well as the classical set-theoretic methods of linear systems. This technique is motivated by the fact that the use of linear matrix inequalities (LMI) and iterative procedures for linear systems has become very popular due to their advantageous properties and the availability of efficient numerical tools to solve LMI problems.

In numerous control applications the parameters of the dynamical models are considered unknown but belong to known intervals. Approximating the robust region of attraction of such models are important in model analysis and controlled design/evaluation. Based on the results of [START_REF] Fiacchini | Convex difference inclusions for systems analysis and design[END_REF][START_REF] Fiacchini | On the computation of convex robust control invariant sets for nonlinear systems[END_REF][START_REF] Riah | Invariance-based analysis of cancer chemotherapy[END_REF], in this work we present computational method for estimating the robust domain of attraction of non-linear systems. This method is based on the robust one-step operator of the parameters-dependent CDI system that approximates the uncertain non-linear system. This is the theoretical contribution of the paper.

Mathematical models for tumor growth have appeared in the last decades and several scientific researchers have been interested by this topic [START_REF] Swan | Role of optimal control theory in cancer chemotherapy[END_REF][START_REF] Martin | Optimal control drug scheduling of cancer chemotherapy[END_REF][START_REF] De Pillis | Mixed immunotherapy and chemotherapy of tumors: modeling, applications and biological interpretations[END_REF][START_REF] Matveev | Application of optimal control theory to analysis of cancer chemotherapy regimens[END_REF]. These models have been used for control-based 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 [START_REF] Chareyron | Mixed immunotherapy and chemotherapy of tumors: Feedback design and model updating schemes[END_REF][START_REF] Alamir | Robust feedback design for combined therapy of cancer[END_REF][START_REF] Alamir | On probabilistic certification of combined cancer therapies using strongly uncertain models[END_REF]. Some of them are based on the evolution of the different populations of cells [START_REF] Afenya | Acute leukemia and chemotherapy: a modeling viewpoint[END_REF][START_REF] Afenya | Some perspectives on modeling leukemia[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][START_REF] De Pillis | Mixed immunotherapy and chemotherapy of tumors: modeling, applications and biological interpretations[END_REF] and may incorporate the effects of external drugs on the tumor growth, for instance chemotherapy and immunotherapy drugs or their combination.

Several objectives for control-based tumoral chemotherapy design were posed in the literature. In some cases, the main goal was to minimize the amount of chemotherapy drug, with the constraint that the tumor size does not exceed the prescribed level at the end of the treatment period. In other cases, the objective was to minimize the number of tumor cells in a fixed therapy period under some constraints [START_REF] Matveev | Application of optimal control theory to analysis of cancer chemotherapy regimens[END_REF][START_REF] Chareyron | Mixed immunotherapy and chemotherapy of tumors: Feedback design and model updating schemes[END_REF][START_REF] Ledzewicz | Optimal control for combination therapy in cancer[END_REF][START_REF] Alamir | On probabilistic certification of combined cancer therapies using strongly uncertain models[END_REF][START_REF] De Pillis | A mathematical tumor model with immune resistance and drug therapy: an optimal control approach[END_REF]. Minimize the amount of the therapy drug and the number of tumor cells at the end of the treatment period simultaneously is also considered in [START_REF] Alamir | Robust feedback design for combined therapy of cancer[END_REF]. As a healthy condition, maintaining a normal cells population above a given level is often used [START_REF] Matveev | Application of optimal control theory to analysis of cancer chemotherapy regimens[END_REF][START_REF] Murray | Optimal drug regimens in cancer chemotherapy for single drugs that block progression through the cell cycle[END_REF].

In this paper, the model based on the works done in [START_REF] Afenya | Acute leukemia and chemotherapy: a modeling viewpoint[END_REF][START_REF] Afenya | Some perspectives on modeling leukemia[END_REF][START_REF] Matveev | Application of optimal control theory to analysis of cancer chemotherapy regimens[END_REF][START_REF] Swan | Role of optimal control theory in cancer chemotherapy[END_REF]] is adopted to reproduce the behaviour of the normal cells population levels in presence of tumor cells and under the effect of chemotherapy drugs. It was given in form of a system of ordinary differential equations. This model originally was 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 tumor growth. The normal cells represent a part of the innate immune system and can be interpreted as the aggregation of NK cells, CD8+ cells and a circulating lymphocytes (or white blood cells) as in [START_REF] Chareyron | Mixed immunotherapy and chemotherapy of tumors: Feedback design and model updating schemes[END_REF]. In this model, Gompertzian functions are employed to describe the tumor and normal cells growths and some assumptions are introduced in order to maintain sufficient simplicity to admit analysis.

The practical contribution of this paper, whose preliminary version is [START_REF] Riah | Domain of attraction estimation of cancer chemotherapy model affected by state proportional uncertainty[END_REF], consists in the development of numerical tool for cancer chemotherapy model considering parametric uncertainties. According to the literature, this assumption is crucial since in general the chemotherapy cancer models involve many unknown and dynamically varying set of parameters [START_REF] Alamir | Robust feedback design for combined therapy of cancer[END_REF]. The presence of these uncertainties is important since, in practice, the model parameters can vary from patient to patient and are usually not exactly known. On the other hand, uncertainties might make the analysis problem much more complex. The practical contribution is done by applying the computational method that will be developed in the first part of this paper, which is based on the properties of invariant sets of parametersdependent CDI systems.

This numerical tool compute the robust domain of attraction of the uncertain cancer chemotherapy model. It contains all the normal and tumor cells states for which a set of appropriate chemotherapy drug administration profiles exist. These profiles drive the states of this region to the safe region. The safe region is defined to be the set for which the number of tumor cells population is small enough and the normal cells population is higher than the minimal admitted level. A minimal level is imposed as a healthy condition for the patients.

In this paper, it will be shown that considering parameters uncertainties in the cancer chemotherapy models is crucial in order to achieve a good tumor contraction without killing the normal cells. It will be proved using simulations that for some patients the chemotherapy aggressive is not the efficiency chemotherapy drug administration. A comparative study between several drug strategies is provided. The characterization of all the robust domains of attraction for each chemotherapy drug administration profile is given, and the suitable strategy for each patient can be inferred by analyzing the different robust domains of attraction.

The paper is organized as follows: In Section II, the problem statement and the theoretic contribution of the paper, which is a computational method to estimate the robust domain of attraction of non-linear systems, are given. In Section III, the mathematical cancer chemotherapy model is first given then the numerical tool for approximating the robust domains of attraction of this model is developed. Simulations results for cancer chemotherapy analysis are illustrated in Section IV. Section V ends the paper by providing conclusions and giving hints for future investigations.

Notation 1. 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.

Robust domains of attraction of non-linear systems

In this section the problem statement and the theoretic contribution of the paper are stated.

We consider discrete-time uncertain non-linear systems of the form

x + = f P (x, u), ∀x ∈ X, ∀u ∈ U, P ∈ Σ, (1)
where x is the state vector, x + is the successor, u is the control input, and P is the vector of uncertain parameters. The function f P is a parameters-dependent non-linear function. Furthermore, the sets X and U are intended to be, hereafter, respectively the set of state constraints and the set of admissible control inputs, and Σ is the bounding set of the parameters vector P. They are assumed to be known convex regions.

The set of all initial conditions from which the trajectories of (1) converge to the safe region (it can be just a point) regardless the parametric uncertainties is called the robust domain of attraction. In this paper, the main goal is to propose an algorithm to compute this domain using set-theoretic methods.

According to set-theoretic methods, the robust domain of attraction of the uncertain non-linear system (1) is the maximal robust invariant set. Since in set-theoretic methods the convexity property of sets is often considered, this domain is convex. In order to compute this domain an algorithm, which is a part of the contributions of this paper, is based on the following procedures:

1. Approximation of the uncertain non-linear system (1) with a parameters-dependent CDI system; 2. Characterization of the robust domain of attraction of the parameters-dependent CDI system using set-theoretic methods.

The CDI systems are characterized by a particular class of set-valued maps as dynamic functions. In particular, the setvalued map determining a CDI system is such that, given a point in the state space, its image through the map is a convex and compact set, for more details about CDI systems, see [START_REF] Fiacchini | Convex difference inclusions for systems analysis and design[END_REF].

The approximation of a non-linear system with a CDI system is useful in order to apply the set-theoretic methods. In cancer chemotherapy analysis, that will be given in the next section, it will be shown that this algorithm can compute the exact robust domain of attraction of cancer chemotherapy model.

In order to approximate a non-linear system, one defines a CDI system which is characterized by a set-valued map whose graph contains the graph of the function determining the nonlinear system. Indeed, an invariant (contractive) set for the CDI system is also invariant (contractive) for the non-linear system. Thus due to the particular properties of the CDI systems, convexity is preserved by the one-step operator and then the numerical methods proper for linear systems are also valid for CDI systems.

Now in order to approximate the uncertain non-linear system [START_REF] Chesi | Domain of attraction: analysis and control via SOS programming[END_REF], which is defined by the parameters-dependent function f P , it is sufficient to determine a parameters-dependent setvalued map that bounds the parameters-dependent function f P . Therefore a robust invariant (contractive) set for the parametersdependent CDI system, which is defined by a parametersdependent set-valued map, is also invariant (contractive) for the uncertain non-linear system [START_REF] Chesi | Domain of attraction: analysis and control via SOS programming[END_REF]. Now let us introduce some useful tools to deal with convex closed sets and CDI systems. First consider the parametersdependent difference inclusions systems

x + ∈ F P (x, u), (2)
where x ∈ X ⊆ R n is the state, x + is the successor, u ∈ U ⊆ R m is the control input, P ∈ Σ ⊆ R r is the vector of parameters and F P (•, •) is a parameters-dependent set-valued map on R n . F P (x, u) represents a function which relates a set to every point (x, u, P) ∈ R n ×R m ×R r , i.e. F P (x, u) ⊆ R n for all x ∈ R n , u ∈ R m and P ∈ R r . We consider that the parameters-dependent difference inclusions system (2) is used to approximate the uncertain non-linear system (1), i.e. f P (x, u) ∈ F P (x, u), ∀x ∈ X, ∀u ∈ U, and P ∈ Σ.

The parameters-dependent difference inclusions system could be obtained by bounding the dynamics of the uncertain non-linear system (1) and this can be done by bounding the parameters-dependent function f P (•, •).

Hereafter, the second step to compute the robust domain of attraction of the uncertain non-linear system (1) is given. First, an important tool to deal with convex closed sets is the support function, it is defined as follows.

Definition 1. Given a set Ω ⊆ R n , the support function of Ω evaluated at η ∈ R n is given by φ Ω (η) such that

φ Ω (η) = sup x∈Ω η T x.
Geometrically the support function of Ω at η is the signed "distance" of the point of the closure of Ω farthest from the origin, along the direction η. See [START_REF] Rockafellar | Convex Analysis[END_REF] for some properties of support functions. Using the support function is helpful to transform a 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. [31] Given the closed convex set Ω ⊆ R n and the set Γ ⊆ R n then x ∈ Ω if and only if η T x ≤ φ Ω (η) for all η ∈ R n , and Γ ⊆ Ω if and only if φ Γ (η) ≤ φ Ω (η) for all η ∈ R n .
Before giving the definition of a controlled robust invariant set let us introduce the following assumption.

Assumption 1. Assume that the parameters-dependent setvalued map F P determining the system dynamics ( 2) is, such that F P (x, u) for all P ∈ Σ is compact and convex for all (x u) T ∈ R n × R m and for every η ∈ R n and P ∈ Σ the function

F P (x, u, η) : R n × R m × R n -→ R defined as F P (x, u, η) = φ F P (x,u) (η) = sup z∈F P (x,u) η T z, (3)
is convex with respect to (x u) T on X × U.

Thus, in practice, given P ∈ Σ, F P (x, u, η) is the support function of the set F P (x, u) evaluated at η ∈ R n and then

F P (x, u) = z ∈ R n : η T z ≤ F P (x, u, η), ∀η ∈ R n , ( 4
)
and it is convex in (x u) T .

The dynamical system (2) for which Assumption 1 holds are tightly related to the 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]. In fact, given a known parameters vector P ∈ Σ then the system (1) can be approximated by the parameters-dependent CDI system (2), i.e.

f P (x, u) ∈ F P (x, u), ∀x ∈ X, and ∀u ∈ U.

However, when the parameters vector P is unknown and is assumed to belong to the known set Σ the parameters-dependent function f P (•, •) has to be approximated by the union of F P (•, •) for all P ∈ Σ. Therefore, the uncertain non-linear system (1) can be approximated by the following parameters-dependent CDI system

x + ∈ P∈Σ F P (x, u), (5)
where F P (•, •) is a parameters-dependent set-valued map on R n defined above. Now let us introduce the standard definition of the controlled robust invariant set for the generic non-linear system, adapted here to the sets that do not necessarily contain the origin and the parameters-dependent 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 system (1) with P ∈ Σ if for all x ∈ Ω there exists u ∈ U such that f P (x, u) ∈ Ω for all P ∈ Σ. It is controlled robust invariant for the system (5) if for all x ∈ Ω there exists u ∈ U such that F P (x, u) ⊆ Ω for all P ∈ Σ.

For every trajectory starting in a controlled robust invariant set Ω there exists a control input such that this trajectory remains inside Ω regardless of the parametric uncertainties. For the system (5), the controlled robust one-step operator is employed to check the controlled robust invariance of Ω ⊆ X. It is defined hereafter. Definition 3. Consider the closed convex set Ω and assume that Assumption 1 holds for the parameters-dependent setvalued map F P (•, •) determining the dynamic system [START_REF] Raković | Invariant approximations of the minimal robust positively invariant set[END_REF]. The controlled robust one-step operator is defined as follows

Q(Ω, U, X) = x ∈ X | ∃u ∈ U : P∈Σ F P (x, u) ⊆ Ω . (6)
Notice that obtaining the controlled robust one-step operator (6) requires sweeping the parameter vector P in the set Σ which might be computationally very demanding. Hence in order to solve this issue one focuses on the vertices of the bounding set Σ. Finally let us introduce the following proposition.

Proposition 1. Consider the closed convex set Ω and assume that Assumption 1 holds for the parameters-dependent setvalued map F P (•, •) determining the dynamic system [START_REF] Raković | Invariant approximations of the minimal robust positively invariant set[END_REF]. Given a polytope Σ and the set V of its vertices, assume that the support function F P (x, u, η) of the set F P (x, u) is convex with respect to P, then the controlled robust one step-operator given by ( 6) is equivalent to

Q(Ω, U, X) = x ∈ X | ∃u ∈ U : F P (x, u) ⊆ Ω, ∀P ∈ V . (7)
Before giving the proof the following property that is helpful for the proof is introduced, see [START_REF] Bertsekas | Convex Optimization Theory[END_REF] for the demonstration and for other details. Property 2. [START_REF] Bertsekas | Convex Optimization Theory[END_REF] Let C be a closed convex subset of R n that has at least one extreme point. A convex function f : C → R n that attains a maximum over C attains the maximum at some extreme point of C.

Proof 1. Applying Property 1, the inclusion P∈Σ F P (x, u) ⊆ Ω in Eq.( 6) is equivalent to

φ P∈Σ F P (x,u) (η) ≤ φ Ω (η), ∀η ∈ R n .
By using the properties of the support function, see [START_REF] Rockafellar | Convex Analysis[END_REF], the latter condition is equivalent to

φ P∈Σ F P (x,u) (η) = max P∈Σ φ F P (x,u) (η) = max P∈Σ F P (x, u, η) ≤ φ Ω (η), ∀η ∈ R n , (8) since F P (x, u, η) for η ∈ R n is the support function of F P (x, u) defined in Assumption 1.
Similarly, applying Property 1 and using the definition of the support function, the inclusion F P (x, u) ⊆ Ω, ∀P ∈ V in Eq.( 7) is equivalent to

F P (x, u, η) ≤ φ Ω (η), ∀P ∈ V, ∀η ∈ R n ,
where V is the set of vertices of the set Σ, such condition is also equivalent to

max P∈V F P (x, u, η) ≤ φ Ω (η), ∀η ∈ R n . ( 9
)
Now one proves that the inequality in Eq.( 8) is equivalent to the inequality in Eq. [START_REF] Blanchini | Set-Theoretic Methods in Control[END_REF]. Necessity is due to V ⊆ Σ since Σ is a polytope. To prove sufficiency, one have to prove that

max P∈Σ F P (x, u, η) ≤ max P∈V F P (x, u, η), ∀η ∈ R n . ( 10
)
Since [START_REF] Aubin | Viability theory[END_REF] would mean that max

P∈V F P (x, u, η) ≤ φ Ω (η) implies that max P∈Σ F P (x, u, η) ≤ φ Ω (η) ∀η ∈ R n ,
resonating by contradiction will be applied and this by supposing that there exists η ∈ R n such that max

P∈Σ F P (x, u, η) > max P∈V F P (x, u, η).
By hypothesis, for a given (x, u, η) ∈ X × U × R n , the function F P (x, u, η) : Σ → R n is convex with respect to P ∈ Σ, and Σ is compact. Thus the function F P (x, u, η) attains his maximum with respect to P over the set Σ.

Therefore, according to the Property 2 the function F P (x, u, η) attains his maximum at some extreme point of Σ. Hence this is a contradiction and the supposition is false then the condition given in ( 10) is true and the sufficiency is proved. Consequently, the condition ( 8) is equivalent to ( 9) and the robust one-step operator given in Definition 3 is equivalent to that given in Proposition 1.

The controlled robust one-step operator associates to every set Ω the set of points for which there exists a set of admissible controls u ∈ U such that these points will be mapped inside Ω through P∈V F P (x, u).

The one-step operator will be used to compute an increasing sequences of nested controlled robust invariant sets, as for linear or non-linear systems. In fact, Algorithm 1, standard for generating increasing controlled robust invariant approximations of the robust domain of attraction, see [START_REF] Blanchini | Set-Theoretic Methods in Control[END_REF], can be also applied in this context.

Algorithm 1 Increasing sequences of controlled robust invariant sets for (5)

Input: Initial convex closed controlled robust invariant set Ω 0 ⊆ X. 1: for k ∈ N N do 2: Compute Ω k+1 = Q(Ω k , U, X) ∩ X 3: if Ω k+1 = Ω k then Ω max = Ω k return 4:
end if 5: end for Thus Ω k , k ∈ N N are the controlled robust invariant sets that converge to the robust domain of attraction of the parametersdependent CDI system [START_REF] Raković | Invariant approximations of the minimal robust positively invariant set[END_REF], and hence of all the non-linear systems approximated by [START_REF] Raković | Invariant approximations of the minimal robust positively invariant set[END_REF]. Since in this paper the robust domain of attraction is needed to be calculated in polyhedral form, the following proposition is introduced, analogous to [START_REF] Fiacchini | Convex difference inclusions for systems analysis and design[END_REF], which is functional for this purpose. Proposition 2. Given the parameters-dependent set-valued map F P (x, u) determining the system dynamics (5), its support function F P (x, u, η) which is assumed to be convex with respect to P, and the state constraints set X; given also a polytope Ω = x ∈ R n : Hx ≤ h with H ∈ R n h ×n ; given also a polytope Σ and V its vertices, the controlled robust one-step operator is given as follows

Q(Ω, U, X) = x ∈ X | ∃u ∈ U : F P (x, u, H i T ) ≤ h i , ∀i ∈ N n h , ∀P ∈ V . (11)
Notice that for every polytope Ω and admissible control polyhedral set U the set Q(Ω, U, X) is closed and convex, as proved in [START_REF] Fiacchini | Convex difference inclusions for systems analysis and design[END_REF]. Moreover, if F P (•, •, η) are piecewise affine functions of (x u) T then Q(Ω, U, X) is a polyhedron, which is the intersection of a finite number of halfspaces for every polytope (or polyhedron) Ω.

The Algorithm 1 is suitable to compute the robust domain of attraction of the parameters-dependent CDI system (5) which is an approximation of the robust domain of attraction of the uncertain non-linear system [START_REF] Chesi | Domain of attraction: analysis and control via SOS programming[END_REF]. A meaningful application of the method developed in this section will be given in the rest of the paper.

In practice, these domains are the stability/save regions of the non-linear systems. In the second part of this paper, this algorithm will be applied for a cancer chemotherapy analysis.

Application to cancer chemotherapy model

Mathematical cancer model under chemotherapy drugs 3.1.1. State variables

The mathematical model considered 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][START_REF] Afenya | Some perspectives on modeling leukemia[END_REF] to describe the behaviour of the cells populations in the presence of tumor and under the chemotherapy treatment is used in this paper. This model involves the following cells populations:

• L, tumor cells population;

• N, normal cells population. These cells are part of the innate immune system and in [START_REF] Afenya | Acute leukemia and chemotherapy: a modeling viewpoint[END_REF][START_REF] Afenya | Some perspectives on modeling leukemia[END_REF] it represents the neutrophil cells population.

Modelling assumptions

The model given below is based on the following assumptions taken from [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][START_REF] Afenya | Some perspectives on modeling leukemia[END_REF]:

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 Gompertzian growth; 3. The chemotherapy drug kills both tumor and normal cells; 4. The chemotherapy drug affects tumor cells more than normal cells; 5. The tumor cells exhibit a negative and inhibiting effect on the growth and development of the normal cells; 6. The drug spreads instantaneously within the body.

Dynamic model of a cancer under chemotherapy drugs

The assumptions given above lead to the following model, for cells populations dynamics, in presence of tumor cells and under chemotherapy drug administration:

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 ,
where α and β are the tumor and normal growth rates respectively, θ L and θ N are the greatest size of tumor cells population and the normal size of normal cells population respectively. The terms αL ln θ L L and βN ln θ N N represent the Gompertizian growth of both tumor and normal cells populations respectively, and the term -Ξ(L)N represents the negative effect of the tumor on the normal cells.

The variable u stands for the concentration of the chemotherapy drug, and Π 1 (u) and Π 2 (u) represent loss functions for both cells populations due to chemotherapy effects.

According to the works [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][START_REF] Afenya | Some perspectives on modeling leukemia[END_REF], we adopt the following approximations:

• 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);

• The loss function Π 1 (u) is considered to be 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;

• The loss function Π 2 (u) is also considered linear, i.e, Π 2 (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);

• The parameters of the model are uncertain and bounded within known intervals. This leads to a polytopic parametric uncertain model.

Thus the resulting model becomes:

L = αL ln θ L L -kuL, L(0) = L 0 , Ṅ = βN ln θ N N -γLN -luN, N(0) = N 0 P LN = α β γ θ N θ L k l T ∈ Σ LN , (12)
where P LN is a vector of the parameters of the uncertain model, and Σ LN is the polytopic bounding set. Since Σ LN is a polytope, then it can be represented as the convex hull of its vertices, i.e. Notice that the non-linear model [START_REF] Fiacchini | Invariant sets computation for convex difference inclusions systems[END_REF] involves two state variables L, N, one manipulated variable u, and n P = 7 uncertain parameters. The consistent nominal values for these parameters are inferred from [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], 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 1.

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, x 1 (0) = x 10 , ẋ2 = -βx 2 + γθ L e -x 1 + lu, x 2 (0) = x 20 , P = α β γ θ L k l T ∈ Σ, (13)
where P, Σ are derived directly from ( 12), x 10 and x 20 are the initial conditions. The discrete-time system modelling the cancer evolution is obtained by sampling the continuous-time system [START_REF] Aubin | Differential inclusions: set-valued maps and viability theory[END_REF], 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 = f P1 (x, u), x 1 (0) = x 10 , x + 2 = (1 -T s β)x 2 + T s γθ L e -x 1 + T s lu = f P2 (x, u), x 2 (0) = x 20 , P = α β γ θ L k l T ∈ Σ. (14)
We notice that this model is more suitable to apply the iterative methods for invariant set computation.

Control problem and constraints

In cancer chemotherapy the main objective is to contract the tumor cells population L while maintaining the normal cells population N above a prescribed level [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 number of the normal cells population N is considered as a measure of the patient health [START_REF] Matveev | Application of optimal control theory to analysis of cancer chemotherapy regimens[END_REF].

In this section, the objective is to apply the algorithm developed in the first part of the paper in order to develop a numerical tool for cancer chemotherapy analysis. This numerical tool leads to determine all the initial tumor and normal states for which there exist appropriate drug injection profiles. These profiles must lead to a substantial regression of the tumor size while avoiding that the health measure reaches dangerous values for the patient. Thus, the substantial regression of the tumor size must be ensured without knowing the values of the parameters of the cancer chemotherapy model.

Explicitly the numerical tool provides the set of the initial points (L, N), such that the dynamics of the model (12) converge to the safe region, when an admissible chemotherapy drug profile is applied. The safe region is defined as the set in the state space L, N, where L is small enough and N ≥ N min , with N min denoting the minimal admitted value of the normal cells population.

Similarly in the domain of coordinates x 1 and x 2 , this numerical tool provides the set of the initial points (x 1 , x 2 ), such that the dynamics of the model ( 14) can be driven to the safe region, by applying admissible chemotherapy drug profiles. The safe region in the state space x 1 , x 2 is 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 are represented by the robust domain of attraction of the cancer chemotherapy model [START_REF] Riah | Invariance-based analysis of cancer chemotherapy[END_REF]. Its calculation is substantially based on the properties of robust positively invariant sets and convex inclusions that are introduced in the first part of this paper.

Robust domain of attraction of cancer chemotherapy model

The robust domain of attraction of the cancer chemotherapy model ( 12) can be computed by following these procedures:

1. Approximating the discrete-time non-linear cancer model ( 14) by a parameters-dependent CDI system (5); 2. Applying Algorithm 1 for the parameters-dependent CDI system ( 5) that approximates the model ( 14); The principal task is to compute the controlled robust one-step operator [START_REF] Fiacchini | Convex difference inclusions for systems analysis and design[END_REF]. The result of this procedure is the robust domain of attraction in the state space (x 1 , x 2 ); 3. Computing the robust domain of attraction in the state space (L, N) by using change of variables.

Controlled robust one-step operator for cancer chemotherapy model

Here the controlled robust one-step operator for the cancer chemotherapy model in the state space (x 1 , x 2 ) ( 14) is determined by using the approximations and considerations given in Section 2.

Firstly, it is necessary to determine the bounding functions F P (x, u, η) for all η ∈ R n that is convex with respect to (x u) T . Its associated parameters-dependent set-valued maps F P (x, u), that defines the parameters-dependent CDI system (5) must satisfy the following inclusion

f P (x, u) ∈ F P (x, u), ∀(x u) T ∈ X × U,
with f P (x, u) defines the cancer model ( 14), and F P (x, u) is defined by the bounding functions F P (x, u, η) as it is given by the equation ( 4).

Thus, according to Property 1, the functions F P (x, u, η) must satisfy the following inequality

η T f P (x, u) ≤ F P (x, u, η), ∀η ∈ R n .
Moreover, since the non-linearity in [START_REF] Riah | Invariance-based analysis of cancer chemotherapy[END_REF] involves only f P2 (x, u) then the bounding functions of η T f P (x, u) are related only to η 2 , since the dynamics of x 1 is linear. Thus a possible choice is

• if η 2 ≥ 0: F P (x, u, η) = η T f P (x, u),
since F P (x, u, η) is already convex in this case;

• if η 2 < 0:

F P (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 .

Remark 1. Notice that, if the non-linear function Ξ(L) in (3.1.3) can be taken strictly increasing, convex and Ξ(0) = 0 as in [START_REF] Matveev | Application of optimal control theory to analysis of cancer chemotherapy regimens[END_REF], the convexity is preserved and in this case, the convex bounding function of η T f P (x, u) is also easily obtainable.

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 ,
is the 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 . Figure 1 illustrates this approximation.

Remark 2. 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] Swan | Role of optimal control theory in cancer chemotherapy[END_REF], since the function ( 15) is already piecewise affine function. Then, we obtain

F P (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 P (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 P (x, u, η).

Notice that the bounding functions F P (x, u, η) are convex upper bounds of η T f P (x, u), for all (x u) T ∈ X × U and for all η ∈ R n . Once convex upper bounds of η T f P (x, u) are determined, the controlled robust one-step operator is given by using the Proposition 2.

Since in [START_REF] Martin | Optimal control drug scheduling of cancer chemotherapy[END_REF] we have the term γθ L then this function is not convex with respect to the vector of model parameters P, and the Proposition 2 can not be applied directly.

In order to derive a representation of F P (x, u, η) with convex functions, we use the following change of variable δ = ∈ Σ, where Σ is determined from Σ by using the change of variable given above. Now the function given by ( 16) with this change of variable becomes

F P(x, u, η)=                      η 1 (1 -T s α)x 1 + T s ku + η 2 (1 -T s β)x 2 +T s δ max i∈N q {c i x 1 + d i } + T s lu if η 2 ≥ 0, η 1 (1 -T s α)x 1 + T s ku + η 2 (1 -T s β)x 2 +T s δ(ax 1 + b) + T s lu if η 2 < 0, (17)
and is convex with respect to P, since it is the sum of convex functions. Now let us assume that the polytope Ω = x ∈ R n : Hx ≤ h with H ∈ R n h ×n , the polyhedron X = x ∈ R n : x 1 ≥ 0; 0 ≤ x 2 ≤ x 2max , with x 2max represents the maximal admitted value of x 2 and the polytope U = u ∈ R m : u min ≤ u ≤ u max } with u max and u min are respectively the maximal and the minimal value of the chemotherapy drug, are given for the cancer chemotherapy analysis.

Then considering that the polytope Σ is also given, the controlled robust one-step operator, introduced in Proposition 2, is defined by using the upper bounds of η T f P (x, u) given by [START_REF] De Pillis | Mixed immunotherapy and chemotherapy of tumors: modeling, applications and biological interpretations[END_REF]. Finally, Algorithm 1 is applied in order to compute the maximal robust invariant polyhedral set for cancer chemotherapy model since F P(x, u, η) is piecewise affine function.

Remark 3. By applying Algorithm 1, we have found that the sequence of polytopes Ω k = {x ∈ R n : H k x ≤ h k } generated by the controlled robust one-step operator have never a facet determined by H k i with H k i,2 < 0 except the case related to the trivial constraints x 2 ≤ 0 which 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 robust 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 2. 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 non-linear cancer model ( 14) by a parameters-dependent CDI system (5) can be reduced increasing the precision as indicated in Remark 2.

Thus according to Remark 3, the piecewise affine bound function of η T f P (x, u) to be employed is that given for η 2 ≥ 0. Therefore the related controlled robust one-step operator for the cancer model ( 14) is given as follow:

Q(Ω, U, X) = x ∈ X : ∃ u ∈ U : ∀i ∈ N n h , ∀ j ∈ N q , ∀ P ∈ V, 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 , ( 18
)
where V is the set of vertices of Σ. Then Q(Ω, U, X) maps polytopes in polytopes and Algorithm 1 generates a sequence of polytopes.

Therefore this one-step operator can be used in Algorithm 1 to compute the robust domain of attraction of the uncertain cancer chemotherapy model ( 14). This domain is given in state space (x 1 , x 2 ).

Tumor-normal cells robust domain of attraction

Once the robust domain of attraction Ω max for the cancer chemotherapy model ( 14) has been obtained by applying Algorithm 1 with the one-step operator [START_REF] Matveev | Application of optimal control theory to analysis of cancer chemotherapy regimens[END_REF], the tumor-normal cells robust domain of attraction of the cancer chemotherapy model ( 12) can be determined using change of variables. Now we assume that Ω max = x ∈ R n : Hx ≥ h with H ∈ R n h ×n . Thus by using the change of variables x 1 = ln θ L L and

x 2 = ln θ N N
the robust domain of attraction in the state space (L, N) is given by the following equation

Υ = (L N) T ∈ R 2 : N ≥ max i∈N n h θ N |H i,2 | θ L |H i,1 | e -h i 1 |H i,2 | L |H i,1 | |H i,2 | . (19)
Since the Eq.( 19) depends on the parameters θ N and θ L , and these parameters are fixed but unknown, i.e. θ N ∈ [θ N θ N ] and θ L ∈ [θ L θ L ], then the tumor-normal cells robust domain of attraction defined by Eq.( 19) becomes

Ῡ = (L N) T ∈ R 2 : N ≥ max i∈N n h θ N |H i,2 | θ L |H i,1 | e -h i 1 |H i,2 | L |H i,1 | |H i,2 | , ∀θ N ∈ [θ N θ N ], ∀θ L ∈ [θ L θ L ] . (20)
In Eq. [START_REF] Ledzewicz | Optimal control for combination therapy in cancer[END_REF], we have to maximize the term

θ N |H i,2 | θ L |H i,1 | e -h i 1 |H i,2 | L |H i,1 | |H i,2 | , ∀i ∈ N n h , ∀θ N ∈ [θ N θ N ] and ∀θ L ∈ [θ L θ L ]. Hence it is equivalent to consider only the terms θ N |H i,2 | θ L |H i,1 | e -h i 1 |H i,2 | L |H i,1 | |H i,2 |
for all i ∈ N n h and Eq.( 20) is equivalent to

Ῡ = (L N) T ∈ R 2 : N ≥ max i∈N n h θ N |H i,2 | θ L |H i,1 | e -h i 1 |H i,2 | L |H i,1 | |H i,2 | . (21)
Consequently the set Ῡ is the maximal non-convex robust domain of attraction of the cancer chemotherapy model [START_REF] Fiacchini | Invariant sets computation for convex difference inclusions systems[END_REF] in the state space (L, N).

Simulation results

In this section, the method for the computation of the robust domain of attraction that is developed in Section 2 is applied to analyse the cancer chemotherapy model given in Section 3. The parameters of this model are considered unknown but assumed to belong into known intervals. According to our knowledge, this assumption is consistent and maybe it is the suitable manner to address the problem of variability of the parameters of cancer models.

For simulation the chemotherapy drug profile is constrained to take value between u min = 0 and u max = 1. The minimal admitted level of the normal cells population is N min = 0.1θ N Cells, which implies that x 2max = 2.30 by using change of variable defined in the previous sections. Since the exact values of the model parameters are not available in the literature the value of N min is also taken approximately as others parameters described and noted in Section 2.

Once the set of state constraints X and the initial robust invariant set Ω 0 are selected, Algorithm 1 is fully automatic and no parameters are needed to be selected.

Robust domains of attraction versus uncertainty degree of model parameters

Applying Algorithm 1 with the controlled robust one-step operator [START_REF] Matveev | Application of optimal control theory to analysis of cancer chemotherapy regimens[END_REF], by starting with an initial controlled robust invariant set Ω 0 , a sequence of controlled robust invariant sets is computed. This sequence of sets converge to the controlled robust invariant set of the cancer chemotherapy model [START_REF] Riah | Invariance-based analysis of cancer chemotherapy[END_REF], in the state space (x 1 , x 2 ).

For the model considered in this paper the set Ω 0 = {x ∈ X : x 1 ≥ a} with a big enough, for which the invariance condition Ω 0 ⊆ Q(Ω 0 , X, U) holds, is used as an initial controlled robust invariant set. The controlled robust one-step operator [START_REF] Ledzewicz | Optimal control for combination therapy in cancer[END_REF] is to be used to check the robust invariance of sets.

In order to illustrate the use of the method developed in the paper, let us assume that the parameters of the cancer chemotherapy model [START_REF] Fiacchini | Invariant sets computation for convex difference inclusions systems[END_REF] can have variations of 30% around their nominal values.

Figure 2 shows a sequence of nested controlled robust invariant sets that converge to the maximal controlled robust invariant set. This sequence of sets is depicted in black line. Notice that Ω k+1 is defined as the set that contains all the states for which there exist a set of admissible chemotherapy drug such that these states can be mapped in Ω k in one step regardless of the parametric uncertainties. In the same figure, the maximal controlled invariant set of the cancer chemotherapy model ( 14) computed taking the nominal values of the parameters is depicted in blue line.

By applying [START_REF] Chareyron | Mixed immunotherapy and chemotherapy of tumors: Feedback design and model updating schemes[END_REF] to the sets given in Figure 2, one gets the tumor-normal cells robust domains of attraction of the cancer chemotherapy model [START_REF] Fiacchini | Invariant sets computation for convex difference inclusions systems[END_REF]. These domains are depicted in Figure 3. In black line, the robust domain of attraction of the cancer chemotherapy model ( 12) is showed. This domain contains all the tumor and normal cells states that can be driven to the set of healthy states considering that the parameters of the cancer chemotherapy model are uncertain. In blue line the domain of attraction of the cancer chemotherapy model ( 12) is shown for comparison purpose.

According to simulation results given in this section, we remark drastic reduction of the domain of attraction and we conclude that this reduction is caused by model parameters uncertainties. Therefore, in order to achieve a successful chemotherapy treatment, the identification of the real values of the model parameters or the assumption that the model parameters are unknown but belong to given known intervals are crucial.

The robust domain of attraction, with arbitrary precision, is the exact domain of all the initial tumor and normal cells states for which there exist appropriate administration chemotherapy profiles. These administration profiles lead to tumor cells contraction without violating the healthy condition of the patients. Consequently, given a patient state one could infer on the existence of successful chemotherapy treatment. Furthermore, the amount of drugs that should be delivered and the therapy length can also be calculated following the method of this paper. 5 show the controlled robust invariant sets of the cancer chemotherapy model [START_REF] Riah | Invariance-based analysis of cancer chemotherapy[END_REF] in the state space (x 1 , x 2 ), and the robust domains of attraction of the cancer chemotherapy model [START_REF] Fiacchini | Invariant sets computation for convex difference inclusions systems[END_REF], computed by considering parameters variations of 0%, 10%, 20%, 30%, 40% and 50% around their nominal values. From these figures, we conclude that the more important the parameters variations around their nominal 14) in the state space (x 1 , x 2 ). The admitted level of x 2 is equivalent to the admitted level of the normal cells population, which is 0.1θ N Cells. The chemotherapy drugs are constrained to take values between u max = 1 and u min = 0 4.2. Analysis of different chemotherapy drug administration profiles For a given chemotherapy drug administration profile the related robust one-step operator can be determined by composing [START_REF] Matveev | Application of optimal control theory to analysis of cancer chemotherapy regimens[END_REF]. For instance, the robust one-step operator for the profile given by one sampling period of full drug injection and one period of null chemotherapy drug is given by the following composition Q Q Ω, {0}, X , {u max }, X .

In order to analyse and compare the robust domains of attraction of appropriate drug administration strategies, three profiles of chemotherapy drugs are taken into account, which are: 12). The admitted level of the normal cells population is 0.1θ N Cells. The chemotherapy drugs are constrained to take values between u max = 1 and u min = 0

0 0.5 1
• chemotherapy drug administration profile 1: Heavy doses of drugs are applied all the time (chemotherapy aggressive);

• chemotherapy drug administration profile 2: Heavy doses of drugs are applied for 12 days then no doses are applied for 12 days;

• chemotherapy drug administration profile 3: Heavy doses of drugs are applied for 4 days then no doses are applied for 4 days.

Model parameters variations of 30% around their nominal values are considered in the cancer chemotherapy model [START_REF] Fiacchini | Invariant sets computation for convex difference inclusions systems[END_REF]. Now applying Algorithm 1, Figure 6 shows the robust domains of attraction of the cancer chemotherapy model [START_REF] Fiacchini | Invariant sets computation for convex difference inclusions systems[END_REF] related to the different chemotherapy drug administration profiles given below. By analyzing the robust domains, one can notice that there are tumor and normal cells populations states that cannot be cured by the chemotherapy drug administration profile 3, whereas they can be healed by applying the chemotherapy drug administration profile 2 and there are tumor and normal cells populations states that can not be cured by chemotherapy drug administration profile 2, however they can be healed by applying the chemotherapy drug administration profile 1. Notice that both profiles 2 and 3 have the same drug delivery rate, i.e. 0.5u max /day. Consequently, it is more beneficial to attack the tumor cells population for long period then release, than for small period then release.

In Figure 6, the robust domain of attraction computed using the controlled robust one-step operator [START_REF] Matveev | Application of optimal control theory to analysis of cancer chemotherapy regimens[END_REF] is depicted in solid line. The chemotherapy drug administration profile is constrained to take value between u max and u min . We remark that this domain is the maximal domain of attraction of the cancer chemotherapy model [START_REF] Fiacchini | Invariant sets computation for convex difference inclusions systems[END_REF]. Consequently, the aggressive chemotherapy is not the recommended chemotherapy treatment for all the patients, and an appropriate treatment for each patient is more beneficial.

Conclusion

In this paper we develop an iterative procedure method to approximate the robust domains of attraction of non-linear systems.This procedure is based on set-theoretic methods and parameters-dependent CDI systems which is used to approximate the uncertain non-linear systems. The good approximation allows us to find the maximal robust domain of attraction as it is shown in the application. In the second part of this paper, we apply the method developed for uncertain nonlinear systems to cancer chemotherapy model considering parametric uncertainties and saturation constraints which limit the chemotherapy drug injections. Thus we have developed numerical tool that leads to compute the robust domain of attraction of cancer chemotherapy model. This domain contains all the initial tumor and normal cells states for which there exists a set of admissible chemotherapy drug administration profiles. It is shown that in order to achieve a successful chemotherapy treatment, either the identification of the real values of the model parameters or the awareness during the treatment synthesis that the parameters of the cancer model are not well known, is crucial.

In future work, alternative modelling frameworks could be employed and other methods for estimating the robust domains of attraction of non-linear systems could follow the work done in this paper.

Σ

  LN = Co P LN (1), ..., P LN (Z) , where Co denotes a convex hull, P LN (z), z = 1, . . . , Z are vertices of the polytope Σ LN . Notice that, for any P LN ∈ Σ LN , there exist non-negative coefficients λ z satisfying Z z=1 λ z = 1, P LN = Z z=1 λ z P LN (z).

1 Figure 1 :

 11 Figure 1: Bounds of the exponential function

Figure 2 :Figure 4

 24 Figure 2: Comparison between the maximal controlled invariant set (blue line) and the maximal controlled robust invariant set computed by considering parameters variations of 30% around their nominal values (black line). The admitted level of x 2 is equivalent to the admitted level of the normal cells, which is 0.1θ N Cells. The chemotherapy drugs are constrained to take values between u max = 1 and u min = 0 Figure 4 and Figure5show the controlled robust invariant sets of the cancer chemotherapy model[START_REF] Riah | Invariance-based analysis of cancer chemotherapy[END_REF] in the state space (x 1 , x 2 ), and the robust domains of attraction of the cancer chemotherapy model[START_REF] Fiacchini | Invariant sets computation for convex difference inclusions systems[END_REF], computed by considering parameters variations of 0%, 10%, 20%, 30%, 40% and 50% around their nominal values. From these figures, we conclude that the more important the parameters variations around their nominal

Figure 3 :

 3 Figure 3: Comparison between the domain of attraction (blue line) and the robust domain of attraction (black line), computed considering parameters variations of 30% around their nominal values, of the cancer model (12). The admitted level of the normal cells is 0.1θ N Cells. The chemotherapy drugs are constrained to take values between u max = 1 and u min = 0

Figure 4 :

 4 Figure 4: Controlled robust invariant sets, computed by considering parameters variations of 0%, 10%, 20%, 30%, 40% and 50% around their nominal values, of the cancer chemotherapy model (14) in the state space (x 1 , x 2 ). The admitted level of x 2 is equivalent to the admitted level of the normal cells population, which is 0.1θ N Cells. The chemotherapy drugs are constrained to take values between u max = 1 and u min = 0

Figure 6 :

 6 Figure 6: Robust domains of attraction of the cancer chemotherapy model (12) related to different chemotherapy drug administration profile considering parameters variations of 30% around their nominal values. The admitted level of the normal cells population is 0.1θ N cells

Table 1 :

 1 Parameters of the cancer chemotherapy model[START_REF] Fiacchini | Invariant sets computation for convex difference inclusions systems[END_REF]

	param value	param value
	α	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

Acknowledgements

This work has been supported by the INSERM (Institut National de la Santé et de la Recherche Médical) projects CATS (Cancer Assisted Therapeutic Strategies).


M Fiacchini
email: mirko.fiacchini@gipsa-lab.fr

G Millérioux
email: gilles.millerioux@univ-lorraine.fr

Dead-beat stabilizability of discrete-time switched linear systems: algorithms and applications Mirko Fiacchini, Gilles Millérioux. Abstract-This paper deals with dead-beat stabilizability of autonomous discrete-time switched linear systems. Based on a constructive necessary and sufficient condition for dead-beat stabilizability, we propose two algorithms. The first one is concerned with the problem of testing dead-beat stabilizability and computing the shorter stabilizing mode sequence, whenever it exists. The other one implements a method to construct a switched system whose shorter dead-beat stabilizing sequence has a prescribed length. Then, we present numerical assessments and possible applications.

Index Terms-Switched systems, dead-beat stabilizability

I. INTRODUCTION

Switching systems are dynamical systems for which the state dynamics vary between different operating modes according to a switching sequence [START_REF] Liberzon | Switching in Systems and Control, ser. Systems and Control: Foundations and Applications[END_REF]. Such systems are found in many practical and theoretical domains. For example they are relevant models in networked control systems [START_REF] Alur | Compositional modeling and analysis of multi-hop control networks[END_REF], [START_REF] Jungers | Feedback stabilization of dynamical systems with switched delays[END_REF], in congestion control for computer networks [START_REF] Shorten | A positive systems model of TCPlike congestion control: asymptotic results[END_REF], in viral mitigation [START_REF] Hernandez-Vargas | Optimal and MPC switching strategies for mitigating viral mutation and escape[END_REF], as abstractions of more complex hybrid systems [START_REF] Liberzon | Basic problems in stability and design of switched systems[END_REF], and other fields (see e.g. [START_REF] Jungers | The joint spectral radius[END_REF]- [START_REF] Shorten | Stability criteria for switched and hybrid systems[END_REF] and references therein).

This papers addresses the problem of finite-time stabilizability of discrete-time switched linear systems. The objective is two-fold. First, we investigate the problem of checking whether there exists a finite sequence of switches so that the product of the dynamical matrices gives the null matrix. Secondly, we are concerned with building a dead-beat stabilizable system with a sequence of switches of prescribed length. The difficulty in obtaining some relevant results stems from the fact that the question of deciding the stabilizability of a switching system is known to be hard in general, see [START_REF] Jungers | The joint spectral radius[END_REF].

First, let us recall some existing results of closely related problems. Dead-beat stabilizability has been addressed for systems with control. In the survey [START_REF] Sun | Analysis and synthesis of switched linear control systems[END_REF], conditions for controllability and observability of switched linear systems based on the concept of A-invariant are provided. Special results concerning the discrete-time case given by [START_REF] Conner | State deadbeat response and observability in multi-modal systems[END_REF], [START_REF]The structure of the controllable set for multimodal systems[END_REF] are recalled. These papers provide bounds on the minimal length of a controlling sequence for completely controllable systems. The results are given for switched systems with nonsingular transition matrices, defined as reversible in [START_REF] Sun | Analysis and synthesis of switched linear control systems[END_REF].

For switched linear autonomous systems, exponential stabilizability has been addressed in the survey [START_REF] Lin | Stability and stabilizability of switched linear systems: a survey of recent results[END_REF]. Therein, it is recalled that the discrete-time counterpart of stabilizability conditions for continuous-time autonomous switched linear systems is not straightforward. Only a sufficient condition, proposed by the authors themselves, is presented. Another criterion, probably less conservative, is given in [START_REF] Geromel | Stability and stabilization of discretetime switched systems[END_REF] but is also sufficient. General necessary and sufficient conditions for exponential stabilizability are provided in [START_REF] Fiacchini | Necessary and sufficient condition for stabilizability of discrete-time linear switched systems: A set-theory approach[END_REF] and their computation-oriented relaxations are provided and analyzed in [START_REF] Fiacchini | On the stabilizability of discrete-time switched linear systems: Novel conditions and comparisons[END_REF] but the conditions are not appropriate for dead-beat stabilizability.

Furthermore, we should note that dead-beat stabilizability is closely related to the problem of mortality of a set of matrices. Indeed, a set of matrices is mortal if the null matrix can be expressed as the product of finite length of matrices. This problem is discussed in the papers [START_REF] Paterson | Unsolvability in 3 × 3 matrices[END_REF]- [START_REF] Bournez | The mortality problem for matrices of low dimensions[END_REF]. It is proved that the decidability on the mortality of a set of matrices is unsolvable except for some particular cases whose mortality problem results NP-complete, e.g. pairs of integer matrices. In the paper [START_REF] Paterson | Unsolvability in 3 × 3 matrices[END_REF], it is proved that the problem is unsolvable for finite sets of 3 × 3 matrices over the integers.

Finally, another closely related issue is the dead-beat stability for autonomous switched linear discrete-time systems. Unlike dead-beat stabilizability, the problems boils down to checking whether every sufficiently long sequence of switches leads to the null matrix. In [START_REF] Parriaux | Nilpotent semigroups for the characterization of flat outputs of switched linear and LPV discrete-time systems[END_REF], dead-beat stability has been used to characterize flatness of discrete-time controlled switched linear systems. The issue has been tackled with the notion of nilpotent semigroups. Dead-beat stability has been examined in [START_REF] Philippe | Deciding the boundedness and dead-beat stability of constrained switching systems[END_REF] as well, in which constraints on the switching sequence are considered.

This literature overview shows that the problem of characterizing dead-beat stabilizability for discrete-time switched linear systems is an open problem. A constructive necessary and sufficient condition has been published in the preliminary paper [START_REF] Fiacchini | Dead-beat stabilizability of autonomous switched linear discrete-time systems[END_REF]. In the present paper, we go further. It is proved that a system is dead-beat stabilizable if and only if there exists one switching sequence that stabilizes the whole state space. Next, such a condition is used for defining a complete algorithm to test if a given switched system is dead-beat stabilizable. This algorithm considerably outperforms the brute force exhaustive search approach. Furthermore, the condition is used to define a second algorithm for constructing dead-beat stabilizable switched systems. Finally, numerical evaluations and possible applications, namely flatness and cryptography, are presented.

The paper is organized as follows. In Section II, the problems of dead-beat stability and stabilizability are presented. Condition for dead-beat stability and dead-beat stabilizability along with existing algorithms to test them are recalled. The equivalence between dead-beat stabilizability and the existence of one switching sequence that stabilizes the whole state space is proved. In Section III, we provide two algorithms for respectively testing and building dead-beat stabilizable systems. Section IV presents a numerical example and two possible practical applications of dead-beat stabilizable systems like flatness and cryptography. Some concluding remarks are finally given in Section V.

Notation: Given n ∈ N, define N n = {j ∈ N : 1 ≤ j ≤ n}. The set of q switching modes is I = N q and the related matrices forms a finite collection A ⊆ R n×n , whose i-th element is denoted with A i , i.e. A = {A i } i∈I , with A i ∈ R n×n for all i ∈ I. All the possible sequences of modes of length N is

I N = N j=1 I, with I N = ∅ if N = 0. The length of a sequence σ is denoted by |σ| and |σ| = N if σ ∈ I N . For all L, M ∈ N such that 0 ≤ L ≤ M , define I [L:M ] = M N =L I N .
Given σ ∈ I N and i, j ∈ N N with i ≤ j, σ i is the i-th element of σ and σ [i:j] the subsequence of σ starting at σ i and terminating at σ j . The subsequence σ [i:j] is the empty sequence if i > j. Given σ, δ sequences of modes, (σ, δ) is their concatenation. Given σ ∈ I N , define

A σ = N j=1 A σj = A σ N • • • A σ1 and n j=m A σj = I if m >
n where I stands for the identity matrix. Given a ∈ R, a is the smallest integer greater than or equal to a.

II. NECESSARY AND SUFFICIENT CONDITION FOR DEAD-BEAT STABILIZABILITY

A. Problem statement

Consider the discrete-time switched linear system

x k+1 = A σ k x k , (1)
where x k ∈ R n is the state at time k ∈ N and σ : N → I is the switching law. With slight abuse of notation, we employ the index of σ, as in σ k for instance, to denote both the time realization of the switching law at instant k and the k-th element of the sequence of modes σ ∈ I N , its meaning being determined by the context. Before addressing the issue of dead-beat stabilizability, let us first recall a closely related problem, that is the robust deadbeat (or finite-time) stability. This property stipulates that for any sequences of switches, every state reaches the origin in finite steps. More formally,

∃N ∈ N s.t. x N = 0 ∀σ ∈ I N ∀x 0 ∈ R n which is equivalent to ∃N ∈ N s.t. A σ = A σ N • • • A σ1 = 0 ∀σ ∈ I N .
The former condition is equivalent to state that the set A generates a nilpotent semigroup whose usual definition is recalled below.

Definition 1: A semigroup S is a set together with an associative internal law. A semigroup S with an absorbing element 0 is said to be nilpotent if there exists a positive integer t ∈ N such that the internal law applied to any t elements of S is always equal to 0. The smallest integer t is called the class of nilpotency of S.

If S is a set of matrices the associative internal law is the matrix multiplication and the absorbing element is the null matrix.

Theorem 1 (Levitsky's Theorem): Any semigroup of nilpotent matrices can be triangularized. This theorem is computationally useful since it gives a way of testing nilpotency, and thus robust dead-beat stability, by means of triangularization routines that are polynomially complex with respect to the dimension n of the matrices. Clearly, a necessary condition for the set {A i } i∈Nq to generate a nilpotent semigroup is that every A i is nilpotent for i ∈ N q .

Similarly to dead-beat stability for which a necessary and sufficient condition and an algorithm to test the condition exist, we aim at providing a condition and an algorithm for characterizing dead-beat stabilizability. We first give its formal definition.

Definition 2: The system (1) is dead-beat stabilizable if for every x ∈ R n , there exist N (x) ∈ N and a finite sequence

σ(x) ∈ I N (x) such that A σ(x) x = 0.
Thus, dead-beat stabilizability implies the existence of a function, associating to every state, a stabilizing switching sequence, that is a closed-loop control, if the switching is considered as the control input. We prove below that a system is dead-beat stabilizable if and only if there exists one switching sequence that stabilizes the whole state space. This would mean that the system is stabilizable through an open-loop switching control. Recall that the existence of a switching sequence exponentially stabilizing the whole space, a condition often referred to as consistent stabilizability, is only sufficient for exponential stabilizability, see [START_REF] Sun | Stabilizability and insensitivity of switched linear systems[END_REF], [START_REF] Sun | Stability Theory of Switched Dynamical Systems[END_REF]. On the other hand, as a peculiarity of dead-beat stabilizability, the following theorem applies.

Theorem 2: The system (1) is dead-beat stabilizable if and only if there exist N ∈ N and γ ∈ I N such that A γ = 0.

Proof: Sufficiency is trivial. We must prove necessity, i.e. that if the system is dead-beat stabilizable, then there exists a finite mode sequence γ stabilizing the whole state space. We proceed by contradiction. Suppose that there exists a deadbeat stabilizable system for which A σ = 0 for all σ ∈ I N and every N ∈ N and recall that A σ = 0 implies that ker(A σ ) < n. Stabilizability implies that x ∈ ker(A σ(x) ) for every x ∈ R n and then

N ∈N σ∈I N ker(A σ ) = R n ,
that is absurd since the space R n cannot be expressed as the union of countably many subspaces of dimensions smaller than n. Indeed, the union of a (finite or infinite) sequence of n-dimensional closed sets is an n-dimensional closed set, see [START_REF] Kuratowski | Introduction to Set Theory and Topology[END_REF].

Hence, from Theorem 2, to characterize dead-beat stabilizability, it suffices to derive a necessary and sufficient condition for the existence of N ∈ N and a finite sequence γ ∈ I N such that for any initial condition x 0 ∈ R n , the state at time N reaches the origin. That is

∃N ∈ N, ∃γ ∈ I N s.t. A γ = 0. (2)
Remark 1: According to the considerations made in the introduction on the notion of mortality, the existence of such a N cannot, in general, be performed in finite time, resulting in an unsolvable problem.

Clearly, robust stability, associated to the notion of nilpotent semigroup, implies dead-beat stabilizability, related to mortality of a set of matrices.

Remark 2: Beyond analysis purposes, let us notice that deadbeat stabilizability can be interesting for control perspectives, in particular flatness-based control for switched linear systems in the form

x k+1 = A σ k x k + B σ k u k , y k = C σ k x k + D σ k u k , (3)
where

u k ∈ R p is the input, y k ∈ R m is the output and A σ k , B σ k , C σ k , D σ k are state space matrices of appropriate dimension.
Recall that a non-autonomous dynamical system (assumed to be square) is flat if there exist an output y k , referred to as flat output, and an integer k 0 , such that all system variables can be expressed, for k ≥ k 0 , as a function of the flat output and a finite number of its backward and/or forward shifts. It turns out that flatness of (3) can be characterized, see [START_REF] Millérioux | Flatness of switched linear discrete-time systems[END_REF], by considering the product of matrices

K i=1 P σ k+i-1
, where every P σ k+i-1 belongs to the finite set of matrices defining the left inverse dynamics of (3) and is expressed in terms of the state space matrices of (3). The switching rule σ is a function of a finite number r of modes σ k , . . . , σ k-r where r is the inherent delay of the system (coinciding with the relative degree for a SISO system). It has been proved that the system (3) is flat if and only if there exists a finite integer K such that

K i=1 P σ k+i-1 = 0 (4)
holds for any arbitrary sequences σ . However, we can consider a more flexible situation when both the modes σ k and the control input u k can be chosen at time k. In such a case, first, we can be interested in an algorithm that finds a sequence of switches which guarantees (4). Then, appropriately applying such a mode sequence ensures that the the system (3) is flat with flat output y k and a flatness-based control u k can be designed separately.

In next subsection, we recall from [START_REF] Fiacchini | Dead-beat stabilizability of autonomous switched linear discrete-time systems[END_REF] the necessary and sufficient conditions to determine N and γ, whenever they exist.

B. Necessary and sufficient condition

Denote by

I s = {i ∈ I : det A i = 0}, I ns = {i ∈ I : det A i = 0}, (5)
the sets of singular and nonsingular matrices in I and by q s and q ns the number of elements of I s and I ns , respectively. Clearly, one has q s + q ns = q where q is the number of modes of the system (1). The following proposition will be instrumental for the sequel.

Proposition 1 ( [21]): The system (1) is dead-beat stabilizable if and only if there exist m ∈ N n finite mode sequences σ j ∈ I n j with n j ∈ N, where j ∈ N m , such that

m j=1 dim im j-1 k=1 A σ k ∩ ker (A σ j ) = n, (6)
and with σ j n j ∈ I s for every j ∈ N m .

Proposition 1 asserts that the system (1) is dead-beat stabilizable if and only if there exists a set of m switching sequences σ j of finite length n j , with j ∈ N m and 1 ≤ m ≤ n, whose last element is related to a singular matrix, i.e. σ j n j ∈ I s , and such that the intersection of the kernel of A σ j and the image of the product j-1 k=1 A σ k of matrices has a dimension strictly greater than zero. Moreover, the dimension of those intersections is equal to n but those intersections do not necessarily span R n . Besides, as shown by Corollary 1 in [START_REF] Fiacchini | Dead-beat stabilizability of autonomous switched linear discrete-time systems[END_REF], if the dead-beat stabilizability condition is satisfied by the set of sequences σ j with j ∈ N m , it also holds if the prefix of σ 1 involving nonsingular matrices is removed. All in all, every stabilizing sequence candidate starts and terminates with a singular matrix. The reader can find in Section IV-A an illustrative example of matrices, images and kernels dimensions of a possible sequence γ for the mortal set of 6 matrices of dimension n = 3 composed by 2 singular matrix and 4 nonsingular ones given in [START_REF] Geromel | Stability and stabilization of discretetime switched systems[END_REF], see Table I.

On one hand, given a finite set of matrices corresponding to the modes of (1), we aim at an approach which is based on the necessary and sufficient condition to find the shorter stabilizing sequence γ that leads to [START_REF] Alur | Compositional modeling and analysis of multi-hop control networks[END_REF]. Clearly, such an approach should perform better than the exhaustive search, that is the brute-force test of all the possible switching sequences that leads to [START_REF] Alur | Compositional modeling and analysis of multi-hop control networks[END_REF]. On the other hand, we are concerned with proposing a method to build a dead-beat stabilizable system with a prescribed length switching sequence. Both issues are addressed in the following section.

III. ALGORITHMS

A. Algorithm to test dead-beat stabilizability

From the explanations given after Proposition 1, the expected algorithm can restrict the search to sequences γ = (σ 1 , . . . , σ m ) ∈ I N such that A γ = 0 is composed of m ∈ N n subsequences σ j terminating with an element of I s and the first subsequence σ 1 is composed of only one singular matrix. Clearly, there might be an infinite number of such sequences. Thus, we propose an algorithm that finds all the stabilizing sequences whose subsequences lengths are bounded by h + 1 with h ∈ N. As a consequence, the algorithm will find the minimal length stabilizing sequence γ, provided that it does not involve subsequences longer than h+1. Algorithm 1 given below meets those requirements.

Algorithm 1 Define the dead-beat stabilizing sequences γ of the system (1) with subsequences length ≤ h + 1. Input: matrices A = {A i } i∈I , sets I and I s , parameter h.

1: C ← ∅, List of candidates 2: M ← ∅;
List of stabilizing sequences 3: for s ∈ I s do insert s in C

Candidates σ 1 4: end for 5: while C is not empty do

6: extract σ from C Candidate σ = (σ 1 , • • • , σ j ) 7:
for s ∈ I s , δ ∈ I [0: h] do 8:

if dim ker (A s A δ A σ ) > dim ker (A σ ) and θ ← (σ, δ, s)

11: if dim ker (A θ ) = n then 12:
insert θ in M θ is stabilizing end for 18: end while Output: M List of stabilizing sequences Some comments on the main steps of Algorithm 1 follow. The structures C and M contain the sequences candidates to dead-beat stabilize (1) and the stabilizing ones, respectively. At Line 3, all the singular matrices are introduced in C, as the first subsequence candidates. At every successive iterative step (starting at Line 5), a candidate σ is considered and extracted from C. The loop starting at Line 7 aims at searching all the subsequence (δ, s) of length h + 1 at most. At iteration p with p ∈ N m , the sequences σ p = (δ, s) must be such that s ∈ I s and

dim im p-1 k=1 A σ k ∩ ker(A σ p ) > 0 (7)
which ensures that every term involved in the sum ( 6) is strictly greater that zero and thus, that the dimension reaches n for a sufficient number of subsequences.

In practice and for the sake of numerical efficiency, the search is performed according to the equivalent test given at Line 8 which is based on the following lemma.

Lemma 1: ( [26]) For every pair of matrices A ∈ R m×n and B ∈ R n×p , we have that dim ker(AB) = dim ker(B) + dim im(B) ∩ ker(A) . Indeed, from Lemma 1, it is clear that ( 7) is equivalent to

dim ker p k=1 A σ k > dim ker p-1 k=1 A σ k . (8)
Hence, at each iteration p, the subsequence σ p = (δ, s) makes drop down the dimension of the image of the product of matrices associated to the new candidate sequence (σ 1 , . . . , σ p ) or equivalently, it makes increase the kernel. The obtained sequence is stored in M if the dimension of the kernel reaches n or in C otherwise. Note that the test at Line 9 aims at disregarding subsequences δ (which may involve singular matrices) leading to dim ker (A δ A σ ) < dim ker (A σ ) and would cause redundancy. Remark 3: Algorithm 1 computes all the mortal sequences with subsequences length r = h + 1 at most and thus can be used to determine the shorter mortal sequence γ by increasing the parameter h. Let us compare the computational complexity of Algorithm 1 with respect to the brute force enumeration. Recall that the number of sequences of q modes, with positive length l ∈ N at most, is given by l i=1 q i = (q l+1 -1)/(q -1) -1. Consider first the case of balanced subsequences σ j , that is such that σ j have length r for i > 1. Thus, γ has length l = (m-1)r+1 and Algorithm 1 has to inspect (m-1)q s (q r -1)/(q-1)+q s sequences, while the exhaustive search needs to test all the (q (m-1)r+2 -1)/(q -1) -1 sequences of lengths l at most. For the general non-balanced case, denote with l ∈ N the length of the shortest mortal sequence and r ∈ N the length of its longer subsequence. Then (l -1)/(m -1) ≤ r ≤ l -m + 1 which means that Algorithm 1 finds γ testing between (m-1)q s (q (l-1)/(m-1) +1 -1)/(q-1)+q s and (m-1)q s (q l-m+2 -1)/(q-1)+q s sequences. Summarizing, the key idea is that Algorithm 1 performs an exhaustive search on m-1 subsequences σ j instead of on the whole γ = (σ 1 , . . . , σ m ), where m ≤ n.

B. Algorithm to build a dead-beat stabilizable system

This section is devoted to the construction of dead-beat stabilizable switched linear systems whose shorter sequence γ such that A γ = 0 has a prescribed length. According to Proposition 1, it is recalled that such a sequence is composed of m subsequences σ j with j ∈ N m . The first subsequence reduces to one element (|σ 1 | = 1) and the corresponding matrix is singular. We deliberately impose the length of the other subsequences such that |σ j | = r j with r j ∈ N, for j ≥ 1. Hence, the length of the expected sequence is 1 + p j=2 r p . The prescribed number of singular and nonsingular matrices are q s and q ns , respectively. We aim at an iterative algorithm which consists in computing successively the subsequences σ j with j > 1.

At Step p ∈ N with 2 ≤ p ≤ n and given σ j with j ∈ N p-1 obtained at Step p -1, the algorithm must compute a subsequence σ p = (δ p , s p ) of given length r p such that:

• the last element A sp is one of the singular matrices I s ;

• the condition

dim ker A sp A δp p-1 k=1 A σ k = dim ker p-1 k=1 A σ k + 1 (9) holds.
Actually, the last condition is a particular case of [START_REF] Lin | Stability and stabilizability of switched linear systems: a survey of recent results[END_REF]. Indeed, it is recalled that, according to Proposition 1, the dimension of the kernel between two consecutive subsequences σ p-1 and σ p increases. Here, we impose that the kernel dimension increases by one between iteration p -1 and p. As a result, at Step p, dim ker p-1 k=1 A σ k = p -1 and the total number m of subsequences is equal to n.

After those considerations, we are able to detail the constructive procedure which is summarized in Algorithm 2. The initialization (Line 4-13) consists in randomly computing the q s singular matrices and q ns -n + 1 nonsingular ones. In particular the singular matrices are given by A s = T -1 s Λ s T s with:

(i) Λ s a matrix with one eigenvalue 0 and the other ones randomly generated but non-null. Without loss of generality, we impose the following structure of Λ s

Λ s =      0 0 • • • 0 0 . . . 0 Λs     
with Λs ∈ R n-1×n-1 a nonsingular matrix. (ii) T s a randomly generated invertible matrix. The columns of T -1 s are the eigenvectors of A s , with the first one related to the null eigenvalue.

The first iteration of Algorithm 2 (Line 15) consists in choosing the only matrix composing σ 1 among the singular ones. Then, the algorithm must build a sequence of subsequences σ p (Line 16-24) causing the kernel dimension to increase of one, see [START_REF] Shorten | Stability criteria for switched and hybrid systems[END_REF]. To this end, consider the following lemma.

Lemma 2:

( [26]) Given A ∈ R m×n and B ∈ R n×p then dim ker (AX) = dim (im (B) ∩ ker(A))
holds, where X is a basis matrix of im (B).

Denoting with X p-1 ∈ R n×n-p+1 a basis matrix of p-1 j=1 A σ j , condition ( 9) is equivalent to dim ker (A σ p X p-1 ) = 1 [START_REF] Sun | Analysis and synthesis of switched linear control systems[END_REF] from Lemmas 1 and 2. This being the case, the computation at iteration p is detailed hereafter.

B p = A -1 αp T -1 sp R p C -1 p C p = A βp X p-1 (A βp X p-1 ) ⊥ (11)
with

R p =      1 0 • • • 0 0 . . . 0 Rp     
where Rp ∈ R n-1×n-1 p is an arbitrary nonsingular matrix and where V ⊥ is a basis of the subspace orthogonal to the one spanned by the columns of V , implying the nonsingularity of C p . The matrix B p is such that

A sp A αp B p A βp X p-1 = T -1 sp Λ sp R p C -1 p A βp X p-1 = T -1 sp Λ sp R p I (:,n-p+1) (12
) where I (:,n-p+1) is the matrix given by the first n -p + 1 columns of the identity matrix (since M -1 M (:,m) = I (:,m) for all nonsingular M and positive m ∈ N). Therefore, since

Λ sp R p =      0 0 • • • 0 0 . . . 0 Λsp Rp      then the kernel of A sp A αp B p A βp X p-1 is {y ∈ R n-p+1 : y i = 0 ∀i ≥ 2}.
Hence, the dimension of the kernel of A sp A αp B p A βp X p-1 increases by 1; (d) define the new mode i p such that A ip = B p . The matrix A ip is introduced in the matrix set of stabilizable system. Define σ p = (β p , i p , α p , s p ) (Line 20-23) Remark 4: The sequences α p and β p implicitly determine the position of the matrix B p within the the p-th subsequence. As particular cases one could either consider an empty α p and |β p | = r p -2, leading to put B p just before A sp , or an empty β p and |α p | = r p -2 yielding B p to be the first matrix related to the p-th subsequence. Note moreover that A αp and A βp are composed by only nonsingular matrices. Also singular ones could be considered, but this could lead to stabilizing sequences whose length is smaller than the desired one, i.e. 1 + p j=2 r p .

Algorithm 2 Build a dead-beat stabilizable system. Input: q s and q ns cardinalities of I and I s , subsequences lengths r p . σ p ← (β p , i p , α p , s p ) 24: end for Output: A Matrix set of the stabilizable system IV. APPLICATIONS A numerical example and a potential application to cryptography of the results presented in this paper are illustrated in this section.

A. Numerical example

This numerical example illustrates the construction of a dead-beat stabilizable system by using Algorithm 2 and then, the search of the shorter mortal sequence γ by using Algorithm 1.

First, by using Algorithm 2, we construct a mortal set of matrices in R 3×3 for which the shorter sequence γ satisfying A γ = 0 has length 11. In particular, the first subsequence has length 1 and the second and third ones have length 5, i.e. r 1 = 1 and r 2 = r 3 = 5. Moreover, we set the number of singular matrices to be q s = 2 and that of nonsingular ones to be q ns = 4. The resulting matrices are (2.1003, -0.1957, 0), (0.3802 + 0.7026i, 0.3802-0.7026i, 0), (-0.5181+0.3747i, -0.5181-0.3747i, 0.5174),

A1 A3 A5 A2 A4 A6 =       -11.
(1.1848, -2.2680, -1.1378), (-1.0175 + 2.0238i, -1.0175 -2.0238i, 0.2518) and (0.2256, 6.9414 + 7.8336i, 6.9414 -7.8336i), respectively. The sequence of modes γ is [START_REF] Alur | Compositional modeling and analysis of multi-hop control networks[END_REF][START_REF] Hernandez-Vargas | Optimal and MPC switching strategies for mitigating viral mutation and escape[END_REF][START_REF] Shorten | A positive systems model of TCPlike congestion control: asymptotic results[END_REF][START_REF] Jungers | Feedback stabilization of dynamical systems with switched delays[END_REF][START_REF] Jungers | Feedback stabilization of dynamical systems with switched delays[END_REF][START_REF] Alur | Compositional modeling and analysis of multi-hop control networks[END_REF][START_REF] Hernandez-Vargas | Optimal and MPC switching strategies for mitigating viral mutation and escape[END_REF][START_REF] Jungers | Feedback stabilization of dynamical systems with switched delays[END_REF][START_REF] Shorten | A positive systems model of TCPlike congestion control: asymptotic results[END_REF][START_REF] Liberzon | Basic problems in stability and design of switched systems[END_REF][START_REF] Liberzon | Switching in Systems and Control, ser. Systems and Control: Foundations and Applications[END_REF], for which the subsequences, image and kernel dimensions are given in Table I.

Aγ A 1 • A 6 •A 4 •A 3 •A 5 • A 2 • A 3 •A 3 •A 4 •A 5 • A 2 dim im(Aσ) 0 1 1 1 1 1 2 2 2 2 2 dim ker(Aσ) 3 2 2 2 2 2 1 1 1 1 1 A σ j A σ 3 A σ 2 A σ 1
The time plot of x i k (i = 1, 2, 3) for 10 different initial conditions (randomly generated in the set x ∞ ≤ 1) are depicted in Figure 1. To illustrate the benefit of using Algorithm 1, the sequence γ has been searched by applying the brute force approach and Algorithm 1 for comparison, using an Intel® Core™ i7-6600U CPU @ 2.60GHz × 4 processor with 16GB of RAM. The brute force approach needed 817s, that is more than 13 minutes, to find the solution. By comparison, Algorithm 1 found the solution in 0.105s. Additional examples highlighting the efficiency of this algorithm can be found in the paper [START_REF] Fiacchini | Dead-beat stabilizability of autonomous switched linear discrete-time systems[END_REF].

B. Cryptography

We illustrate here the potential relevance of dead-beat stabilizability and the use of Algorithm 2 in the context of cryptography. Actually, it is shown that dead-beat stabilizability may allow to break down a barrier in the design of so-called statistical Self-Synchronizing Stream Ciphers, statistical SSSC for short. They are an extension of the class of deterministic Self-Synchronizing Stream Ciphers, often merely named Self-Synchronizing Stream Ciphers. Roughly speaking, SSSC (deterministic or statistical) are based on automata that are dynamical systems operating on finite fields and that must be designed to deliver sequences of symbols, named keystreams, of high complexity from a statistical point of view. Then, those sequences are used to scramble information to be safely transmitted.

A well-admitted approach for the design of SSSC has been first suggested in [START_REF] Maurer | New approaches to the design of self-synchronizing stream cipher[END_REF]. An alphabet denoted by A is considered. It is a finite set of basic elements named symbols. At the ciphering side, the automaton is governed by

z k+1 = g(z k , c k ), w k = h(z k ), (14)
where z k ∈ A is the internal state, g is the next-state transition function, w k ∈ A is the keystream. The ciphertext (cryptogram) is calculated by c k = e(m k , w k ) where e is the ciphering function and m k the plaintext to be encrypted. At the deciphering side, the automaton is governed by

ẑk+1 = g(ẑ k , c k ), ŵk = h(ẑ k ), (15)
where ẑk ∈ A is the internal state, g and h are identical to the cipher functions and ŵk ∈ A is the keystream. The recovered plaintext symbol is calculated by mk = d(c k , ŵk ) where d is the deciphering function fulfilling mk = m k if ŵk = w k .

If the functions g and h are such that there exist a function and an integer M verifying

z k = (c k-1 , . . . , c k-M ), (16)
the automata ( 14) and ( 15) are said to have a finite input memory and M is called the delay of memorization. If such a function exists, it means that after a transient time equal to M , the state z of the cipher does no longer depend on the initial state and so does the state ẑk of the decipher. Actually, the states z k and ẑk coincide each other for any k ≥ M . Hence, synchronization between the cipher and the decipher is automatically achieved and the decryption succeeds. It is the reason why such a cipher is called Self-Synchronizing Stream Cipher. Two classes of automata are defined according to the delay of synchronization:

• Deterministic: The delay of memorization is bounded by the constant M and a priori fixed. • Statistical: The bound of the delay of memorization is not constant but is a random variable with respect to the sequence of ciphertexts or the initial state vector. Automata with a statistical delay of memorization have never been deeply explored so far. We can mention [START_REF] Geraldy | Optimized selfsynchronizing mode of operation[END_REF]- [START_REF] Heys | An analysis of the statistical self-synchronization of stream ciphers[END_REF] for exceptions. Till now, the usual method to obtain statistical Self-Synchronizing Stream Ciphers is to switch between two automata with and without finite input memory property. This hybrid architecture is equipped with a supervisor which is scanning on-line and at both ends (cipher and decipher) a sequence of cryptograms c of prescribed length. This sequence is compared with a cryptogram of reference (called sync pattern) of the same length. Then, when they coincide, the supervisor of both the cipher and the decipher switches from the automaton without input memory to the one with finite input memory, guaranteeing self-synchronization in finite time after the supervisor switch and so proper decryption.

Here, we show that dead-beat stabilizability can be interesting to propose design methodologies leading to statistical self-synchronizing architectures without supervisors. Obtaining lighter architectures can be a central issue in embedded systems. As a clue to tackle the problem, we propose nonlinear automata in the form

z k+1 = Q σ k z k + R σ k c k (17)
where Q σ k belongs to a finite set of matrices and σ k is a switching rule which depends in a nonlinear way on a finite sequence of past ciphertexts c k , . . . , c k-s with s ∈ N. In other words, the switching function is of the form

σ k = ϕ(c k , c k-1 , • • • , c k-s )
where ϕ is a surjective nonlinear function which selects the mode σ k at time k from the knowledge of s past ciphertext symbols. In cryptography, such nonlinear functions are commonly used and are called S-boxes. The equation of the decipher is identical to (17) but involves ẑk instead of z k . The property of finite input memory with a statistical delay of memorization is obtained whenever there exists at least a switching sequence of length K so that K i=1 Q σ k+i-1 = 0. In other words, the auxiliary system defined as x k+1 = Q σ k x k is dead-beat stabilizable. The quantity M is a random variable since it directly depends on the time before a stabilizing switching sequence is expected to appear. The probability law can be simply obtained whenever the ciphertext symbols are uniformly distributed, which is a common feature of the cipher.

We illustrate such a design through a simple example with n = 4 and q = 10 matrices Q i , with i ∈ I. To operate on the field of two elements {0, 1}, we adapted the Algorithm 2 by replacing the matrix product with the product modulo 2. We obtained the following mortal set of matrices

Q1 Q6 Q2 Q7 Q3 Q8 Q4 Q9 Q5 Q10 =    
0 1 0 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 0 0 1 0 0 0 0 0 0 1 0 0 1 1 0 1 0 0 1 0 0 0 0 0 1 0 0 0 1 1 0 1 0 0 0 0 1 1 0 0 1 0 1 0 0 0 0 0 1 0 0 1 0 0 0 0 0 1 1 0 0 1 0 1 1 0 1 0 1 0 0 1 0 0 0 0 0 1 0 1 0 0 0 0 1 0 1 0 1 0 0 0 1 1 1 1 1 1 0 1 0 0 0 0 1 0 0 1 0 1 1 0 0 0 1 1 1 1 0 0 0 1 0 1 0 0 0 1 1 1 1 1 0 0 1 0 0 1 1 0 0 0     for which the shorter sequences γ generating the null matrix have length K = 10. Note that Q 1 is the only singular matrix.

Then, N r = 1000 random switching sequences σ, with σ k uniformly distributed, have been generated. Recall that σ k = ϕ(c k , c k-1 , • • • , c k-s ). Denoting with n i the instant at which the i-th sequence generates the null matrix, Figure 2 (top) shows an histogram depicting the relative frequencies of switching sequences that ensure synchronization within a given range of time. Figure 2 (bottom) shows the obtained approximation of the cumulative probability distribution of synchronization with respect to the time k. As expected, the probability tends towards 1 as k increases. It is clear that the shape of the probability depends on the length of the shorter sequences generating the null matrix and the number of modes. Those parameters should be included in the overall design parameters of the cipher and the decipher but this matter is out of the scope of the paper.

V. CONCLUSION

In this paper, we have considered the problem of characterizing the dead-beat stabilizability for discrete-time switched linear systems. Based on a constructive necessary and sufficient condition, we have proposed an algorithm to check the dead-beat stabilizability of a system and to compute the shorter stabilizing sequence, whenever it exists. It has been shown through a complexity analysis and numerical examples, that the algorithm significantly outperforms an exhaustive search approach. Besides, from the condition, we have derived a second algorithm to build a dead-beat stabilizable systems. Possible applications have been discussed to highlight the interest of the results.

9 :

 9 dim ker (A δ A σ ) = dim ker (A σ ) then 10:

  (a) compute the random sequences α p and β p of modes related to nonsingular matrices and such that |α p |+|β p | = r p -2 (Line 17). Note that A αp and A βp are nonsingular; (b) select s p ∈ I s , then A sp singular (Line 18); (c) define (Line 19)

3 Fig. 1 .

 31 Fig. 1. Evolution in time of x i k (i = 1, 2, 3)

Fig. 2 .

 2 Fig. 2. Relative frequency histogram (top) and cumulative frequency (bottom) of the synchronization instants in function of time, for Nr = 1000 uniformly generated mode sequences.

TABLE I SEQUENCE γ ,

 Iγ SUBSEQUENCES, IMAGE AND KERNEL DIMENSIONS