Replace exact and absolute comparisons in isunstable, baltrunc, delayd_ss and thiran

lib/ControlSystemsBase/src/delay_systems.jl

@@ -42,11 +42,12 @@ end
 Discrete-time statespace realization of a delay ``τ`` sampled with period ``T_s``,
 i.e. of ``z^{-N}`` where ``N = τ / T_s.``
 
-``τ`` must be a multiple of ``T_s``. See [`thiran`](@ref) for approximate discretization of fractional delays.
+``τ`` must be a multiple of ``T_s`` to within a relative tolerance ``\\sqrt{ϵ}``. See [`thiran`](@ref) for approximate discretization of fractional delays.
 """
 function delayd_ss(τ::Number, Ts::Number)
-    n = round(Int, τ / Ts)
-    if !(τ - n*Ts ≈ 0)
+    D = τ / Ts
+    n = round(Int, D)
+    if !isapprox(D, n; rtol = sqrt(eps(float(typeof(D)))))
         error("The delay τ must be a multiple of the sample time Ts, use the function `thiran` to approximately discretize fractional delays.")
     end
     ss(diagm(1 => ones(n-1)), [zeros(n-1,1); 1], [1 zeros(1,n-1)], 0, Ts)
@@ -182,12 +183,15 @@ Discretize a potentially fractional delay ``τ`` as a Thiran all-pass filter wit
 
 The Thiran all-pass filter gives an a maximally flat group delay.
 
-If ``τ`` is an integer multiple of ``Ts``, the Thiran all-pass filter reduces to ``z^{-τ/Ts}``.
+If ``τ`` is an integer multiple of ``Ts`` to within a relative tolerance ``\\sqrt{ϵ}``, the Thiran all-pass filter reduces to ``z^{-N}`` with ``N = \\operatorname{round}(τ/Ts)``.
 
 Ref: T. I. Laakso, V. Valimaki, M. Karjalainen and U. K. Laine, "Splitting the unit delay [FIR/all pass filters design]," in IEEE Signal Processing Magazine, vol. 13, no. 1, 1996.
 """
 function thiran(τ::Real, Ts)
     D = τ/Ts
+    # Rounding errors in τ/Ts must not increase the order of the filter
+    Dr = round(D)
+    isapprox(D, Dr; rtol = sqrt(eps(float(typeof(D))))) && (D = Dr)
     N = ceil(Int, D)
     a = ones(N+1)
     for k = 1:N

lib/ControlSystemsBase/src/matrix_comps.jl

@@ -734,9 +734,9 @@ end
 
 
 """
-    sysr, G, T = baltrunc(sys::StateSpace; atol = √ϵ, rtol=1e-3, n = nothing, residual = false)
+    sysr, G, T = baltrunc(sys::StateSpace; atol = 0, rtol=1e-3, n = nothing, residual = false)
 
-Reduces the state dimension by calculating a balanced realization of the system sys, such that the observability and reachability gramians of the balanced system are equal and diagonal `diagm(G)`, and truncating it to order `n`. If `n` is not provided, it's chosen such that all states corresponding to singular values less than `atol` and less that `rtol σmax` are removed.
+Reduces the state dimension by calculating a balanced realization of the system sys, such that the observability and reachability gramians of the balanced system are equal and diagonal `diagm(G)`, and truncating it to order `n`. If `n` is not provided, it is chosen such that all state variables corresponding to Hankel singular values less than `atol` or less than `rtol σmax` are removed. The Hankel singular values scale with the gain of the system, the default absolute tolerance `atol` is therefore zero and the truncation is determined by the relative tolerance `rtol`.
 
 `T` is the projection matrix between the old state `x` and the newstate `z` such that `z = Tx`. `T` will in general be a non-square matrix.
 
@@ -755,11 +755,11 @@ For more advanced model reduction, see [RobustAndOptimalControl.jl - Model Reduc
 # Extended help
 $(_scaling_notice)
 """
