11institutetext:  Beijing National Laboratory for Condensed Matter Physics and
Institute of Physics, Chinese Academy of Sciences, Beijing 100190, China
School of Physical Sciences, University of Chinese Academy of Sciences,
Beijing 100190, China
School of Physical Science and Technology, ShanghaiTech University, Shanghai
201210, China
Department of Physics and Beijing Key Laboratory of Opto-electronic Functional
Materials & Micro-nano Devices, Renmin University of China, Beijing 100872,
China
Songshan Lake Materials Laboratory, Dongguan, Guangdong 523808, China
School of Physical Science and Technology, Lanzhou University, Lanzhou 730000,
China
Kavli Institute for Theoretical Sciences, Beijing 100190, China

Crystallographic aspects of phase transformations Conductivity phenomena in
semiconductors and insulators Occurrence, potential candidates

# First order transition in trigonal structure
${\textbf{Ca}}{\textbf{Mn}}_{2}{\textbf{P}}_{2}$

Y. J. Li 1122 F. Jin 1122 Z. Y. Mi 1122 J. Guo 1122 W. Wu 1122 Z. H. Yu 33 D.
S. Wu 1122 S. H. Na 1122 C. Mu 1122 X. B. Zhou 11 Z. Li 1122 K. Liu 44 L. L.
Sun 112255 Q. M. Zhang 6611 T. Xiang 112277 G. Li 112255 J. L. Luo E-mail:
112255<EMAIL_ADDRESS>(corresponding author) 11223344556677

###### Abstract

We report structural and physical properties of the single crystalline
${\mathrm{Ca}}{\mathrm{Mn}}_{2}{\mathrm{P}}_{2}$. The X-ray diffraction(XRD)
results show that ${\mathrm{Ca}}{\mathrm{Mn}}_{2}{\mathrm{P}}_{2}$ adopts the
trigonal ${\mathrm{Ca}}{\mathrm{Al}}_{2}{\mathrm{Si}}_{2}$-type structure.
Temperature dependent electrical resistivity $\rho(T)$ measurements indicate
an insulating ground state for
${\mathrm{Ca}}{\mathrm{Mn}}_{2}{\mathrm{P}}_{2}$ with activation energies of
40 meV and 0.64 meV for two distinct regions, respectively. Magnetization
measurements show no apparent magnetic phase transition under 400 K. Different
from other ${\mathrm{A}}{\mathrm{Mn}}_{2}{\mathrm{Pn}}_{2}$ (A = Ca, Sr, and
Ba, and Pn = P, As, and Sb) compounds with the same structure, heat capacity
$C_{\mathrm{p}}(T)$ and $\rho(T)$ reveal that
${\mathrm{Ca}}{\mathrm{Mn}}_{2}{\mathrm{P}}_{2}$ has a first-order transition
at $T$ = 69.5 K and the transition temperature shifts to high temperature upon
increasing pressure. The emergence of plenty of new Raman modes below the
transition, clearly suggests a change in symmetry accompanying the transition.
The combination of the structural, transport, thermal and magnetic
measurements, points to an unusual origin of the transition.

###### pacs:

61.50.Ks

###### pacs:

72.20.-i

###### pacs:

74.10.+v

##  1 Introduction

Exploring new superconductors in the 3d transition-metal-based compounds is
one of the most active research topics in modern condensed matter physics.
Owning to the interplay among charge, spin, and orbital degrees of freedom,
various ground states could be stabilized in these materials.
Superconductivity has been observed in the majority of 3d transition-metal-
based compounds, however it is still rare for Cr- and Mn-based ones. Up to
now, a few Cr-based superconductors has been reported at either ambient or
elevated pressure, such as CrAs [1], quasi one-dimensional
${\mathrm{A}}_{2}{\mathrm{Cr}}_{3}{\mathrm{As}}_{3}$,
${\mathrm{A}}{\mathrm{Cr}}_{3}{\mathrm{As}}_{3}$ (A=Na, K, Rb, Cs) and
${\mathrm{Ln}}_{3}{\mathrm{Cr}}_{10-x}{\mathrm{N}}_{11}$ (Ln=La, Pr)[2, 3, 4,
5, 6, 7, 8]. Meanwhile, to our knowledge, there is only one Mn-based
superconductor, MnP [9]. It is of great interest to explore other possible Mn-
based superconductors[10].

The body-centered tetragonal ternary compounds
${\mathrm{A}}{\mathrm{Fe}}_{2}{\mathrm{As}}_{2}$ (A=Ca, Sr, Ba, Eu) are
metallic and show nearly contiguous antiferromagnetic spin-density wave(SDW)
and structure transitions [11, 12, 13, 14, 15, 16, 17, 18]. The suppression of
these transitions by external pressure or chemical doping leads to
superconductivity with bulk transition temperature, $T_{c}$ up to 45 K [19,
20, 21, 22]. The iron atoms can be fully replaced by other transition metal
atoms, such as Cr, Mn, Co, Ni, and Cu [23, 24, 25, 26, 27]. Among them,
${\mathrm{Sr}}{\mathrm{Ni}}_{2}{\mathrm{As}}_{2}$ and
${\mathrm{Ba}}{\mathrm{Ni}}_{2}{\mathrm{As}}_{2}$ were found to be
superconductors with $T_{c}$ = 0.62 and 0.7 K, respectively [28, 29].

The Mn pnictides ${\mathrm{A}}{\mathrm{Mn}}_{2}{\mathrm{Pn}}_{2}$ (A=alkaline
earth, Pn=P, As, Sb) have been known for decade. Different from
${\mathrm{Ba}}{\mathrm{Mn}}_{2}{\mathrm{Pn}}_{2}$ crystallizing in the
tetragonal ${\mathrm{Th}}{\mathrm{Cr}}_{2}{\mathrm{Si}}_{2}$ structure type
($I4/mmm$) similar to ${\mathrm{A}}{\mathrm{Fe}}_{2}{\mathrm{As}}_{2}$, the
compounds ${\mathrm{Ca}}{\mathrm{Mn}}_{2}{\mathrm{As}}_{2}$,
${\mathrm{Ca}}{\mathrm{Mn}}_{2}{\mathrm{Sb}}_{2}$, and
${\mathrm{Sr}}{\mathrm{Mn}}_{2}{\mathrm{Pn}}_{2}$ crystalize in the trigonal
${\mathrm{Ca}}{\mathrm{Al}}_{2}{\mathrm{Si}}_{2}$ structure ($P\bar{3}m1$).
Previous reports suggest that they are all insulating at low temperature and
have antiferromagnetic (AFM) order [24, 30, 31, 32, 33, 34, 35, 36]. However,
for ${\mathrm{Ca}}{\mathrm{Mn}}_{2}{\mathrm{P}}_{2}$, apart from the
structural data, no detailed physical properties have been reported so far,
which gives us motivation to carry out current study.

In this work, we report the growth, crystal structure, electrical resistivity
$\rho(T)$ in the $ab$ plane at ambient and elevated pressure, Raman spectra,
magnetic susceptibility $\chi(\mathrm{T})$, and heat capacity
$C_{\mathrm{p}}(T)$ measurements of
${\mathrm{Ca}}{\mathrm{Mn}}_{2}{\mathrm{P}}_{2}$ single crystals. It is found
that ${\mathrm{Ca}}{\mathrm{Mn}}_{2}{\mathrm{P}}_{2}$ has a structural type
first order transition at 69.5 K, at the same time there is no long range
magnetic order developed below 400 K, thus it is different from other
isostructural Mn- based pnictides and calls for further study.

##  2 Experiments

