Broken amount of parameters: much more parameters in resulting ODEProblem than symbolic parameters in MTK system

[Datseris]

Describe the bug 🐞

I just updated my MTK version and I find an incredible odd behavior. The number of parameters in the MTK model that generates the ODE Problem is much less than the parameters in the ODE Problem itself.

I have

julia> prob
ODEProblem with uType SVector{5, Float64} and tType Float64. In-place: false
Initialization status: FULLY_DETERMINED
Non-trivial mass matrix: false
timespan: (0.0, Inf)
u0: 5-element SVector{5, Float64} with indices SOneTo(5):
  290.0
    1.0
 1200.0
  290.0
   11.0

All of my parameters are real numbers so by default they become "tunable", and I have:

julia> prob.p.tunable
79-element Vector{Float64}:
 1.0
 0.38
 0.0012
 1.0
 5.0
 0.3
 0.2
 0.0
 ⋮

This is much more than the actual parameters of my model. Indeed, If I check the symbolic parameters in the MTK System I get

julia> mtk = prob.f.sys
Model CTBBL:
Equations (5):
  5 standard: see equations(CTBBL)
Unknowns (5): see unknowns(CTBBL)
  SST(t) [defaults to 290.0]: SCT sea surface temperature
  C(t) [defaults to 1.0]: cloud fraction
  z_b(t) [defaults to 1200.0]: height of the cloud top = top of the boundary layer = at inversion, m
  s_b(t) [defaults to 290.0]: liquid water static energy of boundary layer, K (cₚ units)
  ⋮
Parameters (25): see parameters(CTBBL)
  e_D [defaults to 1.0]: height depression efficiency.
  ⋮

25 parameters is the correct number.

Expected behavior

The number of parameters in the ODE problem should be the same.

There must have been a breaking change because before I updated I didn't have this problem.
Note that I recovered this problem because I build all my models with split = false, so that I do NOT create the specialized MTK container, ecll my parameters are real numbers. Even in that case I have the huge missmatch of 70+ parameters in the created ODEProblem even though my system only has 25.

What was the change that caused this? Is there perhaps a new API way to obtain the "real" parameters of an ODEProb created from an ODESystem?

Minimal Reproducible Example 👇

If you really need a MRE for this, I'll try to make one from scratch. My codebase is now too complicated to modify it into a MRE. But I doubt a MRE is needed, this seems to be a deliberate change.

Environment (please complete the following information):

• Output of using Pkg; Pkg.status()

pkg> st -m ModelingToolkit Symbolics SymbolicUtils SymbolicIndexingInterface
Status `...\Manifest.toml`
  [961ee093] ModelingToolkit v9.64.3
  [2efcf032] SymbolicIndexingInterface v0.3.38
  [d1185830] SymbolicUtils v3.16.0
  [0c5d862f] Symbolics v6.29.2

• Output of versioninfo()

1.10.8

[Datseris]

Well, here is a MRE as well:

using ModelingToolkit

# @variables t # use unitless time
D = Differential(t)
@mtkmodel Roessler begin
    @parameters begin
        a = 0.2
        b = 0.2
        c = 5.7
    end
    @variables begin
        x(t) = 1.0
        y(t) = 0.0
        z(t) = 0.0
        nlt(t) # nonlinear term
    end
    @equations begin
        D(x) ~ -y -z
        D(y) ~ x + a*y
        D(z) ~ b + nlt
        nlt ~ z*(x - c)
    end
end

@mtkbuild model = Roessler()
prob = ODEProblem(model)
prob.p.tunable

7-element Vector{Float64}:
 0.2
 0.2
 5.7
 0.0
 0.0
 1.0
 0.0

this is wrong, clearly only 3 parameters exist. is there a way to ensure that only the existing parameters go into the model? Especially when using split = false?

[Datseris]

Indeed, changing to:

@mtkbuild model = Roessler() split = false
prob = ODEProblem(model)
prob.p

7-element Vector{Float64}:
 0.2
 0.2
 5.7
 0.0
 1.0
 0.0
 0.0

Where are these extra parameters coming from...?

[SebastianM-C]

What happens here is that you also have the initial values for the unknowns of the model, which are also parameters. In previous versions of MTK they were not included.

[SebastianM-C]

And with #3439 the will move to a separate section in the MTKParameters.

[Datseris]

Oh thank god there is a coming fix for that, this is absolutely breaking for me, and may I say it does not make conceptual sense from a dynamical systems point of view to make default values or initial conditions extra parameters.

[Datseris]

also, my initial values are 3, not 4, which is another thing I would consider incorrect. nlt in the MRE does not have a default value and does not need one as it is not a state variable.

[Datseris]

I hope the MTK team will consider at least an option so that when split = true is given these "default values" are not made into extra parameters so that the length of hte parmaeter container matches exactly the symbolic @parameters of hte system.

[SebastianM-C]

and may I say it does not make conceptual sense from a dynamical systems point of view to make default values or initial conditions extra parameters.

Why would they not be parameters? If you consider the dynamical system in general, it can be parameterized by initial conditions and you can optimize those. In previous version if we were optimizing the initial conditions we would have passed the unkonwn to remake/setsym, but rigorously speaking we don't remake the unknown, but the initial value, which is a parameter.

so that the length of hte parmaeter container matches exactly the symbolic @parameters of hte system.

I'm curious why is that needed?

[Datseris]

Why would they not be parameters? If you consider the dynamical system in general, it can be parameterized by initial conditions and you can optimize those. In previous version if we were optimizing the initial conditions we would have passed the unkonwn to remake/setsym, but rigorously speaking we don't remake the unknown, but the initial value, which is a parameter

