Mark Stam
email: mstam@genoscope.cns.fr

Jordan Langlois

Céline Chevalier

Guillaume Reboul

Karine Bastard

Claudine Médigue

David Vallenet

Génomique Labgem

Métabolique

Cea

NetSyn: genomic context exploration of protein families

Background:

The growing availability of large genomic datasets presents an opportunity to discover novel metabolic pathways and enzymatic reactions profitable for industrial or synthetic biological applications. Efforts to identify new enzyme functions in this substantial number of sequences cannot be achieved without the help of bioinformatics tools and the development of new strategies. The classical way to assign a function to a gene uses sequence similarity. However, another way is to mine databases to identify conserved gene clusters (i.e. syntenies) as, in prokaryotic genomes, genes involved in the same pathway are frequently encoded in a single locus with an operonic organisation. This Genomic Context (GC) conservation is considered as a reliable indicator of functional relationships, and thus is a promising approach to improve the gene function prediction. Methods. Here we present NetSyn (Network Synteny), a tool, which aims to cluster protein sequences according to the similarity of their genomic context rather than their sequence similarity. Starting from a set of protein sequences of interest, NetSyn retrieves neighbouring genes from the corresponding genomes as well as their protein sequence. Homologous protein families are then computed to measure synteny conservation between each pair of input sequences using a GC score. A network is then created where nodes represent the input proteins and edges the fact that two proteins share a common GC. The weight of the edges corresponds to the synteny conservation score. The network is then partitioned into clusters of proteins sharing a high degree of synteny conservation. Results. As a proof of concept, we used NetSyn on two different datasets. The first one is made of homologous sequences of an enzyme family (the BKACE family, previously named DUF849)

Introduction

The advance of genome sequencing technology has created an ever increasing gap between protein sequences and annotations [START_REF] Galperin | From complete genome sequence to "complete" understanding?[END_REF]. It is currently estimated that 25% of the families defined in the PFAM database are annotated as "unknown function" [START_REF] Mudgal | De-DUFing the DUFs: Deciphering distant evolutionary relationships of Domains of Unknown Function using sensitive homology detection methods[END_REF]. Moreover, this issue is worsened by current automatic annotation procedures that still have a high false positive rate. This is especially true for methods that rely only on sequence similarity which lead to some protein families having a level of mis-annotation as high as 80% [START_REF] Schnoes | Annotation error in public databases: misannotation of molecular function in enzyme superfamilies[END_REF]. To overcome the lack or error annotation, innovative strategies use a combination of genome analysis and metabolic information [START_REF] Chen | Genome-wide discovery of missing genes in biological pathways of prokaryotes[END_REF][START_REF] Kharchenko | Identifying metabolic enzymes with multiple types of association evidence[END_REF][START_REF] Smith | The CanOE strategy: integrating genomic and metabolic contexts across multiple prokaryote genomes to find candidate genes for orphan enzymes[END_REF][START_REF] Yamanishi | Prediction of missing enzyme genes in a bacterial metabolic network. Reconstruction of the lysine-degradation pathway of Pseudomonas aeruginosa[END_REF], an integrative approach associating sequence family classification and controlled vocabulary [START_REF] Jung | A novel method for functional annotation prediction based on combination of classification methods[END_REF], Protein Protein Interaction (PPI) network [START_REF] Cozzetto | Protein function prediction by massive integration of evolutionary analyses and multiple data sources[END_REF][START_REF] Peng | Improving protein function prediction using domain and protein complexes in PPI networks[END_REF][START_REF] Zhao | An efficient method for protein function annotation based on multilayer protein networks[END_REF].

A way to circumvent the challenge of the gene function annotation, is to compare the genomic context (GC) of genes which is considered as a reliable indicator of functional relationships [START_REF] Huynen | Gene and context: integrative approaches to genome analysis[END_REF][START_REF] Janga | Nebulon: a system for the inference of functional relationships of gene products from the rearrangement of predicted operons[END_REF][START_REF] Mcclean | Synteny mapping between common bean and soybean reveals extensive blocks of shared loci[END_REF][START_REF] Rogozin | Connected gene neighborhoods in prokaryotic genomes[END_REF] and has been used to support functional annotation [START_REF] Ferrer | A systematic study of genome context methods: calibration, normalization and combination[END_REF][START_REF] Lee | Synteny Portal: a web-based application portal for synteny block analysis[END_REF][START_REF] Vallenet | MaGe: a microbial genome annotation system supported by synteny results[END_REF]. GC methods can be classified into four categories [START_REF] Ferrer | A systematic study of genome context methods: calibration, normalization and combination[END_REF]: gene neighbour, gene cluster, gene fusion (or Rosetta Stone) andphylogenetic profile and . Gene cluster (or conserved synteny) method is defined as the conservation of chromosomal proximity between genes [START_REF] Overbeek | Use of contiguity on the chromosome to predict functional coupling[END_REF][START_REF] Tamames | Evolution of gene order conservation in prokaryotes[END_REF]. It is based on conserved genes through different organisms and computes distances, in number of bases or genes, between two adjacent genes that are transcribed in the same direction [START_REF] Overbeek | Use of contiguity on the chromosome to predict functional coupling[END_REF]. Such distances are considered as a good predictor for co-regulated groups of genes like in operon structure [START_REF] Brouwer | The relative value of operon predictions[END_REF] or polysaccharide utilisation locis (PUL) [START_REF] Bjursell | Functional genomic and metabolic studies of the adaptations of a prominent adult human gut symbiont, Bacteroides thetaiotaomicron, to the suckling period[END_REF]. The Phylogenetic profile methods [START_REF] Pellegrini | Assigning protein functions by comparative genome analysis: protein phylogenetic profiles[END_REF] create a presence/absence matrix of homologous genes through different genomes. The length of the matrix depends on the number of considered genomes. It assumes that genes having a similar phylogenetic profile are more likely to share the same fonction. The gene fusion [START_REF] Marcotte | Detecting protein function and protein-protein interactions from genome sequences[END_REF] (or rosetta) method is based on genes fusion events. The idea is that two genes can be fused into some genomes. Pairs of proteins linked by such a fusion relationship are more likely to be functionally linked and occur in the same pathway. The gene neighbour methods [START_REF] Bowers | Prolinks: a database of protein functional linkages derived from coevolution[END_REF]) calculate the distance between two genes in different genomes. If one considers two genes from a genome, the distance between them is the number of genes between them plus one. This distance can be computed into different reference genomes. If the distance increases that means that the two genes are not linked. Two linked genes are more likely to be functionally in interaction. But in order to avoid a taxonomy bias, one must choose carefully the studied genomes which must not be too close in the taxonomy (i.e. organisms from the same genus).

Other efforts have been carried out to incorporate information from gene context to integrative approach. For instance, the STRING database [START_REF] Szklarczyk | The STRING database in 2017: quality-controlled protein-protein association networks, made broadly accessible[END_REF] stores all the known and predicted protein-protein interactions from different sources. These interactions can be physical (i.e. heterodimer) or indirect (e.g. gene regulation, signal relay mechanism) and can be represented by a graph where nodes correspond to the proteins and the edges to the interaction. This network highlights the interactions between different protein families. Another example is the Enzyme Function Initiative (EFI) [START_REF] Zallot | The EFI Web Resource for Genomic Enzymology Tools: Leveraging Protein, Genome, and Metagenome Databases to Discover Novel Enzymes and Metabolic Pathways[END_REF] which has developed a method based on the Genomic Neighborhood Network (GNN). Starting with a list of query sequences, a Sequences Similarity Network (SSN) is created, in which groups of highly connected nodes are considered as a cluster. Then, the genomic context is explored by recovering neighbour genes from the European Nucleotide Archive (ENA) [START_REF] Kanz | The EMBL Nucleotide Sequence Database[END_REF].Finally, neighbour genes are considered related if they belong to the same PFAM family [START_REF] Finn | Pfam: the protein families database[END_REF]. As output, EFI gives a Genome Neighbourhood Diagram (GNDs) for visual representation of the neighbourhoods for the sequences in each SSN cluster. In another study, gene context clustering had been integrated in a mega-clustering approach, which also included active site modeling, sequence similarity and phylogeny. Applied to a family of unknown function PFAM DUF849, forteen potential new activities had been predicted and experimentally validated subsequently [START_REF] Bastard | Revealing the hidden functional diversity of an enzyme family[END_REF].

None of the methods cited above have developed a sequence network that reflects only the GC conservation while avoiding the sequence similarity. A sequence similarity steps limit the input to only one homologous family at once, and make them unable to handle a set of non homologous proteins. Mixing non-homologous proteins as input is necessary when one wants to study genomic structures like operon, in which several genes involved in a specific cellular function (i.e., metabolic pathway, regulation, etc..) are co-localized in the same locus. With a method based only on the synteny conservation, genes belonging to the same structure will be gathered into the same cluster without having any sequence similarity.

Here, we present NetSyn, for Network Synteny, a tool allowing to link sequences based on the conservation of their genomic context in a network. NetSyn also offers several functionalities to explore the synteny network such as clustering the sequences, visualizing the genomic context of each sequence or mapping data from diverse resources. Our method quantifies genomic context conservation of proteins taken from different genomes. The advantage of NetSyn over others methods, is the fact that the sequence comparison is based not on sequence similarity but only on genomic context conservation. This important originality allows the genomic comparison of both a set of homologous sequences like a whole sequence family, as well as a set of non-homologous sequences like several proteomes. NetSyn was tested on two kinds of dataset. One set contains sequences from an homologous family. NetSyn was able to create iso-functional sub-families and pinpoint putative new pathways. The second set contains sequences from three non-homologous families, involved in the same degradation pathway. NetSyn was able to retrieve this known pathway in different organisms. Widely, the NetSyn principle is based on the idea that proteins sharing a similar genomic context are more likely to share similar functions or being involved in the same metabolic pathway. Consequently, our method can be extended to various sets of proteins.

Materials & Methods

Workflow and principles of NetSyn

NetSyn is a workflow developed in Python 3 and is divided into four steps: i) downloading the data, ii) extracting the genomic context, iii) computing the synteny conservation, iiii) computing a synteny network. NetSyn accepts as input a list of UniProt accession numbers (hereafter called target proteins). The workflow makes requests to the UniProt rest API to retrieve the corresponding nucleic accession number and then download the genome assembly file (hereafter called the genome file) from the European Nucleotide Archive (ENA) web site. Genomes files are used to extract the GC for each target protein (Fig. S1). But one can use his own data that are not stored in UniProt and ENA. In that case, NetSyn takes as input a correspondence file (Fig. S2), which is in a tabulation-separated value format with 5 mandatory fields: -protein_AC: the protein accession number in the genome file -protein_AC_field: the name of the field (qualifier) from which the protein accession number can be found in the genome file -nucleic_AC: the accession number of the genome file -nucleic_File_Format: the format of the genome file (EMBL or Genbank format) -nucleic_File_Path: the path to the locally stored genome file.

