[0,0] | T.R

     T.R =
         [0,0,0] | [0,0,0] | [0,0,0]

Here, by assigning a value to dim, the parameters are adjusted to zero
matrices of appropriate dimensions. We note that, this command is not strictly
necessary, as the dimensions of T will also be automatically adjusted once we
specify the values of the parameters.

Next, we assign the interval $[2,3]$ on which the functions $\mathbf{x}_{1}\in
L_{2}^{3}[0,3]$ are defined, by setting the field I as

>> T.I = [2,3];

Since the parameters defining $\mathcal{T}$ also depend on $r,\nu\in[2,3]$, we
have to assign these variables as well. For this, we represent them by pvar
objects r and nu, and set the values of T.var1 and T.var2 as

     >> pvar r nu
     >> T.var1 = r;     T.var2 = nu;

Note that, if the parameters Qi and R.Ri are constant, there is no need to
declare the variables var1 or var2, in which case these fields will default to
var1=s and var2=theta. Having declared the variables, we finally set the
values of the parameters, by assigning them to the appropriate fields of T

     >> T.P = [1,0; 0,2; 3,4];
     >> T.Q1 = [r^2, 0, 0; 3, r^3, 0; 0, r+2*r^2, 0];
     >> T.Q2 = [-5*r, 6];
     >> T.R.R1 = [r, 2*nu, 3*r-nu]
     T =
              [1,0] | [r^2,0,0]
              [0,2] | [3,r^3,0]
              [3,4] | [0,2*r^2+r,0]
           ----------------------
           [-5*r,6] | T.R

     T.R =
         [0,0,0] | [r,2*nu,-nu+3*r] | [0,0,0]

dim is dependent on size of the 6 parameters P, Qi and R.Ri. Modifying those
parameters automatically changes the value stored in dim property. If the
dimensions of the parameters are incompatible, dim will store Nan as its value
to alert the user about the discrepancy.

#####  10.1.2 dpvar objects and the dopvar class

In addition to polynomial functions, polynomial decision variables can also be
declared in PIETOOLS. Such decision variables are used in polynomial
optimization programs, optimizing over the values of the coefficients defining
these polynomials. To declare such polynomial decision variables, we use the
dpvar class, extending the polynomial class to represent polynomials with
decision variables. For example, to represent a variable quadratic polynomial

$\displaystyle p(c_{0},c_{1},c_{2};x)=c_{0}+c_{1}x+c_{2}x^{2},$

with unknown coefficient $\\{c_{0},c_{1},c_{2}\\}$, we use

     >> pvar x
     >> dpvar c0 c1 c2
     >> p = c0 + c1*x + c2*x^2
     p =
       c0 + c1*x + c2*x^2

Here, the second line decision variables $c_{0}$, $c_{1}$ and $c_{2}$, and the
third line uses these decision variables to declare the decision variable
polynomial $p(c_{0},c_{1},c_{2};x)$ as a dpvar object p. Crucially, this
variable is affine in the coefficients $c_{0}$, $c_{1}$ and $c_{2}$, as dpvar
objects can only be used to represent decision variable polynomials that are
affine with respect to their coefficients. This is because PIETOOLS cannot
tackle polynomial optimization programs that are nonlinear in the decision
variables, and accordingly, the dpvar structure has been built to exploit the
linearity of the decision variables to minimize computational effort.

Using polynomial decision variables, we can also define PI operator decision
variables, defining the parameters $Q_{1}$ through $R_{2}$ by polynomial
decision variables rather than polynomials. For example, suppose we have an
operator $\mathcal{D}:L_{2}^{2}[0,1]\rightarrow L_{2}^{2}[0,1]$ defined as

$\displaystyle\bigl{(}\mathcal{D}\mathbf{x}\bigr{)}(s)$
$\displaystyle=\underbrace{\begin{bmatrix}c_{1}&c_{2}s\\\
c_{2}s&c_{3}s^{2}\end{bmatrix}}_{R_{0}(s)}\mathbf{x}(s)+\int_{s}^{1}\underbrace{\begin{bmatrix}c_{4}s^{2}&c_{5}s\theta\\\
-c_{6}s\theta&c_{7}\theta^{2}\end{bmatrix}}_{R_{2}(s,\theta)}\mathbf{x}(\theta)d\theta,$
$\displaystyle s$ $\displaystyle\in[0,1]$

for $\mathbf{x}\in L_{2}^{2}[0,1]$, where $c_{1}$ through $c_{7}$ are unknown
coefficients. Then, we can declare the parameters $R_{1}$ and $R_{2}$ as dpvar
class objects by calling

     >> pvar s th
     >> dpvar c1 c2 c3 c4 c5 c6 c7
     >> R0 = [c1, c2*s; c2*s, c3*s^2];
     >> R2 = [c4, c5*s*th; c6*s*th, c7*th^2];

Then, we can declare $\mathcal{D}$ as a dopvar object D. Such dopvar objects
have a structure identical to that of opvar objects, with the only difference
being that the parameters P through R.R2 can be declared as dpvar objects.
Accordingly, we declare the dopvar object D in a similar manner as we would
and opvar objects:

     >> dopvar D;
     >> D.I = [0,1];
     >> D.var1 = s;     D.var2 = th;
     >> D.R.R0 = R0;    D.R.R1 = R2
     D =
           [] | []
           ---------
           [] | D.R

     D.R =
             [c1,c2*s] | [0,0] |      [c4,c5*s*th]
         [c2*s,c3*s^2] | [0,0] | [c6*s*th,c7*th^2]

####  10.2 Algebraic Operations on opvar Objects

In this section, we go over various methods that help in manipulating and
handling of opvar objects in PIETOOLS. In particular, in Subsection 10.2.1, we
show how the sum of two PI operators can be computed in PIETOOLS, followed by
the composition of PI operators in Subsection 10.2.2. In Subsection 10.2.3, we
then show how to take the adjoint of a PI operator, and finally, in Subsection
10.2.4, we show how a numerical inverse of a PI operator can be computed. To
illustrate each of these operations, we use the PI operators
$\mathcal{A},\mathcal{B}:\left[\begin{smallmatrix}\mathbb{R}^{2}\\\
L_{2}[-1,1]\end{smallmatrix}\right]\rightarrow\left[\begin{smallmatrix}\mathbb{R}^{2}\\\
L_{2}[-1,1]\end{smallmatrix}\right]$ defined as

$\displaystyle\bigl{(}\mathcal{A}\mathbf{x}\bigr{)}(s)$
$\displaystyle=\Biggl{[}\begin{array}[]{ll}\left[\begin{smallmatrix}1&0\\\
2&-1\end{smallmatrix}\right]x_{0}&\\!+\
\int_{-1}^{1}\left[\begin{smallmatrix}1-s\\\
s+1\end{smallmatrix}\right]\mathbf{x}_{1}(s)ds\\\
\left[\begin{smallmatrix}10s&-1\end{smallmatrix}\right]x_{0}&\\!+\
2\mathbf{x}_{1}(s)+\int_{-1}^{s}(s-\theta)\mathbf{x}_{1}(\theta)d\theta+\int_{s}^{1}(s-\theta)\mathbf{x}_{1}(\theta)d\theta\end{array}\Biggr{]},$
(10.6) $\displaystyle\bigl{(}\mathcal{B}\mathbf{x}\bigr{)}(s)$
$\displaystyle=\Biggl{[}\begin{array}[]{ll}\left[\begin{smallmatrix}1&0\\\
0&3\end{smallmatrix}\right]x_{0}&\\\
\left[\begin{smallmatrix}5s&-s\end{smallmatrix}\right]x_{0}&\\!+\
s^{2}\mathbf{x}_{1}(s)+\int_{s}^{1}\theta\mathbf{x}_{1}(\theta)d\theta\end{array}\Biggr{]},$
$\displaystyle s$ $\displaystyle\in[-1,1],$ (10.8)

for $\mathbf{x}=\left[\begin{smallmatrix}x_{0}\\\
\mathbf{x}_{1}\end{smallmatrix}\right]\in\left[\begin{smallmatrix}\mathbb{R}^{2}\\\
L_{2}[-1,1]\end{smallmatrix}\right]$. We declare these operators as

     >> pvar s th
     >> opvar A B;
     >> A.I = [-1,1];         B.I = [-1,1];
     >> A.var1 = s;           B.var1 = s;
     >> A.var2 = th;          B.var2 = th;
     >> A.P = [1,0;2,-1];     B.P = [1,0;0,3];
     >> A.Q1 = [1-s;s+1];
     >> A.Q2 = [10*s,-1];     B.Q2 = [5*s,-s];
     >> A.R.R0 = 2;           B.R.R0 = s^2;
     >> A.R.R1 = (s-th);
     >> A.R.R2 = (s-th);      B.R.R2 = th;

#####  10.2.1 Addition (A+B)

opvar objects, A and B, can be added simply by using the command

>> A+B

For two opvar objects to be added, they must have same dimensions
(A.dim=B.dim), domains (A.I=B.I), and variables (A.var1=B.var1). Furthermore,
if A (or B) is a scalar then PIETOOLS considers that as adding A*I (or B*I)
where I is an identity matrix. Again, this operation is appropriate if and
only if dimensions match. Similarly, if A (or B) is a matrix with matching
dimension, it can be added to opvar B (or A) using the same command.

###### Example

Adding the opvar objects A and B corresponding to operators
$\mathcal{A},\mathcal{B}$ defined as in Equation (10.6), we find

     >> C = A+B
     C =
                [2,0] | [-s+1]
                [2,2] | [s+1]
            ------------------
          [15*s,-s-1] | C.R

     C.R =
         [s^2+2] | [s-th] | [s]

suggesting that, for $\mathbf{x}=\left[\begin{smallmatrix}x_{0}\\\
\mathbf{x}_{1}\end{smallmatrix}\right]\in\left[\begin{smallmatrix}\mathbb{R}^{2}\\\
L_{2}[-1,1]\end{smallmatrix}\right]$,

$\displaystyle\bigl{(}\mathcal{A}\mathbf{x}\bigr{)}(s)+\bigl{(}\mathcal{B}\mathbf{x}\bigr{)}(s)$
$\displaystyle=\Biggl{[}\begin{array}[]{ll}\left[\begin{smallmatrix}2&0\\\
2&2\end{smallmatrix}\right]x_{0}&\\!+\
\int_{-1}^{1}\left[\begin{smallmatrix}1-s\\\
s+1\end{smallmatrix}\right]\mathbf{x}_{1}(s)ds\\\
\left[\begin{smallmatrix}15s&-s-1\end{smallmatrix}\right]x_{0}&\\!+\
(s^{2}+2)\mathbf{x}_{1}(s)+\int_{-1}^{s}(s-\theta)\mathbf{x}_{1}(\theta)d\theta+\int_{s}^{1}s\mathbf{x}_{1}(\theta)d\theta\end{array}\Biggr{]}.$

#####  10.2.2 Multiplication (A*B)

opvar objects, A and B, can be composed simply by using the command

>> A*B

For two opvar objects to be composed, they must have the same domains
(A.I=A.B), the same variables (A.var1=B.var1 and A.var2=B.var2), and the
output dimension of B must match the input dimension of A
(A.dim(:,2)=B.dim(:,1)). Furthermore, if A (or B) is a scalar then PIETOOLS
considers that as a scalar multiplication operation, thus multiplying all
parameters of B (or A) by that value.

###### Example

Composing the opvar objects A and B corresponding to operators
$\mathcal{A},\mathcal{B}$ defined as in Equation (10.6), we find

     >> C = A*B
     C =
                  [-2.3333,0.66667] | [-1.5*s^3+2*s^2+1.5*s]
                   [5.3333,-3.6667] | [1.5*s^3+2*s^2+0.5*s]
          -------------------------------------------
          [20*s-3.3333,-2*s-2.3333] | C.R

     C.R =
         [2*s^2] | [2*s*th^2-1.5*th^3+s*th+0.5*th] | [2*s*th^2-1.5*th^3+s*th+2.5*th]

suggesting that, for $\mathbf{x}=\left[\begin{smallmatrix}x_{0}\\\
\mathbf{x}_{1}\end{smallmatrix}\right]\in\left[\begin{smallmatrix}\mathbb{R}^{2}\\\
L_{2}[-1,1]\end{smallmatrix}\right]$,

$\displaystyle\Bigl{(}\mathcal{A}\bigl{(}\mathcal{B}\mathbf{x}\bigr{)}\Bigr{)}(s)$
$\displaystyle=\left[\begin{array}[]{ll}\left[\begin{smallmatrix}-2\frac{1}{3}&\frac{2}{3}\\\
5\frac{1}{3}&-3\frac{2}{3}\end{smallmatrix}\right]x_{0}&\\!+\
\int_{-1}^{1}\left[\begin{smallmatrix}-1\frac{1}{2}s^{3}+2s^{2}+1\frac{1}{2}s\\\
1\frac{1}{2}s^{3}+2s^{2}+\frac{1}{2}s\end{smallmatrix}\right]\mathbf{x}_{1}(s)ds\\\
\left[\begin{smallmatrix}20s-3\frac{1}{3}&-2s-2\frac{1}{3}\end{smallmatrix}\right]x_{0}&\\!+\
2s^{2}\mathbf{x}_{1}(s)+\int_{-1}^{s}(2s\theta^{2}-1\frac{1}{2}\theta^{3}+s\theta+\frac{1}{2}\theta)\mathbf{x}_{1}(\theta)d\theta\\\
&\quad+\int_{s}^{1}(2s\theta^{2}-1\frac{1}{2}\theta^{3}+s\theta+2\frac{1}{2}\theta)\mathbf{x}_{1}(\theta)d\theta\end{array}\right].$

Although dopvar objects can be multiplied with opvar objects and vice versa,
producing a dopvar object in both cases, it is not possible to compute the
composition of two dopvar objects. This is because dopvar objects depend
linearly (affinely) on the decision variables, and the composition of two
dopvar objects would require taking the product of decision variables.
Similarly for dopvar2d objects.

#####  10.2.3 Adjoint (A’)

The adjoint of an opvar object A can be calculated using the command

>> A’

For an operator $\mathcal{A}:RL^{n_{0},n_{1}}[a,b]\rightarrow
RL^{m_{0},m_{1}}[a,b]$, the adjoint
$\mathcal{A}^{*}:RL^{m_{0},m_{1}}[a,b]\rightarrow RL^{n_{0},n_{1}}[a,b]$ will
be such that, for any $\mathbf{x}\in RL^{n_{0},n_{1}}[a,b]$ and $\mathbf{y}\in
RL^{m_{0},m_{1}}[a,b]$,

$\displaystyle\left\langle{\mathcal{A}\mathbf{x}},{\mathbf{y}}\right\rangle_{RL}=\left\langle{\mathbf{x}},{\mathcal{A}^{*}\mathbf{y}}\right\rangle_{RL},$

where for $\mathbf{x}=\left[\begin{smallmatrix}x_{0}\\\
\mathbf{x}_{1}\end{smallmatrix}\right]\in\left[\begin{smallmatrix}\mathbb{R}^{n_{0}}\\\
L_{2}^{n_{1}}[a,b]\end{smallmatrix}\right]=RL^{n_{0},n_{1}}[a,b]$ and
$\mathbf{y}=\left[\begin{smallmatrix}y_{0}\\\
\mathbf{y}_{1}\end{smallmatrix}\right]\in\left[\begin{smallmatrix}\mathbb{R}^{n_{0}}\\\
L_{2}^{n_{1}}[a,b]\end{smallmatrix}\right]=RL^{n_{0},n_{1}}[a,b]$,

$\displaystyle\left\langle{\mathbf{x}},{\mathbf{y}}\right\rangle_{RL}:=\left\langle{x_{0}},{y_{0}}\right\rangle+\left\langle{\mathbf{x}_{1}},{\mathbf{y}_{1}}\right\rangle_{L_{2}}=x_{0}^{T}y_{0}+\int_{a}^{b}[\mathbf{x}_{1}(s)]^{T}\mathbf{y}_{1}(s)ds$

###### Example

Computing the adjoint of the opvar object A corresponding to operator
$\mathcal{A}$ defined as in Equation (10.6), we find

     >> AT = A’
     AT =
                [1,2] | [10*s]
               [0,-1] | [-1]
           ------------------
           [-s+1,s+1] | AT.R

     AT.R =
         [2] | [-s+th] | [-s+th]

suggesting that, for $\mathbf{x}=\left[\begin{smallmatrix}x_{0}\\\
\mathbf{x}_{1}\end{smallmatrix}\right]\in\left[\begin{smallmatrix}\mathbb{R}^{2}\\\
L_{2}[-1,1]\end{smallmatrix}\right]$,

$\displaystyle\bigl{(}\mathcal{A}^{*}\mathbf{x}\bigr{)}(s)$
$\displaystyle=\left[\begin{array}[]{ll}\left[\begin{smallmatrix}1&2\\\
0&-1\end{smallmatrix}\right]x_{0}&\\!+\
\int_{-1}^{1}\left[\begin{smallmatrix}10s\\\
-1\end{smallmatrix}\right]\mathbf{x}_{1}(s)ds\\\
\left[\begin{smallmatrix}1-s&s+1\end{smallmatrix}\right]x_{0}&\\!+\
2\mathbf{x}_{1}(s)+\int_{-1}^{1}(\theta-s)\mathbf{x}_{1}(\theta)d\theta+\int_{s}^{1}(\theta-s)\mathbf{x}_{1}(\theta)d\theta\end{array}\right].$

#####  10.2.4 Inverse (inv_opvar(A))

The inverse of an opvar object A can be numerically calculated, using the
function

>> inv_opvar(A)

See Lemma 9 for details on Inversion formulae.

###### Example

Computing the inverse of the opvar object A corresponding to operator
$\mathcal{A}$ defined as in Equation (10.6), we find

     >> Ainv = inv_opvar(A)
     Ainv =
                        [-0.2,-0.4] | [  0.3*s + 0.2]
                          [1.2,0.4] | [ -0.3*s - 0.7]
           ------------------------------------------
           [ 0.3*s+0.7,1.35*s+0.65] | AT.R

     Ainv.R =
         [0.5] | [-1.2*s*th-0.675*s-0.3*th-0.575] | [-1.2*s*th-0.675*s-0.3*th-0.575]

suggesting that, for $\mathbf{x}=\left[\begin{smallmatrix}x_{0}\\\
\mathbf{x}_{1}\end{smallmatrix}\right]\in\left[\begin{smallmatrix}\mathbb{R}^{2}\\\
L_{2}[-1,1]\end{smallmatrix}\right]$

$\displaystyle\bigl{(}\mathcal{A}^{-1}\mathbf{x}\bigr{)}(s)$
$\displaystyle=\left[\begin{array}[]{ll}\left[\begin{smallmatrix}-0.2&-0.4\\\
1.2&0.4\end{smallmatrix}\right]x_{0}&\\!+\
\int_{-1}^{1}\left[\begin{smallmatrix}0.3s+0.2\\\
-0.3s-0.7\end{smallmatrix}\right]\mathbf{x}_{1}(s)ds\\\
\left[\begin{smallmatrix}0.3s+0.7&1.35s+0.65\end{smallmatrix}\right]x_{0}&\\!+\
0.5\mathbf{x}_{1}(s)+\int_{-1}^{1}(-1.2s\theta-0.675s-0.3\theta-0.575)\mathbf{x}_{1}(\theta)d\theta\\\
&\quad+\int_{s}^{1}(-1.2s\theta-0.675s-0.3\theta-0.575)\mathbf{x}_{1}(\theta)d\theta\end{array}\right].$

An inverse function has not been defined for opvar2d, dopvar, or dopvar2d
objects.

####  10.3 Matrix Operations on opvar Objects

In this section, we show how matrix operations on opvar objects can be
performed. In particular, in Subsection 10.3.1 we show how opvar] objects can
be concatenated, and in Subsection 10.3.2, we show how desired rows and
columns of opvar objects can be extracted. Although we explain these
operations only for opvar objects, they can also be applied to dopvar,
opvar2d, and dopvar2d objects. For the purpose of illustration, we once more
use the 4-PI operators
$\mathcal{A},\mathcal{B}:\left[\begin{smallmatrix}\mathbb{R}^{2}\\\
L_{2}[-1,1]\end{smallmatrix}\right]\rightarrow\left[\begin{smallmatrix}\mathbb{R}^{2}\\\
L_{2}[-1,1]\end{smallmatrix}\right]$ defined in Equation (10.6), represented
in PIETOOLS by the opvar objects A and B,

     >> A
     A =
              [1,0] | [-s+1]
             [2,-1] | [s+1]
          ----------------
          [10*s,-1] | A.R

     A.R =
        [2] | [s-th] | [s-th]
     >> B
     B =
             [1,0] | [0]
             [0,3] | [0]
          ---------------
          [5*s,-s] | B.R

     B.R =
        [s^2] | [0] | [th]

