11institutetext:  Stellar Astrophysics Centre, Department of Physics and
Astronomy, Aarhus University, Ny Munkegade 120, DK-8000 Aarhus C, Denmark
22institutetext: School of Physics and Astronomy, Sun Yat-sen University, 2
Daxue Road, Tangjia, Zhuhai 519082, Guangdong Province, China 33institutetext:
Instituto de Astrofísica de Canarias (IAC), 38200 La Laguna, Tenerife, Spain
44institutetext: Niels Bohr Institute, University of Copenhagen, Juliane
Maries Vej 30, 2100 København Ø, Denmark 55institutetext: Niels Bohr Institute
& Centre for Star and Planet Formation, University of Copenhagen, Øster
Voldgade 5, DK-1350 – Copenhagen K, Denmark 66institutetext: Centre for
Electronic Imaging, Dept. of Physical Sciences. The Open University, Milton
Keynes MK7 6AA, UK

# The mass and age of the first SONG target: the red giant 46 LMi. ††thanks:
Based on observations made with the Hertzsprung SONG telescope operated at the
Spanish Observatorio del Teide on the island of Tenerife by the Aarhus and
Copenhagen Universities and by the Instituto de Astrofísica de Canarias.

S. Frandsen 11 M. Fredslund Andersen 11 K. Brogaard 11 C. Jiang 22 T. Arentoft
11 F. Grundahl 11 H. Kjeldsen 11 J. Christensen-Dalsgaard 11 E. Weiss 11 P.
Pallé 33 V. Antoci 11 P. Kjærgaard 44 A. N. Sørensen 44 J. Skottfelt 6655 U.
G. Jørgensen 55

(Received ; Accepted )

###### Abstract

Context.  The Stellar Observation Network Group (SONG) is an initiative to
build a worldwide network of 1m telescopes with high-precision radial-velocity
spectrographs. Here we analyse the first radial-velocity time series of a red-
giant star measured by the SONG telescope at Tenerife. The asteroseismic
results demonstrate a major increase in the achievable precision of the
parameters for red-giant stars obtainable from ground-based observations.
Reliable tests of the validity of these results are needed, however, before
the accuracy of the parameters can be trusted.

Aims. We analyse the first SONG time series for the star 46 LMi, which has a
precise parallax and an angular diameter measured from interferometry, and
therefore a good determination of the stellar radius. We use asteroseismic
scaling relations to obtain an accurate mass, and modelling to determine the
age.

Methods.  A 55-day time series of high-resolution, high S/N spectra were
obtained with the first SONG telescope. We derive the asteroseismic parameters
by analysing the power spectrum. To give a best guess on the large separation
of modes in the power spectrum, we have applied a new method which uses the
scaling of Kepler red-giant stars to 46 LMi.

Results. Several methods have been applied: classical estimates, seismic
methods using the observed time series, and model calculations to derive the
fundamental parameters of 46 LMi. Parameters determined using the different
methods are consistent within the uncertainties. We find the following values
for the mass $M$ (scaling), radius $R$ (classical), age (modelling), and
surface gravity (combining mass and radius): $M=1.09\pm 0.04$ ${\rm
M}_{\odot}$, $R=7.95\pm 0.11$ ${\rm R}_{\odot}$ age $t=8.2\pm 1.9$ Gy, and
$\log g=2.674\pm 0.013$.

Conclusions.  The exciting possibilities for ground-based asteroseismology of
solar-like oscillations with a fully robotic network have been illustrated
with the results obtained from just a single site of the SONG network. The
window function is still a severe problem which will be solved when there are
more nodes in the network.

###### Key Words.:

Stars: fundamental parameters – Stars, Individual: HD 94264 – Techniques:
radial velocities – Techniques: telescopes

##  1 Introduction

Asteroseismology has made a big step forward thanks to space-borne photometry
missions. This has been a revolution, particularly for red-giant stars
(Bedding et al. 2011; Beck et al. 2012; Christensen-Dalsgaard 2014), as the
need for very long time series has been solved by the long-duration
observations of CoRoT and Kepler. Ground-based campaigns suffer from much
shorter time coverage, and even for long campaigns (Arentoft et al., 2008) we
are far from the time coverage of space-borne missions. Furthermore, campaigns
must use high-precision radial velocities to detect the stellar oscillations
since ground-based photometry is unable to obtain the required precision
(Stello et al., 2006, 2007).

However, it is still useful to conduct ground-based radial-velocity campaigns,
if we can use dedicated networks of telescopes which can give close to 24-hour
coverage for extended periods as discussed in detail by Arentoft et al.
(2014). By measuring stellar oscillations using radial velocities the effects
of surface granulation are strongly reduced compared to photometric
measurements. This allows us to detect modes at lower frequencies, and modes
of $l\,=\,3$ are easier to detect using radial velocities.

We are particularly interested in nearby stars where an angular diameter can
be measured accurately with an interferometer. The idea of combining
interferometric measurements with observations of solar-type oscillations in
the context of red-giant stars was employed for the star $\epsilon$ Oph
(Mazumdar et al. 2009), where the asteroseismic data were obtained from space
with the photometric MOST mission (Walker et al., 2003). This was a favourable
case because the angular diameter was measured with high accuracy and the
space data do not suffer from atmospheric and alias problems. This is
reflected in the good agreement found between the interferometric radius and
the asteroseismic radius. Another similar study has been presented by Beck et
al. (2015) for the red-giant stars $\gamma$ Per and $\theta^{1}$ Tau, where
the asteroseismic measurements come from multisite radial-velocity
measurements.

The Stellar Observations Network Group (SONG) is an initiative which aims to
establish a ground-based network of automated 1m telescopes to obtain high-
precision radial-velocity measurements of nearby, bright stars and to study
their properties, primarily using the technique of asteroseismology. The
targets being studied will be close enough that accurate parameters, such as
angular diameter, parallax, temperature, and abundances, are well known in
advance or can be easily determined.

The first SONG node consists of a 1m telescope produced by ASTELCO GmBH and is
installed at the Observatorio del Teide, Tenerife, Spain. There are two
instruments: a lucky imaging photometer (Skottfelt et al. 2015) and a high-
resolution spectrograph (R $\sim$ 100,000) installed at the Coudé focus in a
container next to the dome. The spectrograph, which is the instrument relevant
for this paper, has been optimized for high-precision radial-velocity
measurements and employs an iodine cell for wavelength reference. All
operations are fully automatic and no operator is needed on site or remotely.
For more details see Andersen et al. (2014) & Andersen et al. (2016). The SONG
radial-velocity pipeline is described in Grundahl et al. (2017). In this paper
we present the results of the first long test run with the SONG node on
Tenerife.

Table 1: Basic parameters of 46 LMi Parameter | Measured value | Reference
---|---|---
$V$ | 3.83 | SIMBAD Wenger et al. (2000)
Temperature $T_{\rm eff}$ | 4690$\pm 50\,$K | Bubar & King (2010)
Metal abundance [Fe/H] | -0.1$\pm$0.1 | Bubar & King (2010)
Parallax (mas) | $34.38\pm 0.21$ | van Leeuwen (2007)
Angular diameter (mas) | $2.54\pm 0.03$ | Nordgren et al. (1999)
Radius $R$/${\rm R}_{\odot}$ | 7.95$\pm$0.11 | from parallax and angular diameter
Surface gravity $\log g$ | $2.61\pm 0.2$ | Bubar & King (2010)
Luminosity $L$/${\rm L}_{\odot}$ | $27.42\pm 1.38$ | from $T_{\rm eff}$ and $R$ above
Rotation $v\sin(i)$ |  $2.1$ km s-1 | Massarotti et al. (2008)
Possible age (see text) |  $2.7\pm 0.5$ Gy | Bubar & King (2010)

##  2 The first SONG target 46 LMi

The bright ($V=3.83$) red giant 46 LMi (HR 4247, HD 94264, HIP 53229) with
spectral type K0III was selected as a suitable first target for SONG as the
amplitudes and the timescale of the solar-like oscillations in this star are
large. These requirements were necessary to test the telescope and
spectrograph performance.

The star also presents a good test case for asteroseismology as there are
several high-precision measurements of its basic parameters (Table 1). This
includes a well-determined parallax ($\pi=34.38\pm 0.21$ mas, van Leeuwen
2007), and an angular diameter ($\Theta_{L}=2.54\pm 0.03$ mas) measured by
interferometry (Nordgren et al., 1999). Together, these lead to a classical
measurement of the radius of $R=7.95\pm 0.11$ ${\rm R}_{\odot}$.

In Table 1 the selected values from the literature are listed. The estimate of
the age comes from the possible membership of the moving group WOLF 630 (Bubar
& King, 2010) of 46 LMi.

##  3 Observations and data reduction

Figure 1: Time coverage of the 46 LMi campaign. First night: 14 February 2014.

The SONG node permitted a very good time coverage, but with just one node
operational the data still suffer from gaps in the time series during daytime.
Nevertheless, we were able to obtain a very long time series without involving
any observer.

We obtained a total of 3211 exposures, each with an exposure time of 240 s and
a CCD readout time of 2.5 s. The observations were carried out from 14
February to 27 April 2014, with data obtained on 35 of the nights. The main
causes for interruptions were poor weather conditions and technical work being
done on the telescope and control system. All observations were carried out
automatically; in the late afternoon we inserted the observing requests for
the coming night into our queue system, and these requests were then executed
during the night. All calibration files (bias, flats, ThAr reference) were
obtained automatically each afternoon prior to the start of the observations.
Figure 1 shows the distribution of data points over the 55-day period
allocated to this target. The large gap is due to a period of necessary
technical upgrades; the data were obtained while the telescope was being
commissioned. The rest of the time the only interruptions were due to bad
weather.

