# Exploring the Focusing Mechanism of the NuMI Horn Magnets

Katsuya Yonehara<EMAIL_ADDRESS>Fermi National Accelerator Laboratory,
Batavia, IL 60510, USA  Sudeshna Ganguly  Fermi National Accelerator
Laboratory, Batavia, IL 60510, USA  Don Athula Wickremasinghe  Fermi National
Accelerator Laboratory, Batavia, IL 60510, USA  Pavel Snopok  Illinois
Institute of Technology, Chicago, IL 60616, USA  Yiding Yu  Illinois Institute
of Technology, Chicago, IL 60616, USA

###### Abstract

Neutrinos at the Main Injector (NuMI) is a project at Fermilab that provides
an intense beam of neutrinos used by a number of experiments. NuMI creates a
beam of pions that decay into neutrinos, muons, and other particles. Muons are
registered by the muon monitors. Magnetic horns are the key elements of the
NuMI beam line. This paper uses the muon beam profile observed at the muon
monitors to study the NuMI horn focusing mechanism. It is found that the horn
magnet generates dipole and quadrupole fields to focus pions. This suggests
that the optics of the horn magnet are predominantly linear. Our study shows
that the muon beam profile accurately detects the horn current within
$\pm$0.05%.

††preprint: APS/123-QED

##  I Introduction

Fermilab hosts the world’s most powerful accelerator-produced neutrino beam
facility to characterize neutrino properties. The Neutrinos at the Main
Injector (NuMI) facility operates an intense proton beam striking the high-
power target to produce positively and negatively charged pions via hadronic
interactions. The primary pion decay process for neutrino production is as
follows:

$\displaystyle\pi^{+}\rightarrow\nu_{\mu}+\mu^{+},$
$\displaystyle\pi^{-}\rightarrow\bar{\nu}_{\mu}+\mu^{-}.$

A magnetic horn plays a crucial role in the neutrino beamline to select pions
of the correct charge, collect them in a wide energy range of 1–20 GeV, and
transport them to a decay volume to provide the correct desired flavor and
chirality111We omit the discussion of right-handed neutrinos or left-handed
anti-neutrinos in this paper. neutrinos to the neutrino detectors. Given its
pivotal role, the horn needs to produce a high-reliability and high-stability
magnetic field with low systematic uncertainty. Currently, the NuMI neutrino
beam is provided to the NOvA [1, 2], MicroBooNE [3, 4], and ICARUS [5, 6]
detectors. The physics tolerance of the beam-related systematic uncertainty on
the neutrino measurement at NOvA detectors has been estimated [7]. Table 1
shows the range of beam parameters that stay within the physics tolerance. The
beam instrumentation to measure the beam parameters is also shown [8].

Table 1:  Estimated physics tolerance for the beam parameters of the NOvA experiment and the corresponding beam instrumentation. Here CT is a current transformer, BPM is a beam position monitor, and PM is a beam profile monitor. POT stands for the protons on target. Beam parameter | Design value | Tolerance | Beam instrument | Accuracy
---|---|---|---|---
Horn current | 200 kA |  $\pm$2 kA | Horn CT | 0.1%
Horiz. beam position on target | 0 mm |  $\pm$1 mm | BPM | 0.02 mm
Vertical beam position on target | 0 mm |  $\pm$1 mm | BPM | 0.02 mm
RMS beam spot size | 1.3 mm |  $\pm$0.2 mm | Beam PM | 0.1 mm
Beam intensity | 50$\times$1012 POT | 1% | Beam CT | 0.5% (stability 0.1%)

A set of conventional detectors is used as the primary beam instrumentation to
monitor the beam parameters with accuracy shown in Table 1. However, in the
event of a signal fluctuation in the detector, there is no method to check the
signal quality. Therefore, the presence of secondary beam detectors that
measure the beam parameters using a different method are a great advantage in
cross-checking measurements.

This paper shows that the NuMI horn magnet generates dipole and quadrupole
fields to focus pions. This suggests that the horn magnet optics is
predominantly linear. Therefore, the phase space evolution of pions and muons
in the NuMI beamline can be simulated using a linear beam optics model. The
three muon monitors downstream of a hadron absorber after the decay volume
(see Fig. 1) will be used as the secondary beam instrumentation. Each monitor
measures a profile of the muon beam produced as a result of pion decay
anywhere in the NuMI beam line. As Figure 1 shows, the thickness of the
absorber material in front and between the monitors is different. As a result,
each muon monitor detects muons of a different energy spectrum. Consequently,
the beam profile observed by the corresponding monitor can be attributed to
the magnet focusing strength at different pion energies. The observed profile
at each monitor responds to the variation of beam parameters independently. By
analyzing the muon monitor signal, the horn current is detected to $\pm$0.05%
accuracy. The other beam parameters in Table 1 are also detected from the
monitor signal with high accuracy. Cross-checking beam measurements using
different types of instrumentation and methods can effectively reduce beam-
related systematic uncertainties, making it a viable technique for future
neutrino experiments, such as Fermilab’s Long Baseline Neutrino Facility
(LBNF) [9, 10].

The main components of the NuMI beamline are introduced in Sec. II. The NuMI
horn focusing mechanism is investigated by using analytical and numerical
methods in Sec. III. The observed beam response at the muon monitors is
described in Sec. IV. The beam parameters are determined from the muon monitor
signal by using machine learning, which is demonstrated in Sec. V. Finally, we
conclude the study and discuss future prospects in Sec. VI.

##  II NuMI beam components

Figure 1 shows the layout of the NuMI target and horn system, including the
muon monitors. Every 1.2 seconds, a 120-GeV/$c$ proton beam is transported
from the Main Injector to the NuMI target system. The beam spill length is 9.6
$\mu$s. The highest beam intensity is 5.6$\times$1013 protons on target (POT)
per spill demonstrated in Spring 2022. The proton beam structure in the Main
Injector and the NuMI proton beam transport line optics are described in [8].
In the following subsections, we discuss the main beam components and their
relation to the study of the horn focusing mechanism.

Figure 1:  The overall layout of the NuMI beamline.

###  II.1 NuMI target and horn system

Figure 2 shows the layout of the baffle, target canister, and Horn 1 in the
NuMI beamline. The baffle is a 1.5-meter-long graphite cylinder with a 15-mm-
diameter center hole along the beam axis to capture any mis-steered beam. The
target canister is a water-cooled aluminum cylinder filled with atmospheric-
pressure Helium gas that contains the target core. A thin Beryllium foil is
used at both ends of the canister to make upstream and downstream beam
windows. The thickness of the downstream window in which the beam particles
pass through is 0.25 mm.

Figure 2:  The layout of the baffle, target chase, and Horn 1.

Figure 3 shows the target core made of a series of graphite fins (grade: POCO
ZXF-5Q, density: 1.78 g/cm3). The nominal beam position at the target core is
near the upper end of the fin. The first two fins are used for the beam
position measurement. They can be electrically isolated from the ground so
that the beam position can be measured by observing the intensity of the
delta-rays knocked out from the fins by the incident beam. The next four fins
have a thick graphite cylinder around the top of the fin. The cylinder
diameter, 18 mm, is slightly larger than the inner center hole diameter of the
baffle. These fins are called winged fins. A mis-steered proton beam hits the
winged fin to reduce the beam intensity produced by multiple scattering. This
way, the beam makes less thermal impact on the beam elements downstream from
the target, which are not designed to receive the full intensity of the proton
beam. The rest of the target core has 44 rounded rectangular target fins. The
length and width of the fins are 25 and 9 mm, respectively. The dimensions of
the fin are determined from a balance between the pion production yield and
thermal stress on the fin. The edge of the fin is rounded to avoid a shock
wave concentration at a corner. The present target core is capable of
operating with a 1-MW proton beam. The target core length is 1.25 m, and it
has 48+2 fins along the beam direction.

Figure 3:  The layout of the target core.

Figure 4 shows the elevation view of the target hall. There are two horns
called Horn 1 and Horn 2, respectively.

Figure 4:  The elevation view of the target hall. The beam axis is tilted by
58 mrad downward to direct the beam toward the Soudan Mine in Minnesota, where
the MINOS far detector was located.

Each horn is 3 m long and consists of two coaxial aluminum conductors. Figure
5 is the cross-sectional view of Horn 1. The inner conductor has a parabolic
shape at each end of the horn, and there is a narrow straight cylindrical tube
(called the “neck”) to connect the two parabolic shape conductors at both
parabolic vertices; the outer conductor is a straight cylinder. The upstream
end of the horn has a curved conductor that connects the inner and outer
conductors continuously. A stripline conductor is connected at the downstream
end of the horn to inject and extract the current. Horns 1 and 2 are
electrically connected in series.

Figure 5:  The cross-sectional view of Horn 1. A red arrow near the horn
conductor shows a horn current flow for the FHC mode. Trajectories of pions
with different initial momenta. Red: 5 GeV/$c$, green: 10 GeV/$c$, blue: 20
GeV/$c$, black: wrong-sign pion of 10 GeV/$c$ is defocused by the horn running
in the FHC mode.

A Forward Horn Current (FHC) mode or neutrino mode occurs when the horn
current flows in the inner conductor in the same direction as the beam. In
this mode, a toroidal magnetic field is induced between the inner and outer
conductors in the clockwise direction, focusing positively charged pions (Fig.
5). Negatively charged pions are focused in the Reverse Horn Current (RHC)
mode, also known as the antineutrino mode. In this case, the horn current
flows in the inner conductor in the opposite direction of the beam. To check
the quality of the magnetic field in the horn, the radial magnetic field is
measured at the neck before the horn is exposed to the beam. The observed
magnetic field strength agrees with the analytical field calculation to within
0.1% [11]. Because the high-energy pions tend to travel in the forward
direction of the proton beam, they often traverse inner diameter of the neck
where no magnetic field exists. Some of them are wrong-charge pions since they
generate wrong chirality neutrino events in the neutrino detectors. In
addition to that, a small fraction of electron neutrinos produced in the
decays of kaons and muons are also contained in the neutrino flux. The overall
contamination of the neutrino beam is a few percent.

The inner conductor is cooled by a water jet sprayed vertically from the outer
conductor. Atmospheric-pressure Argon gas is contained between the inner and
outer conductors. Argon gas is used because it has a higher electrical
breakdown voltage than Helium gas.

The radial size of Horn 2 is larger than that of Horn 1. Detailed dimensions
are shown in Sec. III.2. The distance between the two horns is 20 m. The
upstream end of the decay pipe is located 25 m downstream of Horn 2. The decay
pipe contains an atmospheric-pressure Helium gas, and most pions decay in the
pipe. The length and inner diameter of the pipe are 675 m and 2 m,
respectively, and there are 1.6-mm-thick Al windows on both ends of the pipe
for sealing.

The first beam monitor, known as the hadron monitor, is located at the
downstream end of the decay pipe. Presently, the hadron monitor is mainly used
for the proton beam-based alignment [8, 12]. The positions of the target and
horns are identified via proton beam tomography. The accuracy of the beam-
based alignment is typically $\pm$0.1 mm. It is consistent with the optical
survey result performed whenever a beam element is installed or replaced in
the beamline.

Figure 6 shows the schematic drawing of the absorber, which is located right
after the hadron monitor. It is essentially a box approximately 5.5 m $\times$
5.6 m $\times$ 8.5 m (W$\times$H$\times$L). The core volume is composed at its
upstream end of eight 1.3 m $\times$ 1.3 m $\times$ 0.3 m
(W$\times$H$\times$L) aluminum blocks followed by ten 1.3 m $\times$ 1.3 m
$\times$ 0.23 m (W$\times$H$\times$L) flame-cut steel blocks. Each aluminum
block has two independent water circuits, while the steel layers are not
water-cooled. The function of the steel blocks is to catch the tails of the
hadronic showers created in the core. Radiation shielding is provided by 88
steel blocks surrounding the absorber core and absorbing thermal and low-
energy neutrons, which iron is relatively transparent to. This double-layer
absorber structure is shown in the simulated muon flux in the muon monitor
(Fig. 7).

