Instantaneous control of interacting particle systems in the mean-field limit

TL;DR

This paper develops an instantaneous control method for large interacting particle systems, using mean-field approximations and adjoint-based optimization to effectively steer particles into desired regions, with demonstrated convergence as particle number grows.

Contribution

It introduces a novel control approach for mean-field particle systems using adjoint optimization, applicable to high-dimensional collective dynamics.

Findings

Optimal controls converge in the mean-field limit

Numerical results validate the control strategy

Effective steering of particles demonstrated

Abstract

Controlling large particle systems in collective dynamics by a few agents is a subject of high practical importance, e.g., in evacuation dynamics. In this paper we study an instantaneous control approach to steer an interacting particle system into a certain spatial region by repulsive forces from a few external agents, which might be interpreted as shepherd dogs leading sheep to their home. We introduce an appropriate mathematical model and the corresponding optimization problem. In particular, we are interested in the interaction of numerous particles, which can be approximated by a mean-field equation. Due to the high-dimensional phase space this will require a tailored optimization strategy. The arising control problems are solved using adjoint information to compute the descent directions. Numerical results on the microscopic and the macroscopic level indicate the convergence of optimal controls and optimal states in the mean-field limit,i.e., for an increasing number of particles.

Tables

Table 1. Table 1 . Illustration of convergence as N → ∞ → 𝑁 N\rightarrow\infty using the norms in ( 26 ) with the respective values of M ​ 100 𝑀 100 M100 as reference. All values are given in %.

ICS1	M25	M50	1000	2000	4000	8000
$\left|J-J_{\text{ref}}\right|/\left|J_{\text{ref}}\right|$	1.2	0.2	0.9	1.1	0.8	0.75
$\parallel u-u_{\text{ref}}\parallel /\parallel u_{\text{ref}}\parallel$	4.21	1.51	1.33	1.33	1.32	1.32
$\left|{\rho }^N-{\rho }_{\text{ref}}\right|/\left|{\rho }_{\text{ref}}\right|$	-	-	0.1	0.1	0.1	0.1
ICS2
$\left|J-J_{\text{ref}}\right|/\left|J_{\text{ref}}\right|$	0.84	0.1	0.1	0.04	0.05	0.05
$\parallel u-u_{\text{ref}}\parallel /\parallel u_{\text{ref}}\parallel$	3.43	1.04	1.29	1.35	1.3	1.32
$\left|{\rho }^N-{\rho }_{\text{ref}}\right|/\left|{\rho }_{\text{ref}}\right|$	-	-	0.01	0.01	0.1	0.01
ICS3
$\left|J-J_{\text{ref}}\right|/\left|J_{\text{ref}}\right|$	0.68	0.12	0.14	0.29	0.13	0.09
$\parallel u-u_{\text{ref}}\parallel /\parallel u_{\text{ref}}\parallel$	2.8	1.01	1.21	1.21	1.22	1.21
$\left|{\rho }^N-{\rho }_{\text{ref}}\right|/\left|{\rho }_{\text{ref}}\right|$	-	-	0.01	0.01	0.01	0.01

Table 2. Table 2 . Illustration of relative computational times and memory consumption for 10 time steps of the optimization algorithm. The reference values are the values of M50, respectively.

computation time	M25	M100	1000	2000	4000	8000
$t/t_{\text{ref}}$	0.15	13	0.16	0.64	2.74	10.75
memory	M25	M100	1000	2000	4000	8000
$m/m_{\text{ref}}$	0.08	15.84	0.15	0.56	2.20	9.26

Equations

$$x_i,v_i,d_m,u_m:I\to {\mathbb{R}}^D,\text{for }i=1,\dots N\text{ and }m=1,\dots ,M.$$

$$\text{x}(t)=(x_i(t){)}_{i=1,\dots ,N}\text{ and }\text{v}(t)=(v_i(t){)}_{i=1,\dots ,N},$$

$$\text{d}(t)=(d_m(t){)}_{m=1,\dots ,M}\text{ and }\text{u}(t)=(u_m(t){)}_{m=1,\dots ,M},$$

$$K_j(x,y)=(\nabla {\Phi }_j)(x-y),j=1,2.$$

$$\frac{\mathrm{d}}{\mathrm{d}t}{x}_i$$

$$\frac{\mathrm{d}}{\mathrm{d}t}{d}_m$$

$$\mathbf{S}(y)=({\mathbf{S}}_i(y){)}_{i=1,\dots ,N},{\mathbf{S}}_i(y)=-\frac{1}{N}\sum _{j\ne i}{K}_1(x_i,x_j)-\frac{1}{M}\sum _{m=1}^M{K}_2(x_i,d_m)-\alpha v_i.$$

$$\frac{\mathrm{d}}{\mathrm{d}t}y=(\frac{\mathrm{d}}{\mathrm{d}t}\text{x},\frac{\mathrm{d}}{\mathrm{d}t}\text{v},\frac{\mathrm{d}}{\mathrm{d}t}\text{d})=(\text{v},\mathbf{S}(y),\text{u})=:F(y,\text{u}),$$

$$\text{x}(0)={\text{x}}_0\in {\mathbb{R}}^{DN},\text{v}(0)={\text{v}}_0\in {\mathbb{R}}^{DN},\text{d}(0)={\text{d}}_0\in {\mathbb{R}}^{DM},$$

$$U:=L^2(I,{\mathbb{R}}^{MD})$$

$$\text{d}(t)={\text{d}}_0+{\int }_0^t\text{u}(s)\mathrm{d}s,$$

$$Y:=H^1(I,{\mathbb{R}}^{ND})\times H^1(I,{\mathbb{R}}^{ND})\times H^1(I,{\mathbb{R}}^{MD}).$$

$${\mu }_t^N(x,v)=\frac{1}{N}\sum _{i=1}^N{\delta }_0(x_i(t)-x)\otimes {\delta }_0(v_i(t)-v).$$

$${\rho }_t(A):={\mu }_t(A\times {\mathbb{R}}^D)={\iint }_{A\times {\mathbb{R}}^D}{f}_t(x,v)\mathrm{d}x\mathrm{d}v,$$

$$\frac{\mathrm{d}}{\mathrm{d}t}⟨{\mu }_t^N,h⟩:=\frac{\mathrm{d}}{\mathrm{d}t}\frac{1}{N}\sum _{i=1}^{N}h(z_i)=\frac{1}{N}\sum _{i=1}^N{\nabla }_{x}h(z_i)\cdot \frac{\mathrm{d}}{\mathrm{d}t}{x}_i+{\nabla }_{v}h(z_i)\cdot \frac{\mathrm{d}}{\mathrm{d}t}{v}_i,$$

