1

A Convex Sum-of-Squares Approach to Analysis,
State Feedback and Output Feedback
Control of Parabolic PDEs

arXiv:1507.05888v3 [cs.SY] 2 Sep 2016

Aditya Gahlawat and Matthew. M. Peet

Abstract—We present an optimization-based framework for
analysis and control of linear parabolic Partial Differential
Equations (PDEs) with spatially varying coefficients without
discretization or numerical approximation. For controller synthesis, we consider both full-state feedback and point observation
(output feedback). The input occurs at the boundary (point
actuation). We use positive definite matrices to parameterize
positive Lyapunov functions and polynomials to parameterize
controller and observer gains. We use duality and an invertible
state variable transformation to convexify the controller synthesis
problem. Finally, we combine our synthesis condition with the
Luenberger observer framework to express the output feedback
controller synthesis problem as a set of LMI/SDP constraints.
We perform an extensive set of numerical experiments to
demonstrate accuracy of the conditions and to prove necessity
of the Lyapunov structures chosen. We provide numerical and
analytical comparisons with alternative approaches to control
including Sturm Liouville theory and backstepping. Finally we
use numerical tests to show that the method retains its accuracy
for alternative boundary conditions.
Index Terms—Distributed parameter systems, partial differential equations (PDEs), control design, sum of squares.

I. I NTRODUCTION
Partial Differential Equations (PDEs) are used to model
quantities which vary in both space and time with early
examples including the D’Alembert wave equation (1746);
the Euler-Bernoulli beam (1750); the Euler equations (1757);
and the Fourier heat equation (1822). Today, the use of PDE
models has expanded to include phenomena such as the magnetohydrodynamics of plasma in a fusion reactor [43], tumour
growth, infectious diseases, and ecological succession [25,
Chapter 11]. However, despite the variety of phenomena
modeled by PDEs, compared to the literature on Ordinary
Differential Equations (ODEs), our knowledge of how to
analyze and control PDEs remains incomplete.
Consider the following class of scalar-valued anisotropic
parabolic PDEs with input u(t) ∈ R,
wt (x, t) = a(x)wxx (x, t) + b(x)wx (x, t) + c(x)w(x, t), (1)
This research was supported by the Chateaubriand program and NSF
CAREER Grant CMMI-1151018.
Aditya Gahlawat is with the Department of Mechanical, Materials and
Aerospace Engineering at the Illinois Institute of Technology, Chicago, IL,
60616 USA e-mail: (agahlawa@hawk.iit.edu).
Matthew. M. Peet is with the School of Engineering of Matter, Transport
and Energy at Arizona State University, Tempe, AZ, 85287-6106 USA e-mail:
(mpeet@asu.edu).
Color versions of one or more of the figures in this paper are available
online at http://ieeexplore.ieee.org.
Digital Object Identifier 10.1109/TAC.2016.2593638

x ∈ [0, 1], t ≥ 0, which has output v(t) = w(1, t) ∈ R and
mixed boundary conditions of the form
w(0, t) = 0,

wx (1, t) = u(t),

(2)

where a, b and c are polynomials with a(x) ≥ α > 0, for
x ∈ [0, 1]. We assume the controller is parameterized
by scalar
R1
R1 and function R2 as u(t) = R1 ŵ(1, t)+ 0 R2 (x)ŵ(x, t)dx
where ŵ is an estimate of w obtained from some set of
observer dynamics. The objective of the paper is to propose
an optimization-based method for determining controller gains
R1 and R2 and observer dynamics which minimize certain
closed-loop gains.
Control of PDE models is a challenging problem in that
slight variations in the type of PDE, boundary conditions, etc.
may dramatically alter properties of the solution [22]. The
model defined above is classified as an anisotropic parabolic
PDE with point inputs and point outputs. The term anisotropic
means that the values of the coefficients a(x), b(x) and
c(x) depend on the spatial variable x ∈ [0, 1]. Examples
of anisotropic systems include heat conduction with nonhomogeneous conductive properties or a wave propagating
through a medium of varying density. The term point input
(boundary actuated) means that the control input determines
one of the boundary values and therefore has no direct measurable effect on Equation (1). This is in contrast to the case
of distributed inputs, wherein the control effort is spread over
some measurable subset of the domain. In a similar manner,
the term point output means that the sensor measures the state
at a single point in the domain and hence the output operator
is unbounded in the L2 induced norm.
Perhaps the most common approach to analysis and control
of PDEs is based on the use of discrete approximation. Such
approximation techniques typically use a model reduction
wherein the PDE is approximated by a set of ODEs. Finitedimensional linear control theory is then used to analyze
stability and design control laws for the finite-dimensional
approximations [24], [23]. Furthermore, results have been
obtained which show that as the order of the discrete approximation increases, stability of the closed-loop approximations
will eventually imply stability of the closed-loop PDE. A
disadvantage of the discrete approach, however, is that the
required order of the approximation cannot be established a
priori. Consequently, the stability of any particular approximation is not guaranteed to imply stability of the actual
PDE. For this reason, among others, there has been some
interest in finding approaches to analysis and control which

2

can be applied directly to the PDE model without the use of
discretization or numerical approximation. Such methods are
sometimes termed direct or infinite-dimensional.
There has been significant progress in the use of direct
methods for control of PDE systems. One approach is to
express the control problem as the solution to a set of
operator-valued Riccati equations. This approach was applied
to distributed input/distributed output optimal control problems
in [41]. The problem of point actuation with full-state feedback
was considered in [21] (and related work) and extended in [20]
to output feedback controller synthesis through the use of a
Luenberger observer. An alternative Riccati-based approach
for static output feedback of a certain class of well-posed
operators can be found in [35], [36], [42]. A limitation of these
Riccati-based methods, however, is that they rely on finitedimensional numerical methods for obtaining the operatorvalued solution. While convergence of these approximations
has been demonstrated [21], for a given level of approximation, it is not possible to determine whether existence of a
solution implies the closed loop is stable when applied to the
original PDE.
Backstepping [18] is a popular and well-developed method
for boundary control of parabolic PDE systems. This approach
is based on the use of a boundary controller to transform
the PDE to a simpler model for which the existence of a
decreasing Lyapunov function has previously been established.
The backstepping approach is commonly used in the literature
and has been extended to many classes of PDE systems - see,
e.g. [17], [33], [34], [32]. A highlight of the backstepping
method is that for certain types of system, stabilizability
guarantees the existence of a backstepping transformation.
However, a drawback of the backstepping approach is that
it is not based on optimization, but rather typically requires
numerical integration of a PDE in order to obtain the stabilizing controller - thereby making extensions to robust and
optimal control more difficult. Although a complete survey
the of the literature on direct control of PDEs is beyond
the scope of this paper, we do note some other significant
results on the use of Lyapunov functions for analysis and
control of infinite dimensional systems including: a rotating
beam [4]; quasilinear hyperbolic systems [3]; and control of
systems governed by conservation laws [5]. As an alternative
to Lyapunov-based methods, a classical spectral approach to
stability and stabilization is based on Sturm-Liouville theory.
In particular, the differential operators which define the PDEs
in this paper can be adapted to the Sturm-Liouville framework,
from whence one can attempt to determine stability and
design static output-feedback controllers. As is demonstrated
in Section XIII, however, the use of dynamic output feedback
offers considerable advantages over this classical framework.
The goal of this paper is to design stabilizing static state
feedback and dynamic output feedback controllers for PDE
systems. Our approach is inspired by the use of Linear Matrix
Inequalities (LMIs) and Semi-Definite Programming (SDP)
in control of ODEs. For stability analysis, as discussed in
Sections V and VI, we use positive definite matrices to create a
linear parametrization of a cone of Lyapunov functions which
are positive on the Hilbert space L2 . Specifically, the Lyapunov

functions have the quadratic form V = hZ(w), P Z(w)iL2
where w ∈ L2 is the infinite-dimensional state, P is a positive
definite matrix and Z is a fixed vector of multiplication and
integral operators with monomial multipliers and kernels. The
derivative of the Lyapunov function is likewise constrained
to be a negative definite quadratic form. If such a Lyapunov
function exists it directly proves stability of the PDE - i.e. there
is no numerical approximation. For state-feedback controller
synthesis, the controller, as defined above, is parameterized by
a scalar R1 and a function R2 . Combining these gains with the
quadratic Lyapunov functions used for stability analysis yields
synthesis conditions which are bilinear in the design variables.
However, as described in Sections VII and VIII, by defining an
invertible state transformation and a variable substitution, we
derive synthesis conditions which are linear in the optimization
variables. Next, in Section IX we introduce a class of infinitedimensional Luenberger observers with observer gains, again
parameterized by the coefficients of polynomials. Again, using
the Lyapunov function from Section V and the invertible state
variable transformation from Section VII, we obtain SDPbased observer synthesis conditions. Finally, in Section XI, we
verify the accuracy of the method with a series of numerical
tests which indicate that the proposed stability conditions are
accurate to several decimal places and suggest that for any
suitably controllable and observable system, the algorithm
will return an observer-based controller. This is followed by
Section XIII, wherein we include numerical and analytical
comparisons with other results in the literature, including
Sturm-Liouville and backstepping.
A significant contribution of the paper, in addition to a
new approach to analysis and control of PDEs, lies in the
flexibility of the optimization-based approach. Specifically, as
the use of LMIs for control of ODEs enabled the field of
robust control, so too does our LMI/Lyapunov-based approach
to control of PDEs allow the extension to analysis and control
of PDEs with parametric uncertainty, PDEs with nonlinearity,
multivariate PDEs and PDEs coupled with ODEs or delays.
Finally, we note that our approach is complementary to several
recent results in the use of LMIs for stability and control of
PDEs, including, e.g. our early work in [27], modeling and
control of nonlinear dynamic systems in [39], stability analysis
of semilinear parabolic and hyperbolic systems in [12] and the
numerous results contained in [26].
II. N OTATION
We denote the vector space of m-by-n real matrices by
Rm×n and the subspace of symmetric matrices by Sn ⊂
Rn×n where the multiplicative and additive identities are
denoted by In ∈ Sn and 0m,n ∈ Rm×n , respectively. For
P ∈ Sn , P > 0 (P ≥ 0) denotes that P is a positive
definite (positive semi-definite) matrix. The spaces of n−times
continuously differentiable and infinitely differentiable functions on an interval W ⊂ R are denoted by C n (W ) and
C ∞ (W ), respectively. In a similar manner, C n,m (W1 , W2 )
represents the space of n and m−times continuously differentiable functions on intervals W1 ⊂ R and W2 ⊂ R,
respectively. The shorthand ux and ut denote the partial

3

derivative of u with respect to independent variables x and
t, respectively. For a bivariate function, f (x, y), we denote
D1 f := fx and D2 f := fy - i.e. D1 is differentiation with
respect to the first variable and D2 is differentiation with
respect to the second. In a similar manner, D12 := fxx and
D22 := fyy . Recall L2 (W ) is the standard Hilbert space
of square Lebesgue integrable functions with standard norm
and inner product. We nuse H n (W ) to denote the Sobolev
o
n
subspace H n (W ) := y ∈ L2 (W ) : ddtny ∈ L2 (W ) with
Pn D dm x dm y E
. We occainner product hx, yiH n =
m=0 dtm , dtm
L2

sionally let L2 (0, 1) := L2 ([0, 1]) and H n (0, 1) := H n ([0, 1]).
For normed spaces X and Y , L(X, Y ) denotes the Banach space of bounded linear operators from X to Y with
induced norm kGkL := supkxkX =1 kGxkY and we denote
L(X) := L(X, X). We define Zd (x) ∈ Rd+1×1 to be the
column vector of all monomials in variables x of degree d
or less arranged in increasing lexicographical order. We often
use the notation Zd (x, y) := Zd ([x; y]) to denote the vector
of monomials in both x and y. For any function T ∈ L2
we use MT : L2 → L2 to denote the multiplier operator
defined by T . i.e. (MT w)(x) = T (x)w(x). For any functions
M, K1 , K2 ∈ C ∞ we define

X{M,K1 ,K2 } w (x)
Z x
Z 1
= M (x)w(x) +
K1 (x, ξ)w(ξ)dξ +
K2 (x, ξ)w(ξ)dξ.
x

0

(3)

III. P ROBLEM S TATEMENT
For the system of Equations (1) - (2), the strict positivity of
a(x) implies that the differential operator defining the PDE
is uniformly elliptic [11, Section 6.1]. This means that w
diffuses from higher density to lower density, a property which
is representative of most physical systems. The choice of
sensor and actuator location is somewhat arbitrary. For the
heat equation, input wx (1, t) = u(t) would represent heat
flow into the rod and the output v(t) = w(1, t) represents the
temperature of the rod at that point. Note that the results of
this paper can be adapted to Dirichlet, Neuman and Robin
boundary conditions with only slight modifications to the
conditions and proofs. These extensions are addressed in
Section XIV.
The goal of this article is to design algorithms which resolve
the following problems:
1) Stability Analysis: Establish global exponential stability
of the trivial solution w ≡ 0 of the autonomous system
u(t) = 0 and determine the exponential rate of decay δ.
2) State feedback control: If the autonomous system is
unstable, construct gains R1 ∈ R and R2 (x) ∈ C ∞ (0, 1)
such that if
Z 1
R2 (x)w(x, t)dx,
(4)
u(t) = R1 w(1, t) +
0

then the trivial solution w ≡ 0 is globally exponentially
stable with some desired rate of decay, µ.
3) Output feedback control: If only output feedback is
available (v(t) = w(1, t)), construct gains L1 ∈ C ∞ (0, 1)

and L2 ∈ R such that for stabilizing gains R1 and R2 , if
Z 1
R2 (x)ŵ(x, t)dx,
u(t) = R1 ŵ(1, t) +
0

where ŵ satisfies

ŵt (x, t) =a(x)ŵxx (x, t) + b(x)ŵx (x, t)
+ c(x)ŵ(x, t) + L1 (x) (v̂(t) − v(t)) ,
(5)
for v(t) = w(1, t) and v̂(t) = ŵ(1, t) with boundary
conditions
ŵ(0, t) = 0,

ŵx (1, t) = u(t) + L2 (v̂(t) − v(t)) ,
(6)
then the trivial solution w ≡ 0 of Equations (1) - (2) is
globally exponentially stable.
Note that if we consider only bounded linear operators, then
the structure of the controller in (4) is not restrictive, as any
bounded linear functional can be represented in this way using
only the integral form (second term). However, we also would
like to consider unbounded operators and hence we include
the term R1 w(1, t) as well. If controllers of this form prove
inadequate, then one can
R 1 generalize the structure further to
include terms such as 0 R3 (x)wx (x, t)dx as in [15].
The choice for the structure of the Luenberger observer was
similarly determined in an ad-hoc manner through inclusion of
terms necessary to achieve separation of controller synthesis
and observer design objectives. That is, the goal of the
observer is to stabilize the dynamics of the estimation error
e = ŵ − w and the terms in Equations (5) - (6) were chosen
as the minimal necessary to achieve this objective. Again,
this structure mirrors the structure of observers found in the
backstepping approach.
A. Existence and Uniqueness
We now briefly discuss the uniqueness and existence of
solutions. Define the operator
d
d2
+ b(x)
+ c(x).
dx2
dx
It is known that the operator A restricted to space
A = a(x)

