Activating muscles in every step of the simulation

[kidzik]

Hi, I am trying to activate muscles one by one every 0.01 second in simulation. It used to work fine, when I was compiling a version from some old commit b3163b22be7bd4cccd8c073e02dd3abb52867031 . In the beta release 399c8d57a779dd5dde2916192f8b92bfc959e269 it seems to only activate muscles in the first step of the simulation and they are deactivated in all the other steps. If I use setExcitation instead, it seems that nothing happens and muscles are deactivated throughout the entire simulation.

import opensim

model_path = "arm26.osim"
model = opensim.Model(model_path)
model.setUseVisualizer(True)
state = model.initSystem()
manager = opensim.Manager(model)
muscleSet = model.getMuscles()
stepsize = 0.01

for i in range(500):
    for j in range(6):
        muscleSet.get(j).setActivation(state, 1.0)
    manager.integrate(state, stepsize * (i + 1))

Is there anything I am missing? I checked arm26 and gait10dof18musc models. Thanks!

[chrisdembia]

The interface for Manager::integrate() changed. That could be it, but it looks like you're using the new interface correctly.

[junja94]

I'm having the same problem with similar code.

[chrisdembia]

@kidzik What do you mean in the first step? You mean for i = 0, or for the first time step the integrator takes within each call to manager.integrate()?

[kidzik]

@chrisdembia I meant i=0 - it basically blinks the muscles to red then and they are turned off throughout the rest of the simulation. I will check it more precisely on Monday (out of office now :). Thanks!

[chrisdembia]

Ah okay. My guess is that the state is not getting updated. That is state.getTime() == 0 even for i > 0.

[carmichaelong]

@chrisdembia I think at least part of the state is being updated. I wrote a simple case directly in C++. The initialization is about the same, and my loop looks like this:

    for (int i = 0; i < 10; ++i)
    {
        for (int j = 0; j < muscleSet.getSize(); ++j)
        {
            muscleSet.get(j).setActivation(state, 1.0);
        }
        cout << state.getTime() << " ";
        manager.integrate(state, stepsize*(i + 1));
        arm.realizeDynamics(state);
        cout << muscleSet.get(0).getActivation(state) << endl;
    }

The output is:

0 0.652489
0.01 0.468143
0.02 0.354615
0.03 0.278174
0.04 0.223565
0.05 0.182883
0.06 0.151619
0.07 0.127015
0.08 0.107286
0.09 0.0912298

[carmichaelong]

Worked on this some with @chrisdembia. Still need more investigation. Initial work and some test cases added in branch manager-initialize

[kidzik]

Just verified that I am getting the same values through the python interface

0.000000 0.652489
0.010000 0.468144
0.020000 0.354616
0.030000 0.278174
0.040000 0.223565
0.050000 0.182883
0.060000 0.151619
0.070000 0.127015
0.080000 0.107286
0.090000 0.091230

Back in b3163b22be7bd4cccd8c073e02dd3abb52867031 I'm getting

0.000000 0.655534
0.010000 0.655546
0.020000 0.655537
0.030000 0.655538
0.040000 0.655534
0.050000 0.655534
0.060000 0.655536
0.070000 0.655536
0.080000 0.655536
0.090000 0.655536

using the same code https://github.com/stanfordnmbl/osim-rl/blob/master/tests/test.manager.py

[kidzik]

Everything works well with @carmichaelong's suggestion (/hack) of recreating the manager

import opensim

model_path = "arm26.osim"
model = opensim.Model(model_path)
model.setUseVisualizer(True)
state = model.initSystem()
muscleSet = model.getMuscles()
stepsize = 0.01

for i in range(500):
    for j in range(6):
        muscleSet.get(j).setActivation(state, 1.0)
    manager = opensim.Manager(model)
    manager.integrate(state, stepsize * (i + 1))

Which is ok for my application. Thanks!

[jenhicks]

@carmichaelong will look into fixing issue.

[chrisdembia]

@carmichaelong I pushed some changes to the manager-initialize branch that uses Integrator::reinitialize().

[jimmyDunne]

Fix the interface so that Lukasz' bug wouldn't happen. Update timestepper in the future to facilitate.

[aymanhab]

@carmichaelong did investigation, will ping @sherm1 on best way to modify state between integration intervals.

[jenhicks]

I am moving to the next release unless you have updates @carmichaelong

[carmichaelong]

