Chiara Gabellieri

Marco Tognon
email: marco.tognon@laas.fr

Lucia Pallottino
email: lucia.pallottino@unipi.it

Antonio Franchi
email: antonio.franchi@laas.fr

A Study on Force-based Collaboration in Flying Swarms

published or not. The documents may come

Introduction

Cooperative transportation without explicit communication, but based only on on the indirect exchange of information through the physical interaction with the environment is a very important feature for social animals. From a scientific point of view, the problem has been addressed as Stigmery theory [START_REF] Grassé | La reconstruction du nid et les coordinations interindividuelles chezbellicositermes natalensis etcubitermes sp. la théorie de la stigmergie: Essai d'interprétation du comportement des termites constructeurs[END_REF]; later, it has been regarded as the main coordinating mechanism in groups of ants for object transportation [START_REF] Kube | Cooperative transport by ants and robots[END_REF] (see Fig. 1 on the left) and has been indeed observed and studied in [START_REF] Berman | Experimental study and modeling of group retrieval in ants as an approach to collective transport in swarm robotic systems[END_REF][START_REF] Franks | Teams in social insects: group retrieval of prey by army ants (eciton burchelli, hymenoptera: Formicidae)[END_REF][START_REF] Franks | Division of labour within teams of new world and old world army ants[END_REF][START_REF] Moffett | Cooperative food transport by an asiatic ant[END_REF].

Such skills observed in nature have inspired researchers to transfer them to swarms of simple robotic agents. In fact, avoiding explicit communication would reduce hardware and software complexity, and overcome possible communication failure issues. So far, the interest has been mainly focused on terrestrial systems [START_REF] Berman | Experimental study and modeling of group retrieval in ants as an approach to collective transport in swarm robotic systems[END_REF][START_REF] Kube | Cooperative transport by ants and robots[END_REF][START_REF] Rubenstein | Collective transport of complex objects by simple robots: theory and experiments[END_REF][START_REF] Wang | Force-amplifying n-robot transport system (force-ants) for cooperative planar manipulation without communication[END_REF], where the possibility of decentralized transportation based on physical interactions has been proved. Instead, in this work, we are interested in the communication-less aerial transportation of objects by swarms. This is particularly interesting not only from a scientific point of view, thanks to the higher complexity of the problem, but also because a simple, robust and scalable object aerial manipulation technique could meet the requirements demanded by many real applications. Aerial transportation can benefit from Fig. 1: Left: red ants in cooperative transportation. Right: the tension in spider webs influences their natural frequency [START_REF] Alam | Mechanics in naturally compliant structures[END_REF]. a larger workspace and the independence from uneven terrains. However, aerial robots, e.g., multi-rotors, though agile and low cost, are typically characterized by a limited payload. Hence, a cooperative approach is a very suitable solution. Some examples of applications can be found in industrial contexts, in agriculture, and in search and rescue missions to carry necessary equipment or first aid.

Considering the lack of results for aerial systems compared to grounded ones, our purpose is to start filling this gap. In particular, we investigate the possibility of a communication-less approach for cooperative aerial manipulation with a swarm of flying agents. Furthermore, we investigate if and how the load and the cable physically connecting the robots may play the role of an implicit communication channel exploiting the forces exchanged. To the best of our knowledge, this is the first work proposing a bio-inspired algorithm for communication-less aerial manipulation by flying swarms, going beyond the two-robot-'only' scenario, recently considered in [START_REF] Gassner | Dynamic collaboration without communication: Vision-based cable-suspended load transport with two quadrotors[END_REF][START_REF] Tagliabue | Collaborative transportation using MAVs via passive force control[END_REF][START_REF] Tognon | Aerial co-manipulation with cables: The role of internal force for equilibria, stability, and passivity[END_REF].

The algorithm that we propose exploits a leader-follower paradigm where the leader agent knows where to go and hence steers the object toward the desired position. On the other hand, the follower agents follow the leader and help to sustain the weight and manipulate the load, exploiting only the implicit information contained in the force received from the load itself. It has been observed that also groups of army ants Eciton burchelli [START_REF] Franks | Teams in social insects: group retrieval of prey by army ants (eciton burchelli, hymenoptera: Formicidae)[END_REF] and Dorylus wilverthi [START_REF] Franks | Division of labour within teams of new world and old world army ants[END_REF] adopt a distinct caste distribution in transportation groups, in the sense that groups have a significant tendency to contain only one submajor, i.e., a particular type of ant. These species of army ants have proved to be very efficient in transporting objects together. Additionally, in the same works, it has been noticed that it is usually a single ant, the submajor, that starts the motion of the object, and then the rest of the group moves accordingly. Such a behavior is actually replicated in the leader-follower paradigm proposed in this paper.

While animals usually deploy items by directly touching them (direct manipulation), in our framework we have chosen an indirect manipulation technique of the load through cables. This choice is motivated by different reasons and it has been proved to be a very effective solution for cooperative aerial transportation [10-12, 19, 23]. First of all, a cable attached to the agent center of mass allows to minimize the coupling between the rotational dynamics of the agent itself from the rest of the system dynamics. This is particularly useful for aerial agents that are underactuated -the most common case -since they need to change their attitude to be able to apply forces in any direction. Furthermore, compared to other possible decoupling gripping mechanisms, as the ones in [START_REF] Nguyen | Aerial tool operation system using quadrotors as rotating thrust generators[END_REF], cables are simpler, low cost, and in general lighter. In this paper we demonstrate through several numerical simulations that a swarm of N flying agents is capable of collaborative manipulation skills based only on implicit communication. We show that a twofold major role is played by the internal force applied to the transported object. Internal forces are forces applied at the contact points on the object that stretch or compress it without producing any movement, since they counterbalance each other. The condition of zero internal force corresponds to the case in which the agents transport the object while keeping the cables vertical and applying only a force that compensates for the gravity. Firstly, we have found that nonzero internal forces allow the swarm to univocally set the attitude of the commonly transported object, and secondly, that larger internal forces reduce the convergence time of the overall system to such unique equilibrium. This creates an interesting analogy with the role of tension in spider webs, see Fig. 1-right. The breadth of analysis covers also the thorough investigation of the leader forces depending on the swarm parameters and the analysis of the benefits of a saturated nonlinear law for the leader force in order to tradeoff compliance/safety and transportation accuracy.

The paper structure is the following: Sec. 2 illustrates the dynamic model of the system. Then, we formulate the hypotheses regarding the properties of the swarm, supported by the numerical results presented in Sec. 3. A thorough discussion follows in Sec. 4. Final conclusions and future developments are presented in Sec 5.

Model

The system is composed by a set of N flying agents attached to a cable suspended load that must be deployed to a particular configuration. In our framework, each agent interacts with the environment, hence we aim at a soft response similar to the behavior of human or of animal during everyday interaction tasks. We model the commonly transported object (the load) as a rigid body. The attitude of the load is parametrized by Euler angles yaw, pitch and roll, indicated with ψ, θ , and φ , respectively. Each agent is attached to the load by means of a cable by means of which it can transfer forces. The cables are attached to the load at the points L i , with i = 1, . . . , N, placed on the same plane, denoted by I . The object center of mass (CoM) is indicated with G.

Each cable is modeled as a linear unidirectional spring, with a dissipative term that damps its longitudinal oscillations. We assume that each flying agent is endowed with a position controller. If the latter is sufficiently precise, we can model the closed loop system as a simple double integrator. In this way we can consider each agent as an actuated point mass capable of exchanging a force with the external world.

By doing so, the proposed method can be applied to different aerial robots. If we consider multidirectional-thrust platforms capable of controlling position and orientation independently (popular in the field of aerial physical interaction) [START_REF] Brescianini | Design, modeling and control of an omni-directional aerial vehicle[END_REF][START_REF] Park | Design, modeling and control of omni-directional aerial robot[END_REF][START_REF] Rajappa | Modeling, control and design optimization for a fully-actuated hexarotor aerial vehicle with tilted propellers[END_REF][START_REF] Ryll | 6D physical interaction with a fully actuated aerial robot[END_REF], the double integrator is an exact model of the position-controlled closed loop system. In the case of underactuated unidirectional-thrust vehicles using a standard position controller, the double integrator is instead a very good approximation. Furthermore, the time-scale separation between the translational and rotational dynamics has been exploited in other works on aerial manipulation like [START_REF] Michael | Cooperative manipulation and transportation with aerial robots[END_REF][START_REF] Nguyen | Aerial tool operation system using quadrotors as rotating thrust generators[END_REF] where the robots are considered as point masses and modeled as double integrators. Denoting by the vector p Ri ∈ R 3 the position of the i-th agent, by M i = m Ai I 3 its inertia matrix (with m Ai ∈ R >0 ) and by f i = f i,x f i,y f i,z the force that the i-th cable exerts on the object (so that -f i is the force exerted by the cable on the agent), the dynamics of the i-th agent is:

M i pRi = -f i + f Ci , (1)
where

f Ci = B i ( ṗd Ri -ṗRi ) + K i (p d Ri -p Ri ) + π i (2)
is the 'control' force of each agent. To better understand how to control a multi-rotor aerial robot in such a way, the reader can refer to [START_REF] Tognon | Aerial co-manipulation with cables: The role of internal force for equilibria, stability, and passivity[END_REF]. Such control force models three simple actions. First, the agent implements a spring like action to move towards a desired position or follow a desired path (the apex d indicates the 'desired' quantities). Secondly, the agent implements a dissipative derivative term proportional to the velocity error. This action damps the oscillations induced by the spring action. Finally, there is a force bias indicated with π i = π i,x π i,y π i,z . This bias is essential to make the flying agents sustain the weight of the load. We shall show that it plays an important role also for shaping the system equilibria. Through π i it is possible to set reference internal forces (forces that do not result in a motion of the object).

