# On Linear Time Invariant Systems Analysis via A Single Trajectory: A Linear
Programming Approach

Hassan Abdelraouf1, Fahad Albalawi2 and Eric Feron3 *This work was supported
by King Abdullah University of Science and Technology (KAUST)1Hassan
Abdelraouf is a Ph.D. student with Mechanical Engineering, (KAUST), Thuwal,
KSA<EMAIL_ADDRESS>Fahad Albalawi is a research scientist
with Electrical and Computer Engineering, KAUST, Thuwal, KSA
<EMAIL_ADDRESS>Eric Feron is with the Faculty of Electrical,
Computer, and Mechanical Engineering at KAUST, Thuwal, KSA
<EMAIL_ADDRESS>
###### Abstract

In this note, a novel methodology that can extract a number of analysis
results for linear time-invariant systems (LTI) given only a single trajectory
of the considered system is proposed. The superiority of the proposed
technique relies on the fact that it provides an automatic and formal way to
obtain a valuable information about the controlled system by only having
access to a single trajectory over a finite period of time (i.e., the system
dynamics is assumed to be unknown). At first, we characterize the stability
region of LTI systems given only a single trajectory dataset by constructing
the associated Lyapunov function of the system. The Lyapunov function is found
by formulating and solving a linear programming (LP) problem. Then, we extend
the same methodology to a variety of essential analysis results for LTI
systems such as deriving bounds on the output energy, deriving bounds on
output peak, deriving $\mathbf{L}_{2}$ and RMS gains. To illustrate the
efficacy of the proposed data-driven paradigm, a comparison analysis between
the learned LTI system metrics and the true ones is provided.

##  I INTRODUCTION

Data-driven based control schemes gained a significant attention in the
control society over the last decade. In fact, such control paradigms become
an attractive alternative to conventional control algorithms [1, 2]. This
paradigm shift from explicit control techniques to learning-based ones has
produced a massive amount of theoretical and practical research works where
the majority of these control techniques were obtained from supervised machine
learning and reinforcement learning [3, 4, 5, 6, 7]. Despite the fact that
theoretical certificates for control performance may be derived in various
ways, the direct relationship between the generated or collected data and the
learning-based control performance is not well-established [8]. Albeit data
quality metrics such as entropy have been utilized substantially to guide
exploration and control strategies [9, 10], they still do not yield direct
conclusions into the influence of data on the provable control performance.

As a remedy, [11] utilized Gaussian process priors to generate Lyapunov-based
measure that can assess the value of data points with respect to a number of
control tasks. The ultimate goal of collecting dataset and assessing their
quality for dynamical system is to derive stability and performance
certificates that can be used later for different control tasks. As customary,
the stability of an equilibrium point of a given dynamic system can be studied
through the Lyapunov function. Lyapunov functions can be seen as stability
certificates for ordinary differential equations (ODEs) [12, 13]. The problem
of finding the Lyapunov function is generally complex, and it has been the
scope of many research papers in the control community [14]. Analytical
techniques are standard tools for constructing Lyapunov functions. Despite the
fact that these techniques are mathematically stable and sound, they require
significant expertise and manual efforts [15]. For LTI systems, which is the
focus class of systems for this work, semi-definite programming is sufficient
to construct Lyapunov functions due to the fact that Lyapunov functions are
inherently quadratic polynomials for LTI systems. However, the knowledge of
the underlying dynamical linear systems, i.e., the matrices (A,B,C,D), is
always assumed to be known in order to derive the Lyapunov function via semi-
definite programming approaches. Such an assumption might not be fulfilled
when first principle models cannot be derived because of the complexity of the
actual system. Alternatively, data-driven algorithms can provide an
approximate mathematical model from real system measurements.

Motivated by the aforementioned observations, we introduce a formal and
automatic methodology for LTI systems analysis where the model of the
controlled system is not available and only a single trajectory dataset is
provided. Specifically, we formulate the problems of finding the Lyapunov
function and observability gramium of LTI systems as a linear programming (LP)
problem given dataset that represents a single system trajectory over a finite
period of time. Many problems in control systems analysis like calculating
system’s RMS gain and output peak are formulated as Semi-definite Programming
(SP) and Linear Matrix Inequalities (LMIs) [16] where a full knowledge of the
system dynamics is assumed. By using our novel approach, these problems can be
solved accurately and efficiently without knowing the system dynamics, and
only a single trajectory over a finite period of time is given. Majority of
the proposed data-driven based Lyapunov function methods in the literature
focus on using Sum of Squares (SOS) methods [17, 18] as well as machine
learning techniques like Neural networks. [19, 20]. Nevertheless, all of these
methods are model-based and computationally expensive where big data sets are
needed for efficient learning. Our approach requires the least data points
along the trajectory to learn the Lyapunov function for LTI systems as well as
other LTI system analysis metrics. Finally, we test our proposed frameworks by
comparing them with the true Lyapunov function as well as the true LTI systems
metrics where superior performance of our proposed algorithms is well-
established.

