Momentum Preserving Simulation of Cooperative Multirotors With Flexible-Cable Suspended Payload

Kaya, Uluhan Cem; Subbarao, Kamesh

doi:10.1115/1.4053343

Abstract

In this paper, a momentum-preserving integration scheme is implemented for the simulation of single and cooperative multirotors with a flexible-cable suspended payload by employing a Lie group–based variational integrator (VI), which provides the preservation of the configuration manifold and geometrical constraints. Due to the desired properties of the implemented VI method, e.g., symplecticity, momentum preservation, and the exact fulfillment of the constraints, exponentially long-term numerical stability, and good energy behavior are obtained for more accurate simulations of aforementioned systems. The effectiveness of Lie group VI method with the corresponding discrete systems are demonstrated by comparing the simulation results of two example scenarios for the single and cooperative systems in terms of the preserved quantities and constraints, where a conventional fixed-step Runge–Kutta 4 (RK4) and variable-step integrator are utilized for the simulation of continuous-time models. It is shown that the implemented VI method successfully performs the simulations with a long-time stable behavior by preserving invariants of the system and the geometrical constraints, whereas the simulation of continuous-time models by RK4 and variable step are incapable of satisfying these desired properties, which inherently results in divergent and unstable behavior in simulations.

1 Introduction

The aerial manipulation of payloads using suspended cables by multirotors has gained a considerable amount of attention for more than a decade among researchers. In addition to the simplistic, cost-efficient, and easily scalable design of multirotor systems, being an inexpensive and effective way of transporting a useful payload in various missions where no other means of transportation is either available or effective such as search and rescue missions or firefighting scenarios in remote/inaccessible areas makes these systems into one of the most tested platforms. Literature has numerous fruitful works done in the various aspects of this very topic. Sreenath et al. [1] proved that quadrotor with a slung load system is actually differentially flat, providing a simpler form of relations between quadrotor and slung load states. Following that, some hybrid approaches modeling the system according to the tautness of the suspension cable are proposed in Refs. 2 and 3 and are successfully used to generate feasible trajectories for tracking, which are experimentally validated as well. As opposed to the taut cable assumptions made for simplified modeling and control purposes by most studies, e.g., Refs. 1–5, Goodarzi et al. [6] adopt a flexible suspension cable consisting of multiple links such that deformations in cable shape can be addressed without an explicit analysis of tension force along the cable, which yields a more accurate approximation of actual continuous cable dynamics. Eventually, previously developed models and control algorithms are extended to cooperative multirotor systems where a rigid-body payload or a lumped-mass payload is suspended by multiple unmanned aerial vehicles (UAV) via either straight, as in Refs. 7–9, or flexible cables, as in Ref. 10.

All of these aforementioned works and design decisions are based on simulations or limited experimental verifications under the assumption that the actual behavior of the real system is captured reasonably well by the developed dynamical models. Leaving aside the accuracy of models capturing the underlying physical phenomena, it is known that the majority of numerical integration methods, including the Runge–Kutta based ODE solvers, do not preserve the symmetries (invariants) and the geometrical properties of the system, where even slight deviations might substantially affect the overall simulation accuracy, especially for the systems under consideration. Inherently, the conclusions drawn from inaccurate simulations of such complex physical systems may cause unexpected results in practice. These are particularly important in systems that involve flexible components, and a simple rigid body assumption does not capture the interactions between the components accurately.

To overcome described integration inaccuracies, a special class of geometric integration scheme, called variational integrators (VIs), is introduced in the literature [11]. VI approach, instead of discretizing continuous equations of motion as the most general-purpose integrators do, directly utilizes the discrete variational principle obtaining the discrete dynamics as a result of the discrete Hamilton's principle. Therefore, the derived discrete dynamics exactly preserves the momentum and the symplectic form of the system, as presented in great detail by Refs. 12 and 13. In addition, Lee and his team in Ref. 14 put forward an enhancement for variational integrators to preserve the geometric structure of the configuration manifold, which is represented as a Lie group, in rigid-body problems. Along with the symplecticity and momentum conservation, the exact preservation of the structure of the manifold ensures an exponentially long-term stability with a good energy behavior, which makes VI methods an ideal candidate for the simulation of complex or highly nonlinear systems.

As the precision and the long-term numerical stability are detrimental for the simulations of space missions, VI methods are one of the natural choices in these missions. The spacecraft dynamics in Ref. 15 is tackled by developing a variant of variational integrators involving quaternions to represent the spacecraft attitude and demonstrates the superiority of VI simulation accuracy among the other continuous solvers. Various other works utilizing VI methods and discrete control systems have been presented by Lee et al. [16,17], who have significant contributions in this area. Indeed, their representation of attitude dynamics via Lie group operations provides a basis for our N-Link cable model in this work as a series of spherical pendulums. Additionally, VI applications on a rigid-body dynamics of a quadrotor UAV in Ref. 18 and an underwater vehicle in Ref. 19 are among many other examples that led to our work. Apart from the aforementioned applications, Ref. 20 reveals another inspiring example of estimation and filtering for mechanical control systems, where multiple structure-preserving methods are incorporated with Kalman filters for higher accuracy in state estimations.

The main contribution of this work is the derivation and the implementation of Lie group variational integrator method to the transportation of a suspended payload via flexible cables by either single or cooperative multirotors. This work is dedicated to accomplish more accurate simulations of cooperative multirotors transporting a slung load via flexible cables by employing a variational integrator approach. The physical system consists of multiple multirotors, serially attached cable segments with concentrated masses and a point-mass payload that is suspended via these cables. Hence, the described system defines a multibody system with a complex geometry where the configuration manifold lies on a highly nonlinear space with geometrical constraints on each link, we benefited from Lie group VI approach in our work. We begin with the derivation of the complete equations of motion in continuous time for a single multirotor with a flexible cable-suspended payload, as in Ref. 6. Following that, the discrete dynamics of the same system is constructed with the Lie group variational integrator approach. Essentially, the steps for the construction of a single-vehicle case establish the foundation for the modeling of the general cooperative system whose continuous and discrete dynamics are both derived following the same steps. Acquiring the models of the single and cooperative systems in continuous and discrete time, the numerical simulations are performed to compare the simulation accuracies of a fixed-step fourth-order Runge–Kutta and ODE45 solvers for the continuous dynamics and the Lie group VI approach for the discrete dynamics in terms of their preservation of the system's total energy, momentum, and the constraints.

The rest of this paper is organized as follows. Section 2 outlines the steps involved in the derivation of the variational integrator procedure for the general Lagrangian systems. In Sec. 3, the derivation of both continuous and discrete time model of the single multirotor with a flexible cable-suspended payload system is presented, where the Lie group variational integration approach is employed to obtain the discrete model. Afterward, the procedure for a single-vehicle system is extended in Sec. 4 for the cooperative multirotors transporting a common payload suspended via flexible cables. Finally, Sec. 5 provides the simulation examples for both single and cooperative system demonstrating the effectiveness of the VI method. Section 6 concludes the paper.

2 Overview of the Variational Integrator Procedure

In this section, a brief introduction to the formulation of variational integrator method is presented for general Lagrangian systems. In the first step, the continuous equations of motion of the system are derived by utilizing the variational principle, and then, the Hamiltonian representation of the system is obtained by Legendre transformation. Afterward, the discrete equivalent of variational principle and Legendre transformation are employed similarly to provide the discrete Hamiltonian formulation, which constitutes the basis of variational integrator. As a result, symplectic and momentum-preserving properties of Hamiltonian systems are benefited from while formulating the variational integration scheme.

2.1 Continuous Dynamics.

General form of the Lagrangian of a dynamical system is expressed by

L (q, \dot{q}) = \frac{1}{2} {\dot{q}}^{T} M \dot{q} - U (q)

(1)

with the generalized coordinates, velocities, mass matrix, and the coordinate dependent potential function are denoted with q, $\dot{q}$ ⁠, M, and U, respectively.

Applying the variational principle on Lagrangian and integrating it along time while holding the endpoints

q (t), t = t_{0}, t_{f}

fixed, the least action principle yields the dynamics of the system, which is called Euler–Lagrange equations.

D_{q} L (q, \dot{q}) - \frac{d}{d t} D_{\dot{q}} L (q, \dot{q}) = 0

(2)

where the partial derivatives are denoted by $D_{q} L (q, \dot{q}) = \partial L (q, \dot{q}) / \partial q$ and $D_{\dot{q}} L (q, \dot{q}) = \partial L (q, \dot{q}) / \partial \dot{q}$ ⁠.

The motion evolves under a configuration-dependent potential in Eq. . However, the general forcing to the system can be included employing Lagrange d'Alembert principle as

D_{q} L (q, \dot{q}) - \frac{d}{d t} D_{\dot{q}} L (q, \dot{q}) + F_{q} = 0

(3)

Given the initial states,

(q, \dot{q})

⁠, the time histories of the states are provided by the integration of Eq. , a second-order ODE, along time. An alternate representation of Eq. can be achieved by Legendre Transformation, which yields Hamiltonian dynamics with the simpler first-order equations of motion.

\dot{p} - D_{q} L (q, p) = F_{q}

(4)

where $p = D_{\dot{q}} L (q, \dot{q})$ is called generalized momentum and it also establishes the symplectic form, (q, p).

2.2 Discrete Dynamics.

Instead of directly discretizing the continuous equations of motion in Eq. (3) or (4), the discrete equivalent of continuous dynamics is derived by applying the discrete version of the calculus of variation and minimum action principle.

The first step is to attain the discrete linear approximation of the Lagrangian as

L_{d} (q_{k}, q_{k + 1}) \approx L ((1 - α) q_{k} + α q_{k + 1}, \frac{q_{k + 1} - q_{k}}{h}) h

(5)

