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,λ)?

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

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?

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.