Merge pull request #1575 from SciML/myb/pss

Add the dummy derivative pass

Author: YingboMa

docs/src/mtkitize_tutorials/modelingtoolkitize_index_reduction.md

@@ -29,8 +29,8 @@ 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(traced_sys)
-prob = ODAEProblem(pendulum_sys, Pair[], tspan)
+pendulum_sys = structural_simplify(dae_index_lowering(traced_sys))
+prob = ODAEProblem(pendulum_sys, [], tspan)
 sol = solve(prob, Tsit5(),abstol=1e-8,reltol=1e-8)
 plot(sol, vars=states(traced_sys))
 ```

src/ModelingToolkit.jl

@@ -175,7 +175,7 @@ export ControlSystem
 export alias_elimination, flatten
 export connect, @connector, Connection, Flow, Stream, instream
 export isinput, isoutput, getbounds, hasbounds, isdisturbance, istunable, getdist, hasdist, tunable_parameters
-export ode_order_lowering, liouville_transform
+export ode_order_lowering, dae_order_lowering, liouville_transform
 export runge_kutta_discretize
 export PDESystem
 export Differential, expand_derivatives, @derivatives

src/bipartite_graph.jl

@@ -234,13 +234,13 @@ Base.in(edge::BipartiteEdge, g::BipartiteGraph) = Graphs.has_edge(g, edge)
 
 ### Maximal matching
 """
-    construct_augmenting_path!(m::Matching, g::BipartiteGraph, vsrc, dstfilter, vcolor=falses(ndsts(g)), ecolor=falses(nsrcs(g))) -> path_found::Bool
+    construct_augmenting_path!(m::Matching, g::BipartiteGraph, vsrc, dstfilter, vcolor=falses(ndsts(g)), ecolor=nothing) -> path_found::Bool
 
 Try to construct an augmenting path in matching and if such a path is found,
 update the matching accordingly.
 """
-function construct_augmenting_path!(matching::Matching, g::BipartiteGraph, vsrc, dstfilter, dcolor=falses(ndsts(g)), scolor=falses(nsrcs(g)))
-    scolor[vsrc] = true
+function construct_augmenting_path!(matching::Matching, g::BipartiteGraph, vsrc, dstfilter, dcolor=falses(ndsts(g)), scolor=nothing)
+    scolor === nothing || (scolor[vsrc] = true)
 
     # if a `vdst` is unassigned and the edge `vsrc <=> vdst` exists
     for vdst in 𝑠neighbors(g, vsrc)

src/inputoutput.jl

@@ -210,7 +210,7 @@ function generate_control_function(
     inputs = map(x->time_varying_as_func(value(x), sys), ctrls)
 
     eqs = [eq for eq in equations(sys) if !isdifferenceeq(eq)]
-    foreach(check_derivative_variables, eqs)
+    check_operator_variables(eqs, Differential)
     # substitute x(t) by just x
     rhss = implicit_dae ? [_iszero(eq.lhs) ? eq.rhs : eq.rhs - eq.lhs for eq in eqs] :
                         [eq.rhs for eq in eqs]
@@ -244,5 +244,5 @@ function toparam(sys, ctrls::AbstractVector)
     eqs = map(eqs) do eq
         substitute(eq.lhs, subs) ~ substitute(eq.rhs, subs)
     end
-    ODESystem(eqs, name=sys.name)
-end
+    ODESystem(eqs, name=nameof(sys))
+end

src/structural_transformation/StructuralTransformations.jl

@@ -27,7 +27,7 @@ using ModelingToolkit.BipartiteGraphs
 import .BipartiteGraphs: invview
 using Graphs
 using ModelingToolkit.SystemStructures
-using ModelingToolkit.SystemStructures: algeqs
+using ModelingToolkit.SystemStructures: algeqs, EquationsView
 
 using ModelingToolkit.DiffEqBase
 using ModelingToolkit.StaticArrays
@@ -40,6 +40,7 @@ using SparseArrays
 using NonlinearSolve
 
 export tearing, partial_state_selection, dae_index_lowering, check_consistency
+export dummy_derivative
 export build_torn_function, build_observed_function, ODAEProblem
 export sorted_incidence_matrix, pantelides!, tearing_reassemble, find_solvables!
 export tearing_assignments, tearing_substitution

src/structural_transformation/bareiss.jl

@@ -222,7 +222,7 @@ function nullspace(A; col_order=nothing)
             @swap(col_order[i],col_order[cp])
         end
     end
-    
+
     N = ModelingToolkit.reduced_echelon_nullspace(rank, B)
     apply_inv_pivot_rows!(N, column_pivots)
 end

src/structural_transformation/bipartite_tearing/modia_tearing.jl

@@ -1,18 +1,6 @@
-# This code is from the Modia project and is licensed as follows:
+# This code is derived from the Modia project and is licensed as follows:
 # https://github.com/ModiaSim/Modia.jl/blob/b61daad643ef7edd0c1ccce6bf462c6acfb4ad1a/LICENSE
 
-################################################
-#
-# Functions to tear systems of equations
-#
-# Author: Martin Otter, DLR-SR (first version: Jan. 14, 2017)
-#
-# Details are described in the paper:
-#   Otter, Elmqvist (2017): Transformation of Differential Algebraic Array Equations to
-#                           Index One Form. Modelica'2017 Conference.
-#
-################################################
-
 function try_assign_eq!(ict::IncrementalCycleTracker, vj::Integer, eq::Integer)
     G = ict.graph
     add_edge_checked!(ict, Iterators.filter(!=(vj), 𝑠neighbors(G.graph, eq)), vj) do G
@@ -21,19 +9,6 @@ function try_assign_eq!(ict::IncrementalCycleTracker, vj::Integer, eq::Integer)
     end
 end
 
-"""
-    (eSolved, vSolved, eResidue, vTear) = tearEquations!(td, Gsolvable, es, vs)
-
-Equations es shall be solved with respect to variables vs. The function returns
-the teared equation so that if vTear is given, vSolved can be computed from eSolved
-in a forward sequence (so solving eSolved[1] for vSolved[1], eSolved[2] for vSolved[2],
-and so on). vTear must be selected, so that the equations eResidues are fulfilled.
-Equations es are the union of eSolved and eResidue.
-Variables vs are the union of vSolved and vTear.
-
-Gsolvable defines the variables that can be explicitly solved in every equation without influencing the solution space
-(= rank preserving operation).
-"""
 function tearEquations!(ict::IncrementalCycleTracker, Gsolvable, es::Vector{Int}, vs::Vector{Int})
     G = ict.graph
     vActive = BitSet(vs)
@@ -59,13 +34,8 @@ function tear_graph_block_modia!(var_eq_matching, graph, solvable_graph, eqs, va
     return nothing
 end
 
-"""
-    tear_graph_modia(sys) -> sys
-
-Tear the bipartite graph in a system. End users are encouraged to call [`structural_simplify`](@ref)
-instead, which calls this function internally.
-"""
-function tear_graph_modia(graph::BipartiteGraph, solvable_graph::BipartiteGraph; varfilter=v->true, eqfilter=eq->true)
+function tear_graph_modia(structure::SystemStructure; varfilter=v->true, eqfilter=eq->true)
+    @unpack graph, solvable_graph = structure
     var_eq_matching = complete(maximal_matching(graph, eqfilter, varfilter))
     var_sccs::Vector{Union{Vector{Int}, Int}} = find_var_sccs(graph, var_eq_matching)
 

src/structural_transformation/codegen.jl

@@ -524,8 +524,11 @@ function ODAEProblem{iip}(
                           parammap = DiffEqBase.NullParameters();
                           callback = nothing,
                           use_union = false,
+                          check = true,
                           kwargs...
                          ) where {iip}
+    eqs = equations(sys)
+    check && ModelingToolkit.check_operator_variables(eqs, Differential)
     fun, dvs = build_torn_function(sys; kwargs...)
     ps = parameters(sys)
     defs = defaults(sys)
@@ -535,7 +538,7 @@ function ODAEProblem{iip}(
     u0 = ModelingToolkit.varmap_to_vars(u0map, dvs; defaults=defs, tofloat=true)
     p = ModelingToolkit.varmap_to_vars(parammap, ps; defaults=defs, tofloat=!use_union, use_union)
 
-    has_difference = any(isdifferenceeq, equations(sys))
+    has_difference = any(isdifferenceeq, eqs)
     if has_continuous_events(sys)
         event_cb = generate_rootfinding_callback(sys; kwargs...)
     else

src/structural_transformation/partial_state_selection.jl

@@ -139,3 +139,98 @@ function partial_state_selection_graph!(structure::SystemStructure, var_eq_match
 
     var_eq_matching
 end
+
+function dummy_derivative_graph!(state::TransformationState, jac=nothing)
+    var_eq_matching = complete(pantelides!(state))
+    complete!(state.structure)
+    dummy_derivative_graph!(state.structure, var_eq_matching, jac)
+end
+
+function dummy_derivative_graph!(structure::SystemStructure, var_eq_matching, jac)
+    @unpack eq_to_diff, var_to_diff, graph = structure
+    diff_to_eq = invview(eq_to_diff)
+    diff_to_var = invview(var_to_diff)
+    invgraph = invview(graph)
+
+    neqs = nsrcs(graph)
+    eqlevel = zeros(Int, neqs)
+    maxlevel = 0
+    for i in 1:neqs
+        level = 0
+        eq = i
+        while diff_to_eq[eq] !== nothing
+            eq = diff_to_eq[eq]
+            level += 1
+        end
+        maxlevel = max(maxlevel, level)
+        eqlevel[i] = level
+    end
+
+    nvars = ndsts(graph)
+    varlevel = zeros(Int, nvars)
+    for i in 1:nvars
+        level = 0
+        var = i
+        while diff_to_var[var] !== nothing
+            var = diff_to_var[var]
+            level += 1
+        end
+        maxlevel = max(maxlevel, level)
+        varlevel[i] = level
+    end
+
+    var_sccs = find_var_sccs(graph, var_eq_matching)
+    eqcolor = falses(neqs)
+    dummy_derivatives = Int[]
+    col_order = Int[]
+    for vars in var_sccs
+        eqs = [var_eq_matching[var] for var in vars if var_eq_matching[var] !== unassigned]
+        isempty(eqs) && continue
+        maxlevel = maximum(map(x->eqlevel[x], eqs))
+        iszero(maxlevel) && continue
+
+        rank_matching = Matching(nvars)
+        for level in maxlevel:-1:1
+            eqs = filter(eq->diff_to_eq[eq] !== nothing, eqs)
+            nrows = length(eqs)
+            iszero(nrows) && break
+            eqs_set = BitSet(eqs)
+
+            # TODO: making the algorithm more robust
+            # 1. If the Jacobian is a integer matrix, use Bareiss to check
+            # linear independence. (done)
+            #
+            # 2. If the Jacobian is a single row, generate pivots. (Dynamic
+            # state selection.)
+            #
+            # 3. If the Jacobian is a polynomial matrix, use Gröbner basis (?)
+            if jac !== nothing && (_J = jac(eqs, vars); all(x->unwrap(x) isa Integer, _J))
+                J = Int.(unwrap.(_J))
+                N = ModelingToolkit.nullspace(J; col_order) # modifies col_order
+                rank = length(col_order)-size(N, 2)
+                for i in 1:rank
+                    push!(dummy_derivatives, vars[col_order[i]])
+                end
+            else
+                rank = 0
+                for var in vars
+                    pathfound = construct_augmenting_path!(rank_matching, invgraph, var, eq->eq in eqs_set, eqcolor)
+                    pathfound || continue
+                    push!(dummy_derivatives, var)
+                    rank += 1
+                    rank == nrows && break
+                end
+                fill!(rank_matching, unassigned)
+            end
+            if rank != nrows
+                @warn "The DAE system is structurally singular!"
+            end
+
+            # prepare the next iteration
+            eqs = map(eq->diff_to_eq[eq], eqs)
+            vars = [diff_to_var[var] for var in vars if diff_to_var[var] !== nothing]
+        end
+    end
+
+    dummy_derivatives
+end

src/structural_transformation/symbolics_tearing.jl

@@ -248,7 +248,7 @@ function tearing(state::TearingState; kwargs...)
     @unpack graph, solvable_graph = state.structure
     algvars = BitSet(findall(v->isalgvar(state.structure, v), 1:ndsts(graph)))
     aeqs = algeqs(state.structure)
-    var_eq_matching = Matching{Union{Unassigned, SelectedState}}(tear_graph_modia(graph, solvable_graph;
+    var_eq_matching = Matching{Union{Unassigned, SelectedState}}(tear_graph_modia(state.structure;
         varfilter=var->var in algvars, eqfilter=eq->eq in aeqs))
     for var in 1:ndsts(graph)
         if isdiffvar(state.structure, var)
@@ -282,3 +282,21 @@ function partial_state_selection(sys; simplify=false)
 
     tearing_reassemble(state, var_eq_matching; simplify=simplify)
 end
+
+"""
+    dummy_derivative(sys)
+
+Perform index reduction and use the dummy derivative techinque to ensure that
+the system is balanced.
+"""
+function dummy_derivative(sys, state=TearingState(sys))
+    function jac(eqs, vars)
+        symeqs = EquationsView(state)[eqs]
+        Symbolics.jacobian((x->x.rhs).(symeqs), state.fullvars[vars])
+    end
+    dds = dummy_derivative_graph!(state, jac)
+    symdds = Symbolics.diff2term.(state.fullvars[dds])
+    subs = Dict(state.fullvars[dd] => symdds[i] for (i, dd) in enumerate(dds))
+    @set! sys.eqs = substitute.(EquationsView(state), (subs,))
+    @set! sys.states = [states(sys); symdds]
+end