Using SOS and Sublevel Set Volume Minimization for Estimation of Forward Reachable Sets

TL;DR

This paper introduces a convex Sum-of-Squares optimization method to compute tight outer approximations of forward reachable sets for nonlinear uncertain ODEs, minimizing volume via a determinant-like objective.

Contribution

It presents a novel convex SOS-based approach for volume minimization in reachable set estimation, applicable to nonlinear systems with bounded disturbances.

Findings

Effective approximation of Lorenz system's reachable set

Accurate estimation for Van der Pol oscillator

Demonstrates tight bounds with volume minimization

Abstract

In this paper we propose a convex Sum-of-Squares optimization problem for finding outer approximations of forward reachable sets for nonlinear uncertain Ordinary Differential Equations (ODE's) with either (or both) L2 or point-wise bounded input disturbances. To make our approximations tight we seek to minimize the volume of our approximation set. Our approach to volume minimization is based on the use of a convex determinant-like objective function. We provide several numerical examples including the Lorenz system and the Van der Pol oscillator.

Figures

Equations

$$\dot{x}(t)=f(x(t),u(t),w(t)),u\in U,w\in W,x(0)=x_0\in {\mathbb{R}}^n,$$

$$\dot{x}(t)=f(x(t),u(t),w(t)),u\in U_Y,w\in W_{\gamma },x(0)=x_0\in {\mathbb{R}}^n,$$

$$FR(X_0,f,T,Y,\gamma ):=\{y\in {\mathbb{R}}^n$$

$$\text{such that }{\varphi }_f(x,u,w,T)=y\},$$

$$\underset{X\in C}{\min }\{$$

$$\text{subject to: }FR(X_0,f,T,Y,\gamma )\subseteq X$$

$$V(x,0)\le 1\text{ for all }x\in X_0.$$

$$\frac{\partial V}{\partial t}(x,t)+{\nabla }_{x}V(x,t{)}^{T}f(x,u,w)\le w^{T}w,$$

$$\text{ for all }x\in X_c,t\in [0,T],u\in Y,w\in {\mathbb{R}}^{m_w}.$$

$${\int }_0^T\frac{d}{dt}$$

$$\text{ for all }{x}_0\in X_0,u\in U,w\in W.$$

$$V(\varphi (x_0,u,w,T),T)$$

$$=V(x_0,0)+\gamma$$

$$\le 1+\gamma .$$

$$\underset{X}{\min }\{$$

$$\text{subject to: }X=\{x\in {\mathbb{R}}^n:V(x,T)\le 1+\gamma \}$$

$$V(x,0)\le 1\text{ for all }x\in X_0.$$

$$\frac{\partial V}{\partial t}(x,t)+{\nabla }_{x}V(x,t{)}^{T}f(x,u,w)\le w^{T}w,$$

$$\text{ for all }x\in X_c,t\in [0,T],u\in Y,w\in {\mathbb{R}}^{m_w}.$$

$$\underset{X}{\min }\{D(X,FR(X_0,f,T,Y,\gamma ))\}$$

$$\text{subject to: }X=\{x\in {\mathbb{R}}^n:V(x,T)\le 1+\gamma \}$$

$$V\in \sum _{SOS}{k}_i\in \sum _{SOS}\text{ for }i=1,2s_i\in \sum _{SOS}\text{ for }i=1,2,3,4,$$

$$\underset{P(T)\in S_{++}^N}{\min }\{-\log \det \{P(T)\}\}$$

$$\text{subject to: }V(x,t)=z_d(x{)}^{T}P(t)z_d(x)$$

$$k_i\in \sum _{SOS}\text{ for }i=1,2s_i\in \sum _{SOS}\text{ for }i=1,2,3,4,$$

$$V(x,0)\le 1\text{ for all }x\in \{x\in {\mathbb{R}}^n:g_X(x)\ge 0\}.$$

$$\frac{\partial V}{\partial t}(x,t)+{\nabla }_{x}V(x,t{)}^{T}f(x,u,w)\le w^{T}w,$$

$$\text{for all }x\in \{y\in {\mathbb{R}}^n:g_C(y)\ge 0\},t\in [0,T],$$

$$u\in \{z\in {\mathbb{R}}^{m_u}:g_U(z)\ge 0\},w\in {\mathbb{R}}^{m_w}.$$

$$\dot{x}(t)=f(x(t)),x(0)=x_0\in {\mathbb{R}}^n,$$

$$FR(X_0,f,T):=\{y\in {\mathbb{R}}^n$$

$$\underset{P(T)\in S_{++}^N}{\min }\{-\log \det \{P(T)\}\}$$

$$\text{subject to: }V(x,t)=z_d(x{)}^{T}P(t)z_d(x)$$

$$k_1\in \sum _{SOS},k_2\in \sum _{SOS},s_i\in \sum _{SOS}\text{ for }i=1,2,3,$$

$$k_1(x)$$

$$k_2(x,t)$$

$$-s_2(x,t)g_C(x)-s_3(x,t)h(t).$$

$$\dot{x}(t)=Ax(t),$$

$${\dot{x}}_1(t)$$

$${\dot{x}}_2(t)$$

$${\dot{x}}_3(t)$$

$${\dot{x}}_1(t)$$

Full text

††thanks: This research was carried out with the financial support of the NSF Grant CNS-1739990.

1 Introduction

In this paper we consider nonlinear dynamical system described by Ordinary Differential Equations (ODE’s) of the form

[TABLE]

where $f:\mathbb{R}^{n}\times\mathbb{R}^{m_{u}}\times\mathbb{R}^{m_{w}}\to\mathbb{R}^{n}$ is a vector field; $u:\mathbb{R}\to\mathbb{R}^{m_{u}}$ and $w:\mathbb{R}\to\mathbb{R}^{m_{w}}$ are inputs; $U\subset L_{2}^{m_{u}}[0,T]$ and $W\subset L_{2}^{m_{w}}[0,T]$ are the sets of admissible inputs; and $x_{0}\in\mathbb{R}^{n}$ is the initial condition. The objective of this paper is to compute the set of states reachable from uncertain initial conditions contained inside some semialgebriac set under bounded $L_{2}$ and point-wise bounded input disturbances.

Computing reachable sets is important practically for certifying systems remain in “safety regions”; regions of the state space that are deemed to have low risks of system failure. Reachable set analysis allows for the possibility of systems to detect impending transitions outside of “safe regions” and then execute control laws to avoid such transitions. Historic examples of system’s transitioning outside “safe regions” include: two of the four reaction wheels on the Kepler Space telescope failing, analyzed in Kampmeier et al. (2018); and the disturbing lateral vibrations of the Millennium footbridge over the River Thames in London on opening day, analyzed in Chen et al. (2018) and Eckhardt et al. (2007). Such examples can potentially be avoided by modeling the system by an ODE of the Form (1) and analyzing under what input disturbances potential failure can occur.

There are many different approaches to reachability analysis. One such approach is simulation based methods; here solution maps are estimated and algorithms are designed such that the numerical error always results in an over approximation of the reachable set. Such simulation methods are explored in the works of Greenstreet and Mitchell (1999) where nonlinear dynamics are approximated by linear dynamics, where solution maps can be analytically found. In Li et al. (2018) solution maps are Taylor expanded to over or under approximate the reachable set of autonomous systems. A further alternative simulation based method can be found in Maidens and Arcak (2015). Another approach to reachability analysis is to construct a particular Hamilton-Jacobi-Isaac PDE that has a viscosity solution such that its zero sublevel set is the reachable set as shown in Mitchell et al. (2005).

Our approach to reachability analysis is to outer approximate the forward reachable set by the sublevel set of a function satisfying energy like disspataion inequalities. To make the outer approximation tight we would like to minimize the volume of our outer set approximation. This problem formulation is also found in Yin et al. (2018). However here a heuristic bisection method is proposed where the sublevel set $\{x\in\mathbb{R}^{n}:q(x)\leq\alpha\}$, where $q:\mathbb{R}^{n}\to\mathbb{R}$ is some hand picked function, is constrained to contain the outer approximation and $\alpha>0$ is iteratively minimized. In our previous work, Jones and Peet (2018), it was argued that the volume of a sublevel set of the form $\{x\in\mathbb{R}^{n}:z_{d}(x)^{T}Pz_{d}(x)\leq 1\}$, where $z_{d}(x)$ is the monomial vector of degree $d\in\mathbb{N}$ and $P$ is a positive matrix, can be minimized by minimizing the convex objective function $-\log\det P$. In this work we will therefore formulate a convex optimization problem who’s solution can construct a tight outer approximation of the forward reachable set of an ODE of Form (1) using a $\log\det$ type objective function. Our reachable set analysis does not require the use of bisection methods or handpicked functions and furthermore is formulated to include both $L_{2}$ bounded type input disturbances and point wise bounded type input disturbances.

The rest of this paper is organized as follows. In Section 3 we propose an optimization whose solution can construct an optimal outer set approximation of a forward reachable set of an ODE. In Section 4 we propose a class of functions that satisfy energy like dissipation inequalities and show how such functions can help us characterize sets that contain forward reachable sets. In Section 5 a convex Sum-of-Squares optimization problem, based upon these energy dissipating functions, is then proposed whose solution can construct an outer approximation of the forward reachable set. Here we justify the tightness of the outer approximation by the minimization of a $\log\det$ objective function. In Section 6 we give several numerical examples of our outer approximations of forward reachable sets. Finally we give our conclusion in 7.

2 Notation

We denote the set $L_{2}^{m}[0,T]:=\{g:\mathbb{R}\to\mathbb{R}^{m}:\int_{0}^{T}g(t)^{T}g(t)dt<\infty\}$. For $Y\subset\mathbb{R}^{m_{u}}$ we denote $U_{Y}=\{u\in L_{2}^{m_{u}}[0,T]:u(t)\in Y\subset\mathbb{R}^{m_{u}}\text{ for all }t\in[0,T]\}$. For $\gamma>0$ we denote $W_{\gamma}=\{w\in L_{2}^{m_{w}}[0,T]:\int_{0}^{T}w(t)^{T}w(t)dt<\gamma\}$. For a set $A\subset\mathbb{R}^{n}$ we define the indicator function $\mathds{1}_{A}:\mathbb{R}^{n}\to\mathbb{R}$ by $\mathds{1}_{A}(x)=\begin{cases}&1\text{ if }x\in A\\ &0\text{ otherwise}\end{cases}$. For a set $A\subset\mathbb{R}^{n}$ we define $vol\{A\}=\int_{\mathbb{R}^{n}}\mathds{1}_{A}(x)dx$. We denote the power set of $\mathbb{R}^{n}$ by $P(\mathbb{R}^{n})=\{X:X\subset\mathbb{R}^{n}\}$. For a differentiable function $V:\mathbb{R}^{n}\times\mathbb{R}\to\mathbb{R}$ we denote $\nabla_{x}V=(\frac{\partial V}{\partial x_{1}},....,\frac{\partial V}{\partial x_{n}})$. For two sets $A,B\subset\mathbb{R}^{n}$ we denote $A/B=\{x\in A:x\notin B\}$. We denote $S^{n}_{++}$ to be the set of positive definite $n\times n$ matrices. For $x\in\mathbb{R}^{n}$ we denote $z_{d}(x)$ to be the vector of monomial basis in $n$-dimensions with maximum degree $d\in\mathbb{N}$. We say the polynomial $p:\mathbb{R}^{n}\to\mathbb{R}$ is Sum-of-Squares (SOS) if there exists polynomials $p_{i}:\mathbb{R}^{n}\to\mathbb{R}$ such that $p(x)=\sum_{i=1}^{k}(p_{i}(x))^{2}$. We denote $\sum_{SOS}$ to be the set of SOS polynomials.

3 An Optimization Problem For Minimal Outer Bounds of Reachable Sets

In this paper we consider systems that can be modeled by nonlinear Ordinary Differential Equations (ODE’s) of the form

[TABLE]

where $f:\mathbb{R}^{n}\times\mathbb{R}^{m_{u}}\times\mathbb{R}^{m_{w}}\to\mathbb{R}^{n}$ is a vector field; $u:\mathbb{R}\to\mathbb{R}^{m_{u}}$ and $w:\mathbb{R}\to\mathbb{R}^{m_{w}}$ are inputs; $U_{Y}=\{u\in L_{2}^{m_{u}}[0,T]:u(t)\in Y\subset\mathbb{R}^{m_{u}}\text{ for all }t\in[0,T]\}$ and $W_{\gamma}=\{w\in L_{2}^{m_{w}}[0,T]:\int_{0}^{T}w(t)^{T}w(t)dt<\gamma\}$ are the sets of admissible inputs; $Y\subset\mathbb{R}^{m_{u}}$; $\gamma>0$; and $x_{0}\in\mathbb{R}^{n}$ is the initial condition. Typically inputs that are members of the set $U_{Y}$ are thought of as uncertainties and inputs that are members of the set $W_{\gamma}$ are thought of as disturbances.

Throughout this paper we will assume the existence and uniqueness of solution maps.

Definition 1

We say $\phi_{f}:\mathbb{R}^{n}\times L_{2}^{m_{u}}[0,T]\times L_{2}^{m_{w}}[0,T]\times\mathbb{R}^{+}\to\mathbb{R}^{n}$ is the solution map for (2) if $\frac{\delta\phi_{f}(x,u,w,t)}{\delta t}=f(\phi_{f}(x,u,w,t),u(t),w(t))$ and $\phi_{f}(x,u,w,0)=x$.

The goal of this paper is to estimate the set of coordinates in $\mathbb{R}^{n}$ that the solution map can attain, at some finite time $T\geq 0$, starting in some set of initial conditions; we call this set the forward reachable set and define it formally next.

Definition 2

For an ODE of Form (2), $X_{0}\subset\mathbb{R}^{n}$, $T>0$, $Y\subset\mathbb{R}^{m_{u}}$, and $\gamma\geq 0$ we define the forward reachable set of $X_{0}$ at time $T>0$ by

[TABLE]

where $U_{Y}=\{u\in L_{2}^{m_{u}}[0,T]:u(t)\in Y\subset\mathbb{R}^{m_{u}}\text{ for all }t\in[0,T]\}$ and $W_{\gamma}=\{w\in L_{2}^{m_{w}}[0,T]:\int_{0}^{T}w(t)^{T}w(t)dt<\gamma\}$.

Lemma 3

Suppose $X_{1},X_{2}\subset\mathbb{R}^{n}$ that are such that $X_{1}\subseteq X_{2}$. Then $FR(X_{1},f,T,Y,\gamma)\subseteq FR(X_{2},f,T,Y,\gamma)$, where $f:\mathbb{R}^{n}\times\mathbb{R}^{m_{u}}\times\mathbb{R}^{m_{w}}\to\mathbb{R}^{n}$, $T>0$, $Y\subset\mathbb{R}^{m_{u}}$, and $\gamma\geq 0$.

{pf}

Suppose $y\in FR(X_{1},f,T,Y,\gamma)$, then there exists $x_{0}\in X_{1}$, $u\in U_{Y}$ and $w\in W_{\gamma}$ such that $\phi_{f}(x_{0},u,w,T)=y$. Since $X_{1}\subseteq X_{2}$ we have $x_{0}\in X_{2}$. Therefore it follows $y\in FR(X_{2},f,T,Y,\gamma)$. Since $y$ was arbitrarily chosen we deduce $FR(X_{1},f,T,Y,\gamma)\subseteq FR(X_{2},f,T,Y,\gamma)$.

For $X_{0}\subset\mathbb{R}^{n}$, $f:\mathbb{R}^{n}\times\mathbb{R}^{m_{u}}\times\mathbb{R}^{m_{w}}\to\mathbb{R}^{n}$, $T>0$, $Y\subset\mathbb{R}^{m_{u}}$, and $\gamma\geq 0$, we now propose the following optimization problem to find the optimal outer set approximation, that is an element of some set $C$, of a reachable set.

[TABLE]

where $C\subset P(\mathbb{R}^{n})$ and $D:P(\mathbb{R}^{n})\times P(\mathbb{R}^{n})\to\mathbb{R}$ is some metric that measures the distance between two subsets of $\mathbb{R}^{n}$.

When solving the above Optimization Problem (3) there are two challenges.

To enforce the constraint $FR(X_{0},f,T,Y,\gamma)\subseteq X$. 2. 2.

To select a metric $D$ that can be tractably minimized.

In the next section we tackle the first of these challenges.

4 Sublevel Sets Of Functions Satisfying Dissipation Like Inequalities Contain Reachable Sets

Reachable sets are implicitly defined using solution maps of ODE’s. Therefore the set containment in Optimization Problem (3) must be indirectly constrained. To enforce the set containment constraint we use energy-like dissipation inequalities. In the next theorem we will show that if there exists a function, that has a rate of change along the solution map less than the magnitude of the $L_{2}$ bounded input, $w^{T}w$, for any point wise admissible input, then it has a sublevel set at time $T>0$ that must contain the forward reachable set at time $T$.

Theorem 4

For some $X_{0}\subset\mathbb{R}^{n}$, $f:\mathbb{R}^{n}\times\mathbb{R}^{m_{u}}\times\mathbb{R}^{m_{w}}\to\mathbb{R}^{n}$, $T>0$, $Y\subset\mathbb{R}^{m_{u}}$, and $\gamma\geq 0$, suppose there exists a function $V:\mathbb{R}^{n}\times\mathbb{R}\to\mathbb{R}$ such that

[TABLE]

Then $FR(X_{0},f,T,Y,\gamma)\subseteq\{x\in\mathbb{R}^{n}:V(x,T)\leq 1+\gamma\}$ and where $X_{c}\subset\mathbb{R}^{n}$ is any set such that $FR(X_{0},f,t,Y,\gamma)\subseteq X_{c}$ for all $t\in[0,T]$ (typically we take $X_{c}=\mathbb{R}^{n}$).

{pf}

Since $FR(X_{0},f,t,Y,\gamma)\subseteq X_{c}$ for all $t\in[0,T]$ we have $\phi_{f}(x_{0},u,w,t)\in X_{c}$ for all $x_{0}\in X_{0}$, $u\in U_{Y}$ and $w\in W_{\gamma}$. Now using (5) and the $L_{2}$ bound on $w(t)$, it follows

[TABLE]

Thus we deduce from rearranging (6) and using (4) that $\text{for all }x_{0}\in X_{0},u\in U,w\in W$

[TABLE]

Now clearly from (7) we have $\phi_{f}(x_{0},u,w,T)\in\{x\in\mathbb{R}^{n}:V(x,T)\leq 1+\gamma\}$ for all $x_{0}\in X_{0}$, $u\in U_{Y}$ and $w\in W_{\gamma}$. Therefore $FR(X_{0},f,T,Y,\gamma)\subseteq\{x\in\mathbb{R}^{n}:V(x,T)\leq 1+\gamma\}$.

The set $X_{c}\subset\mathbb{R}^{n}$ can be thought of as the computation region. In general we can select $X_{c}=\mathbb{R}^{n}$ and $FR(X_{0},f,t,Y,\gamma)\subseteq X_{c}$ will always be satisfied, however setting $X_{c}$ to be some sufficiently large bounded set can result in better numerical results. This is because we will later use Semidefinite Programing (SDP) to find polynomial functions that satisfy the inequalities (4) (5). To the authors knowledge there is no converse theorem that proves the existence of such a polynomial function satisfying these inequalities, however we do know that the Weierstrass approximation theorem states that any continuous function can be uniformly approximated over a closed and bounded set by a polynomial function. In light of this result and from numerical experience the authors recommend the use of bounded computation regions.

We now use Theorem 4 and Optimization Problem (3) to write an optimization problem with a solution that can construct an outer approximation of the reachable set.

[TABLE]

5 Proposing A Convex SOS Optimization Problem For Reachable Set Approximation

Currently the objective function of Optimization Problem (8) is said to be a metric that measures the distance between sets and has not yet been defined exactly; we will tackle this problem later. Firstly we consider the problem of enforcing the constraints of this optimization problem, which currently are not tractable.

To solve the Optimization Problem (8) we must find a function, $V$, that satisfies several inequality constraints. Determining whether a polynomial satisfies an inequality constraint has the same difficulties as proving a polynomial is globally positive ($f(x)>0$ $\forall x\in\mathbb{R}^{n}$); Blum et al. (1998) has shown this problem to be NP-hard. However it can be shown testing if a polynomial is Sum-of-Squares (SOS), and hence positive, is equivalent to solving a semidefinite program (SDP). Although not all positive polynomials are SOS, this gap can be made arbitrarily small, see Hilbert (1888). We thus propose a tightening of the optimization problem and restrict the decision variable, $V$, to be an SOS polynomial.

We now propose a tightened SOS optimization problem of Optimization Problem (8). To do this will assume the existence of polynomial functions, $g_{X}$, $g_{C}$, and $g_{U}$, such that $X_{0}\subseteq\{x\in\mathbb{R}^{n}:g_{X}(x)\geq 0\}$, $X_{c}\subseteq\{x\in\mathbb{R}^{n}:g_{C}(x)\geq 0\}$ and $Y\subseteq\{u\in\mathbb{R}^{m_{u}}:g_{U}(u)\geq 0\}$. To ensure the hypothesis of Theorem 4, $FR(X_{0},f,t,Y,\gamma)\subseteq X_{c}$ for all $t\in[0,T]$, is satisfied we typically select $g_{c}(x)=R^{2}-||x||_{2}^{2}$ where $R>0$ can be made sufficiently large. Lastly we denote the function $h(t)=t[T-t]$ and note the problems time interval can be described as $[0,T]=\{t\in\mathbb{R}:h(t)\geq 0\}$.

[TABLE]

where $k_{1}(x)=(1-V(x,0))-s_{1}(x)g_{X}(x)$ and $k_{2}(x,u,w,t)=-\left(\frac{\partial V}{\partial t}(x,t)+\nabla_{x}V(x,t)^{T}f(x,u,w)-w^{T}w\right)-s_{2}(x,u,w,t)g_{C}(x)$ $-s_{3}(x,u,w,t)h(t)-s_{4}(x,u,w,t)g_{U}(u)$.

To make the above optimization problem (9) tractable we must select a metric $D$ that is convex and hence can be minimized numerically. In our previous work, Jones and Peet (2018), it was shown that in the case where the metric is $D_{V}(X,Y)=vol\{(X/Y)\cup(Y/X)\}$ a heuristic solution to the above optimization problem can be found by using a $\log\det$ type convex objective function. We now therefore propose a convex SOS optimization problem; that as shown in Proposition 5 is solved by a feasible, and in general suboptimal, solution to the intractable Optimization Problem (3). The optimization problem is denoted by $S_{1}(d,T,f,\gamma,g_{X},g_{C},g_{U},h)$:

[TABLE]

where $k_{1}(x)=(1-V(x,0))-s_{1}(x)g_{X}(x)$ and $k_{2}(x,u,w,t)=-\left(\frac{\partial V}{\partial t}(x,t)+\nabla_{x}V(x,t)^{T}f(x,u,w)-w^{T}w\right)-s_{2}(x,u,w,t)g_{C}(x)$ $-s_{3}(x,u,w,t)h(t)-s_{4}(x,u,w,t)g_{U}(u)$.

Proposition 5

For some $X_{0}\subset\mathbb{R}^{n}$, $f:\mathbb{R}^{n}\times\mathbb{R}^{m_{u}}\times\mathbb{R}^{m_{w}}\to\mathbb{R}^{n}$, $T>0$, $Y\subset\mathbb{R}^{m_{u}}$, and $\gamma\geq 0$, suppose there exists functions $g_{X}:\mathbb{R}^{n}\to\mathbb{R}$, $g_{U}:\mathbb{R}^{m_{u}}\to\mathbb{R}$, and $g_{C}:\mathbb{R}^{n}\to\mathbb{R}$ such that $X_{0}\subseteq\{x\in\mathbb{R}^{n}:g_{X}(x)\geq 0\}$, $Y\subseteq\{u\in\mathbb{R}^{m_{u}}:g_{U}(u)\geq 0\}$ and $FR(X_{0},f,t,Y,\gamma)\subseteq\{x\in\mathbb{R}^{n}:g_{C}(x)\geq 0\}$ for all $t\in[0,T]$. Then if $P(t):\mathbb{R}\to\mathbb{R}$ solves the problem $S_{1}(d,T,f,\gamma,g_{X},g_{C},g_{U},h)$, found in (10), for $h(t)=t[T-t]$ and some $d\in\mathbb{N}$ then $FR(X_{0},f,T,Y,\gamma)\subseteq\{x\in\mathbb{R}^{n}:z_{d}(x)P(T)z_{d}(x)\leq 1+\gamma\}$.

{pf}

We will show $V(x,t)=z_{d}(x)P(t)z_{d}(x)$ satisfies inequalities (4) and (5) in order to use Theorem 4.

From the constraints of $S_{1}(d,T,f,\gamma,g_{X},g_{C},g_{U},h)$ we have $k_{1}\in\sum_{SOS}$ and $s_{1}\in\sum_{SOS}$ and thus it follows $k_{1}(x)=(1-V(x,0))-s_{1}(x)g_{X}(x)\geq 0$ and $s_{1}(x)\geq 0$ for all $x\in\mathbb{R}^{n}$. Since a positive function multiplied with a positive function is also a positive function we can now deduce

[TABLE]

Moreover the above inequality also holds for all $x\in X_{0}$ as $X_{0}\subseteq\{x\in\mathbb{R}^{n}:g_{X}(x)\geq 0\}$. Furthermore using a similar argument with the remaining constraints of $S_{1}(d,T,f,\gamma,g_{X},g_{C},g_{U},h)$ we can also deduce

[TABLE]

Moreover the above inequality also holds for all $u\in Y$ as $Y\subseteq\{u\in\mathbb{R}^{m_{u}}:g_{U}(u)\geq 0\}$.

We are now in a position to use Theorem 4 and deduce $FR(X_{0},f,T,Y,\gamma)\subseteq\{x\in\mathbb{R}^{n}:V(x,T)\leq 1+\gamma\}$.

Proposition 5 shows that an outer approximation of the reachable set of an ODE can be constructed from the convex optimization problem found in (10). Furthermore we have argued using the objective function $-\log\det\{P(T)\}$ results in an heuristic optimal representation of the forward reachable set under the volume metric.

5.1 Reachable Sets of ODE’s With No Inputs

We can consider the simpler case of an ODE with no inputs,

[TABLE]

where $f:\mathbb{R}^{n}\to\mathbb{R}^{n}$ is the vector field and $x_{0}\in\mathbb{R}^{n}$ is the initial condition.

We note that the definitions of solution map and forward reachable set can be slightly altered and easily applied to ODE’s with no input of Form (13). In this case for a vector field $f:\mathbb{R}^{n}\to\mathbb{R}^{n}$, a set $X_{0}\subset\mathbb{R}^{n}$, and $T>0$ we denote

[TABLE]

where $\phi_{f}:\mathbb{R}^{n}\times\mathbb{R}\to\mathbb{R}^{n}$ is the solution map of (13).

Following a similar argument we used to derive the Optimization Problem (10) we now propose a convex SOS optimization problem for outer set approximation of forward reachable sets of ODE’s of form (13). The optimization problem is denoted by $S_{2}(d,T,f,g_{X},g_{C},h)$:

[TABLE]

where

[TABLE]

Corollary 6

For some $X_{0}\subset\mathbb{R}^{n}$, $f:\mathbb{R}^{n}\times\mathbb{R}^{m_{u}}\times\mathbb{R}^{m_{w}}\to\mathbb{R}^{n}$, and $T>0$, suppose there exists functions $g_{X}:\mathbb{R}^{n}\to\mathbb{R}$ and $g_{C}:\mathbb{R}^{n}\to\mathbb{R}$ such that $X_{0}\subseteq\{x\in\mathbb{R}^{n}:g_{X}(x)\geq 0\}$ and $FR(X_{0},f,t)\subseteq\{x\in\mathbb{R}^{n}:g_{C}(x)\geq 0\}$ for all $t\in[0,T]$. Then if $P(t):\mathbb{R}\to\mathbb{R}$ solves optimization problem $S_{2}(d,T,f,g_{X},g_{C},h)$, found in (14), for some $d\in\mathbb{N}$ and $h(t)=t(T-t)$, then $FR(X_{0},f,T)\subseteq\{x\in\mathbb{R}^{n}:z_{d}(x)P(T)z_{d}(x)\leq 1\}$.

6 Numerical Examples

In this section we will now compute several forward reachable sets for different dynamical systems. Here constraints from Optimization Problems (10) and (14) were enforced using software such as SOSTOOLS, found in Prajna et al. (2002), that reformulates the problem as an SDP. Using efficient primal-dual interior point methods for SDP’s we are able to solve such proposed problems, see Prajna et al. (1994).

6.1 Computation Of Reachable Sets Of Systems With No Inputs

Example 1

Let us consider the linear ODE:

[TABLE]

where $A=\begin{bmatrix}0&-1\\ 1&0\end{bmatrix}$. Since the eigenvalues of $A$ are $\pm i$ it follows (15) produces non-stable circular trajectories. We now solve the optimization problem $S_{2}(d,T,f,g_{X},g_{C},h)$, found in (14), for the ODE (15), where $d=3$, $T=2$, $f(x)=Ax$, $g_{X}(x)=0.25^{2}-(x_{1}-0.5)^{2}-x_{2}^{2}$, $g_{C}(x)=50^{2}-x_{1}^{2}-x_{2}^{2}$, and $h(t)=t[T-t]$. The results are displayed in Figure 1. Here $24^{2}$ terminal trajectory points, shown in blue, were approximately found by forward-time integrating (15) starting from initial points, shown in red. We find that $FR(X_{0},f,T)\subseteq\{x\in\mathbb{R}^{n}:z_{d}(x)P(T)z_{d}(x)\leq 1\}$; this is demonstrated by the black sublevel set containing the blue circle of points in the figure.

Since the dynamics of the ODE (15) are simple it was numerically tractable to solve the Optimization Problem for 100 evenly spaced points between times [0,5]; Figure 2 shows the results.

Example 2

Let us now consider the Lorenz system defined by the three dimensional second order nonlinear ODE:

[TABLE]

We solved optimization problem $S_{2}(d,T,f,g_{X},g_{C},h)$, found in (14), for the ODE (16) with, $d=2$, $T=0.75$; $f(x)=[\sigma(x_{2}-x_{1}),x_{1}(\rho-x_{3})-x_{2},x_{1}x_{2}-\beta x_{3}]^{T}$ $g_{X}(x)=1-x_{1}^{2}-x_{2}^{2}-x_{3}^{2}$; $g_{C}(x)=50^{2}-x_{1}^{2}-x_{2}^{2}-x_{3}^{2}$; $h(t)=t[T-t]$; and system parameters $\sigma=10$, $\beta=8/3$, $\rho=28$. Using the solution to the optimization problem, $P(t)$, we then constructed the function $V(x,t)=z_{d}(x)^{T}P(t)z_{d}(x)$. In Figure 6 we have plotted our outer approximation of $FR(X_{0},f,T)$, the set $\{x\in\mathbb{R}^{2}:V(x,T)\leq 1\}$ shown as the gray 3D boundary. As expected initial points contained inside the set $X_{0}$, shown as red points, transition to terminal points, shown as blue points, contained in $\{x\in\mathbb{R}^{2}:V(x,T)\leq 1\}$. Moreover the boundary of the set $\{x\in\mathbb{R}^{2}:V(x,T)\leq 1\}$ has a similar shape to the Lorenz attractor.

6.2 Computation Of Reachable Sets Of Systems With Inputs

Example 3

We next consider a third order nonlinear system with bounded L2 inputs, from Jarvis-Wloszek et al. (2005) and Yin et al. (2018), given in the following ODE

[TABLE]

where $w\in W=\{w\in L_{2}^{1}[0,T]:\int_{0}^{T}w(t)^{T}w(t)dt<\gamma\}$.

We solved Optimization Problem $S_{1}(d,T,f,\gamma,g_{X},g_{C},g_{U},h)$, found in (10), for this ODE and the following terms; $d=2$; $T=1$; $f(x)=[-x_{1}(t)+x_{2}(t)-x_{1}x_{2}^{2},-x_{2}(t)-x_{1}^{2}x_{2}+w]^{T}$; $\gamma=2$; $g_{X}(x)=1^{2}-x_{1}^{2}-x_{2}^{2}$; $g_{C}(x)=1.75^{2}-x_{1}^{2}-x_{2}^{2}$; $g_{U}(x)=0$; and $h(t)=t[T-t]$. Terms involving $u$ were ignored from the optimization problem as there is no point-wise input in (17). Using the solution of the optimization problem, $P(t)$, we then constructed the function $V(x,t)=z_{d}(x)^{T}P(t)z_{d}(x)$. In Figure 4 we have plotted the sublevel set $\{x\in\mathbb{R}^{2}:V(x,T)\leq 1+\gamma\}$, shown as the black curve. We have also plotted initial conditions, shown as red points, contained in the set $X_{0}=\{x\in\mathbb{R}^{n}:g_{X}(x)\geq 0\}$. The solution maps at time $T$ generated for randomly generated polynomial inputs of the form $w(t)=c^{T}z_{d}(t)$, shown as the blue points, are also plotted. As expected, since we have shown the set $\{x\in\mathbb{R}^{2}:V(x,T)\leq 1+\gamma\}$ is an outer approximation of $FR(X_{0},f,T,\emptyset,\gamma)$, the blue points are all contained inside the black line.

Example 4

Let us now consider the Van der Pol oscillator with both bounded $L_{2}$ and point-wise input disturbances defined by the nonlinear ODE:

[TABLE]

where $w\in W=\{w\in L_{2}^{m_{w}}[0,T]:\int_{0}^{T}w(t)^{T}w(t)dt<\gamma\}$ and $\mu\in U=\{u\in L_{2}^{1}[0,T]:u(t)\in[\underline{u},\bar{u}]\text{ for all }t\in[0,T]\}$ is a modeling parameter that measures damping strength.

We solved optimization problem $S_{1}(d,T,f,\gamma,g_{X},g_{C},g_{U},h)$, found in (10), for the ODE (18) with, $d=3$; $T=1$; $f(x)=[x_{2}+w,-x_{1}+\mu x_{2}(1-x_{1}^{2})]^{T}$; $\gamma=0.25$; $g_{X}(x)=1-x_{1}^{2}-x_{2}^{2}$; $g_{C}(x)=8^{2}-x_{1}^{2}-x_{2}^{2}$; $g_{U}(u)=(u-\underline{u})(\bar{u}-u)$; $\underline{u}=0.5$; $\bar{u}=1.5$; and $h(t)=t[T-t]$. Using the solution of the optimization problem, $P(t)$, we then constructed the function $V(x,t)=z_{d}(x)^{T}P(t)z_{d}(x)$. To compare our approximation of the forward reachable set with no input disturbances we then also solved $S_{2}(d,T,f,g_{X},g_{C},h)$, found in (14), for the ODE (18) with $w(t)=0$ and $\mu(t)=1$; that is $f(x)=[x_{2},-x_{1}+x_{2}(1-x_{1}^{2})]^{T}$ is now used. Using the solution of this optimization problem, $\tilde{P}(t)$, we then constructed the function $\tilde{V}(x,t)=z_{d}(x)^{T}\tilde{P}(t)z_{d}(x)$. In Figure 5 we have plotted the outer approximation of the forward reachable set for ODE with input disturbances, the set $\{x\in\mathbb{R}^{2}:{V}(x,T)\leq 1+\gamma\}$ shown as the dotted black curve, and outer approximation of the forward reachable set for ODE with no input disturbances, the set $\{x\in\mathbb{R}^{2}:\tilde{V}(x,T)\leq 1\}$ shown as the black curve. As expected our approximation of the forward reachable set for the ODE with input disturbances is much larger than without.

Moreover in Figure 6 we have again solved $S_{2}(d,T,f,g_{X},g_{C},h)$, found in (14), for the ODE (18) with no input disturbances and $\mu(t)=1$ for a later terminal time of $T=5$. Interestingly the boundary of the approximated forward reachable set is very similar to the Van der Pol limit cycle; shown as the red line which was approximately found by forward time integrating (18) at a starting position close to the limit cycle.

7 Conclusion

We have illustrated a method for finding approximations of forward reachable sets by sublevel sets of an SOS polynomials that solve a convex optimization problem. We have used an objective function based on the determinant to heuristically minimizes the volume of these sublevel sets and improve our outer approximations. We have applied our methods to finding reachable of nonlinear systems with both $L_{2}$ or point wise bounded input disturbances. Outer approximations for the reachable sets for the Lorenz system and Van der Pol system show a similar shape to the attractor set and limit cycle respectively.