simtk.org

Dear all,

Background: I am performing parameter optimization on various muscle-tendon and ligament parameters (e.g. optimal_muscle_fibre_length and tendon_slack_length). I am using MocoInverse to prescribe the desired motion. The solver suceeds and provides realistic parameter values.

Problem: However, when I load the motion from the MocoStudy solution it does not resemble the prescribed motion/coordinates.

Question: In what ways can you verify the output from Parameter optimization with Moco? Can the output motion/coordinates from the MocoStudy solution deviate from the input/prescribed motion?

kind regards

Christian

Hi Christian,

Could you be more specific about how the solution doesn't resemble the motion? The input trajectory is splined which allows MocoInverse to prescribe the trajectory; there might be small numerical differences between the original trajectory and splined trajectory but it should be almost identical.

Since you're optimizing muscle-tendon parameters, these values will have no effect on the output kinematics, so this shouldn't be the problem.

-Nick

Hi Nick,

Thanks for the fast reply.

When I load the solution from ParameterOptimization the motion is completely different. Or in other terms, the model "goes crazy". Below I also added the code...maybe I am doing something strange here.

Regards

Christian

###########

import opensim as osim
# from tkinter import filedialog
# import pyplot as plt
import os
#
Path = '/groups/umcg-rehabilitation/tmp01/umcg-cgreve/MocoInverse/'#'C:/Users/Greve\OneDrive - UMCG/UMCG\Onderzoek\OffRoad2021_Clubfoot\HPC_Gearshift\TestData'#
os.chdir(Path)
Scaled_modelFileName = os.path.join(Path, 'CluFO_01_adjustedRRA_ADJUSTED_NoPathActu_HipCoordActu.osim')

GRF_SourceFilename = os.path.join(Path,'grf_CluFO_01_Okt202109_Filt.xml')
Coordinates_Filename = os.path.join(Path,'CluFO_01_IK_Input_Kinematics_q.sto')
Output_Filename = 'MocoInverse_Okt202109_Clubfoot.sto'
# guessFilename = os.path.join(Path,'GuessFile.sto')
EMG_Filename = os.path.join(Path, 'EMG_Clubfoot_CSV.sto')
ReserveActuatorStrength = 25
start_time = 0.59
end_time = 1.49

# This initial block of code is identical to the code above.
inverse = osim.MocoInverse()
modelProcessor = osim.ModelProcessor(Scaled_modelFileName)
modelProcessor.append(osim.ModOpAddExternalLoads(GRF_SourceFilename)) #---------------------------------------------------------------------#
# modelProcessor.append(osim.ModOpIgnoreTendonCompliance())# I should enable tendon compliance, see notes.
modelProcessor.append(osim.ModOpReplaceMusclesWithDeGrooteFregly2016())
modelProcessor.append(osim.ModOpUseImplicitTendonComplianceDynamicsDGF())# Set the tendon compliance dynamics mode to "implicit" for all DeGrooteFregly2016Muscles in the model
# Only valid for DeGrooteFregly2016Muscles.
#modelProcessor.append(osim.ModOpIgnorePassiveFiberForcesDGF())# do I need to enable?
# Only valid for DeGrooteFregly2016Muscles.
modelProcessor.append(osim.ModOpScaleActiveFiberForceCurveWidthDGF(1.5))#shoudl not affect anything?
modelProcessor.append(osim.ModOpAddReserves(ReserveActuatorStrength)) #how high
ClubfootModel = modelProcessor.process()
# #
inverse.setModel(modelProcessor)
inverse.setKinematics(osim.TableProcessor(Coordinates_Filename))
inverse.set_initial_time(start_time)
inverse.set_final_time(end_time)
inverse.set_mesh_interval(0.02)
inverse.set_kinematics_allow_extra_columns(True)
# inverse.set_allow_unused_references(True)

study = inverse.initialize()


MuscleNames = list(['tib_post_r', 'flex_hal_r','flex_dig_r','soleus_r','med_gas_r','lat_gas_r','ext_dig_r', 'per_brev_r',
'per_long_r', 'per_tert_r','ext_hal_r','tib_ant_r',
'tib_post_l', 'flex_hal_l','flex_dig_l','soleus_l','med_gas_l','lat_gas_l','ext_dig_l', 'per_brev_l',
'per_long_l', 'per_tert_l','ext_hal_l','tib_ant_l'])