##  II Preliminaries

###  II-A Class of LTI systems

We consider a class of LTI systems that can be written in the following state
description:

$\displaystyle\begin{split}\dot{x}&=Ax+Bu\\\ z&=Cx+Du\end{split}$ (1)

where $x\in\mathbb{R}^{n}$ represents the state vector, $u\in\mathbb{R}^{m}$
is the control input vector, and $z\in\mathbb{R}^{p}$ is the output state
vector.

###  II-B Notations

The $2$ norm of a vector $x\in\mathbb{R}^{n}$ is defined as follows
$\left(\sum_{i=1}^{n}\left|x_{i}\right|^{2}\right)^{1/2}$. We denote the
vector of unique elements of a symmetric matrix $P\in\mathbb{R}^{n\times n}$
as $\text{vec}(P)$ and it is defined as follows:

$\text{vec}(P)=\begin{bmatrix}[P_{11}\quad\dots\quad P_{nn}]^{T}\\\
[P_{12}\quad\dots\quad P_{1n}]^{T}\\\ [P_{23}\quad\dots\quad P_{2n}]^{T}\\\
\vdots\\\ P_{(n-1)n}\end{bmatrix}$ (2)

where $\text{vec}(P)\in\mathbb{R}^{n(n+1)/2}$. For instance, when
$P\in\mathbb{R}^{2\times 2}$ is symmetric, then $\text{vec}(P)=[P_{11}\quad
P_{22}\quad P_{12}]^{T}$. In addition, for the two vectors
$x,y\in\mathbb{R}^{n}$ , we define $\oplus$ operator as follows:

$x\oplus y=\begin{bmatrix}[x_{1}y_{1}\quad\dots\quad x_{n}y_{n}]^{T}\\\
[x_{1}y_{2}+x_{2}y_{1}\quad\dots\quad x_{1}y_{n}+x_{n}y_{1}]^{T}\\\
[x_{2}y_{3}+x_{3}y_{2}\quad\dots\quad x_{2}y_{n}+x_{n}y_{2}]^{T}\\\ \vdots\\\
x_{n-1}y_{n}+x_{n}y_{n-1}\end{bmatrix}$ (3)

If $z=x\oplus y$, then $z\in\mathbb{R}^{n(n+1)/2}$. For example, if
$x,y\in\mathbb{R}^{2}$, then $x\oplus x=[x_{1}^{2}\quad x_{2}^{2}\quad
2x_{1}x_{2}]^{T}$ and $x\oplus y=[x_{1}y_{1}\quad x_{2}y_{2}\quad
x_{1}y_{2}+x_{2}y_{1}]^{T}$. Using the definitions of both operators $\oplus$
and $\text{vec}(.)$, we can represent the term $x^{T}Py$ as $(x\oplus
y)^{T}\text{vec}(P)$

Finally, the $\mathbf{L}_{2}$ norm of a signal $\zeta$ is defined as
$\left\lVert\zeta\right\rVert_{2}^{2}=\int_{0}^{\infty}\zeta^{T}\zeta\text{dt}$
and its root-mean square gain $\mathbf{RMS}(\xi)$ is defined as follows:

$\mathbf{RMS}(\xi)\triangleq\left(\limsup_{T\rightarrow\infty}\frac{1}{T}\int_{0}^{T}\xi^{T}\xi
dt\right)^{1/2}$ (4)

###  II-C Outline

This paper is organized as follows. Section III presents using our approach to
learn lyapunov function and solve lyapunov equation for LTI system using data
along a single trajectory. Sections IV, V, VI show how the proposed approach
is used to learn bounds on output energy, bounds on output peak and
$\mathbf{L}_{2}$ gain of LTI systems respectively.

##  III Data-driven Construction for Lyapuonv function for LTI systems

Stability analysis of dynamical systems aims to show that a set of initial
states will stay in the neighborhood of an equilibrium point or converge to
it. Based on the construction of Lyapunov functions, the stability of
equilibrium points, which can form positive invariant sets (i.e., regions of
attraction), can be certified. Lyapunov showed that the linear time invariant
system is stable (i.e. all trajectories converge to zero asymptotically) if
and only if there exists a quadratic function $V(x)=x^{T}Px$ which is positive
$V(x)>0$ and its gradient is negative along the system’s trajectories
$\dot{V}(x)<0$. These two conditions can be formulated as the following Linear
matrix inequalities: $P>0$ and $PA+A^{T}P<0$ [16]. Lyapunov inequality in $P$
can be explicitly solved by picking $Q=Q^{T}>0$ , then solving the linear
equation $PA+A^{T}P=-Q$ for the matrix $P$ [16]. This method assumes the
knowledge of the system dynamic equation of Eq. 1. For real systems, it may be
very difficult to obtain the system dynamics and very expensive to estimate
its unknown parameters. Hence, we propose a new method to learn Lyapunov
function for LTI systems given only a single trajectory. To the best of our
knowledge, our proposed approach is very novel and it has never been
introduced in the control society as a tool for LTI system analysis. For
stable LTI systems, Lyapunov function is known to be quadratic ($x^{T}Px$).
So, our approach focuses on deriving a $P$ matrix that makes the Lyapunov
function positive and its gradient negative at all given data points along the
given system trajectory.

