$D_{\mathbf{x}_{0}}\boldsymbol{\phi}\left(\mathbf{F}^{t}(\mathbf{x}_{0})\right)=D\boldsymbol{\phi}\left(\mathbf{F}^{t}(\mathbf{x}_{0})\right)D\mathbf{F}^{t}(\mathbf{x}_{0}),$

and $\mathbf{F}^{t}$ is a diffeomorphism, we can rewrite (79) as

$\frac{d}{dt}\left\\{\text{$\mathcal{\boldsymbol{\mathcal{K}}}^{t}\left[\boldsymbol{\phi}\right]$}(\mathbf{x}_{0})\right\\}=D_{\mathbf{x}_{0}}\boldsymbol{\phi}\left(\mathbf{F}^{t}(\mathbf{x}_{0})\right)\left[D\mathbf{F}^{t}(\mathbf{x}_{0})\right]^{-1}\mathbf{f}\left(\mathbf{F}^{t}(\mathbf{x}_{0})\right)=D_{\mathbf{x}_{0}}\left\\{\mathcal{\boldsymbol{\mathcal{K}}}^{t}\left[\boldsymbol{\phi}\right](\mathbf{x}_{0})\right\\}\left[D\mathbf{F}^{t}(\mathbf{x}_{0})\right]^{-1}\mathbf{f}\left(\mathbf{F}^{t}(\mathbf{x}_{0})\right).$
(80)

Note, however, that

$\left[D\mathbf{F}^{t}(\mathbf{x}_{0})\right]^{-1}\mathbf{f}\left(\mathbf{F}^{t}(\mathbf{x}_{0})\right)=\mathbf{f}(\mathbf{x}_{0}),$
(81)

where we used the fact that
$\mathbf{f}\left(\mathbf{F}^{t}(\mathbf{x}_{0})\right)$ is a solution of the
equation of variations
$\dot{\boldsymbol{\xi}}=D\mathbf{F}^{t}(\mathbf{x}_{0})\boldsymbol{\xi}$ and
hence we have
$\mathbf{f}\left(\mathbf{F}^{t}(\mathbf{x}_{0})\right)=D\mathbf{F}^{t}(\mathbf{x}_{0})\mathbf{f}\left(\mathbf{x}_{0}\right)$.
Therefore, we obtain from eqs. (80) and (81) that

$\frac{d}{dt}\left(\mathcal{\boldsymbol{\mathcal{K}}}^{t}\boldsymbol{\phi}\right)=\boldsymbol{\mathcal{L}}\mathcal{\boldsymbol{\mathcal{K}}}^{t}\boldsymbol{\phi}$
(82)

with the Liouville operator
$\boldsymbol{\mathcal{L}}:\mathcal{G}^{d}\mapsto\mathcal{G}^{d}$ defined as

$\boldsymbol{\mathcal{L}}\colon\boldsymbol{\phi}\mapsto\mathcal{\boldsymbol{L}}\boldsymbol{\phi}:=D\boldsymbol{\phi}\,\mathbf{f}.$

###  A.3 An example: The phase space variable as observable

As an illustration of the solution structure of the functional differential
equation (82), consider the simplest nontrivial case wherein the observable
function $\boldsymbol{\phi}$ is just the identity map on
$\mathcal{P}=\mathbb{R}^{d}$, i.e.,
$\hat{\boldsymbol{\phi}}=\mathbf{I}\colon\mathbb{R}^{d}\to\mathbb{R}^{d}$. In
that case, we specifically have the expressions

$\hat{\boldsymbol{\phi}}(\mathbf{x}_{0})=\mathbf{x}_{0},\qquad\mathcal{\boldsymbol{\mathcal{K}}}^{t}\hat{\boldsymbol{\phi}}(\mathbf{x}_{0})=\mathbf{F}^{t}(\mathbf{x}_{0}),\qquad\frac{d}{dt}\left\\{\text{$\mathcal{\boldsymbol{\mathcal{K}}}^{t}\hat{\boldsymbol{\phi}}$}\right\\}=\dot{\mathbf{F}}^{t}(\mathbf{x}_{0})=\dot{\mathbf{x}},\qquad
D_{\mathbf{x}_{0}}\left\\{\mathcal{\boldsymbol{\mathcal{K}}}^{t}\hat{\boldsymbol{\phi}}(\mathbf{x}_{0})\right\\}=D\mathbf{F}^{t}(\mathbf{x}_{0}).$

Substituting these formulas into the functional differential equation (82), we
obtain

$\dot{\mathbf{x}}=D\mathbf{F}^{t}(\mathbf{x}_{0})\mathbf{f}(\mathbf{x}_{0})=\mathbf{f}(\mathbf{x}).$
(83)

Therefore, for the most commonly used observable in classic dynamical systems,
the phase space variable $\mathbf{x}$, the linear functional differential
equation (82) simplifies to the original nonlinear ODE ((1)).

This again underscores that the linearity of the Koopman operator does not
imply linear dynamics for individual observations
$\mathbf{x}(t;\mathbf{x}_{0})$ of the evolution of the initial condition
$\mathbf{x}_{0}$. Rather, it implies that for any real constant $c$, _the
scaled initial observation
$\tilde{\boldsymbol{\phi}}\,(\mathbf{x}_{0})=c\hat{\boldsymbol{\phi}}(\mathbf{x}_{0})$
of $\mathbf{x}_{0}$ will evolve into the identically scaled current
observation
_$\mathcal{\boldsymbol{\mathcal{K}}}^{t}\hat{\boldsymbol{\phi}}(\mathbf{x}_{0})=c\mathcal{\boldsymbol{\mathcal{K}}}^{t}\hat{\boldsymbol{\phi}}(\mathbf{x}_{0})=c\mathbf{x}(t;\mathbf{x}_{0})$
_of the trajectory at time $t$_, i.e.,
$\mathcal{\boldsymbol{\mathcal{K}}}^{t}$ will return
$\mathbf{x}(t;\mathbf{x}_{0})$ scaled by the same constant $c$. Most
importantly, we will still generally have

$\mathbf{x}(t;c\mathbf{x}_{0})\neq c\mathbf{x}(t;\mathbf{x}_{0}),$

unless the dynamical system (83) is linear. Therefore, the dynamics of the
fixed observable $\mathbf{x}(t;\mathbf{x}_{0})$ is not linear in the usual
sense, i.e., with respect to changes in initial conditions. Rather, it is
linear with respect to changes in the observable in which the evolution of the
same trajectory is observed.

Equation (83) also illustrates that passing from the autonomous nonlinear ODE
(1) to the functional differential equation (82) gives as infinite-dimensional
linear formulation that is at least as complicated to solve as the original
finite-dimensional nonlinear ODE formulation.

###  A.4 Koopman eigenfunctions

There are, nevertheless, non-generic sets of observables restricted to which
e.q. (82) becomes a constant-coefficient linear system of ODEs. Examples are
observables falling in the span of eigenfunctions of the Koopman operator
$\mathcal{\boldsymbol{\mathcal{K}}}^{t}$, if and where such eigenfunctions
exist.

Indeed, if for some constant $\lambda\in\mathbb{C}$, the eigenvalue problem

$\mathcal{\boldsymbol{\mathcal{K}}}^{t}\left[\boldsymbol{\phi}\right](\mathbf{x}_{0})=e^{\lambda
t}\boldsymbol{\phi}(\mathbf{x}_{0}),\quad\mathbf{x}_{0}\in\mathcal{D},$ (84)

has a solution $\boldsymbol{\phi}\colon\mathcal{D}\to\mathbb{C}^{d}$ over some
domain $\mathcal{D}\subset\mathcal{P}$ of the phase space, then we obtain

$\frac{d}{dt}\text{$\mathcal{\boldsymbol{\mathcal{K}}}^{t}\left[\boldsymbol{\phi}\right]$}=\lambda\mathcal{\boldsymbol{\mathcal{K}}}^{t}\left[\boldsymbol{\phi}\right].$

Consequently, if
$\boldsymbol{\phi}_{1}(\boldsymbol{x}_{0}),\ldots,\boldsymbol{\phi}_{k}(\boldsymbol{x}_{0})\in\mathbb{C}^{d}$
are Koopman eigenfunctions with corresponding Koopman eigenvalues
$\lambda_{1},\ldots,\lambda_{k}\in\mathbb{C}$ and domains of definition
$\mathcal{D}_{1},\ldots,\mathcal{D}_{k}\subset\mathcal{P}$, then any
observable
$\boldsymbol{w}\in\mathcal{\Phi}^{d}\left(\cap_{j=1}^{k}\mathcal{D}_{j}\right)$
of the form

$\boldsymbol{w}(\mathbf{x}_{0})=c_{1}\boldsymbol{\phi}_{1}(\mathbf{x}_{0})+\ldots+c_{k}\boldsymbol{\phi}_{k}(\mathbf{x}_{0})$
(85)

satisfies

$\text{$\mathcal{\boldsymbol{\mathcal{K}}}^{t}\left[\boldsymbol{w}\right]$}=\sum_{j=1}^{k}c_{j}\mathcal{\boldsymbol{\mathcal{K}}}^{t}\left[\boldsymbol{\phi}^{j}\right]=\sum_{j=1}^{k}c_{j}e^{\lambda_{j}t}\boldsymbol{\phi}_{j}.$

Therefore, the coordinate representation $\boldsymbol{W}(t)$ of
$\mathcal{\boldsymbol{\mathcal{K}}}^{t}\left[\boldsymbol{w}\right]$ in the
basis $\boldsymbol{\phi}_{1},\ldots,\boldsymbol{\phi}_{j},$ given by

$\mathbf{W}(t)=\left(\begin{array}[]{c}c_{1}e^{\lambda_{1}t}\\\ \vdots\\\
c_{k}e^{\lambda_{k}t}\end{array}\right),$

satisfies the $k$-dimensional autonomous linear ODE

$\frac{d}{dt}\mathbf{W}=\boldsymbol{\Lambda}\mathbf{W},\quad\boldsymbol{\Lambda}=\left(\begin{array}[]{ccc}\lambda_{1}&0&0\\\
0&\ddots&0\\\ 0&0&\lambda_{k}\end{array}\right).$

More generally, if $\boldsymbol{\phi}_{1},\ldots,\boldsymbol{\phi}_{k}$ are
just linearly independent observables in the spectral subspace
$\mathrm{span}\left\\{\boldsymbol{\phi}_{1},\ldots,\boldsymbol{\phi}_{k}\right\\}$,
then a similar ODE holds for $W$ in the basis
$\boldsymbol{\phi}_{1},\ldots,\boldsymbol{\phi}_{k}.$ In that case, $\Lambda$
is not a diagonal matrix but its eigenvalues are still
$\lambda_{1},\ldots,\lambda_{k}$.

It is often forgotten, however, that Koopman eigenfunctions satisfying the
eigenvalue problem (84) generally only exist on a subset $\mathcal{D}$ of the
phase space. This is made precise by the following simple observation:

######  Proposition 1.

All $C^{\infty}$ Koopman eigenfunctions blow up (i.e., become unbounded) at
the boundaries of domains of attraction or repulsion of a compact invariant
set along which the underlying dynamical system admits a $C^{1}$linearization
transformation.

###### Proof.

We only prove the statement for domains of attraction, as domains of repulsion
can be handled in the same fashion in backward time. Let $\mathbf{x}_{0}$ be a
point on the boundary $\mathcal{B}$ of a domain of attraction of a compact
invariant manifold $\mathcal{S}$. Examples of such include attracting fixed
points, periodic orbits and quasiperiodic tori. As the domain of attraction is
open, for any $\rho>0$, there is a closed $\rho$-neighborhood
$B_{\rho}(\mathcal{S})$ around the invariant set $\mathcal{S}$ that lies fully
in the domain of attraction. Consequently, we have
$d(\mathcal{S},\mathbf{F}^{t}(\mathbf{x}_{0}))>\rho/2>0$ for all $t\geq 0$. By
continuous dependence of the flow map on initial conditions, for any $T>0$
there exists $\delta(T)>0$ such that for all $\hat{\mathbf{x}}_{0}\in
B_{\delta(T)}\left(\mathbf{x}_{0}\right)\cap\mathcal{B}$, we have

$d(\mathbf{F}^{t}(\hat{\mathbf{x}}_{0}),\mathbf{F}^{t}(\mathbf{x}_{0}))<\rho/2,\quad
t\in\left[0,T\right].$

This means that for arbitrary large $T$, no initial condition in
$B_{\delta(T)}\left(\mathbf{x}_{0}\right)\cap\mathcal{B}$ will reach the ball
$B_{\rho}(\mathcal{S})$ around the origin for times up to $T$. This in turn
means that for any $\mathbf{x}_{0}\in\mathcal{B},$ there are arbitrarily close
initial conditions $\hat{\mathbf{x}}_{0}$ in the domain of attraction with
arbitrarily long time of flight to the ball $B_{\rho}(\mathcal{S})$ around the
attracting invariant set $\mathcal{S}.$

For the special case of a fixed point $p=\mathcal{S}$, Lan and Mezic [35]
extended the local $C^{1}$ linearization result of Hartman [23] near an
attracting fixed point via the backward-time flow map to a global $C^{1}$
linearizing transformation $\mathbf{x}=\mathbf{h}(\mathbf{y})$ within
$\mathcal{B}$ that maps the original ODE within $\mathcal{B}$ to
$\dot{\mathbf{y}}=\mathbf{A}\mathbf{y}$ defined on $\mathbb{R}^{n}$. This
transformation has also been shown unique by Eldering et al. [17]. By the
construction of this unique linearizing transformation, initial conditions
$\hat{\mathbf{x}}_{0}$ whose times of flight to $B_{\rho}(p)$ is large will
then map those initial conditions $\hat{\mathbf{y}}_{0}$ of the linearized
system that are far from the origin $\mathbf{y}=\mathbf{0}$.

The unique principal Koopman eigenfunctions of the nonlinear system are known
to be

$\boldsymbol{\phi}_{j}(\mathbf{x})=\left\langle\mathbf{h}^{-1}(\mathbf{x}),\mathbf{v}_{j}\right\rangle,\quad
j=1,\ldots,n.$ (86)

As we have just concluded that $\mathbf{y}=\mathbf{h}^{-1}(\mathbf{x})$ will
admit arbitrarily large values arbitrarily close to any point
$\mathbf{x}_{0}\in\mathcal{B},$ the principal Koopman eigenfunctions
$\boldsymbol{\phi}_{j}(\mathbf{x})$ will also have arbitrarily large values
near any $\mathbf{x}_{0}\in\mathcal{B}$ by eq. (86). (We need to exclude
points $\hat{\mathbf{x}}_{0}$ at which
$\mathbf{h}^{-1}(\hat{\mathbf{x}_{0}})\perp\mathbf{v}_{j}$, but those points
form an $n-1$-dimensional plane $\mathcal{V}(\hat{\mathbf{x}}_{0})$ that has a
measure zero intersection with
$B_{\delta(T)}\left(\mathbf{x}_{0}\right)\cap\mathcal{B}$.) Therefore, the set
$B_{\delta(T)}\left(\mathbf{x}_{0}\right)-\mathcal{V}(\hat{\mathbf{x}}_{0})$
still has arbitrarily many points with arbitrarily large values of
$\boldsymbol{\phi}_{j}(\mathbf{x})$. This establishes the blow-up of any
principal Koopman eigenfunction $\boldsymbol{\phi}_{j}(\mathbf{x})$ at any of
the domain boundary points $\mathbf{x}_{0}\in\mathcal{B}$.

