Compute initial states for `u0` if `u0` is parametric

[MarcBerliner]

If u0 is known to be a parametric function, changing the parameters will affect the initial values (e.g., for ODEs) and/or the initial guess for algebraic equations (e.g., DAEs). This problem can be avoided if MTK detects any parametric initial states (when declaring @variables u(t) = p or ODEProblem(sys, [u => p])) and generates a function to recompute the initial states upon declaring new parameters.

It would be helpful if this is included in remake (e.g., when calling remake(prob; p=...)) and as an internal function (e.g., generate_initial_states!(prob) or something).

MWE:

using ModelingToolkit, OrdinaryDiffEq
@parameters p = 1
@variables t u(t)

eqs = [Differential(t)(u) ~ 1]

@named sys = ODESystem(eqs,t,[u],[p],tspan=(0, 10))

prob = ODEProblem(sys, [u => p])

sol = solve(prob, Tsit5())
# sol.u[end] ≈ 11

# remake the problem with new parameters
prob_new = remake(prob; p=[p => 2])
# however, u0 is still 1
# prob_new.u0[1] == 1

sol_new = solve(prob_new, Tsit5())
# sol_new.u[end] ≈ 11, and not 12

An ideal workflow where I can define u(t) = p:

using ModelingToolkit, OrdinaryDiffEq
@parameters p = 1
@variables t u(t) = p

eqs = [Differential(t)(u) ~ 1]

@named sys = ODESystem(eqs,t,[u],[p],tspan=(0, 10))

prob = ODEProblem(sys)

sol = solve(prob, Tsit5())
# sol.u[end] ≈ 11

# remake the problem with new parameters
prob_new = remake(prob; p=[p => 2])
# now, u0 is 2
# prob_new.u0[1] == 2

sol_new = solve(prob_new, Tsit5())
# sol_new.u[end] ≈ 12

[baggepinnen]

Would get_u0_p be helpful (available on master only)? It will at least handle parameters that are functions of other parameters.
https://docs.sciml.ai/ModelingToolkit/dev/systems/ODESystem/#Composition-and-Accessor-Functions

If it does not completely solve your problem, maybe it could be baked into remake somehow?

[MarcBerliner]

Note: I edited the first comment

get_u0_p can help! It looks like it is doing the same thing under the hood.

I thought about this problem more, and the more significant issue is that I cannot easily abstract away the initial states inside components with parametric u0. For complex systems (such as batteries), end users will likely not even know the states inside the battery component and will definitely not know appropriate initial values for them. So, I cannot reasonably expect users to know anything about the states, only the parameters. If sys knew the default u0map, it would be easy to run something like get_u0_p for the new initial states inside remake(prob; p=...)) or generate_initial_states!(prob)

[baggepinnen]

We need to step up our game when it comes to initial conditions across the board, I'll create an issue label for this to make sure we track all relevant issues.

If sys knew the default u0map

It does, or am I not understanding you correctly? It only knows it in it's symbolic form, i.e., it knows u = p in your example, and get_u0_p is required to process it into the numeric array. Or are you asking for remake(prob; p=...)) to automatically handle the fact that u depends on p, and when new p are provided, new u0 are automatically computed? This feature would indeed be nice to have, and this is most likely what the user would expect to happen when they change the parameters.

[MarcBerliner]

It only knows it in it's symbolic form, i.e., it knows u = p in your example

The issue is that I believe it does not store the relationship u = p. When ODEProblem(sys, [u => p]) is evaluated, it computes u0 using the default or specified parameter values, stores u0 as an array of floats, and then discards u => p.

are you asking for remake(prob; p=...)) to automatically handle the fact that u depends on p, and when new p are provided, new u0 are automatically computed? This feature would indeed be nice to have, and this is most likely what the user would expect to happen when they change the parameters.

Yeah, this is precisely what I mean. prob should store the necessary mapping information to make remake(prob; p=...)) the only function a user needs to run.