First, a trajectory for the LTI system is generated from the initial time
$t_{0}=0$ to $T+\text{dt}$ with a time step dt. The states along this
trajectory are $\begin{bmatrix}x(0)&x(1)&\dots&x(N)&x(N+1)\end{bmatrix}$,
where $N=T/\text{dt}$ and $x(i)$ represents the state vector at time
$i\text{dt}$. Then, the state vector derivative w.r.t. time $\dot{x}(i)$ are
calculated numerically using forward finite difference differentiation from
$i=0$ to $i=N$. So, the data points $i=0\dots N$ for state vectors $x$ and
state vectors derivative $\dot{x}$ are only used to obtain the unknown
parameters of $P$ matrix. The objective is to find the unknown parameters of
$P\in\mathbb{R}^{n\times n}$ which satisfies Lyapunov conditions at every
point on the given trajectory:

$\begin{gathered}V(x(i))=x(i)^{T}Px(i)>0\\\
\frac{d}{\text{dt}}V(x(i))=2x(i)^{T}P\dot{x}(i)<0\end{gathered}$ (5)

for all $i=0,\dots,N$. It is known that $P$ matrix is a symmetric positive
definite matrix, so the number of unknown parameters in $P$ is $n(n+1)/2$. Let
$p\in\mathbb{R}^{n(n+1)/2}$ is the vector of unknown parameters. So,
$p=\text{vec}(P)$ using (2). Therefore, the Lyapunov-based conditions of (5)
can be formulated as a set of linear inequalities in $p$ as: $L_{1}p>\epsilon$
and $L_{2}p<-\epsilon$ where,

$L_{1}=\begin{bmatrix}(x(0)\oplus x(0))^{T}\\\ (x(1)\oplus x(1))^{T}\\\
\vdots\\\ (x(N)\oplus x(N))^{T}\end{bmatrix}$ (6)
$L_{2}=2\begin{bmatrix}(x(0)\oplus\dot{x}(0))^{T}\\\
(x(1)\oplus\dot{x}(1))^{T}\\\ \vdots\\\
(x(N)\oplus\dot{x}(N))^{T}\end{bmatrix}$ (7)

such that $L_{1},L_{2}\in\mathbb{R}^{(N+1)\times n(n+1)/2}$. The parameter
$\epsilon$ is a positive small number. These linear inequalities can be solved
trivially by any linear programming solver such as CVX [21] to get $p$ (i.e.,
the vector of unknown parameters of $P$ matrix that defines our desired
Lyapunov function $V(x)$).

Numerical example: We have the following unforced LTI system: $\dot{x}=Ax$
where

$A=\begin{bmatrix}0&1\\\ -1&-3\end{bmatrix}$ (8)

Our approach utilizes only one generated trajectory and then learns the
quadratic Lyapunov function $x^{T}Px$ from data along this trajectory to prove
the system stability. The given trajectory starts at $x(0)=[2,2]^{T}$ to the
final time $T=1\text{sec}$ with time step $dt=0.01\text{sec}$ and $N=100$.
First, the finite difference method is used to get $\dot{x}$ along the given
trajectory. Then, the two sets of linear inequalities that represent the
satisfaction of Lyapunov conditions at every point on the given trajectory
$L_{1}p>\epsilon$ and $L_{2}p<-\epsilon$ are solved. Since the dimension of
the state vector $x$ is $\mathbb{R}^{2}$, $L_{1}$ and $L_{2}$ can be
structured as follows:

$L_{1}=\begin{bmatrix}x_{1}(0)^{2}&x_{2}(0)^{2}&2x_{1}(0)x_{2}(0)\\\
\vdots&\vdots&\vdots\\\
x_{1}(N)^{2}&x_{2}(N)^{2}&2x_{1}(N)x_{2}(N)\end{bmatrix}$ (9)
$L_{2}=2\begin{bmatrix}x_{1}(0)\dot{x_{1}}(0)&x_{2}(0)\dot{x_{2}}(0)&x_{1}(0)\dot{x_{2}}(0)+x_{2}(0)\dot{x_{1}}(0)\\\
\vdots&\vdots&\vdots\\\
x_{1}(N)\dot{x_{1}}(N)&x_{2}(N)\dot{x_{2}}(N)&x_{1}(N)\dot{x_{2}}(N)+x_{2}(N)\dot{x_{1}}(N)\end{bmatrix}$
(10)

These linear inequalities are solved by Matlab CVX solver [21] to get the
unknown $P$ matrix:

$P=\begin{bmatrix}26.4840&5.3151\\\ 5.3151&17.0361\end{bmatrix}$ (11)