In addition to these files, NetSyn can accept a user-supplied metadata file to give attributes to target proteins. These attributes can be displayed on the final synteny network (example is given in Fig. S3). NetSyn parses the different genome files to retrieve the sequence of the target proteins and also the protein sequence of the neighbour genes. The number of retrieved neighbour genes depends on the window size parameter. It is set to 11 by default and means that NetSyn takes into account 5 genes before and 5 genes after the target gene. NetSyn also recovers, from the genome files, the full taxonomic lineage of the organisms having the target proteins. This taxonomic information is reported on the final graph and can be used to remove the redundancy if the graph is too dense (see the paragraph "Merging nodes"). All the retrieved proteins are then clustered by MMseqs2 [START_REF] Steinegger | MMseqs2 enables sensitive protein sequence searching for the analysis of massive data sets[END_REF] based on their sequence identity and alignment coverage percentages. These parameters can be modified and are set to 30% of identity on 80% of coverage by default, respectively. These MMseqs2 clusters allow NetSyn to define homologous protein families that are used to compute the synteny conservation for each pair of target proteins (Fig. 1).

Synteny computation and scoring

The computation of the genomic context conservation (i.e. synteny conservation) between two target proteins is made by an exact graph-theoretical approach described by Boyer et al. [START_REF] Boyer | Syntons, metabolons and interactons: an exact graph-theoretical approach for exploring neighbourhood between genomic and functional data[END_REF]. For each pair of target proteins in two different genomes, namely GA and GB, the graph-theoretical approach defines two networks, namely N1 and N2 where the vertices are the genes and the edges represent the fact that two genes are physically linked (i.e. on the same chromosome) or are homologous. Then, it searches connected components which are a subset of a network where any two vertices are connected to each other by paths but have no edges with the rest of the network. Finally a Common Connected Component (CCC) is computed taking into account what we call "gaps". Here we consider a "gap" in genome A as a gene in genome A which has no homologous gene into the genome B. A CCC is composed of common vertex/edge relationships in connected components in N1 and in N2. Figure 2 (panel A) shows an example of the computation of a CCC. The genome A has a connect component composed with the following relationships:

({G A5 , G B5 }, {T a , T b } {G A6 , G B6 } {G A7 , G B6 } {G A8 , G B8 }
) and the genome B has a connect component composed with these relationships:

({G A1 ,G B2 } {G A5 , G B5 } {T a , T b } {G A6 , G B6 } {G A7 , G B6 } {G A8 , G B8 }). The resulting CCC is composed by relationships present in the two lists: ({G A5 , G B5 } {T a , T b } {G A6 , G B6 } {G A7 , G B6 } {G A8 , G B8 }) (Fig. 2 panel B).
The {G A1 , G B2 } relationship is not taken into account in GA, because the number of genes between G A1 and T A which is superior to the threshold ( gap parameter set to 3 by default). The final CCC corresponds to the synteny considered by NetSyn.

Figure 2

Once a synteny has been defined between two target proteins, a synteny score is then calculated as follow:

𝑆𝑦𝑛𝑡𝑒𝑛𝑦 𝑆𝑐𝑜𝑟𝑒 = (𝐺𝑆/2) * (𝐺𝑆/𝐺𝑇)
Where GS is the number of genes having an homologous gene into the opposite genome, and GT the total number of genes in the considered synteny. In Figure 2 (panel B), the considered synteny contains 9 genes in synteny (or 9 genes having an homologous gene into the opposite genome) and 12 genes in total. Therefore the synteny score is SS = (9/2) * (9/12) = 3.375. In order to consider only significant synteny, by default, NetSyn keeps only synteny with a score equal to or higher than 3. This threshold can be modified by the user.

Synteny network and clustering

Once all syntenies and scores are computed between all the target proteins, NetSyn constructs a synteny network where nodes represent the target proteins. Two target proteins are linked with an edge if a conserved synteny is detected between them. In order to keep only significant syntenies, NetSyn considers only syntenies with score higher than a threshold defined by the user (3 by default) for the generation of the synteny network. Target proteins that are unrelated to any others target proteins are not taken into account for synteny network construction. The weight of the edge is equal to the synteny score which makes it possible to highlight the most conserved synteny.

In order to bring out target proteins sharing the most similar genomics context, several network clustering algorithms are implemented into netsyn : MCL, Walktrap, Louvain and Infomap. The Markov Cluster algorithm [START_REF] Dongen | Graph clustering by flow simulation[END_REF] (MCL) is an unsupervised algorithm based on simulation of flow in graphs and simulates random walks within a graph. The Walktrap algorithm [START_REF] Pons | Computing Communities in Large Networks Using Random Walks[END_REF] is the only implemented algorithm which takes into account the weight of the edge. It is also based on a random walk but unlike MCL, it can handle large graphs (> 1000 vertices). Louvain [START_REF] Blondel | Fast unfolding of communities in large networks[END_REF] is based on the comparison between the edge's density inside and outside a community (i.e. a cluster). Infomap [START_REF] Rosvall | Maps of random walks on complex networks reveal community structure[END_REF] is also based on a randomized walk trajectory and tries to minimize the length of the trajectory into a cluster. All the implemented algorithms can be applied on the final network. The clusters thus created are considered as the NetSyn cluster.

NetSyn network visualization

The final NetSyn's output is a network. We provide two different files format to explore and analyse the network : a graphML file which can be opened with network visualisation tools like Gephi or Cytoscape, and a HTML file, generated with the 3d-force-graph library, and which can be viewed with web browser (Fig. 3).

Figure 3

The HTML view is divided into 4 panels. The synteny network is displayed on the upper central panel. The upper left panel is the legend of the network and displays all attributes linked to each protein target like the NetSyn cluster they belong to, the metadata provided by the user and the different taxonomic rank retrieved during the genomic file parsing step. With the control box on the top, the user can choose which attributes to display. On this panel, if a user clicks on one attribute, all nodes sharing the same attribute are highlighted on the network and a schematic view of the context appears on the upper right panel. Each context is centered on the target protein with the number of genes before and after depending on the window size parameter given in input. In the example Fig. 3, the selected window size was 11, therefore 5 genes before and 5 genes after are shown. In this view, homologous genes are colored with the same color. If the user pass his cursor above the one the gene, a pop up appears and shows information from the NetSyn computation and the genome file: the family id (MMseqs2 family identifier), strain name, protein accession number, Uniprot accession number, gene name, locus tag, product and EC_number (Fig. S4). The lower panel gives information about the computation of homologous protein families by MMSEQ, for the selected NetSyn cluster. For each MMseqs2 family determined, it gives the number of selected proteins belonging to this family, the number of different species, strains and genomes from which the selected proteins come.

Merging nodes and network reduction

In some cases, the final network can be huge and tricky to analyse. To circumvent this challenge, we add an option to reduce the network. For that purpose, nodes identified as belonging to the same synteny cluster (several clustering methods are available) and sharing a given property will be merged into a unique node. The property can be a taxonomic rank, which is recovered by NetSyn when parsing the INSDC file, or a property given into the medatafile ( Fig. S5). The merging can be applied on only one property at a time in a single NetSyn run. The called "merged nodes" are the result of the union of one NetSyn cluster and one property. Merged nodes and the associated network cannot be rebuilt after a NetSyn computation. Only one clustering method must be provided.

Results

Dividing homologous family