D0 = {w ∈ H 2 (0, 1) :

w(0) = wx (1) = 0},

(7)

(8)

generates a strongly-continuous semigroup, or a C0 semigroup, on L2 (0, 1) (see, e.g., [7, Section 2.1]). More
precisely, one can represent A as the negative of a SturmLiouville operator on D0 and hence, using the spectral properties of a Sturm-Liouville operator, it can be proven that A
restricted to D0 generates a C0 -semigroup on L2 (0, 1) [9].
Thus, using Theorems 3.1.3 and 3.1.7 in [7] we conclude that
in the autonomous case (u(t) = 0), for any initial condition
w0 ∈ D0 there exists a unique classical solution of (1) - (2).
For the state-feedback case, using a fixed point argument
similar to the one presented in [2] it can be shown that for
R1 ∈ R and R2 ∈ L∞ (0, 1), the closed loop system (1) - (2)
with
Z 1
u(t) = R1 w(1, t) +
R2 (x)w(x, t)dx,
0

4

admits a unique local in time solution w ∈ C 1,2 ((0, T ), [0, 1]),
for T > 0 sufficiently small, for any initial condition w0 ∈ D,
where
D = {w ∈ H 2 (0, 1) : w(0) = 0 and
wx (1) = R1 w(1) +

Z

1

R2 (x)w(x)dx}.
0

(9)

Thus if we can establish that any solution of the closed
loop system decays exponentially, then this implies the existence and uniqueness of a unique classical solution w ∈
C 1,2 ((0, ∞), [0, 1]) for any w0 ∈ D. The proof of this statement has been omitted, but follows the arguments presented
in [2, Section 6].
Finally, consider the observer-based controller as defined in
Equations (1) - (2) and (5) - (6). Define the estimator error as
e = ŵ − w, which is governed by

with polynomial multipliers and semi-separable kernels. Then,
if P > ǫI,
square root and hence
E
D 1it has a 1symmetric
V (w) = P 2 Z(w), P 2 Z(w) ≥ ǫkwk2 . For the time
derivative, we will similarly require V̇ (w(t)) + µV (w(t)) =
− hZ(w(t)), QZ(w(t))i, for some scalar µ > 0 and Q > 0
where here and throughout the paper we denote by V̇ the
d
V (w(t)) for any
function which satisfies V̇ (w(t)) := dt
solution of the associated PDE - i.e. the derivative along
solutions or time-derivative. Existence of such P, Q > 0
implies exponential stability of the system. As was done for
LMIs in finite-dimensional systems, this approach can then
be extended to controller and observer synthesis, as outlined
below.
a) Controller Synthesis: For controller synthesis, again
consider the LMI approach for the finite-dimensional system:
ẋ(t) = Ax(t) + Bu(t),

n
m
et (x, t) = a(x)exx (x, t)+b(x)ex (x, t)+c(x)e(x, t)+L1 (x)e(1, t),where x(t) ∈ R and u(t) ∈ R . For this system, there
exists
a
stabilizing
state
feedback
controller
of the form u(t) =
(10)
Rx(t)
if
and
only
if
there
exists
a
positive
definite matrix P
with boundary conditions
and Y ∈ Rm×n such that
e(0, t) = 0,
ex (1, t) = L2 e(1, t).
(11)
(AP + BY ) + (AP + BY )T < 0.
It has been established in [12, Section 2] that for L1 ∈
−1
C 1 (0, ∞) and L2 ∈ R, Equations (10) - (11), if expo- If this LMI is feasible, then for R = Y P , the Lyapunov
T −1
nentially stable, admit a unique local in time solution e ∈ function V (x) = x P x is positive definite and has time
C 1,2 ((0, T ), [0, 1]), for T > 0 sufficiently small, for any initial derivative
condition e0 ∈ De , where
V̇ = xT (P −1 A + P −1 BR + AT P −1 + (BR)T P −1 )x
2
De = {w ∈ H (0, 1) : w(0) = 0 and wx (1) = L2 w(1)}.
= (P −1 x)T (AP + BRP + P AT + (BRP )T )(P −1 x)
(12)
= y T (AP + BY + P AT + (BY )T )y < 0,
Therefore, if we can establish that any solution of the
coupled closed-loop dynamics decays exponentially, then the where y = P −1 x. The extension of this LMI approach to
local in time solution can be extended to a classical solution PDEs is to search for a positive definite operator P = Z ∗ P Z
e ∈ C 1,2 ((0, ∞), [0, 1]) for any initial condition e0 ∈ De .
for some RP > 0 and operator Y, defined by (Yz)(z) :=
1
Y1 z(1) + 0 Y2 (x)z(x)dx, such that if u = Rw = YP −1 w,
IV. A F RAMEWORK FOR S TABILITY A NALYSIS AND
Z 1
C ONTROL
R2 (x)w(x, t)dx
u(t) = R1 w(1, t) +
0
Our approach is motivated by the use of LMIs for optimal
Z 1
−1
control of finite-dimensional systems. For example, consider
= Y1 (P w)(1, t) +
Y2 (x)(P −1 w)(x, t)dx,
the autonomous finite-dimensional ODE
0

ẋ(t) = Ax(t),
where x(t) ∈ Rn . This ODE is exponentially stable if and
only if there exists a positive definite matrix P ∈ Sn such that
AT P + P A < 0.

Feasibility of this LMI implies that the Lyapunov function
V (x) = xT P x is positive definite and its derivative along
solutions V̇ (x) = xT (AT P + P A)x is negative definite.
For stability of PDEs, our approach is to use positive matrices to define positive quadratic Lyapunov functions, except that instead of V (x) = xT P x, we will use the form
V (w) = hZ(w), P Z(w)i, where Z : L2 → Rp is a vector
of bases for a subspace of linear operators on L2 (similar
to how x = [x1 , · · · , xn ]T is a vector of bases for the
space of linear functions on Rn ). In our case, however, Z
parameterizes a subspace of multiplier and integral operators

then the Lyapunov function V
=
w, P −1 w
=
−1
−1
Z(P w), PZ(P w) satisfies
V̇
(w(t))
+
2µV
(w(t))
=

− Z P −1 w , QZ P −1 w
for some scalar µ > 0 and
Q > 0, which implies the closed-loop system is exponentially
stable. This is detailed in Section VIII.
b) Observer Synthesis: As mentioned previously, for
observer design, we use a Luenberger observer and a separation principle to decouple the error dynamics as defined
in Equations (10) - (11). For a finite-dimensional Luenberger
observer, where the output is v(t) = Cx(t), the estimator
dynamics are defined using the controller gain F and observer
gain L as
x̂˙ = (A + LC)x̂ − Lv(t) + Bu(t).
If u(t) = F x̂(t), then the error dynamics become
ė(t) = (A + LC)e(t).

5

Existence of an observer gain which renders the error dynamics stable is equivalent to the existence of a P > 0 and T such
that
P A + T C + AT P + C T T T < 0.
If this LMI is feasible, then for L = P −1 T , the Lyapunov
function V (e) = eT P e is positive definite and has derivative
V̇ (e) = eT (P A + P LC + AT P + C T LT P )e
= eT (P A + T C + AT P + C T T T )e < 0.
For the infinite-dimensional PDE, we have two observer gains
which we construct as
L1 (x) = P −1 (T1 (x) + T3 (x))

and L2 = P −1 (T2 ),

for some gains T1 , T2 and T3 and where P = Z ⋆ P Z for
some P > 0. We then use the Lyapunov function V (e) =
hZ(e), P Z(e)i and search for a Q > 0 such that V̇ (e) =
− hZ(e), QZ(e)i ≤ −δV (e) < 0, for some δ > 0. This is
detailed in Section IX.
V. S UM - OF -S QUARES LYAPUNOV F UNCTIONS
S EMI -S EPARABLE K ERNELS

WITH

In this Section, we define the map Z and show how this map
is used to construct Lyapunov functions of the form V (w) =
hZ(w), P Z(w)i. This approach is based on prior work, as
described in [28]. Specifically, we define


Z (x)w(x)
R 1 d1
(Zw)(x) = Rx Zd2 (x, ξ)w(ξ)dξ  ,
x
0 Zd2 (x, ξ)w(ξ)dξ

where recall Zd1 (x) and Zd2 (x, ξ) are the vectors of all
monomials of degree d1 and d2 or less, starting with 1.

Theorem 1. Given d1 , d2 ∈ N and ǫ > 0, ǫ ∈ R, let
Z1 (x) = Zd1 (x) and Z2 (x, ξ) = Zd2 (x, ξ), with n = d1 + 1
and m = 12 (d2 +2)(d2 +1) denoting the length of these vectors,
respectively. Suppose that there exists a matrix P ∈ Sn+2m
such that




ǫ
01,n−1
P12 P13 
 P11 − 0n−1,1 0n−1,n−1
 ≥ 0,
P =
T

P12
P22 P23 
T
T
P13
P23
P33
(13)
where Pij is a partition of P such that P11 ∈ Sn , P22 ∈ Sm
and P33 ∈ Sm . Now let
M (x) = Z1 (x)T P11 Z1 (x),
T

(14)
T

K1 (x, ξ) = Z1 (x) P12 Z2 (x, ξ) + Z2 (ξ, x) P31 Z1 (ξ)
Z ξ
Z2 (η, x)T P33 Z2 (η, ξ)dη
+
Z0 x
+
Z2 (η, x)T P32 Z2 (η, ξ)dη
ξ

+

Z

1

Z2 (η, x)T P22 Z2 (η, ξ)dη,

(15)

x

K2 (x, ξ) = K1 (ξ, x).

(16)

Then
V (w) =

Z

1

w(x)M (x)w(x)dx
Z 1Z x
w(x)K1 (x, ξ)w(ξ)dξdx
+
0
Z0 1 Z 1
+
w(x)K2 (x, ξ)w(ξ)dξdx
(17)
0
x
D 1
E
1
= hZ(w), P Z(w)i = P 2 Z(w), P 2 Z(w) ≥ ǫkwk2 .
0

Proof: The proof follows directly from the definition of
Z and the Sum-of-Squares representation of V .
The form of the Lyapunov function defined by Theorem 1 in
Equation (17) is somewhat atypical for the study of parabolic
PDEs.
A more commonly used version would
be V (w) =
R1
R1
w(x)M
(x)w(x)dx
or
even
yet
V
(w)
=
w(x)M
w(x)dx
0
0
for M > 0. Such forms can be obtained as a special case
of Theorem 1 when Pij = 0 for i 6= j 6= 1. However, as
we discuss in Section XII, neglect of the K1 and K2 terms
results in significantly less accurate conditions for stability and
control.
For polynomials M , K1 and K2 , let X{M,K1 ,K2 } be defined as in (3). If M , K1 and K2 satisfy the conditions
of Theorem 1, then V (w) = w, X{M,K1 ,K2 } w ≥ ǫkwk2 ,
which implies the operator X{M,K1 ,K2 } is positive definite
and furthermore, coercive. Moreover, since M , K1 and K2 are
polynomials, the operator is bounded, which implies that there
exists a θ > 0 such that ǫkwk2 ≤ V (w) ≤ θkwk2 . Finally,
the constraint (16) in Theorem 1 implies that the operator
X{M,K1 ,K2 } is self-adjoint.
As discussed in Section IV, Theorem 1 allows us to use
positive matrices to parameterize positive Lyapunov functions
of the Form (17). By expanding these forms, the coefficients
of the polynomials M , K1 and K2 are linear combinations
of the elements of P > 0. Furthermore, if we can express
the derivative V̇ in the Form (17), where the coefficients are
again linear combinations of the elements of P , then we can
enforce negativity of the derivative along the solutions w by
using V̇ (w) = − hZ(w), QZ(w)i to equate these coefficients
to those defined by Q > 0. Constructing the matrices which
relate the elements of P and Q can be automated using
MATLAB toolboxes for polynomial manipulation such as
MULTIPOLY, contained in the package SOSTOOLS [29] and
further developed in our package DELAYTOOLS [28].
For polynomials M , K1 and K2 , we represent the constraint
w, X{M,K1 ,K2 } w = hZ(w), P Z(w)i for some P > 0 as
{M, K1 , K2 } ∈ Ξ{d1 ,d2 ,ǫ} where
Ξ{d1 ,d2 ,ǫ} := {M, K1 , K2 : M, K1 , K2 satisfy

Theorem 1 for d1 , d2 , ǫ}.

The constraint {M, K1 , K2 } ∈ Ξ{d1 ,d2 ,ǫ} is an LMI constraint
in the coefficients of the polynomials M , K1 and K2 and
the unknown matrix P > 0. In this way, the shorthand
{M, K1 , K2 } ∈ Ξ{d1 ,d2 ,ǫ} allows us to define LMI constraints
implicitly.

6

the classical solution w of (18) - (19) satisfies

VI. A T EST FOR S TABILITY
In this section, we use the results of the previous section
to test the existence of a Lyapunov function which establishes
stability of the scalar parabolic PDE defined in Equations (1) (2). Recall the autonomous (u(t) = 0) form of the PDE
wt (x, t) = a(x)wxx (x, t) + b(x)wx (x, t) + c(x)w(x, t),
(18)
w(0, t) = 0,

wx (1, t) = 0.

(19)

The main technical contribution of this section is reformulating
the derivative of the Lyapunov function V in (17) in the form
of Equation (17). This is achieved in the following theorem
wherein we obtain functions M̂ , K̂1 and K̂2 such that
Z 1
V̇ (w) ≤
w(x)M̂ (x)w(x)dx
0
Z 1Z x
+
w(x)K̂1 (x, ξ)w(ξ)dξdx
0
0
Z 1Z 1
+
w(x)K̂2 (x, ξ)w(ξ)dξdx.
0

x

Note that the inequality in this expression is deliberate, i.e.,
certain negative semidefinite terms have been left out of M̂ ,
K̂1 and K̂2 .
Before giving the main theorem, we define the following
linear map, Ωs , which relates functions M , K1 and K2 to an
upper bound on the time-derivative of the Lyapunov function
defined by these functions. Specifically, we say that
{M̂ , K̂1 , K̂2 } := Ωs (M, K1 , K2 ),

(20)

if


∂
∂
M̂ (x) =
a(x)M (x) − b(x)M (x)
∂x ∂x


∂
[a(x) (K1 (x, ξ) − K2 (x, ξ))]
+2
∂x
ξ=x
+ 2M (x)c(x) −

π2
αǫ,
2

(21)

∂
∂
[a(x)K1 (x, ξ)] − b(x)K1 (x, ξ)
K̂1 (x, ξ) =
∂x ∂x


∂ ∂
[a(ξ)K1 (x, ξ)] − b(ξ)K1 (x, ξ)
+
∂ξ ∂ξ
+ (c(x) + c(ξ)) K1 (x, ξ),
(22)


K̂2 (x, ξ) =K̂1 (ξ, x).

(23)

Theorem 2. Suppose that there exist scalars ǫ, δ > 0,
d1 , d2 , dˆ1 , dˆ2 ∈ N and polynomials M , K1 and K2 such that
{M, K1 , K2 } ∈ Ξd1 ,d2 ,ǫ ,

{−M̂ − 2δM, −K̂1 − 2δK1 , −K̂2 − 2δK2 } ∈ Ξd̂1 ,d̂2 ,0 ,

