Convex Risk Bounded Continuous-Time
Trajectory Planning in Uncertain Nonconvex Environments

∅∅∅Ashkan Jasour*, Weiqiao Han*, and Brian Williams

MIT, Computer Science and Artificial Intelligence Laboratory

This work was partially supported by the Boeing grant 6943358. *These authors contributed equally to the paper.

In this paper, we address the trajectory planning problem in uncertain nonconvex static and dynamic environments that contain obstacles with probabilistic location, size, and geometry. To address this problem, we provide a risk bounded trajectory planning method that looks for continuous-time trajectories with guaranteed bounded risk over the planning time horizon. Risk is defined as the probability of collision with uncertain obstacles.
Existing approaches to address risk bounded trajectory planning problems either are limited to Gaussian uncertainties and convex obstacles or rely on sampling-based methods that need uncertainty samples and time discretization. To address the risk bounded trajectory planning problem,
we leverage the notion of risk contours to transform the risk bounded planning problem into
a deterministic optimization problem. Risk contours are the set of all points in the uncertain environment with guaranteed bounded risk. The obtained deterministic optimization is, in general, nonlinear and nonconvex time-varying optimization. We provide convex methods based on sum-of-squares optimization to efficiently solve the obtained nonconvex time-varying optimization problem and obtain the continuous-time risk bounded trajectories without time discretization. The provided approach deals with arbitrary probabilistic uncertainties, nonconvex and nonlinear, static and dynamic obstacles, and is suitable for online trajectory planning problems.

§ INTRODUCTION

In order for robots to navigate safely in the real world, they need to plan safe trajectories to avoid static and moving obstacles, such as humans and vehicles, under perception uncertainties. The motion planning problem in dynamic environments is known to be computationally hard [1]. In this paper, we address the trajectory planning problem in uncertain nonconvex static and dynamic environments that contain obstacles with probabilistic location, size, and geometry. In this problem, the time-varying, nonconvex, and probabilistic nature of the obstacle-free safe regions makes the trajectory planning problem challenging.

Several approaches have been proposed to address the trajectory planning problems. In the absence of obstacles, one can use standard convex optimization to look for polynomial trajectories that satisfy boundary and way-points conditions [2, 3, 4]. In the presence of obstacles, sampling-based methods, including rapidly exploring random tree (RRT) and probabilistic roadmap (PRM), and virtual potential field methods are widely used to find obstacle-free trajectories [5, 2]. In [6], a mixed-integer optimization is provided for trajectory planning in the presence of convex obstacles. The proposed method first uses convex segmentation to compute convex regions of obstacle-free space. Then, it uses a mixed-integer optimization to assign polynomial trajectories to the computed convex safe regions. Also, [7] provides a moment-sum-of-squares-based convex optimization to obtain piece-wise linear trajectories in the presence of deterministic time-varying polynomial obstacles without the need for time discretization.

Trajectory planning problems under uncertainty look for trajectories with a bounded probability of collision with uncertain obstacles. Existing methods to address trajectory planning problems under uncertainty either are limited to Gaussian uncertainties and convex obstacles [8, 9, 10, 11, 12, 13, 14] or rely on sampling-based methods [15, 16, 17]. For example, chance constrained RRT$^*$ algorithm in [11] assumes Gaussian uncertainties and linear obstacles and performs a probabilistic collision check for the nodes of the search tree. Hence, it can not guarantee to satisfy the probabilistic safety constraints along the edges of the search tree. The Monte Carlo-based motion planning algorithms, e.g., [17], use a large number of uncertainty samples to estimate the probability of collision of a given trajectory. Sampling-based methods do not provide any analytical bounds on the probability of collision and, due to a large number of samples, can be computationally intractable.

Also, [18, 19] use moment-based approaches to address non-Gaussian uncertainties in motion planning problems in the presence of convex obstacles. More precisely, to obtain the risk bounded trajectories, [18] uses first and second-order moments of uncertainties and RRT$^*$ algorithm in the presence of linear obstacles and [19] uses higher-order moments and interior-point nonlinear optimization solvers in the presence of ellipsoidal obstacles.

Statement of Contributions: In this paper, we propose novel convex algorithms for risk bounded continuous-time trajectory planning in uncertain nonconvex static and dynamic environments that contain obstacles with probabilistic location, size, geometry, and trajectories with arbitrary probabilistic distributions. To achieve risk bounded plans:

1)	We provide an analytical method to compute risk contours maps. Risk contours allow us to identify risk bounded regions in uncertain environments and transform nonlinear stochastic planning problems into deterministic standard planning problems, in the presence of arbitrary probabilistic uncertainties. Hence, standard deterministic motion planning algorithms, e.g., RRT*, PRM, can be employed to look for safe (risk-bounded) plans.

2)	To solve the obtained deterministic planning problems, we provide two planners including i) sum-of-squares-based RRT algorithm and ii) sum-of-squares-based convex optimization that allows us to theoretically look for global optimal plans in nonconvex environments.

3)	We also provide continuous-time safety guarantees in stochastic environments. To ensure safety, existing planning under uncertainty algorithms “only” verify the safety of a finite set of waypoints (time-discretization). Unlike the existing planners, the provided planners of this paper ensure the (risk bounded) safety of the continuous-time trajectories without the need for time discretization.

The outline of the paper is as follows: Section II presents the notation adopted in the paper and definitions of polynomials, moments, and sum-of-squares optimization. In Section III, we provide the problem formulation of risk bounded continuous-time trajectory planning. Section IV provides analytical approaches to compute static and dynamic risk contours defined for static and dynamic uncertain obstacles to identify the risk bounded safe regions in uncertain environments. In Section V, using the obtained risk contours, we provide sum-of-squares-based planners to look for continuous-time risk bounded trajectories in uncertain static and dynamic environments. In Section VI, we present experimental results on the risk bounded planning problems of autonomous and robotic systems followed by a discussion section. Finally, concluding remarks and future work are given in Section VII.

§ NOTATION AND DEFINITIONS

This section covers notation and some basic definitions of polynomials, moments of probability distributions, and sum-of-squares optimization.
For a vector $\mathbf{x}\in\R^n$ and multi-index $\alpha\in\mathbb{N}^n$, let $x^\alpha = \prod_{i=1}^nx_i^{\alpha_i}$.

Polynomials: Let $\mathbb{R}[x]$ be the set of real polynomials in the variables $\mathbf{x} \in \mathbb{R}^n$. Given polynomial
$\mathcal{P}(\mathbf{x}):\mathbb{R}^n\rightarrow\mathbb{R}$, we represent $\mathcal{P}$ as $\sum_{\alpha\in\mathbb{N}^n} p_\alpha x^\alpha$ where
$\{x^\alpha\}_{\alpha\in \mathbb{N}^n}$ are the standard monomial basis of $\mathbb{R}[x]$, $\mathbf{p}=\{p_\alpha\}_{\alpha\in\mathbb{N}^n}$ denotes the coefficients, and $\alpha \in \mathbb N ^n$. In this paper, we use polynomials to describe uncertain obstacles and continuous-time trajectories. For example the set $\{ (x_1,x_2): 1-(x_1-\omega_1)^2-(x_2-\omega_2)^2 \geq 0\}$  represents a circle-shaped obstacle whose center is subjected to uncertainty modeled with random variables $(\omega_1,\omega_2)$. Also,
\end{matrix}\right]=\left[\begin{matrix}
-1\end{matrix}\right]t + \left[\begin{matrix}
2\end{matrix}\right]t^2$, $t\in [0,1]$ is an example of polynomial trajectory of order 2 in 2D environment between the points $\mathbf{x}(0)=(0,0)$ and $\mathbf{x}(1)=(1.5,1)$.

Moments of Probability Distributions:
Moments of random variables are the generalization of mean and covariance and are defined as
expected values of monomials of random variables. More precisely, given $(\alpha_1,...,\alpha_n) \in \mathbb{N}^n$ where $\alpha=\sum_{i=1}^{n}\alpha_i$, moment of order $\alpha$ of random vector $\mathbf{\omega}$ is defined as $\mathbb{E}[ \Pi_{i=1}^n \omega_i^{\alpha_i}]$. For example, sequence of the moments of order $\alpha=2$ for $n=3$ is defined as
$\left[\mathbb{E}[w_1^{2}], \mathbb{E}[\omega_1\omega_2],
\mathbb{E}[\omega_1\omega_3],
\mathbb{E}[\omega_2^2] ,
\mathbb{E}[\omega_2\omega_3],
\mathbb{E}[\omega_3^2] \right]$.
Moments of random variables can be easily computed using the characteristic function of probability distributions [20]. We will use a finite sequence of the moments to represent non-Gaussian probability distributions.

Sum of Squares Polynomials and Optimization: In this paper, we will use sum of squares (SOS) techniques to solve nonconvex optimization problems of the risk bounded trajectory planning problems. Polynomial $\mathcal{P}(x)$ is a sum of squares polynomial if it can be written as a sum of finitely many squared polynomials, i.e., $\mathcal{P}(x)= \sum_{j=1}^{m} h_j(x)^2$ for some $m<\infty$ and $h_j(x)\in\mathbb{R}[x]$ for $1\leq j\leq m$.
SOS condition, i.e., $\mathcal{P}(x) \in SOS$, can be represented as a convex constraint of the form of a linear matrix inequality (LMI) in terms of the coefficients of the polynomial, i.e.,
$\mathcal{P}(x) \in SOS \rightarrow \mathcal{P}(x)=\mathbf{x}^TA\mathbf{x}$ where $\mathbf{x}$ is the vector of standard monomial basis and $A$ is a positive semidefinite matrix in terms of the coefficients of the polynomial [21, 22, 23]. One can use different software packages like Yalmip [24] and Spotless [25] to check the SOS condition of the polynomials.

Sum of squares polynomials are used to obtain the convex relaxations of noncovex polynomial optimization problems [21, 22, 23]. More precisely, consider the following noncovex polynomial optimization

\begin{equation} \label{nlp}
\begin{aligned}
& \underset{\mathbf{x} \in \mathbb{R}^n}{\text{minimize}}
& & f(\mathbf{x}) \\
& \text{subject to}
& & g_i(\mathbf{x}) \geq 0, \ i=1,...,m \\
\end{aligned}
\end{equation}

where $f$ and $g_i, i=1,...,m$ are polynomial functions. We can rewrite the polynomial optimization in (<ref>) as the following form:

\begin{equation}\label{nlp_sos}
\begin{aligned}
& \underset{\gamma \in \mathbb{R}}{\text{maximize}}
& & \gamma \\
& \text{subject to}\
& & f(\mathbf{x})-\gamma \geq 0, \  \forall \mathbf{x} \in \{\mathbf{x} \in \mathbb{R}^n: g_i(\mathbf{x}) \geq 0 \  |_{i=1}^{m} \}
\end{aligned}
\end{equation}

where we look for the best lower bound of the function $f(\mathbf{x})$ denoted by $\gamma$ inside the feasible set of the original optimization problem. Hence, if $\mathbf{x}^*$ is the optimal solution of the original optimization problem in (<ref>), then $\gamma^* = f(\mathbf{x}^*)$ is the optimal solution of the optimization problem in (<ref>).

Note that in optimization problem (<ref>), objective function is linear and the constraint is the nonnegativity condition of  polynomial $f(\mathbf{x})-\gamma$. Such nonnegativity condition can be replaced by SOS conditions of polynomials [21, 22, 23]. Hence, we can transform the optimization in (<ref>) into a convex optimization, i.e., semidefinite program, with LMI constraints in terms of the coefficients of the polynomials of the original optimization in (<ref>). Also, we can recover the optimal solution of the original polynomial optimization in (<ref>), i.e., $\mathbf{x}^*$, using the solution of the dual convex optimization problem of (<ref>) as shown in [22, 23]. One can use different software packages like GloptiPoly [26] to solve the polynomial optimization in (<ref>) using SOS-based primal-dual approach.

Recently, SOS optimization techniques have been extended to address time-varying polynomial optimization problems of the form

\begin{equation} \label{nlp2}
\begin{aligned}
& \underset{\mathbf{x} \in \mathbb{R}^n}{\text{minimize}}
& & f(\mathbf{x}) \\
& \text{subject to}
& & g_i(\mathbf{x},t) \geq 0,  \ \forall t\in [t_0,t_f], \ i=1,...,m \\
\end{aligned}
\end{equation}

where one needs to make sure that time-varying constraints are satisfied over the given time horizon $t \in [t_0,t_f]$. SOS-based techniques can be used to solve the time-varying optimization problem in (<ref>) without the need for time discretization by transforming the problem into a convex optimization, i.e., time-varying semidefinite program [27, 7, 28].

§ PROBLEM FORMULATION

Suppose $\mathcal{X}\in \mathbb{R}^{n_x}$ is an uncertain environment and the sets $\mathcal{X}_{obs_i}(\omega_i) \subset \mathcal{X}, \ i=1,...,n_{o_s}$ and
$\mathcal{X}_{obs_i}(\omega_i,t) \subset \mathcal{X}, \ i=n_{o_s}+1,...,n_{o_d}$ are
the static and dynamic uncertain obstacles, respectively, where $\omega_i \in \mathbb{R}^{n_{\omega}}, \ i=1,...,n_{o_s}+n_{o_d}$
are probabilistic uncertain parameters with known probability distributions. We represent static uncertain obstacles in terms of polynomials in $\mathbf{x}\in \mathcal{X}$ and uncertain parameters as follows:

\begin{equation} \label{obs_s}
\mathcal{X}_{obs_i}(\omega_i)= \{ \mathbf{x}\in \mathcal{X}: \mathcal{P}_i(\mathbf{x},\mathbf{\omega}_i) \geq 0  \}, \ i=1,...,{n_{o_s}}
\end{equation}

where $\mathcal{P}_i:  \mathbb{R}^{n_x+n_{\omega}} \rightarrow \mathbb{R}, i=1,...,n_{o_s}$ are the given polynomials. Similarly, we represent dynamic uncertain obstacles in terms of polynomials in $\mathbf{x}\in \mathcal{X}$, time $t$, and uncertain parameters as follows:

