Mean Field Games Flock! The Reinforcement Learning Way

1 Introduction

The term flocking describes the behavior of large populations of birds that fly in compact groups with similar velocities, often exhibiting elegant motion patterns in the sky. Such behavior is pervasive in the animal realm, from fish to birds, bees or ants. This intriguing property has been widely studied in the scientific literature [Shaw, 1975; Reynolds, 1987] and its modeling finds applications in psychology, animation, social science, or swarm robotics. One of the most popular approaches to model flocking was proposed in [Cucker and Smale, 2007] and allows predicting the evolution of each agent's velocity from the speeds of its neighbors.

To go beyond pure description of population behaviours and emphasize on the decentralized aspect of the underlying decision making process, this model has been revisited to integrate an optimal control perspective, see e.g. [Caponigro et al., 2013; Bailo et al., 2018]. Each agent controls its velocity and hence its position by dynamically adapting its acceleration so as to maximize a reward that depends on the others' behavior. An important question is a proper understanding of the nature of the equilibrium reached by the population of agents, emphasizing how a consensus can be reached in a group without centralized decisions. Such question is often studied using the notion of Nash equilibrium and becomes extremely complex when the number of agents grows.

A way to approximate Nash equilibria in large games is to study the limit case of an continuum of identical agents, in which the local effect of each agent becomes negligible. This is the basis of the Mean Field Games (MFGs) paradigm introduced in [Lasry and Lions, 2007]. MFGs have found numerous applications from economics to energy production and engineering. A canonical example is crowd motion modeling in which pedestrians want to move while avoiding congestion effects. In flocking, the purpose is different since the agents intend to remain together as a group, but the mean-field approximation can still be used to mitigate the complexity.

However, finding an equilibrium in MFGs is computationally intractable when the state space exceeds a few dimensions. In traditional flocking models, each agent's state is described by a position and a velocity, while the control is the acceleration. In terms of computational cost, this typically rules out state-of-the-art numerical techniques for MFGs based on finite difference schemes for partial differential equations (PDEs) [Achdou and Capuzzo-Dolcetta, 2010]. In addition, PDEs are in general hard to solve when the geometry is complex and require full knowledge of the model.

For these reasons, Reinforcement Learning (RL) to learn control strategies for MFGs has recently gained in popularity [Guo et al., 2019; Elie et al., 2020; Perrin et al., 2020]. Combined with deep neural nets, RL has been used successfully to tackle problems which are too complex to be solved by exact methods [Silver et al., 2018] or to address learning in multi-agent systems [Lanctot et al., 2017]. Particularly relevant in our context, are works providing techniques to compute an optimal policy [Haarnoja et al., 2018; Lillicrap et al., 2016] and methods to approximate probability distributions in high dimension [Rezende and Mohamed, 2015; Kobyzev et al., 2020].

Our main contributions are: (1) we cast the flocking problem into a MFG and propose variations which allow multigroup flocking as well as flocking in high dimension with complex topologies, (2) we introduce the Flock'n RL algorithm that builds upon the Fictitious Play paradigm and involves deep neural networks and RL to solve the model-free flocking MFG, and (3) we illustrate our approach on several numerical examples and evaluate the solution with approximate performance matrix and exploitability.

2 Background

Three main formalisms will be combined: flocking, Mean Field Games (MFG), Reinforcement Learning (RL). Ez stands for the expectation w.r.t. the random variable z.

2.1 The Model of Flocking

To model a flocking behaviour, we consider the following system of N agents derived in a discrete time setting in [Nourian et al., 2011] from Cucker-Smale flocking modeling. Each agent i has a position and a velocity, each in dimension d and denoted respectively by xi and vi . We assume that it can control its velocity by choosing the acceleration, denoted by ui . The dynamics of agent i ∈ {1, . . . ,N} is:

where ∆t is the time step and ϵit is a random variable playing the role of a random disturbance. We assume that each agent is optimizing for a flocking criterion f that is underlying to the flocking behaviour. For agent i at time t, f is of the form:

where the interactions with other agents are only through the empirical distribution of states and velocities denoted by: µNt = 1N ∑Nj=1 δ(xjt ,vj t ) .

