Problem of activation estimation for Thelen and Millard muscles
Posted: Wed Nov 11, 2020 1:46 am
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).
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
For the Optimizer System
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).
Also for Thelen2003Muscle, there's this errornewControls = ~[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
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.Thelen2003Muscle initialization: tibant_r is at its minimum fiber length of 0.000000
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;
};