tendon compliance runtime error

[Christian Greve]

Dear Nick and Ross,

I now performed the optimization with tendon-compliance on individual muscles (gastroc, peroneus longus and tibialis posterior). If I only allow tendon compliance for one individual muscle at a time the solution does converge. However, once I run the optimization with more than 1 muscle it fails (Restoration_failed). I have two questions right now:

1. Is there any other parameter which I can tune to allow simultaneous optimization with multiple muscles having tendon-compliance turned on? Or is the best solution right now to iteratively repeat the simulation and with each iteration add tendon compliance to one individual muscle.

2. In the Moco output .sto file there is the variable "/forceset/flex_hal_r" and "/forceset/flex_hal_r/activation". They contain the exact same values. Is there a convenient way to also receive muscle force estimations and information on other muscle tendon parameters from the solver?

Thanks and kind regards

Christian

[Nicholas Bianco]

Hi Christian,

What is happening in the solution that fails to converge? Does the normalized tendon force get very small or large, or is the tendon buckling? You could try optimizing optimal fiber lengths and tendon slack lengths use MocoParameter to see if that helps get a solution that converges.

2. In the Moco output .sto file there is the variable "/forceset/flex_hal_r" and "/forceset/flex_hal_r/activation". They contain the exact same values.

The first variable is muscle excitation and the second is muscle activation. Did this solution converge? These values should be different due to the delays in muscle activation dynamics. You can access other muscle quantities from your solution using the "analyze()" function on MocoStudy.

-Nick

[Christian Greve]

Dear Nick,

@MocoParameters: I did a ParameterOptimization but still the solution would not converge....I think the first step would now be to get more grip on the underlying causes as you suggest as well. However, I have some trouble to implement the .analyze() function.

What I do: I use MocoInverse, append the kinematics and then use study= inverse.initialize to modifie the goals. Then to get the required outputs I use the code below.

Code: Select all

    
    solution = study.solve()
    solution.unseal()
    outputs = osim.StdVectorString()
    outputs.append('.*normalized_tendon_force')
    table = study.analyze(solution, outputs)
    solution.write(Output_Filename)

The error message I get is: std::exception in 'OpenSim::TimeSeriesTable OpenSim::MocoStudy::analyze(OpenSim::MocoTrajectory const &,std::vector< std::string,std::allocator< std::string > >) const': The following 60 states from Model 'Clubfoot_HipCoordActu_hipJointConstraints_NoPathActuator_NoHipMuscles' are missing from the data:****here follows a list with values and speeds: e.g. /jointset/ground_pelvis/pelvis_tilt/value****

I guess that the problem is, that the solution does not contain the required states (values and speeds)? Can I add the states to the solution? Or should I use OpenSim's Analyze tool?


Thanks

Christian
////////////////////////////////////////
The entire code is here:

Code: Select all

modelProcessor = osim.ModelProcessor(Scaled_modelFileName)#('Patient_SCALED.osim')
modelProcessor.append(osim.ModOpAddExternalLoads(GRF_SourceFilename))#
modelProcessor.append(osim.ModOpIgnoreTendonCompliance())
modelProcessor.append(osim.ModOpReplaceMusclesWithDeGrooteFregly2016())
modelProcessor.append(osim.ModOpScaleActiveFiberForceCurveWidthDGF(1.5))
modelProcessor.append(osim.ModOpAddReserves(ReserveActuatorStrength))

#%% Construct MocoInverse
inverse = osim.MocoInverse()
inverse.setModel(modelProcessor)
inverse.setKinematics(osim.TableProcessor(Coordinates_Filename))
inverse.set_minimize_sum_squared_activations(True)


# Initial time, final time, and mesh interval.
inverse.set_initial_time(start_time)
inverse.set_final_time(end_time)
inverse.set_mesh_interval(0.02)

# set kinematics
inverse.set_kinematics_allow_extra_columns(True)
inverse.append_output_paths('.*normalized_tendon_force')
inverse.append_output_paths('.*normTendonForce')   
study = inverse.initialize()
problem = study.updProblem()

# %% modify effort goal
controlEffortWeight = 5
# problem = study.updProblem()
# problem.setModel(Clubfoot_Model)
effort = osim.MocoControlGoal("control_effort") #.safeDownCast(problem.updGoal("control_effort"))
effort.setName('effort')
effort.setWeight(controlEffortWeight)
problem.addGoal(effort)


# Tendon velocity equilibrium
#problem = study.updProblem()
TendonVGoal = osim.MocoInitialVelocityEquilibriumDGFGoal()
TendonVGoal.setName("initial_velocity_equilibrium")
#The problem converges in fewer iterations when this goal is in cost mode.
TendonVGoal.setMode("cost")
TendonVGoal.setWeight(0.001)
problem.addGoal(TendonVGoal)

