Merge branch 'master' into myb/obsfix

Author: YingboMa

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.23.0"
+version = "8.25.0"
 
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
@@ -36,6 +36,7 @@ RuntimeGeneratedFunctions = "7e49a35a-f44a-4d26-94aa-eba1b4ca6b47"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 Serialization = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
 Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46"
+SimpleWeightedGraphs = "47aef6b3-ad0c-573a-a1e2-d07658019622"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
@@ -72,6 +73,7 @@ Reexport = "0.2, 1"
 RuntimeGeneratedFunctions = "0.4.3, 0.5"
 SciMLBase = "1.58.0"
 Setfield = "0.7, 0.8, 1"
+SimpleWeightedGraphs = "1"
 SpecialFunctions = "0.7, 0.8, 0.9, 0.10, 1.0, 2"
 StaticArrays = "0.10, 0.11, 0.12, 1.0"
 SymbolicUtils = "0.19"
@@ -93,9 +95,9 @@ 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"
 StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
-SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 SteadyStateDiffEq = "9672c7b4-1e72-59bd-8a11-6ac3964bc41f"
 StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0"
 Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4"

docs/pages.jl

@@ -17,6 +17,7 @@ pages = [
                     "basics/Variable_metadata.md",
                     "basics/Composition.md",
                     "basics/Events.md",
+                    "basics/Linearization.md",
                     "basics/Validation.md",
                     "basics/DependencyGraphs.md",
                     "basics/FAQ.md"],

docs/src/basics/Linearization.md