NetSyn was tested on the β-keto acid cleavage enzymes family (BKACE). This family, initially called DUF849 (Domain of Unknown Function) in the PFAM database, was extensively characterized using high-throughput enzymatic screening [START_REF] Bastard | Revealing the hidden functional diversity of an enzyme family[END_REF]. The profiles of activity issued from the enzymatic screening of 124 representatives on 20 substrates showed a high correlation with the profiles of active sites. Indeed, the 725 proteins of the BKACE family were divided into 7 main groups based on structural classification of active sites by the Active Sites Modeling and Clustering (ASMC) method (de [START_REF] De Melo-Minardi | Identification of subfamily-specific sites based on active sites modeling and clustering[END_REF]. As ASMC groups show a good match with in vitro enzymatic activities, our aim here is to verify NetSyn clustering can identify the ASMC groups, and thus check if NetSyn can cluster iso-functional enzymes.

one CPU to compute. Among the 725 entries, 159 cannot be associated with an ENA file, because they were obsolete UniProt entries. Among the remaining 566 entries, 86 have no relevant conserved genomic context (i.e. with a synteny score >=3). Therefore 480 input entries (85% of the entries that had an corresponding ENA file) are retrieved into the final network. The clustering algorithm which gives about the same number of clusters as ASMC´s grouping was the Louvain algorithm. It gave 33 clusters with sizes varying from 2 to 68 sequences (Table 1). It should be noted that 12 clusters contain only 2 proteins. These clusters are considered as singleton and present low interest to analyse as they are not gathering enough context genomic information to highlight a metabolic pathway. Most of the sequences (391) are clustered into 9 Louvain clusters (0, 1,5,7,12,14,19,22,25) with a number of sequences varying from 23 to 68. All these clusters show a good fit with the ASMC groups (named here after G1 to G7) meaning that between 77 to 100% of their members belong to the same ASMC group. The only exception is the Louvain cluster 12 where only 55% of its sequences belong to the group G1. Cluster 12 is the only cluster that contains sequences from all the ASMC groups (Fig. 4). The network generated by NetSyn on BKACE data is available for visualization at https://doi.org/10.5281/zenodo.5746499 By definition, NetSyn clusters proteins according to their metabolic context while ASMC gathers proteins according to their active sites, which directly reflects the type of substrate they can catalyze (negative, positive, polar or hydrophobic substrates). Group G7 is the perfect example for illustrating this difference. Proteins which lacked at least one catalytic active site residue had been grouped in the G7 [START_REF] Bastard | Revealing the hidden functional diversity of an enzyme family[END_REF] and have been shown not transforming any keto-acid from the enzymatic screening. They have been annotated as non-BKACE enzymes. Thus, the G7 group is a heterogeneous group which functions could not be studied in our previous study [START_REF] Bastard | Revealing the hidden functional diversity of an enzyme family[END_REF]. However, NetSyn provides enough clues to investigate new metabolic pathways for G7 proteins. Proteins from the G7 group are distributed in 4 large Louvain clusters (Table 1) and show a large diversity of genomic context (Table S1). More than half the G7 enzymes are found in cluster Louvain 25 (Fig. S6), with a locus integrating a quine-oxidoreductase, a taurine dioxygenase, an MSF transporter, a porin, a fatty acid-CoA ligase and a transcriptional regulator (Table S1). This protein network has been described in Rhodococcus rhodochrous (actinobacteria class) (https://string-db.org/network/1429046.RR21198_4245). Other examples show that NetSyn performs better at splitting proteins according to the metabolic pathway rather than their active site. G1 was split into two highly populated NetSyn clusters (cluster 12 and 22) (Table 1). According to the genomic context, cluster 12 is involved in the conversion from valine to leucine while cluster 22 is involved in the pathway yet to be discovered, but for which clues are available (Table S1). The G5 group was also divided into high populated Louvain clusters (cluster 0 and 7). Their genomic contexts outline the two routes involved in the carnitine degradation (Table S1).

The benefit of NetSyn over ASMC is the lack of dependency on the molecular modeling threshold. 7 members previously grouped in G3 groups belong to Louvains cluster 25 (Table 1). Because the 3D model of these proteins presented all the catalytic residues, ASMC could not group them with the G7 group (i.e. at least missing one catalytic residue). Consequently, ASMC classified them in the G3 group for which there is not a particular substrate signature. However, according to their metabolic neighbors, these proteins are classified in cluster 25. Same remark occurs for four proteins classified in G1, G3 and G4 now classified in Louvain cluster 19 (Table S2). The genomic neighbouring of these four genes has in common two genes (Acetoacetate:butyrate CoA-transferase α subunit and β subunit) with other members of cluster 19 (Table S2). Contrastly, there are 8 proteins which used to belong to G7 group and now belong to cluster 12 (mostly populated with G1 group members). For instance, C0ZQU4 and A5V766 whose 3D model was of poor quality (sequence identity with the template for homology modeling at 23% and 19 %) were grouped into G7 because one catalytic residue was missing in the 3D model, whereas their metabolic function seems in a good agreement with Louvain cluster 1 and 22 (Table S3).

Louvain cluster 19 contains proteins from organisms with different phylum ranks: Clostridia, Fusobacteriia, Bacteroidia, Beta-proteobacteria and Gamma-proteobacteria, which belong to different phylums (Fig. S6). All the large clusters, as cluster 12, 1, 0, 14 and 25 present a large diversity of taxonomic classes. Only the highly populated cluster 5 with 30 members contains only alpha proteobacteria, mostly with the order Rhodobacterales. All genomes from the organisms of cluster 5 are similar. By consequence, the cluster 5 is not heterogeneous and genomic information repeated a few times. The group seems highly populated but actually the genomic information is repeated a few times. To avoid this skew, it appears necessary to prepare the data by removing the closest organism from the initial set before running NetSyn. Another way of detecting this skew is to use the merging mode of NetSyn with merging on the taxonomy. By comparing the NetSyn graph before and after merging, the user can detect what cluster is populated because of the diversity of taxonomy. We have also shown by looking at all clusters and the redundant genes retrieved, it is possible to highlight new pathways as shown in Table S1.

In our results, NetSyn was able to retrieve the ASMC groups and suggest even more detailed clusters. As ASMC had been shown to reflect in vitro activities, the agreement between these two methods indicates that NetSyn groups iso-functional enzymes. NetSyn clustering is more influenced by the composition of genomic neighborhoods than the taxonomic origin of the target proteins as most of the clusters defined here show a large taxonomic diversity. We have also demonstrated by curating manually the NetSyn clusters with literature that putative pathways can be deduced.

Interconnectivities of families in metabolic pathways

A great advantage of NetSyn is its ability to handle genes with no homology but functionally related. Indeed, half the genes of prokaryotes and a minor fraction of genes in eukaryotes form groups of physically clustered genes on chromosomes that are related in function [START_REF] Nützmann | Metabolic Gene Clusters in Eukaryotes[END_REF]. The most studied gene organization is the operon [START_REF] Jacob | Genetic regulatory mechanisms in the synthesis of proteins[END_REF] in which genes are transcribed as single polycistronic mRNA. Operons are often conserved across species by vertical inheritance, and thus operons can be computationally predicted.

More recently an operon-like structure, grouping CAZymes (for Carbohydrate Active enZymes), transcriptional regulators and transporters, has been highlighted [START_REF] Bjursell | Functional genomic and metabolic studies of the adaptations of a prominent adult human gut symbiont, Bacteroides thetaiotaomicron, to the suckling period[END_REF]. In Bacteroidetes, this system is called Polysaccharide Utilization Locus (PULs) and is composed of gene clusters that encode enzyme and protein ensembles required for the saccharification of complex carbohydrates. The first described PULs was centered around a transport system composed of two proteins : a TonB-dependent receptor, namely susC, and a carbohydrate binding protein, namely susD. Therefore, the different efforts to detect PULs were based on the research of the SusC/SusD couple followed by the inspection of the glycosyl hydrolase (GH) composition nearby [START_REF] Terrapon | Automatic prediction of polysaccharide utilization loci in Bacteroidetes species[END_REF]. But the SusC/SusD paradigme is challenged by the discovery of different sugar systems transport associated to PULs [START_REF] Larsbrink | A complex gene locus enables xyloglucan utilization in the model saprophyte Cellvibrio japonicus[END_REF][START_REF] Sheridan | Polysaccharide utilization loci and nutritional specialization in a dominant group of butyrate-producing human colonic Firmicutes[END_REF] or by the fact that the SusC/SusD system is not necessary to be localized into the gene cluster [START_REF] Ficko-Blean | Carrageenan catabolism is encoded by a complex regulon in marine heterotrophic bacteria[END_REF]. Consequently, in order to detect gene clusters or PULs, a method which is freed from the search for a particular gene and which highlights the colocalization of non homologous genes is necessary. NetSyn is able to handle genes with no homology but having similar GC, like genes involved in the xyloglucan PUL.

Xyloglucans (XyG) are members of the polysaccharides family found in terrestrial plant cell walls. They are composed of a (1→4)-glucan backbone appended with branching (1→6) β α xylosyl moieties (Fig. 5 panel A) [START_REF] Grondin | Polysaccharide Utilization Loci: Fueling Microbial Communities[END_REF]. Larbrink et al. [START_REF] Larsbrink | A complex gene locus enables xyloglucan utilization in the model saprophyte Cellvibrio japonicus[END_REF] have demonstrated that Cellvibrio japonicus strain Ueda107 can degrade the xyloglucan with a gene locus composed by 3 genes belonging to glycosyl hydrolase families GH31, GH35 and GH95 (genes xyl31A, bgl35A and afc95A respectively) and a transporter. Our goal is to focus on the interconnectivity between the GH proteins from different CAZy families that intervene in glycan degradation. For this purpose, we started to use the three key gene (xyl31A, bgl35A and afc95A) of Cellvibrio japonicus strain Ueda107 as reference to search homologous gene among the complete proteomes in UniProt using Blast (default parameter) [START_REF] Altschul | Gapped BLAST and PSI-BLAST: a new generation of protein database search programs[END_REF] We found 20 organisms containing homologs for the 3 key genes present in the Xygul in Cellvibrio Japonicus (Table S4). The complete proteome of these organisms were downloaded from UniProt and we used DBCAN [START_REF] Zhang | dbCAN2: a meta server for automated carbohydrate-active enzyme annotation[END_REF] to retrieve all the GHs. Then we selected all proteins annotated as GH31, GH35 or GH95 in these proteomes, which give a set of 194 UniProt entries and use this set as input for NetSyn. It took 47 second on one CPU to compute the final graph. The final network is available at https://zenodo.org/record/5705596. Only 136 input entries had a conserved genomic context and have been clusterized into 34 clusters with the Walktrap algorithm (Fig 5 It should be noted that the Cellvibrio japonicus Xygul locis lacks a xyloglucanase activity essential for the degradation of the xyloglucan backbone [START_REF] Larsbrink | A complex gene locus enables xyloglucan utilization in the model saprophyte Cellvibrio japonicus[END_REF]. However, this activity is found in the CAZY family GH74. In 7 organisms, the 3 key GHs are associated with a gene annotated as GH74 by DBCAN. Moreover, NetSyn unveils other genes that are involved in xylose metabolism or degradation generally. For instance, in the putative Xygul we have analysed, specific transporters (ABC sugar transporter, xylose ABC transporter) or enzymes (D-xylose-1-dehydrogenase) were found [START_REF] Stephens | Genetic Analysis of a Novel Pathway for d-Xylose Metabolism in Caulobacter crescentus[END_REF]. These findings emphasize that the locus we detected in the 14 organisms, gathered into NetSyn walktrap cluster 0, is dedicated to xyloglucan degradation. Organisms sharing a potential xyloglucan degradation locus belong to different classes as alphaproteobacteria, betaproteobacteria and gammaproteobacteria (Table 2). It shows that the conserved synteny in this Netsyn cluster is more due to a functional conservation rather than taxonomic conservation.

Automatic detection of the Xygul PUL in many different phyla was possible using NetSyn. Whereas in the litterature, this PUL is only described in βand ɣproteobacteria, our analysis demonstrated that this PUL is also potentially present in alphaproteobacteria (Table 2). Here we have shown that NetSyn is able to group together enzymes acting in interaction into the same degrading system. Functional characterisation is needed to confirm our prediction.

Discussion & Conclusions

We have shown that NetSyn is able to divide more precisely homologous sequences classified into families than a method based on actif site modelling as ASMC. NetSyn is also able to group sequences disregarding their origin (i.e taxonomic rank of the organism). NetSyn is able to pinpoint pathways yet to be discovered. These putative pathways, never observed in any organisms, could be guessed by the annotation provided by NetSyn in the genomic context of the target genes. Moreover, NetSyn is able to highlight groups of non homologous proteins but belonging to a pathway not yet observed into some organisms. We have also demonstrated that NetSyn is able to group together enzymes acting in interaction into the same degrading system but without any homologous relationships and sequence similarity. Their grouping is due to their genomic context conservation which reflects their functional conservation. But functional characterisation is needed to confirm our prediction.

NetSyn is a great tool to help functional assignment of protein sequences that are still of unknown function. Classical tools used sequence comparison to assign functions. However, even with recent improvement in homolog search accuracy [START_REF] Buchfink | Fast and sensitive protein alignment using DIAMOND[END_REF], there are still 23 % of protein families with uncharacterized function in databases [START_REF] Mistry | Pfam: The protein families database in 2021[END_REF]. On the other hand, many known activities have not been associated with a protein sequence. In their review, Sorokina et al., show that more than 1000 EC numbers, called orphan reactions, are not linked to a specific gene [START_REF] Sorokina | Profiling the orphan enzymes[END_REF]. To link genes with orphans reactions, methods based on protein sequence similarity are not suitable, but methods able to group proteins according to their function without taking into account sequence similarity is required. NetSyn, which gathers sequences depending on their context conservation, could be the ideal tool for tackling this issue.

Synteny conservation is one promising method that can be used to challenge the issue of linking genes with reaction. Also network generating and clustering allows to extract interesting information and organise results in a meaningful way. Recent efforts based on synteny conservation were published, like EFI-GNT which analysed up to hundred of genomes and generated vast amounts of data [START_REF] Zhao | An efficient method for protein function annotation based on multilayer protein networks[END_REF]. Despite the accuracy of the method, proper data verification is needed. EFI-GNT uses genomic context to predict gene function but such analysis may be biased by the fact that their network is based on sequence similarity. Also the context exploration is made through the PFAM classification which does not cover the entire protein sequence space. NetSyn, on the contrary, disregards database classification, but uses MMseqs2 to create its own homologous families. To our knowledge, NetSyn is the only tool able to handle a set of proteins without any common evolutionary origin whereas the other methods start from an homologous family. NetSyn paves the way to new and innovatives tools to solve the challenging task of predicting new enzymatic activities.

The main limitation of NetSyn is the fact that it may be biased by a selection of sequences from organisms classified at the same low taxonomic range (i.e organisms from the same genus). Indeed the lack of genomic diversity in the input set might lead to detected a synteny only due to closely related organisms rather than functional conservation through evolution. To circumvent this issue, we proposed to reduce the redundancy by merging nodes belonging to the same NetSyn clusters which share the same criteria like a taxonomic rank or a given metadata.

The main guiding idea of NetSyn is that genes sharing a common genomic context are more likely to share the same function. But this comparison must be independent of their sequence similarities. Here, the relationship between two proteins in the network depends only on the number of homologous genes into their corresponding genomic context. Therefore NetSyn is a suitable tool to explore either the functional diversity of enzyme families, either to highlight conserved loci containing genes of different evolutionary origins implied into the same pathway or in the same degradation/biosynthesis process. NetSyn is freely available at https://github.com/labgem/netsyn.

  panel A). Most of them are small clusters with only two or three sequences. The most populated cluster (cluster 0 in green Fig. 5 panel B) is also the only cluster where the three keys GH of the xyloglucan PUL are grouped (Fig. 5 panel C & D) including the GHs of the C. japonicus xyloglucan PUL.

