A Partial Integral Equation (PIE) Representation of Coupled Linear PDEs and Scalable Stability Analysis using LMIs

TL;DR

This paper introduces a novel PIE representation for coupled linear PDEs enabling stability analysis via convex optimization, avoiding discretization, and demonstrating scalability up to 40 PDEs with practical implementation in MATLAB.

Contribution

It provides a standardized PIE representation for coupled linear PDEs and develops convex optimization-based stability analysis methods that are scalable and implementable in PIETOOLS.

Findings

PIE representation is equivalent to well-posed coupled PDEs.

Convex optimization verifies stability without discretization.

Algorithm scales to systems with up to 40 PDEs.

Abstract

We present a new Partial Integral Equation (PIE) representation of Partial Differential Equations (PDEs) in which it is possible to use convex optimization to perform stability analysis with little or no conservatism. The first result gives a standardized representation for coupled linear PDEs in a single spatial variable and shows that any such PDE, suitably well-posed, admits an equivalent PIE representation, defined by the given conversion formulae. This leads to a new prima facie representation of the dynamics without the implicit constraints on system state imposed by boundary conditions. The second result is to show that for systems in this PIE representation, convex optimization may be used to verify stability without discretization. The resulting algorithms are implemented in the Matlab toolbox PIETOOLS, tested on several illustrative examples, compared with previous results, and the code has been posted on Code Ocean. Scalability testing indicates the algorithm can analyze systems of up to 40 coupled PDEs on a desktop computer.

Tables

Table 1. Table 1: Number of PDEs vs. Computation Time for Stability Test

$n$	$1$	$5$	$10$	$20$	$30$	$40$
sec	.504	1.907	71.63	2706	23920	103700

Equations

$$\mathcal{T}\dot{\mathbf{x}}(t)$$

$$\mathbf{x}(0)$$

$$\left({\mathcal{P}}_{\{N_0,N_1,N_2\}}\mathbf{x}\right)(s):=N_0(s)x(s)ds$$

$$+{\int }_a^s{N}_1(s,\theta )x(\theta )d\theta +{\int }_s^b{N}_2(s,\theta )x(\theta )d\theta .\text{vspace}\text{eqnspace}$$

$$\left[\begin{array}{c}{\dot{x}}_0(t,s)\\ {\dot{x}}_1(t,s)\\ {\dot{x}}_2(t,s)\end{array}\right]$$

$$+A_2(s){\left[\begin{array}{c}{x}_2(t,s)\end{array}\right]}_{ss}$$

$$X=\left\{\left[\begin{array}{c}{x}_0\\ x_1\\ x_2\end{array}\right]\in L_2^{n_0}\times H_1^{n_1}\times H_2^{n_2}:B\left[\begin{array}{c}{x}_1(a)\\ x_1(b)\\ x_2(a)\\ x_2(b)\\ x_{2s}(a)\\ x_{2s}(b)\end{array}\right]=0\right\}.\text{vspace}\text{eqnspace}$$

$$\ddot{u}(t,s)$$

$$u(t,0)$$

$$\left[\begin{array}{c}{\dot{u}}_1(t,s)\\ {\dot{u}}_2(t,s)\end{array}\right]=\underset{A_0}{\underbrace{\left[\begin{array}{cc}-2a& -a^2\\ 1& 0\end{array}\right]}}\left[\begin{array}{c}{u}_1(t,s)\\ u_2(t,s)\end{array}\right]+\underset{A_2}{\underbrace{\left[\begin{array}{c}1\\ 0\end{array}\right]}}{u}_{2ss}(t,s)\text{vspace}\text{eqnspace}$$

$$X=\left\{\left[\begin{array}{c}{u}_1\\ u_2\end{array}\right]\in H_1^1\times H_2^1:\underset{B}{\underbrace{\left[\begin{array}{cccccc}0& 0& 1& 0& 0& 0\\ 1& 0& 0& 0& 0& 0\\ 0& k& 0& 0& 0& 1\end{array}\right]}}\left[\begin{array}{c}{u}_1(0)\\ u_1(1)\\ u_2(0)\\ u_2(1)\\ u_{2s}(0)\\ u_{2s}(1)\end{array}\right]=0\right\}.\text{vspace}\text{eqnspace}$$

$${\mathbf{x}}_f=\left[\begin{array}{c}{x}_0\\ x_{1s}\\ x_{2ss}\end{array}\right]=\left[\begin{array}{ccc}I& & \\ & {\partial }_s& \\ & & {\partial }_s^2\end{array}\right]\mathbf{x}\in L_2^{n_0+n_1+n_2}.\text{vspace}\text{eqnspace}$$

$${\begin{array}{c}\text{Find}\end{array}}_{\mathcal{P}={\mathcal{P}}_{\{N_i\}}}:\mathcal{P}\ge ϵI,$$

$${\mathcal{T}}^{\ast }\mathcal{P}\mathcal{A}+{\mathcal{A}}^{\ast }\mathcal{P}\mathcal{T}\le -\delta {\mathcal{T}}^{\ast }\mathcal{T}.$$

$$\dot{\mathbf{x}}(t,s)=A_0(s)\mathbf{x}(t,s)+A_1(s){\mathbf{x}}_s(t,s)+A_2(s){\mathbf{x}}_{ss}(t,s)\text{vspace}\text{eqnspace}$$

$$V(\mathbf{x})={⟨\mathbf{x},{\mathcal{P}}_{\{M,0,0\}}\mathbf{x}⟩}_{L_2}={\int }_a^b\mathbf{x}(s{)}^{T}M(s)\mathbf{x}(s)ds.\text{vspace}\text{eqnspace}$$

$$\dot{V}(\mathbf{x})=\underset{a}{\overset{b}{\int }}{\left[\begin{array}{c}\mathbf{x}(s)\\ {\mathbf{x}}_s(s)\\ {\mathbf{x}}_{ss}(s)\end{array}\right]}^{T}D(s)\left[\begin{array}{c}\mathbf{x}(s)\\ {\mathbf{x}}_s(s)\\ {\mathbf{x}}_{ss}(s)\end{array}\right]ds$$

$$D(s):=\left[\begin{array}{ccc}{A}_0(s{)}^{T}M(s)+M(s)A_0(s)& M(s)A_1(s)& M(s)A_2(s)\\ A_1(s{)}^{T}M(s)& 0& 0\\ A_2(s{)}^{T}M(s)& 0& 0\end{array}\right].\text{vspace}\text{eqnspace}$$

$$\dot{\mathbf{u}}(t,s)=\mathbf{u}(t,s),\mathbf{u}(t,0)=w_1(t),{\mathbf{u}}_s(t,0)=w_2(t).\text{vspace}\text{eqnspace}$$

$$\dot{\mathbf{u}}(t,s)=s{w}_1(t)+w_2(t)+{\int }_0^s(s-\eta ){\mathbf{u}}_{\eta \eta }(t,\eta )d\eta \text{vspace}\text{eqnspace}$$