I think there's a bigger issue going on with caching that should definitely be pushed for next release (I have some tests that try to isolate this problem, but the outcome is inconclusive).

I think for this release, it could be important to do either or both of: 1) Change the signature of Manager::integrate(state, time) to Manager::integrate(time), so that users will not be able to change the state and pass back into the same Manager 2) Document that changes to the state passed back into Manager may not update as expected.

[jenhicks]

@carmichaelong will implement 1

[BeppeInfo]

@carmichaelong So at least in version 4.0 one won't be able to update muscle activations online/in real time ?

[BeppeInfo]

On a related issue (I guess), I'm trying to run a simulation in Python with the attached model, where I update the input for a CoordinateActuator on each pass. Segmentation faults happen after the first iteration, even without the actuator.setOverrideActuation command. what should I do ? right_knee_joint.osim.txt

 import sys
 import opensim

 model = opensim.Model( 'right_knee_joint.osim' )
 model.extendFinalizeFromProperties()

 systemState = model.initSystem()

 actuator = model.getComponent( 'knee_actuator' )
 actuator.overrideActuation( systemState, True )

 musclesList = model.updMuscles()
 for muscle in musclesList:
   muscle.setAppliesForce( systemState, False )

 timeDelta = 0.005
 for timeStep in range( 1, 100 ):

   actuator.setOverrideActuation( systemState, timeStep % 4 - 2 )
   manager = opensim.Manager( model )
   systemState.setTime( ( timeStep - 1 ) * timeDelta )
   manager.initialize( systemState )
   systemState = manager.integrate( timeStep * timeDelta )

 model.setUseVisualizer( False )

[chrisdembia]

Do you know which line is causing the segfault?

[BeppeInfo]

@chrisdembia The segfault happens the second time manager.initialized is called

[kidzik]

I think that with recent updates the initial issue here is not an issue anymore. Right now, I can run everything without hacks, initializing a manager only once:

manager = opensim.Manager(model)
state.setTime(0)
manager.initialize(state)

for i in range(100):
    state = manager.integrate(i * 0.01)

I don't know about @Bitiquinho's issue though

[BeppeInfo]

It seems like with lastest git I'm able to modify the actuator force during simulation without recreating the manager every step:

 import sys
 import opensim

 model = opensim.Model( 'right_knee_joint.osim' )
 model.extendFinalizeFromProperties()

 systemState = model.initSystem()

 actuator = model.getComponent( 'knee_actuator' )
 actuator.overrideActuation( systemState, True )

 musclesList = model.updMuscles()
 for muscle in musclesList:
   muscle.setAppliesForce( systemState, False )

 manager = opensim.Manager( model ) 
 systemState.setTime( 0 )
 manager.initialize( systemState )

 timeDelta = 0.005
 for timeStep in range( 1, 1000 ):
   actuator.setOverrideActuation( systemState, timeStep % 4 - 2 )
   systemState = manager.integrate( timeStep * timeDelta )

Reinitializing it every step still causes segfaults, though

[clnsmith]

@kidzik can you clarify what your final solution was to this problem? I am experiencing a similar issue trying to prescribe the activation for a muscle at each timestep of a forward simulation. Using the timestepper instead of manager so I can use the CPodesIntegrator.

//Set Integrator
SimTK::CPodesIntegrator integrator(_model.getSystem(), SimTK::CPodes::BDF, SimTK::CPodes::Newton);
integrator.setAccuracy(get_integrator_accuracy());

SimTK::TimeStepper timestepper(_model.getSystem(), integrator);
state.setTime(get_start_time());
timestepper.initialize(state);

//Integrate Forward in Time
double dt = get_maximum_time_step();
int nSteps = round((get_stop_time() - get_start_time()) / dt);
AnalysisSet& analysisSet = _model.updAnalysisSet();

SimTK::State& s= timestepper.updIntegrator().updAdvancedState();