We focus on criteria incorporating a term of the form:

where β ≥ 0 is a parameter and (x, v, u, µ) ∈ Rd ×Rd ×Rd × P(Rd ×Rd ), with P(E) denoting the set of probability measures on a set E. This criterion incentivises agents to align their velocities, especially if they are close to each other. Note that β parameterizes the level of interactions between agents and strongly impacts the flocking behavior: if β = 0, each agent tries to align its velocity with all the other agents of the population irrespective of their positions, whereas the larger β > 0, the more importance is given to its closest neighbors (in terms of position).

In the N-agent case, for agent i, it becomes:

The actual criterion will typically include other terms, for instance to discourage agents from using a very large acceleration, or to encourage them to be close to a specific position. We provide such examples in Sec. 4.

Since the agents may be considered as selfish (they try to maximize their own criterion) and may have conflicting goals (e.g. different desired velocities), we consider Nash equilibrium as a notion of solution to this problem and the individual criterion can be seen as the payoff for each agent. The total payoff of agent i given the other agents' strategies u-i = (u1 , . . . , ui-1 , ui+1 , . . . , uN ) is: F iu-i(ui ) =

Exit ,vi t [ ∑t≥0 γ t f i t ], with f i t defined Eq. (1). In this context, a Nash equilibrium is a strategy profile (ˆu1 , ˆu2 , . . . , ˆuN ) such that there's no profitable unilateral deviation, i.e., for every i = 1, . . . ,N, for every control ui , F iˆu-i(ˆui ) ≥ F iˆu-i(ui ).

2.2 Mean Field Games

An MFG describes a game for a continuum of identical agents and is fully characterized by the dynamics and the payoff function of a representative agent. More precisely, denoting by µt the state distribution of the population, and by ξt ∈ Rℓ and αt ∈ Rk the state and the control of an infinitesimal agent, the dynamics of the infinitesimal agent is given by

where b : Rℓ × Rk × P(Rℓ ) → Rℓ is a drift (or transition) function, σ is a ℓ × ℓ matrix and ϵt+1 is a noise term taking values in Rℓ . We assume that the sequence of noises (ϵt)t≥0 is i.i.d. (e.g. Gaussian). The objective of each infinitesimal agent is to maximize its total expected payoff, defined given a flow of distributions µ = (µt)t≥0 and a strategy α (i.e., a stochastic process adapted to the filtration generated by (ϵt)t≥0) as: Jµ(α) = Eξt,αt [ ∑t≥0 γ tϕ(ξt, αt, µt)], where γ ∈ (0, 1) is a discount factor and ϕ : Rℓ ×Rk ×P(Rℓ ) → Rℓ is an instantaneous payoff function. Since this payoff depends on the population's state distribution, and since the other agents would also aim to maximize their payoff, a natural approach is to generalize the notion of Nash equilibrium to this framework. A mean field (Nash) equilibrium is defined as a pair (ˆµ, ˆα) = (ˆµt, ˆαt)t≥0 of a flow of distributions and strategies such that the following two conditions are satisfied: ˆα is a best response against ˆµ (optimality) and ˆµ is the distribution generated by ˆα (consistency), i.e.,

1. ˆα maximizes α ↦→ Jˆµ(α);
2. for every t ≥ 0, ˆµt is the distribution of ξt when it follows the dynamics (4) with (αt, µt) replaced by (ˆαt, ˆµt).

Finding a mean field equilibrium thus amounts to finding a fixed point in the space of (flows of) probability distributions. The existence of equilibria can be proven through classical fixed point theorems [Carmona and Delarue, 2018]. In most mean field games considered in the literature, the equilibrium is unique, which can be proved using either a strict contraction argument or the so-called Lasry-Lions monotonicity condition [Lasry and Lions, 2007]. Computing solutions to MFGs is a challenging task, even when the state is in small dimension, due to the coupling between the optimality and the consistency conditions. This coupling typically implies that one needs to solve a forward-backward system where the forward equation describes the evolution of the distribution and the backward equation characterizes the optimal control. One can not be solved prior to the other one, which leads to numerical difficulties. The basic approach, which consists in iteratively solving each equation, works only in very restrictive settings and is otherwise subject to cycles. A method which does not suffer from this limitation is Fictitious Play, summarized in Alg. 1. It consists in computing the best response against a weighted average of past

