fix: handle derivatives of time-dependent parameters

Author: AayushSabharwal

src/structural_transformation/pantelides.jl

@@ -54,6 +54,7 @@ function pantelides_reassemble(state::TearingState, var_eq_matching)
             D(eq.lhs)
         end
         rhs = ModelingToolkit.expand_derivatives(D(eq.rhs))
+        rhs = fast_substitute(rhs, state.param_derivative_map)
         substitution_dict = Dict(x.lhs => x.rhs
         for x in out_eqs if x !== nothing && x.lhs isa Symbolic)
         sub_rhs = substitute(rhs, substitution_dict)

src/structural_transformation/symbolics_tearing.jl

@@ -65,7 +65,8 @@ function eq_derivative!(ts::TearingState{ODESystem}, ieq::Int; kwargs...)
 
     sys = ts.sys
     eq = equations(ts)[ieq]
-    eq = 0 ~ ModelingToolkit.derivative(eq.rhs - eq.lhs, get_iv(sys))
+    eq = 0 ~ fast_substitute(
+        ModelingToolkit.derivative(eq.rhs - eq.lhs, get_iv(sys)), ts.param_derivative_map)
     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,6 +207,7 @@ mutable struct TearingState{T <: AbstractSystem} <: AbstractTearingState{T}
     fullvars::Vector
     structure::SystemStructure
     extra_eqs::Vector
+    param_derivative_map::Dict{BasicSymbolic, Real}
 end
 
 TransformationState(sys::AbstractSystem) = TearingState(sys)
@@ -264,6 +265,7 @@ function TearingState(sys; quick_cancel = false, check = true)
     var2idx = Dict{Any, Int}()
     symbolic_incidence = []
     fullvars = []
+    param_derivative_map = Dict{BasicSymbolic, Real}()
     var_counter = Ref(0)
     var_types = VariableType[]
     addvar! = let fullvars = fullvars, var_counter = var_counter, var_types = var_types
@@ -295,6 +297,10 @@ 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
+                end
                 continue
             end
             v = scalarize(v)
@@ -439,7 +445,7 @@ function TearingState(sys; quick_cancel = false, check = true)
     ts = TearingState(sys, fullvars,
         SystemStructure(complete(var_to_diff), complete(eq_to_diff),
             complete(graph), nothing, var_types, sys isa AbstractDiscreteSystem),
-        Any[])
+        Any[], param_derivative_map)
     if sys isa DiscreteSystem
         ts = shift_discrete_system(ts)
     end

test/structural_transformation/utils.jl

@@ -282,3 +282,33 @@ end
         @test length(mapping) == 3
     end
 end
+
+@testset "Issue#3480: Derivatives of time-dependent parameters" begin
+    @component function FilteredInput(; 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]
+        return ODESystem(eqs, t, vars, params; systems, name)
+    end
+
+    @mtkbuild sys = FilteredInput()
+    vs = Set()
+    for eq in equations(sys)
+        ModelingToolkit.vars!(vs, eq)
+    end
+    for eq in observed(sys)
+        ModelingToolkit.vars!(vs, eq)
+    end
+
+    @test !(D(sys.k) in vs)
+end