Defaults and tolerances that depend on the time unit, gain, sample time or element type #1081

An audit of the default values and tolerances in ControlSystemsBase and ControlSystems (at 59c065a) identified defaults whose result changes when the time unit, the sample time, the system gain, the state coordinates or the element type is changed. Such defaults should instead be derived from quantities that transform consistently, e.g.,

• a characteristic frequency computed from the poles and zeros (for discrete-time systems after mapping the poles with s = log(z)/Ts), as done for the root locus in Root locus: relative step-size control, improper systems, default gain and plot limits #1079,
• tolerances relative to eps(T) of the element type T, scaled by opnorm(A) or maximum(abs, poles), with the classification of count_eigval_multiplicity wherever a test decides whether a pole is located at the origin, on the imaginary axis or on the unit circle,
• gain crossovers determined from the gain of the system rather than from pole and zero locations alone.

The findings that produce wrong results for realistic systems are addressed in separate pull requests (listed below). This issue collects the remaining findings. Each finding was reproduced experimentally unless noted otherwise; line numbers refer to 59c065a.

Addressed in pull requests

• The discrete-time unstable-pole count in sisomargin (analysis.jl:609), and the margin frequency grid not containing the gain crossover: Fix default margin frequency ranges and missing phase guides #1076
• isunstable for discrete-time systems, the absolute atol of baltrunc, and the exact comparisons in delayd_ss and thiran: Replace exact and absolute comparisons in isunstable, baltrunc, delayd_ss and thiran #1083
• The default frequency vector (absolute cutoff at 1e-4 rad/s, z-plane magnitudes used as frequencies, lightly damped resonances, delays): Improve the default frequency vector: features below 1e-4 rad/s, discrete-time systems, lightly damped resonances and delays #1086 (builds on Fix default margin frequency ranges and missing phase guides #1076)
• The default time vector of step, impulse and lsim (exact integrator test, final time, discrete-time systems): Compute the default simulation time vector from the decay rates of the non-integrator poles #1085
• The non-stiff default solver and the missing return-code check of the ODE-based lsim: Use a stiffness-switching default solver in lsim and check the solver return code #1082
• The default axis limits of nyquistplot (Plots and Makie): Compute the default axis limits of nyquistplot per subplot from the curves and the critical point #1084

Tolerances and stability classification

