Matlab version of example2DWalking with metabolic cost minimization

[Brian Umberger]

Hello All,

First, a big congratulations to Nick, Chris, and the whole team for the 1.0 release of Moco! It is a major contribution to the biomechanics community.

One of the great additions in the new release is a smooth muscle energetics model, contributed by Antoine Falisse. Antoine provided a C++ example for solving a 2-D walking problem with minimization of the metabolic cost of transport (see example2DWalkingMetabolics.cpp). Here is a Matlab implementation for anyone who is interested in cost minimization, but prefers to work in a scripting language.

Best regards,

Brian

Code: Select all

% -------------------------------------------------------------------------- %
% OpenSim Moco: example2DWalkingMetabolicCost.m                              %
% -------------------------------------------------------------------------- %
% Copyright (c) 2021 Stanford University and the Authors                     %
%                                                                            %
% Author(s): Brian Umberger                                                  %
%                                                                            %
% Licensed under the Apache License, Version 2.0 (the "License"); you may    %
% not use this file except in compliance with the License. You may obtain a  %
% copy of the License at http://www.apache.org/licenses/LICENSE-2.0          %
%                                                                            %
% Unless required by applicable law or agreed to in writing, software        %
% distributed under the License is distributed on an "AS IS" BASIS,          %
% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.   %
% See the License for the specific language governing permissions and        %
% limitations under the License.                                             %
% -------------------------------------------------------------------------- %

% This is a Matlab implementation of an example optimal control
% problem (2-D walking) orginally created in C++ by Antoine Falisse
% (see: example2DWalkingMetabolics.cpp).
%
% This example features a tracking simulation of walking that includes
% minimization of the metabolic cost of transport computed using a smooth
% approximation of the metabolic energy model of Bhargava et al (2004).
%
% The code is inspired from Falisse A, Serrancoli G, Dembia C, Gillis J,
% De Groote F: Algorithmic differentiation improves the computational
% efficiency of OpenSim-based trajectory optimization of human movement.
% PLOS One, 2019.
%
% Model
% -----
% The model described in the file '2D_gait.osim' included in this file is a
% modified version of the 'gait10dof18musc.osim' available within OpenSim. We
% replaced the moving knee flexion axis by a fixed flexion axis, replaced the
% Millard2012EquilibriumMuscles by DeGrooteFregly2016Muscles, and added
% SmoothSphereHalfSpaceForces (two contacts per foot) to model the
% contact interactions between the feet and the ground.
%
% Data
% ----
% The coordinate data included in the 'referenceCoordinates.sto' comes from
% predictive simulations generated in Falisse et al. 2019.

clear;

% Load the Moco libraries
import org.opensim.modeling.*;

% Set a coordinate tracking problem where the goal is to minimize the
% difference between provided and simulated coordinate values and speeds
% as well as to minimize an effort cost (squared controls) and a metabolic
% cost (metabolic energy normalized by distance traveled and body mass;
% the metabolics model is based on a smooth approximation of the
% phenomenological model described by Bhargava et al. (2004)). The provided
% data represents half a gait cycle. Endpoint constraints enforce periodicity
% of the coordinate values (except for pelvis tx) and speeds, coordinate
% actuator controls, and muscle activations.

track = MocoTrack();
track.setName('gaitTrackingMetCost');


% Define the optimal control problem
% ==================================
baseModel = Model('2D_gait.osim');


% Add metabolic cost model
metabolics = Bhargava2004SmoothedMuscleMetabolics();
metabolics.setName('metabolic_cost');
metabolics.set_use_smoothing(true);