Figure 5

 5 Figure 5

Table 1 :

 1 Comparison between ASMC groups and Louvain clusters

	ASMC\
	Louvain	G1	G2	G3	G4	G5 G6	G7	Total	max % couverture
	0	5	0	1	0	60	2	0	68	G5: 88%
	1	1	0	29	1	0	0	3	34	G3: 85%
	2	0	0	0	0	1	0	4	5	G7: 80%
	3	0	0	2	0	0	0	0	2	G3: 100%
	4	2	0	0	0	0	0	0	2	G4:100%
	5	0	0	0	30	0	0	1	31	G4: 97%
	6	0	0	3	0	0	0	0	3	G3: 100%
	7	9	0	4	1	47	0	0	61	G5: 77%
	8	0	0	0	0	0	7	0	7	G6: 100%
	9	0	0	0	0	0	0	10	10	G7:100%
	10	0	5	0	0	0	0	2	7	G2: 71%
	11	0	2	0	0	0	0	0	2	G2: 100%
	12	29	2	4	4	4	2	8	53	G1: 55%
	13	4	0	0	0	0	0	0	4	G1: 100%
	14	0	0	0	0	0	23	0	23	G6: 100%
	15	0	0	1	1	0	0	0	2	G3 : 50%
	16	0	0	2	0	0	0	0	2	G3: 100%
	17	0	3	0	0	0	0	0	3	G2: 100%
	18	0	0	0	2	0	0	0	2	G4: 100%
	19	1	27	2	1	0	0	0	31	G2: 87%
	20	1	0	1	0	0	0	0	2	G1: 50%
	21	0	0	0	0	0	6	0	6	G6: 100%
	22	32	0	0	1	0	0	1	34	G1: 94%
	23	2	0	0	0	0	0	0	2	G1: 100%
	24	0	0	0	0	0	6	0	6	G6: 100%

cluster 0. Of these 18 organisms, 14 have their 3 key GHs gathered into a conserved locus. For the 4 organisms lacking at least one key GH, it appears that one have an incomplete context (i.e the locus is interrupt by the end of the genome file), an another have an assembly gap in its assemblage, the two other are false positive due to duplication non essential genes (i.e transporter) (Table 2).


Filiberto Muñoz

Eduardo S Espinoza

Iván González

Sergio Salazar
email: sergio.salazar.cruz@gmail.com

Rogelio Lozano
email: rlozano@hds.utc.fr

E S Espinoza -I

Gonzalez
email: igonzalez@ctrl.cinvestav.mx

S Salazar -R

Lozano

Robust Trajectory Tracking for Unmanned Aircraft Systems Using a Nonsingular Terminal Modified Super-Twisting Sliding Mode Controller

Keywords: Nonsingular Terminal Modified Super Twisting, High Order Sliding Mode Observer, Unmanned Aerial System, RTK-GPS

The precision trajectory tracking problem for Unmanned Aerial Systems (UAS) is addressed in this work. A novel algorithm that combines a Nonsingular Modified Super-Twisting Controller with a High Order Sliding Mode Observer to enable an aerial vehicle tracking a desired trajectory under the assumption that i) its translational velocities are not available and ii) there are unmodeled dynamics and external disturbances. The proposed Sliding Mode Controller is based on a nonlinear sliding mode surface that ensures that the position and velocity tracking errors of all system's state variables converge to zero in finite time. Moreover, the proposed controller generates a continuous control signal eliminating the chattering phenomenon. Finally, simulation results and an extensive set of experiments are presented in order to illustrate the robustness and effectiveness of the proposed control strategy.

Introduction

The growing interest for Quadrotor aircraft development has been accelerated during the last decade, mainly due to the advantages over comparable Vertical Take Off and Landing (VTOL) UASs, such as helicopters. A Quadrotor aircraft has simple hardware structure combined with inexpensive and reliable sensors. Applications of multirotors range from surveillance, aerial photograph, precision agriculture, search and rescue, to acrobatic maneouvers for entertaining. Nowadays, the applications and use of UASs is increasing exponentially but, simultaneously, the challenges related to improving its performance are also increasing, especially in the field of precision navigation.

With the main objective of find a satisfactory performance of Quadrotor vehicles, different control strategies have been proposed, such as Robust Controllers based on Sliding Mode techniques. Sliding Mode Control has been widely explored in unmanned systems due to the property of insensitivity to parametric uncertainties, model uncertainties and external disturbances [START_REF] Parker | Decentralized sliding mode control for flexible link robots[END_REF], [START_REF] Espinoza | Modeling and Sliding Mode Control of a Micro Helicopter-Airplane System[END_REF]. Control strategies based on First Order Sliding Modes have shown to be an effective method for controlling dynamic systems with uncertainties, unfortunately, these First Order strategies induce signals with chattering due to the discontinuous terms in the control law. In order to overcome these drawbacks, techniques known as Higher Order Sliding Modes Control have been developed [START_REF] Emelyanov | Higher order sliding modes in the binary control systems[END_REF]- [START_REF] Edwards | Advances in Variable Structure and Sliding Mode Control[END_REF]. These High Order strategies have improved the controllers' performance with respect to chattering effects.

The Second Order Super-Twisting Controller has been the most popular High Order Sliding Mode Controller used for controlling Quadrotor aircrafts, see for example [START_REF] Luque-Vega | Robust block second order sliding mode control for a quadrotor[END_REF]- [START_REF] Jayakrishnan | Position and Attitude control of a Quadrotor UAV using Super Twisting Sliding Mode[END_REF]. In order to improve the convergence of the sliding variable in the Super Twisting algorithm, two linear correction terms are added to the classical Super Twisting to accelerate the approach of the sliding variable to the sliding surface. These modifications have been reported in [START_REF] Moreno | A Lyapunov approach to output feedback control using second-order sliding modes[END_REF]- [START_REF] Yang | A new modified super-twisting algorithm with double closedloop feedback regulation[END_REF]. The Sliding Mode Controllers based on linear sliding mode surface can only guarantee that closed-loop systems converge to the equilibrium asymptotically. In order to force the state of the system to reach the equilibrium in finite time, a terminal sliding mode control has been proposed.

Comparing the sliding mode control with the linear sliding surface, the terminal sliding mode presents important properties such as finite-time convergence. Although finite-time convergence can be guaranteed employing this method, the discontinuous signal will lead to the chattering phenomena caused by the switching control. Moreover, the terminal sliding mode has a singularity problem, that is, the terminal sliding mode control signal requieres to be infinitely large to maintain the ideal motion in some areas of the state space [START_REF] Fhen | On nonsingular terminal sliding-mode control of nonlinear systems[END_REF]. In [START_REF] Mou Chen | Terminal sliding mode tracking control for a class of SISO uncertain nonlinear systems[END_REF] a terminal sliding mode tracking control for a class of SISO uncertain nonlinear systems with unknown external disturbances is proposed. In [START_REF] Xhao | Distributed adaptive fixed-time consensus tracking for second-order multi-agent systems using modified terminal sliding mode[END_REF] a modified terminal sliding mode is used for solving a fixed-time consensus tracking for second-order multi-agent system.

In order to overcome the classical singularity problem of the conventional terminal sliding mode, the nonsingular terminal sliding mode has been pro-posed. Moreover, nonsingular terminal sliding mode has a better performance against the controllers based on linear sliding mode surfaces [START_REF] Zhu | Space-based line-of-sight tracking control of GEO target using nonsingular terminal sliding mode[END_REF]- [START_REF] Madani | Non-singular terminal sliding mode controller: Application to an actuated exoskeleton[END_REF]. In [START_REF] Fhen | On nonsingular terminal sliding-mode control of nonlinear systems[END_REF] a nonsingular terminal sliding-mode control for nonlinear systems is presented. In [START_REF] Zhu | Space-based line-of-sight tracking control of GEO target using nonsingular terminal sliding mode[END_REF] an adaptive nonsingular terminal mode controller is proposed for spacecraft attitude tracking. A Continuous nonsingular terminal sliding mode control for systems with mismatched disturbances is presented in [START_REF] Yang | Continuous nonsingular terminal sliding mode control for systems with mismatched disturbances[END_REF]. The Nonsingular Terminal Sliding Mode Control of Uncertain Second-Order Nonlinear System was stuied in [START_REF] Tran | Nonsingular Terminal Sliding Mode Control of Uncertain Second-Order Nonlinear Systems[END_REF]. A Non-singular terminal sliding mode controller applied to an actuated exoskeleton was developed in [START_REF] Madani | Non-singular terminal sliding mode controller: Application to an actuated exoskeleton[END_REF].

