Merge pull request #1086 from bolognam/mdb/unbalanced-error-messages

Return Unbalanced Equations/Variables

Author: ChrisRackauckas

src/structural_transformation/StructuralTransformations.jl

@@ -12,11 +12,13 @@ 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, defaults, InvalidSystemException
+using ModelingToolkit: ODESystem, AbstractSystem, var_from_nested_derivative,
+                       Differential, states, equations, vars, Symbolic,
+                       diff2term, value, operation, arguments, Sym, Term,
+                       simplify, solve_for, isdiffeq, isdifferential,
+                       get_structure, defaults, InvalidSystemException,
+                       ExtraEquationsSystemException,
+                       ExtraVariablesSystemException
 
 using ModelingToolkit.BipartiteGraphs
 using ModelingToolkit.SystemStructures

src/structural_transformation/utils.jl

@@ -47,20 +47,61 @@ function matching(g::BipartiteGraph, varwhitelist=nothing, eqwhitelist=nothing)
     return assign
 end
 
+function error_reporting(sys, bad_idxs, n_highest_vars, iseqs)
+    io = IOBuffer()
+    if iseqs
+        error_title = "More equations than variables, here are the potential extra equation(s):\n"
+        out_arr = equations(sys)[bad_idxs]
+    else
+        error_title = "More variables than equations, here are the potential extra variable(s):\n"
+        out_arr = structure(sys).fullvars[bad_idxs]
+    end
+
+    Base.print_array(io, out_arr)
+    msg = String(take!(io))
+    neqs = length(equations(sys))
+    if iseqs
+        throw(ExtraEquationsSystemException(
+            "The system is unbalanced. "
+            * "There are $n_highest_vars highest order derivative variables "
+            * "and $neqs equations.\n"
+            * error_title
+            * msg
+        ))
+    else
+        throw(ExtraVariablesSystemException(
+            "The system is unbalanced. "
+            * "There are $n_highest_vars highest order derivative variables "
+            * "and $neqs equations.\n"
+            * error_title
+            * msg
+        ))
+    end
+end
+
 ###
 ### Structural check
 ###
-function check_consistency(s::SystemStructure)
+function check_consistency(sys::AbstractSystem)
+    s = structure(sys)
     @unpack varmask, graph, varassoc, fullvars = s
     n_highest_vars = count(varmask)
     neqs = nsrcs(graph)
     is_balanced = n_highest_vars == neqs
 
-    (neqs > 0 && !is_balanced) && throw(InvalidSystemException(
-          "The system is unbalanced. "
-        * "There are $n_highest_vars highest order derivative variables "
-        * "and $neqs equations."
-       ))
+    if neqs > 0 && !is_balanced
+        varwhitelist = varassoc .== 0
+        assign = matching(graph, varwhitelist) # not assigned
+        # Just use `error_reporting` to do conditional
+        iseqs = n_highest_vars < neqs
+        if iseqs
+            inv_assign = inverse_mapping(assign) # extra equations
+            bad_idxs = findall(iszero, @view inv_assign[1:nsrcs(graph)])
+        else
+            bad_idxs = findall(isequal(UNASSIGNED), assign)
+        end
+        error_reporting(sys, bad_idxs, n_highest_vars, iseqs)
+    end
 
     # This is defined to check if Pantelides algorithm terminates. For more
     # details, check the equation (15) of the original paper.

src/systems/abstractsystem.jl

