# Photoacoustic Tomography using a Michelson Interferometer with Quadrature
Phase Detection

Rory W. Speirs  School of Physics, The University of Melbourne, Victoria,
3010, Australia  Alexis I. Bishop<EMAIL_ADDRESS>School of Physics,
Monash University, Victoria, 3800, Australia

###### Abstract

We present a pressure sensor based on a Michelson interferometer, for use in
photoacoustic tomography. Quadrature phase detection is employed allowing
measurement at any point on the mirror surface without having to retune the
interferometer, as is typically required by Fabry-Perot type detectors. This
opens the door to rapid full surface detection, which is necessary for
clinical applications. Theory relating acoustic pressure to detected acoustic
particle displacements is used to calculate the detector sensitivity, which is
validated with measurement. Proof-of-concept tomographic images of blood
vessel phantoms have been taken with sub-millimeter resolution at depths of
several millimeters.

Photoacoustic imaging has the potential to become a routinely used medical
imaging modality, combining the superior contrast of optical techniques, with
the penetration depth of ultrasound.Minghua and Lihong V. (2006); Li and Wang
(2009); Wang and Hu (2012) Its inherent ability to distinguish regions of
contrasting optical absorption make it an ideal candidate for imaging vascular
structure, with possible applications in diagnosis of strokeJoshi and Agarwal
(2010) and early stage cancers.Esenaliev  _et al._ (1999); Oraevsky  _et al._
(2001); Yang  _et al._ (2012)

Photoacoustic tomography (PAT) is potentially capable of producing real time,
three-dimensional (3D), high resolution images to depths of several
centimeters.Wang (2008); Xiang  _et al._ (2013) The requirements of the
ultrasonic detection system to achieve this are formidable, and as yet, no
system has emerged which satisfy all criteria simultaneously.

Though much emphasis has been placed on detector sensitivity, there are other,
equally pressing requirements of a high resolution system. These include the
need for a large detection surface, with high spatial resolution. Of
particular importance for a system that can be used in clinical context is the
ability to rapidly capture data, ideally over the whole detector surface
simultaneously. A detector which is transparent to the excitation light is
also favourable, as this allows a large amount of optical energy to be dumped
uniformly on the region being imaged.

Piezoelectric detectors struggle with many of these requirements,Hou  _et al._
(2007) and so a wide variety of optical detectors have been developed.Beard
and Mills (1996); Hamilton and O’Donnell (1998); Paltauf  _et al._ (2007);
Chow  _et al._ (2011) Planar Fabry-Perot based systems show good sensitivity
and bandwidth response, but are typically slow to acquire data because of the
need to tune the probing laser at each point on the detector surface to
achieve peak sensitivity.Zhang  _et al._ (2008) Simultaneous two-dimensional
(2D) data collection has been demonstrated with these systems,Lamont and Beard
(2006) but large detection areas are challenging to produce because of the
difficulty in creating polymer coatings of uniform thickness.

Microring resonators have been made with impressive sensitivity, element size
and can be made transparent. However coupling and addressing a large array of
microrings will be difficult, so it is yet to be seen if simultaneous full
surface detection is achievable.Huang  _et al._ (2008)

Pressure dependent optical reflectance detectors have been demonstrated with
the ability to capture pressure data over a whole surface
simultaneously,Paltauf  _et al._ (1999); Paltauf and Schmidt-Kloiber (1997)
without the need for complicated nanofabrication techniques of some other
methods. Moreover, use of fast-framing, or gated charge-coupled devices (CCDs)
simplifies data collection, and allows for high spatial resolution over a
large detection surface. However, the detection sensitivity of this type of
system has so far been only modest, and may be insufficient for high
resolution imaging of biological tissue.

We have developed a detector based on a Michelson interferometer (MI) with
quadrature phase detection. This detector has comparable sensitivity to other
optical detectors in the literature, but has the potential to perform high
resolution measurements over a full 2D surface simultaneously, without the
need for any position dependent sensitivity tuning.

The MI acts as an ultrasound sensor simply by acoustically coupling an
ultrasound source to a mirror in one of the arms. The acoustic wave of
pressure $p$, has an associated particle displacement $\xi$, which shifts the
position of the mirror as the wave passes through it. This change in position
adjusts the relative phase of the laser beams in the two arms, resulting in a
change in fringe brightness at the output of the interferometer. For a small
amplitude wave travelling in the $x$ direction at time $t$, pressure and
displacement are related by:

$p(x,t)=-E\frac{\partial\xi(x,t)}{\partial x},$ (1)

where $E$ is the appropriate modulus of elasticity for the medium.Blitz (1963)

The intensity, $I$ of the recombined beam in a standard MI varies sinusoidally
with mirror position:

$I=\frac{I_{0}}{2}(1+\cos(\phi)),$ (2)

where $\phi=4\pi n\xi/\lambda$. Here, $I_{0}$ is the input intensity, $\phi$
is the phase, $\lambda$ is the wavelength of the probe beam, and $n$ is the
refractive index of the arm where the mirror position is changing.

Michelson interferometers previously used in ultrasonic detection therefore
suffer from the same problem as Fabry-Perot type detectors, in that the laser
wavelength (or mirror position) must be tuned to a sensitive region at each
point in order to obtain good optical modulation for a given mirror
displacement.Mezrich  _et al._ (1976)

The need for tuning was removed from our system by employing quadrature phase
detection. In quadrature phase detection, two orthogonal linear polarizations
are used to simultaneously obtain two separate interference patterns at the
output of the interferometer, which have a relative phase difference of
$\pi/2$. This phase difference ensures that the interference pattern of at
least one of the polarization components is always in a sensitive region.

The phase shift between the two polarization components is created by first
linearly polarizing the light at 45∘ from the vertical or horizontal axis of
the polarizing beamsplitter. A liquid crystal variable waveplate is placed in
one of the arms, which retards the phase of one polarization component
relative to the other by nominally $\pi/4$ in both the forward and reverse
trips. A variable retarder is used instead of a fixed $\lambda/8$ wave plate
to compensate for small amounts of birefringence present in other optical
components.

The _phase sensitivity_ , $\frac{dI}{d\phi}$ (which is the optical intensity
modulation per radian of phase) of a normal MI varies between $0~{}\rm
rad^{-1}$ and $I_{0}/2~{}\rm rad^{-1}$. For an MI with quadrature detection,
the phase sensitivity for each polarization is simply added together, so is
always between $I_{0}/2~{}\rm rad^{-1}$ and $I_{0}/\sqrt{2}~{}\rm rad^{-1}$.
This ensures that the total sensitivity of the system is always at least as
high as the maximum of a standard MI, irrespective of absolute mirror
position.

A diagram of the detector setup can be seen in Fig. 1. The mirror position is
recovered from the detected intensities of the two polarizations, $I_{1,2}$ by
first scaling them between $-1$ and $1$, then treating them as points on the
unit circle: $\phi=\mathrm{atan2}(I_{1},I_{2})$.

Figure 1: A schematic of the Michelson interferometer detector. Ultrasound
passing through the sensing mirror changes the path length of the light, which
alters the phase of the interference pattern. Quadrature phase detection
allows high sensitivity to be achieved, whatever the mirror positions.

Our setup uses an expanded $5~{}\rm mW$, $633~{}\rm nm$ continuous wave
helium-neon laser as the probe, and the detectors are amplified photodiodes
with a $20~{}\rm MHz$ bandwidth. The photodiodes are apertured, which sets the
spatial resolution of the system. The signal of each photodiode is recorded
using a digital oscilloscope. The sensing mirror is a $150~{}\rm\mu m$ thick
glass substrate with a gold reflective coating. To maintain mechanical
stability, the mirror is mounted on an optical window by bonding it around the
perimeter with epoxy resin. This method of bonding also ensures the presence
of an air gap between the reflective surface and the window, eliminating the
possibility of acoustic waves propagating into the window. The glass surface
of the mirror is acoustically coupled to the medium being imaged, with the
window becoming part of the interferometer arm. The resulting mirror can be
seen in Fig. 2.

Figure 2: (a) Schematic of the sensing mirror showing probe laser. (b) The
back side of the sensing mirror used in this experiment.

It is necessary to use such a thin mirror for two reasons. Firstly, a thick
mirror suffers from the acoustic wave reflecting back and forth off the
boundary of the substrate. This reflected wave interferes with the incoming
wave, making the detected displacement useless. Multiple reflections are
strongly suppressed in a mirror that is much thinner than the acoustic
wavelength, so the displacement on the mirror surface accurately represents
the incoming wave. Secondly, a mirror of thickness significantly greater than
the acoustic wavelength is able to support surface waves. These surface, or
_Rayleigh_ waves can exist whenever there is an impedance mismatch between two
media, and are generated on the surface of the mirror when the photoacoustic
pulse first reaches the boundary. The Rayleigh waves then propagate outwards
across the surface and interfere with incoming photoacoustic waves.

To see the fundamental limitations of the detection system, it is useful to
look at the theoretically achievable sensitivity. Like all Fabry-Perot and
piezoelectric type detectors, the MI detector is sensitive to particle
displacements, rather than directly to pressure. However because the MI
detects the absolute position of a single plane, rather than the relative
position of two planes, its sensitivity to pressure is easier to describe
analytically. For a sinusoidal acoustic wave of pressure amplitude $p_{0}$,
travelling in the positive $x$ direction, the corresponding particle
displacement is given by:

$\xi(x,t)=\frac{p_{0}}{2\pi\nu z}\cos(kx-2\pi\nu t),$ (3)

where $k$ is wavenumber, $\nu$ is frequency, and $z$ is the specific acoustic
impedance of the propagation medium. For optical based pressure sensors the
sensitivity, $S$ can be simply given as a proportion of optical intensity
modulation per unit acoustic pressure: $S=1/I_{0}\frac{dI}{dp}$. Using the
chain rule, this may be expanded to
$\frac{dI}{dp}=\frac{dI}{d\phi}\frac{d\phi}{d\xi}\frac{d\xi}{dp}$. For the
quadrature MI, $\frac{dI}{d\phi}$ is always at least $0.5I_{0}$, and
$\frac{d\phi}{d\xi}$ is simply calculated from the expression for $\phi$.
$\frac{d\xi}{dp}$ can be calculated from Eqn. 3, however it must be modified
to describe the setup employed in our system. Firstly, since the pressure wave
must propagate from the original medium into the glass mirror substrate, the
pressure must be multiplied by the transmission coefficient:
$T_{p}=2z_{2}/(z_{1}+z_{2})$, where $z_{1,2}$ are the specific acoustic
impedances of the first and second media respectively.Shutilov (1988) Also,
the mirror is essentially on a free boundary (since $z_{air}\ll z_{glass}$),
so the particle displacement will be twice as great as in the bulk. Combining
these terms gives the expression for the frequency dependent sensitivity:

$S(\nu)=\frac{4n}{\nu\lambda(z_{1}+z_{2})}.$ (4)

Typical values of acoustic impedance for water and glass are $1.5\times
10^{6}~{}\rm Pa~{}s~{}m^{-1}$ and $13.1\times 10^{6}~{}\rm Pa~{}s~{}m^{-1}$
respectively. Taking $n=1$ as the refractive index of air, and letting
$\lambda=633~{}\rm nm$, the sensitivity of our detector is
$S(\nu)=0.43/\nu~{}\rm Hz~{}Pa^{-1}$. For a photoacoustic wave of frequency
$1~{}\rm MHz$, the sensitivity of the detector should be 4.3% optical
modulation per $100~{}\rm kPa$ of peak acoustic pressure. This value can be
compared directly with the sensitivity for an optical reflectance based
detector of 0.19 to 0.81% reported elsewhere in the literature.Paltauf  _et
al._ (1999)

