# Engineering Phase Competition Between Stripe Order and Superconductivity in
La1.88Sr0.12CuO4

J. Küspert<EMAIL_ADDRESS>Physik-Institut, Universität Zürich,
Winterthurerstrasse 190, CH-8057 Zürich, Switzerland  I. Biało  Physik-
Institut, Universität Zürich, Winterthurerstrasse 190, CH-8057 Zürich,
Switzerland  AGH University of Science and Technology, Faculty of Physics and
Applied Computer Science, 30-059 Kraków, Poland  R. Frison  Physik-Institut,
Universität Zürich, Winterthurerstrasse 190, CH-8057 Zürich, Switzerland  A.
Morawietz  Physik-Institut, Universität Zürich, Winterthurerstrasse 190,
CH-8057 Zürich, Switzerland  L. Martinelli  Physik-Institut, Universität
Zürich, Winterthurerstrasse 190, CH-8057 Zürich, Switzerland  J. Choi  Diamond
Light Source, Harwell Campus, Didcot OX11 0DE, United Kingdom  D. Bucher
Physik-Institut, Universität Zürich, Winterthurerstrasse 190, CH-8057 Zürich,
Switzerland  O. Ivashko  Deutsches Elektronen-Synchrotron DESY, Notkestraße
85, 22607 Hamburg, Germany.  M. v. Zimmermann  Deutsches Elektronen-
Synchrotron DESY, Notkestraße 85, 22607 Hamburg, Germany.  N. B. Christensen
Department of Physics, Technical University of Denmark, DK-2800 Kongens
Lyngby, Denmark  D. G. Mazzone  Laboratory for Neutron Scattering and Imaging,
Paul Scherrer Institut, CH-5232, Villigen-PSI, Switzerland  G. Simutis
Laboratory for Neutron and Muon Instrumentation, Paul Scherrer Institut,
CH-5232 Villigen PSI, Switzerland  A. A. Turrini  Laboratory for Neutron
Scattering and Imaging, Paul Scherrer Institut, CH-5232, Villigen-PSI,
Switzerland  L. Thomarat  Laboratory for Neutron and Muon Instrumentation,
Paul Scherrer Institut, CH-5232 Villigen PSI, Switzerland  École normale
supérieure Paris-Saclay, 91190 Gif-sur-Yvette, France  D. W. Tam  Laboratory
for Neutron Scattering and Imaging, Paul Scherrer Institut, CH-5232, Villigen-
PSI, Switzerland  M. Janoschek  Physik-Institut, Universität Zürich,
Winterthurerstrasse 190, CH-8057 Zürich, Switzerland  Laboratory for Neutron
and Muon Instrumentation, Paul Scherrer Institut, CH-5232 Villigen PSI,
Switzerland  T. Kurosawa  Department of Physics, Hokkaido University, Sapporo
060-0810, Japan  N. Momono  Department of Physics, Hokkaido University,
Sapporo 060-0810, Japan  Department of Applied Sciences, Muroran Institute of
Technology, Muroran 050-8585, Japan  M. Oda  Department of Physics, Hokkaido
University, Sapporo 060-0810, Japan  Qisi Wang<EMAIL_ADDRESS>Physik-
Institut, Universität Zürich, Winterthurerstrasse 190, CH-8057 Zürich,
Switzerland  Department of Physics, The Chinese University of Hong Kong,
Shatin, Hong Kong, China  J. Chang  Physik-Institut, Universität Zürich,
Winterthurerstrasse 190, CH-8057 Zürich, Switzerland

Unconventional superconductivity often couples to other electronic orders in a
cooperative or competing fashion. Identifying external stimuli that tune
between these two limits is of fundamental interest. Here, we show that
$c$-axis strain couples directly to the competing interaction between charge
stripe order and superconductivity in La1.88Sr0.12CuO4 (LSCO). Compressive
$c$-axis pressure amplifies stripe order within the superconducting state,
while having no impact on the normal state. By contrast, strain dramatically
diminishes the magnetic field enhancement of stripe order in the
superconducting state. This observation hints to the existence of a ground
state without phase competition. As such, $c$-axis pressure is a promising
tuning parameter in the search for ground states where superconductivity
adopts broken translational symmetry from charge ordering phenomena.
Electronic phases may coexist microscopically, either in a collaborative or
competing manner. In elementary chromium, for example, spin and charge density
wave orders coexist with commensurate ordering vectors [1, 2, 3]. A similar
spin-charge intertwined order is found in doped lanthanum-based (La-based)
cuprate superconductors [4, 5, 6, 7, 8]. Competing interaction is often found
in the context of unconventional superconductivity. For example, in CsV3Sb5
[9, 10] and BaFe2As2 [11, 12], superconductivity can be optimized through the
suppression of charge or spin density wave orders. However, an interplay
between density waves and superconductivity – at least theoretically – can
lead to a collaborative state. This state would be characterized by a
spatially modulated Cooper pair density with a commensurate wave vector.
Extensive experimental and theoretical efforts have been devoted to study this
novel superconducting state [13]. In the cuprates, pair-density-waves have
been predicted to exist in the limit where superconductivity is sufficiently
weakened [14, 15]. A general challenge is therefore to switch the coupling
between superconductivity and charge order from competing to collaborative.
Ideally, an external stimulus would tune the coupling between these two
phases.

Here, using high-energy x-ray diffraction, we show how compressive $c$-axis
uniaxial pressure enhances stripe order inside the superconducting state of
La1.88Sr0.12CuO4 (LSCO), while charge order remains unchanged in the normal
state. We furthermore discover that the magnetic field enhancement of charge
order inside the superconducting state is dramatically reduced upon
compressive $c$-axis strain application. This observation suggests a
correspondingly reduced phase competition. We thus demonstrate that $c$-axis
pressure acts directly on the coupling between charge stripe order and
superconductivity. Tuning of the phase competition between charge order and
superconductivity may be a route to pair-density-wave ground states.

Figure 1: a Schematic illustration of the scattering geometry for high-energy
x-ray diffraction on LSCO. Uniaxial pressure and magnetic field are applied
along the crystallographic $c$-axis direction as depicted. Detector read-outs
with exemplary region-of-interest for three different rotation angles $\omega$
around the vertical axis are shown. b Lattice parameters, extracted from fits
of (0, 0, $\ell$) Bragg peaks. Solid lines are least-square fits to
$c=\lambda/(2\sin{\omega})\cdot\ell$. c Charge order peak obtained by
integrating the intensities in a region-of-interest (ROI) around
$\textbf{Q}_{co}=(0.231,0,12.5)$ (blue points). A background (bg, grey points)
is estimated by a similar integration of an ROI slightly shifted off the
$(h,0,\ell)$ scattering plane. d Spin flip intensities from polarized neutron
scattering at $\textbf{Q}_{co}=(-1,-1,0)$ on LSCO ($T_{c}=27$ K), normalized
by the monitor and scaled to the integrated Bragg peak intensity of the
unstrained measurement. The drop of the curves yields the superconducting
transition temperature with and without $c$-axis uniaxial pressure (see
Methods for detailed definition of $T_{c}$).

##  I Results

Diffraction and $c$-axis strain application.
Stripe charge order in La-based cuprates manifests itself by weak reflections
at $\bf{Q}_{co}=\boldsymbol{\tau}$+$(\delta,0,0.5)$ where $\boldsymbol{\tau}$
represents fundamental Bragg peaks and $\delta\approx 1/4$ is the stripe
incommensurability [16, 6, 17, 7]. We adopted an x-ray transmission geometry
with crystalline $a$\- and $c$-axis spanning the horizontal scattering plane
as illustrated in Fig. 1a. Magnetic field and uniaxial pressure were applied
along the $c$-axis direction.

The strain as a result of uniaxial $c$-axis pressure can be directly estimated
from lattice parameter measurements. Pressure induced compression of the
$c$-axis lattice parameter is evidenced by a shift of $(0,0,n$) Bragg peaks to
larger scattering angles. Precise strain characterization utilizes multiple
such Bragg peaks $\boldsymbol{\tau}$$=(0,0,\pm 2n)$ with $n$ being an integer
(see Method section and Fig. 1b). The resulting $c$-axis lattice parameter for
strained and unstrained LSCO is exemplified in Fig. 1b. As expected, the
uniaxial c-axis pressure reduces the $c$-axis lattice parameter that in turn
lowers the superconducting transition temperature $T_{c}$ [18, 19]. Using
polarized neutron scattering we confirmed the decrease of $T_{c}$ with
compressive $c$-axis strain [20, 21, 22, 23] – see Fig. 1d and Method section.
Exploiting that $H_{c1}$ is low, we track the excess depolarisation of the
neutron beam due to flux trapped along the $c$-axis. The in-plane polarized
neutrons are depolarized by this trapped flux. Upon crossing $T_{c}$, the flux
is released. The reduced $T_{c}$ stems either from a lower superfluid density
or a weaker Cooper pairing strength.

Uniaxial $c$-axis strain effects on charge order.
X-ray diffraction intensity was collected using a two-dimensional single-
photon detector. Detector regions-of-interest (ROI) are defined such that the
signal or background of interest is covered (see Supplementary Fig. 2). We
constructed standard one-dimensional rocking curves (see Fig. 1c for
$\bf{Q}_{co}\rm=(\delta,0,12.5)$). An advantage of 2D-detectors (over point
detectors) is that a background can be estimated by slightly shifting the ROI
(grey data in Fig. 1c).

Figure 2: Charge order reflection in LSCO upon application of $c$-axis
uniaxial strain. a-d Background subtracted charge order reflections for
temperatures and momenta as indicated. Red (blue) points are recorded with
(without) compressive $c$-axis strain application. Error bars stem from
counting statistics and solid lines are fits with a split-normal distribution
including a linear background.

In Fig. 2, we show scans through $\bf{Q}_{co}$=($\delta$, 0, 12.5) and
$(\delta,0,16.5)$, with and without uniaxial $c$-axis pressure. Data for a
La2-xSrxCuO4 crystal with slightly different doping are shown in Supplementary
Fig. 3. In the normal state ($T>T_{c}$), no pressure effect is observed. By
contrast, we find a significant pressure-induced enhancement of the charge
order reflection inside the superconducting state. The correlation length and
incommensurability $\delta$ remain virtually unaffected by uniaxial pressure.
Observed shifts are within the error bars and thus negligible. From here, we
therefore consider the charge order peak amplitude as a function of
temperature, uniaxial $c$-axis pressure, and magnetic field. The peak
amplitude $I_{co}$ is extracted by fitting intensity profiles with a split-
normal distribution on a linear background – see Fig. 1c, 2, Supplementary
Fig. 5 and Methods. In Fig. 3, the intensities and fits are presented after
subtracting the background.

The temperature dependence of the charge order amplitude is shown in Fig. 3
for strained and unstrained conditions. In the absence of a magnetic field,
compressive $c$-axis pressure enhances the charge order inside the
superconducting state. The charge order peak amplitude, due to phase
competition [17], displays a cusp at $T_{c}$. The cusp is shifted to slightly
lower temperatures upon application of $c$-axis pressure. Assuming phase
competition between charge order and superconductivity, this suggests a
reduction of the superconducting transition temperature, in agreement with the
measurements in Fig. 1d. At our base temperature ($T=10$ K), the relative peak
amplitude, $I_{co}(10\,\text{K})/I_{co}(30\,\text{K})$, scales approximately
linearly with the applied strain $\epsilon_{c}$. Within the examined range of
$\epsilon_{c}$, the charge order peak amplitude increases by about 25 %.