• minreal for tf/zpk (types/TransferFunction.jl:83, SisoTfTypes/SisoRational.jl:41, SisoTfTypes/SisoZpk.jl:221,247): pole-zero cancellation uses the absolute distance sqrt(eps()) with Float64 eps regardless of the element type. minreal(tf([1,1.5e-8],[1,1e-8])) changes the DC gain from 1.5 to 1, zpk([],[-1e-8],1e-8) + 1 has DC gain 1 instead of 2, f+f with f = zpk([-2e6],[-1e6,-3e6,-4e6],1e12) retains five poles instead of three, and for Float32 128/200 random products with a common factor are not cancelled. Proposal: |p - z| ≤ tol*max(|p|, |z|, ρ) with tol = sqrt(eps(real(float(T)))) and a floor ρ relative to the largest root magnitude.
• isapprox for tf/zpk (SisoZpk.jl:216, SisoRational.jl:94): the zpk method uses an absolute atol = sqrt(eps()) on the gain (isapprox(1e-9H, 2e-9H) is true), the tf method compares raw coefficient vectors, which is dominated by the highest powers of s (second-order systems with ω = 1e9 and damping 0.1 vs. 0.2 compare equal). Proposal: normalize s by ω0 = |a₀/aₙ|^(1/n), compare gains and sorted roots relatively, use eps(T).
• isstable (types/Lti.jl:114-115): the exact tests real(p) < 0 and abs(p) < 1 classify an integrating system as stable when the pole at the origin is computed as -1e-17; this happened in 102/200 random similarity transforms of 1/(s(s²+s+1)). Proposal: real(p) < -tol and abs(p) < 1 - tol with the tolerance of count_eigval_multiplicity. Note that this changes the result for integrating systems that are currently classified as stable by rounding error, and callers such as hinfnorm and covar.
• isunstable, continuous time (types/Lti.jl:131-133): the imaginary-axis prefilter abs(real(p)) > sqrt(sqrt(eps(abs(p)))) scales as |p|^(1/4), and the multiplicity is counted around the projection of a single pole of a cluster. A double pole pair on the imaginary axis at ω = 1e6 is missed in 21/50 random coordinate transforms. Two parallel integrators (semisimple, marginally stable) are reported as unstable. Proposal: cluster with the radius of count_eigval_multiplicity, test the cluster centroid, and decide the Jordan structure from rank(A - λI) for state-space input.
• hinfnorm/linfnorm: in discrete time approxcirc = 1e-8 is an absolute bound on ||z| - 1| (matrix_comps.jl:545,549), so hinfnorm(c2d(ss(tf(1,[1,1])), Ts)) is Inf for Ts ≤ 1e-9; in continuous time the floor in approximag*max(1, |z|) (matrix_comps.jl:438) has the unit of frequency, so hinfnorm(ss(tf(1e-11,[1,1e-11]))) is Inf. Proposal: test in the s-domain equivalent log(z) and replace the floor by a magnitude relative to opnorm(A) or the largest pole magnitude.
• place (synthesis.jl:345): @assert all(abs(imag(l)) .< 1e-6) is an absolute test on a dimensioned gain; with B scaled by 1e-9 it fails in 38/200 random cases although the relative imaginary part is below 7e-16. Proposal: abs(imag(l)) <= sqrt(eps(real(T)))*norm(L).
• place_knvd (synthesis.jl:420,492): the rank of B is determined with Float64 eps() for every element type (a rank-deficient Float32 B is treated as having full rank and the assigned eigenvalues are wrong), and the stopping test cond_old - c < 1e-14 is absolute. Proposal: max(size(B)...)*eps(real(T))*R[1] and a relative stopping test.
• relative_gain_array (analysis.jl:466-468): the branch is selected by the exact test det(A) != 0, so relative_gain_array(A) and relative_gain_array(1e-3A) differ by 16 orders of magnitude for a numerically singular A, contradicting the documented scale invariance. tol = 1e-15 is absolute, the iteration has no limit, and the work arrays are Float64. Proposal: select the branch by a rank test relative to the largest singular value, or always use the iterative branch.
• controllability/observability (matrix_comps.jl:185,209): the rank tolerance relative to σ_max([λI-A B]) mixes the units of A and B, so the result changes under diagonal state scaling. This is inherent to a binary rank test; documenting it may suffice.
• _findCrossings (analysis.jl:655-657): grid points at which the phase is exactly -180° are counted as crossings, so margin(tf(1,[1,0,0])) reports 200 phase crossovers and the gain margin equals 1/|L| at the first grid point, which depends on the time unit. filter_th is passed but not used.
• sisomargin (analysis.jl:563): phase crossovers with abs(Giw) > 1e6 && wgm < 0.001 are removed, an absolute test in gain and frequency that is active only for the default grid. For L = 1e6(s+1)²/s³, time_scale(L, 5e-4) reports gm = Inf instead of 5e-7.
• convert_pidparams_* to :series (pid_design.jl:810,839): the exact test Δ < 0 throws a DomainError in 32/1000 round trips with Ti = 4Td due to rounding. Proposal: max(Δ, 0) for Δ ≥ -8eps(T)*Ti².
• ss · improper tf warning (types/Lti.jl:23): the warning threshold 1e7 on matrix entries is absolute.
• common_timeevol (types/TimeEvolution.jl:42): sample times are compared exactly, so ss(1,1,1,1,0.1)*ss(1,1,1,1,0.3/3) throws.

Time and frequency grids

• Time vector of delay systems (delay_systems.jl:88-99, timeresp.jl:439): the default ignores slow poles and the delay lengths; for DemoSystems.woodberry() the final time is 16.6 s while the time constants are 10.9-21 s and the delays up to 7 s, and for an unstable P the final time is 10 regardless of the time unit. Proposal: τ_max + k/min|Re p| over the poles of P.P.
• Frequency grid of delay systems (delay_systems.jl:71-77): log10(τ) is used as a frequency exponent, and only delays τ > 1e4 are considered. Proposal: use the feature frequency 1/τ.
• pidplots (pid_design.jl:270) uses exp10.(range(-3, 3, 500)) and stabregionPID(::Function) (pid_design.jl:447) exp10.(range(-3, 1, 50)). Proposal: _default_freq_vector of the loop transfer functions for pidplots.
• describing_function_plot (nonlinear_components.jl:373): A_range = range(0.1, 10, length=200) is absolute and linear; for Saturation(100) the curve is a single point and for DeadZone(100) no point is finite. Proposal: logarithmic spacing relative to the parameter of the nonlinearity.
• Default grid for BigFloat and MIMO systems (freqresp.jl:522): bode without a frequency vector fails for BigFloat (via tzeros), and the per-channel zpkdata is expensive for MIMO systems. Only the order of magnitude of the features is needed, which can be computed from Float64.(ssdata(sys)) with poles and transmission zeros. The default grid is always Float64.
• Inconsistent grid options: sigmaplot, rgaplot and nicholsplot use the non-adaptive grid (minimum 200 points, 60 per decade), whereas bodeplot, marginplot and nyquistplot use the adaptive grid (minimum 1000 points, 100 per decade). _default_time_data (plotting.jl:899) uses 100dt versus 200dt in _default_time_vector and is unused.