@@ -0,0 +1,64 @@
+# [Linearization](@id linearization)
+A nonlinear dynamical system with state (differential and algebraic) ``x`` and input signals ``u``
+```math
+M \dot x = f(x, u)
+```
+can be linearized using the function [`linearize`](@ref) to produce a linear statespace system on the form
+```math
+\begin{aligned}
+\dot x &= Ax + Bu\\
+y &= Cx + Du
+\end{aligned}
+```
+
+The `linearize` function expects the user to specify the inputs ``u`` and the outputs ``u`` using the syntax shown in the example below:
+
+## Example
+```@example LINEARIZE
+using ModelingToolkit
+@variables t x(t)=0 y(t)=0 u(t)=0 r(t)=0
+@parameters kp = 1
+D = Differential(t)
+
+eqs = [u ~ kp * (r - y) # P controller
+       D(x) ~ -x + u    # First-order plant
+       y ~ x]           # Output equation
+
+@named sys = ODESystem(eqs, t)
+matrices, simplified_sys = linearize(sys, [r], [y]) # Linearize from r to y
+matrices
+```
+The named tuple `matrices` contains the matrices of the linear statespace representation, while `simplified_sys` is an `ODESystem` that, amongst other things, indicates the state order in the linear system through
+```@example LINEARIZE
+using ModelingToolkit: inputs, outputs
+[states(simplified_sys); inputs(simplified_sys); outputs(simplified_sys)]
+```
+
+## Operating point
+The operating point to linearize around can be specified with the keyword argument `op` like this: `op = Dict(x => 1, r => 2)`. 
+
+## Batch linearization and algebraic variables
+If linearization is to be performed around multiple operating points, the simplification of the system has to be carried out a single time only. To facilitate this, the lower-level function [`ModelingToolkit.linearization_function`](@ref) is available. This function further allows you to obtain separate Jacobians for the differential and algebraic parts of the model. For ODE models without algebraic equations, the statespace representation above is available from the output of `linearization_function` as `A, B, C, D = f_x, f_u, h_x, h_u`.
+
+
+## Input derivatives
+Physical systems are always *proper*, i.e., they do not differentiate causal inputs. However, ModelingToolkit allows you to model non-proper systems, such as inverse models, and may sometimes fail to find a realization of a proper system on proper form. In these situations, `linearize` may throw an error mentioning
+```
+Input derivatives appeared in expressions (-g_z\g_u != 0)
+```
+This means that to simulate this system, some order of derivatives of the input is required. To allow `linearize` to proceed in this situation, one may pass the keyword argument `allow_input_derivatives = true`, in which case the resulting model will have twice as many inputs, ``2n_u``, where the last ``n_u`` inputs correspond to ``\dot u``. 
+
+If the modeled system is actually proper (but MTK failed to find a proper realization), further numerical simplification can be applied to the resulting statespace system to obtain a proper form. Such simplification is currently available in the experimental package [ControlSystemsMTK](https://github.com/baggepinnen/ControlSystemsMTK.jl#internals-transformation-of-non-proper-models-to-proper-statespace-form).
+
+
+## Tools for linear analysis
+[ModelingToolkitStandardLibrary](http://mtkstdlib.sciml.ai/dev/API/linear_analysis/) contains a set of [tools for more advanced linear analysis](http://mtkstdlib.sciml.ai/dev/API/linear_analysis/). These can be used to make it easier to work with and analyze causal models, such as control and signal-processing systems. 
+
+```@index
+Pages = ["Linearization.md"]
+```
+
+```@docs
+linearize
+ModelingToolkit.linearization_function
+```

src/ModelingToolkit.jl

@@ -147,6 +147,7 @@ include("systems/dependency_graphs.jl")
 include("systems/systemstructure.jl")
 using .SystemStructures
 
+include("debugging.jl")
 include("systems/alias_elimination.jl")
 include("structural_transformation/StructuralTransformations.jl")
 
@@ -176,7 +177,8 @@ export NonlinearSystem, OptimizationSystem
 export alias_elimination, flatten
 export connect, @connector, Connection, Flow, Stream, instream
 export isinput, isoutput, getbounds, hasbounds, isdisturbance, istunable, getdist, hasdist,
-       tunable_parameters, isirreducible, getdescription, hasdescription
+       tunable_parameters, isirreducible, getdescription, hasdescription, isbinaryvar,
+       isintegervar
 export ode_order_lowering, dae_order_lowering, liouville_transform
 export PDESystem
 export Differential, expand_derivatives, @derivatives
@@ -185,7 +187,7 @@ export Term, Sym
 export SymScope, LocalScope, ParentScope, GlobalScope
 export independent_variables, independent_variable, states, parameters, equations, controls,
        observed, structure, full_equations
-export structural_simplify, expand_connections, linearize, linear_statespace
+export structural_simplify, expand_connections, linearize, linearization_function
 export DiscreteSystem, DiscreteProblem
 
 export calculate_jacobian, generate_jacobian, generate_function
@@ -209,5 +211,6 @@ export build_function
 export modelingtoolkitize
 export @variables, @parameters
 export @named, @nonamespace, @namespace, extend, compose
+export debug_system
 
 end # module

src/bipartite_graph.jl

@@ -63,7 +63,17 @@ end
 function Base.setindex!(m::Matching{U}, v::Union{Integer, U}, i::Integer) where {U}
     if m.inv_match !== nothing
         oldv = m.match[i]
-        isa(oldv, Int) && (m.inv_match[oldv] = unassigned)
+        # TODO: maybe default Matching to always have an `inv_match`?
+
+        # To maintain the invariant that `m.inv_match[m.match[i]] == i`, we need
+        # to unassign the matching at `m.inv_match[v]` if it exists.
+        if v isa Int && (iv = m.inv_match[v]) isa Int
+            m.match[iv] = unassigned
+        end
+        if isa(oldv, Int)
+            @assert m.inv_match[oldv] == i
+            m.inv_match[oldv] = unassigned
+        end
         isa(v, Int) && (m.inv_match[v] = i)
     end
     return m.match[i] = v
@@ -191,14 +201,49 @@ end
 # Matrix whose only purpose is to pretty-print the bipartite graph
 struct BipartiteAdjacencyList
     u::Union{Vector{Int}, Nothing}
+    highligh_u::Union{Set{Int}, Nothing}
+    match::Union{Int, Unassigned}
+end
+function BipartiteAdjacencyList(u::Union{Vector{Int}, Nothing})
+    BipartiteAdjacencyList(u, nothing, unassigned)
+end
+
+struct HighlightInt
+    i::Int
+    highlight::Union{Symbol, Nothing}
+end
+Base.typeinfo_implicit(::Type{HighlightInt}) = true
+
+function Base.show(io::IO, hi::HighlightInt)
+    if hi.highlight !== nothing
+        printstyled(io, hi.i, color = hi.highlight)
+    else
+        print(io, hi.i)
+    end
 end
+
 function Base.show(io::IO, l::BipartiteAdjacencyList)
     if l.u === nothing
         printstyled(io, '⋅', color = :light_black)
     elseif isempty(l.u)
         printstyled(io, '∅', color = :light_black)
-    else
+    elseif l.highligh_u === nothing
         print(io, l.u)
+    else
+        function choose_color(i)
+            i in l.highligh_u ? (i == l.match ? :light_yellow : :green) :
+            (i == l.match ? :yellow : nothing)
+        end
+        if !isempty(setdiff(l.highligh_u, l.u))
+            # Only for debugging, shouldn't happen in practice
+            print(io, map(union(l.u, l.highligh_u)) do i
+                      HighlightInt(i, !(i in l.u) ? :light_red : choose_color(i))
+                  end)
+        else
+            print(io, map(l.u) do i
+                      HighlightInt(i, choose_color(i))
+                  end)
+        end
     end
 end
 
@@ -435,19 +480,23 @@ function set_neighbors!(g::BipartiteGraph, i::Integer, new_neighbors)
         for n in old_neighbors
             @inbounds list = g.badjlist[n]
             index = searchsortedfirst(list, i)
-            deleteat!(list, index)
+            if 1 <= index <= length(list) && list[index] == i
+                deleteat!(list, index)
+            end
         end
         for n in new_neighbors
             @inbounds list = g.badjlist[n]
             index = searchsortedfirst(list, i)
-            insert!(list, index, i)
+            if !(1 <= index <= length(list) && list[index] == i)
+                insert!(list, index, i)
+            end
         end
     end
     if iszero(new_nneighbors) # this handles Tuple as well
         # Warning: Aliases old_neighbors
         empty!(g.fadjlist[i])
     else
-        g.fadjlist[i] = copy(new_neighbors)
+        g.fadjlist[i] = unique!(sort(new_neighbors))
     end
 end
 

src/debugging.jl

@@ -0,0 +1,36 @@
+const LOGGED_FUN = Set([log, sqrt, (^), /, inv])
+is_legal(::typeof(/), a, b) = is_legal(inv, b)
+is_legal(::typeof(inv), a) = !iszero(a)
+is_legal(::Union{typeof(log), typeof(sqrt)}, a) = a isa Complex || a >= zero(a)
+is_legal(::typeof(^), a, b) = a isa Complex || b isa Complex || isinteger(b) || a >= zero(a)
+
+struct LoggedFun{F}
+    f::F
+    args::Any
+end
+Base.nameof(lf::LoggedFun) = nameof(lf.f)
+SymbolicUtils.promote_symtype(::LoggedFun, Ts...) = Real
+function (lf::LoggedFun)(args...)
+    f = lf.f
+    symbolic_args = lf.args
+    if is_legal(f, args...)
+        f(args...)
+    else
+        args_str = join(string.(symbolic_args .=> args), ", ", ", and ")
+        throw(DomainError(args, "$(lf.f) errors with input(s): $args_str"))
+    end
+end
+
+function logged_fun(f, args...)
+    # Currently we don't really support complex numbers
+    term(LoggedFun(f, args), args..., type = Real)
+end
+
+debug_sub(eq::Equation) = debug_sub(eq.lhs) ~ debug_sub(eq.rhs)
+function debug_sub(ex)
+    istree(ex) || return ex
+    f = operation(ex)
+    args = map(debug_sub, arguments(ex))
+    f in LOGGED_FUN ? logged_fun(f, args...) :
+    similarterm(ex, f, args, metadata = metadata(ex))
+end

src/inputoutput.jl

@@ -160,9 +160,9 @@ has_var(ex, x) = x ∈ Set(get_variables(ex))
 # Build control function
 
 """
-    (f_oop, f_ip), dvs, p = generate_control_function(sys::AbstractODESystem, dvs = states(sys), ps = parameters(sys); implicit_dae = false, ddvs = if implicit_dae
+    (f_oop, f_ip), dvs, p = generate_control_function(sys::AbstractODESystem, inputs = unbound_inputs(sys); implicit_dae = false, ddvs = if implicit_dae
 
-For a system `sys` that has unbound inputs (as determined by [`unbound_inputs`](@ref)), generate a function with additional input argument `in`
+For a system `sys` with inputs (as determined by [`unbound_inputs`](@ref) or user specified), generate a function with additional input argument `in`
 ```
 f_oop : (u,in,p,t)      -> rhs
 f_ip  : (uout,u,in,p,t) -> nothing
@@ -187,7 +187,7 @@ function generate_control_function(sys::AbstractODESystem, inputs = unbound_inpu
         error("No unbound inputs were found in system.")
     end
 
-    sys, diff_idxs, alge_idxs = io_preprocessing(sys, inputs, []; simplify, kwargs...)
+    sys, _ = io_preprocessing(sys, inputs, []; simplify, kwargs...)
 
     dvs = states(sys)
     ps = parameters(sys)
@@ -269,8 +269,8 @@ function inputs_to_parameters!(state::TransformationState, io)
         @assert new_v > 0
         new_var_to_diff[new_i] = new_v
     end
-    @set! structure.var_to_diff = new_var_to_diff
-    @set! structure.graph = new_graph
+    @set! structure.var_to_diff = complete(new_var_to_diff)
+    @set! structure.graph = complete(new_graph)
 
     @set! sys.eqs = map(Base.Fix2(substitute, input_to_parameters), equations(sys))
     @set! sys.states = setdiff(states(sys), keys(input_to_parameters))

src/structural_transformation/StructuralTransformations.jl

@@ -21,7 +21,8 @@ using ModelingToolkit: ODESystem, AbstractSystem, var_from_nested_derivative, Di
                        ExtraVariablesSystemException,
                        get_postprocess_fbody, vars!,
                        IncrementalCycleTracker, add_edge_checked!, topological_sort,
-                       invalidate_cache!, Substitutions, get_or_construct_tearing_state
+                       invalidate_cache!, Substitutions, get_or_construct_tearing_state,
+                       AliasGraph
 
 using ModelingToolkit.BipartiteGraphs
 import .BipartiteGraphs: invview

src/structural_transformation/pantelides.jl

@@ -74,16 +74,47 @@ end
 
 Perform Pantelides algorithm.
 """
-function pantelides!(state::TransformationState; maxiters = 8000)
+function pantelides!(state::TransformationState, ag::Union{AliasGraph, Nothing} = nothing;
+                     maxiters = 8000)
     @unpack graph, solvable_graph, var_to_diff, eq_to_diff = state.structure
     neqs = nsrcs(graph)
     nvars = nv(var_to_diff)
     vcolor = falses(nvars)
     ecolor = falses(neqs)
     var_eq_matching = Matching(nvars)
     neqs′ = neqs
+    nnonemptyeqs = count(eq -> !isempty(𝑠neighbors(graph, eq)) && eq_to_diff[eq] === nothing,
+                         1:neqs′)
+
+    # Allow matching for the highest differentiated variable that
+    # currently appears in an equation (or implicit equation in a side ag)
+    varwhitelist = falses(nvars)
+    for var in 1:nvars
+        if var_to_diff[var] === nothing && !varwhitelist[var]
+            while isempty(𝑑neighbors(graph, var)) && (ag === nothing || !haskey(ag, var))
+                var′ = invview(var_to_diff)[var]
+                var′ === nothing && break
+                var = var′
+            end
+            if !isempty(𝑑neighbors(graph, var))
+                if ag !== nothing && haskey(ag, var)
+                    # TODO: remove lower diff vars from whitelist
+                    c, a = ag[var]
+                    iszero(c) || (varwhitelist[a] = true)
+                else
+                    varwhitelist[var] = true
+                end
+            end
+        end
+    end
+
+    if nnonemptyeqs > count(varwhitelist)
+        throw(InvalidSystemException("System is structurally singular"))
+    end
+
     for k in 1:neqs′
         eq′ = k
+        eq_to_diff[eq′] === nothing || continue
         isempty(𝑠neighbors(graph, eq′)) && continue
         pathfound = false
         # In practice, `maxiters=8000` should never be reached, otherwise, the
@@ -93,7 +124,6 @@ function pantelides!(state::TransformationState; maxiters = 8000)
             #
             # the derivatives and algebraic variables are zeros in the variable
             # association list
-            varwhitelist = var_to_diff .== nothing
             resize!(vcolor, nvars)
             fill!(vcolor, false)
             resize!(ecolor, neqs)
@@ -103,11 +133,16 @@ function pantelides!(state::TransformationState; maxiters = 8000)
             pathfound && break # terminating condition
             for var in eachindex(vcolor)
                 vcolor[var] || continue
-                # introduce a new variable
-                nvars += 1
-                var_diff = var_derivative!(state, var)
-                push!(var_eq_matching, unassigned)
-                @assert length(var_eq_matching) == var_diff
+                if var_to_diff[var] === nothing
+                    # introduce a new variable
+                    nvars += 1
+                    var_diff = var_derivative!(state, var)
+                    push!(var_eq_matching, unassigned)
+                    push!(varwhitelist, false)
+                    @assert length(var_eq_matching) == var_diff
+                end
+                varwhitelist[var] = false
+                varwhitelist[var_to_diff[var]] = true
             end
 
             for eq in eachindex(ecolor)