Algorithm 1: Generic Fictitious Play in MFGs

input : MFG = {ξ, ϕ, µ0}; # of iterations J ;

1 Define ¯µ0 = µ0 for j = 1, . . . , J do

• 2 1. Set best response αj = argmaxα J¯µj-1 (α) and let ¯αj be the average of (αi)i=0,...,j
• 3 2. Set µj = γ-stationary distribution induced by αj
• 4 3. Set ¯µj = j-1 j ¯µj-1 + 1j µj
• 5 return ¯αJ, ¯µJ

distributions instead of just the last distribution. This algorithm has been shown to converge for more general MFGs. State-of-the-art numerical methods for MFGs based on partial differential equations can solve such problems with a high precision when the state is in small dimension and the geometry is elementary [Achdou and Capuzzo-Dolcetta, 2010; Carlini and Silva, 2014]. More recently, numerical methods based on machine learning tools have been developed [Carmona and Lauri`ere, 2019; Ruthotto et al., 2020]. These techniques rely on the full knowledge of the model and are restricted to classes of quite simple MFGs.

2.3 Reinforcement Learning

The Reinforcement Learning (RL) paradigm is the machine learning answer to the optimal control problem. It aims at learning an optimal policy for an agent that interacts in an environment composed of states, by performing actions. Formally, the problem is framed under the Markov Decision Processes (MDP) framework. An MDP is a tuple (S,A, p, r, γ) where S is a state space, A is an action space, p : S ×A → P(S) is a transition kernel, r : S × A → R is a reward function and γ is a discount factor (see Eq. (5)). Using action a when the current state is s leads to a new state distributed according to P(s, a) and produces a reward R(s, a). A policy π : S → P(A), s ↦→ π(·|s) provides a distribution over actions for each state. RL aims at learning a policy π∗ which maximizes the total return defined as the expected (discounted) sum of future rewards:

with at ∼ π(·|st) and st+1 ∼ p(·|st, at). Note that if the dynamics (p and r) is known to the agent, the problem can be solved using e.g. dynamic programming. Most of the time, these quantities are unknown and RL is required. A plethora of algorithms exist to address the RL problem. Yet, we need to focus on methods that allow continuous action spaces as we want to control accelerations. One category of such algorithms is based on the Policy Gradient (PG) theorem [Sutton et al., 1999] and makes use of the gradient ascent principle: π ← π + α ∂R(π) ∂π , where α is a learning rate. Yet, PG methods are known to be high-variance because they use Monte Carlo rollouts to estimate the gradient. A vast literature thus addresses the variance reduction problem. Most of the time, it involves an hybrid architecture, namely Actor-Critic, which relies on both a representation of the policy and of the so-called state-action value function

(s, a) ↦→ Qπ (s, a). Qπ (s, a) is the total return conditioned on starting in state s and using action a before using policy π for subsequent time steps. It can be estimated by bootstrapping, using the Markov property, through the Bellman equations. Most recent implementations rely on deep neural networks to approximate π and Q (e.g. [Haarnoja et al., 2018]).

3 Our Approach

In this section, we put together the pieces of the puzzle to numerically solve the flocking model. Based on a mean-field approximation, we first recast the flocking model as an MFG with decentralized decision making, for which we propose a numerical method relying on RL and deep neural networks.

3.1 Flocking as anMFG

Mean field limit. We go back to the model of flocking introduced in Sec. 2.1. When the number of agents grows to infinity, the empirical distribution µNt is expected to converge towards the law µt of (xt, vt), which represents the position and velocity of an infinitesimal agent and have dynamics:

This problem can be viewed as an instance of the MFG framework discussed in Sec. 2.2, by taking the state to be the position-velocity pair and the action to be the acceleration, i.e., the dimensions are ℓ = 2d, k = d, and ξ = (x, v), α = u. To accommodate for the degeneracy of the noise as only the velocities are disturbed, we take σ = (0d 0d 0d 1d ) where 1d is

the d-dimensional identity matrix.

The counterpart of the notion of N-agent Nash equilibrium is an MFG equilibrium as described in Sec. 2.2. We focus here on equilibria which are stationary in time. In other words, the goal is to find a pair (ˆµ, ˆu) where ˆµ ∈ P(Rℓ ×Rℓ ) is a position-velocity distribution and ˆu : Rℓ × Rℓ ↦→ Rℓ is a feedback function to determine the acceleration given the position and velocity, such that: (1) ˆµ is an invariant position-velocity distribution if the whole population uses the acceleration given by ˆu, and (2) ˆu maximizes the rewards when the agent's initial position-velocity distribution is ˆµ and the population distribution is ˆµ at every time step. In mathematical terms, ˆu maximizes Jˆµ(u) = Ext ,vt ,ut [ ∑t≥0 γ tϕ(xt, vt, ut, ˆµ)], where (xt, vt)t≥0 is the trajectory of an infinitesimal agent who starts with distribution ˆµ at time t = 0 and is controlled by the acceleration (ut)t≥0. As the payoff function ϕ we use f flock β fromEq. (2). Moreover, the consistency condition rewrites as: ˆµ is the stationary distribution of (xt, vt)t≥0 if controlled by (ˆut)t≥0.

Theoretical analysis. The analysis of MFG with flocking effects is challenging due to the unusual structure of the dynamics and the payoff, which encourages gathering of the population. This is running counter to the classical LasryLions monotonicity condition [Lasry and Lions, 2007], which typically penalizes the agents for being too close to each other. However, existence and uniqueness have been proved in some cases. If β = 0, every agent has the same influence over the representative agent and it is possible to show that

Algorithm 2: Flock'n RL

input : MFG = {(x, v), f flock β , µ0}; # of iterations J

1 Define ¯µ0 = µ0 for j = 1, . . . , J do

• 2 1. Set best response πj = argmax π J¯µj-1 (π) with

SAC and let ¯πj be the average of (π0, . . . , πj)

• 3 2. Using a Normalizing Flow, compute µj =

γ-stationary distribution induced by πj

• 4 3. Using a Normalizing Flow and samples from (µ0, . . . , µj-1), estimate ¯µj
• 5 return ¯πJ, ¯µJ

the problem reduces to a Linear-Quadratic setting. Th 2 in [Nourian et al., 2011] shows that a mean-consensus in velocity is reached asymptotically with individual asymptotic variance σ2 2 . If β > 0, [Nourian et al., 2011] shows that if the MF problem admits a unique solution, then there exists an ϵN Nash equilibrium for theN agents problem and lim N→+∞ ϵN =

1. Existence has also been proved when β > 0 in [Carmona and Delarue, 2018, Section 4.7.3], with a slight modification of the payoff, namely considering, with ϕ a bounded func-

3.2 The Flock'n RL Algorithm

We propose Flock'n RL, a deep RL version of the Fictitious Play algorithm for MFGs [Elie et al., 2020] adapted to flocking. We consider a γ-discounted setting with continuous state and action spaces and we adapt Alg. 2 from its original tabular formulation [Perrin et al., 2020]. It alternates 3 steps: (1) estimation of a best response (using Deep RL) against the mean distribution ¯µj-1, which is fixed during the process, (2) estimation (with Normalizing Flows [Kobyzev et al., 2020]) of the new induced distribution from trajectories generated by the previous policy, and (3) update of the mean distribution ¯µj .

Computing the Best Response with SAC

The first step in the loop of Alg. 2 is the computation of a best response against ¯µj. In fact, the problem boils down to solving an MDP in which ¯µj enters as a parameter. Following the notation introduced in Sec. 2.3, we take the state and action spaces to be respectively S = R2d (for position-velocity pairs) and A = Rd (for accelerations). Letting s = (x, v) and a = u, the reward is: (x, v, u) = (s, a) ↦→ r(s, a) = f(x, v, u, ¯µj), which depends on the given distribution ¯µj. Remember that r is the reward function of the MDP while f is the optimization criterion in the flocking model.

As we set ourselves in continuous state and action spaces and in possibly high dimensions, we need an algorithm that scales. We choose to use Soft Actor Critic (SAC) [Haarnoja et al., 2018], an off-policy actor-critic deep RL algorithm using entropy regularization. SAC is trained to maximize a trade-off between expected return and entropy, which allows to keep enough exploration during the training. It is designed to work on continuous action spaces, which makes it suited for acceleration controlled problems such as flocking.

The best response is computed against ¯µj, the fixed average distribution at step j of Flock'n RL. SAC maximizes the reward which is a variant of f flock,i β,t from Eq. (3). It needs samples from ¯µj in order to compute the positions and velocities of the fixed population. Note that, in order to measure more easily the progress during the learning at step j, we sample N agents from ¯µj at the beginning of step 1 (i.e. we do not sample new agents from ¯µj every time we need to compute the reward). During the learning, at the beginning of each episode, we sample a starting state s0 ∼ ¯µj .

In the experiments, we will not need ¯πj but only the associated reward (see the exploitability metric in Sec. 4). To this end, it is enough to keep in memory the past policies (π0, . . . , πj) and simply average the induced rewards.

Normalizing Flow for Distribution Embedding

We choose to represent the different distributions using a generative model because the continuous state space prevents us from using a tabular representation. Furthermore, even if we could choose to discretize the state space, we would need a huge amount of data points to estimate the distribution using methods such as kernel density estimators. In dimension 6 (which is the dimension of our state space with 3-dimensional positions and velocities), such methods already suffer from the curse of dimensionality.

Thus, we choose to estimate the second step of Alg. 2 using a Normalizing Flow (NF) [Rezende and Mohamed, 2015; Kobyzev et al., 2020], which is a type of generative model, different from Generative Adversarial Networks (GAN) or Variational Autoencoders (VAE). A flow-based generative model is constructed by a sequence of invertible transformations and allows efficient sampling and distribution approximation. Unlike GANs and VAEs, the model explicitly learns the data distribution and therefore the loss function simply identifies to the negative log-likelihood. An NF transforms a simple distribution (e.g. Gaussian) into a complex one by applying a sequence of invertible transformations. In particular, a single transformation function f of noise z can be written as x = f(z) where z ∼ h(z). Here, h(z) is the noise distribution and will often be in practice a normal distribution.

Using the change of variable theorem, the probability density of x under the flow can be written as: p(x) = h(f -1 (x)) ∣∣∣det ( ∂f -1 ∂x )∣∣∣ . We thus obtain the probability distribution of the final target variable. In practice, the transformations f and f -1 can be approximated by neural networks. Thus, given a dataset of observations (in our case rollouts from the current best response), the flow is trained by maximizing the total log likelihood ∑n log p(x(n) ).

Computation of ¯µj

Due to the above discussion on the difficulty to represent the distribution in continuous space and high dimension, the third step (Line 4 of Alg. 2) can not be implemented easily. We represent every µj as a generative model, so we can not "average" the normalizing flows corresponding to (µi)i=1,...,j in a straightforward way but we can sample data points x ∼ µi for each i = 1, . . . , j. To have access to ¯µj , we keep in memory every model µj, j ∈ {1, . . . , J} and, in order to sample points according to ¯µj for a fixed j, we sample points from µi, i ∈ {1, . . . , j}, with probability 1/j. These points are then used to learn the distribution ¯µj with an NF, as it is needed both for the reward and to sample the starting state of an agent during the process of learning a best response policy.

4 Experiments

Environment. We implemented the environment as a custom OpenAI gym environment to benefit from the powerful gym framework and use the algorithms available in stable baselines [Hill et al., 2018]. We define a state s ∈ S as s = (x, v) where x and v are respectively the vectors of positions and velocities. Each coordinate xi of the position can take any continuous value in the d-dimensional box xi ∈ [-100,+100], while the velocities are also continuous and clipped vi ∈ [-1, 1]. The state space for the positions is a torus, meaning that an agent reaching the box limit reappears at the other side of the box. We chose this setting to allow the agents to perfectly align their velocities (except for the effect of the noise), as we look for a stationary solution.

At the beginning of each iteration j of Fictitious Play, we initialize a new gym environment with the current mean distribution ¯µj, in order to compute the best response.

Model - Normalizing Flows. To model distributions, we use Neural Spline Flows (NSF) with a coupling layer [Durkan et al., 2019]. More details about how coupling layers and NSF work can be found in the appendix.

Model - SAC. To compute the best response at each Flock'n RL iteration, we use Soft Actor Critic (SAC) [Haarnoja et al., 2018] (but other PG algorithms would work). SAC is an off-policy algorithm which, as mentioned above, uses the key idea of regularization: instead of considering the objective to simply be the sum of rewards, an entropy term is added to encourage sufficient randomization of the policy and thus address the exploration-exploitation trade-off. To be specific, in our setting, given a population distribution µ, the objective is to maximize: Jµ(π) = E(st ,ut) [∑+∞ t=0 γ t r(xt, vt, ut, µt) + δH(π(·|st))] , where H denotes the entropy and δ ≥ 0 is a weight.

To implement the optimization, the SAC algorithm follows the philosophy of actor-critic by training parameterized Qfunction and policy. To help convergence, the authors of SAC also train a parameterized value function V . In practice, the three functions are often approximated by neural networks.

In comparison to other successful methods such as Trust Region Policy Optimization (TRPO) [Schulman et al., 2015] or Asynchronous Actor-Critic Agents (A3C), SAC is expected to be more efficient in terms of number of samples required to learn the policy thanks to the use of a replay buffer in the spirit of methods such as Deep Deterministic Policy Gradient (DDPG) [Lillicrap et al., 2016].

Metrics. An issue with studying our flocking model is the absence of a gold standard. Especially, we can not compute the exact exploitability [Perrin et al., 2020] of a policy against a given distribution since we can not compute the exact best response. The exploitability measures how much an agent can gain by replacing its policy π with a best response π′ , when the rest of the population plays π: φ(π) := max π′ J(µ0, π′ , µπ ) - J(µ0, π, µπ ). If φ(¯πj) → 0 as j increases, FP approaches a Nash equilibrium. To cope with these issues, we introduce the following ways to measure progress of the algorithm:

• Performance matrix: we build the matrix M of performance of learned policies versus estimated distributions. The entry Mi,j on the i-th row and the jth column is the total γ-discounted sum of rewards: Mi,j = E [∑Tt=0 γ t rt,i | s0 ∼ ¯µi-1, ut ∼ πj(.|st)] , where rt,i = r(st, ut, ¯µi-1), obtained with πj against ¯µi-1. The diagonal termMj,j corresponds to the value of the best response computed at iteration j.
• Approximate exploitability: We do not have access to the exact best response due to the model-free approach and the continuous spaces. However, we can approximate the first term of φ(¯π) directly in the Flock'n RL algorithm with SAC. The second term, J(µ0, ¯π, µ¯π ), can be approximated by replacing ¯π with the average over past policies, i.e., the policy sampled uniformly from the set {π0, . . . , πj}. At step j, the approximate exploitability is ej = Mj,j - 1 j-1 ∑j-1 k=1 Mj,k. To smoothen the exploitability, we take the best response over the last 5 policies and use a moving average over 10 points. Please note that only relative values are important as it depends on the scale of the reward.

A 4-Dimensional example. We illustrate in a four dimensional setting (i.e. two-dimensional positions and velocities) how the agents learn to adopt similar velocities by controlling their acceleration. We focus on the role of β in the flocking effect. We consider noise ϵit ∼ N(0,∆t) and the following reward: ri t = f flock,i β,t -‖uit ‖22 +‖vi t ‖∞ - min{‖xi2,t ± 50‖}, where xi2,t stands for the second coordinate of the i-th agent's position at time t. The last term attracts the agents' positions towards one of two lines corresponding to the second coordinate of x being either -50 or +50. We added a term regarding the norm of the velocity to prevent agents from stopping. Here we take ‖vi t ‖∞ = max{|vi1,t |, |vi2,t |}. Hence, a possible equilibrium is with two groups of agents, one for each line. When β = 0, the term f flock,i β,t encourages agent i to have the same velocity vector as the rest of the whole population. At equilibrium, the agents in the two groups should thus move in the same direction (to the left or to the right, in order to stay on the two lines of x's). On the other hand, when β > 0 is large enough (e.g. β = 100), agent i gives more importance to its neighbors when choosing its control and it tries to have a velocity similar to the agents that are position-wise close to. This allows the emergence of two groups moving in different directions: one group moves towards the left (overall negative velocity) and the other group moves towards the right (overall positive velocity).

This is confirmed by Fig. 1. In the experiment, we set the initial velocities perpendicular to the desired ones to illustrate the robustness of the algorithm. We observe that the approximate exploitability globally decreases. In the case β = 0, we experimentally verified that there is always a global consen-

(a) Initial positions and velocities

(b) At convergence

• (c) Performance matrix
• (d) Approximate exploitability

Figure 1: Multi-group flocking with noise and β = 100.

sus, i.e., only one line or two lines but moving in the same direction.

Scaling to 6 dimensions and non-smooth topology. We now present an example with arbitrary obstacles (and thus non-smooth topology) in dimension 6 (position and velocity in dimension 3) which would be very hard to address with classical numerical methods. In this setting, we have multiple columns that the agents are trying to avoid. The reward has the following form: ri t = f flock,i β,t - ‖uit ‖22 + ‖vi t ‖∞ - min{‖xi2,t ‖} - c ∗ 1obs, If an agent hits an obstacle, it gets a negative reward and bounces on it like a snooker ball. After a few iterations, the agents finally find their way through the obstacles. This situation can model birds trying to fly in a city with tall buildings. In our experiments, we noticed that different random seeds lead to different solutions. This is not surprising as there are a lot of paths that the agents can take to avoid the obstacles and still maximizing the reward function. The exploitability decreases quicker than in the previous experiment. We believe that this is because agents find a way through the obstacles in the first iterations.

5 RelatedWork

Numerical methods for flocking models. Most work using flocking models focus on the dynamical aspect without optimization. To the best of our knowledge, the only existing numerical approach to tackle a MFG with flocking effects is in [Carmona and Delarue, 2018, Section 4.7.3], but it is restricted to a very special and simpler type of rewards.

Learning inMFGs. MFGs have attracted a surge of interest in the RL community as a possible way to remediate the scalability issues encountered in MARL when the number of agents is large [Perolat et al., 2018]. [Guo et al., 2019] com-

(b) At convergence

(c) Performance matrix

(d) Approximate exploitability

Figure 2: Flocking with noise and many obstacles.

bined a fixed-point method with Q-learning, but the convergence is ensured only under very restrictive Lipschitz conditions and the method can be applied efficiently only to finitestate models. [Subramanian and Mahajan, 2019] solve MFG using a gradient-type approach. The idea of using FP in MFGs has been introduced in [Cardaliaguet and Hadikhanloo, 2017], assuming the agent can compute perfectly the best response. [Elie et al., 2020; Perrin et al., 2020] combined FP with RL methods. However, the numerical techniques used therein do not scale to higher dimensions.

6 Conclusion

In this work we introduced Flock'n RL, a new numerical approach which allows solving MFGs with flocking effects where the agents reach a consensus in a decentralized fashion. Flock'n RL combines Fictitious Play with deep neural networks and reinforcement learning techniques (normalizing flows and soft actor-critic). We illustrated the method on challenging examples, for which no solution was previously known. In the absence of existing benchmark, we demonstrated the success of the method using a new kind of approximate exploitability. Thanks to the efficient representation of the distribution and to the model-free computation of a best response, the techniques developed here could be used to solve other acceleration controlled MFGs [Achdou et al., 2020] or, more generally, other high-dimensional MFGs. Last, the flexibility of RL, which does not require a perfect knowledge of the model, allow us to tackle MFGs with complex topologies (such as boundary conditions or obstacles), which is a difficult problem for traditional methods based on partial differential equations.