\sidecaptionvpos

figuret

¹¹institutetext: Department of Mechanical Engineering
Department of Computer Science
Yale University, New Haven CT 06511, USA
¹¹email: {henry.berger,ian.abraham}@yale.edu

RAnGE: Reachability Analysis
for Guaranteed Ergodicity

Henry Berger Ian Abraham

Abstract

This paper investigates performance guarantees on coverage-based ergodic exploration methods in environments containing disturbances. Ergodic exploration methods generate trajectories for autonomous robots such that time spent in each area of the exploration space is proportional to the utility of exploring in the area. We find that it is possible to use techniques from reachability analysis to solve for optimal controllers that guarantee ergodic coverage and are robust against disturbances. We formulate ergodic search as a differential game between the controller optimizing for ergodicity and an external disturbance, and we derive the reachability equations for ergodic search using an extended-state Bolza-form transform of the ergodic problem. Contributions include the computation of a continuous value function for the ergodic exploration problem and the derivation of a controller that provides guarantees for coverage under disturbances. Our approach leverages neural-network-based methods to solve the reachability equations; we also construct a robust model-predictive controller for comparison. Simulated and experimental results demonstrate the efficacy of our approach for generating robust ergodic trajectories for search and exploration on a 1D system with an external disturbance force.

1 Introduction

Autonomous mobile robots have been used for a number of search and exploration applications, including search and rescue [1, 2], mapping [3, 4], and extraterrestrial exploration [5]. One challenge is that many real-world applications of autonomous search occur in areas with hazards that impede search efforts. Hazards for autonomous mobile robots can be broadly categorized into two classes: static obstacles, such as walls or objects, which do not vary with time; and dynamic disturbances, such as wind and moving obstacles, which can vary over time. It is advantageous to find control solutions that provide guarantees on effective coverage, even in the presence of hazards. ^†^†For our code, see https://github.com/ialab-yale/RAnGE.

Recent advances in ergodic coverage-based exploration algorithms have demonstrated their efficacy for generating trajectories that provide coverage of a search area [6]. Ergodic exploration methods optimize trajectories such that over time, a robot’s trajectory spends time in each area of a search space proportional to the information to be gained in the area [7], providing sparser coverage in low-information areas and more meticulous coverage in high-information areas. However, current ergodic controllers have no formal performance guarantees in the presence of disturbances. Therefore, this paper investigates whether it is possible to provide performance guarantees on coverage-based ergodic control methods under dynamic disturbances.

Refer to caption — Figure 1: Overview of RAnGE. An agent explores an area with distributed information, subject to a disturbance (a), using a precomputed value function (b) to achieve comprehensive ergodic coverage (c). Note that in (b), higher values denote less desirable states.

Methods in reachability analysis offer a path toward providing such performance guarantees on ergodic control problems. Reachability analysis is a method for computing control policies with performance guarantees in worst-case disturbance scenarios [8, 9, 10]. Unfortunately, the prerequisite problem structure for using reachability analysis is not present in the canonical ergodic control problem formulation. In this paper, we investigate suitable formulations to construct reachability-based controllers for ergodic search that guarantee ergodic coverage of an area in the presence of dynamic external disturbances.

Our approach structures the ergodic control problem as a differential game between control and external disturbance force, and we employ a Bolza-form augmented-state transform of the ergodic metric that has also been used in other ergodicity-related work [7, 11, 12]. We show that the problem transform is equivalent to the original ergodic control problem and demonstrate compatibility with existing Hamilton-Jacobi-Isaacs (HJI) reachability methods [9]. We take advantage of computational advancements in HJI methods [13] to solve controllers capable of generating robust ergodic trajectories, although the limitations of these computational methods restrict the complexity of problems for which our controller can currently be computed. In summary, our main contributions are as follows:

•

Derivation of Hamilton-Jacobi-Isaacs reachability conditions for the problem of ergodic exploration;
•

Computation of a robust ergodic controller under dynamic disturbances;
•

Demonstration of the effectiveness of robust reachability-based ergodic controllers for 1D coverage on a double-integrator system (illustrated in Figure 1); and
•

Derivation of a robust ergodic model-predictive controller (reMPC) for ergodic exploration in the presence of disturbance.

This paper is structured as follows: Section 2 outlines the related work, Section 3 provides an overview of preliminary background information on ergodic exploration and reachability analysis, Section 4 derives the main contribution of this work, Section 5 demonstrates performance results of the proposed method, and Section 6 provides a discussion and conclusion. The appendices provide implementation details, proofs, and an additional derivation.

2 Related Work

2.1 Reachability Analysis

Robots operating in the real world are often subject to disturbances, so it is desirable to have guarantees of performance for robotic systems in worst-case conditions. Reachability analysis is a verification technique for robotic systems under such performance scenarios. For a bounded, time-varying disturbance, these analysis methods can determine the worst-case disturbance and provide a lower bound on performance in all possible conditions. Reachability has already been used for a variety of autonomous robotics applications, including control and obstacle avoidance for drones [10, 14], collision avoidance for airplanes [13, 15], and safe navigation for autonomous vehicles [13, 14].

Reachability analysis requires finding the viscosity solution to the Hamilton-Jacobi-Isaacs (HJI) partial differential equation [8, 16]; for most applications, this cannot be done analytically and must be done numerically. For any complex problem where an analytical solution is not feasible, the numerical solution is only an approximation of the value function. In some cases, it is possible to guarantee that all inaccuracies are conservative, which maintains the performance guarantees [10, 17]. In other cases, the reachability-based controllers cannot quite provide the formal guarantees promised by reachability analysis, but as the accuracy of the numerical solution improves, a violation of the guarantees becomes less likely.

Much research on reachability analysis has focused on low/dimensional systems. Additionally, due to the difficulty of learning the value function, reachability is often used modularly to provide safety guarantees (e.g., collision avoidance), without providing performance guarantees relating to the robot’s main objective [18, 19]. In this work, we generate approximate solutions to the HJI problem for the case of ergodic exploration, a high-dimensional problem where reachability is used to guarantee performance, not merely safety.

The Curse of Dimensionality.

The HJI equation is often solved through discretizing the state space on a grid. However, the computational cost of this approach scales exponentially with the problem state dimensionality, and as a result, grid-based solutions to the HJI equation are limited to low/dimensional problems. State decomposition methods [10, 20] can mitigate this scaling, but many high-dimensional problems are still difficult to solve on a grid. Furthermore, decomposition methods compute individual level sets of the value function, not a complete, differentiable value function.

Recent advances using neural networks for solving differential equations [21] have provided new methods for solving the reachability equations [17, 22, 13]. For a neural network, the computational time required to find a satisfactory solution depends on the complexity of the solution, rather than strictly on the dimensionality of the state space. Higher-dimensional problems will often have more complex value functions and therefore take longer to learn, but the scaling is better than exponential [13]. Neural-network-based solvers have demonstrated solutions of the reachability equations for 10D problems with disturbance [13] and 30D problems without disturbance [22].

2.2 Ergodic Exploration

In this paper, we apply reachability to a quintessentially trajectory-based optimization problem: coverage-based exploration. In general, the goal of exploration is to maximize coverage of a search area, and for problems involving a single mobile robot, coverage is inherently a property of an entire trajectory, not of an instantaneous state.

For discretized exploration problems, a variety of methods already exist to provide formal guarantees on performance. In discrete time and space, exploration can be represented as a traveling salesperson problem, enabling the application of fast, locally optimal algorithms [23] or slower, formally guaranteed globally optimal algorithms [24]. Voronoi methods can be used to guarantee maximal coverage of a continuous exploration area by a set of agent locations [25, 26]; these locations can be treated as the static locations of a swarm of agents or as the discretized trajectory of a single agent, but Voronoi methods cannot generate continuous trajectories. The motions of most robotic systems are continuous in time and space, so discrete-time and discrete-space methods necessarily introduce simplifications, which can lead to non-optimal results.

Ergodic exploration is a commonly-used technique for continuous-time, continuous-space exploration. Ergodic controllers seek to minimize the difference between the (expected) spatial distribution of information and the time-averaged spatial distribution of an agent’s trajectory [7]. Dressel and Kochenderfer [6] proved that ergodic exploration is an optimal strategy for gathering information if, as is often the case in real-world applications, the amount of information in a location is depleted as the agent spends time in that location.

Ergodic exploration has been combined with formally guaranteed static obstacle avoidance, using methods such as vector fields [27] and control barrier functions [28]. In this work, we introduce formal performance guarantees with time-varying disturbances into ergodic exploration, broadening the set of situations in which ergodic search can safely and effectively be used.

3 Problem Setup and Definitions

We consider a general system with state $x\in\mathcal{X}\subseteq\mathbb{R}^{n}$ , control input $u\in\mathcal{U}\subseteq\mathbb{R}^{n_{u}}$ , disturbance input $d\in\mathcal{D}\subseteq\mathbb{R}^{n_{d}}$ , and (possibly nonlinear) dynamics

\dot{x}(t)=f(x(t),u(t),d(t)).

(1)

In this work, we consider the system over a finite time interval $t\in[{t_{0}},{t_{f}}]$ . The control input is generated by a control policy $u(\cdot):\mathcal{X}\times[{t_{0}},{t_{f}}]\to\mathcal{U}$ , also called the controller; the disturbance input is generated by a disturbance policy $d(\cdot):\mathcal{X}\times[{t_{0}},{t_{f}}]\times\mathcal{U}\to\mathcal{D}$ , also called the disturbance.