This statement does not make sense to me. How can the initial condition be a parameter of a dynamical system? Does the initial condition influence the form of the state space? Does it influence the nullclines shapes? Does it influence the speed of motion in the state space? Of course not. Does changing an initial condition change when bifurcations happen for example? of course not. The dynamical system is the whole of the state space. This is what it means to consider the dynamical system "in general".

What I think you are calling a "dynamical system" is not the dynamical system in general, but rather a specific trajectory of a it evolved forwards in time and for a finite time window. In this scenario sure you can "optimize" the initial condition so that the trajectory does something else. If your time window is large enough, such an optimization only makes sense in a multistable dynamical system where different initial conditions may converge to different final states. But in a monostable system, what difference does the initial condition make? A parameter change on the other hand has huge impact in both monostable or multistable systems.

---

I'm curious why is that needed?

You are right, it is not needed, if I always and only access parameters by their symbolic reference. If however, for some reason research related or otherwise, I just want to obtain the whole parameter container and e.g., count how many parameters my system has, then you can imagine how confusing it is to have 80 instead of 25. Imagine writing a paper and needing to make a statement "there are 25 free parameters". You can't count the 25 by hand in a complex real-world system like mine.

And as I said above, only for true state variables it makes sense to even consider their starting values. For non-state variables as the nlt above I see no reason why you would even store its starting value.

[Datseris]

In any case, if we have contrasting views on what a "parameter" is, it doesn't make sense to argue about it. Everyone can use the definition that fits their needs.

My point is: since there is an explicit @parameter declarator, all I ask is that my parameter container only has what I explicitly declared as @parameter, or at least there is an API-way to obtain those from the ODEProblem itself, or the final "numeric" container of parameters, whatever this structure is.

[SebastianM-C]

What I think you are calling a "dynamical system" is not the dynamical system in general, but rather a specific trajectory of a it evolved forwards in time and for a finite time window.

The ODEProblem is is precisely that, a specific trajectory with a finite time window, which is what contains the new Initials, so I think we are in agreement 😅

There is a new isinitial function that can help filter these new parameters, that should help (#3433).

For nlt I think the initial value might be stored due to initialization since you can initialize your system by asking for a specific value of that and MTK should solve the nonlinear system to find out how the rest of the variables should start.

[SebastianM-C]

I think that the main distinction here is initial value problem vs dynamical system. So when you have the dynamical system, as you said the initial conditions are less relevant, but when you solve an initial value problem, then the initial conditions become parameters. I think this is illustrated in how you go from a ODESystem to an ODEProblem.

See also #3397 for more discussion on the topic 😅

[Datseris]

I think that the main distinction here is initial value problem vs dynamical system

That's the perfect way to put it I think. For a dynamical system the initial conditions are not parameters but for the IVP they are. I see people in #3397 having similar problems due to this missmatch. I think this is good evidence that some explicit API function that differentiates is useful. It is good to see that there will be a function isinitial, but this is not enough. The discussion in #3397 makes it clear that what is the content of prob.p may change arbitrarily in the future. This is not safe for me, and I'd like to avoid what happened to me today happen again.

So: is there, or will there ever be, a function that returns all parameters generated explicitly via @parameters? Is this function ModelingToolkit.parameters(sys)? If yes, how can I make this function apply to ODEProblem? Will it ever be applicable?
In the prob of the previous MRE I get:

julia> ModelingToolkit.parameters(prob)
Any[]

which is kinda nonsense. It would be better if it erroed with MethodError.

[hersle]

I think that the main distinction here is initial value problem vs dynamical system

That's the perfect way to put it I think.

This viewpoint also generalizes nicely to a boundary-value problem, where it can make sense to refer to both Initial(x(t)) and Final(y(t)) (for example; the latter doesn't exist now) as parameters.

[Datseris]

Good point.

Staying on the topic of the issue, the issue was closed but this was not answered:

is there, or will there ever be, a function that returns all parameters generated explicitly via @parameters? Is this function ModelingToolkit.parameters(sys)? If yes, how can I make this function apply to ODEProblem? Will it ever be applicable?

[ChrisRackauckas]

This viewpoint also generalizes nicely to a boundary-value problem, where it can make sense to refer to both Initial(x(t)) and Final(y(t)) (for example; the latter doesn't exist now) as parameters.

And even further, the way the fitting of parameters is done is to treat the entire thing as a single vector that is then solved via a NonlinearLeastSquaresProblem. Formally, parameters are just initial conditions of state values where the derivative are zero. So a "parameter" is just p(0) where dp/dt = 0. In that sense, there really isn't a difference between them. Everything like the adjoint methods use this equivalence in the derivation BTW, so it's pretty fundamental to a lot of the numerics that initial conditions and parameters are treated similarly. We tend to differentiate between them in the high level interfaces but many of the internals will just treat it as just a single vector.

But anyways, details aside:

1. We document to not rely on the exact formulation of prob.p. It can and will change. SII is recommended for all operations there.
2. We have split initials to a separate vector. This is because it can be specialized in the adjoint (du0/dp is equal to the pullback lambda and thus you can shortcut some computations on their derivative using a few tricks).

[Datseris]

Okay, but can someone please explain: what is the way to obtain all @parameters in an MTK system? Will ModelingToolkit.parameters(sys) guaranteed return only the parameters defined with @parameters, or will it return the full parameter container which would include initial conditions?

[SebastianM-C]

MTK has both parameters and full_parameters. Should Initials only be accessible via full_parameters?

IIRC this was introduced for parameter dependencies and caused some bugs for a while since it had to be replaced/propagated in a lot of places internally.