% This next part can easily be put in a loop for models with more muscles
metabolics.addMuscle('hamstrings_r', Muscle.safeDownCast(baseModel.getComponent('hamstrings_r')));
metabolics.addMuscle('hamstrings_l', Muscle.safeDownCast(baseModel.getComponent('hamstrings_l')));
metabolics.addMuscle('bifemsh_r', Muscle.safeDownCast(baseModel.getComponent('bifemsh_r')));
metabolics.addMuscle('bifemsh_l', Muscle.safeDownCast(baseModel.getComponent('bifemsh_l')));
metabolics.addMuscle('glut_max_r', Muscle.safeDownCast(baseModel.getComponent('glut_max_r')));
metabolics.addMuscle('glut_max_l', Muscle.safeDownCast(baseModel.getComponent('glut_max_l')));
metabolics.addMuscle('iliopsoas_r', Muscle.safeDownCast(baseModel.getComponent('iliopsoas_r')));
metabolics.addMuscle('iliopsoas_l', Muscle.safeDownCast(baseModel.getComponent('iliopsoas_l')));
metabolics.addMuscle('rect_fem_r', Muscle.safeDownCast(baseModel.getComponent('rect_fem_r')));
metabolics.addMuscle('rect_fem_l', Muscle.safeDownCast(baseModel.getComponent('rect_fem_l')));
metabolics.addMuscle('vasti_r', Muscle.safeDownCast(baseModel.getComponent('vasti_r')));
metabolics.addMuscle('vasti_l', Muscle.safeDownCast(baseModel.getComponent('vasti_l')));
metabolics.addMuscle('gastroc_r', Muscle.safeDownCast(baseModel.getComponent('gastroc_r')));
metabolics.addMuscle('gastroc_l', Muscle.safeDownCast(baseModel.getComponent('gastroc_l')));
metabolics.addMuscle('soleus_r', Muscle.safeDownCast(baseModel.getComponent('soleus_r')));
metabolics.addMuscle('soleus_l', Muscle.safeDownCast(baseModel.getComponent('soleus_l')));
metabolics.addMuscle('tib_ant_r', Muscle.safeDownCast(baseModel.getComponent('tib_ant_r')));
metabolics.addMuscle('tib_ant_l', Muscle.safeDownCast(baseModel.getComponent('tib_ant_l')));

baseModel.addComponent(metabolics);
baseModel.finalizeConnections();


% Reference data for tracking problem
tableProcessor = TableProcessor('referenceCoordinates.sto');
tableProcessor.append(TabOpLowPassFilter(6));

% ModelProcessor modelprocessor = ModelProcessor(baseModel);
modelProcessor = ModelProcessor(baseModel);
track.setModel(modelProcessor);
track.setStatesReference(tableProcessor);
track.set_states_global_tracking_weight(30);
track.set_allow_unused_references(true);
track.set_track_reference_position_derivatives(true);
track.set_apply_tracked_states_to_guess(true);
track.set_initial_time(0.0);
track.set_final_time(0.47008941);
study = track.initialize();
problem = study.updProblem();


% Goals
% =====

% Symmetry (to permit simulating only one step)
symmetryGoal = MocoPeriodicityGoal('symmetryGoal');
problem.addGoal(symmetryGoal);
model = modelProcessor.process();
model.initSystem();

% Symmetric coordinate values (except for pelvis_tx) and speeds
for i = 1:model.getNumStateVariables()
    currentStateName = string(model.getStateVariableNames().getitem(i-1));
    if startsWith(currentStateName , '/jointset')
        if contains(currentStateName,'_r')
            pair = MocoPeriodicityGoalPair(currentStateName, ...
                           regexprep(currentStateName,'_r','_l'));
            symmetryGoal.addStatePair(pair);
        end
        if contains(currentStateName,'_l')
            pair = MocoPeriodicityGoalPair(currentStateName, ...
                           regexprep(currentStateName,'_l','_r'));
            symmetryGoal.addStatePair(pair);
        end
        if (~contains(currentStateName,'_r') && ...
            ~contains(currentStateName,'_l') && ...
            ~contains(currentStateName,'pelvis_tx/value') && ...
            ~contains(currentStateName,'/activation'))
            symmetryGoal.addStatePair(MocoPeriodicityGoalPair(currentStateName));
        end
    end
end

% Symmetric muscle activations
for i = 1:model.getNumStateVariables()
    currentStateName = string(model.getStateVariableNames().getitem(i-1));
    if endsWith(currentStateName,'/activation')
        if contains(currentStateName,'_r')
            pair = MocoPeriodicityGoalPair(currentStateName, ...
                         regexprep(currentStateName,'_r','_l'));
            symmetryGoal.addStatePair(pair);
        end
        if contains(currentStateName,'_l')
            pair = MocoPeriodicityGoalPair(currentStateName, ...
                         regexprep(currentStateName,'_l','_r'));
            symmetryGoal.addStatePair(pair);
        end
    end
end

% Symmetric coordinate actuator controls
symmetryGoal.addControlPair(MocoPeriodicityGoalPair('/lumbarAct'));

% Get a reference to the MocoControlGoal that is added to every MocoTrack
% problem by default and change the weight
effort = MocoControlGoal.safeDownCast(problem.updGoal('control_effort'));
effort.setWeight(0.1);


% Metabolic cost; total metabolic rate includes activation heat rate,
% maintenance heat rate, shortening heat rate, mechanical work rate, and
% basal metabolic rate.
metGoal = MocoOutputGoal('met',0.1);
problem.addGoal(metGoal);
metGoal.setOutputPath('/metabolic_cost|total_metabolic_rate');
metGoal.setDivideByDisplacement(true);
metGoal.setDivideByMass(true);


