Local stabilization of an unstable parabolic equation via saturated controls

TL;DR

This paper presents a novel saturated feedback control method that stabilizes unstable linear reaction-diffusion equations without requiring the uncontrolled system to be stable, using Lyapunov techniques and matrix inequalities.

Contribution

It introduces a new control approach for unstable systems, extending stabilization techniques to more general cases with various control types and saturations.

Findings

Effective stabilization of unstable systems demonstrated

Region of attraction estimated via matrix inequalities

Numerical simulations confirm control performance

Abstract

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.

Equations

$$\begin{array}{rl}& w_t(t,x)=w_{xx}(t,x)+c(x)w(t,x)\\ & +\sum _{k=1}^m{b}_k(x)\text{sat}(u_k(t)),t>0,x\in (0,L),\\ & w(t,0)=w(t,L)=0,t>0,\\ & w(0,x)=w^0(x),x\in (0,L).\end{array}$$

$$\text{sat}(v{)}_k:=\{\begin{array}{ll}{v}_k& \text{ if }|v_k|⩽\ell ,\\ \frac{\ell }{|v_k|}{v}_k& \text{ if }|v_k|⩾\ell ,\end{array}$$

$$A={\partial }_{xx}+c(\cdot )\mathrm{id}:X\to X$$

$$w_t(t,\cdot )=Aw(t,\cdot )+\sum _{k=1}^m{b}_k\text{sat}(u_k(t)).$$

$$-\infty <\cdots <{\lambda }_j<\cdots <{\lambda }_1\text{and}{\lambda }_j\underset{j\to +\infty }{⟶}-\infty$$

$$w(t,\cdot )$$

$$w_j(t)$$

$$b_k(\cdot )=\sum _{j=1}^{\infty }{b}_{jk}{e}_j(\cdot ),b_{jk}=⟨b_k(\cdot ),e_j(\cdot ){⟩}_{L_2(0,L)},j\in {\mathrm{I}\mathrm{N}}^{\ast }.$$

$${\dot{w}}_j(t)$$

$$={\lambda }_j{w}_j(t)+{\mathbf{b}}_j\cdot \text{sat}(u(t)),j\in {\mathrm{I}\mathrm{N}}^{\ast },$$

$$\forall j>n:{\lambda }_j<-\eta <0.$$

$$z:=\left(\begin{array}{c}{w}_1\\ ⋮\\ w_n\end{array}\right),\mathbf{A}:=\left(\begin{array}{ccc}{\lambda }_1& \cdots & 0\\ ⋮& \ddots & ⋮\\ 0& \cdots & {\lambda }_n\end{array}\right),\mathbf{B}:=\left(\begin{array}{ccc}{b}_{11}& \cdots & b_{1m}\\ ⋮& & ⋮\\ b_{n1}& \cdots & b_{nm}\end{array}\right)$$

$$\dot{z}(t)=\mathbf{A}z(t)+\mathbf{B}\text{sat}(u(t)).$$

$${\pi }_{n}w(t,\cdot ):=\sum _{j=1}^n{w}_j(t)e_j(\cdot ).$$

$$\ell_{2}(\mathrm{I\kern-2.10007ptN}^{*},\mathrm{I\kern-2.10007ptR}):=\Big{\{}(x_{k})_{k\in\mathrm{I\kern-1.47005ptN}^{*}}\in\mathrm{I\kern-2.10007ptR}^{\mathrm{I\kern-1.47005ptN}^{*}}\ :\ \sum_{k=1}^{\infty}|x_{k}|^{2}<\infty\Big{\}},$$

$$u=$$

$$=$$

$${\dot{w}}_j(t)={\lambda }_j{w}_j(t)+{\mathbf{b}}_j\cdot \text{sat}(\mathbf{K}z(t)),j\in {\mathrm{I}\mathrm{N}}^{\ast }.$$

$$\dot{z}(t)=\mathbf{A}z(t)+\mathbf{B}\text{sat}(\mathbf{K}z(t)).$$

$${\dot{w}}_j(t)$$

$$\dot{z}=\mathbf{A}z+\mathbf{B}u$$

$$|z(t;z_0)|⩽M{e}^{-at}|z_0|\forall t⩾0$$

$$\Vert w(t,\cdot ){\Vert }_X⩽M{e}^{-at}\Vert w(0,\cdot ){\Vert }_X\forall t⩾0.$$

$$|z(t;z_0)|⩽M{e}^{-at}|z_0|,t⩾0.$$

$$w_j(t)=e^{{\lambda }_{j}t}{w}_j(0)+{\mathbf{b}}_j\cdot {\int }_0^t{e}^{{\lambda }_j(t-s)}\text{sat}(\mathbf{K}z(s))ds.$$

$$|\text{sat}(\mathbf{K}z)|⩽|\mathbf{K}z|⩽\Vert \mathbf{K}\Vert |z|.$$

$$\displaystyle|b_{jk}|=\Big{|}\langle b_{k}(\cdot),e_{j}(\cdot)\rangle_{X}\Big{|}\leqslant\|b_{k}\|_{X}\|e_{j}\|_{X}=\|b_{k}\|_{X}.$$

$$|w_j(t)|⩽$$

$$⩽$$

$$=$$

$$|w_j(t)|$$

$$\begin{split}&w_{t}(t,x)=w_{xx}(t,x)+c(x)w(t,x)\\ &+\sum_{k=1}^{m}b_{k}(x)\textup{sat}\big{(}(\mathbf{K}z(t))_{k}+d_{k}(t)\big{)},\ t>0,\;x\in(0,L),\\ &w(t,0)=w(t,L)=0,\quad t>0,\\ &w(0,x)=w^{0}(x),\quad x\in(0,L).\end{split}$$

$$\Vert w(t,\cdot ){\Vert }_X⩽M_2{e}^{-a_{2}t}\Vert w(0,\cdot ){\Vert }_X+{\gamma }_2\underset{s⩾0}{\sup }|d(s)|.$$

Full text

I Introduction

In applications of control technology, physical inputs (like force, torque, thrust, stroke, etc.) are often limited in size [3]. 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. [41, 45, 33] 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 [19] a backstepping approach is applied to design a boundary delayed feedback control for a heat equation (see [20] for further results using the same design methods). In [10] a stable heat partial differential equation (PDE) is controlled by means of a delayed bounded linear control operator (see also [37] for the semilinear case). For the case of unbounded control operators, see [29, 31] 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 [36], where compact and bounded control operators are considered. Another notable early reference is [21] 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 [25] for a discussion.

The definition of the saturation map used in the present paper is aligned to that used in [30] where Lyapunov methods are shown to be useful for the stability analysis of wave equations subject to saturated inputs. See also [25, 23], 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 [36, 21, 30], 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 [42]).*

