Exponential expressions and forward dynamic simulation

[Peter Bishop]

Hello,
I've been looking at the incorporation of passive joint moments into my model. Presently, I have been using the CoordinateLimitForce class to represent these passive moments, but as I understand this only involves a single DOF, whereas ideally I would like to involve all DOFs at a joint (up to three for e.g., the hip) simultaneously. So I have been looking at using the ExpressionBasedBushingForce class, which does allow for multiple DOFs to be handled at once (i.e., the passive moment about one coordinate/axis can be a function of the values of multiple coordinates). Moreover, as passive joint moments have typically been represented with double exponential functions in the literature, I am seeking to implement this kind of function in my ExpressionBasedBushingForce.

The form of the function I have been testing is: Moment = a*exp(b*theta1-c)-d*exp(e*theta1-f), for a single angle. If Moment is a function of multiple angles, then I have additional terms for theta2, theta3, etc.

I've tested how this implementation performs, and it behaves exactly as expected in some basic inverse dynamics test. However, when I try to perform a forward dynamic simulation, I run into problems. Even in a basic 'drop test' (and I've got appropriate passive forces on the pelvis segment which, are a function of vertical displacement and so prevent the pelvis falling through the floor), the simulation will progress so far and then it 'fails' - the model disappears from view and all coordinate values/speeds just become question marks. In the output states file, at these time instances all values are reported as "-1.#IND0000", something I've never come across before. I get the same problem whether I have ground contact forces enabled or not.

If I force the integrator to use smaller timesteps (i.e., reduce the maximum step size value to something like 0.001), this ameliorates the situation somewhat: more of the simulation will now run/solve properly, but sooner or later (if I have set a long enough time range) the problem happens.

If I use a polynomial-based function in my ExpressionBasedBushingForce, I don't run into these problems. I am therefore wondering if the use of the exponential function here is causing the problems, even though it's a 'legal' function, and works fine in inverse dynamics. I am hoping to use exponential functions in some way, because like I said that's how empirical passive joint moment data has been typically modelled previously.

If anybody has any suggestions for what I can do to fix this, I would be most appreciative!

Many thanks,
Peter Bishop.

[Ton van den Bogert]

Peter,

It is well known that exponential force fields make simulation difficult, and you have already found that the time step is the issue.

Your ODE solver will use variable time steps. It will first try one particular step size, evaluate the differential equations, and then reduce the step size if needed. If the angular velocity is large enough, the change in angle could be large enough that exp(b*theta-c) becomes infinity during that initial evaluation. exp(800) is already infinity in a standard double precision calculation. If the ODE solver does not detect this problem, all subsequent calculations based on this initial evaluation will be "not a number" and that is what you are seeing (https://en.wikipedia.org/wiki/NaN#Display). Because of that initial failure, the ODE solver is never able to reduce the step size to avoid infinities and solve the problem.

I suspect that the Matlab ODE solvers will be better at detecting infinities and responding appropriately, but I have not tried.

You could follow up with the authors of the ODE solver to make it work better. However, you will still end up with extremely small time steps. This is slow, and probably not even needed.

You are right that exponential functions fit the passive joint moment data well, but it makes simulation difficult. Hatze (The complete optimization of a human motion, Math Biosci, 1976) used exponential functions for all nonlinear elastic elements. He used a fixed step Runge-Kutta solver and mentions "It was, however, found that the integration interval had to be reduced to 0.0005 sec because of certain stability problems". He did not directly blame the exponentials, and I do not remember how I came to that conclusion, but I have always avoided exponential functions for that reason.

We always use quadratic functions for the passive joint moments and those can get close enough to the experimental data as well (in my opinion), and simulation is much easier. Like the exponential model, there are three parameters for each direction. Parameters are linear stiffness k1, quadratic stiffness k2, and range of motion limit c.

Pseudo-code:

Code: Select all

M = -k1*theta;
if theta > c
   M = M - k2*(theta-c)^2
end     
(and similar code for the other direction) 
 

When the quadratic model and exponential model are both fitted to the same passive moment-angle data. I don't think the difference will be large enough to influence your simulations. The difference between the two models will only be large when theta goes far outside the range of motion limits, but that will never happen during real simulations. The passive moment is there to prevent that.

Ton van den Bogert

[Peter Bishop]

Hello Ton,
thanks for the very detailed and helpful reply!

Thanks also for the quadratic code suggestion. One thing though, if I use something like the ExpressionBasedBushingForce class, I don't think I'm allowed to use if-else statements when specifying the force in the model's .osim file; I think only mathematical expressions are permissible here (but I'll still give it a try).

Peter Bishop.

[Ton van den Bogert]

I have no experience doing forward dynamic simulation with Opensim, because I did this long before Opensim existed. It would be too bad if you can't use code to define a moment-angle relationship. Maybe someone else knows a way to do that with Opensim.

If you can't use if-then, maybe Opensim can use a spline function?

If all else fails, I suppose you can use brute force and keep reducing the max step size until the simulation no longer fails.

Ton