Figure 6:  The layout of the NuMI beam absorber. Figure 7:  Simulated 2D
histogram of the muon spectrum. The horizontal axis is the muon energy at
production, and the vertical axis is the muon energy at Muon Monitor 1. The
thick band is the muons passing through the Aluminum core, while the second
thin band is the muons passing through the steel block outside the Al core.
Muons lose more energy in steel than Al.

There are three beam monitors after the absorber. Most charged particles
produced upstream of the absorber are removed in the absorber volume except
for muons. Because these monitors mostly detect muons, they are called muon
monitors. The first (the most upstream) muon monitor is called Muon Monitor 1
(MM1), the second Muon Monitor 2 (MM2), and the last Muon Monitor 3 (MM3). A
detailed description of the muon monitors is given in Sec. II.2.3. All
residual charged particles, including muons, are stopped in the rock
downstream of the muon monitors. Only neutrinos pass through the rock.

###  II.2 NuMI Beam Instrumentation

The primary beam instrumentation and the muon monitors for this study are
described in this section.

####  II.2.1 Proton beam monitor

A pair of induced-charge pickup electrodes and a multi-wire Secondary Emission
electron Monitor (SEM) are used as a beam position monitor (BPM) and beam
profile monitor (PM), respectively. A special PM, which is permanently placed
upstream of the target, is made of a thin (5$\mu$meter thick $\times$ 25
$\mu$meter wide) Ti foil to mitigate the beam scattering at the foil. There
are $47\times 47$ foils stretched in $x$ and $y$ directions with gaps between
foils of 0.5 mm. The PM provides the proton beam profile for every spill.
Other PMs are periodically inserted in the beam line at different locations
for monitoring the beam optics.

Horizontal and vertical BPMs are paired to measure the beam position on the
transverse plane. Two sets of paired BPMs near the target are used for
measuring the proton beam position and beam orientation at the target. Several
BPMs comprise a network system to automatically correct the proton beam
position and beam orientation at the target, which is called the autotune [8].
The accuracy of the beam position measurement is $\pm$0.02 mm using BPM and PM
[8], and the accuracy of the RMS beam spot size measurement is $\pm$0.1 mm
using PM.

####  II.2.2 Current monitor

The proton beam intensity at the target is measured by using a toroidal beam
Current Transformer (beam CT) located upstream of the target. The signal
calibration is done by using a current source and by comparison to a Direct-
Current Current Transformer (DCCT) at the Main Injector. The accuracy of the
beam CT is 0.5% while the stability of signal pulse-per-pulse is $\pm$0.1% [8,
13].

The horn power supply circuit generates a half-sinusoidal voltage, and the
pulse is propagated through four striplines. The total pulse length of this
half sinusoidal wave is 2.334 msec. A current transformer measures the horn
current for each stripline (horn CT). The current transformer is calibrated
using a pre-calibrated current source. The accuracy of horn CT is $\pm$0.1%
[14].

####  II.2.3 Muon monitors

Figure 8 shows a schematic drawing of the muon monitor. Each muon monitor
consists of a 9$\times$9 grid of ionization chambers.

Figure 8:  Muon monitor engineering drawing showing the nine tubes containing
a row of nine pixels each.

Each chamber has a $7.5\times 7.5$ cm2 sense area Ag-Pt plated electrode on a
1-mm-thick Alumina ceramic electrode base, and the gap between the electrodes
is 3 mm. The distance between the centers of channels is 250 mm to cover a
$2.1\times 2.1$ m2 area. Nine ionization chambers are installed in a
rectangular aluminum tube. There are nine tubes assembled vertically on a
support structure in the beam line. Atmospheric-pressure 99.995% pure Helium
gas is used as an ionization media, injected from the top and exhausted from
the bottom of each tube. The bias voltage is 300 V in which the ionization
chamber is on the plateau where all ionized charged particles in the chamber
are collected at the electrondes without charge multiplication [15, 16].

Figures 9–11 show an example of the observed pixel image per spill. The
pedestal signal, taken when the beam is turned off, is subtracted from the
beam-on signal. The size of the pedestal signal is on the order of 0.1% as
compared to the beam signal. The signal for each channel was calibrated using
a radiation source before the monitor was installed in the beam line. The
integrated signal from all channels is normalized by the beam intensity per
beam spill and is monitored during beam operation. The long-term fluctuation
of the integrated signal normalized by the beam intensity is $\pm 0.2$% for
Muon Monitors 1 and 2 while the fluctuation of normalized Muon Monitor 3
signal is slightly higher at $\pm 0.5$%. It is found that the observed pixel
image pattern is not changed by the beam intensity although the signal size is
proportional to the beam intensity.

Figure 9:  The observed beam profile at Muon Monitor 1. The top and right
subplots show the projection of the beam profile in horizontal and vertical
planes, respectively.  Figure 10:  The observed beam profile at Muon Monitor
2. The top and right subplots show the projection of the beam profile in
horizontal and vertical planes, respectively.  Figure 11:  The observed beam
profile at Muon Monitor 3. The top and right subplots show the projection of
the beam profile in horizontal and vertical planes, respectively.

The muon monitors are located downstream of the hadron absorber and separated
by 12 and 18 m of rock (see Fig. 1). The muons reaching the monitors lose part
of their kinetic energy in the material via energy loss processes. Therefore,
the initial spectrum of muons when they are produced differs from the spectrum
at the monitors. Figure 12 suggests that each monitor detects a different muon
phase space with a specific cutoff energy.

Figure 12:  Simulated muon spectra. A solid red line represents the muon
spectrum when muons are produced while a dashed red line is the muon spectrum
at Muon Monitor 1. The shift is because muons are lost energy in the absorber.
The green and blue lines are for Muon Monitor 2 and 3, respectively.

It is demonstrated in Sec. III that the observed beam profile on the monitors
is the result of different horn magnet focusing strengths for different pion
energies. Figures 13–15 show the simulated phase space of the initial pions
for which the daughter muons reach the corresponding muon monitor. Muon
Monitor 1 detects the initial pion momentum range above 5 GeV/$c$, Muon
Monitor 2 detects the pion momentum range above 12 GeV/$c$, and Muon Monitor 3
detects the pion momentum above 25 GeV/$c$. It should be noted that the
observed muon profile on the monitor is a result of the multiple Coulomb
scattering in matter. As a result, muon signals are not correlated with
specific neutrino events.

Figure 13: Simulated transverse and longitudinal momenta of the pions for
which the corresponding muons reach Muon Monitor 1. It should be noted that
the low statistic bin is cut and set to be zero in the plot. Figure 14:
Simulated transverse and longitudinal momenta of the pions for which the
corresponding muons reach Muon Monitor 2. Figure 15: Simulated transverse and
longitudinal momenta of the pions for which the corresponding muons reach Muon
Monitor 3.

##  III Horn focusing mechanism

To maximize the sensitivity of neutrino oscillation and CP violation
measurements, modern focusing horn systems are designed using numerical
simulations. However, analyzing the horn focusing mechanism is challenging due
to the complexity of the optimized system. For example, the target is not a
point source, and the pion phase space received in the horn is broad.
Additionally, the optimized horn itself has a complicated shape. These issues
make it difficult to accurately model the particle trajectories. To better
understand the mechanism, we start by using an analytical model formulated
with a conventional beam optics method. A second-order polynomial function is
employed to describe the mechanism in this model.

Next, we execute a simple simulation in which pions are ray-traced through the
two horns using an analytical formula. This is therefore a semi-analytical
model. By applying various assumptions in both models, the horn focusing
strength is generalized and we can evaluate the acceptance of horn system
without executing numerical simulations. However, we lose the accuracy of beam
dynamics. The analytical and semi-analytical models are evaluated with a
numerical ray-trace simulation in G4beamline [17], in which pions are tracked
in a realistic horn field and horn conductor material. The mean energy loss is
taken into account in the pion tracking in matter. However, the simulations do
not include stochastic processes such as decay, hadronic interaction, or
multiple scattering.

Finally, the full model simulation is executed by using the modified G4NuMI
framework to simulate all stochastic processes, including decay, hadronic
interactions, and Coulomb interactions in the whole NuMI beamline. The
acceptance of the NuMI horn system is evaluated by comparing the numerical
simulation to the semi-analytical model. It is found that those results are in
good agreement.

###  III.1 Analytical model

A horn conductor generates a toroidal magnetic field between the inner and
outer conductors. As a result, an axially symmetric Lorentz force acts on
charged particles. We use a cylindrical coordinate system ($r$, $z$, $\phi$).
A particle trajectory in a long coaxial transmission line has been studied
[18] in which the particle follows a spiral orbit in the space available
between inner and outer conductors. In this paper, the angular momentum is
ignored ($\phi^{\prime}\sim 0$) in the model because the NuMI horn is too
short to make a periodic motion around the horn inner conductor. A similar
treatment is illustrated in Section 5.2.2 in the textbook Particle Accelerator
Physics by Helmut Wiedemann.[19]

We start from the equation of motion in the transverse direction,

$\frac{dp_{r}}{dt}=\frac{dp_{r}}{dz}\frac{dz}{dt}=\frac{dp_{r}}{dz}v_{z}=-qv_{z}b_{\phi},\\\
$ (1)

where $p_{r}$ is the transverse momentum, $v_{z}$ is the longitudinal
velocity, and $b_{\phi}$ is the toroidal magnetic field. The equation is
simplified by using $\displaystyle{\frac{dp_{r}}{dz}=-qb_{\phi}}$, or
$\displaystyle{dp_{r}=-qb_{\phi}dz}$. By inserting $dp_{r}\equiv
p_{r}-p_{r,0}$ or $p_{r}=p_{r,0}+dp_{r}$ and $\displaystyle{p_{r}=\gamma
m\frac{dr}{dt}=\gamma m\frac{dz}{dt}\frac{dr}{dz}=p_{z}\frac{dr}{dz}}$, where
$p_{z}$ is the longitudinal momentum, the equation is transformed:

$\displaystyle p_{z}\frac{dr}{dz}=p_{r,0}-qb_{\phi}dz,$ (2)
$\displaystyle\frac{dr}{dz}=\frac{p_{r,0}}{p_{z}}-\frac{qb_{\phi}dz}{p_{z}}.$
(3)

Here $p_{r,0}$ is the initial transverse momentum (and we assume $p_{z,0}\sim
p_{z}$). By replacing $\displaystyle{\frac{dr}{dz}\rightarrow\theta}$ and
$\displaystyle{\frac{p_{r,0}}{p_{z}}\rightarrow\theta_{0}}$, the equation
becomes

$\theta=\theta_{0}-\frac{qb_{\phi}dz}{p_{z}}.$ (4)

Eq. (4) shows the evolution of the angle of a charged particle in a horn
field. The horn toroidal field is

$b_{\phi}=\frac{\mu_{0}I}{2\pi r}.$ (5)

Thus, in Eq. (4) with a constant horn current, $\theta$ is a function of
$p_{z}$, $r$, and $dz$:

$\theta(p_{z},r,dz)=\theta_{0}-\frac{q\mu_{0}I}{2\pi}\frac{dz}{p_{z}r}.$ (6)