for (int i = 0; i <= nSteps; ++i) {
    double t = get_start_time() + i * dt;
    std::cout << "Time:" << t << std::endl;

    //Set Prescribed Muscle Forces
    for (int j = 0; j < _prescribed_frc_actuator_paths.size();++j) {
        std::string actuator_path = _prescribed_frc_actuator_paths[j];
        ScalarActuator& actuator = _model.updComponent<ScalarActuator>(actuator_path);
        double value = _frc_functions.get(actuator_path +"_frc").calcValue(SimTK::Vector(1,t));
        actuator.setOverrideActuation(s, value);
    }

    //Set Prescribed Muscle Activations     
    for (int j = 0; j < _prescribed_act_actuator_paths.size();++j) {
        std::string actuator_path = _prescribed_act_actuator_paths[j];
        Muscle& actuator = _model.updComponent<Muscle>(actuator_path);
        double value = _act_functions.get(actuator_path + "_act").calcValue(SimTK::Vector(1,t));
        actuator.setActivation(s, value);
    }
    _model.realizeAcceleration(s);

    timestepper.stepTo(t);

    //Record parameters
    s=timestepper.updIntegrator().updAdvancedState();

    for (int j = 0; j < _prescribed_act_actuator_paths.size();++j) {
        std::string actuator_path = _prescribed_act_actuator_paths[j];
        Muscle& actuator = _model.updComponent<Muscle>(actuator_path);
        std::cout << actuator.getActivation(s) << " " << actuator.getActivation(s) << std::endl;
    }

This works for prescribing the force in a musle, but not for the activation. My work around right now is to use a Prescribed Controller and use ignote_activation_dynamics(true), but this doesnt work for Thelen2003 muscles.

[chrisdembia]

What about just using a sufficiently small time constant? 0.00001?

[kidzik]

@clnsmith In my case recreating the manager in every step of the simulation fixed the issue, by I don't know the underlying logic. @carmichaelong do you remember what was wrong here?

[clnsmith]

@chrisdembia Just tried with a time step of 0.0000001, and it still doesn't work.

As far as I can tell if you make changes to the state from: state=timestepper.updIntegrator().updAdvancedState();

They do not change the internal state used by the timestepper (Prehaps a difference between State and AdvancedState in the integrator? I didn't get that far digging down into Simbody)

So the state used in the call timestepper.initialize(state), is updated internally and will be used throughout a loop with multiple calls to timestepper.stepTo(time), regardless of how timestepper.updIntegrator().updAdvancedState(); is modified.

In my case I got it to work by calling timestepper.initialize(state) at every step in the loop. However, you get performance reduction because of the repeated initialization calls. Furthermore, for my scenario I just want a constant activation throughout the time step, so I would need to take really small steps to compensate for the activation dynamics.

So I will stick with my work around above ^. Looks like https://github.com/opensim-org/opensim-core/issues/366 will eventually address what I am looking for in a cleaner way.

[clnsmith]

hmm I forgot that: actuator.setOverrideActuation(state, value);

Did work without needing to reinitialize, so maybe this is more complex than my explanation. OverrideActuation sets a DiscreteVariable, whereas setActivation sets a StateVariable. So maybe that has something to do with what is going on.

[carmichaelong]

So the state used in the call timestepper.initialize(state), is updated internally and will be used throughout a loop with multiple calls to timestepper.stepTo(time), regardless of how timestepper.updIntegrator().updAdvancedState(); is modified.

In some of the testing while reworking OpenSim::Manager I think I also ran into this problem. I imagine something internal to the TimeStepper keeping track of the state between it in the Integrator makes it difficult to muck around in it as a user.

Furthermore, for my scenario I just want a constant activation throughout the time step, so I would need to take really small steps to compensate for the activation dynamics.

Is this true even if you set both the activation of the muscle AND the PrescribedController to the same constant value? How long are the time steps for each constant value (I assume this is the loop you're referring to that you have to initialize every time?)?

[carmichaelong]

Oh, also @clnsmith, you may have found this but a possible similar use case is in our tests: https://github.com/opensim-org/opensim-core/blob/master/OpenSim/Simulation/Test/testManager.cpp#L278

We similarly have to initialize every step, which seems to fall more in line with how TimeStepper's use is documented. There are some reinitialize methods for Integrators, but I couldn't get that working as intended.

[aseth1]

@clnsmith and @kidzik, if you are updating the activation every time-step then you are not in fact integrating the dynamics for that state variable. Your use of a prescribed controller and ignore_activation_dynamics in this case is the correct approach. My recommendation is to make the change to Thelen2003Muscle so that it correctly ignores activation dynamic. We would be happy to help you out with that. Otherwise, switch to the Millard2012EquilibriumMuscle.

[kidzik]

On my side, as far as I remember, it was just to test setActivations as we wanted to make sure that a certain level is set at the beginning of the simulation. Besides that we use excitations. So I think everything is clear now, thanks!