Skip to content

Commit

Permalink
Setup tests for MTKstdlib v2
Browse files Browse the repository at this point in the history
We don't have a dep on it (it would be circular anyways), but this is all that's needed for mtkstdlib v2. I just did `@test_broken` on the two using Torque because they need to handle the booleans in the `@mtkmodel`, and we can skip that for now.
  • Loading branch information
ChrisRackauckas committed Jul 19, 2023
1 parent 24d0d7c commit 3bbad02
Show file tree
Hide file tree
Showing 2 changed files with 80 additions and 75 deletions.
153 changes: 79 additions & 74 deletions test/input_output_handling.jl
Original file line number Diff line number Diff line change
Expand Up @@ -122,36 +122,39 @@ syss = structural_simplify(sys2)
using ModelingToolkitStandardLibrary
using ModelingToolkitStandardLibrary.Mechanical.Rotational
t = ModelingToolkitStandardLibrary.Mechanical.Rotational.t
@named inertia1 = Inertia(; J = 1)
@named inertia2 = Inertia(; J = 1)
@named spring = Spring(; c = 10)
@named damper = Damper(; d = 3)
@named torque = Torque()
@variables y(t) = 0
eqs = [connect(torque.flange, inertia1.flange_a)
connect(inertia1.flange_b, spring.flange_a, damper.flange_a)
connect(inertia2.flange_a, spring.flange_b, damper.flange_b)
y ~ inertia2.w + torque.tau.u]
model = ODESystem(eqs, t; systems = [torque, inertia1, inertia2, spring, damper],
name = :name)
model_outputs = [inertia1.w, inertia2.w, inertia1.phi, inertia2.phi]
model_inputs = [torque.tau.u]
matrices, ssys = linearize(model, model_inputs, model_outputs)
@test length(ModelingToolkit.outputs(ssys)) == 4

if VERSION >= v"1.8" # :opaque_closure not supported before
matrices, ssys = linearize(model, model_inputs, [y])
A, B, C, D = matrices
obsf = ModelingToolkit.build_explicit_observed_function(ssys,
[y],
inputs = [torque.tau.u],
drop_expr = identity)
x = randn(size(A, 1))
u = randn(size(B, 2))
p = getindex.(Ref(ModelingToolkit.defaults(ssys)), parameters(ssys))
y1 = obsf(x, u, p, 0)
y2 = C * x + D * u
@test y1[] y2[]

@test_broken begin
@named inertia1 = Inertia(; J = 1)
@named inertia2 = Inertia(; J = 1)
@named spring = Spring(; c = 10)
@named damper = Damper(; d = 3)
@named torque = Torque()
@variables y(t) = 0
eqs = [connect(torque.flange, inertia1.flange_a)
connect(inertia1.flange_b, spring.flange_a, damper.flange_a)
connect(inertia2.flange_a, spring.flange_b, damper.flange_b)
y ~ inertia2.w + torque.tau.u]
model = ODESystem(eqs, t; systems = [torque, inertia1, inertia2, spring, damper],
name = :name)
model_outputs = [inertia1.w, inertia2.w, inertia1.phi, inertia2.phi]
model_inputs = [torque.tau.u]
matrices, ssys = linearize(model, model_inputs, model_outputs)
@test length(ModelingToolkit.outputs(ssys)) == 4

if VERSION >= v"1.8" # :opaque_closure not supported before
matrices, ssys = linearize(model, model_inputs, [y])
A, B, C, D = matrices
obsf = ModelingToolkit.build_explicit_observed_function(ssys,
[y],
inputs = [torque.tau.u],
drop_expr = identity)
x = randn(size(A, 1))
u = randn(size(B, 2))
p = getindex.(Ref(ModelingToolkit.defaults(ssys)), parameters(ssys))
y1 = obsf(x, u, p, 0)
y2 = C * x + D * u
@test y1[] y2[]
end
end

## Code generation with unbound inputs
Expand Down Expand Up @@ -260,58 +263,60 @@ m2 = 1
k = 1000 # Spring stiffness
c = 10 # Damping coefficient

@named inertia1 = Rotational.Inertia(; J = m1)
@named inertia2 = Rotational.Inertia(; J = m2)
@named spring = Rotational.Spring(; c = k)
@named damper = Rotational.Damper(; d = c)
@named torque = Rotational.Torque()

function SystemModel(u = nothing; name = :model)
eqs = [connect(torque.flange, inertia1.flange_a)
connect(inertia1.flange_b, spring.flange_a, damper.flange_a)
connect(inertia2.flange_a, spring.flange_b, damper.flange_b)]
if u !== nothing
push!(eqs, connect(torque.tau, u.output))
return @named model = ODESystem(eqs, t;
systems = [
torque,
inertia1,
inertia2,
spring,
damper,
u,
])
@test_broken begin
@named inertia1 = Rotational.Inertia(; J = m1)
@named inertia2 = Rotational.Inertia(; J = m2)
@named spring = Rotational.Spring(; c = k)
@named damper = Rotational.Damper(; d = c)
@named torque = Rotational.Torque()

