Andrii Mironchenko

Christophe Prieur

Fabian Wirth

Local stabilization of an unstable parabolic equation via saturated controls

Keywords: PDE control, reaction-diffusion equation, saturated control, stabilization, attraction region

We derive a saturated feedback control, which locally stabilizes a linear reaction-diffusion equation. In contrast to most other works on this topic, we do not assume the Lyapunov stability of the uncontrolled system and consider general unstable systems. Using Lyapunov methods, we provide estimates for the region of attraction for the closed-loop system, given in terms of linear and bilinear matrix inequalities. We show that our results can be used with distributed as well as scalar boundary control, and with different types of saturations. The efficiency of the proposed method is demonstrated by means of numerical simulations.

I. INTRODUCTION

In applications of control technology, physical inputs (like force, torque, thrust, stroke, etc.) are often limited in size [START_REF] Bernstein | A chronological bibliography on saturating actuators[END_REF]. If such input limitations are neglected, this may result in undesirable oscillations of the closed-loop system, lack of global stabilizability and in a dramatic reduction of the region of attraction of the closed-loop system, see e.g. [START_REF] Tarbouriech | Stability and Stabilization of Linear Systems with Saturating Actuators[END_REF], [START_REF] Zaccarian | Modern Anti-windup Synthesis: Control Augmentation for Actuator Saturation[END_REF], [START_REF] Saberi | Control of Linear Systems with Regulation and Input Constraints[END_REF] for an introduction to the nonlinear behavior induced by input limitations.

This leads to the problem of (local or global) stabilization of control systems with inputs of a norm not exceeding a prescribed value.

In this paper, we study the stabilizability of a onedimensional linear unstable reaction-diffusion equation (also called heat equation) using a saturated control. Many results exist in the literature for the control of this class of equations, using either bounded or unbounded control operators, with or without input delays. More specifically, in [START_REF] Krstic | Control of an unstable reaction-diffusion PDE with long input delay[END_REF] a backstepping approach is applied to design a boundary delayed feedback control for a heat equation (see [START_REF] Krstic | Boundary Control of PDEs: A Course on Backstepping Designs[END_REF] for further results using the same design methods). In [START_REF] Fridman | Exponential stability of linear distributed parameter systems with time-varying delays[END_REF] a stable heat partial differential equation (PDE) is controlled by means of a delayed bounded linear control operator (see also [START_REF] Solomon | Stability and passivity analysis of semilinear diffusion PDEs with time-delays[END_REF] for the semilinear case). For the case of unbounded control operators, see [START_REF] Nicaise | Stability of the heat and wave equations with boundary time-varying delays[END_REF], [START_REF] Prieur | Feedback stabilization of a 1-D linear reactiondiffusion equation with delay boundary control[END_REF] for the computation of delayed control to stabilize the reaction-diffusion equation.

A. Mironchenko is with Faculty of Computer Science and Mathematics, University of Passau, 94030 Passau, Germany. Email: andrii.mironchenko@uni-passau.de. Corresponding author.

C. Prieur is with Univ. Grenoble Alpes, CNRS, Grenoble INP, GIPSA-lab, 38000 Grenoble, France. Email: Christophe.Prieur@gipsa-lab.fr F. Wirth is with Faculty of Computer Science and Mathematics, University of Passau, 94030 Passau, Germany. Email: fabian.(lastname)@uni-passau.de.

A. Mironchenko has been supported by DFG, grant Nr. MI 1886/2-1. C. Prieur has been partially supported by MIAI @ Grenoble Alpes, (ANR-19-P3IA-0003). A preliminary version of this paper was presented at the 11-th IFAC Symposium on Nonlinear Control Systems (NOLCOS 2019), see [START_REF] Mironchenko | Design of saturated controls for an unstable parabolic PDE[END_REF].

To analyze the effect of input saturations and to design saturated controllers for infinite-dimensional systems, different results and techniques are available. One of the first papers in this field is [START_REF] Slemrod | Feedback stabilization of a linear control system in Hilbert space with an a priori bounded control[END_REF], where compact and bounded control operators are considered. Another notable early reference is [START_REF] Lasiecka | Strong stability of elastic control systems with dissipative saturating feedback[END_REF] where an observability assumption is stated for the study of PDEs with constrained controllers. In these results, the control inputs are functions from a certain input space, and the saturation map has to be understood as a limitation on the norm of the input in this space. A more physically relevant type of saturation is given by pointwise saturations, which limit the values of the control input at each point by a prescribed value.

Pointwise saturations are more complex and require further developments, some of which we investigate in the present paper. Other types of saturation functions can also be useful in practice, see [START_REF] Marx | Stability results for infinitedimensional linear control systems subject to saturations[END_REF] for a discussion.

The definition of the saturation map used in the present paper is aligned to that used in [START_REF] Prieur | Wave equation with cone-bounded control laws[END_REF] where Lyapunov methods are shown to be useful for the stability analysis of wave equations subject to saturated inputs. See also [START_REF] Marx | Stability results for infinitedimensional linear control systems subject to saturations[END_REF], [START_REF] Marx | Cone-bounded feedback laws for mdissipative operators on Hilbert spaces[END_REF], where systems in Hilbert spaces with applications to the Kortewegde Vries equation have been addressed.

Our main result is the derivation of a saturated feedback control, which locally stabilizes the unstable heat equation, and for which estimates of the region of attraction for the closed-loop system can be provided. These estimates are formulated in terms of matrix inequalities which can be efficiently solved numerically. We emphasize that, in contrast to the works [START_REF] Slemrod | Feedback stabilization of a linear control system in Hilbert space with an a priori bounded control[END_REF], [START_REF] Lasiecka | Strong stability of elastic control systems with dissipative saturating feedback[END_REF], [START_REF] Prieur | Wave equation with cone-bounded control laws[END_REF], where the stress is on the global stabilization of marginally stable hyperbolic systems (which have infinitely many eigenvalues on the imaginary axis), in our paper we consider parabolic systems, which have finitely many exponentially unstable modes. Therefore the control of the system with a saturating controller only provides local asymptotic stability, paralleling what is known for finitedimensional systems (see in particular [START_REF] Teel | Global stabilization and restricted tracking for multiple integrators with bounded controls[END_REF]).

We would like to mention the related research on model predictive control of parabolic PDEs under constraints on the state and the input [START_REF] Dubljevic | Predictive control of parabolic pdes with state and control constraints[END_REF], and on Lyapunov-based control of parabolic systems by controls of a bounded magnitude [START_REF] El-Farra | Analysis and control of parabolic PDE systems with input constraints[END_REF]. Note however, that in these papers the problem of estimating the region of attraction has not been studied, which is the key objective of this work. In [START_REF] Kang | Boundary control of delayed ODE-heat cascade under actuator saturation[END_REF] stabilizing backstepping-based boundary controllers for coupled heat-ODE systems with timevarying state delays in the presence of actuator saturation are developed, and in [START_REF] Kang | Boundary constrained control of delayed nonlinear Schrödinger equation[END_REF] boundary stabilization of a nonlinear Schrödinger equation with state delay and bounded internal disturbance is performed. Both in [START_REF] Kang | Boundary control of delayed ODE-heat cascade under actuator saturation[END_REF], [START_REF] Kang | Boundary constrained control of delayed nonlinear Schrödinger equation[END_REF] estimates of the region of attraction of the closed-loop system are provided.

Our approach is based on the spectral decomposition of the open-loop dynamics of a heat equation into the finitedimensional unstable part and infinite-dimensional stable dynamics, which is a classical tool used in PDE control [START_REF] Russell | Controllability and Stabilizability Theory for Linear Partial Differential Equations: Recent Progress and Open Questions[END_REF], [START_REF] Barbu | Internal stabilization of Navier-Stokes equations with finite-dimensional controllers[END_REF], [START_REF] Barbu | Tangential boundary stabilization of Navier-Stokes equations[END_REF]. To stabilize the unstable finite-dimensional part using a saturated control we apply techniques, which are well-known for ODE systems (see e.g. [START_REF] Tarbouriech | Stability and Stabilization of Linear Systems with Saturating Actuators[END_REF], [START_REF] Zaccarian | Modern Anti-windup Synthesis: Control Augmentation for Actuator Saturation[END_REF], [START_REF] Saberi | Control of Linear Systems with Regulation and Input Constraints[END_REF] for systems with non-delayed input, and [START_REF] Liu | Delay-dependent methods and the first delay interval[END_REF] for systems with a delayed input).