In real applications, the performance of the UAVs could be affected by external perturbations. In order to overcome this drawback, many research work are devoted to deal with the estimation of perturbations. In this sense, sliding mode observers provide a robust and finite-time convergence observer, while, high-order sliding mode observers preserves the same performance than standard sliding mode but without the undesired chattering effect. In [START_REF] Levant | Higher order sliding modes, differentiation and output-feedback control[END_REF] author proposes an finite-time arbitrary-order robust differentiator using highorder sliding mode, this differentiator could be used like a robust observer. Authors in [START_REF] Benallegue | High order sliding mode observer for a quadrotor UAV[END_REF] use a feedback linearization-based controller with a high-order sliding mode observer for a Quadrotor helicopter, where the high-order sliding mode observer is employed as an observer and estimator of the effect of the external disturbances such as wind and noise. In [START_REF] Héctor Ríos | High-order sliding mode observers for nonlinear autonomous switched systems with unknown inputs[END_REF] the authors proposed observers using high-order sliding mode for obtaining the reconstruction of the exact non-linear system even in presence of disturbances.

This work is motivated by the tracking control in finite time of an UAS type Quadrotor affected by external disturbances, unmodeled dynamics and partial measurement of the state. Compared with most of the existing sliding mode control methods employed to control a UAS type Quadrotor, the main contribution of this paper are summarized as follows:

• A novel algorithm using a Nonsingular Modified Super-Twisting Controller enable to an aerial vehicle to track a desired trajectory that force the tracking error to reach the zero equilibrium in finite time. • An output feedback Nonsingular Modified Super-Twisting Controller based on a High Order Sliding Mode Observer is developed, for case that only the system output is available. Also, it is show that the tracking error converges to zero in finite-time with the Controller-Observer proposed. In addition, a continuous control signal is obtained with the proposed strategy. • Practical implementation is conducted in outdoor environments using a high preceision RTK-GPS system to show the effectiveness of the control strategy proposed.

The remainder of the paper is organized as follows. The dynamic model of the UAS type Quadrotor is presented in Section 2. The development of the Nonsingular Terminal Modified Super-Twisting Controller based on a High Order Sliding Mode Observer is showed in Section 3. The application of the Nonsingular Terminal Controller to the aerial vehicle is carry out in Section 4. A set of simulation and experimental results to illustrate the performance of the control strategie proposed is presented in Section 5, followed by some concluding remarks in Section 6.

Dynamic Model

The quadrotor is an underactuated four-rotor unmanned aerial vehicle having six degrees of freedom and four control inputs including the total thrust and three angular moments. It has four rotors arranged in a cross shape. The front and rear motors rotate counterclockwise (CCW), while the other two motors rotate clockwise (CW). The main thrust is the sum of the thrusts of each motor, the pitch torque is a function of the difference f 1 -f 3 , the roll torque is produced by the difference f 2 -f 4 , and the yaw torque is the sum of τ M1 + τ M2 + τ M3 + τ M4 (see Figure 1), where τ Mi is the reaction torque of motor i due to shaft acceleration and the blades drag. Let Γ I = {i I , j I , k I } be the inertial frame, Γ B = {i B , j B , k B } denotes a set of coordinates fixed to the rigid aircraft.

Let q = (x, y, z, φ, θ, ψ) T ∈ R 6 = (ζ, η) T be the generalized coordinates vector which describes the position and orientation of the flying machine, so the model could be separated in two coordinate subsystems: translational and rotational. These are defined respectively by The dynamic model is obtained using the Euler-Lagrange approach. We can decompose the equations into translational and rotational displacement. The Lagrangian of the aerial vehicle is given by the following equation L(q, q) = (T trasl + T rot ) -U, where -T trasl : Is the translational kinetic energy of the quad-rotor

T trasl = 1 2 m ζT ζ.
-T rot : Is the rotational kinetic energy of the quad-rotor

T rot = 1 2 ηT J η.
-U : Is the potential energy of the quad-rotor

U = mgz.
In compact form, the dynamic model of the Quadrotor is given as:

m ζ = F I + mg Jη = τη + C(η, η) η.
Afterwards, doing a change of coordinates, we propose new variables to the control input τη :

τ =   τ φ τ θ τ ψ   = J -1 (τ η -C (η, η) η) ,
where F B can be expressed in the inertial frame F I as F I = R B→I F B , R B→I represents the rotation matrix which is defined by

R B→I =   cψcθ cψsθ -sψcφsφ cψsθcφ + sψsφ sψcθ sψsθsφ + cψcφ sψsθcφ -cψsφ -sθ cθsφ cθcφ   ,
where s(•) = sin(•) and c(•) = cos(•).

Finally, the dynamic model of the quadrotor can be described by the following equations:

ẍ = 1 m (cos φ sin θ cos ψ + sin φ sin ψ) u 1 - d 1 m ẋ (1a) ÿ = 1 m (cos φ sin θ sin ψ -sin φ cos ψ) u 1 - d 2 m ẏ (1b) z = 1 m (cos φ cos θ) u 1 -g - d 3 m ż (1c) φ = θ ψ I y -I z I x - ld 4 I x φ + l I x u 2 (1d) θ = φ ψ I z -I x I y - ld 5 I y θ + l I y u 3 (1e) ψ = φ θ I x -I y I z - d 6 I z ψ + 1 I z u 4 , (1f)
where x, y, and z represent the Quadrotor's translational position in axis X, Y and Z, respectively. The translational velocities and translational accelerations are represented as ( ẋ, ẏ, ż) and (ẍ, ÿ, z), respectively. The Euler angles roll, pitch and yaw are represented by φ, θ and ψ, respectively. The rotational velocities and rotational accelerations are given by φ, θ, ψ and φ, θ, ψ , respectively. The mass of the vehicle is denoted by m, the term g represents acceleration due to gravity, l is the distance of center of mass to the motors. The inertia moments of the vehicle are represented by I x , I y and I z . The control inputs are given by u i , i = 1, . . . , 4. The terms d i , i = 1 . . . 6, represent drag coefficients.

3 Nonsingular Terminal Modified Super Twisting Controller Based on a High Order Sliding Mode Observer

In this section we present the mathematical structure of the Nonsingular Terminal Modified Super-Twisting Sliding Mode Controller and the High Order Sliding Mode Observer, which are the base for the Controller-Observer proposed in this work. Consider a second order dynamic system defined as

χ1 (t) = χ 2 (t) χ2 (t) = f (χ (t)) + h (χ (t)) u (t) + ξ (t, χ) y (t) = χ 1 (t) , (2)
where

χ (t) = [χ 1 (t) , χ 2 (t)] T ∈ R 2 represents the state vector, y (t) ∈ R is the system output, u (t) is the input control signal, f (χ (t)) ∈ R and h (χ (t)) ∈ R
are known functions. Moreover, the term ξ (t, χ) ∈ R includes both unmodeled dynamics as external disturbances. We assume that disturbance ξ (t, χ) is Lipschitz and satisfies | ξ(t, ξ)| < ξ + for a positive constant ξ + . Hereafter, temporal dependence of the t variable will be omitted for simplicity.

High Order Sliding Mode Observer

The High Order Sliding Mode Observer used to estimate the dynamic χ 2 as well as the disturbance ξ of system ( 2) is given as

χ1 = χ2 + λ 1 | χ1 | 2 3 sign ( χ1 ) χ2 = χ3 + λ 2 | χ1 | 1 3 sign ( χ1 ) + f (χ) + h (χ) u χ3 =λ 3 sign ( χ1 ) , (3)
where χi , i = 1, . . . , 3, represent the estimated variables. In the subsequent, variables with symbol ˆdenote estimated variables. The estimation errors are defined as χ1 = χ 1 -χ1 , χ2 = x 2 -χ2 and χ3 = -χ3 + ξ. In the subsequent, variables with symbol ˜denote estimation error variables. Now, from ( 2) and (3), we can write the estimation error's dynamics as

χ1 = -λ 1 | χ1 | 2 3 sign ( χ1 ) + χ2 χ2 = -λ 2 | χ1 | 1 3 sign ( χ1 ) + χ3 χ3 = -λ 3 sign ( χ1 ) + ξ. ( 4
)
The dynamics of the estimation error (4) have the form of the non-recursive exact robust differentiator presented in [START_REF] Levant | Higher order sliding modes, differentiation and output-feedback control[END_REF]. The convergence proofs for (4) have been obtained by geometric methods in [START_REF] Levant | Higher order sliding modes, differentiation and output-feedback control[END_REF], by using homogeneity properties in [START_REF] Levant | Homogeneity approach to high-order sliding mode design[END_REF] and by using a quadratic and strict Lyapunov function in [START_REF] Moreno | Lyapunov function for Levant's Second Order Differentiator[END_REF]. Therefore, the estimation errors χ1 , χ2 , and χ3 will converge to zero in a finite time t ≥ T 0 if gains λ 1 , λ 2 , and λ 3 are chosen appropriately [START_REF] Levant | Higher order sliding modes, differentiation and output-feedback control[END_REF], [START_REF] Moreno | Lyapunov function for Levant's Second Order Differentiator[END_REF]. From the error variables definition, after finite time t ≥ T 0 , it results that χ 1 = χ1 , χ 2 = χ2 and χ3 = ξ.

Nonsingular Terminal Modified Super-Twisting Sliding Mode Controller

It is well known that the main drawback of the First Order Sliding Mode Control is the chattering phenomenon that appears in the control law. As a solution for solving this issue, authors in [START_REF] Emelyanov | Higher order sliding modes in the binary control systems[END_REF] presented Second Order Sliding algorithms, such as the Twisting algorithm, where the idea of acting on the superior derivatives of the sliding variable was introduced. On the other hand, in the case of the Second Order Sliding Mode Control, the following condition should be verified

s (χ) = ṡ (χ) = 0,
where s (χ) represents the sliding variable. In the remainder of the article, the dependence of χ on the s (χ) sliding variable will be omitted. The Modified Super-Twisting Sliding Mode Controller is a Second Order Sliding Mode (SOSM) algorithm introduced firstly by [START_REF] Moreno | A Lyapunov approach to second-order sliding mode controllers and observers[END_REF], and variations of this algorithm have been presented in [START_REF] Yang | A modified super-twisting sliding mode control with inner feedback and adaptive gain schedule[END_REF], [START_REF] Yang | A new modified super-twisting algorithm with double closedloop feedback regulation[END_REF], [START_REF] Defoort | A Lyapunov-based design of a modified supertwisting algorithm for the Heisenberg system[END_REF].

In the following subsections we introduce the tracking control scheme using the nonsingular terminal modified super-twisting sliding mode controller both for the case where all the state variables are available and when the state variables are measured only partially.

A) Nonsingular Terminal Modified Super-Twisting Sliding Mode Controller

A nonlinear dynamic sliding mode manifold for system (2) is defined as follows

