diff --git a/src/structural_transformation/pantelides.jl b/src/structural_transformation/pantelides.jl index b6877d65f8..585c4a29d1 100644 --- a/src/structural_transformation/pantelides.jl +++ b/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) diff --git a/src/structural_transformation/symbolics_tearing.jl b/src/structural_transformation/symbolics_tearing.jl index f1cdd7ce9c..bda424f2c4 100644 --- a/src/structural_transformation/symbolics_tearing.jl +++ b/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. diff --git a/src/systems/systemstructure.jl b/src/systems/systemstructure.jl index c0c4a5ff4d..f733ba9be0 100644 --- a/src/systems/systemstructure.jl +++ b/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 diff --git a/test/structural_transformation/utils.jl b/test/structural_transformation/utils.jl index 6dfc107cc9..417d81d3da 100644 --- a/test/structural_transformation/utils.jl +++ b/test/structural_transformation/utils.jl @@ -5,6 +5,7 @@ using SparseArrays using UnPack using ModelingToolkit: t_nounits as t, D_nounits as D, default_toterm using Symbolics: unwrap +using DataInterpolations const ST = StructuralTransformations # Define some variables @@ -282,3 +283,65 @@ 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) + + @testset "Called parameter still has derivative" begin + @component function FilteredInput2(; name, x0 = 0, T = 0.1) + ts = collect(0.0:0.1:10.0) + spline = LinearInterpolation(ts .^ 2, ts) + params = @parameters begin + (k::LinearInterpolation)(..) = spline + T = T + end + vars = @variables begin + x(t) = k(t) + dx(t) = 0 + ddx(t) + end + systems = [] + eqs = [D(x) ~ dx + D(dx) ~ ddx + dx ~ (k(t) - x) / T] + return ODESystem(eqs, t, vars, params; systems, name) + end + + @mtkbuild sys = FilteredInput2() + 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(t)) in vs + end +end