post-merge

Author: vyudu

.github/workflows/Downstream.yml

@@ -37,7 +37,6 @@ jobs:
           - {user: SciML, repo: MethodOfLines.jl, group: Interface}
           - {user: SciML, repo: MethodOfLines.jl, group: 2D_Diffusion}
           - {user: SciML, repo: MethodOfLines.jl, group: DAE}
-          - {user: ai4energy, repo: Ai4EComponentLib.jl, group: Downstream}
     steps:
       - uses: actions/checkout@v4
       - uses: julia-actions/setup-julia@v1

.github/workflows/Tests.yml

@@ -38,6 +38,7 @@ jobs:
           - Extensions
           - Downstream
           - RegressionI
+          - FMI
     uses: "SciML/.github/.github/workflows/tests.yml@v1"
     with:
       julia-version: "${{ matrix.version }}"

Project.toml

@@ -1,7 +1,7 @@
 name = "ModelingToolkit"
 uuid = "961ee093-0014-501f-94e3-6117800e7a78"
 authors = ["Yingbo Ma <[email protected]>", "Chris Rackauckas <[email protected]> and contributors"]
-version = "9.61.0"
+version = "9.64.0"
 
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
@@ -64,6 +64,7 @@ Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 BifurcationKit = "0f109fa4-8a5d-4b75-95aa-f515264e7665"
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
 DeepDiffs = "ab62b9b5-e342-54a8-a765-a90f495de1a6"
+FMI = "14a09403-18e3-468f-ad8a-74f8dda2d9ac"
 HomotopyContinuation = "f213a82b-91d6-5c5d-acf7-10f1c761b327"
 InfiniteOpt = "20393b10-9daf-11e9-18c9-8db751c92c57"
 LabelledArrays = "2ee39098-c373-598a-b85f-a56591580800"
@@ -72,6 +73,7 @@ LabelledArrays = "2ee39098-c373-598a-b85f-a56591580800"
 MTKBifurcationKitExt = "BifurcationKit"
 MTKChainRulesCoreExt = "ChainRulesCore"
 MTKDeepDiffsExt = "DeepDiffs"
+MTKFMIExt = "FMI"
 MTKHomotopyContinuationExt = "HomotopyContinuation"
 MTKInfiniteOptExt = "InfiniteOpt"
 MTKLabelledArraysExt = "LabelledArrays"
@@ -106,6 +108,7 @@ FindFirstFunctions = "1"
 ForwardDiff = "0.10.3"
 FunctionWrappers = "1.1"
 FunctionWrappersWrappers = "0.1"
+FMI = "0.14"
 Graphs = "1.5.2"
 HomotopyContinuation = "2.11"
 InfiniteOpt = "0.5"
@@ -132,7 +135,7 @@ RecursiveArrayTools = "3.26"
 Reexport = "0.2, 1"
 RuntimeGeneratedFunctions = "0.5.9"
 SCCNonlinearSolve = "1.0.0"
-SciMLBase = "2.71.1"
+SciMLBase = "2.73"
 SciMLStructures = "1.0"
 Serialization = "1"
 Setfield = "0.7, 0.8, 1"
@@ -142,9 +145,9 @@ SpecialFunctions = "0.7, 0.8, 0.9, 0.10, 1.0, 2"
 StaticArrays = "0.10, 0.11, 0.12, 1.0"
 StochasticDiffEq = "6.72.1"
 StochasticDelayDiffEq = "1.8.1"
-SymbolicIndexingInterface = "0.3.36"
-SymbolicUtils = "3.10"
-Symbolics = "6.22.1"
+SymbolicIndexingInterface = "0.3.37"
+SymbolicUtils = "3.14"
+Symbolics = "6.29"
 URIs = "1"
 UnPack = "0.1, 1.0"
 Unitful = "1.1"

docs/Project.toml

@@ -15,6 +15,7 @@ Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba"
 OptimizationOptimJL = "36348300-93cb-4f02-beb5-3c3902f8871e"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+PreallocationTools = "d236fae5-4411-538c-8e31-a6e3d9e00b46"
 SciMLStructures = "53ae85a6-f571-4167-b2af-e1d143709226"
 Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46"
 StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0"
@@ -38,6 +39,7 @@ Optimization = "3.9, 4"
 OptimizationOptimJL = "0.1, 0.4"
 OrdinaryDiffEq = "6.31"
 Plots = "1.36"