\begin{equation} \label{obs_d}
\mathcal{X}_{obs_i}(\omega_i,t)= \{ \mathbf{x}\in \mathcal{X}: \mathcal{P}_i(\mathbf{x},\mathbf{\omega}_i, t) \geq 0  \} \ |_{i=n_{o_s}+1}^{n_{o_d}}
\end{equation}

where $\mathcal{P}_i:  \mathbb{R}^{n_x+n_{\omega}+1} \rightarrow \mathbb{R}, i=n_{o_s}+1,...,n_{o_d}$ are the given polynomials. Note that, in
general, the sets in (<ref>) and (<ref>) represent nonconvex probabilistic obstacles, e.g., nonconvex obstacles with uncertain size, location, or geometry [29, 30].

Given static and dynamic uncertain obstacles in (<ref>) and (<ref>), we define risk as the probability of collision with uncertain obstacles in the environment. In the risk bounded continuous-time trajectory planning problem, we aim at finding a continuous-time trajectory $\mathbf{x}(t): [t_0, t_f] \rightarrow \mathbb{R}^{n_x}$ defined over the time horizon $t\in[t_0,t_f]$ between the start and final points $\mathbf{x}_0$ and $\mathbf{x}_f$ such that the probability of collision of the trajectory $\mathbf{x}(t)$ with uncertain obstacles is bounded.

More precisely, we define the risk bounded continuous-time trajectory planning problem as the following probabilistic optimization problem:

\begin{align}
& \underset{\mathbf{x}(t): [t_0, t_f] \rightarrow \mathbb{R}^{n_x}}{\text{minimize}}
\int_{t_0}^{t_f}  \left\Vert \dot{\mathbf{x}}(t) \right\Vert_2^2 dt \label{opt_obj} \\
& \text{subject to}  \
\ \ \mathbf{x}(t_0)=\mathbf{x}_0, \  \mathbf{x}(t_f)=\mathbf{x}_f \subeqn \\
&  \ \hbox{Prob}\left( \mathbf{x}(t) \in \mathcal{X}_{obs_i}(\omega_i)    \right) \leq \Delta, \ \forall t\in [t_0,t_f] \ |_{i=1}^{n_{o_s}} \subeqn \label{opt_prob1} \\
& \  \hbox{Prob}\left(  \mathbf{x}(t) \in \mathcal{X}_{obs_i}(\omega_i,t)    \right) \leq \Delta, \ \forall t\in [t_0,t_f] \ |_{i=n_{o_s}+1}^{n_{o_d}} \subeqn \label{opt_prob2}
\end{align}

where objective function (<ref>) is the length of the trajectory $\mathbf{x}(t)$ defined in terms of $\ell_2$ norm $\left\Vert. \right\Vert_2$. Also, constraints (<ref>) and (<ref>) are the defined risks at time $t$ for trajectory $\mathbf{x}(t)$ in terms of the uncertain static and dynamic obstacles, respectively. Moreover, $ 0 \leq \Delta \in \mathbb{R} \leq 1 $ is the given acceptable risk level.
To solve the risk bounded optimization problem in (<ref>), we will look for the following continuous-time trajectories:

i) Polynomial trajectories over the planning horizon $[t_0,t_f]$ of the form

\begin{equation}\label{poly}
\mathbf{x}(t)=\sum_{\alpha=0}^d \mathbf{c}_{\alpha}t^{\alpha}, \  t\in[t_0,t_f]
\end{equation}

where $\mathbf{c}_{\alpha} \in \mathbb{R}^{n_x}, \ \alpha=0,...,d$ are the coefficient vectors,

ii) Piece-wise linear trajectories of the form

\begin{equation}\label{linear}
\mathbf{x}_i(t)= \mathbf{a}_i+ \mathbf{b}_it, \ t \in [t_{i-1}, t_i), \ i=1,...,s
\end{equation}
where $s$ is the number of linear pieces defined over the time intervals $t \in [t_{i-1},t_i] $ of the form $t_{i-1}=t_0+\frac{(i-1)(t_f-t_0)}{s}$ and $t_{i} = t_0+ \frac{i(t_f-t_0)}{s}, i=1,...,s$, and $\mathbf{a}_i,\mathbf{b}_i \in \mathbb{R}^{n_x}$ are the coefficient vectors.

Solving the probabilistic optimization in (<ref>) is challenging,  because i) we need to deal with multivariate integrals of the probabilistic constraints in (<ref>) and (<ref>) defined over the nonconvex sets of the obstacles, ii)
we need to deal with time-varying constraints to ensure that they are all satisfied over the entire planning time horizon $[t_0,t_f]$, and iii) optimization in (<ref>) is, in general, nonconvex optimization; Hence, we cannot guarantee to obtain the global optimal solution.

In this paper, we provide a systematic numerical procedure to efficiently solve the probabilistic optimization problem in (<ref>) in the presence of nonconvex uncertain obstacles with arbitrary probability distributions. For this purpose, we will leverage the notion of risk contours to transform the probabilistic optimization in (<ref>) into a deterministic polynomial optimization problem and use SOS optimization techniques to obtain optimal continuous-time trajectories with guaranteed bounded risk.

§ RISK CONTOURS

In [29], we define the risk contour with respect to the static uncertain obstacle in (<ref>) and the given acceptable risk level
$\Delta$ in (<ref>) as the set of all points in the environment whose probability of collision with the uncertain obstacle is less or equal to $\Delta$.

In this paper, we use static and dynamic risk contours defined for static and dynamic uncertain obstacles, respectively, to identify the safe regions in uncertain environments, i.e., the feasible set of optimization (<ref>).

In [29], to construct the risk contours of uncertain static obstacles, we propose an $(n_x+n_{\omega})$-dimensional convex optimization in the form of a semidefinite program (SDP). Such optimization is not suitable for online computations and is limited to small dimensions $(n_x+n_{\omega})$. In this paper, we propose an optimization-free fast approach, i.e., an analytical method, to construct the risk contours both for static and dynamic uncertain obstacles and show how one can use the obtained risk contours to solve the risk-bounded trajectory planning problem in (<ref>).

§.§ Static Risk Contours

Let $\mathcal{X}_{obs}(\omega)= \{ \mathbf{x}\in \mathcal{X}: \mathcal{P}(\mathbf{x},\mathbf{\omega}) \geq 0  \}$ be the given static uncertain obstacle as defined in (<ref>) and $\Delta \in [0,1]$ be the given acceptable risk level. Then, static $\Delta$-risk contour denoted by $\mathcal{C}^{\Delta}_{r}$ is defined as the set of all points in the environment, i.e., $\mathbf{x} \in \mathcal{X}$, whose risk is less or equal to $\Delta$. More precisely,

\begin{equation}\label{rc_s}
\mathcal{C}^{\Delta}_{r}:= \{\ \mathbf{x} \in \mathcal{X}: \ \hbox{Prob}( \mathbf{x} \in \mathcal{X}_{obs}(\omega)) \leq \Delta  \}
\end{equation}

The main idea to construct the static risk contour in (<ref>) is to replace the probabilistic constraint, i.e.,
$\hbox{Prob}( \mathbf{x} \in \mathcal{X}_{obs}(\omega))=\hbox{Prob}( \mathcal{P}(\mathbf{x},\omega)\geq 0 ) \leq \Delta$, with a deterministic constraint in terms of $\mathbf{x}$. In this paper, we provide an analytical method as follows:

Given the polynomial $\mathcal{P}(\mathbf{x},\mathbf{\omega})$ of the uncertain obstacle $\mathcal{X}_{obs}(\omega)$, we define the set $\hat{\mathcal{C}}_{r}^{\Delta}$ as follows:

\begin{align}\label{rc_s_cheby}
\hat{\mathcal{C}}^{\Delta}_{r}= \left\lbrace \ \mathbf{x} \in \mathcal{X}:
\begin{array}{cc} \frac{\Exp[\mathcal{P}^2(\mathbf{x},\omega)] - \Exp[\mathcal{P}(\mathbf{x},\omega)]^2}{\Exp[\mathcal{P}^2(\mathbf{x},\omega)]} \leq \Delta, \\
\begin{small}	\Exp[\mathcal{P}(\mathbf{x},\omega)]\leq 0  \end{small} \end{array}
\right\rbrace
\end{align}

where the expectation is taken with respect to the distribution of uncertain parameter $\omega$. Note that we can compute polynomials $\Exp[\mathcal{P}^2(\mathbf{x},\omega)]$ and $\Exp[\mathcal{P}(\mathbf{x},\omega)]$ in terms of $\mathbf{x}$ and known moments of $\omega$. More precisely, $\Exp[\mathcal{P}^2(\mathbf{x},\omega)]$ and $\Exp[\mathcal{P}(\mathbf{x},\omega)]$ are polynomials in $\mathbf{x}$ whose coefficients are defined in terms of the moments of $\omega$ and the coefficients of polynomial $\mathcal{P}(\mathbf{x},\omega)$.

The following result holds true.

Theorem 1: The set $\hat{\mathcal{C}}_{r}^{\Delta}$ in (<ref>) is an inner approximation of the static risk contour ${\mathcal{C}}_{r}^{\Delta}$ in (<ref>).

Proof: See Appendix. $\blacksquare$

Remark 1: The set in (<ref>) is a rational polynomial-based inner approximation of the risk contour in (<ref>). It also uses higher order moments of the uncertain parameter $\omega$ up to order $2d$ where $d$ is the order of the polynomial obstacle $\mathcal{P}(\mathbf{x},\omega)$.

Note that since $\hat{\mathcal{C}}^{\Delta}_{r}$ is an inner approximation of ${\mathcal{C}}^{\Delta}_{r}$, any trajectory $\mathbf{x}(t) \in \hat{\mathcal{C}}^{\Delta}_{r} \ \forall t\in [t_0,t_f]$ is guaranteed to have a risk less or equal to $\Delta$. We now provide an illustrative example to show the performance of the proposed method to construct the static $\Delta$-risk contours and benchmark our method against the optimization-based approach in [29].

Illustrative Example 1:

Consider the following illustrative example where $\mathcal{X}=[-1,1]^2$. The set $\mathcal{X}_{obs}(\omega)= \left\lbrace  (x_1,x_2) :  \omega^2-x_1^2-x_2^2 \geq 0 \right\rbrace $ represents a circle-shaped obstacle
whose radius $\omega$ has a uniform probability distribution over $[0.3,0.4]$, [29].
Moment of order $\alpha$ of a uniform distribution defined over $[l,u]$ can be described in a closed-form as $\frac{u^{\alpha+1}-l^{\alpha+1}}{(u-l)(\alpha+1)}$.
To construct the static $\Delta$-risk contour, we compute polynomials $\Exp[\mathcal{P}(\mathbf{x},\omega)]$
and $\Exp[\mathcal{P}^2(\mathbf{x},\omega)]$ using the polynomial obstacle and the moments of $\omega$ as follows:

\begin{align}
\Exp[\mathcal{P}(\mathbf{x},\omega)] &= \Exp[\omega^2-x_1^2-x_2^2 ]= \Exp[\omega^2]-x_1^2-x_2^2 =0.1233-x_1^2-x_2^2 \notag\\
\Exp[\mathcal{P}^2(\mathbf{x},\omega)] &= \Exp[\left(\omega^2-x_1^2-x_2^2 \right)^2]  \notag \\
&=\Exp[\omega^4] - 2\Exp[\omega^2]x_1^2 - 2\Exp[\omega^2]x_2^2 + x_1^4 + 2x_1^2x_2^2 + x_2^4\notag \\
&=0.0156 - 0.2466x_1^2 - 0.2466x_2^2 + x_1^4 + 2x_1^2x_2^2 + x_2^4 \notag
\end{align}

As shown in Figures <ref> and <ref>, we use the sublevels of functions $\frac{\Exp[\mathcal{P}^2(\mathbf{x},\omega)] - \Exp[\mathcal{P}(\mathbf{x},\omega)]^2}{\Exp[\mathcal{P}^2(\mathbf{x},\omega)]} $ and $\Exp[\mathcal{P}(\mathbf{x},\omega)]$ as in (<ref>) to construct the inner approximations of the static $\Delta$-risk contours for different risk levels $\Delta=[0.2,0.1,0.07,0.05]$.

We also compare our proposed method in  (<ref>) with the optimization-based method in [29] as shown in Figure <ref>.
Using the provided SDP in [29], we obtain a polynomial of order 20 to describe the inner approximations of the $\Delta$-risk contours. We note that the proposed analytical method of this paper provides a tight inner approximation of the $\Delta$-risk contours. This is primarily due to the facts that i) the proposed analytical approach results in a rational polynomial representation of the risk-contours as opposed to a standard polynomial representation provided in [29] and ii) with the provided analytical approach, we are able to avoid the numerical issues that arise when solving large-scale SDPs [31]. The proposed method of this paper is also suitable for online large scale planning problems.

Illustrative example 1: Intersection of $\Delta$-sublevel set of function $\frac{\Exp[\mathcal{P}^2(\mathbf{x},\omega)] - \Exp[\mathcal{P}(\mathbf{x},\omega)]^2}{\Exp[\mathcal{P}^2(\mathbf{x},\omega)]}$ and $0$-sublevel set of function $\Exp[\mathcal{P}(\mathbf{x},\omega)]$ describes the static $\Delta$-risk contour as in (<ref>).

Illustrative example 1: True static $\Delta$-risk contour ${\mathcal{C}}^{\Delta}_{r}$ (green) and inner approximation $\hat{\mathcal{C}}^{\Delta}_{r}$ obtained using i) the proposed analytical method in (<ref>) (outside of the solid-line) and ii) the proposed optimization-based method in ([29], Fig.4 ) (outside of the dashed-line). While the optimization-based method in [29] uses a standard polynomial of order 20, the proposed analytical method in (<ref>) uses 4th order rational and 2nd order standard polynomials to describe the risk contours.

§.§ Dynamic Risk Contours

