Caihong Xu

Min Wei

Jianmin Chen
email: jmchen@fudan.edu.cn

Chao Zhu

Jiarong Li

Xianmang Xu

Wenxing Wang

Qingzhu Zhang

Aijun Ding

Haidong Kan

Zhuohui Zhao

Abdelwahid Mellouki

Profile of Inhalable Bacteria in PM 2.5 at Mt. Tai, China: Abundance, Community, and Influence of Air Mass Trajectories

Keywords: Abundance, Community, and Influence of Air Mass Trajectories Bioaerosol, PM, Air mass origin, Diversity

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

The particulate matter with aerodynamic diameters of ≤2.5 μm (PM 2.5 ) is known as a major health risk for the spread of asthma, stroke burden, and chronic diseases such as nasopharyngitis, bronchitis, and pneumonia [START_REF] Bell | Seasonal and regional short-term effects of fine particles on hospital admissions in 202 US Counties, 1999-2005[END_REF][START_REF] Laden | Association of fine particulate matter from different sources with daily mortality in six U.S. cities[END_REF][START_REF] Zhai | A review on airborne microorganisms in particulate matters: composition, characteristics and influence factors[END_REF] . Assessment of air pollution exposure for the global burden of disease revealed that about 29.2% of stroke burden was associated with inhaled particulate matters at a high concentration [START_REF] Brauer | Ambient air pollution exposure estimation for the global burden of disease 2013[END_REF] . The environmental health scholars considered that the toxic effect of PM 2.5 on the respiratory tract and lungs mainly induced by the suspension and reaction of gases, heavy metals, mineral element, and toxic organic matters (cancerogenic substance) in ambient air [START_REF] Kleinman | Toxicity of chemical components of ambient fine particulate matter (PM 2.5 ) inhaled by aged rats[END_REF][START_REF] Maaike | In vitro toxicity of particulate matter (PM) collected at different sites in the Netherlands is associated with PM composition, size fraction and oxidative potential-the RAPTES project[END_REF] . In addition to these components, there is increasing evidence that pathogenic and allergenic microorganisms (bacteria, eukaryotes, archaea, and viruses) attached to the particles, may deposit in the epithelial cells through respiration, and thus causing various ailments in humans and livestock [START_REF] Brodie | Urban aerosols harbor diverse and dynamic bacterial populations[END_REF][START_REF] Xie | Seasonal disparities in airborne bacteria and sssociated antibiotic resistance Genes in PM2.5 between Urban and Rural Sites[END_REF][START_REF] Zhou | Occurrence and distribution of urban dust-associated bacterial antibiotic resistance in Northern China[END_REF][START_REF] Gao | Airborne microbial communities in the atmospheric environment of urban hospitals in China[END_REF] . The pathogenic species (e.g. adenovirus, influenza, SARS, and rhinovirus) and their components (e.g. endotoxins, glucans, and mycotoxins) widely distributed in near-surface atmosphere bring a great threat to the ecosystem balance [START_REF] Liu | Determination and risk assessment of airborne endotoxin concentrations in a university campus[END_REF][START_REF] Carty | Seasonal variability of endotoxin in ambient fine particulate matter[END_REF][START_REF] Chen | Ambient influenza and avian influenza virus during dust storm days and background days[END_REF] .

Bacteria, approximately 10 4 cells m -3 and 0.1 μg m -3 in PM, mainly originated from the emission of soil, vegetation, surface of rock, human and animals (excrements and debris from skin and hair), precipitation (surface snow and ice core), and ocean [START_REF] Despres | Primary biological aerosol particles in the atmosphere: a review[END_REF] . Investigations on the ambient bacteria experienced a recent growth in interest with the aerosol community. To our knowledge, the first study on airborne microorganisms was reported in the mid-19th century [START_REF] Ariatti | Louis Pasteur: the first experimental aerobiologist[END_REF] .

The occurrence of high throughput sequencing offers a unique possibility to examine the bacterial, viral, and eukaryotic diversity in the aerosol. While the understanding of airborne microbe is still being developed due to current challenges, e.g., low density of airborne microbes, variability in airborne microbial community structure, and the lack of unified regulations for the research. However, in spite of this, this method has been applied to study the airborne bacteria and the associated pathogens in Japan, USA, Korea, and China in the last few years [START_REF] Brodie | Urban aerosols harbor diverse and dynamic bacterial populations[END_REF][START_REF] Maki | Assessment of composition and origin of airborne bacteria in the free troposphere over Japan[END_REF][START_REF] Gou | Assessment of microbial communities in PM 1 and PM 10 of Urumqi during winter[END_REF][START_REF] Cha | Alterations in the airborne bacterial community during Asian dust events occurring between February and March 2015 in South Korea[END_REF] . Franzetti et al. (2011) [START_REF] Franzetti | Seasonal variability of bacteria in fine and coarse urban air particulate matter[END_REF] described the significant seasonal variation of bacterial community associated with PM 10 and PM 2.5 based on the Roche 454 sequencing.

Based on the high-throughput sequencing, Bower et al. (2013) [START_REF] Bowers | Seasonal variability in bacterial and fungal diversity of the near-surface atmosphere[END_REF] found the significant seasonal variability of Actinobacteridae and Pseudomonadales across PM 2.5 and PM 2.5-10 in Greeley, CO. Cao et al. (2014) [START_REF] Cao | Inhalable microorganisms in Beijing's PM 2.5 and PM 10 pollutants during a severe smog event[END_REF] developed a comprehensive framework of the microbial community in PM 2.5 and PM 10 under sever smog events by metagenomics. These findings mainly dressed the scanty composition community of bacteria in the urban air, the identification for the inhaled microbes and their pathogenic members at high altitudes were incomplete. Airborne bacteria with active metabolic have been reported in the mountains and the upper troposphere [START_REF] Maki | Assessment of composition and origin of airborne bacteria in the free troposphere over Japan[END_REF][START_REF] Maki | Variations in airborne bacterial communities at high altitudes over the Noto Peninsula (Japan) in response to Asian dust events[END_REF] . These species dissolving into cloud droplets proposed an important impact on the regional precipitation distribution, biogeochemical cycling of elements in ecosystem processes, and cloud chemistry by acting as cloud condensation nucleation or ice nucleation [START_REF] Bowers | Seasonal variability in airborne bacterial communities at a high-elevation site[END_REF][START_REF] Xu | Investigation of diverse bacteria in cloud water at Mt. Tai, China[END_REF] . Thus it is necessary to further the understanding for the profile of airborne bacterial community structure at high elevations.

Here, we expound on our initial efforts to advance the knowledge of airborne bacteria over Mt. Tai, the highest site at the North China Plain. It is located on the transport path of the Asian continent outflow and be affected by air mass from the Asian dust area. Mt. Tai is a tilted fault block mountain with height increasing from the north to the south, facing to the East China Sea and back on the yellow river to the west. In this region, the airborne bacterial populations in cloud system and fungal populations in particles have been reported before [START_REF] Xu | Investigation of diverse bacteria in cloud water at Mt. Tai, China[END_REF][START_REF] Xu | Fungi diversity in PM 2. 5 and PM 1 at the summit of Mt. Tai: abundance, size distribution, and seasonal variation[END_REF][START_REF] Lv | Investigation of new particle formation at the summit of Mt. Tai, China[END_REF] .

The preponderant gram-negative bacteria inhabited in the cold and leaf-surface environments which were involved in the cloud condensation/ice nuclei formation and pathogenicity, demonstrating their importance in ecology, epidemiology and human health [START_REF] Xu | Investigation of diverse bacteria in cloud water at Mt. Tai, China[END_REF][START_REF] Xu | Fungi diversity in PM 2. 5 and PM 1 at the summit of Mt. Tai: abundance, size distribution, and seasonal variation[END_REF] . The pioneering studies were lack of a comprehensive and compound research on the bacteria in particles. The chemical composition in particles have been proved to be related to the variable airborne microbe. Further, they also highlighted the crucial influence of environmental factors humidity, and temperature on the bacterial community in cloud water and fungal community in particles. Generally the chemical composition of particles may be influenced by different air mass origin with the input of local aerosol or long-range transported aerosol. Whether the different air trajectories do play a role in shaping bacterial community of particles still remains a question. In this study, we aim to: (1) profile the bacterial abundance and community of PM 2.5 at Mt. Tai, China, (2) to elucidate the environmental variables (temperature, humidity, PM 2.5 , PM 2.5 chemistry), and the air mass trajectories during sampling time, (3) to estimate whether the different air mass trajectories effect on the bacterial characteristics.

Materials and method

Aerosol collection and DNA Extraction

Fine particulate matter (PM 2.5 ) were sampled on quartz fiber filters (PALL, NY, U.S. 88mm) with a middle-volume sampler (TH-100A, Wuhan Tianhong Instruments Co., Ltd., China) and a flux speed of 100 L min -1 . The blank filters were obtained by placing a sterilized quartz microfiber filter inside of the sampler without operation. The sampling campaigns were carried out in 8~13 day periods in summer, autumn, and winter during July 2014 to August 2015. The sampler was located at the summit of Mt. Tai (36.25 o N and 117.10 o E, 1534 m above sea level), a typical tourist area in the North China Plain (Figure 1). The geographical characteristics of the sampling site have been previously described [START_REF] Xu | Fungi diversity in PM 2. 5 and PM 1 at the summit of Mt. Tai: abundance, size distribution, and seasonal variation[END_REF][START_REF] Lv | Investigation of new particle formation at the summit of Mt. Tai, China[END_REF] . Each collection period lasted 23 h. Before sampling, all the filters were baked at a Muffle furnace at 500°C for 5 h and then placed into sterilized aluminum foil before deposited into a sealed bag.

To avoid the contamination, the sampling filter holder and materials used for changing filters were treated with 75% ethanol before sampling every day. After sampling, the samples were stored at -80 o C until the next analysis. The environmental factors including temperature, relative humidity, visibility, water-soluble inorganic ions (NO 3 -, SO 4 2-, NH 4 + , Ca 2+ , Na + , Mg 2+ , and Cl -), meteorological data, the mass concentration of PM 2.5 , DNA extraction procedure and PCR amplification have been described before [START_REF] Lv | Investigation of new particle formation at the summit of Mt. Tai, China[END_REF] . The final product was ready for resolving for the abundance and community of airborne bacteria by real-time quantitative PCR and Illumina Miseq sequencing.