$$⟨{\partial }_t{\mu }_t^N,h⟩$$

$$\displaystyle=\left<\mu_{t}^{N},\nabla_{x}h\cdot v-\nabla_{v}h\cdot\left(\big{(}K_{1}\ast\rho_{t}^{N}\big{)}(x)+\frac{1}{M}\sum_{m=1}^{M}K_{2}(x,d_{m})+\alpha v\right)\right>,$$

$$\big{(}K_{1}\ast\rho_{t}^{N}\big{)}(x)=\int_{\mathbb{R}^{2D}}K_{1}(x,\bar{x})\,\mathrm{d}\mu_{t}^{N}(\bar{z}).$$

$$0=⟨{\partial }_t{\mu }_t+v\cdot {\nabla }_x{\mu }_t+{\nabla }_v\cdot (S({\mu }_t){\mu }_t),h⟩,$$

$$S(\mu_{t})(x,v,d)=-\big{(}K_{1}\ast\rho_{t}\big{)}(x)-\frac{1}{M}\sum_{m=1}^{M}K_{2}(x,d_{m})-\alpha v$$

$${\partial }_t{\mu }_t+v\cdot {\nabla }_x{\mu }_t+{\nabla }_v\cdot (S({\mu }_t){\mu }_t)=0\text{in the sense of distributions},$$

$$\int_{0}^{T}\int_{\mathbb{R}^{2D}}\big{(}\partial_{t}h_{t}+v\cdot\nabla_{x}h_{t}+S(\mu_{t})\cdot\nabla_{v}h_{t}\big{)}\,\mathrm{d}\mu_{t}\,\mathrm{d}t+\int_{\mathbb{R}^{2D}}h_{0}\,\mathrm{d}\mu_{0}=0.$$

$$\displaystyle\frac{\,\mathrm{d}}{\,\mathrm{d}t}Z=(\frac{\,\mathrm{d}}{\,\mathrm{d}t}x,\frac{\,\mathrm{d}}{\,\mathrm{d}t}v)=\big{(}v,S(\mu_{t})(Z,d)\big{)},\qquad\mu_{t}=\text{law}(Z(t)),$$

$$\partial_{t}p:=(\partial_{t}f_{t},\frac{\,\mathrm{d}}{\,\mathrm{d}t}\textbf{d})=-\big{(}\nabla_{v}\cdot(S(f_{t})f_{t})+v\cdot\nabla_{x}f_{t},{\bf w}\big{)}=:G(p,{\bf w}).$$

$$f_t{|}_{t=0}=f_0\in L^2({\mathbb{R}}^{2D}),\text{d}(0)={\text{d}}_0\in {\mathbb{R}}^{MD}.$$

$$\mathcal{Y}:=H^1(I,L^2(\Omega ))\times H^1(I,{\mathbb{R}}^{MD}).$$

$${\mathbb{E}}_N(\text{x}(t))=\frac{1}{N}\sum _{i=1}^N{x}_i(t),{\mathbb{V}}_N(\text{x}(t))=\frac{1}{N}\sum _{i=1}^N|x_i(t)-{\mathbb{E}}_N(\text{x}(t)){|}^2,$$

$$J_N(y(\text{u}),\text{u})={\int }_0^T\frac{{\sigma }_1}{4T}|{\mathbb{V}}_N(\text{x}(t))-{\bar{V}}_N{|}^2+\frac{{\sigma }_2}{2T}|{\mathbb{E}}_N(\text{x}(t))-E_{\text{des}}{|}^2+\frac{{\sigma }_3}{2M}\Vert \text{u}(t){\Vert }_{{\mathbb{R}}^{MD}}^2\mathrm{d}t,$$

$$J_{\infty }(p(\mathbf{w}),\mathbf{w})={\int }_0^T\frac{{\sigma }_1}{4T}|{\mathbb{V}}_{\infty }(f_t)-{\bar{V}}_{\infty }{|}^2+\frac{{\sigma }_2}{2T}|{\mathbb{E}}_{\infty }(f_t)-E_{\text{des}}{|}^2+\frac{{\sigma }_3}{2M}\Vert {\text{w}}_t{\Vert }_{{\mathbb{R}}^{MD}}^2\mathrm{d}t,$$

$${\mathbb{E}}_{\infty }(f_t)={\iint }_{{\mathbb{R}}^D\times {\mathbb{R}}^D}x{f}_t\mathrm{d}x\mathrm{d}v,{\mathbb{V}}_{\infty }(f_t)={\iint }_{{\mathbb{R}}^D\times {\mathbb{R}}^D}|x-{\mathbb{E}}_{\infty }(f_t){|}^2{f}_t\mathrm{d}x\mathrm{d}v.$$

$${\mathbb{E}}_{\infty }(f_t^N)={\iint }_{{\mathbb{R}}^D\times {\mathbb{R}}^D}x{f}_t^N\mathrm{d}x\mathrm{d}v=\frac{1}{N}\sum _{i=1}^N{x}_i(t),$$

$${\mathbb{V}}_{\infty }(f_t^N)={\iint }_{{\mathbb{R}}^D\times {\mathbb{R}}^D}|x-{\mathbb{E}}_{\infty }(f_t^N){|}^2{f}_t^N\mathrm{d}x\mathrm{d}v=\frac{1}{N}\sum _{i=1}^N(x_i(t)-{\mathbb{E}}_N(\text{x}(t)){)}^2.$$

$$(y_k^{\ast },{\text{u}}_k^{\ast })=\underset{(y,\text{u})\in Y^k\times {\mathcal{U}}_{ad}^k}{\mathrm{argmin}}{J}_N(y,\text{u})$$

Full text

1. Introduction

In the last decades, the behavior of large particle systems and their mean-field limits were intensively investigated on both the theoretical and computational level [25, 2, 12, 11]. Large groups of individuals like flocks of birds and schools of fish, and their attractive and repulsive interaction were considered, which lead to models of different types of collective behavior such as flocking or milling, and a thorough study of their stability, see, e.g., [36, 20, 3, 32] or [12] for a detailed overview. The models were refined using vision cones, self-propulsion and orientation alignment of neighbors[22, 12, 42]. As the behavior of large particle groups of the same type was well understood, the interest in self-organized systems interacting with few external agents arose. This concept and its kinetic limit was first numerically investigated by Albi and Pareschi in [2]. Using interaction potentials introduced by Cucker-Smale [18] or D’Orsogna et al [22], they showed numerically, that the collective behavior of large groups coincides with the behavior of the kinetic model as the number of individuals tends to infinity. While the models for the interaction of particles became more and more realistic, the complexity was significantly increased. At first, interaction potentials were chosen smooth in order to have well-defined derivatives [18, 19, 22]. Later, it was even shown that special classes of singular interaction potentials allow passing to the mean-field limit [26, 5]. Nowadays, particle games are also employed in the field of global optimization [34, 16].