(b(1) − ax (1))K1 (1, x) − a(1)(D1 K1 )(1, x) = 0,
(b(1) − ax (1))M (1) − a(1)Mx (1) ≤ 0,
K2 (0, x) = 0,

where {M̂ , K̂1 , K̂2 } := Ωs (M, K1 , K2 ). Then for any initial
condition w(0) ∈ D0 , there exists a scalar γ > 0 such that

kw(t)k ≤ γkw(0)ke−δt ,

t > 0,

where D0 is defined in Equation (8).
Proof: Recall the operator X{M,K1 ,K2 } is as defined
in (3). As discussed in Section III, for any w(0) ∈ D0 the
autonomous system admits a unique classical solution. By
Theorem 1, if {M, K1 , K2 } ∈ Ξd1 ,d2 ,ǫ , then
V (w) = w, X{M,K1 ,K2 } w = hw, Pwi ,
satisfies ǫkwk2 ≤ V (w) ≤ θkwk2 for some θ > 0. The
calculation of the time derivative V̇ and its reformulation is
lengthy. It involves integration by parts, the Wirtinger inequality and the assumption a(x) ≥ α. For this reason, we have
included this proof in the appendix as Lemma 3. Continuing,
by Lemma 3, for any w which satisfies Equations (18) - (19),
D
E
V̇ (w(t)) ≤ w(t), X{M̂ ,K̂1 ,K̂2 } w(t) .

Now, since {−M̂ − 2δM, −K̂1 − 2δK1 , −K̂2 − 2δK2 } ∈
Ξd̂1 ,d̂2 ,0 , we have that X{M̂ ,K̂1 ,K̂2 } ≤ −2δP and thus
d
−V̇ (w) − 2δV (w) ≥ 0. This implies that dt
V (w(t)) ≤
−2δV (w(t)) for all t ≥ 0. Thus, V (w(t)) ≤ V (w(0))e−2δt .
Concluding, we have that
r
θ
−δt
.
kw(t)k ≤ γkw(0)ke , γ =
ǫ

Note that using the arguments in the proof of [7, Theorem 5.1.3], the above result holds for weak/mild solutions
where the initial condition need only satisfy w0 ∈ L2 (0, 1).
To test the conditions of Theorem 2, the variables are
the coefficients of the polynomials M , K1 and K2 . The
coefficients of M̂ , K̂1 and K̂2 are then linear combinations
of these variables. Finally, the constraints ∈ Ξd1 ,d2 ,ǫ are
LMI constraints, as discussed in Section V. Constructing the
matrices which map these coefficients can be automated using
SOSTOOLS or DelayTOOLs. The algorithm used can be
adapted from the algorithm presented for output feedback
controller in Section X. Application of the conditions of
Theorem 2 to several numerical examples can be found in
Section XI.
VII. I NVERSION

AND

S TATE T RANSFORMATION

As discussed in Section IV, for controller synthesis, we
will use a state variable transformation z = P −1 w so that
Zh, P Z(P −1 w) = hh, wi. Define P = X{M,K1 ,K2 } , where
X{M,K1 ,K2 } is as defined in (3). Then P has the form
Z x
Z 1
(Pz)(x)=M (x)z(x)+ K1 (x, ξ)z(ξ)dξ+ K2 (x, ξ)z(ξ)dξ,
0

x

where if {M, K1 , K2 } ∈ Ξ{d1 ,d2 ,ǫ} , the operator is coercive
with hw, Pwi ≥ ǫkwk2 . Operators of this type are a combination of a multiplier operator and two integral operators.
Furthermore, since K1 and K2 are polynomials, there exist
polynomials Fi and Gi such that K1 (x, ξ) = F1 (x)T G1 (ξ)
and K2 (x, ξ) = F2 (x)T G2 (ξ). This implies that the two

7
−1

Integral operators with semiseparable kernels are used to represent the input-output map of well-posed Linear Time-Varying
(LTV) systems, as explored in [16, Section I.4, Theorem 4.1].
These operators have certain properties which make them wellsuited for use in Lyapunov functions. Specifically, they are not
trace-class, which means that their eigenvalues may not be
summable. Moreover, as discussed in [16, Section II.2], since
M (x) ≥ ǫ > 0, P −1 is a bounded linear operator and can
be calculated explicitly, as in the following theorem, which is
adapted from [16, Section II.3, Theorem 3.1].
Theorem 3. Suppose that {M, K1 , K2 } ∈ Ξ{d1 ,d2 ,ǫ} for
some d1 , d2 , ǫ > 0 with K1 (x, ξ) = F (x)T G(ξ) and
K2 (x, ξ) = G(x)T F (ξ). Let P ∈ L(L2 (0, 1)) be defined as
P = X{M,K1 ,K2 } , where X{M,K1 ,K2 } is as defined in (3).
Define




G(x)
B(x) =
, C(x) = F (x)T −G(x)T ,
F (x)
−1

H = [N1 + N2 U (1)] N2 U (1)




I 0
0 0
N1 =
, N2 =
,
0 0
0 I

and U (x) = limn→∞ Un (x), where
Z x
Un+1 (x) = I −
B(ξ)M (ξ)−1 C(ξ)Un (ξ)dξ,

(24)

10

10

−1
m+1 w

( ·)

−2

w( ·) − P P

