Merge pull request #906 from JuliaControl/ss_inv

handle proper quotient of proper statespace systems

lib/ControlSystemsBase/src/discrete.jl

@@ -21,7 +21,7 @@ ZoH sampling is exact for linear systems with piece-wise constant inputs (step i
 
 FoH sampling is exact for linear systems with piece-wise linear inputs (ramp invariant), this is a good choice for simulation of systems with smooth continuous inputs.
 
-To approximate the behavior of a continuous-time system well in the frequency domain, the `:tustin` (trapezoidal / bilinear) method may be most appropriate. In this case, the pre-warping argument can be used to ensure that the frequency response of the discrete-time system matches the continuous-time system at a given frequency. The tustin transformation alters the meaning of the state components, while ZoH and FoH preserve the meaning of the state components. The Tustin method is commonly used to discretize a continuous-tiem controller.
+To approximate the behavior of a continuous-time system well in the frequency domain, the `:tustin` (trapezoidal / bilinear) method may be most appropriate. In this case, the pre-warping argument can be used to ensure that the frequency response of the discrete-time system matches the continuous-time system at a given frequency. The tustin transformation alters the meaning of the state components, while ZoH and FoH preserve the meaning of the state components. The Tustin method is commonly used to discretize a continuous-time controller.
 
 The forward-Euler method generally requires the sample time to be very small
 relative to the time constants of the system, and its use is generally discouraged.
@@ -99,8 +99,8 @@ function d2c(sys::AbstractStateSpace{<:Discrete}, method::Symbol=:zoh; w_prewarp
     ny, nu = size(sys)
     nx = nstates(sys)
     if method === :zoh
-        M = log([A  B;
-            zeros(nu, nx) I])./sys.Ts
+        M = log(Matrix([A  B;
+            zeros(nu, nx) I]))./sys.Ts
         Ac = M[1:nx, 1:nx]
         Bc = M[1:nx, nx+1:nx+nu]
         if eltype(A) <: Real

lib/ControlSystemsBase/src/types/StateSpace.jl

@@ -378,15 +378,64 @@ end
 
 
 ## DIVISION ##
-/(sys1::AbstractStateSpace, sys2::AbstractStateSpace) = sys1*inv(sys2)
+
+
+"""
+    /(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)))
+
+Compute `sys1 / sys2 = sys1 * inv(sys2)` in a way that tries to handle situations in which the inverse `sys2` is non-proper, but the resulting system `sys1 / sys2` is proper.
+
+See `ControlSystemsBase.MatrixPencils.isregular` for keyword arguments `atol`, `atol1`, `atol2`, and `rtol`.
+"""
+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([sys2; sys1])
+        Ai = [A B; C[1:ny2,:] D[1:ny2,:]]
+        Ei = [I zeros(T,nx1+nx2,ny2); zeros(T,ny2,nx1+nx2+ny2)] |> 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 sys2 is not invertible")
+        Ci = [C[ny2+1:ny1+ny2,:] D[ny2+1:ny1+ny2,:]]
+        Bi = [zeros(T,nx1+nx2,nu1); -I] |> Matrix # TODO: rm call to Matrix when type piracy in https://github.com/JuliaLinearAlgebra/LinearMaps.jl/issues/219 is fixed
+        Di = zeros(T,ny1,nu1)
+        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 sys2 is not invertible"))
+            Ai = luE\Ai
+            Bi = luE\Bi
+        end
+    else
+        D2 = T.(sys2.D)
+        LUD = lu(D2)
+        (norm(D2,Inf) <= atol1 || rcond(LUD.U) <= 10*nu1*eps(real(float(one(T))))) && 
+                error("The system sys2 is not invertible")
+        Ai, Bi, Ci, Di = ssdata(sys1)
+        rdiv!(Bi,LUD); rdiv!(Di,LUD)
+    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)
+    size(D, 1) == size(D, 2) || error("The inverted system must have the same number of inputs and outputs")
     Dinv = try
         inv(D)
     catch
-        error("D isn't invertible")
+        error("D isn't invertible. If you are trying to form a quotient between two systems `N(s) / D(s)` where the quotient is proper but the inverse of `D(s)` isn't, consider calling `N / D` instead of `N * inv(D)")
     end
     return basetype(ST)(A - B*Dinv*C, B*Dinv, -n*Dinv*C, n*Dinv, sys.timeevol)
 end

lib/ControlSystemsBase/test/test_statespace.jl

@@ -159,11 +159,18 @@
 
         # Division
         @test 1/C_222_d == SS([-6 -3; 2 -11],[1 0; 0 2],[-1 0; -0 -1],[1 -0; 0 1])
-        @test C_221/C_222_d == SS([-5 -3 -1 0; 2 -9 -0 -2; 0 0 -6 -3;
-        0 0 2 -11],[1 0; 0 2; 1 0; 0 2],[1 0 0 0],[0 0])
+        @test hinfnorm((C_221/C_222_d) - SS([-5 -3 -1 0; 2 -9 -0 -2; 0 0 -6 -3;
+        0 0 2 -11],[1 0; 0 2; 1 0; 0 2],[1 0 0 0],[0 0]))[1] < 1e-10
         @test 1/D_222_d == SS([-0.8 -0.8; -0.8 -1.93],[1 0; 0 2],[-1 0; -0 -1],
         [1 -0; 0 1],0.005)
 
+        # Division when denominator inverse is non-proper but quotient is proper
+        G1 = tf([1], [1, 1])
+        G2 = tf([1], [1, 1, 1])
+        G1s = ss(G1)
+        G2s = ss(G2)
+        @test tf(G2s / G1s) ≈ G2 / G1
+
         fsys = ss(1,1,1,0)/3 # Int becomes FLoat after division
         @test fsys.B[]*fsys.C[] == 1/3