Single crystal of ${\mathrm{Ca}}{\mathrm{Mn}}_{2}{\mathrm{P}}_{2}$ was grown
using Sn flux. High purity elements Ca (99%), Mn (99.99%), P (99.999%), and Sn
(99.9999%) were weighed in the molar ratio Ca:Mn:P:Sn =1.05:2:2:20 and placed
in an alumina crucible. The crucible was then sealed into a quartz tube under
an Ar pressure of 1/4 ATM. After preheating the mixture at
300${}^{\circ}\mathrm{C}$ for 24h, the assembly was heated to
1100${}^{\circ}\mathrm{C}$, dwelled at the temperature for 24 hours for
homogenization, and subsequently cooled down to 700${}^{\circ}\mathrm{C}$ at
the rate of 2${}^{\circ}\mathrm{C}$/h. Finally the tube was taken out of the
furnace quickly to separate the tin flux by centrifuging. Shiny, plate-like
single crystals in hexagon shape with maximum dimensions 3$\times$ 1.5$\times$
1$\mathrm{mm}^{3}$ were obtained.

Powder XRD experiments were conducted by a PAN-analytical X-ray diffractometer
with Cu $\mathrm{K}_{\alpha}$ radiation, the angle 2$\theta$ range was from
10∘ to 80∘. The element analysis was performed by using an energy dispersive
X-ray spectroscopy (EDX) equipped on a Hitachi S-4800 scanning electron
microscope (SEM). Single crystal XRD structural analysis of
${\mathrm{Ca}}{\mathrm{Mn}}_{2}{\mathrm{P}}_{2}$ was performed at room
temperature and 40 K in a Bruker D8 Venture diffractometer using Mo
$\mathrm{K}_{\alpha}$ radiation, $\lambda$ =0.71073 Å. Structure refinement
was performed by the program SHELXL-2014/7 [37] embedded in the program suite
APEX3.

The electrical resistivity and specific heat were measured in a Quantum Design
physical property measurement system (PPMS-9) by the standard four-probe
method with electrodes along the $ab$ plane and relaxation method,
respectively. High-pressure resistance measurements were performed in a
diamond-anvil cell (DAC), in which diamond anvils with 300 $\mu m$ flats and a
nonmagnetic rhenium gasket with 100 $\mu m$ diameter hole were employed. NaCl
powder was used as the pressure medium to provide a quasi-hydrostatic pressure
environment. The temperature dependence of DC magnetic susceptibility was
measured in a Quantum Design magnetic property measurement system (MPMS).

The polarized Raman spectra were collected with a Jobin Yvon LabRam HR
Evolution spectrometer equipped with a volume Bragg grating low-wave-number
suite, a liquid-nitrogen-cooled back-illuminated charge-coupled device
detector, and a 532-nm laser. The laser was focused into a spot of $\sim$ 5
$\mu m$ in diameter on the sample surface ($ab$ plane), with a power $<$100
$\mu$W, to avoid overheating.

##  3 Results And Discussions

Powder XRD data on crushed ${\mathrm{Ca}}{\mathrm{Mn}}_{2}{\mathrm{P}}_{2}$
single crystals at room temperature is shown in fig. 1(a)

[scale=0.33]Fig1.eps

Figure 1: Structural and elemental analysis of
${\mathrm{Ca}}{\mathrm{Mn}}_{2}{\mathrm{P}}_{2}$. (a) Powder XRD spectrum of
crushed crystals at room temperature with indexing. (b) EDX spectrum. (c) An
illustration of the crystal structure. (d) Single crystal XRD spectrum at room
temperature, inset shows a photo of grown crystal on milimeter scale.

. The reflections can be well indexed by the trigonal unit cell. However, weak
peaks from adventitious Sn flux are also visible and marked by stars. The
reason is that due to a pile of single crystals were used to prepare the
sample for powder XRD, the Sn residual on the surface of
${\mathrm{Ca}}{\mathrm{Mn}}_{2}{\mathrm{P}}_{2}$ single crystals were not
completely removed and entered XRD spectrum. After further cleaning, EDX
analysis was performed to check the chemical composition of the grown single
crystals. The typical EDX spectrum of an individual crystal is shown in fig.
1(b). Only three elements, Ca, Mn and P are detected. The average ratio of the
elements at different locations in the crystals is 19.79:40.12:40.09, which is
close to the 1:2:2 stoichiometry of the compound. The Sn flux and any other
possible impurity were not detected in present analysis within instrumental
sensitivity.

A tiny piece of single crystal was used for the structure refinement at room
temperature. R1 and wR2 are residual factors. The final anisotropic full-
matrix least-squares refinement on $F^{2}$ (F is structural factor) with 10
variables converges at R1 for the observed data and wR2 for all data[37]. The
result confirms that the crystal has the trigonal
${\mathrm{Ca}}{\mathrm{Al}}_{2}{\mathrm{Si}}_{2}$-type structure with space
group $P\bar{3}m1$, and the refined lattice constants are a=4.0973(1) Å,
c=6.8473(3) Å, respectively. The refinement parameters are tabulated in table
1,

Table 1: Refined crystallographic parameters of
${\mathrm{Ca}}{\mathrm{Mn}}_{2}{\mathrm{P}}_{2}$ from single crystal X-ray
structure analysis at 293K and 40K.

Refined parameter | 293K | 40K
---|---|---
Wavelength | 0.71073Å | 0.71073Å
Structure |  ${\mathrm{Ca}}{\mathrm{Al}}_{2}{\mathrm{Si}}_{2}$-type |  ${\mathrm{Ca}}{\mathrm{Al}}_{2}{\mathrm{Si}}_{2}$-type
| trigonal | trigonal
Space group |  $P\bar{3}m1$ (No.164) |  $P\bar{3}m1$ (No.164)
Lattice parameters | a=4.0973(1) Å | a=4.0887(1) Å
| b=4.0973(1) Å | b=4.0887(1) Å
| c=6.8473(3) Å | c=6.8207(3) Å
|  $\alpha$=90∘ |  $\alpha$=90∘
|  $\beta$=90∘ |  $\beta$=90∘
|  $\gamma$=120∘ |  $\gamma$=120∘
Absorption coefficient |  8.187 mm-1 |  8.253 mm-1
R1 | 2.37% | 1.4%
wR2 | 5.7% | 3.3%
The goodness-of fit | 1.412 | 1.167
Ca1 | (0,0,0) | (0,0,0)
Mn1 | (0.667,0.333,0.625) | (0.667,0.333,0.624)
P1 | (0.667,0.333,0.261) | (0.667,0.333,0.262)

which are in accordance to the early report [38].

The crystal structure of trigonal
${\mathrm{Ca}}{\mathrm{Mn}}_{2}{\mathrm{P}}_{2}$ is shown in fig. 1(c).
Similar to other 122-type Mn pnictide compounds, it consists of corrugated
honeycomb [${\mathrm{Mn}}_{2}{\mathrm{P}}_{2}$]-2 layers that are stacked
along the $c$ axis and separated by Ca+2 cations. Alternatively, the Mn
sublattice can be viewed as triangular double layers of Mn stacked along the
$c$ axis and separated by Ca atoms. The XRD results of
${\mathrm{Ca}}{\mathrm{Mn}}_{2}{\mathrm{P}}_{2}$ single crystal are shown in
fig. 1(d). Only four main peaks of (001), (002), (003) and (005) are clearly
visible, whose positions are in good agreement with those in powder XRD data,
indicating that the crystal is grown along the c-axis direction.

Figure 2(a)

[scale=0.6]Fig2.eps

Figure 2: (a) Temperature dependence of in-plane electrical resistivity
$\rho(T)$ of ${\mathrm{Ca}}{\mathrm{Mn}}_{2}{\mathrm{P}}_{2}$ from 30 K to 400
K. The inset shows the hysteresis resistivity for
${\mathrm{Ca}}{\mathrm{Mn}}_{2}{\mathrm{P}}_{2}$ around 69.5 K. (b) The plot
of $\log_{10}\rho$ versus 1000/T. The red dash curves through the data are
linear fittings(See text).

gives the $ab$ plane $\rho(T)$ from 30 to 400 K. It is clear that
${\mathrm{Ca}}{\mathrm{Mn}}_{2}{\mathrm{P}}_{2}$ has an insulating ground
state, on the contrary to the prediction by a semi-local DFT band structure
calculation [39], in which neither magnetic exchange nor on-site Coulomb
repulsion were considered. Whether the insulating gap is due to electron
correlation as in a Mott insulator requires further study. The plot of
$\log_{10}\rho$ versus 1000/T is shown in fig. 2 (b). The $\rho(T)$ data over
confined temperature interval is fitted by the expression