#%% Add electromyography tracking.
emgTracking = osim.MocoControlTrackingGoal('emg_tracking')
emgTracking.setWeight(10.0)
controlsRef = osim.TimeSeriesTable(EMG_Filename)
emgTracking.setReference(osim.TableProcessor(controlsRef))
emgTracking.setReferenceLabel('/forceset/soleus_r', 'soleus_r_activation')
emgTracking.setReferenceLabel('/forceset/med_gas_r', 'med_gas_r_activation')
emgTracking.setReferenceLabel('/forceset/lat_gas_r', 'med_gas_r_activation')
emgTracking.setReferenceLabel('/forceset/tib_ant_r', 'tib_ant_r_activation')
emgTracking.setReferenceLabel('/forceset/per_brev_r', 'per_brev_r_activation')
emgTracking.setReferenceLabel('/forceset/per_long_r', 'per_long_r_activation')# brev en long are same data
#individual weights
emgTracking.setWeightForControl('/forceset/soleus_r',10)
emgTracking.setWeightForControl('/forceset/med_gas_r',10)
emgTracking.setWeightForControl('/forceset/lat_gas_r',10)
emgTracking.setWeightForControl('/forceset/tib_ant_r',10)
emgTracking.setWeightForControl('/forceset/per_brev_r',5)
emgTracking.setWeightForControl('/forceset/per_long_r',5)

problem.addGoal(emgTracking)


#%% Solve the problem and write the solution to a Storage file.
solver = osim.MocoCasADiSolver.safeDownCast(study.updSolver())
solver.resetProblem(problem)
solver.set_multibody_dynamics_mode('implicit')
solver.set_num_mesh_intervals(40)
solver.set_optim_convergence_tolerance(.01)#.....10^-1
solver.set_optim_constraint_tolerance(1e-4)#1e-4 is default
solver.set_optim_max_iterations(10000)
solver.set_parameters_require_initsystem(False)
solver.setGuessFile(guessFilename)
# guess = solver.createGuess()
# solver.setGuess(guess)
solution = study.solve()
solution.unseal()
outputs = osim.StdVectorString()
outputs.append('.*normalized_tendon_force')
table = study.analyze(solution, outputs)
solution.write(Output_Filename)

[Ross Miller]

Hi Christian,

Are those three muscles you mentioned the only three in the model with tendon compliance? If so then I'm not sure why they would work individually but not as a group. If not then I would guess this indicates the problem is in one of the other muscles.

I agree with Nick's suggestion on examining the output of a failed run. My guess is one or more of the values for tendon slack length is too short. It's often necessary in my experience to make small/medium changes to these values for new use cases (e.g. optimal control simulations) even when using previously "validated" models e.g. Rajagopal.

Ross

[Christian Greve]

Thanks a lot Ross,

@investigating failed solution. I do have some trouble to implement the .analyze() function in the MocoStudy. I think the problem is that there are no states in the solution. At least that is what I make out of the error message below. Could you give me some hint how to add the states to the solution output? My entire code is below the error message.

Regards

Christian


///////////////////////////

Error message when using analyze:

std::exception in 'OpenSim::TimeSeriesTable OpenSim::MocoStudy::analyze(OpenSim::MocoTrajectory const &,std::vector< std::string,std::allocator< std::string > >) const': The following 60 states from Model 'Clubfoot_HipCoordActu_hipJointConstraints_NoPathActuator_NoHipMuscles' are missing from the data:****here follows a list with values and speeds: e.g. /jointset/ground_pelvis/pelvis_tilt/value****

Code: Select all

modelProcessor = osim.ModelProcessor(Scaled_modelFileName)#('Patient_SCALED.osim')
modelProcessor.append(osim.ModOpAddExternalLoads(GRF_SourceFilename))#
modelProcessor.append(osim.ModOpIgnoreTendonCompliance())
modelProcessor.append(osim.ModOpReplaceMusclesWithDeGrooteFregly2016())
modelProcessor.append(osim.ModOpScaleActiveFiberForceCurveWidthDGF(1.5))
modelProcessor.append(osim.ModOpAddReserves(ReserveActuatorStrength))

#%% Construct MocoInverse
inverse = osim.MocoInverse()
inverse.setModel(modelProcessor)
inverse.setKinematics(osim.TableProcessor(Coordinates_Filename))
inverse.set_minimize_sum_squared_activations(True)


# Initial time, final time, and mesh interval.
inverse.set_initial_time(start_time)
inverse.set_final_time(end_time)
inverse.set_mesh_interval(0.02)

# set kinematics
inverse.set_kinematics_allow_extra_columns(True)
inverse.append_output_paths('.*normalized_tendon_force')
inverse.append_output_paths('.*normTendonForce')   
study = inverse.initialize()
problem = study.updProblem()

# %% modify effort goal
controlEffortWeight = 5
# problem = study.updProblem()
# problem.setModel(Clubfoot_Model)
effort = osim.MocoControlGoal("control_effort") #.safeDownCast(problem.updGoal("control_effort"))
effort.setName('effort')
effort.setWeight(controlEffortWeight)
problem.addGoal(effort)


# Tendon velocity equilibrium
#problem = study.updProblem()
TendonVGoal = osim.MocoInitialVelocityEquilibriumDGFGoal()
TendonVGoal.setName("initial_velocity_equilibrium")
#The problem converges in fewer iterations when this goal is in cost mode.
TendonVGoal.setMode("cost")
TendonVGoal.setWeight(0.001)
problem.addGoal(TendonVGoal)

