Ground Reaction Forces after MocoContactTracking

[Christian Greve]

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

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');

[Aaron Fox]

Hi Christian,

I think the size of your contact spheres might be a problem here that's limiting finding an appropriate contact tracking solution. From what I've seen the contact spheres produce a force when they intersect the contact half space plane, and if they are not in contact (e.g. completely above or below the plane) then they don't produce any force. The size of the spheres you have therefore have a minimal window of contact where they can produce force, and hence it might be hard to find a solution that satisfies the kinematics you're tracking while generating the GRF curves.

Most examples I've seen have had contact spheres with a radius between 0.02-0.04, I haven't seen any examples with contact spheres this small before.

Aaron

[Christian Greve]

Dear Aaron,

thanks for your reply. However, the size of the contact spheres does not seem to have much of an effect. When increasing the size to 0.04 the GRFs are still nothing close to measured GRFs. I figured that the hertzSmoothing and huntCrossleySmoothing affect the GRF tracking a lot. But also here I cannot generate good solutions.

THis is the list of params I use:
staticFriction = 0.8;
dynamicFriction = 0.8;
viscousFriction = 0.5;
transitionVelocity = 0.2;
hertzSmoothing = 300
huntCrossleySmoothing = 50

All force spheres have the calcaneus as frame. That is also where the GRFs applie.

And the tracking weights are distributed as follows:
GRFTrackingWeight = 10
CoordTrackingWeight = 10
EffortWeight = 1

Thanks

Christian

[Nicholas Bianco]

Hi Christian,

I usually recommend adding residual forces to the model (i.e., CoordinateActuators to the pelvis) in order to generate a solution that tracks the pelvis well. This solution can then be used an initial guess for subsequent problems where you can then reduce the residual force strengths or remove them altogether to get a realistic solution.

-Nick

[Christian Greve]

Thanks Nick,

I do add reserve actuators with the command: modelProcessor.append(osim.ModOpAddReserves(ReserveActuatorStrength))

The reserve actuator strength is set to 250. Should I add in addition coord actuators to the pelvis?

I also changed the effort weight for the pelvis forces in this way:

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)

Regards

Christian

[Ross Miller]

Hi Christian,

A couple things to try/consider:

(1) Try turning the tracking weight on the GRF way high and seeing if the solution then has good tracking of them. This will most likely come at the cost of other terms in the cost function, but it can help with identifying the issue.

(2) From the image (and from the very-low vertical GRF in the simulation result), it looked like your model might only include the lower limbs and pelvis. Is there an upper body? These kinds of problems need a model of the full body, or at least the full mass of the body. An easy way to compute what the mass of your model should be is to take the target GRF over a full stride, average them in the vertical direction, and that average should be the model's weight. This won't work out exactly right using data for a single stride that isn't necessarily periodic but it should be close enough for tracking purposes.

If your model has no upper body and is missing all that mass, a quick way to include the missing mass is to figure out the total body mass from the GRF method above, then increase the mass of the model's pelvis so that the model's mass equals the GRF-based mass.

Hope this helps,
Ross

[Christian Greve]

Thanks for the suggestions Ross,

The GRF tracking weight was set to 10000 with the Coordinate and effort costs being set to 1 and the weight of the torso of the model is centerd at the pelvis. The total model weight corresponds with the subjects weight.

Still the generated GRFs do not correspond to measured GRFs.

Regards

Christian

[Nicholas Bianco]

Hi Christian,

I would suggest adding CoordinateActuators to the pelvis with strengths independent of the reserve actuators (i.e., much stronger). I would also then not penalize these pelvis actuators as much in the control cost term. Basically, make the model as strong as necessary to track the ground reaction force data tightly while still roughly tracking the coordinate trajectories.

If you find that the model can't track the ground reaction forces as well when you decrease the pelvis actuator strength (or penalize them heavier in the cost function), then it might be that there are dynamic inconsistencies between the model, coordinate trajectories, and ground reaction force data. In this case, you could allow Moco to modify the coordinate trajectories in the tracking problem such that you track the GRF data (in the case that you trust the GRF data more than the coordinate data).