Let $\mathcal{X}_{obs}(\omega,t)= \{ \mathbf{x}\in \mathcal{X}: \mathcal{P}(\mathbf{x},\mathbf{\omega},t) \geq 0  \}$ be the given dynamic uncertain obstacle as defined in (<ref>) and $\Delta \in [0,1]$ be the given acceptable risk level in (<ref>). Then, dynamic $\Delta$-risk contour at time $t$ denoted by $\mathcal{C}^{\Delta}_{r}(t)$ is defined as the set of all points in the environment, i.e., $\mathbf{x} \in \mathcal{X}$, whose risk at time $t$ is less or equal to $\Delta$. More precisely,

\begin{equation}\label{rc_d}
\mathcal{C}^{\Delta}_{r}(t):= \{\ \mathbf{x} \in \mathcal{X}: \ \hbox{Prob}( \mathbf{x} \in \mathcal{X}_{obs}(\omega,t)) \leq \Delta  \}
\end{equation}

Similar to the static risk contours, we can replace the probabilistic constraint in (<ref>) with deterministic constraints and construct an inner approximation of the dynamic $\Delta$-risk contour denoted by $\hat{\mathcal{C}}^{\Delta}_{r}(t)$ as follows:

\begin{align}\label{rc_d_cheby}
\hat{\mathcal{C}}^{\Delta}_{r}(t)= \left\lbrace \ \mathbf{x} \in \mathcal{X}:
\begin{array}{cc} \frac{\Exp[\mathcal{P}^2(\mathbf{x},\omega,t)] - \Exp[\mathcal{P}(\mathbf{x},\omega,t)]^2}{\Exp[\mathcal{P}^2(\mathbf{x},\omega,t)]} \leq \Delta, \\
\begin{small}	\Exp[\mathcal{P}(\mathbf{x},\omega,t)]\leq 0  \end{small} \end{array}
\right\rbrace
\end{align}

Note that dynamic $\Delta$-risk contour (<ref>) is described in terms of the time-varying constraints.

Illustrative Example 2:
Consider the following illustrative example where $\mathcal{X}=[-1,1]^2$. The set $\mathcal{X}_{obs}(\omega,t)= \left\lbrace  (x_1,x_2) :  \omega_1^2-(x_1-p_{x_1}(t,\omega_2))^2-(x_2-p_{x_2}(t,\omega_3))^2 \geq 0 \right\rbrace$  represents a moving circle-shaped obstacle with uncertain radius $\omega_1$ and uncertain trajectories $ p_{x_1}(t,\omega_2)=2-t+t^2+0.2\omega_2,\  p_{x_2}(t,\omega_3)=-1+4t-t^2+0.1\omega_3 $ that describe the uncertain motion of the obstacle over the time horizon $t\in[0,1]$. Uncertain parameters have uniform, normal, and Beta distributions as $\omega_1 \sim Uniform[0.3,0.4]$, $\omega_2 \sim \mathcal{N}(0,0.1)$, $\omega_3 \sim Beta(3,3)$. Moment of order $\alpha$ of a Beta distribution with parameters $(a,b)$
and a normal distribution with mean $\mu$ and standard deviation $\sigma$ can be described in closed-forms as $y_{\alpha}=\frac{a+\alpha-1}{a+b+\alpha-1}y_{\alpha-1}, y_0=1$ and $y_{\alpha}=\sigma^{\alpha}(-\sqrt{-1}\sqrt{2})^{\alpha}kummerU(\frac{-\alpha}{2}, \frac{1}{2}, \frac{-\mu^2}{2\sigma^2})$, respectively, where $kummerU$ is "confluent hyper-geometric Kummer U" function. Similar to illustrative example 1, we compute the polynomials $\Exp[\mathcal{P}^2(\mathbf{x},\omega,t)]$ and
$\Exp[\mathcal{P}(\mathbf{x},\omega,t)]$ using the moments of the uncertain parameters $\omega_i, i=1,2,3$ and the polynomial obstacle. We then construct the dynamic $\Delta$-risk contour as a function of time as described in (<ref>).
Figure <ref> shows the obtained dynamic $\Delta$-risk contours for $\Delta=0.1$ at time steps $t=0,0.5,1$ along the given uncertain trajectory $\left( p_{x_1}(t,\omega_2),\  p_{x_2}(t,\omega_3) \right)$.

Illustrative example 2:
a) Function $\frac{\Exp[\mathcal{P}^2(\mathbf{x},\omega,t)] - \Exp[\mathcal{P}(\mathbf{x},\omega,t)]^2}{\Exp[\mathcal{P}^2(\mathbf{x},\omega,t)]}$ at time steps $t=0,0.5,1$, b) Dynamic $\Delta$-risk contours for $\Delta=0.1$ at time steps $t=0,0.5,1$
described in (<ref>). Dashed line shows the expected value of the given uncertain trajectory, i.e., $\Exp[(p_{x_1}(t,\omega_2),p_{x_2}(t,\omega_3))]$. At each time $t$, for any point inside $\hat{\mathcal{C}}^{\Delta}_{r}(t)$ (outside of the closed curve), probability of collision with the moving uncertain obstacle is less or equal to $\Delta=0.1$.

Remark 2: We can use (<ref>)
and (<ref>) to construct static and dynamic risk contours in real-time.
Therefore, standard motion planning algorithms such as RRT$^*$, PRM, and motion primitive-based methods can be used for real-time risk bounded motion planning. To accomplish this, one just needs
to use the safe regions, i.e, risk contours, to construct the trajectories.

In the next section, we provide continuous-time planning algorithms to look for trajectories with guaranteed bounded risk over the entire planning time horizon without the need for time discretization.

§ CONTINUOUS-TIME RISK BOUNDED TRAJECTORY PLANNING USING RISK CONTOURS

In this section, we will use the static and dynamic risk contours to solve the continuous-time risk bounded trajectory planning problem defined in (<ref>). More precisely, we use the static and dynamic risk contours to transform the probabilistic optimization in (<ref>) into a deterministic polynomial optimization.
The obtained deterministic polynomial optimization is, in general, nonconvex and nonlinear.
In addition, we need to ensure that all the obtained deterministic constraints are satisfied over the entire planning time horizon $[t_0,t_f]$. In this section, we provide convex methods based on SOS techniques introduced in Section <ref> to efficiently solve the obtained nonconvex time-varying deterministic planning optimization problem. While the existing planners rely on time discretization to verify the planning safety constraints, the provided SOS-based planners look for continuous-time trajectories with guaranteed bounded risk over the entire planning time horizon without the need for time discretization.

We first begin by addressing the continuous-time risk bounded trajectory planning in static uncertain environments using the static risk contours. We then use the dynamic risk contours to address the planning problems in dynamic uncertain environments.

§.§ Planning in Static Uncertain Environments

In this section, we are concerned with continuous-time risk bounded trajectory planning in the presence of static uncertain obstacles of the form (<ref>). More precisely, we aim at solving
the probabilistic optimization problem in (<ref>) in the presence of static uncertain obstacles $\mathcal{X}_{obs_i}(\omega_i), i=1,...,n_{o_s}$, i.e.,

\begin{align}
& \underset{\mathbf{x}(t): [t_0, t_f] \rightarrow \mathbb{R}^{n_x}}{\text{minimize}}
\int_{t_0}^{t_f}  \left\Vert \dot{\mathbf{x}}(t) \right\Vert_2^2 dt \label{opt_s_obj} \\
& \text{subject to} \
\ \ \mathbf{x}(t_0)=\mathbf{x}_0, \  \mathbf{x}(t_f)=\mathbf{x}_f \subeqn \\
&  \ \hbox{Prob}\left( \mathbf{x}(t) \in \mathcal{X}_{obs_i}(\omega_i)    \right) \leq \Delta, \ \forall t\in [t_0,t_f], \  i=1,...,n_{o_s} \subeqn \label{opt_s_prob1}
\end{align}

To obtain the deterministic polynomial optimization, we replace probabilistic constraints (<ref>) with deterministic constraints in terms of the static $\Delta$-risk contours as follows:

\begin{align}
& \underset{\mathbf{x}(t): [t_0, t_f] \rightarrow \mathbb{R}^{n_x}}{\text{minimize}}
\int_{t_0}^{t_f}  \left\Vert \dot{\mathbf{x}}(t) \right\Vert_2^2 dt \label{opt_s2_obj} \\
& \text{subject to} \
\ \ \mathbf{x}(t_0)=\mathbf{x}_0, \  \mathbf{x}(t_f)=\mathbf{x}_f \subeqn \\
&  \mathbf{x}(t) \in \hat{\mathcal{C}}^{\Delta}_{r_i}   , \ \forall t\in [t_0,t_f], \  i=1,...,n_{o_s} \subeqn \label{opt_s2_prob1}
\end{align}

where $\hat{\mathcal{C}}^{\Delta}_{r_i}$ is the static $\Delta$-risk contour of the static uncertain obstacle $\mathcal{X}_{obs_i}(\omega_i)$. The set of $\hat{\mathcal{C}}^{\Delta}_{r_i}, i=1,...,n_{o_s}$ represents the inner approximation of the feasible set of the probabilistic optimization in (<ref>). The obtained deterministic optimization in (<ref>) is time-varying optimization problem where we need to ensure to satisfy the constraints over the entire planning time horizon $[t_0,t_f]$.

To solve the deterministic optimization in (<ref>), we will look for i) polynomial trajectory defined in (<ref>) and ii) piece-wise linear trajectory defined in (<ref>).

By substituting the polynomial trajectory $\mathbf{x}(t)=\sum_{\alpha=0}^d \mathbf{c}_{\alpha}t^{\alpha}$ in optimization (<ref>), we obtain a deterministic optimization with constraints in terms of the coefficients of the polynomial trajectory $\mathbf{x}(t)$, i.e.,
$(1-\Delta)\Exp[\mathcal{P}^2(\sum_{\alpha=0}^d \mathbf{c}_{\alpha}t^{\alpha},\omega)] - \Exp[\mathcal{P}(\sum_{\alpha=0}^d \mathbf{c}_{\alpha}t^{\alpha},\omega)]^2 \leq 0$, $ \Exp[\mathcal{P}(\sum_{\alpha=0}^d \mathbf{c}_{\alpha}t^{\alpha},\omega)]\leq 0$.

Similarly, by substituting the piece-wise linear trajectory of (<ref>), we obtain the objective function $\sum_{j=1}^s \int_{t_{j-1}}^{t_{j}}  \left\Vert \dot{\mathbf{x}}_j(t) \right\Vert_2^2 dt $
and constraints
$ \mathbf{a}_j+ \mathbf{b}_jt \in \hat{\mathcal{C}}^{\Delta}_{r_i}, \ \forall t\in [t_{j-1},t_j], \ i=1,...,n_{o_s}, j=1,...,s $. To solve the obtained time-varying deterministic polynomial optimization problems, we will provide 3 methods based on SOS optimization techniques as follows:

§.§.§ Time-Varying SOS Optimization

we use time-varying SOS optimization, introduced in Section <ref>, to solve the time-varying deterministic polynomial optimization in (<ref>).
In implementation, we use the heuristic algorithm introduced by [7].

§.§.§ Standard SOS Optimization

In this method, we obtain a standard polynomial optimization of the form (<ref>) by eliminating time $t$. We then use the standard SOS optimization technique in (<ref>) to solve the obtained optimization problem. To achieve this, let $\mathcal{X}_{obs}(\omega)= \{ \mathbf{x}\in \mathcal{X}: \mathcal{P}(\mathbf{x},\mathbf{\omega}) \geq 0  \}$ be the given static uncertain obstacle defined in (<ref>) and $\mathbf{x}(t)$ be the polynomial trajectory in (<ref>). To eliminate time $t$, instead of using the instant risk as in (<ref>), we work with the average risk defined as
$\frac{1}{t_f-t_0} \int \int_{ \{ (t,\mathbf{\omega}): \mathcal{P}(\mathbf{x}(t),\omega)\geq 0\}} pr(\mathbf{\omega})d\mathbf{\omega}dt$ where $pr(\mathbf{\omega})$ is the probability density function of $\omega$. This is equivalent to treating $t$ as a random variable with a uniform probability distribution over the planning time horizon $[t_0,t_f]$.
We should note that the average risk is a weak safety measure, which means that just by bounding the average risk, we cannot guarantee to satisfy the constraints of the probabilistic optimization in (<ref>).

By defining the average risk, we substitute the trajectory $\mathbf{x}(t)$ in the probabilistic constraint and follow the same steps as in Section <ref> to obtain an upper bound of the average risk and construct the set of all coefficients $\mathbf{c}_{\alpha}|_{\alpha =0}^d$ that results in a risk bounded trajectory of the form $\mathbf{x}(t)=\sum_{\alpha=0}^d \mathbf{c}_{\alpha}t^{\alpha}$.
More precisely, we obtain the following set for the coefficients:

\left\lbrace \ \mathbf{c}_{\alpha}|_{\alpha =0}^d:\begin{array}{cc}(1-\Delta)\Exp[\mathcal{P}^2(\sum_{\alpha=0}^d \mathbf{c}_{\alpha}t^{\alpha},\omega)] - \Exp[\mathcal{P}(\sum_{\alpha=0}^d \mathbf{c}_{\alpha}t^{\alpha},\omega)]^2 \leq 0, \\ \Exp[\mathcal{P}(\sum_{\alpha=0}^d \mathbf{c}_{\alpha}t^{\alpha},\omega)]\leq 0\end{array}
\right\rbrace

where the expectation is taken with respect to the probability distributions of $\omega$ and $t$. We can obtain a similar set for the coefficients of the piece-wise linear trajectories, as well. By computing such sets for coefficients of the trajectories, we can transform the probabilistic optimization problem into a deterministic standard polynomial optimization and use the standard SOS optimization technique in (<ref>) to obtain the optimal solution.

§.§.§ RRT-SOS

In this method, we use sampling-based motion planning algorithms like RRT to construct the risk bounded trajectory of the deterministic polynomial optimization in (<ref>). To ensure safety along the edges of the RRT, we use an SOS-based continuous-time technique to verify the constraints of the optimization in (<ref>) as follows:

Let $\mathbf{x}(t)=\sum_{\alpha=0}^d \mathbf{c}_{\alpha}t^{\alpha}, t\in[t_1,t_2]$ be the given trajectory between the two samples $\mathbf{x}_1 \in \mathcal{X}$ and $\mathbf{x}_2 \in \mathcal{X}$ in the uncertain environment. Also, let $\mathcal{S}=\{\mathbf{x}: g_i(x)\geq 0 \ \ i=1,...,\ell \}$ be the feasible set of optimization (<ref>), i.e., the set constructed by the polynomial constraints of all the risk contours $\hat{\mathcal{C}}^{\Delta}_{r_i}, i=1,...,n_{o_s}$. Then, the following result holds true:

Polynomial $\mathbf{x}(t)$ satisfies the safety constraints of the deterministic optimization in (<ref>) over the time interval $[t_1,t_2]$, i.e., $\mathbf{x}(t) \in \mathcal{S}$ for all $t\in[t_1,t_2]$, if and only if polynomials $g_i(\mathbf{x}(t)), i=1,...,\ell$ take the following SOS representation:

\begin{align}\label{SOS_Cond}
g_i(\mathbf{x}(t))={\sigma_0}_i(t)+ {\sigma_1}_i(t)(t-t_1) + {\sigma_2}_i(t)(t_2-t) \ |_{i=1}^{\ell}
\end{align}

where ${\sigma_0}_i(t), {\sigma_1}_i(t), {\sigma_2}_i(t), i=1,...,\ell$ are SOS polynomials with appropriate degrees [32, 22, 23]. Yalmip [24] and Spotless [25] packages can be used to check the SOS condition (<ref>) for the given trajectory $\mathbf{x}(t)$.

Remark 3: The complexity of the safety SOS condition in (<ref>) is independent of the size of the planning time horizon $[t_1,t_2]$ and the length of the polynomial trajectory $\mathbf{x}(t)$. Hence, one can use (<ref>) to easily verify the safety of trajectories in uncertain environments over the long planning time horizons.

The safety condition in (<ref>) can be used in any sampling-based motion planning algorithms to verify the safety of the trajectories between the sample points. In this paper, we will use the following naive RRT algorithm:
To expand the RRT, we use a linear trajectory to connect the given sample point to the closest vertex in the tree if the linear trajectory satisfies the SOS condition in (<ref>), i.e., this implies that the linear trajectory between the two points is inside the risk contours in (<ref>).

As a termination condition, we also check the safety of the linear trajectory
between the selected safe sample and the goal points via the SOS condition in (<ref>). If the linear trajectory satisfies SOS condition (<ref>), we connect the sample point to the goal point; Hence, a feasible trajectory between the start and goal points can be constructed.

To improve the obtained feasible trajectory, we perform the following steps: i) given the obtained tree, we construct a graph, e.g., PRM, whose nodes are the vertex of the tree and edges of the graph are all the linear trajectories between the nodes that satisfy the SOS condition in (<ref>), ii) we then perform a shortest path algorithm, e.g.,  Dijkstra algorithm, to obtain a path from the start to the goal point.
We can also use smart initialization to guide the RRT-SOS algorithm and improve the run-time. For example, we use the straight line between the start and goal points to initialize the  RRT-SOS algorithm. We then perform the sampling in the neighborhood of the given initial path and  incrementally increase the size of the neighborhood until a feasible trajectory is obtained.

Illustrative Example 3: Consider the uncertain obstacle in illustrative example 1. We want to find a risk bounded trajectory between the points $\mathbf{x}(0)=[-1,-1]$ and $\mathbf{x}(1)=[1,1]$ by solving the probabilistic optimization in (<ref>) with $\Delta=0.1$. For this purpose, we solve the deterministic optimization problem in (<ref>) with respect to the $0.1$-risk contour of the uncertain obstacle using the discussed 3 methods as shown in Figure <ref>. More precisely, we use standard SOS optimization to obtain a polynomial trajectory of order 2 and also a piece-wise linear trajectory with 2 pieces. Using GloptiPoly package, we extract two risk bounded trajectories between the given start and goal points as shown in Figure <ref>-a and <ref>-b. We also use RRT-SOS algorithm and time-varying SOS optimization to obtain piece-wise linear trajectories as shown in Figure <ref>-c and <ref>-d, respectively. We note that the standard SOS optimization-based method results in conservative trajectories as shown in Figure <ref>-a and <ref>-b.

We compare our proposed methods with Monte Carlo-based risk bounded RRT algorithm that uses uncertainty samples and time discretization to look for the risk bounded trajectories. To verify the safety of the edges in the RRT, we use only $10$ uncertainty samples and $20$ uniformly sampled way-points on the edges. Such RRT algorithm is significantly slower and also does not provide any guaranteed risk bounded trajectories.

The run-time for the standard SOS optimization to obtain the piece-wise linear and polynomial trajectories are both roughly $1.5$ seconds. Also, the run-time for the time-varying SOS optimization is roughly $2.2$ seconds (for more information see Section <ref>). The run-time for the RRT-SOS algorithm to obtain feasible and optimal trajectories are roughly $9.3$ and $92.1$ seconds, respectively.
The run-time for the Monte Carlo-based risk bounded RRT algorithm to obtain feasible and optimal trajectories are roughly $215.12$ and $3471.4$ seconds, respectively. Also, continuous-time safety verification of each edge in the RRT-SOS algorithm via (<ref>) takes roughly $0.1$ seconds while the sampling-based safety verification in the Monte Carlo-based risk bounded RRT algorithm takes roughly $3$ seconds.

Illustrative Example 3: a) risk bounded polynomial trajectories of order 2 obtained via standard SOS optimization, b) risk bounded piece-wise linear trajectories obtained via standard SOS optimization, c) risk bounded piece-wise linear trajectory obtained via RRT-SOS algorithm, d) risk bounded piece-wise linear trajectory obtained via time-varying SOS optimization.

§.§ Planning in Dynamic Uncertain Environments

In this section, we are concerned with continuous-time risk bounded trajectory planning in the presence of dynamic uncertain obstacles of the form (<ref>). More precisely, we aim at solving
the optimization problem in (<ref>) in the presence of dynamic uncertain obstacles $\mathcal{X}_{obs_i}(\omega_i,t), i=1,...,n_{o_d}$, i.e.,

\begin{align}
& \underset{\mathbf{x}(t): [t_0, t_f] \rightarrow \mathbb{R}^{n_x}}{\text{minimize}}
\int_{t_0}^{t_f}  \left\Vert \dot{\mathbf{x}}(t) \right\Vert_2^2 dt \label{opt_d_obj} \\
& \text{subject to} \
\ \ \mathbf{x}(t_0)=\mathbf{x}_0, \  \mathbf{x}(t_f)=\mathbf{x}_f \subeqn \\
&  \ \hbox{Prob}\left( \mathbf{x}(t) \in \mathcal{X}_{obs_i}(\omega_i,t)    \right) \leq \Delta, \ \forall t\in [t_0,t_f] \ |_{i=n_{o_s}+1}^{n_{o_d}} \subeqn \label{opt_d_prob1}
\end{align}

To obtain the deterministic optimization, we replace probabilistic constraints (<ref>) with deterministic constraints in terms of the dynamic $\Delta$-risk contours as follows:

\begin{align}
& \underset{\mathbf{x}(t): [t_0, t_f] \rightarrow \mathbb{R}^{n_x}}{\text{minimize}}
\int_{t_0}^{t_f}  \left\Vert \dot{\mathbf{x}}(t) \right\Vert_2^2 dt \label{opt_d2_obj} \\
& \text{subject to} \
\ \ \mathbf{x}(t_0)=\mathbf{x}_0, \  \mathbf{x}(t_f)=\mathbf{x}_f \subeqn \\
&  \mathbf{x}(t) \in \hat{\mathcal{C}}^{\Delta}_{r_i}(t)   , \ \forall t\in [t_0,t_f], \ i=n_{o_s}+1,...,n_{o_d} \subeqn \label{opt_d2_prob1}
\end{align}

where $\hat{\mathcal{C}}^{\Delta}_{r_i}(t)$ is the dynamic $\Delta$-risk contour of the dynamic uncertain obstacle $\mathcal{X}_{obs_i}(\omega_i,t)$. The set of $\hat{\mathcal{C}}^{\Delta}_{r_i}(t), i=o_{n_s}+1,...,o_{n_d}$ represents the inner approximation of the feasible set of the probabilistic optimization in (<ref>). In the deterministic polynomial optimization of (<ref>), we need to make sure that the constraints are satisfied over the entire planning time horizon $[t_0,t_f]$. To solve the time-varying deterministic polynomial optimization in (<ref>), we will use i) time-varying SOS optimization described in Section <ref> and also ii) sampling-based motion planning algorithm similar to the one in Section <ref> that uses the SOS-based continuous-time safety verification.

Illustrative Example 4: Risk bounded trajectories between the start (square) and goal (triangle) points obtained via the time-varying SOS optimization (solid line) and RRT-SOS algorithm (dashed line) in the presence of the probabilistic moving obstacle.

In the sampling-based algorithm, given the dynamic nature of the environment, one needs to verify the safety of the trajectory between the given two samples only for the time interval that is needed to traverse between the points. In this paper, we will use the following naive RRT-SOS algorithm:

RRT algorithm looks for piece-wise linear trajectories with $s$ number of linear pieces defined over the time intervals $\Delta t_i = [t_{i-1},t_i), i=1,...,s$ as in (<ref>). To construct the tree, we first fix the number of linear pieces $s$ and build the tree for each time interval incrementally. To expand the tree for the time interval $\Delta t_i$, we connect the given sample point to a vertex in the tree, built for the time interval $\Delta t_{i-1}$, if the linear trajectory between the two points satisfies SOS condition (<ref>) for the time interval $\Delta t_i$.
Moreover, to construct the tree for the time interval $\Delta t_{s-1}$, we verify the safety of the linear trajectory between the given sample and the goal point for the time interval $\Delta t_s$ as well.

Illustrative Example 4: Consider the uncertain moving obstacle in illustrative example 2. We want to find a risk bounded trajectory between the points $\mathbf{x}(0)=[1,-2]$ and $\mathbf{x}(1)=[3,2]$ by solving the probabilistic optimization in (<ref>) with $\Delta=0.1$. For this purpose, we solve the deterministic optimization problem in (<ref>) with respect to the dynamic $0.1$-risk contour of the uncertain obstacle via the time-varying SOS optimization and RRT-SOS algorithm as shown in Figure <ref>. Note that although the obtained RRT-SOS trajectory looks like a straight line, it is a piece-wise linear trajectory consisting of 2 pieces with different velocities. The run-time for the time-varying SOS optimization and RRT-SOS algorithm are roughly $6.9$ and $0.5$ seconds, respectively.

§ EXPERIMENTS

In this section, numerical examples are presented to illustrate the performance of the proposed approaches. To obtain the risk bounded continuous-time trajectories, we use the time-varying SOS optimization and the RRT-SOS algorithm described in Section <ref>[github.com/jasour/Risk-Bounded-Continuous-Time-Trajectory-Planning]. Note that the provided RRT-SOS algorithm uses a naive tree search as explained in Section <ref>.

The main objective of the provided RRT-SOS algorithm is to demonstrate how the risk contours in (<ref>) and (<ref>) and the SOS-based continuous-time safety constraints in (<ref>) can be incorporated into sampling-based motion planning algorithms to look for guaranteed risk bounded continuous-time trajectories in stochastic environments. The computations in this section were performed on a computer with Intel i7 2.7 GHz processors and 8 GB RAM. We use the Spotless MATLAB package [25] to verify the SOS-based continuous-time safety constraints in the RRT-SOS algorithm and the Julia package provided by [7] to solve the time-varying SOS optimization.

§.§ Risk Contours

The purpose of this example is to demonstrate how the provided approach can be used to compute the risk contours in the presence of highly complex uncertain unsafe regions.

§.§.§ 2D Uncertain Obstacle

Static uncertain obstacle of the form (<ref>) is described by the polynomial $\mathcal{P}(\mathbf{x},\omega)=- 0.42x_1^5 - 1.18x_1^4x_2 - 0.47x_1^4 + 0.3x_1^3x_2^2 - 0.57x_1^3x_2 + 0.6x_1^3 - 0.65x_1^2x_2^3 + 0.17x_1^2x_2^2 + 1.87x_1^2x_2 + 0.06x_1^2 + 0.69x_1x_2^4 - 0.14x_1x_2^3 - 0.85x_1x_2^2 + 0.6x_1x_2 - 0.21x_1 + 0.01x_2^5 - 0.06x_2^4 - 0.07x_2^3 - 0.41x_2^2 - 0.08x_2 +0.07 - 0.1w
$ where the uncertain parameter $\omega$ has a Beta distribution with parameters $(9,0.5)$ over $[0,1]$. Figure <ref> shows the uncertain obstacle for different values of the uncertain parameter $\omega$. We obtain the static risk contours defined in (<ref>) for different risk levels $\Delta$ as shown in Figure <ref> and <ref>.

Example A-1: a) Nonconvex obstacle with uncertain parameter $\omega \sim Beta(9,0.5)$, b) Obtained static risk contours $\hat{\mathcal{C}}^{\Delta}_r$ defined in (<ref>) for different risk levels $\Delta=[0.3,0.1,0.05,0.02,0.01]$.

Example A-1: True $\Delta$-risk contour ${\mathcal{C}}^{\Delta}_{r}$ (green) and its inner approximation $\hat{\mathcal{C}}^{\Delta}_{r}$ (outside of the solid-line). Inside the risk contours, probability of collision with the uncertain obstacle is less or equal to the associated risk level $\Delta$.

§.§.§ 3D Uncertain Obstacle

Static uncertain obstacle of the form (<ref>) is described by the polynomial $\mathcal{P}(\mathbf{x},\omega)=
$ where the uncertain parameter $\omega$ has a normal distribution with mean $0.1$ and variance $0.001$, [30]. Figure <ref> shows the uncertain obstacle for different values of the uncertain parameter $\omega$. We obtain the static risk contours defined in (<ref>) for different risk levels $\Delta$ as shown in Figure <ref>.

Example A-2: a) Nonconvex  obstacle with uncertain parameter $\omega \sim \mathcal{N}(0.1,0.001)$, b) Obtained static risk contours $\hat{\mathcal{C}}^{\Delta}_r$ defined in (<ref>) for different risk levels $\Delta=[0.5,0.3,0.1,0.05]$.

