Fix linearize without inputs and name differentiated inputs from the linearization function

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)]

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