GRF not being applied correctly to Inverse Dynamics in Python

[Shivangi Giri]

Hello OpenSim Community,

I am trying to implement inverse dynamics in Python. I am using the full body model developed by Rajagopal and their dataset from which I am taking the kinematics and the GRF.
So below is the code I use to run ID

import opensim as osim
mydata = pd.read_excel("ik_output_walk.xlsx",'Sheet1') #Changed IK mot file into xlsx format to extract data. Also converted degrees into
radians
ID_model = osim.Model('Rajagopal2015.osim') #This model has no muscles

#The GRF mot file was divided into two files containing right and left GRF data
grffile_right="grf_walk_for_python_right.mot"
rGRF = osim.ExternalForce(osim.Storage(grffile_right), 'ground_force_r_v', 'ground_force_r_p','ground_torque_r_', 'calcn_r', 'ground', ground')
ID_model.addForce(rGRF)

# The above was throwing an exception, therefore had to use a fix as mentioned here https://github.com/opensim-org/opensim-core/issues/3204
grf_r_storage = osim.Storage('grf_walk_for_python_right.mot', False); # loading the external force data in a storage object
right_GRF_in_model = ID_model.getForceSet().get(0); # getting the force from the model
rGRF.safeDownCast(right_GRF_in_model).setDataSource(grf_r_storage); # resetting the data source
ID_model.finalizeConnections()
grffile_left="grf_walk_for_python_left.mot"
lGRF = osim.ExternalForce(osim.Storage(grffile_left), 'ground_force_l_v', 'ground_force_l_p','ground_torque_l_', 'calcn_l', 'ground', 'ground')
ID_model.addForce(lGRF)
grf_l_storage = osim.Storage('grf_walk_for_python_left.mot', False); # loading the external force data in a storage object
left_GRF_in_model = ID_model.getForceSet().get(1); # getting the force from the model
lGRF.safeDownCast(left_GRF_in_model).setDataSource(grf_l_storage); # resetting the data source
ID_model.finalizeConnections()
ID_state = ID_model.initSystem()
for m in range(ID_model.getCoordinateSet().getSize()): #Setting the initial pose from first row of ik data
ID_model.getCoordinateSet().get(m).setValue(ID_state,float(mydata.iloc[1,m+1]))

#Simulation time from 0.24s to 1.8s
for i in range(2,len(time)):
timed = time[i+1]-time
ID = osim.InverseDynamicsSolver(ID_model)
for j in range(ID_model.getCoordinateSet().getSize()):
model.getCoordinateSet().get(j).setValue(state, mydata[i,j])
model.getCoordinateSet().get(j).setSpeedValue(state, (mydata(i+1,j)-mydata(i,j))/timed)
u_dot = osim.Vector(39,0)
tau = ID.solve(ID_state,u_dot)

But the graphs being produced report very high moments esp. hip and arm flexion moments, whereas ankle moment is very low, indicating the GRF is not being implemented properly. In the figure starting from the top is ankle moment, arm_flex moment, hip_flex moment and then knee_flex moment of right leg.

Please anyone can tell me how I can check whether GRF is being applied correctly? Or is it something else I should look into?
Thanks in advance.

Regards,
Shivangi

[Mohammadreza Rezaie]

Hi, I think your GRF is applied correctly. Two points:

1. InverseDynamicsSolver requires correct acceleration (not zero u_dot). You can simply use the second derivative of q
2. It also requires coordinates acceleration in multibody tree order not CoordinateSet order.

Please see this post: viewtopicPhpbb.php?f=91&t=16680&p=0&sta ... 4f7dccced1

[Shivangi Giri]

Thank you Mohammadreza, for your reply and suggestions.
I have a question regarding your second point. Is multibody tree order different than coordinate set order?
If yes, can you point me to a documentation link or example which can help me look up the multibody tree order.

So far I have tried assigning accelerations according to the Coordinate Set order, as mentioned below
(I know this is a very crass way of assigning the values)
u_dot.set(0, p_tilt_alpha)
u_dot.set(1, p_list_alpha)
u_dot.set(2, p_rot_alpha)
u_dot.set(3, p_tx_alpha)
u_dot.set(4, p_ty_alpha)
u_dot.set(5, p_tz_alpha)
u_dot.set(6, hip_flex_r_alpha)
u_dot.set(7, hip_add_r_alpha)
u_dot.set(8, hip_rot_r_alpha)
u_dot.set(9, knee_flex_r_alpha)
u_dot.set(10, knee_beta_r_alpha)
u_dot.set(11, ankle_r_alpha)
u_dot.set(12, subtalar_r_alpha)
u_dot.set(13, mtp_r_alpha)
u_dot.set(14, hip_flex_l_alpha)
u_dot.set(15, hip_add_l_alpha)
u_dot.set(16, hip_rot_l_alpha)
u_dot.set(17, knee_flex_l_alpha)
u_dot.set(18, knee_beta_l_alpha)
u_dot.set(19, ankle_l_alpha)
u_dot.set(20, subtalar_l_alpha)
u_dot.set(21, mtp_l_alpha)
u_dot.set(22, lumb_ext_alpha)
u_dot.set(23, lumb_bend_alpha)
u_dot.set(24, lumb_rot_alpha)
u_dot.set(25, arm_flex_r_alpha)
u_dot.set(26, arm_add_r_alpha)
u_dot.set(27, arm_rot_r_alpha)
u_dot.set(28, elbow_flex_r_alpha)
u_dot.set(29, pro_sup_r_alpha)
u_dot.set(30, w_flex_r_alpha)
u_dot.set(31, w_dev_r_alpha)
u_dot.set(32, arm_flex_l_alpha)
u_dot.set(33, arm_add_l_alpha)
u_dot.set(34, arm_rot_l_alpha)
u_dot.set(35, elbow_flex_l_alpha)
u_dot.set(36, pro_sup_l_alpha)
u_dot.set(37, w_flex_l_alpha)
u_dot.set(38, w_dev_l_alpha)

However, the torques still look very bad.
How do I assign accelerations according to the multibody tree order?
Thank you for you help in advance.

Regards,
Shivangi

[Mohammadreza Rezaie]

In this specific model, the coordinates in the multibody tree order would be:

Code: Select all

pelvis_tilt
pelvis_list
pelvis_rotation
pelvis_tx
pelvis_ty
pelvis_tz
hip_flexion_r
hip_adduction_r
hip_rotation_r
hip_flexion_l
hip_adduction_l
hip_rotation_l
lumbar_extension
lumbar_bending
lumbar_rotation
knee_angle_r
knee_angle_r_beta
knee_angle_l
knee_angle_l_beta
arm_flex_r
arm_add_r
arm_rot_r
arm_flex_l
arm_add_l
arm_rot_l
ankle_angle_r
ankle_angle_l
elbow_flex_r
elbow_flex_l
subtalar_angle_r
subtalar_angle_l
pro_sup_r
pro_sup_l
mtp_angle_r
mtp_angle_l
wrist_flex_r
wrist_dev_r
wrist_flex_l
wrist_dev_l

If you check the columns of the ID output, you will figure it out. There is a function for this, getCoordinatesInMultibodyTreeOrder, but it doesn't work in Python/MATLAB:
https://github.com/opensim-org/opensim-core/issues/3569

In this way, you can get the multibody tree order. (if you remove .getName(), you will get OpenSim Coordinate class)

Code: Select all

test = [[c.getBodyIndex(), c.getMobilizerQIndex(), c] for c in model.getCoordinateSet()]
multibodyTreeOrder = [i[2].getName() for i in sorted(test)]

[Shivangi Giri]

Thank you very much Mohammadreza for your detailed answer. I was able to calculate the ID torques using the order you listed .
Just one question, is it mandatory to set the time explicitly to model.setTime(time) to make ID pull GRF at the required instances or does the ID solver take care of it automatically?

Thank you for your help in advance.

Regards,
Shivangi

[Mohammadreza Rezaie]

Yes, according to Ayman's answer in this post: viewtopic.php?f=91&t=16680&p=0&start=0& ... f3a#p46315

4. The time is maintained by the state and thus used to query the GRF for values.

It must be:

Code: Select all

state.setTime(time)

[Shivangi Giri]

Ok got it. Thank you again, Mohammadreza.

Regards,
Shivangi