Using Lyapunov functions, we derive in Section III linear matrix inequalities (LMIs) whose solutions provide estimates of the region of attraction of the closed-loop finite-dimensional system, see [START_REF] Boyd | Linear Matrix Inequalities in System and Control Theory[END_REF], [START_REF] Scherer | Linear matrix inequalities in control[END_REF], [START_REF] Van Antwerp | A tutorial on linear and bilinear matrix inequalities[END_REF] for an introduction to matrix inequalities. We note, that a different Lyapunov function has been used for local ISS stabilization of the diffusion equation by saturated controls in the paper [START_REF] Tanwani | Local input-to-state stabilization of 1-D linear reaction-diffusion equation with bounded feedback[END_REF]. Then we show how asymptotic stability and estimates for the region of attraction can be obtained for the reaction-diffusion equation in closedloop with the nonlinear saturated control.

We show that our results can be extended in a variety of different directions. In Section III-E we demonstrate, how dynamic controllers can be used to enlarge the region of attraction. In Section V-A we show that more complex types of saturation functions can be tackled, such as pointwise (L ∞ ) saturation functions. Although our main results concern the stabilization via distributed controllers, in Section V-B we show how the methods can be applied to unbounded scalar control operators, designed as the output of a finitedimensional boundary control plant, which is fed with the saturated input.

In this light, our approach seems to be useful for any infinite-dimensional systems for which there exist only finitely many unstable modes and which is to be controlled through a bounded input operator (as is the case for systems with input delays considered in e.g. [START_REF] Fiagbedzi | Feedback stabilization of linear autonomous time lag systems[END_REF]), see also Remark 2.

The remaining part of this paper is organized as follows. We first introduce the reaction-diffusion equation with bounded input operator in Section II. The saturation function is defined and the spectral decomposition is provided. A saturated feedback is designed in Section III and an estimate for the region of attraction is provided, first for the open-loop unstable part, and then for both stable and unstable part. Numerical experiments conducted in Section IV illustrate our design method and the obtained estimates of the region of attraction depending on the saturation level. Other saturations are considered in Section V-A. In Section V-B, we show how our results can be applied in the case of unbounded input operators, that is to saturating boundary controllers resulting of a finite-dimensional dynamical system. Concluding remarks and possible future lines of research are collected in Section VI.

Notation: The set of nonnegative reals we denote by IR + . The Euclidean norm on IR n is denoted by | • |, the operator norm induced by this norm on spaces of matrices is denoted by • . The interior of a set S in a topological space is denoted int S and S denotes its closure. In any normed vector space, the ball of radius r around 0 is denoted by B r (0). By IN we denote the nonnegative integers. For convenience,

IN * := IN \ {0}. For k ∈ IN, L > 0, H k (0, L) denotes the
Sobolev space of functions from the space L 2 (0, L), which have weak derivatives of order k, all of which belong to L 2 (0, L). H k 0 (0, L) is the closure of C k 0 (0, L) (the k-times continuously differentiable functions with compact support in (0, L)) in the norm of H k (0, L).

II. PROBLEM FORMULATION

We consider the stabilization problem of a one-dimensional linear reaction-diffusion equation by means of a distributed control u : IR + → IR m . Let L > 0. We are given m functions b k : [0, L] → IR, k = 1, . . . , m which describe at which places the control input u k ∈ IR is acting. The function c models the place-dependent reaction rate. The system model is then

w t (t, x) = w xx (t, x) + c(x)w(t, x) + m k=1 b k (x)sat(u k (t)), t > 0, x ∈ (0, L), w(t, 0) = w(t, L) = 0, t > 0, w(0, x) = w 0 (x), x ∈ (0, L). (1)
We assume that the state space of this system is

X := L 2 (0, L) and that c, b k ∈ X, k = 1, . . . , m.
Here sat is a component-wise saturation function, that is, for all k = 1, . . . , m and for any v

∈ IR m , sat(v) k := v k if |v k | , |v k | v k if |v k | , (2)
where > 0 is the given level of the saturation, which is assumed to be uniform with respect to the index k.

Remark 1: Systems of the form (1) also occur in the problem of stabilizing the linear heat equation by means of boundary control subject to delays or saturation, see e.g. [START_REF] Prieur | Feedback stabilization of a 1-D linear reactiondiffusion equation with delay boundary control[END_REF] as well as Section V-B below.

• For the design of a stabilizing feedback, we use the wellknown eigenfunction decomposition method. Define

A = ∂ xx + c(•)id : X → X (3)
with domain D(A) = H 2 (0, L) ∩ H 1 0 (0, L). Then the control system (1) takes the form

w t (t, •) = Aw(t, •) + m k=1 b k sat(u k (t)). (4)
We note that A is selfadjoint and has compact resolvent, see Appendix B. Hence, the spectrum of A consists of only isolated eigenvalues with finite multiplicity, see [START_REF] Kato | Perturbation Theory for Linear Operators[END_REF]Theorem III.6.29]. Furthermore, there exists a Hilbert basis (e j ) j 1 of X consisting of eigenfunctions of A, associated with the sequence of eigenvalues (λ j ) j 1 . Note that

-∞ < • • • < λ j < • • • < λ 1 and λ j -→ j→+∞ -∞
and that e j (•) ∈ D(A) for every j 1.

We consider (mild) solutions of the system (1) (see [5, Section 3.1]), which exist and are unique for any initial condition in X and for any

u k ∈ L 1,loc ([0, ∞)), for k = 1, . . . , m.
Every solution w(t, •) ∈ D(A) of ( 4) can be expanded as a series in the eigenfunctions e j (•), convergent in H 1 0 (0, L),

w(t, •) = ∞ j=1
w j (t)e j (•),

w j (t) := w(t, •), e j (•) L2(0,L) , j ∈ IN * . (5)
Analogously, we can expand the coefficients b k in the series

b k (•) = ∞ j=1 b jk e j (•), b jk = b k (•), e j (•) L2(0,L) , j ∈ IN * .
As discussed in Appendix A, ( 4) is equivalent to the infinitedimensional control system

ẇj (t) = λ j w j (t) + m k=1 b jk sat(u k (t)) = λ j w j (t) + b j • sat(u(t)), j ∈ IN * , (6)
where "•" is the scalar product in IR m , sat(u(t)) ∈ IR m is the vector with entries sat(u k (t)) and b j is the row vector with entries b jk , k = 1, . . . , m.

Let n ∈ IN * be the number of nonnegative eigenvalues of A and let η > 0 be such that

∀j > n : λ j < -η < 0. (7)
With the matrix notations

z:=    w 1 . . . w n   , A:=    λ 1 • • • 0 . . . . . . . . . 0 • • • λ n   , B:=    b 11 • • • b 1m . . . . . . b n1 • • • b nm    (8)
the n first equations of (6) form the unstable finite-dimensional control system

ż(t) = Az(t) + Bsat(u(t)). (9)
Remark 2: (Possible extensions of our results).

(i) The same method can be used for stabilization with a prescribed decay rate -β. Just define for a suitable η > β the minimal n such that (7) holds, and then use the method described here.

(ii) Since the input operator is bounded and the input space is finite-dimensional, the system can only be stabilized if the unstable spectrum consists of finitely many eigenvalues (counting multiplicities), see [START_REF] Jacob | Equivalent Conditions for Stabilizability of Infinite-Dimensional Systems with Admissible Control Operators[END_REF]Theorem 1.1] and the discussion in that paper relating to the history of the result. If this restriction is satisfied (e.g., this is the case for systems with input delays considered in [START_REF] Fiagbedzi | Feedback stabilization of linear autonomous time lag systems[END_REF]), the methods of this paper apply.

•

III. ESTIMATION OF THE REGION OF ATTRACTION FOR SATURATED INPUTS

A. Decomposition of the system into stable and unstable part

We now introduce a decomposition of the state space into a finite-dimensional space on which the stabilization problem has to be solved and its orthogonal complement, which is invariant under the free dynamics.

Let X n be the subspace of L 2 (0, L) spanned by (e i (•)) n i=1 and π n be the orthogonal projection onto X n , that is

π n w(t, •) := n j=1 w j (t)e j (•). (10)
We define also X ⊥ n as the orthogonal complement of X n in X. Let ι : IR n → X n be the isomorphism defined by ι(e j ) = e j (•), where (e j ) j=1,...,n is the canonical basis of IR n . We will use the isometric representation of L 2 (0, L) as 2 (IN * , IR) obtained by the isomorphism induced by e j (•) → e j , where

2 (IN * , IR) := (x k ) k∈IN * ∈ IR IN * : ∞ k=1 |x k | 2 < ∞ ,
where e j , j ∈ IN * are the standard basis vectors in 2 (IN * , IR) and we use the standard norm on that space. Corresponding to the decomposition L 2 (0, L) = X n X ⊥ n , where " " is the orthogonal sum of subspaces, we denote 2 (IN * , IR) = IR n 2,j>n , where we identify IR n with the sequences with support in {1, . . . , n} and 2,j>n is the set of sequences in 2 (IN * , IR) which are 0 in the first n entries.

Given a linear map K : X n → IR m , consider the feedback

