2 fixes. 1st fix: the m_nr branch which is meant to handle negative reals takes the -1 root and uses the p_nr case causes t_emp to become infinite due to ln|p|=0 causing a division by zero and a RuntimeWarning. FIX exclude |p|=1 so they use the m_w case which uses cycle counting formula which is appropriate for unit circle poles.

UbeenII · UbeenII · commit b59fba934e9f · 2026-05-16T00:12:01.000+01:00
2nd fix: The m_int case used p.real -1 <sqrt_eps instead of np.abs (p.real-1)<sqrt_eps. Without abs any pole with Re(p) < 1 were matched rather than just poles near +1 causing stable poles like 0.95 to be discarded as discrete integrators. I updated the corresponding test to have the expected value to match the correct result. Fixes #1204.
diff --git a/control/tests/timeresp_test.py b/control/tests/timeresp_test.py
@@ -818,14 +818,14 @@ def test_step_robustness(self):
     @pytest.mark.parametrize(
         "tfsys, tfinal",
         [(TransferFunction(1, [1, .5]), 13.81551),        #  pole at 0.5
-         (TransferFunction(1, [1, .5]).sample(.1), 25),  # discrete pole at 0.5
+         (TransferFunction(1, [1, .5]).sample(.1), 13.81551),  # discrete pole at 0.5
          (TransferFunction(1, [1, .5, 0]), 25)])         # poles at 0.5 and 0
     def test_auto_generated_time_vector_tfinal(self, tfsys, tfinal):
         """Confirm a TF with a pole at p simulates for tfinal seconds"""
         ideal_tfinal, ideal_dt = _ideal_tfinal_and_dt(tfsys)
         np.testing.assert_allclose(ideal_tfinal, tfinal, rtol=1e-4)
         T = _default_time_vector(tfsys)
-        np.testing.assert_allclose(T[-1], tfinal, atol=0.5*ideal_dt)
+        np.testing.assert_allclose(T[-1], tfinal, atol=ideal_dt)
 
     @pytest.mark.parametrize("wn, zeta", [(10, 0), (100, 0), (100, .1)])
     def test_auto_generated_time_vector_dt_cont1(self, wn, zeta):
diff --git a/control/timeresp.py b/control/timeresp.py
@@ -2172,13 +2172,13 @@ def _ideal_tfinal_and_dt(sys, is_step=True):
         m_z = np.abs(p) < sqrt_eps
         p = p[~m_z]
         # Negative reals- treated as oscillatory mode
-        m_nr = (p.real < 0) & (np.abs(p.imag) < sqrt_eps)
+        m_nr = (p.real < 0) & (np.abs(p.imag) < sqrt_eps)&(np.abs(p.real+1)>sqrt_eps)
         p_nr, p = p[m_nr], p[~m_nr]
         if p_nr.size > 0:
             t_emp = np.max(log_decay_percent / np.abs((np.log(p_nr)/dt).real))
             tfinal = max(tfinal, t_emp)
         # discrete integrators
-        m_int = (p.real - 1 < sqrt_eps) & (np.abs(p.imag) < sqrt_eps)
+        m_int = (np.abs(p.real - 1) < sqrt_eps) & (np.abs(p.imag) < sqrt_eps)
         p_int, p = p[m_int], p[~m_int]
         # pure oscillatory modes
         m_w = (np.abs(np.abs(p) - 1) < sqrt_eps)