Fig. 3 shows the actual detected displacement caused by a photoacoustic wave
from a single source, positioned $3.3~{}\rm mm$ directly behind the detector.
To recover the acoustic pressure from the displacement, the temporal
derivative must be taken according to Eqn. 1, where $c=dx/dt$ has been used to
change the variable of differentiation. A post processing low pass filter is
applied to the detected displacement before the derivative is taken. This
ensures the calculated pressure does not contain unphysical spikes which are
artefacts of taking the derivative of a noisy signal.

Figure 3: (a) The detected mirror displacement due to a photoacoustic wave
produced from a single source and (b), the corresponding pressure wave. A
comparison between simulated and experimental pressure data shows the
calculated sensitivity for the system is accurate.

The noise equivalent displacement in this trace is approximately $5~{}\rm\AA$,
where the signal has been averaged for 64 pulses. The majority of the noise in
this case is due to the amplified photodiodes, which may be improved using
different equipment.

Fig. 3 also shows how the experimentally detected photoacoustic wave compares
with a simulated one, and the agreement between the two support the values of
sensitivity calculated previously. In the experiment, a long straight silicone
tube with internal and external diameters of $0.5~{}\rm mm$ and $1.3~{}\rm mm$
respectively, was filled with diluted India ink with absorption coefficient
$35~{}\rm cm^{-1}$ at $1064~{}\rm nm$. This value is similar to the optical
absorption coefficient of blood at wavelengths commonly used in PAT. The tube
was illuminated with a $10~{}\rm ns$ pulse of $1064~{}\rm nm$ light with
fluence of $25~{}\rm mJ~{}cm^{-2}$. This fluence is well below the maximum
permissible exposure for human skin of $100~{}\rm mJ~{}cm^{-2}$ based on
American National Standards Institute recommendations.ANS (2007) The pulsed
beam was collimated, and had a 1/e diameter of $1~{}\rm cm$. The simulation
was performed using the method described by Köstli _et al._ Kostli  _et al._
(2001), where Fig. 3b shows a one-dimensional (1D) slice of the full 3D
simulation. The parameters used in the simulation were the same as described
for the experiment, however the simulation assumed an acoustically homogeneous
propagation medium, so any acoustic effects of the silicone tube were ignored.
The large diameter excitation beam minimized intensity variation across the
tube, though due to the Gaussian nature of the beam, the stated fluence at the
sample is only accurate to within 15%. The fluence used in the simulation was
adjusted within this range until the simulated amplitude most closely matched
the detected one. As such, the calculated sensitivity is also only accurate to
within this margin.

For the experimental measurement, the tube was submerged in a water filled
glass cell. One wall of the cell was made of acoustically transparent
polyethylene film, which was coupled to the sensing mirror of the
interferometer with commercial ultrasound coupling gel. The detected pressure
represents the pressure in the bulk of the glass substrate. This value was
divided by the pressure transmission coefficient to give the acoustic pressure
in the water, allowing a direct comparison with simulation.

To create photoacoustic images, the photodiodes could be moved laterally to
build up 1D or 2D scans. However the lower laser intensity at the edges of the
expanded probe beam meant that data collected in these regions had a lower
signal to noise ratio (SNR). Instead, the sample itself was scanned laterally.
This allowed different regions of the generated ultrasonic wavefield to be
sampled in an equivalent manner to scanning the photodiodes, while maintaining
a high SNR. A 1D scan of the wavefield produced by the cylindrical optical
absorber is shown in Fig. 4. The source was the same India ink filled silicone
tube as previously described, illuminated by the same pulsed laser. The
photodiode aperture was set to a diameter of $200~{}\rm\mu m$, which was the
same as the lateral step size.

Figure 4: (a) Detected acoustic pressure from a single photoacoustic source
and (b), the corresponding reconstructed image of the source. The size and
location of the source is indicated by the circle in the magnified inset.

The reconstruction shown in Fig. 4b was performed using the _kspaceLineRecon_
function of the _k-Wave_ photoacoustic package.Treeby and Cox (2010) The
algorithm is based on Fourier transforms, and it is theoretically exact if
pressure is detected over an infinite plane for infinite time.Kostli  _et al._
(2001) In the reconstruction, the source has been positioned correctly, but
has suffered some blurring and distortion which is consistent with other
implementations of this inversion technique.Paltauf  _et al._ (2007) The
blurring is unsurprising given that the diameter of the source was only
$500~{}\rm\mu m$ and the spatial separation of each data point was
$200~{}\rm\mu m$.

The resolution of our system is currently limited by the size of the aperture
in front of the photodiodes, which in turn is limited by the need to get
sufficient laser power to the photodiodes. This is easily improved by
increasing the laser power, and could be achieved using inexpensive diode
lasers. Using a shorter wavelength probe laser would also be a simple way to
boost sensitivity according to Eqn. 4.

Manually scanning the photodetectors (or the source) to build up an image is
too slow for real-time imaging applications, so any useful system must
ultimately be capable of performing simultaneous detection over the whole
surface. This could be simply achieved in our system by replacing the
photodiode and aperture arrangement with fast gated intensified charge-coupled
device cameras. CCDs detect intensity at many points across their surface
simultaneously, eliminating the need to move the sample or the detector. The
potential to use CCDs is a significant advantage over other proposed optical
photoacoustic detectors, which have no easy route to simultaneous measurement
of all elements over a large surface.

The current configuration of our sensing mirror (constructed from a thin glass
substrate) may cause limitations to the detectable bandwidth needed in higher
resolution systems, due to the possible reappearance of Rayleigh waves at
higher frequencies. However this could be addressed by using polymer
substrates impedance matched to water. Also, where deeper imaging is required,
the bandwidth requirements of the detector are much more forgiving, since very
high frequency acoustic waves are strongly attenuated in tissue. As such, the
current sensing mirror is suitable for imaging depths beyond a few millimeters
with no reduction in attainable resolution.

In summary, we have demonstrated the use of a Michelson interferometer as a
photoacoustic detector. Quadrature phase detection removes the need to tune
the sensitivity of the interferometer, allowing for the possibility of
simultaneous full surface detection. We used the detector to produce proof-of-
concept photoacoustic images with sub-millimeter resolution, and suggested
ways that resolution could be improved. Future work will involve demonstrating
simultaneous full surface measurement, with the aim of producing real-time 3D
photoacoustic visualizations of the sub-cranial blood vessel vasculature. This
would allow stroke researchers to visualize the dynamics of reperfusion in
small animal subjects after inducing a stroke, in both high resolution and on
a usefully short timescale.

Please cite the Applied Physics Letters version of this article, available
online at: http://link.aip.org/link/doi/10.1063/1.4816427

## References

  * Minghua and Lihong V. (2006) X. Minghua and W. Lihong V., Rev. Sci. Instrum., 77, 041101 (2006).
  * Li and Wang (2009) C. Li and L. V. Wang, Phys. Med. Biol., 54, R59 (2009).
  * Wang and Hu (2012) L. V. Wang and S. Hu, Science, 335, 1458 (2012).
  * Joshi and Agarwal (2010) S. Joshi and S. Agarwal, Ann. N.Y. Acad. Sci., 1199, 149 (2010).
  * Esenaliev  _et al._ (1999) R. Esenaliev, A. Karabutov, and A. Oraevsky, IEEE J. Sel. Top. Quant. Electron., 5, 981 (1999).
  * Oraevsky  _et al._ (2001) A. A. Oraevsky, A. A. Karabutov, S. V. Solomatin, E. V. Savateeva, V. A. Andreev, Z. Gatalica, H. Singh, and R. D. Fleming, Proc. SPIE, 4256, 6 (2001).
  * Yang  _et al._ (2012) Y. Yang, S. Wang, C. Tao, X. Wang, and X. Liu, Appl. Phys. Lett., 101, 034105 (2012).
  * Wang (2008) L. V. Wang, Med. Phys., 35, 5758 (2008).
  * Xiang  _et al._ (2013) L. Xiang, B. Wang, L. Ji, and H. Jiang, Sci. Rep., 3 (2013).
  * Hou  _et al._ (2007) Y. Hou, J.-S. Kim, S. Ashkenazi, S.-W. Huang, L. J. Guo, and M. ODonnell, Appl. Phys. Lett., 91, 073507 (2007).
  * Beard and Mills (1996) P. C. Beard and T. N. Mills, Appl. Opt., 35, 663 (1996).
  * Hamilton and O’Donnell (1998) J. Hamilton and M. O’Donnell, IEEE Trans. Ultrason., Ferroelectr., Freq. Control, 45, 216 (1998).
  * Paltauf  _et al._ (2007) G. Paltauf, R. Nuster, M. Haltmeier, and P. Burgholzer, Appl. Opt., 46, 3352 (2007).
  * Chow  _et al._ (2011) C. M. Chow, Y. Zhou, Y. Guo, T. B. Norris, X. Wang, C. X. Deng, and J. Y. Ye, J. Biomed. Opt., 16, 017001 (2011).
  * Zhang  _et al._ (2008) E. Zhang, J. Laufer, and P. Beard, Appl. Opt., 47, 561 (2008).
  * Lamont and Beard (2006) M. Lamont and P. Beard, Electron. Lett., 42, 187 (2006).
  * Huang  _et al._ (2008) S.-W. Huang, S.-L. Chen, T. Ling, A. Maxwell, M. ODonnell, L. J. Guo, and S. Ashkenazi, Appl. Phys. Lett., 92, 193509 (2008).
  * Paltauf  _et al._ (1999) G. Paltauf, H. Schmidt-Kloiber, K. P. Köstli, and M. Frenz, Appl. Phys. Lett., 75, 1048 (1999).
  * Paltauf and Schmidt-Kloiber (1997) G. Paltauf and H. Schmidt-Kloiber, J. Appl. Phys., 82, 1525 (1997).
  * Blitz (1963) J. Blitz,  _Fundamentals Of Ultrasonics_ (Butterworths Co., London, 1963).
  * Mezrich  _et al._ (1976) R. Mezrich, D. Vilkomerson, and K. Etzold, Appl. Opt., 15, 1499 (1976).
  * Shutilov (1988) V. A. Shutilov,  _Fundamental Physics of Ultrasound_ (Gordon and Breach Science Publishers, New York, 1988).
  * ANS (2007) _ANSI Z136.1_ , American National Standards Institute (2007).
  * Kostli  _et al._ (2001) K. P. Kostli, M. Frenz, H. Bebie, and H. P. Weber, Phys. Med. Biol., 46, 1863 (2001).
  * Treeby and Cox (2010) B. E. Treeby and B. T. Cox, J. Biomed. Opt., 15, 021314 (2010).


# Mass and Color Dependence of the Hubble Spiral Sequence

Petra Mengistu  Haverford College
370 Lancaster Avenue
Haverford, PA 19041, USA  Karen L. Masters  Haverford College
370 Lancaster Avenue
Haverford, PA 19041, USA

###### Abstract

