A Semi-Definite Programming Approach to Stability Analysis of Linear Partial Differential Equations

TL;DR

This paper introduces a semi-definite programming method for stability analysis of linear 1-D PDEs with polynomial data, enabling systematic construction of Lyapunov certificates for various PDE systems.

Contribution

It develops a novel Lyapunov-based approach that reduces PDE stability verification to solving polynomial problems via semi-definite programming.

Findings

Successfully applied to different PDE types

Provides a systematic computational framework

Enables stability certification for complex PDEs

Abstract

We consider the stability analysis of a large class of linear 1-D PDEs with polynomial data. This class of PDEs contains, as examples, parabolic and hyperbolic PDEs, PDEs with boundary feedback and systems of in-domain/boundary coupled PDEs. Our approach is Lyapunov based which allows us to reduce the stability problem to the verification of integral inequalities on the subspaces of Hilbert spaces. Then, using fundamental theorem of calculus and Green's theorem, we construct a polynomial problem to verify the integral inequalities. Constraining the solution of the polynomial problem to belong to the set of sum-of-squares polynomials subject to affine constraints allows us to use semi-definite programming to algorithmically construct Lyapunov certificates of stability for the systems under consideration. We also provide numerical results of the application of the proposed method on different types of PDEs.

Tables

Table 1. TABLE I: Maximum λ ∈ ℝ 𝜆 ℝ \lambda\in\mathbb{R} for which problem ( 48 ) is feasible for Examples  1 1 1 - 3 3 3 in Equations ( 3 )-( 5 ), respectively, as a function of polynomial degree d 𝑑 d . Here, Inf. denotes infeasibility.

	$d=2$	$4$	$6$	$8$
Example $1$ Eqn. (3)	$\lambda =3.263$	$3.263$	$3.409$	$3.409$
Example $2$ Eqn. (4)	Inf.	Inf.	$\lambda =0.999$	$0.999$
Example $3$ Eqn. (5)	$\lambda =0.999$	$0.999$	$0.999$	-

Table 2. TABLE II: Computer run time (in seconds) for performing the search for variables which solve problem ( 48 ) for Examples  1 1 1 - 3 3 3 in Equations ( 3 )-( 5 ), respectively.

	$d=2$	$d=4$	$d=6$	$d=8$
Example $1$ Eqn. (3)	$4.415$	$8.111$	$20.273$	$48.080$
Example $2$ Eqn. (4)	$7.085$	$12.796$	$49.751$	$141.297$
Example $3$ Eqn. (5)	$25.663$	$113.708$	$360.876$	-

Equations

$$w_{\alpha }(x)=$$

$$w_{\alpha }^b=$$

$${\bar{w}}_{\alpha }(x)=$$

$${\partial }_x{w}_{\alpha -1}(x)=N_{\partial }{w}_{\alpha }(x),$$

$$w_{\alpha -1}(0)=N_0{w}_{\alpha }^b,$$

$$\displaystyle\mathcal{H}^{\alpha}\left([0,1];\mathbb{R}^{\beta}\right)=\Big{\{}w:[0,1]$$

$$\text{absolutely continuous on }[0,1]\text{and}$$

$$\displaystyle\int_{0}^{1}\left(\partial_{x}^{\alpha}w(x)\right)^{T}\left(\partial_{x}^{\alpha}w(x)\right)dx<\infty\Big{\}}.$$

$$\Vert w{\Vert }_{{\mathcal{H}}^{\alpha }}=\sqrt{{⟨w_{\alpha },w_{\alpha }⟩}_{{\mathcal{H}}^{\alpha }}}=\sqrt{{\int }_0^1{w}_{\alpha }(x{)}^T{w}_{\alpha }(x)dx},$$

$$Y_{q(\alpha ,\beta ,d)}(x,y)=I_{\beta (\alpha +1)}\otimes z(x,y)\in {\mathbb{R}}^{q(\alpha ,\beta ,d)\times \beta (\alpha +1)},$$

$$z(x,y)={\left[\begin{array}{cccccc}1& x& y& x^2& xy& y^2\end{array}\right]}^T.$$

