sir_modia.jmd

# Ordinary differential equation model using ModiaSim
Simon Frost (@sdwfrost), 2021-04-15

## Introduction

The classical ODE version of the SIR model is:

- Deterministic
- Continuous in time
- Continuous in state

This tutorial uses [ModiaSim](https://modiasim.github.io/docs/) to specify the model. ModiaSim is an acausal Modelica-like modeling framework composed of a number of Julia packages that allows specification of models as differential-algebraic equations (DAE). ModiaSim expects parameter values to have units; here, we assume that the time variable is in days.

## Libraries

```julia
using TinyModia
using Unitful
using DifferentialEquations
using ModiaPlot
using BenchmarkTools
```

## Transitions

The following function call defines the equations of the model, including the definition of the total population size, as well as the gradients of the variables (using `der`). At this point, the model won't run, as it does not have the parameter values or initial conditions.

```julia
SIR = Model(
    equations = :[
        N = S + I + R
        der(S) = -β*c*I/N*S
        der(I) = β*c*I/N*S - γ*I
        der(R) = γ*I
    ]
);
```

## Time domain

We set the stopping time, `tmax`, for simulations, as well as interval to output results.

```julia
δt = 0.1u"d"
tmax = 40.0u"d"
```

## Initial conditions

Initial conditions are described as a `Map` that defines variable names (here, `S`, `I`, and `R`) as `Var` with an initial value.

```julia
u0 = Map(S = Var(init=990.0),
         I = Var(init=10.0),
         R = Var(init=0.0));
```

## Parameter values

Parameter values are also defined using a `Map`.

```julia
p = Map(β = 0.05,
        c = 10.0u"1/d",
        γ = 0.25u"1/d");
```

## Running the model

To define a full model, we need to merge the equations, initial conditions, and parameter values using the merge operator, `|`.

```julia
model = SIR | u0 | p;
```

The model is then instantiated using a macro.

```julia
sir = @instantiateModel(model);
```

Simulation changes the model instance in-place.

```julia
simulate!(sir,Tsit5(),stopTime=tmax,interval=δt);
```

## Plotting

We can now plot the results using `ModiaPlot`.

```julia
ModiaPlot.plot(sir,[("S","I","R")]);
```

Alternatively, we can extract the result from `sir` and plot using `StatsPlots`. Modia returns the initial condition alongside the simulation, and so we omit the first row of the `DataFrame`, and use the variable names in the `sir` object to name the columns.

```julia
using StatsPlots
using DataFrames
df = DataFrame(sir.result)
dfnames = collect(keys(sir.variables))
dfnames[1] = "t"
rename!(df,dfnames)
df = df[2:end,dfnames]
first(df,6)
```

```julia
@df df StatsPlots.plot(:t,
    [:S :I :R],
    xlabel="Time",
    ylabel="Number")
```

## Benchmarking

```julia
@benchmark simulate!(sir,Tsit5(),stopTime=tmax,interval=δt)
```