$dz$ is the projection of pass length of a charged particle along $z$ axis. In
the following section, it is found that the NuMI magnetic horn works as a beam
focusing element. The length of the horn is relatively short by comparing with
the distance from the horn to the monitors. Therefore, a thin lens
approximation is applied in the analytical model in which $dz$ loses the
information about $z$ position, but is represented as a function of $r$. In
the textbook, $dz\sim ar^{2}$ is used for presenting the desired focusing
properties of the magnetic horn with a constant magnetic field strength, where
$a$ is a quadratic coefficient. To generalize the focusing properties, we
assume that $dz$ is a polynomial function of $r$:

$dz\equiv a_{0}+a_{1}r+a_{2}r^{2}+a_{3}r^{3}+O(r^{4}),$ (7)

where $a_{0},a_{1},a_{2}$ and $a_{3}$ are constants. $a_{0}$ is set to be zero
in the model to avoid a singularity at $r=0$. Hence,

$\frac{dz}{r}=a_{1}+a_{2}r+a_{3}r^{2},$ (8)

where the first, second, and third terms represent the expansion’s dipole,
quadrupole, and sextupole components, respectively. The particle angle after
the transverse kick $\theta$ is

$\theta(p_{z},r,\theta_{0})=\theta_{0}-\frac{q\mu_{0}I}{2\pi
p_{z}}\left(a_{1}+a_{2}r+a_{3}r^{2}\right).$ (9)

Eq. (9) suggests that the horn acts as a multipole magnet. The distance of the
charged particle from the horn to the crossing of the beam axis is

$\lambda(p_{z},r,\theta_{0})=\frac{r}{\tan\theta}\sim\frac{r}{\theta}=\frac{r}{\theta_{0}-\frac{q\mu_{0}I}{2\pi
p_{z}}\left(a_{1}+a_{2}r+a_{3}r^{2}\right)}$ (10)

Thus, the particle direction is parallel to the beam axis when
$\lambda\rightarrow\pm\infty$ (or $\theta\rightarrow\pm 0$). Figure 16 shows
the schematic drawing of the geometry used in the analytical model.

Figure 16:  Schematic drawing of the geometry which is used in the analytical
model. “$\it{O}$” in the drawing is the horn position in the model.

###  III.2 Semi-analytical model

The NuMI beamline has two horns, and each horn’s shape is parabolic at both
ends of the horn (Fig. 5). However, the analytical model only represents a
single horn. A semi-analytical model is proposed to take into account the
geometrical complications. The transverse kick in the magnetic horn is modeled
by using the formula

$\theta_{k}=q\frac{\int b_{\phi}dl}{p_{z}},$ (11)

where $dl$ is the path length of a particle in the magnetic volume of the
horn. Again, we use the assumption that the trajectory of a particle in the
magnetic volume of the horn is parallel to the horn axis. Hence, the
integrated magnetic field along the path length $\int b_{\phi}dl$ is
analytically estimated by calculating the intersection points of a particle
with the horn conductors while the magnetic field strength $b_{\phi}$ is a
constant. The dimensions of the NuMI horns are summarized in Table 2.

Table 2: Magnetic horn profile for Horn 1 and Horn 2.  | z (mm) | profile (mm)
---|---|---
Horn 1 | 0.00 $\leq z<80.00$ | $\displaystyle{\sqrt{\frac{92.85-z}{7.05}}}$
| 80.00 $\leq z<83.98$ (neck) | 1.35
| 83.98 $\leq z<300.00$ | $\displaystyle{\sqrt{\frac{z-80.00}{2.19}}}$
Horn 2 | 0.00 $\leq z<97.62$ | $\displaystyle{\sqrt{\frac{100.00-z}{0.14}}-0.30}$
| 97.62 $\leq z<104.80$ (neck) | 3.90
| 104.80 $\leq z<300.00$ | $\displaystyle{\sqrt{\frac{z-100.00}{0.27}}-0.30}$

The outer radii of Horn 1 and Horn 2 are 174.6 and 395.4 mm, respectively.

###  III.3 Test linear optics model

The particle tracking simulation uses G4beamline [17] to test the linearity of
the NuMI horn optics and compare it with the analytical and semi-analytical
models. In the simulation, a horn has a conductor material, and a mangetic
field (given in Eq.(5)) is applied between inner and outer conductors where
atmospheric-pressure Ar gas is contained. Test particles (pion) are
distributed in front of Horn 1 (“$O$” in Fig. 16 is the upstream end of Horn 1
in the particle track simulation). The initial test particle position is
distributed with a 0.5 mm step in the radial direction (a vertical axis in
Fig. 16). The initial transverse and longitudinal momenta are variable
parameters to fill the possible initial pion phase space escaped from the
target.

Figure 17 shows the estimated $\lambda$ as a function of the initial radial
particle position $r$ in the particle track simulation (a blue point), semi-
analytical model (a solid red curve), and analytical model (a dashed green
curve) for a pion with $p_{z}=12$ GeV/$c$ and $p_{r,0}=0.2$ GeV/$c$. The semi-
analytical model reproduces the particle tracking simulation result at
$r>0.01$ m. Besides, the semi-analytical model reproduces the second
singularity point around $r=0.015$ m. This is caused by the transverse kick in
Horn 2 after crossing the beam axis. It is shown in the ray-trace in both
subplots in Fig. 17, where pions are kicked by Horn 2 around z = 20 m and $r$
= -0.05 $\sim$ -0.10 m in the semi-analytical model and the particle track
simulation. On the other hand, the semi-analytical model is not well
reproduced $\lambda$ at $r<0.01$ m because all test particles below $r<0.01$ m
in the semi-analytical model miss the horns due to the thin lens model, which
is not the case in the particle track simulation.

Figure 17: Simulated $\lambda$ as a function of the initial radial position
$r$ for a pion with $p_{z}=12$ GeV/$c$ and $p_{r,0}=0.2$ GeV/$c$. Blue points:
particle track simulation, solid red line: semi-analytical model, dashed green
line: analytical model. Top subplot: ray trace of pions with the semi-
analytical model, bottom plot: ray trace of pions using the particle tracking
simulation.

The coefficients of the analytical model are found from the fitting of the
particle track simulation result. These are $a_{1}=0.0061$, $a_{2}=0.069$, and
$a_{3}=0.0$. As the fitting is conducted at a singularity point, the fitting
coefficient is determined by minimizing the chi-square value. In this case,
the fitting error is considered to be significantly small and is therefore
omitted. However, we found that the sextupole coefficient has a wide range of
low chi-square values. This suggests that the sextupole coefficient is not
sensitive to changes in the $\lambda$ curve.

The semi-analytical model reproduces the complicated particle tracking in
higher momentum regions while it does not work well at lower momenta. Figure
18 shows the estimated $\lambda$ for a pion with $p_{z}=18$ GeV/$c$ and
$p_{r,0}=0.2$ GeV/$c$. There are two singularity points in this case. The
semi-analytical model precisely reproduces these points as well as the
complicate $\lambda$ curve at $r=0.016\sim 0.02$ m. The best fitting of the
analytical model is made at $a_{1}=0.0063$, $a_{2}=0.12$, and $a_{3}=2.7\times
10^{-4}$. It should be noted that pions at $r<0.01$ m miss Horn 1 while they
are focused by Horn 2 so $\lambda$ has a negative value. The focusing
mechanism is observed in both subplots.

Figure 18: Simulated $\lambda$ as a function of the initial radial position
$r$ for a pion with $p_{z}=18$ GeV/$c$ and $p_{r,0}=0.2$ GeV/$c$. Blue points:
particle track simulation, solid red line: semi-analytical model, dashed green
line: analytical model.

On the other hand, Figure 19 shows the estimated $\lambda$ for a pion with
$p_{z}=10$ GeV/$c$ and $p_{r,0}=0.3$ GeV/$c$. The estimated $\lambda$ in the
semi-analytical model is underestimated by comparison with the particle track
simulation. The transverse kick to a particle in the semi-analytical model is
found to be stronger than in the particle track simulation. As a result, the
estimated $\lambda$ in the semi-analytical model tends to be shorter. It is
also found that the estimated $\lambda$ for low momentum particle is sensitive
to the gaps between the target and Horn 1, and Horn 1 and Horn 2. The fitting
coefficients of the analytical model are $a_{1}=0.0061$, $a_{2}=0.048$, and
$a_{3}=1.5\times 10^{-5}$.

The fitting of the analytical model suggests that the singularity behavior is
reproduced well by the dipole and quadrupole components, while the sextupole
component is negligibly small. Consequently, the linear optics model is a good
approximation for the particle track in the NuMI horn system. It is found that
the coefficient of the quadrupole component is changed by the momentum. It
will represent a chromaticity of the horn. Further study is in progress.

Figure 19: Simulated $\lambda$ as a function of the initial radial position
$r$ for a pion with $p_{z}=10$ GeV/$c$ and $p_{r,0}=0.3$ GeV/$c$. Blue points:
particle track simulation, solid red line: semi-analytical model, dashed green
line: analytical model.

###  III.4 Acceptance study of the NuMI horn system

The acceptance of the NuMI horn system is executed in the full model
simulation. The modified G4NuMI is used. G4NuMI contains the most updated NuMI
beam elements for the neutrino flux simulation. The modified G4NuMI can track
muons and duplicate pion events multiple times to increase the statistics of
muon events. The decay and kinematics of daughter particles in the duplicated
pions are independently executed using the normal Geant4 processes.

The full model simulation starts with the primary protons impinging on the
target. The acceptance is defined by the initial phase space of the pions
right after escaping from the target, for which the decay products (muons)
reach the muon monitors. On the other hand, the initial beam in the semi-
analytical model is comprised of pions. The initial position of pions is
uniformly distributed on the surface of the target volume along the $z$ axis.
The initial transverse and longitudinal momenta are uniformly distributed.
Therefore, the semi-analytical model only evaluates the boundary of the
acceptance. The energy cut is taken into account in the semi-analytical model
for low momentum pions. Figure 20 shows the comparison of the full model
simulation where the acceptance is shown as a 2D histogram with the semi-
analytical model for which the boundary of the acceptance is outlined with a
red dashed line. Both simulations are in good agreement.

There is a mismatch in the high $p_{z}$ region. Most pions which have
$p_{z}=$25 GeV/c or higher can reach the absorber without decay. Most of those
pions are absorbed in the absorber. A small fraction of pions are converted
into high energy muons. The semi-analytical model does not take into account
the pion loss process in the absorber. Indeed, high momentum pions are shown
in the spectrum of Muon Monitor 2 (Fig. 14) which is low statistics by
comparing Muon Monitor 1. The acceptance of the semi-analytical model
represents well the combined spectrum of Muon Monitor 1 and 2. The estimated
$\lambda$ shows that low $p_{z}$ particles have a different tracking between
the semi-analytical model and the numerical simulation in the previous
subsection. Such a discrepancy is too small to affect the acceptance.

Figure 20: Full model simulation result (2D histogram) vs. semi-analytical
model (red dashed boundary outline).

The semi-analytical model is used to further develop our acceptance studies.
Figure 21 demonstrates how the acceptance of pions is influenced by their
initial position upon emission from the target. The blue area shows the
acceptance for which pions are emitted from the first 1/3 of the target
length. The green area shows the acceptance for which pions are emitted from
the middle 1/3 of the target length. The orange area shows the acceptance for
which pions are emitted from the last 1/3 of the target length. Higher-
momentum pions are more likely to be emitted from the first third of the
target length, while lower-momentum pions are predominantly emitted from the
last third of the target length. It suggests that the time domain muon profile
measurement at the muon monitors shows the phase space of pions produced in a
specific target section. The early (late) arriving muons at the monitor will
be the result of the decay of high (low) energy pions which could be produced
at the last (first) third of target. We are considering using a picosecond
detector to measure the time domain muon profile.

