Merge pull request #4402 from SciML/as/better-alias-elim

feat: eliminate perfect aliases in favor of higher priority variables

Author: AayushSabharwal

Project.toml

@@ -90,7 +90,7 @@ LinearSolve = "3.66"
 Logging = "1"
 ModelingToolkitBase = "1.27.0"
 ModelingToolkitStandardLibrary = "2.20"
-ModelingToolkitTearing = "1.9.0"
+ModelingToolkitTearing = "1.12.1"
 Moshi = "0.3.6"
 NonlinearSolve = "4.3"
 OffsetArrays = "1"
@@ -113,7 +113,7 @@ Serialization = "1"
 Setfield = "0.7, 0.8, 1"
 SimpleNonlinearSolve = "0.1.0, 1, 2"
 SparseArrays = "1"
-StateSelection = "1.5"
+StateSelection = "1.9.1"
 StaticArrays = "1.9.14"
 StochasticDelayDiffEq = "1.11"
 StochasticDiffEq = "6.82.0"

lib/ModelingToolkitBase/src/systems/callbacks.jl

@@ -140,7 +140,7 @@ function AffectSystem(
         iv = t_nounits
         @warn "No independent variable specified. Defaulting to t_nounits."
     end
-    append!(affect, extra_eqs)
+    affect = [affect; extra_eqs]
 
     discrete_parameters = SymbolicAffect(affect; discrete_parameters).discrete_parameters
 

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