We would like to mention the related research on model predictive control of parabolic PDEs under constraints on the state and the input [6], and on Lyapunov-based control of parabolic systems by controls of a bounded magnitude [7]. 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 [16] 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 [17] boundary stabilization of a nonlinear Schrödinger equation with state delay and bounded internal disturbance is performed. Both in [16, 17] 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 [32, 2, 1]. To stabilize the unstable finite-dimensional part using a saturated control we apply techniques, which are well-known for ODE systems (see e.g. [41, 45, 33] for systems with non-delayed input, and [22] 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 [4, 35, 43] 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 [39]. 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 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_{\infty}$) 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 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. [9]), 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 $\mathrm{I\kern-2.10007ptR}_{+}$. The Euclidean norm on $\mathrm{I\kern-2.10007ptR}^{n}$ is denoted by $|\cdot|$, the operator norm induced by this norm on spaces of matrices is denoted by $\|\cdot\|$. The interior of a set $S$ in a topological space is denoted ${\rm int}\,S$ and $\overline{S}$ denotes its closure. In any normed vector space, the ball of radius $r$ around [math] is denoted by $B_{r}(0)$. By $\mathrm{I\kern-2.10007ptN}$ we denote the nonnegative integers. For convenience, $\mathrm{I\kern-2.10007ptN}^{*}:=\mathrm{I\kern-2.10007ptN}\setminus\{0\}$. For $k\in\mathrm{I\kern-2.10007ptN},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 $\leqslant 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)$.

II Problem formulation

We consider the stabilization problem of a one-dimensional linear reaction-diffusion equation by means of a distributed control $u:\mathrm{I\kern-2.10007ptR}_{+}\to\mathrm{I\kern-2.10007ptR}^{m}$. Let $L>0$. We are given $m$ functions $b_{k}:[0,L]\to\mathrm{I\kern-2.10007ptR}$, $k=1,\ldots,m$ which describe at which places the control input $u_{k}\in\mathrm{I\kern-2.10007ptR}$ is acting. The function $c$ models the place-dependent reaction rate. The system model is then

[TABLE]

We assume that the state space of this system is $X:=L_{2}(0,L)$ and that $c,b_{k}\in X$, $k=1,\ldots,m$.

Here sat is a component-wise saturation function, that is, for all $k=1,\ldots,m$ and for any $v\in\mathrm{I\kern-2.10007ptR}^{m}$,

[TABLE]

where $\ell>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. [31] as well as Section V-B below.   $\circ$

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

[TABLE]

with domain $D(A)=H^{2}(0,L)\cap H^{1}_{0}(0,L)$. Then the control system (1) takes the form

[TABLE]

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 [18, Theorem III.6.29]. Furthermore, there exists a Hilbert basis $(e_{j})_{j\geqslant 1}$ of $X$ consisting of eigenfunctions of $A$, associated with the sequence of eigenvalues $(\lambda_{j})_{j\geqslant 1}$. Note that

[TABLE]

and that $e_{j}(\cdot)\in D(A)$ for every $j\geqslant 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}\in L_{1,loc}([0,\infty))$, for $k=1,\ldots,m$.

Every solution $w(t,\cdot)\in D(A)$ of (4) can be expanded as a series in the eigenfunctions $e_{j}(\cdot)$, convergent in $H_{0}^{1}(0,L)$,

[TABLE]

Analogously, we can expand the coefficients $b_{k}$ in the series

[TABLE]

As discussed in Appendix -A, (4) is equivalent to the infinite-dimensional control system

[TABLE]

where “$\cdot$” is the scalar product in $\mathrm{I\kern-2.10007ptR}^{m}$, $\textup{sat}(u(t))\in\mathrm{I\kern-2.10007ptR}^{m}$ is the vector with entries $\textup{sat}(u_{k}(t))$ and $\mathbf{b}_{j}$ is the row vector with entries $b_{jk}$, $k=1,\ldots,m$.

Let $n\in\mathrm{I\kern-2.10007ptN}^{*}$ be the number of nonnegative eigenvalues of $A$ and let $\eta>0$ be such that

[TABLE]

With the matrix notations

[TABLE]

the $n$ first equations of (II) form the unstable finite-dimensional control system

[TABLE]

Remark 2

(Possible extensions of our results).

(i) The same method can be used for stabilization with a prescribed decay rate $-\beta$. Just define for a suitable $\eta>\beta$ 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 [14, 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 [9]), the methods of this paper apply.   $\circ$

III Estimation of the region of attraction for saturated inputs

III-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}(\cdot))_{i=1}^{n}$ and $\pi_{n}$ be the orthogonal projection onto $X_{n}$, that is

[TABLE]

We define also $X_{n}^{\bot}$ as the orthogonal complement of $X_{n}$ in $X$. Let $\iota:\mathrm{I\kern-2.10007ptR}^{n}\rightarrow X_{n}$ be the isomorphism defined by $\iota(e^{j})=e_{j}(\cdot)$, where $(e^{j})_{j=1,\ldots,n}$ is the canonical basis of $\mathrm{I\kern-2.10007ptR}^{n}$. We will use the isometric representation of $L_{2}(0,L)$ as $\ell_{2}(\mathrm{I\kern-2.10007ptN}^{*},\mathrm{I\kern-2.10007ptR})$ obtained by the isomorphism induced by $e_{j}(\cdot)\mapsto e^{j}$, where

[TABLE]

where $e^{j},j\in\mathrm{I\kern-2.10007ptN}^{*}$ are the standard basis vectors in $\ell_{2}(\mathrm{I\kern-2.10007ptN}^{*},\mathrm{I\kern-2.10007ptR})$ and we use the standard norm on that space. Corresponding to the decomposition $L_{2}(0,L)=X_{n}\bigoplus X_{n}^{\bot}$, where “$\bigoplus$” is the orthogonal sum of subspaces, we denote $\ell_{2}(\mathrm{I\kern-2.10007ptN}^{*},\mathrm{I\kern-2.10007ptR})=\mathrm{I\kern-2.10007ptR}^{n}\bigoplus\ell_{2,j>n}$, where we identify $\mathrm{I\kern-2.10007ptR}^{n}$ with the sequences with support in $\{1,\ldots,n\}$ and $\ell_{2,j>n}$ is the set of sequences in $\ell_{2}(\mathrm{I\kern-2.10007ptN}^{*},\mathrm{I\kern-2.10007ptR})$ which are [math] in the first $n$ entries.

Given a linear map $K:X_{n}\to\mathrm{I\kern-2.10007ptR}^{m}$, consider the feedback

[TABLE]

where $\mathbf{K}_{j}:=Ke_{j}(\cdot)\in\mathrm{I\kern-2.10007ptR}^{m}$, $j=1,\ldots,n$, and where we use the notation from (8) in the final step and set

$\mathbf{K}:=(\mathbf{K}_{1},\ldots,\mathbf{K}_{n})\in\mathrm{I\kern-2.10007ptR}^{m\times n}$.

Hence the system (II) with the feedback (11) is equivalent to the following set of differential equations:

[TABLE]

Using (8), we rewrite the first $n$ equations of (12) as

[TABLE]

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

[TABLE]

Our general assumption is:

Assumption 1

The pair $(\mathbf{A},\mathbf{B})$ is stabilizable, i.e. there is $\mathbf{K}\in\mathrm{I\kern-2.10007ptR}^{m\times n}$: $\mathbf{A}+\mathbf{B}\mathbf{K}$ is Hurwitz.

Remark 3

Clearly, the closed-loop system (13) with a matrix $\mathbf{K}$ as in Assumption 1 is locally asymptotically stable.

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

[TABLE]

