StackOverFlow Error in SIR model using ModelingTookit

You can add observed variables, which are just scalings of other state variables, and the system will automatically eliminate them from the solve. You then just do sol[r_scaled] or whatnot and it’ll give you those values in the solution.

Hmm, I’m not sure how to implement that…Do you think you could illustrate this with a simple example? I’ve checked this page for more info but this particular topic doesn’t have any code…

I’m not sure how to add them back after the fact, but you can do this:

using ModelingToolkit, DifferentialEquations
using Plots

@parameters t β γ N
@variables S(t) I(t) R(t) s(t) i(t) r(t) x y
D = Differential(t)

SIR_eqs = [D(s) ~ -β*s*i,
           D(i) ~ β*s*i - γ*i,
           D(r) ~ γ*i,
            0 ~ s*N - S,
            0 ~ i*N - I,
            0 ~ r*N - R,
]

SIR = ODESystem(SIR_eqs, name = :SIR)

# After this, the system has only 3 states [s,i,r]
SIR = structural_simplify(SIR)

u0 = [s => 0.99, i => 0.01, r => 0.00]
pars = [β => 0.20, γ => 0.15, N=>1e3]

prob = ODEProblem(SIR, u0, (0.0, 200.0), pars)

sol = solve(prob; saveat=1:1:200)

To plot scaled

plot(sol; vars=[s,i,r])

Or to plot unscaled

plot(sol; vars=[S,I,R])
2 Likes

@contradict Thanks, this helps a lot! Let me see if I understand correctly: by saying

SIR = structural_simplify(SIR)

the variables s,i,r are the variables which are actually being integrated, whereas as the variables S,I,R are the “observables” which are calculated from the solution for s(t),i(t),r(t)?

Yes, in this case that is what happens. structural_simplify is kind of a magic box that applies several tools to try to make a simpler system to solve. As such, it isn’t very controllable so you sometimes have to play with the input equations if you want a particular shape. The fun starts here if you want to see what it is doing.

1 Like

Awesome, this looks to be exactly what I originally wanted!

Another question: How does structural_simplify “know” that s,i,r are the observables? Is it the 0 on the LHS of the equations? Or the fact that the s,i,r equations don’t involve derivatives? Or did we just happen to get lucky in this case and structural_simplify somehow inferred what the observables would be?

…I’ll check out the source code in more detail (though based on a quick glance, it seems a bit over my head, currently :slight_smile: ).