# Impact of large scale climate oscillation on drought in West Africa

Samuel, T. Ogunjo<EMAIL_ADDRESS>+2348026009902 Christiana, F.
Olusegun  Ibiyinka, A. Fuwape  Department of Physics, Federal University of
Technology, Akure, Ondo State, Nigeria  Michael and Cecilia Ibru University,
Ughelli North, Delta State, Nigeria.  Department of Meteorology and Climate
Science, West African Science Service Centre on Climate Change and Adapted
Land Use (WASCAL), FUTA

###### Abstract

Drought poses a significant threat to the delicate economies in subsaharan
Africa. This study investigates the influence of large scale ocean oscillation
on drought in West Africa. Standardized Precipitation Index for the region was
computed using monthly precipitation data from the Climate Research Unit
during the period 1961 -1990. The impact of three ocean oscillation indices -
Southern Ocean Index (SOI), Pacific Decadal Oscillation (PDO) and North
Atlantic Oscillation (NAO) on drought over West Africa was investigated using
linear correlation, co-integration test, mutual information and nonlinear
synchronization methods. SOI showed predominantly positive correlation with
drought over the region while PDO and NAO showed negative correlation. This
was confirmed by the co-integration tests. The nonlinear test revealed more
complex relationship between the indices and drought. PDO has lesser influence
or contribute less to the drought in the coastal region compared to the Sahel
region of West Africa.

††journal: Journal of Hydrology

##  1 Introduction

Variations in local and global climate are influenced by large scale
(atmospheric and ocean) oscillations [1]. Pacific Decadal Oscillation (PDO),
North Atlantic Oscillations (NAO) and El Nino Southern Oscillations (ENSO) are
some of the large scale oscillations whose impacts are felt across different
parts of Africa. The North Atlantic Oscillation is a measure of the variations
in lower atmospheric pressure over Iceland and higher pressure near the Azores
and Iberian Peninsula [2]. NAO has been found to have impact on the
hydroclimatic conditions in the Mediterranean region [3], export of dust from
North Africa [4], precipitation in southeastern Africa [5], European airline
traffic [6], harmful algae blooms [7] and population density of migratory
birds [8]. Over West Africa, strong correlation has been found between
Atlantic Nino Index 1 and drought indices based on CORDEX models [9].
According to Anyamba et al. [10], South Atlantic SST and Global tropics SST
show significant strong correlation with drought in Normalized Difference
Vegetation Index (NDVI) over East Africa. Positive correlation between ENSO
and NDVI have also been reported for the have also been reported for the
African continent [11]. Studies have established coherence between NAO and
ENSO [12], tropical SSTs [13] and East Asian summer monsoon-ENSO relationship
[14]. NAO has stronger influence over agricultural yields of selected crops in
Spain than ENSO [15] as well as drought in different parts of Europe [16].

The Southern Oscillation Index is a large scale oscillation based on the sea
level pressure between Tahiti and Darwin Islands. The SOI is associated with
sustained periods of negative (positive) episodes referred to as El Nino (La
Nina) and occurs every 2-7 years. SOI, though originating in pacific, has been
recognized as a major climatic factor that influences global weather events.
SOI has been shown to have influence on precipitation in Australia with a lag
of 6 months [17], inverse multiscale relationship with winter drought in China
[18], precipitation in South Korea [19], rice production in India [20], wave
climate and beach rotation [21] and extreme temperature [22]. PDO, which
usually occurs every 20 - 30 years, is observed $20^{o}N$ of the North Pacific
ocean. PDO has dominant influence in the North Pacific with secondary
signature in the tropics while the effect of SOI is predominantly in the
tropics [23]. It is related to, affected by but distinct from SOI [24].
Relationship has been found between PDO and fish population in the Northwest
Pacific [25], global precipitation in the extratropics [26] and climate
variability in Mexico [27]. It has been suggested that there exist a
relationship between SOI and PDO [28]. Indirect influence of PDO on
precipitation in continental United States through its influence on SOI has
been proposed [26]. Studies on drought within West Africa include
investigation into the impact of ENSO on precipitation in Nigeria [29], impact
of drought in the Volta Basin [30] and impact of Hardley Circulation on
drought in Africa [31]. Low-frequency signals in global Sea Surface
Temperatures are now recognized as major drivers of the Sahelian drought [32].

Large scale oscillations and climate variables have been shown to be either
stochastic or chaotic [33, 34, 35, 36]. The use of linear analysis in their
investigation will lead to loss of information. For example, Fuwape et al.
[37] showed that the use of correlation overestimate the relationship between
precipitation and radio refractivity in the tropics. It is important to
investigate complex and chaotic tendencies of climate time series for
efficient short and long term predictions, model testing and better
description [35]. Chaotic systems are sensitive to initial conditions and
characterized by positive Lyapunov exponents. Chaoticity and detrending of
these systems tend to introduce spurious correlations and cross-correlations
which will have significant impact on the results obtained in such studies and
their interpretation [38, 39]. Recent development in chaos theory has provided
a veritable tool to effectively study these nonlinear events without loss of
information.

In this study, we aim to investigate the linear and nonlinear relationship
between three large scale oscillation (PDO, NAO and SOI) on drought in the
West African region. A study of linear analysis (correlation and cointegration
tests) and dynamical tests (mutual information and phase synchronization) will
not only give the performance of both tools but a robust relationship between
large scale oscillations and droughts.

##  2 Dataset and Methodology

###  2.1 Study area

The study area for this research is West Africa. It lies within latitude
$0^{o}N$ and $23.5^{o}N$ and longitudes $20^{o}W$ and $15^{o}E$ (Figure 1). It
is bounded in the south and west by the Atlantic Ocean and landed areas in the
North and East. The hydrographic network consist of short seasonal rivers
around the Guinean Highlands, Gulf of Guinea, Niger Delta and Volta rivers.
Omotosho [40] reported that deep convective system is responsible for 70% of
rainfall in the West African region. The tropical diurnal convective systems
that drives the West African Moonson involves the movement of moisture and
heat between the Atlantic Ocean and continental land [41]. Distribution of
rainfall within the region is influenced by the position of the Intertropical
discontinuity (ITD). The ITD is separation between the tropical maritime air
mass and the dry tropical continental air mass. It is a low pressure system
which controls the spatial variation in tropical rainfall [42].

Figure 1: Map of study area showing the mean monthly precipitation (mm) during
the period 1961- 1990. The inset shows the study location within the continent
of Africa. Rivers are shown in thick black lines.

###  2.2 Data