Quantification of 16S rRNA gene copy number

The bacterial abundance was determined by the quantification of the copy numbers of the gene 16S rRNA. This value is not directly related to bacterial number due to the presence of multiple ribosomal operons in the bacterial genomes; however, it can be used to look at relative shifts in microbial biomass as the bias in copy number is likely constant across samples. The V4-V5 hypervariable regions were PCR-amplified with universal primers 515F

(5'-TGCCAGCMGCCGCGGTAA-3') and 926R (5'-CGTCAATTCMTTTGAGTTT-3'). The real-time quantitive PCR was performed in a total volume of 25 μL reaction mixture including 12.5 Μl TransStart Green qPCR SuperMix (Transgene, AQ101-03), 1μL forward primer, 1μL reverse primer, 5 μL DNA, and 5.5 μL double distilled H 2 O. The amplification was run on a BioRad CFX96 Real-Time PCR system under the following conditions: 95 °C for 5 min, followed by 40 cycles of 95 °C for 30 s, 56 °C for 30 s and 72 °C for 30 s, with acquisition of the fluorescence at the end of each 72 °C elongation step. The deionized water was employed as the negative controls and was pipetted into the wells in a 96-well microplate with the DNA extracted from the samples. All the samples and controls were tested in triplicate in each run. Based on the average rRNA gene copy number (3.98), the bacterial cell concentration (cells m -3 ) was calculated using the method described before [START_REF] Xu | Bacterial characterization in ambient submicron particles during severe haze episodes at Ji'nan, China[END_REF][START_REF] An | Development and calibration of real-time PCR for quantification of airborne microorganisms in air samples[END_REF][START_REF] Seunghoon | Identification of airborne bacterial and fungal community structures in an urban area by T-RFLP analysis and quantitative real-time PCR[END_REF] .

Sequence Analyses

A DNA sequencing library was constructed on the basis of the under accession numbers SRR51461. The potential pathogenic sequences corresponding to human health were screened referenced to the previous database [START_REF] Zhou | Continuous observations of water-soluble ions in PM 2.5 at Mount Tai (1534 m a.s.l.) in central-eastern China[END_REF][START_REF] Guo | Profiling bulking and foaming bacteria in activated sludge by high throughput sequencing[END_REF][START_REF] Deng | Chemical characterization of aerosols at the summit of Mountain Tai in Central East China[END_REF][START_REF] Ye | Pathogenic bacteria in sewage treatment plants as revealed by 454 pyrosequencing[END_REF] . The human pathogenic bacteria list, including species names and the diseases caused, was used to find the potential pathogens sequences from total sequences. We selected the operational taxonomic units (OTUs) belonged to these genera mentioned in the reference list. The nearest phylogenetic species were obtained by NCBI BLAST. Finally, the bacterial species belonging to the reference list were determined as the potential pathogens in PM 2.5 .

Statistical Analyses

The difference in bacterial concentration across different air mass trajectories were assessed by an ANOVA test, followed by post-hoc tests. The ANOVAs test was also conducted to determine whether season had a significant impact on the relative abundances of specific microbial groups. The weighted UniFrac analysis was employed to estimate the significant differences of the microbial community among three seasons. The relationship between the airborne microbial concentrations and environmental factors including PM concentrations and chemical compositions was assessed with a Pearson correlation by SPSS 16.0. The results that Bonferroni adjusted P-value for each pairwise comparison higher than 0.05 were considered to be significant. The correlation coefficient between NMDS plot and environmental elements was calculated via metadata analysis with Mothur software.

Statistical values, including Pearson and Spearman correlation coefficients and independent ttests, were calculated using R software with the R package.

Results and discussion

3.1 Bacterial abundance in PM 2.5

During the sampling period, the PM 2.5 mass concentration, ranging from 3 to 111 μg m -3 with a mean value of 40 μg m -3 at Mt. Tai, were much lower than that in summer 2006

(123 μg m -3 ) [START_REF] Deng | Chemical characterization of aerosols at the summit of Mountain Tai in Central East China[END_REF] and comparable with that in summer 2007 (59 μg m -3 ) 30 that conducted at the same site. Table 1 summarizes the corresponding bacterial concentrations in a different size.

The bacterial concentration of PM 2. In general, the airborne bacteria, ranging between 10 and 10 7 cells m -3 , varied in a degree with many aspects [START_REF] Liu | Determination and risk assessment of airborne endotoxin concentrations in a university campus[END_REF][START_REF] Bowers | Seasonal variability in airborne bacterial communities at a high-elevation site[END_REF][START_REF] Seunghoon | Identification of airborne bacterial and fungal community structures in an urban area by T-RFLP analysis and quantitative real-time PCR[END_REF][START_REF] Bertolini | Temporal variability and effect of environmental variables on airborne bacterial communities in an urban area of Northern Italy[END_REF] , like sampling time, site (altitude, temperature, humidity, etc.), sample type (TSP, PM 10 , PM 2.5 , PM 1 , cloud/fog water, etc.), and analysis methods (culture methods, bio-molecular technology, and quick estimation model) [START_REF] Despres | Primary biological aerosol particles in the atmosphere: a review[END_REF][START_REF] John | Microorganisms of the upper atmosphere III. Relationship between altitude and micropopulation[END_REF][START_REF] Li | Comparison of the biological content of air samples collected at ground level and at higher elevation[END_REF] . And these indexes effectively represented the situation of bacterial community on particles under different conditions. [START_REF] Borodulin | Physical characteristics of concentration fields of tropospheric bioaerosols in the South of Western Siberia[END_REF] observed the varied concentration of culturable microorganisms for about more than two orders of magnitude over an altitude range of 0.5~7 km in southwestern Siberia. They found the overall bacterial concentration in the atmosphere where not only dependent on altitude but also related to the cell size of bacteria and particulate matters.

Borodulin et al. (2005)

Different bacteria have respective sizes ranging from 0.1~5 μm, and each bacterium has a particular shape feature that plays an important role in the adhesion on PM [START_REF] Zhai | A review on airborne microorganisms in particulate matters: composition, characteristics and influence factors[END_REF] . Since the study environments were breeding settings in which bacteria had better conditions to grow and propagate, the bacteria concentration may decrease with decreasing particulate size [START_REF] Zheng | Concentrations and size distributions of airborne particulate matter and bacteria in an experimental aviary laying-hen chamber[END_REF] . Hence more high-throughput sequencing researches for the bacterial community distribution were vital to understanding their diversity fully. [START_REF] Despres | Characterization of primary biogenic aerosol particles in urban, rural, and high-alpine air by DNA sequence and restriction fragment analysis of ribosomal RNA genes[END_REF] , the airborne bacteria at Mt. Tai commonly exhibited lower richness and diversity compared with those from urban, marine, sand, and dust sites. The diversity estimators can reflect the state of the bacteria and their relationship to the environment, and microbial resources developments, ecosystem functioning, global change and human health are all closely related to microbial diversity. In addition to the changing sampling sites, this discrepancy may be caused by the different analysis platform (Miseq sequencing, pyrosequencing, and clone library). Using the same method (Miseq Sequencing), TSP seems to have higher OTUs, chao1, and Shannon indexes than PM 2.5 and PM 1 , indicating more and various bacterial species in the larger particulate matter. We attribute this increasement to the larger surface area of TSP for bacterial adhesion. As mentioned above, the size of a single bacterial cell is smaller than particles, there is a mass of nutrients and bacterial slime coat on the bacteria that make the bacterial size bigger but still smaller than particles, so there were more bacteria presenting in particles with a larger size. Actually different bacteria had respective variation tendency for particle sizes, although there were some differences on size distribution, bacteria mainly attached to particles with two peaks (1.1~2.1 μm and >7.1 μm) [START_REF] Raisi | Size distribution of viable, cultivable, airborne microbes and their relationship to particulate matter concentrations and meteorological conditions in a Mediterranean site[END_REF] .

Except for the bacterial size, the structure of bacterial community can also exert a certain extent influence on the diversity characteristics.

At Mt. Tai, the diverse bacterial community including 40 phyla, 78 classes, 154 orders, 286 families, 724 genera was classified in PM 2.5 . About 97.82% of OTUs belongs to the phylum Proteobacteria, Cyanobacteria, Firmicutes, and Bacteroidetes. The co-occurrence network were used to reveal the complex bacterial community structure based on the top 100 OTUs (Figure 2). The species of dominant phyla (Proteobacteria, Actinobacteria, Cyanobacteria, Firmicutes, Acidovacteria, Fusobacteria, and Bacterioidetes) were similar in USA [START_REF] Brodie | Urban aerosols harbor diverse and dynamic bacterial populations[END_REF][START_REF] Bowers | Spatial variability in airborne bacterial communities across land-use types and their relationship to the bacterial communities of potential source environments[END_REF] , Italy [START_REF] Bertolini | Temporal variability and effect of environmental variables on airborne bacterial communities in an urban area of Northern Italy[END_REF] , and East Asian Regions [START_REF] Ariatti | Louis Pasteur: the first experimental aerobiologist[END_REF][START_REF] Cao | Inhalable microorganisms in Beijing's PM 2.5 and PM 10 pollutants during a severe smog event[END_REF][START_REF] Maki | Vertical distribution of airborne bacterial communities in an Asian-dust downwind area, Noto Peninsula[END_REF][START_REF] Katra | Richness and diversity in dust stormborne biomes at the southeast mediterranean[END_REF][START_REF] Seifried | Spatial distribution of marine airborne bacterial communities[END_REF][START_REF] Bowers | Spatial variability in airborne bacterial communities across land-use types and their relationship to the bacterial communities of potential source environments[END_REF] , but the relative abundance were different. Furthermore, the microbiome is a dynamic aspect of the near-surface layer among extreme weather incidents such as dust, haze, rain, and tropical storms. Proteobacteria members were the dominant species in the atmosphere over mountains and from the troposphere. The relative abundance of this bacterium (91.42%) were richer than that in the dust (69.0%), nondust (53.5%) [START_REF] Cha | Alterations in the airborne bacterial community during Asian dust events occurring between February and March 2015 in South Korea[END_REF] , and haze (16.5%) [START_REF] Xu | Bacterial characterization in ambient submicron particles during severe haze episodes at Ji'nan, China[END_REF] days. It was dominated by a few taxa throughout the whole sampling period, with a marked prevalence of Alphaproteobacteria, Betaproteobacteria, and Gammaproteobacteria. At higher taxonomic levels, a bacterium in the nitrogen fixing, gram-negative pathogenic bacteria genus Burkholderia appeared to be the most abundant. [START_REF] Fang | Culturable airborne bacteria in outdoor environments in Beijing, China[END_REF] who reported the 80%~86% culturable bacteria in outdoors environment were gram-positive. They considered the gram-negative bacteria cannot adapt to the higher air pollutants concentration, stronger solar radiation, and complicated meteorological environment outdoors. We speculate this gap attributed to the different analysis platform. The Illumina Miseq sequencing method gained the comprehensive airborne bacteria no matter culturable or unculturable. Actually, the proportion of culturable bacteria is very small.

