test: update tests according to new alias elimination

Author: AayushSabharwal

lib/ModelingToolkitBase/test/analysis_points.jl

@@ -1,6 +1,7 @@
 using ModelingToolkitBase, ModelingToolkitStandardLibrary.Blocks, ControlSystemsBase
 using ModelingToolkitStandardLibrary.Mechanical.Rotational
 using ModelingToolkitStandardLibrary.Blocks
+using SymbolicIndexingInterface
 using OrdinaryDiffEq, LinearAlgebra
 using Test
 using ModelingToolkitBase: t_nounits as t, D_nounits as D, AnalysisPoint, AbstractSystem
@@ -618,6 +619,16 @@ if @isdefined(ModelingToolkit)
     @test SciMLBase.successful_retcode(solve(prob, Rodas5P()))
     matrices_normal, _ = get_sensitivity(sys_normal, sys_normal.normal_inner.ap)
 
+    function compare_matrices(reference, value)
+        @assert size(reference.A) == (2, 2) "This testing function is only valid for `2x2` systems"
+        @test isequal(reference.C, value.C) || isequal(reverse(reference.C), value.C)
+        colorder = isequal(reference.C, value.C) ? [1, 2] : [2, 1]
+        @test isequal(reference.B, value.B) || isequal(reverse(reference.B), value.B)
+        roworder = isequal(reference.B, value.B) ? [1, 2] : [2, 1]
+        @test isequal(reference.D, value.D)
+        @test isequal(reference.A[roworder, colorder], value.A)
+    end
+
     @testset "Analysis point overriding part of connection - normal connect" begin
         @named F1 = FirstOrder(k = 1, T = 1)
         @named F2 = FirstOrder(k = 1, T = 1)
@@ -643,7 +654,7 @@ if @isdefined(ModelingToolkit)
         @test SciMLBase.successful_retcode(solve(prob, Rodas5P()))
 
         matrices, _ = get_sensitivity(sys, sys.ap)
-        @test matrices == matrices_normal
+        compare_matrices(matrices_normal, matrices)
     end
 
     @testset "Analysis point overriding part of connection - variable connect" begin
@@ -671,7 +682,7 @@ if @isdefined(ModelingToolkit)
         @test SciMLBase.successful_retcode(solve(prob, Rodas5P()))
 
         matrices, _ = get_sensitivity(sys, sys.ap)
-        @test matrices == matrices_normal
+        compare_matrices(matrices_normal, matrices)
     end
 
     @testset "Analysis point overriding part of connection - mixed connect" begin
@@ -699,7 +710,7 @@ if @isdefined(ModelingToolkit)
         @test SciMLBase.successful_retcode(solve(prob, Rodas5P()))
 
         matrices, _ = get_sensitivity(sys, sys.ap)
-        @test matrices == matrices_normal
+        compare_matrices(matrices_normal, matrices)
     end
 
     @testset "Ignored analysis points only affect relevant connection sets" begin

lib/ModelingToolkitBase/test/initial_values.jl

@@ -133,7 +133,8 @@ end
 end
 
 @testset "Cyclic dependency checking and substitution limits" begin