$\log_{10}\rho=\mathrm{A}+2.303\Delta/\mathrm{k_{B}T},$ (1)

where A is a constant, $k_{B}$ is Boltzmann,s constant, and $\Delta$ is the
activation energy, respectively. The data between 250 K and 400 K and between
30 K and 60 K are nearly linear in $T$ and were separately fitted by Eq. 1 ,
shown as the dash curves. The high temperature range fit yields an activation
energy $\Delta_{1}$=40 meV, which is of the same order as previously found for
${\mathrm{Ca}}{\mathrm{Mn}}_{2}{\mathrm{As}}_{2}$ (61 meV) and
${\mathrm{Sr}}{\mathrm{Mn}}_{2}{\mathrm{P}}_{2}$ (12.9 meV)[24, 32], and is
the intrinsic activation energy of
${\mathrm{Ca}}{\mathrm{Mn}}_{2}{\mathrm{P}}_{2}$. The low temperature range
fit gives a $\Delta_{2}$= 0.64 meV, which should be extrinsic and originated
from the energy gap between donor (acceptor) energy levels and the conduction
(valence) band [40].

As a feature that has not been observed on other 122-type Mn pnictides, there
is a drop of resistivity at $T$ = 69.5 K, which is enlarged as the inset of
fig. 2(a), the thermal hysteresis is confirmed by repeated measurements with
different cooling/warming speeds, suggesting it corresponds to a first-order
transition.

The existence of a first-order transition at 69.5 K is further revealed by
heat capacity measurement. In the conventional relaxation method, the
temperature versus time curves were recorded as the raw data. Shown as the
lower inset of fig. 3

[scale=0.32]Fig3.eps

Figure 3: Heat capacity $C_{\mathrm{p}}(T)$ of
${\mathrm{Ca}}{\mathrm{Mn}}_{2}{\mathrm{P}}_{2}$ at both 0 T and 9 T, combined
with a Debye lattice contribution simulation. The upper inset shows the
enlarged region around the first order transition. The lower inset shows a raw
temperature relaxation curve demonstrating the existence of latent heat.

, on a typical fitting curve with the starting temperature is slightly lower
than the transition temperature, latent heat is clearly observed on both the
raise and the fall processes. Therefore around that temperature a slope
analysis is employed to extract $C_{p}$ value, resulting in the sharp peak in
the $C_{p}$ versus $T$ curve of fig. 3. By applying 9 T magnetic field
perpendicular to the crystalline $ab$ plane, neither the position nor the
amplitude of the peak is changing, suggesting that the transition is not
affected by magnetic fields up to 9 T. Considering the non-metallic ground
state previously found from the $\rho(T)$ behavior and the absence of magnetic
contribution, the lattice heat capacity alone could account for the total heat
capacity at low temperature. We model the $C_{\mathrm{p}}(T)$ data from 2 K to
60 K by

$C_{p}=C_{lattice}=nC_{Debye}$ (2)

in which the Debye lattice heat capacity per mole is

$C_{Debye}=9R\left(\frac{T}{\Theta_{D}}\right)^{3}\int_{0}^{\Theta_{D}/T}\frac{x^{4}e^{x}}{\left(e^{x}-1\right)^{2}}dx$
(3)

where R=8.314 J/mol K is the molar gas constant, and n is the number of atoms
per formula unit (n=5 for ${\mathrm{Ca}}{\mathrm{Mn}}_{2}{\mathrm{P}}_{2}$)
[41]. The fitting yields Debye temperature $\Theta_{D}$ = 345(1) K and the
calculated lattice heat capacity extended up to 200 K is shown as a solid red
line superimposed to experimental data in fig. 3.

In order to better understand the first-order transition of
${\mathrm{Ca}}{\mathrm{Mn}}_{2}{\mathrm{P}}_{2}$, Raman measurements at both
300 K and 10 K are performed on a well characterized crystal, as shown in fig.
4.

[scale=0.9]Fig4.eps

Figure 4: Raman spectra of ${\mathrm{Ca}}{\mathrm{Mn}}_{2}{\mathrm{P}}_{2}$ at
(a) 300 K (b) 10 K with different polarization.

Four Raman peaks are observed at 300 K, all of which could be assigned
according to the symmetry of the crystal structure. However, another 8 new
Raman peaks (marked with red arrows) appeared at low temperature. Meanwhile
all those high temperature peaks remain present with tiny shift in energy,
thus by decreasing temperature through the 69.5 K first order transition, the
lattice symmetry is not fully altered, but possibly formed super-structure.

Single crystal XRD structural analysis and powder XRD at temperatures below
the first-order transition were also carried out. The refined crystal
parameters obtained at 40 K are listed in table 1. From 293 K to 40 K, the
lattice constants shrink 0.21% and 0.39% for the a- and c-axis, respectively.
Meanwhile the lattice space group remains the same. We did not observe any
distinguishable peak-splitting nor new peaks at 40 K on the XRD spectrum.
Considering the instrumental sensitivity, those results provide further
evidence that ${\mathrm{Ca}}{\mathrm{Mn}}_{2}{\mathrm{P}}_{2}$ forms super-
structure below 69.5 K.

Temperature-dependent resistance measurements on single crystals of
${\mathrm{Ca}}{\mathrm{Mn}}_{2}{\mathrm{P}}_{2}$ were performed in a diamond-
anvil cell (DAC) with the aim to check whether the compound could be
metallized under pressure. However, as shown in fig. 5

[scale=0.29]Fig5.eps

Figure 5: Temperature dependence of electrical resistances at different
pressures for ${\mathrm{Ca}}{\mathrm{Mn}}_{2}{\mathrm{P}}_{2}$. The inset
shows the pressure dependence of the first-order transition temperature.

, for pressures up to 9.1 GPa, the compound remains an insulating state, and
there is no abrupt change in the overall $\rho(T)$ behavior. Meanwhile, the
first order transition is robust. Denoting the resistively drop temperature in
$\rho(T)$ , as the transition temperature $T^{*}$, we find that $T^{*}$
increases almost linearly in the pressure range investigated. It is in
contrast to the suppression of AFM transition by external pressure in the
metallic 122 iron-arsenides [21, 22]. Given that the first order transition is
of structural origin, why the low temperature phase becomes more favorable
with increasing pressure remains unclear.

Figure 6

[scale=0.32]Fig6.eps

Figure 6: Isothermal magnetization of
${\mathrm{Ca}}{\mathrm{Mn}}_{2}{\mathrm{P}}_{2}$ at H=1 T with magnetic field
applied along both the $ab$ plane and $c$ axis.