Since $P$ is positive definite matrix and $PA+A^{T}P$ is negative definite,
the learned quadratic function $V(x)=x^{T}Px$ is a valid Lyapunov function for
the LTI system $\dot{x}=Ax$. Such function was derived from dataset along a
single trajectory of the considered LTI system. Fig. 1 shows one level set of
the learned Lyapunov function ($x^{T}Px<=1000$) which constitutes a forward
invariant set $\Omega_{\rho}$ where $\rho$ is 1000. Several points at the
boundary of the Lyapunov level set are chosen to be initial states for the
system trajectories. All the trajectories remain inside the set which
demonstrates that the learned function is a true Lyapunov function.

Figure 1: Data-driven Lyapunov function with $V(x)<=1000$

Lyapunov functions for stable LTI systems are not necessarily unique. As a
result of such fact, changes on the Lyapunov conditions can be made by the
user to improve the numerical stability of the linear program. Lyapunov
conditions can be modified to be $L_{1}p>c_{1}l$ and $L_{2}p<-c_{2}l$. Where
$c_{1},c_{2}>0$ and $l\in\mathbb{R}^{(N+1)\times 1}$ is:

$l=\begin{bmatrix}\left\lVert x(0)\right\rVert_{2}^{2}\\\ \vdots\\\
\left\lVert x(N)\right\rVert_{2}^{2}\end{bmatrix}$ (12)

To improve the robustness of the learned Lyapunov function $V(x)$, the number
of data points need to be increased and a probing noise to the given
trajectory needs to be augmented.

###  III-A Exact solution of Lyapunov equation from data

The proposed approach can be employed to solve the Lyapunov equation
$PA+A^{T}P=-Q$ accurately by knowing only the states $x$ and their derivatives
w.r.t time evaluated at least at $n(n+1)/2$ points on any given trajectory.
where $n$ is the state space dimension. The solution of the Lyapunov equation
is based on formulating the problem as an LP problem without the need to know
the system dynamics (i.e., $A$ matrix for LTI systems). The matrix
$Q\in\mathbb{R}^{n\times n}$ at the right hand side of the Lyapunov function
is a user defined symmetric positive definite matrix. Hence, the set of linear
inequalities will be: $L_{1}p>\epsilon$, $L_{2}p\leq l_{Q}$ and $L_{2}p\geq
l_{Q}$. The vectors $L_{1}$ and $L_{2}$ are defined in (6) and (7) while
$l_{Q}$ is defined as follows:

$l_{Q}=\begin{bmatrix}x(0)^{T}Qx(0)\\\ x(1)^{T}Qx(1)\\\ \vdots\\\
x(N)^{T}Qx(N)\end{bmatrix}$ (13)

Numerical example: From the same state trajectory used in the previous
example, we will now pick only three data points to derive the Lyapunov
equation of the system of Eq. 8. The three data points are ($x,\dot{x}$) at
the following time instants $0,0.5$ and $1$ sec. When $Q=-I$, then the sets of
linear inequalities in $p$ (i.e., $L_{1}p>\epsilon$, $L_{2}p\leq l_{Q}$ and
$L_{2}p\geq l_{Q}$) are as follows:

$\displaystyle L_{1}$ $\displaystyle=\begin{bmatrix}4&4&8\\\
4.9627&0.4626&-3.0304\\\ 2.6757&2.1311&-4.7759\end{bmatrix}$ (14)
$\displaystyle L_{2}$ $\displaystyle=\begin{bmatrix}8&-32&-24\\\
-3.0304&0.2547&0.0910\\\ -4.7759&-8.0107&13.2384\end{bmatrix}$ (15)
$\displaystyle l_{Q}$ $\displaystyle=\begin{bmatrix}8\\\ 5.4253\\\
4.8068\end{bmatrix}$ (16)

The solution of these linear inequalities is:

$P=\begin{bmatrix}1.8333&0.5000\\\ 0.5000&0.3333\end{bmatrix}$ (17)

which is exactly the same solution of $PA+A^{T}P=-I$, yet derived by only
using three data points along any given system’s trajectory. Such an
interesting result shows the efficacy and accuracy of our proposed approach.

The main idea is that $P$ matrix has $n(n+1)/2$ unknown parameters so, we need
only $n(n+1)/2$ equations to solve for $P$. These equations represents the
satisfaction of Lyapunov second inequality at $n(n+1)/2$ points on the given
trajectory.

##  IV Data-driven evaluation for bounds on the output energy

The maximum output energy for LTI system given the initial state is given by:

$\max\left\\{\int_{0}^{\infty}z^{T}zdt\mid\dot{x}=Ax,\quad z=Cx\right\\}$ (18)

Where the initial state $x(0)$ is given. Suppose that there exists a quadratic
positive function $V(\zeta)=\zeta^{T}P\zeta$ such that:

$P>0\text{ and }\frac{d}{\text{dt}}V(x)\leq-z^{T}z,\quad\text{ for every
}x\text{ and }z$ (19)