$$\Gamma \left[K\right]=\{\begin{array}{ll}K(x,y),& x\ge y\\ K(y,x{)}^T,& y>x\end{array},$$

$${\int }_{\Delta }v(x{)}^T\Gamma \left[K\right]v(y)d\Delta =$$

$$+{\int }_0^1{\int }_x^{1}v(x{)}^{T}K(y,x{)}^{T}v(y)dydx.$$

$${\partial }_{t}w(x,t)=A(x)w_{\alpha }(x,t),$$

$$\mathcal{B}=$$

$${\partial }_{t}v(x,t)=$$

$$v(0,t)=$$

$$w(x,t)=$$

$$A(x)=$$

$${\partial }_t{v}_1(x,t)=\left(\lambda -1\right)\left(x^2+1\right){\partial }_x{v}_1(x,t)+(x-3)v_2(x,t),$$

$${\partial }_t{v}_2(x,t)=\left(x+1\right){\partial }_x{v}_2(x,t),$$

$$v_1(0,t)-3v_2(0,t)=0,v_2(1,t)=0.$$

$$w(x,t)=$$

$$A(x)=$$

$$F=$$

$${\partial }_t^{2}v(x,t)+{\partial }_x^{4}v(x,t)=0,$$

$${\partial }_x^{2}v(0,t)-\frac{1}{1-\lambda }{\partial }_{xt}v(0,t)=0,$$

$${\partial }_x^{3}v(0,t)+(1-\lambda ){\partial }_{t}v(0,t)=0,{\partial }_x^{2}v(1,t)=0,$$

$$v(1,t)=0.$$

$${\partial }_t\left[\begin{array}{c}{\partial }_x^{2}v(x,t)\\ {\partial }_{t}v(x,t)\end{array}\right]=\left[\begin{array}{c}{\partial }_t{\partial }_x^{2}v(x,t)\\ -{\partial }_x^{4}v(x,t)\end{array}\right],$$

$${\partial }_x^{2}v(0,t)-\frac{1}{1-\lambda }{\partial }_{xt}v(0,t)=0,$$

$${\partial }_x^{3}v(0,t)+(1-\lambda ){\partial }_{t}v(0,t)=0,{\partial }_x^{2}v(1,t)=0,$$

$${\partial }_{t}v(1,t)=0.$$

$$w(x,t)=$$

$$A(x)=$$

$$F=$$

$$\text{mcz}V(w)=$$

$$+\frac{1}{2}{\int }_{\Delta }w(x,t{)}^T\Gamma \left[\bar{T}\right]w(y,t)d\Delta ,$$

$$\text{mcz}{V}_d(w)=-{\partial }_t\text{mcz}V(w)-2\delta \text{mcz}V(w),$$

$$\text{mcz}{V}_d(w)=$$

Full text

I INTRODUCTION

Temporal evolution of processes involving spatially distributed physical quantities requires Partial Differential Equations (PDEs) to produce accurate models for analysis and control [19], [3]. Analysis of PDE systems brings forth technical challenges relative to Ordinary Differential Equations (ODEs) due to the infinite dimensional nature of the state-space. Owing to the maturity of research on ODEs, a traditional method of analysis is to reduce PDE systems, using spatial or spectral methods, to a system of ODEs which can then be analyzed with relative ease [6], [9].

Recently, various direct methods, i.e., methods without finite-dimensional approximation, have been developed for analysis and control of PDEs. For controller synthesis, backstepping is an oft used method wherein the problem of stabilizing boundary feedback design is reduced to a hyperbolic PDE whose solution can be obtained either analytically or numerically [10]. Similarly, one may use Lyapunov’s second method for infinite-dimensional systems to establish stability [4], [12]. Lyapunov’s method requires that the structure of a Lyapunov Functional (LF) be chosen a priori. In this paper we choose quadratic LF candidates of the form used in our previous work in [7]. The reason for this choice, as we established in [7], is that such LF candidates are successful in providing stability certificates for parabolic PDEs. Additionally, also in [7], we showed that such LFs are certificates of stability for parabolic PDEs actuated by backstepping boundary feedback.

The LF candidates we consider lead to Lyapunov inequalities defined by integral expressions containing quadratic terms in the state defined on one-dimensional (1-D) and two-dimensional (2-D) domains and their respective boundaries. Moreover, the verification of these inequalities has to be performed on the spaces of functions defined by the boundary conditions of the PDE whose stability we wish to establish. The interplay between the 2-D domain/1-D domain and the boundary values have to be thus accounted for with the Fundamental Theorem of Calculus (FTC) and Green’s theorem, which are used in a similar way to [16].

I-A Contribution

The paper presents a convex optimization based method for the verification of integral inequalities on (subspaces of) Hilbert spaces and its application to stability analysis of PDE systems with polynomial data. The proposed method contains previous approaches by the authors, see, for e.g. [7], [17], thus encompassing a larger class of PDEs.

In order to accomplish our goals we first reduce the problem of stability analysis to the verification of two integral inequalities on (subspaces of) Hilbert spaces. To verify these inequalities, we use FTC and Green’s theorem to construct polynomial equations whose solution, if exists, verifies the Lyapunov integral inequalities. Finally, we show that this polynomial equation may be solved by searching for Sum-of-Squares (SOS) polynomials under affine constraints thus rendering the problem of searching for the solution as a Semi-Definite Program (SDP) [2, Chapter 3], [18], a convex optimization problem.

I-B Outline

The paper is organized as follows. In Section II we present the problem statement and formulate the stability analysis problem as the problem of verification of two integral inequalities. In Section III we use positive semi-definite polynomial matrices to characterize a class of non-negative/positive integral inequalities on Hilbert spaces. In Section IV we use FTC and Green’s theorem to construct integral terms which relate the domains on which the LF candidates are defined to their boundaries, which we call slack integrals. Finally, in Section V we present the main result and numerical experiments.

I-C Notation

We denote by $\overline{\Omega}$ the set $\{(x,y)\in\mathbb{R}^{2}:0\leq y\leq x\leq 1\}$, $\underline{\Omega}=\{(x,y)\in\mathbb{R}^{2}:0\leq x\leq y\leq 1\}$ and $\Delta=[0,1]\times[0,1]$. In the following definitions $\alpha,\beta\in\mathbb{N}$. For $w:[0,1]\rightarrow\mathbb{R}^{\beta}$, $w\in\mathcal{C}^{\alpha}([0,1])$, we denote by $\partial_{x}^{i}w(x)$ the $i$-th derivative of $w$ and

[TABLE]

Thus, $w_{\alpha}:[0,1]\rightarrow\mathbb{R}^{\beta(\alpha+1)}$, $w_{\alpha}^{b}\in\mathbb{R}^{2\beta\alpha}$ and $\bar{w}_{\alpha}:[0,1]\rightarrow\mathbb{R}^{\beta(3\alpha+1)}$. We denote by $N_{\partial},N\in\mathbb{N}^{\beta\alpha\times\beta(\alpha+1)}$ and $N_{0},N_{1}\in\mathbb{N}^{\beta\alpha\times 2\beta\alpha}$ the matrices satisfying

[TABLE]

We denote

[TABLE]

We also write $\mathcal{L}_{2}\left([0,1];\mathbb{R}^{\beta}\right)=\mathcal{H}^{0}([0,1])$. The space $\mathcal{H}^{\alpha}\left([0,1];\mathbb{R}^{\beta}\right)$ is endowed with the standard inner product and norm

[TABLE]

and the space $\mathcal{L}_{2}\left([0,1];\mathbb{R}^{\beta}\right)$ has the norm and inner product $\lVert{\cdot}\rVert_{\mathcal{L}_{2}}=\lVert{\cdot}\rVert_{\mathcal{H}^{0}}$ and $\left\langle\cdot,\cdot\right\rangle_{\mathcal{L}_{2}}=\left\langle\cdot,\cdot\right\rangle_{\mathcal{H}^{0}}$, respectively.

For any $m,n\in\mathbb{N}$, we denote by $0_{m,n}$ the matrix of zeros of dimensions $m$-by-$n$ and $0_{m}$ when $m=n$. Similarly, we denote by $I_{n}$ the identity matrix of dimensions $n$-by-$n$. The set of symmetric matrices of dimension $n$-by$n$ is donated by $\mathbb{S}^{n}$. For any square matrix $Q$, we denote $He(Q)=\frac{1}{2}\left(Q+Q^{T}\right)$.

We denote by $\mathcal{S}^{n}[x]$ and $\mathcal{S}^{n}[(x,y)]$ the sets of symmetric polynomial matrices of size $n$-by-$n$ in variables $x$, and $x$ and $y$, respectively. Similarly, we denote by $\mathcal{R}^{m\times n}[x]$ and $\mathcal{R}^{m\times n}[(x,y)]$ the sets of real polynomial matrices of size $m$-by-$n$ in variables $x$ and $y$. We denote by $\Sigma^{n}[x]$ the set of Sum-of-Squares (SOS) polynomials which, for all $x\in[0,1]$, belong to $\mathbb{S}^{n}$. polynomials in the variable $x$. Note that, by definition, an SOS polynomial is non-negative for all $x\in\mathbb{R}$ [2, Chapter 3]. We say that a polynomial $S\in\mathcal{S}^{n}[x]$ is positive semi-definite on $[0,1]$ if $S(x)\succeq 0$, for all $x\in[0,1]$. Given $\alpha,\beta,d\in\mathbb{N}$, we define

[TABLE]

where $q(\alpha,\beta,d)=\frac{1}{2}\beta(\alpha+1)(d+2)(d+1)$, $z(x,y)\in\mathbb{R}^{\frac{1}{2}(d+2)(d+1)}$ is the vector of monomials in $x$ and $y$ up to degree $d$ and $\otimes$ denotes the Kronecker product. For example, for $d=2$, we have

[TABLE]

For any bivariate function $K:\overline{\Omega}\rightarrow\mathbb{R}^{n\times n}$, we define the linear map

[TABLE]

thus satisfying, for any $v:[0,1]\rightarrow\mathbb{R}^{n}$,

[TABLE]

II PROBLEM STATEMENT

Consider the following class of linear PDEs

[TABLE]

We now provide a few examples of PDEs which can be cast in the form of (2). Each of the following PDEs is parameterized by a positive scalar $\lambda$.

II-1 Example $1$

Consider the following parabolic PDE with distributed coefficients

[TABLE]

This PDE may be set in the form of (2) with $\alpha=2,\beta=1$ and

[TABLE]

II-2 Example $2$

Now consider the following system of hyperbolic PDEs coupled in-domain and at the boundaries

[TABLE]

This PDE may be set in the form of (2) with $\alpha=1,\beta=2$ and

[TABLE]

II-3 Example $3$

Finally, consider the following Euler-Bernoulli beam model

[TABLE]

We may re-write (5) as

[TABLE]

With this representation, we may write (5) in the form of (2) with $\alpha=2,\beta=2$ and

[TABLE]

We wish to establish the stability of (2) by constructing Lyapunov functional (LF) certificates of exponential stability. In particular, we wish to construct LFs of the form

[TABLE]

where $T_{b}\in\mathcal{S}^{\beta}[x]$ and $\bar{T}\in\mathcal{R}^{\beta\times\beta}[(x,y)]$. As stated in the Introduction, parabolic PDE systems with boundary backstepping based control laws admit LF certificates of stability which have the same structure as (6) [7]. The numerical results in [7] also indicate that such LFs are not conservative for spatially distributed scalar parabolic PDEs.

Let us now formulate conditions for the exponential stability of (2). We begin by computing the time derivative of (6) along the trajectories of (2). Consider a scalar $\delta>0$ and define

[TABLE]

which can be put in the form (see the Appendix)

[TABLE]

where

[TABLE]

The following theorem casts the verification of exponential stability of (2) as the verification of two integral inequalities.

Theorem 1

Given the PDE (2), suppose there exist positive scalars $\mu,\delta$ and polynomial matrices $T_{b}\in\mathcal{S}^{\beta}[x]$ and $\bar{T}\in\mathcal{R}^{\beta\times\beta}[(x,y)]$ such that

[TABLE]

where $\mcz{V}(w)$ and $\mcz{V}_{d}(w)$ are defined in (6) and (II-3), respectively.

Then, there exists a positive scalar $\kappa$ such that the solution of (2) satisfies

[TABLE]

Proof:

Let us choose the LF candidate $\mcz{V}(w)$. Since the condition in (9) holds, we have that there exists a positive scalar $\theta$ such that

[TABLE]

where the upper bound is a consequence of $T_{b}$ and $\bar{T}$ being polynomial matrices defined on bounded domains.

Now, for this LF candidate, we have from (7) that

[TABLE]

and since (10) holds for all $w\in\mathcal{B}$ and thus for $w$ that solve (2), we have that

[TABLE]

Integrating this expression in time produces

[TABLE]

Using (12) produces

[TABLE]

Then, taking the square root we conclude that (11) holds with $\kappa=\sqrt{\theta/\mu}$. ∎

We have reduced the problem of stability analysis of (2) to the verification of the integral inequalities in (9) and (10). Thus, the remainder of the work considers the following problem:

Problem: Verify if $\mcz{V}(w)$ is strictly positive, i.e., (9) is verified on $\mathcal{L}_{2}\left([0,1];\mathbb{R}^{\beta}\right)$, and $\mcz{V}_{d}(w)$ is non-negative on $\mathcal{B}\subset\mathcal{H}^{\alpha}\left([0,1];\mathbb{R}^{\beta}\right)$.

We verify (9) by generalizing the methods proposed in [7] and [17]. Namely, we test the positive semi-definiteness of a polynomial matrix associated to $T_{b}$ and $\bar{T}$. Since (10) requires $\mcz{V}_{d}(w)\geq 0$ only on the subset $\mathcal{B}$ of $\mathcal{H}^{\alpha}\left([0,1];\mathbb{R}^{\beta}\right)$, it calls for a different formulation than the one adopted to verify (9). In this case we first construct a set of integral terms $\mcz{S}(w)$ which have the same structure as $\mcz{V}_{d}(w)$ and satisfy $\mcz{S}(w)=0$, for all $w\in\mathcal{B}$. We refer to such expressions as slack integrals. We then check for the positive semi-definiteness of a polynomial matrix associated to $\mcz{V}_{d}(w)+\mcz{S}(w)$ which guarantees $\mcz{V}_{d}(w)+\mcz{S}(w)\geq 0$, for all $w\in\mathcal{H}^{\alpha}\left([0,1];\mathbb{R}^{\beta}\right)$.

The construction of the desired slack integrals is presented in Section IV. These integrals generalize the results presented in [13], [14] for time-delay systems and in [17] for PDE analysis using (6) with $\bar{T}=0$.

III POSITIVE/NON-NEGATIVE INTEGRAL INEQUALITIES ON HILBERT SPACES

In this section we construct a set of positive/non-negative integral inequalities which are parameterized by Positive Semi-Definite (PSD) polynomial matrices. The results provided are a generalization of [7, Theorem 1]. Throughout this section we will use the polynomial matrix $Y_{q(\alpha,\beta,d)}(x,y)$ defined in (1).

We begin by constructing non-negative integral inequalities on $\mathcal{H}^{\alpha}\left([0,1];\mathbb{R}^{\beta}\right)$, given $\alpha,d\in\mathbb{N}$, that have the same form as $\mcz{V}_{d}(w)$ in (II-3). Let us define

[TABLE]

where $R_{b}\in\mathcal{S}^{\beta(3\alpha+1)}[x]$ and

[TABLE]

for some (polynomial) matrices $R_{12},R_{13}\in\mathcal{R}^{\beta(\alpha+1)\times q(\alpha,\beta,d)}[x]$, $R_{22},R_{33}\in\mathbb{S}^{q(\alpha,\beta,d)}$ and $R_{23}\in\mathbb{R}^{q(\alpha,\beta,d)\times q(\alpha,\beta,d)}$, which also define

[TABLE]

The following theorem states the conditions on the polynomial matrix $R(x)$ whose submatrices define $R_{b}$ and $\bar{R}$ such that (14) is non-negative on $\mathcal{H}^{\alpha}\left([0,1];\mathbb{R}^{\beta}\right)$.

Theorem 2

Given (polynomial) matrices which define $\mcz{R}(w)$ in (14), if

[TABLE]

where $R(x)$ is defined in (20), then

[TABLE]

Proof:

Throughout this proof, for notational brevity, we write $q$ in place of $q(\alpha,\beta,d)$.

If we define

[TABLE]

then

[TABLE]

since $R(x)\succeq 0$, for all $x\in[0,1]$. Simplifying the expression we obtain

[TABLE]

where

[TABLE]

Applying Lemma .1 to $\Theta_{1}$ with $K_{1}(x,y)=2R_{12}(x)Y_{q}(x,y)$ and $K_{2}(x,y)=2R_{13}(x)Y_{q}(x,y)$ produces

[TABLE]

Applying Lemma .2 to $\Theta_{2}$ with

[TABLE]

produces

[TABLE]

Similarly, applying Lemma .2 to $\Theta_{3}$ with

[TABLE]

produces

[TABLE]

Substituting (25)-(27) into (24) produces

[TABLE]

Then, (23) completes the proof. ∎

We now present a corollary in which we construct strictly positive integral inequalities on $\mathcal{L}_{2}\left([0,1];\mathbb{R}^{\beta}\right)$ that have the same form as $\mcz{V}(w)$ in (6). This corollary corresponds to setting $\alpha=0$ in Theorem 2.

Given $w\in\mathcal{L}_{2}\left([0,1];\mathbb{R}^{\beta}\right)$ and $\beta,d\in\mathbb{N}$, let us define

[TABLE]

where $T_{b}\in\mathcal{S}^{\beta}[x]$ and

[TABLE]

for any $T_{12},T_{13}\in\mathcal{R}^{\beta\times q(0,\beta,d)}[x]$, $T_{22},T_{33}\in\mathbb{S}^{q(0,\beta,d)}$ and $T_{23}\in\mathbb{R}^{q(0,\beta,d)\times q(0,\beta,d)}$, which also define

[TABLE]

Corollary 1

Given (polynomial) matrices which define $\mcz{T}(w)$ in (28), if there exists a positive scalar $\epsilon$ such that

[TABLE]

where $T(x)$ is defined in (30), then,

[TABLE]

Proof:

Following the same steps as for the proof of Theorem 2, it can be shown that

[TABLE]

where

[TABLE]

Then, (32) holds since (31) holds. ∎

IV SLACK INTEGRALS

In Section II we cast the stability of (2) as a test of positivity and non-negativity of integral inequalities in (9) and (10). In Section V we will show that Corollary 1 is sufficient to verify (9). Theorem 2 is too conservative to verify (10) because it enforces $\mcz{V}_{d}(w)\geq 0$ on the entire space $\mathcal{H}^{\alpha}\left([0,1];\mathbb{R}^{\beta}\right)$, while we are only interested in non-negativity over the subset $\mathcal{B}\subset\mathcal{H}^{\alpha}\left([0,1];\mathbb{R}^{\beta}\right)$. As discussed in Section II, however, to solve this problem we construct slack integrals $\mcz{S}(w)$ which we defined as integral expressions with the same structure as $\mcz{V}_{d}(w)$ and satisfies

[TABLE]

Then, we may use Theorem 2 to test if

[TABLE]

which, if true, would imply that $\mcz{V}_{d}(w)\geq 0$, $\forall w\in\mathcal{B}$, owing to (33).

We will construct slack integrals using quadratic forms of the Fundamental Theorem of Calculus (FTC) and Green’s theorem.

Lemma 1** (FTC quadratic form)**

For any $K_{1}\in\mathcal{R}^{\beta\alpha\times\beta\alpha}[x]$ and $K_{2}\in\mathcal{R}^{\beta\alpha\times 2\beta\alpha}[x]$ the following identity holds

[TABLE]

where

[TABLE]

Proof:

The identity (34) is established by expanding

[TABLE]

with

[TABLE]

∎

Next, we present the quadratic form of the Green’s theorem. The proof of the following lemma is provided in the Appendix.

Lemma 2** (Green’s theorem quadratic form)**

For any $H_{1},H_{2}\in\mathcal{R}^{\beta\alpha\times\beta\alpha}[(x,y)]$, the following identity holds

[TABLE]

where

[TABLE]

In the following lemma we formulate an integral equation that holds on the set $\mathcal{B}$ given in (2b). The proof of the following lemma is provided in the Appendix.

Lemma 3

Given $F\in\mathbb{R}^{\beta\alpha\times 2\beta\alpha}$ the following identity holds true for all $w\in\mathcal{B}$:

[TABLE]

where

[TABLE]

for any $B_{1}\in\mathcal{R}^{\beta(\alpha+1)\times\beta\alpha}[x]$ and $B_{2}\in\mathcal{R}^{2\beta\alpha\times\beta\alpha}[x]$.

We now use the results in Lemmas 1-3 to formulate slack integrals on the set $\mathcal{B}\in\mathcal{H}^{\alpha}\left([0,1];\mathbb{R}^{\beta}\right)$.

Let us define

[TABLE]

where $K_{b}$ is parameterized by $K_{1}$ and $K_{2}$ as in Lemma 1, $H_{b}$ and $\bar{H}$ are parameterized by $H_{1}$ and $H_{2}$ as in Lemma 2 and $B_{b}$ is parameterized by $B_{1}$ and $B_{2}$, and the matrix $F$ which defines the set $\mathcal{B}$ as in Lemma 3.

We now state the main result of this section.

Theorem 3

Given matrix $F$ which defines the set $\mathcal{B}$ in (2b), the following identity holds true

[TABLE]

where $\mcz{S}(w)$ is parameterized by any $K_{i}$, $H_{i}$ and $B_{i}$, $i\in\{1,2\}$, as in in (35).

Proof:

We begin by considering the following decomposition

[TABLE]

where

[TABLE]

From Lemmas 1 and 2 we have that

[TABLE]

From Lemma 3 we have that

[TABLE]

Therefore, from (37)-(39) we conclude that the expression in (36) holds for all $w\in\mathcal{B}$. ∎

V MAIN RESULT

In this section we use the results formulated in Sections III-IV to construct a method of verifying the stability of PDE (2). Let us proceed with the following.

Theorem 4

Consider the PDE (2). For any given $d\in\mathbb{N}$, positive scalars $\epsilon,\delta$ and polynomial $Y_{q(\alpha,\beta,d)}\in\mathcal{R}^{q(\alpha,\beta,d)\times\beta(\alpha+1)}[(x,y)]$ defined in (1), suppose there exist:

[TABLE]

such that

[TABLE]

where $\omega(x)=x(1-x)$, $T(x)$ and $R(x)$ are defined in (30) and (20), respectively, $K_{b}$ is defined using $K_{i}$ as in Lemma 1, $H_{b}$ and $\bar{H}$ are defined using $H_{i}$ as in Lemma 2, $B_{b}$ is defined using $B_{i}$ and $F$ as in Lemma 3, and polynomials $U_{b}$ and $\bar{U}$ are defined using $\delta$ in (II-3) and with $T_{b}$ and $\bar{T}$ defined in (28).

Then, (2) is exponentially stable.

Proof:

Since the polynomial matrix in (41a) is Sum-of-Squares (SOS), $S_{T}$ is SOS in (41b) and $\omega(x)\geq 0$, for all $x\in[0,1]$, using the property that a SOS polynomial is non-negative for all $x\in\mathbb{R}$, we conclude that $T(x)$ satisfies (31). Therefore, $\mcz{T}(w)\geq\epsilon\lVert{w}\rVert_{\mathcal{L}_{2}}$, for all $w\in\mathcal{L}_{2}\left([0,1];\mathbb{R}^{\beta}\right)$. Moreover, from (6) we have that $\mcz{V}(w)=\frac{1}{2}\mcz{T}(w)$ and thus

[TABLE]

with $\mu=\frac{1}{2}\epsilon$. Similarly, since (41c)-(41d) hold, we conclude from Theorem 2 that

[TABLE]

Now, let us define

[TABLE]

where

[TABLE]

Since (41e)-(41g) hold, we deduce that

[TABLE]

Therefore, we get

[TABLE]

From the definition of $\mcz{P}(w)$ in (44) it is clear that $\mcz{P}(w)=\mcz{V}_{d}(w)+\mcz{S}(w)-\mcz{R}(w)$, where $\mcz{V}_{d}(w)$ is defined in (II-3) and $\mcz{S}(w)$ is defined in (35). Thus, using (45) we obtain

[TABLE]

Using (43) the previous expression may be reduced to

[TABLE]

Now, from Theorem 3 we have that $\mcz{S}(w)=0$, for all $w\in\mathcal{B}$, thus, from (46) we deduce that

[TABLE]

Since (42) and (47) hold, we apply Theorem 1 to complete the proof. ∎

Since polynomials are closed under differentiation and integration, we have that $K_{b}$ and $\bar{H}$ containing $\partial_{x}K_{1}(x)$, $\partial_{x}K_{2}(x)$, $\partial_{y}H_{1}(x,y)$ and $\partial_{x}H_{2}(x,y)$ in Lemmas 1-2 are polynomials. Moreover, since the polynomial $Y_{q(\alpha,\beta,d)}(x,y)$ is fixed, the terms $R_{b}$, $\bar{R}$, $T_{b}$ and $\bar{T}$ in Theorem 2 and Corollary 1 are polynomials affine in their respective (polynomial) matrices. Therefore, the conditions in (41) are simply either, $1$) a verification of the membership of polynomial matrices in the set of SOS polynomials as in (41a)-(41f), or, $2$) enforcement of affine constraints on the polynomial variables as in (41g). Indeed, the problem of searching for SOS polynomials subject to affine constraints is a Semi-Definite Program (SDP), [2, Chapter 3], [18]. We are then interested in solving