+PreallocationTools = "0.4"
 SciMLStructures = "1.1"
 Setfield = "1"
 StochasticDiffEq = "6"

docs/pages.jl

@@ -41,7 +41,8 @@ pages = [
         "systems/JumpSystem.md",
         "systems/NonlinearSystem.md",
         "systems/OptimizationSystem.md",
-        "systems/PDESystem.md"],
+        "systems/PDESystem.md",
+        "systems/DiscreteSystem.md"],
     "comparison.md",
     "internals.md"
 ]

docs/src/basics/Debugging.md

@@ -35,6 +35,39 @@ dsol = solve(dprob, Tsit5());
 Now we see that it crashed because `u1` decreased so much that it became negative and outside the domain of the `√` function.
 We could have figured that out ourselves, but it is not always so obvious for more complex models.
 
+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!"])
+```
+
+The assertions must be an iterable of pairs, where the first element is the symbolic condition and
+the second is a message to be logged when the condition fails. All assertions are added to the
+generated code and will cause the solver to reject steps that fail the assertions. For systems such
+as the above where the assertion is guaranteed to eventually fail, the solver will likely exit
+with a `dtmin` failure..
+
+```@example debug
+prob = ODEProblem(sys, [], (0.0, 10.0))
+sol = solve(prob, Tsit5())
+```
+
+We can use `debug_system` to log the failing assertions in each call to the RHS function.
+
+```@repl debug
+dsys = debug_system(sys; functions = []);
+dprob = ODEProblem(dsys, [], (0.0, 10.0));
+dsol = solve(dprob, Tsit5());
+```
+
+Note the logs containing the failed assertion and corresponding message. To temporarily disable
+logging in a system returned from `debug_system`, use `ModelingToolkit.ASSERTION_LOG_VARIABLE`.
+
+```@repl debug
+dprob[ModelingToolkit.ASSERTION_LOG_VARIABLE] = false;
+solve(drob, Tsit5());
+```
+
 ```@docs
 debug_system
 ```

docs/src/basics/MTKLanguage.md

@@ -23,6 +23,7 @@ equations.
 
 `@mtkmodel` definition contains begin blocks of
 
+  - `@description`: for describing the whole system with a human-readable string
   - `@components`: for listing sub-components of the system
   - `@constants`: for declaring constants
   - `@defaults`: for passing `defaults` to ODESystem
@@ -42,20 +43,23 @@ using ModelingToolkit
 using ModelingToolkit: t
 
 @mtkmodel ModelA begin
+    @description "A component with parameters `k` and `k_array`."
     @parameters begin
         k
         k_array[1:2]
     end
 end
 
 @mtkmodel ModelB begin
+    @description "A component with parameters `p1` and `p2`."
     @parameters begin
         p1 = 1.0, [description = "Parameter of ModelB"]
         p2 = 1.0, [description = "Parameter of ModelB"]
     end
 end
 
 @mtkmodel ModelC begin
+    @description "A bigger system that contains many more things."
     @icon "https://github.com/SciML/SciMLDocs/blob/main/docs/src/assets/logo.png"
     @constants begin
         c::Int = 1, [description = "Example constant."]
@@ -91,6 +95,10 @@ end
 end
 ```
 
+#### `@description`
+
+A documenting `String` that summarizes and explains what the model is.
+
 #### `@icon`
 
 An icon can be embedded in 3 ways:

docs/src/basics/Variable_metadata.md

@@ -54,6 +54,11 @@ using ModelingToolkit: t_nounits as t, D_nounits as D
 
 @variables i(t) [connect = Flow]
 @variables k(t) [connect = Stream]
+hasconnect(i)
+```
+
+```@example connect
+getconnect(k)
 ```
 
 ## Input or output
@@ -177,8 +182,46 @@ A variable can be marked `irreducible` to prevent it from being moved to an
 `observed` state. This forces the variable to be computed during solving so that
 it can be accessed in [callbacks](@ref events)
 
-```julia
-@variable important_value [irreducible = true]
+```@example metadata
+@variables important_value [irreducible = true]
+isirreducible(important_value)
+```
+
+## 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.
+
+```@example metadata
+@variables important_dof [state_priority = 10] unimportant_dof [state_priority = -2]
+state_priority(important_dof)
+```
+
+## 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"]
+hasunit(speed)
+```
+
+```@example metadata
+getunit(speed)
+```
+
+## Miscellaneous metadata
+
+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]]
+hasmisc(u)
+```
+
+```@example metadata
+getmisc(y)
 ```
 
 ## Additional functions

