Merge pull request #338 from SebastianM-C/smc/update

Update MTK

Author: ChrisRackauckas

Project.toml

@@ -15,15 +15,15 @@ Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 Aqua = "0.8"
 ChainRulesCore = "1.24"
 ControlSystemsBase = "1.4"
-DataInterpolations = "4.6"
+DataInterpolations = "6"
 DiffEqBase = "6.152"
 IfElse = "0.1"
 LinearAlgebra = "1.10"
-ModelingToolkit = "9.32"
+ModelingToolkit = "9.46.1"
 OrdinaryDiffEq = "6.87"
 OrdinaryDiffEqDefault = "1.1"
 SafeTestsets = "0.1"
-Symbolics = "5.36, 6"
+Symbolics = "6.14"
 Test = "1"
 julia = "1.10"
 

docs/src/tutorials/dc_motor_pi.md

@@ -92,8 +92,8 @@ plot(p1, p2, layout = (2, 1))
 When implementing and tuning a control system in simulation, it is a good practice to analyze the closed-loop properties and verify robustness of the closed-loop with respect to, e.g., modeling errors. To facilitate this, we added two analysis points to the set of connections above, more specifically, we added the analysis points named `:y` and `:u` to the connections (for more details on analysis points, see [Linear Analysis](@ref))
 
 ```julia
-connect(sys.speed_sensor.w, :y, feedback.input2)
-connect(sys.pi_controller.ctr_output, :u, source.V)
+connect(sys.speed_sensor.w, :y, sys.feedback.input2)
+connect(sys.pi_controller.ctr_output, :u, sys.source.V)
 ```
 
 one at the plant output (`:y`) and one at the plant input (`:u`). We may use these analysis points to calculate, e.g., sensitivity functions, illustrated below. Here, we calculate the sensitivity function $S(s)$ and the complimentary sensitivity function $T(s) = I - S(s)$, defined as

test/Blocks/nonlinear.jl

@@ -20,7 +20,7 @@ using OrdinaryDiffEq: ReturnCode.Success
         prob = ODEProblem(sys, [int.x => 1.0], (0.0, 1.0))
 
         sol = solve(prob, Rodas4())
-        @test sol.retcode == Success
+        @test SciMLBase.successful_retcode(sol)
         @test sol[int.output.u][end] ≈ 2
         @test sol[sat.output.u][end] ≈ 0.8
     end
@@ -42,7 +42,7 @@ using OrdinaryDiffEq: ReturnCode.Success
         prob = ODEProblem(sys, unknowns(sys) .=> 0.0, (0.0, 10.0))
 
         sol = solve(prob, Rodas4())
-        @test sol.retcode == Success
+        @test SciMLBase.successful_retcode(sol)
         @test all(abs.(sol[lim.output.u]) .<= 0.5)
         @test all(isapprox.(sol[lim.output.u], _clamp.(sol[source.output.u], y_min, y_max),
             atol = 1e-2))
@@ -69,7 +69,7 @@ end
         prob = ODEProblem(sys, [int.x => 1.0], (0.0, 1.0))
         sol = solve(prob, Rodas4())
 
-        @test sol.retcode == Success
+        @test SciMLBase.successful_retcode(sol)
         @test all(sol[int.output.u][end] .≈ 2)
     end
 
@@ -89,7 +89,7 @@ end
         prob = ODEProblem(sys, [int.x => 1.0], (0.0, 10.0))
         sol = solve(prob, Rodas4())
 