The spectrograph has six possible slits, and for all observations we used the
one providing a resolution of 90,000. For accurate wavelength calibration, we
used an iodine cell. All spectra were bias corrected, flat-fielded and
wavelength calibrated using the REDUCE package implemented in IDL (see
Piskunov & Valenti
2002)111http://www.astro.uu.se/$\sim$piskunov/RESEARCH/REDUCE/. To extract the
precise radial velocities we used the software iSONG (also IDL based) (see
e.g. Corsaro et al. (2012), Antoci et al. (2013), and Grundahl et al. (2017)
for more details). There are a total of 24 useful orders with iodine lines
available for calculating the velocities.

Figure 2: Radial velocity curve of a single night.

In Fig. 2 we present a full night of the radial-velocity measurements. The
typical characteristic timescale of variability of the stochastic oscillations
is clearly seen in the velocity curve as the noise level per data point is
much lower than the periodic signal. Typically, the instrumental noise per
data point is 1.5 ${\rm m\,s}^{-1}$ on a good night, where the noise has been
calculated as the scatter around a smooth curve through the data.

##  4 Analysis of the radial-velocity time series

###  4.1 Frequency of maximum power: $\nu_{\rm max}$

The power spectrum (shown in Fig. 3) was calculated as a fit of sinusoids to
the raw radial-velocity time series, as described in Frandsen et al. (1995).
The oscillation signal is clearly visible as the power excess around
$60\,\mu$Hz. A low-frequency filtering has been applied where ten frequencies
below $6\,\mu$Hz were cleaned from the power spectrum. The effect of single-
site observations is illustrated by the window function in the inset, which
shows strong alias peaks at $\pm 11.57\,\mu$Hz corresponding to the daily gap
in the time series. This significantly complicates the determination of the
frequencies of the $p$ modes.

Figure 3: High-pass filtered power spectrum of 46 LMi. The inset shows the
normalized window function, which shows strong sidelobes.

The frequency of maximum power ($\nu_{\rm max}$) is a widely used value which
relates to the large frequency separation through scaling relations (Stello et
al., 2009; Huber et al., 2010; Handberg et al., 2016a; Arentoft et al., 2017).
Based on the available data we determined $\nu_{\rm max}$ from a Gaussian fit
to the high-pass filtered power spectrum. The uncertainty on $\nu_{\rm max}$
was determined by splitting the time series into seven chunks, each with the
same number of data points. For each chunk the unfiltered power spectrum was
calculated and a value for $\nu_{\rm max}$ was determined by fitting a
Gaussian to the interval showing power excess, $30\,\mu$Hz – $100\,\mu$Hz. The
uncertainty on the final frequency of maximum power was then given as the
standard deviation of the mean error. This procedure yielded $\nu_{\rm
max}=59.4\pm 1.4\,\mu$Hz. It is worth mentioning that the expected mode
lifetime is of the order of the time span of the data set. This can have an
effect on the determined $\nu_{\rm max}$, which will be highly influenced by
the dominant mode in the power spectrum.

Figure 4: Échelle diagram with respect to the daily alias $11.57\,\mu$Hz,
which brings out peaks belonging to the window function. This is the basis for
locating a number of possible modes listed in Table 2. The arrow points to the
five peaks corresponding to the dominant mode in the power spectrum.

###  4.2 Individual frequencies

We have attempted to determine the frequencies of modes in the power spectrum
of 46 LMi using two techniques. In the first method, illustrated in Fig. 4,
peaks found in a smoothed power spectrum are plotted in an échelle diagram
with frequency modulo $11.574\,\mu$Hz on the abscissa and $\nu$ on the
ordinate. Then modes are found by looking for patterns similar to the window
function in the vertical direction. The size of the symbols represents the
amplitude of each detected peak and can be used to select the correct
frequency, but in some cases an alias might be obtained instead. Some of the
peaks shown in the figure are statistically non-significant and are only
included to look for the aliases in the window pattern.

In the second method an iterative frequency determination is performed using
PERIOD04. First, the frequency of the mode with maximum signal is found. A
weighted cosine fit to the time series with the corresponding period is
calculated and subtracted from the time series. This way side bands caused by
the window function are removed in the power spectrum. This is repeated on the
residual spectrum, but the following iterations simultaneously fit multiple
cosines to the time series. This operation is halted when the same number of
frequencies as in the first method has been reached. This method is possible
because the resolution in the power spectrum ($0.2\,\mu$Hz) is larger than the
line width of the stochastic modes in red giants ($\approx 0.1\,\mu$Hz,
Handberg et al. 2016b). This means excited modes will be represented by only
one peak in the power spectrum in most cases and not by a set producing a
Lorentzian profile, which is the case when the observations span many life
cycles.
Finally, we get a list of modes (frequency, amplitude, and phase), where modes
are only included if they are detected by both methods and agree to better
than 0.1 $\mu$Hz. Both results ($\nu_{1},\nu_{2}$) are presented in Table 2 to
give an idea of the precision, and generally good agreement is seen. The
amplitude $A_{2}$ is given for the PERIOD04 method. For two frequencies
PERIOD04 detects what could be an alias, which we have enclosed in
parentheses. All modes included in Table 2 have high signal-to-noise ratios
($S/N$). More modes are likely present; however, for the remaining peaks we
were not able to determine which are the true modes and which are the aliases,
and they have therefore not been included.

[b]  $\nu_{1}$ $\nu_{2}$ $A_{2}$ 1 S/N 2 ($\mu$Hz) ($\mu$Hz) (${\rm
m\,s^{-1}}$) 58.986 58.958 5.1 72.9 53.705 53.703 2.3 32.9 64.838 64.837 1.9
27.1 60.919 60.914 1.9 27.1 51.160 51.157 1.9 27.1 64.148 64.143 1.7 24.3
71.345 (82.904) – ${\sim}25$ 48.796 (60.451) – ${\sim}25$

Table 2: Possible modes in 46 LMi

  * 1

$\nu_{1}$ is the frequency from method 1. $\nu_{2}$ and $A_{2}$ the frequency
and amplitude from method 2 (PERIOD4).

  * 2

Noise level calculated at high frequency to 0.07m/s in PERIOD04.

###  4.3 Solution using $\nu_{\rm max}$

The global asteroseismic observables include $\nu_{\rm max}$ defined above and
the large frequency separation $\Delta\nu$, i.e. the average frequency spacing
between modes of the same degree and adjacent orders.

As presented by Miglio et al. (2012), the asteroseismic scaling relations for
mass can be expressed in the following four different ways if the luminosity
$L$ is included as an independent variable along with $T_{\rm eff}$,
$\Delta\nu$, and $\nu_{\rm max}$:

$\frac{M}{{\rm
M}_{\odot}}\simeq\left({\frac{\Delta\nu}{\Delta\nu_{\odot}}}\right)^{-4}\left(\frac{\nu_{\rm
max}}{\nu_{\rm
max,\odot}}\right)^{3}\left(\frac{{T_{\mathrm{eff}}}}{{T_{\mathrm{eff},\odot}}}\right)^{3/2}\;,$
$\frac{M}{{\rm
M}_{\odot}}\simeq\left(\frac{\Delta\nu}{\Delta\nu_{\odot}}\right)^{2}\left(\frac{L}{{\rm
L}_{\odot}}\right)^{3/2}\left(\frac{{T_{\mathrm{eff}}}}{{T_{\mathrm{eff},\odot}}}\right)^{-6}\;,$
$\frac{M}{{\rm M}_{\odot}}\simeq\left(\frac{\nu_{\rm max}}{\nu_{\rm
max,\odot}}\right)\left(\frac{L}{{\rm
L}_{\odot}}\right)\left(\frac{{T_{\mathrm{eff}}}}{{T_{\mathrm{eff},\odot}}}\right)^{-7/2}\;,$
$\frac{M}{{\rm M}_{\odot}}\simeq\left(\frac{\nu_{\rm max}}{\nu_{\rm
max,\odot}}\right)^{12/5}\left(\frac{\Delta\nu}{\Delta\nu_{\odot}}\right)^{-14/5}\left(\frac{L}{{\rm
L}_{\odot}}\right)^{3/10}\;.$

Using $L$ we exploit the asteroseismic scaling relation for mass in the form
that does not involve $\Delta\nu$. We use the standard asteroseismic values
and $T_{\rm eff}$ for the Sun (Miglio et al., 2012), which are
$\Delta\nu_{\odot}=135\,\mu$Hz, $\nu_{\rm max,\odot}=3100\,\mu$Hz, and $T_{\rm
eff,\odot}=5777\,K$.

We adopted the luminosity ($L=27.42\pm 1.38$ ${\rm L}_{\odot}$) determined by
combining the spectroscopic $T_{\rm eff}$ and radius, $R=7.95\pm 0.11$ ${\rm
R}_{\odot}$ from angular diameter and parallax distance since in that way we
avoid making use of a bolometric correction. We checked that the luminosity
derived in an alternative way from the parallax $\pi=34.38$ mas, the visual
magnitude $V=3.83,$ and the bolometric corrections from Casagrande &
VandenBerg (2014), $BC_{V}=-0.384$ and ${BC}_{V,\odot\ }=-0.068,$ are in
agreement with the adopted values to well within the $1\sigma$ uncertainty.

With this approach we obtained $M=1.09\pm 0.04$ ${\rm M}_{\odot}$, after which
reversing any of the other mass equations leads to $\Delta\nu=6.30\pm
0.08\,\mu$Hz with the uncertainties determined by bootstrapping (Feigelson &
Babu, 2012). With the estimate of the mass and radius we can determine the
surface gravity of 46 LMi to $\log g=2.674\pm 0.013$.