#%% Add electromyography tracking.
emgTracking = osim.MocoControlTrackingGoal('emg_tracking')
emgTracking.setWeight(10.0)
controlsRef = osim.TimeSeriesTable(EMG_Filename)
emgTracking.setReference(osim.TableProcessor(controlsRef))
emgTracking.setReferenceLabel('/forceset/soleus_r', 'soleus_r_activation')
emgTracking.setReferenceLabel('/forceset/med_gas_r', 'med_gas_r_activation')
emgTracking.setReferenceLabel('/forceset/lat_gas_r', 'med_gas_r_activation')
emgTracking.setReferenceLabel('/forceset/tib_ant_r', 'tib_ant_r_activation')
emgTracking.setReferenceLabel('/forceset/per_brev_r', 'per_brev_r_activation')
emgTracking.setReferenceLabel('/forceset/per_long_r', 'per_long_r_activation')# brev en long are same data
#individual weights
emgTracking.setWeightForControl('/forceset/soleus_r',10)
emgTracking.setWeightForControl('/forceset/med_gas_r',10)
emgTracking.setWeightForControl('/forceset/lat_gas_r',10)
emgTracking.setWeightForControl('/forceset/tib_ant_r',10)
emgTracking.setWeightForControl('/forceset/per_brev_r',5)
emgTracking.setWeightForControl('/forceset/per_long_r',5)

problem.addGoal(emgTracking)


#%% Solve the problem and write the solution to a Storage file.
solver = osim.MocoCasADiSolver.safeDownCast(study.updSolver())
solver.resetProblem(problem)
solver.set_multibody_dynamics_mode('implicit')
solver.set_num_mesh_intervals(40)
solver.set_optim_convergence_tolerance(.01)#.....10^-1
solver.set_optim_constraint_tolerance(1e-4)#1e-4 is default
solver.set_optim_max_iterations(10000)
solver.set_parameters_require_initsystem(False)
solver.setGuessFile(guessFilename)
# guess = solver.createGuess()
# solver.setGuess(guess)
solution = study.solve()
solution.unseal()
outputs = osim.StdVectorString()
outputs.append('.*normalized_tendon_force')
table = study.analyze(solution, outputs)
solution.write(Output_Filename)

[Ross Miller]

There's probably a way to do this directly from a MocoSolution, but in my simulations I get the states as Matlab variables by writing the solution to a file then reading the file back in:

Code: Select all

% Solve the problem
soln = study.solve();
    
% Write the solution to a file
soln.write(OutFileName);

% Read solution in
temp = importdata(OutFileName);
states_data = temp.data; % These are the values of the states
states_headers = temp.colheaders; % These are the names of the states

[Christian Greve]

Thanks for the reply.

But the problem is that the coordinate values (e.g. ankle angle) and coordinate speeds are not in the solution. When calling study.solve it only computes and writes the activations, excitations and reserve forces.

Is that because I am using MocoInverse to prescribe the motion and then edit it as a MocoStudy?
study = inverse.initialize()


I can generate muscle outputs when using only MocoInverse and also the solution converges when adding tendon compliance. But when using MocoInverse only I am not tracking EMGs. After parameter optimization, the normalized_fiber_lengths and muscle activations do look fine. Except for one muscle. Please see attachments.

Since MocoInverse converges when adding tendon compliance is there anything from when adding the EMG tracking goal?


Thanks

Christian

[Ross Miller]

The solution file I think will include the time-varying states and controls. With MocoInverse the generalized coordinates and speeds I think are not states, so they won't be in there.

I've never used EMG tracking in Moco but it's good to do things like you described to isolate the problem, e.g. does it converge with/without EMG tracking goal? For "large" problems sometimes there's nothing wrong with the model, it's just a hard problem to solve. E.g. I just did a simulation of transfemoral limb loss that converge on a metabolic cost of walking of ~3.5 J/m/kg (reasonable result). I then added a powered ankle (should reduce cost) and it converged on a metabolic cost of ~11 J/m/kg. The problem here is likely that I need to try a different initial guess, or a different weigh on the ankle actuator.

Ross

[Christian Greve]

Dear Ross and Nick,

I did some debugging on my code and figures that the TendonVGoal is causing trouble. I set the initial weight to 0.001 and it would not converge....but decreasing the weight even further it converges within 90 iterations. I am not sure whether the VGoal still does something but at least its' converging now with tendon compliance....

Thanks for your help and support

Regards

Christian

Code: Select all

TendonVGoal2 = osim.MocoInitialVelocityEquilibriumDGFGoal()
    TendonVGoal2.setName("initial_velocity_equilibrium2")
    # The problem converges in fewer iterations when this goal is in cost mode.
    TendonVGoal2.setMode("cost")
    TendonVGoal2.setWeight(0.0000001)
    problem2.addGoal(TendonVGoal2)

[Nicholas Bianco]

Ah yes, that MocoGoal is somewhat experimental and not documented very well. I'm not too surprised it's causing inconsistent behavior.

Glad the problem is working now!