function SystemModel(u = nothing; name = :model)
eqs = [connect(torque.flange, inertia1.flange_a)
connect(inertia1.flange_b, spring.flange_a, damper.flange_a)
connect(inertia2.flange_a, spring.flange_b, damper.flange_b)]
if u !== nothing
push!(eqs, connect(torque.tau, u.output))
return @named model = ODESystem(eqs, t;
systems = [
torque,
inertia1,
inertia2,
spring,
damper,
u,
])
end
ODESystem(eqs, t; systems = [torque, inertia1, inertia2, spring, damper], name)
end
ODESystem(eqs, t; systems = [torque, inertia1, inertia2, spring, damper], name)
end

model = SystemModel() # Model with load disturbance
model = complete(model)
model_outputs = [model.inertia1.w, model.inertia2.w, model.inertia1.phi, model.inertia2.phi]
model = SystemModel() # Model with load disturbance
model = complete(model)
model_outputs = [model.inertia1.w, model.inertia2.w, model.inertia1.phi, model.inertia2.phi]

@named dmodel = Blocks.StateSpace([0.0], [1.0], [1.0], [0.0]) # An integrating disturbance
@named dmodel = Blocks.StateSpace([0.0], [1.0], [1.0], [0.0]) # An integrating disturbance

@named dist = ModelingToolkit.DisturbanceModel(model.torque.tau.u, dmodel)
(f_oop, f_ip), outersys, dvs, p = ModelingToolkit.add_input_disturbance(model, dist)
@named dist = ModelingToolkit.DisturbanceModel(model.torque.tau.u, dmodel)
(f_oop, f_ip), outersys, dvs, p = ModelingToolkit.add_input_disturbance(model, dist)

@unpack u, d = outersys
matrices, ssys = linearize(outersys, [u, d], model_outputs)
@unpack u, d = outersys
matrices, ssys = linearize(outersys, [u, d], model_outputs)

def = ModelingToolkit.defaults(outersys)
def = ModelingToolkit.defaults(outersys)

# Create a perturbation in the disturbance state
dstate = setdiff(dvs, model_outputs)[]
x_add = ModelingToolkit.varmap_to_vars(merge(Dict(dvs .=> 0), Dict(dstate => 1)), dvs)
# Create a perturbation in the disturbance state
dstate = setdiff(dvs, model_outputs)[]
x_add = ModelingToolkit.varmap_to_vars(merge(Dict(dvs .=> 0), Dict(dstate => 1)), dvs)

x0 = randn(5)
x1 = copy(x0) + x_add # add disturbance state perturbation
u = randn(1)
pn = ModelingToolkit.varmap_to_vars(def, p)
xp0 = f_oop(x0, u, pn, 0)
xp1 = f_oop(x1, u, pn, 0)
x0 = randn(5)
x1 = copy(x0) + x_add # add disturbance state perturbation
u = randn(1)
pn = ModelingToolkit.varmap_to_vars(def, p)
xp0 = f_oop(x0, u, pn, 0)
xp1 = f_oop(x1, u, pn, 0)

@test xp0 matrices.A * x0 + matrices.B * [u; 0]
@test xp1 matrices.A * x1 + matrices.B * [u; 0]
@test xp0 matrices.A * x0 + matrices.B * [u; 0]
@test xp1 matrices.A * x1 + matrices.B * [u; 0]
end

@parameters t
@variables x(t)[1:3] = 0
Expand Down Expand Up @@ -382,4 +387,4 @@ matrices, ssys = linearize(augmented_sys,
# Verify using ControlSystemsBase
# P = ss(A,B,C,0)
# G = ss(matrices...)
# @test sminreal(G[1, 3]) ≈ sminreal(P[1,1])*dist
# @test sminreal(G[1, 3]) ≈ sminreal(P[1,1])*dist
2 changes: 1 addition & 1 deletion test/linearize.jl
Original file line number Diff line number Diff line change
Expand Up @@ -197,7 +197,7 @@ if VERSION >= v"1.8"
D = Differential(t)

@named link1 = Link(; m = 0.2, l = 10, I = 1, g = -9.807)
@named cart = Translational.Mass(; m = 1, s_0 = 0)
@named cart = Translational.Mass(; m = 1, s = 0)
@named fixed = Fixed()
@named force = Force()

Expand Down

0 comments on commit 3bbad02

Please sign in to comment.