The static equilibrium equation of the object subject to the gravity vector g ∈ R 3 is given by [m L g 0 1×3 ] = W f, where W ∈ R 6×3N is the grasp matrix, that maps the forces at the contact points to a wrench applied at the object center of mass and f ∈ R 3N×1 collects the equilibrium cable forces. Resolving the equation for f we obtain:

f = W † [m L g 0 1×3 ] + t, (3)
where † indicates a right (pseudo)inverse, and t ∈ null(W ) is the internal force, which neither influences the object dynamics nor balances any external wrench.

Generally, in leader follower approaches it may occur that only one agent, i.e., the leader, is aware of the desired trajectory (p d Ri , ṗd Ri ). On the other hand, slave agents tend to stay where they are if no external action intervenes. We model such agents by setting K i = 0 N m -1 and ṗd Ri = 0 m s -1 in (2). When the leader starts moving, the followers will perceive a modification of the environment through a change in their cable force. This, in turn, will make the followers move toward the leader agent trying to bring back the force to the initial equilibrium value. It is worthy to note that, once each agent stops, possible load oscillations are damped thanks to the dissipative action modeled in (2).

Numerical Study

For a system of only two aerial agents and a beam/like load, in [START_REF] Tognon | Aerial co-manipulation with cables: The role of internal force for equilibria, stability, and passivity[END_REF] we formally proved, using a Lyapunov-based approach, the stability or instability of all the possible equilibria and the passivity of the overall controlled system. However, it is not trivial to extend those theoretical results to N > 2 and to a more generally shaped object. In particular, it is not straightforward to solve the so called equilibria inverse problem in [START_REF] Tognon | Aerial co-manipulation with cables: The role of internal force for equilibria, stability, and passivity[END_REF], namely

S a t = 0 [N] G ∈ I on a circle none S b t = 0 [N] G / ∈ I on a circle none S c t = 0 [N] G / ∈ I on a circle none Parametric uncertainty S d t = 0 [N] G / ∈ I random m L , N (only bounds)
to find all the possible positions of the agents and forces in the cables for each stable pose of the load. The authors in [START_REF] Abbasnejad | Direct geometrico-static problem of underconstrained cabledriven parallel robots with n cables[END_REF] showed that such problem is very complex to solve even if just less than six cables of assigned length are used. Analytical answers to the problem are therefore difficult to reach. However, a numerical study of whether some of the properties discovered in [START_REF] Tognon | Aerial co-manipulation with cables: The role of internal force for equilibria, stability, and passivity[END_REF] for the two-agent system apply also to the swarm case, where the number of follower agents may be arbitrary large and the object not only a bar, is equally interesting. Thus, in this section we extrapolate some conjectures on the expected behavior of the swarm system and validate these hypotheses through a wide numerical study. Table 1 contains the description of the simulated conditions. The load has mass m L = 5 kg and inertial matrix equal to J L = I 3 kg m 2 . The leader agent, set as the agent 1, is fed with a sufficiently smooth reference trajectory constituted by a 5 th order polynomial in 3D (rest to rest trajectory with zero acceleration at start and end points), which lasts 10 s and covers 2 m.

M 1 = 0.5I 3 kg , M i = 0.01I 3(N-1) kg , B 1 = 100I 3 N s m -1 , B i = 1.5I 3(N-1) N s m -1 , K 1 = 1000I 3 N m -1
, and K i = 0I 3(N-1) N m -1 for i = 1. When N> 30, we changed the apparent mass M i , and the damping B i of the followers, i.e., for i = 2, ..., N reducing both by 90%. This allows the leader agent to drag the system without applying too large forces, similarly to what happens in biological systems as described in Sec. [START_REF] Brescianini | Design, modeling and control of an omni-directional aerial vehicle[END_REF].

Scenarios S a , S b and S c refer to the cases in which the parameters of the system are perfectly known, and so is N. Based on what we demonstrated in [START_REF] Tognon | Aerial co-manipulation with cables: The role of internal force for equilibria, stability, and passivity[END_REF] for the two agents system, we decided to investigate the role of the internal force on the object for the equilibria of the system. We considered the case in which the points L i lie on a circle centered around the projection of G on I (see Figs. 2(a), 2(b), and 2(d)). We ran different simulations where N is a random number between 2 and 100, and the initial pose of the object has a random value of yaw between 0 rad and π/4 rad.

In the scenarios S a and S b the force bias in ( 2) is set to sustain the weight of the load without internal forces, namely π i = [0 0 m L g/N] ∀i = 1, . . . , N, where g is the intensity of g. This implies vertical cables. Once the leader agent stops, the system converges to an equilibrium in which the final attitude of the transported object is not univocally determined but depends on the trajectory resulting from the initial condition and the leader desired trajectory. In particular, the results in the first line of plots of Fig. 3 shows a completely arbitrary final orientation of the load for the scenario S a , and the second line of Fig. 3 an arbitrary final yaw for the scenario S b .

In the scenario S c we set the force bias π i in (2) so that the system reaches an equilibrium with the cables forces applying non zero horizontal forces. More in detail, indicating with li the vector that connects the projection of G on I to the position Preprint version, 11th Swarm Intelligence ANTS 2018, Budapest, Hungary --page 5/13 ). The object is depicted as a grey surface, and the light spot on it coincide with its center of mass G. In case G / ∈ I , the light spot is instead the projection of G on I . The cables are the black lines, while the circles represent the robots (the darker one is the leader). of point L i at the final attitude, we set [π i,x π i,y ] = f int li /|| li ||, where f int is thus the intensity of each agent's planar force bias. In other words, the desired horizontal forces in the cables are oriented radially and outward the object in the final configuration, similarly to what is depicted in Fig. 2(b). Since the object reaches an equilibrium, and the external wrench on the object at the equilibrium is only the vertical force due to the gravity, as in (3), such non vertical components of the cable forces, which do not cause any motion of the object and do not compensate any external wrench, generate an internal force t = 0. In this way, the object reaches always the same attitude: zero pitch, zero roll, and the same yaw. The third line of Fig 3 shows the results of such simulations. Notice that the same results have been obtained even when G ∈ I . Fig. 4 shows the evolution of the forces that the agents apply to the load in four different cases belonging to scenario S c .

We also decided to see whether the value of f int influences the speed of convergence with which the swarm stabilizes the object in scenario S c . Fig. 5 shows the evolution of the load attitude when transported by a group of five agents for different values of f int .

Preprint version, 11th Swarm Intelligence ANTS 2018, Budapest, Hungary --page 6/13 We present now the results concerning scenario S d , where neither the exact real mass of the load, m L , nor the exact total number of agents, N, are known by the agents. Only some upper bounds are given, indicated with m and N respectively. Therefore, we introduced a particular choice of π i .

π i =   0 0 f i,z +U z i   where U z i =      -K Z if -| f i,z | < -f max K Z if -| f i,z | > -f min 0 otherwise . ( 4
)
Namely, we are defining a dead-zone in the sensed vertical force such that if the robots perceive a vertical force in the cable that is inside a certain range defined by f max and f min , they ignore it; otherwise, they apply an upward or downward force trying to restore a vertical force inside the predefined range. K Z is a constant that determines the respon-Preprint version, 11th Swarm Intelligence ANTS 2018, Budapest, Hungary --page 7/13

(a) M i = 0.001I 3(N-1) kg, Fig. 7: The three components of p R1 and p d R1 (dotted lines of the same color) and the intensity of the cable forces. By modeling f C1 as in (5) it is possible, with the same parameters, to have the leader follow its desired trajectory under normal conditions (7(a)), and to limit the forces exerted in case the swarm does not follow it (7(b)). siveness of the robots in trying to maintain the vertical force inside

B i = 0.15I 3(N-1) N s m -1 (b) M i = 0.01I 3(N-1) kg, B i = 0.15I 3(N-1) N s m -1 (c) M i = 0.001I 3(N-1) kg, B i = 1.5I 3(N-1) N s m -1 (d) M i = 0.01I 3(N-1) kg, B i = 1.5I 3(N-1) N s m -1
f int = 2 N f int = 4 N f int = 6 N f int = 0.5 N 12 -10 -8 -6 -4 -2 0 2
[-f max , -f min ]. We choose f max = (m L + m)g N- N and f min = (m L -m)g N+ N .
In this way, we guarantee that the overall force exerted by the agents will be enough to sustain the object weight. In particular, we are not choosing a precise reference force distribution. The actual final force in each cable is induced by the choice of the bias (which does not depend on the exact values of N and m L ) and is not the same for all the agents. Such an implementation allows the robots to successfully cope with a variation of the parameters of the swarm and of the object. Simulations results for this scenario, see Fig. 6, show that the system stabilizes after the transportation. However, due to the uncertain conditions, the final pose of the object cannot be known a priori. The number of agents varies randomly between 7 and 50, m = 0.5 kg, and N = 5.

We conclude this section proposing and testing an alternative version of (2) for the leader agent. Thanks to the model ( 1), (2), the leader agent does not blindly follow the desired trajectory but is aware of the outer world. In the choice of the parameters of (2) one has to face a clear trade-off between compliance (i.e., energy consumption and safety) and tracking error (i.e., performance). A behavior that is too compliant can compromise the reach of the final position, a behavior that is too stiff can require a large amount of force from the leader side and thus a lot of energy consumption and an increased risk of ruptures. Due to space limitations we omit the performed numerical results showing such intuitive trade-off. A possible solution to better deal with such trade-off is to introduce a nonlinear saturation model like, e.g., the following one:

f C1 = B 1 tan -1 ( ṗd R1 -ṗR1 ) + K i tan -1 (p d R1 -p R1 ) + π i . (5)
The relative results are shown in Fig. 7.

Discussion

We first simulated symmetric cables attachment points, lying on a circle around the object center of mass (or its projection on I ). Applying a non-zero internal force that stretches the object the agents are capable of controlling univocally the final orientation of the object without the need of communication. This result actually reflects what we had discovered and formally proven for the two-agent system and a beam-like load in [START_REF] Tognon | Aerial co-manipulation with cables: The role of internal force for equilibria, stability, and passivity[END_REF]. It is possible to assume that also biological agents, transporting an object, may be able to pull the edge of the object towards themselves, applying a force orthogonal to the edge, and thus generating a resulting internal force altogether. For example, applying a constant horizontal force bias that equilibrates the object at the initial configuration, the swarm may be capable of transporting the load and repositioning it with the same attitude without the need of explicit communication. This suggests that the local sensing allows a sort of collective memory, as if the swarm of agents remembered the initial orientation of the object and were capable of recreating it. We were also interested in understanding the role played by the internal forces in the convergence rate. In fact, observing some biological structures, it seems reasonable to assume that a higher internal force in the object might lead to a faster convergence of the system to the equilibrium. For example, in [START_REF] Alam | Mechanics in naturally compliant structures[END_REF] the authors observe the dynamic response of different spider webs. It came out that an initial tension in the net changes the natural frequency in the sense that the net is faster in recovering its equilibrium after being perturbed with respect to the case without initial tension. Similarly, we found that increasing the intensity of the horizontal bias force in f C , and hence in turn increasing the internal force in the object, leads to faster a convergence rate, as shown in Fig. 5.

The apparent mass of the leader agent has been set to a larger value than the one of the followers. Actually, it has been observed in groups of different army ants (Eciton burchelli, Dorylus wilverthi) that the front ant in prey retrieval is larger than the others, and it is characterized by a larger dry weight. Additionally, in our simulations, during the dynamic part, the leader agent has to apply greater forces to the object than the ones in the group. An example is shown in Fig. 4. However, this seems to be true also in super efficient prey retrieval groups of army ants [START_REF] Franks | Teams in social insects: group retrieval of prey by army ants (eciton burchelli, hymenoptera: Formicidae)[END_REF], where a single submajor initiates the motion of the item by itself, and only after the rest of the group helps in the transportation.

The damping parameter of the followers B i and their masses M i have been tuned depending on the number of robots. This has been done to reduce the force that the leader agent needs to apply in order to let each follower start moving. This force increases with the increasing of the damping parameters of the followers that are related to energy dissipation of the overall system. The master force in case of a small N with a large damping is comparable to the force that the master applies with a large N and small damping, a behavior which is reported in Fig. 4(a). To equilibrate as much as possible the force exerted by the leader for different values of N, one solution is to modify the damping parameter of the followers based on N. Another way to modify the effect of the follower inertia on the leader it is to modify their mass. Compare Figs. 4(d) and 4(b) to see the effect of the follower inertia on the master applied force, and Figs. 4(d) and 4(b) to see the effect of the damping parameter. Finally, comparing Figs. 4(d) and 4(a) the benefits of decreasing both the followers mass and damping emerges. It is not immediately clear which is the possible biological meaning of tuning the mass of the agents though. However, it is interesting to highlight that Eciton burchelli ants, very efficient in transportation tasks, tend to maintain a constant ratio between the total weight of porters and the weight of the carried load [START_REF] Franks | Teams in social insects: group retrieval of prey by army ants (eciton burchelli, hymenoptera: Formicidae)[END_REF]. f Concerning the leader spring parameter, it is not trivial to find the right balance in order to have a leader agent capable of both following the desired trajectory not being too much perturbed by small forces (for instance the ones caused by the followers initial resistance to the motion) and being sensitive to large forces on the cable indicating that something is not going as expected. However, that is exactly what an intelligent biological system would do. We were able to mimic this behavior by using a nonlinear control action on the leader agent, see [START_REF] Franks | Teams in social insects: group retrieval of prey by army ants (eciton burchelli, hymenoptera: Formicidae)[END_REF]. In this way, the leader follows the desired trajectory accurately for small deviations, but it changes its behavior consistently if the followers stay still and not follow as expected (see Fig. 7).

Finally, the results in Fig. 6, which refer to scenario S d , show the potential of the algorithm also in conditions were the swarm is not completely aware of the parameters of the system. Even with a limited and imperfect knowledge the agents are capable of commonly carrying an object in a stable way.

Conclusions

This work is a simulative study on communication-less cooperative object transportation by swarms of aerial agents. The main focus concerns the important role of the internal force to make asymptotically stable certain equilibria of the system and to enhance the manipulation capabilities of the swarm. An imperfect knowledge of the system parameters has been also treated. Several parallelisms with biological examples are discussed as well. Validation of the proposed algorithm on real platforms and theoretical proofs is left as an important future work. Of course, in order to realize a practical implementation, additional aspects, such as a collision avoidance technique, are required, especially for very large number of agents. Furthermore, the relaxation of some hypothesis, e.g., on the cables attachment points or on the rigidity of the load, may represent an interesting enlargement of the application domain.

Preprint version, 11th Swarm Intelligence ANTS 2018, Budapest, Hungary --page 11/13

Fig. 2 :

 2 Fig. 2: Simulated system for the case N = 6 (2(a), 2(b), 2(c)) and N = 60 (2(d))). The object is depicted as a grey surface, and the light spot on it coincide with its center of mass G. In case G / ∈ I , the light spot is instead the projection of G on I . The cables are the black lines, while the circles represent the robots (the darker one is the leader).

Fig. 3 :

 3 Fig. 3: Object attitude during transportation. First row: scenario S a : f int = 0 N and G ∈ I . Second row: ??), scenario S b : f int = 0 N and G below I . Third row: ??), scenario S c : f int = 0.8N and G placed below I . The dotted line is the desired value. The first, second, and third columns refer to the trajectories of the yaw (ψ), pitch (θ ), and roll (φ ) angles of the object, respectively.

Fig. 4 :

 4 Fig. 4: Intensities of the vertical force f i,z and of the horizontal force f 2 i,x + f 2 i,y that the agents apply to the object during the transportation for four different parameters of the followers in scenario S c with N = 70. The dotted line refers to the leader agent, while the solid lines refer to the followers. In the subcaptions i = 2, . . . , N.

Fig. 5 : 5 .Fig. 6 :

 556 Fig. 5: Evolution of the object attitude (yaw, pitch and roll, respectively) during transportation for different values of the intensity of the internal force in scenario S c with N = 5.

Table 1 :

 1 The different simulations setups (scenarios).

	Goal of the study	Scenario Internal Forces CoM Points L i Unknown Parameters
	Internal force role

Preprint version, 11th Swarm Intelligence ANTS 2018, Budapest, Hungary --page 2/13

Preprint version, 11th Swarm Intelligence ANTS 2018, Budapest, Hungary --page 3/13

Preprint version, 11th Swarm Intelligence ANTS 2018, Budapest, Hungary --page 4/13

Preprint version, 11th Swarm Intelligence ANTS 2018, Budapest, Hungary --page 9/13

Preprint version, 11th Swarm Intelligence ANTS 2018, Budapest, Hungary --page 10/13

Acknowledgments. This research was partially supported by the ANR, Project ANR-17-CE33-0007 MuRoPhen.


Markus Ryll
email: ryll@mit.edu

Davide Bicego
email: davide.bicego@laas.fr

Antonio Franchi
email: antonio.franchi@laas.fr

A Truly Redundant Aerial Manipulator exploiting a Multi-directional Thrust Base

Keywords: Unmanned aerial robots, Redundant manipulators, Robot kinematics

published or not. The documents may come

INTRODUCTION

Nowadays, Unmanned Aerial Vehicles (UAVs) are extensively employed in different application scenarios like remote monitoring [START_REF] Merino | An unmanned aircraft system for automatic forest fire monitoring and measurement[END_REF] and aerial photography, search and rescue missions like SHERPA (2013)(2014)(2015)(2016)(2017) and to perform many other contact-less operations. In recent years, new research efforts have been made towards the accomplishment of aerial physical interaction tasks which require the UAV, equipped with one (or more) robotic manipulator(s), to get in contact with the surrounding environment. Some examples of such operations can be found in aerial manipulation and grasping [START_REF] Mellinger | Design, modeling, estimation and control for aerial grasping and manipulation[END_REF], peg-in-hole [START_REF] Ryll | 6D physical interaction with a fully actuated aerial robot[END_REF], structure assembly and decommissioning [START_REF] Staub | Towards a flying assistant paradigm: the OTHex[END_REF] and also cooperative transportation [START_REF] Mellinger | Cooperative grasping and transport using multiple quadrotors[END_REF]. Such a direction has also been fostered by European Projects like AeRoArms (2015AeRoArms ( -2019)). This new field of complex tasks yields new challenges in the mechanical structure, in the design of aerial manipulators [START_REF] Orsag | Modeling and control of MM-UAV: Mobile manipulating unmanned aerial vehicle[END_REF], in the modeling methodologies [START_REF] Yang | Dynamics and control of quadrotor with robotic manipulator[END_REF], and specially in the control. In the last years, we saw a development from UAVs equipped with simple gripper, towards highly complex n-degree of freedom (DoF) compliant manipulators [START_REF] Yüksel | Design, identification and experimental testing of a lightweight flexible-joint arm for aerial physical interaction[END_REF], [START_REF] Baizid | Experiments on behavioral coordinated control of an unmanned aerial vehicle manipulator system[END_REF]. As these complex manipulators can be in constant motion, the reaction forces can destabilize the flight and the interaction of the UAV. To circumvent this problem, an applicable control scheme is needed not only to ensure the stability of the platform, but also to allow a precise tracking of a desired trajectory of the manipulator end-effector.

