add option to simplify causality in a different way

Author: baggepinnen

src/ode_system.jl

@@ -154,19 +154,22 @@ end
 
 
 """
-    RobustAndOptimalControl.named_ss(sys::ModelingToolkit.AbstractTimeDependentSystem, inputs, outputs; kwargs...)
+    RobustAndOptimalControl.named_ss(sys::ModelingToolkit.AbstractTimeDependentSystem, inputs, outputs; descriptor=true, kwargs...)
 
 Convert an `ODESystem` to a `NamedStateSpace` using linearization. `inputs, outputs` are vectors of variables determining the inputs and outputs respectively. See docstring of `ModelingToolkit.linearize` for more info on `kwargs`.
 
-This method automatically converts systems that MTK has failed to produce a proper form for into a proper linear statespace system. Learn more about how that is done here:
+If `descriptor = true` (default), this method automatically converts systems that MTK has failed to produce a proper form for into a proper linear statespace system using the method described here:
 https://juliacontrol.github.io/ControlSystemsMTK.jl/dev/#Internals:-Transformation-of-non-proper-models-to-proper-statespace-form
+If `descriptor = false`, the system is instead converted to a statespace realization using `sys[:,uinds] + sys[:,duinds]*tf('s')`, which tends to result in a larger realization on which the user may want to call `minreal(sys, tol)` with a carefully selected tolerance.
+
 
 See also [`ModelingToolkit.linearize`](@ref) which is the lower-level function called internally. The functions [`get_named_sensitivity`](@ref), [`get_named_comp_sensitivity`](@ref), [`get_named_looptransfer`](@ref) similarily provide convenient ways to compute sensitivity functions while retaining signal names in the same way as `named_ss`. The corresponding lower-level functions `get_sensitivity`, `get_comp_sensitivity` and `get_looptransfer` are available in ModelingToolkitStandardLibrary.Blocks and are documented in [MTKstdlib: Linear analysis](https://docs.sciml.ai/ModelingToolkitStandardLibrary/stable/API/linear_analysis/).
 """
 function RobustAndOptimalControl.named_ss(
     sys::ModelingToolkit.AbstractTimeDependentSystem,
     inputs,
     outputs;
+    descriptor = true,
     kwargs...,
 )
 
