Repository: python-control/python-control

Issue:
Different impulse responses for zpk vs tf
### Discussed in https://github.com/python-control/python-control/discussions/1061

<div type='discussions-op-text'>

<sup>Originally posted by **MEnzmann** November 19, 2024</sup>
Hi all! 

I am little puzzled about the declaration of transfer-functions with control.zpk. The step-responses / impulse-responses of a transfer-function declared with control.zpk seem to unstable, although the poles are in the left-half plane and a bode plot does not show unstable poles. Am I missing something in the following MME?

Help would be appreciated! TIA!

_import control as ct
import control.matlab as ml
import matplotlib.pyplot as pl
#### Declare stable system as zpk 
Gzp = ct.zpk([-5],[-1,-10], gain=4)
print(Gzp)
#### Declare same system as tf
Gtf = ct.tf([4,20], [1,11, 10])
print(Gtf)
#### Both systems share the same poles and zeros, the difference of both systems is 0.
#### Yet: the Gzp system is instable in simulations while Gtf is stable (as expected)
[Tzp, Yzp] = ct.impulse_response(Gzp)
pl.plot(Tzp, Yzp)
pl.show()
[Ttf, Ytf] = ct.impulse_response(Gtf)
pl.plot(Ttf, Ytf)
pl.show()__</div>

Hints from issue comments:
@murrayrm: I took a quick look and confirmed that I am getting the same results (and that something is wrong).

The issue appears to be that `Gzp` above has timebase `None` while `Gtf` has timebase `0`.  In the `forced_response` function (which is called by `impulse_response`), the logic is set up so that a timebase of `None` gets interpreted as a discrete time system.  This means that the system is now unstable (since the magnitude of one of the poles is greater than 1).

I'll need to think through how to fix this in the right way.  As a very short term fix, if you set the timebase explicitly, the problem goes away:
```
Gzp = ct.zpk([-5],[-1,-10], gain=4, dt=0)
```



Solution patch:
diff --git a/control/frdata.py b/control/frdata.py
index c5018babb..42ecee0d9 100644
--- a/control/frdata.py
+++ b/control/frdata.py
@@ -124,13 +124,24 @@ def __init__(self, *args, **kwargs):
         To construct frequency response data for an existing LTI
         object, other than an FRD, call FRD(sys, omega).
 
