diff --git a/lib/ModelingToolkitBase/src/problems/nonlinearproblem.jl b/lib/ModelingToolkitBase/src/problems/nonlinearproblem.jl index 1728ed73f3..1a5dc153ca 100644 --- a/lib/ModelingToolkitBase/src/problems/nonlinearproblem.jl +++ b/lib/ModelingToolkitBase/src/problems/nonlinearproblem.jl @@ -90,8 +90,8 @@ function generate_nonlinear_bounds(sys::AbstractSystem, op) left_merge!(lbmap, op) lb = varmap_to_vars(lbmap, dvs; tofloat = false) ubmap = SymmapT() - for (var, b) in zip(dvs, lb) - b === -Inf && continue + for (var, b) in zip(dvs, ub) + b === Inf && continue write_possibly_indexed_array!(ubmap, var, Symbolics.SConst(b), COMMON_NOTHING) end left_merge!(ubmap, op) diff --git a/lib/ModelingToolkitBase/src/systems/index_cache.jl b/lib/ModelingToolkitBase/src/systems/index_cache.jl index c03336f2c6..ed0dcc48a6 100644 --- a/lib/ModelingToolkitBase/src/systems/index_cache.jl +++ b/lib/ModelingToolkitBase/src/systems/index_cache.jl @@ -644,14 +644,20 @@ function reorder_parameters(ic::IndexCache, ps::Vector{SymbolicT}; drop_missing if i isa Int param_buf[i] = p else - param_buf[i] = collect(p) + i = (first(i)::Int):(last(i)::Int) + for (buf_i, p_i) in zip(i, SU.stable_eachindex(p)) + param_buf[buf_i] = p[p_i] + end end elseif haskey(ic.initials_idx, p) i = ic.initials_idx[p] if i isa Int initials_buf[i] = p else - initials_buf[i] = collect(p) + i = (first(i)::Int):(last(i)::Int) + for (buf_i, p_i) in zip(i, SU.stable_eachindex(p)) + initials_buf[buf_i] = p[p_i] + end end elseif haskey(ic.constant_idx, p) i, j = ic.constant_idx[p] diff --git a/lib/ModelingToolkitBase/src/systems/system.jl b/lib/ModelingToolkitBase/src/systems/system.jl index 13b15fd6ea..29eb44c947 100644 --- a/lib/ModelingToolkitBase/src/systems/system.jl +++ b/lib/ModelingToolkitBase/src/systems/system.jl @@ -1176,7 +1176,7 @@ end Get the metadata associated with key `k` in system `sys` or `default` if it does not exist. """ -function SymbolicUtils.getmetadata(sys::AbstractSystem, k::DataType, default) +function SymbolicUtils.getmetadata(sys::AbstractSystem, @nospecialize(k::DataType), default) meta = get_metadata(sys) return get(meta, k, default) end @@ -1188,7 +1188,7 @@ Set the metadata associated with key `k` in system `sys` to value `v`. This is a out-of-place operation, and will return a shallow copy of `sys` with the appropriate metadata values. """ -function SymbolicUtils.setmetadata(sys::AbstractSystem, k::DataType, v) +function SymbolicUtils.setmetadata(sys::AbstractSystem, @nospecialize(k::DataType), @nospecialize(v)) meta = get_metadata(sys) meta = Base.ImmutableDict(meta, k => v)::MetadataT return @set sys.metadata = meta diff --git a/lib/ModelingToolkitBase/src/systems/systems.jl b/lib/ModelingToolkitBase/src/systems/systems.jl index 579bfdec9b..51ce5a689d 100644 --- a/lib/ModelingToolkitBase/src/systems/systems.jl +++ b/lib/ModelingToolkitBase/src/systems/systems.jl @@ -983,6 +983,12 @@ function observed2graph(sys::AbstractSystem, eqs::Vector{Equation}, unknowns::Ve return graph, assigns end +""" +Toggle to control whether `topsort_equations` prints the equations in the +cycle, if present. +""" +TOPSORT_EQS_PRINT_CYCLE::Bool = false + """ $(TYPEDSIGNATURES) @@ -1054,7 +1060,37 @@ function topsort_equations(sys::AbstractSystem, eqs::Vector{Equation}, unknowns: end end - (check && idx != neqs) && throw(ArgumentError("The equations have at least one cycle.")) + if check && idx != neqs + # Build a directed eq→eq subgraph over unsorted equations, find smallest SCC. + if TOPSORT_EQS_PRINT_CYCLE + unsorted = findall(>(0), degrees) + unsorted_set = Set(unsorted) + n_unsorted = length(unsorted) + old_to_new = Dict(old => new for (new, old) in enumerate(unsorted)) + + g = SimpleDiGraph(n_unsorted) + for src_old in unsorted + for dst_old in 𝑑neighbors(graph, assigns[src_old]) + dst_old in unsorted_set || continue + add_edge!(g, old_to_new[src_old], old_to_new[dst_old]) + end + end + + sccs = strongly_connected_components(g) + nontrivial = filter(scc -> length(scc) >= 2, sccs) + smallest_new = isempty(nontrivial) ? collect(1:n_unsorted) : + nontrivial[argmin(length.(nontrivial))] + + println("=== topsort_equations: CYCLE DETECTED ===") + println("Smallest cycle ($(length(smallest_new)) equations):") + for new_idx in smallest_new + old_idx = unsorted[new_idx] + println(" LHS = $(unknowns[assigns[old_idx]])") + println(" EQ = $(eqs[old_idx])") + end + end + throw(ArgumentError("The equations have at least one cycle.")) + end return ordered_eqs end diff --git a/lib/ModelingToolkitBase/test/initializationsystem.jl b/lib/ModelingToolkitBase/test/initializationsystem.jl index c927fd4c1c..a800aa845c 100644 --- a/lib/ModelingToolkitBase/test/initializationsystem.jl +++ b/lib/ModelingToolkitBase/test/initializationsystem.jl @@ -326,8 +326,15 @@ end initprob = ModelingToolkitBase.InitializationProblem(sys, 0.0; missing_guess_value) conditions = getfield.(equations(initprob.f.sys), :rhs) -@test initprob isa SCCNonlinearProblem -@test initprob.probs isa Tuple{<:LinearProblem} +# The init system as-is is linear. Since it is so severely overdetermined, MTK can sometimes +# choose a bad simplification that results in a nonlinear system. +# TODO: Support retaining `observed` equations in proper `mtkcompile` and generate +# `x ~ Initial(x)` as such `observed` equations. This makes `mtkcompile`'s job easier +# on initialization systems, but can lead to parameter-only equations. +if !@isdefined(ModelingToolkit) + @test initprob isa SCCNonlinearProblem + @test initprob.probs isa Tuple{<:LinearProblem} +end if @isdefined(ModelingToolkit) initsol = solve(initprob, reltol = 1.0e-12, abstol = 1.0e-12) @test SciMLBase.successful_retcode(initsol)