Skip to content

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

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
baggepinnen merged 1 commit into from
Oct 9, 2026
Merged

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

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
110 changes: 90 additions & 20 deletions lib/ControlSystemsBase/src/timeresp.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)`,
Expand Down Expand Up @@ -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)`,
Expand Down Expand Up @@ -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
Expand All @@ -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

"""
Expand Down
78 changes: 78 additions & 0 deletions lib/ControlSystemsBase/test/test_timeresp.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down
Loading