From 796cc481ee98e19135baae35eda23b975dc60919 Mon Sep 17 00:00:00 2001 From: Fredrik Bagge Carlson Date: Tue, 29 Sep 2026 10:16:07 +0000 Subject: [PATCH] Fix `linearize` without inputs and name differentiated inputs from the linearization function `linearize(sys, [], outputs)` failed in ForwardDiff with an invalid tag, since the Jacobians with respect to the inputs were prepared for an empty, `Any`-typed input vector. The input vector now has the element type of the state and the input Jacobians are not evaluated when there are no inputs, so that `B` and `D` have zero columns. `LinearizationFunction` stores its input symbols, from which the diagnostic for differentiated inputs takes their names; it previously indexed `inputs(sys)`, which is empty when the inputs of a compiled system are parameters not marked as inputs. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_01QBgBrXDidWG1G2KF64gAte --- src/linearization.jl | 36 ++++++++++++++++++++++++------------ test/linearize.jl | 19 +++++++++++++++++++ 2 files changed, 43 insertions(+), 12 deletions(-) diff --git a/src/linearization.jl b/src/linearization.jl index 747a51b62b0..a3616469abc 100644 --- a/src/linearization.jl +++ b/src/linearization.jl @@ -236,13 +236,15 @@ function linearization_function( p = parameter_values(prob) t0 = current_time(prob) - inputvals = [prob.ps[i] for i in inputs] if u0 === nothing T = typeof(t0) else T = promote_type(eltype(u0), typeof(t0)) end + # Preparing the jacobians with respect to the inputs requires a vector of a concrete + # element type, also when there are no inputs. + inputvals = isempty(inputs) ? T[] : [prob.ps[i] for i in inputs] prob = _linearization_wrap_odeproblem_f(prob, T) ct0 = DI.Constant(T(t0)) u0T = if u0 === nothing @@ -295,7 +297,8 @@ function linearization_function( input_getter = getsym(prob, inputs) lin_fun = LinearizationFunction( - diff_idxs, alge_idxs, input_getter, length(inputs), length(unknowns(sys)), + diff_idxs, alge_idxs, input_getter, length(inputs), collect(SymbolicT, inputs), + length(unknowns(sys)), prob, h, u0 === nothing ? nothing : similar(u0, T), uf_jac, h_jac, pf_jac, hp_jac, initializealg, initialization_kwargs, initial_idxs_for_unknowns, collect(SymbolicT, loop_opening_params) @@ -425,6 +428,11 @@ mutable struct LinearizationFunction{ """ const num_inputs::Int """ + The input variables, or the parameters that represent them, in the order of the columns + of the input jacobians. + """ + const inputs::Vector{SymbolicT} + """ The number of unknowns in the linearized system. """ const num_states::Int @@ -537,20 +545,24 @@ function (linfun::LinearizationFunction)(u, p, t) end fg_xz = linfun.uf_jac(u, DI.Constant(p), DI.Constant(t)) h_xz = linfun.h_jac(u, DI.Constant(p), DI.Constant(t)) - fg_u = linfun.pf_jac( - input_vals, - DI.Constant(u), DI.Constant(p), DI.Constant(t) - ) + # The jacobian with respect to an empty input vector is not defined. + fg_u = if linfun.num_inputs == 0 + zeros(eltype(fg_xz), size(fg_xz, 1), 0) + else + linfun.pf_jac(input_vals, DI.Constant(u), DI.Constant(p), DI.Constant(t)) + end else linfun.num_states == 0 || error("Number of unknown variables (0) does not match the expected number of unknowns ($(linfun.num_states))") fg_xz = zeros(0, 0) - h_xz = fg_u = zeros(0, length(linfun.num_inputs)) + fg_u = zeros(0, linfun.num_inputs) + h_xz = zeros(size(linfun.hp_jac.buf, 1), 0) + end + h_u = if linfun.num_inputs == 0 + zeros(eltype(h_xz), size(h_xz, 1), 0) + else + linfun.hp_jac(input_vals, DI.Constant(u), DI.Constant(p), DI.Constant(t)) end - h_u = linfun.hp_jac( - input_vals, - DI.Constant(u), DI.Constant(p), DI.Constant(t) - ) return ( f_x = fg_xz[linfun.diff_idxs, linfun.diff_idxs], f_z = fg_xz[linfun.diff_idxs, linfun.alge_idxs], @@ -794,7 +806,7 @@ function CommonSolve.solve(prob::LinearizationProblem; allow_input_derivatives = if !iszero(Bs) if !allow_input_derivatives der_inds = findall(vec(any(!=(0), Bs, dims = 1))) - error("Input derivatives appeared in expressions (-g_z\\g_u != 0), the following inputs appeared differentiated: $(inputs(prob.f.prob.f.sys)[der_inds]). Call `linearize` with keyword argument `allow_input_derivatives = true` to allow this and have the returned `B` matrix be of double width ($(2nu)), where the last $nu inputs are the derivatives of the first $nu inputs.") + error("Input derivatives appeared in expressions (-g_z\\g_u != 0), the following inputs appeared differentiated: $(prob.f.inputs[der_inds]). Call `linearize` with keyword argument `allow_input_derivatives = true` to allow this and have the returned `B` matrix be of double width ($(2nu)), where the last $nu inputs are the derivatives of the first $nu inputs.") end B = [B [zeros(nx, nu); Bs]] D = [D zeros(ny, nu)] diff --git a/test/linearize.jl b/test/linearize.jl index 33b56fd18bf..16e44d345e6 100644 --- a/test/linearize.jl +++ b/test/linearize.jl @@ -447,3 +447,22 @@ end @test m1.C ≈ m2.C @test m1.D ≈ m2.D end + +@testset "Linearization without inputs" begin + @variables x(t) y(t) + @named sys0 = System([D(x) ~ -x^2, y ~ 2x], t) + for autodiff in (AutoForwardDiff(), AutoFiniteDiff()) + lsys, _ = linearize(sys0, [], [y]; op = Dict(x => 1.0), autodiff) + @test lsys.A ≈ [-2.0;;] + @test size(lsys.B) == (1, 0) + @test lsys.C == [2.0;;] + @test size(lsys.D) == (1, 0) + end + # A model without unknowns and inputs + @named sys1 = System([y ~ 2.0 + 0 * t], t) + lsys, _ = linearize(sys1, [], [y]) + @test size(lsys.A) == (0, 0) + @test size(lsys.B) == (0, 0) + @test size(lsys.C) == (1, 0) + @test size(lsys.D) == (1, 0) +end