shows the temperature dependent magnetization of single crystalline
${\mathrm{Ca}}{\mathrm{Mn}}_{2}{\mathrm{P}}_{2}$ measured in a magnetic field
$H$=1 T applied in both the $ab$ plane ($H//ab$) and along the
$c$-axis($H//c$). Due to the small mass of tested crystal, the data is noisy
and not normalized to molar unit. However the overall behavior is clear. The
upturn below 30 K is regarded as a contribution from paramagnetic impurities.
Other than that, there is no apparent feature at the first order transition
temperature, neither any sign of an AFM transition below 400 K. In addition,
the anisotropy in magnetization is small. Previous studies on other 122-type
Mn pnictides suggested that due to the competition among different exchange
interaction energies of Mn sites, these compounds tended to have AFM
correlations or develop a long range AFM order. Here for
${\mathrm{Ca}}{\mathrm{Mn}}_{2}{\mathrm{P}}_{2}$, the magnetization increases
with increasing temperature, even if a long range AFM order developed at
temperature higher than 400 K, the small anisotropy with field applied along
two perpendicular directions could rule out a simple AFM type moments
alignment along crystal principle axis. Other possible scenarios would include
wide fluctuation window of AFM in which the long range magnetic order has not
been established or frustration of Mn moment on the hexagonal lattice.

##  4 Conclusions

We have successfully synthesized high quality
${\mathrm{Ca}}{\mathrm{Mn}}_{2}{\mathrm{P}}_{2}$ single crystals by the Sn
flux method. The compound is an insulator and crystalizes in the trigonal
${\mathrm{Ca}}{\mathrm{Al}}_{2}{\mathrm{Si}}_{2}$-type structure with lattice
parameters $a$=4.0973(1) Å, and $c$=6.8473(3) Å at room temperature. Heat
capacity, Raman spectra and electrical resistivity $\rho(T)$ measurements
reveal that ${\mathrm{Ca}}{\mathrm{Mn}}_{2}{\mathrm{P}}_{2}$ undergoes a
structural phase transition at $T$=69.5 K and the transition temperature
increase with increasing pressure.
${\mathrm{Ca}}{\mathrm{Mn}}_{2}{\mathrm{P}}_{2}$ is unique in that it is the
first compound in the hexagonal 122-type Mn pnictides which has a structural
transition at low temperature other than apparent magnetic correlation.

###### Acknowledgements.

We are grateful to Youting Song (Institute of Physics) for his help in the
single crystal X-ray diffraction experiments. This work was supported by the
National Key Research and Development of China (Grant Nos. 2017YFA0302901, and
2017YFA0302903), the National Natural Science Foundation of China (Grant Nos.
11674375, 11634015 and 11921004), the Strategic Priority Research Program of
the Chinese Academy of Sciences (Grant No. XDB33010100), and the Postdoctoral
Science Foundation of China.

## References

  * [1] Wu W., Cheng J. G., Matsubayashi K., Kong P. P., Lin F. K., Jin C. Q., Wang N. L., Uwatoko Y. Luo J. L. Nat. Commun.520145508.
  * [2] Bao J. K., Liu J. Y., Ma C. W., Meng Z. H., Tang Z. T., Sun Y. L., Zhai H. F., Jiang H., Bai H., Feng C. M., Xu Z. A. Cao G. H. Phys. Rev. X52015011013.
  * [3] Tang Z. T., Bao J. K., Liu Y., Sun Y. L., Ablimit A., Zhai H. F., Jiang H., Feng C. M., Xu Z. A. Cao G. H. Phys. Rev. B912015020506.
  * [4] Mu Q. G., Ruan B. B., Pan B. J., Liu T., Yu J., Zhao K., Chen G. F., Ren Z. A. Phys. Rev. B962017140504.
  * [5] Shao Y. T., Wu X. X., Wang L., Shi Y. G., Hu J. P. Luo J. L. Europhys. Lett.123201857001.
  * [6] Wu W., Liu K., Li Y. J., Yu Z.-H., Wu D. S., Shao Y. T., Na S. H., Li G., Huang R. Z., Xiang T. Luo J. L. Nati. Sci. Rev.720201.
  * [7] Li Y. J., Wu W., Liu K., Yu Z. H., Wu D. S., Shao Y. T., Na S. H., Li G., Zheng P., Xiang T Luo J. L. Europhys. Lett.128201967002.
  * [8] Nigro A., Marra P., Autieri C., Wu W., Cheng J. G., Luo J. L. Noce C. Europhys. Lett.125201957002.
  * [9] Cheng J. G., Matsubayashi K., Wu W., Sun J. P., Lin F. K., Luo J. L. Uwatoko Y. Phys. Rev. Lett.1142015117001.
  * [10] Na S. H., Wu W. Luo J. L. Chin. Phys. Lett.372020087301.
  * [11] Rotter M., Tegel M., Johrendt D., Schellenberg I., Hermes W. Pottgen R. Phys. Rev. B782008020503(R).
  * [12] Chen G. F., Li Z., Dong J., Li G., Hu W. Z., Zhang X. D., Song X. H., Zheng P., Wang N. L. Luo J. L. Phys. Rev. B782008224512.
  * [13] Krellner C., Caroca-Canales N., Jesche A., Rosner H., Ormeci A. Geibel C. Phys. Rev. B782008100504(R).
  * [14] Ni N., Budko S. L., Kreyssig A., Nandi S., Rustan G. E., Goldman A. I., Gupta S., Corbett J. D., Kracher A. Canfield P. C. Phys. Rev. B782008014507.
  * [15] Ronning F., Klimczuk T., Bauer E. D., Volz H. Thompson J. D. J. Phys.: Condens. Matter202008322201.
  * [16] Goldman A. I., Argyriou D. N., Ouladdiaf B., Chatterji T., Kreyssig A., Nandi S., Ni N., Bud’ko S. L., Canfield P. C. McQueeney R. J. Phys. Rev. B782008100506(R).
  * [17] Tegel M., Rotter M., Weiss V., Schappacher F. M., Poettgen R. Johrendt D. J. Phys.: Condens. Matter202008452201.
  * [18] Ren Z., Zhu Z., Jiang S., Xu X., Tao Q., Wang C., Feng C., Cao G. Xu Z. Phys. Rev. B782008052501.
  * [19] Chen G. F., Li Z., Li G., Hu W. Z., Dong J., Zhang X. D., Zheng P., Wang N. L. Luo J. L. Chin. Phys. Lett.2520083403.
  * [20] Rotter M., Tegel M. Johrendt D. Phys. Rev. Lett.1012008107006.
  * [21] Kimber S. A. J., Kreyssig A., Zhang Y. Z., Jeschke H. O., Valenti R., Yokaichiya F., Colombier E., Yan J., Hansen T. C., Chatterji T., McQueeney R. J., Canfield P. C., Goldman A. I. Argyriou D. N. Nat. Mater.82009471.
  * [22] Wu J. J., Lin J. F., Wang X. C., Liu Q. Q., Zhu J. L., Xiao Y. M., Chow P. Jin C. Proc. Natl. Acad. Sci. USA110201317263.
  * [23] Sangeetha N. S., Cuervo-Reyes E., Pandey A. Johnston D. C. Phys. Rev. B942016014422.
  * [24] Brock S. L., Greedan J. E. Kauzlarich S. M. J. Solid State Chem.1131994303.
  * [25] Sefat A. S., Singh D. J., Jin R., McGuire M. A., Sales B. C. Mandrus D. Phys. Rev. B792009024512.
  * [26] Li Y., Luo Y., Li L., Chen B., Xu X., Dai J., Yang X., Zhang L., Cao G. Xu Z. J. Phys.: Condens. Matter262014425701.
  * [27] Saparov B. Sefat A. S. J. Solid State Chem.1912012213.
  * [28] Bauer E. D., Ronning F., Scott B. L. Thompson J. D. Phys. Rev. B782008172504.
  * [29] Ronning F., Kurita N., Bauer E. D., Scott B. L., Park T., Klimczuk T., Movshovich R. Thompson J. D. J. Phys.: Condens. Matter202008342203.
  * [30] Wang Z. W., Yang H. X., Tian H. F., Shi H. L., Lu J. B., Qin Y. B., Wang Z. Li J. Q. J. Phys. Chem. Solids722011457.
  * [31] Sangeetha N. S., Pandey A., Benson Z. A. Johnston D. C. Phys. Rev. B942016094417.
  * [32] Das P., Sangeetha N. S., Pandey A., Benson Z. A., Heitmann T. W., Johnston D. C., Goldman A. I. Kreyssig A. J. Phys.: Condens. Matter292017035802.
  * [33] Simonson J. W., Smith G. J., Post K., Pezzoli M., Kistner-Morris J. J., McNally D. E., Hassinger J. E., Nelson C. S., Kotliar G., Basov D. N. Aronson M. C. Phys. Rev. B862012184430.
  * [34] Gibson Q. D., Wu H., Liang T., Ali M. N., Ong N. P., Huang Q. Cava R. J. Phys. Rev. B912015085128.
  * [35] Singh Y, Ellern A. Johnston D. C. Phys. Rev. B792019094519.
  * [36] Saparov B. Sefat A. S. J. Solid State Chem.204201332.
  * [37] Sheldrick G. M. Acta Crystallogr. C7120153.
  * [38] Mewis A. Z. Naturforsch.33b1978606.
  * [39] Jain A., Ong S. P., Hautier G., Chen W., Richards W. D., Dacek S., Cholia S., Gunter D., Skinner D., Ceder G. Persson K. A. APL Materials1(1)2013011002.
  * [40] Pearson G. L. Bardeen J. Phys. Rev.751949865.
  * [41] Kittel C. Introduction to Solid State PhysicsWiley Press, Hoboken, NJ2005.


# Multi-agent Black-box Optimization using a Bayesian Approach to Alternating
Direction Method of Multipliers

Dinesh Krishnamoorthy1 Joel A. Paulson2 1Department of Mechanical Engineering,
Eindhoven University of Technology, Eindhoven, 5600MB, The Netherlands.
2Department of Chemical and Biological Engineering, The Ohio State University,
Columbus, OH 43210, USA.
(email<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Bayesian optimization (BO) is a powerful black-box optimization framework that
looks to efficiently learn the global optimum of an unknown system by
systematically trading-off between exploration and exploitation. However, the
use of BO as a tool for coordinated decision-making in multi-agent systems
with unknown structure has not been widely studied. This paper investigates a
black-box optimization problem over a multi-agent network coupled via shared
variables or constraints, where each subproblem is formulated as a BO that
uses only its local data. The proposed multi-agent BO (MABO) framework adds a
penalty term to traditional BO acquisition functions to account for coupling
between the subsystems without data sharing. We derive a suitable form for
this penalty term using alternating directions method of multipliers (ADMM),
which enables the local decision-making problems to be solved in parallel (and
potentially asynchronously). The effectiveness of the proposed MABO method is
demonstrated on an intelligent transport system for fuel efficient vehicle
platooning.

###### keywords:

Distributed optimization, Bayesian Optimization, Multi-agent systems

††thanks: D.K. acknowledges support from the Research Council of Norway
through the IKTPLUSS program and Eindhoven University of Technology. J.A.P.
acknowledges support from the National Science Foundation Grant 2036982.

##  1 Introduction

Multi-agent decision-making problems (often formulated as optimization tasks)
arise in a wide range of application areas, such as in collaborative robotics,
intelligent transport systems, integrated energy systems, sensor networks, and
modular manufacturing, where a collection of agents work together to achieve
one or more well-defined goals. A common framework for multi-agent decision-
making is the decomposition-coordination approach, where the agents solve
their own local subproblems, and a central coordinator accounts for the
coupling/interactions between the different agents, thus facilitating a
coordinated action (Bertsekas, 2016; Lasdon, 2002). One such popular
distributed optimization approach is the alternating directions method of
multipliers (ADMM) Boyd et al. (2011).

In the traditional optimization framework, each local agent requires
analytical expressions (i.e., models) that relate the local decision variables
to the local objectives and constraints in every subsystem. Obtaining reliable
models is a time consuming and tedious work, and modeling effort is often one
of the biggest bottlenecks in many application areas (Chachuat et al., 2009).
Accurate models may also be challenging to obtain in some applications, due to
lack of domain knowledge. This requirement hinders the use of traditional
distributed optimization methods that rely on exploiting the structure of such
models (in the form of gradients or relaxations) on applications for which
accurate models are not readily available.

A class of optimization methods known as “black box” optimization is a
promising alternative that circumvents the need for detailed mathematical
models. Here, the system is considered to be an unknown “black-box” that can
be queried through a series of simulations or experiments. The data observed
from such evaluations can then be used to find the optimum. Bayesian
optimization (BO) is one such zeroth-order black-box optimization scheme,
where the optimum is learned by sequentially interacting with the system
(Kushner, 1964; Shahriari et al., 2015; Jones et al., 1998).

However, in the context of multi-agent systems, formulating the subproblems
using BO does not account for the coupling and the interactions between the
different agents. In BO, the next sequence of actions are computed by
optimizing a so-called “acquisition function” induced from a probabilistic
surrogate model for the local objective function. Several different
acquisition functions have been proposed in the BO literature, which
fundamentally attempt to trade-off exploration (learning for the future) and
exploitation (immediately advancing toward goal). Interested readers are
referred to the review papers (Shahriari et al., 2015; Frazier, 2018) for
detailed discussions on the different type of acquisition functions and the
relationship between them. In the distributed setting, however, the local
acquisition function in each subsystem does not account for the coupling to
the other subsystems. A truly distributed BO algorithm requires the local
acquisition functions to be modified so that, e.g., a central coordinator can
influence the local BO decision-making in each subsystem. Despite the growing
popularity of BO, decomposition-coordination frameworks for BO has received
very limited attention in the literature.

To address this gap, this paper proposes a multi-agent Bayesian optimization
(MABO) framework that is capable of solving black-box optimization problems
defined over a multi-agent network coupled via shared variables or
constraints. An important and realistic assumption made about the network is
that direct data sharing is not allowed between the local agents due to data
privacy concerns or bandwidth limitations. Such constraints would not be
satisfied by traditional “centralized” BO methods, but is naturally handled by
our proposed MABO approach, which augments the local acquisition functions
with additional penalty terms expressed directly in terms of limited
information received from the central coordinator. MABO can be used in
combination with any existing acquisition function from the BO literature and,
thus, provides a bridge between the developments in BO and distributed
optimization.

The reminder of the paper is organized as follows. Section 2 introduces the
preliminaries of Bayesian optimization. Section 3 describes the proposed MABO
framework and corresponding penalty terms that are derived using ADMM. The
effectiveness of the MABO framework is then demonstrated on a fuel efficient
vehicle platooning problem in Section 4, wherein the goal is to minimize
overall fuel consumption of the platoon by selecting an optimal cruising
speed. Lastly, we conclude the paper and discuss some interesting directions
for future work in Section 5.

##  2 Preliminaries - Bayesian Optimization

Consider the optimization problem

$\min_{x\in\mathcal{X}}f(x)$ (1)

where $x\in\mathcal{X}$ is the decision variable chosen from a known compact
set $\mathcal{X}$ and $f:\mathcal{X}\rightarrow\mathbb{R}$ is some cost black-
box function whose mathematical structure is unknown. Bayesian optimization
(BO) consists of the following two steps:

  1. 1.

Learning: A probabilistic surrogate model, typically a Gaussian process (GP),
is derived for the unknown cost function $f(x)$. The prior distribution for
the surrogate model is updated every time that a noisy observation of the cost
$y=f(x)+w$, where $w$ is a stochastic noise term, is available by conditioning
on the observations (Rasmussen and Williams, 2006).

  2. 2.

Decision-making: The posterior surrogate model is then used to induce a
suitable acquisition function $\alpha:\mathcal{X}\rightarrow\mathbb{R}$, which
provides a measure of how desirable querying a point $x\in\mathcal{X}$ is
expected to be with respect to optimizing the cost. Using this, the next
action to evaluate is computed by solving111Note that we choose to define the
acquisition function $\alpha(\cdot)$ in a way that its a minimization problem
(as opposed to maximization problem) without loss of generality.

$x^{k+1}=\arg\min_{x\in\mathcal{X}}\;\alpha(x|\mathcal{D}^{k}),$ (2)

where $\mathcal{D}^{k}:=\\{(x^{j},y^{j})\\}_{j=1}^{k}$ is the past set of
observations at time $k$ and $\alpha(x|\mathcal{D}^{k})$ denotes the
acquisition function evaluated at $x$ induced by the posterior conditioned on
$\mathcal{D}_{k}$.

Figure 1: Illustration of Bayesian optimization scheme to minimize an unknown
function $f(x)$.

The choice of acquisition function controls the trade-off between exploration
and exploitation when selecting a suitable next action, with the goal being to
find the global optimum in as few iterations as possible. There are several
choices for the acquisition function such as probability of improvement,
expected improvement, upper confidence bound, and Thompson sampling. See
Shahriari et al. (2015); Frazier (2018) for a more detailed description of the
different acquisition functions.

In summary, BO is a black-box decision-making strategy where we sequentially
query the system, observe the cost, update our prior belief and choose the
next action to query. The BO loop is schematically represented in Fig. 1.

##  3 Proposed Multi-agent Bayesian optimization (MABO) Framework

Now, let us consider a multi-agent version of (1) that has the following
additively separable cost function

$\min_{x\in\mathcal{X}}\;f(x)=\textstyle\sum_{i=1}^{N}f_{i}(x),$ (3)

where $f_{i}:\mathcal{X}\to\mathbb{R}$ denotes a local cost function for the
$i^{\text{th}}$ subsystem. We make the following assumption about these local
cost functions.

###### Assumption 1

The local cost functions $f_{i}(x)$ are unknown but can be evaluated at any
desired $x\in\mathcal{X}$ for all $i=1,\ldots,N$.

If this was our only assumption, we could apply traditional BO methods (Fig.
1) on the aggregated set of data from all subsystems. It is worth noting that
an improved version of this algorithm can easily be developed by exploiting
the additive cost structure in (3), as discussed in Kudva et al. (2022).
However, such methods are centralized in the sense that they require the cost
data to be gathered from all local subsystems at every iteration. There are
many situations where this is not desired or even possible, which we summarize
in the following additional assumption.

###### Assumption 2

Data cannot be shared across subsystems such that subsystem $i$ cannot
directly use data from any other subsystem $j\neq i$.

Figure 2: Multi-agent Bayesian optimization (MABO) scheme where each
subproblem minimizes its unknown local cost function $f_{i}(x_{i})$ using a
modified BO method. Each subproblem learns only from the local data set
$\mathcal{D}_{i}$ which is not shared among the different subsystems. The
acquisition function in each subproblem has additional penalty terms as a
function of the master variables $\theta_{i}$ that enables coordinated action
of the different subsystems.

A straightforward approach to tackle (3), under Assumptions 1–2, is to
separately apply BO to each subproblem. That is, learn a GP surrogate model
for each local cost function $f_{i}(x_{i})$ based only on its local
observations $\mathcal{D}_{i}^{k}=\\{(x_{i}^{j},y_{i}^{j})\\}_{j=1}^{k}$. Any
suitable acquisition function $\alpha_{i}(x_{i}|\mathcal{D}^{k}_{i})$ induced
by the GP posterior is then used for local decision-making by each agent.
However, computing the next action by optimizing only the acquisition function
induced locally

$x_{i}^{k+1}=\arg\min_{x_{i}\in\mathcal{X}}\;\alpha_{i}(x_{i}|\mathcal{D}_{i}^{k}),$

results in each of the different subsystems taking independent decisions
without any coordination. Solving the multi-agent version of the problem (3)
requires a coordinated solution in which $x_{1}^{\star}=\cdots=x_{N}^{\star}$,
which is not necessarily achieved with the aforementioned strategy.

To enable the required coordination, we propose to incorporate additional
penalty terms to the local acquisition function. That is, each subsystem now
computes its next action $x_{i}^{k+1}$ by optimizing

$x_{i}^{k+1}=\arg\min_{x_{i}\in\mathcal{X}}\;\alpha_{i}(x_{i}|\mathcal{D}^{k}_{i})+\delta_{i}(x_{i},\theta_{i}^{k}),$
(4)

where $\theta^{k}_{i}\in\mathbb{R}^{n_{\theta}}$ denotes the information that
is communicated to the $i^{\text{th}}$ subsystem at the $k^{\text{th}}$
iteration (e.g. by a central coordinator), $x_{i}\in\mathcal{X}$ denotes the
local decision variables, $\alpha_{i}(\cdot)$ is any suitable acquisition
function used to optimize the unknown local cost $f_{i}(\cdot)$, and
$\delta_{i}:\mathcal{X}\times\mathbb{R}^{n_{\theta}}\rightarrow\mathbb{R}$ is
the augmented penalty term that can influence the local decisions. We propose
to derive an expression for $\delta_{i}(\cdot)$ using the alternating
direction method of multipliers (ADMM) framework, as described below.

The optimization problem (3) can be equivalently rewritten by introducing a
new variable $x_{0}$ as

$\displaystyle\min_{x_{0},\\{x_{i}\\}\in\mathcal{X}}\;\textstyle\sum_{i=1}^{N}f_{i}(x_{i})\;\textup{s.t.}\;x_{i}=x_{0}\quad\forall
i=1,\dots,N.$ (5)

The augmented Lagrangian for (5) is given by

$\min_{x_{0},\\{x_{i}\\}\in\mathcal{X}}\;\sum_{i=1}^{N}\left[f_{i}(x_{i})+\lambda_{i}^{\mathsf{T}}(x_{i}-x_{0})+\frac{\rho}{2}\|x_{i}-x_{0}\|^{2}\right],$
(6)

where $\lambda=(\lambda_{1},\ldots,\lambda_{N})$ denotes the set of dual
variables (or Lagrange multipliers) associated with the equality consraints in
(5) and $\rho>0$ is a penalty parameter. The ADMM scheme, which is derived
from the augmented Lagrangian in (6), then consists of the following
iterations (Boyd et al., 2011, Chapter 7)

$\displaystyle x_{0}^{k+1}$
$\displaystyle=\frac{1}{N}\sum_{i=1}^{N}\left[x_{i}^{k}+\frac{\lambda_{i}^{k}}{\rho}\right],$
(7a) $\displaystyle x_{i}^{k+1}$
$\displaystyle=\arg\min_{x_{i}\in\mathcal{X}}f_{i}(x_{i})+\lambda_{i}^{k\mathsf{T}}x_{i}+\frac{\rho}{2}\|x_{i}-x_{0}^{k+1}\|^{2},$
(7b) $\displaystyle\lambda_{i}^{k+1}$
$\displaystyle=\lambda_{i}^{k}+\rho(x_{i}^{k+1}-x_{0}^{k+1}),$ (7c)

Since $f_{i}(x_{i})$ in each subproblem (7b) is unknown under Assumption 1, we
propose to replace it with a local acquisition function
$\alpha_{i}(x_{i}|\mathcal{D}_{i}^{k})$. In this case, (7b) has the same form
as (4) with the penalty term equal to

$\displaystyle\delta_{i}(x_{i},\theta_{i}^{k})=\lambda_{i}^{k\mathsf{T}}x_{i}+\frac{\rho}{2}\|x_{i}-x_{0}^{k+1}\|^{2},$
(8)

where $\theta^{k}_{i}:=(x^{k+1}_{0},\lambda_{i}^{k})$ denotes the coordinating
variables, which are updated as shown in (7a) and (7c) at every iteration;
these updates are very simple, can be stored with little memory, and ensure
that Assumption 2 is satisfied. We refer to this overall proposed framework as
multi-agent BO (MABO), which is illustrated in Fig. 2. We also provide a
complete algorithmic description of MABO in Algorithm 1, which has several
unique features that we highlight in the series of remarks below.

Algorithm 1 Multi-agent Bayesian optimization.

1:$N$ local GP models for the local cost functions $f_{i}(x_{i})$, initial
values $x_{i}^{0}$ and $\lambda_{i}^{0}$, and penalty $\rho>0$

2:for ADMM iterates $k=0,1,2,\dots$ do

3:
$x_{0}^{k+1}\leftarrow\frac{1}{N}\sum_{i=1}^{N}\left[x_{i}^{k}+\frac{\lambda_{i}^{k}}{\rho}\right]$
$\triangleright$ Central collector

4: for subproblems $i=1,\dots,N$ do (in parallel)

5: $\theta^{k}_{i}\leftarrow[x_{0}^{k+1},\lambda_{i}^{k}]^{\mathsf{T}}$

6:
$\delta_{i}(x_{i},\theta_{i}^{k})\leftarrow\lambda_{i}^{k\mathsf{T}}x_{i}+\frac{\rho}{2}\|x_{i}-x_{0}^{k+1}\|^{2}$

7: for local BO iterates $j=1,2,\dots$ do

8: Update local GP model using $\mathcal{D}_{i}^{k}$

9: Induce any acquisition function $\alpha_{i}(x_{i}|\mathcal{D}_{i}^{k})$

10:
$x_{i}^{k+1,j}\leftarrow\arg\min_{x_{i}\in\mathcal{X}}\alpha_{i}(x_{i}|\mathcal{D}_{i}^{k})+\delta_{i}(x_{i},\theta^{k}_{i})$

11: Query and observe the local cost $f_{i}(x_{i}^{k+1,j})$

12:
$\mathcal{D}_{i}^{k}\leftarrow\mathcal{D}_{i}^{k}\cup\\{(x_{i}^{k+1,j},f_{i}(x_{i}^{k+1,j}))\\}$

13: end for

14:
$\lambda_{i}^{k+1}\leftarrow\lambda_{i}^{k}+\rho(x_{i}^{k+1,j}-x_{0}^{k+1})$

15: end for

16:end for

17:Optimal solution $x_{0}$

######  Remark 1 (Data privacy)

The GP model in each subsystem is conditioned only on the local cost
measurement, which need not be shared with the other subsystems. This enables
a collection of agents to take coordinated actions without sharing its local
data with one another, clearly respecting data privacy concerns or any other
data sharing limitations (such as bandwidth or connection limitations).

######  Remark 2 (Choice of acquisition function)

Note that the choice of the local acquisition function is independent from the
coordinating term $\delta_{i}(\cdot)$. Hence any acquisition function proposed
in the standard BO framework can be directly used by our MABO approach.
Moreover, the choice of the acquisition function can also be different between
the different subsystems. Simply put, in a multi-agent network, each agent is
free to choose its local acquisition function, without affecting the structure
of the MABO strategy.

######  Remark 3 (Multiple BO iterations)

It can be seen that our proposed decomposable MABO framework (Algorithm 1) has
a hierarchical structure with two levels of iterations, namely, the local BO
iterations within each subproblem, and the ADMM iterations between the central
coordinator and the subproblems. As such, MABO can also be implemented with
several BO iterations per ADMM iteration. For example, after every ADMM
iterate that updates $\theta^{k}_{i}$, the local subproblems can run one or
several local Bayesian optimization iterates for a given $\theta^{k}_{i}$.
Furthermore, the different subproblems can also run different number of BO
iterations for each ADMM iterate. This is especially beneficial where the
local evaluation time and budgets are different in the different agents. This
is explicitly shown in Algorithm 1, where the local BO iterations $j$ can be
run for any desired finite number of steps.

######  Remark 4 (Optimal sharing)

Although we have focused our MABO derivation on the consensus problem in (5),
the approach extends to optimal sharing problems that have coupling
constraints of the form $\sum_{i}A_{i}x_{i}=0$. The main modification needed
is to update the form of the penalty term $\delta_{i}(\cdot)$ whose structure
is shown in Boyd et al. (2011, Section 7.3). Many important applications
involving the distribution of a shared resource would fall into this class of
problems (e.g., allocation of energy from a centralized repository in a
chemical process).

###### Remark 5

Although we considered only the ADMM framework in this section, (4) also
allows for other decomposition strategies such as primal decomposition and
Lagrangian decomposition, with $\delta_{i}(\cdot)$ and $\theta_{i}$ chosen
accordingly.

##  4 Illustrative example:
Fuel Efficient Vehicle Platooning

We demonstrate the performance of MABO on the problem of fuel-efficient
vehicle platooning, which is a multi-agent system consisting of several
vehicles that take coordinated action to optimize the overall fuel efficiency.
Here, “platooning” refers to a collection of vehicles that drive in a single-
file fashion at close inter-vehicle distance in highways. Platooning of heavy-
duty vehicles (e.g., trucks used for freight transport) has several benefits
such as reduced traffic congestion, greater fuel economy due to reduced air
drag friction, and fewer traffic collisions (Franke et al., 1995; Browand et
al., 2004).

Early work on truck platooning focused on tracking inter-vehicle distances
using adaptive cruise control (Hao and Barooah, 2013); however, there has been
an increasing interest in optimizing the speed of the platoon to minimize
overall fuel consumption. In fact, a recent experimental study showed that,
although the fuel economy improved at all platooning speeds, certain speeds
resulted in the best overall fuel economy (Lammert et al., 2014). However, the
effect of cruising speed on fuel consumption is not constant and depends on
several factors such as freight load, road and driving conditions, time since
last engine maintenance, etc., to name a few. Therefore, developing accurate
models that can be used to find the optimal platoon speed is challenging.

An attractive alternative is to resort to a black-box optimization strategy.
Black-box formulations are particularly useful when the platoon consists of
vehicles that have different configurations and/or manufacturers since it is
non-obvious how to coordinate different types of vehicles to minimize overall
fuel consumption when each isolated vehicle has a unique optimum speed. Such
problems can be formulated in the form of (3) where $x_{i}$ and
$f_{i}(x_{i})$, respectively, denote the speed and fuel consumption of the
$i^{\text{th}}$ vehicle. Here, we simulate the fuel consumption of a platoon
of $N=7$ vehicles using the approach derived by Sobrino et al. (2016), which
leads to a model of the form

$\displaystyle
f_{i}(x_{i})=a_{i}+\frac{b_{i}}{x_{i}}+c_{i}x_{i}+d_{i}x_{i}^{2},$ (9)

where speed $x_{i}$ has units of [km/h] and fuel consumption $f_{i}(x_{i})$
has units of [g/veh-km]. The nominal model parameters are given by the
experimentally calibrated values in (Sobrino et al., 2016, Section 5). To
emulate $N=7$ completely different trucks, we sample perturbed model
parameters from a uniform distribution defined by $\pm 20\%$ of these nominal
values. We highlight that (9) is used only as a source of data generation, as
the parameter values and structure are unknown in accordance with Assumption
1.

Figure 3: The local Bayesian optimization of one of the subsystems at iterates
$k=1,2,$ and 20. Top subplots show the Gaussian process of the local cost.
Bottom subplots show the acquisition function with (solid green) and without
(dashed green) coordination. The next query point is shown as a red dot.

The proposed MABO approach (Algorithm 1) was implemented in Python using the
GPy package (GPy, since 2012) to construct and train the GP models using a
squared exponential (SE) kernel. We first use a lower confidence bound (LCB)
acquisition function defined as follows

$\displaystyle\alpha_{i}(x_{i}|\mathcal{D}_{i}^{k}):=\mu_{i}^{k}(x_{i})-\sqrt{\beta}\sigma_{i}^{k}(x_{i}),$
(10)

where $\mu_{i}^{k}$ and $\sigma_{i}^{k}$, respectively, denote the posterior
mean and standard deviation predicted by the GP model for the $i^{\text{th}}$
subsystem at the $k^{\text{th}}$ iteration and $\beta$ is a scaling parameter
used to trade-off between exploration and exploitation. Based on previous
work, we select $\beta=4$ throughout this work (Srinivas et al., 2009). We
selected a penalty value of $\rho=10$ and only run $j=1$ local BO iteration
for every ADMM iteration.

We first demonstrate the effect of the coordination term in Fig. 3, which
plots the GP model and corresponding acquisition functions for a single
vehicle at three iterations $k\in\\{1,2,20\\}$. From the bottom subplots, we
see that the local non-perturbed acquisition function
$\alpha_{i}(x_{i}|\mathcal{D}_{i}^{k})$ (green dashed lines) behaves quite
different from the MABO acquisition
$\alpha_{i}(x_{i}|\mathcal{D}_{i}^{k})+\delta_{i}(x_{i},\theta_{i}^{k})$
(green solid line). Already at iteration $k=1$, the effect of the penalty term
$\delta_{i}(x_{i},\theta_{i}^{k})$ is evident, which skews the selection
process toward lower $x_{i}^{k+1}$ values such that the next query point is
not too far from the average speed of the platoon. By iteration $k=20$, the
local BO step converges to around 57 km/hr (near the true optimal value),
which significantly differs from the optimum point of the unconstrained
acquisition function.

Figure 4: Vehicle speeds when the local BO loops use probability of
improvement (leftmost subplot), expected improvement (middle left subplot),
LCB (middle right subplot), and mixed (rightmost subplot) acquisition
functions. The dashed black line represents the true solution.  Figure 5:
Comparison of platoon average speed, primal residual, and dual residual when
local BO loops use probability of improvement (orange), expected improvement
(blue), LCB (green), and mixed (red) acquisition functions. The solid black
lines represent the ideal ADMM iterations when the local functions are
perfectly known.

We repeat the same simulation test using two additional types of acquisition
functions, probability of improvement (PI) and expected improvement (EI), in
addition to LCB. Furthermore, to demonstrate the flexibility of MABO, we also
test a “mixed” case wherein the different vehicles use different acquisition
functions (i.e., vehicles 1 and 2 use EI, vehicles 3 and 4 use PI, vehicles 5
and 6 use LCB, and vehicle 7 using a greedy mean-based policy). The individual
vehicle speeds obtained at each MABO iteration using the different acquisition
function test cases are shown in Fig. 4. We also show the true optimum with a
black dashed line to serve as a benchmark to compare against. The average
speed of the platoon $x_{0}^{k}$, the primal residual
$r^{k}:=\sum_{i}\|x_{i}^{k}-x_{0}^{k}\|_{2}^{2}$, and the dual residual
$s^{k}:=N\rho^{2}\|x_{0}^{k}-x_{0}^{k-1}\|_{2}^{2}$ for the different cases
are shown in Fig. 5 to compare the effect of the choice of the acquisition
function on the progress of the ADMM iterations. It can be seen that the
progress of the ADMM iterates are not heavily affected by the choice of the
acquisition function. We also benchmark this against the traditional model-
based ADMM with full model information in each subsystem, which acts as a
lower bound on the achievable residuals. This is shown in solid black lines.
This empirically quantifies the additional iterations needed to explore the
unknown function when using BO in the subsystems.

##  5 Conclusion and Future work

This paper presents a novel approach for efficiently finding global optimum
solutions of multi-agent systems with unknown (black-box) structure and
limitations on information sharing by combining elements of distributed
optimization with Bayesian optimization (BO). As opposed to independently
applying BO in each of the subsystems, the key idea is to modify the local
acquisition (or expected utility) function with additional penalty terms that
are properly set according to a master coordinator, which requires limited
communication bandwidth and memory storage. Our proposed multi-agent BO (MABO)
framework is broadly applicable in the sense that it works for different
network structures (e.g., consensus and sharing formulations) and allows for
the use of any well-defined acquisition function within the local subproblems.
Furthermore, the local subproblems can be straightforwardly formulated using
constrained BO (Gelbart et al., 2014; Paulson and Lu, 2022) or safe BO methods
(Berkenkamp et al., 2021; Krishnamoorthy and Doyle III, 2022) whenever local
black-box constraints are present in a given subsystem. The proposed MABO
approach is demonstrated on a vehicle platooning problem, where the goal is to
select a common platoon speed for multiple (different configuration) vehicles
such that the overall fuel consumption is minimized. We show that MABO is
capable of quickly identifying the optimal platoon speed, even with no prior
knowledge of the relationship between speed and fuel consumption in any of the
vehicles in the platoon. This work represents a first step toward a
decomposable BO framework and future work should involve better understanding
the theoretical properties (e.g., rate of convergence) of the proposed method.
In addition, it would be interesting to explore the impact of the update rule
chosen for the central coordinator on performance, as there may exist
strategies that accelerate the speed of convergence in certain applications.

## References

  * Berkenkamp et al. (2021) Berkenkamp, F., Krause, A., and Schoellig, A.P. (2021).  Bayesian optimization with safety constraints: safe and automatic parameter tuning in robotics.  _Machine Learning_ , 1–35.
  * Bertsekas (2016) Bertsekas, D.P. (2016).  _Nonlinear Programming_.  Athena Scientific, Belmont, MA, 3 edition.
  * Boyd et al. (2011) Boyd, S., Parikh, N., Chu, E., Peleato, B., and Eckstein, J. (2011).  Distributed optimization and statistical learning via the alternating direction method of multipliers.  _Foundations and Trends® in Machine Learning_ , 3(1), 1–122.
  * Browand et al. (2004) Browand, F., McArthur, J., and Radovich, C. (2004).  Fuel saving achieved in the field test of two tandem trucks.
  * Chachuat et al. (2009) Chachuat, B., Srinivasan, B., and Bonvin, D. (2009).  Adaptation strategies for real-time optimization.  _Computers & Chemical Engineering_, 33(10), 1557–1567.
  * Franke et al. (1995) Franke, U., Bottiger, F., Zomotor, Z., and Seeberger, D. (1995).  Truck platooning in mixed traffic.  In _Proceedings of the Intelligent Vehicles’ 95. Symposium_ , 1–6. IEEE.
  * Frazier (2018) Frazier, P.I. (2018).  A tutorial on bayesian optimization.  _arXiv preprint arXiv:1807.02811_.
  * Gelbart et al. (2014) Gelbart, M.A., Snoek, J., and Adams, R.P. (2014).  Bayesian optimization with unknown constraints.  In _Proceedings of the Thirtieth Conference on Uncertainty in Artificial Intelligence_ , UAI’14, 250–259. AUAI Press.
  * GPy (since 2012) GPy (since 2012).  GPy: A gaussian process framework in python.  http://github.com/SheffieldML/GPy.
  * Hao and Barooah (2013) Hao, H. and Barooah, P. (2013).  Stability and robustness of large platoons of vehicles with double-integrator models and nearest neighbor interaction.  _International Journal of Robust and Nonlinear Control_ , 23(18), 2097–2122.
  * Jones et al. (1998) Jones, D.R., Schonlau, M., and Welch, W.J. (1998).  Efficient global optimization of expensive black-box functions.  _Journal of Global Optimization_ , 13(4), 455–492.
  * Krishnamoorthy and Doyle III (2022) Krishnamoorthy, D. and Doyle III, F.J. (2022).  Safe Bayesian Optimization using Interior-Point Methods - Applied to Personalized Insulin Dose Guidance.  _IEEE Control System Letters (In-Press)_.
  * Kudva et al. (2022) Kudva, A., Sorouifar, F., and Paulson, J.A. (2022).  Efficient robust global optimization for simulation-based problems using decomposed Gaussian processes: Application to MPC calibration.  In _In Proceedings of the American Control Conference_ , 2091–2097. IEEE.
  * Kushner (1964) Kushner, H.J. (1964).  A new method of locating the maximum point of an arbitrary multipeak curve in the presence of noise.  _Journal of Basic Engineering_ , 86, 97–106.
  * Lammert et al. (2014) Lammert, M.P., Duran, A., Diez, J., Burton, K., and Nicholson, A. (2014).  Effect of platooning on fuel consumption of class 8 vehicles over a range of speeds, following distances, and mass.  Technical report, National Renewable Energy Lab.(NREL), Golden, CO (United States).
  * Lasdon (2002) Lasdon, L.S. (2002).  _Optimization theory for large systems_.  Courier Corporation.
  * Paulson and Lu (2022) Paulson, J.A. and Lu, C. (2022).  COBALT: COnstrained Bayesian optimizAtion of computationaLly expensive grey-box models exploiting derivaTive information.  _Computers & Chemical Engineering_, 160, 107700.
  * Rasmussen and Williams (2006) Rasmussen, C.E. and Williams, C.K. (2006).  _Gaussian processes for machine learning_.  MIT press Cambridge, MA.
  * Shahriari et al. (2015) Shahriari, B., Swersky, K., Wang, Z., Adams, R.P., and De Freitas, N. (2015).  Taking the human out of the loop: A review of bayesian optimization.  _Proceedings of the IEEE_ , 104(1), 148–175.
  * Sobrino et al. (2016) Sobrino, N., Monzon, A., and Hernandez, S. (2016).  Reduced carbon and energy footprint in highway operations: the highway energy assessment (hera) methodology.  _Networks and Spatial Economics_ , 16(1), 395–414.
  * Srinivas et al. (2009) Srinivas, N., Krause, A., Kakade, S.M., and Seeger, M. (2009).  Gaussian process optimization in the bandit setting: No regret and experimental design.  _arXiv preprint arXiv:0912.3995_.