$$\mathbf{I}(s)=\{\begin{array}{ll}1,& \text{mbox}ifs>0\\ 0,& \text{mbox}otherwise.\end{array}\text{vspace}\text{eqnspace}\text{vspace}-3mm$$

$$\left[\begin{array}{c}{\dot{x}}_0(t,s)\\ {\dot{x}}_1(t,s)\\ {\dot{x}}_2(t,s)\end{array}\right]$$

$$+A_2(s)\left[\begin{array}{c}{\partial }_s^2{x}_2(t,s)\end{array}\right].$$

$$X:=\left\{\left[\begin{array}{c}{x}_0\\ x_1\\ x_2\end{array}\right]\in L_2^{n_0}\times H_1^{n_1}\times H_2^{n_2}:B\left[\begin{array}{c}{x}_1(a)\\ x_1(b)\\ x_2(a)\\ x_2(b)\\ x_{2s}(a)\\ x_{2s}(b)\end{array}\right]=0\right\}\text{vspace}\text{eqnspace}$$

$$B\left[\begin{array}{ccc}{I}_{n_1}& 0& 0\\ I_{n_1}& 0& 0\\ 0& I_{n_2}& 0\\ 0& I_{n_2}& (b-a)I_{n_2}\\ 0& 0& I_{n_2}\\ 0& 0& I_{n_2}\end{array}\right]\text{is invertible.}\text{vspace}\text{eqnspace}$$

$$T^{\perp }:=\left[\begin{array}{cccccc}{I}_{n_1}& -I_{n_1}& 0& 0& 0& 0\\ 0& 0& I_{n_2}& -I_{n_2}& I_{n_2}(b-a)& 0\\ 0& 0& 0& 0& I_{n_2}& -I_{n_2}\end{array}\right].\text{vspace}\text{eqnspace}$$

$${\int }_a^b{x}_{1s}(s)ds=0.\text{vspace}\text{eqnspace}$$

$${\int }_a^b{x}_{2ss}(s)ds=0,\text{vspace}\text{eqnspace}$$

$${\int }_a^b{\int }_a^s{x}_{2\eta \eta }(\eta )d\eta ds=0.\text{vspace}\text{eqnspace}$$

$$\ddot{u}(t,s)=-c{u}_{ssss}(t,s),\text{where}$$

$$u(0)=u_s(0)=u_{ss}(L)=u_{sss}(L)=0.$$

$${\dot{u}}_1$$

$${\dot{u}}_2$$

$$\dot{\mathbf{x}}(t)=\underset{A_2}{\underbrace{\left[\begin{array}{cc}0& -c\\ 1& 0\end{array}\right]}}{\mathbf{x}}_{ss}(t)\text{vspace}\text{eqnspace}$$

$$u_{ss}(L)=u_2(L)=0\text{and}{u}_{sss}(L)=u_{2s}(L)=0.\text{vspace}\text{eqnspace}$$

Full text

,

††thanks: This work was supported by grants NSF # CNS-1739990, CMMI-1935453 and ONR #N000014-17-1-2117.

1 Introduction

Partial Differential Equations (PDEs) are used to model systems where the state depends continuously on both time and secondary independent variables. Common examples of such secondary dependence include space; as in flexible structures (Bernoulli-Euler beams) and fluid flow (Navier-Stokes); or maturation, as in cell populations and predator-prey dynamics.

The most common method for computational analysis of PDEs is to project the infinite-dimensional state onto a finite-dimensional vector space using, e.g.  [1, 2, 3] and to use the extensive literature on control of Ordinary Differential Equations (ODEs) to test stability of, and design controllers for, the resulting finite-dimensional system. However, such discretization approaches are prone to instability (e.g. in the case of hyperbolic balance laws [4]), numerical ill-conditioning, and large-dimensional state-spaces. Furthermore, representation of PDEs using ODEs inevitably neglects higher-order modes, modes which can be inadvertently excited by feedback control via the well-known spillover effect [5].

Work on computational methods for analysis and control of PDEs which do not rely on discretization has been more limited. Perhaps the most well-known computational method for stabilization of PDEs without discretization is the backstepping approach to controller synthesis [6, 7, 8]. This approach is not optimization-based, however, and not typically used for stability analysis (An exception being [9]). Recently, there has been some work on the use of Linear Matrix Inequalities (LMIs) to find Lyapunov functions for linear and nonlinear PDEs - See [10, 11, 12, 13]. However, because most of these works focus on the nonlinear case, the Lyapunov functions proposed therein are relatively simple and the resulting stability conditions conservative. An extension of the IQC framework to PDEs can be found in [14].

Numerous analytic (non-computational) methods have been proposed over the years for analysis of PDEs, including the well-developed literature on Semigroup theory [15, 16, 17, 4, 18, 19] and the literature on Port-Hamiltonian systems [20] for selecting boundary inputs. However, these methods are typically ad-hoc - relying on the expertise of the user to propose and test energy metrics.

Recently, Sum-of-Squares (SOS) has been used for analysis and control of PDEs and examples can be found in [21, 22, 23, 24] and [25, 26, 27, 28]. While these SOS-based works are relatively accurate, they: 1) Are mostly limited to scalar PDEs; 2) Suffer from high computational complexity; 3) Are mostly ad-hoc, requiring significant effort to extend the results to new PDEs. For example, these methods have never been able to analyze stability of beam or wave equations. The source of the difficulty in using LMIs and SOS for stability analysis of PDEs is that the solution of a PDE is required to satisfy three sets of constraints: the differential equation; the boundary conditions; and continuity constraints. This is in contrast to ODEs, which are defined by bounded linear operators (matrices) and solutions to which need only satisfy a single differential equation.

The Goal of the Paper is to create, for PDEs, an equivalent of the LMI framework developed for ODEs. Historically, PDEs (as old as Newton) are defined by two conflicting sets of equations: the PDE itself, which moves the state; and the Boundary Conditions (BCs), which implicitly constrain the motion of the state. We want to unify these conflicting constraints into a new state-space representation of PDEs, defined by an algebra of bounded linear operators, and which directly incorporates: the PDE, the BCs, and the continuity constraints – thereby eliminating issues of well-posedness and obviating the need to account for implicit constraints on the state.

This approach is fundamentally different than previous work using SOS or LMI-based methods. These previous efforts used SOS or positive matrices to propose candidate Lyapunov functions and then attempted to integrate the effect of boundary conditions into the derivative using, e.g. integration by parts. By contrast, our approach is to integrate the effect of boundary conditions into the dynamics - thereby obviating the need to account for them in the stability analysis. As a result, our algorithms have no obvious source of conservatism and scale to systems of up to 40 coupled PDEs.

Approach: In this paper, we propose the Partial Integral Equation (PIE) representation of PDEs. PIEs are infinite-dimensional state-space systems of the form

[TABLE]

where the state, $\mathbf{x}(t)$ is in $L_{2}^{n}[a,b]$, and the system parameters ($\mathcal{T},\mathcal{A}$) are Partial Integral (3-PI) operators. 3-PI refers to the 3 matrix-valued parameters, $\{N_{0},N_{1},N_{2}\}$ which define every such operator $\mathcal{P}_{\{N_{i}\}}:L_{2}[a,b]\rightarrow L_{2}[a,b]$ as

[TABLE]

As shown in Section 4, all 3-PI operators are $L_{2}$-bounded and together, they form an algebra (closed under addition, composition, scalar multiplication). Because they are algebraic, 3-PI operators inherit many of the properties of matrices and there is now a Matlab toolbox, PIETOOLS (using SOSTOOLS as a model), which allows for manipulation of 3-PI operators using matrix syntax and which can solve Linear PI Inequality (LPI) constrained optimization problems using the YALMIP syntax for LMIs.

Converting PDE state to PIE state: The first contribution of the paper (extending the results in [29]) is to show that the solutions to a broad class of PDEs can be represented using PIEs. For this result, we propose an alternative state-space. To explain this change of state, we consider the following standardized representation of coupled linear PDEs in a single spatial variable, presented in Section 3.

[TABLE]

with associated state-space

[TABLE]

Most 1D PDEs can be formulated using this standardized representation.

For example, if we consider the damped wave equation

[TABLE]

Then, setting $u_{1}=\dot{u}$ and $u_{2}=u$, we have

[TABLE]

with BCs $u_{2}(0)=0$, $u_{1}(0)=0$, and $u_{2s}(1)=-ku_{1}(1)$ so that

[TABLE]

The Fundamental State-Space A reasonable definition of state is the minimal amount of information needed to forward-propagate the solution. By this measure, and referring to the example, defining the state of a PDE as $u(t)\in H_{1}\times H_{2}$ is not minimal, as this function contains redundant information regarding the boundary values. We propose, then, that for a PDE defined by $A_{i}$ and $X$, the correct definition of state is the so-called fundamental state, where for $\mathbf{x}\in X$, we define

[TABLE]

In Section 6, we show that if the PDE is suitably well-posed, there exists a unitary 3-PI operator $\mathcal{T}:L_{2}^{n_{0}+n_{1}+n_{2}}\rightarrow X$ such that $\mathbf{x}=\mathcal{T}\mathbf{x}_{f}$.

Equivalence of PIE and PDE: In Section 7, equipped with the unitary operator, $\mathcal{T}$, we propose a 3-PI operator $\mathcal{A}$ such that for any solution $\mathbf{x}_{f}(t)\in L_{2}$ of the PIE in Eqn. (LABEL:eqn:PIE_0), $\mathbf{x}(t)=\mathcal{T}\mathbf{x}_{f}(t)$ satisfies the PDE defined by $A_{i},X$ and that, conversely, for any solution of the PDE, $\mathbf{x}$, the fundamental state $\mathbf{x}_{f}(t)=\text{diag}(I,\partial_{s},\partial_{s}^{2})\mathbf{x}(t)$ satisfies the PIE. We further show that exponential stability of the PIE and PDE are equivalent.

An Linear PI Inequality (LPI) for Stability: Aside from have a minimal state-space, the advantage of the PIE representation of a PDE is computational. Recall that our goal is to create a framework akin to LMIs which can be used to study PDEs. Because PIEs do not have BCs, all information needed to define the solution is contained in the fundamental state, $\mathbf{x}_{f}\in L_{2}$. Therefore, unlike PDEs, where the effect of the BCs needs to be “brought in” using integration by parts, Poincare inequalities, etc., PIEs are a prima facie representation of the dynamics. This allows us to pose the stability test for PIEs as a convex optimization problem of the following form:

[TABLE]

This LPI, then, is a straightforward generalization of the Lyapunov LMI: $P>0$, $A^{T}P+PA<0$. Furthermore, again motivated by the efficient Matlab parser YALMIP [30], we have recently developed the Matlab toolbox PIETOOLS allows us to manipulate 3-PI operators using matrix syntax and solve Linear PI Inequality (LPI) constrained optimization problems. Thus, if $\{\mathcal{T},\mathcal{A}\}$ are as defined in Section 7, then our stability test reduces to 3 lines of Matlab code

[TABLE]

where the functions sos_posopvar and sos_opineq are defined in Section 9.3 and in [31]. We emphasize that this code applies to any PDE in standardized format and since there is no need to bring in the boundary conditions, there is no obvious source of conservatism. Specifically, in Section 11, we apply the algorithm to beam and wave equations for which there are no previous LMI-based stability conditions. Furthermore, the lack of conservatism is verified in Section 10 by comparing against known stability margins taken from the literature.

Finally, computational complexity depends on the degree of the polynomial parameters in the 3-PI variable $\mathcal{P}$ (corresponding to the complexity of the candidate Lyapunov funtion). However, most problems only require very simple Lyapunov functions. In this case, the scalability of the method is comparable to the complexity of discretization-based analysis. Specifically, if we choose the polynomial parameters to have degree 2, then the algorithm can analyze stability of more than 40 coupled PDEs on a desktop computer.

In the following two illustrations, we attempt to further motivate our BC free PIE representation by illustrating how BCs affect the dynamics and complicate stability analysis in the original PDE state.

1.1 BCs Complicate Stability Analysis

The goal of the paper is to find a representation of PDEs which can be interpreted as we would interpret an ODE - without implicit constraints on the solution imposed by BCs or continuity. To motivate this goal, and to illustrate that BCs cannot be ignored in the PDE framework, let us illustrate what happens when we treat a PDE like an ODE - i.e. without considering continuity constraints or boundary conditions. Obviously, in a more classical analysis, we may account for such boundary conditions using, e.g. integration-by-parts. However, remember that our goal in this paper is to eliminate the need for such ad hoc manipulations - as is possible in the PIE representation. To this end, suppose we are given a PDE, parameterized by $A_{0},A_{1},A_{2}$ as follows.

