fix: require explicitly specifying discrete variable derivatives

Author: AayushSabharwal

src/structural_transformation/symbolics_tearing.jl

@@ -68,6 +68,20 @@ function eq_derivative!(ts::TearingState{ODESystem}, ieq::Int; kwargs...)
     eq = 0 ~ fast_substitute(
         ModelingToolkit.derivative(
             eq.rhs - eq.lhs, get_iv(sys); throw_no_derivative = true), ts.param_derivative_map)
+
+    vs = ModelingToolkit.vars(eq.rhs)
+    for v in vs
+        # parameters with unknown derivatives have a value of `nothing` in the map,
+        # so use `missing` as the default.
+        get(ts.param_derivative_map, v, missing) === nothing || continue
+        _original_eq = equations(ts)[ieq]
+        error("""
+        Encountered derivative of discrete variable `$(only(arguments(v)))` when \
+        differentiating equation `$(_original_eq)`. This may indicate a model error or a \
+        missing equation of the form `$v ~ ...` that defines this derivative.
+        """)
+    end
+
     push!(equations(ts), eq)
     # Analyze the new equation and update the graph/solvable_graph
     # First, copy the previous incidence and add the derivative terms.

src/systems/systemstructure.jl

@@ -207,7 +207,7 @@ mutable struct TearingState{T <: AbstractSystem} <: AbstractTearingState{T}
     fullvars::Vector
     structure::SystemStructure
     extra_eqs::Vector
-    param_derivative_map::Dict{BasicSymbolic, Real}
+    param_derivative_map::Dict{BasicSymbolic, Any}
 end
 
 TransformationState(sys::AbstractSystem) = TearingState(sys)
@@ -254,6 +254,11 @@ function Base.push!(ev::EquationsView, eq)
     push!(ev.ts.extra_eqs, eq)
 end
 
+function is_time_dependent_parameter(p, iv)
+    return iv !== nothing && isparameter(p) && iscall(p) &&
+           (args = arguments(p); length(args)) == 1 && isequal(only(args), iv)
+end
+
 function TearingState(sys; quick_cancel = false, check = true)
     sys = flatten(sys)
     ivs = independent_variables(sys)
@@ -265,7 +270,7 @@ function TearingState(sys; quick_cancel = false, check = true)
     var2idx = Dict{Any, Int}()
     symbolic_incidence = []
     fullvars = []
-    param_derivative_map = Dict{BasicSymbolic, Real}()
+    param_derivative_map = Dict{BasicSymbolic, Any}()
     var_counter = Ref(0)
     var_types = VariableType[]
     addvar! = let fullvars = fullvars, var_counter = var_counter, var_types = var_types
@@ -278,11 +283,17 @@ function TearingState(sys; quick_cancel = false, check = true)
 
     vars = OrderedSet()
     varsvec = []
+    eqs_to_retain = trues(length(eqs))
     for (i, eq′) in enumerate(eqs)
         if eq′.lhs isa Connection
             check ? error("$(nameof(sys)) has unexpanded `connect` statements") :
             return nothing
         end
+        if iscall(eq′.lhs) && (op = operation(eq′.lhs)) isa Differential &&
+           isequal(op.x, iv) && is_time_dependent_parameter(only(arguments(eq′.lhs)), iv)
+            param_derivative_map[eq′.lhs] = eq′.rhs
+            eqs_to_retain[i] = false
+        end
         if _iszero(eq′.lhs)
             rhs = quick_cancel ? quick_cancel_expr(eq′.rhs) : eq′.rhs
             eq = eq′
@@ -297,9 +308,11 @@ function TearingState(sys; quick_cancel = false, check = true)
             any(isequal(_var), ivs) && continue
             if isparameter(_var) ||
                (iscall(_var) && isparameter(operation(_var)) || isconstant(_var))
-                if iv !== nothing && isparameter(_var) && iscall(_var) &&
-                   (args = arguments(_var); length(args)) == 1 && isequal(only(args), iv)
-                    param_derivative_map[Differential(iv)(_var)] = 0.0
+                if is_time_dependent_parameter(_var, iv) &&
+                   !haskey(param_derivative_map, Differential(iv)(_var))
+                    # default to `nothing` since it is ignored during substitution,
+                    # so `D(_var)` is retained in the expression.
+                    param_derivative_map[Differential(iv)(_var)] = nothing
                 end
                 continue
             end
@@ -357,6 +370,9 @@ function TearingState(sys; quick_cancel = false, check = true)
             eqs[i] = eqs[i].lhs ~ rhs
         end
     end
+    eqs = eqs[eqs_to_retain]
+    neqs = length(eqs)
+    symbolic_incidence = symbolic_incidence[eqs_to_retain]
 
     ### Handle discrete variables
     lowest_shift = Dict()

test/structural_transformation/utils.jl

@@ -302,7 +302,28 @@ end
         return ODESystem(eqs, t, vars, params; systems, name)
     end
 
-    @mtkbuild sys = FilteredInput()
+    @component function FilteredInputFix(; name, x0 = 0, T = 0.1)
+        params = @parameters begin
+            k(t) = x0
+            T = T
+        end
+        vars = @variables begin
+            x(t) = k
+            dx(t) = 0
+            ddx(t)
+        end
+        systems = []
+        eqs = [D(x) ~ dx
+               D(dx) ~ ddx
+               dx ~ (k - x) / T
+               D(k) ~ 0]
+        return ODESystem(eqs, t, vars, params; systems, name)
+    end
+
+    @named sys = FilteredInput()
+    @test_throws ["derivative of discrete variable", "k(t)"] structural_simplify(sys)
+
+    @mtkbuild sys = FilteredInputFix()
     vs = Set()
     for eq in equations(sys)
         ModelingToolkit.vars!(vs, eq)