allowed Bode/Nyquist for mixed FRD systems

murrayrm · murrayrm · commit 4b3597fa2157 · 2024-06-16T07:19:44.000-07:00
diff --git a/control/freqplot.py b/control/freqplot.py
@@ -284,7 +284,7 @@ def bode_plot(
     #
 
     # If we were passed a list of systems, convert to data
-    if all([isinstance(
+    if any([isinstance(
             sys, (StateSpace, TransferFunction)) for sys in data]):
         data = frequency_response(
             data, omega=omega, omega_limits=omega_limits,
@@ -1276,7 +1276,11 @@ def nyquist_response(
                 "Nyquist plot currently only supports SISO systems.")
 
         # Figure out the frequency range
-        omega_sys = np.asarray(omega)
+        if isinstance(sys, FrequencyResponseData) and sys.ifunc is None \
+           and not omega_range_given:
+            omega_sys = sys.omega               # use system frequencies
+        else:
+            omega_sys = np.asarray(omega)       # use common omega vector
 
         # Determine the contour used to evaluate the Nyquist curve
         if sys.isdtime(strict=True):
@@ -2477,18 +2481,6 @@ def _determine_omega_vector(syslist, omega_in, omega_limits, omega_num,
         and omega_limits are None.
 
     """
-    # Special processing for FRD systems
-    # TODO: allow different ranges of frequencies
-    if omega_in is None:
-        for sys in syslist:
-            if isinstance(sys, FrequencyResponseData):
-                # FRD already has predetermined frequencies
-                if omega_in is not None and not np.all(omega_in == sys.omega):
-                    raise ValueError(
-                        "List of FrequencyResponseData systems can only have "
-                        "a single frequency range between them")
-                omega_in = sys.omega
-
     # Handle the special case of a range of frequencies
     if omega_in is not None and omega_limits is not None:
         warnings.warn(
@@ -2573,6 +2565,15 @@ def _default_frequency_range(syslist, Hz=None, number_of_samples=None,
         syslist = (syslist,)
 
     for sys in syslist:
+        # For FRD systems, just use the response frequencies
+        if isinstance(sys, FrequencyResponseData):
+            # Add the min and max frequency, minus periphery decades
+            # (keeps frequency ranges from artificially expanding)
+            features = np.concatenate([features, np.array([
+                np.min(sys.omega) * 10**feature_periphery_decades,
+                np.max(sys.omega) / 10**feature_periphery_decades])])
+            continue
+
         try:
             # Add new features to the list
             if sys.isctime():
@@ -2587,7 +2588,8 @@ def _default_frequency_range(syslist, Hz=None, number_of_samples=None,
                 # TODO: What distance to the Nyquist frequency is appropriate?
                 freq_interesting.append(fn * 0.9)
 
-                features_ = np.concatenate((sys.poles(), sys.zeros()))
+                features_ = np.concatenate(
+                    (np.abs(sys.poles()), np.abs(sys.zeros())))
                 # Get rid of poles and zeros on the real axis (imag==0)
                 # * origin and real < 0
                 # * at 1.: would result in omega=0. (logaritmic plot!)
@@ -2602,8 +2604,9 @@ def _default_frequency_range(syslist, Hz=None, number_of_samples=None,
                 # TODO
                 raise NotImplementedError(
                     "type of system in not implemented now")
-            features = np.concatenate((features, features_))
+            features = np.concatenate([features, features_])
         except NotImplementedError:
+            # Don't add any features for anything we don't understand
             pass
 
     # Make sure there is at least one point in the range
diff --git a/control/lti.py b/control/lti.py
@@ -475,6 +475,7 @@ def frequency_response(
         #>>> # s = 0.1i, i, 10i.
 
     """
+    from .frdata import FrequencyResponseData
     from .freqplot import _determine_omega_vector
 
     # Process keyword arguments
@@ -489,13 +490,18 @@ def frequency_response(
 
     responses = []
     for sys_ in syslist:
-        # Add the Nyquist frequency for discrete time systems
-        omega_sys = omega_syslist.copy()
-        if sys_.isdtime(strict=True):
-            nyquistfrq = math.pi / sys_.dt
-            if not omega_range_given:
-                # Limit up to the Nyquist frequency
-                omega_sys = omega_sys[omega_sys < nyquistfrq]
+        if isinstance(sys_, FrequencyResponseData) and sys_.ifunc is None and \
+           not omega_range_given:
+            omega_sys = sys_.omega              # use system properties
+        else:
+            omega_sys = omega_syslist.copy()    # use common omega vector
+
+            # Add the Nyquist frequency for discrete time systems
+            if sys_.isdtime(strict=True):
+                nyquistfrq = math.pi / sys_.dt
+                if not omega_range_given:
+                    # Limit up to the Nyquist frequency
+                    omega_sys = omega_sys[omega_sys < nyquistfrq]
 
         # Compute the frequency response
         responses.append(sys_.frequency_response(omega_sys, squeeze=squeeze))
diff --git a/control/tests/freqplot_test.py b/control/tests/freqplot_test.py
@@ -523,6 +523,34 @@ def test_freqplot_ax_keyword(plt_fcn, ninputs, noutputs):
     np.testing.assert_equal(ct.get_plot_axes(out1), ct.get_plot_axes(out3))
 
 
+def test_mixed_systypes():
+    s = ct.tf('s')
+    sys_tf = ct.tf(
+        (0.02 * s**3 - 0.1 * s) / (s**4 + s**3 + s**2 + 0.25 * s + 0.04),
+        name='tf')
+    sys_ss = ct.ss(sys_tf * 2, name='ss')
+    sys_frd1 = ct.frd(sys_tf / 2, np.logspace(-1, 1, 15), name='frd1')
+    sys_frd2 = ct.frd(sys_tf / 4, np.logspace(-3, 2, 20), name='frd2')
+
+    # Simple case: compute responses separately and plot
+    resp_tf = ct.frequency_response(sys_tf)
+    resp_ss = ct.frequency_response(sys_ss)
+    plt.figure()
+    ct.bode_plot([resp_tf, resp_ss, sys_frd1, sys_frd2], plot_phase=False)
+    ct.suptitle("bode_plot([resp_tf, resp_ss, sys_frd1, sys_frd2])")
+
+    # Same thing, but using frequency response
+    plt.figure()
+    resp = ct.frequency_response([sys_tf, sys_ss, sys_frd1, sys_frd2])
+    resp.plot(plot_phase=False)
+    ct.suptitle("frequency_response([sys_tf, sys_ss, sys_frd1, sys_frd2])")
+
+    # Same thing, but using bode_plot
+    plt.figure()
+    resp = ct.bode_plot([sys_tf, sys_ss, sys_frd1, sys_frd2], plot_phase=False)
+    ct.suptitle("bode_plot([sys_tf, sys_ss, sys_frd1, sys_frd2])")
+
+
 def test_suptitle():
     sys = ct.rss(2, 2, 2)
 
@@ -623,5 +651,4 @@ def test_freqplot_errors(plt_fcn):
     # Run a few more special cases to show off capabilities (and save some
     # of them for use in the documentation).
     #
-
-    pass
+    test_mixed_systypes()
diff --git a/control/tests/nyquist_test.py b/control/tests/nyquist_test.py
@@ -162,35 +162,35 @@ def test_nyquist_fbs_examples():
 
     """Run through various examples from FBS2e to compare plots"""
     plt.figure()
-    plt.title("Figure 10.4: L(s) = 1.4 e^{-s}/(s+1)^2")
+    ct.suptitle("Figure 10.4: L(s) = 1.4 e^{-s}/(s+1)^2")
     sys = ct.tf([1.4], [1, 2, 1]) * ct.tf(*ct.pade(1, 4))
     response = ct.nyquist_response(sys)
     response.plot()
     assert _Z(sys) == response.count + _P(sys)
 
     plt.figure()
-    plt.title("Figure 10.4: L(s) = 1/(s + a)^2 with a = 0.6")
+    ct.suptitle("Figure 10.4: L(s) = 1/(s + a)^2 with a = 0.6")
     sys = 1/(s + 0.6)**3
     response = ct.nyquist_response(sys)
     response.plot()
     assert _Z(sys) == response.count + _P(sys)
 
     plt.figure()
-    plt.title("Figure 10.6: L(s) = 1/(s (s+1)^2) - pole at the origin")
+    ct.suptitle("Figure 10.6: L(s) = 1/(s (s+1)^2) - pole at the origin")
     sys = 1/(s * (s+1)**2)
     response = ct.nyquist_response(sys)
     response.plot()
     assert _Z(sys) == response.count + _P(sys)
 
     plt.figure()
-    plt.title("Figure 10.10: L(s) = 3 (s+6)^2 / (s (s+1)^2)")
+    ct.suptitle("Figure 10.10: L(s) = 3 (s+6)^2 / (s (s+1)^2)")
     sys = 3 * (s+6)**2 / (s * (s+1)**2)
     response = ct.nyquist_response(sys)
     response.plot()
     assert _Z(sys) == response.count + _P(sys)
 
     plt.figure()
-    plt.title("Figure 10.10: L(s) = 3 (s+6)^2 / (s (s+1)^2) [zoom]")
+    ct.suptitle("Figure 10.10: L(s) = 3 (s+6)^2 / (s (s+1)^2) [zoom]")
     with pytest.warns(UserWarning, match="encirclements does not match"):
         response = ct.nyquist_response(sys, omega_limits=[1.5, 1e3])
         response.plot()
@@ -208,7 +208,7 @@ def test_nyquist_fbs_examples():
 def test_nyquist_arrows(arrows):
     sys = ct.tf([1.4], [1, 2, 1]) * ct.tf(*ct.pade(1, 4))
     plt.figure();
-    plt.title("L(s) = 1.4 e^{-s}/(s+1)^2 / arrows = %s" % arrows)
+    ct.suptitle("L(s) = 1.4 e^{-s}/(s+1)^2 / arrows = %s" % arrows)
     response = ct.nyquist_response(sys)
     response.plot(arrows=arrows)
     assert _Z(sys) == response.count + _P(sys)
@@ -222,13 +222,13 @@ def test_nyquist_encirclements():
     plt.figure();
     response = ct.nyquist_response(sys)
     response.plot()
-    plt.title("Stable system; encirclements = %d" % response.count)
+    ct.suptitle("Stable system; encirclements = %d" % response.count)
     assert _Z(sys) == response.count + _P(sys)
 
     plt.figure();
     response = ct.nyquist_response(sys * 3)
     response.plot()
-    plt.title("Unstable system; encirclements = %d" %response.count)
+    ct.suptitle("Unstable system; encirclements = %d" %response.count)
     assert _Z(sys * 3) == response.count + _P(sys * 3)
 
     # System with pole at the origin
@@ -237,7 +237,7 @@ def test_nyquist_encirclements():
     plt.figure();
     response = ct.nyquist_response(sys)
     response.plot()
-    plt.title("Pole at the origin; encirclements = %d" %response.count)
+    ct.suptitle("Pole at the origin; encirclements = %d" %response.count)
     assert _Z(sys) == response.count + _P(sys)
 
     # Non-integer number of encirclements
@@ -251,7 +251,7 @@ def test_nyquist_encirclements():
         response = ct.nyquist_response(
             sys, omega_limits=[0.5, 1e3], encirclement_threshold=0.2)
         response.plot()
-    plt.title("Non-integer number of encirclements [%g]" %response.count)
+    ct.suptitle("Non-integer number of encirclements [%g]" %response.count)
 
 
 @pytest.fixture
@@ -266,7 +266,7 @@ def test_nyquist_indent_default(indentsys):
     plt.figure();
     response = ct.nyquist_response(indentsys)
     response.plot()
-    plt.title("Pole at origin; indent_radius=default")
+    ct.suptitle("Pole at origin; indent_radius=default")
     assert _Z(indentsys) == response.count + _P(indentsys)
 
 
@@ -293,7 +293,7 @@ def test_nyquist_indent_do(indentsys):
         indentsys, indent_radius=0.01, return_contour=True)
     count, contour = response
     response.plot()
-    plt.title("Pole at origin; indent_radius=0.01; encirclements = %d" % count)
+    ct.suptitle("Pole at origin; indent_radius=0.01; encirclements = %d" % count)
     assert _Z(indentsys) == count + _P(indentsys)
     # indent radius is smaller than the start of the default omega vector
     # check that a quarter circle around the pole at origin has been added.
@@ -314,7 +314,7 @@ def test_nyquist_indent_left(indentsys):
     plt.figure();
     response = ct.nyquist_response(indentsys, indent_direction='left')
     response.plot()
-    plt.title(
+    ct.suptitle(
         "Pole at origin; indent_direction='left'; encirclements = %d" %
         response.count)
     assert _Z(indentsys) == response.count + _P(indentsys, indent='left')
@@ -328,14 +328,14 @@ def test_nyquist_indent_im():
     plt.figure();
     response = ct.nyquist_response(sys)
     response.plot()
-    plt.title("Imaginary poles; encirclements = %d" % response.count)
+    ct.suptitle("Imaginary poles; encirclements = %d" % response.count)
     assert _Z(sys) == response.count + _P(sys)
 
     # Imaginary poles with indentation to the left
     plt.figure();
     response = ct.nyquist_response(sys, indent_direction='left')
     response.plot(label_freq=300)
-    plt.title(
+    ct.suptitle(
         "Imaginary poles; indent_direction='left'; encirclements = %d" %
         response.count)
     assert _Z(sys) == response.count + _P(sys, indent='left')
@@ -346,7 +346,7 @@ def test_nyquist_indent_im():
         response = ct.nyquist_response(
             sys, np.linspace(0, 1e3, 1000), indent_direction='none')
         response.plot()
-    plt.title(
+    ct.suptitle(
         "Imaginary poles; indent_direction='none'; encirclements = %d" %
         response.count)
     assert _Z(sys) == response.count + _P(sys)
@@ -465,6 +465,36 @@ def test_freqresp_omega_limits():
         np.array([resp0.contour[1], resp0.contour[-1]]))
 
 
+def test_nyquist_frd():
+    sys = ct.rss(4, 1, 1)
+    sys1 = ct.frd(sys, np.logspace(-1, 1, 10), name='sys1')
+    sys2 = ct.frd(sys, np.logspace(-2, 2, 10), name='sys2')
+    sys3 = ct.frd(sys, np.logspace(-2, 2, 10), smooth=True, name='sys3')
+
+    # Turn off warnings about number of encirclements
+    warnings.filterwarnings(
+        'ignore', message="number of encirclements was a non-integer value",
+        category=UserWarning)
+
+    # OK to specify frequency with FRD sys if frequencies match
+    nyqresp = ct.nyquist_response(sys1, np.logspace(-1, 1, 10))
+    np.testing.assert_allclose(nyqresp.contour, np.logspace(-1, 1, 10) * 1j)
+
+    # If a fixed FRD omega is used, generate an error on mismatch
+    with pytest.raises(ValueError, match="not all frequencies .* in .* list"):
+        nyqresp = ct.nyquist_response(sys2, np.logspace(-1, 1, 10))
+
+    # OK to specify frequency with FRD sys if interpolating FRD is used
+    nyqresp = ct.nyquist_response(sys3, np.logspace(-1, 1, 12))
+    np.testing.assert_allclose(nyqresp.contour, np.logspace(-1, 1, 12) * 1j)
+
+    # Computing Nyquist response w/ different frequencies OK if given as a list
+    nyqresp = ct.nyquist_response([sys1, sys2])
+    out = nyqresp.plot()
+
+    warnings.resetwarnings()
+
+
 if __name__ == "__main__":
     #
     # Interactive mode: generate plots for manual viewing
@@ -508,7 +538,7 @@ def test_freqresp_omega_limits():
     print("Unusual Nyquist plot")
     sys = ct.tf([1], [1, 3, 2]) * ct.tf([1], [1, 0, 1])
     plt.figure()
-    plt.title("Poles: %s" %
+    ct.suptitle("Poles: %s" %
               np.array2string(sys.poles(), precision=2, separator=','))
     response = ct.nyquist_response(sys)
     response.plot()
@@ -517,10 +547,17 @@ def test_freqresp_omega_limits():
     print("Discrete time systems")
     sys = ct.c2d(sys, 0.01)
     plt.figure()
-    plt.title("Discrete-time; poles: %s" %
+    ct.suptitle("Discrete-time; poles: %s" %
               np.array2string(sys.poles(), precision=2, separator=','))
     response = ct.nyquist_response(sys)
     response.plot()
 
-
-
+    print("Frequency response data (FRD) systems")
+    sys = ct.tf(
+        (0.02 * s**3 - 0.1 * s) / (s**4 + s**3 + s**2 + 0.25 * s + 0.04),
+        name='tf')
+    sys1 = ct.frd(sys, np.logspace(-1, 1, 15), name='frd1')
+    sys2 = ct.frd(sys, np.logspace(-2, 2, 20), name='frd2')
+    plt.figure()
+    ct.nyquist_plot([sys, sys1, sys2])
+    ct.suptitle("Mixed FRD, tf data")
