An LMI Framework for Contraction-based Nonlinear Control Design by Derivatives of Gaussian Process Regression

Consider the following discrete-time nonlinear system:

where $x_{k}\in{\mathbb{R}}^{n}$ and $u_{k}\in{\mathbb{R}}$ denote the state and input, respectively; $f:{\mathbb{R}}^{n}\to{\mathbb{R}}^{n}$ is of class $C^{1}$, i.e., continuously differentiable, and $b\in{\mathbb{R}}^{n}$. The $i$th components of $f$ and $b$ are denoted by $f_{i}$ and $b_{i}$, respectively. For the sake of notational simplicity, we consider single-input systems. However, the results can readily be generalized to multiple-input systems as explained below. We also later discuss the case where $b$ is a function of $x$.

In contraction theory [28, 37, 22], we study incremental stability as a property of the pair of trajectories, stated below.

Applying [37, Theorem 15] to control design yields the following IES condition.

By the Schur complement, the set of (2) and (3) is equivalent to $P\succ 0$ and

This is nothing but the definition of uniform contraction with $\Theta=P^{-1/2}$ [37, Definition 6] for the closed-loop system $x_{k+1}=f(x_{k})+bp(x_{k})$. Therefore, [37, Theorem 15] concludes IES of the closed-loop system. ∎

In [37, Theorem 15], it has been shown that a closed IES system admits a state-dependent $P$ satisfying a similar inequality as (2). Moreover, such a $P$ is uniformly lower and upper bounded by constant matrices. It is not yet clear when $P$ becomes constant. If one restricts the class of controllers into linear for a constant $P$, control design can be reduced to linear matrix inequalities (LMIs); see, e.g., [11, 19]. Namely, control design can be sometimes formulated as a practically solvable problem even for nonlinear systems. Indeed, (2) is an LMI with respect to $\varepsilon$, $P$ and $p_{\partial}(x)$ at each $x\in{\mathbb{R}}^{n}$. However, as in Proposition 2.2, designing nonlinear controllers is not fully formulated in the LMI framework, due to the so-called integrability constraint (3), i.e., $p_{\partial}(x)P^{-1}$ needs to be the partial derivative of some function $p(x)$ providing a feedback control law $u=p(x)$. The main objective of this paper is to develop an LMI framework for stabilizing nonlinear control design by tackling the integrability constraint, stated below.

An important feature of contraction theory is to study the convergence between any pair of trajectories. By virtue of this, one can handle observer design [1, 28], tracking control [15, 33], and control design for imposing synchronizations [1] and rich behavior such as limit cycles [11] in the same framework as stabilizing control design. These references mainly focus on continuous-time systems, but similar results can be delivered to discrete-time systems, since incremental stability conditions in contraction analysis have been derived also for discrete-time systems [28, 37]. Moreover, the proposed method in this paper can be generalized to continuous-time systems as will be explained in Section 3.3. Therefore, solving Problem 2.3 can result in LMI frameworks for various design problems.

Integrability constraints sometimes appear in nonlinear adaptive control or observer design. Since directly addressing integrability constraints are difficult, there are techniques for avoiding them by adding the dynamic order of the identifier. However, as pointed out by [27], adding additional dynamics can degenerate control performances, and it has not been validated that such an approach works for contraction-based control design. Therefore, it is worth solving Problem 2.3 directly. In particular, we provide an LMI framework for control design, which is easy-to-implement. The proposed method may be tailored for nonlinear adaptive control or observer design although this is beyond the scope of this paper.

To solve Problem 2.3, we employ Gaussian process regression (GPR) [32, 4]. We, in this subsection, briefly summarize basics of GPR and, in the next section, demonstrate that GPR is a suitable tool for handling problems involving partial derivatives, e.g., integrability conditions.

Let $\{x^{(i)}\}_{i=1}^{N}$, $x^{(i)}\in{\mathbb{R}}^{n}$ be input data, and let $\{y^{(i)}\}_{i=1}^{N}$, $y^{(i)}\in{\mathbb{R}}$ be the corresponding output data given by

where $\omega^{(i)}\sim{\mathcal{N}}(0,\sigma^{2})$ is i.i.d. GPR is a technique to learn an unknown function $p:{\mathbb{R}}^{n}\to{\mathbb{R}}$ from input-output data $\{x^{(i)},y^{(i)}\}_{i=1}^{N}$ by assuming $p$ as a Gaussian process (GP).