§.§ Risk Bounded Lane Changing for Autonomous Vehicles

In this example, we generate a risk bounded trajectory for the lane-change maneuver of an autonomous vehicle in the presence of surrounding vehicles. In this scenario, uncertain locations of the surrounding vehicles are modeled as the following sets:
$\mathcal{X}_{obs_1}(\omega_1,t)=\{(x_1,x_2): 0.3^2 - (x_1-p_{1}(t,\omega_1))^2 - (x_2-1)^2 \geq 0 \}$ and
$\mathcal{X}_{obs_2}(\omega_2,t)=\{(x_1,x_2): 0.3^2 - (x_1-p_{2}(t,\omega_1))^2 - x_2^2 \geq 0 \}$
$p_{1}(t,\omega_1)=t+0.4+\omega_1$ and $p_{2}(t,\omega_2)=2t+0.6+\omega_2$ are the uncertain trajectories of the surrounding vehicles with uncertain parameters $\omega_i\sim Uniform[-0.1,0.1],i=1,2$. For the lane-change maneuver, we look for a risk bounded trajectory between the points $\mathbf{x}(0)=(0,0)$ and $\mathbf{x}(1)=(2,0)$ over the planning time horizon $t\in[0,1]$. Figure <ref> shows the obtained trajectory using the time-varying SOS optimization and RRT-SOS algorithm considering the dynamic $0.1$-risk contours of the surrounding vehicles.
The run-time for the time-varying SOS optimization and RRT-SOS are roughly 1.2 and 6.5 seconds, respectively.

Example B: Risk bounded trajectories for the lane-change maneuver between the start (square) and goal (triangle) points obtained via the time-varying SOS optimization (solid line) and RRT-SOS algorithm (dotted line) in the presence of surrounding vehicles with uncertain trajectories. Dashed lines show the expected values of the uncertain trajectories of the surrounding vehicles.

§.§ Risk Bounded Trajectory Planning for Delivery Robots

In this example, we generate a risk bounded trajectory for a delivery robot in the presence of uncertain moving obstacles. In this scenario, uncertain locations of the moving obstacles are modeled as the following sets: $\mathcal{X}_{obs_1}(\omega_1,t)=\{(x_1,x_2): 0.4^2 - (-x_1+p_{1}(t,\omega_1))^2 - (x_2-1)^2 \geq 0 \}$, $\mathcal{X}_{obs_2}(\omega_2,t)=\{(x_1,x_2): 0.4^2 - (x_1+p_{2}(t,\omega_1))^2 - (x_2-2)^2 \geq 0 \}$, and $\mathcal{X}_{obs_3}(\omega_3,t)=\{(x_1,x_2): 0.4^2 - (-x_1+p_{3}(t,\omega_3))^2 - (x_2-3)^2 \geq 0 \}$ where $p_{1}(t,\omega_1)=t - 0.5+\omega_1$, $p_{2}(t,\omega_1)= 2t-\omega_2-0.8$, and $p_{3}(t,\omega_3)=1.8t - 0.7+\omega_3$ are the uncertain trajectories of the moving obstacles with uncertain parameters $\omega_i\sim Uniform[-0.1,0.1],i=1,2,3$. We look for a risk bounded trajectory between the start and destination points, $\mathbf{x}(0)=(0,0)$ and $\mathbf{x}(1)=(0,4)$, over the planning time horizon $t\in[0,1]$. Figure <ref> shows the obtained trajectories using the time-varying SOS optimization and RRT-SOS algorithm considering the dynamic $0.1$-risk contours of the moving obstacles.
The run-time for the time-varying SOS optimization and RRT-SOS are roughly 7 and 197 seconds, respectively.

Example C: Risk bounded trajectories for a delivery robot between the start (square) and goal (triangle) points obtained via the time-varying SOS optimization (solid line) and RRT-SOS algorithm (dotted line) in the presence of moving uncertain obstacles.

§.§ Planning in Cluttered Uncertain Environments

In this example, we look for risk bounded trajectories in cluttered static and dynamic uncertain environments.

§.§.§ 2D Static Environment

In this scenario, we use the RRT-SOS algorithm to obtain a risk bounded trajectory with $\Delta=0.1$ between the start and goal points, $\mathbf{x}(0)=(0,0)$, $\mathbf{x}(1)=(5,5)$, in the presence of circle-shaped obstacles with uncertain position as shown in Figure <ref>. The position of the obstacles is subjected to an additive zero mean Gaussian noise with 0.001 variance.
The run-time to obtain the feasible and optimal trajectories are roughly 27 and 92 seconds, respectively.

Example D-1: Risk bounded trajectory between the start (square) and goal (triangle) points  obtained via the RRT-SOS algorithm in the static uncertain cluttered environment.

§.§.§ 3D Dynamic Environment

In this scenario, we use the RRT-SOS algorithm to obtain a risk bounded trajectory with $\Delta=0.1$ between the start and goal points, $\mathbf{x}(0)=(-1,-1,0)$, $\mathbf{x}(1)=(1,1,1)$, in the presence of moving sphere-shaped obstacles with uncertain radius and uncertain trajectories as shown in Figure <ref>. Radius of the obstacles has a uniform distribution over $[0.1,0.2]$. Also, trajectories of the obstacles are subjected to an additive zero mean Gaussian noise with 0.001 variance. The run-time to obtain the risk bounded trajectory is roughly 29 seconds.

Example D-2: Risk bounded trajectory between the start (square) and goal (triangle) points  obtained via the RRT-SOS algorithm in the dynamic uncertain cluttered environment.

§.§ Discussion

The provided risk bounded algorithms not only are faster than Monte Carlo-based risk bounded RRT algorithms (illustrative example 3), but also provide continuous-time trajectories with guaranteed bounded risk. In addition, compared with the RRT-SOS algorithm, the time-varying SOS optimization is generally faster. The run-time of the time-varying SOS optimization is a function of the order of polynomials of risk contours, the number of the linear pieces of risk bounded trajectories, and the number of iterations of the heuristic algorithm. We use the heuristic algorithm introduced by [7] that trades off theoretical guarantees for more efficiency to avoid large scale time-varying SOS optimization. Hence, the heuristic algorithm may fail to obtain a feasible trajectory as observed in Experiment D. To ensure safety, we verify the trajectory returned by the heuristic algorithm using the SOS condition in (<ref>). In the time-varying SOS optimization, one can initially start with a small number of iterations and linear pieces and then increase the parameters if the obtained trajectory is infeasible. On the other hand, the RRT-SOS algorithm is more robust and always returns a feasible trajectory. The provided SOS-based RRT algorithm uses a naive tree search algorithm to look for the risk bounded trajectories. One can improve the performance and run-time by considering efficient sampling-based motion planning algorithms and high-order polynomial trajectories between the samples.

§ CONCLUSION

In this paper, we provided continuous-time trajectory planning algorithms to obtain risk bounded polynomial trajectories in uncertain nonconvex environments that contain static and dynamic obstacles with probabilistic location, size, and geometry. The provided algorithms  leverage  the  notion  of  risk  contours to transform the probabilistic trajectory planning problem into a deterministic planning problem and use convex methods to obtain the continuous-time trajectories with guaranteed bounded risk without  the  need  for  time discretization and uncertainty samples. The provided algorithms are suitable for online and large scale planning problems. For the future work, we will incorporate the provided approaches into a model predictive control (MPC) framework to address the online planning problems. Also, we will extend the provided approaches to deal with inaccurate models of uncertain parameters, e.g., inaccurate probability distributions and moments.

§ APPENDIX: PROOF OF THEOREM 1

To compute the static risk contour in (<ref>), we obtain the deterministic constraint as follows:
For a given point $\mathbf{x} \in \mathcal{X}$, the probability of collision with the uncertain static obstacle, i.e., $\hbox{Prob}( \mathbf{x} \in \mathcal{X}_{obs}(\omega))$, is equivalent to the expectation of the indicator function of the superlevel set of $\mathcal{P}(\mathbf{x},\omega)$ as follows [29, 30]:

\begin{equation}
\hbox{Prob}( \mathcal{P}(\mathbf{x},\omega)\geq 0 ) = \int_{ \{ (\mathbf{x},\mathbf{\omega}): \mathcal{P}(\mathbf{x},\omega)\geq 0\}} pr(\mathbf{\omega})d\mathbf{\omega}= \mathbb{E}[\mathbb{I}_{\mathcal{P}\geq 0}]
\end{equation}
where $pr(\mathbf{\omega})$ is the probability density function of $\mathbf{\omega}$ and $\mathbb{I}_{\mathcal{P}\geq 0}$ is the indicator function of the superlevel set of $\mathcal{P}(\mathbf{x},\omega)$ defined as $\mathbb{I}_{\mathcal{P}\geq 0} = 1 $ if $(\mathbf{x},\omega) \in  \{ (\mathbf{x},\omega): \mathcal{P}(\mathbf{x},\omega) \geq 0\} $, and
0 otherwise. The expectation of the indicator function, however, is not necessarily easily computable. To compute the expectation value, we will find a polynomial description of the indicator function that upper bounds the true indicator function $\mathbb{I}_{\mathcal{P}\geq 0}$. If we can find a polynomial $P_{\mathbb{I}}: \R^{n_x+n_{\omega}}\rightarrow\R$ of the order $d$ with coefficients $c_{ij}$ that upper bounds the indicator function, i.e., $  P_{\mathbb{I}}(\mathbf{x},\omega):= \sum_{(i,j)} c_{ij} x^{i} \omega^{j} \geq \mathbb{I}_{\mathcal{P} \geq 0}$,
then we can apply the expectation w.r.t. the probability density function of $\mathbf{\omega}$ to the both sides and describe the upper bound probability in terms of the known moments of $\mathbf{\omega}$ as follows:

\begin{align}\label{risk_ind}\mathbb{E}[P_{\mathbb{I}}(\mathbf{x},\omega)]=\sum_{(i,j)} c_{ij} x^{i} \mathbb{E}[\omega^{j}] \geq \mathbb{E}[\mathbb{I}_{\mathcal{P} \geq 0}]=\hbox{Prob}( \mathcal{P}(\mathbf{x},\omega)\geq 0 ) \end{align}
where $\mathbb{E}[\omega^j]$ is the moment of order $j$ of random vector $\omega$ defined in Section <ref>.
Hence, we can construct an inner approximation of the $\Delta$-risk contour in (<ref>), denoted by $\hat{\mathcal{C}}^{\Delta}_{r}$, using the upper bound probability in (<ref>) as follows:

\begin{equation}\label{rc_s_inner}\hat{\mathcal{C}}^{\Delta}_{r}= \{\ \mathbf{x} \in \mathcal{X}: \ P(\mathbf{x})
\leq \Delta  \}
\end{equation}

where polynomial $P(\mathbf{x})=\mathbb{E}[P_{\mathbb{I}}(\mathbf{x},\omega)]=\sum_{(i,j)} c_{ij} \mathbb{E}[\omega^{j}]x^{i} $ as shown in (<ref>).

According to (<ref>) and (<ref>), the problem of constructing $\Delta$-risk contour (<ref>), reduces to the problem of finding an upper bound probability and an upper bound polynomial indicator function of the superlevel set of $\mathcal{P}(\mathbf{x},\omega)$. In [29], to compute the upper bound polynomial indicator functions, we provide an $(n_x+n_{\omega})$-dimensional convex optimization problem in the form of a semidefinite program. Such optimization is not suitable  for  online  computations and is limited  to small dimensions $(n_x+n_{\omega})$. In this paper, we propose an optimization-free approach to obtain upper bound probability and polynomial indicator function that is suitable for online computations and large scale problems as follows:

Let $\mathcal{X}_{obs}(\omega)= \{ \mathbf{x}\in \mathcal{X}: \mathcal{P}(\mathbf{x},\mathbf{\omega}) \geq 0  \}$ be the given static uncertain obstacle as defined in (<ref>). To obtain an upper bound of the probability  $\hbox{Prob}( \mathbf{x} \in \mathcal{X}_{obs}(\omega))=\hbox{Prob}( \mathcal{P}(\mathbf{x},\omega)\geq 0 )$, we begin by defining a new random variable $z \in \mathbb{R}$ as follows:

\begin{equation}
\end{equation}

Note that $z \in \mathbb{R}$ is a random variable, while $\mathcal{P}(\mathbf{x},\omega): \mathbb{R}^{n_x+n_{\omega}} \rightarrow \mathbb{R}$ is a polynomial in $\mathbf{x} \in \mathbb{R}^{n_x}$ and random vector $\mathbf{\omega} \in \mathbb{R}^{n_{\omega}}$. By doing so, we can transform the $(n_x+n_{\omega})$-dimensional probability assessment problem into a one-dimensional probability assessment problem in terms of the new defined random variable $z$, i.e., $\hbox{Prob}( \mathbf{x} \in \mathcal{X}_{obs}(\omega))=\hbox{Prob}( \mathcal{P}(\mathbf{x},\omega)\geq 0 )=\hbox{Prob}( z \geq 0 )$, [30]. Note that the statistics of the random variable $z$, e.g., moments, are functions of $\mathbf{x}$ and the statistics of the random vector $\omega$.

Now, to compute an upper bound of the probability $\hbox{Prob}( z \geq 0 )$, we just need an upper bound polynomial description of one-dimensional indicator function $\mathbb{I}_{z\geq0}$ defined as $\mathbb{I}_{z\geq0} = 1 $ if $z\geq 0$, and
0 otherwise. In this paper, we will use the upper bound polynomial indicator function and the upper bound probability provided by
Cantelli's inequality defined for scalar random variables as $\hbox{Prob}(z \geq 0) \leq \frac{\Exp[z^2] - \Exp[z]^2}{\Exp[z^2]}$ whenever $\Exp[z] \leq 0 $. In other words, the probability of collision, i.e., $\hbox{Prob}(z\geq0)=\hbox{Prob}( \mathcal{P}(\mathbf{x},\omega)\geq 0 )$, is bounded if the expected value of remaining safe is nonnegative, i.e, $\Exp[z]=\Exp[\mathcal{P}(\mathbf{x},\omega)] \leq 0$. For other different one-dimensional indicator function-based probability bounds see [30, 33, 34].

