A Robust Attitude Tracking Controller for Spacecraft Based on Singularity-Free Quaternion Nonlinear Dynamic Inversion Framework

1. Introduction

The control of spacecraft attitude has evolved significantly, progressing from classical linear techniques to advanced nonlinear architectures capable of meeting the rigorous agility requirements of modern missions. Conventional Proportional-Integral-Derivative (PID) controllers have served as the industry standard due to their implementation simplicity and extensive flight heritage. While effective for set-point regulation, PID designs often struggle with the strong gyroscopic coupling encountered during rapid, large-angle maneuvers. Some researchers have proposed modified PID methods to increase robustness and overcome the phenomenon [1]. To address the multi-variable interactions, Linear Quadratic Regulator (LQR) was introduced, providing a systematic framework for minimizing energy and error for multi-input-multi-output (MIMO) systems [2,3,4,5]. However, both PID and LQR are inherently limited by their linearization around equilibrium points and often need extensive gain scheduling to cover the full flight envelope. To bridge this gap, Linear-Parameter-Varying (LPV) control has emerged as a powerful intermediate solution [6,7], allowing designers to embed nonlinear dynamics into a family of linear models scheduled by state-dependent parameters, thus extending the valid operating range of linear techniques.

Despite these advancements, the nonlinearity of rigid body dynamics - particularly at high angular velocities - has driven research toward purely nonlinear control paradigms. Sliding Mode Control (SMC) has been widely investigated for its robust properties, utilizing discontinuous control laws to force system states onto a predefined sliding manifold [8]. While SMC is theoretically invariant to matched uncertainties, it famously suffers from the “chattering” phenomenon, which can excite high-frequency structural modes and wear down actuators. Another way to address the nonlinearity is nonlinear dynamic inversion (NDI), also known as feedback linearization, which offers a more geometric and physically intuitive alternative. Unlike robust controllers that fight against nonlinearities, NDI exploits the system’s physics to cancel them. By algebraically transforming the nonlinear spacecraft dynamics into a linear time-invariant (LTI) system, NDI enables the application of well-understood linear control laws to the “synthetic” LTI system.

The NDI technique has been widely implemented in aerospace systems, from aircraft flight control [9,10], reentry vehicle [11], reusable rocket [12], and satellites [8,13,14]. A two-loop cascaded flight control architecture is commonly used for attitude control based on the assumption of time-scale separation. In this approach, the angular velocity loop (the inner loop) and the attitude loop (the outer loop) are linearized and controlled separately. This method simplifies the controller design process. Banerjee et al. [14] demonstrate the conventional two-loop structure NDI in the attitude control of a lunar lander, employing Euler angles for attitude representation. However, the use of Euler angles is hindered by a singularity at pitch angles of ±90 degrees, a phenomenon known as gimbal lock. As a result, alternative approaches for attitude representation, such as quaternions, have garnered increasing attention. Zhang et al. [15] design a quaternion-based NDI for fixed-wing UAV flight controller, but still using the conventional two-loop structure.

Unlike the conventional cascaded structure, we can select quaternions as the output and directly linearize them with respect to the input torque using the NDI technique. The angular velocity is implicitly dictated by the quaternion kinematics equation [8]. However, there is a dimension mismatch between the output quaternion, which is represented as a vector in ${\mathbb{R}}^4$, and the input torque, which is a vector in ${\mathbb{R}}^3$. To resolve this issue, previous studies [8,13,16,17] have used only the vector part of the quaternions as the output, while computing the scalar part using the unit-quaternion constraint. It has been noted, however, that a singularity may occur in the matrix inversion during NDI when the scalar part of the quaternions is zero, leading to numerical instability in control processes. Long et al. [16] and Navabi et al. [17] demonstrate that the internal dynamics of the linearized system remain stable as long as the scalar part of the quaternions is non-zero. To address this issue, Bang et al. [8] proposed adding a small value to the scalar part of the quaternions when it approaches zero, thus avoiding singularities. Bhargavapuri et al. [18] focused on quaternion error, defined as the difference between the desired and actual quaternions. They integrated the scalar part of the quaternions into the NDI framework to prevent it from crossing zero.

The NDI can effectively cancel system nonlinearity, but its performance relies heavily on the accuracy of state feedback and the model information. Model discrepancies and external disturbances can impair the effectiveness of the nonlinearity cancellation. To mitigate the sensitivity of conventional exact linearization to model mismatch, recent literature has heavily favored Incremental Nonlinear Dynamic Inversion (INDI) [19,20]. While INDI successfully bypasses deep model dependencies, it strictly requires high-fidelity angular acceleration measurements. In practical attitude control architectures, differentiating noisy IMU signals severely amplifies high-frequency noise, risking actuator saturation. Therefore, returning to an exact, model-based NDI approach augmented with an active disturbance observer remains practically superior for environments where sensor noise is prevalent.

Chen and Yang [21,22] introduced a nonlinear disturbance observer-based control (DOBC) approach to improve the NDI’s ability to manage disturbances and uncertainties resulting from model mismatches. The DOBC consists of two main components: disturbance estimation and compensation. Chen et al. [23] categorized various types of observers based on their application to linear or nonlinear systems, as well as the compensation mechanisms addressing external disturbances or model mismatches.

Consequently, DOBC and Extended State observers (ESO) have seen continued advancement for robust spacecraft and UAV tracking applications [15,24]. These state-of-the-art observer designs provide excellent lumped disturbance rejection. However, integrating these modern observers with quaternion-based NDI presents a distinct challenge. Conventional input-output linearization on the quaternion vector part often induces unwinding or unstable internal zero dynamics, compromising the observer’s stability guarantees.

Previous research indicates that matrix inversion in NDI is related to the control input and the second derivatives of the output function in rotational dynamics. While earlier studies employed Newtonian mechanics to describe these motions, they encountered singularities with the matrix inverse. Therefore, it is worthwhile to explore alternative approaches. Udwadia et al. [25] present a direct method for obtaining Lagrange rotational motion using quaternions. This approach establishes the relationship between torque and the second derivatives of quaternions, while accounting for the unit-quaternion constraints.

In this paper, we derive a quaternion-based NDI using the Lagrange equation of motion, not suffering from the singularity issues observed in previous studies. Crucially, unlike conventional Newtonian approaches that rely on input-output linearization and inadvertently generate internal zero dynamics, our pure Lagrangian derivation achieves an exact input-state linearization. Because the synthetic linear states are chosen to be mathematically identical to the Lagrangian physical states, the required state coordinate transformation constitutes a perfect global identity diffeomorphism. Furthermore, we rigorously prove the complete realization of this LTI framework. By deriving a structurally non-singular input-transformation and satisfying the necessary differential geometry conditions – specifically, proving that the input distribution is inherently involutive and strictly spans the active constraint manifold – we guarantee that the exact input-state linearization is completely realized on the 6-degree-of-freedom holonomic manifold without inducing any hidden zero dynamics.

After eliminating the nonlinearity using NDI, we can implement a linear controller for the input-output linearized system. In this paper, we design a Linear Quadratic Tracking (LQ) controller that incorporates both feedback and feedforward gains to achieve the desired attitude while ensuring the stabilization of the quaternions and their derivatives. Additionally, a first-order filter is introduced before the controller to smooth the reference signal, rather than applying a step signal directly. This approach helps prevent actuator saturation.

To guarantee robustness against exogenous disturbances, we augment the NDI framework with a DOBC. By designing the nonlinear observer gain to respect the geometric properties of the quaternion, we formally prove that the disturbance estimation error is uniformly ultimately bounded (UUB), even in the presence of a non-vanishing disturbance derivative ($\dot{\mathit{d}}\ne 0)$. Furthermore, rather than analyzing the observer and controller in isolation, we construct a composite Lyapunov function to evaluate the complete NDI-DOBC interconnected system. This rigorous stability analysis demonstrates that the total close-loop system remains strictly UUB, mathematically ensuring safe and robust tracking performance without compromising the strictly separated observer-controller software architecture.

Finally, to validate these theoretical contributions, comprehensive numerical simulations are conducted to evaluate the proposed architecture against conventional input-output linearization methods. The comparative analysis explicitly demonstrates that the proposed exact input-state linearization successfully eliminates the discontinuity in the control input derivative at the quaternion singularity (${\mathit{q}}_0=0$), ensuring smooth and physically realizable control efforts during large-angle, shortest-path maneuvers. Furthermore, the simulations are grounded in the physical parameters of an actual spherical air-bearing attitude control testbed. By explicitly modeling the realistic, time-varying external disturbance torque induced by the gravitational imbalance – resulting from the inevitable misalignment between the center of gravity and the center of rotation – we demonstrate that the proposed NDI-DOBC framework achieves exceptional robust tracking performance, fully validating the theoretical UUB stability guarantees.

In summary, the main contributions of this paper are outlined as follows:

• Exact input-state linearization via Lagrangian mechanics: Unlike conventional Newtonian approaches that utilize input-output linearization and inadvertently generate internal zero dynamics, this work derives a quaternion-based NDI strictly utilizing the Udwadia formulation. This achieves an exact input-state linearization on the 6-degree-of-freedom holonomic manifold, eliminating zero dynamics and avoiding the control input derivative discontinuities typically found at ${\mathit{q}}_0=0$ singularity.
• Disturbance observer design: To address dynamic external disturbances, a nonlinear disturbance observer is designed. By selecting an appropriate observer gain that respects the geometric properties of the quaternion, we mathematically prove that the disturbance estimation error is Uniformly Ultimately Bounded (UUB), even in the presence of a non-vanishing disturbance derivative ($\dot{\mathit{d}}\ne 0$).
• Interconnected closed-loop stability proof: Moving beyond isolated observer and controller analysis, a composite Lyapunov function is constructed to evaluate the complete NDI-DOBC architecture. This stability analysis formally guarantees that the interconnected system remains strictly UUB, mathematically validating the separated controller-observer framework.
• Physical testbed simulation and validation: The theoretical framework is validated through comprehensive numerical simulations parameterized by the physical properties of a spherical air-bearing attitude control testbed. The results demonstrate the proposed method’s superiority over the conventional approaches, showcasing exceptional tracking robustness against realistic, time-varying gravitational imbalance torques.