Simulation

• ODE tolerances (src/timeresp.jl:100,336): abstol = reltol = 1e-6 for delay and Hammerstein-Wiener systems is absolute in the signal amplitude; for a gain of 1e-9 the relative output error is 0.23 (delay system) and 0.99 (Hammerstein-Wiener system). Rational lsim uses the solver defaults with dtmax = dt. Proposal: an absolute tolerance relative to an a-priori state scale, and consistent defaults across the methods (delay impulse uses BS3, delay step Tsit5).
• hysteresis (nonlinear_components.jl:189): hardness = 20.0 multiplies the input signal and is not relative to width; the element is bistable only if width*hardness > 1, so for width ≤ 0.04 the hysteresis silently disappears. Tf = 0.001 is an absolute time, which makes the simulation of a slow loop stiff (2.45 s vs. 0.0023 s with Tf = 1 for a plant with time constant 100 s). Proposal: tanh(hardness*y/width) and a required Tf or one derived from the plant.
• dab (discrete.jl:442): println whenever cond(m) > 1e6, an absolute threshold that depends on the scaling of the polynomial coefficients; the early returns at lines 421/430 return nothing.
• Element types: c2d with :zoh/:foh returns Float64 matrices for a Float32 system (discrete.jl:50,57) and fails for BigFloat (exp!), so step and lsim fail for BigFloat; lsim converts the time step with Float64(t[2]-t[1]) (timeresp.jl:194). c2d_x0map returns nx×nx for tustin/euler and nx×(nx+nu) for zoh/foh.

Plotting

• gangoffourplot (plotting.jl:986-997): S and T are plotted with ylims=(1e-2,4), which hides a sensitivity peak of 40.8; under setPlotScale("dB") the limits and the Ms/Mt lines are not converted to dB. Proposal: an upper limit of max(4, 1.1max(Ms, Mt)) and conversion of limits and lines in dB mode.
• Nyquist downsampling (plotting.jl:451-452, ControlSystemsBaseMakieExt.jl:298-299): the tolerance 1/500 is absolute, so of 1000 points 8 are retained for a gain of 1e-3 and 394 for a gain of 100. Proposal: a tolerance relative to the extent of the curve, as in the Bode recipe.
• Makie bodeplot/marginplot downsampling (ControlSystemsBaseMakieExt.jl:194,450): the tolerance is computed on the linear magnitude, whereas Plots uses the log magnitude; for a magnitude spanning 1e-12…0.5 the maximum error is 0.54 decades (Plots: 0.02).
• Makie pzmap (ControlSystemsBaseMakieExt.jl:85,109): the ωn labels use round(ωn, digits=1), so all labels read ωn=0.0 for a system in milliseconds, and the ωn circles are spaced linearly.
• pzmap(...; hz=true) for discrete-time systems (plotting.jl:927-929, Makie 48-70): the z-plane poles are divided by 2π (z = 0.905 is plotted at 0.144); Makie also scales the unit circle.
• nicholsplot (plotting.jl:588,654): w_nyquist = 2π/Ts is twice the Nyquist frequency, and the phase lines are drawn over a fixed range of [-110, 25] dB.
• leadlinkcurve (pid_design.jl:391, ControlSystemsBaseMakieExt.jl:617): 180/pi*atan(sqrt(N)) - atan(1/sqrt(N)) converts only the first term to degrees; for N = 1 the phase advance is shown as 44.2° instead of 0°, for N = 4 as 62.97° instead of 36.87°. M > 1 || @warn "M should be ≥ 1" warns at M = 1.

Other defects found during the audit

• diagonalize (types/StateSpace.jl:631) calls ss(A,B,C,D) and returns a continuous-time system for a discrete-time input; the docstring mentions an unused digits=12.
• Filter parameterization: loopshapingPI/PD/PID use the filter 1/(Tf²s² + √2Tf s + 1), whereas pid(...; Tf) uses 1/(Tf²/2 s² + Tf s + 1) for the default damping, so the same keyword gives filter bandwidths differing by a factor √2.
• Unused utilities (utilities.jl:41,98): issemiposdef and _fix_conjugate_pairs! (Float64-specific rtol = 1e-15) are not called.
• Implicit differentiation extension (ControlSystemsBaseImplicitDifferentiationExt.jl:201-229): hinfnorm and are for dual numbers discard their keyword arguments.

Host: demeter6, Claude Code session session_01HigNnun4aJKMoWmY7wyULn, working directory ~/.julia/dev/ControlSystems/.claude/worktrees/generic-pondering-babbage

🤖 Generated with Claude Code

https://claude.ai/code/session_01HigNnun4aJKMoWmY7wyULn