Make the tolerance of the integrator count independent of the time unit

[baggepinnen]

count_eigval_multiplicity determines how many computed eigenvalues are located at a given point. It is used by count_integrators, integrator_excess, isunstable, and the phase adjustment (adjust_phase_start) in margin, marginplot and bodeplot. The default radius for multiplicity i was eps(ρ)^(1/i), where ρ = maximum(abs, p). For i ≥ 2, this radius is proportional to ρ^(1/i) rather than to ρ, so the result depends on the time unit of the model. For i = 1, the radius ϵρ leaves no margin for the backward error of the eigenvalue computation.

This PR changes the default radius to (100ϵ)^(1/i) * ρ.

Motivation of the tolerance

• A backward-stable eigenvalue algorithm returns the exact eigenvalues of A + E with ‖E‖ ≲ p(n)ϵ‖A‖. An eigenvalue that belongs to a Jordan block of size i is then perturbed by an amount of order ϵ^(1/i)‖A‖, i.e., proportionally to the scale of A. The radius must therefore scale linearly with the magnitude of the eigenvalues. ρ is used in place of ‖A‖ since only the eigenvalues are available (ρ ≤ ‖A‖).
• The constant: 100ϵ and 10√ϵ = (100ϵ)^(1/2) are the tolerances for simple and double imaginary-axis eigenvalues in the Bruinsma–Steinbuch algorithm for the H∞ norm, as implemented in DescriptorSystems.jl (norminfc/norminfd) and in SLICOT AB13DD. DescriptorSystems.jl and MatrixPencils.jl make other structural decisions, such as multiplicities of infinite eigenvalues, by rank decisions with the relative tolerance nϵ; that constant performed worse in the comparison below. MatrixEquations.jl contains no tolerance of this kind.

Further changes

• ϵ is the machine epsilon of the element type of the eigenvalues, so that models with Float32 coefficients are handled.
• Eigenvalues are counted if their distance is <= the radius, so that eigenvalues located exactly at the point are counted also when the scale is zero.
• When no multiplicity is found, the returned tolerance is the radius for multiplicity one instead of the radius for multiplicity n. isunstable and sisomargin use the returned tolerance as the threshold on the real part of unstable poles. With the radius for multiplicity n (0.033 for a system with state dimension 10), isunstable classified a system with a pole at +0.01 as not unstable.
• In integrator_excess, the zeros are counted with the scale max(maximum(abs, p), maximum(abs, z)). The zeros are computed from the system pencil, whose scale is at least that of the poles, and a single zero in the origin does not provide a scale of its own.
• An explicitly provided third argument e retains its meaning (radius e^(1/i)).

Examples

System	master	this PR
1/(s²(s+1)(s+2)) after an orthogonal similarity transformation, time unit scaled by 1e4	0 integrators	2
Poles -1e-3·(1:10), no integrator	10 integrators, margin reports a phase margin of -703.6° (16.4° when the time unit is scaled by 1000)	0, 16.4°
DemoSystems.double_mass_model() (the pole in the origin is computed as -1.1e-14)	0 integrators	1
Discrete-time integrator after a similarity transformation	0 integrators	1
s/(s+1)² after a rotation (the zero is computed as 2.3e-16)	integrator excess 0	-1
isunstable with a pole at +0.01, state dimension 10	false	true

A comparison on 2000 random realizations per case (modal form with one Jordan block of size m at the origin, 2–8 additional stable poles spread over three decades, an orthogonal similarity transformation, and a random time scale in [1e-4, 1e4]) gave the following numbers of wrong counts:

m	master	(ϵ)^(1/i)ρ	(nϵ)^(1/i)ρ	(100ϵ)^(1/i)ρ (this PR)
0	72	0	0	1
1	19	4	0	0
2	763	110	42	5
3	942	151	98	45

Two simple poles closer to the origin than approximately 1e-8 ρ cannot be distinguished from a perturbed double integrator by any tolerance of this form; this ambiguity is inherent. The existing test with pade(OL, k), which compares the count with a BigFloat computation, passes for all k.

The problem was found while reviewing #1076, where a pole of DemoSystems.double_mass_model() that is zero up to rounding error extended the frequency range of marginplot.

Tests

All test suites were run with Julia 1.13.0 on Linux, with this branch of ControlSystemsBase developed into each environment. The downstream packages were tested at their current default branches.

Package	Result with this PR	Result with master (ControlSystemsBase 1.23.0)
ControlSystemsBase	21336 passed, 6 broken	not run
ControlSystems	1151 passed	not run
ControlSystemsMTK	39 passed, 1 broken	not run
ModelPredictiveControl	1606 passed	not run
RobustAndOptimalControl	1497 passed, 1 failed, 6 broken	identical
DyadControlSystems	3483 passed, 2 failed, 1 errored, 20 broken	identical

The failures in RobustAndOptimalControl and DyadControlSystems occur identically with master and do not involve the functions changed here:

• RobustAndOptimalControl test_uncertainty.jl:42: rand(Δ(2, δc), 100) returns a Diagonal rather than a Matrix on Julia 1.13.
• DyadControlSystems test/MPC/test_collocation.jl:140: an allocation test (160 bytes).
• DyadControlSystems test/MPC/test_quadtank.jl:209: a convergence tolerance of the linear MPC; the final state values are bit-identical with master.
• DyadControlSystems, Dyad frequency-response analysis: bodeplot of a ControlSystemIdentification.FRD throws a MethodError, since bodeplot passes the keyword balance (added in Balance the realization in the frequency-response and norm functions #1071) to bode(::FRD, w; unwrap), which does not accept it.

---

Host: demeter2, session 2071354d-b6c5-40db-a9b8-ccbe8a7d8c38, working directory ~/.julia/dev/ControlSystems/.claude/worktrees/joyful-gathering-pumpkin

🤖 Generated with Claude Code

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 91.69%. Comparing base (22ff8c1) to head (3b92a27).

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #1077      +/-   ##
==========================================
+ Coverage   91.64%   91.69%   +0.04%     
==========================================
  Files          42       42              
  Lines        5795     5789       -6     
==========================================
- Hits         5311     5308       -3     
+ Misses        484      481       -3     

☔ 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.