I didn’t see an obvious solution in the docs so I thought I would ask: I have an ODE problem where the state is a different size than the change in state.
Specifically, lets say my state of [pos; vel] includes a quaternion of size 4, for as the orientation of some object, and a 3 vector for the rotational velocity of the object. The change in state would be [vel, accel], both terms being 3 vectors. Thus size(state) == 7, size(state’) == 6. Overall there are 3 degrees of freedom of the system, but I’m representing orientation with a quat.
This is managed if one defines the integration step to specifically integrate the quaternion with the velocity components, but I’m unsure how to do this within DiffEq setting / ODEProblem. Any pointers would be excellent!
The bigger picture is that I’m using a physics engine under the ODE solver, and specifying quats is more stable. I have everything else working (and nicely differentiating) except the quaternions
You will still have a du of the same size. The quaternion update part produces a quaternion and the rotational rate, translational position, and translational velocity update parts produce 3-vectors. So nothing should be changing sizes.
Maybe you were referring to using Euler angles inside your function? If so, you can just convert to Euler angles or (usually better) a direction cosine matrix inside your function for when you need it. But that doesn’t change the fact that you’ll still want to do the derivative update on the quaternion form.
Yes the quaternion update part produces a quaternion, but from a 3d vector. the state’ variable would have the 3d vector while the next state needs to be updated to a 4d quaternion. The nominal form for DifferentialEquations.jl requires a function du = f(u, p, t); u would be the quaternion, and du the 3d vector of angular velocity. The integration step is what maps du to u at then next timestep. From what I can gather, the integration step is generally u = u + dt * du. I would need a way to specify this integration step specifically for the quat situation, I believe.
du for the rotational part should include your angular rates and orientation. So that should be 7 states just for the rotational part and another 6 for the translational part, giving you 13 states in total.
In Systems Theory, the state is the “collection” of minimal length for which the initial value uniquely determines the solution. So if your model is a semi-explicit DAE of index 0 or 1, and of form:
\dot{x} = f(x,y;t) \\
0 = g(x,y;t)
and you collect your unknowns in u = (x,y), then the state is x and not u, since knowledge of x(0) uniquely determines the solution and the length of x is less than the length of u.
Perhaps the concept of “state” is used differently in other fields.
If your problem is a semi-explicit DAE of index 0 or 1, using the DAE solvers of DifferentialEquations.jl should work well.
I haven’t tested the DAE solvers for more general DAEs of form F(\dot{x},x;t)= 0 and of higher index, but I would guess they should work also in that case (using Pantelides’ algorithm for index reduction?).
That’s an accurate summary of the problem setting, and we are using state similarly but the complication is the use of quaternions. As you say,
When rotations are concerned, a particular orientation can be represented in multiple ways with a 3d vector as Euler angles, whereas a quaternion unique represents the given orientation. The underlying system, however, only has three degrees of freedom, so the change in state is still a 3d vector. Said another way, differentiating a quaternion returns a 3d vector.
My system should indeed be representable by a DAE but I’ll have to think a little bit on how to do it here. The complication is that my physics engine is actually a C library (MuJoCo) so I don’t have the algebraic functions off the top of my head.
I don’t have the necessary background in physics to appreciate quaternions, etc. That’s the reason for my somewhat vague response.
However, if your system really is a DAE, it is highly recommended to not differentiate the algebraic equation in order to produce a “fake” ODE. Thus, it is recommended to avoid introducing (assuming index 0 or 1):
@ChrisRackauckas Is there a way to just apply a callback to adjust the state of the solver after every iteration? It would be sufficient to simply renormalize the quaternion after every step of the solver.
After some thinking, I can re-state my question a little more clearly. I can represent my change in (quaternion) state as a quaternion as well; the sizes would be the same.
However integration becomes tricky. Applying a quaternion to another requires a quaternion multiplication (and a re-normalization if necessary). I assume most solvers in DiffEq would do a u = u + dt * du step, but if u is a quat, this step would need to be hijacked to perform a quaternion multiplication on the correct indices of u.
For the project I’m working on, we can’t move to DAEs quite yet, but is there a way to inject a slightly different state update function? I can presumably use a ManifoldProjection to renormalize the quaternion, but that would be after the update. While locally the linear approximation is fine, it would be swell to use the correct quaternion math is possible.