By integrating both sides of the second inequality in (19) from $0$ to $T$, we
obtain,

$V(x(T))-V(x(0))\leq-\int_{0}^{T}z^{T}zdt$ (20)

Given that $V(x(T))\geq 0$, we conclude that $V(x(0))=x(0)^{T}Px(0)$ is an
upper bound for the maximum output energy given the initial condition $x(0)$
[16]. Assume ($A$, $C$) are given. Then, the second inequality can be
formulated as the following LMI: $PA+A^{T}P+C^{T}C\leq 0$. Therefore, we get
the best upper bound on the output energy by solving the following SDP in the
variable $P$:

$\displaystyle\underset{P}{\text{minimize}}$ $\displaystyle x(0)^{T}Px(0)$
(21) subject to $\displaystyle P>0$ $\displaystyle PA+A^{T}P+C^{T}C\leq 0$

This SDP can be solved analytically or using MATLAB CVX. The solution is
exactly equal to the output energy $x(0)^{T}W_{0}x(0)$, where $W_{0}$ is the
observability gramiam of the system:

$W_{\mathrm{o}}\triangleq\int_{0}^{\infty}e^{A^{T}t}C_{z}^{T}C_{z}e^{At}dt$
(22)

Assuming $(A,C)$ are unknown, the problem of finding the observability gramiam
or the maximum output energy given the initial state $x(0)$ of a LTI system
can be solved given only a single trajectory. The proposed method is based on
finding a quadratic function that satisfies conditions in (19) at every point
along the given trajectory. Therefore, the conditions in (19) can be
represented numerically as follows:

$\begin{gathered}x(i)^{T}Px(i)>0\\\
2x(i)^{T}P\dot{x}(i)\leq-z(i)^{T}z(i)\end{gathered}$ (23)

for $i=0,\dots,N$. Hence, instead of formulating the problem as an LMI problem
as in (21) given $(A,C)$, we can reformulate the problem as the following LP
formulation provided that one single trajectory over a finite period of time
is on hand:

$\displaystyle\underset{p}{\text{minimize}}$ $\displaystyle(x(0)\oplus
x(0))^{T}p$ (24) subject to $\displaystyle L_{1}p>0$ $\displaystyle
L_{2}p\leq-l_{z}$

where, $p=\text{vec}(P)$, $L_{1}$ and$L_{2}$ are defined in (6)and (7), and
$l_{z}$ is defined as:

$l_{z}=\begin{bmatrix}z(0)^{T}z(0)\\\ z(1)^{T}z(1)\\\ \vdots\\\
z(N)^{T}z(N)\end{bmatrix}$ (25)

Numerical example: Given the following LTI system: $\dot{x}=Ax$ , $z=Cx$

$A=\begin{bmatrix}0&1\\\ -4&-2\end{bmatrix}\quad\quad
C=\begin{bmatrix}0&1\end{bmatrix}$ (26)

A single trajectory at the initial state
$x(0)=\begin{bmatrix}2&2\end{bmatrix}^{T}$ over a finite time period
$T=5\text{sec}$ with a time step $dt=0.1$ is provided. The output $z$ and
state vector $x$ at each time step can be measured. The derivatives of the
state vector $\dot{x}$ can be calculated numerically. So, $N=50$. The LP
problem of (24) is solved and the following $P$ matrix was the solution:

$\begin{bmatrix}0.5000&0.1250\\\ 0.1250&0.0625\end{bmatrix}$ (27)

Which is exactly the same solution obtained by solving LMI (21) where the full
knowledge of the system dynamics is assumed. Therefore, from a single
trajectory of the LTI system, the observability gramium matrix can be derived,
then the maximum output energy of any given initial condition can be
determined by $x(0)^{T}Px(0)$.

##  V Data-driven Evaluation for bounds on output peak

The problem of deriving bounds on output energy $\left\lVert z(t)\right\rVert$
can be formulated as a SDP [16] when the system dynamics is known and the
initial state $x(0)$ is given. Let $\mathcal{E}=\left\\{\xi\mid\xi^{T}P\xi\leq
1\right\\}$ be an invariant ellipsoid that contains $x(0)$ for the LTI system
$\dot{x}=Ax$ , $z=cx$. Then,

$z(t)^{T}z(t)\leq\max_{\xi\in\mathcal{E}}\xi^{T}C^{T}C\xi$ (28)

In [16], the right hand side of (28) can be expressed as the square root of
the minimum $\delta$ subject to:

$\begin{bmatrix}P&C^{T}\\\ C&\delta I\end{bmatrix}\geq 0$ (29)

Therefore, given the initial state $x(0)$, the minimum bound of the output
peak is the square root of the optimal value of the following SDP where $P$
and $\delta$ are decision variables for the optimization problem.

$\displaystyle\underset{P,\delta}{\text{minimize}}$ $\displaystyle\delta$ (30)
subject to $\displaystyle P>0,\quad PA+A^{T}P\leq 0$ $\displaystyle
x(0)^{T}Px(0)\leq 1,\quad(\ref{bound outpeak cons})$

