Obtaining accelerations of segment CoM

[Michael Woodward]

Hi

We are attempting to extract all the forces (muscle, reaction, inertial and gravitational) acting on a particular body at each frame in time.

We’re currently having issues trying to obtain accurate values for the acceleration components and were wondering if you could point us into the right direction. We’re trying to realize the states to include accelerations. Our understanding is that the acceleration which OpenSim can output (via getSimbodyEngine.getAcceleration) calculates accelerations based on F=ma and therefore the forces acting on the system are required (Muscle forces and Ground Reaction Forces) as well as positions and velocity. Is this correct, or would only GRF suffice (as inverse dynamics does not include muscle forces, but does assume accelerations to be known)?

The current approach we are taking is to update the state frame-by-frame. We can get the model to realize the positions in the form of the joint coordinates(velocities) (state.updQ and state.updU, obtained through state reporter) but I am not sure how to apply an external force to the model. The values attained from just providing position information provide realistic velocity values, but very high acceleration values which can be explained by the lack of forces and the way OpenSim (or Simbody) realizes dynamics through F=ma rather than obtaining accelerations though double differentiation of positions.

Is my thinking correct in terms of having to provide both position and forces to the state in order to calculate the accelerations? And if so, what is the correct process for providing the state with the correct force values in order to be able to execute osimModel.computeStateVariableDerivatives(state)?

Thank you for this, and we’re happy to provide any snippets of code if needed.

Michael

[Dimitar Stanev]

Hi,

We are attempting to extract all the forces (muscle, reaction, inertial and gravitational) acting on a particular body at each frame in time.

For you to calculate the acting forces, the model must be realized to Stage::Dynamics and not acceleration. You can calculate all the body and joint forces as follows:

https://github.com/mitkof6/task-space/b ... el.cpp#L45

which does not include the Coriolis effect.

https://github.com/mitkof6/task-space/b ... l.cpp#L101

Note that in my implementation I disable the muscle forces, before calculating the total forces, you must do otherwise.

We’re currently having issues trying to obtain accurate values for the acceleration components and were wondering if you could point us into the right direction. We’re trying to realize the states to include accelerations. Our understanding is that the acceleration which OpenSim can output (via getSimbodyEngine.getAcceleration) calculates accelerations based on F=ma and therefore the forces acting on the system are required (Muscle forces and Ground Reaction Forces) as well as positions and velocity. Is this correct, or would only GRF suffice (as inverse dynamics does not include muscle forces, but does assume accelerations to be known)?

Can you be more specific why you need the acceleration?

You can add Kinematics and BodyKinematics analyses to you model and use the setRecordAccelerations option.

https://simtk.org/api_docs/opensim/api_ ... c25a8e6095

the getSimbodyEngine.getAcceleration calculates the acceleration of a point on some body.

Is my thinking correct in terms of having to provide both position and forces to the state in order to calculate the accelerations? And if so, what is the correct process for providing the state with the correct force values in order to be able to execute osimModel.computeStateVariableDerivatives(state)?

Please provide more details on what you want to achieve, this is a very general question. In short you provide the muscle controls via a Controller that is added to the model. The forces are computed at a dynamics stage and their effect is accounted in the computation of the state derivatives giving you the current induced acceleration. Calling computeStateVariableDerivatives is a low level manipulation and requires that you know exactly what you are doing. Maybe what you want to achieve can be done without having to call these methods yourself.

Best

[Michael Woodward]

Ultimately, we are attempting to extract all the forces and moments of force and their point of application on the femur for each frame of the gait cycle. This is for the intended purpose of using these specific force vectors as the input to finite element modelling.

We are attempting to verify the values of the extracted force components (FHip, FKnee, FMuscles, FGravity) by also extracting the acceleration, which determines the resultant force (ma) as shown in figure 1. This should equate to the sum of the individually extracted forces, showing that they have been estimated correctly.

figure1.png (12.07 KiB) Viewed 1291 times

For briefness, let’s stick to the linear forces:

• The forces of all the muscles acting on the femur (FMuscles)
• The Joint Reaction Forces acting on the distal and proximal ends of the femur (FHip and FKnee)
• The force of the model’s weight due to gravity (FGravity)
• The resultant force on the model: mass x acceleration (ma)

The equilibrium equation as a result of should therefore be: FHip + FKnee + FMuscles + FGravity - ma = 0

Currently, however, we obtain quite large residual values, indicating we’re not extracting all of the force components correctly.

In an attempt to narrow down what the cause of these residuals is, we have simplified the calculations and have omitted the muscle forces. The free-body diagram therefore now resembles the diagram of figure 2 and the forces in the equation include:

• The intersegmental at the joints (knee and hip) as obtained through inverse dynamics (FInverseDynamics(Hip/Knee))
• The force of the model’s weight due to gravity (FGravity)
• The resultant force (ma)