In the classic Hubble spiral sequence, arm windiness correlates with bulge
size; Sa type spiral galaxies with larger bulges also have the most tightly
wound spirals. Exceptions to this have long been known, and in recent work
using Galaxy Zoo morphologies no strong correlation was seen in a volume
limited sample. In this Research Note, we explore the impact of galaxy mass
and integrated color upon this correlation in the Galaxy Zoo sample, finding
that bluer and lower mass spirals show the “expected” correlation; however, it
becomes slightly negative for redder and/or more massive spiral galaxies.

##  1 Introduction

Hubble’s classification scheme (Hubble, 1926) classifies galaxies into two
major categories: elliptical and spirals. Spirals are ordered in increasing
order of spiral arm windiness, “concentration” of the arms and dominance of
the central bulge size. From this, it has been implicitly assumed that arm
windiness positively correlates with galaxy bulge size and some work has
measured this (e.g. Kennicutt & Hodge, 1982). The correlation motivated, in
part, early models of static density wave spiral arms (i.e. Lin & Shu, 1964)
which predict tighter arms in more centrally concentrated galaxies. However,
it has long been known that the diversity of spiral galaxy morphologies is not
captured by the Hubble sequence (e.g. see Buta 2013), and current
understanding of spiral arms has moved on from static density waves (e.g. see
Sellwood & Masters 2022). Using a volume limited sample of galaxies selected
from the Galaxy Zoo 2 (GZ2) survey, and defining an arm winding and bulge size
score, Masters et al. (2019) found no significant correlation.

In this Research Note, we seek to better understand how spiral arm winding and
bulge size are related in sub-samples of spirals separated by global
properties (namely mass and colour).

##  2 Data and Sample Selection

Our base sample is GZ2 (Willett et al., 2013) – 250,000 nearby galaxies with
Sloan Digital Sky Survey (SDSS) Legacy imaging and spectroscopy (York et al.,
2000).

We limit to $z<0.035$ and $\log(M_{\star}/M_{\odot})>9.0$, for a volume-
limited sample of $N=$21,899. We use debiased vote fractions,
‘p[primaryresponse]’, from GZ2 (Willett et al., 2013)

to select 7,272 spiral galaxies using $p_{\rm features}>0.430$; $p_{\rm
notedge\\_on}>0.715$, and $p_{\rm spiral}>0.5$.

We use a previously defined arm winding score, $w_{\rm avg}=0.5p_{\rm
mediumspiral}+1.0p_{\rm tightspiral}$, and bulge prominence score, $b_{\rm
avg}=0.2p_{\rm justnoticeablebulge}+0.8p_{\rm obviousbulge}+1.0p_{\rm
dominantbulge}$. These values have been previously defined and are discussed
in Hart et al. (2017); Masters et al. (2019). The choice of coefficients is
arbitrary; these values produce scores ranging from 0–1 and were calibrated
against other measures of pitch angle and bulge size in Hart et al. (2017).

We make four subsets of spiral galaxies based on mass and global color as
illustrated in Figure 1. We use NUV-$r$ colours and stellar mass from the NASA
Sloan Atlas (NSA, Blanton et al., 2005, namely ELPETRO$\\_$MASS). We cut at
(NUV-$r$) = 4 and $\log(M_{\star}/M_{\odot})=10.3$ motivated by the typical
colour and mass of the blue cloud–red sequence transition (e.g. Chilingarian &
Zolotukhin 2012; Smethurst et al. 2022).

Overall, we have a total of 7,272 spiral galaxies divided into 5,999 low mass
blue, 758 massive blue, 297 low mass red and 218 massive red spirals.

##  3 Results

Figure 1 shows our results. The histograms show a notable difference in the
distribution of bulge sizes for different mass and color subsets. Lower mass
blue spirals have smaller bulges, while red and/or high mass spirals have
larger bulges, especially massive red spirals. All subsets exhibit a similar
skew in the distribution of arm winding scores towards tightly wound; massive
spirals are slightly more tightly wound than lower mass spirals.

