23-358

# Cislunar Satellite Constellation Design via Integer Linear Programming

Malav Patel Graduate Student, Aerospace Engineering, Georgia Tech Yuri Shimane
PhD Candidate, Aerospace Engineering, Georgia Tech Hang Woon Lee Assistant
Professor, Mechanical and Aerospace Engineering, West Virginia University and
Koki Ho Professor, Aerospace Engineering, Georgia Tech

###### Abstract

Cislunar space awareness is of increasing interest to the international
community as Earth-Moon traffic is projected to increase. This raises the
problem of placing satellites optimally in a constellation to provide
satisfactory coverage for said traffic. The Circular Restricted 3 Body Problem
(CR3BP) provides promising periodic orbits in the Earth-Moon rotating frame
for traffic monitoring. This work converts a spatially and temporally varying
traffic coverage requirement into an integer linear programming problem,
attempting to minimize the number of satellites required for the requested
coverage.

##  1 Introduction

The Earth-Moon system and related orbital trajectories have been of interest
since the early Apollo missions in the 1960s and 1970s. This trend carried
into the early part of the 21st century with research efforts aiming to find
better trajectories for tours in the Earth-Moon system and operations near the
Moon.[1],[2]

The Artemis program has renewed this interest as Earth-Moon traffic is
projected to increase through crewed and uncrewed missions in the coming
decade. This makes the region of space between the Earth and the Moon, coined
cislunar space, an attractive region to monitor. Ground based measurement is
unable to provide sufficient observational power for such a vast region. This
insufficiency motivates space-based monitoring schemes, realized as satellite
constellations with optical payloads on favorable orbits for surveillance.
Technological advancements over the last several decades have equipped
satellites with stronger optical payloads, shrinking the performance gap
between space-based and ground based observation technology. Several
works[3],[4] have begun to investigate optimal architectures for surveillance.
Intuitively, the cost of launching and maintaining a patrol constellation
scales with the number of satellites in the constellation.

Following this demand for surveillance capability, Wilmer et al.[5] and Dahlke
et al.[6] investigated using periodic orbits to monitor activity near the
Lagrange points (L1-L5) in the Earth-Moon system. Zimovan et al.[7]
investigated near-rectilinear halo orbits (NRHOs) around the Moon because of
attractive stability, transfer prospectives, and eclipsing properties. Because
of this, NRHOs are promising orbits for lunar space stations that serve as
outposts for cislunar activity. Klonowski et al.[8] formulated space
situational awareness (SSA) architecture design as a Markov Decision Process
and employed the popular Monte Carlo Tree Search to optimize a constellation
of up to 5 observers to maintain custody of a space object traveling on a
specific orbit between the Earth and Moon. Vendl et al.[9] explored a more
general region of cislunar space, specifically a 20 degree cone emanating from
the Earth towards the Moon. The region consisted of 1071 target points, which
were used to choose an optimal phase of an observer on an L1 (and separately,
L2) Lyapunov orbit to maximize observation capability. However, the number of
satellites was not a design variable in this formulation. In contrast,
Visonneau et al. [10] presented a genetic algorithm for finding optimally
sized constellations for monitoring regions of cislunar space, with no
restriction on number of observers.

Monitoring capability depends on favorable SSA architectures, namely an
attractive choice of orbits to place observers in. Frueh et al. [11]
introduces the 2:1 resonant orbit as an attractive candidate for monitoring
cislunar space because of its close approaches to the Earth and Moon and
because it covers the entire cislunar region in under 20 revolutions. Gupta et
al. [12] investigated how phasing an observer and a data relay satellite on
separate orbits can impact the length of the window in which communication is
possible. This phasing, however, is for the tail of the observation process,
where already collected data needs to be relayed back to Earth.

This work aims to optimize the phasing and number of satellites to track and
observe a specific target moving in the cislunar region.

##  2 Problem Overview

To motivate the solution procedure, we present the results of past work in a
similar domain. Lee et al. [13] presented spatio-temporally varying Earth
observation as an integer linear programming problem (ILP). The formulation
sought to minimize the number of satellites required to satisfy a given
regional coverage requirement on Earth’s surface. This work’s main
contribution is to show that we can extend the formulation presented by Lee et
al. [13] to the cislunar regime.

###  2.1 Accessibility Profile

Suppose there is target point in space between the Earth and Moon, denoted by
$\boldsymbol{p}$ and an observer at position $\boldsymbol{r}$. The target
point’s observability can be quantified by an accessibility measure, $g$. A
simple measure could be the distance between observer and target,
$g(\boldsymbol{p},\boldsymbol{r})=\|\boldsymbol{p}-\boldsymbol{r}\|$. A more
sophisticated measure is the apparent magnitude of the object, which is a
function of the position of the sun as well,
$g=g(\boldsymbol{p},\boldsymbol{r},\boldsymbol{s})$ [14].

Suppose the observer is on an orbit with period $T$, time is discretized for
computational ease, and $\boldsymbol{p}$ is a static target point in space. At
time step $n$, $\boldsymbol{r}_{n}$ is the observer’s position,
$\boldsymbol{s}_{n}$ the sun’s position, and
$a[n]=g(\boldsymbol{p},\boldsymbol{r}_{n},\boldsymbol{s}_{n})$ is the
accessibility measure. If the orbit is discretized into $L$ equal time steps,
then $a$ is a $L\times 1$ vector containing the accessibility measures along
the observer’s orbit. In practice, the target point will be in one of two
states, [observable] or [not observable]. In other words, the continuous
accessibility measure must be transformed into a binary one. To achieve this,
we introduce the accessibility profile, $v$, distinct from the accessibility
measure vector, $a$. Define a threshold value, $M$, for the accessibility
measure. Then,

$v[n]=\left\\{\begin{array}[]{lr}1,&\text{if }a[n]\leq M\\\
0,&\text{otherwise}\end{array}\right\\}$

At time step $n$, $v[n]$ indicates whether the observer can view the target
point.

###  2.2 Constellation Pattern Vector

Instead of a single observer, suppose we have a set (called a constellation)
of observers on the same orbit of period $T$ again discretized into L time
steps. The constellation pattern vector $x$ captures the phasing of the
observers in the constellation relative to a seed satellite. At time step
$n=0$, if the seed satellite is at position $\boldsymbol{r}$, then after
$n_{k}$ time steps, the $k$-th observer will be at position $\boldsymbol{r}$.
The continuous-time time delay for the $k$-th observer is $\Delta
t_{k}=n_{k}*\frac{T}{L}$. This is illustrated below in Figure 1.

Figure 1: Correspondence between constellation pattern vector and
constellation setup[13]

###  2.3 Coverage Timeline and Coverage Requirement

The accessibility profile is a vector of binary variables describing when the
target is in view of the observer as the observer travels along its orbit. The
coverage timeline, denoted by $b$, is a $L\times 1$ vector that extends the
accessibility profile concept to multiple observers. At time step $n$, $b[n]$
describes the number of satellites with view access to the target point.
Instead of a vector of binary variables, the coverage timeline is a vector of
non-negative integers.

It turns out that a circular convolution between the constellation pattern
vector $x$ and the accessibility profile $v$ returns the coverage timeline
$b$. The circular nature of the convolution is a direct consequence of the
periodicity of the observers’ orbit. As such, $b$ can be computed from $x$ and
$v$:

$\displaystyle b[n]$ $\displaystyle=\sum_{m=0}^{L-1}v[m]*x[(n-m)\bmod L]$
$\displaystyle=\sum_{m=0}^{L-1}x[m]*v[(n-m)\bmod L]$

Alternatively, the circular convolution can be conveniently expressed as a
matrix multiplication between an $L\times L$ circulant matrix, $V$, and the
constellation pattern vector, $x$. Here, the $i$-th column of $V$ is produced
by circularly shifting $v$ by $i$ elements.

$Vx=b$

Expanded in matrix form and using 1-based indexing,

$\displaystyle{\begin{bmatrix}v[1]&v[L]&\cdots&v[3]&v[2]\\\
v[2]&v[1]&v[L]&&v[3]\\\ \vdots&v[2]&v[1]&\ddots&\vdots\\\
v[L-1]&&\ddots&\ddots&v[L]\\\ v[L]&v[L-1]&\cdots&v[2]&v[1]\\\
\end{bmatrix}}{\begin{bmatrix}x[1]\\\ x[2]\\\ \vdots\\\ \\\
x[{L}]\end{bmatrix}}={\begin{bmatrix}b[1]\\\ b[2]\\\ \vdots\\\ \\\
\b{[}{L}]\end{bmatrix}}$

Often times, we are interested in satisfying a coverage requirement, $f$,
which demands a different number of observers at different time steps. The
goal, then, is to place a set of observers on the orbit (i.e. find a
constellation pattern $x$) such that the resulting coverage timeline $b$ meets
or exceeds the coverage requirement $f$. Stated mathematically,

Find $\displaystyle x$ such that $\displaystyle b[n]\geq f[n]$
$\displaystyle\text{for}\ n=\\{1,2,\cdots,L\\}$ where $\displaystyle Vx=b$

###  2.4 Integer Linear Program Formulation

In general, there can be many constellation pattern vectors that satisfy the
coverage requirement. However, it is often beneficial to minimize the total
number of observers in the constellation while still satisfying the coverage
requirement, since fewer observers reduce launch and maintenance costs.

The problem above of finding a constellation pattern $x$ is now a problem of
finding the optimal constellation pattern $x^{*}$ that minimizes the number of
observers in the constellation. The following integer linear program presents
this mathematically.

$\displaystyle\underset{x}{\text{minimize}}$
$\displaystyle\boldsymbol{1}^{T}x$ subject to $\displaystyle Vx\geq f,$
$\displaystyle x\in\\{0,1\\}^{L}$

The objective is to minimize the total number of satellites, hence the
objective $\boldsymbol{1}^{T}x$, which sums all the entries in $x$.

###  2.5 Extension to Multiple Orbits & Target Points

The above formulation is linear, and thus allows extensions to multiple target
points and multiple constellations.

####  2.5.1 Multiple Target Points

Suppose we have a set $\mathcal{J}$ of target points and we seek to find an
optimal constellation of observers on a single orbit to monitor these target
points according to a set of coverage requirements
$\mathcal{F}=\\{f_{1},f_{2},\cdots,f_{\lvert\mathcal{J}\rvert}\\}$. Each
target point in $\mathcal{J}$ will have its own accessibility profile, and
thus its own circulant matrix. Denote the circulant matrix for each target
point by $V_{j,1}$, where $j$ indexes the elements in $\mathcal{J}$. Then we
have

