Does OpenSim report state equations?

[Sina Porsa]

Hi every one.
I was wondering if I can have the state equations of a model? either explicit or implicit.
It would be great if I could have the state equations in one of these forms:

dx/dt = f(x,u)

or

f(x,u,dx/dt)=0

where u is the excitation matrix.
to be more clear: these are the same equations that we integrate forward in time when using the "Forward Dynamic tool",

Thanks
Sina

[Sina Porsa]

Hi again,
It seems that I was not clear enough! I'll try to reword my question and I hope this time it will be crystal clear and let's hope someone knows the answer(fingers crossed!)

let's assume I know the states and excitations at a specific time t=t0, in other words, X(t=t0) and U(t=t0) are known. Is there any functions which returns the Xdot vector at that specific time t=t0? In other words: Is there any functions that reports Xdot(t=t0) ?

I searched the Doxygen document and I came up with these 3 functions. But I am not sure if these are reporting the values that I am looking for

state::getQDot () ----> for position variables q
state::getUDot () ----> for velocity variables u
state::getZDot () ----> for auxiliary variables z (Muscle fiber length and muscle activation)

I was wondering if these 3 function are reporting the values that I need? And if yes, how can I input the X and U vectors (states and excitations) into these functions?

Thanks
Sina

[Michael Sherman]

Hi, Sina. Your question is clear, and the answer is yes, but it is somewhat trickier than you might be thinking because you have to make use of the OpenSim API to establish correspondence between entities in the model and individual states and activations. Are you hoping to use the plant model in Matlab? I think if you provide a little more information about what you are trying to do we may be able to find a relevant example.

Regards,
Sherm

[Sina Porsa]

Hi Sherm and thank you for the reply.
I am using the OpenSim API to solve an optimal control problem and I am using an open source optimizer which is written in C++. (Everything is happening in C++). I will try my best to avoid a super technical question!
I am trying to solve a constraint problem which means I want to find an optimum X and U which minimizes this error function:

Err = F( X(t=t0) , U(t=t0) ) - Xdot(t=t0);

To avoid technical terms, I will not explain the F function here. Xdot(t=t0) is described in my previous post.

In other words I have 2 methods to calculate Xdot at time t = t0:

1- using these function: getQDot (), getUDot () and getZDot ()
2- using the F function.

Both of these methods will have the same input (X(t=t0) and U(t=t0)) and I am trying to minimize the error between the outputs.
My question is how can I chang the default values of X and U to (X(t=t0) and U(t=t0)) before using the getQDot() function, or how can I input (X(t=t0) and U(t=t0)) into getQDot()

Cheers
Sina

[Michael Sherman]

Hi, Sina. I'm confused now. Your error equation err=f(x,u)-xdot(t) doesn't specify arguments for xdot except time. Does that mean you are content to let OpenSim choose the state and activations? Or did you intend that OpenSim would calculate xdot based on your chosen values of x and u? Can you be precise about the unknowns, inputs, and outputs of the two functions you are comparing? There is no need to avoid a technical discussion if that would make it more clear.

If what you mean is err=f(t,x,u)-xdot(t,x,u) with the same t, x, and u on each side then the error would be minimized by f being identical to OpenSim's xdot calculation. I'm having trouble seeing that as an optimization problem involving a search for optimum x and u! Alternatively, if you have perhaps a simplified model in f with its own states x and activations u that you would like to find, then it isn't clear what the inputs to OpenSim are nor how the state derivatives could be matched up. So I am sure I don't actually understand what you're trying to do!

Regards,
Sherm

[Sina Porsa]