% Bounds
% ======
problem.setStateInfo('/jointset/groundPelvis/pelvis_tilt/value', [-20*pi/180, -10*pi/180]);
problem.setStateInfo('/jointset/groundPelvis/pelvis_tx/value', [0, 1]);
problem.setStateInfo('/jointset/groundPelvis/pelvis_ty/value', [0.75, 1.25]);
problem.setStateInfo('/jointset/hip_l/hip_flexion_l/value', [-10*pi/180, 60*pi/180]);
problem.setStateInfo('/jointset/hip_r/hip_flexion_r/value', [-10*pi/180, 60*pi/180]);
problem.setStateInfo('/jointset/knee_l/knee_angle_l/value', [-50*pi/180, 0]);
problem.setStateInfo('/jointset/knee_r/knee_angle_r/value', [-50*pi/180, 0]);
problem.setStateInfo('/jointset/ankle_l/ankle_angle_l/value', [-15*pi/180, 25*pi/180]);
problem.setStateInfo('/jointset/ankle_r/ankle_angle_r/value', [-15*pi/180, 25*pi/180]);
problem.setStateInfo('/jointset/lumbar/lumbar/value', [0, 20*pi/180]);


% Configure the solver.
% =====================
solver = MocoCasADiSolver.safeDownCast(study.updSolver());
solver.resetProblem(problem);
solver.set_num_mesh_intervals(50);
solver.set_verbosity(2);
solver.set_optim_solver("ipopt");
solver.set_optim_convergence_tolerance(1e-4);
solver.set_optim_constraint_tolerance(1e-4);
solver.set_optim_max_iterations(10000);


% Solve the problem
% =================
solution = study.solve();

% Create a full stride from the periodic single step solution.
% For details, view the Doxygen documentation for createPeriodicTrajectory().
fullStride = opensimMoco.createPeriodicTrajectory(solution);
fullStride.write('gaitTrackingMetCost_solution_fullcycle.sto');


% To report the COT we multiply the metabolic cost objective term by 10 because
% it had been scaled by 0.1
disp('   ')
disp( ['The metabolic cost of transport is: '  ...
      num2str(10*solution.getObjectiveTerm('met')) ...
      ' J/kg/m'])
disp('   ')

% Uncomment next line to visualize the result
study.visualize(fullStride);

[Nicholas Bianco]

Thanks Brian! This is fantastic.

We'll get this merged into the codebase soon.

-Nick

[Ross Miller]

I installed the new version just now and immediately came here, thinking "I wonder if someone has posted a Matlab example of the metabolic cost minimization?"

The new features/updates in 1.0 all look very nice, looking forward to using this.

Three questions:

(1) When this code finished running, it didn't print the cost function terms in the Command Window, and the output .sto file also doesn't have them. Is that expected? Version 0.4 used to do that.

(2) In the "total metabolic rate", what is the assumed basal rate (typically something like 1.0 W per kg body mass) and is it modifiable?

(3) I wasn't sure the math in the code below works out to give "metabolic cost" (J/m/kg) in the cost function:

Code: Select all

% Metabolic cost; total metabolic rate includes activation heat rate,
% maintenance heat rate, shortening heat rate, mechanical work rate, and
% basal metabolic rate.
metGoal = MocoOutputGoal('met',0.1);
problem.addGoal(metGoal);
metGoal.setOutputPath('/metabolic_cost|total_metabolic_rate');
metGoal.setDivideByDisplacement(true);
metGoal.setDivideByMass(true);

If the output is the "total metabolic rate" (W or J/s) and this is divided by displacement (m) and mass (kg), the output quantity would be J/m/s/kg. Is the output of metGoal being integrated over the movement time?

Ross

[Brian Umberger]

Hi Ross,

Thanks for the questions.

1) I think that is the new behavior in Moco 1.0.0. In my opinion, the output file behavior is an improvement, as it should only write that out if you want, and to the filename that you want (i.e. using solution.write('filename.sto')). The objective and its terms can easily be obtained using some combination of:

solution.getObjective
solution.getObjectiveTerm
solution.getObjectiveTermByIndex
solution.getObjectiveTermNames

2) It looks like the default basal metabolic rate is 1.2 W per kg body mass:
https://github.com/opensim-org/opensim- ... tabolics.h
That should be modifiable but I haven't tried changing it yet. Maybe Nick can confirm.

3) The "total_metabolic_rate" is in W (J/s), but that is not the goal. The rate of metabolic energy consumption is fed into the OutputGoal, which integrates over the duration of the simulation, resulting in J, which is then divided by mass and displacement.

Hope that helps.

Brian

[Brian Umberger]

Here is a follow-up.

If you would like to add a breakdown of the individual terms in the overall objective function to this Matlab code, similar to the output from earlier versions of Moco, here is one way to do it. Add the following code after solution = study.solve().

