diff --git a/docs/Project.toml b/docs/Project.toml index 697757af90..8164c6976b 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -11,6 +11,8 @@ Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba" OptimizationOptimJL = "36348300-93cb-4f02-beb5-3c3902f8871e" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" +StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0" +StructuralIdentifiability = "220ca800-aa68-49bb-acd8-6037fa93a544" SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b" Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" diff --git a/docs/src/basics/Validation.md b/docs/src/basics/Validation.md index bf5a25782c..f86d3a4e4b 100644 --- a/docs/src/basics/Validation.md +++ b/docs/src/basics/Validation.md @@ -6,7 +6,7 @@ ModelingToolkit.jl provides extensive functionality for model validation and uni Units may assigned with the following syntax. -```julia +```@example validation using ModelingToolkit, Unitful @variables t [unit = u"s"] x(t) [unit = u"m"] g(t) w(t) [unit = "Hz"] @@ -45,21 +45,25 @@ ModelingToolkit.get_unit Example usage below. Note that `ModelingToolkit` does not force unit conversions to preferred units in the event of nonstandard combinations -- it merely checks that the equations are consistent. -```julia +```@example validation using ModelingToolkit, Unitful @parameters τ [unit = u"ms"] @variables t [unit = u"ms"] E(t) [unit = u"kJ"] P(t) [unit = u"MW"] D = Differential(t) eqs = eqs = [D(E) ~ P - E/τ, 0 ~ P ] -ModelingToolkit.validate(eqs) #Returns true -ModelingToolkit.validate(eqs[1]) #Returns true -ModelingToolkit.get_unit(eqs[1].rhs) #Returns u"kJ ms^-1" +ModelingToolkit.validate(eqs) +``` +```@example validation +ModelingToolkit.validate(eqs[1]) +``` +```@example validation +ModelingToolkit.get_unit(eqs[1].rhs) ``` An example of an inconsistent system: at present, `ModelingToolkit` requires that the units of all terms in an equation or sum to be equal-valued (`ModelingToolkit.equivalent(u1,u2)`), rather that simply dimensionally consistent. In the future, the validation stage may be upgraded to support the insertion of conversion factors into the equations. -```julia +```@example validation using ModelingToolkit, Unitful @parameters τ [unit = u"ms"] @variables t [unit = u"ms"] E(t) [unit = u"J"] P(t) [unit = u"MW"] diff --git a/docs/src/tutorials/ode_modeling.md b/docs/src/tutorials/ode_modeling.md index 04850add21..13c66cf20e 100644 --- a/docs/src/tutorials/ode_modeling.md +++ b/docs/src/tutorials/ode_modeling.md @@ -9,9 +9,10 @@ illustrate the basic user-facing functionality. A much deeper tutorial with forcing functions and sparse Jacobians is below. But if you want to just see some code and run, here's an example: -```julia +```@example ode using ModelingToolkit + @variables t x(t) # independent and dependent variables @parameters τ # parameters @constants h = 1 # constants have an assigned value @@ -21,15 +22,14 @@ D = Differential(t) # define an operator for the differentiation w.r.t. time @named fol = ODESystem([ D(x) ~ (h - x)/τ]) using DifferentialEquations: solve -using Plots: plot prob = ODEProblem(fol, [x => 0.0], (0.0,10.0), [τ => 3.0]) # parameter `τ` can be assigned a value, but constant `h` cannot sol = solve(prob) + +using Plots plot(sol) ``` - -![Simulation result of first-order lag element, with right-hand side](https://user-images.githubusercontent.com/13935112/111958369-703f2200-8aed-11eb-8bb4-0abe9652e850.png) Now let's start digging into MTK! ## Your very first ODE @@ -46,7 +46,7 @@ variable, ``f(t)`` is an external forcing function, and ``\tau`` is a parameter. In MTK, this system can be modelled as follows. For simplicity, we first set the forcing function to a time-independent value. -```julia +```@example ode2 using ModelingToolkit @variables t x(t) # independent and dependent variables @@ -56,11 +56,6 @@ D = Differential(t) # define an operator for the differentiation w.r.t. time # your first ODE, consisting of a single equation, indicated by ~ @named fol_model = ODESystem(D(x) ~ (h - x)/τ) - # Model fol_model with 1 equations - # States (1): - # x(t) - # Parameters (1): - # τ ``` Note that equations in MTK use the tilde character (`~`) as equality sign. @@ -69,7 +64,7 @@ matches the name in the REPL. If omitted, you can directly set the `name` keywor After construction of the ODE, you can solve it using [DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/): -```julia +```@example ode2 using DifferentialEquations using Plots @@ -77,8 +72,6 @@ prob = ODEProblem(fol_model, [x => 0.0], (0.0,10.0), [τ => 3.0]) plot(solve(prob)) ``` -![Simulation result of first-order lag element](https://user-images.githubusercontent.com/13935112/111958369-703f2200-8aed-11eb-8bb4-0abe9652e850.png) - The initial state and the parameter values are specified using a mapping from the actual symbolic elements to their values, represented as an array of `Pair`s, which are constructed using the `=>` operator. @@ -88,16 +81,10 @@ of `Pair`s, which are constructed using the `=>` operator. You could separate the calculation of the right-hand side, by introducing an intermediate variable `RHS`: -```julia +```@example ode2 @variables RHS(t) @named fol_separate = ODESystem([ RHS ~ (h - x)/τ, D(x) ~ RHS ]) - # Model fol_separate with 2 equations - # States (2): - # x(t) - # RHS(t) - # Parameters (1): - # τ ``` To directly solve this system, you would have to create a Differential-Algebraic @@ -106,15 +93,13 @@ additional algebraic equation now. However, this DAE system can obviously be transformed into the single ODE we used in the first example above. MTK achieves this by structural simplification: -```julia +```@example ode2 fol_simplified = structural_simplify(fol_separate) - equations(fol_simplified) - # 1-element Array{Equation,1}: - # Differential(t)(x(t)) ~ (τ^-1)*(h - x(t)) +``` +```@example ode2 equations(fol_simplified) == equations(fol_model) - # true ``` You can extract the equations from a system using `equations` (and, in the same @@ -128,7 +113,7 @@ in a simulation result. The intermediate variable `RHS` therefore can be plotted along with the state variable. Note that this has to be requested explicitly, through: -```julia +```@example ode2 prob = ODEProblem(fol_simplified, [x => 0.0], (0.0,10.0), [τ => 3.0]) sol = solve(prob) plot(sol, vars=[x, RHS]) @@ -140,8 +125,6 @@ when using `parameters` (e.g., solution of linear equations by dividing out the constant's value, which cannot be done for parameters, since they may be zero). -![Simulation result of first-order lag element, with right-hand side](https://user-images.githubusercontent.com/13935112/111958403-7e8d3e00-8aed-11eb-9d18-08b5180a59f9.png) - Note that the indexing of the solution similarly works via the names, and so `sol[x]` gives the time-series for `x`, `sol[x,2:10]` gives the 2nd through 10th values of `x` matching `sol.t`, etc. Note that this works even for variables @@ -152,7 +135,7 @@ which have been eliminated, and thus `sol[RHS]` retrieves the values of `RHS`. What if the forcing function (the “external input”) ``f(t)`` is not constant? Obviously, one could use an explicit, symbolic function of time: -```julia +```@example ode2 @variables f(t) @named fol_variable_f = ODESystem([f ~ sin(t), D(x) ~ (f - x)/τ]) ``` @@ -166,7 +149,7 @@ So, you could, for example, interpolate given the time-series using we illustrate this option by a simple lookup ("zero-order hold") of a vector of random values: -```julia +```@example ode2 value_vector = randn(10) f_fun(t) = t >= 10 ? value_vector[end] : value_vector[Int(floor(t))+1] @register_symbolic f_fun(t) @@ -178,15 +161,13 @@ sol = solve(prob) plot(sol, vars=[x,f]) ``` -![Simulation result of first-order lag element, step-wise forcing function](https://user-images.githubusercontent.com/13935112/111958424-83ea8880-8aed-11eb-8f42-489f4b44c3bc.png) - ## Building component-based, hierarchical models Working with simple one-equation systems is already fun, but composing more complex systems from simple ones is even more fun. Best practice for such a “modeling framework” could be to use factory functions for model components: -```julia +```@example ode2 function fol_factory(separate=false;name) @parameters τ @variables t x(t) f(t) RHS(t) @@ -202,7 +183,7 @@ end Such a factory can then be used to instantiate the same component multiple times, but allows for customization: -```julia +```@example ode2 @named fol_1 = fol_factory() @named fol_2 = fol_factory(true) # has observable RHS ``` @@ -212,21 +193,11 @@ Now, these two components can be used as subsystems of a parent system, i.e. one level higher in the model hierarchy. The connections between the components again are just algebraic relations: -```julia +```@example ode2 connections = [ fol_1.f ~ 1.5, fol_2.f ~ fol_1.x ] connected = compose(ODESystem(connections,name=:connected), fol_1, fol_2) - # Model connected with 5 equations - # States (5): - # fol_1₊f(t) - # fol_2₊f(t) - # fol_1₊x(t) - # fol_2₊x(t) - # fol_2₊RHS(t) - # Parameters (2): - # fol_1₊τ - # fol_2₊τ ``` All equations, variables, and parameters are collected, but the structure of the @@ -236,26 +207,12 @@ hierarchical model is still preserved. This means you can still get information variables and connection equations from the system using structural simplification: -```julia +```@example ode2 connected_simp = structural_simplify(connected) - # Model connected with 2 equations - # States (2): - # fol_1₊x(t) - # fol_2₊x(t) - # Parameters (2): - # fol_1₊τ - # fol_2₊τ - # Incidence matrix: - # [1, 1] = × - # [2, 1] = × - # [2, 2] = × - # [1, 3] = × - # [2, 4] = × +``` +```@example ode2 full_equations(connected_simp) - # 2-element Array{Equation,1}: - # Differential(t)(fol_1₊x(t)) ~ (fol_1₊τ^-1)*(1.5 - fol_1₊x(t)) - # Differential(t)(fol_2₊x(t)) ~ (fol_2₊τ^-1)*(fol_1₊x(t) - fol_2₊x(t)) ``` As expected, only the two state-derivative equations remain, as if you had manually eliminated as many variables as possible from the equations. @@ -264,7 +221,7 @@ As mentioned above, the hierarchical structure is preserved. So, the initial state and the parameter values can be specified accordingly when building the `ODEProblem`: -```julia +```@example ode2 u0 = [ fol_1.x => -0.5, fol_2.x => 1.0 ] @@ -275,8 +232,6 @@ prob = ODEProblem(connected_simp, u0, (0.0,10.0), p) plot(solve(prob)) ``` -![Simulation of connected system (two first-order lag elements in series)](https://user-images.githubusercontent.com/13935112/111958439-877e0f80-8aed-11eb-9074-9d35458459a4.png) - More on this topic may be found in [Composing Models and Building Reusable Components](@ref acausal). ## Defaults @@ -284,7 +239,7 @@ More on this topic may be found in [Composing Models and Building Reusable Compo It is often a good idea to specify reasonable values for the initial state and the parameters of a model component. Then, these do not have to be explicitly specified when constructing the `ODEProblem`. -```julia +```@example ode2 function unitstep_fol_factory(;name) @parameters τ @variables t x(t) @@ -311,21 +266,25 @@ By default, analytical derivatives and sparse matrices, e.g. for the Jacobian, t matrix of first partial derivatives, are not used. Let's benchmark this (`prob` still is the problem using the `connected_simp` system above): -```julia +```@example ode2 using BenchmarkTools - -@btime solve($prob, Rodas4()); - # 251.300 μs (873 allocations: 31.18 KiB) +@btime solve(prob, Rodas4()); +nothing # hide ``` Now have MTK provide sparse, analytical derivatives to the solver. This has to be specified during the construction of the `ODEProblem`: -```julia -prob_an = ODEProblem(connected_simp, u0, (0.0,10.0), p; jac=true, sparse=true) +```@example ode2 +prob_an = ODEProblem(connected_simp, u0, (0.0,10.0), p; jac=true) +@btime solve($prob_an, Rodas4()); +nothing # hide +``` +```@example ode2 +prob_an = ODEProblem(connected_simp, u0, (0.0,10.0), p; jac=true, sparse=true) @btime solve($prob_an, Rodas4()); - # 142.899 μs (1297 allocations: 83.96 KiB) +nothing # hide ``` The speedup is significant. For this small dense model (3 of 4 entries are diff --git a/docs/src/tutorials/parameter_identifiability.md b/docs/src/tutorials/parameter_identifiability.md index 24c3609f08..26ebc1e7d5 100644 --- a/docs/src/tutorials/parameter_identifiability.md +++ b/docs/src/tutorials/parameter_identifiability.md @@ -4,12 +4,6 @@ Ordinary differential equations are commonly used for modeling real-world proces We will start by illustrating **local identifiability** in which a parameter is known up to _finitely many values_, and then proceed to determining **global identifiability**, that is, which parameters can be identified _uniquely_. -To install `StructuralIdentifiability.jl`, simply run -```julia -using Pkg -Pkg.add("StructuralIdentifiability") -``` - The package has a standalone data structure for ordinary differential equations, but is also compatible with `ODESystem` type from `ModelingToolkit.jl`. ## Local Identifiability @@ -31,7 +25,7 @@ This model describes the biohydrogenation[^1] process[^2] with unknown initial c To define the ode system in Julia, we use `ModelingToolkit.jl`. We first define the parameters, variables, differential equations and the output equations. -```julia +```@example SI using StructuralIdentifiability, ModelingToolkit # define parameters and variables @@ -52,34 +46,20 @@ measured_quantities = [y1 ~ x4, y2 ~ x5] # define the system de = ODESystem(eqs, t, name=:Biohydrogenation) - ``` After that, we are ready to check the system for local identifiability: -```julia +```@example SI # query local identifiability # we pass the ode-system local_id_all = assess_local_identifiability(de, measured_quantities=measured_quantities, p=0.99) - # [ Info: Preproccessing `ModelingToolkit.ODESystem` object - # 6-element Vector{Bool}: - # 1 - # 1 - # 1 - # 1 - # 1 - # 1 ``` We can see that all states (except $x_7$) and all parameters are locally identifiable with probability 0.99. Let's try to check specific parameters and their combinations -```julia +```@example SI to_check = [k5, k7, k10/k9, k5+k6] local_id_some = assess_local_identifiability(de, measured_quantities=measured_quantities, funcs_to_check=to_check, p=0.99) - # 4-element Vector{Bool}: - # 1 - # 1 - # 1 - # 1 ``` Notice that in this case, everything (except the state variable $x_7$) is locally identifiable, including combinations such as $k_{10}/k_9, k_5+k_6$ @@ -106,7 +86,7 @@ Global identifiability needs information about local identifiability first, but __Note__: as of writing this tutorial, UTF-symbols such as Greek characters are not supported by one of the project's dependencies, see [this issue](https://github.com/SciML/StructuralIdentifiability.jl/issues/43). -```julia +```@example SI2 using StructuralIdentifiability, ModelingToolkit @parameters b c a beta g delta sigma @variables t x1(t) x2(t) x3(t) x4(t) y(t) y2(t) @@ -124,25 +104,16 @@ measured_quantities = [y~x1+x2, y2~x2] ode = ODESystem(eqs, t, name=:GoodwinOsc) -@time global_id = assess_identifiability(ode, measured_quantities=measured_quantities) - # 30.672594 seconds (100.97 M allocations: 6.219 GiB, 3.15% gc time, 0.01% compilation time) - # Dict{Num, Symbol} with 7 entries: - # a => :globally - # b => :globally - # beta => :globally - # c => :globally - # sigma => :globally - # g => :nonidentifiable - # delta => :globally +global_id = assess_identifiability(ode, measured_quantities=measured_quantities) ``` We can see that only parameters `a, g` are unidentifiable, and everything else can be uniquely recovered. Let us consider the same system but with two inputs, and we will find out identifiability with probability `0.9` for parameters `c` and `b`: -```julia +```@example SI3 using StructuralIdentifiability, ModelingToolkit @parameters b c a beta g delta sigma -@variables t x1(t) x2(t) x3(t) x4(t) y(t) u1(t) [input=true] u2(t) [input=true] +@variables t x1(t) x2(t) x3(t) x4(t) y(t) y2(t) u1(t) [input=true] u2(t) [input=true] D = Differential(t) eqs = [ @@ -159,9 +130,6 @@ to_check = [b, c] ode = ODESystem(eqs, t, name=:GoodwinOsc) global_id = assess_identifiability(ode, measured_quantities=measured_quantities, funcs_to_check=to_check, p=0.9) - # Dict{Num, Symbol} with 2 entries: - # b => :globally - # c => :globally ``` Both parameters `b, c` are globally identifiable with probability `0.9` in this case. diff --git a/docs/src/tutorials/stochastic_diffeq.md b/docs/src/tutorials/stochastic_diffeq.md index 7db7d6e678..3d9d622a35 100644 --- a/docs/src/tutorials/stochastic_diffeq.md +++ b/docs/src/tutorials/stochastic_diffeq.md @@ -7,7 +7,7 @@ which has a deterministic (drift) component and a stochastic (diffusion) component. Let's take the Lorenz equation from the first tutorial and extend it to have multiplicative noise. -```julia +```@example SDE using ModelingToolkit, StochasticDiffEq # Define some variables