Use a stiffness-switching default solver in lsim and check the solver return code

[baggepinnen]

Part of #1081

Problem

The ODE-based lsim(sys, u::Function, t) for rational continuous-time systems used the explicit solver Tsit5() by default and did not check the return code of the solver. For a stiff system, the solver reaches maxiters and the method then fails with a BoundsError, because the solution contains fewer time points than requested:

using ControlSystems
sys = ss(tf(1, [1, 1e6+1, 1e6])) # poles in -1 and -1e6
lsim(sys, (x, t) -> [1.0], 0:0.1:10)
# Warning: Interrupted. Larger maxiters is needed. ...
# ERROR: BoundsError: attempt to access 2×36 Matrix{Float64} at index [1:2, 37]

The same BoundsError occurs for the diagonal realization ss(diagm([-1.0, -1e6]), [1.0; 1e6;;], [1.0 1.0], 0). The simulation methods for DelayLtiSystem and HammersteinWienerSystem did not check the return code either, and failed with a BoundsError when the solver terminated early, e.g., lsim(delay(1)*tf(1, [1, 1]), (x, t) -> [1.0], 0:0.1:10; maxiters=5) gave BoundsError: attempt to access 2-element Vector{Tuple{...}} at index [3].

Change

• The default solver of the rational method is now AutoTsit5(Rosenbrock23()). It integrates with Tsit5 and switches to Rosenbrock23 when stiffness is detected. Both are available from the existing dependency OrdinaryDiffEq (v6 and v7), so no dependency is added. The keyword argument alg is retained.
• A new internal function check_retcode(sol) throws an error when SciMLBase.successful_retcode(sol) is false. The error states the time at which the simulation terminated and the return code, and suggests selecting another solver with alg or adjusting abstol, reltol and maxiters. It is called after the solver in the rational, DelayLtiSystem and HammersteinWienerSystem methods. The solvers and tolerance defaults of the delay and Hammerstein-Wiener methods are unchanged.
• The lsim docstring (in ControlSystemsBase) states the default solver and the requirement below.

After the change, the reproducer above gives y[end] = 9.9995478e-7 (exact value 1e-6*(1 - e^{-10}) ≈ 9.99954600e-7) in 0.2 ms, and an unsuccessful simulation gives

The simulation terminated at t = 0.1 with return code MaxIters. Consider selecting another solver with the keyword argument `alg`, e.g., a solver for stiff problems, or adjusting the keyword arguments `abstol`, `reltol` and `maxiters`.

Effect on existing users

• Non-stiff problems: AutoTsit5 uses the same Tsit5 steps as long as no stiffness is detected. For ss(tf([1, 2], [1, 0.5, 2, 1])) with input sin(t) on 0:0.01:20, the result is identical to alg = Tsit5() (maximum difference 0.0); the run time is 4.4 ms compared to 3.2 ms with Tsit5().
• Moderately stiff problems, for which the stiffness detection switches to Rosenbrock23, give slightly different results within the solver tolerances. Note that with the default dtmax = t[2]-t[1], the switch already occurs for a fast pole in -100 and dt = 0.1. In the tested examples the results are closer to the exact solution than before: for a system of state dimension 6 with a step input and nonzero x0, the maximum error with respect to the exact zero-order-hold solution is 3.1e-4 (previously 4.9e-4); for the pole pair -1, -1e3 the error of y[end] is 1.4e-10 (previously 2.0e-9).
• Trade-off: Rosenbrock23 computes the Jacobian and the time derivative by forward-mode automatic differentiation, so after a switch u is evaluated with arguments of a ForwardDiff dual-number type. Input functions that do not support this, e.g., (x, t) -> (buf .= -L*x) with a preallocated Float64 buffer, now throw a MethodError for problems that are detected as stiff (for a fast pole in -100 and dt = 0.1, such a function worked with Tsit5). The docstring states this and recommends alg = Tsit5() in that case. A finite-difference Jacobian (Rosenbrock23(autodiff = AutoFiniteDiff())) would avoid this, but the autodiff keyword accepts an ADTypes specifier only in recent OrdinaryDiffEq versions (v7 rejects autodiff = false), so this would require ADTypes as a direct dependency and a higher lower compat bound of OrdinaryDiffEq; finite differences of the time derivative are in addition inaccurate at discontinuities of u such as steps.
• Unsuccessful simulations now throw an ErrorException that states the return code instead of a BoundsError.

Tests

• test/test_timeresp.jl: regression test with ss(tf(1, [1, 1e6+1, 1e6])), compared to the exact zero-order-hold simulation (rtol=1e-3), and a test that alg = Tsit5() with maxiters = 100 throws an error that contains "return code".
• test/test_delay_timeresp.jl, test/test_nonlinear_timeresp.jl: tests that the delay and Hammerstein-Wiener methods throw an error that contains "return code" when maxiters is reached (verified to give a BoundsError without the check).
• test/runtests_controlsystems.jl (test_timeresp.jl, test_delay_timeresp.jl, test_nonlinear_timeresp.jl) with OrdinaryDiffEq v7.8.1 and DelayDiffEq v6.4.2 on Julia 1.12: ControlSystems | 1156 passed, 1156 total.

---

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

🤖 Generated with Claude Code

https://claude.ai/code/session_01HigNnun4aJKMoWmY7wyULn

[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

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 91.75%. Comparing base (59c065a) to head (aa722a4).

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #1082      +/-   ##
==========================================
+ Coverage   91.73%   91.75%   +0.02%     
==========================================
  Files          42       42              
  Lines        5865     5871       +6     
==========================================
+ Hits         5380     5387       +7     
+ Misses        485      484       -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.