integral operators
R 1 can be combined into a single integral
of the form 0 K(x, ξ)z(ξ)dξ where K is a kernel of the
semiseparable type. That is, there exist functions Fi and Gi
such that
(
F1 (x)T G1 (ξ), x ≥ ξ
K(x, ξ) =
.
F2 (x)T G2 (ξ), otherwise

−3

10

−4

10

−5

10

1

2

3
m

4

5

−1
Fig. 1: kw − PPn+1
wk as a function of n.

and integration, which can be performed in MATLAB or
Mathematica. In practice, we have found that Un converges
after only a few iterations. To illustrate, in Figure 1 we have
applied this approach to a given {M, K1 , K2 } ∈ Ξ1,1,1 and
−1
plot kw − PPn+1
wk as a function of n for the arbitrarily
−1
chosen function w(x) = x(x−0.4)(x−1). Here Pn+1
denotes
−1
the construction for P
defined in Theorem 3 with U (x)
replaced by Un+1 (x). In this case, n = 5 yields an L2 norm
error of ≈ 10−5 . In this example, we approximated M (ξ)−1
using the first five terms of its Chebyshev series.
Finally, we emphasize that construction of P −1 is not part
of the optimization algorithm, but rather is performed after the
algorithm has solved the controller synthesis problem (to be
defined in the following section) and returned the polynomial
variables M , K1 , and K2 .

0

VIII. S TATE -F EEDBACK C ONTROLLER S YNTHESIS
and U1 = I. Then, the inverse of the operator P is given by
Z x
Our
approach to controller synthesis is based on the use

P −1 w (x) =M(x)w(x) +
K1 (x, ξ)w(ξ)dξ
of a state variable transformation y = P −1 w which, by
0
Theorem 3, is guaranteed to exist for any P = X{M,K1 ,K2 }
Z 1
+
K2 (x, ξ)w(ξ)dξ, defined by {M, K1 , K2 } ∈ Ξ{d1 ,d2 ,ǫ} . Specifically, we will
x
use the Lyapunov function V (w) = P −1 w, w = hy, Pyi.
M(x) =M (x)−1 ,
Ignoring the input for the moment and using the operator A
K1 (x, ξ) =M (x)−1 C(x)U (x)(H − I)U (ξ)−1 B(ξ)M (ξ)−1 , defined in Equation (7), the time-derivative of this function
yields the dual stability condition
K2 (x, ξ) =M (x)−1 C(x)U (x)HU (ξ)−1 B(ξ)M (ξ)−1 .
V̇ (w) = 2 P −1 w, Aw = 2 hy, APyi ≤ 0,
Note that since M (ξ) ≥ ǫ, M (ξ)−1 is bounded and
continuous and hence the matrix of rational functions which we must enforce for all y ∈ L2 . The critical point
B(ξ)M (ξ)−1 C(ξ) is bounded and continuous. Therefore, it is that the operator P −1 does not appear explicitly in the
follows from [8, Chapter 3] that the uniform limit U (x) stability condition. Rather its existence is only inferred from
exists and is non-singular for x ∈ [0, 1]. Since U (x) is non- the constraint on P that {M, K1 , K2 } ∈ Ξ{d1 ,d2 ,ǫ} . The next
singular on [0, 1], the matrix H is well defined. Therefore, by step in our approach is to combine this dual stability condition
construction M, K1 , K2 ∈ C ∞ . Furthermore, note that since with a variable substitution through the use of a controller of
P satisfies ǫkwk2 ≤ hw, Pwi ≤ θkwk2 for some θ > 0, then the form
Z 1
1/θkwk2 ≤ w, P −1 w ≤ 1/ǫkwk2 .
−1
u(t) = Y1 (P w)(1, t) +
Y2 (x)(P −1 w)(x, t)dx
Theorem 3 not only proves existence, but gives a practi0
cal method for constructing the state variable transformation
Z 1
P −1 for which Zh, P Z(P −1 w) = hh, wi. Specifically, if
R2 (x)w(x, t)dx,
= R1 w(1, t) +
0
we truncate the sequence Un and approximate M (x)−1 by
a Chebyshev series, then construction of the functions M, wherein we have replaced the original controller gains R1 and
K1 and K2 is simply a matter of polynomial multiplication R2 with the new variables Y1 and Y2 . Once Y1 and Y2 are

8

determined by the SOS solver, the actual gains R1 and R2
can be recovered by computing P −1 and applying the formula
listed here.
Before giving the main theorem, we recall that the input
enters the dynamics as
wt (x, t) = a(x)wxx (x, t) + b(x)wx (x, t) + c(x)w(x, t),
(25)
w(0, t) = 0,

wx (1, t) = u(t).

(26)

The goal, then, is to define conditions on P (which defines
M , K1 and K2 ) as well as on Y1 and the polynomial Y2 such
that the closed-loop system is exponentially stable.
To simplify exposition, we now define the following linear
map, Ωc , which relates functions M , K1 and K2 to an upper
bound on the time-derivative of the Lyapunov function defined
by these functions for the controller dynamics. Specifically, we
say that
{M̂ , K̂1 , K̂2 } := Ωc (M, K1 , K2 ),
(27)
if
M̂ (x) = (axx (x) − bx (x)) M (x) + b(x)Mx (x)
π2
αǫ
+ a(x)Mxx (x) + 2c(x)M (x) −
2 

∂
+ a(x) 2
,
[K1 (x, ξ) − K2 (x, ξ)]
∂x
ξ=x
(28)
K̂1 (x, ξ) = a(x)(D12 K1 )(x, ξ) + b(x)(D1 K1 )(x, ξ)
+ a(ξ)(D22 K1 )(x, ξ) + b(ξ)(D2 K1 )(x, ξ)
+ (c(x) + c(ξ)) K1 (x, ξ),
(29)
K̂2 (x, ξ) = K̂1 (ξ, x).

(30)

Theorem 4. Suppose that there exist scalars ǫ, µ > 0,
d1 , d2 , dˆ1 , dˆ2 ∈ N and polynomials M , K1 and K2 such that
{M, K1 , K2 } ∈ Ξd1 ,d2 ,ǫ

and K2 (0, x) = 0.

Further suppose
{−M̂ − 2µM, −K̂1 − 2µK1 , −K̂2 − 2µK2 } ∈ Ξd̂1 ,d̂2 ,0 ,
where {M̂ , K̂1 , K̂2 } = Ωc (M, K1 , K2 ). Let
Y1 <

Mx (1) ax (1) − b(1)
+
M (1),
2
2a(1)

Y2 (x) = (D1 K1 )(1, x).

(31)
If the control input u(t) is defined as
Z 1
−1
Y2 (x)(P −1 w)(x, t)dx
u(t) =Y1 (P w)(1, t) +
0
Z 1
=R1 w(1, t) +
R2 (x)w(x, t)dx,
(32)
0

where P −1 is as defined for P = X{M,K1 ,K2 } in Theorem 3
and X{M,K1 ,K2 } is as defined in (3), then there exists a scalar
γ > 0 such that for any initial condition w(0) ∈ D (where
D is as in Equation (9)) the solution w of (25) - (26) exists,
belongs to C 1,2 ((0, ∞), [0, 1]) and satisfies
kw(t)k ≤ γkw(0)ke−µt ,

t > 0.

Proof: We start the proof by observing that since
{M, K1 , K2 } ∈ Ξd1 ,d2 ,ǫ , as per Theorem 1, these polynomials define a positive operator P = X{M,K1 ,K2 } such that
ǫkwk2 ≤ hw, Pwi ≤ θkwk2 for some θ > 0. Furthermore,
by Theorem 3, there exist bounded and continuously differentiable functions M, K1 and K2 such that P −1 = X{M,K ,K }
1
2
satisfying 1/θkwk2 ≤ w, P −1 w ≤ 1/ǫkwk2. We now
propose the Lyapunov function
V (w) = P −1 w, w = P −1 w, PP −1 w
Z 1
=
(P −1 w)(x)M (x)(P −1 w)(x)dx
0
Z 1Z x
+
(P −1 w)(x)K1 (x, ξ)(P −1 w)(ξ)dξdx
0
0
Z 1Z 1
+
(P −1 w)(x)K2 (x, ξ)(P −1 w)(ξ)dξdx.
x

0

−1

Let y = P w. Note that if w ∈ H 2 (0, 1), then y = P −1 w ∈
H 2 (0, 1). Now, since 1/θkwk2 ≤ w, P −1 w ≤ 1/ǫkwk2, we
have that the Lyapunov function is upper and lower bounded.
Now suppose that
Z 1
−1
Y2 (x)(P −1 w)(x, t)dx
u(t) = Y1 (P w)(1, t) +
0
Z 1
= R1 w(1, t) +
R2 (x)w(x, t)dx.
0

∞

Since M, K1 , K2 ∈ C (0, 1) and Y2 is polynomial, we have
that R2 ∈ C ∞ (0, 1). Therefore, as discussed in Section III, the
closed loop System (25) - (26) admits a solution w ∈ H 2 (0, 1)
which implies y = P −1 w ∈ H 2 (0, 1). Again, the calculation
of the time derivative V̇ and its reformulation is lengthy. It
involves integration by parts, the Wirtinger inequality and
the assumption a(x) ≥ α. This proof is in the appendix as
Lemma 4 which establishes that for any w which satisfies
Equations (25) - (26),
D
E
V̇ (w(t)) ≤ y(t), X{M̂ ,K̂1 ,K̂2 } y(t)

+ y(1, t)N y(1, t) + 2y(1, t)a(1)M (1)yx(1, t),

where N = a(1)Mx (1) + (b(1) − ax (1))M (1). Now,
since {−M̂ − 2µM, −K̂1 − 2µK1 , −K̂2 − 2µK2 } ∈
Ξd1 ,d2 ,0 , we have that X{M̂,K̂1 ,K̂2 } ≤ −2µP and
D
E
hence y(t), X{M̂ ,K̂1 ,K̂2 } y(t)
≤ −2µ hy(t), Py(t)i =
−1
−2µ w(t), P w(t) . Applying this to the inequality, we get
V̇ (w(t)) ≤ − 2µ w(t), P −1 w(t)

+ y(1, t)N y(1, t) + 2y(1, t)a(1)M (1)yx(1, t).

A sufficient condition for stability, then, is that
2y(1)a(1)M (1)yx(1)
≤
−y(1)N y(1). Unfortunately,
however, our control input enters via wx (1) and not yx (1).
To see the relationship between wx (1) and yx (1), we expand
the former and then solve for the latter as follows
wx (1, t) =Mx (1)y(1, t) + M (1)yx (1, t)
Z 1
+
(D1 K1 )(1, x)y(x, t)dx,
0

(33)

9

Z

where solving for M (1)yx (1) yields
M (1)yx (1, t) =wx (1, t) − Mx (1)y(1, t)
Z 1
−
(D1 K1 )(1, x)y(x, t)dx.

(34)

0

0

0

Now, examining the proposed controller, we obtain
Z 1
wx (1, t) = u(t) = R1 w(1, t) +
R2 (x)w(x, t)dx
0
Z 1
= R1 (Py)(1, t) +
R2 (x)(Py)(x, t)dx
0
Z 1
Y2 (x)y(x, t)dx,
= Y1 y(1, t) +
0

which is expressed in the new optimization
R 1 variables Y1 and
Y2 . Now, plugging wx (1) = Y1 y(1) + 0 Y2 (x)y(x)dx into
the time-derivative of the Lyapunov function, we get
V̇ (w(t)) ≤ − 2µV (w(t))

+ y(1, t)2 (N + 2a(1)Y1 − 2a(1)Mx (1))
Z 1
+ 2y(1, t)a(1)
Y2 (x)y(x, t)dx
0
Z 1
(D1 K1 )(1, x)y(x, t)dx. (35)
− 2y(1, t)a(1)
0

By inspection, we see that the stability conditions are now
N + 2a(1)Y1 − 2a(1)Mx (1) < 0
and Y2 (x) = (D1 K1 )(1, x). This then implies that V̇ (w(t)) ≤
−2µV (w(t)) for all t ≥ 0 and hence V (w(t)) ≤
V (w(0))e−2µt . Since kwk2 ≤ θV (w), we have
p
p
kw(t)k ≤ θV (w(t)) ≤ θV (w(0))e−µt
p
≤ θ/ǫkw(0)ke−µt .
At this point, it is significant to note that given values for the
variables Y1 , Y2 , M , K1 and K2 , the controller gains R1 and
R2 can be found by calculating M, K1 and K2 via Theorem
3 and using the formula
Z 1
R1 w(1) +
R2 (x)w(x)dx
0
Z 1
= Y1 y(1) +
Y2 (x)y(x)dx
0
Z 1
= Y1 (P −1 w)(1) +
Y2 (x)(P −1 w)(x)dx
0
Z 1
= Y1 M(1)w(1) +
Y1 K1 (1, x)w(x)dx
Z x 0

Z 1
Z 1
+
Y2 (x)
K2 (x, ξ)w(ξ)dξ dx
K1 (x, ξ)w(ξ)dξdx+
0

0

x

x

0

so that

This implies that the Lyapunov function satisfies
V̇ (w(t)) ≤ − 2µV (w(t)) + y(1, t)N y(1, t)
+ 2y(1, t)a(1) (wx (1, t) − Mx (1)y(1, t))
Z 1
− 2y(1, t)a(1)
(D1 K1 )(1, x)y(x, t)dx.

1

= Y1 M(1)w(1) +
Y1 K1 (1, x)w(x)dx
0

Z 1 Z 1
Z x
Y2 (ξ)K1 (ξ, x)dξ+
+
Y2 (ξ)K2 (ξ, x)dξ w(x)dx,
R1 =Y1 M(1),

(36)
Z 1
Z x
R2 (x)=Y1 K1 (1, x)+ Y2 (ξ)K1 (ξ, x)dξ+ Y2 (ξ)K2 (ξ, x)dξ,
x

0

(37)

where we have used the identity
Z 1Z 1
Z 1Z ξ
Z 1Z x
f (x, ξ)dξdx =
f (x, ξ)dxdξ =
f (ξ, x)dξdx,
0

x

0

0

0

0

and the fact that K1 (x, ξ) = K2 (ξ, x).
IX. O BSERVER S YNTHESIS
In Section VIII, we described LMI conditions under which
one can obtain controller
R 1 gains R1 and R2 (x) such that input
u(t) = R1 w(1, t) + 0 R2 (x)w(x, t)dx ensures exponentially
stability of the closed-loop system. However, this form of
controller requires measurement of the state w(x, t) at every
point x ∈ [0, 1] at all times. Implementation of such a
controller is problematic as such a distributed measurement
is unlikely to be readily available. A more common scenario
is one in which we may only measure the value of w(x, t)
at discrete points in the domain. In particular, we assume
that only a single measurement is available at the boundary
of the domain, and in particular, at v(t) = w(1, t). As
discussed in Section III, in this scenario, we seek to find
an estimator/observer which will yield a real-time estimate
of the state of the system at all points and which, if used in
closed-loop, will ensure exponential stability of the closedloop. Specifically, our observer is a dynamic system with
input v(t) = w(1, t) and output ŵ(x, t), where ŵ(x, t) is
the estimate of the state at time t. We adopt the Luenberger
observer framework discussed previously, which implies that
the dynamics of the observer are given by
ŵt (x, t) = a(x)ŵxx (x, t) + b(x)ŵx (x, t)
+ c(x)ŵ(x, t) + L1 (x) (v̂(t) − v(t)) ,
ŵ(0, t) = 0,

(38)

ŵx (1, t) = u(t) + L2 (v̂(t) − v(t)) , (39)

where v̂(t) = ŵ(1, t) is the predicted output and the scalar L2
and function L1 (x) are gains which map error in this predicted
output to the dynamics of the observer state. In the following
theorem, we seek conditions on L1 and L2 which ensure that
if R1 and R2 are as defined in Theorem 4 and the controller
is defined as
Z 1
R2 (x)ŵ(x, t)dx,
(40)
u(t) = R1 ŵ(1, t) +
0

then Equations (38) - (39) coupled with Equations (25) - (26)
and Equation (40) define an exponentially stable system.
Our approach is based on the separation principle [7,
Chapter 5], [18, Chapter 5]. Specifically, we consider the error
dynamics of the PDE coupled with the observer dynamics

10

in Equations (38) - (39). That is, if we define the error as
e = ŵ − w, then this quantity satisfies
et (x, t) = a(x)exx (x, t) + b(x)ex (x, t) + c(x)e(x, t)
+ z1 (x, t),
(41)
e(0, t) = 0,

ex (1, t) = z2 (t),

(42)

where the feedback signals z1 and z2 are defined as
z1 (x, t) := L1 (x)e(1, t)

and z2 (t) := L2 e(1, t).

(43)

The key point is that the error dynamics do not depend on
the choice of controller gains R1 and R2 . In the following
theorem, this will allow us to choose observer gains L1 and
L2 which stabilize the error dynamics. Then, in Theorem 6
we will show that if the controller gains are chosen as
per Theorem 4 and the observer gains are chosen as per
Theorem 5, then the coupled dynamics are stable in both the
state and state estimate. Unlike for controller synthesis, the
conditions for stabilization of the error dynamics are based
on the use of a simple Lyapunov function V (e) = he, Pei
where the operator P = X{M,K1 ,K2 } is defined by some
{M, K1 , K2 } ∈ Ξd1 ,d2 ,ǫ .
The following theorem is motivated by the LMI approach
as defined in Section IV, wherein as before the variables M ,
K1 and K2 are defined by a positive definite matrix P and the
observer variables are scalar T2 and polynomials T1 and T3
(defined by their vector of coefficients). Referring to the LMI
motivation, these observer variables are similar to the matrix
T and the observer gains are then recovered as L2 = P −1 T2
and L1 = P −1 (T1 + T3 ), which is similar to the LMI observer
gain matrix L = P −1 T .
Theorem 5. Suppose there exist scalars ǫ, δ > 0,
d1 , d2 , dˆ1 , dˆ2 ∈ N and polynomials M , K1 and K2 such that
{M, K1 , K2 } ∈ Ξd1 ,d2 ,ǫ ,

and K2 (0, x) = 0.

Further suppose
{−M̂ − 2δM, −K̂1 − 2δK1 , −K̂2 − 2δK2 } ∈ Ξd̂1 ,d̂2 ,0 ,
where {M̂ , K̂1 , K̂2 } := Ωs (M, K1 , K2 ). Let M, K1 and K2
−1
define X{M,K
= X{M,K ,K } as in Theorem 3 and
1 ,K2 }
1

2

−1

(44)
L2 :=(a(1)M (1)) T2 ,
L1 (x) :=M(x)(T1 (x) + T3 (x))
Z x
+
K1 (x, ξ)(T1 (ξ) + T3 (ξ))dξ
0
Z 1
+
K2 (x, ξ)(T1 (ξ) + T3 (ξ))dξ, (45)
x

where
T1 (x) = −0.5((b(1) − ax (1))K1 (1, x) − a(1)(D1 K1 )(1, x)),
(46)
T2 < −0.5((b(1) − ax (1))M (1) − a(1)Mx (1))
T3 (x) = −L2 a(1)K1 (1, x),

(47)
(48)

and X{M,K1 ,K2 } and X{M,K ,K } are as defined in (3). Then
1
2
for any e which satisfies (41) - (42) with initial condition

e(0) ∈ De (See Equation (12)), there exists a scalar γ > 0
such that
ke(t)k ≤ γke(0)ke−δt , t > 0.
Proof: We start by observing that since {M, K1 , K2 } ∈
Ξd1 ,d2 ,ǫ , as per Theorem 1, these polynomials define a positive
operator P = X{M,K1 ,K2 } such that ǫkwk2 ≤ hw, Pwi ≤
θkwk2 for some θ > 0. Furthermore, by Theorem 3, there exist
bounded and continuously differentiable functions M, K1 and
K2 which define the positive operator P −1 = X{M,K ,K } .
1
2
Therefore, since L1 ∈ C ∞ (0, 1), we have that the closed-loop
error dynamics (41) - (42) admit a local in time solution e for
any e0 ∈ De .
We now propose the Lyapunov function
V (e) = he, Pei
Z 1
=
e(x)M (x)e(x)dx
0
Z 1Z x
+
e(x)K1 (x, ξ)e(ξ)dξdx
0
0
Z 1Z 1
+
e(x)K2 (x, ξ)e(ξ)dξdx.
0

x

The derivative of this Lyapunov function is identical to the
one in Theorem 2 except for the presence of the terms z1 and
z2 defined in (43). Specifically, we have
D
E
V̇ (e) ≤ e, X{M̂,K̂1 ,K̂2 } e
+ 2 hPe, z1 i + 2z2 a(1)(Pe)(1)
+ 2 he, MR1 e(1)i + e(1)R2 e(1),

where R1 (x) = (b(1) − ax (1))K1 (1, x) − a(1)(D1 K1 )(1, x)
and R2 = (b(1) − ax (1))M (1) − a(1)Mx (1). In the proof of
Theorem 2, we had z1 = 0 and z2 = 0 and hence the stability
condition was that R1 = R2 = 0 and that X{M̂,K̂1 ,K̂2 } ≤
−2δP. For the observer, we similarly require X{M̂ ,K̂1 ,K̂2 } ≤
−2δP. However, we now have the observer gains z1 (x) =
L1 (x)e(1) and ex (1) = z2 = L2 e(1) which the algorithm
can choose in order to cancel out R1 and R2 . Unfortunately,
however, these gains depend on M and K1 and the gains
are currently bilinear with the operator variable P (and the
functions M K1 , and K2 which define it). Hence we would
like to perform a variable substitution. This is complicated,
however, by the fact that there are two observer gains - one at
the boundary and one directly injected into the dynamics. Let
us first examine the second gain z2 = L2 e(1) which appears
in the term
z2 a(1)(Pe)(1) = e(1)L2 a(1)(Pe)(1)
= e(1) L2 a(1)M (1) e(1)
{z
}
|
T2
1

+

Z

e(1)L2 a(1)K1 (1, x)e(x)dx
Z 1
= e(1)T2 e(1) +
e(1)L2 a(1)K1 (1, x)e(x)dx,
0

0

where we have made the variable substitution T2 =
L2 a(1)M (1) which implies T2 is a scalar variable. The

11

variable L2 is thereby partially eliminated from the search.
However, since a(x) > 0 and M (x) > 0, given T2 , the gain L2
can later be recovered as L2 = (a(1)M (1))−1 T2 . Of course,
this variable substitution has not completely eliminated the
original variable L2 . To completely eliminate L2 will require
assistance from the second gain L1 . To see how this is done,
we examine the second term in which z1 appears
+
*

hPe, z1 i= e, PML1 e(1) =he, MT1 e(1)i+he, MT3 e(1)i .
| {z }

Theorem 6. Given d1 , d2 , dˆ1 , dˆ2 ∈ N and ǫ, δ, µ > 0, suppose
that there exist polynomials M , N , K1 , K2 , S1 and S2 such
that
{M, K1 , K2 } ∈ Ξd1 ,d2 ,ǫ ,
(49)
{−M̂ − 2µM, −K̂1 − 2µK1 , −K̂2 − 2µK2 } ∈ Ξd̂1 ,d̂2 ,0 ,
(50)
K2 (0, x) = 0,
(51)

MT1 +MT3

Here we have defined a new variable
R xT1 (x) which is defined by T1 (x) := M (x)L1 (x) + 0 K1 (x, ξ)L1 (ξ)dξ +
R1
x K2 (x, ξ)L1 (ξ)dξ − T3 (x) for which MT1 c = PML1 c −
MT3 c for any c ∈ R where T3 will be defined shortly.
Furthermore, for any T3 , the map L1 7→ T1 is invertible with
L1 (x) :=M(x)(T1 (x) + T3 (x))
Z x
+
K1 (x, ξ)(T1 (ξ) + T3 (ξ))dξ
0
Z 1
K2 (x, ξ)(T1 (ξ) + T3 (ξ))dξ
+
x

if P −1 = X{M,K ,K } . In this way, we eliminate the variable
1
2
L1 and replace it with T1 and T3 . The next step, then, is to
choose T3 so as to cancel the remaining term which contains
L2 . This is done using he, MT3 e(1)i, which we expand to get
Z 1
e(x)T3 (x)e(1)dx,
he, MT3 e(1)i =
0

which
we
would
like
to
use
to
eliminate
R1
e(1)L
a(1)K
(1,
x)e(x)dx.
Clearly,
then,
the
appropriate
2
1
0
choice for T3 is
T3 (x) = −L2 a(1)K1 (1, x).
Note that the dependence of T3 on L2 is admissible because
T3 is not a free variable and L2 is computed directly from T2 .
This means that once feasible values for T1 and T2 have been
found, we then calculate L2 from T2 , then use L2 to calculate
T3 and then use T1 and T3 to calculate the gain L1 .
Concluding the proof, the time-derivative of the Lyapunov
function becomes
D
E
V̇ (e) ≤ e, X{M̂ ,K̂1 ,K̂2 } e
+ 2 he, (MT1 + MR1 ) e(1)i + e(1) (2T2 + R2 ) e(1).

Therefore, if T1 = −R1 , 2T2 + R2 < 0 and {−M̂ −
2δM, −K̂1 − 2δK1 , −K̂2 − 2δK2 } ∈ Ξd̂1 ,d̂2 ,0 , we have that
V̇ (e) ≤ −2δV (e),
which, in a similar manner as Theorem 2 establishes exponential stability of the error dynamics with decay rate δ.
X. A N LMI CONDITION FOR O UTPUT-F EEDBACK
S TABILIZATION
In this section we briefly summarize the results of the
paper by giving an LMI formulation of the output-feedback
controller synthesis problem.

{N, S1 , S2 } ∈ Ξd1 ,d2 ,ǫ ,
(52)
{−N̂ − 2δN, −Ŝ1 − 2δS1 , −Ŝ2 − 2δS2 } ∈ Ξd̂1 ,d̂2 ,0 ,
(53)
S2 (0, x) = 0,

(54)

where {M̂ , K̂1 , K̂2 } = Ωc (M, K1 , K2 ), {N̂ , Ŝ1 , Ŝ2 } =
Ωs (N, S1 , S2 ) and 2dˆ1 and 2dˆ2 + 1 are the degrees of M ,
N and K1 , S1 , respectively.
Then, there exist gains R1 , R2 (x), L1 (x) and L2 such that
if
Z 1
u(t) = R1 ŵ(1, t) +
R2 (x)ŵ(x, t)dx,
(55)
0

and w satisfies Equations (25) - (26) and ŵ satisfies Equations (38) - (39) with a zero initial condition then
kw(t)k ≤ γkw(0)ke−κt ,
for some γ > 0 and any κ satisfying 0 < κ < min{µ, δ}.
Proof: If the conditions in (49) - (51) are satisfied,
then the polynomials M , K1 and K2 satisfy the constraints
of Theorem 4. Therefore, we may construct R1 and R2 (x)
using (36) - (37). Similarly, if N , S1 and S2 satisfy (52) - (54),
then the conditions of Theorem 5 are satisfied with M = N ,
K1 = S1 and K2 = S2 . Thus, we can construct observer gains
L1 (x) and L2 using (44) - (45). Now, let Pc = X{M,K1 ,K2 } ,
P̂c = X{M̂ ,K̂1 ,K̂2 } , Po = X{N,S1 ,S2 } and P̂o = X{N̂ ,Ŝ1 ,Ŝ2 } .
Therefore, the theorem conditions imply that P̂c ≤ −2µPc
and P̂o ≤ −2δPo .
Using the proof of Theorem 5, there exists a scalar β1 > 0
such that
V̇o (e) ≤ −2δVo (e) − β1 e(1)2 ,
(56)
where Vo (e) = he, Po ei. Similarly, for the observer dynamics
in (38) - (39) with the input (55), using the proof of Theorem 4
one can prove that there exists a scalar β2 > 0 such that
V̇c (ŵ) ≤ − 2µVc (ŵ) + 2 hŷ, L1 e(1)i + ŷ(1)(2a(1)L2 )e(1)
− β2 ŷ(1)2 ,

(57)

where ŷ = Pc−1 ŵ and Vc (ŵ) = ŵ, Pc−1 ŵ = hŷ, Pc ŷi.
From (56) - (57) we infer that for any r > 0 we have


* ŷ 
ŷ +
rV̇o (e) + V̇c (ŵ) ≤ −2rδVo (e) + ŷ(1) , U ŷ(1) ,
e(1)
e(1)
(58)

12

where

and Mupad. To help with understanding this process, we define
several subroutines which perform specific relevant tasks and
combine them in the pseudo code which would be used to
obtain the observer-based controllers.
and the inner product is defined on L2 (0, 1) × R × R. Now, [M,K1 ,K2 ]=mult_semisep(ǫ)
for any 0 < κ < min{δ, µ}, if we choose r > 0 sufficiently
• Declares polynomial variables M , K1 and K2 and
large, it follows that U ≤ diag(−2κPc , 0, 0). Thus, from (59)
enforces the constraint {M, K1 , K2 } ∈ Ξd1 ,d2 ,ǫ .
we get that
[M̂,K̂1 ,K̂2 ]=omega_primal(M,K1 ,K2 )
• Constructs M̂ , K̂1 and K̂2 as defined by the map Ωs
rV̇o (e) + V̇c (ŵ) ≤ − 2rδVo (e) − 2κVc (ŵ)
in (20).
≤ − 2κ (rVo (e) + Vc (ŵ)) .
[M̂,K̂1 ,K̂2 ]=omega_dual(M,K1 ,K2 )
Therefore defining V (ŵ, e) = rVo (e) + Vc (ŵ), we get that
• Constructs M̂ , K̂1 and K̂2 as defined by the map Ωc
in (27).
V̇ (ŵ, e) ≤ −2κV (ŵ, e).
eq_constr(F)
Integrating in time,
• Given a set F of univariate/bivariate polynomials,
r he, Po ei + ŵ, Pc−1 ŵ ≤ e−2κt r hw(0), Po w(0)i , (59)
declares element wise equality constraint F = 0.
[M,K
1 ,K2 ]=inv_op(M,K1 ,K2 )
where we have used the fact that ŵ(0) = 0 and thus e(0) =
• Given {M, K1 , K2 } = Ξd1 ,d2 ,ǫ calculates the inverse
−w(0). Now, as discussed, there exist scalars θ1 , θ2 > 0 such
multiplier M and kernels K1 and K2 by approximating
that
U
(x) by performing the integration in (24) a finite
1
1
ǫkek2 ≤ he, Po ei ≤ θ1 kek2 ,
kŵk2 ≤ ŵ, Pc−1 ŵ ≤ kŵk2 .
number of times and using a Chebyshev series approxθ2
ǫ
imation of M (x)−1 .
Therefore, using (59) we get
[R1 ,R2 ]=controller_gains(M,K1 ,K1 ,M,K1 ,K2 )
rθ1
2
2
2 −2κt
• The function defines Y1 and Y2 (x) using (31). Consekek + kŵk ≤
kw(0)k e
,
σ
quently, R1 and R2 (x) are defined using (36) and (37),
respectively.
where σ = min(rǫ, 1/θ2 ). Thus,
r
[L1 ,L2 ]=observer_gains(M,K1 ,K1 ,M,K1 ,K2 )
rθ1
kek, kŵk ≤
kw(0)ke−κt .
• The function constructs T2 using (47) and sets L2 usσ
ing (44). Then the function constructs T1 (x) and T3 (x)
Finally, using the fact that kwk ≤ kŵk + kek produces
using (46) and (48) and constructs L1 (x) using (45).
r
A pseudo code for the SOSTOOLS implementation of the SDP
rθ1
kw(t)k ≤ 2
kw(0)ke−κt .
is presented in Algorithm 1.
σ

−2µPc
U = ⋆
⋆

0
−β2
⋆


L1
a(1)L2  ,
−rβ1

The variables in Theorem 6 are polynomials which are
parameterized by vectors of coefficients associated to a predetermined monomial basis. There are two types of constraints
on these variables: equality constraints between polynomials;
and constraints of the form ∈ Ξd1 ,d2 ,ǫ . To test the conditions of Theorem 6, these variables and constraints must be
converted to a form recognized by an SDP solver such as
SeDuMi [37]. Many of these tasks have already been automated in SOSTOOLS [29] and our extended toolbox, DelayTOOLS [28]. Specifically, SOSTOOLS has functionality for
declaring polynomial variables and enforcing scalar equality
constraints. Furthermore, DelayTOOLS [28] allows the user
to declare matrix-valued equality constraints and create new
polynomial variables which satisfy ∈ Ξd1 ,d2 ,ǫ . Furthermore,
the multipoly toolbox allows one to manipulate polynomial
variables in order to construct new dependent polynomials
such as in {M̂ , K̂1 , K̂2 } = Ωc (M, K1 , K2 ). Once all variables
and constraints have been declared, SOSTOOLS converts all
constraints and variables to a format which can be accepted
by SDP solvers such as SeDuMi, SDPT3 or MOSEK. The aposteriori polynomial manipulations such as operator inversion
can be performed using a combination of the multipoly toolbox

XI. N UMERICAL R ESULTS
In this section we test the conditions of Theorems 2, 4 and 5
by applying them to two parameterized instances of scalar
parabolic PDEs. The first instance is a variation of the classical
isotropic heat equation. Because this system is well-studied,
we are able to compare our results with a number of existing
results in the literature. The second system is an anisotropic
PDE with arbitrarily chosen coefficients. Both instances have
an instability term, parameterized by an instability factor,
λ. For both systems, we test stability, find controllers and
construct observer-based controllers.
Example 1: Our first system is defined as follows.
wt (x, t) = wxx (x, t) + λw(x, t),

λ ∈ R,

(60)

with boundary conditions
w(0, t) = 0,

wx (1, t) = u(t).

The output of the PDE is v(t) = w(1, t). For u(t) = 0, the
analytical solution of this PDE is given by
w(x, t) =

∞
X

n=1

eλn t hw0 , φn i φn (x),

13

d=7
8
9
10
11
Algorithm 1: Output-feedback controller synthesis.
λ = 14.3982
17.9626
22.8645
23.3093
27.1179
Declaring optimization variables:
1) [M,K1 ,K2 ]=mult_semisep(ǫ)
TABLE II: Max. λ as a function of d1 = d2 = d for which
2) [N,S1 ,S2 ]=mult_semisep(ǫ)
the conditions of Theorem 4 are feasible, thereby implying
Constructing polynomials:
the existence of an exp. stabilizing state-feedback controller
1) [M̂,K̂1 ,K̂2 ]=omega_dual(M,K1 ,K2 )
for PDE (60).
2) [N̂,Ŝ1 ,Ŝ2 ]=omega_primal(N,S1 ,S2 )
Declaring affine constraints:
d=7
8
9
10
11
λ = 14.5233
17.7643
23.4406
24.7772
27.8820
1) eq_constr((-M̂-2µM,-K̂1 -2µK1 ,-K̂2 -2µK2 )
-mult_semisep(0))
TABLE III: Max. λ as a function of d1 = d2 = d for which
2) eq_constr((-N̂-2δN,-Ŝ1 -2δS1 ,-Ŝ2 -2δS2 )
the
conditions of Theorem 5 are feasible, thereby implying the
-mult_semisep(0))
existence
of an exp. stabilizing output-feedback controller for
3) eq_constr(K2 (0,x))
PDE
(60).
4) eq_constr(S2 (0,x))
if SOS problem is feasible then
Return outputs: M, K1 , K2 , N, S1 , S2 .
a dynamic output feedback controller can be constructed for
Calculating control gains:
any value of λ > 0.
1) [M,K1 ,K2 ]=inv_op(M,K1 ,K2 )
Example 2: To illustrate the versatility of the proposed
2)
method,
we next consider the following arbitrarily chosen
[R1 ,R2 ]=controller_gains(M,K1 ,K1 ,M,K1 ,K2 )
anisotropic
system
Calculating observer gains:
1) [N,S1 ,S2 ]=inv_op(N,S1 ,S2 )
wt (x, t) = a(x)wxx (x, t) + b(x)wx (x, t) + c(x)w(x, t), (61)
2) [L1 ,L2 ]=observer_gains(N,S1 ,S1 ,N,S1 ,S2 )
where a(x) = x3 − x2 + 2, b(x) = 3x2 − 2x and c(x) =
−0.5x3 + 1.3x2 − 1.5x + 0.7 + λ with λ ∈ R. Although the
d=3
4
5
6
7
analytic
analytical solution to this PDE is not readily available, we
λ = 0.59
2.19
2.457
2.46
2.461
2.467
may use a finite-difference scheme to numerically simulate
the system and thereby estimate the range of λ for which
TABLE I: Max. λ as a function of d1 = d2 = d for which the
the PDE (61) is stable. Specifically, we find that the system
exp. stability conditions of Theorem 2 are feasible, implying
is unstable for λ > 4.66. To determine the accuracy of the
stability of PDE (60) with u(t) = 0.
conditions of Theorem 2, we find the largest λ for which the
conditions of Theorem 2 are feasible. Table IV lists the largest
√
such λ using ǫ, δ = 0.001 as a function of polynomial degree
where λn = λ − (2n − 1)2 π 2 /4 and φn = 2 sin((2n −
d1 = d2 = d. The maximum λ for which we can prove the
1)πx/2). This implies that Equation (60) is unstable for λ >
exponential stability for is λ = 4.62, which is 99.14% of
π 2 /4 ≈ 2.467. To test the numerical accuracy of the stability
the predicted stability margin of 4.66. The < 1% discrepancy
conditions in Theorem 2, we found the largest λ > 0 for which
may be due to conservatism or inaccuracy in the predicted
the conditions of Theorem 2 are feasible as a function of the
maximum λ on account of inaccuracy in the discretization or
parameters d1 and d2 which define the degree of the variables
poor choice of initial conditions in the simulation.
M , K1 and K2 . Table I presents these results for ǫ, δ = 0.001.
To test the accuracy of the conditions in Theorem 4, we
For d1 = d2 = 7, we can construct a Lyapunov function which
again find the largest λ for which the conditions of Theorem 4
proves stability for λ = 2.461, which is 99.74% of the stability
are feasible with ǫ = 0.001 and µ = 0.001, thereby implying
2
margin π4 ≈ 2.4674.
the existence of an exponentially stabilizing state-feedback
To test the accuracy of the conditions in Theorem 4, we controller. Table V presents this maximum λ as a function
find the largest λ for which the conditions of Theorem 4 are of the degree d = d = d. The results suggest that for
1
2
feasible with ǫ = 0.001 and µ = 0.001, thereby implying sufficiently high degree, a static state-feedback controller can
the existence of an exponentially stabilizing state-feedback be constructed for any value of λ > 0.
controller. Table II presents this maximum λ as a function
To test the accuracy of the conditions of Theorem 5, we
of the degree d1 = d2 = d. The results suggest that for again find the largest λ for which the conditions of Theorem 5
sufficiently high degree, a static state-feedback controller can are feasible with ǫ = 0.001 and δ = 0.001, thereby implying
be constructed for any value of λ > 0.
To test the accuracy of the conditions of Theorem 5, we
d=3
4
5
6
7
simulation
find the largest λ for which the conditions of Theorem 5 are
λ = 4.37
4.61
4.61
4.62
4.62
4.66
feasible with ǫ = 0.001 and δ = 0.001, thereby implying
the existence of an exponentially stabilizing dynamic outputfeedback controller with output v(t) = w(1, t). Table III TABLE IV: Max. λ as a function of d1 = d2 = d for which the
presents this maximum λ as a function of the degree d1 = exp. stability conditions of Theorem 2 are feasible, implying
d2 = d. The results suggest that for sufficiently high degree, stability of PDE (61) with u(t) = 0.