[TABLE]

The numerical implementation is performed by constructing the underlying SDP for (48) by using the freely available packages SOSTOOLS [1] or YALMIP [11]. Then, the associated SDP is solved, for example, using SeDuMi [15] or SDPA [20] solvers.

V-A Numerical Examples

We now determine the stability of PDEs (3), (4) and (5) presented in Section II by solving (48) for each of these systems. Recall that each of these systems are parameterized by a positive scalar $\lambda\in\mathbb{R}$. The parameter $\lambda$ can alter the stability properties of these systems, and thus, allows us to measure the effectiveness of the proposed methodology.

We run the numerical codes for the example PDEs with $\epsilon=10^{-3}$, $\delta=10^{-4}$ and maximum polynomial degree allowed by the memory resources of $8$ giga bytes of RAM. In the following examples, we search for the maximum $\lambda\in\mathbb{R}$ for which the conditions of problem (48) are feasible using a bisection search with a resolution of $10^{-3}$. These results are provided in Table I. The discussion of the results is provided below.

Example $1$

Using finite-differences with $1500$ uniformly distributed spatial points we approximate that the parabolic PDE in (3) is stable for $\lambda<3.412$. As illustrated in Table I, the maximum $\lambda$ for which the problem (48) is feasible is $\lambda=3.409$ which is $99.91\%$ of the value of $3.412$ obtained via the finite-difference approximation.