+        The timebase for the frequency response can be provided using an
+        optional third argument or the 'dt' keyword.
+
         """
-        # TODO: discrete-time FRD systems?
         smooth = kwargs.pop('smooth', False)
 
         #
         # Process positional arguments
         #
+        if len(args) == 3:
+            # Discrete time transfer function
+            dt = args[-1]
+            if 'dt' in kwargs:
+                warn("received multiple dt arguments, "
+                     "using positional arg dt = %s" % dt)
+            kwargs['dt'] = dt
+            args = args[:-1]
+
         if len(args) == 2:
             if not isinstance(args[0], FRD) and isinstance(args[0], LTI):
                 # not an FRD, but still a system, second argument should be
@@ -200,11 +211,11 @@ def __init__(self, *args, **kwargs):
 
         # Process iosys keywords
         defaults = {
-            'inputs': self.fresp.shape[1], 'outputs': self.fresp.shape[0],
-            'dt': None}
+            'inputs': self.fresp.shape[1], 'outputs': self.fresp.shape[0]}
+        if arg_dt is not None:
+            defaults['dt'] = arg_dt             # choose compatible timebase
         name, inputs, outputs, states, dt = _process_iosys_keywords(
                 kwargs, defaults, end=True)
-        dt = common_timebase(dt, arg_dt)        # choose compatible timebase
 
         # Process signal names
         InputOutputSystem.__init__(
diff --git a/control/phaseplot.py b/control/phaseplot.py
index 6f7ed355f..e6123bc0e 100644
--- a/control/phaseplot.py
+++ b/control/phaseplot.py
@@ -532,8 +532,7 @@ def equilpoints(
     out = []
     for xeq in equilpts:
         with plt.rc_context(rcParams):
-            out.append(
-                ax.plot(xeq[0], xeq[1], marker='o', color=color))
+            out += ax.plot(xeq[0], xeq[1], marker='o', color=color)
     return out
 
 
@@ -649,8 +648,7 @@ def separatrices(
     for i, xeq in enumerate(equilpts):
         # Plot the equilibrium points
         with plt.rc_context(rcParams):
-            out.append(
-                ax.plot(xeq[0], xeq[1], marker='o', color='k'))
+            out += ax.plot(xeq[0], xeq[1], marker='o', color='k')
 
         # Figure out the linearization and eigenvectors
         evals, evecs = np.linalg.eig(sys.linearize(xeq, 0, params=params).A)
@@ -691,13 +689,13 @@ def separatrices(
                 # Plot the trajectory (if there is one)
                 if traj.shape[1] > 1:
                     with plt.rc_context(rcParams):
-                        out.append(ax.plot(
-                            traj[0], traj[1], color=color, linestyle=linestyle))
+                        out += ax.plot(
+                            traj[0], traj[1], color=color, linestyle=linestyle)
 
                     # Add arrows to the lines at specified intervals
                     with plt.rc_context(rcParams):
                         _add_arrows_to_line2D(
-                            ax, out[-1][0], arrow_pos, arrowstyle=arrow_style,
+                            ax, out[-1], arrow_pos, arrowstyle=arrow_style,
                             dir=1)
     return out
 
diff --git a/control/pzmap.py b/control/pzmap.py
index 7b0f8d096..400db2807 100644
--- a/control/pzmap.py
+++ b/control/pzmap.py
@@ -607,14 +607,16 @@ def _compute_root_locus_limits(response):
         # Find the local maxima of root locus curve
         xpeaks = np.where(
             np.diff(np.abs(locus.real)) < 0, locus.real[0:-1], 0)
-        xlim = [
-            min(xlim[0], np.min(xpeaks) * rho),
-            max(xlim[1], np.max(xpeaks) * rho)
-        ]
+        if xpeaks.size > 0:
+            xlim = [
+                min(xlim[0], np.min(xpeaks) * rho),
+                max(xlim[1], np.max(xpeaks) * rho)
+            ]
 
         ypeaks = np.where(
             np.diff(np.abs(locus.imag)) < 0, locus.imag[0:-1], 0)
-        ylim = max(ylim, np.max(ypeaks) * rho)
+        if ypeaks.size > 0:
+            ylim = max(ylim, np.max(ypeaks) * rho)
 
     if isctime(dt=response.dt):
         # Adjust the limits to include some space around features
diff --git a/control/rlocus.py b/control/rlocus.py
index 189c2ccd0..54c672997 100644
--- a/control/rlocus.py
+++ b/control/rlocus.py
@@ -56,10 +56,9 @@ def root_locus_map(sysdata, gains=None):
     Returns
     -------
     rldata : PoleZeroData or list of PoleZeroData
-        Root locus data object(s) corresponding to the .  The loci of
-        the root locus diagram are available in the array
-        `rldata.loci`, indexed by the gain index and the locus index,
-        and the gains are in the array `rldata.gains`.
+        Root locus data object(s).  The loci of the root locus diagram are
+        available in the array `rldata.loci`, indexed by the gain index and
+        the locus index, and the gains are in the array `rldata.gains`.
 
     Notes
     -----
diff --git a/control/xferfcn.py b/control/xferfcn.py
index ee41cbd2b..499359cbc 100644
--- a/control/xferfcn.py
+++ b/control/xferfcn.py
@@ -1677,7 +1677,7 @@ def tf(*args, **kwargs):
         raise ValueError("Needs 1 or 2 arguments; received %i." % len(args))
 
 
-def zpk(zeros, poles, gain, dt=None, **kwargs):
+def zpk(zeros, poles, gain, *args, **kwargs):
     """zpk(zeros, poles, gain[, dt])
 
     Create a transfer function from zeros, poles, gain.