Hi Sherm.
In General, an optimal control problem is to find the optimum control u* which causes the system:
ydot(t) = f (y(t), u(t), t) (eq.1)
To follow the optimum trajectory y* that minimizes/maximizes a cost function.
In (eq. 1) y is the state vector, u is the control vector and t is the independent variable (usually time). In biomechanical studies the states vector (y) includes positions, velocities, muscle fibre lengths and muscle activations while the controls are the muscle excitations. Two well known methods to solve these kind of problems are Direct Shooting (DS) and Direct Collocation (DC).
In DS method, the controls are guessed and the state differential equations (eq. 1) are integrated explicitly. Anderson and Pandy used DS to find the optimal solution for normal gait.
In DC method both controls and states are guessed and the state differential equations (eq. 1) are defined as constraint functions. In other words this method tries to find the optimum y and u which are consistent with each other.
I am trying to use the DC method. As I guess both y and u in this method, I can calculate the left hand side of (eq. 1) based on my guess. I want to input these guesses to OpenSim to calculate the right hand side of (eq. 1). If the error between these two values is less than a threshold, I consider the guessed y and u consistent.
After this brief introduction, I have two questions:

1- Based on your previous post, SimTK::state:: getYDot () should return the right hand side of eq.1. I assume I should use this to input the guessed states
SimTK::state:: setY (const Vector &y)
The question is: how can I input the guessed u (guessed excitations)

2- As I was not sure what the output of SimTK::state:: getYDot () is, I tested it on the TugofWar model. Though the model has 16 states (6 general coordinates + 6 general velocities+ 2 muscle fibre lengths+ 2 muscle activations) the output of this function is a vector with 17 components which is weird.

Cheers
Sina

[Michael Sherman]

Hi, Sina. It's the correspondence between the model and the state and input variables that makes this tricky. There may also be constraint equations that must be satisfied by the states -- you would have to include constraint errors in your cost function or project any guessed states to the constraint manifold prior to evaluating the cost.

Luckily I managed to get hold of Ajay Seth (who is on vacation) -- he said he has an approach he uses for this in OpenSim and will respond.

Regards,
Sherm

[Michael Sherman]

P.S. I'm not certain but I think the reason you got 17 states for the tug of war is that general orientation requires a quaternion to avoid singularities. That is, a rigid body in space requires 7 generalized coordinates and 6 generalized speeds, plus a constraint on the quaternion length.

[Sina Porsa]

Hi again,
I reckon it shouldn't be that tricky, because the forward dynamic tool needs the same ydot vector in order to integrate the state equations forward in time. Something like y(t+h) = y(t) + h*ydot( t, y(t), u(t) ) which is a first order Euler method. Of course the FD tool uses more accurate methods, but all of them will need to calculate ydot( t, y(t), u(t) ) at some point.
Regards,
Sina

[Ajay Seth]

There are multiple ways you can set the individual state variables and the controls and evaluate the model dynamics to yield the system dynamics in first order form: ydot =f(y,u). I use the model component methods like setActivation and setFiberLength to set the state variables for muscles. You can do the same thing for Coordinates. The model component provides the context and components and their states are all accessible by name.

If you don't care about context and want to solve the system of equations directly, you can use Vector& y = state.updY() to get a writable vector. If you want more granularity you can get the vector of coordinates using updQ, mobilities updU, and additional states (typically activations and fiber lengths of muscles) with updZ(). Note that updating these will invalidate the system to different stages (Q position and above, U velocity and above, Z dynamics and above).

You knew this already it seems. The question then is how to update the controls? The controls are not state variables, but they are part of the system and they are owned my the Model (which is also a model component). The model has similar methods to the State access above, such as Vector& controls = model.updControls(State) that provides a writable controls vector. Note, that if you update the controls the dynamics stage is invalidated, since the dynamics of the system are (typically) dependent on the controls.

After you have updated the model state variables and controls you must realize the system to the acceleration stage to access Ydot from the state. E.g. model.getMultibodySystem().realize(state, Stage::Acceleration);

Note, if the system has kinematic constraints, the SimTK integrators employ a projection technique to ensure that the coordinates and speeds remain on the constraint manifold. If you are solving the functional simultaneously across multiple fixed time intervals (as in direct collocation) you need to enforce the kinematic constraints as well. As Sherm suggested, the easiest way is to include the constraint error in your cost function. You can get the Q and U errors from the state as well (once you are realized to the Velocity state or above).

Then you should be all set. Thank you for your patience.