Safe Approximate Dynamic Programming Via Kernelized Lipschitz Estimation

TL;DR

This paper introduces a method to compute safe initial policies for reinforcement learning in uncertain systems using kernelized Lipschitz estimation and semidefinite programming, ensuring safety and stability.

Contribution

It presents a novel approach combining kernelized Lipschitz estimation with semidefinite programming to generate safe, stabilizing initial policies for approximate dynamic programming.

Findings

Provides high-probability guarantees for admissible policies.

Ensures exponential stability of the closed-loop system.

Enables safe initialization in reinforcement learning.

Abstract

We develop a method for obtaining safe initial policies for reinforcement learning via approximate dynamic programming (ADP) techniques for uncertain systems evolving with discrete-time dynamics. We employ kernelized Lipschitz estimation and semidefinite programming for computing admissible initial control policies with provably high probability. Such admissible controllers enable safe initialization and constraint enforcement while providing exponential stability of the equilibrium of the closed-loop system.

Tables

Table 1. TABLE I: Kernelized Lipschitz Learning of benchmark functions

$n$	${\log }_{10}\beta$	Time [s]	${\widehat{𝔏}}_{\varphi }$ (mean $\pm$ stdev)	OE‡?
${\varphi }_1=\left|\cos \left(\pi x\right)\right|$, ${𝔏}_{\varphi }^{\ast }=3.141$ on $[-\pi ,\pi ]$
$100$	$-2$	0.596 $\pm$ 0.127	3.361 $\pm$ 0.093	$✓$
$100$	$-4$	0.610 $\pm$ 0.122	3.528 $\pm$ 0.177	$✓$
$500$	$-2$	2.463 $\pm$ 0.432	3.252 $\pm$ 0.083	$✓$
$500$	$-4$	2.438 $\pm$ 0.427	3.369 $\pm$ 0.158	$✓$
${\varphi }_2=x-x^3/3$, ${𝔏}_{\varphi }^{\ast }=1.000$ on $[-1,1]$
$100$	$-2$	0.455 $\pm$ 0.123	1.018 $\pm$ 0.011	$✓$
$100$	$-4$	0.459 $\pm$ 0.114	1.030 $\pm$ 0.020	$✓$
$500$	$-2$	1.656 $\pm$ 0.218	1.005 $\pm$ 0.004	$✓$
$500$	$-4$	1.585 $\pm$ 0.208	1.010 $\pm$ 0.009	$✓$
${\varphi }_3=\sin \left(x\right)+\sin \left(2x/3\right)$, ${𝔏}_{\varphi }^{\ast }=1.667$ on $[3.1,20.4]$
$100$	$-2$	0.556 $\pm$ 0.121	1.780 $\pm$ 0.074	$✓$
$100$	$-4$	0.547 $\pm$ 0.124	1.923 $\pm$ 0.166	$✓$
$500$	$-2$	1.826 $\pm$ 0.224	1.684 $\pm$ 0.010	$✓$
$500$	$-4$	1.821 $\pm$ 0.221	1.720 $\pm$ 0.002	$✓$
${\varphi }_4=$ Hansen test function from [29], ${𝔏}_{\varphi }^{\ast }=8.378$ on $[0,1]$
$100$	$-2$	0.450 $\pm$ 0.117	8.969 $\pm$ 0.262	$✓$
$100$	$-4$	0.488 $\pm$ 0.123	9.401 $\pm$ 0.476	$✓$
$500$	$-2$	1.921 $\pm$ 0.138	8.507 $\pm$ 0.046	$✓$
$500$	$-4$	1.923 $\pm$ 0.130	8.707 $\pm$ 0.103	$✓$
${\varphi }_5=\max \{1-3\sin \left(x\right),\exp \left(-\sin \left(x\right)\right)\}$, ${𝔏}_{\varphi }^{\ast }=3.0$ on $[-10,10]$
$100$	$-2$	0.546 $\pm$ 0.061	3.139 $\pm$ 0.051	$✓$
$100$	$-4$	0.612 $\pm$ 0.066	3.200 $\pm$ 0.087	$✓$
$500$	$-2$	1.893 $\pm$ 0.245	3.043 $\pm$ 0.014	$✓$
$500$	$-4$	1.989 $\pm$ 0.230	3.104 $\pm$ 0.024	$✓$

Equations

$$x_{t+1}$$

$$q_t$$

$$x_{t+1}$$

$$q_t$$

$$\Vert f(x_1)-f(x_2)\Vert \le {\mathfrak{L}}_f\Vert x_1-x_2\Vert$$

$${\mathfrak{L}}_f^{\ast }=\underset{{\mathbb{R}}_{+}}{\inf }\{{\mathfrak{L}}_f:\text{condition eqrefeq:lipsz holds}\}$$

$$\Vert \varphi (q_1)-\varphi (q_2)\Vert \le {\mathfrak{L}}_{\varphi }^{\ast }\Vert q_1-q_2\Vert$$

$$\mathcal{J}(x_0,u)=\sum _{t=0}^{\infty }{\gamma }^t\mathcal{U}(x_t,u(x_t)),$$

$$x_{t+1}=A{x}_t+Bu(x_t)+G\varphi \left(C_q{x}_t\right).$$

$$\mathcal{J}_{\infty}(x_{0})=\inf_{u\in\mathfrak{U}}\mathcal{J}\big{(}x_{0},u\big{)},$$

$$u_{\infty}=\operatorname*{arg\,inf}_{u\in\mathfrak{U}}\mathcal{J}\big{(}x_{0},u\big{)}.$$

$$\varphi (q_j)=G^{†}\left(x_j^{+}-A{x}_j-B{u}_j\right),$$

$${\ell }_{jk}=\frac{|\varphi (q_j)-\varphi (q_k)|}{\Vert q_j-q_k\Vert },$$

$${\hat{L}}_n(\ell )=\frac{1}{n{h}_n}\sum _{r=1}^n\mathcal{K}\left(\frac{\ell -{\ell }_r}{h_n}\right),$$

$${\hat{S}}_n:=\{\ell \in {\mathbb{R}}_{\ge 0}:{\hat{L}}_n(\ell )\ge {\beta }_n\},$$

$${\hat{\mathfrak{L}}}_{\varphi }:=\max ({\hat{S}}_n).$$

$$\tilde{L}(\lambda )={\mu }_X\left(\left\{(q_1,q_2):q_1\ne q_2,\frac{|\varphi (q_1)-\varphi (q_2)|}{\Vert q_1-q_2\Vert }\le \lambda \right\}\right).$$

$$\phi (t)=\varphi \left(q_1^{\infty }+t(q_2^{\infty }-q_1^{\infty })\right)-\varphi (q_1^{\infty })$$

$$\mathfrak{L}^{*}_{\phi}\|q_{2}^{\infty}-q_{1}^{\infty}\|=\varphi(1)=\int_{0}^{1}\big{(}\mathrm{d}\varphi/\mathrm{d}t\big{)}\,\mathrm{d}s\leq\mathfrak{L}^{*}_{\phi}\|q_{2}^{\infty}-q_{1}^{\infty}\|.$$

$${\varphi }^{\prime }(q_{+},q_{-})\ge (1-\delta ){\mathcal{L}}_{\varphi }^{\ast }-C_0\delta$$

$$\varphi (q^{\star }+q)=\varphi (q^{\star })+\nabla \varphi (q^{\star }{)}^{\top }q+\mathfrak{R}(q),$$

$$\nabla \varphi (q^{\star }{)}^{\top }q=\Vert \nabla \varphi (q^{\star })\Vert \Vert q\Vert \cos \theta ,$$

$$|\varphi (q^{\star }+q_{+})-\varphi (q^{\star }+q_{-})|$$

$$=|\nabla \varphi (q^{\star }{)}^T(q_{+}-q_{-})+\mathfrak{R}(q_{+})-\mathfrak{R}(q_{-})|$$

$$\ge \Vert \nabla \varphi (q^{\star })\Vert \Vert q_{+}-q_{-}\Vert \left(1-\frac{{\chi }^2}{2}\right)-|\mathfrak{R}(q_{+})-\mathfrak{R}(q_{-})|.$$

$${\varphi }^{\prime }(q_{+},q_{-})\ge \left(1-\frac{{\chi }^2}{2}\right){\mathcal{L}}_{\varphi }^{\ast }-\frac{|\mathfrak{R}(q_{+})-\mathfrak{R}(q_{-})|}{\Vert q_{+}-q_{-}\Vert }.$$