[TABLE]

An obvious class of candidate Lyapunov functions for this system is parameterized by $M$ as

[TABLE]

As would be the case for an ODE, $V(\mathbf{x})\geq\epsilon\left\lVert\mathbf{x}\right\rVert^{2}$ if $M(s)\geq\epsilon I$ for all $s$ and some $\epsilon>0$ - a constraint which is easy to enforce using, e.g. SOS. However, if we now take the derivative of this candidate Lyapunov function we obtain

[TABLE]

Now, if we were to treat the system like an ODE (without considering the BCs and continuity constraints), we would constrain $D(s)\leq 0$ and this would imply stability. Unfortunately, however, $D(s)\not\leq 0$ for ANY choice of $M,A_{1},A_{2}\neq 0$. The problem, of course, is that the differentiation operator branches $\mathbf{x}$ into $\mathbf{x}_{s}$ and $\mathbf{x}_{ss}$, neither of which are independent of $\mathbf{x}$. Moreover, the information which determines the relationship between $\mathbf{x}$, $\mathbf{x}_{s}$ and $\mathbf{x}_{ss}$ is not embedded in the differential equation. Rather, this information is implicit in the BCs and continuity constraints. At this point, of course, the user would attempt to “bring in” the boundary and continuity properties to obtain a new stability condition using, e.g. integration by parts or Stokes Theorem [28]. However, such secondary steps require significant expertise, are ad-hoc, and may result in conservatism. In the analysis of ODEs, by contrast, no such secondary analysis is necessary. The goal of this paper, then, is to find a prima facie representation of the PDE with no implicit constraints on continuity and BCs, and wherein we may develop a computational framework which mirrors that of LMIs for ODEs.

1.2 BCs Significantly alter the Dynamics

Having observed that BCs complicate the analysis of PDEs, let us now examine how the BCs can alter the nature of a distributed parameter system. Specifically, in this subsection, we consider an example of how the incorporation of continuity and BCs can take a system (wherein the dynamics do not explicitly depend on the partial derivatives), and transform it into a system with explicit dependence on the partial derivatives. Consider the following non-partial-differential, yet distributed-parameter system.

[TABLE]

To allow for the specified BCs, we restrict continuity of $\mathbf{u}$ as $\mathbf{u}(t)\in H_{2}^{1}$. Note that this system is not, prima facie, a PDE and that the dynamics are identical at every point in the domain. However, if we now combine the fundamental theorem of calculus with integration by parts, we obtain a very different, yet equivalent, set of dynamics - dynamics with explicit dependence on the partial derivatives.

[TABLE]

This formulation of the same system directly incorporates continuity and BCs into the dynamics - which are now expressed using the partial derivative $\mathbf{u}_{ss}$. In addition, while the original formulation was spatially decoupled, with $w_{1},w_{2}$ only acting at the boundary, the new formulation shows that the exogenous inputs $w_{1},w_{2}$ are felt instantaneously at every point in the domain. This observation demonstrates that the continuity properties of the solution (combined with BCs) fundamentally alter the dynamics.

The first goal of this paper, then, is: 1) to show that a broad class of PDEs can be reformulated in a way which specifies precisely how the BCs affect the dynamics and 2) to provide universal formulae for constructing such a representation.

2 Notation

We define $L_{2}[a,b]^{n}$ to be space of $\mathbb{R}^{n}$-valued Lesbegue integrable functions defined on $[a,b]$ and equipped with the standard inner product. We use $W^{k,p}[a,b]^{n}$ to denote the Sobolev subspace of $L_{p}[a,b]^{n}$ defined as $\{u\in L_{p}[a,b]^{n}\,:\,\frac{\partial^{q}}{\partial x^{q}}u\in L_{p}\text{ for all }q\leq k\}$. $H_{k}[a,b]:=W^{k,2}[a,b]$ and $H_{k}^{n}[a,b]=\prod_{i=1}^{n}H_{k}[a,b]$. For efficiency, we typically omit the domain, so that, e.g. $H_{k}^{n}:=H_{k}^{n}[a,b]$ unless otherwise stated. $I_{n}\in\mathbb{R}^{n\times n}$ and $0_{n_{1}\times n_{2}}\in\mathbb{R}^{n_{1}\times n_{2}}$ are used to denote the identity and zero matrices and the subscripts are omitted when the dimension of the matrices is clear from context. $\mathbf{I}$ denotes the indicator function $\mathbf{I}:\mathbb{R}\rightarrow\{0,1\}$, defined as

[TABLE]

3 A Standardized PDE Representation

The two primary contributions of this paper are: a formula for conversion of PDEs to PIEs; and an LPI framework for Lyapunov stability analysis of PIEs (Section 8). The significance of the latter result clearly depends on the set of PDEs which can be converted to PIEs. In this section, we propose a standardized framework for representation of PDEs. In Section 7, we will show that for any such PDE, there exists a PIE for which a solution to the PIE yields a solution to the PDE and vice-versa. The class of PDEs considered here is not exhaustive, however. That is, there exist PDEs not listed here which may be converted to PIEs. Furthermore, there exist PIEs which do not have a coupled PDE representation.

We consider coupled PDEs with infinitesimal generator of the following form

[TABLE]

and with domain

[TABLE]

where

[TABLE]

Specifically, given $\mathbf{x}_{0}\in X$, we say that $\mathbf{x}$ satisfies the PDE defined by $\{A_{i},X\}$ if $\mathbf{x}$ is Frechét differentiable, $\mathbf{x}(0)=\mathbf{x}_{0}$, $\mathbf{x}(t)\in X$ and Equation (3) is satisfied for all $t\geq 0$.

3.1 A Guide to Partition of States

The partition of states in Equation (3) is not overly restrictive and there is no special structure to this generator. The partition is purely organizational and defined by the domain restriction $\mathbf{x}\in L_{2}^{n_{0}}\times H_{1}^{n_{1}}\times H_{2}^{n_{2}}$. States with no restrictions on continuity are assigned to be an element of the vector $x_{0}(t)\in L_{2}[a,b]^{n_{0}}$. If a state has no continuity properties, then these states cannot be differentiated and it is not possible to assign boundary conditions, as the limits $x_{0}(a),x_{0}(b)$ do not exist. States which are continuous, but not continuously differentiable are assigned to $x_{1}(t)\in H_{1}^{n_{1}}$. The continuity property of these states allow for boundary conditions, as $x_{1}(a),x_{1}(b)$ exist. However, since $\partial_{s}x_{1}\in L_{2}[a,b]^{n_{1}}$ is not continuous, we cannot assign boundary conditions which involve $x_{1s}(a)$ or $x_{1s}(b)$, as these limits do not exist. Finally, states which are required to be continuously differentiable are assigned to $x_{2}(t)\in H_{2}^{n_{2}}$ and admit boundary conditions involving $x_{2s}(a)$ or $x_{2s}(b)$ and second-order spatial derivatives, $\partial^{2}_{s}x_{2}$. This standardized representation specifically excludes states in $H^{n}_{k}$ where $k>2$. Although the results of the paper can be extended to such systems, such an extension is not considered here.

3.2 A Guide to Boundary Conditions

In this subsection, we propose restrictions on the matrix $B$ which are equivalent to Equation (5). Specifically, the row rank of $B$ must be $n_{1}+2n_{2}$ and $B$ contains no boundary conditions of a given prohibited form. Note that the rank condition on $B$ is not overly restrictive as, to the best of our knowledge, it is a necessary condition for existence of a unique solution for any PDE in standardized form.

3.2.1 Prohibited Boundary Conditions

A necessary and sufficient condition for $B$ to satisfy Equation (5) is for $B\in\mathbb{R}^{(n_{1}+2n_{2})\times(2n_{1}+4n_{2})}$ to have row rank $n_{1}+2n_{2}$ and to define no boundary conditions consisting of a linear combination of $x_{1}(a)-x_{1}(b)=0$, $x_{2}(a)+(b-a)x_{2s}(a)-x_{2}(b)=0$, or $x_{2s}(a)-x_{2s}(b)=0$.

Lemma 1**.**

Suppose $B\in\mathbb{R}^{(n_{1}+2n_{2})\times(2n_{1}+4n_{2})}$. Define

[TABLE]

Equation (5) is satisfied if and only if $B\in\mathbb{R}^{(n_{1}+2n_{2})\times(2n_{1}+4n_{2})}$, has row rank $n_{1}+2n_{2}$ and the row space of $B$ and the row space of $T^{\perp}$ has trivial intersection.

{pf}

Let the $T$ matrix from Eqn. (5) be as defined in Eqn. (9). Now suppose $BT$ is invertible. Since $T\in\mathbb{R}^{(2n_{1}+4n_{2})\times(n_{1}+2n_{2})}$, we require $B\in\mathbb{R}^{(n_{1}+2n_{2})\times(2n_{1}+4n_{2})}$ in order for $BT$ to exist and be square. Since $BT\in\mathbb{R}^{(n_{1}+2n_{2})\times(n_{1}+2n_{2})}$, $B$ must also have row rank $n_{1}+2n_{2}$. Now, since $T$ has column rank $n_{1}+2n_{2}$, its row rank is also $n_{1}+2n_{2}$. Now $T^{\perp}$ has row rank $n_{1}+2n_{2}$ and $T^{\perp}T=0$. Thus the row space of $T^{\perp}$ lies in $Im(T)^{\perp}$ and since $Im(T)^{\perp}$ is of dimension $n_{1}+2n_{2}$, the row space of $T^{\perp}$ is $Im(T)^{\perp}$. Therefore $x^{T}BT=0$ for some $x\neq 0$, if and only if the intersection of the row space of $B$ and that of $T^{\perp}$ is non-trivial. This establishes necessity. For sufficiency, we assume $B\in\mathbb{R}^{(n_{1}+2n_{2})\times(2n_{1}+4n_{2})}$ and has row rank $n_{1}+2n_{2}$. Furthermore, as shown, the row space of $B$ has trivial intersection with $Im(T)^{\perp}$. Again, we have that $x^{T}BT=0$ implies $x=0$, from which we conclude invertibility. $\square$

