Improve the default frequency vector: features below 1e-4 rad/s, discrete-time systems, lightly damped resonances and delays

[baggepinnen]

Part of #1081.

Problems

All reproducers were checked on top of #1076 (head 1ff5e00), which has since been merged; the branch is rebased on master after the merge of #1076 (f8c5db4).

1. Absolute cutoff at 1e-4 rad/s. Features below 1e-4 rad/s were ignored for bode, nyquist and sigma. For G = ss(tf(1,[1,1])), _default_freq_vector(time_scale(G, a), Val(:bode)) has the extrema [1, 100] for a = 1e-4, 1e-5, 1e-6, i.e., the grid lies on the high-frequency asymptote. For tf(1,[1,1])*tf(1e-5,[1,1e-5]) the grid starts at 0.01 rad/s, where |G| = 1e-3 although the DC gain is 1.
2. Discrete-time systems. log10(abs(z)) of the z-plane roots was used as a log-frequency, and the upper bound log10(π/Ts) was not compared with the lower bound. For c2d(tf(1,[1,1]), Ts) with Ts = 10, 100, 1000 the grid is unsorted and 99 % of its points lie above the Nyquist frequency. For a pole at 1e-3 rad/s sampled with Ts = 0.01 the grid starts at 0.01 rad/s. For tf(0.1,[1,-0.9],60) the grid is [0.01, 0.052] although the equivalent pole is located at 1.8e-3 rad/s.
3. Lightly damped resonances. The grid has a fixed logarithmic density and does not resolve resonances. Peak of tf(1,[1,2ζ,1]), true value vs. bode default grid: ζ = 1e-2: 50 vs. 23.5; ζ = 1e-3: 500 vs. 26.4; ζ = 1e-4: 5000 vs. 26.5. For margin this results in missed crossovers:
   • L = 0.02*tf(1,[1,2e-3,1])*tf(1,[0.1,1]): no gain crossover is found (pm = Inf) and gm = 1.529, whereas a dense grid gives pm = 168.6° at 0.99005 rad/s, pm = 0.059° at 1.00985 rad/s and gm = 1.0102.
   • The example of margin throws an error when the system has delay #941 in test_analysis.jl (C*G*H*delay(1) with a resonance at 15.97 rad/s): the default grid misses the phase crossover at 15.9626 rad/s, so the reported smallest gain margin is 28.6 instead of 0.1143 (dense grid).
4. Delays. _bounds_and_features(::DelayLtiSystem) used log10(τ) as a log-frequency and only considered τ > 1e4, which extended the grid upward to approximately τ rad/s. For τ = 1e-3 and 1e3 the delay did not affect the grid, for τ = 1e5 the grid extended to 1e6 rad/s.

Changes