We assume that the disturbance has an instantaneous information advantage, meaning that the disturbance can respond instantaneously to the controller, although it cannot anticipate the controller’s future actions [9]. More concretely, at any time $t$ , the control is a function of $x$ and $t$ only, whereas the disturbance can be a function of $x$ , $t$ , and $u$ .

3.1 Ergodic Exploration

Ergodicity is a measure of the effectiveness of a trajectory in covering distributed information in a known bounded space. We consider exploration of a $v$ -dimensional space $\mathcal{M}=[0,L_{0}]\times[0,L_{1}]\times\cdots[0,L_{v-1}]\subset\mathbb{R}^{v}$ , where $L_{i}$ are bounds on the exploration space. We define a function $M:\mathcal{X}\to\mathcal{M}$ that maps state $x\in\mathcal{X}$ to a point in the exploration space $m\in\mathcal{M}$ .¹¹1Often, the exploration task only involves a subset of the agent’s state (e.g., pose), though the state includes additional variables (e.g., velocity). The mapping $M$ between $\mathcal{X}$ and $\mathcal{M}$ accounts for this distinction; for example, $M$ could be a selection matrix that selects only the pose components of the state.

Definition 3.1 (Ergodicity).

Given a function $\phi:\mathcal{M}\to\mathbb{R}_{>0}$ , an ergodic trajectory $x(\cdot):[{t_{0}},{t_{f}}]\to\mathcal{X}$ is one that spends time in each region of $\mathcal{M}$ in proportion to the integral of $\phi$ over that region, such that

\lim_{t\to\infty}{\frac{1}{t-{t_{0}}}\int_{{t_{0}}}^{t}{F(M\circ x(\tau))d\tau}}=\int_{\mathcal{M}}{F(m)\phi(m)dm}

(2)

for all Lebesgue integrable functions $F:\mathcal{M}\to\mathbb{R}$ [7].

The definition of $\phi$ is problem-specific. For search and exploration applications, $\phi$ is an information density function that represents the probability of discovering information at a location. Following [7], ergodicity is often measured by a comparison in the Fourier domain between the information density $\phi$ and the time-averaged trajectory statistics, which are defined as follows.

Definition 3.2 (Time-averaged trajectory statistics).

The time-averaged location statistics of the trajectory are given by

c(x(\cdot),t)\coloneqq\frac{1}{t-{t_{0}}}\int_{{t_{0}}}^{t}\delta[M(x(\tau))]d\tau,

(3)

where $\delta$ is the Dirac delta function.

In order to measure whether a trajectory is ergodic, we define a metric on ergodicity.

Definition 3.3 (Spectral ergodic metric).

Following [7], the ergodic metric is given by

\mathcal{E}(x(\cdot),t)\coloneqq\sum_{k\in\mathcal{K}}{\Lambda_{k}\left(c_{k}(x(\cdot),t)-\phi_{k}\right)^{2}},

(4)

where the time-average trajectory statistics $c$ and the information density function $\phi$ are encoded in the coefficients $c_{k}$ and $\phi_{k}$ , respectively, calculated using the Fourier transform

c_{k}(x(\cdot),t)\coloneqq\frac{1}{t-{t_{0}}}\int_{{t_{0}}}^{t}{F_{k}(M\circ x(\tau))d\tau},\ \ \phi_{k}\coloneqq\int_{\mathcal{M}}{F_{k}\left(m\right)\phi(m)dm},

(5)

for modes $k\in\mathcal{K}=\{0,1,2,...,N-1\}^{v}$ . Here, the set $\{\Lambda_{k}\in\mathbb{R}\}$ is a sequence of weighting coefficients defined for each $k\in\mathcal{K}$ as $\Lambda_{k}\coloneqq(1+\lVert k\rVert_{2}^{2})^{-s}$ , where $s\coloneqq(v+1)/2$ .

We use the cosine basis function $F_{k}(m)\coloneqq\frac{1}{h_{k}}\prod_{i=0}^{v-1}{\cos{\left(\frac{k_{i}\pi m_{i}}{L_{i}}\right)}},$ where $\frac{1}{h_{k}}$ is a normalizing constant [7], but other orthonormal bases would also suffice. The ergodic metric measures how far a trajectory is from ergodicity. The canonical ergodic trajectory-optimization problem is the minimization of the ergodic metric with respect to a given information metric $\phi$ [7, 29] (see Problem 2 below).

3.2 Reachability Analysis

Reachability analysis is a technique for solving differential games via the Hamilton-Jacobi-Isaacs (HJI) partial differential equation, which can be used to model safety requirements in control systems. We first consider a problem in which, subject to dynamic constraints, the controller and disturbance seek to minimize and maximize, respectively, a cost function with Bolza form (Problem 1).

Problem 1 (Bolza-form optimal control with disturbance)

Given a terminal cost $h:\mathcal{X}\to\mathbb{R}_{\geq 0}$ and a running cost $g:\mathcal{X}\times\mathcal{U}\to\mathbb{R}_{\geq 0}$ , for some choice of $x_{0}$ and ${t_{0}}$ , find

	$u^{*}(\cdot)\coloneqq\arg\min_{u(\cdot)}\max_{d(\cdot)}\left\{J(x_{0},{t_{0}},u(\cdot),d(\cdot))\right\},$		(6a)
where the cost function $J$ is defined as
	$J\left(x_{0},{t_{0}},u(\cdot),d(\cdot)\right)\coloneqq h(x({t_{f}}))+\int_{{t_{0}}}^{{t_{f}}}{g(x(\tau),u(\tau))d\tau},$		(6b)
subject to the following conditions ( $\forall t\in[{t_{0}},{t_{f}}]$ , where applicable):

	$\displaystyle u(t)$	$\displaystyle=u(x(t),t)$	(6c)
	$\displaystyle d(t)$	$\displaystyle=d(x(t),t,u(t))$	(6d)
	$\displaystyle\dot{x}(t)$	$\displaystyle=f(x(t),u(t),d(t))$	(6e)
	$\displaystyle x({t_{0}})$	$\displaystyle=x_{0}.$	(6f)

Definition 3.4 (Value function).

The value function $V:\mathcal{X}\times[{t_{0}},{t_{f}}]\to\mathbb{R}$ , also called the optimal cost-to-go function, is the cost of the remainder of the trajectory, starting at the given time and state, ending at time ${t_{f}}$ , with optimal control and disturbance.

V(x,t)=\min_{u(\cdot)}\max_{d(\cdot)}\left\{J(x,t,u(\cdot),d(\cdot)\right)\}.

(7)

As proved in Theorem 4.1 of [16], the value function for a Bolza-form optimization problem (Problem 1) is a viscosity solution to the Hamilton-Jacobi-Isaacs optimality conditions [16].

Definition 3.5 (Hamilton-Jacobi-Isaacs optimality conditions).

For the Bolza-form optimization problem (Problem 1), the value function $V:\mathcal{X}\times[{t_{0}},{t_{f}}]\to\mathbb{R}$ satisfies the conditions
	$-\frac{\partial V}{\partial t}=H(x,t),\ \forall t\in[{t_{0}},{t_{f}}],\forall x\in\mathcal{X}$	(8a)
	$V(x,{t_{f}})=h(x),\ \forall x\in\mathcal{X},$	(8b)
where the Hamiltonian $H$ is defined as
	$H(x,t)\coloneqq\min_{u\in\mathcal{U}}\max_{d\in\mathcal{D}}{\Bigr{\{}g(x,u)+\left(\frac{\partial V(x,t)}{\partial x}\right)^{\intercal}f(x,u,d)\Bigl{\}}}.$	(8c)

Once value function $V$ is known, the HJI equation (8a) and the definition of the Hamiltonian (8c) entail that the optimal controller $u^{*}(\cdot)$ and worst-case disturbance $d^{*}(\cdot)$ can be calculated directly as


$\displaystyle u^{*}(x,t)=\arg\min_{u\in\mathcal{U}}\bigg{\{}$	$\displaystyle\max_{d\in\mathcal{D}}{\Bigr{\{}\hat{H}(x,t,u,d)\Bigl{\}}}\bigg{\}}$	(9a)
$\displaystyle d^{*}(x,t,u)=\arg$	$\displaystyle\max_{d\in\mathcal{D}}{\Bigr{\{}\hat{H}(x,t,u,d)\Bigl{\}}},$	(9b)

where $\hat{H}$ represents the target of the minimax condition in the Hamiltonian $H$ [30]. This optimal controller is one of the main results of reachability analysis, since the controller can be used as a robust policy against worst-case disturbance. We seek to find such an optimal controller for exploration problems, in order to provide sufficient coverage guarantees for exploration with disturbance.

4 Reachability Analysis for Guaranteed Ergodic Exploration

We now consider the problem of ergodic control subject to disturbance, as formulated in Problem 2.

Problem 2 (Ergodic control with disturbance)

Solve an optimization problem of the same form as the Bolza optimization problem with disturbance (Problem 1), but with the modified cost function

J\left(x_{0},{t_{0}},u(\cdot),d(\cdot)\right)\coloneqq\mathcal{E}(x(\cdot),{t_{f}})+h(x({t_{f}}))+\int_{{t_{0}}}^{{t_{f}}}{g(x(\tau),u(\tau))d\tau}.

