Automatic differentiation of non-holonomic fast marching for computing most threatening trajectories under sensors surveillance

TL;DR

This paper develops a method combining non-holonomic fast marching and automatic differentiation to compute optimal evasive trajectories and optimize surveillance placement in a monitored region.

Contribution

It introduces a novel approach integrating non-holonomic fast marching with automatic differentiation for strategic surveillance and evasion planning.

Findings

Efficient computation of optimal trajectories considering vehicle maneuverability.

Gradient-based optimization of surveillance sensor placement.

Demonstrated effectiveness in simulated surveillance scenarios.

Abstract

We consider a two player game, where a first player has to install a surveillance system within an admissible region. The second player needs to enter the the monitored area, visit a target region, and then leave the area, while minimizing his overall probability of detection. Both players know the target region, and the second player knows the surveillance installation details.Optimal trajectories for the second player are computed using a recently developed variant of the fast marching algorithm, which takes into account curvature constraints modeling the second player vehicle maneuverability. The surveillance system optimization leverages a reverse-mode semi-automatic differentiation procedure, estimating the gradient of the value function related to the sensor location in time N log N.

Figures

Equations

$$\mathfrak{C}(\Xi ,\Gamma )$$

$${\nabla }_{\xi }\mathfrak{C}(\xi ,\Gamma ),$$

$$\gamma (0)$$

$$\mathfrak{C}(\xi ,\Gamma )$$

$${\mathcal{F}}_{\xi }(x,\dot{x})$$

$$\forall x\in \Omega \setminus {\rm Y},{\mathcal{H}}_{\xi }(x,{\nabla }_x{u}_{\xi }(x))$$

$$\eta (0)$$

$$\mathfrak{C}(\xi ,\Gamma )$$