Kvalheim and Revzen [34] show that all $C^{\infty}$ Koopman eigenfunctions are
products of powers of principal Koopman eigenfunctions. As a consequence, we
obtain that all $C^{\infty}$ Koopman eigenfunctions blow up a the boundary
$\mathcal{B}$ of the domain of attraction of a fixed point. Finally, Eldering
et al. [17] extend the results of Lan and Mezic [35] from attracting fixed
points $p$ to arbitrary compact invariant manifolds $\mathcal{S}$ with a
domain of attraction, as long as a local $C^{1}$ linearizing transformation
exists near $\mathcal{S}.$ ∎

As a consequence of Proposition 1, one cannot simply patch together Koopman
eigenfunctions outside boundaries of domains of attraction and thereby provide
a viable Koopman-linearization for the whole phase space, as is often
suggested. Indeed, even if found somehow, Koopman eigenfunctions would become
unmanageable for the purposes of observer expansions even before reaching
basin boundaries, as shown explicitly by the examples of Page and Kerswell
[44].

##  Appendix B Proof of Theorem 1 and technical remarks

Under assumption (A1) of hyperbolicity of the theorem and under assumption
(A2) on the smoothness class $C^{2}$ of the dynamical system, a refinement of
the Hartman–Grobman theorem (see Guckenheimer and Holmes [19]) by Hartman [23]
guarantees the local existence of a near-identity, linearizing change of
coordinates of the form555Hartman [23] only states that $h\in C^{1}$ in his
main theorem. He adds, however, that his proof of the theorem also shows that
$Dh$ is uniformly Hölder continuous in $U$ with Hölder exponent
$\beta\in\left(0,1\right)$, which implies eq. (87).

$\mathbf{x}=\mathbf{y}+\mathbf{h}(\mathbf{y}),\qquad\mathbf{y}\in
U\subset\mathbb{R}^{n},\qquad\mathbf{h}(\mathbf{y}),\mathbf{h}^{-1}(\mathbf{y})=\mathcal{O}\left(\left|\mathbf{y}\right|^{1+\beta}\right),$
(87)

for some $\beta\in\left(0,1\right)$, under which system (11) takes the exact
linear form

$\dot{\mathbf{y}}=\mathbf{A}\mathbf{y}.$ (88)

As pointed out by Haller et al. [22] in the context of the more restrictive,
$C^{\infty}$ linearization theorem of Sternberg [53], the inverse of the
linearizing transformation (87) maps the $d$-dimensional spectral subspace $E$
of the linearized system (88) into a $d$-dimensional, normally attracting
spectral submanifold (SSM),

$\mathcal{W}(E)=\left(\mathbf{I}+\mathbf{h}\right)^{-1}(E),$ (89)

for system (11). The internal dynamics of this SSM govern the longer-term
behavior of all trajectories near the origin. Generally, there will be
infinitely many such invariant manifolds, but they are all tangent to $E$ at
the origin. We note that, $\mathcal{W}(E)$, as defined in eq. (89), is only
known to be $C^{1}$ under our current set of assumptions. Under further
nonresonance assumption on the spectrum of $\mathbf{A},$ the SSM can be shown
to be as smooth as the underlying dynamical system (see Cabré et al. [10],
Haller and Ponsioen [21]).

Based on these results, we apply the linearizing transformation (87) followed
by a linear change of coordinates to
$\left(\boldsymbol{\xi},\boldsymbol{\eta}\right)$, where $\boldsymbol{\xi}$
are coordinates along the spectral subspace $E$ and $\boldsymbol{\eta}$ are
coordinates along the spectral subspace $F.$ Specifically, let

$\displaystyle\mathbf{x}$
$\displaystyle=\mathbf{T}\mathbf{y}+\mathbf{h}(\mathbf{T}\mathbf{y}),\quad\mathbf{y}=\left(\boldsymbol{\xi},\boldsymbol{\eta}\right)^{\mathrm{T}}\in\mathbb{R}^{d}\times\mathbb{R}^{n-d},$
(90)

with the matrix $\mathbf{T}$ defined in formula (17).

Under the transformation (90), the nonlinear system (11)becomes

$\left(\begin{array}[]{c}\dot{\boldsymbol{\xi}}\\\
\dot{\boldsymbol{\eta}}\end{array}\right)=\left(\begin{array}[]{cc}\boldsymbol{\Lambda}_{E}&0\\\
0&\boldsymbol{\Lambda}_{F}\end{array}\right)\left(\begin{array}[]{c}\boldsymbol{\xi}\\\
\boldsymbol{\eta}\end{array}\right),\qquad\boldsymbol{\Lambda}_{E}=\left(\mathbf{T}^{-1}\mathbf{A}\mathbf{T}\right)|_{E},\qquad\boldsymbol{\Lambda}_{F}=\left(\mathbf{T}^{-1}\mathbf{A}\mathbf{T}\right)|_{F}.$
(91)

In these coordinates, the smoothest invariant manifold tangent to $E$ is

$\mathcal{W}(E)=\left\\{\left(\boldsymbol{\xi},\boldsymbol{\eta}\right)\,:\,\boldsymbol{\eta}=0\right\\},$

within which the dynamics restricted to $\mathcal{W}(E)$ are then given by

$\dot{\boldsymbol{\xi}}=\boldsymbol{\Lambda}_{E}\boldsymbol{\xi}.$ (92)

From eq. (91), we obtain the estimates

$\displaystyle\left|\boldsymbol{\xi}\left(t;\boldsymbol{\xi}_{0}\right)\right|$
$\displaystyle\leq\left|\boldsymbol{\xi}_{0}\right|e^{\left(\mathrm{Re}\text{$\lambda_{1}+\epsilon_{1}$}\right)t},$
$\displaystyle\left|\boldsymbol{\eta}\left(t;\boldsymbol{\eta}_{0}\right)\right|$
$\displaystyle\leq\left|\boldsymbol{\eta}_{0}\right|e^{\left(\mathrm{Re}\text{$\lambda_{d+1}+\epsilon_{d+1}$}\right)t},$

where $\epsilon_{j}\geq 0$ is an arbitrarily small constant that can be chosen
zero if the algebraic multiplicity of $\lambda_{j}$ is equal to its geometric
multiplicity.

With the help of the $\left(\boldsymbol{\xi},\boldsymbol{\eta}\right)$
coordinates and eq. (87), the class $C^{2}$ observable $\boldsymbol{\phi}$ be
locally be written as

$\displaystyle\boldsymbol{\phi}(x)$
$\displaystyle=\boldsymbol{\phi}\left(\mathbf{T}\mathbf{y}+\mathbf{h}(\mathbf{T}\mathbf{y})\right)=D\boldsymbol{\phi}(\mathbf{0})\left[\mathbf{T}\mathbf{y}+\mathbf{h}(\mathbf{Ty})\right]+\mathcal{O}\left(\left|\mathbf{y}+\mathbf{h}(\mathbf{y})\right|^{2}\right)$
$\displaystyle=D\boldsymbol{\phi}(\mathbf{0})\mathbf{T}\left(\begin{array}[]{c}\boldsymbol{\xi}\\\
\boldsymbol{\eta}\end{array}\right)+D\boldsymbol{\phi}(\mathbf{0})\mathbf{h}(\mathbf{T}\mathbf{y})+\mathcal{O}\left(\left|\boldsymbol{\xi}\right|^{2},\left|\boldsymbol{\xi}\right|\left|\boldsymbol{\eta}\right|,\left|\boldsymbol{\eta}\right|^{2}\right)$
(95)
$\displaystyle=D\boldsymbol{\phi}(\mathbf{0})\mathbf{T}_{E}\boldsymbol{\xi}+\mathcal{O}\left(\left|\boldsymbol{\eta}\right|\right)+\mathcal{O}\left(\left|\left(\boldsymbol{\xi},\boldsymbol{\eta}\right)\right|^{1+\beta}\right)+\mathcal{O}\left(\left|\boldsymbol{\xi}\right|^{2},\left|\boldsymbol{\xi}\right|\left|\boldsymbol{\eta}\right|,\left|\boldsymbol{\eta}\right|^{2}\right)$
(96)
$\displaystyle=D\boldsymbol{\phi}(\mathbf{0})\mathbf{T}_{E}\boldsymbol{\xi}+\mathcal{O}\left(\left|\boldsymbol{\eta}\right|\right)+\mathcal{O}\left(\left|\left(\boldsymbol{\xi},\boldsymbol{\eta}\right)\right|^{1+\beta}\right)$

for some $\beta\in\left(0,1\right)$.

We also express the data matrix $\boldsymbol{\Phi}$ defined in eq. (5) in
terms of data matrices with respect to the coordinates
$\mathbf{y}=\left(\boldsymbol{\xi},\boldsymbol{\eta}\right)$ by letting

$\boldsymbol{\Phi}=\boldsymbol{\phi}\left(\mathbf{TY}-\mathbf{h}(\mathbf{TY})\right),\qquad\mathbf{Y}=\left(\begin{array}[]{c}\boldsymbol{\Xi}\\\
\boldsymbol{H}\end{array}\right),$ (97)

where the functions $\boldsymbol{\phi}$ and $\boldsymbol{h}$ are applied to
the matrices involved column by column. From eqs. (96) and (97), we obtain
that

$\displaystyle\Phi$
$\displaystyle=D\boldsymbol{\phi}(\mathbf{0})\mathbf{T}_{E}\boldsymbol{\Xi}+\mathcal{O}\left(\left|\boldsymbol{H}\right|\right)+\mathcal{O}\left(\left|\left(\boldsymbol{\Xi},\boldsymbol{H}\right)\right|^{1+\beta}\right),$
$\displaystyle\hat{\Phi}$
$\displaystyle=D\boldsymbol{\phi}(\mathbf{0})\mathbf{T}_{E}e^{\boldsymbol{\Lambda}_{E}\Delta
t}\boldsymbol{\Xi}+\mathcal{O}\left(\left|\boldsymbol{H}\right|\right)+\mathcal{O}\left(\left|\left(\boldsymbol{\Xi},\boldsymbol{H}\right)\right|^{1+\beta}\right),$
(98)

where we have used the boundedness of the sampling time $\Delta t$ and the
near-identity nature of the linearizing mapping (90), which implies

$\mathcal{O}\left(\left|\boldsymbol{\Xi}\right|^{\beta}\right)=\mathcal{O}\left(\left|\mathbf{P}_{E}\mathbf{X}\right|^{\beta}\right),\qquad\mathcal{O}\left(\left|\boldsymbol{H}\right|^{\beta}\right)=\mathcal{O}\left(\left|\mathbf{P}_{F}\mathbf{X}\right|^{\beta}\right).$
(99)

From the equations (98), we obtain

$\displaystyle\hat{\boldsymbol{\Phi}}\boldsymbol{\Phi}^{\mathrm{T}}$
$\displaystyle=\left[D\boldsymbol{\phi}(0)\mathbf{T}_{E}e^{\boldsymbol{\Lambda}_{E}\Delta
t}\boldsymbol{\Xi}+\mathcal{O}\left(\left|\boldsymbol{H}\right|\right)+\mathcal{O}\left(\left|\left(\boldsymbol{\Xi},\boldsymbol{H}\right)\right|^{1+\beta}\right)\right]\left[D\boldsymbol{\phi}(0)\mathbf{T}_{E}\boldsymbol{\Xi}+\mathcal{O}\left(\left|\boldsymbol{H}\right|\right)+\mathcal{O}\left(\left|\left(\boldsymbol{\Xi},\boldsymbol{H}\right)\right|^{1+\beta}\right)\right]^{\mathrm{T}}$
$\displaystyle=D\boldsymbol{\phi}(0)\mathbf{T}_{E}e^{\boldsymbol{\Lambda}_{E}\Delta
t}\boldsymbol{\Xi}\boldsymbol{\Xi}^{\mathrm{T}}\mathbf{T}_{E}^{\mathrm{T}}\left[D\boldsymbol{\phi}(\mathbf{0})\right]^{\mathrm{T}}+\mathcal{O}\left(\left|\boldsymbol{H}\right|^{2},\left|\boldsymbol{H}\right|\left|\left(\boldsymbol{\Xi},\boldsymbol{H}\right)\right|^{1+\beta},\left|\left(\boldsymbol{\Xi},\boldsymbol{H}\right)\right|^{2+2\beta}\right),$
$\displaystyle\left(\boldsymbol{\Phi}\boldsymbol{\Phi}^{\mathrm{T}}\right)^{\dagger}$
$\displaystyle=\left\\{\left[D\boldsymbol{\phi}(\mathbf{0})\mathbf{T}_{E}\boldsymbol{\Xi}+\mathcal{O}\left(\left|\boldsymbol{H}\right|\right)+\mathcal{O}\left(\left|\left(\boldsymbol{\Xi},\boldsymbol{H}\right)\right|^{1+\beta}\right)\right]\left[D\boldsymbol{\phi}(\mathbf{0})\mathbf{T}_{E}\boldsymbol{\Xi}+\mathcal{O}\left(\left|\boldsymbol{H}\right|\right)+\mathcal{O}\left(\left|\left(\boldsymbol{\Xi},\boldsymbol{H}\right)\right|^{1+\beta}\right)\right]^{\mathrm{T}}\right\\}^{\dagger}$
$\displaystyle=\left[D\boldsymbol{\phi}(\mathbf{0})\mathbf{T}_{E}\boldsymbol{\Xi}\boldsymbol{\Xi}^{\mathrm{T}}\mathbf{T}_{E}^{\mathrm{T}}\left[D\boldsymbol{\phi}(\mathbf{0})\right]^{\mathrm{T}}+\mathcal{O}\left(\left|\boldsymbol{H}\right|^{2},\left|\boldsymbol{H}\right|\left|\left(\boldsymbol{\Xi},\boldsymbol{H}\right)\right|^{1+\beta},\left|\left(\boldsymbol{\Xi},\boldsymbol{H}\right)\right|^{2+2\beta}\right)\right]^{\dagger}.$
(100)

