Problem of activation estimation for Thelen and Millard muscles

[Bach Quoc Hoa]

Hello,

I'm working on the simple model shown in the attached photo. The foot is welded to the ground and connected to the tibia by a pin joint at the ankle. The only two actuators are those representing soleus and the tibialis anterior muscles.

I already have 3 seconds of state and torque profiles with a step time of 1 ms (so each profile consists of 3000 points of reference). With these inputs, I want to estimate the activities of these muscles reproducing the given torque trajectory with the following algorithm

for i = 1 to number of torque trajectory points
_____update state from file (changing state without integration)
_____equilibrate muscle (to calculate muscle fiber length)
_____optimization until cost function value < 1e-6
__________input: torque reference value
__________Output: muscle activations
__________Cost function: minimizing the error between reference and computed torques and the sum of squared muscle forces
end

The above algorithm was tested with three types of actuators: Path Actuator, Millard2012EquilibriumMuscle and Thelen2003Muscle. It worked perfectly fine with Path Actuator. For the two other types of muscle, the optimization function kept trying the same values of control inputs and therefore was stuck right at the first iteration (below are the values printed directly from inside the objectiveFunc function).

newControls = ~[2.71881e-06 0]
f = 18.8932
newControls = ~[-2.71881e-06 0]
f = 18.8932
newControls = ~[0 2.71881e-06]
f = 18.8932
newControls = ~[0 -2.71881e-06]
f = 18.8932
newControls = ~[0 0]
f = 18.8932
newControls = ~[0 0]
f = 18.8932
newControls = ~[0 0]
f = 18.8932
newControls = ~[2.71881e-06 0]
f = 18.8932
newControls = ~[-2.71881e-06 0]
f = 18.8932

Also for Thelen2003Muscle, there's this error

Thelen2003Muscle initialization: tibant_r is at its minimum fiber length of 0.000000

I tried to set the fiber length after initializing the state and before running the algo above, but nothing changed. The default fiber length are 0.1.

From these observations, the algorithm seems not to work with the muscles having a dynamic model. Could you help me to find out what I missed here?

Below are the code snippets

For the algorithm

Code: Select all

for (int nbincrements = 0; nbincrements < nbStateIterations; ++nbincrements) {
        
        // Update state
        StateVector* state = stateRef.getStateVector(nbincrements);
        si.setTime(state->getTime());

        for (int i = 0, j = 0; i < state->getSize(); i +=2, ++j) {
            osimModel.updCoordinateSet().get(j).setValue
                (si, state->getData()[i]);
            osimModel.updCoordinateSet().get(j).
                setSpeedValue(si, state->getData()[i+1]);
        }

        osimModel.equilibrateMuscles(si);
	osimModel.realizeAcceleration(si);
        
        // Torque to be tracked
        int nbcoords = osimModel.getNumCoordinates();
        Vector torqueToTrack(nbcoords, 0.0);
        torqueRef.getDataAtTime(state->getTime(), nbcoords, torqueToTrack);
        
     
        //__________________________________
        // Moment arm matrix
        auto& coordinateSet = osimModel.getCoordinateSet();
        vector<vector<double> > momentArm;

        for (int i = 0; i < muscleSet.getSize(); i++) {
            auto muscle = dynamic_cast<const Millard2012EquilibriumMuscle*>(&muscleSet[i]);
            if (muscle) {
                vector<double> column;
                for (int j = 0; j < coordinateSet.getSize(); j++) {
                    column.push_back(muscle->computeMomentArm(si, coordinateSet[j]));
                }
                momentArm.push_back(column);
            }
        }
        // convert to Matrix
        Matrix R(momentArm[0].size(), momentArm.size());
        for (unsigned int m = 0; m < momentArm.size(); m++) {
            for (unsigned int n = 0; n < momentArm[0].size(); n++) {
                R[n][m] = momentArm[m][n];
            }
        }

	Vector controls = osimModel.updDefaultControls();
     	Real f = 1e3;
        while (f > 1e-6)) {
                TrackingTorqueSystem sys(numControls, si, osimModel, R, torqueToTrack, regweight);

                /* Define initial values and bounds for the controls to optimize */
                Vector lower_bounds(numControls, 0.01);
                Vector upper_bounds(numControls, 0.99);

                sys.setParameterLimits( lower_bounds, upper_bounds );

                // Create an optimizer
               Optimizer opt(sys, BestAvailable);

                // Specify settings for the optimizer
                opt.setConvergenceTolerance(1e-6);
                opt.useNumericalGradient(true);
                opt.setMaxIterations(2000);
                opt.setLimitedMemoryHistory(500);

                // Optimize it!
                f = opt.optimize(controls);
	}
	osimModel.updDefaultControls() = controls;
}