u =Kπ n w(•) = K n j=1 w j e j (•) = n j=1 w j Ke j (•) = n j=1 w j K j = Kz, (11)
where K j := Ke j (•) ∈ IR m , j = 1, . . . , n, and where we use the notation from (8) in the final step and set

K := (K 1 , . . . , K n ) ∈ IR m×n .
Hence the system (6) with the feedback ( 11) is equivalent to the following set of differential equations:

ẇj (t) = λ j w j (t) + b j • sat(Kz(t)), j ∈ IN * . (12)
Using ( 8), we rewrite the first n equations of (12) as

ż(t) = Az(t) + Bsat(Kz(t)). (13)
Now, ( 12) can be considered as a cascade interconnection of an n-dimensional part, described by the equations ( 13) and of an infinite-dimensional part described by the equations

ẇj (t) = λ j w j (t) + b j • sat(Kz(t)), j n + 1. ( 14
)
Our general assumption is: Assumption 1: The pair (A, B) is stabilizable, i.e. there is K ∈ IR m×n : A + BK is Hurwitz.

Remark 3: Clearly, the closed-loop system (13) with a matrix K as in Assumption 1 is locally asymptotically stable.

Conversely, if a feedback K renders the closed-loop system (13) locally asymptotically stable, then also

ż = Az + Bu ( 15
)
is locally and hence globally asymptotically stabilized by means of the feedback u(t) := Kz(t). Thus, local asymptotic stability of ( 13) implies that the pair (A, B) is stabilizable.

We note that in the case m = 1 the situation simplifies further as then ( 15) is a linear diagonal system with scalar control input. The criterion for stabilizability is then that b j1 = 0 for all j = 1, . . . , n and λ k = λ j for all k, j = 1, . . . , n, k = j (which is an easy exercise). That is, the localization function b 1 should not be orthogonal to an unstable eigenfunction and that all unstable eigenvalues need to be simple.

• Next we show that the problem of exponential stabilization of the overall system (6) boils down to the exponential stabilization of the finite-dimensional unstable system [START_REF] Fiagbedzi | Feedback stabilization of linear autonomous time lag systems[END_REF]. This latter problem will be elaborated in Section III-B.

Definition 1: Assume that K is chosen so that 0 is a locally asymptotically stable fixed point of [START_REF] Henry | Geometric Theory of Semilinear Parabolic Equations[END_REF]. We say that S is a region of attraction of 0 if (i) 0 ∈ int S;

(ii) for any z 0 ∈ S the corresponding solution of (13) satisfies z(t; z 0 ) → 0 as t → ∞; (iii) S is forward invariant, i.e. for any z 0 ∈ S it holds that z(t; z 0 ) ∈ S for all t 0. The largest set (with respect to set inclusion) with the properties (i)-(iii) is called the maximal region of attraction.

As unions of regions of attraction are again a region of attraction, it is immediate that the maximal region of attraction is uniquely defined in this way and coincides with what is called domain of attraction in [START_REF] Hahn | Stability of Motion[END_REF].

Definition 2: (13) is called locally exponentially stable in 0 with region of attraction S, if the following conditions hold:

(i) there exist ε, M, a > 0 such that for any initial condition satisfying |z 0 | ε, it holds

|z(t; z 0 )| M e -at |z 0 | ∀t 0 (16)
(ii) B ε (0) ⊂ S and S is a region of attraction of (13).

We call (13) globally exponentially stable in 0 if (16) holds for some a, M > 0 and all z 0 ∈ IR n . Definitions 1 and 2 can be stated analogously for system [START_REF] Hahn | Stability of Motion[END_REF]. We note that if the maximal region of attraction is not IR n , then the system cannot be exponentially stable on the maximal region of attraction (for Lipschitz continuous systems), see [START_REF] Hahn | Stability of Motion[END_REF]. Thus by analyzing regions of attraction with exponential stability we necessarily restrict the region of attraction.

Proposition 1: Assume K is chosen such that the subsystem (13) is locally exponentially stable in 0 with region of attraction S ⊂ IR n . Then: (i) system ( 12) is locally exponentially stable in 0 with region of attraction S × 2,j>n . (ii) system (1) with the feedback (11) is locally exponentially stable in 0 with region of attraction ι(S) × X ⊥ n . In addition, for any closed and bounded set G ⊂ int (ι(S)× X ⊥ n ), there exist two positive values M and a such that for any initial condition w(0, •) in G, the solution w(•) to (1) with the controller [START_REF] Gomes Da | Antiwindup design with guaranteed regions of stability: an LMI-based approach[END_REF] satisfies

w(t, •) X M e -at w(0, •) X ∀t 0. ( 17
)
Proof. Pick a compact subset G of int S. Since we assume that ( 13) is locally exponentially stable in 0 with region of attraction S ⊂ IR n , it follows from a standard compactness argument that there exist M, a > 0 so that for any z 0 ∈ G the solution z(•; z 0 ) to (13) satisfies

|z(t; z 0 )| M e -at |z 0 |, t 0.
From equations ( 6) and [START_REF] Gomes Da | Antiwindup design with guaranteed regions of stability: an LMI-based approach[END_REF], we derive that for j = n + 1, . . . , ∞, for any t 0 and for any (w n+1 (0), w n+2 (0), . . .) ∈ 2,j>n it holds that w j (t) = e λj t w j (0) + b j • t 0 e λj (t-s) sat(Kz(s))ds.

From (2) it follows that for all z ∈ IR n we have

|sat(Kz)| |Kz| K |z|.
Also due to the Cauchy-Bunyakovsky-Schwarz inequality we have that, for all j ∈ IN,

|b jk | = b k (•), e j (•) X b k X e j X = b k X . (18)
Thus, for all j n + 1, we obtain exploiting (18) that We now obtain local exponential stability of ( 12) by choosing G such that 0 ∈ IR n is in the interior of G and noting that then 0 ∈ X is in the interior of ι(G ) × X ⊥ n . For the final statement of the proposition, pick a closed and bounded set

|w j (t)| e -ηt |w j (0)| + |b j | t 0 e -η(t-s) |Kz(s)|ds e -ηt |w j (0)| + |b j | K t 0 e -η(t-s) M e -as |z(0)|ds =e -ηt |w j (0)| + |b j | M K η -a (e -at -e -ηt )|z(0)|.
G ⊂ int (ι(S) × X ⊥ n ). Select G = ι -1 • π n (G)
, then G is a compact subset of int S, the previous computations yield [START_REF] Kang | Boundary constrained control of delayed nonlinear Schrödinger equation[END_REF] for suitable constants M and a and for the superset ι(G ) × X ⊥ n which contains G. Remark 4: It is not hard to see that ( 14) is input-tostate stable (ISS) with respect to the input z. Hence ( 12) is asymptotically stable as a cascade interconnection of a locally asymptotically stable system and an ISS system. However, since the general theorem on cascade interconnections does not guarantee exponential convergence, we needed an extra argument for Proposition 1.

• Remark 5: Our feedback controller is robust w.r.t. additive actuator disturbances. Let the control input to system (1) be u(t) := Kz(t) + d(t), where d ∈ P C(IR + , IR m ) is a piecewise continuous actuator disturbance. Then there are r > 0, M 1 , a 1 , γ 1 > 0 so that for all w(0, •) ∈ B r and all d ∈ P C(IR + , IR m ) with sup s 0 |d(s)| < r the solution of

w t (t, x) = w xx (t, x) + c(x)w(t, x) + m k=1 b k (x)sat (Kz(t)) k + d k (t) , t > 0, x ∈ (0, L), w(t, 0) = w(t, L) = 0, t > 0, w(0, x) = w 0 (x), x ∈ (0, L). (19)
satisfies the local input-to-state stability (LISS) estimate 20) follows as exponential stability of ( 19) without disturbances implies LISS for a system [START_REF] Krstic | Control of an unstable reaction-diffusion PDE with long input delay[END_REF] with disturbances. It is, however, much harder to estimate a region of attraction of system [START_REF] Krstic | Control of an unstable reaction-diffusion PDE with long input delay[END_REF], as in addition to the complexities arising in the undisturbed case the interplay between the size of a region of attraction and the maximal norm of a disturbance has to be analyzed. This can be an interesting topic for a future research. For more on ISS theory of infinite-dimensional systems the reader may consult [START_REF] Mironchenko | Input-to-state stability of infinitedimensional systems: recent results and open questions[END_REF], [START_REF] Mironchenko | Characterizations of input-to-state stability for infinite-dimensional systems[END_REF], [START_REF] Tanwani | Disturbance-to-state stabilization and quantized control for linear hyperbolic systems[END_REF] and references therein. •

w(t, •) X M 2 e -a2t w(0, •) X + γ 2 sup s 0 |d(s)|. (20)

B. Estimate of the region of attraction for the finitedimensional part

