Compute the default simulation time vector from the decay rates of the non-integrator poles

[baggepinnen]

Part of #1081.

Problem

The default time vector of step, impulse and lsim(sys, u::Function, tfinal) (_default_time_vector, _default_dt in timeresp.jl) had the following defects.

1. Integrators were detected with an exact comparison iszero. A multiple pole in the origin is computed with rounding errors of order eps^(1/k)·‖A‖, and the sample interval became very large:

   P = similarity_transform(ss(tf(1, [1, 0, 0])), [1 2; 3 4.5])
   ControlSystemsBase._default_dt(P)  # 4.1e14, tfinal = 8.2e16, step(P).y contains NaN
   P = minreal(ss(tf([1, 1], [1, 1, 0, 0])))  # poles ±2.9e-9im
   ControlSystemsBase._default_dt(P)  # 2.8e7

   A triple integrator in transformed coordinates gave dt ≈ 1e4. The absolute threshold 1e-6 used for the final time was inconsistent with this test.
2. The final time. (a) minimum(w for w in ws if w > 1e-6; init=dt) included a time (dt) in a minimum over frequencies, so the result depended on the time unit: the step response of 1/((s+1)(s+0.01)) ended at y = 0.387 with the time unit s and at y = 1.0 with the time unit ms. (b) 1/(2ω_min) is half of the slowest time constant: for the poles -1 and -1e-3, the step response ended at y(500) = 0.39. (c) The threshold 1e-6 was absolute. (d) |p| was used instead of the decay rate |Re p|: for ζ = 0.01, the envelope was still 0.85 at tfinal.
3. Discrete-time systems. The magnitudes |z| ≈ 1 were treated as frequencies in rad/s, so that tfinal depended only on Ts: step(c2d(ss(tf(1, [100, 1])), Ts)) ended at y = 0.39 for Ts = 0.01, y = 0.86 for Ts = 1 and y = 0.095 for Ts = 1e-4.
4. The limit of 100_000 samples was enforced by silently reducing tfinal.

Change

• _nonintegrator_poles(sys): poles in the origin (z = 1 in discrete time) are classified as integrators with count_eigval_multiplicity, i.e., the tolerance (100ϵ)^(1/k)·scale for a pole of multiplicity k, as in count_integrators/integrator_excess. The scale includes opnorm(A, 1), since the rounding errors of the eigenvalues are proportional to the norm of A, which is much larger than the magnitude of the computed poles if all poles are in the origin. Discrete-time poles are mapped by s = log(z)/Ts; a pole at z = 0 is mapped to s = -∞ and does not affect the final time.
• _default_dt(sys) (continuous time): 1/(12 max|p|) over the non-integrator poles, as before. If all poles are integrators or the system has no poles, dt = 0.05 as before. A system whose poles are all in the origin, such as 1/s^k, has no characteristic time (it is invariant under time scaling up to a gain). I did not use opnorm(A) as a time scale for this case, since the norm of a nilpotent A depends on the scaling of the state variables ([0 c; 0 0] is similar to [0 1; 0 0] for every c ≠ 0).
• _default_tfinal(sys, dt):
  • Stable systems: tfinal = μ + 6s with μ = Σ τᵢ, s = √(Σ τᵢ²) and τᵢ = 1/|Re pᵢ| over the stable non-integrator poles. The impulse response of a series connection of first-order systems with time constants τᵢ equals the probability density of a sum of independent exponentially distributed random variables with means τᵢ, whose mean and standard deviation are μ and s. For a single pole, tfinal = 7τ, after which exp(-7) < 1e-3 of the initial deviation remains; for a pole of multiplicity m, the term 6√m τ accounts for the polynomial factors tᵐ⁻¹e^{-t/τ} (for 1/(s+1)^5, tfinal = 18). Each pole of a complex-conjugate pair is counted, which yields ≈ 10.5τ for a single pair; counting each pair once would make tfinal depend discontinuously on the poles, since rounding errors split a multiple real pole into complex-conjugate pairs (this was observed for (s+1)^5 after a conversion to a transfer function). The time constants are bounded from below by dt, since a discrete-time mode does not decay faster than in one sample interval.
  • Undamped oscillatory poles (|Re p| below the tolerance 100ϵ·scale of a simple eigenvalue): at least five periods of the slowest such mode.
  • Unstable systems: five time constants of the fastest growing mode, i.e., the growth by a factor of about e^5 is visible but the response does not overflow (previously e^16.6 for 1/(s-1)).
  • Integrators only, or no poles: 200dt, as before.
  • tfinal is rounded to two significant digits.
• Limit on the number of samples (100_001): for a continuous-time system, dt is increased (the simulation of a continuous-time system uses a ZoH discretization, so the samples remain exact, only the resolution of the plot is reduced); for a discrete-time system, whose sample interval is fixed, tfinal is reduced and a warning is emitted.
• The result is invariant under ControlSystemsBase.time_scale up to the rounding to two significant digits, and for discrete-time systems it does not depend on Ts (unless the limit on the number of samples is reached).
• The docstrings of step and impulse describe the rule.

