Merge branch 'master' into myb/state_machine

Author: YingboMa

NEWS.md

@@ -20,7 +20,7 @@
       + `t_nounits` and `D_nounits` are unitless.
   - `ODAEProblem` is deprecated in favor of `ODEProblem`.
   - Specifying the independent variable for an `ODESystem` is now mandatory. The `ODESystem(eqs)`
-    constructor is removed.
+    constructor is removed. Use `ODESystem(eqs,t)` instead.
   - Systems must be marked as `complete` before creating `*Function`/`*FunctionExpr`/`*Problem`/
     `*ProblemExpr`. Typically this involved using `@mtkbuild` to create the system or calling
     `structural_simplify` on an existing system.
@@ -45,3 +45,7 @@
     equations. For example, `[p[1] => 1.0, p[2] => 2.0]` is no longer allowed in default equations, use
     `[p => [1.0, 2.0]]` instead. Also, array equations like for `@variables u[1:2]` have `D(u) ~ A*u` as an
     array equation. If the scalarized version is desired, use `scalarize(u)`.
+  - Parameter dependencies are now supported. They can be specified using the syntax
+    `(single_parameter => expression_involving_other_parameters)` and a `Vector` of these can be passed to
+    the `parameter_dependencies` keyword argument of `ODESystem`, `SDESystem` and `JumpSystem`. The dependent
+    parameters are updated whenever other parameters are modified, e.g. in callbacks.

Project.toml

@@ -1,7 +1,7 @@
 name = "ModelingToolkit"
 uuid = "961ee093-0014-501f-94e3-6117800e7a78"
 authors = ["Yingbo Ma <[email protected]>", "Chris Rackauckas <[email protected]> and contributors"]
-version = "8.76.0"
+version = "9.0.1"
 
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
@@ -69,7 +69,7 @@ Compat = "3.42, 4"
 ConstructionBase = "1"
 DataStructures = "0.17, 0.18"
 DiffEqBase = "6.103.0"
-DiffEqCallbacks = "2.16"
+DiffEqCallbacks = "2.16, 3"
 DiffRules = "0.1, 1.0"
 Distributed = "1"
 Distributions = "0.23, 0.24, 0.25"
@@ -97,16 +97,16 @@ RecursiveArrayTools = "2.3, 3"
 Reexport = "0.2, 1"
 RuntimeGeneratedFunctions = "0.5.9"
 SciMLBase = "2.0.1"
+SciMLStructures = "1.0"
 Serialization = "1"
 Setfield = "0.7, 0.8, 1"
 SimpleNonlinearSolve = "0.1.0, 1"
-SciMLStructures = "1.0"
 SparseArrays = "1"
 SpecialFunctions = "0.7, 0.8, 0.9, 0.10, 1.0, 2"
 StaticArrays = "0.10, 0.11, 0.12, 1.0"
 SymbolicIndexingInterface = "0.3.1"
 SymbolicUtils = "1.0"
-Symbolics = "5.7"
+Symbolics = "5.20.1"
 URIs = "1"
 UnPack = "0.1, 1.0"
 Unitful = "1.1"

docs/Project.toml

@@ -25,7 +25,7 @@ BifurcationKit = "0.3"
 DifferentialEquations = "7.6"
 Distributions = "0.25"
 Documenter = "1"
-ModelingToolkit = "8.33"
+ModelingToolkit = "8.33, 9"
 ModelingToolkitDesigner = "1"
 NonlinearSolve = "0.3, 1, 2, 3"
 Optim = "1.7"

docs/src/basics/Composition.md

@@ -16,22 +16,20 @@ of `decay2` is the value of the unknown variable `x`.
 
 ```@example composition
 using ModelingToolkit
+using ModelingToolkit: t_nounits as t, D_nounits as D
 
 function decay(; name)
-    @parameters t a
+    @parameters a
     @variables x(t) f(t)
-    D = Differential(t)
     ODESystem([
             D(x) ~ -a * x + f
-        ];
+        ], t;
         name = name)
 end
 
 @named decay1 = decay()
 @named decay2 = decay()
 
-@parameters t
-D = Differential(t)
 connected = compose(
     ODESystem([decay2.f ~ decay1.x
                D(decay1.f) ~ 0], t; name = :connected), decay1, decay2)
@@ -137,7 +135,7 @@ sys.y = u * 1.1
 In a hierarchical system, variables of the subsystem get namespaced by the name of the system they are in. This prevents naming clashes, but also enforces that every unknown and parameter is local to the subsystem it is used in. In some cases it might be desirable to have variables and parameters that are shared between subsystems, or even global. This can be accomplished as follows.
 
 ```julia