In view of Proposition 1, it is important to study the local exponential stability and to estimate the region of attraction of the finite-dimensional system [START_REF] Henry | Geometric Theory of Semilinear Parabolic Equations[END_REF]. We perform this task in this section.

Defining the deadzone nonlinearity φ : IR m → IR m by

φ(u) = sat(u) -u, u ∈ IR m ,
where sat is defined in (2), the following generalized sector condition holds (see [41, Lemma 1.6, Page 45] for a proof):

Lemma 1: If for some C ∈ R m×n , z ∈ R n and j ∈ {1, . . . , m} it holds that |((K -C)z) j | , then φ j (Kz)(φ j (Kz) + (Cz) j ) 0.
Remark 6: Note that other properties exist for saturation yielding to other representations of the saturation maps (see, e.g., [START_REF] Tarbouriech | Stability and Stabilization of Linear Systems with Saturating Actuators[END_REF]Section 1.7]). In the present work, we insist on using generalized sector conditions because they provide constructive methods to design Lyapunov functions.

• As a consequence of Lemma 1, for all diagonal positive definite matrices D ∈ R m×m , for all C ∈ R m×n , and all z ∈ R n such that |((K -C)z) j | , j = 1 . . . , m, we have:

φ(Kz) D(φ(Kz) + Cz) 0. (21)
We recall the following well-known Schur complement lemma (see e. The following result is only a small variation of [START_REF] Gomes Da | Antiwindup design with guaranteed regions of stability: an LMI-based approach[END_REF]Theorem 1]. Let us give a full proof using the notation used in this paper, for the sake of completeness.

Proposition 2: Under Assumption 1, let K ∈ IR m×n be such that A + BK is Hurwitz. Let P ∈ R n×n be symmetric positive definite, D ∈ R m×m diagonal positive definite and C ∈ R m×n such that

M 2 := P (K -C) K -C 2 I m 0. ( 23
)
Then the finite-dimensional system (13) is locally exponentially stable in 0 with a region of attraction given by

A := {z : z P z 1}. (24)
Moreover, in A, the function V 1 defined by V 1 (z) := z P z, z ∈ IR n , decreases exponentially fast to 0 along the solutions to (13), i.e. there is a constant α > 0 so that

V1 (z) -α|z| 2 , z ∈ A. (25)
Proof. Let us show that ( 22), ( 23) are feasible. As A + BK is Hurwitz, there is a symmetric positive definite matrix P ∈ R n×n such that (A + BK) P + P (A + BK) < 0 and, multiplying if needed P by a positive constant, we may additionally assume that P K K 2 I m 0. Then, letting C = 0 ∈ R m×n , and picking a diagonal positive definite matrix D ∈ R m×m with a norm sufficiently large, we get diagonal positive definite matrix D ∈ R m×m and a matrix C ∈ R m×n such that ( 22) and ( 23) hold. Now pick P, D, C such that 22) and ( 23) hold and consider the Lyapunov function candidate

V 1 (z) = z P z,
where P ∈ IR n×n is the symmetric positive definite matrix given by the assumptions. The time-derivative of V 1 along the solutions to (13) for z ∈ IR n is given by

V1 (z) = z [(A + BK) P + P (A + BK)]z + 2z P Bφ(Kz). Now assuming |((K -C)z) j |
for all j = 1, . . . , m, it follows from (21) that

V1 (z) z [(A + BK) P + P (A + BK)]z +2z P Bφ(Kz) -2φ(Kz) D(φ(Kz) + Cz) = z φ(Kz) M 1 z φ(Kz)
.

In view of [START_REF] Liu | Delay-dependent methods and the first delay interval[END_REF], this implies the estimate (25) selecting -α as the maximal eigenvalue of M 1 .

It remains to ensure that |((K -C)z) j | is satisfied for all j = 1, . . . , m. To do this, we use the Lyapunov function V 1 and we impose that the ellipsoid {z ∈

IR n : z P z 1} is included in the ellipsoid {z ∈ IR n : |(K -C)z| } (this implies that |((K -C)z) j |
holds for all j = 1, . . . , m). This inclusion is equivalent to the inclusion of {z ∈ IR n : z P z 1} in the set {z ∈ IR n :

z (K -C) (K -C)z
2 } which is again equivalent to the matrix inequality P -(K -C) 1 2 (K -C) 0. As 2 I m is positive definite, Lemma 2 ensures that this latter matrix inequality is equivalent to [START_REF] Marx | Cone-bounded feedback laws for mdissipative operators on Hilbert spaces[END_REF].

Remark 7: We note that the formulation of Proposition 2 immediately gives room for a larger estimate for the domain of attraction so that the estimate given by {z : z P z 1} can never be optimal. Indeed, in this region we have V1 (z) -α|z| 2 , so that by a continuity and compactness argument it follows that V1 (z) < 0 on an enlarged region of the form {z : z P z 1 + ε}, for a suitable ε > 0.

• With the previous result, it is also possible to analyze global stability. The region of attraction is global as soon as for all z ∈ R n , it holds that |((K-C)z) j | , j = 1, . . . , m. This is equivalent to K = C. We thus obtain the following corollary:

Corollary 1: If there exist a symmetric positive definite matrix P ∈ R n×n and a diagonal positive matrix D ∈ R m×m such that [START_REF] Liu | Delay-dependent methods and the first delay interval[END_REF] holds with C := K, then the finite-dimensional system (13) is globally exponentially stable in 0. Moreover the Lyapunov function V 1 decreases exponentially fast to 0 along the solutions to [START_REF] Henry | Geometric Theory of Semilinear Parabolic Equations[END_REF]. In our context, the main interest of Proposition 2 lies in the following consequence for system [START_REF] Barbu | Tangential boundary stabilization of Navier-Stokes equations[END_REF].

Theorem 1: Consider system (1) with (A, B) stabilizable. Let K be such that A + BK is Hurwitz and K be as in [START_REF] Gomes Da | Antiwindup design with guaranteed regions of stability: an LMI-based approach[END_REF]. Then the closed-loop system with the feedback K w t (t, x) = w xx (t, x) + c(x)w(t, x)

+ m j=1 b k (x)sat (Kπ n w(t, •)) k , t > 0, x ∈ (0, L), w(t, 0) = w(t, L) = 0, t > 0, (26)
w(0, x) = w 0 (x), x ∈ (0, L). is locally exponentially stable in 0 with region of attraction ı(A) × X ⊥ n . In addition, the constants of decay can be chosen uniformly on ı(A) × X ⊥ n . Proof. Exponential stabilization and the region of attraction are immediate consequences of Proposition 2 in combination with Proposition 1. As the exponential estimate can be chosen uniformly on the set A, the proof of Proposition 1 shows that the uniformity holds on all of ı(A) × X ⊥ n .

C. Numerical implementation of Proposition 2

The inequalities ( 22) and ( 23) are bilinear matrix inequalities, since they have the bilinear cross term DC in the unknowns D and C. The nonlinearity complicates the numerical analysis of ( 22) and ( 23). However, they can be reformulated into equivalent linear matrix inequalities (LMIs), which are easier to solve numerically.

Let S = P -1 and E = D -1 . Since M 1 < 0 if and only if

S 0 0 E M 1 S 0 0 E < 0,
we obtain that ( 22) is equivalent to

S(A + BK) + (A + BK)S BE -SC (BE -SC ) -2E < 0. ( 27
)
Now pre-and post-multiplying ( 23) by S 0 0 I m , we get that ( 23) is equivalent to

S SK -SC (SK -SC ) 2 I m 0. (28)
Letting Y = SC , we obtain (see also [11, Theorem 1] in a different context): Proposition 3: Existence of a symmetric positive definite matrix P ∈ R n×n , a diagonal positive definite matrix D ∈ R m×m and a matrix C ∈ R m×n such that ( 22) and ( 23) hold is equivalent to the existence of a symmetric positive definite matrix S ∈ R n×n , a diagonal positive definite matrix E ∈ R m×m and Y ∈ R m×n for which it holds that

S(A + BK) + (A + BK)S BE -Y (BE -Y ) -2E < 0 (29)
and

S SK -Y (SK -Y ) 2 I m 0. ( 30
)
This reformulation is important as the constraints ( 29) and ( 30) are linear in the variables S, E and Y, in contrast to bilinear inequalities ( 22) and ( 23). Finally, no structure is imposed on both the matrix C ∈ R m×n and on the matrix Y ∈ R m×n , and we let C = (S -1 Y) from the knowledge of S and Y (and respectively we may let Y = P -1 C from the knowledge of P and C).

D. Scalar control inputs

In this section we specialize our results for scalar control inputs, i.e. m = 1, which allows for a numerically more efficient method. We need the following lemma.