Another reason is that those gram-negative bacteria in the atmosphere lose cultivability after aerosolization [START_REF] Amato | Microorganisms isolated from the water phase of tropospheric clouds at the Puy de Dôme: major groups and growth abilities at low temperatures[END_REF] . Hence the majority of airborne bacteria outdoors were gram-negative.

Out of all the gram-negative species, eight opportunistic pathogens were determined based on the human pathogenic bacteria list as described before. Acinetobacter baumannii, primarily associated with hospital-acquired infections, commonly found in water and soil.

Over the past decades, reports on community-acquired A. baumannii infections (wound infection, ventilator-associated pneumonia, and bloodstream infection) have increased and attracted public attention [START_REF] Peleg | Acinetobacter baumannii: emergence of a successful pathogen[END_REF][START_REF] Xie | Analysis of global prevalence of antibiotic resistance in Acinetobacter baumannii infections disclosed a faster increase in OECD countries[END_REF][START_REF] Zhang | Treatment options for multidrug-resistant and extensively drugresistant Acinetobacter baumannii infections[END_REF][START_REF] Hawkey | Evolution of carbapenem resistance in Acinetobacter baumannii during a prolonged infection[END_REF] . An outbreak of soil-borne bacteria Pseudomonas putida was involved in clinic infections in patients with medical devices and immunosuppressed hosts 64- 65 . Besides, five gram-positive potential pathogens were detected including Bacillus subtilis, Rhodococcus equi, Streptococcus anginosus, Bacillus pumilus, and Rhodococcus fascians.

Some strains of B. subtilis cause food poisoning through the ropy spoilage of bread and rots in potatoes [START_REF] Valerio | Use of Lactobacillus plantarum fermentation products in bread-making to prevent Bacillus subtilis ropy spoilage[END_REF][START_REF] Apetroaie-Constantin | Bacillus subtilis and B. mojavensis strains connected to food poisoning produce the heat stable toxin amylosin[END_REF] . R. equi, mainly thrive on water and soils, were connected with pneumonia in immunocompromised patients and animal hosts [START_REF] Prescott | Rhodococcus equi: an animal and human pathogen[END_REF][START_REF] Yamshchikov | Rhodococcus equi infection[END_REF] . S. anginosus were able to cause diseases including bloodstream infections, pyogenic liver abscess, brain and liver abscesses [START_REF] Laupland | Streptococcus anginosus group bloodstream infections in the western interior of British Columbia, Canada[END_REF][START_REF] Murarka | Pyogenic liver abscess secondary to disseminated streptococcus anginosus from sigmoid diverticulitis[END_REF] . These bacterial species suspending in the atmosphere diffuse downwind, yielding a threat to the health of plants and animals in local and downstream area.

Other than the pathogenicity, the diverse bacteria demonstrates a diverse metabolism in the atmospheric process and ecosystem interactions. The bacteria with ice nucleation activity suspending in the atmosphere can take part in the cloud chemistry via act as cloud condensation or ice nuclei. During the same project, we captured some typical cloud events (Event F and G) [START_REF] Xu | Investigation of diverse bacteria in cloud water at Mt. Tai, China[END_REF] , incidentally. The special terrain of Mt. Tai is conducive to the formation of orographic clouds. The tiny droplets could dissolve or adsorb the particles suspending in the atmosphere and then accumulated into the cloud water. Though the same sampling site, the dominant bacterial species varied from Burkholderia in PM 2.5 to the preponderance of Pseudomonas in cloud water. The distant discrepancy of dominating bacterial distribution also found between PM 2.5 of Mt. Tai and PM 1 of Ji'nan (Figure 3A). As to this phenomenon, different bacterial suitable inhabits and growing media (water, soil, atmosphere, etc.) may be important factors. Additionally, Liu et al. ( 2018) [START_REF] Liu | Effect of air pollution on the total bacteria and pathogenic bacteria in different sizes of particulate matter[END_REF] considered the relative humidity, wind scale, CO, and SO 2 were the main factors, which affected the diversity of total bacteria and the proportion of pathogenic bacteria. The PM's chemical components serve as the nutritional ingredient, coupled with surrounding environments such as temperature, humidity, wind direction, and atmospheric pressure seems to shift on the bacterial growth and metabolism in atmosphere [START_REF] Zhai | A review on airborne microorganisms in particulate matters: composition, characteristics and influence factors[END_REF][START_REF] Gao | Airborne microbial communities in the atmospheric environment of urban hospitals in China[END_REF][START_REF] Xie | Characteristics of total airborne microbes at various air quality levels[END_REF] . The chemical and biological data were jointly processed and verified that the varied NH 4 + , NO 3 -, Mg 2+ , Ca 2+ , Cl -and Na + explained about 53.7%~63.1% of the variance of species-environment correlations [START_REF] Innocente | Influence of seasonality, air mass origin and particulate matter chemical composition on airborne bacterial community structure in the Po Valley, Italy[END_REF] . The chemical composition of PM , NH 4 + , and NO3 -) and reactions between particulate species. We also detected the obvious correlation between airborne bacterial community and Ca 2+ , Mg 2+ , and temperature at Mt. Tai, while the influence of atmospheric circulation using the backtrajectory analysis was not evaluated fully. Therefore how does the air mass trajectories vary and whether the air mass trajectories evoked the change of concentration and structure of bacterial community is worth discussing.

Cluster analysis of backward trajectory and PM 2.5 chemistry

At Mt. Tai, the air mass back trajectories have been categorized into four types, typified by the examples in Figure 4, based on trajectories terminating at the site on the day of measurements. One trajectory per day was, in our experience, adequate to represent the flow patterns which tended to change over timescales of a few days. Inevitably, there were a few days on which assignment of trajectory type was difficult, but the predominant type from the four synoptic hours was adopted in these circumstances. The four types of back trajectory were as follows:

Cluster I (25%): The air mass had passed across the marginal seas in the western Pacific Ocean and had been in contact with the land immediately prior to arrival in the study area.

Broadly speaking, the highest mass concentration of PM 2.5 (61 μg m -3 ) and major water inorganic ions were recorded followed the decreasing order: SO 4 2 -> NO 3 -> Cl -> Na + > Ca 2 + > K + > NH 4 + > Mg 2 + . In addition to the secondary inorganic ions, Cl -(4.35 μg m -3 ) and Na + (3.73 μg m -3 ) can be assumed mostly to originate from sea salt aerosols. The Cl -/Na + molar ratios obtained (0.75) was lower than that of bulk seawater (1.18) [START_REF] Yan | Hygroscopic behavior of water-soluble matter in marine aerosols over the East China Sea[END_REF] , suggesting that the chlorine loss of marine aerosols during transport to some extent. These ions are proved to play essential roles in atmospheric processes. e.g., cloud formation, visibility degradation, solar radiation, acidification of cloud, rain, fog and haze formation and served as the nutrients for the growth and development of airborne bacteria [START_REF] Xu | Investigation of diverse bacteria in cloud water at Mt. Tai, China[END_REF][START_REF] Wei | Characteristics of bacterial community in cloud water at Mt Tai: similarity and disparity under polluted and non-polluted cloud episodes[END_REF] .

Cluster II (33%): The air mass had come in from a generally western direction bringing it across Shandong peninsula and its adjacent areas, a typical densely polluted areas of northern China. Influenced by the local terrestrial air flow, the average of PM 2.5 mass was 37 μg m -3 .

The secondary inorganic ions, i.e., SO 4 2-(9.92 μg m -3) , NO 3 -(4.23 μg m -3 ), and NH 4 + (2.29 μg m -3 ), which mostly originated from the transformation of their precursors (SO 2 , NO 2 , and NH 3 )

were found as the dominant species, followed by Cl -(0.75 μg m -3 ), K + (0.63 μg m -3 ), Ca 2+ (0.51 μg m -3 ), Na + (0.31 μg m -3 ), and Mg 2+ (0.07 μg m -3 ). In order to investigate the continental influence, non-sea-salt sulfate (nss-SO 4 2 -) concentrations are estimated as:

. The mass fraction of nss-SO 4 2ranged from 1.17 μg m -3 to 18.89 μg m -3 (mean 9.85 μg m -3 ), accounted for 99.02% of total sulfate. ] mass ratios at Mt. Tai (0.25) was higher or comparable with those reported in other cities, e.g., 0.17 in Hongkong 77 , 0.20 in Guangzhou [START_REF] Han | Chemical characteristics of PM 10 during the summer in the mega-city Guangzhou, China[END_REF] and 0.35 in Qingdao [START_REF] Hu | Seasonal variation of ionic species in fine particles at Qingdao, China[END_REF] . However, the ratio was lower than those cities where the use of coal was limited, for example, Los Angeles (2.0) [START_REF] Kim | Characterization of PM 2.5 and PM 10 in the South Coast Air basin of Southern California: Part 1-Spatial variations[END_REF] and Seoul (0.73) [START_REF] Kim | Characteristics of the major chemical constituents of PM 2.5 and smog events in Seoul, Korea in 2003 and 2004[END_REF] . The result implied that stationary source emissions were likely more important than mobile sources at Mt. Tai under the influence of local terrestrial air flow.

