Modelling Spacecraft Control Systems for On-orbit Manufacturing using Quaternion-based Variational Integration

1. Introduction

The equation of motion for a mechanical system can be solved using variational integrators, which provide advantages over traditional algorithms such as Runge–Kutta methods. Instead of deriving the continuous-time differential equations of motion and then discretizing them, variational approaches begin by discretizing the Lagrangian and the action integral of the system. The tools of variational mechanics are then applied to derive discrete-time equations of motion.Integrators derived in this way retain many of the essential properties of the continuous system, such as momentum and energy conservation [1]. These integrators are symplectic, preserve energy in Hamiltonian conservative systems, and remain highly accurate for long-term integration. In addition, they are computationally efficient and stable, even for relatively large fixed time steps, making them well suited for real-time estimation and control applications.

Examples of integration methods that respect motion integrals of mechanical systems, often called geometric or symplectic integrators, have been known for decades. The classic example is the Stormer-Verlet and Galerkin Variational Integrators (GVI) [2,3]. GVI approach is non-symplectic integrators also reducing simulation time significantly for long integration periods and also noted, errors start growing proportionally to the square root of time [3]. These early methods were often devised in ad hoc ways that do not generalize well to arbitrary mechanical systems or higher orders of accuracy. More recently, systematic methods for deriving symplectic integrators using discrete-time versions of ideas from variational mechanics have been introduced [1]. These methods are straightforward applications of Lagrangian and Hamiltonian dynamics, and they can generate integration schemes of any desired order of accuracy.

Discrete variational mechanics has previously been applied to problems in optimal control and rigid-body dynamic control by using rotation matrices to parameterize attitude [4-6]. Momentum-preserving integrators have also been derived by other (non-variational) means for rigid-body dynamics using quaternions to parameterize attitude [7-11]. The primary contribution of this paper is the application of discrete variational mechanics to spacecraft attitude dynamics using quaternion state variables. The emphasis on quaternions over other attitude parameterizations here is due to both the compact and elegant derivations they enable and their prevalence in the implementation of spacecraft guidance, navigation, and control algorithms. Specifically, they lend themselves to straightforward feedback control and estimation schemes of practical relevance in flight software.

The paper proceeds as follows. Section 2 gives a brief review of quaternions and outlines the notation used throughout the paper. Section 3 derives the classical Euler equation of rigid-body dynamics in continuous time using Hamilton’s principle. Next, Section 4 presents this derivation in discrete time, leading to the variational integrator presented in Section 5. Sections 6-8 incorporate several extensions to the basic rigid-body integrator, including reaction wheels, external torques, and internal damping. Finally, in Section 9, several numerical examples are presented, including an extended Kalman filter for attitude determination.

2. Background

Attitude dynamics and rotations are parameterized with quaternions throughout this paper. A brief review of their properties is presented in this section. A more thorough treatment is given by Altmann [12].

Quaternions form an algebra with a noncommutative binary product operation. It is often convenient to think of them as four-dimensional objects composed of a three-dimensional vector v and a scalar s:

$q\mathrm{<mo>=</mo>}\left[\begin{array}{c}s\\ v\end{array}\right]\text{ (1)}$

This representation allows the quaternion product to be written in terms of scalar and vector products:

$q_1{q}_2\mathrm{<mo>=</mo>}\left[\begin{array}{c}{v}_1\mathrm{<mo>*</mo>}{v}_2+s_1{v}_2+s_2{v}_1\\ s_1{s}_2-v_1.v_2\end{array}\right]\text{ (2)}$

Note that $q_1{q}_2\ne q_2{q}_1$ . Throughout the paper, quaternion products are indicated by juxtaposition, whereas scalar and vector products are indicated in the usual way, with the . and × symbols, respectively.

Rotations can be conveniently represented by unit-length quaternions. If r is a unit vector in ${ℝ}^3$ representing the axis of rotation and θ is the angle of rotation, then the quaternion representing the rotation is as follows:

$q\mathrm{<mo>=</mo>}\left[\begin{array}{c}r\sin \mathrm{<mo\; stretchy="false">(</mo>}\frac{\theta }{2}\mathrm{<mo\; stretchy="false">)</mo>}\\ \cos \mathrm{<mo\; stretchy="false">(</mo>}\frac{\theta }{2}\mathrm{<mo\; stretchy="false">)</mo>}\end{array}\right]\text{ (3)}$

Both q and −q correspond to the same rotation, making the unit quaternions a ``double cover'' of the group of rotations.

The conjugate of a quaternion is denoted with a superscript ${}^{†}$ and represents the rotation about the same axis r by −θ:

$q^{†}\mathrm{<mo>=</mo>}\left[\begin{array}{c}-v\\ s\end{array}\right]\text{ (4)}$

Two rotations can be composed by multiplying their quaternion representations. A quaternion q3 representing a rotation q1 followed by a rotation q2 is simply q3 q2q1. The rotation of a three- dimensional vector x by a unit quaternion q is