To obtain an optimal solution for the SDP of (30), the system dynamics has to
be known a prior. However, this problem can be solved without knowing the
system matrices ($A,C$) where only a single trajectory is assumed to be
available. Using our approach in sections III and IV, this problem can be
formulated as an LP problem given the initial state and $(x,\dot{x},z)$ at
every point along the given system trajectory. First, by Schur complement,
(29) can be rewritten as follows:

$P-C^{T}C/\delta\geq 0$ (31)

this condition should be satisfied along the given trajectory
$x(i)^{T}Px(i)-x(i)C^{T}Cx(i)/\delta\geq 0$. As a result, the constraints of
(30) will be represented as follows:

$\begin{gathered}x(i)^{T}Px(i)>0\\\ 2x(i)^{T}P\dot{x}(i)\leq 0\\\
x(0)^{T}Px(0)\leq 1\\\ x(i)^{T}Px(i)-z(i)^{T}z(i)/\delta\geq 0\end{gathered}$
(32)

The fourth constraint is nonlinear, so let $\lambda=1/\delta$ and instead of
minimizing $\delta$, $\lambda$ is maximized. Therefore, instead of formulating
the problem as a SDP (30), it can be formulated as the following LP problem
given data points along the system trajectory:

$\displaystyle\underset{p,\lambda}{\text{maximize}}$ $\displaystyle\lambda$
(33) subject to $\displaystyle L_{1}p>0,\quad L_{2}p\leq 0$
$\displaystyle(x(0)\oplus x(0))^{T}p\leq 1$ $\displaystyle L_{1}p-\lambda
l_{z}\geq 0$

where, $p=\text{vec}(P)$, $L_{1},L_{2},l_{z}$ are defined in (6),(7) and (25)
respectively. Given the initial state $x(0)$ and a single trajectory, the LP
problem of (33) can be solved to obtain the upper bound of the output peak
$\sqrt{1/\lambda}$ and the invariant ellipsoid $x^{T}Px\leq 1$ at which the
maximum output of any trajectory within the ellipsoid will not exceed that
upper bound.

Numerical example:

Using the same LTI system $\dot{x}=Ax,z=Cx$ such that $A$ and $C$ are defined
in (26). Given the initial state $x(0)=[3,3]$ and a trajectory of the LTI
system starting from $t=0~{}\text{sec}$ to $T=5~{}\text{sec}$ with a time step
$dt=0.1~{}\text{sec}$. Then the LP problem (33) is solved to obtain the upper
bound on the output peak and its corresponding invariant ellipsoid. The upper
bound found to be ($\sqrt{1/\lambda}=3.2901$) and the invariant ellipsoid is
$x^{T}Px\leq 1$ such that

$P=\begin{bmatrix}0.092453&0.001486\\\ 0.001486&0.015684\\\ \end{bmatrix}$
(34)

Fig. 2 shows the invariant ellipsoid $x^{T}Px\leq 1$ that contains the initial
state $x(0)$. As it can be seen from the figure, any trajectory starts at the
boundary or within the defined ellipsoid remains inside it.

Figure 2: Data-driven invariant ellipsoid

Fig. 3 shows the upper bound of the output peak and set of randomly selected
trajectories starting at the boundary of the invariant ellipsoid. The figure
shows that the output $\left\lVert z(t)\right\rVert$ of all trajectories
inside the invariant ellipsoid remains below the obtained upper bound.

Figure 3: Data-driven upper bound for output peak

The data-driven results of our proposed approach are almost identical to the
solution of the SDP (30) assuming the system dynamics is known. The upper
bound for output peak derived from the SDP of (30) is $3.2915$ and the
invariant ellipsoid is $x^{T}Px\leq 1$ such that:

$P_{LMI}=\begin{bmatrix}0.092426&0.001406\\\ 0.001406&0.015873\end{bmatrix}$
(35)

##  VI Data-driven Evaluation for $\mathbf{L}_{2}$ gain and RMS gain

For LTI system(1), $\text{L}_{2}$ gain is defined as:

$\sup_{\left\lVert u\right\rVert_{2}\neq 0}\frac{\left\lVert
z\right\rVert_{2}}{\left\lVert u\right\rVert_{2}}$ (36)

where the supremum is over all nonzero trajectories of the LTI system starting
from $x(0)=0$. In [16], If there exists a quadratic function
$\mathcal{E}=\zeta^{T}P\zeta$, $P>0$ and for all $t\geq 0$,

$\frac{d}{\text{dt}}V(x)+z^{T}z-\gamma^{2}u^{T}u\leq 0\quad\forall x\text{ and
}u\text{ satisfying (\ref{LTI system}) }$ (37)

such that $\gamma\geq 0$, then the $\mathbf{L}_{2}$ gain of the LTI is less
than $\gamma$. To show that, integrate both sides of (37) from $0$ to $T$ with
$x(0)=0$, we obtain

