From c2da3540439ac53be47be983ff569f5b12ba9651 Mon Sep 17 00:00:00 2001 From: Fredrik Bagge Carlson Date: Fri, 9 Oct 2026 04:47:22 +0000 Subject: [PATCH] Replace exact and absolute comparisons in isunstable, baltrunc, delayd_ss and thiran MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - isunstable: in discrete time, compare the pole magnitudes with 1 instead of with the multiplicity tolerance, and check repeated poles on the unit circle instead of on the imaginary axis. - baltrunc: default atol = 0, since the Hankel singular values scale with the gain of the system. Handle the case where all singular values are removed. - delayd_ss: test that τ/Ts is an integer to within a relative tolerance instead of requiring τ - n*Ts == 0. - thiran: round τ/Ts to the nearest integer when it is within a relative tolerance, so that rounding errors do not increase the filter order. Part of #1081 Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01HigNnun4aJKMoWmY7wyULn --- lib/ControlSystemsBase/src/delay_systems.jl | 12 ++++++---- lib/ControlSystemsBase/src/matrix_comps.jl | 8 +++---- lib/ControlSystemsBase/src/types/Lti.jl | 24 +++++++++++++------ lib/ControlSystemsBase/test/test_analysis.jl | 17 ++++++++++++- .../test/test_delayed_systems.jl | 8 +++++++ .../test/test_simplification.jl | 13 ++++++++++ 6 files changed, 66 insertions(+), 16 deletions(-) diff --git a/lib/ControlSystemsBase/src/delay_systems.jl b/lib/ControlSystemsBase/src/delay_systems.jl index a50721234..14257cb48 100644 --- a/lib/ControlSystemsBase/src/delay_systems.jl +++ b/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 diff --git a/lib/ControlSystemsBase/src/matrix_comps.jl b/lib/ControlSystemsBase/src/matrix_comps.jl index 03514ee65..6107a4342 100644 --- a/lib/ControlSystemsBase/src/matrix_comps.jl +++ b/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.") diff --git a/lib/ControlSystemsBase/src/types/Lti.jl b/lib/ControlSystemsBase/src/types/Lti.jl index 719cddb74..6e85674b9 100644 --- a/lib/ControlSystemsBase/src/types/Lti.jl +++ b/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 diff --git a/lib/ControlSystemsBase/test/test_analysis.jl b/lib/ControlSystemsBase/test/test_analysis.jl index f05a4d674..f2bd7d195 100644 --- a/lib/ControlSystemsBase/test/test_analysis.jl +++ b/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 \ No newline at end of file +@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)) \ No newline at end of file diff --git a/lib/ControlSystemsBase/test/test_delayed_systems.jl b/lib/ControlSystemsBase/test/test_delayed_systems.jl index 47d836841..b0c32eba2 100644 --- a/lib/ControlSystemsBase/test/test_delayed_systems.jl +++ b/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) diff --git a/lib/ControlSystemsBase/test/test_simplification.jl b/lib/ControlSystemsBase/test/test_simplification.jl index 33a108d89..40c999cf3 100644 --- a/lib/ControlSystemsBase/test/test_simplification.jl +++ b/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