s = e 2 + c 1 e 1 + c 2 exp -λt e (1-2β) 1 , (5)
where c 1 > 0, c 2 > 0, λ > 0 and 0 < β < 1. The sliding mode manifold (s = 0) guarantees the finite-time convergence of the error e 1 via the following dynamic equation

e 2 = -c 1 e 1 -c 2 exp -λt e (1-2β) 1
, where the tracking errors are defined as e 1 = χ 1 -χ 1 d and e 2 = χ 2 -χ 2 d , and the signal references are defined as χ 1 d and χ 2 d . The finite-time convergence of the sliding mode can be proved by the following theorem.

Theorem 1 [START_REF] Tran | Nonsingular Terminal Sliding Mode Control of Uncertain Second-Order Nonlinear Systems[END_REF] For system (2), considering the sliding manifold (5), the tracking error e 1 will converge to zero in finite time if 2βc 1 -λ > 0 and the finite convergence time is

T s ≤ ln 1 + exp 2βc1t V β (0) /a 2 2βc 1 -λ ,
where

a 2 = 2 (1-β) βc 1 2βc 1 -λ > 0 and V = 1 2 e 2 1 .
The temporal derivative of the sliding variable, with the dynamics of the system (2) is given as

ṡ = f (χ) + h (χ) u + ξ -χd 2 + c 1 ė1 + A, (6)
with A defined as

A = c 2 exp -λt e -2β 1 [(1 -2β) ė1 -λe 1 ] .
Proposition 1. Let the Nonsingular Terminal Modified Super Twisting control input u be defined as

u =h (χ) -1 -f (χ) + χd 2 -c 1 ė1 -A + c 1 χd 1 -k 1 |s| 1 2 sign (s) -k 2 s - t 0 k 3 sign (s) dτ - t 0 k 4 sdτ . (7)
By properly selecting the gains k 1 , k 2 , k 3 and k 4 , then s = ṡ = 0 in finite time, which implies that the tracking errors e 1 and e 2 will converge to zero in finite time.

Proof . Substituting the Nonsingular Terminal Modified Super Twisting Controller [START_REF] Besnard | Quadrotor vehicle control via sliding mode controller driven by sliding mode disturbance observer[END_REF] in ( 6) we obtain

ṡ = -k 1 |s| 1 2 sign (s) -k 2 s - t 0 k 3 sign (s) dτ - t 0 k 4 sdτ + ξ.
Defining new variables ω 1 and ω 2 as

ω 1 =s ω 2 = -k 3 t 0 sign (s) dτ -k 4 t 0 sdτ + ξ,
it can be obtained the following closed loop dynamic corresponding to the sliding surface

ω1 = -k 1 |ω 1 | 1/2 sign (ω 1 ) -k 2 ω 1 + ω 2 ω2 = -k 3 sign (ω 1 ) -k 4 ω 1 + ξ.
Suppose that for system (2) the derivative of the disturbance is globally bounded by | ξ| < ξ 1 + ξ 2 |x 1 | with ξ 1 > 0 and ξ 2 > 0, assuming that the gains k 1 , k 2 , k 3 , k 4 are selected according to

k 1 > ξ 1 k 2 > 1 2 8ξ 2 k 3 > ξ 1 k 4 > k 1 1 2 k 3 1 (2k 2 -ξ 2 ) + 5 2 k 2 2 + ξ 2 p 1 k 1 p 1 -1 2 k 3 1 ,
where

p 1 = k 1 1 4 k 2 1 -ξ 1 + 1 2 k 1 2k 3 + 1 2 k 2 1 .
Then, the Nonsingular Terminal Modified Super Twisting Controller (7) yields finite-time convergence of the sliding surface s = 0 [START_REF] Moreno | Lyapunov function for Levant's Second Order Differentiator[END_REF], [START_REF] Moreno | A Lyapunov approach to second-order sliding mode controllers and observers[END_REF] and from Theorem 1, the tracking errors e 1 and e 2 will converge to zero in finite time.

B) Nonsingular Terminal Modified Super-Twisting Sliding Mode Controller based on High Order Sliding Mode Observer

Assuming that the state variable χ 2 is unmeasured, the design of the Nonsingular Terminal Modified Super-Twisting Sliding Mode Controller will be based on the High Order Sliding Mode Observer presented in Subsection 3.1.

In this sense, a nonlinear dynamic sliding mode manifold for system (2) by using the estimation χ2 is defined as follows

s = ê2 + c 1 e 1 + c 2 exp -λt e (1-2β) 1
, where the tracking error ê2 is defined as ê2 = χ2 -χ 2 d . The temporal derivative of the sliding variable, with the dynamics of the system (2) is given as

ṡ = ė2 + c 1 ė1 + c 2 exp -λt e -2β 1 [(1 -2β) ė1 -λe 1 ] .
Substituting the dynamics of system (2) and observer (3) in the dynamics of the sliding manifold we obtain

ṡ = t 0 λ 3 sign ( χ1 ) dτ + λ 2 | χ1 | 1 3 sign ( χ1 ) + f (χ) + h (χ) u -ẋd 2 + c 1 χ2 + χ2 -ẋd 1 + c 2 exp -λt e 2β 1 (1 -2β) χ2 + χ2 -ẋd 1 -λe 1 . (8)
Proposition 2. Let the Nonsingular Terminal Modified Super Twisting control input u based on High Order Sliding Mode Observer (3) be defined as

u =h (χ) -1 - t 0 λ 3 sign ( χ1 ) dτ -λ 2 | χ1 | 1 3 sign ( χ1 ) -f (χ) + ẋd 2 -c 1 χ2 + c 1 χd 1 -c 2 exp -λt e 2β 1 (1 -2β) χ2 -χd 1 -λe 1 -k 1 |s| 1 2 sign (s) -k 2 s - t 0 k 3 sign (s) dτ - t 0 k 4 sdτ -χ3 . (9)
By properly selecting the gains k 1 , k 2 , k 3 and k 4 , then s = ṡ = 0 in finite time, which implies that by Theorem 1 the tracking errors e 1 and e 2 will converge to zero in finite time.

Proof . Substituting the Nonsingular Terminal Modified Super Twisting Controller based on High Order Sliding Mode Observer ( 9) in ( 8) we obtain

ṡ =c 1 χ2 + c 2 exp -λt e 2β 1 [(1 -2β) χ2 ] -k 1 |s| 1 2 sign (s) -k 2 s - t 0 k 3 sign (s) dτ - t 0 k 4 sdτ + ξ -χ3 .
Therefore, the closed-loop Controller-Observer system is given by the following expressions

Γ =              χ1 =s -c 1 e 1 -c 2 exp -λt e (1-2β) 1 + χ2 + x 2 d ṡ = χ2 c 1 + c 2 exp -λt e 2β 1 (1 -2β) -k 1 |s| 1 2 sign (s) -k 2 s + ξ -χ3 + v v = -k 3 sign (s) -k 4 s, Σ =        χ1 = -λ 1 | χ1 | 2/3 sign ( χ1 ) + χ2 χ2 = -λ 2 | χ1 | 1/3 sign ( χ1 ) + χ3 χ3 = -λ 3 sign ( χ1 ) + ξ. ( 10
)
The estimation error of observer [START_REF] Modirrousta | A novel nonlinear hybrid controller design for an uncertain quadrotor with disturbances[END_REF] converges to zero in finite-time. Consequently, we obtain for the error variables that χ1 = χ2 = χ3 = 0. After time T 0 , the closed-loop system is given by the dynamics

Γ =        χ1 = s -c 1 e 1 -c 2 exp -λt e (1-2β) 1 + x 2 d ṡ = -k 1 |s| 1 2 sign (s) -k 2 s + v v = -k 3 sign (s) -k 4 s. (11)
The last two terms of (11) represent a Modified Super Twisting Sliding Mode Controller introduced in [START_REF] Moreno | A Lyapunov approach to second-order sliding mode controllers and observers[END_REF], where choosing the gains as

k i > 0, i = 1, . . . , 4 4k 3 k 4 > 8k 3 + 9k 2 1 k 2 2 ,
then it is obtained s = ṡ = 0 in finite time. This implies that the tracking errors e 1 and e 2 are asymptotically stable. Finally, we obtain that χ 1 and χ2 converge in finite time to χ 1 d and χ 2 d , respectively.

Robust Controller for the aerial vehicle

In this section, a robust controller for position and attitude tracking of the Quadrotor helicopter will be discussed. Position and attitude robust control are designed by considering a time-scale separation between the translational dynamics and the orientation dynamics. This approach is based on the assumption that the closed-loop attitude dynamics converge faster than the closed-loop translational dynamics [START_REF] Bertrand | A hierarchical controller for miniature VTOL UAVs: Design and stability analysis using singular perturbation theory[END_REF], [START_REF] Liu | Robust Tracking Control of a Quadrotor Helicopter[END_REF].

From (1) it can be seen that the Quadrotor dynamics can be divided into two subsystems: i) an underactuated system composed by dynamics ẍ, ÿ, φ, θ and ii) a fully actuated system composed of dynamics z and ψ.

Control Problem Formulation.

For the Quadrotor dynamics (1), it is necessary to design control inputs u 1 , u 2 , u 3 and u 4 using a Nonsingular Terminal Modified Super-Twisting Sliding Mode Controller, such that:

-The position and velocity tracking errors of all system state variables converge to zero in finite. -For the fully actuated subsystem, it is required that yaw angle ψ and altitude z tracking robustly signals ψ d and z d . -For underactuated subsystem controller, it is necessary to guarantee that positions x and y tracking robustly signals x d and y d , respectively, whilst φ and θ tracking robustly φ des and θ des . -Desired references x d , y d , z d and ψ d are generated by a trajectory generator or an upper level controller, whereas references φ des and θ des are generated by using (1a) and (1b).

-Translational velocities ẋ, ẏ, ż cannot be measured online, so they will be estimated by a High Order Sliding Mode Observer.

Position Control

The position controller generates the desired references for the pitch and roll angles φ des and θ des , respectively, and simultaneously, stabilizes the altitude of the Quadrotor.

Altitude Controller. From (1c), we can express the altitude dynamics as

ż1 =z 2 ż2 = 1 m (cos φ cos θ) u 1 -g + ξ z , (12)
where z 1 = z and z 2 = ż. In this work, we assume that the term -d3 m ż is treated as a unmodeled dynamic and it is included in the term ξ z (external disturbances can also be included in ξ z ). In order to track a desired altitude z 1 d , we employ the Modified Super-Twisting Controller based on the High Order Sliding Mode Observer (3) presented in Section 3. The dynamics of observer (3) applied to the altitude dynamics [START_REF] González-Hernández | Robust trajectory-tracking control design for a small Quad-rotor aircraft via sliding modes[END_REF] are given by