$${\mathcal{F}}_{\xi }^{\pm }((p,n),(\dot{p},\dot{n})):=\{\begin{array}{ll}+\infty & \text{if }\dot{p}\ne \Vert \dot{p}\Vert n,\\ \Vert \dot{p}\Vert {\mathcal{C}}_{\xi }(p,\pm n,\dot{n}/\Vert \dot{p}\Vert )& \text{else},\end{array}$$

$${\mathcal{C}}_{\mathrm{RS}}(\kappa )$$

$${\mathcal{H}}_{\mathrm{RS}}$$

$$X=\Omega \cap (h\mathbb{Z}{)}^d,$$

$${\mathcal{H}}_{\xi }(x,\dot{x})\approx h^{-2}\sum _{y\in X}{c}_{\xi }^2(x,y)⟨x-y,\dot{x}{⟩}_{+}^2,$$

$$\sum _{y\in X}{c}_{\xi }^2(x,y)(U_{\xi }(x)-U_{\xi }(y){)}_{+}^2=h^2/2$$

$$\sum _{y\in X}{\omega }_{\xi }(x,y)\left(\mathrm{d}{U}_{\xi }(x)-\mathrm{d}{U}_{\xi }(y)+(U_{\xi }(x)-U_{\xi }(y))\mathrm{d}\ln c_{\xi }(x,y)\right)=0,$$

$${\mathcal{C}}_{\xi }(p,\dot{p},\ddot{p})={\mathcal{C}}_{\ast }(|\ddot{p}|),$$

$${\mathcal{C}}_{\xi }(p,\dot{p},\ddot{p})=\xi (p){\mathcal{C}}_{\ast }(|\ddot{p}|)$$

$${\mathcal{C}}_{\xi }(p,\dot{p},\ddot{p})={\mathcal{C}}_{\ast }(\kappa )\sum _{\begin{array}{c}q\in \xi \\ [p,q]\subseteq R\end{array}}\frac{1}{\Vert q-p{\Vert }^2},$$

$${\mathcal{C}}_{\xi }(p,\dot{p},\ddot{p})={\mathcal{C}}_{\ast }(|\ddot{p}|)\sqrt{\sum _{q\in \xi }\frac{⟨\dot{p},n_{\overrightarrow{pq}}{⟩}^2+{\delta }^2⟨\dot{p},n_{\overrightarrow{pq}}^{\perp }{⟩}^2}{\Vert p-q{\Vert }^4}}$$

Full text

Optimal trajectories for the second player are computed using a recently developed variant of the fast marching algorithm, which takes into account curvature constraints modeling the second player vehicle maneuverability. The surveillance system optimization leverages a reverse-mode semi-automatic differentiation procedure, estimating the gradient of the value function related to the sensor location in time ${\mathcal{O}}(N\ln N)$.

keywords:

Anisotropic fast marching, motion planning, sensor placement, optimization

1 Introduction

This paper presents a proof of concept numerical implementation of a motion planning algorithm related to a two player game. A first player selects, within a admissible class $\Xi$, an integral cost function on paths, which takes into account their position, orientation, and possibly curvature. The second player selects a path, within an admissible class $\Gamma$, with prescribed endpoints and an intermediate keypoint. The players objective is respectively to maximize and minimize the path cost

[TABLE]

where the path $\gamma$ is parametrized at unit euclidean speed, and the final time $T(\gamma)$ is free. From a game theoretic point of view, this is a non-cooperative zero-sum game, where player $\Xi$ has no information and player $\Gamma$ has full information over the opponent’s strategy.

The game (1) typically models a surveillance problem [15], and $\exp(-{\mathfrak{C}}(\Xi,\Gamma))$ is the probability for player $\Gamma$ to visit the intermediate keypoint without being detected by player $\Xi$ (by convention, players are identified with their sets of strategies). For instance player $\Gamma$ could be a robber trying to steal an object from a museum, and player $\Xi$ would be responsible for the placement of video surveillance cameras. In another instance, player $\Xi$ is responsible for the installation of radar [1] or sonar detection systems [15], and would like to prevent vehicles sent by player $\Gamma$ from spying on some objectives without being detected.

The dependence of the cost ${\mathcal{C}}_{\xi}$ w.r.t. the path tangent $\gamma^{\prime}(t)$ models the variation of a measure of how detectable the target is (radar cross section, directivity index, etc.) w.r.t. the relative positions and orientations of the target and sensor. The dependence of ${\mathcal{C}}_{\xi}$ on the path curvature $\gamma^{\prime\prime}(t)$ models the airplane maneuverability constraints, such as the need to slow down in tight turns [6], or even a hard bound on the path curvature [5].

Strode [15] have shown the interplay of motion planning and game theory in a similar setting, on a multistatic sonar network use case, but using isotropic graph-based path planning. The same year, Barbaresco [2] used fast-marching for computing threatening paths toward a single radar, but without taking into account curvature constraints and without considering a game setting.

The main contributions of this paper are as follows:

Anisotropy and curvature penalization: Strategy optimization for player $\Gamma$ is an optimal motion planning problem, with a known cost function. This is addressed by numerically solving a generalized eikonal PDE posed on a two or three dimensional domain, and which is strongly anisotropic in the presence of a curvature penalty and a detection measurement that depends on orientation. A Fast-Marching algorithm, relying on recent adaptive stencils constructions, based on tools from lattice geometry, is used for that purpose [6, 10, 11]. In contrast, the classical fast marching method [13] used in [4] is limited to cost functions ${\mathcal{C}}_{\xi}(\gamma(t))$ independent of the path orientation $\gamma^{\prime}(t)$ and curvature $\gamma^{\prime\prime}(t)$. 2. 2.

Gradient computation for sensors placement: Strategy optimization for player $\Xi$ is typically a non-convex problem, to which various strategies can be applied, yet gradient information w.r.t. the variable $\xi\in\Xi$ is usually of help. For that purpose, we implement efficient differentiation algorithms, forward and reverse, for estimating the gradient of the value function of player $\Xi$

[TABLE]

Reverse mode differentiation reduced the computation cost of $\nabla_{\xi}{\mathfrak{C}}(\xi,\Gamma)$ from ${\mathcal{O}}(N^{2})$, as used in [4], to ${\mathcal{O}}(N\ln N)$, where $N$ denotes the number of ization points of the domain. As a result, we can reproduce examples from [4] with computation times reduced by several orders of magnitude, and address complex three dimensional problems.

Due to space constraints, this paper is focused on problem modeling, ization and numerical experiments, rather than on mathematical aspects of wellposedness and convergence analysis. Free and open source codes for reproducing (part of) the presented numerical experiments are available on the first author’s webpage333github.com/Mirebeau/HamiltonianFastMarching.

2 Mathematical background of trajectory optimization

We describe in this section the PDE formalism, based on generalized eikonal equations, used to compute the value function $\min_{\gamma\in\Gamma}{\mathfrak{C}}(\xi,\gamma)$ of the second player, where $\xi$ is known and fixed. Their discretization is discussed in §3. We distinguish two cases, depending on wether the path cost function ${\mathcal{C}}_{\xi}(x,\dot{x},\ddot{x})$ appearing in (1) depends on the last entry $\ddot{x}$, i.e. on path curvature.

2.1 Curvature independent cost

Let $\Omega\subseteq{\mathbb{E}}:={\mathbb{R}}^{2}$ be a bounded domain, and let the source set $\Upsilon$ and target set $\Theta$ be subsets of $\Omega$. For each $x\in\Omega$, let $\Gamma_{x}$ denote the set of all paths $\gamma\in C^{1}([0,T],\Omega)$, where $T=T(\gamma)$ is free, such that

[TABLE]

The problem description states that the first player needs to go from $\Upsilon$ to $\Theta$ and back, hence its set of strategies is $\Gamma=\bigcup_{x\in\Theta}\Gamma_{x}\times\Gamma_{x}$. Therefore

[TABLE]

Define the $1$-homogenous metric ${\mathcal{F}}_{\xi}:\Omega\times{\mathbb{E}}\to[0,\infty]$, the Lagrangian ${\mathcal{L}}_{\xi}$ and the Hamiltonian ${\mathcal{H}}_{\xi}$ by

[TABLE]

Under mild assumptions [3], the function $u_{\xi}:\Omega\to{\mathbb{R}}$ is the unique viscosity solution to a generalized eikonal equation

[TABLE]

with outflow boundary conditions on $\partial\Omega$. The discretization of this PDE is discussed in §3. We limit in practice our attention to Isotropic costs ${\mathcal{C}}_{\xi}(x)$, and Riemannian costs ${\mathcal{C}}_{\xi}(x,\dot{x})=\sqrt{\langle\dot{x},M_{\xi}(x)\dot{x}\rangle}$ where $M_{\xi}(x)$ is symmetric positive definite, for which efficient numerical strategies have been developed [8, 10].

2.2 Curvature dependent cost

Let $\Omega\subseteq{\mathbb{R}}^{2}\times{\mathbb{S}}^{1}$ be a bounded domain, within the three dimensional space of all positions and orientations. As before, let $\Upsilon,\Theta\subseteq\Omega$. For all $x\in\Omega$ let $\Gamma^{\pm}_{x}$ be the collection of all $\gamma\in C^{2}([0,T],\Omega)$, such that $\eta:=(\gamma,\pm\gamma^{\prime})$ satisfies

[TABLE]

Here and below, the symbol “$\pm$” must be successively replaced with “$+$” and then “$-$”. Since the first player needs to go from $\Upsilon$ to $\Theta$ and back, its set of strategies is $\Gamma=\bigcup_{x\in\Theta}\Gamma_{x}^{+}\times\Gamma_{x}^{-}$. Therefore

[TABLE]

We denoted ${\mathfrak{C}}^{\pm}$ the path cost defined in terms of the local cost ${\mathcal{C}}_{\xi}(p,\pm\dot{p},\ddot{p})$. Consider the $1$-homogeneous metric ${\mathcal{F}}^{\pm}_{\xi}:{\mathrm{T}}\Omega\to[0,\infty]$, defined on the tangent bundle to $\Omega\subseteq{\mathbb{R}}^{2}\times{\mathbb{S}}^{1}$ by

[TABLE]

where $p\in{\mathbb{R}}^{2}$, $n\in{\mathbb{S}}^{1}$ is a unit vector, and the tangent vector satisfies $\dot{p}\in{\mathbb{R}}^{2}$, $\dot{n}\perp n$.

Introducing the Lagrangian ${\mathcal{L}}^{\pm}_{\xi}=\frac{1}{2}({\mathcal{F}}^{\pm}_{\xi})^{2}$ on ${\mathrm{T}}\Omega$, and its Legendre-Fenchel dual the Hamiltonian ${\mathcal{H}}^{\pm}_{\xi}$, one can again under mild assumptions characterize $u_{\xi}^{\pm}$ as the unique viscosity solution to the generalized eikonal PDE ${\mathcal{H}}^{\pm}_{\xi}(x,\nabla u^{\pm}_{\xi}(x))=1/2$ with appropriate boundary conditions [3]. In practice, we choose cost functions of the form ${\mathcal{C}}_{\xi}(p,\dot{p},\ddot{p})={\mathcal{C}}^{\circ}_{\xi}(p,\dot{p}){\mathcal{C}}_{*}(|\ddot{p}|)$, where ${\mathcal{C}}_{*}$ is the Reeds-Shepp car or Dubins car [5] curvature penalty, namely

[TABLE]

where $\rho>0$ is a parameter which has the dimension of a curvature radius. The Dubins car can only follows path which curvature radius is $\leq\rho$, whereas the Reeds-Shape car, in the sense of [6] and without reverse gear), can rotate into place if needed. The Hamiltonian then has the explicit expression ${\mathcal{H}}((p,n),(\hat{p},\hat{n}))=\frac{1}{2}{\mathcal{C}}^{0}_{\xi}(p,n)^{-2}\allowbreak{\mathcal{H}}_{*}(n,(\hat{p},\hat{n}))$ where

