Merge branch 'master' into myb/ps

Author: YingboMa

CONTRIBUTING.md

@@ -1,3 +1,3 @@
   - This repository follows the [SciMLStyle](https://github.com/SciML/SciMLStyle) and the SciML [ColPrac](https://github.com/SciML/ColPrac).
-  - Please run `using JuliaFormatter, ModelingToolkit; format(joinpath(dirname(pathof(ModelingToolkit)), ".."))` before commiting.
+  - Please run `using JuliaFormatter, ModelingToolkit; format(joinpath(dirname(pathof(ModelingToolkit)), ".."))` before committing.
   - Add tests for any new features.

Project.toml

@@ -1,7 +1,7 @@
 name = "ModelingToolkit"
 uuid = "961ee093-0014-501f-94e3-6117800e7a78"
 authors = ["Yingbo Ma <[email protected]>", "Chris Rackauckas <[email protected]> and contributors"]
-version = "8.64.0"
+version = "8.67.0"
 
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
@@ -31,6 +31,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MLStyle = "d8e11817-5142-5d16-987a-aa16d5891078"
 MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
 NaNMath = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 RuntimeGeneratedFunctions = "7e49a35a-f44a-4d26-94aa-eba1b4ca6b47"
@@ -78,6 +79,7 @@ Latexify = "0.11, 0.12, 0.13, 0.14, 0.15, 0.16"
 MLStyle = "0.4.17"
 MacroTools = "0.5"
 NaNMath = "0.3, 1"
+OrdinaryDiffEq = "6"
 RecursiveArrayTools = "2.3"
 Reexport = "0.2, 1"
 RuntimeGeneratedFunctions = "0.5.9"
@@ -107,7 +109,6 @@ NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
 Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba"
 OptimizationMOI = "fd9f6733-72f4-499f-8506-86b2bdd0dea1"
 OptimizationOptimJL = "36348300-93cb-4f02-beb5-3c3902f8871e"
-OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 ReferenceTests = "324d217c-45ce-50fc-942e-d289b448e8cf"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
@@ -120,4 +121,4 @@ Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["AmplNLWriter", "BenchmarkTools", "ControlSystemsMTK", "NonlinearSolve", "ForwardDiff", "Ipopt", "Ipopt_jll", "ModelingToolkitStandardLibrary", "Optimization", "OptimizationOptimJL", "OptimizationMOI", "OrdinaryDiffEq", "Random", "ReferenceTests", "SafeTestsets", "StableRNGs", "Statistics", "SteadyStateDiffEq", "Test", "StochasticDiffEq", "Sundials", "StochasticDelayDiffEq"]
+test = ["AmplNLWriter", "BenchmarkTools", "ControlSystemsMTK", "NonlinearSolve", "ForwardDiff", "Ipopt", "Ipopt_jll", "ModelingToolkitStandardLibrary", "Optimization", "OptimizationOptimJL", "OptimizationMOI", "Random", "ReferenceTests", "SafeTestsets", "StableRNGs", "Statistics", "SteadyStateDiffEq", "Test", "StochasticDiffEq", "Sundials", "StochasticDelayDiffEq"]

docs/src/basics/FAQ.md

@@ -120,7 +120,7 @@ We can solve this problem by using the `missing_variable_defaults()` function
 prob = ODEProblem(sys, ModelingToolkit.missing_variable_defaults(sys), (0,1))
 ```
 
-This function provides 0 for the default values, which is a safe assumption for dummy derivatives of most models.  However, the 2nd arument allows for a different default value or values to be used if needed.
+This function provides 0 for the default values, which is a safe assumption for dummy derivatives of most models.  However, the 2nd argument allows for a different default value or values to be used if needed.
 
 ```
 julia> ModelingToolkit.missing_variable_defaults(sys, [1,2,3])

src/ModelingToolkit.jl

@@ -106,7 +106,7 @@ abstract type AbstractOptimizationSystem <: AbstractTimeIndependentSystem end
 
 function independent_variable end
 
-# this has to be included early to deal with depency issues
+# this has to be included early to deal with dependency issues
 include("structural_transformation/bareiss.jl")
 function complete end
 function var_derivative! end

src/bipartite_graph.jl

@@ -618,7 +618,7 @@ Essentially the graph adapter performs two largely orthogonal functions
  1. It pairs an undirected bipartite graph with a matching of the destination vertex.
 
     This matching is used to induce an orientation on the otherwise undirected graph:
-    Matched edges pass from destination to source [source to desination], all other edges
+    Matched edges pass from destination to source [source to destination], all other edges
     pass in the opposite direction.
 
  2. It exposes the graph view obtained by contracting the destination [source] vertices
@@ -632,7 +632,7 @@ graph is acyclic.
 # Hypergraph interpretation
 
 Consider the bipartite graph `B` as the incidence graph of some hypergraph `H`.
-Note that a maching `M` on `B` in the above sense is equivalent to determining
+Note that a matching `M` on `B` in the above sense is equivalent to determining
 an (1,n)-orientation on the hypergraph (i.e. each directed hyperedge has exactly
 one head, but any arbitrary number of tails). In this setting, this is simply
 the graph formed by expanding each directed hyperedge into `n` ordinary edges
@@ -685,7 +685,7 @@ function Base.iterate(c::CMONeighbors{false}, (l, state...))
         r === nothing && return nothing
         # If this is a matched edge, skip it, it's reversed in the induced
         # directed graph. Otherwise, if there is no matching for this destination
-        # edge, also skip it, since it got delted in the contraction.
+        # edge, also skip it, since it got deleted in the contraction.
         vsrc = c.g.matching[r[1]]
         if vsrc === c.v || !isa(vsrc, Int)
             state = (r[2],)

src/inputoutput.jl

@@ -112,13 +112,13 @@ end
     same_or_inner_namespace(u, var)
 
 Determine whether `var` is in the same namespace as `u`, or a namespace internal to the namespace of `u`.
-Example: `sys.u ~ sys.inner.u` will bind `sys.inner.u`, but `sys.u` remains an unbound, external signal. The namepsaced signal `sys.inner.u` lives in a namspace internal to `sys`.
+Example: `sys.u ~ sys.inner.u` will bind `sys.inner.u`, but `sys.u` remains an unbound, external signal. The namespaced signal `sys.inner.u` lives in a namespace internal to `sys`.
 """
 function same_or_inner_namespace(u, var)
     nu = get_namespace(u)
     nv = get_namespace(var)
     nu == nv ||           # namespaces are the same
-        startswith(nv, nu) || # or nv starts with nu, i.e., nv is an inner namepsace to nu
+        startswith(nv, nu) || # or nv starts with nu, i.e., nv is an inner namespace to nu
         occursin('₊', string(getname(var))) &&
             !occursin('₊', string(getname(u))) # or u is top level but var is internal
 end
@@ -127,7 +127,7 @@ function inner_namespace(u, var)
     nu = get_namespace(u)
     nv = get_namespace(var)
     nu == nv && return false
-    startswith(nv, nu) || # or nv starts with nu, i.e., nv is an inner namepsace to nu
+    startswith(nv, nu) || # or nv starts with nu, i.e., nv is an inner namespace to nu
         occursin('₊', string(getname(var))) &&
             !occursin('₊', string(getname(u))) # or u is top level but var is internal
 end
@@ -177,7 +177,7 @@ f_ip  : (xout,x,u,p,t) -> nothing
 
 The return values also include the remaining states and parameters, in the order they appear as arguments to `f`.
 
-If `disturbance_inputs` is an array of variables, the generated dynamics function will preserve any state and dynamics associated with distrubance inputs, but the distrubance inputs themselves will not be included as inputs to the generated function. The use case for this is to generate dynamics for state observers that estimate the influence of unmeasured disturbances, and thus require state variables for the disturbance model, but without disturbance inputs since the disturbances are not available for measurement.
+If `disturbance_inputs` is an array of variables, the generated dynamics function will preserve any state and dynamics associated with disturbance inputs, but the disturbance inputs themselves will not be included as inputs to the generated function. The use case for this is to generate dynamics for state observers that estimate the influence of unmeasured disturbances, and thus require state variables for the disturbance model, but without disturbance inputs since the disturbances are not available for measurement.
 See [`add_input_disturbance`](@ref) for a higher-level interface to this functionality.
 
 # Example

src/structural_transformation/bareiss.jl

@@ -34,7 +34,7 @@ else
     function Base.swapcols!(A::AbstractSparseMatrixCSC, i, j)
         i == j && return
 
-        # For simplicitly, let i denote the smaller of the two columns
+        # For simplicity, let i denote the smaller of the two columns
         j < i && @swap(i, j)
 
         colptr = getcolptr(A)
@@ -72,7 +72,7 @@ else
         return nothing
     end
     function swaprows!(A::AbstractSparseMatrixCSC, i, j)
-        # For simplicitly, let i denote the smaller of the two rows
+        # For simplicity, let i denote the smaller of the two rows
         j < i && @swap(i, j)
 
         rows = rowvals(A)
@@ -184,7 +184,7 @@ const bareiss_virtcolswap = ((M, i, j) -> nothing, swaprows!,
 Perform Bareiss's fraction-free row-reduction algorithm on the matrix `M`.
 Optionally, a specific pivoting method may be specified.
 
-swap_strategy is an optional argument that determines how the swapping of rows and coulmns is performed.
+swap_strategy is an optional argument that determines how the swapping of rows and columns is performed.
 bareiss_colswap (the default) swaps the columns and rows normally.
 bareiss_virtcolswap pretends to swap the columns which can be faster for sparse matrices.
 """

src/structural_transformation/partial_state_selection.jl

@@ -209,7 +209,7 @@ function dummy_derivative_graph!(structure::SystemStructure, var_eq_matching, ja
             J = nothing
         else
             _J = jac(eqs, vars)
-            # only accecpt small intergers to avoid overflow
+            # only accept small integers to avoid overflow
             is_all_small_int = all(_J) do x′
                 x = unwrap(x′)
                 x isa Number || return false

src/structural_transformation/symbolics_tearing.jl

@@ -103,7 +103,7 @@ function full_equations(sys::AbstractSystem; simplify = false)
             if rhs isa Symbolic
                 return 0 ~ rhs
             else # a number
-                error("tearing failled because the system is singular")
+                error("tearing failed because the system is singular")
             end
         end
         eq
@@ -194,7 +194,7 @@ function to_mass_matrix_form(neweqs, ieq, graph, fullvars, isdervar::F,
         return (new_lhs ~ new_rhs), invview(var_to_diff)[dervar]
     else # a number
         if abs(rhs) > 100eps(float(rhs))
-            @warn "The equation $eq is not consistent. It simplifed to 0 == $rhs."
+            @warn "The equation $eq is not consistent. It simplified to 0 == $rhs."
         end
         return nothing
     end

src/systems/abstractsystem.jl

@@ -1,3 +1,4 @@
+using OrdinaryDiffEq
 const SYSTEM_COUNT = Threads.Atomic{UInt}(0)
 
 get_component_type(x::AbstractSystem) = get_gui_metadata(x).type
@@ -471,15 +472,18 @@ function namespace_defaults(sys)
          for (k, v) in pairs(defs))
 end
 
-function namespace_equations(sys::AbstractSystem)
+function namespace_equations(sys::AbstractSystem, ivs = independent_variables(sys))
     eqs = equations(sys)
     isempty(eqs) && return Equation[]
-    map(eq -> namespace_equation(eq, sys), eqs)
+    map(eq -> namespace_equation(eq, sys; ivs), eqs)
 end
 
-function namespace_equation(eq::Equation, sys, n = nameof(sys))
-    _lhs = namespace_expr(eq.lhs, sys, n)
-    _rhs = namespace_expr(eq.rhs, sys, n)
+function namespace_equation(eq::Equation,
+    sys,
+    n = nameof(sys);
+    ivs = independent_variables(sys))
+    _lhs = namespace_expr(eq.lhs, sys, n; ivs)
+    _rhs = namespace_expr(eq.rhs, sys, n; ivs)
     _lhs ~ _rhs
 end
 
@@ -489,26 +493,29 @@ function namespace_assignment(eq::Assignment, sys)
     Assignment(_lhs, _rhs)
 end
 
-function namespace_expr(O, sys, n = nameof(sys))
-    ivs = independent_variables(sys)
+function namespace_expr(O, sys, n = nameof(sys); ivs = independent_variables(sys))
     O = unwrap(O)
     if any(isequal(O), ivs)
         return O
-    elseif isvariable(O)
-        renamespace(n, O)
     elseif istree(O)
         T = typeof(O)
-        if symtype(operation(O)) <: FnType
-            renamespace(n, O)::T
+        renamed = let sys = sys, n = n, T = T
+            map(a -> namespace_expr(a, sys, n; ivs)::Any, arguments(O))
+        end
+        if isvariable(O)
+            # Use renamespace so the scope is correct, and make sure to use the
+            # metadata from the rescoped variable
+            rescoped = renamespace(n, O)
+            similarterm(O, operation(rescoped), renamed,
+                metadata = metadata(rescoped))::T
         else
-            renamed = let sys = sys, n = n, T = T
-                map(a -> namespace_expr(a, sys, n)::Any, arguments(O))
-            end
-            similarterm(O, operation(O), renamed)::T
+            similarterm(O, operation(O), renamed, metadata = metadata(O))::T
         end
+    elseif isvariable(O)
+        renamespace(n, O)
     elseif O isa Array
         let sys = sys, n = n
-            map(o -> namespace_expr(o, sys, n), O)
+            map(o -> namespace_expr(o, sys, n; ivs), O)
         end
     else
         O
@@ -656,20 +663,6 @@ function SymbolicIndexingInterface.is_param_sym(sys::AbstractSystem, sym)
     !isnothing(SymbolicIndexingInterface.param_sym_to_index(sys, sym))
 end
 
-struct AbstractSysToExpr
-    sys::AbstractSystem
-    states::Vector
-end
-AbstractSysToExpr(sys) = AbstractSysToExpr(sys, states(sys))
-function (f::AbstractSysToExpr)(O)
-    !istree(O) && return toexpr(O)
-    any(isequal(O), f.states) && return nameof(operation(O))  # variables
-    if issym(operation(O))
-        return build_expr(:call, Any[nameof(operation(O)); f.(arguments(O))])
-    end
-    return build_expr(:call, Any[operation(O); f.(arguments(O))])
-end
-
 ###
 ### System utils
 ###
@@ -1249,7 +1242,7 @@ function io_preprocessing(sys::AbstractSystem, inputs,
 end
 
 """
-    lin_fun, simplified_sys = linearization_function(sys::AbstractSystem, inputs, outputs; simplify = false, kwargs...)
+    lin_fun, simplified_sys = linearization_function(sys::AbstractSystem, inputs, outputs; simplify = false, initialize = true, kwargs...)
 
 Return a function that linearizes the system `sys`. The function [`linearize`](@ref) provides a higher-level and easier to use interface.
 
@@ -1273,12 +1266,14 @@ The `simplified_sys` has undergone [`structural_simplify`](@ref) and had any occ
   - `inputs`: A vector of variables that indicate the inputs of the linearized input-output model.
   - `outputs`: A vector of variables that indicate the outputs of the linearized input-output model.
   - `simplify`: Apply simplification in tearing.
+  - `initialize`: If true, a check is performed to ensure that the operating point is consistent (satisfies algebraic equations). If the op is not consistent, initialization is performed.
   - `kwargs`: Are passed on to `find_solvables!`
 
 See also [`linearize`](@ref) which provides a higher-level interface.
 """
 function linearization_function(sys::AbstractSystem, inputs,
     outputs; simplify = false,
+    initialize = true,
     kwargs...)
     sys, diff_idxs, alge_idxs, input_idxs = io_preprocessing(sys, inputs, outputs; simplify,
         kwargs...)
@@ -1294,6 +1289,14 @@ function linearization_function(sys::AbstractSystem, inputs,
             if u !== nothing # Handle systems without states
                 length(sts) == length(u) ||
                     error("Number of state variables ($(length(sts))) does not match the number of input states ($(length(u)))")
+                if initialize && !isempty(alge_idxs) # This is expensive and can be omitted if the user knows that the system is already initialized
+                    residual = fun(u, p, t)
+                    if norm(residual[alge_idxs]) > √(eps(eltype(residual)))
+                        prob = ODEProblem(fun, u, (t, t + 1), p)
+                        integ = init(prob, Rodas5P())
+                        u = integ.u
+                    end
+                end
                 uf = SciMLBase.UJacobianWrapper(fun, t, p)
                 fg_xz = ForwardDiff.jacobian(uf, u)
                 h_xz = ForwardDiff.jacobian(let p = p, t = t
@@ -1397,7 +1400,7 @@ function linearize_symbolic(sys::AbstractSystem, inputs,
             if !allow_input_derivatives
                 der_inds = findall(vec(any(!iszero, Bs, dims = 1)))
                 @show typeof(der_inds)
-                error("Input derivatives appeared in expressions (-g_z\\g_u != 0), the following inputs appeared differentiated: $(ModelingToolkit.inputs(sys)[der_inds]). Call `linear_statespace` with keyword argument `allow_input_derivatives = true` to allow this and have the returned `B` matrix be of double width ($(2nu)), where the last $nu inputs are the derivatives of the first $nu inputs.")
+                error("Input derivatives appeared in expressions (-g_z\\g_u != 0), the following inputs appeared differentiated: $(ModelingToolkit.inputs(sys)[der_inds]). Call `linearize_symbolic` with keyword argument `allow_input_derivatives = true` to allow this and have the returned `B` matrix be of double width ($(2nu)), where the last $nu inputs are the derivatives of the first $nu inputs.")
             end
             B = [B [zeros(nx, nu); Bs]]
             D = [D zeros(ny, nu)]
@@ -1413,7 +1416,12 @@ function markio!(state, orig_inputs, inputs, outputs; check = true)
     outputset = Dict{Any, Bool}(o => false for o in outputs)
     for (i, v) in enumerate(fullvars)
         if v in keys(inputset)
-            v = setio(v, true, false)
+            if v in keys(outputset)
+                v = setio(v, true, true)
+                outputset[v] = true
+            else
+                v = setio(v, true, false)
+            end
             inputset[v] = true
             fullvars[i] = v
         elseif v in keys(outputset)
@@ -1442,8 +1450,8 @@ function markio!(state, orig_inputs, inputs, outputs; check = true)
 end
 
 """
-    (; A, B, C, D), simplified_sys = linearize(sys, inputs, outputs;    t=0.0, op = Dict(), allow_input_derivatives = false, kwargs...)
-    (; A, B, C, D)                 = linearize(simplified_sys, lin_fun; t=0.0, op = Dict(), allow_input_derivatives = false)
+    (; A, B, C, D), simplified_sys = linearize(sys, inputs, outputs;    t=0.0, op = Dict(), allow_input_derivatives = false, zero_dummy_der=false, kwargs...)
+    (; A, B, C, D)                 = linearize(simplified_sys, lin_fun; t=0.0, op = Dict(), allow_input_derivatives = false, zero_dummy_der=false)
 
 Return a NamedTuple with the matrices of a linear statespace representation
 on the form
@@ -1463,6 +1471,8 @@ the default values of `sys` are used.
 
 If `allow_input_derivatives = false`, an error will be thrown if input derivatives (``u̇``) appear as inputs in the linearized equations. If input derivatives are allowed, the returned `B` matrix will be of double width, corresponding to the input `[u; u̇]`.
 
+`zero_dummy_der` can be set to automatically set the operating point to zero for all dummy derivatives.
+
 See also [`linearization_function`](@ref) which provides a lower-level interface, [`linearize_symbolic`](@ref) and [`ModelingToolkit.reorder_states`](@ref).
 
 See extended help for an example.
@@ -1576,7 +1586,7 @@ function linearize(sys, lin_fun; t = 0.0, op = Dict(), allow_input_derivatives =
         if !iszero(Bs)
             if !allow_input_derivatives
                 der_inds = findall(vec(any(!=(0), Bs, dims = 1)))
-                error("Input derivatives appeared in expressions (-g_z\\g_u != 0), the following inputs appeared differentiated: $(inputs(sys)[der_inds]). Call `linear_statespace` with keyword argument `allow_input_derivatives = true` to allow this and have the returned `B` matrix be of double width ($(2nu)), where the last $nu inputs are the derivatives of the first $nu inputs.")
+                error("Input derivatives appeared in expressions (-g_z\\g_u != 0), the following inputs appeared differentiated: $(inputs(sys)[der_inds]). Call `linearize` with keyword argument `allow_input_derivatives = true` to allow this and have the returned `B` matrix be of double width ($(2nu)), where the last $nu inputs are the derivatives of the first $nu inputs.")
             end
             B = [B [zeros(nx, nu); Bs]]
             D = [D zeros(ny, nu)]