14

d=4
λ = 19.0216

5
36.1359

6
39.7247

7
43.5974

8
44.5219

TABLE V: Max. λ as a function of d1 = d2 = d for which
the conditions of Theorem 4 are feasible, thereby implying
the existence of an exp. stabilizing state-feedback controller
for PDE (61).
d=4
λ = 18.3090

5
36.0199

6
38.0478

7
40.5931

8
44.079

TABLE VI: Max. λ as a function of d1 = d2 = d for which
the conditions of Theorem 5 are feasible, thereby implying the
existence of an exp. stabilizing output-feedback controller for
PDE (61).

the existence of an exponentially stabilizing dynamic outputfeedback controller with output v(t) = w(1, t). Table VI
presents this maximum λ as a function of the degree d1 =
d2 = d. The results suggest that for sufficiently high degree,
a dynamic output feedback controller can be constructed for
any value of λ > 0.
We conclude with the conjecture that the proposed method
is asymptotically accurate in the sense that, for any λ > 0,
if the PDE (1) - (2) is stable in the autonomous sense, then
the conditions of Theorem 2 will be feasible for sufficiently
high d1 and d2 . Moreover, we conjecture that if the system
is observable and controllable for some suitable definition
of controllability and observability, then the conditions of
Theorems 4 and 5 will be feasible for sufficiently high d1 and
d2 . We emphasize, however, that this is only a conjecture and
additional work must be done in order to make this statement
rigorous and determine its veracity. A further caveat to these
results is the observation that the maximum degree d1 and d2
for which the conditions can be tested is a function of the
memory and processing speed of the computational platform
on which the experiments are performed. Specifically, the
number of optimization variables in the underlying SDP problem is determined by the number of polynomial coefficients
which scales as O(d2 ). To illustrate, all numerical experiments
presented in this paper were performed on a machine with
8 gigabytes of random access memory, which limited our
analysis to a maximum degree of d1 = d2 = 11 for PDE (60)
and d1 = d2 = 8 for PDE (61).
In the following subsection, we illustrate the controllers and
observers which result from feasibility of the conditions of
Theorems 4 and 5 using numerical simulation.
A. Numerical Implementation of Observer-Based Controllers
To illustrate the observer-based controllers which result
from feasibility of the conditions of Theorems 4 and 5, we
take the anisotropic PDE (61) with λ = 35. This value of λ
renders the autonomous system unstable. We then synthesize
controller and observer gains using the results of Theorems 4
and 5 for d1 = d2 = 6, along with the inverse state
transformation defined in Theorem 3. For the inverse state
transformation, M (x)−1 is approximated using a sixth order

