Merge pull request #3931 from AayushSabharwal/as/fix-gen-initsys

fix: fix minor bug in `generate_initializesystem`

Author: AayushSabharwal

src/structural_transformation/symbolics_tearing.jl

@@ -889,6 +889,8 @@ function reorder_vars!(state::TearingState, var_eq_matching, var_sccs, eq_orderi
     # the new reality of the system we've just created.
     new_graph = contract_variables(graph, var_eq_matching, varsperm, eqsperm,
         nsolved_eq, nsolved_var)
+    new_solvable_graph = contract_variables(solvable_graph, var_eq_matching, varsperm, eqsperm,
+        nsolved_eq, nsolved_var)
 
     new_var_to_diff = complete(DiffGraph(length(var_ordering)))
     for (v, d) in enumerate(var_to_diff)
@@ -919,6 +921,7 @@ function reorder_vars!(state::TearingState, var_eq_matching, var_sccs, eq_orderi
 
     # Update system structure
     @set! state.structure.graph = complete(new_graph)
+    @set! state.structure.solvable_graph = complete(new_solvable_graph)
     @set! state.structure.var_to_diff = new_var_to_diff
     @set! state.structure.eq_to_diff = new_eq_to_diff
     @set! state.fullvars = new_fullvars

src/systems/connectors.jl

@@ -895,10 +895,12 @@ function expand_connections(sys::AbstractSystem; tol = 1e-10)
     stream_eqs, instream_subs = expand_instream(instream_csets, sys; tol = tol)
 
     eqs = [equations(sys); ceqs; stream_eqs]
-    # substitute `instream(..)` expressions with their new values
-    for i in eachindex(eqs)
-        eqs[i] = fixpoint_sub(
-            eqs[i], instream_subs; maxiters = max(length(instream_subs), 10))
+    if !isempty(instream_subs)
+        # substitute `instream(..)` expressions with their new values
+        for i in eachindex(eqs)
+            eqs[i] = fixpoint_sub(
+                eqs[i], instream_subs; maxiters = max(length(instream_subs), 10))
+        end
     end
     # get the defaults for domain networks
     d_defs = domain_defaults(sys, domain_csets)

src/systems/nonlinear/initializesystem.jl

@@ -154,7 +154,7 @@ function generate_initializesystem_timevarying(sys::AbstractSystem;
     # 3) process other variables
     for var in vars
         if var ∈ keys(op)
-            push!(eqs_ics, var ~ defs[var])
+            push!(eqs_ics, var ~ op[var])
         elseif var ∈ keys(guesses)
             push!(defs, var => guesses[var])
         elseif check_defguess
@@ -824,8 +824,6 @@ Counteracts the CSE/array variable hacks in `symbolics_tearing.jl` so it works w
 initialization.
 """
 function unhack_observed(obseqs::Vector{Equation}, eqs::Vector{Equation})
-    subs = Dict()
-    tempvars = Set()
     rm_idxs = Int[]
     for (i, eq) in enumerate(obseqs)
         iscall(eq.rhs) || continue
@@ -835,20 +833,7 @@ function unhack_observed(obseqs::Vector{Equation}, eqs::Vector{Equation})
         end
     end
 
-    for (i, eq) in enumerate(obseqs)
-        if eq.lhs in tempvars
-            subs[eq.lhs] = eq.rhs
-            push!(rm_idxs, i)
-        end
-    end
-
     obseqs = obseqs[setdiff(eachindex(obseqs), rm_idxs)]
-    obseqs = map(obseqs) do eq
-        fixpoint_sub(eq.lhs, subs) ~ fixpoint_sub(eq.rhs, subs)
-    end
-    eqs = map(eqs) do eq
-        fixpoint_sub(eq.lhs, subs) ~ fixpoint_sub(eq.rhs, subs)
-    end
     return obseqs, eqs
 end
 

src/systems/systemstructure.jl

@@ -635,6 +635,7 @@ function trivial_tearing!(ts::TearingState)
     matched_vars = BitSet()
     # variable to index in fullvars
     var_to_idx = Dict{Any, Int}(ts.fullvars .=> eachindex(ts.fullvars))
+    sys_eqs = equations(ts)
 
     complete!(ts.structure)
     var_to_diff = ts.structure.var_to_diff
@@ -654,6 +655,14 @@ function trivial_tearing!(ts::TearingState)
                 push!(blacklist, i)
                 continue
             end
+            # Edge case for `var ~ var` equations. They don't show up in the incidence
+            # graph because `TearingState` makes them `0 ~ 0`, but they do cause `var`
+            # to show up twice in `original_eqs` which fails the assertion.
+            sys_eq = sys_eqs[i]
+            if isequal(sys_eq.lhs, 0) && isequal(sys_eq.rhs, 0)
+                continue
+            end
+
             # if a variable was the LHS of two trivial observed equations, we wouldn't have
             # included it in the list. Error if somehow it made it through.
             @assert !(vari in matched_vars)

test/extensions/dynamic_optimization.jl

@@ -243,7 +243,7 @@ end
     D = M.D_nounits
 
     @parameters h_c m₀ h₀ g₀ D_c c Tₘ m_c
-    @variables h(..) v(..) m(..) [bounds = (m_c, 1)] T(..) [input = true, bounds = (0, Tₘ)]
+    @variables h(..) v(..) m(..) = m₀ [bounds = (m_c, 1)] T(..) [input = true, bounds = (0, Tₘ)]
     drag(h, v) = D_c * v^2 * exp(-h_c * (h - h₀) / h₀)
     gravity(h) = g₀ * (h₀ / h)