The Model Predictive Control approach, in particular, instantaneous control was applied successfully to traffic flow problems [27] and the Navier-Stokes Equations [28]. In the present work, we use an instantaneous approach to control the interaction of a huge crowd of individuals with a few external agents. The focus is here on the appropriate mathematical models and the tailored numerical approaches to solve these large scale optimization problems. In particular, we are interested in the behavior of microscopic and macroscopic optimization problems in the corresponding mean-field limit.The control parameters are the velocities of the external agents. The cost functionals are designed in such a way that the external agents aim to lead the crowd to a predefined destination while penalizing their kinetic energy. Exemplarily, one might consider a crowd of sheep guided by dogs.

A similar sparse optimization problem was analytically investigated in [24], where the existence of optimal controls on the microscopic level is shown and the concept of $\Gamma$-convergence is used to perform the mean-field limit in a sparse optimal control setting.

Numerically, the computational effort of the state solutions alone is very challenging due to the high-dimensional phase space. On the one hand, we have the pairwise interactions in the microscopic system and on the other hand the mean-field equations are of Vlasov type, which yields a four dimensional problem in two spatial dimensions. The complexity of the optimization problems is even worse, since the state systems need to be solved several times and a huge amount of data needs to be stored. Hence, there is no hope of solving these problems using black-box optimization approaches. Instead, tailored optimization algorithms are required, which are typically based on derivative information [29].

Recently, many researchers have investigated related problems: in [1] an evacuation problem using invisible agents in a crowd was studied and a similar scenario with visible agents in [39], [37] considers a Fokker-Planck feedback control strategy for crowd motion and [23], where optimal strategies for driving a mobile agent in a guidance by repulsion model are studied. Optimal control problems for transport processes are discussed, e.g., in [6]. An extensive overview on the actual mathematical approaches for behavioral social systems can be found in [4].

Here, we are going to construct appropriate numerical approaches which allow us to verify computationally that the controls of the microscopic problem converge to the optimal control of the kinetic problem as the number of individuals increases. The numerical algorithms are based on first-order derivative information using the adjoint variables for interacting particle system and the mean-field equation, respectively.

The manuscript is organized as follows. First, we describe the microscopic interacting particle system and the corresponding kinetic mean-field equation in Section 2. The cost functionals for the constrained optimization problems are introduced in Section 3, where we also specify the instantaneous control approach we use for the numerical simulations. In Section 4 the associated first-order optimality conditions are derived for both the microscopic and the macroscopic optimization problems. The numerical schemes for the state and adjoint equations are presented in Section 5, where also the instantaneous control algorithm is introduced. Numerical results underlining our approach and confirming the desired behavior in the mean-field limit are discussed in Section 6. Concluding remarks are given in Section 7.

2. Microscopic and Mean-Field Optimal Control Problems

First, we describe the agent-based model and its corresponding mean-field limit, which are going to be the respective state systems for the control problems considered later on.

2.1. Microscopic Model

Let $D\geq 1$ denote the dimension of the spatial and the velocity space. The considered particle system consists of $N\in\mathbb{N}$ particles of the same type and $M$ external agents with $N\gg M$. Let $I=[0,T]$ be the time interval of observation. Then, the particles and agents are represented by state vectors

[TABLE]

The vectors

[TABLE]

denote the positions and velocities of the particles and

[TABLE]

the positions and velocities of the external agents, respectively (compare also [2]).

Remark 2.1*.*

Since the time dependence is clear, we often write x or $x_{i}$ instead of $\textbf{x}(t)$ and $x_{i}(t)$, respectively. Note, that we write the vectors x and v in bold to make clear when denoting the positions and the velocities of the individuals and when we refer to the position and velocity space with variables $x$ and $v$ in the mean-field setting.

The interactions of the individuals and external agents are modeled using potentials $\Phi_{j}$, $j=1,2$ which satisfy the following assumption

• (A)

The potentials $\Phi_{j}\colon\mathbb{R}^{D}\to\mathbb{R}$ are radially symmetric and continuously differentiable, with $\nabla\Phi_{j}$ locally Lipschitz and globally bounded.

Remark 2.2*.*

These rather strict assumptions allow for existence and uniqueness of solutions to the adjoint and state systems. The latter are essential for defining the reduced cost functional needed in the algorithms for numerical investigation.

To simplify the presentation, we denote for any $x,y\in\mathbb{R}^{D}$ the interaction forces

[TABLE]

Since $\Phi_{j}$ is symmetric, we have that $K_{j}(x,y)=-K_{j}(y,x)$ for all $x,y\in\mathbb{R}^{D}$.

Altogether, this leads to the particle system

[TABLE]

In addition to the interaction terms, our model includes (linear) friction represented by the parameter $\alpha>0$. The individuals interact pairwise via $K_{1}$ with each other and with the external agents via $K_{2}$. The latter do not interact among themselves nor are influenced by the others. Note that such interactions can easily be considered by a simple transformation of the optimal control variables analogously to the approach in [8]. For system (1) the velocities of the external agents $u$ are assumed to be given. Later, they serve as controls for the optimization problem.

For notational convenience we further define the state vector $y:=(\textbf{x},\textbf{v},\textbf{d})\in\mathbb{R}^{D(2N+M)}$ and the operators

[TABLE]

Using this notation, the ODE system (1) may be written compactly as

[TABLE]

shortly denoted by $y(0)=y_{0}$.

2.1.1. Well-posedness

The well-posedness of (2) follows from standard theory. Here, we consider the Hilbert space parameters, we use

[TABLE]

as the space of control, which is sufficient to prove the well-posedness of (1a). Indeed, d may be explicitly expressed as

[TABLE]

which shows that d is absolutely continuous on $I$. Consequently, we obtain the existence and uniqueness of a global solution due to the theorem of Picard and Lindelöf and we have $\textbf{d}\in H^{1}(I,\mathbb{R}^{MD})$.

Altogether, we are able to define the corresponding control-to-state operator $\mathcal{G}_{N}\colon U\rightarrow Y$, which maps any control parameter $\textbf{u}\in U$ to the unique solution $y=\mathcal{G}_{N}(\textbf{u})$ of (2) in the state space $Y$ defined by

[TABLE]

