MTP Optimizer Crashing

[SpencerTWilliams]

I think the issues with MTP crashing were caused by two bugs. First, there was an issue with calculating model moments producing infinite values. The problem was in passiveForceLengthCurve.m, where an intermediate step in an equation used big exponentials that sometimes went over the limit of what a double could store. I rearranged the log equation so the solution is the same without the intermediate numbers getting too big:

function normalizedPassiveForce = passiveForceLengthCurve(...
    normalizedMuscleFiberLength)

e1 = 0.232000797810576;
e2 = 12.438535493526128;
e3 = 1.329470475731338;
% Result of exp can be greater than a double can store
% normalizedPassiveForce = e1 * log(exp(e2 * (... 
%     normalizedMuscleFiberLength - e3)) + 1);

fiberLengthPower = e2 * (normalizedMuscleFiberLength - e3);
normalizedPassiveForce = e1 * (log(exp(0.5 * fiberLengthPower) + ...
    exp(-0.5 * fiberLengthPower)) + log(exp(0.5 * fiberLengthPower)));
end

The other issue was related to Inf, -Inf, and NaN values for the normalized fiber length penalty. I think the problem was in findCorrectMtpValues.m, in this block:

if (valuesStruct.isIncluded(index))
    [startIndex, endIndex] = findIsIncludedStartAndEndIndex( ...
        valuesStruct.primaryValues, valuesStruct.isIncluded, index);
    % Secondary values are all zero
    output = valuesStruct.secondaryValues(startIndex:endIndex);
else
    output = valuesStruct.primaryValues(index, :);
end

valuesStruct.secondaryValues was defined to contain all zeroes, while primary values are real initial values. I changed the first output line to this:

output = valuesStruct.primaryValues(startIndex:endIndex);

When it was pulling zeroes from the secondary values array, it would cause a later step to divide by zero, so I think my change is correct.

The optimizer now runs, but this is the output, so I'm not sure if everything is right:

Initial point is a local minimum that satisfies the constraints.

Optimization completed because at the initial point, the objective function is non-decreasing 
in feasible directions to within the value of the optimality tolerance, and 
constraints are satisfied to within the value of the constraint tolerance.

The program then crashes at the unfinished results reporting in MusclePersonalizationTool.m.

[SpencerTWilliams]

Committed b08fadcf3583033f782409a785bfee7058d748ba like we talked about in the meeting today to revert findCorrectMtpValues and change the isIncluded vector to zeros instead of ones. The behavior in testing is the same: optimizer does not crash but is completed at the initial point, and the program crashes at results reporting.

[SpencerTWilliams]

I found that calcNormalizedMusceFiberLengthsAndVelocities() sometimes comes up with negative values for fiber length, and it's caused by issues with units. Muscle-tendon length is given in meters, but optimal fiber length and tendon slack length are converted from meters (their units in the .osim model) to cm when imported in getModelInputs() called by parseMuscleTendonPersonalizationSettingsTree(). getModelInputs() also converts the pennation angle from radians to degrees, but calcNormalizedMusceFiberLengthsAndVelocities() uses cos(), which expects radians. The angular units and subtracting a value in cm from a value in meters can both generate negative lengths.

Is there a reason these units need to be converted during the import step, and would it be better to convert them back to their original units while used in calcNormalizedMusceFiberLengthsAndVelocities()?

When I correct the units in calcNormalizedMusceFiberLengthsAndVelocities(), this is the optimizer result:

Local minimum possible. Constraints satisfied.

fmincon stopped because the size of the current step is less than
the value of the step size tolerance and constraints are 
satisfied to within the value of the constraint tolerance.

Optimization stopped because the relative changes in all elements of x are
less than options.StepTolerance = 1.000000e-06, and the relative maximum constraint
violation, 0.000000e+00, is less than options.ConstraintTolerance = 1.000000e-06.

Is there a way I can check how many iterations fmincon took to get this? It doesn't say it converged on the initial point like last time, but it seems weird that the constraint violation would be perfect.

[cvhammond]

The units for the MTP cost function should be (meters, radians) so any issues with getModelInputs() should be rectified to ensure data is loaded into those units appropriately