From 9460ea4397230b154d65b36bcb8073bae6ec6f52 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20Miclu=C8=9Ba-C=C3=A2mpeanu?= Date: Tue, 6 Aug 2024 01:55:11 +0300 Subject: [PATCH 1/5] docs: fix namespacing in DC motor tutorial --- docs/src/tutorials/dc_motor_pi.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/src/tutorials/dc_motor_pi.md b/docs/src/tutorials/dc_motor_pi.md index 420eb43d4..5b5f6982e 100644 --- a/docs/src/tutorials/dc_motor_pi.md +++ b/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 From 9f7f0004b2e01b665ec2df2f8368f62b78030cf8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20Miclu=C8=9Ba-C=C3=A2mpeanu?= Date: Fri, 18 Oct 2024 03:31:29 +0300 Subject: [PATCH 2/5] build: bump MTK, Symbolics & DataInterpolations Since DataInterpolations is only a test dep, we seem to hit an edgecase in CompatHelper and this was holding back MTK. --- Project.toml | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Project.toml b/Project.toml index b604de635..d348cc349 100644 --- a/Project.toml +++ b/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" From fb9e3f2f1c1ec27e5373afad0997300a614211dc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20Miclu=C8=9Ba-C=C3=A2mpeanu?= Date: Fri, 18 Oct 2024 03:34:18 +0300 Subject: [PATCH 3/5] test: use `SciMLBase.successful_retcode(sol)` --- test/Blocks/nonlinear.jl | 10 +++++----- test/Electrical/analog.jl | 20 ++++++++++---------- test/Magnetic/magnetic.jl | 2 +- test/Thermal/demo.jl | 2 +- test/Thermal/motor.jl | 2 +- test/Thermal/piston.jl | 2 +- test/Thermal/thermal.jl | 10 +++++----- 7 files changed, 24 insertions(+), 24 deletions(-) diff --git a/test/Blocks/nonlinear.jl b/test/Blocks/nonlinear.jl index 3b5f58787..8fa2c2ade 100644 --- a/test/Blocks/nonlinear.jl +++ b/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 diff --git a/test/Electrical/analog.jl b/test/Electrical/analog.jl index fbbeb24a9..1ac86ef73 100644 --- a/test/Electrical/analog.jl +++ b/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 diff --git a/test/Magnetic/magnetic.jl b/test/Magnetic/magnetic.jl index c81439fc1..0564aafa9 100644 --- a/test/Magnetic/magnetic.jl +++ b/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 diff --git a/test/Thermal/demo.jl b/test/Thermal/demo.jl index 5403681fc..14b61b215 100644 --- a/test/Thermal/demo.jl +++ b/test/Thermal/demo.jl @@ -22,5 +22,5 @@ using OrdinaryDiffEq: ReturnCode.Success sys = structural_simplify(model) prob = ODEProblem(sys, [mass1.der_T => 1.0, mass2.der_T => 1.0], (0, 3.0)) sol = solve(prob, Tsit5()) - @test sol.retcode == Success + @test SciMLBase.successful_retcode(sol) end diff --git a/test/Thermal/motor.jl b/test/Thermal/motor.jl index c02d43538..a8f951256 100644 --- a/test/Thermal/motor.jl +++ b/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] ≈ diff --git a/test/Thermal/piston.jl b/test/Thermal/piston.jl index d514c5f4f..f1c8e3cb3 100644 --- a/test/Thermal/piston.jl +++ b/test/Thermal/piston.jl @@ -40,7 +40,7 @@ using ModelingToolkitStandardLibrary.Blocks # 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 diff --git a/test/Thermal/thermal.jl b/test/Thermal/thermal.jl index 85348faf5..a62ca5e61 100644 --- a/test/Thermal/thermal.jl +++ b/test/Thermal/thermal.jl @@ -29,7 +29,7 @@ using OrdinaryDiffEq: ReturnCode.Success # Check if Relative temperature sensor reads the temperature of heat capacitor # when connected to a thermal conductor and a fixed temperature source - @test sol.retcode == Success + @test SciMLBase.successful_retcode(sol) @test sol[reltem_sensor.T] + sol[tem_src.port.T] == sol[mass1.T] + sol[th_conductor.dT] @info "Building a two-body system..." @@ -49,7 +49,7 @@ using OrdinaryDiffEq: ReturnCode.Success prob = ODEProblem(sys, u0, (0, 3.0)) sol = solve(prob, Tsit5()) - @test sol.retcode == Success + @test SciMLBase.successful_retcode(sol) m1, m2 = sol.u[end] @test m1≈m2 atol=1e-1 mass_T = reduce(hcat, sol.u) @@ -85,7 +85,7 @@ end prob = ODEProblem(sys, u0, (0, 3.0)) sol = solve(prob, Tsit5()) - @test sol.retcode == Success + @test SciMLBase.successful_retcode(sol) @test sol[th_conductor.dT] .* G == sol[th_conductor.Q_flow] @test sol[th_conductor.Q_flow] ≈ sol[hf_sensor1.Q_flow] + sol[flow_src.port.Q_flow] @@ -127,7 +127,7 @@ end prob = ODEProblem(sys, u0, (0, 3.0)) sol = solve(prob, Rodas4()) - @test sol.retcode == Success + @test SciMLBase.successful_retcode(sol) @test sol[dissipator.dT] == sol[radiator.port_a.T] - sol[radiator.port_b.T] rad_Q_flow = G * σ * (T_gas^4 - T_coolant^4) @test sol[radiator.Q_flow] == fill(rad_Q_flow, length(sol[radiator.Q_flow])) @@ -158,7 +158,7 @@ end prob = ODEProblem(sys, u0, (0, 3.0)) sol = solve(prob, Rodas4()) - @test sol.retcode == Success + @test SciMLBase.successful_retcode(sol) @test sol[collector.port_b.Q_flow] + sol[collector.port_a1.Q_flow] + sol[collector.port_a2.Q_flow] == zeros(length(sol[collector.port_b.Q_flow])) From 574b0fd1001eb1a44474d901b02f4833862a1ba0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20Miclu=C8=9Ba-C=C3=A2mpeanu?= Date: Fri, 18 Oct 2024 03:36:25 +0300 Subject: [PATCH 4/5] test: remove unneeded initial conditions with newer versions of MTK we no longer need to provides these, as initialization can find the appropriate initial conditions. If the extra u0s are provided, they lead to initialization failures. --- test/Blocks/sources.jl | 2 +- test/Thermal/demo.jl | 2 +- test/Thermal/piston.jl | 4 +--- test/Thermal/thermal.jl | 17 ++++------------- 4 files changed, 7 insertions(+), 18 deletions(-) diff --git a/test/Blocks/sources.jl b/test/Blocks/sources.jl index f03a6de4a..c073385b1 100644 --- a/test/Blocks/sources.jl +++ b/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( diff --git a/test/Thermal/demo.jl b/test/Thermal/demo.jl index 14b61b215..001a43840 100644 --- a/test/Thermal/demo.jl +++ b/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 SciMLBase.successful_retcode(sol) end diff --git a/test/Thermal/piston.jl b/test/Thermal/piston.jl index f1c8e3cb3..1d3426902 100644 --- a/test/Thermal/piston.jl +++ b/test/Thermal/piston.jl @@ -33,9 +33,7 @@ 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 diff --git a/test/Thermal/thermal.jl b/test/Thermal/thermal.jl index a62ca5e61..49c5d634c 100644 --- a/test/Thermal/thermal.jl +++ b/test/Thermal/thermal.jl @@ -22,8 +22,7 @@ using OrdinaryDiffEq: ReturnCode.Success @named h1 = ODESystem(eqs, t, systems = [mass1, reltem_sensor, tem_src, th_conductor]) sys = structural_simplify(h1) - u0 = [mass1.T => 2.0 - mass1.der_T => 1.0] + u0 = [mass1.T => 2] prob = ODEProblem(sys, u0, (0, 2.0)) sol = solve(prob, Tsit5()) @@ -42,10 +41,7 @@ using OrdinaryDiffEq: ReturnCode.Success sys = structural_simplify(h2) u0 = [mass1.T => 1.0 - mass2.T => 10.0 - final_T => 12 - mass1.der_T => 1.0 - mass2.der_T => 1.0] + mass2.T => 10.0] prob = ODEProblem(sys, u0, (0, 3.0)) sol = solve(prob, Tsit5()) @@ -79,9 +75,7 @@ end th_resistor, flow_src, th_ground, th_conductor]) sys = structural_simplify(h2) - u0 = [mass1.T => 10.0 - th_resistor.Q_flow => 1.0 - mass1.der_T => 1.0] + u0 = [mass1.T => 10.0] prob = ODEProblem(sys, u0, (0, 3.0)) sol = solve(prob, Tsit5()) @@ -121,9 +115,7 @@ end ]) sys = structural_simplify(rad) - u0 = [base.Q_flow => 10 - dissipator.Q_flow => 10 - mass.T => T_gas] + u0 = [mass.T => T_gas] prob = ODEProblem(sys, u0, (0, 3.0)) sol = solve(prob, Rodas4()) @@ -152,7 +144,6 @@ end sys = structural_simplify(coll) u0 = [ - th_resistor.Q_flow => 1.0, mass.T => 0.0 ] prob = ODEProblem(sys, u0, (0, 3.0)) From 555145ef5bd2cefe68d4a0a44de69d59545e266b Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sat, 19 Oct 2024 09:42:50 -0400 Subject: [PATCH 5/5] Update translational.jl --- test/Mechanical/translational.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/Mechanical/translational.jl b/test/Mechanical/translational.jl index ec077de13..c3c6dc2e5 100644 --- a/test/Mechanical/translational.jl +++ b/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