Example $2$

Using [5, Lemma 3.1] it can be established that the coupled first order hyperbolic PDE in (4) is exponentially stable for $\lambda<1$. As illustrated in Table I the maximum $\lambda$ for which problem (48) is feasible is $\lambda=0.999$ which is $99.9\%$ of the stability value of $1$.

Example $3$

Finally, the Euler-Bernoulli beam model in (5) is stable for $\lambda<1$ (see [10, Exercise 8.3]). From Table I we observe that the maximum $\lambda$ for which problem (48) is feasible is $\lambda=0.999$ which is $99.9\%$ of the stability value of $1$. Note that owing to memory constraints, we could only solve the problem  (48) for a maximum degree of $d=6$.

Note that our method successfully establishes the stability of the example PDEs within $99.9\%$ of the approximated/analytic stability margin for $\lambda$. Finally, the time taken to search for variables which solve the problem (48) is provided in Table II.

VI CONCLUSION

We presented a method to verify stability for a large class of 1-D PDEs with polynomial data. The presented methodology relies on using Lyapunov’s method to reduce the problem of stability to the verification of integral inequalities. An application of the fundamental theorem of calculus and Green’s theorem allows us to formulate a polynomial problem for verifying the integral inequalities. Using the properties of SOS polynomials allowed us to solve the polynomial problem in a computationally efficient manner by casting the problem as an SDP. The numerical experiments indicate that the method can predict the stability of the systems considered up to a high degree of accuracy. We would like to extend this method to consider an even larger class of PDEs, for example, by including Partial (Integro)-Differential Equations (P(I)DEs) and boundary feedback. Furthermore, we would like to formulate this theory for general PDEs, i.e., PDEs not constrained to have polynomial data. Eventually, we would like to extend this framework to in-domain/boundary controller synthesis for PDEs.