where the states are approximated by a linear interpolation as $q (t) \approx (1 - α) q_{k} + α q_{k + 1}$ and $\dot{q} (t) \approx \frac{q_{k + 1} - q_{k}}{h}$ for $α \in [0, 1]$ and t = kh with time-step of h. Many other linearization approaches are available in the literature with differing accuracy levels such as quadrature rules, $n -$ point linear interpolation and so on, as detailed in Refs. 21 and 22. One of the most common approaches is to use the trapezoidal approximation of Lagrangian, midpoint rule, with the selection of $α = 0.5$ ⁠, where we also adopted in the rest of this paper.

Then, applying the variational principle, discrete equivalent of the least action principle is formulated as a sum over time-steps, k.

δ B_{d} = \sum_{k = 0}^{N - 1} [D_{q_{k}} L_{d} (q_{k}, q_{k + 1}) \cdot δ q_{k} + D_{q_{k + 1}} L_{d} (q_{k}, q_{k + 1}) \cdot δ q_{k + 1}]

where $D_{q_{k}} L_{d}$ is the partial differential of the discrete Lagrangian at time-step k with respect to the generalized coordinates at kth time-step, q_k.

Knowing that the endpoints are held fixed,

δ q_{k} = 0

for

k = 0, N

⁠, discrete equations of motion are derived as

D_{q_{k}} L_{d} (q_{k}, q_{k + 1}) + D_{q_{k}} L_{d} (q_{k - 1}, q_{k}) = 0

Employing the discrete version of Lagrange d'Alembert Principle, forced discrete Euler–Lagrange equations are acquired as follows:

D_{q_{k}} L_{d} (q_{k}, q_{k + 1}) + D_{q_{k}} L_{d} (q_{k - 1}, q_{k}) + F_{k}^{-} + F_{k - 1}^{+} = 0

(6)

Discrete virtual work is again approximated by a linear combination of two consecutive states with the forcing as below

δ W = \int_{t_{k}}^{t_{k + 1}} F \cdot δ x d t \approx F_{k}^{-} \cdot δ x_{k} + F_{k}^{+} \cdot δ x_{k + 1} = δ W_{d}

where $F_{k}^{-} = (1 - c) h F_{k}$ and $F_{k}^{+} = c h F_{k + 1}$ with $c \in [0, 1]$ ⁠. For midpoint (trapezoidal) approximation, parameter c is chosen to be 0.5.

Given the initial states of (q₀, q₁), the next state, q₂, can be computed from Eq. and the integration progresses iteratively till the final state q_N, i.e.,

(q_{k - 1}, q_{k}) \to (q_{k}, q_{k + 1}), k = 1, . ., N - 1

⁠. However, instead of working with the consecutive initial states in each integration step, it is more convenient and natural to utilize the initial position and velocity (or momentum) states of the system to define its motion. The conversion between the pair of consecutive positions to position-momentum pair can be established by discrete Legendre Transformation, which also yields the discrete Hamiltonian representation of the system as follows:

p_{k} = - D_{q_{k}} L_{d} (q_{k}, q_{k + 1}) - F_{k}^{-}

(7)

p_{k + 1} = D_{q_{k + 1}} L_{d} (q_{k}, q_{k + 1}) + F_{k}^{+}

(8)

Given a pair of initial conditions, $(q_{0}, {\dot{q}}_{0})$ ⁠, the momentum state, p₀, is assessed from $p_{k} = D_{{\dot{q}}_{k}} L (q_{k}, {\dot{q}}_{k})$ relation. Then, substituting (q₀, p₀) pair into Eq. (7), q₁ can be solved implicitly. In the next step, Eq. (8) provides the next momentum state, $p_{k + 1}$ ⁠, explicitly. Hence, the flow of the integration follows $(q_{k}, p_{k}) \to (q_{k + 1}, p_{k + 1}), k = 0, . ., N - 1$ ⁠, as summarized in Algorithm 1.

Algorithm 1:

Variational integrator algorithm

IC:

(q_{0}, {\dot{q}}_{0})

//Initial Pos and Vel

Input:

(F_{k}), k = 0 : N

//Command

for

k = 0 : N - 1

do

p_{k} \leftarrow ComputeMomentum (q_{k}, {\dot{q}}_{k})

//Implicit Sol For Next Pos

q_{k + 1} \leftarrow SolveForNextPos (q_{k}, p_{k}, F_{k}^{-})

//Propagate Momentum Eq

p_{k + 1} \leftarrow ComputeNextMomentum (q_{k}, q_{k + 1}, F_{k}^{+})

//Find Vel from Momentum

{\dot{q}}_{k + 1} \leftarrow ExtractVelFromMomentum (q_{k + 1}, p_{k + 1})

end

Output: History of States

(q_{k}, {\dot{q}}_{k}), k = 0 : N

3 Single Multi-Rotor With a Flexible Cable Suspended Payload

In this section, the modeling of a multirotor with a flexible cable-suspended payload is performed as a multibody dynamics problem, similar to Refs. 6 and 16. To simplify the modeling of the actual flexible suspension cable, following assumptions are made accordingly.

Assumption 1.

The cable is treated as the serial connection of multiple straight links with concentrated corresponding masses at the center of each link.
The cable segments are inelastic, and they can carry the tension force only.
The cable links do not twist around their own length axis.

Figure 1 illustrates the multirotor with the suspended payload system where the suspension cable is attached to the center of gravity of the multirotor.

Fig. 1

View large Download slide

Illustration of a single multirotor with a flexibly suspended payload

Remark 1. The model shown in Fig. 1 is validated by comparing the results with the analytical catenary equations. It is found that as the number of the cable segment increases, the static curve formed by the N-link cable, once its free end is supported, asymptotically converges to the catenary shape computed by the analytical equations.

The position of the multirotor and the payload with respect to an inertial frame is denoted by

X_{Q}, X_{p} \in ℝ^{3}

⁠, respectively. The position of the mass elements on each link and the multirotor itself are written in terms of the position of the payload and the unit directions along with each link, respectively, as

\begin{matrix} X_{i} = X_{p} + \sum_{a = 1}^{i - 1} l_{a} q_{a} + \frac{1}{2} l_{i} q_{i} \\ X_{Q} = X_{p} + \sum_{i = 1}^{n} l_{i} q_{i} \end{matrix}

where q_i is the unit vector along the ith link through the cable and it lies on a manifold defined by a product of two spheres $S^{2} = {q \in ℝ^{3} | | | q | | = 1}$ ⁠, which is a part of serially connected two-sphere manifolds forming the complete cable. The links are assumed to be identical with length of l_i = l and the mass of m_i for $i = 1, 2, . ., n$ ⁠. Moreover, the payload is modeled as a point mass with the magnitude of m_p.

The rate of change of unit vectors are expressed by a cross product as ${\dot{q}}_{i} = ω_{i} \times q_{i}$ ⁠, where ω_i is the angular velocity of the ith link in its body-fixed frame. Considering the nontwisting cable link assumption made in Assumption 1, the orthogonality constraints on the link attitudes and angular velocities are introduced as $ω_{i} \cdot q_{i} = 0$ and $q_{i} \cdot {\dot{q}}_{i} = 0$ ⁠.

The mass and the moment of inertia of the multirotor are denoted by m_Q and $J_{Q} \in ℝ^{3 \times 3}$ ⁠, respectively, and its attitude is given by $R_{Q} \in SO (3)$ with the angular velocity of Ω_Q defined in the body-fixed coordinate frame. Thus, the configuration manifold of the complete system is $ℝ^{3} \times SO (3) \times {(S^{2})}^{n}$ where ${(S^{2})}^{n} = S^{2} \times \dots \times S^{2}$ ⁠.

3.1 Continuous Dynamics.

To construct the Lagrangian of the system, kinetic and potential energies are found as follows:

T = \frac{1}{2} {\dot{X}}_{Q}^{T} m_{Q} {\dot{X}}_{Q} + \frac{1}{2} Ω_{Q}^{T} J_{Q} Ω_{Q} + \frac{1}{2} \sum_{i = 1}^{n} {\dot{X}}_{i}^{T} m_{i} {\dot{X}}_{i} + \frac{1}{2} {\dot{X}}_{p}^{T} m_{p} {\dot{X}}_{p}

(9)

U = - m_{Q} g X_{Q} \cdot e_{3} - \sum_{i = 1}^{n} m_{i} g X_{i} \cdot e_{3} - m_{p} g X_{p} \cdot e_{3}

(10)

where e₃ is the unit direction along the gravitational acceleration, g, i.e., $e_{3} = {[0 0 1]}^{T}$ ⁠.

Substituting X_i and X_Q into Eqs. and , Lagrangian is obtained as

\begin{matrix} L (X_{p}, {\dot{X}}_{p}, R_{Q}, Ω_{Q}, q_{1}, {\dot{q}}_{1}, . ., q_{n}, {\dot{q}}_{n}) = T - U \\ = \frac{1}{2} {\dot{X}}_{p}^{T} M_{00} {\dot{X}}_{p} + \frac{1}{2} Ω_{Q}^{T} J_{Q} Ω_{Q} + M_{00} g X_{p} \cdot e_{3} \\ + {\dot{X}}_{p}^{T} \sum_{i = 1}^{n} M_{0 i} {\dot{q}}_{i} + \frac{1}{2} \sum_{i, j = 1}^{n} {\dot{q}}_{i}^{T} M_{i j} {\dot{q}}_{j} + \sum_{i = 1}^{n} M_{0 i} g q_{i} \cdot e_{3} \end{matrix}

(11)

where the mass terms are defined by

