Determining the torque needed to rotate a spacecraft to a given rotation/quaternion

Which is the correct algorithm to determine the torque needed to rotate a spacecraft to a given quaternion ? I have a set of quaternions of a spacecraft and the time difference (delta T) between the quaternions.

I've tried estimating the angular velocity and angular acceleration as: 1. Find quaternion q so q = q1/q0

1. Convert quaternion q to axis and angle:

len=sqrt(q.xq.x+q.yq.y+q.z*q.z); angle=2*atan2(len, q.w); axis = q.xyz()/len;

1. angular velocity w = axis * angle / dt.
2. angular acceleration a = w/dt.

Using this algorithm, the computed angular acceleration (a) is not correct.

What I'm doing wrong? Is this the correct way to obtain the angular acceleration and torque from a sequence of two quaternions?

• What are the initial and final angular velocities of the spaceship? Your algorithm doesn't mention torque at all. Are you really trying to compute torque. – Brian Moths Mar 9 '17 at 22:13
• @NowIGetToLearnWhatAHeadIs The inital and final angular velocities are known ( but I don't know the values at this moment ). I didn't specify the torque determination because that part I got it "covered". My issue is finding the angular acceleration between those two quaternions. – Madalin Mar 9 '17 at 22:18
• Well doesn't the angular acceleration depend only on the initial and final angular velocities and and not the orientations? What you are asking right now is like saying something travelled three meters in five seconds, what was the acceleration? You have no idea because acceleration is related to a difference in velocity not position. – Brian Moths Mar 9 '17 at 23:20
• @NowIGetToLearnWhatAHeadIs Actually I'm stucked because I don't have strong knowledge about this. Knowing initial and final angular velocity and the time between these is enough to find the average acceleration ? I don't need to know how the satellite rotated from initial to final state ? – Madalin Mar 9 '17 at 23:29
• @Madalin You need both initial and final orientation and angular velocity, but also how you plan to change angular velocity. For example you could apply a "torque impulses" at start and finish, which first cancels the current angular velocity and adds the average angular velocity, and then cancels this average angular velocity and adds the final angular velocity. But you can also use quaternion interpolation and the corresponding torque to achieve that behavior. – fibonatic Mar 10 '17 at 0:05

To find the average angular acceleration $\boldsymbol{\alpha}$ over the time period $\Delta t$, use the formula $$\boldsymbol{\alpha}=\frac{\boldsymbol{\omega}_f-\boldsymbol{\omega}_i}{\Delta t},$$

where $\boldsymbol{\omega}_i$ is the angular velocity at the beginning of time period and $\boldsymbol{\omega}_f$ is the angular velocity at the end.

Now you say you already know how to compute the torque, but for completeness, I will give the formula for this is well. The average torque $\boldsymbol{\tau}$ is given by

$$\boldsymbol{\tau}=\frac{\mathbf{I}_f\boldsymbol{\omega}_f-\mathbf{I}_i\boldsymbol{\omega}_i}{\Delta t},$$

where $\mathbf{I}_i$ is the moment of inertia tensor at the beginning and $\mathbf{I}_f$ is the moment of inertia tensor at the end.

• Thank you for your explanation. I've calculated the torque using the formulas you've written. Now, after making a propagation with the calculated torque and extract the quaternion from the final state of the propagation (Q'), between the reference quaternion (Q) and the obtained quaternion (Q') I have the following difference in the Euler angles: x:0.0146 deg y:0.4235 deg z:0.0146 deg. Is there any algorithm that I can apply to have a smaller difference between Q and Q' ? – Madalin Mar 10 '17 at 14:58
• For example, I've tried using the angular velocity of Q' state as the "new initial state" and calculated Q'' (same method as for Q') but the differences are not decreasing, they remain constant. – Madalin Mar 10 '17 at 15:01
• @Madalin Your question (as far as I could tell) "Is this the correct way to obtain the angular acceleration and torque from a sequence of two quaternions?" I think I have answered that question. I don't understand how your comments address my answer. Could you make it more clear what you don't understand about my answer. Of course, if you have another, different question, you can feel free to ask that one as well on this site. – Brian Moths Mar 10 '17 at 15:12
• @Madalin You might already know this, but I think the question you are trying to ask is "Given an object whose orientation is given by a quaternion $q$, and given that at some time $t$, $q(t)$ and $\dot{q}(t)$ is known, what constant torque must be applied over a time $\Delta t$ so that $q(t+\Delta t)$ attains some desired value." – Brian Moths Mar 10 '17 at 15:29

From your timeseries of the rotation quaternion $\gamma$, estimate the following derivatives:

$$\Omega = \gamma^{-1} \,\dot{\gamma}$$ $$\dot{\Omega} = \gamma^{-1}\, \ddot{\gamma}-\Omega^2$$

with all the wonted provisos and warnings about calculating numerical derivatives from timeseries. The above two quaternions will be purely imaginary quaternions because $\gamma$ is always a unit quaternion. Then, to calculate the torque, you need to know the inertia tensor of your spacecraft and you use my formula (6) below, at the end of the following analysis.

To work with quaternions ($SU(2)$ representation of the rotation group, with $SU(2)\stackrel{\mathrm{Ad}}{\rightarrow} SO(3)$), you do the following.

As you likely know, a point in 3-space are represented by the purely imaginary quaterion $X=x\,\mathbf{i}+y\,\mathbf{j}+z\,\mathbf{k}$, and the rotation represented by quaternion $\gamma$ acts on $X$ through

$$X\mapsto \gamma\,X\,\gamma^{-1}\tag{1}$$

Differentiate this and you'll find that:

$$\dot{X} = \gamma\,[\Omega,X]\,\gamma^{-1} \tag{2}$$

where $\Omega = \gamma^{-1}\,\dot{\gamma}$ is also a pure quaternion and represents the angular velocity. Indeed, if you write $\Omega = \omega_x\,\mathbf{i}+\omega_y\,\mathbf{j}+\omega_z\,\mathbf{k}$, then the operation represented by the Lie bracket $[\Omega,X] = \Omega\,X-X\,\Omega$ (as matrices) is indeed the cross product $\Omega\times X$ when $\Omega$ and $X$ are thought of as 3-vectors. The moment of the velocity in (2) is:

$$L = \gamma\,[X,\,[\Omega,\,X]]\,\gamma^{-1}\tag{3}$$

At this point, it becomes a little hard to stay in quaternion notation completely because we need to integrate (3) over the whole body to get the inertia tensor. We need to calculate $\int_{body} \rho(X)\,[X,\,[X,\,\Omega]]\mathrm{d} V$; this is simply a $3\times3$ matrix operator that acts on the three components of $\Omega$ and given by:

$$I = \int_{body} \rho(X)\,\mathrm{ad}(X)^2 \mathrm{d} V\tag{4}$$

where $\mathrm{ad}(X)$ is the $3\times 3$ matrix of the linear operation $Y\mapsto [X,\,Y]$ for a pure quaternion $Y$.

Let's write this linear mapping on $\Omega$ by the inertia tensor $I(\Omega)$.

To calculate the torque, we need to calculate the time derivative of (3); repeating the trick we used to get (2), we find:

$$\tau = \gamma\,\left([\Omega,\,I(\Omega)]+\mathrm{d}_t I(\Omega)\right)\,\gamma^{-1}\tag{5}$$

which, when we rotate with the frame attached to the rigid body, becomes (on transforming $\tau\mapsto \gamma^{-1}\,\tau\,\gamma$):

$$\tau = [\Omega,\,I(\Omega)]+ I(\dot{\Omega})\tag{6}$$

which is the equivalent of obtaining the Euler equations through the formula $\mathrm{d}_t = D_t + \Omega\times$.

• Thank you very much for the effort you've done explaining me. I really appreciate it ! I've tried the formula presented below by @NowIGetToLearnWhatAHeadIs and after making a propagation with the calculated torque and extract the quaternion from the final state of the propagation (Q'), between the reference quaternion (Q) and the obtained quaternion (Q') I have the following difference in the Euler angles: x:0.0146 deg y:0.4235 deg z:0.0146 deg. Now I'm trying to implement an iterative algorithm that decreases these errors. – Madalin Mar 10 '17 at 15:07