This is an example that uses the MocoInverse tool and EMG data to create an EMG-driven simulation of walking.
function exampleEMGTracking_answers
clear; close all; clc;
%% Part 0: Load the
OpenSim and Moco libraries.
import org.opensim.modeling.*;
%% Part 1: Muscle redundancy problem: effort minimization.
% Solve the muscle redundancy problem while minimizing muscle excitations
% squared using the MocoInverse tool.
% Part 1a: Load a 19 degree-of-freedom model with 18 lower-limb,
% sagittal-plane muscles and a torque-actuated torso. This includes a set of
% ground reaction forces applied to the model via ExternalLoads, which is
% necessary for the muscle redundancy problem. See the function definition
% at the bottom of this file to see how the model is loaded and constructed.
model = getWalkingModel();
% Part 1b: Create the MocoInverse tool and set the Model.
inverse = MocoInverse();
inverse.setModel(model);
% Part 1c: Create a TableProcessor using the coordinates file from inverse
% kinematics.
coordinates = TableProcessor('coordinates.mot');
coordinates.append(TabOpLowPassFilter(6));
coordinates.append(TabOpUseAbsoluteStateNames());
% Part 1d: Set the kinematics reference for MocoInverse using the
% TableProcessor we just created.
inverse.setKinematics(coordinates);
inverse.set_kinematics_allow_extra_columns(true);
% Part 1e: Provide the solver settings: initial and final time, the mesh
% interval, and the constraint and convergence tolerances.
inverse.set_initial_time(0.83);
inverse.set_final_time(2.0);
inverse.set_mesh_interval(0.04);
inverse.set_constraint_tolerance(1e-3);
inverse.set_convergence_tolerance(1e-3);
if ~exist('effortSolution.sto', 'file')
% Part 1f: Solve the problem!
inverseSolution = inverse.solve();
solution = inverseSolution.getMocoSolution();
solution.write('effortSolution.sto');
end
%% Part 2: Plot the muscle redundancy problem solution.
% Load the experimental electromyography data and compare
% the effort minimization solution against this data. We will also use
% it later for the EMG-tracking problem. Each column in emg.sto is
% normalized so the maximum value for each signal is 1.0.
compareSolutionToEMG(emgReference, 'effortSolution.sto');
%% Part 3: Muscle redundancy problem: EMG-tracking.
% Modify the existing problem we created with the MocoInverse tool to solve
% a new problem where we will track electromyography (EMG) data.
% Part 3a: Call initialize() to get access to the MocoStudy contained within
% the MocoInverse instance. This will allow us to make additional
% modifications to the problem not provided by MocoInverse.
study = inverse.initialize();
problem = study.updProblem();
% Part 3b: Create a MocoControlTrackingGoal, set its weight, and provide
% the EMG data as the tracking reference. We also need to specify the
% reference labels for the four muscles whose EMG we will track.
tracking = MocoControlTrackingGoal('emg_tracking');
tracking.setWeight(5);
tracking.setReference(TableProcessor(emgReference));
tracking.setReferenceLabel('/forceset/gasmed_l', 'gastrocnemius');
tracking.setReferenceLabel('/forceset/tibant_l', 'tibialis_anterior');
tracking.setReferenceLabel('/forceset/bfsh_l', 'biceps_femoris');
tracking.setReferenceLabel('/forceset/glmax2_l', 'gluteus');
% Part 3c: The EMG signals in the tracking are all normalized to have
% a maximum value of 1, but the magnitudes of the excitations from the
% effort minimization solution suggest that these signals should be
% rescaled. Use addScaleFactor() to add a MocoParameter to the problem that
% will scale the reference data for the muscles in the tracking cost.
tracking.addScaleFactor('gastroc_factor', '/forceset/gasmed_l', [0.01 1.0]);
tracking.addScaleFactor('tibant_factor', '/forceset/tibant_l', [0.01 1.0]);
tracking.addScaleFactor('bifem_factor', '/forceset/bfsh_l', [0.01 1.0]);
tracking.addScaleFactor('gluteus_factor', '/forceset/glmax2_l', [0.01 1.0]);
% Part 3d: Add the tracking goal to the problem.
problem.addGoal(tracking)
% Part 3e: Update the MocoCasADiSolver with the updated MocoProblem using
% resetProblem().
solver = MocoCasADiSolver.safeDownCast(study.updSolver());
solver.resetProblem(problem);
% Part 3f: Tell MocoCasADiSolver that the MocoParameters we added to the
% problem via addScaleFactor() above do not require initSystem() calls on
% the model. This provides a large speed-up.
solver.set_parameters_require_initsystem(false);
if ~exist('trackingSolution.sto', 'file')
% Part 3g: Solve the problem!
solution = study.solve();
solution.write('trackingSolution.sto');
end
% Part 3h: Get the values of the optimized scale factors.
trackingSolution = MocoTrajectory('trackingSolution.sto');
gastroc_factor = trackingSolution.getParameter('gastroc_factor');
tibant_factor = trackingSolution.getParameter('tibant_factor');
bifem_factor = trackingSolution.getParameter('bifem_factor');
gluteus_factor = trackingSolution.getParameter('gluteus_factor');
%% Part 4: Plot the EMG-tracking muscle redundancy problem solution.
% Part 4a: Print the scale factor values to the command window.
fprintf('\n')
fprintf('Optimized scale factor values: \n')
fprintf('------------------------------ \n')
fprintf(['gastrocnemius = ' num2str(gastroc_factor) '\n'])
fprintf(['tibialis anterior = ' num2str(tibant_factor) '\n'])
fprintf(['biceps femoris short head = ' num2str(bifem_factor) '\n'])
fprintf(['gluteus = ' num2str(gluteus_factor) '\n'])
fprintf('\n')
% Part 4b: Re-scale the reference data using the optimized scale factors.
gastroc = emgReference.updDependentColumn('gastrocnemius');
tibant = emgReference.updDependentColumn('tibialis_anterior');
bifem = emgReference.updDependentColumn('biceps_femoris');
gluteus = emgReference.updDependentColumn('gluteus');
for t = 0:emgReference.getNumRows() - 1
gastroc.set(t, gastroc_factor * gastroc.get(t));
tibant.set(t, tibant_factor * tibant.get(t));
bifem.set(t, bifem_factor * bifem.get(t));
gluteus.set(t, gluteus_factor * gluteus.get(t));
end
% Part 4c: Generate the plots. Compare results to the effort minimization
% solution.
compareSolutionToEMG(emgReference, 'effortSolution.sto', ...
'trackingSolution.sto');
end
% Add a CoordinateActuator to the provided model with a specified optimal
% force.
function addCoordinateActuator(model, coordinateName, optForce)
import org.opensim.modeling.*;
coordSet = model.getCoordinateSet();
actu = CoordinateActuator();
actu.setName(['torque_' coordinateName]);
actu.setCoordinate(coordSet.get(coordinateName));
actu.setOptimalForce(optForce);
actu.setMinControl(-1);
actu.setMaxControl(1);
% Add to ForceSet
model.addForce(actu);
end
function [modelProcessor] = getWalkingModel()
import org.opensim.modeling.*;
% Load the 19 DOF, 18 muscle model from file.
model = Model('subject_walk_armless_18musc.osim');
% Add actuators representing the pelvis residual actuators and lumbar
% torques. These actuators are only as strong as they need to be for the
% problem to converge.
addCoordinateActuator(model, 'pelvis_tx', 60)
addCoordinateActuator(model, 'pelvis_ty', 300)
addCoordinateActuator(model, 'pelvis_tz', 35)
addCoordinateActuator(model, 'pelvis_tilt', 60)
addCoordinateActuator(model, 'pelvis_list', 35)
addCoordinateActuator(model, 'pelvis_rotation', 25)
addCoordinateActuator(model, 'lumbar_bending', 40)
addCoordinateActuator(model, 'lumbar_extension', 40)
addCoordinateActuator(model, 'lumbar_rotation', 25)
% We need additional actuators for hip rotation and hip adduction since the
% existing muscle act primarily in the sagittal plane.
addCoordinateActuator(model, 'hip_rotation_r', 100)
addCoordinateActuator(model, 'hip_rotation_l', 100)
addCoordinateActuator(model, 'hip_adduction_r', 100)
addCoordinateActuator(model, 'hip_adduction_l', 100)
% Create a ModelProcessor to make additional modifications to the model.
modelProcessor = ModelProcessor(model);
% Weld the subtalar and toe joints.
jointsToWeld = StdVectorString();
jointsToWeld.add('subtalar_r');
jointsToWeld.add('subtalar_l');
jointsToWeld.add('mtp_r');
jointsToWeld.add('mtp_l');
modelProcessor.append(ModOpReplaceJointsWithWelds(jointsToWeld));
% Apply the ground reaction forces to the model.
modelProcessor.append(ModOpAddExternalLoads('external_loads.xml'));
% Update the muscles: ignore tendon compliance and passive forces, replace
% muscle types with the DeGrooteFregly2016Muscle type, and scale the width
% of the active force length curve.
modelProcessor.append(ModOpIgnoreTendonCompliance());
modelProcessor.append(ModOpReplaceMusclesWithDeGrooteFregly2016());
modelProcessor.append(ModOpIgnorePassiveFiberForcesDGF());
modelProcessor.append(ModOpScaleActiveFiberForceCurveWidthDGF(1.5));
% Add a set a weak reserves to the sagittal-plane joints in the model.
modelProcessor.append(ModOpAddReserves(1.0));
end
function compareSolutionToEMG(varargin)
import org.opensim.modeling.*;
% Retrieve the inputs.
emgReference = varargin{1};
effortSolution = MocoTrajectory(varargin{2});
if nargin > 2
trackingSolution = MocoTrajectory(varargin{3});
end
% Create a time vector for the EMG data that is consistent with the problem
% time range.
time = effortSolution.getTimeMat();
startIndex = emgReference.getNearestRowIndexForTime(time(1)) + 1;
endIndex = emgReference.getNearestRowIndexForTime(time(end)) + 1;
emgTime = linspace(time(1), time(end), endIndex - startIndex + 1);
% Plot results from the muscles in the left leg.
titles = {'soleus', 'gastroc.', 'tib. ant.', 'hamstrings', ...
'bi. fem. sh.', 'vastus', 'rec. fem.', 'gluteus', 'psoas'};
emgSignals = {'soleus', 'gastrocnemius', 'tibialis_anterior', ...
'hamstrings', 'biceps_femoris', 'vastus', ...
'rectus_femoris', 'gluteus', 'none'};
muscles = {'soleus_l', 'gasmed_l', 'tibant_l', 'semimem_l', 'bfsh_l', ...
'vasint_l', 'recfem_l', 'glmax2_l', 'psoas_l'};
figure;
for i = 1:length(titles)
subplot(3,3,i)
% If it exists, plot the EMG signal for this muscle.
if ~strcmp(emgSignals{i}, 'none')
col = emgReference.getDependentColumn(emgSignals{i}).getAsMat();
emg = col(startIndex:endIndex)';
patch([emgTime fliplr(emgTime)], [emg fliplr(zeros(size(emg)))], ...
[0.7 0.7 0.7])
alpha(0.5);
end
hold on
% Plot the control from the muscle effort solution.
control = effortSolution.getControlMat(['/forceset/' muscles{i}]);
plot(time, control, '-k', 'linewidth', 2)
% If provided, plot the control from the tracking solution.
if nargin > 2
control = trackingSolution.getControlMat(['/forceset/' muscles{i}]);
plot(time, control, '-r', 'linewidth', 2)
end
% Plot formatting.
title(titles{i})
xlim([0.83 2.0])
ylim([0 1])
if i == 1
if nargin > 2
legend('EMG', 'effort', 'tracking', 'location', 'best')
else
legend('EMG', 'effort', 'location', 'best')
end
end
if i > 6
xlabel('time (s)')
else
xticklabels([])
end
if ~mod(i-1, 3)
ylabel('excitation')
else
yticklabels([])
end
end
end
1.8.14