# Chance-Constrained Multi-Robot Motion Planning under
Gaussian Uncertainties

Anne Theurkauf∗, Justin Kottinger∗, Nisar Ahmed, and Morteza Lahijanian This
work was supported by NASA STTR award 80NSSC20C0314.∗A. Theurkauf and J.
Kottinger have equal contributions in this work.Authors are with the
department of Aerospace Engineering Sciences at the University of Colorado
Boulder, CO, USA<EMAIL_ADDRESS>
###### Abstract

We consider a chance-constrained multi-robot motion planning problem in the
presence of Gaussian motion and sensor noise. Our proposed algorithm, CC-K-
CBS, leverages the scalability of kinodynamic conflict-based search (K-CBS) in
conjunction with the efficiency of Gaussian belief trees as used in the
Belief-$\mathcal{A}$ framework, and inherits the completeness guarantees of
Belief-$\mathcal{A}$’s low-level sampling-based planner. We also develop three
different methods for robot-robot probabilistic collision checking, which
trade off computation with accuracy. Our algorithm generates motion plans
driving each robot from its initial to goal state while accounting for
uncertainty evolution with chance-constrained safety guarantees. Benchmarks
compare computation time to conservatism of the collision checkers, in
addition to characterizing the performance of the planner as a whole. Results
show that CC-K-CBS scales up to 30 robots.

##  I INTRODUCTION

Robotic teams with diverse capabilities provide a powerful tool for
efficiently accomplishing complex tasks. Applications extend from planetary
exploration to automated warehouses and search and rescue missions. However,
planning for many robots is a challenging problem. A major difficulty of
multi-robot motion planning (MRMP) lies in the size of the planning space,
which grows exponentially with the number of robots. Planning must also
respect the robots’ kinodynamic constraints and avoid robot-robot and robot-
obstacle collision, a difficult problem even for a single robot. Additional
complexities arise for realistic robots due to uncertainty in both motion and
sensing, which the planner needs to account for in collision checking. In this
work, we focus on the uncertain MRMP problem, and develop a scalable algorithm
that guarantees collision avoidance and reaching goal with completeness
properties.

Extensive work has considered MRMP problem in deterministic settings. Proposed
methods are either centralized (coupled) [1, 2, 3, 4, 5] or decentralized
(decoupled) [6, 7, 8, 9]. Centralized planners compose all robots’ states
together into a single meta-robot state, to which single-robot planning
methods with completeness guarantees can then be applied. However, the
exponential increase in the state-space dimension makes centralized approaches
scale poorly. Decentralized methods address this by planning for each robot
individually and considering system collisions separately. These algorithms
drastically cut computation time but often sacrifice completeness guarantees.
For example, [10] plans for one robot at a time, where each robot is
constrained by the paths of all previously-planned robots. While fast, this
method is provably incomplete [11]. One recent contribution to scalable MRMP
is Kinodynamic Conflict-Based Search (K-CBS), [12]. K-CBS pairs a low-level
kinodynamic motion planner for individual robots with a high-level tree search
over the space of possible plans. K-CBS offers scalability and includes a
merging technique to uphold probabilistic completeness properties of the low-
level motion planner. Although a great step forward, none of these approaches
are robust to state and measurement uncertainty, which are inherent to
realistic robotic systems.

Figure 1: Sample plans for 30 robots with unicycle dynamics under motion and
sensing uncertainties. Ellipses are $95\%$ safety contours. Star-centered
circles are goal regions.

Motion planning under Gaussian uncertainty has largely been studied for
single-robot systems. The most efficient algorithms are sampling-based
planners, which are analogous to their deterministic counterparts, e.g., the
feedback-based information roadmap [13] and rapidly-exploring random belief
trees [14]. The belief-$\mathcal{A}$ framework introduced in [15] provides a
generic structure for adapting any sampling-based kinodynamic planner
$\mathcal{A}$ into a belief space planner. It shows that fast planning is
possible by using chance constraint formulation of collision avoidance.
Although these approaches can solve uncertain motion planning, they cannot
handle uncertain MRMP out of the box.

One approach to uncertain MRMP is to abstract the problem to a multi-agent
Markov Decision Process (MMDP), which can be reduced to single-agent MDPs with
a joint action space [16]. This relies on abstraction of the robot dynamics to
a discretized transition probability matrix, which is difficult and scales
poorly to many robots. Additionally, MMDPs solely account for uncertain
dynamics, whereas measurement noise requires using a Partially Observable MDP
(POMDP), which can be very difficult to solve and also scales poorly. Other
uncertain MRMP approaches use _online_ distributed planning, e.g., [17, 18,
19, 20], where each robot plans over a short horizon to account for nearby
robots and obstacles. These methods perform well in unknown environments.
However, neither POMDP approaches nor online distributed methods provide any
formal performance or completeness guarantees.

We address this gap by presenting a scalable MRMP algorithm for robots
operating with Gaussian state and measurement noise. We combine the
scalability of K-CBS with the fast belief-space planning of
belief-$\mathcal{A}$ to solve MRMP for robots under Gaussian uncertainty. Our
algorithm, Chance-Constrained K-CBS (CC-K-CBS), is probabilistically complete
via the completeness inheritance properties of K-CBS and belief-$\mathcal{A}$
[12, 15].

We present three main contributions: (i) three probabilistic robot-robot
collision checking algorithms of varying conservatism and computational
complexity; (ii) a scalable decentralized planner that uses novelties from (i)
to solve MRMP for robots operating with Gaussian noise; and (iii) case studies
and benchmarks demonstrating the applicability of the larger MRMP framework,
validity of our collision-checking methods and comparisons between the
algorithms.

##  II PROBLEM FORMULATION

###  II-A Uncertain Multi-Robot System

We consider $N_{A}$ uncertain robots in a shared workspace
$\mathcal{W}\subset\mathbb{R}^{w}$, $w\in\\{2,3\\}$ with a set of static
obstacles $\mathcal{W}_{O}\subset\mathcal{W}$. Each robot
$i\in\\{1,\ldots,N_{A}\\}$ has a rigid body
$\mathcal{B}^{i}\subset\mathcal{W}$ and uncertain dynamics

$x_{k+1}^{i}=A^{i}x_{k}^{i}+B^{i}u_{k}^{i}+w^{i}_{k},\quad
w_{k}^{i}\sim\mathcal{N}(0,Q^{i}),$ (1)

where $x_{k}^{i}\in\mathcal{X}^{i}\subseteq\mathbb{R}^{n_{i}}$ and
$u_{k}^{i}\in\mathcal{U}^{i}\subseteq\mathbb{R}^{p_{i}}$ are the state and
controls, respectively, with the associated $A^{i}\in\mathbb{R}^{n_{i}\times
n_{i}}$ and $B^{i}\in\mathbb{R}^{n_{i}\times p_{i}}$ matrices, and $w^{i}_{k}$
is the motion (process) noise sampled from a Gaussian distribution with zero
mean and covariance $Q^{i}\in\mathbb{R}^{n_{i}\times n_{i}}$.

Note that state $x_{k}^{i}\in\mathcal{X}^{i}$ corresponds to the origin of the
body frame of robot $\mathcal{B}^{i}$ w.r.t. to the global frame. The robot
body $\mathcal{B}^{i}$ is defined as a set of points in its local frame. Since
state $x_{k}^{i}$ includes configuration (position and orientation) as well as
velocity terms, the body of the robot is well-defined w.r.t. global frame by
projection of the state to either the configuration space (C-space) or
workspace ($\mathcal{W}$). We define operator
$\textsc{proj}_{\mathcal{S}}:\cup_{i=1}^{N_{A}}\mathcal{X}^{i}\to\mathcal{S}$
to be this projection of $x_{k}^{i}\in\mathcal{X}^{i}$ to its the lower
dimensional space $\mathcal{S}\subset\mathcal{X}^{i}$, e.g.,
$\textsc{proj}_{\mathcal{W}}(x_{k}^{i})\subset\mathcal{W}$ is the set of
points ($\mathcal{B}^{i}$) that robot $i$ at state $x_{k}^{i}$ occupies in the
workspace.

Each robot $i\in\\{1,\dots,N_{A}\\}$ also has a sensor (measurement) model
described by

$y_{k}^{i}=C^{i}x_{k}^{i}+v^{i}_{k},\quad v_{k}^{i}\sim\mathcal{N}(0,R^{i}),$
(2)

where $y_{k}^{i}\in\mathcal{Y}^{i}\subseteq\mathbb{R}^{m_{i}}$ is the sensor
measurement, $C^{i}\in\mathbb{R}^{m_{i}\times n_{i}}$ is the corresponding
matrix, and $v_{k}^{i}$ is the measurement noise, sampled from Gaussian
distribution with zero mean and covariance $R^{i}$. We additionally assume the
system is both controllable and observable.

###  II-B Estimation, Control, and Motion Plan

Each robot is equipped with a standard Kalman Filer (KF) and feedback
controller. The uncertainty of robot $i$ state, also known as the belief of
$x^{i}_{k}$ and denoted by $b(x^{i}_{k})$, is given by the KF as a Gaussian
distribution at each time step $k$, i.e., $x^{i}_{k}\sim
b(x_{k}^{i})=\mathcal{N}(\hat{x}_{k}^{i},\Sigma_{k}^{i})$, where
$\hat{x}_{k}^{i}\in\mathcal{X}^{i}$ is the state estimate (mean), and
$\Sigma_{k}^{i}\in\mathbb{R}^{n_{i}\times n_{i}}$ is the covariance matrix.

We define _motion plan_ for robot $i$ to be a pair
$(\check{U}^{i},\check{X}^{i})$ of nominal controls
$\check{U}^{i}=(\check{u}_{0}^{i},\check{u}_{1}^{i},\ldots,\check{u}_{T-1}^{i})$
and the corresponding trajectory of nominal states
$\check{X}^{i}=(\check{x}_{0}^{i},\check{x}_{1}^{i},\ldots,\check{x}_{T}^{i})$,
where $\check{X}^{i}$ is obtained via propagation of the dynamics in (1) under
$\check{U}^{i}$ in the absence of process noise. Then, robot $i$ executes
motion plan $(\check{U}^{i},\check{X^{i}})$ by using a proportional
stabilizing feedback controller
$u_{k}^{i}=\check{u}_{k}^{i}-K(\hat{x}_{k}^{i}-\check{x}_{k}^{i}),$ where $K$
is the control gain.

###  II-C Probabilistic Collisions and Task Completion