Magnetic field effect.
Without strain, magnetic field effects on charge and spin order inside the
superconducting state have already been studied [24, 25, 26, 27, 28, 29].
Consistent with previous studies, we find an increase of the charge order
amplitude for $T<T_{c}$ (see Fig. 3b). At base temperature, the charge order
intensity scales approximately linearly with magnetic field. At 10 T, the peak
amplitude is more than doubled. This is in strong contrast to the ground state
reached through application of $c$-axis strain. Here, as depicted in Fig. 3c a
much weaker field effect is observed. The magnetic field effect for strained
and unstrained conditions is shown in Fig. 3e. Upon increasing strain
$\epsilon_{c}$, a decreasing magnetic field effect on the charge order
intensity is observed inside the superconducting state.

Figure 3: Phase competition of charge order and superconductivity. a
Temperature dependence of amplitude $I_{co}$ of the (0.231, 0, 12.5) charge
order peak without and with strain measured in setup without magnetic field.
b, c Temperature dependence for charge order peak at (0.231, 0, 12.5) in 0 T
and 10 T without and with strain, respectively. d Low-temperature (10 K)
charge order amplitude normalized to the normal state (30 K) plotted as a
function of $c$-axis strain. For the black (grey) points, the strains are
directly (indirectly) measured. Indirect measurements use the calibration
curve in Supplementary Fig. 1. Filled (unfilled) points stem from the
experimental setup with (without) a magnet. The grey shaded area is a guide to
the eye. e Magnetic field dependence of the charge order peak amplitude at
(0.231, 0, 12.5) (circles, left $y$-axis) and (0.231, 0, 16.5) (diamonds,
right $y$-axis) for different strains as indicated. Solid lines are fits to
the data – see main text. f Comparison of data obtained with and without
strain at (0.231, 0, 12.5) (circles) and (0.231, 0, 16.5) (diamonds). By
application of uniaxial strain, we roughly double the reachable parameter
space in $H/H_{c}$. Lines in a-c are guides to the eye.

##  II Discussion

There has been an interest to study the stripe order beyond the upper critical
field ($H_{c2}$) for superconductivity [30]. The magnetic field dependence of
the charge order amplitude is often modelled by
$I_{co}(\mathcal{H})=I_{co}(0)+I_{1}\mathcal{H}\ln(1/\mathcal{H})$, where
$\mathcal{H}=H/H_{c}$ and $I_{1}$ is a fitting parameter. The critical field
scale $H_{c}$ has been linked to the upper critical field $H_{c2}$ of
superconductivity [31]. We however stress that these two field scales ($H_{c}$
and $H_{c2}$) are not necessarily identical. As indicated in Fig. 3e, uniaxial
$c$-axis strain reduces $H_{c}$, driving $\mathcal{H}$ closer to the limit
$\mathcal{H}\rightarrow 1$. We show in Fig. 3f that we could double the
reachable parameter space in $\mathcal{H}$ by applying strain along $c$. By
fitting the observed charge order peak amplitude versus magnetic field, we
extract $H_{c}$ at base temperature (see fitted curves in Fig. 3e). In Fig. 4,
we plot $H_{c}$ versus strain $\epsilon_{c}$. The field scale $H_{c}$ seems to
scale approximately linearly with $c$-axis strain:
$H_{c}=H_{c}(0)-\alpha\epsilon_{c}$ with $\alpha$ being a constant. By
extrapolation, our results suggest that there exists a critical $c$-axis
strain at which $H_{c}\rightarrow 0$. At this critical strain
$\epsilon_{crit}=H_{c}(0)/\alpha$ $\approx 0.29$ %, the magnetic field effect
is expected to vanish.

Outlining our findings, we illustrate schematically in Fig. 4b the charge
order peak amplitude versus magnetic field and $c$-axis strain inside the
superconducting state. The region with highest intensities is found at high
magnetic fields and low strains. The strength of the magnetic field effect
depends on the uniaxial strain along $c$. At the highest applied $c$-axis
strain, the magnetic field effect is strongly reduced. Typically, magnetic
field effects inside the superconducting state are interpreted in terms of (i)
vortex physics [32, 33] or (ii) phase competition. The vortex density and
volume increase with increasing magnetic field strength. It is however
difficult to explain the strain-field effects in terms of vortices alone. Even
in the high-magnetic field state where vortices are present, it is not obvious
how to explain the strain effect without also involving phase competition.

Reduced (absence of) magnetic field effect would imply either that one phase
is partially suppressed (absent) or that the competing interaction is
weakened. As no pressure effect is found in the normal state, $c$-axis
pressure is inert to charge order itself. On the other hand, superconductivity
is known to be weakened by increasing compressive $c$-axis pressure [34]. One
commonly accepted explanation is that compressive $c$-axis pressure reduces
the apical oxygen distance to the CuO2 layers [18]. This change of crystal
field environment, in turn, boosts the hybridization of $d_{x^{2}-y^{2}}$ and
$d_{z^{2}}$ orbitals [35]. Theoretically, such hybridization is expected to be
unfavorable for $d$-wave superconductivity [36, 37]. Therefore, $c$-axis
pressure diminishes the superconducting order parameter. However, there is no
evidence of superconductivity being completely suppressed upon moderate
pressure along the $c$-axis [19, 20]. Superconductivity and charge order are
therefore expected to coexist. One possible explanation for the diminished
magnetic field effect is that $c$-axis strain induces a state in which
superconductivity rather cooperates than competes with charge order. This may
pave the way for pair-density-wave states in which superconductivity and
charge order couple in a cooperative fashion [13]. The superconducting order
parameter would have a spatial modulation commensurate to the charge order.
Search for such a superconducting modulation with $c$-axis pressure is an
interesting experimental continuation.

Figure 4: Magnetic field and strain effects. a The magnetic field scale
$H_{c}$ versus compressive $c$-axis strain at 10 K. Vertical error bars are
Gaussian standard deviations from fits to the magnetic field dependence of the
charge order reflection. Horizontal error bars reflect uncertainty of the
strain estimation. The black dashed line is a linear guide to the eye b
Schematic illustration of the magnetic field and strain effect on the charge
order peak amplitude inside the superconducting state. In the grey region, we
expect – by extrapolation – the magnetic field enhancement to fully vanish.
Crosses display the phase space covered by our experiments.

##  III Methods

###  III.1 Crystal Growth and Characterization

La2-xSrxCuO4 crystals were synthesised by the traveling-zone method. We used
samples from three different batches – all with $x\approx 0.12$. All x-ray
diffraction data in the main text stem from a batch with $x=0.12$ and
$T_{c}=27$ K. The neutron experiment was carried out on a batch with $x=0.115$
and $T_{c}=27$ K. X-ray data in the Supplementary Information were recorded on
a $x=0.125$ crystal with $T_{c}=30$ K. Crystal orientations were determined by
x-ray Laue diffraction. In this fashion, samples were cut to enable
application of $c$-axis uniaxial pressure in combination with x-ray
transmission.

###  III.2 X-Ray Diffraction and Pressure Cell

Hard x-ray diffraction experiments were carried out at the P21.1 beamline at
PETRA III (DESY, Hamburg) synchrotron using 101.6 keV photons. All
measurements were performed with a Dectris Pilatus 100K CdTe detector.
Uniaxial pressure and magnetic field were applied along the crystallographic
$c$-axis. In the absence of a magnetic field, a standard (ex-situ) screw cell
was used [25, 38]. An adaption was made to fit into a horizontal magnet [24,
28] with pressure application parallel to the direction of the magnetic field.
The designs of two uniaxial pressure cells are shown in Supplementary Fig. 4a,
b. For both cells, the horizontal scattering plane is spanned by a copper-
oxygen bond direction and the $c$-axis. This gives access to the $(h,0,\ell)$
reciprocal plane. Throughout the manuscript, reciprocal space $(h,k,\ell)$ are
given in units of $(\frac{2\pi}{a},\frac{2\pi}{b},\frac{2\pi}{c})$ with
$a=b=3.78$ Å and $c\approx 13.2$ Å. The magnitude of applied strain is
determined from measurements of Bragg reflections, similar to what is
described in Ref. [25]. We defined the zero-strain, low-temperature lattice
parameter as $c_{0}$ and $c$-axis strain as $\epsilon_{c}=(c_{0}-c)/c_{0}$.
Both $c_{0}$ and $c$ are determined from Bragg peaks of the type $(0,0,\pm
2n)$ with $n$ being an integer (see Fig. 1b). Using this methodology a
calibration curve between screw turn and strain is established (See
Supplementary Fig. 1).

###  III.3 Neutron Diffraction

Neutron diffraction data were measured at the TASP spectrometer at SINQ, PSI
on LSCO $x=0.115$ ($T_{c}\approx 27$ K). The sample in the (ex-situ) pressure
cell was mounted in a coil (see Supplementary Fig. 4c), custom-built to fit
inside a cryostat, which itself was mounted inside the Mu-metal polarization
analysis device MuPAD [39]. The sample was cooled from 40 K to 10 K in a
magnetic field of 30 G, applied along the crystallographic $c$-direction,
perpendicular to the scattering plane. At base temperature, the field was
subsequently turned off. Using MuPAD’s spherical polarimetry capabilities, a
full polarization matrix was measured on the in-plane Bragg peak
$\boldsymbol{\tau}=(-1,-1,0)$. The highest flipping ratio of 6.87 was found
for neutron polarization along $y$ (in the scattering plane but perpendicular
to $\boldsymbol{\tau}$). The temperature dependence of the intensity of
$\boldsymbol{\tau}=(-1,-1,0)$ was tracked in the Py-y spin flip channel upon
slowly heating. The flux trapped along the $c$-direction during the field-
cooling procedure gives rise to an excess depolarization of the neutron beam,
which is reflected in an increased spin-flip intensity. Upon crossing $T_{c}$,
the trapped flux is released. Intensities are normalized to the monitor and
scaled by the integrated Bragg peaks of ($-1$, $-$1, 0) in unstrained LSCO. To
extract the transition temperature, spin flip intensities were fitted with
$p_{0}\cdot\text{tanh}(T-T_{c})-p_{1}$.

Acknowledgements: J.K., R.F., L.M., D.B., and J.C. acknowledge support from
the Swiss National Science Foundation. J.K. is further supported by the PhD
fellowship from the German Academic Scholarship Foundation. I.B. acknowledges
support from the Swiss Government Excellence Scholarship. Q.W. is supported by
the Research Grants Council of Hong Kong (ECS No. 24306223), and the CUHK
Direct Grant (No. 4053613). N.B.C was supported by the Danish National Council
for Research infrastructure (NUFI) through DANSCATT and the ESS-Lighthouse
Q-MAT. This research was carried out at beamline P21.1 at DESY, a member of
the Helmholtz Association (HGF). We would like to thank Philipp Glaevecke for
technical assistance during the experiment. The research leading to this
result has been supported by the project CALIPSOplus under the Grant Agreement
730872 from the EU Framework Programme for Research and Innovation HORIZON
2020. Part of the research have been carried out at the TASP endstation of the
spallation source SINQ (PSI).

Authors Contribution: TK, NM, MO grew and JK prepared the LSCO crystals.
JChoi, DB, MJ designed and tested the pressure cells. X-ray diffraction
experiments were carried out by JK, RF, AM, OI, MvZ, NBC, QW, JChang.
Polarized neutron scattering experiments were conducted by JK, NBC, DGM, GS,
AAT, LT, DWT. JK carried out the data analysis with input from IB, QW, and LM.
The project was conceived by JK and JChang who also wrote the manuscript with
input from all authors.

