Merge branch 'SciML:master' into iss3707

Author: fchen121

.github/workflows/Downstream.yml

@@ -6,6 +6,7 @@ on:
   pull_request:
     paths-ignore:
         - 'docs/**'
+        - 'benchmark/**'
 
 concurrency:
   # Skip intermediate builds: always, but for the master branch and tags.
@@ -26,6 +27,7 @@ jobs:
         os: [ubuntu-latest]
         package:
           - {user: SciML, repo: SciMLBase.jl, group: Downstream}
+          - {user: SciML, repo: SciMLBase.jl, group: SymbolicIndexingInterface}
           - {user: SciML, repo: Catalyst.jl, group: All}
           - {user: SciML, repo: CellMLToolkit.jl, group: Core}
           - {user: SciML, repo: SBMLToolkit.jl, group: All}

.github/workflows/ReleaseTest.yml

@@ -6,6 +6,7 @@ on:
   pull_request:
     paths-ignore:
         - 'docs/**'
+        - 'benchmark/**'
 
 concurrency:
   # Skip intermediate builds: always, but for the master branch and tags.

.github/workflows/Tests.yml

@@ -13,6 +13,7 @@ on:
       - master
     paths-ignore:
       - 'docs/**'
+      - 'benchmark/**'
 
 concurrency:
   # Skip intermediate builds: always, but for the master branch.

.github/workflows/benchmark.yml

@@ -22,3 +22,5 @@ jobs:
       - uses: MilesCranmer/AirspeedVelocity.jl@action-v1
         with:
           julia-version: "${{ matrix.version }}"
+          script: "benchmark/benchmarks.jl"
+          extra-pkgs: "ModelingToolkitStandardLibrary,OrdinaryDiffEqDefault"

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 = "10.2.0"
+version = "10.5.0"
 
 [deps]
 ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
@@ -98,7 +98,7 @@ Combinatorics = "1"
 CommonSolve = "0.2.4"
 Compat = "3.42, 4"
 ConstructionBase = "1"
-DataInterpolations = "6.4"
+DataInterpolations = "7, 8"
 DataStructures = "0.17, 0.18"
 DeepDiffs = "1"
 DelayDiffEq = "5.50"
@@ -129,6 +129,7 @@ LabelledArrays = "1.3"
 Latexify = "0.11, 0.12, 0.13, 0.14, 0.15, 0.16"
 Libdl = "1"
 LinearAlgebra = "1"
+LinearSolve = "3"
 Logging = "1"
 MLStyle = "0.4.17"
 ModelingToolkitStandardLibrary = "2.20"
@@ -148,7 +149,7 @@ RecursiveArrayTools = "3.26"
 Reexport = "0.2, 1"
 RuntimeGeneratedFunctions = "0.5.9"
 SCCNonlinearSolve = "1.0.0"
-SciMLBase = "2.91.1"
+SciMLBase = "2.100.0"
 SciMLPublic = "1.0.0"
 SciMLStructures = "1.7"
 Serialization = "1"
@@ -180,6 +181,7 @@ ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9"
 Ipopt_jll = "9cc047cb-c261-5740-88fc-0cf96f7bdcc7"
 JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b"
+LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 ModelingToolkitStandardLibrary = "16a59e39-deab-5bd0-87e4-056b12336739"
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
@@ -205,4 +207,4 @@ Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["AmplNLWriter", "BenchmarkTools", "BoundaryValueDiffEqMIRK", "BoundaryValueDiffEqAscher", "ControlSystemsBase", "DataInterpolations", "DelayDiffEq", "NonlinearSolve", "ForwardDiff", "Ipopt", "Ipopt_jll", "ModelingToolkitStandardLibrary", "Optimization", "OptimizationOptimJL", "OptimizationMOI", "OrdinaryDiffEq", "OrdinaryDiffEqCore", "OrdinaryDiffEqDefault", "REPL", "Random", "ReferenceTests", "SafeTestsets", "StableRNGs", "Statistics", "SteadyStateDiffEq", "Test", "StochasticDiffEq", "Sundials", "StochasticDelayDiffEq", "Pkg", "JET", "OrdinaryDiffEqNonlinearSolve", "Logging", "OptimizationBase"]
+test = ["AmplNLWriter", "BenchmarkTools", "BoundaryValueDiffEqMIRK", "BoundaryValueDiffEqAscher", "ControlSystemsBase", "DataInterpolations", "DelayDiffEq", "NonlinearSolve", "ForwardDiff", "Ipopt", "Ipopt_jll", "ModelingToolkitStandardLibrary", "Optimization", "OptimizationOptimJL", "OptimizationMOI", "OrdinaryDiffEq", "OrdinaryDiffEqCore", "OrdinaryDiffEqDefault", "REPL", "Random", "ReferenceTests", "SafeTestsets", "StableRNGs", "Statistics", "SteadyStateDiffEq", "Test", "StochasticDiffEq", "Sundials", "StochasticDelayDiffEq", "Pkg", "JET", "OrdinaryDiffEqNonlinearSolve", "Logging", "OptimizationBase", "LinearSolve"]