is locally and hence globally asymptotically stabilized by means of the feedback $u(t):=\mathbf{K}z(t)$. Thus, local asymptotic stability of (13) implies that the pair $(\mathbf{A},\mathbf{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}\neq 0$ for all $j=1,\ldots,n$ and $\lambda_{k}\neq\lambda_{j}$ for all $k,j=1,\ldots,n$, $k\neq 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.   $\circ$

Next we show that the problem of exponential stabilization of the overall system (II) boils down to the exponential stabilization of the finite-dimensional unstable system (9). This latter problem will be elaborated in Section III-B.

Definition 1

Assume that $\mathbf{K}$ is chosen so that [math] is a locally asymptotically stable fixed point of (13). We say that $S$ is a region of attraction of [math] if

• (i)

$0\in{\rm int}\,S$;

• (ii)

for any $z_{0}\in S$ the corresponding solution of (13) satisfies $z(t;z_{0})\to 0$ as $t\to\infty$;

• (iii)

$S$* is forward invariant, i.e. for any $z_{0}\in S$ it holds that $z(t;z_{0})\in S$ for all $t\geqslant 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 [12].

Definition 2

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

• (i)

there exist $\varepsilon,M,a>0$ such that for any initial condition satisfying $|z_{0}|\leqslant\varepsilon$, it holds

[TABLE]

• (ii)

$\overline{B_{\varepsilon}(0)}\subset 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}\in\mathrm{I\kern-2.10007ptR}^{n}$.

Definitions 1 and 2 can be stated analogously for system (12).

We note that if the maximal region of attraction is not $\mathrm{I\kern-2.10007ptR}^{n}$, then the system cannot be exponentially stable on the maximal region of attraction (for Lipschitz continuous systems), see [12]. Thus by analyzing regions of attraction with exponential stability we necessarily restrict the region of attraction.

Proposition 1

Assume $\mathbf{K}$ is chosen such that the subsystem (13) is locally exponentially stable in [math] with region of attraction $S\subset\mathrm{I\kern-2.10007ptR}^{n}$. Then:

• (i)

system (12) is locally exponentially stable in [math] with region of attraction $S\times\ell_{2,j>n}$.

• (ii)

system (1) with the feedback (11) is locally exponentially stable in [math] with region of attraction $\iota(S)\times X_{n}^{\bot}$.

In addition, for any closed and bounded set $G\subset{\rm int}\,(\iota(S)\times X_{n}^{\bot})$, there exist two positive values $M$ and $a$ such that for any initial condition $w(0,\cdot)$ in $G$, the solution $w(\cdot)$ to (1) with the controller (11) satisfies

[TABLE]

Proof:

Pick a compact subset $G^{\prime}$ of ${\rm int}\,S$. Since we assume that (13) is locally exponentially stable in 0 with region of attraction $S\subset\mathrm{I\kern-2.10007ptR}^{n}$, it follows from a standard compactness argument that there exist $M,a>0$ so that for any $z_{0}\in G^{\prime}$ the solution $z(\cdot;z_{0})$ to (13) satisfies

[TABLE]

From equations (II) and (11), we derive that for $j=n+1,\ldots,\infty$, for any $t\geqslant 0$ and for any $(w_{n+1}(0),w_{n+2}(0),\ldots)\in\ell_{2,j>n}$ it holds that

[TABLE]

From (2) it follows that for all $z\in\mathrm{I\kern-2.10007ptR}^{n}$ we have

[TABLE]

Also due to the Cauchy-Bunyakovsky-Schwarz inequality we have that, for all $j\in\mathrm{I\kern-2.10007ptN}$,

[TABLE]

Thus, for all $j\geqslant n+1$, we obtain exploiting (18) that

[TABLE]

The above computations have been performed for the case when $\eta\neq a$. If $a=\eta$, then it holds that

[TABLE]

Now, using the inequality $(a+b)^{2}\leqslant 2(a^{2}+b^{2})$ for any $(a,b)\in\mathrm{I\kern-2.10007ptR}^{2}$, and the square summability of $|w_{j}(0)|$ and $|b_{jk}|$, $k=1,\ldots,m$, it follows that $\sum_{j=n+1}^{\infty}|w_{j}(t)|^{2}$ decays exponentially as well.

We now obtain local exponential stability of (12) by choosing $G^{\prime}$ such that $0\in\mathrm{I\kern-2.10007ptR}^{n}$ is in the interior of $G^{\prime}$ and noting that then $0\in X$ is in the interior of $\iota(G^{\prime})\times X_{n}^{\bot}$.

For the final statement of the proposition, pick a closed and bounded set $G\subset{\rm int}\,(\iota(S)\times X_{n}^{\bot})$. Select $G^{\prime}=\iota^{-1}\circ\pi_{n}(G)$, then $G^{\prime}$ is a compact subset of ${\rm int}\,S$, the previous computations yield (17) for suitable constants $M$ and $a$ and for the superset $\iota(G^{\prime})\times X_{n}^{\bot}$ which contains $G$. ∎

Remark 4

It is not hard to see that (14) is input-to-state 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.   $\circ$

Remark 5

Our feedback controller is robust w.r.t. additive actuator disturbances. Let the control input to system (1) be $u(t):=\mathbf{K}z(t)+d(t)$, where $d\in PC(\mathrm{I\kern-2.10007ptR}_{+},\mathrm{I\kern-2.10007ptR}^{m})$ is a piecewise continuous actuator disturbance. Then there are $r>0$, $M_{1},a_{1},\gamma_{1}>0$ so that for all $w(0,\cdot)\in B_{r}$ and all $d\in PC(\mathrm{I\kern-2.10007ptR}_{+},\mathrm{I\kern-2.10007ptR}^{m})$ with $\sup_{s\geqslant 0}|d(s)|<r$ the solution of

[TABLE]

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

[TABLE]

This property can be obtained quite easily, since for small enough $w$ and $d$ we have $\textup{sat}(\mathbf{K}z(t)+d(t))=\mathbf{K}z(t)+d(t)$, and thus (19) is a linear system with a bounded disturbance operator $d\mapsto\sum_{k=1}^{m}b_{k}(\cdot)d_{k}$, acting on $PC(\mathrm{I\kern-2.10007ptR}_{+},\mathrm{I\kern-2.10007ptR}^{m})$, and (20) follows as exponential stability of (19) without disturbances implies LISS for a system (19) with disturbances. It is, however, much harder to estimate a region of attraction of system (19), 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 [26, 28, 40] and references therein.   $\circ$

III-B Estimate of the region of attraction for the finite-dimensional 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 (13). We perform this task in this section.

Defining the deadzone nonlinearity $\phi:\mathrm{I\kern-2.10007ptR}^{m}\to\mathrm{I\kern-2.10007ptR}^{m}$ by

[TABLE]

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 $\mathbf{C}\in\mathbb{R}^{m\times n}$, $z\in\mathbb{R}^{n}$ and $j\in\{1,\ldots,m\}$ it holds that $|((\mathbf{K}-\mathbf{C})z)_{j}|\leqslant\ell$, then

[TABLE]

Remark 6

Note that other properties exist for saturation yielding to other representations of the saturation maps (see, e.g., [41, Section 1.7]). In the present work, we insist on using generalized sector conditions because they provide constructive methods to design Lyapunov functions.   $\circ$

As a consequence of Lemma 1, for all diagonal positive definite matrices $\mathbf{D}\in\mathbb{R}^{m\times m}$, for all $\mathbf{C}\in\mathbb{R}^{m\times n}$, and all $z\in\mathbb{R}^{n}$ such that $|((\mathbf{K}-\mathbf{C})z)_{j}|\leqslant\ell$, $j=1\ldots,m$, we have:

[TABLE]

We recall the following well-known Schur complement lemma (see e.g. [4, Chapter 2, p. 7]):

Lemma 2

