Moco goals for average speed in cycling model

[Nicholas Bianco]

Hi Ana,

2) Nicholas, is there a significant difference between using MocoScaleFactor-addScaleFactor and simply adding a different constant value as a reference?

The main difference is that the MocoScaleFactor is a parameter that can be optimized during the simulation. So if you want the simulation to discover the optimal constant crank speed as part of the cycling simulation, that could be one way to go.

As for the crank resistance, do you need the resistance force to be active. It might be simpler to use a SpringGeneralizedForce to the coordinate, which models linear stiffness and damping, or an ExpressionBasedCoordinateForce if you want more control over the stiffness/damping model (i.e., non-linear springs/dampers).

Best,
Nick

[Ana de Sousa]

Hi Nicolas!

Thanks a lot for the idea! Using your suggestion, I implemented the expression-based coordinate force, and it worked like a charm!

Code: Select all

expression = osim.ExpressionBasedCoordinateForce("crank_angle", "-200*qdot")
model.addComponent(expression)

Now, it is just a matter of modelling this properly

I have a follow-up question: why the prescribed force doesn't work? You mentioned the need for the force to be active. Does this imply that the prescribed force only operates in forward dynamics and not in other modes utilised by Moco? If there's documentation on this, I'd greatly appreciate a pointer. I want to avoid more stupid questions like this in the future and minimise inconveniences for you guys.

Thank you once again for your guidance!

[Nicholas Bianco]

Hi Ana,

A PrescribedForce can work, it's just more difficult to model correctly. Also, when you said "resistance" I assumed you might want a passive force in the model. I say "passive" in the sense that the force expression is based on the model configuration and not an applied force based on a function w.r.t. time.

PrescribedForce requires OpenSim Function's (i.e., functions w.r.t. time) to define the forces and torques applied to the model. In addition, PrescribedForce applies body forces to the model, which means you need to define the point of force application and the force and torque vectors, either expressed in ground or in a body frame. It takes careful modeling to do this correctly. On the other hand, ExpressionBasedCoordinateForce applies a generalized force to the crank joint (i.e., a portion of the net joint moment for that joint), and doesn't require you to specify how the spatial force is defined.

I want to avoid more stupid questions like this

Definitely not a stupid question! Modeling forces correctly can be one of the trickier aspects of any simulation pipeline. And I guarantee someone else will find this post useful in the future. So keep them coming

Best,
Nick

[Ana de Sousa]

Hi again, Nick, thanks a lot for the motivation!

Your explanation was really great, and I'm starting to wrap my head around it. I reckon what I'm after is a generalised force applied at a joint (crank). But here's the hiccup I've run into:

I managed to model crank resistance based on velocity, where higher velocity means more resistance:

Code: Select all

expression = osim.ExpressionBasedCoordinateForce("crank_angle", "-200*qdot")

It seems to be working. By "working," I mean the forces needed to get the crank moving are higher, and depending on some factors, it also takes longer for the model to move from 0 to 360.

Now, here's where it gets a bit tricky. Physically, the resistance torque on the crank shouldn't really depend on velocity; my first modelling would be more like a constant force (as friction). But when I try this:

Code: Select all

expression = osim.ExpressionBasedCoordinateForce("crank_angle", "-200")

Nothing changes. Whether there's resistance or not, the torque and final time stay the same.

I'm guessing the constant part doesn't play nice in the expression. So, I have two questions:

1) How can this constant part work?
2) Would the constant part (-100) in this expression also be ignored?

Code: Select all

expression = osim.ExpressionBasedCoordinateForce("crank_angle", "-100-200*qdot")

Thanks and cheers!
Ana

[Ana de Sousa]

Oops, my bad. I goofed up in my last post.

The problem I'm solving: the crank model has only the crank (no human) and an actuator. I asked Moco to solve the required torque to move the crank from 0 to 180, minimising the time.

What I get is the following (check out the torque plots here)

Without any ExpressionBasedCoordinateForce: T_final = 0.63 (positive torque)

With 200: T_final = 0.061 (lower than without, negative torque)

With -200: T_final = 0.061 (lower than without, positive torque)

With 200*qdot: T_final = 0.736 (higher than without, no torque)

With -200*qdot: T_final = 0.736 (higher than without, no torque)

Doubts
I have several other questions:

1. Why am I getting zero torque for 200qdot and -200qdot?
2. Is the torque expressed in Newtons? It seems odd that a constant torque of -200 would result in a small difference between with and without resistance...
3. Oddly enough, solving with resistance takes less time than without. That's not what I expected.
4. It is also a bit awkward that the times are similar between 200 and -200, even though the torques are symmetrically inverse. Am I missing something here?


The Moco study for all of them is the same:

Code: Select all

study = osim.MocoStudy()
problem = study.updProblem()
problem.setModel(model)
control_goal = osim.MocoFinalTimeGoal('final_time', 1)
problem.addGoal(control_goal)
problem.setTimeBounds(start_time, [0, end_time])
problem.setStateInfo('/jointset/crank_angle/crank_angle/value',
                         [low_bound, high_gound],
                         start_position,
                         end_position
                         )
