add find_similarity_transform (#1000)

* add `find_similarity_transform`

this was previously implemented in ControlSystemIdentification

* export

* fix

* up docstring

* up docstring

Author: baggepinnen

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.16.1"
+version = "1.17.0"
 
 [deps]
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"

lib/ControlSystemsBase/src/ControlSystemsBase.jl

@@ -44,6 +44,7 @@ export  LTISystem,
         balreal,
         baltrunc,
         similarity_transform,
+        find_similarity_transform,
         time_scale,
         innovation_form,
         observer_predictor,

lib/ControlSystemsBase/src/matrix_comps.jl

@@ -710,7 +710,7 @@ D̃ = D
 ```
 
 If `unitary=true`, `T` is assumed unitary and the matrix adjoint is used instead of the inverse.
-See also [`balance_statespace`](@ref).
+See also [`balance_statespace`](@ref), [`find_similarity_transform`](@ref).
 """
 function similarity_transform(sys::ST, T; unitary=false) where ST <: AbstractStateSpace
     if unitary
@@ -726,6 +726,44 @@ function similarity_transform(sys::ST, T; unitary=false) where ST <: AbstractSta
     ST(A,B,C,D,sys.timeevol)
 end
 
+"""
+    find_similarity_transform(sys1, sys2, method = :obsv)
+
+Find T such that `similarity_transform(sys1, T) == sys2`
+
+Ref: Minimal state-space realization in linear system theory: an overview, B. De Schutter
+
+If `method == :obsv`, the observability matrices of `sys1` and `sys2` are used to find `T`, whereas `method == :ctrb` uses the controllability matrices.
+
+```jldoctest
+julia> using ControlSystemsBase
+
+julia> T = randn(3,3);
+
+julia> sys1 = ssrand(1,1,3);
+
+julia> sys2 = similarity_transform(sys1, T);
+
+julia> T2 = find_similarity_transform(sys1, sys2);
+
+julia> T2 ≈ T
+true
+```
+"""
+function find_similarity_transform(sys1, sys2, method = :obsv)
+    if method === :obsv
+        O1 = obsv(sys1)
+        O2 = obsv(sys2)
+        return O1\O2
+    elseif method === :ctrb
+        C1 = ctrb(sys1)
+        C2 = ctrb(sys2)
+        return C1/C2
+    else
+        error("Unknown method $method")
+    end
+end
+
 """
     time_scale(sys::AbstractStateSpace{Continuous}, a; balanced = false)
     time_scale(G::TransferFunction{Continuous},     a; balanced = true)

lib/ControlSystemsBase/test/test_matrix_comps.jl

@@ -39,6 +39,19 @@ sysdb, _ = balance_statespace(sysd)
 
 @test ControlSystemsBase.balance_transform(A,B,C) ≈ ControlSystemsBase.balance_transform(sys)
 
+@testset "similarity transform" begin
+    @info "Testing similarity transform"
+    T = randn(3,3)
+    sys1 = ssrand(1,1,3)
+    sys2 = ControlSystemsBase.similarity_transform(sys1, T)
+    T2 = find_similarity_transform(sys1, sys2)
+    @test T2 ≈ T atol=1e-8
+
+    T3 = find_similarity_transform(sys1, sys2, :ctrb)
+    @test T3 ≈ T atol=1e-8
+
+end
+
 W = [1 0; 0 1]
 @test covar(sys, W) ≈ [0.002560975609756 0.002439024390244; 0.002439024390244 0.002560975609756]
 D2 = [1 0; 0 1]