It is well known that a correction to the observed $\Delta\nu$ is needed
before its use with the asteroseismic scaling relations (Brogaard et al.,
2016; Handberg et al., 2016a; Miglio et al., 2012; Mosser et al., 2013;
Guggenberger et al., 2016). From figure 3 in Rodrigues et al. (2017) for a
mass of 1.09 ${\rm M}_{\odot}$ and $\nu_{\rm max}=59.4\,\mu$Hz in the two
cases with [Fe/H] = -0.25 and [Fe/H] = 0.00 to match the [Fe/H] = -0.1 of 46
LMi we get a correction factor of 0.966 as the mean of the two solutions.
Applying this correction to $\Delta\nu=6.30\pm 0.08\,\mu$Hz we obtain an
estimate for the observed $\Delta\nu=6.09\pm 0.09\,\mu$Hz. The increase in
uncertainty comes from the metallicity dependence of the correction.

###  4.4 Estimating $\Delta\nu$ using other methods

We tried to find the large frequency separation $\Delta\nu$ by plotting the
detected modes in an échelle diagram with different guesses for the large
separation, but no convincing $\Delta\nu$ could be found in this way due to
the significant number of alias peaks. Another unsuccessful approach was to
calculate the autocorrelation of the power spectrum, and we therefore tried
the method that we now discuss here.

####  4.4.1 Comparison to red giants observed by Kepler

We made use of the wealth of information available from the Kepler
observations (Jenkins et al., 2010). From the APOKASC222APOKASC comes from a
combination of Apache Point Observatory Galactic Evolution Experiment (APOGEE)
and of the Kepler Asteroseismic Science Consortium (KASC) catalogue
(Pinsonneault et al., 2014) we chose five red giants observed by Kepler with
similar APOKASC effective temperatures and radii to 46 LMi (See Table 3). The
asteroseismic masses of the five stars are close to the value we determined in
Sect. 4.3 for 46 LMi. The idea is to apply a scale factor to the power
spectrum of a Kepler red giant to match the power spectrum of 46 LMi, but
first the method was tested on the Kepler stars alone. The frequencies in the
power spectrum of each Kepler star was multiplied by a variable scale factor,
and the cross-correlation between the smoothed power spectrum and a chosen
smoothed reference spectrum (KIC 4672904) was determined for each value of the
scale factor. Each power spectrum was normalized with the power of the highest
peak in the region of the oscillation signal. A first guess on the scale
factor was based on the $\Delta\nu$ given in APOKASC and the cross-correlation
was done for a range around this value. For each star the maximum in the
cross-correlation function was found and this value of the scale factor was
applied to the frequencies in the power spectrum. The result of this test (see
Fig. 5) demonstrates how well all modes align and proves that the method works
for this sample of Kepler stars. We also tested the method on a larger sample
of more than 50 Kepler stars which were not as similar to 46 LMi, but we still
obtained a good agreement.
It is important to mention that applying the determined scale factor for the
individual Kepler stars to the $\Delta\nu$ of the reference does not produce
the exact same $\Delta\nu$ as given in APOKASC but they are well within the
$1\sigma$ uncertainties. The scaled values are given in Table 3.

Figure 5: Four red giants scaled to KIC 4672904. The scale factor applied to
the frequency axis of each power spectrum is shown. The l-mode label marks the
mode identification of the dominant mode in the given spectrum. Figure 6:
Cross-correlation function for each of the Kepler stars correlated with 46
LMi. The abscissa is the scale factor multiplied by the $\Delta\nu$ of KIC
4672904. The dashed lines indicate the different $l$-values coming from the
identified modes in the Kepler stars when scaled to the dominant mode in 46
LMi. See the text for further information.

The method was then applied to 46 LMi which was used as the reference for the
cross-correlation. In Fig. 6 the cross-correlation function for each of the
Kepler stars from Table 3 is shown. The correlation to 46 LMi is almost
entirely dominated by the highest peak in the 46 LMi power spectrum producing
a highly correlated signal when matched with a central peak (around 60 $\mu$Hz
in Fig. 5) in the Kepler red giants. This results in three dominant peaks in
the cross-correlation corresponding to $l=0,1$, and 2 from the Kepler stars
scaled to the dominant peak in 46 LMi. In the case of KIC 8545291 and KIC
3734104 the dominant mode in the power spectrum is outside the central region
(See Fig. 5). Nevertheless, it is still one of the modes in the central region
that produces the highest correlation. This is illustrated in Fig. 7 where the
cross-correlation between KIC 3734104 and 46 LMi is shown with a larger range
of scale factors. The first peak in this figure corresponds to the match when
the dominant mode in KIC 3734104 is scaled to the dominant mode in 46 LMi.
Even though this mode is almost twice as high in the power spectrum the cross-
correlation is still stronger around the central modes around 60 $\mu$Hz.

Figure 7: Cross-correlation between KIC 3734104 and 46LMi for a larger range
of scale factors. The peak at 0.95 corresponds to the match when the dominant
mode in KIC 3734104 (see bottom panel in Fig. 5) is scaled to the dominant
mode in 46 LMi.

The three possible values of $\Delta\nu$ from the correlation corresponding to
the different l-values can be determined with high accuracy (better than
$0.01\,\mu$Hz), but the real uncertainty on these values comes from the
uncertainty in the $\Delta\nu$ from the APOKASC catalogue. We selected KIC
4672904 as our reference for the scaled $\Delta\nu$ in Table 3 because it is
the one closest to 46 LMi in the physical parameters. These values for
$\Delta\nu$ are used to plot the correlation as a function of frequency on the
abscissa on Fig. 6. This is simply done by multiplying the given $\Delta\nu$
with the scale-factor values found from the cross-correlations. This means
that if we selected another star as our reference, the abscissa could shift up
to $0.05\,\mu$Hz (in the case of KIC 6844373).

From this scaling method we get three possible values for the large frequency
separation for 46 LMi, which we take as the mean of the five scaled
$\Delta\nu$ from Table 3 with the uncertainty given as the half of the full
range in these values:

$\Delta\nu_{46LMi}=(5.85\pm 0.03)\,\mu$Hz (Dominant mode in 46 LMi is $l$=0)

$\Delta\nu_{46LMi}=(5.92\pm 0.03)\,\mu$Hz (Dominant mode in 46 LMi is $l$=2)

$\Delta\nu_{46LMi}=(6.12\pm 0.03)\,\mu$Hz (Dominant mode in 46 LMi is $l$=1)

We cannot say which $\Delta\nu$ is the correct one because we do not know the
identification of the dominant mode in 46 LMi, but by looking at each value we
can argue which are more or less likely. Remembering the value of $\nu_{\rm
max}$ found in Sec. 4.3, and using this with the known physical parameters of
46 LMi and the asteroseismic scaling relations, we found $\Delta\nu=6.09\pm
0.09\,\mu$Hz. From this, we can argue that the $\Delta\nu$ corresponding to
the highest peak being a $l$=1 mode is the most likely. This solution has some
implications which do not support it. When the power spectrum from any of the
Kepler stars is scaled using the $l$=1 value, none of their radial modes near
$\nu_{\rm max}$ match peaks in the power spectrum of 46 LMi. This seems
unlikely since no Kepler red giant in our sample shows this behaviour.
However, it is well known that radial modes have less power relative to the
dipole modes in RV compared to intensity (Bedding et al., 1996), which again
could favor this solution. The relative ratio between $l$=0 and $l$=1 modes
for SONG was computed in Handberg & Campante (2011) and was determined to be
1.35. The solution where the dominant mode in 46 LMi is a $l$=0 mode is less
likely (${\sim}3\sigma$) when compared to the $\Delta\nu$ estimated from
$\nu_{\rm max}$ in Sec. 4.3. However, more modes from the power spectrum of a
given Kepler giant match peaks in the power spectrum of 46 LMi. This could be
a result of the $\Delta\nu$ being very close to one half of the daily
splitting produced by the sampling from single-site observations. Every second
mode in the Kepler power spectrum would then roughly match an alias peak if
this were the true value and the high correlation for the $l$=0 mode might
simply be a result of the sampling of our data.
The same is the case for the value related to the dominant mode in 46 LMi
being a $l$=2 mode. This value is also close to one half of the daily alias.
The fact that it is a bit further away from the alias peaks might mirror the
correlation power in Fig. 6, which is in general a bit lower for the $l$=2
peaks than the $l$=0 modes.

The high peaks in Fig. 6 have their origin in the dominant modes in the power
spectra of 46 LMi and the template. To reduce the effect of the dominant modes
we have divided the power spectra by a Gaussian fit to the power spectra
before cross correlation. We have also performed a test with a power law. In
none of the cases did we see any significant change of the results. The
heights of the peaks in Fig. 6 were modified, but not the positions.

All three values of $\Delta\nu$ can be used to scale a given Kepler star to 46
LMi, and all three fail to match all frequencies given in Table 2. This
indicates that we have either incorrectly determined some of the frequencies
given or that none of the stars in our sample matches 46 LMi completely, or a
combination of the two. We know from the Kepler red giants that some show
clear splitting of the dipole modes caused by mixing with g-modes. If this
were the case in 46 LMi it would make the analysis of the power spectrum
significantly more difficult since all split modes would produce alias peaks
due to the window function.

Table 3: Large separation of 46 LMi, $\Delta\nu_{46LMi}$ in $\mu$Hz from
scaling the large separation $\Delta\nu$ of Kepler giants