[TABLE]

3 Discretization of generalized eikonal equations

We construct a discrete domain $X$ by intersecting the computational domain with an orthogonal grid of scale $h>0$

[TABLE]

where $d=2$ for the curvature independent models, $d=3$ for the other models which are posed on ${\mathbb{R}}^{2}\times({\mathbb{R}}/2\pi{\mathbb{Z}})$ —using the angular parametrization ${\mathbb{S}}^{1}\cong{\mathbb{R}}/2\pi{\mathbb{Z}}$ (in the latter periodic case, and $2\pi/h$ must be an integer). We design weights $c_{\xi}(x,y)$, $x,y\in X$ such that for any tangent vector $\dot{x}$ at $x$ one has

[TABLE]

where $a_{+}:=\max\{0,a\}$. (Expression (5), is typical although some models require a slight generalization.) The weights $c_{\xi}(x,y)$ are non-zero for only few $(x,y)\in X$ at distance $\|x-y\|={\mathcal{O}}(h)$, see Figure 1. Their construction exploits the additive structure of the discretization grid $X$ and relies on techniques from lattice geometry [14], see [6, 10, 11] for details. The generalized eikonal PDE ${\mathcal{H}}_{\xi}(x,\nabla_{x}u_{\xi}(x))=1/2$ is discretized as