Cluster III (25%) and Cluster IV (17%): The air mass arrived from a northwesterly sector outer Mongolia and having crossed dust area (Inner Mongolia and Shanxi Province) and high anthropogenic emission area (Henan and Shandong Province). This type air mass moves quickly and the corresponding mass concentration of coarse particulate matter was usually high. Typically, these westerly and northwesterly trajectories had not been in contact with any land for several days prior to arriving at the site.

Influence of air mass trajectories on bacterial populations and community

Obviously, no obvious variation for the bacterial concentration was observed under the influence of air mass trajectories. The terrestrial air masses from the Shandong peninsula and its adjacent areas exhibited a little higher bacterial concentration loading (644 cells m -3 ) and more potential pathogens. While the bacterial concentration of PM influenced by the marine air mass was lower, though the highest chemical composition (Figure 4). In contrast, samples influenced by the long-distance air flow from Siberia and Outer Mongolia were found to have a highest richness and diversity as an average (Table 2). The effects of long-transported air parcels on the aerosol particles and its composition and sources have been discussed in previous studies [START_REF] Castañer | Influence of air mass origins on optical properties and PM concentrations measured at a high mountain station located in the southwestern Mediterranean[END_REF][START_REF] Masiol | A chemometric approach to determine local and regional sources of PM 10 and its geochemical composition in a coastal area[END_REF][START_REF] Moreda-Piñeiro | The influence of oceanic air masses on concentration of major ions and trace metals in PM 2.5 fraction at a coastal European suburban site[END_REF][START_REF] Squizzato | Application of meteorology-based methods to determine local andexternal contributions to particulate matter pollution: A case study in Venice (Italy)[END_REF] . The Asian dust particles blow from these areas in the westerly wind and play an important role in the leeward ecosystem. The bacterial community at phyla, class, and order levels under the influence of four air mass trajectories was shown in Figure 5. However, the bacterial abundance made no significant difference among the top 20 OTUs (Figure 3B), suggesting the stable community structure of airborne bacteria. Yet, some non-dominant bacterial groups that exhibited pronounced changes in relative abundance across the different air mass trajectories (Figure 6) can be linked to specific source environments of PM 2.5 . The long-range transported air parcels seem to be able to harbor diverse and dynamic microbial populations [START_REF] Deleon-Rodriguez | Microbiome of the upper troposphere: species composition and prevalence, effects of tropical storms, and atmospheric implications[END_REF] . Kakikawa et al., (2008) (**p<0.01). Consistent with these findings, the long-range transported westerly air mass increase the diversity and richness of airborne bacterial community, with the increase of some dustborne species such as genus Brevibacillus, Staphylococcus, and Anoxybacillus (Figure 6A). The Bacilli members obtained from mineral dust particles over downwind areas were positively correlated to temperature and Cl -in PM 2.5 (**p<0.01), excepting the obvious effect of humidity and temperature. They showed a stronger ability to resist ultraviolet light and be able to survive in the atmosphere by the formation of resistant endospores. The natural dust particles from dust areas can be transported to the free troposphere and vertically mixed with the anthropogenic particles during the transportation processes. Together, the airborne bacteria community at Mt. Tai contains a mixture of dustborne bacteria, terrestrial bacteria, marine bacteria, and others. The community of dominant bacteria in PM 2.5 was relatively stable, and little affected by the sampling environment. In other words, ambient bacterial community structure were stable in atmosphere. Four air mass trajectories from various origin environments may lead to the discrepant composition in PM 2.5 , and then driving the shift on the pronounced changes observed in some non-dominant bacterial groups.

Conclusions

In this study, 16S rRNA gene sequencing was conducted to study the bacterial communities at Mt. Tai's PM bacterial concentration loading and more potential pathogens at the site. In contrast, samples influenced by the long-distance air flow from Siberia and Outer Mongolia were found to have a highest richness and diversity as an average, which was also marked by the increased abundance of dust-associated bacteria (Brevibacillus and Staphylococcus). We conclude that the air flow from different origin environments may lead to the discrepant chemical composition in PM 2.5 , and then driving the shift on the pronounced changes observed in some non-dominant bacterial groups. We hope our research can provide for some reference value in the global air microbial propagation and air microbial research at higher elevation site.

List of Figures

  5 varied from 40 to 1728 cells m -3 with an average of 490 cells m -3 at Mt. Tai. It was much lower than those TSP from the high-altitude air, urban Japan, Korea, and Italy, PM 2.5 and PM 1 from China. The lower bacterial concentration were found to be correlated with PM 2.5 mass concentration (r = 0.526, p < 0.05), SO 4 2-(r = 0.549, p < 0.01), NO 3 -(r = 0.458, p < 0.05), NH 4 + (r = 0.570, p < 0.01), K + (r = 0.508, p < 0.01), Na + (r = 0.585, p < 0.01), Mg 2+ (r = 0.608, p < 0.01), and Ca 2+ (r = 0.432, p < 0.05) significantly based on the Spearman rank correlations. The ambient pollutants attached to the particles may provide nutrients (sulfur, nitrogen, and ammonia) and thus affect the bacterial characteristics.

  Bowers et al. (2009) 57 also showed the dominant Burkholderiales in particles collected at the top of Mt. Werner near Steamboat Springs CO, USA at 3220 m above sea level via tagged pyrosequencing. The preponderant gram-negative bacteria in PM 2.5 were different with Fang et al. (2007)