Lemma Weyl's inequality (see [38, Chapter IV, Corollary 4.9, p. 203]) implies that all eigenvalues of P a /P 22 are strictly increasing functions of a with a slope bounded below everywhere by λ min (P 11 ) > 0. The claim is now immediate. Clearly, P 0 is not positive definite. As all eigenvalues are strictly increasing with an affine lower bound, eventually at some a * 0 the smallest eigenvalue is equal to 0 by continuity of eigenvalues. For all a > a * all eigenvalues are positive.

The main result of this section is Proposition 4: Under Assumption 1, let K ∈ IR m×n be such that A + BK is Hurwitz. Let P ∈ R n×n be symmetric positive definite and C ∈ R 1×n such that

M1 := (A + BK) P + P (A + BK) P B -C ( P B) -C -2 < 0. ( 31
)
Then the finite-dimensional system (13) is locally exponentially stable in 0 with a region of attraction given by

{z : z P z D -1 }, ( 32
)
where D is the minimal real number such that:

M2 := D P (K -C) K -C 2 0. (33)
Moreover, in this region of attraction, the function V 1 defined by V 1 (z) = z P z := Dz P z, for all z ∈ IR n , decreases exponentially fast to 0 along the solutions to (13), i.e. there is a constant α > 0 so that

V1 (z) -α|z| 2 . ( 34
)
Proof. Define P := D P . Substituting this expression into [START_REF] Liu | Delay-dependent methods and the first delay interval[END_REF] and [START_REF] Marx | Cone-bounded feedback laws for mdissipative operators on Hilbert spaces[END_REF] and multiplying it by D -1 to obtain the inequalities [START_REF] Prieur | Feedback stabilization of a 1-D linear reactiondiffusion equation with delay boundary control[END_REF] and [START_REF] Saberi | Control of Linear Systems with Regulation and Input Constraints[END_REF]. If there exist P and C as in Proposition 4 so that the LMI (31) holds, then according to Lemma 3, one can find D > 0 so that (33) holds as well.

The estimate for the region of attraction follows from Proposition 2. We choose the minimal D to obtain the largest estimate of a region of attraction (for given P and C).

Remark 8: Lemma 3 can be helpful in finding the optimal D as it shows that the problem is to find the unique root of a piecewise analytic function.

•

E. Enlarging the region of attraction using a dynamic controller

Consider a dynamic controller for system [START_REF] Henry | Geometric Theory of Semilinear Parabolic Equations[END_REF], instead of a static controller. By doing so, we add some degrees of freedom and thus we may enlarge the region of attraction. Such approach is useful for many control systems (see e.g. [START_REF] Tarbouriech | Stability and Stabilization of Linear Systems with Saturating Actuators[END_REF]Example 8.1] for a simple 2 × 2 example). To be more specific, we consider an additional finite-dimensional state z c in IR nc (for a given integer n c ) and design matrices K 1 , K 2 , A c and B c of appropriate dimensions so that the (estimation of the) region of attraction of

ż = Az + Bsat(K 1 z + K 2 z c ) żc = A 1 z c + A 2 z ( 35
)
is larger (in the z-direction) than the estimate provided by Proposition 2 for [START_REF] Henry | Geometric Theory of Semilinear Parabolic Equations[END_REF]. We rewrite the dynamics [START_REF] Scherer | Linear matrix inequalities in control[END_REF] by

Ż = ĀZ + Bsat( KZ) (36)
where

Ā = A 0 A 2 A 1 , B = B 0 , K = K 1 K 2 , ( 37
)
with Ā ∈ IR (n+nc)×(n+nc) , B ∈ IR (n+nc)×m and K ∈ IR m×(n+nc) .

Remark 9: For given symmetric positive definite matrices P ∈ IR n×n and P ∈ IR (n+nc)×(n+nc) , the inclusion of the ellipsoid {z ∈ IR n : z P z 1} in the projection (onto IR n ) of the ellipsoid {Z ∈ IR n+nc : Z P Z 1} is equivalent to

I n 0 P I n 0 -P 0 . (38)
• Using this remark and applying Proposition 2 to system [START_REF] Scherer | Linear matrix inequalities in control[END_REF], we have the following:

Proposition 5: Consider system (13) with a stabilizable pair (A, B) ∈ IR n×n × IR n×m , and matrices Ā, B and K defined in [START_REF] Solomon | Stability and passivity analysis of semilinear diffusion PDEs with time-delays[END_REF]. Choose a symmetric positive definite matrix P ∈ R (n+nc)×(n+nc) , a diagonal positive matrix D ∈ R m×m and C ∈ R m×(n+nc) such that (38) holds and

( Ā + B K) P + P ( Ā + B K) P B -( D C) ( P B) -D C -2 D < 0, (39)
P ( K -C) K -C 2 I m 0. ( 40
)
Then the finite-dimensional system ( 35) is locally exponentially stable in 0 with a region of attraction given by A P := {Z ∈ IR n+nc : Z P Z 1}. Moreover, the projection of the ellipsoid {Z : Z P Z 1} onto IR n is larger than the region of attraction A given by Proposition 2 in (24) for system [START_REF] Henry | Geometric Theory of Semilinear Parabolic Equations[END_REF].

Finally, in A P , the function V 1 defined by

V 1 (Z) = Z P Z, Z ∈ IR n+nc ,
decreases exponentially fast to 0 along the solutions to [START_REF] Scherer | Linear matrix inequalities in control[END_REF], i.e. there is a constant α > 0 so that for all

Z ∈ A P V1 (Z) -α|Z| 2 . ( 41
)
As in Section III-C, the inequalities ( 39) and ( 40) can be reformulated as linear matrix inequalities.

F. Lyapunov analysis of the closed-loop system under saturation control

Now under the assumptions of Proposition 2, we provide a Lyapunov function V : 2 (IN * , IR) → IR + for the system [START_REF] Hahn | Stability of Motion[END_REF]. To this end, we use the Lyapunov function P provided by Proposition 2 for the finite-dimensional subsystem [START_REF] Henry | Geometric Theory of Semilinear Parabolic Equations[END_REF].

Proposition 6: Consider system (12) and assume K is chosen such that the subsystem (13) is locally exponentially stable in 0. Assume further that P ∈ IR n×n is symmetric positive definite and such that for Ṽ (z) = z P z, z ∈ IR n we have V (z) -α|z| 2 on the set A := {z ∈ IR n : Ṽ (z) 1} along the solutions of [START_REF] Henry | Geometric Theory of Semilinear Parabolic Equations[END_REF]. Then there exist γ, C > 0 such that for Q :

2 (IN * , IR) → IR defined by Q(w) := ∞ j=n+1 w 2 j the function V (w) := Ṽ (π n w) + γQ(w) = z P z + γQ(w) (42)
(with the identification π n w = z) satisfies for all w ∈ A × 2,j>n that along the solutions of ( 12)

V (w) -C w 2 2 .
In particular, it follows that S × 2,j>n is a region of attraction of 0 for system [START_REF] Hahn | Stability of Motion[END_REF].

Proof. We write w ∈ 2 (IN * , IR) as w = z + z ⊥ , where z ∈ IR n , z ⊥ ∈ 2,j>n . Here we identify IR n with the sequences with support in {1, . . . , n} and 2,j>n is the set of sequences in 2 (IN * , IR) which are 0 in the first n entries.

This decomposition is unique. Due to [START_REF] El-Farra | Analysis and control of parabolic PDE systems with input constraints[END_REF], we have along the solutions to ( 12)

Q(w) =2 ∞ j=n+1 w j (t) ẇj (t) =2 ∞ j=n+1 w j (t) λ j w j (t) + b j • sat(Kz(t)) -2η ∞ j=n+1 w 2 j (t) + 2 ∞ j=n+1 |w j (t)||b j ||sat(Kz(t))| -2η ∞ j=n+1 w 2 j (t) + 2 ∞ j=n+1 |w j (t)||b j | K |z(t)|.
Cauchy-Bunyakovsky-Schwarz inequality implies that

Q(w) -2η z ⊥ (t) 2 2 + 2 z ⊥ (t) 2 b ⊥ (t) 2 K |z(t)|. Using Young's inequality we have for all κ > 0 that 2 z ⊥ (t) 2 |z(t)| κ z ⊥ (t) 2 2 + 1 κ |z(t)| 2 .
We proceed to:

Q(w) -2η -κ b ⊥ 2 K z ⊥ (t) 2 2 + 1 κ b 2 K |z(t)| 2 .
Therefore for any choice of γ in (42), we have by an application of Proposition 2 along the solutions to (12) that

V (w) -α - γ κ b ⊥ 2 K |z(t)| 2 -γ 2η -κ b ⊥ 2 K z ⊥ (t) 2 2 . Therefore selecting first κ > 0 such that 2η-κ b ⊥ 2 K > 0 and then selecting γ > 0 such that α -γ κ b ⊥ 2 K > 0, we get the existence of C > 0 such that V (w) -C|z(t)| 2 -C z ⊥ (t) 2 2 (43)
provided that z(t) lies in the ellipsoid A, whatever the value of z ⊥ (t). As A × 2,j>n is an invariant subset of 2 for a system [START_REF] Hahn | Stability of Motion[END_REF], we obtain that V is a Lyapunov function for [START_REF] Hahn | Stability of Motion[END_REF], with a guaranteed region of attraction containing A × 2,j>n .