benchmark/Project.toml

@@ -1,6 +1,8 @@
 using ModelingToolkit, BenchmarkTools
 using ModelingToolkitStandardLibrary
-using ModelingToolkitStandardLibrary.Thermal
+using ModelingToolkitStandardLibrary.Electrical
+using ModelingToolkitStandardLibrary.Mechanical.Rotational
+using ModelingToolkitStandardLibrary.Blocks
 using OrdinaryDiffEqDefault
 
 const SUITE = BenchmarkGroup()
@@ -51,4 +53,4 @@ tspan = (0.0, 6.0)
 SUITE["ODEProblem"] = @benchmarkable ODEProblem($model, $u0, $tspan)
 
 prob = ODEProblem(model, u0, tspan)
-SUITE["init"] = init($prob)
+SUITE["init"] = @benchmarkable init($prob)

benchmark/benchmarks.jl

@@ -35,7 +35,7 @@ Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 Attractors = "1.24"
 BenchmarkTools = "1.3"
 BifurcationKit = "0.4"
-CairoMakie = "0.13"
+CairoMakie = "0.13, 0.15"
 CommonSolve = "0.2"
 DataInterpolations = "6.5, 8"
 DiffEqDevTools = "2"

docs/Project.toml

@@ -21,6 +21,8 @@ ModelingToolkit.generate_constraint_hessian
 ModelingToolkit.generate_control_jacobian
 ModelingToolkit.build_explicit_observed_function
 ModelingToolkit.generate_control_function
+ModelingToolkit.generate_update_A
+ModelingToolkit.generate_update_b
 ```
 
 For functions such as jacobian calculation which require symbolic computation, there
@@ -42,6 +44,7 @@ ModelingToolkit.cost_hessian_sparsity
 ModelingToolkit.calculate_constraint_jacobian
 ModelingToolkit.calculate_constraint_hessian
 ModelingToolkit.calculate_control_jacobian
+ModelingToolkit.calculate_A_b
 ```
 
 All code generation eventually calls `build_function_wrapper`.

docs/src/API/codegen.md

@@ -1,4 +1,4 @@
-### Solvers
+# [Dynamic Optimization Solvers](@id dynamic_opt_api)
 
 Currently 4 backends are exposed for solving dynamic optimization problems using collocation: JuMP, InfiniteOpt, CasADi, and Pyomo.
 

docs/src/API/dynamic_opt.md

@@ -29,7 +29,7 @@ SciMLBase.DiscreteProblem
 SciMLBase.ImplicitDiscreteProblem
 ```
 
-## Nonlinear systems
+## Linear and Nonlinear systems
 
 ```@docs
 SciMLBase.NonlinearFunction
@@ -41,6 +41,7 @@ SciMLBase.IntervalNonlinearFunction
 SciMLBase.IntervalNonlinearProblem
 ModelingToolkit.HomotopyContinuationProblem
 SciMLBase.HomotopyNonlinearFunction
+SciMLBase.LinearProblem
 ```
 
 ## Optimization and optimal control

docs/src/API/problems.md

@@ -27,7 +27,7 @@ are similarly undocumented. Following is the list of behaviors that should be re
     parameter with the given index.
   - `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
+  - `parameter_index(sys, sym)` will return a `ParameterIndex` object if `sys` has been
     `complete`d (through `mtkcompile`, `complete` or `@mtkcompile`).
   - `copy(::MTKParameters)` is defined and duplicates the parameter object, including the
     memory used by the underlying buffers.

docs/src/basics/FAQ.md

