The goal of this controller was to have an easy to adjust damping ratio set the performance and response characteristics.
You can first set a settling time target and a damping ratio target, which can then be used to compute \(K_p\) and \(K_d\) gains for the system.
1) Start from rigid-body rotational dynamics
Euler’s rotational equation gives the relationship between torque and angular acceleration.
The plant uses spacecraft inertia \(J\), body rate \(\omega\), and total applied torque \(\tau\).
\[
J\dot{\omega} + \omega \times (J\omega) = \tau
\]
For a reaction-wheel controlled spacecraft, the commanded wheel torque drives the body torque with opposite sign.
The controller output is \(u\) and the plant torque is \(\tau = -u\).
\[
J\dot{\omega} + \omega \times (J\omega) = -u
\]
2) Define the quaternion convention and attitude error
The state uses a scalar-first quaternion \(q = [q_0,\ q_v^T]^T\), with \(q_v \in \mathbb{R}^3\).
Quaternion kinematics relate rate to quaternion derivative.
\[
q =
\begin{bmatrix}
q_0 \\
q_v
\end{bmatrix},
\quad
\dot{q}_0 = -\frac{1}{2}q_v^T\omega,
\quad
\dot{q}_v = \frac{1}{2}\left(q_0 I_3 + [q_v]_\times\right)\omega
\]
The attitude error quaternion uses the standard composition rule.
\[
q_e = q_d^{-1} \otimes q =
\begin{bmatrix}
q_{e0}\\
q_{ev}
\end{bmatrix}
\]
As quaternions have a double cover, \(q\) and \(-q\) represent the same attitude.
To enforce the short rotation, I flip the error quaternion when the scalar term is negative. Otherwise, instead of commanding a 45° rotation, the system could command a -315° rotation, which would be significantly less efficient.
\[
q_{e0} < 0 \Rightarrow q_e \leftarrow -q_e
\]
3) Convert quaternion error to a vector error using the log map
A PD law needs an error vector in \(\mathbb{R}^3\).
The log map converts the error quaternion into an axis-angle rotation and then into a 3D error vector.
\[
s = q_{e0}, \quad v = q_{ev}, \quad r = \lVert v \rVert
\]
\[
\theta = 2\ \mathrm{atan2}(r, s)
\]
\[
\hat{a} =
\begin{cases}
v/r, & r > \epsilon \\
\text{any unit vector}, & r \le \epsilon
\end{cases}
\]
\[
e = \theta \hat{a}
\]
This avoids the small-angle issues you hit with raw quaternion vector parts when errors get larger and is able to produces an error direction and magnitude in one vector.
4) Apply the PD control law and linearize for tuning
We define the rate error \(\tilde{\omega} = \omega - \omega_d\). For detumble and hold, typically \(\omega_d = 0\), although this can be set within the scenario command.
For a general PD control law:
\[
u = K_p e + K_d \tilde{\omega}
\]
For gain selection, I chose to use a small-angle approximation and ignore the cross-coupling term to reduce complexity and be able to form a spring-dampener system.
\[
J\dot{\omega} \approx -u
\]
With small attitude errors, the error kinematics reduce to \(\dot{e} \approx -\omega\) for a zero command rate.
We can then substitute \(\omega \approx -\dot{e}\) and \(\dot{\omega} \approx -\ddot{e}\) into the dynamics.
\[
J(-\ddot{e}) = -K_p e - K_d(-\dot{e})
\]
\[
J\ddot{e} + K_d \dot{e} + K_p e = 0
\]
5) Match to a standard second-order system
If \(J\), \(K_p\), and \(K_d\) are diagonal (which they are in this system), each axis behaves like an independent second-order system.
We can then compare the closed-loop form to the standard template.
\[
\ddot{e}_i + 2\zeta_i \omega_{n,i}\dot{e}_i + \omega_{n,i}^2 e_i = 0
\]
\[
\omega_{n,i}^2 = \frac{K_{p,i}}{J_i},
\quad
2\zeta_i\omega_{n,i} = \frac{K_{d,i}}{J_i}
\]
With these terms broken out, we can solve for the response gains in terms of inertia, natural frequency, and damping ratio.
\[
K_{p,i} = J_i \omega_{n,i}^2,
\quad
K_{d,i} = 2\zeta_i J_i \omega_{n,i}
\]
6) Use the 2 percent settling-time approximation
A common design rule for a second-order system is the 2 percent settling-time approximation:
\[
T_{s,i} \approx \frac{4}{\zeta_i \omega_{n,i}}
\]
If we substitute \(\omega_{n,i} = 4/(\zeta_i T_{s,i})\) into the gain expressions, this produces our direct tuning knobs: pick \(T_s\) and \(\zeta\), we can compute \(K_p\) and \(K_d\).
\[
K_{p,i} = J_i\left(\frac{4}{\zeta_i T_{s,i}}\right)^2,
\quad
K_{d,i} = \frac{8J_i}{T_{s,i}}
\]
In matrix form with diagonal \(J\), this matches the MATLAB initialization:
\[
K_p = J\left(\frac{4}{\zeta T_s}\right)^2,
\quad
K_d = J\left(\frac{8}{T_s}\right)
\]
Benefits
With this relationship established, you can get predictable tuning behavior. Smaller \(T_s\) increases \(K_p\) and \(K_d\). Lower \(\zeta\) increases overshoot.
This equation setup allows for easy response adjustments as can be seen in the Demo 2 visualization.
