This question links to my previous post.
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.
If a solver finds a solution depends on a lot of details, e.g.
see: Common Solver Options (Solve Keyword Arguments) · DifferentialEquations.jl
Can you answer these 4 questions for the case that works and the case that doesn’t work?
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.
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?
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:
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!
OK, I get:
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…
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:
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?
This works, but not sure if that is what you want:
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
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.
Thank you for your reply!
This line does not handle tstops at all. The warning from radau5 is
┌ 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:
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
I suggest to create a bug report,
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.
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).
I think with adaptive=true, resetting dt should not be necessary (as far as convergence is concerned)
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
ContinuousCallbackis applied when a given continuous condition function hits zero.
TheDiscreteCallbackis applied when its condition function istrue…
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.
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?
Hi Chris,
Thank you for your reply! Below is the version with DiscreteCallback that fails. If you change DiscreteCallback to ContinuousCallback, it will work.
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.
What if you set dtmin=0?
With dtmin=0, the solver reports a maxiter error.
┌ 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.
There were major improvements to reinitialization of DAEs, and this seems fine now at lower tolerances.
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> 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.