Contact force oscillations/parameter tuning

[Pasha van Bijlert]

Hi all,

I'm doing forward simulations, and I'm trying to tune my contact force parameters so that they "look right". I am simulating a large dinosaur (8000kg), so finding data to compare it is a bit of a challenge, but I've found the following measured GRF data of a running ostrich (Rankin et al 2016):

rankin_ostrich_run.PNG (85.24 KiB) Viewed 2268 times

Rankin's measured data was low pass filtered at 10Hz, so it's bound to be smoother than the "real" value. However, I'm not sure low-pass filtering my simulation data is the way to go (because I want to use mid-shaft joint reaction moments to calculate hindlimb stress). In my own simulations, I get a severe impact peak, and strong oscillations in the x component:

low_stiff.PNG (32.01 KiB) Viewed 2268 times

After lowering the dissipation to 0.2, the impact peak is reduced (but not eliminated):

low_stiffdisp.PNG (34.74 KiB) Viewed 2268 times

I was wondering if anybody has any advice on how to 1) eliminate the impact peak and 2) reduce or eliminate the horizontal oscillations.

I believe setting dissipation to zero will reduce the impact peak further (but not fully eliminate it, I'm currently running an optimization to try this out). I'm not sure how I can tackle the horizontal oscillations. They're implying a skidding motion happening upon contact, but I'm not seeing that in the visualizer. I think changing the friction coefficients is the way to go, but there are four parameters and I'm not 100% sure which one might be causing the oscillations. My friction coefficients are: 0.8 static friction, 0.4 dynamic and viscous friction, and a transition velocity of 0.2 m/s. Has anybody dealt with a similar issue?

Cheers,
Pasha

[Ross Miller]

Hi Pasha,

In predictive simulations, you may have a hard time getting rid of impact peaks. They kind of come with the territory, in my experience.

Ackermann & van den Bogert (2010) reported impact peaks in predictive simulations of human walking were smaller when using "minimum fatigue"-type cost functions (sum of muscle excitations raised to a power) than when using "minimum effort"-type cost functions (sum of muscle excitations, weighted by muscle mass and raised to a power). So you might try different cost functions.

Oscillations in the GRF I think can be artifacts of the transcription scheme. Ton and Anne Koelewijn informed me of this some years ago, maybe one of them will see this and comment. Anne had a solution to it that involved averaging but I forget the details. Not sure if this is an issue with Moco's schemes but you could try other schemes. I think Moco has at least two. I used to see this in simulations of sprinting when I used a Midpoint Euler scheme but not when I used a Backward Euler scheme.

Ross

[Pasha van Bijlert]

Hi Ross,

Thanks for your advice. Regarding fatigue vs. effort cost functions, I'm using an unweighted sum of controls cubed, which I think would be classified as "minimum fatigue", although I remember that paper trying out higher exponents as well. I'll try out different exponents, thanks for the tip! Coincidentally, Anne Koelewijn recently mentioned to me that too high exponents in your cost function can cause convergence issues when using IPOPT. My cost function also contains joint reaction forces and moments from 6 (weld)joints in the model, which I'd hoped might prevent too much "jerkiness" in the gait. I'll see what exponents I can get away with.

About the transcription scheme: I'm using the default (hermite simpson). Do you think it's worth trying out a trapezoidal scheme? Or am I conflating solver settings now?

Cheers,
Pasha

[Ross Miller]

Yes, by "transcription" I'm referring to that temporal collocation grid discretization scheme. I don't know if Hermite-Simpson would necessarily have this issue but it's worth trying out other schemes available in Moco, it would at least help identify the problem or rule out causes.

Ross

[Ross Miller]

Also, if I remember correctly, Anne told me that this artifact is removed by doing some sort of averaging of the GRF in nearby nodes (which is okay since it's an "artifact", not the "real" GRF that's moving the model). So if you identify that this is being caused by the collocation scheme, you may not actually have a problem here.

I don't remember the details of the averaging procedure. I sent Anne this thread to see if she remembers them.

Ross

[Ton van den Bogert]

The impact can certainly be reduced if there is something in the cost function to discourage it. It sounds like you already have that. Maybe it's not weighted high enough.

Ross is right to remind you that higher exponents may cause issues with IPOPT. We went up to 10 (Ackermann 2010) but that was with the SNOPT solver. Also, IPOPT is somewhat sensitive to scaling, so try to scale your objective function (the very first column in the IPOPT output) to be between 0.1 and 10.

The oscillations in the horizontal GRF are annoying but likely an artifact of the discretization scheme and/or the friction model, and not affecting your other results much. It would not hurt to try the Trapezoidal discretization, although theoretically, that should be more sensitive to oscillations than the Hermite-Simpson which has some numerical damping.

It may be better to have the same friction coefficients for static and dynamic friction. The MOCO models also specify a "tanh smoothing" parameter and I am not sure how its value should be chosen. More smoothing may be needed, so you can no longer get large variations in friction force from tiny variations in slip velocity. It is possible that contact models that worked in human models need some tweaking for your dinosaur scale.

And some less elegant solutions:

• discourage the oscillations by penalizing the optimization objective for non-smoothness. The mathematical definition of smoothness is the squared second derivative of a state trajectory. You will find that this penalty can easily be weighted so it prevents the oscillations without affecting any predictions.

• (as Anne suggested via Ross) average two successive force values to get a smooth curve, for presentation. The original and averaged force trajectories are not really mechanically different. The human (and especially dinosaur) body act as low-pass filters and don't respond to high frequencies.

Ton

[Pasha van Bijlert]

Hi Ross, hi Ton,

Thank you for all of the suggestions! You've given me a lot of ideas to think about. I'd also e-mailed Anne, who mentioned that her oscillations were present in all the states, not just the ground contacts. This leads me to believe the issue is the friction parameters.

You are, of course, correct that the animal's inertia acts as an LPF of the contact forces. I'm determining hindlimb stresses by adding in an extra weldjoint at the midshaft location of each segment, and using the joint-reaction moments and forces combined with Euler-Bernoulli beam theory. Dinosaurs effectively have a three-segment leg: femur, shank, and metatarsals. The shank & femurs are sufficiently downstream to see this inertial filtering/dampening effect, but the metatarsals are not, while stress on the metatarsals might be limiting the maximal speed of the animal (as per Sellers 2017). Hence my extra care in trying to minimize the contact oscillations and impacts.

I've weighted the joint reaction moments in the foot more heavily, but I will increase the weight further (my objective is scaled to be below 10).

Regarding the scaling of the contact parameters: I scaled stiffness according to the mass of the animal. I'm using a rigid foot+toes segment, and to compensate for that I'm allowing some deformation in the contact model. This represents the fat pad below the foot, and of course the deformation that would occur from simply being able to bend the toes to absorb the impact. This, combined with setting dissipation to 0, has smoothed out the severe impact peak I plotted earlier (it's still present, but in the order of 2BW during a moderate run).

I'd forgotten about the smoothing parameters of the contact force, I will give that a shot, and also see what effect keeping a constant friction coefficient for static and dynamic situations might have.

Regarding the smoothness/regularization of the accelerations: I believe MOCO only permits this if you use implicit multibody dynamics (and you cannot square the accelerations, only add them in linearly). I will give this another shot, and also see if the different transcription method affects the friction at all. The states in my simulation don't seem to have oscillations, though, so I'm not sure if this will have an effect on the GRF.

Thanks again!
Pasha

[Pasha van Bijlert]

Hi all,

I thought I'd just report back on some of the things I've tried.

Setting all the friction values to the same magnitude reduced or eliminated the horizontal GRF oscillations. This did not eliminate the vertical peak GRF. I could reduce the impact peak to about 2xBW by setting contact dissipation to 0, but it is debatable whether this is still a good representation reality (because your model is now walking on lossless springs). Even if I set it to 0, I didn't get a GRF shape that looks like the perfect bell curve from the ostrich forceplate data, but that is of course because it is low-pass filtered. We did some unfiltered (human) running trials, and while there is an impact peak, the overall vertical GRF shape is dissimilar to the simulations (the simulations preferentially have a higher magnitude at the start than at the end).

I found that minimizing joint accelerations was a more effective way to reduce the contact oscillations (even without setting all the friction parameters to the same magnitude). You can do this in moco using "solver.set_minimize_implicit_multibody_accelerations(true)",which defaults to an exponent of 2 (so you're minimizing squared accelerations). This obviously also requires using implicit multibody dynamics, which takes more time per iteration, but I've found overall optimization times to be roughly the same or slightly lower.

I've not yet tried different collocation schemes. Thanks again for all the suggestions & discussions in this thread!

Cheers,
Pasha

Edit: I incorrectly assumed the exponent of the accelerations was 1 when you're adding them in as a cost, but it's 2

[Nicholas Bianco]

As I mentioned on the other thread, "solver.set_minimize_implicit_multibody_accelerations(true)" does introduce a squared cost term for the accelerations.

-Nick

[Pasha van Bijlert]

HI Nick,

Thanks, I'll edit my post to reflect this.

Cheers,
Pasha