Figure 21: Simulate the boundary of the acceptance of pions in the semi-
analytical model. The acceptance of pions separated by the start position of
pions within the target. Blue: pions start from the first 1/3 of the target
length. Green: pions start from the middle 1/3 of the target length. Orange:
pions start from the last 1/3 of the target length. Higher-momentum pions tend
to start from the first 1/3 of the target length, while lower-momentum pions
are extracted mostly from the last 1/3 of the target length.

##  IV Beam response measurement

The first systematic beam response measurement for studying the horn focusing
mechanism was carried out on December 12, 2019, by changing the beam position
at the target and the horn current. The procedure is as follows:

  * •

Place the horn in the FHC mode.

  * •

Set the beam intensity lower than $4\times 10^{13}$ POT per beam spill.
Limiting the beam intensity protects the target from a mis-steered beam.

  * •

Move the proton beam position at the target in the horizontal and vertical
directions. The beam orientation is kept parallel to the horn axis. The
observed target center from the proton tomography measurement was (horizontal,
vertical) = (0.0, -0.8) mm. The movement range of the beam position is
typically $\pm$0.8 mm from the target center.

  * •

Observe the beam centroid variation on each muon monitor as a function of the
proton beam position at the target.

  * •

Once the beam scan data was taken, vary the horn current. The horn current was
set to 200, 195, 190, and 180 kA.

Figure 22 is the observed horizontal beam centroid at the muon monitors as a
function of the horizontal proton beam position on target at a horn current of
200 kA. In order to establish linearity, two types of fitting are used. The
solid line is a linear fit and the dashed line is a fit with a third-order
polynomial function. The fitting analysis shows that the horizontal beam
centroid variation is linear. The variations of the intercepts at different
monitors are caused by the imperfection of alignment for the beam element.
They are not discussed in this study.

Figure 22: The observed horizontal beam centroid at the muon monitor as a
function of the horizontal proton beam position at the target. The horn
current is 200 kA. Blue, orange, and green points are the observed beam
centroids at Muon Monitors 1, 2, and 3, respectively. The solid line
represents a linear fit, and the dashed line represents a third-order
polynomial fit.

Muon Monitor 2 data are particularly interesting. The slope is almost zero. It
suggests that the beam is in focus at Muon Monitor 2. A negative slope at Muon
Monitor 1 shows that the beam is overfocused at Muon Monitor 1. In contrast, a
positive slope at Muon Monitor 3 shows that the beam is underfocused at Muon
Monitor 3, or the beam diverges. Thus, each muon monitor detects muons focused
differently based on the parent pion energy. These phenomena are consistent
with the linear optics model.

Figure 23 is the observed vertical beam centroid at the muon monitors as a
function of the vertical proton beam position on the target at a horn current
of 200 kA. The fitting analysis suggests that the observed beam centroid
variation is still predominantly linear. However, a nonlinear component
clearly appears, especially in the Muon Monitors 2 and 3. It should be noted
that the target geometry is asymmetric in the vertical direction, so the
variation of the vertical pion phase space should differ from the horizontal.
This is evidenced by the results of the vertical beam scan. The nonlinear
component is reproduced in the full model simulation. Overall, this is a good
validation test of the physics processes in the full model simulation. A
detailed study of validating the physics model will be presented in another
paper.

Figure 23: The observed vertical beam centroid at the muon monitor as a
function of the vertical proton beam position at the target. The horn current
is 200 kA. Blue, orange, and green points are the observed beam centroids at
Muon Monitors 1, 2, and 3, respectively. The solid line represents a linear
fit, and the dashed line represents a third-order polynomial fit.

Figures 24 and 25 show the fitting curves of the observed horizontal and
vertical beam centroid variations at different horn currents, respectively.

Figure 24: The observed vertical beam centroid at the muon monitor as a
function of the vertical proton beam position at the target at horn current
200, 195, 190, and 180 kA. Blue, orange, and green lines are the fits for beam
centroids at Muon Monitors 1, 2, and 3, respectively. The solid line
represents a linear fit, and the dashed line represents a third-order
polynomial fit.  Figure 25: The observed vertical beam centroid at the muon
monitor as a function of the vertical proton beam position at the target at
horn current 200, 195, 190, and 180 kA. Blue, orange, and green lines are the
fits for beam centroids at Muon Monitors 1, 2, and 3, respectively. The solid
line represents a linear fit, and the dashed line represents a third-order
polynomial fit.

Figures 26 and 27 show the measured slope for horizontal and vertical beam
scans from the fitting analysis as a function of the horn current,
respectively. A solid line is a result of linear fitting of the measured slope
while a dashed line is the simulated slope by using the full model simulation.
It is found that the slope is proportional to the horn current in the
measurement and the simulation. The linear horn current dependence is
predicted by the analytical model in Eq. 9.

Overall, the full model simulation reproduces the measured slope well.
However, the simulation does not include two background processes [15]. A
delta-ray which is generated in the material by interacting with high energy
muons affects the muon monitor signal. Elastically scattered protons from the
target produce additional particles in the absorber. Those background
processes are included in a new simulation study. Fine tuning of the full
model simulation does not affect the beam centroid result, but changes the
muon beam profile. Further simulation study of these effects is in progress.

Figure 26: The estimated slope from the horizontal beam scan as a function of
the horn current. Blue, orange, and green points are the observed beam
centroids at Muon Monitors 1, 2, and 3, respectively. Blue, orange, and green
dashed lines are the simulated beam centroid at Muon Monitors 1, 2, and 3 in
the full model simulation, respectively.  Figure 27: The estimated slope from
the vertical beam scan as a function of the horn current. Blue, orange, and
green points are the observed beam centroid at Muon Monitors 1, 2, and 3,
respectively. Blue, orange, and green dashed lines are the simulated beam
centroid at Muon Monitors 1, 2, and 3 in the full model simulation,
respectively.

##  V Machine Learning Algorithm to Detect Beam Parameters using Muon Monitor
Data

As demonstrated in the previous sections, the NuMI horn is a linear beam
optics element. It suggests that the muon monitor response signal should be
linearly related to the beam parameter variations. The muon monitor signal
analysis is carried out by using a Machine Learning (ML) algorithm. An
advantage of using ML for this application is that it finds correlations among
a large number of input parameters and deconvolves them to make a prediction.
For example, it is known that the beam intensity affects the horn current
because the beam heats the horn conductor by the energy deposition and changes
its conductance. By applying ML, the horn current and the beam intensity are
detected independently from the observed muon profile.

An Artificial Neural Network (ANN) with multiple hidden layers is used to
build the ML architecture. The ANN takes 241 observed values from individual
muon monitor channels as inputs. It predicts four beam parameters: horn
current, beam intensity, and horizontal and vertical beam positions at the
target. Two dead channels on Muon Monitor 1, (row, column) = (1,5) and (4,1)
in Fig. 9, are omitted in this analysis. The total number of data points
prepared for the study is $13,563$. The data is divided into training (70%)
and validation (30%) data to optimize the ANN model architecture.

Figure 28 compares the measured horn current by using the horn current
transformer (CT) and the predicted horn current by the trained ANN model.

Figure 28: Comparison of the predicted horn current from the muon monitor
signal vs. the measured horn current by using the horn CT.

On the other hand, Figure 29 compares the beam current measured by the beam CT
and the predicted beam intensity by the trained ANN model.

Figure 29: Comparison of the predicted beam intensity from the muon monitor
signal vs. the measured beam intensity using the beam CT. The units are
$10^{12}$ protons per spill.

These figures show that the horn current and the beam intensity are identified
independently from the same validation data by applying the trained ANN model.
It validates that ML can simultaneously detect the horn current and beam
intensity from the observed muon profile.

Figures 30 and 31 show the measured and detected horn current and beam
intensity as a function of time, respectively. The fluctuations observed in
the signal of the primary beam instrumentation can be reproduced by the values
detected using a trained ANN model. It suggests that the fluctuations are
unlikely to be due to electrical noise but rather are likely to be associated
with variations in the beam intensity and horn current stability. The accuracy
of the estimated beam parameters from the ANN is $\pm$0.05% for the horn
current and $\pm$0.1% for the beam intensity.

Figure 30: Comparison between the observed horn current by using the horn CT
and the predicted horn current from the muon monitor signal by using ML.
Figure 31: comparison between the observed beam intensity by using the beam CT
and the predicted beam intensity from the muon monitor signal by using ML.

Figures 32 and 33 compare the horizontal and vertical beam positions at the
target measured by using the BPM and detected by the trained ANN model,
respectively. Again, the detection values reproduce the measurements well. It
suggests that the signal fluctuation is the stability of the beam position at
the target. The accuracy of ML detection is $\pm$0.018 mm and $\pm$0.013 mm
for the horizontal and vertical positions, respectively.

Figure 32: Comparison between the horizontal beam position at the target
observed by using the BPM and the horizontal proton beam positions predicted
by using the ML. Figure 33: Comparison between the vertical beam position at
the target observed by using the BPM and the vertical proton beam positions
predicted by using the ML.

Table 3 shows that the muon monitor can detect all beam parameters affecting
the NOvA experiment’s physics tolerance. The accuracy of the Machine Learning
(ML) detection from the muon monitor signal is similar to the accuracy of the
primary beam instrumentation, which is shown in Table 1. This implies that the
ML method is able to provide reliable results that are comparable to the
traditional methods of measuring the primary beam instrumentation.

Table 3:  Beam instrumentation and beam parameters monitored by the instrumentation. An open circle marks the parameter that can be detected using ML, while a triangle marks the parameter that will be tested in the future.  Beam parameter | BPM | PM | Beam CT | Horn CT | MM
---|---|---|---|---|---
Beam position | $\bigcirc$ | $\bigcirc$ |  |  | $\bigcirc$
Beam spot size |  | $\bigcirc$ |  |  | $\triangle$
Beam intensity |  |  | $\bigcirc$ |  | $\bigcirc$
Horn current |  |  |  | $\bigcirc$ | $\bigcirc$

From these studies, we confirm that ML has the potential to improve the
ability to diagnose the NuMI target and horn system. One specific application
of ML is to use it as an anomaly event detector, which can identify anomalies
such as target density deterioration, target core damage, or degraded horn
field quality. However, the accuracy of the ML detection is highly dependent
on the quality of the training dataset. Therefore, it is important to optimize
the ML architecture and reduce systematic errors in the training data to
improve the accuracy of the model. These studies are underway.

##  VI Conclusion and Future Prospects

The NuMI horn focusing mechanism has been investigated by analyzing the
observed muon monitor signals in the beam scan. The results suggest that the
horn focusing magnet is a linear optics element, which has been confirmed by
both analytical models and numerical simulations. Moreover, we have employed
machine learning to deconvolve the horn current and beam intensity from the
observed muon profile. The trained ANN model has demonstrated high accuracy in
detecting the horn current and the beam intensity within $\pm$0.05% and
$\pm$0.1% accuracy, respectively. Additionally, the ML approach has enabled us
to detect the beam position at the target with an accuracy of $\pm$0.018 mm
and $\pm$0.013 mm for the horizontal and vertical positions, respectively,
from the observed muon profile. As expected, the detection accuracy is
consistent with the accuracy of the primary beam instrumentation. Our current
focus is on enhancing the ML’s performance by preparing a training dataset
with low systematic uncertainty and developing an algorithm to detect
anomalous events.

The results of this study suggest that the muon monitor can be used as the
secondary beam instrumentation to monitor the beam parameters related to the
NOvA physics tolerance with high accuracy. Muon monitors provide a significant
benefit when validating measurements. Additionally, the study results indicate
that Fermilab’s Long Baseline Neutrino Facility (LBNF) might benefit from this
technique.