$$\frac{|\mathfrak{R}(q_{+})-\mathfrak{R}(q_{-})|}{\Vert q_{+}-q_{-}\Vert }$$

$$\le \frac{2C_{\mathfrak{R}}\max \{\Vert q_{+}\Vert ,\Vert q_{-}\Vert {\}}^2}{\max \{\Vert q_{+}\Vert ,\Vert q_{-}\Vert \}(1-\chi )\left(1-{\chi }^2/2\right)}$$

$$\le \frac{2C_{\mathfrak{R}}\eta (1+{\chi }^2{)}^{1/2}}{(1-\chi )\left(1-{\chi }^2/2\right)}.$$

$${\phi }^{\prime }(q_{+},q_{-})\ge \left(1-\frac{\chi }{2}\right){\mathfrak{L}}_{\varphi }^{\ast }-\frac{2C_{\mathfrak{R}}\eta (1+{\chi }^2{)}^{1/2}}{(1-\chi )\left(1-{\chi }^2/2\right)}.$$

$${\phi }^{\prime }(q_{+},q_{-})\ge \left(1-\delta \right){\mathfrak{L}}_{\varphi }^{\ast }-C_0\delta ,$$

$${\int }_0^{\eta }\kappa (\chi r{)}^{n_q-1}\mathrm{d}r=\frac{\kappa {\chi }^{n_q-1}{\eta }^{n_q}}{n_q}$$

$$1$$

Full text

I Introduction

Recent advances in the field of deep and machine learning has led to a renewed interest in using learning for control of physical systems [1]. Reinforcement learning (RL) is a learning framework that handles sequential decision-making problems, wherein an ‘agent’ or decision maker learns a policy to optimize a long-term reward by interacting with the (unknown) environment. At each step, an RL agent obtains evaluative feedback (called reward or cost) about the performance of its action, allowing it to improve the performance of subsequent actions [2, 3]. While RL has witnessed huge success in recent times [4, 5], there are several unsolved challenges which restricts use of these algorithms for industrial systems. In most practical applications, control policies must be designed to satisfy operational constraints. This leads to the challenge that one has to guarantee constraint satisfaction during learning and policy optimization. Therefore, initializing with an unverified control policy is not ‘safe’ (in terms of stability or constraint handling). In other words, using on-line RL for expensive equipment or safety-critical applications necessitates that the initial policy used for obtaining data for subsequently improved policies must be at least stabilizing, and generally, constraint-enforcing. The work presented in this paper is motivated by this challenge. We present a framework for deriving initial control policies from historical data that can be verified to satisfy constraints and guarantee stability while learning the optimal control policy on-line, from operational data.

A successful RL method needs to balance a fundamental trade-off between exploration and exploitation. One needs to gather data safely (exploration) in order to best extract information from this data for optimal decision-making (exploitation). One way to solve the exploration and exploitation dilemma is to use optimistic initialization [6, 7, 8, 9], but this assumes the optimal policy is available until data is obtained that proves otherwise. Such approaches have been applied to robotics applications, where systems with discrete and continuous state-action spaces [10, 11]. A limitation of these methods is that, before the optimal policy is learned, the agent is quite likely to explore actions that lead to violation of the task-specific constraints as it aims to optimize the cumulative reward for the task. This shortcoming significantly limits such methods to be applicable to industrial applications, since this could lead to irreparable hardware damage or harm human operators due to unexpected dynamics. Consequently, safe learning focuses on learning while enforcing safety constraints. There are primarily two types of approaches to safe RL and approximate/adaptive dynamic programming (ADP). These include: modification of the optimization criterion with a safety component such as barrier functions by transforming the operational constraints into soft constraints [12, 13]; and, modifying the exploration process through the incorporation of external system knowledge or historical data [14]. Our method is amongst the latter class of methods, because our operational constraints are hard constraints and softening them could lead to intermittent failure modes.

High performance model-based control requires precise model knowledge for controller design. However, it is well known that for most applications, accurate model knowledge is practically elusive due to the presence of unmodeled dynamical interactions (e.g., friction, contacts, etc.). Recent efforts tackle this issue by learning control policies from operational (on-line) or archival data (off-line). Since the exact structure of the nonlinearity may be unknown or not amenable for analysis, researchers have proposed ‘indirect’ data-driven controllers that employ non-parametric learning methods such as Gaussian processes to construct models from operational data [15, 16] to improve control policies on-line [17, 18]. Conversely, ‘direct’ methods, such as those proposed in [19, 20, 21, 22], directly compute policies using a combination of archival/legacy and operational input-output data without constructing an intermediate model. For example, in [23], a human expert was introduced into the control loop to conduct initial experiments to ensure safety while generating archival data. A common assumption in many of these approaches is the availability of an initial control policy that is stabilizing and robust to unmodeled dynamics. Designing such safe initial control policies in a computationally tractable manner remains an open challenge.

In this work, we present a formalism for synthesizing safe initial policies for uncertain non-linear systems. We assume the presence of historical/archival/legacy data, with which we estimate Lipschitz constants for the unmodeled system dynamics. The estimation of the Lipschitz constant is done via kernel density estimation (KDE). The estimated Lipschitz constant is used to design control policies via semidefinite programming that can incorporate stability and constraint satisfaction while searching for policies. We show that the proposed approach is able to design feasible policies for different constrained tasks for several systems while respecting all active constraints. Our key insight is that information regarding the structure of classes of unmodeled nonlinearities can be encapsulated using only a few parameters, without knowing the exact form of the nonlinearity. Therefore, it may not be necessary to model the unknown component itself in order to compute a safe control policy. For example, the class of Lipschitz nonlinearities (which constitute a large share of nonlinearities observed in applications) can be described using only a few parameters: the Lipschitz constants of the nonlinear components. Recent work has investigated the utility of Lipschitz properties in constructing controllers when an oracle is available [24] or in designing models for prediction [25] with on-line data used for controller refinement [26] In this paper, we construct control policies that respect constraints and certify stability (with high probability) for applications where only off-line data is available, and no oracle is present. We do so through the systematic use of multiplier matrices that enable the representation of nonlinear dynamics through quadratic constraints [27, 28] without requiring knowledge of the underlying nonlinearity. The control policies can then be obtained by solving semidefinite programs. However, construction of multiplier matrices for Lipschitz systems requires knowledge of the Lipschitz constants, which are not always available, and therefore, must be estimated. We refer to the estimation of Lipschitz constants from data as Lipschitz learning. Historically, methods that estimate the Lipschitz constant [29, 30, 31] do not provide certificates on the quality of the estimate. Herein, we provide conditions that, if satisfied, enable us to estimate the Lipschitz constant of an unknown locally Lipschitz nonlinearity with high probability. To this end, we employ kernel density estimation (KDE), a non-parametric data-driven method that employs kernels to approximate smooth probability density functions to arbitrarily high accuracy. We refer to our proposed KDE-based Lipschitz constant estimation algorithm as kernelized Lipschitz learning.

Contributions

Compared to the existing literature on safe learning, the contributions of the present paper are threefold. First, we formulate an algorithm to construct stabilizing and constraint satisfying policies for nonlinear systems without knowing the exact form of the nonlinearity. Then we leverage a kernelized Lipschitz learning mechanism to estimate Lipschitz constants of the unmodeled dynamics with high probability; and, finally we use a multiplier-matrix based controller design based on Lipschitz learning from legacy data that forces exponential stability on the closed-loop dynamics (with the same probability as the kernelized Lipschitz learner).

Structure

The rest of the paper is structured as follows. We present the formal motivation of our work in Section II. Our kernelized Lipschitz learning algorithm is described in Section III, and benchmarking of the proposed learner on benchmark Lipschitz functions is performed. The utility of Lipschitz learning in policy design via multiplier matrices is elucidated in Section IV, and a numerical example demonstrating the potential of our overall formalism is provided in Section V. We provide concluding remarks and discuss future directions in Section VI.

Notation

