Tendon compliance in MocoInverse example

OpenSim Moco is a software toolkit to solve optimal control problems with musculoskeletal models defined in OpenSim using the direct collocation method.
POST REPLY
User avatar
Israel Luis
Posts: 11
Joined: Thu Oct 24, 2019 3:21 am

Tendon compliance in MocoInverse example

Post by Israel Luis » Tue Jun 16, 2020 7:01 am

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.

User avatar
Christopher Dembia
Posts: 506
Joined: Fri Oct 12, 2012 4:09 pm

Re: Tendon compliance in MocoInverse example

Post by Christopher Dembia » Tue Jun 16, 2020 9:57 am

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?

User avatar
Israel Luis
Posts: 11
Joined: Thu Oct 24, 2019 3:21 am

Re: Tendon compliance in MocoInverse example

Post by Israel Luis » Tue Jun 16, 2020 12:07 pm

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

User avatar
Christopher Dembia
Posts: 506
Joined: Fri Oct 12, 2012 4:09 pm

Re: Tendon compliance in MocoInverse example

Post by Christopher Dembia » 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.

User avatar
Nicholas Bianco
Posts: 1051
Joined: Thu Oct 04, 2012 8:09 pm

Re: Tendon compliance in MocoInverse example

Post by Nicholas Bianco » Tue Jun 16, 2020 2:03 pm

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.

User avatar
Ross Miller
Posts: 375
Joined: Tue Sep 22, 2009 2:02 pm

Re: Tendon compliance in MocoInverse example

Post by Ross Miller » Tue Jun 16, 2020 4:53 pm

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

User avatar
Israel Luis
Posts: 11
Joined: Thu Oct 24, 2019 3:21 am

Re: Tendon compliance in MocoInverse example

Post by Israel Luis » Thu Jun 18, 2020 11:14 am

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 :D 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:

Image

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.

User avatar
Christopher Dembia
Posts: 506
Joined: Fri Oct 12, 2012 4:09 pm

Re: Tendon compliance in MocoInverse example

Post by Christopher Dembia » Thu Jun 18, 2020 11:34 am

Thanks for reporting this :)

The TFL has a very long tendon-to-fiber ratio; that may contribute to any issues you're encountering.

POST REPLY