The trend between spiral arm winding and bulge size is shown for all subsets
in the middle of Figure 1. We observe no strong correlation in any sub-set. In
low mass ($\log(M_{\star}/M_{\odot}<10.3$), blue spirals (top panel of Figure
1), we observe the positive correlation ($p=0.092$) noted by Hubble (Hubble,
1926). For the $\sim$%20 of spiral galaxies that are red and/or massive
($\log(M_{\star}/M_{\odot})>10.3$), the correlation becomes negative
($p$=values from -0.045 to -0.174) – looser arms appear mildly associated with
larger bulges. This is most significant in the subset of massive red spirals
(lower panel of 1), which have the steepest negative slope.

Figure 1: [Middle]: Two subplots showing winding score, wavg against bulge
size, bavg for subsets of different mass and color spirals. The small points
show the raw data, to which we have fitted a trend line; we also show fifteen
equally spaced binned averages - error bars show standard error on the mean.
[Upper middle]: Low mass $\log(M_{\star}/M_{\odot}<10.3)$ sample of galaxies
divided into blue (Ngal=5999), and red (Ngal=297) spirals (dark cyan and pink
data points and lines). [Lower middle]: High mass
$\log(M_{\star}/M_{\odot}>10.3)$ subset subdivided into blue (Ngal=758), and
red (Ngal=218) spirals (blue and orange data points and lines). Histograms of
bulge size and winding score are also shown for all subsets. [Side]: Example
images of tightly and loosely wound spiral galaxies in each subset; blue
spirals on the right, red at left, low mass at upper, and high mass lower.

##  4 Discussion and Conclusion

We find that the observed correlation between arm windiness and bulge size in
spiral galaxies has a significant mass and color dependence. Only low mass
blue spirals show clear evidence for a Hubble sequence like correlation
between bulge size and spiral pitch angle.

The photographic plate images Hubble used would be more sensitive to bluer
wavelengths and typical samples of very nearby galaxies skew towards lower
mass, while SDSS contains both redder and more massive spirals.

The superposition of both low-mass blue and redder/more massive spirals
perhaps explains the lack of correlation noted by Masters et al. (2019). It is
clear that sample selection (along with methods of determining pitch angle)
may play an important role in what correlation is observed.

Our findings suggest that spiral galaxies with different masses and colors,
properties which correlate with varying disc–halo mass fraction, or local
environment (Bamford et al. 2009), may have different mechanisms of spiral arm
formation. Spiral arm formation mechanisms can invoke tidal triggering, links
to galactic bars, or self-excitation linked to shear (Sellwood & Masters,
2022). Redder massive spirals are more likely to be found in higher densities
(e.g. Masters et al., 2010), and show a greater diversity of rotation curves
shapes (i.e. shear; Frosst et al. 2022). Treating all spiral galaxies as
having the same spiral formation mechanism may obscure the physical process.

## References

  * Bamford et al. (2009) Bamford, S. P., Nichol, R. C., Baldry, I. K., et al. 2009, MNRAS, 393, 1324, doi: 10.1111/j.1365-2966.2008.14252.x
  * Blanton et al. (2005) Blanton, M. R., Schlegel, D. J., Strauss, M. A., et al. 2005, AJ, 129, 2562, doi: 10.1086/429803
  * Buta (2013) Buta, R. J. 2013, in Planets, Stars and Stellar Systems. Volume 6: Extragalactic Astronomy and Cosmology, ed. T. D. Oswalt & W. C. Keel, Vol. 6, 1, doi: 10.1007/978-94-007-5609-0_1
  * Chilingarian & Zolotukhin (2012) Chilingarian, I. V., & Zolotukhin, I. Y. 2012, MNRAS, 419, 1727, doi: 10.1111/j.1365-2966.2011.19837.x
  * Frosst et al. (2022) Frosst, M., Courteau, S., Arora, N., et al. 2022, MNRAS, 514, 3510, doi: 10.1093/mnras/stac1497
  * Hart et al. (2017) Hart, R. E., Bamford, S. P., Hayes, W. B., et al. 2017, MNRAS, 472, 2263, doi: 10.1093/mnras/stx2137
  * Hubble (1926) Hubble, E. P. 1926, ApJ, 64, 321, doi: 10.1086/143018
  * Kennicutt & Hodge (1982) Kennicutt, R., J., & Hodge, P. 1982, ApJ, 253, 101, doi: 10.1086/159614
  * Lin & Shu (1964) Lin, C. C., & Shu, F. H. 1964, ApJ, 140, 646, doi: 10.1086/147955
  * Masters et al. (2010) Masters, K. L., Mosleh, M., Romer, A. K., et al. 2010, MNRAS, 405, 783, doi: 10.1111/j.1365-2966.2010.16503.x
  * Masters et al. (2019) Masters, K. L., Lintott, C. J., Hart, R. E., et al. 2019, MNRAS, 487, 1808, doi: 10.1093/mnras/stz1153
  * Sellwood & Masters (2022) Sellwood, J. A., & Masters, K. L. 2022, ARA&A, 60, doi: 10.1146/annurev-astro-052920-104505
  * Smethurst et al. (2022) Smethurst, R. J., Masters, K. L., Simmons, B. D., et al. 2022, MNRAS, 510, 4126, doi: 10.1093/mnras/stab3607
  * Willett et al. (2013) Willett, K. W., Lintott, C. J., Bamford, S. P., et al. 2013, MNRAS, 435, 2835, doi: 10.1093/mnras/stt1458
  * York et al. (2000) York, D. G., Adelman, J., Anderson, Jr., J. E., et al. 2000, AJ, 120, 1579, doi: 10.1086/301513


# Hierarchical Reinforcement Learning with Opponent Modeling for Distributed
Multi-agent Cooperation

Zhixuan Liang, Jiannong Cao, Fellow, IEEE, Shan Jiang§, Divya Saxena, Huafeng
Xu
Department of Computing, The Hong Kong Polytechnic University

###### Abstract

Many real-world applications can be formulated as multi-agent cooperation
problems, such as network packet routing and coordination of autonomous
vehicles. The emergence of deep reinforcement learning (DRL) provides a
promising approach for multi-agent cooperation through the interaction of the
agents and environments. However, traditional DRL solutions suffer from the
high dimensions of multiple agents with continuous action space during policy
search. Besides, the dynamicity of agents’ policies makes the training non-
stationary. To tackle the issues, we propose a hierarchical reinforcement
learning approach with high-level decision-making and low-level individual
control for efficient policy search. In particular, the cooperation of
multiple agents can be learned in high-level discrete action space
efficiently. At the same time, the low-level individual control can be reduced
to single-agent reinforcement learning. In addition to hierarchical
reinforcement learning, we propose an opponent modeling network to model other
agents’ policies during the learning process. In contrast to end-to-end DRL
approaches, our approach reduces the learning complexity by decomposing the
overall task into sub-tasks in a hierarchical way. To evaluate the efficiency
of our approach, we conduct a real-world case study in the cooperative lane
change scenario. Both simulation and real-world experiments show the
superiority of our approach in the collision rate and convergence speed.

###### Index Terms:

Multi-agent Cooperation; Deep Reinforcement Learning; Hierarchical
Reinforcement Learning

§§footnotetext: Corresponding author: Shan Jiang

##  I Introduction

Many complex real-world applications can be modeled as multi-agent
cooperation, such as network packet routing [1], energy distribution [2], and
coordination of autonomous vehicles [3]. In these applications, the agents
need to make individual decisions considering the mutual influence.
Traditional solutions of multi-agent cooperation are mostly model-based, in
which the agents and their interactions are modeled using physical formulas
and prior knowledge [4][5]. These approaches fail to adapt to dynamic,
stochastic, and complex environments. Recently, the development of deep
reinforcement learning (DRL) [6] provides a promising solution through a
trial-and-error process [7]. At each step, the agent observes the environment,
selects the optimal action, and receives rewards as feedback signals related
to the team performance. The goal of each agent is to learn the policies
maximizing the accumulated rewards received from the environment.

Figure 1: Illustration of hierarchical reinforcement learning for distributed
multi-agent cooperation. Each agent maintains a high-level cooperation layer
and a low-level individual control layer.

Existing DRL approaches for multi-agent cooperation can be classified into
three categories: centralized reinforcement learning, centralized training
with decentralized execution, and distributed reinforcement learning.
Centralized reinforcement learning aims to train a centralized value function
that selects actions for all the agents. Such a method is hard to be extended
to large-scale scenarios due to the exponential action spaces when the number
of agents increases. An alternative approach is centralized training with
decentralized execution (CTDE), which trains a critic network to estimate
state-action pairs’ values. At the same time, each agent maintains an
individual actor network for decentralized action selection. By training such
a critic network considering the states and actions of other agents, all the
agents can learn to cooperate in certain states. Some previous works have
applied CTDE for real-world applications, such as flocking [8] and pathfinding
[9]. However, the number of features in the critic network needs to be scaled
up linearly (in the best case) or exponentially (in the worse case) as the
number of agents increases. Besides, the gain of learning a centralized critic
is likely to be minimal in the sparse interaction scenarios [10].

Recent developments in DRL address the above limitations through decentralized
training with decentralized execution (DTDE), where each agent acts and learns
to cooperate in a distributed manner[11][12][13]. However, it is non-trivial
to learn in a distributed way because the training is non-stationary, and it
is hard to represent individual agents’ policies. In particular, standard DRL
approaches relying on end-to-end models suffer from the curse of
dimensionality of multiple agents in continuous action spaces. It demands
large-scale models and a long time for training to learn a mapping function
from the state to the continuous action space. Furthermore, the independence
of individual policy learning leads to poor coordination. Finally, the
dynamicity of the agents’ policies makes it nearly impossible to employ
experience replay [14] that is crucial for stabilizing DRL.

In this paper, we propose a novel hierarchical deep reinforcement learning
approach for distributed multi-agent cooperation. First, we decompose the
policy space of each agent into a high-level cooperation layer and a low-level
individual control layer, which is shown in Fig. 1. The cooperation of
multiple agents can be efficiently learned in high-level discrete action
space, while the low-level individual control can be reduced to independent
reinforcement learning. Besides, we introduce an opponent modeling mechanism
to model other agents’ high-level decisions. The learned opponent models will
encourage cooperation behaviors and stabilize DRL. In contrast to the standard
end-to-end DRL model, we reduce the learning complexity by decomposing the
overall task into sub-tasks that are easier to solve hierarchically.
Furthermore, we conduct a case study on cooperative lane change and present
the hierarchical decision-making model with the specific design of states,
actions, and rewards. Finally, we evaluate our approach in a simulation
environment and a real-world testbed. The main contributions of this paper are
as follows:

  * •

We study the problem of distributed multi-agent cooperation in continuous
action space and propose a hierarchical deep reinforcement learning approach,
which decomposes the overall task into high-level option selection and low-
level individual control.

  * •

To address the non-stationary training issue in distributed multi-agent
learning, we introduce an opponent modeling mechanism to model other agents’
high-level policies.

  * •

We conduct extensive experiments on cooperative lane change and compare our
approach with state-of-art baselines in simulation and a real-world testbed.
The experimental results show the superiority of our approach in both task
performance and convergence speed.

The rest of the paper is organized as follows. Sec. VI summarizes the previous
works of multi-agent reinforcement learning and hierarchical reinforcement
learning. Sec. II introduces the notations and preliminaries of deep
reinforcement learning. Sec. III presents a novel hierarchical reinforcement
learning algorithm with an opponent modeling mechanism. Sec. IV showcases a
real-world case study of hierarchical reinforcement learning on the multi-
vehicle cooperative lance change scenario. Sec. V shows the experimental
results in simulation and a real-world bed. Finally, Sec. VII concludes the
paper and discusses the future directions.

##  II Preliminaries

In this section, we introduce the definitions of Markov decision process
(MDP), reinforcement learning (RL), and hierarchical reinforcement learning
(HRL).

###  II-A Markov Decision Process

Many decision-making problems can be mathematically modeled as a Markov
decision process (MDP). An MDP is defined as a four-tuple $(S,A,r,T)$, in
which $S$ is the set of states, $A$ is the set of available actions,
$r:S\times A\times S\rightarrow R$ is the reward function, and
$\mathcal{T}:S\times A\times S\times\rightarrow[0,1]$ is the probability
function of state-action-state transition. A stochastic policy function is
defined as $\pi:S\times A\rightarrow[0,1]$ and a deterministic policy function
$\mu:S\rightarrow R^{|A|}$ is defined to select an action for current state.
The goal of solving MDP is to find a policy which selects actions maximizing
the cumulative rewards $R=\sum_{t=0}^{N}r_{t}$.

Given a policy $\pi$ and a state $s$, the value function $V:S\rightarrow R$
calculates the expected sum of discounted rewards using the following formula:

$\displaystyle
V^{\pi}(s)=\mathbb{E}_{\pi}\left[\sum\nolimits_{t=0}^{N}\gamma^{t}r_{t}\mid
s\right]$

where $\gamma$ is the discount factor and $r_{t}$ is the reward received at
each time step. In addition, given a policy $\pi$, a state $s$, and an action
$a$, the state-action value function $Q:S\times A\rightarrow R$ calculates the
expected sum of discounted rewards from a given state as follows:

$\displaystyle
Q^{\pi}(s,a)=\sum\nolimits_{\pi}\left[\sum\nolimits_{t=0}^{N}\gamma^{t}r_{t}\mid
s,a\right]$

Solving an MDP requires computing the state-value function and state-action
value function. The reward function and state transition are entirely known,
these functions can be solved by iterating the Bellman equations, and the
decision-making problem can be solved by using dynamic programming [15].
However, the reward function and state transition are unknown to agents in
most cases.

###  II-B Reinforcement learning

Reinforcement learning provides multiple methods for solving MDPs, classified
as value-based and policy-based methods.

In value-based method, the state-action-value function (or state-value
function) is estimated using temporal different (TD) methods as follows:

$\displaystyle\begin{split}&Q(s_{t},a_{t})\\\
=&Q(s_{t},a_{t})+\alpha(r_{t}+\gamma\max_{a^{t+1}}Q(s^{t+1},a^{t+1})-Q(s_{t},a_{t}))\end{split}$

Then the greedy policy can be derived by selecting the action with highest
Q-value at each time step as follows:

$\displaystyle\pi(s_{t})=\arg\max_{a}Q(s_{t},a_{t})$

Such a method is referred to as the Q-learning [16].

Traditional tabular Q-learning suffers from the issue of high dimensionality
in state space. Deep Q-learning (DQN) [17] addresses the issue by employing a
deep neural network to approximate the Q-function. Besides, DQN is a class of
off-policy reinforcement learning methods that uses a replay memory to store
the transition four-tuple $(s,a,r,s^{\prime})$ and these data can be sampled
to train the Q-network. The loss function of training Q-network is as follows:

$\displaystyle\mathcal{L}(\theta)=\mathbb{E}_{s,a,r,s_{{}^{\prime}}^{\prime}}\left[\left(Q\left(s,{a};\theta\right)-y\right)^{2}\right]$

where
$y=r+\gamma\max_{{a}\in\mathcal{A}}\left(Q\left({s}_{{}^{\prime}},{a};\theta^{-}\right)\right)$
and $\gamma\in[0,1)$ is the discount factor. In $y$, the variable $\theta^{-}$
is the parameter of the target network that is periodically copied from
$\theta$ that is kept constant for a number of iterations.

An alternative model-free approach is the policy-based method, which directly
learns a policy $\pi$ parameterized by $\theta^{\pi}$. The objective of these
approaches is to adjust the parameters $\theta^{\pi}$ in order to maximize the
function $J(\theta)=\mathbb{E}_{s\sim p^{\pi},a\sim\pi_{\theta}}[R]$, where
$R$ denotes the expected accumulative rewards which are usually approximated
by $Q(s,a)$. In the actor-critic theory, the policy network $\pi_{\theta}$ can
be considered an actor network, and the state-action value function $Q(s,a)$
can be considered a critic network. According to the policy gradient theorem
[15], the gradient of this objective function is defined as follows:

$\displaystyle\nabla_{\theta}J\left(\pi_{\theta}\right)=\mathbb{E}_{s\sim\rho^{\pi},a\sim\pi_{\theta}}\left[\nabla_{\theta}\log\pi_{\theta}(a|s)Q^{\pi}(s,a)\right]$

where $\rho^{\pi}$ is the state distribution under policy $\pi_{\theta}$ and
$\pi_{\theta}(a|s)$ is the probability of select the action $a$ under the
state $s$. At each time step, the agent samples an action from the
distribution generated from the policy network.

Deep deterministic policy gradient (DDPG) is an actor-critic algorithm
extended from stochastic policy gradient to deterministic policy gradient
[18][19]. It employs a deep neural network parameterized by $\theta^{\mu}$ to
approximate the deterministic policy $\mu_{\theta}:S\rightarrow A$ and a
neural network parameterized by $\theta^{Q}$ to approximate the action-value
function $Q(s,a|\theta^{Q})$. The critic network is learned using the Bellman
equation as in Q-learning, and the actor network is updated by applying the
chain rule to the objective function $J$:

$\displaystyle\nabla_{\theta}J(\theta)=\mathbb{E}_{s\sim\mathcal{D}}\left[\left.\nabla_{\theta}\mu_{\theta}(a\mid
s)\nabla_{a}Q^{\mu}(s,a)\right|_{a=\mu_{\theta}(s)}\right]$

###  II-C Hierarchical Reinforcement learning with Temporal Abstraction

Human decision-making often involves choosing among temporally extended
courses of action over a broad range of time scales [20]. Learning and
operating over different levels of temporal abstraction is a critical
challenge in tasks involving long-term planning. Sutton et al. [21] extended
reinforcement learning framework to include temporally abstract actions,
representations that group together a set of interrelated actions (for
example, moving to block $A$, driving another lane, passing the ball to
another person). These representations can be translated into a series of
individual actions, which are described as temporal abstraction. A recent
extension of this direction is hierarchical deep reinforcement learning with
temporal abstraction proposed by Kulkarni et al. [22]. In their approach, the
agent uses a two-level decision-making model, i.e., a meta-controller and a
controller. The meta-controller receives a state $s_{t}$ and chooses a goal
$g_{t}\in\mathcal{G}$, where $\mathcal{G}$ denotes the set of all possible
current goals. The $g_{t}$ will remain for $T$ time steps, or the terminated
state is reached. The Q-value function for the controller is:

$\displaystyle
Q_{l}(s,a;g)=\mathrm{E}\left[r_{l}+\gamma\max_{a_{{}^{\prime}}}Q_{l}\left(s_{{}^{\prime}},a_{{}^{\prime}};g\right)\right]$

where $g$ is the goal and $r_{l}$ is low-level reward given by a goal-driven
reward function $R(s_{t},a_{t},g_{t})$. This low-level reward is related to
performing the goal (temporal abstraction), which is also regards as intrinsic
reward. The design of intrinsic reward and termination is still an open
question in hierarchical reinforcement learning.

Similarly, the Q-value function of the meta-controller is:

$\displaystyle
Q(s,g)=\mathrm{E}[r_{h}+\gamma\max_{g^{{}^{\prime}}}Q(s^{{}^{\prime}},g^{{}^{\prime}})]$

where $r_{h}=\sum_{t=0}^{T}r$ is the accumulative reward received from the
environment during the $T$ time steps.

##  III HERO: Hierarchical Reinforcement Learning with Opponent Modeling

This section introduces the sequential decision-making problem in distributed
multi-agent systems. Then, we propose HERO, a general hierarchical decision-
making model with high-level cooperative decision-making and low-level
individual control. Note that, in many applications, the primitive actions of
agents are regarded as individual control.

###  III-A Problem Formulation

Typically, multi-agent cooperation can be mathematically modeled as multi-
agent Markov games, which extends MDP to multi-agent setting [23]. A Markov
game is defined as a five-tuple $(I,S,A,r,T)$, where $I$ is the set of $N$
agents, $S$ is the set of states, $A=A_{1}\times A_{2}\ldots\times A_{N}$ is
the set of actions of all agents, and
$r=\left(r_{1},r_{2},\ldots,r_{N}\right)$ where $r_{i}:S\times A\times
S\rightarrow\mathbb{R}$ is the reward function of agent $i$. In the fully
cooperative setting, we have $r_{1}=r_{2},...,=r_{N}$. Given the current state
and the actions of all agents, $T:S\times A\times S\rightarrow[0,1]$ is the
probability distribution over the next states. Each agent $i$ aims to learn a
policy $\pi_{i}:S_{i}\times A_{i}\rightarrow[0,1]$ to select the optimal
action $a_{i}$ maximizing the accumulative rewards
$R_{i}=\sum_{t=0}^{T}\gamma^{t}r_{i}^{t}$, where $\gamma$ is the is the
discount factor and $T$ is the time horizon. We can formulate the joint policy
of other agents as:

$\displaystyle\pi^{-i}\left(a_{t}^{-i}\mid
s_{t}\right)=\prod\nolimits_{j\in\\{-i\\}}\pi^{j}\left(a_{t}^{j}\mid
s_{t}\right)$

Therefore, the objective of each agent $i$ is defined as follows:

$\displaystyle\begin{split}&\max_{\pi^{i}}\eta_{i}\left[\pi^{i},\pi^{-i}\right]\\\
=&\mathbb{E}_{\left(s_{t},a_{t}^{i},a_{t}^{-i}\right)\sim
T,\pi^{i},\pi^{-i}}\left[\sum\nolimits_{t=1}^{\infty}\gamma^{t}R^{i}\left(s_{t},a_{t}^{i},a_{t}^{-i}\right)\right]\end{split}$

Most existing multi-agent reinforcement learning methods assume that each
agent can access all the other agents’ policies during the training. In this
paper, we consider a general scenario in a distributed manner, in which each
agent $i$ has no knowledge of other agents’ policies, and can only observe the
historical states and actions of other agents, i.e.,
$\left\\{s_{1:t-1},a_{1:t-1}^{-i}\right\\}$. Such a setting is also found in
other decentralized MARL works [11][13].

###  III-B Hierarchical Decision-making Model for Multi-agent Cooperation

Figure 2: Two-stage training structure of HERO. (a) Each individual agent
learns different individual control policy with random noise in the first
stage. (b) Multiple agents learn to select options in the second stage.

Hierarchical Policy Representation. In this paper, we introduce a hierarchical
reinforcement learning model for multi-agent cooperation, which decomposes an
overall cooperation task into a hierarchy of discrete sub-tasks. The sub-tasks
are also regarded as options, skills, or temporal abstractions in previous
studies. For simplicity, we use the notation of option in the rest of this
paper. In a two-layer hierarchical reinforcement learning model, the policy
space of each agent is defined as follows:

$\displaystyle\pi^{i}=\left[\pi^{i}_{h},\pi^{i}_{l}\right]$

where $\pi_{h}$ is the high-level policy to select the options and $\pi_{l}$
is the low-level policy to select primitive actions to perform the selected
option.

A high-level option is defined as a three-tuple $o=(I_{o},\pi_{h},\beta_{o})$,
where ${I_{o}}\in\mathcal{S}_{h}$ is the initiation set, $\pi_{h}$: the option
selection policy, and $\beta_{o}:\mathcal{\delta}\rightarrow[0,1]$ specifies
the termination condition while executing $o$.

In the definition of high-level option, $S_{h}$ denotes the state space, and
$A_{h}$ denotes the action space of high-level policy. The option $o$ can also
be interpreted as the high-level action with the state set $S_{h}$. In the
context of multi-agent cooperation, each agent $i$ cooperatively selects the
option based on the current state $s^{i}_{h}$ and the inferred opponent
options $o^{-i}_{h}$, which is shown in Fig. 2.

Asynchronous Option Termination. There are two modes for option termination:
synchronous and asynchronous. Synchronous termination requires all agents to
interrupt the existing option execution and select the next option
synchronously, which is not feasible for fully distributed systems. Thus, we
consider the asynchronous termination mode for each agent. At each time step
$t$, the high-level layer of each agent $i$ will verify whether the current
state $s^{i}_{h,t}$ satisfies the termination condition. In particular, if
$\beta^{i}_{o}$ is equal to $1$, the existing option will terminate, and the
high-level policy will select another option to execute until the overall task
is finished.

Hierarchical Reward. Also, the reward of each agent $i$ can be divided into
two parts for efficient training:

$\displaystyle r^{i}=\left[r^{i}_{h},r^{i}_{l}\right]$

where $r^{i}_{h}$ and $r^{i}_{l}$ are the high-level and low-level rewards
concerning team performance and individual control, respectively.

###  III-C Learning High-level Option Selection with Opponent Modeling

Figure 3: Illustration of the high-level opponent modeling in high-level
layer. Each agent maintain a self policy network for its option selection and
a opponent modeling network for other agents’ option prediction.

Option-value Function. The value function of each agent’s option selection can
be defined as follows:

$\displaystyle Q^{i}_{h}\left(s^{i}_{h,t},o_{t}\right)=r^{i}_{h,t:t+c}+\gamma
Q^{i}_{h}\left(s^{i}_{t+c},o_{t+c}\right)-Q^{i}_{h}\left(s^{i}_{h,t},o_{t}\right)$

where $\mathbf{o_{t}}=[o^{i}_{t},o^{-i}_{t}]$ denotes the option selected by
agent $i$ and the options selected by other agents. Here, we assume that each
agent can observe all the previous option selections of others agents, i.e.,
$\left\\{o^{-i}_{1:t-1}\right\\}$, at each time step. Besides, $r^{i}_{h,t}$
is the high-level reward received from environment which indicates the
progress of the overall task and $r^{i}_{h,t:t+c}=\sum^{t+c}_{t}r^{i}_{h,t}$
denotes the accumulated high-level reward when performs the option
$o^{i}_{t}$. Thereafter, the high-level transition that will be stored in the
replay buffer $\mathcal{D}^{i}_{h}$ is expressed as follows:

$\displaystyle\\{(s^{i}_{h,t+k},o^{i}_{t+k},o^{-i}_{t+k},r^{i}_{h,t+k},s^{i}_{h,t+k+1})\\}_{k=0}^{k=c}$

Algorithm 1 Training High-level Cooperative Strategy with Opponent Modeling

1: initialize the experience replay buffer $\mathcal{D}^{i}_{h}$ and the
parameters $\theta^{Q}_{i,h}$, $\theta^{\pi}_{i,h}$ for each agent i

2: initialize the experience replay buffer $\mathcal{D}^{-i}_{h}$ and the
parameters $\pi^{-i}_{h}$ for opponent modeling

3: for episode $\leftarrow 1$ to $M$ do

4: t $\leftarrow$ 0

5: reset the environment and receive an initial states $s^{i}_{h,t}$ and
$s^{i}_{l,t}$ for each agent $i$

6: for each agent $i$ select high-level action (option) from the actor network
$\pi^{i}_{h}$

7: $o^{i}_{t}\leftarrow\pi^{i}_{h}$

8: while the task is not completed do

9: pass the high-level option $o^{i}_{t}$ to low-level controller and generate
low-level action $a^{i}_{t}$

10: execute action $a^{i}_{t}$ and receive next state $s^{i}_{h,t+1}$

11: compute the high-level team reward $r^{i}_{h}$

12: for agent $i\leftarrow 1$ to $N$ do

13: if $s^{i}_{h,t+1}$ is a terminal state for option $o^{i}_{t}$ then

14: choose next option $o^{i}_{t+1}$ from actor network $\pi^{i}_{h}$

15: store transition ($s^{i}_{h,t},o^{i}_{t},r^{i}_{h},s^{i}_{h,t+1}$) in
$\mathcal{D}^{i}_{h}$

16: $o^{i}_{t}\leftarrow o^{i}_{t+1}$, $s^{i}_{h,t}\leftarrow s^{i}_{h,t+1}$

17: else

18: $o^{i}_{t}\leftarrow o^{i}_{t+1}$

19: end if

20: update the high-level critic network

21: update the high-level actor network

22: end for

23: store transition $(s^{i}_{t},o^{-i}_{t})$ in $\mathcal{D}^{-i}_{h}$

24: t $\leftarrow t+1$, $s_{t}\leftarrow s_{t+1}$

25: update the opponent modeling networks

26: end while

27: update the target network parameters

28: end for

Algorithm 2 Training Low-level Driving Skills with Intrinsic Reward Functions

1: initialize the experience replay buffer $\mathcal{D}_{l}$ and the
parameters $\theta^{l}_{Q}$, $\theta^{l}_{\mu}$ for low-level controller

2: for episode $\leftarrow 1$ to $M$ do

3: reset the environment and receive an initial state $s_{0}$

4: t $\leftarrow$ 0

5: while $s_{t}$ is not a terminal state do

6: select action from low-level actor network

7: $a_{t}\leftarrow\pi_{l}(\cdot\mid\mathbf{s}_{t})$

8: execute action $a_{t}$ and receive next state $s_{t+1}$ and intrinsic
reward $r_{l}$

9: store transition $(s_{t},a_{t},r_{t},s_{t+1})$ in $\mathcal{D}_{l}$

10: t $\leftarrow t+1$

11: $s_{t}\leftarrow s_{t+1}$

12: update the low-level critic network

13: update the low-level actor network

14: end while

15: update the target network parameters

16: end for

Agent Learning. We propose a decentralized actor-critic method for high-level
policy learning. Each agent trains a decentralized critic network
$Q^{i}_{h}\left(s^{i}_{h},o^{i},o^{-i};\theta^{Q}_{i,h}\right)$ to approximate
the option-value function, while maintaining an actor network
$\pi^{i}_{h}(s^{i}_{h},o^{-i};\theta^{\pi}_{i,h})$ for option generation.
Specifically, the critic network is trained by minimizing the loss:

$\displaystyle\mathcal{L}\left(\theta^{Q}_{i,h}\right)=\mathbb{E}_{{s^{i}_{h}},o^{i},o^{-i},r^{i}_{h}\sim\mathcal{D}^{i}_{h}}\left[\left(Q^{i}_{h,t}\left({s^{i}_{h,t}},o^{i}_{t},o^{-i}_{t}\right)-y^{i}_{h}\right)^{2}\right]$
$\displaystyle y^{i}_{h}=r^{i}_{h,t:t+c}+\gamma
Q^{i}_{h,t+c}\left(s^{i},o^{i}_{t+c},\hat{o}^{-i}_{t+c})\right)$

