# Koopman Operator Approximation under Negative Imaginary Constraints

M. A. Mabrok, Ilyasse Aksikas and Nader Meskin M. A. Mabrok and Ilyasse
Aksikas are with Mathematics Program, Department of Mathematics, Statistics
and Physics, College of Arts and Sciences, Qatar University, P.O. box 2713,
Doha, Qatar<EMAIL_ADDRESS><EMAIL_ADDRESS>Nader Meskin is with
Electircal Engineering Department, College of Engineering, Qatar University,
<EMAIL_ADDRESS>
###### Abstract

Nonlinear Negative Imaginary (NI) systems arise in various engineering
applications, such as controlling flexible structures and air vehicles.
However, unlike linear NI systems, their theory is not well-developed. In this
paper, we propose a data-driven method for learning a lifted linear NI
dynamics that approximates a nonlinear dynamical system using the Koopman
theory, which is an operator that captures the evolution of nonlinear systems
in a lifted high-dimensional space. The linear matrix inequality that
characterizes the NI property is embedded in the Koopman framework, which
results in a non-convex optimization problem. To overcome the numerical
challenges of solving a non-convex optimization problem with nonlinear
constraints, the optimization variables are reformatted in order to convert
the optimization problem into a convex one with the new variables. We compare
our method with local linearization techniques and show that our method can
accurately capture the nonlinear dynamics and achieve better control
performance. Our method provides a numerically tractable solution for learning
the Koopman operator under NI constraints for nonlinear NI systems and opens
up new possibilities for applying linear control techniques to nonlinear NI
systems without linearization approximations.

##  I Introduction

Negative imaginary (NI) systems are systems that have a negative imaginary
frequency response when subjected to a positive feedback interconnection [1,
2, 3]. NI systems theory has been widely used for analysis and design of
linear-time-invariant (LTI) control systems, especially for flexible
structures and air vehicles [1, 2, 4, 5]. NI systems have several desirable
properties, such as robust stability under positive feedback, disspativity
with respect to collocated inputs and outputs, and existence of optimal
controllers [3]. A central result in the NI theory is that the positive
feedback interconnection between an NI system with a transfer function matrix
$G(s)$ and a strictly NI system with a transfer function matrix $\bar{G}(s)$,
as shown in Fig. 1, results in a robustly stable feedback interconnection [1,
2]. This implies that if a given system is NI, then, it is a great advantage
to synthesis a strictly NI controller to guarantee robust stability.

There have been various approaches to address the NI control synthesis problem
[6, 7, 8]. The paper [8] proposes a synthesis methodology for non-linear
systems using a library of controllers parametrized and strictly negative
imaginary, optimized through Sequential Quadratic Programming. This
methodology can be applied to different non-linear systems. In [6], a negative
imaginary and strict negative imaginary lemma is used to construct a static
controller that ensures robust stability against strict negative imaginary
uncertainty. In [7], a data-driven controller synthesis methodology is
proposed that uses measured frequency response data to construct the
controller response and transfer function. The data driven linear NI system
identification problem was addressed in [9], where a modified subspace system
identification algorithm that guarantees the negative imaginary property in
the identified model by imposing constraints to ensure stability and negative
imaginariness was proposed.

Another approach is the sufficient conditions developed in [10] for designing
various types of output feedback controllers, including static, dynamic, and
observer-based feedback controllers, based on the above design conditions. The
desired controllers can be found using various numerical algorithms. Discrete-
time negative imaginary systems were introduced in [11, 12] and the authors
investigated the relationship between this concept and discrete-time positive
real system theory. Discrete-time negative imaginary systems are defined in
terms of a sign condition that must be met in a domain of analyticity of the
transfer function. In [13], similar discrete-time results are explored without
the restriction of a symmetric transfer function.

Negative imaginary systems’ theory provides a way to analyze robustness and
design robust controllers. In other words, the theory can be used to design
controllers that are able to maintain stability and performance even in the
presence of uncertainties in the system, such as unmodeled spillover dynamics
or variations in resonant frequencies and damping levels [14].

However, many real-world systems are nonlinear in nature and cannot be
adequately modeled by LTI systems. For example, the dynamics of a mass-spring-
damper system with a nonlinear spring or a nonlinear damper are nonlinear NI
systems [15]. Nonlinear NI systems pose significant challenges for analysis
and control, as the existing NI systems theory does not directly apply to
them. Therefore, there is a need to develop a general framework for nonlinear
NI systems that can capture their essential features and enable their
effective control.

Figure 1: If the transfer function matrix of a plant, $G(s)$, is Negative-
Imaginary (NI) and the transfer function matrix of a controller, $\bar{G}(s)$,
is Strictly Negative-Imaginary (SNI), then the positive-feedback
interconnection will be internally stable if and only if the DC gain condition
$\lambda_{max}(G(0)\bar{G}(0))<1$ is met.

One possible approach to deal with nonlinear systems is to use the Koopman
operator, which is a linear operator that describes the evolution of scalar
observables of nonlinear systems in an infinite-dimensional Hilbert space
[16]. The Koopman operator has attracted considerable attention in recent
years as a powerful tool for nonlinear dynamics modeling and control [17, 16,
18]. By using the Koopman operator, one can lift the nonlinear dynamics to a
high-dimensional space where they become linear, and then apply linear
techniques to analyze and control them. However, learning the Koopman operator
from data is generally a challenging task, as it involves solving a non-convex
optimization problem with nonlinear constraints [16].

In this paper, we propose a data-driven method for learning the Koopman
operator under NI constraints for nonlinear NI systems. We design a data-
driven method for learning a lifted linear NI dynamics that approximates a
nonlinear dynamical system using the Koopman theory. To overcome the numerical
challenges of solving a non-convex optimization problem with nonlinear
constraints, we use a change of variable technique, which allow us to
reformulate the problem as a convex optimization that can be solved
efficiently. We compare our method with local linearization techniques and
show that our method can accurately capture the nonlinear dynamics and achieve
better control performance. Our method provides a numerically tractable
solution for learning the Koopman operator under NI constraints for nonlinear
NI systems and opens up new possibilities for applying linear control
techniques to nonlinear systems without linearization approximations.

