Semidefinite Approximations of Reachable Sets for Discrete-time Polynomial Systems

TL;DR

This paper introduces a method to compute certified outer approximations of the reachable set for discrete-time polynomial systems using a hierarchy of semidefinite relaxations, ensuring convergence to the true set.

Contribution

It develops a novel hierarchy of semidefinite relaxations for approximating reachable sets, with convergence guarantees and practical numerical implementation.

Findings

Hierarchies of super level sets approximate reachable sets with increasing accuracy.

Convex semidefinite programs compute polynomial coefficients efficiently.

Numerical examples demonstrate the method's effectiveness.

Abstract

We consider the problem of approximating the reachable set of a discrete-time polynomial system from a semialgebraic set of initial conditions under general semialgebraic set constraints. Assuming inclusion in a given simple set like a box or an ellipsoid, we provide a method to compute certified outer approximations of the reachable set. The proposed method consists of building a hierarchy of relaxations for an infinite-dimensional moment problem. Under certain assumptions, the optimal value of this problem is the volume of the reachable set and the optimum solution is the restriction of the Lebesgue measure on this set. Then, one can outer approximate the reachable set as closely as desired with a hierarchy of super level sets of increasing degree polynomials. For each fixed degree, finding the coefficients of the polynomial boils down to computing the optimal solution of a convex semidefinite program. When the degree of the polynomial approximation tends to infinity, we provide strong convergence guarantees of the super level sets to the reachable set. We also present some application examples together with numerical results.

Equations

$${\mathbf{X}}^0:=\{\mathbf{x}\in {\mathbb{R}}^n:g_1^0(\mathbf{x})\ge 0,\dots ,g_{m^0}^0(\mathbf{x})\ge 0\}$$

$${\mathbf{X}}^T:={\mathbf{X}}^0\cup f({\mathbf{X}}^0)\cup f(f({\mathbf{X}}^0))\cup \cdots \cup f^T({\mathbf{X}}^0),$$

$${\mathbf{X}}^{\infty }:=\underset{T\to \infty }{\lim }{\mathbf{X}}^T$$

$$\mathbf{X}:=\{\mathbf{x}\in {\mathbb{R}}^n:g_1(\mathbf{x})\ge 0,\dots ,g_m(\mathbf{x})\ge 0\}$$

$$y_{\beta }^{\mathbf{A}}:=\int {\mathbf{x}}^{\beta }{\lambda }_{\mathbf{A}}(d\mathbf{x})\in \mathbb{R},\beta \in {\mathbb{N}}^n$$

$$\mu \le \nu$$