@@ -1732,7 +1732,7 @@ def zpk(zeros, poles, gain, dt=None, **kwargs):
 
     """
     num, den = zpk2tf(zeros, poles, gain)
-    return TransferFunction(num, den, dt=dt, **kwargs)
+    return TransferFunction(num, den, *args, **kwargs)
 
 
 def ss2tf(*args, **kwargs):


Test patch:
diff --git a/control/tests/docstrings_test.py b/control/tests/docstrings_test.py
index 2b1368aeb..27ced105f 100644
--- a/control/tests/docstrings_test.py
+++ b/control/tests/docstrings_test.py
@@ -54,6 +54,7 @@
     control.bode_plot: ['sharex', 'sharey', 'margin_info'], # deprecated
     control.eigensys_realization: ['arg'],                  # quasi-positional
     control.find_operating_point: ['method'],               # internal use
+    control.zpk: ['args']                                   # 'dt' (manual)
 }
 
 # Decide on the level of verbosity (use -rP when running pytest)
diff --git a/control/tests/rlocus_test.py b/control/tests/rlocus_test.py
index a62bc742b..52546a428 100644
--- a/control/tests/rlocus_test.py
+++ b/control/tests/rlocus_test.py
@@ -287,3 +287,18 @@ def test_root_locus_documentation(savefigs=False):
 
     # Run tests that generate plots for the documentation
     test_root_locus_documentation(savefigs=True)
+
+
+# https://github.com/python-control/python-control/issues/1063
+def test_rlocus_singleton():
+    # Generate a root locus map for a singleton
+    L = ct.tf([1, 1], [1, 2, 3])
+    rldata = ct.root_locus_map(L, 1)
+    np.testing.assert_equal(rldata.gains, np.array([1]))
+    assert rldata.loci.shape == (1, 2)
+
+    # Generate the root locus plot (no loci)
+    cplt = rldata.plot()
+    assert len(cplt.lines[0, 0]) == 1      # poles (one set of markers)
+    assert len(cplt.lines[0, 1]) == 1      # zeros
+    assert len(cplt.lines[0, 2]) == 2      # loci (two 0-length lines)
diff --git a/control/tests/timebase_test.py b/control/tests/timebase_test.py
index 79b1492d7..c0e02d3b8 100644
--- a/control/tests/timebase_test.py
+++ b/control/tests/timebase_test.py
@@ -97,3 +97,34 @@ def test_composition_override(dt):
         with pytest.raises(ValueError, match="incompatible timebases"):
             sys3 = ct.interconnect(
                 [sys1, sys2], inputs='u1', outputs='y2', dt=dt)
+
+
+# Make sure all system creation functions treat timebases uniformly
+@pytest.mark.parametrize(
+    "fcn, args", [
+        (ct.ss, [-1, 1, 1, 1]),
+        (ct.tf, [[1, 2], [3, 4, 5]]),
+        (ct.zpk, [[-1], [-2, -3], 1]),
+        (ct.frd, [[1, 1, 1], [1, 2, 3]]),
+        (ct.nlsys, [lambda t, x, u, params: -x, None]),
+    ])
+@pytest.mark.parametrize(
+    "kwargs, expected", [
+        ({}, 0),
+        ({'dt': 0}, 0),
+        ({'dt': 0.1}, 0.1),
+        ({'dt': True}, True),
+        ({'dt': None}, None),
+    ])
+def test_default(fcn, args, kwargs, expected):
+    sys = fcn(*args, **kwargs)
+    assert sys.dt == expected
+
+    # Some commands allow dt via extra argument
+    if fcn in [ct.ss, ct.tf, ct.zpk, ct.frd] and kwargs.get('dt'):
+        sys = fcn(*args, kwargs['dt'])
+        assert sys.dt == expected
+
+        # Make sure an error is generated if dt is redundant
+        with pytest.warns(UserWarning, match="received multiple dt"):
+            sys = fcn(*args, kwargs['dt'], **kwargs)