Solving a system with initial condition of the form $f(0) = g(0)$

[natema]

In some cases when solving a model we have a differential-algebraic system of equations where two variables and have an initial condition of the form .
We created a thread in the Julia forum to ask how to input such initial condition in ModelingToolkit, and it was suggested to open an issue here.

Here's a MWE which returns an error:

using ModelingToolkit, DifferentialEquations

@variables t
D = Differential(t)

function system_f(; name)
    @variables f(t) g(t)
    eqs = [D(f) ~ g - f]
    ODESystem(eqs; name)
end

function system_g(; name)
    @variables g(t)
    eqs = [g ~ t]
    ODESystem(eqs; name)
end

@named f = system_f()
@named g = system_g()

connected_eqs = [
    f.g ~ g.g,
]

@named _model = ODESystem(connected_eqs, t)
@named model = compose(_model, [f, g])

model = structural_simplify(model)
prob = ODEProblem(model, [f.f => g.g], (0.0, 10.0))
sol = solve(prob)

If we replace the last two lines in the above block with the following ones, the code works, but the solution is clearly a workaround that could be improved:

pre_prob = ODEProblem(model, [f.f=>42,g.g=>42], (0.0, 10.0)) 
pre_sol = solve(pre_prob)
prob = ODEProblem(model, [f.f=>pre_sol[g.g, 1],g.g=>pre_sol[g.g, 1]], (0.0, 10.0))
sol = solve(prob)

[natema]

With @piluc we distilled an even simpler MWE:

using ModelingToolkit, DifferentialEquations

@variables t
D = Differential(t)

@variables f(t) g(t)

eqs = [D(f) ~ g - f, g ~ t]

@named sys = ODESystem(eqs)
sys = structural_simplify(sys)
prob = ODEProblem(sys, [f => g], (0.0, 10.0))
sol = solve(prob)

This problem is particularly bothersome when trying to translate Vensim models to MTK, since initializations of the form f => g are often used there (and they allow for more concise formulations).

[ChrisRackauckas]

Why is this an issue? The issue is that it's just underspecified. prob = ODEProblem(sys, [f => g, g => 0.0], (0.0, 10.0)) and you should be good.

[natema]

... Oops!

[natema]

Actually, @ChrisRackauckas, why is the system unspecified?
The above example considers the system $f' = g - f, g = t$ and declares as initial condition $f(0) = g(0)$, where $g(0) = 0$.
In fact, in the proposed workaround of the OP, I do provide an initial condition for $g$, which gets ignored.

AFAIU, it looks like in MTK it is not possible to provide as initial condition the initial value of an algebraic equation, even if the latter can be computed without solving the ODEs.

[paulobruno]

As an example, if we do prob = ODEProblem(sys, [f => g, g => 30.0], (0.0, 10.0)), the initial value of g is 0.0, as expected, but the initial value of f is 30.0. Instead, we would like the initial value of f to be 0.0 as well.

In the MWE, it's clear that $g(0) = 0$, but in reality $g$ is a function that depends on many other values.

[ChrisRackauckas]

No, prob = ODEProblem(sys, [f => g, g => 30.0], (0.0, 10.0)), here the initial value of g is 30.0 and the initial value of f is whatever the initial value of g is, which is 30.0.

but in reality is a function that depends on many other values.

So then do exactly as shown in the example here. Here, the initial value of f is dependent on many other values, which happens to be g. But then you can make the initial value for g just be a function of the parameters, say g => p1 * p2 etc. There's a whole system for functional initial conditions, just use it. But 2 things needs two initial values: saying that one is the initial condition of the other is underspecified until you define the other, either as a constant or a function of parameters.

[natema]

@ChrisRackauckas I don't think this issue should be closed, since the situation remains confusing w.r.t. what a user can understand at this point about how MTK works.
In particular, let me remark again that even if the provided initial value of g is 30, the initial value of g in the solution will be different. Isn't this confusing?
Would it help if we provide a MWE in Vensim and show how we translate it in MTK?

[ChrisRackauckas]

let me remark again that even if the provided initial value of g is 30, the initial value of g in the solution will be different. Isn't this confusing?

It won't be unless it's an inconsistent initialization of a DAE, which then runs the chosen initialization algorithm which will change the differential or algebraic variables as chosen by the user (defaulting to only the algebraic). If that's not the case, provide and MTK MWE of what's not working. There has been no example shown that this doesn't work, and there are tested examples that this does work.