*Let $A\in\mathrm{I\kern-2.10007ptR}^{n\times n}$, $B\in\mathrm{I\kern-2.10007ptR}^{m\times n}$, $C\in\mathrm{I\kern-2.10007ptR}^{m\times m}$ and let $M:=\begin{pmatrix}A&B^{\top}\\ B&C\end{pmatrix}$. If $C$ is positive definite, then $M$ is positive semidefinite if and only if its Schur complement $M/C:=A-B^{\top}C^{-1}B$ is positive semidefinite. *

The following result is only a small variation of [11, 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 $\mathbf{K}\in\mathrm{I\kern-2.10007ptR}^{m\times n}$ be such that $\mathbf{A}+\mathbf{B}\mathbf{K}$ is Hurwitz. Let $P\in\mathbb{R}^{n\times n}$ be symmetric positive definite, $\mathbf{D}\in\mathbb{R}^{m\times m}$ diagonal positive definite and $\mathbf{C}\in\mathbb{R}^{m\times n}$ such that

[TABLE]

and

[TABLE]

Then the finite-dimensional system (13) is locally exponentially stable in [math] with a region of attraction given by

[TABLE]

Moreover, in ${\cal A}$, the function $V_{1}$ defined by $V_{1}(z):=z^{\top}Pz$, $z\in\mathrm{I\kern-2.10007ptR}^{n}$, decreases exponentially fast to [math] along the solutions to (13), i.e. there is a constant $\alpha>0$ so that

[TABLE]

Proof:

Let us show that (22), (25) are feasible. As $\mathbf{A}+\mathbf{B}\mathbf{K}$ is Hurwitz, there is a symmetric positive definite matrix $P\in\mathbb{R}^{n\times n}$ such that $(\mathbf{A}+\mathbf{B}\mathbf{K})^{\top}P+P(\mathbf{A}+\mathbf{B}\mathbf{K})<0$ and, multiplying if needed $P$ by a positive constant, we may additionally assume that \left[\begin{array}[]{cc}P&\mathbf{K}^{\top}\\ \mathbf{K}&\ell^{2}I_{m}\end{array}\right]\geqslant 0. Then, letting $\mathbf{C}=0\in\mathbb{R}^{m\times n}$, and picking a diagonal positive definite matrix $\mathbf{D}\in\mathbb{R}^{m\times m}$ with a norm sufficiently large, we get diagonal positive definite matrix $\mathbf{D}\in\mathbb{R}^{m\times m}$ and a matrix $\mathbf{C}\in\mathbb{R}^{m\times n}$ such that (22) and (25) hold.

Now pick $P,\mathbf{D},\mathbf{C}$ such that 22) and (25) hold and consider the Lyapunov function candidate

[TABLE]

where $P\in\mathrm{I\kern-2.10007ptR}^{n\times n}$ is the symmetric positive definite matrix given by the assumptions. The time-derivative of $V_{1}$ along the solutions to (13) for $z\in\mathrm{I\kern-2.10007ptR}^{n}$ is given by

[TABLE]

Now assuming $|((\mathbf{K}-\mathbf{C})z)_{j}|\leqslant\ell$ for all $j=1,\ldots,m$, it follows from (21) that

[TABLE]

In view of (22), this implies the estimate (27) selecting $-\alpha$ as the maximal eigenvalue of $M_{1}$.

It remains to ensure that $|((\mathbf{K}-\mathbf{C})z)_{j}|\leqslant\ell$ is satisfied for all $j=1,\ldots,m$. To do this, we use the Lyapunov function $V_{1}$ and we impose that the ellipsoid $\{z\in\mathrm{I\kern-2.10007ptR}^{n}:\;z^{\top}Pz\leqslant 1\}$ is included in the ellipsoid $\{z\in\mathrm{I\kern-2.10007ptR}^{n}:\;|(\mathbf{K}-\mathbf{C})z|\leqslant\ell\}$ (this implies that $|((\mathbf{K}-\mathbf{C})z)_{j}|\leqslant\ell$ holds for all $j=1,\ldots,m$). This inclusion is equivalent to the inclusion of $\{z\in\mathrm{I\kern-2.10007ptR}^{n}:\;z^{\top}Pz\leqslant 1\}$ in the set $\{z\in\mathrm{I\kern-2.10007ptR}^{n}:z^{\top}(\mathbf{K}-\mathbf{C})^{\top}(\mathbf{K}-\mathbf{C})z\leqslant\ell^{2}\}$ which is again equivalent to the matrix inequality $P-(\mathbf{K}-\mathbf{C})^{\top}\frac{1}{\ell^{2}}(\mathbf{K}-\mathbf{C})\geqslant 0$. As $\ell^{2}I_{m}$ is positive definite, Lemma 2 ensures that this latter matrix inequality is equivalent to (25). ∎

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^{\top}Pz\leqslant 1\}$ can never be optimal. Indeed, in this region we have $\dot{V}_{1}(z)\leqslant-\alpha|z|^{2}$, so that by a continuity and compactness argument it follows that $\dot{V}_{1}(z)<0$ on an enlarged region of the form $\{z:\;z^{\top}Pz\leqslant 1+\varepsilon\}$, for a suitable $\varepsilon>0$.   $\circ$

With the previous result, it is also possible to analyze global stability. The region of attraction is global as soon as for all $z\in\mathbb{R}^{n}$, it holds that $|((\mathbf{K}-\mathbf{C})z)_{j}|\leqslant\ell$, $j=1,\ldots,m$. This is equivalent to $\mathbf{K}=\mathbf{C}$. We thus obtain the following corollary:

Corollary 1

If there exist a symmetric positive definite matrix $P\in\mathbb{R}^{n\times n}$ and a diagonal positive matrix $\mathbf{D}\in\mathbb{R}^{m\times m}$ such that (22) holds with $\mathbf{C}:=\mathbf{K}$, then the finite-dimensional system (13) is globally exponentially stable in [math]. Moreover the Lyapunov function $V_{1}$ decreases exponentially fast to [math] along the solutions to (13).

In our context, the main interest of Proposition 2 lies in the following consequence for system (1).

Theorem 1

Consider system (1) with $(\mathbf{A},\mathbf{B})$ stabilizable. Let $\mathbf{K}$ be such that $\mathbf{A}+\mathbf{B}\mathbf{K}$ is Hurwitz and $K$ be as in (11). Then the closed-loop system with the feedback $K$

[TABLE]

is locally exponentially stable in [math] with region of attraction $\imath({\cal A})\times X_{n}^{\perp}$. In addition, the constants of decay can be chosen uniformly on $\imath({\cal A})\times X_{n}^{\perp}$.

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 ${\cal A}$, the proof of Proposition 1 shows that the uniformity holds on all of $\imath({\cal A})\times X_{n}^{\perp}$. ∎

III-C Numerical implementation of Proposition 2

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

Let $S=P^{-1}$ and $\mathbf{E}=\mathbf{D}^{-1}$. Since $M_{1}<0$ if and only if

[TABLE]

we obtain that (22) is equivalent to

[TABLE]

Now pre- and post-multiplying (25) by \left[\begin{array}[]{cc}S&0\\ 0&I_{m}\end{array}\right], we get that (25) is equivalent to

[TABLE]

Letting $\mathbf{Y}=S\mathbf{C}^{\top}$, we obtain (see also [11, Theorem 1] in a different context):

Proposition 3