@@ -667,7 +667,7 @@ topological sort of the observed equations.
 """
 function structural_simplify(sys::AbstractSystem)
     sys = initialize_system_structure(alias_elimination(sys))
-    check_consistency(structure(sys))
+    check_consistency(sys)
     if sys isa ODESystem
         sys = dae_index_lowering(sys)
     end
@@ -688,6 +688,16 @@ struct InvalidSystemException <: Exception
 end
 Base.showerror(io::IO, e::InvalidSystemException) = print(io, "InvalidSystemException: ", e.msg)
 
+struct ExtraVariablesSystemException <: Exception
+    msg::String
+end
+Base.showerror(io::IO, e::ExtraVariablesSystemException) = print(io, "ExtraVariablesSystemException: ", e.msg)
+
+struct ExtraEquationsSystemException <: Exception
+    msg::String
+end
+Base.showerror(io::IO, e::ExtraEquationsSystemException) = print(io, "ExtraEquationsSystemException: ", e.msg)
+
 AbstractTrees.children(sys::ModelingToolkit.AbstractSystem) = ModelingToolkit.get_systems(sys)
 AbstractTrees.printnode(io::IO, sys::ModelingToolkit.AbstractSystem) = print(io, nameof(sys))
 AbstractTrees.nodetype(::ModelingToolkit.AbstractSystem) = ModelingToolkit.AbstractSystem

src/systems/systemstructure.jl

@@ -180,7 +180,10 @@ function initialize_system_structure(sys)
     for algvar in algvars
         # it could be that a variable appeared in the states, but never appeared
         # in the equations.
-        algvaridx = var2idx[algvar]
+        algvaridx = get(var2idx, algvar, 0)
+        algvaridx == 0 && throw(InvalidSystemException("The system is missing "
+            * "an equation for $algvar."
+        ))
         vartype[algvaridx] = ALGEBRAIC_VARIABLE
     end
 

test/error_handling.jl

@@ -0,0 +1,60 @@
+using Test
+using ModelingToolkit
+import ModelingToolkit: ExtraVariablesSystemException, ExtraEquationsSystemException
+
+include("../examples/electrical_components.jl")
+
+function UnderdefinedConstantVoltage(;name, V = 1.0)
+    val = V
+    @named p = Pin()
+    @named n = Pin()
+    @parameters V
+    eqs = [
+           V ~ p.v - n.v
+           # Remove equation
+           # 0 ~ p.i + n.i
+          ]
+    ODESystem(eqs, t, [], [V], systems=[p, n], defaults=Dict(V => val), name=name)
+end
+
+function OverdefinedConstantVoltage(;name, V = 1.0, I = 1.0)
+    val = V
+    val2 = I
+    @named p = Pin()
+    @named n = Pin()
+    @parameters V I
+    eqs = [
+           V ~ p.v - n.v
+           # Overdefine p.i and n.i
+           n.i ~ I
+           p.i ~ I
+          ]
+    ODESystem(eqs, t, [], [V], systems=[p, n], defaults=Dict(V => val, I => val2), name=name)
+end
+
+R = 1.0
+C = 1.0
+V = 1.0
+@named resistor = Resistor(R=R)
+@named capacitor = Capacitor(C=C)
+@named source = UnderdefinedConstantVoltage(V=V)
+
+rc_eqs = [
+          connect(source.p, resistor.p)
+          connect(resistor.n, capacitor.p)
+          connect(capacitor.n, source.n)
+         ]
+
+@named rc_model = ODESystem(rc_eqs, t, systems=[resistor, capacitor, source])
+@test_throws ModelingToolkit.ExtraVariablesSystemException structural_simplify(rc_model)
+
+
+@named source2 = OverdefinedConstantVoltage(V=V, I=V/R)
+rc_eqs2 = [
+          connect(source2.p, resistor.p)
+          connect(resistor.n, capacitor.p)
+          connect(capacitor.n, source2.n)
+         ]
+
+@named rc_model2 = ODESystem(rc_eqs2, t, systems=[resistor, capacitor, source2])
+@test_throws ModelingToolkit.ExtraEquationsSystemException structural_simplify(rc_model2)

test/reduction.jl

@@ -233,7 +233,7 @@ eqs = [
       ]
 
 @named sys = ODESystem(eqs, t, [E, C, S, P], [k₁, k₂, k₋₁, E₀])
-@test_throws ModelingToolkit.InvalidSystemException structural_simplify(sys)
+@test_throws ModelingToolkit.ExtraEquationsSystemException structural_simplify(sys)
 
 # Example 5 from Pantelides' original paper
 @parameters t

test/runtests.jl

@@ -38,3 +38,4 @@ println("Last test requires gcc available in the path!")
 @testset "Serialization" begin include("serialization.jl") end
 @safetestset "print_tree" begin include("print_tree.jl") end
 @safetestset "connectors" begin include("connectors.jl") end
+@safetestset "error_handling" begin include("error_handling.jl") end