Kepler target |  $\Delta\nu\leavevmode\nobreak\ (\mu$Hz)a |  $\Delta\nu_{scaled}\leavevmode\nobreak\ (\mu$Hz)b |  $\Delta\nu_{46Lmi}\leavevmode\nobreak\ (\mu$Hz)c |  $T_{\rm eff}$(K)a |  $M$/${\rm M}_{\odot}$ a |  $R$/${\rm R}_{\odot}$ a | [Fe/H]a
---|---|---|---|---|---|---|---
KIC 4672904 | $6.05\pm 0.13$ | - | 5.82, 5.89, 6.09 | 4614 | $1.11\pm 0.11$ | 8.19 | -0.38
KIC 6844373 | $5.78\pm 0.15$ | 5.73 | 5.87, 5.94, 6.14 | 4556 | $1.10\pm 0.15$ | 8.41 | -0.05
KIC 8545291 | $5.81\pm 0.15$ | 5.78 | 5.86, 5.93, 6.13 | 4522 | $1.09\pm 0.13$ | 8.34 | -0.13
KIC 4548564 | $5.98\pm 0.13$ | 5.99 | 5.84, 5.91, 6.11 | 4503 | $1.05\pm 0.12$ | 8.09 | -0.33
KIC 3734104 | $5.85\pm 0.14$ | 5.82 | 5.86, 5.93, 6.13 | 4514 | $1.21\pm 0.17$ | 8.60 | 0.02

  * a

Values from Pinsonneault et al. (2014); $T_{\rm eff}$ & [Fe/H] are the
UNCORRECTED ASPCAP (G1)

  * b

Scaled from KIC 4672904

  * c

Scaled from the given Kepler star based on the $l=0$, 2, and 1 maximum in the
correlation function, respectively.

##  5 Modelling

We can show that a model exists that reproduces our new asteroseismic results
and existing non-seismic measurements. In addition, we can find a ‘best’
solution in terms of all constraints.

We employed the ASTEC evolution code (Christensen-Dalsgaard, 2008b) for
stellar evolution modelling computations, and the ADIPLS oscillation package
(Christensen-Dalsgaard, 2008a) for frequency calculations. ASTEC can evolve
models up to the tip of the red giant branch, which is adequate for 46 LMi.
The input physics of the ASTEC version that we used included the latest OPAL
opacity tables (Iglesias & Rogers, 1996), OPAL 2005 equation of state (Rogers
& Nayfonov, 2002), and NACRE reaction rates (Angulo et al., 1999). At low
temperatures, opacities are obtained from Ferguson et al. (2005). Convection
is treated under the assumption of mixing length theory (Böhm-Vitense, 1958).
We did not take rotation, diffusion, or convective overshoot into
consideration in our calculation.

For model selection, we used a tool called DIAMONDS that performs Bayesian
parameter estimation and model comparison by means of the nested sampling
Monte Carlo (NSMC) algorithm (Corsaro & De Ridder, 2014). The mass, heavy-
element abundance $Z$, and mixing-length parameter $\alpha$ are set as free
parameters. The mass is restricted to the range of 0.95–1.25 ${\rm
M}_{\odot}$, $Z$ to the range 0.001 – 0.03, and $\alpha$ to 1 – 3. Based on a
Galactic chemical-evolution model (Carigi, 2000; Pietrinferni et al., 2004)
the hydrogen abundance $X$ is obtained as $X=-2.4Z+0.748$, which gives $X=0.7$
when $Z=0.02$ and a helium-to-metal enhancement ratio $\Delta Y/\Delta Z=1.4$.
We use the observed parameters (with their uncertainties) $R$, $T_{\rm eff}$
from Table 1, and $\nu_{\rm max}$ from Sect. 4.1 as constraints to obtain the
likelihood of a given model based on the $\chi^{2}$ value. The results are
given in Table 4 where the errors are given as Bayesian errors which only
represent the capability of the model to fit the given observations. We use a
68.3% probability for the Bayesian credible intervals.

Table 4: Most likely model parameters given the constraints from observations. See the text for explanation of the errors given. Parameter | Value |  $1\sigma$ error
---|---|---
$M$/${\rm M}_{\odot}$ | 1.09 | 0.03
$Z$ | 0.017 | 0.003
$\alpha$ | 2.07 | 0.11
$R$/${\rm R}_{\odot}$ | 7.95 | 0.13
$T_{\rm eff}$ (K) | 4688 | 86
$L$/${\rm L}_{\odot}$ | 27.3 | 2.2
Age (Gyr) | 8.2 | 1.9
$\log g$ | 2.67 | 0.01
$\nu_{\rm max}\,(\,\mu{\rm Hz})$ | 59.4 | 2.0
$\Delta\nu\,(\,\mu{\rm Hz})$ | 6.3 | 0.1

If we assume that we know $R$ and $T_{\rm eff}$, we can use the scaling
relation for $\nu_{\rm max}$ to calculate $M$ which we did in Sec. 4.3. Not
surprisingly, the best fitting model therefore has parameters consistent with
the previous results in this paper. The $\chi^{2}$ of the best models is 1.36.
The values given in Table 4 are from the Bayesian and they give
$\chi^{2}=0.87$. Removing the constraint on $\nu_{\rm max}$ changes the
results very little. The age determined here is significantly older than found
by Bubar & King (2010), but it is the expected age for a $\sim 1.09$ ${\rm
M}_{\odot}$ red giant. If 46 LMi were indeed a member of the moving group Wolf
630, this would suggest a higher age for the group. More likely, 46 LMi is not
a member.

##  6 Conclusions

Even if SONG is not yet a network of telescopes, the use of just one dedicated
automatic telescope leads to ground-based results comparable to or better than
previous asteroseismic campaigns. The reason is the easy way SONG can perform
long time observations. The main results are summarized here.

  * •

The first SONG observations demonstrate that we have perhaps obtained the best
ground-based red-giant power spectrum to date. 46 LMi shows a clear solar-like
power excess centred at 59.4$\pm 1.4\,\mu$Hz, where the uncertainty is mainly
due to the stochastic excitation. The power reaches nearly zero at low
frequencies ($\leq 30\,\mu$Hz) for the high-pass filtered radial-velocity time
series. The power spectrum is affected by alias peaks due to the use of
single-site data. Comparison to stars observed by Kepler leads to three
possible values for the large frequency separation. We therefore base our mass
and radius estimates on $\nu_{\rm max}$ as the only asteroseismic measure and
combine it with classical results, as described in Sect. 4.3.

  * •

Our best estimate of the mass of 46 LMi (from the $\nu_{\rm max}$ solution) is
$M=1.09\pm 0.04$ ${\rm M}_{\odot}$ and the radius from classical estimates is
$R=7.95\pm 0.11$ ${\rm R}_{\odot}$. The surface gravity estimate is $\log
g=2.674\pm 0.013$ derived from the mass and radius just listed leading to an
improved accuracy compared to the spectroscopic gravity $\log g=2.61\pm 0.2$.

  * •

The age of $8.2\pm 1.9$ Gyr that we find for 46 LMi is significantly higher
than that reported by Bubar & King (2010). This suggests that 46 LMi is not a
member of the moving group WOLF 630. We thus demonstrate that asteroseismic
age determination can be used to determine membership of stellar
constellations such as moving groups and clusters.

Improvements will be possible from multisite observations, as we will be able
to identify individual mode frequencies since alias peaks will be highly
suppressed. This will lead to a more precise age estimate, and a measurement
of mass and radius from the asteroseismic scaling relations that is
independent of the distance. We will then be able to get a high-accuracy test
of the asteroseismic scaling relations by comparing the radius derived from
asteroseismology to that obtained from interferometry. This clearly shows why
more SONG nodes are needed.

###### Acknowledgements.

This research took advantage of the Simbad and Vizier databases, operated at
the CDS, Strasbourg (France), and NASA’s Astrophysics Data System
Bibliographic Services. We would like to acknowledge the Villum Foundation,
The Danish Council for Independent Research — Natural Science and the
Carlsberg Foundation for the support on building the SONG prototype on
Tenerife. K.B. acknowledges support from the Villum Foundation. The Stellar
Astrophysics Centre is funded by The Danish National Research Foundation
(Grant DNRF106) and research was supported by the ASTERISK project (ASTERo-
seismic Investigations with SONG and Kepler) funded by the European Research
Council (Grant agreement n. 267864). We also gratefully acknowledge the
support from the Spanish Ministry of Economy Competitiveness (MINECO) grant
AYA2016-76378-P.