@@ -206,7 +209,7 @@ function RobustAndOptimalControl.named_ss(
         # This indicates that input derivatives are present
         duinds = findall(any(!iszero, eachcol(matrices.B[:, nu+1:end]))) .+ nu
         u2du = (1:nu) .=> duinds # This maps inputs to their derivatives
-        lsys = causal_simplification(matrices, u2du)
+        lsys = causal_simplification(matrices, u2du; descriptor)
     else
         lsys = ss(matrices...)
     end
@@ -219,26 +222,33 @@ function RobustAndOptimalControl.named_ss(
     )
 end
 
-function causal_simplification(sys, u2duinds::Vector{Pair{Int, Int}})
+function causal_simplification(sys, u2duinds::Vector{Pair{Int, Int}}; descriptor=true)
+    fm(x) = convert(Matrix{Float64}, x)
     nx = size(sys.A, 1)
     ny = size(sys.C, 1)
     ndu = length(u2duinds)
     nu = size(sys.B, 2) - ndu
     u_with_du_inds = first.(u2duinds)
     duinds = last.(u2duinds)
-    B̄ = sys.B[:, duinds]
     B = sys.B[:, 1:nu]
-    Iu = u_with_du_inds .== (1:nu)'
-    E = [I(nx) -B̄; zeros(ndu, nx+ndu)]
-
-    fm(x) = convert(Matrix{Float64}, x)
+    B̄ = sys.B[:, duinds]
+    D = sys.D[:, 1:nu]
+    D̄ = sys.D[:, duinds]
+    iszero(fm(D̄)) || error("Nonzero feedthrough matrix from input derivative not supported")
+    if descriptor
+        Iu = u_with_du_inds .== (1:nu)'
+        E = [I(nx) -B̄; zeros(ndu, nx+ndu)]
 
-    Ae = cat(sys.A, -I(ndu), dims=(1,2))
-    Be = [B; Iu]
-    Ce = [fm(sys.C) zeros(ny, ndu)]
-    De = fm(sys.D[:, 1:nu])
-    dsys = dss(Ae, E, Be, Ce, De)
-    ss(RobustAndOptimalControl.DescriptorSystems.dss2ss(dsys)[1])
+        Ae = cat(sys.A, -I(ndu), dims=(1,2))
+        # Ae[1:nx, nx+1:end] .= B
+        Be = [B; Iu]
+        Ce = [fm(sys.C) zeros(ny, ndu)]
+        De = fm(D)
+        dsys = dss(Ae, E, Be, Ce, De)
+        return ss(RobustAndOptimalControl.DescriptorSystems.dss2ss(dsys)[1])
+    else
+        ss(sys.A, B, sys.C, D) + ss(sys.A, B̄, sys.C, D̄)*tf('s')
+    end
 end
 
 for f in [:sensitivity, :comp_sensitivity, :looptransfer]

test/test_ODESystem.jl

@@ -295,7 +295,7 @@ using ModelingToolkitStandardLibrary.Mechanical.TranslationalPosition
 using Test
 
 using ControlSystemsMTK
-using ControlSystemsMTK.ControlSystemsBase: sminreal, minreal, poles
+using ControlSystemsMTK.ControlSystemsBase: sminreal, minreal, poles, ss, tf
 connect = ModelingToolkit.connect
 
 @independent_variables t
@@ -348,3 +348,38 @@ op2[cart.f] = 0
 @test G.nx == 4
 @test G.nu == length(lin_inputs)
 @test G.ny == length(lin_outputs)
+
+## Test difficult `named_ss` simplification
+using ControlSystemsMTK, ControlSystemsBase, RobustAndOptimalControl
+lsys = (A = [0.0 2.778983834717109e8 1.4122312296634873e6 0.0; 0.0 0.0 0.0 0.037848975765016724; 0.0 24.837541148074962 0.12622006230897712 0.0; -0.0 -4.620724819774693 -0.023481719514324866 -0.6841991610512456], B = [-5.042589978197361e8 0.0; -0.0 0.0; -45.068824982602656 -0.0; 8.384511049369085 54.98555939873381], C = [0.0 0.0 0.954929658551372 0.0], D = [0.0 0.0])
+
+# lsys = (A = [-0.0075449237853825925 1.6716817118020731e-6 0.0; 1864.7356343162514 -0.4131578457122937 0.0; 0.011864343330426718 -2.6287085638214332e-6 0.0], B = [0.0 0.0; 0.0 52566.418015009294; 0.0 0.3284546792274811], C = [1.4683007399899215e8 0.0 0.0], D = [-9.157636303058283e7 0.0])
+
+G = ControlSystemsMTK.causal_simplification(lsys, [1=>2])
+G2 = ControlSystemsMTK.causal_simplification(lsys, [1=>2], descriptor=false)
+G2 = minreal(G2, 1e-12)
+
+@test dcgain(G, 1e-5)[] ≈ dcgain(G2, 1e-5)[] rtol=1e-3
+@test freqresp(G, 1)[] ≈ freqresp(G2, 1)[]
+@test freqresp(G, 10)[] ≈ freqresp(G2, 10)[]
+
+z = 0.462726166562343204837317130554462562
+
+@test minimum(abs, tzeros(G) .- z) < sqrt(eps())
+@test minimum(abs, tzeros(G2) .- z) < sqrt(eps())
+
+# using GenericSchur
+
+Gb = balance_statespace(G)[1]
+Gb = minreal(Gb, 1e-8)
+@test Gb.nx == 2
+@test minimum(abs, tzeros(Gb) .- z) < sqrt(eps())
+
+w = exp10.(LinRange(-12, 2, 2000))
+ControlSystemsBase.bodeplot([G, G2, minreal(G, 1e-8)], w)
+
+
+##
+
+# S = schur(A,B)
+# V = eigvecs(S)