## References

  * Gibbs  _et al._ [1988] D. Gibbs, K. M. Mohanty, and J. Bohr, Phys. Rev. B 37, 562 (1988).
  * Jiang and Fishman [1997] X. W. Jiang and R. S. Fishman, Journal of Physics: Condensed Matter 9, 3417 (1997).
  * Hu  _et al._ [2022] Y. Hu, T. Zhang, D. Zhao, C. Chen, S. Ding, W. Yang, X. Wang, C. Li, H. Wang, D. Feng, and T. Zhang, Nature Communications 13, 445 (2022).
  * Tranquada  _et al._ [1994] J. M. Tranquada, D. J. Buttrey, V. Sachan, and J. E. Lorenzo, Phys. Rev. Lett. 73, 1003 (1994).
  * Tranquada  _et al._ [1996] J. M. Tranquada, J. D. Axe, N. Ichikawa, Y. Nakamura, S. Uchida, and B. Nachumi, Phys. Rev. B 54, 7489 (1996).
  * Hücker  _et al._ [2011] M. Hücker, M. v. Zimmermann, G. D. Gu, Z. J. Xu, J. S. Wen, G. Xu, H. J. Kang, A. Zheludev, and J. M. Tranquada, Phys. Rev. B 83, 104506 (2011).
  * Thampy  _et al._ [2014] V. Thampy, M. P. M. Dean, N. B. Christensen, L. Steinke, Z. Islam, M. Oda, M. Ido, N. Momono, S. B. Wilkins, and J. P. Hill, Phys. Rev. B 90, 100510(R) (2014).
  * Achkar  _et al._ [2016] A. J. Achkar, M. Zwiebler, C. McMahon, F. He, R. Sutarto, I. Djianto, Z. Hao, M. J. P. Gingras, M. Hücker, G. D. Gu, A. Revcolevschi, H. Zhang, Y.-J. Kim, J. Geck, and D. G. Hawthorn, Science 351, 576 (2016).
  * Yu  _et al._ [2021] F. H. Yu, D. H. Ma, W. Z. Zhuo, S. Q. Liu, X. K. Wen, B. Lei, J. J. Ying, and X. H. Chen, Nature Communications 12, 3645 (2021).
  * Song  _et al._ [2021] Y. Song, T. Ying, X. Chen, X. Han, X. Wu, A. P. Schnyder, Y. Huang, J.-g. Guo, and X. Chen, Phys. Rev. Lett. 127, 237001 (2021).
  * Nandi  _et al._ [2010] S. Nandi, M. G. Kim, A. Kreyssig, R. M. Fernandes, D. K. Pratt, A. Thaler, N. Ni, S. L. Bud’ko, P. C. Canfield, J. Schmalian, R. J. McQueeney, and A. I. Goldman, Phys. Rev. Lett. 104, 057006 (2010).
  * Hu  _et al._ [2020] D. Hu, D. W. Tam, W. Zhang, Y. Wei, R. Georgii, B. Pedersen, A. C. Roldan, and P. Dai, Phys. Rev. B 101, 020507 (2020).
  * Agterberg  _et al._ [2020] D. F. Agterberg, J. S. Davis, S. D. Edkins, E. Fradkin, D. J. V. Harlingen, S. A. Kivelson, P. A. Lee, L. Radzihovsky, J. M. Tranquada, and Y. Wang, Annual Review of Condensed Matter Physics 11, 231 (2020).
  * Wang  _et al._ [2018] Y. Wang, S. D. Edkins, M. H. Hamidian, J. C. S. Davis, E. Fradkin, and S. A. Kivelson, Phys. Rev. B 97, 174510 (2018).
  * Dai  _et al._ [2018] Z. Dai, Y.-H. Zhang, T. Senthil, and P. A. Lee, Phys. Rev. B 97, 174511 (2018).
  * Tranquada  _et al._ [1995] J. M. Tranquada, B. J. Sternlieb, J. D. Axe, Y. Nakamura, and S. Uchida, Nature (London) 375, 561 (1995).
  * Croft  _et al._ [2014] T. P. Croft, C. Lester, M. S. Senn, A. Bombardi, and S. M. Hayden, Phys. Rev. B 89, 224513 (2014).
  * Kramer  _et al._ [2019] K. P. Kramer, M. Horio, S. S. Tsirkin, Y. Sassa, K. Hauser, C. E. Matt, D. Sutter, A. Chikina, N. B. M. Schröter, J. A. Krieger, T. Schmitt, V. N. Strocov, N. C. Plumb, M. Shi, S. Pyon, T. Takayama, H. Takagi, T. Adachi, T. Ohgi, T. Kawamata, Y. Koike, T. Kondo, O. J. Lipscombe, S. M. Hayden, M. Ishikado, H. Eisaki, T. Neupert, and J. Chang, Phys. Rev. B 99, 224509 (2019).
  * Pavarini  _et al._ [2001] E. Pavarini, I. Dasgupta, T. Saha-Dasgupta, O. Jepsen, and O. K. Andersen, Phys. Rev. Lett. 87, 047003 (2001).
  * Nakamura  _et al._ [1999] F. Nakamura, J. Hori, T. Goko, Y. Uno, N. Kikugawa, and T. Fujita, Journal of Low Temperature Physics 117, 1145 (1999).
  * Sato and Naito [1997] H. Sato and M. Naito, Physica C: Superconductivity 274, 221 (1997).
  * Locquet and Williams [1997] J. Locquet and E. Williams, Acta Phys. Pol. A 92, 69–84 (1997).
  * Takeshita  _et al._ [2004] N. Takeshita, T. Sasagawa, T. Sugioka, Y. Tokura, and H. Takagi, Journal of the Physical Society of Japan 73, 1123 (2004).
  * Christensen  _et al._ [2014] N. B. Christensen, J. Chang, J. Larsen, M. Fujita, M. Oda, M. Ido, N. Momono, E. M. Forgan, A. T. Holmes, J. Mesot, M. Huecker, and M. v. Zimmermann, arXiv:1404.3192 (2014).
  * Choi  _et al._ [2022] J. Choi, Q. Wang, S. Jöhr, N. B. Christensen, J. Küspert, D. Bucher, D. Biscette, M. H. Fischer, M. Hücker, T. Kurosawa, N. Momono, M. Oda, O. Ivashko, M. v. Zimmermann, M. Janoschek, and J. Chang, Phys. Rev. Lett. 128, 207002 (2022).
  * Chang  _et al._ [2008] J. Chang, C. Niedermayer, R. Gilardi, N. B. Christensen, H. M. Rønnow, D. F. McMorrow, M. Ay, J. Stahn, O. Sobolev, A. Hiess, S. Pailhes, C. Baines, N. Momono, M. Oda, M. Ido, and J. Mesot, Phys. Rev. B 78, 104525 (2008).
  * Lake  _et al._ [2002] B. Lake, H. M. Rønnow, N. B. Christensen, G. Aeppli, K. Lefmann, D. F. McMorrow, P. Vorderwisch, P. Smeibidl, N. Mangkorntong, T. Sasagawa, M. Nohara, H. Takagi, and T. E. Mason, Nature 415, 299 (2002).
  * Hücker  _et al._ [2013] M. Hücker, M. v. Zimmermann, Z. J. Xu, J. S. Wen, G. D. Gu, and J. M. Tranquada, Phys. Rev. B 87, 014501 (2013).
  * Khaykovich  _et al._ [2005] B. Khaykovich, S. Wakimoto, R. J. Birgeneau, M. A. Kastner, Y. S. Lee, P. Smeibidl, P. Vorderwisch, and K. Yamada, Phys. Rev. B 71, 220508 (2005).
  * Wen  _et al._ [2023] J.-J. Wen, W. He, H. Jang, H. Nojiri, S. Matsuzawa, S. Song, M. Chollet, D. Zhu, Y.-J. Liu, M. Fujita, J. M. Jiang, C. R. Rotundu, C.-C. Kao, H.-C. Jiang, J.-S. Lee, and Y. S. Lee, Nature Communications 14, 733 (2023).
  * Demler  _et al._ [2001] E. Demler, S. Sachdev, and Y. Zhang, Phys. Rev. Lett. 87, 067202 (2001).
  * Lake  _et al._ [2001] B. Lake, G. Aeppli, K. N. Clausen, D. F. McMorrow, K. Lefmann, N. E. Hussey, N. Mangkorntong, M. Nohara, H. Takagi, T. E. Mason, and A. Schröder, Science 291, 1759 (2001).
  * Wu  _et al._ [2013] T. Wu, H. Mayaffre, S. Krämer, M. Horvatić, C. Berthier, P. L. Kuhns, A. P. Reyes, R. Liang, W. N. Hardy, D. A. Bonn, and M.-H. Julien, Nature Communications 4, 2113 (2013).
  * Nakamura  _et al._ [2000] F. Nakamura, T. Goko, J. Hori, Y. Uno, N. Kikugawa, and T. Fujita, Phys. Rev. B 61, 107 (2000).
  * Matt  _et al._ [2018] C. E. Matt, D. Sutter, A. M. Cook, Y. Sassa, M. Månsson, O. Tjernberg, L. Das, M. Horio, D. Destraz, C. G. Fatuzzo, K. Hauser, M. Shi, M. Kobayashi, V. N. Strocov, T. Schmitt, P. Dudin, M. Hoesch, S. Pyon, T. Takayama, H. Takagi, O. J. Lipscombe, S. M. Hayden, T. Kurosawa, N. Momono, M. Oda, T. Neupert, and J. Chang, Nat. Commun. 9, 972 (2018).
  * Sakakibara  _et al._ [2010] H. Sakakibara, H. Usui, K. Kuroki, R. Arita, and H. Aoki, Phys. Rev. Lett. 105, 057003 (2010).
  * Sakakibara  _et al._ [2012] H. Sakakibara, H. Usui, K. Kuroki, R. Arita, and H. Aoki, Phys. Rev. B 85, 064501 (2012).
  * Simutis  _et al._ [2022] G. Simutis, J. Küspert, Q. Wang, J. Choi, D. Bucher, M. Boehm, F. Bourdarot, M. Bertelsen, C. N. Wang, T. Kurosawa, N. Momono, M. Oda, M. Månsson, Y. Sassa, M. Janoschek, N. B. Christensen, J. Chang, and D. G. Mazzone, Communications Physics 5, 296 (2022).
  * Janoschek  _et al._ [2007] M. Janoschek, S. Klimko, R. Gähler, B. Roessli, and P. Böni, Physica B: Condensed Matter 397, 125 (2007).

Supplementary Material for
Engineering Phase Competition Between Stripe Order and Superconductivity in
La1.88Sr0.12CuO4

