Repository: python-control/python-control

Issue:
phase_crossover_frequencies generates errors with poles at origin

### Discussed in https://github.com/python-control/python-control/discussions/1089

<div type='discussions-op-text'>

<sup>Originally posted by **alittlepinetree** December 31, 2024</sup>
When I used this function to calculate the phase crossover frequencies, it didn't work.
`import control
import matplotlib.pyplot as plt

G = control.tf([200.0], [1.0, 21.0, 20.0, 0.0])
x_omega, x_gain = control.phase_crossover_frequencies(G)
print(x_omega)`</div>

Discussion:
@murrayrm: Fixed in #1106 

Solution patch:
diff --git a/control/freqplot.py b/control/freqplot.py
index 3bc3a0c24..456431f38 100644
--- a/control/freqplot.py
+++ b/control/freqplot.py
@@ -1447,6 +1447,12 @@ def nyquist_response(
         else:
             contour = np.exp(splane_contour * sys.dt)
 
+        # Make sure we don't try to evaluate at a pole
+        if isinstance(sys, (StateSpace, TransferFunction)):
+            if any([pole in contour for pole in sys.poles()]):
+                raise RuntimeError(
+                    "attempt to evaluate at a pole; indent required")
+
         # Compute the primary curve
         resp = sys(contour)
 
diff --git a/control/margins.py b/control/margins.py
index 301baaf57..019c866be 100644
--- a/control/margins.py
+++ b/control/margins.py
@@ -81,11 +81,16 @@ def _poly_iw_sqr(pol_iw):
 
 def _poly_iw_real_crossing(num_iw, den_iw, epsw):
     # Return w where imag(H(iw)) == 0
+
+    # Compute the imaginary part of H = (num.r + j num.i)/(den.r + j den.i)
     test_w = np.polysub(np.polymul(num_iw.imag, den_iw.real),
                         np.polymul(num_iw.real, den_iw.imag))
+
+    # Find the real-valued w > 0 where imag(H(iw)) = 0
     w = np.roots(test_w)
     w = np.real(w[np.isreal(w)])
     w = w[w >= epsw]
+
     return w
 
 
@@ -471,7 +476,7 @@ def phase_crossover_frequencies(sys):
     omega : ndarray
         1d array of (non-negative) frequencies where Nyquist plot
         intersects the real axis
-    gain : ndarray
+    gains : ndarray
         1d array of corresponding gains
 
     Examples
@@ -493,13 +498,13 @@ def phase_crossover_frequencies(sys):
         omega = _poly_iw_real_crossing(num_iw, den_iw, 0.)
 
         # using real() to avoid rounding errors and results like 1+0j
-        gain = np.real(evalfr(sys, 1J * omega))
+        gains = np.real(sys(omega * 1j, warn_infinite=False))
     else:
         zargs = _poly_z_invz(sys)
         z, omega = _poly_z_real_crossing(*zargs, epsw=0.)
-        gain = np.real(evalfr(sys, z))
+        gains = np.real(sys(z, warn_infinite=False))
 
-    return omega, gain
+    return omega, gains
 
 
 def margin(*args):
diff --git a/control/statefbk.py b/control/statefbk.py
index 57df28d45..7b96c8015 100644
--- a/control/statefbk.py
+++ b/control/statefbk.py
@@ -1074,7 +1074,7 @@ def _control_output(t, states, inputs, params):
 
 
 def ctrb(A, B, t=None):
-    """Controllabilty matrix.
+    """Controllability matrix.
 
     Parameters
     ----------


Test patch:
diff --git a/control/tests/margin_test.py b/control/tests/margin_test.py
index 07e21114f..43cd68ae3 100644
--- a/control/tests/margin_test.py
+++ b/control/tests/margin_test.py
@@ -12,12 +12,9 @@
 from numpy import inf, nan
 from numpy.testing import assert_allclose
 
-from control.frdata import FrequencyResponseData
-from control.margins import (margin, phase_crossover_frequencies,
-                             stability_margins)
-from control.statesp import StateSpace
-from control.xferfcn import TransferFunction
-from control.exception import ControlMIMONotImplemented
+from control import ControlMIMONotImplemented, FrequencyResponseData, \
+    StateSpace, TransferFunction, margin, phase_crossover_frequencies, \
+    stability_margins
 
 s = TransferFunction.s
 
@@ -111,7 +108,6 @@ def test_margin_3input(tsys):
     out = margin((mag, phase*180/np.pi, omega_))
     assert_allclose(out, np.array(refout)[[0, 1, 3, 4]], atol=1.5e-3)
 
-
 @pytest.mark.parametrize(
     'tfargs, omega_ref, gain_ref',
     [(([1], [1, 2, 3, 4]), [1.7325, 0.], [-0.5, 0.25]),
@@ -119,7 +115,10 @@ def test_margin_3input(tsys):
      (([2], [1, 3, 3, 1]), [1.732, 0.], [-0.25, 2.]),
      ((np.array([3, 11, 3]) * 1e-4, [1., -2.7145, 2.4562, -0.7408], .1),
       [1.6235, 0.], [-0.28598, 1.88889]),
+     (([200.0], [1.0, 21.0, 20.0, 0.0]),
+      [4.47213595, 0], [-0.47619048, inf]),
      ])
+@pytest.mark.filterwarnings("error")
 def test_phase_crossover_frequencies(tfargs, omega_ref, gain_ref):
     """Test phase_crossover_frequencies() function"""
     sys = TransferFunction(*tfargs)
diff --git a/control/tests/nyquist_test.py b/control/tests/nyquist_test.py
index 0d6907b64..8d6fd8561 100644
--- a/control/tests/nyquist_test.py
+++ b/control/tests/nyquist_test.py
@@ -517,6 +517,15 @@ def test_nyquist_frd():
     warnings.resetwarnings()
 
 
+def test_no_indent_pole():
+    s = ct.tf('s')
+    sys = ((1 + 5/s)/(1 + 0.5/s))**2   # Double-Lag-Compensator
+
+    with pytest.raises(RuntimeError, match="evaluate at a pole"):
+        resp = ct.nyquist_response(
+            sys, warn_encirclements=False, indent_direction='none')
+
+
 if __name__ == "__main__":
     #
     # Interactive mode: generate plots for manual viewing