problem.setStateInfo('/jointset/crank_angle/crank_angle/speed',
                         [-100, 0],
                         [],
                         []
                         )
problem.setControlInfo('/tau_crank_angle',
                           [-50, 50]
                           )
solver = study.initCasADiSolver()
solver.set_num_mesh_intervals(25)
solver.set_verbosity(2)
solver.set_optim_convergence_tolerance(1e-2)
solver.set_optim_constraint_tolerance(1e-2)

Cheers!

[Nicholas Bianco]

Hi Ana,

I believe that constant expressions should work, but I've never directly verified it myself. It is possible that expressions excluding "q" or "qdot" get ignored. If you want to model friction, then I'm not sure this is the way to go (nor do we have a readily available component for applying generalized friction forces). You'd need something that can model the stick-slip transition and applies friction force as a function of normal force (at least the latter part).

1. Why am I getting zero torque for 200qdot and -200qdot?

If the crank is allowed to have initial velocity, I wonder if the expression force is "taking the crank for a ride" here, where the actuator doesn't need to apply any force at all.

2. Is the torque expressed in Newtons? It seems odd that a constant torque of -200 would result in a small difference between with and without resistance...

Torque is in Newton-meters (N-m).

3. Oddly enough, solving with resistance takes less time than without. That's not what I expected.

Similar to my answer for 1, the constant force might be assisting the actuator move the crank. Have you looked at the difference in the crank actuator force between the no resistance and resistance cases?

4. It is also a bit awkward that the times are similar between 200 and -200, even though the torques are symmetrically inverse. Am I missing something here?

This doesn't seem odd to me. I'd be more concerned if the times were different.

-Nick

[Ross Miller]

Hi Ana,

The following paper might be helpful on modeling the crank mechanics:

https://www.researchgate.net/profile/Be ... dation.pdf

I've never modeled/simulated pedaling myself but I see this paper and BJ Fregly's other related papers cited a lot in pedaling models, specifically on modeling the crank mechanics.

Ross

[Ana de Sousa]

Hey guys!

Just wanted to drop a quick update and a big thank you for all the support. After diving into some simulations, I found out a couple of things that I was doing incorrectly. I'll post it here in case someone else also does the same.


To Ross
Yes, this paper is gold, thank you. It's on my to-do list to work that crank modelling you mentioned. First, I wanted to make sure all tools were working and that I understood all of them.

To Nick

I believe that constant expressions should work, but I've never directly verified it myself.

Turns out, the constant was innocent, my bad. You can see what I was doing wrong below and let me know if that was expected to happen or if is it like a "glitch" in Moco.

Torque is in Newton-meters (N-m)

I needed to confirm this because I'm wrestling with convergence issues. High values do not converge, but a value of 10N-m doesn't seem to make a difference (I thought it should, but I could be wrong). Got more tests in the pipeline before I circle back to this.

Issue alert

If the crank is allowed to have initial velocity, I wonder if the expression force is "taking the crank for a ride" here, where the actuator doesn't need to apply any force at all.

Nick's comment made me realize my initial speed wasn't starting at zero.

OLD) No go
I set the speed boundaries from -100 to 0 (negative to prevent backward movement). Then, I threw in setStateInfoPattern, expecting it to define the initial speed as well. Well, it doesn't.

Code: Select all

problem.setStateInfo('/jointset/crank_angle/crank_angle/speed',
                         [-100, 0],
                         [],
                         []
                         )
problem.setStateInfoPattern('/jointset/.*/.*/speed', [], 0.0, [])

Moco showed me only [-100, 0] was set.

Code: Select all

/jointset/crank_angle/crank_angle/speed. bounds: [-100, 0]

NEW) It's a go
I changed it here:

Code: Select all

problem.setStateInfo('/jointset/crank_angle/crank_angle/speed',
                         [-100, 0],
                         0,
                         []
                         )
problem.setStateInfoPattern('/jointset/.*/.*/speed', [], 0.0, [])

Now Moco starts with the right values.

Code: Select all

/jointset/crank_angle/crank_angle/speed. bounds: [-100, 0] initial: 0

[Carlos Gonçalves]

Hello Ana,

It might be worth checking out if there is a difference when you set the overall coordinate's speed as below.

Code: Select all

problem.setStateInfoPattern('/jointset/.*/speed', [-50, 50], 0)

It usually works for me.

Checking the Moco output for every little bit of state and control is always a good procedure when things stop making sense

Best regards.

[Ana de Sousa]

Hi Carlos, appreciate the quick response!

Do you mean

Code: Select all

/.*/

instead of

Code: Select all

/.*/.*/

? I tested that. It behaves the same way, i.e., defining the initial speed incorrectly.

One thing I missed mentioning, and I suspect it might be the reason: initially, I assumed it would take precedence without prior definitions, but it seems to only come into play when values are left undefined.

Cheers!