Hence, the upper bound of the probability of collision can be described in terms of the polynomial of the uncertain obstacle as follows:

\begin{align}\label{cheby2}
\hbox{Prob}( \mathbf{x} \in \mathcal{X}_{obs}(\omega)) \leq \frac{\Exp[\mathcal{P}^2(\mathbf{x},\omega)] - \Exp[\mathcal{P}(\mathbf{x},\omega)]^2}{\Exp[\mathcal{P}^2(\mathbf{x},\omega)]}
\end{align}

whenever $\Exp[\mathcal{P}(\mathbf{x},\omega)]\leq 0$.
This will results in an inner approximation of the $\Delta$-risk contour as in (<ref>). Note that although the standard Cantelli's inequality uses the first two moments of the scalar random variable $z$, we need higher order moments of random vector $\omega$ to construct the set in (<ref>). In fact, in this paper, we generalize the standard scalar Cantelli probability bound to obtain multivariate probability bound (<ref>) involving nonconvex and nonlinear sets of obstacles.

[1]
J. Reif and M. Sharir, “Motion planning in the presence of moving obstacles,”
Journal of the ACM, vol. 41, no. 4, pp. 764–790, 1994.

[2]
K. M. Lynch and F. C. Park, Modern Robotics.1em plus 0.5em minus
0.4emCambridge University Press, 2017.

[3]
V. Kumar, M. Zefran, and J. Ostrowski, “Motion planning and control of
robots,” New York: Wiley, 1999.

[4]
R. Charles, A. Bry, and N. Roy, “Polynomial trajectory planning for aggressive
quadrotor flight in dense indoor environments,” Robotics Research,
Springer, pp. vol. 114, 649–666, 2016.

[5]
S. M. LaValle, Planning algorithms.1em plus 0.5em minus
0.4emCambridge university press, 2006.

[6]
R. Deits and R. Tedrake, “Efficient mixed-integer planning for uavs in
cluttered environments,” IEEE International Conference on Robotics and
Automation (ICRA), 2015.

[7]
B. E. Khadir, J. B. Lasserre, and V. Sindhwani, “Piecewise-linear motion
planning amidst static, moving, or morphing obstacles,” IEEE
International Conference on Robotics and Automation (ICRA), 2021.

[8]
L. Blackmore and M. Ono, “Convex chance constrained predictive control without
sampling,” AIAA Guidance, Navigation, and Control Conference, 2009.

[9]
L. Blackmore, M. Ono, A. Bektassov, and B. C. Williams, “A probabilistic
particle-control approximation of chance-constrained stochastic predictive
control,” IEEE Transactions on Robotics, vol. 26, no. 3, pp.
502–517, 2010.

[10]
W. Schwarting, J. Alonso-Mora, L. Pauli, S. Karaman, and D. Rus, “Parallel
autonomy in automated vehicles: Safe motion generation with minimal
intervention,” IEEE International Conference on Robotics and
Automation (ICRA), pp. 1928–1935, 2017.

[11]
B. Luders, M. Kothari, and J. How, “Chance constrained rrt for probabilistic
robustness to environmental uncertainty,” AIAA Guidance, Navigation,
and Control Conference, 2010.

[12]
B. Axelrod, L. P. Kaelbling, and T. Lozano-Pérez, “Provably safe robot
navigation with obstacle uncertainty,” The International Journal of
Robotics Research, vol. 37, no. 13-14, pp. 1760–1774, 2018.

[13]
C. Dawson, A. Jasour, A. Hofmann, and B. Williams, “Provably safe trajectory
optimization in the presence of uncertain convex obstacles,” IEEE
International Conference on Intelligent Robots and Systems (IROS), 2020.

[14]
S. Dai, S. Schaffert, A. Jasour, A. Hofmann, and B. Williams, “Chance
constrained motion planning for high-dimensional robots,” IEEE
International Conference on Robotics and Automation (ICRA), 2019.

[15]
M. Cannon, “Chance-constrained optimization with tight confidence bounds,”
arXiv preprint arXiv:1711.03747, 2017.

[16]
G. C. Calafiore and M. C. Campi, “The scenario approach to robust control
design,” IEEE Transactions on Automatic Control, vol. 51, no. 5, pp.
742–753, 2006.

[17]
L. Janson, E. Schmerling, and M. Pavone, “Monte carlo motion planning for
robot trajectory optimization under uncertainty,” Robotics Research,
pp. 343–361, 2018.

[18]
T. Summers, “Distributionally robust sampling-based motion planning under
uncertainty,” IEEE/RSJ International Conference on Intelligent Robots
and Systems (IROS), pp. 6518–6523, 2018.

[19]
A. Wang, A. Jasour, and B. C. Williams, “Non-gaussian chance-constrained
trajectory planning for autonomous vehicles under agent uncertainty,”
IEEE Robotics and Automation Letters, vol. 5, no. 4, pp. 6041–6048,

[20]
J. Jacod and P. Protter, Probability essentials.1em plus 0.5em
minus 0.4emSpringer Science & Business Media, 2012.

[21]
P. A. Parrilo, “Semidefinite programming relaxations for semialgebraic
problems,” Mathematical programming, vol. 96, no. 2, 2003.

[22]
J. B. Lasserre, “Global optimization with polynomials and the problem of
moments,” SIAM Journal on optimization, vol. 11, no. 3, pp. 796–817,

[23]
M. Laurent, “Sums of squares, moment matrices and optimization over
polynomials,” Emerging Applications of Algebraic Geometry, Springer,
pp. 157–270, 2009.

[24]
J. Lofberg, “Yalmip: A toolbox for modeling and optimization in matlab,”
IEEE International Conference on Robotics and Automation (ICRA), 2004.

[25]
M. M. Tobenkin, F. Permenter, and A. Megretski, “spotless: Polynomial and
conic optimization,” 2013. [Online]. Available:

[26]
D. Henrion, J.-B. Lasserre, and J. Löfberg, “Gloptipoly 3: moments,
optimization and semidefinite programming,” Optimization Methods &
Software, vol. 24, no. 4-5, pp. 761–779, 2009.

[27]
A. A. Ahmadi and B. E. Khadir, “Time-varying semidefinite programs,”
Mathematics of Operations Research, arXiv:1808.03994, 2021.

[28]
A. Bellon, D. Henrion, V. Kungurtsev, and J. Marecek, “Time-varying
semidefinite programming: Geometry of the trajectory of solutions,”
arXiv preprint arXiv:2104.05445, 2021.

[29]
A. Jasour and B. Williams, “Risk contours map for risk bounded motion planning
under perception uncertainties,” Robotics: Science and Systems, 2019.

[30]
A. Jasour, A. Hofmann, and B. C. Williams, “Moment-sum-of-squares approach for
fast risk estimation in uncertain environments,” IEEE Conference on
Decision and Control (CDC), pp. 2445–2451, 2018.

[31]
P. Roux, Y.-L. Voronin, and S. Sankaranarayanan, “Validating numerical
semidefinite programming solvers for polynomial invariants,” Formal
Methods in System Design, vol. 53, no. 2, pp. 286–312, 2018.

[32]
M. Putinar, “Positive polynomials on compact semi-algebraic sets,”
Indiana University Mathematics Journal, vol. 42, no. 3, pp. 969–984,

[33]
A. Jasour, “Risk aware and robust nonlinear planning (rarnop),” Course
Notes for MIT 16.S498, 2019. [Online]. Available: <rarnop.mit.edu>

[34]
A. Nemirovski and A. Shapiro, “Convex approximations of chance constrained
programs,” SIAM Journal on Optimization, vol. 17, no. 4, pp.
969–996, 2007.


# Systems of difference equations, symmetries and integrability conditions

Louis Brady111Email<EMAIL_ADDRESS><EMAIL_ADDRESS>and Pavlos
Xenitidis222Email<EMAIL_ADDRESS>
School of Mathematics, Computer Science & Engineering
Liverpool Hope University, L16 9JD Liverpool, UK

###### Abstract

We consider a class of systems of difference equations defined on an
elementary quadrilateral of the ${\mathbb{Z}}^{2}$ lattice, define their
eliminable and dynamical variables, and demonstrate their use. Using the
existence of infinite hierarchies of symmetries as integrability criterion, we
derive necessary integrability conditions and employ them in the construction
of the lowest order symmetries of a given system. These considerations are
demonstrated with the help of three systems from the class of systems under
consideration.

##  1 Introduction

The theory of symmetries of scalar difference equations is well developed and
several methods and approaches have been proposed for their derivation, see
for instance [X1] and references therein. The existence of an infinite
hierarchy of symmetries serves as an integrability criterion and it provides
us with conservation laws via the so-called integrability conditions. The
derivation of the latter relies on the existence of a formal recursion
operator [MWX, MX], and some of these conditions can be interpreted as
determining equations for the symmetries of the equation under consideration
[X1].

The integrability of a system of difference equations is defined in the same
way using the existence of infinite hierarchies of symmetries as a criterion.
However, there is no systematic method to compute the symmetries of a given
system therefore we have to employ other integrability aspects for their
derivation. For instance there exist examples of integrable systems, the
symmetries of which were derived using either the direct linearisation method
[TN] or a Lax pair [FX]. The purpose of this paper is to contribute in the
development of the theory of integrability conditions for systems of
difference equations and the main aim is to propose a systematic method for
the derivation of the lowest order symmetries of such systems using only the
system under consideration.

The systems of partial difference equations we are going to consider are
defined on an elementary quadrilateral of the ${\mathbb{Z}}^{2}$ lattice and
are determined by certain solvability requirements given in Definition 1. We
call the variables involved in a system eliminable and dynamical when the
former can be expressed uniquely in terms of the latter using the system. The
importance of these variables stems from the fact that they play an important
role in functional equations and the solution of the initial value problem,
which is demonstrated with the use of three illustrative examples.

We define formal Laurent and Taylor series of matrix pseudo-difference
operators and discuss how they can be computed systematically. We employ them
in the derivation of integrability conditions and explain how some of the
derived conditions can be employed as determining equations for the symmetries
of a given system. We demonstrate these interpretations using three systems,
namely the coupled system of discrete potential Korteweg-de Vries (KdV)
equations (13), the discrete nonlinear Schrödinger system (18) [KMX], and the
discrete Toda system (21) [A1].

The paper is organised as follows. The next section introduces our notation,
defines the class of systems we are interested in along with eliminable and
dynamical variables, and contains all the necessary definitions on symmetries.
Section LABEL:sec:recop deals with recursion operators, formal series and
integrability conditions, and discusses the implementation of our theory with
the help of the three aforementioned systems. The concluding section contains
an overall evaluation of our results and a discussion of further developments.

##  2 Systems of difference equations and symmetries

In this section we present our notation and define the class $\cal{Q}$ of
systems of difference equations we are going to consider. We introduce the
notion of eliminable and dynamical variables, and give the basic definitions
of symmetries and conservation laws for difference equations. We present three
systems from the class $\cal{Q}$ to exemplify the concepts we introduce and
use them as running examples throughout the paper.

###  2.1 Notation

In what follows we deal with systems of partial difference equations in two
independent discrete variables denoted by $n$ and $m$. The dependence of a
function $u$ on $n$ and $m$ is denoted with indices omitting the independent
variables, i.e.,

$u(n+i,m+j)=u_{i,j},\qquad\forall~{}i,j\in{\mathbb{Z}}.$

The two shift operators related to the two lattice directions are

${\cal{S}}:n\mapsto n+1,\quad{\cal{T}}:m\mapsto m+1,$

and their action on a function of $n$ and $m$ is defined as

${\cal{S}}^{k}{\cal{T}}^{\ell}(f_{i,j})=f_{i+k,j+\ell}.$

Vectors are written in bold face letters and their components are
distinguished by upper indices in parentheses, i.e.,

${\boldsymbol{f}}_{0,0}=\left(f^{(1)}_{0,0},\ldots,f^{(k)}_{0,0}\right).$

Matrices are denoted with upper case Roman letters. In particular, $\rm{I}$
and $\rm{O}$ denote the identity and the zero matrix, respectively. We also
use the notation ${\boldsymbol{f}}([{\boldsymbol{u}}])$ to denote a vector
function depending on a finite, but otherwise unspecified number of shifted
values of $\boldsymbol{u}$, e.g., a vector function like
${\boldsymbol{f}}({\boldsymbol{u}}_{-1,0},{\boldsymbol{u}}_{0,0},{\boldsymbol{u}}_{1,0})$.

###  2.2 A class of systems of difference equations

We consider systems of $k$ equations for $k$ functions
$u^{(i)}:{\mathbb{Z}}^{2}\mapsto{\mathbb{R}}$, $i=1,\ldots,k$, which involve
the values of these functions at four points which form an elementary
quadrilateral on the $\mathbb{Z}^{2}$ lattice. That is systems of the form

${\boldsymbol{Q}}({\boldsymbol{u}}_{0,0},{\boldsymbol{u}}_{1,0},{\boldsymbol{u}}_{0,1},{\boldsymbol{u}}_{1,1})={\boldsymbol{0}},$
(1)

where
${\boldsymbol{u}}_{i,j}=\left(u^{(1)}_{i,j},\ldots,u^{(k)}_{i,j}\right)$. We
can also write system (1) in component form as

$Q^{(i)}({\boldsymbol{u}}_{0,0},{\boldsymbol{u}}_{1,0},{\boldsymbol{u}}_{0,1},{\boldsymbol{u}}_{1,1})=0,\quad
i=1,\ldots,k.$

We refer to system (1) as _a quad system_.

Quad systems form a very wide class of systems, so we restrict our
considerations to quad systems having certain properties which are described
in the following definition. We denote the class of all these systems with
$\cal{Q}$.

######  Definition 1 (Class of quad systems).

Let $V_{i}$, $i=1,\ldots,4$, denote the four vertices of an elementary
quadrilateral of the $\mathbb{Z}^{2}$ lattice where quad system (1) is
defined, and $E_{ij}$, or $E_{ji}$, denote the undirectional edge connecting
two consecutive vertices $V_{i}$ and $V_{j}$ as shown in Figure 1. Let
$\cal{U}$ denote the set of all variables appearing in quad system (1).

