Partial Differential Equations (PDEs) with integral terms (integrated with respect to the spatial variable) typically appear in the following two scenarios: a) when the inherent physics leads to integral terms while modeling (for example, population dynamics , linear thermoelasticity , chemical reaction and transport , etc.); b) when the model involves point sensor measurements which, in practice, are weighted averages of the distributed state; and c) closed-loop dynamics of PDEs that are coupled with a boundary feedback controller or observer designed using the Backstepping methods . Various mathematical analysis methods for such PDEs were presented in [1, 5] whereas  used computational tools like LMIs for analysis. However, these methods did not allow analysis (including stability analysis) of PDEs with integral terms at the boundary and typically can only be used for a fixed type of boundary conditions.
The goal of this paper is develop a computational method to analyze stability of PDEs with integral terms in both the dynamics and the boundary conditions. To describe the type of PDEs considered in this paper, we present the following two examples:
Example 1: Consider a model for population dynamics (called the McKendrick PDE; see ) given by
where is the density of population of age at any given time . The boundary condition can be considered as an approximation of the number of newborns at time and is usually modeled as an integral boundary condition (a weighted average of the population density at different ages).
Example 2: Consider the closed loop dynamics of a reaction-diffusion equation with a boundary sensed observer. The PDE dynamics along with the observer dynamics can be written as
In , the observer gains, , were shown to be modified Bessel functions (denoted by , a convergent but infinite series) and are given by
A practical implementation of such an observer-based controller would involve an approximation of the kernels by a finite sum. While the original observer is provably stable, the stability of finite-dimensional approximation of observer is not guaranteed and must be verified — a non-trivial step. Naturally, we want to simplify this post truncation verification step, possibly, by using a computational method.
Additionally, in practice, point measurements such as (the sensor measurement that drives the observer dynamics) are not exact due to physical limitations of the sensing mechanism. A more accurate representation of such a sensor output is of the form
where is typically a Gaussian function with the center at . Thus, a practical model of the observer dynamics can be written as
Note that Partial integrodifferential equations (PIDEs) such as Volterra integral equation and Fredholm integral equation are not considered since the variable of integration for such systems is the same as the variable in time-differential. To develop computational tests for the stability of PDEs with integral terms, We will use the Partial Integral Equation (PIE) framework, introduced in , which is an approach to analysis/control PDEs by converting PDEs to PIEs and then solving analysis/control problems for the PIE using LMI based methods. This approach for analysis, although straightforward, cannot be directly used since a PIE form for PDEs with integral terms was not proposed in  or any of the subsequent works.
The modified goal of this paper, then, is to extend the class of PDE models for which we have PDE to PIE conversion formulae to include the cases defined above. We approach the incorporation of a new class of PDE models in three steps: (a) we propose a parametric representation of the PDE model class and corresponding state-space; (b) We define an appropriate state-space to be used in the corresponding PIE model; (c) We find a unitary transformation from the PIE state-space to the state space of the PDE model – proving equivalence of solutions and equivalence of stability properties. To summarize the main contribution of this paper, we provide explicit maps from the parameters of a PDE with integral terms to the parameters of the PIE where the PIE obtained will have stability properties identical to that of the PDE. Then we present the stability test for PDE with integral terms as an operator-valued optimization problem.
PIEs are explained in more detail in section II-B, however, here we briefly describe the PIE system and benefits of using the PIE framework to motivate our approach for stability analysis of PDEs with integral terms. A PIE: a) is defined by a single integro-differential equation; b) is parameterized by the algebra of Partial Integral (PI) operators; and c) can be used to represent almost any well-posed PDE model. The general form of a PIE, defined by 3-PI operators, is given by
for all where the state, admits no continuity constraints or boundary conditions. An operator is said to be a 3-PI operator, denoted if there exist and separable functions , all matrix-valued, such that
The benefits of using a PIE for analysis, instead of a PDE, is because PIEs have a state-space form similar to ODEs and are parameterized by the PI operators that share many properties of matrices (specifically, being closed algebraically). Consequently, many numerical methods designed for control, analysis, and simulation of ODEs in state-space form can be extended to PIEs. For example, computational tests for analysis and control of PIEs were presented in [9, 10, 3], etc.
To summarize the organization of the paper, first, we define the algebra of PI operators in Section II-A and define the class of PIE models, along with a definition of a solution in Section II-B. Then we introduce the class of PDE systems in Section III – including a definition of the solution of a PDE. Next, for a given PDE, we propose an admissibility condition that guarantees the existence of an equivalent PIE representation and present the formulae for constructing the PIE representation of the PDE (See Section IV). Finally, in Section VI, we consider several specific PDE models – using the PIE representation and PIETOOLS to simulate the PDE and prove stability.
We use the notation
to represent the zero matrix of dimensionand . Similarly,
is the identity matrix of dimension. We use and for the zero and identity matrix when the dimensions are clear from the context. is the set of non-negative real numbers.
is the Hilbert space of
-dimensional vector-valued Lebesgue square-integrable functions defined on the intervaland is equipped with the standard inner product. is the Banach space of -dimensional essentially bounded measurable matrix-valued functions on
equipped with the essential supremum singular value norm.
Normal font or typically implies that or is a scalar or finite-dimensional vector (e.g. ), whereas the bold font, or , typically implies that or is a scalar or vector-valued function (e.g. ). For a suitably differentiable function, of spatial variable , we use to denote the -th order partial derivative . For a suitably differentiable function of time and possibly space, we denote . We use to denote the Sobolev spaces
with inner product
For brevity, we omit the domain and simply write or when clear from the context.
Ii-a PI Operators: A -algebra of bounded linear operators on
The PI algebras are parameterized classes of operators on . The algebra of 3-PI operators which parameterize maps from .
Definition 1 (Separable Function).
We say a function , is separable if there exist and functions and such that .
Definition 2 (3-PI operators, ).
Given and separable functions , we define the operator for as
Furthermore, we say an operator, , is 3-PI of dimension , denoted , if there exist functions and separable functions such that .
Additionally, we define as the set of positive definite 3-PI operators where the positivity is defined with respect to -inner product, i.e, .
Ii-B Partial Integral Equations
A Partial Integral Equation (PIE) is an extension of the state-space representation of ODEs (vector-valued first-order differential equations on ) to spatially-distributed states on . A PIE model is parameterized by 3-PI operators as
where and are 3-PI operators.
Unlike PDE models, a PIE does not allow for spatial derivatives – only a first-order derivative with respect to time. In particular, the state of the PIE model, is not differentiable and consequently, no boundary conditions are possible in the PIE framework. Finally, given a PIE, we require differentiability of the solution with respect to the -norm which is defined as
Definition 3 (Definition of solution of a PIE).
For given inputs , we say that satisfies the PIE defined by with initial condition if
is Frechét differentiable with respect to the -norm almost everywhere on
Eq. (5) is satisfied for almost all .
Iii A parametric form of PDEs with integral terms
We divide the parameters of a PDE into three groups depending based on the type of constraint they appear in – namely the continuity constraints, the in-domain dynamics, and the boundary conditions. The continuity constraints specify the existence of partial derivatives and boundary values for each state as required by the in-domain dynamics and boundary conditions. The boundary conditions are represented as a real-valued algebraic constraint that maps the distributed state to a vector of boundary values which are then used to constrain the in-domain dynamics. The in-domain dynamics (or generating equation) specify the time derivative of the state at every point in the interior of the domain and allow for both integral and spatial derivative operators. In the following subsections, we highlight the parameters required to unambiguously define the above three types of constraint in a PDE.
Iii-1 The continuity constraint
Given a PDE state, , the continuity constraint can be uniquely defined by the parameter , which partitions and orders the PDE states by increasing differentiability as follows.
Given such an , we can identify all well-defined partial derivatives of .
Notation: For convenience, we define the vector of all continuous partial derivatives of the PDE state as permitted by the continuity constraint as , the vector of all partial derivatives as and the list of all possible boundary values of as , i.e.,
Additionally, we define some notation given the continuity constraint parameter that will repeatedly appear in the subsequent sections. Let be the number of states in , be the number of -times differentiable states, and be the total number of possible partial derivatives of .
Iii-2 Boundary Conditions
Given an , we now parameterize a generalized class of boundary conditions consisting of a combination of boundary values and integrals of the PDE state. Specifically, the boundary conditions are parameterized by square integrable functions for and as
where is the number of specified boundary conditions. For reasons of well-posedness, as discussed in Section IV, we typically require .
These boundary conditions and the continuity constraints collectively define the domain of the infinitesimal generator – which specifies a set of acceptable solutions for the PDE – as
Notation: For convenience, we collect all the parameters which define the boundary-valued constraint in Eq. (8) and introduce the shorthand notation which represents the labelled tuple of system parameters as
When this shorthand notation is used to denote a given set of system parameters, it is presumed that all parameters have appropriate dimensions.
Iii-3 Dynamics of the PDE
Finally, we may now define the dynamics of the PDE which is parameterized by the functions , and as
with the constraint . In this representation, the kernels represent integral operators.
Notation: For convenience, we collect all the parameters from the generating equation (Eq. (III-3)) and introduce the shorthand notation which represents the labelled tuple of system parameters as
When this shorthand notation is used to denote a given set of system parameters, it is presumed that all parameters have appropriate dimensions. Now, we define the notion of a solution to the above described PDE .
Iv Representing a PDE as a PIE
In Section III, we presented the PDE form Eq. (III-3). Recall that our goal is to find a PIE form of such PDEs so that the computational tools developed for PIEs can be used in analysis/control. Specifically, we will express PDEs introduced in previous section as PIEs of the form
However, first, we have to verify that such a PIE form does indeed exist, which can be verified by testing some additional constraints on . Specifically, we can verify that for any well-posed PDE of the form given in Section III with the initial condition , there exists a corresponding PIE with corresponding initial condition whose solution can be used to construct a solution to the PDE — i.e., there exists a invertible map from any initial condition and corresponding solution of the PDE to a corresponding initial condition and solution of the corresponding PIE. In the following subsection, we propose a sufficient condition for the existence of such an invertible map.
Iv-a Admissibility of the Boundary Conditions
The following definition of admissibility imposes a notion of well-posedness on , the domain of the PDE defined by the continuity constraints and the boundary conditions which when satisfied guarantees the existence of a PIE form for the PDE.
Definition 5 (Admissible Boundary Conditions).
Given an and a parameter set, , we say the pair is admissible if is invertible where
and , , and are as defined in LABEL:fig:Gb_definitions.
Since is invertible only when it is a square matrix, naturally, we require , i.e., when we have differentiable states we need boundary conditions for a well-posed solution. Note that the test for invertibility of depends only on the boundary condition parameters and not the dynamics or the initial condition. Hence, we refer to this test as “well-posedness of the boundary conditions”, which is not necessarily the same as “well-posedness of a PDE”.
Iv-B PIE representation of a PDE: and
Assuming that we have a PDE with admissible , in this subsection, we propose a PIE form for any PDE of the form Eq. (III-3) such that the solutions of the PDE and the corresponding PIE are equivalent, in the sense, the existence of the PDE solution guarantees the existence of a PIE solution and vice versa. The conversion of PDE to PIE and the equivalence of the solutions are summarized below.
Given a set of PDE parameters as defined in Equations 10 and 12 with admissible, let the PI operators be as defined in LABEL:fig:Gb_definitions. Then, for any ( is as defined in Section III-2), satisfies the PDE defined by with initial condition if and only if satisfies the PIE defined by with initial condition where .
The proof is similar to the proof of Theorem 4.4 and 4.5 in  with , , and other unused parameters set to empty sets or zeros. ∎
The above result provides a map from the PDE solution, , to the solution of the corresponding PIE, . However, there also exists an inverse map which is given by that maps the solution of the PIE back to the PDE solution. The invertible relation between the solutions can be summarized as follows.
Given an , and with admissible, let be as defined in LABEL:fig:Gb_definitions, as defined in Eq. (III-2) and diag . Then we have the following.
If , then and .
For any , and .
The proof is similar to the proof of Theorem 4.3 in  with , and other unused parameters set to empty sets or zeros. ∎
This bijective mapping ensures the existence of both solutions when one of them exists. Given the PDE parameters we now have a set of formulae to find the corresponding PIE parameters . Furthermore, we know that a solution to the PDE equation provides a unique solution to the PIE and vice versa.
V Equivalence of representations: PDE and PIE
In this section, we show that the solutions to the two representations (PDE and PIE) have equivalent stability properties and present a solvable optimization problem to prove stability. However, we have to first define a notion of stability for the two representations.
V-a Definitions of stability
We define the notion of stability of an PDE based on the stability of that satisfies the PDE for given initial conditions where the stability is defined with respect to the standard Sobolev norm which is defined as
Definition 6 (Exponential Stability of a PDE).
We say a PDE defined by is exponentially stable, if there exists constants , such that for any , if satisfies the PDE defined by with initial condition then
Similar to the stability of a PDE, we can define the stability of a PIE system based on stability that satisfies the PIE for some initial condition and zero inputs. Unlike, the stability of PDE, the stability of a PIE is defined with respect to the -norm since solutions of PIE need not have spatial continuity.
Definition 7 (Exponential Stability of a PIE).
We say a PIE defined by is exponentially stable, if there exists constants , such that for any , if satisfies the PIE defined by with initial condition , then
Based on these definitions, the goal now is to show that for any that satisfies the PDE and that satisfies the PIE, if with invertible, then exponential decay of implies exponential decay of and vice versa. The main hurdle in proving this is the fact that these two notions of stability are defined using different norms ( and ). For any and such that , we know that implies , but the converse is typically not true. Thus, additional steps are needed to prove the converse implication which is accomplished using a new norm on the space , denoted by , to show that:
is a norm-preserving bijection from to (when equipped with ).
Consequently, is closed under (and is unitary)
is equivalent to (and thus also equivalent to ) on the subspace
Finally, the PDE is exponentially stable if and only if the corresponding PIE defined by is exponentially stable (PIE parameters are related to PDE parameters as shown in LABEL:fig:Gb_definitions).
V-B The map are unitary
First, we would like to show that is complete with respect to the norm . Previously, in Theorem 2, we showed that is invertible. Therefore, we just need to show that preserves the inner product. However, we first define the new -inner product as
Suppose is admissible, is as defined in LABEL:fig:Gb_definitions and inner product is as defined in Equation 15. Then, for any
The proof is similar to the proof of Theorem 6.3a in  with , and other unused parameters set to empty sets or zeros. ∎
Now, we can show that norms induced by the inner products and on are equivalent and consequently, notions of stability with respect to these norms will be equivalent.
Suppose pair is admissible. Then, for any , and there exists a constant such that .
The proof is same as the proof of Theorem 6.4 in  with and other unused parameters set to empty sets or zeros. ∎
Now that we have established that -norm can be upper bounded by -norm, we can use that to prove the equivalence of the stability of PDE and PIE because -norm on is isometric to -norm on .
The proof is same as the proof of Theorem 6.6 in  with and other unused parameters set to empty sets or zeros. ∎
Using the above results, we now propose the following optimization problem to test the stability of a PDE with integral terms.
Given a set of PDE parameters , suppose there exist , , such that , , and where , are as defined in LABEL:fig:Gb_definitions. Then, the PDE defined by is exponentially stable.
Suppose and are as stated above. Suppose solves the PDE defined by for some initial condition . Then solves the PIE defined by , for initial condition .
Let the Lyapunov function candidate be . Then for all . Taking the derivative of with respect to time along the solution trajectories of the PIE, we have
Then, from Gronwall-Bellman inequality,
where and . Since the initial condition was an arbitrary function, the above inequality is satisfied for any and hence the PIE defined by is exponentially stable. Consequently, from Theorem 5, the PDE is exponentially stable. ∎
See , for a parametric form of and that allows the use of Linear Matrix Inequalities to enforce positivity constraint. Once the positivity constraint is rewritten as LMI constraints, we can use an SDP solver to find a feasible solution.
Vi Numerical Example
In this section, we present numerical tests for stability of the two PDEs introduced earlier. The steps involved in setting up the computational problem, specifically, defining the PDE parameters, converting to PIE form, setting up the operator-valued optimization problem, and converting the operator-valued optimization problem to an LMI feasibility test are all performed using the PIETOOLS toolbox for MATLAB.
Vi-a Population Dynamics
Recall the McKendrick PDE model for population dynamics given by
where is the density of population of age at any given time . The boundary condition can be considered as an approximation of the number of newborns at time .
We will select the kernel in the integral term of the boundary condition as , which implies that the population outside some normalized limits do not contribute to the birth of newborns. We will employ a constant mortality rate and vary the to find the mortality rate, below which the population would go extinct (i.e, ).
By testing the stability of the above PDE for various values (using the method of bijection), we determined that for mortality rates greater than the population goes extinct. That is the population will survive ( will not converge to zero) when the population growth rate .
Vi-B Observer-based control of reaction-diffusion equation
Consider the reaction-diffusion example presented in the introduction, where the observer is designed based on boundary measurement. We know that the observer gain, , stabilizes the PDE where
Since we are interested in the practical implementation, we will approximate the gains by a polynomial of a fixed order denoted by . Furthermore, the boundary measurements are replaced by an integral. While a typical approach is to replace point measurements by an integral with a Gaussian kernel centered at the point, in this example specifically, we can use the Fundamental Theorem of Calculus to rewrite the point measurement exactly as
Likewise, we replace the point measurement value by an integral to get the closed-loop observer PDE as
Then, using the conversion formulae presented in the appendix we can find the PIE representation for the closed-loop PDE where replaced by which is the order polynomial approximation of .
For we can prove that -order polynomial approximation (a straight line) is sufficient to guarantee the stability of the closed loop PDE system. However, as becomes large higher order polynomial approximation of ( for and so on) were necessary to prove the stability.
To conclude, we presented a computational method to test the stability of PDEs with integral terms appearing in the boundary conditions and dynamics where the integration is performed with respect to the spatial variable. In the process of achieving this, we presented a standard parametric form for such PDEs, provided a sufficient criterion that guarantees the existence of a PIE representation for such PDEs, and showed that both representations have the same stability properties. Finally, using the equivalence in the stability properties, we formulated the test for stability of PDEs as an optimization problem that can be solved using SDP solvers and demonstrated the application using numerical examples.
-  J. Appell, A. Kalitvin, and P. Zabrejko. Partial Integral operators and integro-differential equations: pure and applied mathematics. CRC Press, 2000.
-  J. H. Cushman, B. X. Hu, and F. Deng. Nonlocal reactive transport with physical and chemical heterogeneity: Localization errors. Water Resources Research, 31(9):2219–2237, 1995.
-  Amritam Das, Sachin Shivakumar, Matthew Peet, and Siep Weiland. Robust analysis of uncertain ODE-PDE systems using PI Multipliers, PIEs and LPIs. In 2020 IEEE Conference on Decision and Control (CDC), pages 634–639, 2020.
-  W. A. Day. Heat conduction within linear thermoelasticity, volume 30. Springer Science & Business Media, 2013.
-  M. Gil. On stability of linear Barbashin type integrodifferential equations. Mathematical Problems in Engineering, 2015.
-  B. L. Keyfitz and N. Keyfitz. The McKendrick partial differential equation and its uses in epidemiology and population study. Mathematical and Computer Modelling, 26(6):1–9, 1997.
-  M. Krstic and A. Smyshlyaev. Boundary control of PDEs: A course on backstepping designs, volume 16. SIAM, 2008.
-  M. Peet. A new state-space representation for coupled PDEs and scalable lyapunov stability analysis in the SOS framework. In 2018 IEEE Conference on Decision and Control (CDC), pages 545–550, 2018.
-  M. Peet. A Partial Integral Equation (PIE) representation of coupled linear PDEs and scalable stability analysis using LMIs. Automatica, 125:109473, 2021.
-  S. Shivakumar, A. Das, S. Weiland, and M. Peet. Duality and -optimal control of coupled ODE-PDE systems. In 2020 IEEE Conference on Decision and Control (CDC), pages 5689–5696, 2020.
-  S. Shivakumar, A. Das, S. Weiland, and M. Peet. Extension of the Partial Integral Equation representation to GPDE Input-Output systems. https://rgdoi.net/10.13140/RG.2.2.35592.70401, 2022.
-  C. Yang, A. Zhang, X. Chen, X. Chen, and J. Qiu. Stability and stabilization of a delayed pide system via spid control. Neural Computing and Applications, 28(12):4139–4145, 2017.