Note that this provides an alternative proof of Proposition 1. We can also use Proposition 5 and obtain a similar result for a dynamical controller.

IV. NUMERICAL EXPERIMENTS

In this section we use Proposition 4 to obtain estimates for a region of attraction for the unstable heat equation (1) subject to a saturated scalar feedback controller. Let c(•) in equation (1) be a constant function. With slight abuse of notation we will write c(•) = c = const.

According to [13, pp. 16-17] the eigenvalues of the operator

A := ∂ xx + c id : X → X (44)
on the domain D(A) = H 2 (0, L) ∩ H 1 0 (0, L) are given by

λ j := - π 2 L 2 j 2 + c, j ∈ IN * , (45)
and the eigenfunctions e j , j ∈ IN * of (A, D(A)), which form a basis of L 2 (0, 1) are given by

e j (x) := 2 L 1/2 sin jπx L , j ∈ IN * , x ∈ (0, L). ( 46
)
Next we estimate the region of attraction of system (1) for the following choice of parameters:

c(x) ≡ 10, L = 2, = 2, b = e 1 + e 2 ,
where the e j are defined in (46). For these parameters there are only 2 unstable eigenvalues. Using the eigenfunction decomposition of solutions of (1) as in Section II, we see that the matrices A, B defined in (8) have the form

A ≈ 7.5325989 0 0 0.1303956 , B = 1 1 .
As the diagonal entries in the matrix A are distinct, and all the components of B are nonzero, the pair (A, B) is stabilizable. Different choices of the matrix K for the stabilizing feedback u(t) = Kz(t) lead to different attraction rates and different regions of attraction. We demonstrate this with two examples.

A. Choice 1: Placing the poles at (-1, -1)

First we choose the matrix K so that σ(A + BK) = {-1}, which results in K ≈ -9.835618 0.1726235 .

To estimate a region of attraction of ( 13) by Proposition 4, we proceed in two steps.

(i) First we solve the inequality (31) together with the additional constraints: P -P = 0 and P > 0. Additionally, we impose an optimality condition for C:

(K -C) • (K -C) → min, ( 48
)
where • is the scalar product of vectors. (ii) The idea behind (48) is to minimize the off-diagonal elements of the matrix M2 , which gives us at the second step the possibility to find large D (optimal for the given P , C) satisfying the bilinear matrix inequality [START_REF] Saberi | Control of Linear Systems with Regulation and Input Constraints[END_REF].

This algorithm is implemented in Scilab. For solution of the LMI (31) the LMITOOL package has been used. The resulting matrices P , C and the real number D are (approximately):

P = 2.1277468 -0.0655569 -0.0655569 0.0243008 , (49a)
C =(-2.0635579, 0.0844904), D = 7.359375. (49b)

Figure 1 shows an elliptic region of attraction [START_REF] Marx | Global stabilization of a Korteweg-De Vries equation with saturating distributed control[END_REF], subject to P , C, D given by (49) (in blue). Furthermore, in the same figure some trajectories are depicted obtained through direct simulation of [START_REF] Henry | Geometric Theory of Semilinear Parabolic Equations[END_REF]. The trajectories in black tend to the origin while those in red are diverging. This provides an approximation of the maximal region of attraction of [START_REF] Henry | Geometric Theory of Semilinear Parabolic Equations[END_REF]. It can be seen that in one direction the ellipsoid obtained by our method approximates the actual region of attraction very well, but the results are not tight in the orthogonal direction.

Using Theorem 1, we obtain an approximation of the region of attraction of the PDE (1) subject to the feedback [START_REF] Gomes Da | Antiwindup design with guaranteed regions of stability: an LMI-based approach[END_REF].

Remark 10: (Importance of the optimality conditions) Different solutions P , C, D of the matrix inequalities ( 31), [START_REF] Saberi | Control of Linear Systems with Regulation and Input Constraints[END_REF] lead to very different estimates of a region of attraction [START_REF] Marx | Global stabilization of a Korteweg-De Vries equation with saturating distributed control[END_REF] of the model [START_REF] Henry | Geometric Theory of Semilinear Parabolic Equations[END_REF]. Thus, it is important to pick solutions resulting in as large regions of attraction as possible. Here we chose a solution satisfying the optimality condition (48). Enforcing further optimality conditions may provide different estimates. The union of regions of attractions is again a region of attraction and so the maximal region can be explored further by solving [START_REF] Prieur | Feedback stabilization of a 1-D linear reactiondiffusion equation with delay boundary control[END_REF], [START_REF] Saberi | Control of Linear Systems with Regulation and Input Constraints[END_REF] for different optimality conditions.

In the two dimensional case, one option is to fix the eigendirections of P and to optimize the eigenvalues to get an idea of the extension of the maximal domain of attraction in particular directions. Other size criteria exist for optimization of domains of attraction, as volume maximization and trace minimization. See [41, Section 2.2.5] for more details and a complete introduction of such optimization problems.

• Remark 11: (Computational cost) The time needed to solve the problem was (on a system with the specs: Intel(R) Core(TM) i5-3317U 1.70GHz, 16 GB RAM, Windows 10):

• Finding P , C, D via LMIs: 0.0166561 seconds.

• Plotting the obtained region: 0.0071209 seconds.

• Time for solving the ODE (13) for 31 2 = 961 distinct initial conditions on the time-interval [0, 10] on a grid consisting of 100 points and for plotting the resulting trajectories: 43.359664 seconds.

This shows the computational efficiency of our method. • B. Choice 2: Putting the poles to {-0.1, -0.2}

Now let us choose the matrix K so that σ(A + BK) = {-0.1, -0.2}. This choice makes the attraction rate of the closed-loop system much slower, than in the previous simulation. This has however, some advantages, as we will see next. The resulting matrix K is:

K = -7.9732782 0.0102837 . ( 50
)
As in the previous simulation, we solve the LMI (31) subject to the additional optimality condition (48). The corresponding matrices P , C, D are:

P = 0.3108695 -0.0054849 -0.0054849 0.000195 , (51a)
C =(-0.3053879, 0.0054754), D = 90.625.

As we see, with the choice of the stabilizing feedback (50), the region of attraction becomes significantly larger, although at the cost of reducing the rate of convergence of the trajectories to the origin. Furthermore, the choice of the matrices (51) leads to a better estimate of the region of attraction, in comparison to the situation in Section IV-A. This can be seen by comparing the Figures 1 and2.

Remark 12: (Computational costs) For this problem the elapsed time is (on a system with the specs: Intel(R) Core(TM) i5-3317U 1.70GHz, 16 GB RAM, Windows 10)

• Finding P , C, D via LMIs: 0.018131 seconds.

• Plotting the obtained region: 0.0129341 seconds.

• Time for solving the ODE (13) for 31 2 = 961 distinct initial conditions on the time-interval [0, 60] on a grid consisting of 600 points and for the plotting of the resulting trajectories: 91.802833 seconds. We have chosen a longer time-span for solution of the ODE [START_REF] Henry | Geometric Theory of Semilinear Parabolic Equations[END_REF], since in this simulation the attraction rate of the closed loop system is much slower than in the previous simulation. Again, we obtain a considerable approximation of the region of attraction in a computationally efficient way.

• V. EXTENSIONS

A. Pointwise saturations

The type of the saturation which we have considered until now, i.e. component-wise saturation of finite-dimensional vectors, is not the only type of saturation functions, which appears in engineering practice. A general class of saturation functions has been considered in [START_REF] Marx | Global stabilization of a Korteweg-De Vries equation with saturating distributed control[END_REF].

1) Saturation functions: Various definitions for the saturation map exist. Normwise saturations limit the norm of the input u, i.e. for a given norm • on IR n we may consider sat

• (u) := u , if u u u , if u = min 1 , 1 u u. ( 52
)
Another physically motivated saturation map is the following (pointwise) L ∞ saturation map, defined for all x ∈ (0, L) by

sat ∞ (u)(x) := u(x) , if |u(x)| u(x) |u(x)| , if |u(x)| . = min 1 , 1 |u(x)| u(x). (53)
In this section, we depart from the saturation model in [START_REF] Barbu | Tangential boundary stabilization of Navier-Stokes equations[END_REF] and study instead a heat equation with pointwise saturation in each input channel:

w t (t, x)=w xx (t, x)+c(x)w(t, x)+ m k=1 sat(b k (x)u k (t)). ( 54
)
In a certain sense, in the equation (54) the whole terms u k b k are considered as the input which saturates pointwise. The assumptions on the domain of the problem and the functions c, b k are the same as in Section II.

