# Hybrid control for optimal visiting problems for a single player and for a
crowd

Fabio Bagagiolo  Dipartimento di Matematica
Università di Trento
Via Sommarive, 14, 38123 Povo (TN) Italy<EMAIL_ADDRESS>Adriano
Festa  Dipartimento di Scienze Matematiche “Giuseppe Luigi Lagrange”
Politecnico di Torino
Corso Duca degli Abruzzi, 24, 10129 Torino Italy<EMAIL_ADDRESS>Luciano Marzufero  Dipartimento di Matematica
Università di Trento
Via Sommarive, 14, 38123 Povo (TN) Italy<EMAIL_ADDRESS>
###### Abstract.

In an optimal visiting problem, we want to control a trajectory that has to
pass as close as possible to a collection of target points or regions. We
introduce a hybrid control-based approach for the classic problem where the
trajectory can switch between a group of discrete states related to the
targets of the problem. The model is subsequently adapted to a mean-field game
framework to study viability and crowd fluxes to model a multitude of
indistinguishable players.

###### Key words and phrases:

Hybrid systems, optimal control, Hamilton-Jacobi equations, mean-field games.

###### 1991 Mathematics Subject Classification:

Primary 49L25; Secondary 35Q82; 35Q84

##  1\. Introduction

In this paper, we deal with the problem of optimizing a trajectory to “visit”,
i.e., to touch or at least to pass as close as possible, to a collection of
targets. In the following, we refer to this problem as _optimal visiting_. The
issue presents various inherent difficulties: some related to its high
computational complexity (shared with other well-known optimization problems
as the “Travelling salesman problem” [27]) and other related to its possible
continuous/discontinuous nature.

Let us state the problem more precisely: consider the controlled dynamics

$\begin{cases}y^{\prime}(s)=f(y(s),\alpha(s),s),&s\in]t,T]\\\
y(t)=x,&x\in{\mathbb{R}}^{d}\end{cases},$ (1.1)

where $\alpha:[t,T]\longrightarrow A$ is a measurable control function, and
the dynamics $f:{\mathbb{R}}^{d}\times
A\times[0,+\infty[\longrightarrow{\mathbb{R}}^{d}$ is suitably regular.
Consider a collection of $N$ closed disjoint target sets in
${\mathbb{R}}^{d}$,
$\\{\operatorname{\mathcal{T}}_{1},\operatorname{\mathcal{T}}_{2},\ldots,\operatorname{\mathcal{T}}_{N}\\}$,
and $y_{(x,t)}(\cdot;{\alpha})$ a solution of (1.1) related to a starting
point $x$, a starting time $t$ and a control $\alpha$. We can write the
optimal visiting problem just considering, for example, the minimization of a
cost functional of the form, in a finite horizon feature,

$J(x,t,\alpha)=\int_{t}^{T}\ell(y(s),\alpha(s),s)ds,$

with the running cost $\ell$ suitably designed in order to keep trace of the
distances from the targets of the trajectory. The visiting cost is then
defined as

$v(x,t)=\inf_{\alpha\in\operatorname{\mathcal{A}}}J(x,t,\alpha).$

Actually, the problem requires a particular framework as a standard continuous
optimal control setting fails to describe the problem correctly. Let us
illustrate this difficulty using the following toy example.

Let us consider the $1$-D problem with $A={\mathbb{R}}$, $f(x,a)=a$, and
${\mathcal{T}}_{1}=(-\infty,-1]$, ${\mathcal{T}}_{2}=[1,+\infty)$. We focus on
designing an optimal visit formulation in the interval $[-1,1]$. At first, we
consider the easiest running cost design:

$\ell(y,a,t):=\frac{1}{2}\left(\sum_{j=1}^{N}d(y,\operatorname{\mathcal{T}}_{j})^{2}+\|a\|^{2}\right),$

which penalizes quadratically the distance from the targets and the norm of
the control. It is easy to verify that we have a feedback formula for the
optimal control as

$\alpha(x,t)=-\frac{x}{(1-t)^{2}+1/2},$

which means that the trajectory is led to zero, which is the middle point
between the two targets. Since we want to model a slightly different problem,
i.e., a visit more than a compromise between the distances, we are unsatisfied
with this result.

Therefore, we should include the information about the visit of the targets in
a different way in the model, allowing us to focus on a single target, as well
as on a subfamily of targets, at once. If, for example, we consider the
problem of visiting first target $\operatorname{\mathcal{T}}_{1}$ and then
$\operatorname{\mathcal{T}}_{2}$, we can easily observe that we would obtain a
different problem just swapping the order of the visit. This is a consequence
because no Dynamical Programming Principle would be available for the function
$v(x,t)$, since the only information brought by the state-position $x$ does
not give information about the already visited targets (see also [3]). Hence,
at this level, we can not in general characterize the optimal visiting time as
a solution of a Hamilton-Jacobi-Bellman equation. Consequently, it is quite
challenging to perform a global study of the problem or obtain a feedback
optimal control map.

The argument above suggests that we need to include in the model a “memory” of
the targets already visited. The latter can be done using various tools. Here,
we opt for a hybrid control based construction. In particular, to each target
$\operatorname{\mathcal{T}}_{j}$ we associate a label $p^{j}\in\\{0,1\\}$.
These labels indicate whether the corresponding target has been already
visited $p^{j}=1$, or not $p^{j}=0$, and they have a discontinuous evolution
in time. We then split the optimal visiting problem into several problems,
labeled by the $N$-strings $p\in\\{0,1\\}^{N}$, and we suitably interpret it
as a collection of several optimal stopping/switching problems coupled to each
other by the stopping/switching cost (a switch between $N$-strings corresponds
to the visit of a target or the choice to forgo visiting one or more targets).

In the following, we adapt to our framework some classic results of viscosity
solutions theory that can be found, e.g., in [9, 24]. In particular, the
hybrid framework that we propose is related to hybrid control [15] and somehow
to the mathematical switching hysteresis models [34]. The need of memory
feature, associated with the optimal visiting and dynamic programming and
Hamilton-Jacobi equations has been presented in [3], where a continuous
hysteresis memory was introduced. The use of a switching/discontinuous/hybrid
memory, as in the present paper, was instead used for a one-dimensional
optimal visiting problem on a network in [5]. For switching hybrid control
problems, as well as for differential games, related to the model here
presented, and in connection with Hamilton-Jacobi equations, we refer to [4]
and to [6]. A more general discussion is done in [12] (similar formulations
for the deterministic case have also been proposed in [15, 21]).

The literature concerning the Traveling Salesman Problem, to which our optimal
visiting problem is related, is very large. We only quote an early paper by R.
Bellman [11] devoted to the problem and dynamic programming.

In this paper, we are also interested in the case where a huge population of
agents plays the optimal visiting problem with controlled dynamics and costs
also depending on the distribution of the population. Actually, the study of
the interacting motion of many agents with more than one target seems to be
rather new in the literature, especially for what concerns the corresponding
continuity equation for the mass distribution. In this work, we start such a
kind of study and we furnish some numerical promising results and the analysis
for some continuity equations with a mass-sink (corresponding to the case
where some agents, possibly labeled by the same $N$-strings $p$, forgo
visiting some targets and then pass to another level, labeled by another
$N$-string $p^{\prime}$).

The model for a crowd of indistinguishable players is taken from the framework
of mean-field games [31, 29, 28, 18], while the adaptation of the same
structure to hybrid processes has been only very recently attempted. Some
related works are [13, 14], where the author discusses a mean-field optimal
stopping problem and [25], where a hybrid mean-field game is presented to
model a multi-lane traffic flux of vehicles.

As it is mentioned later in the paper, some complementary results, related to
a single player problem, are referred to [7], while in [8] the same framework
is used to solve a series of applied problem arising from the sport of
orienteering races.

The article is organized as follows. In Section 2, we introduce the optimal
visiting problem of a single agent, shortly reporting all the theoretical
elements that justify the use of a Hamilton-Jacobi formulation. A numerical
scheme to approximate the solution is described in Section 3, with a test to
verify in practice the model. In Section 4, which contains the main results of
the paper, we consider a crowd of indistinguishable players focusing on the
good position of the continuity equation that models the motion of the density
of players. Then, in Section 5, we introduce a numerical approximation for the
continuity equation and we examine the model in action through a collection of
simple tests.

##  2\. The optimal visiting problem

We are given of a collection of $N$ disjoint closed subsets
$\\{{\mathcal{T}}_{j}\\}_{j=1,\ldots,N}\subset{\mathbb{R}}^{d}$ and we
represent the state of the system by
$(x,p)\in{\mathbb{R}}^{d}\times{\mathcal{I}}$, where
$p=(p^{1},p^{2},\ldots,p^{N})\in{\mathcal{I}}=\\{0,1\\}^{N}$. Hence $x$ is the
continuous state variable and $p$ is the switching discrete state variable.
The evolution of the continuous variable is described by the controlled
dynamics

$\begin{cases}y^{\prime}(s)=f(y(s),\alpha(s),q(s)),\ s\in]t,T]\\\ y(t)=x,\
q(t)=p\end{cases},$ (2.1)

where $(x,p)\in{\mathbb{R}}^{d}\times{\mathcal{I}}$ is the initial state,
$t\in[0,T]$ the initial instant, $T>0$ the finite horizon. The measurable
control is (for $A\subset\mathbb{R}^{m}$ compact)