Supplementary Figure 1: Strain calibration curve. Pressure cell screw turns
($st$) versus strain $\epsilon_{c}=(c-c_{0})/c_{0}$, where $c_{0}$ is the
unstrained c-axis lattice parameter and $c$ is the strained lattice parameter.
$c_{0}$ and $c$ are determined as described in the Methods section. The solid
line is a linear fit of the form $st=m\cdot\epsilon_{c}$, with $m$ being a fit
parameter. Supplementary Figure 2: Example of detector image. 2D detector
image recorded on LSCO $x=12$ % ($T_{c}=27$ K) and displayed with false color
scale. Signal region-of-interest (ROI) around $\bf{Q}_{co}\rm=(0.231,0,12.5)$
and background ROIs are indicated by respectively blue and grey rectangles.
The most intense feature running horizontally across the detector is a powder
line stemming from the sample environment. The vertical line is a gap between
two parts of the detector. Supplementary Figure 3: Charge order reflection
recorded on LSCO ($\bf{T_{c}=30}$ K, $\bf{x=12.5}$ %) upon application of
$c$-axis uniaxial strain. a-c Background subtracted intensity along
$\bf{Q}_{co}\rm$ $=(h,0,12.5)$ for temperature and strain conditions as
indicated. Horizontal and vertical error bars represent respectively
systematic errors and counting statistics. Horizontal error bars are indicated
only on selected points.  Supplementary Figure 4: Strain cells for $c$-strain
application. a Strain cell for experimental setup without magnetic field
(displex cryo-environment).We extended the strain cell used in [25] with an
adapter to apply strain in $c$-axis direction while having $a$\- and $c$-axis
in the horizontal scattering plane. b Strain cell for experiment with a cryo-
magnet. The design is more compact in order to fit into the magnet bore. c
Strain cell from a (without adapter) inside the coil to be mounted inside the
MuPad [39] for the neutron experiment. Supplementary Figure 5: Comparison of
fit functions for the charge order peak at $\bf{Q}_{co}=(0.231,0,12.5)$. a
Intensity from rocking scans through the charge order reflection and
corresponding background obtained from integration of ”signal” and
”background” ROIs (see Supplementary Fig. 1). The sample was exposed to an
uniaxial strain of $\epsilon_{c}=(0.17\pm 0.01)$ %, measured by the
experimental setup using a cryo-magnet at 30 K. b Background subtracted
intensities for the rocking scan shown in a. Solid black and red lines are
respectively fits to a Gaussian function and a split-normal distribution. The
resulting amplitudes are $(182.99\pm 6.00)$ counts/s and $(184.94\pm 5.95)$,
respectively. Thus within error bars, the amplitude is independent of the two
fit functions. As the split-normal distribution provides a better fit, this
function is consistently used in our analysis. To be specific, a split-normal
distribution with linear background is defined as:
$I(\omega)=I_{co}\cdot\exp(-(\omega-\omega_{0})^{2}/(2\sigma^{2}))+b_{1}\omega+b_{0}$,
where $\sigma=\sigma_{1}$ if $\omega\leq\omega_{0}$ and $\sigma=\sigma_{2}$ if
$\omega>\omega_{0}$. $\sigma$ is the width of the function, $\omega_{0}$ its
peak position, $I_{co}$ the amplitude, and $b_{1}$ and $b_{0}$ fitting
parameters defining the background.


Duy Nam Bui and Manh Duong Phung

# Radial Basis Function Neural Networks for Formation Control of Unmanned
Aerial Vehicles

Duy-Nam Bui  Vietnam National University, Hanoi, Vietnam  Manh Duong Phung*
Fulbright University Vietnam, Ho Chi Minh City, Vietnam
<EMAIL_ADDRESS>*Corresponding author.

###### Abstract

This paper addresses the problem of controlling multiple unmanned aerial
vehicles (UAVs) cooperating in a formation to carry out a complex task such as
surface inspection. We first use the virtual leader-follower model to
determine the topology and trajectory of the formation. A double-loop control
system combining backstepping and sliding mode control techniques is then
designed for the UAVs to track the trajectory. A radial basis function neural
network (RBFNN) capable of estimating external disturbances is developed to
enhance the robustness of the controller. The stability of the controller is
proven by using the Lyapunov theorem. A number of comparisons and software-in-
the-loop (SIL) tests have been conducted to evaluate the performance of the
proposed controller. The results show that our controller not only outperforms
other state-of-the-art controllers but is also sufficient for complex tasks of
UAVs such as collecting surface data for inspection. The source code of our
controller can be found at https://github.com/duynamrcv/rbf_bsmc.

###### keywords:

Unmanned aerial vehicles (UAVs), formation control, radial basis function
neural network, backstepping sliding mode control

††articletype: RESEARCH ARTICLE

##  1 Introduction

Unmanned aerial vehicles, when combined with computer vision technologies, can
collect visual data of structures to provide valuable information for various
tasks such as inspecting structural surfaces [1, 2], reconstructing 3D models
[3, 4], identifying cracks [5, 6], and detecting corrosion and rust on steel
bridges [7, 8]. However, using a single UAV for these tasks is inefficient due
to the large size of the structures and the limited battery capacity of the
UAV. A group of UAVs flying in a formation can be used to overcome those
limitations [9, 10, 11, 12]. The formation allows the UAVs to perform
collaborative inspection to increase the efficiency and accuracy of data
collection. It also allows for safe operation, as the formation control can
prevent collision among the UAVs.

In formation control, the leader-follower approach is commonly used to provide
flexibility in topology and trajectory selection[13, 14]. In the standard
leader-follower method, one UAV is assigned as the leader, and the others are
followers. The leader plays the role of a reference node for the followers to
determine their locations to form the desired topology. The limitation of this
approach, however, is the dependence of the system on the leader. If the
leader is malfunctioning, the whole system will fail. The virtual leader-
follower model can be used to cope with this problem. In this approach, the
leader is purely a virtual entity, serving as a reference point for the
followers to determine their positions [15]. By decoupling the physical leader
from the model, this method mitigates the risk of complete system failure.

In the leader-follower model, linear controllers are commonly used to control
individual UAVs to form the desired topology [15, 16, 17]. In [16], linear
quadratic and neural networks-based controllers are combined to control a
group of UAVs considering their full dynamics. In [17], the receding horizon
control is employed to yield a fast convergence rate of the formation tracking
control. This controller also considers the orientation between the leader and
the followers for accurate formation. The decentralized $H_{\infty}$-PID
controller is introduced in [18] to maneuver a group of UAVs to deal with the
external disturbance and trailing vortex coupling from their neighbor UAVs. A
leader-follower formation control technique is presented in [19] to address
issues related to backward error and suboptimal dynamic speed tracking in PID
neural network control. Linear controllers, however, have limitations in
handling constraints and parameter variation, especially when applied to
nonlinear systems like UAVs.

In another approach, nonlinear controllers have been used for formation
control [20, 21, 22]. In [23], first and second-order sliding-mode controllers
are deployed to assure the asymptotic stability of the formation, taking into
account modeling uncertainties. In [24], an adaptive controller using the
dynamic estimation of the distance between the leader and the followers is
introduced to address uncertainties related to positioning errors. Two finite-
time observers are used in [25] to deal with bounded external disturbance
force and torque. In [26], a non-uniform vector field that dynamically varies
in magnitude and direction is employed to deal with the influence of wind in
UAV formation control. A distributed model predictive control algorithm is
introduced in [27] to coordinate the operation of a fleet of UAVs considering
their spatial kinematics and unidirectional data transmissions. However, the
convergence of these controllers depends on the characteristics of
disturbances, which are hard to model due to their varying nature. A
sufficient approach would be utilizing neural networks such as the radial
basis function neural network (RBFNN) to estimate disturbances and use it as
the feedback for control [28, 29, 30, 31].

In this work, we present a new controller for a group of UAVs cooperating in a
formation. The UAVs use the virtual leader-follower model to determine their
trajectory and form the desired topology. The controller is developed using
the backstepping and sliding mode control techniques. An RBFNN is then
introduced to estimate external disturbances for better control performance.
Our contributions in this work are as follows:

  1. i.

The proposal of a new controller for UAV formation that is constructed by
combining backstepping and sliding mode control techniques, thereby enabling
the elimination of nonlinear components and enhancing system robustness.
Additionally, the adverse effects associated with these controllers, such as
“explosion of term” and “chattering”, are mitigated through the approximation
of unknown factors by the neural network. As the result, the developed
controller not only addresses the drawbacks of the aforementioned techniques
but also augments the adaptability of the UAV system.

  2. ii.

The design of a radial basis function neural network (RBFNN) that is capable
of estimating external disturbances to compensate for input force control
signals, thereby enabling the controller to maintain the required control
quality.

  3. iii.

The derivation of the stability proof for the designed controller using
Lyapunov’s theorem, which is essential to ensure stable operation of the UAVs
under conditions affected by external forces.

  4. iv.

The comparison of the proposed controller with other popular methods including
model predictive control (MPC), backstepping sliding mode control (BSMC) and
sliding mode control (SMC) in different scenarios to confirm its superior
performance. Software-in-the-loop tests were also conducted with a cooperative
bridge inspection task to verify the validity of the proposed method for
practical applications.

The rest of this paper is structured as follows. Section 2 presents the
dynamic and formation models of the UAVs. Section 3 introduces the proposed
controller. Section 4 shows evaluation results. The paper ends with
conclusions described in 5.

##  2 Problem formulation

Figure 1: The structure of the quadrotor UAV in the global frame

To control a group of UAVs, we first consider their dynamic model and
formation topology with details as follows.

###  2.1 UAV dynamic model

Consider a group of $n$ UAVs, each is a quadrotor with two pairs of propellers
rotating in opposite directions, as described in Figure 1. Frames
$BB_{x}B_{y}B_{z}$ and $OXYZ$ are respectively the body-fixed and inertial
frames. We use Euler angles to represent the attitude of the UAV. The
configuration of the UAV includes its position $\xi=\left[x,y,z\right]^{T}$
and Euler angles $\Xi=\left[\phi,\theta,\psi\right]^{T}$, with
$\left|\phi\right|\leq\pi/2$, $\left|\theta\right|\leq\pi/2$ and
$\left|\psi\right|\leq\pi$. Those angles represent the roll, pitch and yaw
orientation of the UAV, respectively. Control signals of the UAV are defined
as follows:

$\left[\begin{array}[]{c}f_{t}\\\ \tau_{\phi}\\\ \tau_{\theta}\\\
\tau_{\psi}\end{array}\right]=\left[\begin{array}[]{c}f_{1}+f_{2}+f_{3}+f_{4}\\\
l\left(f_{4}-f_{2}\right)\\\ l\left(f_{3}-f_{1}\right)\\\
\tau_{2}+\tau_{4}-\tau_{1}-\tau_{3}\end{array}\right],$ (1)

where $l$ is the arm length; $f_{t}$ is the total thrust of four propellers;
$\tau_{\phi}$, $\tau_{\theta}$, $\tau_{\psi}$ are the torques in three axes;
and $f_{i}$ and $\tau_{i}$, with $i=\left\\{1,2,3,4\right\\}$, are the forces
and torques generated by four propellers, respectively. According to [32], the
dynamic model of the UAV is described as follows:

$\displaystyle\ddot{x}$
$\displaystyle=\left(\cos\phi\sin\theta\cos\psi+\sin\phi\sin\psi\right)\dfrac{f_{t}}{m}+\dfrac{d_{x}}{m}$
(2) $\displaystyle\ddot{y}$
$\displaystyle=\left(\cos\phi\sin\theta\sin\psi-\sin\phi\cos\psi\right)\dfrac{f_{t}}{m}+\dfrac{d_{y}}{m}$
$\displaystyle\ddot{z}$
$\displaystyle=\cos\phi\cos\theta\dfrac{f_{t}}{m}-g+\dfrac{d_{z}}{m}$
$\displaystyle\ddot{\phi}$
$\displaystyle=\dfrac{\dot{\theta}\dot{\psi}\left(I_{y}-I_{z}\right)+\tau_{\phi}}{I_{x}}$
$\displaystyle\ddot{\theta}$
$\displaystyle=\dfrac{\dot{\phi}\dot{\psi}\left(I_{z}-I_{x}\right)+\tau_{\theta}}{I_{y}}$
$\displaystyle\ddot{\psi}$
$\displaystyle=\dfrac{\dot{\phi}\dot{\theta}\left(I_{x}-I_{y}\right)+\tau_{\psi}}{I_{z}}$

where $I_{x}$, $I_{y}$, $I_{z}$ are the moments of inertia, $m$ is the mass of
the UAV, $g$ is the gravitational acceleration, and
$\left[d_{x},d_{y},d_{z}\right]^{T}$ is the disturbance caused by factors such
as wind or turbulent flows.

###  2.2 UAV formation model

The formation model used in this work is the virtual leader-follower model
with two main components:

  * •

Virtual leader: a virtual leader is a non-physical UAV used as a reference for
other UAVs to determine their position. Its trajectory represents the
trajectory of the UAV group.

  * •

Follower: a follower is a UAV that adjusts its position based on the virtual
leader. Given the reference trajectory of the leader and the expected
topology, the followers calculate their trajectories and then track them to
form the desired formation.