Note, that the solution $y\in\mathcal{C}^{1}(I,\mathbb{R}^{ND})\times\mathcal{C}^{1}(I,\mathbb{R}^{ND})\times H^{1}(I,\mathbb{R}^{MD})\subset Y$ is indeed in a subspace of $Y$.

2.2. Mean-Field Model

In order to define the limiting problem for an increasing number of individuals $N$ explicitly, we consider the empirical measure

[TABLE]

By definition, $\mu_{t}^{N}\in\mathcal{P}(\mathbb{R}^{2D})$ is a Borel probability measure that assigns the probability $\mu_{t}^{N}(A)$ of finding particles with states within a Borel measurable set $A\subset\mathbb{R}^{2D}$ in the phase space $\mathbb{R}^{2D}$ at time $t\geq 0$. If a Borel probability measure $\mu_{t}\in\mathcal{P}^{ac}(\mathbb{R}^{2D})$ is absolutely continuous w.r.t. the Lebesgue measure, we denote its density by $f_{t}\in L^{1}(\mathbb{R}^{2D})$. For later use we introduce the macroscopic density $\rho_{t}$ of a Borel probability measure $\mu_{t}\in\mathcal{P}(\mathbb{R}^{2D})$ as its first marginal, i.e.,

[TABLE]

for any Borel measurable set $A\subset\mathbb{R}^{D}$. The last equality holds whenever $\mu_{t}\in\mathcal{P}^{ac}(\mathbb{R}^{2D})$.

The link between the particle system and the mean-field equation is derived formally using ideas from [21, 7, 33]. Let $h\in\mathcal{C}^{\infty}_{c}(\mathbb{R}^{2D})$ be an arbitrary smooth function on $\mathbb{R}^{2D}$ with compact support and $z_{i}=(x_{i},v_{i})\in\mathbb{R}^{2D}$. Then, it holds

[TABLE]

which allows for the formal calculation

[TABLE]

where

[TABLE]

Passing to the limit $N\to\infty$ and integrating by parts, we arrive at the equation

[TABLE]

where we define

[TABLE]

analogous to the particle case. Since $h\in\mathcal{C}^{\infty}_{c}(\mathbb{R}^{2D})$ is arbitrary, we may use the variational lemma to find that

[TABLE]

which is the mean-field single particle distribution, supplemented with the initial condition $\mu(0)=\mu_{0}$.

Remark 2.3*.*

Observe that $\mu_{t}^{N}$ and $\mu_{t}$ satisfy exactly the same equation in the distributional sense. The above limit exists for example in the Wasserstein metric (for a detailed discussion see [10]). Note, that even though $\mu_{t}$ has unit mass by construction, we do not need to explicitly assume that the condition is satisfied, since it is a naturally fulfilled. In particular, the control parameter does not change this property.

Here, we only perform the mean-field limit $N\rightarrow\infty$, while the number of external agents $M$ remains finite.

2.2.1. Well-posedness

Results on the existence and uniqueness of solutions for the Vlasov equation (6) can be found, e.g., in [7, 10, 21, 25], where the notion of solution is established in the Wasserstein space of Borel probability measures.

Definition 2.4**.**

Let $\mathcal{P}_{1}(\mathbb{R}^{2D})$ denote the space of Borel probability measures on $\mathbb{R}^{2D}$ with finite first moment. We say that $\mu\in\mathcal{C}(I,\mathcal{P}_{1}(\mathbb{R}^{2D}))$ is a weak measure solution of (6) with initial condition $\mu_{0}\in\mathcal{P}_{1}(\mathbb{R}^{2D})$ if for any test function $h\in\mathcal{C}_{c}^{\infty}((-\infty,T]\times\mathbb{R}^{2D})$ we have

[TABLE]

Remark 2.5*.*

Equation (6) may be equivalently expressed as a nonlinear flow

[TABLE]

with the initial condition $Z(0)=Z_{0}\in\mathbb{R}^{2D}$, where $Z_{0}$ is a random variable distributed according to $\mu_{0}=\text{law}(Z_{0})$. Assuming the solvability of the nonlinear process (7), a weak measure solution of $\eqref{longPDE}$ may be represented as the push-forward of the measure along the flow $Z(t,Z_{0})$, i.e., $\mu_{t}=Z(t,\cdot)\#\mu_{0}$. If we further assume that the initial measure $\mu_{0}\in\mathcal{P}_{1}(\mathbb{R}^{2D})$ has density $f_{0}\in L^{1}(\mathbb{R}^{2D})$, and that the nonlinear flow given by (7) satisfies $Z\in\mathcal{C}(I,\text{Diff}(\mathbb{R}^{2D}))$, i.e., $Z(t,\cdot)$ is a diffeomorphism for all $t\geq 0$, then $\mu_{t}\in\mathcal{P}_{1}^{ac}(\mathbb{R}^{2D})$ has density $f_{t}\in L^{1}(\mathbb{R}^{2D})$ and $f\in\mathcal{C}(I,L^{1}(\mathbb{R}^{2D}))$.

In order to employ the standard $L_{2}$-calculus we will require more regularity of the states and assume additionally that

• (B)

$\mu_{0}\in\mathcal{P}_{1}(\mathbb{R}^{2D})$ has density $f_{0}\in L^{2}(\mathbb{R}^{2D})$ with compact support and the flow given by (7) satisfies $Z\in\mathcal{C}^{1}(I,\text{Diff}(\mathbb{R}^{2D}))$ such that $f\in H^{1}(I,L^{2}(\mathbb{R}^{2D}))$.

Remark 2.6*.*

Note, that the assumption on the nonlinear flow $Z\in\mathcal{C}^{1}(I,\text{Diff}(\mathbb{R}^{2D}))$ satisfying (7) in (B) is not restrictive due to assumption (A). Indeed, $\mu\in\mathcal{C}(I,\mathcal{P}_{1}(\mathbb{R}^{2D}))$ implies that $S(\mu_{t})$ is continuous in $t$ and locally Lipschitz in $z$. Therefore, standard ODE theory provides the required regularity for the nonlinear flow.

Defining the state variable $p:=(f,\textbf{d})$, the coupled system of (6) and the ODE (1b) that models the movement of the external agents d with control ${\bf w}$ may be written as

[TABLE]

In the following we shall refer to (8b) as $p(0)=p_{0}$.

Remark 2.7*.*

The problem is defined on the whole space $\mathbb{R}^{2D}$, such that no boundary conditions are required. In fact, we expect from (A) and (B) that $f_{t}$ has a compact support which remains within a bounded convex domain $\Omega\subset\mathbb{R}^{2D}$ with smooth boundary for all times $t\in I$ (cf. [24, Thm 4.4]).