-@parameters t a b c d e f
+@parameters a b c d e f
 
 # a is a local variable
 b = ParentScope(b) # b is a variable that belongs to one level up in the hierarchy
@@ -197,19 +195,16 @@ higher level system:
 
 ```@example compose
 using ModelingToolkit, OrdinaryDiffEq, Plots
+using ModelingToolkit: t_nounits as t, D_nounits as D
 
 ## Library code
-
-@parameters t
-D = Differential(t)
-
 @variables S(t), I(t), R(t)
 N = S + I + R
 @parameters β, γ
 
-@named seqn = ODESystem([D(S) ~ -β * S * I / N])
-@named ieqn = ODESystem([D(I) ~ β * S * I / N - γ * I])
-@named reqn = ODESystem([D(R) ~ γ * I])
+@named seqn = ODESystem([D(S) ~ -β * S * I / N], t)
+@named ieqn = ODESystem([D(I) ~ β * S * I / N - γ * I], t)
+@named reqn = ODESystem([D(R) ~ γ * I], t)
 
 sir = compose(
     ODESystem(

docs/src/basics/Events.md

@@ -65,9 +65,10 @@ friction
 
 ```@example events
 using ModelingToolkit, OrdinaryDiffEq, Plots
+using ModelingToolkit: t_nounits as t, D_nounits as D
+
 function UnitMassWithFriction(k; name)
-    @variables t x(t)=0 v(t)=0
-    D = Differential(t)
+    @variables x(t)=0 v(t)=0
     eqs = [D(x) ~ v
            D(v) ~ sin(t) - k * sign(v)]
     ODESystem(eqs, t; continuous_events = [v ~ 0], name) # when v = 0 there is a discontinuity
@@ -87,8 +88,7 @@ an `affect!` on the state. We can model the same system using ModelingToolkit
 like this
 
 ```@example events
-@variables t x(t)=1 v(t)=0
-D = Differential(t)
+@variables x(t)=1 v(t)=0
 
 root_eqs = [x ~ 0]  # the event happens at the ground x(t) = 0
 affect = [v ~ -v] # the effect is that the velocity changes sign
@@ -108,8 +108,7 @@ plot(sol)
 Multiple events? No problem! This example models a bouncing ball in 2D that is enclosed by two walls at $y = \pm 1.5$.
 
 ```@example events
-@variables t x(t)=1 y(t)=0 vx(t)=0 vy(t)=2
-D = Differential(t)
+@variables x(t)=1 y(t)=0 vx(t)=0 vy(t)=2
 
 continuous_events = [[x ~ 0] => [vx ~ -vx]
                      [y ~ -1.5, y ~ 1.5] => [vy ~ -vy]]
@@ -229,7 +228,7 @@ Suppose we have a population of `N(t)` cells that can grow and die, and at time
 
 ```@example events
 @parameters M tinject α
-@variables t N(t)
+@variables N(t)
 Dₜ = Differential(t)
 eqs = [Dₜ(N) ~ α - N]
 

docs/src/basics/FAQ.md

@@ -102,9 +102,10 @@ end
 This error can come up after running `structural_simplify` on a system that generates dummy derivatives (i.e. variables with `ˍt`).  For example, here even though all the variables are defined with initial values, the `ODEProblem` generation will throw an error that defaults are missing from the variable map.
 
 ```
-@variables t
+using ModelingToolkit
+using ModelingToolkit: t_nounits as t, D_nounits as D
+
 sts = @variables x1(t)=0.0 x2(t)=0.0 x3(t)=0.0 x4(t)=0.0
