Merge pull request #208 from JuliaDynamics/hw/implicitoutputs

handle fully implicit outputs in mtk models

Author: hexaeder

ext/MTKExt.jl

@@ -274,27 +274,41 @@ function generate_io_function(_sys, inputss::Tuple, outputss::Tuple;
     alloutputs = reduce(union, outputss)
 
     missing_inputs = Set{Symbolic}()
+    implicit_outputs = Set{Symbolic}() # fully implicit outputs which do not appear in the equations
     sys = if ModelingToolkit.iscomplete(_sys)
         deepcopy(_sys)
     else
         _openinputs = setdiff(allinputs, Set(parameters(_sys)))
         all_eq_vars = mapreduce(get_variables, union, full_equations(_sys), init=Set{Symbolic}())
         if !(_openinputs ⊆ all_eq_vars)
             missing_inputs = setdiff(_openinputs, all_eq_vars)
-            @warn "The specified inputs ($missing_inputs) do not appear in the equations of the system!"
+            verbose && @warn "The specified inputs ($missing_inputs) do not appear in the equations of the system!"
             _openinputs = setdiff(_openinputs, missing_inputs)
         end
-        structural_simplify(_sys, (_openinputs, alloutputs); simplify=false)[1]
+        _definedoutputs = alloutputs ∩ all_eq_vars
+        if !(Set(_definedoutputs) == Set(alloutputs))
+            implicit_outputs = setdiff(alloutputs, _definedoutputs)
+            verbose && @warn "The specified outputs $implicit_outputs do not appear in the equations of the system!"
+        end
+        structural_simplify(_sys, (_openinputs, _definedoutputs); simplify=false)[1]
     end
 
-    states = unknowns(sys)
     allparams = parameters(sys) # contains inputs!
     @argcheck allinputs ⊆ Set(allparams) ∪ missing_inputs
     params = setdiff(allparams, Set(allinputs))
 
     # extract the main equations and observed equations
     eqs::Vector{Equation} = full_equations(sys)
     fix_metadata!(eqs, sys);
+
+    # assert the ordering of states and equations
+    explicit_states = Symbolic[eq_type(eq)[2] for eq in eqs if !isnothing(eq_type(eq)[2])]
+    implicit_states = setdiff(unknowns(sys), explicit_states) ∪ implicit_outputs
+    states = map(eqs) do eq
+        type = eq_type(eq)
+        isnothing(type[2]) ? pop!(implicit_states) : type[2]
+    end
+
     # check hat there are no rhs differentials in the equations
     if !isempty(rhs_differentials(eqs))
         diffs = rhs_differentials(eqs)
@@ -389,7 +403,7 @@ function generate_io_function(_sys, inputss::Tuple, outputss::Tuple;
     # filter out unnecessary parameters
     var_deps = _all_rhs_symbols(eqs)
     unused_params = Set(setdiff(params, (var_deps ∪ out_deps))) # do not exclud obs_deps
-    if verbose && !isempth(unused_params)
+    if verbose && !isempty(unused_params)
         @info "Parameters $(unused_params) do not appear in equations of f and g and will be marked as unused."
     end
 

test/MTK_test.jl

@@ -233,3 +233,27 @@ v = VertexModel(nestedswing, [:i_r, :i_i], [:u_r, :u_i])
 data = NetworkDynamics.rand_inputs_fg(v)
 b = @b $(NetworkDynamics.compfg(v))($data...)
 @test b.allocs == 0
+
+
+# test fully implicit outputs
+@mtkmodel FullyImplicit begin
+    @variables begin
+        u(t), [description = "Input Variable", input=true]
+        x(t), [description = "Explicit Variable"]
+        y(t), [description = "Implicit Variable, present in equations"]
+        z(t), [description = "fully implicit variable, not present but output", output=true]
+    end
+    @equations begin
+        Dt(x) ~ -x
+        0 ~ sqrt(y+x)
+        0 ~ u # implicitly forces z becaus u(z)
+    end
+end
+@named fullyimplicit = FullyImplicit()
+v = VertexModel(fullyimplicit, [:u], [:z])
+@test v.mass_matrix == LinearAlgebra.Diagonal([1,0,0])
+@test v.sym == [:x, :z, :y]
+
+data = NetworkDynamics.rand_inputs_fg(v)
+b = @b $(NetworkDynamics.compfg(v))($data...)
+@test b.allocs == 0