In the Appendix we show how $\mcz{V}_{d}$ is cast as (II-3) and we provide the proofs of the results used in this paper.

The following Lemma is used throughout the exposition.

Lemma .1

For any bivariate matrices $K_{1}(x,y),K_{2}(x,y)\in\mathbb{R}^{\beta(\alpha+1)\times\beta(\alpha+1)}$, the following identity holds

[TABLE]

Proof:

We have

[TABLE]

where we first transposed the integrand, followed by a change of order of integration and finally switched between the variables $x$ and $y$. Thus

[TABLE]

Following the same steps for

[TABLE]

then completes the proof. ∎

We now show how (II-3) is formulated. We begin by writing the integral expression in (6) as111For brevity we have dropped the temporal dependency of $w$.

[TABLE]

where the self-adjoint operator $\Xi$ on $\mathcal{L}_{2}\left([0,1];\mathbb{R}^{\beta}\right)$ is defined as

[TABLE]

Since the operator $\Xi$ is self-adjoint, we may use (2) to obtain

[TABLE]

Therefore,

[TABLE]

where,

[TABLE]

The term $\Phi_{1}$ may be written as

[TABLE]

where we have transposed the two double integrals. Then, applying Lemma .1 to the double integrals and writing the single integral kernel in a symmetric form produces