Each robot $i$ is assigned a goal region
$\mathcal{X}_{G}^{i}\subset\mathcal{X}^{i}$ in its state space. The workspace
contains $|\mathcal{W}_{O}|=N_{O}$ static obstacles and $N_{A}$ robots (moving
obstacles). The motion planning task is to compute a plan for every robot to
reach its goal and avoid collisions with all (static and moving) obstacles. In
a deterministic setting, these objectives (reaching goal and obstacle
avoidance) are binary. In a stochastic setting, they are probabilistic, as
defined below.

The probability that robot $i$ is in the goal region at time $k$ is given by

$P^{i_{k}}_{G}=P(x^{i}_{k}\in\mathcal{X}^{i}_{G})=\int_{\mathcal{X}^{i}_{G}}b(x_{k}^{i})(s)ds,$
(3)

where $b(x^{i}_{k})(s)$ is the distribution $b(x^{i}_{k})$ evaluated at state
$s\in\mathcal{X}^{i}$. For collision probability, we consider static and
moving obstacles separately. Let $\mathcal{X}^{i}_{O}\subseteq\mathcal{X}^{i}$
be the set of static obstacles in the robot $i$’s state space. It is
constructed from $\mathcal{W}_{O}$ (and the limits on the velocity terms).
Then, similar to probability of goal, the probability of collision of robot
$i$ at time step $k$ with the static obstacles is

$P^{i_{k}}_{O}=P(x^{i}_{k}\in\mathcal{X}^{i}_{O})=\int_{\mathcal{X}^{i}_{O}}b(x_{k}^{i})(s)ds.$
(4)

The probability of one robot colliding with another is more difficult to
define because we must account for uncertainty in both robot bodies. To
describe the collision states, it is best to consider the composed space of
the two robots. That is, for robot $i$ and $j$, $i\neq j$, define the composed
state as
$x^{ij}_{k}=(x^{i}_{k},x^{j}_{k})\in\mathcal{X}^{i}\times\mathcal{X}^{j}=\mathcal{X}^{ij}$.
The set of points that correspond to the two robots colliding is described by

$\displaystyle\mathcal{X}^{ij}_{\text{coll}}=\\{(x^{i}_{k},x^{j}_{k})\in\mathcal{X}^{ij}\mid\textsc{proj}_{\mathcal{W}}(x^{i}_{k})\cap\textsc{proj}_{\mathcal{W}}(x^{j}_{k})\neq\emptyset\\}.$

The probability of collision is then given by

$\displaystyle P^{ij_{k}}_{\text{coll}}$
$\displaystyle=P((x^{i}_{k},x^{j}_{k})\in\mathcal{X}^{ij}_{\text{coll}})=\int_{\mathcal{X}^{ij}_{\text{coll}}}b(x_{k}^{ij})(s)ds$
$\displaystyle=\int_{(s^{i},s^{j})\in\mathcal{X}^{ij}_{\text{coll}}}b(x_{k}^{i})(s^{i})b(x_{k}^{j})(s^{j})ds^{i}ds^{j},$
(5)

where the equality in the second line is due to independence of the robots’
uncertainties from each other.

In this work, we use a chance constraint to bound the safety (violating)
probabilities. That is, given a minimum probability of safety
$p_{\text{safe}}$, we want to ensure that probability of collision (with
static and moving obstacles) at every time step is less than
$1-p_{\text{safe}}$, and the probability of ending in the goal is at least
$p_{\text{safe}}$.

###  II-D Chance-Constrained MRMP Problem

Given $N_{A}$ robots with uncertain dynamics in (1), noisy measurements in
(2), and controller from Sec. II-B in workspace $\mathcal{W}$ with static
obstacles $\mathcal{W}_{O}$, a set of initial distributions
$\\{x_{0}^{i}=\mathcal{N}(\hat{x}_{0}^{i},\Sigma^{i}_{0})\\}_{i=1}^{N_{A}}$, a
set of goal regions $\\{\mathcal{X}_{G}^{i}\\}_{i=1}^{N_{A}}$, and a safety
probability threshold $p_{\text{safe}}$, for each robot
$i\in\\{1,\ldots,N_{A}\\}$, compute a motion plan
$(\check{U}^{i},\check{X}^{i})=((\check{u}_{0}^{i},\check{u}_{1}^{i},\ldots,\check{u}_{T-1}^{i}),(\check{x}_{0}^{i},\check{x}_{1}^{i},\ldots,\check{x}_{T}^{i}))$
such that $\check{x}_{0}^{i}=\hat{x}_{0}^{i}$,

$\displaystyle P^{i_{k}}_{O}+\sum_{j=1,j\neq
i}^{N_{A}}P^{ij_{k}}_{\text{coll}}\leq 1-p_{\text{safe}}$
$\displaystyle\forall k\in\\{1,..,T\\}$ (6a) $\displaystyle P^{i_{T}}_{G}\geq
p_{\text{safe}}$ (6b)

where probability terms $P^{i_{k}}_{O}$, $P^{ij_{k}}_{\text{coll}}$, and
$P^{i_{T}}_{G}$ are given in (4), (5), and (3), respectively.

Note that, unlike most related work, e.g., [17, 18, 19, 20], that solve the
uncertain MRMP problem using _online_ (iterative) approaches with finite
horizon predictions, the above problem asks for the entire plan with safety
guarantees. This implies that the plans must be computed _offline_ while the
chance-constraint guarantees must hold at run time. The advantage to this
approach is that the computed plans are ensured to be sound, i.e., hold hard
guarantees on safety and reachability, before deployment, whereas the online
methods may produce incorrect solutions, e.g., run into a deadlock. In this
work, we seek an offline MRMP algorithm that guarantees soundness _and_
exhibits probabilistic completeness, i.e., generates a valid plan almost-
surely if one exists.

##  III Chance-Constrained MRMP Algorithms

We now describe two algorithms for Chance-Constrained MRMP (CC-MRMP). We begin
by formalizing the composed belief space for $N_{A}$ robots with models in (1)
and (2). We then propose a centralized approach to solving CC-MRMP by
employing an existing single-robot belief-space planning algorithm called
_Belief-_ $\mathcal{A}$ [15]. To accommodate larger multi-robot systems, we
introduce a decentralized algorithm called _CC-K-CBS_ which adapts the novel
K-CBS algorithm to plan for systems with state and measurement uncertainty. We
describe the empirical advantages and disadvantages of these methods in Sec.
V.

###  III-A Centralized Chance-Constrained MRMP

Centralized methods for uncertain MRMP require composing the system beliefs
into a single _meta-belief_ , $\textbf{B}_{k}$, with a single dynamical
constraint. We can then leverage Belief-$\mathcal{A}$ on this meta-belief.
Belief-$\mathcal{A}$ provides a framework for adapting any kinodynamic
sampling-based planner $\mathcal{A}$ to efficiently plan in the Gaussian
belief space under chance constraints as described by (6a) and (6b). See [15]
for details.

For centralized CC-MRMP, because all the robots are assumed to be independent,
sampling and propagation of the meta-belief can be performed independently for
each robot, and the meta-agent belief reconstructed from the individual
robots’ _expected beliefs_ , $\textbf{b}(x_{k}^{i})$. These individual beliefs
are obtained through [14], which provides a method for predicting the belief
over states for a linearizable and controllable system that uses a KF for
estimation. This method forecasts the belief with a priori unknown
measurements, yielding an expected belief defined as:

$\displaystyle\textbf{b}(x_{k}^{i})$
$\displaystyle=\mathbb{E}_{Y}[b(x_{k}^{i}|x_{0}^{i},y_{0:k}^{i})]=\int_{y_{0:k}^{i}}b(x_{k}^{i}|x_{0}^{i},y_{0:k}^{i})pr(y_{0:k}^{i})dy.$

We use this expected belief for offline motion planning and evaluation of
chance constraints in (6a) and (6b). For a given initial belief $b(x_{0}^{i})$
and nominal trajectory, $\check{X}^{i}$, the expected belief is
$\textbf{b}(x_{k}^{i})=\mathcal{N}(\check{x}_{k}^{i},\Gamma_{k}^{i})$, where
$\Gamma_{k}^{i}=\Sigma_{k}^{i}+\Lambda_{k}^{i}$ can be recursively calculated
as derived in [14]:

$\displaystyle\Sigma_{k}^{i-}$
$\displaystyle=A^{i}\Sigma^{i}_{k-1}A^{iT}+Q^{i},\quad\Sigma_{k}^{i}=\Sigma_{k}^{i-}-L_{k}^{i}C^{i}\Sigma_{k}^{i-},$
(7) $\displaystyle\Lambda_{k}^{i}$
$\displaystyle=(A^{i}-B^{i}K^{i})\Lambda_{k-1}^{i}(A^{i}-B^{i}K^{i})^{T}+L_{k}^{i}C^{i}\Sigma_{k}^{i-},$
(8)

where $\Sigma_{k}^{i}$ is the online KF uncertainty, and $\Lambda_{k}^{i}$ is
covariance of the expected state estimates $\hat{x}_{k}^{i}$. Intuitively,
this distribution can be thought of as the sum of the online estimation error
and the uncertainty from the a priori unknown measurements made during
execution. Given the expected beliefs for the agents, the meta-belief can be
constructed such that
$\textbf{B}_{k}=\mathcal{N}(\check{X}_{k}^{M},\Gamma_{k}^{M})$, where
$\check{X}_{k}^{M}=(\check{x}_{k}^{iM},...,\check{x}_{k}^{N_{A}M})$, and
$\Gamma_{k}^{M}$ is constructed as a block diagonal with elements
$\Gamma_{k}^{i},i=1,...,N_{A}$. Similarly, the validity checkers can be
iterated for each robot, such that an entire meta-belief node is considered
invalid if a single robot’s belief is deemed invalid.

This centralized algorithm is the simplest of our two approaches, however it
suffers from poor scalability with the number of agents, just as centralized
deterministic planners do. To this end, we introduce a decentralized planner.

###  III-B Decentralized Chance-Constrained K-CBS

Input: $\mathcal{W}$, $\\{x_{0}^{i}\\}_{i=1}^{N_{A}}$, N, B, $p_{\text{safe}}$

Output: Valid motion plans $\\{(\check{U}^{i},\check{X}^{i})\\}_{i=1}^{N_{A}}$

1 $Q,n_{0},\leftarrow\emptyset$; $n_{0}$.plan $\leftarrow$
getInitPlan$()$;$Q$.add($n_{0}$);

2 while  _solution not found_ do

3 $c\leftarrow$ $Q$.top();

4 if  _$c$.plan is empty _ then

5 $Q$.pop(); replan(N, $c$); $Q$.add($c$);

6

7 else

8 K $\leftarrow$ validatePlan($c$.plan, $p_{\text{safe}}$);

9 if  _K is empty_ then

10 return _$c$.plan_

11 else if  _shouldMerge(K, B)_ then