Other secondary detectors, such as beam detectors located near the target, are
highly effective in directly detecting secondary particles escaping from the
target. However, it is worth noting that these detectors may not be reliable
under extreme radiation conditions. Hence, upgrading the muon monitors should
be considered to extend their functionality and reliability. Because the horn
is a linear optics element, the muon profile on the muon monitor results from
the chromaticity of the horn magnets. The chromaticity measurement can be
obtained by observing the time-of-flight (TOF) of muons at various muon
monitor channel locations. The correlation can be replicated in the full model
simulation, and the initial pion energy spectrum will be reconstructed from
the observed time domain muon profile. Ultimately, these are a resource to
reconstruct the initial pion phase space.

It is found that the full-model simulation closely replicates the beam scan
result. However, the observed muon profile is not accurately reproduced by the
simulation. The fine tuning is required in simulation. Two background signals
on the muon monitors, namely the delta rays from high energy muons and excess
particles creating by the elastically scattered protons interacting in the
absorber, must be evaluated to ensure accurate simulation. The hadronic
interaction model needs to be tuned to replicate the measurements. Those
studies are currently underway and will be published in separate papers.

###### Acknowledgements.

We would like to thank the Fermilab Accelerator Operations, Target Systems,
External Beam Delivery, and Main Injector Departments for their contributions
to this study. We also want to thank the NOvA collaboration for providing us
with an opportunity to run the special beam tests. This work is supported by
the Fermi Research Alliance, LLC manages and operates the Fermi National
Accelerator Laboratory pursuant to Contract number DE-AC02-07CH11359 with the
United States Department of Energy. This work is partially supported by the
U.S. Department of Energy grant DE-SC0019264.

## References

  * Acero  _et al._ [2019] M. A. Acero _et al._ (NOvA), First Measurement of Neutrino Oscillation Parameters using Neutrinos and Antineutrinos by NOvA, Phys. Rev. Lett. 123, 151803 (2019), arXiv:1906.04907 [hep-ex] .
  * Adamson  _et al._ [2016a] P. Adamson _et al._ (NOvA), First measurement of electron neutrino appearance in NOvA, Phys. Rev. Lett. 116, 151806 (2016a), arXiv:1601.05022 [hep-ex] .
  * Abratenko  _et al._ [2021] P. Abratenko _et al._ (MicroBooNE), Measurement of the flux-averaged inclusive charged-current electron neutrino and antineutrino cross section on argon using the NuMI beam and the MicroBooNE detector, Phys. Rev. D 104, 052002 (2021), arXiv:2101.04228 [hep-ex] .
  * Acciarri  _et al._ [2017] R. Acciarri _et al._ (MicroBooNE), Design and Construction of the MicroBooNE Detector, JINST 12 (02), P02017, arXiv:1612.05824 [physics.ins-det] .
  * ICARUS Collaboration [2023] ICARUS Collaboration, Icarus at the fermilab short-baseline neutrino program – initial operation (2023).
  * Farnese  _et al._ [2023] C. Farnese _et al._ (ICARUS), Implementation of the trigger system of the ICARUS-T600 detector at Fermilab, Nucl. Instrum. Meth. A 1045, 167498 (2023).
  * Cremonesi [2017] L. Cremonesi, NOvA-DOCDB NOVA Document 17608-v5 (2017) internal technical note.
  * Adamson  _et al._ [2016b] P. Adamson _et al._ , The NuMI Neutrino Beam, Nucl. Instrum. Meth. A 806, 279 (2016b), arXiv:1507.06690 [physics.acc-ph] .
  * Acciarri  _et al._ [2016] R. Acciarri _et al._ ,  _Long-Baseline Neutrino Facility (LBNF) and Deep Underground Neutrino Experiment (DUNE): Conceptual Design Report, Volume 1: The LBNF and DUNE Projects_ , Tech. Rep. (DUNE, 2016) arXiv:1601.05471 [physics.ins-det] .
  * Acciarri  _et al._ [2015] R. Acciarri _et al._ ,  _Long-Baseline Neutrino Facility (LBNF) and Deep Underground Neutrino Experiment (DUNE): Conceptual Design Report, Volume 2: The Physics Program for DUNE at LBNF_ , Tech. Rep. (DUNE, 2015) arXiv:1512.06148 [physics.ins-det] .
  * Hylen [2009] J. Hylen, MINOS-DOCDB MINOS Document 6122-v1 (2009) internal technical note.
  * Zwaska  _et al._ [2006] R. Zwaska _et al._ , Beam-Based Alignment of the NuMI Target Station Components at FNAL, Nucl. Instrum. Meth. A 568, 548 (2006), arXiv:physics/0609106 .
  * Crisp  _et al._ [2011] J. Crisp, B. Fellenz, J. Fitzgerald, D. Heikkinen, and M. A. Ibrahim, Operation of the intensity monitors in beam transport lines at fermilab during run II, Journal of Instrumentation 6 (10), T10001.
  * Hylen and Yonehara [2018] J. Hylen and K. Yonehara, Beams-DOCDB Beams Document 9520-v1 (2018) internal technical note.
  * Loiacono [2010] L. J. Loiacono,  _Measurement of the Muon Neutrino Inclusive Charged Current Cross Section on Iron Using the MINOS Detector_ , Ph.D. thesis, University of Texas Austin (2010).
  * Zwaska [2005] R. M. Zwaska,  _Accelerator Systems and Instrumentation for the NuMI Neutrino Beam_ , Ph.D. thesis, University of Texas Austin (2005).
  * [17] https://www.muonsinc.com/Website1/G4beamline .
  * Regenstreif [1961] E. Regenstreif, Technical Report 64-41 (CERN, 1961) internal technical note.
  * Wiedemann [2015] H. Wiedemann,  _Particle Accelerator Physics_ (Springer Press, 2015).


ifaamas [AAMAS ’22]Proc. of the 21st International Conference on Autonomous
Agents and Multiagent Systems (AAMAS 2022)May 9–13, 2022 Auckland, New
ZealandP. Faliszewski, V. Mascardi, C. Pelachaud, M.E. Taylor (eds.) 2022 2022
University of Dhaka University of Dhaka University of Massachusetts Amherst
University of Dhaka

# Learning Cooperation and Online Planning Through Simulation and Graph
Convolutional Network

Rafid Ameer Mahmud<EMAIL_ADDRESS>,  Fahim Faisal
<EMAIL_ADDRESS>,  Saaduddin Mahmud<EMAIL_ADDRESS>and  Md.
Mosaddek Khan<EMAIL_ADDRESS>
###### Abstract.

Multi-agent Markov Decision Process (MMDP) has been an effective way of
modelling sequential decision making algorithms for multi-agent cooperative
environments. A number of algorithms based on centralized and decentralized
planning have been developed in this domain. However, dynamically changing
environment, coupled with exponential size of the state and joint action
space, make it difficult for these algorithms to provide both efficiency and
scalability. Recently, Centralized planning algorithm FV-MCTS-MP and
decentralized planning algorithm Alternate maximization with Behavioural
Cloning (ABC) have achieved notable performance in solving MMDPs. However,
they are not capable of adapting to dynamically changing environments and
accounting for the lack of communication among agents, respectively. Against
this background, we introduce a simulation based online planning algorithm,
that we call SiCLOP, for multi-agent cooperative environments. Specifically,
SiCLOP tailors Monte Carlo Tree Search (MCTS) and uses Coordination Graph (CG)
and Graph Neural Network (GCN) to learn cooperation and provides real time
solution of a MMDP problem. It also improves scalability through an effective
pruning of action space. Additionally, unlike FV-MCTS-MP and ABC, SiCLOP
supports transfer learning, which enables learned agents to operate in
different environments. We also provide theoretical discussion about the
convergence property of our algorithm within the context of multi-agent
settings. Finally, our extensive empirical results show that SiCLOP
significantly outperforms the state-of-the-art online planning algorithms.

###### Key words and phrases:

Cooperation Learning, Multi-agent Markov Decision Process, Transfer Learning,
Graph Convolutional Network

##  1\. Introduction

Sequential decision making models for multi-agent environments hold the key to
many real life problems such as autonomous vehicles

Shalev-Shwartz et al. (2016), controlling robots Kober et al. (2013), resource
allocation Yang et al. (2018), games with multiple types of agents such as
Starcraft II Rashid et al. (2018) and so on. Cooperation among agents is an
essential part of these problems. The agents need to learn collaboration in
order to work together towards a common goal. This adds to the difficulty of
making decisions because the agents must interact not only with the
environment, but also with one another. As a result, they need to adapt to the
policies of other agents to learn cooperation. Another fundamental challenge
of this domain is the curse of dimensionality. With the increase in the number
of agents, the joint action space and number of states grows exponentially,
rendering single agent planning algorithms inefficient in these cases. In
multi-agent systems, the sequential decision making problem can be modeled as
a variant of Markov Decision Process (MDP) Sutton and Barto (2011), called
Multi-agent Markov Decision Process (MMDP) Boutilier (1996).

Over the years, a number of algorithms have been proposed to solve MMDPs. They
are broadly divided into two categories. Centralized planning and
decentralized planning. To solve MMDPs through centralized planning, the naive
approach is to extend the idea of single agent algorithms by treating all the
agents as a single agent and combining their joint action space to represent
the action space. However, due to the curse of dimensionality, this naive
approach fails to scale up to large environments. Therefore, a different
approach, decentralized planning, has been widely used to solve MMDPs. This is
accomplished by decomposing the multi-agent characteristics through
decentralizing their value function. Best et al. Best et al. (2019) proposed a
decentralized online planning algorithm to solve the dimensionality problem by
using parallel MCTS trees and periodic communication. On the other hand,
Aleksander et. al. Czechowski and Oliehoek (2020) introduced an algorithm ABC,
where agents do not use communication, rather they train their individual
policy prediction functions one at a time and try to induce cooperation by
using the learned agent models. However, ABC experiences poor performance
since, because of the lack of communication, it is unable to account for
penalties generated by agent interaction (see Section 4).

Guestrin et al. Guestrin et al. (2001) showed that communication among agents
can be represented as a coordination graph (CG) and joint value functions can
be estimated from the higher order value factorization. This encouraged other
approaches that use coordination graph to learn value functions or for policy
generation. Choudhury et al. Choudhury et al. (2021) proposed simulation based
anytime online planning algorithm FV-MCTS-MP, where factored value function is
used with CG to incorporate communication among agents. The proposed solution
solves both the cooperation and dimensionality problem but is not applicable
in dynamically changing environments. In such cases, the state space grows
exponentially as the size of the environment rises, limiting FV-MCTS-MP’s
ability to calculate state values and deliver good solutions given limited
computational resources.

Against this background, we introduce a MCTS based multi-agent anytime
planning algorithm, SiCLOP, that can scale for large environments and provide
realtime action selection capabilities. SiCLOP uses Graph Convolutional
Network (GCN) Kipf and Welling (2016) for individual agent policy prediction
and propagates the updated policies to other agents through Coordination
Graphs (CG). The agents use only their local information for policy prediction
and training. This enables them to learn cooperation in a small environments
and transfer that knowledge to a larger setting, making the training easier
and faster. We also introduce an effective pruning method to reduce the size
of the joint action space by sampling only a small subset of all possible
joint actions using Gibbs sampling Geman and Geman (1984) and iteratively
updating the actions of every individual agent. Finally, we evaluate our
algorithm empirically, and the result depicts that our algorithm outperforms
the state-of-the-art online planning algorithms in terms of solution quality,
adaptability, robustness and scalability.

##  2\. Problem Formulation and Background

In this section, we first formulate the MMDP problem. We then review the
general ideas of solving MMDPs with $onlineplanning$. Afterwards, we discuss
the impact of Coordination Graph (CG) and graph based neural networks in
improving the solution quality and adaptability,

### Multi-agent Markov Decision Process, MMDP