We denote by $\mathbb{R}$ the set of real numbers, $\mathbb{R}_{+}$ as the set of positive reals, and $\mathbb{N}$ as the set of natural numbers. The measure-based distance between two measurable subsets $A$ and $B$ of a metric space $\mathbb{R}^{n}$ equipped with the metric $\rho_{\mu}$ is given by $\rho_{\mu}(A,B)=\mu(A\triangle B)$, where $\mu$ is a measure on $\mathbb{R}^{n}$ and $A\triangle B$ is the symmetric difference $(A\setminus B)\cup(B\setminus A)$. We define a ball $\mathcal{B}_{\epsilon}(x):=\{y:\rho(x,y)\leq\epsilon\}$ and the sum $A\oplus\epsilon:=\bigcup_{x\in A}\mathcal{B}_{\epsilon}(x)$. The complement of a set $A$ is denoted by $A^{c}$. The indicator function of the set $A$ is denoted by $\mathbf{1}_{A}$. A block diagonal matrix is denoted by \mathrm{blkdiag}\big{(}\cdot\big{)}. For every $v\in\mathbb{R}^{n}$, we denote $\|v\|=\sqrt{v^{\top}v}$, where $v^{\top}$ is the transpose of $v$. The sup-norm or $\infty$-norm is defined as $\|v\|_{\infty}\triangleq\sup_{t\in\mathbb{R}}\|v(t)\|$. We denote by $\lambda_{\min}(P)$ and $\lambda_{\max}(P)$ as the smallest and largest eigenvalue of a square, symmetric matrix $P$. The symbol $\succ(\prec)$ indicates positive (negative) definiteness and $A\succ B$ implies $A-B\succ 0$ for $A,B$ of appropriate dimensions. Similarly, $\succeq(\preceq)$ implies positive (negative) semi-definiteness. The operator norm is denoted $\|P\|$ and is defined as the maximum singular value of $P$. For a symmetric matrix, we use the $\star$ notation to imply symmetric terms, that is, $\left[\begin{smallmatrix}a&b\\ b^{\top}&c\end{smallmatrix}\right]\equiv\left[\begin{smallmatrix}a&b\\ \star&c\end{smallmatrix}\right]$. The symbol $\mathbf{Pr}$ denotes the probability measure.

II Problem Formulation

II-A Problem statement

Consider the following discrete-time nonlinear system,

[TABLE]

where $x_{t}\in\mathbb{R}^{n_{x}},\ u=u_{t}\in\mathbb{R}^{n_{u}}$ denote the state and the control input of the system respectively.

For simplicity of exposition we will write

[TABLE]

where the system matrices $A$, $B$, $G$ and $C_{q}$ have appropriate dimensions. Denote by $\phi\in\mathbb{R}^{n_{\phi}}$ the system’s uncertainty, or unmodeled nonlinearity, whose argument $q=q_{t}\in\mathbb{R}^{n_{q}}$ is represented by a linear combination of the state. The origin is an equilibrium state for (1); that is, $\phi(0)=0$.

The following assumptions and definition are now needed.

Assumption 1**.**

The matrix $B$ is known. The matrix $G$ has full column rank and is sparse and all entries are 0 or 1. Only the non-zero element locations of $G$ are known. The matrix $C_{q}$ and function $\phi$ are unknown. $\Box$

We require the following definition to describe the class of nonlinearities considered in this paper.

Definition 1**.**

A function $f:\mathbb{X}\to\mathbb{R}^{n_{x}}$ is Lipschitz continuous in the domain $\mathbb{X}\subset\mathbb{R}^{n_{f}}$ if

[TABLE]

for some $\mathfrak{L}_{f}>0$ and all $x_{1},x_{2}\in\mathbb{X}$. We define the scalar

[TABLE]

as the Lipschitz constant of $f$ in $\mathbb{X}$. A function is globally Lipschitz if (2) holds for $\mathbb{X}\equiv\mathbb{R}^{n_{f}}$.  $\Box$

Assumption 2**.**

The nonlinearity $\phi$ is globally Lipschitz continuous. That is,

[TABLE]

for any $q_{1},q_{2}\in\mathbb{R}^{n_{q}}$, and the global Lipschitz constant $\mathfrak{L}^{*}_{\phi}$ is unknown. $\Box$

Assumptions 1 and 2 imply that the linear component of the true system (1) can be assumed, but the rest is unknown. However, we do know the vector space through which the nonlinearity enters the dynamics of (1), since the non-zero locations of $G$ are flagged.

Remark 1*.*

Assumption 1 is mild. For instance, one could relax the assumption on $G$ and take the unknown $\tilde{G}$ to be the identity matrix. Then the nonlinearity would be $\tilde{\phi}(q)=\begin{bmatrix}0&\ldots&\phi_{i}(q)&\ldots\end{bmatrix}^{\top},\ i\in\mathcal{I}$, with $\mathcal{I}$ the index set of non-zero rows of $G$, so that $\tilde{G}\tilde{\phi}=G\phi$.  $\Box$

Given a control policy $u(x)$, we define an infinite horizon cost functional given an initial state $x_{0}\in\mathbb{R}^{n_{x}}$ as

[TABLE]

where $\mathcal{U}$ is a function with non-negative range, $\mathcal{U}(0,0)=0$, and $\{x_{k}\}$ denotes the sequence of states generated by the closed loop system

[TABLE]

The scalar $\gamma\in(0,1]$ is a forgetting/discount factor intended to enable the cost to be emphasized more by current state and control actions and lend less credence to the past.

Before formally stating our objective, we need to introduce the following standard definition [1].

Definition 2**.**

A continuous control policy $u(\cdot):\mathbb{R}^{n_{x}}\to\mathbb{R}^{n_{u}}$ is admissible on $X\subset\mathbb{R}^{n_{x}}$ if it stabilizes the closed loop system (6) on $X$ and $\mathcal{J}(x_{0},u)$ is finite for any $x_{0}\in X$.

We want to design an optimal control policy that achieves the optimal cost

[TABLE]

for any $x_{0}\in\mathbb{R}^{n_{x}}$. Here, $\mathfrak{U}$ denotes the set of all admissible control policies. In other words, we wish to compute an optimal control policy

[TABLE]

Directly constructing such an optimal controller is very challenging for general nonlinear systems; this is further complicated because the system (1) contains unmodeled/uncertain dynamics. Therefore, we shall use adaptive/approximate dynamic programming (ADP): a class of iterative, data-driven algorithms that generate a convergent sequence of control policies whose limit is provably the optimal control policy $u_{\infty}(x)$.

Recall from [32, 33] that a necessary condition for convergence of policy iteration methods (a sub-class of ADP) is the availability of an initial admissible control policy $u_{0}(x)$, which is non-trivial to derive for systems with some unmodeled dynamics. Therefore, our objective in this work is to systematically derive an initial admissible control policy using only partial model information via kernelized Lipschitz learning and semidefinite programming. We also extend this idea to handle the case when the control input is constrained. In such cases, along with an admissible controller, we also derive a domain of attraction of the controller within which the control policy is guaranteed to satisfy input constraints and the closed-loop system remains stable. We refer to the derivation of admissible control policies with guaranteed stabilizability and/or constraint enforcement as safe initialization for ADP: a crucial property required for ADP algorithms to gain traction in expensive industrial applications.

We invoke the assumption in [20, 21] regarding the availability of legacy/archival/historical data generated by the system during prior experiments. That is, at design time, we have a dataset $\mathcal{D}$ consisting of unique triples: state-input pairs along with corresponding state update information. Concretely, we have access to $\mathcal{D}=\{x_{j},u_{j},x^{+}_{j}\}_{j=1}^{N}$. For each $\{x_{j},u_{j},x^{+}_{j}\}\in\mathcal{D}$, we estimate the nonlinear term using (1); that is,

[TABLE]

where $G^{\dagger}$ exists by Assumption 1. Note that we also need to estimate the matrix $C_{q}$ (see (1)) so that $q_{j}$ can be calculated from $x_{j}$. While estimating the exact elements of these matrices is quite challenging, we can estimate the non-zero elements in the matrices, which is enough to design safe initial control policies, because the exact elements of $C_{q}$ will be subsumed within the Lipschitz constant.

Remark 2*.*

The problem of estimating the sparsity pattern of $C_{q}$ is analogous to the problem of feature selection and sparse learning, known as automatic relevance determination (ARD) [34]. The basic idea in ARD is to give feature weights some parametric prior densities; these densities are subsequently refined by maximizing the likelihood of the data [34, 35]. For example, one can define hyperparameters which explicitly represent the relevance of different inputs to a machine learning algorithm w.r.t. the desired output (e.g., a regression problem). These relevance hyperparameters determine the range of variation of parameters relating to a particular input. ARD can then determine these hyperparameters during learning to discover which inputs are relevant. $\Box$

We need the following assumption on the data $\{q_{j},\phi(q_{j})\}$, without which one cannot attain the global Lipschitz constant of the nonlinearity $\phi(\cdot)$ with high accuracy.

Assumption 3**.**

Let $\mathcal{Q}$ denote the convex hull of the samples $\{q_{j}\}$. The Lipschitz constant of $\phi(\cdot)$ in the domain $\mathcal{Q}$ is identical to the global Lipschitz constant $\mathfrak{L}^{*}_{\phi}$. $\Box$