# Parameter optimization
problem = study.updProblem()
problem.setModel(ClubfootModel)
Muscles = ClubfootModel.updMuscles()
for cntrMuscle in range(Muscles.getSize()):
muscle = Muscles.get(cntrMuscle)
if str(muscle.getName()) in MuscleNames:
# print(str(muscle.getName()))
# minimize influence of muscle activity on optimization--
# controlString ='/forceset/'+str(muscle.getName())
# problem.setControlInfo(controlString, [0.001, 0.01])#Nick: problem.setControlInfo('/control/path', [0.001, 0.01])
Min_TendonSlackLength = muscle.getTendonSlackLength()*0.5
Max_TendonSlackLength = muscle.getTendonSlackLength()*1.5
parameter=osim.MocoParameter(str(muscle.getName())+'_tendon_slack_length', ['/forceset/'+str(muscle.getName())], 'tendon_slack_length',osim.MocoBounds(Min_TendonSlackLength,Max_TendonSlackLength))
problem.addParameter(parameter)

Min_OptimLength = muscle.getOptimalFiberLength()*0.5
Max_OptimLength = muscle.getOptimalFiberLength()*1.5
parameter=osim.MocoParameter(str(muscle.getName())+'_optimal_fibre_length', ['/forceset/'+str(muscle.getName())], 'optimal_fiber_length',osim.MocoBounds(Min_OptimLength,Max_OptimLength))
problem.addParameter(parameter)

# optimize contraction velocity?
Min_MaxContractionVelocity = muscle.getMaxContractionVelocity()*0.8
Max_MaxContractionVelocity = muscle.getMaxContractionVelocity()*2.5
parameter=osim.MocoParameter(str(muscle.getName())+'_max_contraction_velocity', ['/forceset/'+str(muscle.getName())], 'max_contraction_velocity',osim.MocoBounds(Min_MaxContractionVelocity,Max_MaxContractionVelocity))
problem.addParameter(parameter)

# pennation angle
Min_PennationAngle = muscle.getPennationAngleAtOptimalFiberLength()*0.5
Max_PennationAngle = muscle.getPennationAngleAtOptimalFiberLength()*1.5
parameter=osim.MocoParameter(str(muscle.getName())+'_pennation_angle', ['/forceset/'+str(muscle.getName())], 'pennation_angle_at_optimal',osim.MocoBounds(Min_PennationAngle,Max_PennationAngle))
problem.addParameter(parameter)

#%% Solve the problem and write the solution to a Storage file.
solver = study.initCasADiSolver()
solver.set_multibody_dynamics_mode('implicit')
solver.set_num_mesh_intervals(20)
solver.set_optim_convergence_tolerance(.01)#.....10^-1
solver.set_optim_constraint_tolerance(1e-3)#1e-4 is default
solver.set_optim_max_iterations(4000)
solver.set_parameters_require_initsystem(False)
# solver.setGuessFile(guessFilename)
solution = study.solve()
solution.unseal()
solution.write('Clubfoot_MocoInverse_ParameterOptimization_Okt202109.sto')

Hi Christian,

I see the problem. When you use the "initialize()" approach, you need to make sure you use the MocoCasADiSolver that was pre-configured within MocoInverse and update the solver with the modified problem.

This line creates a new solver that has no knowledge of the original problem: You should instead do this: The "study.updSolver()" part accesses the pre-configured MocoCasADiSolver, but it is returned with type MocoSolver, so you need "safeDownCast" to cast it to MocoCasADiSolver. Then you need to update the solver with the modified MocoProblem via "resetProblem".

Your current approach is erasing the original MocoProblem which contains the prescribed kinematics constraint created by MocoInverse. That's why your solution is "going crazy".

Let me know if that works!
-Nick

Hi Nick,

I'm tacking this slightly out of field question onto this thread as it reminded me of something I wanted to ask. It seems like you can convert MocoInverse into a MocoStudy where you can add additional goals to the problem, and the spline-based approach of keeping the kinematics consistent in MocoInverse will remain? For example, could I do this if I wanted to keep the kinematics the same and minimise a joint reaction? I understand you could achieve something very similar with state tracking in a MocoStudy, but this still would have the option of allowing the kinematics to fluctuate a little vs. with MocoInverse?

Aaron

Hi Aaron,

I think that's technically possible in general with MocoInverse, I haven't tried MocoInverse with the joint reaction goal personally though.

The usual OpenSim IK/RRA workflow allows the model's kinematics to deviate from the prescribed kinematics a little to reduce the residuals, although I'm not sure if MocoInverse does that too. I think it only deviates them in the situation Nick described (mismatch in kinematic constraints between data vs. problem).

Tangent: I think MocoTrack would be a great alternative to inverse dynamics for motion data analysis because I find tracking errors easier to interpret on questions like "Is this small enough?" than residuals and reserve torques. It also makes the ground contact modeling problem harder / more important I guess, but it would be fast enough in most cases that studying that sensitivity wouldn't be too onerous most of the time.

Ross

Hi Aaron,

Yes, definitely! That is exactly what we have the initialize() interface for. In fact, one of the examples in the Moco paper showcases exactly this application: minimizing knee joint contact forces (Fig 7). The code for the paper example is located here: minimizing joint reactions with MocoInverse.

