As/mtk v9 (#63)

* refactor: update for MTKv9

* collect array broadcast

* run ci on latest release

---------

Co-authored-by: Aayush Sabharwal <[email protected]>

.github/workflows/docs.yml

@@ -16,7 +16,7 @@ jobs:
       - uses: actions/checkout@v2
       - uses: julia-actions/setup-julia@latest
         with:
-          version: '1.8'
+          version: '1'
       - name: Install dependencies
         run: julia --project=docs/ -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()'
       - name: Build and deploy

Project.toml

@@ -17,7 +17,7 @@ UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
 [compat]
 ControlSystemsBase = "1.0.1"
 DataInterpolations = "3, 4"
-ModelingToolkit = "8.49"
+ModelingToolkit = "9.0"
 ModelingToolkitStandardLibrary = "2"
 MonteCarloMeasurements = "1.1"
 OrdinaryDiffEq = "6"

docs/src/index.md

@@ -178,7 +178,7 @@ That's pretty cool, but even nicer is to generate some code for this symbolic sy
 defs = ModelingToolkit.defaults(simplified_sys)
 x, pars = ModelingToolkit.get_u0_p(simplified_sys, defs, defs) # Extract the default state and parameter values
 
-fun = Symbolics.build_function(symbolic_sys, states(simplified_sys), ModelingToolkit.parameters(simplified_sys);
+fun = Symbolics.build_function(symbolic_sys, unknowns(simplified_sys), ModelingToolkit.parameters(simplified_sys);
     expression = Val{false}, # Generate a compiled function rather than a Julia expression
     force_SA   = true,       # Use static arrays instead of regular arrays, for higher performance 🚀
 )

src/ControlSystemsMTK.jl

@@ -21,15 +21,15 @@ using ControlSystemsBase: ssdata, AbstractStateSpace, Continuous, nstates, noutp
 # using ControlSystemIdentification
 using RobustAndOptimalControl, MonteCarloMeasurements
 import ModelingToolkit: ODESystem, FnType, Symbolics
-using ModelingToolkit: states, observed, isdifferential
+using ModelingToolkit: unknowns, observed, isdifferential
 using Symbolics
 using Symbolics: jacobian, solve_for
 using UnPack
 # using Optim, Optim.LineSearches
 
 # using SymbolicControlSystems
 
-export sconnect, feedback, ODESystem, states, observed, named_ss
+export sconnect, feedback, ODESystem, unknowns, observed, named_ss
 export batch_ss, trajectory_ss, GainScheduledStateSpace
 export build_quadratic_cost_matrix
 

src/ode_system.jl

@@ -176,7 +176,7 @@ function RobustAndOptimalControl.named_ss(
         inputs = map(inputs) do inp
             if inp isa ODESystem
                 @variables u(t)
-                if u ∈ Set(states(inp))
+                if u ∈ Set(unknowns(inp))
                     inp.u
                 else
                     error("Input $(inp.name) is an ODESystem and not a variable")
@@ -193,7 +193,7 @@ function RobustAndOptimalControl.named_ss(
         outputs = map(outputs) do out
             if out isa ODESystem
                 @variables u(t)
-                if u ∈ Set(states(out))
+                if u ∈ Set(unknowns(out))
                     out.u
                 else
                     error("Outut $(out.name) is an ODESystem and not a variable")
@@ -231,7 +231,7 @@ function RobustAndOptimalControl.named_ss(
     end
     named_ss(
         lsys;
-        x = symstr.(states(ssys)),
+        x = symstr.(unknowns(ssys)),
         u = unames,
         y = symstr.(outputs),
         name = string(Base.nameof(sys)),
@@ -284,7 +284,7 @@ function named_sensitivity_function(
         inputs = map(inputs) do inp
             if inp isa ODESystem
                 @variables u(t)
-                if u ∈ Set(states(inp))
+                if u ∈ Set(unknowns(inp))
                     inp.u
                 else
                     error("Input $(inp.name) is an ODESystem and not a variable")
@@ -322,7 +322,7 @@ function named_sensitivity_function(
     end
     named_ss(
         lsys;
-        x = symstr.(states(ssys)),
+        x = symstr.(unknowns(ssys)),
         u = unames,
         y = unames, #Symbol.("out_" .* string.(inputs)),
         name = string(Base.nameof(sys)),
@@ -353,7 +353,7 @@ The second problem above, the ordering of the states, can be worked around using
 - `costs`: A vector of pairs
 """
 function build_quadratic_cost_matrix(matrices::NamedTuple, ssys::ODESystem, costs::AbstractVector{<:Pair})
-    x = ModelingToolkit.states(ssys)
+    x = ModelingToolkit.unknowns(ssys)
     y = ModelingToolkit.outputs(ssys)
     nx = length(x)
     new_Cs = map(costs) do (xi, ci)
@@ -388,7 +388,7 @@ The second problem above, the ordering of the states, can be worked around using
 """
 function build_quadratic_cost_matrix(sys::ODESystem, inputs::AbstractVector, costs::AbstractVector{<:Pair}; kwargs...)
     matrices, ssys = ModelingToolkit.linearize(sys, inputs, first.(costs); kwargs...)
-    x = ModelingToolkit.states(ssys)
+    x = ModelingToolkit.unknowns(ssys)
     y = ModelingToolkit.outputs(ssys)
     nx = length(x)
     new_Cs = map(costs) do (xi, ci)
@@ -446,7 +446,7 @@ eqs = [D(x) ~ v
 
 @named duffing = ODESystem(eqs, t)
 
-bounds = getbounds(duffing, states(duffing))
+bounds = getbounds(duffing, unknowns(duffing))
 sample_within_bounds((l, u)) = (u - l) * rand() + l
 # Create a vector of operating points
 ops = map(1:N) do i
@@ -513,7 +513,7 @@ function trajectory_ss(sys, inputs, outputs, sol; t = _max_100(sol.t), allow_inp
     maximum(t) > maximum(sol.t) && @warn("The maximum time in `t`: $(maximum(t)), is larger than the maximum time in `sol.t`: $(maximum(sol.t)).")
     minimum(t) < minimum(sol.t) && @warn("The minimum time in `t`: $(minimum(t)), is smaller than the minimum time in `sol.t`: $(minimum(sol.t)).")
     lin_fun, ssys = linearization_function(sys, inputs, outputs; kwargs...)
-    x = states(ssys)
+    x = unknowns(ssys)
     defs = ModelingToolkit.defaults(sys)
     ops = map(t) do ti
         Dict(x => robust_sol_getindex(sol, ti, x, defs; verbose) for x in x)
@@ -680,7 +680,7 @@ function GainScheduledStateSpace(systems, vt; interpolator, x = zeros(systems[1]
         [Differential(t)(x[i]) ~ sum(A[i, k] * x[k] for k in 1:nx) +
                                  sum(B[i, j] * (input.u[j] - u0[j]) for j in 1:nu)
          for i in 1:nx];
-        output.u .~ C * x .+ D * (input.u .- u0) .+ y0
+        collect(output.u .~ C * x .+ D * (input.u .- u0) .+ y0)
     ]
     compose(ODESystem(eqs, t, name = name), [input, output, scheduling_input])
 end

test/test_ODESystem.jl

@@ -19,8 +19,8 @@ C0 = pid(1, 1) * tf(1, [0.01, 1]) |> ss
 fb              = feedback(loopgain, name = :fb) * ref
 fb              = structural_simplify(fb)
 
-@test length(states(P)) == 3 # 1 + u + y
-@test length(states(C)) == 4 # 2 + u + y
+@test length(unknowns(P)) == 3 # 1 + u + y
+@test length(unknowns(C)) == 4 # 2 + u + y
 
 x0 = Pair[loopgain.P.x[1]=>1]
 
@@ -60,7 +60,7 @@ P1 = ss(Pc, [Pc.input.u], [Pc.output.u])
 
 
 # === Go the other way, ODESystem -> StateSpace ================================
-x = states(P) # I haven't figured out a good way to access states, so this is a bit manual and ugly
+x = unknowns(P) # I haven't figured out a good way to access states, so this is a bit manual and ugly
 @unpack input, output = P
 P02_named = named_ss(P, [input.u], [output.u])
 @test P02_named.x == [Symbol("(x(t))[1]")]
@@ -71,14 +71,14 @@ P02 = ss(P02_named)
 @test P02 == P0 # verify that we get back what we started with
 
 # same for controller
-x = states(C)
+x = unknowns(C)
 @nonamespace C02 = named_ss(C, [C.input], [C.output])
 @test ss(C02) == C0
 
 
 ## Back again for a complete round trip, test that ODESystem get correct names
 @named P2 = ODESystem(P02_named)
-@test Set(states(P2)) == Set(states(P))
+@test Set(unknowns(P2)) == Set(unknowns(P))
 @test Set(ModelingToolkit.inputs(P2)) == Set(ModelingToolkit.inputs(P))
 @test Set(ModelingToolkit.outputs(P2)) == Set(ModelingToolkit.outputs(P))
 
@@ -165,7 +165,7 @@ prob = ODEProblem(simplified_sys, x0, (0.0, 20.0), p, jac = true)
 sol = solve(prob, OrdinaryDiffEq.Rodas5(), saveat = 0:0.01:20)
 if isinteractive()
     @show sol.retcode
-    plot(sol, layout = length(states(simplified_sys)) + 1)
+    plot(sol, layout = length(unknowns(simplified_sys)) + 1)
     plot!(sol.t, sol[P.x[1]] - sol[P.x[3]], sp = 12, lab = "Δq")
 
     ##
@@ -267,4 +267,4 @@ Sn = get_named_sensitivity(sys_outer, [:inner_plant_input, :inner_plant_output])
 
 @test S == Sn.sys
 
-@test Sn.u == Sn.y == [:inner_plant_input, :inner_plant_output]
+@test Sn.u == Sn.y == [:inner_plant_input, :inner_plant_output]

test/test_batchlin.jl

@@ -18,7 +18,7 @@ eqs = [D(x) ~ v
 
 @named duffing = ODESystem(eqs, t, systems=[y, u])
 
-bounds = getbounds(duffing, states(duffing))
+bounds = getbounds(duffing, unknowns(duffing))
 sample_within_bounds((l, u)) = (u - l) * rand() + l
 # Create a vector of operating points
 N = 10