A stochastic process $p$ is said to be GP if any finite set $\{p(x^{(i)})\}_{i=1}^{N}$ has a joint Gaussian distribution [32, Definition 2.1]. A GP is completely specified by its mean function $m_{p}:{\mathbb{R}}^{n}\to{\mathbb{R}}$ and covariance function $k_{p}:{\mathbb{R}}^{n}\times{\mathbb{R}}^{n}\to{\mathbb{R}}$. They are defined by

and we represent the GP by $p\sim\mathcal{GP}(m_{p}(x),k_{p}(x,x^{\prime}))$.

The essence of GPR is to estimate $m_{p}(x)$ and $k_{p}(x,x^{\prime})$ as the posterior mean and covariance given $\{x^{(i)},y^{(i)}\}_{i=1}^{N}$ by the Bayes estimation [32, 4]. Typically, the prior mean $m_{0}(x)$ is selected as zero. The prior covariance $k_{0}(x,x^{\prime})$ needs to be a positive definite kernel [32, Section 6], and we also require smoothness. Standard kernels are linear, polynomial, or squared exponential (SE) [32, 4], which are all smooth and positive definite. For instance, an SE kernel is

where $\beta>0$ and $\Sigma\succ 0$ in addition to $\sigma>0$ in (5) are free parameters, called hyper parameters. The hyper parameters can be selected to maximize the marginal likelihood from observed data; see e.g., [32, Section 5.4].

To use GPR for learning $p(x)$, we need its data as in (5). However, Proposition 2.2 contains no information of $p(x)$ but of $\partial p(x)$ via the integrability constraint (3). To handle this situation, we employ derivatives of GPR as elaborated in the next section.

Taking the partial derivatives of (6) and (7), we have the joint distribution of $p(x)$ and $\partial p(x)$ as in

see, e.g., [30, Equation (2)]. The goal of this subsection is to find $m_{p}(x)$ based on the Bayes estimation such that $p(x)=m_{p}(x)$ is a solution to Problem 2.3. For the sake of notational simplicity, we select the prior mean of $p(x)$ as zero.

A standard procedure of the Bayes estimation is that we first select a class $C^{2}$ positive definite kernel $k_{0}(x,x^{\prime})$ as a prior covariance. Then, we compute the posterior mean $m_{p}(x)$ given data of $\partial p(x)$. Looking at this from a different angle, $m_{p}(x)$ can be viewed as a function of data of $\partial p(x)$. Based on this perspective, we consider generating suitable data such that $m_{p}(x)$ becomes a solution to Problem 2.3 as a non-standard use of GPR.

To this end, let $\{x^{(i)},\hat{p}^{(i)}\}_{i=1}^{N}$ denote a data set to be generated, where

The role of i.i.d. $\omega_{p}^{(i)}\sim{\mathcal{N}}(0,\sigma_{p}^{2}I_{n})$ is explained later; one can simply choose it as zero. Define the vector consisting of $\{\hat{p}^{(i)}\}_{i=1}^{N}$ by

Also, we introduce the following notations:

Then, the joint distribution of the prior distribution of $p(x)$, denoted by $p^{\rm pri}(x)$ and $\partial p(X_{N})$ is

By the Bayes estimation, we can compute the posterior mean $m_{p}(x)$ given $Y_{p}$ as

where $h_{p}^{(i)}$ denotes the $i$th block-component with size $n$ of $h_{p}$ defined by

The partial derivative of $m_{p}(x|Y_{p})$ is easy-to-compute:

Especially if $\sigma_{p}=0$ and $K_{0}\succ 0$, it follows from (14) that

Note that $m_{p}(x|Y_{p})$ and $\partial m_{p}(x|Y_{p})$ are nonlinear functions of $x$ and linear functions of $Y_{p}$. Therefore, control design reduces to generating suitable $Y_{p}$, i.e., $\{x^{(i)},\hat{p}^{(i)}\}_{i=1}^{N}$, which is proceeded based on Proposition 2.2.

