add \ for ss and handle some cases of multiplication with improper tf (#930)

* add \ for ss and handle some cases of multiplication with improper tf

* bump version

* fix stray end

* ch error type

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.10.2"
+version = "1.10.3"
 
 [deps]
 DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2"

lib/ControlSystemsBase/src/types/Lti.jl

@@ -1,7 +1,30 @@
 abstract type LTISystem{TE<:TimeEvolution} <: AbstractSystem end
 +(sys1::LTISystem, sys2::LTISystem) = +(promote(sys1, sys2)...)
 -(sys1::LTISystem, sys2::LTISystem) = -(promote(sys1, sys2)...)
-*(sys1::LTISystem, sys2::LTISystem) = *(promote(sys1, sys2)...)
+function *(sys1::LTISystem, sys2::LTISystem)
+    zeroexcess(sys) = length(numvec(sys)[]) - length(denvec(sys)[])
+    try
+        *(promote(sys1, sys2)...)
+    catch e
+        e isa ImproperException || rethrow()
+        if sys1 isa AbstractStateSpace && sys2 isa TransferFunction
+            issiso(sys2) || rethrow() # Can't invert MIMO tf
+            zeroexcess(sys2) <= sys1.nx || rethrow() # Quotient can't be proper
+            Q = sys1 / ss(inv(sys2))
+        elseif sys1 isa TransferFunction && sys2 isa AbstractStateSpace
+            issiso(sys1) || rethrow() # Can't invert MIMO tf
+            zeroexcess(sys1) <= sys2.nx || rethrow() # Quotient can't be proper
+            Q = ss(inv(sys1)) \ sys2
+        else
+            rethrow()
+        end
+        Q = balance_statespace(Q)[1]
+        if max(maximum(abs, Q.A), maximum(abs, Q.B), maximum(abs, Q.C), maximum(abs, Q.D)) > 1e7
+            @warn "Possible numerical instability detected: Multiplication of a statespace system and a non-proper transfer function may result in numerical inaccuracy. Verify result carefully, and consider making use of DescriptorSystems.jl to represent this product as a DescriptorSystem with non-unit descriptor matrix if result is inaccurate."
+        end
+        return Q
+    end
+end
 /(sys1::LTISystem, sys2::LTISystem) = /(promote(sys1, sys2)...)
 
 # Fallback number

lib/ControlSystemsBase/src/types/StateSpace.jl

@@ -428,6 +428,47 @@ function Base.:(/)(sys1::AbstractStateSpace{TE}, sys2::AbstractStateSpace{TE};
     return StateSpace{TE, T}(Ai, Bi, Ci, Di, timeevol) 
 end
 
+function Base.:(\)(sys1::AbstractStateSpace{TE}, sys2::AbstractStateSpace{TE}; 
+    atol::Real = 0, atol1::Real = atol, atol2::Real = atol, 
+    rtol::Real = max(size(sys1.A,1),size(sys2.A,1))*eps(real(float(one(numeric_type(sys1)))))*iszero(min(atol1,atol2))) where {TE<:ControlSystemsBase.TimeEvolution}
+    T1 = float(numeric_type(sys1))
+    T2 = float(numeric_type(sys2))
+    T = promote_type(T1,T2)
+    timeevol = common_timeevol(sys1, sys2)
+    ny2, nu2 = sys2.ny, sys2.nu
+    nu2 == ny2  || error("The system sys2 must be square")
+    ny1, nu1 = sys1.ny, sys1.nu
+    nu1 == nu2  || error("The systems sys1 and sys2 must have the same number of inputs")
+    nx1 = sys1.nx
+    nx2 = sys2.nx
+    if nx2 > 0
+        A, B, C, D = ssdata([sys1 sys2])
+        Ai = [A B[:,1:nu1]; C D[:,1:nu1]]
+        Ei = [I zeros(T,nx1+nx2,nu1); zeros(T,nu1,nx1+nx2+nu1)] |> Matrix # TODO: rm call to Matrix when type piracy in https://github.com/JuliaLinearAlgebra/LinearMaps.jl/issues/219 is fixed
+        MatrixPencils.isregular(Ai, Ei; atol1, atol2, rtol) || 
+            error("The system sys1 is not invertible")
+        Bi = [B[:,nu1+1:nu1+nu2]; D[:,nu1+1:nu1+nu2]]
+        Ci = [zeros(T,ny1,nx1+nx2) -I] 
+        Di = zeros(T,ny1,nu2)
+        Ai, Ei, Bi, Ci, Di = MatrixPencils.lsminreal(Ai, Ei, Bi, Ci, Di; fast = true, atol1 = 0, atol2, rtol, contr = true, obs = true, noseig = true)
+        if Ei != I
+            luE = lu!(Ei, check=false)
+            issuccess(luE) || throw(ArgumentError("The system sys1 is not invertible"))
+            Ai = luE\Ai
+            Bi = luE\Bi
+        end
+    else
+        D1 = T.(sys1.D)
+        LUD = lu(D1)
+        (norm(D1,Inf) <= atol1 || rcond(LUD.U) <= 10*nu1*eps(real(float(one(T))))) && 
+                error("The system sys1 is not invertible")
+        Ai, Bi, Ci, Di = ssdata(sys2)
+        ldiv!(LUD, Ci); ldiv!(LUD, Di)
+    end
+
+    return StateSpace{TE, T}(Ai, Bi, Ci, Di, timeevol) 
+end
+
 function /(n::Number, sys::ST) where ST <: AbstractStateSpace
     # Ensure s.D is invertible
     A, B, C, D = ssdata(sys)
@@ -445,6 +486,22 @@ Base.inv(sys::AbstractStateSpace) = 1/sys
 Base.:\(n::Number, sys::ST) where ST <: AbstractStateSpace = basetype(ST)(sys.A, sys.B, sys.C/n, sys.D/n, sys.timeevol)
 
 
+function Base.adjoint(sys::ST) where ST <: AbstractStateSpace{Continuous}
+    return basetype(ST)(-sys.A', -sys.C', sys.B', sys.D') 
+end
+
+function Base.adjoint(sys::ST) where ST <: AbstractStateSpace{<:Discrete}
+    nx, ny, nu = sys.nx, sys.ny, sys.nu
+    T = numeric_type(sys)
+    return basetype(ST)(
+        [sys.A' sys.C'; zeros(T,ny,nx+ny)], 
+        [zeros(T,nx,ny) ; -I], 
+        [sys.B' zeros(T,nu,ny)], copy(sys.D'),
+        sys.timeevol
+    ) 
+ end
+
+
 #####################################################################
 ##                       Indexing Functions                        ##
 #####################################################################

lib/ControlSystemsBase/src/types/conversion.jl

@@ -86,10 +86,14 @@ function Base.convert(::Type{StateSpace}, G::TransferFunction{TE,<:SisoTf{T0}};
     convert(StateSpace{TE,T}, G; kwargs...)
 end
 
+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)
-        error("System is improper, a state-space representation is impossible")
+        throw(ImproperException())
     end
 
     ny, nu = size(G)

lib/ControlSystemsBase/test/test_conversion.jl

@@ -177,7 +177,7 @@ Hd = zpk([], [1, 0.5], 1.0, h)
 
 
 # Error, not proper
-@test_throws ErrorException ss(tf([1,0],[1]))
+@test_throws ControlSystemsBase.ImproperException ss(tf([1,0],[1]))
 
 s1 = HeteroStateSpace(zeros(Int, 1,1),zeros(Int, 1,1),zeros(Int, 1,1),zeros(Int, 1,1))
 s2 = HeteroStateSpace(zeros(Float64, 1,1),zeros(Float64, 1,1),zeros(Float64, 1,1),zeros(Float64, 1,1))
@@ -201,4 +201,28 @@ sys.matrix[1,1] = 0
 syszpk = zpk(sys)
 @test syszpk isa TransferFunction{Continuous, ControlSystemsBase.SisoZpk{Float64, ComplexF64}}
 
+## Test multiplication of improper transfer function and statespace system
+P = DemoSystems.double_mass_model()
+C = tf('s')
+@test norm(P*C - C*P) < 1e-10
+mp = minreal(minreal(tf(P)*C) - P*C)
+@test hinfnorm(mp)[1] < 1e-10
+@inferred P*C
+
+@test_logs (:warn,"Possible numerical instability detected: Multiplication of a statespace system and a non-proper transfer function may result in numerical inaccuracy. Verify result carefully, and consider making use of DescriptorSystems.jl to represent this product as a DescriptorSystem with non-unit descriptor matrix if result is inaccurate.") P*tf('s')^3
+
+@test_throws ControlSystemsBase.ImproperException P*tf('s')^5
+# Discrete
+
+P = c2d(P, 0.01, :fwdeuler) # to get poles exactly at 1
+C = tf('z', 0.01)-1
+@test norm(minreal(P*C - C*P)) < 1e-10
+mp = minreal(minreal(tf(P)*C) - P*C)
+@test hinfnorm(mp)[1] < 1e-10
+@inferred P*C
+
+# @test_logs (:warn,"Possible numerical instability detected: Multiplication of a statespace system and a non-proper transfer function may result in numerical inaccuracy. Verify result carefully, and consider making use of DescriptorSystems.jl to represent this product as a DescriptorSystem with non-unit descriptor matrix if result is inaccurate.") P*C^3
+
+@test_throws ControlSystemsBase.ImproperException P*C^5
+
 end

lib/ControlSystemsBase/test/test_statespace.jl

@@ -175,6 +175,8 @@
         G2s = ss(G2)
         @test tf(G2s / G1s) ≈ G2 / G1
 
+        @test tf(G1s \ G2s) ≈ G2 / G1
+
         fsys = ss(1,1,1,0)/3 # Int becomes FLoat after division
         @test fsys.B[]*fsys.C[] == 1/3