Unlike RRA/CMC, MocoInverse does not allow the kinematics to deviate slightly from the prescribed trajectory. We use Simbody's "Motion" class, which prescribes kinematics exactly and then computes the necessary constraint forces to enforce the motion (it basically works like inverse dynamics). We enforce that the constraint forces are zero so that the motion is produced only by the forces in the model. Model kinematic constraints are ignored when using MocoInverse; we assume that the prescribed kinematics already obey all model constraints.

Like you say, MocoTrack can be a great option for producing a dynamically consistent simulation while tracking data. If you don't want to model foot-ground constact, you can still run MocoTrack while applying GRF data via ExternalLoads in the model. One of the problems in exampleMocoTrack demonstrates how to do this.

Best,
Nick

Dear all,

I updated and tested the code again. It performs better but still the solution from Parameter optimization does not come even close to the prescribed coordinates. Even though the simulation succeeded. The entire model is moving in the right direction but there is no change in hip, knee ankle coordinates.

Q1: Do I need to tighten the constraints and convergence tolerance?

Below is the code I used...to speed up the simulation I only included optim_fibre_length as a parameter for now.

Thanks for your help

Christian

#####################################################################
inverse = osim.MocoInverse()
modelProcessor = osim.ModelProcessor(Scaled_modelFileName)
modelProcessor.append(osim.ModOpAddExternalLoads(GRF_SourceFilename)) #---------------------------------------------------------------------#
modelProcessor.append(osim.ModOpIgnoreTendonCompliance())# I should enable tendon compliance, see notes.
modelProcessor.append(osim.ModOpReplaceMusclesWithDeGrooteFregly2016())
modelProcessor.append(osim.ModOpUseImplicitTendonComplianceDynamicsDGF())# Set the tendon compliance dynamics mode to "implicit" for all DeGrooteFregly2016Muscles in the model
#modelProcessor.append(osim.ModOpIgnorePassiveFiberForcesDGF())# do I need to enable?
modelProcessor.append(osim.ModOpScaleActiveFiberForceCurveWidthDGF(1.5))#shoudl not affect anything?
modelProcessor.append(osim.ModOpAddReserves(ReserveActuatorStrength)) #how high
ClubfootModel = modelProcessor.process()
# #
inverse.setModel(modelProcessor)
inverse.setKinematics(osim.TableProcessor(Coordinates_Filename))
inverse.set_initial_time(start_time)
inverse.set_final_time(end_time)
inverse.set_mesh_interval(0.02)
inverse.set_kinematics_allow_extra_columns(True)

study = inverse.initialize()

MuscleNames = list(['tib_post_r', 'flex_hal_r','flex_dig_r','soleus_r','med_gas_r','lat_gas_r','ext_dig_r', 'per_brev_r',
'per_long_r', 'per_tert_r','ext_hal_r','tib_ant_r',
'tib_post_l', 'flex_hal_l','flex_dig_l','soleus_l','med_gas_l','lat_gas_l','ext_dig_l', 'per_brev_l',
'per_long_l', 'per_tert_l','ext_hal_l','tib_ant_l'])

# Ligaments = list(['TibioCalcaneal_r','TibioNavicular_r','PlantarFascia1Ph_r','PlantarFascia2Ph_r','PlantarFascia3Ph_r','LongPlantar1_r',
'LongPlantar2_r','LongPlantar3_r','LongPlantar4_r',
'CalcaneoNavicularPlantar1Mus_r','CalcaneoNavicularPlantar2Mus_r','CalcaneoNavicularPlantar3Mus_r',
'TarsoMetatarsalPlantar1Mus_r','TarsoMetatarsalPlantar2Mus_r','TarsoMetatarsalPlantar3Mus_r',
'CalcaneoCuboidPlantar1Mus_r', 'CalcaneoCuboidPlantar2Mus_r',

'TibioCalcaneal_l','TibioNavicular_l','PlantarFascia1Ph_l','PlantarFascia2Ph_l','PlantarFascia3Ph_l','LongPlantar1_l',
'LongPlantar2_l','LongPlantar3_l','LongPlantar4_l',
'CalcaneoNavicularPlantar1Mus_l','CalcaneoNavicularPlantar2Mus_l','CalcaneoNavicularPlantar3Mus_l',
'TarsoMetatarsalPlantar1Mus_l','TarsoMetatarsalPlantar2Mus_l','TarsoMetatarsalPlantar3Mus_l',
'CalcaneoCuboidPlantar1Mus_l', 'CalcaneoCuboidPlantar2Mus_l',