The rest of the paper is organized as follows. Section II reviews some
preliminaries on NI systems theory and Koopman operator theory. Section III
presents our main results on learning the Koopman operator under NI
constraints for nonlinear NI systems. Section IV illustrates our method on a
numerical example; a mass-spring-damper system. Section V concludes the paper
and discusses some future directions.

##  II Preliminaries and Notation

This section introduces the definitions of negative imaginary systems and
Koopman operator.

###  II-A Negative Imaginary Systems

Consider the following LTI system

$\displaystyle\dot{x}(t)=Ax(t)+Bu(t),$ (1) $\displaystyle y(t)=Cx(t)+Du(t),$
(2)

where $A\in\mathbb{R}^{n\times n},B\in\mathbb{R}^{n\times
m},C\in\mathbb{R}^{m\times n},$ $D\in\mathbb{R}^{m\times m},$ and with the
square transfer function matrix $G(s)=C(sI-A)^{-1}B+D$. The transfer function
matrix $G(s)$ is said to be strictly proper if $G(\infty)=D=0$. We will use
the notation $\begin{bmatrix}\begin{array}[]{c|c}A&B\\\ \hline\cr
C&D\end{array}\end{bmatrix}$ to denote the state space realization (1), (2). A
NI system is defined as follows:

###### Definition 1

[2] A square transfer function matrix $G(s)$ is NI if the following conditions
are satisfied:

  1. 1.

$G(s)$ has no pole in $\textmd{Re}[s]>0$.

  2. 2.

For all $\omega>0$ such that $s=j\omega$ is not a pole of $G(s)$,

$j\left(G(j\omega)-G(j\omega)^{\ast}\right)\geq 0.$ (3)

  3. 3.

If $s=j\omega_{0}$ with $\omega_{0}>0$ is a pole of $G(s)$, then it is a
simple pole and the residue matrix $K=\underset{s\longrightarrow
j\omega_{0}}{\lim}(s-j\omega_{0})jG(s)$ is Hermitian and positive
semidefinite.

  4. 4.

If $s=0$ is a pole of $G(s)$, then $\underset{s\longrightarrow
0}{\lim}s^{k}G(s)=0$ for all $k\geq 3$ and $\underset{s\longrightarrow
0}{\lim}s^{2}G(s)$ is Hermitian and positive semidefinite.

###### Definition 2

[19] A square transfer function matrix $G(s)$ is SNI if the following
conditions are satisfied:

  1. 1.

${G}(s)$ has no pole in $\textmd{Re}[s]\geq 0$.

  2. 2.

For all $\omega>0$, $j\left({G}(j\omega)-{G}(j\omega)^{\ast}\right)>0$.

###  II-B Negative Imaginary lemma

The negative imaginary lemma is a result that describes NI systems using a
pair of LMIs, similar to how the positive-real lemma does [20, 21]. This
result was presented in [22, 1] and it also covers the case where there are
poles on the imaginary axis except at zero [19].

###### Lemma 1

(See [19]) Let $\begin{bmatrix}\begin{array}[]{c|c}A&B\\\ \hline\cr
C&D\end{array}\end{bmatrix}$ be a minimal state space realization of a
transfer function matrix $G(s)$. Then $G(s)$ is NI if and only if
$\textmd{det}(A)\neq 0$, $D=D^{T}$ and there exists a real matrix $P>0$ such
that

$AP+PA^{*}\leq 0,$ (4) $B=-APC^{*}.$ (5)

As shown in Fig. 1, there is a positive feedback interconnection between an NI
system with transfer function matrix $G(s)$ and an SNI system with transfer
function matrix $\bar{G}(s)$. Additionally, let’s assume that the transfer
function matrix $G(s)$ has a minimal state space realization,
$\begin{bmatrix}\begin{array}[]{c|c}A&B\\\ \hline\cr
C&D\end{array}\end{bmatrix},$ and $\bar{G}(s)$ has a minimal state space
realization $\begin{bmatrix}\begin{array}[]{c|c}\bar{A}&\bar{B}\\\
\hline\cr\bar{C}&\bar{D}\end{array}\end{bmatrix}.$

###  II-C Nonlinear Negative Imaginary Systems

The theory of nonlinear NI is not yet well investigated. Preliminary work has
been carried out in [23, 24, 15]. In this subsection, we are presenting the
definition of the nonlinear NI.

Consider the following nonlinear dynamical system

$\displaystyle\dot{x}(t)=f(x(t),u(t)),$ (6) $\displaystyle y(t)=h(x(t)),$ (7)

where $x(t)\in\mathbb{R}^{n}$, $u(t)\in\mathbb{R}^{m}$, and
$y(t)\in\mathbb{R}^{m}$ are the system state, control input, and output,
respectively.

###### Definition 3

The nonlinear system (7) is said to be NI if there exists a positive storage
function $S_{f}:\mathbb{R}^{n}\rightarrow\mathbb{R}$ of a class $C^{1}$ such
that

$S_{f}(x(t))\leq S_{f}(x(0))+\int_{0}^{t}u(s)^{T}\dot{y}(s)ds,$ (8)

for all $t>0$. Alternatively, the integral condition (8) can be represented by
the differential condition

$\dot{S_{f}}(x(t))\leq u(t)^{T}\dot{y}(t).$ (9)

It is important to note that the way output measurement is defined in negative
imaginary systems is different from passive systems. Negative imaginary
systems use position or acceleration measurements, while passive systems use
velocity sensors. Negative imaginary systems also allow for transfer functions
with a relative degree of up to two, while passive systems permit relative
degree zero or one. This is seen in systems that use force actuators and
position sensors, like robotics and nano-positioning systems. Stability
results for negative feedback interconnections between passive and strictly
passive systems do not apply to positive position feedback interconnections
between negative imaginary and strictly negative imaginary systems. This was
shown in [25] for physical systems like spring mass damper systems and RLC
electrical networks by comparing stability results for feedback
interconnections of negative imaginary systems with those for passive systems.

There are several important mechanical systems that satisfy the nonlinear
negative imaginary property. For instance, a mass-spring-damper system with
positive nonlinear damping coefficient is a nonlinear negative imaginary from
force input and displacement output. A more general class of nonlinear
negative imaginary systems is the class of systems that are governed by the
following Euler-Lagrangian dynamics,

