From d4d23a380f60c00227973ce9fb0714390c2f8e94 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Mon, 24 Mar 2025 13:21:28 +0530 Subject: [PATCH 1/8] fix: handle derivatives of time-dependent parameters --- src/structural_transformation/pantelides.jl | 1 + .../symbolics_tearing.jl | 4 ++- src/systems/systemstructure.jl | 8 ++++- test/structural_transformation/utils.jl | 30 +++++++++++++++++++ 4 files changed, 41 insertions(+), 2 deletions(-) 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 552c6d13c3..31307d132f 100644 --- a/src/structural_transformation/symbolics_tearing.jl +++ b/src/structural_transformation/symbolics_tearing.jl @@ -65,7 +65,9 @@ function eq_derivative!(ts::TearingState{ODESystem}, ieq::Int; kwargs...) sys = ts.sys eq = equations(ts)[ieq] - eq = 0 ~ Symbolics.derivative(eq.rhs - eq.lhs, get_iv(sys); throw_no_derivative = true) + eq = 0 ~ fast_substitute( + ModelingToolkit.derivative( + eq.rhs - eq.lhs, get_iv(sys); throw_no_derivative = true), 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 73cdaba255..d31a45a653 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) @@ -438,7 +444,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 8894f0bbc9..aad43f9616 100644 --- a/test/structural_transformation/utils.jl +++ b/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 From b7dbe5138d7e1fed2bc42ae8e4ba2019fc1ecec4 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Mon, 24 Mar 2025 13:39:03 +0530 Subject: [PATCH 2/8] test: test called parameters are still differentiated --- test/structural_transformation/utils.jl | 32 +++++++++++++++++++++++++ 1 file changed, 32 insertions(+) diff --git a/test/structural_transformation/utils.jl b/test/structural_transformation/utils.jl index aad43f9616..2565cb9626 100644 --- a/test/structural_transformation/utils.jl +++ b/test/structural_transformation/utils.jl @@ -311,4 +311,36 @@ end 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 From fdf3cb10a106506f6120f6e345765ac01bb19927 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Mon, 24 Mar 2025 15:21:45 +0530 Subject: [PATCH 3/8] test: import `DataInterpolations` in test --- test/structural_transformation/utils.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/test/structural_transformation/utils.jl b/test/structural_transformation/utils.jl index 2565cb9626..7c0cb633c7 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 From f1c9c164b9eb81847a8f2cf3672a554869215ac0 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Tue, 22 Apr 2025 13:49:24 +0530 Subject: [PATCH 4/8] fix: require explicitly specifying discrete variable derivatives --- .../symbolics_tearing.jl | 14 ++++++++++ src/systems/systemstructure.jl | 26 +++++++++++++++---- test/structural_transformation/utils.jl | 23 +++++++++++++++- 3 files changed, 57 insertions(+), 6 deletions(-) diff --git a/src/structural_transformation/symbolics_tearing.jl b/src/structural_transformation/symbolics_tearing.jl index 31307d132f..548c7da519 100644 --- a/src/structural_transformation/symbolics_tearing.jl +++ b/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. diff --git a/src/systems/systemstructure.jl b/src/systems/systemstructure.jl index d31a45a653..61ae06aab3 100644 --- a/src/systems/systemstructure.jl +++ b/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() diff --git a/test/structural_transformation/utils.jl b/test/structural_transformation/utils.jl index 7c0cb633c7..996f2e7b83 100644 --- a/test/structural_transformation/utils.jl +++ b/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) From 8e3f6e379775d4e10173f17795117c7b36190e6c Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Tue, 22 Apr 2025 18:13:46 +0530 Subject: [PATCH 5/8] test: fix state selection test --- test/state_selection.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/state_selection.jl b/test/state_selection.jl index b8404d1f26..a8d3e57773 100644 --- a/test/state_selection.jl +++ b/test/state_selection.jl @@ -2,7 +2,7 @@ using ModelingToolkit, OrdinaryDiffEq, Test using ModelingToolkit: t_nounits as t, D_nounits as D sts = @variables x1(t) x2(t) x3(t) x4(t) -params = @parameters u1(t) u2(t) u3(t) u4(t) +params = @parameters u1 u2 u3 u4 eqs = [x1 + x2 + u1 ~ 0 x1 + x2 + x3 + u2 ~ 0 x1 + D(x3) + x4 + u3 ~ 0 From 4935f3b32def60b03338ae9baed4967d1e3245ca Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Thu, 24 Apr 2025 17:59:59 +0530 Subject: [PATCH 6/8] fix: handle derivatives of time-dependent array parameters --- src/systems/systemstructure.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/systems/systemstructure.jl b/src/systems/systemstructure.jl index 61ae06aab3..8ec731db5e 100644 --- a/src/systems/systemstructure.jl +++ b/src/systems/systemstructure.jl @@ -256,7 +256,8 @@ 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) + (operation(p) === getindex && is_time_dependent_parameter(arguments(p)[1], iv) || + (args = arguments(p); length(args)) == 1 && isequal(only(args), iv)) end function TearingState(sys; quick_cancel = false, check = true) From 4d8463cbeb9134bcb49a95d8784301ae33be8a40 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Thu, 24 Apr 2025 18:00:23 +0530 Subject: [PATCH 7/8] refactor: default parameter derivatives to zero, allow opting-out --- src/systems/systemstructure.jl | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) diff --git a/src/systems/systemstructure.jl b/src/systems/systemstructure.jl index 8ec731db5e..5927d3ebad 100644 --- a/src/systems/systemstructure.jl +++ b/src/systems/systemstructure.jl @@ -292,8 +292,14 @@ function TearingState(sys; quick_cancel = false, check = true) 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 + # parameter derivatives are opted out by specifying `D(p) ~ missing`, but + # we want to store `nothing` in the map because that means `fast_substitute` + # will ignore the rule. We will this identify the presence of `eq′.lhs` in + # the differentiated expression and error. + param_derivative_map[eq′.lhs] = coalesce(eq′.rhs, nothing) eqs_to_retain[i] = false + # change the equation if the RHS is `missing` so the rest of this loop works + eq′ = eq′.lhs ~ coalesce(eq′.rhs, 0.0) end if _iszero(eq′.lhs) rhs = quick_cancel ? quick_cancel_expr(eq′.rhs) : eq′.rhs @@ -311,9 +317,9 @@ function TearingState(sys; quick_cancel = false, check = true) (iscall(_var) && isparameter(operation(_var)) || isconstant(_var)) 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 + # Parameter derivatives default to zero - they stay constant + # between callbacks + param_derivative_map[Differential(iv)(_var)] = 0.0 end continue end From a2b47451946b5e64c379f9463173654ae632cd7b Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Thu, 24 Apr 2025 18:00:30 +0530 Subject: [PATCH 8/8] test: test new parameter derivative behavior --- test/structural_transformation/utils.jl | 33 ++++++++++++++++++++++--- 1 file changed, 29 insertions(+), 4 deletions(-) diff --git a/test/structural_transformation/utils.jl b/test/structural_transformation/utils.jl index 996f2e7b83..6c5a564e9f 100644 --- a/test/structural_transformation/utils.jl +++ b/test/structural_transformation/utils.jl @@ -302,7 +302,25 @@ end return ODESystem(eqs, t, vars, params; systems, name) end - @component function FilteredInputFix(; name, x0 = 0, T = 0.1) + @component function FilteredInputExplicit(; name, x0 = 0, T = 0.1) + params = @parameters begin + k(t)[1:1] = [x0] + T = T + end + vars = @variables begin + x(t) = k + dx(t) = 0 + ddx(t) + end + systems = [] + eqs = [D(x) ~ dx + D(dx) ~ ddx + D(k[1]) ~ 1.0 + dx ~ (k[1] - x) / T] + return ODESystem(eqs, t, vars, params; systems, name) + end + + @component function FilteredInputErr(; name, x0 = 0, T = 0.1) params = @parameters begin k(t) = x0 T = T @@ -316,14 +334,14 @@ end eqs = [D(x) ~ dx D(dx) ~ ddx dx ~ (k - x) / T - D(k) ~ 0] + D(k) ~ missing] return ODESystem(eqs, t, vars, params; systems, name) end - @named sys = FilteredInput() + @named sys = FilteredInputErr() @test_throws ["derivative of discrete variable", "k(t)"] structural_simplify(sys) - @mtkbuild sys = FilteredInputFix() + @mtkbuild sys = FilteredInput() vs = Set() for eq in equations(sys) ModelingToolkit.vars!(vs, eq) @@ -334,6 +352,13 @@ end @test !(D(sys.k) in vs) + @mtkbuild sys = FilteredInputExplicit() + obsfn1 = ModelingToolkit.build_explicit_observed_function(sys, sys.ddx) + obsfn2 = ModelingToolkit.build_explicit_observed_function(sys, sys.dx) + u = [1.0] + p = MTKParameters(sys, [sys.k => [2.0], sys.T => 3.0]) + @test obsfn1(u, p, 0.0) ≈ (1 - obsfn2(u, p, 0.0)) / 3.0 + @testset "Called parameter still has derivative" begin @component function FilteredInput2(; name, x0 = 0, T = 0.1) ts = collect(0.0:0.1:10.0)