chore: update to ModelingToolkit.jl v11 (#279)

- docs: update code blocks to match new style using `get_u0`, `get_p`
- refactor!: Automatically infer variables and parameters in `System`
construction. No longer need to pass `u0` and `p` arrays explicitly.
- refactor!: remove `defaults` keyword dropped in MTK v11
- refactor!: use `scalarize` operations where appropriate to simplify
model definitions
- refactor: use `AstrodynamicalModels.CartesianSTM()` or `I(6)` in place
of literal matrices
- test: Use `MTK.get_u0`, `MTK.get_p` to correctly handle passing of
variable and parameter values to `ODEFunction`. Closes #270

Author: icweaver

Project.toml

@@ -17,7 +17,7 @@ SPICE = "5bab7191-041a-5c2e-a744-024b9c3a5062"
 
 [compat]
 AstrodynamicalCalculations = "0.5, 1"
-AstrodynamicalModels = "3.4"
+AstrodynamicalModels = "3.4, 3.8"
 AstrodynamicalSolvers = "0.5, 1"
 DocStringExtensions = "0.9"
 EphemerisSources = "1.0.0"

lib/AstrodynamicalModels/Project.toml

@@ -28,9 +28,9 @@ AstrodynamicalCalculations = "1"
 DocStringExtensions = "0.9"
 LinearAlgebra = "1"
 Memoize = "0.4"
-ModelingToolkit = "9"
+ModelingToolkit = "11"
 SPICEBodies = "1"
 SciMLBase = "2"
 StaticArrays = "1"
-Symbolics = "5, 6"
+Symbolics = "7"
 julia = "1.10"

lib/AstrodynamicalModels/docs/src/models/attitude/index.qmd

@@ -7,6 +7,7 @@ _Quaternion kinematics, and dynamics!_
 #| echo: false
 using AstrodynamicalModels
 using ModelingToolkit
+using ModelingToolkit: get_u0, get_p
 ```
 
 ## Overview
@@ -57,7 +58,10 @@ Finally, let's construct a Julia function which implements these dynamics!
 #| echo: true
 #| output: true
 f = AttitudeFunction()
-let u = randn(7), p = randn(15), t = 0
-    f(u, p, t)
+let u = randn(7), p = [randn(9), randn(3), rand(3)], t = 0
+    sys = f.sys
+    u0 = get_u0(sys, ModelingToolkit.unknowns(sys) .=> u)
+    p = get_p(sys, [:J => p[1], :L => p[2], :f => p[3]]) # Or get_p(sys, ModelingToolkit.parameters(sys) .=> p)
+    f(u0, p, t)
 end
 ```

lib/AstrodynamicalModels/docs/src/models/entry/index.qmd

@@ -7,6 +7,7 @@ _Also known as canonical entry dynamics!_
 #| output: false
 using AstrodynamicalModels
 using ModelingToolkit
+using ModelingToolkit: get_u0, get_p
 ```
 
 ## Overview
@@ -44,6 +45,9 @@ Finally, let's construct a Julia function which implements these dynamics!
 #| echo: true
 f = PlanarEntryFunction()
 let u = abs.(randn(4)), p = abs.(randn(7)), t = 0
-    f(u, p, t)
+    sys = f.sys
+    u0 = get_u0(sys, [:γ, :v, :r, :θ] .=> u) # Or get_u0(sys, ModelingToolkit.unknowns(sys) .=> u)
+    p = get_p(sys, ModelingToolkit.parameters(sys) .=> p)
+    f(u0, p, t)
 end
 ```

lib/AstrodynamicalModels/docs/src/models/nb/index.qmd

@@ -7,6 +7,7 @@ _Also known as NBP dynamics!_
 #| output: false
 using AstrodynamicalModels
 using ModelingToolkit
+using ModelingToolkit: get_u0, get_p
 ```
 
 ## Overview
@@ -49,7 +50,10 @@ Finally, let's construct a Julia function which implements these dynamics!
 ```{julia}
 #| echo: true
 f = NBFunction(2)
-let u = randn(12), m = randn(2), G = rand(), t = 0
-    f(u, [G, m...], t)
+let u = randn(12), m = randn(2), G = randn(), t = 0
+    sys = f.sys
+    u0 = get_u0(sys, ModelingToolkit.unknowns(sys) .=> u)
+    p = get_p(sys, [:m => m, :G => G]) # Or get_p(sys, ModelingToolkit.parameters(sys) .=> [m, G])
+    f(u0, p, t)
 end
 ```

lib/AstrodynamicalModels/docs/src/models/r2b/index.qmd

@@ -7,6 +7,7 @@ _Also known as R2BP dynamics!_
 #| output: false
 using AstrodynamicalModels
 using ModelingToolkit
+using ModelingToolkit: get_u0, get_p
 ```
 
 ## Overview
@@ -65,7 +66,10 @@ Finally, let's construct a Julia function which implements these dynamics!
 ```{julia}
 #| echo: true
 f = R2BFunction()
-let u = randn(6), p = [3e6], t = 0
-    f(u, p, t)
+let u = randn(6), p = 3e6, t = 0
+    sys = f.sys
+    u0 = get_u0(sys, [:x, :y, :z, :ẋ, :ẏ, :ż] .=> u) # Or get_u0(sys, ModelingToolkit.unknowns(sys) .=> u)
+    p = get_p(sys, [:μ => p]) # Or get_p(sys, ModelingToolkit.parameters(sys) .=> p)
+    f(u0, p, t)
 end
 ```

lib/AstrodynamicalModels/src/AstrodynamicalModels.jl

@@ -579,6 +579,14 @@ Return the parameter vector for an `Orbit`.
 """
 parameters(orbit::Orbit) = orbit.parameters
 
+"""
+Return the operating point for an `Orbit`.
+"""
+op(orbit::Orbit) = [
+    [:x, :y, :z, :ẋ, :ẏ, :ż] .=> orbit.state
+    :μ => orbit.parameters[1]
+]
+
 Base.getindex(orbit::AstrodynamicalOrbit, args...) = Base.getindex(state(orbit), args...)
 Base.setindex!(orbit::AstrodynamicalOrbit, args...) = Base.setindex!(state(orbit), args...)
 

lib/AstrodynamicalModels/src/Attitude.jl

@@ -167,26 +167,20 @@ model = Attitude()
 @memoize function AttitudeSystem(;
     stm = false,
     name = :Attitude,
-    defaults = Pair{ModelingToolkit.Num,<:Number}[],
     kwargs...,
 )
 
     @variables begin
         q(t)[1:4],   [description = "scalar-last attitude quaternion"]
         ω(t)[1:3],   [description = "body rates in radians per second"]
     end
+
     @parameters begin
         J[1:3, 1:3], [description = "moment of inertial matrix"]
         L[1:3],      [description = "lever arm where input torque is applied"]
         f[1:3],      [description = "input torque"]
     end
 
-    q = Symbolics.scalarize(q)
-    ω = Symbolics.scalarize(ω)
-    J = Symbolics.scalarize(J)
-    L = Symbolics.scalarize(L)
-    f = Symbolics.scalarize(f)
-
     A = [
             0  ω[3] -ω[2] ω[1]
         -ω[3]     0  ω[1] ω[2]
@@ -201,25 +195,15 @@ model = Attitude()
     ]
 
     eqs = [
-        D.(q) .~ (1 // 2) * (A * q) ;
-        D.(ω) .~ (-inv(J) * ωx * J * ω + inv(J) * L + f)
+        Symbolics.scalarize(D(q) ~ (1 // 2) * (A * q)) ;
+        Symbolics.scalarize(D(ω) ~ -inv(J) * ωx * (J * ω) + inv(J) * L + f)
     ]
 
-    u = [q; ω]
-
     if stm
-        @variables (Φ(t))[1:7, 1:7], [description = "state transition matrix estimate"]
-        A = Symbolics.jacobian(map(el -> el.rhs, eqs), u)
-
-        Φ = Symbolics.scalarize(Φ)
-
-        LHS = D.(Φ)
-        RHS = A * Φ
-
-        eqs = [eqs; vec([LHS[i] ~ RHS[i] for i in eachindex(LHS)])]
+        @variables Φ(t)[1:7, 1:7], [description = "state transition matrix estimate"]
 
-        u = [u; vec(Φ)]
-        defaults = [defaults; vec(Φ .=> Float64.(I(7)))]
+        A = Symbolics.jacobian(map(el -> el.rhs, eqs), vcat(q, ω))
+        eqs = vcat(eqs, vec(Symbolics.scalarize(D(Φ) ~ A * Φ)))
     end
 
     if string(name) == "Attitude" && stm
@@ -228,9 +212,8 @@ model = Attitude()
         modelname = name
     end
 
-    return System(eqs, t, u, [vec(J); L; f];
+    return System(eqs, t;
         name = modelname,
-        defaults,
         kwargs...,
     )
 end
@@ -248,15 +231,18 @@ are passed directly to `SciMLBase.ODEFunction`.
 
 ```julia
 f = AttitudeFunction()
-let u = randn(7), p = randn(15), t = NaN # time invariant
-    f(u, p, t)
+let u = randn(7), p = [randn(9), randn(3), rand(3)], t = NaN
+    sys = f.sys
+    u0 = get_u0(sys, ModelingToolkit.unknowns(sys) .=> u)
+    p = get_p(sys, [:J => p[1], :L => p[2], :f => p[3]]) # Or get_p(sys, ModelingToolkit.parameters(sys) .=> p)
+    f(u0, p, t)
 end
 ```
 """
 @memoize function AttitudeFunction(; stm = false, name = :Attitude, kwargs...)
     defaults = (; jac = true)
     options = merge(defaults, kwargs)
-    sys = complete(AttitudeSystem(; stm = stm, name = name); split = false)
+    sys = complete(AttitudeSystem(; stm, name); split = true)
     return ODEFunction{true,SciMLBase.FullSpecialize}(
         sys;
         options...,

lib/AstrodynamicalModels/src/CR3BP.jl

@@ -46,7 +46,6 @@ model = CR3BSystem(; stm=true)
 @memoize function CR3BSystem(;
     stm = false,
     name = :CR3B,
-    defaults = Pair{ModelingToolkit.Num,<:Number}[],
     kwargs...,
 )
     @parameters μ
@@ -73,20 +72,10 @@ model = CR3BSystem(; stm=true)
             );
     ]
 