Note 3.1**.**

The restriction on prohibited boundary condition is subtle. For example, $x_{2}(a)=x_{2}(b)$ is permitted, except if combined with $x_{2s}(a)=0$ (combining with $x_{2s}(b)=0$ is still OK). Meanwhile, $x_{1}(a)=x_{1}(b)$ and $x_{2s}(a)=x_{2s}(b)$ are never OK. Additionally, $x_{2s}(a)=0$ is OK, unless combined with $x_{2}(a)=x_{2}(b)$ or $x_{2s}(b)=0$. Of course, the most reliable way to check if certain boundary conditions are permitted is to simply construct $B$ and check the rank of $BT$. The PIETOOLS implementation described in Section 9.3 will do this automatically and generate an error if $BT$ is not invertible.

3.2.2 A Note on Necessity of Equation (5)

Boundary conditions of the form $x_{1}(a)=x_{1}(b)$ are periodic and imply

[TABLE]

Likewise $x_{2s}(a)=x_{2s}(b)$ implies

[TABLE]

and $x_{2}(a)+(b-a)x_{2s}(a)=x_{2}(b)$ implies

[TABLE]

In this way, the prohibited BCs represent integral constraints on the respective PIE (fundamental) states, $x_{1s}\in L_{2}$ and $x_{2ss}\in L_{2}$, meaning these PIE states are not minimal (dynamics expressed using these states have implicit constraints). One option in these cases may be to redefine the PIE states modulo an integral constraint, however this extension is left for future work.

3.3 Euler-Bernoulli Beam Example

In order to better understand how to write a PDE in the standardized PDE form of Eqn. (3), let us consider the cantilevered Euler-Bernoulli beam:

[TABLE]

We wish to construct a standardized PDE representation of this classic diffusive model. Following the approach in, e.g. [20] (from which we also get the Timoshenko beam model in Section 11), we first introduce the augmented state $u_{1}=\dot{u}$ - a choice which leads to “natural” BCs for which the system is well-posed [19]. This choice also eliminates the second order time-derivative, $\ddot{u}$. Since $u\in H_{4}$, we eliminate the fourth-order spatial derivative by creating the augmented state $u_{2}=u_{ss}$. Taking the time-derivative of these states, we obtain

[TABLE]

These equations are now in the standardized form

[TABLE]

where $A_{0}=A_{1}=0$, $n_{2}=2$, and $n_{0}=n_{1}=0$. We now examine the boundary conditions using these states:

[TABLE]

These boundary conditions are insufficient, as the resulting rank is 2. However, we may differentiate boundary conditions in time to obtain new boundary conditions

[TABLE]

We now have 4 boundary conditions, which we use to construct the $B$ matrix as

[TABLE]

The $B$ matrix is now of rank $4=n_{1}+2n_{2}$ and satisfies Equation (5).

Now, if $u$ satisfies the E-B beam equation for some initial condition, we have that $u_{1}=\dot{u},u_{2}=u_{ss}$ satisfy the standardized PDE model. However, conversely, if $u_{1},u_{2}$ satisfy the standardized PDE for some initial condition, then in order to construct a solution to the original PDE, we must integrate $u_{1}$ in time. However, this requires knowledge of $u(0)$. Thus we find that some information on the system solution has been lost in the standardized representation. Note, however, that we could retain this information by including a third state, $u_{3}=u$, so that $\dot{u}_{3}=u_{1}$ and then the solutions would be equivalent.

3.4 Exponential Stability of Coupled PDE Systems

In this subsection, we define two notions of exponential stability with respect to the standardized PDE representation - stability in the $X$ norm and stability in the $L_{2}$ norm. In Section 5, we will define the notion of exponential stability for PIEs. In Section 8, we will show that exponential stability of a PIE representation of a standardized PDE is equivalent to exponential stability of the original standardized PDE in the $X$ norm.

Definition 2**.**

We say the PDE defined by $\{A_{i},X\}$ is exponentially stable in $X$ if there exist constants $K,\gamma>0$ such that for any $\mathbf{x}_{0}\in X$, any solution $\mathbf{x}$ of the PDE defined by $\{A_{i},X\}$ satisfies

[TABLE]

Definition 3**.**

We say the PDE defined by $\{A_{i},X\}$ is exponentially stable in $L_{2}$ if there exist constants $K,\gamma>0$ such that for any $\mathbf{x}_{0}\in X$, any solution $\mathbf{x}$ of the PDE defined by $\{A_{i},X\}$ satisfies

[TABLE]

Note 3.2**.**

Exponential stability in $X$ implies exponential stability in $L_{2}$, since $\left\lVert\mathbf{x}\right\rVert_{L_{2}}\leq\left\lVert\mathbf{x}\right\rVert_{H_{k}}$ for any $\mathbf{x}\in H_{k}$ and $k\geq 0$. Furthermore, our definitions of exponential stability implies that all states in $\mathbf{x}$ must be exponentially decreasing in the given norm. Because not all standardized PDE representations of a given scalar high-order PDE necessarily use the same set of first-order states (See, e.g. the E-B beam example), this raises the possibility one standardized PIE representation may be exponentially stable, while another may not.

Note that in the case where the $X$ or $L_{2}$ stability definition holds with $\gamma=0$, we say that the system is Lyapunov stable or neutrally stable.

4 3-PI Operators Form an Algebra

In Section 7, we will construct a PIE representation of any PDE in the standardized form described in Section 3. PIEs, as will be defined in Section 5, have the advantage that they are parameterized by the class of 3-PI operators, which are bounded on $L_{2}$ and form an algebra. Furthermore, candidate Lyapunov functions can be parameterized using 3-PI operators. The algebraic nature of 3-PI operators significantly simplifies the problem of analysis and control of PIEs.

Formally, we say that an operator $\mathcal{P}$ is 3-PI if there exist 3 bounded matrix-valued functions $N_{0}:[a,b]\rightarrow\mathbb{R}^{n\times n}$, $N_{1}:[a,b]^{2}\rightarrow\mathbb{R}^{n\times n}$, and $N_{2}:[a,b]^{2}\rightarrow\mathbb{R}^{n\times n}$ such that

[TABLE]

where $N_{0}$ defines a multiplier operator and $N_{1},N_{2}$ define the kernel of an integral operator.

For given $N_{0},N_{1},N_{2}$, we use $\mathcal{P}_{\{N_{0},N_{1},N_{2}\}}:L_{2}^{n}\rightarrow L_{2}^{n}$ to denote the corresponding PI operator. When clear from context, we use the shorthand notation $\mathcal{P}_{\{N_{i}\}}$ to indicate $\mathcal{P}_{\{N_{0},N_{1},N_{2}\}}$.

One may interpret 3-PI operators to be an extension of matrices, wherein $N_{0}$ defines the diagonal of the matrix, $N_{1}$ contains the sub-diagonal terms, and $N_{2}$ contains the terms above the diagonal.

In the following subsections, we show that this class of bounded linear operators is closed under composition and adjoint (closure under addition and scalar multiplication follows immediately from addition and scalar multiplication of parameters). Furthermore, these results show that if we define the set of 3-PI operators with polynomial parameters $N_{0}$, $N_{1}$, and $N_{2}$, then this set forms a sub-algebra.

4.1 Composition of 3-PI operators

In this subsection, we derive an analytic formula for the composition 3-PI operators. Specifically, we have the following.

Lemma 4**.**

For any bounded functions $B_{0},N_{0}:[a,b]\rightarrow\mathbb{R}^{n\times n}$, $B_{1},B_{2},N_{1},N_{2}:[a,b]^{2}\rightarrow\mathbb{R}^{n\times n}$, we have

[TABLE]

where

[TABLE]

{pf}

To prove the theorem, we exploit the linear structure of the operator to decompose

[TABLE]

and

[TABLE]

Then

[TABLE]

We now consider each term separately, starting with the first two, which are trivial. First,

[TABLE]

where

[TABLE]

Similarly,

[TABLE]

where

[TABLE]

We now proceed to the more difficult terms. For these, recall the indicator function

[TABLE]

For the first term, we note that

[TABLE]

This allows us to change the variables of integration as follows.

[TABLE]

where

[TABLE]

Next, we use another identity

[TABLE]

to establish the following.

[TABLE]

where

[TABLE]

The remaining two identities are minor variations on the two we most recently derived.

[TABLE]

where

[TABLE]

Finally,

[TABLE]

where

[TABLE]

Combining all terms, we have

[TABLE]

The last equality holds since

[TABLE]

and

[TABLE]

$\square$

This lemma proves that the class of 3-PI operators is closed under composition.

Corollary 5**.**

Suppose that $\{B_{i}\}$ and $\{N_{i}\}$ are matrices of polynomials. Then if $\mathcal{P}_{\{R_{i}\}}=\mathcal{P}_{\{B_{i}\}}\mathcal{P}_{\{N_{i}\}}$, $\{R_{i}\}$ are matrices of polynomials.

{pf}