@@ -93,12 +93,13 @@ eqs_pert = taylor_coeff(eq_pert, ϵ, 0:2)
 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:
 
 ```@example perturbation
-u0 = Dict([unknowns(sys) .=> 0.0; D(y[0]) => 1.0]) # nonzero initial velocity
+u0 = [y[0] => 0.0, y[1] => 0.0, y[2] => 0.0, D(y[0]) => 1.0, D(y[1]) => 0.0, D(y[2]) => 0.0] # nonzero initial velocity
 prob = ODEProblem(sys, u0, (0.0, 50.0))
 sol = solve(prob)
 plot(sol, idxs = substitute(x_series, ϵ => 0.1); label = "Perturbative (ϵ=0.1)")
 
 x_exact(t, ϵ) = exp(-ϵ * t) * sin(√(1 - ϵ^2) * t) / √(1 - ϵ^2)
+@assert isapprox(sol(π/2; idxs = substitute(x_series, ϵ => 0.1)), x_exact(π/2, 0.1); atol = 1e-2) # compare around 1st peak # hide
 plot!(sol.t, x_exact.(sol.t, 0.1); label = "Exact (ϵ=0.1)")
 ```
 

docs/src/examples/perturbation.md

@@ -161,6 +161,7 @@ include("systems/index_cache.jl")
 include("systems/parameter_buffer.jl")
 include("systems/abstractsystem.jl")
 include("systems/model_parsing.jl")
+include("systems/connectiongraph.jl")
 include("systems/connectors.jl")
 include("systems/state_machines.jl")
 include("systems/analysis_points.jl")
@@ -188,6 +189,7 @@ include("problems/jumpproblem.jl")
 include("problems/initializationproblem.jl")
 include("problems/sccnonlinearproblem.jl")
 include("problems/bvproblem.jl")
+include("problems/linearproblem.jl")
 
 include("modelingtoolkitize/common.jl")
 include("modelingtoolkitize/odeproblem.jl")

src/ModelingToolkit.jl

@@ -175,5 +175,5 @@ macro brownian(xs...)
     return quote
         Base.depwarn("`@brownian` is deprecated. Use `@brownians` instead", :brownian_macro)
         $(@__MODULE__).@brownians $(xs...)
-    end
+    end |> esc
 end

src/deprecations.jl

@@ -1,5 +1,15 @@
 using Symbolics: Operator, Num, Term, value, recursive_hasoperator
 
+"""
+    $(TYPEDSIGNATURES)
+
+Trait to be implemented for operators which determines whether application of the operator
+generates a semantically different variable or not. For example, `Differential` and `Shift`
+are not transparent but `Sample` and `Hold` are. Defaults to `false` if not implemented.
+"""
+is_transparent_operator(x) = is_transparent_operator(typeof(x))
+is_transparent_operator(::Type) = false
+
 """
     function SampleTime()
 
@@ -128,6 +138,8 @@ struct Sample <: Operator
     Sample(clock::Union{TimeDomain, InferredTimeDomain} = InferredDiscrete()) = new(clock)
 end
 
+is_transparent_operator(::Type{Sample}) = true
+
 function Sample(arg::Real)
     arg = unwrap(arg)
     if symbolic_type(arg) == NotSymbolic()
@@ -179,6 +191,9 @@ cont_x = Hold()(disc_x)
 """
 struct Hold <: Operator
 end
+
+is_transparent_operator(::Type{Hold}) = true
+
 (D::Hold)(x) = Term{symtype(x)}(D, Any[x])
 (D::Hold)(x::Num) = Num(D(value(x)))
 SymbolicUtils.promote_symtype(::Hold, x) = x

src/discretedomain.jl

@@ -45,6 +45,7 @@ function linearization_function(sys::AbstractSystem, inputs,
         warn_initialize_determined = true,
         guesses = Dict(),
         warn_empty_op = true,
+        t = 0.0,
         kwargs...)
     op = Dict(op)
     if isempty(op) && warn_empty_op
