Switch in equations using ContinuousCallback is handled by Rodas5() but not QNDF()

Specific Domains Modelling & Simulations
,
1

Here I have a problem where I’d like a constraint equation to switch once a constraint has been met for one of the auxiliary variables. This is a buildup to a
large, stiff system that benefits greatly in performance from using multistep methods. The code works for the Rodas5 solver, but not for QNDF.

The equation relating the two auxiliary variables is: x_3 = x_1 + x_4
The initial constraint equation is: x_4 = 0.5
Once the value of x_3 reaches 50.0, I would like to switch the constraint equation to x_3 = 50.0, as done in the code.

The value of x_3 for QNDF at the condition appears to reach ~53.

Please let me know if I have constructed something erroneously in the callback or how I may be able to contribute to understanding this and fixing what is happening.

Thank you

using ModelingToolkit
import ModelingToolkit: t_nounits as t, D_nounits as D
using Parameters
using OrdinaryDiffEq

function build_system(; name = :sys)
    states = @variables begin
        x_1(t) = 1.0
        x_2(t) = 1.0
        x_3(t) = 1.5, [irreducible=true]
        x_4(t) = 0.5, [irreducible=true]
    end

    params = @parameters begin
        α_1 = 0.5
        α_2 = 0.5
        α_4 = 0.5
        mode::Int = 1
    end

    equations = [
        D(x_1) ~ α_1 * x_1 + α_2 * x_2,
        D(x_2) ~ α_2 * x_2,
        x_3 ~ x_1 + x_4,
        0 ~ ((mode - 1) == 0) * (x_4 - α_4) + ((mode - 2) == 0) * (x_3 - 50.0)
        ]

    return ODESystem(equations, t, [states...], params, name=name)
end

sys = build_system()
simp = structural_simplify(sys)
prob_rodas = ODEProblem(simp, ModelingToolkit.get_defaults(simp), (0, 7.0))

function cb_condition(u,t,integ)
    return u[ModelingToolkit.variable_index(integ.f.sys, :x_3)] - 50.0
end

function cb_affect!(integ)
    integ.ps[:mode] = 2
end

cb = ContinuousCallback(cb_condition, cb_affect!, save_positions=(true, true), interp_points=0)
# solves and the switch in equations is handled expectedly
sol_rodas = solve(prob_rodas, Rodas5(), callback = cb, initializealg=CheckInit())

# set up new problem since mode was changed in previous integrator from callback
prob_qndf = ODEProblem(simp, ModelingToolkit.get_defaults(simp), (0, 7.0))
# reinitialization after affect has been applied is not consistent because x_3 does not equal 50.0
sol_qndf = solve(prob_qndf, QNDF(), callback = cb, initializealg=CheckInit())
2

Would it be because of the interpolant?

3

What happens when you reduce the tolerance?

QNDF does need a better interpolation, that is true.

4

Hello Chris!

With abs/reltol sent up to 1e-2 I am still getting:

ERROR: DAE initialization failed: your u0 did not satisfy the initialization requirements, normresid = 0.47353044373979003 > abstol = 0.01.
If you wish for the system to automatically change the algebraic variables to satisfy the algebraic constraints, please pass `initializealg = BrownBasicInit()` to solve (this option will require `using OrdinaryDiffEqNonlinearSolve`). If you wish to perform an initialization on the complete u0, please pass `initializealg = ShampineCollocationInit()` to `solve`.

I set to CheckInit() because I constructed the problem in such a way that it should be consistent at the event and initial conditions.

The printout of u, which is [x_1, x_2, x_4, x_3] because of reordering, within the callback condition is:

u = [1.0, 1.0, 0.5, 1.5]
u = [1.0000008438823904, 1.000000421941106, 0.5, 1.5000008438823904]
u = [1.0000008438823904, 1.000000421941106, 0.5, 1.5000008438823904]
u = [1.0000106002215667, 1.0000053000993734, 0.5, 1.5000106002215667]
u = [1.0000106002215667, 1.0000053000993734, 0.5, 1.5000106002215667]
u = [1.000110226799142, 1.0000551122077321, 0.5, 1.500110226799142]
u = [1.000110226799142, 1.0000551122077321, 0.5, 1.500110226799142]
u = [1.0011103465192357, 1.0005550530743994, 0.5, 1.5011103465192357]
u = [1.0011103465192357, 1.0005550530743994, 0.5, 1.5011103465192357]
u = [1.0111812496619466, 1.0055784987258456, 0.5, 1.5111812496619466]
u = [1.0111812496619466, 1.0055784987258456, 0.5, 1.5111812496619466]
u = [1.1187526052593677, 1.0580651061192092, 0.5, 1.6187526052593677]
u = [1.1187526052593677, 1.0580651061192092, 0.5, 1.6187526052593677]
u = [1.5510879326571247, 1.2534972603285255, 0.5, 2.0510879326571247]
u = [1.5510879326571247, 1.2534972603285255, 0.5, 2.0510879326571247]
u = [2.073539142920032, 1.4760991255863631, 0.5, 2.573539142920032]
u = [2.073539142920032, 1.4760991255863631, 0.5, 2.573539142920032]
u = [4.241422560753108, 2.275828293513783, 0.5, 4.741422560753108]
u = [4.241422560753108, 2.275828293513783, 0.5, 4.741422560753108]
u = [8.126721191149608, 3.502318797167331, 0.5, 8.626721191149608]
u = [8.126721191149608, 3.502318797167331, 0.5, 8.626721191149608]
u = [14.774118652607996, 5.355805713257135, 0.5, 15.274118652607996]
u = [14.774118652607996, 5.355805713257135, 0.5, 15.274118652607996]
u = [26.033186639636995, 8.164553441537473, 0.5, 26.533186639636995]
u = [26.033186639636995, 8.164553441537473, 0.5, 26.533186639636995]
u = [44.988891465130465, 12.437837366819815, 0.5, 45.488891465130465]
u = [44.988891465130465, 12.437837366819815, 0.5, 45.488891465130465]
u = [76.67675477680977, 18.950320654210778, 0.5, 77.17675477680977]
u = [76.67675477680977, 18.950320654210778, 0.5, 77.17675477680977]
u = [44.988891465130465, 12.437837366819815, 0.5, 45.488891465130465]
u = [76.67675477680977, 18.950320654210778, 0.5, 77.17675477680977]
u = [48.72009756530117, 13.238872149952986, 0.5, 50.20289532511152]
u = [48.54336603491433, 13.201230296816952, 0.5, 49.98835827041415]
u = [48.55303206590665, 13.203289776541785, 0.5, 50.00011302800904]
u = [48.552937920367874, 13.203269717955159, 0.5, 49.99999855038348]
u = [48.55293911366616, 13.203269972198559, 0.5, 50.0000000013932]
u = [48.552939112518786, 13.203269971954098, 0.5, 49.999999999998025]
u = [48.55293911252043, 13.203269971954446, 0.5, 50.000000000000014]
u = [48.55293911252039, 13.203269971954441, 0.5, 49.99999999999997]
u = [48.55293911252039, 13.203269971954441, 0.5, 49.99999999999997]
5

That’s increasing it

6

My apologies. This shows that x_1 approaches the expected value 49.50 as I reduce the tolerance. The differential state still remains outside of tolerance, however.

1e-8:

u = [49.479690908110484, 13.682966651643868, 0.5, 49.99999999999998]

1e-9:

u = [49.4915514356334, 13.685537303140494, 0.5, 49.99999999999999]

1e-10:

u = [49.494069541722695, 13.686082918800729, 0.5, 49.99999999999998]

1e-11:

u = [49.49822477214678, 13.686983072701429, 0.5, 49.999999999999986]

1e-12:

u = [49.49990766890697, 13.687347610721506, 0.5, 49.99999999999997]

I saw Yingbo’s BDF method derivation reference in his blog. Do you recommend any references for understanding of algebraic state interpolation?