handle systems without states in linearize

Author: baggepinnen

src/systems/abstractsystem.jl

@@ -1027,13 +1027,22 @@ function linearization_function(sys::AbstractSystem, inputs, outputs; simplify =
     lin_fun = let fun = fun,
         h = ModelingToolkit.build_explicit_observed_function(sys, outputs)
 
-        (u, p, t) -> begin
-            uf = SciMLBase.UJacobianWrapper(fun, t, p)
-            fg_xz = ForwardDiff.jacobian(uf, u)
-            pf = SciMLBase.ParamJacobianWrapper(fun, t, u)
-            # TODO: this is very inefficient, p contains all parameters of the system
-            fg_u = ForwardDiff.jacobian(pf, p)[:, input_idxs]
-            h_xz = ForwardDiff.jacobian(xz -> h(xz, p, t), u)
+        function (u, p, t)
+            if u !== nothing # Handle systems without states
+                length(sts) == length(u) ||
+                    error("Number of state variables does not match the number of input states")
+                uf = SciMLBase.UJacobianWrapper(fun, t, p)
+                fg_xz = ForwardDiff.jacobian(uf, u)
+                h_xz = ForwardDiff.jacobian(xz -> h(xz, p, t), u)
+                pf = SciMLBase.ParamJacobianWrapper(fun, t, u)
+                # TODO: this is very inefficient, p contains all parameters of the system
+                fg_u = ForwardDiff.jacobian(pf, p)[:, input_idxs]
+            else
+                length(sts) == 0 ||
+                    error("Number of state variables does not match the number of input states")
+                fg_xz = zeros(0,0)
+                h_xz = fg_u = zeros(0, length(inputs))
+            end
             h_u = ForwardDiff.jacobian(p -> h(u, p, t), p)[:, input_idxs]
             (f_x = fg_xz[diff_idxs, diff_idxs],
              f_z = fg_xz[diff_idxs, alge_idxs],

test/linearize.jl

@@ -103,3 +103,26 @@ lsys = ModelingToolkit.reorder_states(lsys, states(ssys), reverse(desired_order)
 @test lsys.B == [10 -10; 2 -2]
 @test lsys.C == [-4000 400]
 @test lsys.D == [4400 -4400]
+
+## Test operating points
+
+# The saturation has no dynamics
+function saturation(; y_max, y_min=y_max > 0 ? -y_max : -Inf, name)
+    @variables u(t)=0 y(t)=0
+    @parameters y_max=y_max y_min=y_min
+    ie = ModelingToolkit.IfElse.ifelse
+    eqs = [
+        # The equation below is equivalent to y ~ clamp(u, y_min, y_max)
+        y ~ ie(u > y_max, y_max, ie( (y_min < u) & (u < y_max), u, y_min))
+    ]
+    ODESystem(eqs, t, name=name)
+end
+
+@named sat = saturation(; y_max=1)
+# inside the linear region, the function is identity
+@unpack u,y = sat
+lsys, ssys = linearize(sat, [u], [y])
+@test isempty(lsys.A) # there are no differential variables in this system
+@test isempty(lsys.B)
+@test isempty(lsys.C)
+@test lsys.D[] == 1