simtk.org

We are currently using python to simulate the model.
-The model I am using is "arm26.osim".
-I am controlling the elbow joint angle by giving excitation as input to the model's biceps brachii muscle (BIClong,BICshort).
-It seems that the internal numerical calculation of the state variables of the model is incorrect when the excitation is given as an input that changes abruptly (and thus the angle changes abruptly).

Here are the details
-The attached image shows the step response when a step input is given as an example.
-"Step_Response_No_Bug" is when no bug occurs, and the simulation tick time is 1e-4 * 5 [sec].
-"Step_Response_Bug" is when a bug occurs, and the simulation is run at 1e-4 * 20 [sec] increments.
-"Step_Response_Comparison" is the comparison of the above two pictures.
-In both cases, a step input (see Axis 2, min0.01->max0.14) is given to the biceps brachii (BIClong,BICshort) after 1 second.
-The simulation is performed by the integrate() method of the Manager class. IntegratorMethod_RungeKuttaMerson" is specified as the integrator.

When I examined all available parameters during the simulation to investigate the cause of the bug, I found that "TendonVelocity" was affected by the bug before "Angle" (1.42[sec] for Angle and 1.36[sec] for TendonVelocity). If other parameters had bugs, they were all derived from "TendonVelocity", so "TendonVelocity" was the fastest parameter that could be output.
I checked the source code in c++ of the opensim API from "OpenSim Project - GitHub" and found that for a muscle with PennationAngle=0[rad], such as the biceps brachii, "TendonVelocity = MuscleVelocity - FiberVelocity". The same waveform for "FiberVelocity" was also checked, and no effect of the bug was observed in the same time range here. The waveform is shown in the attached image. However, it seems that an error occurred first in the numerical calculation of "MuscleVelocity," and that the error affected all variables starting from this error.
Here, I checked the c++ source code of "MuscleVelocity" as well, and it seems that the velocity is calculated directly from the displacement of the position coordinates of the musculoskeletal model, as shown below.

void GeometryPath::computeLengtheningSpeed(const SimTK::State& s) const
{
if (isCacheVariableValid(s, _speedCV)) {
return;
}

const Array<AbstractPathPoint*>& currentPath = getCurrentPath(s);

double speed = 0.0;

for (int i = 0; i < currentPath.getSize() - 1; i++) {
speed += currentPath->calcSpeedBetween(s, *currentPath[i+1]);
}

setLengtheningSpeed(s, speed);
}

In light of this, is there any way around this bug other than to make the time increments even shorter? Or is there another, more fundamental, fixable workaround?
With 1e-4 * 5 [sec] increments, for example, a 20-second control simulation would require about 45 seconds of computation.
I am currently running simulations on a scale of several hundred times in the course of my research, so this time-consuming and buggy process is very stressful for me.
I would appreciate any suggestions on how to solve this problem.
The program I am using is also stored in a zip file. The "python_code.zip\main\control_main.py" is the control program and the "python_code.zip\main\simulator.py" to run the simulation.
Thank you in advance for your help.

If the integrator is being run with fixed time steps, then there is no guarantee that the solution will meet a desired tolerance. You may wish to switch to an integration strategy that uses variable time steps, particularly if the excitation is discontinuous. The integrator can then take more steps when necessary to meet the specified tolerance.

Thank you for your reply, Mr. Uchida.
Indeed I was running the integrator with a fixed time step size. If I make this a variable time step size, might the numerical calculation work better, in which case it would automatically be even finer than the time step size I had set?
Also, I am using a fixed time step size for a calculation that requires a step size, such as PID control, etc. To avoid this bug while using a fixed time step size, is there no other way but to make the integration time even finer?

The integrator step size must be sufficiently small to capture the highest-frequency content in the response. To use fixed steps, you would need to either modify the system so its response doesn't change more quickly than the integrator can capture, or accept lower tolerances in the solution. You could also allow the integrator to take several steps (either fixed or variable) between controller updates.

Thank you, Mr. Uchida, for your numerous advices. I will consider the method again.