$\alpha\in{{\mathcal{A}}}:=\left\\{\alpha:[0,+\infty[\longrightarrow A\
\mbox{measurable}\right\\},$

and the dynamics $q(\cdot)$ of the switching variable is subject to

$\exists\tau\in[t,s],\ y(\tau)\in{\operatorname{\mathcal{T}}_{j}}\Rightarrow\
q^{j}(s)=1;\ \ q^{j}(s)=p^{j}\ \text{otherwise}.$

In particular, $q^{j}(s)=0$ means that the target ${\mathcal{T}}_{j}$ has not
been visited yet in $[t,s]$ and viceversa for $q^{j}(s)=1$. The dynamics
$f:{\mathbb{R}}^{d}\times A\times{\mathcal{I}}\longrightarrow{\mathbb{R}}^{d}$
is continuous, bounded and Lipschitz continuous w.r.t. $x\in{\mathbb{R}}^{d}$
uniformly w.r.t. $(a,p)\in A\times{\mathcal{I}}$, i.e., there exists $L>0$
such that

$\|f(x,a,q)-f(y,a,q)\|\leq L\|x-y\|\ \ \text{for all}\
(x,y)\in{\mathbb{R}}^{d}\ \text{and}\ (a,q)\in A\times{\mathcal{I}}.$

For every initial state $(x,t,p)$ and control $\alpha$, by our hypotheses the
existence of a unique solution $y(s)$ of (2.1) is guaranteed. Note that the
number of switches of the variable $q$ is necessarily finite.

The optimal visiting problem consists then in reaching the discrete state
$\bar{p}=(1,1,\dots,1)$ (i.e., to visit all the targets) at a time
$t\leq\bar{t}\leq T$, minimizing the following cost

$\int_{t}^{\bar{t}}e^{-\lambda(s-t)}\ell(y(s),{\alpha}(s),q(s),s)ds,$

for a given running cost $\ell$ and a discount factor $\lambda>0$.

###  2.1. A hybrid-control relaxation of the problem: optimal switching

The optimal control problem described above requires to “exactly touch” all
the targets. This makes the evolution of the discrete variable $q$ rather
complicated, in particular in view of the corresponding Hamilton-Jacobi
equation. We then relax the problem asking instead for “to pass as close as
possible” to each target. Then we assume that we can definitely get rid of
some targets and take into account only the remaining ones. In doing that, we
also pay an additional cost depending, for instance, on the actual distance
from the discarded targets. In this way, the evolution $q(\cdot)$ of the
discrete variables is no more a solution of system (2.1), but instead, it
becomes a control at our disposal. Obviously, there are some constraints: for
example, for $N=4$, if $p=(1,0,0,0)$, $p^{\prime}=(1,1,0,0)$,
$p^{\prime\prime}=(0,1,1,0)$ and $p^{\prime\prime\prime}=(1,1,1,0)$, then from
$p$ we can not switch to $p^{\prime\prime}$ otherwise we lose the information
about the already visited/discarded target ${\mathcal{T}}_{1}$. However, we
can switch to $p^{\prime\prime\prime}$ directly.

Therefore, for any $p$, we denote by ${\mathcal{I}}_{p}$ the set of all
possible new variables in ${\mathcal{I}}$ after a switch from $p$:

${\mathcal{I}}_{p}=\\{\tilde{p}\in{\mathcal{I}}:p^{i}=1\Rightarrow\tilde{p}^{i}=1\text{
and}\ \exists l=1,\ldots,N:p^{l}=0,\;\tilde{p}^{l}=1\\}.$

Note that in particular ${\mathcal{I}}_{\bar{p}}=\emptyset$, where
$\bar{p}=(1,1,\ldots 1)$.

For a given $p$, the number of the admissible subsequent switches is at most
$N-\sum_{i}p^{i}\leq N$. Given the state $(x,p)$ at the time $t$ with
$p\neq\bar{p}$, the controller chooses: the measurable control
${\alpha}\in{\mathcal{A}}$, and the discrete one
$q:[0,+\infty[\longrightarrow{\mathcal{I}}$ which contains: the number $1\leq
m\leq N-\sum_{i}p^{i}$ of switches to be performed in order to reach
$\bar{p}$, the switching instants $t\leq t_{1}<t_{2}<\ldots<t_{m}\leq T$ and
the switching destinations $p_{1},\ldots,p_{m-1}$, $p_{m}=\bar{p}$. Such
destinations must satisfy $p_{1}\in{\mathcal{I}}_{p}$,
$p_{i+1}\in{\mathcal{I}}_{p_{i}}$, $i=1,\ldots,m-1$. To resume, the control at
disposal is then

$({\alpha},m,t_{1},\ldots,t_{m},p_{1},\ldots,p_{m-1})=:u,$

and note that for any $(x,p,t)$ as above such a string belongs to a set
depending on $p$ and $t$ denoted by ${\mathcal{U}}_{(p,t)}$. The cost to be
minimized is

$J(x,t,p,u)=\sum_{j=1}^{m}\Bigg{(}\int_{t_{j-1}}^{t_{j}}e^{-\lambda(s-t)}\ell(y(s),{\alpha}(s),p_{j-1},s)ds\\\
+e^{-\lambda(t_{j}-t)}C(y(t_{j}),p_{j-1},p_{j})\Bigg{)},$

with $\lambda\geq 0$, $p_{0}=p$, $t_{0}=t$ and $y(s)$ is the solution of (2.1)
where $q(s)=p_{j-1}$ if $s\in[t_{j-1},t_{j}]$.

We assume $\ell:{\mathbb{R}}^{d}\times
A\times{\mathcal{I}}\times[0,T]\longrightarrow[0,+\infty[$ bounded, continuous
and uniformly continuous w.r.t. $x$ uniformly w.r.t. $a\in A$,
$p\in{\mathcal{I}}$ and $t\in[0,T]$. Moreover
$C:{\mathbb{R}}^{d}\times{\mathcal{I}}\times{\mathcal{I}}\longrightarrow[0,+\infty[$
is uniformly continuous w.r.t. $x\in{\mathbb{R}}^{d}$, uniformly w.r.t.
$p,p^{\prime}\in{\mathcal{I}}\times{\mathcal{I}}_{p}$. Note that
$C(x,p,p^{\prime})$ represents the switching cost from $p$ to $p^{\prime}$
when the state position is $x\in{\mathbb{R}}^{d}$. For example, it may depend
on the distance from the discarded targets, that is
$C(x,p,p^{\prime})=\sum_{j}\chi_{j}(p,p^{\prime})d(x,{\mathcal{T}}_{j})$,
where

$\chi_{j}(p,p^{\prime})=\begin{cases}0,&p^{j}=p^{\prime j}\\\
1,&\text{otherwise}\end{cases}.$

The value function of the problem is

$V(x,t,p)=\inf_{u\in{\mathcal{U}}_{(p,t)}}J(x,t,p,u).$ (2.2)

###  2.2. Another possible interpretation: a family of optimal stopping
problems

We briefly present another possible characterization of the above optimal
switching problem. In particular, we divide it into several time-dependent
optimal stopping subproblems, one per every switching variable $p$, suitably
coupled by the stopping costs. To do this, we proceed with a backward
construction, whose details are explained in [7]. Anyway, we give here a short
exposure. We start from the last switching variable, that is, the variable $p$
such that $\sum_{i}p^{i}=N-1$ (i.e., from $p$ we can switch only to
$\bar{p}$). At every step, we get a time-dependent optimal stopping problem in
the state space ${\mathbb{R}}^{d}$ where, for a given $(x,t)$, the admissible
controls are the triples
$u=({\alpha}\in{\mathcal{A}},\tau\in[t,T],p^{\prime}\in{\mathcal{I}}_{p})$ and
the cost to be minimized is

$J_{p}(x,t,{\alpha},\tau,p^{\prime})=\int_{t}^{\tau}e^{-\lambda(s-t)}\ell(y(s),{\alpha}(s),p,s)ds\\\
+e^{-\lambda(\tau-t)}\Big{(}C(y(\tau),p,p^{\prime})+V_{p^{\prime}}(y(\tau),\tau)\Big{)}.$

The value function is

$V_{p}(x,t)=\inf_{({\alpha},\tau,p^{\prime})}J_{p}(x,t,{\alpha},\tau,p^{\prime}).$
(2.3)

Note that, if $p$, for instance, is such that $\sum_{i}p^{i}=N-2$, then from
$p^{\prime}$ we can only switch to the final state $\bar{p}$, and hence
$V_{p^{\prime}}$ can be a priori evaluated as in the previous step of the
procedure. Since when $p=\bar{p}$, the game stops, we set $V_{\bar{p}}\equiv
0$. Therefore, proceeding backwardly in this way, we can at least formally
compute the value functions $V_{p}$ for any $p\in{\mathcal{I}}$.

###  2.3. Equivalence of the two models

We show the equivalence between the optimal switching problem and the family
of the optimal stopping ones, i.e., $V(x,t,p)=V_{p}(x,t)$ for every
$(x,t,p)\in{\mathbb{R}}^{d}\times[0,T]\times{\mathcal{I}}$. Here, and in the
sequel, $V$ is the value function defined in (2.2) and $V_{p}$ is the value
function defined backwardly as in (2.3). For the proofs of all the results in
this section, we refer to [7].

######  Proposition 1.

Under the hypotheses of §2.1, we have that

  * $(i)$

$V$ and $V_{p}$ are bounded and uniformly continuous for every
$p\in{\mathcal{I}}$;

  * $(ii)$

for every $x\in\mathbb{R}^{d}$, $t\in[0,T]$ and $p\in{\mathcal{I}}$,

$V(x,t,p)=\inf_{({\alpha},\tau,p^{\prime}\in{\mathcal{I}}_{p})}\Bigg{(}\int_{t}^{\tau}e^{-\lambda(s-t)}\ell(y(s),{\alpha}(s),p,s)ds\\\
+e^{-\lambda(\tau-t)}\Big{(}C(y(\tau),p,p^{\prime})+V_{p^{\prime}}(y(\tau),\tau)\Big{)}\Bigg{)}.$

As a consequence, $V(x,t,p)=V_{p}(x,t)$ for all $(x,t,p)$.

We are able to obtain a differential characterization of the value function
$V_{p}$ as viscosity solution of an Hamilton-Jacobi-Bellman variational
inequality.

In the sequel, by $(\cdot)_{t}$ and $D_{x}$ we denote the time derivative and
the spatial gradient.

######  Theorem 2.1.

Under the hypotheses of Proposition 1, for every $p\in{\mathcal{I}}$ the value
function $V_{p}$ is the unique bounded and uniformly continuous viscosity
solution $U$ of

$\begin{cases}\max\\{U(x,t)-\psi_{p}(x,t),-U_{t}(x,t)+\lambda
U(x,t)+H^{p}(x,t,D_{x}U(x,t))\\}=0,\\\ \mbox{\hskip
233.3125pt$(x,t)\in\mathbb{R}^{d}\times[0,T[$}\\\
U(x,T)=\psi_{p}(x,T),\mbox{\hskip
193.47882pt$x\in\mathbb{R}^{d}$},\end{cases}$ (2.4)

where

$H^{p}(x,t,\xi)=\sup_{a\in A}\\{-f(x,a,p)\cdot\xi-\ell(x,a,p,t)\\}$ (2.5)

is the $p$-labeled Hamiltonian and

$\psi_{p}(x,t):=\inf_{p^{\prime}\in{\mathcal{I}}_{p}}(C(x,p,p^{\prime})+V_{p^{\prime}}(x,t))$
(2.6)

is the $p$-labeled switching operator. Moreover, the family of functions
$\\{V_{p}:p\in{\mathcal{I}}\\}$ is the unique family of bounded and uniformly
continuous functions $\\{U_{p}:p\in{\mathcal{I}}\\}$ that solves the problem

$\left\\{\begin{array}[]{ll}\displaystyle\text{for any $p\in{\mathcal{I}}$,
$U_{p}$ is the unique viscosity solution of \eqref{V}}\\\
\displaystyle\text{with }\psi_{p}\ \mbox{replaced by
}\psi_{p}^{U}(x,t):=\inf_{q^{\prime}\in{\mathcal{I}}}(C(x,p,p^{\prime})+U_{p^{\prime}}(x,t)),\\\
\displaystyle U_{\bar{p}}=0\end{array}.\right.$

##  3\. Numerical testing

In this section, we briefly describe a technique of approximation for the
equation (2.4). This will be useful to produce some tests to illustrate the
framework and show in practice the advantages of the method. Moreover, the
same approximation scheme will be used in the mean-field game framework.

###  3.1. Numerical approximation

We describe a numerical approximation for (2.4). Let us consider a discrete
grid of nodes $(x_{i},p)$ in the state space that, for simplicity, we consider
as $[-c,c]^{2}\times{\mathcal{I}}$, with $c\in{\mathbb{R}}_{+}$. Thus, a
general $i=(i_{1},i_{2})$, $x_{i}=(-c+i_{1}\Delta x,-c+i_{2}\Delta x)$. We
proceed to define also a discrete time space as $t_{k}:=k\Delta t\in[0,T]$.

Following [23], we write an approximation scheme for (2.4) at $(x_{i},p)$ in a
backward in time explicit form as

$V(x_{i},t_{k-1},p)=\min\left(NV(x_{i},t_{k},p),\Sigma\left(x_{i},t_{k},p,V\right)\right).$
(3.1)

In (3.1), the numerical operator $\Sigma$ is related to the continuous
control, or, in other terms, to the approximation of the Hamiltonian function
(2.5).

We approximate the Hamiltonian part using a semi-Lagrangian approach (see,
i.e., [16]). The main advantage of such an approach is that we build a
monotone scheme unconditional stable to the grid parameters justifying the
choice of an explicit-in-time scheme. A standard semi-Lagrangian
discretization of the Hamiltonian is given by

$\Sigma\left(x_{i},t_{k},p,V\right)=\lambda\Delta tV(x_{i},t_{k},p)+\min_{a\in
A}\left\\{\mathbb{I}\left[V\right](x_{i}-\Delta
t\>f(x_{i},a,p),t_{k},p)\right.\\\ \left.+\Delta
t\,\ell(x_{i},a,p,t_{k})\right\\},$ (3.2)

where we denote by $\mathbb{I}\left[V^{\Delta}\right](x,t,p)$ the values of
the discrete function $V^{\Delta}$ computed at $(x,t,p)$ obtained using a
standard interpolation operator. If the interpolation $\mathbb{I}$ is
monotone, the resulting scheme is consistent, monotone and $L^{\infty}$
stable, and therefore convergent via Barles-Souganidis theorem (see [10]).
Some typical examples of monotone interpolation operators are
$\operatorname{\mathbb{P}}_{1}$ (piecewise linear on triangles/tetrahedra) and
$\operatorname{\mathbb{Q}}_{1}$ (piecewise multilinear on rectangles)
interpolations.

The discrete switch operator $N$ is computed at a node $(x_{i},p)$ as

$NV(x_{i},t_{k},p):=\min_{p^{\prime}\in\bar{\mathcal{I}}^{p}}\left\\{C(x_{i},p,p^{\prime})+V(x_{i},t_{k},p^{\prime})\right\\},$
(3.3)

which is clearly, the same operator introduced before in (2.6) with the only
difference that, being the operator backward in time, we look in the
complementary of the admissible switches
$\bar{\mathcal{I}}^{p}={\mathcal{I}}\setminus{\mathcal{I}}^{p}$.

###  3.2. An optimal visiting test for a single agent

We illustrate the technique on a simple scenario. We consider a classic
quadratic penalization of the control in the running cost and some isotropic
dynamics

$\ell(x,a,p,t)=\frac{\|a\|^{2}}{2},\quad f(x,a,p)=a,$

The control is chosen in the whole $A={\mathbb{R}}^{2}$ (despite the fact that
for classic results it is possible to prove that the optimal control is always
included in a compact subset of ${\mathbb{R}}^{2}$, cf. [9]) and the problem
lives in the domain $\Omega=[-1,1]^{2}$. The target set is composed by the
points

${\mathcal{T}}_{j}=0.6\,\left(\cos\left(j\frac{2\pi}{3}\right),\sin\left(j\frac{2\pi}{3}\right)\right),\quad
j=1,\ldots,3.$

The switching cost is set to

$C(x,p,p^{\prime})=\sum_{j\in{\mathcal{I}}}\chi_{j}(p,p^{\prime})\|x-{\mathcal{T}}_{j}\|,$

where $\chi_{j}(p,p^{\prime})$ is an indicator of the possible switch relative
to the $j$-target (if the player renounces to visit the $j$-target, it pays
the distance from it). We underline that multiple switches are allowed by this
framework:

$\chi_{j}(p,p^{\prime})=\begin{cases}0&\text{if }p_{j}=p^{\prime}_{j}\\\
1&\text{otherwise}\end{cases}.$

We implemented the scheme (3.2), where the optimization is performed by
steepest descent on the control set $A$. The final time is $T=5$ while the
boundary condition penalizes the states different from the final one
$\bar{p}$, i.e.

$\psi_{p}(x,T)=\sum_{j\in{\mathcal{I}}}\chi_{j}(p,\bar{p})\|x-{\mathcal{T}}_{j}\|\quad\hbox{
for any }p\in{\mathcal{I}},\ x\in[-1,1]^{2}.$

The latter means that if a trajectory ends not in the final state $\bar{p}$
where all targets have been considered, it pays the cost of switching to
$\bar{p}$. This choice is made to preserve continuity on the boundary of the
value function for each state.

---
Figure 1. Test 1. Approximated value functions in the various discrete states
of the system

We notice that the optimal control model that we have built tries to make any
trajectory as closely as possible to all, or to the higher number of unvisited
target points on the available time $T-t$, minimizing the effort $\|\alpha\|$.

Figure 1 shows a collection of approximated value functions at time $T=0$ for
any state of the system. The active targets (the ones still labeled as $1$ in
$p$) act as attractors for the trajectory (local minimum in the value
function) while they do not are highlighted if not active.

---
Figure 2. Test 1. Optimal trajectories for various starting points: (top/left)
$x_{0}=(0,0)$ (top/right) $x_{0}=(0,-0.2)$, (bottom) $x_{0}=(0.9,0.9)$.

During the computation of the value function, we found the optimal control and
the switching maps as arguments of the minimization (3.2), (3.3). We call them
respectively $\alpha^{*}(x,t,p)$ and $\sigma^{*}(x,t,p)$. We build an
approximation of an optimal trajectory starting from a point
$x_{0}\in[-1,1]^{2}$ via the standard Euler approximation

$\begin{cases}Y_{x_{0}}(t^{n+1}):=Y_{x_{0}}(t^{n})+\Delta
t\alpha^{*}(Y_{x_{0}}(t^{n}),t^{n},p^{n})\\\
p^{n+1}:=\begin{cases}p^{*}&\hbox{ if
}\sigma^{*}(Y_{x_{0}}(t^{n}),t^{n},p^{n})=p^{*}\\\ p^{n}&\hbox{
otherwise}\end{cases}\\\ Y_{x_{0}}(0)=x_{0}\\\ p^{0}=\hat{p}\end{cases}.$

In Figure 2, we show some of these optimal trajectories for various choices of
the starting points. We observe as, when possible, the path visits all the
targets with an order that is determined by the starting point. When this is
not possible (as in the case of $x_{0}=(0.9,0.9)$) the trajectory gets as
close as possible to each target before renouncing and aiming at the next one.
Outside the switching point, straight lines constitute the trajectories since
the running cost is spatially and temporally homogeneous.

##  4\. Optimal visiting for a crowd: the continuity equation

In a possible study of a mean-field game for a population of agents, each one
of them playing a $p$-labeled optimal stopping problem like the one in Section
2, one would be led to consider the coupling of the system (2.4) of Hamilton-
Jacobi quasi-variational inequalities (coupled by the stopping costs) with a
system of continuity equations (one per each level $p$ and coupled by a
transfer through some sinks and sources). In particular, the sink at the level
$p$ is the region where the agents stop running at level $p$ and pass to a new
subsequent level $p^{\prime}\in{\mathcal{I}}_{p}$, and similarly for the
sources. Such a coupling should provide the optimal flow direction and the
optimal switching maps, for the evolution of the masses of agents, that are
also labeled by the variable $p$.

As said in the Introduction, here, from an analytical point of view, we focus
on a single continuity equation for a given suitably regular field, with
possible sinks and sources and we left further analysis to future studies. In
Section 5, some numerical tests are shown for more general situations.

###  4.1. The continuity equation with a sink

We want to model the evolution on ${\mathbb{R}}^{d}$ of a mass $m$ subject to
a given flow with a presence of a given region of ${\mathbb{R}}^{d}$ acting as
sink: the portion of mass that possibly enters the sink instantaneously
disappears. The next section will give some numerical experiments. The region
representing the sink can be in general moving in time but, here, we regard it
as constant, the generalization to the moving case working with same ideas and
calculations.

For the notation and the construction of the following setting, we mostly rely
on [17]. We consider a flow
$\Phi:{\mathbb{R}}^{d}\times{\mathbb{R}}\times{\mathbb{R}}\longrightarrow{\mathbb{R}}^{d}$
given by the solutions of the ordinary differential system for
$t\in{\mathbb{R}}$ and $x\in{\mathbb{R}}^{d}$,

$\begin{cases}y^{\prime}(s)=b(y(s),s),&s>t\\\ y(t)=x\end{cases},$ (4.1)

that is $\Phi(x,t,s)=y(s)$ solving (4.1). We will be mostly concerned with
$\Phi(\cdot,0,\cdot)$. In (4.1), the field
$b:{\mathbb{R}}^{d}\times{\mathbb{R}}\longrightarrow{\mathbb{R}}^{d}$ is
assumed to be bounded, continuous and Lipschitz continuous w.r.t.
$x\in{\mathbb{R}}^{d}$ uniformly w.r.t. $t\in{\mathbb{R}}$. Then, the flow
$\Phi(\cdot,0,\cdot)$ is Lipschitz continuous.

The sink is represented by a subset ${\mathcal{S}}\subset{\mathbb{R}}^{d}$,
which is assumed to be closed with compact and $C^{1}$ boundary. Then, for a
point $x\in{\mathbb{R}}^{d}$, we define the possible first arrival time (to
the sink) as

$t_{x}:=\inf\\{t\geq 0:\Phi(x,0,t)\in{\mathcal{S}}\\},$ (4.2)

and the set of possible arrival points on the sink, for a given $t$, as

${\mathcal{S}}^{t}:=\\{z\in\partial{\mathcal{S}}:\exists x\in{\mathbb{R}}^{d}\
\text{such that}\ t_{x}=t\ \text{and}\ \Phi(x,0,t)=z\\}.$ (4.3)

We will see in §4.3 that it is possible to characterize the possible first
arrival time $t_{x}$ as the unique (viscosity) solution of an Hamilton-Jacobi-
Bellman equation with suitable boundary conditions.

We work in the set ${\mathcal{G}}$ of positive Radon measures $\mu$ on
${\mathbb{R}}^{d}$ with finite first order moment, bounded by a constant $G$
(i.e., $\int_{{\mathbb{R}}^{d}}d\mu\leq G$ for every $\mu\in{\mathcal{G}}$).
Such a space can be endowed with the generalized Wasserstein distance (see
[32, 33])

${\mathcal{W}}(\mu,\mu^{\prime})=\sup\left\\{\int_{{\mathbb{R}}^{d}}\varphi
d(\mu-\mu^{\prime}):\varphi\in C^{0}_{c}({\mathbb{R}}^{d}),\
\|\varphi\|_{\infty}\leq 1,\ \operatorname{Lip}(\varphi)\leq 1\right\\}.$
(4.4)

Let

$\tilde{m}_{0}(x):=\begin{cases}m_{0}(x),&x\in{\mathbb{R}}^{d}\setminus{\mathcal{S}}\\\
0,&\text{otherwise}\end{cases},\qquad\tilde{m}_{0}\in{\mathcal{G}},$

where $m_{0}\in{\mathcal{G}}$ is given, be the initial distribution on
${\mathbb{R}}^{d}$. We suppose that it is absolutely continuous with a
density, still denoted $\tilde{m}_{0}$, which is bounded and has a compact
support.

Then, the continuity equation with a sink (to be interpreted in a suitable
weak formulation), with a finite horizon $T>0$, is

$\begin{cases}\mu_{t}(x,t)+\operatorname{div}(\mu(x,t)b(x,t))+\mathbbm{1}_{\\{({\mathcal{S}}^{t},t):t\in[0,T]\\}}\mu(x,t)=0,\\\
\mbox{\hskip 227.33766pt$(x,t)\in\mathbb{R}^{d}\times[0,T]$}\\\
\mu(x,0)=\tilde{m}_{0}(x),\mbox{\hskip
204.0062pt$x\in\mathbb{R}^{d}$}\end{cases}.$ (4.5)

######  Remark 1.

In the possible mean field game problem, using the notation of previous
sections, at every level $p$ the sink would be given by the evolutive stopping
set ${\mathcal{S}}_{p}(t)=\\{x\in{\mathbb{R}}^{d}:V_{p}(x,t)=\psi_{p}(x,t)\\}$
and the field $b$ by the gradient of the value function $V_{p}$. Hence, the
regularity assumptions above should be probably adjusted.

In view of the mean field case, in §4.2 we will study a possible dependence of
the field $b$ on the measure and in §5 some numerical tests will be performed
including this possibility.

In the sequel, we denote by $\Psi$ the inverse of the flow $\Phi$ starting
from $\partial{\mathcal{S}}$, i.e., all the states backwardly reached by the
trajectories starting from the points of $\partial{\mathcal{S}}$ in the time
interval $[0,T]$. That is

$(z,\tau)\longmapsto\Psi(z,\tau,\tau),\quad 0\leq\tau\leq T,\ \
z\in\partial{\mathcal{S}}$

with $\Psi(z,\tau,\tau)=\zeta(\tau)$ satisfying

$\begin{cases}\zeta^{\prime}(s)=-b(\zeta(s),\tau-s),&0<s<\tau\\\
\zeta(0)=z\end{cases}.$ (4.6)

By hypotheses, $\Psi$ is Lipschitz continuous as $\Phi$ and it is such that

$\displaystyle\Phi(\Psi(z,\tau,\tau),0,\tau)$ $\displaystyle=z,$ (4.7)
$\displaystyle\Psi(\Phi(x,0,\tau),\tau,\tau)$ $\displaystyle=x.$

Now, fixed $s\in[0,T]$, we define the sink-reaching-points set, in time $s$,
as

${\mathcal{B}}(s):=\\{x\in{\mathbb{R}}^{d}:t_{x}\leq s\\},$

that is the set of all initial points $x\in{\mathbb{R}}^{d}$ from which the
agents enter the sink before $s$. Observe that

${\mathcal{B}}(s)=\bigcup_{\tau\in[0,s]}\underbrace{\\{x\in{\mathbb{R}}^{d}:t_{x}=\tau\\}}_{=:{\mathcal{B}}^{\tau}}=\bigcup_{\tau\in[0,s]}\\{x=\Psi(z,\tau,\tau):z\in{\mathcal{S}}^{\tau}\\}$

and that

$s_{1}\leq s_{2}\ \Rightarrow\
{\mathcal{B}}(s_{1})\subseteq{\mathcal{B}}(s_{2}).$ (4.8)

######  Definition 4.1.

We say that $\mu$ is a weak solution to (4.5) if $\mu\in
L^{1}([0,T],{\mathcal{G}})$ is such that, for any test function $\varphi\in
C^{\infty}_{c}({\mathbb{R}}^{d}\times[0,T[)$, we have

$\int_{{\mathbb{R}}^{d}}\varphi(x,0)d\tilde{m}_{0}(x)\\\
+\int_{0}^{T}\int_{{\mathbb{R}}^{d}\setminus\Phi({\mathcal{B}}(t),0,t)}(\varphi_{t}(x,t)+\langle
D_{x}\varphi(x,t),b(x,t)\rangle)d\mu(t)(x)dt\\\
-\int_{0}^{T}\int_{{\mathbb{R}}^{d}}\mathbbm{1}_{\\{({\mathcal{S}}^{t},t):t\in[0,T]\\}}\varphi(x,t)d\mu^{t}(t)(x)dt=0,$

where $\mu^{t}$ entering the last integral is $\mu^{t}(t)=g(t)\mu^{t}(0)$. The
measure $\mu^{t}(0)$ is the disintegration of $\mu(0)$ on the fibers
${\mathcal{B}}^{\tau}$ that compose ${\mathcal{B}}(t)$ and $g(\cdot)$ is the
density of the measure $\nu$ on the indices $\tau$ of the fibers
${\mathcal{B}}^{\tau}$ such that

$E\subset{\mathcal{B}}(t)\ \Rightarrow\
\mu(0)(E)=\int_{0}^{t}\mu^{\tau}(0)({\mathcal{B}}^{\tau}\cap
E)d\nu=\int_{0}^{t}g(\tau)\mu^{\tau}(0)({\mathcal{B}}^{\tau}\cap E)d\tau.$

For $s\in[0,T]$, we define the following measure on ${\mathbb{R}}^{d}$

$\tilde{\mu}(s)=\begin{cases}\Phi(\cdot,0,s)\\#\tilde{m}_{0}&\text{on}\
{\mathbb{R}}^{d}\setminus\Phi({\mathcal{B}}(s),0,s)\\\
0,&\text{otherwise}\end{cases}.$ (4.9)

Observe that the set $\Phi({\mathcal{B}}(s),0,s)$ takes into account all the
agents who passed through the sink ${\mathcal{S}}$ at least once in the time
interval $[0,T]$ and then disappeared. For simplicity, we set
$\mu(s):=\Phi(\cdot,0,s)\\#\tilde{m}_{0}$.

By the hypotheses on $b$, $\tilde{m}_{0}$ and by (4.9), we have that
$\tilde{\mu}(s)$ is a positive Radon measure on ${\mathbb{R}}^{d}$ with finite
first order moment. Moreover, it certainly satisfies the constraint
$\int_{{\mathbb{R}}^{d}}d\tilde{\mu}(s)\leq G$ because, with respect to
$\mu(s)$, it may only lose mass through the sink. Hence
$\tilde{\mu}(s)\in{\mathcal{G}}$. Furthermore, it is absolutely continuous
with a density which is bounded and has a compact support.

######  Lemma 4.2.

For $s\in[0,T]$, consider the function

$\pi:{\mathcal{B}}(s)\longrightarrow[0,s],\qquad\pi(x):=t_{x}$ (4.10)

and suppose that it is Lipschitz continuous (see Remark 2). Then, the map
$s\longmapsto\tilde{\mu}(s)$ is Lipschitz continuous in ${\mathcal{G}}$ (with
respect to the metrics (4.4)).

###### Proof.

Let $s,s^{\prime}\in[0,T]$, $s\geq s^{\prime}$. Then, recalling (4.8), we have

${\mathcal{W}}(\tilde{\mu}(s^{\prime}),\tilde{\mu}(s))=\sup_{\begin{subarray}{c}\|\varphi\|_{\infty}\leq
1\\\ \operatorname{Lip}(\varphi)\leq
1\end{subarray}}\Bigg{\\{}\int_{{\mathbb{R}}^{d}}\varphi(x)d\tilde{\mu}(s^{\prime})(x)-\int_{{\mathbb{R}}^{d}}\varphi(x)d\tilde{\mu}(s)(x)\Bigg{\\}}\\\
=\sup_{\begin{subarray}{c}\|\varphi\|_{\infty}\leq 1\\\
\operatorname{Lip}(\varphi)\leq
1\end{subarray}}\Bigg{\\{}\int_{{\mathbb{R}}^{d}\setminus\Phi({\mathcal{B}}(s^{\prime}),0,s^{\prime})}\varphi(x)d\mu(s^{\prime})(x)-\int_{{\mathbb{R}}^{d}\setminus\Phi({\mathcal{B}}(s),0,s)}\varphi(x)d\mu(s)(x)\Bigg{\\}}\\\
=\sup_{\begin{subarray}{c}\|\varphi\|_{\infty}\leq 1\\\
\operatorname{Lip}(\varphi)\leq
1\end{subarray}}\Bigg{\\{}\int_{{\mathbb{R}}^{d}}\varphi(x)d(\mu(s^{\prime})-\mu(s))(x)+\int_{\Phi({\mathcal{B}}(s),0,s)\setminus\Phi({\mathcal{B}}(s^{\prime}),0,s)}\varphi(x)d\mu(s)(x)\\\
+\int_{\Phi({\mathcal{B}}(s^{\prime}),0,s)}\varphi(x)d\mu(s)(x)-\int_{\Phi({\mathcal{B}}(s^{\prime}),0,s^{\prime})}\varphi(x)d\mu(s^{\prime})(x)\Bigg{\\}}\\\
=\sup_{\begin{subarray}{c}\|\varphi\|_{\infty}\leq 1\\\
\operatorname{Lip}(\varphi)\leq
1\end{subarray}}\Bigg{\\{}\int_{{\mathbb{R}}^{d}}(\varphi(\Phi(x,0,s^{\prime}))-\varphi(\Phi(x,0,s)))d\tilde{m}_{0}+\\\
\int_{{\mathcal{B}}(s)\setminus{\mathcal{B}}(s^{\prime})}\varphi(\Phi(x,0,s))d\tilde{m}_{0}(x)+\int_{{\mathcal{B}}(s^{\prime})}(\varphi(\Phi(x,0,s))-\varphi(\Phi(x,0,s^{\prime})))d\tilde{m}_{0}(x)\Bigg{\\}}\\\
\leq\sup_{\begin{subarray}{c}\|\varphi\|_{\infty}\leq 1\\\
\operatorname{Lip}(\varphi)\leq
1\end{subarray}}\Bigg{\\{}\int_{{\mathbb{R}}^{d}}\|\Phi(x,0,s^{\prime})-\Phi(x,0,s)\|d\tilde{m}_{0}+\int_{{\mathcal{B}}(s)\setminus{\mathcal{B}}(s^{\prime})}\varphi(\Phi(x,0,s))d\tilde{m}_{0}(x)\\\
+\int_{{\mathcal{B}}(s^{\prime})}\|\Phi(x,0,s)-\Phi(x,0,s^{\prime})\|d\tilde{m}_{0}(x)\Bigg{\\}}\\\
\leq
2GM|s^{\prime}-s|+\|\tilde{m}_{0}\|_{\infty}{\mathcal{L}}^{d}({\mathcal{B}}(s)\setminus{\mathcal{B}}(s^{\prime})),$

where $M$ is the time-Lipschitz constant for $\Phi$ (i.e., the bound for $b$).
Hence we conclude if we estimate
${\mathcal{L}}^{d}({\mathcal{B}}(s)\setminus{\mathcal{B}}(s^{\prime}))$. In
particular, if we prove that the map
$s\longmapsto{\mathcal{L}}^{d}({\mathcal{B}}(s))$ is Lipschitz continuous, we
are done. Observe that the function

$f:\sigma\longmapsto\pi(\Phi(y,0,\sigma))=\pi(y)-\sigma$

is such that

$1=|f^{\prime}(0)|=|\nabla\pi(y)\cdot
b(y,0)|\leq\|\nabla\pi\|\|b\|_{\infty}\leq M\|\nabla\pi\|\quad\text{a.e.}$

and then $\|\nabla\pi\|\geq\frac{1}{M}>0$ almost everywhere. Therefore

$\frac{{\mathcal{L}}^{d}({\mathcal{B}}(s)\setminus{\mathcal{B}}(s^{\prime}))}{M}\leq\int_{{\mathcal{B}}(s)\setminus{\mathcal{B}}(s^{\prime})}\|\nabla\pi\|dx=\int_{s^{\prime}}^{s}{\mathcal{H}}^{d-1}(\pi^{-1}(\tau))d\tau\\\
=\int_{s^{\prime}}^{s}{\mathcal{H}}^{d-1}(\Psi({\mathcal{S}}^{\tau},\tau,\tau))d\tau\leq
K|s-s^{\prime}|,$ (4.11)

where we used the Coarea Formula (see [22]) and the fact that the
$(d-1)$-dimensional Hausdorff measure
${\mathcal{H}}^{d-1}(\Psi({\mathcal{S}}^{\tau},\tau,\tau))$ is bounded by a
constant $K>0$ since, by hypotheses, $\Psi$ is Lipschitz continuous and
${\mathcal{S}}^{\tau}$ is compact. Then, the map
$s\longmapsto{\mathcal{L}}^{d}({\mathcal{B}}(s))$ is Lipschitz continuous and
the thesis follows. ∎

######  Remark 2.

In general, the map $\pi$ is not Lipschitz continuous in ${\mathcal{B}}(s)$.
But, in view of the possible mean field game model, the field $b$ will be the
optimal feedback for an optimal control problem with controlled dynamics from
a system $y^{\prime}={\alpha}$, and hence with total controllability. Then we
expect that such a Lipschitz continuity may hold and, at the moment, it is not
too restrictive to assume it. Anyway, future investigations will be made on
this direction.

######  Theorem 4.3.

The map $s\longmapsto\tilde{\mu}(s)$ is a weak solution of (4.5).

###### Proof.

Let $\varphi\in{\mathcal{C}}^{\infty}_{c}({\mathbb{R}}^{d}\times[0,T[)$. By
Lemma 4.2, the map

$s\longmapsto\int_{{\mathbb{R}}^{d}}\varphi(x,s)d\tilde{\mu}(s)(x)$

is absolutely continuous and then we have

$\frac{d}{ds}\int_{{\mathbb{R}}^{d}}\varphi(x,s)d\tilde{\mu}(s)(x)=\frac{d}{ds}\int_{{\mathbb{R}}^{d}\setminus\Phi({\mathcal{B}}(s),0,s)}\varphi(x,s)d\mu(s)(x)\\\
=\frac{d}{ds}\int_{{\mathbb{R}}^{d}}\varphi(\Phi(x,0,s),s)d\tilde{m}_{0}(x)-\frac{d}{ds}\int_{{\mathcal{B}}(s)}\varphi(\Phi(x,0,s),s)d\tilde{m}_{0}(x)\\\
=\int_{{\mathbb{R}}^{d}}(\varphi_{s}(\Phi(x,0,s),s)+\langle
D_{x}\varphi(\Phi(x,0,s),s),b(\Phi(x,0,s),s)\rangle)d\tilde{m}_{0}(x)\\\
-\frac{d}{ds}\int_{{\mathcal{B}}(s)}\varphi(\Phi(x,0,s),s)d\tilde{m}_{0}(x)\\\
=\int_{{\mathbb{R}}^{d}}(\varphi_{s}(y,s)+\langle
D_{x}\varphi(y,s),b(y,s)\rangle)d\mu(s)(y)\\\
-\frac{d}{ds}\int_{{\mathcal{B}}(s)}\varphi(\Phi(x,0,s),s)d\tilde{m}_{0}(x).$

We have to compute

$\frac{d}{ds}\int_{{\mathcal{B}}(s)}\varphi(\Phi(x,0,s),s)d\tilde{m}_{0}(x).$

By the Disintegration Theorem (see Remark 3), we get

$\frac{d}{ds}\int_{{\mathcal{B}}(s)}\varphi(\Phi(x,0,s),s)d\tilde{m}_{0}(x)\\\
=\frac{d}{ds}\int_{0}^{s}\int_{\\{x\in{\mathbb{R}}^{d}:t_{x}=\tau\\}}\varphi(\Phi(x,0,s),s)d\tilde{m}_{0}^{\tau}(x)d\nu(\tau)\\\
=\frac{d}{ds}\int_{0}^{s}\int_{\\{x\in{\mathbb{R}}^{d}:t_{x}=\tau\\}}\varphi(\Phi(x,0,s),s)g(\tau)d\tilde{m}_{0}^{\tau}(x)d\tau\\\
=\int_{\\{x\in{\mathbb{R}}^{d}:t_{x}=s\\}}\varphi(\Phi(x,0,s),s)g(s)d\tilde{m}_{0}^{s}(x)\\\
+\int_{{\mathcal{B}}(s)}(\varphi_{s}(\Phi(x,0,s),s)+\langle
D_{x}\varphi(\Phi(x,0,s),s),b(\Phi(x,0,s),s)\rangle)d\tilde{m}_{0}(x).$

Now, recalling that
$\\{x\in{\mathbb{R}}^{d}:t_{x}=s\\}=\Psi({\mathcal{S}}^{s},s,s)$ by
definition, we have

$\int_{\\{x\in{\mathbb{R}}^{d}:t_{x}=s\\}}\varphi(\Phi(x,0,s),s)g(s)d\tilde{m}_{0}^{s}(x)\\\
=\int_{\Psi({\mathcal{S}}^{s},s,s)}\varphi(\Phi(x,0,s),s)g(s)d\tilde{m}_{0}^{s}(x)=\int_{{\mathcal{S}}^{s}}\varphi(y,s)d\mu^{s}(s)(y),$

where $\mu^{s}(s):=g(s)(\Phi(\cdot,0,s)\\#\tilde{m}_{0}^{s})$. Finally, we
obtain

$\frac{d}{ds}\int_{{\mathbb{R}}^{d}}\varphi(y,s)d\tilde{\mu}(s)(y)\\\
=\int_{{\mathbb{R}}^{d}\setminus\Phi({\mathcal{B}}(s),0,s)}(\varphi_{s}(y,s)+\langle
D_{x}\varphi(y,s),b(y,s)\rangle)d\mu(s)(y)\\\
-\int_{{\mathcal{S}}^{s}}\varphi(y,s)d\mu^{s}(s)(y).$

Since $\tilde{\mu}(0)=\tilde{m}_{0}$, integrating this between $0$ and $T$ we
get

$\int_{{\mathbb{R}}^{d}}\varphi(y,0)d\tilde{m}_{0}(y)\\\
+\int_{0}^{T}\int_{{\mathbb{R}}^{d}\setminus\Phi({\mathcal{B}}(s),0,s)}(\varphi_{s}(y,s)+\langle
D_{x}\varphi(y,s),b(y,s)\rangle)d\mu(s)(y)ds\\\
-\int_{0}^{T}\int_{{\mathcal{S}}^{s}}\varphi(y,s)d\mu^{s}(s)ds=0,$

that is

$\int_{{\mathbb{R}}^{d}}\varphi(y,0)d\tilde{m}_{0}(y)\\\
+\int_{0}^{T}\int_{{\mathbb{R}}^{d}\setminus\Phi({\mathcal{B}}(s),0,s)}(\varphi_{s}(y,s)+\langle
D_{x}\varphi(y,s),b(y,s)\rangle)d\mu(s)(y)ds\\\
-\int_{0}^{T}\int_{{\mathbb{R}}^{d}}\mathbbm{1}_{\\{({\mathcal{S}}^{s},s):s\in[0,T]\\}}\varphi(y,s)d\mu^{s}(s)(y)ds=0.$

∎

######  Remark 3.

The Disintegration Theorem in the previous proof is applied as follows: we set
$Y={\mathcal{B}}(s)$, $X=[0,s]$ and we consider the map (4.10)

$\pi:Y\longrightarrow X,\qquad\pi(x)=t_{x}$

and $\nu=\pi\\#\tilde{m}_{0}\in{\mathcal{G}}(X)$. In this way
$\pi^{-1}(\tau)=\\{x\in{\mathbb{R}}^{d}:t_{x}=\tau\\}$ for every
$\tau\in[0,s]$. Then, there exists a $\nu$-almost everywhere uniquely
determined family
$\\{\tilde{m}_{0}^{\tau}\\}_{\tau\in[0,s]}\subset{\mathcal{G}}(Y)$ such that
for every Borel measurable function $f:Y\longrightarrow[0,+\infty]$,

$\int_{Y}f(y)d\tilde{m}_{0}(y)=\int_{0}^{s}\int_{\\{x\in{\mathbb{R}}^{d}:t_{x}=\tau\\}}f(y)d\tilde{m}_{0}^{\tau}(y)d\nu(\tau).$

Moreover , in view of (4.11), that is the Lipschitz continuity of the map
$s\longmapsto{\mathcal{L}}^{d}({\mathcal{B}}(s))$, the measure $\nu$ is
absolutely continuous on $X$ with a $L^{\infty}$ density denoted by $g$.

######  Theorem 4.4.

The continuity equation (4.5) has a unique solution given by
$s\longmapsto\tilde{\mu}(s)$.

###### Proof.

Let $\varphi\in{\mathcal{C}}^{\infty}({\mathbb{R}}^{d})$ with
$\operatorname{supp}(\varphi)\subset{\mathbb{R}}^{d}\setminus\Phi({\mathcal{B}}(t),0,t)$
for every $t\leq T$. Fix $t\leq T$ and let us consider the map

$w:{\mathbb{R}}^{d}\times[0,t]\longrightarrow{\mathbb{R}},\qquad
w(x,s):=\varphi(\Phi(x,0,t-s)).$ (4.12)

Then, $w$ is Lipschitz continuous in both variables
$(x,s)\in{\mathbb{R}}^{d}\times[0,t]$ with
$\operatorname{supp}(w)\subset({\mathbb{R}}^{d}\setminus\Phi({\mathcal{B}}(s),0,s))\times[0,t]$.
Moreover, by (4.7) we have

$\varphi(x)=w(\Psi(x,t-s,t),s)=\varphi(\Phi(\Psi(x,t-s,t),0,t-s))$

and, recalling that $\Psi(x,t,t)$ is the solution of (4.6) with $\tau=t$, the
function $w$ satisfies

$0=\frac{d}{ds}\varphi(x)=w_{s}(\Psi(x,t-s,t),s)\\\ +\langle
D_{x}w(\Psi(x,t-s,t),s),b(\Psi(x,t-s,t),s)\rangle\ \ \text{a.e.}$

and hence, in general,

$w_{s}(y,s)+\langle D_{x}w(y,s),b(y,s)\rangle=0\quad\text{a.e. in}\
{\mathbb{R}}^{d}\times(0,t).$

Using $w$ as a test function for a generic $\mu$ satisfying Definition 4.1,
for almost all $s\leq t$ we have

$\frac{d}{ds}\int_{{\mathbb{R}}^{d}}w(y,s)d\mu(s)(y)=\int_{{\mathbb{R}}^{d}}w_{s}(y,s)d\mu(s)(y)+\int_{{\mathbb{R}}^{d}}w(y,s)d\mu_{s}(s)(y)\\\
=\int_{{\mathbb{R}}^{d}}(-\langle D_{x}w(y,s),b(y,s)\rangle+\langle
D_{x}w(y,s),b(y,s)\rangle)d\mu(s)(y)=0$

since
$\operatorname{supp}(w)\subset({\mathbb{R}}^{d}\setminus\Phi({\mathcal{B}}(s),0,s))\times[0,t]$,
which implies

$\int_{{\mathcal{S}}^{s}}w(y,s)d\mu^{s}(s)(y)\\\
=\int_{\Phi({\mathcal{B}}(s),0,s)}(w_{s}(y,s)+\langle
D_{x}w(y,s),b(y,s)\rangle)d\mu(s)(y)=0,\ \ s\leq t.$

Therefore, integrating between $0$ and $t$, we get

$\int_{{\mathbb{R}}^{d}}w(y,t)d\mu(t)(y)=\int_{{\mathbb{R}}^{d}}w(y,0)d\mu(0)(y)$

and then

$\int_{{\mathbb{R}}^{d}}\varphi(y)d\mu(t)(y)=\int_{{\mathbb{R}}^{d}}\varphi(\Phi(y,0,t))d\tilde{m}_{0}(y),$

which shows that $\mu(t)=\Phi(\cdot,0,t)\\#\tilde{m}_{0}$ on
${\mathbb{R}}^{d}\setminus\Phi({\mathcal{B}}(t),0,t)$.

Now we have to prove that $\mu(t)=0$ on $\Phi({\mathcal{B}}(t),0,t)$, that is

$\int_{{\mathbb{R}}^{d}}\varphi(x)d\mu(t)(x)=0$

for any $\varphi\in{\mathcal{C}}^{\infty}({\mathbb{R}}^{d})$ with
$\operatorname{supp}(\varphi)\subset\Phi({\mathcal{B}}(t),0,t)$ for every
$t\leq T$. Again fix $t\leq T$ and let us consider the map (4.12). Then,
proceeding as before, we get

$w_{s}(y,s)+\langle D_{x}w(y,s),b(y,s)\rangle=0\quad\text{a.e. in}\
{\mathbb{R}}^{d}\times(0,t).$

Using $w$ as a test function for a generic $\mu$ satisfying Definition 4.1,
for almost all $s\leq t$ we have

$\frac{d}{ds}\int_{{\mathbb{R}}^{d}}w(y,s)d\mu(s)(y)=\int_{{\mathbb{R}}^{d}}w_{s}(y,s)d\mu(s)(y)+\int_{{\mathbb{R}}^{d}}w(y,s)d\mu_{s}(s)(y)\\\
=\int_{{\mathbb{R}}^{d}}(-\langle D_{x}w(y,s),b(y,s)\rangle+\langle
D_{x}w(y,s),b(y,s)\rangle)d\mu(s)(y)\\\
-\int_{\Phi({\mathcal{B}}(s),0,s)}(w_{s}(y,s)+\langle
D_{x}w(y,s),b(y,s)\rangle)d\mu(s)(y)-\int_{{\mathcal{S}}^{s}}w(y,s)d\mu^{s}(s)(y)\\\
=-\int_{\Phi({\mathcal{B}}(s),0,s)}(w_{s}(y,s)+\langle
D_{x}w(y,s),b(y,s)\rangle)d\mu(s)(y)\\\
-\int_{{\mathcal{S}}^{s}}w(y,s)d\mu^{s}(s)(y).$

Now, observe that

$\int_{{\mathcal{S}}^{s}}w(y,s)d\mu^{s}(s)(y)+\int_{\Phi({\mathcal{B}}(s),0,s)}(w_{s}(y,s)+\langle
D_{x}w(y,s),b(y,s)\rangle)d\mu(s)(y)\\\
=\int_{\\{x\in{\mathbb{R}}^{d}:t_{x}=s\\}}w(\Phi(y,0,s),s)g(s)d\mu^{s}(0)(y)\\\
+\int_{{\mathcal{B}}(s)}(w_{s}(\Phi(y,0,s),s)+\langle
D_{x}w(\Phi(y,0,s),s),b(\Phi(y,0,s),s)\rangle)d\mu(0)(y)\\\
=\frac{d}{ds}\int_{{\mathcal{B}}(s)}w(\Phi(y,0,s),s)d\mu(0)(y).$

Then, by (4.12) and the semigroup property of the flow $\Phi$, since
$\mu(0)=\tilde{m}_{0}$ we get

$\frac{d}{ds}\int_{{\mathbb{R}}^{d}}w(y,s)\mu(s)(y)=-\frac{d}{ds}\int_{{\mathcal{B}}(s)}\varphi(\Phi(y,0,t))d\tilde{m}_{0}(y)$

and hence, integrating between $0$ and $t$,

$\int_{{\mathbb{R}}^{d}}w(y,t)d\mu(t)(y)=\int_{{\mathbb{R}}^{d}}w(y,0)d\mu(0)(y)-\int_{0}^{t}\frac{d}{ds}\int_{{\mathcal{B}}(s)}\varphi(\Phi(y,0,t))d\tilde{m}_{0}(y).$

Therefore

$\int_{{\mathbb{R}}^{d}}\varphi(y)d\mu(t)(y)=\int_{{\mathbb{R}}^{d}}\varphi(\Phi(y,0,t))d\tilde{m}_{0}(y)-\int_{{\mathcal{B}}(t)}\varphi(\Phi(y,0,t))d\tilde{m}_{0}(y)\\\
+\int_{{\mathcal{B}}(0)}\varphi(\Phi(y,0,t))d\tilde{m}_{0}(y).$

Since $\tilde{m}_{0}=0$ in ${\mathcal{S}}={\mathcal{B}}(0)$ and
$\mathrm{supp}\;(\varphi)\subset\Phi({\mathcal{B}}(t),0,t)$, the thesis
follows. ∎

###  4.2. On the measure dependence of the field $b$

We consider now the field $b$ depending on the measure, that is

$\displaystyle b:C^{0}([0,T],{\mathcal{G}})\times{\mathbb{R}}^{d}\times[0,T]$
$\displaystyle\longrightarrow{\mathbb{R}}^{d},$ $\displaystyle(\mu,x,t)$
$\displaystyle\longmapsto b(\mu,x,t).$

We suppose that it is bounded and continuous in the whole entry $(\mu,x,t)$
and Lipschitz continuous w.r.t. $x\in{\mathbb{R}}^{d}$ uniformly w.r.t.
$(\mu,t)\in C^{0}([0,T],{\mathcal{G}})\times[0,T]$, that is, there exists
$L>0$ such that

$\|b(\mu,x,t)-b(\mu,y,t)\|\leq L\|x-y\|,\quad\forall x,y\in{\mathbb{R}}^{d},\
(\mu,t)\in C^{0}([0,T],{\mathcal{G}})\times[0,T].$

In the evolution of the flow given by $b$, the sink is always represented by
${\mathcal{S}}$. In the sequel, for every $\mu\in C^{0}([0,T],{\mathcal{G}})$
fixed, we will also use the notation

$\displaystyle b[\mu]:{\mathbb{R}}^{d}\times[0,T]$
$\displaystyle\longrightarrow{\mathbb{R}}^{d}$ $\displaystyle(x,t)$
$\displaystyle\longmapsto b[\mu](x,t):=b(\mu,x,t)$

and we consider the corresponding flow with sink evolution given by the field
$b[\mu]$. As, for every fixed $\mu$, we denote by $\tilde{\mu}$ the unique
solution of the corresponding problem (4.5), which is, for every $t\in[0,T]$,

$\tilde{\mu}(t)=\begin{cases}\Phi[\mu](\cdot,0,t)\\#\tilde{m}_{0}&\text{on
}{\mathbb{R}}^{d}\setminus\Phi[\mu]({\mathcal{B}}[\mu](t),0,t)\\\
0,&\text{otherwise}\end{cases},$ (4.13)

where $\Phi[\mu]$ is the flow generated by the field $b[\mu]$ and
${\mathcal{B}}[\mu](\cdot)$ is the corresponding sink-reaching-points set. We
also denote by $\pi[\mu]$ the map as in (4.10).

######  Theorem 4.5.

Let us suppose that $\pi[\mu]$ is Lipschitz continuous uniformly in $\mu\in
C^{0}([0,T],{\mathcal{G}})$ (see Remark 4). For every $\mu\in
C^{0}([0,T],{\mathcal{G}})$, we have $\tilde{\mu}\in
C^{0}([0,T],{\mathcal{G}})$. Moreover, the function

$\psi:C^{0}([0,T],{\mathcal{G}})\longrightarrow
C^{0}([0,T],{\mathcal{G}}),\qquad\psi(\mu):=\tilde{\mu}$

has a fixed point in $C^{0}([0,T],{\mathcal{G}})$. This means that the problem
of flow with sink and with field depending on the measure has a solution.

###### Proof.

At first observe that, under the previous hypotheses, by analogous
considerations as in the case with no measure dependence we have
$\tilde{\mu}(t)\in{\mathcal{G}}$ for all $t\in[0,T]$.

Let us prove that $\tilde{\mu}\in C^{0}([0,T],{\mathcal{G}})$. In particular
we have to prove that, whenever $t_{n}\rightarrow t$ in $[0,T]$, then
$\tilde{\mu}(t_{n})\rightarrow\tilde{\mu}(t)$ weakly-star. This comes from
standard regularity results for the push-forward measures (see (4.13)) and
from the fact that ${\mathcal{B}}[\mu](t_{n})\rightarrow{\mathcal{B}}[\mu](t)$
in the Hausdorff metrics and as $d$-dimensional Lebesgue measure, and hence
similarly for $\Phi[\mu]({\mathcal{B}}[\mu](t_{n}),0,t_{n})$ and
$\Phi[\mu]({\mathcal{B}}[\mu](t),0,t)$.

Now, to prove the second statement of the theorem, we have to show that the
function $\psi$ is continuous and compact. In this way, we can conclude by the
Schauder fixed point Theorem.

At first we prove the continuity of $\psi$. Let $\mu_{n}\rightarrow\mu$ in
$C^{0}([0,T],{\mathcal{G}})$. We have to prove that
$\tilde{\mu}_{n}\rightarrow\tilde{\mu}$. Let us consider the two trajectories

$x_{n}(t)=x_{0}+\int_{0}^{t}b[\mu_{n}](x_{n}(s),s)ds,\qquad
x(s)=x_{0}+\int_{0}^{t}b[\mu](x(s),s)ds$

and we prove that $x_{n}$ uniformly converges to $x$ on compact sets. Indeed,
by the continuity and boundedness hypotheses, the field $b$ is bounded and
uniformly continuous on $(\\{\mu_{n}\\}_{n}\cup\\{\mu\\})\times K\times[0,T]$,
where $\\{\mu_{n}\\}_{n}\cup\\{\mu\\}$ is the compact set in
$C^{0}([0,T],{\mathcal{G}})$ given by the whole sequence with its limit, and
$K\subset{\mathbb{R}}^{d}$ is compact. Then, the sequence $b[\mu_{n}]$ is
bounded and equicontinuous on $K\times[0,T]$ and moreover it pointwise
converges to $b[\mu]$. Hence, by the Ascoli-Arzelà Theorem, the convergence is
uniform on $K\times[0,T]$. From this we deduce the desired uniform convergence
of the trajectories. Therefore, also the flows $\Phi[\mu_{n}]$ uniformly
converge to $\Phi[\mu]$ on compact sets. From this we have that
${\mathcal{B}}[\mu_{n}](\cdot)$ converges to ${\mathcal{B}}[\mu](\cdot)$ in
the Hausdorff distance and uniformly in time, and we conclude the convergence
of $\tilde{\mu}_{n}$ to $\tilde{\mu}$ in $C^{0}([0,T],{\mathcal{G}})$, that
is, for every $\varphi\in C_{c}^{0}({\mathbb{R}}^{d})$,

$\sup_{t\in[0,T]}\left|\int_{{\mathbb{R}}^{d}}\varphi(x)d\mu_{n}(t)(x)-\int_{{\mathbb{R}}^{d}}\varphi(x)d\mu(t)(x)\right|\rightarrow
0\quad\text{as}\ n\rightarrow+\infty.$

It remains to prove that $\psi$ has compact image. We can restrict to measures
on a compact set $K\subset{\mathbb{R}}^{d}$ independent of $\mu$, which
contains all the possible compact supports of the measures $\tilde{\mu}$
(because of bounded dynamics and finite horizon). Then, since
$s\longmapsto\tilde{\mu}(s)$ is Lipschitz continuous uniformly in $\mu$
(similarly as in Lemma 4.2 and using the hypothesis on $\pi[\mu]$) with values
in the compact set ${\mathcal{G}}$ of Radon measures on $K$, we get the
desired conclusion by Ascoli-Arzelà Theorem. ∎

######  Remark 4.

In the previous proof, we used the fact that $s\longmapsto\tilde{\mu}(s)$ is
Lipschitz continuous uniformly in $\mu$ and this comes from the Lipschitz
continuity of $\pi[\mu]$ uniformly w.r.t. $\mu$. However, just assuming only
the Lipschitz continuity of $\pi[\mu]$, and not necessarily uniformly in
$\mu$, after some calculations one can prove the equicontinuity w.r.t. $\mu$
of $s\longmapsto\tilde{\mu}(s)$ and then still apply Ascoli-Arzelà Theorem.

###  4.3. A differential characterization of the first arrival time

In this paragraph, we see a characterization of the possible first arrival
time $t_{x}$ (4.2) to the sink ${\mathcal{S}}$ as the unique viscosity
solution of an Hamilton-Jacobi-Bellman equation with proper boundary
conditions. To do this, we define

$\displaystyle\mathrm{T}:{\mathbb{R}}^{d}\times[0,T]\longrightarrow[0,+\infty[,$
$\displaystyle\mathrm{T}(x,t):=t_{(x,t)}=\inf\\{\zeta\geq
t:\Phi(x,t,\zeta)\in{\mathcal{S}}\\}.$

Fixed $s\in[t,T]$, the sink-reaching-points set, that is the set from which
the agents are able to reach the target ${\mathcal{S}}$, is

${\mathcal{B}}(s)=\\{(x,t)\in{\mathbb{R}}^{d}\times[0,T]:\mathrm{T}(x,t)\leq
s\\}.$

Let us set

${\mathcal{B}}:=\bigcup_{s\in[t,T]}{\mathcal{B}}(s)=\\{(x,t)\in{\mathbb{R}}^{d}\times[0,T]:\mathrm{T}(x,t)\leq
T\\}.$

We suppose that $\mathrm{T}$ is continuous in ${\mathcal{B}}$ and that

$\mathrm{T}(x,t)=T\ \ \text{for any}\ (x,t)\in\partial{\mathcal{B}}.$ (4.14)

As with similar calculations as in [9], taking into account that now
$\mathrm{T}$ is time dependent (and then the associated Hamilton-Jacobi-
Bellman equation is evolutive) and that there is no control, we have the
following

######  Proposition 2.

The function $\mathrm{T}$ is the unique viscosity solution of

$\begin{cases}-\mathrm{T}_{t}(x,t)-b(x,t)\cdot
D_{x}\mathrm{T}(x,t)-1=0,&(x,t)\in{\mathcal{B}}\\\
\mathrm{T}(x,t)=t,&(x,t)\in(\partial{\mathcal{S}},t)\\\
\mathrm{T}(x,t)=T,&(x,t)\in\partial{\mathcal{B}}\end{cases}.$

Hence, to conclude, the possible first arrival time $t_{x}$ to the sink
${\mathcal{S}}$ is given by the solution $\mathrm{T}(x,t)$ of the above
problem at the time $t=0$.

######  Remark 5.

Observe that the hypotheses in this section are not too restrictive to be
assumed. Indeed, a priori, $\mathrm{T}$ may be discontinuous in
${\mathcal{B}}$. But, as we noticed in Remark 2 too, in view of the possible
mean field game model, the field $b$ comes as the optimal feedback from a
system like $y^{\prime}={\alpha}$, with total controllability. Then it is not
too strong to assume that $\mathrm{T}$ is continuous in ${\mathcal{B}}$.
Moreover, again in view of the possible mean field model, condition (4.14)
holds only on the boundary which is not viable.

##  5\. Approximation of the continuity equation

We want to build an approximation scheme for the continuity equation (4.5). We
include the possibility of a non linearity, non local dependency of the vector
field $b(\mu,x,t):=b[\mu](x,t)$ on the density $\mu$. In similarly to what
already done for the Hamilton-Jacobi equation, we propose a semi-Lagrangian
scheme to approximate the equation. The advantages of the choice of such
framework are the beneficial stability properties of the scheme without
restrictions on the choice of the parameters and some monotone contained in
the formulation of the scheme. The scheme is an adaptation of the one proposed
in [20, 19]. We also refer the reader to [26], where a similar scheme has been
applied to a non-linear continuity equation modeling a kinetic pedestrian
model.

Starting from equation (4.5) restricted on a finite domain $\Omega$,
$b[\mu]:\Omega\times[0,T]\longrightarrow{\mathbb{R}}$ is a given smooth vector
field, $\mu_{0}:=\tilde{m}_{0}$ a smooth initial datum defined on $\Omega$.
Formally, at time $t\in[0,T]$ the solution of (4.5) is given, implicitly, by
the image of the measure $\mu_{0}dx$ induced by the flow
$x\in\Omega\longmapsto\Phi(x,0,t)$, where, given $0\leq s\leq t\leq T$,
$\Phi(x,s,t)$ denotes the solution of

$\begin{cases}\dot{x}(t)=b[\mu](x,t),&t\in]s,T]\\\ x(s)=x\end{cases}.$

Given the discretization parameters $\Delta x$, $\Delta t$, we construct a
grid on $\Omega\times[0,T]$ as in Section 3.1. We call $M,N$ the number of
elements respectively in the spatial and in the temporal grid.

We discretize (4.5) using its representation formula by means of the flow
$\Phi$. This can be considered an extension of the _characteristic method_ for
these kind of equation (see for details [1]). For $j\in\\{1,\ldots,M\\}$,
$k=0,\ldots,N-1$, we define the discrete characteristics as

$\Phi_{i,k}[\mu]:=x_{i}+\Delta t\,b[\mu](x_{i},t_{k}).$

Let us call $\\{{\beta_{j}}\\}$ the set of base functions, such that
$\beta_{j}(x_{i})=\delta_{i,j}$ (the Kronecker symbol) and
$\sum_{i}\beta_{i}(x)=1$ for each $x\in\overline{\Omega}$. We approximate the
solution $\mu$ of the problem (4.5) by a sequence $\\{\mu_{k}(x_{i})\\}_{k}$,
where each $\mu_{i,k}:=\mu_{k}(x_{i})$ approximates

$\frac{1}{|E_{i}|}\int_{E_{i}}\mu(x_{k},t)dx,$

where $E_{i}=[x_{i}-e_{1}\Delta x/2,x_{i}+e_{1}\Delta
x/2]\times[x_{i}-e_{2}\Delta x/2,x_{i}+e_{2}\Delta x/2]$ with
$e_{1}=(1,0)^{T}$ and $e_{2}=(0,1)^{T}$ and therefore $|E_{i}|=\Delta x^{2}$.
We compute the sequence $\\{\mu_{k}\\}_{k}$ by the following explicit scheme:

$\begin{cases}\mu_{i,k+1}=G(\mu_{k},i,k),&k=0,\ldots,N-1,\ i=1,\ldots,M\\\
\mu_{i,0}=\frac{\int_{E_{i}}\mu_{0}(x)dx}{|E_{i}|},&i=1,\ldots,M,\end{cases},$
(5.1)

where $G$ is defined by

$G(\mu,i,k):=\sum_{j=1}^{M}\mathbbm{1}_{S(t_{k},x_{j})}\mu(x,t)\beta_{i}\left(\Phi_{j,k}\left[\mu\right]\right)\mu_{j}\frac{|E_{j}|}{|E_{i}|}.$

We observe that the function $S(t_{k},x_{j})$ above approximates the action of
the term $\mathbbm{1}_{\\{({\mathcal{S}}^{t},t):t\in[0,T]\\}}\mu(x,t)$ in
(4.5) and it works as a pit canceling the contribution of the density in the
switching point described by the map
$S(t_{k},x_{j})=\\{(\Omega,[0,T])\\}\setminus\\{({\mathcal{S}}^{t_{k}},t_{k})\\}$,
i.e. all the points that are not points of stopping at $(t_{k},x_{j})$.

###  5.1. An optimal stopping test for a crowd

Let us consider the simple scenario of an optimal stopping problem for a
density of agents. In the domain $\Omega=[-1,1]^{2}$ an initial density
distribution of

$\mu_{0}(x)=e^{-8\|x\|^{2}}\quad\hbox{ restricted to }\Omega,$

aims to reach a single target point $\operatorname{\mathcal{T}}=(0,0.6)$. We
use the optimal control vector field of the control problem associated (2.4)
as motion vector field for the continuity equation (4.5) (i.e.
$b(x,t)=-D_{x}V_{p}(x,t)$) and the switching (or stopping) map
${\mathcal{S}}^{t}$ as in (4.3).

For the other operators of the optimal stopping problem, we consider, as
before, quadratic penalization of the control in the running cost and
isotropic dynamics

$\ell(x,a,p,t)=\frac{\|a\|^{2}}{2},\quad f(x,a,p)=a.$

The control is in $A={\mathbb{R}}^{2}$ and the switching cost is set as in
test 3.2 to 5

$C(x,p,p^{\prime})=\sum_{j\in{\mathcal{I}}}\chi_{j}(p,p^{\prime})\|x-{\mathcal{T}}_{j}\|.$

The final time is $T=0.26$ while the boundary condition at time $t=T$ is

$\psi_{p}(x,T)=2\sum_{j\in{\mathcal{I}}}\chi_{j}(p,\bar{p})\|x-{\mathcal{T}}_{j}\|\quad\hbox{
for any }p\in{\mathcal{I}},\;x\in[-1,1]^{2}.$

This means that at any moment an agent can choose to aim to the target
lowering its distance from it (paying the relative running cost) or renounce
to it paying the switching cost and exit from the game.

It is worth to underline that in this case, the coupling between equation
(2.4) and (4.5) is only in one direction. The solution of (2.4) is independent
from the distribution of the density therefore $b(x,t)$ is independent from
$\mu(x,t)$. Hence we can proceed, first of all, to the resolution of (2.4)
with discretization parameters $\Delta x=0.04$ and $\Delta t=0.02$, obtaining
a control map $b(x,t)$ and a switching map ${\mathcal{S}}^{t}$ as reported in
Figure 3.

---
Figure 3. Test 2. Control vector field (left) and switch map (right) of the
optimal stopping problem at $t=0.04,0.24$.

As it was predictable in this easy scenario, we observe that the optimal
solution suggests aiming directly to the target with speed inversely
proportional to the square of the distance from it. Once one reached a region
close enough to the target (approximately a ball of radius 0.2 as it is
deducible from the switching map), it is more convenient to renounce to get
closer to the target and paying the switching cost to end the game. Control
vector field and switching map do not change during the evolution of the
system until its final moment $t\approx T$: here, since it is more convenient
paying the switching cost than the final boundary condition, the switching map
includes all the domain.

---
Figure 4. Test 2. Density evolution at various instants $t=0,0.06,0.16,0.24$.

With this control vector field and switching map, we approximate the
continuity equation (4.5) using the scheme (5.1). On the boundary, we use some
standard homogeneous Neumann conditions as described in [19] to preserve the
total mass of agents.

In Figure 4, we observe the evolution of the density. All agents directly aim
to the target point in $(0,0.6)$ (in light blue) with a higher speed if more
distant. This creates a specific concentration of the density in some areas of
the domain close to where the switching map suggests renouncing to reach the
target paying the switching cost. Clearly, inside this area, no density is
present since it has already left the game. The density still present at
$t=0.24$ will disappear at the next discrete time step since the switching map
includes the whole domain. Therefore at the final step, $t=0.26$, the entire
domain will be empty.

We underline that in this case, since the control problem is independent of
the distribution of mass, no congestion phenomena are considered, and any
player acts as it is alone in the domain.

###  5.2. Two optimal visiting problems for a crowd of agents

We add some complexity to our model. We repeat the previous test adding a
dependency of (2.4) from the distribution of density $\mu(x,t)$. This makes
our model a special kind of mean-field game. We keep all the parameters of
test 2 with the only exception of the final time $T=0.5$, $\Delta x=0.066$,
$\Delta t=0.033$ and the running cost

$\ell(x,a,p,t,\mu)=e^{\sum\limits_{p\in{\mathcal{I}}}\mu(x,t,p)}+\frac{\|a\|^{2}}{2},\quad
f(x,a,p)=a.$

This choice is made to penalize the region of high density where, we suppose,
the circulation is more difficult. All agents contributes to the congestion:
even the ones that have switched already to another state of the system.

Due to the mutual dependency of the two equations (2.4)-(4.5), we have to find
a fixed point between the two equations. The existence of it has been largely
discussed in various situations [30, 28] but never in the hybrid control
framework. Some partial results are described in [13] for an optimal stopping
problem. To derive the fixed point numerically, the literature proposes two
possible options: the use of a global descent method between the two-equation
[2], or a fixed point iteration between the two equations [20, 19]. We opted
for the second one in this test, iterating between the two equations till a
stopping criterium is verified. We use as stop criterium the distance in
discrete infinite norm between two iterations of the density

${\mathcal{E}}(k)=\|\mu^{z}-\mu^{z-1}\|=\max_{k,i}|\mu^{z}(x_{i},t_{k})-\mu^{z-1}(x_{i},t_{k})|.$

As we can see in Figure 5, after few iterations, the criterium
${\mathcal{E}}(k)<\Delta x/2$ is reached.

---
Figure 5. Test 3. Evolution of the stopping criterium in the fixed point
iteration.

We show in Figure 6 the control and switch map that we obtain. Different from
Test 2, in this case, the control field is not constant but changes, trying
actively to avoid the region of congestion. The latter is particularly clear
at time $t=0.26$, where the agents around the center of the domain are
encouraged to split laterally to avoid the congestion between them and the
target. Also, the switching map is constant but moves in the congested areas,
suggesting to pay the switching cost to the agents that would take too much
time to circumnavigate the crowded regions.

---
Figure 6. Test 3. Control vector field (left) and switch map (right) of the
optimal stopping problem at $t=0,0.26,0.4$.

The same phenomenon is shown in Figure 7. In this case, we show the
distribution of the density in state $(0)$ (target not visited) and $(1)$
(target visited). In the state $(1)$, as in all the final states where all
targets are visited, the value function is always identically null, so it is
the control map. For this reason, the density that switches to this state
remains where the switch happened. Initially, at $t=0$, all the density is in
state $(1)$; at $t=0.16$, only the portion of the density very close to the
target does not consider it useful to get closer to the target has switched to
state $(1)$. A part of the rest of the density creates a congested region in
the target’s direction, while the rest tries to avoid such a crowded area
approaching the target laterally. This is evident at time t=0.32 where the
density in $(1)$, i.e., the part of the density that has performed the switch,
is all around the target point (differently from the previous Test 1).
Finally, at $t=T=0.5$, all the density has switched to the final state. The
state $0$ is empty, and state $(1)$ shows the location where all the switches
happened. The density tends to aim to the target point spreading as much as
possible to avoid congestions. A relevant part of the density switches almost
immediately due to the penalization imposed in the congested areas.

---
Figure 7. Test 3. Density evolution at various instants $t=0,0.16,0.32,0.5$.

Finally our last test is performed in the same situation of the single player
Test 1.

Therefore, the target set is composed by the points

${\mathcal{T}}_{j}=0.6\,\left(\cos\left(j\frac{2\pi}{3}\right),\sin\left(j\frac{2\pi}{3}\right)\right),\quad
j=1,\ldots,3.$

The switching cost is again

$C(x,p,p^{\prime})=\sum_{j\in{\mathcal{I}}}\chi_{j}(p,p^{\prime})\|x-{\mathcal{T}}_{j}\|,$

while the final time is $T=0.5$. The discretization parameters are $\Delta
x=0.06$ and $\Delta t=0.03$.

Here the situation is more complex, and we limit our description to an
overview of the behavior of the density reported in Figure 8 and 9. We can see
how the whole density, at time $t=0$, all in state $(0,0,0)$ splits in various
parts aiming to the closer target (Figure 8). At $t=0.08$ some portions of the
density have already switched a first time, appearing in the states $(1,0,0)$,
$(0,1,0)$, $(0,0,1)$, while a part of the density located in the center of the
domain (so then equidistant from the various targets) prefers to switch
directly to the final state renouncing in this way to visit any of the
targets. At this point, each one of these densities will aim (trying to avoid
congestions) to the closest target point not visited yet. At $t=0.25$ (Figure
9) all the states (but the initial state $(0,0,0)$) contain a part of the
density. This means, as we mentioned before, that for some agents, it was more
convenient to renounce to visit one or even all the targets and paying the
relative cost of the switch while in all the intermediate states, a portion of
the density is still actively moving. Finally at time $t=0.5$ the whole
density is in the final state $(1,1,1)$. As it was predictable, most of the
switches happened close to the target points or in the region between them.
All the intermediate states are empty.

---
Figure 8. Test 3. Approximated value functions in the various discrete states
of the system
---
Figure 9. Test 3. Approximated value functions in the various discrete states
of the system

##  6\. Conclusion

The paper discusses both theoretical and applicative aspects of using hybrid
control inside a possible macroscopic multi-agents system for an optimal
visiting problem. The results that are shown here also constitute the first
step toward a general theory of mean-field games in the presence of switches
in the dynamics of the problem. This is an issue of considerable interest for
hybrid-controlled systems but also for impulsive controlled structures. To
complete the theory, various points must be fully developed: first of all, to
consider, besides sinks, the presence of sources in the continuity equation.
Afterwards, a major point would be the introduction of a real coupling between
the continuity equation and the HJ variational inequality. In any case, as
shown in the tests Section §5, we may suppose, due to some promising numerical
evidence, the existence of such an equilibrium. We postpone this study to
future research.

## References

  * [1] Achdou, Y., Camilli, F., Capuzzo Dolcetta I.:  Mean field games: numerical methods for the planning problem.  SIAM J. Control Optim., 50(1) 77–109 (2012)
  * [2] Achdou, Y., Capuzzo Dolcetta, I.:  Mean field games: numerical methods.  SIAM J. Numer. Anal., 48(3) 1136–1162 (2010)
  * [3] Bagagiolo, F., Benetton, M.:  About an optimal visiting problem.  Appl. Math. Optim., 65(1) 31–51 (2012)
  * [4] Bagagiolo, F., Danieli, K.:  Infinite horizon optimal control problems with multiple thermostatic hybrid dynamics.  Nonlinear Anal. Hybrid Syst., 6(2) 824–838 (2012)
  * [5] Bagagiolo, F., Faggian, S., Maggistro, R., Pesenti, R.:  Optimal control of the mean field equilibrium for a pedestrian tourists’ flow model.  Netw. Spat. Econ., (2019)
  * [6] Bagagiolo, F., Maggistro, R., Zoppello, M.:  A hybrid differential game with switching thermostatic–type dynamics and costs.  Minimax Theory Appl., 5(2) 151–180 (2020)
  * [7] Bagagiolo, F., Festa, A., Marzufero, L.:  A hybrid control framework for an optimal visiting problem.  arXiv preprint arXiv:2012.07705, (2020)
  * [8] Bagagiolo, F., Festa, A., Marzufero, L.:  The orienteering problem: a hybrid control formulation.  arXiv preprint arXiv:2012.07712, (2020)
  * [9] Bardi, M., Capuzzo Dolcetta, I.:  Optimal control and viscosity solutions of Hamilton-Jacobi-Bellman equations.  Birkhäuser, 1997
  * [10] Barles, G., Souganidis, P. E.:  Convergence of approximation schemes for fully nonlinear second order equations.  Asymptot. Anal., 4 271–283 (1991)
  * [11] Bellman, R.:  Dynamic programming treatment of the travelling salesman problem.  J. ACM, 9(1) 61–63 (1962)
  * [12] Bensoussan, A., Menaldi, J. L.:  Hybrid control and dynamic programming.  Dynam. Contin. Discrete Impuls. Systems, 3(4) 395–442 (1997)
  * [13] Bertucci, C.:  Optimal stopping in mean field games, an obstacle problem approach.  J. Math. Pures Appl., 120 165–194 (2018)
  * [14] Bertucci, C.:  Fokker-Planck equations of jumping particles and mean field games of impulse control.  Ann. Inst. H. Poincaré Anal. Non Linéaire, 37(5) 1211–1244 (2020)
  * [15] Branicky, M. S., Borkar, V. S., Mitter, S. K.:  A unified framework for hybrid control: Model and optimal control theory.  IEEE Trans. Automat. Contr., 43(1) 31–45 (1998)
  * [16] Camilli, F., Falcone, M.:  An approximation scheme for the optimal control of diffusion processes.  ESAIM-Math. Model. Num., 29(1) 97–122 (1995)
  * [17] Cardaliaguet, P.:  Notes on Mean Field Games.  Lecture notes (from P.-L. Lions’ lectures at Collège de France), 2013
  * [18] Cardaliaguet, P., Graber, P. J.:  Mean field games systems of first order.  ESAIM Control Optim. Calc. Var., 21(3) 690–722 (2015)
  * [19] Carlini, E., Festa, A., Silva, F. J., Wolfram, M.-T.:  A semi-Lagrangian scheme for a modified version of the Hughes’ model for pedestrian flow.  Dyn. Games Appl., 1–23 (2016)
  * [20] Carlini, E., Silva, F. J.:  A semi-Lagrangian scheme for a degenerate second order mean field game system.  Discret. Contin. Dyn. S., 35(9) 4269–4292 (2015)
  * [21] Dharmatti, S., Ramaswamy, M.:  Hybrid control systems and viscosity solutions.  SIAM J. Control Optim., 44(4) 1259–1288 (2005)
  * [22] Evans, L. C., Gariepy, R. F.:  Measure theory and fine properties of functions.  Revised edition, CRC Press, 2015
  * [23] Ferretti, R., Zidani, H.:  Monotone numerical schemes and feedback construction for hybrid control systems.  J. Optim. Theory Appl., 165(2) 507–531 (2015)
  * [24] Festa, A., Guglielmi, R., Hermosilla, R., Picarelli, A., Sahu, S., Sassi, A. Silva, F. J.:  Hamilton-Jacobi-Bellman equations.  Lecture Notes in Mathematics, 2180 127–261 (2017)
  * [25] Festa, A., Göttlich, S.:  A mean field game approach for multi-lane traffic management.  IFAC-PapersOnLine, 51(32) 793–798 (2018)
  * [26] Festa, A., Tosin, A., Wolfram, M.-T.:  Kinetic description of collision avoidance in pedestrian crowds by sidestepping.  Kinet. Relat. Models, 11(3) 491–520 (2018)
  * [27] Gavish, B., Graves, S.-C.:  The travelling salesman problem and related problems.  Massachusetts Institute of Technology, Operations Research Center, 1978.
  * [28] Gomes, D. A., Saúde, J.:  Mean Field Models, A Brief Survey.  Dyn. Games Appl., 1–45 (2013)
  * [29] Huang. M.:  Large-population LQG games involving a major player: the Nash certainty equivalence principle.  SIAM J. Control Optim., 48(5) 3318–3353 (2009)
  * [30] Lasry, J.-M., Lions, P.-L.:  Jeux à champ moyen I. Le cas stationnaire.  C. R. Math. Acad. Sci. Paris, 343 619–625 (2006)
  * [31] Lasry, J.-M., Lions, P.-L.:  Mean field games.  Jpn. J. Math., 2(1) 229–260 (2007)
  * [32] Piccoli, B., Rossi, F.:  Generalized Wasserstein distance and its application to transport equations with source.  Arch. Ration. Mech. Anal., 211(1) 335–358 (2014)
  * [33] Piccoli, B., Rossi, F.:  On properties of the generalized Wasserstein distance.  Arch. Ration. Mech. Anal., 222(3) 1339–1365 (2016)
  * [34] Visintin, A.:  Mathematical models of hysteresis, volume 1.  Elsevier, 2006.