Now, we are ready to develop an LMI framework for contraction-based nonlinear control design. Substituting $p_{\partial}=\partial m_{p}(x|Y_{p})P$ into (2) does not give an LMI because of the coupling between $Y_{p}$ and $P$. This issue can be addressed by using a standard technique of an LMI. According to [34, Theorem 2.3.11], (2), i.e., the set of $P\succ 0$ and (2.1) implies, for all $x\in{\mathbb{R}}^{n}$,

This is an LMI with respect to $P$ and $\varepsilon_{p}>0$ at each $x\in{\mathbb{R}}^{n}$. For a solution $P$ to (18) and $p_{\partial}=\partial m_{p}(x|Y_{p})P$, one only has to solve (2) with respect to $Y_{p}$. The proposed procedure for solving Problem 2.3 is summarized in Algorithm 1.

If Algorithm 1 has a set of solutions, the obtained $p(x)$ is a solution to Problem 2.3 at each data point, stated below.

Using a set of solutions $p(x)$ and $P$, define $p_{\partial}(x)=\partial p(x)P$. Then, (3) holds for all $x\in{\mathbb{R}}^{n}$. It suffices to confirm that (2) holds at $x=x^{(i)}$, $i=1,\dots,N$. Since $p(x)=m_{p}(x|Y_{p})$, we have $\partial p(x^{(i)})=\partial m_{p}(x^{(i)}|Y_{p})$, $i=1,\dots,N$. Therefore, (23) implies (2) for $x=x^{(i)}$, $i=1,\dots,N$. ∎

Note that in Algorithm 1, the number $N$ of training data and data point set $\{x^{(i)}\}_{i=1}^{N}$ can be chosen arbitrarily. As $N$ increases, the number of $Y_{p}\in{\mathbb{R}}^{nN}$ in (23) increase. However, the problem is still convex. As explained in the next subsection, for sufficiently large $N$, one can show that $p(x)=m_{p}(x|Y_{p})$ is a solution to Problem 2.3 other than $\{x^{(i)}\}_{i=1}^{N}$ when $x^{(i)}$ is distributed evenly in the state space. Moreover, such $N$ can be found.

When $\sigma_{p}=0$ and $K_{0}\succ 0$, (15) helps to simplify Algorithm 1. In fact, we only have to solve a finite family of LMIs once, stated below without the proof.

At the end of this subsection, we argue the role of $\omega_{p}^{(i)}\sim{\mathcal{N}}(0,\sigma_{p}^{2}I_{n})$ in (10) by interpreting the procedure of Algorithm 1 as follows. We generate suitable $Y_{p}$ based on the LMIs and then construct $p(x)=m_{p}(x|Y_{p})$ by functional fitting. In functional fitting, overfitting is a common issue, since a constructed function becomes unnecessarily complex. The variance $\sigma_{p}^{2}$ specifies how much a constructed function needs to fit to data, and thus adding the noise $\omega_{p}^{(i)}$ helps to avoid overfitting. However, if $\sigma_{p}^{2}$ is large, a constructed function can fit to $\omega_{p}^{(i)}$ instead of $\partial p(x^{(i)})$. This is well known as the bias-variance tradeoff in functional fitting [16, Section 5.4.4]. It is worth emphasizing that Theorem 3.1 holds for arbitrary $\sigma_{p}>0$.

If $K_{0}$ in (11) is non-singular, $h_{p}$ in (13) is nothing but the minimizer of the following optimization problem:

The first term evaluates the fitting error at each $x^{(i)}$ and the second one is for regularization. Especially when $\sigma_{p}=0$, the optimal value becomes zero for the obtained $h_{p}$, i.e., a complete fitting $\partial p(x^{(i)})=\hat{p}^{(i)}$, $i=1,\dots,N$ is achieved as stated by Corollary 3.2. A reproducing kernel Hilbert space is a formal tool to study a functional fitting problem under regularization as an optimization problem. For the prior variance as a kernel, the pair $(p,\partial p)$ can be understood as a minimizer according to the representer theorem [31, Theorem 2].

In Theorem 3.1, it is not clear whether $p_{\partial}(x)=\partial p(x)P$ satisfies (2) other than $\{x^{(i)}\}_{i=1}^{N}$ in contrast to the integrability condition (3). Applying standard arguments of the polytope approach [5], we improve the LMIs in Algorithm 1 such that its solution satisfies (2) other than $\{x^{(i)}\}_{i=1}^{N}$ when $N$ is sufficiently large.