$$f_{\#}\mu (\mathbf{A}):=\mu (f^{-1}(\mathbf{A}))=\mu (\{\mathbf{x}\in \mathbf{X}:f(\mathbf{x})\in \mathbf{A}\})$$

$${\mathbf{Q}}_r^0$$

$${\mathbf{Q}}_r$$

$$\mathbf{x}\mapsto p(\mathbf{x})=\sum _{\beta \in {\mathbb{N}}^n}{p}_{\beta }{\mathbf{x}}^{\beta },$$

$$({\mathbf{M}}_r(\mathbf{y}){)}_{\beta ,\gamma }:={\ell }_{\mathbf{y}}({\mathbf{x}}^{\beta +\gamma }),\forall \beta ,\gamma \in {\mathbb{N}}_r^n.$$

$$({\mathbf{M}}_r(q\mathbf{y}){)}_{\beta ,\gamma }:={\ell }_{\mathbf{y}}(q(\mathbf{x}){\mathbf{x}}^{\beta +\gamma }),\forall \beta ,\gamma \in {\mathbb{N}}_r^n.$$

$${\mu }_{t+1}$$

$$\nu$$

$${\mu }_T+\nu =f_{\#}\nu +{\mu }_0.$$

$${\mathbf{Y}}^0:={\mathbf{X}}^0,{\mathbf{Y}}^t:=f^t({\mathbf{X}}^0)\backslash {\mathbf{X}}^{t-1},t=1,\dots ,T.$$

$${\mathbf{X}}^T={\cup }_{t=0}^T{\mathbf{Y}}^t.$$

$${\lambda }_{{\mathbf{X}}^T}+{\nu }^T=f_{\#}{\nu }^T+{\mu }_0^T.$$

$${\lambda }_{{\mathbf{X}}^{\infty }}+\nu =f_{\#}\nu +{\mu }_0.$$

$${\lambda }_{{\mathbf{Y}}^t}+{\nu }_t=f_{\#}{\nu }_t+{\mu }_{0,t}.$$

$${\nu }^T:=\sum _{t=0}^T{\nu }_t\in {\mathcal{M}}_{+}(\mathbf{X})$$

$${\mu }_0^T:=\sum _{t=0}^T{\mu }_{0,t}\in {\mathcal{M}}_{+}({\mathbf{X}}^0)$$

$${\lambda }_{{\mathbf{X}}^T}+{\nu }^T=f_{\#}{\nu }^T+{\mu }_0^T.$$

$$p^T:=\underset{{\mu }_0,\mu ,\hat{\mu },\nu ,a}{\sup }$$

$${\int }_{\mathbf{X}}\nu +a=T\mathrm{vol}\mathbf{X},$$

$$\mu +\nu =f_{\#}\nu +{\mu }_0,$$

$$\mu +\hat{\mu }={\lambda }_{\mathbf{X}},$$

$${\mu }_0\in {\mathcal{M}}_{+}({\mathbf{X}}^0),\mu ,\hat{\mu },\nu \in {\mathcal{M}}_{+}(\mathbf{X}),a\in {\mathbb{R}}_{+}.$$

$${\mu }^{\prime }({\mathbf{X}}^{\infty })+{\mu }^{\prime }(\mathbf{X}\backslash {\mathbf{X}}^{\infty })={\mu }^{\prime }(\mathbf{X})={\mu }_0^{\prime }({\mathbf{X}}^0)={\mu }_0^{\prime }({\mathbf{X}}^{\infty })={\mu }^{\prime }({\mathbf{X}}^{\infty })+{\nu }^{\prime }({\mathbf{X}}^{\infty })-f_{\#}{\nu }^{\prime }({\mathbf{X}}^{\infty }),$$

$${\int }_{\mathbf{X}}\nu (d\mathbf{x})+a$$

$${\int }_{\mathbf{X}}v(\mathbf{x})\mu (d\mathbf{x})+{\int }_{\mathbf{X}}v(\mathbf{x})\nu (d\mathbf{x})$$

$${\int }_{\mathbf{X}}w(\mathbf{x})\mu (d\mathbf{x})+{\int }_{\mathbf{X}}w(\mathbf{x})\hat{\mu }(d\mathbf{x})$$

$$p^T=\underset{x}{\sup }$$

$$\mathcal{A}x=b,$$

$$x\in E_1^{+},$$

Full text

The proposed method consists of building a hierarchy of relaxations for an infinite-dimensional moment problem. Under certain assumptions, the optimal value of this problem is the volume of the reachable set and the optimum solution is the restriction of the Lebesgue measure on this set. Then, one can outer approximate the reachable set as closely as desired with a hierarchy of super level sets of increasing degree polynomials. For each fixed degree, finding the coefficients of the polynomial boils down to computing the optimal solution of a convex semidefinite program. When the degree of the polynomial approximation tends to infinity, we provide strong convergence guarantees of the super level sets to the reachable set. We also present some application examples together with numerical results.

Keywords:

reachable set, discrete-time polynomial systems, polynomial optimization, semidefinite programming, moment relaxations, sums of squares, convex optimization.

1 Introduction

Given a dynamical polynomial system described by a discrete-time (difference) equation, the (forward) reachable set (RS) is the set of all states that can be reached from a set of initial conditions under general state constraints. This set appears in different fields such as optimal control, hybrid systems or program analysis. In general, computing or even approximating the RS is a challenge. Note that the RS is typically non-convex and non-connected, even in the case when the set of initial conditions is convex and the dynamics are linear.

Computing or approximating RS has been a topic of intensive research in the last four decades. When the dynamics of the discrete-time system is linear, one can rely on contractive algorithms based on finite LP relaxations combined with polyhedral projections [8]. For more details and historical surveys, we refer the interested reader to [9, 10] as well as [12] for an extension to hybrid systems. A recent approach [15] has extended the scope of problems for which one can construct polyhedral bounds on the reachable set. Again, the method relies on (parametric) LP and numerical integration procedures. The works [5] compares approaches combining LP relaxations with Bernstein decompositions or Krivine/Handelman representations of nonnegative polynomials (also based on sums of squares certificates). See also [6] using template polyhedra and Bernstein form of polynomials. In recent work [13], the author relies on Bernstein expansions to over approximate the RS after in discrete-time after a finite number of time steps. The same technique allows to perform parameter synthesis.

However, all these methods based on linear relaxations often fail to construct tight bounds since they build convex over approximations of possibly nonconvex sets. Furthermore, they usually do not provide convergence guarantees.

Another classical approach relies on Lyapunov theory (see e.g. [39, § 5.7]) in order to approximate from outside. This can be done in a continuous-time setting (with possible extension to discrete-time systems), i.e., when the state variable is constrained from an initial condition to satisfy an ordinary differential equation $\dot{\mathbf{x}}=f(\mathbf{x},t)$. The idea is to search for a Lyapunov function $v$ (also called value or Bellman function in the context of optimal control) which is negative on the set of initial conditions and with negative derivative of states satisfying some general constraints. These inequalities provide sufficient conditions for the RS to be included in the sublevel set of $v$. In the case where the set of initial (resp. general) state constraints are defined by polynomial inequalities, the difficulty of computing such a funcion $v$ can be practically addressed. This is done while reducing the search space to polynomials of bounded degree and by replacing the inequalities satisfied by $v$ (and its total derivative) by stronger equality constraints involving $v$ and (weighted) sums of squares (SOS) of polynomials. Since the weights are the polynomials defining the set of initial and general constraints, computing $v$ together with these SOS polynomials boils down to solving a semidefinite program of fixed size. This general framework has been used in [36] for the safety verification of hybrid systems. In this case, the function $v$ is called a “barrier certificate” and can be constructed by computing an SOS decomposition. The zero level set of $v$ separates a given unsafe region from all possible trajectories starting from a prescribed set of initial conditions.

These dual Lyapunov certificates relying on SOS decompositions also allow to obtain approximations of the (backward) reachable set (also called region of attraction) [11]. In [42], the authors proved the existence of a Lyapunov function, whose sublevel set is the region of attraction of a given equilibrium point of a continuous-time system. When the degree of the approximation $v$ is fixed in advance, one can obtain convergence guarantees by increasing the degree of the SOS polynomials. However, one has no guarantee that when the degree of $v$ goes to infinity, the approximation conservatism asymptotically vanishes. In addition, the conservatism of such approximations relying on dual Lyapunov certificates is not easy to estimate in a systematic way.

In this paper, we propose a characterization of the RS as the solution of an infinite-dimensional linear programming (LP) problem. This characterization is done by considering a hierarchy of converging convex programs through moment relaxations of the LP. Doing so, one can compute tight outer approximations of the RS. Such outer approximations yield invariants for the discrete-time system, which are sets where systems trajectories are confined.

This general methodology is deeply inspired from previous research efforts. The idea of formulation relying on LP optimization over probability measures appears in [24], with a hierarchy of semidefinite programs (SDP) also called moment-sum-of-squares or sometimes Lasserre hierarchy, whose optimal values converge from below to the infimum of a multivariate polynomial. One can see outer approximations of sets as the analogue of lower approximations of real-valued functions. In [19], the authors leverage on these techniques to address the problem of computing outer approximations by single polynomial super level sets of basic compact semialgebraic sets described by the intersection of a finite number of given polynomial super level sets. Further work focused on approximating semialgebraic sets for which such a description is not explicitly known or difficult to compute: in [27], the author derives converging outer approximations of sets defined with existential quantifiers; in [33], the authors approximate the image of a compact semialgebraic set under a polynomial map. The current study can be seen as an extension of [33] where instead of considering only one iteration of the map, we consider infinitely many iterations starting from a set of initial conditions.

This methodology has also been successfully applied for several problems arising in the context of polynomial systems control. Similar convergent hierarchies appear in [17], where the authors approximate the region of attraction (ROA) of a controlled polynomial system subject to compact semialgebraic constraints in continuous time. This framework is extended to hybrid systems in [38]. Note that the ROA is not a semialgebraic set in general. The authors of [21] build upon the infinite-dimensional LP formulation of the ROA problem while providing a similar framework to characterize the maximum controlled invariant (MCI) for discrete and continuous time polynomial dynamical systems. The framework used for ROA and MCI computation both rely on occupation measures. These allow to measure the time spent by solutions of differential or difference equations. As solutions of a linear transport equation called the Liouville Equation, occupation measures also capture the evolution of the semialgebraic set describing the initial conditions. As mentioned in [17], the problem of characterizing the (forward) RS in a continuous setting and finite horizon could be done as ROA computation by using a time-reversal argument.

The modeling power of this approach also extends to the analysis of attractors of dynamical systems, e.g., by approximating the moments and support of invariant measures in both continuous and discrete-time settings [22, 32].

In the present study, we handle the problem in a discrete setting and infinite horizon. Our contribution follows a similar approach but requires to describe the solution set of another Liouville Equation.

In contrast with previous work, our contributions are the following:

• •

we rely on a infinite-dimensional LP formulation to handle the general discrete-time RS problem under semialgebraic state and initial conditions;

• •

we build a hierarchy of finite-dimensional SDP relaxations for this infinite-dimensional LP problem. Under additional assumptions, the optimal value of this LP is the volume of the RS, whose optimum is the restriction of the Lebesgue measure on the RS;

• •

we use the solutions of these SDP relaxations to approximate the RS as closely as desired, with a sequence of certified outer approximations defined as (possibly nonconvex) polynomial super level sets, which is less restrictive than linear or convex approximations used in the literature (such as polyhedra or ellipsoids).

In the sequel, we focus on computation of semidefinite approximations of the forward RS for discrete-time polynomial systems. The problem statement is formalized in Section 2 and reformulated in Section 3 into a primal optimization problem over probability measures satisfying Liouville’s Equation. We explain how to obtain the dual problem in Section 3.3. Then, we show in Section 4 how to solve in practice the primal problem with moment relaxations, as well as the dual with sums of squares strenghtenings. In both cases, this boils down to solving a hierarchy of finite-dimensional SDP problems. We illustrate the method with several numerical experiments in Section 5.

2 Problem Statement and Prerequisite

2.1 Forward Reachable set for Discrete-time Polynomial Systems

Given $r,n\in\mathbb{N}$, let $\mathbb{R}[\mathbf{x}]$ (resp. $\mathbb{R}_{2r}[\mathbf{x}]$) stands for the vector space of real-valued $n$-variate polynomials (resp. of degree at most $2r$) in the variable $\mathbf{x}=(x_{1},\ldots,x_{n})\in\mathbb{R}^{n}$. We are interested in the polynomial discrete-time system defined by

• •

a set of initial constraints assumed to be compact basic semi-algebraic:

[TABLE]

defined by given polynomials $g_{1}^{0},\ldots,g_{m^{0}}^{0}\in\mathbb{R}[\mathbf{x}]$, $m^{0}\in\mathbb{N}$;

• •

a polynomial transition map $f:\mathbb{R}^{n}\to\mathbb{R}^{n}$, $\mathbf{x}\mapsto f(\mathbf{x}):=(f_{1}(\mathbf{x}),\dots,f_{n}(\mathbf{x}))\in\mathbb{R}^{n}[\mathbf{x}]$ of degree $d:=\max\{\deg f_{1},\ldots,\deg f_{n}\}$.

Given $T\in\mathbb{N}$, let us define the set of all admissible trajectories after at most $T$ iterations of the polynomial transition map $f$, starting from any initial condition in $\mathbf{X}^{0}$:

[TABLE]

with $f^{T}$ denoting the $T$-fold composition of $f$. Then, we consider the reachable set (RS) of all admissible trajectories:

[TABLE]

and we make the following assumption in the sequel:

Assumption 2.1**.**

The RS $\mathbf{X}^{\infty}$ is included in a given compact basic semi-algebraic set

[TABLE]

defined by polynomials $g_{1},\ldots,g_{m}\in\mathbb{R}[\mathbf{x}]$, $m\in\mathbb{N}$.

Example 1**.**

Let us consider $\mathbf{X}^{0}=[1/2,1]$ and $f(x)=x/4$. Then $\mathbf{X}^{\infty}=[1/2,1]\cup[1/8,1/4]\cup[1/32,1/16]\ldots$ is included in the basic compact semialgebraic set $\mathbf{X}=[0,1]$, so Assumption 2.1 holds. Note that $\mathbf{X}^{\infty}$ is not connected within $\mathbf{X}$.

We denote the closure of $\mathbf{X}^{\infty}$ by $\bar{\mathbf{X}}^{\infty}$ . Obviously $\mathbf{X}^{\infty}\subseteq\bar{\mathbf{X}}^{\infty}$ and the inclusion can be strict. To circumvent this difficulty later on in Section 3 and Section 4, we make the following assumption in the remainder of the paper.

Assumption 2.2**.**

The volume of the RS is equal to the volume of its closure, i.e. $\operatorname{vol}\mathbf{X}^{\infty}=\operatorname{vol}\bar{\mathbf{X}}^{\infty}$.

Example 2**.**

Let $\mathbf{X}^{0}=[1/2,1]$ and $f(x)=x/2$. Then $\mathbf{X}^{\infty}=[1/2,1]\cup[1/4,1/2]\cup[1/8,1/4]\ldots=(0,1]$ is a half-closed interval within $\mathbf{X}=[0,1]$. Note that $\bar{\mathbf{X}}^{\infty}=\mathbf{X}$, so that $\operatorname{vol}\mathbf{X}^{\infty}=\operatorname{vol}\bar{\mathbf{X}}^{\infty}=1$ and Assumption 2.2 is satisfied.

2.2 Prerequisite and Working Assumptions

Given a compact set $\mathbf{A}\subset\mathbb{R}^{n}$, we denote by $\mathcal{M}(\mathbf{A})$ the vector space of finite signed Borel measures supported on $\mathbf{A}$, namely real-valued functions from the Borel sigma algebra $\mathcal{B}(\mathbf{A})$. The support of a measure $\mu\in\mathcal{M}(\mathbf{A})$ is defined as the closure of the set of all points $\mathbf{x}$ such that $\mu(\mathbf{B})\neq 0$ for any open neighborhood $\mathbf{B}$ of $\mathbf{x}$. We note $\mathcal{C}(\mathbf{A})$ the Banach space of continuous functions on $\mathbf{A}$ equipped with the sup-norm. Let $\mathcal{C}(\mathbf{A})^{\prime}$ stand for the topological dual of $\mathcal{C}(\mathbf{A})$ (equipped with the sup-norm), i.e. the set of continuous linear functionals of $\mathcal{C}(\mathbf{A})$. By a Riesz identification theorem (see for instance [28]), $\mathcal{C}(\mathbf{A})^{\prime}$ is isomorphically identified with $\mathcal{M}(\mathbf{A})$ equipped with the total variation norm denoted by $\|\cdot\|_{\text{TV}}$. Let $\mathcal{C}_{+}(\mathbf{A})$ (resp. $\mathcal{M}_{+}(\mathbf{A})$) stand for the cone of non-negative elements of $\mathcal{C}(\mathbf{A})$ (resp. $\mathcal{M}(\mathbf{A})$). The topology in $\mathcal{C}_{+}(\mathbf{A})$ is the strong topology of uniform convergence in contrast with the weak-star topology in $\mathcal{M}_{+}(\mathbf{A})$. See [37, Section 21.7] and [4, Chapter IV] or [29, Section 5.10] for functional analysis, measure theory and applications in convex optimization.

With $\mathbf{X}$ a basic compact semialgebraic set as in (2), the restriction of the Lebesgue measure on a subset $\mathbf{A}\subseteq\mathbf{X}$ is $\lambda_{\mathbf{A}}(d\mathbf{x}):=\mathbf{1}_{\mathbf{A}}(\mathbf{x})\,d\mathbf{x}$, where $\mathbf{1}_{\mathbf{A}}:\mathbf{X}\to\{0,1\}$ stands for the indicator function of $\mathbf{A}$, namely $\mathbf{1}_{\mathbf{A}}(\mathbf{x})=1$ if $\mathbf{x}\in\mathbf{A}$ and $\mathbf{1}_{\mathbf{A}}(\mathbf{x})=0$ otherwise.

The moments of the Lebesgue measure on $\mathbf{A}$ are denoted by

[TABLE]

where we use the multinomial notation $\mathbf{x}^{\beta}:=x^{\beta_{1}}_{1}x^{\beta_{2}}_{2}\ldots x^{\beta_{n}}_{n}$. The Lebesgue volume of $\mathbf{A}$ is $\operatorname{vol}\mathbf{A}:=y^{\mathbf{A}}_{0}=\int\lambda_{\mathbf{A}}(d\mathbf{x})$.

Given $\mu,\nu\in\mathcal{M}(\mathbf{A})$, the notation

[TABLE]

stands for $\nu-\mu\in\mathcal{M}_{+}(\mathbf{A})$, and we say that $\mu$ is dominated by $\nu$.

Given $\mu\in\mathcal{M}_{+}(\mathbf{X})$, the so-called pushforward measure or image measure, see e.g. [1, Section 1.5], of $\mu$ under $f$ is defined as follows:

[TABLE]

for every set $\mathbf{A}\in\mathcal{B}(\mathbf{X})$. The main property of the pushforward measure is the change-of-variable formula: $\int_{\mathbf{A}}v(\mathbf{x})f_{\#}\mu(d\mathbf{x})=\int_{f^{-1}(\mathbf{A})}v(f(\mathbf{x}))\mu(d\mathbf{x})$, for all $v\in\mathcal{C}(\mathbf{A})$.

With $\mathbf{X}^{0}$ a basic compact semialgebraic set as in (1), we set $r_{j}^{0}:=\lceil(\deg g_{j}^{0})/2\rceil,j=1,\dots,m^{0}$ and with $\mathbf{X}$ a basic compact semialgebraic set as in (2), we set $r_{j}:=\lceil(\deg g_{j})/2\rceil,j=1,\dots,m$. Let $\Sigma[\mathbf{x}]$ stand for the cone of polynomial sums of squares (SOS) and let $\Sigma_{r}[\mathbf{x}]$ denote the cone of SOS polynomials of degree at most $2r$, namely $\Sigma_{r}[\mathbf{x}]:=\Sigma[\mathbf{x}]\cap\mathbb{R}_{2r}[\mathbf{x}]$.

For the ease of further notation, we set $g_{0}^{0}(\mathbf{x}):=1$ and $g_{0}(\mathbf{x}):=1$. For each integer $r$, let $\mathbf{Q}_{r}^{0}$ (resp. $\mathbf{Q}_{r}$) be the $r$-truncated quadratic module generated by $g_{0}^{0},\dots,g_{m}^{m^{0}}$ (resp. $g_{0},\dots,g_{m}$):

[TABLE]

To guarantee the convergence behavior of the relaxations presented in the sequel, we need to ensure that polynomials which are positive on $\mathbf{X}^{0}$ (resp. $\mathbf{X}$) lie in $\mathbf{Q}_{r}^{0}$ (resp. $\mathbf{Q}_{r}$) for some $r\in\mathbb{N}$. The existence of such SOS-based representations is guaranteed by Putinar’s Positivstellensaz (see e.g. [25, Section 2.5]), when the following condition holds:

Assumption 2.3**.**

There exists a large enough integer $N^{0}$ (resp. $N$) such that one of the polynomials describing the set $\mathbf{X}^{0}$ (resp. $\mathbf{X}$) is equal to $g_{0}^{i}:=N^{0}-\|\mathbf{x}\|_{2}^{2}$ (resp. $g^{i}:=N-\|\mathbf{x}\|_{2}^{2}$).

This assumption is slightly stronger than compactness. Indeed, compactness of $\mathbf{X}^{0}$ (resp. $\mathbf{X}$) already ensures that each variable has finite lower and upper bounds. One (easy) way to ensure that Assumption 2.3 holds is to add a redundant constraint involving a well-chosen $N^{0}$ (resp. $N$) depending on these bounds, in the definition of $\mathbf{X}^{0}$ (resp. $\mathbf{X}$).

From now on, the over approximation set $\mathbf{X}$ of the set $\mathbf{X}^{\infty}$ is assumed to be “simple” (e.g. a ball or a box), meaning that $\mathbf{X}$ fulfills the following condition:

Assumption 2.4**.**

The moments (3) of the Lebesgue measure on $\mathbf{X}$ are available analytically.

Remark 1**.**

Since we are interested in characterizing the RS of polynomial systems with bounded trajectories, Assumption 2.1 and Assumption 2.4 are not restrictive. As mentioned above, Assumption 2.3 can be ensured by using Assumption 2.1. While relying on Assumption 2.2, we restrict ourselves to discrete-time systems where the boundary of the RS has zero Lebesgue volume.

Let $\mathbb{N}_{0}$ stands for the set of positive integers. For all $r\in\mathbb{N}$, we set $\mathbb{N}^{n}_{r}:=\{\beta\in\mathbb{N}^{n}:\sum_{j=1}^{n}\beta_{j}\leq r\}$, whose cardinality is $\binom{n+r}{r}$. Then a polynomial $p\in\mathbb{R}[\mathbf{x}]$ is written as follows:

[TABLE]

and $p$ is identified with its vector of coefficients $\mathbf{p}=(p_{\beta})$ in the canonical basis $(\mathbf{x}^{\beta})$, $\beta\in\mathbb{N}^{n}$.

Given a real sequence $\mathbf{y}=(y_{\beta})_{\beta\in\mathbb{N}^{n}}$, let us define the linear functional $\ell_{\mathbf{y}}:\mathbb{R}[\mathbf{x}]\to\mathbb{R}$ by $\ell_{\mathbf{y}}(p):=\sum_{\beta}p_{\beta}y_{\beta}$, for every polynomial $p$.

Then, we associate to $\mathbf{y}$ the so-called moment matrix $\mathbf{M}_{r}(\mathbf{y})$, that is the real symmetric matrix with rows and columns indexed by $\mathbb{N}_{r}^{n}$ and the following entrywise definition:

[TABLE]

Given a polynomial $q\in\mathbb{R}[\mathbf{x}]$, we also associate to $\mathbf{y}$ the so-called localizing matrix, that is the real symmetric matrix $\mathbf{M}_{r}(q\,\mathbf{y})$ with rows and columns indexed by $\mathbb{N}_{r}^{n}$ and the following entrywise definition:

[TABLE]

In the sequel, we propose a method to compute an outer approximation of $\mathbf{X}^{\infty}$ as the super level set of a polynomial. Given its degree $2r$, such a polynomial may be produced using standard numerical optimization tools. Under certain assumptions, the resulting outer approximation can be as tight as desired and converges (with respect to the $L_{1}$ norm on $\mathbf{X}$) to the set $\mathbf{X}^{\infty}$ as $r$ tends to infinity.

3 Primal-Dual Infinite-dimensional LP

3.1 Forward Reachable Sets and Liouville’s Equation

For a given $T\in\mathbb{N}_{0}$ and a measure $\mu_{0}\in\mathcal{M}_{+}(\mathbf{X}^{0})$, let us define the measures $\mu_{1},\dots,\mu_{T},\mu\in\mathcal{M}_{+}(\mathbf{X})$ as follows:

[TABLE]

The measure $\nu$ is a (discrete-time) occupation measure: if $\mu_{0}=\delta_{x_{0}}$ is the Dirac measure at $x_{0}\in\mathbf{X}^{0}$ then $\mu_{t}=\delta_{x_{t}}$ and $\nu=\delta_{x_{0}}+\delta_{x_{1}}+\cdots+\delta_{x_{T-1}}$ , i.e. $\nu$ measures the time spent by the state trajectory in any subset of $\mathbf{X}$ after $T$ iterations, if initialized at $x_{0}$.

Lemma 3.1**.**

For any $T\in\mathbb{N}_{0}$ and $\mu_{0}\in\mathcal{M}(\mathbf{X}^{0})$, there exist $\mu_{T},\nu\in\mathcal{M}(\mathbf{X})$ solving the discrete Liouville Equation:

[TABLE]

Proof.

Follows readily from the definitions given in (4). ∎

Now let

[TABLE]

Note that the RS defined in Section 2.1 is equal to

[TABLE]

Further results involve statements relying on the following technical assumption:

Assumption 3.2**.**

$\lim_{T\to\infty}\sum_{t=0}^{T}t\operatorname{vol}\mathbf{Y}^{t}<\infty$.

This assumption seems to be strong or unjustified at that point. Moreover, if we do not know if the assumption is satisfied a priori, there is an a posteriori validation based on duality theory. Thanks to this validation procedure, we could check that the assumption was satisfied in all the examples we processed. Moreover, we were not able to find a discrete-time polynomial system violating this assumption. We will explain this later on in Section 4, with Theorem 4.1 (see also Remark 4). Before that, we prove in Lemma 3.4 that equation (5) holds when $\mu_{T}=\lambda_{\mathbf{X}^{T}}$, the restriction of the Lebesgue measure over the RS. We rely on the auxiliary result from [33, Lemma 4.1]:

Lemma 3.3**.**

Let $\mathbf{S},\mathbf{B}\subseteq\mathbf{X}$ be such that $f(\mathbf{S})\subseteq\mathbf{B}$. Given a measure $\mu_{1}\in\mathcal{M}_{+}(\mathbf{B})$, there is a measure $\mu_{0}\in\mathcal{M}_{+}(\mathbf{S})$ such that $f_{\#}\mu_{0}=\mu_{1}$ if and only if there is no continuous function $v\in\mathcal{C}(\mathbf{B})$ such that $v(f(\mathbf{x}))\geq 0$ for all $\mathbf{x}\in\mathbf{S}$ and $\int_{\mathbf{B}}v(\mathbf{y})d\mu_{1}(\mathbf{y})<0$.

Lemma 3.4**.**

For any $T\in\mathbb{N}_{0}$, there exist $\mu_{0}^{T}\in\mathcal{M}_{+}(\mathbf{X}^{0})$ and $\nu^{T}\in\mathcal{M}_{+}(\mathbf{X})$ such that the restriction of the Lebesgue measure over $\mathbf{X}^{T}$ solves the discrete Liouville Equation:

[TABLE]

In addition, if Assumption 3.2 holds, then there exist $\mu_{0}\in\mathcal{M}_{+}(\mathbf{X}^{0})$ and $\nu\in\mathcal{M}_{+}(\mathbf{X})$ such that the restriction of the Lebesgue measure over $\mathbf{X}^{\infty}$ solves the discrete Liouville Equation:

[TABLE]

Proof.

We first show that for each $t=1,\dots,T$ there exist a measure $\nu_{t}\in\mathcal{M}_{+}(\mathbf{X})$ and $\mu_{0,t}\in\mathcal{M}_{+}(\mathbf{X}^{0})$ such that

[TABLE]

It follows from Lemma 3.3 that there exists a measure $\mu_{0,t}\in\mathcal{M}_{+}(f^{-t}(\mathbf{Y}^{t}))$ such that $\lambda_{\mathbf{Y}^{t}}=f^{t}_{\#}\mu_{0,t}$ (with the notations $f\leftarrow f^{t}$, $\mathbf{S}\leftarrow f^{-t}(\mathbf{Y}^{t})$, $\mathbf{B}\leftarrow\mathbf{X}$, $\mu_{1}\leftarrow\lambda_{\mathbf{Y}^{t}}$ and $\mu_{0}\leftarrow\mu_{0,t}$). Indeed, by the change-of-variable formula recalled in Section 2.2, one has $\int_{\mathbf{X}}v(\mathbf{x})d\lambda_{\mathbf{Y}^{t}}(\mathbf{x})=\int_{\mathbf{Y}^{t}}v(\mathbf{x})d\lambda_{\mathbf{Y}^{t}}(\mathbf{x})=\int_{f^{-t}(\mathbf{Y}^{t})}v(f^{t}(\mathbf{x}))d\mu_{0,t}(\mathbf{x})$. Thus, it is impossible to find a function $v\in\mathcal{C}(\mathbf{X})$ such that $v(f^{t}(\mathbf{x}))\geq 0$ for all $\mathbf{x}\in f^{-t}(\mathbf{Y}^{t})$ while satisfying $\int_{\mathbf{X}}v(\mathbf{x})d\lambda_{\mathbf{Y}^{t}}(\mathbf{x})<0$. Since $f^{-t}(\mathbf{Y}^{t})\subseteq\mathbf{X}^{0}$, it implies that the measures $\mu_{0,t}$ and $\nu_{t}:=\sum_{i=0}^{t-1}f^{i}_{\#}\mu_{0,t}$ satisfy $\mu_{0,t}\in\mathcal{M}_{+}(\mathbf{X}^{0})$, $\nu_{t}\in\mathcal{M}_{+}(\mathbf{X})$ and Equation (8).

Letting $\mu_{0,0}:=\lambda_{\mathbf{X}^{0}}$, we now prove the first claim of the Lemma by showing that the measures

[TABLE]

and

[TABLE]

satisfy

[TABLE]

Since $\mathbf{X}^{T}=\bigcup_{t=0}^{T}\mathbf{Y}^{t}$ and the $\mathbf{Y}^{t}$ are disjoint, one can write $\lambda_{\mathbf{X}^{T}}=\sum_{t=0}^{T}\lambda_{\mathbf{Y}^{t}}$. This decomposition together with Equation (8) allows to show the first claim.

Now, we prove the second claim of the Lemma by showing that the respective limits of the measure sequences $(\lambda_{\mathbf{X}^{T}})_{T}$, $(\nu^{T})_{T}$ and $(\mu_{0}^{T})_{T}$ exist as $T\to\infty$, and that the respective limits $\lambda_{\mathbf{X}^{\infty}}$, $\nu$ and $\mu_{0}$ satisfy Equation (7):

• •

The limit of $\lambda_{\mathbf{X}^{T}}$ exists and is equal to $\lambda_{\mathbf{X}^{\infty}}$ by definition of $\mathbf{X}^{\infty}$;

• •

The limit of $\mu_{0}^{T}$ exists since $\|\mu_{0}^{T}\|_{\text{TV}}=\sum_{t=0}^{T}\int_{\mathbf{X}^{0}}\mu_{0,t}=\sum_{t=0}^{T}\int_{\mathbf{Y}^{t}}\lambda_{\mathbf{Y}^{t}}=\sum_{t=0}^{T}\operatorname{vol}\mathbf{Y}^{t}\leq\operatorname{vol}\mathbf{X}^{\infty}\leq\operatorname{vol}\mathbf{X}<\infty$. Therefore there is a subsequence which converges to a certain $\mu_{0}\in\mathcal{M}_{+}(\mathbf{X}^{0})$ for the weak-star topology.

• •

The limit of $\nu^{T}$ exists since $\|\nu^{T}\|_{\text{TV}}=\sum_{t=0}^{T}\int_{\mathbf{X}}\nu_{t}=\sum_{t=0}^{T}\int_{\mathbf{Y}_{t}}\sum_{i=0}^{t-1}f^{i}_{\#}\mu_{0,t}=\sum_{t=0}^{T}t\operatorname{vol}\mathbf{Y}^{t}<\infty$ by Assumption 3.2. Therefore, there is a subsequence which converges to a certain $\nu\in\mathcal{M}_{+}(\mathbf{X})$ for the weak-star topology.

Finally, taking the limit to infinity of both sides of Equation (9) yields the initial claim. ∎

Remark 2**.**

In Lemma 3.4, the measure $\mu_{0}^{T}$ (resp. $\mu_{0}$) can be thought as distribution of mass for the initial states of trajectories reaching $\mathbf{X}^{T}$ (resp. $\mathbf{X}^{\infty}$) but it has a total mass which is not required to be normalized to one.

The mass of $\nu^{T}$ measures the volume averaged w.r.t. $\mu_{0}$ occupied by state trajectories reaching $\mathbf{X}^{T}$ after $T$ iterations, by contrast with the mass of $\lambda_{\mathbf{X}^{T}}$ which measures the volume of $\mathbf{X}^{T}$.

The mass of $\nu$ measures the volume averaged w.r.t. $\mu_{0}$ occupied by state trajectories reaching the RS $\mathbf{X}^{\infty}$, by contrast with the mass of $\lambda_{\mathbf{X}^{\infty}}$ which measures the exact RS volume.

3.2 Primal Formulation

To approximate the set $\mathbf{X}^{\infty}$, one considers the infinite-dimensional linear programming (LP) problem, for any $T\in\mathbb{N}_{0}$:

[TABLE]

The first equality constraint ensures that the mass of the occupation measure $\nu$ is bounded (by $T\operatorname{vol}\mathbf{X}$). The second one ensures that Liouville’s equation is satisfied by the measures $\mu_{0}$, $\nu$ and $\mu$, as in Lemma 3.4. The last one ensures that $\mu$ is dominated by the restriction of the Lebesgue measure on $\mathbf{X}$ implying that the mass of $\mu$ (and thus the optimal value $p^{T}$) is bounded by $\operatorname{vol}\mathbf{X}$. The next result explains how the solution of LP (10) relates to $\lambda_{\mathbf{X}^{\infty}}$, the restriction of the Lebesgue measure to the RS.

Lemma 3.5**.**

For any $T\in\mathbb{N}_{0}$, LP (10) admits an optimal solution $(\mu_{0}^{*},\mu^{*},\hat{\mu}^{*},\nu^{*},a^{*})$ such that $\mu^{*}=\lambda_{\mathbf{S}^{T}}$ for some set $\mathbf{S}^{T}$ satisfying $\mathbf{X}^{T}\subseteq\mathbf{S}^{T}\subseteq\bar{\mathbf{X}}^{\infty}$ and $\operatorname{vol}\mathbf{S}^{T}=p^{T}$.

In addition if Assumption 3.2 holds then there exists $T_{0}\in\mathbb{N}$ such that for all $T\geq T_{0}$ one has $\mathbf{S}^{T}=\bar{\mathbf{X}}^{\infty}$, LP (10) has a unique optimal solution with $\mu^{*}=\lambda_{\mathbf{X}^{\infty}}$ and $p^{T}=\operatorname{vol}\mathbf{X}^{\infty}$.

Proof.

Let $T\in\mathbb{N}_{0}$. First we show that the feasible set of LP (10) is nonempty and compact, and that LP (10) has at least one optimal solution $(\mu_{0}^{*},\mu^{*},\hat{\mu}^{*},\nu^{*},a^{*})$. The feasible set of LP (10) is nonempty as $(\mu_{0},\mu,\hat{\mu},\nu,a)=(0,0,\lambda_{\mathbf{X}},0,T\operatorname{vol}\mathbf{X})$ is a feasible solution. Let us consider the sequences of measures $(\mu_{0n})_{n}$, $(\mu_{n})_{n}$, $(\hat{\mu}_{n})_{n}$, $(\nu_{n})_{n}$ and the nonnegative real sequence $(a_{n})_{n}$ such that each $(\mu_{0n},\mu_{n},\hat{\mu}_{n},\nu_{n},a_{n})_{n}$ is feasible for LP (10). One has:

• •

$\|\mu_{n}\|_{\text{TV}}+\|\hat{\mu}_{n}\|_{\text{TV}}=\operatorname{vol}\mathbf{X}<\infty$ (as $\mathbf{X}$ is bounded), which implies that $\|\mu_{n}\|_{\text{TV}}$ and $\|\hat{\mu}_{n}\|_{\text{TV}}$ are both bounded;

• •

$\|\mu_{0n}\|_{\text{TV}}=\|\mu_{n}\|_{\text{TV}}<\infty$, which implies that $\|\mu_{0n}\|_{\text{TV}}$ is bounded;

• •

$\|\nu_{n}\|_{\text{TV}}$ and $a_{n}$ are both bounded by $T\operatorname{vol}\mathbf{X}$.

Thus, the feasible set of the LP (10) is bounded for the weak-star topology. Now we show that this set is closed. Assume that the sequences of measures respectively converge weakly-star to $\mu_{0}$, $\mu$, $\hat{\mu}$, $\nu$ and the sequence of nonnegative real numbers converge to $a$. For each $A\in\mathcal{B}(\mathbf{X})$, one has $\nu(f^{-1}A)\to\nu_{n}(f^{-1}(A))$ and hence $f_{\#}\nu_{n}\to f_{\#}\nu$. Thus, $(\mu_{0},\mu,\hat{\mu},\nu,a)$ is also feasible for LP (10), proving that the feasible set of the LP (10) is closed in the metric including the weak-star topology. The existence of an optimal solution $(\mu_{0}^{*},\mu^{*},\hat{\mu}^{*},\nu^{*},a^{*})$ follows from the fact that LP (10) has a linear objective function with a weak-star compact feasible set.

By Lemma 3.4, there exist $\mu_{0}^{T}\in\mathcal{M}_{+}(\mathbf{X}^{0})$ and $\nu^{T}\in\mathcal{M}_{+}(\mathbf{X})$ such that the restriction of the Lebesgue measure over $\mathbf{X}^{T}$ solves the discrete Liouville Equation, i.e. $\lambda_{\mathbf{X}^{T}}+\nu^{T}=f_{\#}\nu^{T}+\mu_{0}^{T}$ and $\int_{\mathbf{X}}\nu^{T}=\sum_{t=0}^{T}t\operatorname{vol}\mathbf{Y}^{t}\leq T\operatorname{vol}\mathbf{X}$. With $a^{T}=T\operatorname{vol}\mathbf{X}-\sum_{t=0}^{T}t\operatorname{vol}\mathbf{Y}^{t}$ and $\hat{\mu}^{T}=\lambda_{\mathbf{X}}-\lambda_{\mathbf{X}^{T}}$, it follows that $(\mu_{0}^{T},\mu^{T},\hat{\mu}^{T},\nu^{T},a^{T})$ is feasible for LP (10).

Next, given any feasible solution $(\mu_{0}^{\prime},\mu^{\prime},\hat{\mu}^{\prime},\nu^{\prime},a^{\prime})$ of LP (10), we show that the support of $\mu^{\prime}$ is included in $\bar{\mathbf{X}}^{\infty}$. Using the Liouville Equation and the fact that $\mu_{0}^{\prime}$ is supported on $\mathbf{X}^{0}\subseteq\mathbf{X}^{\infty}$, one has:

[TABLE]

which proves that $\mu^{\prime}(\mathbf{X}\backslash\mathbf{X}^{\infty})=\nu^{\prime}(\mathbf{X}^{\infty})-f_{\#}\nu^{\prime}(\mathbf{X}^{\infty})$. Since $\mathbf{X}^{\infty}\subseteq f^{-1}(\mathbf{X}^{\infty})$, one has $\nu^{\prime}(\mathbf{X}^{\infty})\leq\nu^{\prime}(f^{-1}(\mathbf{X}^{\infty}))=f_{\#}\nu^{\prime}(\mathbf{X}^{\infty})$, implying that $\mu^{\prime}(\mathbf{X}\backslash\mathbf{X}^{\infty})=0$. Thus, the support of $\mu^{\prime}$ is included in $\bar{\mathbf{X}}^{\infty}$.

Now we assume that $\sum_{t=0}^{\infty}t\operatorname{vol}\mathbf{Y}^{t}>T\operatorname{vol}\mathbf{X}$. Hence, there exist a minimal integer $T_{2}\in\mathbb{N}_{0}$ and measures $\mu_{0}^{T_{2}}\in\mathcal{M}_{+}(\mathbf{X}^{0})$, $\nu^{T_{2}}\in\mathcal{M}_{+}(\mathbf{X})$ such that $T<T_{2}$, $\lambda_{\mathbf{X}^{T_{2}}}+\nu^{T_{2}}=f_{\#}\nu^{T_{2}}+\mu_{0}^{T_{2}}$ and $T\operatorname{vol}\mathbf{X}<\int_{\mathbf{X}}\nu^{T_{2}}$. Therefore, for some $T_{1}\in\mathbb{N}_{0}$ with $T_{1}\geq T$, there exist some sets $\mathbf{S}_{0},\dots,\mathbf{S}_{T_{1}}\subseteq\mathbf{X}$ such that $\mathbf{S}_{0}=\mathbf{Y}^{0},\dots,\mathbf{S}_{T}=\mathbf{Y}^{T},\mathbf{S}_{T+1}\subseteq\mathbf{Y}^{T+1},\dots,\mathbf{S}_{T_{1}}\subseteq\mathbf{Y}_{T_{1}}$ and $\sum_{t=0}^{T_{1}}t\operatorname{vol}\mathbf{S}_{t}=T\operatorname{vol}\mathbf{X}$. By defining $\mathbf{S}^{T}:=\bigcup_{t=0}^{T_{1}}\mathbf{S}_{t}$, one has $\mathbf{X}^{T}\subseteq\mathbf{S}^{T}\subset\mathbf{X}^{T_{2}}\subseteq\mathbf{X}^{\infty}$.

As for Lemma 3.4, one proves that for all $t=0,\dots,T_{1}$, there exist $\mu^{\prime}_{0,t}\in\mathcal{M}_{+}(\mathbf{X}^{0})$, $\nu^{\prime}_{t}\in\mathcal{M}_{+}(\mathbf{X})$ such that $\lambda_{\mathbf{S}_{t}}+\nu^{\prime}_{t}=f_{\#}\nu^{\prime}_{t}+\mu^{\prime}_{0,t}$.

By superposition, the measures $\lambda_{\mathbf{S}^{T}}$, $\mu_{0}^{T_{1}}:=\sum_{i=0}^{T_{1}}\mu^{\prime}_{0,t}$ and $\nu^{T_{1}}:=\sum_{i=0}^{T_{1}}\nu^{\prime}_{t}$ satisfy $\lambda_{\mathbf{S}^{T}}+\nu^{T_{1}}=f_{\#}\nu^{T_{1}}+\mu_{0}^{T_{1}}$ and $\int_{\mathbf{X}}\nu^{T_{1}}=\sum_{t=0}^{T_{1}}t\operatorname{vol}\mathbf{S}_{t}=T\operatorname{vol}\mathbf{X}$.

This implies that $\lambda_{\mathbf{S}^{T}}$ is feasible for LP (10). Then one proves, exactly as in the proof of [19, Theorem 3.1], that the quintuple $(\mu_{0}^{T_{1}},\lambda_{\mathbf{S}^{T}},\lambda_{\mathbf{X}\backslash\mathbf{S}^{T}},\nu^{T_{1}},0)$ is optimum, yielding the optimal value $p^{T}=\operatorname{vol}\mathbf{S}^{T}$.

Under Assumption 3.2, Lemma 3.4 implies that there exist $\mu_{0}\in\mathcal{M}_{+}(\mathbf{X}^{0})$ and $\nu\in\mathcal{M}_{+}(\mathbf{X})$ such that the restriction of the Lebesgue measure over $\mathbf{X}^{\infty}$ solves the discrete Liouville’s Equation (7). In addition, there exists $T_{0}\in\mathbb{N}$ such that for all $T\geq T_{0}$, $\sum_{t=0}^{T}t\operatorname{vol}\mathbf{Y}^{t}\leq T\operatorname{vol}\mathbf{X}$, thus $\int_{\mathbf{X}}\nu\leq T\operatorname{vol}\mathbf{X}$. This implies that $\lambda_{\mathbf{X}^{\infty}}$ is feasible for LP (10). With $\mu_{0}$ and $\nu$ as in Lemma 3.4, we define $\mu_{0}^{*}:=\mu_{0}$, $\mu^{*}:=\lambda_{\mathbf{X}^{\infty}}$, $\hat{\mu}^{*}:=\lambda_{\mathbf{X}}-\lambda_{\mathbf{X}^{\infty}}$, $\nu^{*}:=\nu$ and $a^{*}:=T\operatorname{vol}\mathbf{X}-\int_{\mathbf{X}}\nu$ and show as in the proof of [19, Theorem 3.1], that $(\mu_{0}^{*},\mu^{*},\hat{\mu}^{*},\nu^{*},a^{*})$ is optimum with unique $\mu^{*}$. This yields the optimal value $p^{T}=\int_{\mathbf{X}}\lambda_{\mathbf{X}^{\infty}}=\operatorname{vol}\bar{\mathbf{X}}^{\infty}=\operatorname{vol}\mathbf{X}^{\infty}$, where the last equality comes from Assumption 2.2. ∎

From now on, we refer to $\mathbf{S}^{T}$ as the support of the optimal solution $\mu^{*}$ of LP (10) which satisfies the condition of Lemma 3.5, i.e. $\mathbf{X}^{T}\subseteq\mathbf{S}^{T}\subseteq\bar{\mathbf{X}}^{\infty}$.

In the sequel, we formulate LP (10) as an infinite-dimensional conic problem on appropriate vector spaces. By construction, a feasible solution of problem (10) satisfies:

[TABLE]

for all continuous test functions $v,w\in\mathcal{C}(\mathbf{X})$.

Then, we cast problem (10) as a particular instance of a primal LP in the canonical form given in [4, 7.1.1]:

[TABLE]

with

• •

the vector space $E_{1}:=\mathcal{M}(\mathbf{X}^{0})\times\mathcal{M}(\mathbf{X})^{3}\times\mathbb{R}$ with its cone $E^{+}_{1}$ of non-negative elements;

• •

the vector space $F_{1}:=\mathcal{C}(\mathbf{X}^{0})\times\mathcal{C}(\mathbf{X})^{3}\times\mathbb{R}$;

• •

the duality $\langle\cdot,\cdot\rangle_{1}:E_{1}\times F_{1}\to\mathbb{R}$, given by the integration of continuous functions against Borel measures, since $E_{1}=F_{1}^{\prime}$ is the dual of $F_{1}$;

• •

the decision variable $x:=(\mu_{0},\mu,\hat{\mu},\nu,a)\in E_{1}$ and the reward $c:=(0,1,0,0,0)\in F_{1}$;

• •

$E_{2}:=\mathbb{R}\times\mathcal{M}(\mathbf{X})^{2}$, $F_{2}:=\mathbb{R}\times\mathcal{C}(\mathbf{X})^{2}$ and the right hand side vector $b:=(T\operatorname{vol}\mathbf{X},0,\lambda)\in E_{2}=F_{2}^{\prime}$;

• •

the linear operator $\mathcal{A}:E_{1}\to E_{2}$ given by

[TABLE]

Note that both spaces $E_{1}$, $E_{2}$ (resp. $F_{1}$, $F_{2}$) are equipped with the weak topologies $\sigma(E_{1},F_{1})$, $\sigma(E_{2},F_{2})$ (resp. $\sigma(F_{1},E_{1})$, $\sigma(F_{2},E_{2})$) and $\sigma(E_{1},F_{1})$ is the weak-star topology (since $E_{1}=F_{1}^{\prime}$). Observe that $\mathcal{A}$ is continuous with respect to the weak topology, as $\mathcal{A}^{\prime}(F_{2})\subset F_{1}$.

3.3 Dual Formulation

Using the same notations, the dual of the primal LP (14) in the canonical form given in [4, 7.1.2] reads:

[TABLE]

with

• •

the dual variable $y:=(u,v,w)\in F_{2}$;

• •

the (pre)-dual cone $F^{+}_{1}:=\mathbb{R}^{+}\times\mathcal{C}_{+}(\mathbf{X})^{2}$, whose dual is $E_{1}^{+}$;

• •

the duality pairing $\langle\cdot,\cdot\rangle_{2}:E_{2}\times F_{2}\to\mathbb{R}$, with $E_{2}=F_{2}^{\prime}$;

• •

the adjoint linear operator $\mathcal{A}^{\prime}:F_{2}\to F_{1}$ given by

[TABLE]

Using our original notations, the dual LP of problem (10) then reads:

[TABLE]

Theorem 3.6**.**

For a fixed $T\in\mathbb{N}_{0}$, there is no duality gap between primal LP (10) and dual LP (16), i.e. $p^{T}=d^{T}$ and there exists a minimizing sequence $(u_{k},v_{k},w_{k})_{k\in\mathbb{N}}$ for the dual LP (16).

In addition, if $u_{k}=0$ for some $k\in\mathbb{N}$, then Assumption 3.2 holds and $p^{T}=d^{T}=\operatorname{vol}\mathbf{X}^{\infty}$.

Proof.

As in [21, Theorem 3], the zero duality gap follows from infinite-dimensional LP duality theory (for more details, see [4]). The feasible set of the LP (10) is nonempty in the metric inducing the weak-star topology on $\mathcal{M}(\mathbf{X}^{0})\times\mathcal{M}(\mathbf{X})^{3}\times\mathbb{R}$ since $(0,0,\lambda_{\mathbf{X}},0,T\operatorname{vol}\mathbf{X})$ is a trivial feasible solution. As shown in Lemma 3.5, the feasible set of the LP (10) is bounded for the same metric. Hence, the first claim follows from nonemptiness and boundedness of the feasible set of the primal LP (10).

Now, let us assume that there exists a minimizing sequence $(u_{k},v_{k},w_{k})_{k\in\mathbb{N}}$ of the dual LP (16) such that the dual variable $u_{k}$ is equal to 0 for some $k\in\mathbb{N}$. This implies that the corresponding constraint in the primal LP (10) is not saturated, i.e. $\int_{\mathbf{X}}\nu<T\operatorname{vol}\mathbf{X}$, for any solution $(\mu_{0},\mu,\hat{\mu},\nu,a)$. Let us show that Assumption 3.2 holds. Otherwise, by contradiction one would have $\sum_{t=0}^{\infty}t\operatorname{vol}\mathbf{Y}^{t}>T\operatorname{vol}\mathbf{X}$. In the proof of Lemma 3.5, we proved that this strict inequality implies the existence of a set $\mathbf{S}^{T}$ and an optimal solution $(\mu_{0}^{T_{1}},\lambda_{\mathbf{S}^{T}},\lambda_{\mathbf{X}\backslash\mathbf{S}^{T}},\nu^{T_{1}},0)$ for the primal LP (10), such that $\int_{\mathbf{X}}\nu^{T_{1}}=T\operatorname{vol}\mathbf{X}$, yielding a contradiction. Eventually, Lemma 3.5 implies that $p^{T}=d^{T}=\operatorname{vol}\mathbf{X}^{\infty}$. ∎

Remark 3**.**

When $u=0$, the first and third inequalities satisfied by $-v$ in dual LP (16) can be seen as a discrete-time analogue of the conditions satisfied by the barrier certificate in [36], that is $-v\leq 0$ on $\mathbf{X}^{0}$ and $-v\circ f\leq-v$ on $\mathbf{X}$, the latter one being similar to the barrier condition $\nabla(-v)\cdot f\leq-v$.

4 Primal-Dual Hierarchies of SDP Approximations

4.1 Primal-Dual Finite-dimensional SDP

For each $r\geq r_{\min}:=\max\{r^{0}_{1},\ldots,r^{0}_{m^{0}},r_{1},\ldots,r_{m}\}$, let $\mathbf{y}_{0}=(y_{0\beta})_{\beta\in\mathbb{N}_{2r}^{n}}$ be the finite sequence of moments up to degree $2r$ of the measure $\mu_{0}$. Similarly, let $\mathbf{y}$, $\hat{\mathbf{y}}$ and $\mathbf{z}$ stand for the sequences of moments up to degree $2r$, respectively associated with $\mu$, $\hat{\mu}$ and $\nu$. The infinite primal LP (10) can be relaxed with the following semidefinite program:

[TABLE]

Consider also the following semidefinite program, which is a strengthening of the infinite dual LP (16) and also the dual of Problem (17):

[TABLE]

where $\mathbf{Q}_{r}^{0}$, $\mathbf{Q}_{r}$ (resp. $\mathbf{Q}_{rd}$) are the $r$-truncated (resp. $rd$) quadratic module respectively generated by $g_{0}^{0},\dots,g_{m}^{m^{0}}$ and $g_{0},\dots,g_{m}$, as defined in Section 2.2.

Theorem 4.1**.**

Let $r\geq r_{\min}$. Suppose that the three sets $\mathbf{X}^{0}$, $\mathbf{S}^{T}$ and $\mathbf{X}\backslash\mathbf{S}^{T}$ have nonempty interior. Then:

$p^{T}_{r}=d^{T}_{r}$, i.e. there is no duality gap between the primal SDP program (17) and the dual SDP program (18). 2. 2.

The dual SDP program (18) has an optimal solution $(u_{r},v_{r},w_{r})\in\mathbb{R}\times\mathbb{R}_{2r}[\mathbf{x}]\times\mathbb{R}_{2r}[\mathbf{x}]$, and the sequence $(w_{r}+u_{r}T)$ converges to $\mathbf{1}_{\mathbf{S}^{T}}$ in $L_{1}$ norm on $\mathbf{X}$:

[TABLE] 3. 3.

Defining the sets

[TABLE]

it holds that

[TABLE] 4. 4.

In addition, if $u_{r}=0$ then Assumption 3.2 holds and the sequence $(w_{r})$ converges to $\mathbf{1}_{{\bar{\mathbf{X}}^{\infty}}}$ in $L_{1}$ norm on $\mathbf{X}$. Defining the sets

[TABLE]

its holds that

[TABLE]

and

[TABLE]

Proof.

Let $\mu^{*}:=\lambda_{\mathbf{S}^{T}}$, $\mu_{0}^{*}$ and $\nu^{*}$ be such that $\mu^{*}+\nu^{*}=f_{\#}\nu^{*}+\mu_{0}^{*}$ as in the proof of Lemma 3.5, and let $\hat{\mu}^{*}=\lambda_{\mathbf{X}}-\mu^{*}$ and $a^{*}:=T\operatorname{vol}\mathbf{X}-\int_{\mathbf{X}}\nu^{*}$ so that $(\mu_{0}^{*},\mu^{*},\hat{\mu}^{*},\nu^{*},a^{*})$ is feasible for LP (10). Given $r\geq r_{\min}$, let $\mathbf{y}_{0}$, $\mathbf{y}$, $\hat{\mathbf{y}}$ and $\mathbf{z}$ be the sequences of moments up to degree $2r$ of $\mu_{0}^{*}$, $\mu^{*}$, $\hat{\mu}^{*}$ and $\nu^{*}$ respectively. Then, as in the proof of the first item of Theorem 4.4 in [33], the optimal set of the primal SDP program (17) is nonempty and bounded and the result of [41] implies that there is no duality gap between the primal SDP program (17) and dual SDP program (18). 2. 2.

The proof follows by the same arguments as in the proof of the second item of Theorem 4.4 in [33]. One first shows that $(\mathbf{y}_{0},\mathbf{y},\hat{\mathbf{y}},\mathbf{z})$ is strictly feasible for program (17). It comes from the fact the three sets $\mathbf{X}^{0}$, $\mathbf{S}^{T}$ and $\mathbf{X}\backslash\mathbf{S}^{T}$ have nonempty interior, respectively yielding the positive definiteness of $\mathbf{M}_{r-r_{j}^{0}}(g_{j}\,\mathbf{y}_{0})$ for each $j=0,\dots,m^{0}$, $\mathbf{M}_{r-r_{j}}(g_{j}\,\mathbf{y})$, $\mathbf{M}_{r-r_{j}}(g_{j}\,\hat{\mathbf{y}})$ and $\mathbf{M}_{r-r_{j}}(g_{j}\,\mathbf{z})$ for each $j=0,\dots,m$.

Then, we conclude that the dual SDP program (18) has an optimal solution $(u_{r},v_{r},w_{r})\in\mathbb{R}\times\mathbb{R}_{2r}[\mathbf{x}]\times\mathbb{R}_{2r}[\mathbf{x}]$.

Now, one proves that there exists a sequence of polynomials $(w_{k})_{k\in\mathbb{N}}\subset\mathbb{R}[\mathbf{x}]$ such that $w_{k}(\mathbf{x})\geq\mathbf{1}_{{\mathbf{S}^{T}}}(\mathbf{x})$, for all $\mathbf{x}\in\mathbf{X}$ and such that

[TABLE]

Closedness of $\mathbf{S}^{T}$ implies that the indicator function $\mathbf{1}_{\mathbf{S}^{T}}$ is upper semi-continuous, so there exists a non-increasing sequence of bounded continuous functions $h_{k}:\mathbf{X}\to\mathbb{R}$ such that $h_{k}(\mathbf{x})\downarrow\mathbf{1}_{\mathbf{S}^{T}}(\mathbf{x})$, for all $\mathbf{x}\in\mathbf{X}$, as $k\to\infty$.

By the Monotone Convergence Theorem [3], $h_{k}\to\mathbf{1}_{\mathbf{S}^{T}}$ for the $L_{1}$ norm. Let $\epsilon>0$ be given. Using the Stone-Weierstrass Theorem [40], there exists a sequence of polynomials $(w_{k}^{\prime})_{k\in\mathbb{N}}\subset\mathbb{R}[\mathbf{x}]$, such that $\sup_{\mathbf{x}\in\mathbf{X}}\mid w_{k}^{\prime}(\mathbf{x})-h_{k}(\mathbf{x})\mid<\epsilon$, thus $\mid w_{k}^{\prime}(\mathbf{x})-\mathbf{1}_{\mathbf{S}^{T}}(\mathbf{x})\mid<2\epsilon$. With $u_{k}:=4\epsilon$, the polynomial $w_{k}:=w_{k}^{\prime}+2\epsilon$ satisfies $0<w_{k}-\mathbf{1}_{\mathbf{S}^{T}}<4\epsilon$ and $4T\epsilon<w_{k}+u_{k}T-\mathbf{1}_{\mathbf{S}^{T}}<(4+4T)\epsilon$, thus (20) holds.

Finally, let us define the polynomials $w_{k}^{\prime\prime}:=w_{k}+\epsilon$, $v_{k}:=w_{k}-1$ and prove that $(u_{k},v_{k}^{\prime\prime},w_{k}^{\prime\prime})$ is a feasible solution of (18) for large enough $r\in\mathbb{N}$. The inequalities $w_{k}^{\prime\prime}>w_{k}>\mathbf{1}_{\mathbf{S}^{T}}$ prove that $w_{k}^{\prime\prime}\in\mathbf{Q}_{r}(\mathbf{X})$, as a consequence of Putinar’s Positivstellensatz [25, Section 2.5]. Similarly, $w_{k}^{\prime\prime}-v_{k}^{\prime\prime}-1=\epsilon>0$ which proves that $w_{k}^{\prime\prime}-v_{k}^{\prime\prime}-1\in\mathbf{Q}_{r}(\mathbf{X})$. For each $\mathbf{x}\in\mathbf{X}$, one has $0<v^{\prime\prime}_{k}(\mathbf{x})<4\epsilon$. The left inequality proves that $v_{k}^{\prime\prime}\in\mathbf{Q}_{r}(\mathbf{X}^{0})$. Using both inequalities, one has for all $\mathbf{x}\in\mathbf{X}$, $u_{k}T+v^{\prime\prime}_{k}(f(\mathbf{x}))-v_{k}^{\prime\prime}(\mathbf{x})>0$, so $u_{k}T+v_{k}^{\prime\prime}\circ f-v_{k}^{\prime\prime}$ lies in $\mathbf{Q}_{rd}(\mathbf{X})$. 3. 3.

Let $\mathbf{x}\in\mathbf{X}^{T}$ and $(u_{r},v_{r},w_{r})\in\mathbb{R}_{2r}[\mathbf{y}]\times\mathbb{R}_{2r}[\mathbf{y}]$ be an optimal solution of (18). There exist $t\in\mathbb{N}$ such that $t\leq T$ and $\mathbf{x}\in\mathbf{Y}^{t}$. Thus, there exists $\mathbf{x}^{0}\in\mathbf{X}^{0}$ such that $\mathbf{x}=f^{t}(\mathbf{x}^{0})$. By feasibility, $v_{r}(\mathbf{x}^{0})\geq 0$ and by induction, one has $v_{r}(\mathbf{x})+u_{r}t\geq v_{r}(\mathbf{x}^{0})\geq 0$. Since $u_{r}\geq 0$, one has $v_{r}(\mathbf{x})+u_{r}T\geq v_{r}(\mathbf{x})+u_{r}t\geq 0$, yielding $\mathbf{X}^{T}_{r}\supseteq\mathbf{X}^{T}$. 4. 4.

If the dual variable $u_{r}$ satisfies $u_{r}=0$, then the corresponding constraint in the primal is not saturated, i.e. $z_{0}<Ty_{0}^{\mathbf{X}}$, which implies as in the proof of Lemma 3.6 that Assumption 3.2 holds. By Lemma 3.5, $\mathbf{S}^{T}=\bar{\mathbf{X}}^{\infty}$ and the sequence $(w_{r})$ converges to $\mathbf{1}_{{\bar{\mathbf{X}}^{\infty}}}$ in $L_{1}$ norm on $\mathbf{X}$. For all $\mathbf{x}\in\mathbf{X}^{\infty}$, there exists $\mathbf{x}^{0}\in\mathbf{X}^{0}$ and $t\in\mathbb{N}$ such that $\mathbf{x}=f^{t}(\mathbf{x}^{0})$. By feasibility, $v_{r}(\mathbf{x}^{0})\geq 0$ implies that $v_{r}(\mathbf{x})\geq v_{r}(\mathbf{x}^{0})\geq 0$. For all $\mathbf{x}\in\bar{\mathbf{X}}^{\infty}$, there exists a sequence $(\mathbf{x}_{n})\subset\mathbf{X}^{\infty}$ converging to $\mathbf{x}$, and as above, one shows that for all $n\in\mathbb{N}$, one has $v_{r}(\mathbf{x}_{n})\geq 0$. Continuity of $v_{r}$ implies that $v_{r}(\mathbf{x})\geq 0$. This proves that $\mathbf{X}^{\infty}_{r}\supseteq\bar{\mathbf{X}}^{\infty}\supseteq\mathbf{X}^{\infty}$. Finally, the proof of the convergence in volume is similar to the proof of Theorem 3.2 (6) in [33].

∎

Remark 4**.**

Theorem 4.1 states that one can over approximate the reachable states of the system after any arbitrary finite number of discrete-time steps (third item). In addition, Theorem 4.1 provides a sufficient condition to obtain a hierarchy of over approximations converging in volume to the RS (fourth item). If $u_{r}$ = 0, then the sequence of optimal values of SDP (18) is nonincreasing and converges to the volume of the RS. If one defines the piecewise polynomial $\overline{v}_{r}:=\min_{k\leq r}v_{k}$, then one shows as in [18, Theorem 1] that we obtain a nonincreasing sequence of functions converging to the indicator function of the RS: one has $\overline{v}_{r}\downarrow\mathbf{1}_{\overline{\mathbf{X}}^{\infty}}$ almost everywhere, almost uniformly and in Lebesgue measure.

4.2 Special Case: linear systems with ellipsoid constraints

Given $A\in\mathbb{R}^{n\times n}$, let us consider a discrete-time linear system $\mathbf{x}_{t+1}=\mathbf{A}\,\mathbf{x}_{t}$ with a set of initial constraints defined by the ellipsoid $\mathbf{X}^{0}:=\{\mathbf{x}\in\mathbb{R}^{n}:1\geq\mathbf{x}^{T}\,\mathbf{V}_{0}\,\mathbf{x}\,\}$ with $\mathbf{V}_{0}\in\mathbb{R}^{n\times n}$ a positive definite matrix.

Similarly the set of state constraints is defined by the ellipsoid $\mathbf{X}:=\{\mathbf{x}\in\mathbb{R}^{n}:1\geq\mathbf{x}^{T}\,\mathbf{G}\,\mathbf{x}\,\}$ with $\mathbf{G}\in\mathbb{R}^{n\times n}$ a positive definite matrix. Since one has $\mathbf{X}^{0}\subseteq\mathbf{X}$, it follows that $\mathbf{V}_{0}\succeq\mathbf{G}$.

Then, one can look for a quadratic function $v(\mathbf{x}):=1-\mathbf{x}^{T}\,\mathbf{V}\,\mathbf{x}$, with $\mathbf{V}\in\mathbb{R}^{n\times n}$ a positive definite matrix solution of the following SDP optimization problem:

[TABLE]

where $\mathbf{M}$ is the second-order moment matrix of the Lebesgue measure on $\mathbf{X}$, i.e. the matrix with entries

[TABLE]

Note that in this special case SDP (21) can be retrieved from SDP (18) and one can over approximate the reachable set with the superlevel set of $v$ or $w-1$:

Lemma 4.2**.**

SDP (21) is equivalent to SDP (18) with $r:=1$, $u_{r}:=0$, $v(\mathbf{x}):=1-\mathbf{x}^{T}\,\mathbf{V}\,\mathbf{x}$ and $w(\mathbf{x})=1+v(\mathbf{x})$. Thus, one has:

[TABLE]

Proof.

The polynomial $v$ is nonnegative over $\mathbf{X}^{0}$ if and only if $\mathbf{V}_{0}\succeq\mathbf{V}$. The “if” part comes from the fact that $\mathbf{V}_{0}\succeq\mathbf{V}$ implies that $1-\mathbf{x}^{T}\,\mathbf{V}\,\mathbf{x}\geq 1-\mathbf{x}^{T}\,\mathbf{V}_{0}\,\mathbf{x}\geq 0$, for all $\mathbf{x}\in\mathbf{X}^{0}$. The other implication is a consequence of the $S$-Lemma [43]: there exists a nonnegative constant $c$ such that $1-\mathbf{x}^{T}\,\mathbf{V}\,\mathbf{x}\geq c\,(1-\mathbf{x}^{T}\,\mathbf{V}_{0}\,\mathbf{x})$ for all $\mathbf{x}\in\mathbb{R}^{n}$. It yields $1\geq c$ when $\mathbf{x}=0$. By defining $\mathbf{y}=(1,\mathbf{x})$, one finds that $\mathbf{y}^{T}\,\begin{pmatrix}1-c&0\\ 0&\mathbf{V}_{0}-\mathbf{V}\end{pmatrix}\,\mathbf{y}\geq 0$ for all $\mathbf{y}\in\mathbb{R}^{n+1}$, which finally gives $\mathbf{V}_{0}\succeq\mathbf{V}$.

In addition, the polynomial $\mathbf{x}\mapsto v(\mathbf{A}\mathbf{x})-v(\mathbf{x})$ is nonnegative over $\mathbf{X}$ if and only if $\mathbf{V}\succeq\mathbf{A}^{T}\mathbf{V}\mathbf{A}$. The “if” part comes from the fact that $\mathbf{V}\succeq\mathbf{A}^{T}\mathbf{V}\mathbf{A}$ implies that $v(\mathbf{A}\mathbf{x})-v(\mathbf{x})=\mathbf{x}^{T}\,(\mathbf{V}-\mathbf{A}^{T}\mathbf{V}\mathbf{A})\,\mathbf{x}\geq 0$, for all $\mathbf{x}\in\mathbf{X}$. The other implication follows from the $S$-Lemma: there exists a constant $c\geq 0$ such that $v(\mathbf{A}\mathbf{x})-v(\mathbf{x})\geq c\,(1-\mathbf{x}^{T}\,\mathbf{G}\,\mathbf{x})$ or equivalently $\mathbf{x}^{T}\,(\mathbf{V}-\mathbf{A}^{t}\,\mathbf{V}\,\mathbf{A})\,\mathbf{x}\geq c\,(1-\mathbf{x}^{T}\,\mathbf{G}\,\mathbf{x})$, for all $\mathbf{x}\in\mathbb{R}^{n}$. As before, this yields $0\geq c$, thus $c=0$ and $\mathbf{V}-\mathbf{A}^{t}\,\mathbf{V}\,\mathbf{A}\succeq 0$, which finally gives $\mathbf{V}_{0}\succeq\mathbf{V}\succeq\mathbf{A}^{T}\mathbf{V}\mathbf{A}$.

Minimizing the integral (w.r.t. the Lebesgue measure) of $w$ over $\mathbf{X}$ is equivalent to maximizing the integral of the trace of the matrix $\mathbf{x}\,\mathbf{x}^{T}\mathbf{V}$ on $\mathbf{X}$.

∎

5 Numerical Experiments

Here, we present experimental benchmarks that illustrate our method. For a given positive integer $r$, we compute the polynomial solution $w_{r}$ of the dual SDP program (18). This dual SDP is modeled using the Yalmip toolbox [30] available within Matlab and interfaced with the SDP solver Mosek [2]. Performance results were obtained with an Intel Core i7-5600U CPU ($2.60\,$GHz) running under Debian 8.

For all experiments, we could find an optimal solution of the dual SDP program (18) either by adding the constraint $u=0$ or by setting $T=100$. In the latter case, the optimal solution is such that $u_{r}\simeq 0$ and the polynomial solution $w_{r}$ is the same than in the former case, up to small numerical errors (in practice the value of $u_{r}$ is less than $1\text{e--}5$). This implies that Assumption 3.2 is satisfied, i.e. the constraint of the mass of the occupation measure is not saturated, and yielding valid outer approximations of ${\mathbf{X}}^{\infty}$. The implementation is freely available on-line111www-verimag.imag.fr/~magron/reachsdp.tar.gz.

5.1 Toy Example

First, let us consider the discrete-time polynomial system defined by

[TABLE]

with initial state constraints $\mathbf{X}^{0}:=\{\mathbf{x}\in\mathbb{R}^{2}:(x_{1}-\frac{1}{2})^{2}+(x_{2}-\frac{1}{2})^{2}\leq 4^{-2}\}$ and general state constraints within the unit ball $\mathbf{X}:=\{\mathbf{x}\in\mathbb{R}^{2}:\|x\|_{2}^{2}\leq 1\}$. On Figure 1, we represent in light gray the outer approximations $\mathbf{X}^{\infty}_{r}$ of $\mathbf{X}^{\infty}$ obtained by our method, for increasing values of the relaxation order $r$ (from $2r=4$ to $14$). On each figure, the colored sets of points are obtained by simulation for the first $7$ iterates. More precisely, each colored set correspond to (under approximations of) the successive image sets $f(\mathbf{X}^{0}),\dots,f^{7}(\mathbf{X}^{0})$ of the points obtained by uniform sampling of $\mathbf{X}^{0}$ under $f,\dots,f^{7}$ respectively. The set $\mathbf{X}^{0}$ is blue and the set $f^{7}(\mathbf{X}^{0})$ is red, while intermediate sets take intermediate colors. The dotted circle represents the boundary of the unit ball $\mathbf{X}$. Figure 1 shows that the over approximations are already quite tight for low degrees.

5.2 Cathala System

Consider the Cathala System (see [21, Section 7.1.2]):

[TABLE]

with initial state constraints $\mathbf{X}^{0}:=\{\mathbf{x}\in\mathbb{R}^{2}:(x_{1}+0.6)^{2}+(x_{2}-0.5)^{2}\leq 0.4^{2}\}$ and state constraints $\mathbf{X}:=\{\mathbf{x}\in\mathbb{R}^{2}:\|x\|_{2}^{2}\leq 1.8^{2}\}$. The value $-0.5952$ corresponds to a parameter for which this system has an attractor (see Figure 2), the Cathala system being known to exhibit chaotic behavior [34].

5.3 FitzHugh-Nagumo Neuron Model

Consider the discretized version (taken from [6, Section 5]) of the FitzHugh-Nagumo model [14], which is originally a continuous-time polynomial system modelling the electrical activity of a neuron:

[TABLE]

with initial state constraints $\mathbf{X}^{0}:=[1,1.25]\times[2.25,2.5]$ and state constraints $\mathbf{X}:=\{\mathbf{x}\in\mathbb{R}^{2}:(\frac{x_{1}-0.1}{3.6})^{2}+(\frac{x_{2}-1.25}{1.75})^{2}\leq 1\}$. Figure 3 illustrates that the outer approximations provide useful indications on the system behavior, in particular for higher values of $r$. Indeed $\mathbf{X}^{\infty}_{5}$ and $\mathbf{X}^{\infty}_{6}$ capture the presence of the central “hole” made by periodic trajectories and $\mathbf{X}^{\infty}_{7}$ shows that there is a gap between the first discrete-time steps and the iterations corresponding to these periodic trajectories.

5.4 Julia Map

Consider the discrete-time system $z^{+}:=z^{2}+c$, with the state variable $z\in\mathbb{C}$ and parameter $c\in\mathbb{C}$. By setting $z=x_{1}+ix_{2}$ and $c=c_{1}+ic_{2}$, with $i$ the imaginary unit, we obtain the following equivalent quadratic two-dimensional formulation:

[TABLE]

with initial state constraints $\mathbf{X}^{0}:=\{\mathbf{x}\in\mathbb{R}^{2}:\|\mathbf{x}\|_{2}^{2}\leq 0.1^{2}\}$ and state constraints $\mathbf{X}\subseteq[-1.2,0.2]\times[-0.5,0.6]$. This recurrence allows to generate the filled Julia set, defined as the set of initial conditions for which the RS of the above quadratic system is bounded. Connectivity of the filled Julia set is ensured when $c$ belongs to the Mandelbrot set. In [21, Section 7.1.3], the authors provide over approximations of the sets of initial condition for different values of the parameter $c$, in particular for the case $c=-0.7+0.2i$ belonging to the Mandelbrot set and $c=-0.9+0.2i$ which lies outside the Mandelbrot set.

By contrast with the experimental results provided in [21, Section 7.1.3], Figure 4 depicts over approximations of the RS for different values of the parameter $c$. In particular, for cases where the parameter $c=-0.7\pm 0.2i$ lies inside the Mandelbrot set (corresponding to Figure 4(a) and Figure 4(b)), the over approximation $\mathbf{X}^{\infty}_{5}$ indicates that the trajectories possibly converge to an attractor. For cases where the parameter $c=-0.9\pm 0.2i$ lies outside the Mandelbrot set (corresponding to Figure 4(c) and Figure 4(d)), the over approximation $\mathbf{X}^{\infty}_{5}$ proves a disconnected behavior, possibly implying the presence of two attractors.

5.5 Phytoplankton Growth Model

Consider the discretized version of the Phytoplankton growth model (also taken from [6, Section 5]). This model is obtained after making assumptions, corroborated experimentally by biologists in order to represent such growth phenomena [7], yielding the following discrete-time polynomial system:

[TABLE]

with initial state constraints $\mathbf{X}^{0}:=[-0.3,-0.2]^{2}\times[-0.05,0.05]$ and state constraints $\mathbf{X}:=[-0.5,1.5]\times[-0.5,0.5]^{2}$. Figure 5 illustrates the system convergence behavior towards an equilibrium point for initial conditions near the origin. One way to obtain more accurate information on such systems would be to design a subdivision procedure (e.g., with branch-and-bound techniques), which boils down to zooming on specific areas of the RS.

6 Conclusion and Perspectives

This paper presented an infinite-dimensional primal-dual LP characterization of the (forward) reachable set (RS) for discrete-time polynomial systems with semialgebraic initial and general state constraints. The problem can be practically handled through solving a hierarchy of finite dimensional primal-dual SDP relaxations.

In particular, the hierarchy of dual SDP problems yields sequences of polynomials of increasing degrees, allowing to construct certified outer approximations of the RS while ensuring convergence guarantees (w.r.t. the $L_{1}$ norm) to the indicator function of the RS when the mass of some occupation measure is bounded. Our approach happens to be not only theoretically consistent but also practically efficient.

In some cases, it is possible to complement the hierarchy of convergent outer approximations of set of interest (ROA, MCI) by providing a sequence of inner approximations. For instance, the work [20] uses similar tools from measure theory to derive such inner approximations of the ROA. Future research perspectives include the study a complementary hierarchy of inner approximations for the RS, in the spirit of [20]. We also intend to investigate the RS problem for continuous time polynomial systems as for the infinite-dimensional convex modeling of the maximum controlled invariant, with infinite horizon [21]. In addition, it would be worth to apply the framework of [17], relying on occupation measures to approximate the region of attraction (ROA). A time-reversal argument would allow to formulate the RS problem in continuous time and finite horizon as ROA characterization. Finally, we also intend to develop a formally certified framework, inspired from [31], in order to guarantee the correctness of the over approximations.