ModelingToolkit.jl icon indicating copy to clipboard operation
ModelingToolkit.jl copied to clipboard

fix linearization t0 and bump di

Open oscardssmith opened this issue 5 months ago • 6 comments

oscardssmith avatar Jun 12 '25 20:06 oscardssmith

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

baggepinnen avatar Jun 13 '25 05:06 baggepinnen

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.

oscardssmith avatar Jun 13 '25 06:06 oscardssmith

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.

ChrisRackauckas avatar Jun 13 '25 09:06 ChrisRackauckas

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

AayushSabharwal avatar Jun 13 '25 09:06 AayushSabharwal

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

ChrisRackauckas avatar Jun 13 '25 09:06 ChrisRackauckas

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 avatar Jun 13 '25 14:06 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}

oscardssmith avatar Jun 16 '25 15:06 oscardssmith

LGTM

oscardssmith avatar Jul 02 '25 12:07 oscardssmith