Ex. 1
Ex. 2

d=1
λ = 3.91
λ = 3.51

2
4.78
7.03

3
4.88
8.59

4...9
4.88
8.59

10
4.88
8.59

K1 , K2 6= 0
27.1179
44.5219

TABLE VII: Re-evaluation of the results of Tables II and V
with added constraint K1 = K2 = 0.
Ex. 1
Ex. 2

d=1
λ = 3.89
λ = 3.51

2
4.79
7.12

3
4.88
8.43

4...9
4.88
8.43

10
4.88
8.43

K1 , K2 6= 0
27.8820
44.079

TABLE VIII: Re-evaluation of the results of Tables III and VI
with added constraint K1 = K2 = 0.
Chebyshev series approximation and 5 iterations are used to
define U∞ ∼
= U5 . The controllers are then applied to the state
and estimator dynamics, which are then discretized using a
trapezoidal approximation. The initial state is set to
2

w0 (x) = e

− (x−0.3)
2(0.07)2

2

−e

− (x−0.7)
2(0.07)2

,

while the initial observer state is set to ŵ(x, 0) = 0. Figures 2(a) - 2(c) illustrate the state evolution of the system,
observer and the control effort respectively. Finally, Figure 3
illustrates the integral control gain R2 (x). Note that its behavior at the boundaries is logical since at x = 0, the boundary
condition w(0, t) = 0 ensures that no control effort is required.
Whereas, at x = 1, the control exerts maximum effort.
XII. N ECESSITY

OF S EMI -S EPARABLE K ERNELS IN THE
LYAPUNOV F UNCTION
Recall that the Lyapunov functions used in Theorems 2, 4,
and 5 all have the form
Z 1
V (w) =
w(x)M (x)w(x)dx
0
Z x

Z 1
Z 1
+
w(x)
K1 (x, ξ)w(ξ)dξ +
K2 (x, ξ)w(ξ)dξ dx.
0

0

x

As mentioned previously, this form is atypical in the study
of parabolic PDEs and the reader may question the necessity of the terms K1 and K2 as their presence significantly
complicates the analysis and increases the complexity of the
stability conditions. Therefore, to illustrate the necessity of
including these terms, in this section we repeat the numerical
examples presented previously with the added restriction that
K1 = K2 = 0 (which translates to Pij = 0 for i 6= j 6= 1 in
Theorem 1). Table VII illustrates these results for the controller
synthesis conditions of Theorem 4 using the same methodology as described in the previous section. These numerical
tests indicate that while inclusion of K1 and K2 allows us
to control the PDE for any λ > 0, when K1 = K2 = 0,
our method will fail for some λ, regardless of the polynomial
degree d1 = d1 = d. As indicated in Table VIII, the results are
similar for the observer synthesis conditions of Theorem 5.
XIII. C OMPARISON W ITH AND R ELATION TO E XISTING
R ESULTS
In this section, we compare our numerical results with
several results in the literature which can be used for stability

15

15

2

2

0

0

10

−2

−4
1

u(t)

ŵ(x, t)

w(x, t)

5

−2

−5
−10

−4
1
1

time

0.5

1

0.5

time

x

0 0

0.5

−15

0.5
0 0

(a) Evolution of closed-loop state w(x, t).

0

x

−20
0

0.2

0.4

time

0.6

0.8

1

(c) Control input wx (1, t) = u(t).

(b) Evolution of closed-loop state estimate ŵ(x, t).

Fig. 2: Evolution of closed loop system for Example 2 with λ = 35 using controller from Theorem 4 and observer from
Theorem 5.

analysis and control, including those based on Sturm-Liouville
theory and backstepping.

q(x) ≤ q1

A. Static Controllers Using Sturm Liouville Theory
The output feedback controllers we construct are dynamic
in that they rely on an auxiliary set of estimator dynamics
which must be simulated in real-time. By contrast, static output
feedback controllers do not use an estimator and instead rely
only on a gain of the form, e.g. u(t) = −κv(t) = −κw(1, t).
Unfortunately, even for finite-dimensional systems the problem
of static output feedback design is unsolved when B 6= I.
That is, there is no LMI or polynomial-time algorithm which
is guaranteed to find a stabilizing output feedback controller
if one exists [38], [13]. However, there are numerous results
which give sufficient conditions for the existence of such a
controller, often based on the use of a fixed Lyapunov function.
For the parabolic PDE which we consider, Sturm-Liouville
theory [10, Chapter 2] can be used to express conditions for
existence of static-output feedback controllers. Specfically, for
u(t) = −κw(1, t), the stability of (1) - (2) depends on the
first eigenvalue of the following Sturm-Liouville eigenvalue
problem


dw(x)
d
p(x)
+ q(x)w(x) = µσ(x)w(x),
(62)
dx
dx
where µ is the eigenvalue and
p(x) = e

R

b(ξ)
dξ
a(ξ)

,

q(x) = c(x)

p(x)
,
a(x)

σ(x) =

p(x)
.
a(x)

The boundary conditions for this eigenvalue problem are
w(0) = 0 and wx (1) + κw(1) = 0. For our system, using
the properties of the coefficients a(x), b(x) and c(x) it can
be established that p is continuously differentiable, q and σ
are continuous and there exist scalars p0 and σ0 such that
p(x) ≥ p0 > 0 and σ(x) ≥ σ0 > 0. If µ1 is the first eigenvalue of (62), then it can be established using the Rayleigh
cc
quotient that µ1 ≤ µcc
1 , where µ1 is the first eigenvalue of
the following constant coefficient Sturm-Liouville eigenvalue
problem
p0

d2 w(x)
+ q1 w(x) = µcc σ1 w(x),
w(x)

subject to the boundary conditions w(0) = 0 and wx (1) +
κw(1) = 0 and where q1 and σ1 are scalars such that

(63)

and σ(x) ≤ σ1 .

Now let us first consider Numerical Example 1, as defined
in Equation (60) in Section XI. In this case, we have that
p0 = 1, q1 = λ and σ1 = 1. Therefore, estimating the first
2
eigenvalue of (63) we get that µcc
1 ≈ λ−π . Since, for stability
cc
we require µ1 < 0, for a large enough κ > 0, a control
input of the form u(t) = −κw(1, t) can stabilize (60) for
λ < π 2 . This result is significantly more conservative than
the results described in Tables II-III which yield a stabilizing
controller for at least λ < 27.1179. Of course this result is not
particularly surprising, as static output feedback controllers are
a subset of dynamic output feedback controllers.
Similarly, for Numerical Example 2 (Equation (61)) we have
p(x) = x3 − x2 + 2, q(x) = −0.5x3 + 1.3x2 − 1.5x + 0.7 + λ
and σ(x) = 1. Thus p0 = 50/27, q1 = 0.7 + λ and σ1 = 1.
Therefore, estimating the first eigenvalue of (63) we get that
cc
µcc
1 ≈ λ − 17.58. As before, we require µ1 < 0. Therefore
for a large enough κ > 0, a control input of the form u(t) =
−κw(1, t) can stabilize (60) for λ < 17.58. Whereas, from
Tables V-VI we see that Theorems 4 and 5 yield a dynamic
output feedback controller for at least λ < 44.079.
B. The Case When A + A⋆ ≤ 0

For some values of the coefficients a(x), b(x) and c(x)
we may have that A + A⋆ ≤ 0 on D0 , where the operator
A is defined in (7) and the set D0 is defined in (8). The
output feedback stabilization of such systems, i.e. systems with
A + A⋆ ≤ 0 and collocated control/observation, is considered
in [6]. The authors in [6] show that for such systems there
exists a scalar κ > 0 (possibly κ = ∞) such that the control
u(t) = −κv(t) exponentially stabilizes the system. We wish
to see if our methodology offers a performance gain over the
controller proposed in [6]. If we choose a(x) = 1, b(x) = 0
and c(x) = π 2 /4, then
π2
d2
+
.
(64)
2
dx
4
Applying integration by parts and Lemma 1, it can be established that A + A⋆ ≤ 0 on D0 . If we apply a controller of
A=

16

d=6
δ = 8.01

7
12.7

8
17.21

9
20.31

10
22.66

11
25.78

TABLE IX: Max. exp. decay rate δ as a function of polynomial
degree, d1 = d2 = d for Equations (1) - (2) with A as in
Equation (64) for which we can construct output feedback
controllers using Theorems 4 and 5.

the form proposed in [6], then u(t) = −κv(t) = −κw(1, t),
for some κ > 0. Using the theory in Subsection XIII-A it is
easily established that even for an arbitrarily large κ > 0, the
closed loop system state will decay with a rate close to, but
less then 3π 2 /4. Whereas, from Table IX we observe that for
d1 = d2 = 11 we can construct an output feedback controller
with a minimum exponential decay rate of 25.78, a significant
improvement over 3π 2 /4.

then for any solution w of Equations (65) - (66),
Z x
E(x, ξ)w(ξ, t)dξ,
z(x, t) = w(x, t) −
0

is a solution of the target system in Equations (67) - (68).
Furthermore, if the map E : w 7→ z is invertible, then stability
of the target system implies stability of the original closedloop PDE. For the example problem given, this E is obtained
as a solution of a kernel-PDE and can be found explicitly
as [18]

p
λ (x2 − ξ 2 )
I1
, 0 ≤ ξ ≤ x ≤ 1, (69)
E(x, ξ) = −λξ p
λ (x2 − ξ 2 )

where I1 is the first order modified Bessel function of the first
kind. Moreover, E has an inverse of the form
Z x

F (x, ξ)z(ξ, t)dξ,
(70)
E −1 z (x) = z(x, t) +
0

20

where

0

J1

R2 (x)

−20

F (x, ξ) = −λξ

−40
−60
−80
−100
0

0.2

0.4

x

0.6

0.8

1

Fig. 3: Control gain R2 (x).

C. Backstepping
Backstepping is a well-known alternative for the construction of stabilizing controllers for parabolic PDEs. Specifically,
the backstepping approach defines a control law which, when
coupled with an invertible state transformation, converts the
controlled parabolic PDE to the form of a desired stable
PDE (the target system). Although backstepping is not an
optimization-based method and does not explicitly search for a
Lyapunov-based stability proof, it turns out that the existence
of a backstepping controller typically implies the existence of
a Lyapunov function of the Form (17), defined by a multiplier
M and semiseparable kernels K1 and K2 . To demonstrate the
existence of this Lyapunov function, let us consider the system
defined by Example 1,
wt (x, t) = wxx (x, t) + λw(x, t),

(65)

w(0, t) = 0,

(66)

wx (1, t) = u(t),

where λ > 0. Now define the target system
zt (x, t) = zxx (x, t),
z(0, t) = 0, zx (1, t) = 0.

(67)
(68)

The key backstepping result is that there exists a function E
such that if
Z 1
(D1 E)(1, x)w(x, t)dx,
u(t) = E(1, 1)w(1, t) +
0

p

λ (x2 − ξ 2 )
p
,
λ (x2 − ξ 2 )

0 ≤ ξ ≤ x ≤ 1, (71)

where J1 is the first order Bessel function of the first kind.
Using properties of Bessel functions, it can be shown that both
kernels E and F are bounded on the domain {(ξ, x) : 0 ≤
ξ ≤ x ≤ 1}. This implies that both E and E −1 are bounded
with induced norms which we denote by kEkL and kE −1 kL .
Now, to understand how this backstepping transformation
implies the existence of a Lyapunov function with semiseparable kernels, we first note that stability of the target
system in Equations (67) - (68) is established using the simple
Lyapunov function
Z 1
Vtarget (z) =
z(x)2 dx = hz, zi ,
0

for which, using (67) - (68), integration by parts and Lemma 1,
we obtain
d
Vtarget (z(t)) ≤ −ǫVtarget (z(t)),
(72)
dt

for any z which satisfies (67) - (68), where ǫ =
implies

π2
2 .

This

ǫ

Vtarget (z(t)) ≤ e−ǫt Vtarget (z(0)) ⇒ kz(x, t)k ≤ e− 2 t kz(x, 0)k.
Now, for the original system we define the Lyapunov
function
Vplant (w) = hEw, Ewi .

(73)

Now, since for any solution, w(t), of the original system, z =
Ew(t) is a solution of the target system, we have that
d
d
Vplant (w(t)) =
hEw(·, t), Ew(·, t)i
dt
dt
d
hz(·, t), z(·, t)i
=
dt
d
= Vtarget (z(t)) ≤ −ǫVtarget (z(t))
dt
= −ǫVtarget (Ew(t)) = −ǫ hEw(t), Ew(t)i

17

= −ǫVplant (w(t)).
Therefore,
Vplant (w(t)) ≤ e−ǫt Vplant (w(0)),
which means
ǫ

kEw(·, t)k ≤ e− 2 t kEw(·, 0)k.

(74)

Boundedness of E and E −1 now implies kw(t)k ≤
kE −1 kL kEw(t)k and kEw(0)k ≤ kEkL kw(0)k, which yields

may be approximated by polynomials, then the existence of a
backstepping controller implies the feasibility of Theorem 4
for some degree. Despite this similarity, there are, of course,
differences between the proposed method and backstepping.
Specifically, our approach is optimization based, whereas the
search for the backstepping transformation is not. Advantages
of the proposed method include the ability to analyze stability
of autonomous PDEs and simple extensions to robust control
of PDEs with parametric uncertainty via Positivstellensatz
results [30].