-    @variables x(t) y(t)
+    # `irreducible` to make sure they are unknowns of the initialization system
+    @variables x(t) [irreducible = true] y(t) [irreducible = true]
     @mtkcompile sys = System(
         [D(x) ~ x, D(y) ~ y], t; initialization_eqs = [x ~ 2y + 3, y ~ 2x],
         guesses = [x => 2y, y => 2x]

lib/ModelingToolkitBase/test/initializationsystem.jl

@@ -331,7 +331,7 @@ conditions = getfield.(equations(initprob.f.sys), :rhs)
 if @isdefined(ModelingToolkit)
     initsol = solve(initprob, reltol = 1.0e-12, abstol = 1.0e-12)
     @test SciMLBase.successful_retcode(initsol)
-    @test maximum(abs.(initsol[conditions])) < 1.0e-8
+    @test maximum(abs.(initsol[conditions])) < 2.0e-8
 end
 
 @test_throws ERRMOD.ExtraEquationsSystemException ModelingToolkitBase.InitializationProblem(
@@ -1343,7 +1343,7 @@ if @isdefined(ModelingToolkit)
         model = dc_motor()
         sys = mtkcompile(model)
 
-        prob = ODEProblem(sys, [sys.L1.i => 0.0], (0, 6.0))
+        prob = ODEProblem(sys, [sys.L1.i => 0.0], (0, 6.0); guesses = [sys.emf.flange.phi => 0])
 
         @test_nowarn remake(prob, p = prob.p)
     end

lib/ModelingToolkitBase/test/input_output_handling.jl

@@ -141,12 +141,16 @@ if @isdefined(ModelingToolkit)
     model_inputs = [torque.tau.u]
     op = Dict(torque.tau.u => 0.0)
     matrices, ssys = linearize(
-        model, model_inputs, model_outputs; op
+        model, model_inputs, model_outputs; op,
+        guesses = [inertia2.flange_a.phi => 0.0, inertia1.flange_b.phi => 0.0]
     )
     @test length(ModelingToolkit.outputs(ssys)) == 4
 
     let # Just to have a local scope for D
-        matrices, ssys = linearize(model, model_inputs, [y]; op)
+        matrices, ssys = linearize(
+            model, model_inputs, [y]; op,
+            guesses = [inertia2.flange_a.phi => 0.0, inertia1.flange_b.phi => 0.0]
+        )
         A, B, C, D = matrices
         obsf = ModelingToolkit.build_explicit_observed_function(
             ssys,

lib/ModelingToolkitBase/test/split_parameters.jl

@@ -164,7 +164,8 @@ if @isdefined(ModelingToolkit)
     inputs = [model.torque.tau.u]
     op = [model.torque.tau.u => 0.0]
     matrices, ssys = ModelingToolkit.linearize(
-        wr(model), inputs, model_outputs; op
+        wr(model), inputs, model_outputs; op,
+        guesses = [model.inertia2.flange_a.phi => 0.0, model.inertia1.flange_b.phi => 0.0]
     )
 
     # Design state-feedback gain using LQR

test/linearize.jl

@@ -122,27 +122,32 @@ connections = [
 
 @named cl = System(connections, t, systems = [f, c, p])
 
-lsys0, ssys = linearize(cl, [f.u], [p.x])
-desired_order = [f.x, p.x]
-lsys = ModelingToolkit.reorder_unknowns(lsys0, unknowns(ssys), desired_order)
-lsys1, ssys = linearize(cl, [f.u], [p.x]; autodiff = AutoFiniteDiff())
-lsys2 = ModelingToolkit.reorder_unknowns(lsys1, unknowns(ssys), desired_order)
+function compare_matrices(reference, value)
+    @assert size(reference.A) == (2, 2) "This testing function is only valid for `2x2` systems"
+    @test isapprox(reference.C, value.C) || isapprox(reverse(reference.C), value.C)
+    colorder = isapprox(reference.C, value.C) ? [1, 2] : [2, 1]
+    @test isapprox(reference.B, value.B) || isapprox(reverse(reference.B), value.B)
+    roworder = isapprox(reference.B, value.B) ? [1, 2] : [2, 1]
+    @test isapprox(reference.D, value.D)
+    @test isapprox(reference.A[roworder, colorder], value.A)
+end
+lsys, ssys = linearize(cl, [f.u], [p.x])
+lsys2, ssys = linearize(cl, [f.u], [p.x]; autodiff = AutoFiniteDiff())
 
-@test lsys.A == lsys2.A == [-2 0; 1 -2]
-@test lsys.B == lsys2.B == reshape([1, 0], 2, 1)
-@test lsys.C == lsys2.C == [0 1]
-@test lsys.D[] == lsys2.D[] == 0
+compare_matrices(lsys, lsys2)
+compare_matrices(lsys, (; A = [-2 0; 1 -2], B = reshape([1, 0], 2, 1), C = [0 1], D = [0;;]))
 
 ## Symbolic linearization
 lsyss_ns, ssys_ns = ModelingToolkit.linearize_symbolic(cl, [f.u], [p.x], split = false)
 lsyss, ssys = ModelingToolkit.linearize_symbolic(cl, [f.u], [p.x])
 @test isequal(lsyss.A, lsyss_ns.A)
 
-lsyss = ModelingToolkit.reorder_unknowns(lsyss, unknowns(ssys), [f.x, p.x])
-@test value.(ModelingToolkit.fixpoint_sub(lsyss.A, ModelingToolkit.initial_conditions(cl))) == lsys.A
-@test value.(ModelingToolkit.fixpoint_sub(lsyss.B, ModelingToolkit.initial_conditions(cl))) == lsys.B
-@test value.(ModelingToolkit.fixpoint_sub(lsyss.C, ModelingToolkit.initial_conditions(cl))) == lsys.C
-@test value.(ModelingToolkit.fixpoint_sub(lsyss.D, ModelingToolkit.initial_conditions(cl))) == lsys.D
+_substituter(M, sys) = value.(ModelingToolkit.fixpoint_sub(M, ModelingToolkit.initial_conditions_and_guesses(sys); fold = Val(true)))
+lsys_m = (;
+    A = _substituter(lsyss.A, cl), B = _substituter(lsyss.B, cl),
+    C = _substituter(lsyss.C, cl), D = _substituter(lsyss.D, cl)
+)
+compare_matrices(lsys, lsys_m)
 ##
 using ModelingToolkitStandardLibrary.Blocks: LimPID
 k = 400
@@ -152,14 +157,12 @@ Nd = 10
 @named pid = LimPID(; k, Ti, Td, Nd)
 
 @unpack reference, measurement, ctr_output = pid
-lsys0,
+lsys,
     ssys = linearize(
     pid, [reference.u, measurement.u], [ctr_output.u];
     op = Dict(reference.u => 0.0, measurement.u => 0.0)
 )
 @unpack int, der = pid
-desired_order = [int.x, der.x]
-lsys = ModelingToolkit.reorder_unknowns(lsys0, unknowns(ssys), desired_order)
 
 @test lsys.A == [0 0; 0 -10]
 @test lsys.B == [2 -2; 10 -10]
@@ -171,31 +174,14 @@ lsyss0,
     pid, [reference.u, measurement.u],
     [ctr_output.u]
 )
-lsyss = ModelingToolkit.reorder_unknowns(lsyss0, unknowns(ssys2), desired_order)
-
-@test value.(
-    ModelingToolkit.fixpoint_sub(
-        lsyss.A, ModelingToolkit.initial_conditions_and_guesses(pid); fold = Val(true)
-    )
-) == lsys.A
-@test value.(
-    ModelingToolkit.fixpoint_sub(
-        lsyss.B, ModelingToolkit.initial_conditions_and_guesses(pid); fold = Val(true)
-    )
-) == lsys.B
-@test value.(
-    ModelingToolkit.fixpoint_sub(
-        lsyss.C, ModelingToolkit.initial_conditions_and_guesses(pid); fold = Val(true)
-    )
-) == lsys.C
-@test value.(
-    ModelingToolkit.fixpoint_sub(
-        lsyss.D, ModelingToolkit.initial_conditions_and_guesses(pid); fold = Val(true)
-    )
-) == lsys.D
+lsyss = (;
+    A = _substituter(lsyss0.A, pid), B = _substituter(lsyss0.B, pid),
+    C = _substituter(lsyss0.C, pid), D = _substituter(lsyss0.D, pid)
+)
+compare_matrices(lsys, lsyss)
 
 # Test with the reverse desired unknown order as well to verify that similarity transform and reoreder_unknowns really works
-lsys = ModelingToolkit.reorder_unknowns(lsys, desired_order, reverse(desired_order))
+lsys = ModelingToolkit.reorder_unknowns(lsys, unknowns(ssys), reverse(unknowns(ssys)))
 
 @test lsys.A == [-10 0; 0 0]
 @test lsys.B == [10 -10; 2 -2]

test/structural_transformation/tearing.jl

@@ -156,7 +156,7 @@ eqs = [
 ]
 @named daesys = System(eqs, t)
 newdaesys = mtkcompile(daesys)
-@test issetequal(equations(newdaesys), [D(x) ~ h * z; 0 ~ y + sin(z) - p * t])
+@test issetequal(equations(newdaesys), [D(x) ~ h * z; 0 ~ x + sin(z) - p * t])
 @test issetequal(
     equations(tearing_substitution(newdaesys)), [D(x) ~ h * z; 0 ~ x + sin(z) - p * t]
 )
@@ -260,7 +260,7 @@ end
     @test length(equations(sys3)) == 1
     @test isequal(only(unknowns(sys3)), sys3.capacitor.v)
 
-    idx = findfirst(isequal(sys3.capacitor.i), observables(sys3))
+    idx = findfirst(isequal(sys3.capacitor.p.i), observables(sys3))
     rhs = observed(sys3)[idx].rhs
     @test operation(rhs) === getindex
     @test operation(arguments(rhs)[1]) === (\)

test/structural_transformation/utils.jl

@@ -194,7 +194,7 @@ end
         @testset "DDEs" begin
             function oscillator(; name, k = 1.0, τ = 0.01)
                 @parameters k = k τ = τ
-                @variables x(..) = 0.1 y(t) = 0.1 jcn(t) = 0.0 delx(t)
+                @variables x(..) = 0.1 y(t) = 0.1 jcn(t) = 0.0 [state_priority = -1] delx(t)
                 eqs = [
                     D(x(t)) ~ y,
                     D(y) ~ -k * x(t - τ) + jcn,
@@ -217,11 +217,11 @@ end
             x1 = operation(unwrap(osc1.x))
             x2 = operation(unwrap(osc2.x))
             @test mapping[sys.osc1.x] == (D(sys.osc1.x) ~ sys.osc1.y)
-            @test mapping[sys.osc1.y] == (D(sys.osc1.y) ~ sys.osc1.jcn - sys.osc1.k * x1(t - sys.osc1.τ))
+            @test mapping[sys.osc1.y] == (D(sys.osc1.y) ~ sys.osc2.delx - sys.osc1.k * x1(t - sys.osc1.τ))
             @test mapping[sys.osc1.delx] == (sys.osc1.delx ~ x1(t - sys.osc1.τ))
             @test mapping[sys.osc1.jcn] == (sys.osc1.jcn ~ sys.osc2.delx)
             @test mapping[sys.osc2.x] == (D(sys.osc2.x) ~ sys.osc2.y)
-            @test mapping[sys.osc2.y] == (D(sys.osc2.y) ~ sys.osc2.jcn - sys.osc2.k * x2(t - sys.osc2.τ))
+            @test mapping[sys.osc2.y] == (D(sys.osc2.y) ~ sys.osc1.delx - sys.osc2.k * x2(t - sys.osc2.τ))
             @test mapping[sys.osc2.delx] == (sys.osc2.delx ~ x2(t - sys.osc2.τ))
             @test mapping[sys.osc2.jcn] == (sys.osc2.jcn ~ sys.osc1.delx)
             @test length(mapping) == 8