From fc815cbffb32577fc86446b36b0970108d1cbe39 Mon Sep 17 00:00:00 2001 From: Yingbo Ma Date: Sat, 27 Feb 2021 21:35:19 -0500 Subject: [PATCH] Move StructuralTransformations to MTK --- LICENSE.md | 26 ++ Project.toml | 2 + src/ModelingToolkit.jl | 4 +- src/bipartite_tearing/modia_tearing.jl | 380 ++++++++++++++++++ .../StructuralTransformations.jl | 44 ++ .../bipartite_tearing/modia_tearing.jl | 380 ++++++++++++++++++ src/structural_transformation/codegen.jl | 351 ++++++++++++++++ src/structural_transformation/pantelides.jl | 174 ++++++++ src/structural_transformation/tearing.jl | 220 ++++++++++ src/structural_transformation/utils.jl | 259 ++++++++++++ test/runtests.jl | 1 + test/structural_transformation/components.jl | 104 +++++ .../index_reduction.jl | 179 +++++++++ test/structural_transformation/runtests.jl | 6 + test/structural_transformation/tearing.jl | 181 +++++++++ test/structural_transformation/utils.jl | 39 ++ 16 files changed, 2349 insertions(+), 1 deletion(-) create mode 100644 src/bipartite_tearing/modia_tearing.jl create mode 100644 src/structural_transformation/StructuralTransformations.jl create mode 100644 src/structural_transformation/bipartite_tearing/modia_tearing.jl create mode 100644 src/structural_transformation/codegen.jl create mode 100644 src/structural_transformation/pantelides.jl create mode 100644 src/structural_transformation/tearing.jl create mode 100644 src/structural_transformation/utils.jl create mode 100644 test/structural_transformation/components.jl create mode 100644 test/structural_transformation/index_reduction.jl create mode 100644 test/structural_transformation/runtests.jl create mode 100644 test/structural_transformation/tearing.jl create mode 100644 test/structural_transformation/utils.jl diff --git a/LICENSE.md b/LICENSE.md index 9697030102..3a7add8ed1 100644 --- a/LICENSE.md +++ b/LICENSE.md @@ -38,3 +38,29 @@ The ModelingToolkit.jl package is licensed under the MIT "Expat" License: > SOFTWARE. > > + +The code in `src/structural_transformation/bipartite_tearing/modia_tearing.jl`, +which is from the [Modia.jl](https://github.com/ModiaSim/Modia.jl) project, is +licensed as follows: + +MIT License + +Copyright (c) 2017-2018 ModiaSim developers + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. diff --git a/Project.toml b/Project.toml index e4c9feb202..139a508a18 100644 --- a/Project.toml +++ b/Project.toml @@ -21,6 +21,7 @@ LightGraphs = "093fc24a-ae57-5d10-9952-331d41423f4d" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09" NaNMath = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3" +NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec" RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" Requires = "ae029012-a4dd-5104-9daa-d747884805df" @@ -51,6 +52,7 @@ Latexify = "0.11, 0.12, 0.13, 0.14" LightGraphs = "1.3" MacroTools = "0.5" NaNMath = "0.3" +NonlinearSolve = "0.3.8" RecursiveArrayTools = "2.3" Reexport = "0.2, 1" Requires = "1.0" diff --git a/src/ModelingToolkit.jl b/src/ModelingToolkit.jl index beed2cadac..4c7fbceaad 100644 --- a/src/ModelingToolkit.jl +++ b/src/ModelingToolkit.jl @@ -1,6 +1,6 @@ module ModelingToolkit -using DiffEqBase, SciMLBase +using DiffEqBase, SciMLBase, Reexport using Distributed using StaticArrays, LinearAlgebra, SparseArrays, LabelledArrays using Latexify, Unitful, ArrayInterface @@ -249,6 +249,8 @@ include("systems/systemstructure.jl") using .SystemStructures include("systems/alias_elimination.jl") +include("structural_transformation/StructuralTransformations.jl") +@reexport using .StructuralTransformations include("latexify_recipes.jl") include("build_function.jl") diff --git a/src/bipartite_tearing/modia_tearing.jl b/src/bipartite_tearing/modia_tearing.jl new file mode 100644 index 0000000000..35b555fdb5 --- /dev/null +++ b/src/bipartite_tearing/modia_tearing.jl @@ -0,0 +1,380 @@ +# This code is 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. +# +# The following utility algorithm is used below to incrementally added edges to a +# DAG (Directed Acyclic Graph). This algorithm leads to an O(n*m) worst time complexity of the +# tearing (instead of O(m*m)) where n is the number of equations and m is the number of +# variable incidences. Note, the traversals of the DAG are not performed with recursive function +# calls but with while loops and an explicit stack, in order to avoid function stack overflow +# for large algebraic loops. +# +# Bender, Fineman, Gilbert, Tarjan: +# A New Approach to Incremental Cycle Detection and Related Problems. +# ACM Transactions on Algorithms, Volume 12, Issue 2, Feb. 2016 +# http://dl.acm.org/citation.cfm?id=2756553 +# +# Text excerpt from this paper (the advantage of Algorithm N is that it +# is simple to implement and needs no sophisticated data structures) +# +# 3. A ONE-WAY-SEARCH ALGORITHM FOR DENSE GRAPHS +# [..] +# To develop the algorithm, we begin with a simple algorithm and then modify it to +# improve its running time. We call the simple algorithm Algorithm N (for “naive”). The +# algorithm initializes k(v) = 1 for each vertex v. The initial vertex levels are a weak +# topological numbering since the initial graph contains no arcs. To add an arc (v,w), +# if k(v) < k(w), merely add the arc. If, on the other hand, k(v) ≥ k(w), add the arc and +# then do a selective forward search that begins by traversing (v,w) and continues until +# the search traverses an arc into v (there is a cycle) or there are no more candidate +# arcs to traverse. To traverse an arc (x, y), if y = v, stop (there is a cycle); otherwise, if +# k(x) ≥ k(y), increase k(y) to k(x)+1 and add all arcs (y, z) to the set of candidate arcs to +# traverse. +# It is easy to show that (1) after each arc addition that does not form a cycle +# the vertex levels are a weak topological numbering, (2) Algorithm N is correct, and +# (3) 1 ≤ k(v) ≤ size(v) ≤ n for all v. Since an arc (x, y) notEqual (v,w) is only traversed as a +# result of k(x) increasing, each arc is traversed at most n times, resulting in a total time +# bound of O(nm) for all arc additions. +# +################################################ + +const Undefined = typemin(Int) + + +""" + td = TraverseDAG(G,nv) + +Generate an object td to traverse a set of equations that are +represented as a DAG (Directed Acyclic Graph). +G is the bipartite graph of all relevant equations +and nv is the largest variable number used in G (or larger). +""" +mutable struct TraverseDAG + minlevel::Int + curlevel::Int + level::Vector{Int} + lastlevel::Vector{Int} + levelStack::Vector{Int} + visitedStack::Vector{Int} + vActive::Vector{Bool} + visited::Vector{Bool} + check::Vector{Bool} + stack::Vector{Int} + eSolved::Vector{Int} + vSolved::Vector{Int} + G # Vector{ Vector{Int} } + assign::Vector{Int} + es::Vector{Int} + vs::Vector{Int} + + function TraverseDAG(G, nv::Int) + visited = fill(false, length(G)) + check = fill(false, length(G)) + vActive = fill(false, nv) + level = fill(Undefined, nv) + lastlevel = fill(Undefined, nv) + levelStack = fill(0, 0) + stack = fill(0, 0) + visitedStack = fill(0, 0) + eSolved = fill(0, 0) + vSolved = fill(0, 0) + assign = fill(0, nv) + + new(0, Undefined, level, lastlevel, levelStack, visitedStack, vActive, + visited, check, stack, eSolved, vSolved, G, assign) + end +end + + +""" + initAlgebraicSystem(td::TraverseDAG,es,vs) + +Define the set of equations and the set variables for which the equations shall be solved for +(equations es shall be solved for variables vs) and re-initialize td. + +eSolvedFixed/vSolvedFixed must be a DAG starting at eSolvedFixed/SolvedFixed[1] +""" +function initAlgebraicSystem(td::TraverseDAG, es::Vector{Int}, vs::Vector{Int}, + eSolvedFixed::Vector{Int}, vSolvedFixed::Vector{Int}, vTearFixed::Vector{Int}) + # check arguments + for i in eachindex(es) + if es[i] <= 0 + error("\n\n... Internal error in Tearing.jl: es[", i, "] = ", es[i], ".\n") + end + end + + for i in eachindex(vs) + if vs[i] <= 0 + error("\n\n... Internal error in Tearing.jl: vs[", i, "] = ", vs[i], ".\n") + end + end + + # check that all elements of eSolvedFixed are in es, vSolvedFixed in vs, + # vTearFixed in vs and that vSolvedFixed and vTearFixed have no variables in common + @assert( length(eSolvedFixed) == length(vSolvedFixed) ) + ediff = setdiff(eSolvedFixed, es) + @assert(length(ediff) == 0) + vdiff = setdiff(vSolvedFixed, vs) + @assert(length(vdiff) == 0) + vdiff2 = setdiff(vTearFixed, vs) + @assert(length(vdiff2) == 0) + vdiff3 = intersect(vSolvedFixed, vTearFixed) + @assert(length(vdiff3) == 0) + + # Re-initialize td + td.minlevel = 0 + td.curlevel = Undefined + for i in eachindex(td.visited) + td.visited[i] = false + td.check[i] = false + end + + for i in eachindex(td.vActive) + td.vActive[i] = false + td.assign[i] = 0 + td.level[i] = Undefined + td.lastlevel[i] = Undefined + end + + for i in eachindex(vs) + td.vActive[ vs[i] ] = true + end + + empty!(td.levelStack) + empty!(td.stack) + empty!(td.visitedStack) + empty!(td.eSolved) + empty!(td.vSolved) + + # Define initial DAG + vs2 = Int[] + for i in eachindex(vSolvedFixed) + vFixed = vSolvedFixed[i] + td.assign[vFixed] = eSolvedFixed[i] + td.level[ vFixed] = i + push!(vs2, vFixed) + end + + for i in eachindex(vTearFixed) + td.vActive[ vTearFixed[i] ] = false # vTearFixed shall not be assigned + end + + # Store es, vs in td + td.es = es + td.vs = vs + + return vs2 +end + + +in_vs(td, v) = td.vActive[v] + +function setlevel(td::TraverseDAG, v::Int, parentLevel::Int) + td.lastlevel[v] = td.level[v] + td.level[v] = parentLevel + 1 + push!(td.visitedStack, v) +end + + +""" + success = visit!(td::TraverseDAG, v) + +Traverse potential DAG starting from new variable node v. +If no cycle is detected return true, otherwise return false. +""" +function visit!(td::TraverseDAG, vcheck::Int) + empty!(td.stack) + empty!(td.levelStack) + empty!(td.visitedStack) + td.curlevel = td.level[vcheck] + push!(td.levelStack, td.curlevel) + push!(td.stack, vcheck) + first = true + + while length(td.stack) > 0 + parentLevel = pop!(td.levelStack) + veq = pop!(td.stack) + eq = td.assign[veq] + if first + first = false + else + if td.level[veq] == td.curlevel + # cycle detected + return false + elseif td.level[veq] == Undefined || td.level[veq] <= parentLevel + setlevel(td, veq, parentLevel) + end + end + + if eq > 0 + # Push all child nodes on stack + parentLevel = td.level[veq] + for v in td.G[eq] + if in_vs(td, v) && v != veq # v is an element of td.vs and is not the variable to solve for + eq2 = td.assign[v] + if eq2 == 0 || td.level[v] <= parentLevel + push!(td.levelStack, parentLevel) + push!(td.stack, v) + end + end + end + end + end + + return true +end + + +""" + visit2!(td::TraverseDAG,v) + +Traverse DAG starting from variable v and store visited equations and variables in stacks +eSolved, vSolved. If a cycle is deteced, raise an error (signals a programming error). +""" +function visit2!(td::TraverseDAG, vVisit::Int) + push!(td.stack, vVisit) + while length(td.stack) > 0 + veq = td.stack[end] + eq = td.assign[veq] + if !td.visited[eq] + td.visited[eq] = true + td.check[eq] = true + for v in td.G[eq] + if in_vs(td, v) && v != veq # v is an element of td.vs and is not the variable to solve for + eq2 = td.assign[v] + if eq2 != 0 + if !td.visited[eq2] # visit eq2 if not yet visited + push!(td.stack, v) + elseif td.check[eq2] # cycle detected + error("... error in Tearing.jl code: \n", + " cycle detected (should not occur): eq = ", eq, ", veq = ", veq, ", eq2 = ", eq2, ", v = ", v) + end + end + end + end + else + td.check[eq] = false + push!(td.eSolved, eq) + push!(td.vSolved, veq) + pop!(td.stack) + end + end + nothing +end + + +""" + (eSolved, vSolved) = sortDAG!(td::TraverseDAG, vs) + +Sort the equations that are assigned by variables vs using object td of type TraverseDAG +and return the sorted equations eSolved and assigned variables vSolved. +""" +function sortDAG!(td::TraverseDAG, vs::Vector{Int}) + # initialize data structure + empty!(td.stack) + empty!(td.eSolved) + empty!(td.vSolved) + + for i in eachindex(td.visited) + td.visited[i] = false + td.check[i] = false + end + + # visit all assigned variables and equations + for veq in vs + if !td.visited[ td.assign[veq] ] + visit2!(td, veq) + end + end + + return (td.eSolved, td.vSolved) +end + + +""" + (eSolved, vSolved, eResidue, vTear) = tearEquations!(td, Gsolvable, es, vs; eSolvedFixed=Int[], vSolvedFixed=Int[], vTearFixed=Int[]) + +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. + +Input argument td is an object of type TraverseDAG. Gsolvable defines the variables +that can be explicitly solved in every equation without influencing the solution space +(= rank preserving operation). + +eSolvedFixed/vSolvedFixed must be a DAG starting at eSolvedFixed/SolvedFixed[1] +""" +function tearEquations!(td::TraverseDAG, Gsolvable, es::Vector{Int}, vs::Vector{Int}; + eSolvedFixed::Vector{Int}=Int[], vSolvedFixed::Vector{Int}=Int[], vTearFixed::Vector{Int}=Int[]) + vs2 = initAlgebraicSystem(td, es, vs, eSolvedFixed, vSolvedFixed, vTearFixed) + # eResidue = fill(0,0) + residue = true + esReduced = setdiff(es, eSolvedFixed) + # println(" es = ", es, ", eSolvedFixed = ", eSolvedFixed, ", esReduced = ", esReduced) + # println(" vs = ", vs, ", vSolvedFixed = ", vSolvedFixed) + for eq in esReduced # iterate only over equations that are not in eSolvedFixed + residue = true + for vj in Gsolvable[eq] + if td.assign[vj] == 0 && in_vs(td, vj) + # vj is an element of vs that is not yet assigned + # Add equation to graph + td.assign[vj] = eq + + # Check for cycles + if td.level[vj] == Undefined + # (eq,vj) cannot introduce a cycle + # Introduce a new level (the smallest level that exists yet) + td.minlevel += -1 + td.level[vj] = td.minlevel + + # Inspect all childs and use level+1, if child has no level yet + for v in td.G[eq] + if in_vs(td, v) && v != vj && + (td.level[v] == Undefined || td.level[v] <= td.level[vj]) # v is an element of td.vs and is not the variable to solve for and no level yet defined + setlevel(td, v, td.level[vj]) + end + end + + push!(vs2, vj) + residue = false + break # continue with next equation + + else # Traverse DAG starting from eq + if visit!(td, vj) + # accept vj + push!(vs2, vj) + residue = false + break # continue with next equation + else + # cycle; remove vj from DAG and undo its changes + for vv in td.visitedStack + td.level[vv] = td.lastlevel[vv] + end + td.assign[vj] = 0 + # continue with next variable in equation eq + end + end + end + end + #if residue + # push!(eResidue, eq) + #end + end + + # Determine solved equations and variables + (eSolved, vSolved) = sortDAG!(td, vs2) + vTear = setdiff(vs, vSolved) + eResidue = setdiff(es, eSolved) + return (eSolved, vSolved, eResidue, vTear) +end diff --git a/src/structural_transformation/StructuralTransformations.jl b/src/structural_transformation/StructuralTransformations.jl new file mode 100644 index 0000000000..bfdb571ea8 --- /dev/null +++ b/src/structural_transformation/StructuralTransformations.jl @@ -0,0 +1,44 @@ +module StructuralTransformations + +const UNVISITED = typemin(Int) +const UNASSIGNED = typemin(Int) + +using Setfield: @set!, @set +using UnPack: @unpack + +using SymbolicUtils +using SymbolicUtils.Code +using SymbolicUtils.Rewriters +using SymbolicUtils: similarterm, istree + +using ModelingToolkit +using ModelingToolkit: ODESystem, var_from_nested_derivative, Differential, + states, equations, vars, Symbolic, diff2term, value, + operation, arguments, Sym, Term, simplify, solve_for, + isdiffeq, isdifferential, + get_structure, default_u0, default_p + +using ModelingToolkit.BipartiteGraphs +using ModelingToolkit.SystemStructures + +using ModelingToolkit.DiffEqBase +using ModelingToolkit.StaticArrays +using ModelingToolkit: @RuntimeGeneratedFunction, RuntimeGeneratedFunctions + +RuntimeGeneratedFunctions.init(@__MODULE__) + +using SparseArrays + +using NonlinearSolve + +export tearing, dae_index_lowering +export build_torn_function, build_observed_function, ODAEProblem +export sorted_incidence_matrix + +include("utils.jl") +include("pantelides.jl") +include("bipartite_tearing/modia_tearing.jl") +include("tearing.jl") +include("codegen.jl") + +end # module diff --git a/src/structural_transformation/bipartite_tearing/modia_tearing.jl b/src/structural_transformation/bipartite_tearing/modia_tearing.jl new file mode 100644 index 0000000000..35b555fdb5 --- /dev/null +++ b/src/structural_transformation/bipartite_tearing/modia_tearing.jl @@ -0,0 +1,380 @@ +# This code is 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. +# +# The following utility algorithm is used below to incrementally added edges to a +# DAG (Directed Acyclic Graph). This algorithm leads to an O(n*m) worst time complexity of the +# tearing (instead of O(m*m)) where n is the number of equations and m is the number of +# variable incidences. Note, the traversals of the DAG are not performed with recursive function +# calls but with while loops and an explicit stack, in order to avoid function stack overflow +# for large algebraic loops. +# +# Bender, Fineman, Gilbert, Tarjan: +# A New Approach to Incremental Cycle Detection and Related Problems. +# ACM Transactions on Algorithms, Volume 12, Issue 2, Feb. 2016 +# http://dl.acm.org/citation.cfm?id=2756553 +# +# Text excerpt from this paper (the advantage of Algorithm N is that it +# is simple to implement and needs no sophisticated data structures) +# +# 3. A ONE-WAY-SEARCH ALGORITHM FOR DENSE GRAPHS +# [..] +# To develop the algorithm, we begin with a simple algorithm and then modify it to +# improve its running time. We call the simple algorithm Algorithm N (for “naive”). The +# algorithm initializes k(v) = 1 for each vertex v. The initial vertex levels are a weak +# topological numbering since the initial graph contains no arcs. To add an arc (v,w), +# if k(v) < k(w), merely add the arc. If, on the other hand, k(v) ≥ k(w), add the arc and +# then do a selective forward search that begins by traversing (v,w) and continues until +# the search traverses an arc into v (there is a cycle) or there are no more candidate +# arcs to traverse. To traverse an arc (x, y), if y = v, stop (there is a cycle); otherwise, if +# k(x) ≥ k(y), increase k(y) to k(x)+1 and add all arcs (y, z) to the set of candidate arcs to +# traverse. +# It is easy to show that (1) after each arc addition that does not form a cycle +# the vertex levels are a weak topological numbering, (2) Algorithm N is correct, and +# (3) 1 ≤ k(v) ≤ size(v) ≤ n for all v. Since an arc (x, y) notEqual (v,w) is only traversed as a +# result of k(x) increasing, each arc is traversed at most n times, resulting in a total time +# bound of O(nm) for all arc additions. +# +################################################ + +const Undefined = typemin(Int) + + +""" + td = TraverseDAG(G,nv) + +Generate an object td to traverse a set of equations that are +represented as a DAG (Directed Acyclic Graph). +G is the bipartite graph of all relevant equations +and nv is the largest variable number used in G (or larger). +""" +mutable struct TraverseDAG + minlevel::Int + curlevel::Int + level::Vector{Int} + lastlevel::Vector{Int} + levelStack::Vector{Int} + visitedStack::Vector{Int} + vActive::Vector{Bool} + visited::Vector{Bool} + check::Vector{Bool} + stack::Vector{Int} + eSolved::Vector{Int} + vSolved::Vector{Int} + G # Vector{ Vector{Int} } + assign::Vector{Int} + es::Vector{Int} + vs::Vector{Int} + + function TraverseDAG(G, nv::Int) + visited = fill(false, length(G)) + check = fill(false, length(G)) + vActive = fill(false, nv) + level = fill(Undefined, nv) + lastlevel = fill(Undefined, nv) + levelStack = fill(0, 0) + stack = fill(0, 0) + visitedStack = fill(0, 0) + eSolved = fill(0, 0) + vSolved = fill(0, 0) + assign = fill(0, nv) + + new(0, Undefined, level, lastlevel, levelStack, visitedStack, vActive, + visited, check, stack, eSolved, vSolved, G, assign) + end +end + + +""" + initAlgebraicSystem(td::TraverseDAG,es,vs) + +Define the set of equations and the set variables for which the equations shall be solved for +(equations es shall be solved for variables vs) and re-initialize td. + +eSolvedFixed/vSolvedFixed must be a DAG starting at eSolvedFixed/SolvedFixed[1] +""" +function initAlgebraicSystem(td::TraverseDAG, es::Vector{Int}, vs::Vector{Int}, + eSolvedFixed::Vector{Int}, vSolvedFixed::Vector{Int}, vTearFixed::Vector{Int}) + # check arguments + for i in eachindex(es) + if es[i] <= 0 + error("\n\n... Internal error in Tearing.jl: es[", i, "] = ", es[i], ".\n") + end + end + + for i in eachindex(vs) + if vs[i] <= 0 + error("\n\n... Internal error in Tearing.jl: vs[", i, "] = ", vs[i], ".\n") + end + end + + # check that all elements of eSolvedFixed are in es, vSolvedFixed in vs, + # vTearFixed in vs and that vSolvedFixed and vTearFixed have no variables in common + @assert( length(eSolvedFixed) == length(vSolvedFixed) ) + ediff = setdiff(eSolvedFixed, es) + @assert(length(ediff) == 0) + vdiff = setdiff(vSolvedFixed, vs) + @assert(length(vdiff) == 0) + vdiff2 = setdiff(vTearFixed, vs) + @assert(length(vdiff2) == 0) + vdiff3 = intersect(vSolvedFixed, vTearFixed) + @assert(length(vdiff3) == 0) + + # Re-initialize td + td.minlevel = 0 + td.curlevel = Undefined + for i in eachindex(td.visited) + td.visited[i] = false + td.check[i] = false + end + + for i in eachindex(td.vActive) + td.vActive[i] = false + td.assign[i] = 0 + td.level[i] = Undefined + td.lastlevel[i] = Undefined + end + + for i in eachindex(vs) + td.vActive[ vs[i] ] = true + end + + empty!(td.levelStack) + empty!(td.stack) + empty!(td.visitedStack) + empty!(td.eSolved) + empty!(td.vSolved) + + # Define initial DAG + vs2 = Int[] + for i in eachindex(vSolvedFixed) + vFixed = vSolvedFixed[i] + td.assign[vFixed] = eSolvedFixed[i] + td.level[ vFixed] = i + push!(vs2, vFixed) + end + + for i in eachindex(vTearFixed) + td.vActive[ vTearFixed[i] ] = false # vTearFixed shall not be assigned + end + + # Store es, vs in td + td.es = es + td.vs = vs + + return vs2 +end + + +in_vs(td, v) = td.vActive[v] + +function setlevel(td::TraverseDAG, v::Int, parentLevel::Int) + td.lastlevel[v] = td.level[v] + td.level[v] = parentLevel + 1 + push!(td.visitedStack, v) +end + + +""" + success = visit!(td::TraverseDAG, v) + +Traverse potential DAG starting from new variable node v. +If no cycle is detected return true, otherwise return false. +""" +function visit!(td::TraverseDAG, vcheck::Int) + empty!(td.stack) + empty!(td.levelStack) + empty!(td.visitedStack) + td.curlevel = td.level[vcheck] + push!(td.levelStack, td.curlevel) + push!(td.stack, vcheck) + first = true + + while length(td.stack) > 0 + parentLevel = pop!(td.levelStack) + veq = pop!(td.stack) + eq = td.assign[veq] + if first + first = false + else + if td.level[veq] == td.curlevel + # cycle detected + return false + elseif td.level[veq] == Undefined || td.level[veq] <= parentLevel + setlevel(td, veq, parentLevel) + end + end + + if eq > 0 + # Push all child nodes on stack + parentLevel = td.level[veq] + for v in td.G[eq] + if in_vs(td, v) && v != veq # v is an element of td.vs and is not the variable to solve for + eq2 = td.assign[v] + if eq2 == 0 || td.level[v] <= parentLevel + push!(td.levelStack, parentLevel) + push!(td.stack, v) + end + end + end + end + end + + return true +end + + +""" + visit2!(td::TraverseDAG,v) + +Traverse DAG starting from variable v and store visited equations and variables in stacks +eSolved, vSolved. If a cycle is deteced, raise an error (signals a programming error). +""" +function visit2!(td::TraverseDAG, vVisit::Int) + push!(td.stack, vVisit) + while length(td.stack) > 0 + veq = td.stack[end] + eq = td.assign[veq] + if !td.visited[eq] + td.visited[eq] = true + td.check[eq] = true + for v in td.G[eq] + if in_vs(td, v) && v != veq # v is an element of td.vs and is not the variable to solve for + eq2 = td.assign[v] + if eq2 != 0 + if !td.visited[eq2] # visit eq2 if not yet visited + push!(td.stack, v) + elseif td.check[eq2] # cycle detected + error("... error in Tearing.jl code: \n", + " cycle detected (should not occur): eq = ", eq, ", veq = ", veq, ", eq2 = ", eq2, ", v = ", v) + end + end + end + end + else + td.check[eq] = false + push!(td.eSolved, eq) + push!(td.vSolved, veq) + pop!(td.stack) + end + end + nothing +end + + +""" + (eSolved, vSolved) = sortDAG!(td::TraverseDAG, vs) + +Sort the equations that are assigned by variables vs using object td of type TraverseDAG +and return the sorted equations eSolved and assigned variables vSolved. +""" +function sortDAG!(td::TraverseDAG, vs::Vector{Int}) + # initialize data structure + empty!(td.stack) + empty!(td.eSolved) + empty!(td.vSolved) + + for i in eachindex(td.visited) + td.visited[i] = false + td.check[i] = false + end + + # visit all assigned variables and equations + for veq in vs + if !td.visited[ td.assign[veq] ] + visit2!(td, veq) + end + end + + return (td.eSolved, td.vSolved) +end + + +""" + (eSolved, vSolved, eResidue, vTear) = tearEquations!(td, Gsolvable, es, vs; eSolvedFixed=Int[], vSolvedFixed=Int[], vTearFixed=Int[]) + +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. + +Input argument td is an object of type TraverseDAG. Gsolvable defines the variables +that can be explicitly solved in every equation without influencing the solution space +(= rank preserving operation). + +eSolvedFixed/vSolvedFixed must be a DAG starting at eSolvedFixed/SolvedFixed[1] +""" +function tearEquations!(td::TraverseDAG, Gsolvable, es::Vector{Int}, vs::Vector{Int}; + eSolvedFixed::Vector{Int}=Int[], vSolvedFixed::Vector{Int}=Int[], vTearFixed::Vector{Int}=Int[]) + vs2 = initAlgebraicSystem(td, es, vs, eSolvedFixed, vSolvedFixed, vTearFixed) + # eResidue = fill(0,0) + residue = true + esReduced = setdiff(es, eSolvedFixed) + # println(" es = ", es, ", eSolvedFixed = ", eSolvedFixed, ", esReduced = ", esReduced) + # println(" vs = ", vs, ", vSolvedFixed = ", vSolvedFixed) + for eq in esReduced # iterate only over equations that are not in eSolvedFixed + residue = true + for vj in Gsolvable[eq] + if td.assign[vj] == 0 && in_vs(td, vj) + # vj is an element of vs that is not yet assigned + # Add equation to graph + td.assign[vj] = eq + + # Check for cycles + if td.level[vj] == Undefined + # (eq,vj) cannot introduce a cycle + # Introduce a new level (the smallest level that exists yet) + td.minlevel += -1 + td.level[vj] = td.minlevel + + # Inspect all childs and use level+1, if child has no level yet + for v in td.G[eq] + if in_vs(td, v) && v != vj && + (td.level[v] == Undefined || td.level[v] <= td.level[vj]) # v is an element of td.vs and is not the variable to solve for and no level yet defined + setlevel(td, v, td.level[vj]) + end + end + + push!(vs2, vj) + residue = false + break # continue with next equation + + else # Traverse DAG starting from eq + if visit!(td, vj) + # accept vj + push!(vs2, vj) + residue = false + break # continue with next equation + else + # cycle; remove vj from DAG and undo its changes + for vv in td.visitedStack + td.level[vv] = td.lastlevel[vv] + end + td.assign[vj] = 0 + # continue with next variable in equation eq + end + end + end + end + #if residue + # push!(eResidue, eq) + #end + end + + # Determine solved equations and variables + (eSolved, vSolved) = sortDAG!(td, vs2) + vTear = setdiff(vs, vSolved) + eResidue = setdiff(es, eSolved) + return (eSolved, vSolved, eResidue, vTear) +end diff --git a/src/structural_transformation/codegen.jl b/src/structural_transformation/codegen.jl new file mode 100644 index 0000000000..26da26c426 --- /dev/null +++ b/src/structural_transformation/codegen.jl @@ -0,0 +1,351 @@ +function torn_system_jacobian_sparsity(sys) + s = structure(sys) + @unpack fullvars, graph, partitions = s + + # The sparsity pattern of `nlsolve(f, u, p)` w.r.t `p` is difficult to + # determine in general. Consider the "simplest" case, a linear system. We + # have + # A u = p. + # Clearly, the sparsity of `u` depends on the sparsity of both `p` and `A` + # in a non-trivial way. However, in the generic case, `u` is dense even when + # `A` and `p` are sparse. For instance + # + # ```julia + # julia> using Random, SparseArrays + # + # julia> A = sprand(MersenneTwister(1234), 100, 100, 0.1); + # + # julia> p = sprand(MersenneTwister(12345), 100, 0.05); + # + # julia> count(x->abs(x) < 1e-5, A \ Vector(p)) + # 0 + # ``` + # + # Let 𝑇 be the set of tearing variables and 𝑉 be the set of all *states* in + # the residual equations. In the following code, we are going to assume the + # connection between 𝑇 (the `u` in from above) and 𝑉 ∖ 𝑇 (the `p` in from + # above) has full incidence. + # + # Note that as we are reducing algebraic equations numerically, it could be + # the case that a later partition (a BLT block) contains tearing variables + # from other partitions. + # + # We know that partitions are BLT ordered. Hence, the tearing variables in + # each partition is unique, and all states in a partition must be + # either differential variables or algebraic tearing variables that are + # from previous partitions. Hence, we can build the dependency chain as we + # traverse the partitions. + + # `avars2dvars` maps a algebraic variable to its differential variable + # dependencies. + avars2dvars = Dict{Int,Set{Int}}() + c = 0 + for (_, _, teqs, tvars) in partitions + # initialization + for tvar in tvars + avars2dvars[tvar] = Set{Int}() + end + for teq in teqs + c += 1 + for var in 𝑠neighbors(graph, teq) + # Skip the tearing variables in the current partition, because + # we are computing them from all the other states. + LightGraphs.insorted(var, tvars) && continue + deps = get(avars2dvars, var, nothing) + if deps === nothing # differential variable + @assert !isalgvar(s, var) + for tvar in tvars + push!(avars2dvars[tvar], var) + end + else # tearing variable from previous partitions + @assert isalgvar(s, var) + for tvar in tvars + union!(avars2dvars[tvar], avars2dvars[var]) + end + end + end + end + end + + dvar2idx(idx) = idx # maps `dvar` to the index of the states + I = Int[]; J = Int[] + eqidx = 0 + for ieq in 𝑠vertices(graph) + isalgeq(s, ieq) && continue + eqidx += 1 + for ivar in 𝑠neighbors(graph, ieq) + if isdiffvar(s, ivar) + push!(I, eqidx) + push!(J, dvar2idx(ivar)) + elseif isalgvar(s, ivar) + for dvar in avars2dvars[ivar] + push!(I, eqidx) + push!(J, dvar2idx(dvar)) + end + end + end + end + sparse(I, J, true) +end + +""" + partitions_dag(s::SystemStructure) + +Return a DAG (sparse matrix) of partitions to use for parallelism. +""" +function partitions_dag(s::SystemStructure) + @unpack partitions, graph = s + + # `partvars[i]` contains all the states that appear in `partitions[i]` + partvars = map(partitions) do (_, _, reqs, tvars) + ipartvars = Set{Int}() + for req in reqs + union!(ipartvars, 𝑠neighbors(graph, req)) + end + ipartvars + end + + I, J = Int[], Int[] + n = length(partitions) + for i in 1:n + for j in i+1:n + # The only way for a later partition `j` to depend on an earlier + # partition `i` is when `partvars[j]` contains one of tearing + # variables of partition `i`. + if !isdisjoint(partvars[j], partitions[i][4]) + # j depends on i + push!(I, i) + push!(J, j) + end + end + end + + sparse(I, J, true, n, n) +end + +function gen_nlsolve(sys, eqs, vars) + @assert !isempty(vars) + @assert length(eqs) == length(vars) + rhss = map(x->x.rhs, eqs) + # We use `vars` instead of `graph` to capture parameters, too. + allvars = unique(collect(Iterators.flatten(map(ModelingToolkit.vars, rhss)))) + params = setdiff(allvars, vars) + + u0map = default_u0(sys) + # splatting to tighten the type + u0 = [map(var->get(u0map, var, 1e-3), vars)...] + # specialize on the scalar case + isscalar = length(u0) == 1 + u0 = isscalar ? u0[1] : SVector(u0...) + + fname = gensym("fun") + f = Func( + [ + DestructuredArgs(vars) + DestructuredArgs(params) + ], + [], + isscalar ? rhss[1] : MakeArray(rhss, SVector) + ) |> SymbolicUtils.Code.toexpr + + solver_call = LiteralExpr(quote + $numerical_nlsolve( + $fname, + # initial guess + $u0, + # "captured variables" + ($(params...),) + ) + end) + + [ + fname ← @RuntimeGeneratedFunction(f) + DestructuredArgs(vars) ← solver_call + ] +end + +function get_torn_eqs_vars(sys) + s = structure(sys) + partitions = s.partitions + vars = s.fullvars + eqs = equations(sys) + + torn_eqs = map(idxs-> eqs[idxs], map(x->x[3], partitions)) + torn_vars = map(idxs->vars[idxs], map(x->x[4], partitions)) + + gen_nlsolve.((sys,), torn_eqs, torn_vars) +end + +function build_torn_function( + sys; + expression=false, + jacobian_sparsity=true, + checkbounds=false, + kw... + ) + + rhss = [] + for eq in equations(sys) + isdiffeq(eq) && push!(rhss, eq.rhs) + end + + out = Sym{Any}(gensym("out")) + odefunbody = SetArray( + checkbounds, + out, + rhss + ) + + s = structure(sys) + states = s.fullvars[diffvars_range(s)] + syms = map(Symbol, states) + + expr = SymbolicUtils.Code.toexpr( + Func( + [ + out + DestructuredArgs(states) + DestructuredArgs(parameters(sys)) + independent_variable(sys) + ], + [], + Let( + collect(Iterators.flatten(get_torn_eqs_vars(sys))), + odefunbody + ) + ) + ) + if expression + expr + else + observedfun = let sys = sys, dict = Dict() + function generated_observed(obsvar, u, p, t) + obs = get!(dict, value(obsvar)) do + build_observed_function(sys, obsvar) + end + obs(u, p, t) + end + end + + ODEFunction{true}( + @RuntimeGeneratedFunction(expr), + sparsity = torn_system_jacobian_sparsity(sys), + syms = syms, + observed = observedfun, + ) + end +end + +""" + find_solve_sequence(partitions, vars) + +given a set of `vars`, find the groups of equations we need to solve for +to obtain the solution to `vars` +""" +function find_solve_sequence(partitions, vars) + subset = filter(x -> !isdisjoint(x[4], vars), partitions) + isempty(subset) && return [] + vars′ = mapreduce(x->x[4], union, subset) + if vars′ == vars + return subset + else + return find_solve_sequence(partitions, vars′) + end +end + +function build_observed_function( + sys, syms; + expression=false, + output_type=Array + ) + + if (isscalar = !(syms isa Vector)) + syms = [syms] + end + syms = value.(syms) + syms_set = Set(syms) + s = structure(sys) + @unpack partitions, fullvars, graph = s + diffvars = fullvars[diffvars_range(s)] + algvars = fullvars[algvars_range(s)] + + required_algvars = Set(intersect(algvars, syms_set)) + obs = observed(sys) + observed_idx = Dict(map(x->x.lhs, obs) .=> 1:length(obs)) + for sym in syms + idx = get(observed_idx, sym, nothing) + idx === nothing && continue + vs = vars(obs[idx].rhs) + union!(required_algvars, intersect(algvars, vs)) + end + + varidxs = findall(x->x in required_algvars, fullvars) + subset = find_solve_sequence(partitions, varidxs) + if !isempty(subset) + eqs = equations(sys) + + torn_eqs = map(idxs-> eqs[idxs[3]], subset) + torn_vars = map(idxs->fullvars[idxs[4]], subset) + + solves = gen_nlsolve.((sys,), torn_eqs, torn_vars) + else + solves = [] + end + + output = map(syms) do sym + if sym in required_algvars + sym + else + obs[observed_idx[sym]].rhs + end + end + + ex = Func( + [ + DestructuredArgs(diffvars) + DestructuredArgs(parameters(sys)) + independent_variable(sys) + ], + [], + Let( + collect(Iterators.flatten(solves)), + isscalar ? output[1] : MakeArray(output, output_type) + ) + ) |> Code.toexpr + + expression ? ex : @RuntimeGeneratedFunction(ex) +end + +struct ODAEProblem{iip} +end + +ODAEProblem(args...; kw...) = ODAEProblem{true}(args...; kw...) +function ODAEProblem{iip}( + sys, + u0map, + tspan, + parammap=DiffEqBase.NullParameters(); + kw... + ) where {iip} + s = structure(sys) + @unpack fullvars = s + dvs = fullvars[diffvars_range(s)] + u0map′ = ModelingToolkit.lower_mapnames(u0map, independent_variable(sys)) + u0 = ModelingToolkit.varmap_to_vars(u0map′, dvs; defaults=default_u0(sys)) + + ps = parameters(sys) + d_p = default_p(sys) + if parammap isa DiffEqBase.NullParameters && isempty(d_p) + isempty(ps) || throw(ArgumentError("The model has non-empty parameters but no parameters are specified in the problem.")) + p = parammap + else + if parammap isa DiffEqBase.NullParameters + pp = Pair[] + else + pp = ModelingToolkit.lower_mapnames(parammap) + end + p = ModelingToolkit.varmap_to_vars(pp, ps; defaults=d_p) + end + + ODEProblem{iip}(build_torn_function(sys; kw...), u0, tspan, p; kw...) +end diff --git a/src/structural_transformation/pantelides.jl b/src/structural_transformation/pantelides.jl new file mode 100644 index 0000000000..09aacdb0c7 --- /dev/null +++ b/src/structural_transformation/pantelides.jl @@ -0,0 +1,174 @@ +### +### Reassemble: structural information -> system +### + +# Naive subtree matching, we can make the dict have levels +# Going forward we should look into storing the depth in Terms +function walk_and_substitute(expr, substitution_dict) + if haskey(substitution_dict, expr) + return substitution_dict[expr] + elseif istree(expr) + return similarterm( + expr, + walk_and_substitute(operation(expr), substitution_dict), + map(x->walk_and_substitute(x, substitution_dict), arguments(expr)) + ) + else + @assert !(expr isa Equation) "RHS cannot contain equations" + return expr + end +end + +function pantelides_reassemble(sys, eqassoc, assign) + s = structure(sys) + @unpack fullvars, varassoc = s + # Step 1: write derivative equations + in_eqs = equations(sys) + out_eqs = Vector{Any}(undef, length(eqassoc)) + fill!(out_eqs, nothing) + out_eqs[1:length(in_eqs)] .= in_eqs + + out_vars = Vector{Any}(undef, length(varassoc)) + fill!(out_vars, nothing) + out_vars[1:length(fullvars)] .= fullvars + + D = Differential(independent_variable(sys)) + + for (i, v) in enumerate(varassoc) + # fullvars[v] = D(fullvars[i]) + v == 0 && continue + vi = out_vars[i] + @assert vi !== nothing "Something went wrong on reconstructing states from variable association list" + # `fullvars[i]` needs to be not a `D(...)`, because we want the DAE to be + # first-order. + if isdifferential(vi) + vi = out_vars[i] = diff2term(vi) + end + out_vars[v] = D(vi) + end + + d_dict = Dict(zip(fullvars, 1:length(fullvars))) + lhss = Set{Any}([x.lhs for x in in_eqs if isdiffeq(x)]) + for (i, e) in enumerate(eqassoc) + if e === 0 + continue + end + # LHS variable is looked up from varassoc + # the varassoc[i]-th variable is the differentiated version of var at i + eq = out_eqs[i] + lhs = if !(eq.lhs isa Symbolic) + 0 + elseif isdiffeq(eq) + # look up the variable that represents D(lhs) + lhsarg1 = arguments(eq.lhs)[1] + @assert !(lhsarg1 isa Differential) "The equation $eq is not first order" + i = get(d_dict, lhsarg1, nothing) + if i !== nothing + lhs = D(out_vars[varassoc[i]]) + if lhs in lhss + # check only trivial equations are removed + @assert isequal(diff2term(D(eq.rhs)), diff2term(lhs)) "The duplicate equation is not trivial: $eq" + lhs = Num(nothing) + end + lhs + else + D(eq.lhs) + end + else + D(eq.lhs) + end + rhs = ModelingToolkit.expand_derivatives(D(eq.rhs)) + substitution_dict = Dict(x.lhs => x.rhs for x in out_eqs if x !== nothing && x.lhs isa Symbolic) + sub_rhs = walk_and_substitute(rhs, substitution_dict) + out_eqs[e] = lhs ~ sub_rhs + end + + final_vars = unique(filter(x->!(operation(x) isa Differential), fullvars)) + final_eqs = map(identity, filter(x->value(x.lhs) !== nothing, out_eqs[sort(filter(x->x != UNASSIGNED, assign))])) + + # remove clashing equations (from order lowering vs index reduction) + return ODESystem(final_eqs, independent_variable(sys), final_vars, parameters(sys)) +end + +""" + pantelides!(sys::ODESystem; kwargs...) + +Perform Pantelides algorithm. +""" +function pantelides!(sys; maxiters = 8000) + s = structure(sys) + # D(j) = assoc[j] + @unpack graph, fullvars, varassoc = s + iv = independent_variable(sys) + neqs = nsrcs(graph) + nvars = length(varassoc) + vcolor = falses(nvars) + ecolor = falses(neqs) + assign = fill(UNASSIGNED, nvars) + eqassoc = fill(0, neqs) + neqs′ = neqs + D = Differential(iv) + for k in 1:neqs′ + eq′ = k + pathfound = false + # In practice, `maxiters=8000` should never be reached, otherwise, the + # index would be on the order of thousands. + for iii in 1:maxiters + # run matching on (dx, y) variables + # + # the derivatives and algebraic variables are zeros in the variable + # association list + varwhitelist = varassoc .== 0 + resize!(vcolor, nvars) + fill!(vcolor, false) + resize!(ecolor, neqs) + fill!(ecolor, false) + pathfound = find_augmenting_path(graph, eq′, assign, varwhitelist, vcolor, ecolor) + pathfound && break # terminating condition + for var in eachindex(vcolor); vcolor[var] || continue + # introduce a new variable + nvars += 1 + add_vertex!(graph, DST) + # the new variable is the derivative of `var` + varassoc[var] = nvars + push!(varassoc, 0) + push!(assign, UNASSIGNED) + end + + for eq in eachindex(ecolor); ecolor[eq] || continue + # introduce a new equation + neqs += 1 + add_vertex!(graph, SRC) + # the new equation is created by differentiating `eq` + eqassoc[eq] = neqs + for var in 𝑠neighbors(graph, eq) + add_edge!(graph, neqs, var) + add_edge!(graph, neqs, varassoc[var]) + end + push!(eqassoc, 0) + end + + for var in eachindex(vcolor); vcolor[var] || continue + # the newly introduced `var`s and `eq`s have the inherits + # assignment + assign[varassoc[var]] = eqassoc[assign[var]] + end + eq′ = eqassoc[eq′] + end # for _ in 1:maxiters + pathfound || error("maxiters=$maxiters reached! File a bug report if your system has a reasonable index (<100), and you are using the default `maxiters`. Try to increase the maxiters by `pantelides(sys::ODESystem; maxiters=1_000_000)` if your system has an incredibly high index and it is truly extremely large.") + end # for k in 1:neqs′ + return sys, assign, eqassoc +end + +""" + dae_index_lowering(sys::ODESystem) -> ODESystem + +Perform the Pantelides algorithm to transform a higher index DAE to an index 1 +DAE. +""" +function dae_index_lowering(sys::ODESystem; kwargs...) + s = get_structure(sys) + (s isa SystemStructure) || (sys = initialize_system_structure(sys)) + sys, assign, eqassoc = pantelides!(sys; kwargs...) + return pantelides_reassemble(sys, eqassoc, assign) +end diff --git a/src/structural_transformation/tearing.jl b/src/structural_transformation/tearing.jl new file mode 100644 index 0000000000..7ee733aa57 --- /dev/null +++ b/src/structural_transformation/tearing.jl @@ -0,0 +1,220 @@ +""" + tear_graph(sys) -> sys + +Tear the bipartite graph in a system. +""" +function tear_graph(sys) + find_solvables!(sys) + s = structure(sys) + @unpack graph, solvable_graph, assign, inv_assign, scc = s + + partitions = map(scc) do c + ieqs = filter(eq->isalgeq(s, eq), c) + vars = inv_assign[ieqs] + + td = TraverseDAG(graph.fadjlist, length(assign)) + e_solved, v_solved, e_residue, v_tear = tearEquations!(td, solvable_graph.fadjlist, ieqs, vars) + end + + @set! sys.structure.partitions = partitions + return sys +end + +function tearing_sub(expr, dict, s) + expr = ModelingToolkit.fixpoint_sub(expr, dict) + s ? simplify(expr) : expr +end + +function tearing_reassemble(sys; simplify=false) + s = structure(sys) + @unpack fullvars, partitions, assign, inv_assign, graph, scc = s + eqs = equations(sys) + + ### extract partition information + rhss = [] + solvars = [] + ns, nd = nsrcs(graph), ndsts(graph) + active_eqs = trues(ns) + active_vars = trues(nd) + rvar2req = Vector{Int}(undef, nd) + for (ith_scc, (e_solved, v_solved, e_residue, v_tear)) in enumerate(partitions) + for ii in eachindex(e_solved) + ieq = e_solved[ii]; ns -= 1 + iv = v_solved[ii]; nd -= 1 + rvar2req[iv] = ieq + + active_eqs[ieq] = false + active_vars[iv] = false + + eq = eqs[ieq] + var = fullvars[iv] + rhs = solve_for(eq, var; simplify=simplify, check=false) + # if we don't simplify the rhs and the `eq` is not solved properly + (!simplify && var in vars(rhs)) && (rhs = SymbolicUtils.polynormalize(rhs)) + # Since we know `eq` is linear wrt `var`, so the round off must be a + # linear term. We can correct the round off error by a linear + # correction. + rhs -= expand_derivatives(Differential(var)(rhs))*var + @assert !(var in vars(rhs)) """ + When solving + $eq + $var remainded in + $rhs. + """ + push!(rhss, rhs) + push!(solvars, var) + end + # DEBUG: + #@show ith_scc solvars .~ rhss + #Main._nlsys[] = eqs[e_solved], fullvars[v_solved] + #ModelingToolkit.topsort_equations(solvars .~ rhss, fullvars) + #empty!(solvars); empty!(rhss) + end + + ### update SCC + eq_reidx = Vector{Int}(undef, nsrcs(graph)) + idx = 0 + for (i, active) in enumerate(active_eqs) + eq_reidx[i] = active ? (idx += 1) : -1 + end + + rmidxs = Int[] + newscc = Vector{Int}[]; sizehint!(newscc, length(scc)) + for component′ in newscc + component = copy(component′) + for (idx, eq) in enumerate(component) + if active_eqs[eq] + component[idx] = eq_reidx[eq] + else + push!(rmidxs, idx) + end + end + push!(newscc, component) + deleteat!(component, rmidxs) + empty!(rmidxs) + end + + ### update graph + var_reidx = Vector{Int}(undef, ndsts(graph)) + idx = 0 + for (i, active) in enumerate(active_vars) + var_reidx[i] = active ? (idx += 1) : -1 + end + + newgraph = BipartiteGraph(ns, nd) + + function visit!(ii, gidx, basecase=true) + ieq = basecase ? ii : rvar2req[ii] + for ivar in 𝑠neighbors(graph, ieq) + # Note that we need to check `ii` against the rhs states to make + # sure we don't run in circles. + (!basecase && ivar === ii) && continue + if active_vars[ivar] + add_edge!(newgraph, gidx, var_reidx[ivar]) + else + # If a state is reduced, then we go to the rhs and collect + # its states. + visit!(ivar, gidx, false) + end + end + return nothing + end + + + ### update equations + newstates = setdiff([fullvars[diffvars_range(s)]; fullvars[algvars_range(s)]], solvars) + varidxmap = Dict(newstates .=> 1:length(newstates)) + neweqs = Vector{Equation}(undef, ns) + newalgeqs = falses(ns) + + dict = Dict(value.(solvars) .=> value.(rhss)) + + for ieq in Iterators.flatten(scc); active_eqs[ieq] || continue + eq = eqs[ieq] + ridx = eq_reidx[ieq] + + visit!(ieq, ridx) + + if isdiffeq(eq) + neweqs[ridx] = eq.lhs ~ tearing_sub(eq.rhs, dict, simplify) + else + newalgeqs[ridx] = true + if !(eq.lhs isa Number && eq.lhs != 0) + eq = 0 ~ eq.rhs - eq.lhs + end + rhs = tearing_sub(eq.rhs, dict, simplify) + if rhs isa Symbolic + neweqs[ridx] = 0 ~ rhs + else # a number + if abs(rhs) > 100eps(float(rhs)) + @warn "The equation $eq is not consistent. It simplifed to 0 == $rhs." + end + neweqs[ridx] = 0 ~ fullvars[inv_assign[ieq]] + end + end + end + + ### update partitions + newpartitions = similar(partitions, 0) + emptyintvec = Int[] + for ii in eachindex(partitions) + _, _, og_e_residue, og_v_tear = partitions[ii] + isempty(og_v_tear) && continue + e_residue = similar(og_e_residue) + v_tear = similar(og_v_tear) + for ii in eachindex(og_e_residue) + e_residue[ii] = eq_reidx[og_e_residue[ii]] + v_tear[ii] = var_reidx[og_v_tear[ii]] + end + # `emptyintvec` is aliased to save memory + # We need them for type stability + push!(newpartitions, (emptyintvec, emptyintvec, e_residue, v_tear)) + end + + obseqs = solvars .~ rhss + + @set! sys.structure.graph = newgraph + @set! sys.structure.scc = newscc + @set! sys.structure.fullvars = fullvars[active_vars] + @set! sys.structure.partitions = newpartitions + @set! sys.structure.algeqs = newalgeqs + @set! sys.eqs = neweqs + @set! sys.states = newstates + @set! sys.observed = vcat(observed(sys), obseqs) + return sys +end + + +""" + algebraic_equations_scc(sys) + +Find strongly connected components of algebraic equations in a system. +""" +function algebraic_equations_scc(sys) + s = get_structure(sys) + if !(s isa SystemStructure) + sys = initialize_system_structure(sys) + s = structure(sys) + end + + # skip over differential equations + algvars = isalgvar.(Ref(s), 1:ndsts(s.graph)) + eqs = equations(sys) + assign = matching(s, algvars, s.algeqs) + + components = find_scc(s.graph, assign) + inv_assign = inverse_mapping(assign) + + @set! sys.structure.assign = assign + @set! sys.structure.inv_assign = inv_assign + @set! sys.structure.scc = components + return sys +end + +""" + tearing(sys; simplify=false) + +Tear the nonlinear equations in system. When `simplify=true`, we simplify the +new residual residual equations after tearing. +""" +tearing(sys; simplify=false) = tearing_reassemble(tear_graph(algebraic_equations_scc(sys)); simplify=simplify) diff --git a/src/structural_transformation/utils.jl b/src/structural_transformation/utils.jl new file mode 100644 index 0000000000..cad871492e --- /dev/null +++ b/src/structural_transformation/utils.jl @@ -0,0 +1,259 @@ +### +### Bipartite graph utilities +### + +""" + find_augmenting_path(g::BipartiteGraph, eq, assign, varwhitelist, vcolor=falses(length(varwhitelist)), ecolor=falses(nsrcs(g))) -> path_found::Bool + +Try to find augmenting paths. +""" +function find_augmenting_path(g, eq, assign, varwhitelist, vcolor=falses(length(varwhitelist)), ecolor=falses(nsrcs(g))) + ecolor[eq] = true + + # if a `var` is unassigned and the edge `eq <=> var` exists + for var in 𝑠neighbors(g, eq) + if varwhitelist[var] && assign[var] == UNASSIGNED + assign[var] = eq + return true + end + end + + # for every `var` such that edge `eq <=> var` exists and `var` is uncolored + for var in 𝑠neighbors(g, eq) + (varwhitelist[var] && !vcolor[var]) || continue + vcolor[var] = true + if find_augmenting_path(g, assign[var], assign, varwhitelist, vcolor, ecolor) + assign[var] = eq + return true + end + end + return false +end + +""" + matching(s::Union{SystemStructure,BipartiteGraph}, varwhitelist=trues(nvars), eqwhitelist=nothing) -> assign + +Find equation-variable bipartite matching. `s.graph` is a bipartite graph. +`nvars` is the number of variables. Matched variables and equations are stored +in like `assign[var] => eq`. +""" +matching(s::SystemStructure, varwhitelist=trues(ndsts(s.graph)), eqwhitelist=nothing) = matching(s.graph, varwhitelist, eqwhitelist) +function matching(g::BipartiteGraph, varwhitelist=trues(ndsts(g)), eqwhitelist=nothing) + assign = fill(UNASSIGNED, ndsts(g)) + for eq in 𝑠vertices(g) + if eqwhitelist !== nothing + eqwhitelist[eq] || continue + end + find_augmenting_path(g, eq, assign, varwhitelist) + end + return assign +end + +### +### BLT ordering +### + +""" + find_scc(g::BipartiteGraph, assign=nothing) + +Find strongly connected components of the equations defined by `g`. `assign` +gives the undirected bipartite graph a direction. When `assign === nothing`, we +assume that the ``i``-th variable is assigned to the ``i``-th equation. +""" +function find_scc(g::BipartiteGraph, assign=nothing) + id = 0 + stack = Int[] + components = Vector{Int}[] + n = nsrcs(g) + onstack = falses(n) + lowlink = zeros(Int, n) + ids = fill(UNVISITED, n) + + for eq in 𝑠vertices(g) + if ids[eq] == UNVISITED + id = strongly_connected!(stack, onstack, components, lowlink, ids, g, assign, eq, id) + end + end + return components +end + +""" + strongly_connected!(stack, onstack, components, lowlink, ids, g, assign, eq, id) + +Use Tarjan's algorithm to find strongly connected components. +""" +function strongly_connected!(stack, onstack, components, lowlink, ids, g, assign, eq, id) + id += 1 + lowlink[eq] = ids[eq] = id + + # add `eq` to the stack + push!(stack, eq) + onstack[eq] = true + + # for `adjeq` in the adjacency list of `eq` + for var in 𝑠neighbors(g, eq) + if assign === nothing + adjeq = var + else + # assign[var] => the equation that's assigned to var + adjeq = assign[var] + # skip equations that are not assigned + adjeq == UNASSIGNED && continue + end + + # if `adjeq` is not yet idsed + if ids[adjeq] == UNVISITED # visit unvisited nodes + id = strongly_connected!(stack, onstack, components, lowlink, ids, g, assign, adjeq, id) + end + # at the callback of the DFS + if onstack[adjeq] + lowlink[eq] = min(lowlink[eq], lowlink[adjeq]) + end + end + + # if we are at a start of a strongly connected component + if lowlink[eq] == ids[eq] + component = Int[] + repeat = true + # pop until we are at the start of the strongly connected component + while repeat + w = pop!(stack) + onstack[w] = false + lowlink[w] = ids[eq] + # put `w` in current component + push!(component, w) + repeat = w != eq + end + push!(components, sort!(component)) + end + return id +end + +function sorted_incidence_matrix(sys, val=true; only_algeqs=false, only_algvars=false) + sys = algebraic_equations_scc(sys) + s = structure(sys) + @unpack assign, inv_assign, fullvars, scc, graph = s + g = graph + varsmap = zeros(Int, ndsts(graph)) + eqsmap = zeros(Int, nsrcs(graph)) + varidx = 0 + eqidx = 0 + for c in scc, eq in c + var = inv_assign[eq] + if var != 0 + eqsmap[eq] = (eqidx += 1) + varsmap[var] = (varidx += 1) + end + end + for i in diffvars_range(s) + varsmap[i] = (varidx += 1) + end + for i in dervars_range(s) + varsmap[i] = (varidx += 1) + end + for i in 1:nsrcs(graph) + if eqsmap[i] == 0 + eqsmap[i] = (eqidx += 1) + end + end + + I = Int[] + J = Int[] + for eq in 𝑠vertices(g) + only_algeqs && (isalgeq(s, eq) || continue) + for var in 𝑠neighbors(g, eq) + only_algvars && (isalgvar(s, var) || continue) + push!(I, eqsmap[eq]) + push!(J, varsmap[var]) + end + end + #sparse(I, J, val, nsrcs(g), ndsts(g)) + sparse(I, J, val) +end + +### +### Structural and symbolic utilities +### + +function find_solvables!(sys) + s = structure(sys) + @unpack fullvars, graph, solvable_graph = s + eqs = equations(sys) + empty!(solvable_graph) + for (i, eq) in enumerate(eqs) + isdiffeq(eq) && continue + term = value(eq.rhs - eq.lhs) + for j in 𝑠neighbors(graph, i) + isalgvar(s, j) || continue + D = Differential(fullvars[j]) + c = expand_derivatives(D(term), false) + if !(c isa Symbolic) && c isa Number && c != 0 + add_edge!(solvable_graph, i, j) + end + end + end + s +end + +### +### Miscellaneous +### + +function inverse_mapping(assign) + invassign = zeros(Int, length(assign)) + for (i, eq) in enumerate(assign) + eq <= 0 && continue + invassign[eq] = i + end + return invassign +end + +# debugging use +function reordered_matrix(sys, partitions=structure(sys).partitions) + s = structure(sys) + @unpack graph = s + eqs = equations(sys) + nvars = ndsts(graph) + I, J = Int[], Int[] + ii = 0 + M = Int[] + for (e_solved, v_solved, e_residue, v_tear) in partitions + append!(M, v_solved) + append!(M, v_tear) + end + M = inverse_mapping(vcat(M, setdiff(1:nvars, M))) + for (e_solved, v_solved, e_residue, v_tear) in partitions + for es in e_solved + isdiffeq(eqs[es]) && continue + ii += 1 + js = [M[x] for x in 𝑠neighbors(graph, es) if isalgvar(s, x)] + append!(I, fill(ii, length(js))) + append!(J, js) + end + + for er in e_residue + isdiffeq(eqs[er]) && continue + ii += 1 + js = [M[x] for x in 𝑠neighbors(graph, er) if isalgvar(s, x)] + append!(I, fill(ii, length(js))) + append!(J, js) + end + end + # only plot algebraic variables and equations + sparse(I, J, true) +end + +### +### Nonlinear equation(s) solving +### + +@noinline nlsolve_failure(rc) = error("The nonlinear solver failed with the return code $rc.") + +function numerical_nlsolve(f, u0, p) + prob = NonlinearProblem{false}(f, u0, p) + sol = solve(prob, NewtonRaphson()) + rc = sol.retcode + rc === :DEFAULT || nlsolve_failure(rc) + # TODO: robust initial guess, better debugging info, and residual check + sol.u +end diff --git a/test/runtests.jl b/test/runtests.jl index 957ef35867..a53e9a79fc 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -36,3 +36,4 @@ using SafeTestsets, Test println("Last test requires gcc available in the path!") @safetestset "C Compilation Test" begin include("ccompile.jl") end @safetestset "Latexify recipes Test" begin include("latexify.jl") end +@safetestset "StructuralTransformations" begin include("structural_transformation/runtests.jl") end diff --git a/test/structural_transformation/components.jl b/test/structural_transformation/components.jl new file mode 100644 index 0000000000..64cb5cadd9 --- /dev/null +++ b/test/structural_transformation/components.jl @@ -0,0 +1,104 @@ +using Test +using ModelingToolkit, OrdinaryDiffEq + +# Basic electric components +const t = Sym{ModelingToolkit.Parameter{Real}}(:t) +function Pin(name) + @variables v(t) i(t) + ODESystem(Equation[], t, [v, i], [], name=name, default_u0=[v=>1.0, i=>1.0]) +end + +function Ground(name) + g = Pin(:g) + eqs = [g.v ~ 0] + ODESystem(eqs, t, [], [], systems=[g], name=name) +end + +function ConstantVoltage(name; V = 1.0) + val = V + p = Pin(:p) + n = Pin(:n) + @parameters V + eqs = [ + V ~ p.v - n.v + 0 ~ p.i + n.i + ] + ODESystem(eqs, t, [], [V], systems=[p, n], default_p=Dict(V => val), name=name) +end + +function Resistor(name; R = 1.0) + val = R + p = Pin(:p) + n = Pin(:n) + @variables v(t) + @parameters R + eqs = [ + v ~ p.v - n.v + 0 ~ p.i + n.i + v ~ p.i * R + ] + ODESystem(eqs, t, [v], [R], systems=[p, n], default_p=Dict(R => val), name=name) +end + +function Capacitor(name; C = 1.0) + val = C + p = Pin(:p) + n = Pin(:n) + @variables v(t) + @parameters C + D = Differential(t) + eqs = [ + v ~ p.v - n.v + 0 ~ p.i + n.i + D(v) ~ p.i / C + ] + ODESystem(eqs, t, [v], [C], systems=[p, n], default_p=Dict(C => val), name=name) +end + +R = 1.0 +C = 1.0 +V = 1.0 +resistor = Resistor(:resistor, R=R) +capacitor = Capacitor(:capacitor, C=C) +source = ConstantVoltage(:source, V=V) +ground = Ground(:ground) + +function connect(ps...) + eqs = [ + 0 ~ sum(p->p.i, ps) # KCL + ] + # KVL + for i in 1:length(ps)-1 + push!(eqs, ps[i].v ~ ps[i+1].v) + end + + return eqs +end +rc_eqs = [ + connect(source.p, resistor.p) + connect(resistor.n, capacitor.p) + connect(capacitor.n, source.n, ground.g) + ] + +rc_model = ODESystem(rc_eqs, t, systems=[resistor, capacitor, source, ground], name=:rc) +sys = tearing(alias_elimination(rc_model)) +@test ModelingToolkit.default_p(sys) == Dict( + capacitor.C => 1.0, + source.V => 1.0, + resistor.R => 1.0, + ) +u0 = [ + capacitor.v => 0.0 + capacitor.p.i => 0.0 + resistor.v => 0.0 + ] +prob = ODAEProblem(sys, u0, (0, 10.0)) +sol = solve(prob, Rodas4()) + +@test sol[resistor.p.i] == sol[capacitor.p.i] +@test sol[resistor.n.i] == -sol[capacitor.p.i] +@test sol[capacitor.n.i] == -sol[capacitor.p.i] +@test iszero(sol[ground.g.i]) +@test iszero(sol[ground.g.v]) +#using Plots +#plot(sol) diff --git a/test/structural_transformation/index_reduction.jl b/test/structural_transformation/index_reduction.jl new file mode 100644 index 0000000000..99f3b7200b --- /dev/null +++ b/test/structural_transformation/index_reduction.jl @@ -0,0 +1,179 @@ +using ModelingToolkit +using LightGraphs +using DiffEqBase +using Test +using UnPack + +# Define some variables +@parameters t L g +@variables x(t) y(t) w(t) z(t) T(t) xˍt(t) yˍt(t) xˍˍt(t) yˍˍt(t) +D = Differential(t) + +eqs2 = [D(D(x)) ~ T*x, + D(D(y)) ~ T*y - g, + 0 ~ x^2 + y^2 - L^2] +pendulum2 = ODESystem(eqs2, t, [x, y, T], [L, g], name=:pendulum) +lowered_sys = ModelingToolkit.ode_order_lowering(pendulum2) + +lowered_eqs = [D(xˍt) ~ T*x, + D(yˍt) ~ T*y - g, + D(x) ~ xˍt, + D(y) ~ yˍt, + 0 ~ x^2 + y^2 - L^2,] +@test ODESystem(lowered_eqs, t, [xˍt, yˍt, x, y, T], [L, g]) == lowered_sys +@test isequal(equations(lowered_sys), lowered_eqs) + +# Simple pendulum in cartesian coordinates +eqs = [D(x) ~ w, + D(y) ~ z, + D(w) ~ T*x, + D(z) ~ T*y - g, + 0 ~ x^2 + y^2 - L^2] +pendulum = ODESystem(eqs, t, [x, y, w, z, T], [L, g], name=:pendulum) + +pendulum = initialize_system_structure(pendulum) +sss = structure(pendulum) +@unpack graph, fullvars, varassoc = sss +@test StructuralTransformations.matching(sss, varassoc .== 0) == map(x -> x == 0 ? StructuralTransformations.UNASSIGNED : x, [0, 0, 0, 0, 1, 2, 3, 4, 0]) + +sys, assign, eqassoc = StructuralTransformations.pantelides!(pendulum) +sss = structure(sys) +@unpack graph, fullvars, varassoc = sss +scc = StructuralTransformations.find_scc(graph, assign) +@test sort(sort.(scc)) == [ + [1], + [2], + [3, 4, 7, 8, 9], + [5], + [6], + ] + +@test graph.fadjlist == sort.([ + [5, 3], # 1 + [6, 4], # 2 + [7, 9, 1], # 3 + [8, 9, 2], # 4 + [2, 1], # 5 + [2, 1, 6, 5], # 6 + [5, 3, 10, 7], # 7 + [6, 4, 11, 8], # 8 + [2, 1, 6, 5, 11, 10], # 9 +]) +# [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11] +# [x, y, w, z, xˍt, yˍt, w', z', T, xˍt', yˍt'] +@test varassoc == [5, 6, 7, 8, 10, 11, 0, 0, 0, 0, 0] +#1: D(x) ~ w +#2: D(y) ~ z +#3: D(w) ~ T*x +#4: D(z) ~ T*y - g +#5: 0 ~ x^2 + y^2 - L^2 +# ---- +#6: D(eq:5) -> 0 ~ 2xx'+ 2yy' +#7: D(eq:1) -> D(D(x)) ~ D(w) -> D(xˍt) ~ D(w) -> D(xˍt) ~ T*x +#8: D(eq:2) -> D(D(y)) ~ D(z) -> D(y_t) ~ T*y - g +#9: D(eq:6) -> 0 ~ 2xx'' + 2x'x' + 2yy'' + 2y'y' +# [1, 2, 3, 4, 5, 6, 7, 8, 9] +@test eqassoc == [7, 8, 0, 0, 6, 9, 0, 0, 0] + +using ModelingToolkit +@parameters t L g +@variables x(t) y(t) w(t) z(t) T(t) xˍt(t) yˍt(t) +D = Differential(t) +idx1_pendulum = [D(x) ~ w, + D(y) ~ z, + #0 ~ x^2 + y^2 - L^2, + D(w) ~ T*x, + D(z) ~ T*y - g, + # intermediate 1: 0 ~ 2x*D(x) + 2y*D(y) - 0, + # intermediate 2(a): 0 ~ 2x*w + 2y*z - 0, (substitute D(x) and D(y)) + #0 ~ 2x*w + 2y*z, + # D(D(x)) ~ D(w) and substitute the rhs + D(xˍt) ~ T*x, + # D(D(y)) ~ D(z) and substitute the rhs + D(yˍt) ~ T*y - g, + # 2x*D(D(x)) + 2*D(x)*D(x) + 2y*D(D(y)) + 2*D(y)*D(y) and + # substitute the rhs + 0 ~ 2x*(T*x) + 2*xˍt*xˍt + 2y*(T*y - g) + 2*yˍt*yˍt] +idx1_pendulum = ODESystem(idx1_pendulum, t, [x, y, w, z, xˍt, yˍt, T], [L, g]) +first_order_idx1_pendulum = ode_order_lowering(idx1_pendulum) + +using OrdinaryDiffEq +using LinearAlgebra +prob = ODEProblem(ODEFunction(first_order_idx1_pendulum), + # [x, y, w, z, xˍt, yˍt, T] + [1, 0, 0, 0, 0, 0, 0.0],# 0, 0, 0, 0], + (0, 100.0), + [1, 9.8], + mass_matrix=calculate_massmatrix(first_order_idx1_pendulum)) +sol = solve(prob, Rodas5()); +#plot(sol, vars=(1, 2)) + +new_sys = dae_index_lowering(ModelingToolkit.ode_order_lowering(pendulum2)) + +prob_auto = ODEProblem(ODEFunction(new_sys), + # [xˍt, yˍt, x, y, T] + [0, 0, 1, 0, 0.0],# 0, 0, 0, 0], + (0, 100.0), + [1, 9.8], + mass_matrix=calculate_massmatrix(new_sys)) +sol = solve(prob_auto, Rodas5()); +#plot(sol, vars=(3, 4)) + +# Define some variables +@parameters t L g +@variables x(t) y(t) T(t) +D = Differential(t) + +eqs2 = [D(D(x)) ~ T*x, + D(D(y)) ~ T*y - g, + 0 ~ x^2 + y^2 - L^2] +pendulum2 = ODESystem(eqs2, t, [x, y, T], [L, g], name=:pendulum) + +# Turn into a first order differential equation system +first_order_sys = ModelingToolkit.ode_order_lowering(pendulum2) + +# Perform index reduction to get an Index 1 DAE +new_sys = dae_index_lowering(first_order_sys) + +u0 = [ + D(x) => 0.0, + D(y) => 0.0, + x => 1.0, + y => 0.0, + T => 0.0 +] + +p = [ + L => 1.0, + g => 9.8 +] + +prob_auto = ODEProblem(new_sys,u0,(0.0,100.0),p) +sol = solve(prob_auto, Rodas5()); +#plot(sol, vars=(D(x), y)) + +### +### More BLT/SCC tests +### + +# Test Tarjan (1972) Fig. 3 +g = [ + [2], + [3,8], + [4,7], + [5], + [3,6], + Int[], + [4,6], + [1,7], + ] +graph = StructuralTransformations.BipartiteGraph(8, 8) +for (eq, vars) in enumerate(g), var in vars + add_edge!(graph, eq, var) +end +scc = StructuralTransformations.find_scc(graph) +@test scc == [ + [6], + [3, 4, 5, 7], + [1, 2, 8], + ] diff --git a/test/structural_transformation/runtests.jl b/test/structural_transformation/runtests.jl new file mode 100644 index 0000000000..beae606172 --- /dev/null +++ b/test/structural_transformation/runtests.jl @@ -0,0 +1,6 @@ +using SafeTestsets + +@safetestset "Utilities" begin include("utils.jl") end +@safetestset "Index Reduction & SCC" begin include("index_reduction.jl") end +@safetestset "Tearing" begin include("tearing.jl") end +@safetestset "Components" begin include("components.jl") end diff --git a/test/structural_transformation/tearing.jl b/test/structural_transformation/tearing.jl new file mode 100644 index 0000000000..9b20dfc6a1 --- /dev/null +++ b/test/structural_transformation/tearing.jl @@ -0,0 +1,181 @@ +using Test +using ModelingToolkit +using ModelingToolkit: Equation +using ModelingToolkit.StructuralTransformations: SystemStructure, find_solvables! +using NLsolve +using LinearAlgebra +using UnPack + +### +### Nonlinear system +### +@parameters t +@variables u1(t) u2(t) u3(t) u4(t) u5(t) +eqs = [ + 0 ~ u1 - sin(u5), + 0 ~ u2 - cos(u1), + 0 ~ u3 - hypot(u1, u2), + 0 ~ u4 - hypot(u2, u3), + 0 ~ u5 - hypot(u4, u1), +] +sys = NonlinearSystem(eqs, [u1, u2, u3, u4, u5], []) +sys = initialize_system_structure(sys) +StructuralTransformations.find_solvables!(sys) +sss = structure(sys) +@unpack graph, solvable_graph, fullvars = sss + +io = IOBuffer() +show(io, sss) +prt = String(take!(io)) + +if VERSION >= v"1.6" +@test prt == """xvars: Any[] +dxvars: Any[] +algvars: Any[u1(t), u2(t), u3(t), u4(t), u5(t)] +Incidence matrix: + × ⋅ ⋅ ⋅ × + × × ⋅ ⋅ ⋅ + × × × ⋅ ⋅ + ⋅ × × × ⋅ + × ⋅ ⋅ × ×""" +end + +# u1 = f1(u5) +# u2 = f2(u1) +# u3 = f3(u1, u2) +# u4 = f4(u2, u3) +# u5 = f5(u4, u1) +sys = initialize_system_structure(sys) +find_solvables!(sys) +sss = structure(sys) +@unpack graph, solvable_graph, assign, partitions = sss +@test graph.fadjlist == [[1, 5], [1, 2], [1, 2, 3], [2, 3, 4], [1, 4, 5]] +@test solvable_graph.fadjlist == map(x->[x], 1:5) + +tornsys = tearing(sys) +sss = structure(tornsys) +@unpack graph, solvable_graph, assign, partitions = sss +@test graph.fadjlist == [[1]] +@test partitions == [([], [], [1], [1])] + +# Before: +# u1 u2 u3 u4 u5 +# e1 [ 1 1 ] +# e2 [ 1 1 ] +# e3 [ 1 1 1 ] +# e4 [ 1 1 1 ] +# e5 [ 1 1 1 ] +# solvable_graphs: +# u1 u2 u3 u4 u5 +# e1 [ 1 ] +# e2 [ 1 ] +# e3 [ 1 ] +# e4 [ 1 ] +# e5 [ 1 ] +# +# Optimal: +# u2 u3 u4 u5 | u1 +# e2 [ 1 | 1 ] +# e3 [ 1 1 | 1 ] +# e4 [ 1 1 1 | ] +# e5 [ 1 1 | 1 ] +# ---------------------|----- +# e1 [ 1 1 | ] +# +# Or: +# u1 u2 u3 u4 | u5 +# e1 [ 1 | 1 ] +# e2 [ 1 1 | ] +# e3 [ 1 1 1 | ] +# e4 [ 1 1 1 | ] +# --------------------|----- +# e5 [ 1 1 | 1 ] + +sys = StructuralTransformations.tear_graph(StructuralTransformations.algebraic_equations_scc(sys)) +sss = structure(sys) +@unpack partitions = sss +S = StructuralTransformations.reordered_matrix(sys, partitions) +@test S == [1 0 0 0 1 + 1 1 0 0 0 + 1 1 1 0 0 + 0 1 1 1 0 + 1 0 0 1 1] + +# unknowns: u5 +# u1 := sin(u5) +# u2 := cos(u1) +# u3 := hypot(u1, u2) +# u4 := hypot(u2, u3) +# solve for +# 0 = u5 - hypot(u1, u4) + +# unknowns: u5 +# solve for +# 0 = u5 - hypot(sin(u5), hypot(cos(sin(u5)), hypot(sin(u5), cos(sin(u5))))) +tornsys = tearing(sys) +@test isequal(equations(tornsys), [0 ~ u5 + (-1 * hypot(hypot(cos(sin(u5)), hypot(sin(u5), cos(sin(u5)))), sin(u5)))]) +nlsys_func = generate_function(tornsys)[2] +f = eval(nlsys_func) +du = zeros(1) +u = ones(1) +p = nothing +sol = nlsolve((out, x) -> f(out, x, p), u) +f(du, sol.zero, p) +@test norm(du) < 1e-10 + +### +### Simple test (edge case) +### +@parameters t +@variables x(t) y(t) z(t) +eqs = [ + 0 ~ x - y, + 0 ~ z + y, + 0 ~ x + z, + ] +nlsys = NonlinearSystem(eqs, [x, y, z], []) +newsys = tearing(nlsys) +@test equations(newsys) == [0 ~ z] +@test isequal(states(newsys), [z]) + +### +### DAE system +### +using ModelingToolkit, OrdinaryDiffEq, BenchmarkTools +@parameters t p +@variables x(t) y(t) z(t) +D = Differential(t) +eqs = [ + D(x) ~ z + 0 ~ x - y + 0 ~ sin(z) + y - p*t + ] +daesys = ODESystem(eqs, t) +newdaesys = tearing(daesys) +@test equations(newdaesys) == [D(x) ~ z; 0 ~ x + sin(z) - p*t] +@test isequal(states(newdaesys), [x, z]) +@test_throws ArgumentError ODAEProblem(newdaesys, [x=>1.0], (0, 1.0)) +prob = ODAEProblem(newdaesys, [x=>1.0], (0, 1.0), [p=>0.2]) +du = [0.0]; u = [1.0]; pr = 0.2; tt = 0.1 +@test (@ballocated $(prob.f)($du, $u, $pr, $tt)) == 0 +@test du ≈ [-asin(u[1] - pr * tt)] atol=1e-5 + +# test the initial guess is respected +infprob = ODAEProblem(tearing(ODESystem(eqs, t, default_u0=Dict(z=>Inf))), [x=>1.0], (0, 1.0), [p=>0.2]) +@test_throws DomainError infprob.f(du, u, pr, tt) + +sol1 = solve(prob, Tsit5()) +sol2 = solve(ODEProblem{false}( + (u,p,t) -> [-asin(u[1] - pr*t)], + [1.0], + (0, 1.0), + 0.2, + ), Tsit5(), tstops=sol1.t, adaptive=false) +@test Array(sol1) ≈ Array(sol2) atol=1e-5 + +obs = build_observed_function(newdaesys, [z, y]) +@test map(u -> u[2], obs.(sol1.u, pr, sol1.t)) == first.(sol1.u) +@test map(u -> sin(u[1]), obs.(sol1.u, pr, sol1.t)) + first.(sol1.u) ≈ pr[1]*sol1.t atol=1e-5 + +@test sol1[y, :] == sol1[x, :] +@test (@. sin(sol1[z, :]) + sol1[y, :]) ≈ pr * sol1.t atol=1e-5 diff --git a/test/structural_transformation/utils.jl b/test/structural_transformation/utils.jl new file mode 100644 index 0000000000..ca4083ddcd --- /dev/null +++ b/test/structural_transformation/utils.jl @@ -0,0 +1,39 @@ +using Test +using ModelingToolkit +using LightGraphs +using SparseArrays +using UnPack + +# Define some variables +@parameters t L g +@variables x(t) y(t) w(t) z(t) T(t) +D = Differential(t) + +# Simple pendulum in cartesian coordinates +eqs = [D(x) ~ w, + D(y) ~ z, + D(w) ~ T*x, + D(z) ~ T*y - g, + 0 ~ x^2 + y^2 - L^2] +pendulum = ODESystem(eqs, t, [x, y, w, z, T], [L, g], name=:pendulum) +sys = initialize_system_structure(pendulum) +StructuralTransformations.find_solvables!(sys) +sss = structure(sys) +@unpack graph, solvable_graph, fullvars, varassoc = sss +@test isequal(fullvars, [x, y, w, z, D(x), D(y), D(w), D(z), T]) +@test graph.fadjlist == sort.([[5, 3], [6, 4], [7, 1, 9], [8, 2, 9], [2, 1]]) +@test graph.badjlist == sort.([[3, 5], [4, 5], [1], [2], [1], [2], [3], [4], [3, 4]]) +@test ne(graph) == nnz(incidence_matrix(graph)) == 12 +@test nv(solvable_graph) == nv(solvable_graph) == 9 + 5 +@test varassoc == [5, 6, 7, 8, 0, 0, 0, 0, 0] + +se = collect(StructuralTransformations.𝑠edges(graph)) +@test se == mapreduce(vcat, enumerate(graph.fadjlist)) do (s, d) + StructuralTransformations.BipartiteEdge.(s, d) +end +de = collect(StructuralTransformations.𝑑edges(graph)) +@test de == mapreduce(vcat, enumerate(graph.badjlist)) do (d, s) + StructuralTransformations.BipartiteEdge.(s, d) +end +ae = collect(StructuralTransformations.edges(graph)) +@test ae == vcat(se, de)