(10)

This problem is not strictly compatible with Bolza form, because the ergodic cost depends on the entire trajectory $x(\cdot)$ , not just on the terminal state $x({t_{f}})$ . Consequently, we reformulate the ergodic control problem in order to apply reachability analysis. To do so, we first derive an equivalent form of the ergodic metric via an augmented-state ergodic variable.

Lemma 4.1 (Augmented-state ergodic metric formulation [7, 11]).

The ergodic metric defined in (4) can be equivalently defined as

	$\displaystyle\mathcal{E}(x(\cdot),{t_{f}})$	$\displaystyle=\sum_{k\in\mathcal{K}}{\left(\Lambda_{k}\left(c_{k}(x(\cdot),t)-\phi_{k}\right)^{2}\right)}$
		$\displaystyle=\frac{1}{t_{f}-t_{0}}\sum_{k\in\mathcal{K}}{\left(\Lambda_{k}z_{k}({t_{f}})^{2}\right)}$		(11)

where $z_{k}(\cdot)$ is the solution to the differential equation


$\displaystyle\dot{z_{k}}(t)$	$\displaystyle=F_{k}(M\circ x(t))-\phi_{k}$	(12a)
$\displaystyle z_{k}({t_{0}})$	$\displaystyle=0\ \forall k\in\mathcal{K}.$	(12b)

Proof 4.2.

See Appendix B.1. ∎

This formulation of the ergodic metric (11) can be used to rewrite the ergodic control problem with disturbance (Problem 2) in Bolza form.

Problem 3 (Bolza-form ergodic control with disturbance)

Given a terminal cost $h:\mathcal{X}\to\mathbb{R}_{\geq 0}$ and a running cost $g:\mathcal{X}\times\mathcal{U}\to\mathbb{R}_{\geq 0}$ , for some choice of $x_{0}$ and ${t_{0}}$ , find

	$u^{*}(\cdot)\coloneqq\arg\min_{u(\cdot)}\max_{d(\cdot)}\left\{J(x_{0},{t_{0}},u(\cdot),d(\cdot))\right\},$			(13a)
where the cost function $J$ is defined as

	$\displaystyle J\left({t_{0}},x_{0},u(\cdot),d(\cdot)\right)=$	$\displaystyle\frac{1}{t_{f}-t_{0}}\sum_{k\in\mathcal{K}}{\left(\Lambda_{k}z_{k}({t_{f}})^{2}\right)}+h(x({t_{f}}))$
		$\displaystyle+\int_{{t_{0}}}^{{t_{f}}}{g(x(\tau),u(\tau))d\tau},$		(13b)
subject to the following conditions ( $\forall t\in[{t_{0}},{t_{f}}]$ and $\forall k\in\mathcal{K}$ , where applicable):

	$\displaystyle u(t)$	$\displaystyle=u(x(t),z(t),t),$	$\displaystyle d(t)=d(x(t),z(t),t,u(t)),$	(13c)
	$\displaystyle\dot{x}(t)$	$\displaystyle=f(x(t),u(t),d(t)),$	$\displaystyle x({t_{0}})=x_{0},$	(13d)
	$\displaystyle\dot{z}_{k}(t)$	$\displaystyle=F_{k}(M\circ x(t))-\phi_{k},\hskip 18.06749pt$	$\displaystyle z_{k}({t_{0}})=0.$	(13e)

Note that $z\in\mathcal{Z}=\mathbb{R}^{N^{v}}$ is a vector of $z_{k},\forall k\in\mathcal{K}$ .

Problem 3 is a canonical Bolza-form problem (Problem 1), which enables the application of reachability analysis.

Theorem 4.3 (Ergodic Hamilton-Jacobi-Isaac PDE).