12 return _mrg &rstrt($\mathcal{W}$, $\\{x_{0}^{i}\\}_{i=1}^{N_{A}}$, N, B,
K)_

13 else

14 $Q$.pop();

15 for  _every robot $i$ in K_ do

16 $c_{new}.\mathcal{C}\leftarrow$ $c.\mathcal{C}\>\cup$ createCSTR(K, $i$);

17 replan(N, $c_{new}$); $Q$.add($c_{new}$);

18

19

Algorithm 1 Chance Constrained K-CBS

Our decentralized CC-MRMP algorithm is based on Kinodynamic Conflict-Based
Search (K-CBS) [12]. K-CBS consists of a high-level conflict-tree search and a
low-level motion planner. At the high-level, K-CBS tracks a constraint-tree
(CT), in which each node represents a suggested plan, which might contain
conflicts (i.e., collisions). At each iteration, K-CBS picks an unexplored
node from the CT based on some heuristic and identifies the conflicts in the
node’s corresponding plan. Then, K-CBS tries to resolve the conflicts by
creating child nodes as follows: if robots $i$ and $j$ are in collision for
some time duration $[t_{s},t_{e}]$ (denoted $\langle
i,j,[t_{s},t_{e}]\rangle$), then two child nodes are created, one with the
constraint that $\mathcal{B}^{i}$ cannot intersect with $\mathcal{B}^{j}$ for
all $t\in[t_{s},t_{e}]$ and the other dually for robot $j$. Then, in each
child node, a low-level sampling-based planner attempts to replan a trajectory
that satisfies the set of all constraints in that branch of the CT within
$N>0$ iterations. This process repeats until a non-colliding plan is found. In
the case where greater than $B>0$ conflicts arise, K-CBS includes a merge and
restart procedure that recursively combines a subset of the robots. We refer
the reader to [12] for details on the original K-CBS. We now describe CC-K-CBS
(Alg. 1), which requires several adaptations of K-CBS (underlined in Alg. 1).

Firstly, note that in K-CBS, collision checking between robots and static
obstacles are performed separately. The low-level planner keeps track of
workspace static obstacles, and the CT node accounts for robot-robot
collisions at the high-level. This separation is key for scalability, but it
introduces a major challenge in chance-constrained planning for uncertain
robots. We expand on this further in Sec. IV. Secondly, since states are now
beliefs, CC-K-CBS requires a Gaussian belief low-level motion planner. We
propose to use Belief-$\mathcal{A}$ for the replan and getInitPlan procedures.
Thirdly, CC-K-CBS requires a definition of constraints; we introduce _belief-
constraints_. This simply replaces the deterministic shape with a non-
deterministic belief. Thus, given a conflict $K=\langle
i,j,[t_{s},t_{e}]\rangle$, one constraint forces robot $i$ to satisfy the
chance-constraint with $b(x_{k}^{j})$ for all $t\in[t_{s},t_{e}]$ and dually
for robot $j$. This enables createCSTR (Line 1) to define meaningful
constraints during replanning with Belief-$\mathcal{A}$.

Lastly, validatePlan (Line 1) requires an efficient way of validating the
chance constraints (Equations 6a and 6b) while simulating prospective MRMP
plans. We propose several techniques for this computation, each of which is
described in detail in Sec. IV. Every method requires comparison with the
maximum allowable probability of collision $p_{coll}=1-p_{\text{safe}}$. We
say validatePlan procedure is _complete_ if it can correctly evaluate the
constraints (6a) and (6b) for all Gaussian beliefs.

The following theorem provides the sufficient conditions for CC-K-CBS to be
probabilistically complete.

######  Theorem 1.

CC-K-CBS is probabilistically complete if the underlying planner
Belief-$\mathcal{A}$ is probabilistically complete and the validatePlan
procedure is complete.

###### Proof.

The proof builds on Theorem 1 in [12], which shows that K-CBS inherits the
probabilistic completeness property of the low-level planner. Thus,
probabilistic completeness of Belief-$\mathcal{A}$ is required for CC-K-CBS.
The completeness of validatePlan ensures that CC-K-CBS does not reject valid
plans. It follows that createCSTR does not pass spurious constraints to the
low-level planner. Therefore, under these conditions, the probabilistic
completeness proof for K-CBS extends to CC-K-CBS. ∎

The completeness condition on validatePlan can be satisfied by exact
evaluation of the probabilities in (3)-(5). However, this operation can be
computationally expensive. To achieve tractability and scalability, we
introduce conservatism into validatePlan by over-approximating the probability
of collision. With this validatePlan, we can still guarantee soundness of the
CC-K-CBS solutions, but completeness is affected since valid plans may be
rejected.

######  Corollary 1.

If validatePlan uses an over-approximation of the probability of collision,
then CC-K-CBS is probabilistically complete up to the accuracy of the over-
approximation.

Below, we introduce several methods for over-approximating the probability of
collision.

##  IV Robot-Robot Collision Checking Algorithms

Recall that the probability of collision between robots $i$ and $j$, as
defined in (5), requires integration of the joint probability distribution of
the two robots over set $\mathcal{X}^{ij}_{\text{coll}}$. This computation is
difficult because (i) $\mathcal{X}^{ij}_{\text{coll}}$ is hard to construct,
and (ii) integration of the joint probability distribution function is
expensive. However, the constraint-checking procedure must be extremely fast
in sampling-based algorithms because it is called in every iteration of tree
extension. In this section, we introduce three approximation methods for this
integration that trade off accuracy with computation effort.

(a)

(b)

Figure 2: (a) 3D view of the c-space obstacle $\mathcal{X}_{coll}^{ij}$ for
robot $i$ with fixed robot $j$ calculated at discrete choices of
$\theta_{k}^{i}$. (b) Top-down view of $\mathcal{X}_{coll}^{ij}$. The
projection of $\mathcal{X}_{coll}^{ij}$ onto x-y plane is bounded by a green
circle.

For clarity of explanation, consider two deterministic (rotating and
translating) robots $i$ and $j$ in a 2D workspace, with configurations
$(\text{x}_{k}^{l},\text{y}_{k}^{l},\theta_{k}^{l})$, $l\in\\{i,j\\}$. We can
construct the 4D c-space obstacles for this pair by taking the Minkowski
difference of the two agents for every pair of orientation angles (see the
c-space construction in Fig. 2(a), which shows the c-space in 3D for a given
orientation of agent $j$). This 4D obstacle is generally difficult to
construct, so deterministic collision checkers typically resort to collision
checking in the workspace, e.g., by bounding the projection of this set, shown
by a green disk in Fig. 2(b).

Now, consider that each robot is described by a random variable (RV), so the
position and orientation of their bodies are no longer deterministic but given
by beliefs, such that each point in the 4D c-space maps to a probability. The
exact probability of collision is the integral of the probability mass over
the c-space obstacle. This results in the two challenges listed above: c-space
obstacle construction and integral computation. Each method proposed below
attempts to address these core problems.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 3:  (a) Circle bound on robots’ bodies. (b) $p_{\text{safe}}$ contour
and corresponding circular over bound on $p_{\text{safe}}$. (c) Inflated
circular bounds intersect, indicating possible collision between robots. (d)
Integration area for probability of collision with circular disk bound on
robots’ bodies. (e) Polygonal bound on circular integration area. (f)
Transformed distribution and integration area, with over-approximated grid.

###  IV-A Method 1 (M.1): Safety Contour

M.1 is inspired by deterministic collision checking methods, where collisions
are detected by intersections of robot bodies in the workspace. We adapt this
framework to uncertain collision checking by inflating bounds on the robot
bodies to encapsulate $p_{\text{safe}}$ probability mass.

Recall that the expected belief is
$\textbf{b}(x_{k}^{i})=\mathcal{N}(\check{x}_{k}^{i},\Gamma_{k}^{i})$. We
define the Gaussian marginal of this distribution in the workspace as
$\textbf{b}_{\mathcal{W}}(x_{k}^{i})$. Let $\bar{x}^{i}_{k}\in\mathcal{W}$ be
the position components of $x_{k}^{i}$, and
$\bar{\Gamma}\in\mathbb{R}^{w\times w}$ be the corresponding covariance matrix
(sub-matrix of $\Gamma^{i}_{k}$). Then
$\textbf{b}_{\mathcal{W}}(x_{k}^{i})=\mathcal{N}(\bar{x}^{i}_{k},\bar{\Gamma}^{i}_{k})$.
The elliptical probability contour containing $p_{\text{safe}}$ probability
mass in the workspace marginal distribution is the region
$\mathcal{E}_{\text{safe}}^{i}$ such that
$\int_{\mathcal{E}_{\text{safe}}^{i}}\textbf{b}_{\mathcal{W}}(x_{k}^{i})dx=p_{\text{safe}}.$
The axes of this ellipsoid are given by $a_{l}^{i}=t_{\chi}\lambda_{l},\quad
l=1,\ldots,w,$ where $\lambda_{l}$ is an eigenvalue of $\bar{\Gamma}_{k}^{i}$,
and $t_{\chi}$ is the inverse $\chi^{2}$ cumulative density function evaluated
at $p_{\text{safe}}$ with $w$ degrees of freedom.

Determining intersection of ellipses is difficult, so we bound the ellipsoid
with a sphere, where checking for intersection only requires comparing the sum
of the radii to the distance between the centers. The sphere containing this
ellipsoid, $\mathcal{C}_{\text{safe}}^{i}$, is defined by the radius
$r^{i}=t_{\chi}\lambda_{max}$, where
$\lambda_{\max}=\max\\{\lambda_{l}\\}_{l=1}^{w}.$ Because
$\mathcal{E}_{\text{safe}}^{i}\subseteq\mathcal{C}_{\text{safe}}^{i}$, it is
assured that the probability mass contained by this sphere is greater than or
equal to $p_{\text{safe}}$, i.e.,
$\int_{\mathcal{C}_{\text{safe}}^{i}}\textbf{b}_{\mathcal{W}}(x_{k}^{i})dx\geq
p_{\text{safe}}.$ An illustration of both the elliptical contour and its
associated spherical bound are shown in Fig. 3(b) for the robots in Fig 3(a).