The algebra of polynomials is closed under multiplication and integration. Therefore, the proof follows from the expressions for $\{R_{i}\}$ given in Lemma 4.$\square$

This corollary implies that the subset of 3-PI operators with polynomial parameters is likewise closed under composition and therefore forms a subalgebra.

4.2 The Adjoint of 3-PI operators

Next, we give a formula for the adjoint of a 3-PI operator.

Lemma 6**.**

For any bounded functions $N_{0}:[a,b]\rightarrow\mathbb{R}^{n\times n}$, $N_{1},N_{2}:[a,b]^{2}\rightarrow\mathbb{R}^{n\times n}$ and any $\mathbf{x},\mathbf{y}\in L_{2}^{n}[a,b]$, we have

[TABLE]

where

[TABLE]

{pf}

Noting that

[TABLE]

we can decompose the adjoint as follows. Clearly

[TABLE]

For the other terms, we use a change of integration to yield

[TABLE]

Likewise,

[TABLE]

Combining these terms, we complete the proof. $\square$

The following Corollary follows immediately from Lemma 6.

Corollary 7**.**

Suppose that $\{N_{i}\}$ are matrices of polynomials. Then, using the adjoint with respect to $L_{2}$, if $\mathcal{P}_{\{\hat{N}_{i}\}}=\mathcal{P}_{\{N_{i}\}}^{*}$, $\{\hat{N}_{i}\}$ are matrices of polynomials.

5 Partial Integral Equations (PIEs)

In this section, we give the autonomous form of a Partial Integral Equation (PIE) and define notions of solution and exponential stability. Specifically, for given 3-PI operators $\mathcal{A},\mathcal{T}$, we say, for an initial condition, $\mathbf{x}_{0}\in L_{2}^{n}$, that $\mathbf{x}:\mathbb{R}^{+}\rightarrow L_{2}^{n}$ satisfies the PIE defined by $\{\mathcal{A},\mathcal{T}\}$ if $\mathbf{x}(0)=\mathbf{x}_{0}$, $\mathbf{x}$ is Frechét differentiable for all $t\geq 0$ and

[TABLE]

for all $t\geq 0$.

Not all PIEs are well-posed in the sense of Hadamard. However, we will show in Section 7 that if a PDE is in standardized form (satisfying Eqn. (5)), and the PIE is generated from that standardized PDE using the formulation in Section 7, then the PIE is well-posed.

5.1 Exponential Stability of PIEs

Having defined PIE’s, we now define the notion of exponential stability we will use.

Definition 8**.**

We say the PIE defined by the 3-PI operators $\{\mathcal{A},\mathcal{T}\}$ is exponentially stable if there exist constants $K$ and $\gamma>0$ such that for $\mathbf{x}(0)\in L_{2}^{n}$, any solution $\mathbf{x}$ satisfies

[TABLE]

In the case where the exponential stability definition holds with $\gamma=0$, we say the PIE is stable in the sense of Lyapunov or neutrally stable.

6 A Unitary map from $X$ to $L_{2}$

In this section, we show equivalence between the Hilbert space $L_{2}^{n_{0}+n_{1}+n_{2}}$ and the space

[TABLE]

where $B$ satisfies Equation (5) and $X$ is equipped with the inner product

[TABLE]

Specifically, in this section, we

• •

Construct a unitary map $\mathcal{T}:L_{2}^{n_{0}+n_{1}+n_{2}}\rightarrow X$ where $\mathcal{T}$ is a 3-PI operator.

• •

Show $\left\langle\cdot,\cdot\right\rangle_{X}$ is an inner product and $X$ is Hilbert with this inner product.

• •

Show that for $\mathbf{x}\in X$, the norm $\left\lVert\cdot\right\rVert_{X}$ is equivalent to the norm $\left\lVert\cdot\right\rVert_{L_{2}\times H_{1}\times H_{2}}$ where recall

[TABLE]

6.1 The Unitary Map, $\mathcal{T}$

In this subsection, we define the 3-PI operator $\mathcal{T}=\mathcal{P}_{\{G_{i}\}}$ such that if

[TABLE]

then

[TABLE]

First, we first show that a PDE state $\mathbf{x}\in H_{2}$ can be represented using the PIE state, $\partial_{s}^{2}\mathbf{x}\in L_{2}$ and a set of ‘core’ boundary conditions ($x(a),x_{s}(a)$).

Lemma 9**.**

Suppose that $\mathbf{x}\in H_{2}^{n}[a,b]$. Then for any $s\in[a,b]$,

[TABLE]

{pf}

The first two identities are the fundamental theorem of calculus. The third identity is a repeated application of the fundamental theorem of calculus, combined with a change of variables. That is, for any $s\in[a,b]$,

[TABLE]

Examining the 3rd term, where recall $\mathbf{I}(s)$ is the indicator function,

[TABLE]

which is the desired result.$\square$ As an obvious corollary, we have

[TABLE]

The implication is that any boundary value can be expressed using two other boundary identities. In the standardized PDE representation, we have a generic set of boundary conditions defined by the matrix $B$. In the following theorem, we generalize Lemma 9 in order to express the PDE state $\mathbf{x}\in X$ in terms of the PIE state, $\hat{\mathbf{x}}\in L_{2}^{n_{0}+n_{1}+n_{2}}$, and generalized BCs (which are equal to zero). This allows us to define the map $\mathcal{T}$.

Theorem 10**.**

Let $\mathcal{T}=\mathcal{P}_{\{G_{0},G_{2},G_{2}\}}$ with $G_{i}$ as defined in Equations (9). Then for any $\mathbf{x}\in X$,

[TABLE]

Furthermore, for any $\hat{\mathbf{x}},\hat{\mathbf{y}}\in L_{2}^{n_{0}+n_{1}+n_{2}}$, $\mathcal{T}\hat{\mathbf{x}},\mathcal{T}\hat{\mathbf{y}}\in X$ and $\left\langle\mathcal{T}\hat{\mathbf{x}},\mathcal{T}\hat{\mathbf{y}}\right\rangle_{X}=\left\langle\hat{\mathbf{x}},\hat{\mathbf{y}}\right\rangle_{L_{2}}$.

{pf}

Suppose $\mathbf{x}=\begin{bmatrix}x_{1}\\ x_{2}\\ x_{3}\end{bmatrix}\in X$. Let us define the PIE state, $\hat{\mathbf{x}}$, and the ‘full’ and ‘core’ sets of BCs as

[TABLE]

Clearly $\hat{\mathbf{x}}\in L_{2}^{n_{0}+n_{1}+n_{2}}$. Using Lemma 9, we can express $x_{bf}$ using $x_{bc}$ and $\hat{\mathbf{x}}$ as

[TABLE]

Likewise, we may express $\mathbf{x}$ in terms of $x_{bc}$ and $\hat{\mathbf{x}}$ as

[TABLE]

where

[TABLE]

We may now express the ‘full’ boundary conditions as

[TABLE]

Now from Equation (5), $BT$ is invertible, and hence

[TABLE]

This yields the following expression for $\mathbf{x}$.

[TABLE]

as desired.

Conversely, suppose that $\hat{\mathbf{x}}=\begin{bmatrix}\hat{x}_{0}\\ \hat{x}_{1}\\ \hat{x}_{2}\end{bmatrix}\in L_{2}^{n_{0}+n_{1}+n_{2}}$ and

[TABLE]

First, we note that using the definition of $\{G_{i}\}$, $x_{0}=\hat{x}_{0}\in L_{2}^{n_{0}}$. Next, we have

[TABLE]

Hence, $\partial_{s}Q(s,\theta)=0$ implies $\partial_{s}x_{1}=\hat{x}_{1}\in L_{2}^{n_{1}}$. Finally,

[TABLE]

Hence, since $\partial_{s}Q(s,\theta)=0$, $\partial_{s}(s-a)=1$ and $\partial_{s}^{2}(s-a)=0$, we have $\partial_{s}^{2}x_{2}=\hat{x}_{2}\in L_{2}^{n_{2}}$. We conclude that $\mathbf{x}\in L_{2}\times H_{1}\times H_{2}$. This now implies the identity

[TABLE]

Using the formulae given above, we have

[TABLE]

and

[TABLE]

which implies

[TABLE]

Hence

[TABLE]

We conclude that $BTx_{bc}=-B\mathcal{P}_{\{0,Q,Q\}}\hat{\mathbf{x}}$ and consequently

[TABLE]

from which we conclude that $\mathbf{x}\in X$.

Finally, let $\hat{\mathbf{x}},\hat{\mathbf{y}}\in L_{2}^{n_{0}+n_{1}+n_{2}}$, $\mathcal{T}\hat{\mathbf{x}},\mathcal{T}\hat{\mathbf{y}}\in X$. Then

[TABLE]

$\square$

Corollary 11**.**

Let $\mathcal{H}=\mathcal{P}_{\{G_{3},G_{4},G_{5}\}}$ with $G_{i}$ as defined in Equations (9). Then for any $\mathbf{x}\in X$,

[TABLE]

{pf}

By Theorem 10,

[TABLE]

Now, for any $\hat{\mathbf{x}}\in L_{2}^{n_{0}+n_{1}+n_{2}}$, it can be readily verified through differentiation that

[TABLE]

which completes the proof. $\square$

Corollary 12**.**

The operator $\mathcal{T}=\mathcal{P}_{\{G_{0},G_{1},G_{2}\}}:L_{2}^{n_{0}+n_{1}+n_{2}}\rightarrow X$ is unitary.

{pf}

