Cannot solve for dependent parameters during initialization

[hersle]

Consider this (trivial) ODE for x and y:

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

@variables x(t) y(t) s(t)
sys1 = structural_simplify(ODESystem([D(x) ~ 0, D(y) ~ 0, s ~ x + y], t, [x, y, s], []; name=:sys))
prob1 = ODEProblem(sys1, [s => 1.0, x => 0.3], (0.0, 1.0))
sol1 = solve(prob1)
@test all(sol1[x] .≈ 0.3) && all(sol1[y] .≈ 0.7) # passes

This works. The constraint s(0) == x(0)+y(0) == 1 correctly initializes y(0) == 0.7. I do not have to specify a guess for y(0).

---

But consider an implementation of its (trivial) analytical solution, using parameters for the initial values x(0) == x0 and y(0) == y0:

@parameters x0 y0
sys2 = structural_simplify(ODESystem([x ~ x0, y ~ y0, s ~ x + y], t, [x, y, s], [x0, y0]; name=:sys))
prob2 = ODEProblem(sys2, [s => 1.0, x0 => 0.3], (0.0, 1.0))
sol2 = solve(prob2)

This does not work. The initialization gives an error that asks for a guess:

ERROR: LoadError: Initial condition underdefined. Some are missing from the variable map.
Please provide a default (`u0`), initialization equation, or guess       
for the following variables:

Any[y0]

But adding e.g. ; guesses = [y0 => 0.0] to ODEProblem(sys2, ... does not change anything.

Is it possible to make this work?

[hersle]

From what I can tell, this happens because InitializationProblem() calls generate_initializesystem(), which does not treat parameters as (possible) unknowns.

For example

isys2 = generate_initializesystem(sys2; u0map = [s => 1.0, x => 0.3])
equations(isys2)

shows that all the 5 necessary equations are there:

5-element Vector{Equation}:
 0 ~ 0.3 - x(t)
 0 ~ 1.0 - s(t)
 0 ~ x0 - x(t)
 0 ~ y0 - y(t)
 0 ~ x(t) + y(t) - s(t)

But

unknowns(isys2)

shows that x0 and y0 are not treated as unknowns:

3-element Vector{SymbolicUtils.BasicSymbolic{Real}}:
 x(t)
 y(t)
 s(t)

[ChrisRackauckas]

This is possible to solve and we should fix this.

[AayushSabharwal]

So for this we need:

• Recognize parameters with no initial value as unknowns
• Respect symbolic parameter defaults as equations (if no explicit value given)
• Create and store a function for getting values from the nonlinear problem and putting them into the ODEProblem

[hersle]

If possible, it would be great to test that everything works seamlessly through autodiff as well, using recent improvements like #2686. In case it can be of use, here is one suggestion for a test which I think should pass:

using Test
using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D
using DifferentialEquations
using ForwardDiff: derivative as ∇

@testset "Solve for unspecified parameters during initialization (issue #2665)" begin

@variables x(t) y(t) s(t)
@parameters x0 y0

# solve using variables
function y_of_x_vars(x0val)
    @named sys = ODESystem([D(x) ~ 0, D(y) ~ 0, s ~ x + y], t, [x, y, s], [])
    sys = structural_simplify(sys)
    prob = ODEProblem(sys, [s => 1.0, x => x0val], (0.0, 1.0))
    sol = solve(prob, Tsit5())
    return sol(0.123, idxs=y)
end

# solve using parameters
function y_of_x_pars(x0val)
    @named sys = ODESystem([x ~ x0, y ~ y0, s ~ x + y], t, [x, y, s], [x0, y0])
    sys = structural_simplify(sys)
    prob = ODEProblem(sys, [s => 1.0], (0.0, 1.0), [x0 => x0val])
    sol = solve(prob, Tsit5())
    return sol(0.123, idxs=y)
end

# solve analytically
function y_of_x_anal(x0val)
    return 1.0 - x0val
end

@test y_of_x_anal(0.3) ≈ y_of_x_vars(0.3) ≈ y_of_x_pars(0.3)
@test ∇(y_of_x_anal, 0.3) ≈ ∇(y_of_x_vars, 0.3) ≈ ∇(y_of_x_pars, 0.3)

end

Thank you! :smiley:

[AayushSabharwal]

~I have something working, but the problem here seems to be that the initialization problem is only solved for DAEs which this example isn't. I'll try and find a solution.~ EDIT: The initialization problem isn't called because the problem has no unknowns.

[ChrisRackauckas]

You need to specify all parameters, that's a given. This now works:

@variables x(t) y(t) s(t)
sys1 = structural_simplify(ODESystem([D(x) ~ 0, D(y) ~ 0, s ~ x + y], t, [x, y, s], []; name=:sys))
prob1 = ODEProblem(sys1, [s => 1.0, x => 0.3], (0.0, 1.0))
sol1 = solve(prob1)
@test all(sol1[x] .≈ 0.3) && all(sol1[y] .≈ 0.7) # passes

@parameters x0 y0
sys2 = structural_simplify(ODESystem([x ~ x0, y ~ y0, s ~ x + y], t, [x, y, s], [x0, y0]; name=:sys))
prob2 = ODEProblem(sys2, [s => 1.0, x0 => 0.3, y0 => 0.3], (0.0, 1.0))
sol2 = solve(prob2)

and was fixed by https://github.com/SciML/ModelingToolkit.jl/pull/2775

[hersle]

But the constraint s => 1.0 is not respected in the latter case. See also #2747.

[AayushSabharwal]

The functionality necessary for this is added in https://github.com/SciML/ModelingToolkit.jl/pull/2928

[hersle]

I am really looking forward to this 😀 Thanks!