Add support to variables with custom type

[ronisbr]

Hi!

If I want to port a satellite simulator to MTK, I will need to have some variables that has custom types. For example, for the attitude, I might want to use quaternions that has a specific algebra for multiplication and propagation given an angular velocity. I know I can model all the algebra inside MTK, but it will not be optimal since we have packages that implement the quaternion algebra already.

I tried with this MWE:

using ModelingToolkit
using OrdinaryDiffEq
using ReferenceFrameRotations
import Base: zero

zero(::Type{Quaternion{T}}) where T = Quaternion{T}(0, 0, 0, 0)

@variables t x(t)
D = Differential(t)
@named quat_model = ODESystem(D(x) ~ Quaternion{Float64}(0, 0, 0, 0))
prob = ODEProblem(quat_model, [x => Quaternion(1.0, 0.0, 0.0, 0.0)], (0.0, 10.0))
sol = solve(prob, Tsit5())

However, this returns the following error:

ERROR: LoadError: setindex! not defined for Quaternion{Float64}

I just cannot set the method setindex! for Quaternion because it is immutable.

Is there any plan to support variables with custom types?

Note: The Quaternion of Quaternions.jl works because it is <: Number whereas Quaternion of ReferenceFrameRotations.jl is <:AbstractVector. For some specificities involving how the attitude is propagated, it is better to treat a quaternion as a 4x1 vector.

[YingboMa]

The support for custom types are there. Quaternion should be <: Number if it's a number.

[ronisbr]

For attitude propagation, we treat quaternions as a 4x1 vector instead of an extension of complex numbers.

[1] https://ntrs.nasa.gov/api/citations/20030093641/downloads/20030093641.pdf

[YingboMa]

Could you point me to the relavent part of the code?

[ronisbr]

Sorry, I think I did not understand. Which code?

[YingboMa]

Could you point me to the relavent part of the code that requires Quaternion <: AbstractVector?

[ronisbr]

oh, https://github.com/JuliaSpace/ReferenceFrameRotations.jl/blob/b6129b3032dc2b92104699d0fe7646e0ecf613bb/src/quaternion.jl#L703

This property is used when you propagate the attitude. In this case, you transform a quaternion into a vector, obtain the propagation matrix here:

https://github.com/JuliaSpace/ReferenceFrameRotations.jl/blob/b6129b3032dc2b92104699d0fe7646e0ecf613bb/src/quaternion.jl#L811

And then convert back to quaternion.

[ronisbr]

The code inside the simulator is:

    # Time derivative of the quaternion.
    dw.q_b_i = dquat(q_b_i, w_bi_b)

where the getters and setters in dw are defined as:

Base.@propagate_inbounds Base.getindex(a::Workspace, s::Val)                 = _getindex(a, s)

Base.@propagate_inbounds _getindex(a::Workspace{T}, ::Val{:q_b_i}) where T   = Quaternion(view(getfield(a, :x), 6:9))

Base.@propagate_inbounds Base.setindex!(a::Workspace, v, ::Val{:q_b_i})    = setindex!(getfield(a, :x), v, 6:9)

[ChrisRackauckas]

I think you're mixing up the mathematical vector with the computational vector. you want to treat the quaternion as a mathematical vector, but that does not mean that a mutable array is the right computational representation of it. An immutable number would be much faster (and safer) for this use case.

[ronisbr]

@ChrisRackauckas in my case, Quaternion <: AbstractVector but it is immutable. That's why I had that problem.

struct Quaternion{T} <: AbstractVector{T}
    q0::T
    q1::T
    q2::T
    q3::T
end

[ronisbr]

I think some part of the code saw that x is AbstractVector and assumed it is mutable, trying to call setindex!.

[ChrisRackauckas]

Yes, make it a StaticArrays.FieldVector?

[ronisbr]

Yes, make it a StaticArrays.FieldVector?

Hum, nice suggestion, I will try :) I think it will not break anything.

[YingboMa]

I have a fix locally.

[ronisbr]

@ChrisRackauckas I changed to Quaternion{T} <: FieldVector{4, T}, but I still get this error when trying to solve the system:

ERROR: setfield!: immutable struct of type Quaternion cannot be changed

meaning that setfield! is being called at some point.

[ronisbr]

I just noticed that the problem also happens if x is a SVector:

using StaticArrays
using ModelingToolkit
using OrdinaryDiffEq

ModelingToolkit.get_unit(x::SVector) = ModelingToolkit.unitless

@variables t x(t)
D = Differential(t)
@named quat_model = ODESystem(D(x) ~ SVector{4, Float64}(0.05, 0, 0, 0))
prob = ODEProblem(quat_model, [x => SVector{4, Float64}(1.0, 0.0, 0.0, 0.0)], (0.0, 10.0))
sol = solve(prob, Tsit5())

ERROR: LoadError: setindex!(::SVector{4, Float64}, value, ::Int) is not defined.