Moco goal for contact parameter optimisation using weights

[Oliver Demuth]

Hi everyone,

I am trying to optimise the contact model, and especially the contact radii of my contact spheres and the stiffness of the SmoothSphereHalfSpaceForce, based on GRF and kinematic data. I use MocoInverse to prescribe my kinematics and a MocoContactTrackingGoal to track my GRF data. I currently use MocoParameters to optimise both the stiffness and the contact radii but there appears to be some local minimum when the contact radii reach the upper bounds. How can I add weights to the parameter optimisation, so that they ideally diverge as little as possible from the initial value?

Here is how I set up the MocoParameters:

Code: Select all

% ==== Add parameters to optimise to the problem ==== %
        
% optimise stiffness across all contact forces
stiffnessParam = MocoParameter('contactStiffness',StdVectorString(HalfSpaceForces),'stiffness',MocoBounds(1e-5,1e6));
problem.addParameter(stiffnessParam); % append contact stiffness to problem
        
% optimise contact geometry radii
contactRadii = MocoParameter('contactRadii',StdVectorString(contactGeometryPathNames),'radius',MocoBounds(0.0025,0.05));
problem.addParameter(contactRadii); % append contact radii to problem

I tried MocoOutputGoal (or MocoInitialOutputGoal or MocoFinalOutputGoal) but I am not sure how to 'connect' setOutputPath of MocoOutputBase with either the MocoParameter or the contact geometry directly.

When I try to connect it directly with the contact geometry like this:

Code: Select all

% ==== Add a MocoOutputGoal to the problem ==== %

radiiGoal = MocoOutputGoal('radii');
problem.addGoal(radiiGoal);
radiiGoal.setOutputPath(append(contactGeometryPathNames(1),'|radius'));
radiiGoal.setWeight(radiiWeight);

I get an error saying that no output called radius exists at that component, i,e.:

java.lang.RuntimeException: no Output 'radius' found for this Component.
In Object 'Digit_I_R_CS' of type ContactSphere.
Thrown at Component.h:1163 in getOutput().

Which makes sense as it is an attribute/property and not an output, so directly linking it to the contact geometry might not be the way forward.
However, I can't seem to define a path to (or find a way to conect) the MocoParameter's output. I assume the optimised value from the MocoParameter is found at contactRadii.getPropertyByIndex(3) (i.e., the porperty_element) but how do I set it up correctly with the MocoOutputGoal?

Many thanks for your help!

Cheers,
Oliver

[Nicholas Bianco]

Hi Oliver,

It is currently not possible to add cost term to minimize deviations from the initial parameter value. The SmoothSphereHalfSpaceForce properties are not available as Outputs, so they cannot be used with MocoOutputGoal.

I've considered adding support for "parameter regularization" as you describe, where you minimize the deviation of the parameter values from either 1) the initial guess or 2) the bounds midpoint value. But unfortunately this do not exist yet.

Good to know that others are interested! Will keep this in mind when prioritizing features in the next release.

Best,
Nick

[Oliver Demuth]

Hi Nick,

Thank you for your clarification!

Good to know that muscle parameter tuning (i.e., optimal fibre length or tendon slack length) in this way is also not yet possible in Moco, which would have been my next step. I assume this means that it is currently only possible in the the Muscle Redundancy Solver of De Groote and colleagues as implemented here by Bishop and colleagues.

Alternatively, is there a way to create a custom MocoGoal for this purpose? How would the output of the optimised MocoParameter best be fed into such a custom goal? How can I access the MocoParameters' output after each itteration?

Cheers,
Oliver

[Nicholas Bianco]

Hi Oliver,

If you need more customization in your parameter optimization workflow, then the Muscle Redundancy Solver is a great alternative (I used it for my research before we created Moco).

You could implement a custom MocoGoal, since parameters are applied to the model and will be available when evaluating the cost. However, a better, more general approach would be to make parameter regularization a solver option. I'll file an issue to track this feature request.

Best,
Nick

[Oliver Demuth]

Hi Nick,

Thank you for your reply. I tried to formulate such a custom MocoGoal to regularise MocoParameters. Generally the problem set up is as follows:

Code: Select all

% Construct the MocoInverse
inverse = MocoInverse();
modelProcessor = ModelProcessor(model); 
inverse.setModel(modelProcessor);

% initialise Moco Study
study = inverse.initialize();
problem = study.updProblem();

% optimise contact geometry radii
contactRadii = MocoParameter('contactRadii',StdVectorString(contactGeometryPathNames),'radius',MocoBounds(0.0025,0.05));
problem.addParameter(contactRadii); % append contact radii to problem