Class $\cal{Q}$ contains quad systems (1) which fulfill the following two
requirements.

  1. (i)

For every ${\boldsymbol{u}}_{p,q}$, $p,q=0,1$, the defining function
$\boldsymbol{Q}$ of system (1) satisfies

${\rm{Q}}_{(p,q)}\neq{\rm{O}}~{}~{}~{}\mbox{for all }~{}p,q=0,1,$ (2)

where ${\rm{Q}}_{(p,q)}$ are the Jacobian matrices of $\boldsymbol{Q}$.

${\rm{Q}}_{(p,q)}=\frac{\partial{\boldsymbol{Q}}}{\partial{\boldsymbol{u}}_{p,q}}=\left(\frac{\partial{\boldsymbol{Q}}}{\partial
u^{(1)}_{p,q}}\cdots\frac{\partial{\boldsymbol{Q}}}{\partial
u^{(k)}_{p,q}}\right),\quad{\mbox{with
}}\left({{\rm{Q}}_{(p,q)}}\right)_{a,b}=\frac{\partial Q^{(a)}}{\partial
u^{(b)}_{p,q}}$ (3)

  2. (ii)

For every edge $E_{ij}$ it can be solved uniquely for $\ell$ variables
$(u^{(a_{1})},\ldots,u^{(a_{\ell})})$ lying on the vertex $V_{i}$ and $k-\ell$
variables $(u^{(b_{1})},\ldots,u^{(b_{k-\ell})})$ lying on vertex $V_{j}$,
where $0\leq\ell\leq k$, and the resulting expressions involve all the
variables in
${\cal{U}}-\\{u^{(a_{1})}_{V_{i}},\ldots,u^{(a_{\ell})}_{V_{i}},u^{(b_{1})}_{V_{j}},\ldots,u^{(b_{k-\ell})}_{V_{j}}\\}$,
where $u^{(a)}_{V_{s}}$ denotes the value of the function $u^{(a)}$ at the
vertex $V_{s}$.

1 0.02 t:F (0 0.5) T (-.25 0) (1.25 0) (0 -.25) (0 1.25) (1 -.25) (1 1.25)
(-.25 1) (1.25 1) (1 0) f:0.0 r:0.07 (0 0) f:0.0 r:0.07 (0 1) f:0.0 r:0.07 (1
1) f:0.0 r:0.07 (-.4 -.27) $(n,m)$ (1.04 -.27) $(n+1,m)$ (-.7 1.1) $(n,m+1)$
(1.05 1.1) $(n+1,m+1)$

1 0.02 t:F (0 0.5) T (-.25 0) (1.25 0) (0 -.25) (0 1.25) (1 -.25) (1 1.25)
(-.25 1) (1.25 1) (1 0) f:0.0 r:0.07 (0 0) f:0.0 r:0.07 (0 1) f:0.0 r:0.07 (1
1) f:0.0 r:0.07 (-.2 -.27) $V_{1}$ (1.04 -.27) $V_{2}$ (-.2 1.1) $V_{4}$ (1.05
1.1) $V_{3}$ (.35 -.27) $E_{12}$ (1.04 .4) $E_{23}$ (-.3 0.4) $E_{41}$ (.35
1.1) $E_{34}$

Figure 1: The elementary quadrilateral of the ${\mathbb{Z}}^{2}$ lattice where
quad system (1) is defined and the naming of its vertices and edges.

######  Definition 2.

We call variables
$\\{u^{(a_{1})}_{V_{i}},\ldots,u^{(a_{\ell})}_{V_{i}},u^{(b_{1})}_{V_{j}},\ldots,u^{(b_{k-\ell})}_{V_{j}}\\}$
_eliminable_ and all the other ones _dynamical_.

We call these variables eliminable because we can eliminate them from
_functional relations and equations_ , like conservation laws and determining
equations for symmetries, which must hold on solutions of the quad system.
They are also related to the initial value problem since the dynamical
variables may be viewed as initial conditions and the eliminable ones as the
values of the solution to be determined in terms of the initial data. These
concepts and interpretations are illustrated below with the help of three
examples.

###  2.3 Symmetries and conservation laws

We define a symmetry of a system of difference equations in the following way.

######  Definition 3.

A vector function

${\boldsymbol{F}}(n,m,[{\boldsymbol{u}}])=\left(F^{(1)}(n,m,[{\boldsymbol{u}}]),\ldots,F^{(k)}(n,m,[{\boldsymbol{u}}])\right)$
(4)

is _a symmetry generator_ , or just _a symmetry_ , of system
${\boldsymbol{Q}}(n,m,[{\boldsymbol{u}}])={\boldsymbol{0}}$ if

$\sum_{p,q}{\rm{Q}}_{(p,q)}{\cal{S}}^{p}{\cal{T}}^{q}({\boldsymbol{F}})={\boldsymbol{0}}~{}{\mbox{
holds on solutions of
}}~{}{\boldsymbol{Q}}(n,m,[{\boldsymbol{u}}])={\boldsymbol{0}},$ (5)

where matrices ${\rm{Q}}_{(p,q)}$ are defined in (3). Equations (5) are
referred to as the determining equations for the symmetries of
${\boldsymbol{Q}}(n,m,[{\boldsymbol{u}}])={\boldsymbol{0}}$.

If the system depends on a parameter $\alpha$ we can also consider symmetries
which act on the parameter.

######  Definition 4.

Let system ${\boldsymbol{Q}}(n,m,[{\boldsymbol{u}}];\alpha)={\boldsymbol{0}}$
depend also on a parameter $\alpha$. Then, a vector function

$\left({\boldsymbol{M}}(n,m,[{\boldsymbol{u}}]),\xi\right)=\left(M^{(1)}(n,m,[{\boldsymbol{u}}]),\ldots,M^{(k)}(n,m,[{\boldsymbol{u}}]),\xi\right)$
(6)

is an _extended symmetry generator_ , or just _an extended symmetry_ , of the
system if

$\sum_{p,q}{\rm{Q}}_{(p,q)}{\cal{S}}^{p}{\cal{T}}^{q}({\boldsymbol{M}})+\xi\partial_{\alpha}{\boldsymbol{Q}}={\boldsymbol{0}}~{}{\mbox{
holds on solutions of
}}~{}{\boldsymbol{Q}}(n,m,[{\boldsymbol{u}}];\alpha)={\boldsymbol{0}}.$ (7)

With every symmetry we can associate a vector field. More precisely, with
symmetry $\boldsymbol{F}$ (4) and extended symmetry $\boldsymbol{M}$ (6) we
can associate the vector fields

$X_{\boldsymbol{F}}=\sum_{i=1}^{k}\sum_{p,q}{\cal{S}}^{p}{\cal{T}}^{q}\left(F^{(i)}\right)\frac{\partial{\phantom{u_{p,q}^{(i)}}}}{\partial
u_{p,q}^{(i)}}\quad{\mbox{and}}\quad
X_{\boldsymbol{M}}=\sum_{i=1}^{k}\sum_{p,q}{\cal{S}}^{p}{\cal{T}}^{q}\left(M^{(i)}\right)\frac{\partial{\phantom{u_{p,q}^{(i)}}}}{\partial
u_{p,q}^{(i)}}+\xi\frac{\partial{\phantom{\alpha}}}{\partial\alpha},$

respectively. The commutator $[~{},~{}]$ of two symmetries $\boldsymbol{F}$
and $\boldsymbol{G}$ results to a new symmetry $\boldsymbol{H}$ which is
defined as

${\boldsymbol{H}}=[{\boldsymbol{F}},{\boldsymbol{G}}]=X_{\boldsymbol{F}}({\boldsymbol{G}})-X_{\boldsymbol{G}}({\boldsymbol{F}}),$
(8a) and its components are given by
$H^{(i)}=\sum_{\ell_{1},\ell_{2}}\sum_{j=1}^{k}{\cal{S}}^{\ell_{1}}{\cal{T}}^{\ell_{2}}\left(F^{(j)}\right)\frac{\partial
G^{(i)}}{\partial
u^{(j)}_{\ell_{1},\ell_{2}}}-{\cal{S}}^{\ell_{1}}{\cal{T}}^{\ell_{2}}\left(G^{(j)}\right)\frac{\partial
F^{(i)}}{\partial u^{(j)}_{\ell_{1},\ell_{2}}}.$ (8b)

If ${\boldsymbol{H}}={\boldsymbol{0}}$ then we say that the two symmetries
commute.

A generalised symmetry is called _master symmetry_ if its commutator with a
generalised symmetry produces another generalised symmetry, and its existence
may be used as an integrability criterion.

######  Remark 1.

Using the relation of symmetries and continuous groups of transformations we
can represent the former as systems of differential-difference equations.
Assuming that $\boldsymbol{u}$ depends on two continuous variables $t$ and
$\tau$ (the group parameters), we can write symmetry generator (4) as

$\frac{\partial{\boldsymbol{u}}_{0,0}}{\partial
t}={\boldsymbol{F}}(n,m,[{\boldsymbol{u}}])\quad\Rightarrow\quad\frac{\partial
u^{(i)}_{0,0}}{\partial
t}=F^{(i)}(n,m,[{\boldsymbol{u}}]),~{}~{}i=1,\ldots,k,$ (9)

and extended symmetry generator (6) as

$\frac{\partial{\boldsymbol{u}}_{0,0}}{\partial\tau}={\boldsymbol{M}}(n,m,[{\boldsymbol{u}}]),~{}~{}\frac{\partial\alpha}{\partial\tau}=\xi\quad\Rightarrow\quad\frac{\partial
u^{(i)}_{0,0}}{\partial\tau}=M^{(i)}(n,m,[{\boldsymbol{u}}]),~{}~{}i=1,\ldots,k,~{}~{}\frac{\partial\alpha}{\partial\tau}=\xi.$
(10)

In what follows we denote symmetries as differential-difference equations.

The sums in all the above definitions are finite. In particular, for quad
systems (1) the determining equations (5) involve only four terms,

${\rm{Q}}_{(0,0)}{\boldsymbol{F}}+{\rm{Q}}_{(1,0)}{\cal{S}}\left({\boldsymbol{F}}\right)+{\rm{Q}}_{(0,1)}{\cal{T}}\left({\boldsymbol{F}}\right)+{\rm{Q}}_{(1,1)}{\cal{ST}}\left({\boldsymbol{F}}\right)={\boldsymbol{0}},$
(11)

and similarly (7) involves only five terms,

${\rm{Q}}_{(0,0)}{\boldsymbol{F}}+{\rm{Q}}_{(1,0)}{\cal{S}}\left({\boldsymbol{F}}\right)+{\rm{Q}}_{(0,1)}{\cal{T}}\left({\boldsymbol{F}}\right)+{\rm{Q}}_{(1,1)}{\cal{ST}}\left({\boldsymbol{F}}\right)+\xi\partial_{\alpha}{\boldsymbol{Q}}={\boldsymbol{0}}.$
(12)

It can also be shown that the symmetries of (1) necessarily have the form

${\boldsymbol{F}}={\boldsymbol{A}}(n,m,{\boldsymbol{u}}_{-M,0},\ldots,{\boldsymbol{u}}_{M,0})+{\boldsymbol{B}}(n,m,{\boldsymbol{u}}_{0,-M},\ldots,{\boldsymbol{u}}_{0,M}),$

where $M$ is a positive integer also known as _order of the symmetry_. This
allows us to consider symmetries involving shifts of $\boldsymbol{u}$ only in
one direction. In what follows we discuss only _first order symmetries in the
first direction_ , i.e. symmetries of the form
${\boldsymbol{F}}={\boldsymbol{F}}(n,m,{\boldsymbol{u}}_{-1,0},{\boldsymbol{u}}_{0,0},{\boldsymbol{u}}_{1,0})$.

######  Definition.

A conservation law of system (5) is a pair of functions
$(\rho,\sigma)=(\rho(n,m,[{\boldsymbol{u}}]),\sigma(n,m,[{\boldsymbol{u}}]))$
such that $({\cal{T}}-1)(\rho)=({\cal{S}}-1)(\sigma)$ on solutions of the
system.
A conservation law is trivial if $\rho=({\cal{S}}-1)(H)$ and
$\sigma=({\cal{T}}-1)(H)$ for some function $H=H(n,m,[{\boldsymbol{u}}])$.

Conservation laws are related to symmetries and integrability as they may be
derived from the integrability conditions as we explain in the next section.

###  2.4 Examples

In order to simplify our notation and avoid the use of many indices in the
three running examples throughout this paper, we use the following notation
for the unknown functions

${\boldsymbol{u}}=(u^{(1)},u^{(2)})=(u,v).$

####  2.4.1 First example

The first example is about the system

$(u_{0,0}-u_{1,1})(v_{1,0}-v_{0,1})-\alpha+\beta=0,\qquad(v_{0,0}-v_{1,1})(u_{1,0}-u_{0,1})-\alpha+\beta=0.$
(13)

It involves all eight components of ${\boldsymbol{u}}$ and its shifts, i.e.
${\cal{U}}=\left\\{u_{0,0},v_{0,0},u_{1,0},v_{1,0},u_{0,1},v_{0,1},u_{1,1},v_{1,1}\right\\}$,
and condition (2) is readily satisfied because none of the Jacobian matrices
are the zero matrix.

${\rm{Q}}_{(0,0)}=-{\rm{Q}}_{(1,1)}=\begin{pmatrix}v_{1,0}-v_{0,1}&0\\\
0&u_{1,0}-u_{0,1}\end{pmatrix},~{}{\rm{Q}}_{(1,0)}=-{\rm{Q}}_{(1,0)}=\begin{pmatrix}0&u_{0,0}-u_{1,1}\\\
v_{0,0}-v_{1,1}&0\end{pmatrix}.$ (14)

For the second requirement in Definition 1, we start by considering edge
$E_{12}$, see Figure 1. It is clear that we can solve system (13) uniquely for
either $(u_{V_{1}},v_{V_{1}})=(u_{0,0},v_{0,0})$ or
$(u_{V_{2}},v_{V_{2}})=(u_{1,0},v_{1,0})$, and in both cases, it is not
difficult to check that the resulting expressions involve all the other
variables appearing in system (13).

