Generating derived variables from ModelingToolkit.jl - variable of interest eliminated during structural_simplify

Hi everyone,

I’m trying to generate variables that are functions of my state variables in an ODE. Here’s my example.

using ModelingToolkit
using OrdinaryDiffEq

@parameters t β c γ
@variables S(t) I(t) R(t) λ(t)
D = Differential(t)
eqs = [D(S) ~ -λ*S,
       D(I) ~ λ*S-γ*I,
       D(R) ~ γ*I,
       λ ~ β*c*I/N];
@named sys = ODESystem(eqs) # 4 equations: S,I,R,λ
simpsys = structural_simplify(sys) # 3 equations: S,I,R

Defining λ not only makes the ODEs easier to read, but it’s also a variable that I want access to, in order to plug it in to other models. λ gets eliminated (as it should) during structural_simplify, but then I lose the ability to access it. Is there any way to run the ODE problem and get back all the variables (S,I,R,λ)?

\lambda has likely been moved to observed(simsys). These variables are still accessible in the solution when you have simulated, but it’s also possible to build a function that computes the observed variables given the states.

The following should do what you want

outputs = [λ]
syso, _ = ModelingToolkit.io_preprocessing(sys, [], outputs;)
obsf = ModelingToolkit.build_explicit_observed_function(syso, outputs)

lambda = obsf[1](x, p, t) # where x is the state of the simplified system

Or perhaps better yet, without using internals, construct an ODEProblem and then

λfun = prob.f.observed(λ)
lambda = λfun(x, p, t) 

In both of these approaches, x and p need to be numerical arrays in the correct order. The order is given by states(simpsys) and parameters(simpsys)

It’s worth mentioning that sol[λ] will still get you the solution for lambda, even if it’s eliminated

Thanks for these tips! Is there any way to add those observed values to the solution, so I can call sol_ode(10.0) (say) and get the states and the derived values from the interpolation?

It’s already in there, I can never remember the syntax to obtain the interpolated value, but I know it is possible :sweat_smile:

λfun = prob.f.observed(λ)
lambda = λfun(x, p, t) 

state_interp = sol(10.0)
complete_interp = [state_interp..., λfun(state_interp, [], 10.0)]

This is the first time I have noticed that the observed aren’t automatically included in the interp, this is probably an oversight. Please post an issue to SciMLBase, and I may get to it soon.

sol_ode(10.0, idxs = λ)

And note you can create new derived variables on the fly and it will compute it.

sol_ode(10.0, idxs = λ + R^2)

Please don’t use that :sweat_smile: . I mean, Fredrik can, but it’s internals most people shouldn’t know and it’s not guaranteed to keep the same syntax.

Thanks Chris!