By the first assumption in (A4) of the theorem,
$\boldsymbol{\Phi}\boldsymbol{\Phi}^{\mathrm{T}}$ is invertible, and hence
$\left(\boldsymbol{\Phi}\boldsymbol{\Phi}^{\mathrm{T}}\right)^{\dagger}=\left(\boldsymbol{\Phi}\boldsymbol{\Phi}^{\mathrm{T}}\right)^{-1}$.
Therefore, by the second equation in (100), for small enough
$\left|\boldsymbol{\Xi}\right|$ and $\left|\boldsymbol{H}\right|$, the matrix
$D\boldsymbol{\phi}(\mathbf{0})\mathbf{T}_{E}\boldsymbol{\Xi}\boldsymbol{\Xi}^{\mathrm{T}}\mathbf{T}_{E}^{\mathrm{T}}\left[D\boldsymbol{\phi}(\mathbf{0})\right]^{\mathrm{T}}$
is also invertible and

$\displaystyle\left(\boldsymbol{\Phi}\boldsymbol{\Phi}^{\mathrm{T}}\right)^{\dagger}$
$\displaystyle=\left[D\boldsymbol{\phi}(\mathbf{0})\mathbf{T}_{E}\boldsymbol{\Xi}\boldsymbol{\Xi}^{\mathrm{T}}\mathbf{T}_{E}^{\mathrm{T}}\left[D\boldsymbol{\phi}(\mathbf{0})\right]^{\mathrm{T}}\right]^{-1}+\mathcal{O}\left(\left|\boldsymbol{H}\right|^{2},\left|\boldsymbol{H}\right|\left|\left(\boldsymbol{\Xi},\boldsymbol{H}\right)\right|^{1+\beta},\left|\left(\boldsymbol{\Xi},\boldsymbol{H}\right)\right|^{2+2\beta}\right).$
(101)

Then, by assumption (A3) of the theorem, we can write

$\left(\boldsymbol{\Phi}\boldsymbol{\Phi}^{\mathrm{T}}\right)^{\dagger}=\left(\mathbf{T}_{E}^{\mathrm{T}}\left[D\boldsymbol{\phi}(\mathbf{0})\right]^{\mathrm{T}}\right)^{-1}\left(\boldsymbol{\Xi}\boldsymbol{\Xi}^{\mathrm{T}}\right)^{-1}\left(D\boldsymbol{\phi}(\mathbf{0})\mathbf{T}_{E}\right)^{-1}+\mathcal{O}\left(\left|\boldsymbol{H}\right|^{2},\left|\boldsymbol{H}\right|\left|\left(\boldsymbol{\Xi},\boldsymbol{H}\right)\right|^{1+\beta},\left|\left(\boldsymbol{\Xi},\boldsymbol{H}\right)\right|^{2+2\beta}\right).$
(102)

Substitution of formulas (101)-(102) into eq. (8) gives

$\displaystyle\boldsymbol{\mathcal{D}}$
$\displaystyle=\hat{\boldsymbol{\Phi}}\boldsymbol{\Phi}^{\mathrm{T}}\left(\boldsymbol{\Phi}\boldsymbol{\Phi}^{\mathrm{T}}\right)^{\dagger}=\hat{\boldsymbol{\Phi}}\boldsymbol{\Phi}^{\mathrm{T}}\left(\boldsymbol{\Phi}\boldsymbol{\Phi}^{\mathrm{T}}\right)^{-1}$
$\displaystyle=\left[D\boldsymbol{\phi}(0)\mathbf{T}_{E}e^{\boldsymbol{\Lambda}_{E}\Delta
t}\boldsymbol{\Xi}\boldsymbol{\Xi}^{\mathrm{T}}\mathbf{T}_{E}^{\mathrm{T}}\left[D\boldsymbol{\phi}(\mathbf{0})\right]^{\mathrm{T}}+\mathcal{O}\left(\left|\boldsymbol{H}\right|^{2},\left|\boldsymbol{H}\right|\left|\left(\boldsymbol{\Xi},\boldsymbol{H}\right)\right|^{1+\beta},\left|\left(\boldsymbol{\Xi},\boldsymbol{H}\right)\right|^{2+2\beta}\right)\right]$
$\displaystyle\times\left[\left(\mathbf{T}_{E}^{\mathrm{T}}\left[D\boldsymbol{\phi}(\mathbf{0})\right]^{\mathrm{T}}\right)^{-1}\left(\boldsymbol{\Xi}\boldsymbol{\Xi}^{\mathrm{T}}\right)^{-1}\left(D\boldsymbol{\phi}(\mathbf{0})\mathbf{T}_{E}\right)^{-1}+\mathcal{O}\left(\left|\boldsymbol{H}\right|^{2},\left|\boldsymbol{H}\right|\left|\left(\boldsymbol{\Xi},\boldsymbol{H}\right)\right|^{1+\beta},\left|\left(\boldsymbol{\Xi},\boldsymbol{H}\right)\right|^{2+2\beta}\right)\right]$
$\displaystyle=D\boldsymbol{\phi}(0)\mathbf{T}_{E}e^{\boldsymbol{\Lambda}_{E}\Delta
t}\left(D\boldsymbol{\phi}(\mathbf{0})\mathbf{T}_{E}\right)^{-1}+\mathcal{O}\left(\left|\boldsymbol{\Xi}\right|^{-2}\left|\boldsymbol{H}\right|^{2},\left|\boldsymbol{\Xi}\right|^{-2}\left|\boldsymbol{H}\right|\left|\left(\boldsymbol{\Xi},\boldsymbol{H}\right)\right|^{1+\beta},\left|\boldsymbol{\Xi}\right|^{-2}\left|\left(\boldsymbol{\Xi},\boldsymbol{H}\right)\right|^{2+2\beta}\right)$
$\displaystyle=D\boldsymbol{\phi}(0)\mathbf{T}_{E}e^{\boldsymbol{\Lambda}_{E}\Delta
t}\left(D\boldsymbol{\phi}(\mathbf{0})\mathbf{T}_{E}\right)^{-1}+\mathcal{O}\left(\left|\boldsymbol{\Xi}\right|^{\beta}\right),$

where we have use the second assumption in (A4) of the Theorem. This completes
the proof of Theorem 1, given the order-of-magnitude relations (99).

######  Remark 4.

In systems where the slowest decaying modes are oscillatory, the slow spectral
subspace $E$ is always even dimensional, as it is spanned by the real and
imaginary parts of the generalized eigenvectors of $\mathbf{A}$ that
correspond to complex conjugate pairs of eigenvalues. For such predominantly
oscillatory systems, therefore, the number $d$ of observables used in DMD has
to be an even number for assumption (A3) to hold.

######  Remark 5.

In most applications the dynamical system will have higher degree of
smoothness, i.e., we will have $\mathbf{f}\in C^{r}$ for some
_$r\in\mathbb{N}^{+}\cup\left\\{\infty,a\right\\}$ , _with $C^{a}$ referring
to the space of analytic functions. Then, under further nonresonance
conditions of the spectrum of $\mathbf{A}$, the spectral submanifold
$\mathcal{W}(E)$ will also be of class $C^{r}$, as we discuss in Section 4. In
those cases, for the approximate topological equivalence holds in the
statement of Theorem 1 hold on the full domain of attraction of the
$\mathbf{\mathbf{x}}=\mathbf{0}$ fixed point within $\mathcal{W}(E)$, as one
deduced from the linearization-in-the-large results of Lan and Mezic [35],
Eldering et al. [17].666Note that these results require at least $C^{2}$
differentiability for the dynamical system. Consequently, one cannot use them
within $\mathcal{W}(E)$ under the general assumptions of Theorem 1, which only
guarantee $C^{1}$ differentiability for the reduced flow in $\mathcal{W}(E)$.
In Theorem 1, our objective was to give a minimal set of conditions under
which DMD can be justified as an approximate, leading-order, $d$-dimensional
reduced model for the original nonlinear system. These minimal conditions only
require hyperbolicity, i.e., robustness of the spectrum under perturbations,
without insisting on the lack of resonances in the spectrum. This is
advantageous in data-driven applications in which details of the spectrum are
generally not known.

######  Remark 6.

As in the case of finite-dimensional systems, linearization theorems
guaranteeing higher-degree of smoothness for the linearizing transformations
are also available in infinite dimensions, but these require various non-
resonance conditions on the spectrum of $\mathbf{A}$ (see, e.g., Elbialy
[16]). In Theorem 3, as in the finite-dimensional case, our objective was to
present a minimal set of assumptions under which DMD can be justified without
the detailed knowledge of the spectrum of $\mathbf{A}$.

######  Remark 7.

We note that the more general Theorem 1.5 of Newhouse [41] also shows
$C^{1,\alpha}$ linearizability for a class of systems that have so called
_$\alpha-$ hyperbolic_ (as opposed to stable hyperbolic) fixed points. This
would enable us to waive the requirement in Theorem 3 that $\mathbf{A}$ is a
contraction and also allow for expanding directions within the spectral
subspace $E$, as we did in our remarks after Theorems 1 and 2. The defining
properties of these $\alpha$-hyperbolic systems are, however, complicated to
verify and do not hold for general hyperbolic fixed points even in finite-
dimensional Banach spaces. For instance, the simple three-dimensional example
of Hartman [23] (eq. (8) in that reference) of a non-$C^{1}$ linearizable
system has a hyperbolic fixed point that is not $\alpha$-hyperbolic in the
sense of Newhouse [41].

##  Appendix C Necessity of assumptions (A2)-(A4) of Theorem 1

Assumption (A2) of Theorem 1 requires $C^{2}$ smoothness for the dynamical
system (11). We now show on an example that if only $C^{1}$ smoothness holds
for the dynamical system, then the DMD approximation can be less accurate.
Consider the 1D system

$\dot{x}=-x+x^{1+\alpha}.$ (103)

The origin $x=0$ is an asymptotically stable fixed point. For $\alpha>0$, the
system is $C^{1}$ at the origin, but the second derivative
$\frac{\partial^{2}\dot{x}}{\partial x^{2}}=\alpha^{2}x^{\alpha-1}$ is
singular at $x=0$ for $\alpha<1$. We compare the accuracy of a simple DMD
prediction for $\alpha=1$, which corresponds to $C^{2}$ smoothness and for
$\alpha=0.9$, which corresponds to only $C^{1}$-smoothness at $x=0$. The DMD
prediction errors can be seen in Fig. 10 as a function of the initial distance
of the trajectory from the origin. Although the $C^{1}$-smooth system is close
to the $C^{2}$system, a large difference can be observed in the prediction
accuracy of DMD asymptotically as $x\to 0$.

Figure 10: Prediction accuracy of DMD on $C^{1}-$ and $C^{2}-$smooth systems
(103). The orange curve shows the accuracy of DMD for $\alpha=0.9$, while the
blue curve corresponds to $\alpha=1$.

Assumption (A3) of Theorem 1 requires a specific nondegeneracy condition to
hold for the observable in order for DMD to give a meaningful approximation.
We now illustrate on an example what happens if this nondegeneracy condition
is not satisfied. Consider the chain of nonlinear oscillators in Section 5.5.
Let us denote the normal modes of the mechanical system (74) linearized around
the stable fixed point $\mathbf{q=0},$ as
$\hat{\mathbf{q}}_{j}(t)=\hat{\mathbf{q}}_{j}(0)e^{\lambda_{j}t}$, where the
mode shapes $\hat{\mathbf{q}}_{j}(0)\in\mathbb{R}^{5}$ satisfy

$\left(\lambda_{j}^{2}\mathbf{M}+\lambda_{j}\mathbf{C}+\mathbf{K}\right)\hat{\mathbf{q}}_{j}(0)=0.$

Therefore, in the phase space spanned by the variable

$\mathbf{x}=\left(q_{1},q_{2},...,q_{5},\dot{q_{1}},\dots,\dot{q}_{5}\right)\in\mathbb{R}^{10},$

the eigenvectors of the linear part of (74)are given as

$\mathbf{e}_{j}=\left(\hat{\mathbf{q}}_{j}(0),\lambda_{j}\hat{\mathbf{q}}_{j}(0)\right)^{T}.$

With our choice of parameters, the slowest pair of eigenvalues are

$\lambda_{1,2}=-0.0012\pm i0.2825,$

and hence there exists a 2D slow SSM tangent to the slow spectral subspace

$E=\text{span}\left([\text{Re}(\mathbf{e}_{1}),\text{Im}(\mathbf{e}_{1})]\right)=\text{span}\left(\left[\hat{\mathbf{q}}_{1}(0),\text{Re}\left(\lambda_{1}\right)\hat{\mathbf{q}}_{1}(0)\right],\left[0,\text{Im}\left(\lambda_{1}\right)\hat{\mathbf{q}}_{1}(0)\right]\right),$

where $\mathbf{e}_{1}$ is the eigenvector corresponding to $\lambda_{1}$.

As the dimension of $\mathcal{W}\left(E\right)$ is $d=2$, we consider three
possible 2D observables with respect to which Theorem 1 could potentially be
applied:

$\boldsymbol{\phi}_{1}(\mathbf{x})=\left(\begin{array}[]{c}\langle\text{Re}(\mathbf{e}_{1}),\mathbf{x}\rangle\\\
\langle\text{Im}(\mathbf{e}_{1}),\mathbf{x}\rangle\end{array}\right)\quad\boldsymbol{\phi}_{2}(\mathbf{x})=\left(\begin{array}[]{c}q_{1}\\\
\dot{q}_{1}\end{array}\right)\quad\boldsymbol{\phi}_{3}(\mathbf{x})=\left(\begin{array}[]{c}q_{1}\\\
q_{2}\end{array}\right).$

The Jacobians of these observables with respect to $x$ can be computed as

$\displaystyle D\boldsymbol{\phi}_{1}(\mathbf{0})$
$\displaystyle=\left(\text{Re}(\mathbf{e}_{1})\,\text{Im}(\mathbf{e}_{1})\right)^{\mathrm{T}}=\left(\begin{array}[]{ccccc}*&*&*&\cdots&*\\\
0&0&0&\cdots&*\end{array}\right),$ $\displaystyle
D\boldsymbol{\phi}_{2}(\mathbf{0})$
$\displaystyle=\left(\begin{array}[]{ccccc}1&0&0&\cdots&0\\\
0&1&0&\cdots&0\end{array}\right),$ $\displaystyle
D\boldsymbol{\phi}_{3}(\mathbf{0})$
$\displaystyle=\left(\begin{array}[]{ccccc}1&0&0&\cdots&0\\\
0&0&1&\cdots&0\end{array}\right),$

where $*$ denotes nonzero entries. Therefore, we require __

$\mathrm{rank\,}\left[D\boldsymbol{\phi}_{1,2,3}\left(\mathbf{0}\right)|_{E}\right]=2.$
(104)

The restriction of the Jacobians of the observables to $E$ is
$D\boldsymbol{\phi}_{i}(\mathbf{0})[\text{Re}(\mathbf{e}_{1}),\text{Im}(\mathbf{e}_{1})]$.
It can be seen from the structure of the matrices, that for
$\boldsymbol{\phi}_{1}$ and $\boldsymbol{\phi}_{2}$, the rank of the
restricted Jacobian is 2, but for $\boldsymbol{\phi}_{3}$, it is only of rank
1. Therefore, the observable vector $\boldsymbol{\phi}_{3}$ cannot be used to
parametrize the SSM and to derive a reduced-order model.