Assumption 3 ensures that the samples obtained from the archival data are contained in a subregion of $\mathbb{R}^{n_{q}}$ where the nonlinearity $\phi(\cdot)$’s local Lipschitz constant is the same as its global Lipschitz constant.

Example 1**.**

Suppose $\phi(q)=1.5\sin(q)$. As long as the convex hull of the samples $\{q\}$ contain zero, the Lipschitz constant of $\phi$ on the convex hull $\mathcal{Q}$ and on $\mathbb{R}$ are identical.  $\Box$

In the following section, we will leverage the dataset $\mathcal{D}$ to estimate the Lipschitz constant of $\phi(\cdot)$ using kernelized Lipschitz learning/estimation, and consequently design an initial admissible linear control policy $u_{0}=K_{0}x$ via semidefinite programs. We will demonstrate how such an initial admissible linear control policy fits into a neural-network based ADP formulation (such as policy iteration) to asymptotically generate the optimal control policy $u_{\infty}(x)$.

Remark 3*.*

The control algorithm proposed in this paper is a direct data-driven controller because no model of $\phi(\cdot)$ is identified in the controller design step. $\Box$

Remark 4*.*

Although we focus only on discrete-time systems, our results hold for continuous-time systems with slight modifications. $\Box$

Remark 5*.*

If $n_{\phi}>1$, our proposed Lipschitz learning algorithm will yield $n_{\phi}$ Lipschitz constant estimates, one for each dimension of $\phi(\cdot)$. To avoid notational complications, we proceed (without loss of generality) with $n_{\phi}=1$. For larger $n_{\phi}$, our algorithm can be used component-wise. $\Box$

III Kernelized Lipschitz Learning

In this section, we provide a brief overview of kernel density estimation (KDE) and provide a methodology for estimating Lipschitz constants from data.

III-A Empirical density of Lipschitz estimates

With the data $\{\phi(q_{j}),q_{j}\}_{j=1}^{N}$, we obtain $n\in\mathbb{N}$ underestimates of the global Lipschitz constant $\mathfrak{L}^{*}_{\phi}$ using

[TABLE]

where $k\in\{1,\ldots,N\}\setminus j$. The sequence $\{\ell_{jk}\}$ are empirical samples drawn from an underlying univariate distribution $L$. Clearly, the true distribution $L$ has finite support; indeed, its left-hand endpoint is a non-negative scalar (zero, if $n_{q}>1$ but may be positive if $n_{q}=1$) and its right-hand endpoint is $\mathfrak{L}^{*}_{\phi}$. This leads us to the key idea of our approach that is to identify the support of the distribution $L$ to yield an estimate of the true Lipschitz constant of $\phi(\cdot)$.

Remark 6*.*

Variants of the estimator (9) such as $\max_{k}\ell_{jk}$ have been widely used in the literature to construct algorithms for determining Lipschitz constants, see for example: [29, 30, 36]. $\Box$

In the literature, common methods of tackling the support estimation problem is by assuming prior knowledge about the exact density of Lipschitz estimates [36] or using Strongin overestimates of the Lipschitz constant [30]. However, we avoid these overestimators because they are provably unreliable, even for globally Lipschitz functions [31, Theorem 3.1]. Instead, we try to fit the density directly from local estimates and the data in a non-parametric manner using KDE and characteristics of the estimated density.

III-B Plug-in support estimation

With a set of $n$ underestimates $\{\ell_{r}\}_{r=1}^{n}$, we generate an estimate $\hat{L}_{n}$ of the true density $L$ using a kernel density estimator

[TABLE]

where $\mathcal{K}:\mathbb{R}\to\mathbb{R}$ is a smooth function called the kernel function and $h_{n}>0$ is the kernel bandwidth. A plug-in estimate of the support $S$ of the true density $L$ is

[TABLE]

where $\beta_{n}$ is an element of a sequence $\{\beta_{n}\}$ that converges to zero as $n\to\infty$; this plug-in estimator was proposed in [37].

III-C Implementation details

Implementing the plug-in estimator involves first constructing a KDE of $L$ with $n$ samples. Then, if one picks $\beta\equiv\beta_{n}$ small enough, one can easily compute $\hat{S}$ from (11). Then

[TABLE]

This is a very straightforward operation with the availability of tools like ksdensity (MATLAB) and the KernelDensity tool in scikit-learn (Python). The pseudocode is detailed herein in Algorithm 1.

Remark 7*.*

Note that the true support $S$ is a subset of $\mathbb{R}_{\geq 0}$. Therefore, when computing the density estimate, this information should be fed into the tool being used. For example, in MATLAB, one has the option $\{$’support’, ’positive’$\}$. Essentially, this subroutine transforms the data into the log-scale and estimates the log-density so that upon returning to linear scale, one preserves positivity.

III-D Theoretical guarantees

We formally describe the density $L$. We consider that the samples $q\in\mathcal{Q}$ are drawn according to some probability distribution $\mu_{0}$ with support $\mathbb{X}$. For any set $S$, suppose that $\mu_{0}$ can be written as $\mu_{0}(S)=\int_{S}\varOmega(q)\,\mathrm{d}\mu(q)$, where $\mu$ is the Lebesgue measure, and $\varOmega$ is continuous and positive on $\mathbb{X}$. Let $\mu_{X}$ denote the product measure $\mu_{0}\times\mu_{0}$ on $\mathbb{X}\times\mathbb{X}$. Since $\mu_{0}$ is absolutely continuous with respect to the Lebesgue measure, $\mu_{X}$ assigns zero mass on the diagonal $\{(q,q):q\in\mathbb{X}\}$. The cumulative distribution function for $L$ is then given by

[TABLE]

Since $\tilde{L}$ is non-decreasing, $L$ exists almost everywhere by Lebesgue’s theorem for differentiability of monotone functions, and $L$’s support is contained within $[0,\mathfrak{L}^{*}_{\phi}]$ because of (9).

We investigate the worst-case sample complexity involved in overestimating $\mathfrak{L}^{*}_{\phi}$ under the following mild assumption.

Assumption 4**.**

The nonlinearity $\phi(\cdot)$ is twice continuously differentiable, that is, $\phi(\cdot)\in\mathcal{C}^{2}$. $\Box$

Lemma 1**.**

Suppose that Assumptions 3 and 4 hold. Then there exists some $q^{\ast}\in\mathcal{Q}$ such that $\|\nabla\phi(q^{\ast})\|=\mathfrak{L}^{*}_{\phi}$.

Proof.

Suppose $\{(q_{1}^{k},q_{2}^{k})\}_{k=1}^{\infty}$ denotes a sequence of paired samples in $\mathcal{Q}$ such that $|\phi(q_{1}^{k})-\phi(q_{2}^{k})|/\|q_{1}^{k}-q_{2}^{k}\|\to\mathfrak{L}^{*}_{\phi}$ as $k\to\infty$. Since $\mathcal{Q}$ is the convex hull of finitely many samples, it is compact, so we can choose a subsequence of $\{q_{1}^{k},q_{2}^{k}\}_{k=1}^{\infty}$ that converges to $(q_{1}^{\infty},q_{2}^{\infty})$ where both limits are in $\mathcal{Q}$. If $q_{1}^{\infty}=q_{2}^{\infty}$, then a Taylor expansion estimate implies $\|\nabla\phi(q_{1}^{\infty})\|\geq\mathfrak{L}^{*}_{\phi}$. Since $\mathfrak{L}^{*}_{\phi}$ is an upper bound of $\|\nabla\phi\|$ at any sample in $\mathcal{Q}$, $\|\nabla\phi(q_{1}^{\infty})\|=\mathfrak{L}^{*}_{\phi}$ and $q^{\ast}=q_{1}^{\infty}$. If $q_{1}^{\infty}\neq q_{2}^{\infty}$, then the result follows by applying the mean value theorem to

[TABLE]

for $t\in[0,1]$, for which $\varphi(0)=0$ and $\varphi(1)=\mathfrak{L}^{*}_{\phi}\|q_{1}^{\infty}-q_{2}^{\infty}\|$. Also, \mathrm{d}\varphi/\mathrm{d}t=\big{(}\nabla\phi(q_{1}^{\infty}+t(q_{2}^{\infty}-q_{1}^{\infty}))\big{)}^{\top}(q_{2}^{\infty}-q_{1}^{\infty}). Since $\|\nabla\phi\|\leq\mathfrak{L}^{*}_{\phi}$, this implies $|\mathrm{d}\varphi/\mathrm{d}t|\leq\mathfrak{L}^{*}_{\phi}\|q_{2}^{\infty}-q_{1}^{\infty}\|$. Reordering $q_{1}^{\infty}$ and $q_{2}^{\infty}$ if needed, we have