Note that this sphere contains only the probability mass associated with
$x_{k}^{i}$, the origin point on the body frame. It does not account for the
collision probability of the other points on the robot body. To address this,
we inflate the sphere by a bounding sphere on the robot’s body (shown in Fig.
3(a)), such that the probability of _any_ point on the robot’s body being
within this inflated sphere is greater than or equal to $p_{\text{safe}}$. The
radius, $R_{i}$, of the bounding sphere on the robot’s body is obtained by
finding the largest possible Euclidean distance between any pair of points
within a robot’s body:
$R_{i}=\frac{1}{2}\max_{x_{1},x_{2}\in\mathcal{B}^{i}}\|x_{1}-x_{2}\|_{2}$.
The inflated bound, $\mathcal{C}_{bound}^{i}$ is a sphere of radius
$r^{i}+R^{i}$. The probability that any point of the robot is outside this
sphere, the complement of $\mathcal{C}_{bound}^{i}$, is guaranteed to be less
than $1-p_{\text{safe}}$. Thus, if any two robot’s bounding spheres do not
intersect, i.e.
$\mathcal{C}_{bound}^{i}\cap\mathcal{C}_{bound}^{j}=\emptyset$, the
probability of collision between them is guaranteed to be less than
$1-p_{\text{safe}}$. The inflated disks and corresponding intersection
checking are shown in Fig. 3(c).

This method has two main sources of conservatism. The first stems from the
bounding sphere on the robot body. This is especially egregious for robots
with high aspect ratios. The second source is associated with checking for
intersections of the probability mass bound. This conservatively assumes all
probability mass not within the sphere is in collision, so any intersection
(no matter how small) is considered a possible constraint violation. Violation
for this method does not imply that the actual probability of collision
exceeds the constraint. Despite these limitations, this method is extremely
fast.

###  IV-B Method 2: Linear Gaussian Transformation

This section introduces a class of methods that reduce conservatism by
simplifying the integration itself, namely by reducing its dimensionality.
This is done by introducing a new RV representing the difference between two
robots’ states: $x_{k}^{ij}=x_{k}^{i}-x_{k}^{j}$. This RV is Gaussian
distributed such that
$x_{k}^{ij}\sim\mathcal{N}(\check{x}_{k}^{i}-\check{x}_{k}^{j},\Gamma_{k}^{i}+\Gamma_{k}^{j})$,
with the workspace marginal $\textbf{b}_{\mathcal{W}}(x_{k}^{ij})$. This new
RV represents a vector connecting the origin points of two robots. Collisions
correspond to specific realizations of $x_{k}^{ij}$, defining the required
integration area.

####  IV-B1 Method 2.1 (M.2.1): Convex Polytopic Bounding

The difference RV somewhat simplifies the high dimensional c-space obstacle
representation, shown in Fig. 2(a), by merging the workspace dimensions for
the two robots. But the orientation of the robots must still be considered. To
obtain a tractable integration area, we use the same bounding spheres
described in Sec. IV-A. Intersection of the bounding spheres is considered a
collision state, corresponding to a spherical integration area
$S_{body}^{ij}$, with radius $R^{i}+R^{j}$, over difference distribution
$\textbf{b}_{\mathcal{W}}(x_{k}^{ij})$ (shown in Fig. 3(d)).

Because exact integration of a Gaussian distribution over a sphere is
difficult, we bound $S_{body}^{ij}$ with a polytope,
$\mathcal{P}_{body}^{ij}$, which can be defined by a set of $N_{s}$ half
planes, $\mathcal{P}_{body}^{ij}=\\{x\mid c_{h,i}^{T}x<d_{h,i},\forall
h\in[1,N_{s}]\\}$, where $c_{h,i}\in\mathbb{R}^{n\times 1}$ and
$d_{h,i}\in\mathbb{R}$. This reduces the robot-robot chance constraint
checking problem to the same problem addressed in [21], where various
(over-)approximations are used to reduce the probabilistic chance constraint
checking to deterministic linear constraints checking. We employ the same
methods to check the chance constraints, with the distribution given by
$\mathcal{N}(\check{x}_{k}^{i}-\check{x}_{k}^{j},\Gamma_{k}^{i}+\Gamma_{k}^{j})$,
and polytope given by $\mathcal{P}_{body}^{ij}$.

A key approximation in this method is the allocation of the acceptable
probability of collision, or risk allocation: the total acceptable probability
of collision $p_{coll}=1-p_{\text{safe}}$ must be distributed across all the
possible sources of collision. In our case, we must distribute the probability
of collision for a robot among the $N_{A}-1$ other robots. We denote the
allowable probability for robot $i$ colliding with any other robot $j$ as
$p_{c}^{ij}$. Two proposed methods of determining $p_{c}^{ij}$ are described
in Sec. IV-C.

This method has two main sources of conservatism. The first occurs with the
spherical bound on the body, the same bound used in M.1 (Sec. IV-A) while the
second arises from the use of deterministic linear constraints for the
integration. This source of conservatism scales with the number of robots and
the number of half planes in the polytope bound. This method is only slightly
more computationally expensive than the first method: while checking linear
constraints is very computationally simple, this method must at worst check
every constraint associated with each pair of robots. This yields slightly
longer computation times than the M.1, which requires checking a single
condition for every robot pair.

####  IV-B2 Method 2.2 (M.2.2): Grid Integration

This method follows directly from M.2.1, but with a direct integration
technique, rather than a check on deterministic linear constraints. This makes
this method more accurate and less conservative, but also more computationally
intensive.

We begin with the same difference distribution over
$\textbf{b}_{\mathcal{W}}(x_{k}^{ij})$ and integration over
$\mathcal{P}_{body}^{ij}$. We then perform a whitening (Mahalanobis)
transformation on the distribution and $\mathcal{P}_{body}^{ij}$. This
transformation decorrelates a Gaussian distribution, yielding the standard
normal distribution, as well as rotating and translating
$\mathcal{P}_{body}^{ij}$. We define this new integration area as
$\mathcal{P}_{white}^{ij}$ (see Fig. 3(f)). Because we have transformed the
probability distribution to a standard normal, the cumulative distribution
function can be easily computed via the error function (erf); we do this over
a grid. This technique is agnostic to the choice of grid cells and only
requires rectangular cells that completely cover the polytope
$\mathcal{P}_{white}^{ij}$, as shown in green in Fig. 3(f). The probability of
collision is the sum of the probability mass contained in each grid cell,
defined as $p_{poly}$. The chance constraint can then be simply checked as
$p_{poly}\leq p_{c}^{ij}$, where $p_{c}^{ij}$ again is the risk allocated to a
pair of robots.

This method inherits the same conservatism associated with the spherical body
bound and polytope integration area bound in M.2.1. However, the approximate
grid calculation is far less conservative than the deterministic linear
constraints. The grid introduces more computational cost, particularly in
transforming the distribution and evaluating the probability of each grid
cell. Finer grid discretization reduces conservatism, but also increases
computation time.

###  IV-C Risk Allocation

Risk allocation is only necessary for M.2.1 and M.2.2 and not required in M.1.
Equally distributing the allowable probability of collision $p_{coll}$ across
the robots and obstacles, as introduced in [21], is the simplest allocation
method. For robot-robot collision this corresponds to
$p_{c}^{ij}=\frac{p_{coll}}{N_{O}+N_{A}-1}$. However, this can be overly
conservative and make it difficult to find viable motion plans. Instead, we
propose allocating $p_{coll}$ more effectively, such that robots with a higher
likelihood of collision are assigned a higher proportion of $p_{coll}$.

We begin by setting $p_{coll}=P_{c}^{A}+P_{c}^{O}$, where $P_{c}^{A}$ and
$P_{c}^{O}$ are the allowable probabilities of collision with all the robots
and all the static obstacles, respectively. This naturally fits the high-level
low-level division in CC-K-CBS because obstacle collision checking occurs only
in the low-level planner. Define the total volume of all the robots as
$V_{A}$, and the volume of all the obstacles as $V_{O}$. Then, we set
$P_{c}^{A}=\frac{V_{A}p_{coll}}{V_{A}+V_{O}}$ and
$P_{c}^{O}=\frac{V_{O}p_{coll}}{V_{A}+V_{O}}$. This allocates more of the
probability of collision to the category that is more likely to cause
collisions. We then divide $P_{c}^{A}$ based on the distance between robots,
with robots that are closer together receiving a larger portion of
$P_{c}^{A}$. Let $d_{k}^{j}=\|\bar{x}_{k}^{i}-\bar{x}_{k}^{j}\|$ be the
workspace distance between robot $i$’s and $j$’s means. Then, the allowable
probability of collision for robot $i$ with $j$ is
$p_{c}^{ij}=\frac{\alpha}{d_{k}^{j}}P_{c}^{A}$, where
$\alpha=1/(\sum_{j}^{N_{A}-1}d_{k}^{j})$ is the normalizing factor.

##  V EVALUATIONS

We evaluate our algorithms in several benchmarks. We first independently test
the robot-robot collision checking algorithms in Sec. IV, allowing us to
isolate the conservatism and computational complexity of the robot-robot
collision checkers. We then analyze the full algorithms from Sec. III with
different collision checking methods. Our implementation of CC-K-CBS is
publicly available [22]. It uses belief-SST as the low-level planner in Open
Motion Planning Library (OMPL) [23]. The benchmarks were performed on AMD 4.5
GHz CPU and 64 GB of RAM.

(a)

(b)

Figure 4: (a) Plot of the conservatism $C_{z}$ for each method in each
environment. (b) Box and whiskers plots of the log(computation time) for each
method.

We consider square-shaped robots of width 0.25 with two types of dynamics: 2D
robotic system taken from [14];
$x_{k+1}^{i}=x_{k}^{i}+u_{k}^{i}+w_{k}^{i},w_{k}^{i}\sim\mathcal{N}(0,0.1^{2}I)$,
and 2nd-order unicycle;
$\dot{\text{x}}^{i}=\text{v}\cos\theta,\dot{\text{y}}^{i}=\text{v}\sin\theta,\dot{\theta}=\omega,\dot{\text{v}}=a$,
where $\omega$ and $a$ are control inputs. We use the feedback linearization
scheme described in [24] to obtain a linear model, which we convert to
discrete time with additive noise. We assume both systems are fully observable
with measurement model in (2) with $C^{i}=I$ and $R^{i}=0.1^{2}I$.

Robot-Robot Collision Checker Benchmark:  We characterize the robot-robot
collision checking algorithms based on computation time and conservatism.
Conservatism impacts the planners ability to find paths in more difficult,
i.e. cluttered, situations. For a given set of beliefs for two robots, we
define the conservatism of method $z$ as $C_{z}=R_{z}-R_{MC}$, where $R_{z}$
is the rejection rate obtained by method $z$, and $R_{MC}$ is the rejection
rate obtained by a Monte Carlo evaluation of the probability of collision.

