handle constants in component models

Author: hexaeder

ext/MTKExt.jl

@@ -298,7 +298,7 @@ function generate_io_function(_sys, inputss::Tuple, outputss::Tuple;
     params = setdiff(allparams, Set(allinputs))
 
     # extract the main equations and observed equations
-    eqs::Vector{Equation} = full_equations(sys)
+    eqs::Vector{Equation} = ModelingToolkit.subs_constants(full_equations(sys))
     fix_metadata!(eqs, sys);
 
     # assert the ordering of states and equations
@@ -325,7 +325,8 @@ function generate_io_function(_sys, inputss::Tuple, outputss::Tuple;
     # extract observed equations. They might depend on eachother so resolve them
     obs_subs = Dict(eq.lhs => eq.rhs for eq in observed(sys))
     obseqs = map(observed(sys)) do eq
-        eq.lhs ~ fixpoint_sub(eq.rhs, obs_subs)
+        expanded_rhs = fixpoint_sub(eq.rhs, obs_subs)
+        eq.lhs ~ ModelingToolkit.subs_constants(expanded_rhs)
     end
     fix_metadata!(obseqs, sys);
     # obs can only depend on parameters (including allinputs) or states

src/utils.jl

@@ -118,20 +118,28 @@ end
 
 function rand_inputs_fg(rng, cf)
     @argcheck hasindim(cf) "ComponentModel has no specified input dimensions/syms"
-    du = rand(rng, dim(cf))
+    du = [NaN for _ in 1:dim(cf)]
     u = rand(rng, dim(cf))
     p = rand(rng, pdim(cf))
     ins = Tuple(rand(rng, l) for l in values(indim(cf)))
     if has_external_input(cf)
         ext = rand(rng, extdim(cf))
         ins = (ins..., ext)
     end
-    outs = Tuple(rand(rng, l) for l in values(outdim(cf)))
+    outs = Tuple([NaN for _ in 1:l] for l in values(outdim(cf)))
     t = NaN
-    (outs, du, u, ins, p, t)
+    (outs, du, u, ins, p, t) # fg extrancs outs and ints internally!
 end
 rand_inputs_fg(cf) = rand_inputs_fg(Random.default_rng(), cf)
 
+function rand_inputs_obsf(rng, cf)
+    (_, u, p, ins, p, t) = rand_inputs_fg(rng, cf)
+    N = length(cf.obssym)
+    outs = [NaN for _ in 1:N]
+    (outs, u, ins..., p, t)
+end
+rand_inputs_obsf(cf) = rand_inputs_obsf(Random.default_rng(), cf)
+
 
 # abstract symbolic index types
 abstract type SymbolicIndex{C,S} end

test/MTK_test.jl

@@ -257,3 +257,38 @@ v = VertexModel(fullyimplicit, [:u], [:z])
 data = NetworkDynamics.rand_inputs_fg(v)
 b = @b $(NetworkDynamics.compfg(v))($data...)
 @test b.allocs == 0
+
+@testset "Test constants in MTK models" begin
+    @mtkmodel DQSwing_Constants begin
+        @extend DQBus()
+        @variables begin
+            ω(t) = 0.0, [description = "Rotor frequency"]
+            θ(t) = 0.0, [description = "Rotor angle"]
+            Pel(t), [description = "Electrical Power"]
+        end
+        @constants begin
+            V = 1, [description = "Voltage magnitude"]
+            useless = 0
+        end
+        @parameters begin
+            M = 1, [description = "Inertia"]
+            D = 0.1, [description = "Damping"]
+            Pmech, [description = "Mechanical Power"]
+        end
+        @equations begin
+            Dt(θ) ~ ω
+            Dt(ω) ~ 1/M * (Pmech - D*ω + Pel)
+            Pel ~ real((u_r + im*u_i) * (i_r - im*i_i)) + useless
+            u_r ~ V*cos(θ)
+            u_i ~ V*sin(θ)
+        end
+    end
+    @named withconst = DQSwing_Constants()
+    v = VertexModel(withconst, [:i_r, :i_i], [:u_r, :u_i])
+
+    data = NetworkDynamics.rand_inputs_fg(v)
+    NetworkDynamics.compfg(v)(data...) # no error
+
+    data = NetworkDynamics.rand_inputs_obsf(v)
+    v.obsf(data...) # no error
+end