Ground Reaction Forces after MocoContactTracking
Posted: Wed Jun 28, 2023 12:26 am
Dear OpenSim Moco Community,
I am trying to perform a MocoTracking study with Contact Tracking. However, the newly generated ground reaction forces are very different from the measured GRF signals (see file attached).
Changing tracking weights between the coordinate tracking, effort and grf tracking cost terms does only marginally affect the results. My contact model consists of osim.SmoothSphereHalfSpaceForce (see also attached file). I am struggling to figure out what the reason for the poor GRFs is. Below is my code used to generate the Tracking solution of one gait cycle.
Thanks for your help and kind regards
Christian
I am trying to perform a MocoTracking study with Contact Tracking. However, the newly generated ground reaction forces are very different from the measured GRF signals (see file attached).
Changing tracking weights between the coordinate tracking, effort and grf tracking cost terms does only marginally affect the results. My contact model consists of osim.SmoothSphereHalfSpaceForce (see also attached file). I am struggling to figure out what the reason for the poor GRFs is. Below is my code used to generate the Tracking solution of one gait cycle.
Thanks for your help and kind regards
Christian
Code: Select all
GRF_tracking = 1
CoordTrack = osim.MocoTrack()
CoordTrack.setName("PiG_CoordTracking")
# Construct a ModelProcessor and set it on the tool. The default
# muscles in the model are replaced with optimization-friendly
# DeGrooteFregly2016Muscles, and adjustments are made to the default muscle
# parameters.
# del modelProcessor
modelProcessor = osim.ModelProcessor(Scaled_modelFileName)
# modelProcessor.append(osim.ModOpAddExternalLoads(GRF_SourceFilename))
# modelProcessor.append(osim.ModOpIgnoreTendonCompliance())
# modelProcessor.append(osim.ModOpReplaceMusclesWithDeGrooteFregly2016())
# modelProcessor.append(osim.ModOpIgnorePassiveFiberForcesDGF())
# modelProcessor.append(osim.ModOpScaleActiveFiberForceCurveWidthDGF(1.5))
modelProcessor.append(osim.ModOpRemoveMuscles())
modelProcessor.append(osim.ModOpAddReserves(ReserveActuatorStrength))
CoordTrack.setModel(modelProcessor)
# Construct a TableProcessor of the coordinate data and pass it to the
# tracking tool. TableProcessors can be used in the same way as
# ModelProcessors by appending TableOperators to modify the base table.
# A TableProcessor with no operators, as we have here, simply returns the
# base table.
tableProcessor = osim.TableProcessor(CoordsFromMarkerTrack_Filename)
# tableProcessor.append(osim.TabOpConvertDegreesToRadians())
CoordTrack.setStatesReference(tableProcessor)
CoordTrack.set_states_global_tracking_weight(CoordTrackingWeight);
# This setting allows extra data columns contained in the states
# reference that don't correspond to model coordinates.
CoordTrack.set_allow_unused_references(True);
# Since there is only coordinate position data in the states references, this
# setting is enabled to fill in the missing coordinate speed data using
# the derivative of splined position data.
CoordTrack.set_track_reference_position_derivatives(True);
# Initial time, final time, and mesh interval.
CoordTrack.set_initial_time(start_time)
CoordTrack.set_final_time(end_time)
# CoordTrack.set_mesh_interval(0.01)
# Instead of calling solve(), call initialize() to receive a pre-configured
# MocoStudy object based on the settings above. Use this to customize the
# problem beyond the MocoTrack interface.
study = CoordTrack.initialize()
problem = study.updProblem()
if GRF_tracking == 1:
#% GRF Tracking
# % Track the right and left vertical and fore-aft ground reaction forces.
contactTracking = osim.MocoContactTrackingGoal('contact', GRFTrackingWeight)
contactTracking.setExternalLoadsFile(GRF_SourceFilename)
forceNamesRightFoot = osim.StdVectorString()
forceNamesRightFoot.append('R_PostHeel_Force')
forceNamesRightFoot.append('R_MedHeel_Force')
forceNamesRightFoot.append('R_LatHeel_Force')
forceNamesRightFoot.append('R_ProxMet1_Force')
forceNamesRightFoot.append('R_ProxMet3_Force')
forceNamesRightFoot.append('R_ProxMet5_Force')
forceNamesRightFoot.append('R_DistMet1_Force')
forceNamesRightFoot.append('R_DistMet3_Force')
forceNamesRightFoot.append('R_DistMet5_Force')
forceNamesRightFoot.append('R_Hallux_Force')
forceNamesRightFoot.append('R_Toes_Force')
#Left
forceNamesLeftFoot = osim.StdVectorString()
forceNamesLeftFoot.append('L_PostHeel_Force')
forceNamesLeftFoot.append('L_MedHeel_Force')
forceNamesLeftFoot.append('L_LatHeel_Force')
forceNamesLeftFoot.append('L_ProxMet1_Force')
forceNamesLeftFoot.append('L_ProxMet3_Force')
forceNamesLeftFoot.append('L_ProxMet5_Force')
forceNamesLeftFoot.append('L_DistMet1_Force')
forceNamesLeftFoot.append('L_DistMet3_Force')
forceNamesLeftFoot.append('L_DistMet5_Force')
forceNamesLeftFoot.append('L_Hallux_Force')
forceNamesLeftFoot.append('L_Toes_Force')
contactTracking.addContactGroup(forceNamesRightFoot, 'externalforce')
contactTracking.addContactGroup(forceNamesLeftFoot, 'externalforce_0')
contactTracking.setProjection('plane')
contactTracking.setProjectionVector(osim.Vec3(0, 0, 1))
problem.addGoal(contactTracking)
# Get a reference to the MocoControlGoal that is added to every MocoTrack
# problem by default.
effort = osim.MocoControlGoal.safeDownCast(problem.updGoal("control_effort"))
effort.setWeight(EffortWeight)
# Put a large weight on the pelvis CoordinateActuators, which act as the
# residual, or 'hand-of-god', forces which we would like to keep as small
# as possible.
model = modelProcessor.process()
model.initSystem()
forceSet = model.getForceSet()
for i in range(forceSet.getSize()-1):
# print(i)
forcePath = forceSet.get(i).getAbsolutePathString()
if 'pelvis' in str(forcePath):
effort.setWeightForControl(forcePath, 2)
#% Define solver and solve.
solver = osim.MocoCasADiSolver.safeDownCast(study.updSolver())
solver.resetProblem(problem)
solver.set_multibody_dynamics_mode('implicit')
solver.set_minimize_implicit_auxiliary_derivatives(True)
waux = 0.0001/forceSet.getSize()
solver.set_implicit_auxiliary_derivatives_weight(waux)
solver.set_num_mesh_intervals(100)
solver.set_optim_convergence_tolerance(1e-3)#.....10^-1
solver.set_optim_constraint_tolerance(1e-3)#1e-4 is default
solver.set_optim_max_iterations(1000)
solver.set_parameters_require_initsystem(False)
guess = solver.createGuess()
guess.randomizeAdd()
guess = osim.MocoTrajectory(Output_FilenameMocoCoordTrack_Guess)
# guess.generateAccelerationsFromValues()
solver.setGuess(guess)
#solve
CoordSolution = study.solve()
CoordSolution.unseal()
CoordSolution.write(Output_FilenameMocoCoordTrack)
#% extract GRF data
contact_r = osim.StdVectorString()
contact_l = osim.StdVectorString()
contact_r.append('R_PostHeel_Force')
# contact_r.add('contactFront_r')
contact_l.append('L_PostHeel_Force')
# contact_l.add('contactFront_l')
externalForcesTableFlat = osim.createExternalLoadsTableForGait(modelProcessor.process(), CoordSolution,contact_r,contact_l)
osim.STOFileAdapter.write(externalForcesTableFlat, 'gaitTracking_solutionGRF.sto');