Running Static Optimization tool of OpenSim in matlab live

[Ratna Sambhav]

Hi all,

Is it possible to use OpenSim API commands in matlab to run static optimization but instead of saving the complete result for all the input instances together in an output directory, getting the result at every time instances.
For example, suppose, in every iteration I am giving the value of moment at joint, state values of the joint and other informations, and get the values of activations and excitations for the muscles of the joint. And do it repeatedly in every iteration?

Thankfully,
Ratna Sambhav

[Dimitar Stanev]

Hi Ratna,

Currently, you cannot do this easily. It is possible, but I am not sure to what extend you can rewrite what the analysis tool is doing in Matlab. In C++ this is possible, however, you have to be familiar with OpenSim's and Simbody's API.

[Ratna Sambhav]

Hi Dimitar,
Thanks for the reply. Surely, I'll shift to C++ but how C++ is more flexible than matlab, I mean, when they both are using the same OpenSim classes as on Doxygen.
Regarding Static Optimization in matlab, for sure, I am not able to use OpenSim API to full extent yet. But combining it with matlab fmincon optimization, I am able to get few results. Following is the procedure I am using:

for i=2:120 (number of rows in my input file containing values of position and velocity of joint coordinates)
#Calculated the value of elbow torque (tau1) using a validated equation (consistent with output from OpenSim GUI ID) and shoulder torque (tau2)
from file.

#Setting the position and speed value obtained from file at each of the joint of the model for the current state.

#Integrated using manager for the time duration equal to time gap between two subsequent states.

#Extracted all the muscle properties from the current state which is used in muscle force calculation as shown below and calculated fm1 and fm2
for shoulder and elbow:
for j = 1:6
f_max(j) = model.getMuscles().get(j-1).getMaxIsometricForce();
f_l(j) = model.getMuscles().get(j-1).getActiveForceLengthMultiplier(state);
f_v(j) = model.getMuscles().get(j-1).getForceVelocityMultiplier(state);
f_pe(j) = model.getMuscles().get(j-1).getPassiveForceMultiplier(state);
cosa(j) = model.getMuscles().get(j-1).getCosPennationAngle(state);
r2(i) = model.getMuscles().get(j-1).computeMomentArm(state, model.getCoordinateSet().get(1));
r1(i) = model.getMuscles().get(j-1).computeMomentArm(state, model.getCoordinateSet().get(0));
fm2(j) = f_max(j) * r2(i)*(f_l(j) * f_v(j) + f_pe(j)) * cosa(j)*10;
fm1(j) = f_max(j) * r1(i)*(f_l(j) * f_v(j) + f_pe(j)) * cosa(j)*10;
end

# Finally performed optimizaion using fmincon optimization of matlab to minimize the objective function of sum of square of activation of each
muscle with the constraint equation as [tau1 tau2]' = [muscle activations] * [ [fm1] [fm2] ]'
end

Using this method, I have been able to get muscle activations which is consistent with OpenSim Static Optimization results to some extent. But I am not sure why I had to multiply the fm1 and fm2 equation with an arbitrary value 10. If I don't multiply with 10, then the activation values are going much higher. Is the equation of fm, I am using correct? Or is there some extra term which has to be added?

Thanks,
Ratna

[Michael Asmussen]

Hi Ratna,

I have done a bit of work similar to you (i.e., SO in Matlab), but I did not have to multiply by a value of 10.

I understand that you may not have wrote out your entire optimization routine in your post, but did you set the bounds of optimization such that activations could not be higher than 1 and less than 0? How close are your activations compared to OpenSim's results with and without multiplying by a value of 10?

The only other solution I can think of is that there may be an issue with the units you extracted from the OpenSim API using Matlab.

All the best,

Mike

[Ratna Sambhav]

Hi Mike,

Regarding my code, I perhaps was using a wrong value of elbow torque, which was resulting in wrong output. Now using InverseDynamicsSolver of OpenSim API, the obtained value of shoulder and elbow torque is giving me almost the same muscle activation plot as that obtained from OpenSim without doing any manipulation anywhere.
I have set lower bound as 0.01 and upper bound as 1 for each muscle in fmincon solver.

Glad to know that you have done a similar work. Can you explain your method, how you are doing SO in matlab. That will be very insightful.

Thanks,
Ratna

[Michael Asmussen]

Hi Ratna,

I use the similar method to what you described in an earlier post. Given the posture of the model, I perform the following:

1. extract the maximum force generating capacity of all the muscles of interest crossing the joints of interest
2. extract the same muscles' moment arms
3. bound the activations from 0 to 1
4. set an equality constraint for the solution from static optimization to match the moment input from inverse dynamics
5. use fmincon to solve for the activations

The equality constraint sometimes makes the solver run quite slow, so another option is to run optimizer with an inequality constraint instead. Just let me know if you have any more questions.

Hope this helps!

Mike