fix linearization t0 and bump di

[oscardssmith]

No description provided.

[baggepinnen]

I'm not sure I understand the implications of t = nothing, but it's very common to want to linearize systems that are not autonomous, any control system that uses an input-signal block falls into this category

[oscardssmith]

The part that's a little weird here is that this api is intended to work for

1. autonomous systems in which case it doesn't want to force the user to give a time
2. non-autonomous systems in which case a time is required.

[ChrisRackauckas]

Okay it's odd to build an ODEProblem with tspan nothing if it's supposed to be a value, so yes we should never do that. But, it looks like tests are still having issues with DI bump.

[AayushSabharwal]

What if we build with tspan = (nothing, nothing) and use promote the eltype of u0 and the tunables for preparation?

[ChrisRackauckas]

That would be fine if most of the ODEProblem is discarded.

[oscardssmith]

I think the only reason we're building an ODE in the first place is to get OrdinaryDiffEq jacobian handling. Now that DI exists, we could probably just use it directly. That said, I think this is a good fix in the short term.

[oscardssmith]

so MWE of the issue for @AayushSabharwal:

julia> using ModelingToolkit, ModelingToolkitStandardLibrary,  Symbolics, Test
julia> using ModelingToolkit: t_nounits as t, D_nounits as D
julia> using ModelingToolkitStandardLibrary.Mechanical.Rotational
julia> @named inertia1 = Inertia(; J = 1);
julia> @named inertia2 = Inertia(; J = 1);
julia> @named spring = Rotational.Spring(; c = 10);
julia> @named damper = Rotational.Damper(; d = 3);
julia> @named torque = Torque(; use_support = false);
julia> @variables y(t) = 0;
julia> eqs = [connect(torque.flange, inertia1.flange_a)
              connect(inertia1.flange_b, spring.flange_a, damper.flange_a)
              connect(inertia2.flange_a, spring.flange_b, damper.flange_b)
              y ~ inertia2.w + torque.tau.u];
julia> model = System(eqs, t; systems = [torque, inertia1, inertia2, spring, damper],
           name = :name, guesses = [spring.flange_a.phi => 0.0]);
julia> model_outputs = [inertia1.w, inertia2.w, inertia1.phi, inertia2.phi];
julia> model_inputs = [torque.tau.u];
julia> op = Dict(torque.tau.u => 0.0);
julia> lin_fun, ssys = linearization_function(model,
           model_inputs,
           model_outputs;
           op);
julia> linearize(ssys, lin_fun; op)

gives

ERROR: Non-concrete element type inside of an `Array` detected.
Arrays with non-concrete element types, such as
`Array{Union{Float32,Float64}}`, are not supported by the
differential equation solvers. Anyways, this is bad for
performance so you don't want to be doing this!

Specifically, the problem comes from SciMLBase.get_initial_values, and specifically, the call to initdata.update_initializeprob!(initprob, valp) which decides that u0 should be a Vector{Union{Nothing, Float64}} instead of a Vector{Float64}

[AayushSabharwal]

Why was this kwarg changed? That's the root cause for the failures

[oscardssmith]

#3733 (comment). The previous problem was that we were preparing the ODE with t=nothing and running with t=0 which meant our jacobian was getting an unexpected t type.

[oscardssmith]

Basically the problem is that linearize has this weird spec where we aren't requiring a t for autonomous systems, but otherwise requires a t to be given.

[AayushSabharwal]

So defaulting to t = 0 for autonomous systems shouldn't be a problem?

[AayushSabharwal]

MTK#master right now accepts t as a kwarg to linearization_function, and forwards it from linearize. So that should just address this?

[oscardssmith]

that's a lot better. The basic issue is that we don't want t=0.0 to be passed by default to non-autonomous systems. For them, we want to error if the user doesn't provide t

[AayushSabharwal]

I still don't get why we don't want to default to t=0.0 for non-autonomous systems 😅 The higher level API for this is linearize and that defaults to t=0.0. And the t provided here is essentially only for typing the DI jacobian prep. The lin_fun::LinearizationFunction returned from linearization_function is called as (u, p, t) which requires passing time, and if the user calls linearize(lin_fun, sys) they can also pass a new t there. This behavior is essentially no different from what we had before. Making t = nothing, defaulting to zero for autonomous and erroring non-autonomous is an orthogonal change and can be done separately.

[oscardssmith]

This PR also removes that as a default unless I messed something up. But also, why should it be Float64 only? If the user passes a Float32 time, do we expect an error here (due to incorrect prep)?

[baggepinnen]

Defaulting to t=0 is perfectly reasonable, the point in which to linearize is also by default the same as the initial condition. Changing this is also breaking since this has been the default before.

[oscardssmith]

is t=0 guaranteed to be an initial condition?