Calculating derivatives of observables returns nonsense

[orebas]

The below MWE is an attempt to programmatically take an ODE system with some observables defined, and add the derivatives of the observables into the system. I construct an array of successive derivatives of the observable, and pass that to ODESystem().
Both ODESystem and the solver call proceed as if nothing is wrong.

Moreover, I'm able to query the solution object for these derivatives: the first derivative seems to be calculated fine, but from the second order derivative and higher, some strange object like "Differential(t)(32.0)" or "Differential(t)(Differential(t)(32.0))" is returned. Moreover, the actual number contained in this is wrong (it's just a copy of the first derivative).

Probably these inputs should be rejected by ODESystem(), but I would love advice on a solution for how to get this from the ODE solver.

Output of MWE:

1×5 Matrix{Num}:
 d_obs₁ˏ₁  d_obs₁ˏ₂  d_obs₁ˏ₃  d_obs₁ˏ₄  d_obs₁ˏ₅
5-element Vector{Vector{Equation}}:
 [d_obs₁ˏ₁ ~ Differential(t)(x1(t))]
 [d_obs₁ˏ₂ ~ Differential(t)(Differential(t)(x1(t)))]
 [d_obs₁ˏ₃ ~ Differential(t)(Differential(t)(Differential(t)(x1(t))))]
 [d_obs₁ˏ₄ ~ Differential(t)(Differential(t)(Differential(t)(Differential(t)(x1(t)))))]
 [d_obs₁ˏ₅ ~ Differential(t)(Differential(t)(Differential(t)(Differential(t)(Differential(t)(x1(t))))))]
t 0.0
0.0
0.0
Differential(t)(0.0)
Differential(t)(Differential(t)(0.0))
Differential(t)(Differential(t)(Differential(t)(0.0)))
Differential(t)(Differential(t)(Differential(t)(Differential(t)(0.0))))
t 1.0
0.9999999999999871
4.0
Differential(t)(4.0)
Differential(t)(Differential(t)(4.0))
Differential(t)(Differential(t)(Differential(t)(4.0)))
Differential(t)(Differential(t)(Differential(t)(Differential(t)(4.0))))
t 2.0
15.999999999966947
32.0
Differential(t)(32.0)
Differential(t)(Differential(t)(32.0))
Differential(t)(Differential(t)(Differential(t)(32.0)))
Differential(t)(Differential(t)(Differential(t)(Differential(t)(32.0))))

MWE:

using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D
using OrdinaryDiffEq
using SymbolicIndexingInterface

function lotka_volterra()
	parameters = @parameters a b c d
	states = @variables x1(t) x2(t) #xp1(t) xp2(t)
	observables = @variables y1(t) y2(t)
	p_true = [0.0000, 0.0]
	ic_true = [0.0]
	time_span = (0.0, 2.0)
	equations = [
		D(x1) ~ x1 * a * b + t * t * t * 4.0,
	]
	measured_quantities = [y1 ~ x1]

	n_observables = length(measured_quantities)
	nderivs = 5
	ObservableDerivatives = Symbolics.variables(:d_obs, 1:n_observables, 1:nderivs)  #T = Symbolics.FnType
	display(ObservableDerivatives)
	# Initialize vector of derivative equations for each order
	SymbolicDerivs = Vector{Vector{Equation}}(undef, nderivs)

	# First derivatives
	SymbolicDerivs[1] = [ObservableDerivatives[i, 1] ~ expand_derivatives(D(measured_quantities[i].rhs)) for i in 1:n_observables]

	# Higher order derivatives
	for j in 2:nderivs
		SymbolicDerivs[j] = [ObservableDerivatives[i, j] ~ expand_derivatives(D(SymbolicDerivs[j-1][i].rhs)) for i in 1:n_observables]
	end

	display(SymbolicDerivs)

	# Add all derivative equations to measured quantities
	append!(measured_quantities, vcat(SymbolicDerivs...))



	# Create the ODESystem
	@named sys = ODESystem(equations, t; observed = measured_quantities)

	# Convert to ODEProblem
	prob = ODEProblem(structural_simplify(sys), ic_true, time_span, p_true)

	# Solve the system
	sol = solve(prob, Vern9())

	for myt in (0.0, 1.0, 2.0)
		y1_val = sol(myt, idxs = y1)
		dy1_val = sol(myt, idxs = ObservableDerivatives[1, 1])
		ddy1_val = sol(myt, idxs = ObservableDerivatives[1, 2])
		dddy1_val = sol(myt, idxs = ObservableDerivatives[1, 3])
		ddddy1_val = sol(myt, idxs = ObservableDerivatives[1, 4])
		dddddy1_val = sol(myt, idxs = ObservableDerivatives[1, 5])


		println("t ", myt)
		display(y1_val)
		display(dy1_val)
		display(ddy1_val)
		display(dddy1_val)
		display(ddddy1_val)
		display(dddddy1_val)
	end