We sampled 2D beliefs in two difference spaces, 5$\times$5 and 10$\times$10,
with $p_{\text{safe}}=0.95$. For M.2.2 (the gridded cdf method), we chose 6
different grid discretizations, denoted by M.2.2($d$). For each
discretization, the maximum range of $\mathcal{P}_{whit}^{ij}$ over each axis,
$D_{w}$, was divided equally, such that the discretization step size was
$D_{w}/d$. Note that each of the proposed methods decreases in conservatism,
shown in Fig. 4(a), and increases in computation time, shown in Fig. 4(b). The
5$\times$5 space is the more ‘cluttered’, with more beliefs in close
proximity, which in turn results in higher collision probabilities. Note that
the 5$\times$5 space generally results in lower conservatism. This is because
for the relatively high choice of $p_{\text{safe}}$, a high proportion of the
sampled beliefs truly violated the chance constraint, making it less
conservative compared to a scenario where more beliefs are valid.

(a) Env8

(b) Env32Obs

Figure 5: (a) Sample plan for 4 robots with 2nd-order unicycle dynamics in
Env8. Ellipses are $95\%$ probability contours. (b) Sample plan for 20 robots
with 2nd-order unicycle dynamics in Env32Obs. Ellipses are $95\%$ probability
contours. TABLE I: CC-K-CBS benchmarks for Env8.

2D Linear
---
$N_{A}$ | Metric | M.1 | M.2.1 | M.2.1* | M.2.2(2) | M.2.2(5) | M.2.2(10)
2 | Succ. Rate | 1.00 | 1.00 | 1.00 | 1.00 | 1.00 | 1.00
Time (s) | 1.08 | 0.68 | 1.04 | 1.11 | 0.69 | 0.74
4 | Succ. Rate | 0.00 | 1.00 | 1.00 | 1.00 | 1.00 | 1.00
Time (s) | - | 18.68 | 35.03 | 4.19 | 3.97 | 7.61
6 | Succ. Rate | - | 0.00 | 0.12 | 1.00 | 1.00 | 1.00
Time (s) | - | - | 66.68 | 4.44 | 6.15 | 12.72
8 | Succ. Rate | - | - | 0.00 | 1.00 | 1.00 | 0.86
Time (s) | - | - | - | 15.38 | 26.54 | 59.08
Unicycle
2 | Succ. Rate | 1.00 | 1.00 | 1.00 | 1.00 | 1.00 | 1.00
Time (s) | 0.95 | 0.63 | 0.61 | 0.37 | 0.57 | 0.83
4 | Succ. Rate | 0.00 | 0.96 | 1.00 | 1.00 | 1.00 | 1.00
Time (s) | - | 17.67 | 9.81 | 2.59 | 2.86 | 5.24
6 | Succ. Rate | - | 0.00 | 0.00 | 1.00 | 1.00 | 1.00
Time (s) | - | - | - | 4.79 | 5.95 | 10.40
8 | Succ. Rate | - | - | - | 1.00 | 0.98 | 0.96
Time (s) | - | - | - | 18.27 | 25.52 | 45.73

Planners Benchmark:  We use benchmarking to evaluate the proposed planners
under each of the collision-checking methods in Sec. IV. We additionally
implement the adaptive risk allocation on M.2.1, denoted by M.2.1*. All other
methods use equal allocation. Benchmarking was performed using 50 runs for
each algorithm with a maximum planning time of 3 minutes and
$P_{\text{safe}}=0.9$. We examine 3 environments: a small 8$\times$8 (Env8),
large 32$\times$32 with 50 random obstacles (Env32Obs), and large 32$\times$32
without obstacles (Env32) shown in Figs. 5(a), 5(b), and 1, respectively. The
considered metrics are _success rate_ (the ratio of the successful planning
instances within 3 minutes to all instances) and _time_ (computation time of
the successful runs).

Env8:  This small environment allows us to study how the different collision
checking methods perform in a cluttered environment with the 2D system. The
results are shown in Table I. Conservatism strongly impacts the planner’s
ability to find paths with additional robots. All methods reliably find plans
for 2 robots, however as the conservatism of the method increases, the success
rate decreases when more agents are added. In this cluttered environment robot
paths are necessarily close together, so more conservative methods are less
likely to find valid paths. Specifically, M.1 performs poorly and M.2.2 has
the best performance.

Env32Obs:  The large environment is far less cluttered, allowing for far more
agents and better insight into the scalability of the CC-K-CBS algorithm. The
results for Env8 are based on the robot-obstacle collision checker from [21].
However, this is very conservative, especially for many obstacles. CC-K-CBS
with this method was unable to find paths in Env32Obs for even 2 robots. To
scale to a larger number of robots, we made the very straightforward
adaptation to M.1 to check intersection of the bound $\mathcal{C}_{bound}^{i}$
with obstacles. As shown in Fig 6, CC-K-CBS is able to scale to $20$ robots
with $68\%$ success rate and runtime of $120.7\pm 5.45$ seconds for the 2nd-
order unicycle dynamics. A sample plan is shown in Fig. 5(b). This illustrates
that for large number of robots, M.1 has the best performance because it is
not affected by risk allocation, whereas the allocated risk gets prohibitively
smaller as the number of robots increases for the other methods, making it
difficult to find valid plans.

We also performed benchmarking on the centralized approach with the simpler 2D
system in Env32Obs, using collision checking M.1. This algorithm was able to
plan for 2 robots with $100\%$ success rate, but the success rate dropped to
$8\%$ for 3 robots, illustrating the classical scalability issue with
centralized approaches.

Env32:  To further test the scalability of CC-K-CBS, we increased the number
of agents in the large open environment Env32. As shown in Fig 6, CC-K-CBS is
able to to scale to $28$ 2nd-order unicycle robots with $56\%$ success rate
and $129\pm 4.97$ seconds computation time. A sample plan is shown in Fig. 1.
CC-K-CBS is also able to plan for 30 robots but at much lower success rate of
$18\%$.

Figure 6: Benchmark results for CC-K-CBS on varying number of 2nd-oder
unicycle robots. For # robots $\leq 20$, Env32Obs is used (dashed line), and
for # robots $\geq 22$, Env32 is used (solid line). The discontinuity is due
to removing obstacles.

We validated the robustness of the motion plans using Monte Carlo simulation.
We randomly selected motion plans, and with 500 simulations for each robot and
observed no more than $2\%$ collision probability over the entire trajectory,
illustrating the plans are indeed robust to uncertainty.

##  VI CONCLUSION

In this work, we consider the MRMP problem under Gaussian uncertainty and
introduce CC-K-CBS, which generates motion plans for each robot with chance-
constraint guarantees. This algorithm is based on extending K-CBS to
accommodate belief space planning. Our main contribution includes three
methods for fast evaluation of collision probabilities between robots. This
offline algorithm successfully generated motion plans across a variety of
scenarios and scaled to 30 robots. Future directions could encompass other
types of noise and nonlinear dynamics as well as new collision checking
methods that could further improve computation time and conservatism.

## References

  * [1] S. M. LaValle, “Rapidly-exploring random trees: A new tool for path planning,” Computer Science Department, Iowa State University, Tech. Rep., 1998\.
  * [2] G. Wagner and H. Choset, “Subdimensional expansion for multirobot path planning,” _Artificial Intelligence_ , vol. 219, pp. 1–24, 2015.
  * [3] R. Shome, K. Solovey, A. Dobson, D. Halperin, and K. E. Bekris, “drrt*: Scalable and informed asymptotically-optimal multi-robot motion planning,” _Autonomous Robots_ , vol. 44, no. 3, pp. 443–467, 2020.
  * [4] J. Kottinger, S. Almagor, and M. Lahijanian, “MAPS-X: Explainable multi-robot motion planning via segmentation,” in _Int’t Conf. on Robotics and Automation_. IEEE, May 2021\.
  * [5] I. A. Sucan and L. E. Kavraki, “A sampling-based tree planner for systems with complex dynamics,” _IEEE Transactions on Robotics_ , vol. 28, no. 1, pp. 116–131, 2011.
  * [6] J. D. Gammell, S. S. Srinivasa, and T. D. Barfoot, “Batch informed trees (bit*): Sampling-based optimal planning via the heuristically guided search of implicit random geometric graphs,” in _Int. Conf. on Robotics and Automation_ , 2015, pp. 3067–3074.
  * [7] S. Tang and V. Kumar, “A complete algorithm for generating safe trajectories for multi-robot teams,” in _Robotics Research_. Springer, 2018, pp. 599–616.
  * [8] D. Le and E. Plaku, “Cooperative multi-robot sampling-based motion planning with dynamics,” _ICAPS_ , vol. 27, no. 1, pp. 513–521, 2017.
  * [9] ——, “Multi-robot motion planning with dynamics via coordinated sampling-based expansion guided by multi-agent search,” _IEEE Robot. and Auto. Letters_ , vol. 4, no. 2, pp. 1868–1875, 2019.
  * [10] J. van den Berg and M. Overmars, “Prioritized motion planning for multiple robots,” in _2005 IEEE/RSJ International Conference on Intelligent Robots and Systems_ , 2005, pp. 430–435.
  * [11] S. M. LaValle, _Planning algorithms_. Cambridge university press, 2006.
  * [12] J. Kottinger, S. Almagor, and M. Lahijanian, “Conflict-based search for multi-robot motion planning with kinodynamic constraints,” in _Int. Conf. on Intelligent Robots and Systems_ , 2022, pp. 13 494–13 499.
  * [13] A.-A. Agha-Mohammadi, S. Chakravorty, and N. M. Amato, “Firm: Sampling-based feedback motion-planning under motion uncertainty and imperfect measurements,” _Int. J. Rob. Res._ , vol. 33, no. 2, p. 268–304, feb 2014\.
  * [14] A. Bry and N. Roy, “Rapidly-exploring random belief trees for motion planning under uncertainty,” in _2011 IEEE International Conference on Robotics and Automation_ , 2011, pp. 723–730.
  * [15] Q. H. Ho, Z. N. Sunberg, and M. Lahijanian, “Gaussian belief trees for chance constrained asymptotically optimal motion planning,” in _Int. Conf. on Robotics and Automation_ , 2022, pp. 11 029–11 035.
  * [16] C. Boutilier, “Planning, learning and coordination in multiagent decision processes,” in _Conf. Theo. Aspects of Rationality and Knowledge_ , San Francisco, CA, USA, 1996, p. 195–210.
  * [17] A. Patwardhan, R. Murai, and A. J. Davison, “Distributing collaborative multi-robot planning with gaussian belief propagation,” _IEEE Robotics and Automation Letters_ , vol. 8, no. 2, pp. 552–559, 2023.
  * [18] J. Alonso-Mora, P. Beardsley, and R. Siegwart, “Cooperative collision avoidance for nonholonomic robots,” _IEEE Transactions on Robotics_ , vol. 34, no. 2, pp. 404–420, 2018.
  * [19] B. Şenbaşlar, W. Hönig, and N. Ayanian, “Robust trajectory execution for multi-robot teams using distributed real-time replanning,” in _Distributed Autonomous Robotic Systems_ , N. Correll, M. Schwager, and M. Otte, Eds. Cham: Springer, 2019, pp. 167–181.
  * [20] R. Van Parys and G. Pipeleers, “Online distributed motion planning for multi-vehicle systems,” in _2016 European Control Conference (ECC)_ , 2016, pp. 1580–1585.
  * [21] L. Blackmore, M. Ono, and B. C. Williams, _Chance-constrained optimal path planning with obstacles_. IEEE, 2011, pp. 1080–1094.
  * [22] J. Kottinger, “Kinodynamic Conflict-Based Search,” 2022. [Online]. Available: https://github.com/aria-systems-group/K-CBS
  * [23] I. A. Şucan, M. Moll, and L. E. Kavraki, “The Open Motion Planning Library,” _IEEE Robotics & Automation Magazine_, vol. 19, no. 4, pp. 72–82, December 2012, https://ompl.kavrakilab.org.
  * [24] A. De Luca, G. Oriolo, and M. Vendittelli, “Stabilization of the unicycle via dynamic feedback linearization,” _IFAC Proceedings Volumes_ , vol. 33, no. 27, pp. 687–692, 2000.