Best,
Nick

[Christian Greve]

Dear community and Moco team,

I wanted to follow up on the Contact tracking task. I made some progress but still do not get the desired outcome. Attached is a .png and below my code.

The solution is not succeeding and restoration fails after 2000 or sometimes 5000 itertions. Wht I noticed is that when visualizing the solution, the GRF point of application remains in the middle of the projected floor (se .png). For now I exclude any muscles and replace them with actuators. When only tracking coordinates the solution converges.

The data is pretty clean from a healthy young adult and two force plates. Model was RRA adjusted.

Thanks for your help

Christian

Code: Select all

Scaled_modelFileName = os.path.join(Path, 'PiG_SphereModel_LargeRadiu_Scaled_RRA.osim')
# Scaled_modelFileName = os.path.join(Path,'Patient_SCALED.osim')
GRF_SourceFilename = os.path.join(Path,'Walk100_04_GRF.xml')
IK_Filename = os.path.join(Path,'MocoInverse.sto')#'#Walk_100_02__Kinematics_q.sto')#'MocoInverseWithEMG_Okt202108_Clubfoot_Final_NoMocoProblem.sto')
# TRC_FileName = os.path.join(Path,'Walk_100_02.trc')
# CoordsFromMarkerTrack_Filename = os.path.join(Path,'MocoTrack_Solution_Marker.sto')
# ModelPath = r'C:\Users\tacko\OneDrive - UMCG\UMCG\Onderzoek\OffRoad2021_Clubfoot\Data\Healthy\CluFO_01\OpenSim\Models'
# Scaled_modelFileName = os.path.join(ModelPath, 'Clubfoot_MuscleParamsOptimized_CluFo01.osim')#'MillardModel.osim')

ReserveActuatorStrength = 250000
Output_FilenameMocoCoordTrack = os.path.join(Path,'MocoTrack_Solution_Coord.sto')
FilenameMocoCoordTrack_Guess = os.path.join(Path,'PiG_CoordTracking_solution_GuessFile.sto')

# Output_Filename = os.path.join(Path,'StateTracking_Solution.sto')
# Guess_Filename = os.path.join(Path,'MocoTrack_Solution.sto')
start_time = 1.2#1.5
end_time = 2.2#2.5
GRFTrackingWeight = 50
CoordTrackingWeight = 10
EffortWeight = 0.1
MarkerTrackingWeight = 1
convergence_tol = 1e-3
constraint_tol = 1e-3



def addCoordinateActuator(model, coordName, optForce):
    coordSet = model.updCoordinateSet()
    actu = osim.CoordinateActuator()
    actu.setName('tau_' + coordName)
    actu.setCoordinate(coordSet.get(coordName))
    actu.setOptimalForce(optForce)
    actu.setMinControl(-1)
    actu.setMaxControl(1)
    model.addComponent(actu)


#%% Build Coordinate Tracking Task
# Create and name an instance of the MocoTrack tool.
GRF_tracking = 1
CoordTrack = osim.MocoTrack()
CoordTrack.setName("PiG_CoordTracking")
 
# Construct a ModelProcessor
modelProcessor = osim.ModelProcessor(Scaled_modelFileName)
#modelProcessor.append(osim.ModOpAddExternalLoads(GRF_SourceFilename))

#modelProcessor.append(osim.ModOpReplaceMusclesWithDeGrooteFregly2016())
#modelProcessor.append(osim.ModOpIgnorePassiveFiberForcesDGF())
#modelProcessor.append(osim.ModOpScaleActiveFiberForceCurveWidthDGF(1.5))
#modelProcessor.append(osim.ModOpIgnoreTendonCompliance())

modelProcessor.append(osim.ModOpRemoveMuscles())
modelProcessor.append(osim.ModOpAddReserves(ReserveActuatorStrength))
CoordTrack.setModel(modelProcessor)
 