To fulfill these requirements, several control schemes have been developed. An adaptive scheme has been presented in [START_REF] Antonelli | Adaptive control of armequippedquadrotors -theory and simulations[END_REF], which compensates the manipulator mass and forces due to its motion. [START_REF] Lippiello | Cartesian impedance control of a uav with a robotic arm[END_REF] follow a similar modeling approach with an underactuated quadrotor design for a Cartesian impedance control. An approach based on the Lagrangian formulation is presented in [START_REF] Yang | Dynamics and control of quadrotor with robotic manipulator[END_REF], while the controller is based on a backstepping-like end-effector tracking law. In [START_REF] Orsag | Valve turning using a dual-arm aerial manipulator[END_REF] a dual arm manipulator is used in a valve turning scenario. A sequential Newton method for unconstrained optimal control is used in [START_REF] Garimella | Towards model-predictive control for aerial pick-and-place[END_REF] for model-predictive control in pick and place tasks with an aerial manipulator and [START_REF] Kim | Operating an unknown drawer using an aerial manipulator[END_REF] presents a work focused on the task of opening an unknown drawer.

In this paper we present a novel approach -the combination of a fully actuated hexarotor platform combined with an arbitrary n DoF manipulator. This new approach comes with the drawback of a more complex system compared to standard planar multirotor systems and higher internal forces but is beneficial in several points. Roll and pitch angles of the aerial platform can be controlled independently from the platform position (e.g., for more advanced obstacle avoidance), [START_REF] Franchi | Full-pose tracking control for aerial robotic systems with laterally-bounded input force[END_REF]. A planar rotor system with in total the same number of DoFs, would need a manipulator with two additional DoFs. Standard aerial manipulating platforms are underactuated in the dynamics of their center of mass, e.g. [START_REF] Yang | Dynamics and control of quadrotor with robotic manipulator[END_REF]. Thus, these systems are only able track trajectories that are smooth in C 4 for the lateral position. For the proposed system a smoothness C 2 is sufficient.

The contributions of our solution are, first of all, a novel combination of fully actuated aerial platform and redundant manipulator. At the best of our knowledge, this arrangement has never been thoroughly investigated so far. Secondly, we present an entire modeled, combining the aerial vehicle and the manipulator, making advantage of the well-known Denavit-Hartenberg parametrization. Thirdly, we propose a novel controller based on full feedback linearization of the dynamics and the exploitation of the system redundancy.

The structure of the paper is as follows. We will derive the dynamical model in Sec. 2. Afterwards, in Sec. 3 we will take advantage of the model to develop a controller based on feedback linearization and present a redundancy Fig. 1. Exemplary schematics of the considered aerial platform and the n-degree robotic arm. For clarity, the platform and the arm are drawn detached and the arm has been cut between link 2 and n -1. The following properties are presented: I li -Moment of inertia of link i, m li -mass of link i, a i -distance between joint i and i + 1, l i -distance between joint i and CoG of link i, q i -generalized coordinate of joint i.

exploitation approach, followed by Sec. 4 where we present a use-case scenario with a 2 and a 4 DoF manipulator.

Based on this scenario, we test the validity of our control framework in practical case simulations in Sec. 5. Finally, we conclude the paper with an outline of our results and a hint of future works.

MODELING

In this work we consider an aerial manipulator composed by a fully-actuated aerial vehicle (in the following denoted as the 'aerial vehicle') and a manipulator mounted onboard of it. The aerial vehicle that we consider is a hexarotor with tilted propellers proposed in [START_REF] Rajappa | Modeling, control and design optimization for a fully-actuated hexarotor aerial vehicle with tilted propellers[END_REF]. Fig. 1 provides a scheme of the aerial manipulator as well as the definition of the main symbols adopted in the paper, which are formally introduced in the following.

We start by defining a fixed world frame, F W , whose axes (unit vectors) are {x W , y W , z W } and origin is O W . Then we define two moving frames. The frame The aerial vehicle configuration in 3D space is fully described by p r and by the rotation matrix R r ∈ SO(3) (where

F R : O R - {x R , y R , z R }
SO(3) = {A ∈ R 3×3 |AA T = I, det(A) = 1}), representing the rotation of F R w.r.t. F W . The angular velocity of F R with respect to F W , expressed in F W , is denoted with ω r ∈ R 3 .
The manipulator has n degrees of freedom and we denote with

q m = [q m 1 • • • q m n ] T ∈ R n its generalized coordi-
nates, where q m i represents the i-th joint variable, with i = 1 . . . n. The end-effector configuration in 3D space is fully described by p e and by the matrix R e ∈ SO(3), representing the rotation of F E w.r.t. F W . The angular velocity of F E with respect to F W , expressed in F W , is denoted with ω e ∈ R 3 .

In order to use the Lagrangian approach we write the configuration of the aerial vehicle with a minimal set of generalized coordinates as q r = [q r T 1 q r T 2 ] T ∈ R 6 , where q r 1 = p r , and q r 2 = η = [φ θ ψ] T is a minimal parameterization of R r , like, e.g., the roll, pitch, and yaw angles. Finally, we denote with T η the matrix such that ω r = T η η. We will show in Sec. 3.2 that using minimal representation is not restrictive in our case as the attitude of the aerial vehicle will not exceed near hovering condition. The complete configuration of the aerial manipulator is then described by q = [q T r q T m ] T ∈ R 6+n . The dynamical model of the aerial manipulator is easily derived applying the Euler-Lagrange equation to the kinetic and potential energy of the system, thus obtaining

M(q)q + c(q, q) + d( q) + g(q) = τ r τ e , (1)
where M(q) ∈ R (6+n)×( 6+n) is the positive definite inertia matrix of the whole aerial manipulator, c(q, q) ∈ R 6+n contains all the centrifugal and Coriolis terms, d( q) ∈ R 6+n contains the friction terms, g(q) contains all the gravity terms, and τ r ∈ R 6 and τ e ∈ R n are the control wrench applied to the aerial vehicle and the control torques applied to each joint of the manipulator, respectively. This approach is thoroughly presented in [START_REF] Lippiello | Cartesian impedance control of a uav with a robotic arm[END_REF].

The wrench τ r can be decomposed in a 3D force and a 3D moment [START_REF] Rajappa | Modeling, control and design optimization for a fully-actuated hexarotor aerial vehicle with tilted propellers[END_REF], i.e.,

τ r = τ r1 τ r2 = τ r1 τ thrust + τ drag , (2)
where the single components are computed as

τ r1 = k f R r 6 i=1 R r pi,α e 3 w i , (3)
where k f is the propeller force coefficient, R r pi,α is the constant rotation matrix expressing the orientation of the i-th propeller with respect to F R , e 3 = [0 0 1] T , and w i is the square of the spinning velocity of the i-th propeller;

τ thrust = k f R r 6 i=1 p R i × R r pi,α e 3 w i , (4)
where p R i ∈ R 3 is the constant position of the center of the i-th propeller expressed in F R ; and

τ drag = k d R r 6 i=1 R r pi,α e 3 (-1) i w i , (5)
where k d is the propeller drag coefficient. The factor (-1) i in ( 5) is used to take into account the contrariwise spinning velocity of every second propeller. The subscript in R r pi,α denotes the dependency on the angle α ∈ (0, π/2], which represents the amount of tilting of the i-th propeller about the direction of p R i . Replacing (3), (4), and (5) in (2) one can compactly write

τ r = G w (η, α)w, (6)
where w = [ ŵ1 ŵ2 ŵ3 ŵ4 ŵ5 ŵ6 ] T , where ŵi = sign(w i )w 2 i and G w (η, α) ∈ R 6×6 . Equation ( 6) can be used to replace τ r in (1) in order to explicitly show the dependence of the dynamics on w.

The matrix G w (η, α) plays a crucial role in the actuation capabilities of the aerial vehicle. If α is 0 then all the propellers would be coplanar and the aerial vehicle would be underactuated, with G w being non invertible and degrading the system to an ordinary coplanar hexarotor. If 0 < α < π 2 then G w results invertible for any η and the aerial vehicle is fully actuated. However if the aerial vehicle is too much tilted (i.e., when φ and θ are too large) the inversion w = G -1 w (η, α)τ r might return non-positive values in one or more entries of w thus resulting in an unfeasible command. It is then important to keep φ and θ within the operating condition that guarantee the fullactuation with positive propeller spinning velocities (see Sec. 3.2).

CONTROL

The In order to achieve such an objective let us start from replacing τ r from ( 6) in (1), thus obtaining

M(q)q + c(q, q) + d( q) + g(q) = K(η, α) w τ e , (7)
where

K(η, α) = G w (η, α) 0 0 I n ∈ R 6+n×6+n .
We then first apply the following inner control loop

w τ e = K -1 Mu q + c + d + g , (8)
where we omitted the dependency from q and q for brevity. Control (8) brings the dynamics of q in the fully linearized and decoupled form q = u q , where u q is an additional virtual input. Notice that this result is not possible for standard aerial manipulators using coplanar multi-rotors as aerial platform.

