From 208ca926894e6e2e71b2426b4eb00eb4afb84564 Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Fri, 14 Feb 2025 07:54:08 +0100
Subject: [PATCH 1/4] improve tf2ss conversion

by using MatrixPencils, which offers possibility of constructing a minimal realization. This option is _much_ faster than the old approach `minreal(ss(P))` for large systems
---
 docs/src/man/creating_systems.md              |   5 +-
 .../src/types/StateSpace.jl                   |   5 +-
 .../src/types/conversion.jl                   | 133 ++++++++----------
 3 files changed, 66 insertions(+), 77 deletions(-)

diff --git a/docs/src/man/creating_systems.md b/docs/src/man/creating_systems.md
index 2767b930b..825dbd368 100644
--- a/docs/src/man/creating_systems.md
+++ b/docs/src/man/creating_systems.md
@@ -102,6 +102,7 @@ in discrete time, is created using
 ```julia
 ss(A,B,C,D)    # Continuous-time system
 ss(A,B,C,D,Ts) # Discrete-time system
+ss(P; balance=true, minimal=false) # Convert transfer function P to state space
 ```
 and they behave similarly to transfer functions.
 
@@ -118,7 +119,7 @@ HeteroStateSpace(sys, to_static)
 ```
 Notice the different matrix types used.
 
-To associate **names** with states, inputs and outputs, see [`named_ss`](https://juliacontrol.github.io/RobustAndOptimalControl.jl/dev/#Named-systems) from RobustAndOptimalControl.jl.
+To associate **names** with state variables, inputs and outputs, see [`named_ss`](https://juliacontrol.github.io/RobustAndOptimalControl.jl/dev/#Named-systems) from RobustAndOptimalControl.jl.
 
 
 ## Converting between types
@@ -137,6 +138,8 @@ Sample Time: 0.1 (seconds)
 Discrete-time transfer function model
 ```
 