## References

  * Andersen et al. (2016) Andersen, M. F., Grundahl, F., Beck, A. H., & Pallé, P. 2016, in Revista Mexicana de Astronomia y Astrofisica Conference Series, Vol. 48, Revista Mexicana de Astronomia y Astrofisica Conference Series, 54–58
  * Andersen et al. (2014) Andersen, M. F., Grundahl, F., Christensen-Dalsgaard, J., et al. 2014, in Revista Mexicana de Astronomia y Astrofisica, vol. 27, Vol. 45, Revista Mexicana de Astronomia y Astrofisica Conference Series, 83–
  * Angulo et al. (1999) Angulo, C., Arnould, M., Rayet, M., et al. 1999, Nuclear Physics A, 656, 3
  * Antoci et al. (2013) Antoci, V., Handler, G., Grundahl, F., et al. 2013, MNRAS, 435, 1563
  * Arentoft et al. (2017) Arentoft, T., Brogaard, K., Jessen-Hansen, J., et al. 2017, ApJ, 838, 115
  * Arentoft et al. (2008) Arentoft, T., Kjeldsen, H., Bedding, T. R., et al. 2008, ApJ, 687, 1180
  * Arentoft et al. (2014) Arentoft, T., Tingley, B., Christensen-Dalsgaard, J., et al. 2014, MNRAS, 437, 1318
  * Beck et al. (2015) Beck, P. G., Kambe, E., Hillen, M., et al. 2015, A&A, 573, A138
  * Beck et al. (2012) Beck, P. G., Montalban, J., Kallinger, T., et al. 2012, Nature, 481, 55
  * Bedding et al. (1996) Bedding, T. R., Kjeldsen, H., Reetz, J., & Barbuy, B. 1996, MNRAS, 280, 1155
  * Bedding et al. (2011) Bedding, T. R., Mosser, B., Huber, D., et al. 2011, Nature, 471, 608
  * Böhm-Vitense (1958) Böhm-Vitense, E. 1958, ZAp, 46, 108
  * Brogaard et al. (2016) Brogaard, K., Jessen-Hansen, J., Handberg, R., et al. 2016, Astronomische Nachrichten, 337, 793
  * Bubar & King (2010) Bubar, E. J. & King, J. R. 2010, AJ, 140, 293
  * Carigi (2000) Carigi, L. 2000, Rev. Mexicana Astron. Astrofis., 36, 171
  * Casagrande & VandenBerg (2014) Casagrande, L. & VandenBerg, D. A. 2014, MNRAS, 444, 392
  * Christensen-Dalsgaard (2008a) Christensen-Dalsgaard, J. 2008a, Ap&SS, 316, 113
  * Christensen-Dalsgaard (2008b) Christensen-Dalsgaard, J. 2008b, Ap&SS, 316, 13
  * Christensen-Dalsgaard (2014) Christensen-Dalsgaard, J. 2014, Asteroseismology of red giants, ed. P. L. Pallé & C. Esteban, 194
  * Corsaro & De Ridder (2014) Corsaro, E. & De Ridder, J. 2014, A&A, 571, A71
  * Corsaro et al. (2012) Corsaro, E., Grundahl, F., Leccia, S., et al. 2012, A&A, 537, A9
  * Feigelson & Babu (2012) Feigelson, E. & Babu, G. 2012, Modern Statistical Methods for Astronomy: With R Applications (Cambridge University Press)
  * Ferguson et al. (2005) Ferguson, J. W., Alexander, D. R., Allard, F., et al. 2005, ApJ, 623, 585
  * Frandsen et al. (1995) Frandsen, S., Jones, A., Kjeldsen, H., et al. 1995, A&A, 301, 123
  * Grundahl et al. (2017) Grundahl, F., Fredslund Andersen, M., Christensen-Dalsgaard, J., et al. 2017, ApJ, 836, 142
  * Guggenberger et al. (2016) Guggenberger, E., Hekker, S., Basu, S., & Bellinger, E. 2016, MNRAS, 460, 4277
  * Handberg et al. (2016a) Handberg, R., Brogaard, K. F., Miglio, A., et al. 2016a, submitted to MNRAS
  * Handberg & Campante (2011) Handberg, R. & Campante, T. L. 2011, A&A, 527, A56
  * Handberg et al. (2016b) Handberg, R., Miglio, A., Brogaard, K., Bossini, D., & Elsworth, Y. P. 2016b, Astronomische Nachrichten, 337, 799
  * Huber et al. (2010) Huber, D., Bedding, T. R., Stello, D., et al. 2010, ApJ, 723, 1607
  * Iglesias & Rogers (1996) Iglesias, C. A. & Rogers, F. J. 1996, ApJ, 464, 943
  * Jenkins et al. (2010) Jenkins, J. M., Caldwell, D. A., Chandrasekaran, H., et al. 2010, ApJ, 713, L120
  * Massarotti et al. (2008) Massarotti, A., Latham, D. W., Stefanik, R. P., & Fogel, J. 2008, AJ, 135, 209
  * Mazumdar et al. (2009) Mazumdar, A., Mérand, A., Demarque, P., et al. 2009, A&A, 503, 521
  * Miglio et al. (2012) Miglio, A., Brogaard, K., Stello, D., et al. 2012, MNRAS, 419, 2077
  * Mosser et al. (2013) Mosser, B., Michel, E., Belkacem, K., et al. 2013, A&A, 550, A126
  * Nordgren et al. (1999) Nordgren, T. E., Germain, M. E., Benson, J. A., et al. 1999, AJ, 118, 3032
  * Pietrinferni et al. (2004) Pietrinferni, A., Cassisi, S., Salaris, M., & Castelli, F. 2004, ApJ, 612, 168
  * Pinsonneault et al. (2014) Pinsonneault, M. H., Elsworth, Y., Epstein, C., et al. 2014, ApJS, 215, 19
  * Piskunov & Valenti (2002) Piskunov, N. E. & Valenti, J. A. 2002, A&A, 385, 1095
  * Rodrigues et al. (2017) Rodrigues, T. S., Bossini, D., Miglio, A., et al. 2017, MNRAS, 467, 1433
  * Rogers & Nayfonov (2002) Rogers, F. J. & Nayfonov, A. 2002, ApJ, 576, 1064
  * Skottfelt et al. (2015) Skottfelt, J., Bramich, D. M., Hundertmark, M., et al. 2015, A&A, 574, A54
  * Stello et al. (2006) Stello, D., Arentoft, T., Bedding, T. R., et al. 2006, MNRAS, 373, 1141
  * Stello et al. (2007) Stello, D., Bruntt, H., Kjeldsen, H., et al. 2007, MNRAS, 377, 584
  * Stello et al. (2009) Stello, D., Chaplin, W. J., Basu, S., Elsworth, Y., & Bedding, T. R. 2009, MNRAS, 400, L80
  * van Leeuwen (2007) van Leeuwen, F. 2007, A&A, 474, 653
  * Walker et al. (2003) Walker, G., Matthews, J., Kuschnig, R., et al. 2003, PASP, 115, 1023
  * Wenger et al. (2000) Wenger, M., Ochsenbein, F., Egret, D., et al. 2000, A&AS, 143, 9


# Mean Field Games Flock! The Reinforcement Learning Way

