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

[baggepinnen]

Part of #1081.

This PR replaces four exact or absolute comparisons in ControlSystemsBase that produce wrong results or errors for realistic systems.

1. isunstable in discrete time

The discrete-time branch compared the pole magnitudes with the multiplicity tolerance tolp (of the order of 100ϵ) instead of with 1, so every discrete-time system with a nonzero pole was reported as unstable. In addition, the loop that detects repeated poles was shared with the continuous-time branch and searched for repeated poles on the imaginary axis, so a double pole at z = 0 (e.g. z^-2) was reported as unstable, while a double pole at z = -1 was detected only through the first error.

using ControlSystemsBase: isunstable
isunstable(c2d(tf(1, [1,1]), 0.1)) # true, expected false
isunstable(tf(1, [1,0,0], 1))      # true, expected false (double pole at z = 0)
isunstable(tf(1, [1,-1], 1))       # true, expected false (simple pole at z = 1)

Change: the discrete-time branch searches for poles on the unit circle (abs(abs(p) - 1) below the same prefilter tolerance as the continuous-time branch) and counts the multiplicity with count_eigval_multiplicity around the projection p/abs(p) onto the unit circle. The final test is abs(p) > 1 + tolp. Simple poles on the unit circle (marginal stability) are classified as stable, repeated poles on the unit circle (at z = 1, z = -1 or a complex pair) as unstable, consistent with the docstring and with the continuous-time branch. The continuous-time branch is unchanged; its known weaknesses are listed in #1081.

2. baltrunc default atol

The default atol = sqrt(eps()) is an absolute threshold on the Hankel singular values, which scale with the gain of the system.

sys = ss(tf(1, [1,2,3,1]))
baltrunc(1e-9*sys) # BoundsError (all singular values removed, then S[1] is read)
baltrunc(ss(-1.0, 0, 1, 0)) # BoundsError (zero gain)

Change: the default is atol = 0. The truncation is then determined by the relative tolerance rtol and by the relative truncation to a minimal realization in balreal (sqrt(eps)*σmax). The relative truncation is skipped when no singular values remain, so a system of zero gain returns a system with a state dimension of zero. The docstring is updated. No other function in the package calls baltrunc.

3. delayd_ss

τ - n*Ts ≈ 0 is an exact test, since isapprox with zero and the default atol = 0 requires equality.

ControlSystemsBase.delayd_ss(0.3, 0.1)       # error: "The delay τ must be a multiple of the sample time Ts"
c2d(tf(1, [1,1])*delay(0.3), 0.1)           # same error

Change: isapprox(τ/Ts, n; rtol = sqrt(eps(T))) with n = round(Int, τ/Ts).

4. thiran

N = ceil(Int, τ/Ts) increases the order when τ/Ts exceeds an integer by a rounding error.

thiran(3*0.1, 0.1) # fourth-order filter with a pole of magnitude 0.9999999999999996 instead of z^-3
thiran(6*0.1, 0.1) # seventh-order filter instead of z^-6

Change: D = τ/Ts is replaced by round(D) when it is within the relative tolerance sqrt(eps(T)) of it, before N = ceil(Int, D). The docstring states the tolerance.

An alternative is N = round(Int, D), i.e. D ∈ [N-0.5, N+0.5), which Laakso et al. recommend for the Thiran filter. This would change the order of the filter for existing users whenever the fractional part of D is below 0.5, and is therefore not included in this PR.

Effect on existing users

• isunstable returns false for discrete-time systems with all poles inside the unit circle or with simple poles on it; previously it returned true for every discrete-time system with a nonzero pole.
• baltrunc with default arguments changes the result only if the largest Hankel singular value is smaller than sqrt(eps())/rtol (about 1.5e-5 with the default rtol), i.e. for systems of small gain, for which it previously removed state variables based on the absolute threshold or threw an error. Users who passed atol explicitly are not affected.
• delayd_ss/c2d of delay systems accept delays that are multiples of the sample time up to rounding errors instead of throwing.
• thiran returns z^-N for delays that are multiples of the sample time up to rounding errors, instead of a filter of order N+1.

Tests

Regression tests are added to test_analysis.jl (discrete-time isunstable: stable, unstable, double pole at z = 0, simple and double poles at z = 1 and z = -1, simple and repeated pole pairs on the unit circle, discretized double integrator), test_simplification.jl (baltrunc of k*sys for k ∈ {1e-12, 1e-9, 1e-6, 1e6} gives the same state dimension and Hankel singular values scaled by k; zero-gain system) and test_delayed_systems.jl (delayd_ss(0.3, 0.1), delayd_ss(0.1+0.2, 0.1), c2d of a delay system, rejection of delayd_ss(0.31, 0.1), thiran(3*0.1, 0.1) == 1/z^3, thiran(6*0.1, 0.1) == 1/z^6).

Run locally (Julia 1.12): test_analysis.jl, test_simplification.jl, test_matrix_comps.jl, test_delayed_systems.jl, test_discrete.jl: 513 passed, 0 failed, 0 errored.

---

Host: demeter6, Claude Code session session_01HigNnun4aJKMoWmY7wyULn, working directory ~/.julia/dev/ControlSystems/.claude/worktrees/agent-a4840bb7eb32e469c

🤖 Generated with Claude Code

https://claude.ai/code/session_01HigNnun4aJKMoWmY7wyULn

[JuliaControlBot]

This is an automated message.
Plots were compared to references. 4/11 images have changed, see differences below.
After pulling this PR, please update the reference images by creating a PR to ControlExamplePlots.jl here.

Difference	Reference Image	New Image
✔️ 0.0
✔️ 0.0
✔️ 0.0
✔️ 0.013

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 91.74%. Comparing base (59c065a) to head (c2da354).

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #1083      +/-   ##
==========================================
+ Coverage   91.73%   91.74%   +0.01%     
==========================================
  Files          42       42              
  Lines        5865     5873       +8     
==========================================
+ Hits         5380     5388       +8     
  Misses        485      485              

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.