# Conditional Proof of the Ergodic Conjecture for Falling Ball Systems

Nandor Simanyi  1402 10th Avenue South
Birmingham AL 35294-1241<EMAIL_ADDRESS>
###### Abstract.

In this paper we present a conditional proof of Wojtkowski’s Ergodicity
Conjecture for the system of $1D$ perfectly elastic balls falling down in a
half line under constant gravitational acceleration, [W85], [W86], [W90a],
[W90b], [W98]. Namely, we prove that almost every such system is (completely
hyperbolic and) ergodic, by assuming the transversality between different
singularities and between singularities and stable (unstable) invariant
manifolds.

###### 2010 Mathematics Subject Classification:

37D05

##  1\. Introduction/Prerequisites

In his paper [W90a] M. Wojtkowski introduced the following Hamiltonian
dynamical system with discontinuities: There is a vertical half line
$\left\\{q|\,q\geq 0\right\\}$ given and $n$ ($\geq 2$) point particles with
masses $m_{1}\geq m_{2}\geq\dots\geq m_{n}>0$ and positions $0\leq q_{1}\leq
q_{2}\leq\dots\leq q_{n}$ are moving on this half line so that they are
subjected to a constant gravitational acceleration $a=-1$ (they fall down),
they collide elastically with each other, and the first (lowest) particle also
collides elastically with the hard floor $q=0$. We fix the total energy

$H=\sum_{i=1}^{n}\left(m_{i}q_{i}+\frac{1}{2m_{i}}p_{i}^{2}\right)$

by taking $H=1$. The arising Hamiltonian flow with collisions
$\left(\mathbf{M},\\{\psi^{t}\\},\mu\right)$ ($\mu$ is the Liouville measure)
is the studied model of this paper.

Before formulating the result of this article, however, it is worth mentioning
here three important facts:

  1. $(1)$

Since the phase space $\mathbf{M}$ is compact, the Liouville measure $\mu$ is
finite.

  2. $(2)$

The phase points $x\in\mathbf{M}$ for which the trajectory
$\\{\psi^{t}(x)|,t\in\mathbb{R}\\}$ hits at least one singularity (i. e. a
multiple collision) are contained in a countable union of proper, smooth
submanifolds of $\mathbf{M}$ and, therefore, such points form a set of $\mu$
measure zero.

  3. $(3)$

For $\mu$-almost every phase point $x\in\mathbf{M}$ the collision times of the
trajectory $\\{\psi^{t}(x)|,t\in\mathbb{R}\\}$ do not have any finite
accumulation point, see Proposition A.1 of [S96].

In the paper [W90a] Wojtkowski formulated his main conjecture pertaining to
the dynamical system $\left(\mathbf{M},\\{\psi^{t}\\},\mu\right)$:

######  Conjecture 1.1 (Wojtkowski’s Conjecture).

If $m_{1}\geq m_{2}\geq\dots\geq m_{n}>0$ and $m_{1}\neq m_{n}$, then all but
one characteristic (Lyapunov) exponents of the flow
$\left(\mathbf{M},\\{\psi^{t}\\},\mu\right)$ are nonzero. Futhermore, the
system is ergodic.

######  Remark 1.1.

1\. The only exceptional exponent zero must correspond to the flow direction.

2\. The condition of nonincreasing masses (as above) is essential for
establishing the invariance of the symplectic cone field — an important
condition for obtaining nonzero characteristic exponents. As Wojtkowski
pointed out in Proposition 4 of [W90a], if $n=2$ and $m_{1}<m_{2}$, then there
exists a linearly stable periodic orbit, thus dimming the chances of proving
ergodicity.

Our paper will very closely follow all the notations and definitions of the
atricle [W90a], so the reader is warmly recommended to be familiar with that
paper.

##  2\. Eventually Strict Cone Invariance

We will be working with the symplectic coordinates $(\delta h,\delta v)$ for
the tangent vectors of the reduced phase space $\mathbf{M}$ satisfying the
usual reduction equations $\sum_{i=1}^{n}\delta h_{i}=0=\sum_{i=1}^{n}\delta
v_{i}$.

######  Remark 2.1.

The coordinates $\delta h_{i}$ and $\delta v_{i}$ ($i=1,2,\dots,n$) serve as
suitable symplectic coordinates in the codimension-one subspace
$\mathcal{T}_{x}$ of the full tangent space $\mathcal{T}_{x}\mathbf{M}$ of
$\mathbf{M}$ at $x$. Recall that the $(2n-2)$-dimensional vector space
$\mathcal{T}_{x}$ is transversal to the flow direction, and the restriction of
the canonical symplectic form

$\omega=\sum_{i=1}^{n}\delta q_{i}\wedge\delta p_{i}=\sum_{i=1}^{n}\delta
h_{i}\wedge\delta v_{i}$

of $\mathbf{M}$ is non-degenerate on $\mathcal{T}_{x}$, see [W90a]. We also
recall that

(2.1) $\delta h_{i}=m_{i}\delta q_{i}+m_{i}v_{i}\delta v_{i}=m_{i}\delta
q_{i}+v_{i}\delta p_{i}.$

Corresponding to the above choice of symplectic coordinates, the considered
monotone Q-form will be

(2.2) $Q_{1}(\delta h,\delta v)=\sum_{i=1}^{n}\delta h_{i}\delta v_{i}.$

It is clear that the evolution of $DS^{t}(\delta h(0),\delta v(0))=(\delta
h(t),\delta v(t))$ between colisions is

(2.3) $\frac{d}{dt}\left(\delta h(t),\delta v(t)\right)=(0,0).$

If a collision of type $(i,i+1)$ ($i=1,2,\dots,n-1$) takes place at time
$t_{k}$, then the derivative of the flow at the collision $\delta
h^{-}(t_{k})\mapsto\delta h^{+}(t_{k})$, $\delta v^{-}(t_{k})\mapsto\delta
v^{+}(t_{k})$ is given by the matrices

(2.4) $\displaystyle\delta h^{+}(t_{k})=R^{*}_{i}\left[\delta
h^{-}(t_{k})+S_{i}\delta v^{-}(t_{k})\right]$ $\displaystyle\delta
v^{+}(t_{k})=R_{i}\delta v^{-}(t_{k}),$

where the matrix $R_{i}$ is the $n\times n$ identity matrix, except that its
$2\times 2$ submatrix at the crossings of the $i$-th and $(i+1)$-st rows and
columns is

$R_{i}^{(i,i+1)}=\left[\begin{array}[]{cc}\gamma_{i}&1-\gamma_{i}\\\
1+\gamma_{i}&-\gamma_{i}\end{array}\right]$

with $\gamma_{i}=\dfrac{m_{i}-m_{i+1}}{m_{i}+m_{i+1}}$. The matrix $S_{i}$ is,
similarly, the $n\times n$ zero matrix, except its $2\times 2$ submatrix at
the crossings of the $i$-th and $(i+1)$-st rows and columns, which takes the
form of

$S_{i}^{(i,i+1)}=\left[\begin{array}[]{cc}\alpha_{i}&-\alpha_{i}\\\
-\alpha_{i}&\alpha_{i}\end{array}\right]$

with

(2.5)
$\alpha_{i}=\frac{2m_{i}m_{i+1}(m_{i}-m_{i+1})}{(m_{i}+m_{i+1})^{2}}(v_{i}^{-}-v_{i+1}^{-})>0.$

These formulas can be found, for example, in Sention 4 of [W90a].

Concerning a floor collision $(0,1)$ at time $t_{k}$, the transformations are

(2.6) $\displaystyle\delta h_{1}^{+}(t_{k})=\delta h_{1}^{-}(t_{k})$
$\displaystyle\delta v_{1}^{+}(t_{k})=\delta v_{1}^{-}(t_{k})+\frac{2\delta
h_{1}^{-}(t_{k})}{m_{1}v_{1}^{+}(t_{k})},$

see, for instance, Section 4 of [W90a] or [W98].

In this section we will be studying _non-singular_ trajectory segments

$S^{[0,T]}x_{0}=\left\\{x_{t}=S^{t}x_{0}\big{|}\;0\leq T\leq T\right\\}$

of the flow $\\{S^{t}\\}$ with the symbolic collision sequence
$\Sigma=(\sigma_{1},\sigma_{2},\dots,\sigma_{N})$, where
$\sigma_{k}=(i_{k},i_{k+1})$, $0\leq i_{k}\leq n-1$, $k=1,2,\dots,N$. The
collision graph $\mathcal{G}=\mathcal{G}(\Sigma)$ of $\Sigma$ has the set
$\\{0,1,\dots,n\\}$ as its vertex set, and the unoriented edges of
$\mathcal{G}$ are the unordered pairs $\\{i_{k},i_{k+1}\\}$, $k=1,2,\dots,N$,
counted without multiplicity.

According to Wojtkowski’s arguments between the Theorem and Proposition I of
Section 5 of [W90a], in order to prove the strict invariance of the cone field
$C_{1}=\\{Q_{1}\geq 0\\}$ along the considered trajectory segment

$S^{[0,T]}x_{0}=\left\\{x_{t}=S^{t}x_{0}\big{|}\;0\leq T\leq T\right\\},$

it is enough to prove that

  1. $(A)$

for every non-zero vector $(0,\delta v)\in\mathcal{T}_{x_{0}}$ it is true that
$Q_{1}(DS^{T}(0,\delta v))>0$, and

  2. $(B)$