2 . 5 correlated

 25 with the air mass back trajectories. Air mass originating from Mongolia region had strong correlations with dust pollution, along with high concentration of Ca 2+ in ambient PM 2.5 . Local emission and air pollutants were more favorable for formation of secondary pollutants (SO 4 2-

SO 4 2 -

 2 mainly originates from anthropogenic and biogenic sources, while nss-SO 4 2originates mainly from stationary sources. We often use the mass ratio of [NO 3 -]/[nss-SO 4 2 -] as an indicator of the relative importance of stationary versus mobile sources of sulfur and nitrogen in the atmosphere. High [NO 3 -]/[nss-SO 4 2 -] mass ratios represent the predominance of mobile sources over stationary sources. Due to the extensive use of sulfur-containing coal, the average ratio [NO 3 -]/[nss-SO 4 2 -] mass ratios in China were ascribed lower value 76 . The measured [NO 3 -]/[nss-SO 4 2 -

  verified the injection of dustborne bacteria (Brevibacillus, Staphylococcus, Rhodococcus, Delftia, Pseudomonas, and Agrobacterium tumefaciens) into the atmosphere during the dust event[START_REF] Kakikawa | 5.1% Bacteroidetes 1.2% Actinobacteria 1.1% Acidobacteria 1.5% Gammaproteobacteria 1.3% Betaproteobacteria 46.4% Alphaproteobacteria 41.7% Clostridia 0.5% Bacilli 0.7% Melainabacteria 5.1% Actinobacteria 1.0% Acidobacteria 1.4% CLUSTER IV PROTEOBACTERIA Proteobacteria 87.6% Firmicutes 5.1% Bacteroidetes 1.5% Actinobacteria 1.9% Acidobacteria 0.9% Cyanobacteria 2.4% Acidobacteria 0.9% Gammaproteobacteria 3.1% Betaproteobacteria 46.0% Alphaproteobacteria 38.2% Clostridia 0.8% Bacilli 4.2% Melainabacteria 2.3% Sphingobacteriia 0.6% Bacteroidia 0.5% Actinobacteria 1[END_REF] . Maki et al., (2017) also observed the variable microbial community in response to the Asian dust events[START_REF] Maki | Variations in airborne bacterial communities at high altitudes over the Noto Peninsula (Japan) in response to Asian dust events[END_REF] . The number of microbial particles at high altitudes (>1200m) increased about 5.5-folds. The corresponding bacterial community was predominantly composed of natural-sand/terrestrial bacteria during dust events instead of anthropogenic/terrestrial bacteria, marine bacteria, and plant-associated bacteria during non-dust events. The prevalent existence of marine bacteria in cloud droplets, the air masses suspended from the Yellow Sea over the marine environment harbored more proteobacterial populations, Pandoraea, and Phyllobacterium associated with plant bodies or surface. The prevalence of Pandoraea had a significant positive relationship with the highest content of SO 4 2-

Figure 1

 1 Figure 1The geographic situation of Mt. Tai in the North China Plain.

Figure 2

 2 Figure 2The co-occurrence network of top 100 OTUs of PM 2.5 at Mt. Tai, China. The bigger the circle, the greater the number of OTUs.

Figure 3

 3 Figure 3The principal component analysis showed the overall distribution pattern of A bacterial communities in different air samples and B profiles among four air mass cluster.

Figure 4 5 Figure 5

 455 Figure 4 Cluster analysis of air-mass trajectories and corresponding chemical composition and bacterial abundance in PM 2.5

Figur 6

 6 Figur 6 Variance analysis of genera species based on the Kruskal-Wallis rank sum test (A) and the Spearman's rank correlations with environmental factors (B; *p <0.05; * * p <0.01).

  Figure 1 708

  Figure 2

  /www.drive5.com/uclust/downloads1_2_22q.html) at different similarity levels to determine the diverse bacteria at phylum, class, order, family, genus, and species levels. The sequences obtained in this study have been deposited in the Sequence Read Archive (SRA)

	the	Silva	database	(http://www.arb-silva.de/)	using	uCLUST
	(http:/
					TruSeq DNA PCR-Free
	Sample Prep Kit (Illumina) and was sequenced by the MiSeq reagent kit v2 (Illumina) on an
	Illumina MiSeq instrument (Illumina, San Diego, CA, USA). The raw sequences were

analyzed using standard Mothur software. Variable sequence numbers (7044 to 32,293) were obtained after the quality control of raw reads by removing the chimeric sequences, primers, barcodes, and low-quality reads (average quality score<20) using FASTX-ToolKit. We normalized the sequences number to 7044 for the subsequent analysis. The Operational Taxonomic Units (OTUs) were defined using a distance-based cutoff method by Usearch (version 7.1; http://drive5.com/uparse/). The alpha-diversity estimators Chao1 and Shannon indexes, which indicated the richness and diversity of the bacterial community, were calculated by QIIME (version 1.17; http://qiime.org/scripts/assign_taxonomy.html) after removing the low-abundance OTUs (<0.01%). The taxonomic classification was aligned to

  2.5 pollutants from July 2014 to August 2015. Results demonstrated the lower bacterial concentration and diverse bacterial community in PM 2.5 compared with some aforementioned studies. We also detected thirteen potential pathogenic species, such as A. baumannii, P. putida, B.subtilis, and R. equi which effect on the immunocompromised populations and livestock. In future studies, we intend to analyze the changing bacterial community structure affected by different air mass trajectories, combined with the meteorological data, PM mass concentration, water soluble ion and other physical and chemical properties. The highest mass concentration of PM 2.5 and major chemical components was recorded during periods when marine southeasterly air masses were dominant. The local

terrestrial air masses from Shandong peninsula and its adjacent areas harbored highest

Table 1 Bacterial abundance in PM with different size fraction.

 1 Copies m -3 ; ** CFU m -3 ; *** Cloud Water

	Type	Site	Abundance Range	Average	Method	Reference
	TSP	Cells m -3 High-altitude air, Caribbean Sea 3.0×10 3 ~9.0×10 4	Cells m -3 5.1×10 3	qPCR	39
	TSP	High-altitude air, Japan	1.0× 10 3 ~2× 10 5		qPCR	40
	TSP	Toyama, Japan	1.1× 10 3 ~1.3×10 5	1.8×10 4	qPCR	41
	TSP	Seoul, Korea	5.2×10~4.3×10 3		qPCR	29
	PM 10	Milan and Venice, Italy	2.4×10 2 ~4.9×10 5*	3.0×10 4*	qPCR	42
	PM 2.5	BTH, China	8.0×10~5.8×10 3**	1.1×10 3** Culture method 43
	PM 2.5	BTH, China	4.8×10 4 ~ 2.6×10 5	1.2×10 5	qPCR	44
	PM 1	Ji'nan, China	1.8×10 4 ~2.6×10 4	2.2×10 4	qPCR	27
	CW *** Mt.Tai, China	5.8×10 4 ~1.6×10 5	1.2×10 5	qPCR	45
	PM 2.5	Mt.Tai, China	4.0×10~1.728×10 3 4.9×10 2	qPCR	This study

*

Table 2 Statistical summary of alpha diversity estimators from different sample type 692 and study environments. 693

 2

	Environment	Sample Type OTUs Chao1 Shannon	Method	Reference
	Urban/Terrestrial PM 2.5	2457	4915	5.89	Miseq Sequencing	46
	Urban/Terrestrial PM 2.5	552	571	3.80	Miseq Sequencing	44
	Urban/Terrestrial PM 0.56-1	222	272	2.44	Miseq Sequencing	27
	Dust/Terrestrial	PM 2.5	62	200	3.51	Clone Library	47
	Dust/Terrestrial	TSP		270		Pyrosequencing	48
	Dust/Terrestrial	Dust	869	1485		Pyrosequencing	49
	Dust/Terrestrial	PM 10	733	1265		Pyrosequencing	18
	Sand/Terrestrial	Sand	1445	1858	5.99	Pyrosequencing	50
	Marine/Coastal	TSP	753			Pyrosequencing	51
	Marine/Coastal	TSP	5015		7.66~9.28 Miseq Sequencing	52
	High altitude	TSP	410			Pyrosequencing	22
	High altitude	TSP	314	241		Miseq Sequencing g 53
	High altitude	CW	*	664	1302	3.80	Miseq Sequencing	24
	High altitude	PM 2.5	580	1282	3.01
	Cluster I	PM 2.5	491	1151	2.67
	Cluster II	PM 2.5	567	1322	2.93	Miseq Sequencing	This Study
	Cluster III	PM 2.5	677	1397	3.09
	Cluster IV	PM 2.5	625	1234	3.63

*

Cloud

Water

Acknowledgements

This work was supported by the Ministry of Science and Technology of China (No. 2016YFC0202700), the National Natural Science Foundation of China (No. 91743202, 21527814). This project has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement No 690958.

Conflict of interest

The authors declare that they have no competing financial interests.


Dario Sanalitro
email: dario.sanalitro@laas.fr

Heitor J Savino
email: heitor.judiss-savino@laas.fr

Marco Tognon
email: marco.tognon@laas.fr

Juan Cortés
email: juan.cortes@laas.fr.

Antonio Franchi
email: a.franchi@utwente.nl

Full-pose Manipulation Control of a Cable-suspended load with Multiple UAVs under Uncertainties

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

I. INTRODUCTION

Unmanned Aerial Vehicles (UAVs) are nowadays used in different domains, from search & rescue and inspection to precise agriculture and transportation of loads [START_REF] Palunko | Agile load transportation: Safe and efficient load manipulation with aerial robots[END_REF]. In the latter case, the object is either rigidly attached to the vehicle [START_REF] Spurnỳ | Cooperative autonomous search, grasping, and delivering in a treasure hunt scenario by a team of unmanned aerial vehicles[END_REF] or attached using a cable [START_REF] Goodarzi | Geometric control of a quadrotor uav transporting a payload connected via flexible cable[END_REF]. Multiple aerial vehicles rigidly attached to the load [START_REF] Ritz | Carrying a flexible payload with multiple flying vehicles[END_REF] or attached to the load's center of mass (CoM) using cables [START_REF] Bernard | Autonomous transportation and deployment with aerial robots for search and rescue missions[END_REF] are typically used to overcome limited payload of the single-vehicle solution.

In the last few years, motivated by manipulation tasks such as assembly, construction, contact-based inspection, etc., the focus moved toward the more complex problem of controlling the load full pose, i.e., both position and orientation. In most of the manipulation tasks, the full pose control is a strong requirement, e.g., transporting objects that cannot enter in limited spaces unless tilting them, or for pegin-hole operations where just controlling the position of the load is not enough.

Among several solutions, the one based on cables is the most suitable for the mentioned applications thanks to the induced decoupling between the rotational dynamics of the vehicles and the one of the load. Full pose control of the load is possible when the robotic system can attain a sixdimensional wrench on the load by controlling the orientation of the cables and the force applied by them (called stress). The minimum setup that possesses such ability is a system with three vehicles connected to three non-collinear points on the load through three cables [START_REF] Prattichizzo | Grasping[END_REF] (see Fig. 1 (a) for a schematic example, and [START_REF] Sreenath | Dynamics, control and planning for cooperative manipulation of payloads suspended by cables from multiple quadrotor robots[END_REF] for additional details). The main drawback of this setup is that the compensation of a time-varying external disturbing wrench (e.g., wind) must be obtained by moving the three aerial vehicles to reorient the three cables. This makes the overall system sluggish and imprecise in real-world conditions. To overcome such drawback one has to employ a statically rigid [START_REF] Connelly | Frameworks, tensegrities and symmetry: Understanding stable structures[END_REF] (or forceclosure [START_REF] Prattichizzo | Grasping[END_REF]) design. Such design allows compensating timevarying external wrenches almost instantaneously by varying only the stresses while keeping all the cables in the same orientation. Six is the minimum number of cables in order to gain such a property. A possible setup is the one using six cable-robot pairs, as in Fig. 1 (b).

However, the use of six aerial vehicles considerably increases the overall system complexity and costs compared to the three aerial vehicle case. A simpler and more efficient choice consists in using only three aerial vehicles connected to a pair of cables each, see Fig. 1 (c). This statically rigid system is called the Fly-Crane and was first presented in [START_REF] Manubens | Motion planning for 6-D manipulation with aerial towed-cable systems[END_REF].

Several algorithms have been proposed to control the pose of a cable-suspended load handled by multiple UAVs. In [START_REF] Geng | Implementation and demonstration of coordinated transport of a slung load by a team of rotorcraft[END_REF], [START_REF] Tan | Cooperative control of multiple unmanned aerial systems for heavy duty carrying[END_REF] and [START_REF] Dhiman | Cooperative load control and transportation[END_REF], position and yaw angle of a slung payload are controlled in real scenarios. In all these works, as in [START_REF] Marina | Flexible collaborative transportation by a team of rotorcraft[END_REF], regulation is achieved using a communicationbased formation control of the UAVs. A communicationless approach, which relies on the sensing of contact forces, is instead proposed in [START_REF] Tognon | Aerial comanipulation with cables: The role of internal force for equilibria, stability, and passivity[END_REF]. Differently, in [START_REF] Sreenath | Dynamics, control and planning for cooperative manipulation of payloads suspended by cables from multiple quadrotor robots[END_REF] the UAVs follow in open loop the nominal trajectories computed from the desired load trajectory using differential flatness. In order to cope with the high sensitivity to model errors and external disturbances of open-loop approaches, other works employ geometric control techniques [START_REF] Lee | Geometric control of quadrotor uavs transporting a cablesuspended rigid body[END_REF] or feedback control laws that require the inversion of the Full Dynamic Model (FDM) [START_REF] Masone | Cooperative transportation of a payload using quadrotors: A reconfigurable cable-driven parallel robot[END_REF]. Although the feedback on the load's pose makes the overall approach more robust to non-ideal scenarios, these methods still require a good knowledge of all the dynamic model parameters. However, it is known that if feedback controlled systems do not take model uncertainties into account, stability may be lost in real world conditions. Few works started to address this problem, e.g., [START_REF] Geng | Implementation and demonstration of coordinated transport of a slung load by a team of rotorcraft[END_REF] showed the robustness of the system to external disturbances, [START_REF] Lee | Geometric control of quadrotor uavs transporting a cablesuspended rigid body[END_REF] suppresses undesired behaviors when facing random uncertainties, [START_REF] Dhiman | Cooperative load control and transportation[END_REF] and [START_REF] Ariyibi | Control of a quadrotor formation carrying a slung load using flexible bars[END_REF] consider the payload itself as a disturbance to be controlled.

Departing from the current state of the art, the main contribution of this work is to propose a novel algorithm aiming at the full-pose control of the load for the Fly- Crane cable-suspended aerial manipulator. The proposed method does not rely on the dynamical model and formally takes into account model parameter uncertainty and external disturbances in its design phase. The advantage of reasoning at the kinematic level rather than at the dynamic one is that we need to identify a smaller set of parameters. Furthermore, the parameters are related only to geometric quantities, which are easier to measure with high accuracy. Furthermore, in our approach, the dynamics of the system is not completely neglected, but all its effects are rather considered as disturbances, and formally addressed in the design phase of the controller. As a drawback, to limit disturbances, and so tracking errors, the acceleration of the desired platform motion should be kept low. On the contrary, an FDM controller theoretically allows acrobatic motions as well. Having quasi-static trajectories is the price to pay to obtain a simpler and more precise identification process of the system parameters. This is not a real problem for the considered manipulation-related applications (e.g., construction and assembly), where accuracy is preferred over agility.

Relying on an inner/outer loop approach, we designed the proposed controller such that it is simple to implement in real world. In the inner loop, the three UAVs aim at regulating their velocities using simple decentralized proportional controllers. The outer loop provides velocity references to the inner loop in order to regulate the system configuration, including the full pose of the load. We demonstrate how the control gain of the outer loop can be designed to: i) cope with the parametric uncertainty of the system model and ensure stability despite the non-perfect knowledge of its kinematic parameters; ii) minimize the effect of external disturbances by optimizing the H ∞ performance index. The proposed framework is tested by real indoor experiments under different conditions.

The paper is organized as follows. Section II describes the Fly-Crane system and its modeling. The proposed robust controller is presented in Sec. III. The experimental results are shown in Sec. IV and on the attached multimedia material. Final discussions and future works are presented in Sec. V.

II. THE FLY-CRANE SYSTEM

For the reasons mentioned above, we believe that the Fly-Crane is one of the most suitable systems for aerial manipulation. The Fly-Crane system is composed by three aerial vehicles and a platform, which is in turn constituted by a load and six cables arranged as in Fig. 2. In order to define the Fly-Crane model upon which the robust controller will be designed, let us define an inertial frame F W = {O W , x W , y W , z W }, where O W is its origin and {x W , y W , z W } are its unit axes. Then we define another frame F L = {O L , x L , y L , z L } rigidly attached to the load. F L is placed such that O L is centered on the CoM of the load. The configuration of the load is then given by the vector Connections between cable-load and cable-robot are done such that no rotational constraints are present. We assume that all the points B i lie on the plane {x L , y L }. As normally done in the state of the art, we assume in the model that the cables have negligible mass and inertia w.r.t. the other bodies of the system. The model also neglects phenomena like sagging and elongation and assumes the cables always taut2 . Thus the generic i-th cable, with i = 1, . . . , 6, has a constant length l i ∈ R >0 in the model. These assumptions are also motivated by the fact that for a Fly-Crane system, the cables are never very long with respect to the dimensions of the system. In fact, differently from cable-driven robots, the Fly-Crane possesses a large positional workspace regardless of the length of the cables. Furthermore, long cables are undesirable because they make the two-cable vehicle-load connection similar to a singlecable connection (see Fig. 1 (a)), at the detriment of the sought static rigidity property.

The configuration of the pair of cables (i, j) ∈ {(1, 2), [START_REF] Goodarzi | Geometric control of a quadrotor uav transporting a payload connected via flexible cable[END_REF][START_REF] Ritz | Carrying a flexible payload with multiple flying vehicles[END_REF], [START_REF] Bernard | Autonomous transportation and deployment with aerial robots for search and rescue missions[END_REF][START_REF] Prattichizzo | Grasping[END_REF]} with respect to F L is represented by the angle α ij ∈ R between the plane {x L , y L } and the plane formed by the cables (i, j). The platform configuration is then entirely described by the variable q

= [p L η L α ] ∈ C ⊂ R 9 , where α = [α 12 α 34 α 56 ] T ∈ R 3 .
The three positions of the robots, p Ri ∈ R 3 , with i = 1, 2, 3, depend on q through the kinematic relation:

p Ri (q) = p L + R L (η L ) P p Ri (α ij ), (1)
where R L (η L ) ∈ SO(3) is the rotation matrix computed from η L ,

F L F W B 1 B 2 B 3 B 4 B 5 B 6 O R1 = O R2 O R3 = O R4 O R5 = O R6 α 12
P p Ri = P b i + l i RP bij (α ij )R z L (β i ) P bij P bij , P b ij is the position of the vector ---→ B i B j expressed in F L ,

and

β i ∈ R is the angle between ---→ B i B j and ----→ B i O Ri .
From (1) we can compactly define the kinematic model of the system as:

p R = f (q), (2)
where

p R = [p R1 p R3 p R5 ] . Differentiating (2)
, we obtain the differential kinematic model: 3) relates velocities of the platform and cable angles rates to the velocities of the robots. The Jacobian matrix J (q) ∈ R 9×9 is a square invertible matrix in C, except for some singular configurations investigated in [START_REF] Six | The kinematics, dynamics and control of a flying parallel robot with three quadrotors[END_REF].

v R := ṗR = ∂f (q) ∂q q = J (q) q, (3) where
, v R = [v R1 v R3 v R5 ] , q = [v L η L α ] ∈ R 9 . Equation (
For configurations in which J (q) is non-singular, the coordinate p R and its time-derivative v R represent an alternative way to describe the platform configuration q and its velocity q, respectively, via the inversion of ( 2) and [START_REF] Goodarzi | Geometric control of a quadrotor uav transporting a payload connected via flexible cable[END_REF].

In this work, we focus on quasi-static operations for which high-order dynamic effects can be considered negligible. Under this assumption, the kinematic model given by ( 2) and ( 3) is enough to describe the system. Nevertheless, possible errors due to this assumption shall be considered in the control design and its stability proof.

III. ROBUST CONTROL OF THE FLY-CRANE

We assume that a motion planner, like the one presented in [START_REF] Manubens | Motion planning for 6-D manipulation with aerial towed-cable systems[END_REF], provides at each time t a desired configuration q d (t) that is far from singularities and ensures taut cables.

The proposed control scheme of the Fly-Crane is shown in Fig. 3. From left to right, the desired trajectory quantities q d , qd and the measured configuration q are used by the outer configuration controller to generate the reference aerial-vehicle velocities v R . Using v R and the measured velocities v R , the velocity controllers of the aerial vehicles compute the commanded accelerations for each vehicle u = [u 1 u 3 u 5 ] . Such commanded accelerations are provided to the low level actuation units of the aerial vehicles, which are in control of the orientation of the vehicles and of the

q d , qd v R u v R v R q q
Fig. 3: The control architecture of the Fly-Crane where the configuration controller of the platform and the velocity controllers of the three aerial vehicles are highlighted.

rotor speeds. In the following we detail the equations used in each controller. To face possible errors of the low-level actuation units, an uncertain model for the aerial vehicles is considered.

We denote with Ĵ (q) ∈ R 9×9 the nominal Jacobian of the system, such that

Ĵ (q) = J (q) + ∆J (q), (4)
where ∆J (q) is the error between the nominal (identified) and real (unkown) Jacobian matrices. Given the desired configuration

q d = [p d L η d L α d ] , the corresponding generalized velocities qd = [v d L ηd L αd ]
, and the measured configuration q, the reference aerial-vehicle velocities are computed:

v R = Ĵ (q) K q e q + qd , (5)
where 3 K q = k q I 9 ∈ R 9×9 >0 is a positive definite matrix and e q = q dq. The commanded accelerations for the vehicles are computed simply as

u = K R (v R -v R ) , (6)
where

K R = k R I 9 ∈ R 9×9
>0 is a positive definite matrix. Each aerial vehicle is equipped with a low-level actuation unit that acts on the orientation and intensity of the total force produced by its propellers in order to let vRi match the given acceleration input u i , for i = 1, 2, 3. However, the non perfect knowledge of the platform dynamics and the presence of external disturbances make such matching inexact. The closed-loop dynamics of the vehicle including the actuation unit is then equivalent to an uncertain double integrator, i.e.,

v Ri = ṗRi , vRi = u i + d i , (7)
where the uncontrollable signal d i accounts for all the aforementioned dynamic uncertainties and external disturbances. Such signal d i can be large, small or even close to zero, depending on the type of platform (e.g., under-or fully-actuated), type of controller (e.g., full dynamic model inversion, vehicle model inversion plus disturbance rejection, etc.) and the desired trajectory. Defining the vector d =

Summarizing, the closed loop system dynamics is given by the equations [START_REF] Connelly | Frameworks, tensegrities and symmetry: Understanding stable structures[END_REF], where u is computed using (32) and ( 5). In such model, the terms ∆J (q) and d account for all the non-idealities arising from modeling errors at the dynamic and kinematic level and the neglected external disturbances. To analyze the system stability under non-ideal conditions, let us write the closed-loop dynamics of the error. Replacing ( 5) and (32) into [START_REF] Connelly | Frameworks, tensegrities and symmetry: Understanding stable structures[END_REF] and applying the differential kinematics (3) yield

J (q) q + J (q, q) q = k R Ĵ (q)k q e q + k R Ĵ (q) qd -k R J (q) q + d. (9)
Considering the previous equation with nominal Jacobian as in (4), we can write the error dynamics as

ëq = -k R k q + k R k q J -1 (q)∆J (q) e q -k R + J -1 (q) J (q, q) ėq + w, (10)
where w contains the terms that do not depend on the error and that are considered as disturbances:

w = qd + J -1 (q) J (q, q) -k R ∆J (q) qd -d . (11)
Normally, the control gain K R is tuned independently to (5) to achieve the best tracking of the velocity reference v R . In this case, we show how to design the configuration control gain K q = k q I 9 , such that to guarantee stability of the closed-loop system with desired H ∞ performance, under some sufficient conditions on the uncertainty and disturbance magnitudes.

In the following, for an arbitrary matrix M , λ min [M ] and λ max [M ] are its minimal and maximal eigenvalues, M = λ max [M M ] is its induced norm, and M > 0 (M < 0) stands for a positive (negative) definite matrix. 0 n ∈ R n×n is the zero matrix of dimension n × n. Finally, with the symbol * we refer to the symmetric term in a matrix.

Theorem 1. Let us consider the control laws (32) and (5) applied to the Fly-Crane under disturbance d and uncertain Jacobian matrix Ĵ (q). k R is given while k q has to be designed. Assume that the following bounds hold true for given positive constants δ, ι, ρ ∈ R >0 :

∆J (q) ∆J (q) ≤ δ 2 , (12)
λ min [J (q) J (q)] ≥ ι 2 , (13) J (q) J (q) ≤ ρ 2 .

(

) 14
The designed control gain k q stabilizes the system with H ∞ performance defined by

e q ≤ γ w , (15)
with a given constant γ > 0, if there exist scalars q and q k and matrices

P = P , F 1 = F 1 , F 2 = F 2 , G 1 = G 1 ,
and G 2 = G 2 of dimension 9 × 9 that satisfy the LMI conditions P qI 9 * qI 9 > 0

and

      Ξ E 2 q E 2 k R q k E 1 q E 1 k R q k * -2F 1 0 9 0 9 0 9 * * -2G 1 0 9 0 9 * * * -2F 2 0 9 * * * * -2G 2       < 0, (17)
with

E 1 = [I 9 0 9 0 9 ] , E 2 = [0 9 I 9 0 9 ] , (18)
and Ξ given as

  Ξ 1,1 1 2 P -1 2 k R q k I 9 -1 2 k R qI 9 qI 9 * ρ 2 2ι 2 (F 1 + F 2 ) -k R qI 9 + qI 9 qI 9 * * -γ 2 I 9   (19) with Ξ 1,1 = I 9 + δ 2 2ι 2 (G 1 + G 2 ) -k R qI 9 .
In particular, the designed controller gain k q is then given by k q = q k /q.

Proof. Provided in the Appendix.

From a design perspective, given the bounds δ, ι, ρ and the velocity controller gain k R , we would like to find the K q that provides the best H ∞ performance. We can then rewrite Theorem 1 as a constrained minimization problem of γ: minimize β γ subject to (16), ( 17),

where β = (P , F 1 , F 2 , G 1 , G 2 , q, q k , δ, ι, ρ, k R , γ). From the solution of (20), and in particular from q k and q, the optimal gain k q is computed as k q = q k /q. The main advantage on this method is to find a feasible controller, with the best H ∞ gain on view of Theorem 1, that formally guarantees stability within bounds of uncertainties in the Jacobian matrix (δ), the proximity to singular configurations (ι), and its speed (ρ).

Remark. In order to solve (20) a first estimation of the bounds δ, ι, and ρ is needed. This can be first computed along the desired trajectory knowing the maximum uncertainty on the Jacobian. The latter comes from the imprecision in the design and manufacturing of the platform. Therefore:

δ 2 > max t∈R ≥0 ∆J (q d (t)) ∆J (q d (t)) , ι 2 < min t∈R ≥0 λ min [J (q d (t)) J (q d (t))], ρ 2 > max t∈R ≥0 J (q d (t)) J (q d (t)) . (21)
Since the bounds strongly depends on the particular trajectory, a motion planner such as the one in [START_REF] Manubens | Motion planning for 6-D manipulation with aerial towed-cable systems[END_REF], could be employed in order to compute the desired trajectory that minimizes the bounds. Notice that also the magnitude of w in (11) depends on the particular motion and especially on the system velocity and acceleration. In fact, w can be bounded as

w ≤ qd + ι (ρ + k R δ) qd + d . (22)
Therefore, the minimization of w can also be considered at the planning level. Equation ( 22) shows that performance is related to desired velocities and accelerations. Since we are more focused on the accuracy rather than on the agility of the maneuvers, we consider that in the domain of interest the motion of the system is done in a quasi-static condition.

Remark. We remark that our controller does not forbid the use a FDM-based method as a low-level controller for the vehicles. This would reduce the velocity tracking error d i , improving the overall performance. On the other hand, this would require a precise estimation of the dynamic parameters, which is hard to get as mentioned before.

IV. EXPERIMENTAL RESULTS

To validate and show the performances of the proposed controller with respect to uncertainties and external disturbances, we designed, implemented and experimentally tested the Fly-Crane in an indoor environment across different conditions.

The Fly-Crane consists of three Quadrotor UAVs connected by six cables to a platform made of carbon fiber bars. The Fly-Crane realization is shown in Fig. 4. Each vehicle weights 1.03 [Kg], the platform weights 0.338 [Kg], and the length of each cable is 1.2 [m]. The aerial vehicles are equipped with a standard flight-controller, four brushless motor controllers regulating the propeller speed in closedloop [START_REF] Franchi | Adaptive closed-loop speed control of BLDC motors with applications to multi-rotor aerial vehicles[END_REF], and an on-board PC that runs the state estimator and the velocity controller 4 . The proposed kinematic controller, implemented in Matlab-Simulink, runs on a desktop PC sending the commanded velocities to the quadrotors through a wifi connection at 100 [Hz]. The control loop is then closed based on the estimated state of the vehicles and of the platform. These estimations are computed onboard at 1 [kHz] by an UKF that fuses the Motion Capture (Mo-Cap) System measurements (at 120 [Hz]) with the IMU measurements (at 1 [kHz]).

In the very first phase of each experiment, the platform is lifted from ground, and the system is brought to a non-singular initial configuration. In this phase, the aerial vehicles are independently controlled by a standard position controller. The proposed controller is activated right after.

The trajectory is designed to take the system from its initial configuration, q d (0), to the final desired configuration q d (T ) where Given the desired trajectory and a maximum model uncertainty of 10% for each parameter, (21) allows computing 4 The full software framework is based on TeleKyb which is open-source and available at https://git.openrobots.org/ projects/telekyb3 the bounds needed for the design of k q . In particular, the following bounds were considered, δ 2 = 0.055, ι 2 = 0.078, and ρ 2 = 0.01. The gain k q can be then designed using Theorem 1. The solution to the minimization problem (20) provides k q = 0.38 and H ∞ gain γ = 1.

Considering this set-up, we performed three experiments under different conditions:

1) Optimal conditions: the model parameters are the one measured with no additional errors (measurement errors and noise are still present). External disturbances were not applied neither. Due to limited space, we show the relative plots in the attached multimedia material only. For experiment 2), in Fig. 5 one can notice that conditions (12), ( 13) and ( 14) are always verified. These conditions are fulfilled for experiments 1) and 3) as well. In Fig. 6 and Fig. 8 the norm and the components of the tracking error are shown for experiments 2) and 3), respectively. We can observe that the tracking error does not go exactly to zero, as expected, due to parameter uncertainties, noisy measurements and external disturbances. These non-idealities affects the exactness of the nominal model and the velocity tracking of the aerial vehicles. However, considering the first row of Fig. 6 and Fig. 8 where e q and γ w are compared, it is clear that condition [START_REF] Lee | Geometric control of quadrotor uavs transporting a cablesuspended rigid body[END_REF] is verified and the maximum H ∞ performance is guaranteed. We can conclude that, under the fulfillment of conditions ( 12), [START_REF] Marina | Flexible collaborative transportation by a team of rotorcraft[END_REF], and ( 14) in Theorem 1, the proposed controller and tuning method guarantee the maximum H ∞ performance. For experiment 2) only, Fig. 7 shows the desired trajectory with respect to the one performed by the platform. Similar plots corresponding to experiments 1) and 3) are reported in the attached multimedia material.