An MMDP can be formalized as a tuple $(\mathcal{N},S,A,\tau,\mathcal{U})$.
$\mathcal{N}$ is the set of all agents, indexed by $1,....,n$. We use $-i$ to
refer to all other agents except $i$. $S$ is the set of all possible states
and $A=A_{1}\times...\times A_{n}$, is a set of action spaces where $A_{i}$ is
the set of actions for agent $i\in\mathcal{N}$. A joint action $a\in A$ is the
combination of individual actions of all agents, $a=\langle
a_{1},...,a_{n}\rangle,a_{i}\in A_{i}$. $\tau:S\times A\rightarrow S$ is a
transition function that takes a joint action and changes the state according
to the deployed environment specifications. $\mathcal{U}=(u_{i},...,u_{n})$ is
the global utility function and for every agent
$i\in\mathcal{N},u_{i}:S\rightarrow\mathbb{R}$ denotes the utility function to
determine the payoff of states in $S$ for that agent. A policy is the task of
deciding which action to choose in any state. The policy for the $i$-th agent
is denoted by $\pi_{i}:S\rightarrow|A_{i}|$, a distribution over the action
space depending on the state the agent is in. The collection of policies of
all the agents is called the joint policy
$\pi=\langle\pi_{1},...,\pi_{n}\rangle$. If the joint policy is fixed, then
all the agents take actions according to that, and the expected reward or the
state value $V_{i}^{\pi}(s)$ of an state $s$ for agent $i$ is:

$V_{i}^{\pi}(s)=u_{i}(s)+\sum_{a\in A}\pi(a|s)V_{i}^{\pi}(\tau(s,a))$ (1)

And the expected reward of an action or action value $Q_{i}^{\pi}(s,a)$ of
state $s$ for agent $i$ in joint action $a$ is:

$Q_{i}^{\pi}(s,a)=u_{i}(s)+V_{i}^{\pi}(\tau(s,a))$ (2)

#### Best Response Policy

This $state\ value$ assumes that the policies of other agents are stationary,
although in multi-agent policy learning contexts, this is rarely the case.
Even if agent $i$ does not change its policy, the expected reward will change
when other agents change their policies.

In a multi-agent setting, the agents have to search the same state space to
create strategies that maximise their payoff. Here, the goal of the agents is
to adapt to the non-stationary policies of other agents. A $best\ response$
policy for agent $i$ is defined by the set $BR(\pi_{-i})$. It is the set of
policies that maximise the reward assuming that the joint policy of other
agents $\pi_{-i}$ is stationary.A policy $\pi^{*}$ belongs to $BR(\pi_{-i})$
if and only if the following inequality holds:

$\forall\pi_{i};\forall s\in S;\pi_{i}^{*}\in BR_{i}(\pi_{-i})\quad
V_{i}^{(\pi_{i}^{*},\pi_{-i})}(s)\geq V_{i}^{(\pi_{i},\pi_{-i})}(s)$ (3)

### Online Planning

Simulation based methods provide a means of policy learning and planning in
non-stationary environments where it is possible for the agents to find their
best response policy according to their current model of other agents. Monte
Carlo Tree Search (MCTS) Kocsis and Szepesvári (2006) is an effective tree
searching algorithm that uses simulation data to balance between exploration
of new solutions and exploitation of solutions with potential. There are four
major steps in the traditional MCTS simulation process. At any given moment,
with an already expanded simulation tree, the first step is to start from the
root and select a new leaf by going down the tree. While traversing, at every
node, scores for its child nodes are calculated and the child with the best
score is selected to go further down until a leaf is reached. The score is
calculated with,

$score_{a}=Q(s,a)+c*\sqrt{(\log N_{i})/n_{i}}$ (4)

where $Q(s,a)$ is the action value, $n_{i}$ is the visit count of action $a$
and $N_{i}$ is the total visit count of state $s$. The constant $c$ determines
the tradeoff between exploration and exploitation. The next step is to expand
the selected leaf and add its child nodes to the tree, called the expansion
step. After selection and expansion, the state value of the child nodes are
estimated, traditionally with random rollouts. The last step is to carry the
newfound information about the expanded node to its ancestors through
backpropagation. In an simulated tree, every node in that tree holds
information about its expected reward, the number of times that node was
visited and the edges to its child nodes.

AlphaZero Silver et al. (2017) takes a different approach with MCTS. In recent
studies, it has been shown that generating policies for large state and action
spaces can be done with deep neural networks rather than random plays Mnih et
al. (2015). Hence, instead of estimating state values with computationally
exhausting rollouts, a neural network is used in AlphaZero. The neural network
also provides an initial policy for the action space to better guide the tree
searching process. This not only makes the process faster, but also helps
focus on better solutions.

### Coordination Graph and Graph Neural Networks

Figure 1. Coordination Graph. The agents have fixed ranges and they are
connected if their communication ranges intersect. The edges determine which
agents share information.

In multi-agent environments, every agent interacts with a subset of other
agents. These interactions can be modeled using a graph structure called
Coordination Graph (CG), where the nodes are the agents and interacting agents
have edges between them. This structure captures the locality of the
dependencies of the agents regarding payoffs. This locality of dependency also
determines the cooperative policies of the agents. In Figure 1, the agents
create a forest of connected graphs based on their range. The edge creation
criteria can be predefined and the generated CG will denote which agents need
to communicate.

New models and architectures have been brought to light that use graphs in
training models to complete tasks regarding graphical structure, called Graph
Neural Networks (GNN). The use of GNNs has been shown to be effective in agent
modeling Böhmer et al. (2020)Liu et al. (2020). Graph Convolutional Networks
(GCN) is a popular GNN for graph embedding. GCN uses a message passing
mechanism that lets the nodes of a graph to aggregate and retain information
about their neighbors of different distances. This greatly helps in
representation learning of graphical structures and use it to generate
practical solutions.

##  3\. Multi-Agent Online Planning

Input: Set of initial states $S$, set of agents $\mathcal{N}$, limited steps
$t$, data storage $d$, transition function $\tau$, utilility dunction
$\mathcal{U}$

1 forall  _$state\in S$_ do

2 for  _$step$ in $0:t$_ do

3 $a\leftarrow$ SiCLOP-S ($state,\mathcal{N}$)

4 $state\leftarrow\tau(state,a)$

5 if  _$state$ is terminal_ then

6 continue

7

8 $reward_{i}\leftarrow\mathcal{U}(state,a)$

9 for  _$step$ in $0:t$_ do

10 store $(state_{step},reward)$ in data storage $d$

11 update SiCLOP-PNN neural network with Algorithm 2

12function $\texttt{SiCLOP-S}{(}s,\mathcal{N}{)}$

13 initialize $root$ with state $s$

14 while  _time limit not exceeded_ do

15 $\Delta\leftarrow$ select leaf using Eq. 4

16 if  _$\Delta$ is terminal_ then

17 break

18 $\Lambda\leftarrow\texttt{SiCLOP-P}(\Delta_{s},\mathcal{N},k)$

19 for  _$a$ in $\Lambda$_ do

20 create child node of $\Delta$ with state $\tau(\Delta_{s},a)$

21 $\rho\leftarrow$ state value predicted by SiCLOP-PNN

22 $reward\leftarrow\mathcal{U}(\Delta_{s},a)+\rho$

23 propagate $reward$ to ancestors up to $root$

24

25 $a\leftarrow$ most visited joint action of $\mathcal{N}^{\prime}$ in $root$

26 return $a$

27

28function $\texttt{SiCLOP-P}{(}s,\mathcal{N},k{)}$

29 initialize set of joint actions $\Lambda=\\{\\}$

30 $\pi\leftarrow$ SiCLOP-PNN($preprocess(s)$)

31 $a\leftarrow$ sample joint action from $\pi$

32 while  _$|\Lambda| <k$_ do

33 for  _$i$ in $1:|\mathcal{N}|$_ do

34 $s^{\prime}\leftarrow\tau(s,a_{-i})$

35 $\pi_{i}\leftarrow$ SiCLOP-NN($preprocess(s^{\prime})$)[i]

36 $a_{i}\leftarrow$ sample joint action from $\pi_{i}$

37

38 $\Lambda\leftarrow\Lambda\cup a$

39

40 return $\Lambda$

41

42function $\texttt{preprocess}{(}s{)}$

43 initialize $s^{*}$, empty set of observation states

44 for  _$i$ in $1:|\mathcal{N}|$_ do

45 $o_{i}\leftarrow$ observed elements of agent $i$ within fixed range
concatenated with environmental details

46 $s^{*}\leftarrow s^{*}\cup o_{i}$

47 $CG\leftarrow$ coordination graph based on fixed rules

48 return $s^{*},CG$

49

Algorithm 1 SiCLOP

Our proposed algorithm Simulation based Cooperation Learner and Online Planner
(SiCLOP) is an MCTS based MMDP solver. SiCLOP has three major components,
simulator SiCLOP-S, pruner SiCLOP-P, and predictor SiCLOP-NN. SiCLOP-S
modifies the traditional MCTS to accommodate the pruning capabilities of
SiCLOP-P and generate a deep simulation tree to find optimal policies.
SiCLOP-P is a novel pruning technique that iterates over the joint policy
space of the agents and uses best response policy predictor SiCLOP-NN to
sample a small number of potential joint actions from the exponential joint
action space. Lastly, SiCLOP-NN is a GCN based neural network that takes the
local information and interaction of the agents (represented by CG) as input
to predict best response policies for the agents and estimate the state value.
Combining these three components, SiCLOP becomes a scalable and adaptive MMDP
solver that is also capable of transfer learning.

The MMDP problems are represented as initial states. The algorithms takes a
set of initial states $S$, set of agents $\mathcal{N}$ and step limit $t$.
SiCLOP is run for every initial state $s$ in $S$ in lines 1-10. For $t$ steps,
the simulator SiCLOP-S finds an optimum joint action in line 3, and the state
transitions to the next step in line 4, and repeats unless terminated in line
5. After $t$ steps, the simulation for that initial state $s$ is finished and
the rewards are calculated in line 7. Then in lines 8-9, the state-action and
state-reward pairs are stored in a database $d$ and then the database is used
to train SiCLOP-NN in line 10. To summarize this section, the goal is to use
simulation (depends on the predictor) to find solutions and then use data from
that simulation to improve the predictor so that next time, SiCLOP can find
even better solutions.

###  3.1. SiCLOP-S, Simulation and Action Selection

The simulator starts with parameters initial state $s$ and set of all the
agents $\mathcal{N}$. The $root$ of the simulation tree is initialized by the
initial state $s$ in line 12.

There are four steps in the simulation process and they are repeated multiple
times within a fixed amount of time in the while loop in lines 13-22. In the
first step, a leaf is reached from the $root$ in line 14 by frequently
calculating the scores of the children of a node with Eq. 4 and choosing the
child with the highest score. Visit counts and state values of visited nodes,
while traversing down the tree, are updated. After choosing a leaf node, its
child nodes are created and the leaf is labeled as expanded. The number of
possible child nodes is exponential in magnitude with respect to the number of
agents. So, in the second step, a small subset of joint actions are sampled
from the joint action space with SiCLOP-P in line 17 and the child nodes are
created from the set of sampled actions to complete step three and four in
lines 18-22. After creating a child node in line 18, its state value is
determined in lines 20-21 by calculating penalties and prediction made by
SiCLOP-NN. Then the state value is backpropagated to its ancestors in line 22.

At the end of generating the simulation tree, the most visited joint action
which denotes the joint action with the most potential, is selected in line 23
and then returned.

###  3.2. SiCLOP-P, Pruning the Joint Action Space