In the same way, for edge $E_{23}$ we can solve system (13) uniquely for
either $(u_{V_{2}},v_{V_{2}})=(u_{1,0},v_{1,0})$ or
$(u_{V_{3}},v_{V_{3}})=(u_{1,1},v_{1,1})$. Working similarly with the
remaining two edges we can clearly conclude that system (13) can be solved
uniquely for any pair $(u_{i,j},v_{i,j})$, $i,j=0,1$, and all the requirements
in Definition 1 are fulfilled. Thus, for every edge we have two different
choices for eliminable variables, and the set of all possible pairs is

${\cal{E}}=\\{(u_{0,0},v_{0,0}),~{}(u_{1,0},v_{1,0}),~{}(u_{1,1},v_{1,1}),~{}(u_{0,1},v_{0,1})\\}.$
(15)

This shows that system (13) belongs in class $\cal{Q}$.

Set $\cal{E}$ contains all the possible pairs of variables which we can
replace in a unique way from a given expression using system (13) and it is
useful in computations with symmetries and conservation laws. For instance, to
show that the pair

$(\rho,\sigma)=(u_{0,0}v_{1,0}-u_{1,0}v_{0,0}\,,~{}u_{0,0}v_{0,1}-u_{0,1}v_{0,0})$
(16)

is a conservation law of system (13), it is sufficient to show that
$({\cal{T}}-1)(\rho)-({\cal{S}}-1)(\sigma)=0$, or in explicit form

$u_{0,1}v_{1,1}-u_{1,1}v_{0,1}-u_{0,0}v_{1,0}+u_{1,0}v_{0,0}-u_{1,0}v_{1,1}+u_{1,1}v_{1,0}+u_{0,0}v_{0,1}-u_{0,1}v_{0,0}=0,$
(17)

modulo system (13). First we choose a pair from (15) which we want to
eliminate from (17), e.g. $(u_{1,1},v_{1,1})$, and then rearrange system (13)
for these two variables.

$u_{1,1}=u_{0,0}-\frac{\alpha-\beta}{v_{1,0}-v_{0,1}},~{}~{}~{}v_{1,1}=v_{0,0}-\frac{\alpha-\beta}{u_{1,0}-u_{0,1}}.$

Finally we use these relations to replace $u_{1,1}$, $v_{1,1}$ in (17) and
turn the latter into an identity.

####  2.4.2 Second example

The second example is about a system derived in [KMX] as a discretization of
the nonlinear Schrödinger equation, namely

$u_{1,0}-u_{0,1}-\frac{\alpha-\beta}{1+u_{0,0}v_{1,1}}u_{0,0}=0,\qquad
v_{1,0}-v_{0,1}+\frac{\alpha-\beta}{1+u_{0,0}v_{1,1}}v_{1,1}=0.$ (18)

This system depends only on six of the eight components of ${\boldsymbol{u}}$
and its shifts, as $u_{1,1}$ and $v_{0,0}$ do not appear in (18), hence
${\cal{U}}=\left\\{u_{0,0},u_{1,0},v_{1,0},u_{0,1},v_{0,1},v_{1,1}\right\\}$.
The Jacobian matrices

${\rm{Q}}_{(0,0)}=f_{0,0}\begin{pmatrix}1&0\\\
v_{1,1}^{2}&0\end{pmatrix},\quad{\rm{Q}}_{(1,0)}={\rm{I}},\quad{\rm{Q}}_{(0,1)}=-{\rm{I}},\quad{\rm{Q}}_{(1,1)}=f_{0,0}\begin{pmatrix}0&-u_{0,0}^{2}\\\
0&-1\end{pmatrix},$ (19)

where $f_{0,0}=\frac{\beta-\alpha}{(1+u_{0,0}v_{1,1})^{2}}$, are non-zero,
thus condition (2) is satisfied.

It is clear that this system can be solved uniquely for the vertices
$(u_{V_{2}},v_{V_{2}})=(u_{1,0},v_{1,0})$ and
$(u_{V_{4}},v_{V_{4}})=(u_{0,1},v_{0,1})$. However there are more choices for
eliminable variables. For edge $E_{12}$, we can solve system (18) uniquely
also for $(u_{V_{1}},v_{V_{2}})=(u_{0,0},v_{1,0})$. For edge $E_{23}$, we have
another pair, namely $(u_{V_{2}},v_{V_{3}})=(u_{1,0},v_{1,1})$. Finally, for
edges $E_{34}$ and $E_{41}$ we find $(u_{V_{4}},v_{V_{3}})=(u_{0,1},v_{1,1})$
and $(u_{V_{1}},v_{V_{4}})=(u_{0,0},v_{0,1})$ respectively. Summarising, for
every edge we have two different choices of eliminable variables, and for
every choice the second requirement in Definition 1 is fully satisfied. Thus
the system belongs in class $\cal{Q}$ and the set of all possible eliminable
variables is

${\cal{E}}=\\{(u_{0,0},v_{1,0}),~{}(u_{1,0},v_{1,0}),~{}(u_{1,0},v_{1,1}),~{}(u_{0,1},v_{1,1}),~{}(u_{0,1},v_{0,1}),~{}(u_{0,0},v_{0,1})\\}.$
(20)

The dynamical variables may be viewed as initial conditions and the eliminable
variables as the values of the solution to be determined in terms of the
initial data. For system (18) the initial value problem along a staircase was
discussed in [KMX] which actually only exploits the eliminable variables
$(u_{1,0},v_{1,0})$ and $(u_{0,1},v_{0,1})$. However we can consider other
initial value problems which exploit more than one set of eliminable
variables. In order to simplify our presentation, let indices in $g_{i,j}$
denote the value of function $g$ at the lattice point $(i,j)$, i.e.,
$g_{i,j}=g(i,j)$. Suppose that initial values $\left\\{u(\pm i,0),~{}v(\pm
i,0),~{}v(0,i),~{}u(0,-i)\,;~{}i\in{\mathbb{Z}}_{\geq 0}\right\\}$ are given,
see Figure 2. Then we can determine the solution everywhere on the lattice
using the system and its shifts. More precisely, using system
${\boldsymbol{Q}}_{0,0}=0$ we uniquely determine $(u_{0,1},v_{1,1})$ in terms
of the initial data $u_{0,0}$, $u_{1,0}$, $v_{1,0}$ (black circles) and
$v_{0,1}$ (grey circle). We continue by using systems
${\boldsymbol{Q}}_{j,0}=0$, $j=1,2,\ldots$, and determine successively the
corresponding pairs $(u_{j,1},v_{j+1,1})$. On the other hand, using systems
${\boldsymbol{Q}}_{-j,0}=0$, $j=1,2,\ldots$, we successively determine the
pairs $(u_{-j,1},v_{-j,1})$. In this way we find the solution at $m=1$. We
work similarly for negative values of $m$. Systems
${\boldsymbol{Q}}_{k,-1}=0$, $k=0,1,\ldots$, yield successively the values
$(u_{k+1,-1},v_{k+1,-1})$, whereas systems ${\boldsymbol{Q}}_{-k,-1}=0$,
$k=1,2,\ldots$, yield values $(u_{-k,-1},v_{-k+1,-1})$. In this way we find
the solution at $m=-1$. Obviously we can repeat this procedure and determine
the solution everywhere on the ${\mathbb{Z}}^{2}$ lattice.

.9 0.01 t:F (-2.5 0) (2.5 0) (-2.5 1) (2.5 1) (-2.5 -1) (2.5 -1) (1 -1.5) (1
1.5) (-1 -1.5) (-1 1.5) (2 -1.5) (2 1.5) (-2 -1.5) (-2 1.5) (0 -1.5) (0 1.5)

(0 0) f:0 r:0.06 (1 0) f:0 r:0.06 (-1 0) f:0 r:0.06 (2 0) f:0 r:0.06 (-2 0)
f:0 r:0.06

(0 1) f:0.8 r:0.06 (0 1) r:0.06

(0 -1) f:1 r:0.06 (0 -1) r:0.06

(1 1) r:0.06 (1 -1) r:0.06 (-1 1) r:0.06 (-1 -1) r:0.06 (2 1) r:0.06 (2 -1)
r:0.06 (-2 1) r:0.06 (-2 -1) r:0.06

(0.4 .4) ${\boldsymbol{Q}}_{0,0}$ (-.8 .4) ${\boldsymbol{Q}}_{-1,0}$ (0.2 -.6)
${\boldsymbol{Q}}_{0,-1}$ (-.8 -.6) ${\boldsymbol{Q}}_{-1,-1}$ (1.4 .4)
${\boldsymbol{Q}}_{1,0}$ (-1.8 .4) ${\boldsymbol{Q}}_{-2,0}$ (1.2 -.6)
${\boldsymbol{Q}}_{1,-1}$ (-1.8 -.6) ${\boldsymbol{Q}}_{-2,-1}$

(2.55 0) $m=0$ (2.55 -1) $m=-1$ (2.55 1) $m=1$ (-0.20 -1.65) $n=0$ (0.05 .85)
$u$ (0.85 .85) $v$ (1.05 .85) $u$ (1.85 .85) $v$ (-.95 .80) $(u,v)$ (.55 -.95)
$(u,v)$ (-1.95 .80) $(u,v)$ (1.55 -.95) $(u,v)$ (-.95 -.95) $u$ (-.15 -.95)
$v$ (-1.95 -.95) $u$ (-1.15 -.95) $v$

Figure 2: The initial value problem: initial values (dynamical variables) are
denoted with disks ($u_{\pm i,0}$, $v_{\pm i,0}$ (black disks), $v_{0,i}$
(grey disk), $u_{0,-i}$ (white disk), $i\in{\mathbb{Z}}_{\geq 0}$), and in
every quadrilateral the eliminable variables $u$, $v$ which can be determined
uniquely in terms of the initial data are given next to the corresponding
vertex indicated with a circle.

####  2.4.3 Third example

The last example is about the following quad system

$u_{0,0}-v_{1,1}=0,\qquad(u_{1,0}-u_{0,0})(v_{1,0}-u_{0,0})-(u_{0,1}-u_{0,0})(v_{0,1}-u_{0,0})=0,$
(21)

which upon the elimination of variable $v$ is equivalent to a discrete Toda
system [A1]. It does not depend on $u_{1,1}$, $v_{0,0}$, i.e.
${\cal{U}}=\left\\{u_{0,0},u_{1,0},v_{1,0},u_{0,1},v_{0,1},v_{1,1}\right\\}$,
and all of its Jacobian matrices are different from the zero matrix.

$\displaystyle{\rm{Q}}_{(0,0)}=\begin{pmatrix}1&0\\\
f_{0,0}&0\end{pmatrix},\quad{\rm{Q}}_{(1,0)}=\begin{pmatrix}0&0\\\
v_{1,0}-u_{0,0}&u_{1,0}-u_{0,0}\end{pmatrix},$ (22a)
$\displaystyle{\rm{Q}}_{(0,1)}=\begin{pmatrix}0&0\\\
u_{0,0}-v_{0,1}&u_{0,0}-u_{0,1}\end{pmatrix},\quad{\rm{Q}}_{(1,1)}=\begin{pmatrix}0&-1\\\
0&0\end{pmatrix},$ (22b)

where $f_{0,0}=u_{0,1}+v_{0,1}-u_{1,0}-v_{1,0}$.

It is not difficult to see that there is only one pair of eliminable variables
for every edge, namely

${\cal{E}}=\\{(u_{0,0},v_{1,0}),~{}~{}(u_{1,0},v_{1,1}),~{}~{}(u_{0,1},v_{1,1}),~{}~{}(u_{0,0},v_{0,1})\\},$
(23)

and then check that the conditions of the second requirement in Definition 1
are satisfied. Thus (21) belongs in class $\cal{Q}$.

To prove that

$G^{(1)}_{0,0}=(u_{1,0}-u_{0,0})(v_{1,0}-u_{0,0}),~{}~{}~{}G^{(2)}_{0,0}=(u_{-1,0}-v_{0,0})(v_{-1,0}-v_{0,0})$

is a symmetry of system (21), we have to show that the determining equations

$\displaystyle G^{(1)}_{0,0}-G^{(2)}_{1,1}=0,$ $\displaystyle
f_{0,0}G^{(1)}_{0,0}+(v_{1,0}-u_{0,0})G^{(1)}_{1,0}+(u_{1,0}-u_{0,0})G^{(2)}_{1,0}+(u_{0,0}-v_{0,1})G^{(1)}_{0,1}+(u_{0,0}-u_{0,1})G^{(2)}_{0,1}=0,$

hold on solutions of (21). These equations involve 12 variables, namely
$u_{-1,1}$, $v_{-1,1}$, $u_{0,0}$, $u_{1,0}$, $u_{0,1}$, $u_{1,1}$, $v_{0,0}$,
$v_{0,1}$, $v_{1,0}$, $v_{1,1}$, $u_{2,0}$ and $v_{2,0}$, which are related
via the system ${\boldsymbol{Q}}={\boldsymbol{0}}$ and its shifts
${\cal{S}}({\boldsymbol{Q}})={\boldsymbol{0}}$,
${\cal{S}}^{-1}({\boldsymbol{Q}})={\boldsymbol{0}}$, see Figure
LABEL:fig:sys3. Using these three systems we can eliminate some of the
variables in a consistent way. First we choose a pair of eliminable variables
in (23), for instance $(u_{0,0},v_{1,0})$, which implies that we have also to
eliminate its shifts $(u_{-1,0},v_{0,0})$ and $(u_{1,0},v_{2,0})$. Then we
eliminate $(u_{-1,0},v_{0,0})$ using
${\cal{S}}^{-1}({\boldsymbol{Q}})={\boldsymbol{0}}$, next $(u_{0,0},v_{1,0})$
using ${\boldsymbol{Q}}={\boldsymbol{0}}$, and finally $(u_{1,0},v_{2,0})$
with the use of ${\cal{S}}({\boldsymbol{Q}})={\boldsymbol{0}}$. Upon the
elimination of these variables, the determining equations become identities.

1 0.01 t:F (1 0) (2 0) (2 1) (1 1) (1 0) (0 0) (1 0) (1 1) (0 1) (0 0)