Merge pull request #2003 from SciML/fb/lind_dummy_der

option to set all dummy derivatives to zero in `linearize`

Author: YingboMa

Project.toml

@@ -87,6 +87,7 @@ julia = "1.6"
 [extras]
 AmplNLWriter = "7c4d4715-977e-5154-bfe0-e096adeac482"
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+ControlSystemsMTK = "687d7614-c7e5-45fc-bfc3-9ee385575c88"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9"
 Ipopt_jll = "9cc047cb-c261-5740-88fc-0cf96f7bdcc7"
@@ -107,4 +108,4 @@ Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["AmplNLWriter", "BenchmarkTools", "NonlinearSolve", "ForwardDiff", "Ipopt", "Ipopt_jll", "ModelingToolkitStandardLibrary", "Optimization", "OptimizationOptimJL", "OptimizationMOI", "OrdinaryDiffEq", "Random", "ReferenceTests", "SafeTestsets", "StableRNGs", "Statistics", "SteadyStateDiffEq", "Test", "StochasticDiffEq", "Sundials"]
+test = ["AmplNLWriter", "BenchmarkTools", "ControlSystemsMTK", "NonlinearSolve", "ForwardDiff", "Ipopt", "Ipopt_jll", "ModelingToolkitStandardLibrary", "Optimization", "OptimizationOptimJL", "OptimizationMOI", "OrdinaryDiffEq", "Random", "ReferenceTests", "SafeTestsets", "StableRNGs", "Statistics", "SteadyStateDiffEq", "Test", "StochasticDiffEq", "Sundials"]

src/systems/abstractsystem.jl

@@ -1362,8 +1362,13 @@ end
 
 function linearize(sys, inputs, outputs; op = Dict(), t = 0.0,
                    allow_input_derivatives = false,
+                   zero_dummy_der = false,
                    kwargs...)
     lin_fun, ssys = linearization_function(sys, inputs, outputs; kwargs...)
+    if zero_dummy_der
+        dummyder = setdiff(states(ssys), states(sys))
+        op = merge(op, Dict(x => 0.0 for x in dummyder))
+    end
     linearize(ssys, lin_fun; op, t, allow_input_derivatives), ssys
 end
 

test/linearize.jl

@@ -157,3 +157,53 @@ lsys, ssys = linearize(sat, [u], [y]; op = Dict(u => 2))
 @test isempty(lsys.B)
 @test isempty(lsys.C)
 @test lsys.D[] == 0
+
+## Test that dummy_derivatives can be set to zero
+if VERSION >= v"1.8"
+    # The call to Link(; m = 0.2, l = 10, I = 1, g = -9.807) hangs forever on Julia v1.6
+    using LinearAlgebra
+    using ModelingToolkit
+    using ModelingToolkitStandardLibrary
+    using ModelingToolkitStandardLibrary.Blocks
+    using ModelingToolkitStandardLibrary.Mechanical.MultiBody2D
+    using ModelingToolkitStandardLibrary.Mechanical.Translational
+
+    using ControlSystemsMTK
+    using ControlSystemsMTK.ControlSystemsBase: sminreal, minreal, poles
+    connect = ModelingToolkit.connect
+
+    @parameters t
+    D = Differential(t)
+
+    @named link1 = Link(; m = 0.2, l = 10, I = 1, g = -9.807)
+    @named cart = Translational.Mass(; m = 1, s_0 = 0)
+    @named fixed = Fixed()
+    @named force = Force()
+
+    eqs = [connect(link1.TX1, cart.flange)
+           connect(cart.flange, force.flange)
+           connect(link1.TY1, fixed.flange)]
+
+    @show @named model = ODESystem(eqs, t, [], []; systems = [link1, cart, force, fixed])
+    def = ModelingToolkit.defaults(model)
+    def[link1.y1] = 0
+    def[link1.x1] = 10
+    def[link1.A] = -pi / 2
+    def[link1.dA] = 0
+    def[cart.s] = 0
+    def[force.flange.v] = 0
+    lin_outputs = [cart.s, cart.v, link1.A, link1.dA]
+    lin_inputs = [force.f.u]
+
+    @info "named_ss"
+    G = named_ss(model, lin_inputs, lin_outputs, allow_symbolic = true, op = def,
+                 allow_input_derivatives = true, zero_dummy_der = true)
+    G = sminreal(G)
+    @info "minreal"
+    G = minreal(G)
+    @info "poles"
+    ps = poles(G)
+
+    @test minimum(abs, ps) < 1e-6
+    @test minimum(abs, complex(0, 1.3777260367206716) .- ps) < 1e-10
+end