% configure the solver and solve the problem
solver = MocoCasADiSolver.safeDownCast(study.updSolver()); % get solver
solver.resetProblem(problem); % reset problem

% get MocoTrajectory
trajectory = solver.createGuess();

The cost function would be relatively straight forward:

Code: Select all

% Cost function for contact radii parameter regularisation
function value = calcMyCustomGoalValue(problem, trajectory)
	Ropt = trajectory.getParameter('contactRadii'); % get optimised parameter
	model = problem.getPhase(0).getModelProcessor().process(); % get model
	Rinit = org.opensim.modeling.ContactSphere.safeDownCast(model.getContactGeometrySet.get(0)).getRadius; % get initial value
	value = (Ropt/Rinit - Rinit/Ropt)^2; % calculate cost value
end

or alternatively, the default/target value could be set manually, which would make it more universally applicable:

Code: Select all

% Cost function for general parameter regularisation
function value = calcMyCustomGoalValue(trajectory, target, targetName)
	opt = trajectory.getParameter(targetName); % get optimised parameter
	value = (opt/target - target/opt)^2; % calculate cost value
end

where target is the reference value (e.g., 0.5), and targetName represents the MocoParameter, in this case 'contactRadii'.

I have tried to translate this into C++, however, I have to say that I am not experienced with C++ at all, and it potentially has some (many?) issues. Here is what I have come up so far. For the MocoParameterRegularisationGoal.h:

Code: Select all

#ifndef OPENSIM_MOCOPARAMETERREGULARISATIONGOAL_H
#define OPENSIM_MOCOPARAMETERREGULARISATIONGOAL_H

#include "MocoGoal.h"

namespace OpenSim {

/** The squared difference between a reference value and an optimised 
MocoParameter.
@ingroup mocogoal */

class OSIMMOCO_API MocoParameterRegularisationGoal : public MocoGoal {
    OpenSim_DECLARE_CONCRETE_OBJECT(MocoParameterRegularisationGoal, MocoGoal);
public:
    MocoParameterRegularisationGoal() { constructProperties(); }
    MocoParameterRegularisationGoal(std::string name) : MocoGoal(std::move(name)) { 
        constructProperties(); 
    }
    MocoParameterRegularisationGoal(std::string name, double weight)
            : MocoGoal(std::move(name), weight) { 
        constructProperties(); 
    }
    
    /// The name of the MocoParameter in the MocoTrajectory which should be 
    /// regularised
    void setParameterName(std::string parameterName) {
        set_parameter_name(std::move(parameterName))
    }

    /// Set desired value of the parameter
    void setReferenceValue(double referenceValue) {
        set_reference_value(std::move(referenceValue))
    }

protected:
    void initializeTrajectoryImpl(const MocoTrajectory&) const override;
    void calcGoalImpl(
            const GoalInput& input, SimTK::Vector& cost) const override; 

private:

    OpenSim_DECLARE_PROPERTY(parameter_name, std::string,
        "The name of the MocoParameter in the MocoTrajectory which should be "
        "regularised.");
    OpenSim_DECLARE_PROPERTY(reference_value, SimTK::Vec3,
        "The desired reference value of the MocoParameter (i.e., initial or "
        "optimal value)");

    void constructProperties();

    mutable SimTK::ReferencePtr<const MocoParameter> m_parameter;
};

} // namespace OpenSim

#endif // OPENSIM_MOCOPARAMETERREGULARISATIONGOAL_H

And the MocoParameterRegularisationGoal.cpp:

Code: Select all

#include "MocoParameterRegularisationGoal.h"

#include <OpenSim/Moco/MocoTrajectory.h>

using namespace OpenSim;

void MocoParameterRegularisationGoal::initializeTrajectoryImpl(const MocoTrajectory& trajectory) const {
    m_parameter.reset(&trajectory.getParameter<MocoParameter>(get_parameter_name()));
}

void MocoParameterRegularisationGoal::calcGoalImpl(
        const GoalInput& input, SimTK::Vector& cost) const override {
    cost[0] = ((m_parameter/get_reference_value())-(get_reference_value()/m_parameter))^2.;
}

void MocoParameterRegularisationGoal::constructProperties() {
    constructProperty_parameter_name("");
    constructProperty_reference_value(double(0));
}

I am especially unsure if I successfully connected the parameter from the MocoTrajectory to the cost function of the goal. I would be extremely appreciative for any comments regarding how to improve it and/or fix it (I am just generally assuming it doesn't work yet).

Many thanks for your help!

Cheers,
Oliver