Figures 11a-b show that observing the dynamics via $\boldsymbol{\phi}_{1}(x)$
and $\boldsymbol{\phi}_{2}(x)$ indeed results in a successful embedding of the
dynamics on $\mathcal{W}\left(E\right)$ near the fixed point. At the same
time, Fig. 11c shows that using the observable $\boldsymbol{\phi}_{3}(x)$ does
not correctly reproduce the dynamics near the origin.

Figure 11: The dynamics near the $\mathbf{x}=\mathbf{0}$ point on the slow 2D
SSM $\mathcal{W}\left(E\right)$ of system (74) represented via different
observables (a) Observable $\boldsymbol{\phi}_{1}(\mathbf{x})$ (b) Observable
$\boldsymbol{\phi}_{2}(\mathbf{x})$ (c) Observable
$\boldsymbol{\phi}_{3}(\mathbf{x})$.

Assumption (A4) of Theorem 1 states that the collected data has to be rich
enough to a neighborhood of the SSM near the fixed point. To illustrate this,
we modify the linear system in (67) to have three real eigenvalues by letting
$\boldsymbol{\Lambda}=\text{diag}\left(a,b,c\right)\in\mathbb{R}^{3}$. We then
collect the resulting trajectories in the phase space in the observable data
matrix

$\boldsymbol{\Phi}=\left(\begin{array}[]{cccc}x_{1}(0)&ax_{1}(0)&&a^{n}x_{1}(0)\\\
x_{2}(0)&bx_{2}(0)&\cdots&b^{n}x_{2}(0)\\\
x_{3}(0)&cx_{3}(0)&&c^{n}x_{3}(0)\end{array}\right).$

Assuming $c<b<a<0$, we seek to identify from the data matrix $\Phi$ the 2D
slow eigenspace $\mathcal{W}\left(E\right)=E$ corresponding to the eigenvalues
$a$ and $b$. Assume now that $x_{2}(0)=x_{3}(0)=0$ holds. In this case, the
second and third rows of $\Phi$ are identically zero and hence its rank is
only 1, which violates assumption (A4) of Theorem 1. Indeed, the trajectory
observed in $\boldsymbol{\Phi}$ is restricted to the slowest 1D eigenspace
$\mathcal{W}\left(E_{1}\right)=E_{1}$ corresponding to the eigenvalue $a$ and
hence cannot be used to approximate the shape and internal dynamics of
$\mathcal{W}\left(E\right)$.

In practice, a generic trajectory of a nonlinear system would not be
restricted to a lower-dimensional invariant subspace and hence the data matrix
$\Phi$ would have full rank with probability one. However, if the available
data is close to a lower-dimensional invariant subspace, the condition number
of the data matrix may become prohibitively large. So a near-failure of
assumption (A4) is already detrimental to the quality of the results obtained
from DMD.

##  Appendix D Proof of Theorem 4

For $\mathbf{f}_{2}\in C^{r}$ and under the nonresonance conditions

$\lambda_{k}\neq\sum_{j=1}^{n}m_{j}\lambda_{j},\qquad
j,k=1,\ldots,n,\qquad\sum_{j=1}^{n}m_{j}\geq 2,$ (105)

the $C^{\infty}$ linearization theorems of Sternberg [52] and Poincaré [46]
guarantee the local existence of a near-identity, linearizing change of
coordinates of the form

$\mathbf{x}=\mathbf{y}+\mathbf{h}(\mathbf{y}),\qquad\mathbf{y}\in
U\subset\mathbb{R}^{n},\qquad\mathbf{h}(\mathbf{y}),\mathbf{h}^{-1}(\mathbf{y})=\mathcal{O}\left(\left|\mathbf{y}\right|^{2}\right),$
(106)

with $\mathbf{h}\in C^{r}$. Additionally, a strengthening of Sternberg’s
linearization theorem by Eldering et al. [17]777This is based on the approach
of Lan and Mezic [35] in the setting of the $C^{1}$ linearization theorem by
Hartman [23]. implies that the transformation in eq. (106) is unique and
globally defined on the full domain of attraction of the fixed point at the
origin. As pointed out by Haller et al. [22], this implies the existence of a
unique $d$-dimensional, normally attracting spectral submanifold
$\mathcal{W}(E)\in C^{\infty}$ for system (11) whose internal dynamics govern
the longer-term behavior of all trajectories near the origin.888Cabré et al.
[10] guarantees the existence and uniqueness of $\mathcal{W}(E)\in C^{\infty}$
under weaker nonresonance conditions that allow resonances inside $E$ but no
resonances of the form (105) with
$\lambda_{k}\in\mathrm{Spect}\,\left(A|_{F}\right)$ and
$\lambda_{j}\in\mathrm{Spect}\,\left(\mathbf{A}|_{E}\right).$ We do not rely
on these results here, as the next step in our construction is local
linearization within $\mathcal{W}(E)$, which cannot be carried out in the
presence of resonances inside $E$.

By the main assumption of the theorem, we have
$\mathrm{Re}\text{$\lambda_{1}<0$}$. Under this assumption, the infinitely
many nonresonance conditions in (105) simplify to finitely many. Indeed, if we
define the spectral quotient as the positive integer $Q$ given by

$Q=\mathrm{Int}\left[\frac{\max_{i}\left|\mathrm{Re}\,\lambda_{i}\right|}{\min_{i}\left|\mathrm{Re}\,\lambda_{i}\right|}\right]+1,$

Then, for all $m_{j}\in\mathbb{N}$ with $\sum_{j=1}^{n}m_{j}>Q$, we have the
estimate

$\displaystyle\sum_{j=1}^{n}m_{j}\mathrm{Re}\,\lambda_{j}$
$\displaystyle\leq-\sum_{j=1}^{n}m_{j}\min_{i}\left|\mathrm{Re}\,\lambda_{i}\right|<-\sum_{j=1}^{n}m_{j}\frac{\max_{i}\left|\mathrm{Re}\,\lambda_{i}\right|}{Q}=-\max_{i}\left|\mathrm{Re}\,\lambda_{i}\right|\frac{\sum_{j=1}^{n}m_{j}}{Q}$
$\displaystyle<-\max_{i}\left|\mathrm{Re}\,\lambda_{i}\right|,$

implying

$\sum_{j=1}^{n}m_{j}\mathrm{Re}\,\lambda_{j}<\mathrm{Re}\,\lambda_{k},\qquad
k=1,\ldots,n,\quad\sum_{j=1}^{n}m_{j}>Q.$ (107)

Consequently, in this case, no integer resonance among the eigenvalues of
$\mathbf{A}$ can occur for $\sum_{j=1}^{n}m_{j}>Q$ and hence the nonresonance
condition (105) can be relaxed to assumption (A1) of Theorem 4.

Based on these results, we apply the $C^{r}$ linearizing transformation (106)
followed by a linear change of coordinates to
$\left(\boldsymbol{\xi},\boldsymbol{\eta}\right)$uses in the proof of Theorem
1 (see eq. (90)). As se have seen there, this coordinate change maps the
spectral submanifold $\mathcal{W}(E)$ into the $\boldsymbol{\eta}=\mathbf{0}$
plane. On $\mathcal{W}(E)$, the reduced dynamics obeys the linear ODE

$\dot{\boldsymbol{\xi}}=\boldsymbol{\Lambda}_{E}\boldsymbol{\xi},\qquad\boldsymbol{\Lambda}_{E}=\left(\mathbf{T}^{-1}\mathbf{A}\mathbf{T}\right)|_{E},$
(108)

where $\boldsymbol{\Lambda}_{E}$, the real Jordan canonical form of
$\mathbf{A}|_{E},$ has eigenvalues $\lambda_{1},\ldots,\lambda_{d}$.

Restricted to $\mathcal{W}(E),$ we can express the $d$-dimensional observable
vector $\boldsymbol{\phi}(\mathbf{x})$ as

$\boldsymbol{\varphi}(\boldsymbol{\xi})=\boldsymbol{\phi}\left(\mathbf{T}\left(\begin{array}[]{c}\boldsymbol{\xi}\\\
\mathbf{0}\end{array}\right)+\mathbf{h}\left(\begin{array}[]{c}\boldsymbol{\xi}\\\
\mathbf{0}\end{array}\right)\right)\in\mathbb{R}^{d}.$ (109)

On the manifold $\mathcal{W}(E)$, we can use $\boldsymbol{\varphi}$ locally
near the origin as a new variable if the coordinate change from
$\boldsymbol{\xi}$ to $\boldsymbol{\varphi}$ is a $C^{r}$ diffeomorphism on
$\mathcal{W}(E)$, i.e., $D_{\boldsymbol{\xi}}\boldsymbol{\varphi}$ is
nonsingular at the origin. To see if this is the case, note that

$D_{\boldsymbol{\xi}}\boldsymbol{\psi}(\boldsymbol{\xi})=D\boldsymbol{\phi}\left(\mathbf{T}\left(\begin{array}[]{c}\boldsymbol{\xi}\\\
\mathbf{0}\end{array}\right)\right)\mathbf{T}\left(\begin{array}[]{c}\mathbf{I}_{d\times
d}\\\ \mathbf{0}\end{array}\right),$

which gives

$D_{\boldsymbol{\xi}}\boldsymbol{\gamma}(\mathbf{0})=D\boldsymbol{\phi}\left(\mathbf{0}\right)\mathbf{T}\left(\begin{array}[]{c}\mathbf{I}_{d\times
d}\\\
\mathbf{0}\end{array}\right)=D\boldsymbol{\phi}\left(\mathbf{0}\right)\mathbf{T}_{E},$

with the operator $\mathbf{T}_{E}$ defined in (17). By assumption (37) of the
theorem, we have

$\det\,\left[D\boldsymbol{\phi}\left(\mathbf{0}\right)\mathbf{T}_{E}\right]\neq
0,$ (110)

and hence, by the inverse function theorem, $\boldsymbol{\xi}$ can be
expressed from the relationship (109) near the origin as a smooth function of
$\boldsymbol{\varphi}$. Specifically, there exists a $C^{r}$ diffeomorphism
$\boldsymbol{\psi}\colon\mathbb{R}^{d}\to\mathbb{R}^{d}$ such that

$\boldsymbol{\xi}=\boldsymbol{\psi}(\boldsymbol{\varphi}),\qquad\boldsymbol{\psi}(\boldsymbol{0})=\boldsymbol{0},\quad
D_{\boldsymbol{\varphi}}\boldsymbol{\psi}(\boldsymbol{\varphi})=\left[D_{\boldsymbol{\xi}}\boldsymbol{\varphi}(\boldsymbol{\xi})\right]^{-1}.$
(111)

Differentiating the first equation in (111), substituting eq. (108), using the
relations (111) and Taylor expanding the result around
$\boldsymbol{\varphi}=\mathbf{0}$, we obtain the reduced dynamics on
$\mathcal{W}(E)$ in the form

$\displaystyle\dot{\boldsymbol{\varphi}}$
$\displaystyle=\left[D_{\boldsymbol{\varphi}}\boldsymbol{\psi}(\boldsymbol{\varphi})\right]^{-1}\boldsymbol{\Lambda}_{E}\boldsymbol{\psi}(\boldsymbol{\varphi})$
$\displaystyle=D_{\boldsymbol{\xi}}\boldsymbol{\varphi}(\mathbf{0})\boldsymbol{\Lambda}_{E}\left[D_{\boldsymbol{\xi}}\boldsymbol{\varphi}(\mathbf{0})\right]^{-1}\boldsymbol{\varphi}+\mathbf{q}\left(\boldsymbol{\varphi}\right),\qquad\mathbf{q}\left(\boldsymbol{\varphi}\right)=o\left(\left|\boldsymbol{\varphi}\right|\right)$
(112)

for some $C^{r}$ function $\mathbf{q}\colon
U\subset\mathbb{R}^{d}\to\mathbb{R}^{d}$ defined in a neighborhood of
$\boldsymbol{\varphi}=\mathbf{0}$. This proves statement (i) of the theorem.

We have seen that $D_{\boldsymbol{\xi}}\boldsymbol{\varphi}(\mathbf{0})$ is
nonsingular and hence the linear part of the nonlinear ODE in (112) has the
same eigenvalues as $\boldsymbol{\Lambda}_{E}$ does.
$\left(\mathbf{T}^{-1}\mathbf{A}\mathbf{T}\right)|_{E}$. Consequently, the
origin of system (113) is a hyperbolic fixed point. If the nonresonance
conditions (36) hold for the full spectrum of $\mathbf{A}$, then they also
hold for any subset of this spectrum, i.e., for the eigenvalues of
$\boldsymbol{\Lambda}_{E}$ as well. Consequently, from the linearization
theorems of Sternberg [52] and Eldering et al. [17], we can also deduce a
similar near-identity, $C^{r}$ linearizing diffeomorphism

$\boldsymbol{\varphi}=\boldsymbol{\gamma}+\boldsymbol{\ell}(\boldsymbol{\gamma}),\qquad\boldsymbol{\ell}(\boldsymbol{\gamma})=o\left(\left|\boldsymbol{\gamma}\right|\right),\quad\boldsymbol{\ell}\in
C^{r},$

for _$r\in\mathbb{N}^{+}\cup\left\\{\infty,a\right\\}$_ under the nonresonance
assumption (36). In the new $\boldsymbol{\gamma}$ coordinates, the dynamics on
the SSM $\mathcal{W}\left(E\right)$ is governed by the linear ODE

$\dot{\boldsymbol{\gamma}}=D_{\boldsymbol{\xi}}\boldsymbol{\gamma}(\mathbf{0})\boldsymbol{\Lambda}_{E}\left[D_{\boldsymbol{\xi}}\boldsymbol{\gamma}(\mathbf{0})\right]^{-1}\boldsymbol{\gamma},$
(113)

which proves statement (ii) of the theorem.

The partial differential equation (41) can then be obtained by differentiating
eq. (39) and substituting eqs. (38) and (40) into the resulting equation. The
expansions (42) and (43) then follow from the guaranteed smoothness class of
the linearizing transformation for general $r$ and for $r=a$, respectively.
This concludes the proof of statement (iii) of the theorem.

##  Appendix E Reduced dynamics on periodically forced SSMs

Let us denote the time-$T$ map for the periodically forced system by
$\mathbf{P}_{\varepsilon}^{T}$. Since $\mathbf{P}_{0}^{T}$ represents the
period-$T$ sampling of the flow map of the unperturbed system, it has a fixed
point at $\mathbf{x}=\mathbf{0}$. By the implicit function theorem,
$\mathbf{P}_{\varepsilon}^{T}$ also has a fixed point $O(\varepsilon)-$close
to $\mathbf{x}=\mathbf{0}$. This corresponds to a periodic orbit of (60) close
to the origin. As explained by Haller et al. [22], the $C^{\infty}$
linearization results of Sternberg [52] and Poincaré [46] also establish the
smoothness of the linearizations in $\varepsilon.$ Applied to the discrete
dynamical system defined by the Poincaré-map, these results guarantee the
existence of a linearizing transformation of the form