-D = Differential(t)
 eqs = [x1 + x2 + 1 ~ 0
        x1 + x2 + x3 + 2 ~ 0
        x1 + D(x3) + x4 + 3 ~ 0
@@ -137,9 +138,9 @@ container type. For example:
 
 ```
 using ModelingToolkit, StaticArrays
-@variables t
+using ModelingToolkit: t_nounits as t, D_nounits as D
+
 sts = @variables x1(t)=0.0
-D = Differential(t)
 eqs = [D(x1) ~ 1.1 * x1]
 @mtkbuild sys = ODESystem(eqs, t)
 prob = ODEProblem{false}(sys, [], (0,1); u0_constructor = x->SVector(x...))

docs/src/basics/Linearization.md

@@ -21,9 +21,9 @@ The `linearize` function expects the user to specify the inputs ``u`` and the ou
 
 ```@example LINEARIZE
 using ModelingToolkit
-@variables t x(t)=0 y(t)=0 u(t)=0 r(t)=0
+using ModelingToolkit: t_nounits as t, D_nounits as D
+@variables x(t)=0 y(t)=0 u(t)=0 r(t)=0
 @parameters kp = 1
-D = Differential(t)
 
 eqs = [u ~ kp * (r - y) # P controller
        D(x) ~ -x + u    # First-order plant

docs/src/basics/Validation.md

@@ -83,7 +83,7 @@ second argument.
 ```@example validation2
 using ModelingToolkit, Unitful
 # Composite type parameter in registered function
-@parameters t
+@variables t
 D = Differential(t)
 struct NewType
     f::Any

docs/src/basics/Variable_metadata.md

@@ -10,14 +10,14 @@ Descriptive strings can be attached to variables using the `[description = "desc
 
 ```@example metadata
 using ModelingToolkit
+using ModelingToolkit: t_nounits as t, D_nounits as D
 @variables u [description = "This is my input"]
 getdescription(u)
 ```
 
 When variables with descriptions are present in systems, they will be printed when the system is shown in the terminal:
 
 ```@example metadata
-@parameters t
 @variables u(t) [description = "A short description of u"]
 @parameters p [description = "A description of p"]
 @named sys = ODESystem([u ~ p], t)
@@ -50,8 +50,9 @@ current in a resistor. These variables sum up to zero in connections.
 
 ```@example connect
 using ModelingToolkit
+using ModelingToolkit: t_nounits as t, D_nounits as D
 
-@variables t, i(t) [connect = Flow]
+@variables i(t) [connect = Flow]
 @variables k(t) [connect = Stream]
 ```
 
@@ -61,6 +62,8 @@ Designate a variable as either an input or an output using the following
 
 ```@example metadata
 using ModelingToolkit
+using ModelingToolkit: t_nounits as t, D_nounits as D
+
 @variables u [input = true]
 isinput(u)
 ```
@@ -136,8 +139,6 @@ For systems that contain parameters with metadata like described above, have som
 In the example below, we define a system with tunable parameters and extract bounds vectors
 
 ```@example metadata
-@parameters t
-Dₜ = Differential(t)
 @variables x(t)=0 u(t)=0 [input = true] y(t)=0 [output = true]
 @parameters T [tunable = true, bounds = (0, Inf)]
 @parameters k [tunable = true, bounds = (0, Inf)]

docs/src/examples/higher_order.md

@@ -13,16 +13,16 @@ We utilize the derivative operator twice here to define the second order:
 
 ```@example orderlowering
 using ModelingToolkit, OrdinaryDiffEq
+using ModelingToolkit: t_nounits as t, D_nounits as D
 
 @parameters σ ρ β
-@variables t x(t) y(t) z(t)
-D = Differential(t)
+@variables x(t) y(t) z(t)
 
 eqs = [D(D(x)) ~ σ * (y - x),
     D(y) ~ x * (ρ - z) - y,
     D(z) ~ x * y - β * z]
 
-@named sys = ODESystem(eqs)
+@named sys = ODESystem(eqs, t)
 ```
 
 Note that we could've used an alternative syntax for 2nd order, i.e.

docs/src/examples/modelingtoolkitize_index_reduction.md

@@ -31,7 +31,7 @@ pendulum_prob = ODEProblem(pendulum_fun!, u0, tspan, p)
 traced_sys = modelingtoolkitize(pendulum_prob)
 pendulum_sys = structural_simplify(dae_index_lowering(traced_sys))
 prob = ODEProblem(pendulum_sys, [], tspan)
-sol = solve(prob, Tsit5(), abstol = 1e-8, reltol = 1e-8)
+sol = solve(prob, Rodas5P(), abstol = 1e-8, reltol = 1e-8)
 plot(sol, idxs = unknowns(traced_sys))
 ```
 
@@ -71,7 +71,7 @@ u0 = [1.0, 0, 0, 0, 0];
 p = [9.8, 1];
 tspan = (0, 10.0);
 pendulum_prob = ODEProblem(pendulum_fun!, u0, tspan, p)
-solve(pendulum_prob, Rodas4())
+solve(pendulum_prob, Rodas5P())
 ```
 
 However, one will quickly be greeted with the unfortunate message:
@@ -151,7 +151,7 @@ numerical solver. Let's try that out:
 traced_sys = modelingtoolkitize(pendulum_prob)
 pendulum_sys = structural_simplify(dae_index_lowering(traced_sys))
 prob = ODEProblem(pendulum_sys, Pair[], tspan)
-sol = solve(prob, Rodas4())
+sol = solve(prob, Rodas5P())
 
 using Plots
 plot(sol, idxs = unknowns(traced_sys))

docs/src/examples/perturbation.md

@@ -45,8 +45,9 @@ As the first ODE example, we have chosen a simple and well-behaved problem, whic
 with the initial conditions $x(0) = 0$, and $\dot{x}(0) = 1$. Note that for $\epsilon = 0$, this equation transforms back to the standard one. Let's start with defining the variables
 
 ```julia
+using ModelingToolkit: t_nounits as t, D_nounits as D
 n = 3
-@variables ϵ t y[1:n](t) ∂∂y[1:n](t)
+@variables ϵ y[1:n](t) ∂∂y[1:n](t)
 ```
 
 Next, we define $x$.
@@ -82,7 +83,6 @@ vals = solve_coef(eqs, ∂∂y)
 Our system of ODEs is forming. Now is the time to convert `∂∂`s to the correct **Symbolics.jl** form by substitution:
 
 ```julia
-D = Differential(t)
 subs = Dict(∂∂y[i] => D(D(y[i])) for i in eachindex(y))
 eqs = [substitute(first(v), subs) ~ substitute(last(v), subs) for v in vals]
 ```
@@ -147,7 +147,6 @@ vals = solve_coef(eqs, ∂∂y)
 Next, we need to replace `∂`s and `∂∂`s with their **Symbolics.jl** counterparts:
 
 ```julia
-D = Differential(t)
 subs1 = Dict(∂y[i] => D(y[i]) for i in eachindex(y))
 subs2 = Dict(∂∂y[i] => D(D(y[i])) for i in eachindex(y))
 subs = subs1 ∪ subs2

docs/src/examples/spring_mass.md

@@ -6,11 +6,9 @@ In this tutorial, we will build a simple component-based model of a spring-mass
 
 ```@example component
 using ModelingToolkit, Plots, DifferentialEquations, LinearAlgebra
+using ModelingToolkit: t_nounits as t, D_nounits as D
 using Symbolics: scalarize
 
-@variables t
-D = Differential(t)
-
 function Mass(; name, m = 1.0, xy = [0.0, 0.0], u = [0.0, 0.0])
     ps = @parameters m = m
     sts = @variables pos(t)[1:2]=xy v(t)[1:2]=u

docs/src/examples/tearing_parallelism.md

@@ -11,10 +11,9 @@ electrical circuits:
 
 ```@example tearing
 using ModelingToolkit, OrdinaryDiffEq
+using ModelingToolkit: t_nounits as t, D_nounits as D
 
 # Basic electric components
-@variables t
-const D = Differential(t)
 @connector function Pin(; name)
     @variables v(t)=1.0 i(t)=1.0 [connect = Flow]
     ODESystem(Equation[], t, [v, i], [], name = name)

docs/src/tutorials/acausal_components.md

@@ -21,8 +21,8 @@ equalities before solving. Let's see this in action.
 
 ```@example acausal
 using ModelingToolkit, Plots, DifferentialEquations
+using ModelingToolkit: t_nounits as t, D_nounits as D
 
-@variables t
 @connector Pin begin
     v(t)
     i(t), [connect = Flow]
@@ -63,8 +63,6 @@ end
     end
 end
 
-D = Differential(t)
-
 @mtkmodel Capacitor begin
     @extend OnePort()
     @parameters begin
@@ -213,8 +211,6 @@ equations and unknowns and extend them with a new equation. Note that `v`, `i` a
 Using our knowledge of circuits, we similarly construct the `Capacitor`:
 
 ```@example acausal
-D = Differential(t)
-
 @mtkmodel Capacitor begin
     @extend OnePort()
     @parameters begin