Finally, we define the corresponding control-to-state operator $\mathcal{G}_{\infty}\colon U\to\mathcal{Y}$, $\textbf{u}\mapsto p=\mathcal{G}_{\infty}(\textbf{u})$, where $p$ satisfies (8) and for the state space we choose

[TABLE]

Here, $\Omega\subset\mathbb{R}^{2D}$ is an a priori given bounded convex domain with smooth boundary which contains the support of $f_{t}$ for all times $t\in I$.

3. Control Problems and Instantaneous Control Approach

Based on the previously introduced state systems we are now in the position to define the respective control problems. The physical task at hand is to guide the crowd of individuals such that it clusters around a certain position, while keeping its variance restricted and utilizing only limited energy for the controlling agents.

Since the question arises on the microscopic and macroscopic level, we need to have complying cost functionals modeling this issue. In fact, they need to be linked in the asymptotic limit.

In the following, we state one possible cost functional using the center of mass and the variance of the crowd, which meets these requirements. Further, we are describing the instantaneous control procedure that we are going to apply to control the crowd.

3.1. Cost Functional in the Microscopic Setting

Again, let x denote the particle positions. Then, the center of mass $\mathbb{E}_{N}$ and variance $\mathbb{V}_{N}$ for the particles are given by

[TABLE]

and we define the cost functional

[TABLE]

where $\bar{V}_{N}>0$ is a given desired variance value.

The first part of the functional penalizes the spread of the particles from the desired variance $\bar{V}_{N}$, while the second term measures the distance of the center of mass from the desired destination $E_{\text{des}}\in\mathbb{R}^{D}$. The third term measures the control costs in terms of the kinetic energy of the agents. Mathematically, it introduces more convexity in the cost functional. The positive parameters $\sigma_{i}$ allow to adjust the influence of the different parts of the cost functional.

3.2. Cost Functional in the Mean-Field Setting

As already discussed, the control problem on the microscopic level has to match the one on the mean-field level as $N\rightarrow\infty$. Thus, we choose the following cost functional

[TABLE]

in which the center of mass $\mathbb{E}_{\infty}$ as well as the variance $\mathbb{V}_{\infty}$ only depend on the macroscopic density and are defined as

[TABLE]

Note that the microscopic cost functional can be derived from the mean-field cost functional by using the empirical measure $\mu^{N}$. In fact, it holds

[TABLE]

3.3. Instantaneous Control

The idea of the instantaneous control approach is to split the time interval of interest into several smaller time intervals and to solve the control problem sequentially on these subintervals. Clearly this reduces the memory consumption of the algorithm significantly. Moreover, it is in our setting more realistic than optimal control over the whole time horizon, because we expect the agents to forecast the movement of the crowd only on short time intervals.

Hence, we split time interval of investigation into $K+1$ subintervals $I_{k}=[t_{k-1},t_{k}]$, where $t_{0}=0$ and $t_{K}=T$, $k=1,\dots,K$ such that $[0,T]=\bigcup_{k=1}^{K}[t_{k-1},t_{k}]$. The corresponding spaces $Y^{k},\mathcal{Y}^{k}$ and $\mathcal{U}_{ad}^{k}$ are precisely defined below. Then, we study the family of optimization problems on each time slice $I_{k}$ on the microscopic level given by

Problem 1**:**

For each subinterval $I_{k}$, $k=1,\dots,K$, find an optimal pair $(y^{*}_{k},\textbf{u}^{*}_{k})\in Y\times\mathcal{U}_{ad}$ such that

$\displaystyle(y_{k}^{*},\textbf{u}_{k}^{*})=\operatorname*{argmin}\nolimits_{(y,\textbf{u})\in Y^{k}\times\mathcal{U}_{ad}^{k}}J_{N}(y,\textbf{u})$

subject to

$\displaystyle\text{ IVP }\eqref{ODEsys}\text{ on }I_{k}$

and initial data $y_{k}(t_{k-1})=y_{k-1}^{*}(t_{k-1})$, where

$\displaystyle J_{N}(y,\textbf{u})=\int_{I_{k}}\frac{\sigma_{1}}{4T}|\mathbb{V}_{N}(\textbf{x}(t))-\bar{V}_{N}|^{2}+\frac{\sigma_{2}}{2T}|\mathbb{E}_{N}(\textbf{x}(t))-E_{\text{des}}|^{2}+\frac{\sigma_{3}}{2M}\lVert\textbf{u}(t)\rVert_{\mathbb{R}^{MD}}^{2}\,\mathrm{d}t.$

On the mean-field level we obtain the following optimization problem:

Problem 2**:**

For each subinterval $I_{k}$, $k=1,\dots,K$, find an optimal pair $(p_{k}^{*},{\bf w}_{k}^{*})\in\mathcal{Y}^{k}\times\mathcal{U}_{ad}^{k}$ such that

$\displaystyle(p_{k}^{*},{\bf w}_{k}^{*})=\operatorname*{argmin}\nolimits_{(p,{\bf w})\in\mathcal{Y}^{k}\times\mathcal{U}_{ad}^{k}}J_{\infty}(p,{\bf w})$

subject to

\displaystyle\text{system \eqref{PDEsys} on }I_{k}\text{ and initial data p_{k}(t_{k-1})=p_{k-1}^{*}(t_{k-1}), },\text{where}

$\displaystyle J_{\infty}(p,{\bf w})=\int_{I_{k}}\frac{\sigma_{1}}{4T}|\mathbb{V}_{\infty}(f_{t})-\bar{V}_{\infty}|^{2}+\frac{\sigma_{2}}{2T}|\mathbb{E}_{\infty}(f_{t})-E_{\text{des}}|^{2}+\frac{\sigma_{3}}{2M}\lVert\textbf{w}_{t}\rVert_{\mathbb{R}^{MD}}^{2}\,\mathrm{d}t.$

Remark 3.1*.*

For the first iterate the initial data is given by $y_{1}^{*}(0)=y_{0}$ and $p_{1}^{*}(0)=p_{0}$, respectively.

Further, we restrict the velocities of the external agents and define the spaces of admissible controls $\mathcal{U}_{ad}^{k}$, for $k=1,\dots,K$, as

[TABLE]

3.3.1. Existence of controls

An existence result for the optimization problem on each time slice may be deduced in a straight-forward way from the results in [24]. Even though the results are based on a sparse control setting, similar arguments may be applied in the present setting as well.

Theorem 3.2**.**