[TABLE]

with adequate boundary conditions. The solution $U_{\varepsilon}:X\to{\mathbb{R}}$ to this system of equations is computed in a single pass with ${\mathcal{O}}(N\ln N)$ complexity using the Fast-Marching algorithm [13].

To be able to use the gradient to solve the problem (1), we need to differentiate the cost ${\mathfrak{C}}(\xi,\Gamma)$ w.r.t. the first player strategy $\xi\in\Xi$. In view of (3) and (4), this only requires the sensitivity of the discrete solution values $U_{\varepsilon}(x_{*})$ at the few points points $x_{*}\in X\cap\Theta$, w.r.t to variations in the weights $c_{\varepsilon}(x,y)$, $x,y\in X$. For that purpose we differentiate (6) w.r.t. $\xi$ and obtain

[TABLE]

where $\omega_{\xi}(x,y):=c^{2}_{\xi}(x,y)(U_{\xi}(x)-U_{\xi}(y))_{+}$. This equation ties ${\mathrm{d}}U_{\varepsilon}(x)$ to those ${\mathrm{d}}U_{\xi}(y)$, $y\in X$, such that $\omega_{\xi}(x,y)>0$, hence satisfying $U_{\xi}(x)<U_{\xi}(y)$ by construction. In the spirit of automatic differentiation by reverse accumulation [16], these relations are back-propagated from any $\Theta$ keypoint $x_{*}$ to the $\Upsilon$ seeds, at a modest cost ${\mathcal{O}}(N)$.

4 Numerical results