where $\mathcal{D}^{i}_{h}$ is the experience replay buffer for offline high-
level policy training. Then, the actor network can be trained with gradient
ascent where the gradient is computed as:

$\displaystyle\begin{split}&\nabla_{\theta^{\pi}_{i,h}}J\left(\theta^{\pi}_{i,h}\right)=\\\
&\mathbb{E}_{s^{i}_{h}\sim
p^{\pi},o^{i}\sim\pi^{i}_{h}}\left[\nabla_{\theta_{i}}\log\pi^{i}_{h}\left(o^{i}\mid
s^{i}_{h},\hat{o}^{-i}\right)Q^{i}_{h}\left(s^{i}_{h},o^{i},o^{-i}\right)\right]\end{split}$

where $p^{\pi}$ denotes the state transition under all agents’ policies
$\pi={\pi_{1},...,\pi_{n}}$ and $\hat{o}^{-i}$ denotes the inferred options of
other agents generated from the opponent modeling network
$\hat{\pi}^{-i}_{h}$.

Opponent-modeling in Option Selection. Opponent modeling is a feasible
solution to the non-stationary issue in distributed reinforcement learning.
Instead of modeling other agents’ primitive actions, we propose to model their
option selections reflecting their temporal abstraction during a few time
steps. Then, two questions arise: how to use the opponent model and how to
train the opponent model.