In order to stabilize (54), we are going to use the same stabilizing control [START_REF] Gomes Da | Antiwindup design with guaranteed regions of stability: an LMI-based approach[END_REF]. We now aim to provide an estimate for a region of attraction for (54). So we assume that K as in [START_REF] Gomes Da | Antiwindup design with guaranteed regions of stability: an LMI-based approach[END_REF] is given and we consider the closed-loop system

w t (t, x) = w xx (t, x)+c(x)w(t, x)+ m k=1 sat b k (x) (Kz(t)) k .
Representing this equation in the basis e j , j ∈ IN * as in Section II, we obtain the equations for the coordinates

ẇj (t) =λ j w j (t) + m k=1 e j , sat ∞ b k (•) (Kz(t)) k , j ∈ IN * . ( 55
)
Let us state the following simple result that will be instrumental to bound the differences

∆ (Kz(t)) k , b k (•) := b k (•)sat Kz(t) k -sat ∞ b k (•)Kz(t) k . Lemma 4: Let r, k ∈ IR and denote ∆(r, k) = rsat(k) - sat(rk). Then |k| and |rk| ⇒ ∆(r, k) = 0. ( 56
) |∆(r, k)| (1 + |r|). (57)
Let k ∈ IR and b ∈ L 2 (0, L) be given. Then

|k| and b(•)k ∞ ⇒ ∆(b(•), k) ≡ 0 a.e. ( 58
) Moreover b(•)sat(k) -sat ∞ (b(•)k) 2 ( 1 + |b| 2 ). (59)
If |k| and χ is the characteristic function of the set The remaining claims follow immediately by applying the pointwise estimates (56) and (57). The claim (60) follows by applying (59) to bχ, after noting that the complement of U does not contribute to the norm of the left hand side.

U := {x ∈ [0, L] : |kb(x)| > } then b(•)sat(k) -sat ∞ (b(•)k) 2 ( χ + |bχ| 2 ). ( 60
) If b(•) is essentially bounded, then b(•)sat(k) -sat ∞ (b(•)k) ∞ (1 + b ∞ ), (61
Proposition 7: Consider system (54) and assume all functions b k (•) ∈ L ∞ ([0, L]), k = 1, . . . , m. Consider also the associated system (13) with a stabilizable pair (A, B) ∈ IR n×n × IR n×m and let K ∈ IR m×n be such that A + BK is Hurwitz. Choose a symmetric positive definite matrix P ∈ R n×n as in Proposition 6. Consider the Lyapunov function candidate V defined in [START_REF] Teel | Global stabilization and restricted tracking for multiple integrators with bounded controls[END_REF].

Then there exist α, β > 0 such that if for w

= z + z ⊥ ∈ 2 (IN * , IR) we have z ∈ A β := {z : z P z β} then V (w) - α 2 |w| 2 . ( 62
)
In particular, {z : z P z β} × X ⊥ n is a region of attraction for system (55).

Proof. First we rewrite (54) by adding and subtracting the term b k (x)sat(u k (t)) to obtain

w t (t, x) = w xx (t, x) + c(x)w(t, x) + m k=1 b k (x)sat(u k (t)) + m k=1 sat b k (x)u k (t) -b k (x)sat u k (t) . (63) Define ∆(t) := m k=1 sat ∞ b k (•) (Kz(t)) k -b k (•)sat(Kz(t) k ). ( 64
)
For j ∈ IN * we define y j (t) := e j , ∆(t) , and let y(t) ∈ IR n be the vector with components y j (t), for j = 1, . . . , n.

Considering the Lyapunov function V defined in (42), we compute its time-derivative along the solutions to (55). Using [START_REF] Van Antwerp | A tutorial on linear and bilinear matrix inequalities[END_REF], we get for all w = z + z ⊥ with z P z 1 that

V (w) -C|z(t)| 2 -C z ⊥ (t) 2 2 + 2z(t) P y(t) + 2γ ∞ j=n+1 w j (t)y j (t). (65)
Using (59) from Lemma 4, we obtain along the solutions to (54) and as long as z P z 1,

V (w) -C|z(t)| 2 -C z ⊥ (t) 2 2 + 2λ max |z(t)| ∆(t) ∞ + 2γ z ⊥ (t) 2 ∆(t) ∞ ,
where λ max denotes the maximal eigenvalue of the matrix P . Therefore for any positive values κ and κ ,

V (w) -C - λ max κ |z(t)| 2 -(C - γ κ ) z ⊥ (t) 2 2 + (λ max κ + γκ ) ∆(t) 2 ∞ . Pick κ > 0 and κ > 0 such that C -λmax κ > 3C 4 and C -γ κ > C 2 .
Due to (58), there exists β > 0, such that for all z in {z :

z P z β}, ∆(t) 2 ∞ C 4(λmaxκ+γκ ) |z(t)| 2 .
We get, for all solutions to (54), as long as z(t) is in {z : z P z β},

V (w) - C 2 |z(t)| 2 - C 2 z ⊥ (t) 2 2 (66)
Moreover, we have, along the solutions to (54), ẇj (t) = λ j w j (t) + sat(bKz(t)), e j , j = 1, 2, . . . (67) therefore, considering V 1 as previously defined, we may check that (62) holds following the same computation as for V along the solutions to (54), and using the fact the dynamics (67) in X n does not depend on the component in X ⊥ n . Therefore the set {z : z P z β} × X ⊥ n is invariant along the dynamics to (54). With (66), we get that {z : z P z β} × X ⊥ n is a region of attraction and V is a Lyapunov function.

B. Applications to boundary control of heat equation subject to control saturations

Let us now start from a heat equation with a dynamical boundary condition y t (t, x) = y xx (t, x) + c(x)y(t, x), t 0, x ∈ (0, L), (68a) y(t, 0) = 0, y(t, L) = y d , t 0, (68b) where y d is the (scalar) output of finite-dimensional dynamical system given by ẋd =A d x d + B d sat(u(t)), (69a)

y d =C d x d . (69b)
Here x d in IR n d is the finite-dimensional state and the dynamics are subject to a saturating control, A d , B d and C d are three matrices of appropriate dimension, and u is the scalar control input for the PDE (68) and the ODE (69) that is subject to a saturation map. Inspired by [START_REF] Prieur | Feedback stabilization of a 1-D linear reactiondiffusion equation with delay boundary control[END_REF], we introduce the following change of variable: w(t, x) = y(t, x) -x L y d (t), t 0, x ∈ (0, L).

The PDE for w then reads as: The boundary conditions for the variable w take the form:

w t (
w(t, 0) = w(t, L) = 0, t 0. ( 71
)
The heat equation ( 70), (71) has to be analyzed along with the ODE (69a).

Performing similar computations as in Section II for the PDE (1) and, using the same notation for w j and λ j , we get ẇj (t) = λ j w j (t) + b j sat(u(t)) + As the control input is scalar, we can apply Proposition 4 to system (72) instead of system [START_REF] Henry | Geometric Theory of Semilinear Parabolic Equations[END_REF], and get sufficient conditions for the estimation of the region of attraction of (72). Coming back to the infinite-dimensional systems (70) and (68), and, applying Proposition 1, we get sufficient conditions for an estimation of attraction region of (68):

Corollary 2: Consider system (13) with a stabilizable pair (A, B) ∈ IR n×n × IR n×1 and let K ∈ IR 1×n be such that A + BK is Hurwitz. Pick a symmetric positive definite P ∈ R (n+n d )×(n+n d ) and a C ∈ R 1×(n+n d ) such that (31) holds.

Then the finite-dimensional system (72) is locally exponentially stable in 0 with a region of attraction given by

A := {z : z P z D -1 }, ( 73
)
where D is the minimal real number such that (33) holds. Moreover, (i) (70) is locally exponentially stable with a region of attraction ı(A) × X ⊥ n , (ii) (68) is locally exponentially stable.

VI. CONCLUSION

A linear unstable reaction-diffusion equation has been considered in this paper. Both boundary control and in-domain control cases have been considered. For this control problem, saturated feedback control laws have been designed so that the origin is a locally asymptotically stable equilibrium. The region of attraction has been estimated by an appropriate Lyapunov function and LMI technique. The interest and the efficiency of our approach have been illustrated by means of numerical simulations.

This work leaves several questions open. In particular it could be useful to consider other classes of Lyapunov functions than the ones considered in this work, and to compare the associated estimations of region of attraction. Moreover, it could be interesting to use this work for the estimation of the region of attraction in presence of disturbance and to study local input-to-state stability (as presented in Remark 5). Finally, let us note that extension of our method to nonlinear systems, other boundary conditions and control input may be considered as a future work.

  The above computations have been performed for the case when η = a. If a = η, then it holds that |w j (t)| e -ηt |w j (0)| + M |b j | K te -ηt |z(0)|. Now, using the inequality (a + b) 2 2(a 2 + b 2 ) for any (a, b) ∈ IR 2 , and the square summability of |w j (0)| and |b jk |, k = 1, . . . , m, it follows that ∞ j=n+1 |w j (t)| 2 decays exponentially as well.

  g. [4, Chapter 2, p. 7]): Lemma 2: Let A ∈ IR n×n , B ∈ IR m×n , C ∈ IR m×m and let M := A B B C . If C is positive definite, then M is positive semidefinite if and only if its Schur complement M/C := A -B C -1 B is positive semidefinite.

M 1

 1 := (A + BK) P + P (A + BK) P B -(DC) (P B) -

Fig. 1 .

 1 Fig. 1. Region of attraction of (13) (in blue) for the choice (47), (49), computed via the LMI technique. Trajectories of (13) are computed by direct solution of the ODE, trajectories attracted to the origin are in black, diverging trajectories are in red.

Fig. 2 .

 2 Fig. 2. Region of attraction of (13) (in blue) for the choice (50), (51), computed via the LMI technique. Trajectories of (13) are computed by direct solution of the ODE, trajectories attracted to the origin are in black, diverging trajectories are in red.

)

  where b ∞ denotes the L ∞ norm of the function b. Proof. The first claim follows as the assumption guarantee that rsat(k) = sat(rk) = rk. The claim in (57) is a direct consequence of the triangle inequality as |rsat(k)sat(rk)| |r||sat(k)| + |sat(rk)|.

  b is a scalar function, and d is a row vector function with d(x) ∈ IR 1×n d , x ∈ [0, L].

  d j x d (t), j = 1, 2, . . . , where b and d are defined for x in [0, L] in (70) and b j = b(•), e j (•) L2(0,L) , d j = d(•), e j (•) L2(0,L) , for j = 1, 2, . . .. Let us consider the first n equations with the ODE (69) and rewrite this finite-dimensional system as follows:z (t) = Az(t) + BKz(t) + Bφ(Kz(t)) = Az(t) + Bsat(Kz(t)), (72)where K in IR 1×(n+n d ) is a row vector to be designed, z(t) := (x d (t), ω 1 (t), . . . , ω n (t)) , t 0B := (B d , b 1 , . . . , b n ) ∈ IR (n+n d )×1 ,and A:= A d 0 D Λ ∈ IR (n+n d )×(n+n d ) , where

  Proof. As P 22 > 0, the matrix P a is positive definite if and only if the Schur complement P a /P 22 = aP 11 -P 12 P -1 22 P 21 is positive definite, see Lemma 2. Note that as (P 12 ) = P 21

	3:	Consider a symmetric matrix P =
	P 11 P 12
	P 21 P 22

