Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

fix: handle derivatives of time-dependent parameters #3493

Open
wants to merge 3 commits into
base: master
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions src/structural_transformation/pantelides.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
3 changes: 2 additions & 1 deletion src/structural_transformation/symbolics_tearing.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down
8 changes: 7 additions & 1 deletion src/systems/systemstructure.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand Down
63 changes: 63 additions & 0 deletions test/structural_transformation/utils.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Loading