Local stabilization of an unstable parabolic equation via saturated controls
Andrii Mironchenko Christophe Prieur Fabian Wirth
math.OC math.AP math.DS

 






fig/

  







matrix,positioning,decorations.pathreplacing




propertyProperty
nnassumptionAssumption
assumption
nntheoremTheorem
theorem
nncorollaryCorollary
corollary
nndefinitionDefinition
definition
nnpropositionProposition
proposition
nnproblemProblem
problem
nnlemmaLemma
lemma
nnremarkRemark
remark ∘
nnexampleExample
example ⋆




































I-0.21emR
I-0.21emN
ℤ
 






















Local stabilization of an unstable parabolic equation via saturated controls
    Andrii Mironchenko, Christophe Prieur and Fabian Wirth

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 <cit.>.


    December 30, 2023
==========================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================









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.






PDE control, reaction-diffusion equation, saturated control, stabilization, attraction region.







§ INTRODUCTION





In applications of control technology, physical inputs (like force,
torque, thrust, stroke, etc.)  are often limited in size <cit.>. 
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. <cit.>
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 one-dimensional 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 <cit.> a backstepping approach is
applied to design a boundary delayed feedback control for a heat equation
(see <cit.> for further results using the same
design methods). In <cit.> a stable heat
partial differential equation (PDE) is controlled by means of a delayed
bounded linear control operator (see also <cit.> for
the semilinear case). 
For the case of unbounded control operators, see <cit.> for the computation of delayed control to stabilize the reaction-diffusion equation. 

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
<cit.>, where compact and bounded control operators are
considered. Another notable early reference is <cit.> 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 <cit.> for a discussion. 

The definition of the saturation map used in the present paper is aligned
to that used in
<cit.> where Lyapunov methods are shown to be useful for the
stability analysis of wave equations subject to saturated inputs. See also
<cit.>, where systems in Hilbert spaces with applications to the Korteweg–de 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 <cit.>, 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 finite-dimensional systems (see in particular <cit.>). 


We would like to mention the related research on model predictive control of parabolic PDEs under constraints on the state and the input  <cit.>, and on Lyapunov-based control of parabolic systems by controls of a bounded magnitude <cit.>.
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 <cit.> stabilizing backstepping-based
boundary controllers for coupled heat-ODE systems with time-varying state delays in the presence of actuator saturation are developed, and in <cit.> boundary stabilization of a nonlinear Schrödinger equation with state delay and bounded internal disturbance
is performed. Both in <cit.> 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 finite-dimensional unstable part and infinite-dimensional stable dynamics, which is a classical tool used in PDE control <cit.>.
To stabilize the unstable finite-dimensional part using a saturated control we apply techniques, which are well-known for ODE systems (see e.g.  <cit.> for systems with non-delayed input, and <cit.> for systems with a delayed input). 


Using Lyapunov functions, we derive in Section <ref> linear matrix inequalities
(LMIs) whose solutions provide estimates of the region of attraction
of the closed-loop finite-dimensional system, see <cit.> 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 <cit.>.
Then we show how asymptotic stability  and estimates for the region of attraction can
be obtained for the reaction-diffusion equation in closed-loop with the nonlinear saturated control. 

We show that our results can be extended in a variety of different directions.
In Section <ref> we demonstrate, how dynamic controllers can be used to enlarge the region of attraction. In Section <ref>
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 <ref> we show
  how the methods can be applied to unbounded scalar control operators, designed as the output of a finite-dimensional 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. <cit.>), see also Remark <ref>.






The remaining part of this paper is organized as follows. We first
introduce the reaction-diffusion equation with bounded input operator in
Section <ref>. The saturation function is defined and the spectral decomposition is provided. 
A saturated feedback is designed in Section
<ref> 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 <ref> illustrate our design method and the obtained estimates of the region of attraction depending on the saturation level.  
Other saturations are considered in Section
<ref>. In Section
<ref>, 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 <ref>.







Notation: The set of nonnegative reals we denote by _+. The Euclidean norm on ^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
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  we denote the nonnegative integers. For convenience,
^*:=∖{0}.  For k ∈, 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_0^k(0,L) 
(the k-times
continuously differentiable functions with compact support in (0,L))
 in
the norm of H^k(0,L).










§ PROBLEM FORMULATION



We consider the stabilization problem of a one-dimensional linear reaction-diffusion equation by means of a
distributed control u:_+→^m. Let L>0. We are given m functions
b_k: [0,L] →, k=1,…,m which describe at which places the control input u_k∈
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)
                  + ∑_k=1^m b_k(x) (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).








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  is a component-wise saturation function, that is, for all k=1,…, m and for any v ∈^m,

    (v)_k
    :=
    
    v_k     if   |v_k|≤ℓ, 
    ℓ/|v_k|v_k     if  | v_k |≥ℓ,

where ℓ>0 is the given level of the saturation, which is assumed to be uniform with respect to the index k.


 
Systems of the form (<ref>) also occur in the problem of stabilizing the linear heat equation by means of boundary control subject to delays or saturation, see e.g. <cit.> as well as Section <ref> below.


For the design of a stabilizing feedback, we use the well-known eigenfunction decomposition method.
Define 

    A=∂_xx+c(·)id: X → X

with domain D(A)=H^2(0,L)∩ H^1_0(0,L). 
Then the control system (<ref>) takes the form

    w_t(t,·)=Aw(t,·)+ ∑_k=1^m  b_k (u_k(t)).


We note that A is selfadjoint and has compact resolvent, see Appendix <ref>. 
Hence, the spectrum of A consists of only isolated eigenvalues with finite multiplicity, 
see <cit.>. 
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  λ_jj→ +∞⟶-∞

and that e_j(·)∈ D(A) for every j≥ 1.

We consider (mild) solutions of the system (<ref>) (see <cit.>), 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 (<ref>) can be
expanded as a series in the eigenfunctions e_j(·), convergent in
H_0^1(0,L),

    w(t,·)   =∑_j=1^∞w_j(t)e_j(·),
    
       w_j(t)   :=⟨ w(t,·),e_j(·)⟩_L_2(0,L), j∈^*.

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(·)⟩_L_2(0,L), j ∈^*.

As discussed in Appendix <ref>, (<ref>) is
equivalent to the infinite-dimensional control system

    ẇ_j(t)    = λ_jw_j(t)+ ∑_k=1^m  b_jk(u_k(t)) 
       =   λ_jw_j(t)+  𝐛_j·(u(t)) ,    
    j∈^*,

where · is the scalar product in ^m, (u(t)) ∈^m is the vector with entries (u_k(t))
and 𝐛_j is the row vector with entries b_jk, k=1,…,m.

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

    ∀ j>n:  λ_j<-η<0.


With the matrix notations

    z:=[ w_1;   ⋮; w_n ] ,
    𝐀:=[ λ_1   ⋯   0;   ⋮   ⋱   ⋮;   0   ⋯ λ_n ] , 
    𝐁:=[ b_11    ⋯ b_1m;    ⋮    ⋮; b_n1    ⋯ b_nm ]

the n first equations of 
(<ref>)
form the unstable finite-dimensional control system 

    ż(t)=𝐀 z(t) + 𝐁(u(t)).



 
 (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 (<ref>) 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 <cit.> 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 <cit.>), the methods of this paper apply.






§ ESTIMATION OF THE REGION OF ATTRACTION FOR SATURATED INPUTS






 §.§ 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(·))_i=1^n and π_n be the orthogonal projection onto X_n, that is

    π_n w(t,·):=∑_j=1^nw_j(t)e_j(·).

We define also X_n^ as the orthogonal complement
of X_n in X.
Let ι: ^n → X_n be the isomorphism defined by
ι(e^j)=e_j(·), where (e^j)_j=1,…, n is the canonical
basis of ^n. 
We will use the isometric representation of
L_2(0,L) as ℓ_2(^*, ) obtained by the isomorphism
induced by e_j(·)↦ e^j, where

    ℓ_2(^*, ):={(x_k)_k∈^*∈^^*  :  ∑_k=1^∞|x_k|^2 <∞},

where e^j, j∈^* are the standard basis vectors in ℓ_2(^*,
) 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(^*, ) = ^n ⊕ℓ_2,j>n, where we identify ^n with the sequences with
support in {1, …,n } and ℓ_2,j>n is the set of sequences
in ℓ_2(^*, ) which are 0 in the first n entries. 


Given a linear map K:X_n →^m, consider the feedback 

    u=   K π_nw(·) = K ( ∑_j=1^nw_je_j(·) ) 
    =  ∑_j=1^n w_jKe_j(·) 
    
    =   ∑_j=1^n w_j𝐊_j =  𝐊 z,

where 𝐊_j:=Ke_j(·) ∈^m, j=1,…,n, and where we
use 
the notation from (<ref>) in the final step and set 


𝐊:=(𝐊_1,…, 𝐊_n)∈^m× n.




Hence the system (<ref>) with the feedback (<ref>) is equivalent to the following set of differential equations:

    ẇ_j(t)= λ_j w_j(t) + 𝐛_j·(𝐊z(t)),   j ∈^*.

Using (<ref>), we rewrite the first n equations of (<ref>) as

    ż(t) = 𝐀z(t) + 𝐁(𝐊z(t)).

Now, (<ref>) can be considered as a cascade interconnection of an n-dimensional part, described by the equations
(<ref>) and of an infinite-dimensional part described by the equations

    ẇ_j(t)   =   λ_j w_j(t) + 𝐛_j·(𝐊z(t)),   j ≥ n+1.


Our general assumption is:

 
The pair (𝐀,𝐁) is stabilizable, i.e. there is 𝐊∈^m× n:
𝐀+𝐁𝐊 is Hurwitz.




 
Clearly, the closed-loop system (<ref>) with a matrix 𝐊 as in Assumption <ref> 
is locally asymptotically stable. 

Conversely, if a feedback 𝐊 renders the closed-loop system (<ref>)
locally asymptotically stable, then also

    ż= 𝐀z + 𝐁 u