The chosen physical domain $R$ is the rectangle $[0,2]\times[0,1]$ minus some obstacles, as illustrated on Figure 2. Source point is $(0.2,0.5)$ and target keypoint $(1.8,0.5)$. The computational domain is thus $\Omega=R$ for curvature independent models and $\Omega=R\times{\mathbb{S}}^{1}$ for curvature dependent models, which is discretized on a $180\times 89$ or $180\times 89\times 60$ grid.

No intervention from the first player.

The cost function is

[TABLE]

where ${\mathcal{C}}_{*}(\kappa)$ is respectively $1$, $\sqrt{1+\rho^{2}\kappa^{2}}$ and $1+\infty\times\chi_{\rho\kappa>1}$, with $\rho:=0.3$. The differences between the three models are apparent: the curvature independent model uses the same path forward and back; the Reeds-Shepp car spreads some curvature along the way but still makes an angle at the target point; the Dubins car maintains the radius of curvature below the bound $\rho$, and its trajectory is a succession of straight and circular segments.

Fresh paint based detection.

In this toy model, see Figure 3, the first player spreads some fresh paint over the domain, and the second player is regarded as detected if he comes back covered in it from his visit to the keypoint. The cost function is

[TABLE]

where $\xi:R\to{\mathbb{R}}_{+}$ is the fresh paint density, decided by the first player, and ${\mathcal{C}}_{*}(\kappa)$ is as above. For wellposedness, we impose a bounds on the density, namely $0.1\leq\xi\leq 1$, and subtract the paint supply cost $\int_{R}\xi(p){\mathrm{d}}p$ to (1). The main interest of this specific game, also considered in [4], is that ${\mathfrak{C}}(\xi,\Gamma)$ is concave w.r.t. $\xi\in\Xi$. The optimal strategy for player $\Xi$ is: in the curvature independent case to make some “fences” of paint between close obstacles, and in the curvature penalized models to deposit paint at the edges of obstacles, as well as along specific circular arcs for the Dubins model. At the optimum for player $\Xi$, a continuum of paths have equal cost for player $\Gamma$, covering a large portion of the domain with a uniform density as illustrated on Figure 3 (bottom), see [4] for details.

Visual detection.

The first player places some cameras, which efficiency at detecting the second player decreases with distance and is blocked by obstacles, see Figure 4. The cost function is

[TABLE]

where $\xi\in\Xi$ is a subset of $R$ with prescribed cardinality, two in out experiments. The green arrows on Fig 4 originate from the current (non optimal) camera position, and point in the direction of greatest growth $\nabla{\mathfrak{C}}(\xi,\Gamma)$ for the first player objective function.

Radar based detection.

The first player places some radars on the domain $R=[0,2]\times[0,1]$, here devoid of obstacles, and the second player has to fly by undetected. The cost function is

[TABLE]

where $n_{\vec{pq}}:=(q-p)/\|q-p\|$, and $\xi$ is a three element subset of $[0.4,1.6]$. The parameter $\delta$ is set to $1$ for an isotropic radar cross section (RCS), or to $0.2$ for an anisotropic RCS. In the latter where a plane showing its side to radar is five times less likely to be detected than a plane showing its nose or back, at the same position. Green arrows on Figure 5 point from the original position to the (locally) optimized position for player $\Xi$. At the equilibrium, several paths are optimal for player $\Gamma$, shown in red on Fig 5.

Computational cost

On a standard Laptop computer (2.7Ghz, 16GB ram), optimizing the second player objective, by solving a generalized eikonal equation, takes $\approx 1$s in the curvature dependent case, and $\approx 60$ times less in the curvature independent case thanks to the absence of angular discretization of the domain. Optimizing the first player objective takes $\approx 100$ $L-BFGS$ iterations, each one taking at most $8$s. For the stability of the minimization procedure, the problems considered were slightly regularized by the use of soft-minimum functions and by “blurring” the target keypoint over the $3\times 3$ box of adjacent pixels.

5 Conclusion

We have modeled a motion planning problem that take into account navigation constraints cand minimizing an anisotropic probability of detection, while computing the gradient of the value function related to the sensors location. This model is thus useful for surveillance applications modeled as a two-player zero-sum game involving a target that tries to avoid detection.