As customary, see, e.g., [START_REF] Siciliano | Robotics: Modelling, Planning and Control[END_REF] the differential kinematics of the aerial manipulator can be expressed as

v e = ṗe ω e = ṗr ω r + J m (q m ) qm = J(q) q, ( 9
)
where J m (q) and J(q) are the geometric Jacobians of the manipulator and of the complete aerial manipulator, respectively. By differentiating (9) we obtain the dynamics of the end-effector configuration ve = J(q, q) q + J(q)q = J(q, q) q + J(q)u q (10) where we applied the inner control loop (8).

Using the fact that J has always full rank for any θ = π(k+ 1 2 ) with k ∈ Z (see Sec. 3.1) as typically done, we choose the virtual input u q as u q = J q) † (u e -J( q, q) q + I 6+n -J(q) † J(q) z (11)

with [] † being the Moore-Penrose inverse of the argument and u e and z two additional virtual inputs. Under the effect of the second inner control loop (11) the dynamics of the end-effector task becomes exactly linearized and fully decoupled ve = u e = u e1 u e2 , regardless of the choice of the second input z, whose use will be defined later. Now we are able to solve the desired control problem stated at the beginning of this section, i.e., the exact asymptotic tracking of a desired trajectory For the end-effector position it is sufficient to use 12) in ( 11) to obtain exponential and decoupled position error convergence where K P1 , K P2 define a Hurwitz polynomial. For the orientation of the end-effector it is instead sufficient to choose u e2 = ωd e + K ω1 (ω d e -ω e ) + K ω2 e R (13) with the orientation error defined as

u e1 = pd + K P1 ( ṗd e -ṗe ) + K P2 (p d e -p e ) (
e R = 1 2 [R T e R d e -R dT e R e ] ∨ (14)
where [] ∨ represents the inverse map from so(3) to R 3 [START_REF] Lee | Geometric tracking control of a quadrotor uav on se(3)[END_REF] and K ω1 , K ω2 define a Hurwitz polynomial as well.

The task Jacobian matrix J possesses an n-dimensional null space onto which the additional control input z is projected, see ( 11). The use of z ∈ R 6+n to fulfill additional tasks besides the tracking of a desired trajectory is discussed in Section 3.2.

Full-rankness of J

In order to demonstrate that J is always full rank, it is sufficient to prove that it exists a [6 × 6] submatrix of J which is invertible. the determinant of A can be found to be det(A) = cos(θ). ( 16)

It follows that rank(J) = 6 for any θ = π(k+ 1 2 ) with k ∈ Z. In fact this corresponds to the singularity of the roll, pitch, yaw representation. Notice that this holds independently from the kind of manipulator mounted onboard of the aerial vehicle.

Redundancy Exploitation

In this section we discuss the exploitation of the ndimensional redundancy of the system by designing a control law for the additional input z in (11) using the projected gradient method [START_REF] De Luca | Robot redundancy resolution at the acceleration level[END_REF]. As optimization criteria to be fulfilled by the redundancy exploitation we opt for (1) Horizontal aerial vehicle orientation: Notice that if q = 0, the most efficient configuration in terms of propeller spinning velocities is with θ = 0 and φ = 0. Therefore the optimization based on H 1 shall automatically drive the system to a state close to θ = φ = 0. By carefully selecting maximum angles for θ and φ we can achieve a bounded behavior of the control output w in (8). (2) Obstacle avoidance: The aerial manipulator is meant to interact with the environment using its endeffector. Consequently any other collision must be avoided. Especially any collision with the propellers would usually result in a highly critical situation. We therefore seek to maximize the distance between the aerial vehicle and surrounding obstacles. (3) Limitation of arm joints: To prevent manipulator self-collisions and collisions between the manipulator and the aerial platform a cost function H 3 is utilized.

Minimization of the three objectives concordant with the output tracking task can be achieved by choosing We are now in the state to compose the individual cost functions. We choose H 1 as

z = ∇ q H 1 + ∇ q H 2 + ∇ q H 3 (
H 1 (q) = 1 2n k 1 5 i=4 tan 2 (γq i ) ( 18
)
with γ = π 2k φ,θ and k 1 as suitable scalar gain and k φ,θ as maximum tilting angle. H 2 is selected as

H 2 (q) = k 2 f ( min oo∈Oo,o R ∈O R o c (q) -o R ) ( 19
)
with o being an obstacle point and O the set of obstacle points. f () is a function mapping the obstacle distance to a force and a torque acting on the aerial robot's CoM and k 2 a suitable scalar gain.

And finally the third cost function is similarly defined as

H 3 (q) = 1 2n k 3 n-1 i=7 tan 2 (γq i ) ( 20
)
with γ = π 2karm and k 3 as suitable scalar gain and k arm as maximum arm angle.

USE-CASE SCENARIO

In order to concretely present the general system described so far, let us now focus on two use case scenarios, first (1) a system with a 2 DoF manipulator depicted in Fig. 2-(a) and second (2) a system with a 4 DoF manipulator depicted in Fig. 2-(b). In both scenarios, all the links of the manipulator lay in the {x R , y R } plane of F R . The first example realizes two elbow joints, while the second example extends the first by a third elbow-joint and a final wrist. The Denavit-Hartenberg parameters of the manipulators can be found in Tab. 1. These configurations have been purposefully chosen driven by the idea that the aerial vehicle can primarily be in charge of the desired endeffector position and its yaw orientation, while the remaining orientation quantities can be primarily controlled by the manipulator. The total number of degrees of freedom of the aerial manipulator, with n + 6 = 8 and n + 6 = 10, still leaves extra maneuvering room for the exploitation of the redundancy.

Added Uncertainties and Simplifications

We decided to test the robustness of the controller against a model with realistic uncertainties. Based on previous experiences, we firstly added to all model parameters a Gaussian distributed uncertainty with different standard deviations depending on the measurability of the property (mass properties: σ 1 = 0.01, length properties: σ 2 = 0.02, inertia properties: σ 3 = 0.05). The controller is fed with the nominal property values while the model contains the altered real values. Secondly, we considered communication delays, control quantization effects and actuation noise on the propellers and joint actuation. Thirdly we overlaid sensor measurements (e.g. joint position, aerial robot position and their derivatives) with a realistic noise profile. Fourthly the dynamics of the brushless controllers and the propellers are modeled with a first order transfer function whereas the controller is based on the presented model. Finally, for simplicity we neglected the effects due to motors and gear boxes driving the manipulator joints.

NUMERICAL VALIDATION

In this section we present the simulation results to validate the presented control framework. Firstly, we present results of two trajectory tracking experiments (I & II) with configuration (1) demonstrating the general tracking capabilities of the system and the benefits of the optimization. Secondly, we will show the robustness of the method with a Monte Carlo experiment (III) exploiting configuration (2). The interested reader is referred to a video (https://youtu.be/9DBKfToHWGM) showing the presented experiments in detail.

General Tracking and Optimization Benefits

In experiment (I) the end-effector of the manipulator shall track a desired trajectory describing a 360 deg rotation about the end-effector axis y e (see Fig. 2-(a) and Fig. 4-(a)) while keeping the position constant. A time-lapse picture of the trajectory is depicted in Fig. 3. The end-effector tracking error remains marginal during the maneuver with a maximum position error max(|e e |) < 0.02 m (see Fig. 4-(b &c)). Thanks to the optimization in Sec. 3.2 the Table 1. Denavit-Hartenberg parameters of the considered manipulator (compare with Fig. 2). In experiment (II) the end-effector of the manipulator shall track a translational trajectory p d e parallel to the world frame axis x W . The desired orientation R d e of the endeffector remains constant. Given the initial state of the robot the aerial vehicle would collide with an obstacle that lies in the way of the trajectory (see Fig. 5). To track the trajectory safely, the aerial robot has to exploit its redundancy and change the position of the aerial robot with respect to the end-effector. It is obvious that the aerial vehicle performs a complex trajectory in order to avoid the obstacle. Thanks to the redundancy exploitation in (19) this capability is achieved.

Link a 1 α i d i υ i 1 a 1 0 0 υ 1 2 a 2 0 0 υ 2 3 0 π 2 0 π 2 + υ 3 4 0 0 a 3 + a 4 υ 4

Monte Carlo Simulation

To further test the stability of our system and to find its limits we conducted a Monte Carlo simulation on configuration (2). In the simulation we follow an eight-shaped trajectory while increasing the uncertainties of the model (see Sec. 4.1) by multiplying the expected uncertainty with a constant parameter k. In every condition of the Monte Carlo simulation, we performed the individual simulations 100 times with randomized modeling error distribution. In total we executed seven conditions with an increasing standard deviation of a truncated Gaussian distribution (truncated at two standard deviations) as follows k[σ 1 σ 2 σ 3 ] with k=[0 1 2 3 5 10 15]. In Fig. 7 the results are presented. It becomes clear that with an increasing uncertainty the performance degrades. Additionally, with a very high uncertainty, trials become unstable. This is expected as the controller cannot deal with extreme modeling errors. Anyway, even with a ten times higher uncertainty than expected no trials are unstable.

CONCLUSION

In this paper we presented a novel aerial manipulator based on a fully actuated aerial vehicle and an n-degree of freedom manipulator. We derived the dynamical model as well as the Denavit-Hartenberg parameters for the system. Then, we developed a controller based on feedback linearization aiming at trajectory tracking with the endeffector. Using an optimization scheme, we exploited the redundancy by optimizing the rotor spinning velocities (by keeping the platform horizontal), the aerial vehicle position and the manipulator joints positions. In an intense simulation and evaluation section we presented the capabilities of the aerial manipulator and as well showed the robustness of the system and its limits in a Monte Carlo simulation.

Fig. 2 .

 2 Fig. 2. Schematic side view of the two considered manipulators. (a) The two joints rotate about the same axis in an elbow manner. (b) The manipulator (a) is extended by a third elbow joint and a last wrist like joint.

  17)where ∇ [] is the gradient with respect to the variables [].

Fig. 3 .Fig. 4 .

 34 Fig. 3. Time-lapse picture of the aerial manipulator during experiment (I). The end-effector performs a full rotation while the position is constant. Thanks to the optimization the aerial vehicle remains almost horizontal. For clarity: the single time instances are presented next to each other, while in reality, the small red dot in all pictures overlaps.

Fig. 5 .Fig. 6 .

 56 Fig. 5. Time-laps picture of the aerial robot avoiding an obstacle. The end-effector is tracking the translational trajectory (blue line) while the aerial vehicle automatically dodges the obstacle.

Fig. 6

 6 Fig.6presents the results of experiment (II). The first plot of Fig.6presents the desired trajectory of experiment (II). While following the desired trajectory the translational error and orientation error remain very small (plot two and three). The orientation of the aerial vehicle remains as well almost horizontal, with a maximum pitch angle of θ = 6 deg. The fifth plot shows the actual spinning velocity of the six propeller. It is clear that all propeller rates stay positive and close to the nominal horizontal spinning velocity. Finally, the last plot depicts the position of the aerial vehicle and the end-effector in the x-z plane.

Fig. 7 .

 7 Fig. 7. Upper figure: Number of trials where instability was reached out 100 trials. With increasing uncertainties the number of unstable trials increases. Lower figure: Box plot of the Monte Carlo simulation showing the position error e p in the different conditions k.

  is rigidly attached to the aerial vehicle and its origin O R coincides with the Center of Mass (CoM) of the aerial vehicle; whereas the frame F E : O E -{x E , y E , z E } is rigidly attached to the end-effector of the manipulator. The positions of O R and O E expressed in F W are denoted with p r ∈ R 3 and p e ∈ R 3 , respectively.

  control problem considered in this work is to let p e and R e exactly track a desired arbitrary trajectory

	expressed as (p d e (t), R d e (t)) ∈ R 3 × SO(3), while exploiting the redundancy to possibly optimize additional require-
	ments. We assume that (p d e (t), R d e (t)) ∈ C3 and that
	( ṗd

e (t), ω d e (t)), and ( ṗd e (t), ωd e (t)) are also provided, as customary.

Preprint version,

12th SYROCO 2018, Budapest, Hungary (2018)

ACKNOWLEDGEMENTS

This work has been funded by the European Union's Horizon 2020 research and innovation programme under grant agreement No 644271 AEROARMS.


Michele Furci
email: michele.furci@laas.fr

Davide Bicego
email: davide.bicego@laas.fr

Antonio Franchi
email: antonio.franchi@laas.fr

Design and Input Allocation for Robots with Saturated Inputs via Genetic Algorithms

Keywords: input allocation, fully actuated/vectored-thrust aerial vehicles, genetic algorithms

In this paper we consider fully-actuated and redundantly-actuated robots, whose saturated inputs can have high bandwidth or can be slowly varying (with dynamics). The slowly varying inputs can be considered as configurations for the system. The proposed strategy allows to find the optimal actuators' configuration to optimize a cost function as the manipulability or the energy consumption. The approach allows for both a static design, which can include actuators' parameters such as position, orientation, saturations, numbers of actuators, and for a dynamic design, where the configurations can be controlled by an input of the system. A generalized solution to the optimal problem is proposed with the use of genetic algorithms. The results are validated in two simulation scenarios: a reconfiguration of the actuators orientation of an redundantly-actuated planar robot for trajectory tracking and the design optimization of the orientation of the motors in a generalized hexa-rotor with arbitrary propeller orientation.

INTRODUCTION

In control theory, an redundantly-actuated system consists of a system with more inputs than the dimension of the configuration space or the D.o.F. A common hierarchical control architecture [START_REF] Johansen | Control allocation-a survey[END_REF]) for such systems consists in having a high-level control that generates a virtual input which stabilizes the system. The control allocation strategy then commands the multiple actuators to obtain the desired virtual input, trying to exploit the higher dimension space of the inputs to optimize some cost functions (for energy, saturations, fault tolerance, etc). Very often in physical-mechanical system we might have actuators with very high bandwidth which can be considered as almost instantaneous, and others with low bandwidth evolving with dynamics [START_REF] Passenbrunner | Optimality-based dynamic allocation with nonlinear first-order redundant actuators[END_REF]). In this paper we consider such systems, where some slow varying parameters (we call configurations) have dynamics, while some high bandwidth inputs (we call instantaneous inputs) can be obtained instantaneously. Moreover, the inputs are subject to saturation constraints. We focus on the optimization of the configurations, while the instantaneous inputs are found using standard techniques as the constrained pseudo-inverse. In particular, we address two problems: the control problem where the configurations can actually be controlled and modified by some inputs, and the design problem where the configurations are static and the goal is to find an optimal static configuration at design time to optimize a cost function. Possible applications include, but are not limited to: ships and vessels [START_REF] Lindegaard | Fuel-efficient rudder and propeller control allocation for marine craft: Experiments with a model ship[END_REF]), aircrafts [START_REF] Oppenheimer | Control allocation for over-actuated systems[END_REF], [START_REF] Boskovic | Control allocation in overactuated aircraft under position and rate limiting[END_REF]), wheeled mobile robots and multirotors unmanned aerial vehicles [START_REF] Ryll | Modeling and control of a quadrotor uav with tilting propellers[END_REF][START_REF] Ryll | Modeling and control of FAST-Hex: a fully-actuated by synchronizedtilting hexarotor[END_REF], [START_REF] Brescianini | Design, modeling and control of an omni-directional aerial vehicle[END_REF], [START_REF] Park | Design, modeling and control of omni-directional aerial robot[END_REF]).