$V(x(T))+\int_{0}^{T}\left(z^{T}z-\gamma^{2}u^{T}u\right)dt\leq 0$ (38)

$V(x(T))\geq 0$, so we conclude that
$\int_{0}^{T}z^{T}zdt\leq\gamma^{2}\int_{0}^{T}u^{T}udt$.

  * •

Take the limit $T\to\infty$, then the system’s $\mathbf{L}_{2}$ gain is less
than $\gamma$

  * •

Divide by $T$ and take the limit $T\to\infty$, then
$\frac{\mathbf{RMS}(z)}{\mathbf{RMS}(u)}\leq\gamma$, Therefore, the system’s
$\mathbf{RMS}$ gain is less than $\gamma$

Given the system dynamics ($A,B,C$), equation (37) can be written as follows:

$x^{T}(PA+A^{T}P+C^{T}C)x+2x^{T}PBu-\gamma^{2}u^{T}u\leq 0$ (39)

that can be formulated as the following LMI:

$\begin{bmatrix}PA+A^{T}P+C^{T}C&PB\\\ B^{T}P&-\gamma^{2}I\end{bmatrix}\leq 0$
(40)

Hence, the smallest upper bound for the LTI system’s $\mathbf{L}_{2}$ or
$\mathbf{RMS}$ gain can be obtained by minimizing $\gamma$ over the variables
$P$ and $\gamma$ while satisfying (40) and $P>0$. This method is based on the
full knowledge of the system’s dynamics ($A,B,C$). Our approach can be also
used here to solve this problem given only a single trajectory of the LTI
system starting from $x(0)=0$. Instead of solving the SDP problem of (40), the
problem can be reformulated as an LP problem by minimizing $\beta=\gamma^{2}$
while the conditions $P>0$ and (37) at any data point $i$ along the given
trajectory are satisfied in the following format:

$\begin{gathered}x(i)^{T}Px(i)>0\\\ 2x(i)^{T}P\dot{x}(i)+z(i)^{T}z(i)-\beta
u(i)^{T}u(i)\leq 0\end{gathered}$ (41)

where $\gamma^{2}$ is replaced by $\beta$ in both the objective and the second
condition to preserve the linearity. Therefore, the Linear program can be
written as:

$\displaystyle\underset{p,\beta}{\text{minimize}}$ $\displaystyle\beta$ (42)
subject to $\displaystyle L_{1}p>0$ $\displaystyle L_{2}p+l_{z}-\beta
l_{u}\leq 0$

where $p=\text{Vec}(P)$ and $L_{1}$,$L_{2}$ and $l_{z}$ are defined by (6),(7)
and (25) respectively. $l_{u}$ is:

$l_{u}=\begin{bmatrix}u(0)^{T}u(0)\\\ u(1)^{T}u(1)\\\ \vdots\\\
u(N)^{T}u(N)\end{bmatrix}$ (43)

Numerical example: Given the following LTI system:

$\displaystyle\dot{x}$ $\displaystyle=\begin{bmatrix}0&1\\\
-1&-2\end{bmatrix}x+\begin{bmatrix}1\\\ 2\end{bmatrix}u$ (44) $\displaystyle
z$ $\displaystyle=\begin{bmatrix}4&1\end{bmatrix}x$

Assuming that ($A,B,C$) are known. So,we can get the the system
$\mathbf{L}_{2}$ gain by solving the following SDP:

$\displaystyle\underset{P,\gamma}{\text{minimize}}$ $\displaystyle\gamma$ (45)
subject to $\displaystyle P>0,\quad(\ref{rbl LMI})$

So, $\mathbf{L}_{2}$ gain of the system is $\mathbf{15}$. Using our approach,
$\gamma$ can be obtained by exciting the system with a unit step input from
$0$ to a final time $T$, then measure the input-output ($u$,$z$) and states
($x,\dot{x}$) data of a single trajectory. Therefore, the linear program (42)
can be solved to get $\gamma=\sqrt{\beta}$. The time step used is $0.01$ sec.
Table I shows that as the length of the trajectory increases, the result’s
accuracy will be improved.

TABLE I: Learned $\mathbf{L}_{2}$ gain Final time $T$ [Seconds] | Learned $\mathbf{L}_{2}$ gain ($\gamma$) [SDP solution=15]
---|---
2 | 7.36147
4 | 12.17092
6 | 14.30461
8 | 14.86621
10 | 14.97684
12 | 14.99619
14 | 14.9994
16 | 14.99990

Note that $\mathbf{L}_{2}$ gain of the LTI system is the $\mathbf{H}_{\infty}$
norm of its transfer function $C(sI-A)^{-1}B$. Therefore, using our approach,
the $\mathbf{H}_{\infty}$ norm of the LTI system can be learned accurately
using only a single trajectory and without need to know its dynamic matrices.

##  VII Conclusion and future work