Assume (A) and (B). Then, there exists on each time interval $I_{K}$, $k=1,\ldots,K$, an optimal pair $(y^{*}_{k},\textbf{u}^{*}_{k})\in Y\times\mathcal{U}_{ad}$ for Problem 1 and an optimal pair $(p_{k}^{*},{\bf w}_{k}^{*})\in\mathcal{Y}^{k}\times\mathcal{U}_{ad}^{k}$ for Problem 2, respectively.

Remark 3.3*.*

The choice of initial conditions for each subinterval $I_{k}$ indicates that we obtain a sub-optimal solution for the problem on the interval $[0,T]$ by cluing the optimal solutions of the subintervals together.

4. First Order Necessary Conditions

We are going to apply adjoint based descent methods to solve the control problems on each time slice. Therefore, we formally derive the adjoints and the optimality conditions for the particle and the mean-field optimal control problem with the help of the extended Lagrange functional. Furthermore, we introduce the reduced cost functional and its gradient. Throughout this section we consider an arbitrary time slice $I$ and denote its right-end point $I_{R}.$

4.1. Adjoint of the microscopic problem

Let the control space $U$ and state space $Y$ be the Hilbert spaces

[TABLE]

with $\mathcal{U}_{ad}\subset U$ defined in (12). Further, we define $X:=[L^{2}(I,\mathbb{R}^{ND})]^{2}\times L^{2}(I,\mathbb{R}^{MD})$ and

[TABLE]

as the space of Lagrange multipliers.

The state operator $e_{N}\colon Y\times U\rightarrow Z^{*}$ of the microscopic problem is given by

[TABLE]

or in weak form

[TABLE]

Let $(\xi,\eta)\in Z$ denote the Lagrange multiplier, which is in fact the adjoint state. Then, the extended Lagrange functional corresponding to Problem 1 reads

[TABLE]

As usual the first-order optimality condition of Problem 1 is derived by solving

[TABLE]

The derivative w.r.t. the adjoint state results in the state equation, while the derivative with respect to the state $y$ yields the adjoint system and the optimality condition is obtained by the derivative w.r.t. the control (for details see, e.g., u [29]).

For the calculations we denote the three parts of the cost functional by $J^{i}$, $i=1,2,3$, with

[TABLE]

For any $h=(h^{y}=(h^{x},h^{v},h^{d}),h^{u})\in Y\times U$, the following Gâteaux derivatives are obtained

[TABLE]

For the second part of the Lagrangian we derive

[TABLE]

Assuming that $\xi\in Y$ one may formally derive the strong formulation of the adjoint system. Indeed, integrating by parts and using the fact that $d_{\textbf{x}}{\bf S}(y)$ and $d_{\textbf{d}}{\bf S}(y)$ are symmetric matrices, we arrive at the following result.

Proposition 4.1**.**

Let $(y_{*},u_{*})$ be an optimal pair for Problem 1. Then, the first-order optimality condition corresponding to Problem 1 reads

[TABLE]

where $\xi=(\xi_{1},\xi_{2},\xi_{3})\in Y$ satisfies the adjoint system given by

[TABLE]

4.2. Adjoint of the Mean-Field Problem

Here, we assume that $p=(f,\textbf{d})$ lies within the state space $\mathcal{Y}$ of the PDE optimization problem, where

[TABLE]

Let $\mathcal{X}:=H^{1}(I,L^{2}(\Omega))\cap L^{2}(I,H^{1}(\Omega))\times L^{2}(I,\mathbb{R}^{MD})$ and set $\mathcal{Z}:=\mathcal{X}\times\left(L^{2}(\Omega)\times\mathbb{R}^{MD}\right)$ to be the space of adjoint states with dual $\mathcal{Z}^{*}$. Note that the control space $U$ is identical to the one of the particle problem.

Now, we define the mapping $e_{\infty}\colon\mathcal{Y}\times U\rightarrow\mathcal{Z}^{*}$ by

[TABLE]

with the adjoint state $(\varphi,\eta)\in\mathcal{X}\times\left(L^{2}(\Omega)\times\mathbb{R}^{MD}\right)$, $\varphi=(g,\varphi_{d})$ and $\eta=(\eta_{f},\eta_{d})$.

Similar to the microscopic case we define the extended Lagrangian corresponding to Problem 2 as

[TABLE]

Again, we introduce the three parts of the cost functional by $J^{i}$, $i=1,2,3$, via

[TABLE]

Analogously to the microscopic case we derive the adjoint system and the optimality condition by calculating the derivatives of $\mathcal{L}_{\infty}$ w.r.t. the state variable $p$ in direction $h^{p}=(h^{f},h^{d})\in\mathcal{Y}$, and the control $w$ in direction $h^{w}\in U$. The standard $L_{2}$-calculus yields

[TABLE]

Let $\varphi=(g,\varphi_{d})$ be the adjoint state corresponding to $p=(f,\textbf{d})$. Then, we obtain for the constraint part in the extended Lagrange functional the following Gâteaux derivatives:

[TABLE]

Assuming again that the adjoint state $\varphi_{d}$ is sufficiently regular, we may integrate by parts to obtain a strong formulation of the adjoint system. For the terms involving derivatives of $S(f_{t})$ we calculate the following representations:

[TABLE]

This yields the following adjoint system and optimality condition.

Proposition 4.2**.**

Let $(p_{*},w_{*})$ be an optimal pair for Problem 2. Then, the optimality condition corresponding to Problem 2 reads

[TABLE]

where $\varphi=(g,\varphi_{d})\in\mathcal{Y}$ satisfies the adjoint system given by

[TABLE]

Remark 4.1*.*

Note, that the optimality conditions of Problem 1 and 2 coincide. This is due to the fact that $J^{3}$ and the set of admissible controls are identical in both problems. The mean-field limit only affects the adjoint system, but formally they can be identified in the following way: along the characteristics of the partial differential equation for $g$, $\xi_{1}$ corresponds to $\frac{1}{N}\nabla_{x}g$ and $\xi_{2}$ to $\frac{1}{N}\nabla_{v}g$. Finally, $\varphi_{d}$ can be identified directly with $\xi_{3}$. [9]

Remark 4.2*.*

The variational inequalities (13) and (16) may be equivalently expressed as fixed point problems in terms of a projection operator $\text{Proj}_{U}\colon U\to\mathcal{U}_{ad}$ which is defined by (see [29])

[TABLE]

Consequently, the variational inequalities (13) and (16) may be expressed as

[TABLE]

where $k(\textbf{u})=\sigma_{3}/M\textbf{u}_{*}-\xi_{3}$ for the microscopic case (13) and $k({\bf w})=\sigma_{3}/M{\bf w}_{*}-\varphi_{d}$ for the mean-field case (16).