$\mathbf{x}=\mathbf{y}+\mathbf{h}^{\varepsilon}(\mathbf{y}),\qquad\mathbf{y}\in
U\subset\mathbb{R}^{n},\qquad\mathbf{h}(\mathbf{y}),\mathbf{h}^{-1}(\mathbf{y})=\mathcal{O}\left(\left|\mathbf{y}\right|^{2},\varepsilon|\mathbf{y}|^{2}\right),$

with $\mathbf{h}\in C^{r}$ in both $\mathbf{y}$ and $\varepsilon$.

We now apply the $\varepsilon$-dependent linearizing transformation followed
by the linear change of coordinates to the coordinates
$\left(\boldsymbol{\xi},\boldsymbol{\eta}\right)$ introduced in the proof of
Theorem 4. Specifically, let

$\mathbf{x}=\mathbf{Ty}+\mathbf{h}^{\varepsilon}(\mathbf{Ty}),\quad\mathbf{y}=\left(\boldsymbol{\xi},\boldsymbol{\eta}\right)^{\mathrm{T}}\in\mathbb{R}^{d}\times\mathbb{R}^{n-d},$
(114)

with the matrix $\mathbf{T}$ defined in formula (17), i.e., containing
eigenvalues of $\mathbf{A}$. As the change of coordinates (114) is smooth in
$\varepsilon$, it can be written as
$\mathbf{x}=\mathbf{T}\mathbf{y}+\mathbf{h}^{0}(\mathbf{T}\mathbf{y})+O(|y|^{2}\varepsilon).$
Differentiation of (114) with respect to time yields

$\dot{\mathbf{x}}=\left(\mathbf{T}+D\mathbf{h}^{\varepsilon}(\mathbf{T}\mathbf{y})\mathbf{T}\right)\dot{\mathbf{y}}=\mathbf{A}(\mathbf{Ty}+\mathbf{h}^{\varepsilon}(\mathbf{Ty}))+\tilde{\mathbf{f}}(\mathbf{Ty}+\mathbf{h}^{\varepsilon}(\mathbf{Ty}))+\varepsilon\mathbf{F}(\mathbf{x}+h^{\varepsilon}(\mathbf{Ty}),t),$

which can be rearranged as

$\displaystyle\dot{\mathbf{y}}$
$\displaystyle=\left(\mathbf{T}+D\mathbf{h}^{0}(\mathbf{Ty})\mathbf{T}+O(|\mathbf{y}|\varepsilon)\right){}^{-1}\left[\mathbf{A}\left(\mathbf{T}\mathbf{y}+\mathbf{h}^{0}(\mathbf{Ty})+O(|\mathbf{y}|^{2}\varepsilon)\right)\right]$
(115)
$\displaystyle+\left(\mathbf{T}+D\mathbf{h}^{0}(\mathbf{Ty})\mathbf{T}+O(|\mathbf{y}|\varepsilon)\right){}^{-1}\left[\tilde{\mathbf{f}}(\mathbf{Ty}+\mathbf{h}^{0}(\mathbf{Ty})+O(|\mathbf{y}|^{2}\varepsilon))+\varepsilon\mathbf{F}(\mathbf{Ty}+\mathbf{h}^{0}(\mathbf{Ty}),t)\right].$

Since $y+h^{0}(y)$ linearizes (60) for $\varepsilon=0$, (115) simplifies to

$\dot{\mathbf{y}}=\boldsymbol{\Lambda}\mathbf{y}+\left(\mathbf{T}+D\mathbf{h}^{0}(\mathbf{Ty})\mathbf{T}\right){}^{-1}\varepsilon\mathbf{F}(\mathbf{0},t)+O(|\mathbf{y}|^{2}\varepsilon,\mathbf{y}\varepsilon^{2}),\quad\mathbf{y}=(\boldsymbol{\xi},\boldsymbol{\eta})^{\mathrm{T}}.$
(116)

The smoothest SSM $\mathcal{W_{\varepsilon}}(E,t)$ of the periodic orbit is
$O(\varepsilon)-$close to that of the fixed point and can be written as

$\mathcal{W_{\varepsilon}}(E,t)=\left(\begin{array}[]{c}\boldsymbol{\xi}\\\
O(\varepsilon)\end{array}\right).$

We now neglect the small, time-dependent correction to the
$\boldsymbol{\eta}$-component of $\mathcal{W_{\varepsilon}}(E,t).$ The reduced
dynamics can then be written as the projection of (116) to the
$\boldsymbol{\eta}=0$ plane, which becomes

$\dot{\boldsymbol{\xi}}=\boldsymbol{\Lambda}_{E}\boldsymbol{\xi}+\left(\mathbf{T}+D\mathbf{h}^{0}(\mathbf{Ty})\mathbf{T}\right){}^{-1}|_{E}\varepsilon\mathbf{F}(0,t)+O(|\mathbf{y}|^{2}\varepsilon,y\varepsilon^{2}),\qquad\boldsymbol{\Lambda}_{E}=\left(\mathbf{T}^{-1}\mathbf{AT}\right)|_{E}.$
(117)

The restriction of the observable vector $\boldsymbol{\phi}$ to the subspace
$E$ can be written as $\boldsymbol{\varphi}(\boldsymbol{\xi}).$ Repeating the
arguments of D, the change of coordinates
$\boldsymbol{\xi}=\boldsymbol{\psi}(\boldsymbol{\varphi})$ can be carried out,
which transforms (117) to

$\displaystyle\dot{\boldsymbol{\varphi}}$
$\displaystyle=D_{\boldsymbol{\xi}}\boldsymbol{\varphi}(\mathbf{0})\boldsymbol{\Lambda}_{E}\left[D_{\boldsymbol{\xi}}\boldsymbol{\varphi}(\mathbf{0})\right]^{-1}\boldsymbol{\varphi}+\mathbf{q}\left(\boldsymbol{\varphi}\right)+\varepsilon\left\\{\left[\mathbf{I}+D\mathbf{h}^{0}\left(D_{\boldsymbol{\xi}}\boldsymbol{\varphi}(\mathbf{0})\right)^{-1}\right]\mathbf{T}_{E}\right\\}^{-1}\mathbf{F}(\mathbf{0},t)+O(|\boldsymbol{\varphi}|^{2}\varepsilon,\left|\boldsymbol{\varphi}\right|\varepsilon^{2}),$
(118) $\displaystyle\mathbf{q}\left(\boldsymbol{\varphi}\right)$
$\displaystyle=o\left(\left|\boldsymbol{\varphi}\right|\right),$

which proves formula (61).

##  Appendix F Forced response from approximate DDL

As stated in Section 4.3.4, the linearizing change of coordinates
$\boldsymbol{\varphi}=\boldsymbol{\kappa}\left(\boldsymbol{\gamma}\right)=\boldsymbol{\gamma}+\boldsymbol{\ell}(\boldsymbol{\gamma})$
transforms the reduced dynamics of the forced system (118) to the form
$\dot{\boldsymbol{\gamma}}=\mathbf{B}\boldsymbol{\gamma}+\varepsilon(\mathbf{I}+D\boldsymbol{\ell}(\boldsymbol{\gamma}))^{-1}\hat{\mathbf{F}}(\mathbf{0},t).$
Therefore, a trajectory of the forced system may be computed by integrating
(63) and transforming to the reduced coordinates $\boldsymbol{\varphi}$ as

$\boldsymbol{\varphi}(t)=\boldsymbol{\gamma}\left(t\right)+\boldsymbol{\ell}\left(\boldsymbol{\gamma}(t)\right).$

An approximate version of (63) may be obtained by assuming (64), i.e.,
$\varepsilon(\mathbf{I}+D\boldsymbol{\ell}(\boldsymbol{\gamma}))^{-1}=\varepsilon\mathbf{I}+O\left(\varepsilon|\boldsymbol{\gamma}|\right)\approx\varepsilon\mathbf{I}.$
This results in the forced linear system

$\dot{\boldsymbol{\gamma}}=\mathbf{B}\boldsymbol{\gamma}+\varepsilon\hat{\mathbf{F}}\left(\mathbf{0},t\right).$
(119)

However, we must keep in mind that the $O\left(\varepsilon|\gamma|\right)$
terms are already retained when we compute
$\boldsymbol{\varphi}(t)=\boldsymbol{\gamma}(t)+\boldsymbol{\ell}\left(\boldsymbol{\gamma}(t)\right)$,
therefore terms of the same order must also be retained in the dynamics (63).
Nevertheless, assumption (64) makes the dynamics (119) an inhomogeneous linear
system of ODEs, which is simple to solve. We refer to DDL under assumption
(64) as Approximate DDL. We compare forced response predictions of the
analytic linearization, DDL, Approximate DDL and DMD on the forced-damped
Duffing system discussed in Section 5.3.

In Fig. 12, we repeat the results already presented in Fig. 6, which should
now be compared to the forced response predicted using Approximate DDL. Panel
(a) shows that although Approximate DDL is accurate only for low amplitudes,
it clearly outperforms the DMD-model. For higher amplitudes, predictions are
not accurate, since the forced linear system (119), by construction, cannot
show nonlinear behavior. In contrast, the analytic linearization and the full
DDL continue to yield accurate predictions for the forced response, even in
the nonlinear regime.

Figure 12: Periodic response of the forced Duffing oscillator (71). Panel (a)
shows DDL-predictions under the simplifying assumption (64). Panels (b)-(d)
are the same as Fig. 6.


# SOS Construction of Compatible Control Lyapunov and Barrier Functions