Existence of a symmetric positive definite matrix $P\in\mathbb{R}^{n\times n}$, a diagonal positive definite matrix $\mathbf{D}\in\mathbb{R}^{m\times m}$ and a matrix $\mathbf{C}\in\mathbb{R}^{m\times n}$ such that (22) and (25) hold is equivalent to the existence of a symmetric positive definite matrix $S\in\mathbb{R}^{n\times n}$, a diagonal positive definite matrix $\mathbf{E}\in\mathbb{R}^{m\times m}$ and $\mathbf{Y}\in\mathbb{R}^{m\times n}$ for which it holds that

[TABLE]

and

[TABLE]

This reformulation is important as the constraints (34) and (37) are linear in the variables $S$, $\mathbf{E}$ and $\mathbf{Y}$, in contrast to bilinear inequalities (22) and (25).

Finally, no structure is imposed on both the matrix $\mathbf{C}\in\mathbb{R}^{m\times n}$ and on the matrix $\mathbf{Y}\in\mathbb{R}^{m\times n}$, and we let $\mathbf{C}=(S^{-1}\mathbf{Y})^{\top}$ from the knowledge of $S$ and $\mathbf{Y}$ (and respectively we may let $\mathbf{Y}=P^{-1}\mathbf{C}^{\top}$ from the knowledge of $P$ and $\mathbf{C}$).

III-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 3

Consider a symmetric matrix $P=\begin{pmatrix}P_{11}&P_{12}\\ P_{21}&P_{22}\end{pmatrix}\in\mathrm{I\kern-2.10007ptR}^{(n_{1}+n_{2})\times(n_{1}+n_{2})}$, where $P_{ij}$ are matrices of appropriate dimension and $P_{11}>0$, $P_{22}>0$. Define $P_{a}:=\begin{pmatrix}aP_{11}&P_{12}\\ P_{21}&P_{22}\end{pmatrix}$ for $a\geqslant 0$. Then there exists $a^{*}\geqslant 0$:

• (i)

$P_{a^{*}}\geqslant 0$, and $P_{a^{*}}$ is not positive definite,

• (ii)

$P_{a}>0$* for all $a>a^{*}$,*

• (iii)

for all $a\in[0,a^{*})$ the matrix $P_{a}$ is not positive semidefinite.

Proof:

As $P_{22}>0$, the matrix $P_{a}$ is positive definite if and only if the Schur complement

[TABLE]

is positive definite, see Lemma 2. Note that as $(P_{12})^{\top}=P_{21}$ and $(P_{22}^{-1})^{\top}=P_{22}^{-1}$, we have that

[TABLE]