docs/src/examples/remake.md

@@ -45,12 +45,22 @@ parameters to optimize.
 
 ```@example Remake
 using SymbolicIndexingInterface: parameter_values, state_values
-using SciMLStructures: Tunable, replace, replace!
+using SciMLStructures: Tunable, canonicalize, replace, replace!
+using PreallocationTools
 
 function loss(x, p)
     odeprob = p[1] # ODEProblem stored as parameters to avoid using global variables
     ps = parameter_values(odeprob) # obtain the parameter object from the problem
-    ps = replace(Tunable(), ps, x) # create a copy with the values passed to the loss function
+    diffcache = p[5]
+    # get an appropriately typed preallocated buffer to store the `x` values in
+    buffer = get_tmp(diffcache, x)
+    # copy the current values to this buffer
+    copyto!(buffer, canonicalize(Tunable(), ps)[1])
+    # create a copy of the parameter object with the buffer
+    ps = replace(Tunable(), ps, buffer)
+    # set the updated values in the parameter object
+    setter = p[4]
+    setter(ps, x)
     # remake the problem, passing in our new parameter object
     newprob = remake(odeprob; p = ps)
     timesteps = p[2]
@@ -81,49 +91,22 @@ We can perform the optimization as below:
 ```@example Remake
 using Optimization
 using OptimizationOptimJL
+using SymbolicIndexingInterface
 
 # manually create an OptimizationFunction to ensure usage of `ForwardDiff`, which will
 # require changing the types of parameters from `Float64` to `ForwardDiff.Dual`
 optfn = OptimizationFunction(loss, Optimization.AutoForwardDiff())
+# function to set the parameters we are optimizing
+setter = setp(odeprob, [α, β, γ, δ])
+# `DiffCache` to avoid allocations
+diffcache = DiffCache(canonicalize(Tunable(), parameter_values(odeprob))[1])
 # parameter object is a tuple, to store differently typed objects together
 optprob = OptimizationProblem(
-    optfn, rand(4), (odeprob, timesteps, data), lb = 0.1zeros(4), ub = 3ones(4))
+    optfn, rand(4), (odeprob, timesteps, data, setter, diffcache),
+    lb = 0.1zeros(4), ub = 3ones(4))
 sol = solve(optprob, BFGS())
 ```
 
-To identify which values correspond to which parameters, we can `replace!` them into the
-`ODEProblem`:
-
-```@example Remake
-replace!(Tunable(), parameter_values(odeprob), sol.u)
-odeprob.ps[[α, β, γ, δ]]
-```
-
-`replace!` operates in-place, so the values being replaced must be of the same type as those
-stored in the parameter object, or convertible to that type. For demonstration purposes, we
-can construct a loss function that uses `replace!`, and calculate gradients using
-`AutoFiniteDiff` rather than `AutoForwardDiff`.
-
-```@example Remake
-function loss2(x, p)
-    odeprob = p[1] # ODEProblem stored as parameters to avoid using global variables
-    newprob = remake(odeprob) # copy the problem with `remake`
-    # update the parameter values in-place
-    replace!(Tunable(), parameter_values(newprob), x)
-    timesteps = p[2]
-    sol = solve(newprob, AutoTsit5(Rosenbrock23()); saveat = timesteps)
-    truth = p[3]
-    data = Array(sol)
-    return sum((truth .- data) .^ 2) / length(truth)
-end
-
-# use finite-differencing to calculate derivatives
-optfn2 = OptimizationFunction(loss2, Optimization.AutoFiniteDiff())
-optprob2 = OptimizationProblem(
-    optfn2, rand(4), (odeprob, timesteps, data), lb = 0.1zeros(4), ub = 3ones(4))
-sol = solve(optprob2, BFGS())
-```
-
 # Re-creating the problem
 
 There are multiple ways to re-create a problem with new state/parameter values. We will go

docs/src/systems/DiscreteSystem.md

@@ -0,0 +1,28 @@
+# 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
+  - `discrete_events(sys)`: The set of discrete events in the discrete system.
+
+## Transformations
+
+```@docs; canonical=false
+structural_simplify
+```
+
+## Problem Constructors
+
+```@docs; canonical=false
+DiscreteProblem(sys::DiscreteSystem, u0map, tspan)
+DiscreteFunction(sys::DiscreteSystem, args...)
+```