Page 1 of 1
Tendon compliance in MocoInverse example
Posted: Tue Jun 16, 2020 7:01 am
by israelluis7
Hi Moco team,
I am not sure why but the MocoInverse matlab code provided as an example: exampleMocoInverse() (as it is), runs only if I ignore tendon Compliance. It is suggested in the forum to changed tendon compliance to implicit formulation to run the simulation. So I did it using either
modelProcessor.append(ModOpUseImplicitTendonComplianceDynamicsDGF());
or
modelProcessor.append(ModOpTendonComplianceDynamicsModeDGF('implicit'));
and allowing/ignoring passive force. But in all cases, my simulation either fails or stops at the third iteration (for a very long time).
Here is the code:
Code: Select all
import org.opensim.modeling.*;
% Construct the MocoInverse tool.
inverse = MocoInverse();
modelProcessor = ModelProcessor('subject_walk_armless.osim');
modelProcessor.append(ModOpAddExternalLoads('grf_walk.xml'));
modelProcessor.append(ModOpTendonComplianceDynamicsModeDGF('implicit'));
modelProcessor.append(ModOpReplaceMusclesWithDeGrooteFregly2016());
% modelProcessor.append(ModOpIgnorePassiveFiberForcesDGF());
modelProcessor.append(ModOpScaleActiveFiberForceCurveWidthDGF(1.5));
modelProcessor.append(ModOpAddReserves(1.0));
inverse.setModel(modelProcessor);
inverse.setKinematics(TableProcessor('coordinates.sto'));
% Initial time, final time, and mesh interval.
inverse.set_initial_time(0.81);
inverse.set_final_time(1.79);
inverse.set_mesh_interval(0.02);
inverse.set_kinematics_allow_extra_columns(true);
% Solve the problem and write the solution to a Storage file.
solution = inverse.solve();
solution.getMocoSolution().write('example3DWalking_MocoInverse_solution.sto');
When it fails, it says the following:
[CasOC] Warning: NaN encountered when detecting sparsity of Jacobian; entry (133, 224). And also
java.lang.RuntimeException: This trajectory is sealed, to force you to acknowledge the
solver failed; call unseal() to gain access.
Thrown at MocoTrajectory.cpp:1347 in ensureUnsealed().
In addition, Dr. Miller shared a code where he used MocoTrack with tendon Compliance
https://simtk.org/projects/umocod, and he initialized tendon force (initial guess). I tried this approach, but I did not work. Maybe I did not implement it correctly.
I would appreciate it if someone could suggest ways to solve this problem.
Re: Tendon compliance in MocoInverse example
Posted: Tue Jun 16, 2020 9:57 am
by chrisdembia
The equations of motion can be sensitive to the state variable used to model tendon compliance (in this case, normalized tendon force), so providing an initial guess can be helpful.
Could you say more about why using Ross Miller's initial guess approach did not work?
Re: Tendon compliance in MocoInverse example
Posted: Tue Jun 16, 2020 12:07 pm
by israelluis7
Hi Christopher,
Yes, I included the configuration of the solver that Dr. Miller implemented in this code and also provided an initial guess to the normalize tendon force as you can see in the following code:
Code: Select all
import org.opensim.modeling.*;
inverse = MocoInverse();
modelProcessor = ModelProcessor('subject_walk_armless.osim');
modelProcessor.append(ModOpAddExternalLoads('grf_walk.xml'));
modelProcessor.append(ModOpTendonComplianceDynamicsModeDGF('implicit'));
modelProcessor.append(ModOpReplaceMusclesWithDeGrooteFregly2016());
% Only valid for DeGrooteFregly2016Muscles.
% modelProcessor.append(ModOpIgnorePassiveFiberForcesDGF());
% Only valid for DeGrooteFregly2016Muscles.
modelProcessor.append(ModOpScaleActiveFiberForceCurveWidthDGF(1.5));
modelProcessor.append(ModOpAddReserves(1.0));
inverse.setModel(modelProcessor);
inverse.setKinematics(TableProcessor('coordinates.sto'));
inverse.set_initial_time(0.81);
inverse.set_final_time(1.79);
inverse.set_mesh_interval(0.02);
inverse.set_kinematics_allow_extra_columns(true);
study = inverse.initialize();
problem = study.updProblem();
problem.setStateInfoPattern('/forceset/.*/normalized_tendon_force', [0, 1.8], [], []);
% Define the solver and set its options
solver = MocoCasADiSolver.safeDownCast(study.updSolver());
%solver.set_multibody_dynamics_mode('implicit')
%solver.set_minimize_implicit_multibody_accelerations(true)
%solver.set_implicit_multibody_accelerations_weight(0.00001)
solver.set_optim_max_iterations(2000);
solver.set_num_mesh_intervals(100);
solver.set_optim_constraint_tolerance(1e-2);
solver.set_optim_convergence_tolerance(1e-2);
solver.set_minimize_implicit_auxiliary_derivatives(true)
solver.set_implicit_auxiliary_derivatives_weight(0.001)
solver.resetProblem(problem);
model = modelProcessor.process();
model.initSystem();
guess = solver.createGuess();
numRows = guess.getNumTimes();
StateNames = model.getStateVariableNames();
for i = 1:model.getNumStateVariables()
currentStateName = string(StateNames.getitem(i-1));
if contains(currentStateName,'normalized_tendon_force')
MusName = currentStateName;
guess.setState(currentStateName, linspace(0.2,0.2,numRows));
end
end
solver.setGuess(guess);
solution = study.solve();
As a result, the optimization did not succeed. This is the output:
Code: Select all
[CasOC] Warning: NaN encountered when detecting sparsity of Jacobian; entry (133, 224).
List of user-set options:
Name Value used
acceptable_compl_inf_tol = 0.01 yes
acceptable_constr_viol_tol = 0.01 yes
acceptable_dual_inf_tol = 0.01 yes
acceptable_tol = 0.01 yes
compl_inf_tol = 0.01 yes
constr_viol_tol = 0.01 yes
dual_inf_tol = 0.01 yes
hessian_approximation = limited-memory yes
max_iter = 2000 yes
print_user_options = yes yes
tol = 0.01 yes
This is Ipopt version 3.12.8, running with linear solver mumps.
NOTE: Other linear solvers might be more efficient (see Ipopt documentation).
Number of nonzeros in equality constraint Jacobian...: 513520
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 0
CasADi - 2020-06-16 14:54:00 WARNING("nlp:nlp_jac_g failed: NaN detected for output jac_g_x, at nonzero index 88 (row 168, col 0).") [D:\a\opensim-moco\opensim-moco\dependencies\casadi\casadi\core\oracle_function.cpp:265]
Error evaluating Jacobian of equality constraints at user provided starting point.
No scaling factors for equality constraints computed!
CasADi - 2020-06-16 14:54:34 WARNING("nlp:nlp_jac_g failed: NaN detected for output jac_g_x, at nonzero index 88 (row 168, col 0).") [D:\a\opensim-moco\opensim-moco\dependencies\casadi\casadi\core\oracle_function.cpp:265]
Number of Iterations....: 0
Number of objective function evaluations = 0
Number of objective gradient evaluations = 0
Number of equality constraint evaluations = 0
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 1
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 0
Total CPU secs in IPOPT (w/o function evaluations) = 37.644
Total CPU secs in NLP function evaluations = 34.682
EXIT: Invalid number in NLP function or derivative detected.
t_proc [s] t_wall [s] n_eval
nlp 72.3 72.3 1
nlp_grad_f 1.74 1.74 1
Breakdown of objective (including weights):
excitation_effort: 3.87638
Active or violated continuous variable bounds
L and U indicate which bound is active; '*' indicates a bound is violated.
The case of lower==upper==value is ignored.
State bounds: no bounds active or violated
Control bounds: no bounds active or violated
Multiplier bounds: no bounds active or violated
Derivative bounds: no bounds active or violated
Active or violated parameter bounds
L and U indicate which bound is active; '*' indicates a bound is violated.
The case of lower==upper==value is ignored.
Time bounds: no bounds active or violated
Parameter bounds: no bounds active or violated
Total number of constraints: 28898.
Differential equation defects:
L2 norm across mesh, max abs value (L1 norm), time of max abs
/forceset/addbrev_r/activation 9.99e-16 2.22e-16 0.819500
/forceset/addlong_r/activation 9.99e-16 2.22e-16 0.819500
/forceset/addlong_r/normalized_tendon_force 4.94e-04 5.17e-05 0.810000
/forceset/addmagDist_r/activation 9.99e-16 2.22e-16 0.819500
/forceset/addmagIsch_r/activation 9.99e-16 2.22e-16 0.819500
/forceset/addmagIsch_r/normalized_tendon_force 5.97e-04 8.60e-05 0.875550
/forceset/addmagMid_r/activation 9.99e-16 2.22e-16 0.819500
/forceset/addmagProx_r/activation 9.99e-16 2.22e-16 0.819500
/forceset/bflh_r/activation 9.99e-16 2.22e-16 0.819500
/forceset/bflh_r/normalized_tendon_force 9.71e-05 1.82e-05 0.810950
/forceset/bfsh_r/activation 9.99e-16 2.22e-16 0.819500
/forceset/edl_r/activation 9.99e-16 2.22e-16 0.819500
/forceset/edl_r/normalized_tendon_force 3.52e-04 5.69e-05 0.861300
/forceset/ehl_r/activation 9.99e-16 2.22e-16 0.819500
/forceset/ehl_r/normalized_tendon_force 4.17e-04 6.76e-05 0.860350
/forceset/fdl_r/activation 9.99e-16 2.22e-16 0.819500
/forceset/fdl_r/normalized_tendon_force 7.84e-05 9.69e-06 0.904050
/forceset/fhl_r/activation 9.99e-16 2.22e-16 0.819500
/forceset/fhl_r/normalized_tendon_force 1.40e-04 1.78e-05 0.904050
/forceset/gaslat_r/activation 9.99e-16 2.22e-16 0.819500
/forceset/gaslat_r/normalized_tendon_force 5.12e-04 3.03e-05 0.904050
/forceset/gasmed_r/activation 9.99e-16 2.22e-16 0.819500
/forceset/gasmed_r/normalized_tendon_force 5.20e-04 1.49e-06 0.810000
/forceset/glmax1_r/activation 9.99e-16 2.22e-16 0.819500
/forceset/glmax2_r/activation 9.99e-16 2.22e-16 0.819500
/forceset/glmax3_r/activation 9.99e-16 2.22e-16 0.819500
/forceset/glmed1_r/activation 9.99e-16 2.22e-16 0.819500
/forceset/glmed2_r/activation 9.99e-16 2.22e-16 0.819500
/forceset/glmed3_r/activation 9.99e-16 2.22e-16 0.819500
/forceset/glmin1_r/activation 9.99e-16 2.22e-16 0.819500
/forceset/glmin2_r/activation 9.99e-16 2.22e-16 0.819500
/forceset/glmin3_r/activation 9.99e-16 2.22e-16 0.819500
/forceset/glmin3_r/normalized_tendon_force 3.50e-04 7.76e-05 0.904050
/forceset/grac_r/activation 9.99e-16 2.22e-16 0.819500
/forceset/iliacus_r/activation 9.99e-16 2.22e-16 0.819500
/forceset/perbrev_r/activation 9.99e-16 2.22e-16 0.819500
/forceset/perbrev_r/normalized_tendon_force 1.51e-04 2.18e-05 0.904050
/forceset/perlong_r/activation 9.99e-16 2.22e-16 0.819500
/forceset/perlong_r/normalized_tendon_force 9.00e-05 1.24e-05 0.904050
/forceset/piri_r/activation 9.99e-16 2.22e-16 0.819500
/forceset/piri_r/normalized_tendon_force 8.24e-05 2.01e-05 0.904050
/forceset/psoas_r/activation 9.99e-16 2.22e-16 0.819500
/forceset/recfem_r/activation 9.99e-16 2.22e-16 0.819500
/forceset/recfem_r/normalized_tendon_force 1.27e-04 1.73e-05 0.810000
/forceset/sart_r/activation 9.99e-16 2.22e-16 0.819500
/forceset/semimem_r/activation 9.99e-16 2.22e-16 0.819500
/forceset/semimem_r/normalized_tendon_force 1.83e+00 -2.33e-05 0.810000
/forceset/semiten_r/activation 9.99e-16 2.22e-16 0.819500
/forceset/semiten_r/normalized_tendon_force 3.04e-04 2.60e-05 0.810950
/forceset/soleus_r/activation 9.99e-16 2.22e-16 0.819500
/forceset/soleus_r/normalized_tendon_force 5.03e-04 6.62e-05 0.904050
/forceset/tfl_r/activation 9.99e-16 2.22e-16 0.819500
/forceset/tfl_r/normalized_tendon_force -nan(ind) -nan(ind) 0.810000
/forceset/tibant_r/activation 9.99e-16 2.22e-16 0.819500
/forceset/tibant_r/normalized_tendon_force 6.24e-04 1.01e-04 0.860350
/forceset/tibpost_r/activation 9.99e-16 2.22e-16 0.819500
/forceset/tibpost_r/normalized_tendon_force 1.17e-04 1.50e-05 0.904050
/forceset/vasint_r/activation 9.99e-16 2.22e-16 0.819500
/forceset/vasint_r/normalized_tendon_force 5.65e-04 8.06e-05 0.864150
/forceset/vaslat_r/activation 9.99e-16 2.22e-16 0.819500
/forceset/vaslat_r/normalized_tendon_force 5.41e-04 7.91e-05 0.861300
/forceset/vasmed_r/activation 9.99e-16 2.22e-16 0.819500
/forceset/vasmed_r/normalized_tendon_force 5.45e-04 7.73e-05 0.862250
/forceset/addbrev_l/activation 9.99e-16 2.22e-16 0.819500
/forceset/addlong_l/activation 9.99e-16 2.22e-16 0.819500
/forceset/addlong_l/normalized_tendon_force 2.26e-04 4.94e-05 0.904050
/forceset/addmagDist_l/activation 9.99e-16 2.22e-16 0.819500
/forceset/addmagIsch_l/activation 9.99e-16 2.22e-16 0.819500
/forceset/addmagIsch_l/normalized_tendon_force 2.58e-04 -2.62e-06 0.810000
/forceset/addmagMid_l/activation 9.99e-16 2.22e-16 0.819500
/forceset/addmagProx_l/activation 9.99e-16 2.22e-16 0.819500
/forceset/bflh_l/activation 9.99e-16 2.22e-16 0.819500
/forceset/bflh_l/normalized_tendon_force 3.34e-04 7.13e-05 0.810950
/forceset/bfsh_l/activation 9.99e-16 2.22e-16 0.819500
/forceset/edl_l/activation 9.99e-16 2.22e-16 0.819500
/forceset/edl_l/normalized_tendon_force 2.46e-04 2.79e-05 0.833750
/forceset/ehl_l/activation 9.99e-16 2.22e-16 0.819500
/forceset/ehl_l/normalized_tendon_force 2.92e-04 3.31e-05 0.833750
/forceset/fdl_l/activation 9.99e-16 2.22e-16 0.819500
/forceset/fdl_l/normalized_tendon_force 5.23e-05 8.52e-06 0.885050
/forceset/fhl_l/activation 9.99e-16 2.22e-16 0.819500
/forceset/fhl_l/normalized_tendon_force 9.79e-05 1.55e-05 0.886000
/forceset/gaslat_l/activation 9.99e-16 2.22e-16 0.819500
/forceset/gaslat_l/normalized_tendon_force 2.92e-04 4.20e-05 0.885050
/forceset/gasmed_l/activation 9.99e-16 2.22e-16 0.819500
/forceset/gasmed_l/normalized_tendon_force 2.59e-04 3.92e-05 0.886000
/forceset/glmax1_l/activation 9.99e-16 2.22e-16 0.819500
/forceset/glmax2_l/activation 9.99e-16 2.22e-16 0.819500
/forceset/glmax3_l/activation 9.99e-16 2.22e-16 0.819500
/forceset/glmed1_l/activation 9.99e-16 2.22e-16 0.819500
/forceset/glmed2_l/activation 9.99e-16 2.22e-16 0.819500
/forceset/glmed3_l/activation 9.99e-16 2.22e-16 0.819500
/forceset/glmin1_l/activation 9.99e-16 2.22e-16 0.819500
/forceset/glmin2_l/activation 9.99e-16 2.22e-16 0.819500
/forceset/glmin3_l/activation 9.99e-16 2.22e-16 0.819500
/forceset/glmin3_l/normalized_tendon_force 2.70e-04 4.78e-05 0.879350
/forceset/grac_l/activation 9.99e-16 2.22e-16 0.819500
/forceset/iliacus_l/activation 9.99e-16 2.22e-16 0.819500
/forceset/perbrev_l/activation 9.99e-16 2.22e-16 0.819500
/forceset/perbrev_l/normalized_tendon_force 1.10e-04 1.74e-05 0.885050
/forceset/perlong_l/activation 9.99e-16 2.22e-16 0.819500
/forceset/perlong_l/normalized_tendon_force 6.53e-05 1.02e-05 0.885050
/forceset/piri_l/activation 9.99e-16 2.22e-16 0.819500
/forceset/piri_l/normalized_tendon_force 9.16e-05 1.50e-05 0.810000
/forceset/psoas_l/activation 9.99e-16 2.22e-16 0.819500
/forceset/recfem_l/activation 9.99e-16 2.22e-16 0.819500
/forceset/recfem_l/normalized_tendon_force 3.83e-04 8.11e-05 0.816650
/forceset/sart_l/activation 9.99e-16 2.22e-16 0.819500
/forceset/semimem_l/activation 9.99e-16 2.22e-16 0.819500
/forceset/semimem_l/normalized_tendon_force 3.44e-04 6.60e-05 0.819500
/forceset/semiten_l/activation 9.99e-16 2.22e-16 0.819500
/forceset/semiten_l/normalized_tendon_force 5.54e-04 1.07e-04 0.810950
/forceset/soleus_l/activation 9.99e-16 2.22e-16 0.819500
/forceset/soleus_l/normalized_tendon_force 3.63e-04 5.67e-05 0.885050
/forceset/tfl_l/activation 9.99e-16 2.22e-16 0.819500
/forceset/tfl_l/normalized_tendon_force -nan(ind) -nan(ind) 0.810000
/forceset/tibant_l/activation 9.99e-16 2.22e-16 0.819500
/forceset/tibant_l/normalized_tendon_force 4.34e-04 4.97e-05 0.833750
/forceset/tibpost_l/activation 9.99e-16 2.22e-16 0.819500
/forceset/tibpost_l/normalized_tendon_force 8.10e-05 1.30e-05 0.885050
/forceset/vasint_l/activation 9.99e-16 2.22e-16 0.819500
/forceset/vasint_l/normalized_tendon_force 5.94e-04 1.64e-04 0.810000
/forceset/vaslat_l/activation 9.99e-16 2.22e-16 0.819500
/forceset/vaslat_l/normalized_tendon_force 5.16e-04 1.37e-04 0.810000
/forceset/vasmed_l/activation 9.99e-16 2.22e-16 0.819500
/forceset/vasmed_l/normalized_tendon_force 5.91e-04 1.61e-04 0.810000
Kinematic constraints: none
Path constraints: none
-------------------------------------------------------------------------------
Elapsed real time: 245 second(s) (4 minute(s), 5 second(s)).
6/16/2020 2:54:36 PM
MocoCasADiSolver did NOT succeed:
Invalid_Number_Detected
Re: Tendon compliance in MocoInverse example
Posted: Tue Jun 16, 2020 1:55 pm
by chrisdembia
Thanks. The way you are setting the initial guess is great.
Which motion are you studying?
You might enable tendon compliance one muscle at a time to determine which muscle is causing the problem.
Re: Tendon compliance in MocoInverse example
Posted: Tue Jun 16, 2020 2:03 pm
by nbianco
It's weird that you are enabling minimizing auxiliary derivatives, but that term is not showing up in the cost breakdown (only 'excitation_effort').
Also, looking defect constraint breakdown in your output, there's NaN's showing up for both TFL muscles. Might want to look at those muscles first.
Re: Tendon compliance in MocoInverse example
Posted: Tue Jun 16, 2020 4:53 pm
by rosshm
Hi Israel,
I would guess this is probably from a mismatch between the motion being prescribed and the muscle model parameters, e.g. the motion is requiring a muscle's tendon to be stretched really far which results in a very large muscle force, or the motion is requiring a muscle to be really short with unusual fiber/tendon kinematics and this is resulting in strange calculations for muscle force.
Even if using a published model like Rajagopal, it's worthwhile to spend some time making sure the model can produce forces of reasonable magnitudes with reasonable fiber/tendon kinematics across the full ranges of motion you're working with. So if you have not already done that I would give that a try.
Ross
Re: Tendon compliance in MocoInverse example
Posted: Thu Jun 18, 2020 11:14 am
by israelluis7
chrisdembia wrote: ↑Tue Jun 16, 2020 1:55 pm
Thanks. The way you are setting the initial guess is great.
Which motion are you studying?
You might enable tendon compliance one muscle at a time to determine which muscle is causing the problem.
Hi Christopher,
I am analyzing walking with MocoInverse. I use the exact motion data ('coordinates.sto'), ground reaction ('grf_walk.xml') and model ('subject_walk_armless.osim') provided by Moco.
I followed the suggestions from you, Nick and Ross, and it worked
So this is my summary:
It seems that there are certain muscles' tendon compliances that are more difficult to handle for casADi. If I only allow the tendon compliance of the soleus or soleus and lateral and medial gastrocnemius it worked with the previous code I shared. Here I show you the tendon strain:
I selected to ignore one specific muscle's tendon compliance by outputting the model after converting it to DeGrooteFregly muscle type, and changing
<ignore_tendon_compliance> to false (then you run the simulation with this model).
There is a curious fact about it, the first time I converted the model I added "modelProcessor.append(ModOpAddExternalLoads('grf_walk.xml'));" before outputting it (accidentally), and then I run the simulation. This simulation had a duration of 1 hour and 30 minutes, which is dramatically higher compared to the 5 min that usually takes. So, I created the model without this, and as it was expected, it took much lesser with tendon compliance.
Finally, the casADi and guess configuration that the code I shared had does not affect the convergence, at least for simulating soleus and lateral and medial gastrocnemius tendon compliance. But it affects the time. Thus, considering both (casADi and guess), it takes 05:44 min; only casADi, 11:31 min; only guess, 07:52 min; and none, 08:06 min. All provide the same results.
Re: Tendon compliance in MocoInverse example
Posted: Thu Jun 18, 2020 11:34 am
by chrisdembia
Thanks for reporting this
The TFL has a very long tendon-to-fiber ratio; that may contribute to any issues you're encountering.