Since the resulting problem might be in high dimension, non-linear, with non-differentiable and non-closed form cost functions, we propose a generalized solution with the use of genetic algorithms [START_REF] Srinivas | Genetic algorithms: a survey[END_REF]).

We propose different cost functions associated to different tasks for the robot, in particular we describe cost functions to optimize: the energy consumption of the instantaneous inputs, the volume of attainable virtual inputs and the manipulability index, that is the biggest ball of obtainable virtual input around a desired virtual input. The contribution of this paper is on the use of Genetic Algorithm to solve complex non-analytic problem. In [START_REF] Passenbrunner | Optimality-based dynamic allocation with nonlinear first-order redundant actuators[END_REF], for example, a gradient based method is employed which is impossible to apply in our case. In [START_REF] Lindegaard | Fuel-efficient rudder and propeller control allocation for marine craft: Experiments with a model ship[END_REF], due to the analytic solution, the approach has limitations on the number of saturated actuators. Many approaches try to use Linear or Quadratic Programming techniques, but those require very specific cost functions. With the proposed approach we are able to optimize very complex cost functions as the manipulability index which to the best of our knowledge was never optimized for an redundantly-actuated system. To conclude, we simulated the proposed approach for both a control and a design problem; the first on a planar robot with actuators consisting in directable thrusters and the second on a generalized Multi Directional Thrust (MDT) hexa-rotor with arbitrary propeller direction.

Organization

The paper is organized as follows: Section 2 briefly presents the notation used in the paper. In section 3 the problem is introduced and some relevant cost functions for different tasks are outlined. In section 4 a generalized solution with the use of Genetic Algorithms is proposed and finally section 5 presents some results on the simulation for both the control and design problems.

NOTATION

A set S ⊂ R n is called polyhedron if it is the intersection of a finite set of closed halfspaces:

S = {x ∈ R n |Cx ≤ b} (1) with C ∈ R q×n , b ∈ R q ,
with q the number of constraints defining the hyper-planes. This is the so called H-representation of polyhedron.

A bounded polyhedron is called a polytope.

The operator S × (a) indicates the skew symmetric matrix such that S × (a)b = a×b with × the cross product operator and a, b two vectors of same dimension.

PROBLEM FORMULATION

Consider a dynamical system: ẋ = f (x, w)

(2) with x ∈ R m is the state of the system and where w ∈ R n are some virtual inputs. The virtual inputs depends linearly from the instantaneous inputs via an allocation matrix, and in particular:

w(t) = A(θ(t))u p (t) (3a) θ(t) = f θ (θ, u θ (t)) (3b) with A ∈ R n×p the allocation matrix, u p ∈ R p a saturated vector of instantaneous inputs, p ≥ n, θ(t) ∈ R n θ some time-variant parameters (the configurations), u θ ∈ R p θ a vector of saturated inputs for the parameters θ, f θ (•) : R p θ → R n θ a locally Lipschitz function. The input u p is such that s i ≤ u p,i ≤ s i for i = 1, • • • , p, with s i , s i ∈ R the saturations limits. The input u θ is such that s i,θ ≤ u θ ≤ s i,θ for i = 1, • • • , p θ , with s i,θ , s i,θ ∈ R the saturations limits.
In this paper we suppose that a reference virtual command w * (t) is generated by a high-level controller and we focus on the optimal design or control of the allocation matrix A(θ) to optimize some cost function V (θ). Given an optimal configuration A * (θ) and a feasible virtual input w * , the instantaneous input u p can be found as:

