From 7f60d0c036dfff2c12cf374bd8c4faa0a6c96b3a Mon Sep 17 00:00:00 2001 From: Fredrik Bagge Carlson Date: Fri, 9 Oct 2026 04:58:15 +0000 Subject: [PATCH] Compute the default simulation time vector from the decay rates of the non-integrator poles The default sample interval and final time of step, impulse and lsim are computed from the poles of the system as follows: - Poles in the origin (z = 1 in discrete time) are classified as integrators with the multiplicity-dependent tolerance of count_eigval_multiplicity, scaled by the norm of A. Previously, an exact comparison with zero was used, and a multiple pole in the origin computed with rounding errors yielded a sample interval of up to 4e14. - Discrete-time poles are mapped to the s-plane by s = log(z)/Ts. Previously, the magnitudes |z| were treated as frequencies, and the final time depended only on Ts. - The final time is computed from the decay rates |Re p| of the stable poles as mu + 6s, where mu and s are the sum of the time constants and the square root of the sum of their squares, which equals seven time constants for a single pole. Previously, the final time was half the slowest time constant, computed from |p| instead of |Re p|, and depended on the time unit through a minimum over both times and frequencies. - Unstable systems are simulated for five time constants of the fastest growing mode, undamped oscillatory modes for at least five periods. - The limit on the number of samples is enforced by increasing the sample interval of a continuous-time system, and by reducing the final time with a warning for a discrete-time system, instead of silently reducing the final time. The resulting time vector is invariant under time_scale up to rounding. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01HigNnun4aJKMoWmY7wyULn --- lib/ControlSystemsBase/src/timeresp.jl | 110 +++++++++++++++---- lib/ControlSystemsBase/test/test_timeresp.jl | 78 +++++++++++++ 2 files changed, 168 insertions(+), 20 deletions(-) diff --git a/lib/ControlSystemsBase/src/timeresp.jl b/lib/ControlSystemsBase/src/timeresp.jl index 806e9a8ac..2af1e41a5 100644 --- a/lib/ControlSystemsBase/src/timeresp.jl +++ b/lib/ControlSystemsBase/src/timeresp.jl @@ -39,7 +39,11 @@ LsimWorkspace(sys::AbstractStateSpace, u::AbstractMatrix) = LsimWorkspace(sys, s Calculate the response of the system `sys` to a unit step at time `t = 0`. If the final time `tfinal` or time vector `t` is not provided, -one is calculated based on the system pole locations. +one is calculated based on the system pole locations: the sample interval +resolves the fastest pole and the final time is chosen such that the slowest +stable mode has decayed, or, for an unstable system, such that the fastest +growing mode has increased by a factor of about ``e^5``. Integrators are +excluded from this computation. The return value is a structure of type `SimResult`. A `SimResul` can be plotted by `plot(result)`, @@ -81,7 +85,11 @@ Calculate the response of the system `sys` to an impulse at time `t = 0`. For continous-time systems, the impulse is a unit Dirac impulse. For discrete-time systems, the impulse lasts one sample and has magnitude `1/Ts`. If the final time `tfinal` or time vector `t` is not provided, -one is calculated based on the system pole locations. +one is calculated based on the system pole locations: the sample interval +resolves the fastest pole and the final time is chosen such that the slowest +stable mode has decayed, or, for an unstable system, such that the fastest +growing mode has increased by a factor of about ``e^5``. Integrators are +excluded from this computation. The return value is a structure of type `SimResult`. A `SimResul` can be plotted by `plot(result)`, @@ -422,20 +430,30 @@ end # HELPERS: -# TODO: This is a poor heuristic to estimate a "good" time vector to use for -# simulation, in cases when one isn't provided. +""" + t = _default_time_vector(sys, tfinal = -1) + +Compute a time vector for simulation of `sys` when none is provided. + +The sample interval `dt` resolves the fastest pole, see `_default_dt`. If `tfinal` is not provided, it is computed from the poles by `_default_tfinal`, such that the slowest stable mode has decayed. The number of samples is limited to 100 001: for a continuous-time system, `dt` is increased if required, for a discrete-time system, whose sample interval is fixed, `tfinal` is reduced and a warning is emitted. + +The time vector is invariant under `time_scale` up to the rounding of `dt` and `tfinal` to two significant digits. +""" function _default_time_vector(sys::LTISystem, tfinal::Real=-1) - dt = _default_dt(sys) # This is set small enough to resolve the fastest dynamics. + dt = _default_dt(sys) # This is set small enough to resolve the fastest dynamics. if tfinal == -1 - # The final simulation time should be chosen to also let the slowest dynamics run its course. - # Integrating systems pose a problem for this heuristic, and we sort out any poles that appear to be integrators. if hasmethod(poles, typeof((sys,))) # DelaySystem does not define poles - ws = abs.(poles(sys)) - ω0_min = minimum(w for w in ws if w > 1e-6; init=dt) - dt_slow = round(1/(2ω0_min), sigdigits=2) - tfinal = max(200dt, dt_slow) - tfinal = min(tfinal, 100_000*dt) + tfinal = _default_tfinal(sys, dt) + nmax = 100_000 + if tfinal > nmax*dt + if isdiscrete(sys) + @warn "The default final time $tfinal of the simulation of the discrete-time system requires more than $nmax samples, it is reduced to $(nmax*dt). Provide the final time or a time vector to simulate for a longer time." + tfinal = nmax*dt + else + dt = round(tfinal/nmax, RoundUp, sigdigits=2) + end + end else tfinal = 200dt end @@ -445,16 +463,68 @@ function _default_time_vector(sys::LTISystem, tfinal::Real=-1) return 0:dt:tfinal end -function _default_dt(sys::LTISystem) +""" + p, tol = _nonintegrator_poles(sys) + +Return the poles of `sys` that are not integrators, mapped to the s-plane by s = log(z)/Ts for discrete-time systems, together with a tolerance `tol` below which the real part of such a pole is considered zero. A pole at z = 0 is mapped to s = -∞. + +Poles in the origin (z = 1 in discrete time) are classified as integrators with the tolerance of `count_eigval_multiplicity`, which accounts for the larger rounding errors of multiple poles. The rounding errors of the eigenvalues of `A` are proportional to the norm of `A`, which may be much larger than the magnitude of the computed poles, e.g., if all poles are in the origin. The norm of `A` is therefore included in the scale of the tolerance. +""" +function _nonintegrator_poles(sys::LTISystem) + p = poles(sys) + location = iscontinuous(sys) ? 0 : 1 + scale = float(maximum(abs, p, init=0.0)) + if sys isa AbstractStateSpace && !isempty(sys.A) + scale = max(scale, opnorm(sys.A, 1)) + end + nint, tol = count_eigval_multiplicity(p, location; scale) + if nint > 0 + p = filter(p -> abs(p - location) > tol, p) + end + tol_axis = 100*eps(float(real(eltype(p))))*scale # Tolerance for a simple pole on the imaginary axis (unit circle) if isdiscrete(sys) - return sys.Ts - elseif all(iszero, poles(sys)) # Static or pure integrators - return 0.05 - else - ω0_max = maximum(abs.(poles(sys))) - dt = round(1/(12*ω0_max), sigdigits=2) - return dt + return log.(complex.(p)) ./ sys.Ts, tol_axis/sys.Ts end + p, tol_axis +end + +""" + tfinal = _default_tfinal(sys, dt) + +Final time of the default simulation time vector with sample interval `dt`, computed from the poles that are not integrators: +- If there are unstable poles, `tfinal = 5τ` where `τ` is the time constant of the fastest growing mode, such that the growth is visible but the response does not overflow. +- Otherwise, `tfinal = μ + 6s` with `μ = Σ τᵢ` and `s = √(Σ τᵢ²)`, where `τᵢ` are the time constants of the stable 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 a fraction `exp(-7) < 1e-3` of the initial deviation remains. For poles of high multiplicity, `s` accounts for the polynomial factors of the response. The time constant of a complex pole is the inverse of the magnitude of its real part, i.e., the time constant of the envelope of the oscillation. Each pole of a complex-conjugate pair is counted, which yields `tfinal ≈ 10.5τ` for a single pair. Counting each pair once would yield a final time that depends discontinuously on the poles, since rounding errors split a multiple real pole into complex-conjugate pairs. The time constants are bounded from below by `dt`, since a discrete-time mode does not decay faster than in one sample interval. +- `tfinal` is at least five periods of the slowest undamped oscillatory mode. +- If all poles are integrators, or the system has no poles, the system has no characteristic time and `tfinal = 200dt`. +""" +function _default_tfinal(sys::LTISystem, dt) + p, tol = _nonintegrator_poles(sys) + isempty(p) && return 200dt + σ = -real.(p) # Decay rates + growth_rate = -minimum(σ) + if growth_rate > tol + return round(5*max(1/growth_rate, dt), sigdigits=2) + end + stable = σ .> tol + τ = [max(1/σ[i], dt) for i in eachindex(p) if stable[i]] + tfinal = sum(τ, init=0.0) + 6*sqrt(sum(abs2, τ, init=0.0)) + ω_undamped = minimum((abs(p[i]) for i in eachindex(p) if !stable[i]), init=Inf) + if isfinite(ω_undamped) + tfinal = max(tfinal, 5*2π/ω_undamped) + end + round(tfinal, sigdigits=2) +end + +""" + dt = _default_dt(sys) + +Default sample interval for simulation of `sys`. For a discrete-time system, `dt = sys.Ts`. For a continuous-time system, `dt = 1/(12 max|pᵢ|)` resolves the fastest pole, where integrators are excluded. If all poles are integrators, or the system has no poles, the system has no characteristic time and `dt = 0.05`. +""" +function _default_dt(sys::LTISystem) + isdiscrete(sys) && return sys.Ts + p, _ = _nonintegrator_poles(sys) + isempty(p) && return 0.05 + round(1/(12*maximum(abs, p)), sigdigits=2) end """ diff --git a/lib/ControlSystemsBase/test/test_timeresp.jl b/lib/ControlSystemsBase/test/test_timeresp.jl index fdac5ea69..6e69ec182 100644 --- a/lib/ControlSystemsBase/test/test_timeresp.jl +++ b/lib/ControlSystemsBase/test/test_timeresp.jl @@ -204,6 +204,84 @@ sys2tau = tf(1, [0.1, 1]) * tf(1, [10, 1]) @test 5 < ControlSystemsBase._default_time_vector(sysint)[end] < 15 @test 20 < ControlSystemsBase._default_time_vector(sys2tau)[end] < 100 +@testset "default time vector" begin + default_tv = ControlSystemsBase._default_time_vector + default_dt = ControlSystemsBase._default_dt + + # Multiple poles in the origin are computed with rounding errors, they are classified as integrators with a tolerance that depends on the multiplicity + P2 = similarity_transform(ss(tf(1, [1, 0, 0])), [1 2; 3 4.5]) + P2m = minreal(ss(tf([1, 1], [1, 1, 0, 0]))) + P3 = similarity_transform(ss(tf(1, [1, 0, 0, 0])), [1 2 0; 3 4.5 1; 0 1 1]) + for P_tv in (P2, P2m, P3) + @test default_dt(P_tv) == 0.05 + @test default_tv(P_tv)[end] ≈ 10 + @test all(isfinite, step(P_tv).y) + end + @test default_tv(tf(1, [1, 0]))[end] ≈ 10 + @test default_tv(ss(1.0))[end] ≈ 10 + + # The step response of a stable system reaches the DC gain, also when the time constants are widely separated or the damping is low + stable_systems = [ + tf(1, [1, 1]), + tf(1, [1, 2*0.01, 1]), + zpk([], [-1, -1e-3], 1e-3), + zpk([], [-1, -0.01], 0.01), + zpk([], fill(-1.0, 5), 1), + tf(1, [1, 0.5]) * tf(4, [1, 0.1, 4]), + ] + for sys_tv in stable_systems, c_tv in (1, 1e-3, 1e3) + sys_tv_c = ControlSystemsBase.time_scale(sys_tv, c_tv) + @test step(sys_tv_c).y[end] ≈ dcgain(sys_tv_c)[] rtol = 1e-2 + # Time-scale invariance + t_tv1 = default_tv(sys_tv) + t_tv2 = default_tv(sys_tv_c) + @test t_tv2[end] ≈ t_tv1[end]/c_tv rtol = 1e-2 + @test step(t_tv2) ≈ step(t_tv1)/c_tv rtol = 1e-2 + end + # The envelope of a lightly damped mode has decayed + @test default_tv(tf(1, [1, 2*0.01, 1]))[end] >= 7/0.01 + + # Discrete-time systems: the final time is determined by the poles mapped to the s-plane, it does not depend on the sample time + G_tv = ss(tf(1, [100, 1])) + for Ts_tv in (0.01, 1, 10) + G_tv_d = c2d(G_tv, Ts_tv) + t_tv = default_tv(G_tv_d) + @test step(t_tv) == Ts_tv + @test t_tv[end] ≈ 700 rtol = 0.02 + @test step(G_tv_d).y[end] ≈ 1 rtol = 1e-2 + end + # The number of samples is limited, a warning is emitted when the final time is reduced + @test_logs (:warn, r"reduced") default_tv(c2d(G_tv, 1e-4)) + @test length(default_tv(c2d(G_tv, 1e-4))) == 100_001 + # Integrators and poles at the origin of the z-plane in discrete time + @test default_tv(tf(1, [1, -1], 0.1))[end] ≈ 20 + @test default_tv(tf(1, [1, 0, 0], 0.1))[end] >= 0.2 + @test step(tf(1, [1, 0, 0], 0.1)).y[end] == 1 + + # For a continuous-time system, the sample interval is increased to limit the number of samples + t_tv = default_tv(tf(1, [1, 2e-4, 1])) + @test length(t_tv) <= 100_001 + @test t_tv[end] >= 7/1e-4 + + # Unstable systems are simulated for a few time constants of the fastest growing mode + for c_tv in (1, 1e-3, 1e3) + sys_tv = ControlSystemsBase.time_scale(tf(1, [1, 0.5]) * tf(1, [1, -1]), c_tv) + t_tv = default_tv(sys_tv) + @test 4/c_tv < t_tv[end] <= 5/c_tv + @test all(isfinite, step(sys_tv).y) + end + @test 0.1 <= default_tv(tf(1, [1, -2], 0.1))[end] <= 1 + + # Undamped oscillatory modes are simulated for several periods + @test default_tv(tf(1, [1, 0, 1]))[end] >= 0.95*5*2π + @test default_tv(tf(1, [1, 0, 1, 0]))[end] >= 0.95*5*2π + + # A user-provided final time + @test default_tv(tf(1, [1, 1]), 3)[end] == 3 + @test step(default_tv(tf(1, [1, 1]), 3)) == 3/200 + @test default_tv(c2d(tf(1, [1, 1]), 0.1), 3) == 0:0.1:3 +end + # Test error hints if VERSION >= v"1.7"