For the model usage, we use the opponent model in the individual option
selection to encourage the coordination behaviors. Each agent $i$ makes a
decision based on both the current state $s^{i}_{h}$ and the inferred opponent
options $\hat{o}^{-i}$, which is shown in Fig. 3. Additionally, we use the
latest opponent model to update the TD-target $y^{i}_{h}$:

$\displaystyle y^{i}_{h}=r^{i}_{h,t:t+c}+\gamma
Q^{i}_{h,t+c}\left(s^{i}_{h},\pi^{i}_{h}(s^{i}_{h,t+c}),\pi^{-i}_{h}(s^{i}_{h,t+c}))\right)$

Note that we input the option log probabilities of other agents directly into
$Q^{i}_{h,t+c}$, rather than sampling.

For the opponent model training, we use a deep neural network to approximate
the option selection policies of other agents and train the network by
maximizing the log probability from the recent observation histories:

$\displaystyle\mathcal{L}\left(\theta^{-i}_{\pi,h}\right)=-\mathbb{E}_{s^{i}_{h},o^{-i}_{h}}\left[\log\pi^{-i}_{h}\left(o^{-i}_{h}\mid
s^{i}_{h}\right)+\lambda H\left(\pi^{-i}_{h}\right)\right]$

where $\pi^{-i}_{h}\left(o^{-i}\mid s^{i}_{h}\right)$ is the predicted
probability of selected action $o^{-i}$ given the observation of agent $i$ and
$H\left(\pi^{-i}_{h}\right)$ is the entropy of the policy distribution which
is used to solve the over-fitting problem in deep learning.

###  III-D Learning Low-level Policies with Intrinsic Reward Functions

Given the option, the objective of the low-level layer is to perform the
option by selecting optimal primitive actions. The low-level individual
control policy $\pi^{i}_{l}$ is represented by an individual deep neural
network $\theta^{i}_{l}$ and trained with independent and parallel deep
reinforcement learning. To successfully learn to perform the option, the
intrinsic reward function $r^{i}_{l}$ is used for efficient training. In this
paper, we adopt the soft actor-critic method. According to the maximum entropy
reinforcement learning theorem [24], the objective function of low-level
policy is defined as follow:

$\displaystyle
J(\pi^{i}_{l})=\sum_{t=0}^{c}\mathbb{E}_{\left(s^{i}_{l},a^{i}\right)}\left[r^{i}_{l}\left(s^{i}_{l},a^{i}\right)+\alpha\mathcal{H}\left(\pi^{i}_{l}\left(\cdot\mid
s^{i}_{l}\right)\right)\right]$

where $\mathcal{H}$ is the regulation term to encourage exploration and learn
diverse policies to perform the specific option and $\alpha$ is the
hyperparameter denote the weight of action exploration. Consequently, the
gradient of each agent is:

$\displaystyle\begin{split}&\nabla_{\theta^{i}_{l}}J\left(\theta^{i}_{l}\right)=\mathbb{E}\left[\nabla_{\theta^{i}_{l}}\epsilon\left[-\alpha\epsilon+Q^{i}_{l}\left({s}^{i}_{l},{a}^{i};o^{i}_{h}\right)\right]\right]\\\
&\text{where }\epsilon=\log\boldsymbol{\pi}^{i}_{l}\left(a^{i}\mid
s^{i}_{l}\right)\end{split}$

In practice, the training of critic can be realized by parameter sharing among
distributed agents. Then, the critic network
$Q^{i}_{l}\left({s^{i}_{l}},{a^{i}};\theta^{i}_{l},o^{i}_{h}\right)$ can be
trained with deep Q-learning and the loss function is:

$\displaystyle\mathcal{L}(\theta^{i}_{l})=\mathbb{E}_{s^{i}_{l},a^{i},r^{i}_{l},s^{i,^{\prime}}_{l}\sim\mathcal{D}^{i}_{l}}\left[\left(Q^{i}_{l}\left(s^{i}_{l},a^{i};\theta^{i}_{l},o^{i}_{h}\right)-y^{i}_{l}\right)^{2}\right]$
$\displaystyle y^{i}_{l}=r^{i}_{l}+\gamma
Q^{i}_{l}\left(s^{i,^{\prime}},a^{{}^{\prime}}_{i};\theta^{i}_{l};o^{i}_{h}\right)$

where $D^{i}_{l}$ is the experience reply buffer for off-policy training and
$r^{i}_{l}$ is the intrinsic reward designed for the option learning.

Figure 4: The left side of this figure shows the reward calculation for the
training of driving in lane, while the right side shows the reward calculation
for training of lane change.

##  IV Case Study: HERO for Multi-vehicle Cooperative Lane Change

In this section, we present a case study of HERO on cooperative lane change to
demonstrate the practicability of HERO in real-world applications.

###  IV-A Cooperative Lane Change

Figure 5: Hierarchical decision-making model for each vehicle.

Lane change is one of the most essential and fundamental behaviors in
autonomous driving, which involves the interaction of multiple vehicles in
different lanes. Successful lane changes require drivers to account for
several safety-related factors, including the road geometries, positions of
ego vehicles, and cooperation with other vehicles. Inaccurate lane change
leads to car accidents, congestion, and waste of energy. This paper aims to
solve decision-making problems under different lane change scenarios. Many
existing works only consider the single-agent driving scenario, while this
work aims to leverage the cooperativeness of multiple vehicles to improve road
safety and efficiency.

In this paper, we propose a two-layer hierarchical decision-making model,
which includes a high-level controller and a low-level controller for
cooperative driving. As shown in Fig. 5, we first decompose the driving task
into several sub-tasks, i.e., lane change, lane-keeping, and slow down. Each
sub-task can be regarded as an option that the low-level controller can
perform. It takes several time steps for the low-level controller to finish
the option successfully. The objective of the high-level controller is to
learn to select different options according to the current state that can
maximize the traffic throughput while avoiding collisions. Specifically, the
high-level controller’s state, option, and reward function are defined in the
following subsection.

###  IV-B High-level State, Option and Reward Function Design

High-level State Space. Here, we separate the state space of each vehicle into
$s_{h}^{i}$ and $s_{l}^{i}$, which are the state for high-level controller and
low-level controller, respectively. The high-level state $s_{h}^{i}$ is
defined as follows:

$\displaystyle
s_{h}^{i}=\left[s_{lidar}^{i},s_{speed}^{i},s_{laneID}^{i}\right]$

where $s_{lidar}^{i}$ is the raw observation from the lidar sensor, which
denotes the distance with other vehicles from $360$ degrees, $s_{speed}^{i}$
is the vehicle’s speed, and $s_{laneID}^{i}$ is the vehicle’s current lane
identifier.

High-level Option Space. The action space of high-level controller is to
choose which option to execute. Here, we consider the discrete action space
defined as follow:

$\displaystyle A_{h}^{i}=\left[\text{keep lane},\text{slow
down},\text{accelerate},\text{lane change}\right]$

We define four types of high-level options according to human driving
behaviors, where each option corresponds to a learned low-level control
policy.

High-level Team Reward. The design of high-level team reward is related to the
overall objective. Typically, safety and efficiency are two main concerns in
cooperative merge scenarios. In terms of safety, each vehicle should not only
avoid collision with the front vehicle when driving in the current lane but
also pay attention to the vehicles in the merged lane. Once the collision
happens, we will assign a negative reward $r_{col}$ to the vehicles as the
penalty, and the current training episode will be ended and restarted. To
avoid traffic congestion, we encourage the vehicle to run as fast as possible
by giving a positive reward $r_{travel}$ to vehicles. Hence, the total reward
of each vehicle is designed as follows:

$\displaystyle r_{h}^{i}=\alpha r_{col}+(1-\alpha)r_{travel}^{i}$

where $r_{col}$ is the penalty of collision, $r_{travel}$ is the reward for
moving forward, and $\alpha$ is a hyperparameter to control the weight of
collision avoidance and move forward.

###  IV-C Low-level State, Action and Intrinsic Reward Function Design

Low-level State Space. Previous works rely on the given waypoint data of the
vehicles and use PID controller to drive. In this paper, we consider learning
driving skills based on vision, which is similar to how human drive. Therefore
the low-level state space $s_{l}^{i}$ is defined as follows:

$\displaystyle
s^{i}_{l}=\left[s^{i}_{img},s^{i}_{speed},s^{i}_{laneID}\right]$

where $s^{i}_{img}$ denotes the captured image, $s^{i}_{speed}$ represents the
speed, and $s^{i}_{laneID}$ is the identifier of current lane.

Low-level Action Space. Because the low-level controller directly determines
the motion of the vehicle, the action space of the low-level controller is the
continuous linear speed and angular speed. If the high-level option selection
is lane-keeping, then the linear and angular speeds will remain the same
compared to the last step. Otherwise, the action spaces of different skills
are defined as follows:

$\displaystyle A^{i}_{\text{linear}}=\begin{cases}0.04:0.08&\text{ slow down
}\\\ 0.08:0.14&\text{ accelerate }\\\ 0.1:0.2&\text{ lane change }\end{cases}$
$\displaystyle A^{i}_{\text{angular}}=\begin{cases}-0.1:0.1&\text{ slow down
}\\\ -0.1:0.1&\text{ accelerate }\\\ 0.12:0.25&\text{ lane change
}\end{cases}$

where $A^{i}_{linear}$ denotes the action space for linear speed and
$A^{i}_{angular}$ is the action space for angular speed. To prevent unsafe
actions and improve exploration efficiency, we set different lower and upper
bounds of the linear speed and angular speed during the training.

