Custom initial condition equations for an ODESystem

[MarcBerliner]

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

This change allows specifying custom initial condition equations for an ODE model which requires metadata for the initial guess. The initializesystem function takes an ODESystem as an argument and returns a NonlinearSystem.

TODO: I added docs for the guess variable metadata, but this still needs docs for initializesystem. Where is an appropriate place for these docs?

https://github.com/JuliaComputing/JuliaSim.jl/issues/747

[ChrisRackauckas]

Yeah let's keep it out of the docs for now until it's tried a bit more but this looks right.

[mtiller-jh]

If I understand these changes, the initial equations are being associated with variables. This doesn't feel right to me. From the test:

        (x(t) = dx ~ dx_ss), [guess = 0.5]
        (y(t) = dy ~ 0), [guess = -0.5]

These equations have nothing to do with these variables. Why associated them like this? Also, what happens when you want to add initial equations to a system where the variables have already been defined? Here is an example of exactly this from my book:

model QuiescentModelWithInheritance "Steady state model with inheritance"
  extends ClassicModel;
initial equation
  der(x) = 0;
  der(y) = 0;
end QuiescentModelWithInheritance;

Is there a reason we are associated these equation directly with variables? I just don't think this would scale very well. Why not simple include an additional field in ODESystem called init_eqs, e.g.,

    # Define the Lotka Volterra system which begins at steady state
    @parameters t
    pars = @parameters a=1.5 b=1.0 c=3.0 d=1.0 dx_ss = 1e-5

    vars = @variables begin
        dx(t),
        dy(t),
        x(t), [guess = 0.5]
        y(t), [guess = -0.5]
    end

    D = Differential(t)

    init_eqs = [dx ~ dx_ss
                      dy ~ 0]
    eqs = [dx ~ a * x - b * x * y
        dy ~ -c * y + d * x * y
        D(x) ~ dx
        D(y) ~ dy]

    @named sys = ODESystem(eqs, t, vars, pars; init_eqs=init_eqs)

[ChrisRackauckas]

Also, what happens when you want to add initial equations to a system where the variables have already been defined?
Is there a reason we are associated these equation directly with variables? I just don't think this would scale very well.
Why not simple include an additional field in ODESystem called init_eqs, e.g.,

Just use the defaults as before. It still does the equivalent thing.

[mtiller-jh]

| Just use the defaults as before. It still does the equivalent thing.

@ChrisRackauckas , I don't really understand this comment. I don't understand what you mean by before. There was not previously a way to specify initial condition equations, right?

Keep in mind, I am not talking about how to programmatically formulate guess values for iteration variables. I'm talking about imposing hard constraints on the values of the initial conditions for the IVP. That is what this change is about, right? What does this have to do with "defaults"?

[ChrisRackauckas]

No, specification of initial equations existed before and this uses the same representation without change. All this does is add the building of a nonlinear system from them.

[mtiller-jh]

Can someone (@MarcBerliner , @ChrisRackauckas , @bradcarman ) show me an example of formulating a non-linear system that combines the differential equations and the initial equations "[using] the defaults as before"?

[MarcBerliner]

Previously, we could only specify a linear system of initial equations for IVPs. For example,

vars = @variables begin
    dx(t),
    dy(t),
    x(t) = 1, [guess = 0.5]
    y(t) = 2, [guess = -0.5]
end

creates a system equivalent to

vars = @variables begin
    dx(t),
    dy(t),
    (x(t) = x ~ 1), [guess = 0.5]
    (y(t) = y ~ 2), [guess = -0.5]
end

where 1 and 2 are the "defaults" for x and y respectively, or in other words, x(t0) = 1 and y(t0) = 2. There was no need for a guess because all unknowns would start at their default.

This PR allows the user to write arbitrary initial condition equations. The line (x(t) = x ~ 1) means "add x(t0) to the list of initial unknowns along with all algebraic unknowns" and "add the equation (x) - (1) = 0 to the list of initial equations along with all algebraic equations."