From e34f3de0ed0edeca69688127df790c94615f0782 Mon Sep 17 00:00:00 2001 From: Keno Fischer Date: Wed, 25 Oct 2023 06:33:03 +0000 Subject: [PATCH 1/2] Give a warning for underconstrained variables forced to 0 MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit @oxinabox noticed that the following MSL test https://github.com/SciML/ModelingToolkitStandardLibrary.jl/blob/0f0fae138fe202140a4862eccaaac734051dc44f/test/Mechanical/translational.jl has a test that depends on an underconstrained variable (either `free.f` or `acc.f` may be set to an arbitrary value without affecting the dynamics of the system). Our structural singularity removal logic detects these situations and chooses one variable arbitrarily to set to 0. However, it does so silently, so users are unaware that they likely have a modeling bug. This PR adds a warning when this happens. For example, the above mentioned test now gives: ``` ┌ Warning: The model is under-constrained. Variable acc₊flange₊f(t) was arbitrarily chosen to be set to 0. This may indicate a model bug! └ @ ModelingToolkit ~/.julia/dev/ModelingToolkit/src/systems/alias_elimination.jl:54 ``` This is styled as an overridable callback, so higher level modeling frameworks that use MTK as a library can hook into this to give more domain-specific errors if desired. --- src/systems/alias_elimination.jl | 44 ++++++++++++++++++++++---------- 1 file changed, 30 insertions(+), 14 deletions(-) diff --git a/src/systems/alias_elimination.jl b/src/systems/alias_elimination.jl index 0437cb72e6..7ff4b2aaa2 100644 --- a/src/systems/alias_elimination.jl +++ b/src/systems/alias_elimination.jl @@ -1,14 +1,15 @@ using SymbolicUtils: Rewriters using Graphs.Experimental.Traversals -function alias_eliminate_graph!(state::TransformationState; kwargs...) +function alias_eliminate_graph!(state::TransformationState; + variable_underconstrained! = force_var_to_zero!, kwargs...) mm = linear_subsys_adjmat!(state; kwargs...) if size(mm, 1) == 0 return mm # No linear subsystems end @unpack graph, var_to_diff, solvable_graph = state.structure - mm = alias_eliminate_graph!(state, mm) + mm = alias_eliminate_graph!(state, mm; variable_underconstrained!) s = state.structure for g in (s.graph, s.solvable_graph) g === nothing && continue @@ -47,9 +48,17 @@ function alias_elimination!(state::TearingState; kwargs...) sys = state.sys complete!(state.structure) graph_orig = copy(state.structure.graph) - mm = alias_eliminate_graph!(state; kwargs...) - fullvars = state.fullvars + + function variable_underconstrained_mtk!(structure::SystemStructure, + ils::SparseMatrixCLIL, + v::Int) + @warn "The model is under-constrained. Variable $(fullvars[v]) was arbitrarily chosen to be set to 0. This may indicate a model bug!" + return force_var_to_zero!(structure, ils, v) + end + mm = alias_eliminate_graph!(state; + variable_underconstrained! = variable_underconstrained_mtk!, kwargs...) + @unpack var_to_diff, graph, solvable_graph = state.structure subs = Dict() @@ -348,7 +357,22 @@ function do_bareiss!(M, Mold, is_linear_variables, is_highest_diff) (rank1, rank2, rank3, pivots) end -function alias_eliminate_graph!(state::TransformationState, ils::SparseMatrixCLIL) +function force_var_to_zero!(structure::SystemStructure, ils::SparseMatrixCLIL, v::Int) + @unpack graph, solvable_graph, eq_to_diff = structure + @set! ils.nparentrows += 1 + push!(ils.nzrows, ils.nparentrows) + push!(ils.row_cols, [v]) + push!(ils.row_vals, [convert(eltype(ils), 1)]) + add_vertex!(graph, SRC) + add_vertex!(solvable_graph, SRC) + add_edge!(graph, ils.nparentrows, v) + add_edge!(solvable_graph, ils.nparentrows, v) + add_vertex!(eq_to_diff) + return ils +end + +function alias_eliminate_graph!(state::TransformationState, ils::SparseMatrixCLIL; + variable_underconstrained! = force_var_to_zero!) @unpack structure = state @unpack graph, solvable_graph, var_to_diff, eq_to_diff = state.structure # Step 1: Perform Bareiss factorization on the adjacency matrix of the linear @@ -360,15 +384,7 @@ function alias_eliminate_graph!(state::TransformationState, ils::SparseMatrixCLI rk1vars = BitSet(@view pivots[1:rank1]) for v in solvable_variables v in rk1vars && continue - @set! ils.nparentrows += 1 - push!(ils.nzrows, ils.nparentrows) - push!(ils.row_cols, [v]) - push!(ils.row_vals, [convert(eltype(ils), 1)]) - add_vertex!(graph, SRC) - add_vertex!(solvable_graph, SRC) - add_edge!(graph, ils.nparentrows, v) - add_edge!(solvable_graph, ils.nparentrows, v) - add_vertex!(eq_to_diff) + ils = variable_underconstrained!(structure, ils, v) end return ils From 09f2122c3b3b1ddbfbe22719556f02435bb85a1d Mon Sep 17 00:00:00 2001 From: Yingbo Ma Date: Tue, 26 Mar 2024 18:40:31 -0400 Subject: [PATCH 2/2] format --- src/systems/alias_elimination.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/systems/alias_elimination.jl b/src/systems/alias_elimination.jl index cc56d996fb..bdb35abd3d 100644 --- a/src/systems/alias_elimination.jl +++ b/src/systems/alias_elimination.jl @@ -2,7 +2,7 @@ using SymbolicUtils: Rewriters using Graphs.Experimental.Traversals function alias_eliminate_graph!(state::TransformationState; - variable_underconstrained! = force_var_to_zero!, kwargs...) + variable_underconstrained! = force_var_to_zero!, kwargs...) mm = linear_subsys_adjmat!(state; kwargs...) if size(mm, 1) == 0 return mm # No linear subsystems @@ -51,8 +51,8 @@ function alias_elimination!(state::TearingState; kwargs...) fullvars = state.fullvars function variable_underconstrained_mtk!(structure::SystemStructure, - ils::SparseMatrixCLIL, - v::Int) + ils::SparseMatrixCLIL, + v::Int) @warn "The model is under-constrained. Variable $(fullvars[v]) was arbitrarily chosen to be set to 0. This may indicate a model bug!" return force_var_to_zero!(structure, ils, v) end @@ -371,7 +371,7 @@ function force_var_to_zero!(structure::SystemStructure, ils::SparseMatrixCLIL, v end function alias_eliminate_graph!(state::TransformationState, ils::SparseMatrixCLIL; - variable_underconstrained! = force_var_to_zero!) + variable_underconstrained! = force_var_to_zero!) @unpack structure = state @unpack graph, solvable_graph, var_to_diff, eq_to_diff = state.structure # Step 1: Perform Bareiss factorization on the adjacency matrix of the linear