ż1 =ẑ 2 + λ 1z |z 1 | 2 3 sign (z 1 ) ż2 =ẑ 3 + λ 2z |z 1 | 1 3 sign (z 1 ) + 1 m (cos φ cos θ) u 1 -g ż3 =λ 3z sign (z 1 ) , (13)
where z1 = z 1 -ẑ1 is the observation error of z 1 , ẑ1 and ẑ2 represent the altitude position and the estimated velocity, respectively. The term ẑ3 is the estimated of the disturbance ξ z . Let us define a nonlinear sliding manifold s z as

s z = ẑ2 -z 2 d + c z1 z 1 -z 1 d + c z2 exp -λzt z 1 -z 1 d (1-2βz) ,
where c z1 > 0, c z2 > 0, λ z > 0 and 0 < β z < 1. The control signal u 1 is given by

u 1 = m cos φ cos θ - t 0 λ 3z sign (z 1 ) dτ -λ 2z |z 1 | 1 3 sign (z 1 ) + g + żd 2 -c z1 ẑ2 + c z1 żd 1 -c z2 exp -λzt z 1 -z 1 d 2βz (1 -2β z ) ẑ2 -żd 1 -λ z z 1 -z 1 d -k 1z |s z | 1 2 sign (s z ) -k 2z s z - t 0 k 3z sign (s z ) dτ - t 0 k 4z s z dτ . (14)
Combining controller ( 14) together with observer (13), ensures that z 1 tracks in finite time to z 1 d under unknown disturbances and estimating on line the altitude velocity z 2 .

Horizontal Controller. The Quadrotor's translational dynamics in X axis represented by (1a) can be rewritten as ẋ1 =x 2 ẋ2 =µ x + ξ x , where µ x = 1 m (cos φ sin θ cos ψ + sin φ sin ψ) u 1 and ξ x = -d1 m ẋ (this term can also include external disturbances). Defining µ xv as a new virtual control input, to track the desired position x 1 d , the following sliding variable is defined

s x = x2 -x 2 d + c x1 x 1 -x 1 d + c x2 exp -λxt x 1 -x 1 d (1-2βx) ,
where c x1 > 0, c x2 > 0, λ x > 0 and 0 < β x < 1. The virtual control input µ xv is given as

µ xv = - t 0 λ 3x sign (x 1 ) dτ -λ 2x |x 1 | 1 3 sign (x 1 ) + ẋd 2 -c x1 x2 + c x1 ẋd 1 -c x2 exp -λxt x 1 -x 1 d 2βx (1 -2β x ) x2 -ẋd 1 -λ x x 1 -x 1 d -k 1x |s x | 1 2 sign (s x ) -k 2x s x - t 0 k 3x sign (s x ) dτ - t 0 k 4x s x dτ . (15)
Following the same procedure, for the dynamic in Y axis, defining µ y = 1 m (cos φ sin θ sin ψ -sin φ cos ψ) u 1 and ξ y = -d2 m ẏ. The dynamics given in (1b) can be rewritten as ẏ1 =y 2 ẏ2 =µ y + ξ y , the sliding variable s y is proposed as

s y = ŷ2 -y 2 d + c y1 y 1 -y 1 d + c y2 exp -λyt y 1 -y 1 d (1-2βy) ,
where c y1 > 0, c y2 > 0, λ y > 0 and 0 < β y < 1. The virtual control input µ yv is given as

µ yv = - t 0 λ 3y sign (ỹ 1 ) dτ -λ 2y |ỹ 1 | 1 3 sign (ỹ 1 ) + ẏd 2 -c y1 ŷ2 + c y1 ẏd 1 -c y2 exp -λyt y 1 -y 1 d 2βy (1 -2β y ) ŷ2 -ẏd 1 -λ y y 1 -y 1 d -k 1y |s y | 1 2 sign (s y ) -k 2y s y - t 0 k 3y sign (s y ) dτ - t 0 k 4y s y dτ . (16)
Virtual controllers ( 15) and ( 16) combined with High Order Sliding Mode Observer (3), guarantee that x 1 and y 1 tracks in finite time to x 1 d and y 1 d , respectively.

In order to force the virtual controllers to match µ x and µ y , i. e. µ xv = µ x and µ yv = µ y , it is necessary to obtain a relation between the virtual controllers and the angles φ and θ. Based on (1a) and (1b), we obtain the following expressions

µ xv cos ψ = 1 m cos φ sin θ cos 2 ψ + sin φ cos ψ sin ψ u 1 + cos ψ ξ x µ yv sin ψ = 1 m cos φ sin θ sin 2 ψ -sin φ cos ψ sin ψ u 1 + sin ψ ξ y .
The desired angles φ des and θ des (in function of virtual controllers µ xv and µ yv given by ( 15) and ( 16), respectively) that enables the system to track desired positions x 1 d and y 1 d are given as

φ des = arcsin m sin ψ d µ xv -cos ψ d µ yv u 1 - sin ψ d ξ xv -cos ψ d ξ yv u 1 θ des = arcsin m cos ψ d µ xv -sin ψ d µ yv u 1 cos φ des - cos ψ d ξ xv + sin ψ d ξ yv u 1 cos φ des .
It is worth mentioning that the disturbance terms ξ x and ξ y will be estimated by using the High Order Sliding Mode Observer.

Attitude Control

To develop the controller for the attitude subsystem, we express (1d)-(1f) as

φ = l I x u 2 + ξ φ θ = l I y u 3 + ξ θ ψ = 1 I z u 4 + ξ ψ ,
where

ξ φ = θ ψ Iy-Iz Ix -ld4 Ix φ, ξ θ = φ ψ Iz-Ix Iy -ld5
Iy θ and ξ ψ = φ θ Ix-Iy

Iz

-d6 Iz ψ. The attitude controller for the Quadrotor is developed using the Nonsingular Terminal Modified Super-Twisting Control, which enables the system to follow the desired orientation given by angles φ des , θ des and ψ d . To accomplish this task, the following sliding manifolds are defined as

s φ = φ 2 -φ 2 des + c φ1 φ 1 -φ 1 des + c φ2 exp -λ φ t φ 1 -φ 1 des (1-2β φ ) s θ = θ 2 -θ 2 des + c θ1 θ 1 -θ 1 des + c θ2 exp -λ θ t θ 1 -θ 1 des (1-2β θ ) s ψ = ψ 2 -ψ 2 des + c ψ1 ψ 1 -ψ 1 des + c ψ2 exp -λ ψ t ψ 1 -ψ 1 des (1-2β ψ )
,

where φ 1 = φ, φ 2 = φ, θ 1 = θ, θ 1 = θ , ψ 1 = ψ and ψ 2 = ψ, the gains c φ1 , c φ2 , c θ1 , c θ2 , c ψ1 , c ψ2 are positive, and 0

< β φ < 1, 0 < β θ < 1, 0 < β ψ < 1.
Based on the assumption that the orientation velocities are available, the following controllers are proposed for Quadrotor's attitude

u 2 = I x l -c φ2 exp -λ φ t φ 1 -φ 1 des -2β φ (1 -2β φ ) φ 2 -φ 2 des -λ φ φ 1 -φ 1 des -θ ψ I y -I z I x + φdes 2 -c φ1 φ 2 + c φ1 φdes 1 -k 1φ |s φ | 1 2 sign (s φ ) -k 2φ s φ - t 0 k 3φ sign (s φ ) dτ - t 0 k 4φ s φ dτ , (17)
u 3 = I x l -c θ2 exp -λ θ t θ 1 -θ 1 des -2β θ (1 -2β θ ) θ 2 -θ 2 des -λ θ θ 1 -θ 1 des -θ ψ I y -I z I x + θdes 2 -c θ1 θ 2 + c θ1 θdes 1 -k 1θ |s θ | 1 2 sign (s θ ) -k 2θ s θ - t 0 k 3θ sign (s θ ) dτ - t 0 k 4θ s θ dτ , (18)
u 4 = I x l -c ψ2 exp -λ ψ t ψ 1 -ψ 1 des -2β ψ (1 -2β ψ ) ψ 2 -ψ 2 des -λ ψ ψ 1 -ψ 1 des -ψ ψ I y -I z I x + ψdes 2 -c ψ1 ψ 2 + c ψ1 ψdes 1 -k 1ψ |s ψ | 1 2 sign (s ψ ) -k 2ψ s ψ - t 0 k 3ψ sign (s ψ ) dτ - t 0 k 4ψ s ψ dτ . ( 19
)

Simulation and Experimental Results

In order to show the right performance of the control strategies proposed, a set of simulation as well as experimental results carry out in outdoor environments is presented in this section.

Simulation Results

We present a set of simulation results where a square trajectory of 10 meters is tracked by the Quadrotor aircraft. The control laws ( 14),( 15),( 16),( 17)-( 19) together with observer (3) are applied to the system given by Equation ( 1). The Quadrotor parameters considered in the simulation test are given in Table I. Table II Super-Twisting Controller and for the High Order Sliding Mode Observer, respectively.

Figure 2 shows the position of the Quadrotor in the X -Y plane when a square reference is tracked. The upper graph shows the performance when disturbances are not estimated by observer (3) and compensated in the control law It can be observed the notable degradation in the performance of the vehicle generated by the external disturbances. The lower graph shows the trajectory of the UAS when the external disturbances are estimated and compensated, a significant improvement in the tracking is appreciated. Finally, Figure 4 depicts the estimation of the disturbances, which are induced to the vehicle's orientation dynamics in the form of an acceleration to simulate a more realistic situation. Notice the effectiveness of the high order sliding mode observer in the estimation of the disturbances.

Experimental Results

In this section we present the real-time results of the proposed controller for trajectory tracking of a desired path subject to external disturbances, as well as a hovering-mode flight in order to demonstrate the good performance of this control strategy when it is applied to the Quadrotor aircraft. For the implementation of the proposed algorithms at outdoors environments we used a Real Time Kinematic-Global Positioning System (RTK-GPS) that provides improved location accuracy, from the 2 -meters nominal GPS accuracy to around 3cm in the XY-plane. We used a Piksi RTK-GPS, which was designed for autonomous vehicle guidance applications, such as, formation flight and autonomous landing, GPS/GNSS research and surveying systems.

We decided to use this kind of GPS because of its effectiveness for locating an object in the 3D space (in our case, the aircraft) and because of its excellent characteristics such as the described below:

-Centimeter-accurate relative positioning (carrier phase RTK) -10Hz position/velocity/time solutions -Open-source software and board design -Low power consumption (-500mW typical) -USB and dual UART connectivity -Integrated patch antenna and external antenna input

The RTK-GPS requires two Piksi receivers to achieve centimeter-level positioning, one called Rover (mounted on top of the Quadrotor aircraft) and one acting as a reference station called the Base, as can be seen in Figure 5.