is locally and hence globally asymptotically stabilized by means of the feedback u(t):=𝐊 z(t). Thus, local asymptotic stability of 
(<ref>) implies that the pair (𝐀,𝐁) is stabilizable.






We note that in the case m=1 the situation simplifies further 
as then (<ref>) 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 (<ref>)
boils down to the exponential stabilization of the finite-dimensional unstable system (<ref>). 
This latter problem will be elaborated in Section <ref>.


 
Assume that 𝐊 is chosen so that 0 is a locally asymptotically
stable fixed point of (<ref>).  We say that S
is a region of attraction of 0 if

	
  (i)  0 ∈ S; 
	
  (ii) for any
z_0 ∈ S the corresponding solution of (<ref>)
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 <cit.>.


 
(<ref>) 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)|≤ Me^-at |z_0|   ∀ t≥ 0

	
  (ii)
B_ε(0)⊂ S and S is a region of attraction of (<ref>).

We call (<ref>) globally exponentially stable in 0 if  (<ref>) holds for some a,M>0 and all z_0 ∈^n.

Definitions <ref> and <ref> can be stated analogously for system (<ref>).

We note that if the maximal region of attraction is not ^n, then the system
cannot be exponentially stable on the maximal region of attraction (for
Lipschitz continuous systems), see <cit.>. Thus by analyzing regions of
attraction with exponential stability we necessarily restrict the region
of attraction.













 
Assume 𝐊 is chosen such that the subsystem (<ref>) is locally
exponentially stable in 0 with region of attraction S ⊂^n.
Then:

	
  (i) system
          (<ref>) is
          locally exponentially stable in 0 with region of attraction