$\widehat{x}\mathrm{<mo>\text{'}</mo><mo>=</mo>}\widehat{q}\times q^{†}\text{ (5)}$

where $x^{\wedge }$ indicates the formation of a quaternion with zero scalar part from the vector x.

$\widehat{x}\mathrm{<mo>=</mo>}\left[\begin{array}{c}x\\ 0\end{array}\right]\text{ (6)}$

Analytic functions can be defined for quaternion arguments in much the same way as for complex numbers and matrices. In particular, the quaternion exponential has a simple closed-form expression in terms of the quaternion’s scalar and vector parts:

$e^q\mathrm{<mo>=</mo>}{\displaystyle \sum _{n\mathrm{<mo>=</mo>0}}^{\infty }}\frac{q^n}{n\mathrm{<mo>!</mo>}}\mathrm{<mo>=</mo>}{e}^s\left[\begin{array}{c}\frac{v}{\mathrm{<mo>|</mo>}v\mathrm{<mo>|</mo>}}sin\mathrm{<mo\; stretchy="false">(</mo><mo>|</mo>}v\mathrm{<mo>|</mo><mo\; stretchy="false">)</mo>}\\ cos\mathrm{<mo\; stretchy="false">(</mo><mo>|</mo>}v\mathrm{<mo>|</mo><mo\; stretchy="false">)</mo>}\end{array}\right]\text{ (7)}$

The formula for a rotation quaternion in Eq. (3) can be compactly written in terms of the exponential:

${\widehat{e}}^{r\theta \mathrm{<mo>/</mo>2}}\mathrm{<mo>=</mo>}\left[\begin{array}{c}r\sin \left(\frac{\theta }{2}\right)\\ \cos \left(\frac{\theta }{2}\right)\end{array}\right]\text{ (8)}$

Finally, in addition to the purely algebraic properties of quaternions outlined so far, the subsequent analysis requires some kinematic identities relating quaternion derivatives to vector quantities more familiar in rigid-body dynamics. First, the time derivative of a body’s attitude quaternion is related to its angular velocity in the following way:

$\widehat{\omega }\mathrm{<mo>=</mo>2}{q}^{†}q\text{ (9)}$

Second, the quaternion generalized force corresponding to a torque on the body is [13,14].

$F\mathrm{<mo>=</mo>2}q\widehat{\pi }\text{ (10)}$

Schaub and Junkins provide a thorough discussion of rigid-body dynamics using quaternions [15].

3. Euler’s equation from Hamilton’s principle

This section presents a detailed derivation of the classical Euler equation using Hamilton’s principle. Although the results are not new, the techniques used provide the foundation for the development of the variational integrators in the following sections. A more in- depth treatment of variational mechanics on Lie groups, including the rotation group SO (3), is given by Holm [16].

The derivation begins with the Lagrangian for a free rigid body, which, in the absence of a potential, is simply its kinetic energy:

$L\mathrm{<mo>=</mo>}\frac{1}{2}\omega \cdot I\cdot \omega \mathrm{<mo>=</mo>}\frac{1}{2}\widehat{\omega }\cdot J\cdot \widehat{\omega }\text{ (11)}$

J is the following augmented inertia matrix:

$J\mathrm{<mo>=</mo>}\left[\begin{array}{cccc}{I}_{11}& I_{12}& I_{13}& 0\\ I_{21}& I_{22}& I_{23}& 0\\ I_{31}& I_{32}& I_{33}& 0\\ 0& 0& 0& 0\end{array}\right]\text{ (12)}$

Following the standard approach in variational mechanics [17,18], an action integral is constructed and its variational derivative is taken:

$\delta S\mathrm{<mo>=</mo>}\delta {\displaystyle {\int }_{t0}^t}f\cdot \frac{1}{2}\widehat{\omega }\mathrm{.}J\mathrm{.}\widehat{\omega }\mathrm{.}dt\mathrm{<mo>=</mo>0}\text{ (13)}$

At this point, care must be taken to vary  in such a way that the quaternion unit-norm constraint is maintained. There are several ways to explicitly enforce the constraint in the action integral [14,19]. An alternative is to incorporate the constraint directly into the variation [4,6,16]. From the fact that the exponential of a quaternion with zero scalar part is always a unit quaternion, a varied unit quaternion can be defined as

${ϵ}_q\mathrm{<mo>=</mo>}q{e}^{ϵ\widehat{\eta }}\text{ (14)}$

Where a left superscript  is used to denote a varied quantity. Next, q is differentiated with respect to time:

${ϵ}_q\mathrm{<mo>=</mo>}q{e}^{ϵ\widehat{\eta }}+ϵ\left(q{e}^{ϵ\widehat{\eta }}\widehat{\eta }\right)\text{ (15)}$

Equations (14) and (15) can be substituted into the identity in Eq. (9) to obtain the desired variation of , keeping in mind that only terms linear in  need to be retained:

$ϵ\widehat{\omega }{\mathrm{<mo>=</mo>2}}^{ϵ}{q}^{†ϵ}q\mathrm{<mo>=</mo>}{e}^{-ϵ\widehat{\eta }}\widehat{\omega }{e}^{ϵ\widehat{\eta }}\approx \widehat{\omega }+ϵ\left(\widehat{\omega }\widehat{\eta }-\widehat{\eta }\widehat{\omega }+2\widehat{\eta }\right)\text{ (16)}$

Using Eq. (16), the variational derivative of the action integral in Eq. (13) is

$\delta S\mathrm{<mo>=</mo>}\frac{d}{dϵ}\mathrm{<mo>|</mo>}ϵ\mathrm{<mo>=</mo>0}{\displaystyle {\int }_{t_0}^{t_f}}\widehat{\omega }\mathrm{.}J\mathrm{.<mo\; stretchy="false">(</mo>2}\widehat{\omega }\mathrm{.}J\mathrm{.<mo\; stretchy="false">(</mo>2}\widehat{\omega }\widehat{\eta }+2\widehat{\eta }\mathrm{<mo\; stretchy="false">)</mo>}dt\mathrm{<mo>=</mo>0}\text{ (17)}$

Following the usual procedure, integration by parts is used to eliminate , noting that variations must be zero at the endpoints of the integration interval:

$\delta S\mathrm{<mo>=</mo>}{\displaystyle {\int }_{t_0}^{t_f}}2\widehat{\omega }\mathrm{.}J\mathrm{.}\widehat{\omega }\widehat{\eta }-2\widehat{\omega }\mathrm{.}J\mathrm{.}\widehat{\eta }dt\mathrm{<mo>=</mo>0}\text{ (18)}$

Because all of the quaternions in Eq. (18) have scalar parts equal to zero, it can be converted to vector form:

$\delta s\mathrm{<mo>=</mo>}{\displaystyle {\int }_{t_0}^{t_f}}\widehat{\omega }\mathrm{.}I\mathrm{.<mo\; stretchy="false">(</mo>}\omega \times n\mathrm{<mo\; stretchy="false">)</mo>}{-}_{\omega }\mathrm{.}I\mathrm{.}\eta dt\mathrm{<mo>=</mo>0}\text{ (19)}$

Using the fact that cyclically permuting the factors in a scalar triple product does not change its value, Eq. (19) can be rewritten as

$\delta s\mathrm{<mo>=</mo>}{\displaystyle {\int }_{t_0}^{t_f}}\eta \cdot \left[\mathrm{<mo\; stretchy="false">(</mo>}I\cdot \omega \mathrm{<mo\; stretchy="false">)</mo>}\times \omega \right]-\omega \cdot I\cdot \eta dt\mathrm{<mo>=</mo>0}\text{ (20)}$

Finally, recognizing that equality must hold for all perturbations results in Euler's equation:

$I\mathrm{.}\omega +\omega \times I\mathrm{.}\omega \mathrm{<mo>=</mo>0}\text{ (21)}$

4. Discrete-time Euler's equation

This section derives an algebraic equation that is a discrete-time analogue of Euler’s equation. The ideas used, collectively known as discrete mechanics, are presented in detail by Marsden and West [1]. The derivation here roughly follows that of Lee et al. [4] but uses quaternions where they have used rotation matrices.

The point of departure from classical mechanics is the action integral in Eq. (13). It is first broken into finite short segments of length h, with $t_k\mathrm{<mo>=</mo>}{t}_0+kh$ :

$S\mathrm{<mo>=</mo>}{\displaystyle {\int }_{t_0}^{t_f}}\frac{1}{2}\widehat{\omega }\mathrm{.}J\mathrm{.}\widehat{\omega }dt\mathrm{<mo>=</mo>}{\displaystyle \sum _{K\mathrm{<mo>=</mo>0}}^{N-1}}{\displaystyle {\int }_{t_k}^{t_k+1}}\frac{1}{2}\widehat{\omega }\mathrm{.}J\mathrm{.}\widehat{\omega }dt\text{ (22)}$

The integral of the Lagrangian over a single time step h inside the summation on the right-hand side of Eq. (22) is known as the exact discrete Lagrangian [1]:

$L_d^E\mathrm{<mo>=</mo>}{\displaystyle {\int }_{t_k}^{t_k+1}}\frac{1}{2}\widehat{\omega }\mathrm{.}J\mathrm{.}\widehat{\omega }dt\text{ (23)}$

The next step is to approximate $L_d^E$ using a quadrature rule. Any quadrature rule for approximating integrals can be used for this purpose, with higher order rules generally leading to higher order variational integrators [1]. The resulting approximation is known as the discrete Lagrangian of the system. In general, different quadrature rules lead to different discrete Lagrangians. For simplicity and clarity, the rectangle rule is used here. First, a finite difference approximation of  is defined:

${\widehat{\omega }}_k\mathrm{<mo>=</mo>2}{q}_k^{†}{q}_k\approx 2q_k^{†}\mathrm{<mo\; stretchy="false">(</mo>}\frac{q_k+\mathrm{1<mo>=</mo>}{q}_k}{h}\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>2<mo\; stretchy="false">(</mo>}\frac{f_k-1}{h}\mathrm{<mo\; stretchy="false">)</mo>}\text{ (24)}$

The quaternion rotation from q_k to q_k+1 is denoted by $f_k\mathrm{<mo>=</mo>}{q}_k^{†}{q}_{k+1}$ . Substituting the approximation for ${\widehat{\omega }}_k$ into Eq. (23), applying the rectangle rule, and simplifying leads to the following discrete Lagrangian:

$L_d\mathrm{<mo>=</mo>}\frac{2}{h}{f}_k\mathrm{.}J\mathrm{.}{f}_k\text{ (25)}$

Using Eq. (25), the discrete action sum for the system can be formed:

$S_d\mathrm{<mo>=</mo>}{\displaystyle \sum _{K\mathrm{<mo>=</mo>0}}^{N-1}}{L}_d\mathrm{<mo>=</mo>}{\displaystyle \sum _{K\mathrm{<mo>=</mo>0}}^{N-1}}\frac{2}{h}{f}_k\mathrm{.}J\mathrm{.}{f}_k\text{ (26)}$

Equation (26) approximates Eq. (22) and serves the same role in discrete mechanics as the action integral does in traditional variational mechanics [1]. Analogously to the continuous case, Hamilton's principle is applied to the action sum. First, a varied f_k that obeys the unit quaternion constraint is needed:

${ϵ}_{f_k}\mathrm{<mo>=</mo>}{ϵ}_{q_k}^{†}{q}_{k+1}\mathrm{<mo>=</mo>}{e}^{-\widehat{ϵ}\eta }{f}_k^{\widehat{ϵ}\eta }\approx f_k+ϵ\mathrm{<mo\; stretchy="false">(</mo>}{f}_k{\widehat{\eta }}_{k+1}-{\widehat{\eta }}_k{f}_k\mathrm{<mo\; stretchy="false">)</mo>}\text{ (27)}$

Using ${ϵ}_{f_k}$ , the variation of the action sum is set equal to zero:

$\delta S_d\mathrm{<mo>=</mo>}\frac{d}{dϵ}{\mathrm{<mo>|</mo>}}_{ϵ\mathrm{<mo>=</mo>0}}{\displaystyle \sum _{K\mathrm{<mo>=</mo>0}}^{N-1}}\frac{2}{h}ϵfk\mathrm{.}J\mathrm{.}{ϵ}_{fk}\mathrm{<mo>=</mo>}{\displaystyle \sum _{K\mathrm{<mo>=</mo>0}}^{N-1}}\frac{4}{h}{f}_k\mathrm{.}J\mathrm{.<mo\; stretchy="false">(</mo>}{f}_k{\widehat{\eta }}_{k+1}-{\widehat{\eta }}_k{f}_k\mathrm{<mo>=</mo>0}\text{ (28)}$

The next step is to eliminate ${\eta }_{k+1}$ from the right-hand side of Eq. (28) by performing the discrete equivalent of integration by parts, which amounts to some simple index manipulation:

$\delta S_d\mathrm{<mo>=</mo>}{f}_{N-1}\mathrm{.}J\mathrm{.}{f}_{N-1}{\widehat{\eta }}_N-f_0\mathrm{.}J\mathrm{.}{f}_0{\widehat{\eta }}_0\mathrm{<mo>=</mo>}{\displaystyle \sum _{K\mathrm{<mo>=</mo>1}}^{N-1}}\frac{4}{h}{f}_{k-1}\mathrm{.}J\mathrm{.<mo\; stretchy="false">(</mo>}{f}_{k-1}{\widehat{\eta }}_k-{\widehat{\eta }}_k{f}_k\mathrm{<mo>=</mo>0}\text{ (29)}$

Using the fact that variations at the endpoints must be zero, just as in the continuous case, the first two terms in Eq. (29) can be eliminated:

$\delta S_d\mathrm{<mo>=</mo>}{\displaystyle \sum _{k\mathrm{<mo>=</mo>1}}^{N-1}}\frac{4}{h}{f}_{k-1}\cdot J\cdot \left(f_{k-1}{\widehat{\eta }}_k-{\widehat{\eta }}_k{f}_k\right)\mathrm{<mo>=</mo>0}\text{ (30)}$

At this point, Eq. (30), which implicitly includes unit-norm constraints on the quaternions, is converted to an unconstrained vector equation by parameterizing f_k in the following way:

$f_k\mathrm{<mo>=</mo>}\frac{{\varphi }_k}{\sqrt{1-{\varphi }_k\mathrm{.}{\varphi }_k}}\text{ (31)}$

This parameterization is only valid for $\parallel {\varphi }_k\parallel \mathrm{<mo><</mo>1}$ . Therefore, h must be chosen small enough to ensure that the incremental rotations between adjacent time steps are less than 1800. A number of other three-parameter attitude representations could be used instead (modified Rodrigues parameters, for example); however, Eq. (31) is a natural choice that leads to simple and elegant expressions. In terms of φ_k, Eq. (30) is

${\displaystyle \sum _{K\mathrm{<mo>=</mo>1}}^{N-1}}\left[\begin{array}{c}\sqrt{1-{\varphi }_k}\mathrm{.}\varphi n_k\mathrm{.}I\mathrm{.}{\varphi }_k-{\varphi }_k\mathrm{.}I\mathrm{.<mo\; stretchy="false">(</mo>}\varphi \times {\eta }_k\mathrm{<mo\; stretchy="false">)</mo>}\\ -\sqrt{1-{\varphi }_{k-1}}\mathrm{.}{\varphi }_{k-1}{\eta }_k\mathrm{.}I\mathrm{.}{\varphi }_{k-1}-{\varphi }_{k-1}\mathrm{.}I\mathrm{.<mo\; stretchy="false">(</mo>}{\varphi }_{k-1}\times {\eta }_k\mathrm{<mo\; stretchy="false">)</mo>}\end{array}\right]\mathrm{<mo>=</mo>0}\text{ (32)}$

Recognizing that Eq. (32) must be true for all ${\eta }_k$ and performing some simple vector algebra reveals an algebraic equation relating the incremental rotation from the last time step to the current time step φ_k to the incremental rotation from the current time step to the next time step φ_k+1:

$\sqrt{1-{\varphi }_k}\mathrm{.}{\varphi }_{k}I\mathrm{.}{\varphi }_k-{\varphi }_k\times I\mathrm{.}{\varphi }_k\mathrm{<mo>=</mo>}\sqrt{1-\varphi -k+1}\mathrm{.}{\varphi }_{k+1}I\mathrm{.}{\varphi }_{k+1}\times I\mathrm{.}\varphi -k+1\text{ (33)}$

As a brief aside, Eq. (33) bears some resemblance to the classical Euler’s equation. Taking its limit as h goes to zero does, in fact, recover Eq. (21). This result confirms that the discrete-time equation converges to the true differential equation for small time steps and establishes consistency with the continuous theory.

5. Variational integrator for the free rigid body

This section uses Eq. (33) as the starting point for the development of a variational integrator for the unforced rigid body. The additional ingredients needed are a way to initialize the integrator given an attitude q₀ and angular velocity ω₀, a way to update the attitude q_k+1 and angular velocity ω_k+1 after solving for φ_k+1, and a way to solve Eq. (33) for φ_k+1 given φ_k. Although it might seem simple enough to approximate φ₀ in any number of ways given ω₀, ad hoc approaches do not maintain the variational integrator's conservation properties. The discrete Legendre transform [1] gives a consistent way to convert between φ_k and ω_k. Similar to the classical Legendre transform [17], it maps from φ_k (which is effectively the discrete-time velocity variable) to p_k, the momentum at time k, which can then be multiplied by I⁻¹ to recover ω_k. Unlike the continuous version, there are actually two discrete Legendre transforms for a given time step [1]:

$p_{{-}_k}\mathrm{<mo>=</mo>}-\frac{\partial L_d\mathrm{<mo\; stretchy="false">(</mo>}{q}_k\mathrm{,}{q}_{k+\mathrm{1<mo\; stretchy="false">)</mo>}}}{\partial q_k}\mathrm{.}\partial q_k\text{ (34)}$

$p_k^{+}\mathrm{<mo>=</mo>}\frac{\partial L_d\mathrm{<mo\; stretchy="false">(</mo>}{q}_k\mathrm{,}{q}_{k+1}\mathrm{<mo\; stretchy="false">)</mo>}}{\partial q_{k+1}}\text{d}{q}_{k+1}\text{ (35)}$

Applying these transformations to the discrete Lagrangian in Eq. (25) reveals that $p_k^{-}$ and p_k correspond to the left and right sides of Eq. (33):

$p_{-k}\mathrm{<mo>=</mo>}\frac{2}{h}\sqrt{1-{\varphi }_k}\mathrm{.}{\varphi }_{k}I\mathrm{.}{\varphi }_k-{\varphi }_k\times I\mathrm{.}{\varphi }_k\text{ (36)}$

$p_{+k}\mathrm{<mo>=</mo>}\frac{2}{h}\sqrt{1-{\varphi }_{k+1}}\mathrm{.}{\varphi }_{k+1}I\mathrm{.}{\varphi }_{k+1}-{\varphi }_{k+1}\times I\mathrm{.}{\varphi }_{k+1}\text{ (37)}$

This result leads to several key conclusions. First, Eqns. (36) and (37) provide a new interpretation of the discrete-time equation of motion as a momentum balance between adjacent time steps. Second, Eq. (36) can be used to initialize the integrator by solving for φ₀ given I and ω₀. Lastly, p_k, and hence ω_k, can be calculated at any point during the integration using either Eqns. (36) or (37).

The final missing piece of the integrator is a method for solving Eq. (33), which is both implicit and nonlinear. Newton’s method, which amounts to solving successive linear approximations of the equation until a desired level of accuracy is achieved [3,20], provides an efficient solution in this case. The necessary linear approximation is the Jacobian matrix of Eq. (37),

$\frac{\partial p_k}{\partial p_{k+1}}\mathrm{<mo>=</mo>}\frac{2}{h}\left[\begin{array}{c}\sqrt{1-{\varphi }_{k+1}^{\top }}\mathrm{.}{\varphi }_{k+1}I\sqrt{I{\varphi }_{k+1}+{\varphi }_{k+1}^{\top }\times I\mathrm{.}\varphi -k+1}\\ +skew\mathrm{<mo\; stretchy="false">(</mo>}{\varphi }_{k+1}\mathrm{<mo\; stretchy="false">)</mo>}I-skew\mathrm{<mo\; stretchy="false">(</mo>}{\varphi }_{k+1}\mathrm{<mo\; stretchy="false">)</mo>}\end{array}\right]\text{ (38)}$

where skew φ indicates the skew-symmetric matrix-multiplication equivalent of the cross-product operation:

$skew\text{}\left[\begin{array}{c}{\varphi }_1\\ {\varphi }_2\\ {\varphi }_3\end{array}\right]\mathrm{<mo>=</mo>}\left[\begin{array}{ccc}0& -{\varphi }_3& {\varphi }_2\\ {\varphi }_3& 0& -{\varphi }_1\\ -{\varphi }_2& {\varphi }_1& 0\end{array}\right]\text{ (39)}$

Three or four Newton iterations are sufficient to reach machine precision in all of the examples presented in Section 9 using standard 64-bit floating-point arithmetic. Once φ_k+1 is computed, the attitude can be updated by simple quaternion multiplication with the previous attitude, $q_{k+1}\mathrm{<mo>=</mo>}{q}_k{q}_{{\varphi }_k}$ .

In summary, given an inertia and initial conditions q₀ and ω₀, the integrator is initialized by computing the momentum $p_0\mathrm{<mo>=</mo>}I\cdot {\omega }_0$ and an initial guess for φ₀. A reasonable guess is

${\varphi }_0\approx \frac{h}{2}{\omega }_0.$

The true value of φ₀ is then calculated to machine precision using Newton’s method with Eqns. (37) and (38). From φ₀, p1 is calculated using Eq. (36), followed by ω1 and q₁. The process is then repeated as desired.

6. Gyrostats

A gyrostat is a system of coupled rigid bodies whose relative motions do not change the total inertia tensor of the system [21]. A practical example is a rigid body with internal rotors or momentum actuators that can spin relative to the carrier body, such as a spacecraft with reaction wheels. Using the variational framework developed in the preceding sections, both the classical equations of motion and a variational integrator are straightforward to derive. The Lagrangian for a gyrostat system is

$L\mathrm{<mo>=</mo>}\frac{1}{2}{\omega }_B\mathrm{.}{I}_B\mathrm{.}{\omega }_B+{\displaystyle \sum _{r\mathrm{<mo>=</mo>0}}^N}\frac{1}{2}\mathrm{<mo\; stretchy="false">(</mo>}{\omega }_B+{\omega }_r\mathrm{<mo\; stretchy="false">)</mo>.}{I}_r\mathrm{.<mo\; stretchy="false">(</mo>}{\omega }_B+{\omega }_r\mathrm{<mo\; stretchy="false">)</mo>}+\frac{1}{2}{m}_r{\mathrm{<mo\; stretchy="false">(</mo>}{\omega }_B\times x_r\mathrm{<mo\; stretchy="false">)</mo>}}^2\text{ (40)}$

where I_B and ω_B are the carrier body’s inertia tensor and angular velocity, the Ir are the rotor inertia tensors, the ωr are the rotor angular velocities relative to the carrier body, and the Xr are the rotor positions relative to the carrier body’s center of mass. The Lagrangian can be simplified by introducing a modified body inertia.

$I_B^0$ , which includes the rotor masses,

$I_B^1\mathrm{<mo>=</mo>}{I}_B-{\displaystyle \sum _{r\mathrm{<mo>=</mo>1}}^N}{m}_{r}skew{\mathrm{<mo\; stretchy="false">(</mo>}{x}_r\mathrm{<mo\; stretchy="false">)</mo>}}^2\text{ (41)}$

Substituting $I_B^0$ into Eq. (40) eliminates the last term, giving the simpler expression

$L\mathrm{<mo>=</mo>}\frac{1}{2}{\omega }_B\mathrm{.}{I}_B^{\mathrm{<mo>\text{'}</mo>}}\mathrm{.}{\omega }_B+{\displaystyle \sum _{r\mathrm{<mo>=</mo>1}}^N}\frac{1}{2}\mathrm{<mo\; stretchy="false">(</mo>}{\omega }_B+{\omega }_r\mathrm{<mo\; stretchy="false">)</mo>.}{I}_r\mathrm{.<mo\; stretchy="false">(</mo>}{\omega }_B+{\omega }_r\mathrm{<mo\; stretchy="false">)</mo>}\text{ (42)}$

The rotor angular velocities ωr are treated as exogenous inputs to the system that can be set arbitrarily (e.g., by a controller), and so variations need only be taken with respect to ω_B. This fact makes the derivation for the gyrostat nearly identical to the free rigid body. The only difference is that a few extra terms involving Ir and ω_r are carried through. Using ω_B to vary the action and following the rest of the steps in Section 3 results in the following differential equation:

$I_B^{\mathrm{<mo>\text{'}</mo>}}\mathrm{.}{\omega }_B+{\omega }_B\times I_B^{\mathrm{<mo>\text{'}</mo>}}\mathrm{.}{\omega }_B+{\displaystyle \sum _{r\mathrm{<mo>=</mo>1}}^N}{I}_B\mathrm{.}{\omega }_B+{\omega }_b\times I_r\mathrm{.}{\omega }_B+I_r\mathrm{.}{\omega }_r+{\omega }_b\times I_r\mathrm{.}{\omega }_r\mathrm{<mo>=</mo>0}\text{ (43)}$

Two new definitions help simplify Eq. (43). First, the gyrostat inertia I_G is

$I_G\mathrm{<mo>=</mo>}{I}_B^{\mathrm{<mo>\text{'}</mo>}}+{\displaystyle \sum _{r\mathrm{<mo>=</mo>1}}^N}{I}_r\mathrm{<mo>=</mo>}{I}_B+{\displaystyle \sum _0^N}{I}_r-m_{r}skew{\mathrm{<mo\; stretchy="false">(</mo>}{x}_r\mathrm{<mo\; stretchy="false">)</mo>}}^2\text{ (44)}$

Second, ρ is the total angular momentum stored in all the rotors:

$\rho \mathrm{<mo>=</mo>}{\displaystyle \sum _{r\mathrm{<mo>=</mo>1}}^N}{I}_r\mathrm{.}{\omega }_r\text{ (45)}$

Substituting I_G and ρ into Eq. (43) results in the classical equation of motion for the gyrostat [21]:

$I_G\mathrm{.}{\omega }_B+{\omega }_B\times \mathrm{<mo\; stretchy="false">(</mo>}{I}_G\mathrm{.}{\omega }_B+\rho \mathrm{<mo\; stretchy="false">)</mo>}+\rho \mathrm{<mo>=</mo>0}\text{ (46)}$

The steps involved in deriving the discrete-time equivalent of Eq. (46) are now briefly highlighted. If ω_B is approximated as in Eq. (24) and the rectangle rule is used, the discrete Lagrangian for the gyrostat is

$L_d\mathrm{<mo>=</mo>}\frac{2}{h}\left[f_k\cdot J_B^{\text{'}}\cdot f_k+{\displaystyle \sum _{r\mathrm{<mo>=</mo>1}}^N}\left(f_k+\frac{h}{2}{\widehat{\omega }}_{r\mathrm{,}k}\right)J_r\left(f_k+\frac{h}{2}{\widehat{\omega }}_{r\mathrm{,}k}\right)\right]\text{ (47)}$

where $J_B^0$ and J_r are the augmented 4 × 4 equivalents formed as in Eq. (26) and its variation taken using ^εf_k of $I_B^0$ and Ir. The discrete action sum Sd can then be formed as in Eq. (26) and its variation taken using ^εf_k from Eq. (27):

$\partial S_d\mathrm{<mo>=</mo>}{\displaystyle \sum _{k\mathrm{<mo>=</mo>1}}^{N-1}}\left[f_k\cdot J_B\mathrm{.}\left(f_k{\widehat{\eta }}_{k+1}-{\widehat{\eta }}_k{f}_k\right)+{\displaystyle \sum _{r\mathrm{<mo>=</mo>1}}^N}\left(f_k+\frac{h}{2}{\widehat{\omega }}_{r\mathrm{,}k}\right)\mathrm{.}{J}_r\mathrm{.}\left(f_k{\widehat{\eta }}_{k+1}-{\widehat{\eta }}_k{f}_k\right)\mathrm{<mo>=</mo>0}\right]\text{ (48)}$

Following the rest of the steps in Section 4 and substituting in I_G and ρ results in the discrete-time gyrostat equation:

$\begin{array}{l}\sqrt{1-{\varphi }_k\mathrm{.}{\varphi }_k}\left[I_G\mathrm{.}{\Phi }_k+\frac{h}{2}{\rho }_k\right]-{\varphi }_k\times \left[I_G\mathrm{.}{\varphi }_k+\frac{h}{2}{\rho }_k\right]\mathrm{<mo>=</mo>}\\ \sqrt{1-{\varphi }_{k+1}\mathrm{.}{\varphi }_{k+1}}\left[I_G\mathrm{.}{\Phi }_{k+1}+\frac{h}{2}{\rho }_{k+1}\right]-{\varphi }_{k+1}\times \left[I_G\mathrm{.}{\varphi }_{k+1}+\frac{h}{2}{\rho }_{k+1}\right]\text{ (49)}\end{array}$

Equation (49) can be directly substituted into the variational integrator developed in Section 5. The only other change necessary is to the Jacobian in the Newton iteration:

$\frac{\partial p^{+}}{\partial {\varphi }_{k+1}}\mathrm{<mo>=</mo>}\frac{2}{h}\left[\sqrt{1-{\varphi }^{\top }\varphi I}-\frac{I\varphi {\varphi }^{\top }}{\sqrt{1-}{\varphi }^{\top }\varphi }+skew\mathrm{<mo\; stretchy="false">(</mo>}\varphi \mathrm{<mo\; stretchy="false">)</mo>}I-skew\mathrm{<mo\; stretchy="false">(</mo>}I\varphi \mathrm{<mo\; stretchy="false">)</mo>}-\frac{h}{2}\frac{\rho {\varphi }^{\top }}{\sqrt{1-{\varphi }^{\top }\varphi }}-\frac{h}{2}skew\mathrm{<mo\; stretchy="false">(</mo>}\rho \mathrm{<mo\; stretchy="false">)</mo>}\right]\text{ (50)}$

7. External torques

External torques can be incorporated into the variational framework using the Lagrange–d’Alembert principle (often known simply as d’Alembert’s principle) [17,18]. In particular, its integral form [1,22]:

$\partial {\displaystyle {\int }_{t0}^{t_f}}Ldt+{\displaystyle {\int }_{t0}^{t_f}}F\mathrm{.}\partial qdt\mathrm{<mo>=</mo>0}\text{ (51)}$

is most readily applied here, where the term on the left is simply the variation of the action δS and the term on the right is the integral of the virtual work done by a generalized force F. To apply the Lagrange–d’Alembert principle to a rigid body, the expression for the quaternion generalized force given in Eq. (10), as well as the variational derivative of the attitude quaternion δq, are substituted into the second term of Eq. (51):

${\displaystyle {\int }_{t_0}^{t_f}}F\delta qdt\mathrm{<mo>=</mo>}{\displaystyle {\int }_{t_0}^{t_f}}2q\widehat{\tau }\cdot \mathrm{<mo\; stretchy="false">(</mo>}q\widehat{\eta }\mathrm{<mo\; stretchy="false">)</mo>}dt\text{ (52)}$

A little algebra reveals that $q\widehat{\tau }\cdot q\widehat{\eta }\tau \cdot \eta \mathrm{,}$ further simplifying the expression. Combining this result with the action term from Eq. (18) leads to the following equation:

${\displaystyle {\int }_{t_0}^{t_f}}\left(2\widehat{\omega }\cdot J\cdot \widehat{\omega }\widehat{\eta }-2\widehat{\omega }\cdot J\cdot \widehat{\eta }+2\widehat{\tau }\cdot \widehat{\eta }\right)dt\mathrm{<mo>=</mo>0}F\mathrm{.}\delta qdt\mathrm{<mo>=</mo>}\text{ (53)}$

It is then straightforward to work through the rest of the steps in Section 3 to arrive at the forced Euler equation:

$I\mathrm{.}\omega +\omega \times I\mathrm{.}\omega \mathrm{<mo>=</mo>}\tau \text{ (54)}$

Incorporating forcing into the discrete variational framework is a bit more subtle than in the continuous case. The discrete version of the Lagrange–D’Alembert principle is

$\delta {\displaystyle \sum _{k\mathrm{<mo>=</mo>1}}^N}{L}_d+\delta {\displaystyle \sum _{k\mathrm{<mo>=</mo>1}}^N}{F}_d^{-}\mathrm{.}\delta q_k+F_d^{+}\mathrm{.}\delta q_{k+1}\mathrm{<mo>=</mo>0}\text{ (55)}$

where $F_d^{-}$ and F_d are known as discrete generalized forces [1]. Similar to what is encountered with the discrete Legendre transform, there are two discrete generalized forces corresponding to the beginning and end of a time step:

$F_d^{-}\mathrm{<mo>=</mo>}{\displaystyle \sum _{t_k}^{t_{k+1}}}\frac{1}{2}F\mathrm{<mo\; stretchy="false">(</mo>}q\mathrm{,}q\mathrm{<mo\; stretchy="false">)</mo>.}\frac{\partial q\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}}{\partial q_k}dt\text{ (56)}$

$F_d^{+}\mathrm{<mo>=</mo>}{\displaystyle \sum _{t_k}^{t_{k+1}}}\frac{1}{2}F\mathrm{<mo\; stretchy="false">(</mo>}q\mathrm{,}q\mathrm{<mo\; stretchy="false">)</mo>.}\frac{\partial q\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}}{\partial q_{k+1}}dt\text{ (57)}$

As with the discrete Lagrangian, the integrals in the definition of the discrete generalized forces are approximated by quadrature. Here, the rectangle rule is used, though more accurate quadrature rules can lead to more accurate integrators at the expense of increased computational burden:

$F_d^{-}\approx h{q}_k{\widehat{\tau }}_k\text{ (58)}$

$F_d^{+}\approx h{q}_{k+1}{\widehat{\tau }}_{k+1}\text{ (59)}$

Substituting the approximations for $F_d^{-}$ and F_d , as well as the discrete Lagrangian for the rigid body from Eq. (25), into Eq. (55) results in

$\delta {\displaystyle \sum _{k\mathrm{<mo>=</mo>0}}^N}\frac{2}{h}{f}_k\mathrm{.}J\mathrm{.}{f}_k+{\displaystyle \sum _{k\mathrm{<mo>=</mo>0}}^N}h{q}_k\widehat{\tau }k\mathrm{.}\delta q_{k+1}{\widehat{\tau }}_{k+1}\mathrm{.}\delta q_{k+1}\mathrm{<mo>=</mo>0}\text{ (60)}$

Carrying out the variational derivatives leads to the following:

$\delta {\displaystyle \sum _{k\mathrm{<mo>=</mo>0}}^N}\frac{4}{h}{f}_k\mathrm{.}J\mathrm{<mo\; stretchy="false">(</mo>}{f}_k{\eta }_{k+1}-{\eta }_k{f}_k\mathrm{<mo\; stretchy="false">)</mo>}+h{\widehat{\tau }}_k\mathrm{.}{\widehat{\eta }}_k+h{\widehat{\tau }}_{k+1}\mathrm{.}{\widehat{\eta }}_{k+1}\mathrm{<mo>=</mo>0}\text{ (61)}$

Eliminating the ${\eta }_{k+1}$ terms again requires a “discrete integration by parts.” Manipulating indices results in

$\delta {\displaystyle \sum _{k\mathrm{<mo>=</mo>0}}^N}\frac{4}{h}{f}_{k-1}\mathrm{.}J\mathrm{<mo\; stretchy="false">(</mo>}{f}_{k-1}{\eta }_k-{\eta }_k{f}_K\mathrm{<mo\; stretchy="false">)</mo>}+2h{\widehat{\tau }}_k\mathrm{.}{\widehat{\eta }}_k\mathrm{<mo>=</mo>0}\text{ (62)}$

Retracing the remaining steps in Section 4 yields the discrete-time equation of motion for the forced rigid body:

$\sqrt{1-{\varphi }_k\cdot I\cdot {\varphi }_k}{\varphi }_k-\left({\varphi }_k\times \mathrm{<mo\; stretchy="false">(</mo>}I\cdot {\varphi }_k\mathrm{<mo\; stretchy="false">)</mo>}\right)\mathrm{<mo>=</mo>}\sqrt{1-{\varphi }_{k+1}\cdot I\cdot {\varphi }_{k+1}}{\varphi }_{k+1}-\left({\varphi }_{k+1}\times \mathrm{<mo\; stretchy="false">(</mo>}I\cdot {\varphi }_{k+1}\mathrm{<mo\; stretchy="false">)</mo>}\right)+\frac{h^2}{2}{\tau }_{k+1}\text{ (63)}$

This result can also be readily applied to the forced gyrostat by adding the same torque term onto the right-hand side of Eq. (49).

8. Gyrostat spacecraft with damping

The tools developed up to this point enable the construction of a variational integrator for a gyrostat with an internal energy-dissipating mechanism. The mechanism considered here is known as a Kane damper and consists of a spherical mass immersed in a viscous fluid inside a spherical cavity in the spacecraft body [23]. The torque exerted on the spacecraft by the damper is

$\tau \mathrm{<mo>=</mo>}C\mathrm{<mo\; stretchy="false">(</mo>}{\omega }_D-{\omega }_B\mathrm{<mo\; stretchy="false">)</mo>}\text{ (64)}$

where C is a damping constant, ω_D is the damper angular velocity, and ω_B is the body angular velocity. The basic approach taken here is to treat the body and damper as separate rigid bodies coupled through the viscous damping force. Equation (64) is approximated by finite differences in the usual way, giving

${\tau }_k\approx \frac{2}{h}C\mathrm{<mo\; stretchy="false">(</mo>}{\gamma }_k-{\varphi }_k\mathrm{<mo\; stretchy="false">)</mo>}\text{ (65)}$

where ${\varphi }_k$ is the incremental body rotation defined in Eq. (31) and γ_k is the analogous incremental rotation for the damper. Substituting this approximation into Eq. (63) yields a set of coupled equations for the gyrostat–damper system:

$\begin{array}{l}\sqrt{1-{\varphi }_k\mathrm{.}{\varphi }_k}\left[I_G\mathrm{.}{\varphi }_k+\frac{h}{2}{\rho }_k\right]-{\varphi }_k\times \left[I_G\mathrm{.}{\varphi }_k+\frac{h}{2}{\rho }_k\right]\mathrm{<mo>=</mo>}\\ \sqrt{1-{\varphi }_{k+1}\mathrm{.}{\varphi }_{k+1}}\left[I_G\mathrm{.}{\varphi }_{k+1}+\frac{h}{2}{\rho }_{k+1}\right]+{\varphi }_{k+1}\times \left[I_G\mathrm{.}{\varphi }_{k+1}+\frac{h}{2}{\rho }_{k+1}\right]+hC\mathrm{<mo\; stretchy="false">(</mo>}{\gamma }_{k+1}-{\varphi }_{k+1}\mathrm{<mo\; stretchy="false">)</mo>}\text{ (66)}\end{array}$

$\sqrt{1-{\gamma }_k\mathrm{.}{\gamma }_k}{I}_D\mathrm{<mo>=</mo>}\sqrt{1-{\gamma }_{k+1}\mathrm{.}{\gamma }_{k+1}}{I}_D\mathrm{.}{\gamma }_{k+1}-hC\mathrm{<mo\; stretchy="false">(</mo>}{\gamma }_{k+1}-{\varphi }_{k+1}\mathrm{<mo\; stretchy="false">)</mo>}\text{ (67)}$

These equations take advantage of the fact that the damper is spherical, and thus has an inertia tensor that is a scalar multiple of the identity, to eliminate the cross-product term in Eq. (67).

Equations (66) and (67) must be solved simultaneously for φ_k+1 and γk+1. Once again, Newton’s method is used, this time with both equations combined to form a single six-dimensional system. The necessary 6×6 Jacobian matrix is easily derived in terms of Eqs. (38) and (50). In addition to the steps outlined in Section 5, a subtlety arises in the implementation of this integrator in that the components of the damper’s angular-momentum vector must be rotated by f_k at the end of each time step to keep them aligned with the spacecraft body frame.

9. Numerical examples

This section presents some numerical examples to demonstrate the performance of the variational integrators derived in Secs. IV–VIII. Comparisons are made to the second-order fixed-step midpoint rule and MATLAB's ODE45 and ODE15s variable-step Runge–Kutta solvers with default error tolerances [24]. The computational cost of the midpoint rule roughly equals that of the variational integrators.In all simulations, the following inertia matrix is used:

$I\mathrm{<mo>=</mo>}\left[\begin{array}{ccc}1& 0& 0\\ 0& 2& 0\\ 0& 0& 3\end{array}\right]\text{ (68)}$

A. Rigid body

The first test compares the energy and momentum behavior of the integrator in Section 5 with the midpoint rule and ODE45 by simulating a free rigid body with an initial angular velocity ${\omega }_0\mathrm{<mo>=</mo>}{\left[\frac{\pi }{4}\mathrm{,}-\frac{\pi }{5}\mathrm{,}\frac{\pi }{6}\right]}^T\text{rad/s}$ . The time steps for the midpoint rule and the variational integrator are chosen to be h = 0.2 s to make the run time for both roughly equal to that of ODE45. Because this system is conservative, both the inertial angular-momentum vector components and the total energy should remain constant throughout the simulation.

The performance of the proposed variational integrator is first evaluated against conventional nonsymplectic schemes. Figure 1 illustrates the applied input torque components in the body frame. The sinusoidal profiles represent externally applied control or disturbance torques driving the rigid-body dynamics. These torque excitations are critical for benchmarking integration schemes since they generate oscillatory responses that challenge numerical stability. The smooth representation of all three components confirms that the variational integrator correctly captures torque transmission without introducing spurious oscillations or numerical artifacts, which are often observed in nonsymplectic methods at larger time steps.

Figure 1: Input torque components applied in the body frame, providing the excitation for rigid-body dynamics simulation.

Figure 2 shows the normalized internal rotor momentum evolution. Here, the momentum components remain bounded and follow physically consistent oscillations over the entire simulation horizon. This behavior highlights the momentum-preserving property of the variational integrator. In contrast, classical Runge–Kutta solvers gradually accumulate drift due to truncation errors, causing artificial energy injection or dissipation into the system. The fact that the variational results align almost exactly with the reference solution demonstrates that the scheme inherently respects conservation laws of rigid-body motion, even under long-duration simulations.

Figure 2: Normalized internal rotor momentum evolution, demonstrating the momentum-preserving property of the variational integrator.

Figure 3 presents the body angular velocity components across an extended propagation interval. The oscillatory motion reflects the natural rotational dynamics of the system under the applied torques. The variational integrator maintains stable angular velocity trajectories for over one million integration steps, showing no evidence of instability or exponential error growth. The slight long-term drift is attributable solely to round-off error from finite-precision floating-point arithmetic. Importantly, this error remains bounded and grows at a much slower rate compared to nonsymplectic integrators, which typically amplify drift due to loss of symplectic structure. This confirms that the variational scheme is particularly well suited for long-duration attitude propagation, where numerical stability is critical.

Figure 3: Body angular velocity components over extended integration, showing stable trajectories with bounded round-off error growth.

Figure 4 depicts the inertial angular momentum components, which remain nearly constant throughout the simulation. This invariance is a direct consequence of the variational integrator’s geometric structure, which enforces conservation of first integrals such as angular momentum. Even over long integration windows, the deviation between the variational and reference solutions remains negligible. By contrast, conventional Runge–Kutta methods introduce gradual dissipation or amplification of momentum because they do not explicitly preserve invariants of motion. The near-perfect overlap of the momentum curves in Figure 4 demonstrates that the variational framework not only stabilizes the numerical solution but also guarantees physical consistency, making it a robust tool for spacecraft dynamics simulations.

Figure 4: Inertial angular momentum components, illustrating near-constant conservation under the variational integrator.

Together, Figures 1-4 confirm that the variational integrator preserves key invariants of motion and provides long-term numerical stability. While higher-order Runge–Kutta schemes can reduce truncation error, they do so at substantially higher computational cost. The variational method strikes a more effective balance, delivering computational efficiency while rigorously maintaining the system’s physical properties. This makes it particularly suitable for spacecraft attitude propagation and real-time guidance, navigation, and control applications.

B. Damped rigid body

The second numerical experiment incorporates the spherical damper described in Section 8 into the rigid-body simulation. The damper inertia is set to I_D = 0.2 and the damping constant C is varied from 0.1 to 100. The classical midpoint rule is excluded from the comparison, since it rapidly diverges as C increases and requires prohibitively small step sizes to maintain stability.

Figure 5 presents the total system energy over time for the case C = 100. Both ODE45 and the proposed variational integrator capture the dissipative behavior of the damped rigid body, showing a monotonic decay of energy. However, the variational integrator achieves this with a fixed time step of h = 0.3s, whereas ODE45 must continually adapt its step size to ensure stability. Despite this, the results remain nearly indistinguishable, confirming that the variational method can reproduce dissipative dynamics with high fidelity while being computationally more efficient.

Figure 5: Total system energy for the damped rigid body at C=100, comparing dissipative behavior captured by ODE45 and the variational integrator.

Figure 6 shows the evolution of attitude quaternion components expressed in terms of system momentum. The variational integrator maintains consistency with ODE45 across all components, accurately reflecting the coupled dynamics between the body and the internal damper. Importantly, the momentum oscillations remain physically interpretable, and no artificial drift is introduced by the variational scheme, even under strong damping conditions.

Figure 6: Evolution of system momentum expressed through quaternion components, highlighting consistency between ODE45 and the variational integrator.

Figure 7 depicts the rotor momentum components. These plots demonstrate the energy exchange between the damper and the primary body dynamics, where the damper dissipates rotational energy and stabilizes the system response. The smooth decay of momentum confirms that the variational framework is able to correctly capture the dissipative interaction between subsystems, a task that often destabilizes traditional integrators unless extremely small step sizes are employed.

Figure 7: Rotor momentum components showing smooth energy exchange and dissipation by the internal damper.

Figure 8 presents the body angular velocity components. While both ODE45 and the variational scheme produce nearly identical results, ODE45 requires adaptive step refinement to avoid instability, whereas the variational method achieves stability with a fixed step size. This robustness highlights a significant advantage for real-time spacecraft applications, where computational budgets are limited and adaptive step solvers may not be practical.

Figure 8: Body angular velocity components under damping, illustrating stability of the variational integrator with a fixed step size.

Together, Figures 5-8 validate that the variational integrator not only conserves momentum and energy in conservative systems (as shown in Section 9A), but also extends naturally to dissipative dynamics without loss of accuracy or stability. This dual capability makes the method particularly attractive for modeling spacecraft equipped with passive damping mechanisms, where long-duration, stable simulations are essential.

C. Extended kalman filter

The final test case demonstrates the advantages of variational integrators in a real-time estimation application. A spacecraft attitude determination problem is simulated where a multiplicative extended Kalman filter (MEKF) [25,26] is used to estimate the attitude quaternion from noisy measurements of two inertial reference vectors. This situation is typical on CubeSats, for example, where magnetometer and sun vector measurements are commonly used for attitude determination.

A simulated truth model is constructed by integrating the rigid-body equations of motion with initial conditions $q_0{\mathrm{0<mo>;</mo>0<mo>;</mo>0<mo>;</mo>1}}^T$ and ${\omega }_0\mathrm{4<mo>;</mo>}-{\mathrm{5<mo>;</mo>6}}^T\text{deg/s}$ using ODE45 in MATLAB. Simulated vector measurements are then generated and Gaussian noise is added. Attitudes for filter initialization are computed from the first pair of noisy measurement vectors using the triad of vectors (TRIAD) algorithm [27]. Figures 8-10 compare a standard EKF to one using a variational integrator to perform its state prediction step. The two filters have nearly identical performance at high sample rates but show very different behaviors as the sample rate decreases. The underlying reason for this performance difference is the quality of the linearization that the variational integrator yields. Equation (38) and the corresponding Jacobian of Eq. (36) lead to the true linearization of the map from φ_k to φ_k+1.

Figure 9: Root-mean-square attitude error from a multiplicative EKF using standard vs. variational integrators, highlighting improved robustness at lower sampling rates.

Figure 10: Total energy response for a damped rigid body, demonstrating efficient capture of dissipative dynamics by the variational integrator.

$\frac{\partial {\varphi }_{k+1}}{\partial {\varphi }_k}\mathrm{<mo>=</mo>}{\left[\frac{\partial {\rho }_k}{\partial {\varphi }_{k+1}}\right]}^{-1}\frac{\partial {\rho }_k}{\partial {\varphi }_k}\text{ (69)}$

This linearization is completely independent of the step size taken. As a result, filters built around variational integrators are highly insensitive to sample rate and can maintain good performance and convergence at much lower rates than standard extended Kalman filters.

10. Conclusion

The integrators developed in this study offer physically realistic momentum and energy behavior while maintaining modest computational costs. This study also highlights the importance of selecting appropriate control points in the quaternion curve and the quadrature points. Such choices lead to increased computational efficiency by enabling larger time steps. The quaternion approach has the potential to become an excellent tool for computing spacecraft orbital trajectories in proximity operations.

Additional problems may also benefit from using this type of integrator for long-term orbital motion propagation, such as the computation of trajectories for small moons around asteroids, orbital trajectories of small landers or non-propelled small satellites, and dust/particle ejection from an asteroid. It is worth noting that research is currently underway in this area, promising a bright future for symplectic implicit methods such as quaternion variational methods [28,29].