and as $P_{11}$ is positive definite, $P_{a}/P_{22}$ is a symmetric matrix.

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 $\lambda_{\mathrm{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^{*}\geqslant 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 $\mathbf{K}\in\mathrm{I\kern-2.10007ptR}^{m\times n}$ be such that $\mathbf{A}+\mathbf{B}\mathbf{K}$ is Hurwitz. Let $\tilde{P}\in\mathbb{R}^{n\times n}$ be symmetric positive definite and $\mathbf{C}\in\mathbb{R}^{1\times n}$ such that

[TABLE]

Then the finite-dimensional system (13) is locally exponentially stable in [math] with a region of attraction given by

[TABLE]

where $\mathbf{D}$ is the minimal real number such that:

[TABLE]

Moreover, in this region of attraction, the function $V_{1}$ defined by $V_{1}(z)=z^{\top}Pz:=\mathbf{D}z^{\top}\tilde{P}z$, for all $z\in\mathrm{I\kern-2.10007ptR}^{n}$, decreases exponentially fast to [math] along the solutions to (13), i.e. there is a constant $\alpha>0$ so that

[TABLE]

Proof:

Define $P:=\mathbf{D}\tilde{P}$. Substituting this expression into (22) and (25) and multiplying it by $\mathbf{D}^{-1}$ to obtain the inequalities (38) and (42). If there exist $\tilde{P}$ and $\mathbf{C}$ as in Proposition 4 so that the LMI (38) holds, then according to Lemma 3, one can find $\mathbf{D}>0$ so that (42) holds as well.

The estimate for the region of attraction follows from Proposition 2. We choose the minimal $\mathbf{D}$ to obtain the largest estimate of a region of attraction (for given $\tilde{P}$ and $\mathbf{C}$). ∎

Remark 8

Lemma 3 can be helpful in finding the optimal $\mathbf{D}$ as it shows that the problem is to find the unique root of a piecewise analytic function.   $\circ$

III-E Enlarging the region of attraction using a dynamic controller

Consider a dynamic controller for system (13), 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. [41, Example 8.1] for a simple $2\times 2$ example). To be more specific, we consider an additional finite-dimensional state $z_{c}$ in $\mathrm{I\kern-2.10007ptR}^{n_{c}}$ (for a given integer $n_{c}$) and design matrices $\mathbf{K}_{1}$, $\mathbf{K}_{2}$, $\mathbf{A}_{c}$ and $\mathbf{B}_{c}$ of appropriate dimensions so that the (estimation of the) region of attraction of

[TABLE]

is larger (in the $z$-direction) than the estimate provided by Proposition 2 for (13). We rewrite the dynamics (44) by

[TABLE]

where

[TABLE]

with $\bar{\mathbf{A}}\in\mathrm{I\kern-2.10007ptR}^{(n+n_{c})\times(n+n_{c})}$, $\bar{\mathbf{B}}\in\mathrm{I\kern-2.10007ptR}^{(n+n_{c})\times m}$ and $\bar{\mathbf{K}}\in\mathrm{I\kern-2.10007ptR}^{m\times(n+n_{c})}$.

Remark 9

For given symmetric positive definite matrices $P\in\mathrm{I\kern-2.10007ptR}^{n\times n}$ and $\bar{P}\in\mathrm{I\kern-2.10007ptR}^{(n+n_{c})\times(n+n_{c})}$, the inclusion of the ellipsoid $\{z\in\mathrm{I\kern-2.10007ptR}^{n}:\;z^{\top}Pz\leqslant 1\}$ in the projection (onto $\mathrm{I\kern-2.10007ptR}^{n}$) of the ellipsoid $\{Z\in\mathrm{I\kern-2.10007ptR}^{n+n_{c}}:\;Z^{\top}\bar{P}Z\leqslant 1\}$ is equivalent to

[TABLE]

$\circ$

Using this remark and applying Proposition 2 to system (44), we have the following:

Proposition 5

Consider system (13) with a stabilizable pair $(\mathbf{A},\mathbf{B})\in\mathrm{I\kern-2.10007ptR}^{n\times n}\times\mathrm{I\kern-2.10007ptR}^{n\times m}$, and matrices $\bar{\mathbf{A}}$, $\bar{\mathbf{B}}$ and $\bar{\mathbf{K}}$ defined in (46). Choose a symmetric positive definite matrix $\bar{P}\in\mathbb{R}^{(n+n_{c})\times(n+n_{c})}$, a diagonal positive matrix $\bar{\mathbf{D}}\in\mathbb{R}^{m\times m}$ and $\bar{\mathbf{C}}\in\mathbb{R}^{m\times(n+n_{c})}$ such that (47) holds and

[TABLE]

[TABLE]

Then the finite-dimensional system (44) is locally exponentially stable in [math] with a region of attraction given by

${\cal A}_{\bar{P}}:=\{Z\in\mathrm{I\kern-2.10007ptR}^{n+n_{c}}:\;Z^{\top}\bar{P}Z\leqslant 1\}$.

Moreover, the projection of the ellipsoid $\{Z:\;Z^{\top}\bar{P}Z\leqslant 1\}$ onto $\mathrm{I\kern-2.10007ptR}^{n}$ is larger than the region of attraction ${\cal A}$ given by Proposition 2 in (26) for system (13).

Finally, in ${\cal A}_{\bar{P}}$, the function $V_{1}$ defined by

[TABLE]

decreases exponentially fast to [math] along the solutions to (44), i.e. there is a constant $\alpha>0$ so that for all $Z\in{\cal A}_{\bar{P}}$

[TABLE]

As in Section III-C, the inequalities (50) and (53) can be reformulated as linear matrix inequalities.

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

Now under the assumptions of Proposition 2, we provide a Lyapunov function $V:\ell_{2}(\mathrm{I\kern-2.10007ptN}^{*},\mathrm{I\kern-2.10007ptR})\to\mathrm{I\kern-2.10007ptR}_{+}$ for the system (12). To this end, we use the Lyapunov function $P$ provided by Proposition 2 for the finite-dimensional subsystem (13).

Proposition 6

Consider system (12) and assume $\mathbf{K}$ is chosen such that the subsystem (13) is locally exponentially stable in [math]. Assume further that $P\in\mathrm{I\kern-2.10007ptR}^{n\times n}$ is symmetric positive definite and such that for $\tilde{V}(z)=z^{\top}Pz,z\in\mathrm{I\kern-2.10007ptR}^{n}$ we have $\dot{\tilde{V}}(z)\leqslant-\alpha|z|^{2}$ on the set ${\cal A}:=\{z\in\mathrm{I\kern-2.10007ptR}^{n}:\tilde{V}(z)\leqslant 1\}$ along the solutions of (13). Then there exist $\gamma,C>0$ such that for $Q:\ell_{2}(\mathrm{I\kern-2.10007ptN}^{*},\mathrm{I\kern-2.10007ptR})\to\mathrm{I\kern-2.10007ptR}$ defined by $Q(w):=\sum_{j=n+1}^{\infty}w_{j}^{2}$ the function

[TABLE]

(with the identification $\pi_{n}w=z$) satisfies for all $w\in{\cal A}\times\ell_{2,j>n}$ that along the solutions of (12)

[TABLE]

In particular, it follows that $S\times\ell_{2,j>n}$ is a region of attraction of [math] for system (12).

Proof:

We write $w\in\ell_{2}(\mathrm{I\kern-2.10007ptN}^{*},\mathrm{I\kern-2.10007ptR})$ as $w=z+z^{\bot}$, where $z\in\mathrm{I\kern-2.10007ptR}^{n},z^{\bot}\in\ell_{2,j>n}$. Here we identify $\mathrm{I\kern-2.10007ptR}^{n}$ with the sequences with support in $\{1,\ldots,n\}$ and $\ell_{2,j>n}$ is the set of sequences in $\ell_{2}(\mathrm{I\kern-2.10007ptN}^{*},\mathrm{I\kern-2.10007ptR})$ which are [math] in the first $n$ entries.

This decomposition is unique. Due to (7), we have along the solutions to (12)

[TABLE]

Cauchy-Bunyakovsky-Schwarz inequality implies that

[TABLE]

Using Young’s inequality we have for all $\kappa>0$ that $2\|z^{\bot}(t)\|_{\ell_{2}}|z(t)|\leqslant\kappa\|z^{\bot}(t)\|_{\ell_{2}}^{2}+\frac{1}{\kappa}|z(t)|^{2}$. We proceed to:

[TABLE]

Therefore for any choice of $\gamma$ in (55), we have by an application of Proposition 2 along the solutions to (12) that

[TABLE]

Therefore selecting first $\kappa>0$ such that $2\eta-\kappa\|b^{\bot}\|_{\ell_{2}}\|\mathbf{K}\|>0$ and then selecting $\gamma>0$ such that $\alpha-\frac{\gamma}{\kappa}\|b^{\bot}\|_{\ell_{2}}\|\mathbf{K}\|>0$, we get the existence of $C>0$ such that

[TABLE]

provided that $z(t)$ lies in the ellipsoid ${\cal A}$, whatever the value of $z^{\bot}(t)$. As ${\cal A}\times\ell_{2,j>n}$ is an invariant subset of $\ell_{2}$ for a system (12), we obtain that $V$ is a Lyapunov function for (12), with a guaranteed region of attraction containing ${\cal A}\times\ell_{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(\cdot)$ in equation (1) be a constant function. With slight abuse of notation we will write $c(\cdot)=c=const$.

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

[TABLE]

on the domain $D(A)=H^{2}(0,L)\cap H^{1}_{0}(0,L)$ are given by

[TABLE]

and the eigenfunctions $e_{j}$, $j\in\mathrm{I\kern-2.10007ptN}^{*}$ of $(A,D(A))$, which form a basis of $L_{2}(0,1)$ are given by

[TABLE]

Next we estimate the region of attraction of system (1) for the following choice of parameters:

[TABLE]

where the $e_{j}$ are defined in (59). 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 $\mathbf{A},\mathbf{B}$ defined in (8) have the form

[TABLE]

As the diagonal entries in the matrix $\mathbf{A}$ are distinct, and all the components of $\mathbf{B}$ are nonzero, the pair $(\mathbf{A},\mathbf{B})$ is stabilizable. Different choices of the matrix $\mathbf{K}$ for the stabilizing feedback $u(t)=\mathbf{K}z(t)$ lead to different attraction rates and different regions of attraction. We demonstrate this with two examples.

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

First we choose the matrix $\mathbf{K}$ so that $\sigma(\mathbf{A}+\mathbf{B}\mathbf{K})=\{-1\}$, which results in

[TABLE]

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

• (i)

First we solve the inequality (38) together with the additional constraints: $\tilde{P}-{\tilde{P}}^{\top}=0$ and ${\tilde{P}}>0$. Additionally, we impose an optimality condition for $\mathbf{C}$:

[TABLE]

where $\cdot$ is the scalar product of vectors.

• (ii)

The idea behind (61) is to minimize the off-diagonal elements of the matrix $\tilde{M}_{2}$, which gives us at the second step the possibility to find large $\mathbf{D}$ (optimal for the given $\tilde{P}$, $\mathbf{C}$) satisfying the bilinear matrix inequality (42).

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

[TABLE]

Figure 1 shows an elliptic region of attraction (26), subject to $\tilde{P},\mathbf{C},\mathbf{D}$ given by (62) (in blue). Furthermore, in the same figure some trajectories are depicted obtained through direct simulation of (13). 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 (13). 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 (11).

Remark 10

(Importance of the optimality conditions) Different solutions $\tilde{P}$, $\mathbf{C}$, $\mathbf{D}$ of the matrix inequalities (38), (42) lead to very different estimates of a region of attraction (26) of the model (13). 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 (61). 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 (38), (42) 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.   $\circ$

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 $\tilde{P}$, $\mathbf{C}$, $\mathbf{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.   $\circ$

IV-B Choice 2: Putting the poles to $\{-0.1,-0.2\}$

Now let us choose the matrix $\mathbf{K}$ so that $\sigma(\mathbf{A}+\mathbf{B}\mathbf{K})=\{-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 $\mathbf{K}$ is:

[TABLE]

As in the previous simulation, we solve the LMI (38) subject to the additional optimality condition (61). The corresponding matrices $\tilde{P},\mathbf{C},\mathbf{D}$ are:

[TABLE]

As we see, with the choice of the stabilizing feedback (63), 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 (64) 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 and 2.

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 $\tilde{P}$, $\mathbf{C}$, $\mathbf{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 (13), 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.   $\circ$

V Extensions

V-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 [24].

V-A1 Saturation functions

Various definitions for the saturation map exist. Normwise saturations limit the norm of the input $u$, i.e. for a given norm $\|\cdot\|$ on $\mathrm{I\kern-2.10007ptR}^{n}$ we may consider

[TABLE]

Another physically motivated saturation map is the following (pointwise) $L_{\infty}$ saturation map, defined for all $x\in(0,L)$ by

[TABLE]

In this section, we depart from the saturation model in (1) and study instead a heat equation with pointwise saturation in each input channel:

[TABLE]

In a certain sense, in the equation (67) 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 (67), we are going to use the same stabilizing control (11). We now aim to provide an estimate for a region of attraction for (67). So we assume that $\mathbf{K}$ as in (11) is given and we consider the closed-loop system

[TABLE]

Representing this equation in the basis $e_{j}$, $j\in\mathrm{I\kern-2.10007ptN}^{*}$ as in Section II, we obtain the equations for the coordinates

[TABLE]

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

[TABLE]

Lemma 4

Let $r,k\in\mathrm{I\kern-2.10007ptR}$ and denote $\Delta(r,k)=r\textup{sat}(k)-\textup{sat}(rk)$. Then

[TABLE]

[TABLE]

Let $k\in\mathrm{I\kern-2.10007ptR}$ and $b\in L_{2}(0,L)$ be given. Then

[TABLE]

Moreover

[TABLE]

If $|k|\leqslant\ell$ and $\chi$ is the characteristic function of the set $U:=\{x\in[0,L]:|kb(x)|>\ell\}$ then

[TABLE]

If $b(\cdot)$ is essentially bounded, then

[TABLE]

where $\|b\|_{\infty}$ denotes the $L_{\infty}$ norm of the function $b$.

Proof:

The first claim follows as the assumption guarantee that $r\textup{sat}(k)=\textup{sat}(rk)=rk$. The claim in (70) is a direct consequence of the triangle inequality as

[TABLE]

The remaining claims follow immediately by applying the pointwise estimates (69) and (70). The claim (73) follows by applying (72) to $b\chi$, after noting that the complement of $U$ does not contribute to the norm of the left hand side. ∎

Proposition 7

Consider system (67) and assume all functions $b_{k}(\cdot)\in L_{\infty}([0,L])$, $k=1,\ldots,m$. Consider also the associated system (13) with a stabilizable pair $(\mathbf{A},\mathbf{B})\in\mathrm{I\kern-2.10007ptR}^{n\times n}\times\mathrm{I\kern-2.10007ptR}^{n\times m}$ and let $\mathbf{K}\in\mathrm{I\kern-2.10007ptR}^{m\times n}$ be such that $\mathbf{A}+\mathbf{B}\mathbf{K}$ is Hurwitz. Choose a symmetric positive definite matrix $P\in\mathbb{R}^{n\times n}$ as in Proposition 6. Consider the Lyapunov function candidate $V$ defined in (55).

Then there exist $\alpha,\beta>0$ such that if for $w=z+z^{\perp}\in\ell_{2}(\mathrm{I\kern-2.10007ptN}^{*},\mathrm{I\kern-2.10007ptR})$ we have $z\in{\cal A}_{\beta}:=\{z:\;z^{\top}Pz\leqslant\beta\}$ then

[TABLE]

In particular, $\{z:\;z^{\top}Pz\leqslant\beta\}\times X_{n}^{\bot}$ is a region of attraction for system (68).

Proof:

First we rewrite (67) by adding and subtracting the term $\sum b_{k}(x)\textup{sat}(u_{k}(t))$ to obtain

[TABLE]

Define

[TABLE]

For $j\in\mathrm{I\kern-2.10007ptN}^{*}$ we define $y_{j}(t):=\langle e_{j},\Delta(t)\rangle$, and let $y(t)\in\mathrm{I\kern-2.10007ptR}^{n}$ be the vector with components $y_{j}(t)$, for $j=1,\ldots,n$.

Considering the Lyapunov function $V$ defined in (55), we compute its time-derivative along the solutions to (68). Using (56), we get for all $w=z+z^{\bot}$ with $z^{\top}Pz\leqslant 1$ that

[TABLE]

Using (72) from Lemma 4, we obtain along the solutions to (67) and as long as $z^{\top}Pz\leqslant 1$,

[TABLE]

where $\lambda_{\max}$ denotes the maximal eigenvalue of the matrix $P$. Therefore for any positive values $\kappa$ and $\kappa^{\prime}$,

[TABLE]

Pick $\kappa>0$ and $\kappa^{\prime}>0$ such that $C-\frac{\lambda_{\max}}{\kappa}>\frac{3C}{4}$ and $C-\frac{\gamma}{\kappa^{\prime}}>\frac{C}{2}$. Due to (71), there exists $\beta>0$, such that for all $z$ in $\{z:\;z^{\top}Pz\leqslant\beta\}$, $\|\Delta(t)\|_{\infty}^{2}\leqslant\frac{C}{4(\lambda_{\max}\kappa+\gamma\kappa^{\prime})}|z(t)|^{2}$. We get, for all solutions to (67), as long as $z(t)$ is in $\{z:\;z^{\top}Pz\leqslant\beta\}$,

[TABLE]

Moreover, we have, along the solutions to (67),

[TABLE]

therefore, considering $V_{1}$ as previously defined, we may check that (75) holds following the same computation as for $\dot{V}$ along the solutions to (67), and using the fact the dynamics (80) in $X_{n}$ does not depend on the component in $X_{n}^{\bot}$. Therefore the set $\{z:\;z^{\top}Pz\leqslant\beta\}\times X_{n}^{\bot}$ is invariant along the dynamics to (67). With (79), we get that $\{z:\;z^{\top}Pz\leqslant\beta\}\times X_{n}^{\bot}$ is a region of attraction and $V$ is a Lyapunov function. ∎

V-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

[TABLE]

where $y_{d}$ is the (scalar) output of finite-dimensional dynamical system given by

[TABLE]

Here $x_{d}$ in $\mathrm{I\kern-2.10007ptR}^{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 (81) and the ODE (82) that is subject to a saturation map. Inspired by [31], we introduce the following change of variable:

[TABLE]

The PDE for $w$ then reads as:

[TABLE]

Please note that $b$ is a scalar function, and $d$ is a row vector function with $d(x)\in\mathrm{I\kern-2.10007ptR}^{1\times n_{d}}$, $x\in[0,L]$.

The boundary conditions for the variable $w$ take the form:

[TABLE]

The heat equation (83), (84) has to be analyzed along with the ODE (82a).

Performing similar computations as in Section II for the PDE (1) and, using the same notation for $w_{j}$ and $\lambda_{j}$, we get

[TABLE]

where $b$ and $d$ are defined for $x$ in $[0,L]$ in (83) and $b_{j}=\langle b(\cdot),e_{j}(\cdot)\rangle_{L_{2}(0,L)}$, $d_{j}=\langle d(\cdot),e_{j}(\cdot)\rangle_{L_{2}(0,L)}$, for $j=1,2,\ldots$.

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

[TABLE]

where $\mathbf{K}$ in $\mathrm{I\kern-2.10007ptR}^{1\times(n+n_{d})}$ is a row vector to be designed,

[TABLE]

and $\mathbf{A}{:=}\begin{pmatrix}A_{d}&0\\ D&\Lambda\end{pmatrix}\in\mathrm{I\kern-2.10007ptR}^{(n+n_{d})\times(n+n_{d})}$, where

[TABLE]

As the control input is scalar, we can apply Proposition 4 to system (85) instead of system (13), and get sufficient conditions for the estimation of the region of attraction of (85). Coming back to the infinite-dimensional systems (83) and (81), and, applying Proposition 1, we get sufficient conditions for an estimation of attraction region of (81):

Corollary 2

Consider system (13) with a stabilizable pair $(\mathbf{A},\mathbf{B})\in\mathrm{I\kern-2.10007ptR}^{n\times n}\times\mathrm{I\kern-2.10007ptR}^{n\times 1}$ and let $\mathbf{K}\in\mathrm{I\kern-2.10007ptR}^{1\times n}$ be such that $\mathbf{A}+\mathbf{B}\mathbf{K}$ is Hurwitz. Pick a symmetric positive definite $\tilde{P}\in\mathbb{R}^{(n+n_{d})\times(n+n_{d})}$ and a $\mathbf{C}\in\mathbb{R}^{1\times(n+n_{d})}$ such that (38) holds.

Then the finite-dimensional system (85) is locally exponentially stable in [math] with a region of attraction given by

[TABLE]

where $\mathbf{D}$ is the minimal real number such that (42) holds.

Moreover,

• (i)

(83) is locally exponentially stable with a region of attraction $\imath({\cal A})\times X_{n}^{\perp}$,

• (ii)

(81) 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.

-A Series expansion of solutions

Denoting the saturated nonlinearity in (1) by $f:\mathrm{I\kern-2.10007ptR}_{+}\to X$, we have for the mild solution of (4) that

[TABLE]

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 [15, Example A.1.13].

Let $\left(e_{j}\right)_{j\geqslant 1}$ be the Hilbert basis of $X$ given by the eigenfunctions of $A$. Then we can define

[TABLE]

and by [44, Corollary V.5.2] we may interchange the linear map $\langle\cdot,e_{j}\rangle$ with the integral to obtain

[TABLE]

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

[TABLE]

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

-B Compactness of the resolvent for Sturm-Liouville operators

The following discussion summarizes some results from [34], 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\leqslant p<\infty$) 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}:=\left(W^{1}_{0,2}\right)^{*}$ (the dual to $W^{1}_{0,2}=H^{1}_{0}$).

Recall that an operator $A\in 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\to X$ with a nonempty resolvent set $\rho(A)$, it is an easy consequence of the resolvent identity that the resolvent $R_{\lambda}(A)$ is compact for some $\lambda$ in the resolvent set, if and only if it is compact on the entire resolvent set, see [18, Theorem III.6.29].

Definition 3

We say that a closed densely defined linear operator $(A,D(A)):X\to X$ has a compact resolvent, if there exists a $\lambda\in\rho(A)$ so that $R_{\lambda}(A)$ is compact.

-B1 Some results from [34]

We note that in [34] the case $L=\pi$ is considered, which requires some rescaling to use their results.

Let $X:=L_{2}(0,L)$, $q\in W^{-1}_{2}(0,L)$, and define (for $y\in W^{1}_{1}(0,L)$) the quasiderivative

[TABLE]

where

[TABLE]

For $q\in W^{-1}_{2}(0,L)$ it holds that $Q\in X$, and for $q\in X$ it holds that $Q\in W^{1}_{2}(0,L)$.

Let $q\in L_{1}(0,L)$ be given. Consider the formal Sturm-Liouville operator $SL:X\to X$ defined by

[TABLE]

where $(0,L)\subset\mathrm{I\kern-2.10007ptR}$. In order to fully define the operator $SL$, we have to introduce its domain of definition. Following [34, Section 1.1] we define the maximal operator $L_{M}$, defined by

[TABLE]

The following result has been shown in [34, Theorem 1.5]:

Theorem 2

Let the operator $A$ be the restriction of $L_{M}$ to the domain

[TABLE]

where for $j=1,2$ it holds that

[TABLE]

where $a_{j1},a_{j2},b_{j1},a_{j2}$ are real numbers, for $j=1,2$.

Let $J_{\alpha\beta}$ be the determinant of the $\alpha$-th and $\beta$-th column of a matrix

[TABLE]

Then the operator $A$ has a nonempty resolvent set $\rho(A)$, has a compact resolvent and discrete spectrum, if one of the following conditions holds:

1. (i)

$J_{42}\neq 0$, 2. (ii)

$J_{42}=0$, $J_{14}+J_{32}\neq 0$, 3. (iii)

$J_{42}=J_{14}=J_{32}=0$, $J_{12}+J_{34}=0$, $J_{13}\neq 0$.

Remark 13

Compactness of the resolvent of $A$ is not mentioned in the formulation of [34, Theorem 1.5], but its compactness was shown in the proof.   $\circ$

Remark 14

If one of conditions (i)–(iii) holds, then the boundary conditions (90) are called Birkhoff-regular.   $\circ$

-B2 Application to the system in Section II

Consider the operator $A$ as in Section II, i.e.

[TABLE]

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\in[1,+\infty)$. Then $f\in W^{1}_{p}(0,L)$ if and only if $f$ is equal a.e. to an absolutely continuous function, the derivative $f^{\prime}$ exists a.e. and is an element of $L_{p}(0,L)$.

First we show

Lemma 6

The domain $D(A)$ of the operator (92) has the form (89) for $U_{1}(y)=y(0)$ and $U_{2}(y)=y(L)$.

Proof:

Let $AC((0,L),\mathrm{I\kern-2.10007ptR})$ be the space of absolutely continuous functions from $(0,L)$ to $\mathrm{I\kern-2.10007ptR}$. In view of Lemma 5 the set $D(A)$ can be rewritten as

[TABLE]

Now consider the domain defined in (89).

As we restrict our attention to the case when $q\in X$, for any $y\in W^{1}_{1}(0,L)$ it holds that $Q(\cdot)y(\cdot)\in W^{1}_{1}(0,L)$, and thus $y^{[1]}\in W^{1}_{1}(0,L)$ if and only if $\frac{dy}{dx}\in W^{1}_{1}(0,L)$.

Hence the domain $D(L_{M})$ can be equivalently written as

[TABLE]

Using again Lemma 5, we see that

[TABLE]

Furthermore, as $f\in D(L_{M})$ is absolutely continuous on $(0,L)$, and $q\in X$, it holds that $qf\in X$ and it holds that $SL(f)\in X$ if and only if $\frac{d^{2}f}{dx^{2}}\in X$. As $X\subset L_{1}(0,L)$, we can finally restate $D(L_{M})$ as

[TABLE]

Using the boundary conditions $U_{1}(y)=y(0)$ and $U_{2}(y)=y(L)$, we see that the set (89) is precisely (93). ∎

In view of Lemma 6, we can use Theorem 2 to study the spectral properties of the Sturm-Liouville operator $A$.

Proposition 8

The operator (92) has a compact resolvent.

Proof:

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 2 have the form: $a_{11}=1$, $b_{21}=1$, and all other entries of a matrix in (91) are zeros. By item (iii) of Theorem 2 the boundary conditions are Birkhoff-regular and the operator $A$ has a compact resolvent. ∎