∈ IR (n1+n2)×(n1+n2) , where P ij are matrices of appropriate dimension and P 11 > 0, P 22 > 0. Define P a := aP 11 P 12 P 21 P 22 for a 0. Then there exists a * 0:

(i) P a * 0, and P a * is not positive definite, (ii) P a > 0 for all a > a * , (iii) for all a ∈ [0, a * ) the matrix P a is not positive semidefinite. and (P -1 22 ) = P -1 22 , we have that (P 12 P -1 22 P 21 ) = P 12 (P -1 22 ) (P 21 ) = P 12 P -1 22 P 21 , and as P 11 is positive definite, P a /P 22 is a symmetric matrix.

APPENDIX

A. Series expansion of solutions

Denoting the saturated nonlinearity in (1) by f : IR + → X, we have for the mild solution of (4) that

where T (t) is the strongly continuous semigroup generated by A. Here the integral on the right-hand side is well-defined as the Bochner integral, see [START_REF] Jacob | Linear Port-Hamiltonian Systems on Infinite-Dimensional Spaces[END_REF]Example A.1.13].

Let (e j ) j 1 be the Hilbert basis of X given by the eigenfunctions of A. Then we can define Here the integral on the right is a standard Lebesque integral and so w j solves an integral equation. Thus w j is absolutely continuous and satisfies, for almost all t, the Carathéodory equation

This justifies the consideration of (6) as an equivalent system for (4).

B. Compactness of the resolvent for Sturm-Liouville operators

The following discussion summarizes some results from [START_REF] Savchuk | Sturm-Liouville operators with distribution potentials[END_REF], which provide the necessary arguments to show the compactness of the resolvent of the operator A introduced in Section II.

Here we use the following notation. All function spaces are considered on the interval [0, L]. The Sobolev space W k p is the space of L p -functions (1 p < ∞) such that the function is ktimes weakly differentiable and the corresponding derivatives are again in L p . Note in particular that W k 2 = H k . For negative indices, we set W -1 2 := W 1 0,2 * (the dual to W 1 0,2 = H 1 0 ). Recall that an operator A ∈ L(X) is said to be compact, if A maps bounded sets into precompact sets. For a densely defined linear operator (A, D(A)) : X → X with a nonempty resolvent set ρ(A), it is an easy consequence of the resolvent identity that the resolvent R λ (A) is compact for some λ in the resolvent set, if and only if it is compact on the entire resolvent set, see [START_REF] Kato | Perturbation Theory for Linear Operators[END_REF]Theorem III.6.29].

Definition 3: We say that a closed densely defined linear operator (A, D(A)) : X → X has a compact resolvent, if there exists a λ ∈ ρ(A) so that R λ (A) is compact.

1) Some results from [START_REF] Savchuk | Sturm-Liouville operators with distribution potentials[END_REF]: We note that in [START_REF] Savchuk | Sturm-Liouville operators with distribution potentials[END_REF] the case L = π is considered, which requires some rescaling to use their results.

Let X := L 2 (0, L), q ∈ W -1 2 (0, L), and define (for y ∈ W 1 1 (0, L)) the quasiderivative

where

For q ∈ W -1 2 (0, L) it holds that Q ∈ X, and for q ∈ X it holds that Q ∈ W 1 2 (0, L). Let q ∈ L 1 (0, L) be given. Consider the formal Sturm-Liouville operator SL : X → X defined by SL(y) := -d 2 y dx 2 + q(x)y(x), x ∈ (0, L),

where (0, L) ⊂ IR. In order to fully define the operator SL, we have to introduce domain of definition. Following [34, Section 1.1] we define the maximal operator L M , defined by

D(L M ) :={y : y, y [1] ∈ W 1 1 (0, L), SL(y) ∈ X}.

The following result has been shown in [START_REF] Savchuk | Sturm-Liouville operators with distribution potentials[END_REF]Theorem 1.5]: Theorem 2: Let the operator A be the restriction of L M to the domain

where for j = 1, 2 it holds that U j (y) = a j1 y(0) + a j2 y [1] (0) + b j1 y(L) + a j2 y [1] (L), (77

where a j1 , a j2 , b j1 , a j2 are real numbers, for j = 1, 2.

Let J αβ be the determinant of the α-th and β-th column of a matrix

Then the operator A has a nonempty resolvent set ρ(A), has a compact resolvent and discrete spectrum, if one of the following conditions holds: (i) J 42 = 0, (ii) J 42 = 0, J 14 + J 32 = 0, (iii) J 42 = J 14 = J 32 = 0, J 12 + J 34 = 0, J 13 = 0.

Remark 13: Compactness of the resolvent of A is not mentioned in the formulation of [START_REF] Savchuk | Sturm-Liouville operators with distribution potentials[END_REF]Theorem 1.5], but its compactness was shown in the proof.

• Remark 14: If one of conditions (i)-(iii) holds, then the boundary conditions (77) are called Birkhoff-regular.

• 2) Application to the system in Section II: Consider the operator A as in Section II, i.e.

Let us show that Theorem 2 can be used for our operator A. We need the following lemma, see [8, Exercise 4 at p. 306].

Lemma 5: Let p ∈ [1, +∞). Then f ∈ W 1 p (0, L) if and only if f is equal a.e. to an absolutely continuous function, the derivative f exists a.e. and is an element of L p (0, L).

First we show Lemma 6: The domain D(A) of the operator (79) has the form (76) for U 1 (y) = y(0) and U 2 (y) = y(L).

Proof. Let AC((0, L), IR) be the space of absolutely continuous functions from (0, L) to IR. In view of Lemma 5 the set D(A) can be rewritten as

(80) Now consider the domain defined in (76). As we restrict our attention to the case when q ∈ X, for any y ∈ W 1 1 (0, L) it holds that Q(•)y(•) ∈ W 1 1 (0, L), and thus y [1] ∈ W 1 1 (0, L) if and only if dy dx ∈ W 1 1 (0, L). Hence the domain D(L M ) can be equivalently written as

Using again Lemma 5, we see that

Furthermore, as f ∈ D(L M ) is absolutely continuous on (0, L), and q ∈ X, it holds that qf ∈ X and it holds that SL(f ) ∈ X if and only if