The SiCLOP-P component of the algorithm described in lines 25-35, focuses on
pruning the joint action space for SiCLOP-S. We propose a novel joint action
sampling method inspired by Gibbs sampling to effectively search over the
joint policy space. First, we elaborate the sampling process. The predictor
SiCLOP-NN generates the best response policy for the agents given a state.
Using this predictor, the best response policy of every agent is sampled
considering the agent policies as static periodically. The sampling runs for
multiple cycle. At cycle $k$, the policy of agent $i$ is updated with the
conditional best response policy

$\pi^{k}_{i}=\Psi(\pi_{i}|\pi^{k}_{1},...,\pi^{k}_{i-1},\pi^{k-1}_{i+1},...,\pi^{k-1}_{n})$
(5)

In SiCLOP-P, the initial policy is represented by initial joint action
selection in line 28. Then, in lines 29-34, a fixed number of joint actions
are sampled. The agents are selected periodically in lines 30-32 where the
current policy of agent $i$ is replaced by the best response policy. To update
the policy, all other agents’ policies are considered to be static. The
current adopted policies of the agents are represented by $s^{\prime}$ in line
31. Then the best response policy is updated in line 32 and new action is
sampled from that policy in line 33. At the end of the cycles, in line 34, the
updated joint action is added to the sampled joint actions and after a fixed
number of cycles, the set of sampled joint actions are returned.

The policy generator SiCLOP-PNN predicts the best response policy for an
agent. So, the generated policy for agent $i$ belongs to $BR_{i}(\pi_{-i})$.
In Lemma 3.2, we show that SiCLOP-P eventually converges to an equilibrium
joint policy for a group of agents.

###### Lemma

In any iteration of sampling for agent $i$ in group $\mathcal{N}$, the new
sampled policy $\pi^{*}_{i}\in BR_{i}(\pi_{-i})$ will lead to non-decreasing
reward and iterating through all the agents of $\mathcal{N}^{\prime}$ will
lead to an equilibrium for that group.

###### Proof.

If the current joint policy of all the agents is $\pi$ and the updated policy
for agent $i$ is $\pi_{i}^{*}$, then for all $s\in S:$

$\displaystyle\forall i\in\mathcal{N}^{\prime}\quad
V_{i}^{(\pi_{i}^{*},\pi_{-i})}(s)\geq V_{i}^{(\pi_{i},\pi_{-i})}(s)$ (6)

If all the agents in the group acquire non-decreasing rewards, the combined
reward will also be non-decreasing. As both the joint action space and reward
are finite, iterating through the policy space will converge to a joint policy
where no individual policy can be updated to get a better result, hence
converging to an equilibrium for the group $\mathcal{N}^{\prime}$. ∎

This shows that with a well trained policy prediction function, the joint
policy of every group will converge to a best response policy equilibrium. But
in non-stationary environments, there might not a single joint policy
equilibrium for the agents to converge to. In that case, the agents
continuously create new best response policies to maximise their respective
payoffs.

The characteristics of the sampled joint actions change gradually throughout
the learning process. In the early stages of learning, the policies of the
agents are initialized to be random or according to prior knowledge. The set
sampled joint actions represent that initial joint policy. If it is random,
then the agents explore different strategies through simulation as a group and
develop their policies by finding strategies that lead to better payoffs. At
the same time, agents adapt to other agents’ policies and search for best
responses. And if it is initialized with handcrafted policy, then the agents
evolve by finding counter strategies and updating their existing policies. In
both cases, they initially search randomly for new strategies and select the
optimal ones in simulation. After repeating the learning process, the agents
start to generate more definitive policies, meaning some actions are more
preferred than others in their respective policies. This causes the sampling
method to sample joint actions that have higher probability of being chosen to
maximize rewards for all the agents. Directly predicting the optimal policy
can concentrate on single policies and might lose robustness. Rather,
simulating a search tree tests different ways of solving a MMDP with a set of
best response policies to choose from. This reduces the dependency on single
strategies and provides stable robust solutions.

Figure 2. Neural network architecture of SiCLOP-NN

Input: set of agents $\mathcal{N}$, prediction neural network SiCLOP-PNN and
data storage $d$

Output: Updated prediction neural network SiCLOP-NN

1 for  _$epoch$ in $1:nEpochs$_ do

2 for  _$i$ in $1:|\mathcal{N}|$_ do

3 prepare batch $B_{i}$ of state-action and state-reward pairs from $d_{i}$ of
recent episodes by random selection

4 $loss\leftarrow$ policy_error + value_error (Eq. 7 and 8)

5 minimize $loss$

6

Algorithm 2 Training algorithm for SiCLOP-NN

###  3.3. SiCLOP-NN, Prediction with Neural Network

The neural network SiCLOP-NN is a GCN based policy predictor and state value
estimator. It has two parts, feature extraction and prediction, as shown in
Figure 2. The first part is used for graph embedding as our algorithm focuses
on CGs. This feature extraction part consists of multiple GCN layers. The
number of GCN layers determine the number of hops of feature information among
agents. Two GCN layers would mean that the features of the agents will reach
other agents up to two hops away. The second part takes the extracted feature
vectors of the agents from the GCN block and generates predicted best response
policies and state value estimations for every agent. The state values are
aggregated to construct a combined state value. In Figure 2, The red outlined
pipeline is policy predictor and the blue outlined pipeline is state value
estimator. SiCLOP-NN returns a boltzman probability distribution over the
actions to enable weighted sampling of actions in SiCLOP-NN.

SiCLOP-NN takes the local information of the agents and a CG as input. The
setting of the local information can be defined, for example, the cells within
a limited range of the agent. The CG is also constructed according to
predefined rules. In SiCLOP, the input is created in the preprocess function
where the local information of the agents are collected into a set in lines
39-40. Then the set is returned with the constructed CG.

The training of the neural network occurs after a fixed number of simulation
episodes. In Algorithm 2, a fixed-sized batch $B_{i}$ of state-policy and
state-value pairs are randomly selected from the data from recent episodes for
every agent. The reason behind selecting data from recent episodes is to train
the neural network on recently found policies of the agents. This allows the
agents to adapt to new policies faster. Then predictions are made on that
batch of states to calculate the cross-entropy error between the real policy
and predicted policy for agent $i$. This is calculated using

$policy\\_error=\sum_{(s^{b},\pi^{b}_{i})\in B_{i}}-\pi^{b}\cdot ln(SiCLOP-
NN(s^{b},i))$ (7)

And the mean squared error between real value and predicted value is
calculated with

$value\\_error=\sum_{(s^{b},v^{b}_{i})\in B_{i}}||v^{b}_{i}-SiCLOP-
NN(s^{b},i)||^{2}$ (8)

The neural network is updated to minimize the combined loss.

The neural network trains and predicts on states consisting of any number of
agents, with their internal interactions passed as CG into the GCN block. This
adds the flexibility of performing predictions on variable number of agents.
On the other hand, the input is a collection of observations of the agents
dictated by a constraint of locality, i.e. communication range. SiCLOP-NN’s
architecture allows making prediction in environments of any size and shape.
If the agents are trained in a smaller environment with fewer agents, the
learned agents can still apply the learned policies on much larger
environments with more agents as their policy prediction only depends on the
locality of information.

##  4\. Empirical Results

Figure 3. In this environment, green cells are destinations, blue cells are
agents and red cells are obstacles. The arrows from blue cells to green denote
the destination of the agents.

To test and compare our algorithm, we used a modified version of the drone
delivery environment Choudhury et al. (2021). The modification is the
inclusion of obstacles to make the environment more realistic, making the
movement of the agents difficult. The environment is shown in Figure 3. Here,
agents need to reach their designated objectives but agents can collide with
the obstacles or other agents if they try to enter the same cell (can be
objective cell) at the same time. They can also choose actions that take them
out of the environment area. These actions are heavily penalized. The goal of
the environment design is to test if the agents can learn to avoid penalizing
actions by using cooperation and complete the tasks, as seen in most realistic
problems. Every agent gets $1.0$ point for reaching an objective, $-1.0$ for
collision and $-0.5$ for going out of bounds. There are also penalties if
multiple agents stick too close to each other. Agents are rewarded if they get
closer to their objectives. The agents have a total of $9$ possible actions,
going to the $8$ surrounding cells and staying in place. The MMDP episode ends
when all agents reach their goals or the pre-defined number of steps are used
up. The results are shown as average score per agent as the number of agents
vary in multiple environmental setting.

Next, we will present an analysis of how the algorithm performs, how it can
scale up and how it performs compared to FV-MCTS-MP and ABC planning and
learning algorithms.

###  4.1. Agent Training and Transfer Learning

The first experiment is to observe how the agents improve their policies to
get better results. In this experiment, there are $8$ agents in a $15x15$ grid
with $15$ obstacles. The agents have maximum 40 steps to complete the tasks. A
total of 500 randomly generated problem states were generated and SiCLOP was
tested on these 500 solution episodes. Agent policy predictor was trained
every 10 episodes. For every step in the episode, the simulator can either run
for a maximum amount of time or can expand limited number of simulation nodes.
The algorithm was tested on three such settings. To test the scaling up
ability of the algorithm, the trained policy predictors were used to test if
it can generate good quality solutions on new environments of different sizes
and difficulties.

Figure 4. Evolution of the agents in different settings Environment List
---
Grid Size |  Agents |  Obstacles |  Average Score
8x8 |  5 |  4 |  1.004
20x20 |  20 |  10 |  1.06
50x50 |  100 |  150 |  1.35
100x100 |  150 |  200 |  1.52
Table 1. Transferring and testing the learned agents in new environments

From the results shown in Figure 4, we can see that at the beginning of the
experiment, the agents adopt random policies and as the episodes progress they
start to adapt to the rules of the environment and converge to a high score
around $1.0$. An average score around $1.0$ denotes that the agents are
avoiding penalties and succeeding at their given tasks. If the simulation
limitations are tighter, the algorithm struggles to find good solutions at the
beginning, even yielding to negative scores. But as SiCLOP is allowed to do
more simulation before taking actions, it generates better solutions from the
start, training the agents on quality solutions. This helps converging to
optimal scores much faster, as we can see in Figure 4. It is also shown that,
in general, more simulation leads to better score on average.

After training the agents in a $15x15$ grid, the policy predictor is
transferred to randomly generated environments. These environments differ in
size, number of agents and obstacles. In Table 1, we can see that the average
score per agent in these environments are above $1.0$, which implies that the
algorithm is capable of adapting to a varying environment and avoid penalties
by using knowledge from a different environment. But in case of other planning
algorithms, FV-MCTS-MP stores state values of visited states so it cannot use
any prior knowledge on new environments and ABC uses the entire state to
predict policies, so it can not be applied to an environment of different
size. Transfer learning lets the algorithm to train on different environments,
enabling agents to initially train in a small environment and then train in
larger environments. As seen from the experiment, the prediction model can
transfer the knowledge allowing the training in large environments to converge
faster. This can save time and computational resources in case of complex
environments and large single agent action spaces, making the algorithm
applicable to most realistic MMDP problems.

Figure 5. Comparison of adaptability to collision (top) and out of bounds
(bottom)

###  4.2. Learning Cooperation

In the second experiment, the cooperation learning and adaptability of the
algorithms were tested. SiCLOP and ABC are both model learning algorithms and
FV-MCTS-MP is a centralized tree searching algorithm. While ABC adopts a
decentralized approach, SiCLOP chooses and learns its policies in a
centralized manner. To test the cooperation learning ability of the agents,
both algorithms were run on the same 500 episodes with $5$ seconds time limit
per action in a $15x15$ grid with $5$ obstacles and $8$ agents. In Figure 5,
we can see how the algorithms learn cooperation and adapt to handle
$collision$ and $outofbounds$ penalties.

