UNIVERSIDADE ESTADUAL PAULISTA “JÚLIO DE MESQUITA FILHO” FACULDADE DE ENGENHARIA DE ILHA SOLTEIRA DEPARTAMENTO DE ENGENHARIA MECÂNICA Renan Cavenaghi Silva Modeling and Control of Multirotor UAVs Formation Ilha Solteira 2023 UNIVERSIDADE ESTADUAL PAULISTA “JÚLIO DE MESQUITA FILHO” Campus de Ilha Solteira Modeling and Control of Multirotor UAVs Formation Renan Cavenaghi Silva Advisor: Prof. Dr. Douglas Domingues Bueno Co-advisor: Prof. Dr. Rodrigo Borges Santos Dissertation presented to the School of Engineer- ing of Ilha Solteira - UNESP - as part of the requirements to obtain the title of Master in Me- chanical Engineering. Knowledge Area: Solid Mechanics Ilha Solteira - SP February/2023 Cavenaghi SilvaModeling and Control of Multirotor UAVs FormationIlha Solteira2023 67 Sim Dissertação (mestrado)Engenharia MecânicaMecânica dos SólidosNão . FICHA CATALOGRÁFICA Desenvolvido pelo Serviço Técnico de Biblioteca e Documentação Silva, Renan Cavenaghi. Modeling and control of multirotor UAVs formation / Renan Cavenaghi Silva. -- Ilha Solteira: [s.n.], 2023 65 f. : il. Dissertação (mestrado) - Universidade Estadual Paulista. Faculdade de Engenharia de Ilha Solteira. Área de conhecimento: Mecânica dos Sólidos, 2023 Orientador: Douglas Domingues Bueno Coorientador: Rodrigo Borges Santos Inclui bibliografia 1. UAV. 2. Dynamics. 3. Control. 4. Formation flight. S586m ACKNOWLEDGMENTS The authors thank the São Paulo State University (UNESP), the team of Intelligent Mate- rials and Structures (GMSINT), and the National Council for Scientific and Technological Development (CNPq) for the financial support through the grant nº152165/2021-5. This study was financed in part by the Coordenação de Aperfeiçoamento de Pessoal de Nı́vel Superior - Brasil (CAPES) - Finance Code 001. ABSTRACT The Unmanned Aerial Vehicles (UAV) have been extensively researched due to their under- lying characteristics that suit them for different applications such as border surveillance, object recognition, crop spraying and aerial transportation. Despite the enthusiasm from the general community, practical tasks are restricted since the specialized equipment that the UAV have to transport in a mission have limited weight. In this sense, a possible alternative consists of using multiple UAVs in a formation flight to collaboratively transport a heavy-payload. Then, in this dissertation, the dynamic model from an UAV formation is obtained using an approach to parameterize the equation of motion in terms of the number of vehicles, such that different formation configurations can be considered. Moreover, the presented controller involves a leader-follower strategy with virtual-constraints for the vehicles to keep each particular formation. The results show that the employed control strategy is feasible for transporting a heavy-payload despite showing a persistent error due to the unmodeled dynamic from the suspended payload in the formation controller. Keywords: UAV. Dynamics. Control. Formation Flight. RESUMO Os Véıculos Aéreos não Tripulados (VANT) têm sido amplamente pesquisados devido às suas caracteŕısticas que permitem seu uso em aplicações como monitoramento de fronteiras, reconhecimento de objetos, pulverização de defensivos agŕıcolas e transporte aéreo de cargas. Apesar do entusiasmo da comunidade, as aplicações são restritas no sentido de que o equipamento especializado para cada missão que o VANT deve transportar tem peso limitado. Nesse sentido, uma posśıvel alternativa consiste em utilizar múltiplos véıculos em um voo em formação para colaborativamente transportar uma carga pesada. Assim, neste trabalho apresenta-se o modelo dinâmico de uma formação de VANTs usando uma abordagem que permite parameterizar as equações do movimento no número de véıculos, de modo que diferentes configurações de formações podem ser consideradas. Ainda, apresenta- se um controlador implementando a estratégia ĺıder-seguidor com v́ınculos virtuais para definir e manter os seguidores em suas respectivas posições na formação. Os resultados mostram que a estratégia empregada é viável para o transporte de cargas pesadas, apesar de apresentar um erro persistente decorrente da dinâmica da carga desconsiderada no controlador da formação. Palavras-chave: VANT. Dinâmica. Controle. Voo em Formação. List of Figures 1.1 Illustration of the steps of a typical multi-lift UAV mission. . . . . . . . . . 14 2.1 Representation of a generic rigic body . . . . . . . . . . . . . . . . . . . . . 18 2.2 Sequential transformation from the coordinate system through the yaw (a), pitch (b), and roll (c). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 2.3 Quadrotor configuration of a multirotor UAV. . . . . . . . . . . . . . . . . 23 2.4 Generic trajectories connecting the configurations at the instants t1 and t2. The system takes the trajectory that satisfies the Euler-Lagrange equation indicated by the solid line. . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 2.5 Representation of the pendulum. . . . . . . . . . . . . . . . . . . . . . . . 27 2.6 Single-lift configuration with an UAV transporting a suspended payload. . 28 2.7 Multi-lift UAV formation transporting the suspended payload. . . . . . . . 30 2.8 Representation from the normalized 1-cosine gust equation at (x0 y0). . . . 33 2.9 Illustration from the constraint drift and the stabilization method. . . . . . 34 2.10 Schematic representation of the equilibrium configuration for the single-lift. 38 2.11 Schematic representation of the equilibrium configuration for the twin-lift. 38 2.12 Eigenvalues from the twin-lift configuration. . . . . . . . . . . . . . . . . . 39 3.1 Diagram illustrating the process for the obtaining the formation trajectories. 43 3.2 Schematic representation of the “V” formation. . . . . . . . . . . . . . . . 45 4.1 Trajectories from the single-lift system. . . . . . . . . . . . . . . . . . . . . 48 4.2 Constraint equation from the single-lift system . . . . . . . . . . . . . . . . 49 4.3 Position error from the UAV and payload for αf = 100. . . . . . . . . . . . 50 4.4 Constraints over time from the twin-lift system. . . . . . . . . . . . . . . . 52 4.5 Top-view from the twin-lift simulation. . . . . . . . . . . . . . . . . . . . . 53 4.6 Attitude from the UAVs considering the twin-lift configuration. . . . . . . 53 4.7 Virtual constraints from the twin-lift formation control. . . . . . . . . . . . 54 4.8 Thrust force required from the UAVs in the twin-lift configuration. . . . . 54 4.9 UAVs trajectories from the twin-lift configuration. . . . . . . . . . . . . . . 55 4.10 Constraints from the multi-lift system. . . . . . . . . . . . . . . . . . . . . 55 4.11 Virtual constraints from the multi-lift formation control. . . . . . . . . . . 56 4.12 Thrust forces from the UAVs in the multi-lift configuration. . . . . . . . . 56 4.13 Trajectories (Top-view) for the multi-lift configuration. . . . . . . . . . . . 57 4.14 Trajectories from the multi-lift system. . . . . . . . . . . . . . . . . . . . . 57 4.15 Comparison between the payload position when considering the twin-lift and multi-lift configurations. . . . . . . . . . . . . . . . . . . . . . . . . . . 58 4.16 Comparison between the virtual constraints when considering the twin-lift and multi-lift configurations. . . . . . . . . . . . . . . . . . . . . . . . . . . 58 4.17 Comparison from the rotor speed from the leader in the twin-lift and multi-lift formations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 List of Tables 2.1 Multirrotors configurations and the associated control allocation matrix. . 22 2.3 Comparative of the kinematic constraints from the pendulum geometry using Cartesian and polar coordinates. . . . . . . . . . . . . . . . . . . . . 27 2.4 Dimensions of the matrices and vectors for the multi-lift equation of motion. 32 3.1 Different formation geometries and their corresponding virtual constraint equations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 4.1 Single-lift system parameters. . . . . . . . . . . . . . . . . . . . . . . . . . 47 4.2 Aerial transportation mission parameters. . . . . . . . . . . . . . . . . . . 51 Contents List of Figures 8 List of Tables 10 1 INTRODUCTION 12 1.1 PROBLEM DESCRIPTION . . . . . . . . . . . . . . . . . . . . . . . . . . 13 1.2 OBJECTIVE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 1.3 PUBLISHED WORKS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 1.4 OUTLINE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2 MODELING THE UAV FORMATION 17 2.1 NOMENCLATURE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 2.2 MODELING PRELIMINARIES . . . . . . . . . . . . . . . . . . . . . . . . 17 2.2.1 Rotation Matrix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 2.2.2 Euler Angles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 2.3 MULTIROTOR UAV . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 2.4 MODEL DERIVATION . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 2.4.1 Translational Dynamics . . . . . . . . . . . . . . . . . . . . . . . . 23 2.4.2 Attitude Dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.4.3 Dynamic Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 2.5 EULER-LAGRANGE EQUATIONS FOR CONSTRAINED SYSTEMS . . 26 2.6 SINGLE-LIFT DYNAMIC . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 2.7 MULTI-LIFT FORMATION DYNAMIC . . . . . . . . . . . . . . . . . . . 29 2.7.1 Aerodynamic Drag Model . . . . . . . . . . . . . . . . . . . . . . . 32 2.7.2 Gust Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 2.7.3 Index Reduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 2.7.4 Constraint Stabilization . . . . . . . . . . . . . . . . . . . . . . . . 33 2.7.5 Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 2.7.6 Null Space Formulation . . . . . . . . . . . . . . . . . . . . . . . . 36 2.8 FORMATION MODAL ANALYSIS . . . . . . . . . . . . . . . . . . . . . . 36 2.8.1 Single-Lift Modal . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 2.8.2 Twin-Lift Modal . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 3 CONTROL OF UAV FORMATIONS 40 3.1 RELATED WORK . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 3.1.1 Propulsion Saturation . . . . . . . . . . . . . . . . . . . . . . . . . 41 3.2 FORMATION CONTROL . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 3.2.1 “V” formation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 4 RESULTS AND DISCUSSION 47 4.1 SINGLE-LIFT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 4.2 MISSION DEFINITION . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 4.3 TWIN-LIFT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 4.4 MULTI-LIFT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 4.5 RESULTS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 5 FINAL REMARKS 60 5.1 CONCLUSIONS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 5.2 SUGGESTIONS FOR FUTURE WORK . . . . . . . . . . . . . . . . . . . 61 Chapter 1 INTRODUCTION Unmanned Aerial Vehicles (UAV) have attracted attention from both researchers and the community due to their potential to perform different applications. The advantages include the lower cost of operation and acquisition, relatively simple training, feasible to use in dangerous environments, updated regulation from certification agents, besides the autonomous operation. In Brazil, UAVs have been extensively employed in the agriculture, with emphasis to modernize the analysis of crop health using aerial photography (e.g., with the Normalized Difference Vegetation Index). Then, farmers can act in localized regions and reduce the use of water and chemicals. The logistic industry can also benefit from UAVs in regions with higher population density, low infrastructure, or difficult terrain, where UAVs can surpass obstacles for ground vehicles and allow fast transportation of critical supplies such as medical products. The use of UAVs in challenging missions, such as those ones involving heavy payload transportation increases the UAV complexity due to the increasing maximum take-off weight (MTOW). Then, their costs increase and, in practice, there is a very limited number of these aerial vehicles for this type of application. An interesting solution that is also the main focus of this research consists of using a formation to transport heavy payloads, where multiple UAVs collaboratively lift a payload heavier than the MTOW from each single vehicle. This chapter is divided as follows: Section 1.1 describes the context and technical challenges that motivates this work. The objectives are outlined in Section 1.2. The published works are presented in Section 1.3. A description from the discussed topics are presented in Section 1.4. 13 1.1 PROBLEM DESCRIPTION Formations are solutions that benefits from the collective mainly to perform tasks that are difficult, time-consuming, and even impossible for single individuals. Occurrences of formations in nature can be observed from migrating birds and transportation of heavy weights by ants. Portugal et al. (2014) show that the northern bald ibises, Geronticus eremita, position themselves in an aerodynamically efficient “V” formation to exploit the upwash vortices from the birds ahead. This biological behavior inspires research in commercial aviation to reduce the carbon-footprint by smart scheduling flight routes to allow aircrafts to fly in formation (DURANGO; LAWSON; SHAHNEH, 2016). Other species such as the Humpback whales, Megaptera novaeangliae, can swim in formations to manage their prey to form clusters, then the whales efficiently eat large amounts of food (HAIN et al., 1981). Cooperative transportation using UAVs is a strategy under development by research groups with the perspective to transport heavy-payloads (RASTGOFTAR; ATKINS, 2019). The main idea consists of distributing the payload weight such that the lift force of each one is in a feasible domain. A cable is used for constraining the payload to each UAV, and it transfers a pulling force which is used to guide the payload through the desired path. A typical mission performed by UAVs has common features with manned vehicles (RAYMER, 1989), indicated in Fig. 1.1 with a sequence of well-defined steps. � Phase 1 - take-off: consists of the UAVs increasing lift to compensate the weight; � Phase 2 - cruise: consists of moving the vehicles to the designated mission area. This motion need to occur in an power efficiently manner to increase the system endurance; � Phase 3 - mission: consist of operating the equipment on the desired areas to carry out a particular task. It can be searching a specific area, following a target, and hovering at certain spots; � Phase 4 - returning cruise: the UAVs returns towards the landing area; � Phase 5 - landing: the UAVs lands. 14 1 2 4 5 Figure 1.1: Illustration of the steps of a typical multi-lift UAV mission. Aerial transportation with a single vehicle, denoted as single-lift system, with a cable suspension allow one to transport cargo that do not fit inside its cargo compartment (RAST- GOFTAR; ATKINS, 2019). However, this transporting configuration also poses additional challenges due to the introduction of pendulum modes (SHI; WU; CHOU, 2018). By suspending the payload with two UAVs, i.e., twin-lift, the pendulum motion occurs per- pendicular to the plane formed by the payload and attachments points (BERRIOS et al., 2014). In contrast, by using three or more UAVs, i.e., multi-lift formation, the pendulum motion from the suspended payload is reduced. Despite the dynamic advantages, remotely piloting an UAV formation poses additional challenges due to closeness from the vehicles requiring coordination to avoid collisions. The twin-lift and multi-lift formations show a change from the equilibrium configuration which requires a change of the attitude from each vehicle to compensate the additional weight and maintain a separation. This effect increases the complexity of manually piloting a multi-lift formation (RAZ; ROSEN, 2005). However, UAVs can autonomously perform multi-lift missions using feedback control. The autonomous operation from an UAV requires the use of onboard controllers to manage the mission path and compute the adequate set of control inputs to guide the vehicle towards the desired trajectory. Modeling the UAV is an essential step to design a suitable controller. Although modeling the dynamic from a single UAV is widely presented in the specialized literature, this step increases in complexity for the aerial transportation mission, mainly due to kinematic constraints that connects the UAV to the payload. For the single-lift system, the common approach consists of applying a coordinate transformation to yield a set of independent generalized coordinates that maps directly to the constraint hypersurface in the original coordinate configuration (SHI; WU; CHOU, 2018; GUERRERO-SÁNCHEZ et al., 2017). Nonetheless, this approach introduces singular configurations and there is an important difficulty to scale it for an arbitrary number of vehicles. An alternative approach consists in obtaining the dynamics in the original coordinate space using the Euler-Lagrange formulation for constrained systems (SALETAN; CROMER, 1970). 15 Expressing the dynamic from the multi-lift system in the original coordinate space requires satisfying both the obtained differential equations, modeled from the Euler- Lagrange formulation for constrained systems (SALETAN; CROMER, 1970), and the algebraic equations from the imposed constraints, resulting in a set of differential and algebraic equations (DAE). To solve them for the trajectories in time, the resulting system of DAE are transformed into a set of ODE through the time-differentiation from the algebraic constraints (ASCHER et al., 1995), and this process is known in the literature by index-reduction. This procedure, however, introduces numerical instabilities because the constraint equations are solved only on acceleration level, such that the error is accumulated in the velocity and position levels. This phenomenon is widely recognized in the literature and it is named constraint-drift. This problem is solved by using the technique proposed by Baumgarte (1972), which is interpreted by Bisgaard, Bendtsen, and Cour-Harbo (2009) as artificially introducing springs and dampers in parallel with each one of the rods. Then, the constraints equations are approximately satisfied during the solution. Based on the dynamic model describing the multi-lift formation, it is possible to evaluate different strategies to control the formation. There is an extensive amount of strategies for controlling multi-agent systems such as the Potential Field Method (PAUL; KROGSTAD; GRAVDAHL, 2008; ZE-SU; JIE; JIAN, 2012) and Leader-Follower Architecture (YUN et al., 2010; HE et al., 2018; LIANG; DONG; ZHAO, 2020). However, there is narrow discussion in the literature regarding the effects of a formation controller for multi-lift system. 1.2 OBJECTIVE The main objective of this research is to investigate the dynamic and control of an UAV formation transporting a heavy-payload. 1.3 PUBLISHED WORKS The work entitled “Disturbance Observer of an UAV with a Suspended Payload” was presented in the 26th International Congress of Mechanical Engineering. The work entitled “Quaternion-Based Attitude Control of a Multirotor UAV ” was presented in the workshop to celebrate the 25th anniversary from the Graduate Program of Mechanical Engineering from the School of Engineering of Ilha Solteira (PPGEM). The following items summarize the contributions from this work: • A dynamic model is parameterized in the number of vehicles, which allows one to obtain the EOM from the single-lift, twin-lift, and multi-lift configurations; • A controller based on the virtual constraints is developed to maintain a formation during all flight and the implication from the suspended payload on the “V” formation 16 geometry during the multi-lift transportation mission. 1.4 OUTLINE The remaining of this work is divided in chapters. Chapter 2 presents the modeling procedure from first principles, where the Euler-Lagrange formalism for constrained systems is applied to describe the dynamics from a multi-lift UAV formation. Chapter 3 presents a literature review on techniques to control UAVs. The developed method to control a UAV formation is also presented. Chapter 4 presents the results for a single-lift, twin-lift and the multi-lift UAV formation. Chapter 5 presents the final conclusions from this work along suggestions for future work. Chapter 2 MODELING THE UAV FORMATION This chapter presents the dynamic modeling of the multi-lift UAV formation. The requirements are introduced, then, the model from the UAV formation is derived and technical challenges are discussed. This chapter is divided into sections described by the following chapter content. Section 2.1 defines the mathematical symbols. Section 2.2 introduces the parameters used in the modeling development throughout this text. Section 2.3 presents particularities from the construction of multirotors UAVs and the major differences from other configurations. Section 2.4 introduces the dynamic modeling of a single multirotor UAV using the Euler- Lagrange formulation. Additional considerations for system with constraints are presented in Section 2.5. The modeling of the dynamic from the multi-lift system is presented in Section 2.7. Section 2.8 presents the procedure to analyze the multi-lift system dynamic without solving the EOM. 2.1 NOMENCLATURE For the remaining of this text, scalar variables are defined by a lower-case letter (e.g., a), vectors by a lower-case bold letter (e.g., a), and matrices by an upper-case bold letter (e.g.. A). The transpose is denoted by A>. The identity n× n matrix is defined by In×n, and the n×m zero matrix is defined by 0n×m. 2.2 MODELING PRELIMINARIES Rigid bodies are representations of objects whose dimensions are relevant for the problem. The underlying assumption in a rigid body is the distance between any two of its particles is constant (NETO, 2013). Figure 2.1 illustrates a generic rigid body and an arbitrary particle p. Consider an inertial frame of reference denoted by the tuple I : (O, î, ĵ, k̂), where O denotes an arbitrary origin, and î = {1 0 0}>, ĵ = {0 1 0}> and k̂ = {0 0 1}> the orthonormal basis that spans the Euclidean space R3. Moreover, there is a body frame 18 O î ĵ k̂ Ob îb k̂b r Bρ p ĵb rρ Figure 2.1: Representation of a generic rigid body, where O is the coordinate system origin; î, ĵ and k̂ are the basis from the inertial frame of reference and îb, ĵb and k̂b the basis from the body frame of reference, r and rρ are the position from the center of mass and the arbitrary particle respectively. of reference B : (Ob, îb, ĵb, k̂b) with its origin Ob attached to the body center of mass at a distance r(t) with respect to the inertial frame, and with the same orientation from the body. Then, the position from any particle p from the rigid body relative to the center of mass is constant in the body reference frame B, i.e., Bρ ≡ constant (2.1) The position from this arbitrary particle with respect to I can be parameterized by the sum of the position from the center of mass and its distance relative to the body center of mass expressed in the inertial frame by means of the rotation matrix R(t),i.e., rρ(t) = r(t) +R(t) Bρ. (2.2) 2.2.1 Rotation Matrix The rotation matrix R(t) is a transformation from the special orthogonal group of transfor- mations SO (3) : { R(t)|R(t) ∈ R3×3,R(t)>R(t) = R(t)R(t)> = I } that has the property to rotate a vector around its origin while preserving the norm and handness (FOSSEN, 1994). 2.2.2 Euler Angles The rotation matrix can be parameterized using Euler angles η(t) = {φ(t) θ(t) ψ(t)}> (FOS- SEN, 1994). By defining the rotations of yaw Tψ : ( î′ ĵ ′ k̂′ ) −→ ( î ĵ k̂ ) , pitch Tθ :( î′′ ĵ ′′ k̂′′ ) −→ ( î′ ĵ ′ k̂′ ) , and rolling Tφ : ( îb ĵb k̂b ) −→ ( î′′ ĵ ′′ k̂′′ ) , the components from the rotation matrix can be obtained by R = TψTθTφ. (2.3) Figure 2.2 illustrates each consecutive rotation of yaw, pitch and roll, which allows one to 19 O î ĵ ψ ψ î′ ĵ ′ k̂//k̂′ (a) O î′ ĵ ′//ĵ ′′ k̂′ θ θ î′′ k̂′′ (b) O φφ î′′//îb ĵ ′′ k̂′′ĵb k̂b (c) Figure 2.2: Sequential transformation from the coordinate system through the yaw (a), pitch (b), and roll (c). obtain the following equations î =cψ î ′ − sψĵ ′ ĵ =sψ î ′ + cψĵ ′ k̂ =k̂′ î′ =cθ î ′′ + sψk̂ ′′ ĵ ′ =ĵ ′′ k̂′ =− sθ î′′ + cθk̂ ′′ î′′ =îb ĵ ′′ =cφĵb − sφk̂b k̂′′ =sφĵb + cφk̂b (2.4) where c(·) and s(·) denotes cos (·) and sin (·) respectively. The transformation matrices T(·) are obtained from Eq. (2.4) and given by Eq. (2.5) Tφ = 1 0 0 0 cφ −sφ 0 sφ cφ  Tθ =  cθ 0 sθ 0 1 0 −sθ 0 cθ  Tψ = cψ −sψ 0 sψ cψ 0 0 0 1  (2.5) R = cψcθ −sψcφ + cψsθsφ sψsφ + cψcφsθ sψcθ cψcφ + sφsθsψ −cψsφ + sθsψcφ −sθ cθsφ cθcφ  . (2.6) In this case, the transformation R : ( îb ĵb k̂b ) −→ ( î ĵ k̂ ) maps a vector from the body reference frame to the inertial reference frame. The inverse transformation, i.e., that maps a vector from the inertial reference system to the body reference system is then R> = T>φ T > θ T > ψ . The angular velocity vector from the body Bω = Bω(t) = {ωx(t), ωy(t), ωz(t)}> is given in terms of the Euler angles, such that (FOSSEN, 1994): Bω =T>φT>θ  0 0 ψ̇ + T>φ  0 θ̇ 0 +  φ̇ 0 0  (2.7) 20 that results in Bω = Wηη̇ , with Wη denoting the transformation from the angular velocity in the inertial reference frame to the body reference frame (RAFFO; ORTEGA; RUBIO, 2010) Wη = 1 0 −sθ 0 cφ sφcθ 0 −sφ cφcθ . The Euler angles rate are obtained in terms of the body angular velocity through the inverse transformation defined by (WANG et al., 2016) η̇ = W−1 η Bω, (2.8) such that W−1 η = 1 sinφ tan θ cosφ tan θ 0 cosφ − sinφ 0 sinφ sec θ cosφ sec θ  . Euler angles are intuitive for understanding the geometric properties from rotations when considering isolated maneuvers. Nonetheless, once det(Wη) = cθ, the inverse transformation is singular for θ = ±π/2. 2.3 MULTIROTOR UAV A distinct feature from multirrotor UAVs is that the propulsion and control is performed by varying the angular velocity from the rotors composed of an electric motor and a fixed-pitch propeller. This characteristics allows one to achieve significant cost reduction in acquisition and maintenance since, in contrast with the helicopter configuration, the swash-plate mechanism is absent. However, the benefit of controlling the UAV by varying the rotor velocity is also challenging to scale for transporting heavier payloads, since this task usually requires bigger propellers which imply to increase the rotational inertia that slows the rotor dynamics and may render the vehicle uncontrollable (POUNDS; MAHONY, 2009; PORTER; SHIRINZADEH; CHOI, 2015). The thrust force and drag moment from the ith rotor can be reasonable modeled respectively by the following set of equations if small size propellers are considered. Bf i(ωr,i) = −kTω2 r,ik̂b (2.9) Bmi(ωr,i) = sign(ωr,i)kDω 2 r,ik̂b (2.10) with ωr,i denotes the angular velocity from the corresponding ith rotor, kT and kD are respectively linear coefficients of thrust and moment, commonly determined from the propeller performance data. Note the negative sign in Eq. (2.9) indicates the thrust pointing to the negative k̂b axis. Moreover, the drag moment direction depends if the rotor is rotating clockwise (CW) or counterclockwise (CCW). Thus it is necessary to arrange the 21 rotors to balance the yaw moment. This is one from the explanations of why multirotors are commonly built with an even number of rotors (e.g., quadrotors, hexarotors, and octorotors). In contrast, the yaw moment in a helicopter configuration is balanced using a tail-rotor. In the planar multirotor configuration, the force from the ith rotor is applied at Bri = {lrcξi − lrsξi 0}> in the body coordinate system, where ξi denotes the angle from the line connecting the UAV center of mass to the propeller rotation axis to the îb axis. Considering a symmetric arrangement of rotors, this angle is computed for the ith rotor by ξi = 2π nr i+ ξc, (2.11) where the two distinct + and × configurations can be obtained by setting the offset angle ξc to 0 or π nr , respectively. Then, the net lift and moment exerted on the UAV center of mass are obtained respectively by Bf = nr∑ i=1 Bf i (2.12) Bm = nr∑ i=1 Bri × Bf i + Bmi (2.13) The control input is defined in terms of the angular speed from each rotor and the control allocation matrix Kr (KOTARSKI et al., 2021), representing a map from the quadratic angular velocity, denoted by the vector ω2 r = { ω2 r,1 ω 2 r,2 . . . ω2 r,nr }> , of a given configuration to the control forces u = { uT uîb uĵb uk̂b }> , with uT denoting the k̂b component from the force vector in Eq. (2.12), and uîb , uĵb and uk̂b the moments on the UAV body axes given by Eq. (2.13) (PILJEK; KOTARSKI; KRZNAR, 2020). u = Krω 2 r , (2.14) Equation (2.14) allows generalizing the forces for any configuration of multirotor. Then, the rotors speed for a particular configuration are obtained through the following equation ω2 r = K−1r u. (2.15) Table 2.1 shows the matrix Kr for different multirrotor configurations. The different multirrotor configurations are named by Q, H, and O, where the letter Q denotes the quadrirrotor, H denotes the hexarrotor, and O the octorrotor configurations, respectively. The Q× configuration is widely used in practice since a camera aligned with the body coordinate system usually has its field of view clear from the rotor arm (RAO et al., 2022). Other configurations with more rotors, such as the hexacopter (H+) and octocopter (O+), are mainly used to provide more thrust and fail-safe features. Table 2.1: Multirrotors configurations and the associated control allocation matrix. Propellers indicated in red and blue represent a CW and CCW rotation respectively. representation conf. Kr ref. îb ĵb Q +  kT kT kT kT 0 −lrkT 0 lrkT lrkT 0 −lrkT 0 −kD kD −kD kD  (LUO; DU; YU, 2019) îb ĵb Q ×  kT kT kT kT √ 2 2 lrkT − √ 2 2 lrkT − √ 2 2 lrkT √ 2 2 lrkT √ 2 2 lrkT √ 2 2 lrkT − √ 2 2 lrkT − √ 2 2 lrkT −kD kD −kD kD  (RAO et al., 2022) îb ĵb H+  kT kT kT kT kT kT 0 √ 3 2 lrkT √ 3 2 lrkT 0 − √ 3 2 lrkT − √ 3 2 lrkT lrkT 1 2 lrkT −1 2 lrkT −lrkT −1 2 lrkT 1 2 lrkT −kD kD −kD kD −kD kD  îb ĵb O+  kT kT kT kT kT kT kT kT 0 √ 2 2 lrkT lrkT √ 2 2 lrkT 0 − √ 2 2 lrkT −lrkT − √ 2 2 lrkT −lrkT − √ 2 2 lrkT 0 √ 2 2 lrkT lrkT √ 2 2 lrkT 0 − √ 2 2 lrkT −kD kD −kD kD −kD kD −kD kD  23 2.4 MODEL DERIVATION The quadrotor is illustrated in Fig. 2.3. Two reference systems are defined according to Sec. 2.2 and the generalized coordinates that uniquely describes the UAV configuration is given by the vector q(t) = {r(t) η(t)}>, with r(t) = {x(t) y(t) z(t)}> denoting the position from the UAV center of mass and η(t) = {φ(t) θ(t) ψ(t)}> describing the attitude parameterized using Euler angles. O î ĵ k̂ îb ĵb k̂b f1 m1 f2 m2 f3 m3 f4 m4 ωz ωy ωx Figure 2.3: Quadrotor configuration of a multirotor UAV. The kinetic energy from the UAV is composed of the translational energy due the motion of the center of mass, and the kinetic energy due to body angular velocity, such that T = Tr + Tη. The potential energy due to the gravitational field is computed by U = Ur. Then, the Lagrangian is given by: L = T −U (GOLDSTEIN, 2002), which allows one to obtain the equation of motion for each generalized coordinate qi ∈ q through the Euler-Lagrange equation (LEMOS, 2007), i.e., d dt ( ∂L ∂qi ) − ∂L ∂qi = fi, (2.16) with fi denoting the generalized force vector. 2.4.1 Translational Dynamics The kinetic energy from the translating rigid body is computed by Tr = 1 2 ṙ>Mrṙ. (2.17) Due to the gravitational potential field, the UAV is also subjected to a potential energy Ur = −mgz. The partial derivative of Eq. (2.17) with respect to the velocity is obtained by ∂Tr ∂ṙ = mṙ. Differentiating with respect to the time, the resulting equation is d dt ( ∂Tr ∂ṙ ) = mr̈. (2.18) 24 The potential energy depends only on the position vector. Nonetheless, its partial derivative is computed by ∂Ur ∂r = −mgk̂. (2.19) Then, the equations for the translational dynamics in matrix form are given by Mrr̈ + fg = BruT + dr (2.20) with Mr = mI3×3 representing the translation inertia, fg = −mgk̂ is the gravitational force, and dr = {dx dy dz}> is the disturbance forces vector acting on the UAV. The components of matrix Br are the projections from the thrust acting on the UAV k̂b axis to the inertial reference frame, such that Br = cφsθcψ + sφsψ cφsθsψ − sφcψ cφcθ  . Note that for the multirotor there is only one control term to affect the three translational coordinates. This render the multirotor UAV as an underactuated vehicle (BRANDÃO; FILHO; CARELLI, 2013). 2.4.2 Attitude Dynamics The kinetic energy due to the rotational motion is computed by Tη = 1 2 Bω>J Bω. (2.21) Substituting Eq. (2.8) in the expression above, the kinetic energy is written in terms of the Euler angles Tη = 1 2 η̇>W> η JWηη̇. (2.22) The EOM from the attitude dynamic is decoupled from the translational dynamic, such that ∂Tη ∂η̇ = ∂ ∂η̇ ( 1 2 η̇>W> η JWηη̇ ) = W> η JWηη̇ = Mηη̇ (2.23) and the symmetric and positive definite inertia matrix Mη = W> η JWη is obtained, such that Mη =  Ixx 0 −Ixxsθ 0 Iyyc 2 φ + Izzs 2 φ sφcφcθ(Iyy − Izz) −Ixxsθ sφcφcθ(Iyy − Izz) Ixxs 2 θ + c2θ(Iyys 2 φ + Izzc 2 φ)  . (2.24) Differentiating Eq. (2.23) with respect to the time using the chain rule, the following equation is obtained d dt ( ∂Tη ∂η̇ ) = Ṁηη̇ +Mηη̈. (2.25) 25 The second term of the Euler-Lagrange equation corresponds to the partial derivative from the kinetic energy with respect to the coordinates. Then, ∂Tη ∂η = ∂ ∂η ( 1 2 η̇>W> η JWηη̇ ) = 1 2 η̇> ∂Mη ∂η η̇ (2.26) By grouping the angular velocity dependent terms, the Coriolis matrix is obtained as the following (BRANDÃO; FILHO; CARELLI, 2013) Cη = Ṁη − 1 2 η̇> ∂Mη ∂η . (2.27) The first term on the right-hand side of Eq. (2.27) corresponds to the derivative of a matrix with respect to a scalar, which consists of differentiating every component of Mη with respect to time. The second term consists of differentiating a matrix with respect to a vector, which is a tensor. The result is a 3× 3 matrix Cη with its elements cij given by c11 = 0 c12 = sφcφ(Iyy − Izz)θ̇ + s2φcθ(Iyy − Izz)ψ̇ − Ixxψ̇cθ c13 = sφcφc 2 θ(Izz − Iyy)ψ̇ + c2φcθ(Izz − Iyy)θ̇ − Ixxcθθ̇ c21 = Izzψ̇cθ + (Iyy − Ixx)(−θ̇sθcφ + ψ̇cθc 2 φ − ψ̇cθs2φ), c22 = −(Iyy − Ixx)sφcφφ̇ c22 = (Izz − Iyy)sφcφφ̇ c23 = ψ̇sθcθ(−Izz + Iyys 2 φ + Ixxc 2 φ) c31 = −(Ixxθ̇cθ − (Iyy − Ixx)(ψ̇c2θsφcφ)) c32 = Izzψ̇sθcθ − (Iyy − Ixx)(θ̇sθcφsφ + φ̇cθs 2 φ)− (Iyy + Ixx)(ψ̇sθcθc 2 φ − φ̇cθc2φ) c33 = Izz θ̇sθcθ − (Iyy + Ixx)(θ̇sθcθs 2 φ) + (Iyy − Ixx)φ̇c2θsφcφ. Substituting the inertia and the Coriolis matrices, the model for the UAV attitude dynamic is given by Mηη̈ +Cηη̇ = uη + dη. (2.28) with uη = {uφ uθ uψ} representing the moments with respect to the inertial frame of reference, and dη = {dφ dθ dψ}> the disturbance moments vector action on the UAV. 2.4.3 Dynamic Model Considering the equations from the translational dynamics ( Sec. 2.4.1), and rotational dynamics (Sec. 2.4.2), the dynamic model from the UAV is presented by the following system of ODEs Mq̈ +Cq̇ + fg = Bu+ d, (2.29) where the term d = { d>r d > η }> represents the disturbances forces acting on the model. Moreover, the model matrices for the full UAV dynamic are defined by the following 26 equations M = [ mI3×3 03×3 03×3 Mη ] C = [ 03×3 03×3 03×3 Cη ] B = [ Br 03×3 03×1 W> η ] (2.30) Note that the rotational dynamic is independent of the translational coordinates, however, the translational dynamic depends on the attitude, such as the control forces in these direction are projections from the thrust force. Equation (2.29) allows one to compute the UAV trajectory depending on the com- manded input history u(t) over time. It is also used to design the controller that stabilizes the system in a required position. 2.5 EULER-LAGRANGE EQUATIONS FOR CONSTRAINED SYSTEMS The typical problem of dynamics consists of finding the trajectory that a system takes from an instant t1 to an instant t2. For a system with independent coordinates, the solution is obtained by solving the Euler-Lagrange equations (Sec. 2.4), which minimizes the action computed through S = ∫ t2 t1 L(q, q̇)dt. This process is illustrated in Fig. 2.4, where the trajectory that satisfies the Euler-Lagrange equations is represented by a solid line and other trajectories are indicated by dashed lines. t î ĵ t1 t2 Figure 2.4: Generic trajectories connecting the configurations at the instants t1 and t2. The system takes the trajectory that satisfies the Euler-Lagrange equation indicated by the solid line. Mechanical systems can be subject to restrictions on their motion, which consists to kinematic restraints imposed on them, such as the rod connecting a pendulum to its pivot joint (LEMOS, 2007). Considering the illustrative pendulum from Fig. 2.5, the mass attached to the rod is constrained to the perimeter of a circle with radius equal to the length l, which is equivalent to state that the mass position (xp yp) is constrained by the equation Θ = x2p + y2p − l2 = 0. A convenient approach used to model the constrained system consists of reparame- terizing the coordinates to obtain a new set of independent coordinates, such that the 27 O î ĵ l (xp yp) constraint perimeter Figure 2.5: Representation of the pendulum. constraint equations are identically satisfied in terms of the new coordinates. In the case of the illustrated pendulum, this strategy is alternatively achieved by parameterizing the pendulum mass using polar coordinates. Table 2.3 shows the constraint equations in terms of the Cartesian and polar coordinates. The first row indicates the constraint equation in the position level, whereas second and third rows represent the constraint equations in the velocity and acceleration levels, respectively. By using Cartesian coordinates, the solution consists to find the values of xp and yp such that the left-hand side of those equations are satisfied. On the other hand, using polar coordinates by defining xp = lcθ and yp = lsθ, the constraint equations are satisfied for any value of θ. Table 2.3: Comparative of the kinematic constraints from the pendulum geometry using Cartesian and polar coordinates. Cartesian coordinates polar coordinates Θ x2p + y2p − l2 = 0 (lsθ) 2 + (lcθ) 2 − l2 = 0 l2 (s2θ + c2θ)− l2 = 0 0 = 0 Θ̇ 2 (xpẋp + ypẏp) = 0 2 [ (lsθ) lθ̇cθ − (lcθ) lθ̇sθ ] = 0 2 [ l2θ̇sθcθ − l2θ̇sθcθ ] = 0 0 = 0 Θ̈ 2 ( ẋ2p + ẏ2p + xpẍp + ypÿp ) = 0 2 [ l2θ̇2c2θ + l2θ̇2s2θ + l2θ̈cθsθ − l2θ̇2s2θ ] = 0 2 ( l2θ̇2 − l2θ̇2 ) = 0 0 = 0 Although this illustrative system is described above, the change of coordinates to obtain an independent set of equation usually is not a trivial task. For these cases, the Lagrangian is modified to include the holonomic constraints, such that the associated Euler-Lagrange equations for such system are given by (SALETAN; CROMER, 1970; NETO, 2013) d dt ( ∂L ∂q̇i ) − ∂L ∂qi + ∂Θ ∂qi > Λ = fi (2.31) with Λ = {λ1 λ2 . . . λm} denotes the m Lagrange multipliers associated to each constraint equations Θ = {Θ1 Θ2 . . . Θm}>. The result is a system of both differential and algebraic 28 equations (DAE) given by Mq̈ + ∂Θ ∂q > Λ =f Θ =0 (2.32) The first row from Eq (2.31) corresponds to the n second order differential equations to solve for the n + m variables q and Λ. The second row consists of the remaining m algebraic equations to be solved simultaneously, which introduces additional challenges to find an accurate solution (BRAUN; GOLDFARB, 2009). 2.6 SINGLE-LIFT DYNAMIC The EOM from the single-lift system using the reparameterized coordinates is obtained in Silva, Bueno, and Santos (2021) and reproduced herein. Figure 2.6 illustrates the single-lift model, where the payload position is parameterized with spherical coordinates using the rod length l, and angles α and β. O î ĵ k̂ rp îb ĵb k̂b O′ î1, î2 k̂1 ĵ2 k̂2r Figure 2.6: Single-lift configuration with an UAV transporting a suspended payload. The position in terms of the reparameterized coordinates q̃ = {x y z φ θ ψ α β}> is computed through Eq. (2.33). rp (q̃) =  x+ lcαsβ y − lsα z + lcαcβ  (2.33) The equation of motion for from the single-lift system with reparameterized coordinates is given by the following equation M̃ ¨̃q + C̃ ˙̃q + f̃g = B̃uu (2.34) 29 such that M̃ = (m+mp) I3×3 03×3 M̃r 03×3 Mη 03×2 M̃> r 02×3 M̃p  C̃ = 03×3 03×3 C̃r 03×3 Cη 03×3 02×3 02×3 C̃p  (2.35) M̃r = −mplsαsβ mplcαcβ −mplcα 0 −mplsαcβ −mplcαsβ  (2.36) C̃r =  mpl(−cαsβα̇− sαcββ̇) mpl(−sαcβα̇− cαsββ̇) mplsαα̇ 0 mpl ( −cαcβα̇ + sαsββ̇ ) mpl ( sαsβα̇− cαcββ̇ )  (2.37) M̃p = [ mpl 2 0 0 mpl 2c2α ] C̃p = [ 0 mpl 2sαcαβ̇ −mpl 2sαcαβ̇ −mpl 2sαcαα̇ ] (2.38) f̃g =  02×1 − (m+mp) g 03×1 mpglsαcβ mpglcαsβ  (2.39) B̃u = [ B 02×4 ] (2.40) Note that in the single-lift with inelastic suspension model, the suspended payload exert no influence on the UAV rotational dynamic. On the other hand, rigidly coupling the payload to the UAV body increases the rotational inertia. A disadvantage is the introduction of a pendulum mode in the system dynamic which is discussed in Section 2.8. 2.7 MULTI-LIFT FORMATION DYNAMIC Figure 2.7 shows the multi-lift formation in which the ith UAV is denoted by Vi, with i = 1, . . . , n. The system consists of an arbitrary number n ≥ 1 of UAVs transporting a single suspended payload of mass mp. Attached to each UAV center of mass there is an inelastic suspension cable (i.e., a no mass rigid rod) of length li connecting the payload to each ith center of mass. 30 î ĵ k̂ O rp dp,x dp,y dp,z li îb,i ĵb,i k̂b,i Vi l1 îb,1 ĵb,1 k̂b,1 V1 l2 îb,2 ĵb,2 k̂b,2 V2 Figure 2.7: Multi-lift UAV formation with three illustrative UAVs transporting the suspended payload. Note the first, second and ith UAVs illustrated. The inertial reference frame I is located at the origin O and it is defined by the orthonormal basis {î ĵ k̂}. A body reference frame Bi with orthonormal basis {îb,i ĵb,i îk,i} attached to each ith UAV center of mass Oi and rotating solidarity to the vehicle is defined. The coordinates that unique define the multi-lift configuration is denoted by, qf (t) = { q1(t) > q2(t) > . . . qn(t)> qp(t) >}> (2.41) where the vector qi(t) = { ri(t) > ηi(t) >}> corresponds to the position and attitude from the ith UAV. Moreover, qp(t) = rp(t) = {xp yp zp}> denotes the position of the payload with respect to the inertial frame. Due to the rigid rods assumption, the multi-lift formation coordinates are constrained to satisfy the constant length equation, which is defined for each ith rod by the equation gi(qf ) = (xp − xi)2 + (yp − yi)2 + (zp − zi)2 − l2i = 0. The constraint equation from the n rods are combined such that Eq. (2.42) represents the set of holonomic constraints from the multi-lift formation. gf = {g1 g2 . . . gn}> = 0 (2.42) The dynamics of the multi-lift formation is obtained by using the Euler-Lagrange equations and the holonomic constraints handled using Lagrange multipliers (GOLDSTEIN, 2002). This approach allows one to generalize the equation of motion for an arbitrary number of UAVs. The Lagrangian from the multi-lift formation system is composed of the 31 kinetic and potential energies from the UAVs and payload. Tf = 1 2 mpq̇ > p q̇p + 1 2 n∑ i=1 mṙ>i ṙi + ω>i Jiωi, (2.43) Uf = −mpgzp − n∑ i=1 mgzi. (2.44) Additionally, there are n Lagrange multipliers denoted by Λf = {λ1 λ2 . . . λn} from the holonomic constraints that are added such that the Euler-Lagrange equations yield the correct equation of motion for the multi-lift system (SALETAN; CROMER, 1970). Then, the Lagrangian from the multi-lift system is given by Lf = Tf − Uf + g>f Λf , (2.45) Applying the Euler-Lagrange equation for each coordinate, and considering Eq. (2.42), the equation of motion from the multi-lift formation is given by the following system of second order differential and algebraic equations Mf q̈f +Cf q̇f + fg +G>f Λf =Bfuf + df gf =0 (2.46) where Mf is the formation inertia, Cf is the Coriolis matrix, fg is the gravitational force vector. The matrix Gf = ∂gi ∂qf is the Jacobian from the constraint equations and the term G>f Λf corresponds to the constraint forces. Bf is the formation input matrix and uf is the input vector for all UAVs. In addition, df represents the disturbance force vector acting on the system. The matrices from Eq. (2.46) are computed as Mf = diag (M1 M2 . . . Mn Mp) (2.47) Cf = diag (C1 C2 . . . Cn 03×3) (2.48) Bf = diag (B1 B2 . . . Bn) (2.49) fg = { f>g,1 f > g,2 . . . f>g,n f > g,p }> , (2.50) uf = { u>1 u > 2 . . . u>n }> (2.51) df = { d>1 d > 2 . . . d>n d > p }> (2.52) The formation gravitational forces from the ith UAV and the payload are computed by fg,i = −migk̂ and fg,p = −mpgk̂, respectively. Table 2.4 summarizes the dimensions of these involved matrices. 32 Table 2.4: Dimensions of the matrices and vectors for the multi-lift equation of motion. matrix dimension Mf 6n+ 3× 6n+ 3 Cf 6n+ 3× 6n+ 3 Gf n× 6n+ 3 Bf 6n+ 3× 4n fg 6n+ 3 uf 4n df 6n+ 3 2.7.1 Aerodynamic Drag Model The disturbance force vector acting on the payload is included in the model from Eq. (2.46) through dp = {dp,x dp,y dp,z}>. The drag force due to for the wind resistance on the payload is given by (ANDERSON, 1991) dp = −cpdṙp|ṙp|, (2.53) where cpd is the coefficient of drag, whose value is given in terms of the Reynolds number and the shape from the payload. The Reynolds number is an adimensional parameters which expresses the ratio between inertial and viscous forces, and it is given by (ANDERSON, 1991, p.32) Re = ρ∞v∞lref µ∞ , (2.54) where ρ∞ is the specific mass from the fluid, v∞ is the freestream velocity, µ the dynamic viscosity from the fluid, and lref a reference dimension, typically the chord for wings or the diameter in the case of a sphere. 2.7.2 Gust Model The 1-cosine gust model presented in Eq. 2.55 is a discrete gust for evaluating the impact of atmospheric turbulence on the system dynamics. It is parameterized by the gust gradient distance H, a reference gust velocity Uref and the parameter r representing the position with respect to the gust. The gust acts either vertically or laterally and it is consist of a pulse gradually increasing from 0 to Uref at the center of the gust. Figure 2.8 illustrates the normalized gust profile acting vertically in a position (xgo ygo) on the grid. The model is incorporated to the equation of motion by defining 1-cosine gust in specific positions on the environment. Considering a gust centered at (xgo ygo) and a distance r = √ (x− xgo)2 + (y − ygo)2 from the gust origin, then the following equation model the gust velocity U Uref =  1 2 [ 1 + cos (πr H )] for r ≤ H 0 for r > H (2.55) 33 x0 y0 1 î ĵ Figure 2.8: Representation from the normalized 1-cosine gust equation at (x0 y0). 2.7.3 Index Reduction The DAE system requires an especial strategy to be solved (BRAUN; GOLDFARB, 2009). The algebraic constraint imposes a restriction on the configuration manifold and violate this constraint implies changing the length of the suspension cables. A typical strategy consists of differentiating the algebraic equations to change the problem from solving a system of DAE to solve a set of ODEs, which is known in the literature as an index reduction. However, the resulting system is mildly numerically unstable and an accurate solution requires the use of a constraint stabilization technique (BAUMGARTE, 1972). Differentiating the constraint equation twice in respect to the time t result in the constraint equations on the acceleration level (ASCHER et al., 1995), i.e., dgf (qf ) dt = dgf dqf dqf dt = Gf q̇f = 0, (2.56) d2gf (qf ) dt2 = Gf q̈f + Ġf q̇f = 0. (2.57) Then, considering ff = Bfuf + df − fg − Cf q̇f , the resulting model is the following system of differential equations[ Mf G>f Gf 0 ]{ q̈f Λf } = { ff −Ġf q̇f } , (2.58) 2.7.4 Constraint Stabilization Equation (2.58) usually is numerically unstable. This numerical issue occurs since the problem of solving the algebraic equations is transformed in an equivalent equation in terms of acceleration g̈f = 0. However, numerical approaches commonly produce residuals that are integrated and accumulated on the velocity and position levels, resulting in the problem know as constraint drift (MASARATI, 2011), which is violation from the constraints equation on position and velocity levels, i.e., Eq. (2.42) and Eq. (2.56) respectively. 34 Baumgarte (1972) suggests the addition of the numerically stabilizing terms −2αf ġf and −β2 fgf to Eq. (2.57) to modify the constraint hypersurface to an attractor. Bisgaard, Bendtsen, and Cour-Harbo (2009) interpret this modification as an artificial introduction of spring and dampers in parallel with the rigid rods. This method is illustrated by Figure 2.9, where the non-stabilized solution diverges from the exact solution over time, whereas the stabilized case, the solution remains close to accurate value. However, note qi time exact α = β = 0 α, β 6= 0 Figure 2.9: Illustration from the constraint drift and the stabilization method. that both new terms are equal to zero. Then, the model from the multi-lift including the stabilizing terms assumes the following form[ Mf G>f Gf 0 ]{ q̈f Λf } = { ff −Ġf q̇f − 2αf ġf − β2 fgf } . (2.59) Then, the Lagrange multiplier with the numerically stabilizing terms can be written by Λf = − ( GfM −1 f G > f )−1 (−Ġf q̇f − 2αf ġf − β2 fgf −GfM −1 f ff ) , (2.60) The matrix GfM −1 f G > f is invertible since rank(Gf) = n, ∀ qf |gf = 0, and this property is demonstrated by the following Theorem. Theorem 1. the constraint Jacobian Gf has rank n for all qf such that gf (qf ) = 0. Proof. The rows from Gf are linearly indepedent if, and only if, the trivial solution α1 = α2 = · · · = αn = 0 is the unique solution to this following equation α1 dg1 dqf + α2 dg2 dqf + · · ·+ αn dgn dqf = 0. (2.61) The equation from the ith rod depends only on the coordinates from the ith UAV and payload, i.e., gi ≡ gi(qi qp), such that the partial differentiation with respect to the coordinates from the other UAVs is zero, i.e., ∂gi ∂qj = 0, ∀ i 6= j. (2.62) This implies that the ith row ∂gi ∂qi from Gf can not be expressed as a linear combination of gj ∂qi , i 6= j. Then, the trivial solution is the only possible solution for Eq. (2.61). An 35 additional condition is that the rows from Gf do not correspond to be the null-vector, i.e., dgi dqf 6= 0, i = 1 2 . . . n. Consider, by contradiction that dgi(qf ) dqf = 0, (2.63) then, if Eq. (2.63) is satisfied, the dot product is such that dgi(qf ) dqf dgi(qf ) dqf > = 0. (2.64) Writing dgi(qf ) dqf = [ ∂gi ∂q1 ∂gi ∂q2 . . . ∂gi ∂qn ∂gi ∂qp ] , (2.65) and from the holonomic constraint written in terms of the corresponding UAV and payload coordinates, i.e., gi(qf ) ≡ gi(qi, qp), the non-null terms are: ∂gi(qi, qp) ∂qi = [ −2(xp − xi) −2(yp − yi) −2(zp − zi) 0 0 0 ] (2.66) and ∂gi(qi, qp) ∂qp = [ 2(xp − xi) 2(yp − yi) 2(zp − zi) ] . (2.67) which allows one to rewrite Eq. (2.64) by dgi(qf ) dqf dgi(qf ) dqf > = 4 (xp − xi)2 + 4 (yp − yi)2 + 4 (zp − zi)2 . (2.68) Substituting the holonomic constraint from the ith rod result in: dgi(qf ) dqf dgi(qf ) dqf > = 4l2i , (2.69) which is zero only if li = 0. However, li 6= 0, and then, the rows Gf(qf) are non null for all qf that satisfies gf (qf ) = 0, which demonstrates that the matrix Gf is invertible quod erat demonstrandum. Substituting the stabilized Lagrange multiplier from Eq.(2.60) in the first row of Eq. (2.58), and defining the variable vf = q̇f leads to a system of first order differential equations: q̇f = vf v̇f = M−1 f ( ff +G>f ( GfM −1 f G > f )−1(− Ġf q̇f −GfM −1 f ff − 2αf ġf − β2 fgf )) (2.70) which represents the EOM for the multi-lift formation. Integrating Eq. (2.70) using a numerical method, such as the Runge-Kutta method, the trajectories of the multi-lift formation system can be obtained over time. 36 2.7.5 Algorithm Algorithm 1 summarizes the procedure for obtaining the trajectories from the multi-lift system. The algorithm consists of defining the number n of UAVs in the formation and their initial configuration, such that: the UAVs are on the constraint manifold. Then, the state from the formation at instant (t+ ∆t) is obtained by integrating Eq. (2.70) using the Runge-Kutta algorithm. The process is repeated until t achieves tmax. Algorithm 1 Compute the trajectories from the multi-lift system 1: n← Z+. 2: define the acceptable error in position εgf and velocity εġf levels. 3: define the simulation and stabilization parameters tmax, ∆t, αf , and βf . 4: define the initial system configuration qf (0) satisfying gf (qf ) = 0 and their derivatives according to Eq. (2.42), Eq. (2.56), Eq. (2.57). 5: main loop: 6: for t ≤ tmax do 7: qf (t+ ∆t)← integrate Eq. (2.70) (qf (t),uf (t), n) 8: uf (t+ ∆t)← update the control signal. 9: t← t+ ∆t. 10: check if the constraints are satisfied for all t. 11: if |gf (t)| < εgf and |ġf (t)| < εġf then 12: accept result and plot trajectories. 13: else 14: redefine the simulation and stabilization parameters. 15: goto main loop. 2.7.6 Null Space Formulation The null-space from Gf , i.e., ker (Gf ) is one subspace that maps any vector from it to the null vector. If one select a matrix Nf whose columns forms the basis from the null-space of Gf , then the following equation holds GfNf = 0 (2.71) where Nf is a 6n+ 3× 5n+ 3 matrix that spans the null-space from Gf . Applying the null-space in the EOM from the multi-lift dynamics (Eq.(2.46)), the result is the projection such that the constraint forces are zero, as represented by the following system of equations N>f Mf q̈f = N>f ff (2.72) 2.8 FORMATION MODAL ANALYSIS The geometric coupling of the multi-lift formation implies to a dynamics for which an equivalent modal analysis can be performed. This analysis can be carried out for the 37 single-lift configuration, since the coordinates are reparameterized analytically and the EOM linearized. However, considering the multi-lift system from Eq. (2.46), there are n constraints, such that only (5n+ 3) from those (6n+ 3) coordinates are independent from each other. This increases the complexity since a choice of reparameterization is not trivial. Yang (1992) noted that a linear system under linear and time-invariant constraints is stiffer than the unconstrained system. For the nonlinear constraints, such as Eq. (2.57), a stable and linear system remains stable for small vibrations around a equilibrium condition. In addition, the constraints introduce additional excitation sources and modify the disturbances forces, such that their effects on the system dynamics must be evaluated. Berrios et al. (2014) present such analysis for the twin-lift system and find unstable formation modes. Sohoni and Whitesell (1986) linearizes the EOM around an operating point and assumes a linear and time-invariant map from the independent and dependent coordinates. Then, they eliminate both the dependent coordinates and Lagrange multipliers to analyze the modal properties. Liang and Lance (1987) use the Gram-Schmidt algorithm to generate an orthonormal basis using the constraint Jacobian such that both the null-space and its time-derivative are efficiently obtained. They obtain the EOM using the velocity transformation in terms of the independent coordinates which allows one to compute the trajectories without using a constraint stabilization method González et al. (2017). The main disadvantage from using independent coordinates is the extra-step required to solve the non-linear algebraic equations to find the dependent coordinates at each time-step, which is prone to numerical convergence issues. 2.8.1 Single-Lift Modal The linearized EOM for the single-lift system with independent coordinates is obtained in Silva, Bueno, and Santos (2021). The modal properties from the single-lift system are obtained from the eigenvalues of matrix A defined by A = [ 08×8 I8×8 Ap 08×8 ] where Ap =  0 0 0 0 −g(m+mp) m 0 0 gmp m 0 0 0 g(m+mp) m 0 0 −gmp m 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 g(m+mp) ml 0 0 −g(m+mp) ml 0 0 0 0 0 g(m+mp) ml 0 0 −g(m+mp) ml  38 Two of the eigenvalues from A are identical purely imaginary pairs λ1 = λ2 = ±2.13j and they correspond to the pendulum modes from the single-lift system with a period of 2.94 s. The coordinate vector qe = {x y z 0 0 ψ x y z + l}> is a configuration in equilibrium as illustrated in Fig. 2.10. O î ĵ k̂ îb ĵb k̂b Figure 2.10: Schematic representation of the equilibrium configuration for the single-lift. 2.8.2 Twin-Lift Modal The formation modal analysis for the twin-lift configuration is realized by identifying a subset of dependent coordinates from the formation coordinates, i.e., qf , such that the constraint Jacobian with respect to these dependent coordinates has full rank. For the twin-lift system, the altitude from the leader and follower vehicles are choosen as dependent coordinates. The multi-lift equilibrium configuration is the solution from Eq. 2.46 with q̇f = 0 fg +G>f Λf =Bfuf gf =0 (2.73) For the twin-lift system, the equilibrium presented in Figure 2.11 is a solution. O î ĵ k̂ Figure 2.11: Schematic representation of the equilibrium configuration for the twin-lift. The dynamic equations from the twin-lift system are linearized with respect to the equilibrium and the modal properties obtained. Figure 2.12 shows the eigenvalues from 39 the linearized twin-lift configuration. −8 −6 −4 −2 0 2 4 6 8 ·10−2 −2 −1 0 1 2 < = Figure 2.12: Eigenvalues from the twin-lift configuration. Chapter 3 CONTROL OF UAV FORMATIONS The use of electronic controllers is fundamental for multirotors UAVs due to their unstable dynamics, nonlinearities, and underactuated system configuration, which render this type of vehicle challenging for piloting without computer assistance (EMRAN; NAJJARAN, 2018). A common requirement from electronic controllers is a reference model from the plant to be controlled, which contains the information in the form of mathematical equations or nonparametric functions to predict the future states depending on the applied control. An assumption usually employed in literature is that the model is time-invariant and the differences in relation to its real dynamic are regarded as disturbances, which are classified into two categories: a) internal disturbances from unmodeled dynamics (e.g., simplifications of aerodynamics (SANZ et al., 2016), time delays in control (FRIDMAN; SEURET; RICHARD, 2004), and model uncertainties (LU; REN; PARAMESWARAN, 2020)); b) external disturbances due to different other sources, such as due to wind gusts (GUO et al., 2020). In addition to the disturbances that a single UAV is subjected to, in aerial transportation mission using a multi-lift UAV formation, the coupled payload introduces additional forces on the system, since the UAV needs to provide extra thrust to balance the constraint force, which may saturate the UAV actuators. This chapter presents a discussion from control of UAVs, with reference to the control literature in Section 3.1. Section 3.2 introduces the developed strategy to control a multi-lift UAV formation. 3.1 RELATED WORK The PID controller is maybe the first controller that one encounters while researching for control of UAVs. This occurs due to its extensive use both in the hobbyist and scientific communities due to its ease of implementation, low computational requirements, and simple understanding from the controller parameters. To design a PID controller, the multirotor UAV dynamic is linearized in a hover condition, then one typically choose the PID gains 41 such that the closed-loop error dynamic from the linearized model is asymptotically stable. Sliding Mode Control is also a strategy present in diverse literature to control UAVs. One of its features is their excellent capability of rejecting disturbances by rapidly switching the actuator to stabilize the system. However, the switching also introduces high frequency vibrations on the UAV which is denominated as “chattering”. Artificial Intelligence (AI) is also in research to overcome the limitation from classical controllers, that include the dependence of models and online adapting to a changing environment. Waslander et al. (2005) model the UAV dynamic using input and output data with a non-parametric model. Then, they optimize the controller gains by adding a Gaussian noise to it in each iteration. Hwangbo et al. (2017) uses a neural network as the controller whose gains are updated through reinforcement learning. Due to the underactuated characteristics of multirotor UAV, the control is designed in a cascade fashion, where an inner-loop is responsible of stabilizing the actuated degrees of freedom, and an outer-loop responsible of finding adequate set-points to allow the vehicle track the underactuated degrees of freedom. Raffo, Ortega, and Rubio (2011) and Brandão, Filho, and Carelli (2013) uses the partial-feedback linearization to design an inner-loop controller (SPONG, 1994). 3.1.1 Propulsion Saturation Multirrotor UAVs commonly uses multiple BLDC (brushless DC) eletric motors for both propulsion and control. The strategy is to keep the mechanical simplicity and reduce overall costs. The motor uses an electronic controller, such as an ESC (Electronic Speed Controller), to switch the DC voltage from the batteries, resulting in a changing magnetic field that the permanent magnets follows. Then, the motor speed can be controlled by conditioning the signal sent to the ESC. 3.2 FORMATION CONTROL Formation control consists of coordinating the motion from each UAV to achieve a common goal. In the multi-lift transportation mission, the proximity from each vehicle also introduces a collision hazard that requires implementation to detect and avoid capability. There are different ways to achieve this goal, which requires that an UAV acknowledge the position from its peers by: a) using the onboard cameras to estimate the obstacle path, then plan a route to avoid a collision (PEDRO et al., 2021); b) creating a communication network to share the information from the inertial navigation system; c) using a dedicated sensor such as ultrasonic or lidar (TAHIR et al., 2019). The individuals position are known over time, and they are used to develop a formation controller based on leader-follower strategy.The formation control consists of considering two types of UAVs: i) the leader, which is responsible of communicating with ground 42 control station and tracking the desired mission path, and ii) the followers, which are responsible to establish the formation geometry. The leader UAV has its control structure unaltered such that the stability and robustness properties are preserved. The followers UAVs have their outer-loop modified to use the leader position and determine their attitude set-points based on their current position within the formation. The desired formation geometry is defined using virtual constraints. If a distance d from the ith to the jth UAVs is desired, the corresponding constraint equation is simply the Euclidian distance given by (xi − xj)2 + (yi − yj)2 + (zi − zj)2 − d2 = 0. (3.1) If a linear separation in the î direction of dx units is desired, then the equation is (xi − xj)− dx = 0. (3.2) Similarly, for the ĵ and k̂ directions, the constraint equations are given by the following equations (yi − yj)− dy = 0, (3.3) (zi − zj)− dz = 0. (3.4) Table 3.1 presents different formation configurations and their corresponding virtual constraints. The first row presents the linear separation among the vehicles such that they keep the distance in both the î and ĵ axes. In the second row of the table, the UAVs are separated from each other by the Euclidian distance. The “V” and the circular formation geometries are obtained by combining the constraint equations such that every vehicle is constrained in the formation. The “V” formation is parameterized through the vertex angle 2γ and the distance d between UAVs. The formation control consists of keeping each UAV in its position, such that the virtual constraint equations are satisfied. If the virtual constraints exist, then a constraint force is considered to keep each UAV on its desired position according to the formation. Since there is no such force, it is artificially introduced by the followers, by projecting the thrust such that its components provide the computed virtual constraint forces. The virtual constraint equations are grouped in the following virtual constraints vector gvc = gvc(rf ) = 0, (3.5) which is written in terms of the formation configuration, described by the formation position vector rf = {r1 r2 . . . rn}>. This equation is linear in terms of the coordinates in “V” formations, and nonlinear for the circular formation. Figure 3.1 illustrates the formation controller. The formation dynamics is used in the integrator block, which computes the future configuration based on the previous control. This future state is validated to check if the constraints are approximately satisfied within 43 a tolerance, if not, the solution is interrupted for adjusting the constraint stabilization parameters. Then, the information is processed and the state from the leader is updated. The leader is updated with the mission trajectory, then, it can compute its control for tracking the desired trajectory. The formation controller receives the leader state, then, the followers can compute their control through the virtual constraints. The control from the leader and the followers are validated through the saturation block. The saturation is applied if a control that result in angular velocity outside the operational range from the multirotors, . The saturated formation control is informed to the integration block and the process is repeated. Integration Formation Dynamic gf ≈ 0 Processing Invalid Simulation Leader Formation Controller Mission Follower 2 ... Follower n Saturation uf (t+ ∆t) q̇f = f(qf ) qf (t + ∆t) q1(t + ∆t) No q1(t + ∆t) u1(t + ∆t) u2(t + ∆t) un(t + ∆t) q1,d(t + ∆t) Figure 3.1: Diagram illustrating the process for the obtaining the formation trajectories. 44 Table 3.1: Different formation geometries and their corresponding virtual constraint equations. representation geometry virtual constraints dy dx V1 V2 linear separation x1 − x2 − dx y1 − y2 − dy z1 − z2 − dz V1 V2 d Euclidian distance (x1 − x2)2 + (y1 − y2)2 + (z1 − z2)2 − d V1 V2 V3 2γ d “V” x1 − x2 − dcγ y1 − y2 − dsγ z1 − z2 x1 − x3 − dcγ y1 − y3 + dsγ z1 − z3 + V2 V3 V1 d circular (x1 − x2)2 + (y1 − y2)2 + (z1 − z2)2 − d (x1 − x3)2 + (y1 − y3)2 + (z1 − z3)2 − d (x2 − x3)2 + (y2 − y3)2 + (z2 − z3)2 − d 3.2.1 “V” formation The geometry of the desired “V” formation is obtained by setting the coefficients v1i using the geometry parameters of the vertex angle 2γ and the desired distance d between each two UAVs on the “V”, such as illustrated in Fig. 3.2. The “V-formation” geometry is obtained by considering a linear separation of v1j units from Vi to V1, such that the following system of equations is obtained gvc = Jvcrf + vc (3.6) where gvc = {g1i . . . g1n}>, i = 2, 3, . . . n, denotes virtual constraint equations relating the desired position from each ith follower to the leader. Vector vc = { v>12 . . . v>1n }> , i = 2, 3, . . . n, introduces the desired distance from the follower to the leader, which for 45 the “V” formation it is parameterized by v1i = i 2 {−dcγ − dsγ 0}> for an even number of vehicles, and v1i = i−1 2 {−dcγ dsγ 0}> otherwise. î ĵ V1 V2 V3 V4 V5 2γ d Figure 3.2: Schematic representation of the “V” formation with a vertex angle 2γ and the relative distance d for every two consecutive UAVs represented by Vi, i = 1, . . . , 5. Jvc =  I3×3 −I3×3 03×3 . . . 03×3 I3×3 03×3 −I3×3 . . . 03×3 ... ... ... . . . ... I3×3 03×3 0 . . . −I3×3  . (3.7) Differentiating twice Eq. (3.6) with respect to the time, the equation obtained is Jvcr̈f = 0. The terms 2αc (Jvcṙf ) and β2 c (Jvcrf + vc), such that αc > 0 and βc > 0, are added to stabilize the controller (TARRAF; ASADA, 2002). Jvcr̈f + 2αc (Jvcṙf ) + β2 c (Jvcrf + vc) = 0 (3.8) The translational dynamics from the UAV formation is given by Mf,rr̈f + ff,r + J>vcΛc = uf,r, (3.9) where Mf,r = diag(m1I3×3 m2I3×3 . . .mnI3×3) and ff,r = { f>1,r f > 2,r . . . f > n,r }> with fi,r = −migk̂, and uf,r = {u1,r . . . un,r}> with ui,r = Br,iuT,i denoting the lift from the ith UAV projected in the x, y and z directions. The formation acceleration is substituted in Eq. (3.8) resulting the following equations JvcM −1 f,r ( uf,r − ff,r − J>c Λc ) + 2αc (Jvcṙf ) + β2 c (Jvcrf + vc) = 0 (3.10) − ( JvcM −1 f,rJ > vc ) Λc = −JvcM−1 c (uf,r − ff,r)− 2αcJvcṙf − β2 c (Jvcrf − cc) (3.11) Λc = ( JvcM −1 f,rJ > vc )−1 ( JvcM −1 c (uf,r − ff,r) + 2αcJvcṙf + β2 c (Jvcrf − cc) ) . (3.12) The term uvf = −J>vcΛc = { u>1,vf u > 2,vf . . . u>n,vf }> from translational dynamics acts as the virtual forces and they are artificially introduced by the followers to keep the UAVs 46 on the formation constraint hypersurface defined by Eq.(3.6). The virtual force developed by the ith UAV is written in terms of its components in the x, y and z directions, such that ui,vf = { ui,x,vf ui,y,vf ui,z,vf }> . Due to the underactuated configuration from multirotors UAVs, these computed virtual forces are obtained by considering the following attitude set-points: φd,i = sin−1 ( ui,x,vf sψ,i − ui,y,vf cψ,i uT,i ) , θd,i = sin−1 ( ui,x,vf cψ,i + ui,y,vf sψ,i uT,icφd,i ) . (3.13) These attitude set-points in addition to the leader altitude zd,i = z1 are then transmitted to the inner-loop controller. Chapter 4 RESULTS AND DISCUSSION This chapter presents the results to discuss the introduced approach. Section 4.1 presents a time-domain analysis from the trajectories of the single-lift system using the method of independent coordinates and the method of Lagrange multipliers for constrained systems. Section 4.2 presents the parameters employed to fully define the mission environment. Section 4.3 presents numerical data from the twin-lift configuration. Section 4.4 presents the trajectories from the multi-lift system, and Section 4.5 introduces a comparative study. 4.1 SINGLE-LIFT The proposed methodology employed for solving the system of constrained differential equations is evaluated in the time-domain by integrating the equation of motion, by comparing the trajectories obtained using the methods of independent coordinates and Lagrange multipliers. The simulations start with identical initial conditions and a constant thrust is imposed to statically balance the combined UAV and payload weight. This specific control input is motivated to avoid compromising the trajectories due to the closed-loop control. Then, any difference in the results is attributed to the solving methodology. Table 4.1 presents the physical parameters which define the single-lift system. The trajectories from the single-lift with reparameterized coordinates (see Eq. 2.34) is obtained by integrating the EOM using the 5th order Runge-Kutta algorithm and it is a reference solution used for comparison. The trajectory from the single-lift system using the Lagrange multiplier method is obtained by setting n = 1 and the resulting EOM is solved through the Algorithm 1. Table 4.1: Single-lift system parameters. parameter description value unit m UAV mass 2.5 kg mp payload mass 0.4 kg J inertia tensor diag (0.957 0.486 0.355) kg.m2 l suspension length 2.5 m 48 Figure 4.1 shows the UAV trajectories and a snapshot from the configuration at t = 0 s. The payload is initially placed at α|t=0 = 15◦ and β|t=0 = 25◦ and these initial conditions correspond to rp|t=0 = {1.02 − 0.65 2.19}> m. The oscillation exhibits a period equal to 2.94 s, corresponding to the eigenvalue obtained from the formation modal analysis. −1.5 −1 −0.5 0 0.5 1 1.5 −1 0 1 0 2 x - [m] y - [m] z - [m ] r rp t = 0s Figure 4.1: Trajectories from the single-lift system. Figure 4.2 presents the constraint equation from the single-lift system in position level, i.e., g = (x− xp)2+(y − yp)2+(z − zp)2− l2 for different values of stabilization parameters. The solution starts from a configuration consistent with the imposed constraints (BRAUN; GOLDFARB, 2009). The constraint rapidly increases when computing the non-stabilized case, which is represented by the curve g(q)|αf=0, and then, the resulting trajectory becomes inaccurate. Including the constraint stabilization, represented by the curves g(q)|αf=20 and g(q)|αf=100, the constraint equation is satisfied within εgf < 10−1. The UAV and payload positions obtained by integrating the dynamic model using the Euler-Lagrange formulation for constrained system are approximately equal to the positions obtained by using reparameterized coordinates. As noted by Neto and Ambrósio (2003), the constraint stabilization method simply keeps the constraint violations under control, such that the error shown in Fig. 4.2 is integrated resulting in the increasing differences over time from the trajectories shown in Fig. 4.3. However, the UAV and payload system is designed to operate under a closed-loop controller, and then, the observed differences related to these different approaches to model the dynamics are insignificant. 49 0 5 10 15 20 25 30 35 40 45 50 0 1 2 3 4 5 6 ·10−2 time - [s] co n st ra in t eq u at io n - [m 2 ] g (q) |αf=0 g (q) |αf=20 g (q) |αf=100 Figure 4.2: Constraint equation from the single-lift system for different stabilization parameters. 50 0 10 20 30 40 50 −1 0 1 ·10−3 time - [s] x (q̃ ) − x (q ) - [m ] 0 10 20 30 40 50 −1 −0.5 0 0.5 1 ·10−2 time - [s] x p (q̃ ) − x p (q ) - [m ] 0 10 20 30 40 50 −1 −0.5 0 0.5 1 ·10−3 time - [s] y (q̃ ) − y (q ) - [m ] 0 10 20 30 40 50 −5 0 5 ·10−3 time - [s] y p (q̃ ) − y p (q ) - [m ] 0 10 20 30 40 50 −4 −2 0 2 4 ·10−4 time - [s] z( q̃ ) − z( q ) - [m ] 0 10 20 30 40 50 −2 0 2 ·10−3 time - [s] z p (q̃ ) − z p (q ) - [m ] Figure 4.3: Position error from the UAV and payload for αf = 100. 51 4.2 MISSION DEFINITION The transportation mission considered herein consists of assigning a desired path for the leader UAV, whereas the followers collaboratively lift the suspended payload and maintain the relative distances through virtual constraints. The mission trajectory assigned to the leader is composed by four different parts: i) ascending with a climb rate of −0.3 m.s−1 for 10 seconds, ii) cruise in the î direction with velocity 0.5 m.s−1 for 15 seconds, iii) cruise in the ĵ direction with velocity 0.5 m.s−1 by 15 seconds, and iv) cruise in the −î direction with velocity 0.5 m.s−1 until t = 100 s, until the formation achieves the final flight time, i.e., tmax = 120 s. The atmospheric disturbance is a lateral gust acting on the payload in the +ĵ direction. This gust effect is positioned at xgo,1 = −4 m, ygo,1 = 6 m, with a gradient of 2 m and reference velocity equal to 3 m.s−1. A gust load with similar characteristics action in the +î direction is defined at xgo,2 = −12 m and ygo,2 = 6 m. The system parameters are further presented in Tab. 4.2, and the UAVs for both the following investigations are assumed to be identical. The maximum thrust due to the rotor parameter is 32.28 N, which limit the maximum payload to mp = 0.79 kg. However, this upper limit is not practical because the attitude control needs to balance the rotors speed to rotate the UAV by reducing the lift and the UAV altitude control can fail. Then, a payload with weight exceeding such limit is defined, and the formation controller parameters αvc and βvc are arbitrary chosen. Table 4.2: Aerial transportation mission parameters. parameter description value unit m UAV mass 2.5 kg J inertia tensor diag (0.957 0.486 0.355) kg.m2 lr arm length 0.2 m kT thrust coefficient 1.83 · 10−6 N.s2/rad2 kD moment coefficient 1.85 · 10−7 N.m.s2/rad2 ωr,min minimum rotor speed 0 rad/s ωr,max maximum rotor speed 2100 rad/s mp payload mass 0.8 kg l suspension length 10.26 m cpd payload drag coefficient 0.12 - αvc formation derivative gain 2.1 - βvc formation proportional gain 1.3 - 4.3 TWIN-LIFT The dynamic model for describing the twin-lift system is obtained by setting n = 2 in Algorithm 1, and it consists of a formation with two identical UAVs coordinated by using 52 the leader-follower strategy. The formation is parameterized with a fixed separation of dx = 3.46 m, dy = −2.00 m and dz = 0.00 m between the leader and follower aerial vehicles (see Sec. 3.2). Figure 4.4 presents the resulting constraints in position level (see Eq. 2.42) from the leader and follower, denoted by the curves g1 and g2, respectively. The constraints are satisfied within an error of 2 · 10−2 m2, which is considered as an acceptable result for practical applications. During the hover phase, the constraints asymptotically converges to zero due to the stabilization method, such as previously noted by Baumgarte (1972). 0 20 40 60 80 100 120 −0.5 0 0.5 1 ·10−2 climb cruise+î cruise+ĵ cruise−î hover time - [s] co n st ra in ts - [m 2 ] g1 g2 Figure 4.4: Constraints over time from the twin-lift system. Figure 4.5 shows the top-view from the trajectories. The payload is set to an initial angle with respect to the oscillation plane formed by the UAVs and the vertical direction. This condition introduces an initial pendulum motion that is gradually suppressed due to the included drag force (see Eq. 2.53). According to these results, this oscillatory motion is further excited when the UAVs changes directions and when the payload moves through regions with gusts. 53 −25 −20 −15 −10 −5 0 5 10 −2 0 2 4 6 8 x - [m] y - [m ] r1 r2 rp gust gust Figure 4.5: Top-view from the twin-lift simulation. Figure 4.6 shows the roll and pitch angles from the UAVs. Note that the UAVs keeps its non-zero attitude. This behavior is required to balance the forces to constraint the payload, which are at an angle with respect to the vertical direction resulting in forces in the î and ĵ directions. These constraint forces are not included in the control design architecture, and they are considered as disturbances instead. They pull the UAVs towards each other, and this behavior increases the errors from the virtual constraints for the follower controller, and from the leader outer-loop. Then, a compensating control is considered (see Eq. 3.13), which results in that observed non-zero attitude. 0 50 100 −5 0 5 time - [s] ro ll an gl e φ - [◦ ] φ1 φ2 0 50 100 −5 0 5 time - [s] p it ch an gl e θ - [◦ ] θ1 θ2 Figure 4.6: Attitude from the UAVs considering the twin-lift configuration. Figure 4.7 shows the virtual constraint from the twin-lift formation control. The UAVs start at the desired formation, and during the climb phase the follower UAV exhibits a time delay to keep at the same altitude as the leader vehicle, such as indicated by the curve gvc|k̂. This behavior occurs since the reference for this parameter is updated in the outer-loop instead of the inner-loop, such as verified to the leader UAV. The follower 54 UAV accurately keeps the leader altitude for the remaining of the mission, even with the additional weight from the payload, and such performance is achieved due to the inner-loop controller robustness. In addition, regarding to the virtual constraints, i.e., gvc|̂i and gvc|ĵ , the payload pulls the UAVs close to each other, as shown by the difference from the follower UAV with respect to its position in the formation. 0 20 40 60 80 100 120 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 climb cruise+î cruise+ĵ cruise−î hover time - [s] fo rm at io n v ir tu al co n st ra in ts - [m ] gvc|̂i gvc|ĵ gvc|k̂ Figure 4.7: Virtual constraints from the twin-lift formation control. Figure 4.8 shows the thrust forces for both UAVs. During the climb phase, since the follower vehicle is flying at a lower altitude than the leader, the required thrust from the leader is 6.14% higher than that force required from the follower. The required thrust force is approximately distributed equally for the remaining of the mission. 0 20 40 60 80 100 120 −30 −25 −20 climb cruise+î cruise+ĵ cruise−î hover −29.33 N −27.63 N time - [s] th ru st fo rc e - [N ] uT,1 uT,2 Figure 4.8: Thrust force required from the UAVs in the twin-lift configuration. 55 −20 −15 −10 −5 0 5 −2 0 2 4 6 0 10 x - [m] y - [m] z - [m ] r1 r2 rp Figure 4.9: UAVs trajectories from the twin-lift configuration. 4.4 MULTI-LIFT The multi-lift system configuration is evaluated by considering n = 3, i.e., three UAVs. Two of them are defined as followers and their outer-loop controllers are replaced by a single controller designed according to Sec. 3.2. The desired formation geometry for this case is the “V” formation, and the geometric parameters are d = 4 m and γ = 30◦. Figure 4.10 shows the corresponding constraint equation. Similarly to the twin-lift simulation, the constraint stabilization parameters reduces the constraint error over time until they stabilize at zero. 0 20 40 60 80 100 120 −0.5 0 0.5 1 ·10−2 climb cruise+î cruise+ĵ cruise−î hover time - [s] co n st ra in ts - [m 2 ] g1 g2 g3 Figure 4.10: Constraints from the multi-lift system. Figure 4.11 presents the virtual constraints from the multi-lift system. There are 6 56 virtual constraints to define the “V” formation, such that 3 of them constrain the position of the first follower to the leader, and the other ones constrain the position of the second follower to the leader. Note that, there is a persistent error in this multi-lift configuration, but its value decreases in comparison with the twin-lift system. This behavior occurs because the payload weight is also distributed to the additional UAV. 0 20 40 60 80 100 120 −1 −0.5 0 0.5 climb cruise+î cruise+ĵ cruise−î hover time - [s] fo rm at io n v ir tu al co n st ra in ts - [m ] gvc|̂i, 1 gvc|ĵ, 1 gvc|k̂, 1 gvc|̂i, 2 gvc|ĵ, 2 gvc|k̂, 2 Figure 4.11: Virtual constraints from the multi-lift formation control. Figure 4.12 presents the thrust forces computed for this multi-lift configuration. Their magnitudes change more substantially when the formation changes the phase of flight. 0 20 40 60 80 100 120 −30 −25 −20 −15 climb cruise+î cruise+ĵ cruise−î hover time - [s] th ru st fo rc e - [N ] uT,1 uT,2 uT,3 Figure 4.12: Thrust forces from the UAVs in the multi-lift configuration. Figure 4.13 illustrates the top-view of the trajectories. The payload weight changes the formation geometry, resulting in the UAVs flying closer to each other. Figure 4.14 presents 57 the multi-lift system trajectories of the UAVs simultaneously lifting the heavy payload. −25 −20 −15 −10 −5 0 5 10 −2 0 2 4 6 8 10 x - [m] y - [m ] r1 r2 r3 rp gust gust Figure 4.13: Trajectories (Top-view) for the multi-lift configuration. −20 −15 −10 −5 0 5 0 5 0 10 position x - [m] position y - [m] p os it io n z - [m ] r1 r2 r3 rp Figure 4.14: Trajectories from the multi-lift system. 4.5 RESULTS The multi-lift system in general can reduce the payload oscillations. Figure 4.15 presents the payload position computed when considering the twin-lift configuration (xp|TL), and multi-lift system (xp|ML). Note that the twin-lift system exhibits a decaying oscillation over time, due to the effect of the drag force. On the oher hand, no significant oscillation is observed when considering multi-lift systems. In addition, note that the gust load 58 generates higher oscillations when considering the first gust, which shows that payload dynamic depends on formation geometry. 0 20 40 60 80 100 120 −20 −10 0 time - [s] p ay lo ad p os it io n x - [m ] xp|TL xp|ML 0 20 40 60 80 100 120 0 5 time - [s] p ay lo ad p os it io n y - [m ] yp|TL yp|ML Figure 4.15: Comparison between the payload position when considering the twin-lift and multi-lift configurations. Figure 4.16 shows a comparison of the virtual constraint error in the twin-lift and multi-lift configurations. Additional vehicle in the multi-lift configuration results in the decrease of the disturbance force exerted to each one from the vehicles. The effect is that the virtual constraint error is decreased in comparison with the twin-lift configuration. 0 20 40 60 80 100 120 −1 −0.5 0 0.5 1 time - [s] v ir tu al co n st ra in ts - [m ] gvc |̂i, 1 gvc|ĵ, 1 gvc|k̂, 1 0 20 40 60 80 100 120 −1 −0.5 0 0.5 1 time - [s] v ir tu al co n st ra in ts - [m ] gvc |̂i, 1 gvc|ĵ, 1 gvc|k̂, 1 gvc|̂i, 2 gvc|ĵ, 2 gvc|k̂, 2 Figure 4.16: Comparison between the virtual constraints when considering the twin-lift and multi-lift configurations. Figure 4.17 shows the comparison from the rotor speed from the leader vehicle for the twin-lift configuration (ωr,1|TL), and for the multi-lift configuration (ωr,1|ML). Note that the rotor speed for both configurations is below the saturation region, as indicated by the gray region. The multi-lift configuration allows one to reduce 2.37% from the rotor speed, for the illustrated time-instant. The immediate advantage is a reduction from the required power, as the rotor speed directly correlates with the energy consumption, thus increasing the multi-lift flight time. 59 0 20 40 60 80 100 120 1600 1800 2000 2200 rotor saturation climb cruise+î cruise+ĵ cruise−î hover 1975 rad 1928 rad time - [s] ro to r sp ee d - [r ad /s ] ωr,1|TL ωr,1|ML Figure 4.17: Comparison from the rotor speed from the leader in the twin-lift and multi-lift formations. Chapter 5 FINAL REMARKS The current technology readness level of UAVs allows one to use multirotor aerial vehicles for applications in which the cost and complexity of conventional aviation are prohibitive. However, due to the characteristics from the multirotor configuration, the applications are limited to those ones involving small weight and short-range flights. Formation flight can address the first limitation by collaboratively transporting heavy-payloads, i.e., distributing the additional weight on multiple vehicles, such that the force exerted by each one is within its operational range. Formation flight also provides means for including redundancy, since a failure of an individual UAV occurs, the failed vehicle can decouple from the formation, and the mission can be continued. However, the design of such features require a robust computational environment to evaluate the controller performance and stability before experimental tests with prototypes. This type of tool allows one to collect data and analyze the controller in different scenarios. This present work introduces a methodology to investigate the dynamics of a multi- lift UAV formation, by using the Euler-Lagrange formalism combined with Lagrange multipliers for dealing with the geometric constraints. The result contribute to future works involving designing controllers by providing a complete modeling for evaluating the system dynamics for a formation with an arbitrary number of UAVs. A leader-follower formation controller is investigated using virtual constraints to guide the UAVs in a formation through the desired trajectory. 5.1 CONCLUSIONS The strategy of using multi-lift systems offers advantages in comparison with the use of a single high-lift vehicle such as a lower acquisition and operational costs of a smaller UAV. The use of multiple vehicles is also more convenient to define a redundant configuration to increase the mission robustness, in relation to failure events for an instance, because the failed UAV can decouple from the formation and the mission can continue with the remaining UAVs. In addition, the use of three or more vehicles can reduce the pendulum 61 dynamics commonly verified in both single and twin-lift configurations. This reduction is benefit since the periodic disturbance introduced by the pendulum motion may adverse affect the system performance. The payload pendulum motion verified in both single-lift and twin-lift configurations is introduced when the UAVs change their directions or due to the wind gusts effects. The obtained dynamic model for the multi-lift configuration is written in terms of an arbitrary number of vehicles which allows one to investigate formations with any particular aerial transportation configuration, such as the single-lift, twin-lift, and multi-lift. The model is obtained by transforming a system of DAE into a system of ODE and involving an index-reduction with a constraint stabilization procedure. The approach introduces an extra-step of analyzing and tuning the stabilization parameters to obtain the constraints satisfied within a specified tolerance. The twin-lift and multi-lift configurations require a non-zero attitude from each UAV to stabilize the system, due to the additional weight from the suspended payload. This condition usually is very difficult to be achieved by a human remotely controlled formation. However, the results herein show that a multi-lift configuration using a formation controller using virtual constraints can be successfully employed to transport heavy-payloads. It is assumed that each follower vehicle knows its relative distance from the leader, and then, the followers can move to their designed position in the “V” formation geometry by artificially introducing virtual constraint forces. The formation controller reference model consists of disconnected UAVs moving without the suspended payload. This design consideration in a multi-lift mission considers the payload as a disturbance source. The result show that the payload pulls the UAVs close to each other. This approximation also increases the formation virtual constraints errors, and then a compensating control force is included to limit the pulling effect. 5.2 SUGGESTIONS FOR FUTURE WORK The following suggestions are proposed for future works in this field. � extend the modeling for rigid-body payloads and cable-suspension attachment off the UAVs center of mass; � design of a payload-based controller to reduce pendulum oscillations and keep the formation geometry; � design a controller to account for the propulsion-saturation and reorganize the formation to improve the payload-weight distribution; � implement formation fault-tolerant controllers to increase mission robustness. References ANDERSON, J. D. Fundamentals of aerodynamics. 2. ed. Maidenhead, England: McGraw Hill Higher Education, Mar. 1991. (McGraw-Hill series in aeronautical and aerospace engineering). ASCHER, U. M. et al. Stabilization of Constrained Mechanical Systems with DAEs and Invariant Manifolds. Mechanics of Structures and Machines, Taylor & Francis, New York, v. 23, n. 2, p. 135–157, 1995. DOI: 10.1080/08905459508905232. BAUMGARTE, J. Stabilization of constraints and integrals of motion in dynamical systems. Computer Methods in Applied Mechanics and Engineering, Elsevier BV, Amsterdam, v. 1, n. 1, p. 1–16, June 1972. DOI: 10.1016/0045-7825(72)90018-7. BERRIOS, M. et al. Stability, control, and simulation of a dual lift system using autonomous R-MAX helicopters. In: AHS 70th ANNUAL FORUM. Montréal: [s.n.], May 2014. v. 3, p. 2195–2213. BISGAARD, M.; BENDTSEN, J.; COUR-HARBO, A. la. Modeling of Generic Slung Load System. Journal of Guidance, Control, and Dynamics, Reston, v. 32, n. 2, p. 573– 585, Mar. 2009. DOI: 10.2514/1.36539. eprint: https://doi.org/10.2514/1.36539. BRANDÃO, A. S.; FILHO, M. S.; CARELLI, R. High-level underactuated nonlinear control for rotorcraft machines. In: 2013 IEEE INTERNATIONAL CONFERENCE ON MECHATRONICS (ICM). Vicenza: [s.n.], 2013. p. 279–285. DOI: 10.1109/ICMECH.2013. 6518549. BRAUN, D. J.; GOLDFARB, M. Eliminating constraint drift in the numerical simulation of constrained dynamical systems. Computer Methods in Applied Mechanics and Engineering, Amsterdam, v. 198, n. 37, p. 3151–3160, 2009. ISSN 0045-7825. DOI: https://doi.org/10.1016/j.cma.2009.05.013. DURANGO, G.; LAWSON, C.; SHAHNEH, A. Formation flight investigation for highly efficient future civil transport aircraft. The Aeronautical Journal, Cambridge University Press (CUP), Cambridge, v. 120, n. 1229, p. 1081–1100, June 2016. DOI: 10.1017/aer. 2016.52. EMRAN, B. J.; NAJJARAN, H. A review of quadrotor: An underactuated mechanical system. Annual Reviews in Control, Oxford, v. 46, p. 165–180, 2018. ISSN 1367-5788. DOI: https://doi.org/10.1016/j.arcontrol.2018.10.009. FOSSEN, T. Guidance and Control of Ocean Vehicles. [S.l.]: Wiley, 1994. ISBN 9780471941132. FRIDMAN, E.; SEURET, A.; RICHARD, J.-P. Robust sampled-data stabilization of linear systems: an input delay approach. Automatica, Elsevier BV, Oxford, v. 40, n. 8, p. 1441–1446, Aug. 2004. ISSN 0005-1098. DOI: 10.1016/j.automatica.2004.03.003. https://doi.org/10.1080/08905459508905232 https://doi.org/10.1016/0045-7825(72)90018-7 https://doi.org/10.2514/1.36539 https://doi.org/10.2514/1.36539 https://doi.org/10.1109/ICMECH.2013.6518549 https://doi.org/10.1109/ICMECH.2013.6518549 https://doi.org/https://doi.org/10.1016/j.cma.2009.05.013 https://doi.org/10.1017/aer.2016.52 https://doi.org/10.1017/aer.2016.52 https://doi.org/https://doi.org/10.1016/j.arcontrol.2018.10.009 https://doi.org/10.1016/j.automatica.2004.03.003 63 GOLDSTEIN, H. Classical mechanics. San Francisco: Addison Wesley, 2002. ISBN 978-0-201-65702-9. GONZÁLEZ, F. et al. Assessment of Linearization Approaches for Multibody Dynamics Formulations. Journal of Computational and Nonlinear Dynamics, New York, v. 12, n. 4, Jan. 2017. 041009. ISSN 1555-1415. DOI: 10.1115/1.4035410. GUERRERO-SÁNCHEZ, M. E. et al. Swing-attenuation for a quadrotor transporting a cable-suspended payload. ISA Transactions, Elsevier BV, New York, v. 68, p. 433–449, May 2017. DOI: 10.1016/j.isatra.2017.01.027. GUO, K. et al. Multiple observers based anti-disturbance control for a quadrotor UAV against payload and wind disturbances. Control Engineering Practice, Elsevier BV, Oxford, v. 102, p. 104560, Sept. 2020. ISSN 0967-0661. DOI: 10.1016/j.conengprac. 2020.104560. HAIN, J. et al. Feeding behavior of the Humpback whale, Megaptera novaeangliae, in the western North Atlantic. Fishery Bulletin, U.S. National Marine Fisheries Service * Scientific Publications Office, Seattle, v. 80, Jan. 1981. ISSN 0090-0656. HE, L. et al. Feedback formation control of UAV swarm with multiple implicit leaders. Aerospace Science and Technology, Elsevier BV, Issy-les-Moulineaux, v. 72, p. 327– 334, Jan. 2018. ISSN 1270-9638. DOI: 10.1016/j.ast.2017.11.020. HWANGBO, J. et al. Control of a Quadrotor With Reinforcement Learning. IEEE Robotics and Automation Letters, Institute of Electrical and Electronics Engineers (IEEE), Piscataway, v. 2, n. 4, p. 2096–2103, Oct. 2017. ISSN 2377-3766. DOI: 10.1109/ lra.2017.2720851. KOTARSKI, D. et al. Performance Analysis of Fully Actuated Multirotor Unmanned Aerial Vehicle Configurations with Passively Tilted Rotors. Applied Sciences, Basel, v. 11, n. 18, 2021. ISSN 2076-3417. DOI: 10.3390/app11188786. LEMOS, N. Mecânica anaĺıtica. São Paulo: Livraria da F́ısica, 2007. ISBN 9788588325241. LIANG, C. G.; LANCE, G. M. A Differentiable Null Space Method for Constrained Dynamic Analysis. Journal of Mechanisms, Transmissions, and Automation in Design, The American Society of Mechanical Engineers, New York, v. 109, n. 3, p. 405–411, Sept. 1987. ISSN 0738-0666. DOI: 10.1115/1.3258810. LIANG, Y.; DONG, Q.; ZHAO, Y. Adaptive leader–follower formation control for swarms of unmanned aerial vehicles with motion constraints and unknown disturbances. Chinese Journal of Aeronautics, Zhongguo Hangkong Xuehui, Beijing, v. 33, n. 11, p. 2972–2988, Nov. 2020. ISSN 1000-9361. DOI: 10.1016/j.cja.2020.03.020. LU, Q.; REN, B.; PARAMESWARAN, S. Uncertainty and Disturbance Estimator-Based Global Trajectory Tracking Control for a Quadrotor. IEEE/ASME Transactions on Mechatronics, Institute of Electrical and Electronics Engineers, Piscataway, v. 25, n. 3, p. 1519–1530, 2020. ISSN 1083-4435. DOI: 10.1109/TMECH.2020.2978529. LUO, C.; DU, Z.; YU, L. Neural Network Control Design for an Unmanned Aerial Vehicle with a Suspended Payload. Electronics, MDPI AG, Basel, v. 8, n. 9, p. 931, Aug. 2019. ISSN 2079-9292. DOI: 10.3390/electronics8090931. MASARATI, P. Constraint Stabilization of Mechanical Systems in Ordinary Differential Equations Form. Proceedings of the Institution of Mechanical Engineers, Part K: Journal of Multi-body Dynamics, SAGE Publications, London, v. 225, n. 1, p. 12–33, Feb. 2011. ISSN 1464-4193. DOI: 10.1177/2041306810392117. https://doi.org/10.1115/1.4035410 https://doi.org/10.1016/j.isatra.2017.01.027 https://doi.org/10.1016/j.conengprac.2020.104560 https://doi.org/10.1016/j.conengprac.2020.104560 https://doi.org/10.1016/j.ast.2017.11.020 https://doi.org/10.1109/lra.2017.2720851 https://doi.org/10.1109/lra.2017.2720851 https://doi.org/10.3390/app11188786 https://doi.org/10.1115/1.3258810 https://