In our particular case, $\text{Proj}_{U}$ has the explicit representation given by

[TABLE]

Due to the high dimensionality of the problem, it is not recommended solving the system consisting of state equation, adjoint equation and optimality condition all at once. Therefore, we employ a gradient descent method. In order to define this method we employ the gradient of the reduced cost functional which is derived in the following for both the microscopic problem and its mean-field counterpart.

4.3. Gradient of the Reduced Cost Functional

In this section we introduce the reduced cost functionals $\hat{J}_{N}(\textbf{u})$ and $\hat{J}_{\infty}({\bf w})$ and formally calculate their gradients $\nabla\hat{J}_{N}(\textbf{u})$ and $\nabla\hat{J}_{\infty}({\bf w})$ which we need for the descent algorithms.

By means of the control-to-state operators $\mathcal{G}_{N}\colon U\to Y$ and $\mathcal{G}_{\infty}\colon U\to\mathcal{Y}$ introduced in Section 2, we define the reduced cost functionals as

[TABLE]

Assuming sufficient regularity for $\mathcal{G}^{N}$ and $\mathcal{G}^{\infty}$ we further derive the gradients of the reduced cost functionals. Making use of the state equation $e_{N}(y,\textbf{u})=0$ and $e_{\infty}(p,{\bf w})=0$ we implicitly obtain $d\mathcal{G}_{N}(\textbf{u})$ and $d\mathcal{G}_{\infty}({\bf w})$ via

[TABLE]

With the help of the adjoint equations

[TABLE]

we may calculate the Gâteaux derivative of $\hat{J}$ in the direction $h\in U$, which gives

[TABLE]

Since $U$ is a Hilbert space, we may use the Riesz representation theorem to identify the gradients

[TABLE]

5. Numerical Schemes and Instantaneous Control

For the numerical solution of the two control problems we use the instantaneous control algorithm. In this case the time slices are equally sized and the size coincides with the time step size of the numerical scheme, i.e., the gradient is based on the information of one time step only. Further, we discuss the numerical schemes for the solution of the respective forward and backward differential equations.

5.1. Numerics for the Forward and Adjoint IVP

We compute the state variable with one step of the explicit fourth order Runge–Kutta solver for the IVP (2). For the mean-field limit we need a numerical scheme which is independent of the number of particles involved. Hence, we rescale the adjoint ODE by multiplying with $N$. This has the effect that the $N$-dependence of the terms in (14b) emerging from the cost functional vanishes.

In fact, for $i=1,\dots,N$, we set $r_{i}(t)=N\xi_{i}^{1}(t)$ and $s_{i}(t)=N\xi_{i}^{2}(t)$, and obtain the rescaled adjoint ODE system

[TABLE]

We apply one step of the explicit fourth order Runge-Kutta scheme to (20) for computing the adjoint.

Note, that the main computational effort for a large particle number $N$ comes from the interaction potentials.

5.2. Numerics for the Mean-field Equation and its Adjoint

The forward and backward solves for the mean-field optimization are realized using a Strang splitting scheme[17]. This scheme applies a semi-Lagrangian solver in the spatial direction and a semi-implicit finite-volume scheme in the velocity space. Using the following short hand notation for (8)

[TABLE]

we define the splitting

[TABLE]

A Semi-Lagrangian method [38, 30] is used to solve (21b), that means the computations rely on a fixed grid and yet we make use of the Lagrangian ansatz and consider the transport along characteristics. To obtain the characteristic curves we solve ODEs using a second order Runge-Kutta scheme. The inital point for each transport step is a grid point. Since we cannot ensure that also the endpoint of each transport step is again grid point, we need to interpolate the data to the grid. This interpolation is realized by a polynomial reconstruction based on cubic Bezier curves, which is again of second order.

For the discretization in velocity space, (21a) and (21c), we adopt a second order finite volume scheme where the advection is approximated by a Lax-Wendroff flux [31, 35]. Oscillations caused by non-smooth solutions are intercepted with the help of a van-Leer limiter [41]. More details on this second order scheme can be found in [14].

Basically the same code is used for the adjoint system: we rewrite

[TABLE]

Further, we need to add the term resulting from the linearization of the non-linear interaction and the cost functional. Altogether we obtain the splitting for the adjoint system

[TABLE]

with $D_{f}(f_{t})[g_{t}]$ as defined in (15) and $d_{x}J_{\infty}$ as in (17c).

Remark 5.1*.*

Note, that the forward as well as the backward system are high-dimensional, since for our example in $\mathbb{R}^{2}$ we need in fact to solve equations in $\mathbb{R}^{4}$. This yields an immense computational effort and is also challenging for the data handling. The last is the essential bottleneck for an optimal control approach, since then one needs to store the whole time evolution of the forward problem to do one backward solve for the evaluation of the gradient. Nevertheless, numerical results in the optimal control setting on a short time horizon can be found in [40].

5.3. Instantaneous Control Algorithm

We use the steepest descent steps to update the control once on every time slice, i.e.,

[TABLE]

where $c_{k}$ denotes the current control and $q_{k}=\nabla\hat{J}(u_{k})$. Due to the uniform bound on the velocity we need to ensure that the computed the controls are feasible. Therefore, we project the controls onto the feasible set in every iteration using the operator $\text{Proj}_{U}$ defined in (18). The step sizes $\omega_{k}$ are obtained with help of the projected Armijo step size rule [29] (see Algorithm 1).

The control on time slice $k+1$ is initialized with

[TABLE]

which ensures that the cost term of the cost functional is initially non-zero. Hence, changing the direction of the controls while fixing the speed does not increase the third part of the cost functional. Therefore, it is easier to find descent directions.

We combine the solvers for the state system, the adjoint system and the gradient update to obtain the Instantaneous Control algorithm, see Algorithm 2.

Remark 5.2*.*

Note that we perform only one gradient step on each time slice. Numerical tests have shown that more gradient steps would only marginally improve the results, but drastically increase the computational effort. Further, we use the time slices of the instantaneous approach as time discretization for the simulation. There are no intermediate steps computed on one time slice.

6. Numerical Results

The numerical simulations for the mean-field equation are performed on the computational domain $\Omega=[-100,100]^{2}\times[-5,5]^{2}\subset\mathbb{R}^{4}$, i.e., we set $D=2$. Further, we use the scaled CFL condition

[TABLE]

with $L=200,|V|=5,dt=0.02$ and thus $K=T/dt=500.$