2. System Architecture

This section introduces the schematic of the NDI-DOBC architecture and the attitude control testbed. First, we present a system overview in the form of a block diagram, providing a visual summary of the framework.

2.1. System Block Diagram

The NDI-DOBC block diagram is shown in Figure 1. For clarity, we delineate the signal propagation from the reference input to the system output. The design process is divided into signal filtering, linear quadratic (LQ) tracking control, nonlinear-dynamic-inversion (NDI) linearization, and disturbance observation. This graphical mapping connects the general overview to each subsystem described below.

A first-order filter, acting as a signal smoother, is placed before the controller to smooth the input signal. For a large-angle maneuver, this smoother outputs a continuous curve to the controller instead of a step function. This prevents excitation of the nonlinearity and actuator saturation from rapid movement. The linear-quadratic tracking controller is designed using the linearized system. It stabilizes the system with state-feedback gains and provides feed-forward gains to improve tracking accuracy. The outputs of the LQ tracking controller, known as “virtual inputs,” serve as the inputs to the linearized system. The NDI then converts these “virtual inputs” into real, physically meaningful control torques for rotational motion.

The disturbance observer, shown as an orange block, estimates external disturbance using information from the original nonlinear system. The nonlinear function $\mathit{p}(\mathit{x})$ and the observer gain $\mathit{L}(\mathit{x})$ are to be designed to ensure the disturbance estimate error stays asymptotically stable regardless of the states $\mathit{x}$. The estimated disturbance is then multiplied by the compensation gain $\mathit{\beta }$ and added to the NDI control torques. These serve as inputs to the nonlinear system. Details are given in Section IV.

2.2. Physical Disturbance Modeling and Theoretical Bounds

The primary exogenous disturbance acting on the attitude testbed arises from physical misalignment between the system’s center of gravity (CG) and its center of rotation, as illustrated in Figure 2. Let ${\mathit{r}}_{cg}$ denote the position vector of the CG relative to the center of rotation, and $m$ represent the total mass of the testbed. The resulting gravitational imbalance torque, $\mathit{d}$, is formulated as:

$$\mathit{d}={\mathit{r}}_{cg}\times (m\mathit{g})$$

(1)

Where $\mathit{g}$ represents the local gravitational acceleration vector. Because the direction of ${\mathit{r}}_{cg}$ relative to the inertial gravity vector changes continuously as the testbed rotates, this disturbance torque is inherently state-dependent and time-varying.

However, because the testbed is mounted on a spherical air-bearing, its maximum tilt angle is strictly constrained by the physical limits of the mechanical structure. Consequently, the magnitude of the gravitational imbalance torque is fundamentally upper bounded by an absolute physical maximum, yielding the inequality $||\mathit{d}||\le D_{max}$. Furthermore, since the angular velocity of the system is physically constrained by the finite saturation limits of the control actuators, the time derivative of the disturbance torque is mathematically guaranteed to be strictly bounded, satisfying $||\dot{\mathit{d}}||\le \delta$.

3. Quaternion Rotation Dynamics

This section introduces the Lagrange equation of rotational motion using quaternions, which is essential for deriving NDI in the next section. Udwadia first derived this equation in 2010 using a direct Lagrange approach without relying on the concept of Lagrange multipliers or Newtonian mechanics [25]. The derivation constructs the Lagrangian from energy terms and expresses the unconstrained equation of motion in terms of quaternions, their derivatives, angular velocity, and generalized torque. The unit-quaternion constraint is incorporated by applying the fundamental equation of constrained motion [26], resulting in a constrained equation of motion. The relationship between generalized and physical torque is then established, providing an equation that describes rotational dynamics in terms of quaternions while satisfying the unit-quaternion constraint for a given input torque.

The derivation begins with the Lagrangian, defined as kinetic energy minus potential energy. For rotational motion, only kinetic energy is present, so the Lagrangian is $T=1/2{\mathit{\omega }}^{\top }\widehat{\mathit{J}}\mathit{\omega }$. The body-fixed coordinate axes are assumed to align with the principal axes of the rigid body, resulting in a diagonal moment of inertia matrix $\widehat{\mathit{J}}=diag(J_1,J_2,J_3)$. The Lagrange equation of motion is:

$${\displaystyle \frac{d}{dt}}\left({\displaystyle \frac{\partial \mathit{T}}{\partial \dot{\mathit{q}}}}\right)-{\displaystyle \frac{\partial \mathit{T}}{\partial \mathit{q}}}={\mathsf{\Gamma }}_{\mathit{q}}$$

(2)

Equation (2) is an unconstrained equation of motion. The quaternion vector $\mathit{q}=\left[q_0 q_1 q_2 q_3\right]$ serves as the generalized coordinate, with each element considered independent. Quaternions are divided into the scalar part $q_0$ and the vector part ${\mathit{q}}_v=\left[q_1 q_2 q_3\right]$. The input to Equation (2) is a $4\times 1$ generalized quaternion torque vector ${\mathsf{\Gamma }}_{\mathit{q}}$. Several properties of quaternions are useful for expanding and rearranging Equation (2):

Quaternion time derivatives: The quaternion time derivative is defined as the “quaternion product” of the quaternion and the augmented angular velocity:

$$\dot{\mathit{q}}={\displaystyle \frac{1}{2}}\mathit{q}\otimes \mathit{w}$$

(3)

Here, $\otimes$ denotes the quaternion product, and the augmented angular velocity is $\mathit{w}={\left[0 {\omega }_1 {\omega }_2 {\omega }_3\right]}^{\top }$. Using the properties of quaternions, Equation (3) can be rearranged into two matrix forms [27]:

$$\dot{\mathit{q}}={\displaystyle \frac{1}{2}}{\mathit{E}}^{\top }\mathit{w}$$

(4)

$\mathit{E}$ and $\dot{\mathit{E}}$ represent the orthogonal matrix of quaternions and its time derivative:

$$\mathit{E}=\left[\begin{array}{cccc}{q}_0& q_1& q_2& q_3\\ -q_1& q_0& q_3& -q_2\\ -q_2& -q_3& q_0& q_1\\ -q_3& q_2& -q_1& q_0\end{array}\right]=\left[\begin{array}{c}{\mathit{q}}^{\top }\\ \\ {\mathit{E}}_1\\ \end{array}\right]$$

(5)

From Equation (4), the augmented angular velocity can be expressed as:

$$\mathit{w}=2\mathit{E}\dot{\mathit{q}}$$

(6)

Properties of matrix$\mathit{E}$ and ${\mathit{E}}_1$: ${\mathit{E}}_1$ is a submatrix of $\mathit{E}$ shown in Equation (5). $\mathit{E}$ is an orthogonal matrix, and ${\mathit{E}}_1$ is a $3\times 4$ matrix with orthogonal rows, hence they have the following property:

$$\begin{array}{ccc}\mathit{E}{\mathit{E}}^{\top }={\mathit{I}}_{4\times 4}& ,& {\mathit{E}}_1{\mathit{E}}_1^{\top }={\mathit{I}}_{3\times 3}\end{array}$$

(7)

$\mathit{q}$ lies in its null space because each row of ${\mathit{E}}_1$ is orthogonal to $\mathit{q}$. Multiplying ${\mathit{E}}_1$ by the quaternion vector yields a zero vector:

$${\mathit{E}}_1\mathit{q}=0$$

(8)

From Equations (5) and (8), it follows that:

$$\mathit{E}\mathit{q}=\left[\begin{array}{c}1\\ 0\\ 0\\ 0\end{array}\right]$$

(9)

Substituting the quaternions with their time derivatives in Equation (8), the ${\mathit{E}}_1$ matrix transforms into ${\dot{\mathit{E}}}_1$. Each row of ${\dot{\mathit{E}}}_1$ is orthogonal to $\dot{\mathit{q}}$:

$${\dot{\mathit{E}}}_1\dot{\mathit{q}}=0$$

(10)

Equation (10) will appear later for the observer gain design.

Unit-quaternion constraint: When quaternions represent attitude, they must satisfy the unit-norm constraint: ${∥\mathit{q}∥}_2=1$. This leads to the following expression:

$$N(\mathit{q})={\mathit{q}}^{\top }\mathit{q}=1$$

(11)

Differentiate Equation (9) with respect to time, and substitute the result $\dot{\mathit{E}}\mathit{q}=-\mathit{E}\dot{\mathit{q}}$ into Equation (6), this yields another expression of augmented angular velocity:

$$\mathit{w}=-2\dot{\mathit{E}}\mathit{q}$$

(12)

Substituting Equations (6) and (12) into the kinetic energy equation, the rotational kinetic energy can be expressed in terms of quaternions and their derivatives:

$$T={\displaystyle \frac{1}{2}}{\mathit{\omega }}^{\top }\widehat{\mathit{J}}\mathit{\omega }={\displaystyle \frac{1}{2}}{\mathit{w}}^{\top }\mathit{J}\mathit{w}=2{\dot{\mathit{q}}}^T{\mathit{E}}^{\top }\mathit{J}\mathit{E}\dot{\mathit{q}}=2{\mathit{q}}^{\top }{\dot{\mathit{E}}}^{\top }\mathit{J}\dot{\mathit{E}}\mathit{q}$$

(13)

$\widehat{\mathit{J}}$ is a $3\times 3$ moment of inertia matrix, and $\mathit{J}=diag(J_0,J_1,J_2,J_3)$ is a $4\times 4$ extended moment of inertia matrix matching the dimension of $\mathit{w}$, where $J_0$ is an arbitrary positive number. The unconstrained Lagrange equation of motion (2) can also be rearranged as follows:

$$4{\mathit{E}}^T\mathit{J}\dot{\mathit{E}}\ddot{\mathit{q}}+8{\dot{\mathit{E}}}^{\top }\mathit{J}\mathit{E}\dot{\mathit{q}}+4J_{0}N(\dot{\mathit{q}})\mathit{q}={\mathsf{\Gamma }}_{\mathit{q}}$$

(14)

Applying the fundamental equation of constrained motion [26] incorporates the unit-quaternion constraint into the unconstrained system. After algebraic manipulation, the arbitrary number $J_0$ is eliminated, and the generalized acceleration $\ddot{\mathit{q}}$ is given by Equation (14) as

$$\ddot{\mathit{q}}=-{\displaystyle \frac{1}{2}}{\mathit{E}}_1^{\top }{\widehat{\mathit{J}}}^{-1}{\mathit{\omega }}^{\times }\widehat{\mathit{J}}\mathit{\omega }-N(\dot{\mathit{q}})\mathit{q}+{\displaystyle \frac{{\mathit{E}}_1^{\top }{\widehat{\mathit{J}}}^{-1}{\mathit{E}}_1{\mathsf{\Gamma }}_{\mathit{q}}}{4}}$$

(15)

${\mathit{\omega }}^{\times }$ in Equation (15) denotes the skew-symmetric matrix of angular velocity, defined as:

$${\mathit{\omega }}^{\times }=\left[\begin{array}{ccc}0& -{\omega }_3& {\omega }_2\\ {\omega }_3& 0& -{\omega }_1\\ -{\omega }_2& {\omega }_1& 0\end{array}\right]$$

(16)

For practical applications, the generalized quaternion torque vector must be converted to a physical torque ${\mathsf{\Gamma }}_B$. Udwadia [25] established the relationship between these torques using the “virtual work method,” resulting in:

$${\mathsf{\Gamma }}_{\mathit{q}}=2{\mathit{E}}_1^{\top }{\mathsf{\Gamma }}_B$$

(17)

By substituting Equation (17) into Equation (15), the Lagrange equation of motion for rigid body rotation using quaternions and body-fixed frame torque is obtained:

$$\ddot{q}=-{\displaystyle \frac{1}{2}}{\mathit{E}}_1^{\top }{\widehat{\mathit{J}}}^{-1}{\mathit{\omega }}^{\times }\widehat{\mathit{J}}\mathit{\omega }-N(\dot{\mathit{q}})\mathit{q}+{\displaystyle \frac{{\mathit{E}}_1^{\top }{\widehat{\mathit{J}}}^{-1}{\mathsf{\Gamma }}_B}{2}}$$

(18)

4. Quaternion NDI and LQ Tracking Controller

In this section, we apply nonlinear dynamic inversion (NDI) to the rotational dynamics – based on the Lagrange equations of motion presented in Equation (18) – to obtain a linearized dynamic model. Following this, a linear quadratic (LQ) tracking controller is developed for the linearized system, utilizing both state feedback and feedforward gains to ensure precise tracking performance.

To simplify controller design, we apply exact input-state linearization to transform the nonlinear rotational dynamics into synthetic Linear Time-Invariant (LTI) system. By comparing the system equations with the Lagrange equation of motion from Section 3, we define the complete system state vector using the quaternions and their time derivatives, $\mathit{x}={\left[\begin{array}{cc}{\mathit{q}}^{\mathit{\top }}& {\dot{\mathit{q}}}^{\mathit{\top }}\end{array}\right]}^{\mathit{\top }}\in {\mathbb{R}}^8$.

Equations (4) and (18) can be rewritten as follows:

$$\dot{\mathit{x}}(t)=\mathit{f}(\mathit{x}(t))+\mathit{g}(\mathit{x}(t))\mathit{u}$$

(19)

where the drift vector field $\mathit{f}(\mathit{x})$ and the input matrix $\mathit{g}(\mathit{x})$ are partitioned as:

$$\begin{array}{c}f(x)=\left[\begin{array}{c}{\mathit{f}}_u\\ {\mathit{f}}_l\end{array}\right]=\left[\begin{array}{c}{\displaystyle \frac{1}{2}}{\mathit{E}}^{\top }\mathit{w}\\ -{\displaystyle \frac{1}{2}}{\mathit{E}}_1^{\top }{\widehat{\mathit{J}}}^{-1}{\mathit{\omega }}^{\times }\widehat{\mathit{J}}\mathit{\omega }-N(\dot{\mathit{q}})\mathit{q}\end{array}\right],g(x)=\left[\begin{array}{c}{\mathit{g}}_u\\ {\mathit{g}}_l\end{array}\right]=\left[\begin{array}{c}{0}_{4\times 3}\\ {\displaystyle \frac{1}{2}}{\mathit{E}}_1^{\top }{\widehat{\mathit{J}}}^{-1}\end{array}\right]\end{array}$$

(20)

Here, the angular velocity vector $\mathit{\omega }$ is treated as a parameter, with its relationship to $\dot{\mathit{q}}$ given in Equation (4). If the quaternion time derivatives are bounded, the angular velocity will also remain bounded.

The control inputs are the torques applied to the body axes, $\mathit{u}={\mathsf{\Gamma }}_B$. To perform the input-state linearization, we examine the lower partition of the state equations, which governs the angular acceleration ($\ddot{\mathit{q}}$):

$$\ddot{\mathit{q}}={\mathit{f}}_l\left(\mathit{x}\right)+{\mathit{g}}_l\left(\mathit{x}\right){\mathsf{\Gamma }}_B$$

(21)

By defining the second derivative of the quaternions as the “virtual input”, $\mathit{\nu }=\ddot{\mathit{q}}\in {\mathbb{R}}^4$, we can algebraically cancel the system nonlinearities. Let $\mathit{F}\left(\mathit{x}\right)={\mathit{f}}_l(\mathit{x})$ and $\mathit{G}\left(\mathit{x}\right)={\mathit{g}}_l\left(\mathit{x}\right)$. The acceleration dynamics become:

$$\ddot{\mathit{q}}=\mathit{F}\left(\mathit{x}\right)+\mathit{G}\left(\mathit{x}\right){\mathsf{\Gamma }}_B=\mathit{\nu }$$

(22)

Utilizing the Moore-Penrose pseudoinverse of the decoupling matrix $\mathit{G}(\mathit{x})$, we derive the input transformation mapping the virtual control $\mathit{\nu }$ to the physical torque ${\mathsf{\Gamma }}_{\mathrm{B}}$:

$${\mathsf{\Gamma }}_B={\mathit{G}}^{†}\left(\mathit{x}\right)\left(\mathit{\nu }-\mathit{F}\left(\mathit{x}\right)\right)$$

(23)

Expanding the Equation (23) yields:

$${\mathsf{\Gamma }}_B=({\mathit{E}}_1^{\top }{\widehat{\mathit{J}}}^{-1}{)}^{†}(2\mathit{\nu }+{\mathit{E}}_1^{\top }{\widehat{\mathit{J}}}^{-1}{\mathit{\omega }}^{\times }\widehat{\mathit{J}}\mathit{\omega }+2({\dot{\mathit{q}}}^{\top }\dot{\mathit{q}})\mathit{q})$$

(24)

Since ${\mathit{E}}_1^{\top }{\widehat{\mathit{J}}}^{-1}$ is not a square matrix, we use the pseudoinverse of $\mathit{G}(\mathit{x})$. In the next subsection, we examine whether a singularity exists in this NDI formulation, as observed in previous studies [13,16,17,18].

Finally, the linearization system is shown below:

$${\displaystyle \frac{d}{dt}}\mathit{\xi }=\left[\begin{array}{cc}{0}_{4\times 4}& {\mathit{I}}_{4\times 4}\\ 0_{4\times 4}& 0_{4\times 4}\end{array}\right]\mathit{\xi }+\left[\begin{array}{c}{0}_{4\times 4}\\ {\mathit{I}}_{4\times 4}\end{array}\right]\mathit{\nu }$$

(25)

where $\mathit{\xi }=[{\mathit{q}}^{\top } {\dot{\mathit{q}}}^{\top }{]}^{\top }$ is the $8\times 1$ states vector and $\mathit{\nu }$ is the input vector of the LTI system.

4.1. Singularity Analysis of the Matrix Inversion

The matrix ${\mathit{E}}_1^{\top }{\widehat{\mathit{J}}}^{-1}$ in Equation (24) is a $4\times 3$ matrix. To determine whether a pseudoinverse exists, we examine its column rank. If the rank is 3, the pseudoinverse exists, and no singularity occurs during inversion. To show that ${\mathit{E}}_1^{\top }{\widehat{\mathit{J}}}^{-1}$ is full column rank, it is equivalent showing that the matrix $({\mathit{E}}_1^{\top }{\widehat{\mathit{J}}}^{-1}{)}^T({\mathit{E}}_1^{\top }{\widehat{\mathit{J}}}^{-1})$ is full rank and invertible.

$$\begin{array}{c}{\left({\mathit{E}}_1^{\top }{\widehat{\mathit{J}}}^{-1}\right)}^{\top }\left({\mathit{E}}_1^{\top }{\widehat{\mathit{J}}}^{-1}\right)={\widehat{\mathit{J}}}^{-1}{\mathit{E}}_1{\mathit{E}}_1^{\top }{\widehat{\mathit{J}}}^{-1}={\widehat{\mathit{J}}}^{-2}\\ \\ =\left[\begin{array}{ccc}{\displaystyle \frac{1}{J_1^2}}& 0& 0\\ 0& {\displaystyle \frac{1}{J_2^2}}& 0\\ 0& 0& {\displaystyle \frac{1}{J_3^2}}\end{array}\right]\end{array}$$

(26)

Since the result in equation (26) is a diagonal matrix and the moments of inertia of the principal axes are always positive, ${\mathit{E}}_1^{\top }{\widehat{\mathit{J}}}^{-1}$ is full column rank and its pseudoinverse exists.