Low-level Intrinsic Reward Functions. In contrast to the high-level team
reward which directly related to the success of the overall task, the
principle of intrinsic reward is to provide feedback for learning different
skills. Here, we categorize the intrinsic rewards into two types:
$r^{i}_{\text{driving in lane}}$ and $r^{i}_{\text{lane change}}$. An
illustration of reward calculation is shown in Fig. 4. In particular,
$r^{i}_{\text{driving in lane}}$ is used for $o^{i}_{h}=\left[\text{slow
down},\text{accelerate},\text{keep lane}\right]$ calculated as follows:

$\displaystyle r^{i}_{\text{driving in lane}}=\beta\cdot
r^{i}_{deviate}+(1-\beta)\cdot r^{i}_{travel}$

where $r^{i}_{deviate}$ denotes the distance deviated from the center of
current lane, $r^{i}_{travel}$ is the travel distance used to encourage the
movement of each vehicle, and $\beta$ is a hyper parameter controlling the
weights of deviated distance and travel distance.

Similarly, the intrinsic reward for lane change skill is designed as follows:

$\displaystyle r^{i}_{\text{lane change}}=\begin{cases}20&\text{ success }\\\
-20&\text{ fail }\\\ r^{i}_{travel}&\text{ else }\end{cases}$

For the lane change skill, we assign a positive reward of $20$ when the
vehicle successfully changes the lane and assign a negative penalty if the
vehicle fails to change lane within a predefined maximum time step.

Figure 6: Illustration of the cooperative lane change scenario, where the
vehicle $1$ should coordinate with vehicle $2$ to avoid the collision when
vehicle $2$ is performing lane change

##  V Performance Evaluation

This section presents the results of applying our approach in cooperative lane
change scenarios, which is shown in Fig. 6. We first create a simulation
environment for training and then build a real-world testbed for further
evaluation. We compare the performance of our approach with several well-known
reinforcement learning approaches.

TABLE I: Hyperparameters for Training Hyperparameter | Value
---|---
Training episode | $14,000$
Episode length | $30$
Buffer capacity | $100,000$
Batch size | $1024$
Learning rate | $0.01$
Discount factor $\gamma$ | $0.95$
Dimension of the hidden layer | $32$
Target network update rate | $0.01$
Figure 7: Comparison of the learning curve of different approaches in the
cooperative lane change scenarios.

###  V-A Baseline Approaches

In this paper, we apply different MARL approaches to cooperative lane change
scenarios, including distributed approaches (DTDE) and centralized training
with decentralized execution (CTDE) approaches. Specifically, the considered
methods for comparison include:

  * •

Independent Deep Q-learning (DQN). It is a distributed learning approach, in
which each agent trains a Q-network using its local observation and shared
team reward. Each agent applies the $\epsilon$-greedy strategy for action
exploration during the training.

  * •

Counterfactual policy gradient (COMA). It is a standard CTDE approach where
the centralized critic is trained with Q-learning. The actor network is
trained with the counterfactual policy gradient theorem.

  * •

Multi-agent mixed actor-critic (MADDPG). It is another approach adopting the
framework of CTDE. Unlike COMA, it considers the environment with explicit
communication and trains a centralized critic for each agent, allowing each
agent to have different reward functions.

  * •

Multi-Actor-Attention-Critic (MAAC). It is the state-of-the-art approach that
extends the attention mechanism, which is widely used in image processing and
natural language processing tasks, to multi-agent reinforcement learning. It
trains an actor-attention-critic network for each agent and allows parameter
sharing to improve the learning efficiency. Noted that, MAAC uses
decentralized critics with a decentralized actor with parameter sharing.

###  V-B Training Environment and Evaluation Metrics

To evaluate our approach, we first create a simulation environment based on
Gazebo, a popular physical simulation engine in robotics [25]. The simulation
environment includes the lane scenario and multiple vehicles equipped with
cameras and lidar sensors. Besides, we provide several RL-friendly APIs that
are convenient for algorithm implementation. Fig. 12 shows the simulation
scenarios of the cooperative lane change. When the vehicle in front of vehicle
$2$ stops or moves slowly, vehicle $2$ needs to control itself to perform lane
change which needs the coordination among vehicle $1$ and vehicle $2$. We use
a conventional neural network to encode the image data and a multi-layer
fully-connected neural network as the critic network for policy
implementation. The dimension of the hidden layer is set to $32$ in all the
algorithms.

At the beginning of each episode, we randomly initialize several positions of
the vehicles and set the same initial speed to them. Besides, the episode
length is set to $18$ time steps. The hyper-parameters of our network
architecture and learning algorithms are presented in Tab. I. When the
vehicles get collision, or the max episode length is reached, the program will
call the reset function provided by our testbed to reset the environment and
continue the training.

In this paper, we consider four evaluation metrics to measure the performance
of different approaches, including:

  * •

Mean episode reward: the average reward of each time step when sampling the
replay buffer from the different episodes.

  * •

Collision rate: the collision ratio of the vehicles in each episode.

  * •

Lane merge successful rate: the proportion of vehicles that successfully merge
in each episode.

  * •

Mean speed: the average of the vehicle’s speeds at each time step.

Figure 8: Episode reward of learning different low-level skills.

###  V-C Learning Low-level Skills

Figure 9: Real-world evaluation

Before training cooperative strategies, we first train the low-level skills,
including lane tracking and lane change. We create parallel training
environments with different intrinsic reward functions so that the low-level
learned skills can be further shared. To improve learning efficiency, we set a
lower angular speed in training the lane-keeping skill and set a higher
angular speed when training the lane change skill. Fig. 8 shows the episode
rewards during the training process. The experimental results imply that the
soft actor-critic algorithm can successfully converge in training the lane
tracking and lane change tasks. To be noticed, the episode reward of learning
lane change policy remains a low value before $5,000$ episodes. The reason is
that the agent will explore the action space at the beginning to maximize the
entropy of action probability. The exploration strategy makes the low-level
skill model robust.

###  V-D Learning High-level Cooperative Policy

Figure 10: Learning loss of modeling other partners’ policies. Figure 11:
Learning high-level cooperative policy in the simulation environment we
created. Figure 12: Training in the Gazebo simulator.

After training the low-level skills in the single-vehicle environment, we
gradually increase the number of vehicles to learn the coordination policy. As
shown in Fig. 9, we set up four vehicles in a double-lane track scenario.
Vehicle $4$ is set with a plodding speed to simulate traffic congestion or
traffic accident. The other vehicles are initialized with an average speed and
random positions. Once the collision happens during lane-changing or lane-
keeping, each vehicle will receive a negative reward of $-20$, and the episode
will end and restart.

During the training, each vehicle will train opponent models to predict
others’ actions and use the opponent model to stabilize the offline training
of the Q-value network. Fig. 10 shows the loss of modeling other vehicles’
policies from vehicle 2’s perspective. As we can see, the vehicle $1$’s action
prediction model become converged at a fast speed while the vehicle $3$’
action prediction model become converged after 12000 episodes. The difference
in convergence speed also illustrates the different interactions among
vehicles.

Figure 13: We develop two prototypes of vehicles in the real-world testbed. In this experiment, we use the first prototype in the testing. TABLE II: Performance evaluation in real-world platform Method Metrics | Collision Rate |  Successful Rate |  Mean Speed
---|---|---|---
COMA | $0.35$ | $0.65$ | $0.06344$
Independent DQN | $1.0$ | $0.0$ | $0.05395$
MAAC | $0.25$ | $0.65$ | $0.0625$
MADDPG | $0.95$ | $0.5$ | $0.07029$
Ours | $0.2$ | $0.8$ | $0.072$

Fig. 7 shows the performance of different RL approaches in terms of mean
episode reward, collision rate and lane change successful rate. As shown in
Fig. 7(a), our method achieved the highest episode reward than other
baselines. Besides, the lower bound of our method is highest than others which
also implicitly shows the stability of our approach. Fig. 7(b) and Fig. 7(c)
shows the collision rate and lane change successful rate among different
approaches. Almost all the RL approaches can reduce the collision rate after
14000 training episodes except MADDPG. We found that MADDPG has a lower
learning efficiency as it still maintains a higher collision rate than others.
Both DQN and our method reach the lowest collision rate. However, the lane
change success rate is near $0$. We found that vehicle $2$ learn a policy to
drive slowly instead of changing the lane during the training. Instead, the
vehicle trained with our method can successfully perform cooperative lane
change without collision with others. Moreover, Fig. 11 shows that our method
achieves the highest mean speed, which is $0.08$, and the vehicles trained
with MAAC received the lowest speed of $0.048$. In general, our method shows
advantages in terms of safety and efficiency.

###  V-E Real-world Experiments

To investigate the gap between simulation and reality, we deploy the learned
policies in simulation to a real-world testbed, which consists of a two-lane
track and multiple vehicles shown in Fig. 13. Each vehicle is equipped with a
camera, lidar, and edge server. We run $20$ episodes for each MRRL method with
several random initial positions. Besides, we run a master node on the server
to monitor the state of each vehicle and calculate the collision rate, lane-
merging successful rate, and mean speed.

As shown in Tab. II, our approach reaches a low collision rate and high speed
than others, which is $0.2$ and $0.07$, respectively. Besides, the collision
rate of MAAC is lower than other baselines such as COMA and MADDPG. However,
the performance of Independent DQN is different in simulation. We argue that
the diversity and robustness of DQN are poor than other policy-based
approaches.

##  VI Related Work

This section summarizes the previous works related to multi-agent
reinforcement learning and hierarchical reinforcement learning in Sec. VI-A
and Sec. VI-B, respectively.

###  VI-A Multi-agent Reinforcement Learning

Recent works on reinforcement learning (RL) have shifted from single-agent
reinforcement learning to multi-agent reinforcement learning (MARL). MARL
corresponds to the learning problem in multi-agent systems in which multiple
agents learn simultaneously in a shared environment. Existing MARL approaches
can be classified into centralized RL, centralized training with decentralized
execution (CTDE), and distributed RL. Centralized RL trains a centralized
network mapping from the state space to all the agents’ joint action space.
Such an approach suffers from the high dimension issue of joint state and
action spaces, making it hard to be extended to large-scale scenarios.

Some previous works adopt the paradigm of centralized training with
decentralized execution (CTDE). Lowe et al. proposed a multi-agent actor-
critic (MADDPG) approach and applied it for cooperative navigation [26].
Foerster et al. proposed counterfactual policy gradient (COMA), using a
counterfactual baseline to address the credit assignment problem in MARL [27].
Similar works of CTDE can also be found in flocking [8] and pathfinding [9].
However, the number of features of the centralized critic network needs to be
scaled up linearly (in the best case) or exponentially (in the worse case)
with the increases of the number of agents [10]. Besides, different agents may
have different influences, and decentralized critics are expected to perform
better in these scenarios.

Recent developments in DRL can be further divided into two types: independent
RL and decentralized critics with decentralized actors. Independent RL is
self-interested and does not consider the states and actions of other agents
[28]. Independent RL has good scalability but poor cooperation performance. An
alternative method is to use decentralized critics with decentralized
actors(DTDE). For example, Zhang et al. proposed a general decentralized MARL
framework [11] and Iqbal et al. proposed an actor-attention-critic approach
for multi-agent cooperation [12]. However, these methods either rely on
explicit information sharing or use end-to-end RL models, which are not
adaptive to complex and complicated tasks. Moreover, these end-to-end RL
models lack interpretability on the interaction among agents.

###  VI-B Hierarchical Reinforcement Learning