Michael Schneeberger  Florian Dörfler  Silvia Mastellone  Institute of
Electrical Engineering, FHNW, Windisch, Switzerland (e-mails:
<EMAIL_ADDRESS><EMAIL_ADDRESS> Department of
Information Technology and Electrical Engineering, ETH Zürich, Zürich,
Switzerland, (e-mail<EMAIL_ADDRESS>
###### Abstract

We propose a novel approach to certify closed-loop stability and safety of a
constrained polynomial system based on the combination of Control Lyapunov
Functions (CLFs) and Control Barrier Functions (CBFs). For polynomial systems
that are affine in the control input, both classes of functions can be
constructed via Sum Of Squares (SOS) programming. Using two versions of the
Positivstellensatz we derive an SOS formulation seeking a rational controller
that — if feasible — results in compatible CLF and multiple CBFs.

##  1 Introduction

When dealing with systems that have state constraints, it is crucial to have a
controller that ensures both stability and compliance with the constraints. In
most cases, feedback control design focuses mainly on achieving stability,
while protection functions only engage when constraints are violated.
Unfortunately, this approach results in downtime and the need to investigate
faults. By contrast, a controller that is both stable and safe can trade off
control performance for the ability to prevent unsafe states.

The focus of this paper is on constructing compatible CLFs and CBFs that can
certify both stability and safety in control systems. A CLF establishes
conditions for the existence of a stabilizing controller for a given control
system. Similarly, a CBF guarantees the existence of a controller that can
render the control system safe. According to Wieland and Allgöwer (2007), a
system is considered safe if any state trajectory starting from a safe set of
states remains within an allowable region defined by the state constraints.
For systems that are affine in the input, a controller that meets the CLF and
CBF conditions can be implemented by solving an online Quadratic Programming
(QP), see Ames et al. (2019). However, for some states, these conditions may
conflict with each other. For both to hold jointly for all states, the
control-sharing property is additionally required, see Grammatico et al.
(2013); Xu (2016). Finding such compatible CLF and CBF is generally difficult.
For polynomial systems, however, this can be achieved by formulating SOS
constraints and solving them via Semidefinite Programming (SDP).

The contributions of this paper is twofold: first, we derive SOS constraints
on compatible CLF and multiple CBFs using two versions of the
Positivstellensatz, and second, an algorithm is developed that efficiently
finds solutions to these SOS constraints by maximizing a surrogate of the
volume of the safe set. The conditions on compatible CBFs are formulated in
Isaly et al. (2022); Tan and Dimarogonas (2022) but without giving a method to
construct them. A constructive approach described by Clark (2021) is based on
the introduction of additional SOS constraints to enforce compatibility. In
this paper, we reveal a correspondence between the SOS constraint derived from
the CLF, resp. CBF, condition, and the existence of a rational controller that
renders the closed-loop system stable, resp. safe. By restricting such
controllers to be identical, we derive a new set of SOS constraints that
guarantee compatibility between a CLF and multiple CBFs without the
introduction of additional SOS constraints.

These SOS constraints contain bilinear terms, and hence, cannot be directly
converted to an SDP. We therefore present an alternating algorithm that
searches simultaneously for a CLF and multiple CBFs by repeatedly solving two
SDPs. In particular, the algorithm seeks to maximize the volume of the safe
set, which is given by the intersection of the invariant sets defined by each
CBF. Multiple CBFs offer additional flexibility to increase the volume of the
safe set when a single CBF does not suffice. Similar approaches were presented
in Anghel et al. (2013); Kundu et al. (2019); Wang et al. (2022) but they
either only searched for a single CLF or a single CBF. Korda et al. (2014)
proposed another intriguing approach for identifying a safe set using an
infinite-dimensional linear programming problem. We demonstrate the utility of
our approach with a power converter control example for which safety is of
paramount importance.

The paper is structured as follows: In Section 2, we introduce the notation
adopted in the paper and review some preliminaries. Then, we recall the
definition of a CLF and CBF and derive a rational controller in Section 3. In
Section 4, we combine the CLF and CBF by introducing the control-sharing
property. Section 5 defines the SOS program encoding the CLF and CBF
conditions, and Section 6 presents the algorithm that solves the SOS program.
Numerical simulation are given in Section 7. Finally, Section 8 is dedicated
to concluding remarks and future work.

##  2 Preliminaries & Notation

###  2.1 Notation

Across this paper, we adopt the following notation. The shorthand
$[t]:=\\{1,2,...,t\\}$ is used to denote a range of numbers. A scalar function
$V:\mathbb{R}^{n}\to\mathbb{R}$ is positive definite w.r.t. $x^{*}$, if
$V(x^{*})=0$ and $V(x)>0$ for $x\neq x^{*}$. $R[x]$ refers to the set of
scalar polynomials in variables $x\in\mathbb{R}^{n}$, and $\Sigma[x]$ refers
to the set of scalar SOS polynomials in $x$. A polynomial $p(x)$ is an SOS
polynomial if it can be written as $p(x)=\sum_{i=1}^{k}g_{i}(x)^{2}$ for
$g_{i}(x)\in R[x]$ and $i\in[k]$. If $p(x)\in\Sigma[x]$, it can be expanded to

$\displaystyle p(x)=Z(x)^{T}QZ(x),$

where $Z(x)$ is a vector of monomials in $x$, and $Q$ is a square positive
semidefinite matrix. In the following, we consider a polynomial control system
described as

$\displaystyle\dot{x}=\mathcal{F}_{c}(x,u)=f(x)+G(x)u_{c},$ (1)

where $f(x)\in\left(R[x]\right)^{n}$ and $G(x)\in\left(R[x]\right)^{n\times
m}$ are polynomial matrices, and $u_{c}\in\mathbb{R}^{m}$ is the control input
vector. System (1) with polynomial state feedback control policy
$u(x)\in\left(R[x]\right)^{m}$ results in a closed-loop polynomial system of
the form

$\displaystyle\dot{x}=\mathcal{F}_{a}(x)=f(x)+G(x)u(x)$ (2)

The state $x^{*}$ is called an equilibrium of (2), if
$\mathcal{F}_{a}(x^{*})=0$.

A set $\mathcal{X}$ is called forward invariant Khalil (2002) with respect to
(2) if for every $x_{0}\in\mathcal{X}$, $x(t)\in\mathcal{X}$ for all
$t\in\mathbb{R}_{+}$. We say that a system (2) is safe Wieland and Allgöwer
(2007) w.r.t. an allowable set of states
$\mathcal{X}_{a}\subseteq\mathbb{R}^{n}$ and the safe states
$\mathcal{X}_{s}\subseteq\mathbb{R}^{n}$, if $\mathcal{X}_{s}$ is forward
invariant and $\mathcal{X}_{s}\subseteq\mathcal{X}_{a}$.

###  2.2 SOS Programming

An SOS program minimizing a quadratic cost function subject to SOS constraints
is defined as follows:

$\begin{array}[]{lll}\underset{v}{\mbox{minimize}}&v^{T}Pv+c^{T}v\\\
\mbox{subject to}&p_{i}(x,v)\in\Sigma[x]&i=1,...,N.\end{array}$

Here $P\in\mathbb{R}^{l\times l}$ is a positive semi-definite matrix,
$c\in\mathbb{R}^{l}$ is a vector, and $p_{i}(x,v)\in\Sigma[x]$ is an SOS
polynomial in $x\in\mathbb{R}^{n}$ parameterized by $v\in\mathbb{R}^{l}$, i.e.
it can be expanded as

$\displaystyle p_{i}(x,v)=Z_{i}(x)^{T}Q_{i}(v)Z_{i}(x),$ (3)

where $Q_{i}(v)=Q_{0,i}+v_{1}Q_{1,i}+...+v_{l}Q_{1,l}$ is a square positive
semidefinite matrix linearly parameterized by $v$ encoding the coefficients of
the polynomial.

When deriving the SOS constraints in the following subsections, the parameter
$v$ is omitted for simplicity, and the polynomial is simply written as
$p_{i}(x)$ instead.

###  2.3 Positivstellensätze

In the following, we present a version of the Positivstellensatz that can be
regarded as a specialization of the weak Positivstellensatz. This version will
be relevant for our analysis later on.

###### Theorem 1

Given polynomials $f_{1}(x)$, $f_{2}(x)$, $g_{1}(x)$, $h_{1}(x)$, …,
$h_{m+1}(x)\in R[x]$ such that

$\displaystyle\begin{split}\Bigl{\\{}x\in\mathbb{R}^{n}\mid&f_{1}(x)\geq
0,f_{2}(x)\geq 0,g_{1}(x)\neq 0,\\\
&h_{1}(x)=0,...,h_{m+1}(x)=0\Bigr{\\}}=\emptyset,\end{split}$ (4)

then there exist polynomials

  * •

$\begin{aligned} f_{\text{cone}}(x)=s_{1}(x)f_{1}(x)&+s_{2}(x)f_{2}(x)\\\
&+s_{3}(x)f_{1}(x)f_{2}(x)\end{aligned}$

  * •

$h_{\text{ideal}}(x)=p_{1}(x)h_{1}(x)+...+p_{m+1}(x)h_{m+1}(x)$

such that

$\displaystyle-
f_{\text{cone}}(x)-h_{\text{ideal}}(x)-g_{1}(x)^{2k}\in\Sigma[x]$ (5)

for all $x\in\mathbb{R}^{n}$, where $p_{1}(x),...,p_{m+1}(x)\in R[x]$ are
polynomials, $s_{1}(x),s_{2}(x),s_{3}(x)\in\Sigma[x]$ are SOS polynomials, and
$k\in\mathbb{N}_{+}$.

Using (Bochnak et al., 2013, Theorem 4.4.2), we note that the cone $P$
generated by $f_{1}(x)$ and $f_{2}(x)$ is contained in
$\Bigl{\\{}s_{0}(x)+s_{1}(x)f_{1}(x)+s_{2}(x)f_{2}(x)+s_{3}(x)f_{1}(x)f_{2}(x)\mid
s_{0},s_{1},s_{2},s_{3}\in\Sigma[x]\Bigr{\\}}$. $\blacksquare$

If the set (4) is empty, the Positivstellensatz ensures that the polynomials
$s_{1}(x),s_{2}(x),...,p_{m+1}(x)$ exist without specifying their degree. To
find these polynomials computationally will require to iteratively increase
the degree of the polynomials until a solution can be found for (5). An
increase in the degree of the polynomials will, however, deteriorate
computation time. For many practical examples of interest, however, low-degree
polynomials suffice.

The following theorem states a version of the Positivstellensatz that can be
seen as specializations of Putinar’s Positivstellensatz.

###### Theorem 2

Given polynomials $f_{1}(x)$, $f_{2}(x)$, $h_{1}(x)$, …, $h_{m+1}(x)\in R[x]$
such that the set $\Bigl{\\{}x\in\mathbb{R}^{n}\mid f_{2}(x)\geq 0\Bigr{\\}}$
is compact and

$\displaystyle\begin{split}\Bigl{\\{}x\in\mathbb{R}^{n}\mid&f_{1}(x)\geq
0,f_{2}(x)\geq 0,\\\
&h_{1}(x)=0,...,h_{m+1}(x)=0\Bigr{\\}}=\emptyset,\end{split}$ (6)

then there exist polynomials

  * •

$f_{\text{cone}}(x)=f_{1}(x)+s_{2}(x)f_{2}(x)$

  * •

$h_{\text{ideal}}(x)=p_{1}(x)h_{1}(x)+...+p_{m+1}(x)h_{m+1}(x)$

such that

$\displaystyle-f_{\text{cone}}(x)-h_{\text{ideal}}(x)\in\Sigma[x]$ (7)

for all $x\in\mathbb{R}^{n}$, where $p_{1}(x),...,p_{m+1}(x)\in R[x]$ are
polynomials and $s_{2}(x)\in\Sigma[x]$ is an SOS polynomial.

Let’s consider the closed semialgebraic set

$\displaystyle\begin{split}K:=\bigl{\\{}x\in\mathbb{R}^{n}\mid&f_{2}(x)\geq
0,\\\ &h_{k}(x)\geq 0\quad k=1,...,m+1,\\\ -&h_{k}(x)\geq 0\quad
k=1,...,m+1\bigr{\\}},\end{split}$

and the quadratic module

$\displaystyle M(f_{2}(x),h_{1}(x),-h_{1}(x),...,-h_{m+1}(x)),$

then $M$ is Archimedean according to (Laurent, 2009, Theorem 3.17) using the
fact that $\bigl{\\{}x\in\mathbb{R}^{n}\mid f_{2}(x)\geq 0\bigr{\\}}$ is
compact. The empty set condition (6) is equivalent to the condition that
$-f_{1}(x)>0$ on K, and, therefore, $K$ is equivalent to

$\displaystyle\begin{split}-f_{1}(x)=&s_{0}(x)+s_{2}(x)f_{2}(x)\\\
&+\sum_{k=1}^{m+1}\left(s_{h,k}^{p}(x)-s_{h,k}^{n}(x)\right)h_{k}(x),\end{split}$
(8)

where $s_{0}(x),s_{2}(x),s_{h,1}^{p}(x),...,s_{h,m+1}^{n}(x)\in\Sigma[x]$
(Laurent, 2009, Theorem 3.20). We note that $s_{h,k}^{p}(x)-s_{h,k}^{n}(x)$
can be replaced by a polynomial $p_{k}(x)$.

$\blacksquare$

##  3 Closed-loop Stability and Safety

In this section, we first review the concept of CLFs and CBFs. Given a
polynomial system (1), we then propose an approach to construct a rational
controller resulting from the SOS formulation of CLF and CBF conditions. We
prove global asymptotic stability and safety of the resulting closed-loop
system (2) given such a controller.

###  3.1 Control Lyapunov Function

For a control system (1) to be stabilizable around the equilibrium point
$x^{*}$, we require the existence of a control input $u_{c}\in\mathbb{R}^{m}$
at every state $x\in\mathbb{R}^{n}$ that renders the sublevel sets of a scalar
polynomial $V(x)\in R[x]$ forward invariant. This motivates the following
definition of a CLF:

###### Definition 1

(Isidori (1995)) Consider a differentiable function
$V:\mathbb{R}^{n}\to\mathbb{R}$ such that $V(x^{*})=0$ and $V(x)>0$ for all
$x\neq x^{*}$. Such a scalar function is called a Control Lyapunov Function
(CLF) for the control system (1) if

$\displaystyle\nabla V(x)^{T}f(x)<0$ (9)

for all $x\in\Bigl{\\{}x\in\mathbb{R}^{n}\mid\nabla V(x)^{T}G(x)=0,x\neq
x^{*}\Bigr{\\}}$.

The inequality (9) can also be formulated as the empty set condition

$\displaystyle\begin{split}\Bigl{\\{}x\in\mathbb{R}^{n}\mid&\nabla
V(x)^{T}f(x)\geq 0,\nabla V(x)^{T}G(x)=0,\\\ &x\neq
x^{*}\Bigr{\\}}=\emptyset.\end{split}$ (10)

Condition (10) when restricted to _polynomial_ CLFs can be solved via SOS
programming. Hence, we restrict our attention to polynomial scalar functions
$V\in R[x]$ for the rest of the paper.

Next, we replace the inequality constraints
$x=(x_{1},...,x_{n})\allowbreak\neq x^{*}$ in (10) by a single inequality
constraint $l(x)\neq 0$ (c.f. Tan and Packard (2004)), where $l(x^{*})=0$ and
$l(x)\neq 0$ elsewhere. A single inequality constraint has the advantage that
it translates into a simpler SOS constraint. The resulting empty set condition
equivalent to (10) is then given by:

$\displaystyle\begin{split}\Bigl{\\{}x\in\mathbb{R}^{n}\mid&\nabla
V(x)^{T}f(x)\geq 0,\nabla V(x)^{T}G(x)=0,\\\ &l(x)\neq
0\Bigr{\\}}=\emptyset.\end{split}$ (11)

According to Theorem 1, by choosing $k=1$, the empty set condition (11)
becomes

$\displaystyle\begin{split}-s_{1}(x)\nabla V(x)^{T}f(x)-\nabla
V(x)^{T}G(x)p(x)-l(x)^{2}\in\Sigma[x],\end{split}$ (12)

where $s_{1}(x)\in\Sigma[x]$ is an SOS polynomial, and
$p(x)=\begin{bmatrix}p_{1}(x)&...&p_{m}(x)\end{bmatrix}^{T}\in(R[x])^{m}$ is a
vector of polynomials.

If $s_{1}(x)$ is strictly positive w.r.t. $x^{*}$, the following controller —
in the form of a rational function — naturally results from the SOS constraint
(12):

$\displaystyle u_{CLF}(x):=p(x)/s_{1}(x).$ (13)

###### Remark 1

Fixing the controller in (1) to $u_{c}=u_{CLF}(x)$ results in the closed-loop
system (2). Hence, stability is asserted by the existence of a Lyapunov
function $V(x)$. For the sake of readability, however, we keep denoting $V(x)$
a CLF.

###### Lemma 1

Given a CLF $V(x)$ for the control system (1) with $s_{1}(x)$ in (12) strictly
positive, the closed-loop system (2) using controller $u(x)=u_{CLF}(x)$
defined in (13) is globally asymptotically stable (GAS).

The polynomial CLF $V(x)$ is strictly positive definite w.r.t. $x^{*}$. Hence,
$V(x)$ is radially unbounded. From (12), we derive $s_{1}(x)\nabla
V(x)^{T}f(x)+\nabla V(x)^{T}G(x)p(x)<0$ for all $x\neq x^{*}$. Since
$s_{1}(x)$ is strictly positive w.r.t. $x^{*}$, the division by $s_{1}(x)$
results in $\nabla V(x)^{T}\left(f(x)+G(x)u_{CLF}(x)\right)\allowbreak<0$ for
all $x\neq x^{*}$. GAS follows from Theorem 4.2 in Khalil (2002).
$\blacksquare$

The strict positivity condition on $s_{1}(x)$ can be enforced by the SOS
constraint

$\displaystyle s_{1}(x)-\epsilon_{s1}\in\Sigma[x],$ (14)

for some $\epsilon_{s1}>0$.

###  3.2 Control Barrier Function

Similar to the stability argument, safety can be asserted with the existence
of a scalar function $B_{1}:\mathbb{R}^{n}\to\mathbb{R}$. Specifically, the
control system (1) is safe w.r.t. the set of safe states

$\displaystyle\mathcal{X}_{s,1}=\Bigl{\\{}x\in\mathbb{R}^{n}\mid B_{1}(x)\leq
0\Bigr{\\}},$ (15)

if there exists a control input $u_{c}\in\mathbb{R}^{m}$ for every state
$x\in\mathbb{R}^{n}$ such that $\mathcal{X}_{s,1}$ is forward invariant. This
motivates the following definition of a CBF (cf. Wang et al. (2022)):

###### Definition 3.3

Consider a differentiable function $B_{1}:\mathbb{R}^{n}\to\mathbb{R}$ such
that $\mathcal{X}_{s,1}$ is non-empty. Such a scalar function is called a
Control Barrier Function (CBF) for the control system (1) if

$\displaystyle\nabla B_{1}(x)^{T}f(x)<0$ (16)

for all $x\in\Bigl{\\{}x\in\mathbb{R}^{n}\mid\nabla
B_{1}(x)^{T}G(x)=0,B_{1}(x)=0\Bigr{\\}}$.

An alternative definition of a CBF (cf. Ames et al. (2019)) involves a
supremum and a class K function. This definition, however, is not well suited
for polynomial optimization since the resulting functions cannot be directly
translated to polynomial inequalities. Similarly to CLF, we restrict our focus
on _polynomial_ CBFs $B_{1}(x)\in R[x]$ for the rest of the paper.

Inequality (16) can also be formulated as the empty set condition

$\displaystyle\begin{split}\Bigl{\\{}x\in\mathbb{R}^{n}\mid&\nabla
B_{1}(x)^{T}f(x)\geq 0,\nabla B_{1}(x)^{T}G(x)=0,\\\
&B_{1}(x)=0\Bigr{\\}}=\emptyset.\end{split}$ (17)

###### Assumption 1

The set $\bigl{\\{}x\in\mathbb{R}\mid B_{1}(x)=0\bigr{\\}}$ is compact.

According to Theorem 2 and Assumption 1, the empty set condition (17) is
equivalent to the SOS constraint

$\displaystyle\begin{split}&-s_{1}(x)\nabla B_{1}(x)^{T}f(x)-\nabla
B_{1}(x)^{T}G(x)p(x)\\\ &\qquad-p_{m+1}(x)B_{1}(x)\in\Sigma[x],\end{split}$
(18)

where $s_{1}(x)\in\Sigma[x]$ is an SOS polynomial,
$p(x)=\begin{bmatrix}p_{1}(x)&...&p_{m}(x)\end{bmatrix}^{T}\in(R[x])^{m}$ is a
vector of polynomials, and $p_{m+1}(x)\in R[x]$ is a scalar polynomial.

If $s_{1}(x)$ in (18) is strictly positive w.r.t. $x^{*}$, there exists a
controller defined by the rational function

$\displaystyle u_{\text{CBF},1}(x):=p(x)/s_{1}(x).$ (19)

###### Remark 3.4

As in Remark 1, we denote $B(x)$ as CBF for both the control and closed-loop
system.

###### Lemma 2

Given a CBF $B_{1}(x)$ for the control system (1) with $s_{1}(x)$ in (18)
strictly positive, the set $\mathcal{X}_{s,1}$ in (15) — if compact (cf.
Assumption 1) — is forward invariant w.r.t. system (2) using controller
$u(x)=u_{\text{CBF},1}(x)$ defined in (19).

A compact zero sublevel set of a scalar polynomial $B_{1}(x)$ is forward
invariant w.r.t. closed-loop system (2) with control $u(x)$, if $\nabla
B_{1}(x)^{T}\left(f(x)+G(x)u(x)\right)<0$ for all
$x\in\\{x\in\mathbb{R}^{n}\mid B_{1}(x)=0\\}$ (Blanchini, 1999, Nagumo’s
Theorem 3.1). From (18), we derive $s_{1}(x)\nabla B_{1}(x)^{T}f(x)+\nabla
B_{1}(x)^{T}G(x)p(x)<0$ for all $x\in\\{x\in\mathbb{R}^{n}\mid B_{1}(x)=0\\}$.
Since $s_{1}(x)$ is strictly positive w.r.t. $x^{*}$, the division by
$s_{1}(x)$ results in $s_{1}(x)\nabla
B_{1}(x)^{T}\left(f(x)+G(x)u_{\text{CBF},1}(x)\right)<0$ for all
$x\in\\{x\in\mathbb{R}^{n}\mid B_{1}(x)=0\\}$.

$\blacksquare$

##  4 SOS Construction of compatible CLF and multiple CBFs

In this section, we are interested in finding a CLF and multiple CBFs that are
compatible with each other. This is achieved by deriving a new set of SOS
constraints that result from restricting the rational controllers (13) and
(19) to be identical.

Multiple CBFs provide an additional flexibility to increase the volume of the
safe set $\mathcal{X}_{s}$, which is defined by the intersection of the zero
sublevel sets of all the CBFs. Intuitively, this makes sense since the
solution of a single CBF can be recovered by equating all CBFs.

###  4.1 Multiple Control Barrier Functions

Consider a set of CBFs $\\{B_{i}(x)\\}_{i\in[t]}$ for some
$t\in\mathbb{N}_{+}$, each defining a forward invariant set

$\displaystyle\mathcal{X}_{s,i}=\Bigl{\\{}x\in\mathbb{R}^{n}\mid B_{i}(x)\leq
0\Bigr{\\}}.$ (20)

Multiple CBFs are compatible with each other if they have the control-sharing
property (cf. Grammatico et al. (2013)).

###### Definition 4.5

Consider the closed-loop system (2) with feedback controller $u(x)$. A set of
CBFs $\\{B_{i}(x)\\}_{i\in[t]}$ with associated controllers
$u_{\text{CBF},i}(x)$ rendering (20) invariant has the control-sharing
property, if $u(x)=u_{\text{CBF},i}(x)$ for all $i\in[t]$.

The sets of compatible CBFs $\\{B_{i}(x)\\}_{i\in[t]}$ form a new forward
invariant set $\mathcal{X}_{s}$, called safe set, defined as the intersection
of the invariant sets from each $B_{i}(x)$:

$\displaystyle\mathcal{X}_{s}:=\bigcap_{i\in[t]}\mathcal{X}_{s,i}=\Bigl{\\{}x\in\mathbb{R}^{n}\mid
B_{i}(x)\leq 0\text{ }\forall i\in[t]\Bigr{\\}}.$ (21)

Consider a set of allowable states of the form

$\displaystyle\mathcal{X}_{a}=\Bigl{\\{}x\in\mathbb{R}^{n}\mid w_{i}(x)\leq
0\text{ }\forall i\in[t]\Bigr{\\}},$ (22)

where $w_{1}(x),...,w_{t}\in R[x]$ are polynomials. Note that for each
polynomial $w_{i}(x)$, that defines the allowable set $\mathcal{X}_{a}$, a
corresponding CBF can be assigned.

###### Assumption 2

The sets $\bigl{\\{}x\in\mathbb{R}\mid w_{i}(x)\leq 0\bigr{\\}}_{i\in[t]}$ are
compact.

To conclude safety, the resulting invariant set $\mathcal{X}_{s}$ needs to be
contained in the allowable set, i.e.
$\mathcal{X}_{s}\subseteq\mathcal{X}_{a}$. This can be achieved by restricting
the invariant set of each CBF to the zero sublevel set of the corresponding
$w_{i}(x)$ as stated by the following lemma.

###### Lemma 3

Consider a set of CBFs $\\{B_{i}(x)\\}_{i\in[t]}$ that have the control-
sharing property with input $u(x)$. If for every $i\in[t]$

$\displaystyle\mathcal{X}_{s,i}\subseteq\Bigl{\\{}x\in\mathbb{R}^{n}\mid
w_{i}(x)\leq 0\Bigr{\\}},$ (23)

then the closed-loop system (2) using controller $u(x)$ is safe w.r.t. the
safe set $\mathcal{X}_{s}$ (21) and the allowable set $\mathcal{X}_{a}$ (22).

For a closed-loop system to be safe, we need to show that $\mathcal{X}_{s}$ is
(i) forward invariant and (ii) $\mathcal{X}_{s}\subseteq\mathcal{X}_{a}$. Note
that the sets $\mathcal{X}_{s,i}$, $i\in[t]$ are compact due to (23). By Lemma
2, every $\mathcal{X}_{s,i}$, $i\in[t]$, is forward invariant w.r.t. the
closed-loop system, and the intersection of forward invariant sets is also
forward invariant. This proves the condition (i). The condition (ii) can be
seen from

$\displaystyle\mathcal{X}_{s}=\bigcap_{i\in[t]}\mathcal{X}_{s,i}\subseteq\bigcap_{i\in[t]}\Bigl{\\{}x\in\mathbb{R}^{n}\mid
w_{i}(x)\leq 0\Bigr{\\}}=\mathcal{X}_{a}.$

$\blacksquare$ Finally, we formalize (23) as the empty set condition

$\displaystyle\Bigl{\\{}x\in\mathbb{R}^{n}\mid w_{i}(x)\geq 0,-B_{i}(x)\geq
0\Bigr{\\}}=\emptyset$ (24)

for all $i\in[t]$. Under Theorem 2 and Assumption 2, condition (24) is
equivalent to the following SOS condition:

$\displaystyle B_{i}(x)-s_{4,i}(x)w_{i}(x)\in\Sigma[x]$ (25)

where $s_{4,i}\in\Sigma[x]$ is an SOS polynomial for all $i\in[t]$.

###  4.2 Compatible CLF and CBFs

The notion of stability and safety can be unified by the existence of both a
CLF and multiple CBFs that are compatible with each other via the following
control-sharing property.

###### Definition 4.6

Consider a closed-loop system (2) with feedback controller $u(x)$. A CLF with
associated controller $u_{\text{CLF}}(x)$ stabilizing the system to $x^{*}$
and set of CBFs $\\{B_{i}(x)\\}_{i\in[t]}$ with associated controllers
$u_{\text{CBF},i}(x)$ rendering (20) invariant have the control-sharing
property, if $u(x)=u_{\text{CLF}}(x)$ and $u(x)=u_{\text{CBF},i}(x)$ for all
$i\in[t]$.

With Definition 4.6 in place, we can state our main proposition:

###### Proposition 4

Given a CLF $V(x)$ and a set of CBFs $\\{B_{i}(x)\\}_{i\in[t]}$ that have the
control-sharing property with input $u(x)$. If condition (23) holds for every
$i\in[t]$, then the closed-loop system (2) with input $u(x)$ is asymptotically
stable and safe w.r.t. the allowable set $\mathcal{X}_{a}$ (22).

Lemma 1 and Lemma 3 prove stability and safety for the closed-loop system
using input $u(x)$.

$\blacksquare$

Based on conditions (12), (14), (18) and (25), the stability and safety
condition in Lemma 4 can be restated as the SOS constraints:

$\displaystyle\begin{array}[]{lll}-\nabla
V(x)^{T}\left(s_{1}(x)f(x)+G(x)p(x)\right)-l(x)^{2}\\\ \qquad\in\Sigma[x],\\\
-\nabla B_{i}(x)^{T}\left(s_{1}(x)f(x)+G(x)p(x)\right)\\\ \qquad-
p_{m+1,i}(x)B_{i}(x)\in\Sigma[x],&&\forall i\in[t]\\\
B_{i}(x)-s_{4,i}(x)w_{i}(x)\in\Sigma[x],&&\forall i\in[t]\\\
s_{1}(x)-\epsilon_{s1}\in\Sigma[x],\\\ \end{array}$ (32)

where $s_{1}(x),s_{4,i}(x)\in\Sigma[x]$ are SOS polynomials,
$p(x)=\begin{bmatrix}p_{1}(x)&...&p_{m}(x)\end{bmatrix}^{T}\in(R[x])^{m}$ is a
vector of polynomials, and $p_{m+1,i}\in R[x]$ is a scalar polynomial.

The resulting control input $u(x)=p(x)/s_{1}(x)$ renders the system
asymptotically stable on $\mathcal{X}_{a}$ and safe w.r.t. the allowable set
$\mathcal{X}_{a}$ and safe set $\mathcal{X}_{s}$.

##  5 SOS Program

In this section, we state the SOS program that finds a CLF and multiple CBFs
for the control system (1). The SOS program maximizes a surrogate of the
volume of $\mathcal{X}_{s}$ subject to the SOS condition (32).

###  5.1 Cost function

Ideally, the SOS program would optimize the CLF and CBFs over the volume of
the safe set $\mathcal{X}_{s}$:

$\displaystyle\begin{array}[]{ll}\underset{v}{\mbox{minimize}}&-vol\left(\Bigl{\\{}x\in\mathbb{R}^{n}\mid
B_{i}(x,v)\leq 0\text{ for }i\in[t]\Bigr{\\}}\right)\\\ \mbox{subject
to}&\text{SOS constraints~{}(\ref{eq:sos_condition})}\end{array}$ (35)

where $B_{i}(x,v)=Z_{B,i}(x)^{T}Q_{B,i}(v)Z_{B,i}(x)$ (cf. (3)). To make (35)
computationally tractable, the volume of ${\mathcal{X}}_{s}$ is surrogated by
the traces of the quadratic matrices encoding the $B_{i}(x)$. The cost then
becomes:

$\displaystyle c_{1}(v)=\sum_{i=1}^{t}tr(Q_{B,i}(v)).$ (36)

Furthermore, we have empirically observed that the outcome is improved by
restricting the CBF to a predefined center point. The justification for this
lies in the fact that the CBF condition is only valid when $B(x)=0$ and,
therefore, it does not restrict the CBF from taking arbitrary large values
when $B(x)\neq 0$.

Let us define a center point $(x_{c,i},B_{c,i})\in(\mathbb{R}^{n},\mathbb{R})$
for each CBF $B_{i}(x),\,i\in[t]$. Then, we define a cost function by the
deviation from that point:

$\displaystyle c_{2}(v)=\sum_{i=1}^{t}(B_{i}(x_{c,i})-B_{c,i})^{2}$ (37)

The cost function of the SOS program is set to

$\displaystyle c(v)=c_{1}(v)+c_{2}(v).$ (38)

###  5.2 SOS Program

Given the cost (38) and the SOS constraint (32), we summarize the SOS program
by:

$\displaystyle\begin{array}[]{lll}\mbox{find}&V,B_{i},s_{1},s_{4,i},p_{m+1,i}\in
R[x],\\\ &p\in(\mathbb{R}[x])^{m}\\\
\mbox{minimize}&\sum_{i=1}^{t}tr(Q_{B,i})+(B_{i}(x_{c,i})-B_{c,i})^{2}\\\
\mbox{subject to}&-\nabla V^{T}\left(s_{1}f+Gp\right)-l^{2}\in\Sigma[x]\\\
&-\nabla B^{T}_{i}\left(s_{1}f+Gp\right)-p_{m+1,i}B_{i}\\\
&\qquad\in\Sigma[x]&\forall i\in[t]\\\ &B_{i}-s_{4,i}w_{i}\in\Sigma[x]&\forall
i\in[t]\\\ &s_{1}-\epsilon_{s1}\in\Sigma[x]\\\ &s_{4,i}\in\Sigma[x]&\forall
i\in[t]\\\ &V\in\Sigma[x],V(x^{*})=0\end{array}$

where $(x_{c,i},B_{c,i})\in(\mathbb{R}^{n},\mathbb{R})$ are the predefined
center points for each CBF.

##  6 Algorithm

The SOS problem in Section 5.2 cannot directly be converted to an SDP problem
due to bilinear terms in its decision variables. However, an alternating
algorithm can find adequate solutions to this non-convex problem.

###  6.1 Alternating Algorithm

Consider an abstract optimization problem with $n_{c}$ SOS inequality
constraints defined by $f_{1}(v_{1},v_{2})\in(R[x])^{n_{c}}$ that are bilinear
in the decision variables $v_{1}\in\mathbb{R}^{n_{1}}$ and
$v_{2}\in\mathbb{R}^{n_{2}}$, and a cost function $f_{0}(v_{1})\in\mathbb{R}$
that only depends on $v_{1}$:

$\displaystyle\begin{array}[]{ll}\underset{v_{1},v_{2}}{\mbox{minimize}}&f_{0}(v_{1})\\\
\mbox{subject to}&f_{1}(v_{1},v_{2})\in(\Sigma[x])^{n_{c}}\end{array}$ (41)

The optimization problem (41) cannot be solved directly with an SOS program.
Instead, we propose a method that alternates between searching over one
variable while holding fixed the other.

Starting from an iteration $v_{1}^{(0)}$ and $v_{2}^{(0)}$, the algorithm at
iteration $k$ is defined by:

  * •

Step 1: Substitute $v_{2}=v_{2}^{(k-1)}$ and solve for $v_{1}$

$\displaystyle\begin{array}[]{ll}\underset{v_{1}}{\mbox{minimize}}&f_{0}(v_{1})\\\
\mbox{subject to}&f_{1}(v_{1},v_{2}^{(k-1)})\in(\Sigma[x])^{n_{c}}\end{array}$
(44)

  * •

Step 2: Substitute $v_{1}=v_{1}^{(k)}$ and solve for $\epsilon$ and $v_{2}$

$\displaystyle\begin{array}[]{ll}\underset{\epsilon,v_{2}}{\mbox{minimize}}&\epsilon\\\
\mbox{subject
to}&f_{1}(v_{1}^{(k)},v_{2})+s(v_{1}^{(k)})\epsilon\in(\Sigma[x])^{n_{c}}\end{array}$
(47)

The algorithm terminates when $|f_{0}(v_{1}^{(k)})-f_{0}(v_{1}^{(k-1)})|$ is
lower than a given threshold.

###### Lemma 5

If the initial pair $(v_{1}^{(0)},v_{2}^{(0)})$ is feasible for (41), then
feasibility is maintained along the alternating algorithms and the cost
$f_{0}(v_{1}^{(k)})$ in (44) is non-increasing with iteration $k$.

First, we show that if Step 1 in iteration $k$ is feasible for (44) then also
Step 2 in iteration $k$ must be feasible (47). For feasibility, we only need
to find a point $(v_{2},\epsilon)$, for which the constraint in Step 2 holds.
From Step 1, we know that $f_{1}(v_{1}^{(k)},v_{2}^{(k-1)})$ is an SOS
polynomial. But then, $(v_{2},\epsilon)=(v_{2}^{(k-1)},0)$ is a feasible point
of (47). As a consequence, $\epsilon$ must be equal or smaller than zero, i.e.
$\epsilon\leq 0$.

Next, we show that the solution of iteration $k$ is a feasible point of Step 1
in iteration $k+1$. First, note that $-s(v_{1}^{(k)})\epsilon$ is an SOS
polynomial vector since $\epsilon\leq 0$. But then also
$f_{1}(v_{1}^{(k)},v_{2}^{(k)})$ is an SOS polynomial. Therefore,
$(v_{1}^{(k)},v_{2}^{(k)})$ is a feasible point of Step 1 in iteration $k$.
Taken together, we showed that $(v_{1}^{(k)},v_{2}^{(k)})$ remain feasible for
(41). Therefore , $f_{0}(v_{1}^{(k+1)})$ is either equal or smaller than
$f_{0}(v_{1}^{(k)})$.

$\blacksquare$

###  6.2 Main algorithm

The main algorithm improves feasibility by introducing an operating region
$\mathcal{X}_{op}$. Under Assumption 2, there exists a scalar polynomial
$r(x)$ whose zero sublevel set is compact and strictly contains the allowable
set $\mathcal{X}_{a}$, i.e.

$\displaystyle\mathcal{X}_{s}\subseteq\mathcal{X}_{a}\subseteq{\mathcal{X}}_{op}:=\Bigl{\\{}x\in\mathbb{R}^{n}\mid
r(x)\leq 0\Bigr{\\}}.$ (48)

Note that all state trajectories of interest are contained in the operating
region $\mathcal{X}_{op}$. We will use this fact to slightly alter the empty
set conditions and the corresponding SOS constraints developed in the previous
chapters.

Given an allowable set $\mathcal{X}_{a}$ defined by $w_{i}(x)$ for $i\in[t]$
(22), positive definite SOS polynomials $l(x)\in\Sigma[x]$, a positive scalar
$\epsilon_{s1}>0$, center points
$(x_{c,i},B_{c,i})\in(\mathbb{R}^{n},\mathbb{R})$ for each $i\in[t]$, a region
of operation defined by $r(x)$ (48) and initial values $s_{1}(x)$, $p(x)$, and
$p_{m+1,i}(x)$ for $i\in[t]$, the main SOS algorithm is summarized below.

  * •

Step 1: Given a controller $u=p/s_{1}$ find a CLF and multiple CBFs.

$\displaystyle\begin{array}[]{lll}\mbox{find}&V,B_{i},s_{2},s_{3,i},s_{4,i}\in
R[x]\\\ \mbox{hold fixed}&s_{1},p_{m+1,i}\in R[x],p\in(\mathbb{R}[x])^{m}\\\
\mbox{minimize}&\sum_{i=1}^{t}tr(Q_{B,i})+(B_{i}(x_{c,i})-B_{c,i})^{2}\\\
\mbox{subject to}&-\nabla
V^{T}\left(s_{1}f+Gp\right)-l^{2}+s_{2}r\in\Sigma[x]\\\ &-\nabla
B^{T}_{i}\left(s_{1}f+Gp\right)-p_{m+1,i}B_{i}\\\
&\qquad+s_{3,i}r\in\Sigma[x]&\forall i\in[t]\\\
&B_{i}-s_{4,i}w_{i}\in\Sigma[x]&\forall i\in[t]\\\ &s_{2}\in\Sigma[x]\\\
&s_{3,i}\in\Sigma[x]&\forall i\in[t]\\\ &s_{4,i}\in\Sigma[x]&\forall
i\in[t]\\\ &V\in\Sigma[x],V(x^{*})=0\\\ \\\ \end{array}$

  * •

Step 2: Given a CLF and multiple CBFs find a controller $u=p/s_{1}$.

$\displaystyle\begin{array}[]{lll}\mbox{find}&\epsilon_{i}\in\mathbb{R},s_{1},s_{2},s_{3,i},p_{m+1,i}\in
R[x],\\\ &p\in(\mathbb{R}[x])^{m}\\\ \mbox{hold fixed}&V,B_{i}\in R[x]\\\
\mbox{minimize}&\sum_{i=1}^{t}\epsilon_{i}\\\ \mbox{subject to}&-\nabla
V^{T}\left(s_{1}f+Gp\right)-l^{2}+s_{2}r\in\Sigma[x]\\\ &-\nabla
B^{T}_{i}\left(s_{1}f+Gp\right)-p_{m+1,i}B_{i}\\\
&\qquad+s_{3,i}r+\epsilon_{i}\in\Sigma[x]&\forall i\in[t]\\\
&s_{1}-\epsilon_{s1}\in\Sigma[x]\\\ &s_{2}\in\Sigma[x]\\\
&s_{3,i}\in\Sigma[x]&\forall i\in[t]\\\ \end{array}$

##  7 Simulation

The alternating algorithm from Section 6.2 is implemented for a three-
dimensional vector field realization of a dc/ac power converter model. A
feedback controller that avoids unsafe states is of paramount importance for
this application since all electrical variables (e.g. voltages and currents)
need to be constrained at all times. With the dc voltage and ac currents as
states, the control system (1) is defined by

$\displaystyle f(x)=\begin{bmatrix}-0.05x_{1}-57.9x_{2}+0.00919x_{3}\\\
1710x_{1}+314x_{3}\\\ -0.271x_{1}-314x_{2}\end{bmatrix}$ (49)

and

$\displaystyle G(x)=\begin{bmatrix}0.05-57.9x_{2}&-57.9x_{3}\\\
1710+1710x_{1}&0\\\ 0&1710+1710x_{1}\end{bmatrix}.$ (50)

The allowable set $\mathcal{X}_{a}$ is defined by two polynomials $w_{1}(x)$
and $w_{2}(x)$ encoding the state constraints for current and voltage
respectively:

$\displaystyle\begin{split}w_{1}(x)&=(x_{1}+0.3)^{2}+(x_{2}/20)^{2}+(x_{3}/20)^{2}-0.5^{2}\\\
w_{2}(x)&=(x_{1}/20)^{2}+x_{2}^{2}+x_{3}^{2}-1.2^{2}.\end{split}$ (51)

An operating region $\mathcal{X}_{op}$ in (48) can be parameterized by:

$\displaystyle r(x)=(x_{1}/0.8)^{2}+(x_{2}/1.2)^{2}+(x_{3}/1.2)^{2}-1.8.$ (52)

A $r(x)$ that tighter fits the allowable set results in a faster convergence
of the alternating algorithm. For illustration purposes, however, we selected
a suboptimal operation region. The center points of the CBFs are given by:
$x_{c,i}=\begin{bmatrix}-0.3&0&0\end{bmatrix}^{T}$ and $B_{c,i}=-10$ for
$i\in[2]$. The degrees of the polynomials involved in the algorithm are
summarized as follows:

$V(x)$ | $B_{i}(x)$ | $s_{1}(x)$ | $s_{2}(x)$ | $s_{3,i}(x)$ | $s_{4,i}(x)$ | $p(x)$ | $p_{m+1,i}(x)$
---|---|---|---|---|---|---|---
4 | 4 | 2 | 6 | 6 | 3 | 4 | 2

For $V(x)$ and $B_{i}(x)$, the odd degree coefficients are set to zero:

  * •

$deg\left(V(x)\right)=(2,4)$

  * •

$deg\left(B_{1}(x-x_{c,1})\right)=(0,2,4)$

  * •

$deg\left(B_{2}(x-x_{c,2})\right)=(0,2,4)$.

Figure 1 illustrates the sequence of safe sets calculated at each iteration of
the algorithm. The computed CLF and CBFs not only prove that there exists a
stabilizing feedback controller but also guarantee that such controllers
renders the system safe w.r.t. a safe set $\mathcal{X}_{s}$. By using two
CBFs, the algorithm terminates with a safe set that tightly fits the allowable
set $\mathcal{X}_{a}$.

The number of variables and constraints involved in the resulting SDP are
summarized as follows:

| Step 1 | Step 2
---|---|---
Variables | 337 | 372
Constraints | 5188 | 4892

The algorithm was completed in 90 seconds. For each iteration, the SDP
algorithm took on average 33 inner iterations to solve step 1 and on average
21 inner iterations to solve step 2.

Figure 1: shows the evolution of the sequence of safe sets $\mathcal{X}_{s}$
projected to the $x_{1}$ and $x_{2}$ coordinates over 6 iterations. The
boundary of the allowable set $\mathcal{X}_{a}$ is shown in red. The sequence
of zero-level sets of the CBF $B_{1}(x)$ and $B_{2}(x)$ are indicated in gray.

##  8 Conclusions and future work

This paper presented a framework to combine stability and safety conditions
using CLF and multiple CBFs. By employing two versions of the
Positivstellensatz, we synthesized a controller used to prove compatibility
between CLF and CBFs. We then formalized SOS constraints that encode
compatible CLF and CBF conditions. Finally, we proposed an algorithm that
solves the resulting SOS program by iteratively solving two SDPs.

For future work, we plan to study the computational complexity of our
algorithm and also to explore a unified framework, that proves stability and
safety with weaker CLF and CBF conditions. Moreover, we intend to incorporate
noise-robustness into our SOS formulation, as demonstrated by Kang et al.
(2023).

## References

  * Ames et al. (2019) Ames, A.D., Coogan, S., Egerstedt, M., Notomista, G., Sreenath, K., and Tabuada, P. (2019).  Control barrier functions: Theory and applications.  In _2019 18th European control conference (ECC)_ , 3420–3431. IEEE.
  * Anghel et al. (2013) Anghel, M., Milano, F., and Papachristodoulou, A. (2013).  Algorithmic construction of Lyapunov functions for power system stability analysis.  _IEEE Transactions on Circuits and Systems I: Regular Papers_ , 60(9), 2533–2546.
  * Blanchini (1999) Blanchini, F. (1999).  Set invariance in control.  _Automatica_ , 35(11), 1747–1767.
  * Bochnak et al. (2013) Bochnak, J., Coste, M., and Roy, M.F. (2013).  _Real algebraic geometry_ , volume 36.  Springer Science & Business Media.
  * Clark (2021) Clark, A. (2021).  Verification and synthesis of control barrier functions.  In _2021 60th IEEE Conference on Decision and Control (CDC)_ , 6105–6112. IEEE.
  * Grammatico et al. (2013) Grammatico, S., Blanchini, F., and Caiti, A. (2013).  Control-sharing and merging control Lyapunov functions.  _IEEE Transactions on Automatic Control_ , 59(1), 107–119.
  * Isaly et al. (2022) Isaly, A., Ghanbarpour, M., Sanfelice, R.G., and Dixon, W.E. (2022).  On the feasibility and continuity of feedback controllers defined by multiple control barrier functions for constrained differential inclusions.  In _2022 American Control Conference (ACC)_ , 5160–5165. IEEE.
  * Isidori (1995) Isidori, A. (1995).  _Nonlinear Control Systems_.  Springer.
  * Kang et al. (2023) Kang, S., Chen, Y., Yang, H., and Pavone, M. (2023).  Verification and synthesis of robust control barrier functions: Multilevel polynomial optimization and semidefinite relaxation.  _arXiv preprint arXiv:2303.10081_.
  * Khalil (2002) Khalil, H. (2002).  _Nonlinear Systems_.  Pearson Education. Prentice Hall.
  * Korda et al. (2014) Korda, M., Henrion, D., and Jones, C.N. (2014).  Convex computation of the maximum controlled invariant set for polynomial control systems.  _SIAM Journal on Control and Optimization_ , 52(5), 2944–2969.
  * Kundu et al. (2019) Kundu, S., Geng, S., Nandanoori, S., Hiskens, I.A., and Kalsi, K. (2019).  Distributed barrier certificates for safe operation of inverter-based microgrids.  1042–1047.
  * Laurent (2009) Laurent, M. (2009).  Sums of squares, moment matrices and optimization over polynomials.  _Emerging applications of algebraic geometry_ , 157–270.
  * Tan and Packard (2004) Tan, W. and Packard, A. (2004).  Searching for control Lyapunov functions using sums of squares programming.  _sibi_ , 1(1).
  * Tan and Dimarogonas (2022) Tan, X. and Dimarogonas, D.V. (2022).  Compatibility checking of multiple control barrier functions for input constrained systems.  _arXiv preprint arXiv:2209.02284_.
  * Wang et al. (2022) Wang, H., Margellos, K., and Papachristodoulou, A. (2022).  Safety verification and controller synthesis for systems with input constraints.  _arXiv preprint arXiv:2204.09386_.
  * Wieland and Allgöwer (2007) Wieland, P. and Allgöwer, F. (2007).  Constructive safety using control barrier functions.  _IFAC Proceedings Volumes_ , 40(12), 462–467.
  * Xu (2016) Xu, X. (2016).  Control sharing barrier functions with application to constrained control.  In _2016 IEEE 55th Conference on Decision and Control (CDC)_ , 4880–4885. IEEE.