Figure 2: Illustration of the virtual leader-follower formation structure

Consider virtual leader $B_{L}$ having position
$\xi_{L}=\left[x_{L},y_{L},z_{L}\right]^{T}$ and heading angle $\psi_{L}$ and
follower $B_{i}$ having position $\xi_{L}=\left[x_{i},y_{i},z_{i}\right]^{T}$
and yaw angle $\psi_{i}$. Let $B^{\prime}_{L}$ and $B^{\prime}_{i}$ be their
projection on the $OXY$ plane, respectively, $B^{\prime\prime}_{L}$ be the
projection of $B_{i}$ on $B_{L}B^{\prime}_{L}$, and $\Delta_{i}=\left[\Delta
x_{i},\Delta y_{i},\Delta z_{i}\right]^{T}$ be the desired distance between
follower $B_{i}$ and the virtual leader, as depicted in Figure 2. Since
$\Delta z_{i}=B_{L}B^{\prime\prime}_{L}$, the desired position of follower
$B_{i}$ can be computed as:

$\displaystyle{}^{d}\xi_{i}$
$\displaystyle=\text{Rot}_{z}\left(\psi_{L}\right)\Delta_{i}+\xi_{L}$ (3)
$\displaystyle{}^{d}\psi_{i}$ $\displaystyle=\psi_{L},$

where $\text{Rot}_{z}(\cdot)\in\mathbb{R}^{3\times 3}$ is the rotation matrix
around z-axis. Equation (3) allows the followers to compute their trajectory
based on the trajectory of the virtual leader and the desired formation
topology.

##  3 Controller design for UAV formation

Figure 3: The proposed controller

Given the trajectory of the virtual leader, denoted as $(\xi_{L},\psi_{L})$,
the desired trajectory of follower $i$, $({}^{d}\xi_{i},{}^{d}\psi_{i})$, in
the formation can be computed based on (3). To track this trajectory, we
design a dual-loop control system for each follower as shown in Figure 3. The
outer loop is a position controller that regulates the altitude and horizontal
position, while the inner loop is a backstepping sliding mode controller
(BSMC) that handles the UAV’s attitude, including its roll, pitch, and yaw
angles. To account for external disturbances, the position controller was
designed with a radial basis function neural network (RBFNN). A converter
block is also included to convert the desired translational control forces
into roll and pitch angles. Details of each controller are described as
follows.

###  3.1 Position controller design

The position controller aims to keep the UAV’s position aligned with the
desired trajectory. It is designed based on BSMC with the use of RBFNN for
disturbance estimation. According to (2), dynamic equations for the
translational motion of UAV $i$ are given as follows:

$\displaystyle\ddot{x}_{i}$ $\displaystyle=u_{xi}+\dfrac{d_{xi}}{m},$ (4)
$\displaystyle\ddot{y}_{i}$ $\displaystyle=u_{yi}+\dfrac{d_{yi}}{m},$
$\displaystyle\ddot{z}_{i}$ $\displaystyle=u_{zi}+\dfrac{d_{zi}}{m},$

where

$\displaystyle u_{xi}$
$\displaystyle=\left(\cos\phi_{i}\sin\theta_{i}\cos\psi_{i}+\sin\phi_{i}\sin\psi_{i}\right)\dfrac{f_{ti}}{m},$
(5) $\displaystyle u_{yi}$
$\displaystyle=\left(\cos\phi_{i}\sin\theta_{i}\sin\psi_{i}-\sin\phi_{i}\cos\psi_{i}\right)\dfrac{f_{ti}}{m},$
$\displaystyle u_{zi}$
$\displaystyle=\cos\phi_{i}\cos\theta_{i}\dfrac{f_{ti}}{m}-g.$

Let $\xi_{i}=\left[x_{i},y_{i},z_{i}\right]^{T}$ and
$\dot{\xi}_{i}=\left[\dot{x}_{i},\dot{y}_{i},\dot{z}_{i}\right]^{T}$
respectively be the position and velocity of the translational motion,
$U_{i}=\left[u_{xi},u_{yi},u_{zi}\right]^{T}$ be the control signal, and
$D_{i}=\left[\dfrac{d_{xi}}{m},\dfrac{d_{yi}}{m},\dfrac{d_{zi}}{m}\right]^{T}$,
with $\left\|D_{i}\right\|\leq\bar{d}$, be the external disturbance affecting
UAV $i$. Equation (4) can be rewritten as:

$\ddot{\xi}_{i}=U_{i}+D_{i}.$ (6)

The BSMC is then designed as follows.

####  3.1.1 Backstepping sliding mode controller (BSMC) design

Let ${}^{e}\xi_{i}$ be the translational error,
${}^{e}\xi_{i}=\xi_{i}-{}^{d}\xi_{i}$. The virtual velocity ,${}^{v}\xi_{i}$,
of the subsystem is designed as:

${}^{v}\xi_{i}={}^{d}\dot{\xi}_{i}-\lambda_{\xi}{}^{e}\xi_{i},$ (7)

where $\lambda_{\xi}>0$ is a positive definite gain. The first candidate
Lyapunov function is chosen as

${}^{1}V_{\xi_{i}}=\dfrac{1}{2}{}^{e}\xi_{i}^{T}{}^{e}\xi_{i}.$ (8)

Its derivative is given by

${}^{1}\dot{V}_{\xi_{i}}={}^{e}\xi_{i}^{T}{}^{e}\dot{\xi}_{i}={}^{e}\xi_{i}^{T}\left(\dot{\xi}_{i}-{}^{d}\dot{\xi}_{i}\right).$
(9)

Substituting $\dot{\xi}_{i}={}^{v}{\xi}_{i}$ into (9) gives

${}^{1}\dot{V}_{\xi_{i}}=-\lambda_{\xi}\left\|{}^{e}\xi_{i}\right\|^{2}\leq
0.$ (10)

Hence, the system is stable with the virtual velocity chosen in (7). The
sliding mode control (SMC) algorithm is then utilized to design the input
control signal for the position system. The sliding surface is chosen as
follows:

$s_{\xi_{i}}=\gamma_{\xi}{}^{e}\xi_{i}+\left(\dot{\xi}_{i}-{}^{v}{\xi}_{i}\right),$
(11)

where $\gamma_{\xi}>0$ is a positive definite gain. Denote $\hat{D}_{i}$ as
the disturbance estimated via an estimator such as the RBFNN in Section 3.1.2.
The derivative of $s_{\xi_{i}}$ then can be obtained by using $\hat{D}_{i}$
instead of $D_{i}$ as follows:

$\dot{s}_{\xi_{i}}=\gamma_{\xi}\left(\dot{\xi}_{i}-{}^{d}\dot{\xi}_{i}\right)+U_{i}+\hat{D}_{i}-{}^{v}\dot{\xi}_{i}.$
(12)

The second Lyapunov function of the subsystem is chosen as follows:

${}^{2}V_{\xi_{i}}=\dfrac{1}{2}s_{\xi_{i}}^{T}s_{\xi_{i}}$ (13)

The control signals are designed as follows:

$\displaystyle U_{ieq}$
$\displaystyle={}^{v}\dot{\xi}_{i}-\gamma_{\xi}\left(\dot{\xi}_{i}-{}^{d}\dot{\xi}_{i}\right)-\hat{D}_{i}$
(14) $\displaystyle U_{isw}$ $\displaystyle=-\left(c_{\xi
1}\text{sg}\left(s_{\xi_{i}}\right)+c_{\xi 2}s_{\xi_{i}}\right),$

where $c_{\xi 1}$ and $c_{\xi 2}$ are positive gains, $U_{ieq}$ is the
equivalent control signal that maintains the position variables on the sliding
manifold, $U_{isw}$ is the signal that leads the subsystem to the sliding
surface $s_{\xi_{i}}$, and $\text{sg}\left(\cdot\right)$ is the piece-wise
continuous function defined as

