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.
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