[TABLE]

Hence, the rightmost inequality must be an equality, which implies that $\textrm{d}\varphi(s)/\textrm{d}t=\mathfrak{L}^{*}_{\phi}\|q_{1}^{\infty}-q_{2}^{\infty}\|$ for all $s$. That is, if the Lipschitz constant is attained with $q_{1}^{\infty}\neq q_{2}^{\infty}$, then $\phi(\cdot)$ restricted to the segment connecting $q_{1}^{\infty}$ and $q_{2}^{\infty}$ is linear with slope $\mathfrak{L}^{*}_{\phi}$. This concludes the proof.

Lemma 1 enables the worst-case complexity result described in the following theorem.

Theorem 1**.**

Let $\varphi^{\prime}(q_{1},q_{-1})=|\phi(q_{1})-\phi(q_{-1})|/\|q_{1}-q_{-1}\|$, and suppose that Assumptions 3 and 4 hold. There exists $C_{0}>0$ such that for all sufficiently small $\delta>0$ and and any set $\{q_{j}\}_{j=1}^{n}$ of $n$ uniform random samples in ${\mathbb{X}}$, the probability that some pair $q_{+},q_{-}\in\{q_{j}\}$ gives the Lipschitz estimate

[TABLE]

is at least $1-\epsilon(n,\delta)$. Here $\epsilon(n,\delta)\leq 3\exp(-n\kappa\delta^{2n_{q}-1})$, where $\kappa$ is a constant depending on $n_{q}$.

Proof.

By Lemma 1, there exists at least one $q^{\ast}$ such that $\|\nabla\phi(q^{\star})\|=\mathfrak{L}^{*}_{\phi}$. For the worst-case analysis, suppose this occurs only at a single sample, $q^{\star}$. A Taylor expansion at $q^{\star}$ yields

[TABLE]

where $\mathfrak{R}$ is a remainder term with $|\mathfrak{R}(q)|\leq C_{\mathfrak{R}}\|q\|^{2}$ when $\|q\|\leq\eta$, for some $C_{\mathfrak{R}}>0$ and $\eta>0$. Note that

[TABLE]

where $\theta$ is the angle between $\nabla\phi(q^{\star})$ and $q$. To obtain a good estimate of $\|\nabla\phi(q^{\star})\|$, one needs to sample two points in this neighborhood, one point $q_{+}$ with $\cos\theta\approx 1$ and a second point $q_{-}$ with $\cos\theta\approx-1$. Each of these conditions defines a cone. Regarding one of these cones in cylindrical coordinates $0\leq\ell\leq\eta$ and $\|y\|\leq\chi\ell$ for $\chi=\tan\theta$, we can integrate the $n_{q}-1$ dimensional volume to get the volume of this cone as $C_{0}\chi^{n_{q}-1}\eta^{n_{q}}$ for some dimension-dependent constant $C_{0}$. A calculation shows that $\cos\theta\geq 1-\chi^{2}/2$ for small $\theta$. With $q_{+}$ and $q_{-}$ as one sample from each cone, we have that $q_{+}-q_{-}$ is contained in the cone $\|y\|\leq\chi\ell$, and so we can use (13) to approximate

[TABLE]

Dividing by $\|q_{+}-q_{-}\|$ and using the defining property of $q^{\star}$ gives

[TABLE]

Subsequently, one can decompose $q_{+}=q_{+}^{\parallel}+q_{+}^{\perp}$, with $q_{+}^{\parallel}$ parallel to $\nabla\phi(q^{\star})$ and $q_{+}^{\perp}$ perpendicular and satisfying $\|q_{+}^{\perp}\|\leq\chi|q_{+}^{\parallel}|$. Identical arguments can be used to infer $\|q_{-}^{\perp}\|\leq\chi|q_{-}^{\parallel}|$, and hence, $\|q_{+}-q_{-}\|\geq\|q_{+}^{\|}-q_{-}^{\|}\|-\|q_{+}^{\perp}-q_{-}^{\perp}\|.$ Since $q_{+}$ and $q_{-}$ are chosen from opposite cones, we have $\|q_{+}^{\|}-q_{-}^{\|}\|=\|q_{+}^{\|}\|+\|q_{-}^{\|}\|$. Using $\|q^{\perp}\|\leq\chi\|q^{\|}\|$ and $\cos\theta\geq 1-\chi^{2}/2$, we have $\|q_{+}-q_{-}\|\geq(\|q_{+}\|+\|q_{-}\|)(1-\chi)\left(1-\chi^{2}/2\right)$. Hence,

[TABLE]

By combining the aforementioned inequality with  (14), one obtains

[TABLE]

Set $\delta=\chi/2$ and take $\eta=\delta$. Then there exists $C_{0}>0$ such that for all sufficiently small $\delta$,

[TABLE]

which implies that $\varphi^{\prime}\to\mathfrak{L}^{*}_{\phi}$ as $\delta\to 0$.

From the assumption on uniformly drawn samples, the probability of sampling in one of the $(\chi,\eta)$ cones is,

[TABLE]

for some $\kappa>0$ that depends on $n_{q}$. Using $\delta=\eta=\chi/2$ and absorbing the factors of 2 and $1/n_{q}$ into $\kappa$ yields $\kappa\delta^{2n_{q}-1}$.

Let $\mathfrak{X}_{1\pm}$ be the event of sampling at least one point in the $(\chi,\eta)$ cone as above, and let $\mathfrak{X}_{0\pm}$ be the event of sampling nothing in the $(\chi,\eta)$ cone. The probability of sampling at least one of each of the points $q_{+}$ and $q_{-}$ just described is,

[TABLE]

where the factor 2 before $\kappa$ in the second term comes from the fact that both cones are excluded and they are disjoint, and the 2 before the third term comes by combining the final two terms in the first expression.

By using the fact that for any $\epsilon^{\prime}\in(0,1)$ and $n>0$ the inequality $(1-\epsilon^{\prime})^{n}\leq\exp(-n\epsilon^{\prime})$ holds, we can conclude that,

[TABLE]

for any given $\epsilon>0$. The latter can be ensured by choosing $n$ large enough. This gives a lower bound on the probability of obtaining (15) and hence, the desired result.

III-E Benchmarking the Lipschitz estimator

Our Lipschitz estimator is tested on well-studied benchmark examples studied previously in [29, 25]: the benchmark functions are described in Table I along with their domains and true local Lipschitz constants. Note that all the functions are not globally Lipschitz (e.g. $\phi_{2}$), not differentiable everywhere (e.g. $\phi_{1}$, $\phi_{4}$), and, in the special case of $\phi_{4}$, specifically constructed to ensure that naive overestimation of $\mathfrak{L}^{*}_{\phi}$ using Strongin methods provably fails [31]. To evaluate the proposed Lipschitz estimator, we vary the number of data points $N$ and the confidence parameter $\beta$. Over 100 runs, we report mean $\pm$ one standard deviation of the following quantities: the time required by our learning algorithm, the estimated Lipschitz constant $\hat{\mathfrak{L}}_{\phi}$, and the error $\hat{\mathfrak{L}}_{\phi}-\mathfrak{L}^{*}_{\phi}$ (which should be positive when we overestimate $\mathfrak{L}^{*}_{\phi}$).

The final column of Table I reveals an important empirical detail: all our estimates of $\mathfrak{L}^{*}_{\phi}$ are overestimates for $\beta\leq 0.01$ and $n\geq 100$. This is a critical advantage of our proposed approach, because an overestimate will enable us to provide stability and constraint satisfaction guarantees about the data-driven controller, as we will discuss in subsequent sections. Furthermore, the estimation error is small for every test run, and as expected, the error increases as $\beta$ decreases, because a smaller value of $\beta$ indicates the need for greater confidence, which results in more conservative estimates.

IV Safe Initialization in ADP

In this section, we begin by reviewing a general ADP procedure, and then explain how to safely initialize unconstrained, as well as input-constrained ADP.

IV-A Unconstrained ADP

Recall the optimal value function given by (7) and the optimal control policy (8). From the Bellman optimality principle, we know that the discrete-time Hamilton-Jacobi-Bellman equations are given by

[TABLE]