Theorem 10 shows that for any $\mathbf{x}\in X$, there exists some $\hat{\mathbf{x}}\in L_{2}$ such that $\mathbf{x}=\mathcal{T}\hat{\mathbf{x}}$ (surjective). Furthermore, for any $\hat{\mathbf{x}},\hat{\mathbf{y}}\in L_{2}$, $\left\langle\mathcal{T}\hat{\mathbf{x}},\mathcal{T}\hat{\mathbf{y}}\right\rangle_{X}=\left\langle\hat{\mathbf{x}},\hat{\mathbf{y}}\right\rangle_{L_{2}}$, which concludes the proof. $\square$

Because $L_{2}^{n_{0}+n_{1}+n_{2}}$ is a Hilbert space and $\mathcal{T}$ is unitary, Corollary 12 implies $X$ is a Hilbert space.

6.2 Equivalence of Norms

In this subsection, we briefly show that the norms $\left\lVert\cdot\right\rVert_{X}$ and $\left\lVert\cdot\right\rVert_{L_{2}\times H_{1}\times H_{2}}$ are equivalent.

Lemma 13**.**

For any $\mathbf{x}\in X$, $\left\lVert\mathbf{x}\right\rVert_{X}\leq\left\lVert\mathbf{x}\right\rVert_{L_{2}\times H_{1}\times H_{2}}$ and there exists a constant $c>0$ such that $\left\lVert\mathbf{x}\right\rVert_{L_{2}\times H_{1}\times H_{2}}\leq c\left\lVert\mathbf{x}\right\rVert_{X}$.

{pf}

First, we note that

[TABLE]

and hence $\left\lVert\mathbf{x}\right\rVert_{X}\leq\left\lVert\mathbf{x}\right\rVert_{L_{2}\times H_{1}\times H_{2}}$. Now, since $G_{i}\in L_{\infty}[a,b]$, there exist $c_{1},c_{2}>0$ such that

[TABLE]

and

[TABLE]

Therefore, we conclude that

[TABLE]

as desired.$\square$

This result shows that for PDE systems in standardized form, stability in $\left\lVert\cdot\right\rVert_{X}$ and $\left\lVert\cdot\right\rVert_{L_{2}\times H_{1}\times H_{2}}$ are equivalent.

7 Converting PDEs to PIEs

In this section, we show that for any PDE in standardized form, there exists a PIE for which any solution to the PDE defines a solution to the PIE and any solution to the PIE defines a solution to the PDE. We further show that this result implies that exponential stability of the PIE is equivalent to exponential stability of the PDE in $X$.

7.1 Equivalence of Solutions for PDEs and PIEs

Now that we have the unitary 3-PI operator $\mathcal{T}:=\mathcal{P}_{\{G_{0},G_{1},G_{2}\}}$ where

[TABLE]

for any $\mathbf{x}\in X$, conversion of the PDE to a PIE (Eqn. (8)) is direct.

Lemma 14**.**

Given $\hat{\mathbf{x}}_{0}(t)\in L_{2}^{n_{0}+n_{1}+n_{2}}$, the function $\hat{\mathbf{x}}(t)\in L_{2}^{n_{0}+n_{1}+n_{2}}$ satisfies the PIE defined by $\{\mathcal{T},\mathcal{A}\}$ if and only if for $\mathbf{x}_{0}=\mathcal{T}\hat{\mathbf{x}}_{0}$, the function $\mathbf{x}(t)=\mathcal{T}\hat{\mathbf{x}}(t)$ satisfies the PDE defined by $\{A_{i},X\}$ where

[TABLE]

where the $G_{i}$ are as defined in Eqns. (9).

{pf}

Define $\mathcal{H}:=\mathcal{P}_{\{G_{3},G_{4},G_{5}\}}$. Suppose $\mathbf{x}$ satisfies the PDE. Define

[TABLE]

and $\hat{\mathbf{x}}(t)=\mathcal{D}_{1}\mathbf{x}(t)$. By Theorem 10 and Lemma 4 and the definition of the $G_{i}$, we have

[TABLE]

Finally, $\dot{{\mathbf{x}}}(t)=\mathcal{T}\dot{\hat{\mathbf{x}}}(t)$ and $\hat{\mathbf{x}}(0)=\mathcal{D}_{1}{\mathbf{x}}(0)=\mathcal{D}_{1}{\mathbf{x}}_{0}=\mathcal{D}_{1}\mathcal{T}{\hat{\mathbf{x}}}_{0}=\hat{\mathbf{x}}_{0}$.

Conversely, suppose $\hat{\mathbf{x}}(t)$ solves the PIE. Define $\mathbf{x}(t)=\mathcal{T}\hat{\mathbf{x}}(t)$. Then by Theorem 10, $\mathbf{x}(t)\in X$ and

[TABLE]

as desired. Furthermore, ${\mathbf{x}}(0)=\mathcal{T}\hat{\mathbf{x}}(0)=\mathcal{T}\hat{\mathbf{x}}_{0}={\mathbf{x}}_{0}$. $\square$

Note 7.1**.**

While the conversion formulae in Eqns. 9 are relatively complex, this is because they encompass a very large class of PDEs and must account for every case. Individual PIE representations of specific PDEs, by contrast are typically rather simple. In the following subsection, we demonstrate one such representation.

7.2 PIE Representation of the E-B Beam

To illustrate the PIE representation, we again consider the Euler-Bernoulli beam model, using the standardized PDE representation of Subsection 3.3. Applying the formulae in Eqns. (9), we obtain the PIE $\{\mathcal{T},\mathcal{A}\}$ where

[TABLE]

7.3 Stability Equivalence for PDEs and PIEs

Lemma 15**.**

The PDE defined by $\{A_{i},X\}$ is exponentially stable in $X$ with constants $K,\gamma>0$ if and only if the PIE defined by $\{\mathcal{T},\mathcal{A}\}$, where $\{\mathcal{T},\mathcal{A}\}$ are as defined in Eqn. (10), is exponentially stable with constants $K,\gamma>0$.

{pf}

Suppose the PDE defined by $\{A_{i},X\}$ is exponentially stable with constants $K,\gamma>0$. Then for any ${\mathbf{x}}_{0}\in X$, any solution ${\mathbf{x}}$ of the PDE defined by $\{A_{i},X\}$ satisfies $\left\lVert{\mathbf{x}}(t)\right\rVert_{X}\leq K\left\lVert{\mathbf{x}}_{0}\right\rVert_{X}e^{-\gamma t}$. Now for $\hat{\mathbf{x}}_{0}\in L_{2}^{n_{0}+n_{1}+n_{2}}$, let $\hat{\mathbf{x}}$ be a solution of the PIE defined by $\{\mathcal{T},\mathcal{A}\}$. Define ${\mathbf{x}}_{0}:=\mathcal{T}\hat{\mathbf{x}}_{0}\in X$ and ${\mathbf{x}}(t):=\mathcal{T}\hat{\mathbf{x}}(t)$. Then by Lemma 14, ${\mathbf{x}}(t)$ satisfies the PDE defined by $\{A_{i},X\}$ with initial condition ${\mathbf{x}_{0}}$. Therefore, by Theorem 10,

[TABLE]

Conversely, suppose the PIE defined by $\{\mathcal{T},\mathcal{A}\}$ is exponentially stable with constants $K,\gamma>0$. Then for any $\hat{\mathbf{x}}_{0}\in L_{2}$, any solution $\hat{\mathbf{x}}$ of the PIE defined by $\{\mathcal{T},\mathcal{A}\}$ satisfies $\left\lVert\hat{\mathbf{x}}(t)\right\rVert_{L_{2}}\leq K\left\lVert\hat{\mathbf{x}}_{0}\right\rVert_{L_{2}}e^{-\gamma t}$. Now for $\mathbf{x}_{0}\in X$, let $\mathbf{x}$ be a solution of the PDE defined by $\{A_{i},X\}$. Define

[TABLE]

Then by Lemma 14, $\mathbf{x}(t)=\mathcal{T}\hat{\mathbf{x}}(t)$ and $\hat{\mathbf{x}}(t)$ satisfies the PIE defined by $\{\mathcal{T},\mathcal{A}\}$ with initial condition $\hat{\mathbf{x}}_{0}$. Therefore, by Theorem 10,

[TABLE]

$\square$

Having shown that PIEs are equivalent to PDEs, we now proceed to define a Linear PI Inequality (LPI), whose feasibility guarantees exponential stability a PDE in standardized form.

8 Lyapunov Stability as an LPI

Using the 3-PI algebra, we may now succinctly represent our Lyapunov stability conditions. The procedure is relatively straightforward.

Theorem 16**.**

Suppose there exist $\epsilon,\delta>0$, $N_{0}:[a,b]\rightarrow\mathbb{R}^{n\times n}$, $N_{1},N_{2}:[a,b]^{2}\rightarrow\mathbb{R}^{n\times n}$ such that for $\mathcal{P}:=\mathcal{P}_{\{N_{0},N_{1},N_{2}\}}$, $\mathcal{P}=\mathcal{P}^{*}\geq\alpha I$ and

[TABLE]

where $\mathcal{T}$ and $\mathcal{A}$ are as defined in Eqn. (10). Then any solution, $\mathbf{x}(t)$ of the PDE defined by $\{A_{i},X\}$ satisfies

[TABLE]

where $\zeta=\left\lVert\mathcal{P}\right\rVert_{\mathcal{L}(L_{2})}$.

{pf}

Suppose $\hat{\mathbf{x}}$ solves the PIE defined by $\{\mathcal{T},\mathcal{A}\}$ for some $\hat{\mathbf{x}}_{0}$. Consider the candidate Lyapunov function defined as

[TABLE]

The derivative of $V$ along solution $\hat{\mathbf{x}}$ is

[TABLE]

Recall $\left\lVert\mathcal{P}\right\rVert_{\mathcal{L}(L_{2})}=\zeta$. Then by a standard application of Gronwall-Bellman, we have

