{'fraction_non_alphanumeric': 0.04163069184278689, 'fraction_numerical': 0.01009797665077518, 'mean_word_length': 5.2623413258110014, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Quaternion variational integration for inertial maneuvering in a biomimetic UAV

Biological flying, gliding, and falling creatures are capable of extraordinary forms of inertial maneuvering: free-space maneuvering based on fine control of their multibody dynamics, as typified by the self-righting reflexes of cats. However, designing inertial maneuvering capability into biomimetic robots, such as biomimetic unmanned aerial vehicles (UAVs) is challenging. Accurately simulating this maneuvering requires numerical integrators that can ensure both singularity-free integration, and momentum and energy conservation, in a strongly coupled system - properties unavailable in existing conventional integrators. In this work, we develop a pair of novel quaternion variational integrators (QVI) showing these properties, and demonstrate their capability for simulating inertial maneuvering in a biomimetic UAV showing complex multibody-dynamics coupling. Being quaternion-valued, these QVIs are innately singularity-free; and being variational, they can show excellent energy and momentum conservation properties. We explore the effect of variational integration order (left-rectangle vs. midpoint) on the conservation properties of integrator, and conclude that, in complex coupled systems in which canonical momenta may be time-varying, the midpoint integrator is required. The resulting midpoint QVI is well-suited to the analysis of inertial maneuvering in a biomimetic UAV - a feature that we demonstrate in simulation - and of other complex dynamical systems.


Introduction
Inertial maneuvering is a form of free-space orientation control that is common in biological creatures: it involves the use of biological multi-body dynamics (e.g., the effects of wing or limb inertia) to achieve orientation changes independent of aerodynamic forces. As a free-fall response, it is found in a wide range of flightless creaturesfor instance, in cats, where it is responsible for their extraordinary capability to land on their feet [1,2]; but also, in species of gecko [3,4] and insect [5]. As a flight control mechanism, it is observed in a range of flying creatures as a mechanism for ultra-low-airspeed controlfor instance, during flight initiation from an inverted state [6], slow turning [7] and vertical landing [8,9]. Both these forms of inertial maneuvering are highly relevant to the design of biomimetic robots: self-righting and self-stabilizing biomimetic terrestrial robots [10][11][12]; and to biomimetic unmanned aerial vehicles (UAVs) undergoing complex maneuvering and/or landing operations [13][14][15].
Designing and simulating complex inertial maneuvering in such systems poses several challenges. For instance, in a vertical landing maneuver in a biomimetic UAV, the UAV may pass through any given orientation, or attitude, state: for robust simulation and analysis, without the limitations of the gimbal lock phenomenon [16][17][18], a singularity-free representation of the UAV orientation is required. A singularity-free representation of choice for many aerospace systems involves the quaternions: a system of four-parameter hypercomplex numbers [19,20]. The set of normalized quaternions provide a smooth and compact method of representing orientation, and thus, tracking and simulating rotating systems. They allow singularity-free analysis of the dynamics of UAVs [21,22], manned aircraft [18,23], spacecraft and space objects [24][25][26], and a range of other structures. Unlike alternative three-parameter representations using Euler angles or Rodriguez parameters, the quaternions are singularity-free [16]; and also, as a non-singular higher-dimensional representation, they are more compact than rotation matrices, and the constraint that generates dimensional reduction is more numerically accessible: quaternion normalization [27]; instead of rotation matrix special orthogonality [28,29]. However, when utilizing quaternions for orientation representation, conventional numerical integration techniques are not directly applicable. The quaternion normalization constraintor, constraint to unit hyperspherenecessitates that numerical integration be carried out not in unconstrained space, but on this hypersphere. Formally, the set of orientation quaternions is a Lie group [30,31]. Neglect of this property, via the use of the quaternion pseudoderivative, involving finite-difference integration, with periodic normalization of the output; or via naïve treatment as a differential-algebraic equation (DAE), both can lead to a degradation of integrator quality and generality [32][33][34][35].
This degradation in quality has spurred the development of dedicated quaternion numerical integrators.
Variational integration represent a recent development in computational mechanics, showing several advantages over non-variational integration: in particular, favorable energy and momentum conservation properties [44,45].
The application to Lie groups, including quaternions, also allows an efficient, embedded, treatment of Lie group constraints, via appropriate definition of the perturbation defining the variational derivative. However, existing formulations, particularly of QVIs, have been limited to systems that are not directly applicable to the complex dynamics of biomimetic UAVsinstead, being focused on uncoupled rigid-body motion [24,43,46,47] and/or problems in orbital and spacecraft mechanics [40][41][42]. The dynamics of biomimetic UAVs (and, other robotic systems) undergoing inertial maneuvering are differ significantly: multibody dynamic systems showing complex coupling, and subjected to complex aerodynamic forces [13,48]. In addition, the degradation in integrator quality arising from simplistic treatments of quaternion normalization is particularly concerning in the context of inertial maneuvering: given that inertial maneuvering is a phenomenon that is founded on conservative inertial dynamics, rather than, e.g., non-conservative aerodynamics, it is of key importance to ensure accurate energy and momentum conservation during integration. Despite ongoing generalization efforts [24,30,41], there remain no quaternion variational integrators available for complex inertial-maneuvering systems.
In this work, we seek to remedy this deficit, with the application of quaternion variational integration to a complex and general multibody-dynamic model of inertial maneuvering in a biomimetic UAV. We develop a pair of novel QVI formulations for the complex multibody dynamic behavior of this inertial maneuvering behavior: a conventional left-rectangle variational integrator, and a higher-order midpoint variational integrator.
We study the effects of canonical momenta, finding that retaining favorable energy and momentum conservation requires midpoint integration, when canonical momenta are varying. The application of quaternion variational integration to this system relieves an analysis bottleneck associated with integration availability, and enables an analysis of complex forms of inertial maneuvering in biomimetic aircraft. The midpoint integrator that is developed is applicable not only to analogous models, e.g., of inertial maneuvering biological creatures; but to more general formulations of multibody-dynamic systems with strong translation-rotation coupling. The application of variational and higher-order quaternion integrators to the dynamics of complex aerospace systems shows significant potential for accelerating simulation-based design processes in such systems.

Multibody modelling of a biomimetic UAV
Pons and Cirak [13][14][15] studied the flight dynamics of a biomimetic morphing-wing unmanned aerial vehicle (UAV), in order to explore the complex forms of maneuverability available to such systems. The analysis of these maneuvers lends itself to a quaternion-based flight dynamic model: a wide range of orientation states needs to be modelled; along with the possibility for rapid orientation changes. In some flight maneuvers, ensuring accurate energy and momentum conservation via quaternion variational integration may not be a top concern; but in the case of biomimetic inertial maneuvering, it is paramount. Fig. 1 shows the biomimetic UAV under consideration, alongside one hypothetical form of inertial maneuvering applicable to such a system: inverted takeoff and landing, from a fixed inverted landing state, to an arbitrary heading.  [8,9] and pigeons [6]: wing inertial control enables takeoff from, and landing on, an inverted surface.
The multibody flight dynamic model of this UAV can be expressed in Lagrangian form, via the definitions of the UAV kinetic energy ( ) [13][14][15]: where ̇( ) is the UAV velocity in the body-fixed frame; ( ) is its angular velocity in the body-fixed frame; is a vector description of the UAV morphing configuration parameters; and a , a 0 , A , A , , and are coefficient functions, given in [13]. Our analysis will assume nothing regarding the form of these coefficients, except for the fact that they are continuous functions of . The form Eq. 1 is well-suited for variational integration, because the coefficients functions a , etc., are independent of the state variables, ̇( ) and ( ) , and are only dependent on the morphing configuration parameters,this reduces the extent of kinematic coupling in the model. The notation (⋅) ( ) denotes resolution in the UAV body-fixed reference frame, and (⋅) ( ) resolution in the earth reference frame. If the orientation of the UAV is represented by an orientation quaternion, , defined under the Hamilton convention [27], then the relationship between resolution in earth and body-fixed reference frames, for any 3-vector (i.e., imaginary quaternion) , may be represented: The dynamic behavior of this UAV model can be expressed via the principle of least action, asserting that the system action functional is stationary with respect to first-order perturbations [49,50]. In the context of simple inertial maneuvering in the absence of gravitational acceleration, cf. [8], the UAV Lagrangian is , and we may express the UAV dynamics as: for generalized forces , generalized coordinates r, time , and simulation time interval [ 0 , ]. The primary objective for quaternion variational integration is to accurately model the inertial maneuvering behavior: the changes in energy and momentum that arise from changes in morphing configuration ( ), in a context with strong translation-rotation coupling.

Choice of dimensional-reduction map
To undertake a variational analysis of this the UAV system, Eq. 1-3, a set of generalized coordinates and associated time-derivative variables are required. In a conventional variational analysis, e.g., with purely translational dynamics, these would be related directly via the time differentiation operator [44,45]. In the case of a quaternion generalized coordinatedefined under the Hamilton convention [27] direct time differentiation is unsuitable. The derivative ̇ is an underconstrained parameterization of the orientation rate, requiring an additional constraint in the integrator inter-step equations derived from the hypersphere constraint, ‖ ‖ = 1, for quaternion norm ‖⋅‖. Moreover the -̇ pairing introduces an additional dependency of the kinetic energy on -(̇, )via the quaternion representation of ( ) : for quaternion multiplication ⊗ and transpose (⋅) † . For properties of quaternion algebra, see [51]. The (̇, ) dependency significantly complicates the variational analysis by introducing further terms in the chain rule expansion of .
Variational treatments of the constraint on ̇ utilize the following three-step process [30]. (1) In discrete, or infinitesimal continuous, definitions of the derivative ̇, the rotation of on the hypersphere between adjacent time-steps will be small. (2) Under a small-rotation approximation, these rotations can be parameterized via a local dimensional-reduction map (diffeomorphism) generating a local, unconstrained, 3-parameter representation.
(3) Local discrete or infinitesimal rotations can be solved in this local unconstrained representation, and then mapped back to full quaternion space. (end). For unit quaternions at adjacent time-steps, +1 and ; some map = ( ) with mapping parameter ; and dynamical-system solution process (⋅), this process may be expressed: Note that a distinction exists between right ( ⊗ Δ ) and left (Δ ⊗ ) multiplication, as ⊗ is noncommutative [51].
Two diffeomorphisms, ( ), are available: the exponential map, and the Cayley map [30]. Existing QVIs have commonly used the Cayley map [24,30]: for ∈ ℝ 3 , i.e., an imaginary quaternion. However, as an alternative for systems with complex dynamics, such as our inertial maneuvering UAV, we utilize the exponential map: for ∈ ℝ 3 and quaternion exponential exp(⋅) [27]. The elegance of this approach is that the physical angular velocity variable ( ) = ⁄formally, a logarithmic derivativeis itself used as a dimensionally-reduced proxy derivative for the orientation quaternion. The continuous relation between ̇ and ( ) is given directly by Eq. 3; and the discrete relation, by any one of a range of quaternion kinematic integrators, as per [46]. The simplest is the first-order Crouch-Grossman (CG) method, simply: for timestep ℎ. Other integrators can be utilized or derived, depending on the available information regarding the variation in ( ) over the integration step.

Left-rectangle integrator
In the discrete mechanics framework of [44,45], the time integral in the principle of least action, Eq. 3, may be approximated with left-rectangle integration: The variational derivative of , , is defined with reference to first-order perturbations in the orientation ( ) and position ( ( ) ) of the UAV. The process of computing involves first defining these perturbations in and ( ) ; then propagating these perturbations to other relevant kinematic variables (̇( ) and ( ) ); and finally propagating the complete kinematic perturbations to the perturbed Lagrangian, . First, defining the perturbations. Following [24], the orientation, , is subjected to a continuous norm-preserving perturbation defined via the quaternion exponential: where ( ) represents a perturbative angular velocity axis in the body-fixed frame; that is, the axis around which the system will be perturbed by a small angle. This perturbation can be propagated directly to a variational derivative in the angular velocity, ( ) . Under a first-order (two step) approximation, the discrete analogue of Eq. 4 can be expressed: with associated variational derivative: Turning to the translation, ( ) , the perturbation definition is more straightforward, as there is no normalization constraint. We may utilize simply: Combining this translational perturbation with the angular velocity perturbation of Eq. 12, we can compute the perturbation on ̇( ) , as: Finally, for the complete system of Eq. 1, the variational derivative of the Lagrangian, , can be expressed via the chain rule as: The discrete principle of least action, Eq. 9, is then directly expressible as two coupled inter-step equations; in translation and rotation, respectively; and coupling the step at − 1 to the step at : where the generalized force is may be split into translational force ( ) , in the earth frame, and rotational moment ( ) , in the body-fixed frame.
It is worth observing these inter-step equations more closely. The equations take the form of momentum balances in the system's canonical momenta [52]; here: The translational equation is symmetric, with momentum term , . In the absence of external force, ( ) , Eq. 16 ensures that the translational momentum of the UAV is conserved in a numerical-integration context. This conservation law is of significant utility in preserving the accuracy of inertial maneuvering simulation. The rotational equation is nearly symmetric, with momentum term , : the exception is the symmetry-breaking term, ℎ̇( ) × 1, . The presence of this symmetry-breaking term means that rotational momentum, , , is not conserved within the inter-step equation: there is a problematic potential for rotational momentum drift in simulation. This potential for momentum drift cannot be eliminated by trying to choose an alternative generalized coordinate formulation (e.g., with ̇( ) , ( ) ) so as to remove the symmetry-breaking term. The indirect dependency of the system kinetic energy on generalized coordinatescontained in the proxy derivative relation, Eq. 11 -means that at least one of the system's canonical momenta will always not be conserved, as the associated generalized coordinate will not be ignorable [53]. In an uncoupled system, a particle-based axis-angle representation ( -, the generalised coordinate(s) corresponding to ( ) ) leads to conservation of angular momentum [54] though the quaternion canonical momentum is not conserved in any case [55]. However, in Eq. Finally, as is the case even in simpler dynamical systems [24], this inter-step equation, Eq. 16, does not permit an analytical solution, but must instead be solved numerically. The Newton iteration for the solution of the inter-step relation may itself be expressed in an asymmetric momentum-balance form: This completes the derivation of this left-rectangle integrator.

Midpoint integrator
The implication of the canonical-momentum analysis in Section 3.1 is that left-rectangle quaternion variational integrators, while suitable for uncoupled rigid-body dynamics [24]; are likely to be unsuitable for the complex coupling present in an inertial maneuvering simulation. In an uncoupled simulation, when the system canonical momenta are conserved, a left rectangle discretization is an exact representation of the system conservation properties. When the system canonical momenta may vary, left rectangle integration becomes a loworder approximation of the system canonical momenta behavior. A more accurate form of integration is thus one key avenue to re-attaining good conservation properties in this situation. Here we apply midpoint integration, a small, but significant increase in integration accuracyand one not utilized before in quaternion variational integration. Applying discrete midpoint integration to the principle of least action, Eq. 2, we obtain: where the tilde (̃) denotes evaluation at an inter-step midpoint. Heuristically, we would expect midpoint integration to be able to capture linear trends in the Lagrangian . The midpoint location in all the non-quaternion system variables ( , ( ) , etc.) can be computed via linear interpolation, but to compute the midpoint quaternion (̃) a different scheme is appropriate, as linear interpolation does not preserve the orientation quaternion normalization [56]. Other schemes available include normalized linear interpolation (NLERP), spherical linear interpolation (SLERP), spherical spline interpolation (SQUAD) [56] and eigenvector quaternion averaging [57].
For any choice of interpolation, the variational derivative of ̃ is: This change simplifies the manipulation, and is not fundamentalit may be verified that the equivalence between the earth and body perturbations is exact: = exp( ( ) ) = exp( ( ) † ) = exp( ( ) ).
By parameterizing the perturbed midpoint in a local perturbative direction ̃( ) ; ̃= exp(̃( ) )̃≅̃+̃( )̃+ ( 2 ), results are obtained for the generalized derivatives ̃( ) and ̇( ) : This has transformed the problem of computing ̃ into one of computing ̃( ) and ̇( ) , as a function of ( ) +1 ( ) . Again, this cannot be done exactly, but discrete approximations are significantly easier to obtain: Both these estimates are accurate to first order in ℎ. More formally, this approach has approximated the quaternion interpolation problem by an interpolation of the quaternion generators (cf. [58]). This leads to results for the variational derivatives of the proxy derivatives: and thus, with some manipulation, the integrator inter-step equations, in translation and rotation respectively: Notably, the rotational equation shows a symmetry after the manner of Eq. 16, though with cross products in ̇( ) .
The midpoint quaternion ̃ contains a dependency on ̃( ) via SLERP interpolation (Eq. 23) and +1 , integrated under the assumption of a constant angular velocity of ( ) =̃( ) over the step interval ( , + 1), corresponding again to midpoint integration: These inter-step equations are solved via Newton's method; though the dependency of ̃ on ̃( ) precludes the computation of an analytical Jacobian. A numerical Jacobian is implemented instead. This completes the derivation of the integrator.

Integrator conservation properties under test conditions
As an initial test of the conservation properties of our integrators, we first study their behavior in a UAV model with fixed wings, i.e., no wing morphing. Note that, while this test condition corresponds to the simulation of a single rigid body, it still contains translation-rotation coupling because of the use of a body-fixed reference point for the UAV dynamics [13,14]. As such, we can test the conservation behavior of our integrator in where tot is the total system mass, ̇c .o.m.
In addition, the kinetic energy conservation error can be estimated directly from Eq. 1: Kinetic energy error: Fig. 2A shows the results from the left-rectangle QVI are compared to those from midpoint QVI, and also, for output validation, an adaptive RK45 integrator in Euler angles with pole-switching, developed for this biomimetic UAV [13][14][15]. Fig. 2B illustrates the behavior of these integrators in terms of their momentum and energy conservation properties. Several points may be noted. The midpoint QVI significantly outperforms the left-rectangle QVI in energy and rotational momentum conservation: errors are consistently lower by a factor of 10 2 . However, notably, both integrators are effectively exact in translational momentum conservationthey do not reach errors above 10 −13 over 100 s. This is representative of the fact that the QVI inter-step equations in translational DOFs function as a momentum balance, whereas the inter-step equations in angular DOFs do not (Eq. 30). However, a notable feature of the midpoint QVI's angular momentum conservation is that the angular momentum error is oscillatory at near constant amplitude: Fig. 2C illustrates this error over the 100 s integration window. In Fig. 2B, this corresponds to the constant maximum angular momentum error for > 10 0 . The constant-amplitude nature of this error suggests that there may be a method of accounting for this dynamical effect exactly (e.g., via more detailed analysis of Eq. 30) and thereby reducing all momentum errors to the level of machine precision. It also provides additional assurance regarding the long-term behavior of the midpoint QVI: the only conserved property in the midpoint QVI which shows any steady trending deviation is the kinetic energy itself. Comparing the QVIs to the adaptive RK45 integrator in Euler angles: the latter allows significantly larger step sizes, and at these step sizes shows conservation properties that lie between the left-rectangle QVI and midpoint QVIexcept in translational momentum, in which these latter integrators are effectively exact. For very low-order integrators, these QVI perform well.

Integrator performance under inertial maneuvering behavior
Finally, as a test of the midpoint QVI for true inertial maneuvering, Fig. 3 shows a maneuver arising from an oscillatory pattern of UAV wingbeat motion. The symmetric wing control inputs are dihedral ( ) = sin( ) rad and incidence ( ) = −0.5 cos( ) rad. Note that, as a result of wing motion, the system coefficients are now time-varying. Following per Bergou et al. [8], we simulate a bio-inspired pitching maneuver, including aerodynamic forces, but neglecting gravitational ones. This motion could be of utility in the design of UAV vertical landing maneuvers, analogous to those carried out by bats [8,9]. Fig. 3 shows the kinematic history of the UAV over the simulation, alongside rendering of the maneuver. Results are presented for the timestep, ℎ = ℎ RK , where the mean stepsize of the adaptive RK45 integrator in Euler angles is ℎ RK = 183 ms (over a 60 s simulation). The midpoint QVI performs excellently, to the extent that the difference between QVI and RK45 results may be attributable to inaccuracy in the RK45 integrator. This demonstrates the effectiveness of the QVI for the simulation of some forms of flapping-wing flight. It additionally demonstrates the differences between the QVI and the RK45 integrator in terms of the suitable step size for a given system. In this flapping-wing flight simulation, the QVI is capable of matching and exceeding the adaptive RK45 step size; whereas in the free rotation and translation test case these integrators require smaller step sizes. The mean adaptive RK45 step size is thus not necessary a good metric for comparing the QVI performance. Further numerical analysis techniques may be able to shed light on these effects, but given the complexity of the integrators system-specific testing is likely to be a more practical approach for assessing integrator suitability. Overall, this midpoint QVI is well-suited to the analysis of complex inertial maneuvering in UAVs; and shows potential for applications in other systems showing complex dynamics.

Conclusions
In this work, we developed a pair of QVIs for simulating the strongly coupled dynamics that govern inertial maneuvering in biomimetic UAV. Formulating the inertial maneuvering dynamics in terms of a general quaternion-parameterized system with strong translation-rotation coupling, we developed integrators under two variational integration conditions: left-rectangle variational integration, as per existing QVIs, and a novel midpoint variational integration. Resolving aspects such as the appropriate definitions of quaternion and perturbation interpolation, we applied these two integrators to numerical tests of inertial maneuvering behavior in a biomimetic UAV showing complex coupling behavior, and studied their behavior. The midpoint QVI is best suited to the analysis of inertial maneuvering behavior: with the canonical momenta of the system varying due to coupling, a midpoint approximation of these momenta is significantly more precise than a conventional left-rectangle approximation. The midpoint QVI performs excellently in energy and momentum conservation, and is a candidate for precise and singularity-free simulation of practical inertial maneuvering behavior in biomimetic UAVs and other complex systems.