{'fraction_non_alphanumeric': 0.06094400376973219, 'fraction_numerical': 0.019869630095028665, 'mean_word_length': 3.849200304645849, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Stability Analysis of Parabolic Linear PDEs with Two Spatial Dimensions Using Lyapunov Method and SOS

In this paper, we address stability of parabolic linear Partial Differential Equations (PDEs). We consider PDEs with two spatial variables and spatially dependent polynomial coefficients. We parameterize a class of Lyapunov functionals and their time derivatives by polynomials and express stability as optimization over polynomials. We use Sum-of-Squares and Positivstellensatz results to numerically search for a solution to the optimization over polynomials. We also show that our algorithm can be used to estimate the rate of decay of the solution to PDE in the L2 norm. Finally, we validate the technique by applying our conditions to the 2D biological KISS PDE model of population growth and an additional example.


I. INTRODUCTION
Stability analysis and controller design for Partial Differential Equations (PDEs) is an active area of research [1], [2]. One approach to stability analysis of PDEs is to approximate the PDEs with Ordinary Differential Equations (ODEs) using, e.g. Galerkin's method or finite difference, and then apply finite-dimensional optimal control methods, e.g. [3]- [6]. In this paper we consider stability analysis without discretization. Specifically, we use Sum-of-Squares (SOS) optimization to construct Lyapunov functionals for parabolic PDEs with two spatial variables and spatially dependent coefficients.
It is well-known that existence of a Lyapunov function for a system of ODEs or PDEs is a sufficient condition for stability. For example, [7] uses a Lyapunov approach and Linear Operator Inequalities (LOIs) to provide sufficient conditions for exponential stability of a controlled heat and delayed wave equations. Another method based on Lyapunov and semigroup theories was applied in [8] for analysis of wave and beam PDEs with constant coefficients and delayed boundary control. In [9] a Lyapunov based analysis of semilinear diffusion equations with delays gave stability conditions in terms of Linear Matrix Inequalities (LMIs). Extensive examples of applying the backstepping method to the boundary control of PDEs can be found in [10]- [14]. Briefly speaking, backstepping uses a Volterra operator to search for an invertible mapping from the original PDE to a chosen "target" PDE, known to be stable. In order to find such mapping, one has to solve analytically or numerically a PDE for Volterra operator's kernel. If the mapping is found, it provides the boundary control law. Applications for twodimensional cases were discussed in [15] and [16]. However, backstepping requires us to guess on the target PDE and solve the PDE for kernel, which may be a challenging task for PDEs with two spatial variables and spatially dependent coefficients. Moreover, backstepping cannot be used for stability analysis in the absence of a control input.
Note that SOS has been previously applied in [17] and [18] to find Lyapunov functionals for 1D parabolic PDE. Input-Output analysis of PDE systems with SOS implementation is discussed in [19]. Examples of using SOS in controller design for one-dimensional PDEs can be found in [20]- [22].
The main contribution of this paper is to provide an algorithm for stability analysis of 2D parabolic PDEs with spatially varying coefficients. We search for polynomial s which defines a Lyapunov functional of the form where u is the solution of the PDE, t represents time, x and Ω are the spatial variable and domain. We introduce stability conditions in the form of parameter dependent LMIs, which certify Lyapunov inequalities, i.e.
V (u(t, ·)) > 0 and d dt [V (u(t, ·))] < 0 for all t > 0. (1) There exist several algorithms for optimization over polynomials which can be applied to (1) in order to find V . See e.g., [24] for a survey on algorithms for optimization of polynomials. In this paper, we apply SOS and the Positivstellensatz results to find a solution to parameter dependent LMIs. For details on Positivstellensatz see [23] or [25]. We use the MATLAB toolbox SOSTOOLS (see [26] for details on SOSTOOLS) and SeDuMi, a well-known semidefinite programming solver, to solve the SOS optimization problem. The paper is organized as follows. In Section 2 we specify the notation and provide some background. Lyapunov stability conditions for parabolic PDEs are given in Section 3. We demonstrate the proposed method in Section 4 and present SOS optimization problems in Section 5. Finally, we discuss numerical implementations in Section 6 and conclude the paper in Section 7 with ideas about future research steps.

A. Notation
N is the set of natural numbers and N 0 := N ∪ {0}. R n and S n are the n-dimensional Euclidean space and space of n × n real symmetric matrices. For x ∈ R n , let x T denote transposed x and x i ∈ R is the i-th component of x. · 1 is a norm on R n , defined as x 1 := n i=1 |x i |. For X ∈ S n , X ≤ 0 means that X is negative semidefinite. The symbol * will denote the symmetric elements of a symmetric matrix.
For Ω ⊆ R n and f : Ω → R let f (x) stand for f (x 1 , ..., x n ) and Ω f (x) dx represent an integral of f over Ω with dx := dx 1 dx 2 ...dx n .
Let N n 0 := {α ∈ R n : α i ∈ N 0 }. A vector α ∈ N n 0 is called multi-index. For l ∈ N define the set For α ∈ N n 0 , x ∈ R n and g : R n → R partial derivative Note that ∂ 0 for any i ∈ {1, ..., n}. We will also use classical notation, u x1x2 (t, If for a function f : Ω → R and some α ∈ N n 0 derivative Further we use W 2,2 (Ω). It is one of Banach spaces W k,p (Ω) which denote Sobolev spaces of functions u : Ω → R with D α [u] ∈ L p (Ω) for all α ∈ Q n k , where Q n k is defined as in (2) and norm u k,p := where L p (Ω) stands for the space of Lebesgue-measurable functions g : Ω → R with norm, for p ∈ N It is known that for continuous functions u : [0, ∞) → W 2,2 (Ω) and V : W 2,2 (Ω) → R their composition (V • u) : [0, ∞) → R is also continuous and the upper right-hand derivative D + t V (u(t)) is defined by

B. Polynomials
For a multi-index α ∈ N n 0 and x ∈ R n , let .. x αn n . Then x α is a monomial of degree α 1 ∈ N 0 . A polynomial is a finite linear combination of monomials p(x) := α p α x α , where the summation is applied over a given finite set of multi-indexes α and p α ∈ R denotes the corresponding coefficient. Let R[x] stand for polynomial ring in n variables x 1 , ..., x n with coefficients in R. The degree of a polynomial p ∈ R[x] is the largest degree among all monomials, and is denoted by deg(p) Polynomial matrices and SOS matrices are defined in a similar manner (See, e.g. [27]). We use [S(x)] to denote the set of real symmetric SOS polynomial matrices. If M ∈ [S m (x)] then for all x ∈ R n , M (x) ∈ S m and M (x) ≥ 0.

C. Comparison Lemma
Recall the comparison principle from, e.g. [28]. Lemma 1: Consider the scalar differential equation x) is continuous in t and locally Lipschitz in x, for all t ≥ 0 and all x ∈ J ⊂ R. Let [t 0 , T ) (T could be infinity) be the maximal interval of existence of the solution u, and suppose u(t) ∈ J for all t ∈ [t 0 , T ). Let v be a continuous function whose

III. LYAPUNOV STABILITY FOR PARABOLIC PDE
First consider the following general form of parabolic PDEs. For all t ∈ (0, ∞) and x ∈ Ω ⊂ R n , are partial derivatives as denoted in (3). Assume that solutions to (4) exist, are unique and depend continuously on initial conditions. This implies that the function is continuously differentiable in t and for each t ∈ [0, ∞), u(t, ·) ∈ W 2,2 (Ω).
Theorem 1: Let there exist continuous V : L 2 (Ω) → R, l, m ∈ N and b, a > 0 such that for all v ∈ L 2 (Ω). Furthermore, suppose that there exists c ≥ 0 such that for all t ≥ 0 the upper right-hand derivative where u satisfies (4). Then Proof: Let conditions of Theorem 1 be satisfied. Since Dividing both sides of the second inequality in (7) After multiplying both sides of (8) by −c, we have From (6) and (9) it follows that To use the comparison principle, consider the ODE where t ∈ (0, ∞) and function φ : [0, ∞) → R is continuous. The well-known solution for (11) is for all t ≥ 0. Applying Lemma 1 for (10) and (11) for all t ≥ 0. Substituting t = 0 in the second inequality of (7) implies Combining the first inequality of (7) with (13) and (12) gives Dividing (14) by a and taking the l th root results in

IV. STABILITY TEST EXPRESSED AS OPTIMIZATION OVER POLYNOMIALS
In this paper, we focus on PDEs of the following form. For all t > 0 and x ∈ Ω : where a, b, c, d, e, f ∈ R[x]. Assume that solution to (15) exists, is unique and depends continuously on initial conditions. For each t ≥ 0 suppose u(t, ·) ∈ W 2,2 (Ω). Furthermore, let u satisfy zero Dirichlet boundary conditions, i.e.
where s ∈ R[x]. Using u(t, ·) for v in (17) and differentiating the result with respect to t gives Substituting for u t (t, x) from (15) into (18) implies where Note that, based on section 5.2.3 of [29], we have for all x ∈ Ω. We used property (20) to define I 2 and I 3 .
Alternatively, I 5 can be formulated as where for all x ∈ Ω and Using integration by parts and boundary conditions (16), I 1 can be rewritten as follows.
where for all x ∈ Ω Following steps of (23) for I 2 , I 3 and I 4 , we get where U is defined as in (22) and for all x ∈ Ω By combining (19), (21), (23) and (24) it follows that d dt where for all x ∈ Ω

A. Spacing functions
If Q(x) ≤ 0 for all x ∈ Ω, then the time derivative in (25) is clearly non-positive for all t > 0. However, such a condition on Q is conservative. To decrease that conservatism, we introduce matrix valued functions Υ i such that and, therefore, can be added to Q without altering the integral. Υ i are called spacing functions. We parameterize Υ i by polynomials p i . The idea came from [27].
The following holds for any p 1 ∈ R[x], because of the boundary conditions (16).
Using the chain rule, we have where for all x ∈ Ω Combining (27) and (28) results in Likewise in (27), because of the boundary conditions (16), the following holds for any .
Following steps of (28) Similarly to (27), the following is true for any .
Note that the left-hand side of (33) can be written as follows.
where we need property (20). Applying integration by parts to the second integral of the right-hand side of the last equation in (34) results in From (34) and (35) the following holds.
Proof: If conditions of Theorem 2 are satisfied, let Since s(x) ≥ θ > 0 for all x ∈ Ω, then b, a > 0 and the following holds for all v ∈ L 2 (Ω).
Then for all t > 0 solution to (15) satisfies where a, b are defined as in (46). Proof: Under the assumptions of Theorem 3, (47) holds. With (48) and (22) we can write Since M (x) + γS(x) ≤ 0 for all x ∈ Ω, it holds that Since γ is a scalar, (51) can be easily satisfied as follows.

VI. NUMERICAL VERIFICATION OF THE PROPOSED METHOD A. KISS model
We applied the method described in Sections (IV) and (V) to stability of the biological KISS PDE named after Kierstead, Slobodkin and Skelam, which describes population growth on a finite area. For more details see [30]. The system is modeled by the following PDE.
where h, r > 0, x ∈ Ω ⊂ R 2 and scalar function u satisfies zero Dirichlet boundary conditions.
It is claimed in [30] that if Ω is a square with edge of length l, then    Table (I). Now we choose l = 1, h = 2 and r = 4. Using a bisection search over γ, we determine the maximum γ for which the SOS problem in Theorem 5, with (58), may be shown to be feasible. Results for different deg(s) are presented in Table  (II). Using finite difference scheme, we numerically solve (55) with u(0, x) = 10 3 x 1 x 2 (1 − x 1 )(1 − x 2 ). Plots of log 10 ( u(t, ·) L2 ) versus t, using a numerical solution, and bounds on log 10 ( u(t, ·) L2 ), given by the proposed method for different deg(s), are presented in Fig. (1). These plots allow us to determine γ by examining the rate of decrease in the L 2 norm. Plots are aligned at t = 0 in order to better compare our SOS estimates of γ to the estimate of γ derived from numerical simulation as a function of increasing deg(s).

VII. CONCLUSION AND FUTURE WORK
In this paper we have presented a method which allows us to search for Lyapunov functionals for parabolic linear PDEs with two spatial variables and spatially dependent coefficients. We demonstrated accuracy of the method in estimating critical diffusion for the KISS PDE and the rate of decay for the KISS PDE and for a randomly chosen PDE. In future work, we will extend the method by considering Lyapunov functionals of more complicated types, e.g. functionals with semi-separable kernels, as in [31]. In addition, we will consider alternative boundary conditions, as well as semi-linear and nonlinear 2D and 3D parabolic PDEs.