$\displaystyle\underset{x}{\text{minimize}}$
$\displaystyle\boldsymbol{1}^{T}x$ subject to $\displaystyle V_{1,1}\ x\geq
f_{1},$ $\displaystyle V_{2,1}\ x\geq f_{2},$ $\displaystyle\vdots$
$\displaystyle V_{\lvert\mathcal{J}\rvert,1}\ x\geq
f_{\lvert\mathcal{J}\rvert}$ $\displaystyle x\in\\{0,1\\}^{L}$

The linear constraint can be condensed into a single matrix equation by
vertically concatenating the circulant matrices $V_{j,1}$ and coverage
requirements $f_{j}$

$\displaystyle\underset{x}{\text{minimize}}$
$\displaystyle\boldsymbol{1}^{T}x$ subject to
$\displaystyle{\begin{bmatrix}V_{1,1}\\\ V_{2,1}\\\ \vdots\\\
V_{\lvert\mathcal{J}\rvert,1}\end{bmatrix}}x\geq{\begin{bmatrix}f_{1}\\\
f_{2}\\\ \vdots\\\ f_{\lvert\mathcal{J}\rvert}\end{bmatrix}}$ $\displaystyle
x\in\\{0,1\\}^{L}$

####  2.5.2 Multiple Orbits

Suppose we have a set $\mathcal{Z}$ of orbits and we seek to an optimal
constellation of observers on each orbit to monitor a single target point.
Each orbit in $\mathcal{Z}$ will have its own accessibility profile, and thus
its own circulant matrix. Denote the circulant matrix for each orbit by
$V_{1,z}$ where $z$ indexes the elements in $\mathcal{Z}$. Furthermore, each
orbit will have its own constellation pattern vector, $x_{z}$. Then we have

$\displaystyle\underset{x}{\text{minimize}}$
$\displaystyle\sum_{z=1}^{\lvert\mathcal{Z}\rvert}\boldsymbol{1}^{T}x_{z}$
subject to $\displaystyle V_{1,1}\ x_{1}\geq f,\ V_{1,2}\ x_{2}\geq f,\
\cdots,V_{1,\lvert\mathcal{Z}\rvert}\ x_{\lvert\mathcal{Z}\rvert}\geq f$
$\displaystyle x\in\\{0,1\\}^{L}$

The linear constraint can be condensed into a single matrix equation by
horizontally concatenating the circulant matrices $V_{1,z}$ and vertically
concatenating the constellation pattern vectors $x_{z}$.

$\displaystyle\underset{\boldsymbol{x}}{\text{minimize}}$
$\displaystyle\sum_{z=1}^{\lvert\mathcal{Z}\rvert}\boldsymbol{1}^{T}x_{z}$
subject to
$\displaystyle{\begin{bmatrix}V_{1,1}&V_{1,2}&\cdots&V_{1,\lvert\mathcal{Z}\rvert}\end{bmatrix}}{\begin{bmatrix}x_{1}\\\
x_{2}\\\ \vdots\\\ x_{\lvert\mathcal{Z}\rvert}\end{bmatrix}}\geq f$
$\displaystyle x_{z}\in\\{0,1\\}^{L}$

We can combine the linearities above to provide a single comprehensive integer
linear program. Let $\boldsymbol{x}$ be the augmented constellation pattern
vector obtained by vertically concatenating the constellation pattern vectors
$x_{z}$, where $z$ indexes the elements in $\mathcal{Z}$. Let $\boldsymbol{f}$
be the augmented coverage requirement obtained by vertically concatenating the
coverage requirements $f_{j}$ for each target point, where $j$ indexes the
elements in $\mathcal{J}$. Let $\boldsymbol{V}$ be the augmented circulant
matrix obtained by vertically concatenating the circulant matrices for each
target point and horizontally concatenating the circulant matrices for each
orbit. For clarity, $\boldsymbol{V}$, $\boldsymbol{x}$, and $\boldsymbol{f}$
are shown in matrix form,

$\displaystyle\boldsymbol{V}={\begin{bmatrix}V_{1,1}&V_{1,2}&\cdots&\cdots&V_{1,\lvert\mathcal{Z}\rvert}\\\
V_{2,1}&V_{2,2}&\cdots&\cdots&V_{2,\lvert\mathcal{Z}\rvert}\\\
\vdots&\vdots&\ddots&\ddots&\vdots\\\ \vdots&\vdots&\ddots&\ddots&\vdots\\\
V_{\lvert\mathcal{J}\rvert,1}&V_{\lvert\mathcal{J}\rvert,2}&\cdots&\cdots&V_{\lvert\mathcal{J}\rvert,\lvert\mathcal{Z}\rvert}\\\
\end{bmatrix}}\boldsymbol{x}={\begin{bmatrix}x_{1}\\\ x_{2}\\\ \vdots\\\
x_{\lvert\mathcal{Z}\rvert}\end{bmatrix}}\boldsymbol{f}={\begin{bmatrix}f_{1}\\\
f_{2}\\\ \vdots\\\ f_{\lvert\mathcal{J}\rvert}\end{bmatrix}}$

Thus, our full integer linear program is written as

$\displaystyle\underset{x}{\text{minimize}}$
$\displaystyle\boldsymbol{1}^{T}\boldsymbol{x}$ subject to
$\displaystyle\boldsymbol{Vx}\geq\boldsymbol{f},$
$\displaystyle\boldsymbol{x}\in\\{0,1\\}^{\lvert\mathcal{Z}\rvert L}$

##  3 Circular Restricted 3 Body Problem

To leverage this formulation in cislunar space, a gravitational model of the
region is needed. The circular restricted three body problem (CR3BP) offers an
accurate model for the dynamics of a small satellite as it moves under the
influence of two much larger bodies (in our case, the Earth and Moon). The
equations are given by,

$\displaystyle\ddot{x}-2\dot{y}$ $\displaystyle=\frac{dU}{dx}$
$\displaystyle\ddot{y}+2\dot{x}$ $\displaystyle=\frac{dU}{dy}$
$\displaystyle\ddot{z}$ $\displaystyle=\frac{dU}{dz}$

where U is an effective potential of the rotating system, given by

$U=\frac{x^{2}+y^{2}}{2}+\frac{1-\mu}{r_{1}}+\frac{\mu}{r_{2}}$

Here, $\mu=\frac{m_{2}}{m_{1}+m_{2}}$ is the mass ratio of the primaries and
$r_{1}$ and $r_{2}$ are the distances of the object to each primary.

$\displaystyle r_{1}$ $\displaystyle=\sqrt{(x+\mu)^{2}+y^{2}+z^{2}}$
$\displaystyle r_{2}$ $\displaystyle=\sqrt{(x-(1-\mu))^{2}+y^{2}+z^{2}}$

This formulation of the CR3BP uses normalized mass, distance, and time units
which are aggregated in Table 1 below. It is important to note that the CR3BP
admits periodic orbits as solutions, which is convenient for the integer
linear program developed. Furthermore, the CR3BP admits 5 stationary points,
called Lagrange points, where the acceleration and velocity of a satellite is
zero. They are commonly referrred to as L1 through L5.

For certain accessibility measures such as the apparent magnitude of the
object, the position of the sun is required in addition to that of the
observer and target. In the rotating frame of the CR3BP, the sun’s motion is
clockwise with period equal to the synodic period of the Earth-Moon system. A
simplified model of the sun’s motion is a coplanar clockwise circular orbit of
radius 1 AU around the Earth-Moon barycenter. Note that the sun is only used
as an illumination source in this model. We do not consider the sun’s
gravitational perturbation on the system. Parameters for the dynamics of this
system are gathered in Table 1 in the appendix.

##  4 Utilizing the Formulation in Cislunar Space

A flow graph of the proposed formulation is shown in Figure 2 below. The first
stage consists of precomputation, where we determine the accessibility for
each (target point, observer location) pair in the entire constellation for
the simulation period. Second, we construct the augmented circulant matrix
$\boldsymbol{V}$ by phasing the accessibility measures. This primarily
consists of a simple circular shift operation. Lastly, we provide the
necessary variables to an optimizer in the form of an integer linear program.
The algorithm is outlined in pseudocode in the appendix. All orbits are
discretized such that successive position states are 0.015 TU $\approx$ 93
minutes apart.

Figure 2: Optimizer Illustration

###  4.1 Target Point Set $\mathcal{J}$

As traffic between the Earth and the Moon is anticipated to increase, an
attractive choice is a transfer orbit from L1 to geosynchronous Earth orbit
(GEO). Figure 3 illustrates this transfer. It consists of a manifold insertion
from L1 followed by a high impulse manuever into a connecting orbit to GEO.
The transfer takes $\approx$ 23.12 days. The discretized position history
(i.e. the blue circles in Figure 3) are target points comprising the set
$\mathcal{J}$. Note that although we define $\mathcal{J}$ as a set of points,
it is actually a sequence. Target point 1 is followed by target point 2 and so
on, tracking the position of a hypothetical target on this path from L1 to
GEO.

Figure 3: L1 to GEO transfer

###  4.2 Observer Orbit Set $\mathcal{Z}$

There are a variety of periodic orbits in the CR3BP. Of these, the most
attractive options offer close approaches to L1 and GEO while moving through
the majority of cislunar space. Furthermore, as per an AFRL objective to
monitor a 30 degree cone emanating from the Earth, orbits in this region are
also attractive. We choose a set of 6 orbits inspired by Gupta et al.[12],
which are plotted in Figure 4 below. The initial conditions are shown in Table
3 in the appendix. All orbits were generated via a single-shooting
differential correction algorithm with initial guesses sampled from data
provided by NASA JPL[15] and Broucke. [16]. All orbits have period $T$ = 6.45
TU $\approx$ 28.008 days. It should be noted that the L1 Lyapunov and L2 Halo
have a fundamental period of $\frac{T}{2}$ = 3.225 TU which means they are
also $T$-periodic. The circular restricted problem does not admit many orbits
with period equal to the synodic period of the Earth-Moon system
(approximately 29.05 days). As a consequence, while the orbits themselves
repeat, the illumination conditions will not. As such, the results of the
formalism apply on a case-by-case basis depending on the sun’s initial phase
(in radians) relative to the Earth-Moon axis, denoted as $\phi_{0}$. We stress
that the candidate orbits are chosen to have periods as close as possible to
the synodic period and are still resonant in the CR3BP sense. This ensures the
difference in illumination conditions is small after $t=6.45$ TU compared to a
randomly chosen constellation period. We aimed to approximately minimize this
difference while preserving candidate orbit diversity. In this work we choose
$\phi_{0}$ = 0 for optimization (Earth-Moon-Sun are collinear at simulation
start). However, we investigate the robustness of the constellation to
different initial phase angles later in the Results section.

Figure 4: Candidate Orbits

###  4.3 Accessibility Metric and Threshold

To quantify the the observability of a target from an observer’s position, we
utilize the target point’s apparent magnitude.[14],[9]Consider an observer at
position $\boldsymbol{r}$, a spherical target with diameter $d$ at position
$\boldsymbol{p}$, and the sun at position $\boldsymbol{s}$. Define the solar
phase angle $\psi$,