@@ -73,7 +74,7 @@ function linearization_function(sys::AbstractSystem, inputs,
     end
 
     prob = ODEProblem{true, SciMLBase.FullSpecialize}(
-        sys, merge(op, anydict(p)), (nothing, nothing); allow_incomplete = true,
+        sys, merge(op, anydict(p)), (t, t); allow_incomplete = true,
         algebraic_only = true, guesses)
     u0 = state_values(prob)
 
@@ -753,7 +754,7 @@ function linearize(sys, inputs, outputs; op = Dict(), t = 0.0,
         inputs,
         outputs;
         zero_dummy_der,
-        op,
+        op, t,
         kwargs...)
     mats, extras = linearize(ssys, lin_fun; op, t, allow_input_derivatives)
     mats, ssys, extras

src/linearization.jl

@@ -46,7 +46,7 @@ end
 """
     $(TYPEDSIGNATURES)
 
-Return a symbolic state for the given proble `prob.`. `t` is the independent variable.
+Return a symbolic state for the given problem `prob.`. `t` is the independent variable.
 `u_names` optionally contains the names to use for the created symbolic variables.
 """
 function construct_vars(prob, t, u_names = nothing)
@@ -287,7 +287,7 @@ end
 """
     $(TYPEDSIGNATURES)
 
-Return a symbolic parameter object for the given proble `prob.`. `t` is the independent
+Return a symbolic parameter object for the given problem `prob.`. `t` is the independent
 variable. `p_names` optionally contains the names to use for the created symbolic
 variables.
 """
@@ -319,7 +319,7 @@ end
     $(TYPEDSIGNATURES)
 
 Given the differential operator `D`, mass matrix `mm` and ordered list of unknowns `vars`,
-return the list of 
+return the list of
 """
 function lhs_from_mass_matrix(D, mm, vars)
     var_set = Set(vars)

src/modelingtoolkitize/common.jl

@@ -169,3 +169,12 @@ function check_no_equations(sys::System, T)
         """))
     end
 end
+
+function check_affine(sys::System, T)
+    if !isaffine(sys)
+        throw(SystemCompatibilityError("""
+        A non-affine system cannot be used to construct a `$T`. Consider a
+        `NonlinearProblem` instead.
+        """))
+    end
+end

src/problems/compatibility.jl

@@ -72,7 +72,7 @@ end
         eval_module, check_compatibility, implicit_dae = true, expression, kwargs...)
 
     kwargs = process_kwargs(sys; expression, callback, eval_expression, eval_module,
-        kwargs...)
+        op, kwargs...)
 
     diffvars = collect_differential_variables(sys)
     sts = unknowns(sys)

src/problems/daeproblem.jl

@@ -66,7 +66,7 @@ end
     end
 
     kwargs = process_kwargs(
-        sys; expression, callback, eval_expression, eval_module, kwargs...)
+        sys; expression, callback, eval_expression, eval_module, op, kwargs...)
     args = (; f, u0, h, tspan, p)
 
     return maybe_codegen_scimlproblem(expression, DDEProblem{iip}, args; kwargs...)

src/problems/ddeproblem.jl

@@ -391,3 +391,32 @@ $PROBLEM_INTERNALS_HEADER
 
 $PROBLEM_INTERNAL_KWARGS
 """ SciMLBase.IntervalNonlinearProblem
+
+@doc """
+    SciMLBase.LinearProblem(sys::System, op; kwargs...)
+    SciMLBase.LinearProblem{iip}(sys::System, op; kwargs...)
+
+Build a `LinearProblem` given a system `sys` and operating point `op`. `iip` is a boolean
+indicating whether the problem should be in-place. 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 `LinearProblem`. Any values not provided will
+fallback to the corresponding default (if present).
+
+Note that since `u0` is optional for `LinearProblem`, values of unknowns do not need to be
+specified in `op` to create a `LinearProblem`. In such a case, `prob.u0` will be `nothing`
+and attempting to symbolically index the problem with an unknown, observable, or expression
+depending on unknowns/observables will error.
+
+Updating the parameters automatically updates the `A` and `b` arrays.
+
+# Keyword arguments
+
+$PROBLEM_KWARGS
+$(prob_fun_common_kwargs(LinearProblem, false))
+
+All other keyword arguments are forwarded to the $func constructor.
+
+$PROBLEM_INTERNALS_HEADER
+
+$PROBLEM_INTERNAL_KWARGS
+""" SciMLBase.LinearProblem

src/problems/docs.jl

@@ -80,7 +80,8 @@
     end
 
     # handle events, making sure to reset aggregators in the generated affect functions
-    cbs = process_events(sys; callback, eval_expression, eval_module, reset_jumps = true)
+    cbs = process_events(
+        sys; callback, eval_expression, eval_module, op, reset_jumps = true)
 
     if rng !== nothing
         kwargs = (; kwargs..., rng)