• _remove_origin_roots (renamed from _margin_nonintegrators of Fix default margin frequency ranges and missing phase guides #1076, since it is now used by all plot types) removes only the poles and zeros located in the origin (at z = 1 in discrete time), using the classification of integrator_excess_with_tol via count_eigval_multiplicity. All other features are kept, so the bounds are equivariant under time_scale up to the rounding to whole decades. The 1e-4 cutoff is removed for all plot types.
• The radius of count_eigval_multiplicity is (100ϵ)^(1/n)*scale, which approaches scale for large n. For high-order systems with poles distributed over several decades, it therefore reports a large multiplicity at the origin: count_integrators of a sum of 30 lightly damped modes with ωn in [0.1, 100] returns 54, and integrator_excess of prod(tf(1,[1/10.0^k,1]) for k in -4:5) returns 8. With the classification of Fix default margin frequency ranges and missing phase guides #1076 alone, the grid of the first example started at 1 rad/s (Fix default margin frequency ranges and missing phase guides #1076 has the same behaviour for its margin grid). _remove_origin_roots therefore accepts a cluster only if the mean of its roots is located at the origin up to rounding error: the computed roots of an n-fold root are spread over a radius proportional to ϵ^(1/n), but their mean is perturbed by an amount proportional to ϵ only. For Jordan blocks of size 1 to 5 in 200 random coordinate transformations each, the mean test never rejected a cluster that count_eigval_multiplicity had found (largest ratio of the mean to the tolerance: 0.86 for n = 1, ≤ 2e-3 for n ≥ 2). This change modifies the margin grid of Fix default margin frequency ranges and missing phase guides #1076 only for such systems. The defect of count_eigval_multiplicity itself, which also affects integrator_excess and therefore the phase adjustment in sisomargin, is not changed here.
• Discrete-time roots are mapped to the s-plane by s = log(z)/Ts (a root at z = 0 has no finite equivalent and is excluded). Features above the Nyquist frequency are treated as features at the Nyquist frequency, the upper bound is π/Ts, and the lower bound is at least two decades below π/Ts. This replaces the w1 = min(w1, w2 - 1) rule that Fix default margin frequency ranges and missing phase guides #1076 introduced for discrete-time margin grids; the margin-specific extension of Fix default margin frequency ranges and missing phase guides #1076 below the lower bound is unchanged.
• _resonance_frequencies adds frequencies close to each lightly damped pole or zero pair with natural frequency ωn and damping ratio ζ: ωp(1 ± ζ sinh(u)) for u = 0, 1/8, 2/8, ..., where ωp = ωn sqrt(1 - 2ζ²) is the frequency of the extremum of the magnitude. Close to ωp the spacing is ζωp/8, far from it the frequencies are logarithmically spaced with ratio approximately exp(1/8). Frequencies are added only as long as this spacing is smaller than that of the logarithmic grid, so no frequencies are added for well-damped pairs (the grid of tf(1,[1,1,1]) is unchanged). For an undamped pair, ωp itself is omitted, so the response is not evaluated at a pole. For ζ = 1e-3 this adds approximately 100 frequencies per pair, for 30 such modes (60 pole and zero pairs) the grid has 6497 frequencies. The returned vector is sorted and unique but not logarithmically equidistant.
• Callers that assumed logarithmically uniform spacing: downsample approximated the data linearly in the index, which is equivalent to linear approximation in log(w) only for logarithmically equidistant grids. It now approximates linearly in log(t) when all t are positive and in the index otherwise (Nyquist curves, where the index is a common parameter of the real and imaginary parts and log(w) is equivalent for logarithmically spaced grids). The linear approximation is written as a convex combination, so that it equals the data at the segment end points without rounding error; otherwise the loop does not terminate when the rounding error exceeds the tolerance (this occurred in test_pid_design.jl). The resonance frequencies are added after the grid extension of Fix default margin frequency ranges and missing phase guides #1076, which computes a logarithmic step from the uniform Bode grid. _findCrossings/interpolate in sisomargin, unwrap and stabregionPID do not assume uniform spacing.
• _bounds_and_features(::DelayLtiSystem) adds a real feature at -1/τ for each delay τ > 0 and recomputes the bounds with the same rules as for poles and zeros.
• The docstrings of bode, nyquist and sigma describe the default frequency vector.

Results after the change

	before	after
time_scale(G, a), a = 1e-4 / 1e-5 / 1e-6	[1, 100]	[1e-6, 1e-2] / [1e-7, 1e-3] / [1e-8, 1e-4]
poles 1 and 1e-5: |G(w[1])|	1e-3	0.99995 (grid starts at 1e-7)
c2d(tf(1,[1,1]), Ts), Ts = 10 / 100 / 1000	unsorted, 99 % above π/Ts	sorted, [π/Ts/100, π/Ts]
pole 1e-3 rad/s, Ts = 0.01	[0.01, 314]	[1e-5, 314]
tf(0.1,[1,-0.9],60)	[0.01, 0.052]	[1e-4, 0.052]
resonance peak ζ = 1e-2 / 1e-3 / 1e-4	23.5 / 26.4 / 26.5	50 / 500 / 5000
margin of L above	pm = Inf, gm = 1.529	pm = [168.64, 0.049]° at [0.99004, 1.00987], gm = 1.0109
#941 example, smallest gm	28.6 at 20.47 rad/s	0.11436 at 15.9626 rad/s
delay τ = 1e-3 / 1e3 / 1e5	[0.01, 100] / [0.01, 100] / [0.01, 1e6]	[0.01, 1e5] / [1e-5, 100] / [1e-7, 100]

The phase margin 0.049° of L differs by 0.01° from the dense-grid value 0.059°, since the gain crossover is interpolated linearly; the phase margin is the difference of two terms of approximately 5.7°.

Effect on existing users

• Default grids of systems with features below 1e-4 rad/s, of discrete-time systems, of systems with lightly damped poles or zeros, and of delay systems change. The number of points increases by approximately 100 per lightly damped pair with ζ = 1e-3.
• bode(sys), nyquist(sys) and sigma(sys) return a frequency vector that is not logarithmically equidistant if the system has lightly damped poles or zeros.
• margin without a frequency vector finds crossovers close to resonances that were previously missed (see the margin throws an error when the system has delay #941 example, whose reported smallest gain margin changes from 28.6 to 0.114).
• downsample (internal, used by the plot recipes) approximates the data linearly in log(t) for positive t.

Tests

• New testset "Default frequency vector" in test_freqresp.jl: equivariance under time_scale (a = 1e-6, 1e-4, 1e3) for all plot types, slow poles, integrators, classification of poles distributed over several decades, discrete-time grids for Ts from 1e-3 to 1000 (sorted, unique, ending at π/Ts, containing the equivalent pole frequency), roots at z = 0 and z = 1, resonance peaks and anti-resonances within 1 % for ζ = 1e-2, 1e-3, 1e-4 and ωn = 1e-3, 1, 1e3, a discrete-time resonance, an undamped resonance, and an unchanged grid for a well-damped system.
• test_analysis.jl: the margin example L above; the margin throws an error when the system has delay #941 test asserted imag(L(iω)) ≈ 0 with the absolute tolerance 1e-6 at the gain-margin frequency. The smallest gain margin is now located at the resonance, where |L| = 8.7, and the crossover is interpolated, so the tolerance is made relative (1e-3|L|) and the gain margin and its frequency are compared with dense-grid reference values.
• test_delayed_systems.jl: grid bounds for delays 1e-3, 1e3, 1e5.
• test_plots.jl: downsample retains the resonance peak of a non-equidistant grid and terminates with zero tolerance.

Results with Julia 1.12.7 on Linux:

• test_freqresp.jl: 415/415 passed. The allocation test at line 69 passes both with this PR and on the head of Fix default margin frequency ranges and missing phase guides #1076 (102/102), so the allocation-limit failure reported in Fix default margin frequency ranges and missing phase guides #1076 (Julia 1.10 on Windows) does not occur in this environment.
• test_analysis.jl: 248/248 passed.
• test_delayed_systems.jl: 147/147 passed.
• test_hammerstein_wiener.jl: 70/70 passed.
• test_pid_design.jl: 152/152 passed.
• test_discrete.jl: 60/60 passed.
• test_plots.jl: 34/34 passed.
• test_makie_plots.jl: 103/103 passed.

The full package test suite was not run.

---

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

🤖 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 92.22%. Comparing base (ef5f104) to head (e782e2a).

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #1086      +/-   ##
==========================================
+ Coverage   92.17%   92.22%   +0.05%     
==========================================
  Files          42       42              
  Lines        6119     6146      +27     
==========================================
+ Hits         5640     5668      +28     
+ Misses        479      478       -1     

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

[JuliaControlBot]

This is an automated message.
Plots were compared to references. 5/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.0
✔️ 0.013