Use the dummy derivative to lower derivatives

Author: ChrisRackauckas

src/structural_transformation/StructuralTransformations.jl

@@ -54,7 +54,7 @@ export tearing, partial_state_selection, dae_index_lowering, check_consistency
 export dummy_derivative
 export build_torn_function, build_observed_function, ODAEProblem
 export sorted_incidence_matrix,
-       pantelides!, tearing_reassemble, find_solvables!,
+       pantelides!, pantelides_reassemble, tearing_reassemble, find_solvables!,
        linear_subsys_adjmat!
 export tearing_assignments, tearing_substitution
 export torn_system_jacobian_sparsity

src/systems/systemstructure.jl

@@ -631,6 +631,11 @@ function _structural_simplify!(state::TearingState, io; simplify = false,
     end
     if fully_determined && dummy_derivative
         sys = ModelingToolkit.dummy_derivative(sys, state; simplify, mm, check_consistency)
+    elseif fully_determined
+        var_eq_matching = pantelides!(state; finalize = false, kwargs...)
+        sys = pantelides_reassemble(state, var_eq_matching)
+        state = TearingState(sys)
+        sys = ModelingToolkit.dummy_derivative(sys, state; simplify, mm, check_consistency)
     else
         sys = ModelingToolkit.tearing(sys, state; simplify, mm, check_consistency)
     end

test/structural_transformation/index_reduction.jl

@@ -155,3 +155,20 @@ let sys = structural_simplify(pendulum2)
     @test sol.retcode == ReturnCode.Success
     @test norm(sol[x] .^ 2 + sol[y] .^ 2 .- 1) < 1e-2
 end
+
+let
+    @parameters g
+    @variables x(t) [state_priority = 10] y(t) λ(t)
+
+    eqs = [
+        D(D(x)) ~ λ * x
+        D(D(y)) ~ λ * y - g
+        x^2 + y^2 ~ 1
+        ]
+    @named pend = ODESystem(eqs,t)
+    sys = complete(structural_simplify(pend; dummy_derivative = false))
+    prob = ODEProblem(sys, [x => 1, y => 0, D(x) => 0.0], (0.0, 10.0), [g => 1], guesses = [λ => 0.0])
+    sol = solve(prob,Rodas5P())
+    @test SciMLBase.successful_retcode(sol)
+    @test sol[x^2 + y^2][end] < 1.1
+end