\begin{matrix} M_{00} = m_{Q} + \sum_{i = 1}^{n} m_{i} + m_{p} \\ M_{0 i} = (\frac{2 (n - i) + 1}{2} l m_{i} + l m_{Q}) \\ M_{i j} = {\begin{matrix} \frac{4 (n - i) + 1}{4} l^{2} m_{i} + l^{2} m_{Q}, if i = j \\ \frac{2 (n - a) + 1}{2} l^{2} m_{i} + l^{2} m_{Q}, if i \neq j (a = \max (i, j)) \end{matrix} \end{matrix}

Using the variational principle on the Lagrangian, we obtain

\begin{matrix} δ L = D_{X_{p}} L \cdot δ X_{p} + D_{{\dot{X}}_{p}} L \cdot δ {\dot{X}}_{p} + D_{R_{Q}} L \cdot δ R_{Q} \\ + D_{Ω_{Q}} L \cdot δ Ω_{Q} + \sum_{i = 1}^{n} (D_{q_{i}} L \cdot δ q_{i} + D_{{\dot{q}}_{i}} L \cdot δ {\dot{q}}_{i}) \end{matrix}

(12)

Variations on the states are derived as follows [16]:

\begin{matrix} δ R = R \hat{η} \\ δ \hat{Ω} = {(δ R)}^{T} \dot{R} + R^{T} δ \dot{R} = \hat{\dot{η}} + \hat{Ω} \hat{η} - \hat{η} \hat{Ω} \\ = \hat{\dot{η}} + \hat{Ω \times η} = \hat{\dot{η}} + {ad}_{\hat{Ω}} \hat{η} \\ δ q_{i} = ξ_{i} \times q_{i} \\ δ {\dot{q}}_{i} = {\dot{ξ}}_{i} \times q_{i} + ξ_{i} \times {\dot{q}}_{i} \end{matrix}

where

η, ξ \in ℝ^{3}

⁠. The hat operator

\hat{(\cdot)} : ℝ^{3} \to s o (3)

is the Lie algebra isomorphism between a vector

\in ℝ^{3}

and its skew-symmetric form

\in s o (3)

as

\hat{a} = \hat{[\begin{matrix} a_{x} \\ a_{y} \\ a_{z} \end{matrix}]} = [\begin{matrix} 0 & - a_{z} & a_{y} \\ a_{z} & 0 & - a_{x} \\ - a_{y} & a_{x} & 0 \end{matrix}]

The adjoint action ad $_{x} y : g \to g$ ⁠, where $g$ is a Lie algebra over some arbitrary field and $x, y \in g$ ⁠, is the group homomorphism providing a linear map as ad $_{x} y = [x, y] = x y - y x$ ⁠. The coadjoint action ad $_{x}^{*} y : g^{*} \to g^{*}$ is the dual of the adjoint action defined as ad $_{x}^{*} y = [y, x] = y x - x y$ ⁠. Hence, for the $s o (3)$ ⁠, ad $_{\hat{a}} \hat{b} = \hat{a} \hat{b} - \hat{b} \hat{a} = \hat{a \times b}$ and ad $_{\hat{a}}^{*} \hat{b} = \hat{b} \hat{a} - \hat{a} \hat{b} = \hat{b \times a}$ ⁠.

Equations below express the variations on the attitude states of multirotor and the cable links.

D_{R} L \cdot δ R = D_{R} L \cdot (R \hat{η}) = R^{T} D_{R} L \cdot \hat{η}

(13)

\begin{matrix} D_{Ω} L \cdot δ Ω = D_{Ω} L \cdot (\dot{η} + Ω \times η) \\ = D_{Ω} L \cdot \dot{η} - (Ω \times D_{Ω} L) \cdot η \\ D_{q_{i}} L \cdot δ q_{i} = (q_{i} \times D_{q_{i}} L) \cdot ξ_{i} \\ D_{{\dot{q}}_{i}} L \cdot δ {\dot{q}}_{i} = (q_{i} \times D_{{\dot{q}}_{i}} L) \cdot {\dot{ξ}}_{i} + ({\dot{q}}_{i} \times D_{{\dot{q}}_{i}} L) \cdot ξ_{i} \end{matrix}

(14)

The minimum action principle is formulated for the given system as

δ B = \int_{0}^{T} [δ L + f_{I} \cdot δ X_{Q} + τ \cdot η] d t = 0

(15)

where the variation in multirotor position is $δ X_{Q} = δ X_{p} + \sum_{i = 1}^{n} l δ q_{i}$ ⁠. The total net force generated by the propellers is represented in inertial frame with $f_{I}$ ⁠, which can be obtained by utilizing the body frame representation $f_{B}$ as $f_{I} = R_{B I}^{T} f_{B}$ ⁠. $τ$ is the net torque propellers generated, which is written in multirotor's body-fixed frame.

Continuous equations of motion for this single multirotor with a flexible cable suspended payload are derived below.

D_{X_{p}} L - \frac{d}{d t} D_{{\dot{X}}_{p}} L + f_{I} = 0

(16)

q_{i} \times q_{i} \times (D_{q_{i}} L - \frac{d}{d t} D_{{\dot{q}}_{i}} L + l f_{I}) = 0

(17)

R_{Q}^{T} D_{R_{Q}} L - \frac{d}{d t} D_{{\hat{Ω}}_{Q}} L + ({ad}_{{\hat{Ω}}_{Q}}^{*} \cdot D_{{\hat{Ω}}_{Q}} L) + \hat{τ} = 0

(18)

It is worth noting that the cable twist constraint $q_{i} \cdot ξ_{i} = 0$ is used in Eq. (17) as $(q_{i} \times (D_{q_{i}} L - \frac{d}{d t} D_{{\dot{q}}_{i}} L)) \cdot ξ_{i} = (c q_{i}) \cdot ξ_{i}$ for any arbitrary constant c.

Partial derivatives of Lagrangian are given below

\begin{matrix} D_{X_{p}} L = M_{00} g e_{3} \\ D_{{\dot{X}}_{p}} L = M_{00} {\dot{X}}_{p} + \sum_{i = 1}^{n} M_{0 i} {\dot{q}}_{i} \\ D_{R_{Q}} L = 0 \\ D_{Ω_{Q}} L = J_{Q} Ω_{Q} \\ D_{q_{i}} L = M_{0 i} g e_{3} \\ D_{{\dot{q}}_{i}} L = M_{0 i} {\dot{X}}_{p} + \sum_{j = 1}^{n} M_{i j} {\dot{q}}_{j} \end{matrix}

Substituting the partial derivatives and using the relation

{\ddot{q}}_{i} = {\dot{ω}}_{i} \times q_{i} + ω_{i} \times (ω_{i} \times q_{i}) = {\dot{ω}}_{i} \times q_{i} - (ω_{i} \cdot ω_{i}) q_{i}

⁠, continuous equations of motion become

\begin{matrix} M_{00} {\dot{V}}_{p} - \sum_{i = 1}^{n} M_{0 i} {\hat{q}}_{i} {\dot{ω}}_{i} = \sum_{i = 1}^{n} (ω_{i} \cdot ω_{i}) M_{0 i} q_{i} + M_{00} g e_{3} + f_{I} \\ M_{i i} {\dot{ω}}_{i} = q_{i} \times (\sum_{j = 1 (j \neq i)}^{n} (M_{i j} (q_{j} \times {\dot{ω}}_{j}) + (ω_{j} \cdot ω_{j}) M_{i j} q_{j})) \\ + q_{i} \times (M_{0 i} g e_{3} - M_{0 i} {\dot{V}}_{p} + l f_{I}) \\ J_{Q} {\dot{Ω}}_{Q} + {\hat{Ω}}_{Q} (J_{Q} Ω_{Q}) = τ \end{matrix}

Generalized momenta of the system are found by Legendre Transformation as

P_{V} = D_{{\dot{X}}_{p}} L

(19)

\begin{matrix} P_{Ω} = D_{Ω_{Q}} L \\ P_{ω_{i}} = q_{i} \times D_{{\dot{q}}_{i}} L \end{matrix}

(20)

3.2 Discrete Dynamics.

In the first step, Lagrangian of the system is approximated by employing the method described in Eq. to the kinetic and potential energies given in Eqs. and .

L_{d} (X_{p_{k}}, X_{p_{k + 1}}, R_{k}, Ω_{k}, q_{1_{k}}, q_{1_{k + 1}}, . ., q_{n_{k}}, q_{n_{k + 1}}) = (T_{k} - U_{k}) h

One of the key advantages of the discrete dynamics over its continuous opponent is rooted from the way that the attitude dynamics of a body is represented. The attitude of a rigid body evolves in a nonlinear manifold such as $SO (3)$ for the attitude of the multirotor and $S^{2}$ for the attitude of the cable link for the cable suspended payload system. It is crucial for any numerical integration scheme to preserve the structure of the manifold to capture the actual behavior of the model. For this reason, the discrete system takes advantage of the structure preserving Lie group actions to approximate the attitude dynamics of the multirotor and the cable links.

The group action for $SO (3)$ manifold is defined as the right matrix multiplication. The discrete update map of the multirotor attitude is given by $R_{k + 1} = R_{k} F_{k}$ where $F_{k} \in SO (3)$ is an infinitesimal rotation from R_k to $R_{k + 1}$ ⁠. The group action for the configuration space in product of two spheres, S², is similarly defined but with the left matrix multiplication. The discrete update map of the cable link attitude is given by $q_{k + 1} = F_{q_{k}} q_{k}$ where $F_{q_{k}} \in SO (3)$ is an infinitesimal rotation from q_k to $q_{k + 1}$ ⁠.

It is desired to approximate the continuous rotational kinematics of the multirotor,

{\dot{R}}_{k} = R_{k} {\hat{Ω}}_{k}

⁠, in terms of discrete updates of infinitesimal rotations to preserve orthogonality and the structure of the manifold

SO (3)

during the evolution of the states. Hence, the need for additional manipulations to fix the deviations and accumulated errors on the orthogonality of the attitude representations can be eliminated since the group action in

SO (3)

preserves the manifold, i.e.,

R_{k + 1}^{T} R_{k + 1} = I_{3 \times 3} = F_{k}^{T} R_{k}^{T} R_{k} F_{k} = F_{k}^{T} F_{k}

⁠. Now, the approximation of the vehicle's angular velocity can be written as

{\hat{Ω}}_{k} = R_{k}^{T} {\dot{R}}_{k} \approx R_{k}^{T} \frac{R_{k + 1} - R_{k}}{h} = \frac{F_{k} - I_{3 \times 3}}{h}

Using the trace identity,

Ω^{T} Ω = Tr [\hat{Ω} {\hat{Ω}}^{T}]

⁠, angular kinetic energy is rewritten by

\frac{1}{2} Ω^{T} J_{Q} Ω = \frac{Tr [(I - F_{k}) J_{Q_{d}}]}{h^{2}}

(21)

where $J_{Q_{d}}$ is the nonstandard moment of inertia, which can be obtained by $J_{Q_{d}} = \frac{1}{2} Tr [J_{Q}] I_{3 \times 3} - J_{Q}$ ⁠.

Substituting Eq. , discrete Lagrangian becomes

\begin{matrix} L_{d_{k}} = L_{d} (X_{p_{k}}, X_{p_{k + 1}}, R_{k}, F_{k}, q_{1_{k}}, q_{1_{k + 1}}, . ., q_{n_{k}}, q_{n_{k + 1}}) \\ = \frac{1}{2 h} {(X_{p_{k + 1}} - X_{p_{k}})}^{T} M_{00} (X_{p_{k + 1}} - X_{p_{k}}) + \frac{1}{h} Tr [(I - F_{k}) J_{Q_{d}}] \\ + \frac{1}{2 h} \sum_{i, j = 1}^{n} {(q_{i_{k + 1}} - q_{i_{k}})}^{T} M_{i j} (q_{j_{k + 1}} - q_{j_{k}}) \\ + \frac{{(X_{p_{k + 1}} - X_{p_{k}})}^{T}}{h} \sum_{i = 1}^{n} M_{0 i} (q_{i_{k + 1}} - q_{i_{k}}) \\ + \frac{h}{2} (M_{00} g (X_{p_{k}} + X_{p_{k + 1}}) + \sum_{i = 1}^{n} M_{0 i} g (q_{i_{k}} + q_{i_{k + 1}})) \cdot e_{3} \end{matrix}

Secondly, the variation on the discrete Lagrangian is assessed with

\begin{matrix} δ L_{d_{k}} = D_{X_{p_{k}}} L_{d_{k}} \cdot δ X_{p_{k}} + D_{X_{p_{k + 1}}} L_{d_{k}} \cdot δ X_{p_{k + 1}} + D_{R_{k}} L_{d_{k}} \cdot δ R_{k} \\ + D_{F_{k}} L_{d_{k}} \cdot δ F_{k} + \sum_{i = 1}^{n} (D_{q_{i_{k}}} L_{d_{k}} \cdot δ q_{i_{k}} + D_{q_{i_{k + 1}}} L_{d_{k}} \cdot δ q_{i_{k + 1}}) \end{matrix}

and the variation on the states are expressed below

\begin{matrix} δ R_{k} = R_{k} {\hat{η}}_{k} \\ δ F_{k} = {(δ R_{k})}^{T} R_{k + 1} + R_{k}^{T} R_{k + 1} {\hat{η}}_{k + 1} \\ = - {\hat{η}}_{k} F_{k} + F_{k} {\hat{η}}_{k + 1} = F_{k} (- F_{k}^{T} {\hat{η}}_{k} F_{k} + {\hat{η}}_{k + 1}) \\ = F_{k} (- {Ad}_{F_{k}^{T}} \cdot {\hat{η}}_{k} + {\hat{η}}_{k + 1}) \\ δ q_{i_{k}} = ξ_{i_{k}} \times q_{i_{k}} \end{matrix}

where the Adjoint group operations in $SO (3)$ are defined as ${Ad}_{F_{k}^{T}} \cdot {\hat{η}}_{k} = F_{k}^{T} {\hat{η}}_{k} F_{k}$ and ${Ad}_{F_{k}^{T}}^{*} \cdot {\hat{η}}_{k} = F_{k} {\hat{η}}_{k} F_{k}^{T}$ ⁠.

The terms regarding the variation on the attitude of the multirotor can be expressed as follows:

\begin{matrix} D_{R_{k}} L_{d_{k}} \cdot δ R_{k} = D_{R_{k}} L_{d_{k}} \cdot (R_{k} R_{k}^{T} δ R_{k}) = R_{k}^{T} D_{R_{k}} L_{d_{k}} \cdot (R_{k}^{T} δ R_{k}) \\ = R_{k}^{T} D_{R_{k}} L_{d_{k}} \cdot {\hat{η}}_{k} \\ D_{F_{k}} L_{d_{k}} \cdot δ F_{k} = D_{F_{k}} L_{d_{k}} \cdot (F_{k} F_{k}^{T} δ F_{k}) = F_{k}^{T} D_{F_{k}} L_{d_{k}} \cdot (F_{k}^{T} δ F_{k}) \\ = F_{k}^{T} D_{F_{k}} L_{d_{k}} \cdot (- {Ad}_{F_{k}^{T}} \cdot {\hat{η}}_{k} + {\hat{η}}_{k + 1}) \\ = - {Ad}_{F_{k}^{T}}^{*} \cdot F_{k}^{T} D_{F_{k}} L_{d_{k}} \cdot {\hat{η}}_{k} + F_{k}^{T} D_{F_{k}} L_{d_{k}} \cdot {\hat{η}}_{k + 1} \end{matrix}

Minimum action principle is formulated with an action sum, including the forces and moments propellers exerted to the system.

\begin{matrix} δ B_{k} = \sum_{k = 0}^{N - 1} [δ L_{d_{k}} + δ W_{d_{k}}] = 0 \\ δ W_{d_{k}} = f_{I k}^{-} \cdot (δ X_{p_{k}} + \sum_{i = 1}^{n} l δ q_{i_{k}}) \\ + f_{I k}^{+} \cdot (δ X_{p_{k + 1}} + \sum_{i = 1}^{n} l δ q_{i_{k + 1}}) + τ_{k}^{-} \cdot η_{k} + τ_{k}^{+} \cdot η_{k + 1} \end{matrix}

Finally, considering the vanishing end-point conditions, forced discrete Euler–Lagrange equations are obtained as

D_{X_{p_{k}}} L_{d_{k}} + D_{X_{p_{k}}} L_{d_{k - 1}} + f_{I k - 1}^{+} + f_{I k}^{-} = 0

(22)

\begin{matrix} R_{k}^{T} D_{R_{k}} L_{d_{k}} - {Ad}_{F_{k}^{T}}^{*} \cdot F_{k}^{T} D_{F_{k}} L_{d_{k}} \\ + F_{k - 1}^{T} D_{F_{k - 1}} L_{d_{k - 1}} + {\hat{τ}}_{k - 1}^{+} + {\hat{τ}}_{k}^{-} = 0 \\ q_{i_{k}} \times q_{i_{k}} \times (D_{q_{i_{k}}} L_{d_{k}} + D_{q_{i_{k}}} L_{d_{k - 1}} + l f_{I k - 1}^{+} + l f_{I k}^{-}) = 0 \end{matrix}

(23)

where the partial derivatives of the discrete Lagrangian are derived as below.

For the position and attitude states of the multirotor

\begin{matrix} D_{X_{p_{k}}} L_{d_{k}} = - M_{00} \frac{(X_{p_{k + 1}} - X_{p_{k}})}{h} + \frac{h}{2} M_{00} g e_{3} - \sum_{i = 1}^{n} M_{0 i} \frac{(q_{i_{k + 1}} - q_{i_{k}})}{h} \\ D_{X_{p_{k}}} L_{d_{k - 1}} = M_{00} \frac{(X_{p_{k}} - X_{p_{k - 1}})}{h} + \frac{h}{2} M_{00} g e_{3} + \sum_{i = 1}^{n} M_{0 i} \frac{(q_{i_{k}} - q_{i_{k - 1}})}{h} \\ F_{k}^{T} D_{F_{k}} L_{d_{k}} = \frac{1}{h} [J_{Q_{d}} F_{k} - F_{k}^{T} J_{Q_{d}}] \\ R_{k}^{T} D_{R_{k}} L_{d_{k}} = 0 \end{matrix}

(24)

For the attitude states of the suspension cable

\begin{matrix} D_{q_{i_{k}}} L_{d_{k}} = - \frac{(X_{p_{k + 1}} - X_{p_{k}})}{h} M_{0 i} - \sum_{j = 1}^{n} M_{i j} \frac{(q_{j_{k + 1}} - q_{j_{k}})}{h} + \frac{h}{2} M_{0 i} g e_{3} \\ D_{q_{i_{k}}} L_{d_{k - 1}} = \frac{(X_{p_{k}} - X_{p_{k - 1}})}{h} M_{0 i} + \sum_{j = 1}^{n} M_{i j} \frac{(q_{j_{k}} - q_{j_{k - 1}})}{h} + \frac{h}{2} M_{0 i} g e_{3} \end{matrix}

Remark 2. Equation (24) is derived utilizing the identities $Tr [A B] = Tr [B A] = Tr [A^{T} B^{T}] = Tr [B^{T} A^{T}], A, B \in ℝ^{n \times n}$ and $y^{T} \hat{x} z = - Tr [y z^{T} \hat{x}] = < y z^{T} - z y^{T}, \hat{x} >$ ⁠.

Solution of the discrete Euler–Lagrange equations in Eqs. and requires the initial conditions to include the states at consecutive steps, i.e.,

(X_{p_{0}}, X_{p_{1}}, R_{0}, F_{0}, q_{1_{0}}, q_{1_{1}}, \dots, q_{n_{0}}, q_{n_{1}})

⁠. However, specifying consecutive initial positions is not practical in general for most dynamical simulations, and will be avoided using the discrete Legendre Transformation, which yields position and velocity pair as defined below for the derived discrete system.

P_{V_{k}} = - D_{X_{p_{k}}} L_{d_{k}} - f_{I k}^{-}

(25)

\begin{matrix} P_{{\hat{Ω}}_{k}} = - R_{k}^{T} D_{R_{k}} L_{d_{k}} + {Ad}_{F_{k}^{T}}^{*} \cdot (F_{k}^{T} D_{F_{k}} L_{d_{k}}) - {\hat{τ}}_{k}^{-} \\ P_{ω_{i_{k}}} = - q_{i_{k}} \times (D_{q_{i_{k}}} L_{d_{k}} + l f_{I k}^{-}) \end{matrix}

(26)

P_{V_{k + 1}} = D_{X_{p_{k + 1}}} L_{d_{k}} + f_{I k}^{+}

(27)

\begin{matrix} P_{{\hat{Ω}}_{k + 1}} = F_{k}^{T} D_{F_{k}} L_{d_{k}} + {\hat{τ}}_{k}^{+} \\ P_{ω_{i_{k + 1}}} = q_{i_{k + 1}} \times (D_{q_{i_{k + 1}}} L_{d_{k}} + l f_{I k}^{+}) \end{matrix}

(28)

3.3 Solution Procedure for the Discrete Dynamics.

The solution procedure consists of four sequential stages and these stages are iterated till the final time-step, k = N. The first stage involves computation of the momenta, $(P_{V_{0}}, P_{Ω_{0}}, P_{ω_{i_{0}}})$ ⁠, from the given initial states as $(X_{p_{0}}, V_{p_{0}}, R_{0}, Ω_{0}, q_{i_{0}}, ω_{i_{0}})$ ⁠. This is achieved by employing the relations derived in continuous Legendre Transformation (19) and (20). In the second stage, the computed momenta are used along with the states at the current time-step to provide $(X_{p_{1}}, R_{1}, q_{i_{1}})$ by implicitly solving Eqs. (25) and (26). Afterwards, substituting the consecutive states at both time steps into Eqs. (27) and (28), the momenta at the next time-step, $(P_{V_{1}}, P_{Ω_{1}}, P_{ω_{i_{1}}})$ ⁠, are computed in the third stage. Finally, the last stage returns the velocity states, $(V_{p_{1}}, Ω_{1}, ω_{i_{1}})$ ⁠, utilizing the same relations in Eqs. (19) and (20) inversely, but by substituting $(P_{V_{1}}, P_{Ω_{1}}, P_{ω_{i_{1}}})$ ⁠. Thereby, the iteration repeats these sequential stages to propagate the states of the system in discrete time steps.

4 Cooperative Multi-Rotors With a Suspended Payload Via Flexible Cables

In this section, the procedure explained in Sec. 3 is extended to involve cooperative multirotors transporting a common payload that are suspended via flexible cables. Figure 2 illustrates the proposed configuration.

Fig. 2

View large Download slide

Illustration of cooperative multirotors with a flexibly suspended payload

In this configuration, the position vectors of each multirotor and their corresponding cable segments are represented in terms of the payload position, X_p, and the unit vectors along with cable links.

\begin{matrix} X_{i}^{a} = X_{p} + \sum_{r = 1}^{i - 1} l q_{r}^{a} + \frac{1}{2} l q_{i}^{a} \\ X_{Q}^{a} = X_{p} + \sum_{i = 1}^{n} l q_{i}^{a} \end{matrix}

where $q_{i}^{a}$ expresses the unit vector along the ith link of the ath multirotor. Without losing generality, suspension cables are assumed to be identical and consisting of the same number of segments.

4.1 Continuous Dynamics.

Lagrangian of the cooperative multirotors with a slung-load system is constructed as

\begin{matrix} L (X_{p}, {\dot{X}}_{p}, R_{Q}^{1}, Ω_{Q}^{1}, R_{Q}^{2}, Ω_{Q}^{2}, . ., R_{Q}^{m}, Ω_{Q}^{m}, q^{1}, q^{2}, . ., q^{m}) \\ = \frac{1}{2} {\dot{X}}_{p}^{T} M_{00} {\dot{X}}_{p} + M_{00} g X_{p} \cdot e_{3} + \sum_{a = 1}^{m} \sum_{i = 1}^{n} M_{0 i}^{a} g q_{i}^{a} \cdot e_{3} \\ + \sum_{a = 1}^{m} [\frac{1}{2} Ω_{Q}^{a}^{T} J_{Q}^{a} Ω_{Q}^{a} + {\dot{X}}_{p}^{T} \sum_{i = 1}^{n} M_{0 i}^{a} {\dot{q}}_{i}^{a} + \frac{1}{2} \sum_{i, j = 1}^{n} {\dot{q}}_{i}^{a T} M_{i j}^{a} {\dot{q}}_{j}^{a}] \end{matrix}

where

q^{a} = [q_{1}^{a}, {\dot{q}}_{1}^{a}, . ., q_{n}^{a}, {\dot{q}}_{n}^{a}]

and the mass terms are

\begin{matrix} M_{00} = \sum_{a = 1}^{m} (m_{Q}^{a} + \sum_{i = 1}^{n} m_{i}^{a}) + m_{p} \\ M_{0 i}^{a} = (\frac{2 (n - i) + 1}{2} l m_{i}^{a} + l m_{Q}^{a}) \\ M_{i j}^{a} = {\begin{matrix} \frac{4 (n - i) + 1}{4} l^{2} m_{i}^{a} + l^{2} m_{Q}^{a}, if i = j \\ \frac{2 (n - a) + 1}{2} l^{2} m_{i}^{a} + l^{2} m_{Q}^{a}, if i \neq j (a = \max (i, j)) \end{matrix} \end{matrix}

The variational principle on Lagrangian gives us

\begin{matrix} δ L = D_{X_{p}} L \cdot δ X_{p} + D_{{\dot{X}}_{p}} L \cdot δ {\dot{X}}_{p} + \sum_{a = 1}^{m} [D_{R_{Q}^{a}} L \cdot δ R_{Q}^{a} + D_{Ω_{Q}^{a}} L \cdot δ Ω_{Q}^{a}] \\ + \sum_{a = 1}^{m} [D_{Ω_{Q}^{a}} L \cdot δ Ω_{Q}^{a} + \sum_{i = 1}^{n} (D_{q_{i}^{a}} L \cdot δ q_{i}^{a} + D_{{\dot{q}}_{i}^{a}} L \cdot δ {\dot{q}}_{i}^{a})] \end{matrix}

The variation on the work done by the forces and moments on each multirotor is

δ W = \sum_{a = 1}^{m} [f_{I}^{a} \cdot (δ X_{p} + \sum_{i = 1}^{n} l δ q_{i}^{a}) + τ^{a} \cdot η_{Q}^{a}]

(29)

Using the relations given in Eqs. and , the variation on the action integral,

δ B = \int_{0}^{T} (δ L + δ W) d t = 0

⁠, is achieved as below.

\begin{matrix} δ B = \int_{0}^{T} [(D_{X_{p}} L - \frac{d}{d t} D_{{\dot{X}}_{p}} L + \sum_{a = 1}^{m} f_{I}^{a}) \cdot δ X_{p} \\ + \sum_{a = 1}^{m} \sum_{i = 1}^{n} q_{i}^{a} \times (D_{q_{i}^{a}} L - \frac{d}{d t} D_{{\dot{q}}_{i}^{a}} L + l f_{I}^{a}) \cdot ξ_{i}^{a} \\ + \sum_{a = 1}^{m} [R_{Q}^{a}^{T} D_{R_{Q}^{a}} L - \frac{d}{d t} D_{{\hat{Ω}}_{Q}^{a}} L + ({ad}_{{\hat{Ω}}_{Q}^{a}}^{*} \cdot D_{{\hat{Ω}}_{Q}^{a}} L) + {\hat{τ}}^{a}] \cdot {\hat{η}}_{Q}^{a}] d t \end{matrix}

which provides the continuous equations of motion for this cooperative system as follows:

\begin{matrix} M_{00} {\dot{V}}_{p} - \sum_{a = 1}^{m} \sum_{i = 1}^{n} M_{0 i}^{a} (q_{i}^{a} \times {\dot{ω}}_{i}^{a}) \\ = \sum_{a = 1}^{m} \sum_{i = 1}^{n} (ω_{i}^{a} \cdot ω_{i}^{a}) M_{0 i}^{a} q_{i}^{a} + M_{00} g e_{3} + \sum_{a = 1}^{m} f_{I}^{a} \end{matrix}

(30)

\begin{matrix} M_{i i}^{a} {\dot{ω}}_{i}^{a} = q_{i}^{a} \times \sum_{j = 1 (j \neq i)}^{n} (M_{i j}^{a} (q_{j}^{a} \times {\dot{ω}}_{j}^{a}) + (ω_{j}^{a} \cdot ω_{j}^{a}) M_{i j}^{a} q_{j}^{a}) \\ + q_{i}^{a} \times (M_{0 i}^{a} g e_{3} - M_{0 i}^{a} {\dot{V}}_{p} + l f_{I}^{a}) \end{matrix}

(31)

J_{Q}^{a} {\dot{Ω}}_{Q}^{a} + {\hat{Ω}}_{Q}^{a} (J_{Q}^{a} Ω_{Q}^{a}) = τ^{a}

(32)

One can deduce from Eqs. (30)–(32) that the motion of the payload is affected by both the forces applied by multirotors and the dynamics of each flexible cable connected to it. On the other hand, Eq. (31) indicates that each distinct suspension cable and the multirotor subsystem has its own separate dynamics, which are coupled by the whole system through the motion of the payload itself.

Generalized momentum states of the systems are derived by Legendre Transformation as

P_{V_{p}} = D_{{\dot{X}}_{p}} L

(33)

\begin{matrix} P_{Ω_{Q}^{a}} = D_{Ω_{Q}^{a}} L \\ P_{ω_{i}^{a}} = q_{i}^{a} \times D_{{\dot{q}}_{i}^{a}} L \end{matrix}

(34)

These momentum states will be used in discrete Hamiltonian formalism in Sec. 4.2.

4.2 Discrete Dynamics.

Discrete Lagrangian for the cooperative system is constructed as below

\begin{matrix} L_{d_{k}} = L_{d} (X_{p_{k}}, X_{p_{k + 1}}, ξ_{k}^{1}, ξ_{k}^{2}, . ., ξ_{k}^{m}) \\ = \frac{1}{2 h} {(X_{p_{k + 1}} - X_{p_{k}})}^{T} M_{00} (X_{p_{k + 1}} - X_{p_{k}}) \\ + \frac{h}{2} M_{00} g (X_{p_{k}} + X_{p_{k + 1}}) \cdot e_{3} \\ + \sum_{a = 1}^{m} [\frac{1}{h} Tr [(I - F_{k}^{a}) J_{Q_{d}}^{a}] + \frac{h}{2} \sum_{i = 1}^{n} M_{0 i}^{a} g (q_{i_{k}}^{a} + q_{i_{k + 1}}^{a}) \cdot e_{3} \\ + \frac{1}{2 h} \sum_{i, j = 1}^{n} {(q_{i_{k + 1}}^{a} - q_{i_{k}}^{a})}^{T} M_{i j}^{a} (q_{j_{k + 1}}^{a} - q_{j_{k}}^{a}) \\ + \frac{{(X_{p_{k + 1}} - X_{p_{k}})}^{T}}{h} \sum_{i = 1}^{n} M_{0 i}^{a} (q_{i_{k + 1}}^{a} - q_{i_{k}}^{a})] \end{matrix}

where $ξ^{a} = [R_{k}^{a}, F_{k}^{a}, q_{1_{k}}^{a}, q_{1_{k + 1}}^{a}, . ., q_{n_{k}}^{a}, q_{n_{k + 1}}^{a}]$ ⁠.

Applying the variation principle on the discrete Lagrangian yields

\begin{matrix} δ L_{d_{k}} = D_{X_{p_{k}}} L_{d_{k}} \cdot δ X_{p_{k}} + D_{X_{p_{k + 1}}} L_{d_{k}} \cdot δ X_{p_{k + 1}} \\ + \sum_{a = 1}^{m} [D_{R_{k}^{a}} L_{d_{k}} \cdot δ R_{k}^{a} + D_{F_{k}^{a}} L_{d_{k}} \cdot δ F_{k}^{a} \\ + \sum_{i = 1}^{n} (D_{q_{i_{k}}^{a}} L_{d_{k}} \cdot δ q_{i_{k}}^{a} + D_{q_{i_{k + 1}}^{a}} L_{d_{k}} \cdot δ q_{i_{k + 1}}^{a})] \end{matrix}

Discrete equivalent of Eq. for the work done on the system by each multirotor is expressed as follows:

\begin{matrix} δ W_{d_{k}} = \sum_{a = 1}^{m} [f_{I k}^{a -} \cdot (δ X_{p_{k}} + \sum_{i = 1}^{n} l δ q_{i_{k}}^{a}) + τ_{k}^{a^{-}} \cdot η_{k}^{a} \\ + f_{I k}^{a +} \cdot (δ X_{p_{k + 1}} + \sum_{i = 1}^{n} l δ q_{i_{k + 1}}^{a}) + τ_{k}^{a^{+}} \cdot η_{k + 1}^{a}] \end{matrix}

Ultimately, the minimization of the action provides the discrete Euler–Lagrange equations for the cooperative system as

\begin{matrix} D_{X_{p_{k}}} L_{d_{k}} + D_{X_{p_{k}}} L_{d_{k - 1}} + \sum_{a = 1}^{m} [f_{I k}^{a -} + f_{I k - 1}^{a^{+}}] = 0 \\ R_{k}^{a^{T}} D_{R_{k}^{a}} L_{d_{k}} - {Ad}_{F_{k}^{a^{T}}}^{*} \cdot F_{k}^{a^{T}} D_{F_{k}^{a}} L_{d_{k}} \\ + F_{k - 1}^{a^{T}} D_{F_{k - 1}^{a}} L_{d_{k - 1}} + {\hat{τ}}_{k - 1}^{a^{+}} + {\hat{τ}}_{k}^{a^{-}} = 0 \\ q_{i_{k}}^{a} \times q_{i_{k}}^{a} \times (D_{q_{i_{k}}^{a}} L_{d_{k}} + D_{q_{i_{k}}^{a}} L_{d_{k - 1}} + l f_{I k}^{a -} + l f_{I k - 1}^{a +}) = 0 \end{matrix}

for $a = 1, . ., m$ ⁠.

Substituting the partial derivatives of the discrete Lagrangian, which can be derived exactly as in the single multirotor case for the given cooperative system, discrete equations of motion become

\begin{matrix} M_{00} \frac{(X_{k + 1} - 2 X_{k} - X_{k - 1})}{h} = h M_{00} g e_{3} \\ + \sum_{a = 1}^{m} (h f_{k}^{a} + \sum_{i = 1}^{n} M_{0 i}^{a} \frac{(- q_{i_{k + 1}}^{a} + 2 q_{i_{k}}^{a} - q_{i_{k - 1}}^{a})}{h}) \\ \frac{1}{h} {[F_{k - 1}^{a}^{T} J_{Q_{d}}^{a} - J_{Q_{d}}^{a} F_{k - 1}^{a} + F_{k}^{a} J_{Q_{d}}^{a} - J_{Q_{d}}^{a} F_{k}^{a}^{T}]}^{\lor} = h τ_{k}^{a} \\ q_{i_{k}}^{a} \times \sum_{j = 1}^{n} M_{i j}^{a} \frac{(- q_{j_{k + 1}}^{a} + 2 q_{j_{k}}^{a} - q_{j_{k - 1}}^{a})}{h} = - q_{i_{k}}^{a} \times l h f_{k}^{a} \\ + q_{i_{k}}^{a} \times [\frac{(X_{k + 1} - 2 X_{k} + X_{k - 1})}{h} M_{0 i}^{a} - h M_{0 i}^{a} g e_{3} - l h f_{k}^{a}] \end{matrix}

By applying discrete Legendre Transformation, discrete Hamiltonian form of the derived equations are obtained with discrete momenta as follows.

P_{V_{p_{k}}} = - D_{X_{p_{k}}} L_{d_{k}} - \sum_{a = 1}^{m} f_{I k}^{a -}

(35)

\begin{matrix} P_{{\hat{Ω}}_{k}^{a}} = - R_{k}^{a^{T}} D_{R_{k}^{a}} L_{d_{k}} + {Ad}_{F_{k}^{a^{T}}}^{*} \cdot (F_{k}^{a^{T}} D_{F_{k}^{a}} L_{d_{k}}) - {\hat{τ}}_{k}^{a^{-}} \\ P_{ω_{i_{k}}^{a}} = - q_{i_{k}}^{a} \times (D_{q_{i_{k}}^{a}} L_{d_{k}} + l f_{I k}^{a -}) \end{matrix}

(36)

P_{V_{p_{k + 1}}} = D_{X_{p_{k + 1}}} L_{d_{k}} + \sum_{a = 1}^{m} f_{I k}^{a +}

(37)

\begin{matrix} P_{{\hat{Ω}}_{k + 1}^{a}} = F_{k}^{a^{T}} D_{F_{k}^{a}} L_{d_{k}} + {\hat{τ}}_{k}^{a^{+}} \\ P_{ω_{i_{k + 1}}^{a}} = q_{i_{k + 1}}^{a} \times (D_{q_{i_{k + 1}}^{a}} L_{d_{k}} + l f_{I k}^{a +}) \end{matrix}

(38)

4.3 Solution Procedure for the Discrete Dynamics.

The solution procedure follows the same stages as described in the discrete dynamics of a single multirotor with a slung load system. The momenta, $(P_{V_{p_{0}}}, P_{Ω_{0}^{a}}, P_{ω_{i_{0}}^{a}})$ ⁠, are calculated from the given initial states, $(X_{p_{0}}, V_{p_{0}}, R_{0}^{a}, Ω_{0}^{a}, q_{i_{0}}^{a}, ω_{i_{0}}^{a})$ using Eqs. (33) and (34) in the first stage. Following that, the implicit solution of Eqs. (35) and (36) supplies the states at the next time-step, $(X_{p_{1}}, R_{1}^{a}, q_{i_{1}}^{a})$ ⁠. Then, the states at the consecutive steps are substituted into Eqs. (37) and (38) to get the momenta at the next time-step, $(P_{V_{p_{1}}}, P_{Ω_{1}^{a}}, P_{ω_{i_{1}}^{a}})$ ⁠, which are also used to compute the velocity states as $(V_{p_{1}}, Ω_{1}^{a}, ω_{i_{1}}^{a})$ by the use of Eqs. (33) and (34). These sequential stages form an integration flow to solve the discrete equations of the derived Hamiltonian system.

5 Numerical Simulations

In this section, numerical simulations for the systems developed in Secs. 3 and 4 are provided to illustrate the effectiveness of the variational integration scheme in these complex systems. Simulations are carried out in each scenario for both continuous and discrete dynamics, where the fourth-order fixed-step Runge–Kutta (RK4) and variable-step ODE45 solvers are selected for the integration of continuous dynamics separately. Simulation results obtained from the implementation of variational integrator are compared with the results from the selected continuous-time integrators, RK4 and ODE45, in terms of their preservation of momentum, total energy and the geometrical constraints of the systems.

5.1 Simulation of a Single Multi-Rotor With a Flexible Cable Suspended Payload System.

In this example, the response of the continuous and discrete time models of the single multirotor with a flexible cable-suspended payload system to an initial deflection of the suspension cable is analyzed to demonstrate the effectiveness of variational integrator solution comparing the fixed-step RK4 and the variable-step ODE45 integrator solutions of the continuous-time model.

The payload is initially located at the origin and it is suspended via a cable, which is constructed from 10 rigid links and it connects payload to the center of mass of the multirotor, such that the initial deflection of the cable makes 30-degree angle with the horizontal plane as illustrated in Fig. 3. Multirotor applies a vertical force to compensate for the total weight of the whole system. The system does not have an initial velocity and an acceleration except the gravity.

Fig. 3

View large Download slide

Initial configuration of the single multirotor with a flexibly suspended payload

Remark 3. It is known that in the case where the suspension cable is not straight initially, i.e., having relative displacements among cable links, simulated systems suffer from numerical instabilities due to the chaotic and unstable behavior of the undamped dynamics regardless of the integration method utilized. A common approach to prevent such numerical instabilities is to introduce a damping into the system, such as torsional damping between the cable links in the form of $τ_{d_{i}} = - C_{d} (ω_{i - 1 / i} - ω_{i / i + 1})$ ⁠. However, since the damping would prevent the preservation of the energy, the scenarios are simulated without any damping to illustrate the effectiveness of VI method.

Parameters used in this simulation are chosen as below.

\begin{matrix} m_{Q} = 1.2 kg, m_{p} = 0.5 kg, m_{i} = 0.005 kg \\ J_{Q} = diag (1.367 e - 2, 1.367 e - 2, 2.586 e - 2) \\ h = 0.001 s (for both RK 4 and VI) \\ l_{cable} = 1 m (total cable length) \end{matrix}

Simulation results are plotted in Fig. 4. The left and middle columns correspond to the continuous-time model solution by RK4 and ODE45 integrators, respectively, whereas the right column shows the discrete-time model solution by Lie group VI method. Since the system has zero net force in all axes, it is expected to undergo a stable oscillatory motion like a dumbbell due to the restoring moment around the center of gravity, which is resulting from the vertically applied force by the propellers. The time histories of the systems' total energies are plotted in Figs. 4(a)–4(c). Although the total energies follow a sinusoidal wave in all methods, as expected, an artificial change in system's total energy is noticed over time with the RK4 and ODE45 solutions, while the variational integrator successfully preserves the total energy over a thousand second without any sign of artificial variations.

Fig. 4

View large Download slide

Comparison between the simulations by RK4, ODE45, and VI methods in terms of simulation accuracy for a single multirotor with a flexible cable-suspended payload: (a) Total energy of the system (RK4). (b) Total energy of the system (ODE45). (c) Total energy of the system (VI). (d) Linear momentum of Payload (RK4). (e) Linear momentum of Payload (ODE45). (f) Linear momentum of Payload (VI). (g) Angular momentum of the cable links (RK4). (h) Angular momentum of the cable links (ODE45). (i) Angular momentum of the cable links (VI). (j) Unit length constraints of the cable links (RK4). (k) Unit length constraints of the cable links (ODE45). (l) Unit length constraints of the cable links (VI).

Figures 4(d)–4(f) show the generalized linear momentum of the payload during the motion. As the system rotates around the center of gravity, linear momentum should be preserved, which corresponds to zero linear momentum since the system initially starts from the rest. The VI method solution on the right graph displays the anticipated behavior with a relatively small numerical tolerance, whereas RK4 solution has a clear separation on the z component of the linear momentum from the origin. Similarly, the ODE45 solution demonstrates a constant drift apart from the numerical inaccuracy although it is significantly smaller than the deviation in RK4. It can be inferred from these graphs that the conservation of linear momentum of the payload is achieved only in variational integrator solution. Likewise, the angular momentum of the cable links are plotted in Figs. 4(g)–4(i). In contrast to linear momentum results, the angular momentum of the links are preserved by each method. Since the relative motion between the links are not substantial for the desired initial configuration of the system, the error in the angular motion is expected to be small.

Finally, Figs. 4(j)–4(l) indicate how well the constraint of the constant cable length, which is enforced by the norm of the unit vectors, is satisfied over time. The deviation on the unit vector norm would also cause the warping on the configuration manifold, which is essentially what we see in the RK4 and ODE45 results. The consistent drift from the constraints affects the symplecticity and the other invariant properties on the system for the continuous-time model solution. The use of Lie group operations proves to be a substantial improvement for preserving the structure and constraints of the system as seen in VI method results, which seems to be affected by only the numerical truncation errors.

The resultant trajectories of multirotor and payload obtained from the simulations are plotted in Fig. 5. Figure 5(a) evidently illustrates the effect of accumulation of the errors, especially in vertical position, over time on the simulated trajectories during the numerical integration process. The adaptive variable-step solution of ODE45 noticeably improves the solution in Fig. 5(b) even though it does not eliminate the drift completely. Conversely, Fig. 5(c) exhibits an achievement of a long-term stable response of the single multirotor with a flexible cable suspended payload system.

Fig. 5

View large Download slide

Comparison between the simulated trajectories of the single multirotor with a flexible cable suspended payload system by RK4, ODE45, and VI integration methods: (a) Multi-rotor and payload trajectories (RK4). (b) Multi-rotor and payload trajectories (ODE45). (c) Multi-rotor and payload trajectories (VI).

5.2 Simulation of a Cooperative Multi-Rotors With a Suspended Payload System Via Flexible Cables.

This example involves three multirotors cooperatively suspending a common payload via flexible cables as formulated in Sec. 4. The simulation scenario is extending the previous example for the three vehicles where each multirotor is 120 degrees apart from the others and supports a common payload located at the origin initially as illustrated in Fig. 6.

Fig. 6

View large Download slide

Initial configuration of the cooperative multirotors with a flexibly suspended payload

The parameters that are used in this simulation are identical to the previous example, except for the cable links. The continuous flexible cables of 1 meter in length are represented with 5 links with a mass of 0.01 kg for each link. Similarly, the step-size of 0.001 s is chosen for both continuous-time and discrete-time solvers.

The expected behavior from the simulation of the cooperative system slightly differs from the results of the single multirotor and payload system. In this configuration, since the multirotors are evenly spread, the motion of the payload is constrained to the z-axis only. Considering each multirotor and the corresponding suspension cable as an individual subsystem, one can see that each subsystem again has a dumbbell-like oscillatory motion around the center of gravity of the whole system, which is a constant location on the z axis above the payload.

Remark 4. It is worth noting here that due to the coupled constraints between the individual subsystems, the numerical inaccuracies of simulation are significantly amplified regardless of the integration method.

Figure 7 presents the result of the simulation for the cooperative system. Figures on the left column correspond to the solution by RK4 integrator, where we see a quick divergent behavior after around 400 s. It can be deduced that for the given parameters and the configuration, the continuous-time model solution by RK4 does not preserve any of the system's invariants, which inherently yields inaccurate simulation results. In the middle column, the adaptive variable-step ODE45 solver slows down the rate of deviation from the actual states, but it cannot eliminate the main factor contributing the fast divergence since the continuous-time equations of motion do not strongly enforce the constraint and structure preservation. However, the figures on the right column, corresponding to the Lie group VI method, demonstrate the predicted behavior for over a long time. Despite that the VI method also suffers from the chaotic nature of the model and the amplified numerical inaccuracies, it manages to maintain the acceptable level of accuracy on the preservation of total energy, the momentum and the constraints. It can also be inferred from these results that one can minimize the artificial control command effort to mitigate the numerical errors introduced by the solver instabilities by just adopting a suitable numerical solver that complies with the physical properties and geometrical constraints of the actual system such as the family of variational integrators.

Fig. 7

View large Download slide

Comparison between the simulations by RK4 and VI methods in terms of simulation accuracy for three cooperative multirotors with a suspended payload via flexible cables. (a) Total energy of the system (RK4). (b) Total energy of the system (ODE45). (c) Total energy of the system (VI). (d) Linear momentum of Payload (RK4). (e) Linear momentum of Payload (ODE45). (f) Linear momentum of Payload (VI). (g) Angular momentum of the cable links (RK4). (h) Angular momentum of the cable links (ODE45). (i) Angular momentum of the cable links (VI). (j) Unit length constraints of the cable links (RK4). (k) Unit length constraints of the cable links (ODE45). (l) Unit length constraints of the cable links (VI).

The payload trajectories are plotted in Fig. 8. Figure 8(c) clearly concludes in favor of the VI method for the solution of cooperative multirotors with a flexibly suspended payload system.

Fig. 8

View large Download slide

Comparison between the simulated trajectories of the cooperative multirotors with a flexible cable-suspended payload system by RK4, ODE45, and VI integration methods: (a) Payload trajectory (RK4). (b) Payload trajectory (ODE45). (c) Payload trajectory (VI).

6 Conclusion

In this work, exponentially long-time numerically stable simulations of a single and cooperative multirotors transporting a suspended payload via flexible cables have been put forward by implementing a momentum and structure-preserving Lie group VI. The example simulations of both single and cooperative systems, where a fixed-step RK4 and an adaptive variable-step integrator are used for the simulation of continuous-time models, were compared in terms of their accuracy on the preservation of momentum, energy and the geometrical constraints of the system. It is shown that the implemented VI method is capable of maintaining the invariants and the constraints of the system over a long time, whereas, the continuous-time model simulation by RK4 and variable step method fell short on the same task with a consistent divergent behavior. Hence, for more consistent and accurate simulations of aforementioned physical systems, the effectiveness of the variational integrator approach has been demonstrated.

Funding Data

Office of Naval Research (Grant No. N00014-18-1-2215; Funder ID: 10.13039/100000006).
National Science Foundation (S&AS) (Grant No. 1724248; Funder ID: 10.13039/100000084).

References

1.

Sreenath

,

K.

,

Lee

,

T.

, and

Kumar

,

V.

,

2013

, “

Geometric Control and Differential Flatness of a Quadrotor Uav With a Cable-Suspended Load

,”

52nd IEEE Conference on Decision and Control

, Florence, Italy, Dec. 10–13, pp.

2269

–

2274

.10.1109/CDC.2013.6760219

Google Scholar

Crossref

2.

Sreenath

,

K.

,

Michael

,

N.

, and

Kumar

,

V.

,

2013

, “

Trajectory Generation and Control of a Quadrotor With a Cable-Suspended Load—A Differentially-Flat Hybrid System

,”

IEEE International Conference on Robotics and Automation

, Karlsruhe, Germany, May 6-10, pp.

4888

–

4895

.10.1109/ICRA.2013.6631275

Google Scholar

Crossref

3.

Tang

,

S.

, and

Kumar

,

V.

,

2015

, “

Mixed Integer Quadratic Program Trajectory Generation for a Quadrotor With a Cable-Suspended Payload

,”

IEEE International Conference on Robotics and Automation

, Seattle, WA, May 26–30, pp.

2216

–

2222

.10.1109/ICRA.2015.7139492

Google Scholar

Crossref

4.

Godbole

,

A. R.

, and

Subbarao

,

K.

,

2019

, “

Nonlinear Control of Unmanned Aerial Vehicles With Cable Suspended Payloads

,”

Aerosp. Sci. Technol.

,

93

, p.

105299

.10.1016/j.ast.2019.07.032

Google Scholar

Crossref

5.

Guo

,

M.

,

Gu

,

D.

,

Zha

,

W.

,

Zhu

,

X.

, and

Su

,

Y.

,

2020

, “

Controlling a Quadrotor Carrying a Cable-Suspended Load to Pass Through a Window

,”

J. Intell. Rob. Syst.

,

98

(

2

), pp.

387

–

401

.10.1007/s10846-019-01038-6

Google Scholar

Crossref

6.

Goodarzi

,

F. A.

,

Lee

,

D.

, and

Lee

,

T.

,

2014

, “

Geometric Stabilization of a Quadrotor UAV With a Payload Connected by Flexible Cable

,”

American Control Conference

, Portland, OR, June 4–6, pp.

4925

–

4930

.10.1109/ACC.2014.6859419

Google Scholar

Crossref

7.

Wu

,

G.

, and

Sreenath

,

K.

,

2014

, “

Geometric Control of Multiple Quadrotors Transporting a Rigid-Body Load

,”

53rd IEEE Conference on Decision and Control

, Los Angeles, CA, Dec. 15–17, pp.

6141

–

6148

.10.1109/CDC.2014.7040351

Google Scholar

Crossref

8.

Sreenath

,

K.

, and

Kumar

,

V.

,

2013

, “

Dynamics, Control and Planning for Cooperative Manipulation of Payloads Suspended by Cables From Multiple Quadrotor Robots

,”

rn

,

1

(

r2

), p.

r3

.

9.

Lee

,

T.

,

2014

, “

Geometric Control of Multiple Quadrotor UAVs Transporting a Cable-Suspended Rigid Body

,”

53rd IEEE Conference on Decision and Control

, Los Angeles, CA, Dec. 15–17, pp.

6155

–

6160

.10.1109/CDC.2014.7040353

Google Scholar

Crossref

10.

Goodarzi

,

F. A.

, and

Lee

,

T.

,

2015

, “

Dynamics and Control of Quadrotor UAVs Transporting a Rigid Body Connected Via Flexible Cables

,”

American Control Conference (ACC)

, Chicago, IL, July 1–3, pp.

4677

–

4682

.10.1109/ACC.2015.7172066

Google Scholar

Crossref

11.

Simo

,

J.

,

Tarnow

,

N.

, and

Wong

,

K.

,

1992

, “

Exact Energy-Momentum Conserving Algorithms and Symplectic Schemes for Nonlinear Dynamics

,”

Comput. Methods Appl. Mech. Eng.

,

100

(

1

), pp.

63

–

116

.10.1016/0045-7825(92)90115-Z

Google Scholar

Crossref

12.

Marsden

,

J. E.

, and

West

,

M.

,

2001

, “

Discrete Mechanics and Variational Integrators

,”

Acta Numer.

,

10

(

1

), pp.

357

–

514

.10.1017/S096249290100006X

13.

Hairer

,

E.

,

Wanner

,

G.

, and

Lubich

,

C.

,

2006

,

Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations

,

Springer

,

Berlin/Heidelberg

.

14.

Lee

,

T.

,

Leok

,

M.

, and

McClamroch

,

N. H.

,

2007

, “

Lie Group Variational Integrators for the Full Body Problem

,”

Comput. Methods Appl. Mech. Eng.

,

196

(

29–30

), pp.

2907

–

2924

.10.1016/j.cma.2007.01.017

15.

Manchester

,

Z. R.

, and

Peck

,

M. A.

,

2016

, “

Quaternion Variational Integrators for Spacecraft Dynamics

,”

J. Guid. Control Dyn.

,

39

(

1

), pp.

69

–

76

.10.2514/1.G001176

Google Scholar

Crossref

16.

Lee

,

T.

,

McClamroch

,

N. H.

, and

Leok

,

M.

,

2005

, “

A Lie Group Variational Integrator for the Attitude Dynamics of a Rigid Body With Applications to the 3D Pendulum

,”

Proceedings of 2005 IEEE Conference on Control Applications

, Toronto, ON, Canada, Aug. 28–31, pp.

962

–

967

.10.1109/CCA.2005.1507254

17.

Lee

,

T.

,

Leok

,

M.

, and

McClamroch

,

N. H.

,

2007

,

Discrete Control Systems

, arXiv preprint arXiv:0705.3868.

18.

Nikhilraj

,

A.

,

Simha

,

H.

, and

Priyadarshan

,

H.

,

2019

, “

Optimal Energy Trajectory Generation for a Quadrotor UAV Using Geometrically Exact Computations on se(3)

,”

IEEE Control Syst. Lett.

,

3

(

1

), pp.

216

–

221

.10.1109/LCSYS.2018.2874103

Google Scholar

Crossref

19.

Nordkvist

,

N.

, and

Sanyal

,

A. K.

,

2010

, “

A Lie Group Variational Integrator for Rigid Body Motion in SE (3) With Applications to Underwater Vehicle Dynamics

,”

49th IEEE Conference on Decision and Control (CDC)

, Atlanta, GA, Dec. 15–17, pp.

5414

–

5419

.10.1109/CDC.2010.5717622

Google Scholar

Crossref

20.

Schultz

,

J.

,

Flaßkamp

,

K.

, and

Murphey

,

T. D.

,

2017

, “

Variational Integrators for Structure-Preserving Filtering

,”

ASME J. Comput. Nonlinear Dyn.

,

12

(

2

), p.

021005

.10.1115/1.4034728

Google Scholar

Crossref

21.

Leok

,

M.

, and

Shingel

,

T.

,

2012

, “

General Techniques for Constructing Variational Integrators

,”

Front. Math. China

,

7

(

2

), pp.

273

–

303

.10.1007/s11464-012-0190-9

Google Scholar

Crossref

22.

Davis

,

P. J.

, and

Rabinowitz

,

P.

,

1984

, “

Chapter 2—Approximate Integration Over a Finite Interval

,”

Methods of Numerical Integration

, 2nd ed.,

P. J.

Davis

and

P.

Rabinowitz

, eds.,

Academic Press

, Cambridge, MA, pp.

51

–

198

.

Momentum Preserving Simulation of Cooperative Multirotors With Flexible-Cable Suspended Payload

Abstract

1 Introduction

2 Overview of the Variational Integrator Procedure

2.1 Continuous Dynamics.

2.2 Discrete Dynamics.

3 Single Multi-Rotor With a Flexible Cable Suspended Payload

3.1 Continuous Dynamics.

3.2 Discrete Dynamics.

3.3 Solution Procedure for the Discrete Dynamics.

4 Cooperative Multi-Rotors With a Suspended Payload Via Flexible Cables

4.1 Continuous Dynamics.

4.2 Discrete Dynamics.

4.3 Solution Procedure for the Discrete Dynamics.

5 Numerical Simulations

5.1 Simulation of a Single Multi-Rotor With a Flexible Cable Suspended Payload System.

5.2 Simulation of a Cooperative Multi-Rotors With a Suspended Payload System Via Flexible Cables.

6 Conclusion

Funding Data

References

Get Email Alerts

Cited By

ASME Journals

ASME Conference Proceedings

ASME eBooks

Resources

Opportunities

Momentum Preserving Simulation of Cooperative Multirotors With Flexible-Cable Suspended Payload

Abstract

1 Introduction

2 Overview of the Variational Integrator Procedure

2.1 Continuous Dynamics.

2.2 Discrete Dynamics.

3 Single Multi-Rotor With a Flexible Cable Suspended Payload

3.1 Continuous Dynamics.

3.2 Discrete Dynamics.

3.3 Solution Procedure for the Discrete Dynamics.

4 Cooperative Multi-Rotors With a Suspended Payload Via Flexible Cables

4.1 Continuous Dynamics.

4.2 Discrete Dynamics.

4.3 Solution Procedure for the Discrete Dynamics.

5 Numerical Simulations

5.1 Simulation of a Single Multi-Rotor With a Flexible Cable Suspended Payload System.

5.2 Simulation of a Cooperative Multi-Rotors With a Suspended Payload System Via Flexible Cables.

6 Conclusion

Funding Data

References

Get Email Alerts

Cited By

Related Articles

Related Proceedings Papers

Related Chapters

ASME Journals

ASME Conference Proceedings

ASME eBooks

Resources

Opportunities

This Feature Is Available To Subscribers Only