-        @test sol.retcode == Success
+        @test SciMLBase.successful_retcode(sol)
         @test all(sol[dz.output.u] .<= 2)
         @test all(sol[dz.output.u] .>= -1)
         @test all(isapprox.(sol[dz.output.u],
@@ -115,7 +115,7 @@ end
 
     tS = 0.01
     sol = solve(prob, Rodas4(), saveat = tS, abstol = 1e-10, reltol = 1e-10)
-    @test sol.retcode == Success
+    @test SciMLBase.successful_retcode(sol)
     @test all(abs.(sol[rl.output.u]) .<= 0.51)
     @test all(-1 - 1e-5 .<= diff(sol[rl.output.u]) ./ tS .<= 1 + 1e-5) # just an approximation
 end

test/Blocks/sources.jl

@@ -25,7 +25,7 @@ end
 
 @testset "TimeVaryingFunction" begin
     f(t) = t^2 + 1
-    vars = @variables y(t)=1 dy(t)=0 ddy(t)=0
+    vars = @variables y(t) dy(t) ddy(t)
     @named src = TimeVaryingFunction(f)
     @named int = Integrator()
     @named iosys = ODESystem(

test/Electrical/analog.jl

@@ -48,7 +48,7 @@ using OrdinaryDiffEq: ReturnCode.Success
     # Plots.plot(sol; vars=[capacitor.v, voltage_sensor.v])
     # Plots.plot(sol; vars=[power_sensor.power, capacitor.i * capacitor.v])
     # Plots.plot(sol; vars=[resistor.i, current_sensor.i])
-    @test sol.retcode == Success
+    @test SciMLBase.successful_retcode(sol)
     @test sol[capacitor.v]≈sol[voltage_sensor.v] atol=1e-3
     @test sol[power_sensor.power]≈sol[capacitor.i * capacitor.v] atol=1e-3
     @test sol[resistor.i]≈sol[current_sensor.i] atol=1e-3
@@ -75,7 +75,7 @@ end
     sys = structural_simplify(model)
     prob = ODEProblem(sys, Pair[R2.i => 0.0], (0, 2.0))
     sol = solve(prob, Rodas4()) # has no state; does not work with Tsit5
-    @test sol.retcode == Success
+    @test SciMLBase.successful_retcode(sol)
     @test sol[short.v] == sol[R0.v] == zeros(length(sol.t))
     @test sol[R0.i] == zeros(length(sol.t))
     @test sol[R1.p.v][end]≈10 atol=1e-3
@@ -103,7 +103,7 @@ end
     sol = solve(prob, Tsit5())
 
     # Plots.plot(sol; vars=[source.v, capacitor.v])
-    @test sol.retcode == Success
+    @test SciMLBase.successful_retcode(sol)
     @test sol[capacitor.v][end]≈10 atol=1e-3
 end
 
@@ -127,7 +127,7 @@ end
     sol = solve(prob, Tsit5())
 
     # Plots.plot(sol; vars=[inductor.i, inductor.i])
-    @test sol.retcode == Success
+    @test SciMLBase.successful_retcode(sol)
     @test sol[inductor.i][end]≈10 atol=1e-3
 end
 
@@ -160,9 +160,9 @@ end
         sys = structural_simplify(model)
         prob = ODEProblem(sys, Pair[], (0.0, 10.0))
         sol = solve(prob, Tsit5())
-        @test sol.retcode == Success
+        @test SciMLBase.successful_retcode(sol)
         sol = solve(prob, Rodas4())
-        @test sol.retcode == Success
+        @test SciMLBase.successful_retcode(sol)
 
         # Plots.plot(sol; vars=[voltage.v, capacitor.v])
     end
@@ -188,7 +188,7 @@ end
     prob = ODEProblem(sys, Pair[], (0.0, 10.0))
     sol = solve(prob, Tsit5())
     y(x, st) = (x .> st) .* abs.(collect(x) .- st)
-    @test sol.retcode == Success
+    @test SciMLBase.successful_retcode(sol)
     @test sum(reduce(vcat, sol[capacitor.v]) .- y(sol.t, start_time))≈0 atol=1e-2
 end
 
@@ -228,7 +228,7 @@ end
           R1.v => 0.0]
     prob = ODEProblem(sys, u0, (0, 100.0))
     sol = solve(prob, Rodas4())
-    @test sol.retcode == Success
+    @test SciMLBase.successful_retcode(sol)
     @test sol[opamp.v2] == sol[C1.v] # Not a great one however. Rely on the plot
     @test sol[opamp.p2.v] == sol[sensor.v]
 
@@ -298,7 +298,7 @@ _damped_sine_wave(x, f, A, st, ϕ, d) = exp((st - x) * d) * A * sin(2 * π * f *
         prob = ODEProblem(vsys, u0, (0, 10.0))
         sol = solve(prob, dt = 0.1, Tsit5())
 
-        @test sol.retcode == Success
+        @test SciMLBase.successful_retcode(sol)
         @test sol[voltage.V.u]≈waveforms(i, sol.t) atol=1e-1
         @test sol[voltage.p.v] ≈ sol[voltage.V.u]
         # For visual inspection
@@ -364,7 +364,7 @@ end
         prob = ODEProblem(isys, u0, (0, 10.0))
         sol = solve(prob, dt = 0.1, Tsit5())
 
-        @test sol.retcode == Success
+        @test SciMLBase.successful_retcode(sol)
         @test sol[current.I.u]≈waveforms(i, sol.t) atol=1e-1
         @test sol[current.I.u]≈sol[current.p.i] atol=1e-1
         # For visual inspection

test/Magnetic/magnetic.jl

@@ -54,7 +54,7 @@ using OrdinaryDiffEq: ReturnCode.Success
     # Plots.plot(sol; vars=[r.i])
     # Plots.plot(sol; vars=[r_mFe.V_m, r_mFe.Phi])
 
-    @test sol.retcode == Success
+    @test SciMLBase.successful_retcode(sol)
     @test sol[r_mFe.Phi] == sol[r_mAirPar.Phi]
     @test all(sol[coil.port_p.Phi] + sol[r_mLeak.Phi] + sol[r_mAirPar.Phi] .== 0)
 end

test/Mechanical/translational.jl

@@ -155,7 +155,7 @@ end
             delta_s = 1 / 1000
             s_b = 2 - delta_s + 1
 
-            @test sol[s.pos_value.u][end]≈s_b atol=1e-3
+            @test sol[s.pos_value.u][end]≈1.0 atol=1e-3
             @test all(sol[s.spring.flange_a.f] .== sol[s.force_output.u])
         end
 

test/Thermal/demo.jl

@@ -20,7 +20,7 @@ using OrdinaryDiffEq: ReturnCode.Success
     @named model = ODESystem(connections, t,
         systems = [mass1, mass2, conduction, Tsensor1, Tsensor2])
     sys = structural_simplify(model)
-    prob = ODEProblem(sys, [mass1.der_T => 1.0, mass2.der_T => 1.0], (0, 3.0))
+    prob = ODEProblem(sys, [], (0, 3.0))
     sol = solve(prob, Tsit5())
-    @test sol.retcode == Success
+    @test SciMLBase.successful_retcode(sol)
 end

test/Thermal/motor.jl

@@ -46,7 +46,7 @@ using ModelingToolkitStandardLibrary.Blocks
     sol = solve(prob)
 
     # plot(sol; vars=[T_winding.T, T_core.T])
-    @test sol.retcode == Success
+    @test SciMLBase.successful_retcode(sol)
     @test sol[motor.T_winding.T] == sol[motor.winding.T]
     @test sol[motor.T_core.T] == sol[motor.core.T]
     @test sol[-motor.core.port.Q_flow] ≈

test/Thermal/piston.jl

@@ -33,14 +33,12 @@ using ModelingToolkitStandardLibrary.Blocks
 
     @mtkbuild piston = Piston()
 
-    u0 = [piston.coolant.dT => 5.0
-          piston.wall.Q_flow => 10.0]
-    prob = ODEProblem(piston, u0, (0, 3.0))
+    prob = ODEProblem(piston, [], (0, 3.0))
     sol = solve(prob)
 
     # Heat-flow-rate is equal in magnitude
     # and opposite in direction
-    @test sol.retcode == Success
+    @test SciMLBase.successful_retcode(sol)
     # The initial value doesn't add up to absolute zero, while the rest do. To avoid
     # tolerance on the latter, the test is split in two parts.
     @test sol[piston.gas.Q_flow][1] + sol[piston.coolant.Q_flow][1]≈0 atol=1e-6