$\psi=\frac{(\boldsymbol{p}-\boldsymbol{r})\cdot(\boldsymbol{p}-\boldsymbol{s})}{\|\boldsymbol{p}-\boldsymbol{r}\|\|\boldsymbol{p}-\boldsymbol{s}\|}$

It should be noted that $\psi$ is the familiar angle used in the dot product
formula between two vectors, taking on values between 0 and $\pi$. The diffuse
phase angle function can be written as a function of $\psi$,

$p_{\text{diff}}(\psi)=\frac{2}{3\pi}(\sin{\psi}+(\pi-\psi)\cos{\psi})$

On the interval (0, $\pi$), $p_{\text{diff}}$ is a strictly decreasing
function. The apparent magnitude of a target is given by,

$g(\boldsymbol{p},\boldsymbol{r},\boldsymbol{s})\triangleq
m_{\text{sun}}-2.5*\log_{10}\Big{[}\frac{d^{2}}{\zeta^{2}}*(\frac{a_{\text{spec}}}{4}+a_{\text{diff}}*p_{\text{diff}}(\psi))\Big{]}$

Where $\zeta=\|\boldsymbol{p}-\boldsymbol{r}\|$ is the distance between the
oberver and target, and $a_{\text{spec}},a_{\text{diff}}$ are the specular and
diffuse reflection coefficients of said target, respectively. Values can be
found in Table 2 in the appendix. It should be noted that the solar phase
angle $\psi$ is distinct from the initial phase angle $\phi_{0}$ described in
the Observer Orbits section above. $\phi_{0}$ descibes the initial orientation
of the Sun relative to the Earth-Moon axis. $\psi$ describes the orientation
of an observer, a target point, and the Sun. The value of $g$ increases with
$\psi$, meaning solar phase angles closer to $\pi$ result in an unfavorably
dim target. Every unit increase in apparent magnitude corresponds to a 2.5
times dimmer target. In addition to this, we consider the possibility that the
Earth and Moon may occlude the target. If this is the case, the accessibility
measure is set to $\infty$, indicating an invisible target point. In practice,
the measure will be set to the largest floating point number (often denoted
REALMAX) available on the machine. We choose a threshold value of $M=17$ as
our access measure cutoff. Together, $g$ and $M$ define our accessibility
metric and threshold, respectively.

###  4.4 Coverage Requirement $\boldsymbol{f}$

We seek a set
$\mathcal{F}=\\{f_{1},f_{2},\cdots,f_{\lvert\mathcal{J}\rvert}\\}$ of coverage
requirements, one for each target point in the set $\mathcal{J}$. As our
object of interest travels along the transfer, it travels through each target
point in order. Thus, the coverage requirement for target point $j+1$ should
be that of target point $j$, but shifted (in fact, circularly shifted due to
the periodicity of observer orbits) forward by one index. Mathematically,

$f_{j+1}[i]=f_{j}[(i-1)\bmod L]$ (1)

Where $[\ ]$ symbolizes 0-based indexing into the coverage requirement array
$f_{j}$. This recursive relation is convenient because we are only required to
define the coverage requirement $f_{1}$, with subsequent $f_{j}$ defined by
Eqn. 1. An example set $\mathcal{F}=\\{f_{1},f_{2},f_{3}\\}$ is illustrated
below in Figure 5.

Figure 5: Illustration of Equation 1

##  5 Results

In the above image, there are two departure windows (shown as two impulses) in
each coverage requirement. We choose to investigate how this number of
departure windows affects the satellite constellation size. The number of
requested departure windows reflects the amount of uncertainty in launch time.
With more departure windows, the constellation is optimized for robustness
against target departure delays. However, this will come at the cost of
increasing constellation size. We specifically investigate
$N\in\\{1,2,4,8,16\\}$, where the departure windows are uniformly separated
amongst the $L$ time steps. The powers of two ensure that the indices where
$f_{j}=1$ for $k$ departure windows are a subset of the indices where
$f_{j}=1$ for $2k$ departure windows, all while keeping the windows uniformly
distributed in the interval $[0,L]$. Figure 7 illustrates the result of the
optimizer for $N=16$ departure windows. The constellation consists of 9
satellites distributed on 3:1 and 2:1 resonant orbits. These orbits were
likely chosen by the optimizer because of their close approaches to the moon,
where a significant portion of the spatio-temporal coverage requirement is
located. Figure 14 in the appendix illustrates the results of the optimizer
for every element $N\in\\{1,2,4,8,16\\}$. Figure 7 illustrates the optimal
number of satellites versus the requested number of departure windows. The
constellation size grows with the coverage demand. This is intuitive, since
coverage potential grows with the number of observers in the constellation. We
see linear growth does not hold as the coverage requirement grows. At some
finite constellation size, it is expected that coverage ability saturates.
Significantly stronger compute is necessary to analyze $>$ 16 departure
windows . A continuous coverage requirement, i.e.
$\boldsymbol{f}=\boldsymbol{1}$, would require parallel computation for
reasonable evaluation time.

Figure 6: Optimal Constellation for 16 Departure Windows

Figure 7: Constellation Size vs. Number of Departure Windows

###  5.1 Robustness to Delays

The constellation shows impressive robustness to delays in departure. Figures
9 and 9 show binary indicator maps for constellations optimized for $N=1$ and
$N=16$ departure windows. Figure 15 in the appendix shows the maps for each
$N\in\\{1,2,4,8,16\\}$ requested departure windows. Areas shaded in green
indicate a pair (Time Step, Target Point) which is observable by at least one
satellite in the constellation. Areas shaded in black are not observable. The
red lines (shown as diagonal streaks across Figures 9 and 9) illustrate the
spatio-temporal coverage requirement sent to the optimizer. These regions are
always covered, as required by the formulation. We see that the majority of
the optimizer’s work attempts to cover the target points closer to GEO (target
point No. 300 and onwards). As seen in Figure 3, these points are spaced
further apart (the target is moving quickly in this region), naturally making
it more difficult to observe. With increasing $N$, we make progress in closing
the gap, with most of the binary map displaying green for $N=16$.

Figure 8: Coverage Robustness for 1 Departure Window

Figure 9: Coverage Robustness for 16 Departure Windows

###  5.2 Robustness to Initial Sun Phase Angle, $\phi_{0}$

We acknowledged earlier that because our candidate orbits are not resonant
with the Earth-Moon synodic period, after one simulation period, the
illumination conditions will differ slightly at the start of the second epoch
($t=6.45$ TU). The orbits repeat but now with a slightly different initial sun
phase angle, $\phi_{0}$. As discussed above, we chose an initial sun phase
angle of 0 rad, meaning that at the start of the simulation, the Sun is on the
positive x-axis at distance 1 AU. The constellations are optimized for this
phase angle, but they show impressive robustness against various initial sun
phase angles. We investigated numerous $\phi_{0}$ between 0 and 2$\pi$, given
a constellation already optimized for $N=2$ departure windows (see Figure 15
in the appendix for an illustration of the constellation). Figure 11 shows
that despite changes in $\phi_{0}$, the constellation still satisfies a
majority of the coverage requirement. As such, the constellation’s
effectiveness extends beyond just one simulation period. For example, if
mission requirements mandate a 70$\%$ coverage requirement satisfaction, then
our constellation meets the mandate.

Figure 10: Coverage Requirement Satisfaction vs. $\phi_{0}$ for 2 Departure
Windows

Figure 11: Coverage Robustness for 2 Departure Windows with
$\phi_{0}=162^{\circ}$

Furthermore, Figure 11 shows that there exist several values of $\phi_{0}$
still satisfying $>$ 90$\%$ of the coverage requirement. Figure 11 shows the
coverage map for one such angle, $\phi_{0}=162^{\circ}$. Not only is 97$\%$ of
the coverage requirement satisfied, the impressive ratio of green to black
shows that robustness to delay is still present.

Of all the investigated $\phi_{0}$ in the range [0, 2$\pi$],
$\phi_{0}=250^{\circ}$ was the worst performing constellation. However, it
still satisfied 71.48 $\%$ of the coverage requirement. Figure 13 shows the
coverage map for this choice of $\phi_{0}$, illustrating that despite only
satisfying 71.48 $\%$ of the coverage requirement, the constellation still
boasts impressive robustness to delays (shown via the ratio of green to
black). There is a band of invisibility (appearing as an hourglass figure in
black) in the approximate time step range 70 to 130. This is due to
unfavorable solar phase angles $\psi$ that are incurred in this region. Figure
13 shows the position of the observers as well as the direction of the Sun’s
rays at time step 100. We see that the positions of the Sun, target points,
and observer on the Lyapunov orbit result in solar phase angle
$>\frac{\pi}{2}$. In this scenario the Sun’s rays are too direct, resulting a
poorly illuminated target, making it difficult for optical sensors to monitor.

Figure 12: Coverage Robustness for 2 Departure Windows with
$\phi_{0}=250^{\circ}$

Figure 13: Satellite Positions at Time Step 100

##  6 Conclusions and Future Work

In this work we show the extension of the formulation provided by Lee et al.
[13] to the cislunar regime. We show that the extension is viable in this
regime due to periodic orbits present in the CR3BP. Monitoring specific
transfers (such as L1 to GEO presented in this work) is possible with small
constellation sizes. With more satellites, we get increasing coverage and
thus, robustness to delays. This formulation can be used with any set of
target points, not a trajectory specifically. This method is also fast, as
illustrated by its implementation on personal computers. Significant compute
is not required until the coverage requirement becomes close to continuous.

Rerunning the simulation with higher fidelity models is left to future work.
Instead of the Sun on a coplanar orbit, we may incline it by $\approx
5^{\circ}$ to reflect the Earth-Moon system’s tilt with respect to the
ecliptic. Furthermore, the Sun’s gravitational perturbation can be
incorporated to determine its effect on the stability and periodicity of
candidate orbits in the CR3BP. Lastly, the objective of optimization can
incorporate more complex costs. We can penalize constellations that place
satellites on separate orbits and include a weighting term to reflect the
degree of its importance.

##  7 Acknowledgement

The authors gratefully acknowledge support for this research from the Air
Force Office of Scientific Research (AFOSR), as part of the Space University
Research Initiative (SURI), grant FA9550-22-1-0092 (grant principal
investigator: J. Crassidisfrom University at Buffalo, The State University of
New York).

