Handle overdetermined systems gracefully when fully_determined = false

Author: YingboMa

src/bipartite_graph.jl

@@ -422,7 +422,7 @@ return `false` may not be matched.
 """
 function maximal_matching(g::BipartiteGraph, srcfilter = vsrc -> true,
         dstfilter = vdst -> true, ::Type{U} = Unassigned) where {U}
-    matching = Matching{U}(ndsts(g))
+    matching = Matching{U}(max(nsrcs(g), ndsts(g)))
     foreach(Iterators.filter(srcfilter, 𝑠vertices(g))) do vsrc
         construct_augmenting_path!(matching, g, vsrc, dstfilter)
     end

src/parameters.jl

@@ -102,7 +102,7 @@ function split_parameters_by_type(ps)
         split_ps = tighten_types.(Base.Fix1(getindex, ps).(split_idxs))
 
         if ps isa StaticArray
-            parrs = map(x-> SArray{Tuple{size(x)...}}(x), split_ps)
+            parrs = map(x -> SArray{Tuple{size(x)...}}(x), split_ps)
             split_ps = SArray{Tuple{size(parrs)...}}(parrs)
         end
         if length(split_ps) == 1  #Tuple not needed, only 1 type

src/systems/abstractsystem.jl

@@ -213,7 +213,8 @@ function SymbolicIndexingInterface.variable_index(sys::AbstractSystem, sym::Symb
     if idx !== nothing
         return idx
     elseif count('₊', string(sym)) == 1
-        return findfirst(isequal(sym), Symbol.(sys.name, :₊, getname.(variable_symbols(sys))))
+        return findfirst(isequal(sym),
+            Symbol.(sys.name, :₊, getname.(variable_symbols(sys))))
     end
     return nothing
 end
@@ -255,7 +256,8 @@ function SymbolicIndexingInterface.parameter_index(sys::AbstractSystem, sym::Sym
     if idx !== nothing
         return idx
     elseif count('₊', string(sym)) == 1
-        return findfirst(isequal(sym), Symbol.(sys.name, :₊, getname.(parameter_symbols(sys))))
+        return findfirst(isequal(sym),
+            Symbol.(sys.name, :₊, getname.(parameter_symbols(sys))))
     end
     return nothing
 end

src/systems/diffeqs/abstractodesystem.jl

@@ -912,7 +912,10 @@ function DiffEqBase.ODEProblem(sys::AbstractODESystem, args...; kwargs...)
     ODEProblem{true}(sys, args...; kwargs...)
 end
 
-function DiffEqBase.ODEProblem(sys::AbstractODESystem, u0map::StaticArray, args...; kwargs...)
+function DiffEqBase.ODEProblem(sys::AbstractODESystem,
+        u0map::StaticArray,
+        args...;
+        kwargs...)
     ODEProblem{false, SciMLBase.FullSpecialize}(sys, u0map, args...; kwargs...)
 end
 

src/systems/nonlinear/initializesystem.jl

@@ -14,7 +14,7 @@ function initializesystem(sys::ODESystem; name = nameof(sys), guesses = Dict(),
     u0 = Vector{Pair}(undef, 0)
     defs = ModelingToolkit.defaults(sys)
 
-    full_states = [sts;getfield.((observed(sys)),:lhs)]
+    full_states = [sts; getfield.((observed(sys)), :lhs)]
 
     # Refactor to ODESystem construction
     # should be ModelingToolkit.guesses(sys)
@@ -28,26 +28,28 @@ function initializesystem(sys::ODESystem; name = nameof(sys), guesses = Dict(),
             def = defs[st]
 
             if def isa Equation
-                st ∉ keys(guesses) && error("Invalid setup: unknown $(st) has an initial condition equation with no guess.")
-                push!(eqs_ics,def)
-                push!(u0,st => guesses[st])
+                st ∉ keys(guesses) &&
+                    error("Invalid setup: unknown $(st) has an initial condition equation with no guess.")
+                push!(eqs_ics, def)
+                push!(u0, st => guesses[st])
             else
-                push!(eqs_ics,st ~ def)
+                push!(eqs_ics, st ~ def)
                 push!(u0, st => def)
             end
         elseif st ∈ keys(guesses)
-            push!(u0,st => guesses[st])
+            push!(u0, st => guesses[st])
         else
             error("Invalid setup: unknown $(st) has no default value or initial guess")
         end
     end
 
     pars = parameters(sys)
+    nleqs = [eqs_ics; observed(sys)]
 
-    sys_nl = NonlinearSystem([eqs_ics; observed(sys)],
+    sys_nl = NonlinearSystem(nleqs,
         full_states,
         pars;
-        defaults = merge(ModelingToolkit.defaults(sys),todict(u0)),
+        defaults = merge(ModelingToolkit.defaults(sys), todict(u0)),
         name,
         kwargs...)
 

test/inversemodel.jl

@@ -145,7 +145,7 @@ sol = solve(prob, Rodas5P())
 # plot(sol, idxs=[model.tank.xc, model.tank.xT, model.controller.ctr_output.u], layout=3, sp=[1 2 3])
 # hline!([prob[cm.ref.k]], label="ref", sp=1)
 
-@test sol(tspan[2], idxs = cm.tank.xc)≈ getp(prob, cm.ref.k)(prob) atol=1e-2 # Test that the inverse model led to the correct reference
+@test sol(tspan[2], idxs = cm.tank.xc)≈getp(prob, cm.ref.k)(prob) atol=1e-2 # Test that the inverse model led to the correct reference
 
 Sf, simplified_sys = Blocks.get_sensitivity_function(model, :y) # This should work without providing an operating opint containing a dummy derivative
 x, p = ModelingToolkit.get_u0_p(simplified_sys, op)

test/symbolic_indexing_interface.jl

@@ -8,11 +8,13 @@ eqs = [D(x) ~ a * y + t, D(y) ~ b * t]
 
 @test all(is_variable.((odesys,), [x, y, 1, 2, :x, :y]))
 @test all(.!is_variable.((odesys,), [a, b, t, 3, 0, :a, :b]))
-@test variable_index.((odesys,), [x, y, a, b, t, 1, 2, :x, :y, :a, :b]) == [1, 2, nothing, nothing, nothing, 1, 2, 1, 2, nothing, nothing]
+@test variable_index.((odesys,), [x, y, a, b, t, 1, 2, :x, :y, :a, :b]) ==
+      [1, 2, nothing, nothing, nothing, 1, 2, 1, 2, nothing, nothing]
 @test isequal(variable_symbols(odesys), [x, y])
 @test all(is_parameter.((odesys,), [a, b, 1, 2, :a, :b]))
 @test all(.!is_parameter.((odesys,), [x, y, t, 3, 0, :x, :y]))
-@test parameter_index.((odesys,), [x, y, a, b, t, 1, 2, :x, :y, :a, :b]) == [nothing, nothing, 1, 2, nothing, 1, 2, nothing, nothing, 1, 2]
+@test parameter_index.((odesys,), [x, y, a, b, t, 1, 2, :x, :y, :a, :b]) ==
+      [nothing, nothing, 1, 2, nothing, 1, 2, nothing, nothing, 1, 2]
 @test isequal(parameter_symbols(odesys), [a, b])
 @test all(is_independent_variable.((odesys,), [t, :t]))
 @test all(.!is_independent_variable.((odesys,), [x, y, a, :x, :y, :a]))
@@ -23,10 +25,10 @@ eqs = [D(x) ~ a * y + t, D(y) ~ b * t]
 @variables x y z
 @parameters σ ρ β
 
-eqs = [0 ~ σ*(y-x),
-       0 ~ x*(ρ-z)-y,
-       0 ~ x*y - β*z]
-@named ns = NonlinearSystem(eqs, [x,y,z],[σ,ρ,β])
+eqs = [0 ~ σ * (y - x),
+    0 ~ x * (ρ - z) - y,
+    0 ~ x * y - β * z]
+@named ns = NonlinearSystem(eqs, [x, y, z], [σ, ρ, β])
 
 @test !is_time_dependent(ns)
 
@@ -37,20 +39,20 @@ Dtt = Differential(t)^2
 Dt = Differential(t)
 
 #2D PDE
-C=1
-eq  = Dtt(u(t,x)) ~ C^2*Dxx(u(t,x))
+C = 1
+eq = Dtt(u(t, x)) ~ C^2 * Dxx(u(t, x))
 
 # Initial and boundary conditions
-bcs = [u(t,0) ~ 0.,# for all t > 0
-       u(t,1) ~ 0.,# for all t > 0
-       u(0,x) ~ x*(1. - x), #for all 0 < x < 1
-       Dt(u(0,x)) ~ 0. ] #for all  0 < x < 1]
+bcs = [u(t, 0) ~ 0.0,# for all t > 0
+    u(t, 1) ~ 0.0,# for all t > 0
+    u(0, x) ~ x * (1.0 - x), #for all 0 < x < 1
+    Dt(u(0, x)) ~ 0.0] #for all  0 < x < 1]
 
 # Space and time domains
-domains = [t ∈ (0.0,1.0),
-           x ∈ (0.0,1.0)]
+domains = [t ∈ (0.0, 1.0),
+    x ∈ (0.0, 1.0)]
 
-@named pde_system = PDESystem(eq,bcs,domains,[t,x],[u])
+@named pde_system = PDESystem(eq, bcs, domains, [t, x], [u])
 
 @test pde_system.ps == SciMLBase.NullParameters()
 @test parameter_symbols(pde_system) == []