From 4d9a4f14d4dc9fd36bac0a62fc722d1a054a3110 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Mon, 24 Mar 2025 13:21:28 +0530 Subject: [PATCH 1/7] 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 d27e5c93a1..7ad2ceb587 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 6dfc107cc9..fd7963626d 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 18e17959a903ad07bee307bf75eb2fb06f5933f8 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Mon, 24 Mar 2025 13:39:03 +0530 Subject: [PATCH 2/7] 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 fd7963626d..cfee3c260a 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 4dbd08f11022a37e2f6e958d1991d38972e6bc76 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Mon, 24 Mar 2025 15:21:45 +0530 Subject: [PATCH 3/7] 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 cfee3c260a..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 From 105dd7295a451c2d98f0e24c2ed109b5d8ce1f55 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Tue, 22 Apr 2025 13:49:24 +0530 Subject: [PATCH 4/7] 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 7ad2ceb587..7704e85938 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 417d81d3da..ded5581b94 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 e660b0b840f806e8a7683853f1fa164704c6ab18 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Tue, 22 Apr 2025 18:13:46 +0530 Subject: [PATCH 5/7] 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 7561632666da62901d1d8df03da77a1f9c707aa5 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Wed, 23 Apr 2025 00:42:06 +0530 Subject: [PATCH 6/7] fix: handle scalarized time-dependent array parameter derivatives --- 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 7704e85938..91681984d7 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 046d51c97bb88b1a954a4a56e9518f78b7153531 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Wed, 23 Apr 2025 00:42:27 +0530 Subject: [PATCH 7/7] feat: add support for obtaining derivatives of CoSimulation FMU outputs --- ext/MTKFMIExt.jl | 87 +++++++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 82 insertions(+), 5 deletions(-) diff --git a/ext/MTKFMIExt.jl b/ext/MTKFMIExt.jl index 5cfe9a82ef..04cf96a211 100644 --- a/ext/MTKFMIExt.jl +++ b/ext/MTKFMIExt.jl @@ -249,15 +249,21 @@ function MTK.FMIComponent(::Val{Ver}; fmu = nothing, tolerance = 1e-6, FMI3CSFunctor(state_value_references, output_value_references) end @parameters (functor::(typeof(_functor)))(..)[1:(length(__mtk_internal_u) + length(__mtk_internal_o))] = _functor - # for co-simulation, we need to ensure the output buffer is solved for - # during initialization + + diffeqs = Equation[] for (i, x) in enumerate(collect(__mtk_internal_o)) + # for co-simulation, we need to ensure the output buffer is solved for + # during initialization push!(initialization_eqs, x ~ functor( - wrapper, __mtk_internal_u, __mtk_internal_x, __mtk_internal_p, t)[i]) - end + wrapper, (__mtk_internal_u), __mtk_internal_x, __mtk_internal_p, t)[i]) - diffeqs = Equation[] + # also add equations for output derivatives + push!(diffeqs, + D(x) ~ term( + getOutputDerivative, functor, wrapper, i, 1, collect(__mtk_internal_u), + __mtk_internal_x, __mtk_internal_p, t; type = Real)) + end # use `ImperativeAffect` for instance management here cb_observed = (; inputs = __mtk_internal_x, params = copy(params), @@ -739,6 +745,15 @@ struct FMI2CSFunctor The value references of output variables in the FMU. """ output_value_references::Vector{FMI.fmi2ValueReference} + """ + Simply a buffer to store the order of output derivative required from + `getRealOutputderivatives` and avoid unnecessary allocations. + """ + output_derivative_order_buffer::Vector{FMI.fmi2Integer} +end + +function FMI2CSFunctor(svref, ovref) + FMI2CSFunctor(svref, ovref, FMI.fmi2Integer[1]) end function (fn::FMI2CSFunctor)(wrapper::FMI2InstanceWrapper, states, inputs, params, t) @@ -764,6 +779,41 @@ end ndims = 1 end +""" + $(TYPEDSIGNATURES) + +Calculate the `order` order derivative of the `var`th output of the FMU. +""" +function getOutputDerivative(fn::FMI2CSFunctor, wrapper::FMI2InstanceWrapper, var::Int, + order::FMI.fmi2Integer, states, inputs, params, t) + states = states isa SubArray ? copy(states) : states + inputs = inputs isa SubArray ? copy(inputs) : inputs + params = params isa SubArray ? copy(params) : params + instance = get_instance_CS!(wrapper, states, inputs, params, t) + fn.output_derivative_order_buffer[] = order + return FMI.fmi2GetRealOutputDerivatives( + instance, fn.output_value_references[var], fn.output_derivative_order_buffer) +end + +# @register_symbolic getOutputDerivative(fn::FMI2CSFunctor, w::FMI2InstanceWrapper, var::Int, order::FMI.fmi2Integer, states::Vector{<:Real}, inputs::Vector{<:Real}, params::Vector{<:Real}, t::Real) + +# HACK-ish for allowing higher order output derivatives +# The first `D(output)` will result in a `getOutputDerivatives` expression. +# Subsequent differentiations of this expression will expand to +# `Σ_{i = 1:8} Differential(args[i])(getOutputDerivative(args...)) * D(args[i])` +# using the chain rule. `i = 1:4` are not time-dependent (or real). We define +# the derivatives for `i = 5:7` to be zero, and the derivative for `i = 8` (w.r.t `t`) +# to be the same `getOutputDerivative` call but with the order increased. This results +# in `D(output) = getOutputDerivative(fn, w, var, order + 1, states, inputs, params, t) * 1` +# which is exactly what we want. +for i in 1:7 + @eval Symbolics.derivative(::typeof(getOutputDerivative), args::NTuple{8, Any}, ::Val{$i}) = 0 +end +function Symbolics.derivative(::typeof(getOutputDerivative), args::NTuple{8, Any}, ::Val{8}) + term(getOutputDerivative, args[1], args[2], args[3], + args[4] + 1, args[5], args[6], args[7], args[8]) +end + """ $(TYPEDSIGNATURES) @@ -848,6 +898,15 @@ struct FMI3CSFunctor The value references of output variables in the FMU. """ output_value_references::Vector{FMI.fmi3ValueReference} + """ + Simply a buffer to store the order of output derivative required from + `getRealOutputderivatives` and avoid unnecessary allocations. + """ + output_derivative_order_buffer::Vector{FMI.fmi3Int32} +end + +function FMI3CSFunctor(svref, ovref) + FMI3CSFunctor(svref, ovref, FMI.fmi3Int32[1]) end function (fn::FMI3CSFunctor)(wrapper::FMI3InstanceWrapper, states, inputs, params, t) @@ -871,6 +930,24 @@ end ndims = 1 end +""" + $(TYPEDSIGNATURES) + +Calculate the `order` order derivative of the `var`th output of the FMU. +""" +function getOutputDerivative(fn::FMI3CSFunctor, wrapper::FMI3InstanceWrapper, var::Int, + order::FMI.fmi3Int32, states, inputs, params, t) + states = states isa SubArray ? copy(states) : states + inputs = inputs isa SubArray ? copy(inputs) : inputs + params = params isa SubArray ? copy(params) : params + instance = get_instance_CS!(wrapper, states, inputs, params, t) + fn.output_derivative_order_buffer[] = order + return FMI.fmi3GetOutputDerivatives( + instance, fn.output_value_references[var], fn.output_derivative_order_buffer) +end + +# @register_symbolic getOutputDerivative(fn::FMI3CSFunctor, w::FMI3InstanceWrapper, var::Int, order::FMI.fmi3Int32, states::Vector{<:Real}, inputs::Vector{<:Real}, params::Vector{<:Real}, t::Real) false + """ $(TYPEDSIGNATURES) """