$\text{sg}\left(x\right)=\left\\{\begin{array}[]{c}1\\\ -1\\\
\dfrac{x}{\epsilon}\end{array}\right.\begin{array}[]{c}x>\epsilon\\\
x<-\epsilon\\\ \text{otherwise}\end{array}$ (15)

where $0<\epsilon<1$ is a pre-defined constant.

######  Theorem 3.1.

Consider the position control system of the UAV. If the control signal is
chosen as

$U_{i}=U_{ieq}+U_{isw},$ (16)

the system is stable.

######  Proof 3.2.

Taking the first derivative of ${}^{2}V_{\xi_{i}}$ gives

${}^{2}\dot{V}_{\xi_{i}}=s_{\xi_{i}}^{T}\dot{s}_{\xi_{i}}.$ (17)

Substituting (12) into (17) gives

${}^{2}\dot{V}_{\xi_{i}}=s_{\xi_{i}}^{T}\left(\gamma_{\xi}\left(\dot{\xi}_{i}-{}^{d}\dot{\xi}_{i}\right)+U_{i}+\hat{D}_{i}-{}^{v}\dot{\xi}_{i}\right).$
(18)

By substituting (14) and (16) into (18), $\dot{V}_{\xi_{i}}$ becomes

${}^{2}\dot{V}_{\xi_{i}}=-c_{\xi
1}s_{\xi_{i}}^{T}\text{sg}\left(s_{\xi_{i}}\right)-c_{\xi
2}\left\|s_{\xi_{i}}\right\|^{2}\leq 0.$ (19)

According to Lyapunov’s stability theorem, the system is stable.

####  3.1.2 Radial basis function neural network (RBFNN) design

Figure 4: The RBFNN structure

During operation, UAVs are subject to inevitable disturbances such as wind or
turbulent flows. Those disturbances affect the system performance, but are
complex to model and analyze. We address this problem by exploiting the online
learning capability of neural networks to estimate the disturbances. Previous
studies on universal approximation theorems for RBFNN show that RBFNN can
approximate any nonlinear function on a bounded set with an arbitrary level of
accuracy [31]. In this work, we design a disturbance estimator using a neural
network with the radial basis function (RBF). The network has three layers
including an input layer, a hidden layer, and an output layer, as shown in
Figure 4. Position vector $\xi_{i}$ and its derivation $\dot{\xi}_{i}$ are the
input of the network. At the hidden layer, neurons are activated by a radial
basis function. The output of neuron $j$ is computed as:

$h_{j}=\exp{\left(-\dfrac{\left\|\xi_{i}-\mu_{1j}\right\|^{2}+\left\|\dot{\xi}_{i}-\mu_{2j}\right\|^{2}}{b^{2}}\right)},$
(20)

where $b$ is a parameter controlling the width of the Gaussian function,
$\mu_{1j}$ and $\mu_{2j}$ are predefined center points, and
$j\in\\{1,2,...,m\\}$ is the neuron index with $m$ being the number of neurons
in the hidden layer. The output layer is a weighted sum. Let $W_{i}$ be the
optimal weight matrix, $H_{i}$ be the output of the hidden layer, and
$\sigma_{i}$ be the approximation error. Disturbance $D_{i}$ affecting UAV $i$
then can be expressed by:

$D_{i}=W_{i}^{T}H_{i}+\sigma_{i}$ (21)

The output $\hat{D}_{i}$ of the RBFNN approximates $D_{i}$ as:

$\hat{D}_{i}=\hat{W}_{i}^{T}H_{i},$ (22)

where $\hat{W}_{i}$ is a trained weight matrix. This matrix is updated based
on the following rule:

$\dot{\hat{W}}_{i}=a\left(H_{i}s_{\xi_{i}}^{T}-\eta\left\|s_{\xi_{i}}\right\|\hat{W}_{i}\right),$
(23)

where $a$ is a positive definite gain matrix. With this structure, the
estimation of disturbance $D_{i}$ can be described as in Algorithm 1, where
$\beta\in(0,1)$ is the momentum factor.

/* Initialization only one time */

1 Initialize parameters $m$, $a$, $b$, $\eta$, $\beta$;

2 Initialize predefined center points $\mu_{1}$, $\mu_{2}$;

3 Create random weight matrix $\hat{W}_{i}$;

4 Initialize stored weight matrices $\hat{W}_{1i}=\hat{W}_{2i}=\hat{W}_{i}$;

5

/* Online training */

6 foreach  _computation step_ do

7 Get the current values of $\xi_{i}$, $\dot{\xi}_{i}$, $s_{\xi_{i}}$;

8 $H_{i}=\text{zeros}(m,1)$;

9 for  _$j=1$ to $m$_ do

Compute sub-hidden layer $h_{j}$;

/* Equation 20 */

10

Compute the value of $\dot{\hat{W}}_{i}$;

/* Equation 23 */

11 Update the weight matrix
$\hat{W}_{i}=\hat{W}_{1i}+\dot{\hat{W}}_{i}+\beta(\hat{W}_{1i}-\hat{W}_{2i})$;

12 Update stored weight matrices $\hat{W}_{2i}=\hat{W}_{1i}$ and
$\hat{W}_{1i}=\hat{W}_{i}$;

Compute the current estimated disturbance $\hat{D}_{i}$;

/* Equation 22 */

13

Algorithm 1 Pseudocode to estimate disturbances by RBFNN for UAV $i$

####  3.1.3 Stability of the position controller

The stability of this controller is addressed in Theorem 3.3 as follows.

######  Theorem 3.3.

Consider UAV $i$ affected by external disturbance $D_{i}$ as described in (6),
the control signals designed in (14), the bounded weight
$\left\|W_{i}\right\|\leq\bar{W}_{i}$, and the update rule for RBFNN in (23).
If the following inequality condition is satisfied:

$\left\|s_{\xi_{i}}\right\|\geq\dfrac{\bar{\sigma}+\dfrac{1}{4}\eta\bar{W}_{i}^{2}}{c_{\xi
2}},$ (24)

the position control system is stable.

######  Proof 3.4.

Choose the candidate Lyapunov function as follows:

$V_{\xi_{i}}={}^{2}V_{\xi_{i}}+\dfrac{1}{2}\text{tr}\left(\tilde{W}_{i}^{T}a^{-1}\tilde{W}_{i}\right),$
(25)

where $\tilde{W}_{i}=W_{i}-\hat{W}_{i}$ is the error weight matrix. Taking the
first derivative of $V_{\xi_{i}}$ gives:

$\displaystyle\dot{V}_{\xi_{i}}$
$\displaystyle=s_{\xi_{i}}^{T}\dot{s}_{\xi_{i}}+\text{tr}\left(\tilde{W}_{i}^{T}a^{-1}\dot{\tilde{W}_{i}}\right)$
(26) $\displaystyle=-c_{\xi
1}s_{\xi_{i}}^{T}\text{sg}\left(s_{\xi_{i}}\right)-c_{\xi
2}\left\|s_{\xi_{i}}\right\|^{2}-s_{\xi_{i}}^{T}\left(\hat{D}_{i}-D_{i}\right)-\text{tr}\left(\tilde{W}_{i}^{T}a^{-1}\dot{\hat{W}}_{i}\right)$
$\displaystyle=-c_{\xi
1}s_{\xi_{i}}^{T}\text{sg}\left(s_{\xi_{i}}\right)-c_{\xi
2}\left\|s_{\xi_{i}}\right\|^{2}+s_{\xi_{i}}^{T}\sigma_{i}+s_{\xi_{i}}^{T}\tilde{W}_{i}^{T}H_{i}-\text{tr}\left(\tilde{W}_{i}^{T}a^{-1}\dot{\hat{W}}_{i}\right)$
$\displaystyle=-c_{\xi
1}s_{\xi_{i}}^{T}\text{sg}\left(s_{\xi_{i}}\right)-c_{\xi
2}\left\|s_{\xi_{i}}\right\|^{2}+s_{\xi_{i}}^{T}\sigma_{i}+\text{tr}\left(-\tilde{W}_{i}^{T}\left(a^{-1}\dot{\hat{W}}_{i}-H_{i}s_{\xi_{i}}^{T}\right)\right)$

With the updated rule of the neural network, $\dot{V}_{\xi_{i}}$ can be
rewritten as follows:

$\displaystyle\dot{V}_{\xi_{i}}$ $\displaystyle=-c_{\xi
1}s_{\xi_{i}}^{T}\text{sg}\left(s_{\xi_{i}}\right)-c_{\xi
2}\left\|s_{\xi_{i}}\right\|^{2}+s_{\xi_{i}}^{T}\sigma_{i}+\text{tr}\left(\tilde{W}_{i}^{T}\eta\left\|s_{\xi_{i}}\right\|\hat{w}_{i}\right)$
(27) $\displaystyle=-c_{\xi
1}s_{\xi_{i}}^{T}\text{sg}\left(s_{\xi_{i}}\right)-c_{\xi
2}\left\|s_{\xi_{i}}\right\|^{2}+s_{\xi_{i}}^{T}\sigma_{i}+\eta\left\|s_{\xi_{i}}\right\|\text{tr}\left(\tilde{W}_{i}^{T}\left(W_{i}-\tilde{W}_{i}\right)\right)$

According to the Cauchy-Schwarz inequality, the following inequality equation
can be satisfied:

$\text{tr}\left(\tilde{W}_{i}^{T}\left(W_{i}-\tilde{W}_{i}\right)\right)\leq\left\|\tilde{W}_{i}\right\|\left\|W_{i}\right\|-\left\|\tilde{W}_{i}\right\|^{2}$
(28)

Thus,

$\displaystyle\dot{V}_{\xi_{i}}$ $\displaystyle\leq-c_{\xi
1}s_{\xi_{i}}^{T}\text{sg}\left(s_{\xi_{i}}\right)-c_{\xi
2}\left\|s_{\xi_{i}}\right\|^{2}+\left\|s_{\xi_{i}}\right\|\bar{\sigma}+\eta\left\|s_{\xi_{i}}\right\|\left(\left\|\tilde{W}_{i}\right\|\left\|W_{i}\right\|-\left\|\tilde{W}_{i}\right\|^{2}\right)$
(29) $\displaystyle\leq-c_{\xi
1}s_{\xi_{i}}^{T}\text{sg}\left(s_{\xi_{i}}\right)-c_{\xi
2}\left\|s_{\xi_{i}}\right\|^{2}+\left\|s_{\xi_{i}}\right\|\bar{\sigma}+\dfrac{1}{4}\eta\left\|s_{\xi_{i}}\right\|\bar{W}_{i}^{2}-\eta\left\|s_{\xi_{i}}\right\|\left(\dfrac{1}{2}\left\|W_{i}\right\|-\left\|\tilde{W}_{i}\right\|\right)^{2}$

Based on an extension of the Lyapunov theorem [33],
$\left\|s_{\xi_{i}}\right\|$ is bounded. Moreover, the control gain $c_{\xi
2}$ can be selected large enough so that

$\left[{\bar{\sigma}+\eta\bar{W}_{i}^{2}/{4}}\right]/{c_{\xi 2}}\leq b_{\xi}$
(30)

Therefore, with the inequality condition (24), $\dot{V}_{p}$ can be rewritten
as follows:

$\dot{V}_{\xi_{i}}\leq-c_{\xi
1}s_{\xi_{i}}^{T}\text{sg}\left(s_{\xi_{i}}\right)-\eta\left\|s_{\xi_{i}}\right\|\left(\dfrac{1}{2}\left\|W_{i}\right\|-\left\|\tilde{W}_{i}\right\|\right)^{2}\leq
0$ (31)

The Lyapunov stability condition is satisfied.

###  3.2 Attitude controller design

In our system, the position controller is the outer loop of the UAV control
system, as depicted in Figure 3. Its control signal is then fed to the
converter block to calculate the desired angles and translational forces based
on (5) as:

$\displaystyle{}^{d}\theta_{i}$
$\displaystyle=\arctan\left(\dfrac{u_{xi}\cos{}^{d}\psi_{i}+u_{yi}\sin{}^{d}\psi_{i}}{u_{zi}+g}\right)$
(32) $\displaystyle{}^{d}\phi_{i}$
$\displaystyle=\arctan\left(\cos{}^{d}\theta_{i}\dfrac{u_{xi}\sin{}^{d}\psi_{i}-u_{yi}\cos{}^{d}\psi_{i}}{u_{zi}+g}\right)$
$\displaystyle f_{ti}$
$\displaystyle=\dfrac{u_{zi}+g}{\cos{}^{d}\phi_{i}\cos{}^{d}\theta_{i}}$

They are used as the reference for the attitude controller, which is designed
based on the BSMC. From (2), the dynamic equation for the roll angle is given
by:

$\displaystyle\ddot{\phi}_{i}$
$\displaystyle=\dfrac{\dot{\theta}_{i}\dot{\psi}_{i}\left(I_{y}-I_{z}\right)+\tau_{\phi_{i}}}{I_{x}}.$
(33)

Denote ${}^{e}\phi_{i}=\phi_{i}-{}^{d}\phi_{i}$ as the roll angle error. The
virtual velocity, ${}^{v}\phi_{i}$, is defined as:

${}^{v}\phi_{i}={}^{d}\dot{\phi}_{i}-\lambda_{\phi}\phi_{1e},$ (34)

where $\lambda_{\phi}>0$ is a positive gain. The first candidate Lyapunov
function for subsystem $\phi_{1e}$ is chosen as:

${}^{1}V_{\phi_{i}}=\dfrac{1}{2}{}^{e}\phi_{i}^{2}.$ (35)

Taking the first derivative of ${}^{1}V_{\phi}$ gives:

${}^{1}\dot{V}_{\phi_{i}}={}^{e}\phi_{i}{}^{e}\dot{\phi}_{i}={}^{e}\phi_{i}\left(\dot{\phi}_{i}-{}^{d}\dot{\phi}_{i}\right).$
(36)

Substituting $\dot{\phi}_{i}={}^{v}\phi_{i}$ into (36) gives

${}^{1}\dot{V}_{\phi_{i}}=-\lambda_{\phi}{}^{e}\phi_{i}^{2}\leq 0.$ (37)

Thus, the Lyapunov stability is guaranteed. The sliding surface of the roll
angle subsystem is expressed as:

$s_{\phi_{i}}=\gamma_{\phi}{}^{e}\phi_{i}+\left(\dot{\phi}_{i}-{}^{e}\phi_{i}\right),$
(38)

where $\gamma_{\phi}>0$ is a positive gain. The first derivative of
$s_{\phi_{i}}$ is given by:

$\dot{s}_{\phi_{i}}=\gamma_{\phi}\left(\dot{\phi}_{i}-{}^{d}\dot{\phi}_{i}\right)+\dfrac{\dot{\theta}_{i}\dot{\psi}_{i}\left(I_{y}-I_{z}\right)+\tau_{\phi_{i}}}{I_{x}}-{}^{v}\dot{\phi}_{i}.$
(39)

The control signal is then designed with two sub-control signals,
${}^{eq}\tau_{\phi_{i}}$ and ${}^{sw}\tau_{\phi_{i}}$.
${}^{eq}\tau_{\phi_{i}}$ is the equivalent control signal that maintains the
roll angle on the sliding manifold and ${}^{sw}\tau_{\phi_{i}}$ is the signal
that leads the subsystem to the sliding surface $s_{\phi_{i}}$. They are
chosen as follows:

$\displaystyle{}^{eq}\tau_{\phi_{i}}$
$\displaystyle=I_{x}\left({}^{v}\dot{\phi}_{i}-\gamma_{\phi}\left(\dot{\phi}_{i}-{}^{d}\dot{\phi}_{i}\right)\right)-\dot{\theta}_{i}\dot{\psi}_{i}\left(I_{y}-I_{z}\right)$
(40) $\displaystyle{}^{sw}\tau_{\phi_{i}}$ $\displaystyle=-I_{x}\left(c_{\phi
1}\text{sg}\left(s_{\phi_{i}}\right)+c_{\phi 2}s_{\phi_{i}}\right),$

where $c_{\phi 1}$ and $c_{\phi 2}$ are positive gains.

######  Theorem 3.5.

Consider the roll angle subsystem (33). If the control signal is designed as:

$\tau_{\phi_{i}}={}^{eq}\tau_{\phi_{i}}+{}^{sw}\tau_{\phi_{i}},$ (41)

the roll angle control system is stable.

######  Proof 3.6.

The candidate Lyapunov function of the roll angle subsystem is chosen as
follows:

$V_{\phi_{i}}=\dfrac{1}{2}s_{\phi_{i}}^{2}.$ (42)

Taking the first derivative of $V_{\phi}$ gives

$\dot{V}_{\phi_{i}}=s_{\phi_{i}}\dot{s}_{\phi_{i}}.$ (43)

By substituting (39) into (43), we have

$\dot{V}_{\phi_{i}}=s_{\phi_{i}}\left(\gamma_{\phi}\left(\dot{\phi}_{i}-{}^{d}\dot{\phi}_{i}\right)+\dfrac{\dot{\theta}_{i}\dot{\psi}_{i}\left(I_{y}-I_{z}\right)+\tau_{\phi_{i}}}{I_{x}}-{}^{v}\dot{\phi}_{i}\right).$
(44)

Finally, substituting (40) and (41) into (44) gives

$\dot{V}_{\phi_{i}}=-c_{\phi
1}s_{\phi_{i}}\text{sg}\left(s_{\phi_{i}}\right)-c_{\phi
2}s_{\phi_{i}}^{2}\leq 0.$ (45)

Thus, the Lyapunov stability of the roll angle control system is guaranteed.

The control signals for the pitch and yaw angles can be obtained by applying
the design process similar to the roll angle. As a result, the pitch control
signals are obtained as:

$\displaystyle\tau_{\theta_{i}}$
$\displaystyle={}^{eq}\tau_{\theta_{i}}+{}^{sw}\tau_{\theta_{i}}$ (46)
$\displaystyle{}^{eq}\tau_{\theta_{i}}$
$\displaystyle=I_{y}\left({}^{v}\dot{\theta}_{i}-\gamma_{\theta}\left(\dot{\theta}_{i}-{}^{d}\dot{\theta}_{i}\right)\right)-\dot{\phi}_{i}\dot{\psi}_{i}\left(I_{z}-I_{x}\right)$
$\displaystyle{}^{sw}\tau_{\theta_{i}}$ $\displaystyle=-I_{y}\left(c_{\theta
1}\text{sg}\left(s_{\theta_{i}}\right)+c_{\theta 2}s_{\theta_{i}}\right),$

and the yaw control signals are given by:

$\displaystyle\tau_{\psi_{i}}$
$\displaystyle={}^{eq}\tau_{\psi_{i}}+{}^{sw}\tau_{\psi_{i}}$ (47)
$\displaystyle{}^{eq}\tau_{\psi_{i}}$
$\displaystyle=I_{z}\left({}^{v}\dot{\psi}_{i}-\gamma_{\psi}\left(\dot{\psi}_{i}-{}^{d}\dot{\psi}_{i}\right)\right)-\dot{\phi}_{i}\dot{\theta}_{i}\left(I_{x}-I_{y}\right)$
$\displaystyle{}^{sw}\tau_{\psi_{i}}$ $\displaystyle=-I_{z}\left(c_{\psi
1}\text{sg}\left(s_{\psi_{i}}\right)+c_{\psi 2}s_{\psi_{i}}\right).$

##  4 Results

Table 1: Parameters of the Hummingbird quadrotor Description |  Notation |  Value
---|---|---
The mass of UAV (kg) |  $m$ |  0.68
The gravitational acceleration (m/s2) |  $g$ |  9.81
Moment of inertia about x and y axes (kg.m2) |  $I_{x}$, $I_{y}$ |  0.007
Moment of inertia about z-axis (kg.m2) |  $I_{z}$ |  0.012
Length of UAV arm (m) |  $l$ |  0.17
Thrust coefficient |  $k_{t}$ |  29$\times$10-6
Drag coefficient |  $k_{d}$ |  1.1$\times$10-6
Table 2: Parameters of the proposed controller Parameters of the BSMC
---
Attitude controller |  |  $\lambda_{\phi}=\lambda_{\theta}=\lambda_{\psi}=5$;
---
$\gamma_{\phi}=\gamma_{\theta}=\gamma_{\psi}=5$;
$c_{\phi 1}=c_{\theta 1}=c_{\psi 1}=2$;
$c_{\phi 2}=c_{\theta 2}=c_{\psi 2}=2$;
Position controller |  $\lambda_{\xi}=3$; $\gamma_{\xi}=2$; $c_{\xi 1}=c_{\xi 2}=2$
Parameters of the RBFNN
Position controller |  | $m=50;a=0.1;b=10;\eta=0.2;$
---
$\left\\{\begin{array}[]{c}\mu_{x}=\text{linspace}\left(-r_{x},r_{x},n\right)\\\
\mu_{y}=\text{linspace}\left(-r_{y},r_{y},n\right)\\\
\mu_{z}=\text{linspace}\left(-r_{z},r_{z},n\right)\end{array}\right.;r_{x}=r_{y}=r_{z}=15;$
$\mu_{1}=\left[\mu_{x},\mu_{y},\mu_{z}\right]^{T};\mu_{2}=\mu_{1}$;
Figure 5: The desired topology Figure 6: The external disturbance acting on
the formation generated based on the rectangle and full wavelength “1-cosine”
wind model in Scenario 1

(a) RBF-BSMC

(b) MPC

(c) BSMC

(d) SMC

Figure 7: Trajectories of the UAV formation generated by the four controllers
in Scenario 1

(a) The tracking errors

(b) The estimated disturbances

Figure 8: The tracking errors and estimated disturbances of the controllers in
Scenario 1 Figure 9: The external disturbance acting on the formation
generated based on the rectangle wind model in Scenario 2

(a) RBF-BSMC

(b) MPC

(c) BSMC

(d) SMC

Figure 10: Trajectories of the UAV formation generated by the four controllers
in Scenario 2

(a) The tracking errors

(b) The estimated disturbances

Figure 11: The tracking errors and estimated disturbances of the controllers
in Scenario 2

(a) The real Hummingbird drone

(b) The Hummingbird model

Figure 12: The Hummingbird drone model used in validation [32, 34]

(a) Vertical formation

(b) Triangular formation

Figure 13: The vertical and triangular formation topologies used in SIL tests

(a) 3D view

(b) Side view

(c) Top view

Figure 14: The planned paths to inspect the bridge Figure 15: The tracking
errors of UAV formation in the Gazebo SIL

To evaluate the performance of the proposed control system, we have conducted
a number of evaluations and comparisons111Evaluation results in Scenario 1 and
Scenario 2 - https://youtu.be/LYD7269n1-c. The UAV model used is the
Hummingbird quadrotors [34], whose parameters are shown in Table 1. Parameters
of the position and attitude controllers are chosen as shown in Table 2. The
desired formation is a triangular shape with
$\Delta_{1}=\left[2,0,0\right]^{T}$, $\Delta_{2}=\left[0,0,2\right]^{T}$, and
$\Delta_{3}=\left[0,0,-2\right]^{T}$, as depicted in Figure 5. Comparisons are
conducted between the proposed controller (RBF-BSMC) and three other
controllers namely model predictive control (MPC) [35, 27], backstepping
sliding mode control (BSMC) [36, 37] and sliding mode control (SMC) [23, 38]
in different scenarios.

###  4.1 Scenario 1

In this scenario, external disturbances acting on the formation are generated
based on the combination of the rectangle and full wavelength “1-cosine” wind
model [39], as shown in Figure 6. The initial positions of the UAVs are set as
$\xi_{1}=\left[3,2,4\right]^{T}$, $\xi_{2}=\left[2,1,4\right]^{T}$, and
$\xi_{3}=\left[0,0,4\right]^{T}$. The desired trajectory of the virtual leader
is a spiral with the $z$ coordinate increasing over time as expressed in (48).

$\displaystyle x_{L}$ $\displaystyle=5\cos\left(\dfrac{2\pi}{20}t\right)\text{
(m)},$ (48) $\displaystyle y_{L}$
$\displaystyle=5\sin\left(\dfrac{2\pi}{20}t\right)\text{ (m)},$ $\displaystyle
z_{L}$ $\displaystyle=0.5t+5\text{ (m)},$ $\displaystyle\psi_{L}$
$\displaystyle=\dfrac{2\pi}{20}t+\dfrac{\pi}{2}\text{ (rad)}.$

Figure 7 shows the 3D views of the trajectory tracking results of the UAV
formation. It can be seen that the UAVs quickly reach the initial positions to
form the desired shape. They then maintain the shape while following the
reference trajectory. However, the trajectory of the proposed method is
smoother and more accurate than the others due to its capability to estimate
the disturbance via the RBFNN and use it as feedback to adjust the control
signals. This result can be further verified via the tracking errors as shown
in Figure 8(a). It can be seen that the average tracking errors of the
proposed controller quickly converge to zeros, whereas those errors of the
other controllers largely fluctuate due to disturbances. In addition, the
maximum and minimum tracking errors of the UAVs are also very small with our
method, which confirm its stability for formation control.

Figure 8(b) shows the disturbances estimated for each UAV by the proposed
controller. After the transition period, the estimation starts to converge to
the real disturbance. This provides feedback for the controller to adjust the
control signal for better tracking performance.

###  4.2 Scenario 2

In this scenario, the rectangle wind model is used to generate external
disturbances as shown in Figure 9. The initial positions of the UAVs are set
as $\xi_{1}=\left[2,0,9\right]^{T}$, $\xi_{2}=\left[5,2,7\right]^{T}$, and
$\xi_{3}=\left[-3,-2,8\right]^{T}$. The desired trajectory of the virtual
leader is generated based on an inspection path used to collect surface data
[40].

The 3D of the UAV formation tracking results are shown in Figure 10. It can be
seen that all controllers are able to drive the UAVs to reach their reference
positions and then track them to form the desired shape. The proposed
controller, however, introduces smaller tracking errors than the other
controllers, as shown in Figure 11(a), due to its disturbance estimator. As
shown in Figure 11(b), the estimation closely follows the actual disturbances
except for the positions where step changes happen. At those positions, a
transition period is needed for the estimator to converge to the new steady
state. However, the settling time of the estimator is short allowing it to
provide timely feedback to the controller.

###  4.3 Validation with software-in-the-loop tests

To further validate the proposed control system, we have carried out software-
in-the-loop (SIL) tests that involve the inspection of a scaled-down 3D model
of a real bridge with 5 columns, as shown in Figure 14. The UAV model used is
a Hummingbird quadrotor222Source code used for Gazebo validation -
https://github.com/duynamrcv/hummingbird_simulator developed based on Gazebo-
based RotorS simulator [32], as depicted in Figure 12. The formation used
includes two topologies, vertical and triangular shapes, as shown in Figure
13. According to our previous work [41], the generated path to inspect the
bridge includes two stages. The first stage covers all columns of the bridge
using a vertical formation. The second one uses a triangular formation to
cover the side and top surfaces of the bridge, as depicted in Figure 14.

Given the planned paths, the UAVs start to fly from positions
$\left[-8,-8,0\right]^{T}$, $\left[-4,-8,0\right]^{T}$,
$\left[4,-8,0\right]^{T}$, and reach their initial positions to form a
vertical formation as shown in Figure 14(a). The formation then tracks the
planned path to acquire surface images of the bridge333SIL validation -
https://youtu.be/1yUCzWRDcp0. Figure 15 shows the tracking errors of the UAVs
during operation. It can be seen that the errors quickly converge to small
values in both inspection stages except between time steps 1070 and 1120,
where there is a change in the formation topology. The errors in the first
stage are slightly larger than in the second one as the UAVs frequently
changes their direction to navigate around each column of the bridge.
Nevertheless, the average tracking error of less than 5 cm is sufficient for
most UAV-related applications and thus confirms the validity of our approach.

##  5 Conclusion

In this paper, we have presented a robust control system using RBFNN for a
group of UAVs flying in a formation. By combining BSC with SMC, the controller
can handle nonlinearity to increase its control performance. The use of RBFNN
enables the system to estimate external disturbances to enhance its control
robustness. By using Lyapunov’s theorem, we proved that the control system is
stable and the proposed controller can track the reference trajectory.
Evaluation results show that the proposed controller outperforms the state-of-
the-art BSMC in terms of accuracy and robustness and is sufficient for most
UAV applications.

Author Contributions Duy-Nam Bui: Conceptualization, Methodology,
Implementation, Writing - original draft. Manh Duong Phung: Investigation,
Conceptualization, Supervision, Writing - review and editing.

Financial Support Duy-Nam Bui was funded by the Master, PhD Scholarship
Programme of Vingroup Innovation Foundation (VINIF), code VINIF.2022.ThS.057.

Conflicts of Interest The authors declare no conflicts of interest exist.

## References

  * [1] J. Zeng, H. Zhong, Y. Wang, S. Fan, and H. Zhang, “Autonomous control design of an unmanned aerial manipulator for contact inspection,” Robotica, vol. 41, no. 4, p. 1145–1158, 2023.
  * [2] M. Rizia, J. A. Reyes-Munoz, A. G. Ortega, A. Choudhuri, and A. Flores-Abad, “Autonomous aerial flight path inspection using advanced manufacturing techniques,” Robotica, vol. 40, no. 7, p. 2128–2151, 2022.
  * [3] L. Inzerillo, G. Di Mino, and R. Roberts, “Image-based 3D reconstruction using traditional and UAV datasets for analysis of road pavement distress,” Automation in Construction, vol. 96, pp. 457–469, 2018.
  * [4] S. Zhao, F. Kang, J. Li, and C. Ma, “Structural health monitoring and inspection of dams based on UAV photogrammetry with image 3D reconstruction,” Automation in Construction, vol. 130, p. 103832, 2021.
  * [5] K. Chen, G. Reichard, X. Xu, and A. Akanmu, “Automated crack segmentation in close-range building façade inspection images using deep learning techniques,” Journal of Building Engineering, vol. 43, p. 102913, 2021.
  * [6] X. Peng, X. Zhong, C. Zhao, A. Chen, and T. Zhang, “A UAV-based machine vision method for bridge crack recognition and width quantification through hybrid feature learning,” Construction and Building Materials, vol. 299, p. 123896, 2021.
  * [7] H. M. La, T. H. Dinh, N. H. Pham, Q. P. Ha, and A. Q. Pham, “Automated robotic monitoring and inspection of steel structures and bridges,” Robotica, vol. 37, no. 5, p. 947–967, 2019.
  * [8] Y. Tian, G. Zhang, K. Morimoto, and S. Ma, “Automated rust removal: Rust detection and visual servo control,” Automation in Construction, vol. 134, p. 104043, 2022.
  * [9] W. Jing, D. Deng, Y. Wu, and K. Shimada, “Multi-UAV coverage path planning for the inspection of large and complex structures,” in 2020 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pp. 1480–1486, 2020.
  * [10] G. Silano, T. Baca, R. Penicka, D. Liuzza, and M. Saska, “Power line inspection tasks with multi-aerial robot systems via signal temporal logic specifications,” IEEE Robotics and Automation Letters, vol. 6, no. 2, pp. 4169–4176, 2021.
  * [11] V. Hoang, M. Phung, T. Dinh, and Q. Ha, “Angle-encoded swarm optimization for UAV formation path planning,” in 2018 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pp. 5239–5244, 2018.
  * [12] K.-K. Oh, M.-C. Park, and H.-S. Ahn, “A survey of multi-agent formation control,” Automatica, vol. 53, pp. 424–440, 2015.
  * [13] Y. Liu and R. Bucknall, “A survey of formation control and motion planning of multiple unmanned vehicles,” Robotica, vol. 36, pp. 1019–1047, Mar. 2018.
  * [14] V. T. Hoang, M. D. Phung, T. H. Dinh, and Q. P. Ha, “System architecture for real-time surface inspection using multiple UAVs,” IEEE Systems Journal, vol. 14, no. 2, pp. 2925–2936, 2020.
  * [15] J. Zheng, X. Zong, H. Ge, Z. Zheng, and M. C. Makuwatsine, “Virtual leader-follower synchronization controller design for distributed parameter multi-agent systems with time-varying disturbances,” Neurocomputing, vol. 450, pp. 389–398, 2021.
  * [16] F. Rinaldi, S. Chiesa, and F. Quagliotti, “Linear quadratic control for quadrotors UAVs dynamics and formation flight,” Journal of Intelligent & Robotic Systems, vol. 70, pp. 203–220, Apr 2013.
  * [17] J. Chen, D. Sun, J. Yang, and H. Chen, “Leader-follower formation control of multiple non-holonomic mobile robots incorporating a receding-horizon scheme,” The International Journal of Robotics Research, vol. 29, no. 6, pp. 727–747, 2010.
  * [18] B.-S. Chen, Y.-C. Liu, M.-Y. Lee, and C.-L. Hwang, “Decentralized h PID team formation tracking control of large-scale quadrotor UAVs under external disturbance and vortex coupling,” IEEE Access, vol. 10, pp. 108169–108184, 2022.
  * [19] Y. Chen and T. Deng, “Leader-follower UAV formation flight control based on feature modelling,” Systems Science & Control Engineering, vol. 11, Oct. 2023.
  * [20] F. Fahimi, “Full formation control for autonomous helicopter groups,” Robotica, vol. 26, pp. 143–156, Mar. 2008.
  * [21] A. K. Khalaji and R. Zahedifar, “Lyapunov-based formation control of underwater robots,” Robotica, vol. 38, pp. 1105–1122, Aug. 2019.
  * [22] M. A. Dehghani and M. B. Menhaj, “Integral sliding mode formation control of fixed-wing unmanned aircraft using seeker as a relative measurement system,” Aerospace Science and Technology, vol. 58, pp. 318–327, 2016.
  * [23] M. Defoort, T. Floquet, A. Kokosy, and W. Perruquetti, “Sliding-mode formation control for cooperative autonomous mobile robots,” IEEE Transactions on Industrial Electronics, vol. 55, no. 11, pp. 3944–3953, 2008.
  * [24] R. Li, L. Zhang, L. Han, and J. Wang, “Multiple vehicle formation control based on robust adaptive control algorithm,” IEEE Intelligent Transportation Systems Magazine, vol. 9, no. 2, pp. 41–51, 2017.
  * [25] Y. Huang, W. Liu, B. Li, Y. Yang, and B. Xiao, “Finite-time formation tracking control with collision avoidance for quadrotor UAVs,” Journal of the Franklin Institute, vol. 357, no. 7, pp. 4034–4058, 2020.
  * [26] X. Wang, S. Baldi, X. Feng, C. Wu, H. Xie, and B. De Schutter, “A fixed-wing UAV formation algorithm based on vector field guidance,” IEEE Transactions on Automation Science and Engineering, vol. 20, no. 1, pp. 179–192, 2023.
  * [27] Q. Yuan and X. Li, “Distributed model predictive formation control for a group of UAVs with spatial kinematics and unidirectional data transmissions,” IEEE Transactions on Network Science and Engineering, vol. 10, no. 6, pp. 3209–3222, 2023.
  * [28] S. Yang, W. Bai, T. Li, Q. Shi, Y. Yang, Y. Wu, and C. L. P. Chen, “Neural-network-based formation control with collision, obstacle avoidance and connectivity maintenance for a class of second-order nonlinear multi-agent systems,” Neurocomputing, vol. 439, pp. 243–255, 2021.
  * [29] K. Shojaei, “Neural network formation control of Underactuated Autonomous Underwater vehicles with saturating actuators,” Neurocomputing, vol. 194, pp. 372–384, 2016.
  * [30] C.-W. Kuo, C.-C. Tsai, and C.-T. Lee, “Intelligent leader-following consensus formation control using recurrent neural networks for small-size unmanned helicopters,” IEEE Transactions on Systems, Man, and Cybernetics: Systems, vol. 51, no. 2, pp. 1288–1301, 2021.
  * [31] E. J. Hartman, J. D. Keeler, and J. M. Kowalski, “Layered neural networks with gaussian hidden units as universal approximations,” Neural Computation, vol. 2, no. 2, pp. 210–215, 1990.
  * [32] F. Furrer, M. Burri, M. Achtelik, and R. Siegwart, Robot Operating System (ROS): The Complete Reference (Volume 1), ch. RotorS—A Modular Gazebo MAV Simulator Framework, pp. 595–625.  Cham: Springer International Publishing, 2016.
  * [33] K. S. NARENDRA and A. M. ANNASWAMY, “Persistent excitation in adaptive systems,” International Journal of Control, vol. 45, pp. 127–160, Jan. 1987.
  * [34] D. N. Bui, T. T. Van Nguyen, and M. D. Phung, “Lyapunov-based nonlinear model predictive control for attitude trajectory tracking of unmanned aerial vehicles,” International Journal of Aeronautical and Space Sciences, Oct 2022.
  * [35] D. Wang, Q. Pan, Y. Shi, J. Hu, and C. Zhao, “Efficient nonlinear model predictive control for quadrotor trajectory tracking: Algorithms and experiment,” IEEE Transactions on Cybernetics, vol. 51, no. 10, pp. 5057–5068, 2021.
  * [36] L.-X. Xu, H.-J. Ma, D. Guo, A.-H. Xie, and D.-L. Song, “Backstepping sliding-mode and cascade active disturbance rejection control for a quadrotor UAV,” IEEE/ASME Transactions on Mechatronics, vol. 25, no. 6, pp. 2743–2753, 2020.
  * [37] D. J. Almakhles, “Robust backstepping sliding mode control for a quadrotor trajectory tracking application,” IEEE Access, vol. 8, pp. 5515–5525, 2020.
  * [38] I. Ahmad, M. Liaquat, F. M. Malik, H. Ullah, and U. Ali, “Variants of the sliding mode control in presence of external disturbance for quadrotor,” IEEE Access, vol. 8, pp. 227810–227824, 2020.
  * [39] B. H. Wang, D. B. Wang, Z. A. Ali, B. T. Ting, and H. Wang, “An overview of various kinds of wind effects on unmanned aerial vehicle,” Measurement and Control, vol. 52, no. 7-8, pp. 731–739, 2019.
  * [40] M. D. Phung, C. H. Quach, T. H. Dinh, and Q. Ha, “Enhanced discrete particle swarm optimization path planning for UAV vision-based surface inspection,” Automation in Construction, vol. 81, pp. 25–33, 2017.
  * [41] D. N. Bui, T. N. Duong, and M. D. Phung, “Ant colony optimization for cooperative inspection path planning using multiple unmanned aerial vehicles,” in 2024 IEEE/SICE International Symposium on System Integration (SII), pp. 675–680, 2024.