Minimizing duty factor in a running simulation using MocoOutputGoal?

[Pasha van Bijlert]

Hi all,

I've been trying to think of ways to increase or decrease duty factors in a running simulation. The release notes for OpenSim 4.5 say

Foot-ground contact forces can now be optimized (e.g., maximize contact impulse) via MocoOutputGoal using the new Outputs 'sphere_force' and 'half_space_force' added to SmoothSphereHalfSpaceForce.

, prompting this thread.

I think this goal would allow you to add a cost that minimizes duty factors by adding the vertical component of sphere_force to the goal, setting the exponent higher than 1 so the goal is focused on the peak values, and then setting the weight negative.

1) I've read in earlier forum posts that setting a negative weight results in the maximization of a goal, instead of minimization. I don't fully understand how it works. I thought the optimizer reduces costs to zero, but then the sign of the weight shouldn't matter. Does the optimizer attempt to reduce costs to -inf?

2) If my model has 10 contact-forces per foot, would I need to add a MocoOutputGoal for each force element?

I'm also open to trying other approaches - I was considering is adding differently shaped GRF curves as a tracking cost, but maybe there's a more direct way to constrain duty factors that I'm not thinking of.

Cheers,
Pasha

[Nicholas Bianco]

Hi Pasha,

I thought the optimizer reduces costs to zero, but then the sign of the weight shouldn't matter. Does the optimizer attempt to reduce costs to -inf?

Yes, the optimizer can minimize the objective function below zero. We typically have non-negative cost terms so this is less common, but it is totally valid.

2) If my model has 10 contact-forces per foot, would I need to add a MocoOutputGoal for each force element?

Yes, you will need one for each foot.

Best,
Nick

[Pasha van Bijlert]

Hi Nick,

Adding the following goal throws no errors.

Code: Select all

%% Add peakGRFy maximization goal
forceset =  model.getForceSet()

contactforcenames = {}; #to be reused in a different section of the script
weight = -0.01;

for f = 1:(forceset.getSize)
    force = forceset.get(f-1);
    CCN = force.getConcreteClassName().char;
    if contains(CCN,'SmoothSphereHalfSpaceForce')


        force = SmoothSphereHalfSpaceForce.safeDownCast(force);
        forcename = force.getName.char;
        contactforcenames{end+1} = forcename;

        GRFGoal = MocoOutputGoal([forcename '_peakGRFy'],weight);
        GRFGoal.setOutputPath(['/forceset/'  forcename  '|sphere_force'])
        GRFGoal.setOutputIndex(1) %is this GRF_y? Or should the index be 4?
        GRFGoal.setExponent(2);
        problem.addGoal(GRFGoal);
        
    end
end

I'm unsure whether the OutputIndex should be 1 or 4, if I want the y Component of the force minimized. The documentation says the first three components of a SpatialVec are the rotational components, but when using force.getRecordValues you get the force components first, then the torque components.

Also, the above code adds 20 extra goals to my problem, because I have 10 contact forces per foot, which substantially slows things down.

Yes, you will need one for each foot.

Is there some way to group all the force outputs per foot, and then minimize/maximize the grouped/totaled output?

Thanks!
Pasha

[Nicholas Bianco]

Hi Pasha,

The documentation says the first three components of a SpatialVec are the rotational components, but when using force.getRecordValues you get the force components first, then the torque components.

MocoOutputGoal uses the SpatialVec convention, so you'll want to use "GRFGoal.setOutputIndex(4)" to pick out the GRFy value.

Is there some way to group all the force outputs per foot, and then minimize/maximize the grouped/totaled output?

Unfortunately no, at the moment. It would be a good idea to find a way to combine multiple Outputs into a single goal to avoid the overhead. One possible workaround is to create an additional component that sums all the Output values from other components and provides it as another Output value, but this still requires C++ coding.

Best,
Nick