The present paper is concerned with the dynamic modeling and design of control laws for a small non-rigid multi-rotor airship constituted of an oblate-spheroid helium balloon coupled with an electric-powered hexa-rotor airframe. The vehicle is assumed to operate in windless and low- speed conditions. A six-degree-of-freedom nonlinear dynamic model is derived for it using the Newton-Euler approach and considering, among other common e orts, a restoring torque due to the displacement of the balloon's center of buoyancy above the vehicle's center of mass and the added-mass e ect resulting from the air-structure interaction. Using the derived model and as- suming a time-scale separation between the translational and rotational dynamics, the attitude and position control laws are designed separately from each other. Both laws are formulated using feedback linearization combined with control input saturation within appropriate parallelepipedal sets, which are carefully chosen to respect pre-de ned bounds on the control torque, control force and inclination angle. The e ect of temperature and pressure uctuations is taken into account through a parametric probabilistic approach, where Maximum Entropy Principle is used to con- struct a physically consistent stochastic model and Monte Carlo method is used as stochastic solver to propagate the uncertainties through the system. Extensive simulation results show the e ectiveness of the proposed control system and quanti es the uncertainty of its performance over a wide range of local temperature and pressure. Keywords: multi-rotor airship, multi-rotor aerial vehicle, ight control system, uncertainty quanti cation. Preprint submitted to ISA Transactions January 21, 2019 1. Introduction The traditional airships are lighter-than-air (LTA) aircrafts constituted of a streamlined or elliptical envelop (or hull) lled with a low-density gas (generally, helium), a pair of propellers installed through the gondola on its bottom part and a tail equipped with rudder and elevator. They mainly rely on the aerostatic lift, resulting from the di erence between the air and inner gas densities, for opposing the gravity and keep ying, which provides them with better endurance, payload capacity, and energy consumption characteristics compared to heavier-than-air (HTA) vehicles (such as airplanes and helicopters) [1]. In the last three decades, thanks to advances (and cost reduction) in sensor technology, em- bedded computers, and materials, the airships have resurged in academia, mainly in non-rigid (blimp), small, and unmanned versions [2]. Such interest renewal in turn gives rise to many works concerned with the dynamic modeling and identi cation [3, 4] as well as the ight control [5, 6] of such aerial vehicles. Despite of the aforementioned advantages of the traditional airships over the HTA vehicles, they su er from a weak controllability in the vertical and lateral translational degrees of freedom (DOF), besides of being more vulnerable to undesirable aerodynamic e ects, such as those due to the drag and the added mass [1]. The interplay between these drawbacks has motivated the design of a number of small and electric hybrid airships [7, 8, 9, 10] provided with an augmented maneuverability by an unconventional combination of actuators. In particular, the reference [7] presents a n-less ellipsoidal blimp equipped with four vectoring rotors; the reference [8] analyzes a small streamlined blimp actuated by one vertical and two horizontal xed rotors; the reference [9] is concerned with a ying-wing-shaped hull equipped with four vectoring double rotors; and the reference [10] comes up with a cylindrical blimp with conventional actuators, whose rotor axis can be controlled to slide along a longitudinal keel. The present paper investigates the dynamic modeling and ight control of a hexa-rotor airship, which is a new hybrid aircraft consisting of an oblate-spheroid balloon lled with helium, attached to a hexa-rotor airframe (see a picture in Figure 1). Their two parts are strapped down by ropes at six contact points. The balloon is made on polyurethane and has 2 :5 m of horizontal diameter and 2 1:6 m of height. With a payload of 4 :42 kg, the vehicle has a buoyancy ratio of 0 :7. In a qualitative pre-analysis, one can say that, in low speed and windless condition, the hexa-rotor airship has potentially a better ight duration and payload capacity than a pure hexa-rotor vehicle, since its aerostatic lift is sized to counteract part of the weight aboard it. On the other hand, compared to other airships, the one investigated here has a full maneuverability in the position and heading DOFs, which is a characteristic inherited from conventional multi-rotor vehicles. Figure 1: The hexa-rotor airship. It is patented according to reference [11]. In particular, the present paper is concerned with the dynamic modeling and design of ight control laws for the hexa-rotor airship depicted in Figure 1. A nonlinear six-DOF dynamic model is derived for this vehicle using the Newton-Euler approach, considering its operation in a low- speed and windless condition. On the basis of this model and considering a time-scale separation between the translational and rotational dynamics, a hierarchical ight control architecture with the attitude control loop put inside the position control loop is adopted here. This architecture is commonly found in the literature of control applied to multi-rotor aerial vehicles (MAV) [12, 13, 14]. Here, on the basis of the derived model, the attitude and position control laws are separately designed using feedback linearization and considering the saturation of the control vector within appropriate parallelepipedal sets that ensure the satisfaction of design bounds on the control torque 3 and force, as well as on the vehicle's inclination angle. The design of position and attitude control laws for the hexa-rotor airship, even if it is intended to operate indoors and in low speed, must face two paramount dynamic aspects not found in the MAV control problem. The rst one is the restoring torque, which appears as a consequence of the displacement between the balloon's center of buoyancy (CB) and the vehicle's center of mass (CM) [15]. For a xed CB-CM displacement, the larger the inclination angle of the vehicle with respect to the local vertical, the larger the magnitude of the restoring torque. Therefore, for a given sizing of the rotor set, the attitude controller must respect a maximum bound on the inclination angle for the control system to maintain its e ectiveness to counteract the restoring torque (and more). The side e ect of such inclination constraint is a low lateral acceleration capability. The second aspect is the added-mass e ect [16], which consists in the reaction e orts (torque and force) of the surrounding air on the vehicle's body. It turns out that the added-mass e ect causes an additional coupling between the translational and rotational dynamics, besides an apparent increasing in the vehicle's mass and inertia. The above two aspects are considered here. Moreover, an uncertainty quanti cation analysis is carried out to assess the variation in the closed-loop performance as a result of the deviation of the operating temperature and pressure with respect to their nominal values; note that such deviation causes a change in the air and helium densities and, as a consequence, in the aerostatic lift and restoring torque. From the above justi cation, we can assert that the main contributions of this paper are: 1) the dynamic modeling of the hexa-rotor airship for windless and low-speed conditions; 2) the design of a simple but e ective control scheme for this vehicle, considering the aerostatic e orts; and 3) the uncertainty quanti cation study for assessing its closed-loop performance over a speci ed range of possible operating temperatures and pressures. The remaining text is organized in the following manner. Section 2 derives the six-DOF dynamic model. Section 3 is concerned with the design of the nonlinear attitude and position control laws, as well as the stability analysis of the corresponding closed-loop dynamics. Sections 4 and 5 evaluates the overall system numerically. Particularly, Section 4 adopts a nominal scenario, whereas Section 5 carries out a stochastic uncertainty quanti cation analysis considering a random variation of both the local temperature and pressure conditions over given ranges. Finally, Section 6 concludes 4 60𝑜 30𝑜 𝐵 𝑧 B 𝑥 B 𝑦 B 𝐺 𝑧 G 𝑥 G 𝑦 G 1 2 3 4 5 6 Figure 2: Illustration of the hexa-rotor airship and representaiton of the Cartesian coordinate systems. the paper. 2. Dynamic Modeling This section derives rotational and translational equations of motion, as well as actuator models for the hexa-rotor airship, considering the most relevant e orts in low-speed indoor ight scenarios. We start with preliminary de nitions in Subsection 2.1, then we model the rotor dynamics and control e orts in Subsection 2.2, the restoring torque and aerostatic lift generated by the balloon in Subsection 2.3, the vehicle's translational dynamics in Subsection 2.5, and its rotational dynamics in Subsection 2.4. Finally, the added-mass e ects are treated in Subsection 2.6. 2.1. Preliminary De nitions, Notation, and Nomenclature De ne two Cartesian coordinate systems (CCS) as illustrated in Figure 2. The body CCS, SB,fB; ^xB;^yB;^zBgis attached to the vehicle's body with the origin at the vehicle's center of massB, the ^xBaxis is pointing forward, aligned with the bisectrix of the separation angle between 5 the rotors 1 and 2, the ^ zBaxis is pointing upward, normal to the airframe, and the ^ yBaxis completes a right-handed coordinate system. The ground CCS, SG,fG; ^xG;^yG;^zGgis xed to the ground at a known point G, with the ^zGaxis pointing upward vertically. For our purposes, SGcan be considered as an inertial frame and will be adopted as such. Denote the set of real numbers by Rand the set of positive and non-negative real numbers by R>0 and R0, respectively. Moreover, denote the special orthogonal group by SO(3). Denote arbitrary sets by uppercase calligraphic letters, e.g.,AandB. The Minkowski sum of two sets AandBis de ned byAB,fa+b:a2A;b2Bg . To clarify the vector notation adopted here, consider two di erent kinds of vectors: physical (geometric) vectors and algebraic vectors. Physical vectors are denoted by italic letters with a right arrow superscript, e.g.,~ a. The corresponding algebraic vector resulting from the projection of ~ aonto an arbitrary CCS SAis denoted by the same letter, but in a bold format, with the subscript A, i.e.,aA2 R3. Matrices are denoted by bold uppercase letters, e.g.,A. Let~ abe an arbitrary physical vector. The relation between its representations aA andaBisaA=DA=BaB, where DA=B2SO(3) is the attitude matrix of SAw.r.t.SB. The inverse ofDA=B(which coincides with its transpose) is sometimes denoted by DB=A. The Euclidean norm ofaAis denoted bykaAk. Now, let~ aA=Brepresent a physical quantity of SArelative toSB,e.g., ~ rB=Gdenotes the relative position of SBw.r.t.SG. Let eidenote ann-dimensional (for some n) standard unit vector. The ith component of ais denoted by aiand sometimes by eT ia. Now, consider two physical vectors ~ aand~bas well as the respective algebraic representations in SA, aA= [a1a2a3]TandbA, respectively. Denote the SArepresentation of the vector product ~ a~b by the matrix multiplication [ aA]bA, where [ aA] is the following skew-symmetric matrix: [aA],2 666640a3a2 a3 0a1 a2a1 03 77775: (1) The text also distinguishes between a variable representing an e ective quantity, e.g.,a, and the corresponding command (or desired value). The latter is denoted with the same letter as the e ective quantity, but including an overbar; e.g., a. 6 2.2. Rotor Dynamics and E orts The set of six rotors equipping the airframe is responsible for generating the control forces and torques as described here. The ith rotor individually produces a thrust force and a reaction torque on the airframe along the ^ zBaxis with magnitudes denoted by fiand i, respectively. We describe these e orts by the following aerodynamic models [17]: fi=kf!2 i; (2) i=k !2 i; (3) i= 1;:::;4, wherekfis the thrust force coecient, k is the reaction torque coecient, and !iis the rotation speed of the ith rotor. The rotor dynamics can be modeled by the following rst-order linear model: _!i=1 !!i+k! !!i; (4) where !i2[0;!max] is the rotation speed command of the ith rotor,k!is the speed coecient, and !is the rotor time constant. The rotation upper bound  !maxis assumed to be known. Consider that all the six thrusts fipoint upward. Moreover, consider that the reaction torque 1is positive, 2is negative, 3is positive, and so on. Therefore, one can show that the magnitude Fcof the resulting control force and the SBrepresentation Tc Bof the resulting control torque are given by 2 4Fc Tc B3 5=f; (5) where f,[f1f2::: f 6]Tand 7 𝐵 𝑥 B 𝐺 𝑧 G 𝑥 G 𝑦 G 𝑧 B 𝐻 𝐹 𝑏 𝑑 Figure 3: The restoring torque. ,2 666666641 1 1 1 1 1 l=2l=2ll=2l=2l lp 3=2lp 3=2 0lp 3=2lp 3=2 0 kk kk kk3 777777752 R46; (6) withldenoting the length of each vehicle's arm, and k,k =kf. 2.3. Aerostatic Lift and Restoring Torque This subsection models two crucial e orts generated by the balloon. One is an aerostatic lift force~Fband the other one is a restoring torque ~Tb. The force ~Fbis explained by the Archimedes' Principle, which says that the buoyancy force always points upwards parallel to the local vertical and its magnitude is equal to the weight of the air volume displaced by the balloon. Therefore, Fb G=Vg aire3; (7) 8 whereVis the volume of the balloon (and of the gas as well), gis the gravitational acceleration, andairis the air density. Denote the magnitude of ~FbbyFb,Vg air. On the other hand, the restoring torque ~Tbis an e ort acting about the vehicle's CM, which appears as a consequence of the displacement dbetween the balloon's CB and the vehicle's CM. Figure 3 illustrates the balloon-airframe connection. Assume that this connection is rigid. In this case, the CB, which is placed at point H, is xed w.r.t.SB. From the illustration, one can immediately write ~Tb= (d^zB)~Fb: (8) By representing equation (8) in SB, we nally have Tb B=d[e3]DB=Ge3Fb; (9) where DB=G2SO(3) is the attitude matrix of SBw.r.t.SG. 2.4. Translational Motion By invoking the Second Newton's Law considering all the vectors represented in SG, one can immediately write 5  mt+mhrB=G G=Fg G+Fb G+Fc G+Fd G; (10) wheremtis the total mass of the vehicle without lifting gas and including the payload, mh, heliumVis the helium mass, with helium representing the helium density, rB=G G2 R3is theSG representation of the position of the vehicle's center of mass Bw.r.t.G,Fg Gis theSGrepresentation of the gravitational force, Fb Gis theSGrepresentation of the aerostatic lift force, Fc Gis theSG representation of the control force, and Fd Gis theSGrepresentation of the (unknown) disturbance force. The force Fb Gis given by equation (7), while Fg GandFc Gare modeled by Fg G= mt+mh ge3; (11) Fc G=Fc DB=GT e3: (12) 9 By replacing equations (11)-(12) into (10), one can obtain rB=G G=Fc mt+mhnG+Vg air mt+mhg e3+1 mt+mhFd G; (13) where nG2 R3is the transpose of the third line of DB=G, which corresponds to the SGrepresen- tation of the unit (geometric) vector ^ zB. 2.5. Rotational Motion The kinematic equation of the rotational motion of SBw.r.t.SGis given in SO(3) by _DB=G=h B=G Bi DB=G; (14) where B=G B2 R3is theSBrepresentation of the vehicle's angular velocity w.r.t. SG. Assume that the vehicle has a rigid structure and consider SGas an inertial frame. Therefore, the Euler's Law gives _HB+h B=G Bi HB=Tc B+Tb B+Td B; (15) where HBis theSBrepresentation of the total angular momentum of the vehicle, Tc Bis the SBrepresentation of the control torque (see equation (5)), Td Bis theSBrepresentation of the (unknown) disturbance torque, and Tb Bis theSBrepresentation of the balloon restoring torque (see equation (9)). Considering the rotation of both the body and the propellers and noting that the latter rotates much faster, the total angular momentum HBcan be written as HB=J B=G B+Jr6X i=1(1)i!ie3; (16) where J2 R33is the inertia matrix of the vehicle and Jr2 Ris the moment of inertia of each rotor about ^ zB. Therefore, by replacing equation (16) into equation (15), one can obtain the dynamic equation of the rotational motion of SBw.r.t.SG: 10 _ B=G B=J1" J B=G B+Jr6X i=1(1)i!ie3! # B=G BJ1Jr6X i=1(1)i_!ie3+ J1Tb B+J1Tc B+J1Td B:(17) 2.6. Added Mass The added-mass (also called virtual or apparent-mass) forces and torques are aerodynamic reaction e orts from a surrounding uid on a body which tries to accelerate or decelerate inside it; note that the displaced mass of uid is accelerated or decelerated together with the body interacting with it (see, e.g., [18, 19]). This e ect is only relevant if the mass of the displaced uid is signi cant compared to the body's mass, which is generally the case for conventional or hybrid airships [20, 9, 16]. It is also worth noting that it appears even in low-speed and windless conditions. From the Kirchho equation (for a body moving in an ideal uid) one can show that the added-mass force and torque are, respectively, given by [18] 2 4Fa B Ta B3 5=Ma2 4_vB=G B _ B=G B3 5 2 4h B=G Bi M11vB=G B+M12 B=G B h vB=G Bi M11vB=G B+M12 B=G B +h B=G Bi M21vB=G B+M22 B=G B3 5;(18) where Ma,2 4M11M12 M21M223 52 R66; (19) is the added mass in SB,Mij2 R33,8i;j2f1;2g. By writing the added-mass force as well as the vehicle's linear velocity in SG, equation (18) 11 yields Fa G=DG=BM11DB=G_vB=G GDG=BM12_ B=G B+DG=BM11h B=G Bi DB=GvB=G G DG=Bh B=G Bi M11DB=GvB=G G+M12 B=G B ;(20) Ta B=M21DB=G_vB=G GM22_ B=G B+M21h B=G Bi DB=GvB=G G h DB=GvB=G G i M11DB=GvB=G G+M12 B=G B h B=G Bi M21DB=GvB=G G+M22 B=G B :(21) Before computing Mafor the oblate-spheroid balloon under consideration, it is convenient to obtain the added mass Ma,diag( M11;M22) with respect to a CCS which is parallel to SBand has its origin at CB, using the procedure of [21], resulting in M11=2 66664m0 0 0 m 0 0 0 m33 77775;and M22=2 66664J0 0 0J0 0 0 03 77775; (22) where 12 m,4 3a2bair 0 2 0; (23) m3,4 3a2bair 0 2 0; (24) J,4 15a2bair a2+b2 0 0 2; (25) 0,p 1e2 e3 sin1eep 1e2 ; (26) 0,2p 1e2 e3ep 1e2sin1e ; (27) e,r a2b2 a2: (28) Finally, to obtain Ma(inSB), one can use the coordinate transformation (see [22], p.27) Ma=UTMaU (29) with U,2 4I3d[e3] 033 I33 5: (30) Therefore, to account for the added mass, the new external inputs represented by Fa GandTa B have to be considered in equations (13) and (17), respectively. 3. Flight Control System A hierarchical control strategy is adopted here, as illustrated in Figure 4. In this strategy, the ight control is realized by two nested control loops, where the inner loop is responsible for the attitude control, while the outer loop performs position control. The position controller receives an external position command rB=G G2 R3as well as feedbacks of the vehicle position rB=G Gand velocity vB=G G. On the other hand, it produces the SGrepresentation of the thrust command vector Fc G,FcnG2 R3, where Fc2 Ris the total thrust magnitude command and nG2 R3 (a unit vector) can be seen as a two-DOF attitude command that, together with the external 13 Flight Control System Position Controller Attitude Command Generator Attitude Controller Control Allocation Hexa -Rotor Airship 𝑓 𝑖 𝑖=1,…,6 𝐫GB/G ,𝐯GB/G 𝐫 GB/G 𝜓 𝐧 G 𝐹 𝑐 𝐃 B/G 𝐓 Bc 𝐃B/G,𝛀BB/G Figure 4: Architecture of the hexa-rotor airship control system. heading command  2 R, composes the three-DOF command input DB=G2SO(3) of the attitude controller. The latter also receives feedback of the vehicle's attitude DB=Gand angular velocity B=G B, and computes the torque command Tc B2 R3. The control allocation block is responsible for generating the individual thrust commands fi2 R0,i= 1;:::;6, from Tc Band Fc. This architecture is based on the assumption that there is a time-scale separation between the (slower) closed-loop translational and (faster) rotational vehicle dynamics [12]. This assumption is ensured in practice by tuning the attitude control loop to converge much faster than the position control loop. Under such assumptions, when designing the position control law, one can assume that the actual attitude DB=Gconverges to the corresponding command DB=Ginstantaneously, i.e., DB=GDB=G. On the other hand, when designing the attitude control law, the attitude command DB=Gis assumed to be constant, or equivalently, the angular velocity command  B=G B2 R3is assumed to be zero. On the basis of these conditions, the position and attitude control laws can be designed separately. 3.1. Position Controller First, the design model is obtained from equation (13), by assuming that: 1) the disturbance force Fd Gis negligible, 2) the actual control force magnitude Fcis identical to the corresponding 14 command Fc, and 3) DB=G=DB=G(time-scale separation). The resulting design model is rB=G G=1 mt+mhFc G+Vg air mt+mhg e3: (31) Suppose that Fc Gis bounded within a parallelepipedal set from Fmin,[Fmin 1Fmin 2Fmin 3]T2 R3 toFmax,[Fmax 1Fmax 2Fmax 3]T2 R3. The position controller proposed here is the saturated law Fc G=[Fmin;Fmax]( p); (32) where p[ p 1 p 2 p 3]T2 R3is given by p, mt+mh gVg air e3+ mt+mh K1 rB=G GrB=G G  mt+mh K2_rB=G G;(33) with the matrices K1;K22 R33representing the controller gains (assumed to be diagonal), and [Fmin;Fmax]( p),2 66664[Fmin 1;Fmax 1]( p 1) [Fmin 2;Fmax 2]( p 2) [Fmin 3;Fmax 3]( p 3)3 77775; (34) [Fmax l;Fmax l]( p l),8 >>>>< >>>>:Fmin l; p l<Fmin l p l; p l2[Fmin l;Fmax l] Fmax l; p l>Fmax l; l= 1;2;3: (35) Note that the proposed controller (32)-(33) is such that, if no saturation is active, it cancels the second term on the right-hand side of equation (33), remaining a double-integrator closed- loop dynamics controlled by the proportional-derivative actions appearing in the last two terms of equation (33). In saturation-free conditions, it is straightforward to show asymptotic stability of the proposed translational control loop using linear time-invariant control methods. Subsection 3.5 presents a stability analysis considering the saturation. By constraining the control input Fc Gto keep inside a parallelepipedal set, one can deal with 15 𝑧 R 𝑦 R 𝑥 R 𝐅 G𝑐 𝐹3min 𝐹3max 𝐵 𝐹12max 𝑥 R 𝑦 R 𝐅 G𝑐 𝐵 𝑥 R 𝑧 R 𝐵 𝜑max 𝐹3min (a) (b) (c) Figure 5: Constraint set of the position control input. (a) Perspective view showing the control force command inside the parallelepipedal set. (b) Superior view showing the bounds Fmax 1 =Fmax 1 =Fmax 12. (c) Lateral view showing the maximum inclination angle 'max. bounds on its magnitude Fcas well as on its inclination angle w.r.t. the local vertical. To illustrate it, de ne a reference CCS SR,fB; ^xR;^yR;^zRgwith origin at Band axes parallel to the SGones, as depicted in Figure 5. From the geometry of the constraint set, one can see that Fmin 1=Fmax 1, Fmin 2=Fmax 2, andFmax 1=Fmax 2=Fmax 12=Fmin 3tan'max, where'maxis the maximum allowed inclination angle of Fc Gw.r.t. ^zR. 3.2. Attitude Controller Here the design model is obtained from equation (17) by: 1) neglecting the disturbance torque Td B; 2) replacing the actual control torque Tc Bby the corresponding command Tc B; 3) considering that the rotor dynamics are so fast that one can assume !i= !i,8i; and 4) assuming that _ !i0, 8i. The resulting model is _ B=G B =J1" J B=G B+Jr6X i=1(1)i!ie3! # B=G B+dFbJ1[e3]DB=Ge3+J1Tc B:(36) Suppose that the torque command Tc Bmust be bounded from Tmax2 R3toTmax, [Tmax 1Tmax 2Tmax 3]T. The attitude controller proposed here is the saturated law Tc B=[Tmax;Tmax]( a); (37) 16 where a,dFb[e3]DB=Ge3+h B=G Bi J B=G B+Jrh B=G Bi e36X i=1(1)i!i JK3"JK4 B=G B;(38) with"2 R3containing the Euler angles (in the 1-2-3 sequence) corresponding to the attitude control error ~D=DB=G DB=GT ,K3;K42 R33are controller (diagonal-matrix) gains, and the saturation function [Tmax;Tmax](:) is as de ned in equations (34)-(35). One can note that the proposed attitude control law (37)-(38) is such that, if no saturation is active, it cancels the rst three terms on the right-hand side of equation (38), remaining a feedback-linearized closed-loop dynamic model controlled by the proportional-derivative actions appearing in the last two terms of (38). Subsection 3.5 presents a stability analysis. 3.3. Attitude Command Generation In order to generate the three-DOF attitude command DB=R, note that, by de nition, its third line is the transpose of nG. Then, consider the formula to convert from Euler angles (in the 1-2-3 sequence) to the attitude matrix (see [23], p. 52): DB=R=2 66664      scs cc 3 77775; (39) where s and c stands for sin and cos, respectively, and compute  andfrom  =ataneT 2nG=eT 3nG; (40) = asin eT 1nG: (41) 17 Finally, considering the external heading command  , one can compute DB=R=2 66664c c c ss + s c c sc + s s s cs ss + c c s sc + c s scs cc 3 77775: (42) 3.4. Control Allocation Consider equation (5) again; it relates the e ective resultant e orts FcandTc Bwith the six individual rotor forces (stored in f). We can consider that the respective e ort commands are interrelated in the same way, i.e.,2 4Fc Tc B3 5=f; (43) where f,f1f2:::f6T. The control allocation pursued here is to nd a solution fto the system of linear equations (43). Since this is an undetermined system (with in nitely many solutions), we choose here the one that minimizes kfk2, considering (43) as a set of equality constraints. The explicit solution to such a problem is well known; by using Lagrange multipliers, it can be found to be f=T T12 4Fc Tc B3 5: (44) The above control allocation can be seen as an energy-ecient solution, since it minimizes the squared norm of f. Moreover, considering that the commands FcandTc Bbelong to the admissible set of the above optimization problem (meaning that the problem has at least one solution), the resulting control allocation is exact, i.e., it satis es equation (43). 3.5. Stability Analysis Still considering the time-scale separation assumption, this subsection presents stability proofs for both the position and attitude control loops. 18 3.5.1. Position control Denote the admissible-force set by F, F2 R3:FminFFmax (45) and de ne the tracking error ~r,rB=G GrB=G G, the viable set Lp,~r;_~r 2 R6: p2F ; (46) as well as the maximal 0-centered Euclidean ball inside Lp L p,max Bs:t:BLp; (47) whereB,fx2 R6:kxkg. The following theorem gives a stability result for the closed-loop system (31)-(33). Theorem 1. Consider the closed-loop system (31)-(33). Assume that the position command rB=G G is such that rB=G G(t) = 0;8t2 R0. In this context, if the initial state ( rB=G G(0);_rB=G G(0))2 L p(rB=G G(0);_rB=G G(0)), then rB=G Gconverges exponentially to rB=G Gand (rB=G G(t);_rB=G G(t))2L p (rB=G G(t);_rB=G G(t));8t2 R>0. Proof. Under the condition that p2F or, equivalently, ( ~r;_~r)2Lp, the closed-loop transla- tional dynamics represented by equations (31)-(33) can be re-modeled in terms of ~rand its time derivatives, resulting in ~r+K2_~r+K1~r=rB=G G: (48) Since, by assumption, rB=G G(t) = 0;8t2 R0and the gain matrices K3andK4are positive-de nite and diagonal, from the linear control theory, one see that the state ( ~r;_~r) = (031;031) of system (48), without considering the control saturation, is globally exponentially stable. This implies that, for any 2 R>0, if (~r(0);_~r(0))2B, then ( ~r(t);_~r(t))2B;8t2 R>0, which means that Bis an invariant set of system (48). In particular, consider the ball L pde ned in (47). Since, by construction, L pLp, it holds that 19 p2F, from the above argumentation, one can see that L pis an invariant set of system (31)-(33), which, in this case, behaves as (48) (with rB=G G(t) = 0;8t2 R0). Therefore, if ( ~r(0);_~r(0))2L p, then ( ~r(t);_~r(t))2L p;8t2 R>0and~rconverges to 031exponentially. From the de nition of ~r, this is equivalent to say that ( rB=G G(t);_rB=G G(t))2L p(rB=G G(t);_rB=G G(t));8t2 R0andrB=G Gconverges torB=G Gexponentially, thus completing the proof.  3.5.2. Attitude Control Denote the admissible-torque set T, T2 R3:TmaxTTmax (49) and de ne the corresponding viable set La,n "; B=G B 2 R6: a2To ; (50) as well as the maximal 0-centered Euclidean ball inside La L a,max Bs:t:BLa: (51) From the assumption that DB=G(t) is constant (time-scale separation), it follows that  B=G B(t) = 0. In this case, the attitude error kinematics can be described by (see [23], p. 73) _"=A(") B=G B; (52) where A("),2 66664c"3=c"2s"3=c"20 s"3 c"3 0 c"3s"2=c"2s"3s"2=c"213 77775; (53) with"i;i= 1;2;3;denoting the components of ". In this section, we assume that "has small components. Such small-angle approximation is 20 reasonable, since "represents an attitude error rather than the total attitude itself and, from the time-scale separation assumption, DB=Gis supposed to rapidly converge to DB=G(which means that"rapidly converges to zero). In this case, equation (52) can be approximated by _"= B=G B: (54) The following theorem gives a stability result for the closed-loop system (52), (37)-(38). Theorem 2. Consider the closed-loop system (52), (37)-(38). Assume that DB=G(t) is constant at anyt2 R0(time-scale separation) and that "is so small that its kinematics can be described by (54). In this context, if the initial state ( "(0); B=G B(0))2L a, then DB=Gconverges exponentially toDB=Gand ("(t); B=G B(t))2L a;8t2 R>0. Proof. Under the condition that a2Tor, equivalently, ( "; B=G B)2La, and considering the small- angle kinematics described by (54), the closed-loop rotational dynamics represented by equations (54),(37)-(38) can be re-modeled in terms of "and its time derivatives, resulting in "+K4_"+K3"=0: (55) Since the gain matrices K3andK4are positive-de nite and diagonal, from the linear control theory, one see that the state ( ";_") = (031;031) of system (55), without considering the control saturation, is globally exponentially stable. This implies that, for any 2 R>0, if ("(0);_"(0))2B, then ("(t);_"(t))2B;8t2 R>0, which means that Bis an invariant set of system (55). In particular, consider the ball L ade ned in (51). Since, by construction, L aLa, it holds that a2T, from the above argumentation and assumptions, one can see that L ais an invariant set of system (54), (37)-(38), which, in this case, behaves as (55). Therefore, if ( "(0);_"(0))2L a, then ("(t);_"(t))2L a;8t2 R>0and"converges to 031exponentially. From the de nition of "and considering (54), this is equivalent to say that ( DB=G; B=G B) converges to ( DB=G;031) exponentially, thus completing the proof.  21 4. Nominal Simulation As a rst step to evaluate the proposed method, this section presents a numerical simulation of the the proposed ight control system in nominal conditions. It is worth noting at this point that both the aerostatic lift ~Fband the restoring torque ~Tb depend on the value of airandhelium . These densities in turn depend on the local temperature and pressure, as expressed by the Ideal Gas Law: q=p RqT; (56) where \q" can be replaced by either \air" or \helium", pis the local pressure in Pa, Tis the local temperature in K, and Rqis the speci c gas constant ( Rair= 286:9 J/Kg/K and Rhelium = 2077 J/Kg/K). In particular, in this section, the surrounding air and contained gas are suposed to have the same temperature and pressure; their values are 20oC and 101325 Pa, respectively. Moreover, the air is assumed to be dry. Under such conditions, the air and helium densities are air1:204 kg/m3andhelium0:1664 kg/m3. These values are used here to simulate the ground truth as well as to compute the control laws. Subsection 4.1 describes the simulation and shows the adopted nominal parameters, while Subsection 4.2 presents and analyzes the simulation results. 4.1. Plant and Controller Parameters For simulating the overall closed-loop ight control system represented in Figure 4, we use the models formulated in Section 2 as well as the control laws and control allocation proposed in Section 3. The simulator is implemented in MATLAB, using the forth order Runge-Kutta solver, with an integration step of 0.001 s. Table 1 shows the nominal values of the hexa-rotor airship parameters. On the other hand, Table 2 presents the adopted controller parameters. The controller gains are tuned by trial and error, taking into account their proportional or derivative e ect and considering the time-scale separation assumption as well. In this paper, for obtaining simulation data that are consistent with a typical operation of MAVs, the proposed ight control system is commanded to follow a waypoint-based position 22 Table 1: Physical parameters of the hexa-rotor airship. Description Symbol Value Force coecient kf 1:2838105N/(rad/s)2 Torque coecient k 3:0811107Nm/(rad/s)2 Maximum rotor speed  !max 906:66 rad/s Motor speed coecient k! 1 Motor time constant ! 0:01 s Arm length l 1 m Volume of the balloon V 5:3 m3 CB-CM displacement d 0:85 m Semi-major axis a 1.25 m Semi-minor axis b 0.8 m Total inertia matrix J diag(2:0633;2:0651;1:9556) kgm2 Moment of inertia of the rotors Jr0:001 kgm2 Total empty mass mt9:392 Kg Table 2: Parameters of the attitude and position control laws. Description Symbol Value Proportional gain of the position controller K1 diag(0:5;0:2;0:7) Derivative gain of the position controller K2 diag(2;1;3) Proportional gain of the attitude controller K3 diag(20;50;1) Derivative gain of the attitude controller K4 diag(10;20;1) Maximum torque command Tmax[16:3 14:1 0:58]TNm Minimum force command Fmin[5:85:8 2:7]TN Maximum force command Fmax[5:8 5:8 54:6]TN Maximum inclination angle 'max12 degrees trajectory. In this trajectory, the waypoints are connected by straight lines with length of 5 m and a constant desired velocity of  v= 0:5 m/s. Moreover, the heading angle command  is set to zero. 4.2. Nominal Simulation Results Figures 6{8 are the results of the nominal simulation using the parameters of Tables 1{2. Figure 6 shows the e ective position and the corresponding position command. In the ramp parts of the component trajectories, one can verify a steady-state error of about 2 m. After nishing the ramps, there are no overshoots and the accommodation times (within 5 cm of the nal value) are between 10 s and 14 s. In the at parts of the trajectory, the steady-state errors are zero. The large errors during the ramps and the long accommodation times re ect the vehicle's inherent low acceleration 23 capability, which can be attributed to both the restricted inclination bound 'maxand the large inertia parameters. Figure 6: Nominal performance of the position control. It shows the controlled position components together with the respective command. Figure 7 shows the e ective and commanded attitude along the time using Euler angles in the 1-2-3 sequence. The attitude pitch and roll commands deviate from zero at the waypoints and the corresponding e ective angles follow them with an apparently fast but attenuated response. In the yaw angle, one can see a very small and short transient starting at 20 s; it is due to the gyroscopic coupling of the rotational dynamics, which becomes apparent when there exist a simultaneous motion in roll and pitch. Figure 7: Nominal performance of the attitude control. It shows the controlled attitude variables together with the respective command. Finally, Figure 8 presents the force magnitude command Fccomputed by the position controller and Figure 9 shows the components of the torque command Tc Bgenerated by the attitude controller. 24 There are transients at the instants the vehicle is commanded to ascend or to stop ascending, after which, Fcconverges and stays around 42 :1 N. This force value is equivalent to the total weight of the vehicle, payload, and helium, minus the aerostatic lift. The force command Fckeeps far from its maximum bound, even during the ascending acceleration. On the other hand, the components of the torque command Tc Bdeviate from zero at the waypoints, but stay inside their bounds Tmax with a large margin. Figure 8: Force magnitude command provided by the position controller in the nominal simulation. Figure 9: Torque command provided by the attitude controller in the nominal simulation. 25 5. Uncertainty Quanti cation The present section conducts an uncertainty quanti cation analysis of the ight control system under consideration using a parametric probabilistic approach [24, 25]. In this analysis, ambient the temperature and pressure are modeled as independent random variables and, based on such probabilistic models, for each ight realization, one new outcome is generated for both variables. For simulating the plant, the surrounding air and contained helium are assumed to have the same (varying) temperature and pressure. On the other hand, in the implementation of the control laws, their densities are maintained constant (and equal the respective nominal values; see Section 4). 5.1. Probabilistic Model Let (; ; P) be a probability space, where  is the sample space, is a-algebra over , and P: ![0;1] a probability measure. In this stochastic framework, the local temperature Tland pressurePlare modeled by (independent) random variables X: ! Rand Y: ! R, respectively. We assume that the joint probability distribution of these random variables, P XY(dxdy ), admit a probability density function (PDF) ( x;y)7!p XY(x;y), with respect to dxdy . Similarly, the marginal distributions P X(dx) andP Y(dy) are assumed as having marginal PDFs x7!p X(x), with respect todx, andy7!p Y(y), with respect to dy. The maximum entropy principle [26] is adopted here as a way to specify the aforementioned distributions. This approach provides the least biased distributions in a scenario with little (or no) experimental information on a particular random variable [24]. In the present investigation, the random local temperature Xand local pressure Yare both completely unknown beforehand. However, one should better specify for them resonable operational ranges, from which it is possible to immediately de ne the underlying supports. In this way, we assume from now on that Xand Y take values in known positive nite intervals [ x1;x2](0;+1) and [y1;y2](0;+1), respectively. Therefore, the maximum entropy principle is formulated here as the maximization of the Shannon entropy S(p XY) =Z Z R2lnp XY(x;y)p XY(x;y)dxdy; (57) 26 subject to the normalization constraint (which brings information about the supports) Zy2 y1Zx2 x1p XY(x;y)dxdy = 1: (58) By using the Lagrange multiplier method, one can solve the above optimization problem to obtain the following PDF p XY(x;y) =p X(x)p Y(y); (59) with p X(x) = 1[x1;x2](x)1 x2x1andp Y(y) = 1[y1;y2](y)1 y2y1; (60) where 1X(x) denotes the indicator function of the set X. Note that the above marginal PDFs correspond to uniform distributions on the intervals [ x1;x2] and [y1;y2], respectively. Informally, the above result makes all sense, because since the minimum and maximum values are the only known informaion about the random parameters, the most reasonable choice for their distributions is the one that assigns equal weight to all their possible values. Due to the randomness of Xand Y, the air density airand the helium density helium are also random variables [see equation (56)]. Further, since the aerostatic lift [see equation (7)] and the restoring torque [see equation (9)] depend on airandhelium , then the translational (14) and rotational (17) dynamic models are rather di erential equations with random parameters and, consequently, their variables are stochastic processes. In particular, we choose here the minimum temperature x1= 0oC, the maximum temperature x2= 40oC, the minimum pressure y1= 0:7739 atm (at 2000 m and 0oC), and the maximum pressure y2= 1 atm. 5.2. Propagation of Uncertainties To compute the propagation of uncertainties of the random parameters Xand Ythrough the control system dynamics, we employ the Monte Carlo (MC) method [27]. It consists in a stochastic solver that, based on nsrealizations of the random parameters Xand Y, generates nsdi erent real- izations of the nominal closed-loop dynamics and computes the sample statistics of the respective system responses. These statistics are then used to approximate the solution (or response) of the original stochastic system. 27 Figure 10: Convergence metrics for choosing the number of Monte Carlo realizations. In order to guarantee the convergence of the MC solution to the true stochastic behavior of the system, a sucient number nsof realizations must be considered. The reference [28] shows that such a convergence occurs in the mean-square sense if the following metrics converge in the point-wise sense: p= 1 nsnsX n=1Ztf t=0 jjrB=G G(t;n)jj2 dt!1=2 ; (61) a= 1 nsnsX n=1Ztf t=0 jj B=G(t;n)jj2 dt!1=2 ; (62) wherendenotes an arbitrary MC realization, tdenotes the continuous time, tfis the nal simula- tion instant, rB=G G(t;n) is thenth realization of the vehicle's position at instant t, and B=G(t;n) is thenth realization of the vehicle's attitude (in Euler angles 1-2-3) at instant t. Figure 10 shows a plot of the position (61) and attitude (62) convergence metrics versus ns. One can see that after 50 realizations, both metrics show to converge. Therefore, in the problem under consideration, we considerns= 100, which is sucient for the convergence of the MC solution. Figure 11 shows the position statistics (mean and con dence region). One can see that, as a consequence of the randomness in the temperature and pressure conditions, there is a very small variation in the r1andr2performance (less the 5 cm), and a considerable variation in the r3 performance (about 70 cm). Figure 12 shows the statistics (mean and con dence region) of the attitude. The largest variations, which are very small (less than 0 :3 degree), occur during the 28 transients. Figure 11: Position statistics (mean and con dence region) estimated by Monte Carlo method. Figure 12: Attitude control error statistics (mean and con dence region) estimated by Monte Carlo method. Figures 13 and 14 show the statistics (mean and con dence region) of the force magnitude command and torque command, respectively. One can observe a small random variation in the torque command components. On the other hand, there is a large variation in the force magnitude command; this is a direct consequence of the aerostatic lift variation (see Figure 15), which is proportional to the variation of the air density (see equation (7)). The extreme values of both commands keep inside their bounds with a good margin. Figure 13: Torque command statistics (mean and con dence region) estimated by Monte Carlo method. Figure 14: Force magnitude command statistics of the Monte Carlo simulation. 29 Figure 15: Aerostatic lift statistics of the Monte Carlo simulation. 6. Concluding Remarks The paper proposed a ight control system for a new micro aerial vehicle (MAV) resulting from the combination of a hexa-rotor airframe with a balloon lled with helium. This MAV is a straightforward attempt to extend the ight duration and load capacity of conventional multi-rotor vehicles by means of an aerostatic lift. On one hand, such a combination results in a very simple vehicle. On the other hand, even for indoor ight, its control system must be carefully designed to overcome the e ects of a restoring torque (which does not occur in conventional multi-rotor vehicles), while being robust to variations in the local temperature and pressure operation con- ditions. The proposed ight control system is evaluated on the basis of both deterministic and stochastic simulations. In particular, the deterministic simulation shows that in nominal condi- tions, it is possible to control the vehicle's position to follow a desired waypoint-based trajectory with reference speed of 0 :5 m/s. The obtained performance is sucient accurate and fast for many MAV applications. On the other hand, in the stochastic simulation, the local temperature and pressure are modeled by uniformly-distributed random variables and the resulting performance is rigorously evaluated by Monte Carlo (MC) simulation. The variability veri ed in the position control performance can be considered small enough to dispense a special design of the control laws in order to suppress it. However, the variability shown in the attitude control error motivates a special robust control method for a better performance of the system in a stochastic sense. Acknowledgment The authors would like to thank all those people who indirectly contributed or motivated us to conduct the present investigation. In special, we are grateful to Mr. Luciano S. Araujo (from Elio Tecnologia), who designed and built the rst prototype of the balloon-hexacopter and donated it to the Laboratorio de Robotica Aerea do ITA. The authors are thankful for the nancial support provided by S~ ao Paulo Research Foundation (FAPESP) grant number XXX, Carlos Chagas Filho Research Foundation of Rio de Janeiro State (FAPERJ) under grants E-26/010.002.178/2015 and E-26/010.000.805/2018. 30