Differential Equation solve problem

Specific Domains Modelling & Simulations
,
1

I have a dae problem to solve ,but always have singualar error, could someone help?

using DifferentialEquations, Plots

function ga!(du,u,p,t)
    t2=0.02;
    τ, kp, ki, tc=p
          
    du[1] = ki * (u[5] - u[2])     
    du[2] = 1/t2 * (u[3]-u[2])
    du[3] = u[4] / tc
    du[4] = u[4] - (u[1] + kp * (u[5] - u[2]))
    du[5] = 0.0
end

u0 = [0, 0.1, 0.1, 0.0, 0.1]

p = [0.1, 1, 1, 1]


function condition(u, t, integrator)
    t - integrator.p[1]
end

function affect!(integrator)
    integrator.u[5] = 0.9
end
cb = ContinuousCallback(condition,affect!)

discrete_condition1 = (u, t, integrator) -> integrator.u[4] > 1.0 
discrete_affect1! = integrator -> integrator.u[4] = 1.0
db1 = DiscreteCallback(discrete_condition1, discrete_affect1!)

discrete_condition2 = (u, t, integrator) -> integrator.u[4] < -1.0 
discrete_affect2! = integrator -> integrator.u[4] = -1.0
db2 = DiscreteCallback(discrete_condition2, discrete_affect2!)

M = [1.0 0.0 0.0 0.0 0.0;
     0.0 1.0 0.0 0.0 0.0;
     0.0 0.0 1.0 0.0 0.0; 
     0.0 0.0 0.0 0.0 0.0;
     0.0 0.0 0.0 0.0 0.0]

f = ODEFunction((du,u,p,t)->ga!(du,u,p,t),mass_matrix = M)

tspan = (0.0, 10.0)

prob=ODEProblem(f,u0,tspan,p;callback = CallbackSet(cb, db1, db2))

sol=solve(prob,Rosenbrock23();verbose = true)
2

That doesn’t look like an Index-1 DAE? Try using the tricks in Automated Index Reduction of DAEs · ModelingToolkit.jl

3

my system like the figure:

model
model787×317 10.7 KB

u1-u5 is odefunction state, u5 is step implement by callback, the saturation of u4 is implement by discretecallback.
as your subjust,I try the modelingtoolkitize, structural_simplify,dae_index_lowering, it dont work.

using DifferentialEquations, Plots

using ModelingToolkit
using LinearAlgebra
function ga!(du,u,p,t)
    t2=0.02;
    τ, kp, ki, tc=p
    du[1] = ki * (u[4] - u[2])     
    du[2] = 1/t2 * (u[3]-u[2])

    du[3] = u[4] / tc
    du[4] = u[4] - (u[1] + kp * (u[5] - u[2]))
    du[5] = 0.0
end

u0 = [0, 0.1, 0.1, 0.0, 0.1]
p = [0.1, 1, 1, 1]

function condition(u, t, integrator)
    t - integrator.p[1]
end

function affect!(integrator)
    integrator.u[5] = 0.9
end

cb = ContinuousCallback(condition,affect!)
discrete_condition1 = (u, t, integrator) -> integrator.u[4] > 1.0 
discrete_affect1! = integrator -> integrator.u[4] = 1.0
db1 = DiscreteCallback(discrete_condition1, discrete_affect1!)

discrete_condition2 = (u, t, integrator) -> integrator.u[4] < -1.0 
discrete_affect2! = integrator -> integrator.u[4] = -1.0
db2 = DiscreteCallback(discrete_condition2, discrete_affect2!)

f = ODEFunction((du,u,p,t)->ga!(du,u,p,t); mass_matrix = Diagonal([1,1,1,0,0]))

tspan = (0.0, 10.0)

prob=ODEProblem(f,u0,tspan,p;callback = CallbackSet(cb,db1,db2))

mtk_prob = modelingtoolkitize(prob)

sys = structural_simplify(dae_index_lowering(mtk_prob))

prob = structural_simplify(dae_index_lowering(sys))

@time sol=solve(prob,Rosenbrock23();verbose = true)

plot(sol)

so I change the code the saturation of u4 implement in code ,it work

using DifferentialEquations, Plots
using ModelingToolkit
using LinearAlgebra
function ga!(du,u,p,t)
    t2=0.02;
    τ, kp, ki, tc=p
    po =  (u[1] + kp * (u[4] - u[2]))
    po > 1 ? po = 1 : po
    po < -1 ? po = -1 : po   
    du[1] = ki * (u[4] - u[2])     
    du[2] = 1/t2 * (u[3]-u[2])
    du[3] = po / tc
end

u0 = [0, 0.1, 0.1, 0.1]
p = [0.1, 1, 1, 1]

function condition(u, t, integrator)
    t - integrator.p[1]
end

function affect!(integrator)
    integrator.u[4] = 0.9
end
cb = ContinuousCallback(condition,affect!)

f = ODEFunction((du,u,p,t)->ga!(du,u,p,t))
tspan = (0.0, 10.0)

prob=ODEProblem(f,u0,tspan,p;callback = cb)

@time sol=solve(prob,Rosenbrock23();verbose = true)

plot(sol)

so for common nonlinear part in ode or dae like Blacklash,dead zone,rate limiter,saturation,relay, which is the recommend method to implement by code it in odefunction or using callback, how about performance or recommend tutorial code for nonlinear part implement, thank @ChrisRackauckas