for every non-zero vector $(\delta h,0)\in\mathcal{T}_{x_{0}}$ it is true that
$Q_{1}(DS^{T}(\delta h,0))>0$.

Moreover, Wojtkowski’s mentioned arguments in Section 5 of [W90a] actually
prove property (A) above in the case when the collision graph
$\mathcal{G}(\Sigma)$ restricted to the vertex set $\\{1.2.\dots.n\\}$ is
connected, i. e. all collisions $(i,i+1)$ with $i>0$ occur. Here we briefly
rephrase his ideas:

Formula 2.6 shows that a tangent vector of the form $(0,\delta v)$ is
unchanged at any floor collision. Suppose now that a collision $(i,i+1)$
($1\leq i\leq n$) occurs at time $t_{k}$ and $\delta h^{-}(t_{k})=0$. Equation
2.4 shows that, after pushing the tangent vector $(0,\delta v^{-}(t_{k}))$
through the collision $(i,i+1)$, either

(2.7) $Q_{1}(\delta h^{+}(t_{k}),\delta v^{+}(t_{k}))=\alpha_{i}(\delta
v_{i}^{-}(t_{k})-\delta v_{i+1}^{-}(t_{k}))^{2}>0,$

or

(2.8) $\delta v_{i}^{-}(t_{k})=\delta v_{i+1}^{-}(t_{k})=\delta
v_{i}^{+}(t_{k})=\delta v_{i+1}^{+}(t_{k}).$

Thus, as long as $Q_{1}(DS^{T}(0,\delta v(0)))=0$ and the restriction of the
collision graph $\mathcal{G}$ on the vertex set $\\{1,2,\dots,n\\}$ is
connected, we have that

(2.9) $\delta v_{i}(t)=\delta v_{1}(0)$

for all $t$, $0\leq t\leq T$ and all $i=1,2,\dots,n$. Therefore, $\delta
v_{i}(t)=0$ for all $t$, $0\leq t\leq T$ and all $i=1,2,\dots,n$, thanks to
the usual reduction equation $\sum_{i=1}^{n}\delta v_{i}=0$. This finishes the
proof of Property (A) in the case when the collision graph
$\mathcal{G}(\Sigma)$ contains all edges $(i,i+1)$ with $1\leq i\leq n-1$.

######  Remark 2.2.

From now on, we will always assume that the collision graph $\mathcal{G}$
contains all collisions $(i,i+1)$ with $i>0$.

Checking Property (B) for the non-singular trajectory segment $S^{[0,T]}x_{0}$
is, however, a bit harder and it requires a bit more on the side of the
collision graph $\mathcal{G}(\Sigma)$ than simply containing all edges
$(i,i+1)$, $i=1,2,\dots,n-1$.

First of all, we define the linear space of all _neutral vectors_ as follows:

(2.10) $\mathcal{N}_{x_{0}}^{T}=\left\\{(\delta h(0),\delta
v(0))\in\mathcal{T}_{x_{0}}\big{|}\;Q_{1}(\delta h(t),\delta v(t))=0,\text{
for }0\leq t\leq T\right\\}.$

Recall that $(\delta h(t),\delta v(t))=DS^{t}\left((\delta h(0),\delta
v(0))\right)$.

It floows from the time-evolution equations 2.3, 2.4, and 2.6 that

(2.11) $DS^{t}\left((\delta h(0),\delta v(0))\right)=(\delta h(t),0)$

for $0\leq t\leq T$, i.e. $\delta v(t)=0$. Indeed, at any floor collision the
form $Q_{1}$ increases by the amount of

(2.12) $\frac{2(\delta h_{1}^{-}(t_{k}))^{2}}{m_{1}v_{1}^{+}(t_{k})},$

which forces $\delta h_{1}^{-}(t_{k})=\delta h_{1}^{+}(t_{k})=0$ and, as a
corollary, $\delta v_{1}^{-}(t_{k})=\delta v_{1}^{+}(t_{k})$. According to
2.4, at an $(i,i+1)$ collision ($i\geq 1$), occuring at time $t_{k}$, the form
$Q_{1}$ increases by the amount

(2.13) $\alpha_{i}\left(\delta v_{i}^{-}(t_{k})-\delta
v_{i+1}^{-}(t_{k})\right)^{2},$

and this forces $\delta v_{i}^{-}(t_{k})=\delta v_{i+1}^{-}(t_{k})=\delta
v_{i}^{+}(t_{k})=\delta v_{i+1}^{+}(t_{k})$, thus $\delta v(t)=0$ for $0\leq
t\leq T$, as claimed.

So, the neutral vector takes the form $DS^{t}(\delta h(0),0)=(\delta h(t),0)$
for $0\leq t\leq T$, and $\delta h(t)$ only changes at ball-to-ball collisions
$(i,i+1)$ (at time $t_{k}$) according to the law

(2.14) $\delta h^{+}(t_{k})=R_{i}^{*}\delta h^{-}(t_{k}),$

whereas at a floor collision, occuring at time $t_{k}$, the law

(2.15) $\delta h^{-}(t_{k})=\delta h^{+}(t_{k})=0$

is forced by the neutrality requirement $Q_{1}(\delta h(t),\delta v(t))=0$.
Accordingly, we define $R_{0}=I=R_{0}^{*}$, the $n\times n$ identity matrix.

Out of the considered collisions $\sigma_{1},\sigma_{2},\dots,\sigma_{N}$, let
the floor collisions be

$\sigma_{k_{1}},\dots,\sigma_{k_{m}},$

$1\leq k_{1}<\dots<k_{m}\leq N$. The above argument shows that the neutral
space $\mathcal{N}_{x_{0}}^{T}$ can be defined by the following system of
homogeneous linear equations:

(2.16) $\mathcal{N}_{x_{0}}^{T}=\mathcal{N}(\Sigma,\vec{m})\\\
=\left\\{(\delta
h(0),0)\in\mathcal{T}_{x_{0}}\big{|}\;\Pi_{1}R^{*}_{i_{k_{l}}}\dots
R^{*}_{i_{2}}R^{*}_{i_{1}}\delta h(0)=0,\\\ \text{ for }l=1,\dots,m\right\\},$

where $\Pi_{1}(\delta h_{1},\delta h_{2},\dots,\delta h_{n})=\delta h_{1}$ is
the projection onto the first component.

We observe that the neutral space
$\mathcal{N}_{x_{0}}^{T}=\mathcal{N}(\Sigma,\vec{m})$ only depends on the
n-tuple os masses $m_{1}>m_{2}>\dots>m_{n}>0$ and on the symbolic collision
sequence $\Sigma=(\sigma_{1},\sigma_{2},\dots,\sigma_{N})$.

######  Definition 2.3.

The non-singular trajectory segment $S^{[0,T]}x_{0}$ or, equivalently, the
corresponding pair $(\Sigma,\vec{m})$ is said to be _sufficient_ iff

$\mathcal{N}_{x_{0}}^{T}=\mathcal{N}(\Sigma,\vec{m})=\\{0\\}.$

Otherwise, these objects are said to be insufficient.

######  Remark 2.4.

A concequence of the previous definition is that the sufficiency of
$(\Sigma,\vec{m})$ implies that $\Sigma$ must contain all $n$ types of
collisions. Indeed, a missing $(0,1)$ collision would mean that the system of
homogeneous linear equations 2.16 contains no equations, whereas, a missing
collision $(i_{0},i_{0}+1)$ ($i_{0}=1,2,\dots,n-1$) would establish no
connection between the subsystems of particles $\\{1,2,\dots,i_{0}\\}$ and
$\\{i_{0}+1,i_{0}+2,\dots,n\\}$, thus preventing sufficiency.

Summarizing what we have seen so far, we have

######  Proposition 2.17.

The sufficiency of the trajectory segment $S^{[0,T]}x_{0}$ (or, equivalently,
the sufficiency of $(\Sigma,\vec{m})$) is equivalent to the _strict
invariance_ of the cone field

$C_{1}(x_{t})=\left\\{(\delta h,\delta
v)\in\mathcal{T}_{x_{t}}\big{|}\;\sum_{i=1}^{n}\delta h_{i}\delta v_{i}\geq
0\right\\}$

along

$S^{[0,T]}x_{0}=\left\\{x_{t}=S^{t}x_{0}\big{|}\;0\leq t\leq T\right\\},$

that is, this sufficiency exactly means that

$DS^{T}(C_{1}(x_{0}))\subset\text{int}\left(C_{1}(x_{T})\right).$

According to the time evolution equations 2.3, 2.4, 2.6, in order for the
derivative $DS^{T}$ of the flow to be defined on the tangent space
$\mathcal{T}_{x_{0}}$, it is not enough to know the pair $(\Sigma,\vec{m})$,
but one needs to know the relative velocities

(2.18)
$\rho_{k}=\rho_{k}(\sigma_{k})=v^{-}_{i_{k}}(t_{k})-v^{-}_{i_{k}+1}(t_{k})>0$

for all ball-to-ball collisions $\sigma_{k}=(i_{k},i_{k}+1)$ (they play a role
in 2.4 as a part of $\alpha_{i_{k}}$) and the velocities

(2.19) $\rho_{k}=\rho_{k}(\sigma_{k})=v_{1}^{+}(t_{k})$

for any floor collision $\sigma_{k}=(0,1)$ that play a role in 2.6,
$k=1,2,\dots,N$.

Therefore, the _strict cone invariance_ formulated in 2.17 is a property of
the _extended symbolic sequence_

(2.20)
$\tilde{\Sigma}=\left((\sigma_{1},\rho_{1}),(\sigma_{2},\rho_{2}),\dots(\sigma_{N},\rho_{N})\right)$

and the vector of the masses $\vec{m}$.

The characterization of sufficiency with the strict cone invariance 2.17 has
the big advantage that it appears to be _combinatorially monotone_ , i.e. if
one inserts an additional collision $(\sigma^{*},\rho^{*})$ into the extended
symbolic sequence

$\tilde{\Sigma}=\left((\sigma_{1},\rho_{1}),(\sigma_{2},\rho_{2}),\dots,(\sigma_{N},\rho_{N})\right)$

between $(\sigma_{k},\rho_{k})$ and $(\sigma_{k+1},\rho_{k+1})$, then the
sufficiency will not be lost by this insertion:

######  Proposition 2.21 (Combinatorial Monotonicity Property, CMP).

Assume that, for a given mass distribution $m_{1}>m_{2}>\dots>m_{n}$, the
extended symbolic sequence

$\tilde{\Sigma}=\left((\sigma_{1},\rho_{1}),(\sigma_{2},\rho_{2}),\dots,(\sigma_{N},\rho_{N})\right)$

is sufficient. Let the extended symbolic sequence

