feat: allow parameters in ODESystem to be unknowns in initialization system

[AayushSabharwal]

Close #2665

Requires:

• https://github.com/SciML/OrdinaryDiffEq.jl/pull/2228
• https://github.com/SciML/DiffEqBase.jl/pull/1051
• https://github.com/SciML/SciMLBase.jl/pull/698

Checklist

• [ ] Appropriate tests were added
• [ ] Any code changes were done in a way that does not break public API
• [ ] All documentation related to code changes were updated
• [ ] The new code follows the contributor guidelines, in particular the SciML Style Guide and COLPRAC.
• [ ] Any new documentation only uses public API

Additional context

Add any other context about the problem here.

[ChrisRackauckas]

Why would parameters be unknowns?

[hersle]

This is the way I see it. Consider a system of ODEs where some can be solved analytically. Suppose also that their initial values are constrained. For example (with "alternative 1"):

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

@parameters y0 # represents y(0) in alternative 2 below
@variables x(t) y(t)
@mtkbuild sys = ODESystem([
    D(x) ~ sin(y^2)  # hard to solve analytically for arbitrary y
    D(y) ~ t         # alternative 1: use ODE (works)
    # y ~ t^2/2 + y0 # alternative 2: use analytical solution of ODE (fails)
], t; initialization_eqs = [y * tan(y) ~ 1.0]) # non-trivial initialization constraint
prob = ODEProblem(sys, [x => 0.0], (0.0, 1.0), []; guesses = [y, y0] .=> 1.0)
sol = solve(prob, Tsit5()) # expect sol[y][begin] ≈ 0.8603335890193797

This works as expected: y is correctly initialized from the constraint.

But what if I want to implement the same system, but with the analytical solution for y instead? This would reduce the number of ODEs and improve performance. Then I need a parameter y0 to represent the initial value of y(0), and this would have to be treated as unknown during initialization. I expect the same code to work, only with D(y) ~ ... commented and # y ~ ... uncommented. But this fails precisely because y0 is not treated as unknown.

Note that if the initialization equation was something simpler, like y ~ 1.0 - x, this could be handled in both alternative 1 and 2 with defaults for both [y, y0] .=> 1.0 - x. In this case, one can interpret y0 as an "unknown" parameter. This issue can thus be seen as a matter of generalizing unknown parameters to work with arbitrary initialization equations (harder than substitution).

In other words, I think allowing for mixed "numerical/analytical" ODEs in a consistent way requires unknowns and parameters to be treated on equal footing during initialization. Looking (very) far ahead, I also think overcoming this is a prerequisite to solve e.g. #1586.

[AayushSabharwal]

This PR is too old and will be reimplemented