The grid parameter $h$ is varied throughout the simulations to investigate the grid convergence of the scheme. In fact, we use $25,50$ or $100$ grid points in each of the two directions leading to $h=0.04,0.02$ or $0.01$.

Our particular choice of the interaction potentials are the Morse potentials as proposed in [22, 15]. For fixed positive parameters $A_{j},a_{j},R_{j},r_{j}$ we have

[TABLE]

The parameters $A_{j},R_{j}$ denote the attraction and repulsion strengths and $a_{j},r_{j}$ the radii of interaction. The case $j=1$ refers to the interaction of the individuals, the interaction of individuals with external agents is denoted by $j=2$. Inspired by [13], we use the values

[TABLE]

Remark 6.1*.*

Note, that the coefficients of the interaction potentials are scaled with respect to the spatial size of the domain $\Omega$. As mentioned in the introduction we had the model problem of dogs herding sheep to the destination $E_{\text{des}}$ in mind. The parameters for the sheep-sheep interaction have long-range attraction and repulsion on a very short range. The parameters modeling the sheep-dog interaction have a larger repulsive influence in order to reflect the guiding property correctly.

Further parameters are fixed as follows: the destination $E_{\text{des}}=(-20,-20)$ is depicted in green. The initial spatial support of the sheep is set to $\Omega_{0}=[-10,55]\times[-20,55]$, i.e., $f_{0}(x,v)$ is the uniform distribution and the initial positions and velocities of the sheep are chosen by realizations of i.i.d. random variables with $\text{law}(f_{0})$. In every figure, the sheep and dogs are represented by blue markers and red triangles, respectively. The initial configurations of the microscopic and the mean-field case are shown in Figure 1.

The parameters for the following calculations are set to

[TABLE]

Thus, the desired variance is $10\%$ less than the initial variance given by the initial distribution of the crowd. The Armijo parameter is set to $\omega_{0}=1\cdot 10^{3}.$

Remark 6.2*.*

Note, that in the given time interval $[0,T]$ the task of steering the crowd to the destination $E_{\text{des}}$ is impossible to realize. In fact, our focus lies on the comparison of the behavior of the controls as $N$ increases.

6.1. Numerical Results Using the IC Algorithm

In this section we discuss the results obtained by the IC algorithm. First, we illustrate the influence of the cost functional parameters. Therefore, we set up three test cases: S1 stresses the variance part $J_{1}$ of the cost functional, we expect the dogs to move towards the corners of the crowd in order to reduce the spread of the crowd. In S2 the second part $J_{2}$ is emphasized. Thus, the dogs would be more inclined to push the crowd towards the destination. For test case S3 we choose the weights of the cost functional such that the focus lies on steering the crowd to the destination $E_{\text{des}}$ while the variance term has minor influence but cannot be neglected. In particular, we use the following parameter values

[TABLE]

These choices assure that $J_{1}$ of test case S1 has the same order of magnitude as $J_{2}$ of test case S2 and the other way around. The following notation is used in the Figures:

[TABLE]

The spatial and the velocity grid are discretized with the same number of grid points ($25,50$ or $100$) in each of the two directions. The graphs corresponding to mean-field solutions are labeled M$\#$ where $\#$ denotes the number of grid points. The graphs corresponding to microscopic simulations are denoted by their respective number of particles. In Figure 2 the influence of the cost functional parameters is illustrated.

Considering the application we expect that S3 is the most realistic setting, since the focus lies on $J_{2}$ which measures the distance to the destination, while the influence of $J_{1}$ is not negligible. In Figure 3(left) the trajectories of the dogs and the crowd at $T=10$ are shown. Note, that the dogs are not splitting the crowd as much as in test case S1.

Comparing the values of $J$, $J_{1}$ and $J_{2}$ in Figure 3(right) and Figure 4 we conclude that the part measuring the crowd’s distance to the destination is dominating the others.

Remark 6.3*.*

Note, that the graphs of the cost functional values are given with respect to time, i.e., it may happen that the values increase. Since the initial positions of the dogs is chosen arbitrarily, they move to appropriate positions first, causing a slight increase of the cost functional. Afterwards the cost functional is decreasing as expected.

The simulation on the coarse mean-field discretization $M25$ overestimates the variance term significantly as can be seen in Figure 4(left). The evolution of $J_{1}$ is in good resemblance for all other discretizations. Similarly, there is a good agreement in the graphs showing the evolution of $J_{2}$. Furthermore, we see the convergence of the graphs to the solution for $M100$ as $N$ increases and the convergence of the mean-field simulations as the grid is refined.

In Table 1 we illustrate the convergence for $N\rightarrow\infty$ by the relative errors of

[TABLE]

The reference values are the results of the simulation $M100$. To compute the norms the empirical density $\rho^{N}$ of the microscopic simulations is approximated by a histogram which is based on the grid of the corresponding mean-field simulation.

In the first two columns we see the convergence of the mean-field scheme, as expected the values are decreasing as the grid is refined. The four columns on the right illustrate the behavior as $N$ increases. All norm values are decreasing for increasing number of particles from $N=1000$ to $N=8000$.

Remark 6.4*.*

Since mean-field quantities are averaged it is very common to compare integral values like $J$ and $\rho^{N}$ in this setting. Note that the velocities and the trajectories of the agents are no such quantities. This explains why the convergence is more transparent for the integral quantities than for the velocities.

In Figure 5 the evolution of the crowd and the dogs is illustrated for the mean-field case $M50$ on the time interval $I=[0,55]$ with $\bar{V}_{\infty}=0.1\mathbb{V}_{\infty}(f_{0})$. The red lines show the trajectories of the dogs as above. Additionally, to the information in Figure 3, the red markers indicate the current positions of the dogs at different times. We see that the variance of the crowd is reduced first and then the center of mass is transported to the destination.

In Table 2 we compare the relative computation times and memory needed to compute 10 time steps of the optimization procedure. Replacing the 8000 particle simulation with the mean-field simulation with 50 grid points in each direction reduces the computational time by factor 10. Similar savings can be expected regarding the memory consumption.

7. Conclusions

Altogether, the numerical tests with different cost functional parameters show that the cost functional yields in the optimization the expected behaviour. Further, the results underline the convergence of the states and the convergence of the controls as $N$ tends to infinity. The convergence will be investigated from an analytical point of view in future work. Finally, comparing the first order necessary conditions, one realizes that there appear derivatives of $S$ on the microscopic level but not in the corresponding the mean-field equations. This relation will be discussed in some upcoming work as well.

Acknowledgements

MB acknowledges support by ERC via Grant EU FP 7 - ERC Consolidator Grant 615216 LifeInverse.