figure2.png (15.17 KiB) Viewed 1291 times

The equilibrium equation therefore becomes: ∑FInverseDynamics + FGravity - m∙a = 0

Realizing State

Below is the current method used to realize the state. This method is called for each frame of the gait and is ran each time one of the force components is calculated.

Code: Select all

%% Realizing the state %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
% Qframe = Coordinate Orientations at specified frame / Uframe = Velocities at specified frame
Qframe = Q.data(n_frame,:), 
Uframe = U.data(n_frame,:)
 
% Initialize state
state = osimModel.initSystem;
 
% Get number of values in Q
numQ = length(Qframe);
 
% Pre-allocate Q and U in OpenSim vector
Qvec = Vector(numQ, 0);
Uvec = Vector(numQ, 0);
 
% Set entries of Q and U in OpenSim vector
for i = 1:numQ
    Qvec.set(i-1, Qframe(i));
    Uvec.set(i-1, Uframe(i));
end
 
% Set state
state.setQ(Qvec)
state.setU(Uvec)
 
 
% Realize velocities and accelerations
osimModel.computeStateVariableDerivatives(state);

Current methods for obtaining forces:

Joint reaction force:
The intersegmental loads, or joint reaction forces (excluding muscle force contribution), are obtained through running a static optimisation and a Joint Reaction Analysis with a model where the muscle actuators are replaced by a single coordinate actuator for each DoF.

Gravitational force:

Code: Select all

%% Obtaining gravitational direction %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
% get body of interest and its mass
bodyOfInterest = osimModel.getBodySet.get(bodyOfInterest_name);
bodyMass = bodyOfInterest.getMass;
 
% gravity direction in global frame
g_vec = osimModel.getGravity;
ground_body = osimModel.getGroundBody;
 
% transform gravity vector to local reference system 
g_local    = Vec3(0);
 
osimModel.getSimbodyEngine.transform(state, ground_body, g_vec, bodyOfInterest, g_local);
g_local_double(n_frame, 1:3) = [g_local.get(0), g_local.get(1), g_local.get(2)];
 
% computing gravitational force in local reference system at each frame.
Gravity_term = g_local_double * bodyMass;

Acceleration component:
For the resultant inertial force (m∙a), we are attempting to obtain values of acceleration for the femur CoM. Obtaining these values through this method should allow us to multiply by the model mass (presented as an inertial matrix given from the OpenSim API: bodyOfInterest.getInertia(m);) and thus equate them to the other force components to verify them.

Currently, we are struggling to obtain the accelerations of the CoM of the femur through realizing the state. We are currently trying various methods of obtaining the accelerations:

1. Through the function osimModel.getSimbodyEngine.getAcceleration after realising the state at each frame in time. Currently this method overestimates the acceleration, resulting in accelerations that are far larger than would even be physically possible. We believe this is due to not realising the state to the appropriate level as probably not all of the required forces are applied to the model. (Code Below)

Code: Select all

%% Obtaining acceleration (method 1) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
% COM in local
COM_local = Vec3(0);
bodyOfInterest.getMassCenter(COM_local);
 
% generate vector for COM in global
COM_global  = Vec3(0);
acc         = Vec3(0);
acc_local   = Vec3(0);
          
% Transform COM to global 
osimModel.getSimbodyEngine.transformPosition(state, bodyOfInterest, COM_local, ground_body, COM_global);
COM_global_double(n_frame, 1:3) = [COM_global.get(0), COM_global.get(1), COM_global.get(2)];
    
% Get linear acceleration in global
osimModel.getSimbodyEngine.getAcceleration(state, bodyOfInterest, COM_global, acc);
    
% Transform acceleration to local 
osimModel.getSimbodyEngine.transformPosition(state, ground_body, acc, bodyOfInterest, acc_local);
    
acc_local_double(n_frame, 1:3) = [acc_local.get(0), acc_local.get(1), acc_local.get(2)];
 
Resultant_force_term = bodyMass * acc_local_double;

2. Through running a Body Kinematics Analysis and using the accelerations outputted from this process to calculate F = ma. Attached is a copy of the settings file for this analysis.

BKsettings.xml

(4.79 KiB) Downloaded 24 times

In short, our question is what the correct way of obtaining values for the CoM acceleration would be in this case? If this involves realising the model to a specific state, what stage would this be and what would be the required steps to do so? If realizing the state is not the advised way, what would be an alternative?

[Bart van Veen]

Is there anyone who could help with this? I would also be very interested in the answer to this topic.

Thanks in advance.

Bart

[Ryan Byrne]

Hi Michael,

I realize it's been a while - did you ever find a solution to your problem? I also have the same question.

Ryan

[Giulia Ghielmi]

Hi all!
Have you found a solution to this problem? I'm trying to do the equilibrium for each time step and I don't understand how I can calculate the inertial force.
Thank you in advance

Giulia