u p (t) = A * -1 w * (t) p = n u p (t) = A * † w * (t) p > n (4)
where A † is the constrained pseudoinverse operator [START_REF] Sabharwal | Convexly constrained linear inverse problems: iterative least-squares and regularization[END_REF] to guarantee that the solution lies in the feasible saturated u p .

We define two problems:

Design Problem: we suppose that f θ (θ, u θ ) = 0 and p θ = 0, i.e., the configurations θ are constant and there are no inputs to steer θ. The goal is to find a static optimal configuration θ * that maximizes/minimizes a cost function V (θ) ∈ R. The design problem is useful to design some intrinsic properties of the vehicle such as orientation and position of actuators, maximum power of the actuators, etc.. The optimization problem is defined as: max/min θ V (θ) subject to (3)

s i ≤ u p,i ≤ s i (5)
Control Problem: the goal is to find an optimal configuration and configuration input trajectory {θ * (t), u θ (t) * } that maximize/minimize a cost function V (θ) ∈ R. The optimization problem is defined as: max/min

θ(t),u θ (t) t f 0 V (θ(t))dt subject to (3) s i ≤ u p,i ≤ s i s i,θ ≤ u θ ≤ s i,θ (6)
To describe the problem we have to define the sets W ⊆ R n and U ⊆ R p as the virtual input set and the instantaneous input set respectively, i.e., where the virtual input w and the instantaneous input u p vectors can lie. Due to the nature of the saturated input u p , U is a p-orthotope (hyper-rectangle of dimension p) polytope defined in Hrepresentation as:

U = {x ∈ R p |C U x ≤ b U } (7)
with C U ∈ R 2p×p defined as:

C U = I p -I p (8)
with I p an identity matrix of dimension p, and b U defined as:

b U = s 1 • • • s p -s 1 • • • -s p T (9) Due to (3a) W is the image of U subject to linear transformation A(θ)
, which is again a polytope [START_REF] Zhang | Polyhedron under linear transformations[END_REF] defined in H-representation as:

W = {x ∈ W p |C W x ≤ b W } .
(10)

Cost Function Definition

In the following we define suitable cost functions V to be optimized, based on different tasks for the robot.

Manipulability task: a desired virtual input trajectory is given by w * (t). The goal is to find the optimal configuration θ * (t) and the configuration input u θ such that the manipulability index around w * (t) is the maximum, considering as manipulability index the radius of the largest sphere inscibed into the polytope W and centered in w * (t). This goal can be useful to provide robustness to the high level control. The higher the manipulability index, the higher the high level control has margin around the nominal control input w * (t) to compensate external disturbances or disturbances due to non perfect modeling.

If we define as d i the distance of w * (t) to the i-th hyperplane that constraints W, we define the cost function as:

V M (θ) = min(d i (θ)), i = 1, . . . , t, t ≤ 2p (11)
In the evaluation of the cost functions using W, we might have different measurements units (for example [N] and [N/m] for forces and torques, respectively) or preferred directions in W (for example we want higher manipulability in the forces of x-axis). In this case we can simply scale the matrix A by pre-multiplying it with a scale matrix A λ ∈ R n×n defined as:

A λ =   λ 1 0 . . . 0 λ n   (12)
Remark 1. Optimizing for the manipulability task, it is not guaranteed that the desired w * , around which we maximize the radius of the sphere, is feasible w.r.t. the constraints on u p .

Energy task: a desired virtual input trajectory is given by w * (t). From (4) we can compute u p (t). The cost function has to be chosen to minimize the norm of the instantaneous inputs. The cost function to maximize in case of (p = n) is defined as:

V E (θ) = A -1 (θ)w * (t) -1 w * (t) ∈ W V E (θ) = A -1 (θ)w * (t) -1 w * (t) / ∈ W (13)
while in case of (p > n) it is defined as:

V E (θ) = A † (θ)w * (t) -1 w * (t) ∈ W V E (θ) = A † (θ)w * (t) -1 w * (t) / ∈ W (14)
where w * (t) is the obtainable virtual input, closest to w * (t) (can be found with quadratic programming [START_REF] Wolfe | Finding the nearest point in a polytope[END_REF]).

Remark 2. Optimizing for the energy task, it is not guaranteed that the desired w * is feasible using the configuration obtained.

Wrench maximum volume task: the goal is to find the optimal θ * to maximize the volume of W. The cost function to maximize is defined as:

V V (θ) = vol(W) (15) where vol(W) is the volume of W.
The cost function can be defined as weighted sum of different tasks, i.e., we can define

V as V = β 1 V M +β 2 V E + β 3 V V ,
with β i some positive scalar.

OPTIMIZATION PROBLEM SOLUTION

In the following, we present a solution to the optimal problems (5) and ( 6) for the mentioned tasks, and we also propose a generalized solution to consider other possible tasks or a cost function mixing the different tasks.

Manipulability

The solution lies into describing W as function of U and A(θ). The solution is different in case of fully actuation (p = n) or in case of redundant actuation (p > n):

• Fully Actuation (p = n): if A(θ) is invertible, we can compute b W = b U and C W = C U A(θ) -1 (closed form solution). If A(θ) is non-invertible, the cost function V M = 0. • Redundant Actuation (p > n): since A is a projection
on a lower dimension space, the problem is known in literature for not having a closed form solution.

A numerical solution consists in expressing U in its V-representation (the vertex representation), transform the verteces of U via A(θ), compute the Vrepresentation of W as convex hull of the transformed points using [START_REF] Chazelle | An optimal convex hull algorithm in any fixed dimension[END_REF]. The problem is known as the facet enumeration problem [START_REF] Fukuda | Frequently asked questions in polyhedral computation[END_REF].

Once W is defined, we can compute d i as

d i = |C i W w * -b i W | C i W (16)
where the i index refers to the i-th row or the i-th hyperplane defining W.

Energy

• Fully Actuation (p = n): the cost function requires the computation of W which can be done in closed form and the computation of A -1 which again can be done in closed form from A(θ). Some gradient-based algorithms can be used to compute the optimal θ * .

• Redundant Actuation (p > n): both the computation of W and A † cannot be done in closed form, so the resulting cost function is non-differentiable and noncontinuous.

Wrench Maximum Volume

The cost function is directly dependent on the volume of W which can be computed in a numerical fashion by one of the algorithms in [START_REF] Bueler | Exact volume computation for convex polytopes: a practical study[END_REF].

Other Tasks

The solution of the generalized problem, as described in the following, requires only to be able to evaluate the cost function for a particular configuration θ. Hence complex cost function can be built to be optimized.

Optimization Problem Solution

In general, the cost function can be non-linear, nondifferentiable (for the computation of W and/or A † ) and in high dimensions, so standard optimization algorithms based on differentiation [START_REF] Ruder | An overview of gradient descent optimization algorithms[END_REF] or complete solutions cannot be applied. A general solution can instead be derived using genetic algorithms (GA) [START_REF] Srinivas | Genetic algorithms: a survey[END_REF]. In case of the Design problem a single GA can be executed to find the optimal configuration θ * . In the Control problem, instead, we discretize w * (t) with a given time step ∆t, and solve a GA for each step. Moreover, we need to consider the dynamic equation (3b) in the solutions of the GA. We have then to introduce the concept of Initial Genetic Algorithm (IGA) and on Trajectory Genetic Algorithm (TGA). The IGA is used in the Design problem and to find the initial configuration θ * (0) for the Control problem. The TGA is instead used along the trajectory to find the θ * (t > 0). Since the TGA has to obey (3b) we need to find a suitable mechanism to find feasible solutions of θ * (t). For sake of compactness we call TGA(k∆t) as TGA k and we use the subscript k to indicate a k∆t. The choice of the initial population, selection and mutation of each TGA k should guarantee that θ * k = θ * k-1 + f θ (u θ,k ). Suitable strategies to pick a feasible θ * can be derived from kino-dynamic planning LaValle and J. J. [START_REF] Lavalle | Randomized kinodynamic planning[END_REF] and are model specific and GA specific.

SIMULATIONS

In the following, we present the results of two case studies, one for the control problem and one for the design problem.

Control Case: redundantly-actuated planar robot

The platform described in [START_REF] Nainer | Hierarchical control of the overactuated rospo platform via static input allocation[END_REF] and depicted in Figure 1, is composed by a rigid rectangular frame sustained by omnidirectional passive spherical wheels and n actuating modules. Each module consists of a turret, which is orientable by means of a servo motor. Each turret carries a propeller driven by a BLDC motor. We consider the platform to have 4 actuator modules.

We can define equation (3) for this particular robot by defining w, θ, A, u p , f θ , u θ . In particular w = [F x , F y , τ ] T is the wrench vector, where F x is the body force along xaxis, F y is the body force along y-axis and τ the torque around the vertical axis;

θ = [θ 1 , • • • , θ 4 ] with θ i the angle of the i-th turret; u p = [k ω ω 2 1 , k ω ω 2 2 , k ω ω 2 3 , k ω ω 2 4 ] T ,
where k ω is an aerodynamic constant and ω i the angular velocity of the i-th propeller. f (θ, u θ ) = u θ , u θ = [u θ,1 , u θ,2 , u θ,3 , u θ,4 ] T where u θ,i is the control input for the turret's angle. The turret model is a single integrator system. Finally the allocation matrix is defined as:

A(θ) =    c(θ 1 ) • • • c(θ 4 ) s(θ 1 ) • • • s(θ 4 ) (Πr 1 ) T c(θ 1 ) s(θ 1 ) • • • (Πr 4 ) T c(θ 4 ) s(θ 4 )    ( 17
)
where Π = 0 1 -1 0 , r i ∈ R 2 the position of the turret w.r.t. the center of gravity of the platform and c(•), s(•) indicating the cos and sin function respectively.