For a set of matrices ${\mathbb{A}}:=\{A_{1},\dots,A_{L}\}$, let ${\rm ConvexHull}({\mathbb{A}})$ denote its convex hull, i.e.,

Now, we choose $L=N$ and

From the standard discussion of the polytope approach [5], (2) and (3) hold for all $x\in{\mathbb{R}}^{n}$ belonging to

This set $D$ can be made larger by increasing the number $N$ of data used for control design.

Increasing $N$ is not the only approach to enlarge a set of $x$ in which (2) holds. Another approach is to improve the LMIs in Algorithm 1. First, we decide ${\mathbb{A}}^{(i)}:=\{A^{(i)}_{1},\dots,A^{(i)}_{L^{(i)}}\}\subset{\mathbb{R}}^{n\times n}$ which is arbitrary as long as

Instead of (21), we consider the following finite family of LMIs with respect to $P\in{\mathbb{R}}^{n\times n}$ and $\varepsilon_{p}>0$:

Using a solution $P$, we consider the following finite family of LMIs with respect to $Y_{p}\in{\mathbb{R}}^{nN}$ and $\varepsilon>0$ instead of (23):

The newly obtained $p(x):=m_{p}(x|Y_{p})$ satisfies (2) for all $x\in{\mathbb{R}}^{n}$ belonging to

Now, we are ready to state a control design procedure such that the IES conditions (2) and (3) hold on a given bounded set.

First, we show (31). Since a continuous function is bounded on a bounded set, $\{\partial f(x)\in{\mathbb{R}}^{n\times n}:x\in D^{(i)}\}$ is bounded for each $i=1,\dots,N$. Each bounded set admits its convex hull. Therefore, there exists ${\mathbb{A}}^{(i)}$ satisfying (31).

Next, we consider the latter statement. We choose $p_{\partial}(x)=\partial p(x)P=\partial m_{p}(x|Y_{p})P$. Then, (3) holds on ${\mathbb{R}}^{n}$. It remains to show that (2) holds on $D$. Items 1) and 2) and (31) imply that for each $i=1,\dots,N$,

The strict inequality holds for any positive $\bar{\varepsilon}<\varepsilon$. Since $p_{\partial}(x)=\partial m_{p}(x|Y_{p})P$ is continuous (with respect to $x$), there exists a sufficiently small $\bar{D}_{\bar{\varepsilon}}^{(i)}\subset{\mathbb{R}}^{n}$ centered at $x^{(i)}$ such that

From item 3), the continuity of $\partial f$ and $p_{\partial}$, and the boundedness of $D^{(i)}$, there exists a sufficiently large $r\in{\mathbb{Z}}_{>0}$ such that $D^{(i)}\subset\bar{D}_{\bar{\varepsilon}}^{(i)}$ for all $i=1,\dots,N$. Consequently, (2) hold on $D$, a union of $D^{(i)}$, $i=1,\dots,N$. ∎

In the above theorem, $D$ and $D^{(i)}$ are not necessarily to be hypercubes or closed. Essential requirements are that $D$ is covered by $\{D^{(i)}\}_{i=1}^{N}$, and each $D^{(i)}$ shrinks as $N$ increases. In the proof, we show that a switching controller $u=\partial p(x^{(i)})x$, $x\in D^{(i)}$ also satisfies (2) and (3). However, a continuous controller $u=p(x)$ is more easy-to-implement.

In this subsection, we discuss how to generalize our results to the cases where the input vector field $b$ is a function of $x$. Also, we mention the continuous-time case.

First, we consider the system with non-constant $b$:

where $b:{\mathbb{R}}^{n}\to{\mathbb{R}}^{n}$ is of class $C^{1}$. A modification of Proposition 2.2 implies that a controller $u=p(x)$ achieves IES if there exist $\varepsilon>0$, $P\in{\mathbb{R}}^{n\times n}$, and $p:{\mathbb{R}}^{n}\to{\mathbb{R}}$ of class $C^{1}$ such that for all $x\in{\mathbb{R}}^{n}$,

The difference from (2) is the additional term $p(x)\partial b(x)$.

For finding $P$ first, one can utilize a modification of (18):