Solve does not solve initialization problem of remade problems

[hersle]

I want to solve an ODESystem for many values of a parameter P, where the initial value for an unknown x depends on P (here I write a nonlinear initialization equation to emphasize the general case where defaults do not suffice):

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

@variables x(t)
@parameters P
@named sys = ODESystem([D(x) ~ 0], t, [x], [P]; initialization_eqs = [x^2 ~ P], guesses = [x => +1.0])
ssys = structural_simplify(sys)

prob1 = ODEProblem(ssys, [], (0.0, 1.0), [P => 1.0])
prob2 = remake(prob1, p = [P => 4.0])
sol2 = solve(prob2)
@test sol2[x][begin] == 2.0

This fails because remake (or the following solve) does not resolve the initialization problem. Currently, I am resorting to reconstructing the ODEProblem with new parameters every time, which is very slow.

I think there should be an efficient and well-supported way to remake and then (re)solve an ODEProblem, where the full initialization system is also solved before the ODE is integrated.

[hersle]

Is a fix to this already underway with the PRs linked in #2747?
If so I apologize for the spam 😅 Anyway, I think the example above motivates this feature very well.

[bradcarman]

I believe that PR is for another feature. I'm also finding this bug, I added some color to your MWE and provided a proof...

@variables x(t)
@parameters P
@mtkbuild sys = ODESystem([D(x) ~ 0], t, [x], [P]; initialization_eqs = [x^2 ~ P], guesses = [x => +1.0])

prob1 = ODEProblem(sys, [], (0.0, 1.0), [P => 1.0])
prob2 = remake(prob1, p = [P => 4.0])
sol1 = solve(prob1)
sol2 = solve(prob2)

sol1(0; idxs=x) #1.0 OK
sol2(0; idxs=x) #1.0 <-- ERROR, should be 2.0

# Building and solving the initialization system explicitly gives the correct result
initsys = ModelingToolkit.generate_initializesystem(sys)
initsys = structural_simplify(initsys)
initprob = NonlinearProblem(initsys, [t=>0],[P => 1.0])
initsol = solve(initprob)
initsol[x] #1.0 OK

initprob = remake(initprob; p = [P => 4.0])
initsol = solve(initprob)
initsol[x] #2.0 OK

[bradcarman]

New PR that hopefully fixes this is #2928

[ChrisRackauckas]

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

@variables x(t)
@parameters P
@mtkbuild sys = ODESystem([D(x) ~ 0], t, [x], [P]; initialization_eqs = [x^2 ~ P], guesses = [x => +1.0])

prob1 = ODEProblem(sys, [], (0.0, 1.0), [P => 1.0])
prob2 = remake(prob1, p = [P => 4.0])
sol1 = solve(prob1)
sol2 = solve(prob2)

sol1(0; idxs=x) #1.0 OK
sol2(0; idxs=x) #2.0

# Building and solving the initialization system explicitly gives the correct result
initsys = ModelingToolkit.generate_initializesystem(sys)
initsys = structural_simplify(initsys)
initprob = NonlinearProblem(initsys, [t=>0],[P => 1.0])
initsol = solve(initprob)
initsol[x] #1.0

initprob = remake(initprob; p = [P => 4.0])
initsol = solve(initprob)
initsol[x] #2.0 OK

fixed on the latest version by that PR.