where $J_{\infty}(x_{t})$ is the optimal value function and $u_{\infty}(x_{t})$ is the optimal control policy. The key operations in ADP methods [32] involve setting an admissible control policy $u_{0}(x)$ and then iterating the policy evaluation step

[TABLE]

until convergence.

IV-A1 Semidefinite programming for safe initial control policy

Recall the following definition.

Definition 3**.**

The equilibrium point $x=0$ of the closed-loop system (6) is globally exponentially stable with a decay rate $\alpha$ if there exist scalars $C_{0}>0$ and $\alpha\in(0,1)$ such that $\|x_{t}\|\leq C_{0}\alpha^{(t-t_{0})}\|x_{0}\|$ for any $x_{0}\in\mathbb{R}^{n_{x}}$. $\Box$

Conditions for global exponential stability (GES) of the equilibrium state, adopted from [38], is provided next.

Lemma 2**.**

Let $V(\cdot,\cdot):[0,\infty)\times\mathbb{R}^{n_{x}}\to\mathbb{R}$ be a continuously differentiable function such that

[TABLE]

for any $t\geq t_{0}$ and $x\in\mathbb{R}^{n_{x}}$ along the trajectories of the system

[TABLE]

where $\gamma_{1}$, $\gamma_{2}$, and $\alpha$ are positive scalars, and $\varphi(\cdot)$ is a nonlinear function. Then the equilibrium state $x=0$ for the system (20) is GES with decay rate $\alpha$. $\Box$

The following design theorem provides a method to construct an initial linear stabilizing policy $u_{0}(x)=K_{0}x$ such that the origin is a GES equilibrium state of the closed-loop system (6).

Theorem 2**.**

Suppose that Assumptions 1–2 hold, and that there exist matrices $P=P^{\top}\succ 0\in\mathbb{R}^{n_{x}\times n_{x}}$, $K_{0}\in\mathbb{R}^{n_{u}\times n_{x}}$, and scalars $\alpha\in(0,1)$, $\nu>0$ such that

[TABLE]

where

[TABLE]

Then the equilibrium $x=0$ of the closed-loop system (6) is GES with decay rate $\alpha$.

Proof.

Let $V=x^{\top}Px$. Then (19a) in Lemma 2 is satisfied with $\gamma_{1}=\lambda_{\min}(P)$ and $\gamma_{2}=\lambda_{\max}(P)$. Let $\Delta V=V^{+}-V$. Note that

[TABLE]

Therefore,

[TABLE]

and

[TABLE]

Thus, pre- and post-multiplying (21) with $\begin{bmatrix}x&\phi\end{bmatrix}^{\top}$ and its transpose, respectively, we get

[TABLE]

By inequality (4) in Assumption 2 and recalling that $\phi(0)=0$, we get $(\mathfrak{L}^{*}_{\phi})^{2}q^{\top}q-\phi^{\top}\phi\geq 0.$ Since $\nu>0$, this implies $\Delta V+(1-\alpha^{2})V\leq 0,$ which is identical to (19b).

Note that we do not need to know $\phi(\cdot)$ to satisfy conditions (21). Instead, Theorem 2 provides conditions that leverage matrix multipliers similar to those described in [27].

We shall now provide LMI-based conditions for computing the initial control policy $K_{0}$, the initial domain of attraction $P$ and $\nu$ via convex programming.

Theorem 3**.**

*Fix $\alpha\in(0,1)$ and $\hat{\mathfrak{L}}_{\phi}$ obtained via (12). If there exist matrices $S=S^{\top}\succ 0$, $Y$, and a scalar $\nu>0$ such that the LMI conditions *

[TABLE]

are satisfied, then the matrices $K_{0}=YS^{-1}$, $P=S^{-1}$ and scalar $\nu$ satisfy the conditions (21) with the same $\alpha$ and $\hat{\mathfrak{L}}_{\phi}$.

Proof.

A congruence transformation of (30) with the matrix $\mathrm{blkdiag}\left(\begin{bmatrix}P&\nu^{-1}\,I&P&I\end{bmatrix}\right)$ and substituting $S$ with $P^{-1}$ and $Y$ with $K_{0}P^{-1}$ yields

[TABLE]

Taking the Schur complement with the submatrices shown by the guidelines in the above inequality, we get (28). Since $\nu>0$, taking the Schur complement again yields (29) which can be rewritten as

[TABLE]

which is exactly (21). Thus, the conditions (21) and (30) are equivalent.

Empirically, we observe that our proposed kernelized Lipschitz learner typically provides overestimates of $\mathfrak{L}^{*}_{\phi}$ (see Appendix). A benefit of overestimating $\mathfrak{L}^{*}_{\phi}$ is that admissibility of the control policy is ensured. This is demonstrated by the following result.

Theorem 4**.**

Let $(P,K_{0},\nu,\alpha)$ be a feasible solution to the conditions (21) with an overestimate of the Lipschitz constant $\hat{\mathfrak{L}}_{\phi}>\mathfrak{L}^{*}_{\phi}$. Then $(P,K_{0},\nu,\alpha)$ is also a feasible solution to the conditions (21).

Proof.

Let $\delta L=\hat{\mathfrak{L}}_{\phi}-\mathfrak{L}^{*}_{\phi}$. Since $\hat{\mathfrak{L}}_{\phi}$ is an overestimator of $\mathfrak{L}^{*}_{\phi}$, $\delta L>0$. Since $(P,K,\nu,\alpha)$ is a feasible solution to (21) with $\hat{\mathfrak{L}}_{\phi}$, it satisfies

[TABLE]

which can be written as

[TABLE]

As $\nu>0$, we infer that $\delta\mathcal{M}\preceq 0$, hence $\Gamma^{\top}\delta\mathcal{M}\Gamma\preceq 0$. Therefore, $\Psi+\Gamma^{\top}\mathcal{M}\Gamma\preceq 0$. Since the other conditions in (21) are independent of $\mathfrak{L}^{*}_{\phi}$, the other conditions are automatically satisfied. This concludes the proof.

Theorem 4 indicates that if our learned $\hat{\mathfrak{L}}_{\phi}$ is an overestimate of $\mathfrak{L}^{*}_{\phi}$, and we use $\hat{\mathfrak{L}}_{\phi}$ to obtain a safe stabilizing control policy, then this is also a safe stabilizing control policy for the true system (1). Having a feasible solution to (21) with an underestimator of $\mathfrak{L}^{*}_{\phi}$ is not sufficient to guarantee a feasible solution for the true Lipschitz constant, because $\delta\mathcal{M}$ may not be negative semi-definite in that case. Of course, extremely conservative overestimates of $\hat{\mathfrak{L}}_{\phi}$ will result in conservative control policies or result in infeasibility. In our proposed approach, we have observed that the confidence parameter $\beta$ dictates the conservativeness of the overestimate; that is $\beta\to 1$ makes the estimate $\hat{\mathfrak{L}}_{\phi}$ more conservative.

IV-A2 Safely initialized PI

We begin by proving the following critical result.

Theorem 5**.**

Let $\mathcal{U}(x,u)$ be defined as in (5). If $K_{0}$ is obtained by solving (30) for $\hat{\mathfrak{L}}_{\phi}\geq\mathfrak{L}^{*}_{\phi}$, then the initial control policy $u_{0}=K_{0}x$ is an admissible control policy on $\mathbb{R}^{n_{x}}$.

Proof.

Clearly, $u_{0}$ is continuous, and (by Theorem 2 and 3) is a stabilizing control policy for (1). It remains to show that the cost induced by $u_{0}$ is finite. Since $u_{0}$ is stabilizing and $\hat{\mathfrak{L}}_{\phi}\geq\mathfrak{L}^{*}_{\phi}$, we know that $\|x_{t}\|\to 0$ as $t\to\infty$, which implies $u_{0}\to 0$ and, by therefore, $\mathcal{U}(x_{t},u_{t})\to 0$ as $t\to\infty$. Since $\mathcal{U}(x_{t},u_{t})$ converges to a finite limit, $\mathcal{U}(x_{t},u_{t})$ is bounded for all $t\geq 0$. Therefore, any partial sum $\sum_{t=0}^{t^{\prime}}\mathcal{U}(x_{t},u_{t})$ is bounded and monotonic; that is, $\mathcal{J}$ converges to a finite limit.

Admissibility of $u_{0}$ for the specific linear quadratic regulator (LQR) cost function follows directly from Theorem 5.

Corollary 1**.**

Let

[TABLE]

for some matrices $Q=Q^{\top}\succeq 0$ and $R=R^{\top}\succ 0$. Then the initial control policy $u_{0}=K_{0}x$ obtained by solving (30) is an admissible control policy on $\mathbb{R}^{n_{x}}$.