# Construct a TableProcessor
tableProcessor = osim.TableProcessor(IK_Filename)
#tableProcessor.append(osim.TabOpConvertDegreesToRadians())
CoordTrack.setStatesReference(tableProcessor)
CoordTrack.set_states_global_tracking_weight(CoordTrackingWeight);
CoordTrack.set_allow_unused_references(True);
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.
effort = osim.MocoControlGoal("control")
effort.setWeight(EffortWeight)

model = modelProcessor.process()
model.initSystem()
forceSet = model.updForceSet()
# 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.
#for i in range(forceSet.getSize()-1):
    # print(i)
    #forcePath = forceSet.get(i).getAbsolutePathString()
    #if 'pelvis' in str(forcePath):
        #effort.setWeightForControl(forcePath, 20)
       

problem.addGoal(effort)

#/////////////////////////
#5 add coord actuators for the pelvis

#optForce = 20000
#CoordSet = model.updCoordinateSet()
#for cntrCoord in range(CoordSet.getSize()):
#    coord = CoordSet.get(cntrCoord)
#    coordName = coord.getName()
#    coord.setWeightForControl(coordName,2)
   
#    if 'pelvis' in coordName:
#        print(coordName)
#        actu = osim.CoordinateActuator()
#        actu.setName('tau_' + coordName)
#        actu.setCoordinate(CoordSet.get(coordName))
#        actu.setOptimalForce(optForce)
#        actu.setMinControl(-1)
#        actu.setMaxControl(1)
#        model.addComponent(actu)
                 
#model.finalizeConnections()

#///////////////////////////////////////////
 
#%% 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)
forceSet = model.getForceSet()
waux = 0.0001/forceSet.getSize()
solver.set_implicit_auxiliary_derivatives_weight(waux)
solver.set_num_mesh_intervals(100)
solver.set_optim_convergence_tolerance(convergence_tol)#.....10^-1
solver.set_optim_constraint_tolerance(constraint_tol)#1e-4 is default
solver.set_optim_max_iterations(50000)
solver.set_parameters_require_initsystem(False)

# create guess file
guess = solver.createGuess()
guess.randomizeAdd()
prevSolution = osim.MocoTrajectory(FilenameMocoCoordTrack_Guess)
StatesTableGuess = guess.exportToStatesTable()
columnLabels_to_keep = StatesTableGuess.getColumnLabels()
StatesTablePrevSol = prevSolution.exportToStatesTable()
columnLabels_all = StatesTablePrevSol.getColumnLabels()
for column in columnLabels_all:
    if not column in columnLabels_to_keep:
        StatesTablePrevSol.removeColumn(column)
guess.insertStatesTrajectory(StatesTablePrevSol, True);
ControlsTable = prevSolution.exportToControlsTable();
guess.insertControlsTrajectory(ControlsTable, True);

#set guess file
solver.setGuess(guess)
#solve
CoordSolution = study.solve()
CoordSolution.unseal()
CoordSolution.write(Output_FilenameMocoCoordTrack)

[Nicholas Bianco]

Hi Christian,

The GRFs at the origin are just how they are visualized, assuming you computed them using "createExternalLoadsTableForGait()". If you download the latest version of the API, there's a new addition to that utility that computes centers-of-pressure in addition to the GRFs, which will give you the correct visualization.

As for your convergence issues, I recommend simplifying the problem as much as you can to get a set of kinematics that tracks the GRFs, and then add problem complexity back in. For example, start with a torque-driven model with strong actuators at all degrees-of-freedom to get a valid trajectory for the kinematics which you can then use in a subsequent initial guess. Depending on the problem, this step could take a few iterations, but it's much better iterate with a simple model.

Also, I noticed you're using implicit multibody dynamics mode. Most users have found more consistent results when using explicit multibody dynamics, so I recommend you switch over to that.

Best,
Nick