end

lotka_volterra()

Versions: MTK 9.60.0, Symbolics 6.22.1, SymbolicUtils 3.10.0, OrdinaryDiffEq 6.90.1

[hersle]

One workaround I have used is if you can "anticipate" beforehand that derivatives are wanted from the solver output, you can define observed equations for them like this:

using ModelingToolkit, OrdinaryDiffEq, Plots
using ModelingToolkit: t_nounits as t, D_nounits as D

@parameters a b
@variables x(t) ẋ(t) ẍ(t) x⃛(t) x⃜(t)
eqs = [
    D(x) ~ a * b * x + 4 * t^3
    ẋ ~ D(x)
    ẍ ~ D(ẋ)
    x⃛ ~ D(ẍ)
    x⃜ ~ D(x⃛)
]
@named sys = ODESystem(eqs, t)
sys = structural_simplify(sys)
prob = ODEProblem(sys, [x => 0.0], (0.0, 2.0), [a => 1.0, b => 2.0])
sol = solve(prob, Tsit5())
plot(sol, idxs = [x, ẋ, ẍ, x⃛, x⃜])

This avoids using the D operator when indexing the solution (which creates the D(0.0) nonsense), and relies on internal derivative expansion and dummy variables (see observed(sys)). It works well for this simple example, and I have found it semi-reliable for more complex problems. But it requires manual intervention, generally for every single quantity, so I don't think it's a good general solution. I think it would be better if it was somehow possible to automatically calculate derivatives automatically from the solver without "anticipating" them when defining the system.

[orebas]

I solved my problem in a slightly different way; you can repeatedly substitute yourself using the state equations. Below is a moderately generic function to do that. Thank you for your suggestion, it seems not too hard to programmatically do either.

In any case, it would be nice if MTK either had a built-in for this, or rejected the inputs in OP with an error, or accepted them and "did the right thing".

"""
	calculate_observable_derivatives(equations, measured_quantities, nderivs=5)

Calculate symbolic derivatives of observables up to the specified order using ModelingToolkit.
Returns the expanded measured quantities with derivatives and the derivative variables.
"""
function calculate_observable_derivatives(equations, measured_quantities, nderivs = 5)
	# Create equation dictionary for substitution
	equation_dict = Dict(eq.lhs => eq.rhs for eq in equations)

	n_observables = length(measured_quantities)

	# Create symbolic variables for derivatives
	ObservableDerivatives = Symbolics.variables(:d_obs, 1:n_observables, 1:nderivs)

	# Initialize vector to store derivative equations
	SymbolicDerivs = Vector{Vector{Equation}}(undef, nderivs)

	# Calculate first derivatives
	SymbolicDerivs[1] = [ObservableDerivatives[i, 1] ~ substitute(expand_derivatives(D(measured_quantities[i].rhs)), equation_dict) for i in 1:n_observables]

	# Calculate higher order derivatives
	for j in 2:nderivs
		SymbolicDerivs[j] = [ObservableDerivatives[i, j] ~ substitute(expand_derivatives(D(SymbolicDerivs[j-1][i].rhs)), equation_dict) for i in 1:n_observables]
	end

	# Create new measured quantities with derivatives
	expanded_measured_quantities = copy(measured_quantities)
	append!(expanded_measured_quantities, vcat(SymbolicDerivs...))

	return expanded_measured_quantities, ObservableDerivatives
end