From Figure 5, it is clear that ABC fails to learn how to avoid collisions but
can learn to avoid going out of bounds. The penalty of going out of bounds
depends on individual agents and collision mostly depends on the interaction
of the agents. ABC learns the agent models one at a time and tries every
available action. This helps the agent models to learn individual actions
effectively. But in case of collision, it is very hard to avoid without
communication among agents, causing the ABC agent models to not converge to
effective policies. As there is no communication, there is no differentiation
between successful and unsuccessful joint actions before the actions are
executed. FV-MCTS-MP consistently outperforms ABC in cooperation but due to
time limit and small number of iterations, there were no significant increase
in performance. SiCLOP outperforms both algorithms as it searches over the
joint policy space and needs lesser number of iterations to find cooperative
joint policies to avoid the penalties. Throughout the training process, SiCLOP
shows a stable performance and converges to almost no average penalty.

Contrary to the collision penalty, all three algorithms manage to adapt to the
out of bounds penalty. Before training, FV-MCTS-MP outperforms both SiCLOP and
ABC but after training these two perform better. Similar to previous
experiment, FV-MCTS-MP’s performance does not get better because of the
limitations. The reason behind the similar performances of SiCLOP and ABC is
that avoiding out of bounds penalty does not require any cooperation.

###  4.3. Performance Comparison

In this experiment, SiCLOP was compared to FV-MCTS-MP and ABC by average
scores accumulated on gradually increasing size of environments and number of
agents. All of the algorithms had $10$ seconds time limit per action. SiCLOP
and ABC algorithms were trained on 500 episodes. In Figure 6, we can see the
comparison in scores.

Figure 6. Performance comparison

The algorithms were tested in three different environments of varying sizes.
As the size of the environment gets larger, accumulating a high score becomes
harder as there are more obstacles and the objectives are farther. From this
experiment, it is clear that as the environment gets larger, ABC is more prone
to penalties and accumulates lower average scores. FV-MCTS-MP, with limited
time and number of iterations, is unable to find sequence of actions that lead
to high score. As the simulation tree of FV-MCTS-MP is wide, it calculates the
expected action value over a wide and shallow tree, causing it to have average
negative values. So, FV-MCTS-MP prefers not to move the agents to minimize
penalties. But, as SiCLOP searches over a deeper search tree, it can find
better solutions to the MMDP problems. Also, the size of the environment
hardly affects the performance of SiCLOP as the agents only use their local
information and interactions. These attributes and better cooperation learning
allow SiCLOP to outperform both ABC and FV-MCTS-MP by a significant margin.

##  5\. Conclusions and Future Work

In this paper, we introduced a scalable and adaptive algorithm to solve MMDP.
Our algorithm SiCLOP uses local information of the agents to construct state
dependant dynamic CGs and uses it to find optimum policies through simulation
and effective pruning. We have shown that our algorithm can adapt to larger,
more realistic environments and outperform existing online MMDP solvers.

SiCLOP is an algorithm that uses policy predictions as a tool to find better
policies, which implies that the quality of the initial prediction holds great
significance on the final outcome. More work is required for finding better
ways to train the prediction models, which may include newer ways of selecting
self-play data on which the model is trained on and better neural network
architectures that suit the process better.

## References

  * (1)
  * Bernstein et al. (2009) Daniel S Bernstein, Christopher Amato, Eric A Hansen, and Shlomo Zilberstein. 2009.  Policy iteration for decentralized control of Markov decision processes.  _Journal of Artificial Intelligence Research_ 34 (2009), 89–132.
  * Best et al. (2019) Graeme Best, Oliver M Cliff, Timothy Patten, Ramgopal R Mettu, and Robert Fitch. 2019.  Dec-MCTS: Decentralized planning for multi-robot active perception.  _The International Journal of Robotics Research_ 38, 2-3 (2019), 316–337.
  * Böhmer et al. (2020) Wendelin Böhmer, Vitaly Kurin, and Shimon Whiteson. 2020\.  Deep coordination graphs. In _International Conference on Machine Learning_. PMLR, 980–991.
  * Boutilier (1996) Craig Boutilier. 1996\.  Planning, learning and coordination in multiagent decision processes. In _TARK_ , Vol. 96. Citeseer, 195–210.
  * Buşoniu et al. (2010) Lucian Buşoniu, Robert Babuška, and Bart De Schutter. 2010.  Multi-agent reinforcement learning: An overview.  _Innovations in multi-agent systems and applications-1_ (2010), 183–221.
  * Choudhury et al. (2021) Shushman Choudhury, Jayesh K Gupta, Peter Morales, and Mykel J Kochenderfer. 2021.  Scalable Anytime Planning for Multi-Agent MDPs.  _arXiv preprint arXiv:2101.04788_ (2021).
  * Claes et al. (2017) Daniel Claes, Frans Oliehoek, Hendrik Baier, Karl Tuyls, et al. 2017\.  Decentralised online planning for multi-robot warehouse commissioning. In _AAMAS’17: PROCEEDINGS OF THE 16TH INTERNATIONAL CONFERENCE ON AUTONOMOUS AGENTS AND MULTIAGENT SYSTEMS_. 492–500.
  * Claes et al. (2015) Daniel Claes, Philipp Robbel, Frans Oliehoek, Karl Tuyls, Daniel Hennes, and Wiebe Van der Hoek. 2015\.  Effective approximations for multi-robot coordination in spatially distributed tasks. In _Proceedings of the 14th international conference on autonomous agents and multiagent systems (AAMAS 2015)_. International Foundation for Autonomous Agents and Multiagent Systems, 881–890.
  * Czechowski and Oliehoek (2020) Aleksander Czechowski and Frans Oliehoek. 2020.  Decentralized MCTS via Learned Teammate Models.  _arXiv preprint arXiv:2003.08727_ (2020).
  * Foerster et al. (2016) Jakob N Foerster, Yannis M Assael, Nando De Freitas, and Shimon Whiteson. 2016.  Learning to communicate with deep multi-agent reinforcement learning.  _arXiv preprint arXiv:1605.06676_ (2016).
  * Geman and Geman (1984) Stuart Geman and Donald Geman. 1984.  Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images.  _IEEE Transactions on pattern analysis and machine intelligence_ 6 (1984), 721–741.
  * Guestrin et al. (2001) Carlos Guestrin, Daphne Koller, and Ronald Parr. 2001\.  Multiagent Planning with Factored MDPs.. In _NIPS_ , Vol. 1. 1523–1530.
  * Guo et al. (2014) Xiaoxiao Guo, Satinder Singh, Honglak Lee, Richard L Lewis, and Xiaoshi Wang. 2014.  Deep learning for real-time Atari game play using offline Monte-Carlo tree search planning. In _Advances in neural information processing systems_. 3338–3346.
  * Gupta et al. (2017) Jayesh K Gupta, Maxim Egorov, and Mykel Kochenderfer. 2017\.  Cooperative multi-agent control using deep reinforcement learning. In _International Conference on Autonomous Agents and Multiagent Systems_. Springer, 66–83.
  * Kipf and Welling (2016) Thomas N Kipf and Max Welling. 2016.  Semi-supervised classification with graph convolutional networks.  _arXiv preprint arXiv:1609.02907_ (2016).
  * Kober et al. (2013) Jens Kober, J Andrew Bagnell, and Jan Peters. 2013\.  Reinforcement learning in robotics: A survey.  _The International Journal of Robotics Research_ 32, 11 (2013), 1238–1274.
  * Kocsis and Szepesvári (2006) Levente Kocsis and Csaba Szepesvári. 2006.  Bandit based monte-carlo planning. In _European conference on machine learning_. Springer, 282–293.
  * Liu et al. (2020) Yong Liu, Weixun Wang, Yujing Hu, Jianye Hao, Xingguo Chen, and Yang Gao. 2020\.  Multi-agent game abstraction via graph attention neural network. In _Proceedings of the AAAI Conference on Artificial Intelligence_ , Vol. 34. 7211–7218.
  * Matignon et al. (2007) Laëtitia Matignon, Guillaume J Laurent, and Nadine Le Fort-Piat. 2007.  Hysteretic q-learning: an algorithm for decentralized reinforcement learning in cooperative multi-agent teams. In _2007 IEEE/RSJ International Conference on Intelligent Robots and Systems_. IEEE, 64–69.
  * Mnih et al. (2015) Volodymyr Mnih, Koray Kavukcuoglu, David Silver, Andrei A Rusu, Joel Veness, Marc G Bellemare, Alex Graves, Martin Riedmiller, Andreas K Fidjeland, Georg Ostrovski, et al. 2015\.  Human-level control through deep reinforcement learning.  _nature_ 518, 7540 (2015), 529–533.
  * OroojlooyJadid and Hajinezhad (2019) Afshin OroojlooyJadid and Davood Hajinezhad. 2019.  A review of cooperative multi-agent deep reinforcement learning.  _arXiv preprint arXiv:1908.03963_ (2019).
  * Panait and Luke (2005) Liviu Panait and Sean Luke. 2005.  Cooperative multi-agent learning: The state of the art.  _Autonomous agents and multi-agent systems_ 11, 3 (2005), 387–434.
  * Rashid et al. (2018) Tabish Rashid, Mikayel Samvelyan, Christian Schroeder, Gregory Farquhar, Jakob Foerster, and Shimon Whiteson. 2018.  Qmix: Monotonic value function factorisation for deep multi-agent reinforcement learning. In _International Conference on Machine Learning_. PMLR, 4295–4304.
  * Shalev-Shwartz et al. (2016) Shai Shalev-Shwartz, Shaked Shammah, and Amnon Shashua. 2016\.  Safe, Multi-Agent, Reinforcement Learning for Autonomous Driving.  _arXiv preprint arXiv:1610.03295_ (2016).
  * Silver et al. (2017) David Silver, Thomas Hubert, Julian Schrittwieser, Ioannis Antonoglou, Matthew Lai, Arthur Guez, Marc Lanctot, Laurent Sifre, Dharshan Kumaran, Thore Graepel, et al. 2017\.  Mastering chess and shogi by self-play with a general reinforcement learning algorithm.  _arXiv preprint arXiv:1712.01815_ (2017).
  * Silver and Veness (2010) David Silver and Joel Veness. 2010.  Monte-Carlo planning in large POMDPs. Neural Information Processing Systems.
  * Sunehag et al. (2017) Peter Sunehag, Guy Lever, Audrunas Gruslys, Wojciech Marian Czarnecki, Vinicius Zambaldi, Max Jaderberg, Marc Lanctot, Nicolas Sonnerat, Joel Z Leibo, Karl Tuyls, et al. 2017\.  Value-decomposition networks for cooperative multi-agent learning.  _arXiv preprint arXiv:1706.05296_ (2017).
  * Sutton and Barto (2011) Richard S Sutton and Andrew G Barto. 2011.  Reinforcement learning: An introduction.
  * Tampuu et al. (2017) Ardi Tampuu, Tambet Matiisen, Dorian Kodelja, Ilya Kuzovkin, Kristjan Korjus, Juhan Aru, Jaan Aru, and Raul Vicente. 2017\.  Multiagent cooperation and competition with deep reinforcement learning.  _PloS one_ 12, 4 (2017), e0172395.
  * Watkins and Dayan (1992) Christopher JCH Watkins and Peter Dayan. 1992.  Q-learning.  _Machine learning_ 8, 3-4 (1992), 279–292.
  * Yang et al. (2018) Yaodong Yang, Rui Luo, Minne Li, Ming Zhou, Weinan Zhang, and Jun Wang. 2018\.  Mean field multi-agent reinforcement learning. In _International Conference on Machine Learning_. PMLR, 5571–5580.
  * Zhang and Lesser (2010) Chongjie Zhang and Victor Lesser. 2010.  Multi-agent learning with policy prediction. In _Twenty-fourth AAAI conference on artificial intelligence_.