Parameter optimization verification

[Christian Greve]

I think I found the bug and I also understand your question now....

I removed the line:
problem.setModel(ClubfootModel)

now the output only contains the new optimal_fibre_length.....

Thanks

Christian

[Nicholas Bianco]

Ah yes, that should be it!

You can't reset the model here because the prescribed motion constraint that MocoInverse implements is applied to the model directly.

Do your results look more reasonable now?

[Christian Greve]

@reasonable results: the proposed changes in muscle parameters look reasonable but this is still hard to verify I think. When I run muscle analysis with the adjusted model and plot the normalized fibre length it does still look odd.....I attached a .png of the muscle analysis output. This is after optimizing tendon slack length, optimal muscle fibre length and max contraction velocity.

Would I get more realistic results if I enable tendon compliance during the Parameter Optimization?

Thanks

Christian

[Nicholas Bianco]

Hi Christian,

Some of the normalized fiber lengths here seem very low, but others seem reasonable (i.e., near 1).

What are you optimizing for in this problem? Are you tracking EMG data? The problem might choose optimal fiber lengths near the extreme ends of the bounds unless you design a cost function that directly influences these parameters.

-Nick

[Christian Greve]

I use EMG tracking and the effort goal. When optimizing on both optimal_fibre_length and tendon slack I get the attached results on normalized muscle fibre length. Please see also the code for clarification.

Thanks

Christian

Code: Select all

def solveMocoInverse_AND_OptimizeParameter(Scaled_modelFileName,GRF_SourceFilename, Coordinates_Filename, Output_Filename, ReserveActuatorStrength, start_time, end_time):
    ''' 28_01_2022: worked for opt_fibre_length right side--> try tendon-slack and left side'''

    # 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() 
    
    
    #%% effort goal
    controlEffortWeight = 10
    problem = study.updProblem()
    effort = osim.MocoControlGoal("control_effort") #.safeDownCast(problem.updGoal("control_effort"))
    effort.setName('effort')
    effort.setWeight(controlEffortWeight)
    problem.addGoal(effort)
    
    #%% EMG Tracking Goal
    problem = study.updProblem()
    # Add electromyography tracking.
    emgTracking = osim.MocoControlTrackingGoal('emg_tracking')
    emgTracking.setWeight(10.0)
    # Each column in electromyography.sto is normalized so the maximum value in
    # each column is 1.0.
    controlsRef = osim.TimeSeriesTable(EMG_Filename)
    emgTracking.setReference(osim.TableProcessor(controlsRef))
    # Associate actuators in the model with columns in electromyography.sto.
    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
    problem.addGoal(emgTracking)
    
    #%% Muscle and Lig names
    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:
            # 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)
            
    # 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)



        #%% adjust model and add ligament optimization (as in MocoStudy)
        
        
    
    #%% 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-4)#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')
    
    return solution

[Nicholas Bianco]

Hi Christian,

I scanned the code briefly, but I didn't see anything abnormal. It might be that the ranges on the muscle parameters are too wide, and the optimization is finding a local minimum that leads to short fiber lengths.

You might also try only optimizing one parameter at a time, to limit their interaction during the optimization and get a better sense of what matters to your EMG optimization.

-Nick