simtk.org

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

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

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

In this specific model, the coordinates in the multibody tree order would be:
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)

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

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:

Ok got it. Thank you again, Mohammadreza.

Regards,
Shivangi