Although out of the scope of this work, we tested our system for far from quasi-static trajectories. In this 'exper-IEEE Robotics and Automation Letters 2020 12)-( 14) for experiment 2). iment 4)', we asked the system to accelerate the platform up to 3.2 [m/s 2 ] reaching a velocity up to 1.6 [m/s]. As expected, the tracking error increases but the bounds and the H ∞ gain of Theorem 1 are still verified. The corresponding plots are reported in the attached multimedia material.

V. CONCLUSIONS

A controller capable of dealing with uncertainties in the parameters of the Fly-Crane system has been provided given bounds that can be assessed using the planned trajectory. These bounds can be guaranteed in a real scenario when the Fly-Crane is asked to perform specific trajectories translating and rotating the platform in different directions. The structure of the controller gain K q was limited to a diagonal matrix with equal values to allow design conditions. In the future, we intend to relax this structure and to extend the controller to a second order kinematic approach. The inclusion of a robust dynamic controller is left as future works. Consider the Lyapunov candidate function V (e q , ėq ) =

V 1 + V 2 + V 3 , with: V 1 =
1 2 e q P e q , V 2 = 1 2 ė q q ėq , V 3 = e q q ėq .

To have V as a Lyapunov function we need V > 0, ∀ e q = 0, ėq = 0. Re-writing V , this is verified if