$\tilde{\Sigma}^{*}=\left((\sigma_{1},\rho_{1}),\dots,(\sigma_{k},\rho_{k}),(\sigma^{*},\rho^{*}),(\sigma_{k+1},\rho_{k+1}),\dots,(\sigma_{N},\rho_{N})\right)$

be obtained from $\tilde{\Sigma}$ by the insertion of $(\sigma^{*},\rho^{*})$,
as indicated above. Then $\tilde{\Sigma}^{*}$ also satifies the strict cone
invariance property formulated in 2.17.

###### Proof.

For any number $l$, $l=1,2,\dots,N$, let $D_{l}$ be the derivative of the flow
(in terms of $(\delta h,\delta v)$, as always) resulting from the collision
$(\sigma_{l},\rho_{l})$, and $D^{*}$ be the derivative of the flow resulting
from the collision $(\sigma^{*},\rho^{*})$. We have

$\displaystyle D_{N}D_{N-1}\dots
D_{1}(C_{1}(x_{0}))\subset\rm{int}(C_{1}(x_{T})),$ $\displaystyle
D^{*}D_{k}D_{k-1}\dots D_{1}(C_{1}(x_{0}))\subset D_{k}D_{k-1}\dots
D_{1}(C_{1}(x_{0})),$

thus

$D_{N}\dots D_{k+1}D^{*}D_{k}\dots D_{1}(C_{1}(x_{0}))\subset
D_{N}D_{N-1}\dots D_{1}(C_{1}(x_{0}))\subset\rm{int}(C_{1}(x_{T})).$

∎

Since the coefficients of the system of homogeneous linear equations 2.16 are
given rational functions of the masses $m_{1},m_{2},\dots,m_{n}$, we
immediately obtain

######  Proposition 2.22 (Dichotomy).

For any given symbolic sequence

$\Sigma=(\sigma_{1},\dots,\sigma_{N})$

either

  1. $(D_{1})$

for almost every n-tuple of masses $\vec{m}$ ($m_{1}>m_{2}>\dots>m_{n}>0$, the
exceptional set being a proper algebraic subset of $\mathbb{R}^{n}$) it is
true that $\mathcal{N}(\Sigma,\vec{m})=\\{0\\}$, or

  2. $(D_{2})$

for every mass vector $\vec{m}$ we have
$\mathcal{N}(\Sigma,\vec{m})\neq\\{0\\}$.

######  Definition 2.5.

In the case $(D_{1})$ above the symbolic sequence $\Sigma$ itself is said to
be sufficient, otherwise it is said to be insufficient.

Thanks to the characterization result 2.17 and the Combinatorial Monotonicity
Property 2.21 we get

######  Corollary 2.23.

For any given mass distribution $\vec{m}$, for any symbolic collision sequence
$\Sigma$, and for any (not necessarily contiguous) subsequence $\Sigma_{1}$ of
$\Sigma$, the sufficiency of the pair $(\Sigma_{1},\vec{m})$ implies the
sufficiency of $(\Sigma,\vec{m})$.

Thanks to 2.22, we obtain the typical, i.e. the mass distribution-free version
of 2.23:

######  Corollary 2.24.

For any symbolic collision sequence $\Sigma$ and for any (not necessarily
contiguous) subsequence $\Sigma_{1}$ of $\Sigma$, the sufficiency of
$\Sigma_{1}$ implies the sufficiency of $\Sigma$.

Next, we take a quick look at the limiting case $m_{1}=m_{2}=\dots=m_{n}$. In
this case, when a ball-to-ball collision $(i,i+1)$ occurs, we can simply swap
the labels $i$ and $i+1$ of the involved particles. In this way, by performing
all these dynamic label changes, the entire flow becomes the independent
motion of the $n$ particles, each of them falling at unit acceleration and
bouncing back elastically from the floor. (Doing so, they freely pass through
each other.)

Observe that, in the case of equal masses, the transformation matrices
$R_{i}^{*}$ in 2.14 are the reflections

$R_{i}=R_{i}^{*}=\left[\begin{array}[]{cc}0&1\\\ 1&0\end{array}\right],$

thus the components $\delta h_{i}(t)$ will be independent of the time $t$
after the dynamic re-labeling.

The above makes it clear that the trajectory segment $S^{[0,T]}x_{0}$ is
sufficient in the sense of 2.3, provided that all particles hit the floor.
(After the dynamic re-labeling of the particles, of course.) This means that,
for the limiting system $m_{1}=m_{2}=\dots=m_{n}$ there exists a uniform time
treshold $\tau_{0}>0$ such that $S^{[0,T]}x_{0}$ is sufficient, provided that
$T\geq\tau_{0}$. By continuity, this property extends to every falling ball
system $m_{1}>m_{2}>\dots>m_{n}>0$ with $m_{n}/m_{1}\geq 1-\epsilon_{0}$ for
some fixed $\epsilon_{0}>0$. However, this property, along with the CMP 2.24,
implies that there is a large enough $K\in\mathbb{Z}_{+}$ such that every
symbolic sequence $\Sigma=(\sigma_{1},\sigma_{2},\dots,\sigma_{N})$ is
sufficient if $\Sigma$ contains at least $K$ consecutive, contiguous,
connected subsequences. (Recall that a subsequence is connected when its
collision graph is connected or, equivalently, the subsequence contains all
types of collisions.) However, the system of falling balls obviously possesses
the property that every trajectory $S^{[0,\infty)}x_{0}$ contains infinitely
many appearances of each collision type. Thus we obtain the main result of
this section:

######  Theorem 2.6 (Main Result).

For almost every n-tuple of masses $m_{1}>m_{2}>\dots>m_{1}>0$ (the
exceptional set being a countable union of proper algebraic subsets of
$\mathbb{R}^{n}$) the following is true: Every non-singular positive orbit
$S^{[0,\infty)}x_{0}$ is sufficient, hence possesses the property of
eventually strict cone invariance.

##  3\. Conclusions

Here we summarize the corollaries of our Theorem 2.6. These corollaries will
be _conditional_ on the condition that the following property is possessed by
the falling ball system:

######  Claim 3.1 (Transversality Condition).

Singularities of different order are transversal to each other. Analogously,
stable and unstable invariant manifolds are transversal to all singularities.

Our goal is to prove the Ergodic Conjecture of Wojtkowski [W90a] for almost
every falling ball system, assuming the above transversality condition. The
exceptional systems are the ones for which our main result 2.6 is false, so
they form a countable union a proper algebraic sets in terms of the masses
$\vec{m}$. For this purpose we are going to use the celebrated Local
Ergodicity Theorem (LET, for short) of [L-W95], so we investigate the first
return map $T:\;U_{0}\to U_{0}$ of the billiard map to a suitably small, open
neighborhood of a hyperbolic phase point $x_{0}$ with at most one singularity
on its trajectory.

First of all, the second part of 3.1 is the so called Chernov-Sinai Ansatz, a
condition of the LET of [L-W95].

Second, the transversality of singularites guarantees that the set of double
singularities is a countable union of smooth, codimension-two ($\geq 2$)
submanifolds, hence a slim set, negligible in dynamical considerations, see
Section 2 of [K-S-Sz92], esp. Lemma 2.12 there. Thus, we may safely assume
that the considered hyperbolic phase point $x_{0}$ has at most one singularity
on its trajectory. Hyperbolicity means, of course, that the cone field $C_{1}$
is strictly invariant along the trajectory of $x_{0}$.

Third, the second part of 3.1 is Property 5’ of [Ch93], thus the result of
[Ch93] guarantees the Proper Alignment of Singularities, another condition of
the LET in [L-W95].

Since the first return map $T:\;U_{0}\to U_{0}$ enjoys strict $C_{1}$-cone
invariance, according to Proposition 6.1 of [L-W95], the minimum
$Q_{1}$-expansion rate $\sigma$ of $T$ is uniformly bigger than $1$. This
establishes the Non-Contraction Property and the Strict Unboundedness Property
for the first returm map $T:\;U_{0}\to U_{0}$, two more conditions of the LET.

Fourth, the Regularity of the Singular Set condition of [L-W95] is obtained as
follows:

  1. $(1)$

The phase points with two singularities on their trajectory form a _slim_ set,
as stated before.

  2. $(2)$

The smooth components of different singularities do not accumulate at any
point of the phase space, since there are no infinitely many collisions in
finite time (Proposition A.1 of the Appendix of [S96]), and the horizon is
finite.

  3. $(3)$

The derivatives of the singularities do not “explode”, since there are no
grazing (tangential) singularities in the falling ball system, only of corner
type.

Finally, the Local Ergodicity Theorem of [L-W95] states that the entire open
neighborhood of the hyperbolic phase point $x_{0}$ (with at most one
singularity on its trajectory) belongs to a single ergodic component of the
falling ball system. For almost every mass vector $\vec{m}$, the set of phase
points $x_{0}$ not having the two properties required above (strict cone
invariance and at most one singularity along its trajectory) is a slim set.
According to Lemma 2.12 of [K-S-Sz92], such slim sets are unable to cut the
phase space into different, open ergodic components. Thus, we conclude:

######  Theorem 3.1 (Main Corollary).

Assuming the condition 3.1, for almost every mass vector
$m_{1}>m_{2}>\dots>m_{n}>0$ the falling ball system is ergodic.

## References

  * [Ch93] Chernov, N.; Local ergodicity in hyperbolic systems with singularities, Functional Analysis & Applications, 27:1 (1993)
  * [L-W95] Liverani, C., Wojtkowski, M.; Ergodicity of Hamiltonian systems, Dynamics Reported 27:1, 130-202 (1995)
  * [K-S-Sz92] Krámli A., Simányi, N., Szász D.; The K-Property of Four Billiard Balls, Commun. Math, Phys. 144, 107-148 (1992)
  * [S96] Simányi, N.; The Characteristic Exponents of the Falling Ball Model, Commun. Math. Phys. 182, 457-468 (1996)
  * [W85] Wojtkowski, M.; Invariant families of cones and Lyapunov exponents, Ergodic Theory and Dynamical Systems, 5, 145-161 (1985)
  * [W86] Wojtkowski, M.; Principles for the Design of Billiards with Nonvanishing Lyapunov Exponents, Commun. Math. Phys. 105, 391-414 (1986)
  * [W90a] Wojtkowski, M.; A system of one dimensional balls with gravity, Commun. Math. Phys. 126, 507-533 (1990)
  * [W90b] Wojtkowski, M.; The system of one-dimensional balls in an external field. II, Commun. Math. Phys. 127, 425-432 (1990)
  * [W98] Wojtkowski, M.; Hamiltonian systems with linear potential and elastic constraints, Fundamenta Mathematicae, Volume 157, 305-341 (1998)