A desired wrench trajectory w * (t), derived from the model inversion to track a 8-shaped trajectory, is given by:

w * (t) = F * x F * y τ * = -M f 2 1 R cos(f 1 t) -M f 2 2 R sin(f 2 t) 0 ( 18
)
with R, f 1 , f 2 some parameters to define the size and the frequency of the trajectory and M the mass of the robot.

The following results are for the manipulability and energy tasks, whose cost functions are given in ( 11) and ( 14).

Because we are in the redundantly-actuated problem, the Genetic Algorithm solution is employed.

The population is defined as N p individuals p i defined by a series of chromosomes. In particular the individuals p i are defined as the sets {θ 1 , θ 2 , θ 3 , θ 4 }, i.e., the chromosomes are the parameters θ to be optimized. The fitness function is defined by the cost function ( 11) or ( 14), depending on the task. The steps of the genetic algorithm are defined as the following:

• Initial Population (IGA): the initial population is selected randomizing

θ i ∈ [-π, π]. • Initial Population (TGA): the initial population is selected randomizing θ i = θ o i ± δ θ,i
, where θ o i is the optimal value of θ i from the previous iteration and δ θ,i is the maximum variation of θ i in ∆t given a limitation of the u θ,i . Since (3b) is a single integrator model δ θ,i = u θ,i .

• Selection: the selection mechanism used is the fitness proportionate selection [START_REF] Blickle | A comparison of selection schemes used in evolutionary algorithms[END_REF]. • Crossover: the crossover mechanism consists in taking 2 chromosomes (two θ i ) from each parent. • Mutation (IGA): in the mutation process at the initial iteration, each chromosome has a probability P M of being randomized again in [-π, π]. • Mutation (TGA): in the mutation process along the trajectory, each chromosome has a probability P M of being increased by ±δ θ,i .

The parameters used in the simulations are gathered in Table 1, where we indicate with N G the number of generations for each GA solution. Fig. 3. The obtainable wrench polytope after the solution of IGA (t = 0). The inscribed ball is the so called manipulability index, which measure the possibility of generating virtual inputs around a reference virtual input for robustness purpose.

In Figure 2 the fitness value for the IGA is depicted.

The fitness converges to a steady value in around 500 generations.

In Figure 3 we can see the obtainable virtual input w polytope (in light blue) with the optimal configuration solution given by the IGA (t = 0), for both the manipulability and the energy tasks. Inside the polytope is depicted the sphere indicating the manipulability index centered around the initial desired wrench. We can notice that in the configuration optimized for the manipulability we can obtain a much wider range of wrenches in all direction, while in the configuration optimized for the energy we can obtain a large range in x-axis (the direction of w * (t = 0)) and very small in the other directions.

In Figure 4 the system is tested with the high level controller developed in [START_REF] Nainer | Hierarchical control of the overactuated rospo platform via static input allocation[END_REF]. We can see the different behavior of the turrets. In the energy task the turrets are pointing mainly in the direction of the requested wrench (the blue line), but the manipulability index is very low (the mean V M = 0.41), such that the platform cannot compensate big external disturbances. In the manipulability task the turrets are more distributed, the manipulability index much higher (the mean V M = 4), but the energy consumption is almost double (mean V E = 0.221 against V E = 0.39).

Design Case: Fully Actuated Hexa-rotor

For the design case we consider a hexa-rotor unmanned aerial vehicle (UAV). Typical hexa-rotors have 6 equidistant motor-propeller actuators which point in the vertical z-axis. Some particular hexa-rotors instead have tilted actuators to generate forces in 3 axes and torques in 3 axes. Inspired by those designs, the goal is to define the optimal orientation of each independent propeller to optimize some cost functions. We define r i ∈ R 3 the position vector of the i-th actuator w.r.t. the center of gravity of the vehicle. We define v i ∈ S 2 the orientation of the i-th actuator in body frame, where the group S n defines the unit sphere group S n = x ∈ R n+1 | x = 1 . The vector v i describes where the propeller is pointing. We define σ i ∈ {-1, +1} the spinning direction of the i-th propeller. We can define equation (3) for this particular robot by defining w, θ, A, u p , f θ , u θ . In particular w = [F, τ ] T is the wrench vector, where F ∈ R 3 is the body force in three dimensions x, y, z and τ ∈ R 3 the torque in three dimensions.

θ = [v 1 , • • • , v 6 ]. u p [k ω ω 2 1 , • • • , k ω ω 2 6 ] T ,
where k ω is aerodynamic constant and ω i the angular velocity of the i-th propeller. We consider the design problem, so θ = 0. Finally the allocation matrix is defined as:

A(θ) = θ γ 1 • • • γ 6 (19) with γ i = kθ i σ i + S × (r i )v i , i = 1, • • • , 6
and k an aerodynamic constant that relates the thrust of the propeller with the torque produced around the rotation axis.

In the following, we present the results of the optimization problem to optimize the manipulability and energy task for the hovering condition (the typical working point), that is w * = [0, 0, M g, 0, 0, 0] T , with M the mass of the vehicle and g the gravitational acceleration. To simplify the problem, only the directions of the propellers are optimized, while r i and σ i are decided a priori. Remark 3. The generalized approach with GA is so powerful that one can optimize almost any parameter, including the number of actuators, the position r i , the spinning direction σ i , the saturation levels s i , s i , etc. The steps of the genetic algorithm are defined as the following:

• Initial Population (IGA): the initial population is selected randomizing v i ∈ S 2 . • Selection: the selection mechanism used is the fitness proportionate selection [START_REF] Blickle | A comparison of selection schemes used in evolutionary algorithms[END_REF]. • Crossover: the crossover mechanism consists in taking 3 chromosomes (three v i ) from each parent. • Mutation (IGA): in the mutation process each chromosome has a probability P M of being randomized such that v i ∈ S 2 .

The parameters in Table . 2 were used in the simulations.

We considered 3 optimization problems:

• a) the cost function (11) for manipulability.

• b) the cost function ( 11) for manipulability but with a scale matrix A s = diag {1, 1, 1, 3, 3, 8} to scale the torques w.r.t. the forces. • c) the cost function ( 13) for energy and w * a feasible wrench (i.e., V M > 0).

In Fig. 5 the fitness function for the problems a) and c), while in Fig. 6 the visual representation of the optimal configurations. Finally, in Table 3 the optimized configurations for the problems a), b), c).

It is easy to note how the orientations of the configurations optimized for the manipulability are much more tilted, while the ones optimized for the energy are almost vertical. This corresponds to the capability of resisting an external wrench coming from any direction (in the first two cases), or just the one produced by gravity.

Furthermore, it is worth to notice that the configurations shaped for the manipulability use around twice and four times the energy of the configuration optimized for the energy. On the other hand, the configuration optimized for the energy can generate a very small wrench around the hovering condition.

Given the consistent computational time required by the solving of the optimization problem with Genetic Algorithms w.r.t. other state-of-the-art techniques, we limited its usage to offline optimization, with a feasible w * precomputed for all the trajectory.

CONCLUSION AND FUTURE WORKS

In this paper we considered fully-actuated or redundantlyactuated system, where a virtual input can be considered as a linear combination of the instantaneous inputs of the system via an allocation matrix that depends on some configuration parameters. An optimization problem was defined, to optimize the configurations of to maximize or minimize a function, as manipulability, energy consumption and volume of the obtainable virtual input. Due to the non-differentiable nature of the costs functions, the general problem was approached with Genetic Algorithms. The proposed solution allowed to optimize static parameters of the system for a design problem, as well as to find the optimal trajectory of the configurations in the control problem, where the configurations can be controlled by mean of some inputs. The approach was tested on two scenarios: finding the reconfiguration of the orientation of some actuators in an redundantly-actuated mobile robot, and on the design of the optimal configuration for a generalized hexa-rotor with arbitrary propeller orientation. The preliminary results were very promising.

Future works include: the definition of cost functions for other tasks, the application of the method to other robotic platforms for the design problem and the implementation of the allocation strategy on the real robotic platform in [START_REF] Nainer | Hierarchical control of the overactuated rospo platform via static input allocation[END_REF].

Table 3. Optimal configuration for the different tasks.

Fig. 1 .

 1 Fig. 1. The ROSPO platform.

Fig. 2 .

 2 Fig. 2. Fitness functions for the IGA problem in both manipulability and energy tasks.

Fig. 4 .

 4 Fig. 4. The trajectory of the platform with the high level controller with the solution of the TGA in manipulability and energy task. The direction of the turrets along the trajectory is indicated by orange arrows starting from the platform.

Fig. 5 .

 5 Fig. 5. The Fitness evolution for the manipulability and energy tasks.

Fig. 6 .

 6 Fig. 6. The configurations of the hexarotor propellers for the three optimization problems.

Table 1 .

 1 [START_REF] Ryll | Modeling and control of FAST-Hex: a fully-actuated by synchronizedtilting hexarotor[END_REF] Parameters for GCM weight functions.

	R	f 1	f 2	∆t	Np	N G	P M	s i	s i
	8	0.3 0.3	0.1(s)	90	500	0.05	5.8 (N)	0(N)

Table 2 .

 2 Parameters for the hexa-rotor simulation.

	N p	N G	P M	s i	s i	r i	σ i
	500 10000 0.05	6.5 (N)	0(N)	0.7(m)	-1 i

Preprint version,

12th SYROCO 2018, Budapest, Hungary (2018)

Manipulability

Manipulability Scaled Energy