-    u = [r; v]
-
     if stm
-        @variables (Φ(t))[1:6, 1:6], [description = "state transition matrix estimate"]
-        A = Symbolics.jacobian(map(el -> el.rhs, eqs), u)
-
-        Φ = Symbolics.scalarize(Φ)
-        LHS = D.(Φ)
-        RHS = A * Φ
-
-        eqs = [eqs; vec([LHS[i] ~ RHS[i] for i in eachindex(LHS)])]
-
-        u = [u; vec(Φ)]
-        defaults = [defaults; vec(Φ .=> Float64.(I(6)))]
+        @variables Φ(t)[1:6, 1:6], [description = "state transition matrix estimate"]
+        A = Symbolics.jacobian(map(el -> el.rhs, eqs), [r; v])
+        eqs = [eqs; vec(Symbolics.scalarize(D(Φ) ~ A * Φ))]
     end
 
     if string(name) == "CR3B" && stm
@@ -95,9 +84,8 @@ model = CR3BSystem(; stm=true)
         modelname = name
     end
 
-    return System(eqs, t, u, [μ];
+    return System(eqs, t;
         name = modelname,
-        defaults,
         kwargs...,
     )
 end
@@ -119,15 +107,18 @@ directly to `SciMLBase.ODEFunction`.
 
 ```julia
 f = CR3BFunction(; stm=false, jac=true)
-let u = randn(6), p = randn(1), t = 0
-    f(u, p, t)
+let u = randn(6), p = randn(), t = 0
+    sys = f.sys
+    u0 = get_u0(sys, [:x, :y, :z, :ẋ, :ẏ, :ż] .=> u) # Or get_u0(sys, ModelingToolkit.unknowns(sys) .=> u)
+    p = get_p(sys, [:μ => p]) # Or get_p(sys, ModelingToolkit.parameters(sys) .=> p)
+    f(u0, p, t)
 end
 ```
 """
 @memoize function CR3BFunction(; stm = false, name = :CR3B, kwargs...)
     defaults = (; jac = true)
     options = merge(defaults, kwargs)
-    sys = complete(CR3BSystem(; stm = stm, name = name); split = false)
+    sys = complete(CR3BSystem(; stm, name); split = true)
     return ODEFunction{true,SciMLBase.FullSpecialize}(
         sys;
         options...,
@@ -146,13 +137,12 @@ AstrodynamicalModels.CR3BOrbit(; state::AbstractVector, parameters::AbstractVect
 """
 Return an `ODEProblem` for the provided CR3B system.
 """
-CR3BProblem(u0, tspan, p; kwargs...) = ODEProblem(CR3BFunction(), u0, tspan, p; kwargs...)
+CR3BProblem(op, tspan; kwargs...) = ODEProblem(complete(CR3BSystem()), op, tspan; kwargs...)
 CR3BProblem(orbit::AstrodynamicalOrbit, tspan::Union{<:Tuple,<:AbstractArray}; kwargs...) =
     ODEProblem(
-        CR3BFunction(),
-        AstrodynamicalModels.state(orbit),
-        tspan,
-        AstrodynamicalModels.parameters(orbit);
+        complete(CR3BSystem()),
+        AstrodynamicalModels.op(orbit),
+        tspan;
         kwargs...,
     )
 CR3BProblem(orbit::AstrodynamicalOrbit, Δt; kwargs...) =

lib/AstrodynamicalModels/src/Entry.jl

@@ -77,7 +77,6 @@ model = PlanarEntrySystem()
 """
 @memoize function PlanarEntrySystem(;
     name = :PlanarEntry,
-    defaults = Pair{ModelingToolkit.Num,<:Number}[],
     kwargs...,
 )
 
@@ -114,9 +113,8 @@ model = PlanarEntrySystem()
         D(θ) ~ (v / r) * cos(γ),
     ]
 
-    model = System(eqs, t, [γ, v, r, θ], [R, P, H, m, A, C, μ];
-        name,
-        defaults,
+    model = System(eqs, t;
+        name = name,
         kwargs...,
     )
 
@@ -142,14 +140,17 @@ passed directly to `SciMLBase.ODEFunction`.
 ```julia
 f = PlanarEntryFunction()
 let u = randn(4), p = randn(7), t = NaN # time invariant
-    f(u, p, t)
+    sys = f.sys
+    u0 = get_u0(sys, [:γ, :v, :r, :θ] .=> u) # Or get_u0(sys, ModelingToolkit.unknowns(sys) .=> u)
+    p = get_p(sys, ModelingToolkit.parameters(sys) .=> p)
+    f(u0, p, t)
 end
 ```
 """
 @memoize function PlanarEntryFunction(; name = :PlanarEntry, kwargs...)
     defaults = (; jac = true)
     options = merge(defaults, kwargs)
-    sys = complete(PlanarEntrySystem(; name = name); split = false)
+    sys = complete(PlanarEntrySystem(; name); split = true)
     return ODEFunction{true,SciMLBase.FullSpecialize}(
         sys;
         options...,