We use a Quadrotor aircraft propelled by four motors fitted with 10-inch propellers. The RTK-Piksi module is enclosed within a protective carbon fiber triangle that ensures adequate operation of the RTK-GPS. The base antenna of the RTK-GPS is placed on the top of a tripod in order to obtain a better satellite reception (see Figure 5). The rover's position is sent by the Piksi module through a serial port by using the Swift Binary Protocol (SBP) which is a protocol employed for communicating with Swift devices. The SBP is read by an embedded computer (Odroid XU4) which sends to the on-board Pixhawk the 3D position measurements (X,Y,Z) obtained from the RTK module, allowing the Quadrotor aircraft to execute the trajectory tracking or performing hover-mode flight over a specific point depending on the specific escenario. Finally, we have a laptop computer as a ground station which helps us to visualize the current state of all the variables involved in the control of the aerial vehicle when the programmed task is executed. It also allows to see if there is good reception of the GPS signal in the RTK module.

In the following sections we present the obtained results of four experiments performed In order to prove the effectiveness of the proposed controller, -Experiment 1: Hovering-mode test -Experiment 2: Line trajectory-tracking test -Experiment 3: Square trajectory-tracking test -Experiment 4: Circular trajectory-tracking test

Hovering-mode test

This experiment was conducted for performing flight in stationary mode in the cartesian X-Y plane. Initial conditions were set for the Quadrotor aircraft positions as x d = 0 m and y d = 0 m. Figure 6 and Figure 7 show the performance of the controller while flying at stationary mode over a particular point under external disturbances. Figure 6 shows the response in hover-mode subject to small wind gusts as external perturbation. As can be seen, the control successfully takes the helicopter to the desired coordinates (0, 0) and then performs flight in hover mode. Moreover, from this figure, it can be appreciated the accuracy of the hovering flight since the Quadrotor position is always inside a circle of radius of 25cm.

On the other hand, Figure 7 illustrates the behavior of the positions on the X and Y axes, respectively. It can be seen that these variables are maintained satisfactorily close to the given references (in our case the coordinates (x, y) = (0, 0)). For example, at time instants t = 2s and t = 3s on the X-Y axes (respectively) the Quadrotor suffers a disturbance due to a gust of wind, then the control compensates this disturbance in both axes returning the aerial vehicle to the given reference. We can conclude the proper operation of the RTK-GPS for this type of applications under moderated weeather conditions. Fig. 6: Best tuning for hover at (0, 0)

Line trajectory-tracking test

In this experiment, the Quadrotor aircraft performs two specific paths (lines) during the trajectory-tracking stage. Figure 8 and Figure 9 depict the position and velocity of the aerial vehicle in the X-axis in order to demonstrate the effectiveness controller to follow the desired trajectory. We used different trajectories for each experiment described by the following parametric equations

x (t) = p f 10 t t f 3 -15 t t f 4 + 6 t t f 5 ẋ (t) = p f 30 t 2 t 3 f -60 t 3 t 4 f + 30 t 4 t 5 f
, where x(t) is the position and ẋ(t) is the velocity of the Quadrotor aircraft at time t, respectively. Whereas p f is the final position (length of the path) and t f is the total route time. The following initial conditions are considered to generate the trajectory of the line: (x, ẋ) = (0, 0). As can be seen in Figure 8, the Quadrotor aircraft starts the track-line at t = 0s, and completes the trackline 20 seconds later with a total distance of 20 meters travelled. Also in the same figure you can observe the behavior of the velocity profile programmed in the aerial vehicle in order to properly perform the desired path. In both graphs of the Figure 8, it is possible to clearly observe the good performance of the controller to execute the tracking of the desired path (line) under external disturbances such as wind gusts.

On the other hand, Figure 9 shows the results of position and velocity obtained from the experiment in real-time for trajectory-tracking of a line of 30 meters. As in the previous case, we can observe an excellent performance of the proposed control strategy due to the fact that even if there are disturbances such as wind gusting, the vehicle manages to follow in an adequate way the trajectory given by the parametric equation described above As well as the given velocity profile. Notice that the position and velocity obtained from the performed experiments are very similar.

Square trajectory-tracking test

For this experiment, the aerial vehicle has the task of performing a squareshaped trajectory. Figure 10 and Figure 11 show in detail the trajectory tracking of a square reference with disturbances induced by the gust of wind present in the moment of carrying out the testing at outdoors environments.

From Figure 10 a good performance of the controller can be observed since the vehicle satisfactorily followed the desired trajectory (square) under external perturbations. Also, it can be clearly seen that in the upper-left and lower-right corners of the square-trajectory there exists small deviations due to the induction of perturbations caused by the gust of wind that occurred at that moment of the testing. While the efficiency of the control algorithm implemented in the Quadrotor aircraft for following the desired square-trajectory can be seen in Figure 10 and Figure 11, one can observer the excellent tracking performance of the proposed controller based on RTK-GPS system for X-axis and Y-axis, respectively.

In addition, the different instants in which the control signal (PWM) efficiently compensates for disturbances affecting the aerial vehicle while executing the desired square-trajectory are appreciated. A clear example of the above occurs on the Y-axis at the instant of time t = 30 s, shown in Figure 11, where the control signal reacts quickly to reject the effect of the disturbance on the Quadrotor aircraft. Finally, the control signals in roll and pitch for the square trajectory tracking experiment are presented in Figure 12.

Circular trajectory-tracking test

In Figure 13 and Figure 14 the experimental results for tracking a circular path are presented. The circular trajectory has a radius of 5 meters with center in the origin, with constant altitude (2 meters) and constant yaw angle (0 rad). The tracking performance in the X-Y plane can be appreciated in Figure 13 where it can be noted the behavior of the robust control strategies presented in this chapter. The references for X-axis and Y-axis that generate the circular path are shown in Figure 14. In this figure it can be appreciated the real-time position of the Quadrotor aerial vehicle in each one of the X and Y axes, where the system's behavior despite the external disturbances indicates the satisfactory performance of the designed controller. The control signals for roll and pitch angle generated to track the circular trajectory are depicted in Figure 15.

Conclusions

In this paper, the Nonsingular Terminal Modified Super-Twisting Controller has been proposed for the trajectory tracking problem of an UAS type Quadrotor. To estimate the traslational velocities of the system and to improve the ability of system performance robustness, a High Order Sliding Mode Observer is employed together with the Nonsingular Terminal Controller. The Controller-Observer proposed ensures that the tracking errors converge to zero in finite-time. Finally, experimental results were obtained to show the performance of the control strategy and the experimental setup.

-

  ζ = (x, y, z) T ∈ R 3 : denotes the position of the aerial vehicle's center of mass relative to the inertial frame Γ I . η = (φ, θ, ψ) T ∈ R 3 : describe the orientation of the aerial vehicle and (φ,θ,ψ) are the three Euler angles: roll, pitch and yaw, respectively. These angles are bounded as follows: roll angle by (-π/2 < φ < π/2), pitch angle by (-π/2 < θ < π/2) and yaw angle by (-π < ψ < π).

Fig. 1 :

 1 Fig. 1: Inertial and Body frame of the mini Quadrotor aircraft.

Fig. 2 :

 2 Fig. 2: Square trajectory tracking. Upper graph shows the performance without disturbance compensation. Lower graph presents the performance with disturbance compensation.

Fig. 3 :

 3 Fig. 3: Position in X-axis (upper graph) and Y-axis (bottom graph) without disturbances compensation, when the Quadrotor tracks a square reference.

Fig. 4 :

 4 Fig. 4: Estimation of the disturbances induced on the acceleration dynamics of roll and pitch angles.

Fig. 5 :

 5 Fig. 5: The RTK-GPS setup implemented on Quadrotor aircraft.

Fig. 7 :

 7 Fig. 7: Quadrotor's position in X-axis and Y-axis for Hover experiment.

Fig. 8 :

 8 Fig. 8: Position and velocity of the UAS in X-axis. Tracking a signal reference in a distance of -20 meters from the origin.

Fig. 9 :

 9 Fig. 9: Position and velocity of the UAS in X-axis. Tracking a signal reference in a distance of 30 meters from the origin.

Fig. 10 :

 10 Fig. 10: Real-time experiment of trajectory tracking of a square of 10-meter .

Fig. 11 :

 11 Fig. 11: Reference signal and real-time position of the Quadrotor in X-axis and Y-axis for the square trajectory tracking experiment

  Control signal in roll.

Fig. 12 :

 12 Fig. 12: Quadrotor's control signals for X and Y axes for Square Experiment.

Fig. 13 :Fig. 14 :

 1314 Fig. 13: Tracking performance for a circular reference trajectory in the X-Y plane

  Control signal in pitch.

Fig. 15 :

 15 Fig. 15: Quadrotor's control signals for X and Y axes for Square Experiment.

Table 1 :

 1 and Table III present the parameters employed for the Modified Quadrotor model parameters.

	Parameter	Value	Unit
	m	0.8	kg
	l	0.2	m
	g	9.81	m/s 2
	Ix = Iy	0.0075	Ns 2 /rad
	Iz	0.013	Ns 2 /rad
	k 1 = k 2 = k 3	0.01	Ns/m
	k 4 = k 5 = k 6	0.005	Ns/m

Table 2 :

 2 Modified Super-Twisting Controller Parameters.

	Orientation gains	Value	Position gains Value
	k 1φ	10.0	k 1x	1.5
	k 2φ	7.5	k 2x	3.0
	k 3φ	5.0	k 3x	1.5
	k 4φ	3.5	k 4x	1.5
	k 1θ	10.0	k 1y	2.0
	k 2θ	7.5	k 2y	4.0
	k 3θ	5.0	k 3y	1.5
	k 4θ	3.5	k 4y	1.5
	k 1ψ	10.0	k 1z	5.0
	k 2ψ	7.5	k 2z	3.5
	k 3ψ	5.0	k 3z	2.5
	k 4ψ	3.5	k 4z	1.5

Table 3 :

 3 High Order Sliding Mode Observer parameters.

	Parameter	Value
	λ 1x	6.75
	λ 2x	30.00
	λ 3x	6.75
	λ 4x	1.50
	λ 1y	6.75
	λ 2y	30.00
	λ 3y	6.75
	λ 1z	5.00
	λ 2z	1.00
	λ 3z	10.00
	λ 4z	1.30

Acknowledgements A short version of this paper, entitled: Robust Trajectory Tracking for Unmanned Aircraft Systems Using High Order Sliding Mode Controllers-Observers, was presented in ICUAS 2017. This work was partially supported by the Mexican Secretariat of Public Education (SEP) and by the Mexican National Council for Science and Technology Project Laboratorio Nacional en Vehículos Autónomos y Exoesqueletos 281920.