ǫ

kw(t)k ≤ kE −1 kL kEkL e− 2 t kw(0)k,

which proves that Vplant (w) = kEwk2 establishes exponential
stability of the original system.
We now show that Vplant (w) has a form consistent with
Theorem 4. Expanding
Vplant (w) = hEw, Ewi ,
we get
Vplant (w) =

Z

Z 1Z x
w(x)2 −
w(x)E(x, ξ)w(ξ)dξdx
0
0
0
Z 1Z 1
w(x)E(ξ, x)w(ξ)dξdx
−
0
x
Z 1Z 1Z ξ
+
w(x)E(ξ, x)E(ξ, η)w(η)dηdξdx.
1

0

x

0

Changing the order of integration twice in the last integral and
collecting like terms, we obtain
Z 1
Vplant (w) =
w(x)2 dx
0
Z 1Z x
+
w(x)H1 (x, ξ)w(ξ)dξdx
0
0
Z 1Z 1
w(x)H2 (x, ξ)w(ξ)dξdx,
+
0

x

= w, X{I,H1 ,H2 } w ,

where
H1 (x, ξ) =
H2 (x, ξ) =

Z

1

Zx1
ξ

w(0, t) = 0,

E(η, x)E(η, ξ)dη − E(ξ, x),

ξ

then P

= X{I,H1 ,H2 } and hence
Vplant (w) = P

−1

w, PP

In this subsection we consider the merits of the SOS
approach with respect to finite-dimensional approximation.
That is, we consider whether there are advantages over model
reduction techniques wherein the PDE is reduced to a set of
coupled ODEs - as in, e.g. [1].
Before continuing, we note that establishing a suitable
metric for comparison of finite-dimensional and infinitedimensional approaches is complicated by the fact that that
the methods proposed in this paper are suboptimal. That is, we
are not seeking observer-based controllers which are optimal
in any sense. Rather, we simply seek observer-based controllers which establish closed-loop stability. In this sense, our
methods are roughly equivalent to existing finite-dimensional
approaches in that for all numerical examples considered, we
are able to construct observer-based controllers for suitably
high polynomial degree. In a sense, then, one could argue
that finite-dimensional approaches are superior in that they are
able to go beyond stabilization and construct optimal observerbased controllers using a suitably high level of discretization.
In practice, however, our experience has shown that there are
disadvantages to discretization-based methods such as poleplacement. Specifically, we have seen that if the reduction
scheme is not carefully chosen, discretization may result in
loss of controllability or poorly conditioned controllability
matrices. To illustrate, consider the following model:
wt (x, t) = wxx (x, t) + 15w(x, t),

E(η, x)E(η, ξ)dη − E(x, ξ),

which has the form of a Lyapunov function consistent with
Equation (17) using a semi-separable kernel where we have
M (x) = 1, K1 = H1 and K2 = H2 . In a similar manner, if
we define P = X{I,G1 ,G2 } where
Z 1
F (η, x)F (η, ξ)dη + F (x, ξ),
G1 (x, ξ) =
Zx1
F (η, x)F (η, ξ)dη + F (ξ, x),
G2 (x, ξ) =
−1

D. Finite-Dimensional Approximations

wx (1, t) = u(t).

One approach to reduction of this PDE to a system of ODEs
is to use a finite difference method to approximate the spatial
derivative as


w(x + ∆x2 , t) − w(x, t)
2
wxx (x, t) ≈
∆x1 + ∆x2
∆x2


w(x, t) − w(x − ∆x1 , t)
2
,
−
∆x1 + ∆x2
∆x1
where ∆x1 is the step size to the left of x and ∆x2 is the
step size to the right. Using this scheme we obtain an ODE
model of the form
ẇm (t) = Am wm (t) + B m u(t),

−1

w ,

(75)

which is a form consistent with Theorem 4. Thus we conclude
that for this class of systems, if we assume the function F

(76)

where wm (t), B m ∈ Rm×1 and Am ∈ Rm×m and m ∈ N is
the order of reduction. While relatively straightforward, this
approach creates significant technical challenges. For example:

18

m
cond(C(Am , B m )) ≈

5
107

10
1025
m

20
1063

Dirichlet

m

TABLE X: Condition number of C(A , B ) as a function of
order of reduction m.

Neumann
Robin

a) Controllability of the Reduced Model: The reducedorder model must be chosen so as to maintain the properties
of controllability and observability. In most cases, however,
there is no guarantee that a finite-difference approximation
scheme will preserve these properties. For example, for the
finite difference scheme defined above, it is known that if
the original system is controllable and a uniform grid size is
chosen, then the reduced system is also controllable. However,
if one were to chose a non-uniform grid, then controllability
is no longer guaranteed. For example if one were to chose
a logarithmic grid, for m > 13 the reduced model is not
controllable (although it is still stabilizable). In such a case,
the performance of the closed loop system will be limited by
the location of the uncontrollable eigenvalues.
b) Ill-conditioned Controllability Matrix: Now suppose we
wish to perform pole placement by applying Ackermann’s
formula to the reduced order model. As mentioned, it can be
shown that the reduced order model in (76) is controllable
for any m ∈ N when derived using uniform step sizes
(∆x1 = ∆x2 ) as established by the Hautus test. However,
the pole placement problem (which is similar to our condition
for exponential stabilization with desired decay rate) relies on
inversion of the controllability matrix C(Am , B m ) - a step
which is numerically sensitive to conditioning of C(Am , B m ).
This is problematic since, as seen in Table X, the controllability matrix for this system is ill-conditioned and the condition
number worsens as the level of disretization m increases. This
implies that as the level of discretization increases, numerical
errors may dominate - potentially resulting in unstable or
unpredictable controllers. Naturally, these issues are wellknown and have been addressed in the literature through
methods such as robust place placement [40] or Galerkin
schemes [19]. The advantage of the SOS approach, however, is
that the controllers are provably stable at the pre-lumping stage
and thus the only numerical concern is implementation, which
does not appear to be sensitive to issues such as condition
number.
XIV. A LTERNATIVE B OUNDARY C ONDITIONS
The results of this paper may be readily adapted to other
types of boundary conditions. Specifically, the conditions of
Theorems 2, 4 and 5 can be easily modified to consider
alternative boundary conditions. Although economy of space
prohibits us from presenting these conditions in full, in this
section we give the results of numerical tests performed
using Dirichlet, Neumann and Robin boundary conditions.
Specifically, for the two PDEs (60) and (61) which define
Examples 1 and 2, respectively, in Section XI, we consider
the boundary conditions and the outputs as listed in Table XI.

Output v(t)

Boundary Condition
w(0, t) = 0

wx (1, t)

w(1, t) = u(t)
wx (0, t) = 0
wx (1, t) = u(t)
w(0, t) + wx (0, t) = 0

w(1, t)

w(1, t) + wx (1, t) = u(t)

w(1, t)

TABLE XI: Alternative boundary conditions and outputs for
PDEs (60) and (61).

Dirichlet
Neumann
Robin

d=8
λ = 17.7634
14.8163
13.8367

9
22.8645
17.1814
16.6565

10
23.3093
21.8781
18.6050

11
27.1179
21.8781
18.9758

TABLE XII: Max. λ as a function of polynomial degree,
d1 = d2 = d for PDE (60) with boundary conditions as in
Table XI for which we can construct output-feedback boundary
controllers.

Tables XII and XIII illustrate the maximum λ for which we
can construct output-feedback based controllers as a function
of d1 = d2 = d for PDEs (60) and (61), respectively, for the
boundary conditions listed in Table XI using exponential decay
rates of δ = µ = 0.001. Similar to the observation made in
Section XI, the numerical results in this section suggest that
our methodology is asymptotically accurate for the considered
alternative boundary conditions, that is, given any λ > 0, we
can construct controllers/observers by choosing a large enough
d1 = d2 = d. A more detailed study of alternative boundary
conditions can be found in the thesis work of [14].
XV. C ONCLUSION

AND

F UTURE W ORK

We have defined an algorithmic, polynomial-time approach
to the design of observer-based controllers for a general
class of scalar parabolic partial differential equations using
measurements and feedback at the boundary. The results use
polynomials and semidefinite programming to parameterize
a convex set of positive Lyapunov functions on the Hilbert
space L2 . By combining these Lyapunov functions with an
invertible state transformation, we obtain convex conditions
for stability, controller synthesis and Luenberger observer
design. Furthermore, we have tested our results using parameterized numerical examples in order to show that the
stability conditions are accurate to several significant figures
and the synthesis conditions yield controllers for a large class
of controllable and observable systems. Furthermore, we have

Dirichlet
Neumann
Robin

d=5
λ = 36.0199
29.8492
24.6490

6
38.0478
31.1447
27.8503

7
40.5930
31.1447
27.8503

8
44.079
34.1584
29.4373

TABLE XIII: Max. λ as a function of polynomial degree,
d1 = d2 = d for PDE (61) with boundary conditions as in
Table XI for which we can construct output-feedback boundary
controllers.

19

adapted the approach to three alternative classes of boundary
measurements and actuators. Finally, we have performed a
series of comparisons with existing results in the literature,
showing, e.g. that the method is analytically equivalent to
backstepping for controller synthesis and furthermore is numerically competitive for the examples considered. By using
an optimization-based algorithm defined by polynomials, the
results presented here have the advantage that they may be
further extended to the problem of nonlinear stability analysis,
robust control, and control of coupled, multivariate, hyperbolic
and elliptic PDEs - topics of ongoing research.

is self-adjoint, we can write V̇ (w(t)) = 2 hwt , Pwi. Using
Equation (18) we expand this out to get
V̇ (w(t)) = 2 hwt , Pwi = 2

Γ1 =
Γ2 =

To facilitate presentation in this appendix, we use the
following lemmas. The first is simply a restatement of the
Wirtinger inequality

Γ4 =

Lemma 2. For any bivariate polynomials K and P the
following identity holds for any w ∈ L2 (0, 1)
Z x

Z 1
Z 1
w(x)
K(x, ξ)w(ξ)dξ +
P (x, ξ)w(ξ)dξ dx
0
0
x
Z 1Z x
1
=
w(x) [K(x, ξ) + P (ξ, x)] w(ξ)dξdx
2
0
0
Z 1Z 1
1
+
w(x) [P (x, ξ) + K(ξ, x)] w(ξ)dξdx.
2
0
x
Lemma 3 (Analysis). Given polynomials a, b and c with
a(x) ≥ α > 0, for all x ∈ [0, 1], suppose that there exists
a scalar ǫ > 0 and polynomials M , K1 and K2 such that
{M, K1 , K2 } ∈ Ξd1 ,d2 ,ǫ ,

(b(1) − ax (1))K1 (1, x) − a(1)(D1 K1 )(1, x) = 0,
(b(1) − ax (1))M (1) − a(1)Mx (1) ≤ 0,

K2 (0, x) = 0.
Let

V (w) = w, X{M,K1 ,K2 } w ,
where X{M,K1 ,K2 } is as defined in (3). Then, for any w which
satisfies Equations (18) - (19),
D
E
V̇ (w(t)) ≤ w(t), X{M̂ ,K̂1 ,K̂2 } w(t) ,

where {M̂ , K̂1 , K̂2 } = Ωs {M, K1 , K2 }.

Proof: Let P = X{M,K1 ,K2 } so that V (w) = hw, Pwi.
If w satisfies (18) - (19), then taking the time derivative
of V (w(t)) and since {M, K1 , K2 } ∈ Ξd1 ,d2 ,ǫ implies P

(77)

Z

1

wxx (x, t)a(x)M (x)w(x, t)dx,

0

Γ3 =

The second lemma is accomplished by splitting the integral
in two parts and applying a change in the variable of integration to the second part.

Γn ,

n=1

where

A PPENDIX

Lemma 1 ([31]). Let z ∈ H 2 (0, 1) be a scalar function. Then
Z 1
Z 1
4
zx (x)2 dx.
(z(x) − z(0))2 dx ≤ 2
π 0
0

5
X

Γ5 =

Z

1

wx (x, t)b(x)M (x)w(x, t)dx,

0
2 Z
X

wxx (x, t)a(x)Ki (x, ξ)w(ξ, t)dξdx,

i=1 ∆i
2 Z
X
i=1
Z 1

wx (x, t)b(x)Ki (x, ξ)w(ξ, t)dξdx

∆i

w(x, t)2 M (x)c(x)dx

0

+

2 Z
X
i=1

w(x, t)c(x)Ki (x, ξ)w(ξ, t)dξdx,

∆i

where ∆1 = {(ξ, x) : 0 ≤ ξ ≤ x ≤ 1} and ∆2 = {(ξ, x) :
0 ≤ x ≤ ξ ≤ 1}. Applying integration by parts twice and
using the boundary condition w(0, t) = wx (1, t) = 0 yields
Z 1
Γ1 = −
wx (x, t)2 a(x)M (x)dx
0
Z
1 1 ∂2
[a(x)M (x)] w(x, t)2 dx
+
2 0 ∂x2
1
− (ax (1)M (1) + a(1)Mx (1)) w(1, t)2 .
2
Since a(x) ≥ α > 0 and {M, K1 , K2 } ∈ Ξd1 ,d2 ,ǫ , we have
a(x)M (x) ≥ αǫ. Thus, by application of Lemma 1 we get
Z 1
Z 1
π2
2
−
wx (x, t) a(x)M (x)dx ≤ − αǫ
w(x, t)2 dx.
4
0
0
Therefore, we conclude that
 2

Z
π2
1 1
∂
2
[a(x)M (x)] −
αǫ dx
Γ1 ≤
w(x, t)
2 0
∂x2
2
1
− (ax (1)M (1) + a(1)Mx (1)) w(1, t)2 .
2

(78)

Again, applying integration by parts once and using w(0, t) =
0,
Z
1
∂
1 1
[b(x)M (x)] dx+ b(1)M (1)w(1, t)2 .
w(x, t)2
Γ2 =−
2 0
∂x
2
(79)
Since {M, K1 , K2 } ∈ Ξd1 ,d2 ,ǫ , we have K1 (x, ξ) = K2 (ξ, x)
and thus K1 (x, x) = K2 (x, x). Exploiting this property,
the constraint K2 (0, x) = 0, and the boundary conditions
w(0, t) = wx (1, t) = 0, we apply integration by parts twice
to obtain
 !

Z 1
∂
2
[a(x)(K1 (x, ξ)−K2 (x, ξ))]
dx
Γ3 = w(x, t)
∂x
0
ξ=x

20

+

2 Z
X
i=1

w(x, t)

∆i

− w(1, t)
− w(1, t)

Z

Z

1




∂2
[a(x)K
(x,
ξ)]
w(ξ, t)dξdx
i
∂x2

ax (1)K1 (1, x)w(x, t)dx

where y = P −1 w and {M̂ , K̂1 , K̂2 } ∈ Ωc {M, K1 , K2 }.
Proof: Taking the time derivative of V (w(t)) and since
P −1 is self-adjoint, we obtain
V̇ (w(t)) = 2 wt , P −1 w

0
1

= 2 a(·)wxx + b(·)wx + c(·)w, P −1 w