The pseudoinverse of ${\mathit{E}}_1^{\top }{\widehat{\mathit{J}}}^{-1}$ is referred to as a “left inverse” because it is a tall matrix. The inverse is provided below:

$$\begin{array}{c}({\mathit{E}}_1^{\top }{\widehat{\mathit{J}}}^{-1}{)}^{†}={\left[({\mathit{E}}_1^{\top }{\widehat{\mathit{J}}}^{-1}{)}^T({\mathit{E}}_1^{\top }{\widehat{\mathit{J}}}^{-1})\right]}^{-1}({\mathit{E}}_1^{\top }{\widehat{\mathit{J}}}^{-1}{)}^{\top }\\ =\left[\begin{array}{cccc}-J_1{q}_1& J_1{q}_0& J_1{q}_3& -J_1{q}_2\\ -J_2{q}_2& -J_2{q}_3& J_2{q}_0& J_2{q}_1\\ -J_3{q}_3& J_3{q}_2& -J_3{q}_1& J_3{q}_0\end{array}\right]\end{array}$$

(27)

In previous works [16,17] that derive the NDI formulation using a Newtonian mechanics approach, the counterpart matrix $\mathit{G}(\mathit{x})$ to be inverted has the form:

$$\mathit{G}(\mathit{x})={\displaystyle \frac{1}{2}}\left[\begin{array}{ccc}{\displaystyle \frac{q_0}{J_1}}& -{\displaystyle \frac{q_3}{J_2}}& {\displaystyle \frac{q_2}{J_3}}\\ {\displaystyle \frac{q_3}{J_1}}& {\displaystyle \frac{q_0}{J_2}}& -{\displaystyle \frac{q_1}{J_3}}\\ -{\displaystyle \frac{q_2}{J_1}}& {\displaystyle \frac{q_1}{J_2}}& {\displaystyle \frac{q_0}{J_3}}\end{array}\right]$$

(28)

With its determinant

$$\det \left(\mathit{G}\left(\mathit{x}\right)\right)={\displaystyle \frac{q_0}{8J_1{J}_2{J}_3}}$$

(29)

The result shows a singularity at $q_0=0$, which leads to numerical instability when calculating the control torque.

In summary, the NDI formulation derived from the Lagrange mechanics approach can include both the vector and scalar parts of quaternions as outputs. Importantly, this approach does not exhibit the singularity reported in previous work, providing a robust foundation for controller design.

4.2. Proof of Exact Input-State Linearization

To formally establish that the proposed Lagrangian NDI constitutes an exact input-state linearization rather than a conventional input-output linearization, the geometric condition of the state transformation must be rigorously verified. In standard Newtonian formulations, linearizing the rotational dynamics typically involves selecting the vector part of the quaternion as the system output. This approach dictates a coordinate transformation that inevitably splits the state space, leaving the scalar part of the quaternion as internal zero dynamics.

Conversely, the Udwadia formulation operates directly on the generalized coordinates while intrinsically preserving the holonomic constraint (${\mathit{q}}^{\top }\mathit{q}=1$) without requiring state partitioning. Because the Lagrangian derivation comprehensively maps the physical constraints into system’s kinetic energy, the nonlinear coordinate transformation $\mathit{\xi }=T(\mathit{x})$ required to map the original dynamics into the synthetic LTI space is simply the identity mapping:

$$\mathit{\xi }=\mathit{x}$$

(30)

For exact input-state linearization to be theoretically valid, the coordinate transformation $T(\mathit{x})$ must be a diffeomorphism. The Jacobian matrix of the proposed transformation evaluates trivially to the identity matrix:

$${\displaystyle \frac{\partial T\left(\mathit{x}\right)}{\partial \mathit{x}}}=\mathit{I}$$

(31)

Since the identity matrix is globally non-singular and its inverse is smooth, the mapping $\mathit{\xi }=\mathit{x}$ constitutes a global diffeomorphism. This mathematical property is highly significant: it guarantees that the entire mathematical state vector is seamlessly and invertibly mapped into the synthetic linear domain. Consequently, the entire state space is fully observable and controllable within the linear framework, rigorously proving the complete absence of internal zero dynamics.

With the global diffeomorphism established, the second requirement for input-state linearization is a well-defined, non-singular input mapping. As derived in the beginning of Section 4, the algebraic mapping between the virtual control input $\mathit{\nu }$ and the physical torque ${\mathsf{\Gamma }}_B$ is given by ${\mathsf{\Gamma }}_B={\mathit{G}}^{†}\left(\mathit{x}\right)\left(\mathit{\nu }-\mathit{F}\left(\mathit{x}\right)\right)$. Because the decoupling matrix $\mathit{G}\left(\mathit{x}\right)$ is analytically proven to maintain full column rank across the entire operational manifold, its pseudoinverse ${\mathit{G}}^{†}\left(\mathit{x}\right)$ is globally well-defined. Therefore, this input transformation strictly guarantees the exact algebraic cancellation of the nonlinear dynamics without encountering the mathematical singularities inherent to conventional Newtonian formulations.

Finally, to eliminate any theoretical ambiguity regarding the realization of this LTI system, the proposed transformation must satisfy the necessary and sufficient differential geometry conditions established by the Frobenius Theorem [28]. Let the $8\times 3$ system input matrix $\mathit{g}(\mathit{x})$ be partitioned into its three constituent column vector fields, such that $\mathit{g}\left(x\right)=\left[\begin{array}{ccc}{\mathit{g}}_1\left(\mathit{x}\right)& {\mathit{g}}_2\left(\mathit{x}\right)& {\mathit{g}}_3\left(\mathit{x}\right)\end{array}\right]$, where each ${\mathit{g}}_i\left(\mathit{x}\right)\in {\mathbb{R}}^8$ corresponds to a distinct control input channel. For the feedback linearization of our system, the input distribution is defined by the span of these column vector fields:

$${\mathsf{\Delta }}_0=span\left\{{\mathit{g}}_1,{\mathit{g}}_2,{\mathit{g}}_3\right\}$$

(32)

For exact linearization to be mathematically realized, this distribution ${\mathsf{\Delta }}_0$ must be involutive. Because the input vector fields ${\mathit{g}}_i\left(\mathit{x}\right)$ are structurally independent of the angular velocities $\dot{\mathit{q}}$, the Jacobian cross-multiplications strictly evaluate to zero. Consequently, the Lie brackets of any two input vector fields commute and vanish ($\left[g_i,g_j\right]=0$), proving ${\mathsf{\Delta }}_0$ is completely involutive.

Furthermore, to evaluate the multi-input-multi-output exact linearization, the “controllability-like” block matrix $\mathcal{g}\left(\mathit{x}\right)$ must be constructed using the Lie brackets between the drift vector field $\mathit{f}(\mathit{x})$ and the input distribution. Because the physical control inputs are torques acting on rotational dynamics, the relative degree of each input channel to the attitude state is $r_i=2$ for $i=\mathrm{1,2},3$. According to feedback linearization theory, the maximum order of the required Lie brackets in $\mathcal{g}(\mathit{x})$ matrix is $r_i-1=1$. Thus the required block matrix is formulated as:

$$\mathcal{g}\left(\mathit{x}\right)=\left[\begin{array}{cccccc}{\mathit{g}}_1& {\mathit{g}}_2& {\mathit{g}}_3& a{d}_{\mathit{f}}{\mathit{g}}_1& a{d}_{\mathit{f}}{\mathit{g}}_2& a{d}_{\mathit{f}}{\mathit{g}}_3\end{array}\right]$$

(33)

Let $D={\mathbb{R}}^8$ define the unconstrained domain of the mathematical vector $\mathit{x}={\left[\begin{array}{cc}{\mathit{q}}^{\top }& {\dot{\mathit{q}}}^{\top }\end{array}\right]}^{\top }$. In conventional unconstrained systems, exact linearization would require $\mathcal{g}(\mathit{x})$ to have a rank of 8. However, governed by the Udwadia formulation, the system’s actual physical trajectories are permanently restricted to the active holonomic constraint manifold defined by ${\mathit{q}}^{\top }\mathit{q}=1$ and ${\mathit{q}}^{\top }\dot{\mathit{q}}=0$. Let this 6-degree-of-freedom constraint manifold be defined as the specific operational subdomain $D_0\subset D$.

As established by the full column rank of the decoupling matrix $\mathit{G}(\mathit{x})$, the block matrix $\mathcal{g}(\mathit{x})$ strictly possesses a rank of 6 for all $\mathit{x}\in D_0$. Because the sum of the relative degrees ($\sum r_i=6$) precisely equals the topological dimension of the active manifold $D_0$, the exact full input-state linearization is completely realized on this constraint manifold. This rigorously proves that within the physically valid domain $D_0$, the transformation achieves perfect LTI realization without generating any internal zero dynamics.

4.3. Linear Quadratic Tracking Controller

After the NDI formulation, a linear quadratic controller is designed to determine the virtual input in Equation (25). In Figure 1, the blue area represents the linearized equivalent system, and the LQ controller provides the virtual input $\mathit{\nu }$ to control the linearized state $\mathit{\xi }$ and track the reference signal ${\mathit{\xi }}_{ref}$.

A linear quadratic controller for attitude tracking is formulated to minimize the cost functional:

$$J={\int }_0^{\mathrm{\infty }}\left((\mathit{\xi }-{\mathit{\xi }}_{\mathit{r}\mathit{e}\mathit{f}}{)}^{\top }\mathit{Q}(\mathit{\xi }-{\mathit{\xi }}_{\mathit{r}\mathit{e}\mathit{f}})+{\mathit{\nu }}^{\top }\mathit{R}\mathit{\nu }\right)dt$$

(34)