Hierarchical reinforcement learning (HRL) is a type of reinforcement learning
(RL) that leverages the hierarchical structure of a given task learns a
hierarchical policy [29][30]. Dayan et al. developed the first HRL work and
proposed feudal reinforcement learning (FRL), where the high-level managers
learn how to set the tasks to the sub-managers while the sub-managers learn
how to perform the sub-tasks [31]. In this way, a complex task can be divided
into small sub-tasks that are easier to be solved [31]. Kulkarni et al.
introduced intrinsic reward design for low-level tasks and used a deep neural
network to approximate the value function [22]. The proposed intrinsic reward
mechanism motivates the agents to explore new behavior for its own sake.

Some previous works have applied hierarchical reinforcement learning for
multi-agent cooperation. For example, Dietterich proposed the MAXQ framework
to learn the hierarchical policy of each agent [32]. Besides, Ghavamzadeh et
al. proposed a COM-Cooperative HRL for multi-agent communication [33].
Recently, Ahilan et al. proposed a deep MAHRL approach, in which a central
manager policy chooses subtasks for other workers simultaneously [34]. Similar
deep MAHRL approaches can also be found in [35]. However, all of these works
either only consider the discrete action space of each agent or rely on
centralized training.

##  VII Conclusion and Future Directions

This paper studies how distributed multiple agents learn to cooperate in
continuous action space and the non-stationarity issue. Unlike the traditional
end-to-end RL methods, we propose a hierarchical reinforcement learning
approach for distributed multi-agent cooperation. The agent cooperation is
efficiently learned in high-level discrete action space, while the low-level
individual control is handled by independent reinforcement learning. We
incorporate opponent modeling in the high-level layer to encourage cooperation
and stabilize Q-learning. Then, we present a case study of cooperative lane
change and conduct extensive experiments via simulation and a real-world
testbed. In the future, we will investigate deeper hierarchy discovery and the
theoretical guarantee. Currently, the design of task decomposition and state-
action space in different layers still requires human knowledge. It is still
an open question how to discover the hierarchical architecture automatically.
Besides, most existing MARL approaches are trained and evaluated in
simulation. The gap between simulation and reality also deserves more
investigation.

##  VIII Acknowledgment

The research is supported by Hong Kong RGC TRS (project T41-603), RIF
(projects R5009-21 and R5060-19), CRF (projects C5026-18G and C5018-20GF), and
GRF (project 15204921).

## References

  * [1] M. A. L. Silva, S. R. de Souza, M. J. F. Souza, and A. L. C. Bazzan, “A reinforcement learning-based multi-agent framework applied for solving routing and scheduling problems,” _Expert Systems with Applications_ , vol. 131, pp. 148–171, 2019.
  * [2] P. Zhao, S. Suryanarayanan, and M. G. Simoes, “An energy management system for building structures using a multi-agent decision-making control methodology,” _IEEE Transactions on Industry Applications_ , vol. 49, no. 1, pp. 322–330, 2012.
  * [3] S. Jiang, J. Liang, J. Cao, and R. Liu, “An ensemble-level programming model with real-time support for multi-robot systems,” in _IEEE International Conference on Pervasive Computing and Communication Workshops (PerCom Workshops)_ , 2016, pp. 1–3.
  * [4] L. Ma, Z. Wang, and H.-K. Lam, “Event-triggered mean-square consensus control for time-varying stochastic multi-agent system with sensor saturations,” _IEEE Transactions on Automatic Control_ , vol. 62, no. 7, pp. 3524–3531, 2016.
  * [5] Y. Sahni, J. Cao, and S. Jiang, “Middleware for multi-robot systems,” in _Mission-oriented Sensor Networks and Systems: Art and Science_. Springer, 2019, pp. 633–673.
  * [6] J. Wang, J. Cao, M. Stojmenovic, M. Zhao, J. Chen, and S. Jiang, “Pattern-rl: Multi-robot cooperative pattern formation via deep reinforcement learning,” in _18th IEEE International Conference on Machine Learning And Applications (ICMLA)_ , 2019, pp. 210–215.
  * [7] P. Hernandez-Leal, B. Kartal, and M. E. Taylor, “A survey and critique of multiagent deep reinforcement learning,” _Autonomous Agents and Multi-Agent Systems_ , vol. 33, no. 6, pp. 750–797, 2019.
  * [8] P. Zhu, W. Dai, W. Yao, J. Ma, Z. Zeng, and H. Lu, “Multi-robot flocking control based on deep reinforcement learning,” _IEEE Access_ , vol. 8, pp. 150 397–150 406, 2020.
  * [9] M. Damani, Z. Luo, E. Wenzel, and G. Sartoretti, “Primal $\\_2$: Pathfinding via reinforcement and imitation multi-agent learning-lifelong,” _IEEE Robotics and Automation Letters_ , vol. 6, no. 2, pp. 2666–2673, 2021.
  * [10] X. Lyu, Y. Xiao, B. Daley, and C. Amato, “Contrasting centralized and decentralized critics in multi-agent reinforcement learning,” in _International Conference on Autonomous Agents and Multiagent Systems (AAMAS)_. ACM, 2021, pp. 844–852.
  * [11] K. Zhang, Z. Yang, H. Liu, T. Zhang, and T. Basar, “Fully decentralized multi-agent reinforcement learning with networked agents,” in _International Conference on Machine Learning (ICML)_ , 2018, pp. 5872–5881.
  * [12] S. Iqbal and F. Sha, “Actor-attention-critic for multi-agent reinforcement learning,” in _International Conference on Machine Learning (ICML)_ , 2019, pp. 2961–2970.
  * [13] C. Qu, S. Mannor, H. Xu, Y. Qi, L. Song, and J. Xiong, “Value propagation for decentralized networked deep multi-agent reinforcement learning,” in _Advances in Neural Information Processing Systems (NeurIPS)_ , 2019, pp. 1182–1191.
  * [14] T. Schaul, J. Quan, I. Antonoglou, and D. Silver, “Prioritized experience replay,” in _International Conference on Learning Representations (ICLR) Posters_ , 2016.
  * [15] R. S. Sutton, D. A. McAllester, S. P. Singh, and Y. Mansour, “Policy gradient methods for reinforcement learning with function approximation,” in _Advances in Neural Information Processing Systems (NeurIPS)_ , 2000, pp. 1057–1063.
  * [16] L. Busoniu, R. Babuska, and B. De Schutter, “A comprehensive survey of multiagent reinforcement learning,” _IEEE Transactions on Systems, Man, and Cybernetics, Part C (Applications and Reviews)_ , vol. 38, no. 2, pp. 156–172, 2008.
  * [17] V. Mnih, K. Kavukcuoglu, D. Silver, A. A. Rusu, J. Veness, M. G. Bellemare, A. Graves, M. Riedmiller, A. K. Fidjeland, G. Ostrovski _et al._ , “Human-level control through deep reinforcement learning,” _Nature_ , vol. 518, no. 7540, pp. 529–533, 2015.
  * [18] D. Silver, G. Lever, N. Heess, T. Degris, D. Wierstra, and M. Riedmiller, “Deterministic policy gradient algorithms,” in _International Conference on Machine Learning (ICML)_ , 2014, pp. 387–395.
  * [19] T. P. Lillicrap, J. J. Hunt, A. Pritzel, N. Heess, T. Erez, Y. Tassa, D. Silver, and D. Wierstra, “Continuous control with deep reinforcement learning,” in _International Conference on Learning Representations (ICLR) Posters_ , 2016.
  * [20] S. Pateria, B. Subagdja, A.-h. Tan, and C. Quek, “Hierarchical reinforcement learning: A comprehensive survey,” _ACM Computing Surveys_ , vol. 54, no. 5, pp. 1–35, 2021.
  * [21] R. S. Sutton, D. Precup, and S. Singh, “Between mdps and semi-mdps: A framework for temporal abstraction in reinforcement learning,” _Artificial Intelligence_ , vol. 112, no. 1-2, pp. 181–211, 1999.
  * [22] T. D. Kulkarni, K. Narasimhan, A. Saeedi, and J. Tenenbaum, “Hierarchical deep reinforcement learning: Integrating temporal abstraction and intrinsic motivation,” in _Advances in Neural Information Processing Systems (NeurIPS)_ , vol. 29, 2016, pp. 3675–3683.
  * [23] M. L. Littman, “Markov games as a framework for multi-agent reinforcement learning,” in _Machine Learning Proceedings (ICML)_. Elsevier, 1994, pp. 157–163.
  * [24] T. Haarnoja, A. Zhou, P. Abbeel, and S. Levine, “Soft actor-critic: Off-policy maximum entropy deep reinforcement learning with a stochastic actor,” in _International Conference on Machine Learning (ICML)_ , 2018, pp. 1861–1870.
  * [25] N. Koenig and A. Howard, “Design and use paradigms for gazebo, an open-source multi-robot simulator,” in _IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , vol. 3, 2004, pp. 2149–2154.
  * [26] R. Lowe, Y. Wu, A. Tamar, J. Harb, P. Abbeel, and I. Mordatch, “Multi-agent actor-critic for mixed cooperative-competitive environments,” in _Advances in Neural Information Processing Systems (NeurIPS)_ , 2017, pp. 6379–6390.
  * [27] J. Foerster, G. Farquhar, T. Afouras, N. Nardelli, and S. Whiteson, “Counterfactual multi-agent policy gradients,” in _AAAI Conference on Artificial Intelligence (AAAI)_ , 2018, pp. 2974–2982.
  * [28] K. Sivanathan, B. Vinayagam, T. Samak, and C. Samak, “Decentralized motion planning for multi-robot navigation using deep reinforcement learning,” in _International Conference on Intelligent Sustainable Systems (ICISS)_. IEEE, 2020, pp. 709–716.
  * [29] T. G. Dietterich, “The maxq method for hierarchical reinforcement learning,” in _International Conference on Machine Learning (ICML)_ , vol. 98. Citeseer, 1998, pp. 118–126.
  * [30] N. K. Jong and P. Stone, “Hierarchical model-based reinforcement learning: R-max+ maxq,” in _International Conference on Machine Learning (ICML)_ , 2008, pp. 432–439.
  * [31] P. Dayan and G. E. Hinton, “Feudal reinforcement learning,” in _Advances in Neural Information Processing Systems (NeurIPS)_. Morgan Kaufmann, 1992, pp. 271–278.
  * [32] T. G. Dietterich, “Hierarchical reinforcement learning with the maxq value function decomposition,” _Journal of Artificial Intelligence Research_ , vol. 13, pp. 227–303, 2000.
  * [33] M. Ghavamzadeh and S. Mahadevan, “Learning to communicate and act using hierarchical reinforcement learning,” in _International Joint Conference on Autonomous Agents and Multiagent Systems (AAMAS)_ , 2004, pp. 1114–1121.
  * [34] S. Ahilan and P. Dayan, “Feudal multi-agent hierarchies for cooperative reinforcement learning,” in _Multidisciplinary Conference on Reinforcement Learning and Decision Making (RLDM)_ , 2019, p. 57.
  * [35] J. Chakravorty, P. N. Ward, J. Roy, M. Chevalier-Boisvert, S. Basu, A. Lupu, and D. Precup, “Option-critic in cooperative multi-agent systems,” in _International Conference on Autonomous Agents and Multiagent Systems (AAMAS)_ , 2020, pp. 1792–1794.