[TABLE]

Now for any solution, $\mathbf{x}$ of the PDE defined by $\{A_{i},X\}$ with initial condition $\mathbf{x}_{0}$, $\mathbf{x}(t)=\mathcal{T}\hat{x}(t)$ for solution of the PIE with initial condition $\hat{\mathbf{x}}_{0}$ where $\mathbf{x}_{0}=\mathcal{T}\hat{\mathbf{x}}_{0}$. Thus

[TABLE]

$\square$

Note 8.1**.**

Theorem 16 proves exponential stability of the PDE with respect to the $L_{2}$ norm and not the $X$-norm. While it is possible to formulate a PIE for stability in the $X$-norm, this would differ from most existing results and the literature and hence is omitted. Note, however, that for any $\mathbf{x}\in L_{2}$, $\mathcal{T}\mathbf{x}=0$ if and only if $\mathbf{x}=0$ (modulo a set of zero measure).

Note 8.2**.**

Theorem 16 is equivalent to the Lyapunov inequality for PDEs with the restriction that the Lyapunov operator be a PI operator. This, in turn, may be interpreted as a dissipativity condition on the generator. Such conditions are sometimes enforced using multiplier approaches as in, e.g. [19], and have been shown to be necessary and sufficient for stability of infinite-dimensional systems, as in [32, 16]. Note that the constraint that the operator $\mathcal{P}$ be self-adjoint is not conservative as any Lyapunov function defined by a non-self-adjoint operator admits a representation using a self-adjoint operator.

Theorem 16 poses a convex optimization problem, whose feasibility implies stability of solutions of a coupled linear PDE. We refer to such optimization problems as Linear PI Inequalities (LPIs). Solving an LPI requires parameterizing the 3-PI operator $\mathcal{P}$ using polynomials and enforcing the inequalities using LMIs. In the following section, we briefly introduce a method of enforcing positivity of a 3-PI operator using LMI constraints.

9 Solving the Stability LPI via PIETOOLS

In Section 8, we formulated the question of Lyapunov stability as an LPI. In this section, we will we propose a general form of LPI and show how these convex optimization problems can be solved under the assumption that all 3-PI operators are parameterized by polynomials.

For given 3-PI operators $\{\mathcal{E}_{ij},\mathcal{F}_{ij},\mathcal{G}_{i}\}$ and linear operator $\mathcal{L}$, a Linear PI Inequality (LPI) is a convex optimization of the form

[TABLE]

LPIs of the form of Eqn. (12) can be solved directly using PIETOOLS [31]. Composition and adjoint are algebraic operations on the 3-PI parameters and are computed using the formulae in Section 4. Positivity is enforced using an LMI constraint as described in the following subsection.

9.1 Enforcing Positivity of 3-PI Operators

In this subsection, for a given self-adjoint 3-PI operator with polynomial parameters ($\{N_{i}\}$), we given an LMI constraint on the coefficients of the polynomials $\{N_{i}\}$ which enforces a constraint of the form $\mathcal{P}_{\{N_{i}\}}\geq 0$. Specifically, the following proposition (a slight modification of the result in [33]) gives necessary and sufficient conditions for a 3-PI operator to have a 3-PI square root.

Proposition 17**.**

For any bounded functions $Z(s)$, $Z(s,\theta)$, and $g$, where $g$ is scalar and $g(s)\geq 0$ for all $s\in[a,b]$ and

[TABLE]

where

[TABLE]

we have $\mathcal{P}_{\{N_{i}\}}^{*}=\mathcal{P}_{\{N_{i}\}}$ and $\left\langle\mathbf{x},\mathcal{P}_{\{N_{i}\}}\mathbf{x}\right\rangle_{L_{2}}\geq 0$ for all $\mathbf{x}\in L_{2}[a,b]$.

{pf}

It is relatively easy to show that $\{N_{i}\}$ satisfy Equation (7) with $\{\hat{N}_{i}\}=\{N_{i}\}$. Therefore, by Lemma 6 $\mathcal{P}_{\{N_{i}\}}$ is self adjoint. Now define the operator

[TABLE]

Then by expanding out the composition formulae, we find $\mathcal{P}_{\{N_{i}\}}=\mathcal{Z}^{*}P\mathcal{Z}$ and since $P\geq 0$, $P=(P^{\half})^{T}P^{\half}$ for some $P^{\half}$. Thus

[TABLE]

$\square$

Note that Prop. 17 does not ensure that $\mathcal{P}_{\{N_{i}\}}$ is coercive. To obtain a coercive operator, one must add a coercive term of the form $\mathcal{P}_{\{\epsilon I,0,0\}}$.

When we desire the $\{N_{i}\}$ to be polynomial, we may choose $Z$ to be the vector of monomials of bounded degree, $d$. For $g(s)=1$, the operators are positive on any domain. However, for $g(s)=(s-a)(b-s)$ the operator is only positive on the given domain $[a,b]$. For the most accurate results, we combine both choices of $g$. For notational convenience, we define the set of functions which satisfy Prop. 17 in this way. Specifically, we denote $Z_{d}(x)$ as the matrix whose rows are a vector monomial basis for the vector-valued polynomials of degree $d$ or less and define the cone of positive operators with polynomial multipliers and kernels associated with degree $d$ as

[TABLE]

The dimension of the matrices $M_{i}$ and $N_{i}$ should be clear from context. The constraint $\mathcal{P}_{\{R_{i}\}}\in\Omega_{d}$ is then an LMI constraint on the coefficients of the polynomials $\{R_{i}\}$ and guarantees that $P_{\{R_{i}\}}\geq 0$. A Matlab toolbox (PIETOOLS) for setting up and solving LPIs based on Prop. 17 has recently been proposed and is discussed in Subsection 9.3.

9.2 The Degree-Bounded Stability Test

By restricting the degree of the polynomial parameters, $\{N_{i}\}$, we obtain a PIETOOLS-based LMI which enforces the LPI conditions of Theorem 16.

Theorem 18**.**

For any $d\in\mathbb{N}$, suppose there exist $\epsilon,\delta>0$, and polynomials $N_{0}:[a,b]\rightarrow\mathbb{R}^{n\times n}$, $N_{1},N_{2}:[a,b]^{2}\rightarrow\mathbb{R}^{n\times n}$ such that

[TABLE]

and

[TABLE]

where $\mathcal{T}$ and $\mathcal{A}$ are as defined in Eqn. (10). Then any solution the PDE defined by $\{A_{i},X\}$ is exponentially stable in $L_{2}$.

Note that as mentioned in the previous subsection, the constraint $\in\Omega_{d}$ is an LMI constraint.

9.3 PIETOOLS Implementation

In this subsection, we give sample code using the PIETOOLS toolbox which verifies that the conditions of Theorem 18 are satisfied.

A detailed manual for the PIETOOLS toolbox can be found in [31]. This toolbox allows for declaration and manipulation of 3-PI operators and 3-PI decision variables and enforcement of LPI constraints. PIETOOLS uses aspects of the SOSTOOLS LMI conversion process and pvar polynomial objects as defined in MULTIPOLY. PIETOOLS defines the opvar class of PI operators and overloads the multiplication (*), addition (+) and adjoint (’) operations using the formulae in Lemma 4 and Lemma 6. Concatenation, and scalar multiplication are likewise defined so that 3-PI operators can be treated in a similar manner to matrices.

To facilitate implementation of the conditions of Theorem 18, we have created the script PIETOOLS_PDE, which is distributed with the PIETOOLS toolbox. To use this script only requires the user to define the standardized form of the PDE - as illustrated in Step (3) below. Specifically, the user must define n0,n1,n2,A0,A1,A2,B, although A2 may be omitted if n2=0. The user specifies that a stability test is desired by setting stability=1 and can specify the desired accuracy through settings_PIETOOLS scripts, although by default we use settings_PIETOOLS_light script, which corresponds to $d=1$. An overview of the steps included in the PIETOOLS_PDE script are provided below along with a brief description of each step.

Define independent polynomial variables. These are the spatial variables in the PDE.

[TABLE] 2. 2.

Initialize an optimization problem structure, $X$.

[TABLE] 3. 3.

Define the standardized PDE representation and use the provided script to construct the PIE. The script convert_PIETOOLS_pde then converts the standardized PDE to a PIE.

[TABLE] 4. 4.

Declare the positive 3-PI operator $\mathcal{P}$ and add inequality constraints. $n$ is state dimension, $I$ is the interval $[a,b]$, and $d$ is the degree of the polynomial parameters in $\mathcal{P}$. This step is encoded in the script executive_PIETOOLS_stability and is executed automatically if the user has declared the option stability=1.

[TABLE] 5. 5.

Call the SDP solver.

[TABLE] 6. 6.

Get the solution. P_s is the 3-PI operator, $\mathcal{P}$.

[TABLE]

Not that the degree, $d$, enters at Step (4) and is defined in the settings script, which defaults to settings_PIETOOLS_light ($d=1$). If higher degree is required, the setting may be changed manually or using the settings_PIETOOLS_heavy ($d=2$) script. Instructions for declaring the PDE are included in the header to PIETOOLS_PDE.

10 Numerical Tests of Accuracy and Scalability

In this section, we examine the accuracy and computational complexity (scalability) of the proposed stability analysis algorithm by applying Theorem 18 to several well-studied and relatively trivial test cases. The algorithms are implemented using the PIETOOLS toolbox described in the previous section, and use the settings_PIETOOLS_light ($d=1$) option. All computation times are listed for an Intel i7-6950x processor with 64GB RAM and only account for time taken to solve the resulting LMI using Sedumi, excluding time taken for problem setup and polynomial manipulations. In cases where the limiting value of a parameter is listed for which the system is stable, the limiting value was determined using a bisection on that parameter.