#####  10.3.1 Concatenation ([A,B])

Just like matrices, multiple opvar objects can be concatenated, provided the
dimensions match. In particular, two opvar objects A and B can be horizontally
or vertically concatenated by respectively using the command

>> [A B] % for horizontal concatenation
>> [A; B] % for vertical concatenation

Note that concatenation of opvar objects is allowed only if their spatial
domain is the same (A.I=B.I), and the variables involved in each are identical
(A.var1=B.var1 and A.var2=B.var2). Moreover, A and B can be horizontally
concatenated only if they have the same row dimensions (A.dim(:,1)=B.dim(:,2),
and they can be concatenated vertically only if they have the same column
dimensions (A.dim(:,2)=B.dim(:,2)).

###### Example

Horizontally concatenating the opvar objects A and B corresponding to
operators $\mathcal{A},\mathcal{B}$ defined as in Equation (10.6), we find

     >> C = [A, B]
     C =
                 [1,0,1,0] | [-s+1,0]
                [2,-1,0,3] | [s+1,0]
          -----------------------
          [10*s,-1,5*s,-s] | C.R

     C.R =
         [2,s^2] | [s-th,0] | [s-th,th]

Note that, since
$\mathcal{A},\mathcal{B}:\left[\begin{smallmatrix}\mathbb{R}^{2}\\\
L_{2}[-1,1]\end{smallmatrix}\right]\rightarrow\left[\begin{smallmatrix}\mathbb{R}^{2}\\\
L_{2}[-1,1]\end{smallmatrix}\right]$, we have
$\mathcal{C}:\left[\begin{smallmatrix}\mathbb{R}^{4}\\\
L_{2}^{2}[-1,1]\end{smallmatrix}\right]\rightarrow\left[\begin{smallmatrix}\mathbb{R}^{2}\\\
L_{2}[-1,1]\end{smallmatrix}\right]$, not
$\mathcal{C}:\left[\begin{smallmatrix}\mathbb{R}\\\ L_{2}[-1,1]\\\
\mathbb{R}\\\
L_{2}[-1,1]\end{smallmatrix}\right]\rightarrow\left[\begin{smallmatrix}\mathbb{R}^{2}\\\
L_{2}[-1,1]\end{smallmatrix}\right]$. Similarly, taking the vertical
concatenation,

     >> C = [A; B]
     C =
              [1,0] | [-s+1]
             [2,-1] | [s+1]
              [1,0] | [0]
              [0,3] | [0]
          ----------------
          [10*s,-1] | C.R
           [5*s,-s] |

     C.R =
           [2] | [s-th] | [s-th]
         [s^2] |    [0] |   [th]

the resulting operator $\mathcal{C}$ will map
$\left[\begin{smallmatrix}\mathbb{R}^{2}\\\
L_{2}[-1,1]\end{smallmatrix}\right]\rightarrow\left[\begin{smallmatrix}\mathbb{R}^{4}\\\
L_{2}^{2}[-1,1]\end{smallmatrix}\right]$, as PIETOOLS cannot represent
operators mapping e.g. $\left[\begin{smallmatrix}\mathbb{R}^{2}\\\
L_{2}[-1,1]\end{smallmatrix}\right]\rightarrow\left[\begin{smallmatrix}\mathbb{R}\\\
L_{2}[-1,1]\\\ \mathbb{R}\\\ L_{2}[-1,1]\end{smallmatrix}\right]$.

#####  10.3.2 Subs-indexing (A(i,j))

4-PI operators can also be sliced the way matrices are sliced in matrices. The
index slicing is performed in the same manner as matrices.

>> T(row_ind, col_ind)

Indexing 4-PI operators is slightly different from matrix indexing due to
presence of multiple components. These components can be visualized as being
stacked as in a matrix:

$\displaystyle\texttt{B}=\left[\begin{array}[]{c|c}\texttt{P}&\texttt{Q1}\\\
\hline\cr\texttt{Q2}&\texttt{Ri}\end{array}\right]$

Then, row indices specified in row_ind correspond to the rows in this big
matrix. Column indices, col_ind, are associated with the columns of this big
matrix in similar manner. The retrieved slices are put in appropriate
components and a 4-PI operator is returned. Note the bottom-right block of the
big matrix B has 3 components in Ri. If indices in the slice correspond to
rows and columns in this block, then the slice is extracted from all three
components and stored in a Ri part of the new sliced PI operator.

###### Example

Extracting row 3 and columns 1, and 3 of the opvar object A corresponding to
the operator $\mathcal{A},\mathcal{B}$ defined as in Equation (10.6), we find

     >> C = A(3,[1,3])
     C =
              [] | []
          -------------
          [10*s] | C.R
     C.R =
         [2] | [s-th] | [s-th]

####  10.4 Additional Methods for opvar Objects

There are some additional functions included in PIETOOLS that can be used in
debugging or as the user sees fit. In this section, we compile the list of
those functions, without going into details or explanation. However, users can
find additional information by using help command in MATLAB.

Function Name |  Description
---|---
A==B |  The function checks whether the variables, domain, dimensions and parameters of the operators A and B are equal, and returns a binary value 1 if this is the case, or 0 if it is not. The function can also be used to check if an operator A has all parameters equal to zero operator by calling A==0.
isvalid(P) |  The function returns a logical value. 0 is everything is in order, 1 if the object has incompatible dimensions, 2 if property P is not a matrix, 3 if properties Q1, Q2 or R0 are not polynomials in $s$, 4 if properties R1 or R2 are not polynomials in $s$ and $\theta$. For opvar2d objects, returns boolean value true or false to indicate if P is appropriate or not.
degbalance(T) |  Estimates polynomial degrees needed to create an opvar object Q associated to a positive PI operator $\mathcal{Q}\succeq 0$ in poslpivar, such that T=Q has at least one solution.
getdeg(T) |  Returns highest and lowest degree of $s$ and $\theta$ in the components of the opvar object T.
rand_opvar(dim, deg) |  Creates a random opvar object of specified dimensions dim and polynomial degrees deg.
show(T,opts) |  Alternative display format for opvar objects with optional argument to omit selected properties from display output. Not defined for opvar2d objects.
opvar_postest(T) |  Numerically test for sign definiteness of T. Returns -1 if negative definite, 0 if indefinite and 1 if positive definite. Use opvar_postest_2d for opvar2d objects.
diff_opvar(T) |  Returns composition of derivative operator with opvar T as described in Lem. 10. Use diff(T) for opvar2d objects.

##  Part III Examples and Applications

###  Chapter 11 PIETOOLS Demonstrations

In this Chapter, we illustrate several applications of PIE simulation and LPI
programming, and how each of these problems can be implemented in PIETOOLS.
Each of these problems has also been implemented as a DEMO file in PIETOOLS,
which can be found in the PIETOOLS_demos directory.

####  11.1 DEMO 1: Simple Stability, Simulation and Control Problem

See Chapter 2 for a description.

####  11.2 DEMO 2: Estimating the Volterra Operator Norm

The Volterra integral operator $\mathcal{T}:L_{2}[0,1]\rightarrow L_{2}[0,1]$
is perhaps the simplest example of a 3-PI operator, defined as

$\displaystyle\bigl{(}\mathcal{T}\mathbf{x}\bigr{)}(s)$
$\displaystyle=\int_{0}^{s}\mathbf{x}(\theta)d\theta,$ $\displaystyle s$
$\displaystyle\in[0,1],$

for any $\mathbf{x}\in L_{2}[0,1]$. In PIETOOLS, this operator can be easily
declared as

     a=0;    b=1;
     opvar Top;
     Top.R.R1 = 1;   Top.I = [a,b];

Then, an upper bound on the norm of this operator can be computed by solving
the LPI

$\displaystyle\min_{\gamma\geq 0}\ \gamma,$ s.t.
$\displaystyle\mathcal{T}^{*}\mathcal{T}\leq\gamma,$

so that case $C=\sqrt{\gamma}$ for any feasible value $\gamma$ is an upper
bound on the norm of $\mathcal{T}$. This optimization problem is an LPI, that
can be declared and solved in PIETOOLS as

     % First, define dpvar gam and set up an optimization problem
     vars = [Top.var1;Top.var2];     % Free vars in optimization problem
     dpvar gam;                      % Decision var in optimization problem
     prob = sosprogram(vars,gam);

     % Next, set gam as objective function min{gam}
     prob = sossetobj(prob, gam);

     % Then, enforce the constraint Top’*Top-gam<=0
     opts.psatz = 1;                             % Allow Top’*Top-gam>0 outside of [a,b]
     prob = lpi_ineq(prob,-(Top’*Top-gam),opts); % lpi_ineq(prob,Q) enforces Q>=0

     % Finally, solve and retrieve the solution
     prob = sossolve(prob);
     operator_norm = sqrt(double(sosgetsol(prob,gam)));

This code can also be run by calling “volterra_operator_norm_DEMO”. We obtain
an upper bound $C=0.68698$ on the induced norm $\|\mathcal{T}\|$ of the
Volterra operator. The exact value of the induced norm of this operator is
known to be equal to $\|\mathcal{T}\|=\frac{2}{\pi}=0.6366...$.

####  11.3 DEMO 3: Solving the Poincaré Inequality

In a one dimensional domain $\Omega=[a,b]$, the Poincaré inequality imposes a
bound on the norm of a function $x(s)$ in terms of the spatial derivative
$\partial_{s}x(s)$ of this function,

$\displaystyle\|x\|_{L_{2}}$ $\displaystyle\leq
C\|\partial_{s}x\|_{L_{2}},\hskip 56.9055pt\forall x\in W_{1}[a,b],$ where
$\displaystyle W_{1}[a,b]$ $\displaystyle:=\\{x\in
L_{2}[a,b]\mid\partial_{s}x\in L_{2}[a,b],x(a)=x(b)=0\\}.$

In PIETOOLS, we can find a value $C$ that satisfies this inequality, by
solving the optimization problem

$\displaystyle\min_{\gamma\geq 0}\ \gamma,$ s.t.
$\displaystyle\left\langle{x},{x}\right\rangle_{L_{2}}-\gamma\left\langle{\partial_{s}x},{\partial_{s}x}\right\rangle_{L_{2}}\leq
0,$ $\displaystyle\forall x$ $\displaystyle\in W_{1}[a,b]$

in which case $C=\sqrt{\gamma}$ satisfies the Poincaré inequality. To declare
this problem, we first note that we can represent $x$ and $\partial_{s}x$ in
terms of a fundamental state $\partial_{s}^{2}x\in L_{2}[a,b]$, which is free
of the boundary conditions and continuity constraints imposed upon $x$ and
$\partial_{s}x$. In particular, we first declare a PDE

$\displaystyle\dot{x}(t,s)$ $\displaystyle=\partial_{s}x(t,s),$ $\displaystyle
s$ $\displaystyle\in[a,b],$ $\displaystyle x(t,a)$ $\displaystyle=x(t,b)=0,$

as a pde_struct object by calling

     % % Initialize the PDE structure and spatial variable s in [a,b]
     pvar s theta;
     pde_struct PDE;
     a = 0;      b = 1;

     % % Declare the state variables x(t,s)
     PDE.x{1}.vars = s;
     PDE.x{1}.dom = [a,b];
     PDE.x{1}.diff = 2;          % Let x be second order differentiable wrt s.

     % % Declare the PDE \dot{x}(t,s) = \partial_{s} x(t,s)
     PDE.x{1}.term{1}.D = 1;     % Order of the derivative wrt s

     % % Declare the boundary conditions x(t,a) = x(t,b) = 0
     PDE.BC{1}.term{1}.loc = a;      % Evaluate x at s=a
     PDE.BC{2}.term{1}.loc = b;      % Evaluate x at s=b

     % % Initialize the system
     PDE = initialize(PDE);

Here, we explicitly indicate that the function $x$ must be second order
differentiable, as we wish to derive a PIE representation in terms of the
fundamental state $x_{\text{f}}:=\partial_{s}^{2}x$. To obtain this PIE
representation, we call convert,

     PIE = convert(PDE,’pie’);
     H2 = PIE.T;     % H2 x_{ss} = x
     H1 = PIE.A;     % H1 x_{ss} = x_{s}

arriving at an equivalent representation of the PDE as

$\displaystyle\mathcal{H}_{2}\dot{x}_{\text{f}}(t,s)$
$\displaystyle=\mathcal{H}_{1}x_{\text{f}}(t,s),$ $\displaystyle s$
$\displaystyle\in[a,b].$

In this representation, the fundamental state
$x_{\text{f}}:=\partial_{s}^{2}x$ is free of any boundary conditions and
continuity constraints. Moreover, we note that

$\displaystyle\mathcal{H}_{2}x_{\text{f}}(t,s)$
$\displaystyle=x(t,s),\qquad\text{and},\qquad\mathcal{H}_{1}x_{\text{f}}(t,s)$
$\displaystyle=\partial_{s}x(t,s).$

As such, the Poincaré inequality optimization problem can be equivalently
represented as

$\displaystyle\min_{\gamma\geq 0}\
\gamma,\qquad\text{s.t.}\quad\left\langle{\mathcal{H}_{2}x_{\text{f}}},{\mathcal{H}_{2}x_{\text{f}}}\right\rangle-\gamma\left\langle{\mathcal{H}_{1}x_{\text{f}}},{\mathcal{H}_{1}x_{\text{f}}}\right\rangle\leq
0,\qquad\forall x_{\text{f}}\in L_{2}[a,b]$

giving rise to an LPI

$\displaystyle\min_{\gamma\geq 0}\
\gamma,\qquad\text{s.t.}\quad\mathcal{H}_{2}^{*}\mathcal{H}_{2}-\gamma\mathcal{H}_{1}^{*}\mathcal{H}_{1}\leq
0.$

We declare and solve this LPI in PIETOOLS as

     % % First, define dpvar gam and set up an optimization problem
     dpvar gam;
     vars = [H2.var1; H2.var2];
     prob = sosprogram(vars,gam);

     % % Set gam as objective function to minimize
     prob = sossetobj(prob, gam);

     % % Set up the constraint H2’*H2-gam H1’*H1<=0
     opts.psatz = 1;     % allow H2’*H2 > gam H1’*H1 outside of [a,b]
     prob = lpi_ineq(prob,-(H2’*H2-gam*H1’*H1),opts);

     % Solve and retrieve the solution
     prob = sossolve(prob);
     poincare_constant = sqrt(double(sosgetsol(prob,gam)));

This code can also be run calling the function “poincare_inequality_DEMO”,
arriving at a constant $C=0.42664$ that satisfies the Poincaré inequality on
the domain $[a,b]=[0,1]$. On this domain, a minimal value of $C$ is known to
be $C=\frac{1}{\pi}=0.31831...$.

####  11.4 DEMO 4: Finding an Optimal Stability Parameter

We consider a reaction-diffusion equation on an interval $[a,b]$,

$\displaystyle\dot{\mathbf{x}}(t,s)$
$\displaystyle=\lambda\mathbf{x}(t,s)+\partial_{s}^{2}\mathbf{x}(t,s),$
$\displaystyle s$ $\displaystyle\in[a,b],\mathbf{x}(t,a)$
$\displaystyle=\mathbf{x}(t,b)=0,$

for some $\lambda\in\mathbb{R}$. On the interval $[a,b]=[0,1]$, this system is
know to be stable whenever $\lambda\leq\pi^{2}=9.8696...$. In PIETOOLS, we can
numerically estimate this limit using the function stability_PIETOOLS. In
particular, we may equivalent represent the PDE as a PIE of the form

$\displaystyle\bigl{(}\mathcal{T}\dot{\mathbf{x}}_{\text{f}}\bigr{)}(t,s$
$\displaystyle=\bigl{(}\mathcal{A}(\lambda)\mathbf{x}_{\text{f}}\bigr{)}(t,s),$
$\displaystyle s$ $\displaystyle\in[a,b],$

where $\mathbf{x}_{\text{f}}(t,s):=\partial_{s}^{2}\mathbf{x}(t,s)$. Then, a
maximal value of $\lambda$ for which the PIE is stable can be estimated by
solving the optimization problem

$\displaystyle\max_{\lambda\in\mathbb{R},\mathcal{P}\in\Pi}\ \lambda,$ s.t.
$\displaystyle\mathcal{P}\succ 0,$
$\displaystyle\mathcal{T}^{*}\mathcal{P}\mathcal{A}(\lambda)+\mathcal{A}^{*}(\lambda)\mathcal{P}\mathcal{T}\preccurlyeq
0.$

Unfortunately, both $\lambda$ and $\mathcal{P}$ are decision variables in this
optimization program, and so the product $\mathcal{P}\mathcal{A}(\lambda)$ is
not linear in the decision variables. As such, this problem cannot be directly
implemented as a convex optimization program. However, for any fixed value of
$\lambda$, stability of the PIE can be verified by testing feasibility of the
LPI

$\displaystyle\mathcal{P}\succ 0,$
$\displaystyle\mathcal{T}^{*}\mathcal{P}\mathcal{A}(\lambda)+\mathcal{A}^{*}(\lambda)\mathcal{P}\mathcal{T}\preccurlyeq
0,$

an optimization program that has already been implemented in
stability_PIETOOLS. Therefore, to estimate an upper bound on the value of
$\lambda$ for which the PDE is stable, we can test stability for given values
of $\lambda$, and perform bisection over some domain
$\lambda\in[\lambda_{\min},\lambda_{\max}]$ to find an optimal value. For our
demonstration, we use $\lambda_{\min}=0$ and $\lambda_{\max}=20$, testing
stability for $8$ values of $\lambda$ between these upper and lower bounds:

     %%% Set bisection limits for lam.
     lam_min = 0;        lam_max = 20;
     lam = 0.5*(lam_min + lam_max);
     n_iters = 8;

We initialize a PDE structure in Terms-format, so that we may easily update
the value of $\lambda$:

     %%% Initialize a PDE structure.
     a = 0;  b = 1;
     pvar s
     pde_struct PDE;
     PDE.x{1}.vars = s;
     PDE.x{1}.dom = [a,b];

     % Set the PDE \dot{x}(t,s) = lam*x(t,s) + x_{ss}(t,s);
     PDE.x{1}.term{1}.C = lam;
     PDE.x{1}.term{2}.D = 2;

     % Set the BCs x(t,a) = x(t,b) = 0;
     PDE.BC{1}.term{1}.loc = a;
     PDE.BC{2}.term{1}.loc = b;

Finally, we fix settings for the LPI program, using veryheavy settings to
achieve relatively high accuracy

     %%% Initialize settings for solving the LPI
     settings = lpisettings(’veryheavy’);
     if strcmp(settings.sos_opts.solver,’sedumi’)
         settings.sos_opts.params.fid = 0;   % Suppress output in command window
     end

Having initialized everything, we perform bisection on the value of $\lambda$.
In particular, for a given value of
$\lambda\in[\lambda_{\min},\lambda_{\max}]$, we update the PDE structure PDE,
compute the associated PIE structure PIE, and test stability of this PIE using
stability_PIETOOLS. Then, we check the value of feasratio in the output
optimization program structure, which should be close to 1 if the LPI was
successfully solved, and thus the system was found to be stable.

  * •

If the system is stable, then it is stable for any value of $\lambda$ smaller
than the value lam used in the test. As such, we update the value of
$\lambda_{\min}\leftarrow\lambda$, and repeat the test with a greater value
$\lambda=\frac{1}{2}(\lambda_{\min}+\lambda_{\max})$.

  * •

If stability could not be verified, then stability can also not be verified
for any value of $\lambda$ greater than the value lam used in the test. As
such, we update the value of $\lambda_{\max}\leftarrow\lambda$, and repeat the
test with a greater value
$\lambda=\frac{1}{2}(\lambda_{\min}+\lambda_{\max})$.

This algorithm can be implemeted as

     %%% Perform bisection on the value of lam
     for iter = 1:n_iters
         % Update the value of lam in the PDE.
         PDE.x{1}.term{1}.C = lam;

         % Update the PIE.
         PIE = convert(PDE,’pie’);

         % Re-run the stability test.
         prog = PIETOOLS_stability(PIE,settings);

         % Check if the system is stable
         if prog.solinfo.info.dinf || prog.solinfo.info.pinf ...
                || abs(prog.solinfo.info.feasratio - 1)>0.3
             % Stability cannot be verified, decreasing the value of lam...
             lam_max = lam;
             lam = 0.5*(lam_min + lam_max);
         else
             % The system is stable, trying a larger value of lam...
             lam_success = lam;
             lam_min = lam;
             lam = 0.5*(lam_min + lam_max);
         end
     end

This code can also be called using “stability_parameter_bisection_DEMO”.
Running this demo, we find that stability can be verified whenever
$\lambda\leq 9.8438$.

####  11.5 DEMO 5: Constructing and Simulating an Optimal Estimator

We consider a reaction-diffusion PDE, with an observed output $y$,

$\displaystyle\dot{\mathbf{x}}(t,s)$
$\displaystyle=\partial_{s}^{2}\mathbf{x}(t,s)+4\mathbf{x}(t,s)+w(t),$
$\displaystyle s$ $\displaystyle\in[0,1],$ with BCs $\displaystyle 0$
$\displaystyle=\mathbf{x}(t,0)=\partial_{s}\mathbf{x}(t,1),$ and outputs
$\displaystyle z(t)$ $\displaystyle=\int_{0}^{1}\mathbf{x}(t,s)ds+w(t),$
$\displaystyle y(t)$ $\displaystyle=\mathbf{x}(t,1).$ (11.1)

This PDE can be easily declared in PIETOOLS as

     pvar s t
     PDE = sys();
     x = state(’pde’);    w = state(’in’);
     y = state(’out’);    z = state(’out’);
     eqs = [diff(x,t) == diff(x,s,2) + 4*x + w;
            z == int(x,s,[0,1]) + w;
            y == subs(x,s,1);
            subs(x,s,0) == 0;
            subs(diff(x,s),s,1) == 0];
     PDE = addequation(PDE,eqs);
     PDE = setObserve(PDE,y);

at which point an equivalent PIE representation can be derived by calling
convert:

     PIE = convert(PDE,’pie’);        PIE = PIE.params;
     T = PIE.T;
     A = PIE.A;      C1 = PIE.C1;      C2 = PIE.C2;
     B1 = PIE.B1;    D11 = PIE.D11;    D21 = PIE.D21;

Then, the PDE (11.5) can be equivalently represented by a PIE

$\displaystyle\bigl{(}\mathcal{T}\dot{\mathbf{v}}\bigr{)}(t,s)$
$\displaystyle=\bigl{(}\mathcal{A}\mathbf{v}\bigr{)}(t,s)+\bigl{(}\mathcal{B}_{1}w\bigr{)}(t,s),$
$\displaystyle s$ $\displaystyle\in[0,1]$ $\displaystyle z(t)$
$\displaystyle=\bigl{(}\mathcal{C}_{1}\mathbf{v}\bigr{)}(t)+\bigl{(}\mathcal{D}_{11}w)(t),$
$\displaystyle y(t)$
$\displaystyle=\bigl{(}\mathcal{C}_{2}\mathbf{v}\bigr{)}(t)+\bigl{(}\mathcal{D}_{12}w)(t),$
(11.2)

where we define $\mathbf{v}:=\partial_{s}^{2}\mathbf{x}$, and where
$\mathbf{x}=\mathcal{T}\mathbf{v}$.

We consider the problem of designing an optimal estimator for the PIE (11.5).
In particular, we construct an estimator of the form

$\displaystyle\mathcal{T}\dot{\hat{\mathbf{v}}}(t)$
$\displaystyle=\mathcal{A}{\mathbf{\hat{v}}}(t)+\mathcal{L}\bigl{(}y(t)-\hat{y}(t)\bigr{)},$
$\displaystyle\hat{z}(t)$ $\displaystyle=\mathcal{C}_{1}\mathbf{\hat{v}}(t),$
$\displaystyle\hat{y}(t)$ $\displaystyle=\mathcal{C}_{2}\mathbf{\hat{v}}(t),$
(11.3)

so that the error $\mathbf{e}(t,s):=\hat{\mathbf{v}}(t,s)-\mathbf{v}(t,s)$ in
the state and $\tilde{z}(t):=\hat{z}(t)-z(t)$ in the output satisfy

$\displaystyle\mathcal{T}\dot{\mathbf{e}}(t)$
$\displaystyle=(\mathcal{A}-\mathcal{L}\mathcal{C}_{2})\mathbf{e}(t)-(\mathcal{B}_{1}+\mathcal{L}\mathcal{D}_{21})w(t)$
$\displaystyle\tilde{z}(t)$
$\displaystyle=\mathcal{C}_{1}\mathbf{e}(t)-\mathcal{D}_{11}w(t).$ (11.4)

The goal of $H_{\infty}$-optimal estimation, then, is to determine a value for
the observer operator $\mathcal{L}$ that minimizes the gain
$\gamma:=\frac{\|\tilde{z}\|_{L_{2}}}{\|w\|_{L_{2}}}$ from disturbances $w$ to
errors $\tilde{z}$ in the output. See also Section 13.2. To construct such an
operator, we solve the LPI,

$\displaystyle\min\limits_{\gamma,\mathcal{P},\mathcal{Z}}\leavevmode\nobreak\
\leavevmode\nobreak\ \gamma$ $\displaystyle\mathcal{P}\succ 0,$ $\displaystyle
Q:=\begin{bmatrix}-\gamma
I&-\mathcal{D}_{11}^{\top}&-(\mathcal{P}\mathcal{B}_{1}+\mathcal{Z}\mathcal{D}_{21})^{*}\mathcal{T}\\\
(\cdot)^{*}&-\gamma I&\mathcal{C}_{1}\\\
(\cdot)^{*}&(\cdot)^{*}&(\mathcal{P}\mathcal{A}+\mathcal{Z}\mathcal{C}_{2})^{*}\mathcal{T}+(\cdot)^{*}\end{bmatrix}\preccurlyeq
0$ (11.5)

so that, for any solution $(\gamma,\mathcal{P},\mathcal{Z})$ to this problem,
letting $\mathcal{L}:=\mathcal{P}^{-1}\mathcal{Z}$, the estimation error will
satisfy $\|\hat{z}-z\|\leq\gamma\|w\|$.

In Chapter 7, we showed in detail how the LPI (11.5) can be declared for a PIE
of the form (11.5). The code declaring and solving this LPI can also be found
in the demo file “Hinf_Optimal_Estimator_DEMO”. We will therefore not discuss
this implementation here. Instead, we compute an optimal value of the observer
operator using one of the pre-defined executive files, as

     settings = lpisettings(’heavy’);
     [prog, Lval, gam_val] = PIETOOLS_Hinf_estimator(PIE, settings);

returning a value Lval of the operator $\mathcal{L}$ such that the
$L_{2}$-gain $\frac{\|\tilde{z}\|_{L_{2}}}{\|w\|_{L_{2}}}$ is bounded from
above by gam_val.

Given the observer operator $\mathcal{L}$, we can construct the Estimator
(11.5), obtaining a PIE

$\displaystyle\left(\begin{bmatrix}\mathcal{T}&0\\\
0&\mathcal{T}\end{bmatrix}\begin{bmatrix}\dot{\mathbf{v}}\\\
\dot{\hat{\mathbf{v}}}\end{bmatrix}\right)(t,s)$
$\displaystyle=\left(\begin{bmatrix}\mathcal{A}&0\\\
-\mathcal{L}\mathcal{C}_{2}&\mathcal{A}+\mathcal{L}\mathcal{C}_{2}\end{bmatrix}\begin{bmatrix}\mathbf{v}\\\
\hat{\mathbf{v}}\end{bmatrix}\right)(t,s)+\left(\begin{bmatrix}\mathcal{B}_{1}\\\
\mathcal{L}\mathcal{D}_{21}\end{bmatrix}w\right)(t)$
$\displaystyle\begin{bmatrix}z\\\ \hat{z}\end{bmatrix}(t)$
$\displaystyle=\left(\begin{bmatrix}\mathcal{C}_{1}&0\\\
0&\mathcal{C}_{1}\end{bmatrix}\begin{bmatrix}\mathbf{v}\\\
\hat{\mathbf{v}}\end{bmatrix}\right)(t)+\left(\begin{bmatrix}\mathcal{D}_{11}\\\
0\end{bmatrix}w\right)(t)$

We can declare this PIE in PIETOOLS as

     % Construct the operators defining the PIE.
     T_CL = [T, 0*T; 0*T, T];
     A_CL = [A, 0*A; -Lval*C2, A+Lval*C2];   B_CL = [B1; Lval*D21];
     C_CL = [C1, 0*C1; 0*C1, C1];            D_CL = [D11; 0*D11];

     % Declare the PIE.
     PIE_CL = pie_struct();
     PIE_CL.vars = PIE.vars;
     PIE_CL.dom = PIE.dom;
     PIE_CL.T = T_CL;
     PIE_CL.A = A_CL;              PIE_CL.B1 = B_CL;
     PIE_CL.C1 = C_CL;             PIE_CL.D11 = D_CL;
     PIE_CL = initialize(PIE_CL);

Then, we can use PIESIM to simulate the evolution of the PDE state
$\mathbf{x}(t)=\bigl{(}\mathcal{T}\mathbf{v}\bigr{)}(t)$ and its estimate
$\hat{\mathbf{x}}(t)=\bigl{(}\mathcal{T}\hat{\mathbf{v}}\bigr{)}(t)$
associated to the system defined by PIE_CL. For this, we first declare initial
conditions $\left[\begin{smallmatrix}\mathbf{v}(0,s)\\\
\hat{\mathbf{v}}(0,s)\end{smallmatrix}\right]$ and values of the disturbance
$w(t)$ as

     % Declare initial conditions for the state components of the PIE
     syms st sx real
     uinput.ic.PDE = [-10*sx;    % Actual initial PIE state value
                      0];        % Estimated initial PIE state value

     % Declare the value of the disturbance w(t)
     uinput.w = 2*sin(pi*st);

Here we use symbolic objects st and sx to represent a temporal variable $t$
and spatial variable $s$ respectively. Note that, since we will be simulating
the PIE directly, the initial conditions uinput.ic.PDE will also correspond to
the initial values of the PIE state
$\left[\begin{smallmatrix}\mathbf{v}(0,s)\\\
\hat{\mathbf{v}}(0,s)\end{smallmatrix}\right]$, not to those of the PDE state.
Here, we let the initial PIE state $\mathbf{v}(0)$ be parabolic, and start
with an estimate of this state that is just $\hat{\mathbf{v}}(0)=0$. Given
these initial conditions, we then simulate the evolution of the PDE state
$\left[\begin{smallmatrix}\mathbf{x}\\\
\mathbf{\hat{x}}\end{smallmatrix}\right]=\left[\begin{smallmatrix}\mathcal{T}&0\\\
0&\mathcal{T}\end{smallmatrix}\right]\left[\begin{smallmatrix}\dot{\mathbf{v}}\\\
\dot{\hat{\mathbf{v}}}\end{smallmatrix}\right]$ using PIESIM as

     % Set options for the discretization and simulation:
     opts.plot = ’no’;   % Do not plot the final solution
     opts.N = 8;         % Expand using 8 Chebyshev polynomials
     opts.tf = 1;        % Simulate up to t = 1;
     opts.dt = 1e-3;     % Use time step of 10^-3
     opts.intScheme = 1; % Time-step using Backward Differentiation Formula (BDF)
     ndiff = [0,0,2];    % The PDE state involves 2 twice differentiable states

     % Simulate the solution to the PIE with estimator.
     [solution,grid] = PIESIM(PIE_CL,opts,uinput,ndiff);

     % Extract actual and estimated solution at each time step.
     x_act = reshape(solution.timedep.pde(:,1,:),opts.N+1,[]);
     x_est = reshape(solution.timedep.pde(:,2,:),opts.N+1,[]);
     tval = solution.timedep.dtime;

Here, we use opts.N=8 Chebyshev polynomials in the expansion of the solution,
and simulate up to opts.tf=1, taking time steps of opts.dt=1e-3. Having
obtained the solution, we then collect the actual values $\mathbf{x}(t,s)$ of
the PDE state at each time step and each grid point in x_act, and the
estimated values $\mathbf{\hat{x}}(t,s)$ of the PDE state at each time step
and each grid point in x_est. Plotting the obtained values at several grid
points, as well as the error in estimated state, we obtain a graph as in
Figure 11.1.

Figure 11.1: Simulated value of PDE state $\mathbf{x}(t,s)$ and estimated
state $\mathbf{\hat{x}}(t,s)$, along with the error
$\mathbf{e}=\hat{\mathbf{x}}(t,s)-mbf{x}(t,s)$ associated to the PDE (11.5) at
several grid points $s\in[0,1]$, using the Estimator (11.5) with operator
$\mathcal{L}$ computed by solving the $H_{\infty}$-optimal estimator LPI. See
the demo file “Hinf_Optimal_Estimator_DEMO”.

As Figure 11.1 shows, the value of the estimate state $\mathbf{\hat{x}}$ at
each of the grid points quickly converge to those of the actual state
$\mathbf{x}$, despite starting off with a rather poor estimate of
$\mathbf{\hat{x}}=0$. Within a time of 1, the values of the actual and
estimated state become indistinguishable (at the displayed scale), with the
error converging to zero. This is despite the fact that the value of the state
itself continues to increase, as the PDE (11.5) is unstable.

The full code constructing the optimal estimator, simulating the PDE state and
its estimate, and plotting the results has been included in PIETOOLS as the
demo file
“Hinf_Optimal_Estimator_DEMO”.

####  11.6 DEMO 6: $H_{\infty}$-optimal Controller synthesis for PDEs

We consider an unstable reaction-diffusion PDE, with output $z$, disturbance
$w$ and control input $u$

$\displaystyle\dot{\mathbf{x}}(t,s)$
$\displaystyle=\partial_{s}^{2}\mathbf{x}(t,s)+4\mathbf{x}(t,s)+sw(t)+su(t),$
$\displaystyle s$ $\displaystyle\in[0,1],$ with BCs $\displaystyle 0$
$\displaystyle=\mathbf{x}(t,0)=\partial_{s}\mathbf{x}(t,1),$ and outputs
$\displaystyle z(t)$
$\displaystyle=\begin{bmatrix}\int_{0}^{1}\mathbf{x}(t,s)ds+w(t)\\\
u(t)\end{bmatrix}.$ (11.6)

This PDE can be easily declared in PIETOOLS as

    pvar s t
    lam = 4;
    PDE = sys();
    x = state(’pde’);   w = state(’in’);
    z = state(’out’, 2);   u = state(’in’);
    eqs = [diff(x,t) == diff(x,s,2) + lam*x + s*w + s*u;
        z == [int(x,s,[0,1]) + w; u];
        subs(x,s,0)==0;
        subs(diff(x,s),s,1)==0];
    PDE = addequation(PDE,eqs);
    PDE = setControl(PDE,u);
    display_PDE(PDE);

Next, we convert the PDE system to a PIE by using the following code and
extract relevant PI operators for easier access:

    PIE = convert(PDE,’pie’);       PIE = PIE.params;
    T = PIE.T;
    A = PIE.A;      C1 = PIE.C1;    B2 = PIE.B2;
    B1 = PIE.B1;    D11 = PIE.D11;  D12 = PIE.D12;

Then, the PDE (11.6) can be equivalently represented by a PIE

$\displaystyle\bigl{(}\mathcal{T}\dot{\mathbf{v}}\bigr{)}(t,s)$
$\displaystyle=\bigl{(}\mathcal{A}\mathbf{v}\bigr{)}(t,s)+\bigl{(}\mathcal{B}_{1}w\bigr{)}(t,s)+\bigl{(}\mathcal{B}_{2}u\bigr{)}(t,s),$
$\displaystyle s$ $\displaystyle\in[0,1]$ $\displaystyle z(t)$
$\displaystyle=\bigl{(}\mathcal{C}_{1}\mathbf{v}\bigr{)}(t)+\bigl{(}\mathcal{D}_{11}w)(t)+\bigl{(}\mathcal{D}_{12}u)(t),$
(11.7)

where we define $\mathbf{v}:=\partial_{s}^{2}\mathbf{x}$ and
$\mathbf{x}=\mathcal{T}\mathbf{v}$.

We now attempt to find a state-feedback $H_{\infty}$-optimal controller for
the PIE (11.6). In particular, we use $u(t)=\mathcal{K}\mathbf{v}(t)$, where
$\mathcal{K}:L_{2}\to\mathbb{R}$ is a PI operator of the form

$\mathcal{K}\mathbf{v}=\int_{a}^{b}K(s)\mathbf{v}(s)ds.$

Then we get the closed loop system

$\displaystyle\mathcal{T}\dot{\mathbf{v}}(t)$
$\displaystyle=(\mathcal{A}+\mathcal{B}_{2}\mathcal{K}){\mathbf{v}}(t)+\mathcal{B}_{1}w(t),$
$\displaystyle z(t)$
$\displaystyle=(\mathcal{C}_{1}+\mathcal{D}_{12}\mathcal{K})\mathbf{z}(t)+\mathcal{D}_{11}w(t).$
(11.8)

Next, we solve the LPI for $H_{\infty}$-optimal control to find a
$\mathcal{K}$ that minimizes the gain
$\gamma:=\frac{\|z\|_{L_{2}}}{\|w\|_{L_{2}}}$ from disturbances $w$ to errors
$\tilde{z}$ in the output. See also Section 13.3. To find such an operator, we
solve the LPI,

$\displaystyle\min\limits_{\gamma,\mathcal{P},\mathcal{Z}}\leavevmode\nobreak\
\leavevmode\nobreak\ \gamma\leavevmode\nobreak\ s.t.\quad$
$\displaystyle\mathcal{P}\succ 0$ $\displaystyle\begin{bmatrix}-\gamma
I&\mathcal{D}_{11}&\mathcal{T}(\mathcal{P}\mathcal{C}_{1}+\mathcal{Z}\mathcal{D}_{12})\\\
(\cdot)^{*}&-\gamma I&\mathcal{B}_{1}^{*}\\\
(\cdot)^{*}&(\cdot)^{*}&\mathcal{T}(\mathcal{A}\mathcal{P}+\mathcal{B}_{2}\mathcal{Z})^{*}+(\mathcal{A}\mathcal{P}+\mathcal{B}_{2}\mathcal{Z})\mathcal{T}^{*}\end{bmatrix}\preccurlyeq
0$ (11.9)

so that, for any solution $(\gamma,\mathcal{P},\mathcal{Z})$ to this problem,
letting $\mathcal{K}:=\mathcal{Z}\mathcal{P}^{-1}$, we have
$\|z\|\leq\gamma\|w\|$.

In Chapter 7, we showed in detail how the LPI (11.6) can be declared for a PIE
of the form (11.6). The code declaring and solving this LPI can also be found
in the demo file “Hinf_optimal_control_DEMO”. Here we use pre-defined
executive files to solve the above LPI, which is a simple function call as
shown below

     settings = lpisettings(’heavy’);
     [prog, Kval, gam_val] = PIETOOLS_Hinf_control(PIE, settings);

The executive, in this particular case, returns a value Kval of the operator
$\mathcal{K}$ such that the $L_{2}$-gain $\frac{\|z\|_{L_{2}}}{\|w\|_{L_{2}}}$
is bounded from above by gam_val. Next, we construct the closed loop PIE using
the function closedLoop_PIE,

        PIE_CL = closedLoopPIE(PIE,Kval);

Note, closed loop system can be constructed manually, similar to the method
used in 11.5.

Then, we can use PIESIM to simulate the evolution of the PDE state
$\mathbf{x}(t)=\bigl{(}\mathcal{T}\mathbf{v}\bigr{)}(t)$ associated to the
system defined by PIE_CL. As described before we can set initial condition and
disturbance using the commands below

     % Declare initial conditions for the state components of the PIE
     syms st sx real
     uinput.ic.PDE = -10*sx;    % IC PIE
     uinput.w = exp(-st);   % disturbance

Note that, since we will be simulating the PIE directly, the initial
conditions uinput.ic.PDE will also correspond to the initial values of the PIE
state $\mathbf{v}(0,s)$, not to those of the PDE state. In this example, we
have used multiple initial conditions and disturbance inputs. For brevity, the
code corresponding to all the different initial conditions and disturbances is
not presented.

Next, we set the parameters related to the numerical scheme and simulation as
shown below

     % Set options for the discretization and simulation:
     opts.plot = ’no’;   % Do not plot the final solution
     opts.N = 8;         % Expand using 8 Chebyshev polynomials
     opts.tf = 1;        % Simulate up to t = 1;
     opts.dt = 1e-3;     % Use time step of 10^-3
     opts.intScheme = 1; % Time-step using Backward Differentiation Formula (BDF)
     ndiff = [0,0,1];    % The PDE state involves 2 twice differentiable states

These parameters are same as the ones described in 11.5, and hence the
description is omitted here.

We run multiple simulations for both open (without controller) and closed loop
PIEs (with controller) the results of which are plotted in Figure 11.2 and
Figure 11.3, respectively.

Figure 11.2: Simulated value of PDE state $\mathbf{x}(t,s)$ at several grid
points $s\in[0,1]$, without controller input, for different initial conditions
and a bounded $L_{2}$ disturbance $w$. Figure 11.3: Simulated value of PDE
state $\mathbf{x}(t,s)$ at several grid points $s\in[0,1]$, with
$H_{\infty}$-optimal state feedback controller, for different initial
conditions and a bounded $L_{2}$ disturbance $w$. $L_{2}$-gain $=1.6205$.

Lastly, for the initial condition $\mathbf{v}(0,s)=-10s$, we plot the
regulated output $z(t)$ when subjected to two different disturbance inputs,
namely, $w(t)=\exp(-t)$ and $w(t)=\frac{\sin(\pi t)}{t}$. This can be obtained
by numerical integration of the solution as shown in the code below, which can
they be plotted to obtain Figure 11.4.

    tval = linspace(0,1,2000); % grid points in time
    XX = linspace(0,1,20); % grid points in space
    w1_tval = subs(sin(pi*st)./(st+eps),tval);
    w2_tval = subs(exp(-st),tval);
    z_quadrature = double(subs(C1.Q1(1,1),s,XX));
    k_quadrature = double(subs(Kval.Q1,s,XX));
    ZZ1 = trapz(z_quadrature,x_CL_a)+double(w1_tval); %z1 for w1 disturbance
    ZZ2 = trapz(z_quadrature,x_CL_aa)+double(w2_tval);%z1 for w2 disturbance
    ZZ3 = trapz(k_quadrature,x_CL_a); %u(t) for w1 disturbance
    ZZ4 = trapz(k_quadrature,x_CL_aa); %u(t) for w1 disturbance

Figure 11.4: Simulated value of $Z(t)$, with $H_{\infty}$-optimal state
feedback controller, for different bounded $L_{2}$ disturbances, $w$.
$L_{2}$-gain $=1.6205$.

The full code for designing the optimal controller, simulating the PDE state,
and plotting the results has been included in PIETOOLS as the demo file
“Hinf_optimal_control_DEMO”

####  11.7 DEMO 7: Observer-based Controller design and simulation for PDEs

In Sections 11.5 and 11.6, we showed the procedure to find optimal estimator
and controller for a PDE. The controller was designed based on the information
of the complete PDE state, which, in practice, is unrealistic. So, in this
section, we use the observer and controller obtained in the previous sections
and couple them by changing the control law to be dependent on the estimated
state instead of PDE state, i.e., we set $u=\mathcal{K}\hat{\mathbf{v}}$,
where $\mathcal{K}$ is the control gains obtained in 11.6 and
$\hat{\mathbf{v}}$ is the state estimate from the observer found in 11.5.

For demonstration, we reuse the unstable reaction-diffusion PDE, with output
$z$, disturbance $w$ and control input $u$

$\displaystyle\dot{\mathbf{x}}(t,s)$
$\displaystyle=\partial_{s}^{2}\mathbf{x}(t,s)+4\mathbf{x}(t,s)+sw(t)+su(t),$
$\displaystyle s$ $\displaystyle\in[0,1],$ with BCs $\displaystyle 0$
$\displaystyle=\mathbf{x}(t,0)=\partial_{s}\mathbf{x}(t,1),$ and outputs
$\displaystyle z(t)$
$\displaystyle=\begin{bmatrix}\int_{0}^{1}\mathbf{x}(t,s)ds+w(t)\\\
u(t)\end{bmatrix}$ $\displaystyle y(t)$ $\displaystyle=\mathbf{x}(t,1).$
(11.10)

Recall, this PDE can be defined using the code

    pvar s t;
    lam = 4;
    PDE = sys();
    x = state(’pde’);   w = state(’in’);
    z = state(’out’, 2);   u = state(’in’); y = state(’out);
    eqs = [diff(x,t) == diff(x,s,2) + lam*x + s*w + s*u;
        z == [int(x,s,[0,1]) + w; u];
        y == subs(x,s,1);
        subs(x,s,0)==0;
        subs(diff(x,s),s,1)==0];
    PDE = addequation(PDE,eqs);
    PDE = setControl(PDE,u);
    PDE = setObserve(PDE,u);
    display_PDE(PDE);

Next, we convert the PDE system to a PIE by using the following code and
extract relevant PI operators for easier access:

    PIE = convert(PDE,’pie’);       PIE = PIE.params;
    T = PIE.T;
    A = PIE.A;      C1 = PIE.C1;    B2 = PIE.B2;
    B1 = PIE.B1;    D11 = PIE.D11;  D12 = PIE.D12;
    C2 = PIE.C2;    D21 = PIE.D21;  D22 = PIE.D22;

Details on observer design and controller design can be found in Sections 11.5
and 11.6 respectively, and therefore, will not be repeated here. For now, we
proceed assuming the observer gains $\mathcal{L}$ and the controller gains
$\mathcal{K}$ are known.

Recall, from Equations (11.5) and (11.6), the combined system is given by

$\displaystyle\bigl{(}\mathcal{T}\dot{\mathbf{v}}\bigr{)}(t)$
$\displaystyle=\bigl{(}\mathcal{A}\mathbf{v}\bigr{)}(t)+\bigl{(}\mathcal{B}_{1}w\bigr{)}(t)+\bigl{(}\mathcal{B}_{2}u\bigr{)}(t),$
$\displaystyle z(t)$
$\displaystyle=\bigl{(}\mathcal{C}_{1}\mathbf{v}\bigr{)}(t)+\bigl{(}\mathcal{D}_{11}w)(t)+\bigl{(}\mathcal{D}_{12}u)(t),$
$\displaystyle\mathcal{T}\dot{\hat{\mathbf{v}}}(t)$
$\displaystyle=\mathcal{A}{\mathbf{\hat{v}}}(t)+\mathcal{L}\bigl{(}\mathcal{C}_{2}\mathbf{v}(t)-\mathcal{C}_{2}\mathbf{\hat{v}}(t)\bigr{)},$
$\displaystyle\hat{z}(t)$ $\displaystyle=\mathcal{C}_{1}\mathbf{\hat{v}}(t).$
(11.11)

Using the observer-controller coupling, $u=\mathcal{K}\hat{\mathbf{v}}$, we
get the closed loop PIE

$\displaystyle\left(\begin{bmatrix}\mathcal{T}&0\\\
0&\mathcal{T}\end{bmatrix}\begin{bmatrix}\dot{\mathbf{v}}\\\
\dot{\hat{\mathbf{v}}}\end{bmatrix}\right)(t,s)$
$\displaystyle=\left(\begin{bmatrix}\mathcal{A}&\mathcal{B}_{2}\mathcal{K}\\\
-\mathcal{L}\mathcal{C}_{2}&\mathcal{A}+\mathcal{L}\mathcal{C}_{2}\end{bmatrix}\begin{bmatrix}\mathbf{v}\\\
\hat{\mathbf{v}}\end{bmatrix}\right)(t,s)+\left(\begin{bmatrix}\mathcal{B}_{1}\\\
\mathcal{L}\mathcal{D}_{21}\end{bmatrix}w\right)(t)$
$\displaystyle\begin{bmatrix}z\\\ \hat{z}\end{bmatrix}(t)$
$\displaystyle=\left(\begin{bmatrix}\mathcal{C}_{1}&\mathcal{D}_{12}\mathcal{K}\\\
0&\mathcal{C}_{1}\end{bmatrix}\begin{bmatrix}\mathbf{v}\\\
\hat{\mathbf{v}}\end{bmatrix}\right)(t)+\left(\begin{bmatrix}\mathcal{D}_{11}\\\
0\end{bmatrix}w\right)(t).$ (11.12)

We can construct the above closed loop system using the code

    T_CL = [T, 0*T; 0*T, T];
    A_CL = [A, B2*Kval; -Lval*C2, A+Lval*C2];   B_CL = [B1; Lval*D21];
    C_CL = [C1, D12*Kval; 0*C1, C1];            D_CL = [D11; 0*D11];
    PIE_CL = pie_struct();
    PIE_CL.vars = PIE.vars;
    PIE_CL.dom = PIE.dom;
    PIE_CL.T = T_CL;
    PIE_CL.A = A_CL;        PIE_CL.B1 = B_CL;
    PIE_CL.C1 = C_CL;       PIE_CL.D11 = D_CL;
    PIE_CL = initialize(PIE_CL);

Once, we have the closed loop PIE object, we can use PIESIM to simulate the
system with different initial conditions and disturbances. Following the
process presented in Section 11.6 and 11.5, we use the initial condition
$\mathbf{v}(0,s)=-10s$ and $\hat{\mathbf{v}}(0,s)=0$. Using the same options
and disturbance as before, we simulate the open loop system (without control
or observer) using the code

    syms st sx real
    opts.plot = ’no’;   % Do not plot the final solution
    opts.N = 8;         % Expand using 8 Chebyshev polynomials
    opts.tf = 1;        % Simulate up to t = 1;
    opts.dt = 1e-3;     % Use time step of 10^-3
    opts.intScheme=1;   % Time-step using Backward Differentiation Formula (BDF)
    ndiff = [0,0,1];    % The PDE state involves 1 second order differentiable state variables

    % Simulate the solution to the PIE without controller for different IC.
    uinput.ic.PDE = [-10*sx];   % IC PIE
    uinput.w = exp(-st); % disturbance
    [solution_OL,grid] = PIESIM(PIE,opts,uinput,ndiff);

Then, we can access the time-varying PDE solution from solution_OL for post-
processing. In Figure 11.5 we post the solution at different time values to
show that the open loop system is unstable.

Figure 11.5: Simulated value of PDE state $\mathbf{x}(t,s)$ at several grid
points $s\in[0,1]$, without controller input, with some initial condition and
a bounded $L_{2}$ disturbance $w$.

Likewise, we can simulate and extract the closed loop system solution using
the code

    % Simulate the solution to the PIE with controller for different IC and disturbance.
    ndiff = [0,0,2];
    uinput.ic.PDE = [-10*sx; 0];    % IC PIE and observed state
    uinput.w = exp(-st);   % disturbance
    [solution_CL_a,grid] = PIESIM(PIE_CL,opts,uinput,ndiff);
    uinput.ic.PDE = [sin(sx*pi/2); 0];    % IC PIE
    uinput.w = 5*sin(pi*st)./(st+eps);   % disturbance
    [solution_CL_b,grid] = PIESIM(PIE_CL,opts,uinput,ndiff);
    % Extract actual solution at each time step.
    tval = solution_CL_a.timedep.dtime;
    x_CL_a = reshape(solution_CL_a.timedep.pde(:,1,:),opts.N+1,[]);
    hatx_CL_a = reshape(solution_CL_a.timedep.pde(:,2,:),opts.N+1,[]);
    x_CL_b = reshape(solution_CL_b.timedep.pde(:,1,:),opts.N+1,[]);
    hatx_CL_b = reshape(solution_CL_b.timedep.pde(:,2,:),opts.N+1,[]);

where, note that ndiff is changed to $[0,0,2]$ because the closed loop PIE
system has two distributed states corresponding to two PDE states (PDE state
and observer state) that are twice differentiable in space. Using the
solutions stored in x_CL_a(b), we can plot the response of the closed loop
system for different initial conditions and disturbances as shown below in
Figure 11.6.

Figure 11.6: Simulated value of PDE state $\mathbf{x}(t,s)$ at several grid
points $s\in[0,1]$, with $H_{\infty}$-optimal state feedback controller, for
different initial conditions and bounded $L_{2}$ disturbances $w$.
$L_{2}$-gain $=1.6205$ for controller, $L_{2}$-gain $=1.005$ for observer.

Lastly, for the initial condition $\mathbf{v}(0,s)=-10s$ and
$\hat{\mathbf{v}}(0)=0$, we plot (see Figure 11.7 the regulated output $z(t)$
when subjected to two different disturbance inputs, namely, $w(t)=\exp(-t)$
and $w(t)=5\frac{\sin(\pi t)}{t}$, where $z$ is obtained by numerical
integration of the solution as described in Section 11.6.

Figure 11.7: Simulated value of $Z(t)$, with $H_{\infty}$-optimal state
feedback controller, for different bounded $L_{2}$ disturbances, $w$.

###  Chapter 12 Libraries of PDE and TDS Examples in PIETOOLS

In Chapters 4, 8 and 9, we have shown how partial differential equations and
time-delay systems can be declared in PIETOOLS through different input
formats. To help get started with each of these input formats, PIETOOLS
includes a variety of example pre-defined PDE and TDS systems, including
common examples and particular models from the literature. These examples are
collected in the PIETOOLS_examples folder, and can be accessed calling the
function examples_PDE_library_PIETOOLS and the scripts
examples_DDE_library_PIETOOLS and examples_NDSDDF_library_PIETOOLS. In this
Chapter, we illustrate how this works, focusing on PDE examples in Section
12.1, and TDS examples in Section 12.2.

####  12.1 A Library of PDE Example Problems

To help get started with analysing and simulating PDEs in PIETOOLS, a variety
of PDE models have been included in the directory
PIETOOLS_examples/Examples_Library. These systems include common PDE models,
as well as examples from the literature, and are defined through separate
MATLAB functions. Each of these functions takes two arguments

  1. 1.

GUI: A binary index (0 or 1) indicating whether the example should be opened
in the graphical user interface;

  2. 2.

params: A string for specifying allowed parameters in the system.

For example, calling help PIETOOLS_PDE_Ex_Reaction_Diffusion_Eq, PIETOOLS
indicates that this function declares a reaction-diffusion PDE

$\displaystyle\dot{\mathbf{x}}(t,s)$
$\displaystyle=\lambda\mathbf{x}(t,s)+\partial_{s}^{2}\mathbf{x}(t,s),$
$\displaystyle s$ $\displaystyle\in[0,1]\mathbf{x}(t,0)$
$\displaystyle=\mathbf{x}(t,1)=0,$

where the value of the parameter $\lambda$ can be set. Then, calling

     >> PDE = PIETOOLS_PDE_Ex_Reaction_Diffusion_Eq(0,{’lam=10;’})

we obtain a pde_struct object PDE representing the reaction-diffusion equation
with $\lambda=10$. Calling

     >> PDE = PIETOOLS_PDE_Ex_Reaction_Diffusion_Eq(1,{’lam=10;’})

The PDE will also be loaded in the GUI, though a default value of
$\lambda=9.86$ will be used, as the GUI will always load a pre-defined file,
which cannot be adjusted from the command line.

To simplify the process of extracting PDE examples, PIETOOLS includes a
function examples _PDE_library_PIETOOLS(). In this function, each of the pre-
defined PDEs is assigned an index, allowing desired PDEs to be extracted by
calling examples_PDE_library_PIETOOLS with the associated index. For example,
scrolling through this function we find that the reaction-diffusiong equation
is the fifth system in the list, and therefore, we can obtain a pde_struct
object defining this system by calling the library function with argument “5”,
returning

     >> PDE = examples_PDE_library_PIETOOLS(5);
      --- Extracting ODE-PDE example 5 ---

      (d/dt) x(t,s) = 9.86 * x(t,s) + (d^2/ds^2) x(t,s);
      0 = x(t,0);
      0 = x(t,1);

     Would you like to run the executive associated to this problem? (y/n)
     -->

We note that the function asks whether an executive should be run for the
considered PDE. This is because, for each of the PDE examples, an associated
LPI problem has also been declared, matching one of the executive files (see
also Chapter 13). For the reaction-diffusion equation, the proposed executive
is the PIETOOLS_stability function, testing stability of the PDE. Entering yes
in the command line window, this executive will be automatically run, whilst
entering no will stop the function, and just return the pde_struct object PDE.

Using the examples_PDE_library_PIETOOLS function, parameters in the PDE can
also be adjusted, calling e.g.

     >> PDE = examples_PDE_library_PIETOOLS(5,’lam=10;’);
      --- Extracting ODE-PDE example 5 ---

      (d/dt) x(t,s) = 10 * x(t,s) + (d^2/ds^2) x(t,s);
      0 = x(t,0);
      0 = x(t,1);

Similarly, if multiple parameters can be specified, we specify each of these
parameters separately. For example, we note that example 7 corresponds to a
PDE

$\displaystyle\dot{\mathbf{x}}(t,s)$
$\displaystyle=c(s)\mathbf{x}(t,s)+b(s)\partial_{s}\mathbf{x}(t,s)+a(s)\partial_{s}^{2}\mathbf{x}(t,s),$
$\displaystyle s$ $\displaystyle\in[0,1],\mathbf{x}(t,0)$
$\displaystyle=\partial_{s}\mathbf{x}(t,1)=0,$

where the values of the functions $a(s)$, $b(s)$ and $c(s)$ can be specified.
As such, we can declare this PDE for $a=1$, $b=2$ and $c=3$ by calling

     >> PDE = examples_PDE_library_PIETOOLS(7,’a=1;’,’b=2;’,’c=3;’);
      --- Extracting ODE-PDE example 7 ---

      (d/dt) x(t,s) = 3 * x(t,s) + 2 * (d/ds) x(t,s) + (d^2/ds^2) x(t,s);
      0 = x(t,0);
      0 = (d/ds) x(t,1);

Finally, we can also open the PDE in the GUI by calling

     >> examples_PDE_library_PIETOOLS(7,’GUI’);
      --- Extracting ODE-PDE example 7 ---

or extract the PDE as a pde_struct (terms-format) and open it in the GUI by
calling

     >> PDE = examples_PDE_library_PIETOOLS(7,’TERM’,’GUI’);
      --- Extracting ODE-PDE example 7 ---

In each case, the function will still ask whether the executive associated
with the PDE should be run as well. Of course, you can also convert the PDE to
a PIE yourself using convert, and then run any desired executive manually,
assuming this executive makes sense (e.g. there’s no sense in computing an
$H_{\infty}$-gain if your system has no outputs).

####  12.2 Libraries of DDE, NDS, and DDF Examples

Aside from the PDE examples, a list of TDS examples is also included in
PIETOOLS, in DDE, NDS, and DDF format. Unlike the PDE problems, however, these
TDS examples are all not separated into distinct functions, but are divided
over two scripts: examples_DDE_library_PIETOOLS and
examples_NDSDDF_library_PIETOOLS. In each of these scripts, most examples are
commented, and only one example should be uncommented at any time. This
example can then be extracted by calling the script, adding a structure DDE,
NDS or DDF to the MATLAB workspace. To extract a different example, the
desired example must be uncommented, and all other examples must be commented,
at which point the script can be called again to obtain a structure
representing the desired system.

We expect to update the DDE and NDS/DDF example files in a future release to
match the format used for the PDE example library.

#####  12.2.1 DDE Examples

We have compiled a list of 23 DDE numerical examples, grouped into: stability
analysis problems; input-output systems; estimator design problems; and
feedback control problems. These examples are drawn from the literature and
citations are used to indicate the source of each example. For each group, the
relevant flags have been included to indicate which executive mode should be
called after the example has been loaded.

#####  12.2.2 NDS and DDF Examples

There are relatively few DDFs which do not arise from a DDE or NDS. Hence, we
have combined the DDF and NDS example libraries into the script
examples_NDSDDF_library_PIETOOLS. The Neutral Type systems are listed first,
and currently consist only of stability analysis problems - of which we
include 13. As for the DDE case, the library is a script, so the user must
uncomment the desired example and call the script from the root file or
command window. For the NDS problems, after calling the example library, in
order to convert the NDS to a DDF or PIE, the user can use the following
commands:

     >> NDS = initialize\_PIETOOLS\_NDS(NDS);
     >> DDF = convert\_PIETOOLS\_NDS(NDS,’ddf’);
     >> DDF = minimize\_PIETOOLS\_DDF(DDF);
     >> PIE = convert\_PIETOOLS\_DDF(DDF,’pie’);

In contrast to the NDS case, we only include 3 DDF examples. The first two are
difference equations which cannot be represented in either the NDS or DDE
format. The third is a network control problem, which is also included in the
DDE library in DDE format.

###  Chapter 13 Standard Applications of LPI Programming

In Chapter 7, we showed how general LPI optimization programs can be declared
and solved in PIETOOLS. In this chapter, we provide several applications of
LPI programming for analysis, estimation, and control of PIEs. Recall that
such PIEs take the form

$\displaystyle\mathcal{T}\dot{\mathbf{v}}(t)+\mathcal{T}_{w}\dot{w}(t)+\mathcal{T}_{u}\dot{u}(t)$
$\displaystyle=\mathcal{A}\mathbf{v}(t)+\mathcal{B}_{1}w(t)+\mathcal{B}_{2}u(t),\quad\mathbf{v}(0)=\mathbf{v}_{\text{I}}$
$\displaystyle z(t)$
$\displaystyle=\mathcal{C}_{1}\mathbf{v}(t)+\mathcal{D}_{11}w(t)+\mathcal{D}_{12}u(t),$
$\displaystyle y(t)$
$\displaystyle=\mathcal{C}_{2}\mathbf{v}(t)+\mathcal{D}_{21}w(t)+\mathcal{D}_{22}u(t),$
(13.1)

where $\mathbf{v}=\left[\begin{smallmatrix}v_{0}\\\ \mathbf{v}_{1}\\\
\mathbf{v}_{2}\\\
\mathbf{v}_{3}\end{smallmatrix}\right]\in\left[\begin{smallmatrix}\mathbb{R}^{n_{0}}\\\
L_{2}^{n_{1}}[a,b]\\\ L_{2}^{n_{2}}[c,d]\\\
L_{2}^{n_{3}}[[a,b]\times[c,d]]\end{smallmatrix}\right]$, and where
$\mathcal{T}$ through $\mathcal{D}_{22}$ are all PI operators. In Section
13.1, we provide several LPIs for stability analysis and $H_{\infty}$-gain
estimation of such PIEs, setting $u=0$. Then, in Section 13.2, we present an
LPI for $H_{\infty}$-optimal estimation of PIEs of the form (13), followed by
an LPI for $H_{\infty}$-optimal full-state feedback control in Section 13.3.
We note that, almost all of these LPIs have already been implemented as
executive functions in PIETOOLS, and we will refer to these executives when
applicable.

####  13.1 LPIs for Analysis of PIEs

Using LPIs, several properties of a PIE as in Equation (13) may be tested, as
listed in this section. In particular, the LPIs listed below are extensions of
classical results used in analysis of ODEs using LMIs. For most of these LPIs,
PIETOOLS includes an executive function that may be run to solve it for a
given PIE.

#####  13.1.1 Operator Norm

For a PI operator $\mathcal{T}$, an upper bound $\sqrt{\gamma}$ on the
operator norm $\|\mathcal{T}\|$ can be computed by solving the LPI

$\displaystyle\min\limits_{\gamma,\mathcal{P}}\leavevmode\nobreak\
\leavevmode\nobreak\ \gamma$
$\displaystyle\mathcal{T}^{*}\mathcal{T}-\gamma\preccurlyeq 0$ (13.2)

This LPI has not been implemented as an executive in PIETOOLS, but has been
implemented in the demo file volterra_operator_norm_DEMO (see also Section
11.2).

#####  13.1.2 Stability

For given PI operators $\mathcal{T}$ and $\mathcal{A}$, stability of the PIE

$\displaystyle\mathcal{T}\dot{\mathbf{v}}(t)=\mathcal{A}\mathbf{v}(t)$

can be tested by solving the LPI

$\displaystyle\mathcal{P}\succ 0$
$\displaystyle\mathcal{T}^{*}\mathcal{P}\mathcal{A}+\mathcal{A}^{*}\mathcal{P}\mathcal{T}\preccurlyeq
0$ . (13.3)

If there exists a PI operator $\mathcal{P}$ such that this LPI is feasible,
then the PIE is stable. Given a structure PIE, this LPI may be solved for the
associated PIE by calling

     >> [prog, Pop] = PIETOOLS_stability(PIE, settings);

Here prog will be an SOS program structure describing the solved problem and
Pop will be a dopvar object describing the (unsolved) decision operator
$\mathcal{P}$ from which the solved operator can be derived using

     >> Pop = getsol_lpivar(prog,Pop);

See Chapter 7 for more information on the operation of this function and the
settings input.

#####  13.1.3 Dual Stability

For given PI operators $\mathcal{T}$ and $\mathcal{A}$, stability of the PIE

$\displaystyle\mathcal{T}\dot{\mathbf{v}}(t)=\mathcal{A}\mathbf{v}(t)$

can also be tested by solving the LPI

$\displaystyle\mathcal{P}\succ 0$
$\displaystyle\mathcal{T}\mathcal{P}\mathcal{A}^{*}+\mathcal{A}\mathcal{P}\mathcal{T}^{*}\preccurlyeq
0$ (13.4)

If there exists a PI operator $\mathcal{P}$ such that this LPI is feasible,
then the PIE is stable. Given a structure PIE, this LPI may be solved for the
associated PIE by calling

     >> [prog, Pop] = PIETOOLS_stability_dual(PIE, settings);

Here prog will be an SOS program structure describing the solved problem and
Pop will be a dopvar object describing the (unsolved) decision operator
$\mathcal{P}$.

#####  13.1.4 Input-Output Gain

Consider a system of the form

$\displaystyle\mathcal{T}\dot{\mathbf{v}}(t)$
$\displaystyle=\mathcal{A}\mathbf{v}(t)+\mathcal{B}_{1}w(t),\quad\mathbf{v}(0)=\mathbf{0}$
$\displaystyle z(t)$
$\displaystyle=\mathcal{C}_{1}\mathbf{v}(t)+\mathcal{D}_{11}w(t),$ (13.5)

where $w\in L_{2}^{n_{w}}[0,\infty)$. Then, $z\in L_{2}^{n_{z}}[0,\infty)$,
and $\|z\|_{L_{2}}\leq{\gamma}\|w\|_{L_{2}}$, if the following LPI is feasible

$\displaystyle\min\limits_{\gamma,\mathcal{P}}\leavevmode\nobreak\
\leavevmode\nobreak\ \gamma$ $\displaystyle\mathcal{P}\succ 0$
$\displaystyle\begin{bmatrix}-\gamma
I&\mathcal{D}_{11}^{*}&\mathcal{B}_{1}^{*}\mathcal{PT}\\\ (\cdot)^{*}&-\gamma
I&\mathcal{C}_{1}\\\
(\cdot)^{*}&(\cdot)^{*}&\mathcal{T}^{*}\mathcal{P}\mathcal{A}+\mathcal{A}^{*}\mathcal{P}\mathcal{T}\end{bmatrix}\preccurlyeq
0$ (13.6)

Given a structure PIE, this LPI may be solved for the associated PIE by
calling

     >> [prog, Pop, gam] = PIETOOLS_Hinf_gain(PIE, settings);

Here prog will be an SOS program structure describing the solved problem, and
gam will be the smallest value of $\gamma$ for which the LPI was found to be
feasible, offering a bound on the $L_{2}$-gain from $w$ to $z$ of the system.
The output Pop will be a dopvar object describing the (unsolved) decision
operator $\mathcal{P}$.

#####  13.1.5 Dual Input-Output Gain

For a System (13.1.4) with $w\in L_{2}^{n_{w}}[0,\infty)$, an upper bound
$\gamma$ on the $L_{2}$-gain from $w$ to $z$ can also be obtained by solving
the LPI

$\displaystyle\min\limits_{\gamma,\mathcal{P}}\leavevmode\nobreak\
\leavevmode\nobreak\ \gamma$ $\displaystyle\mathcal{P}\succ 0$
$\displaystyle\begin{bmatrix}-\gamma
I&\mathcal{D}_{11}&\mathcal{T}\mathcal{P}\mathcal{C}_{1}\\\
(\cdot)^{*}&-\gamma I&\mathcal{B}_{1}^{*}\\\
(\cdot)^{*}&(\cdot)^{*}&\mathcal{T}\mathcal{P}\mathcal{A}^{*}+\mathcal{A}\mathcal{P}\mathcal{T}^{*}\end{bmatrix}\preccurlyeq
0$ (13.7)

Given a structure PIE, this LPI may be solved for the associated PIE by
calling

     >> [prog, Pop, gam] = PIETOOLS_Hinf_gain_dual(PIE, settings);

Here prog will be an SOS program structure describing the solved problem, and
gam will be the found optimal value for $\gamma$. The output Pop will be a
dopvar object describing the (unsolved) decision operator $\mathcal{P}$.

#####  13.1.6 Positive Real Lemma

For a PIE of the form

$\displaystyle\mathcal{T}\dot{\mathbf{v}}(t)$
$\displaystyle=\mathcal{A}\mathbf{v}(t)+\mathcal{B}_{1}w(t),\quad\mathbf{v}(0)=\mathbf{0}$
$\displaystyle z(t)$
$\displaystyle=\mathcal{C}_{1}\mathbf{v}(t)+\mathcal{D}_{11}w(t),$

we can test whether the system is passive by solving the LPI

$\displaystyle\mathcal{P}\succ 0$
$\displaystyle\begin{bmatrix}-\mathcal{D}_{11}^{*}-\mathcal{D}_{11}&\mathcal{B}_{1}^{*}\mathcal{PT}-\mathcal{C}_{1}\\\
(\cdot)^{*}&\mathcal{T}^{*}\mathcal{P}\mathcal{A}+\mathcal{A}^{*}\mathcal{P}\mathcal{T}\end{bmatrix}\preccurlyeq
0$ (13.8)

If there exists a PI operator $\mathcal{P}$ such that this LPI is feasible,
then the system is passive. Note that this LPI has not been implemented as an
executive in PIETOOLS.

####  13.2 LPIs for Optimal Estimation of PIEs

For the following PIE

$\displaystyle\mathcal{T}\dot{\mathbf{v}}(t)+\mathcal{T}_{w}\dot{w}(t)$
$\displaystyle=\mathcal{A}\mathbf{v}(t)+\mathcal{B}_{1}w(t),$ $\displaystyle
z(t)$ $\displaystyle=\mathcal{C}_{1}\mathbf{v}(t)+\mathcal{D}_{11}w(t),$
$\displaystyle y(t)$
$\displaystyle=\mathcal{C}_{2}\mathbf{v}(t)+\mathcal{D}_{21}w(t),$ (13.9)

a state estimator has the following structure:

$\displaystyle\mathcal{T}\dot{\mathbf{\hat{v}}}(t)$
$\displaystyle=\mathcal{A}\mathbf{\hat{v}}(t)+\mathcal{L}(\hat{y}(t)-y(t)),$
$\displaystyle\hat{z}(t)$ $\displaystyle=\mathcal{C}_{1}\mathbf{\hat{v}}(t)$
$\displaystyle\hat{y}(t)$ $\displaystyle=\mathcal{C}_{2}\mathbf{\hat{v}}(t).$
(13.10)

so that the errors $\mathbf{e}:=\hat{\mathbf{v}}-\mathbf{v}$ and
$\tilde{z}:=\hat{z}-z$ in respectively the state and regulated output
estimates satisfy

$\displaystyle\mathcal{T}\dot{\mathbf{e}}-\mathcal{T}_{w}w(t)$
$\displaystyle=(\mathcal{A}+\mathcal{L}\mathcal{C}_{2})\mathbf{e}(t)-(\mathcal{B}_{1}+\mathcal{L}\mathcal{D}_{21})w(t),$
$\displaystyle\tilde{z}(t)$
$\displaystyle=\mathcal{C}_{1}\mathbf{e}(t)-D_{11}w(t)$

The $H_{\infty}$-optimal estimation problem amounts to synthesizing
$\mathcal{L}$ such that the estimation error $\tilde{z}:={\hat{z}-z}$ admits
$\|\tilde{z}\|\leq\gamma\|w\|$ for a particular $\gamma>0$. To establish such
an estimator, we can solve the LPI

$\displaystyle\min\limits_{\gamma,\mathcal{P},\mathcal{Z}}\leavevmode\nobreak\
\leavevmode\nobreak\ \gamma$ $\displaystyle\mathcal{P}\succ 0$
$\displaystyle\begin{bmatrix}\mathcal{T}_{w}^{*}(\mathcal{P}\mathcal{B}_{1}+\mathcal{Z}\mathcal{D}_{21})+(\cdot)^{*}&0&(\cdot)^{*}\\\
0&0&0\\\
-(\mathcal{P}\mathcal{A}+\mathcal{Z}\mathcal{C}_{2})^{*}\mathcal{T}_{w}&0&0\end{bmatrix}+\begin{bmatrix}-\gamma
I&-\mathcal{D}_{11}^{\top}&-(\mathcal{P}\mathcal{B}_{1}+\mathcal{Z}\mathcal{D}_{21})^{*}\mathcal{T}\\\
(\cdot)^{*}&-\gamma I&\mathcal{C}_{1}\\\
(\cdot)^{*}&(\cdot)^{*}&(\mathcal{P}\mathcal{A}+\mathcal{Z}\mathcal{C}_{2})^{*}\mathcal{T}+(\cdot)^{*}\end{bmatrix}\preccurlyeq
0$ (13.11)

Then, if this LPI is feasible for some $\gamma>0$ and PI operators
$\mathcal{P}$ and $\mathcal{Z}$, then, letting
$\mathcal{L}:=\mathcal{P}^{-1}\mathcal{Z}$, the estimation error will satisfy
$\|z-\hat{z}\|\leq\gamma\|w\|$. Given a structure PIE, this LPI may be solved
for the associated PIE by calling

     >> [prog, Lop, gam, Pop, Zop] = PIETOOLS_Hinf_estimator(PIE, settings);

Here prog will be an SOS program structure describing the solved problem, gam
will be the found optimal value for $\gamma$, and Lop will be an opvar object
describing the optimal estimator $\mathcal{L}$. Outputs Pop and Zop will be
dopvar objects describing the (unsolved) decision operators $\mathcal{P}$ and
$\mathcal{Z}$, from which the solved operators can be derived using

     >> Pop = getsol_lpivar(prog,Pop);        Zop = getsol_lpivar(prog,Zop);

See Chapter 7 for more information on the operation of this function and the
settings input.

####  13.3 LPIs for Optimal Control of PIEs

In this section, we discuss the synthesis of $H_{\infty}$ optimal control of a
PIE of the form

$\displaystyle\mathcal{T}\dot{\mathbf{v}}(t)$
$\displaystyle=\mathcal{A}\mathbf{v}(t)+\mathcal{B}_{1}w(t)+\mathcal{B}_{2}u(t),\qquad\mathbf{v}(0)=\mathbf{v}_{0}$
$\displaystyle z(t)$
$\displaystyle=\mathcal{C}_{1}\mathbf{v}(t)+\mathcal{D}_{11}w(t)+\mathcal{D}_{12}u(t),$
(13.12)

where $w\in L_{2}[0,\infty)$. The problem of synthesizing an
$H_{\infty}$-optimal controller amounts to determining a PIE operator
$\mathcal{K}$ such that, using the full-state feedback law
$u(t)=\mathcal{K}\mathbf{v}(t)$, the regulated output ${z}$ admits
$\|z\|_{L_{2}}\leq\gamma\|w\|_{L_{2}}$ for a particular $\gamma>0$. To
establish such a controller, we can solve the LPI

$\displaystyle\min\limits_{\gamma,\mathcal{P},\mathcal{Z}}\leavevmode\nobreak\
\leavevmode\nobreak\ \gamma$ $\displaystyle\mathcal{P}\succ 0$ $\displaystyle
dj$ $\displaystyle\begin{bmatrix}-\gamma
I&\mathcal{D}_{11}&(\mathcal{C}_{1}\mathcal{P}+D_{12}\mathcal{Z})\mathcal{T}^{*}\\\
\mathcal{D}_{11}^{*}&-\gamma I&\mathcal{B}_{1}^{*}\\\
()^{*}&\mathcal{B}_{1}&()^{*}+\left(\mathcal{AP}+\mathcal{B}_{2}\mathcal{Z}\right)\mathcal{T}^{*}\end{bmatrix}\preccurlyeq
0$ (13.13)

If this LPI is feasible for some $\gamma>0$ and PI operators $\mathcal{P}$ and
$\mathcal{Z}$, then, letting $\mathcal{K}:=\mathcal{Z}\mathcal{P}^{-1}$, the
$L_{2}$-gain for the controlled system with $u=\mathcal{K}\mathbf{v}$ will be
such that $\|z\|\leq\gamma\|w\|$. Given a structure PIE, this LPI may be
solved for the associated PIE by calling

     >> [prog, Kop, gam, Pop, Zop] = PIETOOLS_Hinf_control(PIE, settings);

Here prog will be an SOS program structure describing the solved problem, gam
will be the found optimal value for $\gamma$, and Kop will be an opvar object
describing the optimal feedback $\mathcal{K}$. Outputs Pop and Zop will be
dopvar objects describing the (unsolved) decision operators $\mathcal{P}$ and
$\mathcal{Z}$, from which the solved operators can be derived using

     >> Pop = getsol_lpivar(prog,Pop);        Zop = getsol_lpivar(prog,Zop);

See Chapter 7 for more information on the operation of this function and the
settings input.

##  Part IV Appendices

###  Appendix A PI Operators and their Properties

In this appendix, we discuss in a bit more detail the crucial properties of PI
operators that PIETOOLS relies on for implementation and analysis of PIEs. For
this, in Section A.1, we first recap the definitions of PI operators as
presented in Chapter 5, also introducing some notation that we will continue
to use throughout the appendix. In Section A.2, A.3, and A.4, we show that
respectively the sum, composition, and adjoint of PI operators can be
expressed as PI operators. In Section A.5, and A.6, we then show how the
respectively the inverse of a PI operator, and the composition of a PI
operator with a differential operator can be computed. Finally, in Section
A.7, we show how a cone of positive PI operators can be parameterized by
positive matrices, allowing an LPI constraint $\mathcal{P}\succcurlyeq 0$ to
be posed as an LMI $P\succcurlyeq 0$.

For more information on PI operators, and full proofs of each of the results,
we refer to e.g. [5] [8] (1D) and [4] (2D).

####  A.1 PI Operators on Different Function Spaces

Recall that we denote the space of square integrable functions on a domain
$\Omega$ as $L_{2}[\Omega]$, with inner product

$\displaystyle\left\langle{\mathbf{x}},{\mathbf{y}}\right\rangle_{L_{2}}=\int_{\Omega}\bigl{[}\mathbf{x}(s)\bigr{]}^{T}\mathbf{y}(s)ds.$

In defining the different PI operators, we will restrict ourselves to domains
of 1D or 2D hypercubes. In 1D, such a hypercube is simply an interval $[a,b]$,
for which we define the following PI operator:

######  Definition 1 (3-PI Operator).

For given parameters

$\displaystyle R:=\\{R_{0},R_{1},R_{2}\\}\in\left\\{L_{2}^{m\times
n}[a,b],L_{2}^{m\times n}\bigl{[}[a,b]^{2}\bigr{]},L_{2}^{m\times
n}\bigl{[}[a,b]^{2}\bigr{]}\right\\}=:\mathcal{N}_{1D}^{m\times n}[a,b],$

we define the associated 3-PI operator
$\mathcal{P}[R]:=\mathcal{P}_{\\{R_{0},R_{1},R_{2}\\}}:L_{2}^{n}[a,b]\rightarrow
L_{2}^{m}[a,b]$ as

$\displaystyle\bigl{(}\mathcal{P}[R]\mathbf{x}\bigr{)}(s):=R_{0}(s)\mathbf{x}(s)+\int_{a}^{s}R_{1}(s,\theta)\mathbf{x}(\theta)d\theta+\int_{s}^{b}R_{2}(s,\theta)\mathbf{x}(\theta)d\theta,$
(A.1)

for any $\mathbf{x}\in L_{2}^{n}[a,b]$.

We note that 3-PI operators can be seen as (one possible) generalization of
matrices to infinite dimensional vector spaces. In particular, suppose we have
a matrix $P\in\mathbb{R}^{m\times n}$, which we decompose as $P=D+L+U$, where
$D$ is diagonal, $L$ is strictly lower triangular, and $U$ is strictly upper-
triangular. Then, for any $x\in\mathbb{R}^{n}$, the $i$th element of the
product $Px$ is given by:

$\displaystyle\bigl{(}Px\bigr{)}_{i}$
$\displaystyle=D_{ii}x_{i}+\sum_{j=1}^{i-1}L_{ij}x_{j}+\sum_{j=i+1}^{n}U_{ij}x_{j}.$
Compare this to the value of $\mathcal{P}[R]\mathbf{x}$ at a position
$s\in[a,b]$ for some $\mathbf{x}\in L_{2}[a,b]$ and 3-PI parameters $R$:
$\displaystyle\bigl{(}\mathcal{P}[R]\mathbf{x}\bigr{)}(s)$
$\displaystyle=R_{0}(s)\mathbf{x}(s)+\int_{a}^{s}R_{1}(s,\theta)\mathbf{x}(\theta)d\theta+\int_{s}^{b}R_{2}(s,\theta)\mathbf{x}(\theta)d\theta.$

Replacing row and column indices $(i,j)$ by primary and dummy variables
$(s,\theta)$, and performing integration instead of summation, 3-PI operators
have a structure very similar to that of matrices, wherein we can recognize a
diagonal, lower-triangular, and upper-triangular part. Accordingly, we will
occasionally refer to a PI operator of the form
$\mathcal{P}_{\\{R_{0},0,0\\}}$ as a diagonal 3-PI operator, and to PI
operators of the forms $\mathcal{P}_{\\{0,R_{1},0\\}}$ and
$\mathcal{P}_{\\{0,0,R_{2}\\}}$ as lower- and upper-triangular PI operators
respectively. The similar structure between matrices and PI operator also
ensures that matrix operations such as addition and multiplication are valid
for PI operators as well, as we will discuss in more detail in the next
sections.

To map functions on a domain $[a,b]\times[c,d]\subset\mathbb{R}^{2}$, we also
define the $9$-PI operator:

######  Definition 2 (9-PI Operator).

For given parameters

$\displaystyle R$
$\displaystyle:={\small\left[\\!\\!\begin{array}[]{lll}R_{00}&R_{01}&R_{02}\\\
R_{10}&R_{11}&R_{12}\\\ R_{20}&R_{21}&R_{22}\end{array}\\!\\!\right]}$
$\displaystyle\hskip
14.22636pt\in{\small\left[\begin{array}[]{lll}L_{2}^{m\times
n}\bigl{[}[a,b]\times[c,d]\bigr{]}&L_{2}^{m\times
n}\bigl{[}[a,b]\times[c,d]^{2}\bigr{]}&L_{2}^{m\times
n}\bigl{[}[a,b]\times[c,d]^{2}\bigr{]}\\\ L_{2}^{m\times
n}\bigl{[}[a,b]^{2}\times[c,d]\bigr{]}&L_{2}^{m\times
n}\bigl{[}[a,b]^{2}\times[c,d]^{2}\bigr{]}&L_{2}^{m\times
n}\bigl{[}[a,b]^{2}\times[c,d]^{2}\bigr{]}\\\ L_{2}^{m\times
n}\bigl{[}[a,b]^{2}\times[c,d]\bigr{]}&L_{2}^{m\times
n}\bigl{[}[a,b]^{2}\times[c,d]^{2}\bigr{]}&L_{2}^{m\times
n}\bigl{[}[a,b]^{2}\times[c,d]^{2}\bigr{]}\end{array}\right]}=:\mathcal{N}_{2D}^{m\times
n}\bigl{[}[a,b]\\!\times\\![c,d]\bigr{]}$

we define the associated 9-PI operator
$\mathcal{P}[R]:=\mathcal{P}\left[\begin{smallmatrix}R_{00}&R_{01}&R_{02}\\\
R_{10}&R_{11}&R_{12}\\\
R_{20}&R_{21}&R_{22}\end{smallmatrix}\right]:L_{2}^{n}\bigl{[}[a_{1},b_{1}]\times[a_{2},b_{2}]\bigr{]}\rightarrow
L_{2}^{m}\bigl{[}[a_{1},b_{1}]\times[a_{2},b_{2}]\bigr{]}$ as

$\displaystyle\left(\mathcal{P}[R]\mathbf{x}\right)(s,r)=R_{00}(s,r)\mathbf{x}(s,r)$
$\displaystyle+\int_{c}^{r}\\!R_{01}(s,r,\nu)\mathbf{x}(s,\nu)d\nu+\int_{r}^{d}\\!R_{02}(s,r,\nu)\mathbf{x}(s,\nu)d\nu$
$\displaystyle+\int_{a}^{s}\\!R_{10}(s,r,\theta)\mathbf{x}(\theta,r)d\theta$
$\displaystyle+\int_{a}^{s}\\!\int_{c}^{r}\\!R_{11}(s,r,\theta,\nu)\mathbf{x}(\theta,\nu)d\nu
d\theta+\int_{a}^{s}\\!\int_{r}^{d}\\!R_{12}(s,r,\theta,\nu)\mathbf{x}(\theta,\nu)d\nu
d\theta$
$\displaystyle+\int_{s}^{b}R_{20}(s,r,\theta)\mathbf{x}(\theta,r)d\theta$
$\displaystyle+\int_{s}^{b}\\!\int_{c}^{r}\\!R_{21}(s,r,\theta,\nu)\mathbf{x}(\theta,\nu)d\nu
d\theta+\int_{s}^{b}\\!\int_{r}^{d}\\!R_{22}(s,r,\theta,\nu)\mathbf{x}(\theta,\nu)d\nu
d\theta$ (A.2)

for any $\mathbf{x}\in
L_{2}^{n}\bigl{[}[a_{1},b_{2}]\times[a_{2},b_{2}]\bigr{]}$.

Note that, similar to how 3-PI operators can be seen as a generalization of
matrices, operating on infinite-dimensional states $\mathbf{x}(s)$ instead of
a finite-dimensional vectors $x_{i}$, 9-PI operators are a generalization of
(a particular class of) tensors, operating on infinite-dimensional states
$\mathbf{x}(s,r)$ instead of matrix-valued states $x_{ij}$. However, this
comparison is not quite as easy to visualize as that between 3-PI operators
and matrices, so we will mostly use 3-PI operators to illustrate the different
properties of PI operators in the remaining sections.

Finally we define a general class of PI operators, encapsulating 3-PI
operators and 9-PI operators, as well as matrices and “cross-operators”. In
particular, we consider operators defined on the set
$Z^{\textnormal{n}}\bigr{[}[a,b],[c,d]\bigl{]}:={\scriptsize\left[\begin{array}[]{l}\mathbb{R}^{n_{0}}\\\
L_{2}^{n_{s}}[a,b]\\\ L_{2}^{n_{r}}[c,d]\\\
L_{2}^{n_{2}}\bigl{[}[a,b]\times[c,d]\bigr{]}\end{array}\right]}$, where
$\text{n}:=\\{n_{0},n_{s},n_{r},n_{2}\\}$, with each element being a coupled
state of finite-dimensional variables $x_{0}\in\mathbb{R}^{n_{0}}$, 1D
functions $\mathbf{x}_{s}\in L_{2}^{n_{s}}[a,b]$ and $\mathbf{x}_{r}\in
L_{2}^{n_{r}}[a,b]$, and 2D functions $\mathbf{x}_{2}\in
L_{2}^{n_{2}}\bigl{[}[a,b]\times[c,d]\bigr{]}$.

######  Definition 3 (PI Operator).

For any operator
$\mathcal{R}:Z^{\textnormal{n}}\bigr{[}[a,b],[c,d]\bigl{]}\rightarrow
Z^{\textnormal{m}}\bigr{[}[a,b],[c,d]\bigl{]}$ with
$\textnormal{m}:=\\{m_{0},m_{s},m_{r},m_{2}\\}$ and
$\textnormal{n}:=\\{n_{0},n_{s},n_{r},n_{2}\\}$, we say that $\mathcal{R}$ is
a PI operator, denoted by
$\mathcal{R}\in\Pi^{\textnormal{m}\times\textnormal{m}}$ if there exist
parameters

$\displaystyle
R:=\left[\\!\\!\begin{array}[]{llll}R_{00}&R_{0s}&R_{0r}&R_{02}\\\
R_{s0}&R_{ss}&R_{sr}&R_{s2}\\\ R_{r0}&R_{rs}&R_{rr}&R_{s2}\\\
R_{20}&R_{2s}&R_{2r}&R_{22}\end{array}\\!\\!\right]$
$\displaystyle\qquad\in\left[\\!\\!\begin{array}[]{llll}\mathbb{R}^{m_{0}\times
n_{0}}&L_{2}^{m_{0}\times n_{s}}[a,b]&L_{2}^{m_{0}\times
n_{r}}[c,d]&L_{2}^{m_{0}\times n_{2}}\bigl{[}[a,b]\times[c,d]\bigr{]}\\\
L_{2}^{m_{s}\times n_{0}}[a,b]&\mathcal{N}_{1D}^{m_{s}\times
n_{s}}[a,b]&L_{2}^{m_{s}\times
n_{r}}\bigl{[}[a,b]\times[c,d]\bigr{]}&\mathcal{N}_{1D\leftarrow
2D}^{m_{s}\times n_{2}}\bigl{[}[a,b],[c,d]\bigr{]}\\\ L_{2}^{m_{r}\times
n_{0}}[c,d]&L_{2}^{m_{r}\times
n_{s}}\bigl{[}[a,b]\times[c,d]\bigr{]}&\mathcal{N}_{1D}^{m_{r}\times
n_{r}}[c,d]&\mathcal{N}_{1D\leftarrow 2D}^{m_{r}\times
n_{2}}\bigl{[}[c,d],[a,b]\bigr{]}\\\ L_{2}^{m_{2}\times
n_{0}}\bigl{[}[a,b]\times[c,d]\bigr{]}&\mathcal{N}_{2D\leftarrow
1D}^{m_{2}\times n_{s}}\bigl{[}[a,b],[c,d]\bigr{]}&\mathcal{N}_{2D\leftarrow
1D}^{m_{2}\times
n_{r}}\bigl{[}[c,d],[a,b]\bigr{]}&\mathcal{N}_{2D}^{m_{2}\times
n_{r}}\bigl{[}[a,b]\times[c,d]\bigr{]}\end{array}\\!\\!\right]=:\mathcal{N}^{\textnormal{m}\times\textnormal{n}}\bigl{[}[a,b]\times[c,d]\bigr{]}$

such that

$\displaystyle\mathcal{R}$
$\displaystyle=\left(\mathcal{P}[R]\mathbf{x}\right)(s,r)$
$\displaystyle=\mathcal{P}\left[\begin{smallmatrix}R_{00}&R_{0s}&R_{0r}&R_{02}\\\
R_{s0}&R_{ss}&R_{sr}&R_{s2}\\\ R_{r0}&R_{rs}&R_{rr}&R_{s2}\\\
R_{20}&R_{2s}&R_{2r}&R_{22}\end{smallmatrix}\right]\left[\begin{smallmatrix}x_{0}\\\
\mathbf{x}_{s}\\\ \mathbf{x}_{r}\\\ \mathbf{x}_{2}\end{smallmatrix}\right]$
$\displaystyle:=\left[\\!\begin{array}[]{llll}R_{00}x_{0}&\\!\\!+\
\int_{a}^{b}\\!R_{0s}(s)\mathbf{x}_{s}(s)ds&\\!\\!+\
\int_{c}^{d}\\!R_{0r}(r)\mathbf{x}_{r}(r)dr&\\!\\!+\
\int_{a}^{b}\\!\int_{c}^{d}R_{02}\mathbf{x}_{2}(s,r)drds\\\
R_{s0}(s)x_{0}&\\!\\!+\
\bigl{(}\mathcal{P}[R_{ss}]\mathbf{x}_{s}\bigr{)}(s)&\\!\\!+\
\int_{c}^{d}\\!R_{sr}(s,r)\mathbf{x}_{r}(r)dr&\\!\\!+\
\bigl{(}\mathcal{P}[R_{s2}]\mathbf{x}_{2}\bigr{)}(s)\\\
R_{r0}(r)x_{0}&\\!\\!+\ \int_{a}^{b}\\!R_{rs}(s,r)\mathbf{x}_{s}(s)ds&\\!\\!+\
\bigl{(}\mathcal{P}[R_{rr}]\mathbf{x}_{r}\bigr{)}(r)&\\!\\!+\
\bigl{(}\mathcal{P}[R_{r2}]\mathbf{x}_{2}\bigr{)}(r)\\\
R_{20}(s,r)x_{0}&\\!\\!+\
\bigl{(}\mathcal{P}[R_{2s}]\mathbf{x}_{s}\bigr{)}(s,r)&\\!\\!+\
\bigl{(}\mathcal{P}[R_{2r}]\mathbf{x}_{r}\bigr{)}(s,r)&\\!\\!+\
\bigl{(}\mathcal{P}[R_{22}]\mathbf{x}_{2}\bigr{)}(s,r)\end{array}\\!\right],$
(A.7)

for any $\mathbf{x}=\left[\begin{smallmatrix}x_{0}\\\ \mathbf{x}_{s}\\\
\mathbf{x}_{r}\\\
\mathbf{x}_{2}\end{smallmatrix}\right]\in{\scriptsize\left[\begin{array}[]{l}\mathbb{R}^{n_{0}}\\\
L_{2}^{n_{s}}[a,b]\\\ L_{2}^{n_{r}}[c,d]\\\
L_{2}^{n_{2}}\bigl{[}[a,b]\times[c,d]\bigr{]}\end{array}\right]}=:Z^{\textnormal{n}}$,
where for given parameters

$\displaystyle P$ $\displaystyle:=\\{P_{0},P_{1},P_{2}\\}$
$\displaystyle\qquad\in\left\\{L_{2}^{m_{s}\times
n_{2}}\bigl{[}[a,b]\\!\times\\![c,d]\bigr{]},L_{2}^{m_{s}\times
n_{2}}\bigl{[}[a,b]^{2}\\!\times\\![c,d]\bigr{]},L_{2}^{m_{s}\times
n_{2}}\bigl{[}[a,b]^{2}\\!\times\\![c,d]\bigr{]}\right\\}=:\mathcal{N}^{m_{s}\times
n_{2}}_{1D\leftarrow 2D}\bigl{[}[a,b],[c,d]\bigr{]},$ $\displaystyle Q$
$\displaystyle:=\\{Q_{0},Q_{1},Q_{2}\\}$
$\displaystyle\qquad\in\left\\{L_{2}^{m_{s}\times
n_{2}}\bigl{[}[a,b]\\!\times\\![c,d]\bigr{]},L_{2}^{m_{s}\times
n_{2}}\bigl{[}[a,b]^{2}\\!\times\\![c,d]\bigr{]},L_{2}^{m_{s}\times
n_{2}}\bigl{[}[a,b]^{2}\\!\times\\![c,d]\bigr{]}\right\\}=:\mathcal{N}^{m_{s}\times
n_{2}}_{2D\leftarrow 1D}\bigl{[}[a,b],[c,d]\bigr{]},$

we define

$\displaystyle\bigl{(}\mathcal{P}[P]\mathbf{x}_{2}\bigr{)}(s)$
$\displaystyle:=\int_{c}^{d}\biggl{[}P_{0}(s,r)\mathbf{x}_{2}(s,r)+\int_{a}^{s}P_{1}(s,r,\theta)\mathbf{x}_{2}(\theta,r)d\theta+\int_{s}^{b}P_{2}(s,r,\theta)\mathbf{x}(\theta,r)d\theta\biggr{]}dr,$
$\displaystyle\bigl{(}\mathcal{P}[Q]\mathbf{x}_{s}\bigr{)}(s,r)$
$\displaystyle:=Q_{0}(s,r)\mathbf{x}_{s}(s)+\int_{a}^{s}Q_{1}(s,r,\theta)\mathbf{x}_{s}(\theta)d\theta+\int_{s}^{b}Q_{2}(s,r,\theta)\mathbf{x}_{s}(\theta)d\theta,$

for any $\mathbf{x}_{2}\in L_{2}^{n_{2}}\bigl{[}[a,b]\times[c,d]\bigr{]}$ and
$\mathbf{x}_{s}\in L_{2}^{n_{s}}[a,b]$.

####  A.2 Addition of PI Operators

An obvious but crucial property of PI operators is that the sum of two PI
operators (of appropriate dimensions) is again a PI operator.

######  Lemma 4.

For any PI parameters
$Q,R\in\mathcal{N}^{\text{m}\times\text{n}}\bigl{[}[a,b],[c,d]\bigr{]}$, there
exist unique parameters
$P\in\mathcal{N}^{\text{m}\times\text{n}}\bigl{[}[a,b],[c,d]\bigr{]}$ such
that

$\displaystyle\mathcal{P}[R]+\mathcal{P}[Q]=\mathcal{P}[P].$

That is, for any $\mathbf{x}\in
Z^{\textnormal{n}}\bigl{[}[a,b],[c,d]\bigr{]}$,

$\displaystyle\bigl{(}(\mathcal{P}[Q]+\mathcal{P}[R])\mathbf{x}\bigr{)}(s)=\bigl{(}\mathcal{P}[P]\mathbf{x})(s).$

###### Proof.

We outline the proof for 3-PI operators, for which it is easy to see that, by
linearity of the integral,

$\displaystyle\bigl{(}\mathcal{P}_{\\{R_{0},R_{1},R_{2}\\}}\mathbf{x}\bigr{)}(s)+\bigl{(}\mathcal{P}_{\\{Q_{0},Q_{1},Q_{2}\\}}\mathbf{x}\bigr{)}(s)$
$\displaystyle\qquad=[R_{0}(s)+Q_{0}(s)]\mathbf{x}(s)+\int_{a}^{s}[R_{1}(s,\theta)+Q_{1}(s,\theta)]\mathbf{x}(\theta)d\theta+\int_{s}^{b}[R_{2}(s,\theta)+Q_{2}(s,\theta)]\mathbf{x}(\theta)d\theta$
$\displaystyle\qquad\qquad=\bigl{(}\mathcal{P}_{\\{R_{0}+Q_{0},R_{1}+Q_{1},R_{2}+Q_{2}\\}}\mathbf{x}\bigr{)}(s).$

Similar results can be easily derived for more general PI operators. For a
full proof, we refer to [5]. ∎

Comparing the addition operation for $3$-PI operators to that for matrices
$A,B\in\mathbb{R}^{m\times n}$, we can draw direct parallels. In particular,
where the sum $C=A+B$ of two matrices is simply computed by adding the
elements $[C]_{ij}=[A]_{ij}+[B]_{ij}$ for each row $i$ and column $j$, the sum
of two 3-PI operators is computed by simply adding the values of the
parameters $P(s,\theta)=Q(s,\theta)+R(s,\theta)$ at each position $s$ and
$\theta$ within the domain.

####  A.3 Composition of PI Operators

In addition to the sum of two PI operators being a PI operator, the
composition of two PI operators can also be shown to be a PI operator, as
stated in the following lemma:

######  Lemma 5.

For any PI parameters
$R_{1}\in\mathcal{N}^{\text{m}\times\text{p}}\bigl{[}[a,b],[c,d]\bigr{]}$ and
$R_{2}\in\mathcal{N}^{\text{p}\times\text{n}}\bigl{[}[a,b],[c,d]\bigr{]}$,
there exist unique parameters
$R_{3}\in\mathcal{N}^{\text{m}\times\text{n}}\bigl{[}[a,b],[c,d]\bigr{]}$ such
that

$\displaystyle\mathcal{P}[R_{1}]\circ\mathcal{P}[R_{2}]=\mathcal{P}[R_{3}].$

That is, for any $\mathbf{x}\in
Z^{\textnormal{n}}\bigl{[}[a,b],[c,d]\bigr{]}$,

$\displaystyle\Bigl{(}\mathcal{P}[R_{1}]\bigl{(}\mathcal{P}[R_{2}]\mathbf{x}\bigr{)}\Bigr{)}(s)=\bigl{(}\mathcal{P}[R_{3}]\mathbf{x})(s).$

###### Proof.

We once again outline the proof only for 3-PI operators. For this, we define
the indicator function

$\displaystyle\mathbf{I}(s-\theta)=\begin{cases}1,&\text{if }s\geq\theta\\\
0.&\text{else}\end{cases}$

allowing us to write, e.g.

$\displaystyle\bigl{(}\mathcal{P}_{\\{R_{0},R_{1},R_{2}\\}}\mathbf{x}\bigl{)}(s)=R_{0}(s)\mathbf{x}+\int_{a}^{b}\Bigl{[}\mathbf{I}(s-\theta)R_{1}(s,\theta)+\mathbf{I}(\theta-s)R_{2}(s,\theta)\Bigr{]}\mathbf{x}(\theta)d\theta.$

Furthermore, we have the following relations for the indicator function

$\displaystyle\mathbf{I}(s-\eta)\mathbf{I}(\eta-\theta)$
$\displaystyle=\begin{cases}\mathbf{I}(s-\theta),&\text{if
}\eta\in[\theta,s],\\\ 0,&\text{else,}\end{cases}$
$\displaystyle\mathbf{I}(s-\eta)\mathbf{I}(\theta-\eta)$
$\displaystyle=\mathbf{I}(s-\theta)\mathbf{I}(\theta-\eta)+\mathbf{I}(\theta-s)\mathbf{I}(s-\eta)$

Using the first relation, it follows that for any $R_{1},Q_{1}\in
L_{2}^{m_{s}\times n_{s}}\bigl{[}[a,b]^{2}\bigr{]}$,

$\displaystyle\Bigl{(}\mathcal{P}_{\\{0,R_{1},0\\}}\bigl{(}\mathcal{P}_{\\{0,Q_{1},0\\}}\mathbf{x}\bigr{)}\Bigr{)}(s)$
$\displaystyle=\int_{a}^{s}R_{1}(s,\eta)\int_{a}^{\eta}Q_{1}(\eta,\theta)\mathbf{x}(\theta)d\theta
d\eta$
$\displaystyle=\int_{a}^{b}\int_{a}^{b}\mathbf{I}(s-\eta)\mathbf{I}(\eta-\theta)R_{1}(s,\eta)Q_{1}(\eta,\theta)\mathbf{x}(\theta)d\theta
d\eta$
$\displaystyle=\int_{a}^{s}\left[\int_{\theta}^{s}R_{1}(s,\eta)Q_{1}(\eta,\theta)d\eta\right]\mathbf{x}(\theta)d\theta=\bigl{(}\mathcal{P}_{\\{0,P_{11},0\\}}\mathbf{x}\bigr{)}(s),$

where
$P_{11}(s,\theta):=\int_{\theta}^{s}R_{1}(s,\eta)Q_{1}(\theta,\eta)d\eta$.
Similarly, we can show that

$\displaystyle\Bigl{(}\mathcal{P}_{\\{0,R_{1},0\\}}\bigl{(}\mathcal{P}_{\\{0,0,Q_{2}\\}}\mathbf{x}\bigr{)}\Bigr{)}(s)$
$\displaystyle=\int_{a}^{s}\left[\int_{a}^{\theta}R_{1}(s,\eta)Q_{2}(\eta,\theta)d\eta\right]\mathbf{x}(\theta)d\theta+\int_{s}^{b}\left[\int_{a}^{s}R_{1}(s,\eta)Q_{2}(\eta,\theta)d\eta\right]\mathbf{x}(\theta)d\theta$
$\displaystyle\Bigl{(}\mathcal{P}_{\\{0,0,R_{2}\\}}\bigl{(}\mathcal{P}_{\\{0,Q_{1},0\\}}\mathbf{x}\bigr{)}\Bigr{)}(s)$
$\displaystyle=\int_{a}^{s}\left[\int_{s}^{b}R_{2}(s,\eta)Q_{1}(\eta,\theta)d\eta\right]\mathbf{x}(\theta)d\theta+\int_{s}^{b}\left[\int_{\theta}^{b}R_{2}(s,\eta)Q_{1}(\eta,\theta)d\eta\right]\mathbf{x}(\theta)d\theta$
$\displaystyle\Bigl{(}\mathcal{P}_{\\{0,0,R_{2}\\}}\bigl{(}\mathcal{P}_{\\{0,0,Q_{2}\\}}\mathbf{x}\bigr{)}\Bigr{)}(s)$
$\displaystyle=\int_{s}^{b}\left[\int_{s}^{\theta}R_{2}(s,\eta)Q_{2}(\theta,\eta)d\eta\right]\mathbf{x}(\theta)d\theta=\bigl{(}\mathcal{P}_{\\{0,0,P_{22}\\}}\mathbf{x}\bigr{)}(s),$

proving that the composition of lower-triangular and upper-triangular partial
integrals can always be expressed as partial integrals as well. It is also
easy easy to see that

$\displaystyle\Bigl{(}\mathcal{P}_{\\{R_{0},0,0\\}}\bigl{(}\mathcal{P}_{\\{0,Q_{1},Q_{2}\\}}\mathbf{x}\bigr{)}\Bigr{)}(s)$
$\displaystyle=\int_{a}^{s}R_{0}(s)Q_{1}(s,\theta)\mathbf{x}(\theta)d\theta+\int_{s}^{b}R_{0}(s)Q_{2}(s,\theta)\mathbf{x}(\theta)d\theta=\bigl{(}\mathcal{P}_{\\{0,P_{01},P_{02}\\}}\mathbf{x}\bigr{)}(s),$
$\displaystyle\Bigl{(}\mathcal{P}_{\\{0,R_{1},R_{2}\\}}\bigl{(}\mathcal{P}_{\\{Q_{0},0,0\\}}\mathbf{x}\bigr{)}\Bigr{)}(s)$
$\displaystyle=\int_{a}^{s}R_{1}(s,\theta)Q_{0}(\theta)\mathbf{x}(\theta)d\theta+\int_{s}^{b}R_{2}(s,\theta)Q_{2}(\theta)\mathbf{x}(\theta)d\theta=\bigl{(}\mathcal{P}_{\\{0,P_{10},P_{20}\\}}\mathbf{x}\bigr{)}(s),$

from which it follows that the composition of any 3-PI operators can be
expressed as a 3-PI operator as well. Moreover, since we can repeat these
steps along any spatial directions, this result extends to more general (2D)
PI operators as well. For a full proof, we refer to .

∎

We note again the similarity to matrices: just like the product of two lower
triangular matrices $L_{1},L_{2}$ is a lower triangular matrix $L_{3}$, the
composition of two lower-triangular 3-PI operators
$\mathcal{P}_{\\{0,R_{1},0\\}},\mathcal{P}_{\\{0,Q_{1},0\\}}$ is also a lower-
triangular 3-PI operator $\mathcal{P}_{\\{0,P_{11},0\\}}$. Similarly, the
product of two upper-triangular 3-PI operators
$\mathcal{P}_{\\{0,0,R_{2}\\}},\mathcal{P}_{\\{0,0,Q_{2}\\}}$ is also an
upper-triangular 3-PI operator $\mathcal{P}_{\\{0,0,P_{22}\\}}$, but the
composition of lower- and upper-triangular PI operators need not be lower- or
upper-triangular – just as with matrices. Finally, the composition of a
diagonal operator $\mathcal{P}_{\\{R_{0},0,0\\}}$ with a lower- or upper-
diagonal PI operator is also respectively lower- or upper-diagonal.

####  A.4 Adjoint of PI Operators

To define the adjoint of a PI operator
$\mathcal{R}\in\Pi^{\textnormal{m}\times\text{normal{n}}}$, we first recall
the definition of the function space that these operators map:
$Z^{\textnormal{n}}\bigr{[}[a,b],[c,d]\bigl{]}:={\scriptsize\left[\begin{array}[]{l}\mathbb{R}^{n_{0}}\\\
L_{2}^{n_{s}}[a,b]\\\ L_{2}^{n_{r}}[c,d]\\\
L_{2}^{n_{2}}\bigl{[}[a,b]\times[c,d]\bigr{]}\end{array}\right]}$, where
$\text{n}:=\\{n_{0},n_{s},n_{r},n_{2}\\}$, with each element being a coupled
state of finite-dimensional variables $x_{0}\in\mathbb{R}^{n_{0}}$, 1D
functions $\mathbf{x}_{s}\in L_{2}^{n_{s}}[a,b]$ and $\mathbf{x}_{r}\in
L_{2}^{n_{r}}[a,b]$, and 2D functions $\mathbf{x}_{2}\in
L_{2}^{n_{2}}\bigl{[}[a,b]\times[c,d]\bigr{]}$. We endow this space with the
inner product

$\displaystyle\left\langle{\mathbf{x}},{\mathbf{y}}\right\rangle_{Z}$
$\displaystyle=\left\langle{x_{0}},{y_{0}}\right\rangle+\left\langle{\mathbf{x}_{s}},{\mathbf{y}_{s}}\right\rangle_{L_{2}}+\left\langle{\mathbf{x}_{r}},{\mathbf{y}_{r}}\right\rangle_{L_{2}}+\left\langle{\mathbf{x}_{2}},{\mathbf{y}_{2}}\right\rangle_{L_{2}}$
$\displaystyle=x_{0}^{T}y_{0}+\int_{a}^{b}[\mathbf{x}_{s}(s)]^{T}\mathbf{y}(s)ds+\int_{c}^{d}[\mathbf{x}_{r}(r)]^{T}\mathbf{y}(r)dr+\int_{a}^{b}\int_{c}^{d}[\mathbf{x}_{2}(s,r)]^{T}\mathbf{y}(s,r)drds$

Defining this inner product, we can also define the adjoint of PI operators.

######  Lemma 6.

For any PI parameters
$R\in\mathcal{N}^{\text{m}\times\text{n}}\bigl{[}[a,b],[c,d]\bigr{]}$, there
exist unique parameters
$Q\in\mathcal{N}^{\text{n}\times\text{m}}\bigl{[}[a,b],[c,d]\bigr{]}$ such
that

$\displaystyle\bigl{(}\mathcal{P}[R]\bigr{)}^{*}=\mathcal{P}[Q],$

where $\mathcal{P}^{*}$ denotes the adjoint of a PI operator $\mathcal{P}$.
That is, for any $\mathbf{x}\in Z^{\textnormal{n}}\bigl{[}[a,b],[c,d]\bigr{]}$
and $\mathbf{y}\in Z^{\textnormal{m}}\bigl{[}[a,b],[c,d]\bigr{]}$,

$\displaystyle\left\langle{\mathcal{P}[R]\mathbf{x}},{\mathbf{y}}\right\rangle_{Z}=\left\langle{\mathbf{x}},{\mathcal{P}[Q]\mathbf{y}}\right\rangle_{Z}.$

###### Proof.

We outline the proof only for 4-PI operators. In particular, let
$\textnormal{n}=\\{n_{0},n_{1},0,0\\}$ and
$\textnormal{m}=\\{m_{0},m_{1},0,0\\}$, and let
$B=\left[\begin{smallmatrix}P&Q_{1}\\\
Q_{2}&\\{R_{0},R_{1},R_{2}\\}\end{smallmatrix}\right]$ define a 4-PI operator
$\mathcal{P}[B]\in\Pi^{\textnormal{n}\times\textnormal{m}}$. Define
$\hat{B}=\left[\begin{smallmatrix}\hat{P}&\hat{Q}_{1}\\\
\hat{Q}_{2}&\\{\hat{R}_{0},\hat{R}_{1},\hat{R}_{2}\\}\end{smallmatrix}\right]$,
where

$\displaystyle\hat{P}$ $\displaystyle=P^{T},$ $\displaystyle\hat{Q}_{1}(s)$
$\displaystyle=Q_{2}^{T}(s),$ $\displaystyle\hat{R}_{1}(s,\theta)$
$\displaystyle=R_{2}^{T}(\theta,s),$ $\displaystyle\hat{Q}_{2}(s)$
$\displaystyle=Q_{1}^{T}(s),$ $\displaystyle\hat{R}_{0}(s)$
$\displaystyle=R_{0}^{T}(s),$ $\displaystyle\hat{R}_{2}(s,\theta)$
$\displaystyle=R_{1}^{T}(\theta,s),$

Then, for arbitrary $\mathbf{x}=\left[\begin{smallmatrix}x_{0}\\\
\mathbf{x}_{1}\end{smallmatrix}\right]\in Z^{\textnormal{n}}$ and
$\mathbf{y}=\left[\begin{smallmatrix}y_{0}\\\
\mathbf{y}_{1}\end{smallmatrix}\right]\in Z^{\textnormal{m}}$, we note that

$\displaystyle\left\langle{\mathcal{P}[B]\mathbf{x}},{\mathbf{y}}\right\rangle_{Z}$
$\displaystyle=[Px_{0}]^{T}y_{0}+\left[\int_{a}^{b}Q_{1}(s)\mathbf{x}_{1}(s)ds\right]^{T}y_{0}+\int_{a}^{s}[Q_{2}(s)x_{0}]^{T}\mathbf{y}_{1}(s)ds+\int_{a}^{b}[R_{0}(s)\mathbf{x}_{1}(s)]^{T}\mathbf{y}_{1}(s)ds$
$\displaystyle\qquad+\int_{a}^{b}\left[\left(\int_{a}^{b}\mathbf{I}(s-\theta)R_{1}(s,\theta)+\mathbf{I}(\theta-s)R_{2}(s,\theta)\right)\mathbf{x}_{1}(\theta)d\theta\right]^{T}\mathbf{y}_{1}(s)ds$
$\displaystyle=x_{0}[P^{T}y_{0}]+x_{0}^{T}\left[\int_{a}^{b}Q_{2}^{T}(s)\mathbf{y}_{1}(s)ds\right]+\int_{a}^{b}\mathbf{x}_{1}^{T}(s)\left[Q_{1}^{T}(s)y_{0}\right]ds+\int_{a}^{b}\mathbf{x}_{1}^{T}(s)\left[R_{0}^{T}(s)\mathbf{y}_{1}(s)\right]ds$
$\displaystyle\qquad+\int_{a}^{b}\mathbf{x}_{1}^{T}(s)\left[\left(\int_{a}^{b}\mathbf{I}(-\theta-s)R_{1}(\theta,s)+\mathbf{I}(s-\theta)R_{2}(\theta,s)\right)\mathbf{y}_{1}(\theta)d\theta\right]ds$
$\displaystyle=\left\langle{\mathbf{x}},{\mathcal{P}[\hat{B}]\mathbf{y}}\right\rangle_{Z}.$

∎

We note again the similarities with matrices: Just like the adjoint of a
matrix can be determined by switching the rows and columns, the adjoint of a
3-PI operator is determined by switching the primary and dual variables
$(s,\theta)$, as well as switching the lower- and upper-triangular parts.

####  A.5 Inversion of PI operators

In this section, we focus on constructing the controller from the feasible
solutions $\mathcal{P}$ and $\mathcal{Z}$ to the LPIs described in previous
sections. The controller gains $\mathcal{K}$ is given by the relation
$\mathcal{K}=\mathcal{Z}\mathcal{P}^{-1}$. However, $\mathcal{P}^{-1}$,
although is a 4-PI operator, may not have polynomial parameters. Hence the
inverse is calculated numerically. First, we propose a numerical method to
find the inverse of $\mathcal{P}$. Then, we use the numerical inverse to
construct the controller gains. To find the inverse of a 4-PI operator,
$\mathcal{P}{\footnotesize\begin{bmatrix}P,&Q\\\
Q^{T},&\left\\{R_{i}\right\\}\end{bmatrix}}$, we first find the inverse of
3-PI operators of the form $\mathcal{P}_{\\{I,H_{1},R_{2}\\}}$ and then
express the inverse of $\mathcal{P}{\footnotesize\begin{bmatrix}P,&Q\\\
Q^{T},&\left\\{R_{i}\right\\}\end{bmatrix}}$ in terms of
$\mathcal{P}_{\\{I,H_{1},H_{2}\\}}^{-1}$ where $H_{i}$ are dependent on the
parameters $\\{P,Q,R_{i}\\}$.

#####  A.5.1 Inversion of 3-PI operators

First, we note that any matrix-valued polynomial $H(s,\theta)$ can be factored
as $F(s)G(\theta)$. Then, for any given 3-PI operator
$\mathcal{P}_{\\{I,H_{1},H_{2}\\}}$ with matrix-valued polynomial parameters
$H_{1}$ and $H_{2}$, we have

$\displaystyle\mathcal{P}_{\\{I,H_{1},H_{2}\\}}$
$\displaystyle=\mathcal{P}_{\\{I,-F_{1}G_{1},-F_{2}G_{2}\\}},\leavevmode\nobreak\
\text{where}$ $\displaystyle H_{i}(s,\theta)$
$\displaystyle=-F_{i}(s)G_{i}(\theta),$

for some matrix-valued polynomials $F_{i}$ and $G_{i}$. We can now find an
inverse for $\mathcal{P}_{\\{I,H_{1},H_{2}\\}}$ using the following result.

######  Lemma 7.

Suppose $F_{1}:[a,b]\to\mathbb{R}^{n\times p}$,
$G_{1}:[a,b]\to\mathbb{R}^{p\times n}$, $F_{2}:[a,b]\to\mathbb{R}^{n\times
q}$, $G_{2}:[a,b]\to\mathbb{R}^{q\times n}$ and $U$ is the unique function
that satisfies the equation

$\displaystyle U(s)$
$\displaystyle=I_{(p+q)}+\int_{a}^{s}U(t)\begin{bmatrix}G_{1}(t)F_{1}(t)&G_{1}(t)F_{2}(t)\\\
-G_{2}(t)F_{1}(t)&-G_{2}(t)F_{2}(t)\end{bmatrix}dt,$

where $U$ is partitioned as

$U=\begin{bmatrix}U_{11}&U_{12}\\\ U_{21}&U_{22}\end{bmatrix},\quad
U_{11}(s)\in\mathbb{R}^{p\times p},U_{22}(s)\in\mathbb{R}^{q\times q}.$

Then, the 3-PI operator $\mathcal{P}_{\\{I,-F_{1}G_{1},-F_{2}G_{2}\\}}$ is
invertible if and only if $U_{22}(b)$ is invertible and

$\displaystyle(\mathcal{P}_{\\{I,-F_{1}G_{1},-F_{2}G_{2}\\}})^{-1}$
$\displaystyle=\mathcal{P}_{\\{I,L_{1},L_{2}\\}},$

where

$\displaystyle L_{1}(s,t)$
$\displaystyle=\begin{bmatrix}F_{1}(s)&F_{2}(s)\end{bmatrix}U(s)V(t)\begin{bmatrix}G_{1}(t)\\\
-G_{2}(t)\end{bmatrix}-L_{2}(s,t),$ (A.8a) $\displaystyle L_{2}(s,t)$
$\displaystyle=-\begin{bmatrix}F_{1}(s)&F_{2}(s)\end{bmatrix}U(s)PV(t)\begin{bmatrix}G_{1}(t)\\\
-G_{2}(t)\end{bmatrix},$ (A.8b) $\displaystyle P$
$\displaystyle=\begin{bmatrix}0_{p\times p}&0_{p\times q}\\\
U_{22}(b)^{-1}U_{21}(b)&I_{q}\end{bmatrix},$

and $V$ is the unique function satisfying the equation

$V(t)=I_{(p+q)}-\int_{a}^{t}\begin{bmatrix}G_{1}(s)F_{1}(s)&G_{1}(s)F_{2}(s)\\\
-G_{2}(s)F_{1}(s)&-G_{2}(s)F_{2}(s)\end{bmatrix}V(s)ds.$

###### Proof.

Proof can be found in [3, Chapter IX.2]. ∎

Using Lemma 2.2. of [3, Chapter IX.2], we can use an iterative process and
numerical integration to approximate $U$ and $V$ functions in the above result
to find an inverse for the 3-PI operators of the form
$\mathcal{P}_{\\{I,R_{1},R_{2}\\}}$ where $R_{i}$ are matrix-valued
polynomials. By extension, given an $R_{0}$ invertible, we can obtain the
inverse of a general 3-PI operator as shown below.

######  Corollary 8.

Suppose $R_{0}:[a,b]\to\mathbb{R}^{n\times n}$,
$R_{1},R_{2}:[a,b]^{2}\to\mathbb{R}^{n\times n}$, with $R_{0}$ invertible on
$[a,b]$. Then, the inverse of the 3-PI operator, $\mathcal{P}_{\\{R_{i}\\}}$,
is given by $\mathcal{P}_{\\{\hat{R}_{0},\hat{R}_{1},\hat{R}_{2}\\}}$ where

$\displaystyle\hat{R}_{0}(s)=R_{0}(s)^{-1},\quad\hat{R}_{i}(s,\theta)=L_{i}(s,\theta)\hat{R}_{0}(\theta),\quad
i\in\\{1,2\\},$

where $L_{1}$ and $L_{2}$ are as defined in (A.8) for functions $F_{i}$ and
$G_{i}$ such that $F_{i}(s)G_{i}(\theta)=R_{0}(s)^{-1}R_{i}(s,\theta)$.

###### Proof.

Let $R_{i}$ be as stated above.

$\displaystyle(\mathcal{P}_{\\{R_{0}(s),R_{1}(s,\theta),R_{2}(s,\theta)\\}})^{-1}$
$\displaystyle=(\mathcal{P}_{\\{R_{0}(s),0,0\\}}\mathcal{P}_{\\{I,R_{0}(s)^{-1}R_{1}(s,\theta),R_{0}(s)^{-1}R_{2}(s,\theta)\\}})^{-1}$
$\displaystyle=(\mathcal{P}_{\\{I,R_{0}(s)^{-1}R_{1}(s,\theta),R_{0}(s)^{-1}R_{2}(s,\theta)\\}})^{-1}(\mathcal{P}_{\\{R_{0}(s),0,0\\}})^{-1}$
$\displaystyle=\mathcal{P}_{\\{I,L_{1}(s,\theta),L_{2}(s,\theta)\\}}\mathcal{P}_{\\{R_{0}(s)^{-1},0,0\\}}$
$\displaystyle=\mathcal{P}_{\\{R_{0}(s)^{-1},L_{1}(s,\theta)R_{0}(\theta)^{-1},L_{2}(s,\theta)R_{0}(\theta)^{-1}\\}}.$

where $L_{i}$ are obtained from the Lemma 7 and the composition of PI
operators is performed using the formulae in [8].

∎

Note the above expressions for the inverse are exact, however, in practice,
$R_{0}^{-1}$ may not have an analytical expression (or very hard to
determine). Thus, finding $F_{i}$ and $G_{i}$ such that
$F_{i}(s)G_{i}(\theta)=R_{0}(s)^{-1}R_{i}(s,\theta)$ may not be possible. To
overcome this problem, we approximate $R_{0}^{-1}$ by a polynomial which
guarantees that $R_{0}^{-1}R_{i}$ are polynomials and can be factorized into
$F_{i}$ and $G_{i}$. Using this approach, we can find an approximate inverse
for $\mathcal{P}_{\\{R_{i}\\}}$ using Lemma 7.

#####  A.5.2 Inversion of 4-PI operators

Given, $R_{0},R_{1},R_{2}$ with $R_{0}$ invertible, we proposed a way to find
the inverse of the operator $\mathcal{P}_{\\{R_{i}\\}}$. Now, we use this
method to find the inverse of a 4-PI operator
$\mathcal{P}{\footnotesize\begin{bmatrix}P,&Q\\\
Q^{T},&\left\\{R_{i}\right\\}\end{bmatrix}}$. Given $P$, $Q$ and $R_{i}$ with
invertible $P$ and $R_{0}$, define the 3-PI operator
$\mathcal{P}_{\\{H_{i}\\}}$ with parameters $H_{i}$

$\displaystyle H_{0}(s)=R_{0}(s),\quad H_{i}(s,\theta)$
$\displaystyle=R_{i}(s,\theta)-Q(s)^{T}P^{-1}Q(\theta),$
$\displaystyle\qquad\qquad\quad i\in\\{1,2\\}.$

Next we suppose that $\mathcal{P}_{\\{H_{i}\\}}$ is invertible. We will use
the inverse of this 3-PI operator to find the inverse of the 4-PI operator as
follows.

######  Corollary 9.

Suppose $P\in\mathbb{R}^{m\times m}$, $Q:[a,b]\to\mathbb{R}^{m\times n}$,
$R_{0}:[a,b]\to\mathbb{R}^{n\times n}$ and
$R_{1},R_{2}:[a,b]^{2}\to\mathbb{R}^{n\times n}$ are matrices and matrix-
valued polynomials with $P$ and $\mathcal{P}_{\\{H_{i}\\}}$ invertible such
that $(\mathcal{P}_{\\{H_{i}\\}})^{-1}=\mathcal{P}_{\\{L_{i}\\}}$, where

$\displaystyle H_{0}(s)$ $\displaystyle=R_{0}(s),\quad
H_{i}(s,\theta)=R_{i}(s,\theta)-Q(s)^{T}P^{-1}Q(\theta),$

for $\leavevmode\nobreak\ i\in\\{1,2\\}$. Then, the inverse of
$\mathcal{P}{\footnotesize\begin{bmatrix}P,&Q\\\
Q^{T},&\left\\{R_{i}\right\\}\end{bmatrix}}$ is given by

$\displaystyle\left(\mathcal{P}{\footnotesize\begin{bmatrix}P,&Q\\\
Q^{T},&\left\\{R_{i}\right\\}\end{bmatrix}}\right)^{-1}$
$\displaystyle=\mathcal{P}{\footnotesize\begin{bmatrix}\hat{P},&\hat{Q}\\\
\hat{Q}^{T},&\left\\{\hat{R}_{i}\right\\}\end{bmatrix}},$

where $\hat{P},\hat{Q},\hat{R}_{i}$ are defined as

$\displaystyle\hat{R}_{0}(s)$
$\displaystyle=R_{0}(s)^{-1},\quad\hat{R}_{i}(s,\theta)=L_{i}(s,\theta),\quad\text{for}\leavevmode\nobreak\
i\in\\{1,2\\},$ $\displaystyle\hat{P}$
$\displaystyle=P^{-1}+P^{-1}\left(\int_{a}^{b}(Q(s)\hat{R}_{0}(s)Q(s)^{T}\right.$
$\displaystyle\qquad\left.+\int_{s}^{b}Q(\theta)L_{1}(\theta,s)Q(\theta)^{T}d\theta\right.$
$\displaystyle\quad\qquad\left.+\int_{a}^{s}Q(\theta)L_{2}(\theta,s)Q(\theta)^{T}d\theta\right)dsP^{-1}$
$\displaystyle\hat{Q}(s)$
$\displaystyle=-P^{-1}\left(Q(s)\hat{R}_{0}(s)+\int_{s}^{b}Q(\theta)L_{1}(\theta,s)d\theta\right.$
$\displaystyle\qquad\left.+\int_{a}^{s}Q(\theta)L_{2}(\theta,s)d\theta\right).$
(A.9)

###### Proof.

This result can be proved by performing the composition of $\mathcal{P}$ and
its inverse $\mathcal{P}^{-1}$ where the inverse is provided by the formulae
stated above. Using composition formulae for 4-PI operators (See [8]), we show
that the resulting operator is an identity. ∎

####  A.6 Composition of Differential and PI operator

Given the well-known relationship between integrals and derivatives (think
e.g. the fundamental theorem of calculus), it is natural to assume that the
composition of a differential operator and a PI operator may be expressed as a
PI operator as well. Unfortunately, this is not true in general, as e.g. the
operator $\mathcal{P}$ defined as
$\bigl{(}\mathcal{P}\mathbf{v}\bigr{)}(s)=P(s)\mathbf{v}(s)$ is a PI operator,
but there clearly does not exist a PI operator $\mathcal{Q}$ such that
$\partial_{s}\bigl{(}\mathcal{P}\mathbf{v}\bigr{)}(s)=\bigl{(}\mathcal{Q}\mathbf{v}\bigr{)}(s)$.
Nevertheless, if the function $\mathbf{v}$ is differentiable, i.e.
$\mathbf{v}\in H_{1}$, then we can always express the derivative
$\bigl{(}\mathcal{P}\mathbf{v}\bigr{)}(s)$ for a PI operator $\mathcal{P}$ in
terms of $\mathbf{v}(s)$ and $\partial_{s}\mathbf{v}(s)$ as
$\partial_{s}\bigl{(}\mathcal{P}\mathbf{v}\bigr{)}(s)=\bigl{(}\mathcal{Q}\left[\begin{smallmatrix}\mathbf{v}\\\
\partial_{s}\mathbf{v}\end{smallmatrix}\right]\bigr{)}(s)$, as we show in the
next lemma.

######  Lemma 10.

Suppose $\mathcal{P}{\footnotesize\begin{bmatrix}P,&Q_{1}\\\
Q_{2},&\left\\{R_{i}\right\\}\end{bmatrix}}:\mathbb{R}^{m}\times
H_{1}^{n}\to\mathbb{R}^{p}\times L_{2}^{q}$, and define
$\partial_{s}\mathcal{P}{\footnotesize\begin{bmatrix}P,&Q_{1}\\\
Q_{2},&\left\\{R_{i}\right\\}\end{bmatrix}}:\mathbb{R}^{m}\times
H_{1}^{n}\times L_{2}^{n}\to\mathbb{R}^{p}\times L_{2}^{q}$ as

$\displaystyle\partial_{s}\mathcal{P}{\footnotesize\begin{bmatrix}P,&Q_{1}\\\
Q_{2},&\left\\{R_{i}\right\\}\end{bmatrix}}=\mathcal{P}{\footnotesize\begin{bmatrix}0,&0\\\
\bar{Q}_{2},&\left\\{\bar{R}_{i}\right\\}\end{bmatrix}}$ (A.10)

where $\bar{Q}_{2}(s)=\partial_{s}Q_{2}(s)$,
$\bar{R}_{0}(s)=\begin{bmatrix}\partial_{s}R_{0}(s)+R_{1}(s,s)-R_{2}(s,s)&R_{0}(s)\end{bmatrix}$
and
$\bar{R}_{i}(s,\theta)=\begin{bmatrix}\partial_{s}R_{i}(s,\theta)&0\end{bmatrix}$
for $i\in\\{1,2\\}$. Then, for any $x\in\mathbb{R}^{m}$, $\mathbf{x}\in
H_{1}^{n}$,

$\displaystyle\partial_{s}\left(\mathcal{P}{\footnotesize\begin{bmatrix}P,&Q_{1}\\\
Q_{2},&\left\\{R_{i}\right\\}\end{bmatrix}}\begin{bmatrix}x\\\
\mathbf{x}\end{bmatrix}\right)=\left(\partial_{s}\mathcal{P}{\footnotesize\begin{bmatrix}P,&Q_{1}\\\
Q_{2},&\left\\{R_{i}\right\\}\end{bmatrix}}\right)\begin{bmatrix}x\\\
\mathbf{x}\\\ \partial_{s}\mathbf{x}\end{bmatrix}$

###### Proof.

The result can be easily derived using the Leibniz integral rule, stating that
for any $F\in H_{1}[a,b]^{2}$,

$\displaystyle\frac{d}{ds}\left(\int_{L(s)}^{U(s)}F(s,\theta)d\theta\right)=F(s,U(s))\frac{d}{ds}U(s)-F(s,L(s))\frac{d}{ds}+\int_{L(s)}^{U(s)}\frac{d}{ds}F(s,\theta)d\theta.$

For a similar result for PI operators in 2D, we refer to [4]. ∎

####  A.7 Matrix Parametrization of Positive Definite PI Operators

In order to be able to solve optimization programs involving PI operator
$\mathcal{P}$, we need to be able to enforce positivity constraints
$\mathcal{P}\succcurlyeq 0$. For this, we parameterize PI operators by
positive matrices, expanding them as $\mathcal{P}=\mathcal{Z}^{*}P\mathcal{Z}$
for a fixed operator $\mathcal{Z}$, and positive semidefinite matrix
$P\succcurlyeq 0$. The following theorem provides a sufficient condition for
positivity of a 4-PI operator defined as

$\displaystyle\bigl{(}\mathcal{P}\left[\begin{smallmatrix}P,&Q_{1}\\\
Q_{2},&\\{R_{i}\\}\end{smallmatrix}\right]\mathbf{x}\bigr{)}(s)=\left[\begin{array}[]{ll}Px_{0}\nobreak\leavevmode&+\
\int_{a}^{b}Q_{1}(s)\mathbf{x}_{1}(s)ds\\\ Q_{2}(s)x_{0}&+\
R_{0}(s)\mathbf{x}_{1}(s)+\int_{a}^{s}R_{1}(s,\theta)\mathbf{x}_{1}(\theta)d\theta+\int_{s}^{b}R_{2}(s,\theta)\mathbf{x}_{1}(\theta)d\theta\end{array}\right]$
(A.13)

for $\mathbf{x}=\begin{bmatrix}x_{0}\\\
\mathbf{x}_{1}\end{bmatrix}\in\begin{bmatrix}\mathbb{R}^{n_{0}}\\\
L_{2}^{n_{1}}[a,b]\end{bmatrix}$. This result allows us to parameterize a cone
of positive PI operators as positive matrices, implement LPI constraints as
LMI constraints, allowing us to solve optimization problems with PI operators
using semi-definite programming solvers.

######  Theorem 11.

For any functions $Z_{1}:[a,b]\to\mathbb{R}^{d_{1}\times n}$,
$Z_{2}:[a,b]\times[a,b]\rightarrow\mathbb{R}^{d_{2}\times n}$, if $g(s)\geq 0$
for all $s\in[a,b]$ and

$\displaystyle P$ $\displaystyle=T_{11}\int_{a}^{b}g(s)ds,$ $\displaystyle
Q(\eta)$
$\displaystyle=g(\eta)T_{12}Z_{1}(\eta)+\int_{\eta}^{b}g(s)T_{13}Z_{2}(s,\eta)\text{d}s+\int_{a}^{\eta}g(s)T_{14}Z_{2}(s,\eta)\text{d}s,$
$\displaystyle R_{1}(s,\eta)$
$\displaystyle=g(s)Z_{1}(s)^{\top}T_{23}Z_{2}(s,\eta)+g(\eta)Z_{2}(\eta,s)^{\top}T_{42}Z_{1}(\eta)+\int_{s}^{b}g(\theta)Z_{2}(\theta,s)^{\top}T_{33}Z_{2}(\theta,\eta)\text{d}\theta$
$\displaystyle\hskip 36.98866pt\leavevmode\nobreak\ \leavevmode\nobreak\
+\int_{\eta}^{s}g(\theta)Z_{2}(\theta,s)^{\top}T_{43}Z_{2}(\theta,\eta)\text{d}\theta+\int_{a}^{\eta}g(\theta)Z_{2}(\theta,s)^{\top}T_{44}Z_{2}(\theta,\eta)\text{d}\theta,$
$\displaystyle R_{2}(s,\eta)$
$\displaystyle=g(s)Z_{1}(s)^{\top}T_{32}Z_{2}(s,\eta)+g(\eta)Z_{2}(\eta,s)^{\top}T_{24}Z_{1}(\eta)+\int_{\eta}^{b}g(\theta)Z_{2}(\theta,s)^{\top}T_{33}Z_{2}(\theta,\eta)\text{d}\theta$
$\displaystyle\hskip 36.98866pt\leavevmode\nobreak\ \leavevmode\nobreak\
+\int_{s}^{\eta}g(\theta)Z_{2}(\theta,s)^{\top}T_{34}Z_{2}(\theta,\eta)\text{d}\theta+\int_{a}^{s}g(\theta)Z_{2}(\theta,s)^{\top}T_{44}Z_{2}(\theta,\eta)\text{d}\theta,$
$\displaystyle R_{0}(s)$ $\displaystyle=g(s)Z_{1}(s)^{\top}T_{22}Z_{1}(s).$
(A.14)

where

$\displaystyle T=\begin{bmatrix}T_{11}&T_{12}&T_{13}&T_{14}\\\
T_{21}&T_{22}&T_{23}&T_{24}\\\ T_{31}&T_{32}&T_{33}&T_{34}\\\
T_{41}&T_{42}&T_{43}&T_{44}\end{bmatrix}\succcurlyeq 0,$

then the operator $\mathcal{P}{\footnotesize\begin{bmatrix}P,&Q_{1}\\\
Q_{2},&\left\\{R_{i}\right\\}\end{bmatrix}}$ as defined in (A.13) is positive
semidefinite, i.e.
$\left\langle{\mathbf{x}},{\mathcal{P}{\footnotesize\begin{bmatrix}P,&Q_{1}\\\
Q_{2},&\left\\{R_{i}\right\\}\end{bmatrix}}\mathbf{x}}\right\rangle\geq 0$ for
all $\mathbf{x}\in\mathbb{R}^{m}\times L_{2}^{n}[a,b]$.

To see the PIETOOLS implementation, check Section 7.2.2. For an extension of
this result to 2D PI operators, we refer to [4].

###  Appendix B Troubleshooting

This appendix is dedicated to tackling issues related to installation and
setting up of PIETOOLS 2022. Additionally, we discuss some common issues users
may run into while solving LPIs.

####  B.1 Troubleshooting: Installation

When installing the PIETOOLS toolbox using the install script, in rare
circumstances, the user may run into one of following errors. In case of an
error, a message is displayed explaining the issue, which can be one of the
following.

  1. 1.

An error appeared when trying to create the folder …
Check if the folder already has a folder named PIETOOLS_2022. If that does not
resolve the issue try running MATLAB from an administrator account which has
the authority to create or modify folders. As a last resort, try installing in
a different folder. If that does not fix this, contact us.

  2. 2.

The installation directory “ ” already exists…
Obvious error. However, if there is no folder with the name PIETOOLS_2022,
check for a hidden folder.

  3. 3.

‘tbxmanager’ or ‘SeDuMi’ were not downloaded or installed…
Check your internet connection. Verify that MATLAB is allowed to download
files using the internet connection. Check if the websites for tbxmanager and
SeDuMi are operational.

  4. 4.

‘PIETOOLS’ was not downloaded or installed…
Check the suggestions for the previous error. If that does not fix the issue,
contact us.

  5. 5.

Could not modify the initialization file “startup.m”…
Try running MATLAB from an administrator account which has the authority to
create or modify folders. If that does not fix the issue, manually add SeDuMi
or the relevant SDP solver (like MOSEK, sdpt3, sdpnal) to your MATLAB path,
and extract the files in PIETOOLS_2022.zip. Also add PIETOOLS to your MATLAB
path.

  6. 6.

Could not save the path to a default location…
Try running MATLAB from an administrator account which has the authority to
create or modify folders. If that does not fix it, just add PIETOOLS and
SeDuMi to your MATLAB path and skip this step.

####  B.2 Troubleshooting: Solving LPIs

PIETOOLS can be used for solving LPI optimization problems and users may run
into errors while setting up and solving them. For any errors in setting up an
LPI, refer to the function and script headers to ensure that input-output
formats are correct. Ensure that PI objects, defined in the LPI problem, are
well-defined and valid PI operators. You can use the isvalid function to check
if a PI operator is well defined. If that does not fix the problem, feel free
to contact us.
PIETOOLS 2022 relies on recently released SOSTOOLS 4.00 with solver SeDuMi to
solve optimization problems. If you are unfamiliar with SOSTOOLS, and are
unsure how to intepret the results of a solved optimization problem, please
check the SOSTOOLS manual available at http://www.cds.caltech.edu/sostools/,
or the SeDuMi manual avaialable at https://sedumi.ie.lehigh.edu/?page_id=58
for more information. A brief overview of how to interpret results, and what
to do in case of error, is provided below:

  1. 1.

How do I interpret the results of a solved optimization problem?
A general rule of thumb is to look at: pinf, dinf, feasratio and Residual
norm. pinf and dinf should be 0, while feasratio is in between -1 and 1
(preferably closer to 1). The lower the residual norm the better. Refer to
SeDuMi manual to interpret other output parameters and more details [10].

  2. 2.

What if pinf is 1?
Verify if the LPI constraints are in fact feasible. Verify if the sign-
definiteness of the PI operator is on a compact interval (use psatz term if
local sign-definite is needed) or the entire real line. Use ’getdeg’ function
to check if your LPI constraint has high degree polynomials. If yes, make sure
that all opvar variables used in lpi_eq function have high enough degrees to
match it. If this does not resolve the issue contact us and attach the files
that you are trying to run along with a snapshot of the error/output.

  3. 3.

What if dinf is 1 and feasratio is -1?
This issue typically occurs when the objective function is unbounded from
below and becomes $-\infty$. Check if the objective function is bounded below.
If this does not resolve the issue contact us by email and attach the files
that you are trying to run along with a snapshot of the error/output.

####  B.3 Contact Details

To resolve issues, report bugs or to collaborate on any development work
regarding PIETOOLS, please contact us through email and we will get back to
you as quickly as possible. In case of issues with installation, solving
problems or bugs identified, please include the script file that generates the
error along with images of the error generated in MATLAB. You can reach us
through email at<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>and
<EMAIL_ADDRESS>
### Bibliography

  * [1] A. Das, S. Shivakumar, S. Weiland, and M. Peet.  ${H}_{\infty}$ optimal estimation for linear coupled PDE systems.  In Proceedings of the IEEE Conference on Decision and Control, 2019\.
  * [2] John Doyle, Andy Packard, and Kemin Zhou.  Review of LFTs, LMIs, and $\mu$.  In [1991] Proceedings of the 30th IEEE Conference on Decision and Control, pages 1227–1232 vol.2. IEEE, 1991.
  * [3] Israel Gohberg, Seymour Goldberg, and Marius A Kaashoek.  Classes of linear operators, volume 63.  Birkhäuser, 2013.
  * [4] Declan S. Jagt and Matthew M. Peet.  A PIE representation of coupled 2D PDEs and stability analysis using LPIs.  In Proceedings of the American Control Conference, 2021.
  * [5] M. Peet.  A partial integral equation representation of coupled linear PDEs and scalable stability analysis using LMIs.  Automatica, 2020.  Submitted.
  * [6] Stephen Prajna, Antonis Papachristodoulou, and Pablo A Parrilo.  Introducing SOSTOOLS: A general purpose sum of squares programming solver.  In Proceedings of the 41st IEEE Conference on Decision and Control, 2002., volume 1, pages 741–746. IEEE, 2002.
  * [7] Peter Seiler.  SOSOPT: A toolbox for polynomial optimization.  arXiv preprint arXiv:1308.1889, 2013.
  * [8] S. Shivakumar, A. Das, S. Weiland, and M. Peet.  A generalized LMI formulation for input-output analysis of linear systems of ODEs coupled with PDEs.  In Proceedings of the IEEE Conference on Decision and Control, 2019\.
  * [9] Sachin Shivakumar, Amritam Das, Siep Weiland, and Matthew M Peet.  Duality and ${H}_{\infty}$ -optimal control of coupled ODE-PDE systems.  arXiv preprint arXiv:2004.03638, 2020.
  * [10] Jos F Sturm.  Using SeDuMi 1.02, a MATLAB toolbox for optimization over symmetric cones.  In Optimization methods and software, volume 11, pages 625–653. Taylor & Francis, 1999.