diff --git a/control/freqplot.py b/control/freqplot.py index c53e83f31..a71d44cce 100644 --- 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 index c71f5cd49..bff14d241 100644 --- 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 index aa3670d8c..9e8dd44ce 100644 --- 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 index ec35f6626..c12936f42 100644 --- 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 index 934f4b5ca..e02d86352 100644 --- 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);