Skip to content

Commit 04ca326

Browse files
AayushSabharwalAayushSabharwal
committed
test: add tests from MTKStdlib
1 parent f0bab4a commit 04ca326

File tree

1 file changed

+239
-3
lines changed

1 file changed

+239
-3
lines changed

‎test/analysis_points.jl

Lines changed: 239 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1,11 +1,14 @@
1-
using ModelingToolkit, ModelingToolkitStandardLibrary
1+
using ModelingToolkit, ModelingToolkitStandardLibrary.Blocks
2+
using OrdinaryDiffEq, LinearAlgebra
23
using Test
3-
using ModelingToolkit: t_nounits as t, D_nounits as D
4+
using ModelingToolkit: t_nounits as t, D_nounits as D, AnalysisPoint, get_sensitivity,
5+
get_comp_sensitivity, get_looptransfer
6+
using Symbolics: NAMESPACE_SEPARATOR
47

58
@testset "AnalysisPoint is lowered to `connect`" begin
69
@named P = FirstOrder(k = 1, T = 1)
710
@named C = Gain(; k = -1)
8-
11+
912
ap = AnalysisPoint(:plant_input)
1013
eqs = [connect(P.output, C.input)
1114
connect(C.output, ap, P.input)]
@@ -21,3 +24,236 @@ using ModelingToolkit: t_nounits as t, D_nounits as D
2124