## References

  * [1] D. C. Folta, T. A. Pavlak, A. F. Haapala, and K. C. Howell, “Preliminary design considerations for access and operations in Earth-Moon L1/L2 orbits,” AAS/AIAA Spaceflight Mechanics Meeting, No. GSFC-500-13-D-0037, 2013.
  * [2] A. Haapala, M. Vaquero, T. Pavlak, K. C. Howell, and D. C. Folta, “Trajectory selection strategy for tours in the earth-moon system,” AAS/AIAA Astrodynamics Specialist Conference, Hilton Head, South Carolina, 2013.
  * [3] P. M. Cunio, M. J. Bever, and B. R. Flewelling, “Payload and Constellation Design for a Solar Exclusion-Avoiding Cislunar SSA Fleet,” Advanced Maui Optical and Space Surveillance Technologies Conference (AMOS), 2020.
  * [4] M. R. Thompson, N. P. Ré, C. Meek, and B. Cheetham, “Cislunar Orbit Determination and Tracking via Simulated Space-Based Measurements,” Advanced Maui Optical and Space Surveillance Technologies Conference, Maui, HI, 2021.
  * [5] A. P. Wilmer, R. A. Bettinger, and B. D. Little, “Cislunar Periodic Orbits for Earth–Moon L1 and L2 Lagrange Point Surveillance,” Journal of Spacecraft and Rockets, Vol. 59, No. 6, 2022, pp. 1809–1820, 10.2514/1.A35337.
  * [6] J. A. Dahlke, A. P. Wilmer, and R. A. Bettinger, “Preliminary comparative assessment of l2 and l3 surveillance using select cislunar periodic orbits,” AAS/AIAA Astrodynamics Specialist Conference, 2022, pp. 1–19.
  * [7] E. M. Zimovan, K. C. Howell, and D. C. Davis, “Near rectilinear halo orbits and their application in cis-lunar space,” 3rd IAA Conference on Dynamics and Control of Space Systems, Moscow, Russia, Vol. 20, 2017, p. 40.
  * [8] M. Klonowski, M. Holzinger, and N. Owens Fahrner, “Optimal Cislunar Architecture Design Using Monte Carlo Tree Search Methods,” 11 2022.
  * [9] J. K. Vendl and M. J. Holzinger, “Cislunar periodic orbit analysis for persistent space object detection capability,” Journal of Spacecraft and Rockets, Vol. 58, No. 4, 2021, pp. 1174–1185.
  * [10] L. Visonneau, Y. Shimane, and K. Ho, “Optimizing Multi-Spacecraft Cislunar Space Domain Awareness Systems via Hidden-Genes Genetic Algorithm,” The Journal of the Astronautical Sciences, Vol. 70, No. 22, 2023.
  * [11] C. Frueh, K. Howell, K. DeMars, S. Bhadauria, and M. Gupta, “Cislunar space traffic management: Surveillance through earth-moon resonance orbits,” 8th European Conference on Space Debris, 2021.
  * [12] M. Gupta, K. C. Howell, and C. Frueh, “CONSTELLATION DESIGN TO SUPPORT CISLUNAR SURVEILLANCE LEVERAGING SIDEREAL RESONANT ORBITS,”
  * [13] H. W. Lee, S. Shimizu, S. Yoshikawa, and K. Ho, “Satellite Constellation Pattern Optimization for Complex Regional Coverage,” Journal of Spacecraft and Rockets, Vol. 57, Nov. 2020, pp. 1309–1327, 10.2514/1.A34657.
  * [14] W. E. Krag, “Visible magnitude of typical satellites in synchronous orbits,” Sept. 1974.
  * [15] M. Vaquero and J. Senent, “Poincare: A Multi-Body, Multi-System Trajectory Design Tool,” 2018, 2014/48975.
  * [16] R. A. Broucke, Periodic orbits in the restricted three-body problem with earth-moon masses.  1968\.
  * [17] C. Rackauckas and Q. Nie, “DifferentialEquations.jl–a performant and feature-rich ecosystem for solving differential equations in Julia,” Journal of Open Research Software, Vol. 5, No. 1, 2017.

##  Appendix A Appendix

###  A.1 CR3BP System Parameters and Candidate Orbit Initial Conditions

System parameters for the Earth-Moon CR3BP with the Sun as illumination source
are presented in Table 1. Initial conditions for candidate orbits are shown
below in Table 3. Orbits were generated via a single shooting differential
corrections algorithm written in Julia. We use the Vern7 ODE propagator
(available in the Julia DifferentialEquations[17] module) in this algorithm.

Table 1: CR3BP Parameters
---
Parameter | Value
Mass parameter, $\mu$ | 1.215058560962404e-02
Length scale, DU (km) | 384400.0
Time scale, TU (s) | 3.751902619517228e+05
Relative distance of Sun (DU) | 389.17794
Angular velocity of Sun (rad/TU) | -0.9253018261815922
Earth Radius, (km) | 6371
Moon Radius, (km) | 1737.4
Table 2: Illumination Parameters
---
Parameter | Value
Sun apparent magnitude, $m_{\text{sun}}$ | -26.74
Target specular reflection coef., $a_{\text{spec}}$ | 0.0
Target diffuse reflection coef., $a_{\text{diff}}$ | 0.2
Target Diameter, $d$ (km) | 0.001
Table 3: Candidate Orbit Initial Conditions
---
| 3:1 Resonant | 2:1 Resonant | 1:1 L1 Lyapunov
$x_{0}$ | 0.13603399956670137 | 0.9519486347314083 | 0.65457084231188
$y_{0}$ | 0 | 0 | 0
$z_{0}$ | 0 | 0 | 0
$\dot{x}_{0}$ | 1.9130717669166003e-12 | 0 | 3.887957091335523e-13
$\dot{y}_{0}$ | 3.202418276067991 | -0.952445273435512 | 0.7413347560791179
$\dot{z}_{0}$ | 0 | 0 | 0
T | 6.45 | 6.45 | 6.45
| 1:1 L2 Lyapunov | L1 Lyapunov | L2 Halo
$x_{0}$ | 0.9982702689023665 | 0.8027692908754149 | 1.1540242813087864
$y_{0}$ | 0 | 0 | 0
$z_{0}$ | 0 | 0 | -0.1384196144071876
$\dot{x}_{0}$ | -2.5322340091977996e-14 | -1.1309830924549648e-14 | 4.06530060663289e-15
$\dot{y}_{0}$ | 1.5325475708886613 | 0.33765564334938736 | -0.21493019200956867
$\dot{z}_{0}$ | 0 | 0 | 8.48098638414804e-15
T | 6.45 | 3.225 | 3.225

###  A.2 Optimizing the Constellation Pattern Vector

Algorithm 1 presents pseudocode for the formulation. In this work, we utilize
Gurobi MATLAB for computation. Computation is done with an M1 Pro CPU using 8
high-performance cores. Algorithm run time is $<$ 300 seconds for each
$N\in\\{1,2,4,8,16\\}$.

Algorithm 1

$\mathcal{Z}$, $\mathcal{J}$, $\boldsymbol{f}$, $g$, $M$, $T$, tstep

tspan $\leftarrow$ 0:tstep:T

L $\leftarrow$ length(tspan)

$\mathcal{S}\leftarrow$ suntrajectory(tspan)

V $\leftarrow$ zeros(L, $\lvert\mathcal{Z}\rvert$, $\lvert\mathcal{J}\rvert$)
$\triangleright$ Precompute

for each $z_{j}\in\mathcal{Z}$ do

$[\mathcal{R},\mathcal{V}]\leftarrow$ cr3bp($z_{j},\text{tspan}$)

for each $(r_{i},s_{i})\in\\{(r_{m},s_{n})\in\mathcal{R}\times\mathcal{S}\
\lvert\ n=m\\}$ do

for each $p_{k}\in\mathcal{J}$ do

m $\leftarrow g(p_{k},r_{i},s_{i})$

V[i, j, k] $\leftarrow$ m $\leq M\ $? 1 : 0

end for

end for

end for

$\mathbb{V}$ $\leftarrow$ [ ] $\triangleright$ Phase Orbits

for $k=1$ to $\lvert\mathcal{J}\rvert$ do

C $\leftarrow$ [ ]

for $j=1$ to $\lvert\mathcal{Z}\rvert$ do

A $\leftarrow$ V[: , j , k]

B $\leftarrow$ zeros(L, L)

for $i=1$ to L do

B[: , $i$] $\leftarrow$ circshift(A, $i-1$)

end for

C $\leftarrow$ horzcat(C, B)

end for

$\mathbb{V}$ $\leftarrow$ vertcat($\mathbb{V}$, C)

end for

c $\leftarrow$ ones($\lvert\mathcal{Z}\rvert$ * L) $\triangleright$ Optimize

result $\leftarrow$ integer-linear-program($\mathbb{V}$, $\boldsymbol{f}$, c)

###  A.3 Extended Results

Figures 14 and 15 show the constellations and indicator maps of all the
optimization runs for each $N\in\\{1,2,4,8,16\\}$.

|
---|---
|
|
Figure 14: Optimized Constellations for Given Number of Departure Windows |
---|---
|
|
Figure 15: Constellation Robustness Binary Maps


# Policy Search with Rare Significant Events:
Choosing the Right Partner to Cooperate with

Paul Ecoffet  Institut des Systèmes Intelligents et de Robotique, Sorbonne
Université4 place JussieuParisFrance75005<EMAIL_ADDRESS>0000-0002-7394-6134 ,  Nicolas Fontbonne  Institut des Systèmes Intelligents
et de Robotique, Sorbonne Université4 place JussieuParisFrance75005
<EMAIL_ADDRESS>,  Jean-Baptiste André  Institut Jean
Nicod, Département d’Études Cognitives, École Normale Supérieure29 rue
d’UlmParisFrance75005<EMAIL_ADDRESS>and  Nicolas Bredeche
Institut des Systèmes Intelligents et de Robotique, Sorbonne Université4 place
JussieuParisFrance75005<EMAIL_ADDRESS>
(2021)

###### Abstract.

This paper focuses on a class of reinforcement learning problems where
significant events are rare and limited to a single positive reward per
episode. A typical example is that of an agent who has to choose a partner to
cooperate with, while a large number of partners are simply not interested in
cooperating, regardless of what the agent has to offer. We address this
problem in a continuous state and action space with two different kinds of
search methods: a gradient policy search method and a direct policy search
method using an evolution strategy. We show that when significant events are
rare, gradient information is also scarce, making it difficult for policy
gradient search methods to find an optimal policy, with or without a deep
neural architecture. On the other hand, we show that direct policy search
methods are invariant to the rarity of significant events, which is yet
another confirmation of the unique role evolutionary algorithms has to play as
a reinforcement learning method.

reinforcement learning, rare significant events, on-policy, on-line,
continuous state and action spaces, cooperation and partner choice, gradient
policy search, direct policy search, evolutionary algorithms, PPO, CMAES

††copyright: none††journalyear: 2021††ccs: Theory of computation Evolutionary
algorithms††ccs: Computing methodologies Reinforcement learning

##  1\. Introduction

