python-control · murrayrm · Feb 5, 2018 · Jan 13, 2018 · Jan 15, 2018 · Jan 15, 2018
diff --git a/control/freqplot.py b/control/freqplot.py
@@ -171,22 +171,47 @@ def bode_plot(syslist, omega=None, dB=None, Hz=None, deg=None,
             #! TODO: Not current implemented; just use subplot for now
 
             if (Plot):
+                # Set up the axes with labels so that multiple calls to
+                # bode_plot will superimpose the data.  This was implicit
+                # before matplotlib 2.1, but changed after that (See
+                # https://github.com/matplotlib/matplotlib/issues/9024). 
+                # The code below should work on all cases.
+
+                # Get the current figure 
+                fig = plt.gcf()
+                ax_mag = None
+                ax_phase = None
+
+                # Get the current axes if they already exist
+                for ax in fig.axes:
+                    if ax.get_label() == 'control-bode-magnitude':
+                        ax_mag = ax
+                    elif ax.get_label() == 'control-bode-phase':
+                        ax_phase = ax
+
+                # If no axes present, create them from scratch
+                if ax_mag is None or ax_phase is None:
+                    plt.clf()
+                    ax_mag = plt.subplot(211, label = 'control-bode-magnitude')
+                    ax_phase = plt.subplot(212, label = 'control-bode-phase',
+                                           sharex=ax_mag)
+
                 # Magnitude plot
-                ax_mag = plt.subplot(211);
                 if dB:
-                    pltline = ax_mag.semilogx(omega_plot, 20 * np.log10(mag), *args, **kwargs)
+                    pltline = ax_mag.semilogx(omega_plot, 20 * np.log10(mag),
+                                              *args, **kwargs)
                 else:
                     pltline = ax_mag.loglog(omega_plot, mag, *args, **kwargs)
 
                 if nyquistfrq_plot:
-                    ax_mag.axvline(nyquistfrq_plot, color=pltline[0].get_color())
+                    ax_mag.axvline(nyquistfrq_plot,
+                                   color=pltline[0].get_color())
 
                 # Add a grid to the plot + labeling
                 ax_mag.grid(True, which='both')
                 ax_mag.set_ylabel("Magnitude (dB)" if dB else "Magnitude")
 
                 # Phase plot
-                ax_phase = plt.subplot(212, sharex=ax_mag);
                 if deg:
                     phase_plot = phase * 180. / math.pi
                 else:
@@ -353,28 +378,50 @@ def gangof4_plot(P, C, omega=None):
         L = P * C;
         S = feedback(1, L);
         T = L * S;
-
+
+        # Set up the axes with labels so that multiple calls to
+        # gangof4_plot will superimpose the data.  See details in bode_plot.
+        plot_axes = {'t' : None, 's' : None, 'ps' : None, 'cs' : None}
+        for ax in plt.gcf().axes:
+            label = ax.get_label()
+            if label.startswith('control-gangof4-'):
+                key = label[len('control-gangof4-'):]
+                if key not in plot_axes:
+                    raise RuntimeError("unknown gangof4 axis type '{}'".format(label))
+                plot_axes[key] = ax
+
+        # if any of the axes are missing, start from scratch
+        if any((ax is None for ax in plot_axes.values())):
+            plt.clf()
+            plot_axes = {'t' : plt.subplot(221,label='control-gangof4-t'),
+                         'ps' : plt.subplot(222,label='control-gangof4-ps'),
+                         'cs' : plt.subplot(223,label='control-gangof4-cs'),
+                         's' : plt.subplot(224,label='control-gangof4-s')}
+
+        #
         # Plot the four sensitivity functions
+        #
+
         #! TODO: Need to add in the mag = 1 lines
         mag_tmp, phase_tmp, omega = T.freqresp(omega);
         mag = np.squeeze(mag_tmp)
         phase = np.squeeze(phase_tmp)
-        plt.subplot(221); plt.loglog(omega, mag);
+        plot_axes['t'].loglog(omega, mag);
 
         mag_tmp, phase_tmp, omega = (P * S).freqresp(omega);
         mag = np.squeeze(mag_tmp)
         phase = np.squeeze(phase_tmp)
-        plt.subplot(222); plt.loglog(omega, mag);
+        plot_axes['ps'].loglog(omega, mag);
 
         mag_tmp, phase_tmp, omega = (C * S).freqresp(omega);
         mag = np.squeeze(mag_tmp)
         phase = np.squeeze(phase_tmp)
-        plt.subplot(223); plt.loglog(omega, mag);
+        plot_axes['cs'].loglog(omega, mag);
 
         mag_tmp, phase_tmp, omega = S.freqresp(omega);
         mag = np.squeeze(mag_tmp)
         phase = np.squeeze(phase_tmp)
-        plt.subplot(224); plt.loglog(omega, mag);
+        plot_axes['s'].loglog(omega, mag);
 
 #
 # Utility functions

diff --git a/control/statesp.py b/control/statesp.py
@@ -1078,18 +1078,18 @@ def ss(*args):
         output equations:
 
         .. math::
-            \dot x = A \cdot x + B \cdot u
+            \\dot x = A \\cdot x + B \\cdot u
 
-            y = C \cdot x + D \cdot u
+            y = C \\cdot x + D \\cdot u
 
     ``ss(A, B, C, D, dt)``
         Create a discrete-time state space system from the matrices of
         its state and output equations:
 
         .. math::
-            x[k+1] = A \cdot x[k] + B \cdot u[k]
+            x[k+1] = A \\cdot x[k] + B \\cdot u[k]
 
-            y[k] = C \cdot x[k] + D \cdot u[ki]
+            y[k] = C \\cdot x[k] + D \\cdot u[ki]
 
         The matrices can be given as *array like* data types or strings.
         Everything that the constructor of :class:`numpy.matrix` accepts is

diff --git a/control/tests/freqresp_test.py b/control/tests/freqresp_test.py
@@ -8,6 +8,7 @@
 
 import unittest
 import numpy as np
+import control as ctrl
 from control.statesp import StateSpace
 from control.matlab import ss, tf, bode
 from control.exception import slycot_check
@@ -34,6 +35,50 @@ def test_siso(self):
       systf = tf(sys)
       bode(systf)
 
+   def test_superimpose(self):
+      # Test to make sure that multiple calls to plots superimpose their
+      # data on the same axes unless told to do otherwise
+
+      # Generate two plots in a row; should be on the same axes
+      plt.figure(1); plt.clf()
+      ctrl.bode_plot(ctrl.tf([1], [1,2,1]))
+      ctrl.bode_plot(ctrl.tf([5], [1, 1]))
+
+      # Check to make sure there are two axes and that each axes has two lines
+      assert len(plt.gcf().axes) == 2
+      for ax in plt.gcf().axes:
+         # Make sure there are 2 lines in each subplot
+         assert len(ax.get_lines()) == 2
+
+      # Generate two plots as a list; should be on the same axes
+      plt.figure(2); plt.clf();
+      ctrl.bode_plot([ctrl.tf([1], [1,2,1]), ctrl.tf([5], [1, 1])])
+
+      # Check to make sure there are two axes and that each axes has two lines
+      assert len(plt.gcf().axes) == 2
+      for ax in plt.gcf().axes:
+         # Make sure there are 2 lines in each subplot
+         assert len(ax.get_lines()) == 2
+
+      # Generate two separate plots; only the second should appear
+      plt.figure(3); plt.clf();
+      ctrl.bode_plot(ctrl.tf([1], [1,2,1]))
+      plt.clf()
+      ctrl.bode_plot(ctrl.tf([5], [1, 1]))
+
+      # Check to make sure there are two axes and that each axes has one line
+      assert len(plt.gcf().axes) == 2
+      for ax in plt.gcf().axes:
+         # Make sure there is only 1 line in the subplot
+         assert len(ax.get_lines()) == 1
+
+      # Now add a line to the magnitude plot and make sure if is there
+      for ax in plt.gcf().axes:
+         if ax.get_label() == 'control-bode-magnitude':
+            break
+      ax.semilogx([1e-2, 1e1], 20 * np.log10([1, 1]), 'k-')
+      assert len(ax.get_lines()) == 2
+
    def test_doubleint(self):
       # 30 May 2016, RMM: added to replicate typecast bug in freqresp.py
       A = np.matrix('0, 1; 0, 0');

diff --git a/doc/control.rst b/doc/control.rst
@@ -44,6 +44,15 @@ Frequency domain plotting
     gangof4_plot
     nichols_plot
 
+Note: For plotting commands that create multiple axes on the same plot, the
+individual axes can be retrieved using the axes label (retrieved using the
+`get_label` method for the matplotliib axes object).  The following labels
+are currently defined:
+
+* Bode plots: `control-bode-magnitude`, `control-bode-phase`
+* Gang of 4 plots: `control-gangof4-s`, `control-gangof4-cs`,
+  `control-gangof4-ps`, `control-gangof4-t`
+
 Time domain simulation
 ======================
 

diff --git a/examples/pvtol-nested.py b/examples/pvtol-nested.py
@@ -81,26 +81,45 @@
 T = feedback(L, 1);
 
 # Compute stability margins
-#! Not yet implemented
-# (gm, pm, wgc, wpc) = margin(L); 
+(gm, pm, wgc, wpc) = margin(L);
+print("Gain margin: %g at %g" % (gm, wgc))
+print("Phase margin: %g at %g" % (pm, wpc))
 
-#! TODO: this figure has something wrong; axis limits mismatch
 figure(6); clf; 
-bode(L);
+bode(L, logspace(-4, 3));
 
-# Add crossover line
-subplot(211); hold(True);
-loglog([1e-4, 1e3], [1, 1], 'k-')
+# Add crossover line to the magnitude plot
+#
+# Note: in matplotlib before v2.1, the following code worked:
+#
+#   subplot(211); hold(True);
+#   loglog([1e-4, 1e3], [1, 1], 'k-')
+#
+# In later versions of matplotlib the call to subplot will clear the
+# axes and so we have to extract the axes that we want to use by hand.
+# In addition, hold() is deprecated so we no longer require it.
+#
+for ax in gcf().axes:
+    if ax.get_label() == 'control-bode-magnitude':
+        break
+ax.semilogx([1e-4, 1e3], 20 * np.log10([1, 1]), 'k-')
 
+#
 # Replot phase starting at -90 degrees
-bode(L, logspace(-4, 3));
+#
+# Get the phase plot axes
+for ax in gcf().axes:
+    if ax.get_label() == 'control-bode-phase':
+        break
+
+# Recreate the frequency response and shift the phase
 (mag, phase, w) = freqresp(L, logspace(-4, 3));
 phase = phase - 360;
-subplot(212);
-semilogx([1e-4, 1e3], [-180, -180], 'k-')
-hold(True);
-semilogx(w, np.squeeze(phase), 'b-')
-axis([1e-4, 1e3, -360, 0]);
+
+# Replot the phase by hand
+ax.semilogx([1e-4, 1e3], [-180, -180], 'k-')
+ax.semilogx(w, np.squeeze(phase), 'b-')
+ax.axis([1e-4, 1e3, -360, 0]);
 xlabel('Frequency [deg]'); ylabel('Phase [deg]');
 # set(gca, 'YTick', [-360, -270, -180, -90, 0]);
 # set(gca, 'XTick', [10^-4, 10^-2, 1, 100]);
@@ -109,7 +128,6 @@
 # Nyquist plot for complete design
 #
 figure(7); clf;
-axis([-700, 5300, -3000, 3000]); hold(True);
 nyquist(L, (0.0001, 1000));
 axis([-700, 5300, -3000, 3000]);
 
@@ -118,7 +136,6 @@
 
 # Expanded region  
 figure(8); clf; subplot(231); 
-axis([-10, 5, -20, 20]); hold(True);
 nyquist(L);
 axis([-10, 5, -20, 20]);
 
@@ -136,7 +153,7 @@
 
 figure(9); 
 (Yvec, Tvec) = step(T, linspace(0, 20));
-plot(Tvec.T, Yvec.T); hold(True);
+plot(Tvec.T, Yvec.T); 
 
 (Yvec, Tvec) = step(Co*S, linspace(0, 20));
 plot(Tvec.T, Yvec.T);