ModelingToolkit.jl
ModelingToolkit.jl copied to clipboard
fix linearization t0 and bump di
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
The part that's a little weird here is that this api is intended to work for
- autonomous systems in which case it doesn't want to force the user to give a time
- non-autonomous systems in which case a time is required.
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.
What if we build with tspan = (nothing, nothing) and use promote the eltype of u0 and the tunables for preparation?
That would be fine if most of the ODEProblem is discarded.
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.
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}
LGTM