-function baltrunc(sys::ST; atol = sqrt(eps(numeric_type(sys))), rtol = 1e-3, n = nothing, residual=false) where ST <: AbstractStateSpace
+function baltrunc(sys::ST; atol = 0, rtol = 1e-3, n = nothing, residual=false) where ST <: AbstractStateSpace
     sysbal, S, T = balreal(sys)
     if n === nothing
         S = S[S .>= atol]
-        S = S[S .>= S[1]*rtol]
+        isempty(S) || (S = S[S .>= S[1]*rtol])
         n = length(S)
     else
         n > sys.nx && error("n too large. A state dimension of n = $n was requested, but the original system has a $(sys.nx)-dimensional state.")

lib/ControlSystemsBase/src/types/Lti.jl

@@ -119,20 +119,30 @@ isstable(sys::LTISystem{<:Discrete}) = all(abs.(poles(sys)) .< 1)
     isunstable(sys)
 
 Returns `true` if `sys` is exponentially unstable, else returns `false`.
-Marginally stable systems (systems with a simple poles on the imaginary axis) are considered stable by this function, see [`isstable`](@ref) for a function that returns true only for exponentially stable systems.
+Marginally stable systems (systems with simple poles on the imaginary axis, or on the unit circle in discrete time) are considered stable by this function, see [`isstable`](@ref) for a function that returns true only for exponentially stable systems.
 """
 function isunstable(sys::LTISystem)
     inte, p, z, tolp, tolz = integrator_excess_with_tol(sys)
     if inte > 1
         return true
     end
-    # Go through all poles on the imaginary axis and check if they are duplicated
-    for pi in p
-        abs(real(pi)) > sqrt(sqrt(eps(abs(pi)))) && continue # The pole is too far away to be on the imaginary axis
-        # check if there are multiple poles with this imaginary part
-        count_eigval_multiplicity(p, complex(0.0, imag(pi)))[1] > 1 && return true
+    if iscontinuous(sys)
+        # Go through all poles on the imaginary axis and check if they are duplicated
+        for pi in p
+            abs(real(pi)) > sqrt(sqrt(eps(abs(pi)))) && continue # The pole is too far away to be on the imaginary axis
+            # check if there are multiple poles with this imaginary part
+            count_eigval_multiplicity(p, complex(0.0, imag(pi)))[1] > 1 && return true
+        end
+        return any(real.(p) .> tolp)
+    else
+        # Go through all poles on the unit circle and check if they are duplicated
+        for pi in p
+            abs(abs(pi) - 1) > sqrt(sqrt(eps(abs(pi)))) && continue # The pole is too far away to be on the unit circle
+            # check if there are multiple poles with this argument
+            count_eigval_multiplicity(p, pi/abs(pi))[1] > 1 && return true
+        end
+        return any(abs.(p) .> 1 + tolp)
     end
-    return iscontinuous(sys) ? any(real.(p) .> tolp) : any(abs.(p) .> tolp)
 end
 
 # Fallback since LTISystem not AbstractArray

lib/ControlSystemsBase/test/test_analysis.jl

@@ -505,4 +505,19 @@ using ControlSystemsBase: isunstable
 @test isunstable(tf(1, [1,0,0])) # Double integrator
 @test isunstable(zpk([1], [im, im, -im, -im], 1)) # Repeated pole on imaginary axis not in origin
 @test !isunstable(zpk([1], [im+1e-8im, im, -im-1e-8im, -im], 1)) # Almost repeated pole on imaginary axis not in origin
-@test isunstable(zpk(Float64[], [0.01; -collect(1.0:9)], 1.0)) # Slightly unstable pole in a system with a large state dimension
+@test isunstable(zpk(Float64[], [0.01; -collect(1.0:9)], 1.0)) # Slightly unstable pole in a system with a large state dimension
+
+# Discrete time
+@test !isunstable(c2d(tf(1, [1,1]), 0.1))
+@test !isunstable(c2d(ss(tf(1, [1,1])), 0.1))
+@test isunstable(c2d(tf(1, [1,-1]), 0.1))
+@test isunstable(zpk(Float64[], [1.001, 0.5], 1.0, 1))
+@test !isunstable(tf(1, [1,0,0], 1)) # Double pole in the origin z = 0
+@test !isunstable(tf(1, [1,-1], 1)) # Simple pole at z = 1
+@test !isunstable(tf(1, [1,1], 1)) # Simple pole at z = -1
+@test isunstable(tf(1, [1,-2,1], 1)) # Double pole at z = 1
+@test isunstable(c2d(ss(tf(1, [1,0,0])), 0.1)) # Discretized double integrator
+@test isunstable(tf(1, [1,2,1], 1)) # Double pole at z = -1
+@test !isunstable(ss(c2d(tf(1, [1,0,1]), 0.1))) # Simple pole pair on the unit circle
+@test isunstable(ss(c2d(tf(1, [1,0,1])^2, 0.1))) # Repeated pole pair on the unit circle
+@test isunstable(zpk(Float64[], [cis(0.5), cis(0.5), cis(-0.5), cis(-0.5)], 1.0, 1))

lib/ControlSystemsBase/test/test_delayed_systems.jl

@@ -138,6 +138,10 @@ d = exp(-2*s)
 @test freqresp(ControlSystemsBase.delayd_ss(1.0, 0.2), Ω)[:] ≈ exp.(-im*Ω) atol=1e-14
 @test freqresp(ControlSystemsBase.delayd_ss(3.2, 0.4), Ω)[:] ≈ exp.(-3.2*im*Ω) atol=1e-14
 @test_throws ErrorException ControlSystemsBase.delayd_ss(3.2, 0.5)
+@test ControlSystemsBase.delayd_ss(0.3, 0.1).nx == 3 # 0.3/0.1 is not exactly 3 in floating-point arithmetic
+@test ControlSystemsBase.delayd_ss(0.1+0.2, 0.1).nx == 3
+@test_throws ErrorException ControlSystemsBase.delayd_ss(0.31, 0.1)
+@test c2d(tf(1, [1,1])*delay(0.3), 0.1).nx == 4
 
 # Simple tests for c2d of DelayLtiSystems
 @test freqrespv(c2d(feedback(ss(0,1,1,0), delay(1.5)), 0.5), Ω) ≈ [0.5/((z - 1) + 0.5*z^-3) for z in exp.(im*Ω*0.5)]
@@ -248,6 +252,10 @@ for Ts = [1, 1.1]
     z = tf('z', Ts)
     @test thiran(2Ts, Ts) == 1/z^2
 end
+let z = tf('z', 0.1) # τ/Ts is not exactly an integer in floating-point arithmetic
+    @test thiran(3*0.1, 0.1) == 1/z^3
+    @test thiran(6*0.1, 0.1) == 1/z^6
+end
 
 @test thiran(pi, 1) ≈ tf([-0.00031815668236122736, 0.0042438423976339556, -0.03424682772398137, 0.8290601401044773, 1.0], [1.0, 0.8290601401044773, -0.03424682772398137, 0.0042438423976339556, -0.00031815668236122736], 1)
 

lib/ControlSystemsBase/test/test_simplification.jl

@@ -63,6 +63,19 @@ y2,x2 = step(sysmin,t, method=:zoh)[[1,3]]
     sysr,_ = baltrunc(sys)
     @test norm(sys-sysr) < 0.02
 
+    # The truncation is invariant to the gain of the system
+    for sys0 in (sys, ss(tf(1, [1,2,3,1])))
+        sysr0, Σ0 = baltrunc(sys0)
+        for k in (1e-12, 1e-9, 1e-6, 1e6)
+            sysrk, Σk = baltrunc(k*sys0)
+            @test sysrk.nx == sysr0.nx
+            @test Σk ≈ k*Σ0 rtol=1e-8
+        end
+    end
+    sysr3, Σ3 = baltrunc(ss(-1.0, 0, 1, 0)) # Zero gain
+    @test sysr3.nx == 0
+    @test isempty(Σ3)
+
 
     sysr,Σ = baltrunc(sys, n=3, residual=true)
     @test sysr.nx == 3