$M(q)\ddot{q}+C(q,\dot{q})\dot{q}+G(q)=F.$ (10)

It has been shown in [15] that the dynamical systems governed by Equation (10)
are indeed nonlinear negative imaginary systems.

###  II-D Koopman theory

In this section, a brief introduction for the Koopman operator given in [26]
is presented. The Koopman operator is a mathematical tool used in the analysis
of dynamical systems. It was introduced by Bernard Koopman in 1931, and has
since become an important tool in the study of nonlinear systems. In
particular, in the context of dynamical systems, the Koopman operator provides
a powerful tool for studying the evolution of the system over time. By mapping
the system’s state space to a higher dimensional space, the operator allows
for the analysis of nonlinear dynamics through linear methods. This is
particularly useful in cases where traditional nonlinear analysis techniques,
such as numerical simulations or bifurcation analysis, become computationally
expensive or infeasible. One key application of the Koopman operator is in
data-driven modeling of complex systems. By using data from the system, the
Koopman operator can be estimated and used to approximate the system’s
behavior. This approach has been particularly successful in fields such as
fluid dynamics, where traditional modeling techniques are often impractical
due to the complexity of the system.

Consider the following nonlinear difference equation

$\displaystyle\left\\{\begin{array}[]{ll}x(j+1)&=f\left(x(j),u(j)\right),\\\
y(j)&=g\left(x(j)\right),\end{array}\right.$ (13)

where $j$ represents the discrete time, $u\in\mathbb{R}^{m}$ represents the
input, $x\in\mathbb{R}^{n}$ represents the state of the dynamical system,
$y\in\mathbb{R}^{l}$ represents the output, and
$f(x,u):\mathbb{R}^{n+m}\to\mathbb{R}^{n}$ and
$g(x):\mathbb{R}^{n}\to\mathbb{R}^{l}$ are the nonlinear functions. Let $\xi$
to be defined as, $\xi:=\left[\begin{array}[]{c}x\\\
u\end{array}\right]\in\mathbb{R}^{n+m}.$ Also, define $\mathbb{F}$ to be the
nonlinear operator given by
$\mathbb{F}(\xi):=\left[\begin{array}[]{c}f(x,u)\\\ \Xi(u)\end{array}\right],$
where $\Xi$ is the time-shift operator, i.e $\Xi(u(j)):=u(j+1).$ This implies
that the time evolution of $\xi$ can be given as
$\xi(j+1)=\mathbb{F}(\xi(j)).$ Also, define $\phi_{\mathrm{inf}}(\xi)$ to be
the infinite-dimensional lifting function in the new state space
representation,

$\displaystyle\phi_{\mathrm{inf}}(\xi)=\left[\begin{array}[]{c}\phi_{1}(\xi)\\\
\phi_{2}(\xi)\\\ \vdots\end{array}\right].$ (17)

Now, we can define the Koopman operator $\mathcal{K}$ as

$\displaystyle\mathcal{K}(\phi_{\mathrm{inf}}(\xi)):=\phi_{\mathrm{inf}}(\mathbb{F}(\xi)),$

where the time evolution of the lifting function $\phi_{\mathrm{inf}}(\xi)$ is
given as follows;

$\displaystyle\phi_{\mathrm{inf}}(\xi(j+1))=\mathcal{K}(\phi_{\mathrm{inf}}(\xi(j))).$
(18)

which is a linear mapping defined on the infinite-dimensional state space.

In order to approximate the infinite-dimensional Koopman operator by a finite-
dimensional approximation, we define the $N_{\phi}$-dimensional lifting
function $\phi(\xi):\mathbb{R}^{n+m}\to\mathbb{R}^{N_{\phi}}$ as

$\displaystyle\phi(\xi)=\left[\begin{array}[]{c}\phi_{1}(\xi)\\\ \vdots\\\
\phi_{N_{\phi}}(\xi)\end{array}\right]\in\mathbb{R}^{N_{\phi}}.$ (22)

Now, define $U\in\mathbb{R}^{N_{\phi}\times N_{\phi}}$ to be a finite-
dimensional matrix that approximates the Koopman operator $\mathcal{K}$, i.e.,
$U$ minimizes the following norm

$\displaystyle\|U\phi(\xi)-\phi(\mathcal{F}(\xi))\|.$ (23)

This implies that

$\displaystyle\phi(\xi(j+1))\approx U\phi(\xi(j)),$ (24)

which describes the behavior of $\phi_{\mathrm{inf}}(\xi)$, defined in (18).
Similar to the approach in [26], we consider a class of the lifting function

$\displaystyle\phi(\xi)=\left[\begin{array}[]{c}\psi(x)\\\
u\end{array}\right]\in\mathbb{R}^{N+m},$ (27)

where $\psi(x):\mathbb{R}^{n}\to\mathbb{R}^{N}$ is $N$-dimensional lifting
function given by

$\displaystyle\psi(x)=\left[\begin{array}[]{c}\psi_{1}(x)\\\ \vdots\\\
\psi_{N}(x)\end{array}\right]\in\mathbb{R}^{N},$ (31)

and $N+m=N_{\phi}$ holds. The matrix $U$ can be given as

$\displaystyle
U=\left[\begin{array}[]{cc}A_{d}&B_{d}\end{array}\right]\in\mathbb{R}^{(N)\times(N+m)},$
(33)

where $A_{d}\in\mathbb{R}^{N\times N}$ and $B_{d}\in\mathbb{R}^{N\times m}$.
This implies that (24) can be rewritten as

$\displaystyle\psi(x(j+1))\approx A_{d}\psi(x(j))+B_{d}u(j).$ (34)

In order to formulate the output equation in Koopman operator, the
approximation of the output equation in (13) can be approximated as follows;

$\displaystyle y(j)\approx C_{d}\psi(x(j)),$

where $C_{d}\in\mathbb{R}^{l\times N}$. The final linear state-space model can
be defined as follows;

$\displaystyle\left\\{\begin{array}[]{ll}\psi(j+1)&=A_{d}\psi(j)+B_{d}u(j),\\\
y(j)&=C_{d}\psi(j).\end{array}\right.$ (37)

##  III Formulating the Koopman NI constraints

In this section, the NI constraints for discrete dynamical systems are
formulated and then imposed into the learning problem for obtaining the
approximated Koopman operator given in (24). Consider the following discrete
time system

$\displaystyle x(j+1)=A_{d}x(j)+B_{d}u(j),$ (38) $\displaystyle
y(j)=C_{d}x(j)+D_{d}u(j),$ (39)

where $A_{d}\in\mathbb{R}^{n\times n},B_{d}\in\mathbb{R}^{n\times
m},C_{d}\in\mathbb{R}^{m\times n}$ and $D_{d}\in\mathbb{R}^{m\times m}.$ The
NI lemma given in Lemma 1, which is formulated in continuous-time and
therefore, a bilinear transformation, as presented in [9], in the following
form,

$\displaystyle A$ $\displaystyle=\frac{1}{T}(I+A_{d})^{-1}(A_{d}-I),\ \ \
B=\frac{1}{\sqrt{T}}(I+A_{d})^{-1}B_{d},$ $\displaystyle C$
$\displaystyle=\frac{1}{\sqrt{T}}C_{d}(I+A_{d})^{-1},\ \ \ \
D=D_{d}-C_{d}(I+A_{d})^{-1}B_{d},$

with $T$ denoting a sampling time is used to transform the conditions (4) into
corresponding discrete time conditions. The LMI (4) is transformed from a
continuous-time to a discrete-time form as follows:

$\displaystyle AP+PA^{T}\leq 0$ (40) $\displaystyle\Leftrightarrow$
$\displaystyle(I+A_{d})^{-1}(A_{d}-1)P+P(A_{d}^{T}-1)(I+A_{d}^{T})^{-1}\leq
0,$ $\displaystyle\Leftrightarrow$
$\displaystyle(A_{d}P-P)(I+A_{d}^{T})+(P+A_{d}P)(A_{d}^{T}-I)\leq 0,$
$\displaystyle\Leftrightarrow$ $\displaystyle
A_{d}P+A_{d}PA_{d}^{T}-P-PA_{d}^{T}+PA_{d}^{T}+A_{d}PA_{d}^{T}$
$\displaystyle\;\;\;-P-A_{d}P\leq 0,$ $\displaystyle\Leftrightarrow$
$\displaystyle A_{d}PA_{d}^{T}-P+A_{d}PA_{d}^{T}-P\leq 0,$
$\displaystyle\Leftrightarrow$ $\displaystyle A_{d}PA_{d}^{T}-P\leq 0.$ (41)

Also,

$\displaystyle B$ $\displaystyle=-APC^{T}$ $\displaystyle\Leftrightarrow
B_{d}$ $\displaystyle=-\frac{1}{T}(A_{d}-I)P(I+A_{d}^{T})^{-1}C_{d}^{T}.$ (42)

###  III-A Koopman operator learning problem

Similar approach in [26], define the following data matrices, which are
generated from $L$ measurements,

$\displaystyle\Omega$
$\displaystyle:=[u(j),u(j+1),\dots,u(j+L-1)]\in\mathbb{R}^{m\times L},$ (43)
$\displaystyle Y$
$\displaystyle:=[y(j),y(j+1),\dots,y(j+L-1)]\in\mathbb{R}^{m\times L},$ (44)
$\displaystyle\Theta$
$\displaystyle:=[\psi(j),\psi(j+1),\dots,\psi(j+L-1)]\in\mathbb{R}^{N\times
L},$ (45) $\displaystyle\Theta_{+}$
$\displaystyle:=[\psi(j+1),\psi(j+2),\dots,\psi(j+L)]\in\mathbb{R}^{N\times
L}$ (46)

This implies that the Koopman operator learning problem can be formulated as
follows:

###### Problem 1

—s—[3] A_d,B_d , C_dJ_1(A_d, B_d)+J_2(C_d) P ≥0 A_dPA_d^T-P ≤0 B_d
=-1T(A_d-I)P(I+A_d^T)^-1C_d^T where

$\displaystyle J_{1}(A_{d},B_{d})$
$\displaystyle:=\left\|\Theta_{+}-\left[\begin{array}[]{ll}A_{d}&B_{d}\end{array}\right]\left[\begin{array}[]{c}\Theta\\\
\Omega\end{array}\right]\right\|_{F}^{2},$ (50) $\displaystyle J_{2}(C_{d})$
$\displaystyle:=\left\|Y_{j}-C_{d}\Theta\right\|_{F}^{2}.$ (51)

The optimization problem given in (1) is non-convex due to the non-convexity
of the constraints. In order to overcome this issue, we replace the matrix
inequality (41) by a strict inequality with a positive $\alpha>0$ such that
the following holds,

$\displaystyle P$ $\displaystyle\geq\alpha I$ (52) $\displaystyle
A_{d}PA_{d}^{T}-P$ $\displaystyle\leq-\alpha I.$ (53)

Now, using Schur complements, the LMIs given in (52) and (53) can be written
as follows;

$\displaystyle\left[\begin{array}[]{cc}P-\alpha I&A_{d}P\\\
PA_{d}^{T}&P\end{array}\right]\geq 0\ \
\Rightarrow\left[\begin{array}[]{cc}P-\alpha I&Q\\\ Q^{T}&P\end{array}\right]$
$\displaystyle\geq 0,$ (58)

where $Q=A_{d}P$.

The LMI (58) is now convex in the variables $Q$ and $P$. Furthermore, the
matrix $A_{d}$ can be recovered from $A_{d}=QP^{-1}$. However, the overall
optimization problem is still non-convex because the variables in the cost
function still ${A_{d},B_{d},C_{d}}$. To tackle this issue, the identification
of the matrix $C_{d}$ is considered as a separate optimization problem by
minimizing the cost function given in (51), which has an optimal solution that
is given as

$C_{d}=Y\Theta^{\dagger},$ (59)

where ${(\cdot)}^{\dagger}$ donates the Moore-Penrose pseudoinverse of
$\Theta$.

Second, two weighting matrices $W$ and $\hat{W}$ are introduced into the cost
function (50) as follows;

$\displaystyle J_{1}(A_{d},B_{d})$
$\displaystyle:=\left\|W\left(\Theta_{+}-\left[\begin{array}[]{ll}A_{d}&B_{d}\end{array}\right]\left[\begin{array}[]{c}\Theta\\\
\Omega\end{array}\right]\right)\hat{W}\right\|_{F}^{2},$ (63)

where $\hat{W}$ is selected as:

$\displaystyle\hat{W}=$ $\displaystyle\begin{bmatrix}\Theta\\\
\Omega\end{bmatrix}^{T}\times\left(\begin{bmatrix}\Theta\\\
\Omega\end{bmatrix}\begin{bmatrix}\Theta\\\
\Omega\end{bmatrix}^{T}\right)^{{\dagger}}\times\begin{bmatrix}P&0\\\ 0&I\\\
\end{bmatrix}.$ (64)

It is assumed that there exists sufficient data to be collected, which allows
us to assume that $\begin{bmatrix}\Theta\\\ \Omega\\\ \end{bmatrix}$ to be
full row rank. This implies that the cost function in (63) can be written as
follows:

$\displaystyle J_{1}(P,Q,B_{d})$
$\displaystyle:=\left\|W\Theta_{+}\hat{W}-W\left[\begin{array}[]{ll}Q&B_{d}\end{array}\right]\right\|_{F}^{2}.$
(66)

Therefore, the optimization problem can now be written as follows;

###### Problem 2

—s—[3] P, Q, B_d∥ WΘ_+ ^W\- W[QB_d]∥_F^2. [ P-αIQQ^TP] ≥0.

where $\hat{W}$ is defined as in (64), and the matrix $W$ is a weighting
matrix, which can be treated as a design parameter to tune the optimization
problem.

The optimization problem in (2) is now convex in the variables ${P,Q,B_{d}}$.
Using any convex optimization solver, such as CVX, can accurately compute the
variables ${P,Q,B_{d}}$, and therefore, we can compute the matrix $A_{d}$ from
the fact that $Q=A_{d}P$ as $A_{d}=QP^{-1}$. Note that the matrix $C_{d}$ is
computed separately in (59).

##  IV Example: Mass-Spring-Damper System

Consider a nonlinear mass-spring-damper system

$\ddot{z}(t)+\beta(z(t),\dot{z}(t))\dot{z}(t)+K(z(t))=u(t),\quad y(t)=z(t).$
(67)

where $\beta(z(t),\dot{z}(t))$ and $k(z(t))$ donate the friction coefficient
and the spring stiffness, respectively.

$m$$u$ Figure 2: Spring-damper-mass system

It can be easily shown that the above system is a nonlinear NI system by
choosing a storage function as follows;

$V(z(t),\dot{z}(t))=\frac{1}{2}m\dot{z}^{2}(t)+\int_{0}^{z}(t)k(\xi)d\xi.$

The time derivative of $V(z(t),\dot{z}(t))$ gives,

$\displaystyle\dot{V}(z(t),\dot{z}(t))$
$\displaystyle=m\dot{z}(t)\ddot{z}(t)+k(z(t))\dot{z}(t)$
$\displaystyle=\dot{z}(t)u(t)-\beta(z(t),\dot{z}(t))\dot{z}^{2}(t)\leq\dot{z}(t)u(t),$

which implies that the spring-damper-mass system given in (67) is indeed NI
system.

Now, by defining the states as ${x}_{1}(t)=z(t)$ and ${x}_{2}(t)=\dot{z}(t)$,
the spring-damper-mass system, given in (67), can be written as follows;

$\displaystyle\left[\begin{array}[]{l}\dot{x}_{1}(t)\\\
\dot{x}_{2}(t)\end{array}\right]=\left[\begin{array}[]{c}x_{2}(t)\\\
\frac{-K(x_{1}(t))}{m}-\frac{\beta(x_{1}(t),x_{2}(t))}{m}x_{2}(t)+\frac{u(t)}{m}\end{array}\right],$
(72) $\displaystyle
y(t)=\left[\begin{array}[]{ll}1&0\end{array}\right]\left[\begin{array}[]{l}x_{1}(t)\\\
x_{2}(t)\end{array}\right].$ (76)

Next, the optimization problem in (2) is used to learn a linear system that
represents the nonlinear dynamical system given in (72).

The lifting functions $\phi(x)$ is chosen as a class of radial basis functions
in the following form:

$\psi_{i}(x(j))=\left\|x(j)-r_{i}\right\|_{2}^{2}\ln\left\|x(j)-r_{i}\right\|_{2},$

where the center $r_{i}$ is selected from the uniform distribution randomly.

The new state space representation in terms of the lifting functions is given
as follows,

$\psi(x(j))=\left[x_{1}(j),x_{2}(j),\psi_{1}(x(j)),\dots,\psi_{N_{\textmd{rbf}}}(x(j))\right]^{\mathrm{T}}$
(77)

We assumed that $K(x_{1}(t))=\left(x_{1}(t)+x_{1}^{3}(t)\right)$ and the
nonlinearity in the damper is considered as
$\beta(x_{1}(t),x_{2}(t))=x_{1}^{2}(t)+x_{2}^{2}(t)$. Also, let the number of
the lifted states to be 6 states, i.e., $N_{\textmd{rbf}}=6$. The resulting
model is compared to the non-constrained Koopman model with the same number of
the states, the linearized model around the initial point
$x_{0}=[0,0]^{\mathrm{T}}$, and at the linearized model around the initial
point $x_{0}=[0.5,0.5]^{\mathrm{T}}$. A random forcing input signal was used
for generating data for the optimization problem.

Figure 3 presents the state $x_{2}$ evolutions for different models compared
with the true data generated from the nonlinear system. Also, Figure 4
presents the state $x_{1}$ evolutions for different models compared with the
true data generated from the nonlinear system. As it is shown in Figures 4 and
3, the linear models fail to approximate the nonlinear terms, whereas the
Koopman model was able accurately to approximate the nonlinearity.

Figure 3: The comparison of the state $x_{2}$ for different models compared
with the true data generated from the nonlinear system Figure 4: The
comparison of the state $x_{1}$ for different models compared with the true
data generated from the nonlinear system Figure 5: Bode plots for the
approximated linear system using Koopman operator. The blue lines represent
the NI constrained model developed in this paper. While the orange lines
represents the unconstrained model.

Figure 5 displays the Bode plots corresponding to the approximated linear
system employing the Koopman operator. The blue lines depict the model
constrained by NI that has been developed in this manuscript, whereas the
unconstrained model is represented by the orange lines. It is clear from the
phase plot that the NI constrained model satisfy the NI property, i.e., the
phase between $(0,-\pi)$, however, the non-constrained model violates the NI
property.

The mean squared error (MSE) between the true values of the states and the
generated constrained model is presented in the Table I. The table shows that
the raise in the MSE in the constrained model is very small.

TABLE I: The mean squared error for different models. MSE between: | $x_{1}$ | $x_{2}$
---|---|---
True and the constrained model | $0.0078$ | $0.0042$
True and the unconstrained model | $0.0030$ | $0.0021$
Constrained and unconstrained model | $0.0095$ | $0.0059$

In order to emphasize on the importance of imposing the NI constraints on the
Koopman operator, we consider a feedback control design for the linearized
model. Suppose that a strictly NI controller in a form of a positive position
feedback controller (PPF) in (78), which is a common controller for such
systems, is designed for the spring-damper-mass.

$\displaystyle C(s)=\frac{K}{s^{2}+2\zeta\omega\ s+\omega^{2}}.$ (78)

When the controller given in (78) is connected in a positive feedback
connection with both, the linear NI constrained model and the unconstrained
model, a stable closed-loop is obtained in the case of the NI constrained
model while in the case of the unconstrained model, the closed-loop system is
unstable. This is shown in the step response of the closed-loop in Figure 6
and Nyquist plot in Figure 7. It indeed shows the importance of obtaining a NI
linear model approximation of a NI nonlinear system for controller synthesis.

Figure 6: A positive position feedback controller connection with the
nonlinear model (dashed line), the linear NI constrained model (blue solid
line) and the unconstrained model Figure 7: Nyquist plots for the closed loop.

##  V Conclusion

In conclusion, this paper proposes a data-driven approach to learning a lifted
linear Negative Imaginary (NI) dynamics that approximate nonlinear dynamical
systems using Koopman theory. The proposed method embeds the linear matrix
inequality that characterizes the NI property in the Koopman framework,
resulting in a non-convex optimization problem. To overcome the numerical
challenges of solving the problem, the optimization variables are reformatted
to convert the problem into a convex one. The comparison of the proposed
method with local linearization techniques shows that the proposed method can
accurately capture the nonlinear dynamics and achieve better control
performance. This method provides a numerically tractable solution for
learning the Koopman operator under NI constraints for nonlinear NI systems
and opens up new possibilities for applying linear control techniques to
nonlinear systems without linearization approximations. Overall, the proposed
method is a significant contribution to the field of control engineering, as
it provides a promising approach for controlling nonlinear NI systems in
various engineering applications.

## References

  * [1] A. Lanzon and I. R. Petersen, “Stability robustness of a feedback interconnection of systems with negative imaginary frequency response,” _IEEE Transactions on Automatic Control_ , vol. 53, no. 4, pp. 1042–1046, 2008.
  * [2] M. A. Mabrok, A. G. Kallapur, I. R. Petersen, and A. Lanzon, “Generalizing negative imaginary systems theory to include free body dynamics: Control of highly resonant structures with free body motion,” _IEEE Transactions on Automatic Control_ , vol. 59, no. 10, pp. 2692–2707, 2014.
  * [3] I. R. Petersen, “Negative imaginary systems theory and applications,” _Annual Reviews in Control_ , vol. 42, pp. 309–318, 2016.
  * [4] V. P. Tran, F. Santoso, M. Garratt, and S. Anavatti, “Neural network-based self-learning of an adaptive strictly negative imaginary tracking controller for a quadrotor transporting a cable-suspended payload with minimum swing,” _IEEE Transactions on Industrial Electronics_ , vol. 68, no. 10, pp. 10 258–10 268, 2020.
  * [5] A. Preumont, _Vibration Control of Active Structures: An Introduction_ , N. Jacobs, Ed. Springer, 2011.
  * [6] M. A. Mabrok, A. G. Kallapur, I. R. Petersen, and A. Lanzon, “Stabilization of conditional uncertain negative-imaginary systems using riccati equation approach,” in _Proceedings of the 20th International Symposium on Mathematical Theory of Networks and Systems (MTNS)_ , 2012, pp. 9–13.
  * [7] M. A. Mabrok and I. R. Petersen, “Controller synthesis for negative imaginary systems: a data driven approach,” _IET Control Theory & Applications_, vol. 10, no. 12, pp. 1480–1486, 2016.
  * [8] M. Mabrok, “Controller synthesis for negative imaginary systems using nonlinear optimisation and h 2 performance measure,” _International Journal of Control_ , vol. 94, no. 3, pp. 579–587, 2021.
  * [9] M. A. Mabrok, M. A. Haggag, and I. R. Petersen, “System identification algorithm for negative imaginary systems,” _Appl. Comput. Math_ , vol. 14, no. 3, pp. 336–348, 2015.
  * [10] J. Xiong, J. Lam, and I. R. Petersen, “Output feedback negative imaginary synthesis under structural constraints,” _Automatica_ , vol. 71, pp. 222–228, 2016.
  * [11] A. Ferrante, A. Lanzon, and L. Ntogramatzidis, “Foundations of negative imaginary systems theory and relations with positive real systems,” _arXiv preprint arXiv:1412.5709_ , 2014.
  * [12] ——, “Discrete-time negative imaginary systems,” _Automatica_ , vol. 79, pp. 1–10, 2017.
  * [13] M. Liu and J. Xiong, “Properties and stability analysis of discrete-time negative imaginary systems,” _Automatica_ , vol. 83, pp. 58–64, 2017.
  * [14] I. R. Petersen, A. Lanzon, and Z. Song, “Stabilization of uncertain negative-imaginary systems via state-feedback control,” in _Proceedings of the European Control Conference_ , Budapest, Hungary, August 2009.
  * [15] A. G. Ghallab and I. R. Petersen, “Negative imaginary systems theory for nonlinear systems: A dissipativity approach,” _arXiv preprint arXiv:2201.00144_ , 2022.
  * [16] A. Mauroy, Y. Susuki, and I. Mezić, _Koopman operator in systems and control_. Springer, 2020.
  * [17] M. Korda and I. Mezić, “Linear predictors for nonlinear dynamical systems: Koopman operator meets model predictive control,” _Automatica_ , vol. 93, pp. 149–160, 2018.
  * [18] P. Bevanda, S. Sosnowski, and S. Hirche, “Koopman operator dynamical models: Learning, analysis and control,” _Annual Reviews in Control_ , vol. 52, pp. 197–212, 2021.
  * [19] J. Xiong, I. R. Petersen, and A. Lanzon, “A negative imaginary lemma and the stability of interconnections of linear negative imaginary systems,” _IEEE Transactions on Automatic Control_ , vol. 55, no. 10, pp. 2342–2347, 2010.
  * [20] B. D. O. Anderson and S. Vongpanitlerd, _Network Analysis and Synthesis: A Modern Systems Approach_. Englewood Cliffs, N.J., USA: Prentice-Hall, 1973.
  * [21] B. Brogliato, R. Lozano, B. Maschke, and O. Egeland, _Dissipative Systems Analysis and Control_ , 2nd ed., ser. Communications and Control Engineering. London, UK: Springer, 2007\.
  * [22] A. Lanzon and I. R. Petersen, “A modified positive-real type stability condition,” in _Proceedings of the European Control Conference_ , Kos, Greece, Jul. 2007, pp. 3912–3918.
  * [23] A. G. Ghallab, M. A. Mabrok, and I. R. Petersen, “Extending negative imaginary systems theory to nonlinear systems,” in _2018 IEEE Conference on Decision and Control (CDC)_. IEEE, 2018, pp. 2348–2353.
  * [24] M. A. Mabrok, M. A. Alyami, and E. E. Mahmoud, “On the dissipativity property of negative imaginary systems,” _Alexandria Engineering Journal_ , vol. 60, no. 1, pp. 1403–1410, 2021.
  * [25] I. R. Petersen, “Physical interpretations of negative imaginary systems theory,” in _2015 10th Asian Control Conference (ASCC)_ , May 2015, pp. 1–6.
  * [26] K. Hara, M. Inoue, and N. Sebe, “Learning koopman operator under dissipativity constraints,” _IFAC-PapersOnLine_ , vol. 53, no. 2, pp. 1169–1174, 2020\.


# Sums of Consecutive Prime Squares

Janyarak Tongsomporn, Saeree Wananiyakul, Jörn Steuding

(Date: January 2021)

###### Abstract.

We prove explicit bounds for the number of sums of consecutive prime squares
below a given magnitude

Keywords: prime numbers, sums of squares

MSC Numbers: 11A41, 00A08

##  1\. Motivation and the Main Result

Early last year the authors learned that 2020 can be represented as a sum of
squares of consecutive prime numbers, namely

$2020=17^{2}+19^{2}+23^{2}+29^{2}.$

It is a natural question to ask what the next year with this property will be.
We shall show that such a representation is a rare event.

Indeed, if ${\rm{scp}}(x)$ counts the number of sums of squares of consecutive
primes below $x$, i.e.,

${\rm{scp}}(x)=\sharp\left\\{p_{n}^{2}+p_{n+1}^{2}+\ldots+p_{n-1+m}^{2}\leq
x\,:\,m\in\hbox{{\dubl N}}\right\\},$

where $p_{j}$ denotes the $j$-th prime number in ascending order, then
$\lim_{x\to\infty}{\rm{scp}}(x)/x=0$. The following theorem provides more
precise bounds.

######  Theorem 1.

We have

$2\,{x^{1/2}\over\log
x}<\pi(\sqrt{x})\leq{\rm{scp}}(x)<10.9558\,{x^{2/3}\over(\log x)^{4/3}}\,,$

where the inequality on the far left is valid for $x\geq 289$ and all those to
the right for $x>1$.

Here, as usual, $\pi(N)$ is counting the number of primes $p\leq N$ and
explicit bounds for this prime counting function are the main tool for proving
the inequalities above; we have chosen a recent paper [1] by Pierre Dusart.
The dear reader is invited to improve upon the bounds of our theorem; maybe it
is even possible to prove an asymptotic formula for the number of sums of
consecutive prime squares below a given magnitude. Note that we do not
consider here the question whether or not an integer can have two or even more
such representations or how many of these exist.

Using a computer algebra package one can verify that the next sum of squares
of consecutive primes is given by the expected suspect, namely

$2189=13^{2}+17^{2}+19^{2}+23^{2}+29^{2}.$

A list with all integers below $5000$ that can be written as a sum of
consecutive prime squares can be found in the third and final section.

This year’s prime factorization is $2021=43\cdot 47$ which is a product of two
consecutive primes. Following our approach one can also discuss sums of
products of two consecutive primes; the corresponding bounds should be close
to those found here for sums of consecutive prime squares.

##  2\. Proof of the Theorem

It is convenient to define, for fixed $m\in\hbox{{\dubl N}}$, the counting
function for sums of $m$ consecutive prime squares, i.e.,

${\rm{scp}}_{m}(x)=\sharp\left\\{p_{n}^{2}+p_{n+1}^{2}+\ldots+p_{n-1+m}^{2}\leq
x\right\\}.$

For the lower bound we first observe that the number of squares of prime
numbers $p^{2}$ below or equal to $x$ is given by $\pi(\sqrt{x})$.

In the sequel we shall use the explicit bounds

(1) ${N\over\log N}<\pi(N)<1.2551\,{N\over\log N},$

where the left inequality is valid for $N\geq 17$ and the one on the right for
$N>1$ (see Corollary 5.2 in [1]); of course, the celebrated prime number
theorem provides an asymptotic formula for $\pi(N)$ with main term $N/\log N$,
however, for excluding the related error term for our analysis, we prefer the
version above with the factor $1.2551$. The corresponding range for these
inequalities (resp. the range for $x$ in our theorem) is also useful with
respect to computer experiments.

It follows from (1) that

${\rm{scp}}(x)\geq{\rm{scp}}_{1}(x)=\pi(\sqrt{x})>{\sqrt{x}\over{1\over 2}\log
x},$

which is valid for $x\geq 17^{2}=289$. This proves the lower bound.

The reasoning for the upper bound is a little more advanced. First we note
that for $n={\rm{scp}}(x)$ we have

$mp_{n}^{2}\leq p_{n}^{2}+p_{n+1}^{2}+\ldots+p_{n-1+m}^{2}\leq
x<p_{n}^{2}+p_{n+1}^{2}+\ldots+p_{n-1+m}^{2}+p_{n+m}^{2}.$

Hence, by the inequality in (1),

(2) ${\rm{scp}}_{m}(x)=n\leq\pi(\sqrt{x/m})<1.2551\,{\sqrt{x/m}\over{1\over
2}\log(x/m)}\leq 2.5102\,{(x/m)^{1/2}\over\log x},$

which is valid for $x>m$, which, obviously, is no severe restriction (since
the largest integer $\leq x$ is a trivial upper bound for the length of a sum
of consecutive primes squares $\leq x$).

To continue we shall next bound the length of possible sums of consecutive
prime squares below $x$. For this purpose we shall use an old result due to
Barkley Rosser [2] which has been improved several times, in particular by
Dusart [1], however, we prefer the more simple inequality

$p_{n}>n\log n,$

valid for all $n\in\hbox{{\dubl N}}$; this lower bound is trivial for $n=1$.
We shall use this so-called Rosser theorem for the sum of the squares of the
first primes:

$p_{1}^{2}+p_{2}^{2}+\ldots+p_{M}^{2}>\sum_{2\leq n\leq M}(n\log n)^{2}.$

If we can show that the right hand side is larger than $x$, then the least sum
of $M$ consecutive prime squares already exceeds the given magnitude. Assuming
that this $M$ is the least positive integer with this property, this leads to
a bound for $M$ depending on $x$. This estimate in combination with the
previous one allows us to derive the upper bound of the theorem.
Alternatively, one could also use partial summation here together with the
prime number theorem, however, it is our intention to circumvent error terms.

Obviously, for $M\geq 4$,

$\displaystyle\sum_{2\leq n\leq M}(n\log n)^{2}$ $\displaystyle\geq$
$\displaystyle\sum_{\sqrt{M}\leq n\leq M}n^{2}(\log n)^{2}$
$\displaystyle\geq$ $\displaystyle({\textstyle{1\over 2}}\log
M)^{2}\sum_{\sqrt{M}\leq n\leq M}n^{2}\geq{\textstyle{1\over 12}}M^{3}(\log
M)^{2},$

where we have used in the final step the well-known formula

$1+2^{2}+3^{2}+\ldots+M^{2}={\textstyle{1\over 6}}\,M(M+1)(2M+1)$

and some pen and paper. It thus follows that every sum of consecutive prime
squares below $x$ has less than roughly $x^{1/3}$ summands. For a more precise
bound we observe that substituting

(3) $M=\lfloor 108^{1/3}\,x^{1/3}(\log x)^{-2/3}\rfloor$

into the lower bound above yields a quantity slightly larger than $x$; here
$\lfloor z\rfloor$ denotes the largest integer $\leq z$.

To use this for an upper bound we first observe that (2) implies

${\rm{scp}}(x)=\sum_{1\leq m\leq M}{\rm{scp}}_{m}(x)<2.5102\,{x^{1/2}\over\log
x}\sum_{1\leq m\leq M}m^{-1/2}.$

In general, we have, for $\alpha\in(0,1)$,

$\sum_{1\leq m\leq M}m^{-\alpha}<1+\sum_{2\leq m\leq
M}\int_{m-1}^{m}u^{-\alpha}{\rm{d}}u=1+\int_{1}^{M}u^{-\alpha}{\rm{d}}u={M^{1-\alpha}-\alpha\over
1-\alpha}.$

This in combination with (3) leads to

${\rm{scp}}(x)<5.0204\,{(xM)^{1/2}\over\log
x}\leq{5.0204\cdot(108)^{1/6}}\,{x^{2/3}\over(\log x)^{4/3}}.$

This proves the upper bound of the theorem.

##  3\. Explicit Sums of Consecutive Prime Squares

We conclude with a list of all integers below $5000$ that can be written as a
sum of consecutive prime squares:

4 |  9 |  13 |  25 |  34 |  38 |  49
---|---|---|---|---|---|---
74 |  83 |  87 |  121 |  169 |  170 |  195
204 |  208 |  289 |  290 |  339 |  361 |  364
373 |  377 |  458 |  529 |  579 |  628 |  650
653 |  662 |  666 |  819 |  841 |  890 |  940
961 |  989 |  1014 |  1023 |  1027 |  1179 |  1348
1369 |  1370 |  1469 |  1518 |  1543 |  1552 |  1556
1681 |  1731 |  1802 |  1849 |  2020 |  2189 |  2209
2310 |  2330 |  2331 |  2359 |  2384 |  2393 |  2397
2692 |  2809 |  2981 |  3050 |  3150 |  3171 |  3271
3320 |  3345 |  3354 |  3358 |  3481 |  3530 |  3700
3721 |  4011 |  4058 |  4061 |  4350 |  4489 |  4519
4640 |  4689 |  4714 |  4723 |  4727 |  4852 |  4899

## References

  * [1] P. Dusart, Explicit estimates of some functions over primes, Ramanujan J. 45 (2018), 227–251
  * [2] B. Rosser, The $n$-th prime is greater than $n\log n$, Proc. London math. Soc. (2) 45 (1938), 21-44

Janyarak Tongsomporn, Saeree Wananiyakul, ${\mathcal{W}}$alailak University,
School of Science, Nakhon Si Thammarat 80 160, Thailand<EMAIL_ADDRESS>
Jörn Steuding, Department of Mathematics, ${\mathcal{W}}$ürzburg University,
Am Hubland, 97 218 Würzburg, Germany<EMAIL_ADDRESS>