We consider a particular class of reinforcement learning problems where only
rare events can result in non-zero rewards and when the agent can experience
at most one positive reward in a limited time. This problem is closely related
to the problem of learning with rare significant events in reinforcement
learning (Bhatnagar et al., 2006; Frank et al., 2008; Ciosek and Whiteson,
2017), where rare events can significantly affect performance (e.g. in network
and communication systems or control problems where failure can be
catastrophic). In this paper, we consider that while significant events occur
independently of the agent’s actions, the agent’s policy determines if a
positive reward should be obtained when such an event occurs. Significant
events are thus defined as unique opportunities to obtain a positive reward
and stop the game. Each opportunity can either be seized for an immediate
reward or ignored if the agent hopes to get a better reward in the future.

We address this problem in the context of an independent, on-line and on-
policy episodic learning task with continuous state and action spaces. The
practical application addressed in this paper is that of an agent learning to
choose a partner for a task that requires cooperation (e.g., predators hunting
a large prey or individuals selecting a lifelong mate). The agent can choose
to cooperate or not with a potential partner, based on the effort this partner
is willing to invest in the cooperation. At the same time, the agent must
invest enough so that its partner also accepts to cooperate. In this setup,
the agent may face partners willing to invest various amount of energy in
cooperation (i.e., a possibly significant event), or even refuse to cooperate
whatever the agent is ready to invest (i.e. a non-significant event).

Results from theoretical biology (McNamara et al., 2008; Campennì and Schino,
2014; Debove et al., 2015; Ecoffet et al., 2020b) have shown that cooperation
with partner choice is optimal only under certain conditions. First, the
number of cooperation opportunities must be large enough that an agent can
refuse to cooperate with a potential partner and still have the opportunity to
meet a more interesting partner. Second, if an agent and its partner both
decide to cooperate, the actual duration of this cooperation must be long
enough to make cooperation with an uninteresting partner significantly costly
(which is the case when there can be only one single successful cooperation
event). Under these conditions, the optimal strategy for an agent is to be
very demanding in choosing its partner.

The question raised in this paper is whether reinforcement learning algorithms
actually succeed in learning an optimal strategy when the necessary conditions
are met. We are particularly interested in how the rarity of significant
events influences convergence speed and performance of policy learning.
Indeed, it is not clear how gradient-based policy search method can deal with
a possibly large number of non-significant events that provide zero-reward.

We use two state-of-the-art methods for on-policy reinforcement learning with
continuous state and action spaces: (1) a deep learning method (PPO (Schulman
et al., 2017)) for gradient policy search and (2) an evolutionary method
(CMAES (Hansen and Ostermeier, 2001)) for direct policy search. While both
methods provide similar results when the agent is always presented with
significant events, policy search methods are not equals when such events
become rarer. While the direct policy method is oblivious to rarity of
significant events, the gradient policy search method suffers significantly
from rarity.

The paper is structured as follows: the reinforcement learning problem with
significant rare events and single reward per episode is formalized, and the
partner choice learning problem is presented as a variation of a continuous
prisoner’s dilemma. Algorithms and results are then presented, and learned
policies are analysed and compared.

##  2\. Methods

###  2.1. Learning with Rare Significant Events

Formally, we consider an independent learner $x_{\bullet}$, called the focal
agent, which is placed in an aspatial environment. At each time step,
$x_{\bullet}$ is presented with either a cooperative partner $x_{i}^{+}\in
X^{+}$ or a non-cooperative partner $x_{j}^{-}\in X^{-}$. $X^{+}$ (resp.
$X^{-}$) is the finite set of all cooperative (resp. non-cooperative) agents,
with both $i$ and $j\in\mathbb{N}$ and $i>0,j\geq 0$. When presented with a
non-cooperative partner $x_{j}^{-}$, the focal agent’s reward will always be
zero. When presented with a cooperative partner $x_{i}^{+}$, the focal agent’s
reward will depend on its own action and that of its partner. (see Section 2.2
for details).

Our objective is to endow the focal agent $x_{\bullet}$ with the ability to
learn how to best cooperate, which implies to negotiate with its potential
partners and decide whether cooperation is worth investing energy in, or not
(see Section 2.3 for details). The focal agent faces an individual learning
problem as it must optimize its own gain over time in a competitive setup,
whether its partners are also learning agents or not. For cooperation to occur
between the focal agent and a partner, the partner must willing to cooperate
(ie. be one of $x_{i}^{+}$) and both the focal agent and the cooperative
partner must estimate that one’s own energy invested in cooperation is worth
the benefits.

We use the standard reinforcement learning framework proposed by Sutton and
Barto (2018) to formalize the learning task from the focal agent’s viewpoint,
which is essentially a single agent reinforcement learning problem.

The focal agent $x_{\bullet}$ interacts with the environment in a discrete
time manner. At each time step $t=0,1,2,...$, $x_{\bullet}$ is in a state
$s\in\mathbb{R}$ which describes its current partner’s investment value, and
plays a continuous value $a\in\mathbb{R}$ which represents its decision to
cooperate ($a>0$) or not ($a<=0$).

Let $\pi_{\theta}$ be the parametrised policy of the focal agent, with
$\theta\in\mathbb{R}^{n}$. The learning task is to search for $\theta^{*}$,
such as:

(1) $\theta^{*}=\underset{\theta}{argmax}J(\theta)$

With $J$ the global function to be optimized, defined as:

(2) $J(\theta)=\mathbb{E}\underset{t}{\sum}r_{t}$

with reward $r_{t}$ at time $t$. Rewards are defined such that
$r\in\mathbb{R}$ and depends on the current state $s$ and action $a$, and are
produced according to the probability generator defined as follow:

(3) $r(s,a)=\left\\{\begin{array}[]{ll}\mathit{payoff}(s,a)&\mbox{with
probability }p\\\ 0&\mbox{otherwise.}\end{array}\right.$

The probability $p\in[0,1]$ determines the probability to encounter a
cooperative agent (i.e. one of $x_{i}^{+}$). The value of $p$ depends on the
setup, and determines how rare significant events occur when $p<1.0$. A
probability of $p=1.0$ means the focal agent $x_{\bullet}$ encounters a
cooperative partner at each time step $t$, with a possible positive reward (if
cooperation is accepted by both agents) that depends on the $\mathit{payoff}$
function. Non-zero rewards become rarer (but still possible) as $p\rightarrow
0$. Note that $\mathit{payoff}(s,a)$ is non-zero only if both the focal agent
and its cooperative partner accept to cooperate. Cf. Section 2.3 for details
on the negotiation process.

The problem presented here is very similar to that of Rare Significant Events
as formulated by Frank et al. (2008). However, our problem differs on two
aspects. Firstly, we consider on-line on-policy search of a parametrised
policy, where the frequency of significant events cannot be controlled.
Secondly, and even more importantly, a learning episode stops right after the
focal agent and one cooperative agent have reached a consensus to cooperate.
If no cooperation is triggered, an episode stops after a maximum number of
iterations $T$, defined as:

(4) $T=\frac{100}{p}\mbox{ time steps}$

It results that the expected number of meetings $M$ is held constant
independently from the value of $p$ (i.e. $\mathbb{E}(M)=100$). It is
therefore possible to obtain episodes of different lengths but with the same
number of significant events.

The situation that is modelled here corresponds to many collective tasks
observed in nature (Bshary and Noe, 2002; Simms and Taylor, 2002; Wilkinson et
al., 2016), where each agent has to balance between looking for partners and
cooperating with the current partner, the latter possibly taking significant
time. As a matter of fact, it has been shown elsewhere (McNamara et al., 2008;
Campennì and Schino, 2014; Debove et al., 2015; Ecoffet et al., 2020b, a) that
optimal partner choice strategies can be reached only when the cost of
cooperation is large (ie. the duration of cooperation is long with regards to
looking for cooperative partners).

###  2.2. Partner Choice and Payoff Function

Whenever the focal agent $x_{\bullet}$ and a cooperative partner $x_{i}^{+}$
interact together, they play a variation of a continuous Prisoner’s Dilemma.
Cooperation actually takes place if both agents deem it worthwhile. The two-
step procedure for partner choice is the following:

  1. (1)

each agent simultaneously announce the investment they are willing to pay to
cooperate;

  2. (2)

each agent then chooses to continue the cooperation based on the investment
announced by its partner and its own.

To simplify notations, we use $x_{\bullet}$ and $x_{i}^{+}$ to represent both
the agents and the investment values they play, i.e. $x_{\bullet}$ (resp.
$x_{i}^{+}$) plays $x_{\bullet}$ (resp. $x_{i}^{+}$). The gain received by the
focal agent $x_{\bullet}$ is defined as:

(5) $P(x_{\bullet},x_{i}^{+})=a\times x_{\bullet}+b\times
x_{i}^{+}-\frac{1}{2}x_{\bullet}^{2}$

With $a,b\geq 0$ and $a+b>0$. This payoff function combines both a prisoner’s
dilemma and a public good game, and was first introduced in Ecoffet et al.
(2020b). Two different equilibria111These are actually Nash equilibria, when
all agents are learning. can be reached for $x_{\bullet}$:

  * •

$x_{d}=a$. This is a sub-optimal equilibrium, which corresponds to an agent
cheating, a typical outcome in the prisoner’s dilemma where an agent maximizes
its own gain, but also minimizes its exposure to defection. This ensure the
best payoff for the agent if it is unable to distinguish a cheater from a
cooperator.

  * •

$x_{c}=a+b$. This is the optimal equilibrium, where both agents cooperate to
maximize their long-term gain.

The public good game is included in the payoff function to help distinguish
between agents that are simply ignoring the cooperation game
($x_{\bullet}=0$), from those who takes part in it, even if they defect
($x_{\bullet}\geq x_{d}$).

The focal agent can get the optimal payoff if it plays $x_{\bullet}=x_{c}$ and
its partner plays $x_{i}^{+}\geq x_{c}$, which can occur if particular
conditions are met when partner choice is enabled. Partner choice can lead to
optimal individual gain whenever a successful cooperation removes the
possibility for further gain with other partners. In other words: the focal
agent can meet with any number of possible partners but will take the gain of
the first and single mutually accepted cooperation offer.

In this paper, we set $a=5$ and $b=5$, therefore $x_{d}=5$ and $x_{c}=10$. The
maximum payoff the agent can obtain is to invest $x_{\bullet}=x_{c}$ with its
partner investing equally $x_{i}^{+}=x_{c}$. In this context,
$P(x_{\bullet},x_{i}^{+})=50$. The focal agent’s investment is bounded as
$0.0\leq x_{\bullet}\leq 15.0$. This is similar for $x_{i}^{+}$.

$P(x_{\bullet},x_{i}^{+})$ and $\mathit{payoff}(s,a)$ (introduced in Equation
3) differs as the $P$ function relates to the game theoretical setting while
the $\mathit{payoff}$ function relates to the reinforcement learning problem.
On the one hand, the $\mathit{payoff}$ function computes the focal
individual’s reward whether or not cooperation was initiated. On the other
hand, $P$ computes the focal individual’s gain that results from a cooperation
game between two agents that accepted to cooperate. However, both functions
are linked. From a notational standpoint, $s$ represents the investment value
of the focal individual $x_{\bullet}$, and $a$ represents the decision to
cooperate and depends on both $s$ and that of its partner $s_{i}^{+}$ (which
is implicit). The return value of $\mathit{payoff}(s,a)$ depends on whether
cooperation was initiated or not. If both agents decided to cooperate, then
the focal agent’s payoff is $\mathit{payoff}(s,a)=P(x_{\bullet},x_{i}^{+})$,
with $P(x_{\bullet},x_{i}^{+})\leq 50$ in this case. If cooperation fails, the
focal agent’s payoff is $\mathit{payoff}(s,a)=0$ (which is obtained without
having to compute $P$). The $payoff$ function in Equation 3 can be written as
follow, with updated notations and assuming $a_{\bullet}>0$ (resp.
$a_{i}^{+}>0$) means the focal agent (resp. partner) is willing to cooperate:

(6)
$\mathit{payoff}(s_{\bullet},a_{\bullet})=\left\\{\begin{array}[]{ll}P(x_{\bullet},x_{i}^{+})&\mbox{if
$a_{\bullet}>0$ and $a_{i}^{+}>0$}\\\ 0&\mbox{otherwise.}\end{array}\right.$

###  2.3. Behavioural Strategies

For each interaction, the focal agent’s investment value
$x_{\bullet}\in[0,15]$ is computed, and when the investment value of its
partner is known, its decision to cooperate $a_{\bullet}\in\mathbb{R}$ is
computed to determine if cooperation should be pursued or not. Each value is
provided by a dedicated decision module:

  * •

the investment module which provides the cost $x_{\bullet}$ that the focal
agent is willing to invest to cooperate. This module takes no input as it is
endogenous to the agent (i.e. the proposed cost $x_{\bullet}$ is fixed
throughout an episode);

  * •

the choice module takes both the focal agent’s own investment value
($x_{\bullet}$) and that of its partner ($x_{i}^{+}$ or $x_{j}^{-}$), and
computes $a_{\bullet}$, which is used to determine if cooperation is an
interesting choice ($a_{\bullet}>0$) or not ($a_{\bullet}\leq 0$). The choice
module is essentially a function
$f_{choice}(x_{\bullet},x_{partner})\rightarrow a_{\bullet}$ with
$x_{partner}\in X^{+}\cup X^{-}$. The parameters of the function are learned,
and the decision to cooperate is computed (as the decision to cooperate is
conditioned by the partner’s investment).

With respect to the focal individual, Section 3 describes how the investment
and choice modules are defined and how learning is performed depending on the
learning algorithm used.

Cooperative partners $x_{i}^{+}$ and non-cooperative partners $x_{j}^{-}$ also
use similar decision modules, providing investment and choice values. However,
all use deterministic fixed strategies, which may differ from one partner to
another. Firstly, non-cooperative partners $x_{j}^{-}$ all follow the same
strategy. Both the investment value $x_{j}^{-}$ and the decision to cooperate
$a_{j}^{-}$ are always $0$, $\forall j$.

Secondly, cooperative partners $x_{i}^{+}$ each follows a stereotypical
cooperative strategy depending on the value $i$. Each cooperating partner
invests a fixed value $x_{i}^{+}\in[0,15]$ defined as:

(7) $x_{i}^{+}=\frac{i-1}{i_{max}}\times 15,\,i\in\\{1,\ldots,i_{max}\\}$

Each cooperative partner then accepts to cooperate if the focal agent’s
investment value $x_{\bullet}$ is greater or equal to their investment, which
is written as follow:

(8) $a_{i}^{+}=\left\\{\begin{array}[]{ll}1&\mbox{if $x_{\bullet}\geq
x_{i}^{+}$}\\\ -1&\mbox{otherwise.}\end{array}\right.$

In the following, there are $i_{max}=31$ cooperating partners ($x_{i}^{+}\in
X^{+},i\in\\{1,\ldots,31\\}$). Following Eq.8, this means cooperating partner
$x_{1}^{+}$ (resp. $x_{2}^{+}$, …, $x_{31}^{+}$) plays $0$ (resp. $0.5$, …,
$15$).

##  3\. Parameter Settings and Algorithms

We use two reinforcement learning algorithms: a gradient policy search
algorithm (PPO) and a direct policy search algorithm (CMAES). Both algorithms
are used to learn the parameters of the focal agent’s decision modules.

For both algorithms, the performance of a policy (i.e. the return or the
fitness, depending on the vocabulary used) during one episode is computed as
the sum of rewards during the episode (cf. Section 2.1), which is either zero,
or the value of the unique non-zero reward obtained before the episode ends.

###  3.1. Proximal Policy Optimization

The deep reinforcement learning Proximal Policy Optimisation (PPO) (Schulman
et al., 2017) is a variation of the Policy Gradient algorithm (Sutton and
Barto, 2018). Policy gradient algorithms maximize the global performance by
updating the parameters $\theta$ of the policy $\pi$ (cf. Eq. 2).

Though, as the expected value of a certain state-action pair varies according
to the policy itself, updating a new policy from samples acquired from an old
policy may cause inaccurate predictions, as the expected value of an action-
state pair may be wrong with respect to the new policy. PPO ensures that the
policy generated from the samples of the new policy does remain in a so-called
trust region at each learning step.

As we are dealing with episodes and do not want to encourage the focal agent
to act in the least amount of time steps as possible, the discount factor is
set to $\gamma=1.0$, as recommended by Sutton and Barto (2018, p.68). The PPO
hyper-parameters used are reported in Table 1.

Parameters | Values
---|---
Learning rate | $0.005$
Optimiser Algorithm | SGD
Number of optimisation epochs | $10$
Minibatch size | $128$
Batch size | $4000$
Discount factor $\gamma$ | 1.0
Search space PPO-MLP ($\theta_{MLP}$) | $\mathbb{R}^{33}$
Search space PPO-DEEP ($\theta_{DEEP}$) | $\mathbb{R}^{133894}$
Table 1. Parameters for the PPO algorithm

The investment and choice modules are both represented as Artificial Neural
Networks (ANN). A module is composed of both a decision network and a Value
function, as PPO runs as an actor-critic algorithm. The Value function network
has the same layout as the decision network, but only output the (continuous)
value of the state.

The decision network for the investment module is a simple neural network with
one single input set to $1.0$, no hidden layer and two outputs: the investment
mean $m$ and standard deviation $\sigma$. The investment $x_{\bullet}$ is
picked along the distribution $\mathcal{N}(m,\sigma^{2})$ and clipped between
0 and 15. The continuous stochastic action selection is essential to the PPO
search algorithm.

The decision network for the choice module is a multilayer perceptron with two
input neurons and two output neurons (for accepting or refusing cooperation).
The output neurons use a linear activation function, and a softmax
probabilistic choice is done to choose which action to make (accept or
decline). Hidden units use an hyperbolic tangent activation function. A bias
node is used, that projects on both the hidden layer(s) and output neurons.
The Value Function estimator use the same architecture as the choice neural
networks, with only one output.

In Section 4, two different architectures are evaluated, which we refer to as
PPO-MLP and PPO-DEEP. While both use the decision network for the investment
module described before, they differ with respect to the architecture used for
the choice module. PPO-MLP implements a single hidden layer with 3 neurons,
and PPO-DEEP implements a deep architecture with two hidden layers, each with
$256$ neurons. While PPO-DEEP may seem overpowered at first sight, over-
parametrization has been shown to be very effective in deep learning as
multiple gradients can be followed in wide neural networks (Du and Lee, 2018;
Neyshabur et al., 2018; Allen-Zhu et al., 2019).

All parameter values and module architecture result from an extensive search
(summarised in the Supplementary Materials). In particular, a grid search was
performed to select the best values for each parameters, including the
learning rate ($lr$). The number of Simple Gradient Descent iterations, the
batch size and the mini-batch size had little impact on neither performance
nor convergence. In addition, we performed additional experiments to evaluate
the impact of using (1) a discount factor $\gamma<1.0$ (i.e. $0.9$, $0.99$ and
$0.999$) and (2) PPO without actor-critic. None of these settings provided
better (or even comparable) results to those obtained with the parameters used
in Table 1.

###  3.2. Covariance Matrix Adaptation Evolution Strategy

The Covariance Matrix Adaptation Evolution Strategy (CMAES) is an optimisation
algorithm that does black box optimisation and is derivative-free (Hansen and
Ostermeier, 2001). The goal of CMAES is to find $\theta^{*}$ that maximizes
(or minimizes) a continuous function $f$. CMAES does not require the function
to be convex or differentiable, and relies on stochastic sampling around the
current estimate of the solution. CMAES creates a population of size $\lambda$
using a multivariate Gaussian distribution. Each individual of the population
is evaluated and CMAES then updates its distribution estimation based on the
average of the sampled agents weighted by their evaluation rank. Furthermore,
the covariance matrix of the multivariate Gaussian distribution is updated so
that the distribution is biased toward the most promising direction.

The investment module is represented as a single real value (the investment),
which is clipped between 0 and 15 when used. The partner choice module is a
neural network with 2 inputs, one hidden layer with three neurons and two
neurons on the output layer used to compute the probability to accept or
refuse cooperation. A softmax probabilistic choice is made to choose which
action to make. A bias node is also used, neurons from the hidden layer use an
hyperbolic tangent activation function, and the output units use a linear
activation function. There are $17$ neural network weights.

The parameters for both modules are compiled into a single vector of real
values. To make the search space similar to that of PPO, dummy parameters are
added to the vector (i.e. values which can be modified by the algorithm, but
with no impact on the outcome) to reach a total number of $34$ real values
(i.e. $\Theta\in\mathbb{R}^{34}$).

Table 2 summarizes the parameters used for the CMAES algorithm. As CMAES is
mostly parameter-free, there were no need to perform extensive preliminary
search, and we used the default values. We choose $\sigma_{init}=1.0$ for the
initial standard deviation and a vector of zeros as initial guess. The
population size $\lambda$ is the default population size in the python CMAES
implementation (Hansen et al., 2020), i.e. $\lambda=4+\lfloor
3\times\ln(N)\rfloor=14$ with $N$ the number of dimensions in the model. Once
the $\lambda$ candidate solutions are evaluated, a new population is generated
according to their performance. A new population is generated every 14
episodes, and so forth until the evaluation budget is consumed.

A candidate solution for the focal agent is evaluated on one episode only,
which length may vary depending on when the focal agent and its partner both
accepts to cooperate (maximal duration defined in Eq. 4).

Parameter | Value
---|---
Population size | 14
Number of episode per evaluation | 1
$\sigma_{init}$ | 1.0
Search space ($\theta_{CMAES}$) | $\mathbb{R}^{34}$
Table 2. Parameters for the CMAES algorithm

##  4\. Results

The environment, the models and the learning algorithms are implemented with
ray222https://docs.ray.io/en/master/,
rllib333https://docs.ray.io/en/master/rllib.html and
pytorch444https://pytorch.org/. We use the cma555https://pypi.org/project/cma/
package in python for the CMAES implementation. Source code is available at
https://github.com/PaulEcoffet/RLCoopExp/releases/tag/v1.1.

For a given value of probability of rare significant events $p$, we performed
$24$ independent runs for each algorithm. A run lasts $200\,000$ episodes. The
maximum duration of an episode is fixed as described in Section 2.1 so the
expected number of significant events remains identical independently from the
actual rarity throughout one episode (cf. equation 4). In practical, an
episode lasts at most $100$ (resp. $200$, $500$, $1000$) iterations for
$p=1.0$ (resp. $0.5$, $0.2$, $0.1$).

Performance of the current policy is plotted every $4000$ iterations, which
corresponds to the batch size used by both PPO instances for learning. As
episodes last significantly shorter than $4000$ iterations this means the
policy’s performance is averaged. For CMAES, we extract the best policy of the
current generation and re-evaluate it $10$ times (i.e. for $10$ episodes) to
get a similarly averaged performance. Results are shown on figures with a data
point every $1000$ episodes.

###  4.1. Learning when All Events are Significant

Figure 1 shows the performance throughout learning for CMAES, PPO-DEEP and
PPO-MLP when $p=1.0$ (i.e. the focal agent faces only cooperative partners).
Each Figure shows $24$ curves corresponding the $24$ independent runs. Both
PPO versions and CMAES are shown to learn near optimal policies
($performance\rightarrow 50$) in almost all runs. CMAES is the fastest to
converge, and PPO-DEEP (despite the huge number of dimensions) is faster than
PPO-MLP. On the other hand, CMAES offers less robustness as $20$ (out of $24$)
runs with CMAES reach a performance above $40$, to be compared to $23$ (out of
$24$) runs with PPO-MLP and $24$ runs with PPO-DEEP.

In order to better compare the quality of the policies learned by each
algorithm, the best policy from the end of each run is selected and re-
evaluated for $1000$ extra episodes without learning. Results are shown in
Figure 2 with all three methods faring similar performance. The median value
for CMAES ($47.64$) is only slightly more than that of PPO-DEEP ($46.99$) and
PPO-MLP ($45.58$).

Therefore, we conclude that all three algorithms provide excellent and
comparable results when only significant events are experienced ($p=1.0$).

Figure 1. Performance of the best policy throughout learning with CMAES (top),
PPO-DEEP (center) and PPO-MLP (bottom), with 24 independent runs per method,
for $200*10^{3}$ episodes. There are $20/24$ runs that produced a policy where
performance above $40$ with CMAES, $20/24$ for PPO-DEEP and $23/24$ for PPO-
MLP. Note that PPO-DEEP produces $24/24$ runs with performance above $40$
around episode $80*10^{3}$, with performance occasionally degrading and
immediately recovering for some runs afterwards due to the learning step size
(see Annex for further analysis). Figure 2. Performance of the best policies
from CMAES, PPO-DEEP and PPO-MLP with $p=1.0$ after re-evaluating policies for
$1000$ episodes without learning. Two-tailed Mann-Whitney U-test, $n=24$,
gives $p-value=0.12$ (CMAES vs. PPO-DEEP), $p-value=0.019$ (CMAES vs. PPO-
MLP), $p-value=0.018$ (PPO-DEEP vs. PPO-MLP). Median values and Median
Absolute differences are: CMAES (median=$47.64$, MAD=$4.72$) is only slightly
more than that of PPO-DEEP (median=$46.99$, MAD=$13.04$) and PPO-MLP
(median=$45.58$, MAD=$1.86$).

###  4.2. Learning when Significant Events are Rare

Figure 3 show the performance of the agent throughout its learning with both
PPO algorithms and the CMAES algorithm for different conditions of rare
significant events ($p\in\\{0.1,0.2,0.5\\}$), as well as with the control
condition when all events are significant ($p=1.0$, taken from the previous
Section). Each figure shows the mean performance of $24$ independent runs per
conditions, compiling each setup by tracing the median performance and 95%
confidence interval from the $24$ runs.

CMAES is only marginally impacted when significant events become rarer (i.e.
$p<1.0$), with all setups showing convergence towards a similar performance
value close to the optimal (above $40$). While PPO-DEEP fares better than PPO-
MLP for $p<1.0$, both are largely affected. In the extreme case where $p=0.1$,
the average performance of $35.7\pm 5.2$ for PPO-DEEP and $24.9\pm 4.2$ of
PPO-MLP, to be compared to $46.2\pm 3.2$ for CMAES.

Figure 4 shows the results for the additional analysis where the best policy
from each run for each condition $p\in\\{0.1,0.2,0.5,1.0\\}$ is selected and
re-evaluated for $1000$ extra episodes without learning and with the condition
$p=1.0$ (i.e. only significant events matter here). Results confirm that the
difference in the performance of policies obtained with CMAES compared to
either versions of PPO widens as significant events become rarer ($p<1.0$)
with both PPO-MLP and PPO-DEEP faring significantly worse than CMAES (p-value
$<0.0001$, Mann-Whitney U-test).

Figure 3. Performance of the best policies (median and $95\%$ confidence
interval) throughout learning with CMAES, PPO-DEEP and PPO-MLP for the 3
conditions with rare significant events $(p\in\\{0.1,0.2,0.5\\})$ and 1
control condition ($p=1.0$, same data as shown in Fig.1), for the first
$75*10^{3}$ episodes (out of $200*10^{3}$). Figure 4. Performance of the best
policies (medians and quartiles) from CMAES, PPO-DEEP and PPO-MLP with
$p\in\\{0.1,0.2,0.5,1.0\\}$ after re-evaluating policies for $1000$ episodes
without learning. Two-tailed Mann-Whitney U-test, $n=24$ marked as: * for
$p-value<0.05$, ** for $p-value<0.01$, *** for $p-value<0.001$ and **** for
$p-value<0.0001$.

###  4.3. Analysing the Best Policies for Partner Choice

In order to better understand why policies’ performance differ among learning
algorithms and conditions, the agent’s policy obtained at the end of each run
is extracted and analysed (i.e. $24$ policies per algorithm per condition).

Figure 5 illustrates the outcome of the Investment Module ($x_{\bullet}$),
i.e. the investment value offered by the focal agent when faced with a
potential partner. It is obtained by measuring the investment value of the
focal agent666Note that for CMAES the Investment Module follows a
deterministic policy (but not the Choice Module). Therefore, it would have
been equivalent to take the investment value from the policy parameters in
that particular case. from $1000$ episodes with $p=1.0$ and without learning.
Policies learned with CMAES play close to $x_{c}=10$, which is the optimal
play for the payoff function (Section 2.2), whatever the frequency of
significant events. As expected, this is different for policies learned with
PPO, as the outcome values of the Investment Module are significantly lower
when the frequency of significant events decreases ($p<1.0$).

Figure 6 illustrates the investment values played by cooperative partners,
when the focal agent accepts to cooperate (whether or not cooperation will
actually take place, as it also depends on the partner’s acceptance). In other
words, it represents how demanding is the focal agent with respects to its
partners’ intention to invest in cooperation. The probability to accept
cooperation is computed for the policies of each run. Each policy is presented
with all $31$ possible cooperative partners, $100$ times each, to estimate the
focal agent strategy. While CMAES produced consistent policies that follow
quasi-identical strategies for all conditions (ie. accepting partners that
invest close to the optimal $x_{c}=10$ or above), this is not the case for PPO
policies which are less demanding for lower value of $p$, with many of the
policies learned by PPO-MLP with condition $p=0.1$ actually accepting any
partners). PPO-DEEP policies fare better than PPO-MLP policies, but still
worse than policies learned with CMAES when significant events are rarer.

Figure 7 takes a detailed look at the results shown in Figure 6. It shows the
strategy profile for partner choice by the best policy obtained with each
algorithm in each condition. Focal agents obtained with CMAES follow an
efficient and clear-cut strategy: they play the optimal investment value
($x_{\bullet}=x_{c}=10$, green vertical line) and accept partners only when
those play a similar or better value ($x_{i}^{+}\geq 10$, blue line). Policies
obtained with PPO-DEEP and PPO-MLP either follow roughly the same profile with
a more stochastic behaviour (PPO-MLP policies for $p=1.0$ and $0.5$, PPO-DEEP
policies for $p=1.0$, $0.2$ and $0.1$) or display a selective strategy,
choosing partners only when they play close to the optimal investment value
$x_{i}^{+}\approx x_{c}$. Only PPO-MLP produced policies which are clearly
sub-optimal for $p=0.2$ and $p=0.1$, with a mean investment below the optimal
investment value $x_{\bullet}<x_{c}$.

Figure 5. Investment value of the focal agent given by the Investment Module
for the best learned policies with CMAES (blue), PPO-DEEP (orange) and PPO-MLP
(green) algorithms, for each condition $p$. Each violin graph represents the
results of the outcome of the $24$ best policies for a given algorithm and
condition after being re-evaluate for $1000$ episodes without learning. Figure
6. Decision to accept to cooperate taken by the focal agent, when facing a
cooperative partner with a particular investment value. Results for CMAES
(blue), PPO-DEEP (orange) and PPO-MLP (green) are shown as violin graph.
X-axis: algorithms and conditions, Y-axis: partner’s investment value for
which the focal agent accept to cooperate.  Figure 7.  Analysis of the Partner
Choice module for all conditions (by columns: $p\in\\{0.1,0.2,0.5,1.0\\}$) and
all algorithms (top: CMAES, center: PPO-DEEP, bottom: PPO-MLP). For each
setup, only the best policy is shown. Each graph plots the probability to
accept cooperation for the focal agent following the best policy (y-axis)
depending on its partner’s proposed investment (x-axis). Data are computed by
presenting each of the $31$ possible cooperative partners to the focal agent
for $100$ iterations as policies are stochastic. The green vertical line
represent the mean investment of the focal agent.

##  5\. Concluding Remarks

In this article, we focused on an on-policy reinforcement learning problem of
an autonomous agent that needs to maximize its gain when interacting with
other agents, with whom our agent may or may not decide to cooperate. The
peculiarity of this problem is to present a (very) small number of significant
events during which the agent can obtain only one single positive reward. The
challenge is therefore to learn how to best choose a partner, by making a
compromise between the chances of finding a better partner, and the cost of an
interaction.

We studied the dynamics of two reinforcement learning methods: a gradient
policy search algorithm and a direct policy search algorithm with an evolution
strategy. Both algorithms succeeded in learning policies that make an optimal
use of partner choice when interaction opportunities are frequent. However,
the two algorithms differ fundamentally when interaction opportunities are
rare. The direct policy search algorithm shows total robustness, while the
gradient policy search algorithm collapses, resulting in sub-optimal policies.

The robustness of the direct policy search method can be expected as the
sequential and temporal aspects of the task is lost within one evaluation. As
long as the evaluation time is long enough to sample the whole population of
relevant partners, there is no cost nor change in the algorithm dynamics to
deal with a situation where significant events are lost within a longer
sequence, but still of the same number. Such independence to action frequency
and delayed rewards have actually been observed elsewhere, though for
different problems (e.g.: robotic control problem (Salimans et al., 2017)).
This is of course different for the gradient policy search method, where
increased rarity means that many learning steps will be performed with zero-
reward, resulting in poor gradient information most of the time. Not only this
slows down learning, even with a similar number of iterations, but it also
prevents learning from converging towards a truly optimal partner choice
strategy. This remains true even when a large search space is considered, in
which over-parametrization in deep neural networks help gradient search (Du
and Lee, 2018; Neyshabur et al., 2018).

The broader motivation behind this work is to identify reinforcement learning
problems for which evolutionary algorithms as a direct policy search method
offer a competitive advantage over gradient policy research methods (cf. also
(Igel, 2003; Morse and Stanley, 2016; Such et al., 2017; Conti et al., 2017;
Salimans et al., 2017; Sigaud and Stulp, 2019; Pagliuca et al., 2020)). The
take-home message that emerges from this paper is that one of these problems
occurs when important events are rare, for which direct policy search shows an
invariance to rarity.

As a final remark, it may be tempting to relate the problem of rare
significant events with that of sparse rewards, which has gain a lot of
attention recently (Konidaris and Barto, 2006; Jaderberg et al., 2017;
Riedmiller et al., 2018). However, they differ fundamentally as significant
events may be rare, but eventually occur. This is not the case with sparse
rewards, which occurrences are conditioned by the policy itself (e.g. a
robotic arm must be within the length of a target to trigger a reward) and may
never be obtained. We also argue that problems where significant events are
rare rather than sparse may be more numerous than expected: a complex
environment offers multiple learning opportunities, as long as one is able to
seize them as they arise.

###### Acknowledgements.

This work is funded by the Agence Nationale pour la Recherche under Grant No
ANR-18-CE33-0006. We would like to thank Yann Chevaleyre, Olivier Sigaud and
Mathieu Seurin for feedbacks and comments.

## References

  * (1)
  * Allen-Zhu et al. (2019) Zeyuan Allen-Zhu, Yuanzhi Li, and Zhao Song. 2019\.  A Convergence Theory for Deep Learning via Over-Parameterization. In Proceedings of the 36th International Conference on Machine Learning, Kamalika Chaudhuri and Ruslan Salakhutdinov (Eds.), Vol. 97. PMLR, 242–252.
  * Bhatnagar et al. (2006) Shalabh Bhatnagar, Vivek S Borkar, and Madhukar Akarapu. 2006\.  A simulation-based algorithm for ergodic control of Markov chains conditioned on rare events.  Journal of Machine Learning Research 7, Oct (2006), 1937–1962.
  * Bshary and Noe (2002) Redouan Bshary and Ronald Noe. 2002.  Biological markets - The ubiquitous influence of partner choice on the dynamics of cleaner fish - Client reef fish interactions. In 90th Dahlem Workshop on Genetic and Cultural Evolution of Cooperation, Peter Hammerstein (Ed.). MIT Press, Berlin, Germany, 167\.
  * Campennì and Schino (2014) Marco Campennì and Gabriele Schino. 2014.  Partner choice promotes cooperation: The two faces of testing with agent-based models.  Journal of Theoretical Biology 344 (2014), 49–55.
  * Ciosek and Whiteson (2017) Kamil Andrzej Ciosek and Shimon Whiteson. 2017.  OFFER: Off-Environment Reinforcement Learning. In Proceedings of the Thirty-First AAAI Conference on Artificial Intelligence, February 4-9, 2017, San Francisco, California, USA, Satinder P. Singh and Shaul Markovitch (Eds.). AAAI Press, 1819–1825.
  * Conti et al. (2017) Edoardo Conti, Vashisht Madhavan, Felipe Petroski Such, Joel Lehman, Kenneth O. Stanley, and Jeff Clune. 2017.  Improving Exploration in Evolution Strategies for Deep Reinforcement Learning via a Population of Novelty-Seeking Agents.  CoRR abs/1712.06560 (2017).  arXiv:1712.06560
  * Debove et al. (2015) Stéphane Debove, Jean Baptiste André, and Nicolas Baumard. 2015.  Partner choice creates fairness in humans.  Proceedings of the Royal Society B: Biological Sciences 282, 1808 (2015), 1–7.
  * Du and Lee (2018) Simon Du and Jason Lee. 2018.  On the Power of Over-parametrization in Neural Networks with Quadratic Activation. In Proceedings of the 35th International Conference on Machine Learning, Jennifer Dy and Andreas Krause (Eds.), Vol. 80. PMLR, 1329–1338.
  * Ecoffet et al. (2020a) Paul Ecoffet, Jean-Baptiste André, and Nicolas Bredeche. 2020a.  Learning to Cooperate in a Socially Optimal Way in Swarm Robotics. In Artificial Life Conference Proceedings. MIT Press, 251–259.
  * Ecoffet et al. (2020b) Paul Ecoffet, Nicolas Bredeche, and Jean-Baptiste André. 2020b.  Nothing better to do? Environment quality and the evolution of cooperation by partner choice.  bioRxiv (2020).
  * Frank et al. (2008) Jordan Frank, Shie Mannor, and Doina Precup. 2008\.  Reinforcement Learning in the Presence of Rare Events. In Proceedings of the 25th International Conference on Machine Learning. ACM, New York, NY, USA, 336–343.
  * Hansen et al. (2020) Nikolaus Hansen, Youhei Akimoto, and Petr Baudis. 2020\.  CMA-ES/pycma: r3.0.3.  (April 2020).  https://doi.org/10.5281/zenodo.3764210
  * Hansen and Ostermeier (2001) Nikolaus Hansen and Andreas Ostermeier. 2001.  Completely derandomized self-adaptation in evolution strategies.  Evolutionary computation 9, 2 (2001), 159–195.
  * Igel (2003) Christian Igel. 2003\.  Neuroevolution for reinforcement learning using evolution strategies. In The 2003 Congress on Evolutionary Computation, 2003. CEC’03., Vol. 4. IEEE, 2588–2595.
  * Jaderberg et al. (2017) Max Jaderberg, Volodymyr Mnih, Wojciech Marian Czarnecki, Tom Schaul, Joel Z. Leibo, David Silver, and Koray Kavukcuoglu. 2017.  Reinforcement Learning with Unsupervised Auxiliary Tasks. In Proceedings of the 5th International Conference on Learning Representations, ICLR.
  * Konidaris and Barto (2006) George Konidaris and Andrew Barto. 2006.  Autonomous shaping: Knowledge transfer in reinforcement learning. In Proceedings of the 23rd international conference on Machine learning. 489–496.
  * McNamara et al. (2008) John M. McNamara, Zoltan Barta, Lutz Fromhage, and Alasdair I. Houston. 2008.  The coevolution of choosiness and cooperation.  Nature 451, 7175 (2008), 189–192.
  * Morse and Stanley (2016) Gregory Morse and Kenneth O. Stanley. 2016.  Simple Evolutionary Optimization Can Rival Stochastic Gradient Descent in Neural Networks. In Proceedings of the Genetic and Evolutionary Computation Conference 2016 (GECCO ’16). Association for Computing Machinery, New York, NY, USA, 477–484.
  * Neyshabur et al. (2018) Behnam Neyshabur, Zhiyuan Li, Srinadh Bhojanapalli, Yann LeCun, and Nathan Srebro. 2018.  Towards Understanding the Role of Over-Parametrization in Generalization of Neural Networks.  CoRR abs/1805.12076 (2018).  arXiv:1805.12076
  * Pagliuca et al. (2020) Paolo Pagliuca, Nicola Milano, and Stefano Nolfi. 2020\.  Efficacy of Modern Neuro-Evolutionary Strategies for Continuous Control Optimization.  Frontiers in Robotics and AI 7 (2020), 98.  https://doi.org/10.3389/frobt.2020.00098
  * Riedmiller et al. (2018) Martin Riedmiller, Roland Hafner, Thomas Lampe, Michael Neunert, Jonas Degrave, Tom van de Wiele, Vlad Mnih, Nicolas Heess, and Jost Tobias Springenberg. 2018.  Learning by Playing Solving Sparse Reward Tasks from Scratch. In Proceedings of the 35th International Conference on Machine Learning, Jennifer Dy and Andreas Krause (Eds.), Vol. 80. PMLR, 4344–4353.
  * Salimans et al. (2017) Tim Salimans, Jonathan Ho, Xi Chen, and Ilya Sutskever. 2017\.  Evolution Strategies as a Scalable Alternative to Reinforcement Learning.  CoRR abs/1703.03864 (2017).  arXiv:1703.03864
  * Schulman et al. (2017) John Schulman, Filip Wolski, Prafulla Dhariwal, Alec Radford, and Oleg Klimov. 2017.  Proximal Policy Optimization Algorithms.  CoRR abs/1707.06347 (2017).  arXiv:1707.06347
  * Sigaud and Stulp (2019) Olivier Sigaud and Freek Stulp. 2019.  Policy search in continuous action domains: an overview.  Neural Networks 113 (2019), 28–40.
  * Simms and Taylor (2002) Ellen L. Simms and D. Lee Taylor. 2002.  Partner Choice in Nitrogen-Fixation Mutualisms of Legumes and Rhizobia.  Integrative and Comparative Biology 42, 2 (04 2002), 369–380.
  * Such et al. (2017) Felipe Petroski Such, Vashisht Madhavan, Edoardo Conti, Joel Lehman, Kenneth O. Stanley, and Jeff Clune. 2017.  Deep Neuroevolution: Genetic Algorithms Are a Competitive Alternative for Training Deep Neural Networks for Reinforcement Learning.  CoRR abs/1712.06567 (2017).  arXiv:1712.06567
  * Sutton and Barto (2018) Richard S. Sutton and Andrew G. Barto. 2018.  Reinforcement learning: An introduction.  MIT press, Cambridge, MA.
  * Wilkinson et al. (2016) Gerald S. Wilkinson, Gerald G. Carter, Kirsten M. Bohn, and Danielle M. Adams. 2016.  Non-kin cooperation in bats.  Philosophical Transactions of the Royal Society B: Biological Sciences 371, 1687 (2016).