V = 1 2
e q ėq P qI 9 * qI 9 e q ėq > 0,

which is satisfied only if LMI ( 16) holds. Next, we consider the time-derivative of V . For V 1 , V1 = e q P ėq . Concerning V 2 , after replacing [START_REF] Geng | Implementation and demonstration of coordinated transport of a slung load by a team of rotorcraft[END_REF] we obtain VQ ≤ -ė q qk R k q J -1 (q)∆J (q)e q -ė q qk R k q e q ė q qk R ėqė q qJ -1 (q) J (q) ėq + ė q qw. (25)

Given the relation -2a b ≤ a X -1 a + b Xb, X > 0, choosing X as F 1 > 0, we can use ( 13) and ( 14) to bound:

ė q qJ -1 (q) J (q) ėq ≤ 1 2 ė q qF -1 1 q ėq + 1 2 ė q (J -1 (q) J (q)) F 1 J -1 (q) J (q) ėq

≤ 1 2 ė q qF -1 1 q ėq + ρ 2 2ι 2 ė q F 1 ėq . (26)
In the same way, choose X as G 1 , with [START_REF] Dhiman | Cooperative load control and transportation[END_REF], to bound:

-ė q qk R k q J -1 (q)∆J (q)e q ≤ 1 2 ė q (qk R k q )G -1 1 (qk R k q ) ėq + δ 2
2ι 2 e q G 1 e q . Finally, we have that the derivative of V 2 is bounded by