The solution to the ergodic control problem with disturbance (Problem 3) is given by the value function $V:\mathcal{X}\times\mathcal{Z}\times[{t_{0}},{t_{f}}]\to\mathbb{R}$ , optimal controller $u^{*}(\cdot):\mathcal{X}\times\mathcal{Z}\times[{t_{0}},{t_{f}}]\to\mathcal{U}$ , and optimal disturbance $d^{*}(\cdot):\mathcal{X}\times\mathcal{Z}\times[{t_{0}},{t_{f}}]\times\mathcal{U}\to\mathcal{D}$ such that $V$ satisfies


	$\displaystyle-\frac{\partial V(x,z,t)}{\partial t}$	$\displaystyle=\min_{u\in\mathcal{U}}\max_{d\in\mathcal{D}}\left\{\hat{H}(x,z,t,u,d)\right\},\text{where}$	(14a)
	$\displaystyle\hat{H}(x,z,t,u,d)$	$\displaystyle=\bigg{\{}g(x,u)+\left(\frac{\partial V(x,z,t)}{\partial x}\right)^{\intercal}f(x,u,d)$
		$\displaystyle\phantom{=\bigg{\{}}+\sum_{k\in\mathcal{K}}{\left(\left(\frac{\partial V(x,z,t)}{\partial z_{k}}\right)^{\intercal}\left(F_{k}(M(x))-\phi_{k}\right)\right)}\bigg{\}},$	(14b)
with terminal condition
	$V(x,z,{t_{f}})=\frac{1}{{t_{f}}-{t_{0}}}\sum_{k\in\mathcal{K}}{\left(\Lambda_{k}z_{k}({t_{f}})^{2}\right)}+h(x),$		(14c)

and such that $u^{*}(\cdot)$ and $d^{*}(\cdot)$ satisfy


$\displaystyle u^{*}(x,z,t)$	$\displaystyle=\arg\min_{u\in\mathcal{U}}\max_{d\in\mathcal{D}}\left\{\hat{H}(x,z,t,u,d)\right\},$	(15a)
$\displaystyle d^{*}(x,z,t,u)$	$\displaystyle=\arg\max_{d\in\mathcal{D}}\left\{\hat{H}(x,z,t,u,d)\right\}.$	(15b)

Proof 4.4.

(Informal) The proof follows [31] and Section 4.2 of [30]. Using the definition of the value function (Def. 3.4) with the ergodic cost function (13b), the value function at the terminal time can be shown to satisfy (14c). We then break the value function into Isaacs’ canonical principle of optimality (Eq. 2.7 of [31]), which is a minimax variant of the Bellman principle of optimality. From there, following the standard reachability derivation [30, 31] yields the remaining conditions for $V,u^{*},d^{*}$ . A more detailed proof is provided in Appendix B.2. ∎

We additionally prove in Appendix C that the canonical and extended-state formulations of the ergodic control problem (Problems 2 and 3, respectively) yield the same HJI equation (14). This demonstrates that the presence of extended states in the HJI equation (14) is a feature of the problem structure, not an artifact of the choice of problem formulation.

The optimal controller $u^{*}$ from Theorem 4.3 is robust in the sense that it satisfies the following lemmas, which are direct consequences of (13a).

Remark 4.5 (Cost guarantees).

Applying $u^{*}(\cdot)$ starting at $(x,z,{t_{0}})$ guarantees that despite any possible disturbance $d\in\mathcal{D}$ , the trajectory cost will be upper-bounded by $V(x,z,{t_{0}})$ ; furthermore, no controller can guarantee an upper bound lower than $V(x,z,{t_{0}})$ .

Remark 4.6 (Reachable coverage set).

If $g(x,t)=0$ and $h(x)=0$ , then the set $\{(x,z)\in\mathcal{X}\times\mathcal{Z}\ \rvert\ V(x,z,{t_{0}})\leq\epsilon\}$ is precisely the set of starting states from which the optimal controller can guarantee achieving an ergodic metric of at most $\epsilon$ .

We solve the HJI equation using a deep neural network (DNN), following [13], and we call the resulting controller RAnGE. Implementation details are provided in Appendix A.

4.1 Practical Considerations Regarding Scalability and Robustness

The HJI PDE (14) can only be solved numerically, which raises two concerns. Firstly, the solved value function is only an approximation, and errors in the approximation could invalidate the guarantees of Theorem 4.3. This is a common issue for reachability methods; one common solution is to make the value function overly conservative, which is possible even when using neural networks [17], but this compromises the accuracy of the value function. We do not provide such conservative guarantees here, but as numerical HJI solvers improve, any approximation errors and their consequences should diminish.

Secondly, scalability is a major limitation for HJI solvers, and the extra dimensionality from the augmented state $z$ limits the scope of ergodic problems that can be solved. The number of added dimensions is $N^{v}$ , which scales polynomially with respect to the resolution $N$ of the basis functions and exponentially with respect to the dimensionality $v$ of the exploration space $\mathcal{M}$ . The curse of dimensionality is common to all reachability applications, and developing new solvers with better scalability is an active area of research.

Alternative formulations of reachability for ergodic search may be possible, but we do not anticipate that any such alternative formulation could avoid the issue of increased dimensionality. The ergodic metric depends on the history of the trajectory, so an optimal controller must have a way to respond to information about the history of the trajectory up to the current time. Whether this information is stored in the augmented state variable $z$ (as in this paper, as well as [7, 11, 12]) or some other way, it will necessarily constitute an additional input to the value function, thus increasing the dimensionality of the problem. Indeed, we show in Appendix C that extended state variables emerge from the application of reachachability analysis to ergodic exploration, even when using an ergodic formulation without extended states (Problem 2).

4.2 An Alternative Robust Ergodic Model-Predictive Controller

We additionally develop a robust ergodic model-predictive controller, which we call reMPC. This controller uses first-order gradient descent to solve a discretized version of the minimax condition of the Bolza-form ergodic optimization problem (13):²²2 Here, the dynamics are discretized using the explicit Euler method, but any numerical integration method could be used.

Problem 4 (Discretized Bolza-form ergodic control with disturbance)

Given a start state $(x\in\mathcal{X},z\in\mathcal{Z})$ , for some timestep $\Delta t\in\mathbb{R}_{>0}$ and time horizon $T\in\mathbb{N}$ , solve for the optimal control $u^{*}\in\mathcal{U}^{T}$ such that

\displaystyle\mathbf{u}^{*}=\arg\min_{\mathbf{u}\in\mathcal{U}^{T}}\max_{\mathbf{d}\in\mathcal{D}^{T}}\bigg{\{}\sum_{\tau=0}^{T}{g(\mathbf{x}_{\tau},\mathbf{u}_{\tau})\Delta t}+\frac{1}{T\Delta t}\sum_{k\in\mathcal{K}}{\left(\Lambda_{k}(\mathbf{z}_{T})_{k}^{2}\right)}+h(\mathbf{x}_{T})\bigg{\}},

(16)

subject to


$\displaystyle\mathbf{x}_{\tau}$	$\displaystyle=\mathbf{x}_{\tau-1}+f(\mathbf{x}_{\tau-1},\mathbf{u}_{\tau-1},\mathbf{d}_{\tau-1})\Delta t,\ \forall\tau\in[1,2,\dots,T]$	(17a)
$\displaystyle(\mathbf{z}_{\tau})_{k}$	$\displaystyle=(\mathbf{z}_{\tau-1})_{k}+(F_{k}(M(\mathbf{x}_{\tau-1})-\phi_{k})\Delta t,\ \forall\tau\in[1,2,\dots,T],\forall k\in\mathcal{K}$	(17b)
$\displaystyle(\mathbf{x}_{0},\mathbf{z}_{0})$	$\displaystyle=(x,z).$	(17c)

The complete reMPC algorithm is given as Alg. 1. The first-order gradient descent method (lines 1-1) produces only local optima, so reMPC has no performance guarantees. However, given a long enough time horizon, if reMPC finds solutions close to the global optimum, then its performance should be similar to that of RAnGE.

Input: Starting state

(x_{0}\in\mathcal{X},z_{0}\in\mathcal{Z})

Input: Dynamics

f:\mathcal{X}\times\mathcal{U}\times\mathcal{D}\to\mathcal{T}_{\mathcal{X}}

Input: Basis functions

\{F_{k}:\mathcal{M}\to\mathbb{R},\forall k\in\mathcal{K}\}

, information

\{\phi_{k}\in\mathbb{R},\forall k\in\mathcal{K}\}

Input: Running cost

g:\mathcal{X}\times\mathcal{U}\to\mathbb{R}_{\geq 0}

, terminal cost

h:\mathcal{X}\to\mathbb{R}_{\geq 0}

1 function reMPC(

x_{0}\in\mathcal{X},z_{0}\in\mathcal{Z}

)

\mathbf{u},\mathbf{z}\leftarrow\mathbf{0}

3 while $\mathrm{controller\ is\ running}$ do

\mathbf{u}_{0,\dots,T-2}\leftarrow\mathbf{u}_{1,\dots,T-1}

\mathbf{u}_{T-1}\leftarrow\mathbf{0}

\mathbf{d}_{0,\dots,T-2}\leftarrow\mathbf{d}_{1,\dots,T-1}

\mathbf{d}_{T-1}\leftarrow\mathbf{0}

6 for N $\ \mathrm{iterations}$ do

\mathbf{u}\leftarrow\mathbf{u}-\lambda_{u}\cdot\left(\partial\text{\sc{PredRollout}}/\partial\mathbf{u}\right)

\mathbf{d}\leftarrow\text{clip}\left(\mathbf{d}+\lambda_{d}\cdot\left(\partial\text{\sc{PredRollout}}/\partial\mathbf{d}\right),\mathcal{D}\right)

\triangleright

Constrain

d\in\mathcal{D}

10 Apply

\mathbf{u}_{0}

13function

\textsc{PredRollout}{(}x_{0}\in\mathcal{X},z_{0}\in\mathcal{Z},\mathbf{u}\in\mathcal{U}^{T},\mathbf{d}\in\mathcal{D}^{T}{)}

\mathbf{x}_{0}\leftarrow x_{0},\ \mathbf{z}_{0}\leftarrow z_{0}

15 for $\tau\mathrm{\ in\ }[1,2,\dots,T]$ do

\mathbf{x}_{\tau}\leftarrow\mathbf{x}_{\tau-1}+f(\mathbf{x}_{\tau-1},\mathbf{u}_{\tau-1},\mathbf{d}_{\tau-1})\Delta t

(\mathbf{z}_{\tau})_{k}\leftarrow(\mathbf{z}_{\tau-1})_{k}+(F_{k}(M(\mathbf{x}_{\tau-1})-\phi_{k})\Delta t,\ \forall k\in\mathcal{K}

19 return

\sum_{\tau=0}^{T}{g(\mathbf{x}_{\tau},\mathbf{u}_{\tau})\Delta t}+\frac{1}{T\Delta t}\sum_{k\in\mathcal{K}}{\left(\Lambda_{k}(\mathbf{z}_{T})_{k}^{2}\right)}+h(\mathbf{x}_{T})

Algorithm 1 Robust Ergodic Model-Predictive Control (reMPC)

5 Results

In this section, we demonstrate the RAnGE controller for the double-integrator system moving along one axis.

Due to the computational complexities of the value function, we found that training the value function was only possible for a one-second time horizon. The results in this section were achieved by running the controller with a receding one-second time horizon, always inputting $t\approx{t_{0}}$ into the controller.³³3Because $t={t_{0}}$ was on the boundary of the time interval over which the DNN was trained, a slightly larger time ( $t=0.02s$ ) in the interior of the time interval was found to provide better behavior. This modification violates the preconditions for the guarantees of Theorem 4.3 to apply, but if the HJI equation could be solved over a longer time horizon, the guarantees would apply over that time horizon.

We generate trajectories from RAnGE and compare against reMPC (Section 4.2). We also compare against a standard model-predictive controller (MPC), a type of controller that does not explicitly account for disturbance but has often been shown to perform well under external disturbance [27, 32]. We implement a conventional MPC that minimizes the Bolza-form ergodic cost function (13b) without considering disturbance.

The goal for this section is to answer the following questions:

Q1:

Can RAnGE generate robust ergodic trajectories in the presence of a disturbance?
Q2:

How does the value function evolve over the course of a trajectory?
Q3:

How does the performance of RAnGE compare to the performance of MPC and reMPC?
Q4:

Regarding its effects on controller performances, how does the worst-case disturbance generated by RAnGE compare to other types of disturbance?

The following subsections provide an overview of the problem setup and answers to the questions stated above; each subsection A $n$ in 5.2 below answers the corresponding question Q $n$ above.

5.1 Problem Setup: One-Axis Double Integrator

We use a single-axis double-integrator dynamical system to study the proposed method. The system has state $x=[x_{1}\ x_{2}]^{\intercal}$ , control $u=\left[u\right]$ , disturbance $d=\left[d\right]$ , and dynamics $\dot{x}=[x_{2},u+d]^{\intercal}$ . Here, we subject the system to the control and disturbance constraints $u\in\left[-u_{\text{max}},u_{\text{max}}\right]\subset\mathbb{R},\ d\in\left[-d_{\text{max}},d_{\text{max}}\right]\subset\mathbb{R}$ . For all results in this paper, we use $d_{\text{max}}=0.4u_{\text{max}}$ ; performance was found to decline significantly for larger disturbance magnitudes.

We set $N=6$ , leading to an 8-dimensional system: one dimension for time, two dimensions for the double integrator state, and five dimensions for the augmented state.⁴⁴4 $N=6$ only necessitates five augmented state variables, because $z_{k=\textbf{0}}(t)=0$ for all $t$ , so there is no need to compute and track this term. Note that while the plant dynamics are linear, the dynamics of the augmented state variables (12) are nonlinear, so the augmented system as a whole is nonlinear.

The ergodic problem setup was that of (14), with running cost

g(x,z,u)\coloneqq q\sum_{k\in\mathcal{K}}{\left(\Lambda_{k}z_{k}^{2}\right)}+u^{\intercal}Ru+\text{barr}(x),

(18)

where $q$ is a nonnegative scalar weighting factor, $R\in\mathbb{R}^{1\times 1}$ is a positive semi-definite matrix representing a weight on the control input, and $\text{barr}(x)$ is a barrier function that penalizes the trajectory for leaving the search area bounds on position. We introduce the augmented state $z$ into the running cost with the running ergodic cost term $q\sum_{k\in\mathcal{K}}{\left(\Lambda_{k}z_{k}^{2}\right)}$ , which was found to improve convergence by helping the solver avoid a local minimum where $V=0$ everywhere except close to $t={t_{f}}$ . Other than the ergodic metric, there is no additional terminal cost $h(x)$ .

5.2 Simulated Results

5.2.1 A1: Robust ergodic trajectories.

A robust controller must produce ergodic trajectories even with worst-case disturbance, and with the controller from RAnGE, it is possible to calculate this worst-case disturbance (15b). In this particular problem setup, the control and disturbance are both affine, and they enter the dynamics in exactly the same way; it follows that the worst-case disturbance is always in the opposite direction of the optimal control. Thus, the worst-case disturbance can be calculated as

d^{*}=(-\mathrm{sign}(u^{*}))d_{\mathrm{max}}.

(19)

We demonstrate that the controller can generate robust ergodic trajectories with worst-case disturbance, as shown in Figure 2. More trajectory results are summarized in Figure 4, discussed below.

5.2.2 A2: Value function evolution.

Figure 3 shows the evolution of a trajectory, as well as a cross-section of the value function at times during the trajectory. Qualitatively, properties of the value function can be related to the problem setup. Recall that the controller seeks states where the value function is lowest. At the top right and lower left corners, the value function is high, indicating that the controller should avoid these states; this is logical, because the agent is about to leave the exploration bounds and incur a boundary cost. Elsewhere, the value function is lowest in states where the agent is at or moving toward positions it has not adequately explored, and highest in states where the agent is at or moving toward positions it has already over-explored.

5.2.3 A3: Comparison of RAnGE, MPC and reMPC.

The performance of the three controllers for two information distributions and four disturbance types is summarized in Figure 4. There is considerable variation in the performances within each trial, due in part to the limitations of the controllers’ short time horizons, but some trends emerge.

As expected, the worst performances by RAnGE are better than the worst performances by either of the other two methods. This is the primary theoretical advantage of RAnGE: it should have the highest floor for performance of any controller. By contrast, the best performance by the other controllers rival or exceed the performance of RAnGE, particularly with non-worst-case noise. This is consistent with the formulation of RAnGE, which is designed to optimize its worst-case performance, not its best-case performance. In situations where average performance with non-worst-case disturbance is most important, the MPC could be a better choice.

Interestingly, the reMPC controller performs its best—both absolutely and relative to the other methods—with worst-case disturbance. This is because the reMPC controller optimizes its controls under the assumption of a near-worst-case disturbance, so it is not optimized for non-worst-case disturbances. Consequently, reMPC outperforms MPC with worst-case disturbance but underperforms it otherwise.

5.2.4 A4: Comparison of Disturbances

Just as the RAnGE controller optimizes its worst-case performance, the RAnGE-generated disturbance is optimally designed to impede best-case performance by any controller. Across both information distributions, the best performance by any controller with worst-case disturbance is worse than the best performance by any controller with each of the other types of noise.

The effects of the two types of random disturbance on controller performance are very similar to the effect of no disturbance at all. This is unsurprising, because the random disturbances average to zero and thus tend to cancel out.

5.3 Experimental Results

We tested the controller on a Bitcraze Crazyflie 2.1 quadrotor drone [33] with fans for disturbance; the experimental setup is shown in Figure 5a. The drone was tested with no disturbance; with the fans at constant low speed; and with the fans varying between low, medium, and high speeds. Trajectories from these tests are shown in subfigures b-d of Figure 5. We were unable to implement a worst-case disturbance in the experiment because the worst-case disturbance must react instantaneously to the trajectory of the drone.

In all cases, the drone executed ergodic trajectories, as shown with the reconstructed time-averaged statistics in Figure 5. This demonstrates that the RAnGE controller is robust to disturbance and can be implemented on physical systems. The trajectory was most ergodic in the case of no disturbance, which is consistent with the results of Figure 4.

6 Discussion

In conclusion, we showed that it is possible to guarantee coverage using reachability analysis on ergodic search problems despite dynamic disturbances. We demonstrated a Bolza-form formulation for ergodic exploration using an augmented state that allows us to derive the ergodic HJI conditions of optimality. To compute solutions to the ergodic HJI conditions, we leveraged recent deep learning methods to approximate the value function. The controller was demonstrated to provide robust ergodic trajectories in the presence of dynamic disturbances, both in simulation and in experiments.

One limitation of this work is that due to the value function complexity, the HJI equation could not be solved for more complex problems (e.g., time horizons much greater than one second and exploration spaces with multiple dimensions). This limitation is a result of numerical implementation and not the fundamental theory. With advances in HJI solving techniques, future work could extend these results to longer time horizons, enabling longer guarantees on ergodic coverage trajectories. Furthermore, we seek to extend this work to higher, multidimensional search spaces and introduce dynamic obstacles in the formulation. This paper provides the theoretical foundation for performing reachability analysis on ergodic exploration, and improved solving techniques will broaden its applicability.

Appendix A Implementation for Solving the HJI Equation

Solving the Hamilton-Jacobi-Isaacs problem (14) for the value function $V$ is not feasible analytically, but it can be done numerically. Following [13], we employ a deep neural network (DNN) to compute a numerical solution for $V$ , with the sine function as the activation function. The sine activation function provides smooth gradients of the output of the neural network [21], which is necessary for computing the controller and disturbance (15). The input to the DNN is the full state $(x,z,t)$ , and the output is trained to be the scalar value function $V$ .

The training loss $\ell$ is defined as $\ell(\theta)=\ell_{B}(\theta)+\lambda\ell_{D}(\theta)$ , where


$\displaystyle\ell_{B}(\theta)$	$\displaystyle\coloneqq\frac{1}{\|S_{{t_{f}}}\|}\sum_{(x,z,{t_{f}})\in S_{{t_{f}}}}{\lVert V_{\theta}(x,z,{t_{f}})-\frac{1}{{t_{f}}-{t_{0}}}\sum_{k\in\mathcal{K}}{\left(\Lambda_{k}z_{k}^{2}\right)}-h(x)\rVert},$	(20a)
$\displaystyle\ell_{D}(\theta)$	$\displaystyle\coloneqq\frac{1}{\|S\|}\sum_{(x,z,t)\in S}{\Biggr{\lVert}\frac{\partial V_{\theta}}{\partial t}+\min_{u\in\mathcal{U}}\max_{d\in\mathcal{D}}\biggr{\{}\hat{H}(x,z,t,u,d)\biggl{\}}\Biggl{\rVert}},$	(20b)

where $V_{\theta}$ is the output of the DNN with parameters $\theta$ , $S$ is a sampled set of state data points of the form $(x,z,t)$ , and $S_{t_{f}}\subseteq S$ is the subset of $S$ of all state data points with $t={t_{f}}$ , and $\hat{H}$ is defined in (14b).

In the loss function, $\ell_{B}$ is the loss associated with the terminal boundary condition (14c), so it is only considered for state data points where $t={t_{f}}$ . By contrast, $\ell_{D}$ is the loss associated with the HJI differential equation (14a), so it is considered for all state points sampled. The tuning parameter $\lambda$ is used to adjust the weighting of the terminal and differential loss components.

The training of the DNN was performed as in [13], such that the terminal condition (14c) was learned first, and then the solution was learned expanding backwards from $t={t_{f}}$ to $t={t_{0}}$ , following the differential equation (14a). The DNN was trained for a total of $1.5\times 10^{6}$ iterations, which took an average of 247 minutes on an NVIDIA GeForce RTX 3080 GPU.

Appendix B Proofs

B.1 Proof of Lemma 4.1

Proof B.1.

Integrating the differential equation for $z_{k}$ (12), recognizing that $\phi_{k}$ is constant, produces

z_{k}(x(\cdot),{t_{f}})=\int_{{t_{0}}}^{{t_{f}}}{F_{k}(M\circ x(t))dt}-({t_{f}}-{t_{0}})\phi_{k}.

(21)

The definition of $c_{k}$ (5) can be used to simplify to

z_{k}(x,{t_{f}})=({t_{f}}-{t_{0}})\left(c_{k}(x(\cdot),{t_{f}})-\phi_{k}\right).

(22)

Finally, substituting this into (11) yields

\frac{1}{{t_{f}}-{t_{0}}}\sum_{k\in\mathcal{K}}{\left(\Lambda_{k}z_{k}^{2}\right)}=\sum_{k\in\mathcal{K}}{\Lambda_{k}\left(c_{k}(x(\cdot),{t_{f}})-\phi_{k}\right)^{2}}.

(23)

This completes the proof. ∎

B.2 Proof of Theorem 4.3

Proof B.2.

For an alternative proof, consult [16]. In (7), repeated here, we define the value function as the optimal cost of the remainder of the trajectory, given by

\displaystyle V(x,z,{t_{1}})

\displaystyle=\min_{u(\cdot)}\max_{d(\cdot)}\left\{J(x,{t_{1}},u(\cdot),d(\cdot)\right)\}.

(24)

Using the ergodic cost function (13b), this becomes

	$\displaystyle V(x,z,{t_{1}})=$	$\displaystyle\min_{u(\cdot)}\max_{d(\cdot)}\bigg{\{}\int_{{t_{1}}}^{{t_{f}}}{g(x(\tau),u(\tau))d\tau}$
		$\displaystyle\phantom{\min_{u(\cdot)}\max_{d(\cdot)}\bigg{\{}}+\frac{1}{{t_{f}}-{t_{0}}}\sum_{k\in\mathcal{K}}{\left(\Lambda_{k}z_{k}({t_{f}})^{2}\right)}+h(x({t_{f}}))\bigg{\}}.$		(25)

Note that in this and following equations, $x(t),z(t),u(t)$ are governed by the constraints (13c)-(13e). Substituting ${t_{1}}={t_{f}}$ into the definition of $V$ (25) directly yields the HJI boundary condition

V(x,z,{t_{f}})=\frac{1}{{t_{f}}-{t_{0}}}\sum_{k\in\mathcal{K}}{\left(\Lambda_{k}z_{k}({t_{f}})^{2}\right)}+h(x({t_{f}})).

(26)

Next, we define $t\in({t_{1}},{t_{f}})$ , then split the integral in (25), yielding

$\displaystyle V(x,z,{t_{1}})=$	$\displaystyle\min_{u(\cdot)}\max_{d(\cdot)}\bigg{\{}\int_{{t_{1}}}^{t}{g(x(\tau),u(\tau))d\tau}+\int_{t}^{{t_{f}}}{g(x(\tau),u(\tau))d\tau}$
	$\displaystyle\phantom{\min_{u(\cdot)}\max_{d(\cdot)}\bigg{\{}}+\frac{1}{{t_{f}}-{t_{0}}}\sum_{k\in\mathcal{K}}{\left(\Lambda_{k}z_{k}({t_{f}})^{2}\right)}+h(x({t_{f}}))\bigg{\}}$	(27)
$\displaystyle=$	$\displaystyle\min_{u(\cdot)}\max_{d(\cdot)}\bigg{\{}\int_{{t_{1}}}^{t}{g(x(\tau),u(\tau))d\tau}+V(x(t),z(t),t)\bigg{\}}$	(28)

Equation (28) is the Isaacs recursive principle of optimality (c.f. Eq. 2.7 of [31]). Taking the derivative with respect to $t$ yields

	$\displaystyle 0=$	$\displaystyle\min_{u(t)}\max_{d(t)}\left\{g(x(t),u(t))+\frac{\partial V}{\partial x}\frac{dx(t)}{dt}+\frac{\partial V}{\partial z}\frac{dz(t)}{dt}+\frac{\partial V}{\partial t}\right\}$		(29)
	$\displaystyle-\frac{\partial V}{\partial t}=$	$\displaystyle\min_{u(t)}\max_{d(t)}\left\{g(x(t),u(t))+\frac{\partial V}{\partial x}\frac{dx(t)}{dt}+\frac{\partial V}{\partial z}\frac{dz(t)}{dt}\right\}$		(30)

Since all terms are evaluated at time $t$ , we can replace $x(t)$ with $x$ , and likewise for $z,u,d$ . We can now substitute the equations for $\frac{dx}{dt}$ (1) and $\frac{dz}{dt}$ (12), yielding the HJI differential condition

	$\displaystyle-\frac{\partial V(x,z,t)}{\partial t}$	$\displaystyle=\min_{u\in\mathcal{U}}\max_{d\in\mathcal{D}}\bigg{\{}g(x,u)+\frac{\partial V(x,z,t)}{\partial x}f(x,u,d)$
		$\displaystyle\phantom{\min_{u\in\mathcal{U}}\max_{d\in\mathcal{D}}\bigg{\{}abc}+\sum_{k\in\mathcal{K}}{\left(\frac{\partial V(x,z,t)}{\partial z_{k}}\left(F_{k}(M\circ x)-\phi_{k}\right)\right)}\bigg{\}}.$		(31)

By inspection, it is clear that for any $(x,z,t)$ , the optimal control $u$ and disturbance $d$ are precisely those that satisfy the minimax condition in (31). ∎

Appendix C Derivation of the Ergodic HJI Equation

This section presents an alternative derivation of the HJI equation for ergodic exploration, presented above as (14). Notably, although we begin from the canonical formulation of the ergodic metric, which does not include extended state variables, the resulting HJI equation will nevertheless contain extended state variables.

C.1 Problem Setup

We start with the simple cost function

$\displaystyle J(x(\cdot),u(\cdot))\coloneqq$	$\displaystyle\smallint_{{t_{0}}}^{{t_{f}}}{g(x(t),u(t))dt}+h(x({t_{f}}))+\mathcal{E}(x(\cdot),{t_{f}})$	(32)
$\displaystyle=$	$\displaystyle\smallint_{{t_{0}}}^{{t_{f}}}{g(x(t),u(t))dt}+h(x({t_{f}}))$
	$\displaystyle+\sum_{k}{\Lambda_{k}\left(\frac{1}{{t_{f}}-{t_{0}}}\int_{t_{0}}^{t_{f}}F_{k}(M\circ x(t))dt-\phi_{k}\right)^{2}}$	(33)

First, we bring the $\phi_{k}$ into the integral of the ergodic metric, yielding

	$\displaystyle J(x(\cdot),u(\cdot))=$	$\displaystyle\smallint_{{t_{0}}}^{{t_{f}}}{g(x(t),u(t))dt}+h(x({t_{f}}))$
		$\displaystyle+\sum_{k}{\Lambda_{k}\left(\frac{1}{{t_{f}}-{t_{0}}}\int_{t_{0}}^{t_{f}}\left(F_{k}(M\circ x(t))-\phi_{k}\right)dt\right)^{2}}$		(34)

For notational ease, we define $\epsilon_{k}(t)\coloneqq F_{k}(M\circ x(t))-\phi_{k}$ . The ergodic metric is the sum of a single squared term for each value of $k$ , so it can be rewritten as

\mathcal{E}(x,{t_{f}})=\sum_{k}\left(\frac{\Lambda_{k}}{({t_{f}}-{t_{0}})^{2}}\left(\smallint_{{t_{0}}}^{{t_{f}}}{\epsilon_{k}(t)dt}\right)^{2}\right),

(35)

C.2 Solving for the HJI Equation

Define times ${t_{1}},{t_{2}}$ such that ${t_{0}}<{t_{1}}<{t_{2}}<{t_{f}}$ . We expand the expression for the cost function (34) by substituting in the definition of the ergodic metric (35) and recognizing that $\smallint_{{t_{0}}}^{{t_{f}}}{\epsilon_{k}(t)dt}=\smallint_{{t_{0}}}^{{t_{1}}}{\epsilon_{k}(t)dt}+\smallint_{{t_{1}}}^{{t_{2}}}{\epsilon_{k}(t)dt}+\smallint_{{t_{2}}}^{{t_{f}}}{\epsilon_{k}(t)dt}$ , yielding

	$\displaystyle J(x(\cdot),u(\cdot))=$	$\displaystyle\smallint_{{t_{0}}}^{{t_{f}}}{g(x(t),u(t))dt}+h(x({t_{f}}))$
		$\displaystyle+\sum_{k\in\mathcal{K}}\Bigg{(}\frac{\Lambda_{k}}{({t_{f}}-{t_{0}})^{2}}\bigg{(}\smallint_{{t_{0}}}^{{t_{1}}}{\epsilon_{k}(t)dt}+\smallint_{{t_{1}}}^{{t_{2}}}{\epsilon_{k}(t)dt}+\smallint_{{t_{2}}}^{{t_{f}}}{\epsilon_{k}(t)dt}\bigg{)}^{2}\Bigg{)}.$		(36)

Expanding yields

$\displaystyle J(x(\cdot),u(\cdot))=$	$\displaystyle\smallint_{{t_{0}}}^{{t_{1}}}{g(x(t),u(t))dt}+\smallint_{{t_{1}}}^{{t_{2}}}{g(x(t),u(t))dt}+\smallint_{{t_{2}}}^{{t_{f}}}{g(x(t),u(t))dt}+h(x({t_{f}}))$
	$\displaystyle+\sum_{k\in\mathcal{K}}\Bigg{(}\frac{\Lambda_{k}}{({t_{f}}-{t_{0}})^{2}}\Bigg{[}\left(\smallint_{{t_{0}}}^{{t_{1}}}{\epsilon_{k}(t)dt}\right)^{2}+\left(\smallint_{{t_{1}}}^{{t_{2}}}{\epsilon_{k}(t)dt}\right)^{2}$
	$\displaystyle\hskip 72.26999pt+2\left(\smallint_{{t_{0}}}^{{t_{1}}}{\epsilon_{k}(t)dt}\right)\left(\smallint_{{t_{1}}}^{{t_{2}}}{\epsilon_{k}(t)dt}\right)$
	$\displaystyle\hskip 72.26999pt+2\left[\smallint_{{t_{0}}}^{{t_{1}}}{\epsilon_{k}(t)dt}+\smallint_{{t_{1}}}^{{t_{2}}}{\epsilon_{k}(t)dt}\right]\left(\smallint_{{t_{2}}}^{{t_{f}}}{\epsilon_{k}(t)dt}\right)$
	$\displaystyle\hskip 72.26999pt+\left(\smallint_{{t_{2}}}^{{t_{f}}}{\epsilon_{k}(t)dt}\right)^{2}\Bigg{]}\Bigg{)}\Bigg{\}}.$	(37)

We define the cost-to-go $V(x_{[{t_{0}},\tau]},\tau)$ as the optimal value of the remaining running cost and the ergodic metric of entire trajectory, where the portion of the trajectory prior to $t=\tau$ is held fixed. We use the subscript ${\boxed{\cdot}}_{[t_{a},t_{b}]}$ to denote the portion of the trajectory, control, or disturbance function limited to the time range $t\in[t_{a},t_{b}]$ . Applying this definition of $V$ to (37) produces

$\displaystyle V(x_{[{t_{0}},{t_{2}}]},{t_{2}})=$	$\displaystyle\left(\smallint_{{t_{0}}}^{{t_{2}}}{\epsilon_{k}(t)dt}\right)^{2}$
	$\displaystyle+\min_{u_{[{t_{2}},{t_{f}}]}}{{\max_{d_{[{t_{2}},{t_{f}}]}}}}\Bigg{\{}\smallint_{{t_{2}}}^{{t_{f}}}{g(x(t),u(t))dt}+h(x({t_{f}}))$
	$\displaystyle\hskip 36.135pt+\sum_{k\in\mathcal{K}}\Bigg{(}\frac{\Lambda_{k}}{({t_{f}}-{t_{0}})^{2}}\Bigg{[}2\left(\smallint_{{t_{0}}}^{{t_{2}}}{\epsilon_{k}(t)dt}\right)\left(\smallint_{{t_{2}}}^{{t_{f}}}{\epsilon_{k}(t)dt}\right)$
	$\displaystyle\hskip 108.405pt+\left(\smallint_{{t_{2}}}^{{t_{f}}}{\epsilon_{k}(t)dt}\right)^{2}\Bigg{]}\Bigg{)}\Bigg{\}}$	(38)

and

$\displaystyle V(x_{[{t_{0}},{t_{1}}]},{t_{1}})=$	$\displaystyle\left(\smallint_{{t_{0}}}^{{t_{1}}}{\epsilon_{k}(t)dt}\right)^{2}$
	$\displaystyle+\min_{u_{[{t_{1}},{t_{f}}]}}{{\max_{d_{[{t_{1}},{t_{f}}]}}}}\Bigg{\{}\smallint_{{t_{1}}}^{{t_{2}}}{g(x(t),u(t))dt}+\smallint_{{t_{2}}}^{{t_{f}}}{g(x(t),u(t))dt}+h(x({t_{f}}))$
	$\displaystyle\hskip 36.135pt+\sum_{k\in\mathcal{K}}\Bigg{(}\frac{\Lambda_{k}}{({t_{f}}-{t_{0}})^{2}}\Bigg{[}\left(\smallint_{{t_{1}}}^{{t_{2}}}{\epsilon_{k}(t)dt}\right)^{2}$
	$\displaystyle\hskip 90.3375pt+2\left(\smallint_{{t_{0}}}^{{t_{1}}}{\epsilon_{k}(t)dt}\right)\left(\smallint_{{t_{1}}}^{{t_{2}}}{\epsilon_{k}(t)dt}\right)$
	$\displaystyle\hskip 90.3375pt+2\left[\smallint_{{t_{0}}}^{{t_{1}}}{\epsilon_{k}(t)dt}+\smallint_{{t_{1}}}^{{t_{2}}}{\epsilon_{k}(t)dt}\right]\left(\smallint_{{t_{2}}}^{{t_{f}}}{\epsilon_{k}(t)dt}\right)$
	$\displaystyle\hskip 90.3375pt+\left(\smallint_{{t_{2}}}^{{t_{f}}}{\epsilon_{k}(t)dt}\right)^{2}\Bigg{]}\Bigg{)}\Bigg{\}}$	(39)
$\displaystyle=$	$\displaystyle\min_{u_{[{t_{1}},{t_{2}}]}}{{\max_{d_{[{t_{1}},{t_{2}}]}}}}\Bigg{\{}\smallint_{{t_{1}}}^{{t_{2}}}{g(x(t),u(t))dt}+\left(\smallint_{{t_{0}}}^{{t_{2}}}{\epsilon_{k}(t)dt}\right)^{2}$
	$\displaystyle\hskip 18.06749pt+\min_{u_{[{t_{2}},{t_{f}}]}}{{\max_{d_{[{t_{2}},{t_{f}}]}}}}\Bigg{\{}\smallint_{{t_{2}}}^{{t_{f}}}{g(x(t),u(t))dt}+h(x({t_{f}}))$
	$\displaystyle\hskip 36.135pt+\sum_{k\in\mathcal{K}}\Bigg{(}\frac{\Lambda_{k}}{({t_{f}}-{t_{0}})^{2}}\Bigg{[}2\left(\smallint_{{t_{0}}}^{{t_{2}}}{\epsilon_{k}(t)dt}\right)\left(\smallint_{{t_{2}}}^{{t_{f}}}{\epsilon_{k}(t)dt}\right)$
	$\displaystyle\hskip 108.405pt+\left(\smallint_{{t_{2}}}^{{t_{f}}}{\epsilon_{k}(t)dt}\right)^{2}\Bigg{]}\Bigg{)}\Bigg{\}}\Bigg{\}}$	(40)
$\displaystyle=$	$\displaystyle\min_{u_{[{t_{1}},{t_{2}}]}}{{\max_{d_{[{t_{1}},{t_{2}}]}}}}\Big{\{}\smallint_{{t_{1}}}^{{t_{2}}}{g(x(t),u(t))dt}+V(x_{[{t_{0}},{t_{2}}]},{t_{2}})\Big{\}}.$	(41)

The final result (41) is the recursive Isaacs principle of optimality for a problem with running cost. Note that aside from the different notations for the arguments of $V$ , this is equivalent to (28), from the proof of Theorem 4.3 above.

As written in (C.2), the cost-to-go $V$ at time $\tau$ depends on the entire trajectory up to that time, $x_{[{t_{0}},\tau]}$ . However, inspection shows that $x_{[{t_{0}},\tau]}$ is only used in two ways:

•

The entire history $x_{[{t_{0}},\tau]}$ is used in the calculation of $\smallint_{{t_{0}}}^{\tau}{\epsilon_{k}(t)dt},\forall k\in\mathcal{K}$ .
•

The state $x(\tau)$ is necessary for rolling out the dynamics at times greater than $\tau$ .

Thus, $V$ at time $\tau$ can be written as $V(x(\tau),\left\{\smallint_{{t_{0}}}^{\tau}{\epsilon_{k}(t)dt},\forall k\in\mathcal{K}\right\},\tau)$ . Taking the limit of (41) as $\Delta t\coloneqq{t_{2}}-{t_{1}}\to 0$ results in

$\displaystyle V$	$\displaystyle\left(x({t_{1}}),\left\{\smallint_{{t_{0}}}^{{t_{1}}}{\epsilon_{k}(t)dt},\forall k\in\mathcal{K}\right\},{t_{1}}\right)$	(42)
	$\displaystyle=\min_{u\in\mathcal{U}}\max_{d\in\mathcal{D}}\Bigg{\{}g(x({t_{1}}),u({t_{1}}))\Delta t$
	$\displaystyle\hskip 18.06749pt+V\left(x({t_{1}})+f(x,u,d)\Delta t,\left\{\smallint_{{t_{0}}}^{{t_{1}}}{\epsilon_{k}(t)dt}+\epsilon_{k}({t_{1}})\Delta t,\forall k\in\mathcal{K}\right\},{t_{1}}+\Delta t\right)\Bigg{\}}.$

We use the Taylor expansion of the value function, resulting in

$\displaystyle V$	$\displaystyle\left(x({t_{1}}),\left\{\smallint_{{t_{0}}}^{{t_{1}}}{\epsilon_{k}(t)dt},\forall k\in\mathcal{K}\right\},{t_{1}}\right)$	(43)
	$\displaystyle=\min_{u\in\mathcal{U}}\max_{d\in\mathcal{D}}\Bigg{\{}g(x({t_{1}}),u({t_{1}}))\Delta t+V\left(x({t_{1}}),\left\{\smallint_{{t_{0}}}^{{t_{1}}}{\epsilon_{k}(t)dt},\forall k\in\mathcal{K}\right\},{t_{1}}\right)$
	$\displaystyle\hskip 18.06749pt+\frac{\partial V}{\partial x}f(x,u,d)\Delta t+\sum_{k\in\mathcal{K}}{\left(\frac{\partial V}{\partial\left(\smallint_{{t_{0}}}^{t}{\epsilon_{k}(t)dt}\right)}\epsilon_{k}\right)}\Delta t+\frac{\partial V}{\partial t}\Delta t+o(\Delta t^{2})\Bigg{\}},$

where $o(\Delta t^{2})$ is a term that approaches 0 much faster than $\Delta t$ as $\Delta t\to 0$ . Dividing by $\Delta t$ and then taking the limit as $\Delta t\to 0$ yields

0=\min_{u\in\mathcal{U}}\max_{d\in\mathcal{D}}\Bigg{\{}g(x,u)+\frac{\partial V}{\partial x}f(x,u,d)+\sum_{k\in\mathcal{K}}{\left(\frac{\partial V}{\partial\left(\smallint_{{t_{0}}}^{t}{\epsilon_{k}(t)dt}\right)}\epsilon_{k}\right)}+\frac{\partial V}{\partial t}\Bigg{\}}.

(44)

We recognize now that $\smallint_{{t_{0}}}^{\tau}{\epsilon_{k}(t)dt}$ is precisely equivalent to $z_{k}(\tau)$ . We refrained from making this substitution earlier in order to emphasize that this derivation does not depend on any properties of $z_{k}$ proved elsewhere in the paper. Making this substitution, expanding $\epsilon_{k}$ , and rearranging produces the differential condition

	$\displaystyle-\frac{\partial V}{\partial t}=\min_{u\in\mathcal{U}}\max_{d\in\mathcal{D}}\Bigg{\{}$	$\displaystyle g(x,u)+\frac{\partial V}{\partial x}f(x,u,d)$
		$\displaystyle+\sum_{k\in\mathcal{K}}{\left(\frac{\partial V}{\partial z_{k}}\left(F_{k}(M(x))-\phi_{k}\right)\right)}\Bigg{\}}.$		(45)

Finally, evaluating the non-recursive formulation of the value function (C.2) at $t={t_{f}}$ yields the terminal condition

V(x,z,{t_{f}})=\sum_{k\in\mathcal{K}}{\left(\frac{\Lambda_{k}}{({t_{f}}-{t_{0}})^{2}}z_{k}^{2}\right)}+h(x),\ \forall x\in\mathcal{X},\forall z\in\mathcal{Z}

(46)

Taken together, (45) and (46) constitute a PDE that defines the value function, and they are identical to the HJI PDE (14) derived from the extended-state formulation of the ergodic control problem (Problem 3).

References

[1] T. Tomic, K. Schmid, P. Lutz, A. Domel, M. Kassecker, E. Mair, I. L. Grixa, F. Ruess, M. Suppa, and D. Burschka, “Toward a fully autonomous UAV: Research platform for indoor and outdoor urban search and rescue,” IEEE Robotics & Automation Magazine, vol. 19, no. 3, pp. 46–56, 2012.
[2] D. C. Schedl, I. Kurmi, and O. Bimber, “An autonomous drone for search and rescue in forests using airborne optical sectioning,” Science Robotics, vol. 6, no. 55, p. eabg1188, 2021.
[3] S. Nuske, S. Choudhury, S. Jain, A. Chambers, L. Yoder, S. Scherer, L. Chamberlain, H. Cover, and S. Singh, “Autonomous exploration and motion planning for an unmanned aerial vehicle navigating rivers,” Journal of Field Robotics, vol. 32, no. 8, pp. 1141–1162, 2015.
[4] H. Qin, Z. Meng, W. Meng, X. Chen, H. Sun, F. Lin, and M. H. Ang, “Autonomous exploration and mapping system using heterogeneous UAVs and UGVs in GPS-denied environments,” IEEE Transactions on Vehicular Technology, vol. 68, no. 2, pp. 1339–1350, 2019.
[5] J. G. Serna, F. Vanegas, F. Gonzalez, and D. Flannery, “A review of current approaches for UAV autonomous mission planning for Mars biosignatures detection,” in 2020 IEEE Aerospace Conference, pp. 1–15, 2020.
[6] L. Dressel and M. J. Kochenderfer, “On the optimality of ergodic trajectories for information gathering tasks,” in 2018 Annual American Control Conference (ACC), pp. 1855–1861, IEEE, 2018.
[7] G. Mathew and I. Mezić, “Metrics for ergodicity and design of ergodic dynamics for multi-agent systems,” Physica D: Nonlinear Phenomena, vol. 240, no. 4, pp. 432–442, 2011.
[8] J. Lygeros, “On reachability and minimum cost optimal control,” Automatica, vol. 40, no. 6, pp. 917–927, 2004.
[9] S. Bansal, M. Chen, S. Herbert, and C. J. Tomlin, “Hamilton-Jacobi reachability: A brief overview and recent advances,” in 2017 IEEE 56th Annual Conference on Decision and Control (CDC), pp. 2242–2253, IEEE, 2017.
[10] M. Chen, S. L. Herbert, M. S. Vashishtha, S. Bansal, and C. J. Tomlin, “Decomposition of reachable sets and tubes for a class of nonlinear systems,” IEEE Transactions on Automatic Control, vol. 63, no. 11, pp. 3675–3688, 2018.
[11] G. De La Torre, K. Flasskamp, A. Prabhakar, and T. D. Murphey, “Ergodic exploration with stochastic sensor dynamics,” in 2016 American Control Conference (ACC), pp. 2971–2976, IEEE, 2016.
[12] D. Dong, H. Berger, and I. Abraham, “Time optimal ergodic search,” in 2023 Robotics: Science and Systems (RSS), 2023.
[13] S. Bansal and C. J. Tomlin, “DeepReach: A deep learning approach to high-dimensional reachability,” in 2021 IEEE International Conference on Robotics and Automation (ICRA), pp. 1817–1824, IEEE, 2021.
[14] Y. S. Shao, C. Chen, S. Kousik, and R. Vasudevan, “Reachability-based trajectory safeguard (RTS): A safe and fast reinforcement learning safety layer for continuous control,” IEEE Robotics and Automation Letters, vol. 6, no. 2, pp. 3663–3670, 2021.
[15] K. D. Julian and M. J. Kochenderfer, “Guaranteeing safety for neural network-based aircraft collision avoidance systems,” in 2019 IEEE/AIAA 38th Digital Avionics Systems Conference (DASC), pp. 1–10, IEEE, 2019.
[16] L. C. Evans and P. E. Souganidis, “Differential games and representation formulas for solutions of Hamilton-Jacobi-Isaacs equations,” Indiana University Mathematics Journal, vol. 33, no. 5, pp. 773–797, 1984.
[17] F. Jiang, G. Chou, M. Chen, and C. J. Tomlin, “Using neural networks to compute approximate and guaranteed feasible Hamilton-Jacobi-Bellman PDE solutions,” 2017.
[18] S. L. Herbert, M. Chen, S. Han, S. Bansal, J. F. Fisac, and C. J. Tomlin, “FaSTrack: A modular framework for fast and guaranteed safe motion planning,” in 2017 IEEE 56th Annual Conference on Decision and Control (CDC), pp. 1517–1522, IEEE, 2017.
[19] A. Majumdar and R. Tedrake, “Funnel libraries for real-time robust feedback motion planning,” The International Journal of Robotics Research, vol. 36, no. 8, pp. 947–982, 2017.
[20] C. He, Z. Gong, M. Chen, and S. Herbert, “Efficient and guaranteed Hamilton–Jacobi reachability via self-contained subsystem decomposition and admissible control sets,” IEEE Control Systems Letters, vol. 7, pp. 3824–3829, 2023.
[21] V. Sitzmann, J. N. P. Martel, A. W. Bergman, D. B. Lindell, and G. Wetzstein, “Implicit neural representations with periodic activation functions,” 2020.
[22] T. Nakamura-Zimmerer, Q. Gong, and W. Kang, “Adaptive deep learning for high-dimensional Hamilton-Jacobi-Bellman equations,” SIAM Journal on Scientific Computing, vol. 43, no. 2, pp. A1221–A1247, 2021.
[23] S. Lin, “Computer solutions of the traveling salesman problem,” The Bell System Technical Journal, vol. 44, no. 10, pp. 2245–2269, 1965.
[24] D. L. Miller and J. F. Pekny, “Exact solution of large asymmetric traveling salesman problems,” Science, vol. 251, no. 4995, pp. 754–761, 1991.
[25] J. Cortes, S. Martinez, T. Karatas, and F. Bullo, “Coverage control for mobile sensing networks,” IEEE Transactions on Robotics and Automation, vol. 20, no. 2, pp. 243–255, 2004.
[26] J. Chen and P. Dames, “Distributed and collision-free coverage control of a team of mobile sensors using the convex uncertain Voronoi diagram,” in 2020 American Control Conference (ACC), pp. 5307–5313, 2020.
[27] H. Salman, E. Ayvali, and H. Choset, “Multi-agent ergodic coverage with obstacle avoidance,” Proceedings of the International Conference on Automated Planning and Scheduling, vol. 27, pp. 242–249, 2017.
[28] C. Lerch, D. Dong, and I. Abraham, “Safety-critical ergodic exploration in cluttered environments via control barrier functions,” in 2023 IEEE International Conference on Robotics and Automation (ICRA), pp. 10205–10211, 2023.
[29] S. E. Scott, T. C. Redd, L. Kuznetsov, I. Mezić, and C. K. Jones, “Capturing deviation from ergodicity at different scales,” Physica D: Nonlinear Phenomena, vol. 238, no. 16, pp. 1668–1679, 2009.
[30] R. Isaacs, Differential games: a mathematical theory with applications to warfare and pursuit, control and optimization. Wiley, 1965.
[31] E. Galperin, “The Isaacs equation for differential games, totally optimal fields of trajectories and related problems,” Computers & Mathematics with Applications, vol. 55, no. 6, pp. 1333–1362, 2008.
[32] L. M. Miller, Y. Silverman, M. A. MacIver, and T. D. Murphey, “Ergodic exploration of distributed information,” IEEE Transactions on Robotics, vol. 32, no. 1, pp. 36–52, 2016.
[33] J. A. Preiss*, W. Hönig*, G. S. Sukhatme, and N. Ayanian, “Crazyswarm: A large nano-quadcopter swarm,” in IEEE International Conference on Robotics and Automation (ICRA), pp. 3299–3304, IEEE, 2017. Software available at https://github.com/USC-ACTLab/crazyswarm.

RAnGE: Reachability Analysis for Guaranteed Ergodicity