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
Parameter optimization verification
- Christian Greve
- Posts: 41
- Joined: Mon Jun 13, 2016 11:14 pm
- Nicholas Bianco
- Posts: 1003
- Joined: Thu Oct 04, 2012 8:09 pm
Re: Parameter optimization verification
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?
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
- Posts: 41
- Joined: Mon Jun 13, 2016 11:14 pm
Re: Parameter optimization verification
@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
Would I get more realistic results if I enable tendon compliance during the Parameter Optimization?
Thanks
Christian
- Attachments
-
- MuscleAnalysis_ParamsOptimized.png (83.68 KiB) Viewed 326 times
- Nicholas Bianco
- Posts: 1003
- Joined: Thu Oct 04, 2012 8:09 pm
Re: Parameter optimization verification
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
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
- Posts: 41
- Joined: Mon Jun 13, 2016 11:14 pm
Re: Parameter optimization verification
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
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
- Attachments
-
- ParamOptimization_EMGEffort.PNG (145.57 KiB) Viewed 303 times
- Nicholas Bianco
- Posts: 1003
- Joined: Thu Oct 04, 2012 8:09 pm
Re: Parameter optimization verification
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
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