2225
@test isequal(sys_ap2, sys_normal2)
2326
end
27+
28+
# also tests `connect(input, name::Symbol, outputs...)` syntax
29+
@testset "AnalysisPoint is accessible via `getproperty`" begin
30+
@named P = FirstOrder(k = 1, T = 1)
31+
@named C = Gain(; k = -1)
32+
33+
eqs = [connect(P.output, C.input), connect(C.output, :plant_input, P.input)]
34+
sys_ap = ODESystem(eqs, t, systems = [P, C], name = :hej)
35+
ap2 = @test_nowarn sys_ap.plant_input
36+
@test nameof(ap2) == Symbol(join(["hej", "plant_input"], NAMESPACE_SEPARATOR))
37+
@named sys = ODESystem(Equation[], t; systems = [sys_ap])
38+
ap3 = @test_nowarn sys.hej.plant_input
39+
@test nameof(ap3) == Symbol(join(["sys", "hej", "plant_input"], NAMESPACE_SEPARATOR))
40+
sys = complete(sys)
41+
ap4 = sys.hej.plant_input
42+
@test nameof(ap4) == Symbol(join(["hej", "plant_input"], NAMESPACE_SEPARATOR))
43+
end
44+
45+
### Ported from MTKStdlib
46+
47+
@named P = FirstOrder(k = 1, T = 1)
48+
@named C = Gain(; k = -1)
49+
50+
ap = AnalysisPoint(:plant_input)
51+
eqs = [connect(P.output, C.input), connect(C.output, ap, P.input)]
52+
sys = ODESystem(eqs, t, systems = [P, C], name = :hej)
53+
@named nested_sys = ODESystem(Equation[], t; systems = [sys])
54+
55+
@testset "simplifies and solves" begin
56+
ssys = structural_simplify(sys)
57+
prob = ODEProblem(ssys, [P.x => 1], (0, 10))
58+
sol = solve(prob, Rodas5())
59+
@test norm(sol.u[1]) >= 1
60+
@test norm(sol.u[end]) < 1e-6 # This fails without the feedback through C
61+
end
62+
63+
@testset "get_sensitivity - $name" for (name, sys, ap) in [
64+
("inner", sys, sys.plant_input),
65+
("nested", nested_sys, nested_sys.hej.plant_input),
66+
("inner - nonamespace", sys, :plant_input),
67+
("inner - Symbol", sys, nameof(sys.plant_input)),
68+
("nested - Symbol", nested_sys, nameof(nested_sys.hej.plant_input))
69+
]
70+
matrices, _ = get_sensitivity(sys, ap)
71+
@test matrices.A[] == -2
72+
@test matrices.B[] * matrices.C[] == -1 # either one negative
73+
@test matrices.D[] == 1
74+
end
75+
76+
@testset "get_comp_sensitivity - $name" for (name, sys, ap) in [
77+
("inner", sys, sys.plant_input),
78+
("nested", nested_sys, nested_sys.hej.plant_input),
79+
("inner - nonamespace", sys, :plant_input),
80+
("inner - Symbol", sys, nameof(sys.plant_input)),
81+
("nested - Symbol", nested_sys, nameof(nested_sys.hej.plant_input))
82+
]
83+
matrices, _ = get_comp_sensitivity(sys, ap)
84+
@test matrices.A[] == -2
85+
@test matrices.B[] * matrices.C[] == 1 # both positive or negative
86+
@test matrices.D[] == 0
87+
end
88+
89+
#=
90+
# Equivalent code using ControlSystems. This can be used to verify the expected results tested for above.
91+
using ControlSystemsBase
92+
P = tf(1.0, [1, 1])
93+
C = 1 # Negative feedback assumed in ControlSystems
94+
S = sensitivity(P, C) # or feedback(1, P*C)
95+
T = comp_sensitivity(P, C) # or feedback(P*C)
96+
=#
97+
98+
@testset "get_looptransfer - $name" for (name, sys, ap) in [
99+
("inner", sys, sys.plant_input),
100+
("nested", nested_sys, nested_sys.hej.plant_input),
101+
("inner - nonamespace", sys, :plant_input),
102+
("inner - Symbol", sys, nameof(sys.plant_input)),
103+
("nested - Symbol", nested_sys, nameof(nested_sys.hej.plant_input))
104+
]
105+
matrices, _ = get_looptransfer(sys, ap)
106+
@test matrices.A[] == -1
107+
@test matrices.B[] * matrices.C[] == -1 # either one negative
108+
@test matrices.D[] == 0
109+
end
110+
111+
#=
112+
# Equivalent code using ControlSystems. This can be used to verify the expected results tested for above.
113+
using ControlSystemsBase
114+
P = tf(1.0, [1, 1])
115+
C = -1
116+
L = P*C
117+
=#
118+
119+
@testset "open_loop - $name" for (name, sys, ap) in [
120+
("inner", sys, sys.plant_input),
121+
("nested", nested_sys, nested_sys.hej.plant_input),
122+
("inner - nonamespace", sys, :plant_input),
123+
("inner - Symbol", sys, nameof(sys.plant_input)),
124+
("nested - Symbol", nested_sys, nameof(nested_sys.hej.plant_input))
125+
]
126+
open_sys, (du, u) = open_loop(sys, ap)
127+
matrices, _ = linearize(open_sys, [du], [u])
128+
@test matrices.A[] == -1
129+
@test matrices.B[] * matrices.C[] == -1 # either one negative
130+
@test matrices.D[] == 0
131+
end
132+
133+
# Multiple analysis points
134+
135+
eqs = [connect(P.output, :plant_output, C.input)
136+
connect(C.output, :plant_input, P.input)]
137+
sys = ODESystem(eqs, t, systems = [P, C], name = :hej)
138+
@named nested_sys = ODESystem(Equation[], t; systems = [sys])
139+
140+
@testset "get_sensitivity - $name" for (name, sys, ap) in [
141+
("inner", sys, sys.plant_input),
142+
("nested", nested_sys, nested_sys.hej.plant_input),
143+
("inner - nonamespace", sys, :plant_input),
144+
("inner - Symbol", sys, nameof(sys.plant_input)),
145+
("nested - Symbol", nested_sys, nameof(nested_sys.hej.plant_input))
146+
]
147+
matrices, _ = get_sensitivity(sys, ap)
148+
@test matrices.A[] == -2
149+
@test matrices.B[] * matrices.C[] == -1 # either one negative
150+
@test matrices.D[] == 1
151+
end
152+
153+
@testset "linearize - $name" for (name, sys, inputap, outputap) in [
154+
("inner", sys, sys.plant_input, sys.plant_output),
155+
("nested", nested_sys, nested_sys.hej.plant_input, nested_sys.hej.plant_output)
156+
]
157+
@testset "input - $(typeof(input)), output - $(typeof(output))" for (input, output) in [
158+
(inputap, outputap),
159+
(nameof(inputap), outputap),
160+
(inputap, nameof(outputap)),
161+
(nameof(inputap), nameof(outputap)),
162+
(inputap, [outputap]),
163+
(nameof(inputap), [outputap]),
164+
(inputap, [nameof(outputap)]),
165+
(nameof(inputap), [nameof(outputap)])
166+
]
167+
if input isa Symbol
168+
# broken because MTKStdlib defines the method for
169+
# `input::Union{Symbol, Vector{Symbol}}` which we can't directly call
170+
@test_broken linearize(sys, input, output)
171+
linfun, ssys = @invoke linearization_function(sys::AbstractSystem,
172+
input::Union{Symbol, Vector{Symbol}, AnalysisPoint, Vector{AnalysisPoint}},
173+
output::Any)
174+
matrices = linearize(ssys, linfun)
175+
else
176+
matrices, _ = linearize(sys, input, output)
177+
end
178+
# Result should be the same as feedpack(P, 1), i.e., the closed-loop transfer function from plant input to plant output
179+
@test matrices.A[] == -2
180+
@test matrices.B[] * matrices.C[] == 1 # both positive
181+
@test matrices.D[] == 0
182+
end
183+
end
184+
185+
@testset "linearize - variable output - $name" for (name, sys, input, output) in [
186+
("inner", sys, sys.plant_input, P.output.u),
187+
("nested", nested_sys, nested_sys.hej.plant_input, sys.P.output.u)
188+
]
189+
matrices, _ = linearize(sys, input, [output])
190+
@test matrices.A[] == -2
191+
@test matrices.B[] * matrices.C[] == 1 # both positive
192+
@test matrices.D[] == 0
193+
end
194+
195+
using ModelingToolkit, OrdinaryDiffEq, LinearAlgebra
196+
using ModelingToolkitStandardLibrary.Mechanical.Rotational
197+
using ModelingToolkitStandardLibrary.Blocks: Sine, PID, SecondOrder, Step, RealOutput
198+
using ModelingToolkit: connect
199+
200+
@testset "Complicated model" begin
201+
# Parameters
202+
m1 = 1
203+
m2 = 1
204+
k = 1000 # Spring stiffness
205+
c = 10 # Damping coefficient
206+
@named inertia1 = Inertia(; J = m1)
207+
@named inertia2 = Inertia(; J = m2)
208+
@named spring = Spring(; c = k)
209+
@named damper = Damper(; d = c)
210+
@named torque = Torque()
211+
212+
function SystemModel(u = nothing; name = :model)
213+
eqs = [connect(torque.flange, inertia1.flange_a)
214+
connect(inertia1.flange_b, spring.flange_a, damper.flange_a)
215+
connect(inertia2.flange_a, spring.flange_b, damper.flange_b)]
216+
if u !== nothing
217+
push!(eqs, connect(torque.tau, u.output))
218+
return ODESystem(eqs, t;
219+
systems = [
220+
torque,
221+
inertia1,
222+
inertia2,
223+
spring,
224+
damper,
225+
u
226+
],
227+
name)
228+
end
229+
ODESystem(eqs, t; systems = [torque, inertia1, inertia2, spring, damper], name)
230+
end
231+
232+
@named r = Step(start_time = 0)
233+
model = SystemModel()
234+
@named pid = PID(k = 100, Ti = 0.5, Td = 1)
235+
@named filt = SecondOrder(d = 0.9, w = 10)
236+
@named sensor = AngleSensor()
237+
@named er = Add(k2 = -1)
238+
239+
connections = [connect(r.output, :r, filt.input)
240+
connect(filt.output, er.input1)
241+
connect(pid.ctr_output, :u, model.torque.tau)
242+
connect(model.inertia2.flange_b, sensor.flange)
243+
connect(sensor.phi, :y, er.input2)
244+
connect(er.output, :e, pid.err_input)]
245+
246+
closed_loop = ODESystem(connections, t, systems = [model, pid, filt, sensor, r, er],
247+
name = :closed_loop, defaults = [
248+
model.inertia1.phi => 0.0,
249+
model.inertia2.phi => 0.0,
250+
model.inertia1.w => 0.0,
251+
model.inertia2.w => 0.0,
252+
filt.x => 0.0,
253+
filt.xd => 0.0
254+
])
255+
256+
sys = structural_simplify(closed_loop)
257+
prob = ODEProblem(sys, unknowns(sys) .=> 0.0, (0.0, 4.0))
258+
sol = solve(prob, Rodas5P(), reltol = 1e-6, abstol = 1e-9)
259+
end

0 commit comments

Comments
 (0)