VQ ≤ 1 2 ė q (qk R k q )G -1 1 (qk R k q ) ėq + δ 2
2ι 2 e q G 1 e q ė q qk R k q e q -ė q qk R ėq

+ 1 2 ė q qF -1 1 q ėq + ρ 2 2ι 2 ė q F 1 ėq + ė q qw. (27)
Taking the time-derivative of term V 3 and proceeding analogously, for some F 2 > 0 and G 2 > 0, we can obtain

VR ≤ 1 2 e q (qk R k q )G -1 2 (qk R k q ) e q + δ 2
2ι 2 e q G 2 e q e q qk R k q e q -e q qk R ėq + ė q q ėq + 1 2 e q qF -1 2 qe q + ρ 2 2ι 2 ė q F 2 ėq + e q qw.

(28)

To show that stability with the H ∞ performance is satisfied we use the following condition:

e q e q -γ 2 w w

+ V < 0, (29)
given the Lyapunov candidate function V and a scalar γ > 0. Integrating both sides in any time-interval [0, T ] and assuming zero initial conditions e q (0) = 0 and ėq (0) = 0, we have T 0 e q e q -γ 2 w w dt + V (e q (T ), ėq (T )) < 0, (30) such that, if condition (29) is satisfied then (15) holds and the system is stable with H ∞ performance.

To verify condition (29), since V = V1 + V2 + V3 , we use the upper bounds in ( 27) and ( 28) into (29) to write ξ Aξ < 0 (31) with ξ = [e q ė q w ],

A = Ξ + E 2 q 1 2 F -1 1 (E 2 q) + E 2 qk R k q 1 2 G -1 1 (E 2 qk R k q ) +E 1 q 1 2 F -1 2 (E 1 q) + E 1 qk R k q 1 2 G -1 2 (E 1 qk R k q ) ,
Ξ given in [START_REF] Franchi | Adaptive closed-loop speed control of BLDC motors with applications to multi-rotor aerial vehicles[END_REF], and E 1 and E 2 in (18). Thus, for the condition (31) to hold, A has to be imposed negative definite.

Applying the Schur's Complement on this term, and writing q k = qk q to allow for control design, yields the LMI condition in [START_REF] Ariyibi | Control of a quadrotor formation carrying a slung load using flexible bars[END_REF]. If the LMI conditions are satisfied, k q can be obtained by taking k q = q k /q. This completes the proof.

B. Additional experimental results

In Sec. IV we mentioned some experimental results concerning the demonstration of the robustness of our controller in different possible scenarios and tracking different trajectories. Here we provide additional results and details of the presented experimental. Considering the experimental set-up described in Sec. IV, in the following we provide additional results for the four experiments:

1) Optimal conditions: the model parameters are the one measured with no additional errors (non-zero measurement errors and noise are still present). External disturbances were not applied neither. 2) Model error: an additional model error in the length of the cables equal to the 10% of their nominal value is added. 3) External disturbance: an artificial force equal to [1.7 1.7 1.7] [N] is added to the low level controller of the aerial vehicles. This simulates the effects of external disturbances acting on the vehicles, like wind. 4) Track of a far from quasi-static trajectory: we asked the system to accelerate the platform up to 3.2 [m/s 2 ] reaching a velocity up to 1.6 [m/s]. Our intent is to provide, for each experiment, the following results related to:

1) comparison between q d and q; 2) tracking errors;

3) desired and actual vehicle velocities, v R and v R , respectively. 1) Optimal Conditions: The model parameters are the one measured with no additional errors (measurement errors and noise are still present). External disturbances were not applied neither. From Fig. 9 one can observe that the controller is able to steer the system to the desired configuration with a error less than 0.02[m] in position, 2 • in the attitude and 2 • in the cables angle, α.

2) Model error: From Fig. 2 we can observe that the tracking error does not go perfectly to zero due to parameter uncertainties, noisy measurements and external disturbances. These non-idealities affects the exactness of the nominal model and the velocity tracking of the aerial vehicles. However, looking at e q and γ w , it is clear that condition of Theorem 1 is verified and the maximum H ∞ performance is guaranteed.

3) Simulated Constant Disturbance: We tested our system when an an artificial force equal to f ext = [1.7 1.7 1.7] [N] is added to the low level controller of the aerial vehicles:

u i = k R (v Ri -v Ri ) + f ext m Ri (32)
where m Ri is the mass of the i-th vehicle.

As a consequence of this external disturbance, one can observe from Fig. 11 that the position error of the platform increases. Nevertheless, the bounds are largely respected and ||e q || remains bounded by γ||w||.

4) Fast Trajectory Experiment: With this experiment, we tested the robustness of the proposed method when the assumption of quasi-static trajectory is largely violated, i.e., when the system is asked to follow a dynamic trajectories. In doing so, we aim at understanding the limits of our method. Nevertheless, we remark that our focus on accurate tasks where slow trajectories are preferred with respect to fast ones. The trajectory is designed to take the system from the starting configuration to a successive configurations of interest, q d (T 1 ) = [0 -0 0.4 -1.5 • 6 • 0.5 • 45 • 51 • 70 • ], and q d (T 2 ) = [0.28 -0.06 0.7 8

• 13 • 53 • 32 • 20 • 62 • ]
where T 1 = 42 and T 2 = 48. As the reader can see in Fig. 12 the developed controller is able to quickly follow, once again, the given reference trajectory trying to minimize the error. Even though the ||e q || increases as expected during the dynamic part of the trajectory, it still remains bounded by γ||w|| and goes back to smaller values when the trajectory is static.

Fig. 1 :

 1 Fig. 1: Cable-suspended load configurations: (a) Three robots and three cables [7]; (b) Six robots and six cables [16]; (c) Three robots and six cables (the Fly-Crane).

1 .

 1 W p L = [p Lx , p Ly , p Lz ] ∈ R 3 , describing the position of O L with respect to (w.r.t.) F W , and by the Euler angles W η L = [φ, θ, ψ] ∈ R 3 describing the orientation of F L w.r.t. F W We complete the description of the load state with its linear velocity v L = dp L /dt ∈ R 3 w.r.t. F W , and the Euler angles velocities ηL ∈ R 3 . The i-th cable is attached at one end to the load at point B i and at the other end to an aerial vehicle at point O Ri . The position of B i and O Ri w.r.t. F W are described by the vectors b i ∈ R 3 and p Ri ∈ R 3 , respectively. Since each pair of cables (1, 2), (3, 4) and (5, 6) is attached to the same point of the same aerial vehicle, we have that O R1 ≡ O R2 , O R3 ≡ O R4 and O R5 ≡ O R6 .

1 Fig. 2 :

 12 Fig. 2: Schematic representation of the Fly-Crane and its main variables.

[d 1 d 3 d 5 ]

 5 we can write the dynamics of the three aerial vehicles as v R = ṗR , vR = u + d.

  p d L = [0.28 -0.06 0.7] [m], η d L = [8 13 53] [deg], α d = [32 20 62] [deg], and T = 20 [s]. The gain of the velocity controller is set to k R = 15 (this value provides the smallest velocity error considering our low-level controller).

Fig. 4 :

 4 Fig. 4: Experimental screenshot of the Fly-Crane flying while controlling the platform in full pose (position plus orientation).

2 ) 3 )

 23 Model error: an additional model error in the length of the cables equal to the 10% of their nominal value is added. The corresponding results are shown in Fig. 5, Fig. 6 and Fig. 7. External disturbance: an artificial force equal to [1.7 1.7 1.7] [N] is added to the low level controller of the aerial vehicles. This simulates the effects of external disturbances acting on the vehicles, like wind. The main corresponding results are shown in Fig. 8. Additional plots are also available in the attached multimedia material together with a video of the conducted experiments.

Fig. 5 :

 5 Fig.5: Verification of conditions (12)-(14) for experiment 2).

Fig. 6 :

 6 Fig. 6: Tracking error for experiment 2).

Fig. 7 :

 7 Fig. 7: Tracking of the desired trajectory during system's evolution for experiment 2).

Fig. 8 :

 8 Fig.8: Tracking error for experiment 3). At time 52 [s] we introduced the virtual disturbance. This cause the increasing of the tracking error in position.

  Velocity tracking for the aerial vehicles.

Fig. 9 :

 9 Fig. 9: Plots for experiment 1)

  Velocity tracking for the aerial vehicles.

Fig. 10 :

 10 Fig. 10: Plots for experiment 2)

The left superscript indicates the reference frame. From now on, F W is considered as reference frame when the superscript is omitted.

All the non-idealities and parametric uncertainties in the model will be dealt with in the robust control design, and validated in real experiments. Furthermore, the tautness is ensured using a motion planner, like the one in[START_REF] Manubens | Motion planning for 6-D manipulation with aerial towed-cable systems[END_REF], that provides configurations of reference that always ensure taut cables.Preprint version, final version at http://ieeexplore.ieee.org/

In ∈ R n×n is the identity matrix of dimension n.Preprint version, final version at http://ieeexplore.ieee.org/ 3 IEEE Robotics and Automation Letters 2020

Preprint version, final version at http://ieeexplore.ieee.org/ 4 IEEE Robotics and Automation Letters 2020

Preprint version, final version at http://ieeexplore.ieee.org/

IEEE Robotics and Automation Letters 2020

Preprint version, final version at http://ieeexplore.ieee.org/ 8 IEEE Robotics and Automation Letters 2020

Preprint version, final version at http://ieeexplore.ieee.org/ 10 IEEE Robotics and Automation Letters 2020

Preprint version, final version at http://ieeexplore.ieee.org/ 12 IEEE Robotics and Automation Letters 2020

This research was partially supported by the ANR, Project ANR-17-CE33-0007 MuRoPhen, and by the European Union's Horizon 2020 research and innovation programme under grant agreement ID: 871479 AERIAL-CORE