Errors and wrong solutions when solving `NonlinearSystem` with `initialization_eqs`

[TorkelE]

From the latest understanding this is what one should use with structural_simplify? Although, admittedly, I am not really sure anymore.

Either case, there is problem with the solution, generating a mix of errors and wrong outputs when you check the solution.

using ModelingToolkit, OrdinaryDiffEqDefault, Test
using ModelingToolkit: t_nounits as t, D_nounits as D

# Defines the model.
@parameters k1 k2
@variables X1(t) X2(t)
@parameters Γ[1:1] = missing [guess = [1.0]]
eqs = [
    0 ~ k1 * (Γ[1] - X1) - k2 * X1
]
initialization_eqs = [
    X2 ~ Γ[1] - X1
]
@mtkbuild nlsys = NonlinearSystem(eqs, [X1, X2], [k1, k2, Γ]; initialization_eqs)

# Creates the problem.
u0 = [X1 => 1.0, X2 => 2.0]
ps = [k1 => 0.1, k2 => 0.2]
nlprob = NonlinearProblem(nlsys, u0, ps)

# Solves the problem and check the solution.
sol = solve(nlprob)
sol[X1] # -0.0 (should be 1.0).
sol[X2] # ERROR: ArgumentError: Symbol X2(t) is not present in the system. (should be 2.0).
sol.ps[k1] # 0.1 (Correct.)
sol.ps[k2] # 0.2 (Correct.)
sol.ps[Γ[1]] # -0.0 (should be 3.0)

Either this is a bug, or I still do not know how you are supposed to do this.

[ChrisRackauckas]

sol[X2] # ERROR: ArgumentError: Symbol X2(t) is not present in the system. (should be 2.0).

There is no equation for X2, so I don't see how or why this should be 2.0. You're missing an equation?

sol.ps[Γ[1]] # -0.0 (should be 3.0)

Since X2 is not defined, Γ[1] can be anything.

[TorkelE]

Ok, that might make a bit of sense. My understanding was that

initialization_eqs = [
    X2 ~ Γ[1] - X1
]

is the equation that allow us to compute X2 from Γ[1] (parameter) and X1 (unknown).

Previously I instead did

observed = [
    X2 ~ Γ[1] - X1
]

which still would be my preferred way, but I am also trying to explore alternatives.

[ChrisRackauckas]

is the equation that allow us to compute X2 from Γ[1] (parameter) and X1 (unknown).

It was pointed out in that conversation that it's a bit ambiguous since X2 ~ Γ[1] - X1 is not a parameter expression. The unambiguous form would be:

using ModelingToolkit, OrdinaryDiffEqDefault, Test
using ModelingToolkit: t_nounits as t, D_nounits as D

# Defines the model.
@parameters k1 k2
@variables X1(t) X2(t)
@parameters Γ[1:1] = missing [guess = [1.0]]
eqs = [
    0 ~ k1 * (Γ[1] - X1) - k2 * X1
    X2 ~ Γ[1] - X1
]
initialization_eqs = [
    Initial(X2) ~ Γ[1] - Initial(X1)
]
@mtkbuild nlsys = NonlinearSystem(eqs, [X1, X2], [k1, k2, Γ]; initialization_eqs)

# Creates the problem.
u0 = [X1 => 1.0, X2 => 2.0]
ps = [k1 => 0.1, k2 => 0.2]
nlprob = NonlinearProblem(nlsys, u0, ps)

# Solves the problem and check the solution.
sol = solve(nlprob)
sol[X1]
sol[X2]
sol.ps[k1]
sol.ps[k2]
sol.ps[Γ[1]]

This won't work quite yet as we need to setup Initial(X1) to refer correctly, but this is at least the plan.