Code: Select all

% Breakdown of terms in objective function
disp('   ')
disp('Breakdown of objective (including weights):')
for i = 1:solution.getNumObjectiveTerms()
   termName = solution.getObjectiveTermNames().get(i-1);
   termNamestr = termName.toCharArray';
   disp(['  ' termNamestr ': ' num2str(solution.getObjectiveTermByIndex(i-1))])
end

[Nicholas Bianco]

When this code finished running, it didn't print the cost function terms in the Command Window, and the output .sto file also doesn't have them. Is that expected? Version 0.4 used to do that.

I can't remember if we disabled this feature or if it got lost somewhere when moving Moco into OpenSim. Either way, Brian's suggestions will work. You can also open the MocoSolution file and the cost values are in the metadata in the file header.

In the "total metabolic rate", what is the assumed basal rate (typically something like 1.0 W per kg body mass) and is it modifiable?

Yes, it's modifiable via 'set_basal_coefficient()'.

-Nick

[Burak Kula]

Hello Brian,
Firstly, I want to thank you for implementing "example2DWalkingMetabolics.cpp" into MATLAB script.

But I had some problems, when I run your code on MATLAB. Currently I am trying to minimize metabolic energy expenditure of upper body and it seems like your "example2DWalkingMetabolics.m" is quite useful for me.
When I run your script without any change. I had an error that says "Bhargava2004SmoothedMuscleMetabolics() cannot be recognized as a function or parameter".

Do you have any idea about how can I run it properly ?

Regards,
Burak KULA.

[Nicholas Bianco]

Hi Burak,

You're probably using an older version of Moco. We renamed Bhargava2004Metabolics to Bhargava2004SmoothedMuscleMetabolics in the Moco version that comes with OpenSim 4.2.

-Nick

[Pasha van Bijlert]

Hi all,

I've got some questions related to this tutorial & Bhargava2004SmoothedMuscleMetabolics

1) In the above tutorial, the muscles are SafeDownCast - is that necessary in this situation? I can see the metabolics probe added to the componentset if I print the model without downcasting the muscles.

2) Does the smoothed version of this probe not count as the class Probe, only as a Component? If I call addProbe() instead of addComponent(), I get an error.

3) The smoothed version, is recommended if you want to use a metabolics as a cost in your optimization. Will the non-smoothed version (also the Umberger2010 model) give any problems if you add it to the model, but don't try to minimize the output (i.e. just using it to check MCOT post-hoc)?

4) I've been trying to compute the metabolic cost using the analyze tool, I was able to do so using:

Code: Select all

outputPaths_metabolics = StdVectorString();
outputPaths_metabolics.add('.*total_metabolic_rate');
metabolics_table=study.analyze(trajectory, outputPaths_metabolics);

This works, but I was initially having some trouble assigning the output Path more directly. I tried (outputPaths_metabolics.add('/metabolic_cost|total_metabolic_rate')) and also ('/components/metabolic_cost|total_metabolic_rate'). Here, 'metabolic_cost' is the probe as named in the tutorial. The analyze function didn't recognize the outputs when assigning them in this way, am I doing something wrong here?

Cheers,
Pasha

[Nicholas Bianco]

Hi Pasha,

1) In the above tutorial, the muscles are SafeDownCast - is that necessary in this situation? I can see the metabolics probe added to the componentset if I print the model without downcasting the muscles.

Do you get the same result after running the optimization? It's possible to add a metabolics component that has no muscles associated with it. The argument type in addMuscle() is the base class "Muscle", and in Matlab, funny stuff might happen if you don't use safeDownCast().

2) Does the smoothed version of this probe not count as the class Probe, only as a Component? If I call addProbe() instead of addComponent(), I get an error.

Bhargava2004SmoothedMuscleMetabolics is not a Probe, it's just a regular Component.

3) The smoothed version, is recommended if you want to use a metabolics as a cost in your optimization. Will the non-smoothed version (also the Umberger2010 model) give any problems if you add it to the model, but don't try to minimize the output (i.e. just using it to check MCOT post-hoc)?

No, it will not give you any problems.

4) I've been trying to compute the metabolic cost using the analyze tool

We use regular expressions (regex) by default to parse the output paths you provide, and the vertical bar character, "|", has a specific function in regex. To avoid the regex parser from interpreting the bar character, you can "escape" it by putting a forward slash before it:

Code: Select all

outputPaths_metabolics.add('/metabolic_cost/|total_metabolic_rate');

Sometimes you might need two forward slashes (I think only in Python on Windows):

Code: Select all

outputPaths_metabolics.add('/metabolic_cost//|total_metabolic_rate');

Best,
Nick