cuhk1]Department of Mechanical and Automation Engineering, The Chinese University of Hong Kong, HKSAR, China cuhk2]Department of Mechanical and Automation Engineering, The Chinese University of Hong Kong, HKSAR, China

On Estimating the Probabilistic Region of Attraction for Partially Unknown Nonlinear Systems: An Sum-of-Squares Approach

Hejun Huang\arefcuhk1 Dongkun Han\arefcuhk2 [ hjhuang@mae.cuhk.edu.hk [ dkhan@mae.cuhk.edu.hk

Abstract

Estimating the region of attraction for partially unknown nonlinear systems is a challenging issue. In this paper, we propose a tractable method to generate an estimated region of attraction with probability bounds, by searching an optimal polynomial barrier function. Techniques of Chebyshev interpolants, Gaussian processes and sum-of-squares programmings are used in this paper. To approximate the unknown non-polynomial dynamics, a polynomial mean function of Gaussian processes model is computed to represent the exact dynamics based on the Chebyshev interpolants. Furthermore, probabilistic conditions are given such that all the estimates are located in certain probability bounds. Numerical examples are provided to demonstrate the effectiveness of the proposed method.

keywords:

Region of attraction, Sum-of-squares programming, Chebyshev interpolants, Gaussian processes.

1 Introduction

Tracking the performance of uncertain nonlinear systems is an essential problem of significant research interests. In engineering practices, the concepts of safety and stability are central to these uncertain nonlinear systems in most scenarios, e.g., flight dynamics, bipedal robotics and power systems [1, 2, 3, 4]. Consequently, scientists across multiple disciplines have summarized this type of problem into the analysis of region of attraction (ROA, also called domain of attraction). Estimating the ROA can directly obtain the safety or the stability margin in many practical implementations [5].

Fruitful results have been obtained to estimate the ROA of fixed nonlinear systems. Exploiting the sublevel set of Lyapunov functions is one of the useful methods [6]. Barrier functions are another powerful tool to guarantee safety such that states avoid entering a specified unsafe region [7], [8]. To meet the safety and the stability conditions simultaneously, the quadratic programs have been investigated to find qualified Lyapunov barrier functions [9]. Meanwhile, sum-of-squares programmings (SOSPs) are proposed to find more permissive results in polynomial systems [10, 11, 12, 13]. However, in many engineering practices, there exist model inaccuracies and unknown disturbances that might influence the system dynamics or even result in operation failures in the worst case.

To handle with the above issue, effective methods have been proposed regarding partially unknown systems. First, Lyapunov-based methods are developed and extended to uncertain systems, where a Lyapunov certified ROA (LCROA) estimation is conducted for uncertain, polynomial systems [14], [15]. Furthermore, learning-based methods have been studied in the literature. Among these, the use of Gaussian processes (GP) is shown to be a promising approach to quantify the uncertainty in the stochastic process [16, 17]. The idea of unifying SOSP and GP naturally arises to compute the LCROA [18, 19, 20].

Inspired by the work in [21] that firstly uses GP to compute the barrier certified ROA (BCROA), based on our previous work in [13, 22], this paper combines SOSPs and GP to estimate the optimal BCROA of partially unknown nonlinear systems. Different from the learned polynomial systems in [20], we use Chebyshev interpolants and polynomial mean functions of GP models to find the optimal BCROA rather than LCROA with relaxed assumptions. The polynomial mean function is proven to be more flexible to match other nonlinear kernels instead of polynomial kernel. To the best of our knowledge, this paper is the first to compute barrier functions via SOSPs in partially unknown nonlinear systems.

The main contributions of this paper are threefold. First, a learned polynomial system is built with probability bounds. Second, a positive but safe sample policy is developed to prepare appropriate prior information for higher accuracy of the GP model. Third, a tractable method based on SOSP is proposed to compute the probabilistic optimal BCROA.

2 Preliminary

Consider an autonomous system as follows,

\dot{x}=\underbrace{f(x)+g(x)}_{Known}+\underbrace{d(x)}_{\textit{Unknown}},

(1)

where $x\in\mathcal{X}\subseteq\mathbb{R}^{n}$ denotes the state, $f,g:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ denote the polynomial and the non-polynomial term respectively, and $d:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ denotes the unknown term. All the terms in (1) are Lipschitz continuous.

Without losing generality, systems with a single equilibrium are considered, and this equilibrium could be transformed to the origin via variables substitution [23].

Assumption 1.

The origin ( $x=0$ ) is a single stable equilibrium of (1), that is $f(0)=g(0)=d(0)=0$ . $\hfill\square$

To model the system dynamics in a Bayesian framework as developed in [24], we define a prior distribution of noise $[\epsilon_{1},\epsilon_{2},\dots,\epsilon_{k}]^{\mathrm{T}}$ over $k$ measurements $[\dot{x_{1}},\dot{x_{2}},\dots,\dot{x_{k}}]^{\mathrm{T}}$ .

Assumption 2.

The noise $[\epsilon_{1},\epsilon_{2},\dots,\epsilon_{k}]^{\mathrm{T}}$ over the system measurements $[\dot{x}_{1},\dot{x}_{2},\dots,\dot{x}_{k}]^{\mathrm{T}}$ of (1) is uniformly bound by $\sigma_{n}$ , i.e., $\{\epsilon_{i}\sim\mathcal{N}(0,\sigma_{n}^{2})\}_{i=1}^{k}$ . $\hfill\square$

The known dynamics $f(x)+g(x)$ can be computed directly. Thus, $d(x)$ can be obtained by subtracting $f(x)+g(x)$ from $\dot{x}$ . In this work, we restrict our attention to the system (1) with bounded $d(x)$ :

Assumption 3.

The unknown term $d(x)$ in (1) exists a bounded norm in the reproducing kernel Hilbert space (RKHS), i.e., $\|d(x)\|\leq c_{g}$ , where $c_{g}$ is a constant. $\hfill\square$

As [21] introduced, the RKHS is a Hilbert space of square integrable functions that contains functions of the form $l(x)=\sum_{i}\alpha_{i}k(x,x_{i})$ , where $\alpha_{i}\in\mathbb{R}$ denotes coefficient and $k:\mathcal{X}\times\mathcal{X}\rightarrow\mathbb{R}^{+}_{0}$ denotes a symmetric positive definite kernel function of states $x,x_{i}$ . For more details about the RKHS norm, we kindly refer interested readers to [25].

2.1 Barrier Functions

To introduce the barrier function, let us first consider a simple autonomous system as follows,

\displaystyle\dot{x}=f(x),

(2)

where $f(x)$ is locally Lipschitz continuous with a single stable equilibrium point at the origin. As we mentioned, the safety and stability of (2) can be guaranteed by ROA.

State trajectories starting inside the barrier function certified ROA (BCROA) $\mathcal{L}=\{x\in\mathcal{X}|h(x)\geq 0\}$ will never enter into the unsafe region $\mathcal{L}_{U}=\mathcal{X}\backslash\mathcal{L}$ as defined,

$\displaystyle\forall x\in\mathcal{L}\;\;\;\;$	$\displaystyle h(x)\geq 0,$	(3)
$\displaystyle\forall x\in\mathcal{L}$	$\displaystyle\frac{\partial h}{\partial x}f(x)\geq 0,$
$\displaystyle\forall x\in\mathcal{L}_{U}$	$\displaystyle h(x)<0.$

Meanwhile, the Lyapunov function $V(x)$ certified ROA (LCROA) $\mathcal{L}_{V}=\{x\in\mathcal{X}\,|V(x)\leq c\}$ denotes a sublevel set of $V(x)$ . The state trajectories inside $\mathcal{L}_{V}$ will always converge to the origin. The following lemma demonstrates a relationship between $\mathcal{L}$ and $\mathcal{L}_{V}$ .

Lemma 1.

(Lemma 3.2 of [13]) Given an autonomous dynamical system (2) that is asymptotically stable at the origin, the estimate of ROA by barrier function is no smaller than the estimate by Lyapunov function, i.e., $\mathcal{L}_{v}\subseteq\mathcal{L}$ . $\hfill\square$

The details about how to obtain a Lyapunov maximum sublevel $c^{*}$ and an optimal barrier function $h^{*}(x)$ are introduced [13]. In Example 1, we will illustrate the comparison of LCROA and BCROA.

Example $1$ : Consider an autonomous system

\displaystyle\begin{bmatrix}\dot{x}_{1}\\ \dot{x}_{2}\end{bmatrix}=\begin{bmatrix}-x_{1}^{3}-x_{1}x_{2}^{2}\\ -x_{2}-x_{1}^{2}x_{2}\end{bmatrix}.

(4)

Let $\mathcal{A}_{1}$ denote the LCROA and $\mathcal{A}_{2}$ denote the BCROA. In (4), $\mathcal{A}_{1}=\{x\,|\,V_{1}(x)\leq 1.063\}$ is established by a $4^{th}$ degree polynomial $V_{1}(x)=x_{1}^{4}+x_{2}^{4}+x_{1}^{2}x_{2}^{2}+x_{1}^{2}+x_{2}^{2}+x_{1}x_{2}$ , and $\mathcal{A}_{2}=\{x\,|\,h_{1}^{*}(x)\geq 0\}$ is obtained by another $4^{th}$ degree polynomial barrier function $h_{1}^{*}(x)$ based on $\mathcal{A}_{1}$ .

Refer to caption — Figure 1: Comparison of $\mathcal{A}_{1}$ and $\mathcal{A}_{2}$ . The red region enclosed by a dashed blue ellipse depicts $\mathcal{A}_{1}$ , while the light blue region enclosed by the solid black line depicts $\mathcal{A}_{2}$ .

In Figure 1, $\mathcal{A}_{2}$ is significantly larger than $\mathcal{A}_{1}$ , which means the former allows more states to be considered in (4). In this paper, we will leverage the superiority of the barrier function to obtain a permissively safe and stable guarantee.

2.2 Gaussian Processes

Gaussian processes (GP) provide a non-parametric regression method to capture the unknown dynamics in Bayesian framework. In the GP model, every variable is associated with a normal distribution, where the prior and the posterior GP model of these variables obey a joint Gaussian distribution [24]. A GP model is characterized by the mean function $m(x)$ and the covariance function $k(x,x^{\prime})$ , where $k(x,x^{\prime})$ computes the similarity between any two states $x,\,x^{\prime}\in\mathcal{X}$ . In most cases, the covariance function $k(x,x^{\prime})$ is also called a kernel. The measurements of $d(x)$ in (1) can be used to construct a GP model with polynomial mean functions as follows.

Proposition 1.

Suppose there exist $k$ measurements of $d(x)$ in (1) that satisfy Assumption 2, 3. Then, the following GP model of $d(x)$ can be established with polynomial mean function $m(x_{*})$ and covariance function $k(x,x_{*})$ as,

	$\displaystyle m(x_{*})$	$\displaystyle=\varphi(x_{*})^{\mathrm{T}}w,$		(5)
	$\displaystyle k(x,x_{*})$	$\displaystyle=k(x_{},x_{})-k_{}^{\mathrm{T}}(K+\sigma_{n}^{2}I)^{-1}k_{}.$		(5)

where $x_{*}$ is a query state, $\varphi(x_{*})$ is a monomial vector, $w$ is a coefficient vector, $[K]_{(i,j)}=k(x_{i},x_{j})$ is a kernel Gramian matrix and $k_{*}=[k(x_{1},x_{*}),k(x_{2},x_{*}),\dots,k(x_{w},x_{*})]^{\mathrm{T}}$ .

Proof.

See Appendix. ∎

Remark 1.

Proposition 1 supports us to define the degree of polynomial $m(x_{*})$ in (5). Besides, it establishes a feasible link between polynomial mean functions and other flexible kernels, including polynomial kernel [20]. $\hfill\square$

3 Learning the Partially Unknown System

To illustrate our work about learning the partially unknown dynamics (1) in polynomial form, we propose the approximated probabilistic model in (12), which combines Chebyshev interpolants and a GP model as shown in the first two subsections. We also introduce a covariance oriented safe sample policy based on GP to learn the dynamics consistently in the third subsection.

3.1 Chebyshev Interpolants

Chebyshev interpolants provide a useful way to approximate a class of nonlinear functions with a bounded remainder [26]. The term $g(x)$ in (1) can be approximated by the Chebyshev interpolants $P_{k}$ of degree $k$ in $[-1,1]$ as,

\displaystyle g(x)\approx P_{k}(x)

\displaystyle=\sum_{i=0}^{k}\alpha_{i}\tau_{i}(x)=\left\langle A,T\right\rangle,\quad i\in[0,k],

(6)

where $i$ denotes the subscript of the coefficients vector $A=[\alpha_{0},\alpha_{1},\dots,\alpha_{k}]$ and the Chebyshev polynomials vector $T=[\tau_{0},\tau_{1},\dots,\tau_{k}]$ that satisfies,

	$\displaystyle\alpha_{i}(x)$	$\displaystyle=\left\{\begin{aligned} &\frac{1}{\pi}\int_{-1}^{1}\frac{{g(x)\tau_{i}(x)}}{{\sqrt{1-x^{2}}}}dx,\quad i=0,\\ &\frac{2}{\pi}\int_{-1}^{1}\frac{{g(x)\tau_{i}(x)}}{{\sqrt{1-x^{2}}}}dx,\quad i\in[1,k],\\ \end{aligned}\right.$		(7)
	$\displaystyle\tau_{i}(x)$	$\displaystyle=\cos(i\operatorname{arccos}(x)),\qquad\quad i\in[0,k].$		(7)

Note that, Chebyshev interpolants are applicable to any arbitrary interval $I=[a,b]$ by the following transformation,

I(x)=\frac{2x-(b+a)}{b-a}.

(8)

We define the remainder $\xi=g(x)-P_{k}(x)$ based on (6). The following inequality from [26] declares that the remainder $\xi$ of Chebyshev interpolants is bounded over the domain.

Lemma 2.

(Theorem 8.2 of [26]) Let an analytic function $g(x)$ in $[-1,1]$ be analytically containable to the open Bernstein ellipse $E$ , where it satisfies $|g(x)|\leq m$ for some $m$ . Then, its Chebyshev interpolants $P_{k}(x)$ satisfies $\begin{aligned} \|g(x)-P_{k}(x)\|\leq\frac{4m\rho^{-k}}{\rho-1},\quad k\geq 0.\end{aligned}$ $\square$

The Bernstein ellipse $E$ has foci $\pm 1$ and major radius $1+\rho$ for all $\rho\geq 0$ . For more details about the Bernstein ellipses, we kindly refer to the book by Trefethe [26].

Based on Lemma 2, the system (1) can be expressed as,

\displaystyle\dot{x}

\displaystyle=\;f(x)+P_{k}(x)+d_{\xi}(x),

(9)

where the unknown term $d_{\xi}(x)$ satisfies Assumption 2, 3,

\displaystyle d_{\xi}(x)=d(x)+\xi.

(10)

The approximation accuracy of Chebyshev interpolants is linearly dependent on the degree. We will illustrate a comparison in the following example.

Example $2$ : Consider a nonlinear function $y=\sqrt{|x-3|}$ approximated by the Chebyshev interpolants of degree $4$ and $80$ in $[0,6]$ , the higher degree approximation is more accurate with less bounded remainder $\xi$ in Figure 2.

3.2 Probability Bounds on the GP model

The following result provides the probability bounds in the distribution of a GP model.

Lemma 3.

(Lemma 2 of [17]) Suppose there exist $k$ measurements $\{(x_{1},y_{1}),(x_{2},y_{2}),\dots,(x_{k},y_{k})\}$ of a bounded function $h(x)\in\mathcal{H}$ , where $x$ denotes the state, $y$ denotes the measurement of $h(x)$ that corrupted with noise $\{\epsilon_{1},\epsilon_{2},\cdots\epsilon_{k}\}\sim\mathcal{N}(0,\sigma_{n}^{2})$ , and $\mathcal{H}$ denotes the RKHS. Let $\delta\in(0,1)$ . For $h(x)$ and its inferred GP model $(m_{h}(x),\sigma^{2}_{h}(x))$ holds,

\mathbb{P}\{\|h(x)-m_{h}(x)\|\leq\|\beta\|\|\sigma_{h}(x)\|,\forall x\in\mathcal{X}\}\geq(1-\delta)^{k},

(11)

where $[\beta]_{k}=\sqrt{2\lVert h_{k}\rVert^{2}+300\gamma_{k}\log^{3}(k/\delta)}$ is a discounting factor, $\|\cdot\|$ is a norm of $\mathcal{H}$ , and $\gamma_{k}$ is a factor denoting the maximum mutual information over the measurements. $\hfill\square$

Remark 2.

The value of $\gamma_{k}$ is related to the type of kernels, e.g., RBF kernel, Matérn kernel and linear kernel from [27]. Thus, its sublinear dependent term $\beta_{k}$ can be regarded as a constant [18]. In [17], a more correlated statement yields that with more appropriate prior information, the value of $\sigma_{h}(x)$ will decrease such that the dynamics $h(x)$ can be represented by $m_{h}(x)$ with fewer differences. $\hfill\square$

The unknown term $d_{\xi}(x)$ in (10) can be learned by using a GP model. Combining previous results of the system (9), a probabilistic statement toward the exact dynamics,

\displaystyle\dot{x}=f(x)+P_{k}(x)+m_{d_{\xi}}(x),

(12)

holds with probability greater or equal to $(1-\delta)^{k},\delta\in(0,1)$ . Therefore, we can obtain a high probability statement (12) of (1) based on Lemma 3, which will be needed in the analysis of the proposed ROA estimate.

3.3 Covariance Oriented Safe Sample Policy

Computing a high-confidence GP model is closely related to the appropriate prior information. If the mean function $m(x)$ of a GP model is closer to the prior information, $k(x,x_{*})$ will decrease to get close to the exact dynamics [24]. Therefore, handling the high covariance data will accelerate the learning process. Meanwhile, to maintain the prior information is always safe, we propose a positive sample policy inside an existing ROA as follows,

\displaystyle x^{*}=\mathop{\arg\max}\limits_{x,\psi(t;d_{\xi}(x))\in\mathcal{L}_{k}}\quad\sum\limits_{k=0}^{K}{{k}}(x,\psi(t;d_{\xi}(x))),

(13)

where $\sum k(x,\psi(t;d_{\xi}(x)))$ denotes the covariance of a trajectory $\psi(t;d_{\xi}(x))$ starting from $x\in\mathcal{L}_{k}$ , $\mathcal{L}_{k}$ denotes the estimate ROA and $k$ denotes the index number of the total episode $K$ . The trajectory $\psi(t;d_{\xi}(x^{*}))$ of the result in (13) can be used to enrich the GP prior model recursively. The example below shows a comparison of the dynamics between the original system (1) and a learned system (12).

Example $3$ : Consider a partially unknown nonlinear system (14) with unknown term $d(x)$ as

\displaystyle\begin{aligned} \begin{bmatrix}\dot{x}_{1}\\ \dot{x}_{2}\end{bmatrix}=\begin{bmatrix}-x_{1}+x_{2}\\ x_{1}^{2}x_{2}+1-\sqrt{|\exp{(x_{1})}\cos{(x_{1})}|}\end{bmatrix}+\begin{bmatrix}0\\ d(x)\end{bmatrix}.\end{aligned}

(14)

Figure 3(b) shows the learned polynomial system dynamics of (14). The non-polynomial term in (14) is approximated by $4^{th}$ degree Chebyshev interpolants in $[-2,2]^{2}$ , while the state trajectory $\psi(t;d_{\xi}(-0.05,-0.05))$ is used to construct a GP model with $3^{rd}$ degree polynomial mean function. Figure 4 shows the obtained sum of covariance of the learned polynomial system based on (13). Compared to the points around the equilibrium point, those nearby the ROA edge have larger sum of covariance and can be considered with high uncertainty. Thus, the distribution of the sum of covariance in Figure 4 poses a direction to reduce the uncertainty during the learning process.

4 Estimating the Region of Attraction

In this section, we will estimate the ROA of the learned polynomial system (12) via SOSPs . The relationship of estimates between the learned polynomial system and the true dynamics will also be discussed.

4.1 Parameterize the ROA by Barrier Function

The parameterization of a barrier function in (12) via SOSPs is developed in [13]. To deal with the non-negativity constraints, we will introduce the Positivestellensatz (P-satz) here for further SOSPs implementation.

Let $\mathcal{P}$ be the set of polynomials and $\mathcal{P}^{\text{SOS}}$ be the set of sum of squares polynomials, e.g., $P(x)=\sum_{i=1}^{k}p_{i}^{2}(x),$ where $p_{i}(x)\in\mathcal{P}$ and $P(x)\in\mathcal{P}^{SOS}$ .

Lemma 4.

([28]) For polynomials $\{a_{i}\}_{i=1}^{m}$ , $\{b_{j}\}_{j=1}^{n}$ and $p$ , define a set $\mathcal{B}=\{x\in\mathcal{X}:\{a_{i}(x)\}_{i=1}^{m}=0,\{b_{j}(x)\}_{j=1}^{n}\geq 0\}$ . Let $\mathcal{B}$ be compact, then $\forall x\in\mathcal{X},p(x)\geq 0$ holds if,

$\left\{\begin{array}[]{l}\exists r_{1},\dots,r_{m}\in\mathcal{P},\leavevmode\nobreak\ s_{1},\dots,s_{n}\in\mathcal{P}^{\text{SOS}},\\ p-\sum^{m}_{i=1}r_{i}a_{i}-\sum^{n}_{j=1}s_{j}b_{j}\in\mathcal{P}^{\text{SOS}}.\end{array}\right.$ $\square$

This lemma declared that any strictly positive polynomial $p$ is in the cone that generated by polynomials $\{a_{i}\}_{i=1}^{m}$ and $\{b_{j}\}_{j=1}^{n}$ . Using Lemma 4, the optimal barrier function $h^{*}(x)$ of the learned polynomial system (12) can be searched in the state space via 3 SOSPs.

1. Obtain the maximum sublevel set $\{x|V(x)\leq c^{*}\}$ of a specified Lyapunov function $V(x)$ by a bilinear search,

		$\displaystyle c^{*}=$	$\displaystyle\underset{c\in\mathbb{R}^{+},\,L_{c}(x)\in\mathcal{P}^{\text{SOS}}}{\text{max}}$	$\displaystyle\>\>c$		(15)
		s.t.	$\displaystyle-\frac{\partial V(x)}{\partial x}\dot{x}-L_{c}(x)(c-V(x))$	$\displaystyle\in\mathcal{P}^{\text{SOS}},$		(15)

where $L_{c}(x)$ is an auxiliary factor to relax the non-negativity constraint for the initial barrier function $h(x)$ .

2. Search another two auxiliary factors $L_{1,2}(x)$ for $h(x)$ ,

$\displaystyle\exists\;L_{1}(x),L_{2}(x)\in\mathcal{P}^{\text{SOS}}$				(16)
	s.t.	$\displaystyle-\frac{\partial V(x)}{\partial x}\dot{x}-L_{1}(x)h(x)$	$\displaystyle\in\mathcal{P}^{\text{SOS}},$
$\displaystyle\frac{\partial h(x)}{\partial x}\dot{x}-L_{2}(x)h(x)$	$\displaystyle\in\mathcal{P}^{\text{SOS}}.$

Meanwhile, $h(x)$ can be re-written into the square matrix representation form as $h(x)=Z(x)^{\mathrm{T}}QZ(x)$ , where $Q$ is a semi-definite coefficient matrix and $Z(x)$ is a monomial vector [13]. The trace $\textit{Tr}(\cdot)$ of $Q$ is used to approximate the volume of barrier certified ROA (BCROA).

3. Enlarge $\textit{Tr}(Q)$ to parameterize a permissive $h(x)$ with fixed $L_{1}(x)$ and $L_{2}(x)$ ,

$\displaystyle\underset{\begin{subarray}{c}h(x)\in\mathcal{P}\end{subarray}}{\text{max}}\quad\textit{Tr}(Q)$				(17)
	s.t.	$\displaystyle-\frac{\partial V(x)}{\partial x}\dot{x}-L_{1}(x)h(x)$	$\displaystyle\in\mathcal{P}^{\text{SOS}},$
$\displaystyle\frac{\partial h(x)}{\partial x}\dot{x}-L_{2}(x)h(x)$	$\displaystyle\in\mathcal{P}^{\text{SOS}}.$

The optimal barrier function $h^{*}(x)$ can be found with a permissive BCROA if the increase of $\textit{Tr}(Q)$ is less than a threshold, otherwise repeat 2 and 3 for a long term.

These 3 SOSPs demonstrate the process to obtain an optimal BCROA directly. Thus, if a ROA exists in the learned polynomial system (9), a permissive BCROA would be computed by these SOSPs directly.

4.2 Existence of Probabilistic BCROA

The relationship of a BCROA in the learned system (12) toward the true dynamics is given in the following theorem.

Theorem 1.

Given $k$ measurements of a partially unknown system (9) such that we can obtain a learned system (12) with a probability greater or equal to $(1-\delta)^{k},\delta\in(0,1)$ . Based on the SOSPs (15), (16) and (17), if there exists a BCROA $\mathcal{L}=\{x\in\mathcal{X}|h(x)\geq 0\}$ such that

\displaystyle\frac{\partial h(x)}{x}(f(x)+P_{k}(x)+m_{d_{\xi}}(x))\geq 0

(18)

holds for the states $x\in\mathcal{L}$ , then $\mathcal{L}$ can be regarded as a BCROA in (9) with probability bounds $((1-\delta)^{k},1)$ .

Proof.

Given an arbitrary initial state $x_{0}\in\mathcal{L}$ in the learned system (12), its state trajectory $\psi(t;x_{0}),t\in[0,\infty]$ is guaranteed to converge to the origin due to the Lyapunov constraints in (15), (16) and (17). Since the derivative of the barrier function $\frac{\partial h(\psi(t;x_{0}))}{\partial t}$ is non-negative over the trajectory $\psi(t;x_{0})$ , $h(\psi(t;x_{0}))$ will strictly increase along $\psi(t;x_{0}),t\in[0,\infty)$ . Then, the BCROA $\mathcal{L}$ established a region with the safety and stability guarantee in (12).

Because $\mathcal{L}$ denotes a certain region in (12), while (12) is an approximation result of (9) within probability bounds $((1-\delta)^{k},1)$ . Thus, $\mathcal{L}$ also denotes a potential BCROA toward the true dynamics in (9) with the same probability bounds, which ends the proof. ∎

The BCROA in Theorem 1 that satisfies (15), (16) and (17) encodes the safety of the learned polynomial system (12). It further shows a probabilistic relationship toward the real ROA in the concerned system (9) such that we can make high-probability statements about system safety.

4.3 Probabilistic Optimal BCROA Estimate

To illustrate our procedures of finding probabilistic optimal BCROAs, we present an algorithm below and a corresponding flowchart in Figure 5, which both consist of a learned polynomial system and SOSPs. Besides, the following theorem proposes the confidence range of the result from Algorithm 1 over episodes.

Input: Original system (1); initial sublevel set of a given Lyapunov function

\mathcal{L}_{0}=\{x|V(x)\leq c_{0}\}

; SOSP termination

\epsilon

and episodes

K

Output: Optimal BCROAs

R=[\mathcal{L}^{*}_{1},\mathcal{L}^{*}_{2},\dots,\mathcal{L}^{*}_{K}]^{\mathrm{T}}

; probability bounds

P=[\delta_{1},\delta_{2},\dots,\delta_{K}]^{\mathrm{T}}

1 Preprocess (1) into (9) by Chebyshev interpolants. Initialize the prior information set

D_{\xi}=\emptyset

2 Store the measurements of

d_{\xi}(x)

in (9) as

D_{0}

3 for $i\in\{1,2,\dots,K\}$ do

4 Update prior information

D_{\xi}=\cup D_{i-1}

5 Construct (12) with polynomial

m(x)

6 Execute SOSP (15) and set the initial barrier function as

h_{i}(x)=c^{*}-V(x)

7 Execute SOSPs (16) and (17) to obtain the optimal BCROA

\mathcal{L}^{*}_{i}=\{x|h^{*}_{i}(x)\geq 0\}

8 Record the probability bounds

\delta_{i}

and update the measurement of

d_{\xi}(x)

according to (13).

9return

R,P

Algorithm 1 Optimal BCROAs Estimate

Figure 5: Flowchart of Algorithm 1.

Theorem 2.

Algorithm 1 establishes $K$ probabilistic ROA $[\mathcal{L}^{*}_{1},\mathcal{L}^{*}_{2},\dots,\mathcal{L}^{*}_{K}]^{\mathrm{T}}$ with probabilities greater or equal to $[(1-\delta_{1})^{n_{1}},(1-\delta_{2})^{n_{2}},\dots,(1-\delta_{K})^{n_{K}}]^{\mathrm{T}}$ , where $n_{i}$ denotes the sample number with index $i=1,2\cdots,K$ . Then, the following probability bounds are given:

	$\displaystyle\mathbb{P}\left\{\forall i\in[1,K],\forall\mathcal{L}^{}_{i}\neq\emptyset,\mathcal{L}^{}_{\cup}\subset D\right\}$	$\displaystyle=\min\{(1-\delta_{i})^{n_{i}}\}$		(19)
	$\displaystyle\mathbb{P}\left\{\forall i\in[1,K],\forall\mathcal{L}^{}_{i}\neq\emptyset,\mathcal{L}^{}_{\cap}\subset D\right\}$	$\displaystyle=\max\{(1-\delta_{i})^{n_{i}}\},$		(19)

where $\mathcal{L}^{*}_{\cup}=\cup_{i=1}^{K}\mathcal{L}^{*}_{i}$ and $\mathcal{L}^{*}_{\cap}=\cap_{i=1}^{K}\mathcal{L}^{*}_{i}$ denote the union set and the intersection set of these optimal BCROAs, respectively, and $D$ denotes the exact ROA of the system (9).

Proof.

After executing the algorithm $K$ times, we could get optimal BCROAs $[\mathcal{L}^{*}_{1},\mathcal{L}^{*}_{2},\dots,\mathcal{L}^{*}_{K}]^{\mathrm{T}}$ with probability bounds $[[(1-\delta_{1})^{n_{1}},1),[(1-\delta_{2})^{n_{2}},1),\dots,[(1-\delta_{K})^{n_{K}},1)]^{\mathrm{T}}$ toward the dynamical system (9). Since the increasing $K$ generates larger information set $D_{\xi}$ , when $K\rightarrow\infty$ , the GP model can be inferred with the richest $D_{\xi}$ such that the lower bound of this probabilistic range is approximating $1$ . Let $\mathcal{L}^{*}_{\cap}$ be the intersection set of $\mathcal{L}^{*}_{i}$ , which represents the most common part of these optimal BCROAs toward the exact dynamics of (9). Then, a conservative statement about $\mathcal{L}^{*}_{\cap}$ toward the real ROA $D$ in (9) establishes with a minimum value of $\{(1-\delta_{i})^{n_{K}}\}$ . Therefore, the second equation in (19) about $\mathcal{L}^{*}_{\cup}$ toward $D$ establishes in a similar way, which completes the proof. ∎

Theorem 2 presents a probabilistic statement of $\mathcal{L}^{*}_{i}$ generated from Algorithm 1. Obviously, a higher degree mean function $m_{d_{\xi}}(x)$ in (12) could approximate the real dynamics much more precisely. Note that the shape of $\mathcal{L}^{*}_{i}$ is various with uncertainty. It is worthy noting that the corresponding Lyapunov function $V(x)$ and the barrier function $h(x)$ could increase their degrees simultaneously for a better estimate. This is because a higher degree polynomial certified ROA can explore more areas, so the probabilistic bound narrows naturally.

5 Numerical Examples

Based on the Matlab toolbox of Chebfun, GPML, SOSOPT and Mosek solver, we estimate the BCROA of the autonomous partially unknown system by two examples.

5.1 Example 4: A 2D Nonlinear System

Extended to Example 3, we consider the system:

\displaystyle\begin{aligned} \begin{bmatrix}\dot{x}_{1}\\ \dot{x}_{2}\end{bmatrix}=\begin{bmatrix}-x_{1}+x_{2}\\ x_{1}^{2}x_{2}+1-\sqrt{|\exp{(x_{1})}\cos{(x_{1})}|}\end{bmatrix}+\begin{bmatrix}0\\ d(x)\end{bmatrix}.\end{aligned}

(20)

Example 3 shows the construction of the first learned dynamical system in Figure 3(b) and the noise over the system measurement is bounded by $0.01$ . The learned GP models are consisted by a $2^{nd}$ degree polynomial mean function zero, and a RBF kernel function $k(x,x^{\prime})$ with signal variances $\exp(0.1)$ and length scale $\exp(0.2)$ .

At each episode, a given LCROA $\{x|\;V(x)\leq c_{0}=0.1\}$ is considered to the next step, where $V(x)=2x_{1}^{4}+0.5x_{2}^{4}-x_{1}^{2}x_{2}^{2}+x_{1}^{2}+x_{2}^{2}-0.5x_{1}x_{2}$ is a $4^{th}$ degree polynomial. Besides, we obtained the probability bounds with a fixed discounting factor $\sqrt{\beta}=2$ according to the method in [21]. After $10$ episodes, the probability bounds of these optimal BCROAs is $\mathbb{P}\{\forall\mathcal{L}_{i\in[1,10]}^{*},\mathcal{L}^{*}_{i}\subset D\}\in[0.9500,0.9700]$ , where $D$ denotes the exact ROA in (20).

In Figure 6, each computed optimal BCROA is displayed with a fixed transparency. The intersection region of these optimal BCROAs is established with the probability $[0.9700,1.0000]$ , while the union region of these optimal BCROAs is established with the probability $[0.9500,1.0000]$ . Therefore, the result of these probabilistic optimal BCROAs in Figure 6 is not only obviously larger than the LCROA, but also guaranteeing the system safety by a clear probability range.

5.2 Example 5: A 3D Nonlinear System

Consider a 3D nonlinear system corrupted with unknown terms $d_{1}(x),d_{3}(x)$ as the form of (1),

\displaystyle\begin{aligned} \begin{bmatrix}\dot{x}_{1}\\ \dot{x}_{2}\\ \dot{x}_{3}\end{bmatrix}=\begin{bmatrix}-x_{1}^{2}-\cos{(x_{1}^{2})}\sin{(x_{1})}\\ -x_{2}-x_{1}^{3}x_{2}\\ -x_{1}^{2}x_{3}+1-\sqrt{|\exp{(x_{1})}\cos{(x_{1})}|}\end{bmatrix}+\begin{bmatrix}d_{1}(x)\\ 0\\ d_{3}(x)\end{bmatrix}.\end{aligned}

(21)

The non-polynomial terms in (21) are approximated by the $4^{th}$ degree Chebyshev interpolants in $[-2,2]^{3}$ . To learn the system dynamics, we first collect the trajectory starting at $(-0.1,-0.1,0.1)$ , where the noises over $d_{1}(x)$ and $d_{3}(x)$ are bounded by $0.01$ . In this case, we use a $3^{rd}$ degree polynomial mean function and a RBF kernel with signal variances $\exp(0.01)$ and length scale $\exp(0.2)$ .

A sublevel set $\{x|\;V(x)\leq c_{0}=5\times 10^{-3}\}$ is given to construct the LCROA with a $4^{th}$ order $V(x)=10x_{1}^{4}+x_{2}^{4}+20x_{3}^{4}+2x_{1}^{2}x_{2}^{2}-4x_{3}^{2}x_{2}^{2}+3x_{1}^{2}x_{3}^{2}$ . Set $\sqrt{\beta}=2$ such that the probability bounds of $10$ optimal BCROAs can be computed as $\mathbb{P}\{\forall\mathcal{L}_{i=1}^{10},\;\mathcal{L}_{i}\subset D\}\in[0.8395,0.9567]$ .

These optimal BCROAs are displayed in Figure 7 with a fixed transparency. From Figure 7, we see that there are various differences of these optimal BCROAs. Furthermore, as seen in Table 1, with the same hyperparameters, the computation time of Example 5 is greater than Example 4 significantly. However, at the episode 10, Example 5 will first stop learning due to the high computation cost. Although higher dimensional functions find a comparable BCROA, but this is often not the case due to the fragile and changeable dynamics under uncertainty.

Table 1: Computation Time

t_{c}

[sec] of Algorithm 1.

\hhline	Episode 1	Episode 5	Episode 10
Example 4	$35.6$	$168.3$	$1752.6$
Example 5	$406.4$	$837.8$	$1329.4$
\hhline

6 Conclusion

In this work, a method is proposed to compute a barrier certified region of attraction such that the stability and safety can both be guaranteed for a class of partially unknown systems. The proposed method is built based on Chebyshev interpolants, Gaussian processes, sum-of-squares programming and a safe sample policy. The effectiveness of the proposed algorithm has been demonstrated via two numerical examples.

On the other hand, we would like to admit that a large amount of unknown terms may affect the accuracy and efficiency of proposed algorithm. Possible solutions to address this issue could be to construct an efficient sample policy [29] or to exploit local GP regression [30]. Thus, future efforts will be devoted to prior data selection in further enlarging the estimated region of attraction for partially unknown systems [13, 17].

References

[1] E. Glassman, A. L. Desbiens, M. Tobenkin, M. Cutkosky, and R. Tedrake, “Region of attraction estimation for a perching aircraft: A Lyapunov method exploiting barrier certificates,” in Proceedings of the International Conference on Robotics and Automation, pp. 2235–2242, 2012.
[2] L. Wang, E. A. Theodorou, and M. Egerstedt, “Safe learning of quadrotor dynamics using barrier certificates,” in Proceedings of the International Conference on Robotics and Automation, pp. 2460–2465, 2018.
[3] S.-C. Hsu, X. Xu, and A. D. Ames, “Control barrier function based quadratic programs with application to bipedal robotic walking,” in Proceedings of the American Control Conference, pp. 4542–4548, 2015.
[4] P. Magne, D. Marx, B. Nahid-Mobarakeh, and S. Pierfederici, “Large-Signal Stabilization of a DC-Link Supplying a Constant Power Load Using a Virtual Capacitor: Impact on the Domain of Attraction,” IEEE Transactions on Industry Applications, vol. 48, no. 3, pp. 878–887, 2012.
[5] G. Chesi, Domain of attraction: Analysis and control via SOS programming, vol. 415. Springer Science & Business Media, 2011.
[6] V. I. Zubov, Methods of AM Lyapunov and their application. P. Noordhoff, 1964.
[7] F. Blanchini, “Set invariance in control,” Automatica, vol. 35, no. 11, pp. 1747–1767, 1999.
[8] A. D. Ames, S. Coogan, M. Egerstedt, G. Notomista, K. Sreenath, and P. Tabuada, “Control barrier functions: Theory and applications,” in Proceedings of the European Control Conference, pp. 3420–3431, 2019.
[9] A. D. Ames, X. Xu, J. W. Grizzle, and P. Tabuada, “Control barrier function based quadratic programs for safety critical systems,” IEEE Transactions on Automatic Control, vol. 62, no. 8, pp. 3861–3876, 2016.
[10] P. A. Parrilo, Structured semidefinite programs and semialgebraic geometry methods in robustness and optimization. PhD thesis, California Institute of Technology, 2000.
[11] I. R. Manchester, M. M. Tobenkin, M. Levashov, and R. Tedrake, “Regions of attraction for hybrid limit cycles of walking robots,” IFAC Proceedings Volumes, vol. 44, no. 1, pp. 5801–5806, 2011.
[12] A. El-Guindy, D. Han, and M. Althoff, “Estimating the region of attraction via forward reachable sets,” in Proceedings of the American Control Conference, pp. 1263–1270, 2017.
[13] L. Wang, D. Han, and M. Egerstedt, “Permissive barrier certificates for safe stabilization using sum-of-squares,” in Proceedings of the American Control Conference, pp. 585–590, 2018.
[14] D. Han, and M. Althoff, “On estimating the robust domain of attraction for uncertain non-polynomial systems: An LMI approach,” in Proceedings of the Conference on Decision and Control, pp. 2176–2183, 2016.
[15] A. Iannelli, A. Marcos, and P. Seiler, “An equilibrium-independent region of attraction formulation for systems with uncertainty-dependent equilibria,” in Proceedings of the Conference on Decision and Control, pp. 725–730, 2018.
[16] J. Vinogradska, B. Bischoff, D. Nguyen-Tuong, and J. Peters, “Stability of controllers for Gaussian process dynamics,” The Journal of Machine Learning Research, vol. 18, no. 1, pp. 3483–3519, 2017.
[17] J. Umlauft, L. Pöhler, and S. Hirche, “An uncertainty-based control Lyapunov approach for control-affine systems modeled by Gaussian process,” IEEE Control Systems Letters, vol. 2, no. 3, pp. 483–488, 2018.
[18] F. Berkenkamp, R. Moriconi, A. P. Schoellig, and A. Krause, “Safe learning of regions of attraction for uncertain, nonlinear systems with Gaussian processes,” in Proceedings of the Conference on Decision and Control, pp. 4661–4666, 2016.
[19] J. Umlauft, A. Lederer, and S. Hirche, “Learning stable Gaussian process state space models,” in Proceedings of the American Control Conference, pp. 1499–1504, 2017.
[20] A. Devonport, H. Yin, and M. Arcak, “Bayesian safe learning and control with sum-of-squares analysis and polynomial kernels,” in Proceedings of the Conference on Decision and Control, pp. 3159–3165, 2020.
[21] P. Jagtap, G. J. Pappas, and M. Zamani, “Control barrier functions for unknown nonlinear systems using Gaussian processes,” in Proceedings of the Conference on Decision and Control, pp. 3699–3704, 2020.
[22] D. Han, and M. Althoff, “Estimating the domain of attraction based on the invariance principle,” in Proceedings of the Conference on Decision and Control, pp. 5569–5576, 2016.
[23] H. K. Khalil, Nonlinear systems; 3rd ed. Prentice-Hall, 2002.
[24] C. K. Williams and C. E. Rasmussen, Gaussian processes for machine learning, vol. 2. MIT press, 2006.
[25] V. I. Paulsen and M. Raghupathi, An introduction to the theory of reproducing kernel Hilbert spaces, vol. 152. Cambridge university press, 2016.
[26] L. N. Trefethen, Approximation Theory and Approximation Practice. SIAM, 2019.
[27] N. Srinivas, A. Krause, S. M. Kakade, and M. Seeger, “Gaussian process optimization in the bandit setting: No regret and experimental design,” arXiv preprint: 0912.3995, 2009.
[28] M. Putinar, “Positive polynomials on compact semi-algebraic sets,” Indiana University Mathematics Journal, vol. 42, no. 3, pp. 969–984, 1993.
[29] J. Umlauft, T. Beckers, A. Capone, A. Lederer, and S. Hirche, “Smart forgetting for safe online learning with Gaussian Processes,” in Learning for Dynamics and Control, pp. 160–169, PMLR, 2020.
[30] D. Nguyen-Tuong, M. Seeger, and J. Peters, “Model learning with local Gaussian Process regression,” Advanced Robotics, vol. 23, no. 15, pp. 2015–2034, 2009.
[31] C. Santoyo, M. Dutreix, and S. Coogan, “A barrier function approach to finite-time stochastic system verification and control,” Automatica, vol. 125, p. 109439, 2021.
[32] H. König, Eigenvalue distribution of compact operators, vol. 16. Birkhäuser, 2013.

Appendix: Proof of Proposition 1

Proof.

First, Mercer’s Theorem [32] allows us to represent a kernel by a kernel basis vector $\phi(\cdot)$ . For this paper, we restrict our attention to the RBF kernel $k(x,x^{\prime})$ , though our methodology can support more general kernels,

$\displaystyle k(x,x^{\prime})$	$\displaystyle=\sigma^{2}\exp(\frac{(x-x^{\prime})^{2}}{-2l^{2}})$	(22)
	$\displaystyle=\sum\limits_{i=0}^{\infty}\frac{\sigma^{2}2^{i}}{2l^{2}i!}\exp(\frac{(x)^{2}}{-2l^{2}})(\frac{x}{2l^{2}})^{i}\exp(\frac{(x^{\prime})^{2}}{-2l^{2}})(\frac{x^{\prime}}{2l^{2}})^{i}$
	$\displaystyle=\phi(x)^{\mathrm{T}}\phi(x^{\prime}),$

where $\sigma$ denotes the signal variance, $l$ denotes the length scale and $\phi(x)$ denotes an infinite dimensional vector,

\displaystyle\phi(x)=\exp(\frac{(x)^{2}}{-2l^{2}})\cdot(\sqrt{\frac{\sigma^{2}2^{0}}{2l^{2}0!}}(\frac{x}{2l^{2}})^{0},\sqrt{\frac{\sigma^{2}2^{1}}{2l^{2}1!}}(\frac{x}{2l^{2}})^{1},\dots).

(23)

Then, since $d(x)$ and its measurement $y(x)$ in (1) are bounded in the RKHS, we can formulate,

\displaystyle d(x)=\phi(x)^{\mathrm{T}}w,\quad y(x)

\displaystyle=d(x)+\epsilon,

(24)

where $w$ denotes the weight vector and $\epsilon$ denotes the noise. Let $p(w)\sim\mathcal{N}(0,\sigma_{p}^{2})$ be the prior distribution of $w$ and let $p(w|y,x)$ be the posterior distribution of $w$ based on $y$ and $x$ , we can infer their relationship as follows,

\small\small\text{posterior}=\frac{\text{likelihood}\times\text{prior}}{\text{marginal\;likelihood}},p{(w|y,x)}=\frac{p(y|x,w)\cdot p(w)}{p(y|x)}.

(25)

Next, let $K_{1}$ be the weight-independent marginal likelihood $p(y|x)$ . GP is relocating likelihood $p(y|x,w)$ with measurements $[(x_{1},y_{1}),(x_{2},y_{2}),\dots,(x_{k},y_{k})]^{\mathrm{T}}$ as below,

\displaystyle p(y|x,w)

\displaystyle=\prod_{i=1}^{k}p\left(y_{i},x_{i},w\right)=\prod^{k}_{i=1}\frac{\exp{(\frac{(y_{i}-\phi(x_{i})^{\mathrm{T}}w)^{2}}{-2\sigma_{n}^{2}})}}{\sigma_{n}\sqrt{2\pi}},

(26)

to compute the $p{(w|y,x)}$ based on (25) and (26),

	$\displaystyle p(w\|x,y)$	$\displaystyle=\frac{\prod\limits_{i=1}^{k}p\left(y_{i},x_{i},w\right)\cdot\frac{1}{\sqrt{2\pi\sigma^{2}_{p}}}\exp{(\frac{w^{\mathrm{T}}\sigma_{p}^{-2}w}{-2})}}{K_{1}}$		(27)
		$\displaystyle=\!K_{2}\cdot\exp{(\frac{w^{\mathrm{T}}\sigma_{p}^{-2}w}{-2})}\cdot\exp{(\frac{\sum\limits_{i=1}^{k}(y_{i}-\phi(x_{i})^{\mathrm{T}}w)^{2}}{-2\sigma_{n}^{2}})},$		(27)

where $K_{2}=\frac{1}{K_{1}\cdot\sqrt{2\pi\sigma_{p}^{2}}\cdot(\sigma_{n}\sqrt{2\pi})^{k}}$ .

By using the Maximum a Posterior estimation, the exact value of the optimal weight $\hat{w}$ can be computed as,

$\displaystyle\hat{w}$	$\displaystyle=\arg\max_{w}\;p(w\|x,y)=\arg\max_{w}\log{p(w\|x,y)}$	(28)
	$\displaystyle=\arg\min_{w}\;-\log{p(w\|x,y)}$
	$\displaystyle=\arg\min_{w}\;-\log{K_{2}}+\frac{w^{\mathrm{T}}\sigma_{p}^{-2}w}{2}+\frac{\sum\limits_{i=1}^{k}(y_{i}-\phi(x_{i})^{\mathrm{T}}w)^{2}}{2\sigma_{n}^{2}}$
	$\displaystyle=\arg\min_{w}\;w^{\mathrm{T}}(\sigma_{n}^{2}\sigma_{p}^{-2})w+\sum\limits_{i=1}^{k}(y_{i}-\phi(x_{i})^{\mathrm{T}}w)^{2}$
	$\displaystyle=\arg\min_{w}\;w^{\mathrm{T}}(\phi(x)^{\mathrm{T}}\phi(x)+\sigma_{n}^{2}\sigma_{p}^{-2})w-2w^{\mathrm{T}}\phi(x)^{\mathrm{T}}y+y^{\mathrm{T}}y$
	$\displaystyle=\arg\min_{w}\;L(w).$

The derivative of $L(w)$ in (28) can be computed as,

\displaystyle\frac{\partial L(w)}{\partial w}

\displaystyle=2(\phi(x)^{\mathrm{T}}\phi(x)+\sigma_{n}^{2}\sigma_{p}^{-2})w-2\phi(x)^{\mathrm{T}}y,

(29)

and the optimal weight $\hat{w}$ can be parameterized directly,

$\displaystyle\hat{w}$	$\displaystyle=(\phi(x)^{\mathrm{T}}\phi(x)+\sigma_{n}^{2}\sigma_{p}^{-2})^{-1}\phi(x)^{\mathrm{T}}y$	(30)
	$\displaystyle=\phi(x)^{\mathrm{T}}(\phi(x)\phi(x)^{\mathrm{T}}+\sigma_{n}^{2}\sigma_{p}^{-2})^{-1}y$
	$\displaystyle=\sigma_{p}^{2}\phi(x)^{\mathrm{T}}(\phi(x)\sigma_{p}^{2}\phi(x)^{\mathrm{T}}+\sigma_{n}^{2})^{-1}y,$

Thus, we can obtain the optimal $p(\hat{w}|x,y)$ as (27) shows,

	$\displaystyle p(\hat{w}\|x,y)\sim($	$\displaystyle\sigma_{p}^{2}\phi(x)^{\mathrm{T}}(\phi(x)\sigma_{p}^{2}\phi(x)^{\mathrm{T}}+\sigma_{n}^{2})^{-1}y,$		(31)
		$\displaystyle(\phi(x)\sigma_{p}^{2}\phi(x)^{\mathrm{T}}+\sigma_{n}^{2})^{-1}),$		(31)

where the predictive output $d(x_{*})$ of the query states $x_{*}$ is,

\displaystyle d(x_{*})

\displaystyle=\phi(x_{*})^{\mathrm{T}}\hat{w}.

(32)

Because $\phi(\cdot)$ is an infinite vector, which is used to approximate $y$ with optimal weight $\hat{w}$ . To avoid the loss of generality, we can truncate $\phi(x)$ by a finite monomial vector $\varphi(x)$ without any impact of the inference processes from (26) to (31). This will reshape the mean function in the polynomial form and adjust the second term in (31) as,

	$\displaystyle m(x_{*})$	$\displaystyle=\varphi(x_{*})^{\mathrm{T}}\hat{w},$		(33)
	$\displaystyle\sigma^{2}(x_{*})$	$\displaystyle=k(x_{},x_{})-k_{}^{\mathrm{T}}(K+\sigma_{n}^{2}I)^{-1}k_{}.$		(33)

where $[K]_{(i,j)}=k(x_{i},x_{j})$ is a kernel Gramian matrix and $k_{*}=[k(x_{1},x_{*}),k(x_{2},x_{*}),\dots,k(x_{w},x_{*})]^{\mathrm{T}}$ . Thus, Proposition 1 is obtained, which completes the proof. ∎