Sarah Perrin1 Mathieu Laurière2 Julien Pérolat3 Matthieu Geist4 Romuald Élie3
&Olivier Pietquin4 1Univ. Lille, CNRS, Inria, Centrale Lille, UMR 9189 CRIStAL
2Princeton University, ORFE
3DeepMind Paris
4Google Research, Brain Team<EMAIL_ADDRESS><EMAIL_ADDRESS>{perolat, mfgeist, relie<EMAIL_ADDRESS>
###### Abstract

We present a method enabling a large number of agents to learn how to flock,
which is a natural behavior observed in large populations of animals. This
problem has drawn a lot of interest but requires many structural assumptions
and is tractable only in small dimensions. We phrase this problem as a Mean
Field Game (MFG), where each individual chooses its acceleration depending on
the population behavior. Combining Deep Reinforcement Learning (RL) and
Normalizing Flows (NF), we obtain a tractable solution requiring only very
weak assumptions. Our algorithm finds a Nash Equilibrium and the agents adapt
their velocity to match the neighboring flock’s average one. We use Fictitious
Play and alternate: (1) computing an approximate best response with Deep RL,
and (2) estimating the next population distribution with NF. We show
numerically that our algorithm learn multi-group or high-dimensional flocking
with obstacles.

##  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 multi-group 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).

$\mathbb{E}_{z}$ 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 $x^{i}$ and $v^{i}$. We assume that
it can control its velocity by choosing the acceleration, denoted by $u^{i}$.
The dynamics of agent $i\in\left\\{1,\dots,N\right\\}$ is:

$x^{i}_{t+1}=x^{i}_{t}+v^{i}_{t}\Delta t,\qquad
v^{i}_{t+1}=v^{i}_{t}+u^{i}_{t}\Delta t+\epsilon^{i}_{t+1},$

where $\Delta t$ is the time step and $\epsilon^{i}_{t}$ 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:

$f^{i}_{t}=f(x^{i}_{t},v^{i}_{t},u^{i}_{t},\mu^{N}_{t}),$ (1)

where the interactions with other agents are only through the empirical
distribution of states and velocities denoted by:
$\mu^{N}_{t}=\frac{1}{N}\sum_{j=1}^{N}\delta_{(x^{j}_{t},v^{j}_{t})}.$

We focus on criteria incorporating a term of the form:

$\footnotesize
f^{\mathrm{flock}}_{\beta}(x,v,u,\mu)=\text{-}\left\|\int_{\mathbb{R}^{2d}}\frac{\left(v-v^{\prime}\right)}{(1+\left\|x-x^{\prime}\right\|^{2})^{\beta}}\,d\mu(x^{\prime},v^{\prime})\right\|^{2},$
(2)

where $\beta\geq 0$ is a parameter and
$(x,v,u,\mu)\in\mathbb{R}^{d}\times\mathbb{R}^{d}\times\mathbb{R}^{d}\times\mathcal{P}(\mathbb{R}^{d}\times\mathbb{R}^{d})$,
with $\mathcal{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 $\beta$ parameterizes the level of
interactions between agents and strongly impacts the flocking behavior: if
$\beta=0$, each agent tries to align its velocity with all the other agents of
the population irrespective of their positions, whereas the larger $\beta>0$,
the more importance is given to its closest neighbors (in terms of position).

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

$f^{\mathrm{flock},i}_{\beta,t}=-\left\|\frac{1}{N}\sum^{N}_{j=1}\frac{(v^{i}_{t}-v^{j}_{t})}{(1+\|x^{i}_{t}-x^{j}_{t}\|^{2})^{\beta}}\right\|^{2}.$
(3)

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}=(u^{1},\dots,u^{i-1},u^{i+1},\dots,u^{N})$ is:
$F^{i}_{u^{-i}}(u^{i})=\mathbb{E}_{x^{i}_{t},v^{i}_{t}}\Big{[}\sum_{t\geq
0}\gamma^{t}f^{i}_{t}\Big{]},$ with $f^{i}_{t}$ defined Eq. (1). In this
context, a _Nash equilibrium_ is a strategy profile
$(\hat{u}^{1},\hat{u}^{2},\dots,\hat{u}^{N})$ such that there’s no profitable
unilateral deviation, i.e., for every $i=1,\dots,N$, for every control
$u^{i}$, $F^{i}_{\hat{u}^{-i}}(\hat{u}^{i})\geq F^{i}_{\hat{u}^{-i}}(u^{i}).$

###  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 $\mu_{t}$ the state distribution of the
population, and by $\xi_{t}\in\mathbb{R}^{\ell}$ and
$\alpha_{t}\in\mathbb{R}^{k}$ the state and the control of an infinitesimal
agent, the dynamics of the infinitesimal agent is given by

$\xi_{t+1}=\xi_{t}+b(\xi_{t},\alpha_{t},\mu_{t})+\sigma\epsilon_{t+1},$ (4)

where
$b:\mathbb{R}^{\ell}\times\mathbb{R}^{k}\times\mathcal{P}(\mathbb{R}^{\ell})\to\mathbb{R}^{\ell}$
is a drift (or transition) function, $\sigma$ is a $\ell\times\ell$ matrix and
$\epsilon_{t+1}$ is a noise term taking values in $\mathbb{R}^{\ell}$. We
assume that the sequence of noises $(\epsilon_{t})_{t\geq 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 $\mu=(\mu_{t})_{t\geq
0}$ and a strategy $\alpha$ (i.e., a stochastic process adapted to the
filtration generated by $(\epsilon_{t})_{t\geq 0}$) as:
$J_{\mu}(\alpha)=\mathbb{E}_{\xi_{t},\alpha_{t}}\Big{[}\sum_{t\geq
0}\gamma^{t}\varphi(\xi_{t},\alpha_{t},\mu_{t})\Big{]},$ where
$\gamma\in(0,1)$ is a discount factor and
$\varphi:\mathbb{R}^{\ell}\times\mathbb{R}^{k}\times\mathcal{P}(\mathbb{R}^{\ell})\to\mathbb{R}^{\ell}$
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 $(\hat{\mu},\hat{\alpha})=(\hat{\mu}_{t},\hat{\alpha}_{t})_{t\geq 0}$
of a flow of distributions and strategies such that the following two
conditions are satisfied: $\hat{\alpha}$ is a best response against
$\hat{\mu}$ (optimality) and $\hat{\mu}$ is the distribution generated by
$\hat{\alpha}$ (consistency), i.e.,

  1. 1.

$\hat{\alpha}$ maximizes $\alpha\mapsto J_{\hat{\mu}}(\alpha)$;

  2. 2.

for every $t\geq 0$, $\hat{\mu}_{t}$ is the distribution of $\xi_{t}$ when it
follows the dynamics (4) with $(\alpha_{t},\mu_{t})$ replaced by
$(\hat{\alpha}_{t},\hat{\mu}_{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 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ère (2019);
Ruthotto et al. (2020). These techniques rely on the full knowledge of the
model and are restricted to classes of quite simple MFGs.

input : MFG = $\\{\xi,\varphi,\mu_{0}\\}$; # of iterations $J$ ;

1 Define $\bar{\mu}_{0}=\mu_{0}$ for  _$j=1,\dots,J$_ do

2 1. Set best response $\alpha_{j}=$
$\operatorname*{arg\,max}_{\alpha}J_{\bar{\mu}_{j-1}}(\alpha)$ and let
$\bar{\alpha}_{j}$ be the average of $(\alpha_{i})_{i=0,\dots,j}$

3 2. Set $\mu_{j}=$ $\gamma$-stationary distribution induced by $\alpha_{j}$

4 3. Set $\bar{\mu}_{j}=\frac{j-1}{j}\bar{\mu}_{j-1}+\frac{1}{j}\mu_{j}$

return _$\bar{\alpha}_{J},\bar{\mu}_{J}$_

Algorithm 1 Generic Fictitious Play in 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 $(\mathcal{S},\mathcal{A},p,r,\gamma)$ where
$\mathcal{S}$ is a state space, $\mathcal{A}$ is an action space,
$p:\mathcal{S}\times\mathcal{A}\to\mathcal{P}(\mathcal{S})$ is a transition
kernel, $r:\mathcal{S}\times\mathcal{A}\to\mathbb{R}$ is a reward function and
$\gamma$ 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
$\pi:\mathcal{S}\to\mathcal{P}(\mathcal{A})$, $s\mapsto\pi(\cdot|s)$ provides
a distribution over actions for each state. RL aims at learning a policy
$\pi^{*}$ which maximizes the total return defined as the expected
(discounted) sum of future rewards:

$R(\pi)=\mathbb{E}_{a_{t},s_{t+1}}\Big{[}\sum_{t\geq
0}\gamma^{t}r(s_{t},a_{t})\Big{]},$ (5)

with $a_{t}\sim\pi(\cdot|s_{t})$ and $s_{t+1}\sim p(\cdot|s_{t},a_{t})$. 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: $\pi\leftarrow\pi+\alpha\frac{\partial
R(\pi)}{\partial\pi},$ where $\alpha$ 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)\mapsto Q^{\pi}(s,a)$. $Q^{\pi}(s,a)$ is the
total return conditioned on starting in state $s$ and using action $a$ before
using policy $\pi$ 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 $\pi$ 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 an MFG

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
$\mu^{N}_{t}$ is expected to converge towards the law $\mu_{t}$ of
$(x_{t},v_{t})$, which represents the position and velocity of an
infinitesimal agent and have dynamics:

$x_{t+1}=x_{t}+v_{t}\Delta t,\qquad v_{t+1}=v_{t}+u_{t}\Delta
t+\epsilon_{t+1}.$

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 $\ell=2d$, $k=d$, and
$\xi=(x,v)$, $\alpha=u$. To accommodate for the degeneracy of the noise as
only the velocities are disturbed, we take
$\sigma=\begin{pmatrix}0_{d}&0_{d}\\\ 0_{d}&1_{d}\end{pmatrix}$ where $1_{d}$
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
$(\hat{\mu},\hat{u})$ where
$\hat{\mu}\in\mathcal{P}(\mathbb{R}^{\ell}\times\mathbb{R}^{\ell})$ is a
position-velocity distribution and
$\hat{u}:\mathbb{R}^{\ell}\times\mathbb{R}^{\ell}\mapsto\mathbb{R}^{\ell}$ is
a feedback function to determine the acceleration given the position and
velocity, such that: (1) $\hat{\mu}$ is an invariant position-velocity
distribution if the whole population uses the acceleration given by $\hat{u}$,
and (2) $\hat{u}$ maximizes the rewards when the agent’s initial position-
velocity distribution is $\hat{\mu}$ and the population distribution is
$\hat{\mu}$ at every time step. In mathematical terms, $\hat{u}$ maximizes
$J_{\hat{\mu}}(u)=\mathbb{E}_{x_{t},v_{t},u_{t}}\Big{[}\sum_{t\geq
0}\gamma^{t}\varphi(x_{t},v_{t},u_{t},\hat{\mu})\Big{]},$ where
$(x_{t},v_{t})_{t\geq 0}$ is the trajectory of an infinitesimal agent who
starts with distribution $\hat{\mu}$ at time $t=0$ and is controlled by the
acceleration $(u_{t})_{t\geq 0}$. As the payoff function $\varphi$ we use
$f^{\mathrm{flock}}_{\beta}$ from Eq. (2). Moreover, the consistency condition
rewrites as: $\hat{\mu}$ is the stationary distribution of
$(x_{t},v_{t})_{t\geq 0}$ if controlled by $(\hat{u}_{t})_{t\geq 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 Lasry-Lions 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 $\beta=0$, every
agent has the same influence over the representative agent and it is possible
to show that 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 $\frac{\sigma^{2}}{2}$. If
$\beta>0$, Nourian et al. (2011) shows that if the MF problem admits a unique
solution, then there exists an $\epsilon_{N}$ Nash equilibrium for the $N$
agents problem and $\lim\limits_{N\rightarrow+\infty}\epsilon_{N}=0$.
Existence has also been proved when $\beta>0$ in (Carmona and Delarue, 2018,
Section 4.7.3), with a slight modification of the payoff, namely considering,
with $\varphi$ a bounded function,
$r_{t}=-\varphi\Big{(}\Big{\|}\int_{\mathbb{R}^{2d}}\frac{\left(v_{t}-v^{\prime}\right)}{\left(1+\left\|x_{t}-x^{\prime}\right\|^{2}\right)^{\beta}}\mu_{t}(dx^{\prime},dv^{\prime})\Big{\|}^{2}\Big{)}$.

###  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 $\gamma$-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
$\bar{\mu}_{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 $\bar{\mu}_{j}$.

input : MFG = $\\{(x,v),f^{\mathrm{flock}}_{\beta},\mu_{0}\\}$; # of
iterations $J$

1 Define $\bar{\mu}_{0}=\mu_{0}$ for  _$j=1,\dots,J$_ do

2 1. Set best response
$\displaystyle\pi_{j}=\operatorname*{arg\,max}_{\pi}J_{\bar{\mu}_{j-1}}(\pi)$
with SAC and let $\bar{\pi}_{j}$ be the average of $(\pi_{0},\dots,\pi_{j})$

3 2. Using a Normalizing Flow, compute $\mu_{j}=$ $\gamma$-stationary
distribution induced by $\pi_{j}$

4 3. Using a Normalizing Flow and samples from $(\mu_{0},\dots,\mu_{j-1})$,
estimate $\bar{\mu}_{j}$

return _$\bar{\pi}_{J},\bar{\mu}_{J}$_

Algorithm 2 Flock’n RL

####  3.2.1 Computing the Best Response with SAC

The first step in the loop of Alg. 2 is the computation of a best response
against $\bar{\mu}_{j}$. In fact, the problem boils down to solving an MDP in
which $\bar{\mu}_{j}$ enters as a parameter. Following the notation introduced
in Sec. 2.3, we take the state and action spaces to be respectively
$\mathcal{S}=\mathbb{R}^{2d}$ (for position-velocity pairs) and
$\mathcal{A}=\mathbb{R}^{d}$ (for accelerations). Letting $s=(x,v)$ and $a=u$,
the reward is: $(x,v,u)=(s,a)\mapsto r(s,a)=f(x,v,u,\bar{\mu}_{j}),$ which
depends on the given distribution $\bar{\mu}_{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 $\bar{\mu}_{j}$, the fixed average
distribution at step $j$ of Flock’n RL. SAC maximizes the reward which is a
variant of $f^{\mathrm{flock},i}_{\beta,t}$ from Eq. (3). It needs samples
from $\bar{\mu}_{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 $\bar{\mu}_{j}$ at
the beginning of step 1 (i.e. we do not sample new agents from $\bar{\mu}_{j}$
every time we need to compute the reward). During the learning, at the
beginning of each episode, we sample a starting state
$s_{0}\sim\bar{\mu}_{j}$.

In the experiments, we will not need $\bar{\pi}_{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 $(\pi_{0},\dots,\pi_{j})$ and simply average
the induced rewards.

####  3.2.2 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\sim 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))\left|\text{det}\left(\frac{\partial
f^{-1}}{\partial x}\right)\right|.$ 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 $\sum_{n}\log
p(x^{(n)})$.

####  3.2.3 Computation of $\bar{\mu}_{j}$

Due to the above discussion on the difficulty to represent the distribution in
continuous space and high dimension, the third step (Line 2 of Alg. 2) can not
be implemented easily. We represent every $\mu_{j}$ as a generative model, so
we can not “average” the normalizing flows corresponding to
$(\mu_{i})_{i=1,\dots,j}$ in a straightforward way but we can sample data
points $x\sim\mu_{i}$ for each $i=1,\dots,j$. To have access to
$\bar{\mu}_{j}$, we keep in memory every model $\mu_{j}$,
$j\in\left\\{1,\dots,J\right\\}$ and, in order to sample points according to
$\bar{\mu}_{j}$ for a fixed $j$, we sample points from
$\mu_{i},i\in\\{1,\dots,j\\},$ with probability $1/j$. These points are then
used to learn the distribution $\bar{\mu}_{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\in S$
as $s=(x,v)$ where $x$ and $v$ are respectively the vectors of positions and
velocities. Each coordinate $x_{i}$ of the position can take any continuous
value in the $d$-dimensional box $x_{i}\in[-100,+100]$, while the velocities
are also continuous and clipped $v_{i}\in[-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 $\bar{\mu}_{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
$\mu$, the objective is to maximize:
$J_{\mu}(\pi)=\mathbb{E}_{(s_{t},u_{t})}\left[\sum_{t=0}^{+\infty}\gamma^{t}r(x_{t},v_{t},u_{t},\mu_{t})+\delta\mathcal{H}(\pi(\cdot|s_{t}))\right],$
where $\mathcal{H}$ denotes the entropy and $\delta\geq 0$ is a weight.

To implement the optimization, the SAC algorithm follows the philosophy of
actor-critic by training parameterized $Q$-function 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 $\pi$ with a best response $\pi^{\prime}$, when
the rest of the population plays $\pi$:
$\phi(\pi):=\max\limits_{\pi^{\prime}}J(\mu_{0},\pi^{\prime},\mu^{\pi})-J(\mu_{0},\pi,\mu^{\pi})$.
If $\phi(\bar{\pi}_{j})\to 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 $\mathcal{M}$ of performance of
learned policies versus estimated distributions. The entry $\mathcal{M}_{i,j}$
on the $i$-th row and the $j$-th column is the total $\gamma$-discounted sum
of rewards:
$\mathcal{M}_{i,j}=\mathbb{E}\left[\sum_{t=0}^{T}\gamma^{t}r_{t,i}\;|\;s_{0}\sim\bar{\mu}_{i-1},\;u_{t}\sim\pi_{j}(.|s_{t})\right],$
where $r_{t,i}=r(s_{t},u_{t},\bar{\mu}_{i-1})$, obtained with $\pi_{j}$
against $\bar{\mu}_{i-1}$. The diagonal term $\mathcal{M}_{j,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 $\phi(\bar{\pi})$ directly in the Flock’n RL
algorithm with SAC. The second term, $J(\mu_{0},\bar{\pi},\mu^{\bar{\pi}})$,
can be approximated by replacing $\bar{\pi}$ with the average over past
policies, i.e., the policy sampled uniformly from the set
$\\{\pi_{0},\dots,\pi_{j}\\}$. At step $j$, the approximate exploitability is
$e_{j}=\mathcal{M}_{j,j}-\frac{1}{j-1}\sum_{k=1}^{j-1}\mathcal{M}_{j,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 $\beta$ in the
flocking effect. We consider noise $\epsilon^{i}_{t}\sim\mathcal{N}(0,\Delta
t)$ and the following reward:
$r^{i}_{t}=f^{\mathrm{flock},i}_{\beta,t}-\|u^{i}_{t}\|_{2}^{2}+\|v^{i}_{t}\|_{\infty}-\min\\{\|x^{i}_{2,t}\pm
50\|\\},$ where $x^{i}_{2,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
$\|v^{i}_{t}\|_{\infty}=\max\\{|v^{i}_{1,t}|,|v^{i}_{2,t}|\\}$. Hence, a
possible equilibrium is with two groups of agents, one for each line. When
$\beta=0$, the term $f^{\mathrm{flock},i}_{\beta,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 $\beta>0$ is large enough (e.g. $\beta=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 $\beta=0$, we experimentally verified that there is always a
global consensus, i.e., only one line or two lines but moving in the same
direction.

(a) Initial positions and velocities

(b) At convergence

(c) Performance matrix

(d) Approximate exploitability

Figure 1: Multi-group flocking with noise and $\beta=100$.

##### 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:
$r^{i}_{t}=f^{\mathrm{flock},i}_{\beta,t}-\|u^{i}_{t}\|_{2}^{2}+\|v^{i}_{t}\|_{\infty}-\min\\{\|x^{i}_{2,t}\|\\}-c*\mathbf{1}_{obs},$
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.

(a) Initial positions and velocities

(b) At convergence

(c) Performance matrix

(d) Approximate exploitability

Figure 2: Flocking with noise and many obstacles.

##  5 Related Work

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 in MFGs.  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)
combined 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 finite-state 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.

## References

  * Achdou and Capuzzo-Dolcetta [2010] Yves Achdou and Italo Capuzzo-Dolcetta.  Mean field games: numerical methods.  SIAM J. Numer. Anal., 2010.
  * Achdou et al. [2020] Yves Achdou, Paola Mannucci, Claudio Marchi, and Nicoletta Tchou.  Deterministic mean field games with control on the acceleration.  NoDEA, 2020.
  * Bailo et al. [2018] Rafael Bailo, Mattia Bongini, José A Carrillo, and Dante Kalise.  Optimal consensus control of the Cucker-Smale model.  IFAC, 2018.
  * Caponigro et al. [2013] Marco Caponigro, Massimo Fornasier, Benedetto Piccoli, and Emmanuel Trélat.  Sparse stabilization and optimal control of the Cucker-Smale model.  Mathematical Control and Related Fields, 2013.
  * Cardaliaguet and Hadikhanloo [2017] Pierre Cardaliaguet and Saeed Hadikhanloo.  Learning in mean field games: the fictitious play.  ESAIM Cont. Optim. Calc. Var., 2017.
  * Carlini and Silva [2014] Elisabetta Carlini and Francisco J. Silva.  A fully discrete semi-Lagrangian scheme for a first order mean field game problem.  SIAM J. Numer. Anal., 2014.
  * Carmona and Delarue [2018] René Carmona and François Delarue.  Probabilistic theory of mean field games with applications. I.  2018\.
  * Carmona and Laurière [2019] René Carmona and Mathieu Laurière.  Convergence analysis of machine learning algorithms for the numerical solution of mean field control and games: II - the finite horizon case.  2019\.
  * Cucker and Smale [2007] Felipe Cucker and Steve Smale.  Emergent behavior in flocks.  IEEE Transactions on automatic control, 2007.
  * Durkan et al. [2019] Conor Durkan, Artur Bekasov, Iain Murray, and George Papamakarios.  Neural spline flows, 2019.
  * Elie et al. [2020] Romuald Elie, Julien Perolat, Mathieu Laurière, Matthieu Geist, and Olivier Pietquin.  On the convergence of model free learning in mean field games.  In proc. of AAAI, 2020.
  * Guo et al. [2019] Xin Guo, Anran Hu, Renyuan Xu, and Junzi Zhang.  Learning mean-field games.  In proc. of NeurIPS, 2019.
  * Haarnoja et al. [2018] Tuomas Haarnoja, Aurick Zhou, Pieter Abbeel, and Sergey Levine.  Soft actor-critic: Off-policy maximum entropy deep reinforcement learning with a stochastic actor.  CoRR, abs/1801.01290, 2018.
  * Hill et al. [2018] Ashley Hill, Antonin Raffin, Maximilian Ernestus, Adam Gleave, Anssi Kanervisto, Rene Traore, Prafulla Dhariwal, Christopher Hesse, Oleg Klimov, Alex Nichol, Matthias Plappert, Alec Radford, John Schulman, Szymon Sidor, and Yuhuai Wu.  Stable baselines.  https://github.com/hill-a/stable-baselines, 2018.
  * Kobyzev et al. [2020] Ivan Kobyzev, Simon Prince, and Marcus Brubaker.  Normalizing flows: An introduction and review of current methods.  PAMI, 2020.
  * Lanctot et al. [2017] Marc Lanctot, Vinicius Zambaldi, Audrunas Gruslys, Angeliki Lazaridou, Karl Tuyls, Julien Perolat, David Silver, and Thore Graepel.  A unified game-theoretic approach to multiagent reinforcement learning.  In proc. of NeurIPS, 2017.
  * Lasry and Lions [2007] Jean-Michel Lasry and Pierre-Louis Lions.  Mean field games.  Jpn. J. Math., 2007.
  * Lillicrap et al. [2016] Timothy P Lillicrap, Jonathan J Hunt, Alexander Pritzel, Nicolas Heess, Tom Erez, Yuval Tassa, David Silver, and Daan Wierstra.  Continuous control with deep reinforcement learning.  In proc. of ICLR, 2016.
  * Nourian et al. [2011] Mojtaba Nourian, Peter E Caines, and Roland P Malhamé.  Mean field analysis of controlled Cucker-Smale type flocking: Linear analysis and perturbation equations.  IFAC, 2011.
  * Perolat et al. [2018] Julien Perolat, Bilal Piot, and Olivier Pietquin.  Actor-critic fictitious play in simultaneous move multistage games.  In Proc. of AISTATS, 2018.
  * Perrin et al. [2020] Sarah Perrin, Julien Pérolat, Mathieu Laurière, Matthieu Geist, Romuald Elie, and Olivier Pietquin.  Fictitious play for mean field games: Continuous time analysis and applications.  In proc. of NeurIPS, 2020.
  * Reynolds [1987] Craig W Reynolds.  Flocks, herds and schools: A distributed behavioral model.  In Proc. of SIGGRAPH, 1987.
  * Rezende and Mohamed [2015] Danilo Rezende and Shakir Mohamed.  Variational inference with normalizing flows.  In proc. of ICML, 2015.
  * Ruthotto et al. [2020] Lars Ruthotto, Stanley J Osher, Wuchen Li, Levon Nurbekyan, and Samy Wu Fung.  A machine learning framework for solving high-dimensional mean field game and mean field control problems.  PNAS, 2020.
  * Schulman et al. [2015] John Schulman, Sergey Levine, Philipp Moritz, Michael I. Jordan, and Pieter Abbeel.  Trust region policy optimization.  CoRR, abs/1502.05477, 2015.
  * Shaw [1975] E Shaw.  Naturalist at large-fish in schools.  Natural History, 1975.
  * Silver et al. [2018] David Silver, Thomas Hubert, Julian Schrittwieser, Ioannis Antonoglou, Matthew Lai, Arthur Guez, Marc Lanctot, Laurent Sifre, Dharshan Kumaran, Thore Graepel, Timothy Lillicrap, Karen Simonyan, and Demis Hassabis.  A general reinforcement learning algorithm that masters chess, shogi, and Go through self-play.  Science, 2018.
  * Subramanian and Mahajan [2019] Jayakumar Subramanian and Aditya Mahajan.  Reinforcement learning in stationary mean-field games.  In proc. of AAMAS, 2019.
  * Sutton et al. [1999] R. S. Sutton, D. Mcallester, S. Singh, and Y. Mansour.  Policy gradient methods for reinforcement learning with function approximation.  In proc. of NeurIPS. MIT Press, 1999.

##  Appendix A More numerical tests

###  A.1 A Simple Example in Four dimensions

We illustrate in a simple four dimensional setting (i.e. two-dimensional
positions and velocities hence the total dimension is $4$) how the agents
learn to adopt similar velocities by controlling their acceleration.

Here, we take $\epsilon^{i}_{t}\equiv 0$ (no noise), and we define the reward
as:

$r^{i}_{t}=f^{\mathrm{flock},i}_{\beta=0,t}-\|u^{i}_{t}\|_{2}^{2}+\|v^{i}_{t}\|_{2}^{2}.$
(6)

We set $\beta=0$ to have a intuitive example. The first term encourages the
agent to adapt its velocity to the crowd’s one, giving equal importance to all
the agents irrespective of their distance. The second term penalizes a strong
acceleration and we added the last term (not present in the original flocking
model) to reduce the number of Nash equilibria and prevent the agent to
converge to a degenerate solution which consists in putting the whole crowd at
a common position with a null velocity.

As the velocity is bounded by 1 in our experiments, there are at least four
obvious Nash equilibria in terms of velocity that remain from reward (6):
$\hat{v}^{i}_{t}\equiv(-1,-1)$, $(-1,1)$, $(1,-1)$ and $(1,1)$, while the
position is not important anymore and any distribution for the positions is
valid. Experimentally, we observe that if we start from a normal initial
distribution $v^{i}_{0}\sim\mathcal{N}(0,1)$, then the equilibrium found by
Flock’n RL is randomly one of the four previous velocities. However, if we set
the initial velocities with a positive or negative bias, then we observe
experimentally that the corresponding equilibrium is reached. Thus, as we
could expect, the Nash equilibria found by the algorithm depends on the
initial distribution.

In Fig. 3, we can see that the agents have adopted the same velocities and the
performance matrix indicates that the policy learned at each step of Flock’n
RL tends to perform better than previous policies. The exploitablity decreases
quickly because Flock’n RL algorithm learns during the first iterations an
approximation of a Nash equilibrium.

(a) Initial positions and velocities

(b) At convergence

(c) Performance matrix

(d) Approximate exploitability

Figure 3: Flocking with an intuitive example

###  A.2 A Simple Example in Six Dimensions

We consider the simple example defined with reward from Eq. (6) but now we add
an additional dimension for the position and the velocity, making the problem
six-dimensional. As before, the agents are still encouraged to maximize their
velocities. We set $\beta=0$ and we do not put any noise on the velocity
dynamics. We can notice that experimentally a consensus is reached by the
agents as they learn to adopt the same velocities (Fig. 4(b)), even if the
agents start from a random distribution in terms of positions and velocities
(Fig. 4(a)). The performance matrix (Fig. 4(c)) highlights that the best
response improves until iteration 40, which explains why there is a bump in
performance in the exploitability before the 40th iteration.

(a) Initial positions and velocities

(b) At convergence

(c) Performance matrix

(d) Approximate exploitability

Figure 4: Flocking for a simple example in dimension 6.

###  A.3 Examples with an obstacle

We present an example with a single obstacle and noise in the dynamics, both
in dimension 4 and 6. In dimension 4, the obstacle is a square located at the
middle of the environment, whereas in dimension 6 it is a cube. In dimension
4, we can see in Fig. 5 that the agents that are initially spawned in the
environment with random positions and velocities manage to learn to adopt the
same velocities while avoiding the obstacle (Fig. 5(b)). The same behavior is
observed in dimension 6. We also notice that the exploitability is slower to
decrease in dimension 6, similarly to what we observed with the simple example
without the obstacle.

(a) Initial positions and velocities

(b) At convergence

(c) Performance matrix

(d) Approximate exploitability

Figure 5: Flocking in 4D with one obstacle.

(a) Initial positions and velocities

(b) At convergence

(c) Performance matrix

(d) Approximate exploitability

Figure 6: Flocking in 6D with one obstacle.

###  A.4 Many obstacles in 4D

Finally, we present an example in 4D with many obstacles, in the same fashion
than the 6-dimensional example with the columns located in the main part of
the article.

(a) Initial positions and velocities

(b) At convergence

(c) Performance matrix

(d) Approximate exploitability

Figure 7: Flocking in 4D with many obstacles.

##  Appendix B Normalizing Flows

In this section we provide some background on Normalizing Flows, which are an
important part of our approach.

##### Coupling layers

A coupling layer applies a coupling transform, which maps an input $x$ to an
output $y$ by first splitting $x$ into two parts $x=[x_{1:d-1},x_{d:D}]$,
computing parameters $\theta=NN(x_{1:d-1})$ of an arbitrary neural network (in
our case a fully connected neural network with 8 hidden units), applying $g$
to the last coordinates $y_{i}=g_{\theta_{i}}(x_{i})$ where $i\in[d,\dots,D]$
and $g_{\theta_{i}}$ is an invertible function parameterized by $\theta_{i}$.
Finally, we set $y_{1:d-1}=x_{1:d-1}$. Coupling transforms offer the benefit
of having a tractable Jacobian determinant, and they can be inverted exactly
in a single pass.

##### Neural Spline Flows

NSFs are based on monotonic rational-quadratic splines. These splines are used
to model the function $g_{\theta_{i}}$. A rational-quadratic function takes
the form of a quotient of two quadratic polynomials, and a spline uses $K$
different rational-quadratic functions.

Following the implementation described in Durkan et al. [2019], we detail how
a NSF is computed.

  1. 1.

A neural network NN takes $x_{1:d-1}$ as inputs and outputs $\theta_{i}$ of
length $3K-1$ for each $i\in[1,\dots,D]$.

  2. 2.

$\theta_{i}$ is partitioned as
$\theta_{i}=[\theta_{i}^{w},\theta_{i}^{h},\theta_{i}^{d}]$, of respective of
sizes $K,K,$ and $K-1$.

  3. 3.

$\theta_{i}^{w}$ and $\theta_{i}^{h}$ are passed through a softmax and
multiplied by $2B$, the outputs are interpreted as the widths and heights of
the $K$ bins. Cumulative sums of the $K$ bin widths and heights yields the
$K+1$ knots ${(x^{k},y^{k})}^{K}_{k=0}$.

  4. 4.

$\theta_{i}^{d}$ is passed through a softplus function and is interpreted as
the values of the derivatives of the internal knots.

##  Appendix C Visual Rendering with Unity

Once we have trained the policy with the Flock’n RL algorithm, we generate
trajectories of many agents and stock them in a csv file. We have coded an
integration in Unity, making it possible to load these trajectories and
visualize the flock in movement, interacting with its environment. We can then
easily load prefab models of fishes, birds, or any animal that we want our
agents to be. We can also load the obstacles and assign them any texture we
want. Examples of rendering are available in Fig. 8.

Figure 8: Visual rendering with Unity