S ×ℓ_2,j>n.
	
  (ii) system (<ref>) with the feedback
          (<ref>) is locally exponentially
          stable in 0 with region of attraction ι(S) × X_n^.


  In addition, for any closed and bounded set G⊂(ι(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 

	(<ref>) with the controller (<ref>) satisfies

    w(t,·)_X ≤ Me^-atw(0,·)_X   ∀ t≥ 0.
 



Pick a compact subset G' of S. Since we assume that
(<ref>) is locally exponentially stable in 0 with
region of attraction S ⊂^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 (<ref>) satisfies

    |z(t;z_0)|≤ Me^-at |z_0|,    t≥ 0.



From equations (<ref>) and
(<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) +  𝐛_j·∫_0^t e^λ_j (t-s)(𝐊z(s)) ds.

From (<ref>) it follows that for all z ∈^n we have

    |(𝐊z)| ≤ |𝐊z| ≤𝐊 |z|.

Also due to the Cauchy-Bunyakovsky-Schwarz inequality we have that, for all j∈,

    |b_jk|=|⟨ b_k(·),e_j(·)⟩_X| ≤b_k_X e_j_X = b_k_X.

Thus, for all j≥ n+1,
we obtain exploiting (<ref>) that

    |w_j(t)| 
    				 ≤     e^-η t |w_j(0)| + |𝐛_j| ∫_0^t e^ - η (t-s) |𝐊z(s)| ds
    ≤     e^-η t |w_j(0)| + |𝐛_j|𝐊∫_0^t e^ - η (t-s) Me^-as |z(0)| ds
    
    
    				 =     e^-η t |w_j(0)| +  |𝐛_j| M𝐊/η - a (e^-at - e^-η t) |z(0)|.

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  |𝐛_j|𝐊 t e^-η t|z(0)|.


Now, using the inequality (a+b)^2≤ 2(a^2 + b^2) for any (a,b)∈^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. 

We now obtain local exponential stability of (<ref>) by choosing G' such that 0∈^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 G⊂(ι(S) × X_n^). Select G'=ι^-1∘π_n (G), then G' is a compact subset of S, the previous computations yield
(<ref>)
for suitable constants M and a and for the superset ι(G')× X_n^ which contains G.





 
It is not hard to see that
(<ref>) is
input-to-state stable (ISS) with respect to the input z. Hence
(<ref>) 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 <ref>.






 
Our feedback controller is robust w.r.t. additive actuator disturbances. 
Let the control input to system (<ref>) be u(t):=𝐊z(t) + d(t), where 
d ∈ PC(_+,^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 ∈ PC(_+,^m) with sup_s≥ 0|d(s)|<r the solution of 

    w_t(t,x)=w_xx(t,x)+c(x)w(t,x)
       + ∑_k=1^m b_k(x) ((𝐊z(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).

satisfies the local input-to-state stability (LISS) estimate

    w(t,·)_X ≤ M_2e^-a_2tw(0,·)_X + γ_2 sup_s≥ 0|d(s)|.
 
This property can be obtained quite easily, since for small enough w and d we have 
(𝐊z(t)+d(t)) = 𝐊z(t)+d(t), and thus
(<ref>) is a linear system with a bounded
disturbance operator d ↦∑_k=1^m b_k(·)d_k, acting on PC(_+,^m), and (<ref>) follows as exponential stability of
(<ref>) without disturbances implies LISS for a system (<ref>)
with disturbances.
It is, however, much harder to estimate a region of attraction of system (<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 <cit.> and references therein.

















 §.§ Estimate of the region of attraction for the finite-dimensional part






































In view of Proposition <ref>, it is important to study the local exponential stability and to estimate the region of attraction of the finite-dimensional system (<ref>).
We perform this task in this section.







Defining the deadzone nonlinearity ϕ:^m→^m by 

    ϕ(u)=  (u)- u,   u∈^m,

where  is defined in (<ref>), the following generalized sector condition holds (see <cit.> for a proof):

 
If for some 𝐂∈ℝ^m× n,  z ∈ℝ^n and j∈{ 1,…, m } it holds that
|((𝐊-𝐂) z)_j| ≤ℓ, then

    ϕ_j( 𝐊z)  (ϕ_j(𝐊z)+(𝐂 z) _j)≤ 0.




Note that other properties exist for saturation yielding to other representations of the saturation maps (see, e.g.,  <cit.>). 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 <ref>, for all diagonal positive definite matrices
𝐃∈ℝ^m× m, for all
𝐂∈ℝ^m× n, and all z ∈ℝ^n such
that  |((𝐊-𝐂) z)_j| ≤ℓ, j=1…, m, we have:




    ϕ( 𝐊z) ^⊤𝐃   (ϕ(𝐊z)+ 𝐂z )≤ 0.


We recall the following well-known Schur complement lemma (see e.g. <cit.>):

 
Let A∈^n× n, B∈^m× n, C∈^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^-1B is positive semidefinite.
 

	

	







The following result is only a small variation of <cit.>. Let us give a full proof using the notation used in this paper, for the sake of completeness.









 
Under Assumption <ref>, let 𝐊∈^m× n be such that 𝐀+𝐁𝐊 is Hurwitz.

Let 
P∈ℝ^n× n be symmetric positive definite, 𝐃∈ℝ^m× m diagonal positive definite and 𝐂∈ℝ^m× n such that 

    M_1:=[[ (𝐀+ 𝐁𝐊)^⊤ P + P (𝐀+ 𝐁𝐊)            P 𝐁 - (𝐃𝐂)^⊤;            (P 𝐁)^⊤ - 𝐃𝐂                    -2 𝐃 ]]
    < 0

and 

    M_2:=[[         P (𝐊-𝐂 ) ^⊤;       𝐊-𝐂   ℓ^2 I_m ]]
    ≥ 0.

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

    𝒜:= {z:   z ^⊤ P z ≤ 1}.

Moreover, in A, the function V_1 defined by V_1(z):= z ^⊤ P z, z∈^n, decreases exponentially fast to 0 along the solutions to (<ref>),
i.e. there is a constant α>0 so that 

    V̇_1 (z) ≤ - α |z| ^2,    z ∈ A.






Let us show that (<ref>), (<ref>) are feasible. As 𝐀+𝐁𝐊 is Hurwitz, there is a symmetric positive definite matrix
P∈ℝ^n× n such that (𝐀+ 𝐁𝐊)^⊤ P + P (𝐀+ 𝐁𝐊) <0 and, multiplying if needed P by a positive constant, we may additionally assume that [[       P     𝐊^⊤;       𝐊 ℓ^2 I_m ]]
≥ 0. Then, letting 𝐂=0∈ℝ^m× n, and picking a diagonal positive definite matrix 𝐃∈ℝ^m× m with a norm sufficiently large, we get diagonal positive definite matrix 𝐃∈ℝ^m× m and a matrix 𝐂∈ℝ^m× n such that (<ref>) and (<ref>) hold. 

Now pick P, 𝐃, 𝐂 such that <ref>) and (<ref>) hold and consider the Lyapunov function candidate

    V_1(z)= z ^⊤ P z,

where P∈^n× n is the symmetric positive definite matrix given
by the assumptions. The time-derivative of V_1 along the solutions to
(<ref>) for z ∈^n is given by

    V̇_1(z) =  z ^⊤ [(𝐀+ 𝐁𝐊)^⊤ P + P (𝐀+ 𝐁𝐊)]   z + 2  z ^⊤ P𝐁ϕ(𝐊 z).


Now assuming |((𝐊-𝐂) z)_j| ≤ℓ for all j=1,…,m, it follows from (<ref>) that 

    V̇_1(z)    ≤     
     z ^⊤ [(𝐀+ 𝐁𝐊)^⊤ P + P (𝐀+ 𝐁𝐊)]   z 
             + 2  z ^⊤ P𝐁ϕ( 𝐊 z) - 2 ϕ( 𝐊z)^⊤𝐃 (ϕ(𝐊z)+𝐂 z)
       =   [[     z; ϕ(𝐊z) ]]^⊤ 
    M_1    
    [[     z; ϕ(𝐊z) ]]
    .








In view of (<ref>), this implies the estimate (<ref>) selecting -α as the maximal eigenvalue of M_1.


It remains to ensure that |((𝐊-𝐂) 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∈^n:   z ^⊤ P  z ≤ 1} is included in the ellipsoid {z∈^n:  |(𝐊-𝐂) z| ≤ℓ} (this implies that |((𝐊-𝐂) z)_j| ≤ℓ holds for all j=1, …, m). This inclusion is equivalent to the inclusion of {z∈^n:   z ^⊤ P  z ≤ 1} in the set {z∈^n:  z ^⊤ (𝐊-𝐂)^⊤ (𝐊-𝐂) z ≤ℓ^2}
which is again equivalent to the matrix inequality P-(𝐊-𝐂 ) ^⊤1/ℓ^2 (𝐊-𝐂 )≥ 0. 
As ℓ^2 I_m is positive definite, 
Lemma <ref>

ensures that this latter matrix inequality is equivalent to (<ref>).




















 







We note that the formulation of
Proposition <ref>
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 V̇_1 (z) ≤ - α |z|
^2, so that by a continuity and compactness argument it follows that
V̇_1 (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 ∈ℝ^ n, it holds that
|((𝐊-𝐂) z)_j| ≤ℓ, j=1,…,m. This is
equivalent to 𝐊=𝐂. We thus obtain the following corollary:


 
If there exist a symmetric positive definite matrix P∈ℝ^n× n
and a diagonal positive matrix 𝐃∈ℝ^m× m such that  (<ref>)
holds with 𝐂:=𝐊, 





then the finite-dimensional system (<ref>) is globally exponentially stable in 0. Moreover the Lyapunov function V_1 decreases exponentially fast to 0 along the solutions to (<ref>).

In our context, the main interest of
Proposition <ref> lies in
the following consequence for system (<ref>).



    
Consider system (<ref>) with (𝐀,𝐁) stabilizable. Let 𝐊 be such that 
𝐀 + 𝐁𝐊 is Hurwitz and K be as in (<ref>).





Then the closed-loop system with the feedback K

    w_t(t,x)=w_xx(t,x)+c(x)w(t,x) 
                + ∑_j=1^m b_k(x) ((K π_n w(t,·) )_k),   t>0,   x ∈ (0,L), 
        w(t,0)=w(t,L)=0,    t>0, 
        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^⊥.


    






















    Exponential stabilization and the region of attraction are immediate
    consequences of
    Proposition <ref> in
    combination with
    Proposition <ref>. As
    the exponential estimate can be chosen uniformly on the set A,
    the proof of
    Proposition <ref>
    shows that the uniformity holds on all of ( A) ×
    X_n^⊥.








 §.§ Numerical implementation of Proposition <ref>



The inequalities (<ref>) and (<ref>) are bilinear matrix inequalities, since they have 
the bilinear cross term 𝐃𝐂 in the unknowns 𝐃 and 𝐂. The nonlinearity complicates the numerical analysis of (<ref>) and (<ref>). However, they can be reformulated into equivalent linear matrix inequalities (LMIs), which are easier to solve numerically.

Let S= P^-1 and 𝐄 = 𝐃 ^-1. 
Since M_1<0 if and only if 

    [[ S 0; 0 𝐄 ]]^⊤
    M_1
    [[ S 0; 0 𝐄 ]]<0,


we obtain that (<ref>)  is equivalent to

    [[ S (𝐀+ 𝐁𝐊)^⊤ +  (𝐀+ 𝐁𝐊)S             𝐁𝐄  - S 𝐂^⊤;       ( 𝐁𝐄  - S 𝐂^⊤ )^⊤                    -2 𝐄 ]]
    < 0.

Now pre- and post-multiplying (<ref>) by [[   S   0;   0 I_m ]], we get that (<ref>)  is equivalent to 

    [[                     S      S  𝐊 ^⊤- S 𝐂  ^⊤; (S  𝐊 ^⊤- S 𝐂  ^⊤
    )^⊤               ℓ^2 I_m ]]
    ≥ 0.


Letting 𝐘= S 𝐂  ^⊤, we obtain (see also <cit.> in a different context):

 
Existence of a symmetric positive definite matrix
P∈ℝ^n× n, a diagonal positive definite matrix 𝐃∈ℝ^m× m and a matrix 𝐂∈ℝ^m× n such that (<ref>) and (<ref>) hold is equivalent to 
 the existence of a symmetric positive definite matrix
S∈ℝ^n× n, a diagonal positive definite matrix 𝐄∈ℝ^m× m and 𝐘∈ℝ^m× n for which it holds that

    [[ S (𝐀+ 𝐁𝐊)^⊤ +  (𝐀+ 𝐁𝐊)S                 𝐁𝐄  - 𝐘;         ( 𝐁𝐄  - 𝐘^⊤ )^⊤                    -2 𝐄 ]]
    < 0

and 

    [[                    S      S  𝐊 ^⊤-  𝐘  ^⊤; (S  𝐊 ^⊤-  𝐘  ^⊤
    )^⊤              ℓ^2 I_m ]]
    ≥ 0.


This reformulation is important as the constraints (<ref>) and  (<ref>) are linear in the variables S, 𝐄 and 𝐘, in contrast to bilinear inequalities (<ref>) and  (<ref>).

Finally, no structure is imposed on both the matrix 𝐂∈ℝ^m× n and on the matrix 𝐘∈ℝ^m× n, and we let 𝐂= (S^-1𝐘)  ^⊤ from the knowledge of S and 𝐘 (and respectively we may let 𝐘= P^-1𝐂  ^⊤ from the knowledge of P and 𝐂).


















 

















 





 













 

















 §.§ 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.

 
Consider a symmetric matrix P = 
[ P_11 P_12; P_21 P_22 ]∈^(n_1+n_2)×(n_1+n_2), where P_ij are matrices of appropriate dimension and P_11>0, P_22>0.
Define
P_a:=
[ a P_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.		



	







 
	
	




    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_22^-1 P_21

    is positive definite, see Lemma <ref>. Note that as (P_12)^⊤=P_21 and (P_22^-1)^⊤ = P_22^-1, we have that

    (P_12 P_22^-1 P_21)^⊤ = P^⊤_12 (P_22^-1)^⊤ (P_21)^⊤ = P_12 P_22^-1 P_21,

		and as P_11 is positive definite, P_a/P_22 is a symmetric matrix.
    		
		Weyl's inequality (see <cit.>) 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

 
Under Assumption <ref>, let 𝐊∈^m× n be such that 𝐀+𝐁𝐊 is Hurwitz.

Let 
P̃∈ℝ^n× n be symmetric positive definite and 𝐂∈ℝ^1× n such that 



 


    M̃_1:=[[ (𝐀+ 𝐁𝐊)^⊤P̃ + P̃ (𝐀+ 𝐁𝐊)                P̃𝐁 - 𝐂^⊤;              (P̃𝐁)^⊤ - 𝐂                       -2 ]]
    < 0.

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

    {z:   z ^⊤P̃ z ≤𝐃^ -1},

where 𝐃 is the minimal real number such that:

    M̃_2:=[[       𝐃P̃ (𝐊-𝐂 ) ^⊤;       𝐊-𝐂       ℓ^2 ]]
    ≥ 0.

Moreover, in this region of attraction, the function V_1 defined by
V_1(z)= z ^⊤ P z := 𝐃  z ^⊤P̃ z, for all z∈^n, decreases exponentially fast to 0 along the solutions to (<ref>),
i.e. there is a constant α>0 so that 

    V̇_1 (z) ≤ - α |z| ^2.





Define P:=𝐃P̃. Substituting this expression into 
(<ref>) and (<ref>) and multiplying it by 𝐃^-1 to obtain the inequalities
(<ref>) and (<ref>).
If there exist P̃ and 𝐂 as in Proposition <ref> so that the LMI (<ref>) holds, then according to Lemma <ref>, one can find 𝐃>0 so that 
(<ref>) holds as well.

The estimate for the region of attraction follows from Proposition <ref>.
We choose the minimal 𝐃 to obtain the largest estimate of a region of attraction (for given P̃ and 𝐂).



 





Lemma <ref> can be helpful in finding the
optimal 𝐃 as it shows that the problem is to find the unique
root of a piecewise analytic function.









 §.§ Enlarging the region of attraction using a dynamic controller



Consider a
dynamic controller for system (<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. <cit.> for a simple 2× 2 example). 

To be more specific, we consider 
an additional finite-dimensional state z_c in ^n_c (for a given integer n_c) and design matrices 𝐊_1, 𝐊_2, 𝐀_c and 𝐁_c 
of appropriate dimensions so that the (estimation of the) region of attraction of 

    [                       ż                       = 𝐀z + 𝐁 (𝐊_1 z+𝐊_2 z _c);                     ż_c                       =         𝐀_1z_c + 𝐀_2  z ]

is larger (in the z-direction) than the estimate provided by Proposition <ref> for (<ref>). We rewrite the dynamics (<ref>) by 

    [                  Ż                  = 𝐀̅Z +  𝐁̅ (𝐊̅ Z  ) ]

where 

    𝐀̅= [   𝐀   0; 𝐀_2 𝐀_1 ], 𝐁̅=[ 𝐁; 0 ],  𝐊̅=[ 𝐊_1 𝐊_2 ],

with 𝐀̅∈^(n+n_c)× (n+n_c), 𝐁̅∈^(n+n_c)× m and 𝐊̅∈^m× (n+n_c).




 
For given symmetric positive definite matrices P ∈^n× n and P̅∈^(n+n_c)× (n+n_c), the inclusion
of the ellipsoid {z∈^n :   z^⊤ Pz ≤ 1} in
the projection (onto ^n) of the ellipsoid {Z ∈^n+n_c:   Z^⊤P̅Z ≤ 1} is equivalent to 

    [ I_n   0 ]P̅[ I_n;   0 ] - P ≤ 0 .






Using this remark  and applying Proposition <ref> to system (<ref>), we have the following:



 
Consider system (<ref>) with a stabilizable pair (𝐀,𝐁) ∈^n× n×^n× m, and matrices 𝐀̅, 𝐁̅ and 𝐊̅ defined in (<ref>).





Choose a symmetric positive definite matrix P̅∈ℝ^(n+n_c)× (n+n_c), a diagonal positive matrix 𝐃̅∈ℝ^m× m and 𝐂̅∈ℝ^m× (n+n_c) such that (<ref>) holds and

    [[ (𝐀̅+ 𝐁̅𝐊̅)^⊤P̅ + P̅ (𝐀̅+ 𝐁̅𝐊̅)                P̅𝐁̅ - (𝐃̅𝐂̅)^⊤;                (P̅𝐁̅)^⊤ - 𝐃̅𝐂̅                          -2 𝐃̅ ]]
    < 0,




    [[          P̅ (𝐊̅-𝐂̅ ) ^⊤;       𝐊̅-𝐂̅     ℓ^2 I_m ]]
    ≥ 0.

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


A_P̅ :={Z ∈^n + n_c :   Z ^⊤P̅ Z ≤ 1}. 



Moreover, the projection of the ellipsoid {Z :   Z ^⊤P̅ Z ≤
1} onto ^n is larger than the region of attraction A given
by Proposition <ref>
in (<ref>) for  system (<ref>).


Finally, in A_P̅, the function V_1 defined by 

    V_1(Z)= Z ^⊤P̅ Z,     Z∈^n+n_c,
 
decreases exponentially fast to 0 along the solutions to (<ref>),
i.e. there is a constant α>0 so that for all Z ∈ A_P̅

    V̇_̇1̇(Z) ≤ - α |Z| ^2.













 






 




























As in Section <ref>, the inequalities (<ref>) and (<ref>) can be reformulated as linear matrix inequalities.




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








Now under the assumptions of Proposition <ref>, we provide a
Lyapunov function V:ℓ_2(^*,)→_+ for the system (<ref>). To this end,
we use the Lyapunov function P provided by Proposition
<ref> for the
finite-dimensional subsystem (<ref>).

 
Consider system (<ref>) and
assume 𝐊 is chosen such that the subsystem (<ref>) is locally
exponentially stable in 0. Assume further that P ∈^n × n is
symmetric positive definite and such that for Ṽ(z) = z^⊤ P z, z ∈^n we have
Ṽ̇ (z)≤ - α | z |^2 on the set A:={ z ∈^n: Ṽ(z) ≤ 1 }
along the solutions of (<ref>). 
Then there exist γ, C > 0 such that for Q: ℓ_2 (^*, ) → defined by Q(w) := ∑_j=n+1^∞ w_j^2 the function

    V(w) := Ṽ(π_n w) + γ Q(w) = z^⊤ P z + γ Q(w)

(with the identification π_n w =z) satisfies for all w ∈ A×ℓ_2,j>n that along the solutions of (<ref>)

    V̇(w)≤ - C  w ^2_ℓ_2.

In particular, it follows that S ×ℓ_2,j>n is a region of attraction of
0 for system (<ref>).



   We write w∈ℓ_2 (^*, ) as w=z + z^, where z∈^n,
   z^∈ℓ_2,j>n.
	Here we identify ^n with the sequences with
support in {1, …,n } and ℓ_2,j>n is the set of sequences
in ℓ_2(^*, ) which are 0 in the first n entries. 
	
	This decomposition is unique.
   Due to (<ref>), we have along the solutions to (<ref>)

    Q̇(w) 
     =     2∑_j=n+1^∞ w_j(t)ẇ_j(t)
    
    =    2∑_j=n+1^∞ w_j(t)(λ_j w_j(t) + 𝐛_j ·(𝐊z(t)))
    ≤    -2η∑_j=n+1^∞ w^2_j(t)  + 2∑_j=n+1^∞ |w_j(t)| |𝐛_j| |(𝐊z(t))|
    ≤    -2η∑_j=n+1^∞ w^2_j(t)  + 2∑_j=n+1^∞ |w_j(t)| |𝐛_j| 𝐊|z(t)|.

Cauchy-Bunyakovsky-Schwarz inequality implies that

    Q̇(w) ≤    -2ηz^(t)_ℓ_2^2  + 2 z^(t)_ℓ_2b^(t)_ℓ_2𝐊|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𝐊)  z^(t)^2_ℓ_2 +1/κb_ℓ_2𝐊 | z(t)|^2.

Therefore for any choice of γ in (<ref>), we have by an application of Proposition
<ref> along the solutions to
(<ref>) that

    V̇(w) ≤    -(α - γ/κb^_ℓ_2𝐊) | z(t)|^2
               - γ(2η-κb^_ℓ_2𝐊)  z^(t)^2_ℓ_2.

Therefore selecting first κ>0 such that 2η-κb^_ℓ_2𝐊>0 and then selecting γ>0 such that 
α - γ/κb^_ℓ_2𝐊>0, we get the existence of C>0 such that

    [                        V̇(w)                            ≤ -C | z(t)|^2 - C z^(t)^2_ℓ_2 ]

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 (<ref>), we obtain that V is a Lyapunov function for (<ref>), with a guaranteed region of attraction containing
A×ℓ_2,j>n. 





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

 




§ NUMERICAL EXPERIMENTS



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

According to <cit.> the eigenvalues of the operator

    A:=∂_xx+c id: X → X

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∈^*,

and the eigenfunctions e_j, j∈^* of (A,D(A)), which form a basis of L_2(0,1) are given by

    e_j(x):= ( 2/L)^1/2sinjπ x/L,   j∈^*,    x∈(0,L).





















Next we estimate the region of attraction of system (<ref>) for the following choice of parameters:

    c(x) ≡ 10,   L = 2,   ℓ = 2,   b = e_1 + e_2,

where the e_j are defined in (<ref>).
For these parameters there are only 2 unstable eigenvalues. Using the eigenfunction decomposition of solutions of (<ref>) as in Section <ref>, we see that the matrices 𝐀, 𝐁 defined in (<ref>) have the form

    𝐀≈[ 7.5325989         0;         0 0.1303956 ]
    ,
        𝐁 = 
    [ 1; 1 ]
    .

As the diagonal entries in the matrix 𝐀 are distinct, and all
the components of 𝐁 are nonzero, the 
pair (𝐀,𝐁) is stabilizable.



Different choices of the matrix 𝐊 for the stabilizing feedback
u(t)=𝐊z(t) lead to different attraction rates and different
regions of attraction. We demonstrate this with two examples.












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



First we choose the matrix 𝐊 so that σ(𝐀 + 𝐁𝐊) = {-1}, which results in 

    𝐊≈[ -9.835618 0.1726235 ]
    .

To estimate a region of attraction of (<ref>) by Proposition <ref>, we proceed in two steps.

	
  (i) First we solve the inequality
          (<ref>) together with the additional constraints: P̃-P̃^⊤=0 and P̃>0.
Additionally, we impose an optimality condition for 𝐂:

    (𝐊-𝐂)·(𝐊-𝐂)^⊤→min,

where · is the scalar product of vectors.

	
  (ii) The idea behind (<ref>)
 is to minimize the off-diagonal elements of the matrix M̃_2, which gives us at the second step 
the possibility to find large 𝐃 (optimal for the given P̃, 𝐂) satisfying the bilinear matrix inequality (<ref>).




This algorithm is implemented in Scilab. For solution of the LMI (<ref>) the LMITOOL package has been used.
The resulting matrices P̃, 𝐂 and the real number
𝐃 are (approximately):


    P̃ =   [  2.1277468 -0.0655569; -0.0655569  0.0243008 ],
    𝐂 =    (  -2.0635579,  0.0844904)
    ,  𝐃 =  7.359375.









Figure <ref> shows an elliptic region of attraction (<ref>),
subject to P̃, 𝐂, 𝐃 given by (<ref>) (in
blue).
Furthermore, in the same figure some trajectories are depicted
obtained through direct simulation of (<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 (<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 <ref>, we obtain an approximation of the region of attraction of the PDE (<ref>) subject to the feedback (<ref>).












 
(Importance of the optimality conditions) Different solutions
P̃, 𝐂, 𝐃 of the matrix inequalities
(<ref>), (<ref>) lead to very different
estimates of a region of attraction (<ref>) of the
model (<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
(<ref>).  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
(<ref>), (<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 <cit.> for more details and a complete introduction of such optimization problems.




 
(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̃, 𝐂, 𝐃 via LMIs: 0.0166561 seconds.
	
  * Plotting the obtained region: 0.0071209 seconds.
	
  * Time for solving the ODE (<ref>) 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.





 §.§ Choice 2: Putting the poles to {-0.1,-0.2}



Now let us choose the matrix 𝐊 so that σ(𝐀 + 𝐁𝐊) = {-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 𝐊 is:

    𝐊 = 
    [ -7.9732782  0.0102837 ]
    .


As in the previous simulation, we solve the LMI (<ref>) subject to the additional optimality condition
(<ref>). The corresponding matrices P̃, 𝐂, 𝐃 are:


    P̃ =   [  0.3108695 -0.0054849; -0.0054849   0.000195 ],	
    𝐂 =    (-0.3053879, 0.0054754)
    ,  𝐃 =  90.625.



As we see, with the choice of the stabilizing feedback (<ref>), 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 (<ref>) leads to a
better estimate of the region of attraction, in comparison to the
situation in Section <ref>. This can be seen by comparing the Figures <ref>
and <ref>.



 
(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̃, 𝐂, 𝐃 via LMIs: 0.018131 seconds.
	
  * Plotting the obtained region: 0.0129341 seconds.
	
  * Time for solving the ODE (<ref>) 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 (<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.













	
	






 


§ EXTENSIONS





 §.§ 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 <cit.>.




  §.§.§ 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 ^n we may consider

    _·(u) :=
    u    ,  if u≤ℓ
    ℓu/ u   ,  if  u ≥ℓ
     = 
    ℓmin{1/ℓ,1/u }u.

Another physically motivated saturation map is the following (pointwise) L_∞  saturation map, defined for all x ∈ (0,L) by

    _∞(u)(x)
    :=   
    
    u(x)    ,  if  |u(x)|≤ℓ
    ℓu(x)/|u(x)|   ,  if  | u(x) |≥ℓ.
    
    
     =   ℓmin{1/ℓ,1/|u(x)|}u(x).

In this section, we depart from the saturation model in
(<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)+ ∑_k=1^m (b_k(x)u_k(t)).

In a certain sense, in the equation (<ref>) the whole terms
u_kb_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 <ref>. 

In order to stabilize (<ref>), we are going to use the same
stabilizing control (<ref>). We now aim to
provide an estimate for a region of attraction for (<ref>). 
So we assume that 𝐊 as in (<ref>) is
given and we consider the closed-loop system

    w_t(t,x)= w_xx(t,x)+c(x)w(t,x)+ 
    ∑_k=1^m
     (b_k(x)(𝐊z(t))_k).



Representing this equation in the basis e_j, j∈^*
as in Section <ref>, we obtain the equations for the coordinates

    ẇ_j(t)=   λ_j w_j(t)  + ∑_k=1^m⟨ e_j,
    _∞(b_k(·)(𝐊z(t))_k)⟩
    ,  j∈^*.


Let us state the following simple result that will be instrumental to bound the differences

    Δ( (𝐊z(t))_k, b_k(·)):= b_k(·)(𝐊z(t)_k)-_∞(b_k(·)𝐊z(t)_k).



Let r,k ∈ and denote Δ(r,k) = r(k)- (rk). Then

    |k| ≤ℓ |rk | ≤ℓ  ⇒  Δ(r,k) = 0.


    |Δ(r,k) | ≤ℓ ( 1 + |r|).



Let k ∈ and b ∈ L_2(0,L) be given. Then

    |k| ≤ℓb(·) k _∞≤ℓ  ⇒  Δ (b(·),k) ≡ 0

Moreover

    b(·)(k)-_∞(b(·)k)_2≤ℓ ( 1+  |b|_2).

If |k| ≤ℓ and χ is the characteristic function of the set U:={ x ∈ [ 0, L ]: |kb(x) | > ℓ} then 

    b(·)(k)-_∞(b(·)k)_2≤ℓ ( χ+  |b χ|_2).

If b(·) is essentially bounded, then

    b(·)(k)-_∞(b(·)k)_∞≤ℓ ( 1+  b_∞),

where b_∞ denotes the L_∞ norm of the function b. 



The first claim follows as the assumption guarantee that r(k) =
(rk) = rk. The claim in (<ref>) is a direct consequence of the
triangle inequality as

    |r(k)- (rk)| ≤ |r||(k)| + |(rk)|.

The remaining claims follow immediately by applying the pointwise
estimates (<ref>) and (<ref>). The claim (<ref>)
follows by applying (<ref>) to bχ, after noting that
the complement of U does not contribute to the norm of the left hand side.

























    
    Consider system (<ref>) and assume all functions
    b_k(·) ∈ L_∞([ 0, L ]), k=1,…,m.  		
Consider also the
    associated system (<ref>) with a stabilizable pair (𝐀,𝐁) ∈^n× n×^n× m and let 𝐊∈^m× n be such that 𝐀+𝐁𝐊 is Hurwitz.
Choose a symmetric positive definite matrix
P∈ℝ^n× n as in Proposition <ref>.
		
    
    
    
    
    
    
    
    

    
		Consider the Lyapunov function candidate V defined in
    (<ref>).



   
   
	




Then there exist α, β >0 such that if for w = z + z^⊥∈ℓ_2 (^*, ) we have z∈ A_β := {z:   z^⊤ Pz≤β} then

    V̇(w) ≤ -α/2 |w|^2.
    
In particular, {z:   z^⊤ Pz≤β}× X_n^ is a region
of attraction for system (<ref>). 



First we rewrite (<ref>) by adding and subtracting the term
∑ b_k(x) (u_k(t)) to obtain

    w_t(t,x)=w_xx(t,x)+c(x)w(t,x)+ ∑_k=1^m b_k(x) (u_k(t)) 
    
     + ∑_k=1^m((b_k(x)u_k(t)) - b_k(x)(u_k(t)) ).


Define 

    Δ(t):= ∑_k=1^m _∞(b_k(·)(𝐊z(t))_k) - b_k(·)  (𝐊z(t)_k).

For j ∈^* we define
    y_j(t):=  ⟨ e_j, Δ(t) ⟩, and let  y(t)∈^n be the vector with components y_j(t), for j=1,…,n.

Considering the Lyapunov function V defined in (<ref>), we compute its time-derivative along the solutions to 
(<ref>).
Using (<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+ 2 z(t)^⊤ P y(t)  
    
    
    + 2 γ∑_j=n+1^∞ w_j(t) y_j(t).









Using (<ref>) from Lemma <ref>, we obtain along the solutions to (<ref>) and as long as z^⊤ Pz ≤ 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 (<ref>), there exists β>0, such that for all z in {z:   z^⊤ Pz≤β},  Δ(t)_∞ ^2≤C/4(λ_maxκ   +γκ' )| z(t)|^2. We get, for all solutions to (<ref>), as long as z(t) is in {z:   z^⊤ Pz≤β}, 

    V̇ (w) ≤ -C/2 |z(t)|^2-C/2 z^(t)^2_ℓ_2

Moreover, we have, along the solutions to (<ref>),

    ẇ_j(t)=λ_jw_j(t)+⟨( b K z(t)),e_j⟩ ,     j=1,  2, …

therefore, considering V_1 as previously defined, we may check that
(<ref>) holds
following the same computation as for V̇ along the solutions to
(<ref>), and using the fact the dynamics
(<ref>) in X_n does not depend on the
component in X_n^. 
Therefore the set {z:   z^⊤ Pz≤β}× X_n^ is invariant along the dynamics to (<ref>). With (<ref>), we get that {z:   z^⊤ Pz≤β}× X_n^ is a region of attraction and V is a Lyapunov function.    
































































 §.§ 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),
    
    y(t,0)=    0 ,  y(t,L)= y_d ,   t≥ 0,


where y_d is the (scalar) output of finite-dimensional dynamical system
given by



    ẋ_d =    A_d x_ d + B_d   (u(t)), 
    
    y_d =    C_d x_d.











Here x_d in ^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 (<ref>) and the ODE
(<ref>) that is subject to a saturation map. 
Inspired by <cit.>, 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:

    w_t(t,x)    =    y_t(t,x)- x/Lẏ_d(t)
       =    y_xx(t,x) + c(x)  y(t,x) - x/L C_d ẋ_d(t)
       =    w_xx(t,x) + c(x) (w(t,x) + x/L y_d(t)) 
              - x/L C_d (A_d x_d(t) + B_d (u(t)))
       =    w_xx(t,x) + c(x) (w(t,x) + x/L C_d x_d(t)) 
              - x/L C_d (A_d x_d(t) + B_d (u(t)))
       =    w_xx(t,x) + c(x) w(t,x) 
              + (c(x) x/L C_d - x/L C_d A_d_=:d(x)) x_d(t)
              +(- x/L C_d B_d_=:b(x)) (u(t)).


Please note that b is a scalar function, and d is a row vector function with d(x) ∈^1× n_d, x∈[0,L].

The boundary conditions for the variable w take the form:

    w(t,0)=w(t,L)=0 ,   t≥ 0.

The heat equation (<ref>),
(<ref>) has to be analyzed along with the ODE (<ref>).

Performing similar computations as in Section <ref> for the PDE (<ref>) and, using the same notation for w_j and λ_j, we get

    ẇ_ j (t) = λ_ j w_ j(t) + b_j  (u(t)) + d_ j x_d (t), j=1,  2, …,

where b and d are defined for x in [0,L] in (<ref>) and b_j=⟨ b(·),e_j(·)⟩_L_2(0,L), d_j=⟨ d(·),e_j(·)⟩_L_2(0,L), for j=1,  2, …. 

Let us consider the first n equations with the ODE (<ref>) and rewrite this finite-dimensional system as follows:

    z'(t)    =   𝐀z(t) + 𝐁𝐊z(t) + 𝐁ϕ(𝐊z(t))
       =   𝐀 z(t) + 𝐁(𝐊z(t)),

where 𝐊 in ^1× (n+n_d) is a row vector to be designed,

    z(t)   :=   (x_d^⊤(t),ω_1(t),…,ω_n(t))^⊤,    t≥ 0
    𝐁   :=    (B_d^⊤,b_1,…,b_n)^⊤∈^(n+n_d)× 1,

and 
𝐀:=[ A_d   0;   D   Λ ]∈^(n+n_d)× (n+n_d), where

    D:=[   d_11   d_12      ⋯ d_1n_d;      ⋮      ⋮      ⋮      ⋮;   d_n1   d_n2      ⋯ d_nn_d ],
      Λ:=[ λ_1       0;       ⋱    ;   0     λ_n ]
    .









 





 














 
















As the control input is scalar, we can apply
Proposition <ref>
to system (<ref>) instead of system
(<ref>), and get sufficient conditions for the
estimation of the region of attraction of (<ref>). Coming back to the infinite-dimensional systems (<ref>) and
(<ref>), and, applying Proposition <ref>, we get sufficient conditions for an estimation of attraction region of (<ref>):

Consider system (<ref>) with a stabilizable pair (𝐀,𝐁) ∈^n× n×^n× 1 and let 𝐊∈^1× n be such that 𝐀+𝐁𝐊 is Hurwitz.
Pick a symmetric positive definite P̃∈ℝ^(n+n_d)× (n+n_d) and a 𝐂∈ℝ^1× (n+n_d) such that (<ref>) holds.

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

    A:={z:   z ^⊤P̃ z ≤𝐃^ -1},

where 𝐃 is the minimal real number such that (<ref>) holds.


Moreover, 

	
  (i) (<ref>) is locally exponentially stable with a region of attraction
( A) × X_n^⊥,
	
  (ii) (<ref>) is locally exponentially stable.

















 








 








 







	

	








§ 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 <ref>). 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.







abbrv







[
    < g r a p h i c s >
]Andrii Mironchenko
received his MSc at the I.I. Mechnikov Odessa National University in 2008 and his PhD at the University of Bremen in 2012. 
He has held a research position at the University of Würzburg and was a Postdoctoral JSPS fellow at the Kyushu Institute of Technology (2013–2014). 
In 2014 he joined the Faculty of Mathematics and Computer Science at the University of Passau. 
His research interests include infinite-dimensional systems, stability theory, hybrid systems and applications of control theory to biological systems. 




[
    < g r a p h i c s >
]Christophe Prieur
 was born in 1974. He is currently a senior researcher of the CNRS at the Gipsa-lab, Grenoble, France. He is currently a member of the EUCA-CEB, an associate editor for the  AIMS Evolution Equations and Control Theory and IEEE Trans. on Control Systems Technology, a senior editor for the IEEE Control Systems Letters, and an editor for the IMA Journal of Mathematical Control and Information. He was the Program Chair of the 9th IFAC Symposium on Nonlinear Control Systems (NOLCOS 2013) and of the 14th European Control Conference (ECC 2015). His current research interests include nonlinear control theory, hybrid systems, and control of partial differential equations. 
 

[
    < g r a p h i c s >
]Fabian Wirth
received his PhD from the Institute of Dynamical Systems at
the University of Bremen in 1995. He has since held positions at the Centre Automatique et Systèmes of Ecole des Mines, the
Hamilton Institute at NUI Maynooth, Ireland, the University of Würzburg and IBM Research Ireland. 
He now holds the chair for Dynamical Systems at the University
of Passau. His current interests include stability
theory, switched systems and large scale networks with applications to networked systems and in the domain of smart cities.



 



 §.§ Series expansion of solutions



Denoting the saturated nonlinearity
in (<ref>) by f:_+ → X, we have for the mild solution of (<ref>) that

    w(t)=T(t)w(0)+∫_0^tT(t-s)f(s)ds ,

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 <cit.>.

Let (e_j)_j≥ 1 be the Hilbert basis of X given by the eigenfunctions of A. Then we can define

    w_j(t)    :=⟨ w(t), e_j⟩
       =⟨ T(t) w(0), e_j⟩ + ⟨∫_0^tT(t-s)f(s)ds, e_j⟩

and by <cit.> we may interchange the linear map
⟨·, e_j ⟩ with the integral to obtain

    w_j(t)    = ⟨ T(t)w(0),e_j⟩ +∫_0^t⟨ T(t-s)f(s), e_j⟩ ds
       =⟨ w(0), T(t)^∗e_j⟩ + ∫_0^t⟨ f(s), T(t-s)^∗e_j⟩ ds
       =e^λ_jt⟨ w(0), e_j⟩ + ∫_0^te^λ_j(t-s)⟨ f(s), e_j⟩ ds .

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

    ẇ_j(t)    =λ_jw_j(t)+⟨∑_k=1^mb_ksat(u_k(t)), e_j⟩
       =λ_j w_j(t)+∑_k=1^mb_kjsat(u_k(t)) .

This justifies the consideration of (<ref>) as an
equivalent system for (<ref>).



 §.§ Compactness of the resolvent for Sturm-Liouville operators



The following discussion summarizes some results from <cit.>, which
provide the necessary arguments to show the compactness of the resolvent
of the operator
A introduced in Section <ref>.

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
k-times 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 <cit.>.












 
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.






  §.§.§ Some results from <cit.>


We note that in <cit.> 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

    y^[1](x):= dy/dx - Q(x)y(x),

where 

    Q(x):= ∫_0^xq(s)ds.



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^2y/dx^2 + q(x)y(x),    x ∈ (0,L),

where (0,L) ⊂.
In order to fully define the operator SL, we have to introduce its domain of definition.
Following <cit.> we define the maximal operator L_M, defined by



    L_My =    SL y, 
    
    D(L_M):=   {y : y,y^[1]∈ W^1_1(0,L), SL(y) ∈ X}.



The following result has been shown in <cit.>:

 
Let the operator A be the restriction of L_M to the domain

    D(A):={y∈ D(L_M): U_1(y)=U_2(y) = 0},

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),

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 

    [ a_11 a_12 b_11 b_12; a_21 a_22 b_21 b_22 ].

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.
 



 
Compactness of the resolvent of A is not mentioned in the formulation of <cit.>, but 
its compactness was shown in the proof.





 
If one of conditions (i)–(iii) holds, then the boundary conditions (<ref>) are called
Birkhoff-regular.







  §.§.§ Application to the system in Section <ref>



Consider the operator A as in Section <ref>, i.e.



    A   :=   ∂_xx+c(·)id: X → X,
    
    D(A)   =   H^2(0,L)∩ H^1_0(0,L)



Let us show that Theorem <ref> can be used for our operator A.
We need the following lemma, see <cit.>.

 
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

 
The domain D(A) of the operator (<ref>) has the form (<ref>) for U_1(y) = y(0) and U_2(y)=y(L).



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

    D(A)    =   {f∈ AC((0,L),):df/dx∈ AC((0,L),), 
              d^2f/dx^2∈ L_2(0,L), f(0)=f(L)=0}.


Now consider the domain defined in  (<ref>).

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 

    D(L_M)={f : f, df/dx∈ W^1_1(0,L), SL(f) ∈ X}.

Using again Lemma <ref>, we see that

    D(L_M)   =   {f∈ AC((0,L),): df/dx∈ AC((0,L),), 
              d^2f/dx^2∈ L_1(0,L),  SL(f) ∈ X}.

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 d^2f/dx^2∈ X.
As X ⊂ L_1(0,L), we can finally restate D(L_M) as

    D(L_M)   =   {f∈ AC((0,L),): df/dx∈ AC((0,L),),
              d^2f/dx^2∈ L_2(0,L)}.

Using the boundary conditions U_1(y) = y(0) and U_2(y)=y(L), we see that 
the set (<ref>) is precisely (<ref>).


In view of Lemma <ref>, we can use Theorem <ref> 
to study the spectral properties of the Sturm-Liouville operator A.


 
The operator (<ref>) has a compact resolvent.



For the boundary conditions U_1(y) = y(0) and U_2(y)=y(L)
the coefficients a_ij and b_ij from the formulation of Theorem <ref>
have the form:
a_11 = 1, b_21 = 1, and all other entries of a matrix in (<ref>) are zeros.

By item (iii) of Theorem <ref> the boundary conditions are Birkhoff-regular and the operator A has a compact resolvent.















§ TEMPORARY MATERIAL





 §.§ Some questions



	
  * Has one considered the global stabilization of nonlinear unstable systems by means of saturated controllers?
	For example, if the system has an asymptotically stable bounded attractor (for u≡ 0), then such a design would be possible. E.g. for a system like this:

    ẋ = x-x^3 +(u)
	
In this case the problem would be to construct a saturated feedback so that the region of attraction of the closed loop system includes the global attractor for a zero input.

Yes you are right. Our conditions for the local stability analysis may be particularized to conditions for the global stability analysis. Such sufficient conditions for the global stability analysis may be satisfied or not.



  * Just a newbie-question. One of the motivations for the problem which we consider was that there are no results for local stabilization of PDEs by means of the saturated controllers.

However, it seems to me that if a certain system can be locally uniformly asymptotically stabilized by a certain controller u = k(x), where k: X → U can be linear or nonlinear function, satisfying 

    k(x)_U ≤ψ (x_X),

then the controller with saturation ũ= (k(x)) will necessarily locally stabilize the same system.

Indeed, since there is r>0 so that the trajectories of the closed-loop system satisfy the estimate 	

    ϕ(t,x)_X≤β(x_X,t),  x_X≤ r, t≥ 0

for a certain β∈, then there is c>0 so that β(c,0)≤ψ^-1(1).
Hence if x∈ B_c(0), then ϕ(t,x)_X≤ 1 for all t≥ 0, and thus

    k(ϕ(t,x))_U ≤ψ (ψ(t,x)_X) ≤ 1,   t≥ 0,

and since for y∈ U: y_U≤ 1 we have (y) = y, for the states emanating from B_c(0) the saturated controller acts as the controller without the saturation.

Do I miss some important points? Or is the main idea to give a good estimate for the region of attraction?

Yes, you are right, but in that case you have (y) = y, thus all dynamics is linear. You seem to impose a linear dynamics and thus the estimation of region of attraction would be worse than by allowing some saturations of the inputs. In the next computations, the input may saturate and thus I guess that the region of attraction would be larger.










§ STABILIZATION OF THE LINEAR PARABOLIC SYSTEM




Let X_n be the subspace of L_2(0,L) spanned by
e_1(·),…,e_n(·) and π_n be the orthogonal projection onto X_n, that is

    π_n w(t,·):=∑_j=1^nw_j(t)e_j(·).



Consider the following feedback control 

    u(t)=K π_nw(t,·) = K ( ∑_j=1^nw_j(t)e_j(·) ) 
    =  ∑_j=1^n w_j(t)Ke_j(·) 
    =  ∑_j=1^n w_j(t)𝐊_j
    =  𝐊 z(t),

where K:X_n →, 𝐊_j:=Ke_j(·), j=1,…,n and 𝐊:=(𝐊_1,…, 𝐊_n)∈^n is a row vector.

Using this controller for the system (<ref>), we obtain the closed-loop system of the form

    z'(t)= (𝐀 + 𝐁𝐊) z(t).

Furthermore, substituting this control to the system (<ref>), we see that the corresponding closed-loop system is a 
cascade interconnection of the system
(<ref>)
with a system 

    w_j'(t)=λ_jw_j(t)+b_j 𝐊 z(t),     j=n+1,…, ∞.

which is exponentially ISS if the ℓ_2 norm for the state is chosen and z is treated as an input. 

Assuming that there is 𝐊, so that (<ref>) exponentially stabilizes the system (<ref>),
the closed-loop system for (<ref>) will be a cascade interconnection of a UGAS system and ISS system, and hence UGAS.
Since the general theorem on cascade interconnections does not guarantee exponential UGAS, we perform some trivial computations next to ensure that the overall closed-loop system is exponentially UGAS.

Since 𝐀 + 𝐁𝐊 is a Hurwitz matrix, there are M,a>0 so that 

    |z(t)|≤ Me^-at |z(0)|.

Also due to Cauchy-Schwarz inequality we have that for j=1,…, ∞ it holds that

    |b_j|=|⟨ b(·),e_j(·)⟩_L_2(0,L)| ≤b_X e_j_X = b_X.

From equations (<ref>) we derive that for j=n+1,…, ∞ it holds that

    w_j(t) = e^λ_j t w_j(0) + b_j ∫_0^t e^λ_j (t-s)𝐊z(s) ds

and thus

    |w_j(t)|    ≤     e^λ_j t |w_j(0)| + |b_j| ∫_0^t e^λ_j (t-s) |𝐊z(s)| ds
       ≤     e^-η t |w_j(0)| + |b_j| 𝐊∫_0^t e^ - η (t-s) Me^-as |z(0)| ds
       ≤     e^-η t |w_j(0)| + |b_j| 𝐊 e^-η t|z(0)| ∫_0^t e^(η-a) s  ds
       =     e^-η t |w_j(0)| + |b_j| 𝐊 e^-η t|z(0)| 1/η - a (e^(η-a)t - 1)
       =     e^-η t |w_j(0)| + b_X 𝐊 |z(0)| 1/η - a (e^-at - e^-η t).

Above computations have been performed for the case when η≠ a. If a= η, then
it holds that

    |w_j(t)|    ≤     e^-η t |w_j(0)| + b_X t e^-η t|z(0)|.


Now, using formula w(t)^2_X = ∑_j=1^∞ |w_j(t)|^2, one can show that w(t)_X decays exponentially as well.

Overall, if the controller (<ref>) stabilizes (<ref>) exponentially, then it stabilizes the system
(<ref>) as well.









§ OLD STUFF (PROBABLY FOR DELETION)


OLDER STUFF from Section <ref>

In particular for C=K it holds, for all j=1,…, n and for all x in ℝ^n

    ϕ_j(KX)_j(KX ) ≤ 0

and also 

    ϕ_j(KX) (KX)_j ≤ 0

just by noting that  _j(KX) and (KX)_j have the same sign. This latter inequality is a property that we needed in our computations.

Another fact that may be useful is the following. Under some conditions (that CP needs to check), there exists a positive definite matrix P such that the function V̂ (X)=X ^⊤ PX̂ satisfies, along the solutions of (<ref>) issuing from V(X)≤ 1,

    V̇_1≤ - α X ^2



Moreover there exists l≤ 1 such that for all along the solutions of (<ref>) issuing from V(X)≤ l, it holds additionally

    ( K X) =  KX

(that is the closed-loop system is linear).







Should be rechecked after we resolve above questions.



    [                                                                                                 d/dtṼ (w)                                                                                                         ≤           - α̃Ṽ(w) 
    + 2 ∑_j=1^∞w_jb_j(ϕ(K π_nw(t)))_j
    + 2  (π_nw) ^⊤ P (∑_j=1^n b_j(ϕ(K π_nw(t)))_j e_j );                                                                                                                                                                                                                   = - α̃Ṽ(w) + 2 ∑_j=1^∞w_jb_j( _ j (K π_nw) -(K π_nw)_j)
    + 2  (π_nw) ^⊤ P (b( (K π_nw) - K π_nw)) _j=1, ...n;                                                                                                                                                                                                                   =                                                    - α̃Ṽ(w) + 2 ∑_j=1^n w_jb_j( _ j (K π_nw) -(K π_nw)_j);                                                                                                                               + 2 ∑_j=n+1^∞w_jb_j _ j (K π_nw)
    + 2  (π_nw) ^⊤ P (b( (K π_nw) - K π_nw)) _j=1, ...n. ]





!!! Previous Parts !!!

Then it could be checked that along the dynamics (<ref>), it holds

    [                                                                                          d/dtṼ = 2 ∑_j=1^∞w_iAw_i + b u                                                                                                                       ≤                                         - αṼ + 2 ∑ _j w_j ( _ j (K X) -(K X)_j) + 2  X ^⊤ P  ( (K X) - K X ) _j=1, ...n;                                                                                                                                                                                                                                               ≤ - αṼ + 2 ∑ _j=1,…, n  w_j ( _ j (K X) -(K X)_j)+ 2 ∑ _j=1,…, n  w_j  _ j (K X)  + 2  X ^⊤ P  ( (K X) - K X ) _j=1, ...n ]

let us come back to this equation to finish the proof and to estimate the region of attraction for both definition of the saturation map. 










§ OLD CRAZY VERSION



Old version of the Section <ref>.

Now we come back to the problem (<ref>), which we rewrite using deadzone function ϕ.

    w_t=w_xx+cw+bu + bϕ(u),    ∀ t>0,  ∀ x ∈ (0,L), 
        w(t,0)=w(t,L)=0,   ∀ t>0, 
        w(0,x)=w^0(x) ,   ∀ x ∈ (0,L).

Here we continue to assume that b, c and u are scalar functions.

We would like to apply to this equation the finite-dimensional control, which has been designed in 
Section <ref>.
To this end we set:

    u(t):=𝐊 w(t) = K π_nw(t).

Since K ∈ L(X_n, X_n), we have the following decomposition:

    u(t)=𝐊 w(t) = ∑_j=1^n (K π_nw(t))_j e_j,

where 

    (K π_nw(t))_j := K π_nw(t)e_j.


Using this control, the spectral decomposition of the state and the fact that Ae_j=λ_j e_j, the first equation can be rewritten as 

    d/dt(∑_j=1^∞ w_j(t) e_j )    =    A (∑_j=1^∞ w_j(t) e_j ) + 
    ∑_j=1^n (bK π_nw(t))_j e_j + ∑_j=1^∞(bϕ(K π_nw(t)))_j e_j
       =   ∑_j=1^∞( λ_j w_j(t) + (bK π_nw(t))_j + (bϕ(K π_nw(t)))_j )e_j.

Note that the saturation function as well as the multiplication by b do not necessarily leave the set X_n invariant, therefore the term ϕ(u) has to be decomposed over the whole Hilbert basis of L_2(0,L).

Hence we have the following set of differential equations:

    d/dtw_j(t)= λ_j w_j(t) + (bK π_nw(t))_j + (bϕ(K π_nw(t)))_j,   j=1,…, +∞.



Now the problem is that even if there is no saturation (i.e. (u)= u for all u and hence ϕ≡ 0),
then the decomposition seemingly does not divide the problem in the bad finite-dimensional and good infinite-dimensional part. Instead, all equations look the same (unless b is a constant, in which case the decomposition indeed works).



















We define z(t):=(w_1(t),…,w_n(t))^⊤. 
With this notation one formally write π_nw= z^⊤(t) (e_1,…,e_n)^⊤ = ∑_j=1^n  w_j(t) e_j.

First let us consider the following Lyapunov function candidate:

    Ṽ(w(t)):= w(t)^2_X +  z^⊤(t) P z(t),    w(t) ∈ X,

where P is computed as in the previous section.


I have some problems with the understanding the term (π_nw) ^⊤ P (π_nw) 
(which we previously wrote as X ^⊤ P X).
We have that π_nw∈ X_n ⊂ X. Thus, π_nw is merely a function.
It seems that then P is just a constant.

Or am I mistaken?




Using the spectral decomposition of w∈ X, we can rewrite the definition of a Lyapunov function Ṽ as:

    Ṽ(w) =Ṽ(∑_j=1^∞  w_j e_j)= ∑_j=1^∞ w^2_j +  (π_nw) ^⊤ P (π_nw).

Let us compute the Lie derivative of Ṽ along the dynamics (<ref>).
We have:

    [                                                                                                d/dtṼ (w)                                                                                                        =                                        2 ∑_j=1^∞ w_j (λ_j w_j(t) +  (bK π_nw(t))_j + (bϕ(K π_nw(t)))_j );                                                                                                                                  + 2  (π_nw)^⊤ P (∑_j=1^n (λ_j w_j(t) + (bK π_nw(t))_j + (bϕ(K π_nw(t)))_j) e_j );                                                                                                                                                                                                                 = 2 ∑_j=1^∞  w_j (λ_j w_j(t) + (bK π_nw(t))_j)+ 2  (π_nw)^⊤ P (∑_j=1^n (λ_j w_j(t) + (bK π_nw(t))_j) e_j );                                                                                                                              + 2 ∑_j=1^∞  w_j (bϕ(K π_nw(t)))_j
    + 2  (π_nw)^⊤ P (∑_j=1^n (bϕ(K π_nw(t)))_j e_j ). ]

The terms in the first line of the previous expression are a derivative of Ṽ with respect to the systems
(<ref>), where the control is not saturated.
Since we have assumed that the controller u which we have constructed achieves exponential stabilization of this system,
there is some constant α̃ so that 



Should be rechecked after we resolve above questions.



    [                                                                                                 d/dtṼ (w)                                                                                                         ≤           - α̃Ṽ(w) 
    + 2 ∑_j=1^∞w_jb_j(ϕ(K π_nw(t)))_j
    + 2  (π_nw) ^⊤ P (∑_j=1^n b_j(ϕ(K π_nw(t)))_j e_j );                                                                                                                                                                                                                   = - α̃Ṽ(w) + 2 ∑_j=1^∞w_jb_j( _ j (K π_nw) -(K π_nw)_j)
    + 2  (π_nw) ^⊤ P (b( (K π_nw) - K π_nw)) _j=1, ...n;                                                                                                                                                                                                                   =                                                    - α̃Ṽ(w) + 2 ∑_j=1^n w_jb_j( _ j (K π_nw) -(K π_nw)_j);                                                                                                                               + 2 ∑_j=n+1^∞w_jb_j _ j (K π_nw)
    + 2  (π_nw) ^⊤ P (b( (K π_nw) - K π_nw)) _j=1, ...n. ]






!!! Previous Parts !!!

Then it could be checked that along the dynamics (<ref>), it holds

    [                                                                                          d/dtṼ = 2 ∑_j=1^∞w_iAw_i + b u                                                                                                                       ≤                                         - αṼ + 2 ∑ _j w_j ( _ j (K X) -(K X)_j) + 2  X ^⊤ P  ( (K X) - K X ) _j=1, ...n;                                                                                                                                                                                                                                               ≤ - αṼ + 2 ∑ _j=1,…, n  w_j ( _ j (K X) -(K X)_j)+ 2 ∑ _j=1,…, n  w_j  _ j (K X)  + 2  X ^⊤ P  ( (K X) - K X ) _j=1, ...n ]

let us come back to this equation to finish the proof and to estimate the region of attraction for both definition of the saturation map. 













———————————————-   

What follows are parts which may be useful in the future, but which are not connected to the previous notes.

———————————————-





§ LOCAL STABILIZATION OF LINEAR PARABOLIC PDES VIA CONTROLS WITH SATURATIONS




To stabilize the unstable finite-dimensional part we need to assume that the pair (𝐀,𝐁) is controllable which is equivalent to each b_j for  j=1, ..., n different to 0.










Now we come back to the problem (<ref>), which we rewrite using deadzone function ϕ.

    w_t(t,x)=w_xx(t,x)+c(x)w(t,x)+b(x) u(t) + b(x) ϕ(u(t)),    ∀ t>0,  ∀ x ∈ (0,L), 
        w(t,0)=w(t,L)=0,   ∀ t>0, 
        w(0,x)=w^0(x) ,   ∀ x ∈ (0,L).

Here we continue to assume that b, c and u are scalar functions.

We would like to apply to this equation the finite-dimensional control, which has been designed in 
Section <ref>.

As in Section <ref>, we consider the following scalar controller, see (<ref>):

    u(t):=K̂ w(t) = 𝐊π_nw(t) = 𝐊z(t),

where K ∈ L(X_n, ) and 𝐊∈^n.



Using this control, the spectral decompositions (<ref>), (<ref>)
and the fact that Ae_j=λ_j e_j, the first equation can be rewritten as 

    d/dt(∑_j=1^∞ w_j(t) e_j )    =    A (∑_j=1^∞ w_j(t) e_j ) + 
    ∑_j=1^n b_j 𝐊z(t) e_j + ∑_j=1^∞ b_j ϕ(𝐊z(t)) e_j
       =   ∑_j=1^∞( λ_j w_j(t) + b_j 𝐊z(t) + b_j ϕ(𝐊z(t)) )e_j.

Hence the system (<ref>) with a control (<ref>) is equivalent to the following set of differential equations:

    d/dtw_j(t)= λ_j w_j(t) + b_j 𝐊z(t) + b_j ϕ(𝐊z(t)),   j=1,…, +∞.

Using the notation (<ref>), we can rewrite the first n equations of this system as

    z'(t) = 𝐀z(t) + 𝐁𝐊z(t) + 𝐁ϕ(𝐊z(t))
    = 𝐀 z(t) + 𝐁(𝐊z(t)).


 
The system (<ref>) will be analyzed via two approaches. The first is via construction of the LF for the whole system, and the second one in Section <ref> 
is based on the decoupling of the system (<ref>), and seems to be easier in terms of computations and is better related to Section <ref>.





 §.§ Approach with a construction of LF for the whole system


First let us consider the following Lyapunov function candidate for (<ref>)

    Ṽ(w(t)):= w(t)^2_X +  z^⊤(t) P z(t),    w(t) ∈ X,

where P will be specified later.



Using the spectral decomposition of w∈ X, Ṽ can be written as:

    Ṽ(w) =Ṽ(∑_j=1^∞  w_j e_j)= ∑_j=1^∞ w^2_j +  z(t)^⊤ P z(t).

Let us compute the Lie derivative of Ṽ along the dynamics (<ref>).
We have:

    d/dtṼ (w(t))    =     2 ∑_j=1^∞ w_j(t) (λ_j w_j(t) + b_j 𝐊z(t) + b_j ϕ(𝐊z(t)) ) 
                     + 2  z(t)^⊤ P(𝐀z(t) + 𝐁𝐊z(t) + 𝐁ϕ(𝐊z(t))) 
       =     2 ∑_j=n+1^∞ w_j(t) (λ_j w_j(t) + b_j 𝐊z(t) + b_j ϕ(𝐊z(t)) ) 
          +2 ∑_j=1^n w_j(t) (λ_j w_j(t) + b_j 𝐊z(t) + b_j ϕ(𝐊z(t)) ) 
                     + 2  z(t)^⊤ P(𝐀z(t) + 𝐁𝐊z(t) + 𝐁ϕ(𝐊z(t)))

Now we would like to transform the terms in (<ref>) the second and third line of the previous expression:

    2 ∑_j=1^n w_j(t)   (λ_j w_j(t) + b_j 𝐊z(t) + b_j ϕ(𝐊z(t)) ) + 2  z(t)^⊤ P(𝐀z(t) + 𝐁𝐊z(t) + 𝐁ϕ(𝐊z(t))) 
     
    =    2 z^⊤(t) 𝐀 z(t) + 2 z^⊤(t) 𝐁𝐊z(t) + 2 z^⊤(t) 𝐁ϕ(𝐊z(t)) 
                  + 2  z(t)^⊤ P(𝐀z(t) + 𝐁𝐊z(t) + 𝐁ϕ(𝐊z(t))) 
     
    =     2 z^⊤ (t)(𝐀 + 𝐁𝐊)z(t) + 2 z^⊤ (t)P(𝐀 + 𝐁𝐊) z(t) + 2 z^⊤ (t)(P+I)𝐁ϕ(𝐊z(t)).


Now assume that the matrices C,D are chosen as in Lemma <ref>.
Then we continue estimates (<ref>) as follows:

    2z^⊤ (t)   (𝐀 + 𝐁𝐊)z(t) + 2 z^⊤ (t)P(𝐀 + 𝐁𝐊) z(t) + 2 z^⊤(t)(P+I)𝐁ϕ(𝐊z(t)) 
    ≤   
    2z^⊤ (t)(𝐀 + 𝐁𝐊)z(t) + 2 z^⊤ (t)P(𝐀 + 𝐁𝐊) z(t) + 2 z^⊤ (t)(P+I)𝐁ϕ(𝐊z(t))
    - 2 ϕ( 𝐊z)^⊤𝐃 (ϕ(𝐊z)+𝐂 z)
    
    =    
    z^⊤ (t)(𝐀 + 𝐁𝐊)z(t)
    +
    [[     z; ϕ(𝐊z) ]]^⊤[[ (𝐀+ 𝐁𝐊)^⊤ P + P (𝐀+ 𝐁𝐊)      (P + I) 𝐁 - (𝐃C)^⊤;      ( (P + I)𝐁)^⊤ - 𝐃C                    -2 𝐃 ]]
    [[     z; ϕ(𝐊z) ]]
    
    =   [[     z; ϕ(𝐊z) ]]^⊤[[ (𝐀+ 𝐁𝐊)^⊤ (P+I) + (P+I) (𝐀+ 𝐁𝐊)              (P + I) 𝐁 - (𝐃C)^⊤;              ( (P + I)𝐁)^⊤ - 𝐃C                            -2 𝐃 ]]
    [[     z; ϕ(𝐊z) ]]


If we assume that C, D, P and 𝐊 can be chosen so that 
the matrix

    [[ (𝐀+ 𝐁𝐊)^⊤ (P+I) + (P+I) (𝐀+ 𝐁𝐊)              (P + I) 𝐁 - (𝐃C)^⊤;              ( (P + I)𝐁)^⊤ - 𝐃C                            -2 𝐃 ]]

is negative definite, then we have to proceed as in Proposition <ref>
to obtain some estimates for the stability region.

Clearly, we have also to estimate somehow the term
2 ∑_j=n+1^∞ w_j(t) (λ_j w_j(t) + b_j 𝐊z(t) +
b_j ϕ(𝐊z(t)) )
in (<ref>).






§ AN APPROACH WITHOUT FULL DECOMPOSITION


Also we can take the notation x=(x_1,…,x_n) = ∑_i=1^n x_i e_i.


Alternatively, we can compute the Lie derivative of Ṽ without using decomposition of the state space.
According to (<ref>) we have:

    d/dtṼ (w)    =    2 ∫_0^L w(z)ẇ(z) dz +  d/dt ( X ^⊤ P X)
       =    2 ∫_0^L w(z)(w_xx(z)+c(z)w+b(z)KX + b(z)ϕ(KX)) dz +  d/dt ( X ^⊤ P X)
       =    2 ∫_0^L w(z)(w_xx(z)+c(z)w(z)+b(z)KX)dz +  d/dt ( X ^⊤ P X)
           + 2 ∫_0^L w(z) b(z)ϕ(KX) dz.










§ A SECOND LYAPUNOV FUNCTION FOR THE NONLINEAR PDE



Now let us consider the following Lyapunov function candidate:

    V̂ (∑_j  w_j e_j)= ∑_j> n -1/λ _j w_j ^2+  X ^⊤ P X


And we have, along the solutions to (<ref>),

    [                                                                                   d/dtV̂                                                                                        ≤ - αV̂ + 2 ∑ _j -1/λ _j w_j ( _ j (K X) -(K X)_j) + 2  X ^⊤ P  ( (K X) - K X ) _j=1, ...n ]