Now that we know $u_{0}=K_{0}x$ is an admissible control policy, we are ready to proceed with the policy iteration steps (18). Typically, an analytical form of $\mathcal{J}_{k}$ is not known a priori, so we resort to a shallow neural approximator/truncated basis expansion for fitting this function, assuming $\mathcal{J}_{k}$ is smooth for every $k\in\mathbb{N}\cup\{\infty\}$. Concretely, we represent the value function and cost functions as:

[TABLE]

where $\psi_{0}(\cdot):\mathbb{R}^{n_{x}}\to\mathbb{R}^{n_{0}}$ denotes the set of differentiable basis functions (equivalently, hidden layer neuron activations) and $\omega:\mathbb{R}^{n_{0}}$ is the corresponding column vector of basis coefficients (equivalently, neural weights).

It is not always clear how to initialize the weights of the neural approximators (32). Commonly, small random numbers drawn from a uniform distribution are used [39], but there is no safety guarantee associated with random initialization. We propose initializing the weights as follows. Since our initial Lyapunov function is quadratic, we include the quadratic terms of the components of $x$ to be in the basis $\psi(x)$. Then we can express the initial Lyapunov function $x^{\top}Px$ obtained by solving (30) with appropriate weights in the $\psi(x)$, respectively, setting all other weights to be zero. With the approximator initialized as above, the policy evaluation step (18a) is replaced by

[TABLE]

This can be evaluated as $R$ and $B$ are known, and $\psi$ is differentiable and chosen by the user, so $\nabla\psi$ is computable.

Since we prove that $u_{0}$ is an admissible control policy, we can use arguments identical to [39, Theorem 3.2 and Theorem 4.1] to claim that if the optimal value function and the optimal control policy are dense in the space of functions induced by the basis function expansions (32), then the weights of the neural approximator employed in the PI steps (33) converges to the optimal weights; that is, the optimal value function $\mathcal{J}_{\infty}$ and the optimal control policy $u_{\infty}$ are achieved asymptotically. A pseudocode for implementation is provided next.

IV-B Input-constrained ADP with safety

Herein, we tackle the case when the control input is to be constrained, which is very common in practical applications. We make the following assumption on the constraints.

Assumption 5**.**

The control input $u\in\mathbb{U}$, where

[TABLE]

for $i=1,\ldots,n_{c}$, where $n_{c}$ is the number of input constraints, and $\xi_{i}\in\mathbb{R}^{n_{u}}$ for every $i$. $\Box$

Remark 8*.*

The matrix inequality (34) defines a polytopic input constraint set. Clearly, constraints of the form $|u|\leq\bar{u}$ can be written as

[TABLE]

which is of the form (34). $\Box$

Note that with any control policy $u_{0}=K_{0}x$, the constraint set described in (34) is equivalent to the set

[TABLE]

for $i=1,\ldots,n_{c}$. Before we state the main design theorem, we require the following result from [40, pp. 69].

Lemma 3**.**

The ellipsoid

[TABLE]

for $i=1,\ldots,n_{c}$. $\Box$

IV-B1 Constrained admissible initial control policy and invariant set estimation

Since the control input is constrained, we need to characterize an invariant set of the form $\mathcal{E}_{P}$ within which all control actions satisfy (34) and the following stability certificate holds.

Definition 4**.**

The equilibrium point $x=0$ of the closed-loop system (6) is locally exponentially stable with a decay rate $\alpha$ and a domain of attraction $\mathcal{E}_{P}$ if there exist scalars $C_{0}>0$ and $\alpha\in(0,1)$ such that $\|x_{t}\|\leq C_{0}\alpha^{(t-t_{0})}\|x_{0}\|$ for any $x_{0}\in\mathcal{E}_{P}$. $\Box$

A standard result for testing local exponential stability of the equilibrium point adopted from [41] is provided next.

Lemma 4**.**

Let $V:[0,\infty)\times\mathcal{E}_{P}\to\mathbb{R}$ be a continuously differentiable function such that the inequalities (19) hold for any $t\geq t_{0}$ and $x\in\mathcal{E}_{P}$ along the trajectories of the system (6) where $\gamma_{1}$, $\gamma_{2}$, and $\alpha$ are positive scalars. Then the equilibrium $x=0$ for the system (6) is locally exponentially stable with a decay rate $\alpha$ and a domain of attraction $\mathcal{E}_{P}$. $\Box$

The following design theorem provides a method to construct a stabilizing policy such that the origin is a locally exponentially stable equilibrium of the closed-loop system and constraint satisfaction is guaranteed within a prescribed ellipsoid of attraction $\mathcal{E}_{P}\subset\mathcal{X}$ without knowing the nonlinearity $\phi(\cdot)$.

Theorem 6**.**

Fix $\alpha\in(0,1)$ and $\hat{\mathfrak{L}}_{\phi}$. Suppose $\hat{\mathfrak{L}}_{\phi}\geq\mathfrak{L}^{*}_{\phi}$, and there exist matrices $S=S^{\top}\succ 0$, $Y$, and a scalar $\nu>0$ such that the LMI conditions (30) and

[TABLE]

for every $i=1,\ldots,n_{c}$. Then, the equilibrium $x=0$ of the closed-loop system (6) is locally exponentially stable with a decay rate $\alpha$ and a domain of attraction $\mathcal{E}_{P}$ defined in (36a). Furthermore, given that the initial state $x_{0}\in\mathcal{E}_{P}$, then the control actions $u_{t}$ satisfy the constraints (34) for all $t\geq 0$.

Proof.

From Theorem 2, we know that (19) holds. Taking Schur complements of (37) yields (36b) which, by Lemma 3, implies that $\mathcal{E}_{P}\subset\mathcal{X}$ and hence, the input constraints are satisfied for all $t\geq 0$ by the closed-loop system with policy $u=K_{0}x$, because $x_{0}\in\mathcal{E}_{P}$. Thus, all the conditions of Lemma 4 are satisfied, which concludes the proof.

Remark 9*.*

Note that the conditions (30) and (37) are LMIs in $S$, $Y$, and $\nu$ for a fixed $\hat{\mathfrak{L}}_{\phi}$. Therefore one can maximize the volume of $\mathcal{E}_{P}$ by solving a constrained convex program with cost function $-\log|S|$ (the log-determinant of $S$) subject to the constraints (30) and (37) while line searching for $\alpha$. This will reduce the conservativeness of the domain of attraction.

IV-B2 Safely initialized input-constrained PI

By adopting the work of [42, 43, 44, 45, 46] for input-constrained/actuator saturated ADP we choose a cost function of the form

[TABLE]

where $Q(x):\mathbb{R}^{n_{x}}\to\mathbb{R}$ is a positive definite function satisfying $Q(0)=0$ and $R\succ 0$.

We begin by demonstrating that the constrained policy is an admissible policy on its domain of attraction.

Theorem 7**.**

Let $\mathcal{U}$ be defined as in (38). Then the initial control policy $u_{0}=K_{0}x$ obtained by solving (30) and (37) is an admissible control policy on $\mathcal{E}_{P}$.

Proof.

By definition $Q(0)=0$. Also, the integrand in (38) is zero when the upper limit is zero. Therefore, $\mathcal{U}(0,0)=0$. For any $x_{0}\in\mathcal{E}_{P}$, $u_{0}$ is a stabilizing constrained control policy, therefore, $\|x_{t}\|\to 0$ and $\|u_{t}\|\to 0$ as $t\to\infty$. Hence, $\mathcal{U}\to 0$ as $t\to\infty$. The rest of the proof follows identically as in the proof of Theorem 5.

Since the control policy is constrained, we can initialize ADP safely using the neural approximator (32) as discussed in the previous subsection. The policy evaluation step is given by

[TABLE]

where $u_{1},\ u_{2},\ u_{3},\ \cdots,u_{n_{u}}$ are the individual components of the vector $u$. Subsequently, the policy improvement step is given by

[TABLE]

which satisfies the control constraints, since $\|\tanh(\cdot)\|_{\infty}\leq 1$.

Since the initial control policy is constrained and admissible, one can use [43, Theorem 2] to prove convergence of the value function and the control policy to the optimal using the constrained policy iteration steps (39) and (40).

IV-C Remarks on on-policy VI and Q-learning

