# DAE mass matrix formulation fails but implicit formulation works

**URL:** <https://discourse.julialang.org/t/dae-mass-matrix-formulation-fails-but-implicit-formulation-works/104940>\
**Category:** Modelling & Simulations\
**Tags:** question\
**Created:** [October 13, 2023, 10:05pm UTC](https://discourse.julialang.org/t/dae-mass-matrix-formulation-fails-but-implicit-formulation-works/104940 "2023-10-13T22:05:41Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![aaliqureshi](https://avatars.discourse-cdn.com/v4/letter/a/c5a1d2/32.png) [@aaliqureshi](https://discourse.julialang.org/u/aaliqureshi)\
**Post date:** [October 13, 2023, 10:05pm UTC](https://discourse.julialang.org/t/dae-mass-matrix-formulation-fails-but-implicit-formulation-works/104940/1 "2023-10-13T22:05:41Z")

</div>

This question links to my previous post.

> [@DiffEq: Discrete callback causes dt\<=dtmin](https://discourse.julialang.org/t/diffeq-discrete-callback-causes-dt-dtmin/104679/5):
>
> Closing this thread.

I am simulating a DAE system and noticed that the mass matrix based formulation fails but the implicit DAE formulation works fine.

Is this a bug or is there something else going on? I would appreciate if someone could provide some guidance on this. Thanks!

P.S. A MWE is given in my previous post linked above.

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [October 14, 2023, 5:09am UTC](https://discourse.julialang.org/t/dae-mass-matrix-formulation-fails-but-implicit-formulation-works/104940/2 "2023-10-14T05:09:43Z")

</div>

If a solver finds a solution depends on a lot of details, e.g.

- which solver you use
- max time step (dtmax)
- precision (abstol, reltol)
- do you provide values for the tstop argument

see: [Common Solver Options (Solve Keyword Arguments) · DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/basics/common_solver_opts/#CommonSolve.solve-Tuple%7BSciMLBase.AbstractDEProblem,%20Vararg%7BAny%7D%7D)

Can you answer these 4 questions for the case that works and the case that doesn’t work?

---

<div class="post-metadata">

**Author:** ![aaliqureshi](https://avatars.discourse-cdn.com/v4/letter/a/c5a1d2/32.png) [@aaliqureshi](https://discourse.julialang.org/u/aaliqureshi)\
**Post date:** [October 14, 2023, 4:57pm UTC](https://discourse.julialang.org/t/dae-mass-matrix-formulation-fails-but-implicit-formulation-works/104940/3 "2023-10-14T16:57:44Z")

</div>

Working Case:  
1- I used IDA solver from Sundials  
2- `dtmax = 1e-2`  
3- Precision set to default values until the `DiscreteCallback` is called. At this time, I explicitly set `abstol = 1e-2` and `reltol = 1e-2`.  
4- Yes, I have provided `tstops` vector.

Not Working Case:  
1- I have tried with different solvers including Trapezoid, Rodas4, Rodas5, Rodas5P. I also tried without specifying any particular solver and letting Diffeq choose the solver, but the same problem persists.  
2- I have tried with both adaptive as well as fixed time stepping. With adaptive time stepping, `dtmax = 1e-2`. With fixed time stepping, pre-event `dt` was set to `1e-3` and once `DiscreteCallback` is called, I use `set_proposed_dt!(integrator,integrator.opts.dtmin)` to explicitly set `dt` to its minimum value.  
3- Precision set to default values until the `DiscreteCallback` is called. At this time, I explicitly set `abstol = 1e-2` and `reltol = 1e-2`.  
4- Yes, I do provide ‘tstops’ argument. In addition to that, I also provide `d_discontinuities` argument.

Thank you very much for your time and consideration.

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [October 14, 2023, 5:01pm UTC](https://discourse.julialang.org/t/dae-mass-matrix-formulation-fails-but-implicit-formulation-works/104940/4 "2023-10-14T17:01:20Z")

</div>

> [@aaliqureshi](#):
>
> IDA solver from Sundials

Well, IDA is a very nice solver. I also used that a lot. But in theory other Julia solvers should be better. But making them work can be tricky…

I use Rodas5 now, which works very well for me. So we need to have a closer look at your not-working solution. Can you post the code again?

---

<div class="post-metadata">

**Author:** ![aaliqureshi](https://avatars.discourse-cdn.com/v4/letter/a/c5a1d2/32.png) [@aaliqureshi](https://discourse.julialang.org/u/aaliqureshi)\
**Post date:** [October 14, 2023, 5:16pm UTC](https://discourse.julialang.org/t/dae-mass-matrix-formulation-fails-but-implicit-formulation-works/104940/5 "2023-10-14T17:16:34Z")

</div>

Sure.

A little background on my problem:  
I am simulating a DAE based model of power systems. What’s very particular in power system dynamics is the sudden very large jump of algebraic variables at the time of event. The MWE is a bit large but its the simplest model I could come up with.

MWE for Not Working case:

```julia
using DifferentialEquations

function fx!(du,u,p,t)
    δ,ω,i_d,i_q,v_1,θ_1,v_3,θ_3=u
    k,e_q_prime0=p
    
    e_q_prime=e_q_prime0
    x_d_prime=0.24
    p_m=0.944
    H=3.5

    y_12=0.0
    y_13=6.66666667
    y_11=y_12+y_13
    y_21=0.0
    y_23=5.05050505
    y_22=y_21+y_23
    y_31=6.66666667
    y_32=5.05050505
    y_33=y_31+y_32

    ϕ_11=-1.57079633
    ϕ_12=0.0
    ϕ_13= 1.57079633
    ϕ_21=0.0
    ϕ_22=-1.57079633
    ϕ_23=1.57079633
    ϕ_31= 1.57079633
    ϕ_32=1.57079633
    ϕ_33=-1.57079633
    
    v_2=1.0
    θ_2=0.0
    
    Ω=2π*60
    d=1

    r_f=0
    x_f=1e-4im
    
    y_f = 1/(r_f+x_f)

    p_fault = (v_3^2)*abs(y_f)*cos(-1*angle(y_f))
    q_fault = (v_3^2)*abs(y_f)*sin(-1*angle(y_f))

    p_l3=0 - k*p_fault
    q_l3=0 - k*q_fault

    du[1] = Ω*(ω-1)
    du[2] = (p_m - (i_d*v_1*sin(δ-θ_1) + i_q*v_1*cos(δ-θ_1)) - d*(ω-1))/(2H)
    du[3] = v_1*sin(δ-θ_1) - x_d_prime*i_q
    du[4] = -e_q_prime + v_1*cos(δ-θ_1) + x_d_prime*i_d
    du[5] = i_d*v_1*sin(δ-θ_1) + i_q*v_1*cos(δ-θ_1) - (v_1*v_1*y_11*cos(θ_1 - θ_1 -ϕ_11) + v_1*v_2*y_12*cos(θ_1 - θ_2 -ϕ_12) + v_1*v_3*y_13*cos(θ_1 - θ_3 -ϕ_13))
    du[6] = i_d*v_1*cos(δ-θ_1) - i_q*v_1*sin(δ-θ_1) - (v_1*v_1*y_11*sin(θ_1 - θ_1 -ϕ_11) + v_1*v_2*y_12*sin(θ_1 - θ_2 -ϕ_12)+ v_1*v_3*y_13*sin(θ_1 - θ_3 -ϕ_13))
    du[7] = p_l3 - (v_3*v_1*y_31*cos(θ_3 - θ_1 -ϕ_31) + v_3*v_2*y_32*cos(θ_3 - θ_2 -ϕ_32) + v_3*v_3*y_33*cos(θ_3 - θ_3 -ϕ_33))
    du[8] = q_l3 - (v_3*v_1*y_31*sin(θ_3 - θ_1 -ϕ_31) + v_3*v_2*y_32*sin(θ_3 - θ_2 -ϕ_32) + v_3*v_3*y_33*sin(θ_3 - θ_3 -ϕ_33))
end

M = zeros(8,8)
M[1,1],M[2,2]=1,1

t_span=(0.0,5)

v_10=1.03
v_20=1.0
θ_20=0.0
δ0 = 0.5243271427959093
ω0 = 1.0
i_d0 = 0.41340918211007344
i_q0 = 0.8514126244761779
θ_10 = 0.32461479155703327
v_30 = 1.004014153353707
θ_30 = 0.18725718301004352
e_q_prime0 = 1.1087455985400045

u0=[δ0,ω0,i_d0,i_q0,v_10,θ_10,v_30,θ_30];

f = ODEFunction(fx!,mass_matrix=M)
k=0
p=[k,e_q_prime0]
prob = ODEProblem(f,u0,t_span,p);

const event_time1=[0.10]

function condition1(u,t,integrator,save_positions=(true,true))
    t in event_time1
end

function affect1!(integrator)
    integrator.p[1] = 1
    set_proposed_dt!(integrator,1e-15)
end

cb1=DiscreteCallback(condition1,affect1!)

const tstop = [0.10]

sol=solve(prob,Rodas5P(),callback=cb1,tstops=tstop,d_discontinuities=[0.10])

```

Thank you very much for your time!

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [October 14, 2023, 5:46pm UTC](https://discourse.julialang.org/t/dae-mass-matrix-formulation-fails-but-implicit-formulation-works/104940/6 "2023-10-14T17:46:20Z")

</div>

OK, I get:

```julia
julia> include("mwe.jl")
┌ Warning: dt(8.881784197001252e-16) <= dtmin(8.881784197001252e-16) at t=0.1, and step error estimate = 867.4048433300105. Aborting. There is either an error in your model specification or the true solution is unstable.
└ @ SciMLBase ~/.julia/packages/SciMLBase/McEqc/src/integrator_interface.jl:599
retcode: DtLessThanMin
Interpolation: specialized 4rd order "free" stiffness-aware interpolation
t: 8-element Vector{Float64}:
 0.0
 1.0e-6
 1.1e-5
 0.00011099999999999999
 0.0011109999999999998
 0.011110999999999996
 0.1
 0.1
u: 8-element Vector{Vector{Float64}}:
 [0.5243271427959093, 1.0, 0.41340918211007344, 0.8514126244761779, 1.03, 0.32461479155703327, 1.004014153353707, 0.18725718301004352]
 [0.5243271427959093, 1.0, 0.41340918221034595, 0.8514126245145126, 1.030000000002244, 0.3246147915462685, 1.0040141533558353, 0.18725718300399358]
 [0.5243271427959093, 1.0, 0.4134091822103452, 0.8514126245145128, 1.030000000002244, 0.32461479154626843, 1.0040141533558355, 0.1872571830039936]
 [0.5243271427959093, 1.0000000000000002, 0.41340918221034634, 0.8514126245145125, 1.0300000000022438, 0.32461479154626843, 1.004014153355835, 0.1872571830039936]
 [0.5243271427959101, 1.0000000000000042, 0.41340918221034667, 0.8514126245145138, 1.0300000000022438, 0.32461479154626893, 1.004014153355835, 0.18725718300399388]
 [0.524327142796002, 1.0000000000000442, 0.4134091822104251, 0.8514126245146492, 1.0300000000022318, 0.3246147915463263, 1.004014153355824, 0.1872571830040259]
 [0.5243271428028935, 1.000000000000343, 0.4134091822162912, 0.8514126245247954, 1.030000000001335, 0.3246147915506294, 1.0040141533549731, 0.18725718300644448]
 [0.5243271428028935, 1.000000000000343, 0.4134091822162912, 0.8514126245247954, 1.030000000001335, 0.3246147915506294, 1.0040141533549731, 0.18725718300644448]

```

Lets see why…

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [October 14, 2023, 6:41pm UTC](https://discourse.julialang.org/t/dae-mass-matrix-formulation-fails-but-implicit-formulation-works/104940/7 "2023-10-14T18:41:20Z")

</div>

There seams to be a principle problem. Best results so far with the Rosenbrock23() solver. I can set the integrator in the callback to 8e-5, but not higher.  
Code:

```julia
using DifferentialEquations

function fx!(du,u,p,t)
    δ,ω,i_d,i_q,v_1,θ_1,v_3,θ_3=u
    k,e_q_prime0=p
    
    e_q_prime=e_q_prime0
    x_d_prime=0.24
    p_m=0.944
    H=3.5

    y_12=0.0
    y_13=6.66666667
    y_11=y_12+y_13
    y_21=0.0
    y_23=5.05050505
    y_22=y_21+y_23
    y_31=6.66666667
    y_32=5.05050505
    y_33=y_31+y_32

    ϕ_11=-1.57079633
    ϕ_12=0.0
    ϕ_13= 1.57079633
    ϕ_21=0.0
    ϕ_22=-1.57079633
    ϕ_23=1.57079633
    ϕ_31= 1.57079633
    ϕ_32=1.57079633
    ϕ_33=-1.57079633
    
    v_2=1.0
    θ_2=0.0
    
    Ω=2π*60
    d=1

    r_f=0
    x_f=1e-4im
    
    y_f = 1/(r_f+x_f)

    p_fault = (v_3^2)*abs(y_f)*cos(-1*angle(y_f))
    q_fault = (v_3^2)*abs(y_f)*sin(-1*angle(y_f))

    p_l3=0 - k*p_fault
    q_l3=0 - k*q_fault

    du[1] = Ω*(ω-1)
    du[2] = (p_m - (i_d*v_1*sin(δ-θ_1) + i_q*v_1*cos(δ-θ_1)) - d*(ω-1))/(2H)
    du[3] = v_1*sin(δ-θ_1) - x_d_prime*i_q
    du[4] = -e_q_prime + v_1*cos(δ-θ_1) + x_d_prime*i_d
    du[5] = i_d*v_1*sin(δ-θ_1) + i_q*v_1*cos(δ-θ_1) - (v_1*v_1*y_11*cos(θ_1 - θ_1 -ϕ_11) + v_1*v_2*y_12*cos(θ_1 - θ_2 -ϕ_12) + v_1*v_3*y_13*cos(θ_1 - θ_3 -ϕ_13))
    du[6] = i_d*v_1*cos(δ-θ_1) - i_q*v_1*sin(δ-θ_1) - (v_1*v_1*y_11*sin(θ_1 - θ_1 -ϕ_11) + v_1*v_2*y_12*sin(θ_1 - θ_2 -ϕ_12)+ v_1*v_3*y_13*sin(θ_1 - θ_3 -ϕ_13))
    du[7] = p_l3 - (v_3*v_1*y_31*cos(θ_3 - θ_1 -ϕ_31) + v_3*v_2*y_32*cos(θ_3 - θ_2 -ϕ_32) + v_3*v_3*y_33*cos(θ_3 - θ_3 -ϕ_33))
    du[8] = q_l3 - (v_3*v_1*y_31*sin(θ_3 - θ_1 -ϕ_31) + v_3*v_2*y_32*sin(θ_3 - θ_2 -ϕ_32) + v_3*v_3*y_33*sin(θ_3 - θ_3 -ϕ_33))
end

M = zeros(8,8)
M[1,1], M[2,2] = 1, 1

t_span = (0.0, 0.5)

v_10 = 1.03
v_20 = 1.0
θ_20 = 0.0
δ0 = 0.5243271427959093
ω0 = 1.0
i_d0 = 0.41340918211007344
i_q0 = 0.8514126244761779
θ_10 = 0.32461479155703327
v_30 = 1.004014153353707
θ_30 = 0.18725718301004352
e_q_prime0 = 1.1087455985400045

u0 = [δ0, ω0, i_d0, i_q0, v_10, θ_10, v_30, θ_30];

f = ODEFunction(fx!,mass_matrix=M)
k = 0
p = [k,e_q_prime0]
prob = ODEProblem(f,u0,t_span,p);

const event_time1 = [0.10]

function condition1(u, t, integrator, save_positions = (true, true))
    t in event_time1
end

function affect1!(integrator)
    println("integrator: ", integrator.p[1])
    integrator.p[1] = 8e-5
    # set_proposed_dt!(integrator,1e-15)
end

cb1 = DiscreteCallback(condition1,affect1!)

const tstop = [0.10]

# abstol=tol, dtmax=1, reltol = tol,
tol = 1e-5
# Rodas4 Rodas5P Rosenbrock23
sol = solve(prob, Rosenbrock23(), callback=cb1, tstops=tstop, abstol=tol, reltol=tol, dtmax=1e-6, maxiters=1e7,
            #d_discontinuities=[0.10])
            )
sol.retcode

```

Is it possible that you need to change not only integrator.p[1], but also other states in the callback to avoid inconsistencies, and IDA does some automatic initial state solving after such a callback and the mass matrix solvers not?

Is your mass matrix constant?

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [October 14, 2023, 6:51pm UTC](https://discourse.julialang.org/t/dae-mass-matrix-formulation-fails-but-implicit-formulation-works/104940/8 "2023-10-14T18:51:13Z")

</div>

This works, but not sure if that is what you want:

```julia
using DifferentialEquations, ODEInterfaceDiffEq

function fx!(du,u,p,t)
    δ,ω,i_d,i_q,v_1,θ_1,v_3,θ_3=u
    k,e_q_prime0=p
    
    e_q_prime=e_q_prime0
    x_d_prime=0.24
    p_m=0.944
    H=3.5

    y_12=0.0
    y_13=6.66666667
    y_11=y_12+y_13
    y_21=0.0
    y_23=5.05050505
    y_22=y_21+y_23
    y_31=6.66666667
    y_32=5.05050505
    y_33=y_31+y_32

    ϕ_11=-1.57079633
    ϕ_12=0.0
    ϕ_13= 1.57079633
    ϕ_21=0.0
    ϕ_22=-1.57079633
    ϕ_23=1.57079633
    ϕ_31= 1.57079633
    ϕ_32=1.57079633
    ϕ_33=-1.57079633
    
    v_2=1.0
    θ_2=0.0
    
    Ω=2π*60
    d=1

    r_f=0
    x_f=1e-4im
    
    y_f = 1/(r_f+x_f)

    p_fault = (v_3^2)*abs(y_f)*cos(-1*angle(y_f))
    q_fault = (v_3^2)*abs(y_f)*sin(-1*angle(y_f))

    p_l3=0 - k*p_fault
    q_l3=0 - k*q_fault

    du[1] = Ω*(ω-1)
    du[2] = (p_m - (i_d*v_1*sin(δ-θ_1) + i_q*v_1*cos(δ-θ_1)) - d*(ω-1))/(2H)
    du[3] = v_1*sin(δ-θ_1) - x_d_prime*i_q
    du[4] = -e_q_prime + v_1*cos(δ-θ_1) + x_d_prime*i_d
    du[5] = i_d*v_1*sin(δ-θ_1) + i_q*v_1*cos(δ-θ_1) - (v_1*v_1*y_11*cos(θ_1 - θ_1 -ϕ_11) + v_1*v_2*y_12*cos(θ_1 - θ_2 -ϕ_12) + v_1*v_3*y_13*cos(θ_1 - θ_3 -ϕ_13))
    du[6] = i_d*v_1*cos(δ-θ_1) - i_q*v_1*sin(δ-θ_1) - (v_1*v_1*y_11*sin(θ_1 - θ_1 -ϕ_11) + v_1*v_2*y_12*sin(θ_1 - θ_2 -ϕ_12)+ v_1*v_3*y_13*sin(θ_1 - θ_3 -ϕ_13))
    du[7] = p_l3 - (v_3*v_1*y_31*cos(θ_3 - θ_1 -ϕ_31) + v_3*v_2*y_32*cos(θ_3 - θ_2 -ϕ_32) + v_3*v_3*y_33*cos(θ_3 - θ_3 -ϕ_33))
    du[8] = q_l3 - (v_3*v_1*y_31*sin(θ_3 - θ_1 -ϕ_31) + v_3*v_2*y_32*sin(θ_3 - θ_2 -ϕ_32) + v_3*v_3*y_33*sin(θ_3 - θ_3 -ϕ_33))
end

M = zeros(8,8)
M[1,1], M[2,2] = 1, 1

t_span = (0.0, 0.5)

v_10 = 1.03
v_20 = 1.0
θ_20 = 0.0
δ0 = 0.5243271427959093
ω0 = 1.0
i_d0 = 0.41340918211007344
i_q0 = 0.8514126244761779
θ_10 = 0.32461479155703327
v_30 = 1.004014153353707
θ_30 = 0.18725718301004352
e_q_prime0 = 1.1087455985400045

u0 = [δ0, ω0, i_d0, i_q0, v_10, θ_10, v_30, θ_30];

f = ODEFunction(fx!,mass_matrix=M)
k = 0
p = [k,e_q_prime0]
prob = ODEProblem(f,u0,t_span,p);

const event_time1 = [0.10]

function condition1(u, t, integrator, save_positions = (true, true))
    t in event_time1
end

function affect1!(integrator)
    println("integrator: ", integrator.p[1])
    integrator.p[1] = 1
    # set_proposed_dt!(integrator,1e-15)
end

cb1 = DiscreteCallback(condition1,affect1!)

const tstop = [0.10]

# abstol=tol, dtmax=1, reltol = tol,
tol = 1e-5
# Rodas4 Rodas5P Rosenbrock23 ROS34PW2
sol = solve(prob, radau5(), callback=cb1, abstol=tol, reltol=tol, dtmax=1e-6, maxiters=1e7)
sol.retcode

```

---

<div class="post-metadata">

**Author:** ![cuihantao](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cuihantao/32/8783_2.png) [@cuihantao](https://discourse.julialang.org/u/cuihantao)\
**Post date:** [October 14, 2023, 10:58pm UTC](https://discourse.julialang.org/t/dae-mass-matrix-formulation-fails-but-implicit-formulation-works/104940/9 "2023-10-14T22:58:46Z")

</div>

The change in `integrator.p[1]` is to connect a reactive power load to the bus. That may not require manually modifying states, as that should be the result of solving for the new condition.

It is possible that IDA has addition steps… not sure what they are though.

The mass matrix is constant. The first two are differential equations, and the rest are always algebraic equations.

---

<div class="post-metadata">

**Author:** ![cuihantao](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cuihantao/32/8783_2.png) [@cuihantao](https://discourse.julialang.org/u/cuihantao)\
**Post date:** [October 14, 2023, 11:13pm UTC](https://discourse.julialang.org/t/dae-mass-matrix-formulation-fails-but-implicit-formulation-works/104940/10 "2023-10-14T23:13:03Z")

</div>

> [@ufechner7](#):
>
> ```julia
> sol = solve(prob, radau5(), callback=cb1, abstol=tol, reltol=tol, dtmax=1e-6, maxiters=1e7)
> 
> ```

Thank you for your reply!

This line does not handle `tstops` at all. The warning from `radau5` is

```julia
┌ Warning: The tstops argument is ignored by radau5{Nothing}(nothing, nothing).
└ @ SciMLBase ~/.julia/packages/SciMLBase/emOoT/src/utils.jl:345
┌ Warning: https://docs.sciml.ai/DiffEqDocs/stable/basics/compatibility_chart/
└ @ SciMLBase ~/.julia/packages/SciMLBase/emOoT/src/utils.jl:357

```

`radau5` does not seem compatible with `tstops`, so the event is not triggered.

That can also be seen here:

```julia
t: 9-element Vector{Float64}:
 0.0
 1.0e-6
 9.0e-6
 7.3e-5
 0.000585
 0.004680999999999999
 0.037448999999999996
 0.299593
 0.5

```

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [October 15, 2023, 1:56am UTC](https://discourse.julialang.org/t/dae-mass-matrix-formulation-fails-but-implicit-formulation-works/104940/11 "2023-10-15T01:56:28Z")

</div>

I suggest to create a bug report,

---

<div class="post-metadata">

**Author:** ![aaliqureshi](https://avatars.discourse-cdn.com/v4/letter/a/c5a1d2/32.png) [@aaliqureshi](https://discourse.julialang.org/u/aaliqureshi)\
**Post date:** [October 15, 2023, 5:08am UTC](https://discourse.julialang.org/t/dae-mass-matrix-formulation-fails-but-implicit-formulation-works/104940/12 "2023-10-15T05:08:34Z")

</div>

Thank you for your time.

I would like to mention here that if I do not use `DiscreteCallback` and start the simulation with `integrator.p[1]` set to `1`, the solver works perfectly fine.

This difference in behavior when using the callback to change `integrator.p[1]` to 1, even with identical initial conditions, makes me curious.

---

<div class="post-metadata">

**Author:** ![hexaeder](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/hexaeder/32/24403_2.png) [@hexaeder](https://discourse.julialang.org/u/hexaeder)\
**Post date:** [October 15, 2023, 5:56am UTC](https://discourse.julialang.org/t/dae-mass-matrix-formulation-fails-but-implicit-formulation-works/104940/13 "2023-10-15T05:56:30Z")

</div>

For events which disturb a steadystate I found that `auto_dt_reset!(integrator)` at the end of the callback is important (since you normally have quite large timesteps before that).

---

<div class="post-metadata">

**Author:** ![aaliqureshi](https://avatars.discourse-cdn.com/v4/letter/a/c5a1d2/32.png) [@aaliqureshi](https://discourse.julialang.org/u/aaliqureshi)\
**Post date:** [October 15, 2023, 4:30pm UTC](https://discourse.julialang.org/t/dae-mass-matrix-formulation-fails-but-implicit-formulation-works/104940/14 "2023-10-15T16:30:11Z")

</div>

I think with `adaptive=true`, resetting `dt` should not be necessary (as far as convergence is concerned)

---

<div class="post-metadata">

**Author:** ![aaliqureshi](https://avatars.discourse-cdn.com/v4/letter/a/c5a1d2/32.png) [@aaliqureshi](https://discourse.julialang.org/u/aaliqureshi)\
**Post date:** [October 15, 2023, 8:17pm UTC](https://discourse.julialang.org/t/dae-mass-matrix-formulation-fails-but-implicit-formulation-works/104940/15 "2023-10-15T20:17:34Z")

</div>

Update: I was able to simulate the system successfully with using `ContinousCallback` instead of `DiscreteCallback`. I am yet to figure out why, but it seems that actions that are taken on `DiscreteCallback` might be missing something…

As per documentation,

> The `ContinuousCallback` is applied when a given continuous _condition function_ hits zero.  
> The `DiscreteCallback` is applied when its _condition function_ is `true` …

This justifies my earlier use of `DiscreteCallback`.

I would appreciate if someone from the developer community could comment on this observation.

Thank you all for your kind responses.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [October 15, 2023, 8:46pm UTC](https://discourse.julialang.org/t/dae-mass-matrix-formulation-fails-but-implicit-formulation-works/104940/16 "2023-10-15T20:46:33Z")

</div>

What happens if you set `u_modified!(integrator,true)`?

Can you share the DiscreteCallback and ContinuousCallback versions so it’s straightforward to compare them side by side?

---

<div class="post-metadata">

**Author:** ![cuihantao](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cuihantao/32/8783_2.png) [@cuihantao](https://discourse.julialang.org/u/cuihantao)\
**Post date:** [October 15, 2023, 9:27pm UTC](https://discourse.julialang.org/t/dae-mass-matrix-formulation-fails-but-implicit-formulation-works/104940/17 "2023-10-15T21:27:12Z")

</div>

> [@ChrisRackauckas](#):
>
> Can you share the DiscreteCallback and ContinuousCallback versions so it’s straightforward to compare them side by side?

Hi Chris,

Thank you for your reply! Below is the version with `DiscreteCallback` that fails. If you change `DiscreteCallback` to `ContinuousCallback`, it will work.

```julia
using OrdinaryDiffEq

function fx!(du,u,p,t)
    δ,ω,i_d,i_q,v_1,θ_1,v_3,θ_3=u
    k,e_q_prime0=p
    
    e_q_prime=e_q_prime0
    x_d_prime=0.24
    p_m=0.944
    H=3.5

    y_12=0.0
    y_13=6.66666667
    y_11=y_12+y_13
    y_21=0.0
    y_23=5.05050505
    y_22=y_21+y_23
    y_31=6.66666667
    y_32=5.05050505
    y_33=y_31+y_32

    ϕ_11=-1.57079633
    ϕ_12=0.0
    ϕ_13= 1.57079633
    ϕ_21=0.0
    ϕ_22=-1.57079633
    ϕ_23=1.57079633
    ϕ_31= 1.57079633
    ϕ_32=1.57079633
    ϕ_33=-1.57079633
    
    v_2=1.0
    θ_2=0.0
    
    Ω=2π*60
    d=1

    r_f=0
    x_f=1e-4im
    
    y_f = 1/(r_f+x_f)

    p_fault = (v_3^2)*abs(y_f)*cos(-1*angle(y_f))
    q_fault = (v_3^2)*abs(y_f)*sin(-1*angle(y_f))

    p_l3=0 - k*p_fault
    q_l3=0 - k*q_fault

    du[1] = Ω*(ω-1)
    du[2] = (p_m - (i_d*v_1*sin(δ-θ_1) + i_q*v_1*cos(δ-θ_1)) - d*(ω-1))/(2H)
    du[3] = v_1*sin(δ-θ_1) - x_d_prime*i_q
    du[4] = -e_q_prime + v_1*cos(δ-θ_1) + x_d_prime*i_d
    du[5] = i_d*v_1*sin(δ-θ_1) + i_q*v_1*cos(δ-θ_1) - (v_1*v_1*y_11*cos(θ_1 - θ_1 -ϕ_11) + v_1*v_2*y_12*cos(θ_1 - θ_2 -ϕ_12) + v_1*v_3*y_13*cos(θ_1 - θ_3 -ϕ_13))
    du[6] = i_d*v_1*cos(δ-θ_1) - i_q*v_1*sin(δ-θ_1) - (v_1*v_1*y_11*sin(θ_1 - θ_1 -ϕ_11) + v_1*v_2*y_12*sin(θ_1 - θ_2 -ϕ_12)+ v_1*v_3*y_13*sin(θ_1 - θ_3 -ϕ_13))
    du[7] = p_l3 - (v_3*v_1*y_31*cos(θ_3 - θ_1 -ϕ_31) + v_3*v_2*y_32*cos(θ_3 - θ_2 -ϕ_32) + v_3*v_3*y_33*cos(θ_3 - θ_3 -ϕ_33))
    du[8] = q_l3 - (v_3*v_1*y_31*sin(θ_3 - θ_1 -ϕ_31) + v_3*v_2*y_32*sin(θ_3 - θ_2 -ϕ_32) + v_3*v_3*y_33*sin(θ_3 - θ_3 -ϕ_33))
end

M = zeros(8,8)
M[1,1],M[2,2]=1,1

t_span=(0.0,5)

v_10=1.03
v_20=1.0
θ_20=0.0
δ0 = 0.5243271427959093
ω0 = 1.0
i_d0 = 0.41340918211007344
i_q0 = 0.8514126244761779
θ_10 = 0.32461479155703327
v_30 = 1.004014153353707
θ_30 = 0.18725718301004352
e_q_prime0 = 1.1087455985400045

u0=[δ0,ω0,i_d0,i_q0,v_10,θ_10,v_30,θ_30];

f = ODEFunction(fx!,mass_matrix=M)
k=0
p=[k,e_q_prime0]
prob = ODEProblem(f,u0,t_span,p);

function condition1(u,t,integrator,save_positions=(true,true))
    t == 0.1
end

function affect1!(integrator)
    @show "entering affect1"
    integrator.p[1] = 1
    u_modified!(integrator,true)
end

cb1=DiscreteCallback(condition1,affect1!)

sol=solve(prob,Rodas5P(),callback=cb1,tstops=[0.1], dtmax=0.02)

```

Edit: removed an undef var.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [October 15, 2023, 9:29pm UTC](https://discourse.julialang.org/t/dae-mass-matrix-formulation-fails-but-implicit-formulation-works/104940/18 "2023-10-15T21:29:14Z")

</div>

What if you set `dtmin=0`?

---

<div class="post-metadata">

**Author:** ![cuihantao](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cuihantao/32/8783_2.png) [@cuihantao](https://discourse.julialang.org/u/cuihantao)\
**Post date:** [October 15, 2023, 9:32pm UTC](https://discourse.julialang.org/t/dae-mass-matrix-formulation-fails-but-implicit-formulation-works/104940/19 "2023-10-15T21:32:22Z")

</div>

With `dtmin=0`, the solver reports a maxiter error.

```julia
┌ Warning: Interrupted. Larger maxiters is needed. If you are using an integrator for non-stiff ODEs or an automatic switching algorithm (the default), you may want to consider using a method for stiff equations. See the solver pages for more details (e.g. https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/#Stiff-Problems).
└ @ SciMLBase ~/.julia/packages/SciMLBase/McEqc/src/integrator_interface.jl:580
retcode: MaxIters
Interpolation: specialized 4rd order "free" stiffness-aware interpolation

```

@aaliqureshi also observed a chattering issue of an algebraic variable when using Trapezoid method. That chattering causes algebraic variables to jump between two values, which does not look like the numerical oscillation of Trapezoid method. I’m not sure though if that’s related to this reaching maxiter.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [April 7, 2024, 3:32pm UTC](https://discourse.julialang.org/t/dae-mass-matrix-formulation-fails-but-implicit-formulation-works/104940/20 "2024-04-07T15:32:27Z")

</div>

There were major improvements to reinitialization of DAEs, and this seems fine now at lower tolerances.

```julia
using OrdinaryDiffEq

function fx!(du,u,p,t)
    δ,ω,i_d,i_q,v_1,θ_1,v_3,θ_3=u
    k,e_q_prime0=p

    e_q_prime=e_q_prime0
    x_d_prime=0.24
    p_m=0.944
    H=3.5

    y_12=0.0
    y_13=6.66666667
    y_11=y_12+y_13
    y_21=0.0
    y_23=5.05050505
    y_22=y_21+y_23
    y_31=6.66666667
    y_32=5.05050505
    y_33=y_31+y_32

    ϕ_11=-1.57079633
    ϕ_12=0.0
    ϕ_13= 1.57079633
    ϕ_21=0.0
    ϕ_22=-1.57079633
    ϕ_23=1.57079633
    ϕ_31= 1.57079633
    ϕ_32=1.57079633
    ϕ_33=-1.57079633

    v_2=1.0
    θ_2=0.0

    Ω=2π*60
    d=1

    r_f=0
    x_f=1e-4im

    y_f = 1/(r_f+x_f)

    p_fault = (v_3^2)*abs(y_f)*cos(-1*angle(y_f))
    q_fault = (v_3^2)*abs(y_f)*sin(-1*angle(y_f))

    p_l3=0 - k*p_fault
    q_l3=0 - k*q_fault

    du[1] = Ω*(ω-1)
    du[2] = (p_m - (i_d*v_1*sin(δ-θ_1) + i_q*v_1*cos(δ-θ_1)) - d*(ω-1))/(2H)
    du[3] = v_1*sin(δ-θ_1) - x_d_prime*i_q
    du[4] = -e_q_prime + v_1*cos(δ-θ_1) + x_d_prime*i_d
    du[5] = i_d*v_1*sin(δ-θ_1) + i_q*v_1*cos(δ-θ_1) - (v_1*v_1*y_11*cos(θ_1 - θ_1 -ϕ_11) + v_1*v_2*y_12*cos(θ_1 - θ_2 -ϕ_12) + v_1*v_3*y_13*cos(θ_1 - θ_3 -ϕ_13))
    du[6] = i_d*v_1*cos(δ-θ_1) - i_q*v_1*sin(δ-θ_1) - (v_1*v_1*y_11*sin(θ_1 - θ_1 -ϕ_11) + v_1*v_2*y_12*sin(θ_1 - θ_2 -ϕ_12)+ v_1*v_3*y_13*sin(θ_1 - θ_3 -ϕ_13))
    du[7] = p_l3 - (v_3*v_1*y_31*cos(θ_3 - θ_1 -ϕ_31) + v_3*v_2*y_32*cos(θ_3 - θ_2 -ϕ_32) + v_3*v_3*y_33*cos(θ_3 - θ_3 -ϕ_33))
    du[8] = q_l3 - (v_3*v_1*y_31*sin(θ_3 - θ_1 -ϕ_31) + v_3*v_2*y_32*sin(θ_3 - θ_2 -ϕ_32) + v_3*v_3*y_33*sin(θ_3 - θ_3 -ϕ_33))
end

M = zeros(8,8)
M[1,1],M[2,2]=1,1

t_span=(0.0,5)

v_10=1.03
v_20=1.0
θ_20=0.0
δ0 = 0.5243271427959093
ω0 = 1.0
i_d0 = 0.41340918211007344
i_q0 = 0.8514126244761779
θ_10 = 0.32461479155703327
v_30 = 1.004014153353707
θ_30 = 0.18725718301004352
e_q_prime0 = 1.1087455985400045

u0=[δ0,ω0,i_d0,i_q0,v_10,θ_10,v_30,θ_30];

f = ODEFunction(fx!,mass_matrix=M)
k=0
p=[k,e_q_prime0]
prob = ODEProblem(f,u0,t_span,p);

function condition1(u,t,integrator,save_positions=(true,true))
    t == 0.1
end

function affect1!(integrator)
    @show "entering affect1"
    integrator.p[1] = 1
    u_modified!(integrator,true)
end

cb1=DiscreteCallback(condition1,affect1!)

p=[k,e_q_prime0]
prob = ODEProblem(f,u0,t_span,p);
sol1=solve(prob,Rodas5P(),callback=cb1,tstops=[0.1], saveat = 0.1, abstol=1e-10, reltol=1e-10)

p=[k,e_q_prime0]
prob = ODEProblem(f,u0,t_span,p);
sol2=solve(prob,FBDF(),callback=cb1,tstops=[0.1], saveat = 0.1, abstol=1e-10, reltol=1e-10)

```

Both solve, checking the difference:

```julia
julia> sol1 - sol2
VectorOfArray{Float64,2}:
52-element Vector{Vector{Float64}}:
 [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
 [-9.848788451449764e-13, 1.2589929099249275e-13, -8.392730954653871e-13, -1.4493961586481419e-12, 1.283417816466681e-13, -6.151745779447992e-13, 1.2168044349891716e-13, -3.456956942926581e-13]
 [-9.848788451449764e-13, 1.2589929099249275e-13, -1.787725523172412e-11, 3.288296284570391e-12, 4.9669712787192566e-12, 1.4911405443740477e-12, 8.300534305033525e-12, -4.545008702727671e-9]
 [-2.3977719809664677e-10, 1.1612932837579137e-13, -7.131387880221496e-8, -2.4595117536249482e-8, 1.7109302419271444e-8, -5.209714711240654e-5, 1.1407411998299796e-8, -1.6107900885575788e-5]
 [-2.400699639082404e-10, 8.548717289613705e-14, -1.3591581637939498e-8, 2.545654042445701e-8, 3.2664175253493966e-9, -9.81995196513985e-6, 4.586860387497359e-9, 1.7506057102822936e-6]
 [-2.3677015903444953e-10, 4.773959005888173e-14, 5.39947775379801e-9, 3.534371990226732e-9, -1.296616924939542e-9, -1.068550248284339e-6, -1.5788576581115404e-9, 1.0352748133934941e-5]
 [-2.198099480210658e-10, 2.3092638912203256e-14, 5.8237463740340445e-9, -1.7031891305294145e-8, -1.3949167376736682e-9, -1.6947831351643572e-6, 6.419038963562107e-10, -4.72278221397282e-6]
 ⋮
 [6.565073817910161e-9, -4.0678571622265736e-13, 1.8942347246309055e-8, -2.7873452780689117e-8, -4.543984410965152e-9, -4.644618911697762e-9, 7.285772300139742e-12, -1.7673110691007743e-5]
 [6.442348876589676e-9, 2.156053113822054e-13, 6.684017339608772e-9, 4.116655970623495e-8, -1.5972529410213099e-9, -4.8331173729820875e-8, 2.77156804451284e-9, 9.480809935347168e-6]
 [6.3961351770558394e-9, -1.3427037259816643e-12, -4.3778739300393e-8, -1.6818814049313676e-8, 1.0500807001712076e-8, -1.8222806374978973e-8, 6.449213563902191e-9, -8.993695660708401e-7]
 [6.219920578587335e-9, -1.261879489788953e-12, 7.422707337667589e-9, -9.85985547526623e-9, -1.7808509111816306e-9, 3.201819254172733e-9, -1.3493299235107017e-10, -7.1553641446397975e-6]
 [6.038987976353383e-9, -1.4672707493446069e-12, 4.871099079650776e-9, 3.7149643253124587e-9, -1.169345953400125e-9, -1.0965095498249866e-10, -6.293475783916702e-10, 6.202034353530728e-6]
 [3.4690970096562523e-9, -6.326250634458574e-11, -9.49850464948554e-9, 2.320217818438733e-8, 2.2827632562005817e-9, -2.543339405747247e-8, 2.7335558761215975e-9, 2.267419438817253e-6]
 [1.0565486263658386e-9, -6.265277185946161e-11, 1.5987211554602254e-13, 1.3588889941681254e-12, -3.824718319833664e-14, 1.055809661920648e-9, -2.2072404667894396e-14, 3.808920956416273e-10]

```

Both give the same solution.

I would recommend FBDF for this kind of case.