For the Optimizer System

Code: Select all

class TrackingTorqueSystem : public OptimizerSystem {

    public:

        TrackingTorqueSystem(int numParameters, State& s, Model& aModel, Matrix& momentarms, Vector& torqueRef, double weight):
            OptimizerSystem(numParameters), 
            _state(s), 
            _osimModel(aModel),
            _momentarms(momentarms),
            _regWeight(weight),
            _torqueRef(torqueRef)
{}

        int objectiveFunc(const Vector &newControls, bool new_coefficients, Real& f) const override {

            State s = _state;

	    _osimModel.updDefaultControls() = newControls;	//(edited)

             _osimModel.equilibrateMuscles(s);
            _osimModel.realizeAcceleration(s);

            Vector muscleForces(muscleSet.getSize(), 0.0);

            for (int i = 0; i < muscleSet.getSize(); i++) {
                auto muscle = dynamic_cast<const Millard2012EquilibriumMuscle*>(&muscleSet[i]);
                muscleForces(i) = muscle->getTendonForce(s);
            }


            Vector diff = _momentarms * muscleForces - _torqueRef;
            f =  ~diff * diff  + ~muscleForces * muscleForces * _regWeight;

            //cout << "newControls = " << newControls << endl;
            //cout << "f = " << f << endl;

            return (0);
        }


    private:

        State& _state;
        Model& _osimModel;
        Vector& _torqueRef;
        Matrix& _momentarms;
        double _regWeight;

}; 

[Bach Quoc Hoa]

Is there anyone who can help me?

[Thomas Uchida]

Please confirm that the problem you are asking the optimizer to solve actually has a solution. It looks like you are specifying the kinematics and, therefore, the length of the muscle--tendon actuators. For each muscle, does a muscle fiber length exist where the muscle is generating the necessary force and the fiber is in equilibrium with the tendon? You can try solving this by hand for the first time point that the optimizer cannot solve. You may need to adjust the model parameters (optimal fiber lengths, tendon slack lengths, maximum isometric forces, etc.) to make the problem solvable.

[Bach Quoc Hoa]

Thank you for getting to me!

My solution for the previous mentioned problem is to replace updDefaultControls() in my personalized class OptimizerSystem with setActivation. However, I still can use updDefaultControls outside this class (in the main function) to update the muscle controls.

So it changes from

Code: Select all

        int objectiveFunc(const Vector &newControls, bool new_coefficients, Real& f) const override {

            State s = _state;

	    _osimModel.updDefaultControls() = newControls;
	    
            _osimModel.equilibrateMuscles(s);
            _osimModel.realizeAcceleration(s);

to

Code: Select all

        int objectiveFunc(const Vector &newControls, bool new_coefficients, Real& f) const override {

            State s = _state;

	    auto& muscleSet = _osimModel.getMuscles();
	    for(int i=0; i< muscleSet.getSize(); i++ )
	    	muscleSet[i].setActivation(s, newControls(i));
	    
            _osimModel.equilibrateMuscles(s);
            _osimModel.realizeAcceleration(s);

I have some more questions though.

1. osimModel.updDefaultControls() and muscle.setActivation (state, controls): what are the differences between these two functions?
I just luckily found the solution but don't understand at all why it worked. Plus, the old code had worked for PathActuator.

2. osimModel.equilibrateMuscles() : what does this function do exactly?
Specifically, does this function compute the muscle forces by taking into account a static pose (e., zero muscle-tendon velocity) ?
Should this function be called every time when at least one muscle has its activation updated?
If realizeAcceleration is called, is it necessary to call equilibrateMuscles?

Thank you in advance!