DelayLtiSystem and HammersteinWienerSystem do not define poles, they keep tfinal = 200dt with their own _default_dt (HammersteinWienerSystem uses the new _default_dt of its linear part). _default_time_data in plotting.jl has no callers and is unchanged.

Effect on users

The default time vector changes for most systems. Systems with only integrators and static systems are unchanged. Representative systems (y(tfinal) is the last value of the step response, the DC gain is 1 for all stable systems below; dcgain of 1/(s-1) is -1):

system	tfinal before	dt before	y(tfinal) before	tfinal after	dt after	y(tfinal) after
1/(s+1)	16.6	0.083	1.0	6.97	0.083	0.999
1/(s+1), time unit ms (time_scale(·, 1e3))	8.3	8.3e-5	1.0	0.00697	8.3e-5	0.999
1/(s+1)^5	16.6	0.083	1.0	17.9	0.083	1.0
1/(s²+0.02s+1) (ζ = 0.01)	16.6	0.083	1.54	1000	0.083	1.0
1/(s²+2e-4s+1) (ζ = 1e-4)	16.6	0.083	1.63	1e5	1.0	1.0
1/((s+1)(s+1e-3))	500	0.083	0.393	7000	0.083	0.999
1/((s+1)(s+0.01))	50	0.083	0.387	700	0.083	0.999
1/((s+1)(s+0.01)), time unit ms	8.3	8.3e-5	1.0	0.7	8.3e-5	0.999
1/s	10	0.05	10	10	0.05	10
1/(s(s+1))	16.6	0.083	15.6	6.97	0.083	5.97
1/s², transformed coordinates	8.2e16	4.1e14	NaN	10	0.05	50
minreal(ss((s+1)/(s²(s+1))))	5.6e9	2.8e7	NaN	10	0.05	50
1/s³, transformed coordinates	2.4e6	1.2e4	NaN	10	0.05	167
1/(s²+1)	16.6	0.083	1.63	31	0.083	0.103
1/(s-1)	16.6	0.083	1.6e7	4.98	0.083	144
1/((s+0.5)(s-1))	16.6	0.083	1.1e7	4.98	0.083	95
static gain	10	0.05	1	10	0.05	1
1/(100s+1)	1660	8.3	1.0	697	8.3	0.999
c2d(1/(100s+1), 1e-4)	10	1e-4	0.095	10 (warning)	1e-4	0.095
c2d(1/(100s+1), 0.01)	50	0.01	0.393	700	0.01	0.999
c2d(1/(100s+1), 1)	200	1	0.865	700	1	0.999
c2d(1/(100s+1), 10)	2000	10	1.0	700	10	0.999
1/(z-1), Ts = 0.1	20	0.1	200	20	0.1	200
1/z, Ts = 0.1	20	0.1	1	0.7	0.1	1

Notable changes:

• Simple stable systems are simulated for a shorter time than before (7 instead of 16.6 time constants for a first-order system, with 85 instead of 201 samples), systems with slow or lightly damped modes for a longer time.
• For ζ = 1e-4, the limit on the number of samples increases dt to 1.0, i.e., about six samples per period of the oscillation.
• Deadbeat discrete-time systems (poles at z = 0) are simulated for a few samples only.
• step of a discrete-time system whose slowest mode requires more than 100_000 samples now emits a warning.

Tests

Added a testset default time vector to lib/ControlSystemsBase/test/test_timeresp.jl: the reproducers above (multiple integrators in transformed coordinates), the final value of the step response compared with the DC gain for stable systems (first order, ζ = 0.01, widely separated poles, a pole of multiplicity five, a third-order system) in three time units, time-scale invariance of tfinal and dt, independence of Ts for a discrete-time system, the warning for a discrete-time system that exceeds the limit on the number of samples, the coarsening of dt for a continuous-time system, unstable systems in three time units, undamped oscillatory systems, and user-provided tfinal. The existing tests of _default_dt and _default_time_vector pass unchanged.

Results (Julia 1.12):

• lib/ControlSystemsBase/test: test_timeresp, test_discrete, test_delayed_systems, test_hammerstein_wiener: 457/457 passed.
• lib/ControlSystemsBase/test: test_simplification, test_plots, test_staticsystems, test_analysis, test_demo_systems, test_complex, test_partitioned_statespace, test_dsp, test_implicit_diff, test_rootlocus: 3494/3494 passed.
• test/runtests_controlsystems.jl (ControlSystems: test_timeresp, including the test of the length of the default time vector, test_delay_timeresp, test_nonlinear_timeresp): 1151/1151 passed.

---

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

🤖 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.76%. Comparing base (59c065a) to head (7f60d0c).

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #1085      +/-   ##
==========================================
+ Coverage   91.73%   91.76%   +0.03%     
==========================================
  Files          42       42              
  Lines        5865     5892      +27     
==========================================
+ Hits         5380     5407      +27     
  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.