'CalcaneoFibular_r','TaloNavicularDorsal1Mus_r','TaloNavicularDorsal2Mus_r','CalcaneoCuboidBifurcateMus_r',
'CalcaneoNavicularBifurcateMus_r','TarsoMetatarsalDorsal1Mus_r','TarsoMetatarsalDorsal3Mus_r',
'TarsoMetatarsalDorsal5Mus_r','TarsoMetatarsalDorsal7Mus_r',

'CalcaneoFibular_l','TaloNavicularDorsal1Mus_l','TaloNavicularDorsal2Mus_l','CalcaneoCuboidBifurcateMus_l',
'CalcaneoNavicularBifurcateMus_l','TarsoMetatarsalDorsal1Mus_l','TarsoMetatarsalDorsal3Mus_l',
'TarsoMetatarsalDorsal5Mus_l','TarsoMetatarsalDorsal7Mus_l'])

# Parameter optimization
problem = study.updProblem()
problem.setModel(ClubfootModel)
Muscles = ClubfootModel.updMuscles()
Forces = ClubfootModel.updForceSet()
for cntrMuscle in range(Muscles.getSize()):
muscle = Muscles.get(cntrMuscle)
if str(muscle.getName()) in MuscleNames:

Min_TendonSlackLength = muscle.getTendonSlackLength()*0.5
Max_TendonSlackLength = muscle.getTendonSlackLength()*1.5
parameter=osim.MocoParameter(str(muscle.getName())+'_tendon_slack_length', ['/forceset/'+str(muscle.getName())], 'tendon_slack_length',osim.MocoBounds(Min_TendonSlackLength,Max_TendonSlackLength))
#problem.addParameter(parameter)

Min_OptimLength = muscle.getOptimalFiberLength()*0.5
Max_OptimLength = muscle.getOptimalFiberLength()*1.5
parameter=osim.MocoParameter(str(muscle.getName())+'_optimal_fibre_length', ['/forceset/'+str(muscle.getName())], 'optimal_fiber_length',osim.MocoBounds(Min_OptimLength,Max_OptimLength))
problem.addParameter(parameter)

# optimize contraction velocity?
Min_MaxContractionVelocity = muscle.getMaxContractionVelocity()*0.8
Max_MaxContractionVelocity = muscle.getMaxContractionVelocity()*2.5
parameter=osim.MocoParameter(str(muscle.getName())+'_max_contraction_velocity', ['/forceset/'+str(muscle.getName())], 'max_contraction_velocity',osim.MocoBounds(Min_MaxContractionVelocity,Max_MaxContractionVelocity))
#problem.addParameter(parameter)

# pennation angle
Min_PennationAngle = muscle.getPennationAngleAtOptimalFiberLength()*0.5
Max_PennationAngle = muscle.getPennationAngleAtOptimalFiberLength()*1.5
parameter=osim.MocoParameter(str(muscle.getName())+'_pennation_angle', ['/forceset/'+str(muscle.getName())], 'pennation_angle_at_optimal',osim.MocoBounds(Min_PennationAngle,Max_PennationAngle))
#problem.addParameter(parameter)

# # Ligaments
# for cntrForce in range(Forces.getSize()):
# force = Forces.get(cntrForce)
# Ligament_Force = osim.Ligament.safeDownCast(force)
# if str(force.getName()) in Ligaments:
# #Resting Length
# Min_RestingLength = Ligament_Force.getRestingLength()*0.6
# Max_RestingLength = Ligament_Force.getRestingLength()*1.4
# parameter_resting_length=osim.MocoParameter(str(force.getName())+'_resting_length', ['/forceset/'+str(force.getName())], 'resting_length',osim.MocoBounds(Min_RestingLength,Max_RestingLength))
# problem.addParameter(parameter_resting_length)

# # PCSA force
# Min_PCSA_force = Ligament_Force.get_pcsa_force()*0.6
# Max_PCSA_force = Ligament_Force.get_pcsa_force()*1.4
# parameter_pcsa_force=osim.MocoParameter(str(force.getName())+'_pcsa_force', ['/forceset/'+str(force.getName())], 'pcsa_force',osim.MocoBounds(Min_PCSA_force,Max_PCSA_force))
# problem.addParameter(parameter_pcsa_force)



#%% 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(20)
solver.set_optim_convergence_tolerance(.01)#.....10^-1
solver.set_optim_constraint_tolerance(1e-3)#1e-4 is default
solver.set_optim_max_iterations(4000)
solver.set_parameters_require_initsystem(False)
solution = study.solve()
solution.unseal()
solution.write('Clubfoot_MocoInverse_ParameterOptimization_Okt202109.sto')

Hi Christian,

What does Moco output to the console when you run the problem? It should print a list of the states in the problem, and none of the skeletal kinematic states should not appear in that list.

-Nick

I am not sure whether I understand your question correctly.....Moco is not printing anything in the console except for "Warning: intermediate_callback is disfunctional ......" . There is no list of the states.

Christian