[TABLE]

The term $\Phi_{2}$ may be written as

[TABLE]

Substituting (51) and (52) into (50) gives (II-3).

We now provide the proofs of Lemmas 2 and 3.

Proof:

Consider the vector field

[TABLE]

Then, by Green’s theorem

[TABLE]

where $\partial\overline{\Omega}$ denotes the boundary of the domain $\overline{\Omega}$. Therefore, we obtain

[TABLE]

Using the definitions of the projection matrices described in Section I-C, we obtain

[TABLE]

Similarly, we also obtain the following identity

[TABLE]

Then, applying Lemma .1 with $K_{1}=2\bar{H}$ and $K_{2}=0$ produces

[TABLE]

Substituting (54) and (55) into (53) completes the proof. ∎

Proof:

Since for all $w\in\mathcal{B}$,

[TABLE]

we get for all $w\in\mathcal{B}$

[TABLE]

Then, writing the kernel in symmetric form completes the proof. ∎

Finally, we provide the following result which we will use in the proof of Theorem 2 and Corollary 1.

Lemma .2

For any $v\in\mathcal{L}_{2}\left([0,1];\mathbb{R}^{\beta}\right)$ and polynomial matrices $F_{1},F_{2},G_{1},G_{2}\in\mathcal{R}^{\beta\times\beta}[x,y]$, the following identity holds

[TABLE]

where

[TABLE]

The same result holds for any $v\in\mathcal{H}^{\alpha}\left([0,1];\mathbb{R}^{\beta}\right)$, $F_{1},F_{2},G_{1},G_{2}\in\mathcal{R}^{\beta(\alpha+1)\times\beta(\alpha+1)}[x,y]$ and $v(x)$ replaced by $v_{\alpha}(x)$ in (56).

Proof:

We begin by observing that the left hand side of (56) may be written as

[TABLE]

where the linear bounded operators on $\mathcal{L}_{2}\left([0,1];\mathbb{R}^{\beta}\right)$ are defined as

[TABLE]

and where the adjoint operator is given by

[TABLE]

Now, using the fact that

[TABLE]

we get

[TABLE]

where

[TABLE]

In each of the $\Theta_{i}(x)$ we change the order of integration and switch between the variables $y$ and $z$ to obtain

[TABLE]

Substituting into (58) and consequently in (56) we get

[TABLE]

where

[TABLE]

Then, applying Lemma .1 to (59) completes the proof. ∎