and subject to the linearized dynamic system (25). $\mathit{Q}$ and $\mathit{R}$ are positive definite matrices representing the weighting matrices for the linear system’s states and virtual inputs. The control input vector computed by the linear quadratic method takes the form: $\mathit{\nu }=-\mathit{K}\left(\mathit{\xi }-{\mathit{\xi }}_{\mathit{r}\mathit{e}\mathit{f}}\right)=-\mathit{K}\mathit{\xi }+\mathit{r}$, where $\mathit{K}$ is the state feedback gain and $\mathit{r}$ is the feedforward term [29].

The state feedback gain is obtained by solving the Riccati equation. In our case, it can be partitioned into two submatrices, $\mathit{K}=\left[{\mathit{K}}_q {\mathit{K}}_{\dot{\mathrm{q}}}\right]$, corresponding to their respective state variables. When tracking a desired attitude, ${\mathit{q}}_d$, we set the desired state vector as ${\mathit{x}}_d={\left[\begin{array}{cc}{\mathit{q}}_d^{\top }& 0_{4\times 1}^{\mathit{\top }}\end{array}\right]}^{\top }$, and the 1st-order filter generates smooth sequence of equilibrium points ${\mathit{\xi }}_{ref}=[{\mathit{q}}_{ref}^{\top } 0_{4\times 1}^{\top }{]}^{\top }$ as reference signal for the controller instead of a step signal. The feedforward term is then $\mathit{r}={\mathit{K}}_q{\mathit{q}}_{ref}$. To summarize, the virtual control input consists of both feedback and feedforward components, as shown below.

$$\mathit{\nu }=-\left[{\mathit{K}}_q {\mathit{K}}_{\dot{q}}\right]\mathit{\xi }+{\mathit{K}}_q{q}_{ref}$$

(35)

5. Nonlinear Disturbance Observer-Based Controller

Nonlinear systems often encounter unmodeled disturbances, which can significantly impact the performance of model-based controllers. These disturbances may originate from the external environment, parameter uncertainty, or unmodeled dynamics. Chen et al. [23] classified such disturbances as either matched or mismatched. In our application, external disturbances result from shifts in the center of gravity and the center of the spherical air-bearing, producing additional torque on the attitude control testbed. Because this disturbance acts through the same channel as the control torque, it is considered as a matched disturbance. The rotation dynamics with external disturbances can be represented as follows:

$$\dot{\mathit{x}}=\mathit{f}(\mathit{x})+{\mathit{g}}_1(\mathit{x})\mathit{u}+{\mathit{g}}_2(\mathit{x})\mathit{d}$$

(36)

The state vector, $\mathit{x}=\left[\mathit{q} \dot{\mathit{q}}\right]$, matches the nonlinear system in equation (19). ${\mathit{g}}_1(\mathit{x})$ corresponds to the $\mathit{g}(\mathit{x})$ matrix in equation (20), while ${\mathit{g}}_2(\mathit{x})$ maps disturbances $\mathit{d}$ into the state space. For matched disturbances, we have ${\mathit{g}}_1\left(\mathit{x}\right)={\mathit{g}}_2\left(\mathit{x}\right)=\mathit{g}(\mathit{x})$.

In the DOBC framework, the controller and disturbance observer loops are designed independently. This separation principle enables the direct application of the NDI and LQ controller results from the previous section. We first introduce the DOBC framework to propose a nonlinear observer gain function based on the Lagrange equation of motion using quaternions, and then combine the control torque from the LQ controller with the estimated disturbance to obtain the final control inputs for the nonlinear system (38). The orange block in Figure 1 illustrates the disturbance observer structure. The estimated disturbance is fed forward to the control input through the triangular gain block $\beta$, located between the blue and orange blocks.

Chen et al. [23] reviewed various disturbance observers and their corresponding DOBCs, considering the characteristics of the target system and disturbances. For nonlinear systems with matched disturbances, they propose a nonlinear disturbance observer:

$$\left\{\begin{array}{c}\dot{\mathit{z}}=-L(x)g(x)z-L(x)\left[\mathit{g}(\mathit{x})\mathit{p}(\mathit{x})+\mathit{f}(\mathit{x})+\mathit{g}(\mathit{x})\mathit{u}\right]\\ \widehat{\mathit{d}}=z+p(x)\end{array}\right.$$

(37)

The observer output, $\widehat{\mathit{d}}$ in the lower part of Equation (37), represents the estimated disturbance. It is calculated as the sum of the observer’s internal states, $\mathit{z}$, and a nonlinear function, $\mathit{p}(\mathit{x})$, which must be designed. The upper part of Equation (37) describes the dynamics of the internal states, which depend on the nonlinear system (38). The relationship between observer gain function $\mathit{L}(\mathit{x})$ and $\mathit{p}(\mathit{x})$ is given below:

$$\mathit{L}(\mathit{x})={\displaystyle \frac{\partial \mathit{p}(\mathit{x})}{\partial \mathit{x}}}$$

(38)

The observer gain function $\mathit{L}(\mathit{x})$ is synthesized based on the nominal disturbance assumption ($\dot{\mathit{d}}=0$). Under this condition, the nominal disturbance estimation error dynamics are defined as:

$${\dot{\mathit{e}}}_{\mathit{d}}(t)+\mathit{L}(\mathit{x})\mathit{g}(\mathit{x}){\mathit{e}}_d(t)=0$$

(39)

Equation (39) is a first-order ODE. If the designed gain ensures that $\mathit{L}\left(\mathit{x}\right)\mathit{g}(\mathit{x})$ is positive definite, the nominal error dynamics are asymptotically stable. The disturbance estimation error is defined as the difference between the actual disturbance and its estimate, ${\mathit{e}}_d=\mathit{d}-\widehat{\mathit{d}}$. Therefore, when designing the disturbance observer, we first determine $\mathit{L}(\mathit{x})$, then obtain $\mathit{p}(\mathit{x})$ by integrating $\mathit{L}(\mathit{x})$. The following subsection discusses the details of designing the observer gain function.

Equation (39) is a first-order ODE. If $\mathit{L}(\mathit{x})\mathit{g}(\mathit{x})$ is positive definite, the ODE is stable. The gain function can be defined as $\mathit{L}(\mathit{x})=\left[{\mathit{L}}_1 {\mathit{L}}_2\right]$, so $\mathit{L}\left(\mathit{x}\right)\mathit{g}\left(\mathit{x}\right)={\mathit{L}}_1{\mathit{g}}_u+{\mathit{L}}_2{\mathit{g}}_l$. Since $\mathit{g}(\mathit{x})=[0_{4\times 3} {\displaystyle \frac{1}{2}}{\mathit{E}}_1^{\top }{\widehat{\mathit{J}}}^{-1}{]}^{\top }$, and ${\mathit{L}}_1{\mathit{g}}_u=0$, equation (39) is equivalent to:

$${\dot{\mathit{e}}}_{\mathit{d}}(t)+{\mathit{L}}_2\left(\mathit{x}\right){\mathit{g}}_l(\mathit{x}){\mathit{e}}_d(t)=0$$

(40)

Therefore, we can design ${\mathit{L}}_2$ first by setting ${\mathit{L}}_2{\mathit{g}}_l\succ 0$. To guarantee that ${\mathit{L}}_2{\mathit{g}}_l$ be positive definite, we can simply let it be an identity matrix multiplied by a positive-scalar gain $K_{obs}$. Learning that ${\mathit{g}}_l={\displaystyle \frac{1}{2}}{\mathit{E}}_1^{\top }{\widehat{\mathit{J}}}^{-1}$, and the property in Equation (7),we can set ${\mathit{L}}_2=K_{obs}\widehat{\mathit{J}}{\mathit{E}}_1$ so that

$${\mathit{L}}_2\left(\mathit{x}\right){\mathit{g}}_l\left(\mathit{x}\right)={\displaystyle \frac{1}{2}}{K}_{obs}\widehat{\mathit{J}}{\mathit{E}}_1{\mathit{E}}_1^{\top }{\widehat{\mathit{J}}}^{-1}={\displaystyle \frac{1}{2}}{K}_{obs}{\mathit{I}}_{3\times 3}$$

(41)

It is clear from Equation (38) that ${\mathit{L}}_2$ is the partial derivative of $\mathit{p}(\mathit{x})$ with respect to $\dot{\mathit{q}}$. Expanding ${\mathit{L}}_2$ yields:

$${\mathit{L}}_2\left(\mathit{x}\right)={\displaystyle \frac{\partial \mathit{p}(\mathit{x})}{\dot{\mathit{q}}}}=K_{obs}\left[\begin{array}{cccc}-J_1{q}_1& J_1{q}_0& J_1{q}_3& -J_1{q}_2\\ -J_2{q}_2& -J_2{q}_3& J_2{q}_0& J_2{q}_1\\ -J_3{q}_3& J_3{q}_2& -J_3{q}_1& J_3{q}_0\end{array}\right]$$

(42)

The nonlinear function $\mathit{p}(\mathit{x})$ can be obtained by integrating ${\mathit{L}}_2$ with respect to $\dot{\mathit{q}}$. Assuming a zero constant of integration for simplicity, which leads to:

$$\mathit{p}(\mathit{x})=K_{obs}\left[\begin{array}{c}{J}_1(-q_1\dot{q_0}+q_0\dot{q_1}+q_3\dot{q_2}-q_2\dot{q_3})\\ J_2(-q_2\dot{q_0}-q_3\dot{q_1}+q_0\dot{q_2}+q_1\dot{q_3})\\ J_3(-q_3\dot{q_0}+q_2\dot{q_1}-q_1\dot{q_2}+q_0\dot{q_3})\end{array}\right]$$

(43)

And ${\mathit{L}}_1$ can be derived by partially differentiating $\mathit{p}(\mathit{x})$ with respect to $\mathit{q}$.

$${\mathit{L}}_1\left(\mathit{x}\right)={\displaystyle \frac{\partial \mathit{p}(\mathit{x})}{\mathit{q}}}=K_{obs}\left[\begin{array}{cccc}{J}_1{\dot{q}}_1& {-J}_1{\dot{q}}_0& -J_1{\dot{q}}_3& J_1{\dot{q}}_2\\ J_2{\dot{q}}_2& J_2{\dot{q}}_3& {-J}_2{\dot{q}}_0& -J_2{\dot{q}}_1\\ J_3{\dot{q}}_3& -J_3{\dot{q}}_2& J_3{\dot{q}}_1& -J_3{\dot{q}}_0\end{array}\right]=-K_{obs}\widehat{\mathit{J}}{\dot{\mathit{E}}}_1$$

(44)

While the condition ${\mathit{L}}_1{\mathit{g}}_u=0$ does not affect Equation (39) and allows us to design ${\mathit{L}}_2$ first, the $\dot{\mathit{z}}$ dynamics in Equation (37) includes the term $\mathit{L}\left(\mathit{x}\right)\mathit{f}\left(\mathit{x}\right)={\mathit{L}}_1{\mathit{f}}_u+{\mathit{L}}_2{\mathit{f}}_l$. Since ${\mathit{f}}_u\ne 0$, it is necessary to investigate the influence of ${\mathit{L}}_1$ on $\dot{\mathit{z}}$. From the property described by Equation (10), it follows that:

$${\mathit{L}}_1{\mathit{f}}_u=-K_{obs}\widehat{\mathit{J}}{\dot{\mathit{E}}}_1\left({\displaystyle \frac{1}{2}}{\mathit{E}}^{\top }\mathit{w}\right)=-K_{obs}\widehat{\mathit{J}}{\dot{\mathit{E}}}_1\dot{\mathit{q}}=0_{3\times 1}$$

(45)

Given that ${\mathit{L}}_1$ has no influence on the dynamics of $\mathit{z}$, as shown in Equation (47), the gain $K_{obs}$ becomes negligible. Thus, the observer gain function $\mathit{L}(\mathit{x})$ is constructed as:

$$\mathit{L}\left(\mathit{x}\right)=\left[\begin{array}{cc}-\widehat{\mathit{J}}{\dot{\mathit{E}}}_1& K_{obs}\widehat{\mathit{J}}{\mathit{E}}_1\end{array}\right]$$

(46)

After determining $\mathit{L}(\mathit{x})$ and $\mathit{p}(\mathit{x})$, the estimated disturbance can be obtained from Equation (37). The system input $\mathit{u}$ should not be taken directly from the control torque calculated by the NDI and LQ controller. Disturbance compensation is applied before the control torque is used in the nonlinear system:

$$\mathit{u}={\mathit{\Gamma }}_{\mathit{B}}+\mathit{\beta }\widehat{\mathit{d}}$$

(47)

The $\beta$ in the equation above is the disturbance compensation gain. Chen et al. [23] provided various methods for determining this gain. In the case of matched disturbances, where ${\mathit{g}}_1(\mathit{x})={\mathit{g}}_2(\mathit{x})$, the gain reduces to -1. Therefore, the control input $\mathit{u}$ for the nonlinear system becomes:

$$\mathit{u}={\mathit{\Gamma }}_{\mathit{B}}-\widehat{\mathit{d}}$$

(48)

5.1. Lyapunov Stability Analysis of the Disturbance Observer

Having established the observer gain matrix $\mathit{L}(\mathit{x})$ to stabilize the nominal system, we now evaluate the actual disturbance estimation error in the presence of physical, time-varying disturbances. In reality, the gravitational imbalance torque possesses a non-vanishing derivative ($\dot{\mathit{d}}\ne 0$). Therefore, the true error dynamics are governed by:

$${\dot{\mathit{e}}}_{\mathit{d}}(t)+\mathit{L}(\mathit{x})\mathit{g}(\mathit{x}){\mathit{e}}_d(t)=\dot{\mathit{d}}(t)$$

(49)

Consider the following positive-definite Lyapunov candidate function for the observer error dynamics:

$$V_o\left({\mathit{e}}_{\mathit{d}}\right)={\displaystyle \frac{1}{2}}{\mathit{e}}_{\mathit{d}}^{\mathit{\top }}{\mathit{e}}_{\mathit{d}}$$

(50)

Taking the time derivative of $V_0$ along the trajectories of the error dynamic yields:

$${\dot{V}}_o\left({\mathit{e}}_{\mathit{d}}\right)={\mathit{e}}_{\mathit{d}}^{\mathit{\top }}{\dot{\mathit{e}}}_{\mathit{d}}={\mathit{e}}_{\mathit{d}}^{\mathit{\top }}\left(\dot{\mathit{d}}-\mathit{L}\left(\mathit{x}\right)\mathit{g}\left(\mathit{x}\right){\mathit{e}}_{\mathit{d}}\right)$$

(51)

Substituting the derived relation $\mathit{L}\left(\mathit{x}\right)\mathit{g}\left(\mathit{x}\right)=1/2K_{obs}{\mathit{I}}_{3\times 3}$, the derivative becomes:

$${\dot{V}}_o\left({\mathit{e}}_{\mathit{d}}\right)=-{\displaystyle \frac{1}{2}}{K}_{obs}{\mathit{e}}_{\mathit{d}}^{\mathit{\top }}{\mathit{e}}_{\mathit{d}}+{\mathit{e}}_{\mathit{d}}^{\mathit{\top }}\dot{\mathit{d}}$$

(52)

Utilizing the Rayleigh quotient and the Cauchy-Schwarz inequality, we can bound the Lyapunov derivative as follows:

$${\dot{V}}_o\left({\mathit{e}}_{\mathit{d}}\right)\le -{\displaystyle \frac{1}{2}}{K}_{obs}||{\mathit{e}}_{\mathit{d}}|{|}^2+||{\mathit{e}}_{\mathit{d}}||||\dot{\mathit{d}}||$$

(53)

As established in Section 2.2, the time derivative of the gravitational imbalance disturbance is strictly bounded by $||\dot{\mathit{d}}||<\delta$. Applying this physical bound:

$${\dot{V}}_o\left({\mathit{e}}_{\mathit{d}}\right)\le -||{\mathit{e}}_{\mathit{d}}||\left({\displaystyle \frac{1}{2}}{K}_{obs}||{\mathit{e}}_{\mathit{d}}||-\delta \right)$$

(54)

For the Lyapunov derivative to be strictly negative ($V_0<0$), the norm of the estimation error must satisfy:

$$||{\mathit{e}}_{\mathit{d}}||>{\displaystyle \frac{2\delta }{K_{obs}}}$$

(55)

This result mathematically guarantees that the disturbance estimation error is Uniformly Ultimately Bounded (UUB). The error will converge exponentially to an ultimate bound defined by $b_d=2\delta /K_{obs}$, and remain within this bounded region thereafter.

The theoretical bound $b_d$ reveals a critical engineering design trade-off. While mathematically minimizing this bound requires maximizing the convergence speed via a very large observer gain, doing so practically risks amplifying unmodeled high-frequency noise and inducing undesirable transient behaviors. Prioritizing system stability and smooth transitions dictates the selection of a more moderate gain. The specific numerical selection of this gain, alongside a comprehensive analysis of its impact on the transient behavior, is detailed in the simulation results presented in Section 7.

6. Stability Analysis of the Interconnected NDI-DOBC System

To guarantee the robustness of the complete attitude control architecture, it is imperative to analyze the coupled stability of the linear quadratic tracking controller and the nonlinear disturbance observer. While the previous sections established the stability of each subsystem in isolation, the estimated disturbance $\widehat{\mathit{d}}$ introduces a cross-coupling effect into the closed-loop tracking dynamics must be bounded.

6.1. Perturbed Tracking Error Dynamics

Let the state tracking error of the linearized system be defined as $\mathit{e}\left(t\right)=\mathit{\xi }\left(t\right)-{\mathit{\xi }}_{ref}\left(t\right)$. Without environment disturbances, the closed-loop LTI dynamics under the state-feedback LQ controller are governed by $\dot{\mathit{e}}={\mathit{A}}_c\mathit{e}$, where ${\mathit{A}}_c=(\mathit{A}-\mathit{B}\mathit{K})$ is a Hurwitz system matrix.

However, under the influence of the physical disturbance $\mathit{d}$ and the composite control law $\mathit{u}={\mathsf{\Gamma }}_B-\widehat{\mathit{d}}$, the input transformation is purturbed by the disturbance estimation error, ${\mathit{e}}_{\mathit{d}}=\mathit{d}=\widehat{\mathit{d}}$. Consequently, the closed-loop tracking error dynamics are formulated as:

$$\dot{\mathit{e}}={\mathit{A}}_c\mathit{e}+\mathit{B}\mathit{G}\left(\mathit{x}\right){\mathit{e}}_{\mathit{d}}$$

(56)

Where $B={\left[\begin{array}{cc}{0}_{4\times 4}& {\mathit{I}}_{4\times 4}\end{array}\right]}^{\top }$ is the linear input matrix mapping the virtual acceleration.

6.2. Composite Lyapunov Function and Stability Proof

To evaluate the interconnected system, we construct a composite Lyapunov candidate function. A weighting constant $\gamma >0$ is introduced to scale the observer’s estimation error dynamics:

$$V_{total}\left(\mathit{e},{\mathit{e}}_{\mathit{d}}\right)={\mathit{e}}^{\mathit{\top }}\mathit{P}\mathit{e}+\gamma \left({\displaystyle \frac{1}{2}}{\mathit{e}}_{\mathit{d}}^{\mathit{\top }}{\mathit{e}}_{\mathit{d}}\right)$$

(57)

For the LTI tracking system, because ${\mathit{A}}_c$ is Hurwitz, for any given symmetric positive-definite matrix $\mathit{Q}$, there exists a unique symmetric positive-definite matrix $\mathit{P}$ that satisfies the continuous-time Lyapunov equation:

$${\mathit{A}}_c^{\mathit{\top }}\mathit{P}+\mathit{P}{\mathit{A}}_c=-\mathit{Q}$$

(58)

Taking the time derivative of $V_{total}$ along the system trajectories yield:

$${\dot{V}}_{total}=\left({\dot{\mathit{e}}}^{\top }\mathit{P}\mathit{e}+{\mathit{e}}^{\mathit{\top }}\mathit{P}\dot{\mathit{e}}\right)+\gamma {\mathit{e}}_{\mathit{d}}^{\mathit{\top }}{\dot{\mathit{e}}}_{\mathit{d}}$$

(59)

Substituting Equation (56) and the observer error dynamics (${\dot{\mathit{e}}}_{\mathit{d}}=-\left(1/2\right)K_{obs}{\mathit{e}}_{\mathit{d}}+\dot{\mathit{d}}$) results in:

$${\dot{V}}_{total}={\mathit{e}}^{\mathit{\top }}\left({\mathit{A}}_c^{\mathit{\top }}\mathit{P}+\mathit{P}{\mathit{A}}_c\right)\mathit{e}+2{\mathit{e}}^{\mathit{\top }}\mathit{P}\mathit{B}\mathit{G}\left(\mathit{x}\right){\mathit{e}}_{\mathit{d}}-{\displaystyle \frac{\gamma }{2}}{K}_{obs}||{\mathit{e}}_{\mathit{d}}|{|}^2+\gamma {\mathit{e}}_{\mathit{d}}^{\mathit{\top }}\dot{\mathit{d}}$$

(60)

Applying Equation (58) and utilizing the Rayleigh quotient, the first term on the right-hand side, the nominal tracking term, is bounded by the minimum eigenvalue of $\mathit{Q}$:

$$-{\mathit{e}}^{\mathit{\top }}\mathit{Q}\mathit{e}\le -{\lambda }_{min}\left(\mathit{Q}\right)||\mathit{e}|{|}^2$$

(61)

To separate the cross-coupling term $2{\mathit{e}}^{\mathit{\top }}\mathit{P}\mathit{B}\mathit{G}\left(\mathit{x}\right){\mathit{e}}_{\mathit{d}}$, we apply Young’s inequality. For any vectors $\mathit{a},\mathit{b}$ and any strict positive scalar $ϵ>0$, the inequality states $2{\mathit{a}}^{\mathit{\top }}\mathit{b}\le \left(1/ϵ\right)||\mathit{a}|{|}^2+ϵ||\mathit{b}|{|}^2$.

Let $\mathit{a}={\mathit{G}}^{\mathit{\top }}\left(\mathit{x}\right){\mathit{B}}^{\mathit{\top }}\mathit{P}\mathit{e}$ and $\mathit{b}={\mathit{e}}_{\mathit{d}}$:

$$\begin{array}{cc}\hfill 2{\mathit{e}}^{\top }\mathit{P}\mathit{B}\mathit{G}\left(\mathit{x}\right){\mathit{e}}_{\mathit{d}}& \le {\displaystyle \frac{1}{ϵ}}||{\mathit{G}}^{\mathit{\top }}\left(\mathit{x}\right){\mathit{B}}^{\mathit{\top }}\mathit{P}\mathit{e}|{|}^2+ϵ||{\mathit{e}}_{\mathit{d}}{|}^2\hfill \\ & \le {\displaystyle \frac{1}{ϵ}}||\mathit{P}\mathit{B}\mathit{G}\left(\mathit{x}\right)|{|}^2||\mathit{e}|{|}^2+ϵ||{\mathit{e}}_{\mathit{d}}|{|}^2\hfill \end{array}$$

(62)

Substituting these bounds into Equation (60) and applying the physical disturbance bound $||\dot{\mathit{d}}||\le \delta$, we obtain:

$${\dot{V}}_{total}\le -{\lambda }_{min}\left(\mathit{Q}\right)||\mathit{e}|{|}^2+{\displaystyle \frac{1}{ϵ}}||\mathit{P}\mathit{B}\mathit{G}\left(\mathit{x}\right)|{|}^2||\mathit{e}|{|}^2+ϵ||{\mathit{e}}_{\mathit{d}}|{|}^2-{\displaystyle \frac{\gamma }{2}}{K}_{obs}||{\mathit{e}}_{\mathit{d}}|{|}^2+\gamma \delta ||{\mathit{e}}_{\mathit{d}}||$$

(63)

Rearranging the terms to group the state norms:

$${\dot{V}}_{total}\le -\left({\lambda }_{min}\left(\mathit{Q}\right)-{\displaystyle \frac{1}{ϵ}}||\mathit{P}\mathit{B}\mathit{G}\left(\mathit{x}\right)|{|}^2\right)||\mathit{e}|{|}^2-\left({\displaystyle \frac{\gamma }{2}}{K}_{obs}-ϵ\right)||{\mathit{e}}_{\mathit{d}}|{|}^2+\gamma \delta ||{\mathit{e}}_{\mathit{d}}||$$

(64)

For the interconnected system to maintain bounded stability, the coefficients of the quadratic error terms must be strictly positive. This establishes two critical design conditions that govern the selection of the Young’s inequality scaling parameter $ϵ$, the Lyapunov weighting constant $\gamma$, and the observer gain $K_{obs}$:

$$\begin{array}{cc}1.& {\lambda }_{min}\left(\mathit{Q}\right)-{\displaystyle \frac{1}{ϵ}}||\mathit{P}\mathit{B}\mathit{G}\left(\mathit{x}\right)|{|}^2>0\Rightarrow ϵ>{\displaystyle \frac{||\mathit{P}\mathit{B}\mathit{G}\left(\mathit{x}\right)|{|}^2}{{\lambda }_{min}\left(\mathit{Q}\right)}}\\ 2.& {\displaystyle \frac{\gamma }{2}}{K}_{obs}-ϵ>0\Rightarrow \gamma >{\displaystyle \frac{2ϵ}{K_{obs}}}\end{array}$$

(65)

By first selecting a sufficiently large $ϵ$ to satisfy the first condition, we can subsequently define a weighting constant $\gamma$ and an observer gain $K_{obs}$ that satisfy the second condition. With these parameters established, let the strictly positive constants be defined as $C_1={\lambda }_{min}\left(\mathit{Q}\right)-{\displaystyle \frac{1}{ϵ}}||\mathit{P}\mathit{B}\mathit{G}\left(\mathit{x}\right)|{|}^2$ and $C_2={\displaystyle \frac{\gamma }{2}}{K}_{obs}-ϵ$.

The Lyapunov derivative simplifies to:

$${\dot{V}}_{total}\le -C_1||\mathit{e}|{|}^2-C_2||{\mathit{e}}_{\mathit{d}}|{|}^2+\gamma \delta ||{\mathit{e}}_{\mathit{d}}||$$

(66)

Rearranging Equation (66) yields:

$${\dot{V}}_{total}\le -C_1||\mathit{e}|{|}^2-||{\mathit{e}}_{\mathit{d}}||\left(C_2||{\mathit{e}}_{\mathit{d}}||-\gamma \delta \right)$$

(67)

For the overall Lyapunov derivative to be strictly negative (${\dot{V}}_{total}<0$), it is sufficient that the norm of the disturbance estimation error satisfies:

$$||{\mathit{e}}_{\mathit{d}}||>{\displaystyle \frac{\gamma \delta }{C_2}}$$

(68)

Since $-C_1||\mathit{e}|{|}^2$ is always semi-negative, whenever this condition holds, ${\dot{V}}_{total}$ is negative definite outside a bounded compact set.

This mathematical inequality rigorously guarantees that the complete interconnected NDI-DOBC system is Uniformly Ultimately Bounded (UUB). The state tracking error $||\mathit{e}||$ and the disturbance estimation error $||{\mathit{e}}_{\mathit{d}}||$ will concurrently converge to a compact region bounded by the physical limits of the testbed’s disturbance derivative ($\delta$). Once the errors converge into this boundary, they will remain there indefinitely, formally proving the robust tracking stability of the proposed dual-loop architecture.

7. Simulation Results

This section presents the simulation results utilized to validate the efficacy of the proposed NDI-DOBC methodology. We categorize the simulations into four cases. First, the attitude tracking control with the zero-crossing of the scalar component of quaternion ($q_0=0$) is examined under ideal conditions, characterized by a diagonal moment of inertia matrix and the absence of external disturbance. To underscore the singular-free attributes of the proposed NDI framework, we adopt the moment of inertia matrix and initial, final conditions from previous work [8] and compare the simulation results. In the subsequent case, we adjust the moment of inertia matrix and add off-diagonal elements to test the robustness of the linear quadratic tracking controller in cooperation with the NDI framework.

Before detailing the specific simulation cases, the physical parameters of the attitude control testbed – which define the theoretical disturbance bounds established in Section 2.2 – are quantified. The modeled testbed has a total mass of 19.899 kg and a center of gravity (CG) displacement of [0.02, 0.02, 0.05] m from the center of the spherical air-bearing. Evaluated under standard local gravitational acceleration, this yields an absolute maximum gravitational imbalance torque bound of $D_{max}=3.83N\cdot m$. Furthermore, the physical control torque saturation limit of the actuator is strictly set to 5.25 $N\cdot m$. Constrained by this maximum control effort and the resulting rotational dynamics, the maximum allowable rate of change of the disturbance torque is calculated and bounded by $\delta =11.21N\cdot m/s$.

The first case focuses on the general spacecraft static-to-static attitude control scenario, where the initial and target attitudes are denoted by ${\mathit{q}}_I=[-0.742 0.2 0.4 0.5{]}^{\top }$ and ${\mathit{q}}_F=[0.447 0.6 -0.2 -0.632{]}^{\top }$, respectively. To further demonstrate the performance optimization capabilities of the proposed architecture, Case 1 also includes an additional comparative result utilizing a smaller system time constant. This specific tuning is designed to significantly accelerate the convergence rate while remaining strictly within the physical control saturation limits, highlighting an optimal balance between aggressive tracking and actuator constraints.

In contrast, the second case considers a dynamic transition from an initially rotating state with an angular velocity of $\mathit{\omega }={\left[\begin{array}{ccc}0.5& 0.5& 0.5\end{array}\right]}^{T}rad/s$ to a static target attitude ${\mathit{q}}_{\mathit{F}}$. The moment of inertia matrices for simulation Case 1 and 2 are defined as follows. In Case 1, a purely diagonal moment of inertia matrix, ${\widehat{\mathit{J}}}_d$, is used for both the derivation of the NDI control law and the plant simulation. In case 2, ${\widehat{\mathit{J}}}_d$ is retained for the NDI controller derivation to introduce a deliberate model mismatch, while the actual plant simulation is conducted using the full matrices featuring off-diagonal coupled inertias (${\widehat{\mathit{J}}}_1$, ${\widehat{\mathit{J}}}_2$ and ${\widehat{\mathit{J}}}_3$):

$$\begin{array}{ccc}{\widehat{\mathit{J}}}_{\mathit{d}}=\left[\begin{array}{ccc}300& 0& 0\\ 0& 320& 0\\ 0& 0& 250\end{array}\right]& {\widehat{\mathit{J}}}_1=\left[\begin{array}{ccc}300& -32& 25\\ -32& 320& -30\\ 25& -30& 250\end{array}\right]& \\ & & (\mathrm{kg}\cdot {\mathrm{m}}^2)\\ {\widehat{\mathit{J}}}_2=\left[\begin{array}{ccc}150& -32& 25\\ -32& 160& -30\\ 25& -30& 125\end{array}\right]& {\widehat{\mathit{J}}}_3=\left[\begin{array}{ccc}90& -32& 25\\ -32& 96& -30\\ 25& -30& 75\end{array}\right]& \end{array}$$

(69)

The third simulation case introduces external disturbances, focusing on the attitude control testbed as the physical target. As illustrated in Figure 2, the disturbance is induced by the misalignment between the testbed’s center of gravity and the support center of the air bearing, generating a state-dependent gravitational moment. The proposed DOBC is utilized to estimate and dynamically counteract this specific disturbance impact. Given the constrained roll and pitch angles of the physical testbed, new operational initial and target conditions are defined. Represented in Euler angles, the testbed transitions from an initial roll, pitch, and yaw of -5°, 10°, and 0° to a final target of 0°, 0°, and 60°.

The fourth case considers the simultaneous presence of both model mismatch and external disturbance, reflecting the realistic scenario of the attitude control testbed. The true moment of inertia of the testbed, ${\widehat{\mathit{J}}}_{true}$, is calculated using CAD software and applied as the plant baseline in simulation Cases 3 and 4. To test the framework’s robustness, the NDI control law is deliberately handicapped by only being supplied with the diagonal elements of ${\widehat{\mathit{J}}}_{true}$.

$${\widehat{\mathit{J}}}_{true}=\left[\begin{array}{ccc}1.1932& -0.0474& 0\\ -0.0474& 1.116& 0\\ 0& 0& 1.0798\end{array}\right] (\mathrm{kg}\cdot {\mathrm{m}}^2)$$

(70)

7.1. Case 1: Zero Crossing of Quaternions

The simulation of quaternions is shown in Figure 3. The results are compared with those from Bang’s work [8], which utilizes a conventional linear controller alongside Newtonian NDI. The scalar component of quaternion, $q_0$, crosses zero at approximately t=12.7 seconds in our simulation, indicated by the green dashed line; In Bang’s case, this crossing occurs at t=7.9 seconds, marked by the blue dashed line. We can see that Bang’s controller experiences significant oscillations before reaching the desired attitude, whereas the proposed method yields a remarkably smooth convergence.

The most significant difference between the two controllers lies in the derivatives of the control input, as shown in Figure 4. Because Bang’s controller relies on an input-output linearization that possesses a mathematical singularity at $q_0=0$, the system experiences severe numerical instability. This singularity causes massive, rapid fluctuations in the control input, with the derivative exceeding 500 $N\cdot m/s$, even with their proposed modified gain intended to mitigate the singularity. In contrast, the change in the control input for the proposed Lagrangian NDI remains completely smooth and bounded throughout the exact moment $q_0$ crosses zero.

To further demonstrate the tuning flexibility of the proposed LQ controller and the robustness of the singularity-free NDI framework, an additional simulation is conducted with a reduced system time constant ($\tau =3.5s$). The quaternion trajectories and the corresponding control input derivative for the aggressive tuning are presented in Figure 5 and Figure 6, respectively. As depicted in Figure 5, decreasing the time constant significantly accelerates the system’s convergence rate, allowing the spacecraft to reach the target attitude much faster than the nominal case.

However, as confirmed by the control input derivatives in Figure 6, the proposed architecture successfully accommodates the accelerated tracking without numerical failure. It is important to clarify the distinct roles of the control architecture components: while actuator saturation is prevented by the integrated first-order filter, the fundamental contribution of the proposed NDI is the complete elimination of the mid-maneuver control derivative peak. As seen in Figure 6, the control input derivative remains entirely smooth and continuous as $q_0$ crosses zero, proving that the Lagrangian formulation is fundamentally free from mathematical singularities regardless of the desired convergence speed.

Furthermore, an initial peak in the control input derivative – reaching approximately 400 $N\cdot m/s$ – is observable at $t=0$ in Figure 6. This is a standard start-up transient caused by the instantaneous step in the commanded control torque from a zero initial state. Because the reduced time constant dictates a more aggressive initial angular acceleration, this start-up derivative is naturally larger than the case presented in Figure 4. Crucially, this initial transient is an expected physical response, completely distinct from the mathematically induced singularity spike that plagues conventional methods during the transient tracking phase.

7.2. Case 2: Robustness to Mismatched Moment of Inertia

In this case, we use the moment of inertia matrices ${\widehat{\mathit{J}}}_1$, ${\widehat{\mathit{J}}}_2$, and ${\widehat{\mathit{J}}}_3$ shown in Equation (69) for numerical simulation, but still adopt ${\widehat{\mathit{J}}}_d$ for deriving the NDI formulation. The cross-coupling terms and the decreasing diagonal terms in the moment of inertia matrices simulate the misalignment between principal axes and body-fixed axes and decreasing mass. The spacecraft is controlled by the LQ controller defined in Equation (34) with the weighting matrices $\mathit{Q}=10\mathit{I}$ and $\mathit{R}=\mathit{I}$.

According to Newtonian rotational dynamics, the moment of inertia predominantly governs a system’s transient rotational behavior. Consequently, to explicitly evaluate this effect, Case 2 – unlike Case 1 – is initialized with a non-zero angular velocity. As illustrated in Figure 7, although variations in the moment of inertia yield distinct transient paths, all trajectories successfully converge to the reference quaternions. These results effectively demonstrate the robustness of the combined NDI and LQ control scheme in handling model mismatches.

7.3. Case 3: Robustness to External Disturbances

We introduce the DOBC derived in Section 5 to compensate for the influence of disturbances. The disturbance observer gain, $K_{obs}$, is tested with values of 1, 5, 10, and 20. Figure 8 shows the disturbance estimation error for the observer with different gains. Since the disturbances are related to the testbed’s attitude, we use the error between the disturbance and the estimated one as the performance index for the observer. Figure 8 indicates that the larger the observer gain, the faster the estimation error converges.

Figure 9 and Figure 10 also show that the observer gain will significantly influence the performance of attitude control. It appears that all the estimation errors converge to zero in Figure 8, but biases still exist at the steady state. Take y-axis for example, the value of $K_{obs}$ equals 1 (green line) and 20 (red line) are in the order of magnitude of ${10}^{-2}$ and ${10}^{-4}$. These biases cause the poor performance of the green line in Figure 9 and Figure 10.

Although the red lines ($K_{obs}=20$) have the fastest convergence time, this large gain not only induces oscillations in disturbance estimation but also affects attitude control. Therefore, we select $K_{obs}=10$ as the observer gain for implementation.

7.4. Case 4: Attitude Control with Both Model Mismatch and Disturbance

Finally, we adopt controller and observer gains designed from the two cases above and implement them in the scenario that contains both model mismatch and disturbance. The result is shown in Figure 11, where the yaw angle smoothly aligns with the desired heading. Meanwhile, the pitch and roll axes struggle to compensate for the disturbances caused by the initial pitch and roll angles. However, the DOBC rapidly takes effect to suppress these perturbations. Consequently, both axes achieve a smooth convergence to zero, while the yaw axis perfectly tracks the 60° target.

8. Conclusion

This paper presents a robust, singularity-free attitude tracking architecture for spacecraft, integrating a novel quaternion-based nonlinear dynamic inversion (NDI) framework with a disturbance observer-based controller (DOBC). By deriving the rotational dynamics strictly through Udwadia’s Lagrangian formulation, the proposed NDI achieves exact input-state linearization directly on the active holonomic constraint manifold. The geometric realization of this state transformation is mathematically proven to eliminate internal zero dynamics and entirely avoid the severe control derivative discontinuities historically associated with the $q_0=0$ singularity. To ensure robust performance against physical model uncertainties and time-varying environmental disturbances, a nonlinear DOBC is synthesized to augment a baseline linear quadratic (LQ) tracking controller. Rigorous Lyapunov stability analysis formally guarantee that both the isolated disturbance estimation error and fully interconnected dual-loop NDI-DOBC system are Uniformly Ultimately Bounded (UUB), even in the presence of non-vanishing disturbance derivatives ($\dot{\mathit{d}}\ne 0$). Comprehensive numerical simulations, parameterized by a physical spherical air-bearing testbed, validate the architecture’s exceptional precision, smooth transient response, and robust disturbance rejection without inducing actuator saturation. Because the proposed Lagrangian NDI seamlessly maps the complex nonlinear rotational dynamics into a globally valid synthetic linear time-invariant domain, future research will focus on leveraging this exact linearization to deploy advanced optimal control strategies, such as Model Predictive Control (MPC), and exploiting the linear state constraints to significantly accelerate computational convergence for real-time trajectory optimization problems.