The value iteration algorithm (see Algorithm 3) does not generally require an admissible control policy in order to converge optimally using data. Although this is true in off-policy implementations (that is, when the updated control policy is not used on-line), in on-policy implementations, a lack of stabilizing initial policies could result in unsafe transient behavior unless the underlying system is open-loop stable, leading to unsafe exploration during the initial data collection phase.

Q-learning is a provably convergent direct optimal adaptive control algorithm and model-free reinforcement learning technique [47, 48, 49, 50]. Q-learning can be used to find an optimal action-selection policy based on measurements of previous state and action observations controlled using a sub-optimal policy. In most of the existing work the reward/cost function is manipulated to guarantee correction of the unsafe actions in the learning phase. Our proposed method does not require a corrective modification of the reward/cost function on-line for safety. Instead, historical data and solving SDPs based on Lipschitz estimation is used to generate safe control policies that enables safe data collection during on-policy Q-learning implementation, because the states are guaranteed not to diverge with the initial policy (this divergence could happen if the initial policy was unsafe).

V Numerical Examples

V-A Nonlinear torsional pendulum

We demonstrate our proposed approach using the torsional pendulum which is modeled by discretizing the system

[TABLE]

with mass $M=0.333$ Kg, length $l=0.667$ m, acceleration due to gravity $g=0.981$ m/s2, friction factor $f_{d}=0.2$, and moment of inertia $J=0.1975$ Kg-m2. With Euler discretization and a sampling time of $\tau=0.01$ s, we get a discrete-time model of the form (1) with

[TABLE]

We assume that the nonlinearity $\phi=Mgl\sin\theta/J$ is completely unknown; clearly $\phi(\cdot)$ has a Lipschitz constant $\mathfrak{L}^{*}_{\phi}=Mgl/J=11.038$, which is also unknown to us.

In the data collection phase, we initialize the system (41) from ten different initial conditions in the space $[-\pi,\pi]\times[-2,2]$ and collect data each $0.1$ s, leading to a total dataset of $N=50$ samples. Note that the initialization procedure of [43] requires 400 data points, which is considerably more than ours, and in that procedure, the original policy in the pre-training phase is not guaranteed to be admissible. Automatic relevance determination reveals that the nonlinearity only acts through the second state, and the argument of the nonlinearity is $q=\theta$. Proceeding as in Algorithm 1, we perform cross-validation using an Epanechnikov kernel with bandwidth $h_{n}=0.05$ and choose $\beta=0.01$. This yields the overestimate $\hat{\mathfrak{L}}_{\phi}=11.511>\mathfrak{L}^{*}_{\phi}$. Using this Lipschitz estimate, we solve (30) with $\alpha=0.95$ and $\nu=1$ for an initial value function $x^{\top}Px$ and control policy estimate $K_{0}x$.

We construct a $2-11-1$ value function neural approximator with a set of polynomial basis functions

[TABLE]

where $x_{1}$, $x_{2}$ denote the first and second components of $x$, respectively. Our initial weight vector is set to

[TABLE]

where $P_{ij}$ is the $(i,j)$th element of $P$. We fix the learning rate at $\eta=10^{-4}$ and the forgetting factor $\gamma=0.95$.

Unconstrained scenario

We first test the unconstrained scenario, where the cost function is $\sum\|Qx\|_{1}+u^{\top}Ru$, with $Q=I_{2}$ and $R=0.5$, and compare four initial policies and value functions obtained via: (i) kernelized Lipschitz learning with $\beta=0.1$; (ii) kernelized Lipschitz learning with $\beta=0.001$; (iii) solving an algebraic Riccati equation and ignoring the nonlinearity; and (iv) randomly initializing with small weights from a normal distribution with small variance and zero mean (which is by far the most common initializer). The comparison study results are shown in Figure 1. We observe that all of the methods (i)–(iv) listed above work, and result in stabilizing control policies that result in the state of the torsional pendulum to converge to its equilibrium. Interestingly, both the Lipschitz constant estimates result in similar trajectories implying that the SDPs (30) are not extremely sensitive to the Lipschitz estimate. However, based on the $\mathcal{J}$ subplot which shows the variation of $\sum x^{\top}Qx+u^{\top}Ru$ with time, there is a slight improvement of performance in the $\beta=0.1$ (continuous red line) case compared to the $\beta=10^{-3}$ (dashed blue line) case since the Lipschitz estimate in the former is closer to the true Lipschitz constant. As expected, the cost incurred by the control policy ignoring the nonlinearity (dotted black line) is by far the worst, since the control actions required early on are of larger magnitude and the tracking performance is severely compromised. Randomly selecting weights also results in worse performance than our proposed method, as the cost incurred is increased due to oscillatory behaviour in the states and poor tracking in the initial time frame. Summarily, this experiment demonstrates the effectiveness of the proposed approach and its robustness to Lipschitz estimate conservatism.

In Fig. 2, we demonstrate the on-policy value iteration algorithm with safe initialization. All initial conditions converge to the origin using our proposed approach. In constrast, randomly initializing a policy and value as is typical in on-policy value iteration results in the states initially diverging (not shown in the plot) and poor performance before the rank condition is reached for determining a least-squares solution to update the neural weights.

Constrained scenario

We also test the scenario where the control actions are constrained by $|u|\leq 1$. In this case, we use the cost functional defined in (38) with $\bar{u}=1$, $Q=I_{2}$ and $R=0.5$. We begin by solving (30) and (37) with $\nu=1$ and $\alpha=0.95$ to get $P$ and $K_{0}$, as in the previous subsection. We also select the same basis functions (42). We randomly initialize (using 20 random initial conditions) the system (41) from within the domain of attraction of the initial control policy, that is, from within the set $\{x^{\top}Px\leq 1\}$. We know from Theorem 6 that this ensures that the initial control policy will satisfy input constraints. Consequently, because the policy improvement step is also guaranteeed to satisfy input constraints and the initial policy is stabilizing, there are no constraint violations, and the initialization is deemed safe. The performance of the proposed algorithm is provided in Figure 3. The convergence of $\|x_{t}\|$ to zero and the satisfaction of input bounds are illustrated. Finally, we demonstrate the convergence of the neural weights $\omega_{t}$, noting that learning did occur, that is, the weights were not static throughout the simulation (which would indicate that the initial policy was optimal).

V-B Large randomized linear system

In order to study the scalability of the proposed approach on higher state-space dimensions, we randomly generate a 20-state, 10-input linear system of the form $x_{t+1}=Ax_{t}+Bu_{t}$, where $B$ is known, and $A$ is unknown. We randomly choose $A_{0}\neq A$ to be a known matrix such that $(A_{0},B)$ is controllable. Clearly, the unmodeled component is $\phi(x)=(A-A_{0})x$. We assume we do not know $G$ and $C_{q}$, so both are set to identity matrices of appropriate dimensions. The initial dataset is generated using a small random perturbation signal on the unknown system and 500 data points are stored, from which we compute a Lipschitz constant for $\phi$ as $\hat{\mathfrak{L}}_{\phi}=0.43$ (the true $\mathfrak{L}^{*}_{\phi}=0.34$) and a safe initial policy with $\beta=0.001$; this is illustrated in Figure 4’s top left subplot; the shaded blue area is the subgraph of $\hat{L}_{n}$ within the $\beta$-confident support. The unknown system is simulated from 100 randomly selected initial conditions in $\mathbb{R}^{20}$ following a normal distribution with variance 2 and zero-mean. The cost function has the form (31) with $Q=5I_{20}$, $R=2I_{10}$, and $\gamma=0.95$. In all cases, the weights of the $20-210-10$ neural approximator with learning rate $\eta=0.1$ converges within a few seconds, and the system stabilizes, in spite of $A$ being unstable, as shown in Figure 4.

VI Conclusions

This work provides a methodology for constructing admissible initial control policies for ADP methods using Lipschitz learning by using kernel density estimation and semi-definite programming. Such admissible controllers enable safe initialization, that is, with constraint satisfaction using only historical data, which is necessary not only in policy iteration methods, but also in value iteration and Q-learning for safely obtaining initial data on-line for on-policy learning when the underlying system is not open-loop stable. Simulations on a discretized torsional pendulum model and a high-dimensional linear system are provided to show the efficiency of our approach. Future research efforts will focus on more general costs and uncertain nonlinear safety constraints while ensuring feasibility with a high probability in terms of regret.

Acknowledgments

We would like to thank Drs. Mouhacine Benosman and Piyush Grover at Mitsubishi Electric Research Laboratories, Cambridge, MA, USA, for their time and helpful insights. Kyriakos Vamvoudakis was supported in part by NSF under grant Nos. CPS-1851588 and S&AS-1849198.