+When a transfer function `P` is converted to a state-space model using `ss(P; balance=true, minimal=false)`, the user may choose whether to balance the state-space model (default=true) and/or to return a minimal realization (default=false).
+
 
 ## Converting between continuous and discrete time
 A continuous-time system represents differential equations or a transfer function in the [Laplace domain](https://en.wikipedia.org/wiki/Laplace_transform), while a discrete-time system represents difference equations or a transfer function in the [Z-domain](https://en.wikipedia.org/wiki/Z-transform).
diff --git a/lib/ControlSystemsBase/src/types/StateSpace.jl b/lib/ControlSystemsBase/src/types/StateSpace.jl
index 2cd66eb09..23b55dd42 100644
--- a/lib/ControlSystemsBase/src/types/StateSpace.jl
+++ b/lib/ControlSystemsBase/src/types/StateSpace.jl
@@ -128,6 +128,7 @@ StateSpace(sys::LTISystem; kwargs...) = convert(StateSpace, sys; kwargs...)
 """
     sys = ss(A, B, C, D)      # Continuous
     sys = ss(A, B, C, D, Ts)  # Discrete
+    sys = ss(P::TransferFunction; balance=true, minimal=false) # Convert TransferFunction to StateSpace
 
 Create a state-space model `sys::StateSpace{TE, T}`
 with matrix element type `T` and TE is `Continuous` or `<:Discrete`.
@@ -138,7 +139,9 @@ Otherwise, this is a discrete-time model with sampling period `Ts`.
 `D` may be specified as `0` in which case a zero matrix of appropriate size is constructed automatically. 
 `sys = ss(D [, Ts])` specifies a static gain matrix `D`.
 
-To associate names with states, inputs and outputs, see [`named_ss`](https://juliacontrol.github.io/RobustAndOptimalControl.jl/dev/#Named-systems).
+When a transfer function `P` is converted to a state-space model, the user may choose whether to balance the state-space model (default=true) and/or to return a minimal realization (default=false). This conversion has more keyword argument options, see `?ControlSystemsBase.MatrixPencils.rm2ls` for additional details.
+
+To associate names with state variables, inputs and outputs, see [`named_ss`](https://juliacontrol.github.io/RobustAndOptimalControl.jl/dev/#Named-systems).
 """
 ss(args...;kwargs...) = StateSpace(args...;kwargs...)
 
diff --git a/lib/ControlSystemsBase/src/types/conversion.jl b/lib/ControlSystemsBase/src/types/conversion.jl
index 8c8b9e39e..b3777659a 100644
--- a/lib/ControlSystemsBase/src/types/conversion.jl
+++ b/lib/ControlSystemsBase/src/types/conversion.jl
@@ -1,16 +1,3 @@
-# Base.convert(::Type{<:SisoTf}, b::Real) = Base.convert(SisoRational, b)
-# Base.convert{T<:Real}(::Type{<:SisoZpk}, b::T) = SisoZpk(T[], T[], b)
-# Base.convert{T<:Real}(::Type{<:SisoRational}, b::T) = SisoRational([b], [one(T)])
-# Base.convert{T1}(::Type{SisoRational{Vector{T1}}}, t::SisoRational) =  SisoRational(Polynomial(T1.(t.num.coeffs$1),Polynomial(T1.(t.den.coeffs$1))
-# Base.convert(::Type{<:StateSpace}, t::Real) = ss(t)
-#
-
-#
-# function Base.convert{T<:AbstractMatrix{<:Number}}(::Type{StateSpace{T}}, s::StateSpace)
-#     AT = promote_type(T, arraytype(s))
-#     StateSpace{AT}(AT(s.A),AT(s.B),AT(s.C),AT(s.D), s.timeevol, s.statenames, s.inputnames, s.outputnames)
-# end
-
 # TODO Fix these to use proper constructors
 # NOTE: no real need to convert numbers to transfer functions, have addition methods..
 # How to convert a number to either Continuous or Discrete transfer function
@@ -42,19 +29,7 @@ Base.convert(::Type{TransferFunction{TE,SisoZpk{T,TR}}}, d::Number) where {TE, T
 
 Base.convert(::Type{StateSpace{TE,T}}, d::Number) where {TE, T} =
     convert(StateSpace{TE,T}, fill(d, (1,1)))
-#
-# Base.convert(::Type{TransferFunction{Continuous,<:SisoRational{T}}}, b::Number) where {T} = tf(T(b), Continuous())
-# Base.convert(::Type{TransferFunction{Continuous,<:SisoZpk{T, TR}}}, b::Number) where {T, TR} = zpk(T(b), Continuous())
-# Base.convert(::Type{StateSpace{Continuous,T, MT}}, b::Number) where {T, MT} = convert(StateSpace{Continuous,T, MT}, fill(b, (1,1)))
-
-# function Base.convert{T<:Real,S<:TransferFunction}(::Type{S}, b::VecOrMat{T})
-#     r = Matrix{S}(size(b,2),1)
-#     for j=1:size(b,2)
-#         r[j] = vcat(map(k->convert(S,k),b[:,j])...)
-#     end
-#     hcat(r...)
-# end
-#
+
 
 function convert(::Type{TransferFunction{TE,S}}, G::TransferFunction) where {TE,S}
     Gnew_matrix = convert.(S, G.matrix)
@@ -92,77 +67,85 @@ struct ImproperException <: Exception end
 Base.showerror(io::IO, e::ImproperException) = print(io, "System is improper, a state-space representation is impossible")
 
 # Note: balancing is only applied by default for floating point types, integer systems are not balanced since that would change the type. 
-function Base.convert(::Type{StateSpace{TE,T}}, G::TransferFunction; balance=!(T <: Union{Integer, Rational, ForwardDiff.Dual})) where {TE,T<:Number}
-    if !isproper(G)
-        throw(ImproperException())
-    end
+function Base.convert(::Type{StateSpace{TE,T}}, G::TransferFunction; balance=!(T <: Union{Integer, Rational, ForwardDiff.Dual}), minimal=false, contr=minimal, obs=minimal, noseig=minimal, kwargs...) where {TE,T<:Number}
 
-    ny, nu = size(G)
 
-    # A, B, C, D matrices for each element of the transfer function matrix
-    abcd_vec = [siso_tf_to_ss(T, g) for g in G.matrix[:]]
+    NUM = numpoly(G)
+    DEN = denpoly(G)
+    (A, E, B, C, D, blkdims) = MatrixPencils.rm2ls(NUM, DEN; contr, obs, noseig, minimal, kwargs...)
+    E == I || throw(ImproperException())
 
-    # Number of states for each transfer function element realization
-    nvec = [size(abcd[1], 1) for abcd in abcd_vec]
-    ntot = sum(nvec)
 
-    A = zeros(T, (ntot, ntot))
-    B = zeros(T, (ntot, nu))
-    C = zeros(T, (ny, ntot))
-    D = zeros(T, (ny, nu))
+    # if !isproper(G)
+    #     throw(ImproperException())
+    # end
 
-    inds = -1:0
-    for j=1:nu
-        for i=1:ny
-            k = (j-1)*ny + i
+    # ny, nu = size(G)
 
-            # states corresponding to the transfer function element (i,j)
-            inds = (inds.stop+1):(inds.stop+nvec[k])
+    # # A, B, C, D matrices for each element of the transfer function matrix
+    # abcd_vec = [siso_tf_to_ss(T, g) for g in G.matrix[:]]
 
-            A[inds,inds], B[inds,j:j], C[i:i,inds], D[i:i,j:j] = abcd_vec[k]
-        end
-    end
+    # # Number of states for each transfer function element realization
+    # nvec = [size(abcd[1], 1) for abcd in abcd_vec]
+    # ntot = sum(nvec)
+
+    # A = zeros(T, (ntot, ntot))
+    # B = zeros(T, (ntot, nu))
+    # C = zeros(T, (ny, ntot))
+    # D = zeros(T, (ny, nu))
+
+    # inds = -1:0
+    # for j=1:nu
+    #     for i=1:ny
+    #         k = (j-1)*ny + i
+
+    #         # states corresponding to the transfer function element (i,j)
+    #         inds = (inds.stop+1):(inds.stop+nvec[k])
+
+    #         A[inds,inds], B[inds,j:j], C[i:i,inds], D[i:i,j:j] = abcd_vec[k]
+    #     end
+    # end
     if balance
         A, B, C = balance_statespace(A, B, C)
     end
     return StateSpace{TE,T}(A, B, C, D, TE(G.timeevol))
 end
 
-siso_tf_to_ss(T::Type, f::SisoTf) = siso_tf_to_ss(T, convert(SisoRational, f))
+# siso_tf_to_ss(T::Type, f::SisoTf) = siso_tf_to_ss(T, convert(SisoRational, f))
 
 # Conversion to statespace on controllable canonical form
-function siso_tf_to_ss(T::Type, f::SisoRational)
+# function siso_tf_to_ss(T::Type, f::SisoRational)
 
-    num0, den0 = numvec(f), denvec(f)
-    # Normalize the numerator and denominator to allow realization of transfer functions
-    # that are proper, but not strictly proper
-    num = num0 ./ den0[1]
-    den = den0 ./ den0[1]
+#     num0, den0 = numvec(f), denvec(f)
+#     # Normalize the numerator and denominator to allow realization of transfer functions
+#     # that are proper, but not strictly proper
+#     num = num0 ./ den0[1]
+#     den = den0 ./ den0[1]
 
-    N = length(den) - 1 # The order of the rational function f
+#     N = length(den) - 1 # The order of the rational function f
 
-    # Get numerator coefficient of the same order as the denominator
-    bN = length(num) == N+1 ? num[1] : zero(eltype(num))
+#     # Get numerator coefficient of the same order as the denominator
+#     bN = length(num) == N+1 ? num[1] : zero(eltype(num))
 
-    @views if N == 0 #|| num == zero(Polynomial{T})
-        A = zeros(T, 0, 0)
-        B = zeros(T, 0, 1)
-        C = zeros(T, 1, 0)
-    else
-        A = diagm(1 => ones(T, N-1))
-        A[end, :] .= .-reverse(den)[1:end-1]
+#     @views if N == 0 #|| num == zero(Polynomial{T})
+#         A = zeros(T, 0, 0)
+#         B = zeros(T, 0, 1)
+#         C = zeros(T, 1, 0)
+#     else
+#         A = diagm(1 => ones(T, N-1))
+#         A[end, :] .= .-reverse(den)[1:end-1]
 
-        B = zeros(T, N, 1)
-        B[end] = one(T)
+#         B = zeros(T, N, 1)
+#         B[end] = one(T)
 
-        C = zeros(T, 1, N)
-        C[1:min(N, length(num))] = reverse(num)[1:min(N, length(num))]
-        C[:] .-= bN .* reverse(den)[1:end-1] # Can index into polynomials at greater inddices than their length
-    end
-    D = fill(bN, 1, 1)
+#         C = zeros(T, 1, N)
+#         C[1:min(N, length(num))] = reverse(num)[1:min(N, length(num))]
+#         C[:] .-= bN .* reverse(den)[1:end-1] # Can index into polynomials at greater inddices than their length
+#     end
+#     D = fill(bN, 1, 1)
 
-    return A, B, C, D
-end
+#     return A, B, C, D
+# end
 
 """
     A, B, C, T = balance_statespace{S}(A::Matrix{S}, B::Matrix{S}, C::Matrix{S}, perm::Bool=false)

From 455e14b2b0ee84b1a9aa99f2d36d766a8fe27d2a Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Fri, 14 Feb 2025 08:18:17 +0100
Subject: [PATCH 2/4] make time responses of tf use minimal realization

---
 lib/ControlSystemsBase/src/timeresp.jl | 6 +++---
 test/test_timeresp.jl                  | 4 ++--
 2 files changed, 5 insertions(+), 5 deletions(-)

diff --git a/lib/ControlSystemsBase/src/timeresp.jl b/lib/ControlSystemsBase/src/timeresp.jl
index 3f8dc411c..3dcaef736 100644
--- a/lib/ControlSystemsBase/src/timeresp.jl
+++ b/lib/ControlSystemsBase/src/timeresp.jl
@@ -71,7 +71,7 @@ end
 
 Base.step(sys::LTISystem, tfinal::Real; kwargs...) = step(sys, _default_time_vector(sys, tfinal); kwargs...)
 Base.step(sys::LTISystem; kwargs...) = step(sys, _default_time_vector(sys); kwargs...)
-Base.step(sys::TransferFunction, t::AbstractVector; kwargs...) = step(ss(sys), t::AbstractVector; kwargs...)
+Base.step(sys::TransferFunction, t::AbstractVector; kwargs...) = step(ss(sys, balance=true, minimal=true), t::AbstractVector; kwargs...)
 
 """
     y, t, x = impulse(sys[, tfinal])
@@ -119,7 +119,7 @@ end
 
 impulse(sys::LTISystem, tfinal::Real; kwargs...) = impulse(sys, _default_time_vector(sys, tfinal); kwargs...)
 impulse(sys::LTISystem; kwargs...) = impulse(sys, _default_time_vector(sys); kwargs...)
-impulse(sys::TransferFunction, t::AbstractVector; kwargs...) = impulse(ss(sys), t; kwargs...)
+impulse(sys::TransferFunction, t::AbstractVector; kwargs...) = impulse(ss(sys, balance=true, minimal=true), t; kwargs...)
 
 """
     result = lsim(sys, u[, t]; x0, method])
@@ -301,7 +301,7 @@ function lsim(sys::AbstractStateSpace, u::Function, t::AbstractVector;
 end
 
 
-lsim(sys::TransferFunction, args...; kwargs...) = lsim(ss(sys), args...; kwargs...)
+lsim(sys::TransferFunction, args...; kwargs...) = lsim(ss(sys, balance=true, minimal=true), args...; kwargs...)
 
 
 """
diff --git a/test/test_timeresp.jl b/test/test_timeresp.jl
index 77314f4bc..ddf0e95c1 100644
--- a/test/test_timeresp.jl
+++ b/test/test_timeresp.jl
@@ -75,7 +75,7 @@ xreal[3,:,2] = exp.(-t).*t
 y, t, x = step(systf, t0, method=:zoh)
 yreal[1,:,1] = 1 .- exp.(-t)
 yreal[2,:,2] = -1 .+ exp.(-t) + 2*exp.(-t).*t
-@test y ≈ yreal atol=1e-14
+@test y ≈ yreal atol=1e-13
 #Step ss
 y, t, x = step(sysss, t, method=:zoh)
 @test y ≈ yreal atol=1e-13
@@ -88,7 +88,7 @@ xreal[3,:,2] = exp.(-t).*(-t .- 1) .+ 1
 y, t, x = impulse(systf, t, method=:zoh)
 yreal[1,:,1] = exp.(-t)
 yreal[2,:,2] = exp.(-t).*(1 .- 2*t)
-@test y ≈ yreal atol=1e-14
+@test y ≈ yreal atol=1e-13
 #Impulse ss
 y, t, x = impulse(1.0sysss, t, method=:zoh)
 @test y ≈ yreal atol=1e-13

From e69d2fcf6b76f1f40965266855c57942a7626a77 Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Fri, 14 Feb 2025 08:31:35 +0100
Subject: [PATCH 3/4] test differentiability of tf2ss conversion

---
 lib/ControlSystemsBase/src/timeresp.jl            |  6 +++---
 lib/ControlSystemsBase/test/test_implicit_diff.jl | 14 ++++++++++++++
 2 files changed, 17 insertions(+), 3 deletions(-)

diff --git a/lib/ControlSystemsBase/src/timeresp.jl b/lib/ControlSystemsBase/src/timeresp.jl
index 3dcaef736..e7d91ea1d 100644
--- a/lib/ControlSystemsBase/src/timeresp.jl
+++ b/lib/ControlSystemsBase/src/timeresp.jl
@@ -71,7 +71,7 @@ end
 
 Base.step(sys::LTISystem, tfinal::Real; kwargs...) = step(sys, _default_time_vector(sys, tfinal); kwargs...)
 Base.step(sys::LTISystem; kwargs...) = step(sys, _default_time_vector(sys); kwargs...)
-Base.step(sys::TransferFunction, t::AbstractVector; kwargs...) = step(ss(sys, balance=true, minimal=true), t::AbstractVector; kwargs...)
+Base.step(sys::TransferFunction, t::AbstractVector; kwargs...) = step(ss(sys, minimal=numeric_type(sys) isa BlasFloat), t::AbstractVector; kwargs...)
 
 """
     y, t, x = impulse(sys[, tfinal])
@@ -119,7 +119,7 @@ end
 
 impulse(sys::LTISystem, tfinal::Real; kwargs...) = impulse(sys, _default_time_vector(sys, tfinal); kwargs...)
 impulse(sys::LTISystem; kwargs...) = impulse(sys, _default_time_vector(sys); kwargs...)
-impulse(sys::TransferFunction, t::AbstractVector; kwargs...) = impulse(ss(sys, balance=true, minimal=true), t; kwargs...)
+impulse(sys::TransferFunction, t::AbstractVector; kwargs...) = impulse(ss(sys, minimal=numeric_type(sys) isa BlasFloat), t; kwargs...)
 
 """
     result = lsim(sys, u[, t]; x0, method])
@@ -301,7 +301,7 @@ function lsim(sys::AbstractStateSpace, u::Function, t::AbstractVector;
 end
 
 
-lsim(sys::TransferFunction, args...; kwargs...) = lsim(ss(sys, balance=true, minimal=true), args...; kwargs...)
+lsim(sys::TransferFunction, args...; kwargs...) = lsim(ss(sys, minimal=numeric_type(sys) isa BlasFloat), args...; kwargs...)
 
 
 """
diff --git a/lib/ControlSystemsBase/test/test_implicit_diff.jl b/lib/ControlSystemsBase/test/test_implicit_diff.jl
index 7e95f41cb..7a4cb513e 100644
--- a/lib/ControlSystemsBase/test/test_implicit_diff.jl
+++ b/lib/ControlSystemsBase/test/test_implicit_diff.jl
@@ -212,3 +212,17 @@ J2 = fdgrad(difffun, pars)[:]
 # @show norm(J1-J2)
 @test J1 ≈ J2 rtol = 1e-5
 
+## Test differentiation of tf 2 ss conversion
+
+function difffun(pars)
+    P = tf(1, pars)
+    sum(abs2, step(P, 0:0.1:10, method=:tustin).y + 
+        impulse(P, 0:0.1:10, method=:tustin).y
+    )
+end
+
+pars = [1.0, 2, 1]
+
+g1 = ForwardDiff.gradient(difffun, pars)
+g2 = fdgrad(difffun, pars)
+@test g1 ≈ g2 rtol = 1e-5
\ No newline at end of file

From aba649d9e351e4bd8f2fd2d7b79873d3913ae565 Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Fri, 14 Feb 2025 10:05:46 +0100
Subject: [PATCH 4/4] adjust tests to account for changes to internals

---
 lib/ControlSystemsBase/test/test_complex.jl    | 2 +-
 lib/ControlSystemsBase/test/test_conversion.jl | 4 ++--
 lib/ControlSystemsBase/test/test_timeresp.jl   | 4 ++--
 3 files changed, 5 insertions(+), 5 deletions(-)

diff --git a/lib/ControlSystemsBase/test/test_complex.jl b/lib/ControlSystemsBase/test/test_complex.jl
index 605256e17..0a44493ac 100644
--- a/lib/ControlSystemsBase/test/test_complex.jl
+++ b/lib/ControlSystemsBase/test/test_complex.jl
@@ -37,6 +37,6 @@ s = tf("s");
 @test tzeros(ss([-1.0-im 1-im; 2 0], [2; 0], [-1+1im -0.5-1.25im], 1)) ≈ [-1-2im, 2-im]
 
 @test tzeros(ss((s-2.0-1.5im)^3/(s+1+im)/(s+2)^3)) ≈ fill(2.0 + 1.5im, 3) rtol=1e-4
-@test tzeros(ss((s-2.0-1.5im)*(s-3.0)/(s+1+im)/(s+2)^2)) ≈ [3.0, 2.0 + 1.5im] rtol=1e-14
+@test tzeros(ss((s-2.0-1.5im)*(s-3.0)/(s+1+im)/(s+2)^2)) ≈ [2.0 + 1.5im, 3.0] rtol=1e-14
 
 end
diff --git a/lib/ControlSystemsBase/test/test_conversion.jl b/lib/ControlSystemsBase/test/test_conversion.jl
index 36576e46a..344d4ec43 100644
--- a/lib/ControlSystemsBase/test/test_conversion.jl
+++ b/lib/ControlSystemsBase/test/test_conversion.jl
@@ -8,8 +8,8 @@ A = randn(ComplexF64,3,3)
 
 G = tf(1.0,[1,1])
 H = zpk([0.0], [1.0], 1.0)
-@inferred ControlSystemsBase.siso_tf_to_ss(Float64, G.matrix[1,1])
-@inferred ControlSystemsBase.siso_tf_to_ss(Float64, H.matrix[1,1])
+# @inferred ControlSystemsBase.siso_tf_to_ss(Float64, G.matrix[1,1])
+# @inferred ControlSystemsBase.siso_tf_to_ss(Float64, H.matrix[1,1])
 
 # Easy second order system
 sys1 = ss([-1 0;1 1],[1;0],[1 1],0)
diff --git a/lib/ControlSystemsBase/test/test_timeresp.jl b/lib/ControlSystemsBase/test/test_timeresp.jl
index e078d4080..fdac5ea69 100644
--- a/lib/ControlSystemsBase/test/test_timeresp.jl
+++ b/lib/ControlSystemsBase/test/test_timeresp.jl
@@ -147,7 +147,7 @@ xreal = zeros(3, length(t0), 2)
 y, t, x = step(systf, t0, method=:zoh) # Method is :zoh so discretization is applied
 yreal[1,:,1] = 1 .- exp.(-t)
 yreal[2,:,2] = -1 .+ exp.(-t) + 2*exp.(-t).*t
-@test y ≈ yreal atol=1e-14
+@test y ≈ yreal atol=1e-13
 #Step ss
 y, t, x = step(sysss, t, method=:zoh)
 @test y ≈ yreal atol=1e-13
@@ -160,7 +160,7 @@ xreal[3,:,2] = exp.(-t).*(-t .- 1) .+ 1
 y, t, x = impulse(systf, t, method=:zoh)
 yreal[1,:,1] = exp.(-t)
 yreal[2,:,2] = exp.(-t).*(1 .- 2*t)
-@test y ≈ yreal atol=1e-14
+@test y ≈ yreal atol=1e-13
 #Impulse ss
 y, t, x = impulse(1.0sysss, t, method=:zoh)
 @test y ≈ yreal atol=1e-13