Add function for stable/unstable decomposition (#947)

* Add function for stable/unstable decomposition

* fix another DSP deprecation

* fix

* test tweak

Author: baggepinnen

docs/src/man/creating_systems.md

@@ -497,9 +497,10 @@ Filters can be designed using [DSP.jl](https://docs.juliadsp.org/stable/filters/
 using DSP, ControlSystemsBase, Plots
 
 fs = 100
-df = digitalfilter(Bandpass(5, 10; fs), Butterworth(2))
+df = digitalfilter(Bandpass(5, 10), Butterworth(2); fs)
 G = tf(df, 1/fs) # Sample time must be provided in the conversion to get the correct frequency scale in the Bode plot
 bodeplot(G, xscale=:identity, yscale=:identity, hz=true)
+vline!([5 10], l=(:black, :dash), label="Band-pass limits", sp=1)
 ```
 
 See also

lib/ControlSystemsBase/Project.toml

@@ -2,7 +2,7 @@ name = "ControlSystemsBase"
 uuid = "aaaaaaaa-a6ca-5380-bf3e-84a91bcd477e"
 authors = ["Dept. Automatic Control, Lund University"]
 repo = "https://github.com/JuliaControl/ControlSystems.jl.git"
-version = "1.11.3"
+version = "1.12.0"
 
 [deps]
 DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2"

lib/ControlSystemsBase/src/ControlSystemsBase.jl

@@ -22,6 +22,7 @@ export  LTISystem,
         # Linear Algebra
         balance,
         balance_statespace,
+        stab_unstab,
         are,
         lqr,
         kalman,

lib/ControlSystemsBase/src/analysis.jl

@@ -1,7 +1,10 @@
 """
     poles(sys)
 
-Compute the poles of system `sys`."""
+Compute the poles of system `sys`.
+
+Note: Poles with multiplicity `n > 1` may suffer numerical inaccuracies on the order `eps(numeric_type(sys))^(1/n)`, i.e., a double pole in a system with `Float64` coefficients may be computed with an error of about `√(eps(Float64)) ≈ 1.5e-8`.
+"""
 poles(sys::AbstractStateSpace) = eigvalsnosort(sys.A)
 
 # Seems to have a lot of rounding problems if we run the full thing with sisorational,

lib/ControlSystemsBase/src/matrix_comps.jl

@@ -671,6 +671,27 @@ function baltrunc(sys::ST; atol = sqrt(eps(numeric_type(sys))), rtol = 1e-3, n =
     return ST(A,B,C,D,sys.timeevol), S, T
 end
 
+"""
+    stab, unstab, sep = stab_unstab(sys; kwargs...)
+
+Decompose `sys` into `sys = stab + unstab` where `stab` contains all stable poles and `unstab` contains unstable poles. 
+
+`0 ≤ sep ≤ 1` is the estimated separation between the stable and unstable spectra.
+
+The docstring of `MatrixPencils.ssblkdiag`, reproduced below, provides more information on the keyword arguments:
+$(@doc(MatrixPencils.ssblkdiag))
+"""
+function stab_unstab(sys::AbstractStateSpace; kwargs...)
+    stable_unstable = true
+    disc = isdiscrete(sys)
+    A, B, C, _, _, blkdims, sep = MatrixPencils.ssblkdiag(sys.A, sys.B, sys.C; disc, stable_unstable, withQ = false, withZ = false, kwargs...)
+    n1 = blkdims[1];
+    i1 = 1:n1; i2 = n1+1:sys.nx 
+    return (; stab=ss(A[i1,i1], B[i1,:], C[:,i1], sys.D, timeevol(sys)), 
+            unstab=ss(A[i2,i2], B[i2,:], C[:,i2], 0, timeevol(sys)),
+            sep)
+end
+
 """
     syst = similarity_transform(sys, T; unitary=false)
 Perform a similarity transform `T : Tx̃ = x` on `sys` such that

lib/ControlSystemsBase/test/test_matrix_comps.jl

@@ -88,6 +88,20 @@ syst = similarity_transform(sys, Tr)
 @test sys.C*Tr ≈ syst.C
 
 
+## stab_unstab
+sys = ssrand(2,3,40, stable=false)
+stab, unstab = stab_unstab(sys)
+@test all(real(poles(stab)) .< 0)
+@test all(real(poles(unstab)) .>= 0)
+@test linfnorm(stab + unstab - sys)[1] < 1e-5
+
+sys = ssrand(2,3,40, stable=false, Ts=1)
+stab, unstab = stab_unstab(sys)
+@test all(abs.(poles(stab)) .< 1)
+@test all(abs.(poles(unstab)) .>= 1)
+@test linfnorm(stab + unstab - sys)[1] < 1e-5
+
+
 sys = ss([1 0.1; 0 1], ones(2), [1. 0], 0)
 sysi = ControlSystemsBase.innovation_form(sys, I, I)
 @test sysi.A ≈ sysi.A