In this work, a number of data-driven techniques that aim to evaluate various
metrics of LTI systems such as the upper bound on the output ’s energy and
peak as well as the $\mathbf{L}_{2}$ and $\mathbf{RMS}$ gains were introduced.
In addition, the exact and approximate construction of the Lyapunov function
of LTI systems were proposed. To demonstrate the proposed methodologies, a
number of numerical examples were given and thoroughly discussed. As for
future work, we are extending the proposed data-driven construction for the
Lyapunov function to a non-quadratic Lyapunov function for hybrid linear time-
invariant systems and specific classes of nonlinear systems. Finally, we are
extending the proposed data-driven methodologies to linear systems where
uncertainties happen and the measurements are noisy.

## References

  * [1] M. P. Deisenroth, D. Fox, and C. E. Rasmussen, “Gaussian processes for data-efficient learning in robotics and control,” IEEE transactions on pattern analysis and machine intelligence, vol. 37, pp. 408–423, 2013.
  * [2] K. Chua, R. Calandra, R. McAllister, and S. Levine, “Deep reinforcement learning in a handful of trials using probabilistic dynamics models,” arXiv preprint arXiv:1805.12114, 2018.
  * [3] A. Aswani, H. Gonzalez, S. S. Sastry, and C. Tomlin, “Provably safe and robust learning-based model predictive control,” Automatica, vol. 49, pp. 1216–1226, 2013.
  * [4] F. Berkenkamp, A. P. Schoellig, and A. Krause, “Safe controller optimization for quadrotors with gaussian processes,” in 2016 IEEE International Conference on Robotics and Automation (ICRA), pp. 491–496, 2016.
  * [5] T. Beckers, D. Kulić, and S. Hirche, “Stable gaussian process based tracking control of euler-lagrange systems,” Automatica, vol. 103, pp. 390–397, 2019.
  * [6] A. Lederer, A. Capone, and S. Hirche, “Parameter optimization for learning-based control of control-affine systems,” in Learning for Dynamics and Control, pp. 465–475, 2020.
  * [7] J. F. Fisac, A. K. Akametalu, M. N. Zeilinger, S. Kaynama, J. Gillula, and C. J. Tomlin, “A general safety framework for learning-based control in uncertain robotic systems,” IEEE Transactions on Automatic Control, vol. 64, pp. 2737–2752, 2018.
  * [8] A. Lederer, A. Capone, T. Beckers, J. Umlauft, and S. Hirche, “The impact of data on the stability of learning-based control,” in Learning for Dynamics and Control, pp. 623–635, 2021.
  * [9] F. Pukelsheim, Optimal design of experiments.  SIAM, 2006.
  * [10] P. Hennig and C. J. Schuler, “Entropy search for information-efficient global optimization.,” Journal of Machine Learning Research, vol. 13, 2012.
  * [11] A. Lederer, A. Capone, J. Umlauft, and S. Hirche, “How training data impacts performance in learning-based control,” IEEE Control Systems Letters, vol. 5, pp. 905–910, 2020.
  * [12] K. S. Narendra and J. Balakrishnan, “A common Lyapunov function for stable lti systems with commuting a-matrices,” IEEE Transactions on automatic control, vol. 39, pp. 2469–2471, 1994.
  * [13] O. Mason and R. Shorten, “On linear copositive Lyapunov functions and the stability of switched positive linear systems,” IEEE Transactions on Automatic Control, vol. 52, pp. 1346–1349, 2007.
  * [14] H. Ravanbakhsh and S. Sankaranarayanan, “Learning Lyapunov (potential) functions from counterexamples and demonstrations,” arXiv preprint arXiv:1705.09619, 2017.
  * [15] A. Abate, D. Ahmed, M. Giacobbe, and A. Peruffo, “Formal synthesis of Lyapunov neural networks,” IEEE Control Systems Letters, vol. 5, pp. 773–778, 2020.
  * [16] S. Boyd, L. El Ghaoui, E. Feron, and V. Balakrishnan, Linear matrix inequalities in system and control theory.  SIAM, 1994.
  * [17] A. Papachristodoulou and S. Prajna, “On the construction of Lyapunov functions using the sum of squares decomposition,” in Proceedings of the 41st IEEE Conference on Decision and Control, 2002., vol. 3, pp. 3482–3487, IEEE, 2002.
  * [18] A. Papachristodoulou and S. Prajna, “A tutorial on sum of squares techniques for systems analysis,” in Proceedings of the 2005, American Control Conference, 2005., pp. 2686–2700, IEEE, 2005.
  * [19] Y.-C. Chang, N. Roohi, and S. Gao, “Neural Lyapunov control,” arXiv preprint arXiv:2005.00611, 2020.
  * [20] S. Chen, M. Fazlyab, M. Morari, G. J. Pappas, and V. M. Preciado, “Learning Lyapunov functions for hybrid systems,” in Proceedings of the 24th International Conference on Hybrid Systems: Computation and Control, pp. 1–11, 2021.
  * [21] M. Grant and S. Boyd, “CVX: Matlab software for disciplined convex programming, version 2.1,” 2014.