From c0052f200f5435bb765151e39fa75e7469c209d8 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Fri, 23 May 2025 15:26:01 +0530 Subject: [PATCH 01/49] docs: update `basics` docs to new syntax --- docs/src/basics/AbstractSystem.md | 2 +- docs/src/basics/Composition.md | 10 +++++----- docs/src/basics/Debugging.md | 4 ++-- docs/src/basics/Events.md | 29 +++++++++++++++-------------- docs/src/basics/FAQ.md | 8 ++++---- docs/src/basics/InputOutput.md | 2 +- docs/src/basics/Linearization.md | 4 ++-- docs/src/basics/MTKLanguage.md | 18 +++++++++--------- docs/src/basics/Precompilation.md | 2 +- docs/src/basics/Validation.md | 2 +- 10 files changed, 41 insertions(+), 40 deletions(-) diff --git a/docs/src/basics/AbstractSystem.md b/docs/src/basics/AbstractSystem.md index d1707f822f..62f842b2b2 100644 --- a/docs/src/basics/AbstractSystem.md +++ b/docs/src/basics/AbstractSystem.md @@ -152,7 +152,7 @@ a lower level in the system. ## Namespacing By default, unsimplified systems will namespace variables accessed via `getproperty`. -Systems created via `@mtkbuild`, or ones passed through `structural_simplify` or +Systems created via `@mtkcompile`, or ones passed through `mtkcompile` or `complete` will not perform this namespacing. However, all of these processes modify the system in a variety of ways. To toggle namespacing without transforming any other property of the system, use `toggle_namespacing`. diff --git a/docs/src/basics/Composition.md b/docs/src/basics/Composition.md index d3de71d696..f6abe97d38 100644 --- a/docs/src/basics/Composition.md +++ b/docs/src/basics/Composition.md @@ -42,7 +42,7 @@ equations(connected) # Differential(t)(decay1₊x(t)) ~ decay1₊f(t) - (decay1₊a*(decay1₊x(t))) # Differential(t)(decay2₊x(t)) ~ decay2₊f(t) - (decay2₊a*(decay2₊x(t))) -simplified_sys = structural_simplify(connected) +simplified_sys = mtkcompile(connected) equations(simplified_sys) ``` @@ -84,7 +84,7 @@ example, let's say there is a variable `x` in `unknowns` and a variable `x` in `subsys`. We can declare that these two variables are the same by specifying their equality: `x ~ subsys.x` in the `eqs` for `sys`. This algebraic relationship can then be simplified by transformations -like `structural_simplify` which will be described later. +like `mtkcompile` which will be described later. ### Numerics with Composed Models @@ -169,7 +169,7 @@ parameters(level3) In many cases, the nicest way to build a model may leave a lot of unnecessary variables. Thus one may want to remove these equations -before numerically solving. The `structural_simplify` function removes +before numerically solving. The `mtkcompile` function removes these trivial equality relationships and trivial singularity equations, i.e. equations which result in `0~0` expressions, in over-specified systems. @@ -227,7 +227,7 @@ values. The user of this model can then solve this model simply by specifying the values at the highest level: ```@example compose -sireqn_simple = structural_simplify(sir) +sireqn_simple = mtkcompile(sir) equations(sireqn_simple) ``` @@ -251,7 +251,7 @@ sol[reqn.R] ## Tearing Problem Construction Some system types (specifically `NonlinearSystem`) can be further -reduced if `structural_simplify` has already been applied to them. This is done +reduced if `mtkcompile` has already been applied to them. This is done by using the alternative problem constructors (`BlockNonlinearProblem`). In these cases, the constructor uses the knowledge of the strongly connected components calculated during the process of simplification diff --git a/docs/src/basics/Debugging.md b/docs/src/basics/Debugging.md index 6e2d471461..5bc509acfd 100644 --- a/docs/src/basics/Debugging.md +++ b/docs/src/basics/Debugging.md @@ -13,7 +13,7 @@ using ModelingToolkit: t_nounits as t, D_nounits as D eqs = [D(u1) ~ -√(u1), D(u2) ~ -√(u2), D(u3) ~ -√(u3)] defaults = [u1 => 1.0, u2 => 2.0, u3 => 3.0] @named sys = ODESystem(eqs, t; defaults) -sys = structural_simplify(sys) +sys = mtkcompile(sys) ``` This problem causes the ODE solver to crash: @@ -38,7 +38,7 @@ We could have figured that out ourselves, but it is not always so obvious for mo Suppose we also want to validate that `u1 + u2 >= 2.0`. We can do this via the assertions functionality. ```@example debug -@mtkbuild sys = ODESystem(eqs, t; defaults, assertions = [(u1 + u2 >= 2.0) => "Oh no!"]) +@mtkcompile sys = ODESystem(eqs, t; defaults, assertions = [(u1 + u2 >= 2.0) => "Oh no!"]) ``` The assertions must be an iterable of pairs, where the first element is the symbolic condition and diff --git a/docs/src/basics/Events.md b/docs/src/basics/Events.md index 2808204d61..5d5df0a377 100644 --- a/docs/src/basics/Events.md +++ b/docs/src/basics/Events.md @@ -59,7 +59,7 @@ For example, consider the following system. ```julia @variables x(t) y(t) @parameters p(t) -@mtkbuild sys = ODESystem([x * y ~ p, D(x) ~ 0], t) +@mtkcompile sys = ODESystem([x * y ~ p, D(x) ~ 0], t) event = [t == 1] => [x ~ Pre(x) + 1] ``` @@ -134,7 +134,7 @@ function UnitMassWithFriction(k; name) D(v) ~ sin(t) - k * sign(v)] ODESystem(eqs, t; continuous_events = [v ~ 0], name) # when v = 0 there is a discontinuity end -@mtkbuild m = UnitMassWithFriction(0.7) +@mtkcompile m = UnitMassWithFriction(0.7) prob = ODEProblem(m, Pair[], (0, 10pi)) sol = solve(prob, Tsit5()) plot(sol) @@ -154,8 +154,9 @@ like this root_eqs = [x ~ 0] # the event happens at the ground x(t) = 0 affect = [v ~ -Pre(v)] # the effect is that the velocity changes sign -@mtkbuild ball = ODESystem([D(x) ~ v - D(v) ~ -9.8], t; continuous_events = root_eqs => affect) # equation => affect +@mtkcompile ball = ODESystem( + [D(x) ~ v + D(v) ~ -9.8], t; continuous_events = root_eqs => affect) # equation => affect tspan = (0.0, 5.0) prob = ODEProblem(ball, Pair[], tspan) @@ -174,7 +175,7 @@ Multiple events? No problem! This example models a bouncing ball in 2D that is e continuous_events = [[x ~ 0] => [vx ~ -Pre(vx)] [y ~ -1.5, y ~ 1.5] => [vy ~ -Pre(vy)]] -@mtkbuild ball = ODESystem( +@mtkcompile ball = ODESystem( [ D(x) ~ vx, D(y) ~ vy, @@ -254,7 +255,7 @@ end reflect = [x ~ 0] => (bb_affect!, [v], [], [], nothing) -@mtkbuild bb_sys = ODESystem(bb_eqs, t, sts, par, +@mtkcompile bb_sys = ODESystem(bb_eqs, t, sts, par, continuous_events = reflect) u0 = [v => 0.0, x => 1.0] @@ -299,7 +300,7 @@ injection = (t == tinject) => [N ~ Pre(N) + M] u0 = [N => 0.0] tspan = (0.0, 20.0) p = [α => 100.0, tinject => 10.0, M => 50] -@mtkbuild osys = ODESystem(eqs, t, [N], [α, M, tinject]; discrete_events = injection) +@mtkcompile osys = ODESystem(eqs, t, [N], [α, M, tinject]; discrete_events = injection) oprob = ODEProblem(osys, u0, tspan, p) sol = solve(oprob, Tsit5(); tstops = 10.0) plot(sol) @@ -318,7 +319,7 @@ to ```@example events injection = ((t == tinject) & (N < 50)) => [N ~ Pre(N) + M] -@mtkbuild osys = ODESystem(eqs, t, [N], [M, tinject, α]; discrete_events = injection) +@mtkcompile osys = ODESystem(eqs, t, [N], [M, tinject, α]; discrete_events = injection) oprob = ODEProblem(osys, u0, tspan, p) sol = solve(oprob, Tsit5(); tstops = 10.0) plot(sol) @@ -345,7 +346,7 @@ killing = ModelingToolkit.SymbolicDiscreteCallback( tspan = (0.0, 30.0) p = [α => 100.0, tinject => 10.0, M => 50, tkill => 20.0] -@mtkbuild osys = ODESystem(eqs, t, [N], [α, M, tinject, tkill]; +@mtkcompile osys = ODESystem(eqs, t, [N], [α, M, tinject, tkill]; discrete_events = [injection, killing]) oprob = ODEProblem(osys, u0, tspan, p) sol = solve(oprob, Tsit5(); tstops = [10.0, 20.0]) @@ -374,7 +375,7 @@ killing = ModelingToolkit.SymbolicDiscreteCallback( [20.0] => [α ~ 0.0], discrete_parameters = α, iv = t) p = [α => 100.0, M => 50] -@mtkbuild osys = ODESystem(eqs, t, [N], [α, M]; +@mtkcompile osys = ODESystem(eqs, t, [N], [α, M]; discrete_events = [injection, killing]) oprob = ODEProblem(osys, u0, tspan, p) sol = solve(oprob, Tsit5()) @@ -414,7 +415,7 @@ example: ev = ModelingToolkit.SymbolicDiscreteCallback( 1.0 => [c ~ Pre(c) + 1], discrete_parameters = c, iv = t) -@mtkbuild sys = ODESystem( +@mtkcompile sys = ODESystem( D(x) ~ c * cos(x), t, [x], [c]; discrete_events = [ev]) prob = ODEProblem(sys, [x => 0.0], (0.0, 2pi), [c => 1.0]) @@ -435,7 +436,7 @@ will be saved. If we repeat the above example with `c` not a `discrete_parameter @variables x(t) @parameters c(t) -@mtkbuild sys = ODESystem( +@mtkcompile sys = ODESystem( D(x) ~ c * cos(x), t, [x], [c]; discrete_events = [1.0 => [c ~ Pre(c) + 1]]) prob = ODEProblem(sys, [x => 0.0], (0.0, 2pi), [c => 1.0]) @@ -538,7 +539,7 @@ to the system. ```@example events @named sys = ODESystem( eqs, t, [temp], params; continuous_events = [furnace_disable, furnace_enable]) -ss = structural_simplify(sys) +ss = mtkcompile(sys) prob = ODEProblem(ss, [temp => 0.0, furnace_on => true], (0.0, 10.0)) sol = solve(prob, Tsit5()) plot(sol) @@ -651,7 +652,7 @@ We can now simulate the encoder. ```@example events @named sys = ODESystem( eqs, t, [theta, omega], params; continuous_events = [qAevt, qBevt]) -ss = structural_simplify(sys) +ss = mtkcompile(sys) prob = ODEProblem(ss, [theta => 0.0], (0.0, pi)) sol = solve(prob, Tsit5(); dtmax = 0.01) sol.ps[cnt] diff --git a/docs/src/basics/FAQ.md b/docs/src/basics/FAQ.md index 10671299c6..e83b1f1336 100644 --- a/docs/src/basics/FAQ.md +++ b/docs/src/basics/FAQ.md @@ -28,7 +28,7 @@ are similarly undocumented. Following is the list of behaviors that should be re - `setindex!(::MTKParameters, value, ::ParameterIndex)` can be used to set the value of a parameter with the given index. - `parameter_values(sys, sym)` will return a `ParameterIndex` object if `sys` has been - `complete`d (through `structural_simplify`, `complete` or `@mtkbuild`). + `complete`d (through `mtkcompile`, `complete` or `@mtkcompile`). - `copy(::MTKParameters)` is defined and duplicates the parameter object, including the memory used by the underlying buffers. @@ -194,7 +194,7 @@ p, replace, alias = SciMLStructures.canonicalize(Tunable(), prob.p) # ERROR: ArgumentError: SymbolicUtils.BasicSymbolic{Real}[xˍt(t)] are missing from the variable map. -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. +This error can come up after running `mtkcompile` 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. ```julia using ModelingToolkit @@ -206,7 +206,7 @@ eqs = [x1 + x2 + 1 ~ 0 x1 + D(x3) + x4 + 3 ~ 0 2 * D(D(x1)) + D(D(x2)) + D(D(x3)) + D(x4) + 4 ~ 0] @named sys = ODESystem(eqs, t) -sys = structural_simplify(sys) +sys = mtkcompile(sys) prob = ODEProblem(sys, [], (0, 1)) ``` @@ -237,7 +237,7 @@ using ModelingToolkit: t_nounits as t, D_nounits as D sts = @variables x1(t) = 0.0 eqs = [D(x1) ~ 1.1 * x1] -@mtkbuild sys = ODESystem(eqs, t) +@mtkcompile sys = ODESystem(eqs, t) prob = ODEProblem{false}(sys, [], (0, 1); u0_constructor = x -> SVector(x...)) ``` diff --git a/docs/src/basics/InputOutput.md b/docs/src/basics/InputOutput.md index 4dc5a3d50f..35a3885dbd 100644 --- a/docs/src/basics/InputOutput.md +++ b/docs/src/basics/InputOutput.md @@ -30,7 +30,7 @@ This function takes a vector of variables that are to be considered inputs, i.e. !!! note "Un-simplified system" - This function expects `sys` to be un-simplified, i.e., `structural_simplify` or `@mtkbuild` should not be called on the system before passing it into this function. `generate_control_function` calls a special version of `structural_simplify` internally. + This function expects `sys` to be un-simplified, i.e., `mtkcompile` or `@mtkcompile` should not be called on the system before passing it into this function. `generate_control_function` calls a special version of `mtkcompile` internally. ### Example: diff --git a/docs/src/basics/Linearization.md b/docs/src/basics/Linearization.md index 1c06ce72d4..27b8ec9903 100644 --- a/docs/src/basics/Linearization.md +++ b/docs/src/basics/Linearization.md @@ -29,7 +29,7 @@ eqs = [u ~ kp * (r - y) # P controller D(x) ~ -x + u # First-order plant y ~ x] # Output equation -@named sys = ODESystem(eqs, t) # Do not call @mtkbuild when linearizing +@named sys = ODESystem(eqs, t) # Do not call @mtkcompile when linearizing matrices, simplified_sys = linearize(sys, [r], [y]) # Linearize from r to y matrices ``` @@ -47,7 +47,7 @@ using ModelingToolkit: inputs, outputs !!! note "Un-simplified system" - Linearization expects `sys` to be un-simplified, i.e., `structural_simplify` or `@mtkbuild` should not be called on the system before linearizing. + Linearization expects `sys` to be un-simplified, i.e., `mtkcompile` or `@mtkcompile` should not be called on the system before linearizing. ## Operating point diff --git a/docs/src/basics/MTKLanguage.md b/docs/src/basics/MTKLanguage.md index ba6d2c34b5..05847581f6 100644 --- a/docs/src/basics/MTKLanguage.md +++ b/docs/src/basics/MTKLanguage.md @@ -1,4 +1,4 @@ -# [ModelingToolkit Language: Modeling with `@mtkmodel`, `@connectors` and `@mtkbuild`](@id mtk_language) +# [ModelingToolkit Language: Modeling with `@mtkmodel`, `@connectors` and `@mtkcompile`](@id mtk_language) ## MTK Model @@ -150,7 +150,7 @@ end - Whenever a parameter or variable has initial value, for example `v(t) = 0.0`, a symbolic variable named `v` with initial value 0.0 and a keyword argument `v`, with default value `nothing` are created.
This way, users can optionally pass new value of `v` while creating a component. ```julia -julia> @mtkbuild model_c1 = ModelC(; v = 2.0); +julia> @mtkcompile model_c1 = ModelC(; v = 2.0); julia> ModelingToolkit.getdefault(model_c1.v) 2.0 @@ -182,7 +182,7 @@ One or more partial systems can be extended in a higher system with `@extend` st ```@example mtkmodel-example using ModelingToolkit: getdefault -@mtkbuild model_c3 = ModelC(; model_a.k_array = [1.0, 2.0]) +@mtkcompile model_c3 = ModelC(; model_a.k_array = [1.0, 2.0]) getdefault(model_c3.model_a.k_array[1]) # 1.0 @@ -529,28 +529,28 @@ end ## Build structurally simplified models: -`@mtkbuild` builds an instance of a component and returns a structurally simplied system. +`@mtkcompile` builds an instance of a component and returns a structurally simplied system. ```julia -@mtkbuild sys = CustomModel() +@mtkcompile sys = CustomModel() ``` This is equivalent to: ```julia @named model = CustomModel() -sys = structural_simplify(model) +sys = mtkcompile(model) ``` -Pass keyword arguments to `structural_simplify` using the following syntax: +Pass keyword arguments to `mtkcompile` using the following syntax: ```julia -@mtkbuild sys=CustomModel() fully_determined=false +@mtkcompile sys=CustomModel() fully_determined=false ``` This is equivalent to: ```julia @named model = CustomModel() -sys = structural_simplify(model; fully_determined = false) +sys = mtkcompile(model; fully_determined = false) ``` diff --git a/docs/src/basics/Precompilation.md b/docs/src/basics/Precompilation.md index 97111f0d6b..0bf9a86653 100644 --- a/docs/src/basics/Precompilation.md +++ b/docs/src/basics/Precompilation.md @@ -22,7 +22,7 @@ using ModelingToolkit @variables x(ModelingToolkit.t_nounits) @named sys = ODESystem([ModelingToolkit.D_nounits(x) ~ -x + 1], ModelingToolkit.t_nounits) -prob = ODEProblem(structural_simplify(sys), [x => 30.0], (0, 100), [], +prob = ODEProblem(mtkcompile(sys), [x => 30.0], (0, 100), [], eval_expression = true, eval_module = @__MODULE__) end diff --git a/docs/src/basics/Validation.md b/docs/src/basics/Validation.md index 79c5d0d214..74715d351e 100644 --- a/docs/src/basics/Validation.md +++ b/docs/src/basics/Validation.md @@ -112,7 +112,7 @@ ps = @parameters s=-1 [unit = u"cm"] c=c [unit = u"cm"] eqs = [D(a) ~ dummycomplex(c, s);] sys = ODESystem( eqs, t, [sts...;], [ps...;], name = :sys, checks = ~ModelingToolkit.CheckUnits) -sys_simple = structural_simplify(sys) +sys_simple = mtkcompile(sys) ``` ## `DynamicQuantities` Literals From 750db43807fa40cb734b1253c2560a7dbe272ed0 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Fri, 23 May 2025 15:26:12 +0530 Subject: [PATCH 02/49] docs: update `examples` docs to new syntax --- docs/src/examples/higher_order.md | 8 +++---- .../modelingtoolkitize_index_reduction.md | 4 ++-- docs/src/examples/perturbation.md | 4 ++-- docs/src/examples/remake.md | 4 ++-- docs/src/examples/sparse_jacobians.md | 2 +- docs/src/examples/spring_mass.md | 24 +++++++++---------- docs/src/examples/tearing_parallelism.md | 24 +++++++++---------- 7 files changed, 35 insertions(+), 35 deletions(-) diff --git a/docs/src/examples/higher_order.md b/docs/src/examples/higher_order.md index fac707525f..95480e283b 100644 --- a/docs/src/examples/higher_order.md +++ b/docs/src/examples/higher_order.md @@ -4,7 +4,7 @@ ModelingToolkit has a system for transformations of mathematical systems. These transformations allow for symbolically changing the representation of the model to problems that are easier to numerically solve. One simple to demonstrate transformation, is -`structural_simplify`, which does a lot of tricks, one being the +`mtkcompile`, which does a lot of tricks, one being the transformation that turns an Nth order ODE into N coupled 1st order ODEs. @@ -32,7 +32,7 @@ using ModelingToolkit: t_nounits as t, D_nounits as D D(z) ~ x * y - β * z end end -@mtkbuild sys = SECOND_ORDER() +@mtkcompile sys = SECOND_ORDER() ``` The second order ODE has been automatically transformed to two first order ODEs. @@ -43,7 +43,7 @@ and this syntax extends to `N`-th order. Also, we can use `*` or `∘` to compos `Differential`s, like `Differential(t) * Differential(x)`. Now let's transform this into the `ODESystem` of first order components. -We do this by calling `structural_simplify`: +We do this by calling `mtkcompile`: Now we can directly numerically solve the lowered system. Note that, following the original problem, the solution requires knowing the @@ -54,7 +54,7 @@ but we still have to provide a value for the latter. ```@example orderlowering u0 = [D(sys.x) => 2.0] tspan = (0.0, 100.0) -prob = ODEProblem(sys, u0, tspan, [], jac = true) +prob = ODEProblem(sys, u0, tspan, jac = true) sol = solve(prob, Tsit5()) using Plots plot(sol, idxs = (sys.x, sys.y)) diff --git a/docs/src/examples/modelingtoolkitize_index_reduction.md b/docs/src/examples/modelingtoolkitize_index_reduction.md index b19ea46701..9759dc2081 100644 --- a/docs/src/examples/modelingtoolkitize_index_reduction.md +++ b/docs/src/examples/modelingtoolkitize_index_reduction.md @@ -29,7 +29,7 @@ p = [9.8, 1] tspan = (0, 10.0) pendulum_prob = ODEProblem(pendulum_fun!, u0, tspan, p) traced_sys = modelingtoolkitize(pendulum_prob) -pendulum_sys = structural_simplify(dae_index_lowering(traced_sys)) +pendulum_sys = mtkcompile(dae_index_lowering(traced_sys)) prob = ODEProblem(pendulum_sys, [], tspan) sol = solve(prob, Rodas5P(), abstol = 1e-8, reltol = 1e-8) plot(sol, idxs = unknowns(traced_sys)) @@ -157,7 +157,7 @@ numerical solver. Let's try that out: ```@example indexred traced_sys = modelingtoolkitize(pendulum_prob) -pendulum_sys = structural_simplify(dae_index_lowering(traced_sys)) +pendulum_sys = mtkcompile(dae_index_lowering(traced_sys)) prob = ODEProblem(pendulum_sys, Pair[], tspan) sol = solve(prob, Rodas5P()) diff --git a/docs/src/examples/perturbation.md b/docs/src/examples/perturbation.md index 0d1a493cb4..20cef4067c 100644 --- a/docs/src/examples/perturbation.md +++ b/docs/src/examples/perturbation.md @@ -43,7 +43,7 @@ eqs_pert = taylor_coeff(eq_pert, ϵ, 0:2) These are the ODEs we want to solve. Now construct an `ODESystem`, which automatically inserts dummy derivatives for the velocities: ```@example perturbation -@mtkbuild sys = ODESystem(eqs_pert, t) +@mtkcompile sys = System(eqs_pert, t) ``` To solve the `ODESystem`, we generate an `ODEProblem` with initial conditions $x(0) = 0$, and $ẋ(0) = 1$, and solve it: @@ -87,7 +87,7 @@ We follow the same steps as in the previous example to construct the `ODESystem` ```@example perturbation eq_pert = substitute(eq, x => x_series) eqs_pert = taylor_coeff(eq_pert, ϵ, 0:2) -@mtkbuild sys = ODESystem(eqs_pert, t) +@mtkcompile sys = System(eqs_pert, t) ``` We solve and plot it as in the previous example, and compare the solution with $ϵ=0.1$ to the exact solution $x(t, ϵ) = e^{-ϵ t} \sin(\sqrt{(1-ϵ^2)}\,t) / \sqrt{1-ϵ^2}$ of the unperturbed equation: diff --git a/docs/src/examples/remake.md b/docs/src/examples/remake.md index 91dba4d7ae..9d40a83dbf 100644 --- a/docs/src/examples/remake.md +++ b/docs/src/examples/remake.md @@ -14,7 +14,7 @@ using ModelingToolkit: t_nounits as t, D_nounits as D @variables x(t) y(t) eqs = [D(x) ~ (α - β * y) * x D(y) ~ (δ * x - γ) * y] -@mtkbuild odesys = ODESystem(eqs, t) +@mtkcompile odesys = System(eqs, t) ``` To create the "data" for optimization, we will solve the system with a known set of @@ -24,7 +24,7 @@ parameters. using OrdinaryDiffEq odeprob = ODEProblem( - odesys, [x => 1.0, y => 1.0], (0.0, 10.0), [α => 1.5, β => 1.0, γ => 3.0, δ => 1.0]) + odesys, [x => 1.0, y => 1.0, α => 1.5, β => 1.0, γ => 3.0, δ => 1.0], (0.0, 10.0)) timesteps = 0.0:0.1:10.0 sol = solve(odeprob, Tsit5(); saveat = timesteps) data = Array(sol) diff --git a/docs/src/examples/sparse_jacobians.md b/docs/src/examples/sparse_jacobians.md index 03bd80d432..79a93c3471 100644 --- a/docs/src/examples/sparse_jacobians.md +++ b/docs/src/examples/sparse_jacobians.md @@ -54,7 +54,7 @@ prob = ODEProblem(brusselator_2d_loop, u0, (0.0, 11.5), p) Now let's use `modelingtoolkitize` to generate the symbolic version: ```@example sparsejac -@mtkbuild sys = modelingtoolkitize(prob); +@mtkcompile sys = modelingtoolkitize(prob); nothing # hide ``` diff --git a/docs/src/examples/spring_mass.md b/docs/src/examples/spring_mass.md index 355e5c20b2..8d42592796 100644 --- a/docs/src/examples/spring_mass.md +++ b/docs/src/examples/spring_mass.md @@ -13,13 +13,13 @@ 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 eqs = scalarize(D.(pos) .~ v) - ODESystem(eqs, t, [pos..., v...], ps; name) + System(eqs, t, [pos..., v...], ps; name) end function Spring(; name, k = 1e4, l = 1.0) ps = @parameters k=k l=l @variables x(t), dir(t)[1:2] - ODESystem(Equation[], t, [x, dir...], ps; name) + System(Equation[], t, [x, dir...], ps; name) end function connect_spring(spring, a, b) @@ -43,9 +43,9 @@ g = [0.0, -9.81] eqs = [connect_spring(spring, mass.pos, center) scalarize(D.(mass.v) .~ spring_force(spring) / mass.m .+ g)] -@named _model = ODESystem(eqs, t, [spring.x; spring.dir; mass.pos], []) +@named _model = System(eqs, t, [spring.x; spring.dir; mass.pos], []) @named model = compose(_model, mass, spring) -sys = structural_simplify(model) +sys = mtkcompile(model) prob = ODEProblem(sys, [], (0.0, 3.0)) sol = solve(prob, Rosenbrock23()) @@ -56,18 +56,18 @@ plot(sol) ### Building the components -For each component, we use a Julia function that returns an `ODESystem`. At the top, we define the fundamental properties of a `Mass`: it has a mass `m`, a position `pos` and a velocity `v`. We also define that the velocity is the rate of change of position with respect to time. +For each component, we use a Julia function that returns an `System`. At the top, we define the fundamental properties of a `Mass`: it has a mass `m`, a position `pos` and a velocity `v`. We also define that the velocity is the rate of change of position with respect to time. ```@example component 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 eqs = scalarize(D.(pos) .~ v) - ODESystem(eqs, t, [pos..., v...], ps; name) + System(eqs, t, [pos..., v...], ps; name) end ``` -Note that this is an incompletely specified `ODESystem`. It cannot be simulated on its own, since the equations for the velocity `v[1:2](t)` are unknown. Notice the addition of a `name` keyword. This allows us to generate different masses with different names. A `Mass` can now be constructed as: +Note that this is an incompletely specified `System`. It cannot be simulated on its own, since the equations for the velocity `v[1:2](t)` are unknown. Notice the addition of a `name` keyword. This allows us to generate different masses with different names. A `Mass` can now be constructed as: ```@example component Mass(name = :mass1) @@ -85,7 +85,7 @@ Next, we build the spring component. It is characterized by the spring constant function Spring(; name, k = 1e4, l = 1.0) ps = @parameters k=k l=l @variables x(t), dir(t)[1:2] - ODESystem(Equation[], t, [x, dir...], ps; name) + System(Equation[], t, [x, dir...], ps; name) end ``` @@ -129,7 +129,7 @@ eqs = [connect_spring(spring, mass.pos, center) Finally, we can build the model using these equations and components. ```@example component -@named _model = ODESystem(eqs, t, [spring.x; spring.dir; mass.pos], []) +@named _model = System(eqs, t, [spring.x; spring.dir; mass.pos], []) @named model = compose(_model, mass, spring) ``` @@ -153,10 +153,10 @@ parameters(model) ### Simplifying and solving this system -This system can be solved directly as a DAE using [one of the DAE solvers from DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/solvers/dae_solve/). However, we can symbolically simplify the system first beforehand. Running `structural_simplify` eliminates unnecessary variables from the model to give the leanest numerical representation of the system. +This system can be solved directly as a DAE using [one of the DAE solvers from DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/solvers/dae_solve/). However, we can symbolically simplify the system first beforehand. Running `mtkcompile` eliminates unnecessary variables from the model to give the leanest numerical representation of the system. ```@example component -sys = structural_simplify(model) +sys = mtkcompile(model) equations(sys) ``` @@ -177,7 +177,7 @@ sol = solve(prob, Rosenbrock23()) plot(sol) ``` -What if we want the timeseries of a different variable? That information is not lost! Instead, `structural_simplify` simply changes unknown variables into `observed` variables. +What if we want the timeseries of a different variable? That information is not lost! Instead, `mtkcompile` simply changes unknown variables into `observed` variables. ```@example component observed(sys) diff --git a/docs/src/examples/tearing_parallelism.md b/docs/src/examples/tearing_parallelism.md index 9540e610bd..4ecc8d2a45 100644 --- a/docs/src/examples/tearing_parallelism.md +++ b/docs/src/examples/tearing_parallelism.md @@ -1,7 +1,7 @@ # Exposing More Parallelism By Tearing Algebraic Equations in ODESystems Sometimes it can be very non-trivial to parallelize a system. In this tutorial, -we will demonstrate how to make use of `structural_simplify` to expose more +we will demonstrate how to make use of `mtkcompile` to expose more parallelism in the solution process and parallelize the resulting simulation. ## The Component Library @@ -16,13 +16,13 @@ using ModelingToolkit: t_nounits as t, D_nounits as D # Basic electric components @connector function Pin(; name) @variables v(t)=1.0 i(t)=1.0 [connect = Flow] - ODESystem(Equation[], t, [v, i], [], name = name) + System(Equation[], t, [v, i], [], name = name) end function Ground(; name) @named g = Pin() eqs = [g.v ~ 0] - compose(ODESystem(eqs, t, [], [], name = name), g) + compose(System(eqs, t, [], [], name = name), g) end function ConstantVoltage(; name, V = 1.0) @@ -32,12 +32,12 @@ function ConstantVoltage(; name, V = 1.0) @parameters V = V eqs = [V ~ p.v - n.v 0 ~ p.i + n.i] - compose(ODESystem(eqs, t, [], [V], name = name), p, n) + compose(System(eqs, t, [], [V], name = name), p, n) end @connector function HeatPort(; name) @variables T(t)=293.15 Q_flow(t)=0.0 [connect = Flow] - ODESystem(Equation[], t, [T, Q_flow], [], name = name) + System(Equation[], t, [T, Q_flow], [], name = name) end function HeatingResistor(; name, R = 1.0, TAmbient = 293.15, alpha = 1.0) @@ -51,7 +51,7 @@ function HeatingResistor(; name, R = 1.0, TAmbient = 293.15, alpha = 1.0) h.Q_flow ~ -v * p.i # -LossPower v ~ p.v - n.v 0 ~ p.i + n.i] - compose(ODESystem(eqs, t, [v, RTherm], [R, TAmbient, alpha], + compose(System(eqs, t, [v, RTherm], [R, TAmbient, alpha], name = name), p, n, h) end @@ -62,7 +62,7 @@ function HeatCapacitor(; name, rho = 8050, V = 1, cp = 460, TAmbient = 293.15) eqs = [ D(h.T) ~ h.Q_flow / C ] - compose(ODESystem(eqs, t, [], [rho, V, cp], + compose(System(eqs, t, [], [rho, V, cp], name = name), h) end @@ -74,7 +74,7 @@ function Capacitor(; name, C = 1.0) eqs = [v ~ p.v - n.v 0 ~ p.i + n.i D(v) ~ p.i / C] - compose(ODESystem(eqs, t, [v], [C], + compose(System(eqs, t, [v], [C], name = name), p, n) end @@ -88,7 +88,7 @@ function parallel_rc_model(i; name, source, ground, R, C) connect(capacitor.n, source.n, ground.g) connect(resistor.h, heat_capacitor.h)] - compose(ODESystem(rc_eqs, t, name = Symbol(name, i)), + compose(System(rc_eqs, t, name = Symbol(name, i)), [resistor, capacitor, source, ground, heat_capacitor]) end ``` @@ -114,7 +114,7 @@ eqs = [ D(E) ~ sum(((i, sys),) -> getproperty(sys, Symbol(:resistor, i)).h.Q_flow, enumerate(rc_systems)) ] -@named _big_rc = ODESystem(eqs, t, [E], []) +@named _big_rc = System(eqs, t, [E], []) @named big_rc = compose(_big_rc, rc_systems) ``` @@ -122,7 +122,7 @@ Now let's say we want to expose a bit more parallelism via running tearing. How do we do that? ```@example tearing -sys = structural_simplify(big_rc) +sys = mtkcompile(big_rc) ``` Done, that's it. There's no more to it. @@ -175,5 +175,5 @@ so this is better than trying to do it by hand. After performing this, you can construct the `ODEProblem` and set `parallel_form` to use the exposed parallelism in multithreaded function -constructions, but this showcases why `structural_simplify` is so important +constructions, but this showcases why `mtkcompile` is so important to that process. From 554305c249b4b43002067feb1cb9e7a05c8dd2db Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Fri, 23 May 2025 15:28:52 +0530 Subject: [PATCH 03/49] docs: update `tutorials` docs to new syntax --- docs/src/tutorials/SampledData.md | 9 ++--- docs/src/tutorials/acausal_components.md | 6 +-- docs/src/tutorials/attractors.md | 6 +-- .../bifurcation_diagram_computation.md | 4 +- docs/src/tutorials/callable_params.md | 8 ++-- .../tutorials/change_independent_variable.md | 16 ++++---- docs/src/tutorials/discrete_system.md | 9 ++--- docs/src/tutorials/disturbance_modeling.md | 8 ++-- docs/src/tutorials/domain_connections.md | 38 +++++++++---------- docs/src/tutorials/fmi.md | 22 ++++++----- docs/src/tutorials/initialization.md | 32 ++++++++-------- docs/src/tutorials/linear_analysis.md | 6 ++- docs/src/tutorials/modelingtoolkitize.md | 2 +- docs/src/tutorials/nonlinear.md | 6 +-- docs/src/tutorials/ode_modeling.md | 32 ++++++++-------- docs/src/tutorials/optimization.md | 6 +-- .../tutorials/parameter_identifiability.md | 8 ++-- .../tutorials/programmatically_generating.md | 21 +++++----- docs/src/tutorials/stochastic_diffeq.md | 8 ++-- 19 files changed, 125 insertions(+), 122 deletions(-) diff --git a/docs/src/tutorials/SampledData.md b/docs/src/tutorials/SampledData.md index c700bae5c2..c2d6c9308d 100644 --- a/docs/src/tutorials/SampledData.md +++ b/docs/src/tutorials/SampledData.md @@ -155,10 +155,9 @@ k = ShiftIndex(clock) function plant(; name) @variables x(t)=1 u(t)=0 y(t)=0 - D = Differential(t) eqs = [D(x) ~ -x + u y ~ x] - ODESystem(eqs, t; name = name) + System(eqs, t; name = name) end function filt(; name) # Reference filter @@ -166,7 +165,7 @@ function filt(; name) # Reference filter a = 1 / exp(dt) eqs = [x(k) ~ a * x(k - 1) + (1 - a) * u(k) y ~ x] - ODESystem(eqs, t, name = name) + System(eqs, t, name = name) end function controller(kp; name) @@ -174,7 +173,7 @@ function controller(kp; name) @parameters kp = kp eqs = [yd ~ Sample(y) ud ~ kp * (r - yd)] - ODESystem(eqs, t; name = name) + System(eqs, t; name = name) end @named f = filt() @@ -187,7 +186,7 @@ connections = [r ~ sin(t) # reference signal Hold(c.ud) ~ p.u # controller output to plant input (zero-order-hold) p.y ~ c.y] # plant output to controller feedback -@named cl = ODESystem(connections, t, systems = [f, c, p]) +@named cl = System(connections, t, systems = [f, c, p]) ``` ```@docs diff --git a/docs/src/tutorials/acausal_components.md b/docs/src/tutorials/acausal_components.md index 751b678dae..b364be4012 100644 --- a/docs/src/tutorials/acausal_components.md +++ b/docs/src/tutorials/acausal_components.md @@ -99,7 +99,7 @@ end end end -@mtkbuild rc_model = RCModel(resistor.R = 2.0) +@mtkcompile rc_model = RCModel(resistor.R = 2.0) u0 = [ rc_model.capacitor.v => 0.0 ] @@ -272,7 +272,7 @@ We can create a RCModel component with `@named`. And using `subcomponent_name.pa the parameters or defaults values of variables of subcomponents. ```@example acausal -@mtkbuild rc_model = RCModel(resistor.R = 2.0) +@mtkcompile rc_model = RCModel(resistor.R = 2.0) ``` This model is acausal because we have not specified anything about the causality of the model. We have @@ -323,7 +323,7 @@ plot(sol) By default, this plots only the unknown variables that had to be solved for. However, what if we wanted to plot the timeseries of a different variable? Do not worry, that information was not thrown away! Instead, transformations -like `structural_simplify` simply change unknown variables into observables which are +like `mtkcompile` simply change unknown variables into observables which are defined by `observed` equations. ```@example acausal diff --git a/docs/src/tutorials/attractors.md b/docs/src/tutorials/attractors.md index 317384b01a..426551b017 100644 --- a/docs/src/tutorials/attractors.md +++ b/docs/src/tutorials/attractors.md @@ -41,10 +41,10 @@ eqs = [ Because our dynamical system is super simple, we will directly make an `ODESystem` and cast it in an `ODEProblem` as in the [`ODESystems` tutorial](@ref programmatically). Since all state variables and parameters have a default value we can immediately write ```@example Attractors -@named modlorenz = ODESystem(eqs, t) -ssys = structural_simplify(modlorenz) +@named modlorenz = System(eqs, t) +ssys = mtkcompile(modlorenz) # The timespan given to the problem is irrelevant for DynamicalSystems.jl -prob = ODEProblem(ssys, [], (0.0, 1.0), []) +prob = ODEProblem(ssys, [], (0.0, 1.0)) ``` This `prob` can be turned to a dynamical system as simply as diff --git a/docs/src/tutorials/bifurcation_diagram_computation.md b/docs/src/tutorials/bifurcation_diagram_computation.md index f15d46e1e4..811ea89e83 100644 --- a/docs/src/tutorials/bifurcation_diagram_computation.md +++ b/docs/src/tutorials/bifurcation_diagram_computation.md @@ -14,7 +14,7 @@ using ModelingToolkit: t_nounits as t, D_nounits as D @parameters μ α eqs = [0 ~ μ * x - x^3 + α * y, 0 ~ -y] -@mtkbuild nsys = NonlinearSystem(eqs, [x, y], [μ, α]) +@mtkcompile nsys = System(eqs, [x, y], [μ, α]) ``` we wish to compute a bifurcation diagram for this system as we vary the parameter `μ`. For this, we need to provide the following information: @@ -95,7 +95,7 @@ using ModelingToolkit: t_nounits as t, D_nounits as D @parameters μ eqs = [D(x) ~ μ * x - y - x * (x^2 + y^2), D(y) ~ x + μ * y - y * (x^2 + y^2)] -@mtkbuild osys = ODESystem(eqs, t) +@mtkcompile osys = System(eqs, t) bif_par = μ plot_var = x diff --git a/docs/src/tutorials/callable_params.md b/docs/src/tutorials/callable_params.md index 74e7ea87fa..2500c27015 100644 --- a/docs/src/tutorials/callable_params.md +++ b/docs/src/tutorials/callable_params.md @@ -2,7 +2,7 @@ ModelingToolkit.jl allows creating parameters that represent functions to be called. This is especially useful for including interpolants and/or lookup tables inside ODEs. In this -tutorial we will create an `ODESystem` which employs callable parameters to interpolate data +tutorial we will create an `System` which employs callable parameters to interpolate data inside an ODE and go over the various syntax options and their implications. ## Callable parameter syntax @@ -69,7 +69,7 @@ Tspline = typeof(spline) @variables x(t) @parameters (interp::Tspline)(..) -@mtkbuild sys = ODESystem(D(x) ~ interp(t), t) +@mtkcompile sys = System(D(x) ~ interp(t), t) ``` The derivative of `x` is obtained via an interpolation from DataInterpolations.jl. Note @@ -77,7 +77,7 @@ the parameter syntax. The `(..)` marks the parameter as callable. `(interp::Tspl indicates that the parameter is of type `Tspline`. ```@example callable -prob = ODEProblem(sys, [x => 0.0], (0.0, 1.0), [interp => spline]) +prob = ODEProblem(sys, [x => 0.0, interp => spline], (0.0, 1.0)) solve(prob) ``` @@ -85,7 +85,7 @@ Note that the the following will not work: ```julia ODEProblem( - sys; [x => 0.0], (0.0, 1.0), [interp => LinearInterpolation(0.0:0.1:1.0, 0.0:0.1:1.0)]) + sys; [x => 0.0, interp => LinearInterpolation(0.0:0.1:1.0, 0.0:0.1:1.0)], (0.0, 1.0)) ``` Since the type of the spline doesn't match. diff --git a/docs/src/tutorials/change_independent_variable.md b/docs/src/tutorials/change_independent_variable.md index d55639a669..b5b548a73b 100644 --- a/docs/src/tutorials/change_independent_variable.md +++ b/docs/src/tutorials/change_independent_variable.md @@ -24,8 +24,8 @@ using ModelingToolkit: t_nounits as t, D_nounits as D @parameters g=9.81 v # gravitational acceleration and horizontal velocity eqs = [D(D(y)) ~ -g, D(x) ~ v] initialization_eqs = [D(x) ~ D(y)] # 45° initial angle -M1 = ODESystem(eqs, t; initialization_eqs, name = :M) -M1s = structural_simplify(M1) +M1 = System(eqs, t; initialization_eqs, name = :M) +M1s = mtkcompile(M1) @assert length(equations(M1s)) == 3 # hide M1s # hide ``` @@ -44,7 +44,7 @@ This transformation is well-defined for any non-zero horizontal velocity $v$, so ```@example changeivar M2 = change_independent_variable(M1, x) -M2s = structural_simplify(M2; allow_symbolic = true) +M2s = mtkcompile(M2; allow_symbolic = true) # a sanity test on the 10 x/y variables that are accessible to the user # hide @assert allequal([x, M1s.x]) # hide @assert allequal([M2.x, M2s.x]) # hide @@ -67,7 +67,7 @@ It is straightforward to evolve the ODE for 10 meters and plot the resulting tra ```@example changeivar using OrdinaryDiffEq, Plots -prob = ODEProblem(M2s, [M2s.y => 0.0], [0.0, 10.0], [v => 8.0]) # throw 10 meters +prob = ODEProblem(M2s, [M2s.y => 0.0, v => 8.0], [0.0, 10.0]) # throw 10 meters sol = solve(prob, Tsit5()) plot(sol; idxs = M2.y) # must index by M2.y = y(x); not M1.y = y(t)! ``` @@ -88,16 +88,16 @@ In terms of conformal time $t$, they can be written a = GlobalScope(a) # global var needed by all species function species(w; kw...) eqs = [D(Ω) ~ -3(1 + w) * D(a) / a * Ω] - return ODESystem(eqs, t, [Ω], []; kw...) + return System(eqs, t, [Ω], []; kw...) end @named r = species(1 // 3) # radiation @named m = species(0) # matter @named Λ = species(-1) # dark energy / cosmological constant eqs = [Ω ~ r.Ω + m.Ω + Λ.Ω, D(a) ~ √(Ω) * a^2] initialization_eqs = [Λ.Ω + r.Ω + m.Ω ~ 1] -M1 = ODESystem(eqs, t, [Ω, a], []; initialization_eqs, name = :M) +M1 = System(eqs, t, [Ω, a], []; initialization_eqs, name = :M) M1 = compose(M1, r, m, Λ) -M1s = structural_simplify(M1) +M1s = mtkcompile(M1) ``` Of course, we can solve this ODE as it is: @@ -137,7 +137,7 @@ M3 = change_independent_variable(M2, b, [Da(b) ~ exp(-b), a ~ exp(b)]) We can now solve and plot the ODE in terms of $b$: ```@example changeivar -M3s = structural_simplify(M3; allow_symbolic = true) +M3s = mtkcompile(M3; allow_symbolic = true) prob = ODEProblem(M3s, [M3s.r.Ω => 5e-5, M3s.m.Ω => 0.3], (0, -15), []) sol = solve(prob, Tsit5()) @assert Symbol(sol.retcode) == :Success # surrounding text assumes the solution was successful # hide diff --git a/docs/src/tutorials/discrete_system.md b/docs/src/tutorials/discrete_system.md index 8f6828fde7..d070156076 100644 --- a/docs/src/tutorials/discrete_system.md +++ b/docs/src/tutorials/discrete_system.md @@ -1,7 +1,6 @@ # (Experimental) Modeling Discrete Systems -In this example, we will use the new [`DiscreteSystem`](@ref) API -to create an SIR model. +In this example, we will use the [`System`](@ref) API to create an SIR model. ```@example discrete using ModelingToolkit @@ -23,12 +22,12 @@ recovery = rate_to_proportion(γ * h, δt) * I(k - 1) eqs = [S(k) ~ S(k - 1) - infection * h, I(k) ~ I(k - 1) + infection - recovery, R(k) ~ R(k - 1) + recovery] -@mtkbuild sys = DiscreteSystem(eqs, t) +@mtkcompile sys = System(eqs, t) u0 = [S(k - 1) => 990.0, I(k - 1) => 10.0, R(k - 1) => 0.0] p = [β => 0.05, c => 10.0, γ => 0.25, δt => 0.1] tspan = (0.0, 100.0) -prob = DiscreteProblem(sys, u0, tspan, p) +prob = DiscreteProblem(sys, vcat(u0, p), tspan) sol = solve(prob, FunctionMap()) ``` @@ -39,7 +38,7 @@ the Fibonacci series: ```@example discrete @variables x(t) = 1.0 -@mtkbuild sys = DiscreteSystem([x ~ x(k - 1) + x(k - 2)], t) +@mtkcompile sys = System([x ~ x(k - 1) + x(k - 2)], t) ``` The "default value" here should be interpreted as the value of `x` at all past timesteps. diff --git a/docs/src/tutorials/disturbance_modeling.md b/docs/src/tutorials/disturbance_modeling.md index db8d926498..41c7bd86ee 100644 --- a/docs/src/tutorials/disturbance_modeling.md +++ b/docs/src/tutorials/disturbance_modeling.md @@ -67,7 +67,7 @@ This outer model, `ModelWithInputs`, contains two disturbance inputs, both of ty ```@example DISTURBANCE_MODELING @named model = ModelWithInputs() # Model with load disturbance -ssys = structural_simplify(model) +ssys = mtkcompile(model) prob = ODEProblem(ssys, [], (0.0, 6.0)) sol = solve(prob, Tsit5()) using Plots @@ -108,7 +108,7 @@ y &= g(x, p, t) \end{aligned} ``` -where ``x`` is the state, ``y`` are observed variables, ``p`` are parameters, and ``t`` is time. When using MTK, which variables constitute ``x`` and which are considered part of the output, ``y``, is up to the tool rather than the user, this choice is made by MTK during the call to `@mtkbuild` or the lower-level function `structural_simplify`. +where ``x`` is the state, ``y`` are observed variables, ``p`` are parameters, and ``t`` is time. When using MTK, which variables constitute ``x`` and which are considered part of the output, ``y``, is up to the tool rather than the user, this choice is made by MTK during the call to `@mtkcompile` or the lower-level function `mtkcompile`. If we further consider external inputs to the system, such as controlled input signals ``u`` and disturbance inputs ``w``, we can write the system as @@ -169,7 +169,7 @@ end We demonstrate that this model is complete and can be simulated: ```@example DISTURBANCE_MODELING -ssys = structural_simplify(model_with_disturbance) +ssys = mtkcompile(model_with_disturbance) prob = ODEProblem(ssys, [], (0.0, 10.0)) sol = solve(prob, Tsit5()) using Test @@ -191,7 +191,7 @@ g = ModelingToolkit.build_explicit_observed_function( io_sys, outputs; inputs) op = ModelingToolkit.inputs(io_sys) .=> 0 -x0, _ = ModelingToolkit.get_u0_p(io_sys, op, op) +x0 = ModelingToolkit.get_u0(io_sys, op) p = MTKParameters(io_sys, op) u = zeros(1) # Control input w = zeros(length(disturbance_inputs)) # Disturbance input diff --git a/docs/src/tutorials/domain_connections.md b/docs/src/tutorials/domain_connections.md index d6dc2d8781..f9cd3040d2 100644 --- a/docs/src/tutorials/domain_connections.md +++ b/docs/src/tutorials/domain_connections.md @@ -20,7 +20,7 @@ using ModelingToolkit: t_nounits as t, D_nounits as D dm(t), [connect = Flow] end - ODESystem(Equation[], t, vars, pars; name, defaults = [dm => 0]) + System(Equation[], t, vars, pars; name, defaults = [dm => 0]) end nothing #hide ``` @@ -47,7 +47,7 @@ The fluid medium setter for `HydralicPort` may be defined as `HydraulicFluid` wi dm ~ 0 ] - ODESystem(eqs, t, vars, pars; name, defaults = [dm => 0]) + System(eqs, t, vars, pars; name, defaults = [dm => 0]) end nothing #hide ``` @@ -63,7 +63,7 @@ Now, we can connect a `HydraulicFluid` component to any `HydraulicPort` connecto eqs = [port.p ~ p] - ODESystem(eqs, t, [], pars; name, systems) + System(eqs, t, [], pars; name, systems) end @component function FixedVolume(; vol, p_int, name) @@ -89,7 +89,7 @@ end rho ~ port.ρ * (1 + p / port.β) dm ~ drho * vol] - ODESystem(eqs, t, vars, pars; name, systems) + System(eqs, t, vars, pars; name, systems) end nothing #hide ``` @@ -97,7 +97,7 @@ nothing #hide When the system is defined we can generate a fluid component and connect it to the system. Here `fluid` is connected to the `src.port`, but it could also be connected to `vol.port`, any connection in the network is fine. ```@example domain -@component function System(; name) +@component function HydraulicSystem(; name) systems = @named begin src = FixedPressure(; p = 200e5) vol = FixedVolume(; vol = 0.1, p_int = 200e5) @@ -108,17 +108,17 @@ When the system is defined we can generate a fluid component and connect it to t eqs = [connect(fluid, src.port) connect(src.port, vol.port)] - ODESystem(eqs, t, [], []; systems, name) + System(eqs, t, [], []; systems, name) end -@named odesys = System() +@named odesys = HydraulicSystem() nothing #hide ``` -To see how the domain works, we can examine the set parameter values for each of the ports `src.port` and `vol.port`. First we assemble the system using `structural_simplify()` and then check the default value of `vol.port.ρ`, whichs points to the setter value `fluid₊ρ`. Likewise, `src.port.ρ`, will also point to the setter value `fluid₊ρ`. Therefore, there is now only 1 defined density value `fluid₊ρ` which sets the density for the connected network. +To see how the domain works, we can examine the set parameter values for each of the ports `src.port` and `vol.port`. First we assemble the system using `mtkcompile()` and then check the default value of `vol.port.ρ`, whichs points to the setter value `fluid₊ρ`. Likewise, `src.port.ρ`, will also point to the setter value `fluid₊ρ`. Therefore, there is now only 1 defined density value `fluid₊ρ` which sets the density for the connected network. ```@repl domain -sys = structural_simplify(odesys) +sys = mtkcompile(odesys) ModelingToolkit.defaults(sys)[odesys.vol.port.ρ] ``` @@ -151,7 +151,7 @@ If we have a more complicated system, for example a hydraulic actuator, with a s port_a.dm ~ +(port_a.ρ) * dx * area port_b.dm ~ -(port_b.ρ) * dx * area] - ODESystem(eqs, t, vars, pars; name, systems) + System(eqs, t, vars, pars; name, systems) end nothing #hide ``` @@ -174,14 +174,14 @@ A system with 2 different fluids is defined and connected to each separate domai connect(src_a.port, act.port_a) connect(src_b.port, act.port_b)] - ODESystem(eqs, t, [], []; systems, name) + System(eqs, t, [], []; systems, name) end @named actsys2 = ActuatorSystem2() nothing #hide ``` -After running `structural_simplify()` on `actsys2`, the defaults will show that `act.port_a.ρ` points to `fluid_a₊ρ` and `act.port_b.ρ` points to `fluid_b₊ρ`. This is a special case, in most cases a hydraulic system will have only 1 fluid, however this simple system has 2 separate domain networks. Therefore, we can connect a single fluid to both networks. This does not interfere with the mathematical equations of the system, since no unknown variables are connected. +After running `mtkcompile()` on `actsys2`, the defaults will show that `act.port_a.ρ` points to `fluid_a₊ρ` and `act.port_b.ρ` points to `fluid_b₊ρ`. This is a special case, in most cases a hydraulic system will have only 1 fluid, however this simple system has 2 separate domain networks. Therefore, we can connect a single fluid to both networks. This does not interfere with the mathematical equations of the system, since no unknown variables are connected. ```@example domain @component function ActuatorSystem1(; name) @@ -198,7 +198,7 @@ After running `structural_simplify()` on `actsys2`, the defaults will show that connect(src_a.port, act.port_a) connect(src_b.port, act.port_b)] - ODESystem(eqs, t, [], []; systems, name) + System(eqs, t, [], []; systems, name) end @named actsys1 = ActuatorSystem1() @@ -224,7 +224,7 @@ In some cases a component will be defined with 2 connectors of the same domain, eqs = [port_a.dm ~ (port_a.p - port_b.p) * K 0 ~ port_a.dm + port_b.dm] - ODESystem(eqs, t, [], pars; systems, name) + System(eqs, t, [], pars; systems, name) end nothing #hide ``` @@ -245,14 +245,14 @@ Adding the `Restrictor` to the original system example will cause a break in the connect(src.port, res.port_a) connect(res.port_b, vol.port)] - ODESystem(eqs, t, [], []; systems, name) + System(eqs, t, [], []; systems, name) end -@mtkbuild ressys = RestrictorSystem() +@mtkcompile ressys = RestrictorSystem() nothing #hide ``` -When `structural_simplify()` is applied to this system it can be seen that the defaults are missing for `res.port_b` and `vol.port`. +When `mtkcompile()` is applied to this system it can be seen that the defaults are missing for `res.port_b` and `vol.port`. ```@repl domain ModelingToolkit.defaults(ressys)[ressys.res.port_a.ρ] @@ -278,10 +278,10 @@ To ensure that the `Restrictor` component does not disrupt the domain network, t port_a.dm ~ (port_a.p - port_b.p) * K 0 ~ port_a.dm + port_b.dm] - ODESystem(eqs, t, [], pars; systems, name) + System(eqs, t, [], pars; systems, name) end -@mtkbuild ressys = RestrictorSystem() +@mtkcompile ressys = RestrictorSystem() nothing #hide ``` diff --git a/docs/src/tutorials/fmi.md b/docs/src/tutorials/fmi.md index 0e01393652..d17452f612 100644 --- a/docs/src/tutorials/fmi.md +++ b/docs/src/tutorials/fmi.md @@ -76,7 +76,7 @@ initialization semantics. We can simulate this model like any other ModelingToolkit system. ```@repl fmi -sys = structural_simplify(model) +sys = mtkcompile(model) prob = ODEProblem(sys, [sys.mass__s => 0.5, sys.mass__v => 0.0], (0.0, 5.0)) sol = solve(prob, Tsit5()) ``` @@ -105,11 +105,11 @@ constant until the next time the callback triggers. The periodic interval must b `communication_step_size` keyword argument. A smaller step size typically leads to less error but is more computationally expensive. -This model alone does not have any differential variables, and calling `structural_simplify` will lead +This model alone does not have any differential variables, and calling `mtkcompile` will lead to an `ODESystem` with no unknowns. ```@example fmi -structural_simplify(inner) +mtkcompile(inner) ``` Simulating this model will cause the OrdinaryDiffEq integrator to immediately finish, and will not @@ -117,7 +117,7 @@ trigger the callback. Thus, we wrap this system in a trivial system with a diffe ```@example fmi @variables x(t) = 1.0 -@mtkbuild sys = ODESystem([D(x) ~ x], t; systems = [inner]) +@mtkcompile sys = System([D(x) ~ x], t; systems = [inner]) ``` We can now simulate `sys`. @@ -184,7 +184,7 @@ metadata. We can now use this component as a subcomponent of a larger system. ```@repl fmi @variables a(t) b(t) c(t) [guess = 1.0]; -@mtkbuild sys = ODESystem( +@mtkcompile sys = System( [adder.a ~ a, adder.b ~ b, D(a) ~ t, D(b) ~ adder.out + adder.c, c^2 ~ adder.out + adder.value], t; @@ -198,8 +198,9 @@ FMUs in initialization to solve for initial states. As mentioned earlier, we can through an FMU. Thus, automatic differentiation has to be disabled for the solver. ```@example fmi -prob = ODEProblem(sys, [sys.adder.c => 2.0, sys.a => 1.0, sys.b => 1.0], - (0.0, 1.0), [sys.adder.value => 2.0]) +prob = ODEProblem( + sys, [sys.adder.c => 2.0, sys.a => 1.0, sys.b => 1.0, sys.adder.value => 2.0], + (0.0, 1.0)) solve(prob, Rodas5P(autodiff = false)) ``` @@ -217,12 +218,13 @@ fmu = loadFMU( @named adder = ModelingToolkit.FMIComponent( Val(2); fmu, type = :CS, communication_step_size = 1e-3, reinitializealg = BrownFullBasicInit()) -@mtkbuild sys = ODESystem( +@mtkcompile sys = System( [adder.a ~ a, adder.b ~ b, D(a) ~ t, D(b) ~ adder.out + adder.c, c^2 ~ adder.out + adder.value], t; systems = [adder]) -prob = ODEProblem(sys, [sys.adder.c => 2.0, sys.a => 1.0, sys.b => 1.0], - (0.0, 1.0), [sys.adder.value => 2.0]) +prob = ODEProblem( + sys, [sys.adder.c => 2.0, sys.a => 1.0, sys.b => 1.0, sys.adder.value => 2.0], + (0.0, 1.0)) solve(prob, Rodas5P()) ``` diff --git a/docs/src/tutorials/initialization.md b/docs/src/tutorials/initialization.md index ba733e0bfb..cdf6e7e813 100644 --- a/docs/src/tutorials/initialization.md +++ b/docs/src/tutorials/initialization.md @@ -44,7 +44,7 @@ using ModelingToolkit: t_nounits as t, D_nounits as D eqs = [D(D(x)) ~ λ * x D(D(y)) ~ λ * y - g x^2 + y^2 ~ 1] -@mtkbuild pend = ODESystem(eqs, t) +@mtkcompile pend = System(eqs, t) ``` While we defined the system using second derivatives and a length constraint, @@ -59,7 +59,7 @@ To see how the system works, let's start the pendulum in the far right position, i.e. `x(0) = 1` and `y(0) = 0`. We can do this by: ```@example init -prob = ODEProblem(pend, [x => 1, y => 0], (0.0, 1.5), [g => 1], guesses = [λ => 1]) +prob = ODEProblem(pend, [x => 1, y => 0, g => 1], (0.0, 1.5), guesses = [λ => 1]) ``` This solves via: @@ -93,7 +93,7 @@ sol[λ][1] We can similarly choose `λ = 0` and solve for `y` to start the system: ```@example init -prob = ODEProblem(pend, [x => 1, λ => 0], (0.0, 1.5), [g => 1], guesses = [y => 1]) +prob = ODEProblem(pend, [x => 1, λ => 0, g => 1], (0.0, 1.5); guesses = [y => 1]) sol = solve(prob, Rodas5P()) plot(sol, idxs = (x, y)) ``` @@ -102,7 +102,7 @@ or choose to satisfy derivative conditions: ```@example init prob = ODEProblem( - pend, [x => 1, D(y) => 0], (0.0, 1.5), [g => 1], guesses = [λ => 0, y => 1]) + pend, [x => 1, D(y) => 0, g => 1], (0.0, 1.5); guesses = [λ => 0, y => 1]) sol = solve(prob, Rodas5P()) plot(sol, idxs = (x, y)) ``` @@ -114,7 +114,7 @@ We can also directly give equations to be satisfied at the initial point by usin the `initialization_eqs` keyword argument, for example: ```@example init -prob = ODEProblem(pend, [x => 1], (0.0, 1.5), [g => 1], guesses = [λ => 0, y => 1], +prob = ODEProblem(pend, [x => 1, g => 1], (0.0, 1.5); guesses = [λ => 0, y => 1], initialization_eqs = [y ~ 0]) sol = solve(prob, Rodas5P()) plot(sol, idxs = (x, y)) @@ -125,7 +125,7 @@ variables and parameters: ```@example init prob = ODEProblem( - pend, [x => 1, D(y) => g], (0.0, 3.0), [g => 1], guesses = [λ => 0, y => 1]) + pend, [x => 1, D(y) => g, g => 1], (0.0, 3.0); guesses = [λ => 0, y => 1]) sol = solve(prob, Rodas5P()) plot(sol, idxs = (x, y)) ``` @@ -141,7 +141,7 @@ conditions = getfield.(equations(pend)[3:end], :rhs) when we initialize with ```@example init -prob = ODEProblem(pend, [x => 1, y => 0], (0.0, 1.5), [g => 1], guesses = [y => 0, λ => 1]) +prob = ODEProblem(pend, [x => 1, y => 0, g => 1], (0.0, 1.5); guesses = [y => 0, λ => 1]) ``` we have two extra conditions to satisfy, `x ~ 1` and `y ~ 0` at the initial point. That gives @@ -149,7 +149,7 @@ we have two extra conditions to satisfy, `x ~ 1` and `y ~ 0` at the initial poin case? ```@example init -prob = ODEProblem(pend, [x => 1], (0.0, 1.5), [g => 1], guesses = [y => 0, λ => 1]) +prob = ODEProblem(pend, [x => 1, g => 1], (0.0, 1.5); guesses = [y => 0, λ => 1]) ``` Here we have 4 equations for 5 unknowns (note: the warning is post-simplification of the @@ -167,7 +167,7 @@ direction, we may have an overdetermined system: ```@example init prob = ODEProblem( - pend, [x => 1, y => 0.0, D(y) => 0], (0.0, 1.5), [g => 1], guesses = [λ => 1]) + pend, [x => 1, y => 0.0, D(y) => 0, g => 1], (0.0, 1.5); guesses = [λ => 1]) ``` Can that be solved? @@ -184,7 +184,7 @@ a `SciMLBase.ReturnCode.InitialFailure`: ```@example init prob = ODEProblem( - pend, [x => 1, y => 0.0, D(y) => 2.0, λ => 1], (0.0, 1.5), [g => 1], guesses = [λ => 1]) + pend, [x => 1, y => 0.0, D(y) => 2.0, λ => 1, g => 1], (0.0, 1.5); guesses = [λ => 1]) sol = solve(prob, Rodas5P()) ``` @@ -217,7 +217,7 @@ y-velocities from a given position. ```@example init prob = ODEProblem( - pend, [x => 1, D(y) => 0], (0.0, 1.5), [g => 1], guesses = [λ => 0, y => 1, x => 1]) + pend, [x => 1, D(y) => 0, g => 1], (0.0, 1.5); guesses = [λ => 0, y => 1, x => 1]) sol1 = solve(prob, Rodas5P()) ``` @@ -312,7 +312,7 @@ using ModelingToolkit: t_nounits as t, D_nounits as D # hidden @variables x(t) y(t) @parameters total -@mtkbuild sys = ODESystem([D(x) ~ -x, total ~ x + y], t; +@mtkcompile sys = System([D(x) ~ -x, total ~ x + y], t; defaults = [total => missing], guesses = [total => 1.0]) ``` @@ -380,7 +380,7 @@ with observables, those observables are too treated as initial equations. We can resulting simplified system via the command: ```@example init -isys = structural_simplify(isys; fully_determined = false) +isys = mtkcompile(isys; fully_determined = false) ``` Note `fully_determined=false` allows for the simplification to occur when the number of equations @@ -392,7 +392,7 @@ isys = ModelingToolkit.generate_initializesystem( ``` ```@example init -isys = structural_simplify(isys; fully_determined = false) +isys = mtkcompile(isys; fully_determined = false) ``` ```@example init @@ -503,8 +503,8 @@ eqs = [D(x) ~ α * x - β * x * y D(y) ~ -γ * y + δ * x * y z ~ x + y] -@named sys = ODESystem(eqs, t) -simpsys = structural_simplify(sys) +@named sys = System(eqs, t) +simpsys = mtkcompile(sys) tspan = (0.0, 10.0) ``` diff --git a/docs/src/tutorials/linear_analysis.md b/docs/src/tutorials/linear_analysis.md index 7fcb67f9fb..250dbaa6a7 100644 --- a/docs/src/tutorials/linear_analysis.md +++ b/docs/src/tutorials/linear_analysis.md @@ -58,12 +58,14 @@ The following example builds a simple closed-loop system with a plant $P$ and a ```@example LINEAR_ANALYSIS using ModelingToolkitStandardLibrary.Blocks, ModelingToolkit +using ModelingToolkit: t_nounits as t + @named P = FirstOrder(k = 1, T = 1) # A first-order system with pole in -1 @named C = Gain(-1) # A P controller -t = ModelingToolkit.get_iv(P) + eqs = [connect(P.output, :plant_output, C.input) # Connect with an automatically created analysis point called :plant_output connect(C.output, :plant_input, P.input)] -sys = ODESystem(eqs, t, systems = [P, C], name = :feedback_system) +sys = System(eqs, t, systems = [P, C], name = :feedback_system) matrices_S = get_sensitivity(sys, :plant_input)[1] # Compute the matrices of a state-space representation of the (input)sensitivity function. matrices_T = get_comp_sensitivity(sys, :plant_input)[1] diff --git a/docs/src/tutorials/modelingtoolkitize.md b/docs/src/tutorials/modelingtoolkitize.md index 545879f842..fc74578ea4 100644 --- a/docs/src/tutorials/modelingtoolkitize.md +++ b/docs/src/tutorials/modelingtoolkitize.md @@ -52,7 +52,7 @@ If we want to get a symbolic representation, we can simply call `modelingtoolkit on the `prob`, which will return an `ODESystem`: ```@example mtkize -@mtkbuild sys = modelingtoolkitize(prob) +@mtkcompile sys = modelingtoolkitize(prob) ``` Using this, we can symbolically build the Jacobian and then rebuild the ODEProblem: diff --git a/docs/src/tutorials/nonlinear.md b/docs/src/tutorials/nonlinear.md index 057e856229..13caf96231 100644 --- a/docs/src/tutorials/nonlinear.md +++ b/docs/src/tutorials/nonlinear.md @@ -23,12 +23,12 @@ using ModelingToolkit, NonlinearSolve eqs = [0 ~ σ * (y - x) 0 ~ x * (ρ - z) - y 0 ~ x * y - β * z] -@mtkbuild ns = NonlinearSystem(eqs) +@mtkcompile ns = System(eqs) guesses = [x => 1.0, y => 0.0, z => 0.0] ps = [σ => 10.0, ρ => 26.0, β => 8 / 3] -prob = NonlinearProblem(ns, guesses, ps) +prob = NonlinearProblem(ns, vcat(guesses, ps)) sol = solve(prob, NewtonRaphson()) ``` @@ -38,6 +38,6 @@ Just like with `ODEProblem`s we can generate the `NonlinearProblem` with its ana Jacobian function: ```@example nonlinear -prob = NonlinearProblem(ns, guesses, ps, jac = true) +prob = NonlinearProblem(ns, vcat(guesses, ps), jac = true) sol = solve(prob, NewtonRaphson()) ``` diff --git a/docs/src/tutorials/ode_modeling.md b/docs/src/tutorials/ode_modeling.md index 0a4cd80803..3113ab0fc2 100644 --- a/docs/src/tutorials/ode_modeling.md +++ b/docs/src/tutorials/ode_modeling.md @@ -35,7 +35,7 @@ using ModelingToolkit: t_nounits as t, D_nounits as D end using OrdinaryDiffEq -@mtkbuild fol = FOL() +@mtkcompile fol = FOL() prob = ODEProblem(fol, [], (0.0, 10.0), []) sol = solve(prob) @@ -77,12 +77,12 @@ using ModelingToolkit: t_nounits as t, D_nounits as D end end -@mtkbuild fol = FOL() +@mtkcompile fol = FOL() ``` Note that equations in MTK use the tilde character (`~`) as equality sign. -`@mtkbuild` creates an instance of `FOL` named as `fol`. +`@mtkcompile` creates an instance of `FOL` named as `fol`. After construction of the ODE, you can solve it using [OrdinaryDiffEq.jl](https://docs.sciml.ai/DiffEqDocs/stable/): @@ -90,7 +90,7 @@ After construction of the ODE, you can solve it using [OrdinaryDiffEq.jl](https: using OrdinaryDiffEq using Plots -prob = ODEProblem(fol, [], (0.0, 10.0), []) +prob = ODEProblem(fol, [], (0.0, 10.0)) plot(solve(prob)) ``` @@ -105,15 +105,15 @@ you likely do not want to write an entirely new `@mtkmodel`. ModelingToolkit supports overwriting the default values: ```@example ode2 -@mtkbuild fol_different_values = FOL(; τ = 1 / 3, x = 0.5) -prob = ODEProblem(fol_different_values, [], (0.0, 10.0), []) +@mtkcompile fol_different_values = FOL(; τ = 1 / 3, x = 0.5) +prob = ODEProblem(fol_different_values, [], (0.0, 10.0)) plot(solve(prob)) ``` Alternatively, this overwriting could also have occurred at the `ODEProblem` level. ```@example ode2 -prob = ODEProblem(fol, [fol.τ => 1 / 3], (0.0, 10.0), [fol.x => 0.5]) +prob = ODEProblem(fol, [fol.x => 0.5, fol.τ => 1 / 3], (0.0, 10.0)) plot(solve(prob)) ``` @@ -147,7 +147,7 @@ using ModelingToolkit: t_nounits as t, D_nounits as D end end -@mtkbuild fol = FOL() +@mtkcompile fol = FOL() ``` If you copy this block of code to your REPL, you will not see the above LaTeX equations. @@ -175,7 +175,7 @@ in a simulation result. The intermediate variable `RHS` therefore can be plotted along with the unknown variable. Note that this has to be requested explicitly: ```@example ode2 -prob = ODEProblem(fol, [], (0.0, 10.0), []) +prob = ODEProblem(fol, [], (0.0, 10.0)) sol = solve(prob) plot(sol, idxs = [fol.x, fol.RHS]) ``` @@ -221,7 +221,7 @@ Obviously, one could use an explicit, symbolic function of time: end end -@mtkbuild fol_variable_f = FOL() +@mtkcompile fol_variable_f = FOL() ``` However, this function might not be available in an explicit form. @@ -252,11 +252,11 @@ f_fun(t) = t >= 10 ? value_vector[end] : value_vector[Int(floor(t)) + 1] end end -@mtkbuild fol_external_f = FOLExternalFunction() +@mtkcompile fol_external_f = FOLExternalFunction() ``` ```@example ode2 -prob = ODEProblem(fol_external_f, [], (0.0, 10.0), []) +prob = ODEProblem(fol_external_f, [], (0.0, 10.0)) sol = solve(prob) plot(sol, idxs = [fol_external_f.x, fol_external_f.f]) ``` @@ -292,7 +292,7 @@ end fol_2.f ~ fol_1.x end end -@mtkbuild connected = FOLConnected() +@mtkcompile connected = FOLConnected() ``` Here the total model consists of two of the same submodels (components), @@ -316,7 +316,7 @@ initial unknown and the parameter values can be specified accordingly when building the `ODEProblem`: ```@example ode2 -prob = ODEProblem(connected, [], (0.0, 10.0), []) +prob = ODEProblem(connected, [], (0.0, 10.0)) plot(solve(prob)) ``` @@ -344,13 +344,13 @@ Now have MTK provide sparse, analytical derivatives to the solver. This has to be specified during the construction of the `ODEProblem`: ```@example ode2 -prob_an = ODEProblem(connected, [], (0.0, 10.0), []; jac = true) +prob_an = ODEProblem(connected, [], (0.0, 10.0); jac = true) @btime solve(prob_an, Rodas4()); nothing # hide ``` ```@example ode2 -prob_sparse = ODEProblem(connected, [], (0.0, 10.0), []; jac = true, sparse = true) +prob_sparse = ODEProblem(connected, [], (0.0, 10.0); jac = true, sparse = true) @btime solve(prob_sparse, Rodas4()); nothing # hide ``` diff --git a/docs/src/tutorials/optimization.md b/docs/src/tutorials/optimization.md index 20f1079dcd..9eb72b36ea 100644 --- a/docs/src/tutorials/optimization.md +++ b/docs/src/tutorials/optimization.md @@ -24,7 +24,7 @@ using ModelingToolkit, Optimization, OptimizationOptimJL end @parameters a=1.0 b=1.0 rosenbrock = (a - x)^2 + b * (y - x^2)^2 -@mtkbuild sys = OptimizationSystem(rosenbrock, [x, y], [a, b]) +@mtkcompile sys = OptimizationSystem(rosenbrock, [x, y], [a, b]) ``` Every optimization problem consists of a set of optimization variables. @@ -52,7 +52,7 @@ ModelingToolkit is also capable of constructing analytical gradients and Hessian u0 = [y => 2.0] p = [b => 100.0] -prob = OptimizationProblem(sys, u0, p, grad = true, hess = true) +prob = OptimizationProblem(sys, vcat(u0, p), grad = true, hess = true) u_opt = solve(prob, GradientDescent()) ``` @@ -86,7 +86,7 @@ rosenbrock = (a - x)^2 + b * (y - x^2)^2 cons = [ x^2 + y^2 ≲ 1 ] -@mtkbuild sys = OptimizationSystem(rosenbrock, [x, y], [a, b], constraints = cons) +@mtkcompile sys = OptimizationSystem(rosenbrock, [x, y], [a, b], constraints = cons) prob = OptimizationProblem(sys, [], grad = true, hess = true, cons_j = true, cons_h = true) u_opt = solve(prob, IPNewton()) ``` diff --git a/docs/src/tutorials/parameter_identifiability.md b/docs/src/tutorials/parameter_identifiability.md index 0875a7698c..1731d5795f 100644 --- a/docs/src/tutorials/parameter_identifiability.md +++ b/docs/src/tutorials/parameter_identifiability.md @@ -4,7 +4,7 @@ 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_. -The package has a standalone data structure for ordinary differential equations, but is also compatible with `ODESystem` type from `ModelingToolkit.jl`. +The package has a standalone data structure for ordinary differential equations, but is also compatible with `System` type from `ModelingToolkit.jl`. ## Local Identifiability @@ -22,7 +22,7 @@ y_2 = x_5\end{cases}$$ This model describes the biohydrogenation[^1] process[^2] with unknown initial conditions. -### Using the `ODESystem` object +### Using the `System` object To define the ode system in Julia, we use `ModelingToolkit.jl`. @@ -61,7 +61,7 @@ using ModelingToolkit: t_nounits as t, D_nounits as D end # define the system -@mtkbuild de = Biohydrogenation() +@mtkcompile de = Biohydrogenation() ``` After that, we are ready to check the system for local identifiability: @@ -179,7 +179,7 @@ using ModelingToolkit: t_nounits as t, D_nounits as D end end -@mtkbuild ode = GoodwinOscillator() +@mtkcompile ode = GoodwinOscillator() # check only 2 parameters to_check = [ode.b, ode.c] diff --git a/docs/src/tutorials/programmatically_generating.md b/docs/src/tutorials/programmatically_generating.md index 9fc1db1834..406f65d8d9 100644 --- a/docs/src/tutorials/programmatically_generating.md +++ b/docs/src/tutorials/programmatically_generating.md @@ -18,7 +18,8 @@ as demonstrated in the Symbolics.jl documentation. This looks like: ```@example scripting using ModelingToolkit # reexports Symbolics -@variables t x(t) y(t) # Define variables +@independent_variables t +@variables x(t) y(t) # Define variables D = Differential(t) eqs = [D(x) ~ y D(y) ~ x] # Define an array of equations @@ -43,11 +44,11 @@ using ModelingToolkit: t_nounits as t, D_nounits as D eqs = [D(x) ~ (h - x) / τ] # create an array of equations # your first ODE, consisting of a single equation, indicated by ~ -@named model = ODESystem(eqs, t) +@named model = System(eqs, t) # Perform the standard transformations and mark the model complete # Note: Complete models cannot be subsystems of other models! -fol = structural_simplify(model) +fol = mtkcompile(model) prob = ODEProblem(fol, [], (0.0, 10.0), []) using OrdinaryDiffEq sol = solve(prob) @@ -56,22 +57,22 @@ using Plots plot(sol) ``` -As you can see, generating an ODESystem is as simple as creating an array of equations -and passing it to the `ODESystem` constructor. +As you can see, generating an `System` is as simple as creating an array of equations +and passing it to the `System` constructor. -`@named` automatically gives a name to the `ODESystem`, and is shorthand for +`@named` automatically gives a name to the `System`, and is shorthand for ```@example scripting -fol_model = ODESystem(eqs, t; name = :fol_model) # @named fol_model = ODESystem(eqs, t) +fol_model = System(eqs, t; name = :fol_model) # @named fol_model = ODESystem(eqs, t) ``` -Thus, if we had read a name from a file and wish to populate an `ODESystem` with said name, we could do: +Thus, if we had read a name from a file and wish to populate an `System` with said name, we could do: ```@example scripting namesym = :name_from_file -fol_model = ODESystem(eqs, t; name = namesym) +fol_model = System(eqs, t; name = namesym) ``` ## Warning About Mutation -Be advsied that it's never a good idea to mutate an `ODESystem`, or any `AbstractSystem`. +Be advsied that it's never a good idea to mutate an `System`, or any `AbstractSystem`. diff --git a/docs/src/tutorials/stochastic_diffeq.md b/docs/src/tutorials/stochastic_diffeq.md index 79c71a2e8d..9bc0086d35 100644 --- a/docs/src/tutorials/stochastic_diffeq.md +++ b/docs/src/tutorials/stochastic_diffeq.md @@ -53,7 +53,7 @@ eqs = [D(x) ~ σ * (y - x) + 0.3x * B, D(y) ~ x * (ρ - z) - y + 0.3y * B, D(z) ~ x * y - β * z + 0.3z * B] -@mtkbuild de = System(eqs, t) +@mtkcompile de = System(eqs, t) ``` Even though we did not explicitly use `SDESystem`, ModelingToolkit can still infer this from the equations. @@ -65,7 +65,7 @@ typeof(de) We continue by solving and plotting the SDE. ```@example SDE -prob = SDEProblem(de, [], (0.0, 100.0), []) +prob = SDEProblem(de, [], (0.0, 100.0)) sol = solve(prob, SRIW1()) plot(sol, idxs = [(1, 2, 3)]) ``` @@ -87,8 +87,8 @@ multiple `@brownian` variables have to be declared. eqs = [D(x) ~ σ * (y - x) + 0.3x * Bx, D(y) ~ x * (ρ - z) - y + 0.3y * By, D(z) ~ x * y - β * z + 0.3z * Bz] -@mtkbuild de = System(eqs, t) -prob = SDEProblem(de, [], (0.0, 100.0), []) +@mtkcompile de = System(eqs, t) +prob = SDEProblem(de, [], (0.0, 100.0)) sol = solve(prob, SRIW1()) plot(sol) ``` From 073403976e93a3d4bfc57ee309fdc4dce5d935a8 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Fri, 23 May 2025 15:30:16 +0530 Subject: [PATCH 04/49] docs: remove old system docs --- docs/src/systems/DiscreteSystem.md | 33 ---------- docs/src/systems/ImplicitDiscreteSystem.md | 33 ---------- docs/src/systems/JumpSystem.md | 33 ---------- docs/src/systems/NonlinearSystem.md | 57 ----------------- docs/src/systems/ODESystem.md | 74 ---------------------- docs/src/systems/OptimizationSystem.md | 40 ------------ docs/src/systems/SDESystem.md | 70 -------------------- 7 files changed, 340 deletions(-) delete mode 100644 docs/src/systems/DiscreteSystem.md delete mode 100644 docs/src/systems/ImplicitDiscreteSystem.md delete mode 100644 docs/src/systems/JumpSystem.md delete mode 100644 docs/src/systems/NonlinearSystem.md delete mode 100644 docs/src/systems/ODESystem.md delete mode 100644 docs/src/systems/OptimizationSystem.md delete mode 100644 docs/src/systems/SDESystem.md diff --git a/docs/src/systems/DiscreteSystem.md b/docs/src/systems/DiscreteSystem.md deleted file mode 100644 index 55a02e5714..0000000000 --- a/docs/src/systems/DiscreteSystem.md +++ /dev/null @@ -1,33 +0,0 @@ -# DiscreteSystem - -## System Constructors - -```@docs -DiscreteSystem -``` - -## Composition and Accessor Functions - - - `get_eqs(sys)` or `equations(sys)`: The equations that define the discrete system. - - `get_unknowns(sys)` or `unknowns(sys)`: The set of unknowns in the discrete system. - - `get_ps(sys)` or `parameters(sys)`: The parameters of the discrete system. - - `get_iv(sys)`: The independent variable of the discrete system - -## Transformations - -```@docs; canonical=false -structural_simplify -``` - -## Problem Constructors - -```@docs; canonical=false -DiscreteProblem(sys::DiscreteSystem, u0map, tspan) -DiscreteFunction(sys::DiscreteSystem, args...) -``` - -## Discrete Domain - -```@docs; canonical=false -Shift -``` diff --git a/docs/src/systems/ImplicitDiscreteSystem.md b/docs/src/systems/ImplicitDiscreteSystem.md deleted file mode 100644 index d687502b49..0000000000 --- a/docs/src/systems/ImplicitDiscreteSystem.md +++ /dev/null @@ -1,33 +0,0 @@ -# ImplicitDiscreteSystem - -## System Constructors - -```@docs -ImplicitDiscreteSystem -``` - -## Composition and Accessor Functions - - - `get_eqs(sys)` or `equations(sys)`: The equations that define the implicit discrete system. - - `get_unknowns(sys)` or `unknowns(sys)`: The set of unknowns in the implicit discrete system. - - `get_ps(sys)` or `parameters(sys)`: The parameters of the implicit discrete system. - - `get_iv(sys)`: The independent variable of the implicit discrete system - -## Transformations - -```@docs; canonical=false -structural_simplify -``` - -## Problem Constructors - -```@docs; canonical=false -ImplicitDiscreteProblem(sys::ImplicitDiscreteSystem, u0map, tspan) -ImplicitDiscreteFunction(sys::ImplicitDiscreteSystem, args...) -``` - -## Discrete Domain - -```@docs; canonical=false -Shift -``` diff --git a/docs/src/systems/JumpSystem.md b/docs/src/systems/JumpSystem.md deleted file mode 100644 index 5bd0d50602..0000000000 --- a/docs/src/systems/JumpSystem.md +++ /dev/null @@ -1,33 +0,0 @@ -# JumpSystem - -## System Constructors - -```@docs -JumpSystem -``` - -## Composition and Accessor Functions - - - `get_eqs(sys)` or `equations(sys)`: The equations that define the jump system. - - `get_unknowns(sys)` or `unknowns(sys)`: The set of unknowns in the jump system. - - `get_ps(sys)` or `parameters(sys)`: The parameters of the jump system. - - `get_iv(sys)`: The independent variable of the jump system. - - `discrete_events(sys)`: The set of discrete events in the jump system. - -## Transformations - -```@docs; canonical=false -structural_simplify -``` - -## Analyses - -## Problem Constructors - -```@docs; canonical=false -DiscreteProblem(sys::JumpSystem, u0map, tspan) -``` - -```@docs -JumpProblem(sys::JumpSystem, prob, aggregator) -``` diff --git a/docs/src/systems/NonlinearSystem.md b/docs/src/systems/NonlinearSystem.md deleted file mode 100644 index 06d587b1b9..0000000000 --- a/docs/src/systems/NonlinearSystem.md +++ /dev/null @@ -1,57 +0,0 @@ -# NonlinearSystem - -## System Constructors - -```@docs -NonlinearSystem -``` - -## Composition and Accessor Functions - - - `get_eqs(sys)` or `equations(sys)`: The equations that define the nonlinear system. - - `get_unknowns(sys)` or `unknowns(sys)`: The set of unknowns in the nonlinear system. - - `get_ps(sys)` or `parameters(sys)`: The parameters of the nonlinear system. - - `get_u0_p(sys, u0map, parammap)` Numeric arrays for the initial condition and parameters given `var => value` maps. - -## Transformations - -```@docs; canonical=false -structural_simplify -alias_elimination -tearing -``` - -## Analyses - -```@docs; canonical=false -ModelingToolkit.isaffine -ModelingToolkit.islinear -``` - -## Applicable Calculation and Generation Functions - -```@docs; canonical=false -calculate_jacobian -generate_jacobian -jacobian_sparsity -``` - -## Problem Constructors - -```@docs -NonlinearFunction(sys::ModelingToolkit.NonlinearSystem, args...) -NonlinearProblem(sys::ModelingToolkit.NonlinearSystem, args...) -``` - -## Torn Problem Constructors - -```@docs -BlockNonlinearProblem -``` - -## Expression Constructors - -```@docs -NonlinearFunctionExpr -NonlinearProblemExpr -``` diff --git a/docs/src/systems/ODESystem.md b/docs/src/systems/ODESystem.md deleted file mode 100644 index 24e2952fc5..0000000000 --- a/docs/src/systems/ODESystem.md +++ /dev/null @@ -1,74 +0,0 @@ -# ODESystem - -## System Constructors - -```@docs -ODESystem -``` - -## Composition and Accessor Functions - - - `get_eqs(sys)` or `equations(sys)`: The equations that define the ODE. - - `get_unknowns(sys)` or `unknowns(sys)`: The set of unknowns in the ODE. - - `get_ps(sys)` or `parameters(sys)`: The parameters of the ODE. - - `get_iv(sys)`: The independent variable of the ODE. - - `get_u0_p(sys, u0map, parammap)` Numeric arrays for the initial condition and parameters given `var => value` maps. - - `continuous_events(sys)`: The set of continuous events in the ODE. - - `discrete_events(sys)`: The set of discrete events in the ODE. - - `alg_equations(sys)`: The algebraic equations (i.e. that does not contain a differential) that defines the ODE. - - `get_alg_eqs(sys)`: The algebraic equations (i.e. that does not contain a differential) that defines the ODE. Only returns equations of the current-level system. - - `diff_equations(sys)`: The differential equations (i.e. that contain a differential) that defines the ODE. - - `get_diff_eqs(sys)`: The differential equations (i.e. that contain a differential) that defines the ODE. Only returns equations of the current-level system. - - `has_alg_equations(sys)`: Returns `true` if the ODE contains any algebraic equations (i.e. that does not contain a differential). - - `has_alg_eqs(sys)`: Returns `true` if the ODE contains any algebraic equations (i.e. that does not contain a differential). Only considers the current-level system. - - `has_diff_equations(sys)`: Returns `true` if the ODE contains any differential equations (i.e. that does contain a differential). - - `has_diff_eqs(sys)`: Returns `true` if the ODE contains any differential equations (i.e. that does contain a differential). Only considers the current-level system. - -## Transformations - -```@docs -structural_simplify -ode_order_lowering -dae_index_lowering -change_independent_variable -liouville_transform -alias_elimination -tearing -``` - -## Analyses - -```@docs -ModelingToolkit.islinear -ModelingToolkit.isautonomous -ModelingToolkit.isaffine -``` - -## Applicable Calculation and Generation Functions - -```@docs; canonical=false -calculate_jacobian -calculate_tgrad -calculate_factorized_W -generate_jacobian -generate_tgrad -generate_factorized_W -jacobian_sparsity -``` - -## Standard Problem Constructors - -```@docs -ODEFunction(sys::ModelingToolkit.AbstractODESystem, args...) -ODEProblem(sys::ModelingToolkit.AbstractODESystem, args...) -SteadyStateProblem(sys::ModelingToolkit.AbstractODESystem, args...) -DAEProblem(sys::ModelingToolkit.AbstractODESystem, args...) -``` - -## Expression Constructors - -```@docs -ODEFunctionExpr -DAEFunctionExpr -SteadyStateProblemExpr -``` diff --git a/docs/src/systems/OptimizationSystem.md b/docs/src/systems/OptimizationSystem.md deleted file mode 100644 index bcc8b21de7..0000000000 --- a/docs/src/systems/OptimizationSystem.md +++ /dev/null @@ -1,40 +0,0 @@ -# OptimizationSystem - -## System Constructors - -```@docs -OptimizationSystem -``` - -## Composition and Accessor Functions - - - `get_op(sys)`: The objective to be minimized. - - `get_unknowns(sys)` or `unknowns(sys)`: The set of unknowns for the optimization. - - `get_ps(sys)` or `parameters(sys)`: The parameters for the optimization. - - `get_constraints(sys)` or `constraints(sys)`: The constraints for the optimization. - -## Transformations - -## Analyses - -## Applicable Calculation and Generation Functions - -```julia -calculate_gradient -calculate_hessian -generate_gradient -generate_hessian -hessian_sparsity -``` - -## Problem Constructors - -```@docs -OptimizationProblem(sys::ModelingToolkit.OptimizationSystem, args...) -``` - -## Expression Constructors - -```@docs -OptimizationProblemExpr -``` diff --git a/docs/src/systems/SDESystem.md b/docs/src/systems/SDESystem.md deleted file mode 100644 index 5789d2d9cb..0000000000 --- a/docs/src/systems/SDESystem.md +++ /dev/null @@ -1,70 +0,0 @@ -# SDESystem - -## System Constructors - -```@docs -SDESystem -``` - -To convert an `ODESystem` to an `SDESystem` directly: - -``` -ode = ODESystem(eqs,t,[x,y,z],[σ,ρ,β]) -sde = SDESystem(ode, noiseeqs) -``` - -## Composition and Accessor Functions - - - `get_eqs(sys)` or `equations(sys)`: The equations that define the SDE. - - `get_unknowns(sys)` or `unknowns(sys)`: The set of unknowns in the SDE. - - `get_ps(sys)` or `parameters(sys)`: The parameters of the SDE. - - `get_iv(sys)`: The independent variable of the SDE. - - `continuous_events(sys)`: The set of continuous events in the SDE. - - `discrete_events(sys)`: The set of discrete events in the SDE. - - `alg_equations(sys)`: The algebraic equations (i.e. that does not contain a differential) that defines the ODE. - - `get_alg_eqs(sys)`: The algebraic equations (i.e. that does not contain a differential) that defines the ODE. Only returns equations of the current-level system. - - `diff_equations(sys)`: The differential equations (i.e. that contain a differential) that defines the ODE. - - `get_diff_eqs(sys)`: The differential equations (i.e. that contain a differential) that defines the ODE. Only returns equations of the current-level system. - - `has_alg_equations(sys)`: Returns `true` if the ODE contains any algebraic equations (i.e. that does not contain a differential). - - `has_alg_eqs(sys)`: Returns `true` if the ODE contains any algebraic equations (i.e. that does not contain a differential). Only considers the current-level system. - - `has_diff_equations(sys)`: Returns `true` if the ODE contains any differential equations (i.e. that does contain a differential). - - `has_diff_eqs(sys)`: Returns `true` if the ODE contains any differential equations (i.e. that does contain a differential). Only considers the current-level system. - -## Transformations - -```@docs; canonical=false -structural_simplify -alias_elimination -``` - -```@docs -ModelingToolkit.Girsanov_transform -``` - -## Analyses - -## Applicable Calculation and Generation Functions - -```@docs; canonical=false -calculate_jacobian -calculate_tgrad -calculate_factorized_W -generate_jacobian -generate_tgrad -generate_factorized_W -jacobian_sparsity -``` - -## Problem Constructors - -```@docs -SDEFunction(sys::ModelingToolkit.SDESystem, args...) -SDEProblem(sys::ModelingToolkit.SDESystem, args...) -``` - -## Expression Constructors - -```@docs -SDEFunctionExpr -SDEProblemExpr -``` From e61981103e91c33a5af6d71343560ba90d811bd6 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Fri, 23 May 2025 15:30:37 +0530 Subject: [PATCH 05/49] docs: update more docs to new syntax --- docs/src/comparison.md | 4 ++-- docs/src/internals.md | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/src/comparison.md b/docs/src/comparison.md index 52d5ab2f70..6691b0023a 100644 --- a/docs/src/comparison.md +++ b/docs/src/comparison.md @@ -12,7 +12,7 @@ - All current Modelica compiler implementations are fixed and not extendable by the users from the Modelica language itself. For example, the Dymola compiler [shares its symbolic processing pipeline](https://www.claytex.com/tech-blog/model-translation-and-symbolic-manipulation/), - which is roughly equivalent to the `dae_index_lowering` and `structural_simplify` + which is roughly equivalent to the `dae_index_lowering` and `mtkcompile` of ModelingToolkit.jl. ModelingToolkit.jl is an open and hackable transformation system which allows users to add new non-standard transformations and control the order of application. @@ -90,7 +90,7 @@ [Dymola symbolic processing pipeline](https://www.claytex.com/tech-blog/model-translation-and-symbolic-manipulation/) with some improvements. ModelingToolkit.jl has an open transformation pipeline that allows for users to extend and reorder transformation passes, where - `structural_simplify` is an adaptation of the Modia.jl-improved alias elimination + `mtkcompile` is an adaptation of the Modia.jl-improved alias elimination and tearing algorithms. - Both Modia and ModelingToolkit generate `DAEProblem` and `ODEProblem` forms for solving with [DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/). diff --git a/docs/src/internals.md b/docs/src/internals.md index 00b29f1a64..ed83192f21 100644 --- a/docs/src/internals.md +++ b/docs/src/internals.md @@ -18,7 +18,7 @@ and are then used to generate the `observed` equation found in the variable when necessary. In this sense, there is an equivalence between observables and the variable elimination system. -The procedure for variable elimination inside [`structural_simplify`](@ref) is +The procedure for variable elimination inside [`mtkcompile`](@ref) is 1. [`ModelingToolkit.initialize_system_structure`](@ref). 2. [`ModelingToolkit.alias_elimination`](@ref). This step moves equations into `observed(sys)`. From b4335a7191e4a375abda92d8eeacc1c8624021f1 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Fri, 23 May 2025 15:30:53 +0530 Subject: [PATCH 06/49] docs: move `PDESystem` docs --- docs/src/{systems => API}/PDESystem.md | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename docs/src/{systems => API}/PDESystem.md (100%) diff --git a/docs/src/systems/PDESystem.md b/docs/src/API/PDESystem.md similarity index 100% rename from docs/src/systems/PDESystem.md rename to docs/src/API/PDESystem.md From 55ea020116bbdcb3e7df405056e73e2af2c12b6f Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Fri, 23 May 2025 15:31:03 +0530 Subject: [PATCH 07/49] feat: add new API docs --- docs/src/API/System.md | 174 +++++++++++++++++++++++ docs/src/API/model_building.md | 253 +++++++++++++++++++++++++++++++++ docs/src/API/problems.md | 121 ++++++++++++++++ 3 files changed, 548 insertions(+) create mode 100644 docs/src/API/System.md create mode 100644 docs/src/API/model_building.md create mode 100644 docs/src/API/problems.md diff --git a/docs/src/API/System.md b/docs/src/API/System.md new file mode 100644 index 0000000000..33a4d064ab --- /dev/null +++ b/docs/src/API/System.md @@ -0,0 +1,174 @@ +# [The `System` type](@id System_type) + +ModelingToolkit.jl uses `System` to symbolically represent all types of numerical problems. +Users create `System`s representing the problem they want to solve and `mtkcompile` transforms +them into a format ModelingToolkit.jl can generate code for (alongside performing other +optimizations). + +```@docs +System +``` + +## Utility constructors + +Several utility constructors also exist to easily construct alternative system formulations. + +```@docs +NonlinearSystem +SDESystem +JumpSystem +OptimizationSystem +``` + +## Accessor functions + +Several accessor functions exist to query systems for the information they contain. In general, +for every field `x` there exists a `has_x` function which checks if the system contains the +field and a `get_x` function for obtaining the value in the field. Note that fields of a system +cannot be accessed via `getproperty` - that is reserved for accessing variables, subsystems +or analysis points of the hierarchical system. + +```@docs +ModelingToolkit.has_eqs +ModelingToolkit.get_eqs +equations +equations_toplevel +full_equations +ModelingToolkit.has_noise_eqs +ModelingToolkit.get_noise_eqs +ModelingToolkit.has_jumps +ModelingToolkit.get_jumps +jumps +ModelingToolkit.has_constraints +ModelingToolkit.get_constraints +constraints +ModelingToolkit.has_costs +ModelingToolkit.get_costs +costs +ModelingToolkit.has_consolidate +ModelingToolkit.get_consolidate +ModelingToolkit.has_unknowns +ModelingToolkit.get_unknowns +unknowns +unknowns_toplevel +ModelingToolkit.has_ps +ModelingToolkit.get_ps +parameters +parameters_toplevel +tunable_parameters +ModelingToolkit.has_brownians +ModelingToolkit.get_brownians +brownians +ModelingToolkit.has_iv +ModelingToolkit.get_iv +ModelingToolkit.has_observed +ModelingToolkit.get_observed +observed +observables +ModelingToolkit.has_name +ModelingToolkit.get_name +nameof +ModelingToolkit.has_description +ModelingToolkit.get_description +ModelingToolkit.description +ModelingToolkit.has_defaults +ModelingToolkit.get_defaults +defaults +ModelingToolkit.has_guesses +ModelingToolkit.get_guesses +guesses +ModelingToolkit.has_initialization_eqs +ModelingToolkit.get_initialization_eqs +initialization_equations +ModelingToolkit.has_continuous_events +ModelingToolkit.get_continuous_events +continuous_events +continuous_events_toplevel +ModelingToolkit.has_discrete_events +ModelingToolkit.get_discrete_events +discrete_events_toplevel +ModelingToolkit.has_assertions +ModelingToolkit.get_assertions +assertions +ModelingToolkit.has_metadata +ModelingToolkit.get_metadata +SymbolicUtils.getmetadata(::ModelingToolkit.AbstractSystem, ::DataType, ::Any) +SymbolicUtils.setmetadata(::ModelingToolkit.AbstractSystem, ::DataType, ::Any) +ModelingToolkit.has_is_dde +ModelingToolkit.get_is_dde +ModelingToolkit.is_dde +ModelingToolkit.has_tstops +ModelingToolkit.get_tstops +ModelingToolkit.symbolic_tstops +ModelingToolkit.has_tearing_state +ModelingToolkit.get_tearing_state +ModelingToolkit.does_namespacing +toggle_namespacing +ModelingToolkit.iscomplete +ModelingToolkit.has_preface +ModelingToolkit.get_preface +ModelingToolkit.preface +ModelingToolkit.has_parent +ModelingToolkit.get_parent +ModelingToolkit.has_initializesystem +ModelingToolkit.get_initializesystem +ModelingToolkit.is_initializesystem +``` + +## `getproperty` syntax + +ModelingToolkit allows obtaining in a system using `getproperty`. For a system `sys` with a +subcomponent `inner` containing variable `var`, `sys.inner.var` will obtain the appropriately +namespaced version of `var`. Note that this can also be used to access subsystems (`sys.inner`) +or analysis points. + +!!! note + + By default, top-level systems not marked as `complete` will apply their namespace. Systems + marked as `complete` will not do this namespacing. This namespacing behavior can be toggled + independently of whether the system is completed using [`toggle_namespacing`](@ref) and the + current namespacing behavior can be queried via [`ModelingToolkit.does_namespacing`](@ref). + +```@docs +Base.getproperty(::ModelingToolkit.AbstractSystem, ::Symbol) +``` + +## Functions for querying system equations + +```@docs +has_diff_eqs +has_alg_eqs +get_diff_eqs +get_alg_eqs +has_diff_equations +has_alg_equations +diff_equations +alg_equations +is_alg_equation +is_diff_equation +``` + +## String parsing + +ModelingToolkit can parse system variables from strings. + +```@docs +ModelingToolkit.parse_variable +``` + +## Dumping system data + +```@docs +ModelingToolkit.dump_unknowns +ModelingToolkit.dump_parameters +``` + +```@docs; canonical = false +ModelingToolkit.dump_variable_metadata +``` + +## Debugging utilities + +```@docs +debug_system +``` diff --git a/docs/src/API/model_building.md b/docs/src/API/model_building.md new file mode 100644 index 0000000000..40072bc2c4 --- /dev/null +++ b/docs/src/API/model_building.md @@ -0,0 +1,253 @@ +# [Model building reference](@id model_building_api) + +This page lists functionality and utilities related to building hierarchical models. It is +recommended to read the page on the [`System`](@ref System_type) before this. + +## Hierarchical model composition + +The `System` data structure can represent a tree-like hierarchy of systems for building models +from composable blocks. The [`ModelingToolkit.get_systems`](@ref) function can be used for +querying the subsystems of a system. The `@component` macro should be used when writing +building blocks for model composition. + +```@docs +@component +``` + +Every constructor function should build either a component or a connector. Components define +the dynamics of the system. Connectors are used to connect components together and propagate +information between them. See also [`@connector`](@ref). + +### Scoping of variables + +When building hierarchical systems, is is often necessary to pass variables from a parent system +to the subsystems. If done naively, this will result in the child system assuming it "owns" the +variables passed to it and any occurrences of those variables in the child system will be +namespaced. To prevent this, ModelingToolkit has the concept of variable scope. The scope allows +specifying which system a variable belongs to relative to the system in which it is used. + +```@docs +LocalScope +ParentScope +GlobalScope +``` + +Note that the scopes must be applied to _individual variables_ and not expressions. For example, +`ParentScope(x + y)` is incorrect. Instead, `ParentScope(x) + ParentScope(y)` is the correct usage. +Applying the same scope (more generally, the same function) to all variables in an expression is a +common task, and ModelingToolkit exposes a utility for the same: + +```@docs +ModelingToolkit.apply_to_variables +``` + +It is still tedious to manually use `apply_to_variables` on any symbolic expression passed to a +subsystem. The `@named` macro automatically wraps all symbolic arguments in `ParentScope` and +uses the identifier being assigned as the name of the system. + +```@docs +@named +``` + +### Exploring the tree structure + +The `System` type implements the `AbstractTrees` interface. This can be used to explore the +hierarchical structure. + +```@docs +hierarchy +``` + +### [Connection semantics](@id connect_semantics) + +ModelingToolkit implements connection semantics similar to those in the [Modelica specification](https://specification.modelica.org/maint/3.6/connectors-and-connections.html). +We do not support the concept of `inner` and `outer` elements or `expandable` connectors. +Connectors in ModelingToolkit are systems with the appropriate metadata added via the `@connector` +macro. + +```@docs +connect +domain_connect +@connector +``` + +Connections can be expanded using `expand_connections`. + +```@docs +expand_connections +``` + +Similar to the `stream` and `flow` keyword arguments in the specification, ModelingToolkit +allows specifying how variables in a connector behave in a connection. + +```@docs +Equality +Flow +Stream +``` + +These are specified using the `connect` metadata. ModelingToolkit also supports `instream`. +Refer to the Modelica specification on [Stream connectors](https://specification.modelica.org/maint/3.6/stream-connectors.html) +for more information. + +```@docs +instream +``` + +### System composition utilities + +```@docs +extend +compose +substitute_component +``` + +### Flattening systems + +The hierarchical structure can be flattened. This operation is performed during simplification. + +```@docs +flatten +``` + +## System simplification + +`System`s can be simplified to reformulate them in a way that enables it to be solved numerically, +and also perform other optimizations. This is done via the `mtkcompile` function. Connection expansion +and flattening are preprocessing steps of simplification. + +```@docs +mtkcompile +@mtkcompile +``` + +It is also possible (though not always advisable) to build numerical problems from systems without +passing them through `mtkcompile`. To do this, the system must first be marked as "complete" via +the `complete` function. This process is used to indicate that a system will not be modified +further and allows ModelingToolkit to perform any necessary preprocessing to it. `mtkcompile` +calls `complete` internally. + +```@docs +complete +``` + +### Exploring the results of simplification + +Similar to how [`full_equations`](@ref) returns the equations of a system with all variables +eliminated during `mtkcompile` substituted, we can perform this substitution on an arbitrary +expression. + +```@docs +ModelingToolkit.substitute_observed +ModelingToolkit.empty_substitutions +ModelingToolkit.get_substitutions +``` + +### Experimental simplification + +ModelingToolkit may have a variety of experimental simplification passes. These are not +enabled by default, but can be used by passing to the `additional_passes` keyword argument +of `mtkcompile`. + +```@docs +ModelingToolkit.IfLifting +``` + +## Event handling + +Time-dependent systems may have several events. These are used to trigger discontinuities +in the model. They compile to standard callbacks from `DiffEqCallbacks.jl`. + +```@docs +ModelingToolkit.SymbolicContinuousCallback +ModelingToolkit.SymbolicDiscreteCallback +``` + +The affect functions for the above callbacks can be symbolic or user-defined functions. +Symbolic affects are handled using equations as described in the [Events](@ref events) +section of the documentation. User-defined functions can be used via `ImperativeAffect`. + +```@docs +ModelingToolkit.AffectSystem +ModelingToolkit.ImperativeAffect +``` + +## Modelingtoolkitize + +ModelingToolkit can take some numerical problems created non-symbolically and build a +symbolic representation from them. + +```@docs +modelingtoolkitize +``` + +## Using FMUs + +ModelingToolkit is capable of importing FMUs as black-box symbolic models. Currently only +a subset of FMU features are supported. This functionality requires importing `FMI.jl`. + +```@docs +ModelingToolkit.FMIComponent +``` + +## Model transformations + +ModelingToolkit exposes a variety of transformations that can be applied to models to aid in +symbolic analysis. + +```@docs +liouville_transform +stochastic_integral_transform +Girsanov_transform +change_independent_variable +add_accumulations +noise_to_brownians +convert_system_indepvar +``` + +## Hybrid systems + +Hybrid systems are dynamical systems involving one or more discrete-time subsystems. These +discrete time systems follow clock semantics - they are synchronous systems and the relevant +variables are only defined at points where the clock ticks. + +While ModelingToolkit is unable to simplify, compile and solve such systems on its own, it +has the ability to represent them. Compilation strategies can be implemented independently +on top of [`mtkcompile`](@ref) using the `additional_passes` functionality. + +!!! warn + + These operators are considered experimental API. + +```@docs; canonical = false +Sample +Hold +SampleTime +``` + +ModelingToolkit uses the clock definition in SciMLBase + +```@docs +SciMLBase.TimeDomain +SciMLBase.Clock +SciMLBase.SolverStepClock +SciMLBase.Continuous +``` + +### State machines + +While ModelingToolkit has the capability to represent state machines, it lacks the ability +to compile and simulate them. + +!!! warn + + This functionality is considered experimental API + +```@docs +initial_state +transition +activeState +entry +ticksInState +timeInState +``` diff --git a/docs/src/API/problems.md b/docs/src/API/problems.md new file mode 100644 index 0000000000..5549224a09 --- /dev/null +++ b/docs/src/API/problems.md @@ -0,0 +1,121 @@ +# Building and solving numerical problems + +Systems are numerically solved by building and solving the appropriate problem type. +Numerical solvers expect to receive functions taking a predefeined set of arguments +and returning specific values. This format of argument and return value depends on +the function and the problem. ModelingToolkit is capable of compiling and generating +code for a variety of such numerical problems. + +## Dynamical systems + +```@docs +ODEFunction(::System, args...) +ODEProblem(::System, args...) +DAEFunction(::System, args...) +DAEProblem(::System, args...) +SDEFunction(::System, args...) +SDEProblem(::System, args...) +DDEFunction(::System, args...) +DDEProblem(::System, args...) +SDDEFunction(::System, args...) +SDDEProblem(::System, args...) +JumpProblem(::System, args...) +BVProblem(::System, args...) +DiscreteProblem(::System, args...) +ImplicitDiscreteProblem(::System, args...) +``` + +## Nonlinear systems + +```@docs +NonlinearFunction(::System, args...) +NonlinearProblem(::System, args...) +SCCNonlinearProblem(::System, args...) +NonlinearLeastSquaresProblem(::System, args...) +SteadyStateProblem(::System, args...) +IntervalNonlinearFunction(::System, args...) +IntervalNonlinearProblem(::System, args...) +ModelingToolkit.HomotopyContinuationProblem +HomotopyNonlinearFunction(::System, args...) +``` + +## Optimization and optimal control + +```@docs +OptimizationFunction(::System, args...) +OptimizationProblem(::System, args...) +ODEInputFunction(::System, args...) +JuMPDynamicOptProblem(::System, args...) +InfiniteOptDynamicOptProblem,(::System, args...) +CasADiDynamicOptProblem(::System, args...) +DynamicOptSolution +``` + +## The state vector and parameter object + +Typically the unknowns of the system are present as a `Vector` of the appropriate length +in the numerical problem. The state vector can also be constructed manually without building +a problem. + +```@docs +ModelingToolkit.get_u0 +ModelingToolkit.varmap_to_vars +``` + +By default, the parameters of the system are stored in a custom data structure called +`MTKParameters`. The internals of this data structure are undocumented, and it should +only be interacted with through defined public API. SymbolicIndexingInterface.jl contains +functionality useful for this purpose. + +```@docs +MTKParameters +ModelingToolkit.get_p +``` + +The following functions are useful when working with `MTKParameters` objects, and especially +the `Tunables` portion. For more information about the "portions" of `MTKParameters`, refer +to the [`SciMLStructures.jl`](https://docs.sciml.ai/SciMLStructures/stable/) documentation. + +```@docs +reorder_dimension_by_tunables! +reorder_dimension_by_tunables +``` + +## Initialization + +```@docs +generate_initializesystem +InitializationProblem +``` + +## Linear analysis + +```@docs +linearization_function +LinearizationProblem +linearize +CommonSolve.solve(::LinearizationProblem) +linearize_symbolic +``` + +There are also utilities for manipulating the results of these analyses in a symbolic context. + +```@docs +ModelingToolkit.similarity_transform +ModelingToolkit.reorder_unknnowns +``` + +### Analysis point transformations + +Linear analysis can also be done using analysis points to perform several common +workflows. + +```@docs +get_sensitivity_function +get_sensitivity +get_comp_sensitivity_function +get_comp_sensitivity +get_looptransfer_function +get_looptransfer +open_loop +``` From b84017c7835fb308ad0d06e7a74badc6d5c37f91 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Fri, 23 May 2025 15:31:28 +0530 Subject: [PATCH 08/49] docs: document `getproperty(::AbstractSystem, ::Symbol)` --- src/systems/abstractsystem.jl | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl index 41cc96bf57..b997219eba 100644 --- a/src/systems/abstractsystem.jl +++ b/src/systems/abstractsystem.jl @@ -955,6 +955,14 @@ function Base.propertynames(sys::AbstractSystem; private = false) end end +""" + Base.getproperty(sys::AbstractSystem, name::Symbol) + +Access the subsystem, variable or analysis point of `sys` named `name`. To check if `sys` +will namespace the returned value, use `ModelingToolkit.does_namespacing(sys)`. + +See also: [`ModelingToolkit.does_namespacing`](@ref). +""" function Base.getproperty( sys::AbstractSystem, name::Symbol; namespace = does_namespacing(sys)) if has_parent(sys) && (parent = get_parent(sys); parent !== nothing) From b85a193e49863365b451c85392bd248c6fc54c1c Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Fri, 23 May 2025 15:32:50 +0530 Subject: [PATCH 09/49] docs: add `INTERNAL_FIELD_WARNING` and `INTERNAL_ARGS_WARNING` --- src/ModelingToolkit.jl | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/src/ModelingToolkit.jl b/src/ModelingToolkit.jl index 4b0286f218..9bb96e1e14 100644 --- a/src/ModelingToolkit.jl +++ b/src/ModelingToolkit.jl @@ -122,6 +122,16 @@ for fun in [:toexpr] end end +const INTERNAL_FIELD_WARNING = """ +This field is internal API. It may be removed or changed without notice in a non-breaking \ +release. Usage of this field is not advised. +""" + +const INTERNAL_ARGS_WARNING = """ +The following arguments are internal API. They may be removed or changed without notice \ +in a non-breaking release. Usage of these arguments is not advised. +""" + """ $(TYPEDEF) From d5dab58bab11ab7aaa62753b2dab3f77992fb5a1 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Fri, 23 May 2025 15:31:49 +0530 Subject: [PATCH 10/49] docs: document `System` and its fields --- src/systems/system.jl | 187 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 187 insertions(+) diff --git a/src/systems/system.jl b/src/systems/system.jl index 5da5d9db79..f327a9ca9a 100644 --- a/src/systems/system.jl +++ b/src/systems/system.jl @@ -11,50 +11,237 @@ end const MetadataT = Base.ImmutableDict{DataType, Any} +""" + $(TYPEDEF) + +A symbolic representation of a numerical system to be solved. This is a recursive +tree-like data structure - each system can contain additional subsystems. As such, +it implements the `AbstractTrees.jl` interface to enable exploring the hierarchical +structure. + +# Fields + +$(TYPEDFIELDS) +""" struct System <: AbstractSystem + """ + $INTERNAL_FIELD_WARNING + A unique integer tag for the system. + """ tag::UInt + """ + The equations of the system. + """ eqs::Vector{Equation} # nothing - no noise # vector - diagonal noise # matrix - generic form # column matrix - scalar noise + """ + The noise terms for each equation of the system. This field is only used for flattened + systems. To represent noise in a hierarchical system, use brownians. In a system with + `N` equations and `K` independent brownian variables, this should be an `N x K` + matrix. In the special case where `N == K` and each equation has independent noise, + this noise matrix is diagonal. Diagonal noise can be specified by providing an `N` + length vector. If this field is `nothing`, the system does not have noise. + """ noise_eqs::Union{Nothing, AbstractVector, AbstractMatrix} + """ + Jumps associated with the system. Each jump can be a `VariableRateJump`, + `ConstantRateJump` or `MassActionJump`. See `JumpProcesses.jl` for more information. + """ jumps::Vector{JumpType} + """ + The constraints of the system. This can be used to represent the constraints in an + optimal-control problem or boundary-value differential equation, or the constraints + in a constrained optimization. + """ constraints::Vector{Union{Equation, Inequality}} + """ + The costs of the system. This can be the cost in an optimal-control problem, or the + loss of an optimization problem. Scalar loss values must also be provided as a single- + element vector. + """ costs::Vector{<:Union{BasicSymbolic, Real}} + """ + A function which combines costs into a scalar value. This should take two arguments, + the `costs` of this system and the consolidated costs of all subsystems in the order + they are present in the `systems` field. It should return a scalar cost that combines + all of the individual values. This defaults to a function that simply sums all cost + values. + """ consolidate::Any + """ + The variables being solved for by this system. For example, in a differential equation + system, this contains the dependent variables. + """ unknowns::Vector + """ + The parameters of the system. Parameters can either be variables that parameterize the + problem being solved for (e.g. the spring constant of a mass-spring system) or + additional unknowns not part of the main dynamics of the system (e.g. discrete/clocked + variables in a hybrid ODE). + """ ps::Vector + """ + The brownian variables of the system, created via `@brownians`. Each brownian variable + represents an independent noise. A system with brownians cannot be simulated directly. + It needs to be compiled using `mtkcompile` into `noise_eqs`. + """ brownians::Vector + """ + The independent variable for a time-dependent system, or `nothing` for a time-independent + system. + """ iv::Union{Nothing, BasicSymbolic{Real}} + """ + Equations that compute variables of a system that have been eliminated from the set of + unknowns by `mtkcompile`. More generally, this contains all variables that can be + computed from the unknowns and parameters and do not need to be solved for. Such + variables are termed as "observables". Each equation must be of the form + `observable ~ expression` and observables cannot appear on the LHS of multiple + equations. Equations must be sorted such that every observable appears on + the left hand side of an equation before it appears on the right hand side of any other + equation. + """ observed::Vector{Equation} parameter_dependencies::Vector{Equation} + """ + $INTERNAL_FIELD_WARNING + A mapping from the name of a variable to the actual symbolic variable in the system. + This is used to enable `getproperty` syntax to access variables of a system. + """ var_to_name::Dict{Symbol, Any} + """ + The name of the system. + """ name::Symbol + """ + An optional description for the system. + """ description::String + """ + Default values that variables (unknowns/observables/parameters) should take when + constructing a numerical problem from the system. These values can be overridden + by initial values provided to the problem constructor. Defaults of parent systems + take priority over those in child systems. + """ defaults::Dict + """ + Guess values for variables of a system that are solved for during initialization. + """ guesses::Dict + """ + A list of subsystems of this system. Used for hierarchically building models. + """ systems::Vector{System} + """ + Equations that must be satisfied during initialization of the numerical problem created + from this system. For time-dependent systems, these equations are not valid after the + initial time. + """ initialization_eqs::Vector{Equation} + """ + Symbolic representation of continuous events in a dynamical system. See + [`SymbolicContinuousCallback`](@ref). + """ continuous_events::Vector{SymbolicContinuousCallback} + """ + Symbolic representation of discrete events in a dynamica system. See + [`SymbolicDiscreteCallback`](@ref). + """ discrete_events::Vector{SymbolicDiscreteCallback} + """ + $INTERNAL_FIELD_WARNING + If this system is a connector, the type of connector it is. + """ connector_type::Any + """ + A map from expressions that must be through throughout the solution process to an + associated error message. By default these assertions cause the generated code to + output `NaN`s if violated, but can be made to error using `debug_system`. + """ assertions::Dict{BasicSymbolic, String} + """ + The metadata associated with this system, as a `Base.ImmutableDict`. This follows + the same interface as SymbolicUtils.jl. Metadata can be queried and updated using + `SymbolicUtils.getmetadata` and `SymbolicUtils.setmetadata` respectively. + """ metadata::MetadataT + """ + $INTERNAL_FIELD_WARNING + Metadata added by the `@mtkmodel` macro. + """ gui_metadata::Any # ? + """ + Whether the system contains delay terms. This is inferred from the equations, but + can also be provided explicitly. + """ is_dde::Bool + """ + Extra time points for the integrator to stop at. These can be numeric values, + or expressions of parameters and time. + """ tstops::Vector{Any} + """ + The `TearingState` of the system post-simplification with `mtkcompile`. + """ tearing_state::Any + """ + Whether the system namespaces variables accessed via `getproperty`. `complete`d systems + do not namespace, but this flag can be toggled independently of `complete` using + `toggle_namespacing`. + """ namespacing::Bool + """ + Whether the system is marked as "complete". Completed systems cannot be used as + subsystems. + """ complete::Bool + """ + $INTERNAL_FIELD_WARNING + For systems simplified or completed with `split = true` (the default) this contains an + `IndexCache` which aids in symbolic indexing. If this field is `nothing`, the system is + either not completed, or completed with `split = false`. + """ index_cache::Union{Nothing, IndexCache} + """ + $INTERNAL_FIELD_WARNING + Connections that should be ignored because they were removed by an analysis point + transformation. The first element of the tuple contains all such "standard" connections + (ones between connector systems) and the second contains all such causal variable + connections. + """ ignored_connections::Union{ Nothing, Tuple{Vector{IgnoredAnalysisPoint}, Vector{IgnoredAnalysisPoint}}} + """ + `SymbolicUtils.Code.Assignment`s to prepend to all code generated from this system. + """ preface::Any + """ + After simplification with `mtkcompile`, this field contains the unsimplified system + with the hierarchical structure. There may be multiple levels of `parent`s. The root + parent is used for accessing variables via `getproperty` syntax. + """ parent::Union{Nothing, System} + """ + A custom initialization system to use if no initial conditions are provided for the + unknowns or observables of this system. + """ initializesystem::Union{Nothing, System} + """ + Whether the current system is an initialization system. + """ is_initializesystem::Bool + """ + $INTERNAL_FIELD_WARNING + Whether the system has been simplified by `mtkcompile`. + """ isscheduled::Bool + """ + $INTERNAL_FIELD_WARNING + The `Schedule` containing additional information about the simplified system. + """ schedule::Union{Schedule, Nothing} function System( From 3cf5248b7d63323852b403d173d4b57533dc848a Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Fri, 23 May 2025 15:32:01 +0530 Subject: [PATCH 11/49] docs: document `System` constructors --- src/systems/system.jl | 41 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 41 insertions(+) diff --git a/src/systems/system.jl b/src/systems/system.jl index f327a9ca9a..4909f5a9d3 100644 --- a/src/systems/system.jl +++ b/src/systems/system.jl @@ -301,6 +301,22 @@ function default_consolidate(costs, subcosts) return reduce(+, costs; init = 0.0) + reduce(+, subcosts; init = 0.0) end +""" + $(TYPEDSIGNATURES) + +Construct a system using the given equations `eqs`, independent variable `iv` (`nothing`) +for time-independent systems, unknowns `dvs`, parameters `ps` and brownian variables +`brownians`. + +## Keyword Arguments + +- `discover_from_metadata`: Whether to parse metadata of unknowns and parameters of the + system to obtain defaults and/or guesses. +- `checks`: Whether to perform sanity checks on the passed values. + +All other keyword arguments are named identically to the corresponding fields in +[`System`](@ref). +""" function System(eqs::Vector{Equation}, iv, dvs, ps, brownians = []; constraints = Union{Equation, Inequality}[], noise_eqs = nothing, jumps = [], costs = BasicSymbolic[], consolidate = default_consolidate, @@ -394,10 +410,23 @@ function System(eqs::Vector{Equation}, iv, dvs, ps, brownians = []; initializesystem, is_initializesystem; checks) end +""" + $(TYPEDSIGNATURES) + +Create a time-independent [`System`](@ref) with the given equations `eqs`, unknowns `dvs` +and parameters `ps`. +""" function System(eqs::Vector{Equation}, dvs, ps; kwargs...) System(eqs, nothing, dvs, ps; kwargs...) end +""" + $(TYPEDSIGNATURES) + +Create a time-dependent system with the given equations `eqs` and independent variable `iv`. +Discover variables, parameters and brownians in the system by parsing the equations and +other symbolic expressions passed to the system. +""" function System(eqs::Vector{Equation}, iv; kwargs...) iv === nothing && return System(eqs; kwargs...) @@ -486,6 +515,13 @@ function System(eqs::Vector{Equation}, iv; kwargs...) eqs, iv, collect(allunknowns), collect(new_ps), collect(brownians); kwargs...) end +""" + $(TYPEDSIGNATURES) + +Create a time-independent system with the given equations `eqs`. Discover variables and +parameters in the system by parsing the equations and other symbolic expressions passed to +the system. +""" function System(eqs::Vector{Equation}; kwargs...) eqs = collect(eqs) @@ -516,6 +552,11 @@ function System(eqs::Vector{Equation}; kwargs...) return System(eqs, nothing, collect(allunknowns), collect(new_ps); kwargs...) end +""" + $(TYPEDSIGNATURES) + +Create a `System` with a single equation `eq`. +""" System(eq::Equation, args...; kwargs...) = System([eq], args...; kwargs...) function gather_array_params(ps) From f0075c105bebea53d83ed0bf20da541fd384d8bc Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Fri, 23 May 2025 15:32:25 +0530 Subject: [PATCH 12/49] docs: document `flatten`, `getmetadata` and `setmetadata` --- src/systems/system.jl | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/src/systems/system.jl b/src/systems/system.jl index 4909f5a9d3..9c48375e86 100644 --- a/src/systems/system.jl +++ b/src/systems/system.jl @@ -693,6 +693,12 @@ function _check_if_dde(eqs, iv, subsystems) return is_dde end +""" + $(TYPEDSIGNATURES) + +Flatten the hierarchical structure of a system, collecting all equations, unknowns, etc. +into one top-level system after namespacing appropriately. +""" function flatten(sys::System, noeqs = false) systems = get_systems(sys) isempty(systems) && return sys @@ -863,11 +869,23 @@ function Base.hash(sys::System, h::UInt) return h end +""" + $(TYPEDSIGNATURES) + +Get the metadata associated with key `k` in system `sys` or `default` if it does not exist. +""" function SymbolicUtils.getmetadata(sys::AbstractSystem, k::DataType, default) meta = get_metadata(sys) return get(meta, k, default) end +""" + $(TYPEDSIGNATURES) + +Set the metadata associated with key `k` in system `sys` to value `v`. This is an +out-of-place operation, and will return a shallow copy of `sys` with the appropriate +metadata values. +""" function SymbolicUtils.setmetadata(sys::AbstractSystem, k::DataType, v) meta = get_metadata(sys) meta = Base.ImmutableDict(meta, k => v)::MetadataT From dd432517867891437665eb8022adc36e47c47b0b Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Fri, 23 May 2025 15:32:37 +0530 Subject: [PATCH 13/49] docs: document `XSystem` utility constructors --- src/systems/system.jl | 22 ++++++++++++++++++++++ 1 file changed, 22 insertions(+) diff --git a/src/systems/system.jl b/src/systems/system.jl index 9c48375e86..872cf5fc95 100644 --- a/src/systems/system.jl +++ b/src/systems/system.jl @@ -899,6 +899,13 @@ function check_complete(sys::System, obj) iscomplete(sys) || throw(SystemNotCompleteError(obj)) end +""" + $(TYPEDSIGNATURES) + +Convert a time-dependent system `sys` to a time-independent system of nonlinear +equations that solve for the steady state of the system where `D(x)` is zero for +each continuous variable `x`. +""" function NonlinearSystem(sys::System) if !is_time_dependent(sys) throw(ArgumentError("`NonlinearSystem` constructor expects a time-dependent `System`")) @@ -924,6 +931,11 @@ end # Utility constructors ######## +""" + $(METHODLIST) + +Construct a [`System`](@ref) to solve an optimization problem with the given scalar cost. +""" function OptimizationSystem(cost; kwargs...) return System(Equation[]; costs = [cost], kwargs...) end @@ -940,6 +952,11 @@ function OptimizationSystem(cost::Array, dvs, ps; kwargs...) return System(Equation[], nothing, dvs, ps; costs = vec(cost), kwargs...) end +""" + $(METHODLIST) + +Construct a [`System`](@ref) to solve a system of jump equations. +""" function JumpSystem(jumps, iv; kwargs...) mask = isa.(jumps, Equation) eqs = Vector{Equation}(jumps[mask]) @@ -954,6 +971,11 @@ function JumpSystem(jumps, iv, dvs, ps; kwargs...) return System(eqs, iv, dvs, ps; jumps, kwargs...) end +""" + $(METHODLIST) + +Construct a system of equations with associated noise terms. +""" function SDESystem(eqs::Vector{Equation}, noise, iv; is_scalar_noise = false, kwargs...) if is_scalar_noise if !(noise isa Vector) From 9ff28e46c412d6688c36095d6b2176ccfc26cad3 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Fri, 23 May 2025 17:23:27 +0530 Subject: [PATCH 14/49] docs: move variable metadata docs to `API` section --- .../Variable_metadata.md => API/variables.md} | 218 +++++++++++++----- 1 file changed, 154 insertions(+), 64 deletions(-) rename docs/src/{basics/Variable_metadata.md => API/variables.md} (63%) diff --git a/docs/src/basics/Variable_metadata.md b/docs/src/API/variables.md similarity index 63% rename from docs/src/basics/Variable_metadata.md rename to docs/src/API/variables.md index 43d0a7f2ea..0ff2e7799d 100644 --- a/docs/src/basics/Variable_metadata.md +++ b/docs/src/API/variables.md @@ -1,8 +1,22 @@ -# [Symbolic Metadata](@id symbolic_metadata) +# Symbolic variables and variable metadata -It is possible to add metadata to symbolic variables, the metadata will be displayed when calling help on a variable. +ModelingToolkit uses [Symbolics.jl](https://docs.sciml.ai/Symbolics/stable/) for the symbolic +manipulation infrastructure. In fact, the `@variables` macro is defined in Symbolics.jl. In +addition to `@variables`, ModelingToolkit defines `@parameters`, `@independent_variables`, +`@constants` and `@brownians`. These macros function identically to `@variables` but allow +ModelingToolkit to attach additional metadata. -The following information can be added (note, it's possible to extend this to user-defined metadata as well) +```@docs +Symbolics.@variables +@independent_variables +@parameters +@constants +@brownians +``` + +Symbolic variables can have metadata attached to them. The defaults and guesses assigned +at variable construction time are examples of this metadata. ModelingToolkit also defines +additional types of metadata. ## Variable descriptions @@ -39,6 +53,11 @@ help?> u Symbolics.VariableSource: (:variables, :u) ``` +```@docs +hasdescription +getdescription +``` + ## Connect Variables in connectors can have `connect` metadata which describes the type of connections. @@ -61,6 +80,16 @@ hasconnect(i) getconnect(k) ``` +```@docs +hasconnect +getconnect +``` + +```@docs; canonical = false +Flow +Stream +``` + ## Input or output Designate a variable as either an input or an output using the following @@ -78,16 +107,20 @@ isinput(u) isoutput(y) ``` +```@docs +isinput +isoutput +ModelingToolkit.setinput +ModelingToolkit.setoutput +``` + ## Bounds Bounds are useful when parameters are to be optimized, or to express intervals of uncertainty. -```@example metadata -@variables u [bounds = (-1, 1)] +```@repl metadata +@variables u [bounds = (-1, 1)]; hasbounds(u) -``` - -```@example metadata getbounds(u) ``` @@ -95,51 +128,45 @@ Bounds can also be specified for array variables. A scalar array bound is applie element of the array. A bound may also be specified as an array, in which case the size of the array must match the size of the symbolic variable. -```@example metadata -@variables x[1:2, 1:2] [bounds = (-1, 1)] +```@repl metadata +@variables x[1:2, 1:2] [bounds = (-1, 1)]; hasbounds(x) -``` - -```@example metadata getbounds(x) -``` - -```@example metadata getbounds(x[1, 1]) -``` - -```@example metadata getbounds(x[1:2, 1]) -``` - -```@example metadata -@variables x[1:2] [bounds = (-Inf, [1.0, Inf])] +@variables x[1:2] [bounds = (-Inf, [1.0, Inf])]; hasbounds(x) -``` - -```@example metadata getbounds(x) -``` - -```@example metadata getbounds(x[2]) +hasbounds(x[2]) ``` -```@example metadata -hasbounds(x[2]) +```@docs +hasbounds +getbounds ``` ## Guess -Specify an initial guess for custom initial conditions of an `ODESystem`. +Specify an initial guess for variables of a `System`. This is used when building the +[`InitializationProblem`](@ref). -```@example metadata -@variables u [guess = 1] +```@repl metadata +@variables u [guess = 1]; hasguess(u) +getguess(u) ``` -```@example metadata -getguess(u) +```@docs +hasguess +getguess +``` + +When a system is constructed, the guesses of the involved variables are stored in a `Dict` +in the system. After this point, the guess metadata of the variable is irrelevant. + +```@docs; canonical = false +guesses ``` ## Mark input as a disturbance @@ -151,6 +178,10 @@ Indicate that an input is not available for control, i.e., it's a disturbance in isdisturbance(u) ``` +```@docs +isdisturbance +``` + ## Mark parameter as tunable Indicate that a parameter can be automatically tuned by parameter optimization or automatic control tuning apps. @@ -160,20 +191,32 @@ Indicate that a parameter can be automatically tuned by parameter optimization o istunable(Kp) ``` +```@docs +istunable +ModelingToolkit.isconstant +``` + +!!! note + + [`@constants`](@ref) is a convenient way to create `@parameters` with `tunable = false` + metadata + ## Probability distributions A probability distribution may be associated with a parameter to indicate either uncertainty about its value, or as a prior distribution for Bayesian optimization. -```julia -using Distributions -d = Normal(10, 1) -@parameters m [dist = d] +```@repl metadata +using Distributions; +d = Normal(10, 1); +@parameters m [dist = d]; hasdist(m) +getdist(m) ``` -```julia -getdist(m) +```@docs +hasdist +getdist ``` ## Irreducible @@ -187,6 +230,10 @@ it can be accessed in [callbacks](@ref events) isirreducible(important_value) ``` +```@docs +isirreducible +``` + ## State Priority When a model is structurally simplified, the algorithm will try to ensure that the variables with higher state priority become states of the system. A variable's state priority is a number set using the `state_priority` metadata. @@ -196,18 +243,24 @@ When a model is structurally simplified, the algorithm will try to ensure that t state_priority(important_dof) ``` +```@docs +state_priority +``` + ## Units Units for variables can be designated using symbolic metadata. For more information, please see the [model validation and units](@ref units) section of the docs. Note that `getunit` is not equivalent to `get_unit` - the former is a metadata getter for individual variables (and is provided so the same interface function for `unit` exists like other metadata), while the latter is used to handle more general symbolic expressions. -```@example metadata -using DynamicQuantities -@variables speed [unit = u"m/s"] +```@repl metadata +using DynamicQuantities; +@variables speed [unit = u"m/s"]; hasunit(speed) +getunit(speed) ``` -```@example metadata -getunit(speed) +```@docs +hasunit +getunit ``` ## Miscellaneous metadata @@ -215,13 +268,23 @@ getunit(speed) User-defined metadata can be added using the `misc` metadata. This can be queried using the `hasmisc` and `getmisc` functions. -```@example metadata -@variables u [misc = :conserved_parameter] y [misc = [2, 4, 6]] +```@repl metadata +@variables u [misc = :conserved_parameter] y [misc = [2, 4, 6]]; hasmisc(u) +getmisc(y) ``` -```@example metadata -getmisc(y) +```@docs +hasmisc +getmisc +``` + +## Dumping metadata + +ModelingToolkit allows dumping the metadata of a variable as a `NamedTuple`. + +```@docs +ModelingToolkit.dump_variable_metadata ``` ## Additional functions @@ -250,25 +313,52 @@ lb, ub = getbounds(p) # operating on a vector, we get lower and upper bound vect b = getbounds(sys) # Operating on the system, we get a dict ``` -See also: [`ModelingToolkit.dump_variable_metadata`](@ref), [`ModelingToolkit.dump_parameters`](@ref), -[`ModelingToolkit.dump_unknowns`](@ref). +See also: + +```@docs; canonical = false +tunable_parameters +ModelingToolkit.dump_unknowns +ModelingToolkit.dump_parameters +``` + +## Symbolic operators -## Index +ModelingToolkit makes heavy use of "operators". These are custom functions that are applied +to symbolic variables. The most common operator is the `Differential` operator, defined in +Symbolics.jl. -```@index -Pages = ["Variable_metadata.md"] +```@docs +Symbolics.Differential ``` -## Docstrings +ModelingToolkit also defines a plethora of custom operators. -```@autodocs -Modules = [ModelingToolkit] -Pages = ["variables.jl"] -Private = false +```@docs +Pre +Initial +Shift ``` +While not an operator, `ShiftIndex` is commonly used to use `Shift` operators in a more +convenient way when writing discrete systems. + ```@docs -ModelingToolkit.dump_variable_metadata -ModelingToolkit.dump_parameters -ModelingToolkit.dump_unknowns +ShiftIndex +``` + +### Sampled time operators + +The following operators are used in hybrid ODE systems, where part of the dynamics of the +system happen at discrete intervals on a clock. While ModelingToolkit cannot yet simulate +such systems, it has the capability to represent them. + +!!! warn + + These operators are considered experimental API. + +```@docs +Sample +Hold +SampleTime +sampletime ``` From ebde6245b99fc5e26f8321f0e36eb8f44c06e329 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Fri, 23 May 2025 18:18:54 +0530 Subject: [PATCH 15/49] docs: add codegen page to `API` docs --- docs/src/API/codegen.md | 45 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 45 insertions(+) create mode 100644 docs/src/API/codegen.md diff --git a/docs/src/API/codegen.md b/docs/src/API/codegen.md new file mode 100644 index 0000000000..b3a9c01e58 --- /dev/null +++ b/docs/src/API/codegen.md @@ -0,0 +1,45 @@ +# Code generation utilities + +These are lower-level functions that ModelingToolkit leverages to generate code for +building numerical problems. + +```@docs +ModelingToolkit.generate_rhs +ModelingToolkit.generate_diffusion_function +ModelingToolkit.generate_jacobian +ModelingToolkit.generate_tgrad +ModelingToolkit.generate_hessian +ModelingToolkit.generate_W +ModelingToolkit.generate_dae_jacobian +ModelingToolkit.generate_history +ModelingToolkit.generate_boundary_conditions +ModelingToolkit.generate_cost +ModelingToolkit.generate_cost_gradient +ModelingToolkit.generate_cost_hessian +ModelingToolkit.generate_cons +ModelingToolkit.generate_constraint_jacobian +ModelingToolkit.generate_constraint_hessian +ModelingToolkit.generate_control_jacobian +ModelingToolkit.build_explicit_observed_function +``` + +For functions such as jacobian calculation which require symbolic computation, there +are `calculate_*` equivalents to obtain the symbolic result without building a function. + +```@docs +ModelingToolkit.calculate_tgrad +ModelingToolkit.calculate_jacobian +ModelingToolkit.jacobian_sparsity +ModelingToolkit.jacobian_dae_sparsity +ModelingToolkit.calculate_hessian +ModelingToolkit.hessian_sparsity +ModelingToolkit.calculate_massmatrix +ModelingToolkit.W_sparsity +ModelingToolkit.calculate_W_prototype +ModelingToolkit.calculate_cost_gradient +ModelingToolkit.calculate_cost_hessian +ModelingToolkit.cost_hessian_sparsity +ModelingToolkit.calculate_constraint_jacobian +ModelingToolkit.calculate_constraint_hessian +ModelingToolkit.calculate_control_jacobian +``` From 5b64a6b7d2aefd3d1b352151a160c2ca3c36f9a0 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Sat, 24 May 2025 11:17:50 +0530 Subject: [PATCH 16/49] docs: document `nameof` and `brownian` --- demo.jl | 107 +++++++++++++++++++++++ docs/src/getting_started/odes.md | 145 +++++++++++++++++++++++++++++++ src/systems/abstractsystem.jl | 11 +++ 3 files changed, 263 insertions(+) create mode 100644 demo.jl create mode 100644 docs/src/getting_started/odes.md diff --git a/demo.jl b/demo.jl new file mode 100644 index 0000000000..0a1bea54db --- /dev/null +++ b/demo.jl @@ -0,0 +1,107 @@ +macro mtkcompile(ex...) + quote + @mtkbuild $(ex...) + end +end + +function mtkcompile(args...; kwargs...) + structural_simplify(args...; kwargs...) +end + +################################# + +using ModelingToolkit, OrdinaryDiffEq, StochasticDiffEq +using ModelingToolkit: t_nounits as t, D_nounits as D + +## ODEs + +@parameters g +@variables x(t) y(t) λ(t) +eqs = [D(D(x)) ~ λ * x + D(D(y)) ~ λ * y - g + x^2 + y^2 ~ 1] +@mtkbuild pend = System(eqs, t) +prob = ODEProblem(pend, [x => -1, y => 0], (0.0, 10.0), [g => 1], guesses = [λ => 1]) + +sol = solve(prob, FBDF()) + +## SDEs and unified `System` + +@variables x(t) y(t) z(t) +@parameters σ ρ β +@brownian a + +eqs = [ + D(x) ~ σ * (y - x) + 0.1x * a, + D(y) ~ x * (ρ - z) - y + 0.1y * a, + D(z) ~ x * y - β * z + 0.1z * a +] + +@mtkbuild sys1 = System(eqs, t) + +eqs = [ + D(x) ~ σ * (y - x), + D(y) ~ x * (ρ - z) - y, + D(z) ~ x * y - β * z +] + +noiseeqs = [0.1*x; + 0.1*y; + 0.1*z;;] + +@mtkbuild sys2 = SDESystem(eqs, noiseeqs, t) + +u0 = [ + x => 1.0, + y => 0.0, + z => 0.0] + +p = [σ => 28.0, + ρ => 10.0, + β => 8 / 3] + +sdeprob = SDEProblem(sys1, u0, (0.0, 10.0), p) +sdesol = solve(sdeprob, ImplicitEM()) + +odeprob = ODEProblem(sys1, u0, (0.0, 10.0), p) # error! +odeprob = ODEProblem(sys1, u0, (0.0, 10.0), p; check_compatibility = false) + +@variables x y z +@parameters σ ρ β + +# Define a nonlinear system +eqs = [0 ~ σ * (y - x), + y ~ x * (ρ - z), + β * z ~ x * y] +@mtkbuild sys = System(eqs) + +## ImplicitDiscrete Affects + +@parameters g +@variables x(t) y(t) λ(t) +eqs = [D(D(x)) ~ λ * x + D(D(y)) ~ λ * y - g + x^2 + y^2 ~ 1] +c_evt = [t ~ 5.0] => [x ~ Pre(x) + 0.1] +@mtkbuild pend = System(eqs, t, continuous_events = c_evt) +prob = ODEProblem(pend, [x => -1, y => 0], (0.0, 10.0), [g => 1], guesses = [λ => 1]) + +sol = solve(prob, FBDF()) + +## `@named` and `ParentScope` + +function SysA(; name, var1) + @variables x(t) + return System([D(x) ~ var1], t; name) +end +function SysB(; name, var1) + @variables x(t) + @named subsys = SysA(; var1) + return System([D(x) ~ x], t; systems = [subsys], name) +end +function SysC(; name) + @variables x(t) + @named subsys = SysB(; var1 = x) + return System([D(x) ~ x], t; systems = [subsys], name) +end +@mtkbuild sys = SysC() diff --git a/docs/src/getting_started/odes.md b/docs/src/getting_started/odes.md new file mode 100644 index 0000000000..7f51accc16 --- /dev/null +++ b/docs/src/getting_started/odes.md @@ -0,0 +1,145 @@ +# [Building ODEs and DAEs with ModelingToolkit.jl](@id getting_started_ode) + +This is an introductory tutorial for ModelingToolkit.jl (MTK). We will demonstrate the +basics of the package by demontrating how to build systems of Ordinary Differential +Equations (ODEs) and Differential-Algebraic Equations (DAEs). + +## Installing ModelingToolkit + +To install ModelingToolkit, use the Julia package manager. This can be done as follows: + +```julia +using Pkg +Pkg.add("ModelingToolkit") +``` + +## The end goal + +TODO + +## Basics of MTK + +ModelingToolkit.jl is a symbolic-numeric system. This means it allows specifying a model +(such as an ODE) in a similar way to how it would be written on paper. Let's start with a +simple example. The system to be modeled is a first-order lag element: + +```math +\dot{x} = \frac{f(t) - x(t)}{\tau} +``` + +Here, ``t`` is the independent variable (time), ``x(t)`` is the (scalar) unknown +variable, ``f(t)`` is an external forcing function, and ``\tau`` is a +parameter. + +For simplicity, we will start off by setting the forcing function to a constant `1`. Every +ODE has a single independent variable. MTK has a common definition for time `t` and the +derivative with respect to it. + +```@example ode2 +using ModelingToolkit +using ModelingToolkit: t_nounits as t, D_nounits as D +``` + +Next, we declare the (dependent) variables and the parameters of our model: + +```@example ode2 +@variables x(t) +@parameters τ +``` + +Note the syntax `x(t)`. We must declare that the variable `x` is a function of the independent +variable `t`. Next, we define the equations of the system: + +```@example ode2 +eqs = [D(x) ~ (1 - x) / τ] +``` + +Since `=` is reserved as the assignment operator, MTK uses `~` to denote equality between +expressions. Now we must consolidate all of this information about our system of ODEs into +ModelingToolkit's `System` type. + +```@example ode2 +sys = System(eqs, t, [x], [τ]; name = :sys) +``` + +The `System` constructor accepts a `Vector{Equation}` as the first argument, followed by the +independent variable, a list of dependent variables, and a list of parameters. Every system +must be given a name via the `name` keyword argument. Most of the time, we want to name our +system the same as the variable it is assigned to. The `@named` macro helps with this: + +```@example ode2 +@named sys = System(eqs, t, [x], [τ]) +``` + +Additionally, it may become inconvenient to specify all variables and parameters every time +a system is created. MTK allows omitting these arguments, and will automatically infer them +from the equations. + +```@example ode2 +@named sys = System(eqs, t) +``` + +Our system is not quite ready for simulation yet. First, we must use the `mtkcompile` +function which transforms the system into a form that MTK can handle. For our trivial +system, this does not do much. + +```@example ode2 +simp_sys = mtkcompile(sys) +``` + +Since building and simplifying a system is a common workflow, MTK provides the `@mtkcompile` +macro for convenience. + +```@example ode2 +@mtkcompile sys = System(eqs, t) +``` + +We can now build an `ODEProblem` from the system. ModelingToolkit generates the necessary +code for numerical ODE solvers to solve this system. We need to provide an initial +value for the variable `x` and a value for the parameter `p`, as well as the time span +for which to simulate the system. + +```@example ode2 +prob = ODEProblem(sys, [x => 0.0, τ => 3.0], (0.0, 10.0)) +``` + +Here, we are saying that `x` should start at `0.0`, `τ` should be `3.0` and the system +should be simulated from `t = 0.0` to `t = 10.0`. To solve the system, we must import a +solver. + +```@example ode2 +using OrdinaryDiffEq + +sol = solve(prob) +``` + +[OrdinaryDiffEq.jl](https://docs.sciml.ai/DiffEqDocs/stable/) contains a large number of +numerical solvers. It also comes with a default solver which is used when calling +`solve(prob)` and is capable of handling a large variety of systems. + +We can obtain the timeseries of `x` by indexing the solution with the symbolic variable: + +```@example ode2 +sol[x] +``` + +We can even obtain timeseries of complicated expressions involving the symbolic variables in +the model + +```@example ode2 +sol[(1 - x) / τ] +``` + +Perhaps more interesting is a plot of the solution. This can easily be achieved using Plots.jl. + +```@example ode2 +using Plots + +plot(sol) +``` + +Similarly, we can plot different expressions: + +```@example ode2 +plot(sol; idxs = (1 - x) / τ) +``` diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl index b997219eba..549a86afda 100644 --- a/src/systems/abstractsystem.jl +++ b/src/systems/abstractsystem.jl @@ -197,7 +197,18 @@ end const MTKPARAMETERS_ARG = Sym{Vector{Vector}}(:___mtkparameters___) +""" + $(TYPEDSIGNATURES) + +Obtain the name of `sys`. +""" Base.nameof(sys::AbstractSystem) = getfield(sys, :name) +""" + $(TYPEDSIGNATURES) + +Obtain the description associated with `sys` if present, and an empty +string otherwise. +""" description(sys::AbstractSystem) = has_description(sys) ? get_description(sys) : "" """ From 292ff836b6d7accb1b54e4105d5f2d8fb58b059e Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Sat, 24 May 2025 11:19:07 +0530 Subject: [PATCH 17/49] docs: add docs for `jumps`, `brownians`, `cost`, `constraints`, `symbolic_tstops`, `preface` --- src/systems/abstractsystem.jl | 36 ++++++++++++++++++++++++++++++++++- 1 file changed, 35 insertions(+), 1 deletion(-) diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl index 549a86afda..8f43d2f66f 100644 --- a/src/systems/abstractsystem.jl +++ b/src/systems/abstractsystem.jl @@ -1679,6 +1679,12 @@ function equations_toplevel(sys::AbstractSystem) return get_eqs(sys) end +""" + $(TYPEDSIGNATURES) + +Get the flattened jumps of the system. In other words, obtain all of the jumps in `sys` and +all the subsystems of `sys` (appropriately namespaced). +""" function jumps(sys::AbstractSystem) js = get_jumps(sys) systems = get_systems(sys) @@ -1688,6 +1694,12 @@ function jumps(sys::AbstractSystem) return [js; reduce(vcat, namespace_jumps.(systems); init = [])] end +""" + $(TYPEDSIGNATURES) + +Get all of the brownian variables involved in the system `sys` and all subsystems, +appropriately namespaced. +""" function brownians(sys::AbstractSystem) bs = get_brownians(sys) systems = get_systems(sys) @@ -1697,6 +1709,12 @@ function brownians(sys::AbstractSystem) return [bs; reduce(vcat, namespace_brownians.(systems); init = [])] end +""" + $(TYPEDSIGNATURES) + +Recursively consolidate the cost vector of `sys` and all subsystems of `sys`, returning the +final scalar cost function. +""" function cost(sys::AbstractSystem) cs = get_costs(sys) consolidate = get_consolidate(sys) @@ -1726,7 +1744,12 @@ function namespace_constraints(sys) map(cstr -> namespace_constraint(cstr, sys), cstrs) end -function constraints(sys) +""" + $(TYPEDSIGNATURES) + +Get all constraints in the system `sys` and all of its subsystems, appropriately namespaced. +""" +function constraints(sys::AbstractSystem) cs = get_constraints(sys) systems = get_systems(sys) isempty(systems) ? cs : [cs; reduce(vcat, namespace_constraints.(systems))] @@ -1753,6 +1776,11 @@ function initialization_equations(sys::AbstractSystem) end end +""" + $(TYPEDSIGNATURES) + +Get the tstops present in `sys` and its subsystems, appropriately namespaced. +""" function symbolic_tstops(sys::AbstractSystem) tstops = get_tstops(sys) systems = get_systems(sys) @@ -1761,6 +1789,12 @@ function symbolic_tstops(sys::AbstractSystem) return tstops end +""" + $(TYPEDSIGNATURES) + +Obtain the preface associated with `sys` and all of its subsystems, appropriately +namespaced. +""" function preface(sys::AbstractSystem) has_preface(sys) || return nothing pre = get_preface(sys) From 2f5f27fdf6ca153870952ddffcfd091ef611f9c9 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Sat, 24 May 2025 11:19:47 +0530 Subject: [PATCH 18/49] refactor: move `full_equations` to `ModelingToolkit`, don't error on singular systems --- .../symbolics_tearing.jl | 42 --------------- src/systems/abstractsystem.jl | 51 +++++++++++++++++++ 2 files changed, 51 insertions(+), 42 deletions(-) diff --git a/src/structural_transformation/symbolics_tearing.jl b/src/structural_transformation/symbolics_tearing.jl index a608b8b5ee..9eccc309bd 100644 --- a/src/structural_transformation/symbolics_tearing.jl +++ b/src/structural_transformation/symbolics_tearing.jl @@ -99,48 +99,6 @@ function eq_derivative!(ts::TearingState, ieq::Int; kwargs...) return eq_diff end -function tearing_sub(expr, dict, s) - expr = Symbolics.fixpoint_sub(expr, dict; operator = ModelingToolkit.Initial) - s ? simplify(expr) : expr -end - -function tearing_substitute_expr(sys::AbstractSystem, expr; simplify = false) - empty_substitutions(sys) && return expr - substitutions = get_substitutions(sys) - return tearing_sub(expr, substitutions, simplify) -end - -""" -$(TYPEDSIGNATURES) - -Like `equations(sys)`, but includes substitutions done by the tearing process. -These equations matches generated numerical code. - -See also [`equations`](@ref) and [`ModelingToolkit.get_eqs`](@ref). -""" -function full_equations(sys::AbstractSystem; simplify = false) - isempty(observed(sys)) && return equations(sys) - subs = Dict([eq.lhs => eq.rhs for eq in observed(sys)]) - neweqs = map(equations(sys)) do eq - if iscall(eq.lhs) && operation(eq.lhs) isa Union{Shift, Differential} - return tearing_sub(eq.lhs, subs, simplify) ~ tearing_sub(eq.rhs, subs, - simplify) - else - if !(eq.lhs isa Number && eq.lhs == 0) - eq = 0 ~ eq.rhs - eq.lhs - end - rhs = tearing_sub(eq.rhs, subs, simplify) - if rhs isa Symbolic - return 0 ~ rhs - else # a number - error("tearing failed because the system is singular") - end - end - eq - end - return neweqs -end - function tearing_substitution(sys::AbstractSystem; kwargs...) neweqs = full_equations(sys::AbstractSystem; kwargs...) @set! sys.eqs = neweqs diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl index 8f43d2f66f..d736edf614 100644 --- a/src/systems/abstractsystem.jl +++ b/src/systems/abstractsystem.jl @@ -1679,6 +1679,57 @@ function equations_toplevel(sys::AbstractSystem) return get_eqs(sys) end +""" + $(TYPEDSIGNATURES) + +Recursively substitute `dict` into `expr`. Use `Symbolics.simplify` on the expression +if `simplify == true`. +""" +function substitute_and_simplify(expr, dict::AbstractDict, simplify::Bool) + expr = Symbolics.fixpoint_sub(expr, dict; operator = ModelingToolkit.Initial) + simplify ? Symbolics.simplify(expr) : expr +end + +""" + $(TYPEDSIGNATURES) + +Recursively substitute the observed equations of `sys` into `expr`. If `simplify`, call +`Symbolics.simplify` on the resultant expression. +""" +function substitute_observed(sys::AbstractSystem, expr; simplify = false) + empty_substitutions(sys) && return expr + substitutions = get_substitutions(sys) + return substitute_and_simplify(expr, substitutions, simplify) +end + +""" +$(TYPEDSIGNATURES) + +Like `equations(sys)`, but also substitutes the observed equations eliminated from the +equations during `mtkcompile`. These equations matches generated numerical code. + +See also [`equations`](@ref) and [`ModelingToolkit.get_eqs`](@ref). +""" +function full_equations(sys::AbstractSystem; simplify = false) + empty_substitutions(sys) && return equations(sys) + subs = get_substitutions(sys) + neweqs = map(equations(sys)) do eq + if iscall(eq.lhs) && operation(eq.lhs) isa Union{Shift, Differential} + return substitute_and_simplify(eq.lhs, subs, simplify) ~ substitute_and_simplify( + eq.rhs, subs, + simplify) + else + if !_iszero(eq.lhs) + eq = 0 ~ eq.rhs - eq.lhs + end + rhs = substitute_and_simplify(eq.rhs, subs, simplify) + return 0 ~ rhs + end + eq + end + return neweqs +end + """ $(TYPEDSIGNATURES) From 04c02bc43903a0737a537850ab6a93f3109299da Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Sat, 24 May 2025 11:20:06 +0530 Subject: [PATCH 19/49] refactor: use new `substitute_observed` instead of `tearing_substitute_expr` --- src/systems/systems.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/systems/systems.jl b/src/systems/systems.jl index 13a62f7915..ceab3c6c2f 100644 --- a/src/systems/systems.jl +++ b/src/systems/systems.jl @@ -147,7 +147,7 @@ function __mtkcompile(sys::AbstractSystem; simplify = false, is_scalar_noise = false end - noise_eqs = StructuralTransformations.tearing_substitute_expr(ode_sys, noise_eqs) + noise_eqs = substitute_observed(ode_sys, noise_eqs) ssys = System(Vector{Equation}(full_equations(ode_sys)), get_iv(ode_sys), unknowns(ode_sys), parameters(ode_sys); noise_eqs, name = nameof(ode_sys), observed = observed(ode_sys), defaults = defaults(sys), From d507ff969d568aadcb0d925bec92d0bdfe327f52 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Sat, 24 May 2025 11:20:25 +0530 Subject: [PATCH 20/49] docs: add docs for `empty_substitutions` and `get_substitutions` --- src/utils.jl | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/src/utils.jl b/src/utils.jl index e2bc68e03f..5397b23077 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -693,10 +693,21 @@ end isarray(x) = x isa AbstractArray || x isa Symbolics.Arr +""" + $(TYPEDSIGNATURES) + +Check if any variables were eliminated from the system as part of `mtkcompile`. +""" function empty_substitutions(sys) isempty(observed(sys)) end +""" + $(TYPEDSIGNATURES) + +Get a dictionary mapping variables eliminated from the system during `mtkcompile` to the +expressions used to calculate them. +""" function get_substitutions(sys) Dict([eq.lhs => eq.rhs for eq in observed(sys)]) end From c5cc07af1217ee56afc175fb33aaa72751814f38 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Sat, 24 May 2025 12:40:50 +0530 Subject: [PATCH 21/49] docs: add docs for `apply_to_variables` --- src/systems/abstractsystem.jl | 14 ++++++++++++-- 1 file changed, 12 insertions(+), 2 deletions(-) diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl index d736edf614..e2bd39ca5c 100644 --- a/src/systems/abstractsystem.jl +++ b/src/systems/abstractsystem.jl @@ -1048,8 +1048,18 @@ function Base.setproperty!(sys::AbstractSystem, prop::Symbol, val) """) end -apply_to_variables(f::F, ex) where {F} = _apply_to_variables(f, ex) -apply_to_variables(f::F, ex::Num) where {F} = wrap(_apply_to_variables(f, unwrap(ex))) +""" + $(TYPEDSIGNATURES) + +Apply function `f` to each variable in expression `ex`. `f` should be a function that takes +a variable and returns the replacement to use. A "variable" in this context refers to a +symbolic quantity created directly from a variable creation macro such as +[`Symbolics.@variables`](@ref), [`@independent_variables`](@ref), [`@parameters`](@ref), +[`@constants`](@ref) or [`@brownians`](@ref). +""" +apply_to_variables(f, ex) = _apply_to_variables(f, ex) +apply_to_variables(f, ex::Num) = wrap(_apply_to_variables(f, unwrap(ex))) +apply_to_variables(f, ex::Symbolics.Arr) = wrap(_apply_to_variables(f, unwrap(ex))) function _apply_to_variables(f::F, ex) where {F} if isvariable(ex) return f(ex) From e1ca3682cc243b1aa578cf1db6c3ef2bcfe0c82c Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Sat, 24 May 2025 12:41:12 +0530 Subject: [PATCH 22/49] docs: add docs for `@component` and `@connector` --- src/systems/abstractsystem.jl | 20 ++++++++++++++++++++ src/systems/connectors.jl | 31 +++++++++++++++++++++++++++++++ 2 files changed, 51 insertions(+) diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl index e2bd39ca5c..1d97417d93 100644 --- a/src/systems/abstractsystem.jl +++ b/src/systems/abstractsystem.jl @@ -2471,6 +2471,26 @@ function component_post_processing(expr, isconnector) end end +""" + $(TYPEDSIGNATURES) + +Mark a system constructor function as building a component. For example, + +```julia +@component function AddOne(; name) + @variables in(t) out(t) + eqs = [out ~ in + 1] + return System(eqs, t, [in, out], []; name) +end +``` + +ModelingToolkit systems are either components or connectors. Components define dynamics of +the model. Connectors are used to connect components together. See the +[Model building reference](@ref model_building_api) section of the documentation for more +information. + +See also: [`@connector`](@ref). +""" macro component(expr) esc(component_post_processing(expr, false)) end diff --git a/src/systems/connectors.jl b/src/systems/connectors.jl index ba41b2b011..b4a3d47711 100644 --- a/src/systems/connectors.jl +++ b/src/systems/connectors.jl @@ -20,6 +20,37 @@ function get_connection_type(s) getmetadata(s, VariableConnectType, Equality) end +""" + $(TYPEDSIGNATURES) + +Mark a system constructor function as building a connector. For example, + +```julia +@connector function ElectricalPin(; name, v = nothing, i = nothing) + @variables begin + v(t) = v, [description = "Potential at the pin [V]"] + i(t) = i, [connect = Flow, description = "Current flowing into the pin [A]"] + end + return System(Equation[], t, [v, i], []; name) +end +``` + +Since connectors only declare variables, the equivalent shorthand syntax can also be used: + +```julia +@connector Pin begin + v(t), [description = "Potential at the pin [V]"] + i(t), [connect = Flow, description = "Current flowing into the pin [A]"] +end +``` + +ModelingToolkit systems are either components or connectors. Components define dynamics of +the model. Connectors are used to connect components together. See the +[Model building reference](@ref model_building_api) section of the documentation for more +information. + +See also: [`@component`](@ref). +""" macro connector(expr) esc(component_post_processing(expr, true)) end From 81a2cd64a0e051cac224467a968821d4261d0378 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Sat, 24 May 2025 12:41:42 +0530 Subject: [PATCH 23/49] docs: add docs for `hierarchy`, `extend`, `compose` --- src/systems/abstractsystem.jl | 49 ++++++++++++++++++++++++++++++++--- 1 file changed, 46 insertions(+), 3 deletions(-) diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl index 1d97417d93..3ad7c35923 100644 --- a/src/systems/abstractsystem.jl +++ b/src/systems/abstractsystem.jl @@ -2619,6 +2619,11 @@ end hierarchy(sys::AbstractSystem; describe = false, bold = describe, kwargs...) Print a tree of a system's hierarchy of subsystems. + +# Keyword arguments + +- `describe`: Whether to also print the description of each subsystem, if present. +- `bold`: Whether to print the name of the system in **bold** font. """ function hierarchy(sys::AbstractSystem; describe = false, bold = describe, kwargs...) print_tree(sys; printnode_kw = (describe = describe, bold = bold), kwargs...) @@ -2661,9 +2666,14 @@ end """ $(TYPEDSIGNATURES) -Extend `basesys` with `sys`. +Extend `basesys` with `sys`. This can be thought of as the `merge` operation on systems. +Values in `sys` take priority over duplicates in `basesys` (for example, defaults). + By default, the resulting system inherits `sys`'s name and description. +The `&` operator can also be used for this purpose. `sys & basesys` is equivalent to +`extend(sys, basesys)`. + See also [`compose`](@ref). """ function extend(sys::AbstractSystem, basesys::AbstractSystem; @@ -2710,14 +2720,31 @@ function extend(sys::AbstractSystem, basesys::AbstractSystem; return T(args...; kwargs...) end +""" + $(TYPEDSIGNATURES) + +Extend `sys` with all systems in `basesys` in order. +""" function extend(sys, basesys::Vector{T}) where {T <: AbstractSystem} foldl(extend, basesys, init = sys) end +""" + $(TYPEDSIGNATURES) + +Syntactic sugar for `extend(sys, basesys)`. + +See also: [`extend`](@ref). +""" function Base.:(&)(sys::AbstractSystem, basesys::AbstractSystem; kwargs...) extend(sys, basesys; kwargs...) end +""" + $(TYPEDSIGNATURES) + +Syntactic sugar for `extend(sys, basesys)`. +""" function Base.:(&)( sys::AbstractSystem, basesys::Vector{T}; kwargs...) where {T <: AbstractSystem} extend(sys, basesys; kwargs...) @@ -2726,8 +2753,11 @@ end """ $(SIGNATURES) -Compose multiple systems together. The resulting system would inherit the first -system's name. +Compose multiple systems together. This adds all of `systems` as subsystems of `sys`. +The resulting system inherits the name of `sys` by default. + +The `∘` operator can also be used for this purpose. `sys ∘ basesys` is equivalent to +`compose(sys, basesys)`. See also [`extend`](@ref). """ @@ -2749,9 +2779,22 @@ function compose(sys::AbstractSystem, systems::AbstractArray; name = nameof(sys) @set! sys.ps = unique!(vcat(get_ps(sys), collect(newparams))) return sys end +""" + $(TYPEDSIGNATURES) + +Syntactic sugar for adding all systems in `syss` as the subsystems of `first(syss)`. +""" function compose(syss...; name = nameof(first(syss))) compose(first(syss), collect(syss[2:end]); name = name) end + +""" + $(TYPEDSIGNATURES) + +Syntactic sugar for `compose(sys1, sys2)`. + +See also: [`compose`](@ref). +""" Base.:(∘)(sys1::AbstractSystem, sys2::AbstractSystem) = compose(sys1, sys2) function UnPack.unpack(sys::ModelingToolkit.AbstractSystem, ::Val{p}) where {p} From 51bd5e8038b9ff4c60d97fe56d97d459477edd96 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Sat, 24 May 2025 12:42:01 +0530 Subject: [PATCH 24/49] docs: add docs for `Equality`, `Flow`, `Stream` and `instream` --- src/systems/connectors.jl | 9 +++++++++ src/variables.jl | 42 +++++++++++++++++++++++++++++++++++++++ 2 files changed, 51 insertions(+) diff --git a/src/systems/connectors.jl b/src/systems/connectors.jl index b4a3d47711..473b1bd680 100644 --- a/src/systems/connectors.jl +++ b/src/systems/connectors.jl @@ -94,6 +94,15 @@ is_domain_connector(s) = isconnector(s) && get_connector_type(s) === DomainConne get_systems(c::Connection) = c.systems +""" + $(TYPEDSIGNATURES) + +`instream` is used when modeling stream connections. It is only allowed to be used on +`Stream` variables. + +Refer to the [Connection semantics](@ref connect_semantics) section of the docs for more +information. +""" instream(a) = term(instream, unwrap(a), type = symtype(a)) SymbolicUtils.promote_symtype(::typeof(instream), _) = Real diff --git a/src/variables.jl b/src/variables.jl index 230c98a985..746548b432 100644 --- a/src/variables.jl +++ b/src/variables.jl @@ -91,8 +91,50 @@ end ### Connect abstract type AbstractConnectType end +""" + $(TYPEDEF) + +Flag which is meant to be passed to the `connect` metadata of a variable to affect how it +behaves when the connector it is in is part of a `connect` equation. `Equality` is the +default value and such variables when connected are made equal. For example, electric +potential is equated at a junction. + +For more information, refer to the [Connection semantics](@ref connect_semantics) section +of the docs. + +See also: [`Symbolics.connect`](@ref), [`@connector`](@ref), [`Flow`](@ref), +[`Stream`](@ref). +""" struct Equality <: AbstractConnectType end # Equality connection +""" + $(TYPEDEF) + +Flag which is meant to be passed to the `connect` metadata of a variable to affect how it +behaves when the connector it is in is part of a `connect` equation. `Flow` denotes that +the sum of marked variable in all connectors in the connection set must sum to zero. For +example, electric current sums to zero at a junction (assuming appropriate signs are used +for current flowing in and out of the function). + +For more information, refer to the [Connection semantics](@ref connect_semantics) section +of the docs. + +See also: [`Symbolics.connect`](@ref), [`@connector`](@ref), [`Equality`](@ref), +[`Stream`](@ref). +""" struct Flow <: AbstractConnectType end # sum to 0 +""" + $(TYPEDEF) + +Flag which is meant to be passed to the `connect` metadata of a variable to affect how it +behaves when the connector it is in is part of a `connect` equation. `Stream` denotes that +the variable is part of a special stream connector. + +For more information, refer to the [Connection semantics](@ref connect_semantics) section +of the docs. + +See also: [`Symbolics.connect`](@ref), [`@connector`](@ref), [`Equality`](@ref), +[`Flow`](@ref). +""" struct Stream <: AbstractConnectType end # special stream connector """ From 950cede02b8d2aba10083f178748218ac20b9f1f Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Sat, 24 May 2025 12:42:28 +0530 Subject: [PATCH 25/49] chore: add SciMLPublic.jl as a dependency --- Project.toml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/Project.toml b/Project.toml index 41568022ab..527451bca2 100644 --- a/Project.toml +++ b/Project.toml @@ -51,6 +51,7 @@ Reexport = "189a3867-3050-52da-a836-e630ba90ab69" RuntimeGeneratedFunctions = "7e49a35a-f44a-4d26-94aa-eba1b4ca6b47" SCCNonlinearSolve = "9dfe8606-65a1-4bb3-9748-cb89d1561431" SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" +SciMLPublic = "431bcebd-1456-4ced-9d72-93c2757fff0b" SciMLStructures = "53ae85a6-f571-4167-b2af-e1d143709226" Serialization = "9e88b42a-f829-5b0c-bbe9-9e923198166b" Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46" @@ -145,6 +146,7 @@ Reexport = "0.2, 1" RuntimeGeneratedFunctions = "0.5.9" SCCNonlinearSolve = "1.0.0" SciMLBase = "2.91.1" +SciMLPublic = "1.0.0" SciMLStructures = "1.7" Serialization = "1" Setfield = "0.7, 0.8, 1" From ee8a461a0aabf58dcb05f628b964daa497c95adf Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Sat, 24 May 2025 12:42:39 +0530 Subject: [PATCH 26/49] chore: mark `apply_to_variables` as public --- src/ModelingToolkit.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/ModelingToolkit.jl b/src/ModelingToolkit.jl index 9bb96e1e14..79d2afcd67 100644 --- a/src/ModelingToolkit.jl +++ b/src/ModelingToolkit.jl @@ -99,6 +99,7 @@ const DQ = DynamicQuantities import DifferentiationInterface as DI using ADTypes: AutoForwardDiff +import SciMLPublic: @public export @derivatives @@ -358,4 +359,6 @@ export AbstractDynamicOptProblem, JuMPDynamicOptProblem, InfiniteOptDynamicOptPr CasADiDynamicOptProblem export DynamicOptSolution +@public apply_to_variables + end # module From 6367ed389ca8b408d92089983977d8750d8bee3b Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Sat, 24 May 2025 15:29:56 +0530 Subject: [PATCH 27/49] docs: improve docstring for `complete` --- src/systems/abstractsystem.jl | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl index 3ad7c35923..40029086a1 100644 --- a/src/systems/abstractsystem.jl +++ b/src/systems/abstractsystem.jl @@ -732,6 +732,9 @@ the global structure of the system. One property to note is that if a system is complete, the system will no longer namespace its subsystems or variables, i.e. `isequal(complete(sys).v.i, v.i)`. + +This namespacing functionality can also be toggled independently of `complete` +using [`toggle_namespacing`](@ref). """ function complete( sys::AbstractSystem; split = true, flatten = true, add_initial_parameters = true) @@ -2495,6 +2498,22 @@ macro component(expr) esc(component_post_processing(expr, false)) end +""" + $(TYPEDSIGNATURES) + +Macro shorthand for building and compiling a system in one step. + +```julia +@mtkcompile sys = Constructor(args...; kwargs....) +``` + +Is shorthand for + +```julia +@named sys = Constructor(args...; kwargs...) +sys = mtkcompile(sys) +``` +""" macro mtkcompile(exprs...) expr = exprs[1] named_expr = ModelingToolkit.named_expr(expr) From f4783b027983e9cf8a31781f8a9ac8dc7f5a2919 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Sat, 24 May 2025 15:30:17 +0530 Subject: [PATCH 28/49] docs: document `IfLifting` as experimental --- src/systems/if_lifting.jl | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/systems/if_lifting.jl b/src/systems/if_lifting.jl index 130a4bdab4..c2eba035ff 100644 --- a/src/systems/if_lifting.jl +++ b/src/systems/if_lifting.jl @@ -410,6 +410,11 @@ If lifting converts (nested) if statements into a series of continuous events + Lifting proceeds through the following process: * rewrite comparisons to be of the form eqn [op] 0; subtract the RHS from the LHS * replace comparisons with generated parameters; for each comparison eqn [op] 0, generate an event (dependent on op) that sets the parameter + +!!! warn + + This is an experimental simplification pass. It may have bugs. Please open issues with + MWEs for any bugs encountered while using this. """ function IfLifting(sys::System) if !is_time_dependent(sys) From 7769b5ba0d6167a396f34acc4eff8f03672bc4fc Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Sat, 24 May 2025 15:30:33 +0530 Subject: [PATCH 29/49] docs: improve docstring for `mtkcompile` --- src/systems/systems.jl | 14 +++++++++++--- 1 file changed, 11 insertions(+), 3 deletions(-) diff --git a/src/systems/systems.jl b/src/systems/systems.jl index ceab3c6c2f..63585541e0 100644 --- a/src/systems/systems.jl +++ b/src/systems/systems.jl @@ -9,10 +9,18 @@ end """ $(SIGNATURES) -Structurally simplify algebraic equations in a system and compute the -topological sort of the observed equations in `sys`. +Compile the given system into a form that ModelingToolkit can generate code for. Also +performs a variety of symbolic-numeric enhancements. For ODEs, this includes processes +such as order reduction, index reduction, alias elimination and tearing. A subset of the +unknowns of the system may be eliminated as observables, eliminating the need for the +numerical solver to solve for these variables. + +Does not rely on metadata to identify variables/parameters/brownians. Instead, queries +the system for which symbolic quantites belong to which category. Any variables not +present in the equations of the system will be removed in this process. + +# Keyword Arguments -### Optional Keyword Arguments: + When `simplify=true`, the `simplify` function will be applied during the tearing process. + `allow_symbolic=false`, `allow_parameter=true`, and `conservative=false` limit the coefficient types during tearing. In particular, `conservative=true` limits tearing to only solve for trivial linear systems where the coefficient has the absolute value of ``1``. + `fully_determined=true` controls whether or not an error will be thrown if the number of equations don't match the number of inputs, outputs, and equations. From ac0ad247a2519fc9c25be5367a5964ac5bc7dc90 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Sat, 24 May 2025 15:30:53 +0530 Subject: [PATCH 30/49] refactor: export `Girsanov_transform` --- src/ModelingToolkit.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ModelingToolkit.jl b/src/ModelingToolkit.jl index 79d2afcd67..08a1b8e03e 100644 --- a/src/ModelingToolkit.jl +++ b/src/ModelingToolkit.jl @@ -296,7 +296,7 @@ export isinput, isoutput, getbounds, hasbounds, getguess, hasguess, isdisturbanc hasunit, getunit, hasconnect, getconnect, hasmisc, getmisc, state_priority export liouville_transform, change_independent_variable, substitute_component, - add_accumulations, noise_to_brownians + add_accumulations, noise_to_brownians, Girsanov_transform export PDESystem export Differential, expand_derivatives, @derivatives export Equation, ConstrainedEquation From 62a2bfaff784570cefecfbfeed5f6ce6f41edaf5 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Sat, 24 May 2025 15:31:36 +0530 Subject: [PATCH 31/49] refactor: move `Connection` and `connect` here instead of Symbolics --- src/ModelingToolkit.jl | 2 +- src/systems/analysis_points.jl | 6 ++--- src/systems/connectors.jl | 46 +++++++++++++++++++++++++++++++++- src/systems/systemstructure.jl | 4 +-- src/variables.jl | 6 ++--- 5 files changed, 54 insertions(+), 10 deletions(-) diff --git a/src/ModelingToolkit.jl b/src/ModelingToolkit.jl index 08a1b8e03e..c99ac17c0c 100644 --- a/src/ModelingToolkit.jl +++ b/src/ModelingToolkit.jl @@ -72,7 +72,7 @@ using RuntimeGeneratedFunctions: drop_expr using Symbolics: degree using Symbolics: _parse_vars, value, @derivatives, get_variables, exprs_occur_in, symbolic_linear_solve, build_expr, unwrap, wrap, - VariableSource, getname, variable, Connection, connect, + VariableSource, getname, variable, NAMESPACE_SEPARATOR, set_scalar_metadata, setdefaultval, initial_state, transition, activeState, entry, hasnode, ticksInState, timeInState, fixpoint_sub, fast_substitute, diff --git a/src/systems/analysis_points.jl b/src/systems/analysis_points.jl index 2ce6356aa7..12185c943f 100644 --- a/src/systems/analysis_points.jl +++ b/src/systems/analysis_points.jl @@ -176,7 +176,7 @@ function renamespace(sys, ap::AnalysisPoint) end # create analysis points via `connect` -function Symbolics.connect(in, ap::AnalysisPoint, outs...; verbose = true) +function connect(in, ap::AnalysisPoint, outs...; verbose = true) return AnalysisPoint() ~ AnalysisPoint(in, ap.name, collect(outs); verbose) end @@ -217,11 +217,11 @@ typically is not (unless the model is an inverse model). - `verbose`: Warn if an input is connected to an output (reverse causality). Silence this warning if you are analyzing an inverse model. """ -function Symbolics.connect(in::AbstractSystem, name::Symbol, out, outs...; verbose = true) +function connect(in::AbstractSystem, name::Symbol, out, outs...; verbose = true) return AnalysisPoint() ~ AnalysisPoint(in, name, [out; collect(outs)]; verbose) end -function Symbolics.connect( +function connect( in::ConnectableSymbolicT, name::Symbol, out::ConnectableSymbolicT, outs::ConnectableSymbolicT...; verbose = true) allvars = (in, out, outs...) diff --git a/src/systems/connectors.jl b/src/systems/connectors.jl index 473b1bd680..af5be85096 100644 --- a/src/systems/connectors.jl +++ b/src/systems/connectors.jl @@ -1,4 +1,48 @@ using Symbolics: StateMachineOperator + +""" + $(TYPEDEF) + +Struct used to represent a connection equation. A connection equation is an `Equation` +where the LHS is an empty `Connection(nothing)` and the RHS is a `Connection` containing +the connected connectors. + +For special types of connections, the LHS `Connection` can contain relevant metadata. +""" +struct Connection + systems::Any +end + +Base.broadcastable(x::Connection) = Ref(x) +Connection() = Connection(nothing) +Base.hash(c::Connection, seed::UInt) = hash(c.systems, (0xc80093537bdc1311 % UInt) ⊻ seed) +Symbolics.hide_lhs(_::Connection) = true + +""" + $(TYPEDSIGNATURES) + +Connect multiple connectors created via `@connector`. All connected connectors +must be unique. +""" +function connect(sys1::AbstractSystem, sys2::AbstractSystem, syss::AbstractSystem...) + syss = (sys1, sys2, syss...) + length(unique(nameof, syss)) == length(syss) || error("connect takes distinct systems!") + Equation(Connection(), Connection(syss)) # the RHS are connected systems +end + +function Base.show(io::IO, c::Connection) + print(io, "connect(") + if c.systems isa AbstractArray || c.systems isa Tuple + n = length(c.systems) + for (i, s) in enumerate(c.systems) + str = join(split(string(nameof(s)), NAMESPACE_SEPARATOR), '.') + print(io, str) + i != n && print(io, ", ") + end + end + print(io, ")") +end + isconnection(_) = false isconnection(_::Connection) = true """ @@ -188,7 +232,7 @@ var1 ~ var3 # ... ``` """ -function Symbolics.connect(var1::ConnectableSymbolicT, var2::ConnectableSymbolicT, +function connect(var1::ConnectableSymbolicT, var2::ConnectableSymbolicT, vars::ConnectableSymbolicT...) allvars = (var1, var2, vars...) validate_causal_variables_connection(allvars) diff --git a/src/systems/systemstructure.jl b/src/systems/systemstructure.jl index 16551008e9..4b01917bcf 100644 --- a/src/systems/systemstructure.jl +++ b/src/systems/systemstructure.jl @@ -1,11 +1,11 @@ using DataStructures -using Symbolics: linear_expansion, unwrap, Connection +using Symbolics: linear_expansion, unwrap using SymbolicUtils: iscall, operation, arguments, Symbolic using SymbolicUtils: quick_cancel, maketerm using ..ModelingToolkit import ..ModelingToolkit: isdiffeq, var_from_nested_derivative, vars!, flatten, value, InvalidSystemException, isdifferential, _iszero, - isparameter, + isparameter, Connection, independent_variables, SparseMatrixCLIL, AbstractSystem, equations, isirreducible, input_timedomain, TimeDomain, InferredTimeDomain, diff --git a/src/variables.jl b/src/variables.jl index 746548b432..104616f6e4 100644 --- a/src/variables.jl +++ b/src/variables.jl @@ -102,7 +102,7 @@ potential is equated at a junction. For more information, refer to the [Connection semantics](@ref connect_semantics) section of the docs. -See also: [`Symbolics.connect`](@ref), [`@connector`](@ref), [`Flow`](@ref), +See also: [`connect`](@ref), [`@connector`](@ref), [`Flow`](@ref), [`Stream`](@ref). """ struct Equality <: AbstractConnectType end # Equality connection @@ -118,7 +118,7 @@ for current flowing in and out of the function). For more information, refer to the [Connection semantics](@ref connect_semantics) section of the docs. -See also: [`Symbolics.connect`](@ref), [`@connector`](@ref), [`Equality`](@ref), +See also: [`connect`](@ref), [`@connector`](@ref), [`Equality`](@ref), [`Stream`](@ref). """ struct Flow <: AbstractConnectType end # sum to 0 @@ -132,7 +132,7 @@ the variable is part of a special stream connector. For more information, refer to the [Connection semantics](@ref connect_semantics) section of the docs. -See also: [`Symbolics.connect`](@ref), [`@connector`](@ref), [`Equality`](@ref), +See also: [`connect`](@ref), [`@connector`](@ref), [`Equality`](@ref), [`Flow`](@ref). """ struct Stream <: AbstractConnectType end # special stream connector From afc33dfb57cf24a04d1949af14cf0f329d056cba Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Sat, 24 May 2025 19:50:51 +0530 Subject: [PATCH 32/49] feat: move statemachine code to MTK --- src/ModelingToolkit.jl | 4 +- src/systems/connectors.jl | 2 - src/systems/state_machines.jl | 155 ++++++++++++++++++++++++++++++++++ 3 files changed, 157 insertions(+), 4 deletions(-) create mode 100644 src/systems/state_machines.jl diff --git a/src/ModelingToolkit.jl b/src/ModelingToolkit.jl index c99ac17c0c..033d83d20b 100644 --- a/src/ModelingToolkit.jl +++ b/src/ModelingToolkit.jl @@ -74,8 +74,7 @@ using Symbolics: _parse_vars, value, @derivatives, get_variables, exprs_occur_in, symbolic_linear_solve, build_expr, unwrap, wrap, VariableSource, getname, variable, NAMESPACE_SEPARATOR, set_scalar_metadata, setdefaultval, - initial_state, transition, activeState, entry, hasnode, - ticksInState, timeInState, fixpoint_sub, fast_substitute, + hasnode, fixpoint_sub, fast_substitute, CallWithMetadata, CallWithParent const NAMESPACE_SEPARATOR_SYMBOL = Symbol(NAMESPACE_SEPARATOR) import Symbolics: rename, get_variables!, _solve, hessian_sparsity, @@ -163,6 +162,7 @@ include("systems/parameter_buffer.jl") include("systems/abstractsystem.jl") include("systems/model_parsing.jl") include("systems/connectors.jl") +include("systems/state_machines.jl") include("systems/analysis_points.jl") include("systems/imperative_affect.jl") include("systems/callbacks.jl") diff --git a/src/systems/connectors.jl b/src/systems/connectors.jl index af5be85096..9761d18493 100644 --- a/src/systems/connectors.jl +++ b/src/systems/connectors.jl @@ -1,5 +1,3 @@ -using Symbolics: StateMachineOperator - """ $(TYPEDEF) diff --git a/src/systems/state_machines.jl b/src/systems/state_machines.jl new file mode 100644 index 0000000000..347f92e6f8 --- /dev/null +++ b/src/systems/state_machines.jl @@ -0,0 +1,155 @@ +_nameof(s) = nameof(s) +_nameof(s::Union{Int, Symbol}) = s +abstract type StateMachineOperator end +Base.broadcastable(x::StateMachineOperator) = Ref(x) +Symbolics.hide_lhs(_::StateMachineOperator) = true +struct InitialState <: StateMachineOperator + s::Any +end +Base.show(io::IO, s::InitialState) = print(io, "initial_state(", _nameof(s.s), ")") +initial_state(s) = Equation(InitialState(nothing), InitialState(s)) + +Base.@kwdef struct Transition{A, B, C} <: StateMachineOperator + from::A = nothing + to::B = nothing + cond::C = nothing + immediate::Bool = true + reset::Bool = true + synchronize::Bool = false + priority::Int = 1 + function Transition(from, to, cond, immediate, reset, synchronize, priority) + cond = unwrap(cond) + new{typeof(from), typeof(to), typeof(cond)}(from, to, cond, immediate, + reset, synchronize, + priority) + end +end +function Base.:(==)(transition1::Transition, transition2::Transition) + transition1.from == transition2.from && + transition1.to == transition2.to && + isequal(transition1.cond, transition2.cond) && + transition1.immediate == transition2.immediate && + transition1.reset == transition2.reset && + transition1.synchronize == transition2.synchronize && + transition1.priority == transition2.priority +end + +""" + transition(from, to, cond; immediate::Bool = true, reset::Bool = true, synchronize::Bool = false, priority::Int = 1) + +Create a transition from state `from` to state `to` that is enabled when transitioncondition `cond` evaluates to `true`. + +# Arguments: +- `from`: The source state of the transition. +- `to`: The target state of the transition. +- `cond`: A transition condition that evaluates to a Bool, such as `ticksInState() >= 2`. +- `immediate`: If `true`, the transition will fire at the same tick as it becomes true, if `false`, the actions of the state are evaluated first, and the transition fires during the next tick. +- `reset`: If true, the destination state `to` is reset to its initial condition when the transition fires. +- `synchronize`: If true, the transition will only fire if all sub-state machines in the source state are in their final (terminal) state. A final state is one that has no outgoing transitions. +- `priority`: If a state has more than one outgoing transition, all outgoing transitions must have a unique priority. The transitions are evaluated in priority order, i.e., the transition with priority 1 is evaluated first. +""" +function transition(from, to, cond; + immediate::Bool = true, reset::Bool = true, synchronize::Bool = false, + priority::Int = 1) + Equation( + Transition(), Transition(; from, to, cond, immediate, reset, + synchronize, priority)) +end +function Base.show(io::IO, s::Transition) + print(io, _nameof(s.from), " → ", _nameof(s.to), " if (", s.cond, ") [") + print(io, "immediate: ", Int(s.immediate), ", ") + print(io, "reset: ", Int(s.reset), ", ") + print(io, "sync: ", Int(s.synchronize), ", ") + print(io, "prio: ", s.priority, "]") +end + +function activeState end +function entry end +function ticksInState end +function timeInState end + +for (s, T) in [(:timeInState, :Real), + (:ticksInState, :Integer), + (:entry, :Bool), + (:activeState, :Bool)] + seed = hash(s) + @eval begin + $s(x) = wrap(term($s, x)) + SymbolicUtils.promote_symtype(::typeof($s), _...) = $T + function SymbolicUtils.show_call(io, ::typeof($s), args) + if isempty(args) + print(io, $s, "()") + else + arg = only(args) + print(io, $s, "(", arg isa Number ? arg : nameof(arg), ")") + end + end + end + if s != :activeState + @eval $s() = wrap(term($s)) + end +end + +@doc """ + timeInState() + timeInState(state) + +Get the time (in seconds) spent in a state in a finite state machine. + +When used to query the time spent in the enclosing state, the method without arguments is used, i.e., +```julia +@mtkmodel FSM begin + ... + @equations begin + var(k+1) ~ timeInState() >= 2 ? 0.0 : var(k) + end +end +``` + +If used to query the residence time of another state, the state is passed as an argument. + +This operator can be used in both equations and transition conditions. + +See also [`ticksInState`](@ref) and [`entry`](@ref) +""" timeInState + +@doc """ + ticksInState() + ticksInState(state) + +Get the number of ticks spent in a state in a finite state machine. + +When used to query the number of ticks spent in the enclosing state, the method without arguments is used, i.e., +```julia +@mtkmodel FSM begin + ... + @equations begin + var(k+1) ~ ticksInState() >= 2 ? 0.0 : var(k) + end +end +``` + +If used to query the number of ticks in another state, the state is passed as an argument. + +This operator can be used in both equations and transition conditions. + +See also [`timeInState`](@ref) and [`entry`](@ref) +""" ticksInState + +@doc """ + entry() + entry(state) + +When used in a finite-state machine, this operator returns true at the first tick when the state is active, and false otherwise. + +When used to query the entry of the enclosing state, the method without arguments is used, when used to query the entry of another state, the state is passed as an argument. + +This can be used to perform a unique action when entering a state. +""" +entry + +@doc """ + activeState(state) + +When used in a finite state machine, this operator returns `true` if the queried state is active and false otherwise. +""" activeState From a6235ea5dfdbd5eba36aaf957534c2ccc7aa1e80 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Mon, 26 May 2025 11:47:16 +0530 Subject: [PATCH 33/49] docs: remove old `BVProblem`, `SteadyStateProblem` docstrings --- src/problems/bvproblem.jl | 45 -------------------------------------- src/problems/odeproblem.jl | 14 ------------ 2 files changed, 59 deletions(-) diff --git a/src/problems/bvproblem.jl b/src/problems/bvproblem.jl index 30bd7e61c3..7040b32ffd 100644 --- a/src/problems/bvproblem.jl +++ b/src/problems/bvproblem.jl @@ -1,48 +1,3 @@ -""" -```julia -SciMLBase.BVProblem{iip}(sys::AbstractSystem, u0map, tspan, - parammap = DiffEqBase.NullParameters(); - constraints = nothing, guesses = nothing, - version = nothing, tgrad = false, - jac = true, sparse = true, - simplify = false, - kwargs...) where {iip} -``` - -Create a boundary value problem from the [`System`](@ref). - -`u0map` is used to specify fixed initial values for the states. Every variable -must have either an initial guess supplied using `guesses` or a fixed initial -value specified using `u0map`. - -Boundary value conditions are supplied to Systems in the form of a list of constraints. -These equations should specify values that state variables should take at specific points, -as in `x(0.5) ~ 1`). More general constraints that should hold over the entire solution, -such as `x(t)^2 + y(t)^2`, should be specified as one of the equations used to build the -`System`. - -If a `System` without `constraints` is specified, it will be treated as an initial value problem. - -```julia - @parameters g t_c = 0.5 - @variables x(..) y(t) λ(t) - eqs = [D(D(x(t))) ~ λ * x(t) - D(D(y)) ~ λ * y - g - x(t)^2 + y^2 ~ 1] - cstr = [x(0.5) ~ 1] - @mtkcompile pend = System(eqs, t; constraints = cstrs) - - tspan = (0.0, 1.5) - u0map = [x(t) => 0.6, y => 0.8] - parammap = [g => 1] - guesses = [λ => 1] - - bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(pend, u0map, tspan, parammap; guesses, check_length = false) -``` - -If the `System` has algebraic equations, like `x(t)^2 + y(t)^2`, the resulting -`BVProblem` must be solved using BVDAE solvers, such as Ascher. -""" @fallback_iip_specialize function SciMLBase.BVProblem{iip, spec}( sys::System, op, tspan; check_compatibility = true, cse = true, diff --git a/src/problems/odeproblem.jl b/src/problems/odeproblem.jl index e2901292ee..fbe7af5dec 100644 --- a/src/problems/odeproblem.jl +++ b/src/problems/odeproblem.jl @@ -83,20 +83,6 @@ end maybe_codegen_scimlproblem(expression, ODEProblem{iip}, args; kwargs...) end -""" -```julia -SciMLBase.SteadyStateProblem(sys::System, u0map, - parammap = DiffEqBase.NullParameters(); - version = nothing, tgrad = false, - jac = false, - checkbounds = false, sparse = false, - linenumbers = true, parallel = SerialForm(), - kwargs...) where {iip} -``` - -Generates an SteadyStateProblem from a `System` of ODEs and allows for automatically -symbolically calculating numerical enhancements. -""" @fallback_iip_specialize function DiffEqBase.SteadyStateProblem{iip, spec}( sys::System, op; check_length = true, check_compatibility = true, expression = Val{false}, kwargs...) where {iip, spec} From 228b2629f90b2e45dc451fbe3c3dd2a8dbbfde44 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Mon, 26 May 2025 11:47:34 +0530 Subject: [PATCH 34/49] fix: fix bugs in `HomotopyNonlinearFunction` --- .../nonlinear/homotopy_continuation.jl | 20 +++---------------- 1 file changed, 3 insertions(+), 17 deletions(-) diff --git a/src/systems/nonlinear/homotopy_continuation.jl b/src/systems/nonlinear/homotopy_continuation.jl index d74e4bd226..41ec50052d 100644 --- a/src/systems/nonlinear/homotopy_continuation.jl +++ b/src/systems/nonlinear/homotopy_continuation.jl @@ -472,22 +472,8 @@ function handle_rational_polynomials(x, wrt; fraction_cancel_fn = simplify_fract return num, den end -function SciMLBase.HomotopyNonlinearFunction(sys::System, args...; kwargs...) - ODEFunction{true}(sys, args...; kwargs...) -end - -function SciMLBase.HomotopyNonlinearFunction{true}(sys::System, args...; - kwargs...) - ODEFunction{true, SciMLBase.AutoSpecialize}(sys, args...; kwargs...) -end - -function SciMLBase.HomotopyNonlinearFunction{false}(sys::System, args...; - kwargs...) - ODEFunction{false, SciMLBase.FullSpecialize}(sys, args...; kwargs...) -end - -function SciMLBase.HomotopyNonlinearFunction{iip, specialize}( - sys::System, args...; eval_expression = false, eval_module = @__MODULE__, +@fallback_iip_specialize function SciMLBase.HomotopyNonlinearFunction{iip, specialize}( + sys::System; eval_expression = false, eval_module = @__MODULE__, p = nothing, fraction_cancel_fn = SymbolicUtils.simplify_fractions, cse = true, kwargs...) where {iip, specialize} if !iscomplete(sys) @@ -510,7 +496,7 @@ function SciMLBase.HomotopyNonlinearFunction{iip, specialize}( # we want to create f, jac etc. according to `sys2` since that will do the solving # but the `sys` inside for symbolic indexing should be the non-polynomial system - fn = NonlinearFunction{iip}(sys2; eval_expression, eval_module, cse, kwargs...) + fn = NonlinearFunction{iip}(sys2; p, eval_expression, eval_module, cse, kwargs...) obsfn = ObservedFunctionCache( sys; eval_expression, eval_module, checkbounds = get(kwargs, :checkbounds, false), cse) fn = remake(fn; sys = sys, observed = obsfn) From 298b694aa7b423f9bdd43e471eb24c9a7b54d0ac Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Mon, 26 May 2025 11:48:19 +0530 Subject: [PATCH 35/49] docs: add common documentation functionality for problems and functions --- src/ModelingToolkit.jl | 1 + src/problems/docs.jl | 393 +++++++++++++++++++++++++++++++++++++++++ 2 files changed, 394 insertions(+) create mode 100644 src/problems/docs.jl diff --git a/src/ModelingToolkit.jl b/src/ModelingToolkit.jl index 033d83d20b..306a1f5bc7 100644 --- a/src/ModelingToolkit.jl +++ b/src/ModelingToolkit.jl @@ -168,6 +168,7 @@ include("systems/imperative_affect.jl") include("systems/callbacks.jl") include("systems/system.jl") include("systems/codegen_utils.jl") +include("problems/docs.jl") include("systems/codegen.jl") include("systems/problem_utils.jl") include("linearization.jl") diff --git a/src/problems/docs.jl b/src/problems/docs.jl new file mode 100644 index 0000000000..94b77b8772 --- /dev/null +++ b/src/problems/docs.jl @@ -0,0 +1,393 @@ +const U0_P_DOCS = """ +The order of unknowns is determined by `unknowns(sys)`. If the system is split +[`is_split`](@ref) create an [`MTKParameters`](@ref) object. Otherwise, a parameter vector. +Initial values provided in terms of other variables will be symbolically evaluated. +The type of `op` will be used to determine the type of the containers. For example, if +given as an `SArray` of key-value pairs, `u0` will be an appropriately sized `SVector` +and the parameter object will be an `MTKParameters` object with `SArray`s inside. +""" + +const EVAL_EXPR_MOD_KWARGS = """ +- `eval_expression`: Whether to compile any functions via `eval` or + `RuntimeGeneratedFunctions`. +- `eval_module`: If `eval_expression == true`, the module to `eval` into. Otherwise, the + module in which to generate the `RuntimeGeneratedFunction`. +""" + +const INITIALIZEPROB_KWARGS = """ +- `guesses`: The guesses for variables in the system, used as initial values for the + initialization problem. +- `warn_initialize_determined`: Warn if the initialization system is under/over-determined. +- `initialization_eqs`: Extra equations to use in the initialization problem. +- `fully_determined`: Override whether the initialization system is fully determined. +- `use_scc`: Whether to use `SCCNonlinearProblem` for initialization if the system is fully + determined. +""" + +const PROBLEM_KWARGS = """ +$EVAL_EXPR_MOD_KWARGS +$INITIALIZEPROB_KWARGS +- `check_initialization_units`: Enable or disable unit checks when constructing the + initialization problem. +- `tofloat`: Passed to [`varmap_to_vars`](@ref) when building the parameter vector of + a non-split system. +- `u0_eltype`: The `eltype` of the `u0` vector. If `nothing`, finds the promoted floating point + type from `op`. +- `u0_constructor`: A function to apply to the `u0` value returned from + [`varmap_to_vars`](@ref). + to construct the final `u0` value. +- `p_constructor`: A function to apply to each array buffer created when constructing the + parameter object. +- `warn_cyclic_dependency`: Whether to emit a warning listing out cycles in initial + conditions provided for unknowns and parameters. +- `circular_dependency_max_cycle_length`: Maximum length of cycle to check for. Only + applicable if `warn_cyclic_dependency == true`. +- `circular_dependency_max_cycles`: Maximum number of cycles to check for. Only applicable + if `warn_cyclic_dependency == true`. +- `substitution_limit`: The number times to substitute initial conditions into each other + to attempt to arrive at a numeric value. +""" + +const TIME_DEPENDENT_PROBLEM_KWARGS = """ +- `callback`: An extra callback or `CallbackSet` to add to the problem, in addition to the + ones defined symbolically in the system. +""" + +const PROBLEM_INTERNALS_HEADER = """ +# Extended docs + +The following API is internal and may change or be removed without notice. Its usage is +highly discouraged. +""" + +const INTERNAL_INITIALIZEPROB_KWARGS = """ +- `time_dependent_init`: Whether to build a time-dependent initialization for the problem. A + time-dependent initialization solves for a consistent `u0`, whereas a time-independent one + only runs parameter initialization. +- `algebraic_only`: Whether to build the initialization problem using only algebraic equations. +- `allow_incomplete`: Whether to allow incomplete initialization problems. +""" + +const PROBLEM_INTERNAL_KWARGS = """ +- `build_initializeprob`: If `false`, avoids building the initialization problem. +- `check_length`: Whether to check the number of equations along with number of unknowns and + length of `u0` vector for consistency. If `false`, do not check with equations. This is + forwarded to `check_eqs_u0`. +$INTERNAL_INITIALIZEPROB_KWARGS +""" + +function problem_ctors(prob, istd) + if istd + """ + SciMLBase.$prob(sys::System, op, tspan::NTuple{2}; kwargs...) + SciMLBase.$prob{iip}(sys::System, op, tspan::NTuple{2}; kwargs...) + SciMLBase.$prob{iip, specialize}(sys::System, op, tspan::NTuple{2}; kwargs...) + """ + else + """ + SciMLBase.$prob(sys::System, op; kwargs...) + SciMLBase.$prob{iip}(sys::System, op; kwargs...) + SciMLBase.$prob{iip, specialize}(sys::System, op; kwargs...) + """ + end +end + +function prob_fun_common_kwargs(T, istd) + return """ + - `check_compatibility`: Whether to check if the given system `sys` contains all the + information necessary to create a `$T` and no more. If disabled, assumes that `sys` + at least contains the necessary information. + - `expression`: `Val{true}` to return an `Expr` that constructs the corresponding + problem instead of the problem itself. `Val{false}` otherwise. + $(istd ? " Constructing the expression does not support callbacks" : "") + """ +end + +function problem_docstring(prob, func, istd; init = true, extra_body = "") + if func isa DataType + func = "`$func`" + end + return """ + $(problem_ctors(prob, istd)) + + Build a `$prob` given a system `sys` and operating point `op` + $(istd ? " and timespan `tspan`" : ""). `iip` is a boolean indicating whether the + problem should be in-place. `specialization` is a `SciMLBase.AbstractSpecalize` subtype + indicating the level of specialization of the $func. The operating point should be an + iterable collection of key-value pairs mapping variables/parameters in the system to the + (initial) values they should take in `$prob`. Any values not provided will fallback to + the corresponding default (if present). + + $(init ? istd ? TIME_DEPENDENT_INIT : TIME_INDEPENDENT_INIT : "") + + $extra_body + + # Keyword arguments + + $PROBLEM_KWARGS + $(istd ? TIME_DEPENDENT_PROBLEM_KWARGS : "") + $(prob_fun_common_kwargs(prob, istd)) + + All other keyword arguments are forwarded to the $func constructor. + + $PROBLEM_INTERNALS_HEADER + + $PROBLEM_INTERNAL_KWARGS + """ +end + +const TIME_DEPENDENT_INIT = """ +ModelingToolkit will build an initialization problem where all initial values for +unknowns or observables of `sys` (either explicitly provided or in defaults) will +be constraints. To remove an initial condition in the defaults (without providing +a replacement) give the corresponding variable a value of `nothing` in the operating +point. The initialization problem will also run parameter initialization. See the +[Initialization](@ref initialization) documentation for more information. +""" + +const TIME_INDEPENDENT_INIT = """ +ModelingToolkit will build an initialization problem that will run parameter +initialization. Since it does not solve for initial values of unknowns, observed +equations will not be initialization constraints. If an initialization equation +of the system must involve the initial value of an unknown `x`, it must be used as +`Initial(x)` in the equation. For example, an equation to be used to solve for parameter +`p` in terms of unknowns `x` and `y` must be provided as `Initial(x) + Initial(y) ~ p` +instead of `x + y ~ p`. See the [Initialization](@ref initialization) documentation +for more information. +""" + +const BV_EXTRA_BODY = """ +Boundary value conditions are supplied to Systems in the form of a list of constraints. +These equations should specify values that state variables should take at specific points, +as in `x(0.5) ~ 1`). More general constraints that should hold over the entire solution, +such as `x(t)^2 + y(t)^2`, should be specified as one of the equations used to build the +`System`. + +If a `System` without `constraints` is specified, it will be treated as an initial value problem. + +```julia + @parameters g t_c = 0.5 + @variables x(..) y(t) λ(t) + eqs = [D(D(x(t))) ~ λ * x(t) + D(D(y)) ~ λ * y - g + x(t)^2 + y^2 ~ 1] + cstr = [x(0.5) ~ 1] + @mtkcompile pend = System(eqs, t; constraints = cstrs) + + tspan = (0.0, 1.5) + u0map = [x(t) => 0.6, y => 0.8] + parammap = [g => 1] + guesses = [λ => 1] + + bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(pend, u0map, tspan, parammap; guesses, check_length = false) +``` + +If the `System` has algebraic equations, like `x(t)^2 + y(t)^2`, the resulting +`BVProblem` must be solved using BVDAE solvers, such as Ascher. +""" + +for (mod, prob, func, istd, kws) in [ + (SciMLBase, :ODEProblem, ODEFunction, true, (;)), + (SciMLBase, :SteadyStateProblem, ODEFunction, false, (;)), + (SciMLBase, :BVProblem, ODEFunction, true, + (; init = false, extra_body = BV_EXTRA_BODY)), + (SciMLBase, :DAEProblem, DAEFunction, true, (;)), + (SciMLBase, :DDEProblem, DDEFunction, true, (;)), + (SciMLBase, :SDEProblem, SDEFunction, true, (;)), + (SciMLBase, :SDDEProblem, SDDEFunction, true, (;)), + (JumpProcesses, :JumpProblem, "inner SciMLFunction", true, (; init = false)), + (SciMLBase, :DiscreteProblem, DiscreteFunction, true, (;)), + (SciMLBase, :ImplicitDiscreteProblem, ImplicitDiscreteFunction, true, (;)), + (SciMLBase, :NonlinearProblem, NonlinearFunction, false, (;)), + (SciMLBase, :NonlinearLeastSquaresProblem, NonlinearFunction, false, (;)), + (SciMLBase, :SCCNonlinearProblem, NonlinearFunction, false, (; init = false)), + (SciMLBase, :OptimizationProblem, OptimizationFunction, false, (; init = false)) +] + @eval @doc problem_docstring($mod.$prob, $func, $istd) $mod.$prob +end + +function function_docstring(func, istd, optionals) + return """ + $func(sys::System; kwargs...) + $func{iip}(sys::System; kwargs...) + $func{iip, specialize}(sys::System; kwargs...) + + Create a `$func` from the given `sys`. `iip` is a boolean indicating whether the + function should be in-place. `specialization` is a `SciMLBase.AbstractSpecalize` + subtype indicating the level of specialization of the $func. + + # Keyword arguments + + - `u0`: The `u0` vector for the corresponding problem, if available. Can be obtained + using [`ModelingToolkit.get_u0`](@ref). + - `p`: The parameter object for the corresponding problem, if available. Can be obtained + using [`ModelingToolkit.get_p`](@ref). + $(istd ? TIME_DEPENDENT_FUNCTION_KWARGS : "") + $EVAL_EXPR_MOD_KWARGS + - `checkbounds`: Whether to enable bounds checking in the generated code. + - `simplify`: Whether to `simplify` any symbolically computed jacobians/hessians/etc. + - `cse`: Whether to enable Common Subexpression Elimination (CSE) on the generated code. + This typically improves performance of the generated code but reduces readability. + - `sparse`: Whether to generate jacobian/hessian/etc. functions that return/operate on + sparse matrices. Also controls whether the mass matrix is sparse, wherever applicable. + $(prob_fun_common_kwargs(func, istd)) + $(process_optional_function_kwargs(optionals)) + + All other keyword arguments are forwarded to the `$func` struct constructor. + """ +end + +const TIME_DEPENDENT_FUNCTION_KWARGS = """ +- `t`: The initial time for the corresponding problem, if available. +""" + +const JAC_KWARGS = """ +- `jac`: Whether to symbolically compute and generate code for the jacobian function. +""" + +const TGRAD_KWARGS = """ +- `tgrad`: Whether to symbolically compute and generate code for the `tgrad` function. +""" + +const SPARSITY_KWARGS = """ +- `sparsity`: Whether to provide symbolically compute and provide sparsity patterns for the + jacobian/hessian/etc. +""" + +const RESID_PROTOTYPE_KWARGS = """ +- `resid_prototype`: The prototype of the residual function `f` for a problem involving a + nonlinear solve where the residual and `u0` have different sizes. +""" + +const GRAD_KWARGS = """ +- `grad`: Whether the symbolically compute and generate code for the gradient of the cost + function with respect to unknowns. +""" + +const HESS_KWARGS = """ +- `hess`: Whether to symbolically compute and generate code for the hessian function. +""" + +const CONSH_KWARGS = """ +- `cons_h`: Whether to symbolically compute and generate code for the hessian function of + constraints. Since the constraint function is vector-valued, the hessian is a vector + of hessian matrices. +""" + +const CONSJ_KWARGS = """ +- `cons_j`: Whether to symbolically compute and generate code for the jacobian function of + constraints. +""" + +const CONSSPARSE_KWARGS = """ +- `cons_sparse`: Identical to the `sparse` keyword, but specifically for jacobian/hessian + functions of the constraints. +""" + +const INPUTFN_KWARGS = """ +- `inputs`: The variables in the input vector. The system must have been simplified using + `mtkcompile` with these variables passed as `inputs`. +- `disturbance_inputs`: The disturbance input variables. The system must have been + simplified using `mtkcompile` with these variables passed as `disturbance_inputs`. +""" + +const CONTROLJAC_KWARGS = """ +- `controljac`: Whether to symbolically compute and generate code for the jacobian of + the ODE with respect to the inputs. +""" + +const OPTIONAL_FN_KWARGS_DICT = Dict( + :jac => JAC_KWARGS, + :tgrad => TGRAD_KWARGS, + :sparsity => SPARSITY_KWARGS, + :resid_prototype => RESID_PROTOTYPE_KWARGS, + :grad => GRAD_KWARGS, + :hess => HESS_KWARGS, + :cons_h => CONSH_KWARGS, + :cons_j => CONSJ_KWARGS, + :cons_sparse => CONSSPARSE_KWARGS, + :inputfn => INPUTFN_KWARGS, + :controljac => CONTROLJAC_KWARGS +) + +const SPARSITY_OPTIONALS = Set([:jac, :hess, :cons_h, :cons_j, :controljac]) + +const CONS_SPARSITY_OPTIONALS = Set([:cons_h, :cons_j]) + +function process_optional_function_kwargs(choices::Vector{Symbol}) + if !isdisjoint(choices, SPARSITY_OPTIONALS) + push!(choices, :sparsity) + end + if !isdisjoint(choices, CONS_SPARSITY_OPTIONALS) + push!(choices, :cons_sparse) + end + join(map(Base.Fix1(getindex, OPTIONAL_FN_KWARGS_DICT), choices), "\n") +end + +for (mod, func, istd, optionals) in [ + (SciMLBase, :ODEFunction, true, [:jac, :tgrad]), + (SciMLBase, :ODEInputFunction, true, [:inputfn, :jac, :tgrad, :controljac]), + (SciMLBase, :DAEFunction, true, [:jac, :tgrad]), + (SciMLBase, :DDEFunction, true, Symbol[]), + (SciMLBase, :SDEFunction, true, [:jac, :tgrad]), + (SciMLBase, :SDDEFunction, true, Symbol[]), + (SciMLBase, :DiscreteFunction, true, Symbol[]), + (SciMLBase, :ImplicitDiscreteFunction, true, Symbol[]), + (SciMLBase, :NonlinearFunction, false, [:resid_prototype, :jac]), + (SciMLBase, :IntervalNonlinearFunction, false, Symbol[]), + (SciMLBase, :OptimizationFunction, false, [:jac, :grad, :hess, :cons_h, :cons_j]) +] + @eval @doc function_docstring($mod.$func, $istd, $optionals) $mod.$func +end + +@doc """ + SciMLBase.HomotopyNonlinearFunction(sys::System; kwargs...) + SciMLBase.HomotopyNonlinearFunction{iip}(sys::System; kwargs...) + SciMLBase.HomotopyNonlinearFunction{iip, specialize}(sys::System; kwargs...) + +Create a `HomotopyNonlinearFunction` from the given `sys`. `iip` is a boolean indicating +whether the function should be in-place. `specialization` is a `SciMLBase.AbstractSpecalize` +subtype indicating the level of specialization of the $func. + +# Keyword arguments + +- `u0`: The `u0` vector for the corresponding problem, if available. Can be obtained + using [`ModelingToolkit.get_u0`](@ref). +- `p`: The parameter object for the corresponding problem, if available. Can be obtained + using [`ModelingToolkit.get_p`](@ref). +$EVAL_EXPR_MOD_KWARGS +- `checkbounds`: Whether to enable bounds checking in the generated code. +- `simplify`: Whether to `simplify` any symbolically computed jacobians/hessians/etc. +- `cse`: Whether to enable Common Subexpression Elimination (CSE) on the generated code. + This typically improves performance of the generated code but reduces readability. +- `fraction_cancel_fn`: The function to use to simplify fractions in the polynomial + expression. A more powerful function can increase processing time but be able to + eliminate more rational functions, thus improving solve time. Should be a function that + takes a symbolic expression containing zero or more fraction expressions and returns the + simplified expression. While this defaults to `SymbolicUtils.simplify_fractions`, a viable + alternative is `SymbolicUtils.quick_cancel` + +All keyword arguments are forwarded to the wrapped `NonlinearFunction` constructor. +""" SciMLBase.HomotopyNonlinearFunction + +@doc """ + SciMLBase.IntervalNonlinearProblem(sys::System, uspan::NTuple{2}, parammap = SciMLBase.NullParameters(); kwargs...) + +Create an `IntervalNonlinearProblem` from the given `sys`. This is only valid for a system +of nonlinear equations with a single equation and unknown. `uspan` is the interval in which +the root is to be found, and `parammap` is an iterable collection of key-value pairs +providing values for the parameters in the system. + +$TIME_INDEPENDENT_INIT + +# Keyword arguments + +$PROBLEM_KWARGS +$(prob_fun_common_kwargs(IntervalNonlinearProblem, false)) + +All other keyword arguments are forwarded to the `IntervalNonlinearFunction` constructor. + +$PROBLEM_INTERNALS_HEADER + +$PROBLEM_INTERNAL_KWARGS +""" SciMLBase.IntervalNonlinearProblem From 36b49434a6ff49ccd5dfcd889c1d8029e3a2894a Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Mon, 26 May 2025 11:48:35 +0530 Subject: [PATCH 36/49] docs: use common docs in `process_SciMLProblem` docstring --- src/systems/problem_utils.jl | 52 +++++------------------------------- 1 file changed, 7 insertions(+), 45 deletions(-) diff --git a/src/systems/problem_utils.jl b/src/systems/problem_utils.jl index 76ade8e910..c1ca718b0a 100644 --- a/src/systems/problem_utils.jl +++ b/src/systems/problem_utils.jl @@ -1188,58 +1188,20 @@ Return the SciMLFunction created via calling `constructor`, the initial conditio and parameter object `p` given the system `sys`, and user-provided initial values `u0map` and `pmap`. `u0map` and `pmap` are converted into variable maps via [`to_varmap`](@ref). -The order of unknowns is determined by `unknowns(sys)`. If the system is split -[`is_split`](@ref) create an [`MTKParameters`](@ref) object. Otherwise, a parameter vector. -Initial values provided in terms of other variables will be symbolically evaluated using -[`evaluate_varmap!`](@ref). The type of `u0map` and `pmap` will be used to determine the -type of the containers (if parameters are not in an `MTKParameters` object). `Dict`s will be -turned into `Array`s. +$U0_P_DOCS This will also build the initialization problem and related objects and pass them to the SciMLFunction as keyword arguments. Keyword arguments: -- `build_initializeprob`: If `false`, avoids building the initialization problem. -- `t`: The initial time of the `ODEProblem`. If this is not provided, the initialization - problem cannot be built. +$PROBLEM_KWARGS +$PROBLEM_INTERNAL_KWARGS +- `t`: The initial time of the `SciMLProblem`. This does not need to be provided for time- + independent problems. If not provided for time-dependent problems, will be assumed as + zero. - `implicit_dae`: Also build a mapping of derivatives of states to values for implicit DAEs. - Changes the return value of this function to `(f, du0, u0, p)` instead of - `(f, u0, p)`. -- `guesses`: The guesses for variables in the system, used as initial values for the - initialization problem. -- `warn_initialize_determined`: Warn if the initialization system is under/over-determined. -- `initialization_eqs`: Extra equations to use in the initialization problem. -- `eval_expression`: Whether to compile any functions via `eval` or `RuntimeGeneratedFunctions`. -- `eval_module`: If `eval_expression == true`, the module to `eval` into. Otherwise, the module - in which to generate the `RuntimeGeneratedFunction`. -- `fully_determined`: Override whether the initialization system is fully determined. -- `check_initialization_units`: Enable or disable unit checks when constructing the - initialization problem. -- `tofloat`: Passed to [`varmap_to_vars`](@ref) when building the parameter vector of - a non-split system. -- `u0_eltype`: The `eltype` of the `u0` vector. If `nothing`, finds the promoted floating point - type from `op`. -- `u0_constructor`: A function to apply to the `u0` value returned from `varmap_to_vars` - to construct the final `u0` value. -- `p_constructor`: A function to apply to each array buffer created when constructing the parameter object. -- `check_length`: Whether to check the number of equations along with number of unknowns and - length of `u0` vector for consistency. If `false`, do not check with equations. This is - forwarded to `check_eqs_u0` + Changes the return value of this function to `(f, du0, u0, p)` instead of `(f, u0, p)`. - `symbolic_u0` allows the returned `u0` to be an array of symbolics. -- `warn_cyclic_dependency`: Whether to emit a warning listing out cycles in initial - conditions provided for unknowns and parameters. -- `circular_dependency_max_cycle_length`: Maximum length of cycle to check for. - Only applicable if `warn_cyclic_dependency == true`. -- `circular_dependency_max_cycles`: Maximum number of cycles to check for. - Only applicable if `warn_cyclic_dependency == true`. -- `substitution_limit`: The number times to substitute initial conditions into each - other to attempt to arrive at a numeric value. -- `use_scc`: Whether to use `SCCNonlinearProblem` for initialization if the system is fully - determined. -- `force_initialization_time_independent`: Whether to force the initialization to not use - the independent variable of `sys`. -- `algebraic_only`: Whether to build the initialization problem using only algebraic equations. -- `allow_incomplete`: Whether to allow incomplete initialization problems. All other keyword arguments are passed as-is to `constructor`. """ From 1897db8ad30f741a18cd627022f32657c1e919c0 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Mon, 26 May 2025 11:48:46 +0530 Subject: [PATCH 37/49] docs: tag the initialization docs --- docs/src/tutorials/initialization.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/tutorials/initialization.md b/docs/src/tutorials/initialization.md index cdf6e7e813..fedb8b9a22 100644 --- a/docs/src/tutorials/initialization.md +++ b/docs/src/tutorials/initialization.md @@ -1,4 +1,4 @@ -# Initialization of ODESystems +# [Initialization of Systems](@id initialization) While for simple numerical ODEs choosing an initial condition can be an easy affair, with ModelingToolkit's more general differential-algebraic equation From 1ff060aefa34aec087dc07129a5687075aba2b3b Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Mon, 26 May 2025 12:28:44 +0530 Subject: [PATCH 38/49] refactor: modernize `ODEInputFunction` --- src/systems/optimal_control_interface.jl | 59 ++++++------------------ 1 file changed, 13 insertions(+), 46 deletions(-) diff --git a/src/systems/optimal_control_interface.jl b/src/systems/optimal_control_interface.jl index c88aceabff..2b460f77a2 100644 --- a/src/systems/optimal_control_interface.jl +++ b/src/systems/optimal_control_interface.jl @@ -47,16 +47,10 @@ end is_explicit(tableau) = tableau isa DiffEqBase.ExplicitRKTableau -""" -Generate the control function f(x, u, p, t) from the ODESystem. -Input variables are automatically inferred but can be manually specified. -""" -function SciMLBase.ODEInputFunction{iip, specialize}(sys::System, - dvs = unknowns(sys), - ps = parameters(sys), u0 = nothing, +@fallback_iip_specialize function SciMLBase.ODEInputFunction{iip, specialize}(sys::System; inputs = unbound_inputs(sys), - disturbance_inputs = disturbances(sys); - version = nothing, tgrad = false, + disturbance_inputs = disturbances(sys), + u0 = nothing, tgrad = false, jac = false, controljac = false, p = nothing, t = nothing, eval_expression = false, @@ -66,7 +60,6 @@ function SciMLBase.ODEInputFunction{iip, specialize}(sys::System, checkbounds = false, sparsity = false, analytic = nothing, - split_idxs = nothing, initialization_data = nothing, cse = true, kwargs...) where {iip, specialize} @@ -75,61 +68,49 @@ function SciMLBase.ODEInputFunction{iip, specialize}(sys::System, f = f[1] if tgrad - tgrad_gen = generate_tgrad(sys, dvs, ps; + _tgrad = generate_tgrad(sys; simplify = simplify, expression = Val{true}, + wrap_gfw = Val{true}, expression_module = eval_module, cse, checkbounds = checkbounds, kwargs...) - tgrad_oop, tgrad_iip = eval_or_rgf.(tgrad_gen; eval_expression, eval_module) - _tgrad = GeneratedFunctionWrapper{(2, 3, is_split(sys))}(tgrad_oop, tgrad_iip) else _tgrad = nothing end if jac - jac_gen = generate_jacobian(sys, dvs, ps; + _jac = generate_jacobian(sys; simplify = simplify, sparse = sparse, expression = Val{true}, + wrap_gfw = Val{true}, expression_module = eval_module, cse, checkbounds = checkbounds, kwargs...) - jac_oop, jac_iip = eval_or_rgf.(jac_gen; eval_expression, eval_module) - - _jac = GeneratedFunctionWrapper{(2, 3, is_split(sys))}(jac_oop, jac_iip) else _jac = nothing end if controljac - cjac_gen = generate_control_jacobian(sys, dvs, ps; + _cjac = generate_control_jacobian(sys; simplify = simplify, sparse = sparse, - expression = Val{true}, + expression = Val{true}, wrap_gfw = Val{true}, expression_module = eval_module, cse, checkbounds = checkbounds, kwargs...) - cjac_oop, cjac_iip = eval_or_rgf.(cjac_gen; eval_expression, eval_module) - - _cjac = GeneratedFunctionWrapper{(2, 3, is_split(sys))}(cjac_oop, cjac_iip) else _cjac = nothing end M = calculate_massmatrix(sys) - _M = if sparse && !(u0 === nothing || M === I) - SparseArrays.sparse(M) - elseif u0 === nothing || M === I - M - else - ArrayInterface.restructure(u0 .* u0', M) - end + _M = concrete_massmatrix(M; sparse, u0) observedfun = ObservedFunctionCache( sys; steady_state, eval_expression, eval_module, checkbounds, cse) + _W_sparsity = W_sparsity(sys) + W_prototype = calculate_W_prototype(_W_sparsity; u0, sparse) if sparse uElType = u0 === nothing ? Float64 : eltype(u0) - W_prototype = similar(W_sparsity(sys), uElType) controljac_prototype = similar(calculate_control_jacobian(sys), uElType) else - W_prototype = nothing controljac_prototype = nothing end @@ -142,25 +123,11 @@ function SciMLBase.ODEInputFunction{iip, specialize}(sys::System, jac_prototype = W_prototype, controljac_prototype = controljac_prototype, observed = observedfun, - sparsity = sparsity ? W_sparsity(sys) : nothing, + sparsity = sparsity ? _W_sparsity : nothing, analytic = analytic, initialization_data) end -function SciMLBase.ODEInputFunction(sys::System, args...; kwargs...) - ODEInputFunction{true}(sys, args...; kwargs...) -end - -function SciMLBase.ODEInputFunction{true}(sys::System, args...; - kwargs...) - ODEInputFunction{true, SciMLBase.AutoSpecialize}(sys, args...; kwargs...) -end - -function SciMLBase.ODEInputFunction{false}(sys::System, args...; - kwargs...) - ODEInputFunction{false, SciMLBase.FullSpecialize}(sys, args...; kwargs...) -end - # returns the JuMP timespan, the number of steps, and whether it is a free time problem. function process_tspan(tspan, dt, steps) is_free_time = false From 6518c996928c831332f9d69b822b7550e1890f28 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Mon, 26 May 2025 12:52:19 +0530 Subject: [PATCH 39/49] refactor: remove old `force_initialization_time_independent` kwarg --- src/problems/initializationproblem.jl | 3 +-- src/problems/odeproblem.jl | 2 +- src/systems/problem_utils.jl | 6 +++--- 3 files changed, 5 insertions(+), 6 deletions(-) diff --git a/src/problems/initializationproblem.jl b/src/problems/initializationproblem.jl index e957fa7652..5d6c77b9a6 100644 --- a/src/problems/initializationproblem.jl +++ b/src/problems/initializationproblem.jl @@ -28,7 +28,6 @@ initial conditions for the given DAE. check_units = true, use_scc = true, allow_incomplete = false, - force_time_independent = false, algebraic_only = false, time_dependent_init = is_time_dependent(sys), kwargs...) where {iip, specialize} @@ -58,7 +57,7 @@ initial conditions for the given DAE. # useful for `SteadyStateProblem` since `f` has to be autonomous and the # initialization should be too - if force_time_independent + if !time_dependent_init idx = findfirst(isequal(get_iv(sys)), get_ps(isys)) idx === nothing || deleteat!(get_ps(isys), idx) end diff --git a/src/problems/odeproblem.jl b/src/problems/odeproblem.jl index fbe7af5dec..e2c8f9cd8a 100644 --- a/src/problems/odeproblem.jl +++ b/src/problems/odeproblem.jl @@ -91,7 +91,7 @@ end f, u0, p = process_SciMLProblem(ODEFunction{iip}, sys, op; steady_state = true, check_length, check_compatibility, expression, - force_initialization_time_independent = true, kwargs...) + time_dependent_init = false, kwargs...) kwargs = process_kwargs(sys; expression, kwargs...) args = (; f, u0, p) diff --git a/src/systems/problem_utils.jl b/src/systems/problem_utils.jl index c1ca718b0a..6eb44893e9 100644 --- a/src/systems/problem_utils.jl +++ b/src/systems/problem_utils.jl @@ -1217,7 +1217,7 @@ function process_SciMLProblem( circular_dependency_max_cycle_length = length(all_symbols(sys)), circular_dependency_max_cycles = 10, substitution_limit = 100, use_scc = true, time_dependent_init = is_time_dependent(sys), - force_initialization_time_independent = false, algebraic_only = false, + algebraic_only = false, allow_incomplete = false, is_initializeprob = false, kwargs...) dvs = unknowns(sys) ps = parameters(sys; initial_parameters = true) @@ -1276,8 +1276,8 @@ function process_SciMLProblem( eval_expression, eval_module, fully_determined, warn_cyclic_dependency, check_units = check_initialization_units, circular_dependency_max_cycle_length, circular_dependency_max_cycles, use_scc, - force_time_independent = force_initialization_time_independent, algebraic_only, allow_incomplete, - u0_constructor, p_constructor, floatT, time_dependent_init) + algebraic_only, allow_incomplete, u0_constructor, p_constructor, floatT, + time_dependent_init) kwargs = merge(kwargs, kws) end From bf24fdf8528296d94e5ef045fcd0b1b037ed714b Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Mon, 26 May 2025 12:52:37 +0530 Subject: [PATCH 40/49] docs: add docstring for `HomotopyContinuationProblem` --- src/systems/nonlinear/homotopy_continuation.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/systems/nonlinear/homotopy_continuation.jl b/src/systems/nonlinear/homotopy_continuation.jl index 41ec50052d..77e39e4839 100644 --- a/src/systems/nonlinear/homotopy_continuation.jl +++ b/src/systems/nonlinear/homotopy_continuation.jl @@ -514,6 +514,9 @@ end struct HomotopyContinuationProblem{iip, specialization} end +@doc problem_docstring( + HomotopyContinuationProblem, HomotopyNonlinearFunction, false; init = false) HomotopyContinuationProblem + function HomotopyContinuationProblem(sys::System, args...; kwargs...) HomotopyContinuationProblem{true}(sys, args...; kwargs...) end From 613a9b48285333dbd661fbc452dca8d06526c396 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Mon, 26 May 2025 12:52:54 +0530 Subject: [PATCH 41/49] docs: improve docstring for `InitializationProblem` --- src/problems/initializationproblem.jl | 34 +++++++++++++-------------- 1 file changed, 17 insertions(+), 17 deletions(-) diff --git a/src/problems/initializationproblem.jl b/src/problems/initializationproblem.jl index 5d6c77b9a6..8b30e449b0 100644 --- a/src/problems/initializationproblem.jl +++ b/src/problems/initializationproblem.jl @@ -1,22 +1,22 @@ struct InitializationProblem{iip, specialization} end -""" -```julia -InitializationProblem{iip}(sys::AbstractSystem, t, op; - version = nothing, tgrad = false, - jac = false, - checkbounds = false, sparse = false, - simplify = false, - linenumbers = true, parallel = SerialForm(), - initialization_eqs = [], - fully_determined = false, - kwargs...) where {iip} -``` - -Generates a NonlinearProblem or NonlinearLeastSquaresProblem from a System -which represents the initialization, i.e. the calculation of the consistent -initial conditions for the given DAE. -""" +@doc """ + InitializationProblem(sys::AbstractSystem, t, op = Dict(); kwargs...) + InitializationProblem{iip}(sys::AbstractSystem, t, op = Dict(); kwargs...) + InitializationProblem{iip, specialize}(sys::AbstractSystem, t, op = Dict(); kwargs...) + +Generate a `NonlinearProblem`, `SCCNonlinearProblem` or `NonlinearLeastSquaresProblem` to +represent a consistent initialization of `sys` given the initial time `t` and operating +point `op`. The initial time can be `nothing` for time-independent systems. + +# Keyword arguments + +$INITIALIZEPROB_KWARGS +$INTERNAL_INITIALIZEPROB_KWARGS + +All other keyword arguments are forwarded to the wrapped nonlinear problem constructor. +""" InitializationProblem + @fallback_iip_specialize function InitializationProblem{iip, specialize}( sys::AbstractSystem, t, op = Dict(); From 170cf53302518f2c28cf221ecfaa58430c0f2afd Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Mon, 26 May 2025 14:14:55 +0530 Subject: [PATCH 42/49] refactor: disallow passing `dvs`, `ps` to `generate_*` functions, add docstrings --- src/systems/codegen.jl | 353 ++++++++++++++++++++++++++++++++++++++--- 1 file changed, 327 insertions(+), 26 deletions(-) diff --git a/src/systems/codegen.jl b/src/systems/codegen.jl index 064be6fdfb..1c4e780e64 100644 --- a/src/systems/codegen.jl +++ b/src/systems/codegen.jl @@ -1,18 +1,38 @@ +const GENERATE_X_KWARGS = """ +- `expression`: `Val{true}` if this should return an `Expr` (or tuple of `Expr`s) of the + generated code. `Val{false}` otherwise. +- `wrap_gfw`: `Val{true}` if the returned functions should be wrapped in a callable + struct to make them callable using the expected syntax. The callable struct itself is + internal API. If `expression == Val{true}`, the returned expression will construct the + callable struct. If this function returns a tuple of functions/expressions, both will + be identical if `wrap_gfw == Val{true}`. +$EVAL_EXPR_MOD_KWARGS """ - $(TYPEDSIGNATURES) -Generate the RHS function for the `equations` of a `System`. +""" + $(TYPEDSIGNATURES) -# Arguments +Generate the RHS function for the [`equations`](@ref) of a [`System`](@ref). # Keyword Arguments +$GENERATE_X_KWARGS +- `implicit_dae`: Whether the generated function should be in the implicit form. Applicable + only for ODEs/DAEs or discrete systems. Instead of `f(u, p, t)` (`f(du, u, p, t)` for the + in-place form) the function is `f(du, u, p, t)` (respectively `f(resid, du, u, p, t)`). +- `override_discrete`: Whether to assume the system is discrete regardless of + `is_discrete_system(sys)`. +- `scalar`: Whether to generate a single-out-of-place function that returns a scalar for + the only equation in the system. + +All other keyword arguments are forwarded to [`build_function_wrapper`](@ref). """ -function generate_rhs(sys::System, dvs = unknowns(sys), - ps = parameters(sys; initial_parameters = true); implicit_dae = false, +function generate_rhs(sys::System; implicit_dae = false, scalar = false, expression = Val{true}, wrap_gfw = Val{false}, eval_expression = false, eval_module = @__MODULE__, override_discrete = false, kwargs...) + dvs = unknowns(sys) + ps = parameters(sys; initial_parameters = true) eqs = equations(sys) obs = observed(sys) u = dvs @@ -81,10 +101,22 @@ function generate_rhs(sys::System, dvs = unknowns(sys), res; eval_expression, eval_module) end -function generate_diffusion_function(sys::System, dvs = unknowns(sys), - ps = parameters(sys; initial_parameters = true); expression = Val{true}, +""" + $(TYPEDSIGNATURES) + +Generate the diffusion function for the noise equations of a [`System`](@ref). + +# Keyword Arguments + +$GENERATE_X_KWARGS + +All other keyword arguments are forwarded to [`build_function_wrapper`](@ref). +""" +function generate_diffusion_function(sys::System; expression = Val{true}, wrap_gfw = Val{false}, eval_expression = false, eval_module = @__MODULE__, kwargs...) + dvs = unknowns(sys) + ps = parameters(sys; initial_parameters = true) eqs = get_noise_eqs(sys) if ndims(eqs) == 2 && size(eqs, 2) == 1 # scalar noise @@ -106,6 +138,12 @@ function generate_diffusion_function(sys::System, dvs = unknowns(sys), expression, wrap_gfw, (p_start, nargs, is_split(sys)), res; eval_expression, eval_module) end +""" + $(TYPEDSIGNATURES) + +Calculate the gradient of the equations of `sys` with respect to the independent variable. +`simplify` is forwarded to `Symbolics.expand_derivatives`. +""" function calculate_tgrad(sys::System; simplify = false) # We need to remove explicit time dependence on the unknown because when we # have `u(t) * t` we want to have the tgrad to be `u(t)` instead of `u'(t) * @@ -121,6 +159,16 @@ function calculate_tgrad(sys::System; simplify = false) return tgrad end +""" + $(TYPEDSIGNATURES) + +Calculate the jacobian of the equations of `sys`. + +# Keyword arguments + +- `simplify`, `sparse`: Forwarded to `Symbolics.jacobian`. +- `dvs`: The variables with respect to which the jacobian should be computed. +""" function calculate_jacobian(sys::System; sparse = false, simplify = false, dvs = unknowns(sys)) obs = Dict(eq.lhs => eq.rhs for eq in observed(sys)) @@ -147,6 +195,18 @@ function calculate_jacobian(sys::System; return jac end +""" + $(TYPEDSIGNATURES) + +Generate the jacobian function for the equations of a [`System`](@ref). + +# Keyword Arguments + +$GENERATE_X_KWARGS +- `simplify`, `sparse`: Forwarded to [`calculate_jacobian`](@ref). + +All other keyword arguments are forwarded to [`build_function_wrapper`](@ref). +""" function generate_jacobian(sys::System; simplify = false, sparse = false, eval_expression = false, eval_module = @__MODULE__, expression = Val{true}, wrap_gfw = Val{false}, @@ -182,11 +242,24 @@ function assert_jac_length_header(sys) end end +""" + $(TYPEDSIGNATURES) + +Generate the tgrad function for the equations of a [`System`](@ref). + +# Keyword Arguments + +$GENERATE_X_KWARGS +- `simplify`: Forwarded to [`calculate_tgrad`](@ref). + +All other keyword arguments are forwarded to [`build_function_wrapper`](@ref). +""" function generate_tgrad( - sys::System, dvs = unknowns(sys), ps = parameters( - sys; initial_parameters = true); + sys::System; simplify = false, eval_expression = false, eval_module = @__MODULE__, expression = Val{true}, wrap_gfw = Val{false}, kwargs...) + dvs = unknowns(sys) + ps = parameters(sys; initial_parameters = true) tgrad = calculate_tgrad(sys, simplify = simplify) p = reorder_parameters(sys, ps) res = build_function_wrapper(sys, tgrad, @@ -215,6 +288,11 @@ function calculate_hessian(sys::System; simplify = false, sparse = false) return hess end +""" + $(TYPEDSIGNATURES) + +Return the sparsity pattern of the hessian of the equations of `sys`. +""" function Symbolics.hessian_sparsity(sys::System) hess = calculate_hessian(sys; sparse = true) return similar.(hess, Float64) @@ -222,10 +300,23 @@ end const W_GAMMA = only(@variables ˍ₋gamma) -function generate_W(sys::System, γ = 1.0, dvs = unknowns(sys), - ps = parameters(sys; initial_parameters = true); +""" + $(TYPEDSIGNATURES) + +Generate the `W = γ * M + J` function for the equations of a [`System`](@ref). + +# Keyword Arguments + +$GENERATE_X_KWARGS +- `simplify`, `sparse`: Forwarded to [`calculate_jacobian`](@ref). + +All other keyword arguments are forwarded to [`build_function_wrapper`](@ref). +""" +function generate_W(sys::System; simplify = false, sparse = false, expression = Val{true}, wrap_gfw = Val{false}, eval_expression = false, eval_module = @__MODULE__, kwargs...) + dvs = unknowns(sys) + ps = parameters(sys; initial_parameters = true) M = calculate_massmatrix(sys; simplify) if sparse M = SparseArrays.sparse(M) @@ -244,10 +335,25 @@ function generate_W(sys::System, γ = 1.0, dvs = unknowns(sys), expression, wrap_gfw, (2, 4, is_split(sys)), res; eval_expression, eval_module) end -function generate_dae_jacobian(sys::System, dvs = unknowns(sys), - ps = parameters(sys; initial_parameters = true); simplify = false, sparse = false, +""" + $(TYPEDSIGNATURES) + +Generate the DAE jacobian `γ * J′ + J` function for the equations of a [`System`](@ref). +`J′` is the jacobian of the equations with respect to the `du` vector, and `J` is the +standard jacobian. + +# Keyword Arguments + +$GENERATE_X_KWARGS +- `simplify`, `sparse`: Forwarded to [`calculate_jacobian`](@ref). + +All other keyword arguments are forwarded to [`build_function_wrapper`](@ref). +""" +function generate_dae_jacobian(sys::System; simplify = false, sparse = false, expression = Val{true}, wrap_gfw = Val{false}, eval_expression = false, eval_module = @__MODULE__, kwargs...) + dvs = unknowns(sys) + ps = parameters(sys; initial_parameters = true) jac_u = calculate_jacobian(sys; simplify = simplify, sparse = sparse) t = get_iv(sys) derivatives = Differential(t).(unknowns(sys)) @@ -262,6 +368,18 @@ function generate_dae_jacobian(sys::System, dvs = unknowns(sys), expression, wrap_gfw, (3, 5, is_split(sys)), res; eval_expression, eval_module) end +""" + $(TYPEDSIGNATURES) + +Generate the history function for a [`System`](@ref), given a symbolic representation of +the `u0` vector prior to the initial time. + +# Keyword Arguments + +$GENERATE_X_KWARGS + +All other keyword arguments are forwarded to [`build_function_wrapper`](@ref). +""" function generate_history(sys::System, u0; expression = Val{true}, wrap_gfw = Val{false}, eval_expression = false, eval_module = @__MODULE__, kwargs...) p = reorder_parameters(sys) @@ -272,6 +390,13 @@ function generate_history(sys::System, u0; expression = Val{true}, wrap_gfw = Va expression, wrap_gfw, (1, 2, is_split(sys)), res; eval_expression, eval_module) end +""" + $(TYPEDSIGNATURES) + +Calculate the mass matrix of `sys`. `simplify` controls whether `Symbolics.simplify` is +applied to the symbolic mass matrix. Returns a `Diagonal` or `LinearAlgebra.I` wherever +possible. +""" function calculate_massmatrix(sys::System; simplify = false) eqs = [eq for eq in equations(sys)] M = zeros(length(eqs), length(eqs)) @@ -285,7 +410,7 @@ function calculate_massmatrix(sys::System; simplify = false) error("Only semi-explicit constant mass matrices are currently supported. Faulty equation: $eq.") end end - M = simplify ? simplify.(M) : M + M = simplify ? Symbolics.simplify.(M) : M if isdiag(M) M = Diagonal(M) end @@ -293,6 +418,12 @@ function calculate_massmatrix(sys::System; simplify = false) M == I ? I : M end +""" + $(TYPEDSIGNATURES) + +Return a modified version of mass matrix `M` which is of a similar type to `u0`. `sparse` +controls whether the mass matrix should be a sparse matrix. +""" function concrete_massmatrix(M; sparse = false, u0 = nothing) if sparse && !(u0 === nothing || M === I) SparseArrays.sparse(M) @@ -305,6 +436,11 @@ function concrete_massmatrix(M; sparse = false, u0 = nothing) end end +""" + $(TYPEDSIGNATURES) + +Return the sparsity pattern of the jacobian of `sys` as a matrix. +""" function jacobian_sparsity(sys::System) sparsity = torn_system_jacobian_sparsity(sys) sparsity === nothing || return sparsity @@ -313,6 +449,13 @@ function jacobian_sparsity(sys::System) [dv for dv in unknowns(sys)]) end +""" + $(TYPEDSIGNATURES) + +Return the sparsity pattern of the DAE jacobian of `sys` as a matrix. + +See also: [`generate_dae_jacobian`](@ref). +""" function jacobian_dae_sparsity(sys::System) J1 = jacobian_sparsity([eq.rhs for eq in full_equations(sys)], [dv for dv in unknowns(sys)]) @@ -322,6 +465,13 @@ function jacobian_dae_sparsity(sys::System) J1 + J2 end +""" + $(TYPEDSIGNATURES) + +Return the sparsity pattern of the `W` matrix of `sys`. + +See also: [`generate_W`](@ref). +""" function W_sparsity(sys::System) jac_sparsity = jacobian_sparsity(sys) (n, n) = size(jac_sparsity) @@ -359,6 +509,18 @@ function get_constraint_unknown_subs!(subs::Dict, cons::Vector, stidxmap::Dict, end end +""" + $(TYPEDSIGNATURES) + +Generate the boundary condition function for a [`System`](@ref) given the state vector `u0`, +the indexes of `u0` to consider as hard constraints `u0_idxs` and the initial time `t0`. + +# Keyword Arguments + +$GENERATE_X_KWARGS + +All other keyword arguments are forwarded to [`build_function_wrapper`](@ref). +""" function generate_boundary_conditions(sys::System, u0, u0_idxs, t0; expression = Val{true}, wrap_gfw = Val{false}, eval_expression = false, eval_module = @__MODULE__, kwargs...) @@ -393,6 +555,17 @@ function generate_boundary_conditions(sys::System, u0, u0_idxs, t0; expression = expression, wrap_gfw, (2, 3, is_split(sys)), res; eval_expression, eval_module) end +""" + $(TYPEDSIGNATURES) + +Generate the cost function for a [`System`](@ref). + +# Keyword Arguments + +$GENERATE_X_KWARGS + +All other keyword arguments are forwarded to [`build_function_wrapper`](@ref). +""" function generate_cost(sys::System; expression = Val{true}, wrap_gfw = Val{false}, eval_expression = false, eval_module = @__MODULE__, kwargs...) obj = cost(sys) @@ -429,6 +602,18 @@ function calculate_cost_gradient(sys::System; simplify = false) return Symbolics.gradient(obj, dvs; simplify) end +""" + $(TYPEDSIGNATURES) + +Generate the gradient of the cost function with respect to unknowns for a [`System`](@ref). + +# Keyword Arguments + +$GENERATE_X_KWARGS +- `simplify`: Forwarded to [`calculate_cost_gradient`](@ref). + +All other keyword arguments are forwarded to [`build_function_wrapper`](@ref). +""" function generate_cost_gradient( sys::System; expression = Val{true}, wrap_gfw = Val{false}, eval_expression = false, eval_module = @__MODULE__, simplify = false, kwargs...) @@ -452,10 +637,29 @@ function calculate_cost_hessian(sys::System; sparse = false, simplify = false) end end +""" + $(TYPEDSIGNATURES) + +Return the sparsity pattern for the hessian of the cost function of `sys`. +""" function cost_hessian_sparsity(sys::System) return similar(calculate_cost_hessian(sys; sparse = true), Float64) end +""" + $(TYPEDSIGNATURES) + +Generate the hessian of the cost function for a [`System`](@ref). + +# Keyword Arguments + +$GENERATE_X_KWARGS +- `simplify`, `sparse`: Forwarded to [`calculate_cost_hessian`](@ref). +- `return_sparsity`: Whether to also return the sparsity pattern of the hessian as the + second return value. + +All other keyword arguments are forwarded to [`build_function_wrapper`](@ref). +""" function generate_cost_hessian( sys::System; expression = Val{true}, wrap_gfw = Val{false}, eval_expression = false, eval_module = @__MODULE__, simplify = false, @@ -481,6 +685,17 @@ function canonical_constraints(sys::System) end end +""" + $(TYPEDSIGNATURES) + +Generate the constraint function for a [`System`](@ref). + +# Keyword Arguments + +$GENERATE_X_KWARGS + +All other keyword arguments are forwarded to [`build_function_wrapper`](@ref). +""" function generate_cons(sys::System; expression = Val{true}, wrap_gfw = Val{false}, eval_expression = false, eval_module = @__MODULE__, kwargs...) cons = canonical_constraints(sys) @@ -491,13 +706,20 @@ function generate_cons(sys::System; expression = Val{true}, wrap_gfw = Val{false expression, wrap_gfw, (2, 2, is_split(sys)), res; eval_expression, eval_module) end -function generate_constraint_jacobian( - sys::System; expression = Val{true}, wrap_gfw = Val{false}, - eval_expression = false, eval_module = @__MODULE__, return_sparsity = false, - simplify = false, sparse = false, kwargs...) +""" + $(TYPEDSIGNATURES) + +Return the jacobian of the constraints of `sys` with respect to unknowns. + +# Keyword arguments + +- `simplify`, `sparse`: Forwarded to `Symbolics.jacobian`. +- `return_sparsity`: Whether to also return the sparsity pattern of the jacobian. +""" +function calculate_constraint_jacobian(sys::System; simplify = false, sparse = false, + return_sparsity = false) cons = canonical_constraints(sys) dvs = unknowns(sys) - ps = reorder_parameters(sys) sparsity = nothing if sparse jac = Symbolics.sparsejacobian(cons, dvs; simplify)::AbstractSparseArray @@ -505,19 +727,51 @@ function generate_constraint_jacobian( else jac = Symbolics.jacobian(cons, dvs; simplify) end + return return_sparsity ? (jac, sparsity) : jac +end + +""" + $(TYPEDSIGNATURES) + +Generate the jacobian of the constraint function for a [`System`](@ref). + +# Keyword Arguments + +$GENERATE_X_KWARGS +- `simplify`, `sparse`: Forwarded to [`calculate_constraint_jacobian`](@ref). +- `return_sparsity`: Whether to also return the sparsity pattern of the jacobian as the + second return value. + +All other keyword arguments are forwarded to [`build_function_wrapper`](@ref). +""" +function generate_constraint_jacobian( + sys::System; expression = Val{true}, wrap_gfw = Val{false}, + eval_expression = false, eval_module = @__MODULE__, return_sparsity = false, + simplify = false, sparse = false, kwargs...) + dvs = unknowns(sys) + ps = reorder_parameters(sys) + jac, sparsity = calculate_constraint_jacobian( + sys; simplify, sparse, return_sparsity = true) res = build_function_wrapper(sys, jac, dvs, ps...; expression = Val{true}, kwargs...) fn = maybe_compile_function( expression, wrap_gfw, (2, 2, is_split(sys)), res; eval_expression, eval_module) return return_sparsity ? (fn, sparsity) : fn end -function generate_constraint_hessian( - sys::System; expression = Val{true}, wrap_gfw = Val{false}, - eval_expression = false, eval_module = @__MODULE__, return_sparsity = false, - simplify = false, sparse = false, kwargs...) +""" + $(TYPEDSIGNATURES) + +Return the hessian of the constraints of `sys` with respect to unknowns. + +# Keyword arguments + +- `simplify`, `sparse`: Forwarded to `Symbolics.hessian`. +- `return_sparsity`: Whether to also return the sparsity pattern of the hessian. +""" +function calculate_constraint_hessian( + sys::System; simplify = false, sparse = false, return_sparsity = false) cons = canonical_constraints(sys) dvs = unknowns(sys) - ps = reorder_parameters(sys) sparsity = nothing if sparse hess = map(cons) do cstr @@ -527,12 +781,46 @@ function generate_constraint_hessian( else hess = [Symbolics.hessian(cstr, dvs; simplify) for cstr in cons] end + return return_sparsity ? (hess, sparsity) : hess +end + +""" + $(TYPEDSIGNATURES) + +Generate the hessian of the constraint function for a [`System`](@ref). + +# Keyword Arguments + +$GENERATE_X_KWARGS +- `simplify`, `sparse`: Forwarded to [`calculate_constraint_hessian`](@ref). +- `return_sparsity`: Whether to also return the sparsity pattern of the hessian as the + second return value. + +All other keyword arguments are forwarded to [`build_function_wrapper`](@ref). +""" +function generate_constraint_hessian( + sys::System; expression = Val{true}, wrap_gfw = Val{false}, + eval_expression = false, eval_module = @__MODULE__, return_sparsity = false, + simplify = false, sparse = false, kwargs...) + dvs = unknowns(sys) + ps = reorder_parameters(sys) + hess, sparsity = calculate_constraint_hessian( + sys; simplify, sparse, return_sparsity = true) res = build_function_wrapper(sys, hess, dvs, ps...; expression = Val{true}, kwargs...) fn = maybe_compile_function( expression, wrap_gfw, (2, 2, is_split(sys)), res; eval_expression, eval_module) return return_sparsity ? (fn, sparsity) : fn end +""" + $(TYPEDSIGNATURES) + +Calculate the jacobian of the equations of `sys` with respect to the inputs. + +# Keyword arguments + +- `simplify`, `sparse`: Forwarded to `Symbolics.jacobian`. +""" function calculate_control_jacobian(sys::AbstractSystem; sparse = false, simplify = false) rhs = [eq.rhs for eq in full_equations(sys)] @@ -547,10 +835,23 @@ function calculate_control_jacobian(sys::AbstractSystem; return jac end -function generate_control_jacobian(sys::AbstractSystem, dvs = unknowns(sys), - ps = parameters(sys; initial_parameters = true); +""" + $(TYPEDSIGNATURES) + +Generate the jacobian function of the equations of `sys` with respect to the inputs. + +# Keyword arguments + +$GENERATE_X_KWARGS +- `simplify`, `sparse`: Forwarded to [`calculate_constraint_hessian`](@ref). + +All other keyword arguments are forwarded to [`build_function_wrapper`](@ref). +""" +function generate_control_jacobian(sys::AbstractSystem; expression = Val{true}, wrap_gfw = Val{false}, eval_expression = false, eval_module = @__MODULE__, simplify = false, sparse = false, kwargs...) + dvs = unknowns(sys) + ps = parameters(sys; initial_parameters = true) jac = calculate_control_jacobian(sys; simplify = simplify, sparse = sparse) p = reorder_parameters(sys, ps) res = build_function_wrapper(sys, jac, dvs, p..., get_iv(sys); kwargs...) From 2f0acc1ff8b5d0c033bd68f9f4ff98e3f46e46f9 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Mon, 26 May 2025 14:15:25 +0530 Subject: [PATCH 43/49] refactor: remove old `calculate_*` and `generate_*` function stubs --- src/systems/abstractsystem.jl | 126 ---------------------------------- 1 file changed, 126 deletions(-) diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl index 40029086a1..33c4d72c00 100644 --- a/src/systems/abstractsystem.jl +++ b/src/systems/abstractsystem.jl @@ -8,132 +8,6 @@ end GUIMetadata(type) = GUIMetadata(type, nothing) -""" -```julia -calculate_tgrad(sys::AbstractTimeDependentSystem) -``` - -Calculate the time gradient of a system. - -Returns a vector of [`Num`](@ref) instances. The result from the first -call will be cached in the system object. -""" -function calculate_tgrad end - -""" -```julia -calculate_gradient(sys::AbstractSystem) -``` - -Calculate the gradient of a scalar system. - -Returns a vector of [`Num`](@ref) instances. The result from the first -call will be cached in the system object. -""" -function calculate_gradient end - -""" -```julia -calculate_jacobian(sys::AbstractSystem) -``` - -Calculate the Jacobian matrix of a system. - -Returns a matrix of [`Num`](@ref) instances. The result from the first -call will be cached in the system object. -""" -function calculate_jacobian end - -""" -```julia -calculate_control_jacobian(sys::AbstractSystem) -``` - -Calculate the Jacobian matrix of a system with respect to the system's controls. - -Returns a matrix of [`Num`](@ref) instances. The result from the first -call will be cached in the system object. -""" -function calculate_control_jacobian end - -""" -```julia -calculate_factorized_W(sys::AbstractSystem) -``` - -Calculate the factorized W-matrix of a system. - -Returns a matrix of [`Num`](@ref) instances. The result from the first -call will be cached in the system object. -""" -function calculate_factorized_W end - -""" -```julia -calculate_hessian(sys::AbstractSystem) -``` - -Calculate the hessian matrix of a scalar system. - -Returns a matrix of [`Num`](@ref) instances. The result from the first -call will be cached in the system object. -""" -function calculate_hessian end - -""" -```julia -generate_tgrad(sys::AbstractTimeDependentSystem, dvs = unknowns(sys), ps = parameters(sys), - expression = Val{true}; kwargs...) -``` - -Generates a function for the time gradient of a system. Extra arguments control -the arguments to the internal [`build_function`](@ref) call. -""" -function generate_tgrad end - -""" -```julia -generate_gradient(sys::AbstractSystem, dvs = unknowns(sys), ps = parameters(sys), - expression = Val{true}; kwargs...) -``` - -Generates a function for the gradient of a system. Extra arguments control -the arguments to the internal [`build_function`](@ref) call. -""" -function generate_gradient end - -""" -```julia -generate_jacobian(sys::AbstractSystem, dvs = unknowns(sys), ps = parameters(sys), - expression = Val{true}; sparse = false, kwargs...) -``` - -Generates a function for the Jacobian matrix of a system. Extra arguments control -the arguments to the internal [`build_function`](@ref) call. -""" -function generate_jacobian end - -""" -```julia -generate_hessian(sys::AbstractSystem, dvs = unknowns(sys), ps = parameters(sys), - expression = Val{true}; sparse = false, kwargs...) -``` - -Generates a function for the hessian matrix of a system. Extra arguments control -the arguments to the internal [`build_function`](@ref) call. -""" -function generate_hessian end - -""" -```julia -generate_function(sys::AbstractSystem, dvs = unknowns(sys), ps = parameters(sys), - expression = Val{true}; kwargs...) -``` - -Generate a function to evaluate the system's equations. -""" -function generate_rhs end - """ ```julia generate_custom_function(sys::AbstractSystem, exprs, dvs = unknowns(sys), From 8ed2e02d93b5caf6014119275a825ef6713c352e Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Mon, 26 May 2025 14:16:02 +0530 Subject: [PATCH 44/49] refactor: update from old `generate_*` syntax --- src/problems/daeproblem.jl | 6 ++---- src/problems/ddeproblem.jl | 5 +---- src/problems/discreteproblem.jl | 6 ++---- src/problems/implicitdiscreteproblem.jl | 3 +-- src/problems/intervalnonlinearproblem.jl | 4 +--- src/problems/nonlinearproblem.jl | 4 +--- src/problems/odeproblem.jl | 6 ++---- src/problems/optimizationproblem.jl | 3 +-- src/problems/sddeproblem.jl | 7 ++----- src/problems/sdeproblem.jl | 8 +++----- test/nonlinearsystem.jl | 4 ++-- test/odesystem.jl | 4 ++-- 12 files changed, 20 insertions(+), 40 deletions(-) diff --git a/src/problems/daeproblem.jl b/src/problems/daeproblem.jl index 860f2d2da2..6134923d4d 100644 --- a/src/problems/daeproblem.jl +++ b/src/problems/daeproblem.jl @@ -7,9 +7,7 @@ check_complete(sys, DAEFunction) check_compatibility && check_compatible_system(DAEFunction, sys) - dvs = unknowns(sys) - ps = parameters(sys) - f = generate_rhs(sys, dvs, ps; expression, wrap_gfw = Val{true}, + f = generate_rhs(sys; expression, wrap_gfw = Val{true}, implicit_dae = true, eval_expression, eval_module, checkbounds = checkbounds, cse, kwargs...) @@ -25,7 +23,7 @@ end if jac - _jac = generate_dae_jacobian(sys, dvs, ps; expression, + _jac = generate_dae_jacobian(sys; expression, wrap_gfw = Val{true}, simplify, sparse, cse, eval_expression, eval_module, checkbounds, kwargs...) else diff --git a/src/problems/ddeproblem.jl b/src/problems/ddeproblem.jl index 6f072b3d51..45af3edddb 100644 --- a/src/problems/ddeproblem.jl +++ b/src/problems/ddeproblem.jl @@ -6,10 +6,7 @@ check_complete(sys, DDEFunction) check_compatibility && check_compatible_system(DDEFunction, sys) - dvs = unknowns(sys) - ps = parameters(sys) - - f = generate_rhs(sys, dvs, ps; expression, wrap_gfw = Val{true}, + f = generate_rhs(sys; expression, wrap_gfw = Val{true}, eval_expression, eval_module, checkbounds = checkbounds, cse, kwargs...) diff --git a/src/problems/discreteproblem.jl b/src/problems/discreteproblem.jl index f640ee9ff5..ff0500ccf3 100644 --- a/src/problems/discreteproblem.jl +++ b/src/problems/discreteproblem.jl @@ -7,9 +7,7 @@ check_complete(sys, DiscreteFunction) check_compatibility && check_compatible_system(DiscreteFunction, sys) - dvs = unknowns(sys) - ps = parameters(sys) - f = generate_rhs(sys, dvs, ps; expression, wrap_gfw = Val{true}, + f = generate_rhs(sys; expression, wrap_gfw = Val{true}, eval_expression, eval_module, checkbounds = checkbounds, cse, kwargs...) @@ -45,7 +43,7 @@ end check_compatibility && check_compatible_system(DiscreteProblem, sys) dvs = unknowns(sys) - u0map = to_varmap(u0map, dvs) + u0map = to_varmap(op, dvs) add_toterms!(u0map; replace = true) f, u0, p = process_SciMLProblem(DiscreteFunction{iip, spec}, sys, op; t = tspan !== nothing ? tspan[1] : tspan, check_compatibility, expression, diff --git a/src/problems/implicitdiscreteproblem.jl b/src/problems/implicitdiscreteproblem.jl index 5b61e34df5..9e500d7354 100644 --- a/src/problems/implicitdiscreteproblem.jl +++ b/src/problems/implicitdiscreteproblem.jl @@ -9,8 +9,7 @@ iv = get_iv(sys) dvs = unknowns(sys) - ps = parameters(sys) - f = generate_rhs(sys, dvs, ps; expression, wrap_gfw = Val{true}, + f = generate_rhs(sys; expression, wrap_gfw = Val{true}, implicit_dae = true, eval_expression, eval_module, checkbounds = checkbounds, cse, override_discrete = true, kwargs...) diff --git a/src/problems/intervalnonlinearproblem.jl b/src/problems/intervalnonlinearproblem.jl index 65e7f7cb62..2dc929fe25 100644 --- a/src/problems/intervalnonlinearproblem.jl +++ b/src/problems/intervalnonlinearproblem.jl @@ -6,9 +6,7 @@ function SciMLBase.IntervalNonlinearFunction( check_complete(sys, IntervalNonlinearFunction) check_compatibility && check_compatible_system(IntervalNonlinearFunction, sys) - dvs = unknowns(sys) - ps = parameters(sys) - f = generate_rhs(sys, dvs, ps; expression, wrap_gfw = Val{true}, + f = generate_rhs(sys; expression, wrap_gfw = Val{true}, scalar = true, eval_expression, eval_module, checkbounds, cse, kwargs...) observedfun = ObservedFunctionCache( diff --git a/src/problems/nonlinearproblem.jl b/src/problems/nonlinearproblem.jl index ef36179e06..9986ee665b 100644 --- a/src/problems/nonlinearproblem.jl +++ b/src/problems/nonlinearproblem.jl @@ -8,9 +8,7 @@ check_complete(sys, NonlinearFunction) check_compatibility && check_compatible_system(NonlinearFunction, sys) - dvs = unknowns(sys) - ps = parameters(sys) - f = generate_rhs(sys, dvs, ps; expression, wrap_gfw = Val{true}, + f = generate_rhs(sys; expression, wrap_gfw = Val{true}, eval_expression, eval_module, checkbounds = checkbounds, cse, kwargs...) diff --git a/src/problems/odeproblem.jl b/src/problems/odeproblem.jl index e2c8f9cd8a..e78ff00525 100644 --- a/src/problems/odeproblem.jl +++ b/src/problems/odeproblem.jl @@ -7,9 +7,7 @@ check_complete(sys, ODEFunction) check_compatibility && check_compatible_system(ODEFunction, sys) - dvs = unknowns(sys) - ps = parameters(sys) - f = generate_rhs(sys, dvs, ps; expression, wrap_gfw = Val{true}, + f = generate_rhs(sys; expression, wrap_gfw = Val{true}, eval_expression, eval_module, checkbounds = checkbounds, cse, kwargs...) @@ -26,7 +24,7 @@ if tgrad _tgrad = generate_tgrad( - sys, dvs, ps; expression, wrap_gfw = Val{true}, + sys; expression, wrap_gfw = Val{true}, simplify, cse, eval_expression, eval_module, checkbounds, kwargs...) else _tgrad = nothing diff --git a/src/problems/optimizationproblem.jl b/src/problems/optimizationproblem.jl index 6d5fb1d79e..36b34a8867 100644 --- a/src/problems/optimizationproblem.jl +++ b/src/problems/optimizationproblem.jl @@ -10,8 +10,7 @@ function SciMLBase.OptimizationFunction{iip}(sys::System; expression = Val{false}, kwargs...) where {iip} check_complete(sys, OptimizationFunction) check_compatibility && check_compatible_system(OptimizationFunction, sys) - dvs = unknowns(sys) - ps = parameters(sys) + cstr = constraints(sys) f = generate_cost(sys; expression, wrap_gfw = Val{true}, eval_expression, diff --git a/src/problems/sddeproblem.jl b/src/problems/sddeproblem.jl index 0daa04fcf7..e1cc00b2a7 100644 --- a/src/problems/sddeproblem.jl +++ b/src/problems/sddeproblem.jl @@ -6,12 +6,9 @@ check_complete(sys, SDDEFunction) check_compatibility && check_compatible_system(SDDEFunction, sys) - dvs = unknowns(sys) - ps = parameters(sys) - - f = generate_rhs(sys, dvs, ps; expression, wrap_gfw = Val{true}, + f = generate_rhs(sys; expression, wrap_gfw = Val{true}, eval_expression, eval_module, checkbounds = checkbounds, cse, kwargs...) - g = generate_diffusion_function(sys, dvs, ps; expression, + g = generate_diffusion_function(sys; expression, wrap_gfw = Val{true}, eval_expression, eval_module, checkbounds, cse, kwargs...) if spec === SciMLBase.FunctionWrapperSpecialize && iip diff --git a/src/problems/sdeproblem.jl b/src/problems/sdeproblem.jl index 57ca120641..1bc47118ff 100644 --- a/src/problems/sdeproblem.jl +++ b/src/problems/sdeproblem.jl @@ -7,12 +7,10 @@ check_complete(sys, SDEFunction) check_compatibility && check_compatible_system(SDEFunction, sys) - dvs = unknowns(sys) - ps = parameters(sys) - f = generate_rhs(sys, dvs, ps; expression, wrap_gfw = Val{true}, + f = generate_rhs(sys; expression, wrap_gfw = Val{true}, eval_expression, eval_module, checkbounds = checkbounds, cse, kwargs...) - g = generate_diffusion_function(sys, dvs, ps; expression, + g = generate_diffusion_function(sys; expression, wrap_gfw = Val{true}, eval_expression, eval_module, checkbounds, cse, kwargs...) if spec === SciMLBase.FunctionWrapperSpecialize && iip @@ -27,7 +25,7 @@ end if tgrad - _tgrad = generate_tgrad(sys, dvs, ps; expression, + _tgrad = generate_tgrad(sys; expression, wrap_gfw = Val{true}, simplify, cse, eval_expression, eval_module, checkbounds, kwargs...) else diff --git a/test/nonlinearsystem.jl b/test/nonlinearsystem.jl index bf05e317ce..e4a5b418b9 100644 --- a/test/nonlinearsystem.jl +++ b/test/nonlinearsystem.jl @@ -29,7 +29,7 @@ eqs = [0 ~ σ * (y - x) * h, @test eval(toexpr(ns)) == ns test_nlsys_inference("standard", ns, (x, y, z), (σ, ρ, β, h)) @test begin - f = generate_rhs(ns, [x, y, z], [σ, ρ, β, h], expression = Val{false})[2] + f = generate_rhs(ns, expression = Val{false})[2] du = [0.0, 0.0, 0.0] f(du, [1, 2, 3], [1, 2, 3, 1]) du ≈ [1, -3, -7] @@ -65,7 +65,7 @@ eqs = [0 ~ σ * a * h, 0 ~ x * y - β * z] @named ns = System(eqs, [x, y, z], [σ, ρ, β, h]) ns = complete(ns) -nlsys_func = generate_rhs(ns, [x, y, z], [σ, ρ, β, h]) +nlsys_func = generate_rhs(ns) nf = NonlinearFunction(ns) jac = calculate_jacobian(ns) diff --git a/test/odesystem.jl b/test/odesystem.jl index 6dc8b06eb9..7bf4dfebb0 100644 --- a/test/odesystem.jl +++ b/test/odesystem.jl @@ -148,7 +148,7 @@ f(du, [1.0, 2.0, 3.0], [x -> x + 7, 2, 3, 1], 5.0) eqs = [D(x) ~ x + 10σ(t - 1) + 100σ(t - 2) + 1000σ(t^2)] @named de = System(eqs, t) test_diffeq_inference("many internal iv-varying", de, t, (x,), (σ,)) -f = generate_rhs(de, [x], [σ], expression = Val{false}, wrap_gfw = Val{true}) +f = generate_rhs(de, expression = Val{false}, wrap_gfw = Val{true}) du = [0.0] f(du, [1.0], [t -> t + 2], 5.0) @test du ≈ [27561] @@ -173,7 +173,7 @@ eqs = [D(x) ~ -A * x, @named de = System(eqs, t) @test begin local f, du - f = generate_rhs(de, [x, y], [A, B, C], expression = Val{false}, wrap_gfw = Val{true}) + f = generate_rhs(de, expression = Val{false}, wrap_gfw = Val{true}) du = [0.0, 0.0] f(du, [1.0, 2.0], [1, 2, 3], 0.0) du ≈ [-1, -1 / 3] From 6874e2b7b2757f4f0532d4f92962ffc37078a1d3 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Mon, 26 May 2025 14:23:46 +0530 Subject: [PATCH 45/49] docs: add new doc pages, remove outdated pages --- docs/pages.jl | 18 +++++++----------- 1 file changed, 7 insertions(+), 11 deletions(-) diff --git a/docs/pages.jl b/docs/pages.jl index 0177947ed2..6c0806b15a 100644 --- a/docs/pages.jl +++ b/docs/pages.jl @@ -27,9 +27,13 @@ pages = [ "Advanced Examples" => Any["examples/tearing_parallelism.md", "examples/sparse_jacobians.md", "examples/perturbation.md"]], - "Basics" => Any["basics/AbstractSystem.md", - "basics/ContextualVariables.md", - "basics/Variable_metadata.md", + "API" => Any["API/System.md", + "API/variables.md", + "API/model_building.md", + "API/problems.md", + "API/codegen.md", + "API/PDESystem.md"], + "Basics" => Any[ "basics/Composition.md", "basics/Events.md", "basics/Linearization.md", @@ -40,14 +44,6 @@ pages = [ "basics/DependencyGraphs.md", "basics/Precompilation.md", "basics/FAQ.md"], - "System Types" => Any["systems/ODESystem.md", - "systems/SDESystem.md", - "systems/JumpSystem.md", - "systems/NonlinearSystem.md", - "systems/OptimizationSystem.md", - "systems/PDESystem.md", - "systems/DiscreteSystem.md", - "systems/ImplicitDiscreteSystem.md"], "comparison.md", "internals.md" ] From 175b9cd4ce10a01a4ba098ed3189a23afd480746 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Mon, 26 May 2025 14:24:18 +0530 Subject: [PATCH 46/49] docs: add sources to `docs/Project.toml` to enable doc builds --- docs/Project.toml | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/docs/Project.toml b/docs/Project.toml index ae0d89dd3a..f9cf903ed7 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -17,6 +17,7 @@ ModelingToolkitStandardLibrary = "16a59e39-deab-5bd0-87e4-056b12336739" NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec" Optim = "429524aa-4258-5aef-a3af-852621145aeb" Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba" +OptimizationBase = "bca83a33-5cc9-4baa-983d-23429ab6bcbb" OptimizationOptimJL = "36348300-93cb-4f02-beb5-3c3902f8871e" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" @@ -29,6 +30,10 @@ SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b" Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" +[sources] +ModelingToolkitStandardLibrary = {rev = "mtk-v10", url = "https://github.com/SciML/ModelingToolkitStandardLibrary.jl/"} +OptimizationBase = {rev = "as/mtk-v10", url = "https://github.com/AayushSabharwal/OptimizationBase.jl"} + [compat] Attractors = "1.24" BenchmarkTools = "1.3" From 276aef366b83c6d22f53356bc47eb2194a44ea91 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Tue, 27 May 2025 12:14:11 +0530 Subject: [PATCH 47/49] fix: use `time_dependent_init` instead of `is_time_dependent(sys)` where applicable --- src/systems/problem_utils.jl | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/src/systems/problem_utils.jl b/src/systems/problem_utils.jl index 6eb44893e9..2641d92108 100644 --- a/src/systems/problem_utils.jl +++ b/src/systems/problem_utils.jl @@ -1080,7 +1080,7 @@ function maybe_build_initialization_problem( end initializeprob = remake(initializeprob; p = initp) - get_initial_unknowns = if is_time_dependent(sys) + get_initial_unknowns = if time_dependent_init GetUpdatedU0(sys, initializeprob.f.sys, op) else nothing @@ -1092,9 +1092,6 @@ function maybe_build_initialization_problem( sys, initializeprob.f.sys; u0_constructor, p_constructor), get_initial_unknowns, SetInitialUnknowns(sys)) - if time_dependent_init === nothing - time_dependent_init = is_time_dependent(sys) - end if time_dependent_init all_init_syms = Set(all_symbols(initializeprob)) solved_unknowns = filter(var -> var in all_init_syms, unknowns(sys)) From 03dc7992ded37b48d07f08e04a1cf98ff8bfda27 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Tue, 27 May 2025 12:14:41 +0530 Subject: [PATCH 48/49] feat: implement latexification for `Connector` --- src/systems/analysis_points.jl | 2 ++ src/systems/connectors.jl | 9 +++++++++ 2 files changed, 11 insertions(+) diff --git a/src/systems/analysis_points.jl b/src/systems/analysis_points.jl index 12185c943f..bbd23b9a50 100644 --- a/src/systems/analysis_points.jl +++ b/src/systems/analysis_points.jl @@ -140,6 +140,8 @@ function Base.show(io::IO, ::MIME"text/plain", ap::AnalysisPoint) end end +Symbolics.hide_lhs(::AnalysisPoint) = true + @latexrecipe function f(ap::AnalysisPoint) index --> :subscript snakecase --> true diff --git a/src/systems/connectors.jl b/src/systems/connectors.jl index 9761d18493..d771445074 100644 --- a/src/systems/connectors.jl +++ b/src/systems/connectors.jl @@ -41,6 +41,15 @@ function Base.show(io::IO, c::Connection) print(io, ")") end +@latexrecipe function f(c::Connection) + index --> :subscript + return Expr(:call, :connect, map(nameof, c.systems)...) +end + +function Base.show(io::IO, ::MIME"text/latex", ap::Connection) + print(io, latexify(ap)) +end + isconnection(_) = false isconnection(_::Connection) = true """ From 1daaa32e2c9d1f2936fae30bbc769b123cffc614 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Tue, 27 May 2025 22:34:13 +0530 Subject: [PATCH 49/49] build: bump Symbolics compat --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 527451bca2..e77aa6377c 100644 --- a/Project.toml +++ b/Project.toml @@ -158,7 +158,7 @@ StochasticDelayDiffEq = "1.10" StochasticDiffEq = "6.72.1" SymbolicIndexingInterface = "0.3.39" SymbolicUtils = "3.26.1" -Symbolics = "6.37" +Symbolics = "6.40" URIs = "1" UnPack = "0.1, 1.0" Unitful = "1.1"