Example 1: We begin with several variations of the diffusion equation. The first is adapted from [26].

[TABLE]

where $x(0)=x(1)=0$ and which is known to be stable if and only if $\lambda<\pi^{2}=9.869604\cdots$. For the choice of $d=1$ in Thm. 18, the algorithm is able to prove stability for $\lambda\leq 9.8696$ with a computation time of .54s.

Example 2: The second example from [27] is the same, but changes the boundary conditions to $x(0)=0$ and $x_{s}(1)=0$ and is unstable for $\lambda>2.467$. For $d=1$, the algorithm is able to prove stability for $\lambda\leq 2.467$ with identical computation time.

Example 3: The third example from [22] is not homogeneous

[TABLE]

where $x(0)=0$ and $x_{s}(1)=0$ and was estimated numerically to be unstable for $\lambda>4.65$. For $d=1$, the algorithm is able to prove stability for $\lambda\leq 4.65$ with similar computation time (compare to $\lambda=4.62$ in [22]).

Example 4: In this example from [27], we have

[TABLE]

with $x(0)=0$ and $x_{s}(1)=0$. In this case, using $d=1$, we were able to prove stability for any tested value of $R$ (vs. $R\leq 21$ in [27]) with a computation time of $4.06s$. No upper limit was found.

Example 5: For our last numerical comparison, we consider some of the recent literature on coupled linear hyperbolic systems [34, 35, 9], often representing conservation or balance laws. Although there are several variations of the problem formulation, we consider the recent work of [9], as it seems to contain the most accurate results.

[TABLE]

with boundary conditions $x_{1}(0)=qx_{2}(0)$ and $x_{2}(1)=\rho x_{1}(1)$. In this case, we have

[TABLE]

Using $d=1$, $r_{1}=0$, $r_{2}=1.1$, $\sigma_{1}=1$, $q=1.2$, by gridding the parameters $\sigma_{2}$ and $\rho$, we are able to verify stability for all stable parameter values indicated in Figure 5 in [9]. For example, at $\rho=-.4$, we were able to prove stability for $\sigma_{2}\leq 1.048$.

Example 6 (Scalability): Finally, we explore computational complexity using a simple $n$-dimensional diffusion equation

[TABLE]

where $x(t,s)\in\mathbb{R}^{n}$. We then evaluate the computation time to perform the feasibility test for different size problems, from $n=1$ to $n=40$, choosing $d=1$ - See Fig. 1. Note that no factors other than $d$ influence computation time and the result is always stability.

11 Illustrations of Beam and Wave Equations

In Section 10, we demonstrated that the proposed stability test has no obvious conservatism by finding parameter values corresponding to the stability limit for several well-studied examples. However, the representation of these PDEs in the generalized PDE for of Eqn. (3) was straightforward. In this section, we provide some guidance on how the user might identify the $A_{i}$ and $B$ matrices in Eqn. (3) for several less obvious examples - focusing on four well-known wave and beam examples. The beam examples are particularly interesting in that (to the best of our knowledge) they have not previously been analyzed using LMI-based methods. As we proceed, we call particular attention to the following two questions.

• •

What are the states?

• •

What are the boundary conditions?

Choice of State: Prior to the introduction of state-space, ODEs would often be represented using scalar equations. For example, the spring-mass:

[TABLE]

is a scalar ODE. To represent this in the vector-valued state-space framework, we use $x_{1}$ and define an auxiliary state $x_{2}=\dot{x}$. Similarly, PDEs are often represented as scalar equations using higher-order time derivatives (e.g. The wave equation is $\ddot{w}=w_{xx}$). The standardized PDE representation in Eqn. (3), however, uses only first-order time derivatives. Furthermore, as discussed in Subsection 3.3, the use of the standardized representation occasionally involves loss of some state information and may affect the question of stability. Specifically, the exponential stability criterion in Theorem 18 implies all states decay exponentially. For example, If a PDE is $L_{2}$-stable in $u$, but not $u_{s}$, then if $u_{s}$ is included in the standardized representation, the PIE stability analysis will not be able to verify stability.

Boundary Conditions: Identification of a sufficient number of boundary conditions in the universal framework is particularly important. For the $B$ matrix to have sufficient rank, the solution must be uniquely defined. One consideration to be aware of is that when we introduce additional variables to eliminate higher-order time-derivatives, these new variables must also have associated boundary conditions. This is typically solved by differentiating the original boundary conditions in time to obtain boundary conditions for the new variables.

In the following examples, we illustrate the process of choosing state and constructing the $A_{i}$ and $B$ matrices.

11.1 Beam Equation Examples

We first both the Euler-Bernoulli (E-B) and Timoshenko (T) models beam equations. This case is particularly interesting, as the E-B model is fundamentally diffusive and the T model has hyperbolic character. Furthermore, both these models are known to be energy-conserving[19], meaning that they are stable, but not exponentially stable.

Euler-Bernoulli: In this first case, we simply recall our formulation of the cantilevered E-B beam from Subsection 3.3:

[TABLE]

where $A_{0}=A_{1}=0$, $n_{2}=2$, and $n_{0}=n_{1}=0$. The boundary conditions take the form

[TABLE]

Entering $\{A_{i},B\}$ into the script PIETOOLS_PDE, we find the E-B beam is neutrally stable (using $\delta=0$ in Thm. 18) for any tested value of $c>0$. However, when $\delta>0$, the code is unable to find a Lyapunov function, indicating this formulation is not exponentially stable (as expected). Note that to set $\delta=0$, we modify the script using the command epneg=0.

Timoshenko Beam We now consider the Timoshenko beam model where, for simplicity, we set $\rho=E=I=\kappa=G=1$:

[TABLE]

with boundary conditions of the form

[TABLE]

Our first step is to eliminate the second-order time-derivatives, and hence we choose $u_{1}=\dot{w}$ and $u_{3}=\dot{\phi}$. Using the boundary conditions as a guide, we choose the remaining states as $u_{2}=w_{s}-\phi$ and $u_{4}=\phi_{s}$. Note that this choice of states is a natural set of coordinates as the Timoschenko beam is known to be energy conserving with respect to these states [19]. In summary, we have

[TABLE]

This gives us 4 first order boundary conditions

[TABLE]

Reconstructing the dynamics, we now have

[TABLE]

Expressing this in our standard form we have the purely hyperbolic construction

[TABLE]

where $A_{2}=[\,]$ and $n_{0}=n_{2}=0$ and $n_{1}=4$. The $B$ matrix is then

[TABLE]

where $B$ has row rank $n_{1}=4$ and satisfies Eqn. (5). The script PIETOOLS_PDE indicates this system is neutrally stable (using $\delta=0$ in Thm. 18). However, when $\delta>0$, the code is unable to find a Lyapunov function, indicating this formulation is not exponentially stable (as expected).

11.2 Wave Equation with Boundary Feedback Examples

In this subsection, we consider wave equations attached at one end and free at the other with damping at the free end. This is a well-studied problem for which numerous stability results are available in the literature [36, 37]. The simplest formulation is

[TABLE]

As with the beam examples, this has a purely hyperbolic formulation. Guided by the boundary conditions, we choose

[TABLE]

This yields

[TABLE]

where $A_{0}=0$, $A_{2}=[\,]$ $n_{1}=n_{2}=0$ and $n_{1}=2$. The boundary conditions are now

[TABLE]

This formulation is computed to be exponentially stable in the given state $u_{t},u_{s}$ for any tested value of $k>0$. We now consider a variation on this formulation.

Diffusive Formulation As a variation, we consider a non-diffusive formulation from [36] which was shown to be asymptotically stable in the state $u$ for $a^{2}+k^{2}>0$.

[TABLE]

In this case, we are forced to choose the variables $u_{1}=u_{t}$ and $u_{2}=u$ yielding the diffusive formulation

[TABLE]

where $A_{1}=0$, $n_{0}=0$, $n_{1}=1$, and $n_{2}=1$. Note in this case that the boundary conditions on $u_{1}$ force us to consider this a hyperbolic state and the boundary conditions on $u_{2}$ make this a diffusive state! These boundary conditions are now expressed as

[TABLE]

Computation indicates this model is neutrally stable, but not exponentially stable in the given state - a result confirmed in [36, 37].

12 Conclusion

In this paper, we have shown how to use LMIs to accurately test stability of a large class of coupled linear PDEs. To achieve this result, we have shown how to convert well-posed coupled linear PDEs - defined on state $\mathbf{x}_{p}\in X$, with associated boundary conditions and continuity constraints - to Partial-Integral Equations (PIEs) with state $\mathbf{x}_{f}\in L_{2}$ - a formulation which is defined using the algebra of 3-PI partial-integral operators and which does not require boundary conditions or continuity constraints on $\mathbf{x}_{f}$. We have shown that stability of PDEs can be reformulated as a Linear PI Inequality (LPI) expressed using 3-PI operators and operator positivity constraints. We have shown how to parameterize 3-PI operators using polynomials and how to enforce positivity of 3-PI operators using LMI constraints on the coefficients of these polynomials. We have used the Matlab toolbox PIETOOLS to solve the resulting LPIs and applied the results to a variety of numerical examples. The numerical results indicate little or no conservatism in the resulting stability conditions to several significant figures even for low polynomial degree. By conversion of LMIs developed for ODEs to LPIs, it is possible that these results can be extended to: PDEs with uncertainty; $H_{\infty}$-gain analysis of PDEs; $H_{\infty}$-optimal observer synthesis for PDEs; and $H_{\infty}$-optimal control of PDEs. In addition, it is possible that the framework may be extended to multiple spatial dimensions using the multivariate representation proposed in [38].