∂
∂2
= 2 a(·) 2 (Py) + b(·) (Py) + c(·)Py, y
∂x
∂x
5
X
=2
Γn ,
(83)

a(1)(D1 K1 )(1, x)w(x, t)dx.
0

Applying Lemma 2 and using K1 (x, ξ) = K2 (ξ, x), we get
!


Z 1
∂
2
Γ3 =
w(x, t)
[a(x)(K1 (x, ξ) − K2 (x, ξ))]
dx
∂x
0
n=1
ξ=x
Z 1
−1
where y = P w and
−w(1, t) (ax (1)K1 (1, x)+a(1)(D1 K1 )(1, x)) w(x, t)dx
Z 1
0
∂2
"
#
Γ
=
a(x)
(M (x)y(x, t))y(x, t)dx,
T 
1

Z
2
2
∂
X1
∂x2
0
a(x)Ki (x, ξ)
2
∂x
Z
+
w(x, t) ∂ 2
w(ξ, t)dξdx.
1
∂
a(ξ)Ki (x, ξ)
2 ∆i
∂ξ 2
i=1
b(x) (M (x)y(x, t))y(x, t)dx,
Γ2 =
∂x
(80)
0
Z

2 Z 1
X
∂2
Applying integration by parts once and following the same
Ki (x, ξ)y(ξ, t)dξ y(x, t)dx,
Γ3 =
a(x) 2
∂x
procedure as for Γ3 , we get
βi
i=1 0
Z

Z
2

 ∂ T 
Z
2
X 1
X
∂
1
b(x)Ki (x, ξ)
∂x
Γ4 =
b(x)
Ki (x, ξ)y(ξ, t)dξ y(x, t)dx,
w(ξ, t)dξdx
w(x, t) ∂
Γ4 = −
∂x
b(ξ)Ki (x, ξ)
βi
2 ∆i
∂ξ
i=1 0
i=1
Z 1
Z 1
Γ
=
c(x)M (x)y(x, t)2 dx
5
+ w(1, t)
b(1)K1 (1, x)w(x, t)dx.
(81)
0
0
2 Z
X
Finally, employing Lemma 2 produces
y(x, t)c(x)Ki (x, ξ)y(ξ, t)dξdx,
+
Z 1
i=1 ∆i
Γ5 =
w(x, t)2 M (x)c(x)dx
where ∆ = {(ξ, x) : 0 ≤ ξ ≤ x ≤ 1}, ∆ = {(ξ, x) : 0 ≤
1

0

+

Z
2
X
1
i=1

2

w(x, t) ([c(x) + c(ξ)] Ki (x, ξ)) w(ξ, t)dξdx.

∆i

(82)

Finally, we combine the terms (78) - (82) into the derivative (77) and use the constraints
(b(1) − ax (1))K1 (1, x) − a(1)(D1 K1 )(1, x) = 0,

(b(1) − ax (1))M (1) − a(1)Mx (1) ≤ 0,

to eliminate extraneous terms, thereby completing the proof.
Lemma 4 (Controller Synthesis). Given polynomials a, b and
c with a(x) ≥ α > 0, for all x ∈ [0, 1], suppose that there
exists a scalar ǫ > 0 and polynomials M , K1 and K2 such
that
{M, K1 , K2 } ∈ Ξd1 ,d2 ,ǫ ,

K2 (0, x) = 0.

Let
V (w) = w, P −1 w ,
where P = X{M,K1 ,K2 } and X{M,K1 ,K2 } is as defined in (3).
Then, for any w which satisfies (25) - (26)
D
E
V̇ (w(t)) ≤ y(t), X{M̂ ,K̂1 ,K̂2 } y(t)

+ [a(1)Mx (1) + (b(1) − ax (1))M (1)] y(1, t)2
+ 2a(1)M (1)yx (1, t)y(1, t),

2

x ≤ ξ ≤ 1}, β1 = [0, x] and β2 = [x, 1]. Before proceeding
we calculate y(0, t). The definition y = P −1 w implies
Z 1
w(0, t) = M (0)y(0, t) +
K2 (0, x)y(x, t)dx.
0

Therefore, since w(0, t) = 0 and K2 (0, x) = 0, we get
y(0, t) = 0. Now, since M (x)a(x) ≥ αǫ and y(0, t) = 0,
applying integration by parts twice and using Lemma 1
produces

Z 
π2
1 1
axx (x)M (x)+a(x)Mxx (x)− αǫ y(x, t)2 dx
Γ1 ≤
2 0
2

1
+ a(1)Mx (1) − ax (1)M (1) y(1, t)2
2
+ a(1)M (1)yx (1, t)y(1, t). (84)

Similarly, applying integration by parts once yields
Z
1 1
Γ2 =
(b(x)Mx (x) − bx (x)M (x)) y(x, t)2 dx
2 0
1
(85)
+ b(1)M (1)y(1, t)2 .
2
Applying integration by parts twice and Lemma 2 yields
!


Z 1
∂
[K1 (x, ξ) − K2 (x, ξ)]
y(x, t)2 dx
a(x)
Γ3 =
∂x
0
ξ=x
#T 
"

2
Z
2
∂
X1
a(x) ∂x2
Ki (x, ξ)
+
y(x, t)
y(ξ, t)dξdx.
∂2
Ki (x, ξ)
2 ∆i
a(ξ) ∂ξ
2
i=1
(86)

21

In a similar manner as for Γ3 , we obtain


Z
2
∂ T 
X
1
b(x) ∂x
Ki (x, ξ)
Γ4 =
y(ξ, t)dξdx.
y(x, t)
∂
b(ξ) ∂ξ
Ki (x, ξ)
2 ∆i
i=1
(87)
Finally, applying Lemma 2 to Γ5 produces
Z 1
Γ5 =
c(x)M (x)y(x, t)2 dx
0

2

+

1X
2 i=1

Z

y(x, t)(c(x) + c(ξ))Ki (x, ξ)y(ξ, t)dξdx.

∆i

(88)

Substituting Equations (84) - (88) into (83) completes the
proof.
R EFERENCES
[1] M. Balas. Feedback control of linear diffusion processes. International
Journal of Control, 29:523–534, 1979.
[2] A. Balogh and M. Krstic. Stability of partial difference equations
governing control gains in infinite-dimensional backstepping. Systems
and Control Letters, 51:151–164, 2004.
[3] J. M. Coron, G. Bastin, and B. d’Andrea-Novel. Dissipative boundary
conditions for one-dimensional nonlinear hyperbolic systems. SIAM
Journal on Control and Optimization, 47:1460–1498, 2008.
[4] J. M. Coron and B. d’Andrea-Novel. Stabilization of a rotating body
beam without damping. IEEE Transactions on Automatic Control,
43:608–618, 1998.
[5] J. M. Coron, B. d’Andrea-Novel, and G. Bastin. A strict Lyapunov
function for boundary control of hyperbolic systems of conservation
laws. IEEE Transactions on Automatic Control, 52:2–11, 2007.
[6] R. Curtain and G. Weiss. Exponential stabilization of well-posed systems
by colocated feedback. SIAM Journal on Control and Optimization,
45:273–297, 2006.
[7] R. F. Curtain and H. J. Zwart. An introduction to infinite-dimensional
linear systems theory. Springer, 1995.
[8] J. L. Daleckii and M. J. Krejn. Stability of solutions of differential
equations in Banach space. American Mathematical Society, 2002.
[9] C. Delattre, D. Dochain, and J. Winkin. Sturm-Liouville systems are
Riesz-spectral systems. International Journal of Applied Mathematics
and Computer Science, 13:481–484, 2003.
[10] Y. Egorov and V. Kondratiev. On spectral theory of elliptic operators,
Volume 89 of Operator Theory: Advances and Applications. Birkhäuser
Verlag Basel, 1996.
[11] L. C. Evans. Partial Differential Equations. American Mathematical
Society, 1998.
[12] E. Fridman and Y. Orlov. An LMI approach to H∞ boundary control
of semilinear parabolic and hyperbolic systems. Automatica, 45:2060–
2066, 2009.
[13] M. Fu. Pole placement via static output feedback is NP-hard. IEEE
Transactions on Automatic Control, 49:855–857, 2004.
[14] A. Gahlawat. Analysis and control of parabolic partial differential equations with application to Tokamaks using sum-of-squares polynomials.
PhD thesis, Illinois Institute of Technology, Université de Grenoble,
2015.
[15] A. Gahlawat and M. M. Peet. Designing observer-based controllers
for PDE systems: A heat-conducting rod with point observation and
boundary control. In In Proc. of IEEE Conference on Decision and
Control and European Control Conference, pages 6985–6990, 2011.
[16] I. Gohberg and M. A. Kaashoek. Time varying linear systems with
boundary conditions and integral operators. I. The transfer operator and
its properties. Integral Equations and Operator Theory, 7:325–391,
1984.
[17] M. Krstic and A. Smyshlyaev. Adaptive boundary control for unstable
parabolic PDEs, Part I: Lyapunov design. IEEE Transactions on
Automatic Control, 53:1575–1591, 2008.
[18] M. Krstic and A. Smyshlyaev. Boundary control of PDEs: A course on
backstepping designs. Society for Industrial Mathematics, 2008.

[19] K. Kunisch and S. Volkwein. Galerkin proper orthogonal decomposition
methods for parabolic problems. Numerische mathematik, 90:117–148,
2001.
[20] I. Lasiecka and R. Triggiani. Control and stabilization of distributed
parameter systems; theoretical and computational aspects. Technical
report, DTIC Document, 1994.
[21] I. Lasiecka and R. Triggiani. Control theory for partial differential
equations: Volume 1, Abstract parabolic systems: Continuous and approximation theories. Cambridge University Press, 2000.
[22] J. L. Lions and E. Magenes. Non-Homogeneous Boundary Value
Problems and Applications. Springer-Verlag, Berlin/New York, 1972.
[23] K. Morris and C. Navasca. Approximation of low rank solutions for
linear quadratic control of partial differential equations. Computational
Optimization and Applications, 46:93–111, 2010.
[24] K. A. Morris. Design of finite-dimensional controllers for infinitedimensional systems by approximation. Journal of Mathematical Systems, Estimation, and Control, 4:30, 1994.
[25] J. D. Murray. Mathematical biology. Springer, 2002.
[26] Y. V. Orlov and L. T. Aguilar. Advanced H∞ control: Towards
nonsmooth theory and applications. Springer Science & Business
Media, 2014.
[27] A. Papachristodoulou and M. M. Peet. On the analysis of systems
described by classes of partial differential equations. In Proc. of IEEE
Conference on Decision and Control, pages 747–752, 2006.
[28] M. M. Peet. LMI parametrization of Lyapunov functions for infinitedimensional systems: A framework. In Proc. of American Control
Conference, pages 359–366, 2014.
[29] S. Prajna, A. Papachristodoulou, and P. A. Parrilo. Introducing SOSTOOLS: A general purpose sum of squares programming solver. In
Proc. of IEEE Conference on Decision and Control, pages 741–746,
2002.
[30] M. Putinar. Positive polynomials on compact semi-algebraic sets.
Indiana University Mathematics Journal, 42:969–984, 1993.
[31] A. Seuret and F. Gouaisbaut. On the use of the Wirtinger inequalities
for time-delay systems. In Proc. of 10th IFAC Workshop on Time Delay
Systems, pages 260–265, 2012.
[32] A. Smyshlyaev. Lyapunov adaptive boundary control for parabolic
PDEs with spatially varying coefficients. In Proc. of American Control
Conference, pages 41–48, 2006.
[33] A. Smyshlyaev and M. Krstic. Adaptive boundary control for unstable parabolic PDEs, Part II: Estimation-based designs. Automatica,
43:1543–1556, 2007.
[34] A. Smyshlyaev and M. Krstic. Adaptive boundary control for unstable
parabolic PDEs, Part III: Output feedback examples with swapping
identifiers. Automatica, 43:1557–1564, 2007.
[35] O. Staffans. Quadratic optimal control of stable well-posed linear
systems. Transactions of the American Mathematical Society, 349:3679–
3715, 1997.
[36] O. Staffans. Quadratic optimal control of well-posed linear systems.
SIAM Journal on Control and Optimization, 37:131–164, 1998.
[37] J. F. Sturm. Using SeDuMi 1.02, a MATLAB toolbox for optimization
over symmetric cones. Optimization Methods and Software, 11:625–653,
1999.
[38] V. L. Syrmos, C. T. Abdallah, P. R. Dorato, and K. Grigoriadis. Static
output feedback - A survey. Automatica, 33:125–137, 1997.
[39] K. Tanaka, H. Yoshida, H. Ohtake, and H. O. Wang. A sum-of-squares
approach to modeling and control of nonlinear dynamical systems
with polynomial fuzzy systems. IEEE Transactions on Fuzzy Systems,
17:911–922, 2009.
[40] A. Tits and Y. Yang. Globally convergent algorithms for robust pole
assignment by state feedback. IEEE Transactions on Automatic Control,
41:1432–1452, 1996.
[41] B. Van Keulen. H∞ -control for distributed parameter systems: a state
space approach. Birkhauser, 1993.
[42] M. Weiss and G. Weiss. Optimal control of stable weakly regular linear
systems. Mathematics of Control, Signals and Systems, 10:287–330,
1997.
[43] E. Witrant, E. Joffrin, S. Brémond, G. Giruzzi, D. Mazon, O. Barana, and
P. Moreau. A control-oriented model of the current profile in tokamak
plasma. Plasma Physics and Controlled Fusion, 49:1075, 2007.

Aditya Gahlawat received the B.Tech degree in mechanical engineering from
Punjabi University, Patiala, India in 2007, the M.S. degree in mechanical and

22

aerospace engineering from Illinois Institute of Technology, Chicago, USA
in 2009 and the Ph.D. degree in automatique-productique from Université
Grenoble Alpes, St. Martin d’Heres, France, in 2015. He is currently a
Ph.D. candidate in mechanical and aerospace engineering at Illinois Institute
of Technology, Chicago, USA His research focuses on the application of
convex optimization based methods for the analysis and control of systems
governed by partial differential equations with application to thermonuclear
fusion. Aditya Gahlawat was awarded the Chateaubriand fellowship in 2011
and 2012.

Matthew M. Peet received the B.Sc. degree in physics and in aerospace
engineering from the University of Texas, Austin, TX, USA, in 1999 and
the M.S. and Ph.D. degrees in aeronautics and astronautics from Stanford
University in 2001 and 2006, respectively. He was a Postdoctoral Fellow
at the National Institute for Research in Computer Science and Control
(INRIA), Paris, France, from 2006 to 2008. He was an Assistant Professor
of Aerospace Engineering in the Mechanical, Materials, and Aerospace
Engineering Department at the Illinois Institute of Technology in Chicago, IL,
USA, from 2008 to 2012. Currently, he is an Assistant Professor of Aerospace
Engineering in the School for the Engineering of Matter, Transport, and
Energy (SEMTE) at Arizona State University, Tempe, AZ, USA, and director
of the Cybernetic Systems and Controls Laboratory (CSCL). His research
interests are in the role of computation as it is applied to the understanding
and control of complex and large-scale systems with an emphasis on methods
such as SOS for the optimization of polynomials. Dr. Peet received a National
Science Foundation CAREER award in 2011.

