simtk.org

Hi all,

My model comprises two bodies (A parent, B child) joined together by a custom joint with 1 rotational DOF and I would like to express the angular velocity in B's frame. After setting the states and using the computeStateDerivatives function I then can get the inertial Angular velocity between the two bodies as follows:

w_velG = Vec3(0);
simbodyEngine.getAngularVelocity(osimState,bodyB,w_velG);

My intuition then tells me that I need to determine the rotation matrix between the two bodies:

tform = simbodyEngine.getTransform(osimState,bodyB);
rotMat_B = tform.R();

Afterwards I should then multiply the transpose of the rotation matrix with the inertial angular velocity to determine the angular velocity in the local frame of B:

w_velB = rotMat_B' * w_velG

I am, however, very uncertain as to whether this is correct, as when I perform:

rotMat_B * w_velG

the vector components are identical and this leads me to thinking I am making a mistake. Could anyone please clarify whether my approach is correct?

Thanks,

Nicos

Nicos,

What you are doing is fundamentally correct, but it is indeed odd that you get the same result, whether you transpose the rotation matrix or not.

In general, this can only happen if the rotation matrix is the identity matrix.

But in this special case, the angular velocity vector is oriented along the axis of rotation of the joint. This axis may be a coordinate axis that is common to both bodies. In that case, the global (inertial) and local coordinates of the angular velocity vector are the same.

Does that make sense?

Inspecting the rotation matrix may help understand this.

Ton van den Bogert

Hi Nicos,

Just a suggestion. You can use the new OpenSim API to transform quantities into different frames easily, without having to specify the transformations explicitly:

https://simtk.org/api_docs/opensim/api_ ... Frame.html

Each body is considered an OpenSim::Frame.

In your case, it might be that the transformation is not correct. Is the state properly realized? As Ton mentioned you can check the values of the rotation matrix.

Thank you both for your prompt replies.

Ton,

If I set the DOF to be 0 zero then the rotation matrix returned is the identity matrix as expected. If I set the the DOF to be pi/3 and the generalised speed to be 1 rad/s, for example, the rotation matrix and angular velocity returned are:
0.5017 -0.8646 -0.0274
0.8645 0.5000 0.0522
-0.0314 -0.0499 0.9983

~[-0.058898,0.0023,0.998261]

I would have expected it to look something more like the following for the rotation matrix:
0.5000 -0.8660 0
0.8660 0.5000 0
0 0 1.0000

But I believe this is because the rotation takes place about an axis that is defined by the following vector <axis>-0.05889802 0.0023 0.99826136</axis>.

If the rotation had taken place about <axis>0 0 1</axis> your point regarding the 'special case' makes sense to me, global angular velocity == local angular velocity. But I still cannot seem to grasp why the multiplication of the transposed and non-transposed rotation matrix with the angular velocity vector return the same value.

Dimitar,

I haven't transitioned yet over to the new API for Matlab - perhaps it is something I should try. I am still using 3.3 as I currently have all my code compatible with that version. I am confident I have realised the State correctly as when I use calculateStateDerivatives the first returned values correspond to the generalised velocities I input to the State - which I would expect to see.

Thanks,

Nicos

The rotation matrix seems fine. If you transpose it, then it will yield a different matrix. Maybe the MATLAB transpose operator ' is not interpreted correctly, because the matrix is of type SimTK::Matrix and not MATLAB array. Try using the conventional function transpose() instead:

w_velB = rotMat_B.transpose() * w_velG

Nicos,

Thanks for posting the numerical values.

This confirms my suspicion. In your system there is only one rotation, so the angular velocity vector of the moving body has the same direction as the axis of rotation.

The angular velocity vector is an eigenvector of the rotation matrix (and of the inverse/transpose rotation matrix),
and the corresponding eigenvalue is 1. So the vector components are not affected by the rotation.

See Matlab below.

Ton

>> R = [0.5017 -0.8646 -0.0274; 0.8645 0.5000 0.0522; -0.0314 -0.0499 0.9983]

R =

0.5017 -0.8646 -0.0274
0.8645 0.5000 0.0522
-0.0314 -0.0499 0.9983

>> w = [-0.058898,0.0023,0.998261]'

w =

-0.0589
0.0023
0.9983

>> R*w

ans =

-0.0589
0.0023
0.9983

>> R'*w

ans =

-0.0589
0.0023
0.9983

>> [v,d] = eig(R)

v =

0.0001 + 0.7059i 0.0001 - 0.7059i -0.0589 + 0.0000i
0.7071 + 0.0000i 0.7071 + 0.0000i 0.0023 + 0.0000i
-0.0016 + 0.0417i -0.0016 - 0.0417i 0.9983 + 0.0000i


d =

0.5000 + 0.8661i 0.0000 + 0.0000i 0.0000 + 0.0000i
0.0000 + 0.0000i 0.5000 - 0.8661i 0.0000 + 0.0000i
0.0000 + 0.0000i 0.0000 + 0.0000i 1.0000 + 0.0000i

[ note the third eigenvector is 1 and the third eigenvector is the same as w ]

If you had tried an angular velocity of 2 rad/s, it still would be invariant under this rotation (or inverse rotation):

>> w = 2*w

w =

-0.1178
0.0046
1.9965

>> R*w

ans =

-0.1178
0.0047
1.9966

[ the small error in the final digit is likely due to precision of the numerical values you posted]

Hi Ton,

Thanks for your detailed reply. I needed a few days to remember/learn a few new properties! (i.e. one of the eigenvalues of a rotation matrix is 1). I tried to then spend sometime figuring out a few things, and it all seemed to come back to Euler's theorem and the eigenvector associated with the eigenvalue of 1. The eigenvector associated with the eigenvalue of 1 is the axis about which the rotation in my scenario takes place, and a rotation about this axis leads to the components of the angular velocity been identical in both systems.

Chapter 6 (Euler parameters) in the book by Nikravesh explains this concept well in my opinion:
http://www.u.arizona.edu/~pen/ame553/

The full book is also very useful for multibody dynamics and the full copy can be found on the link above.

Kind regards,

Nicos Haralabidis