# Neural Lyapunov Control for Nonlinear Systems
with Unstructured Uncertainties

Shiqing Wei, Prashanth Krishnamurthy, and Farshad Khorrami The authors are
with Control/Robotics Research Laboratory, Dept. of ECE, NYU Tandon School of
Engineering, 5 Metrotech Center, Brooklyn, NY 11201, USA. {shiqing.wei,
prashanth.krishnamurthy<EMAIL_ADDRESS>work was supported in part by
ARO grant W911NF-22-1-0028 and in part by the New York University Abu Dhabi
(NYUAD) Center for Artificial Intelligence and Robotics, funded by Tamkeen
under the NYUAD Research Institute Award CG010.

###### Abstract

Stabilizing controller design and region of attraction (RoA) estimation are
essential in nonlinear control. Moreover, it is challenging to implement a
control Lyapunov function (CLF) in practice when only partial knowledge of the
system is available. We propose a learning framework that can synthesize
state-feedback controllers and a CLF for control-affine nonlinear systems with
unstructured uncertainties. Based on a regularity condition on these
uncertainties, we model them as bounded disturbances and prove that a CLF for
the nominal system (estimate of the true system) is an input-to-state stable
control Lyapunov function (ISS-CLF) for the true system when the CLF’s
gradient is bounded. We integrate the robust Lyapunov analysis with the
learning of both the control law and CLF. We demonstrate the effectiveness of
our learning framework on several examples, such as an inverted pendulum
system, a strict-feedback system, and a cart-pole system.

##  I Introduction

While the knowledge of a control Lyapunov function (CLF) for a system can
enable the implementation of universal controllers [1], designing a CLF for a
given real-world nonlinear system can be highly challenging, especially when
the system’s dynamics are uncertain. Furthermore, when faced with a real-world
uncertain system, proving that a controller provides global stabilization
might be infeasible or even ill-defined, and one instead needs to estimate the
region of attraction (RoA) of the closed-loop system using online data.
Motivated by the several advances in learning-based methods for various
control design tasks (_e.g_., [2, 3, 4, 5, 6]), we consider the problems of
designing a controller for uncertain systems and estimating the achieved RoA
from a learning-based or data-driven approach in this paper.

Specifically, we address the problem of synthesizing state-feedback
controllers and estimating the RoA for control-affine nonlinear systems with
unstructured uncertainties. We model these uncertainties as a bounded
disturbance but with no assumptions about the source of these uncertainties.
They may result from incorrect model parameters or dynamics that are not
reflected in the model. Our work takes advantage of the rich literature on
Lyapunov theory [7, 8, 9, 10], and a data-driven approach is adopted.

Several approaches have been explored in the existing literature to enhance
the robustness of Lyapunov analysis and control designs to uncertainties in
the underlying dynamics of real-world systems, such as an unknown constant
parameter [11], an unknown Gaussian process [12], or an unknown linearly
parameterized control-affine system [13]. In particular, the authors of [14]
consider the modeling error as a disturbance and propose the projection-to-
state stability approach to characterize the tracking error. However, one
common point among the above methods is that they either derive controllers
[11, 14, 13] or refine the RoA [12] based on a given CLF. In practice, such a
CLF may not be readily available for complex nonlinear systems and may provide
a conservative estimation of the actual RoA.

Another area closely related to this paper is the automated formulation of
Lyapunov functions. Lyapunov functions for polynomial systems can be found by
solving linear matrix inequalities (LMIs) [15]. Approximation of Lyapunov
functions by sum-of-squares (SOS) polynomials can be found through the
solution of a semidefinite programming (SDP) problem [16]. Computational
methods for Lyapunov functions have been reviewed in [17]. More recently,
neural networks have been used to approximate a Lyapunov function with SMT
(Satisfiability Modulo Theories) solvers being a verification tool [18, 19].
However, these learning-based approaches (_e.g_., [18, 5]) have been developed
assuming exact knowledge of the system. In [2], knowledge of closed-loop
dynamics is needed when computing the Lipschitz constant. The method proposed
in [5], Neural Lyapunov Redesign (NLR), is an offline method that finds a
neural Lyapunov function and a stabilizing controller assuming known system
dynamics.

Our contributions: We consider control-affine nonlinear systems with
unstructured uncertainties and propose a learning framework that
simultaneously learns the following: a state-feedback controller that seeks to
enlarge the RoA, a CLF that can be used to estimate the RoA, and an improved
model of the uncertain system dynamics (starting with an initial nominal
approximate system model). Inspired by [14], we model the unstructured
uncertainties as a bounded disturbance. We prove that when the gradient of the
CLF is bounded, and the disturbance is bounded by a particular quantity, a CLF
for the nominal system is an input-to-state stable control Lyapunov function
(ISS-CLF) for the true system and thus can correctly estimate the RoA for the
true system. We apply a machine learning-based data-driven approach for
learning the controller, CLF, and system dynamics and demonstrate the
effectiveness of our learning framework on three different examples.

##  II Preliminaries

Let $\mathcal{X}\subset\mathbb{R}^{n}$ and $\mathcal{U}\subset\mathbb{R}^{m}$
be the state and control input spaces (in general, subsets of $\mathbb{R}^{n}$
and $\mathbb{R}^{m}$ to model physical constraints of real-world systems).
Consider the following control-affine system:

$\dot{x}=f(x)+g(x)u$ (1)

with drift dynamics $f:\mathcal{X}\to\mathbb{R}^{n}$ and actuation matrix
$g:\mathcal{X}\to\mathbb{R}^{n\times m}$. To ensure the existence and
uniqueness of the solution, we assume that $f$ and $g$ are Lipschitz
continuous on $\mathcal{X}$. We further assume $0\in\mathcal{X}$ and
$f(0)+g(0)u_{0}=0$ for a certain $u_{0}\in\mathcal{U}$ and also assume
controllability of the system. Introducing a disturbance in (1), we consider
the perturbed system:

$\dot{x}=f(x)+g(x)u+d$ (2)

where $d\in\mathcal{D}$ is the disturbance assumed to be essentially bounded
in time (i.e., bounded everywhere except possibly on a set of measure zero)
and $\mathcal{D}\subset\mathbb{R}^{n}$ is the disturbance space. A natural
framework for modeling the effects of perturbations is given by the widely
used notions of input-to-state stability (ISS) [20] and input-to-state stable
control Lyapunov functions (ISS-CLF) [9].

######  Definition 1 (CLF and ISS-CLF).

A class $C^{1}$ function $V:\mathcal{X}\to\mathbb{R}_{+}$ is a CLF for
$(\ref{eq:true_sys})$ on $\mathcal{X}$ if there exist
$\alpha_{1},\alpha_{2},\alpha_{3}\in\mathcal{K}_{\infty}$ such that for all
$x\in\mathcal{X}$:

$\displaystyle\alpha_{1}(\lVert x\rVert)\leq V(x)$
$\displaystyle\leq\alpha_{2}(\lVert x\rVert),$ (3)
$\displaystyle\inf_{u\in\mathcal{U}}\dot{V}(x,u)$
$\displaystyle\leq-\alpha_{3}(\lVert x\rVert).$ (4)

$V$ is an ISS-CLF for (2) on $\mathcal{X}$ if it satisfies (3) and
additionally, there exist $\alpha_{4},\rho\in\mathcal{K}_{\infty}$ such that

$\lVert x\rVert\geq\rho\left(\operatorname*{ess\,sup}_{\tau\geq t_{0}}\lVert
d(\tau)\rVert\right)\Rightarrow\inf_{u\in\mathcal{U}}\dot{V}(x,u,d)\leq-\alpha_{4}(\lVert
x\rVert)$ (5)

for all $x\in\mathcal{X}$ and $d\in\mathcal{D}$.

The existence of a CLF implies the existence of a state-feedback controller
$k:\mathcal{X}\rightarrow\mathcal{U}$ [8] such that
$\dot{V}\leq-\alpha_{3}(||x||)$ in the closed-loop system (i.e., when control
inputs are generated using $k$). We refer to controllers satisfying this
inequality as admissible. We note that for the definition of a CLF,
$\alpha_{1}$ and $\alpha_{2}$ need only to be class $\mathcal{K}$ functions,
and condition (4) can be reduced to $\inf_{u\in\mathcal{U}}\dot{V}(x,u)<0$ to
guarantee local asymptotic stability of the origin [8].

##  III Uncertain Dynamics

Let $\hat{f}:\mathcal{X}\to\mathbb{R}^{n}$ and
$\hat{g}:\mathcal{X}\to\mathbb{R}^{n\times m}$ be Lipschitz continuous
functions denoting the estimates of $f$ and $g$ in (1). Then, system (1),
called true system, can be written as

$\dot{x}=\hat{f}(x)+\hat{g}(x)u+\underbrace{(f(x)-\hat{f}(x))+(g(x)-\hat{g}(x))u}_{d}$
(6)

where the estimation errors are written as the disturbance signal $d$. If $d$
is essentially bounded in time, system (6) can be seen as the perturbed system
of the nominal system

$\dot{x}=\hat{f}(x)+\hat{g}(x)u.$ (7)

Finding a CLF for (1) valid on the entire state space $\mathcal{X}$ is
difficult for nonlinear systems, and the conditions (3) and (4) are usually
satisfied only on a compact subset $\mathcal{C}$ of $\mathcal{X}$. Then, given
a continuous state-feedback controller $k:\mathcal{X}\rightarrow\mathcal{U}$,
the estimation error $d$ is indeed bounded on $\mathcal{C}$ as a result of the
continuity of $f,\hat{f},g,\hat{g}$, and $k$ on $\mathcal{X}$.

Given a CLF for the nominal system (7) on $\mathcal{C}$, by (6), the time
derivative of $V$ can be written as

$\dot{V}(x,k(x),d)=\underbrace{\nabla{V}(x)^{\top}(\hat{f}(x)+\hat{g}(x)k(x))}_{\hat{\dot{V}}(x,k(x))}+\underbrace{\nabla
V(x)^{\top}d}_{\delta}$ (8)

where $\nabla V$ is the gradient of $V$ and $\delta=\nabla V(x)^{\top}d$ can
be seen as a disturbance in the dynamics of $V$. Since $V$ is $C^{1}$ on
$\mathcal{X}$, $\lVert\nabla V\rVert$ is bounded on $\mathcal{C}$ (with bound
$L_{V}>0$), and the disturbance term $\delta$ is also bounded.

######  Property 1.

Let $V$ be a CLF with
$\alpha_{1},\alpha_{2},\alpha_{3}\in\mathcal{K}_{\infty}$ for system (7)
associated with an admissible continuous controller
$k:\mathcal{X}\rightarrow\mathcal{U}$ over a compact set
$\mathcal{C}\subset\mathcal{X}$. Let $\Omega\subset\mathcal{C}$ be a sublevel
set of $V$, and $\alpha_{p}$ and $\alpha_{q}$ be class $\mathcal{K}_{\infty}$
functions such that $\alpha_{p}+\alpha_{q}=\alpha_{3}$. If $\lVert
x\rVert\geq\alpha_{q}^{-1}\left(L_{V}\max_{x\in\Omega}\lVert d\rVert\right)$
for all $x\in\partial\Omega$ (the boundary of set $\Omega$), then $V$ is an
ISS-CLF for system (6) with controller $k$ over $\Omega$.

###### Proof.

Let $c=V(x)$ for any $x\in\partial\Omega$. Then, for all $x\in\partial\Omega$,

$\dot{V}(x,k(x),d)\leq-\alpha_{p}(\lVert x\rVert)-\alpha_{q}(\lVert
x\rVert)+\delta\leq-\alpha_{p}(\lVert x\rVert)<0,$

and by Nagumo’s Theorem, $V(x(t))\in[0,c]$ for $t\geq t_{0}$ if
$V(x(t_{0}))\in[0,c]$. Therefore, the set $\Omega$ is forward invariant, which
means that for any trajectory starting within $\Omega$, the disturbance along
the trajectory will be bounded by $\max_{x\in\Omega}\lVert d\rVert$. Note that
the inverse of a class $\mathcal{K}_{\infty}$ function is defined on
$\mathbb{R}_{+}$ and belongs to class $\mathcal{K}_{\infty}$,
$\rho:r\mapsto\alpha_{q}^{-1}(L_{V}r)$ is therefore a class
$\mathcal{K}_{\infty}$ function. $V$ is an ISS-CLF on $\Omega$ since if
$\lVert x\rVert\geq\rho\left(\sup_{\tau\geq t_{0}}\lVert
d(\tau)\rVert\right)$, by (4) and (8),

$\dot{V}(x,k(x),d)\leq-\alpha_{p}(\lVert x\rVert)-\alpha_{q}(\lVert
x\rVert)+\delta\leq-\alpha_{p}(\lVert x\rVert)$

where $\alpha_{p}\in\mathcal{K}_{\infty}$, and condition (5) is satisfied. ∎

The next property is a direct result of Nagumo’s Theorem.

######  Property 2.

Let $V$ be an ISS-CLF for (6) under a state-feedback controller
$k:\mathcal{X}\rightarrow\mathcal{U}$ on a compact set $\mathcal{C}$ and
$\Omega\subseteq\mathcal{C}$ be a sublevel set of $V$. $\Omega$ is forward
invariant if $\lVert x\rVert\geq\rho\left(\lVert d\rVert\right)$ for all
$x\in\partial\Omega$.

Finally, based on [9, Theorem 1], we have:

######  Property 3.

If $V$ is an ISS-CLF for system (6) under an admissible state-feedback
controller $k:\mathcal{X}\rightarrow\mathcal{U}$ on $\mathcal{C}$, then (6) is
ISS with the controller $k$ on $\mathcal{C}$.

##  IV Learning the Controller and ISS-CLF

###  IV-A Neural Network Structures

As part of our learning framework, we update the nominal dynamics during
training. For this purpose, we decompose the estimates $\hat{f}$ and $\hat{g}$
into two parts:

$\hat{f}=f_{0}+f_{\theta_{1}},\quad\hat{g}=g_{0}+g_{\theta_{2}}$ (9)

where $f_{0}$ and $g_{0}$ depend on the initial knowledge of the system, and
the residual parts $f_{\theta_{1}}$ and $g_{\theta_{2}}$ are neural networks
(parameterized by ${\theta_{1}}$ and ${\theta_{2}}$, respectively). If the
system is completely unknown, we may start with $f_{0}=g_{0}=0$.

Let $V_{\theta_{3}}:\mathcal{X}\rightarrow\mathbb{R}_{+}$, called Lyapunov
candidate, be a neural network parameterized by $\theta_{3}$. To ensure that
$V_{\theta_{3}}$ satisfies (3) in Definition 1, its structure is defined as:

$V_{\theta_{3}}(x)=x^{\top}(MM^{\top}+\gamma I)x+\phi(x)^{\top}\phi(x)$ (10)

where $M\in\mathbb{R}^{n\times n}$ is a trainable lower triangular matrix,
$I\in\mathbb{R}^{n\times n}$ is the identity matrix, $\gamma>0$ is a constant
and $\phi:\mathcal{X}\rightarrow\mathbb{R}^{d}$ is a neural network.
$\theta_{3}$ includes both $M$ and weights in $\phi$. Following [2], $\phi$ is
a composition of linear layers (with no bias term) and activation functions
(Lipschitz continuous functions with a trivial null space), and its layer
dimensions are increasing. Let $d_{\ell}$ and $W_{\ell}$ be the output
dimension and layer weights of layer $\ell$, and

$W_{\ell}=\begin{bmatrix}G_{\ell 1}^{\top}G_{\ell 1}+\varepsilon
I_{d_{\ell-1}}\\\ G_{\ell 2}\end{bmatrix}$ (11)

where $G_{\ell 1}\in\mathbb{R}^{q_{\ell}\times d_{\ell-1}}$ for some integer
$q_{\ell}\geq 1$, $G_{\ell 2}\in\mathbb{R}^{(d_{\ell}-d_{\ell-1})\times
d_{\ell-1}}$, $I_{d_{\ell-1}}\in\mathbb{R}^{d_{\ell-1}\times d_{\ell-1}}$ is
the identity matrix, and $\varepsilon>0$ is a constant. Note that (11) is
enforced during training to ensure that $\phi$ has a trivial null space,
implying that $\phi(x)^{\top}\phi(x)$ is positive definite. Our Lyapunov
candidate differs from that of [2] in two aspects. Firstly, the
$\phi(x)^{\top}\phi(x)$ term is not necessarily lower bounded by a class
$\mathcal{K}_{\infty}$ function; so we have added the $\gamma x^{\top}x$ term
to achieve this property. Secondly, we added the $x^{\top}MM^{\top}x$ term to
better capture quadratic behaviors and expect the $\phi(x)^{\top}\phi(x)$ term
to capture other nonlinear behaviors. The upper bounding function
$\alpha_{2}\in\mathcal{K}_{\infty}$ can be constructed based on
$V_{\theta_{3}}$.

The state-feedback controller
$u_{\theta_{4}}:\mathcal{X}\rightarrow\mathcal{U}$ is also a neural network
(parameterized by $\theta_{4}$). The activation functions of $u_{\theta_{4}}$
are chosen such that it is continuous.

###  IV-B Learning Algorithm

Let $\mathcal{X}_{\tau}\subset\mathcal{X}$ be a discretization of
$\mathcal{X}$ with $\lVert x-[x]_{\tau}\rVert_{2}\leq\tau/2$, where
$[x]_{\tau}$ denotes the closest point in $\mathcal{X}_{\tau}$ to
$x\in\mathcal{X}$. Based on the structure and our choice of activation
functions of the neural networks, the estimated time derivative of
$V_{\theta_{3}}$

$\hat{\dot{V}}_{\theta_{3}}(x,u_{\theta_{4}}(x))=\nabla
V_{\theta_{3}}(x)^{\top}(\hat{f}(x)+\hat{g}(x)u_{\theta_{4}}(x))$ (12)

is locally Lipschitz continuous on $\mathcal{X}$. For $V_{\theta_{3}}$ to be a
CLF for the nominal system (7), condition (4) in Definition 1 has to be
satisfied on a certain set $\mathcal{C}$ (to be determined), i.e.,

$\hat{\dot{V}}_{\theta_{3}}(x,u_{\theta_{4}}(x))\leq-\kappa\lVert
x\rVert^{2}\text{ on }\mathcal{C}$ (13)

where we have chosen $\alpha_{3}(\lVert x\rVert)=\kappa\lVert x\rVert^{2}$ and
$\kappa$ is a positive constant. The Lyapunov loss is defined as

$\displaystyle\mathcal{L}_{\theta_{3},\theta_{4}}$
$\displaystyle=\frac{\lambda_{\text{RoA}}}{N_{i}}\sum_{x\in\mathcal{S}_{i}}\text{ReLU}[\hat{\dot{V}}_{\theta_{3}}(x,u_{\theta_{4}}(x))+\kappa\lVert
x\rVert^{2}+\epsilon]$
$\displaystyle+\frac{\lambda_{\text{Lip}}}{N_{i}}\sum_{x\in\mathcal{S}_{i}}\lVert\nabla
V_{\theta_{3}}(x)\rVert$ (14)

where $\lambda_{\text{RoA}},\lambda_{\text{Lip}},$ and $\epsilon$ are positive
constants, $\mathcal{S}_{i}$ is the training set at iteration $i$, $N_{i}$ is
the number of training samples in $\mathcal{S}_{i}$, and
$\text{ReLU}(x)=\max(0,x)$ stands for the Rectified Linear Unit. The first
term of (IV-B) accounts for condition (13), and $\epsilon$ is a positive
constant offset. The second term aims to limit the norm of $\nabla
V_{\theta_{3}}$. Compared with [2], our loss function has a more general
structure with a regularization term that is designed to limit
$L_{V_{\theta_{3}}}$ (upper bound of $\lVert\nabla V_{\theta_{3}}\rVert$), as
required by part of our methodology.

Input: Dynamics network $f_{\theta_{1}}$ and $g_{\theta_{2}}$, Lyapunov
candidate $V_{\theta_{3}}$, state-feedback controller $u_{\theta_{4}}$, loss
function $\mathcal{L}$, mesh $\mathcal{X}_{\tau}$, initial stable set
$\mathcal{X}_{0}$222We assume that the algorithm starts with an initial
locally stable controller., learning rate scheduler StepLR, level multiplier
$\eta_{0}$ and step $k_{\eta}$, and number of iterations $N$.

1

2Pretrain $\hat{\dot{V}}_{\theta_{3}}$ on $\mathcal{X}_{0}$ by minimizing
$\operatorname{MSE}(\hat{\dot{V}}_{\theta_{3}},\tilde{\dot{V}}_{\theta_{3}})$
via stochastic gradient descent (SGD);

3

4for  _$i=1,2,...,N$_ do

5 Determine $\mathcal{X}_{\text{stable}}$ from sampled trajectories;

6

7 $c_{i}=\max_{c>0}c,$ s.t.
$\hat{\dot{V}}_{\theta_{3}}(x,u_{\theta_{4}}(x))\leq-\kappa\lVert x\rVert^{2}$
for all $x\in\mathcal{X}_{\text{stable}}\cap\mathcal{V}(c)$ and
$\mathcal{V}(c)\subset\mathcal{X}$;

8

9 $\eta_{i}=1+\frac{\eta_{0}}{1+\left\lfloor i/k_{\eta}\right\rfloor}$;

10

11 $\mathcal{S}_{i}=\mathcal{V}(\eta_{i}c_{i})$;

12

13 Update $f_{\theta_{1}}$ and $g_{\theta_{2}}$ by minimizing
$\operatorname{MSE}(\hat{\dot{V}}_{\theta_{3}},\tilde{\dot{V}}_{\theta_{3}})$
on $\mathcal{S}_{i}$ via SGD;

14

15 Determine learning rate for $\theta_{3}$ and $\theta_{4}$ by StepLR;

16

17 Update $V_{\theta_{3}}$ and $u_{\theta_{4}}$ by minimizing the loss
function (IV-B) on $\mathcal{S}_{i}$ via SGD;

18 end for

Algorithm 1 Learning to maximize the RoA

Our learning framework learns the nominal dynamics, the Lyapunov candidate,
and the controller at the same time. Instead of directly working on full-state
dynamics, we update the nominal dynamics using the dynamics of the CLF during
learning. During each iteration, our approach first updates $f_{\theta_{1}}$
and $g_{\theta_{2}}$ by minimizing the mean squared error (MSE) between
$\hat{\dot{V}}_{\theta_{3}}$ (by (12)) and $\tilde{\dot{V}}_{\theta_{3}}$
(approximated value of $\dot{V}_{\theta_{3}}$ obtained by numerical
differentiation), and then updates $V_{\theta_{3}}$ and $u_{\theta_{4}}$ by
minimizing the Lyapunov loss (IV-B).

Let $\mathcal{V}(c)=\\{\,x\in\mathbb{R}^{n}\;|\;V_{\theta_{3}}(x)\leq c\,\\}$
denote the $c$-sublevel set of $V_{\theta_{3}}$. At iteration $i$, we collect
from the system the trajectories initialized at $\mathcal{X}_{\tau}$,
determine the set of stable initial states***The stable initial states are the
initial states from which the trajectory asymptotically converges to the
origin., and note them as
$\mathcal{X}_{\text{stable}}\subset\mathcal{X}_{\tau}$. The estimation of the
RoA $R_{i}$ is given by largest sublevel set of $V_{\theta_{3}}$ contained in
$\mathcal{X}$ where (13) is satisfied, i.e., $R_{i}=\mathcal{V}(c_{i})$, where
$c_{i}=\max_{c>0}c$ subject to

$\hat{\dot{V}}_{\theta_{3}}(x,u_{\theta_{4}}(x))\leq-\kappa\lVert
x\rVert^{2}\text{ for all }x\in\mathcal{X}_{\text{stable}}\cap\mathcal{V}(c)$
(15)

and $\mathcal{V}(c)\subset\mathcal{X}$. To enlarge the estimated RoA, we
identify an exploration region by defining a level multiplier $\eta_{i}>1$ to
include more states in the training set. The training set at iteration $i+1$
is $\mathcal{S}_{i+1}=\mathcal{V}(\eta_{i}c_{i})$. The details are provided in
Algorithm 1.

###  IV-C Verification Condition

The Lyapunov candidate $V_{\theta_{3}}$ minimizing (IV-B) is called Lyapunov-
like. Define the following first-order logic formula

$\displaystyle\Phi_{\zeta}(x)$ $\displaystyle=\left(\lVert
x\rVert\geq\zeta,V_{\theta_{3}}(x)\leq c\right)$
$\displaystyle\wedge\left(\hat{\dot{V}}_{\theta_{3}}(x,u_{\theta_{4}}(x))+\kappa\lVert
x\rVert^{2}\geq 0\right)$ (16)

to check the violation of the Lyapunov condition (13) where $\zeta>0$ is a
small constant that rules out a small region around the origin to avoid
numerical instabilities, and $c$ is found from Algorithm 1. We use dReal [21],
an SMT solver for nonlinear constraints, to solve (IV-C). dReal runs a delta-
complete algorithm whose numerical error bound is specified by the user. If
dReal cannot find any solution satisfying (IV-C), then condition (13) is
certified on $R=\mathcal{V}(c)$. In this case, $V_{\theta_{3}}$ is a CLF on
$R$ for the nominal system (7) under controller $u_{\theta_{4}}$ by Definition
1 and is an ISS-CLF for (6) by Property 1. Further, $R$, as a sublevel set of
$V_{\theta_{3}}$, is forward invariant by Property 2, and system (6) is ISS
with $u_{\theta_{4}}$ by Property 3. However, for computational tractability
(specifically, due to the limit on the number of parameters in the dReal SMT
solver that we use), we can only run the SMT solver on a simplified setting,
an example for which is provided in Section V-A.

Note that even if condition (13) is violated on $R$, it is still possible for
states to remain bounded under $u_{\theta_{4}}$. Let
$\delta^{\prime}=\max_{x\in
R}\left[\hat{\dot{V}}_{\theta_{3}}(x,u_{\theta_{4}}(x))+\kappa\lVert
x\rVert^{2}\right]$ denote upper bound of the violation on $R$. With $d$
defined in (6) and $\delta=\nabla V_{\theta_{3}}(x)^{\top}d$, we have
$\dot{V}_{\theta_{3}}(x,u_{\theta_{4}}(x))\leq-\kappa\lVert
x\rVert^{2}+\delta+\delta^{\prime}$ on $R$. Let $\alpha_{p}$ and $\alpha_{q}$
be class $\mathcal{K}_{\infty}$ functions such that
$\alpha_{p}+\alpha_{q}=\kappa\lVert x\rVert^{2}$. If $\lVert
x\rVert\geq\alpha_{q}^{-1}\left(\lVert\delta+\delta^{\prime}\rVert\right)$ for
all $x\in\partial R$, we have
$\dot{V}_{\theta_{3}}(x,u_{\theta_{4}}(x))\leq-\alpha_{p}(\lVert x\rVert)$ for
all $x\in\partial R$. Let $y(x)=V_{\theta_{3}}(x,u_{\theta_{4}}(x))$ and $c$
be the corresponding level value of $R$. Then, by Nagumo’s Theorem,
$y(x(t))\in[0,c]$ for $t\geq t_{0}$ if $y(x(t_{0}))\in[0,c]$. Hence, states
will stay bounded.

##  V Experiments

In this section, we study the efficacy of our proposed learning framework. The
Lyapunov candidate $V_{\theta_{3}}$ is randomly initialized and then
pretrained to

$V_{0}(x)=0.1x^{\top}x$ (17)

on $\mathcal{X}$. The controller $u_{\theta_{4}}$ has two parts:

$u_{\theta_{4}}(x)=\textrm{LS}_{a,b,m_{a},m_{b}}(u_{0}(x)+\psi(x))$ (18)

where $u_{0}$ is an initial locally stabilizing controller (fixed during
training), $\psi$ is the trainable part, and $\textrm{LS}_{a,b,m_{a},m_{b}}$
is the loose saturation filter (see Fig. 1). The thresholds $a$ and $b$ are
fixed, but the slopes $m_{a}$ and $m_{b}$ (initialized to zero) are trainable
and included in the controller parameters $\theta_{4}$. The motivation for
using a loose saturation filter is to more clearly demonstrate that our
algorithm is indeed learning to enlarge the RoA, as for some systems, a
locally stable controller could automatically be a globally stable controller
in the absence of such a filter. Also, a similar setup is used in reference
[5], and we use the loose saturation filter in our work for a fair comparison.
In (18), $u_{0}$ is the linear controller given by the Linear–Quadratic
Regulator (LQR) solution using the nominal dynamics (linearized at $x=0$).
Furthermore, it is to be noted that our algorithm can still work without any
initial stabilizing controller. In this case, we need to ensure that the
training set always has a sufficient amount of data during training, and it
could take a longer time to stabilize the learning in the beginning.

Figure 1: Loose saturation function.

We test our approach on three examples: an inverted pendulum, a third-order
strict-feedback system, and a cart-pole system. The hyperparameters and
structures of neural networks are listed in Tables I and II. The experiment
results are reported in Table III along with the baseline given by the
estimated RoA by the LQR solution. The LQR Lyapunov function
$V_{\textrm{LQR}}=x^{\top}Px$ is obtained using the nominal dynamics before
training. Then, the estimation RoA given by $V_{\textrm{LQR}}$ is
$\mathcal{V}(c^{\prime})$, where $c^{\prime}=\max_{c>0}c$ subject to
$\hat{\dot{V}}_{\textrm{LQR}}(x,u_{\theta_{4}}(x))\leq-\kappa\lVert
x\rVert^{2}$ for all $x\in\mathcal{X}_{\text{stable}}\cap\mathcal{V}(c)$ and
$\mathcal{V}(c)\subset\mathcal{X}$. Our method enlarges the true RoA by 740%,
219%, and 173% on the three examples and increases the estimated RoA by at
least 200% compared with the baseline. For examples 1 to 3, the average
training times for each iteration (as in Algorithm 1) are 147, 90, and 109
$\mathrm{s}$ respectively on an i7-5930K CPU. The sample sizes are 10000,
15625, and 10000, respectively.

TABLE I: Hyperparameters used for training (a dash means that the level multiplier $\eta_{i}$ is kept fixed during training). For all examples, $\gamma=10^{-6}$, $\kappa=0.1$, $\epsilon=0.01$. |  Inverted Pendulum |  Strict Feedback Form |  Cart-pole
---|---|---|---
$\lambda_{\text{RoA}}$ | 1000 | 500 | 500
$\lambda_{\text{Lip}}$ | 0.1 | 0.01 | 0.01
$\eta_{0}$ | 5 | 2 | 9
$k_{\eta}$ | 15 | - | -
$a$ | -2 | -1 | -5
$b$ | 2 | 1 | 5
TABLE II: Network structures and activation functions (id stands for the identity mapping). |  Inverted Pendulum |  Strict Feedback Form |  Cart-pole
---|---|---|---
$f_{\theta_{1}}$ |  [16,16,16,1] [tanh,tanh,tanh,id ] |  [16,16,16,3] [tanh,tanh,tanh,id] |  [16,16,16,2] [tanh,tanh,tanh,id]
$g_{\theta_{2}}$ | scalar | scalar |  [16,16,16,2] [tanh,tanh,tanh,id]
$\phi$ (in $V_{\theta_{3}}$) |  [64,64,64] [tanh,tanh,tanh] |  [64,64,64] [tanh,tanh,tanh] |  [64,64,64] [tanh,tanh,tanh]
$\psi$ (in $u_{\theta_{4}}$) |  [16,16,16,1] [tanh,tanh,tanh,id] |  [16,16,16,1] [tanh,tanh,tanh,id] |  [16,16,16,1] [tanh,tanh,tanh,id]
TABLE III: Percentages of RoA, forward invariant RoA, and estimated RoAs over
the state space.

|  Inverted Pendulum |  Strict Feedback Form |  Cart-pole
---|---|---|---
True RoA* (before/after training) | 11.9/100 | 31.1/99.0 | 27.8/76.0
Forward invariant RoA** (before/after training) | 11.9/81.7 | 30.6/70.8 | 16.7/19.2
Estimated RoA (ours) | 37.4 | 16.3 | 1.2
Estimated RoA (LQR) | 9.6 | 0.8 | 0.4

  * *

The percentages presented in Table III are approximate estimations based on a
mesh. For example, the percentage of the true RoA is defined as the ratio of
the number of stable mesh points over the total number of mesh points.

  * **

Forward invariant RoA refers to the set of initial states from which the
trajectory never leaves $\mathcal{X}$ before converging to the origin.

###  V-A Stationary Inverted Pendulum

The stationary inverted pendulum is a second-order nonlinear system with
dynamics:

$ml^{2}\ddot{\theta}-mgl\sin\theta=u$

where $\theta=0$ is the upward position of the pendulum, and $\theta$ is
positive in the counter-clockwise direction. The states are $(\theta,\omega)$
with $\omega=\dot{\theta}$ and the state space is

$\mathcal{X}=[-\pi,\pi]\times[-\pi,\pi].$

The true parameters of the system are $m=1$ $\mathrm{k}\mathrm{g}$ and $l=0.5$
$\mathrm{m}$, while the nominal parameters are $m^{\prime}=0.8$
$\mathrm{k}\mathrm{g}$ and $l^{\prime}=0.4$ $\mathrm{m}$. Thus, our initial
knowledge of the system is not accurate. When determining the output
dimensions of $f_{\theta_{1}}$ and $g_{\theta_{2}}$, we have taken some known
relationships into account, _e.g_., $\dot{\theta}=\omega$ for the inverted
pendulum. Therefore, for this example, the output dimensions of
$f_{\theta_{1}}$ and $g_{\theta_{2}}$ should both be one. The same reasoning
has been applied to the other example.

Figure 2: Estimated RoA by NLR [5] without accounting for the dynamic
uncertainty during training. Left: true RoA (green), false RoA based on the
nominal model (gray), and estimated RoA (blue contour). Right: phase plot and
four divergent trajectories starting from within the estimated RoA.

Figure 3: Left: true RoA and estimated RoA ratios. Right: phase plot and 20
randomly sampled trajectories starting from the boundary of the estimated RoA.

Figure 4: Evolution of the estimated RoA, sampling area $\mathcal{S}_{i}$,
forward invariant RoA, and true RoA (see legend in the first figure). The
number of iterations is 0, 25, 50, 75, 100, and 200 from top left to bottom
right.

Figure 2 demonstrates an overestimation of the true RoA generated by the
offline method called Neural Lyapunov Redesign (NLR, [5]), which does not
consider uncertainties in the system dynamics. Figure 3 shows the percentages
of RoA, forward invariant RoA, and estimated RoA over the entire state space.
The enlargement of the true RoA and the estimated RoA based on our method are
reflected in Fig. 4.

Next, we consider using an SMT solver to verify the CLF. Due to the SMT
solver’s limit on the number of parameters, we consider a simpler setting for
(IV-C). We reduce $f_{\theta_{1}}$, $\phi$, and $\psi$ to neural networks of
three, two, and three layers, respectively. The nonlinear constraint (IV-C) is
solved with $\zeta=0.3$ and a precision of $10^{-3}$, and no counter-example
is found. This means that (13) is verified on the estimated RoA, which is
around $34.6\%$.

###  V-B A Third-Order Strict Feedback Form

Consider a third-order system of the strict feedback form:

$\dot{x}_{1}=e_{1}x_{2},\quad\dot{x}_{2}=e_{2}x_{3},\quad\dot{x}_{3}=e_{3}x_{1}^{2}+e_{4}u.$

The states are $(x_{1},x_{2},x_{3})$ and the state space is

$\mathcal{X}=\\{\,(x_{1},x_{2},x_{3})\;|\;\lvert x_{1}\rvert\leq 1.5,\lvert
x_{2}\rvert\leq 1.5,\lvert x_{3}\rvert\leq 2\,\\}.$

Again, we assume that our initial knowledge of the system is not accurate: the
true parameters are $e_{1}=1$, $e_{2}=1$, $e_{3}=1$, and $e_{4}=1$, while the
nominal parameters are $e_{1}^{\prime}=0.9$, $e_{2}^{\prime}=0.8$,
$e_{3}^{\prime}=0.9$, and $e_{4}^{\prime}=0.8$. Figure 5 shows the evolution
of the percentages of RoA, forward invariant RoA, and estimated RoA along with
a 3D visualization. We note that as our learning framework tries to stabilize
more states, there can be transition periods where the size of the RoA drops,
which might result from the sensitivity of the nonlinear controller to its
weights. Figure 6 demonstrates ten randomly sampled trajectories starting from
the boundary of the estimated RoA.

Figure 5: Left: true RoA and estimated RoA ratios. Right: true RoA (green),
forward invariant RoA (red), and estimated RoA (blue).

Figure 6: Ten randomly sampled trajectories starting from the boundary of the
estimated RoA.

###  V-C Cart-Pole System

Finally, our method is tested on the cart-pole system:

$\displaystyle(M+m)\ddot{x}-ml\ddot{\theta}\cos{\theta}+ml\dot{\theta}^{2}\sin{\theta}+b_{c}\dot{x}$
$\displaystyle=u,$ $\displaystyle ml^{2}\ddot{\theta}-mgl\sin{\theta}$
$\displaystyle=ml\ddot{x}\cos{\theta}.$

The states are $(\theta,\omega,x,v)$, where $\omega=\dot{\theta}$ and
$v=\dot{x}$. The state space is defined as

$\mathcal{X}=\\{\,(\theta,\omega,x,v)\;|\;\lvert\theta\rvert\leq\pi/6,\lvert\omega\rvert\leq
1,\lvert x\rvert\leq 1,\lvert v\rvert\leq 1.5\,\\}.$

Note that $\theta=0$ corresponds to the upward position of the pole, and
$\theta$ is positive in the counter-clockwise direction. The true parameters
are $M=1$ $\mathrm{k}\mathrm{g}$, $m=0.3$ $\mathrm{k}\mathrm{g}$, $l=1$
$\mathrm{m}$, and $b_{c}=0$ $\mathrm{k}\mathrm{g}\,\mathrm{s}^{-1}$, while the
nominal parameters are $M^{\prime}=0.8$ $\mathrm{k}\mathrm{g}$,
$m^{\prime}=0.27$ $\mathrm{k}\mathrm{g}$, $l^{\prime}=0.8$ $\mathrm{m}$, and
$b_{c}^{\prime}=0$ $\mathrm{k}\mathrm{g}\,\mathrm{s}^{-1}$. Figure 7 shows the
evolution of the percentages of RoA, forward invariant RoA, and estimated RoA.
Ten randomly sampled trajectories starting from the boundary of the estimated
RoA are plotted in Fig. 8. The offline method NLR ([5]) generates an estimated
forward invariant RoA ratio of 1.7% using the true parameters, which
demonstrates the difficulty of this task. Finally, we test the robustness of
the RoA estimation to sudden perturbations, where $b_{c}$ is increased to
$9.1$ $\mathrm{k}\mathrm{g}\,\mathrm{s}^{-1}$. As seen in Fig. 9, our
estimated RoA remains valid while the NLR estimated RoA no longer holds.

Figure 7: evolution of the percentages of RoA, forward invariant RoA, and
estimated RoA.

Figure 8: Ten randomly sampled trajectories starting from the boundary of the
estimated RoA.

Figure 9: Ten randomly sampled trajectories (some visually overlapped in the
plots) starting from boundaries of RoAs estimated by our method (left) and NLR
(right), respectively.

##  VI Conclusion

We proposed a learning framework that can synthesize state-feedback
controllers and CLF for control-affine nonlinear systems with unstructured
uncertainties. Our approach initializes the system at different initial
conditions and observes the system trajectories. Exact knowledge of system
dynamics is not required. Based on a regularity condition, we model
uncertainties as bounded and structured disturbances. Experiments show that
our method can find a controller that enlarges the RoA and a CLF that
estimates the RoA. In addition, our work produces RoA estimations with the
uncertainties considered to avoid overestimation (Section V-A) and has better
robustness to sudden changes in dynamics (Section V-C). Future work includes
relaxing assumptions, generalizing our method to higher-dimension systems, and
exploring further combinations of the proposed approach and neural network
verifying tools (such as SMT solvers).

## References

  * [1] E. D. Sontag, “A ‘universal’construction of Artstein’s theorem on nonlinear stabilization,” _Systems & Control Lett._, vol. 13, no. 2, pp. 117–123, 1989.
  * [2] S. M. Richards, F. Berkenkamp, and A. Krause, “The Lyapunov neural network: adaptive stability certification for safe learning of dynamical systems,” in _Proc. Conf. on Robot Learning_ , (Zürich, Switzerland), Oct. 2018, pp. 466–476.
  * [3] B. Dai, P. Krishnamurthy, A. Papanicolaou, and F. Khorrami, “State constrained stochastic optimal control using LSTMs,” in _Proc. American Control Conf._ , (Virtual), May 2021, pp. 1294–1299.
  * [4] B. Dai, P. Krishnamurthy, and F. Khorrami, “Learning a better control barrier function,” in _Proc. IEEE Conf. on Decision and Control_ , (Cancun, Mexico), Dec. 2022, pp. 945–950.
  * [5] A. Mehrjou, M. Ghavamzadeh, and B. Schölkopf, “Neural Lyapunov redesign,” in _Proc. Annual Learning for Dynamics & Control Conf._, (Virtual), June 2021, pp. 459–470.
  * [6] S. Wei, X. Chen, X. Zhang, and C. Qi, “Towards safe and socially compliant map-less navigation by leveraging prior demonstrations,” in _Proc. International Conf. on Intelligent Robotics and Applications_ , (Virtual), Nov 2020, pp. 133–145.
  * [7] A. M. Lyapunov, “The general problem of the stability of motion,” _International J. Control_ , vol. 55, no. 3, pp. 531–534, 1992.
  * [8] Z. Artstein, “Stabilization with relaxed controls,” _Nonlinear Analysis: Theory Methods & Appl._, vol. 7, no. 11, pp. 1163–1173, 1983.
  * [9] E. D. Sontag and Y. Wang, “On characterizations of input-to-state stability with respect to compact sets,” in _Nonlinear Control Systems Design 1995_. Pergamon, 1995, pp. 203–208.
  * [10] P. Krishnamurthy and F. Khorrami, “On uniform solvability of parameter‐dependent Lyapunov inequalities and applications to various problems,” _SIAM Journal on Control and Optimization_ , vol. 45, no. 4, pp. 1147–1164, 2006.
  * [11] M. Krstić and P. V. Kokotović, “Control Lyapunov functions for adaptive nonlinear stabilization,” _Systems & Control Lett._, vol. 26, no. 1, pp. 17–23, 1995.
  * [12] F. Berkenkamp, R. Moriconi, A. P. Schoellig, and A. Krause, “Safe learning of regions of attraction for uncertain, nonlinear systems with Gaussian processes,” in _Proc. IEEE Conf. on Decision and Control_ , (Las Vegas, NV), Dec. 2016, pp. 4661–4666.
  * [13] R. Sun, M. L. Greene, D. M. Le, Z. I. Bell, G. Chowdhary, and W. E. Dixon, “Lyapunov-based real-time and iterative adjustment of deep neural networks,” _IEEE Control Syst. Lett._ , vol. 6, pp. 193–198, 2021.
  * [14] A. J. Taylor, V. D. Dorobantu, M. Krishnamoorthy, H. M. Le, Y. Yue, and A. D. Ames, “A control Lyapunov perspective on episodic learning via projection to state stability,” in _Proc. IEEE Conf. on Decision and Control_ , (Nice, France), Dec. 2019, pp. 1448–1455.
  * [15] B. Tibken, “Estimation of the domain of attraction for polynomial systems via LMIs,” in _Proc. IEEE Conf. on Decision and Control_ , (Sydney, Australia), Dec. 2000, pp. 3860–3864.
  * [16] A. Papachristodoulou and S. Prajna, “On the construction of Lyapunov functions using the sum of squares decomposition,” in _Proc. IEEE Conf. on Decision and Control_ , (Las Vegas, Nevada), Dec. 2002, pp. 3482–3487.
  * [17] P. Giesl and S. Hafstein, “Review on computational methods for Lyapunov functions,” _Discrete and Continuous Dynamical Systems - B_ , vol. 20, no. 8, pp. 2291–2331, 2015.
  * [18] Y.-C. Chang, N. Roohi, and S. Gao, “Neural Lyapunov control,” in _Proc. Conf. on Neural Information Processing Systems_ , (Vancouver, Canada), Dec. 2019.
  * [19] T. Zhao, J. Wang, X. Lu, and Y. Du, “Neural Lyapunov control for power system transient stability: A deep learning-based approach,” _IEEE Trans. on Power Systems_ , vol. 37, no. 2, pp. 955–966, 2021.
  * [20] E. D. Sontag _et al._ , “Smooth stabilization implies coprime factorization,” _IEEE Trans. Autom. Control_ , vol. 34, no. 4, pp. 435–443, 1989.
  * [21] S. Gao, J. Avigad, and E. M. Clarke, “$\delta$-complete decision procedures for satisfiability over the reals,” in _International Joint Conf. on Automated Reasoning_ , (Manchester, UK), June 2012, pp. 286–300.