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

Specific Domains Modelling & Simulations
1

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)
N=S+I+R
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,λ)?

2

\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)

1 Like
3

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

4

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?

5

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

6 
λ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.

7

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)

2 Likes
8

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.

9

Thanks Chris!