Monthly time series for North Atlantic Oscillation (NAO) Index, Southern
Oscillation (SOI) Index, and Pacific Decadal Oscillations (PDO) were obtained
for 1961 - 1990 from The Earth System Research Laboratory, NOAA
https://www.esrl.noaa.gov/psd/data/climateindices/list/. Monthly precipitation
data used was the version 4.01 of the Climate Research Unit Time Series (CRU
TS - http://data.ceda.ac.uk//badc/cru/data/cru_ts/cru_ts_4.01/data/pre/) at a
resolution of $0.5^{o}$ resolution from 1961 - 1990 [43]. There were no
missing data in both the climatic indices and precipitation data set.
Anomalies of the climatic indices were used to remove annual variations. The
temporal variation of the three climate indices under the period of study is
shown in Figure 2.

Figure 2: Monthly time series of climate indices SOI (upper panel), PDO
(middle panel), NAO (lower panel).

###  2.3 Standard Precipitation Index

SPI was introduced by McKee et al. [44] to provide a drought indicator at
different timescales based on probability. In this approach, the precipitation
data for a location is fitted to a gamma or Pearson Type III distribution
before transforming to a normal distribution. For precipitation values equal
zero, the cumulative probability of the form

$H(x)=1+(1-q)G(x)$ (1)

where $q$ is the probability of zero precipitation [45, 46]. The values of SPI
is interpreted as

$SPI=\left\\{\begin{array}[]{ll}SPI\geq 2.00,&\hbox{Extremely wet;}\\\
1.50\leq SPI<2.00,&\hbox{Very wet;}\\\ 1.00\leq SPI<1.50,&\hbox{Moderately
wet;}\\\ -1.00\leq SPI<1.00,&\hbox{Normal;}\\\ -1.50\leq
SPI<-1.00,&\hbox{Moderately drought;}\\\ -2.00\leq SPI<-1.50,&\hbox{Severely
drought;}\\\ SPI<-2.00,&\hbox{Extreme drought.}\end{array}\right.$

The use of SPI was recommend for drought monitoring by WMO [47]. SPI is easy
to compute, can be calculated for different ranges of time (1, 3, 6, 13, 24
months), capture multitemporal nature of rainfall deficiency and uses readily
available precipitation data [28]. Typical values of SPI values during dry and
wet seasons is shown in Figure 3.

Figure 3: Spatial variation of SPI index over West Africa in (a) February 1965
and (b) June 1965 using a scale of 12 months.

###  2.4 Linear correlation

The interrelation between any climate index $x$ and SPI (represented by $y$)
such that a change in one variable results in a positive or negative change in
the other variable can be investigated using the correlation coefficient
$\rho$. The Pearson method of computing $\rho$ is described as

$\rho=\frac{\sum_{i=1}^{n}(x_{i}-\overline{x})(y_{i}-\overline{y})}{\sqrt{\sum_{i=1}^{n}(x_{i}-\overline{x})^{2}\sum_{i=1}^{n}(y_{i}-\overline{y})^{2}}}$
(2)

The values of $\rho$ in the ranges $-1\leq\rho\leq 0$ and $0\leq\rho\leq 1$
represents positive and negative correlation respectively. If $\rho=0$, the
two time series are said to be uncorrelated [48].

###  2.5 Co-integration

Correlation (positive or negative) does not infer causality between two time
series. To determine causality between two time series the method of
cointegration is used. If the linear combination of two non-stationary time
series, $x$ and $y$ are stationary after differencing, we say, the time series
are cointegrated. In this work, the Engle-Granger approach will be used [49].
The two time series are tested for stationarity, if not, they are made
stationary by differencing with order 1 (ie I(1). Ordinary least square is
used to investigate the relationship $y_{t}=ax_{t}+e_{t}$. The residual
$e_{t}$ is checked for stationarity. In this work, the linear relationship
with constant trend is used. Stationarity is tested using the Augumented
Dickey-Fuller test.

###  2.6 Mutual Information and sychronization

Mutual information (I) quantifies linear and nonlinear interdependence between
two systems and quantity of information shared [50]. It is computed as

$I(X;Y)=\sum_{x,y}P_{xy}(x,y)\log\frac{P_{xy(x,y)}}{P_{x}(x)P_{y}(y)}$ (3)

where $P_{x,y}$ and $P(x)$ are the joint and marginal distribution
respectively.

An interdepence measure, H, was introduced by Quiroga et al. [51, 52] which is
easier to compute, robust to noise, and more sensitive to weak dependencies
than mutual information. It is expressed as

$H^{(k)}(X|Y)=\frac{1}{N}\sum_{n=1}^{N}\ln\frac{R_{n}(X)}{R^{(k)}_{n}(X|Y)}$
(4)

$H$ is zero if X and Y are completely independent. A new measure similar to
$H$ but computed with a different averaging method is used in this study.

$N^{(k)}(X|Y)=\frac{1}{N}\sum_{n=1}^{N}\ln\frac{R_{n}(X)-R^{(k)}_{n}(X|Y)}{R_{n}(X)}$
(5)

$N$ and $H$ have positive values if nearness in Y implies also nearness in X
for equal time partner while negative values are obtained if close pairs in Y
correspond to distant pairs in X [51, 52].

##  3 Results and Discussion

Figure 4 shows the correlation coefficient and its significance, between the
teleconnection patterns considered and SPI drought for the West African
region. NAO shows a predominantly small but significant correlation with SPI
drought index. Positive but non-significant correlation could be observed
along the Eastern coast of West Africa. The Pacific Decadal Oscillation showed
similar but stronger results over the same region (Figure 4e). However,
negative trends were not found along the Eastern Coast as in the case of NAO.
The Southern Oscillation Index showed the opposite trends observed in both NAO
and PDO with Significant negative correlation distributed over the region.

Figure 4: Zero lag correlation coefficient between SPI drought index and (a)
NAO, (c) SOI and (e) PDO and corresponding significance for (b) NAO, (d) SOI
and (f) PDO

The result of the Engle-Granger cointegration test between large scale
oscillations (PDO, NAO and SOI) and SPI index is shown in Figure 5. Both NAO
and PDO exhibit entirely negative linear relationship over the region. In the
linear relationship between NAO and SPI, the eastern seaboard of West Africa
showed the strongest negative correlation. However, although SOI showed
patches of negative relationship, the region boasts of predominantly positive
linear relationship.

The nonlinear relationships between the large scale oscillations (PDO, NAO,
SOI) and SPI drought index over West Africa were investigated using mutual
information and phase synchronization. NAO showed weak mutual information with
SPI in sparse locations (Figure 6). This implies that not much information can
be obtained concerning SPI by knowing the values of NAO. Similar patterns
could be observed for SOI in the region. Large amount of information could,
however, be obtained concerning SPI by having PDO values. This is much more
significant in the Sahel region of West Africa.

Nonlinear dependence measure, $N$ was found to have predominantly negative
values between NAO and SPI values (7). Positive values could be observed off
the coast of Nigeria to Liberia which indicate nonlinear dependence. SPI also
showed positive nonlinear dependence with SOI along the coast of West Africa.
PDO showed patches of nonlinear dependence in the region.

Comparing the five different approaches to quantifying the relationship
between drought in the West African region and large scale oscillation in both
the Pacific and Atlantic Ocean, the mechanism of influence could be deduced.
It has been posited that sea surface temperature contribute significantly to
drought regimes in West Africa but the specific contributions of major
teleconnection patterns have not been investigated [53].

Regions with highest negative correlation (dark blue) were found to correspond
to regions with lowest co-integration (yellow).

Figure 5: Engel-Granger co-integration test between SPI drought index and (a)
NAO, (b) SOI and (c) PDO
Figure 6: Mutual Information between SPI drought index and (a) NAO, (b) SOI
and (c) PDO
Figure 7: Nonlinear dependence measure, $N$, between SPI drought index and (a)
NAO, (b) SOI and (c) PDO

##  4 Conclusion

Several factors such as land-atmosphere interaction, sea surface temperature
and greenhouse gases have been identified as drivers of drought in the West
African region [53]. These factors have spatio-temporal intensity and
dominance for different parts of the region. In this paper, the linear and
nonlinear dependency between three large scale oscillation and drought index
(SPI) in the West African region has been investigated. This was achieved
through the use of linear analysis (correlation coefficient) and nonlinear
analysis (phase synchronization, nonlinear dependence test, and mutual
information).

## References

  * Zhang et al. [2007] K. Zhang, J. S. Kimball, K. C. McDonald, J. J. Cassano, S. W. Running, Impacts of large-scale oscillations on pan-Arctic terrestrial net primary production, Geophysical Research Letters 34 (21) (2007) n/a–n/a, ISSN 1944-8007, URL http://dx.doi.org/10.1029/2007GL031605, l21403.
  * Jones et al. [1997] P. D. Jones, Jónsson, D. Wheeler, Extension to the North Atlantic Oscillation using early instrumental pressure observations from Gibraltar and South-West Iceland., International Journal of Climatology 17 (1997) 1433–1450, ISSN 1097-0088, URL http://dx.doi.org/10.1002/joc.3711.
  * Brandimarte et al. [2011] L. Brandimarte, G. Di Baldassarre, G. Bruni, P. D’Odorico, A. Montanari, Relation Between the North-Atlantic Oscillation and Hydroclimatic Conditions in Mediterranean Areas, Water Resources Management 25 (5) (2011) 1269–1279.
  * Moulin et al. [1997] C. Moulin, C. E. Lambert, F. Dulac, U. Dayan, Control of atmospheric export of dust from North Africa by the North Atlantic Oscillation, Nature 387 (6634) (1997) 691\.
  * McHugh and Rogers [2001] M. J. McHugh, J. C. Rogers, North Atlantic Oscillation Influence on Precipitation Variability around the Southeast African Convergence Zone, Journal of Climate 14 (17) (2001) 3631–3642, URL https://doi.org/10.1175/1520-0442(2001)014<3631:NAOIOP>2.0.CO;2.
  * Rossello [2011] J. Rossello, North Atlantic Oscillation influences on European airline traffic, Transportation Research Part D 16 (2011) 183–187, URL https://doi.org/10.1016/j.trd.2010.09.003.
  * Báez et al. [2014] J. C. Báez, R. Real, V. López-Rodas, E. Costas, A. E. Salvo, C. García-Soto, A. Flores-Moya, The North Atlantic Oscillation and the Arctic Oscillation favour harmful algal blooms in SW Europe, Harmful Algae 39 (2014) 121–126, URL https://doi.org/10.1016/j.hal.2014.07.008.
  * ANDERS and POST [2006] A. D. ANDERS, E. POST, Distribution-wide effects of climate on population densities of a declining migratory landbird, Journal of Animal Ecology 75 (1) (2006) 221–227, ISSN 1365-2656, URL http://dx.doi.org/10.1111/j.1365-2656.2006.01034.x.
  * Adeniyi and Dilau [2018] M. O. Adeniyi, K. A. Dilau, Assessing the link between Atlantic Niño 1 and drought over West Africa using CORDEX regional climate models, Theoretical and Applied Climatology 131 (3-4) (2018) 937–949, URL http://link.springer.com/10.1007/s00704-016-2018-0.
  * Anyamba et al. [2001] A. Anyamba, C. Tucker, J. Eastman, NDVI anomaly patterns over Africa during the 1997/98 ENSO warm event, International journal of remote sensing 22 (10) (2001) 1847–1860.
  * Anyamba and Eastman [1996] A. Anyamba, J. Eastman, Interannual variability of NDVI over Africa and its relation to El Niño/Southern Oscillation, International Journal of Remote Sensing 17 (1996) 2533–2548.
  * Huang et al. [1998] J. Huang, K. Higuchi, A. Shabbar, The relationship between the North Atlantic Oscillation and El Ni?o-Southern Oscillation, Geophysical Research Letters 25 (14) (1998) 2707–2710, ISSN 1944-8007, URL http://dx.doi.org/10.1029/98GL01936.
  * Rajagopalan et al. [1998] B. Rajagopalan, Y. Kushnir, Y. M. Tourre, Observed decadal midlatitude and tropical Atlantic climate variability, Geophysical Research Letters 25 (21) (1998) 3967–3970, ISSN 1944-8007, URL http://dx.doi.org/10.1029/1998GL900065.
  * Wu et al. [2012] Z. Wu, J. Li, Z. Jiang, J. He, X. Zhu, Possible effects of the North Atlantic Oscillation on the strengthening relationship between the East Asian Summer monsoon and ENSO, International Journal of Climatology 32 (5) (2012) 794–800, ISSN 1097-0088, URL http://dx.doi.org/10.1002/joc.2309.
  * Gimeno et al. [2002] L. Gimeno, P. Ribera, R. Iglesias, L. de la Torre, R. García, E. Hernández, Identification of empirical relationships between indices of ENSO and NAO and agricultural yields in Spain, Climate Research 21 (2) (2002) 165–172.
  * Bonaccorso et al. [2015] B. Bonaccorso, A. Cancelliere, G. Rossi, Probabilistic forecasting of drought class transitions in Sicily (Italy) using Standardized Precipitation Index and North Atlantic Oscillation Index, Journal of Hydrology 526 (2015) 136–150, URL https://doi.org/10.1016/j.jhydrol.2015.01.070.
  * McBride and Nicholls [1983] J. L. McBride, N. Nicholls, Seasonal Relationships between Australian Rainfall and the Southern Oscillation, Monthly Weather Review 111 (10) (1983) 1998–2004, URL https://doi.org/10.1175/1520-0493(1983)111<1998:SRBARA>2.0.CO;2.
  * Liu et al. [2018] Z. Liu, X. Zhang, R. Fang, Multi-scale linkages of winter drought variability to ENSO and the Arctic Oscillation: A case study in Shaanxi, North China, Atmospheric Research 200 (March 2016) (2018) 117–125, URL http://dx.doi.org/10.1016/j.atmosres.2017.10.012.
  * Jin et al. [2005] Y. H. Jin, A. Kawamura, K. Jinno, R. Berndtsson, Quantitative relationship between SOI and observed precipitation in southern Korea and Japan by nonparametric approaches, Journal of Hydrology 301 (1-4) (2005) 54–65, URL https://doi.org/10.1016/j.jhydrol.2004.06.026.
  * Bhuvaneswari et al. [2013] K. Bhuvaneswari, V. Geethalakshmi, A. Lakshmanan, R. Srinivasan, N. U. Sekhar, The Impact of El Niño/Southern Oscillation on Hydrology and Rice Productivity in the Cauvery Basin, India: Application of the Soil and Water Assessment Tool, Weather and Climate Extremes 2 (2013) 39–47, URL https://doi.org/10.1016/j.wace.2013.10.003.
  * Ranasinghe et al. [2004] R. Ranasinghe, R. McLoughlin, A. Short, G. Symonds, The Southern Oscillation Index, wave climate, and beach rotation, Marine Geology 204 (3-4) (2004) 273–287, URL https://doi.org/10.1016/S0025-3227(04)00002-7.
  * Zhong et al. [2017] K. Zhong, F. Zheng, H. Wu, C. Qin, X. Xu, Dynamic changes in temperature extremes and their association with atmospheric circulation patterns in the Songhua River Basin, China, Atmospheric Research 190 (2016) (2017) 77–88, URL http://dx.doi.org/10.1016/j.atmosres.2017.02.012.
  * Mantua [2002] N. J. Mantua, Pacific ? Decadal Oscillation ( PDO ), in: Encyclopedia of Global Environmental Change, vol. 1, John Wiley & Sons, Ltd, ISBN 0471977969, 592–594, 2002\.
  * Henley [2017] B. J. Henley, Pacific decadal climate variability: Indices, patterns and tropical-extratropical interactions, Global and Planetary Change 155 (2017) 42 – 55, URL http://www.sciencedirect.com/science/article/pii/S0921818116304398.
  * Zhou et al. [2015] X. Zhou, Y. Sun, W. Huang, J. P. Smol, Q. Tang, L. Sun, The Pacific decadal oscillation and changes in anchovy populations in the Northwest Pacific, Journal of Asian Earth Sciences 114 (2015) 504–511, URL https://doi.org/10.1016/j.jseaes.2015.06.027.
  * Shi et al. [2018] J. Shi, L. Cui, K. Wen, Z. Tian, P. Wei, B. Zhang, Trends in the consecutive days of temperature and precipitation extremes in China during 1961?2015, Environmental Research 161 (July 2017) (2018) 381–391, URL https://doi.org/10.1016/j.envres.2017.11.037.
  * Park et al. [2017] J. Park, R. Byrne, H. Böhnel, The combined influence of Pacific decadal oscillation and Atlantic multidecadal oscillation on central Mexico since the early 1600s, Earth and Planetary Science Letters 464 (2017) 1–9, URL http://dx.doi.org/10.1016/j.epsl.2017.02.013.
  * Shi et al. [2017] P. Shi, T. Yang, C. Y. Xu, B. Yong, Q. Shao, Z. Li, X. Wang, X. Zhou, S. Li, How do the multiple large-scale climate oscillations trigger extreme precipitation?, Global and Planetary Change 157 (April) (2017) 48–58, URL http://dx.doi.org/10.1016/j.gloplacha.2017.08.014.
  * Okeke et al. [2006] F. N. Okeke, J. A. Marengo, C. Nobre, Investigation of El Nino and La Nina effects and the impact of Atlantic sea surface temperatures (SSTs) on precipitation in Nigeria from 1950 to 1992, Surv. Geophys 27 (2006) 545–555, URL 10.1007/s10712-006-9009-0.
  * Oguntunde et al. [2017] P. G. Oguntunde, B. J. Abiodun, G. Lischeid, Impacts of climate change on hydro-meteorological drought over the Volta Basin, West Africa, Global and Planetary Change 155 (March) (2017) 121–132, URL http://dx.doi.org/10.1016/j.gloplacha.2017.07.003.
  * Nicholson [1981] S. E. Nicholson, Rainfall and atmospheric circulation during drought periods and wetter years in West Africa, Monthly Weather Review 109 (10) (1981) 2191\.
  * Giannini et al. [2003] A. Giannini, R. Saravanan, P. Chang, Oceanic Forcing of Sahel Rainfall on Interannual to Interdecadal Time Scales, Science 302 (5647) (2003) 1027–1030, URL http://science.sciencemag.org/content/302/5647/1027.
  * Ahn and Kim [2005] J. H. Ahn, H. S. Kim, Nonlinear Modeling of El Nino/Southern Oscillation Index, Hydrological Engineering 10 (2005) 8–15.
  * Fuwape et al. [2017] I. Fuwape, S. Ogunjo, S. Oluyamo, A. Rabiu, Spatial variation of deterministic chaos in mean daily temperature and rainfall over Nigeria, Theoretical and Applied Climatology 130 (1-2).
  * Sivakumar [2004] B. Sivakumar, Chaos theory in geophysics: past, present and future, Chaos Solitons & Fractals 19 (2) (2004) 441–462, URL https://doi.org/10.1016/s0960-0779(03)00055-9.
  * Fuwape and Ogunjo [2016] I. A. Fuwape, S. T. Ogunjo, Quantification of scaling exponents and dynamical complexity of microwave refractivity in a tropical climate, Journal of Atmospheric and Solar-Terrestrial Physics 150-151 (2016) 61–68, URL https://doi.org/10.1016/j.jastp.2016.10.010.
  * Fuwape et al. [2016] I. A. Fuwape, S. T. Ogunjo, J. B. Dada, G. A. Ashidi, I. Emmanuel, Phase synchronization between tropospheric radio refractivity and rainfall amount in a tropical region, Journal of Atmospheric and Solar-Terrestrial Physics 149 (2016) 46–51, URL http://linkinghub.elsevier.com/retrieve/pii/S1364682616302802.
  * Builes-Jaramillo et al. [2017] A. Builes-Jaramillo, N. Marwan, G. Poveda, J. Kurths, Nonlinear interactions between the Amazon River basin and the Tropical North Atlantic at interannual timescales, Climate Dynamics (2017) 1–19URL http://link.springer.com/10.1007/s00382-017-3785-8.
  * Ogunjo [2015] S. Ogunjo, Effect of Data Transformation on Long Term Memory of Chaotic Time Series, The African Review of Physics 10 (2015) 219–224.
  * Omotosho [1990] J. B. Omotosho, Onset of thunderstorms and precipitation over northern nigeria, International Journal of Climatology 10 (1990) 849–860.
  * Vondou et al. [2010] D. A. Vondou, A. Nzeukou, A. Lenouo, F. Mkankam Kamga, Seasonal variations in the diurnal patterns of convection in Cameroon-Nigeria and their neighboring areas, Atmospheric Science Letters 11 (4) (2010) 290–300, URL https://doi.org/10.1002/asl.297.
  * Joseph O. Adejuwon [2012] Joseph O. Adejuwon, Rainfall seasonality in the Niger Delta Belt, Nigeria, URL https://doi.org/10.5897/JGRP11.096, 2012.
  * Harris et al. [2014] I. Harris, P. Jones, T. Osborn, D. Lister, Updated high-resolution grids of monthly climatic observations ? the CRU TS3.10 Dataset, International Journal of Climatology 34 (3) (2014) 623–642, ISSN 1097-0088, URL http://dx.doi.org/10.1002/joc.3711.
  * McKee et al. [1993] T. McKee, N. Doesken, J. Kleist, The relationship of drought frequency and duration to time scales, Proceedings of the Eight Conference on Applied Climatology (1993) 179–184.
  * Costa [2011] A. C. Costa, Local patterns and trends of the standard precipitation index in southern Portugal (1940-1999), Advances in Geosciences 30 (2011) 11–16.
  * Oloruntade et al. [2017] A. J. Oloruntade, T. A. Mohammad, A. H. Ghazali, A. Wayayok, Analysis of meteorological and hydrological droughts in the Niger-South Basin, Nigeria, Global and Planetary Change 155 (May) (2017) 225–233, URL http://dx.doi.org/10.1016/j.gloplacha.2017.05.002.
  * Svoboda et al. [2012] M. Svoboda, M. Hayes, D. Wood, Standardized precipitation index user guide, World Meteorological Organization Geneva, Switzerland .
  * Arora et al. [2014] P. N. Arora, S. Arora, S. Arora, A. Arora, Comprehensive Statistical Methods, S. Chand and Company Pvt, Ltd, 4th edn., 2014\.
  * Engle and Granger [1987] R. F. Engle, C. W. Granger, Co-integration and error correction: representation, estimation, and testing, Econometrica: journal of the Econometric Society (1987) 251–276.
  * Baptista et al. [2012] M. S. Baptista, R. M. Rubinger, E. R. Viana, J. C. Sartorelli, U. Parlitz, C. Grebogi, Mutual Information Rate and Bounds for It, PLoS ONE 7 (10) (2012) 1–11.
  * Quiroga et al. [2000a] R. Q. Quiroga, J. Arnhold, P. Grassberger, Learning driver-response relationships from synchronization patterns, Physical Review E 61 (5) (2000a) 5142–5148.
  * Quiroga et al. [2000b] R. Q. Quiroga, A. Kraskov, T. Kreuz, P. Grassberger, Performance of different synchronization measures in real data: A case study of electroencephalogrpahic signals, Physical Review E 65 (2000b) 041903.
  * VIS [2019] Chapter 4 - Precipitation Extremes in the West African Sahel: Recent Evolution and Physical Mechanisms, in: V. Venugopal, J. Sukhatme, R. Murtugudde, R. Roca (Eds.), Tropical Extremes, Elsevier, 95 – 138, URL http://www.sciencedirect.com/science/article/pii/B9780128092484000042, 2019\.


11institutetext: Department of Electrical and Computer Engineering
22institutetext: Department of Mechanical Engineering
The University of New Mexico, Albuquerque NM 87131, USA.
22email<EMAIL_ADDRESS>33institutetext:
Department of Computer Science & Engineering
Lehigh University, Bethlehem PA 18015, USA.
33email<EMAIL_ADDRESS>
# Decentralized Adaptive Aerospace Transportation of Unknown Loads Using A
Team of Robots

Longsen Gao 11 Kevin Aubert 11 David Saldaña 33 Claus Danielson 22
Rafael Fierro 11

###### Abstract

Transportation missions in aerospace are limited to the capability of each
aerospace robot and the properties of the target transported object, such as
mass, inertia, and grasping locations. We present a novel decentralized
adaptive controller design for multiple robots that can be implemented in
different kinds of aerospace robots. Our controller adapts to unknown objects
in different gravity environments. We validate our method in an aerial
scenario using multiple fully actuated hexarotors with grasping capabilities,
and a space scenario using a group of space tugs. In both scenarios, the
robots transport a payload cooperatively through desired three-dimensional
trajectories. We show that our method can adapt to unexpected changes that
include the loss of robots during the transportation mission.
\hrefhttp://tiny.cc/dars2024unmhttp://tiny.cc/dars2024unm

###### keywords:

Multi-robot System, Adaptive Control, Aerospace Robot

##  1 Introduction

The payload capacity of an aerospace robot, such as an Unmanned Aerial Vehicle
(UAV) or a space tug, presents a fundamental constraint that limits its
utility across various domains, notably in the construction and transportation
of substantial payloads. Augmenting the payload capacity of an aerospace robot
involves intricate mechanical redesigns, which are often resource-intensive.
An alternative approach that leverages the collaborative efforts of multiple
aerospace robots for transporting and manipulating target payload emerges as a
cost-effective and promising solution. Multiple aerospace robots for
cooperative manipulation and transportation considering the dynamic coupling
between the individual subsystems have been widely studied in both space [1,
2, 3, 4] and aerial environments [5, 6] in recent decade.

The domain of cooperative aerial transportation facilitated through tethers,
has been the subject of extensive research [7, 8, 9, 10]. These investigations
underscore the potential of leveraging contact forces in cooperative
transportation tasks, highlighting the significance of designing aerial
vehicles capable of synchronized complex maneuvers. However, manipulating a
tether-suspended object using multiple aerospace robots requires considering
the feasibility of each robot’s position with the contact force direction to
avoid damage to the whole system, increasing the complexity of internal
dynamics analysis. In this case, recent research works in [11, 12, 13, 14]
learn rigid attachments for transportation purposes using a group of drones
that simplify internal dynamics analysis and improve maneuverability. However,
the rigid attachment scheme in those studies assumes the payload and aerospace
robots are firmly connected through a non-detachable mechanic structure,
decreasing the aerospace robots’ dexterity. The work in [15] studies a group
of drones installed grippers in cooperative payload stabilization task. The
gripper aims to penetrate the surface to hold the object, introducing the
hidden risk of damage and limiting the scope of the object material. An
inflight self-disassembly gripper design in [16] uses a self-adaptation scheme
by self-reconfiguring during the transportation mission process to grasp the
object. Still, the structural design of each module limits the application of
the whole system to the shape of the target operating object, especially for
the thin panel which is hard to hold from a flat plane during the
transportation mission. In [17], a multi-drone system with a self-adaptive
gripper is developed for grasping and transportation tasks. In this paper, we
choose to install a rigid rod underneath the drone’s base with a gripper on
its end as shown in Fig.1. This structure design not only introduces the
control complexity for the whole dynamics of the vehicle but also increases
the detachability for different shapes of the object. Additionally, all
aforementioned works limit their drone’s structure in the directions of all
propellers along $z$-axis in its body frame, which can only generate the
wrench
$\mathbf{w}=\left[f_{z},\boldsymbol{\tau}\right]^{\intercal}\in\mathbb{R}^{4}$
in four dimensions where
$\boldsymbol{\tau}=\left[\tau_{x},\tau_{y},\tau_{z}\right]^{\intercal}\in\mathbb{R}^{3}$.
The limited structure lost $f_{x}$ and $f_{y}$ generated by the 4 motors of
the quadrotor, which would also limit the performance of the transportation
mission in 3D space. In this paper, we consider using a fully actuated
hexarotor UAV [18, 19] in which each rotor has been tilted to a specific angle
along 2 axes to give the drone full actuation in 6 DoF.

Manipulating and interacting with partially known environments are challenging
as they involve adapting to different operational conditions, especially for
missions using multiple aerospace robots in [20, 21, 22]. One of the main
drawbacks of the common methods is their lack of adaptability when the
object’s parameters to service or transport are unknown [23]. Similarly, robot
manipulators performing a contact-based inspection of a variable friction
surface would have difficulty maintaining contact with the surface during the
inspection task [24]. In this case, we design a novel decentralized adaptive
controller in paper to implement on two different aerospace-manipulated
systems: one composed of fully actuated hexarotor UAVs and another one
composed of a group of space tugs, in the same transportation tasks under
gravity and zero-gravity environment to examine the feasibility of our
controller design.

The contribution of this paper is threefold: (i) We combine a grasping system
and a fully actuated hexarotor UAV to output desired wrench
$\boldsymbol{w}\in\mathbb{R}^{6}$ in 6 DoF on grasping system to hold and
transport the object, which improves the detachability and actuation
performance during our transportation task in 3D space, and (ii) introduce a
decentralized adaptive method implemented on multiple aerospace systems that
operate effectively under uncertainties, including unknown mass, inertia and
the location of grasping points, and (iii) we extend our work on multiple
space tugs under a zero-gravity environment by testing the feasibility of our
decentralized adaptive controller design. This method represents a robust
solution to the longstanding challenges associated with satellite detumbling.

##  2 Decentralized Adaptive Transportation Problem

This section describes the transportation problem for a free-floating rigid
panel in the shape of a convex polygon by implementing our decentralized
adaptive controller design, including the dynamics of the two different
aerospace systems and the grasping dynamics from each agent to the payload.

###  2.1 Robots and Reference Frames

This paper proposes a novel methodology for transporting a rigid polygon panel
using multiple drones and then extends the work by using a group of space tugs
in space. We consider a team of $n$ vehicles, indexed by $i=1,...,n$, that can
be either drones or space tugs.

Figure 1: Fully actuated hexarotor UAV with a grasping mechanism.

######  Definition 2.1 (Drone).

A drone is a fully actuated hexarotor [19, 18] with a grasping mechanism
composed of a rigid rod and a gripper as shown in Fig. 1. Once the drone
grasps the object, the connection is considered rigid.

It comprises six rotors, each rotor $j$ generates a wrench denoted as
$\mathbf{w}_{i_{j}}=\left[\boldsymbol{f}_{i_{j}}^{\intercal},\boldsymbol{\tau}_{i_{j}}^{\intercal}\right]^{\intercal}$
$\in\mathbb{R}^{6}$ where
$\boldsymbol{f}_{i_{j}}=\left[0,0,k_{f}\,\omega_{i_{j}}^{2}\right]^{\intercal}$
is the force and
$\boldsymbol{\tau}_{i_{j}}=\left[0,0,k_{m}\,\omega_{i_{j}}^{2}\right]^{\intercal}$
is the torque that the rotor generates when it rotates at an angular velocity
$\omega_{i_{j}}$. The coefficients $k_{f}$ and $k_{m}$ can be obtained
experimentally. The rotors are tilted as pictured in Fig. 1.

######  Definition 2.2 (Space Tug).

A space tug is an independently operating, self-propelled spacecraft for
assembly, maintenance, repair, and contingency operations in $SE(3)$ [2] (see
Fig. 2).

Figure 2: Space tug with 6 rockets around its sides. The red arrow shows each
rocket propulsion in its direction.

The tug is propelled using 6 rockets in different directions to push or pull
in 3D space, forming a fully actuated vehicle that generates a wrench
$\mathbf{w}_{i}\in\mathbb{R}^{6}$. Each vehicle $i$, either drone or space
tug, has a mass $m_{i}$ and moment-of-inertia $\boldsymbol{J}_{i}$. Note that
we will use drones and space tugs to process the same transportation task at
different gravity levels to prove the feasibility of our controller design.

The coordinate systems are illustrated in Fig. 3. The world frame
$\mathcal{W}$ is fixed and has its $z$-axis pointing upwards. We consider $n$
vehicles using grippers to grasp a rigid object. The body frame $\mathcal{B}$
is located on the object’s center of mass. It is assumed that the body frame
axes are chosen as the principal axes of the entire system. Each vehicle has
an individual body frame $\mathcal{U}_{i}$, attached to its center of mass
with $z_{\mathcal{U}_{i}}$ perpendicular to the plane of the base and pointing
up. We require the $z_{\mathcal{U}_{i}}$ axes and $z_{\mathcal{B}}$ to be
parallel. We use Row-Pitch-Yaw to model the rotation of frame $\mathcal{B}$
and frame $\mathcal{U}_{i}$ in the world frame $\mathcal{W}$. We set up one
estimated point $\boldsymbol{p}_{s}$ on the target object that each vehicle
can share their common measurement information $\boldsymbol{d}_{i}$ which is
from each grasping point to $\boldsymbol{p}_{s}$ to estimate the object’s
unknown properties concerning the unknown information $\boldsymbol{\bar{d}}$
which is from $\boldsymbol{p}_{s}$ to the origin of the frame $\mathcal{B}$.
Let $\boldsymbol{q}\in\mathbb{R}^{6}$ be the position and orientation for the
manipulated object.
$\boldsymbol{q}=[\boldsymbol{p}^{\intercal},\boldsymbol{\theta}^{\intercal}]^{\intercal}\in\mathbb{R}^{6\times
1}$, where $\boldsymbol{p}\in\mathbb{R}^{3}$ and
$\boldsymbol{\theta}\in\mathbb{R}^{3}$ denotes the pose of $\mathcal{B}$ with
respect to $\mathcal{W}$. The rotation matrix
$\mathbf{{}^{\mathcal{B}}_{\mathcal{W}}R}\in SO(3)$ denotes the orientation of
$\mathcal{B}$ with respect to $\mathcal{W}$. We use
$\mathbf{{}^{\mathcal{U}_{i}}_{\mathcal{W}}R}$ to denote the rotation matrix
of $\mathcal{U}_{i}$ with respect to $\mathcal{W}$.

###  2.2 Dynamics

Figure 3: A group of fully actuated hexarotor UAVs grasping a free-floating
object during the transportation task in 3D space.

The dynamics of the $n$ vehicles attached to the object are modeled using the
Euler-Largrian equation,

$\displaystyle\mathbf{w}=\boldsymbol{M}(\boldsymbol{q})\ddot{\boldsymbol{q}}+\boldsymbol{C}(\boldsymbol{q},\dot{\boldsymbol{q}})\dot{\boldsymbol{q}}+\boldsymbol{g},$
(1)

where
$\mathbf{w}=\begin{bmatrix}\boldsymbol{f},\boldsymbol{\tau}\end{bmatrix}^{\intercal}\in\mathbb{R}^{6\times
1}$ is the total wrench generated by the $n$ vehicles, where $\boldsymbol{f}$
is the total force, and $\boldsymbol{\tau}$ is the total torque with respect
to the point ${\boldsymbol{p}}_{s}$ in the frame $\mathcal{B}$. The positive
semi-definite inertia matrix of the object is

$\displaystyle\boldsymbol{M}(\boldsymbol{q})=\begin{bmatrix}m\boldsymbol{\mathrm{I}}&-m(\mathbf{{}^{\mathcal{B}}_{\mathcal{W}}R}\boldsymbol{\bar{d}})^{\times}\\\
-m(\mathbf{{}^{\mathcal{B}}_{\mathcal{W}}R}\boldsymbol{\bar{d}})^{\times}&\mathbf{{}^{\mathcal{B}}_{\mathcal{W}}R}\boldsymbol{J_{d}}\mathbf{{}^{\mathcal{B}}_{\mathcal{W}}R^{\intercal}}\end{bmatrix},$

where $\boldsymbol{\mathrm{I}}\in\mathbb{R}^{3\times 3}$ is identity matrix,
$\boldsymbol{\bar{d}}\in\mathbb{R}^{3}$ and $\boldsymbol{J}_{d}$ are
respectively the position vector of the estimated point concerning the center-
of-mass and moment-of-inertia of the client in frame $\mathcal{B}$.
$(\cdot)^{\times}$ denotes the map from $\mathbb{R}^{3}$ to $SO(3)$ which is
skew-symmetric [25]. The moment of inertia is obtained using the parallel axis
theorem[26],

$\displaystyle\boldsymbol{J_{d}}=\boldsymbol{J_{cm}}+\boldsymbol{m}((\boldsymbol{\bar{d}^{\intercal}\bar{d}})\boldsymbol{\mathrm{I}}-\boldsymbol{\bar{d}\bar{d}^{\intercal}}),$

where $\boldsymbol{J_{cm}}$ is the moment-of-inertia about the center-of-mass.
The centripetal and Coriolis matrix is

$\leavevmode\resizebox{432.17372pt}{}{$\boldsymbol{C}(\boldsymbol{q},\boldsymbol{\dot{q}})=\begin{bmatrix}\mathbf{0}_{3\times
3}&-m\boldsymbol{\omega}^{\times}\left(\mathbf{\mathbf{{}^{\mathcal{B}}_{\mathcal{W}}R}\bar{d}}\right)^{\times}\\\
-m\boldsymbol{\omega}^{\times}\left(\mathbf{\mathbf{{}^{\mathcal{B}}_{\mathcal{W}}R}\bar{d}}\right)^{\times}&\boldsymbol{\omega}^{\times}\mathbf{{}^{\mathcal{B}}_{\mathcal{W}}R}\boldsymbol{J}_{d}\mathbf{{}^{\mathcal{B}}_{\mathcal{W}}R}-m\left(\left(\mathbf{{}^{\mathcal{B}}_{\mathcal{W}}R}\mathbf{\bar{d}}\right)^{\times}\dot{\mathbf{q}}\right)^{\times}\end{bmatrix}$}.$

The gravity term $\boldsymbol{g}$ is a $6\times 1$ vector of the gravitational
wrench. Note that this term only applies to the scenario using aerial
vehicles. Otherwise, this term is assumed to be zero for space tugs.

###  2.3 Vehicle’s Wrench

The rotors in the vehicle generate a wrench
$\mathbf{w}_{p_{i}}=\left[\boldsymbol{f}_{p_{i}},\boldsymbol{\tau}_{p_{i}}\right]^{\intercal}$
with respect to the the center of mass of the $i^{\text{th}}$ vehicle in frame
$\mathcal{U}_{i}$, in which $\boldsymbol{f}_{p_{i}}$ and
$\boldsymbol{\tau}_{p_{i}}$ are defined as

$\displaystyle\boldsymbol{f}_{p_{i}}$
$\displaystyle=\sum_{j=0}^{6}\mathbf{{}^{\mathcal{M}_{i_{j}}}_{\mathcal{U}_{i}}R}\boldsymbol{f}_{i_{j}}$
(2) $\displaystyle\boldsymbol{\tau}_{p_{i}}$
$\displaystyle=\sum_{j=0}^{6}\left(\mathbf{{}^{\mathcal{M}_{i_{j}}}_{\mathcal{U}_{i}}R}\boldsymbol{f}_{i_{j}}\right)^{\times}\boldsymbol{d}_{i_{j}}+\sum_{j=0}^{6}\mathbf{{}^{\mathcal{M}_{i_{j}}}_{\mathcal{U}_{i}}R}\boldsymbol{\tau}_{i_{j}},$

where $\boldsymbol{d}_{i_{j}}\in\mathbb{R}^{3}$ denotes the position of
$j^{\text{th}}$ motor in frame $\mathcal{U}_{i}$. The rotation matrix
${}^{\mathcal{M}_{i_{j}}}_{\mathcal{U}_{i}}\mathbf{R}\in SO(3)$ defines the
rotor’s orientation with respect to $\mathcal{U}_{i}$. It is assumed the
rotors are in different orientations to form a fully actuated vehicle
following the conditions in [19].

The wrench
$\mathbf{w}_{i}=\left[\boldsymbol{f}_{i},\boldsymbol{\tau}_{i}\right]$
generated by $i^{\text{th}}$ vehicle in frame $\mathcal{W}$ can be calculated
using the Newton-Euler’s equation,

$\displaystyle\boldsymbol{f}_{i}$
$\displaystyle=\mathbf{{}^{\mathcal{U}_{i}}_{\mathcal{W}}R}\boldsymbol{f}_{p_{i}}-m_{i}\dot{\boldsymbol{v}}_{i}+\boldsymbol{g}_{i}$
(3) $\displaystyle\boldsymbol{\tau}_{i}$
$\displaystyle=\mathbf{{}^{\mathcal{U}_{i}}_{\mathcal{W}}R}(-\left(\boldsymbol{J}_{i}\boldsymbol{\omega}_{i}\right)^{\times}\boldsymbol{\omega}_{i}+\boldsymbol{\tau}_{p_{i}}-\boldsymbol{J}_{i}\dot{\boldsymbol{\omega}}_{i}),$

where $\boldsymbol{g}_{i}$ is the vector of the gravitational force of the
$i^{\text{th}}$ vehicle, $\boldsymbol{v}_{i}\in\mathbb{R}^{3}$ denote the
linear velocity of the $i^{\text{th}}$ vehicle in frame $\mathcal{W}$, and
$\boldsymbol{\omega}_{i}\in\mathbb{R}^{3}$ represent the angular velocity of
the $i^{\text{th}}$ vehicle in frame $\mathcal{U}_{i}$.

###  2.4 Grasping Dynamics of Manipulator System

Next, we model the wrench $\mathbf{w}$ applied to the target object by an
aerospace robot. The total wrench applied on the object is the sum of the
wrenches $\mathbf{w}_{i}$ applied by each robot

$\displaystyle\mathbf{w}=\sum_{i=1}^{N}\boldsymbol{G}\left(\boldsymbol{q},\boldsymbol{d}_{i}\right)\mathbf{w}_{i},$
(4)

where the matrix $\boldsymbol{G}$ maps the forces and torques produced by each
robot into wrenches on the object

$\displaystyle\boldsymbol{G}\left(\boldsymbol{q},\boldsymbol{d}_{i}\right)=\begin{bmatrix}\boldsymbol{\mathrm{I}}&\mathbf{0}_{3\times
3}\\\
(\mathbf{{}^{\mathcal{B}}_{\mathcal{W}}R}\,\mathbf{d}_{i})^{\times}(\mathbf{{}^{\mathcal{E}_{i}}_{\mathcal{W}}R}\,\mathbf{l}_{g})^{\times}&\boldsymbol{\mathrm{I}}\end{bmatrix}.$

where
$\mathbf{{}^{\mathcal{E}_{i}}_{\mathcal{W}}R}=\mathbf{{}^{\mathcal{U}_{i}}_{\mathcal{W}}R}\mathbf{{}^{\mathcal{E}_{i}}_{\mathcal{U}_{i}}R}$
map the transforming coordinates from frame $\mathcal{E}_{i}$ to frame
$\mathcal{W}$. $\mathbf{l}_{g}\in\mathbb{R}^{3}$ is the position of the origin
of frame $\mathcal{E}_{i}$ relative to frame $\mathcal{U}_{i}$. Combined with
the dynamics in (1), we obtain

$\displaystyle\boldsymbol{M}(\boldsymbol{q})\ddot{\boldsymbol{q}}+\boldsymbol{C}(\boldsymbol{q},\boldsymbol{\dot{q}})+\boldsymbol{g}=\sum_{i=1}^{N}\boldsymbol{G}(\boldsymbol{q},\boldsymbol{d_{i}})\mathbf{w}_{i}.$
(5)

Our objective is to design a decentralized control law for the wrench
$\mathbf{w}_{i}$ from each robot to the client that provides reference
tracking despite parametric uncertainty about the dynamics of the client,
$\boldsymbol{d}_{i}$ and $\boldsymbol{\bar{d}}$.

###  2.5 Decentralized Adaptive Transportation

In this section, we will design a decentralized adaptive controller for our
transportation task to let the group of vehicles manipulate the object to
track the desired position trajectory. First, we introduce a composite error
$\boldsymbol{s}$:

$\begin{aligned} \boldsymbol{s}=\begin{bmatrix}\boldsymbol{\epsilon}\\\
\boldsymbol{o}\end{bmatrix}=\begin{bmatrix}\dot{\boldsymbol{e}}_{p}+\beta\boldsymbol{e}_{p}(t)\\\
\boldsymbol{e}_{w}+\dfrac{1}{2}\beta\mathbf{{}^{\mathcal{B}}_{\mathcal{W}}R}_{d}(\mathbf{{}^{\mathcal{B}}_{\mathcal{W}}R}_{e}+\mathbf{{}^{\mathcal{B}}_{\mathcal{W}}R}_{e}^{\intercal})\end{bmatrix}\end{aligned},$
(6)

where $\boldsymbol{e}_{p}(t)=\boldsymbol{p_{d}}(t)-\boldsymbol{p}(t)$ is the
position tracking error, $\beta<1$ is a constant positive number.
$\boldsymbol{\dot{e}}_{p}(t)$ is the velocity tracking error,
$\boldsymbol{e}_{w}$ is the angular velocity error, and
$\mathbf{R}_{e}=\mathbf{R}_{d}^{\intercal}\mathbf{R}$ is the rotation error.
To solve Task-1, we use the “reference rotation error”
$\mathbb{Rot_{err}}=tr(\mathbf{R}_{e}-\boldsymbol{\mathrm{I}})$ to denote the
rotation error. we need $\mathbb{Rot_{err}}\rightarrow 0$ and
$\boldsymbol{e}_{p},\boldsymbol{e}_{w}\rightarrow 0$ as $t\rightarrow\infty$.

The decentralized adaptive controller for Task-1 uses the regressor matrix
$\boldsymbol{Y}_{g}(\boldsymbol{q},\boldsymbol{\dot{q}},\boldsymbol{\dot{q}_{d}},\boldsymbol{\ddot{q}_{d}})$,

$\displaystyle\sum_{i=1}^{N}\boldsymbol{Y}_{\varphi}\boldsymbol{\widetilde{\varphi}}_{i}=\sum_{i=1}^{N}\gamma(\boldsymbol{\widetilde{M}}\ddot{\boldsymbol{q}}_{d}+\boldsymbol{\widetilde{C}\dot{q}}_{d}+\boldsymbol{\widetilde{g}}),$
(7)

where $\gamma=\tfrac{1}{N}$ divides the control workload evenly among the
space tugs, and
$\boldsymbol{\widetilde{\varphi}}_{i}=\boldsymbol{\hat{\varphi}}_{i}-\boldsymbol{\varphi}_{i}$
denotes the error of parameters estimation related to the physical properties
of the service module, scaled by $\gamma$, and
$\boldsymbol{\widetilde{M}(u)}=\boldsymbol{\hat{M}(q)}-\boldsymbol{M(q)}$ is
the inertia error of the stiffness system,
$\boldsymbol{\widetilde{C}(\dot{q},q)}=\boldsymbol{\hat{C}(\dot{q},q)}-\boldsymbol{C(\dot{q},q)}$
is the Coriolis and centripetal error and
$\boldsymbol{\widetilde{g}}=\boldsymbol{\hat{g}}-\boldsymbol{g}$ is the
gravity error.

To provide asymptotic tracking despite parametric uncertainty, we use the
proportional-derivative control law with adaptive feedback linearization

$\displaystyle\mathbf{\hat{w}}_{i}=\boldsymbol{\hat{G}}_{i}\mathbf{w}_{i}=\boldsymbol{Y}_{\varphi}\boldsymbol{\hat{\varphi}}_{i}-\boldsymbol{K}_{\scriptstyle\text{PD}}\boldsymbol{s},$
(8)

where the term $\boldsymbol{Y}_{\varphi}\boldsymbol{\hat{\varphi}}_{i}$
feedback linearizes the service satellite dynamics (1) using the estimated
parameters $\boldsymbol{\hat{\varphi}}_{i}$ and the term
$-\boldsymbol{K}_{\scriptstyle\text{PD}}\boldsymbol{s}$ provides proportional-
derivative control. Note that the composite error (6) contains both position
and orientation errors as well as velocity and angular velocity errors,
producing a proportional-derivative controller.

Now let’s introduce another regressor
$\boldsymbol{Y}_{d}(\mathbf{{\hat{w}}}_{i},\boldsymbol{q})$ and then we can
get

$\displaystyle\boldsymbol{Y}_{d}(\mathbf{{\hat{w}}}_{i},\boldsymbol{q})\boldsymbol{\hat{d}}_{i}=-\boldsymbol{\widetilde{G}}_{i}\mathbf{{\hat{w}}}_{i}\;.$
(9)

Thus, we can get the adaptation laws as

$\displaystyle\dot{\hat{\boldsymbol{\varphi}}}_{i}$
$\displaystyle=-\Gamma_{\varphi}\boldsymbol{Y}_{\varphi}\left(\boldsymbol{q},\dot{\boldsymbol{q}},\dot{\boldsymbol{q}}_{d},\ddot{\boldsymbol{q}}_{d}\right)^{\intercal}\boldsymbol{s},$
(10a) $\displaystyle\dot{\hat{\boldsymbol{d}}}_{i}$
$\displaystyle=-\Gamma_{d}\boldsymbol{Y}_{d}\left(\boldsymbol{{\hat{W}}}_{i},\boldsymbol{q}\right)^{\intercal}\boldsymbol{s}.$
(10b)

We use a Lyapunov-like function to prove our controller design makes the whole
system exponentially stable. Consider the Lyapunov-like function

$\displaystyle
V(t)=\frac{1}{2}\left[\boldsymbol{s}^{\intercal}\boldsymbol{M}\boldsymbol{s}+\sum_{i=1}^{N}\widetilde{\boldsymbol{\varphi}}_{i}^{\intercal}\boldsymbol{\Gamma}\widetilde{\boldsymbol{\varphi}}_{i}+\widetilde{\boldsymbol{d}}_{i}^{\intercal}\boldsymbol{K}\widetilde{\boldsymbol{d}}_{i}\right],$
(11)

where $\widetilde{\boldsymbol{d}}_{i}$ is the position error that is the
difference between the actual position of $i$-th tug $\boldsymbol{d}_{i}$
concerning the measurement point $\boldsymbol{P}_{d}$ and the tug’s estimated
$\boldsymbol{\hat{d}}_{i}$. $\boldsymbol{\Gamma}$ and $\boldsymbol{K}$ are
symmetric positive definite matrices related the adaptive gain, usually
diagonal.

Then, by taking the derivative, we obtain

$\displaystyle\dot{V}(t)$
$\displaystyle=\boldsymbol{s}^{\intercal}(\sum_{i=1}^{N}\boldsymbol{G}(\boldsymbol{q,d_{i}})\mathbf{{w}}_{i}-\boldsymbol{M}\ddot{\boldsymbol{q}}_{d}-\boldsymbol{C\dot{q}}_{d}-\boldsymbol{g})$
(12)
$\displaystyle+\sum_{i=1}^{N}\left[\widetilde{\boldsymbol{\varphi}}_{i}^{\intercal}\boldsymbol{\Gamma}\dot{\hat{\boldsymbol{\varphi}}}+\widetilde{\boldsymbol{d}}_{i}^{\intercal}\boldsymbol{K}\dot{\hat{\boldsymbol{d}}}_{i}\right],$

where
$\dot{\hat{\boldsymbol{\varphi}}}_{i}=\dot{\widetilde{\boldsymbol{\varphi}}}_{i}$
and $\dot{\hat{\boldsymbol{d}}}_{i}=\dot{\widetilde{\boldsymbol{d}}}_{i}$
since the parameter $\boldsymbol{\varphi}_{i}$ and $\boldsymbol{d}_{i}$ are
constant. Then we use the properties of skew-symmetry [27] to eliminate the
term
$\boldsymbol{s}^{\intercal}(\dfrac{1}{2}\dot{\boldsymbol{M}}-\boldsymbol{C})\boldsymbol{s}=0$
and $\boldsymbol{\dot{q}}_{d}=\boldsymbol{s}-\boldsymbol{\dot{q}}$. Let’s
define the control-law as

$\displaystyle\mathbf{w}=\boldsymbol{\hat{M}}\ddot{\boldsymbol{q}}_{d}+\boldsymbol{\hat{C}\dot{q}}_{d}+\boldsymbol{\hat{g}}\;.$
(13)

We substitute (12) to (13) and then rewrite $\dot{V}(t)$ using the regressors
(7) and (9) and then substitute adaptation-laws (10) to obtain

$\displaystyle\begin{split}\dot{V}(t)&=\sum_{i=1}^{N}[-\boldsymbol{s}^{\intercal}\boldsymbol{K}_{{\scriptstyle\text{PD}}}\boldsymbol{s}+\boldsymbol{\widetilde{\varphi}}_{i}^{\intercal}(\boldsymbol{Y_{\varphi}\widetilde{\varphi}_{i}}\\\
&+\boldsymbol{\Gamma}_{\varphi}\dot{\hat{\boldsymbol{\varphi}}}_{i})+\boldsymbol{\widetilde{d}}^{\intercal}(\boldsymbol{Y_{d}\widetilde{d}_{i}}+\boldsymbol{\Gamma}_{d}\dot{\hat{\boldsymbol{d}}}_{i})]\\\
&=N\left(-\boldsymbol{s}^{\intercal}\boldsymbol{K}_{\scriptstyle\text{PD}}\boldsymbol{s}\right)\leq
0.\end{split}$ (14)

Thus, we get $V(t)>0$ and $\dot{V}(t)\leq 0$, the system is stable in the
sense of Lyapunov.

##  3 Simulation Results

This section presents two simulation results for transportation tasks using
two different aerospace robots: $N=4$ fully actuated hexarotor vehicles with
grasping systems in Section 3.1 and $N=4$ space tugs in Section 3.2. The
parameters of the object used in the two simulations are summarized in Table
1. Two simulations are implemented on MuJoCo [28] platform as shown in Fig.4.

Assume all aerospace robots in both simulations are distributed around the
object. Each aerospace robot uses its grasping scheme to firmly grasp the edge
to ensure the module can be transported in any direction and orientation in 3D
space. For the desired position trajectory
$\boldsymbol{p}_{d}=[p_{d_{x}},p_{d_{y}},p_{d_{z}}]^{\intercal}$, we take
$\boldsymbol{p}_{x}=\sin{(\omega_{x}t)}$,
$\boldsymbol{p}_{y}=\cos{(\omega_{y}t)}$, $\boldsymbol{p}_{z}=0$. Here we
adopt the iterative method [29] to simulate the trajectory of the robot which
we take the derivation form $\dot{\boldsymbol{p}}$ to
$\boldsymbol{p}(t)=\boldsymbol{p}(0)+\dot{\boldsymbol{p}}(t)\cdot t$ in which
$\boldsymbol{p}(0)$ is the initial state of the service module. For rotation
trajectory, here we assume $\boldsymbol{q}_{d}\in\mathbb{R}^{4}$ in both
simulations are $\boldsymbol{q}_{d}=[1,0,0,0]$, which means the object should
be kept in a fixed orientation during the transportation task.

Definition | Parameter | Value
---|---|---
Object | $m_{b}$ | 5 kg
| $I_{xx}$ |  $1.4255$ kg$\cdot$m2
| $I_{yy}$ |  $1.4255$ kg$\cdot$m2
| $I_{zz}$ |  $0.8411$ kg$\cdot$m2
Position | $\boldsymbol{l}_{g}$ |  $[0.1,0.0,-0.3]$ m
| $\mathbf{\bar{d}}$ |  $[0.74,0.0.1,-0.2]$ m
| $\mathbf{d}_{1}$ |  $[-0.8,1.2,0.1]$ m
| $\mathbf{d}_{2}$ |  $[1,1,0.1]$ m
| $\mathbf{d}_{3}$ |  $[1,-0.7,0.1]$ m
| $\mathbf{d}_{4}$ |  $[-0.7,-1.1,0.1]$ m
Table 1: Simulation parameters for transportation task.

###  3.1 Simulation-1: Transportation Task in Earth

((a))

((b))

Figure 4: Transportation task implemented on different aerospace robots on
MuJoCo platform. (a) shows the task by multiple fully actuated hexarotor UAVs
with a grasping system in the earth (with gravity); (b) shows the task by a
group of space tugs in space (without gravity).

Fig. 5(a) shows the tracking effect over time in 3D space. Fig. 5(b) shows the
composite error (6) for the client manipulated by a group of vehicles in the
space environment. The results show at $10$ seconds, we disable drone-1, which
introduces the disturbance to the system.

Fig. 5(c) and Fig. 5(d) show the parameter estimates for the physical
properties of the client and the contact point for each tug to the measurement
point $\boldsymbol{P}_{d}$. Based on the results, we can see that all signals
are bounded, and the final values are close to the actual values.

Fig. 5(e) shows the Lyapunov-like function $V(t)$ versus time for the
simulation results. The Lyapunov-like function is positive $V(t)>0$ and
decreasing $\dot{V}(t)\leq 0$ as time increases. This verifies that our system
is asymptotically stable.

((a))

((b))

((c))

((d))

((e))

Figure 5: Simulation results for using $N$ UAVs collaboratively to transport
the payload in $SE(3)$ under gravity environment.

###  3.2 Simulation-2: Transportation Task in Space

Fig. 6(a) shows the tracking effect over time in 3D space. Fig. 6(b) shows the
composite error (6) for the client manipulated by a group of space tugs in the
space environment. The results show at $10$ seconds, we disable tug-1, which
introduces the disturbance to the system.

Fig. 6(c) and Fig. 6(d) show the parameter estimates for the physical
properties of the client and the contact point for each tug to the measurement
point $\boldsymbol{P}_{d}$. Based on the results we can see that all signals
are bounded and the final of values are close to the actual values.

Fig. 6(e) shows the Lyapunov-like function $V(t)$ versus time for the
simulation results. The Lyapunov-like function is positive $V(t)>0$ and
decreasing $\dot{V}(t)\leq 0$ as time increases. This verifies that our system
is asymptotically stable.

((a))

((b))

((c))

((d))

((e))

Figure 6: Simulation results for using $N$ space tugs collaboratively to
transport the payload in $SE(3)$ under zero-gravity environment.

##  4 Conclusion and Future Work

We addressed the problem of centralized control in transportation tasks by
using multiple aerospace vehicles. All vehicles corporately hold and transport
a rigid object without prior knowledge, such as mass, inertia, and position
information relative to the estimated point and center of mass in frame
$\mathcal{B}$. We presented our decentralized adaptive controller design and
then implemented it on multiple fully actuated hexarotor UAVs in a gravity
environment first, and extended it to multiple space tugs for the same
transportation task, following the desired position and rotation in a zero-
gravity environment. We proved the feasibility of our controller in two
simulation environments. During the transportation mission, we deactivated one
of the vehicles to introduce disturbance to the whole system and examine the
robustness of our method. Two simulations using $N=4$ fully actuated hexarotor
UAVs in an Earth environment and $N=4$ space tugs in a space environment,
along with related results, demonstrated the efficacy of the presented
controller design.

We believe that the improvements made in designing the decentralized adaptive
control can be applied to other critical robotic servicing tasks where how
agents are assigned and distributed is crucial. These tasks include having
multiple robots working collaboratively in manufacturing processes. We are
currently working on learning-based robust adaptive control implemented on
multiple drones for transporting and dislodging non-rigid objects with unknown
physical properties, including mass, inertia, and material flexibility, to
handle servicing tasks using multiple drones in a wide range of applications.

##  5 Acknowledgement

This material is based on research sponsored by Air Force Research Laboratory
(AFRL) under agreements FA9453-18-2-0022 and FA9550-22-1-0093. We would like
to acknowledge Zhi Zheng from Zhejiang University for useful discussion in
dynamics analysis of the drone and also acknowledge Giovanni Cordova for his
assistance with Solidworks modeling.

## References

  * [1] Ian Down and Manoranjan Majji.  Adaptive detumbling of uncontrolled planar spacecraft using finite module deposition.  In AIAA SCITECH 2023 Forum, page 0158, 2023.
  * [2] Longsen Gao, Giovanni Cordova, Claus Danielson, and Rafael Fierro.  Autonomous multi-robot servicing for spacecraft operation extension.  In 2023 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pages 10729–10735. IEEE, 2023.
  * [3] Hamed Farivarnejad, Amir Salimi Lafmejani, and Spring Berman.  Fully decentralized controller for multi-robot collective transport in space applications.  In 2021 IEEE Aerospace Conference (50100), pages 1–9. IEEE, 2021.
  * [4] Isuru Basnayake, Hyeongjun Park, Jonathan Kohler, Jennifer Hudson, and Marcello Romano.  In-space demonstration of model predictive control approaches for space towing of uncertain loads.  In AIAA SCITECH 2024 Forum, page 1069, 2024.
  * [5] Daniel Hert, Tomas Baca, Pavel Petracek, Vit Kratky, Robert Penicka, Vojtech Spurny, Matej Petrlik, Matous Vrba, David Zaitlik, Pavel Stoudek, et al.  Mrs drone: A modular platform for real-world deployment of aerial multi-robot systems.  Journal of Intelligent & Robotic Systems, 108(4):64, 2023.
  * [6] Lucio R Salinas, Javier Gimenez, Daniel C Gandolfo, Claudio D Rosales, and Ricardo Carelli.  Unified motion control for multilift unmanned rotorcraft systems in forward flight.  IEEE Transactions on Control Systems Technology, 2023.
  * [7] Gustavo A Cardona, Diego S D’Antonio, Rafael Fierro, and David Saldaña.  Adaptive control for cooperative aerial transportation using catenary robots.  In 2021 Aerial Robotic Systems Physically Interacting with the Environment (AIRPHARO), pages 1–8. IEEE, 2021.
  * [8] Guanrui Li, Xinyang Liu, and Giuseppe Loianno.  Rotortm: A flexible simulator for aerial transportation and manipulation.  IEEE Transactions on Robotics, 2023.
  * [9] Sihao Sun and Antonio Franchi.  Nonlinear mpc for full-pose manipulation of a cable-suspended load using multiple uavs.  In 2023 International Conference on Unmanned Aircraft Systems (ICUAS), pages 969–975. IEEE, 2023.
  * [10] Antonio E Jiménez-Cano, Dario Sanalitro, Marco Tognon, Antonio Franchi, and Juan Cortés.  Precise cable-suspended pick-and-place with an aerial multi-robot system: A proof of concept for novel robotics-based construction techniques.  Journal of Intelligent & Robotic Systems, 105(3):68, 2022.
  * [11] Carlo Bosio, Jerry Tang, Ting-Hao Wang, and Mark W Mueller.  Automated layout design and control of robust cooperative grasped-load aerial transportation systems.  arXiv preprint arXiv:2310.07649, 2023.
  * [12] Bingguo Mu and Pakpong Chirarattananon.  Universal flying objects: Modular multirotor system for flight of rigid objects.  IEEE Transactions on Robotics, 36(2):458–471, 2019.
  * [13] Shraddha Barawkar, Manish Kumar, and Michael Bolender.  Decentralized adaptive controller for multi-drone cooperative transport with offset and moving center of gravity.  Aerospace Science and Technology, 145:108960, 2024.
  * [14] Dimitris Chaikalis, Nikolaos Evangeliou, Anthony Tzes, and Farshad Khorrami.  Modular multi-copter structure control for cooperative aerial cargo transportation.  Journal of Intelligent & Robotic Systems, 108(2):31, 2023.
  * [15] Daniel Mellinger, Michael Shomin, Nathan Michael, and Vijay Kumar.  Cooperative grasping and transport using multiple quadrotors.  In Distributed Autonomous Robotic Systems: The 10th International Symposium, pages 545–558. Springer, 2013.
  * [16] Bruno Gabrich, David Saldaña, Vijay Kumar, and Mark Yim.  A flying gripper based on cuboid modular robots.  In 2018 IEEE International Conference on Robotics and Automation (ICRA), pages 7024–7030, 2018.
  * [17] Zhongmou Li, Vincent Bégoc, Abdelhamid Chriette, and Isabelle Fantoni.  Wrench capability analysis and control allocation of a collaborative multi-drone grasping robot.  Journal of Mechanisms and Robotics, 15(2):021003, 2023.
  * [18] Alejandro Flores and Gerardo Flores.  Fully actuated hexa-rotor uav: Design, construction, and control. simulation and experimental validation.  In 2022 International Conference on Unmanned Aircraft Systems (ICUAS), pages 1497–1503. IEEE, 2022.
  * [19] S. Rajappa, M. Ryll, H. H. Bülthoff, and A. Franchi.  Modeling, control and design optimization for a fully-actuated hexarotor aerial vehicle with tilted propellers.  In 2015 IEEE International Conference on Robotics and Automation (ICRA), pages 4006–4013, 2015.
  * [20] Anibal Ollero, Marco Tognon, Alejandro Suarez, Dongjun Lee, and Antonio Franchi.  Past, present, and future of aerial robotic manipulators.  IEEE Transactions on Robotics, 38(1):626–645, 2022.
  * [21] Hossein Bonyan Khamseh, Farrokh Janabi-Sharifi, and Abdelkader Abdessameud.  Aerial manipulation—a literature survey.  Robotics and Autonomous Systems, 107:221–235, 2018.
  * [22] Chenyu Yang, Guo Ning Sue, Zhongyu Li, Lizhi Yang, Haotian Shen, Yufeng Chi, Akshara Rai, Jun Zeng, and Koushil Sreenath.  Collaborative navigation and manipulation of a cable-towed load by multiple quadrupedal robots.  IEEE Robotics and Automation Letters, 7(4):10041–10048, 2022.
  * [23] Preston Culbertson, Jean-Jacques Slotine, and Mac Schwager.  Decentralized adaptive control for collaborative manipulation of rigid bodies.  IEEE Transactions on Robotics, 37(6):1906–1920, 2021.
  * [24] Farhad Aghili.  Robust impedance-matching of manipulators interacting with uncertain environments: Application to task verification of the space station’s dexterous manipulator.  IEEE/ASME Transactions on Mechatronics, 24(4):1565–1576, 2019.
  * [25] Nikolay Atanasov.  Sensing & estimation in robotics.  Accessed: 4-2-2024.
  * [26] K Hassan Khalil.  Nonlinear systems.  Prentice Hall, 2002.
  * [27] Jean-Jacques E Slotine and Weiping Li.  On the adaptive control of robot manipulators.  The international journal of robotics research, 6(3):49–59, 1987.
  * [28] Emanuel Todorov, Tom Erez, and Yuval Tassa.  Mujoco: A physics engine for model-based control.  In 2012 IEEE/RSJ International Conference on Intelligent Robots and Systems, pages 5026–5033. IEEE, 2012.
  * [29] Josef Stoer and Roland Bulirsch.  Introduction to Numerical Analysis, volume 12.  Texts in Applied Mathematics, Springer, New York, 1996.