non-acurate pelvis velocities in local frame. Is this due to Gimbal lock?

[Sina Porsa]

Hi everyone.
I am trying to solve a problem where I need to calculate Angular velocities in local frames. I just realized that pelvis's orientation can affect it's local angular velocities which does not make sense to me. Please take a look at thesimplified MATLAB script. In this case my model has 1 body only (pelvis) which has 6DOF (3 rotations followed by 3 translations)

Code: Select all

pos = [deg2rad([0 0 0]) 0 0 0];
vel = [1 2 3 0 0 0];
for coord = 1:CoordSet.getSize()
    CoordSet.get(coord-1).setValue(states,pos(coord));
    CoordSet.get(coord-1).setSpeedValue(states,vel(coord));
end

model.computeStateVariableDerivatives(states);
v_global = BodySet.get('pelvis').getAngularVelocityInGround(states);
v_local = Ground.expressVectorInAnotherFrame(states, v_global, BodySet.get('pelvis'))

I am expecting v_local to be exactly equal to the first 3 components of vel array. I just realized v_local is a function of 2nd and 3rd orientations which does not mae sense to me. how come pelvis rotation value can affect pelvis list angular velocity? I dont's think this is a Gimbal lock as even pelvis_list=5 degrees can cause inaccurate results. (see attached jpeg )
file

Capture.JPG (98.4 KiB) Viewed 495 times

[Renee Degutis]

Hi everyone.
I am trying to solve a problem where I need to calculate Angular velocities in local frames. I just realized that pelvis's orientation can affect it's local angular velocities. Please take a look at thesimplified MATLAB script. In this case my model has 1 body only (pelvis) which has 6DOF (3 rotations followed by write my essay in 3 hours translations)

Code: Select all

pos = [deg2rad([0 0 0]) 0 0 0];
vel = [1 2 3 0 0 0];
for coord = 1:CoordSet.getSize()
    CoordSet.get(coord-1).setValue(states,pos(coord));
    CoordSet.get(coord-1).setSpeedValue(states,vel(coord));
end

model.computeStateVariableDerivatives(states);
v_global = BodySet.get('pelvis').getAngularVelocityInGround(states);
v_local = Ground.expressVectorInAnotherFrame(states, v_global, BodySet.get('pelvis'))

I am expecting v_local to be exactly equal to the first 3 components of vel array. I just realized v_local is a function of 2nd and 3rd orientations which does not mae sense to me. how come pelvis rotation value can affect pelvis list angular velocity? I dont's think this is a Gimbal lock as even pelvis_list=5 degrees can cause inaccurate results. (see attached jpeg )
fileCapture.JPG

Hello Sina,

I don't think it's due to Gimbal lock either. It does not really lock anything, only for certain operations (from a mathematical point if view).

Regards,
Renee

[Ton van den Bogert]

If I understand correctly what you asked Opensim to do, this is as expected. The time derivatives of the three orientation angles are not the same as the cartesian components of the angular velocity vector (in either reference frame).

See, for example, https://en.wikipedia.org/wiki/Angular_v ... ler_angles

Below is Matlab code to find the relationship between angle derivatives and angular velocity in symbolic form. This is for one particular rotation sequence, but it is easily changed to another sequence. It finds the matrix B such that omega = B*qdot, where qdot are the time derivatives of the orientation angles.

Ton van den Bogert

Code: Select all

syms R Rx Ry Rz 
q    = sym('q',[3,1], 'real');            % q is a symbolic 3x1 matrix
qdot = sym('qdot',[3,1], 'real');         % another symbolic 3x1 matrix
Rx = [  1        0           0              ;     ...
        0    cos(q(1))  -sin(q(1))      ;     ...
        0    sin(q(1))   cos(q(1))      ]
    
Ry = [  cos(q(2))   0    sin(q(2))      ;     ...
             0        1         0           ;     ...
       -sin(q(2))   0    cos(q(2))      ]  
   
Rz = [  cos(q(3))  -sin(q(3))   0       ;     ...
        sin(q(3))   cos(q(3))   0       ;     ...
        0             0        1       ] 
R = Rx * Ry * Rz

% find Rdot using the multivariate chain rule
Rdot = diff(R,q(1)) * qdot(1) + diff(R,q(2)) * qdot(2) + diff(R,q(3)) * qdot(3)

% find omega using R and Rdot
A = R' * Rdot
omega = [-A(2,3) ; A(1,3) ; -A(1,2)]

% omega is linear in qdot, so find the coefficient matrix B
B = jacobian(omega,qdot)
B = simplify(B)