From 54583d52a223430c79998f8d09ebf0cfb5952b7f Mon Sep 17 00:00:00 2001 From: David de Jong Date: Wed, 23 May 2018 15:02:36 +0200 Subject: [PATCH 01/21] first sisotool chechpoint --- control/__init__.py | 1 + control/rlocus.py | 48 ++++++++++++++++++++++++++++----------------- control/sisotool.py | 10 ++++++++++ 3 files changed, 41 insertions(+), 18 deletions(-) create mode 100644 control/sisotool.py diff --git a/control/__init__.py b/control/__init__.py index 99cf9a73d..4746d28a3 100644 --- a/control/__init__.py +++ b/control/__init__.py @@ -67,6 +67,7 @@ from .canonical import * from .robust import * from .config import * +from .sisotool import * # Exceptions from .exception import * diff --git a/control/rlocus.py b/control/rlocus.py index dea9f183b..4ee55d7dc 100644 --- a/control/rlocus.py +++ b/control/rlocus.py @@ -59,7 +59,7 @@ # Main function: compute a root locus diagram def root_locus(sys, kvect=None, xlim=None, ylim=None, plotstr='-', Plot=True, - PrintGain=True, grid=False): + PrintGain=True, grid=False, sisotool = False, f=None, ax =None): """Root locus plot Calculate the root locus by finding the roots of 1+k*TF(s) @@ -103,21 +103,21 @@ def root_locus(sys, kvect=None, xlim=None, ylim=None, plotstr='-', Plot=True, # Create the Plot if Plot: - figure_number = pylab.get_fignums() - figure_title = [pylab.figure(numb).canvas.get_window_title() for numb in figure_number] - new_figure_name = "Root Locus" - rloc_num = 1 - while new_figure_name in figure_title: - new_figure_name = "Root Locus " + str(rloc_num) - rloc_num += 1 - f = pylab.figure(new_figure_name) + if sisotool == False: + figure_number = pylab.get_fignums() + figure_title = [pylab.figure(numb).canvas.get_window_title() for numb in figure_number] + new_figure_name = "Root Locus" + rloc_num = 1 + while new_figure_name in figure_title: + new_figure_name = "Root Locus " + str(rloc_num) + rloc_num += 1 + f = pylab.figure(new_figure_name) + ax = pylab.axes() - ax = pylab.axes() if PrintGain: - click_point, = ax.plot([0], [0],color='k',markersize = 0,marker='s',zorder=20) f.canvas.mpl_connect( - 'button_release_event', partial(_RLFeedbackClicks, sys=sys,fig=f,point=click_point)) + 'button_release_event', partial(_RLFeedbackClicks,sys=sys,fig=f,sisotool=sisotool)) # plot open loop poles poles = array(denp.r) @@ -344,21 +344,33 @@ def _RLSortRoots(mymat): return sorted -def _RLFeedbackClicks(event, sys,fig,point): +def _RLFeedbackClicks(event,sys,fig,sisotool): """Print root-locus gain feedback for clicks on the root-locus plot """ + (nump, denp) = _systopoly1d(sys) s = complex(event.xdata, event.ydata) K = -1./sys.horner(s) if abs(K.real) > 1e-8 and abs(K.imag/K.real) < 0.04: print("Clicked at %10.4g%+10.4gj gain %10.4g damp %10.4g" % (s.real, s.imag, K.real, -1 * s.real / abs(s))) - point.set_ydata(s.imag) - point.set_xdata(s.real) - point.set_markersize(8) fig.suptitle("Clicked at: %10.4g%+10.4gj gain: %10.4g damp: %10.4g" % - (s.real, s.imag, K.real, -1 * s.real / abs(s))) - fig.canvas.draw() + (s.real, s.imag, K.real, -1 * s.real / abs(s))) + + ax = fig.axes[0] + for line in reversed(ax.lines): + if len(line.get_xdata()) == 1: + line.remove() + del line + else: + break + + if sisotool: + mymat = _RLFindRoots(nump, denp, K.real) + ax.plot([root.real for root in mymat],[root.imag for root in mymat],'m.',marker='s',markersize=8, zorder=20) + else: + ax.plot(s.real, s.imag, 'k.', marker='s', markersize=8, zorder=20) + fig.canvas.draw() def _sgrid_func(fig=None, zeta=None, wn=None): if fig is None: diff --git a/control/sisotool.py b/control/sisotool.py new file mode 100644 index 000000000..fc5cf437b --- /dev/null +++ b/control/sisotool.py @@ -0,0 +1,10 @@ +__all__ = ['sisotool'] + +from .freqplot import bode_plot +from .rlocus import root_locus +import matplotlib.pyplot as plt + +def sisotool(sys, kvect=None,PrintGain=True,grid=False,dB=None,Hz=None,deg=None): + f, (ax1, ax2) = plt.subplots(1, 2) + root_locus(sys,ax=ax1,f=f,sisotool=True) + #bode_plot(sys,ax=ax2) \ No newline at end of file From 2e1bf7d6dc2933e95ea7d097fd7af7f74c3cc732 Mon Sep 17 00:00:00 2001 From: David de Jong Date: Thu, 24 May 2018 22:20:31 +0200 Subject: [PATCH 02/21] second sisotool checkpoint --- control/freqplot.py | 44 +++++++++++++++++------------- control/rlocus.py | 47 +++++++++++++++++++++++--------- control/sisotool.py | 65 +++++++++++++++++++++++++++++++++++++++++---- control/timeresp.py | 2 +- 4 files changed, 121 insertions(+), 37 deletions(-) diff --git a/control/freqplot.py b/control/freqplot.py index a71d44cce..73f8968c4 100644 --- a/control/freqplot.py +++ b/control/freqplot.py @@ -177,24 +177,32 @@ def bode_plot(syslist, omega=None, dB=None, Hz=None, deg=None, # 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) + # Get the current figure + if 'sisotool' in kwargs: + fig = kwargs['fig'] + ax_mag = fig.axes[0] + ax_phase = fig.axes[2] + sisotool = kwargs['sisotool'] + del kwargs['fig'] + del kwargs['sisotool'] + else: + 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 if dB: diff --git a/control/rlocus.py b/control/rlocus.py index 4ee55d7dc..7a80837dd 100644 --- a/control/rlocus.py +++ b/control/rlocus.py @@ -52,6 +52,7 @@ import pylab # plotting routines from .xferfcn import _convertToTransferFunction from .exception import ControlMIMONotImplemented +from .sisotool import _SisotoolUpdate from functools import partial __all__ = ['root_locus', 'rlocus'] @@ -59,7 +60,8 @@ # Main function: compute a root locus diagram def root_locus(sys, kvect=None, xlim=None, ylim=None, plotstr='-', Plot=True, - PrintGain=True, grid=False, sisotool = False, f=None, ax =None): + PrintGain=True, grid=False, **kwargs): + """Root locus plot Calculate the root locus by finding the roots of 1+k*TF(s) @@ -96,14 +98,16 @@ def root_locus(sys, kvect=None, xlim=None, ylim=None, plotstr='-', Plot=True, (nump, denp) = _systopoly1d(sys) if kvect is None: + start_mat = _RLFindRoots(nump, denp, [1]) kvect, mymat, xlim, ylim = _default_gains(nump, denp, xlim, ylim) else: + start_mat = _RLFindRoots(nump, denp, kvect[0]) mymat = _RLFindRoots(nump, denp, kvect) mymat = _RLSortRoots(mymat) # Create the Plot if Plot: - if sisotool == False: + if 'sisotool' not in kwargs: figure_number = pylab.get_fignums() figure_title = [pylab.figure(numb).canvas.get_window_title() for numb in figure_number] new_figure_name = "Root Locus" @@ -113,11 +117,19 @@ def root_locus(sys, kvect=None, xlim=None, ylim=None, plotstr='-', Plot=True, rloc_num += 1 f = pylab.figure(new_figure_name) ax = pylab.axes() - + else: + sisotool = kwargs['sisotool'] + f = kwargs['fig'] + ax = f.axes[1] + bode_plot_params = kwargs['bode_plot_params'] + tvect = kwargs['tvect'] + event = type('event', (object,), {'xdata': start_mat[0][0].real,'ydata':start_mat[0][0].imag})() + _RLFeedbackClicks(event, sys, f, sisotool,bode_plot_params,tvect) if PrintGain: + f.canvas.mpl_connect( - 'button_release_event', partial(_RLFeedbackClicks,sys=sys,fig=f,sisotool=sisotool)) + 'button_release_event', partial(_RLFeedbackClicks,sys=sys,fig=f,sisotool=sisotool,bode_plot_params=bode_plot_params,tvect=tvect)) # plot open loop poles poles = array(denp.r) @@ -139,7 +151,9 @@ def root_locus(sys, kvect=None, xlim=None, ylim=None, plotstr='-', Plot=True, ax.set_ylim(ylim) ax.set_xlabel('Real') ax.set_ylabel('Imaginary') - if grid: + if grid and sisotool: + _sgrid_func(f) + elif grid: _sgrid_func() return mymat, kvect @@ -344,7 +358,7 @@ def _RLSortRoots(mymat): return sorted -def _RLFeedbackClicks(event,sys,fig,sisotool): +def _RLFeedbackClicks(event,sys,fig,sisotool,bode_plot_params,tvect): """Print root-locus gain feedback for clicks on the root-locus plot """ (nump, denp) = _systopoly1d(sys) @@ -356,26 +370,33 @@ def _RLFeedbackClicks(event,sys,fig,sisotool): fig.suptitle("Clicked at: %10.4g%+10.4gj gain: %10.4g damp: %10.4g" % (s.real, s.imag, K.real, -1 * s.real / abs(s))) - ax = fig.axes[0] - for line in reversed(ax.lines): + if sisotool: + ax_rlocus = fig.axes[1] + else: + ax_rlocus = fig.axes[0] + + for line in reversed(ax_rlocus.lines): if len(line.get_xdata()) == 1: line.remove() del line - else: - break + #else: + # break if sisotool: mymat = _RLFindRoots(nump, denp, K.real) - ax.plot([root.real for root in mymat],[root.imag for root in mymat],'m.',marker='s',markersize=8, zorder=20) + ax_rlocus.plot([root.real for root in mymat],[root.imag for root in mymat],'m.',marker='s',markersize=8, zorder=20) + _SisotoolUpdate(sys,fig,K.real[0][0],bode_plot_params,tvect) else: - ax.plot(s.real, s.imag, 'k.', marker='s', markersize=8, zorder=20) + ax_rlocus.plot(s.real, s.imag, 'k.', marker='s', markersize=8, zorder=20) fig.canvas.draw() def _sgrid_func(fig=None, zeta=None, wn=None): if fig is None: fig = pylab.gcf() - ax = fig.gca() + ax = fig.gca() + else: + ax = fig.axes[1] xlocator = ax.get_xaxis().get_major_locator() ylim = ax.get_ylim() diff --git a/control/sisotool.py b/control/sisotool.py index fc5cf437b..73cba8e48 100644 --- a/control/sisotool.py +++ b/control/sisotool.py @@ -1,10 +1,65 @@ __all__ = ['sisotool'] from .freqplot import bode_plot -from .rlocus import root_locus +from .timeresp import step_response +from .lti import issiso import matplotlib.pyplot as plt -def sisotool(sys, kvect=None,PrintGain=True,grid=False,dB=None,Hz=None,deg=None): - f, (ax1, ax2) = plt.subplots(1, 2) - root_locus(sys,ax=ax1,f=f,sisotool=True) - #bode_plot(sys,ax=ax2) \ No newline at end of file +def sisotool(sys, kvect = None, xlim = None, ylim = None, plotstr_rlocus = '-',rlocus_grid = False, omega = None, dB = None, Hz = None, deg = None, omega_limits = None, omega_num = None, tvect=None): + + from .rlocus import root_locus + + # Check if it is a single SISO system + issiso(sys,strict=True) + + # Setup sisotool figure or superimpose if one is already present + fig = plt.gcf() + if fig.canvas.get_window_title() != 'Sisotool': + plt.close(fig) + fig,axes = plt.subplots(2, 2) + fig.canvas.set_window_title('Sisotool') + + # Extract bode plot parameters + bode_plot_params = { + 'omega': omega, + 'dB': dB, + 'Hz': Hz, + 'deg': deg, + 'omega_limits': omega_limits, + 'omega_num' : omega_num, + 'sisotool': True, + 'fig': fig, + } + + #To-do find out clever way to pass correct settings to other plots + _SisotoolUpdate(sys, fig, 1,bode_plot_params,tvect) + + # Setup the root-locus plot window + root_locus(sys,kvect=kvect,xlim=xlim,ylim = ylim,plotstr=plotstr_rlocus,grid = rlocus_grid,fig=fig,bode_plot_params=bode_plot_params,tvect=tvect,sisotool=True) + +def _SisotoolUpdate(sys,fig,K,bode_plot_params,tvect): + + # Get the subaxes and clear them + ax_mag,ax_rlocus,ax_phase,ax_step = fig.axes[0],fig.axes[1],fig.axes[2],fig.axes[3] + ax_mag.cla(),ax_phase.cla(),ax_step.cla() + + # Set the titles and labels + ax_mag.set_title('Bode magnitude') + ax_phase.set_title('Bode phase') + ax_rlocus.set_title('Root locus') + ax_step.set_title('Step response') + ax_step.set_xlabel('Time (seconds)') + ax_step.set_ylabel('Amplitude') + + # Update the bodeplot + bode_plot_params['syslist'] = sys*K.real + bode_plot(**bode_plot_params) + + # Generate the step response + sys_closed = (K*sys).feedback(1) + if tvect is None: + tvect, yout = step_response(sys_closed) + else: + tvect, yout = step_response(sys_closed,tvect) + ax_step.plot(tvect, yout) + diff --git a/control/timeresp.py b/control/timeresp.py index 4c137c933..82ce4c101 100644 --- a/control/timeresp.py +++ b/control/timeresp.py @@ -339,7 +339,7 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False): def _get_ss_simo(sys, input=None, output=None): """Return a SISO or SIMO state-space version of sys - If input is not specified, select first input and issue warning + If input is not sfpecified, select first input and issue warning """ sys_ss = _convertToStateSpace(sys) if sys_ss.issiso(): From f4545dd76a96098867a330b8194abb70649778ff Mon Sep 17 00:00:00 2001 From: David de Jong Date: Fri, 25 May 2018 13:14:19 +0200 Subject: [PATCH 03/21] changed matplotlib event handler --- control/rlocus.py | 49 +++++++++++++++++++++++++++++++-------------- control/sisotool.py | 8 ++++---- 2 files changed, 38 insertions(+), 19 deletions(-) diff --git a/control/rlocus.py b/control/rlocus.py index 7a80837dd..8ca60267c 100644 --- a/control/rlocus.py +++ b/control/rlocus.py @@ -124,12 +124,15 @@ def root_locus(sys, kvect=None, xlim=None, ylim=None, plotstr='-', Plot=True, bode_plot_params = kwargs['bode_plot_params'] tvect = kwargs['tvect'] event = type('event', (object,), {'xdata': start_mat[0][0].real,'ydata':start_mat[0][0].imag})() - _RLFeedbackClicks(event, sys, f, sisotool,bode_plot_params,tvect) - - if PrintGain: + _RLFeedbackClicksSisotool(event, sys, f,bode_plot_params,tvect) + if PrintGain and sisotool == False: + f.canvas.mpl_connect( + 'button_release_event', partial(_RLFeedbackClicksPoint,sys=sys,fig=f)) + elif sisotool == True: + print('this is run') f.canvas.mpl_connect( - 'button_release_event', partial(_RLFeedbackClicks,sys=sys,fig=f,sisotool=sisotool,bode_plot_params=bode_plot_params,tvect=tvect)) + 'button_release_event',partial(_RLFeedbackClicksSisotool,sys=sys, fig=f, bode_plot_params=bode_plot_params,tvect=tvect)) # plot open loop poles poles = array(denp.r) @@ -357,24 +360,35 @@ def _RLSortRoots(mymat): prevrow = sorted[n, :] return sorted +def _RLFeedbackClicksSisotool(event,sys,fig,bode_plot_params,tvect): + """Update Sisotool plots if a new point on the root locus plot is clicked + """ + s = complex(event.xdata, event.ydata) + K = -1./sys.horner(s) + ax_rlocus = fig.axes[1] + if _RLFeedbackClicksPoint(event,sys,fig,ax_rlocus): + _SisotoolUpdate(sys,fig,K.real[0][0],bode_plot_params,tvect) + -def _RLFeedbackClicks(event,sys,fig,sisotool,bode_plot_params,tvect): - """Print root-locus gain feedback for clicks on the root-locus plot +def _RLFeedbackClicksPoint(event,sys,fig,ax_rlocus=None): + """Display root-locus gain feedback point for clicks on the root-locus plot """ + if ax_rlocus is None: + ax_rlocus = fig.axes[0] + sisotool = False + else: + sisotool = True + (nump, denp) = _systopoly1d(sys) s = complex(event.xdata, event.ydata) - K = -1./sys.horner(s) - if abs(K.real) > 1e-8 and abs(K.imag/K.real) < 0.04: + K = -1. / sys.horner(s) + if abs(K.real) > 1e-8 and abs(K.imag / K.real) < 0.04: print("Clicked at %10.4g%+10.4gj gain %10.4g damp %10.4g" % (s.real, s.imag, K.real, -1 * s.real / abs(s))) fig.suptitle("Clicked at: %10.4g%+10.4gj gain: %10.4g damp: %10.4g" % (s.real, s.imag, K.real, -1 * s.real / abs(s))) - if sisotool: - ax_rlocus = fig.axes[1] - else: - ax_rlocus = fig.axes[0] - + # Remove the previous points for line in reversed(ax_rlocus.lines): if len(line.get_xdata()) == 1: line.remove() @@ -382,15 +396,20 @@ def _RLFeedbackClicks(event,sys,fig,sisotool,bode_plot_params,tvect): #else: # break + # Visualise clicked point, display all roots for sisotool mode if sisotool: mymat = _RLFindRoots(nump, denp, K.real) - ax_rlocus.plot([root.real for root in mymat],[root.imag for root in mymat],'m.',marker='s',markersize=8, zorder=20) - _SisotoolUpdate(sys,fig,K.real[0][0],bode_plot_params,tvect) + ax_rlocus.plot([root.real for root in mymat], [root.imag for root in mymat], 'm.', marker='s', markersize=8, + zorder=20) else: ax_rlocus.plot(s.real, s.imag, 'k.', marker='s', markersize=8, zorder=20) fig.canvas.draw() + return True + + + def _sgrid_func(fig=None, zeta=None, wn=None): if fig is None: fig = pylab.gcf() diff --git a/control/sisotool.py b/control/sisotool.py index 73cba8e48..bb22983f9 100644 --- a/control/sisotool.py +++ b/control/sisotool.py @@ -31,13 +31,13 @@ def sisotool(sys, kvect = None, xlim = None, ylim = None, plotstr_rlocus = '-',r 'fig': fig, } - #To-do find out clever way to pass correct settings to other plots - _SisotoolUpdate(sys, fig, 1,bode_plot_params,tvect) + # First time call to setup the bode and step response plots + _SisotoolUpdate(sys, fig,1 if kvect is None else kvect[0],bode_plot_params) # Setup the root-locus plot window root_locus(sys,kvect=kvect,xlim=xlim,ylim = ylim,plotstr=plotstr_rlocus,grid = rlocus_grid,fig=fig,bode_plot_params=bode_plot_params,tvect=tvect,sisotool=True) -def _SisotoolUpdate(sys,fig,K,bode_plot_params,tvect): +def _SisotoolUpdate(sys,fig,K,bode_plot_params,tvect=None): # Get the subaxes and clear them ax_mag,ax_rlocus,ax_phase,ax_step = fig.axes[0],fig.axes[1],fig.axes[2],fig.axes[3] @@ -55,7 +55,7 @@ def _SisotoolUpdate(sys,fig,K,bode_plot_params,tvect): bode_plot_params['syslist'] = sys*K.real bode_plot(**bode_plot_params) - # Generate the step response + # Generate the step response and plot it sys_closed = (K*sys).feedback(1) if tvect is None: tvect, yout = step_response(sys_closed) From e3c1f31de3f0715676b9b30356f327fcf4a114bf Mon Sep 17 00:00:00 2001 From: David de Jong Date: Fri, 25 May 2018 16:46:48 +0200 Subject: [PATCH 04/21] improved sisotool readability --- control/rlocus.py | 44 +++++++++++++++++++++++--------------------- control/sisotool.py | 1 + 2 files changed, 24 insertions(+), 21 deletions(-) diff --git a/control/rlocus.py b/control/rlocus.py index 8ca60267c..161da743c 100644 --- a/control/rlocus.py +++ b/control/rlocus.py @@ -105,9 +105,16 @@ def root_locus(sys, kvect=None, xlim=None, ylim=None, plotstr='-', Plot=True, mymat = _RLFindRoots(nump, denp, kvect) mymat = _RLSortRoots(mymat) + # Check for sisotool mode + sisotool = False if 'sisotool' not in kwargs else True + # Create the Plot if Plot: - if 'sisotool' not in kwargs: + if sisotool: + f = kwargs['fig'] + ax = f.axes[1] + + else: figure_number = pylab.get_fignums() figure_title = [pylab.figure(numb).canvas.get_window_title() for numb in figure_number] new_figure_name = "Root Locus" @@ -117,22 +124,16 @@ def root_locus(sys, kvect=None, xlim=None, ylim=None, plotstr='-', Plot=True, rloc_num += 1 f = pylab.figure(new_figure_name) ax = pylab.axes() - else: - sisotool = kwargs['sisotool'] - f = kwargs['fig'] - ax = f.axes[1] - bode_plot_params = kwargs['bode_plot_params'] - tvect = kwargs['tvect'] - event = type('event', (object,), {'xdata': start_mat[0][0].real,'ydata':start_mat[0][0].imag})() - _RLFeedbackClicksSisotool(event, sys, f,bode_plot_params,tvect) if PrintGain and sisotool == False: f.canvas.mpl_connect( 'button_release_event', partial(_RLFeedbackClicksPoint,sys=sys,fig=f)) + elif sisotool == True: - print('this is run') + f.axes[1].plot([root.real for root in start_mat], [root.imag for root in start_mat], 'm.', marker='s', markersize=8,zorder=20) + f.suptitle("Clicked at: %10.4g%+10.4gj gain: %10.4g damp: %10.4g" % (start_mat[0][0].real, start_mat[0][0].imag, 1, -1 * start_mat[0][0].real / abs(start_mat[0][0]))) f.canvas.mpl_connect( - 'button_release_event',partial(_RLFeedbackClicksSisotool,sys=sys, fig=f, bode_plot_params=bode_plot_params,tvect=tvect)) + 'button_release_event',partial(_RLFeedbackClicksSisotool,sys=sys, fig=f, bode_plot_params=kwargs['bode_plot_params'],tvect=kwargs['tvect'])) # plot open loop poles poles = array(denp.r) @@ -158,6 +159,9 @@ def root_locus(sys, kvect=None, xlim=None, ylim=None, plotstr='-', Plot=True, _sgrid_func(f) elif grid: _sgrid_func() + else: + ax.axhline(0., linestyle=':', color='k') + ax.axvline(0., linestyle=':', color='k') return mymat, kvect @@ -366,23 +370,22 @@ def _RLFeedbackClicksSisotool(event,sys,fig,bode_plot_params,tvect): s = complex(event.xdata, event.ydata) K = -1./sys.horner(s) ax_rlocus = fig.axes[1] - if _RLFeedbackClicksPoint(event,sys,fig,ax_rlocus): + if _RLFeedbackClicksPoint(event,sys,fig,ax_rlocus,sisotool=True): _SisotoolUpdate(sys,fig,K.real[0][0],bode_plot_params,tvect) - -def _RLFeedbackClicksPoint(event,sys,fig,ax_rlocus=None): +def _RLFeedbackClicksPoint(event,sys,fig,ax_rlocus=None,sisotool=False): """Display root-locus gain feedback point for clicks on the root-locus plot """ - if ax_rlocus is None: + if sisotool == False: ax_rlocus = fig.axes[0] - sisotool = False - else: - sisotool = True (nump, denp) = _systopoly1d(sys) s = complex(event.xdata, event.ydata) K = -1. / sys.horner(s) - if abs(K.real) > 1e-8 and abs(K.imag / K.real) < 0.04: + print(K) + if abs(K.real) > 1e-8 and abs(K.imag / K.real) < 0.04 and event.inaxes == ax_rlocus.axes: + + # Display the parameters in the output window and figure print("Clicked at %10.4g%+10.4gj gain %10.4g damp %10.4g" % (s.real, s.imag, K.real, -1 * s.real / abs(s))) fig.suptitle("Clicked at: %10.4g%+10.4gj gain: %10.4g damp: %10.4g" % @@ -393,8 +396,6 @@ def _RLFeedbackClicksPoint(event,sys,fig,ax_rlocus=None): if len(line.get_xdata()) == 1: line.remove() del line - #else: - # break # Visualise clicked point, display all roots for sisotool mode if sisotool: @@ -404,6 +405,7 @@ def _RLFeedbackClicksPoint(event,sys,fig,ax_rlocus=None): else: ax_rlocus.plot(s.real, s.imag, 'k.', marker='s', markersize=8, zorder=20) + # Update the canvas fig.canvas.draw() return True diff --git a/control/sisotool.py b/control/sisotool.py index bb22983f9..e7736a5ad 100644 --- a/control/sisotool.py +++ b/control/sisotool.py @@ -62,4 +62,5 @@ def _SisotoolUpdate(sys,fig,K,bode_plot_params,tvect=None): else: tvect, yout = step_response(sys_closed,tvect) ax_step.plot(tvect, yout) + ax_step.axhline(1.,linestyle=':',color='k') From bde60f58833e8b4696054227c67e46e02f3aa4d4 Mon Sep 17 00:00:00 2001 From: David de Jong Date: Mon, 28 May 2018 15:06:27 +0200 Subject: [PATCH 05/21] Added display of GM and PM in bode plot. Fixed Hz bug support. Fixed phase limit bug. --- control/freqplot.py | 26 +++++++++++++++++++++----- control/rlocus.py | 1 - control/sisotool.py | 3 ++- 3 files changed, 23 insertions(+), 7 deletions(-) diff --git a/control/freqplot.py b/control/freqplot.py index ad791487e..390b3be6d 100644 --- a/control/freqplot.py +++ b/control/freqplot.py @@ -86,7 +86,7 @@ def bode_plot(syslist, omega=None, dB=None, Hz=None, deg=None, omega_num: int number of samples margins : boolean - if True, plot gain and phase margin + If True, plot gain and phase margin *args, **kwargs: Additional options to matplotlib (color, linestyle, etc) @@ -233,10 +233,14 @@ def bode_plot(syslist, omega=None, dB=None, Hz=None, deg=None, if margins: margin = stability_margins(sys) gm, pm, Wcg, Wcp = margin[0], margin[1], margin[3], margin[4] - if pm >= 0.: - phase_limit = -180. - else: + phase_at_cp = phases[0][(np.abs(omegas[0] - Wcp)).argmin()] + if phase_at_cp >= 0.: phase_limit = 180. + else: + phase_limit = -180. + + if Hz: + Wcg, Wcp = Wcg/(2*math.pi),Wcp/(2*math.pi) ax_mag.axhline(y=0 if dB else 1, color='k', linestyle=':') ax_phase.axhline(y=phase_limit if deg else math.radians(phase_limit), color='k', linestyle=':') @@ -271,7 +275,14 @@ def bode_plot(syslist, omega=None, dB=None, Hz=None, deg=None, ax_mag.set_ylim(mag_ylim) ax_phase.set_ylim(phase_ylim) - plt.suptitle('Gm = %.2f %s(at %.2f rad/s), Pm = %.2f %s (at %.2f rad/s)'%(20*np.log10(gm) if dB else gm,'dB ' if dB else '\b',Wcg,pm if deg else math.radians(pm),'deg' if deg else 'rad',Wcp)) + + if sisotool: + ax_mag.text(0.04, 0.06, 'G.M.: %.2f %s\nFreq: %.2f %s'%(20*np.log10(gm) if dB else gm,'dB ' if dB else '\b',Wcg,'Hz' if Hz else 'rad/s'), horizontalalignment='left', verticalalignment='bottom', + transform=ax_mag.transAxes,fontsize=10) + ax_phase.text(0.04, 0.06, 'P.M.: %.2f %s\nFreq: %.2f %s'%(pm if deg else math.radians(pm),'deg' if deg else 'rad',Wcp,'Hz' if Hz else 'rad/s'), horizontalalignment='left', verticalalignment='bottom', + transform=ax_phase.transAxes,fontsize=10) + else: + plt.suptitle('Gm = %.2f %s(at %.2f %s), Pm = %.2f %s (at %.2f %s)'%(20*np.log10(gm) if dB else gm,'dB ' if dB else '\b','Hz' if Hz else 'rad/s',Wcg,pm if deg else math.radians(pm),'deg' if deg else 'rad',Wcp,'Hz' if Hz else 'rad/s')) if nyquistfrq_plot: ax_phase.axvline(nyquistfrq_plot, color=pltline[0].get_color()) @@ -644,6 +655,11 @@ def gen_prefix(pow1000): 'z', # zepto (10^-21) 'y'][8 - pow1000] # yocto (10^-24) +def find_nearest_omega(omega_list, omega): + omega_list = np.asarray(omega_list) + idx = (np.abs(omega_list - omega)).argmin() + return omega_list[(np.abs(omega_list - omega)).argmin()] + # Function aliases bode = bode_plot nyquist = nyquist_plot diff --git a/control/rlocus.py b/control/rlocus.py index 161da743c..f32287ced 100644 --- a/control/rlocus.py +++ b/control/rlocus.py @@ -382,7 +382,6 @@ def _RLFeedbackClicksPoint(event,sys,fig,ax_rlocus=None,sisotool=False): (nump, denp) = _systopoly1d(sys) s = complex(event.xdata, event.ydata) K = -1. / sys.horner(s) - print(K) if abs(K.real) > 1e-8 and abs(K.imag / K.real) < 0.04 and event.inaxes == ax_rlocus.axes: # Display the parameters in the output window and figure diff --git a/control/sisotool.py b/control/sisotool.py index e7736a5ad..558837712 100644 --- a/control/sisotool.py +++ b/control/sisotool.py @@ -5,7 +5,7 @@ from .lti import issiso import matplotlib.pyplot as plt -def sisotool(sys, kvect = None, xlim = None, ylim = None, plotstr_rlocus = '-',rlocus_grid = False, omega = None, dB = None, Hz = None, deg = None, omega_limits = None, omega_num = None, tvect=None): +def sisotool(sys, kvect = None, xlim = None, ylim = None, plotstr_rlocus = '-',rlocus_grid = False, omega = None, dB = None, Hz = None, deg = None, omega_limits = None, omega_num = None,margins = True, tvect=None): from .rlocus import root_locus @@ -29,6 +29,7 @@ def sisotool(sys, kvect = None, xlim = None, ylim = None, plotstr_rlocus = '-',r 'omega_num' : omega_num, 'sisotool': True, 'fig': fig, + 'margins': margins, } # First time call to setup the bode and step response plots From bc74a0a01e2eb075e83211e5a85a08c19a86df81 Mon Sep 17 00:00:00 2001 From: David de Jong Date: Tue, 29 May 2018 15:47:36 +0200 Subject: [PATCH 06/21] Added support for both matplotlib 1.x.x. and 2.x.x and fixed redrawing issue. --- control/freqplot.py | 39 ++++++++++++++++++------------------ control/rlocus.py | 20 +++++++++++-------- control/sisotool.py | 48 ++++++++++++++++++++++++++++++++++----------- 3 files changed, 69 insertions(+), 38 deletions(-) diff --git a/control/freqplot.py b/control/freqplot.py index 390b3be6d..9112b7abf 100644 --- a/control/freqplot.py +++ b/control/freqplot.py @@ -40,7 +40,7 @@ # SUCH DAMAGE. # # $Id$ - +import matplotlib import matplotlib.pyplot as plt import scipy as sp import numpy as np @@ -192,6 +192,7 @@ def bode_plot(syslist, omega=None, dB=None, Hz=None, deg=None, fig = plt.gcf() ax_mag = None ax_phase = None + sisotool = False # Get the current axes if they already exist for ax in fig.axes: @@ -242,47 +243,47 @@ def bode_plot(syslist, omega=None, dB=None, Hz=None, deg=None, if Hz: Wcg, Wcp = Wcg/(2*math.pi),Wcp/(2*math.pi) - ax_mag.axhline(y=0 if dB else 1, color='k', linestyle=':') - ax_phase.axhline(y=phase_limit if deg else math.radians(phase_limit), color='k', linestyle=':') + ax_mag.axhline(y=0 if dB else 1, color='k', linestyle=':',zorder=-20) + ax_phase.axhline(y=phase_limit if deg else math.radians(phase_limit), color='k', linestyle=':',zorder=-20) mag_ylim = ax_mag.get_ylim() phase_ylim = ax_phase.get_ylim() if pm != float('inf') and Wcp != float('nan'): if dB: - ax_mag.semilogx([Wcp, Wcp], [0.,-1e5],color='k', linestyle=':') + ax_mag.semilogx([Wcp, Wcp], [0.,-1e5],color='k', linestyle=':',zorder=-20) else: - ax_mag.loglog([Wcp,Wcp], [1.,1e-8],color='k',linestyle=':') + ax_mag.loglog([Wcp,Wcp], [1.,1e-8],color='k',linestyle=':',zorder=-20) if deg: - ax_phase.semilogx([Wcp, Wcp], [1e5, phase_limit+pm],color='k', linestyle=':') - ax_phase.semilogx([Wcp, Wcp], [phase_limit + pm, phase_limit],color='k') + ax_phase.semilogx([Wcp, Wcp], [1e5, phase_limit+pm],color='k', linestyle=':',zorder=-20) + ax_phase.semilogx([Wcp, Wcp], [phase_limit + pm, phase_limit],color='k',zorder=-20) else: - ax_phase.semilogx([Wcp, Wcp], [1e5, math.radians(phase_limit)+math.radians(pm)],color='k', linestyle=':') - ax_phase.semilogx([Wcp, Wcp], [math.radians(phase_limit) +math.radians(pm), math.radians(phase_limit)],color='k') + ax_phase.semilogx([Wcp, Wcp], [1e5, math.radians(phase_limit)+math.radians(pm)],color='k', linestyle=':',zorder=-20) + ax_phase.semilogx([Wcp, Wcp], [math.radians(phase_limit) +math.radians(pm), math.radians(phase_limit)],color='k',zorder=-20) if gm != float('inf') and Wcg != float('nan'): if dB: - ax_mag.semilogx([Wcg, Wcg], [-20.*np.log10(gm), -1e5],color='k', linestyle=':') - ax_mag.semilogx([Wcg, Wcg], [0,-20*np.log10(gm)],color='k') + ax_mag.semilogx([Wcg, Wcg], [-20.*np.log10(gm), -1e5],color='k', linestyle=':',zorder=-20) + ax_mag.semilogx([Wcg, Wcg], [0,-20*np.log10(gm)],color='k',zorder=-20) else: - ax_mag.loglog([Wcg, Wcg], [1./gm,1e-8],color='k', linestyle=':') - ax_mag.loglog([Wcg, Wcg], [1.,1./gm],color='k') + ax_mag.loglog([Wcg, Wcg], [1./gm,1e-8],color='k', linestyle=':',zorder=-20) + ax_mag.loglog([Wcg, Wcg], [1.,1./gm],color='k',zorder=-20) if deg: - ax_phase.semilogx([Wcg, Wcg], [1e-8, phase_limit],color='k', linestyle=':') + ax_phase.semilogx([Wcg, Wcg], [1e-8, phase_limit],color='k', linestyle=':',zorder=-20) else: - ax_phase.semilogx([Wcg, Wcg], [1e-8, math.radians(phase_limit)],color='k', linestyle=':') + ax_phase.semilogx([Wcg, Wcg], [1e-8, math.radians(phase_limit)],color='k', linestyle=':',zorder=-20) ax_mag.set_ylim(mag_ylim) ax_phase.set_ylim(phase_ylim) if sisotool: - ax_mag.text(0.04, 0.06, 'G.M.: %.2f %s\nFreq: %.2f %s'%(20*np.log10(gm) if dB else gm,'dB ' if dB else '\b',Wcg,'Hz' if Hz else 'rad/s'), horizontalalignment='left', verticalalignment='bottom', - transform=ax_mag.transAxes,fontsize=10) + ax_mag.text(0.04, 0.06, 'G.M.: %.2f %s\nFreq: %.2f %s'%(20*np.log10(gm) if dB else gm,'dB ' if dB else '',Wcg,'Hz' if Hz else 'rad/s'), horizontalalignment='left', verticalalignment='bottom', + transform=ax_mag.transAxes,fontsize=8 if int(matplotlib.__version__[0]) == 1 else 6) ax_phase.text(0.04, 0.06, 'P.M.: %.2f %s\nFreq: %.2f %s'%(pm if deg else math.radians(pm),'deg' if deg else 'rad',Wcp,'Hz' if Hz else 'rad/s'), horizontalalignment='left', verticalalignment='bottom', - transform=ax_phase.transAxes,fontsize=10) + transform=ax_phase.transAxes,fontsize=8 if int(matplotlib.__version__[0]) == 1 else 6) else: - plt.suptitle('Gm = %.2f %s(at %.2f %s), Pm = %.2f %s (at %.2f %s)'%(20*np.log10(gm) if dB else gm,'dB ' if dB else '\b','Hz' if Hz else 'rad/s',Wcg,pm if deg else math.radians(pm),'deg' if deg else 'rad',Wcp,'Hz' if Hz else 'rad/s')) + plt.suptitle('Gm = %.2f %s(at %.2f %s), Pm = %.2f %s (at %.2f %s)'%(20*np.log10(gm) if dB else gm,'dB ' if dB else '\b',Wcg,'Hz' if Hz else 'rad/s',pm if deg else math.radians(pm),'deg' if deg else 'rad',Wcp,'Hz' if Hz else 'rad/s')) if nyquistfrq_plot: ax_phase.axvline(nyquistfrq_plot, color=pltline[0].get_color()) diff --git a/control/rlocus.py b/control/rlocus.py index f32287ced..dbd778c93 100644 --- a/control/rlocus.py +++ b/control/rlocus.py @@ -47,6 +47,7 @@ # Packages used by this module import numpy as np +import matplotlib from scipy import array, poly1d, row_stack, zeros_like, real, imag import scipy.signal # signal processing toolbox import pylab # plotting routines @@ -131,7 +132,7 @@ def root_locus(sys, kvect=None, xlim=None, ylim=None, plotstr='-', Plot=True, elif sisotool == True: f.axes[1].plot([root.real for root in start_mat], [root.imag for root in start_mat], 'm.', marker='s', markersize=8,zorder=20) - f.suptitle("Clicked at: %10.4g%+10.4gj gain: %10.4g damp: %10.4g" % (start_mat[0][0].real, start_mat[0][0].imag, 1, -1 * start_mat[0][0].real / abs(start_mat[0][0]))) + f.suptitle("Clicked at: %10.4g%+10.4gj gain: %10.4g damp: %10.4g" % (start_mat[0][0].real, start_mat[0][0].imag, 1, -1 * start_mat[0][0].real / abs(start_mat[0][0])),fontsize = 12 if int(matplotlib.__version__[0]) == 1 else 10) f.canvas.mpl_connect( 'button_release_event',partial(_RLFeedbackClicksSisotool,sys=sys, fig=f, bode_plot_params=kwargs['bode_plot_params'],tvect=kwargs['tvect'])) @@ -160,7 +161,7 @@ def root_locus(sys, kvect=None, xlim=None, ylim=None, plotstr='-', Plot=True, elif grid: _sgrid_func() else: - ax.axhline(0., linestyle=':', color='k') + ax.axhline(0., linestyle=':', color='k',zorder=-20) ax.axvline(0., linestyle=':', color='k') return mymat, kvect @@ -367,11 +368,14 @@ def _RLSortRoots(mymat): def _RLFeedbackClicksSisotool(event,sys,fig,bode_plot_params,tvect): """Update Sisotool plots if a new point on the root locus plot is clicked """ - s = complex(event.xdata, event.ydata) - K = -1./sys.horner(s) - ax_rlocus = fig.axes[1] - if _RLFeedbackClicksPoint(event,sys,fig,ax_rlocus,sisotool=True): - _SisotoolUpdate(sys,fig,K.real[0][0],bode_plot_params,tvect) + try: + s = complex(event.xdata, event.ydata) + K = -1. / sys.horner(s) + ax_rlocus = fig.axes[1] + if _RLFeedbackClicksPoint(event, sys, fig, ax_rlocus, sisotool=True): + _SisotoolUpdate(sys, fig, K.real[0][0], bode_plot_params, tvect) + except TypeError: + pass def _RLFeedbackClicksPoint(event,sys,fig,ax_rlocus=None,sisotool=False): """Display root-locus gain feedback point for clicks on the root-locus plot @@ -388,7 +392,7 @@ def _RLFeedbackClicksPoint(event,sys,fig,ax_rlocus=None,sisotool=False): print("Clicked at %10.4g%+10.4gj gain %10.4g damp %10.4g" % (s.real, s.imag, K.real, -1 * s.real / abs(s))) fig.suptitle("Clicked at: %10.4g%+10.4gj gain: %10.4g damp: %10.4g" % - (s.real, s.imag, K.real, -1 * s.real / abs(s))) + (s.real, s.imag, K.real, -1 * s.real / abs(s)),fontsize = 12 if int(matplotlib.__version__[0]) == 1 else 10) # Remove the previous points for line in reversed(ax_rlocus.lines): diff --git a/control/sisotool.py b/control/sisotool.py index 558837712..4de675054 100644 --- a/control/sisotool.py +++ b/control/sisotool.py @@ -3,9 +3,10 @@ from .freqplot import bode_plot from .timeresp import step_response from .lti import issiso +import matplotlib import matplotlib.pyplot as plt -def sisotool(sys, kvect = None, xlim = None, ylim = None, plotstr_rlocus = '-',rlocus_grid = False, omega = None, dB = None, Hz = None, deg = None, omega_limits = None, omega_num = None,margins = True, tvect=None): +def sisotool(sys, kvect = None, xlim = None, ylim = None, plotstr_rlocus = 'b' if int(matplotlib.__version__[0]) == 1 else 'C0',rlocus_grid = False, omega = None, dB = None, Hz = None, deg = None, omega_limits = None, omega_num = None,margins = True, tvect=None): from .rlocus import root_locus @@ -29,7 +30,7 @@ def sisotool(sys, kvect = None, xlim = None, ylim = None, plotstr_rlocus = '-',r 'omega_num' : omega_num, 'sisotool': True, 'fig': fig, - 'margins': margins, + 'margins': margins } # First time call to setup the bode and step response plots @@ -40,22 +41,43 @@ def sisotool(sys, kvect = None, xlim = None, ylim = None, plotstr_rlocus = '-',r def _SisotoolUpdate(sys,fig,K,bode_plot_params,tvect=None): + if int(matplotlib.__version__[0]) == 1: + title_font_size = 12 + label_font_size = 10 + else: + title_font_size = 10 + label_font_size = 8 + # Get the subaxes and clear them ax_mag,ax_rlocus,ax_phase,ax_step = fig.axes[0],fig.axes[1],fig.axes[2],fig.axes[3] ax_mag.cla(),ax_phase.cla(),ax_step.cla() - # Set the titles and labels - ax_mag.set_title('Bode magnitude') - ax_phase.set_title('Bode phase') - ax_rlocus.set_title('Root locus') - ax_step.set_title('Step response') - ax_step.set_xlabel('Time (seconds)') - ax_step.set_ylabel('Amplitude') - # Update the bodeplot bode_plot_params['syslist'] = sys*K.real bode_plot(**bode_plot_params) + # Set the titles and labels + ax_mag.set_title('Bode magnitude',fontsize = title_font_size) + ax_mag.set_ylabel(ax_mag.get_ylabel(), fontsize=label_font_size) + + ax_phase.set_title('Bode phase',fontsize=title_font_size) + ax_phase.set_xlabel(ax_phase.get_xlabel(),fontsize=label_font_size) + ax_phase.set_ylabel(ax_phase.get_ylabel(),fontsize=label_font_size) + ax_phase.get_xaxis().set_label_coords(0.5, -0.15) + ax_phase.get_shared_x_axes().join(ax_phase, ax_mag) + + ax_step.set_title('Step response',fontsize = title_font_size) + ax_step.set_xlabel('Time (seconds)',fontsize=label_font_size) + ax_step.set_ylabel('Amplitude',fontsize=label_font_size) + ax_step.get_xaxis().set_label_coords(0.5, -0.15) + ax_step.get_yaxis().set_label_coords(-0.15, 0.5) + + ax_rlocus.set_title('Root locus',fontsize = title_font_size) + ax_rlocus.set_ylabel('Imag', fontsize=label_font_size) + ax_rlocus.set_xlabel('Real', fontsize=label_font_size) + ax_rlocus.get_xaxis().set_label_coords(0.5, -0.15) + ax_rlocus.get_yaxis().set_label_coords(-0.15, 0.5) + # Generate the step response and plot it sys_closed = (K*sys).feedback(1) if tvect is None: @@ -63,5 +85,9 @@ def _SisotoolUpdate(sys,fig,K,bode_plot_params,tvect=None): else: tvect, yout = step_response(sys_closed,tvect) ax_step.plot(tvect, yout) - ax_step.axhline(1.,linestyle=':',color='k') + ax_step.axhline(1.,linestyle=':',color='k',zorder=-20) + + # Manually adjust the spacing and draw the canvas + fig.subplots_adjust(top=0.9,wspace = 0.3,hspace=0.35) + fig.canvas.draw() From 1e3a925ccc804348c52a1f50fda460df30d2240d Mon Sep 17 00:00:00 2001 From: David de Jong Date: Tue, 29 May 2018 22:06:34 +0200 Subject: [PATCH 07/21] Added docstring and unittest --- control/freqplot.py | 2 + control/rlocus.py | 3 +- control/sisotool.py | 63 ++++++++++++++++++++++++--- control/tests/sisotool_test.py | 78 ++++++++++++++++++++++++++++++++++ 4 files changed, 138 insertions(+), 8 deletions(-) create mode 100644 control/tests/sisotool_test.py diff --git a/control/freqplot.py b/control/freqplot.py index 9112b7abf..44b22084d 100644 --- a/control/freqplot.py +++ b/control/freqplot.py @@ -213,8 +213,10 @@ def bode_plot(syslist, omega=None, dB=None, Hz=None, deg=None, 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()) diff --git a/control/rlocus.py b/control/rlocus.py index dbd778c93..8ba9bca8a 100644 --- a/control/rlocus.py +++ b/control/rlocus.py @@ -380,6 +380,7 @@ def _RLFeedbackClicksSisotool(event,sys,fig,bode_plot_params,tvect): def _RLFeedbackClicksPoint(event,sys,fig,ax_rlocus=None,sisotool=False): """Display root-locus gain feedback point for clicks on the root-locus plot """ + print(event) if sisotool == False: ax_rlocus = fig.axes[0] @@ -413,8 +414,6 @@ def _RLFeedbackClicksPoint(event,sys,fig,ax_rlocus=None,sisotool=False): return True - - def _sgrid_func(fig=None, zeta=None, wn=None): if fig is None: fig = pylab.gcf() diff --git a/control/sisotool.py b/control/sisotool.py index 4de675054..916eef1cd 100644 --- a/control/sisotool.py +++ b/control/sisotool.py @@ -5,9 +5,56 @@ from .lti import issiso import matplotlib import matplotlib.pyplot as plt - -def sisotool(sys, kvect = None, xlim = None, ylim = None, plotstr_rlocus = 'b' if int(matplotlib.__version__[0]) == 1 else 'C0',rlocus_grid = False, omega = None, dB = None, Hz = None, deg = None, omega_limits = None, omega_num = None,margins = True, tvect=None): - +import warnings + +def sisotool(sys, kvect = None, xlim_rlocus = None, ylim_rlocus = None, plotstr_rlocus = 'b' if int(matplotlib.__version__[0]) == 1 else 'C0',rlocus_grid = False, omega = None, dB = None, Hz = None, deg = None, omega_limits = None, omega_num = None,margins_bode = True, tvect=None): + + """Sisotool + + Sisotool style collection of plots inspired by the matlab sisotool. + The left two plots contain the bode magnitude and phase diagrams. + The top right plot is a clickable root locus plot, clicking on the + root locus will change the gain of the system. The bottom left plot + shows a closed loop time response. + + Parameters + ---------- + sys : LTI object + Linear input/output systems (SISO only) + kvect : list or ndarray, optional + List of gains to use for plotting root locus + xlim_rlocus : tuple or list, optional + control of x-axis range, normally with tuple (see matplotlib.axes) + ylim_rlocus : tuple or list, optional + control of y-axis range + plotstr_rlocus : Additional options to matplotlib + plotting style for the root locus plot(color, linestyle, etc) + rlocus_grid: boolean (default = False) + If True plot s-plane grid. + omega : freq_range + Range of frequencies in rad/sec for the bode plot + dB : boolean + If True, plot result in dB for the bode plot + Hz : boolean + If True, plot frequency in Hz for the bode plot (omega must be provided in rad/sec) + deg : boolean + If True, plot phase in degrees for the bode plot (else radians) + omega_limits: tuple, list, ... of two values + Limits of the to generate frequency vector. + If Hz=True the limits are in Hz otherwise in rad/s. + omega_num: int + number of samples + margins_bode : boolean + If True, plot gain and phase margin in the bode plot + tvect : list or ndarray, optional + List of timesteps to use for closed loop step response + + Examples + -------- + >>> sys = tf([1000], [1,25,100,0]) + >>> sisotool(sys) + + """ from .rlocus import root_locus # Check if it is a single SISO system @@ -30,14 +77,14 @@ def sisotool(sys, kvect = None, xlim = None, ylim = None, plotstr_rlocus = 'b' i 'omega_num' : omega_num, 'sisotool': True, 'fig': fig, - 'margins': margins + 'margins': margins_bode } # First time call to setup the bode and step response plots _SisotoolUpdate(sys, fig,1 if kvect is None else kvect[0],bode_plot_params) # Setup the root-locus plot window - root_locus(sys,kvect=kvect,xlim=xlim,ylim = ylim,plotstr=plotstr_rlocus,grid = rlocus_grid,fig=fig,bode_plot_params=bode_plot_params,tvect=tvect,sisotool=True) + root_locus(sys,kvect=kvect,xlim=xlim_rlocus,ylim = ylim_rlocus,plotstr=plotstr_rlocus,grid = rlocus_grid,fig=fig,bode_plot_params=bode_plot_params,tvect=tvect,sisotool=True) def _SisotoolUpdate(sys,fig,K,bode_plot_params,tvect=None): @@ -50,7 +97,11 @@ def _SisotoolUpdate(sys,fig,K,bode_plot_params,tvect=None): # Get the subaxes and clear them ax_mag,ax_rlocus,ax_phase,ax_step = fig.axes[0],fig.axes[1],fig.axes[2],fig.axes[3] - ax_mag.cla(),ax_phase.cla(),ax_step.cla() + + # Catch matplotlib 2.1.x and higher userwarnings when clearing a log axis + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + ax_step.clear(), ax_mag.clear(), ax_phase.clear() # Update the bodeplot bode_plot_params['syslist'] = sys*K.real diff --git a/control/tests/sisotool_test.py b/control/tests/sisotool_test.py new file mode 100644 index 000000000..69a4c0b3c --- /dev/null +++ b/control/tests/sisotool_test.py @@ -0,0 +1,78 @@ +import unittest +import numpy as np +from control.sisotool import sisotool +from control.tests.margin_test import assert_array_almost_equal +from control.rlocus import _RLFeedbackClicksSisotool +from control.xferfcn import TransferFunction +import matplotlib.pyplot as plt + +class TestSisotool(unittest.TestCase): + """These are tests for the sisotool in sisotool.py.""" + + def setUp(self): + """This contains some random LTI systems and scalars for testing.""" + + # Two random SISO systems. + sys1 = TransferFunction([1000],[1,25,100,0]) + sys2 = TransferFunction([1,1,1],[1,1]) + self.systems = (sys1, sys2) + + def test_sisotool(self): + sisotool(self.systems[0]) + fig = plt.gcf() + ax_mag,ax_rlocus,ax_phase,ax_step = fig.axes[0],fig.axes[1],fig.axes[2],fig.axes[3] + + # Check the initial root locus plot points + initial_point_0 = (np.array([-22.53155977]),np.array([0.])) + initial_point_1 = (np.array([-1.23422011]), np.array([-6.54667031])) + initial_point_2 = (np.array([-1.23422011]), np.array([06.54667031])) + assert_array_almost_equal(ax_rlocus.lines[0].get_data(),initial_point_0) + assert_array_almost_equal(ax_rlocus.lines[1].get_data(),initial_point_1) + assert_array_almost_equal(ax_rlocus.lines[2].get_data(),initial_point_2) + + # Check the bode plot magnitude line + bode_mag_original = np.array([ 15.12670407, 12.48358602, 10.28367521, 8.44961576, 6.91740918, 5.63436301, 4.55750665, 3.65237575, 2.89200358, 2.25587634]) + assert_array_almost_equal(ax_mag.lines[0].get_data()[1][10:20],bode_mag_original) + + # Check the step response before moving the point + print(ax_step.lines[0].get_data()[1][:10]) + step_response_original = np.array([ 0., 0.02233651, 0.13118374, 0.33078542, 0.5907113, 0.87041549, 1.13038536, 1.33851053, 1.47374666, 1.52757114]) + assert_array_almost_equal(ax_step.lines[0].get_data()[1][:10],step_response_original) + + bode_plot_params = { + 'omega': None, + 'dB': False, + 'Hz': False, + 'deg': True, + 'omega_limits': None, + 'omega_num': None, + 'sisotool': True, + 'fig': fig, + 'margins': True + } + + # Move the rootlocus to another point + event = type('test', (object,), {'xdata': 2.31206868287,'ydata':15.5983051046, 'inaxes':ax_rlocus.axes})() + _RLFeedbackClicksSisotool(event=event, sys=self.systems[0], fig=fig, bode_plot_params=bode_plot_params, tvect=None) + + # Check the moved root locus plot points + moved_point_0 = (np.array([-29.91742755]), np.array([0.])) + moved_point_1 = (np.array([2.45871378]), np.array([-15.52647768])) + moved_point_2 = (np.array([2.45871378]), np.array([15.52647768])) + assert_array_almost_equal(ax_rlocus.lines[-3].get_data(),moved_point_0) + assert_array_almost_equal(ax_rlocus.lines[-2].get_data(),moved_point_1) + assert_array_almost_equal(ax_rlocus.lines[-1].get_data(),moved_point_2) + + # Check if the bode_mag line has moved + bode_mag_moved = np.array([ 111.83321224, 92.29238035, 76.02822315, 62.46884113, 51.14108703, 41.6554004, 33.69409534, 27.00237344, 21.38086717, 16.67791585]) + assert_array_almost_equal(ax_mag.lines[0].get_data()[1][10:20],bode_mag_moved) + + # Check if the step response has changed + step_response_moved = np.array([[ 0., 0.02458187, 0.16529784 , 0.46602716 , 0.91012035 , 1.43364313, 1.93996334 , 2.3190105 , 2.47041552 , 2.32724853] ]) + assert_array_almost_equal(ax_step.lines[0].get_data()[1][:10],step_response_moved) + +def test_suite(): + return unittest.TestLoader().loadTestsFromTestCase(TestSisotool) + +if __name__ == "__main__": + unittest.main() From 94a0ea154a7c0b6a3a0b576615aa89f73c37f5be Mon Sep 17 00:00:00 2001 From: David de Jong Date: Wed, 30 May 2018 13:04:06 +0200 Subject: [PATCH 08/21] removed bug where a pole or zero would be deleted if only 1 of them was present --- control/rlocus.py | 9 ++++----- control/tests/sisotool_test.py | 12 ++++-------- 2 files changed, 8 insertions(+), 13 deletions(-) diff --git a/control/rlocus.py b/control/rlocus.py index 8ba9bca8a..3762c12e2 100644 --- a/control/rlocus.py +++ b/control/rlocus.py @@ -131,7 +131,7 @@ def root_locus(sys, kvect=None, xlim=None, ylim=None, plotstr='-', Plot=True, 'button_release_event', partial(_RLFeedbackClicksPoint,sys=sys,fig=f)) elif sisotool == True: - f.axes[1].plot([root.real for root in start_mat], [root.imag for root in start_mat], 'm.', marker='s', markersize=8,zorder=20) + f.axes[1].plot([root.real for root in start_mat], [root.imag for root in start_mat], 'm.', marker='s', markersize=8,zorder=20,label='gain_point') f.suptitle("Clicked at: %10.4g%+10.4gj gain: %10.4g damp: %10.4g" % (start_mat[0][0].real, start_mat[0][0].imag, 1, -1 * start_mat[0][0].real / abs(start_mat[0][0])),fontsize = 12 if int(matplotlib.__version__[0]) == 1 else 10) f.canvas.mpl_connect( 'button_release_event',partial(_RLFeedbackClicksSisotool,sys=sys, fig=f, bode_plot_params=kwargs['bode_plot_params'],tvect=kwargs['tvect'])) @@ -380,7 +380,6 @@ def _RLFeedbackClicksSisotool(event,sys,fig,bode_plot_params,tvect): def _RLFeedbackClicksPoint(event,sys,fig,ax_rlocus=None,sisotool=False): """Display root-locus gain feedback point for clicks on the root-locus plot """ - print(event) if sisotool == False: ax_rlocus = fig.axes[0] @@ -397,7 +396,7 @@ def _RLFeedbackClicksPoint(event,sys,fig,ax_rlocus=None,sisotool=False): # Remove the previous points for line in reversed(ax_rlocus.lines): - if len(line.get_xdata()) == 1: + if len(line.get_xdata()) == 1 and line.get_label()=='gain_point': line.remove() del line @@ -405,9 +404,9 @@ def _RLFeedbackClicksPoint(event,sys,fig,ax_rlocus=None,sisotool=False): if sisotool: mymat = _RLFindRoots(nump, denp, K.real) ax_rlocus.plot([root.real for root in mymat], [root.imag for root in mymat], 'm.', marker='s', markersize=8, - zorder=20) + zorder=20,label='gain_point') else: - ax_rlocus.plot(s.real, s.imag, 'k.', marker='s', markersize=8, zorder=20) + ax_rlocus.plot(s.real, s.imag, 'k.', marker='s', markersize=8, zorder=20,label='gain_point') # Update the canvas fig.canvas.draw() diff --git a/control/tests/sisotool_test.py b/control/tests/sisotool_test.py index 69a4c0b3c..8db6fe9a6 100644 --- a/control/tests/sisotool_test.py +++ b/control/tests/sisotool_test.py @@ -10,15 +10,11 @@ class TestSisotool(unittest.TestCase): """These are tests for the sisotool in sisotool.py.""" def setUp(self): - """This contains some random LTI systems and scalars for testing.""" - - # Two random SISO systems. - sys1 = TransferFunction([1000],[1,25,100,0]) - sys2 = TransferFunction([1,1,1],[1,1]) - self.systems = (sys1, sys2) + # One random SISO system. + self.system = TransferFunction([1000],[1,25,100,0]) def test_sisotool(self): - sisotool(self.systems[0]) + sisotool(self.system) fig = plt.gcf() ax_mag,ax_rlocus,ax_phase,ax_step = fig.axes[0],fig.axes[1],fig.axes[2],fig.axes[3] @@ -53,7 +49,7 @@ def test_sisotool(self): # Move the rootlocus to another point event = type('test', (object,), {'xdata': 2.31206868287,'ydata':15.5983051046, 'inaxes':ax_rlocus.axes})() - _RLFeedbackClicksSisotool(event=event, sys=self.systems[0], fig=fig, bode_plot_params=bode_plot_params, tvect=None) + _RLFeedbackClicksSisotool(event=event, sys=self.system, fig=fig, bode_plot_params=bode_plot_params, tvect=None) # Check the moved root locus plot points moved_point_0 = (np.array([-29.91742755]), np.array([0.])) From c0e533e219f9ee3ddcfa7ffaaea61ad549bb4465 Mon Sep 17 00:00:00 2001 From: David de Jong Date: Sat, 2 Jun 2018 13:16:15 +0200 Subject: [PATCH 09/21] root locus smooth with zoom bug --- control/rlocus.py | 210 ++++++++++++++++++++++++++++++++++------------ 1 file changed, 155 insertions(+), 55 deletions(-) diff --git a/control/rlocus.py b/control/rlocus.py index 3762c12e2..41c17445f 100644 --- a/control/rlocus.py +++ b/control/rlocus.py @@ -48,6 +48,7 @@ # Packages used by this module import numpy as np import matplotlib +import matplotlib.pyplot as plt from scipy import array, poly1d, row_stack, zeros_like, real, imag import scipy.signal # signal processing toolbox import pylab # plotting routines @@ -58,9 +59,8 @@ __all__ = ['root_locus', 'rlocus'] - # Main function: compute a root locus diagram -def root_locus(sys, kvect=None, xlim=None, ylim=None, plotstr='-', Plot=True, +def root_locus(sys, kvect=None, xlim=None, ylim=None, plotstr='b' if int(matplotlib.__version__[0]) == 1 else 'C0', Plot=True, PrintGain=True, grid=False, **kwargs): """Root locus plot @@ -128,13 +128,13 @@ def root_locus(sys, kvect=None, xlim=None, ylim=None, plotstr='-', Plot=True, if PrintGain and sisotool == False: f.canvas.mpl_connect( - 'button_release_event', partial(_RLFeedbackClicksPoint,sys=sys,fig=f)) + 'button_release_event', partial(_RLClickDispatcher,sys=sys, fig=f,ax_rlocus=f.axes[0],plotstr=plotstr)) elif sisotool == True: f.axes[1].plot([root.real for root in start_mat], [root.imag for root in start_mat], 'm.', marker='s', markersize=8,zorder=20,label='gain_point') f.suptitle("Clicked at: %10.4g%+10.4gj gain: %10.4g damp: %10.4g" % (start_mat[0][0].real, start_mat[0][0].imag, 1, -1 * start_mat[0][0].real / abs(start_mat[0][0])),fontsize = 12 if int(matplotlib.__version__[0]) == 1 else 10) f.canvas.mpl_connect( - 'button_release_event',partial(_RLFeedbackClicksSisotool,sys=sys, fig=f, bode_plot_params=kwargs['bode_plot_params'],tvect=kwargs['tvect'])) + 'button_release_event',partial(_RLClickDispatcher,sys=sys, fig=f,ax_rlocus=f.axes[1],plotstr=plotstr, sisotool=sisotool, bode_plot_params=kwargs['bode_plot_params'],tvect=kwargs['tvect'])) # plot open loop poles poles = array(denp.r) @@ -146,8 +146,8 @@ def root_locus(sys, kvect=None, xlim=None, ylim=None, plotstr='-', Plot=True, ax.plot(real(zeros), imag(zeros), 'o') # Now plot the loci - for col in mymat.T: - ax.plot(real(col), imag(col), plotstr) + for index,col in enumerate(mymat.T): + ax.plot(real(col), imag(col), plotstr,label='rootlocus') # Set up plot axes and labels if xlim: @@ -165,8 +165,7 @@ def root_locus(sys, kvect=None, xlim=None, ylim=None, plotstr='-', Plot=True, ax.axvline(0., linestyle=':', color='k') return mymat, kvect - -def _default_gains(num, den, xlim, ylim): +def _default_gains(num, den, xlim, ylim, point_tolerance=5e-2,zoom_xlim=None,zoom_ylim=None): """Unsupervised gains calculation for root locus plot. References: @@ -198,7 +197,7 @@ def _default_gains(num, den, xlim, ylim): "locus with equal order of numerator and denominator.") if xlim is None and false_gain > 0: - x_tolerance = 0.05 * (np.max(np.real(mymat_xl)) - np.min(np.real(mymat_xl))) + x_tolerance = point_tolerance * (np.max(np.real(mymat_xl)) - np.min(np.real(mymat_xl))) xlim = _ax_lim(mymat_xl) elif xlim is None and false_gain < 0: axmin = np.min(np.real(important_points))-(np.max(np.real(important_points))-np.min(np.real(important_points))) @@ -206,39 +205,119 @@ def _default_gains(num, den, xlim, ylim): axmax = np.max(np.real(important_points))+np.max(np.real(important_points))-np.min(np.real(important_points)) axmax = np.max(np.array([axmax, np.max(np.real(mymat_xl))])) xlim = [axmin, axmax] - x_tolerance = 0.05 * (axmax - axmin) + x_tolerance = 5e-2 * (axmax - axmin) else: - x_tolerance = 0.05 * (xlim[1] - xlim[0]) + x_tolerance = 5e-2 * (xlim[1] - xlim[0]) if ylim is None: - y_tolerance = 0.05 * (np.max(np.imag(mymat_xl)) - np.min(np.imag(mymat_xl))) + y_tolerance = 5e-2 * (np.max(np.imag(mymat_xl)) - np.min(np.imag(mymat_xl))) ylim = _ax_lim(mymat_xl * 1j) else: - y_tolerance = 0.05 * (ylim[1] - ylim[0]) - + y_tolerance = 5e-2 * (ylim[1] - ylim[0]) + print(len(mymat)) tolerance = np.max([x_tolerance, y_tolerance]) - distance_points = np.abs(np.diff(mymat, axis=0)) - indexes_too_far = np.where(distance_points > tolerance) - - while (indexes_too_far[0].size > 0) and (kvect.size < 5000): - for index in indexes_too_far[0]: - new_gains = np.linspace(kvect[index], kvect[index+1], 5) - new_points = _RLFindRoots(num, den, new_gains[1:4]) - kvect = np.insert(kvect, index+1, new_gains[1:4]) - mymat = np.insert(mymat, index+1, new_points, axis=0) - + # print('normal smoothing start') + # mymat,kvect =_smooth_rootlocus(num, den, tolerance,mymat,kvect, xlim=None, ylim=None) + # print('normal smoothing end') + mymat = _RLSortRoots(mymat) + print(len(mymat)) + + # If a zoom on the plot is used insert points on this interval and use a smaller tolerance + if zoom_xlim != None and zoom_ylim != None: + print('zoom smoothing start') + y_tolerance = 5e-2 * (zoom_ylim[1] - zoom_ylim[0]) + x_tolerance = 5e-2 * (zoom_xlim[1] - zoom_xlim[0]) + tolerance = np.max([x_tolerance, y_tolerance]) + tolerance = np.max([x_tolerance, y_tolerance]) + #print(mymat) + mymat,kvect = _smooth_rootlocus(num,den,tolerance,mymat,kvect,zoom_xlim,zoom_ylim) mymat = _RLSortRoots(mymat) - distance_points = np.abs(np.diff(mymat, axis=0)) > tolerance # distance between points - indexes_too_far = np.where(distance_points) + #print(mymat) + print('zoom smoothing end') + print(len(mymat)) + kvect = np.sort(kvect) - new_gains = kvect[-1] * np.hstack((np.logspace(0, 3, 4))) - new_points = _RLFindRoots(num, den, new_gains[1:4]) - kvect = np.append(kvect, new_gains[1:4]) - mymat = np.concatenate((mymat, new_points), axis=0) - mymat = _RLSortRoots(mymat) + mymat = _RLFindRoots(num, den, kvect) return kvect, mymat, xlim, ylim +def _smooth_rootlocus(num,den,tolerance,mymat,kvect, xlim,ylim): + """Smooth the rootlocus plot by inserting points at locations where the distance between + two points exceeds a certain tolerance.""" + + distance_points = np.abs(np.diff(mymat, axis=0)) + indexes_too_far = np.where(distance_points > tolerance) + indexes_too_far_filtered = _indexes_filter(indexes_too_far,mymat,xlim,ylim) + + print('length of indexes too_far_filtered') + print(len(indexes_too_far_filtered)) + print(indexes_too_far_filtered) + + if len(indexes_too_far_filtered) > 0: + print('points are added') + # if indexes_too_far_filtered[0] != 0: + # point_before_start = indexes_too_far_filtered[0] -1 + # indexes_too_far_filtered.insert(0,point_before_start) + # + # if indexes_too_far_filtered[-1] != len(mymat): + # point_at_end = indexes_too_far_filtered[-1]+1 + # indexes_too_far_filtered.append(point_at_end) + + print('before while loop') + print(len(indexes_too_far_filtered)) + print(indexes_too_far_filtered) + while (len(indexes_too_far_filtered) > 0) and (kvect.size < 5000): + counter = 0 + for list_index, index in enumerate(indexes_too_far_filtered): + index+= counter*3 + new_gains = np.linspace(kvect[index], kvect[index+1], 5) + new_points = _RLFindRoots(num, den, new_gains[1:4]) + + print('inside while loop') + print(index) + print(kvect[index],kvect[index+1]) + print(new_gains[1:4]) + + kvect = np.insert(kvect, index+1, new_gains[1:4]) + + mymat = np.insert(mymat, index+1, new_points, axis=0) + counter +=1 + + + + mymat = _RLSortRoots(mymat) + distance_points = np.abs(np.diff(mymat, axis=0)) > tolerance + indexes_too_far = np.where(distance_points) + indexes_too_far_filtered = _indexes_filter(indexes_too_far, mymat, xlim, ylim) + print('new indexes too far') + print(indexes_too_far_filtered) + print('K SORTED?') + print(all(kvect[i] <= kvect[i+1] for i in range(len(kvect)-1))) + + #print(kvect) + + new_gains = kvect[-1] * np.hstack((np.logspace(0, 3, 4))) + new_points = _RLFindRoots(num, den, new_gains[1:4]) + kvect = np.append(kvect, new_gains[1:4]) + mymat = np.concatenate((mymat, new_points), axis=0) + + + return mymat,kvect + +def _indexes_filter(indexes_too_far,mymat,xlim,ylim): + """Filter the indexes so only the resolution of points within the xlim and ylim is improved""" + if xlim== None and ylim == None: + indexes_too_far_filtered = list(np.unique(indexes_too_far[0])) + else: + indexes_too_far_filtered = [] + for index in np.unique(indexes_too_far[0]): + for point in mymat[index]: + if (xlim[0] <= point.real <= xlim[1]) and (ylim[0] <= point.imag <=ylim[1]): + indexes_too_far_filtered.append(index) + break + + return indexes_too_far_filtered + def _break_points(num, den): """Extract break points over real axis and the gains give these location""" # type: (np.poly1d, np.poly1d) -> (np.array, np.array) @@ -365,27 +444,46 @@ def _RLSortRoots(mymat): prevrow = sorted[n, :] return sorted -def _RLFeedbackClicksSisotool(event,sys,fig,bode_plot_params,tvect): - """Update Sisotool plots if a new point on the root locus plot is clicked - """ +def _RLClickDispatcher(event,sys,fig,ax_rlocus,plotstr,sisotool=False,bode_plot_params=None,tvect=None): + """Rootlocus plot click dispatcher""" + + # If zoom is used on the rootlocus plot smooth and update it + if plt.get_current_fig_manager().toolbar.mode == 'zoom rect' and event.inaxes == ax_rlocus.axes: + + (nump, denp) = _systopoly1d(sys) + xlim = ax_rlocus.get_xlim() + ylim = ax_rlocus.get_ylim() + kvect,mymat, xlim,ylim = _default_gains(nump, denp,xlim=None,ylim=None, zoom_xlim=xlim,zoom_ylim=ylim) + _removeLine('rootlocus', ax_rlocus) + + for index,col in enumerate(mymat.T): + ax_rlocus.plot(real(col), imag(col), plotstr,label=index) + + fig.canvas.draw() + + # if a point is clicked on the rootlocus plot visually emphasize it + else: + K = _RLFeedbackClicksPoint(event, sys, fig,ax_rlocus,sisotool) + if sisotool and K is not None: + _SisotoolUpdate(sys, fig, K, bode_plot_params, tvect) + + # Update the canvas + fig.canvas.draw() + + +def _RLFeedbackClicksPoint(event,sys,fig,ax_rlocus,sisotool=False): + """Display root-locus gain feedback point for clicks on the root-locus plot""" + + (nump, denp) = _systopoly1d(sys) + + # Catch type error when event click is in the figure but not in an axis try: s = complex(event.xdata, event.ydata) K = -1. / sys.horner(s) - ax_rlocus = fig.axes[1] - if _RLFeedbackClicksPoint(event, sys, fig, ax_rlocus, sisotool=True): - _SisotoolUpdate(sys, fig, K.real[0][0], bode_plot_params, tvect) - except TypeError: - pass -def _RLFeedbackClicksPoint(event,sys,fig,ax_rlocus=None,sisotool=False): - """Display root-locus gain feedback point for clicks on the root-locus plot - """ - if sisotool == False: - ax_rlocus = fig.axes[0] + except TypeError: + K = float('inf') - (nump, denp) = _systopoly1d(sys) - s = complex(event.xdata, event.ydata) - K = -1. / sys.horner(s) if abs(K.real) > 1e-8 and abs(K.imag / K.real) < 0.04 and event.inaxes == ax_rlocus.axes: # Display the parameters in the output window and figure @@ -394,11 +492,8 @@ def _RLFeedbackClicksPoint(event,sys,fig,ax_rlocus=None,sisotool=False): fig.suptitle("Clicked at: %10.4g%+10.4gj gain: %10.4g damp: %10.4g" % (s.real, s.imag, K.real, -1 * s.real / abs(s)),fontsize = 12 if int(matplotlib.__version__[0]) == 1 else 10) - # Remove the previous points - for line in reversed(ax_rlocus.lines): - if len(line.get_xdata()) == 1 and line.get_label()=='gain_point': - line.remove() - del line + # Remove the previous line + _removeLine(label='gain_point',ax=ax_rlocus) # Visualise clicked point, display all roots for sisotool mode if sisotool: @@ -408,10 +503,15 @@ def _RLFeedbackClicksPoint(event,sys,fig,ax_rlocus=None,sisotool=False): else: ax_rlocus.plot(s.real, s.imag, 'k.', marker='s', markersize=8, zorder=20,label='gain_point') - # Update the canvas - fig.canvas.draw() + return K.real[0][0] + - return True +def _removeLine(label,ax): + """Remove a line from the ax when a label is specified""" + for line in reversed(ax.lines): + if line.get_label() == label: + line.remove() + del line def _sgrid_func(fig=None, zeta=None, wn=None): if fig is None: From fe13ae19fdf9e25679038dda6f8927daf48982e3 Mon Sep 17 00:00:00 2001 From: David de Jong Date: Sat, 2 Jun 2018 18:56:51 +0200 Subject: [PATCH 10/21] Fixed zooming on plot. Now displays smoother root locus plot. --- control/rlocus.py | 182 ++++++++++++++++++---------------------------- 1 file changed, 72 insertions(+), 110 deletions(-) diff --git a/control/rlocus.py b/control/rlocus.py index 41c17445f..65b95641f 100644 --- a/control/rlocus.py +++ b/control/rlocus.py @@ -165,9 +165,10 @@ def root_locus(sys, kvect=None, xlim=None, ylim=None, plotstr='b' if int(matplot ax.axvline(0., linestyle=':', color='k') return mymat, kvect -def _default_gains(num, den, xlim, ylim, point_tolerance=5e-2,zoom_xlim=None,zoom_ylim=None): - """Unsupervised gains calculation for root locus plot. - + +def _default_gains(num, den, xlim, ylim,zoom_xlim=None,zoom_ylim=None): + """Unsupervised gains calculation for root locus plot. + References: Ogata, K. (2002). Modern control engineering (4th ed.). Upper Saddle River, NJ : New Delhi: Prentice Hall..""" @@ -197,126 +198,92 @@ def _default_gains(num, den, xlim, ylim, point_tolerance=5e-2,zoom_xlim=None,zoo "locus with equal order of numerator and denominator.") if xlim is None and false_gain > 0: - x_tolerance = point_tolerance * (np.max(np.real(mymat_xl)) - np.min(np.real(mymat_xl))) + x_tolerance = 0.05 * (np.max(np.real(mymat_xl)) - np.min(np.real(mymat_xl))) xlim = _ax_lim(mymat_xl) elif xlim is None and false_gain < 0: - axmin = np.min(np.real(important_points))-(np.max(np.real(important_points))-np.min(np.real(important_points))) + axmin = np.min(np.real(important_points)) - ( + np.max(np.real(important_points)) - np.min(np.real(important_points))) axmin = np.min(np.array([axmin, np.min(np.real(mymat_xl))])) - axmax = np.max(np.real(important_points))+np.max(np.real(important_points))-np.min(np.real(important_points)) + axmax = np.max(np.real(important_points)) + np.max(np.real(important_points)) - np.min( + np.real(important_points)) axmax = np.max(np.array([axmax, np.max(np.real(mymat_xl))])) xlim = [axmin, axmax] - x_tolerance = 5e-2 * (axmax - axmin) + x_tolerance = 0.05 * (axmax - axmin) else: - x_tolerance = 5e-2 * (xlim[1] - xlim[0]) + x_tolerance = 0.05 * (xlim[1] - xlim[0]) if ylim is None: - y_tolerance = 5e-2 * (np.max(np.imag(mymat_xl)) - np.min(np.imag(mymat_xl))) + y_tolerance = 0.05 * (np.max(np.imag(mymat_xl)) - np.min(np.imag(mymat_xl))) ylim = _ax_lim(mymat_xl * 1j) else: - y_tolerance = 5e-2 * (ylim[1] - ylim[0]) - print(len(mymat)) - tolerance = np.max([x_tolerance, y_tolerance]) - # print('normal smoothing start') - # mymat,kvect =_smooth_rootlocus(num, den, tolerance,mymat,kvect, xlim=None, ylim=None) - # print('normal smoothing end') - mymat = _RLSortRoots(mymat) - print(len(mymat)) + y_tolerance = 0.05 * (ylim[1] - ylim[0]) + + tolerance = np.min([x_tolerance, y_tolerance]) + indexes_too_far = _indexes_filt(mymat,tolerance,zoom_xlim,zoom_ylim) + + while (len(indexes_too_far) > 0) and (kvect.size < 5000): + for counter,index in enumerate(indexes_too_far): + index = index + counter*3 + new_gains = np.linspace(kvect[index], kvect[index + 1], 5) + new_points = _RLFindRoots(num, den, new_gains[1:4]) + kvect = np.insert(kvect, index + 1, new_gains[1:4]) + mymat = np.insert(mymat, index + 1, new_points, axis=0) - # If a zoom on the plot is used insert points on this interval and use a smaller tolerance - if zoom_xlim != None and zoom_ylim != None: - print('zoom smoothing start') - y_tolerance = 5e-2 * (zoom_ylim[1] - zoom_ylim[0]) - x_tolerance = 5e-2 * (zoom_xlim[1] - zoom_xlim[0]) - tolerance = np.max([x_tolerance, y_tolerance]) - tolerance = np.max([x_tolerance, y_tolerance]) - #print(mymat) - mymat,kvect = _smooth_rootlocus(num,den,tolerance,mymat,kvect,zoom_xlim,zoom_ylim) mymat = _RLSortRoots(mymat) - #print(mymat) - print('zoom smoothing end') - print(len(mymat)) - kvect = np.sort(kvect) + indexes_too_far = _indexes_filt(mymat,tolerance,zoom_xlim,zoom_ylim) - mymat = _RLFindRoots(num, den, kvect) + new_gains = kvect[-1] * np.hstack((np.logspace(0, 3, 4))) + new_points = _RLFindRoots(num, den, new_gains[1:4]) + kvect = np.append(kvect, new_gains[1:4]) + mymat = np.concatenate((mymat, new_points), axis=0) + mymat = _RLSortRoots(mymat) return kvect, mymat, xlim, ylim - -def _smooth_rootlocus(num,den,tolerance,mymat,kvect, xlim,ylim): - """Smooth the rootlocus plot by inserting points at locations where the distance between - two points exceeds a certain tolerance.""" +def _indexes_filt(mymat,tolerance,zoom_xlim=None,zoom_ylim=None): + """Calculate the distance between points and return the indexes. + Filter the indexes so only the resolution of points within the xlim and ylim is improved when zoom is used""" distance_points = np.abs(np.diff(mymat, axis=0)) - indexes_too_far = np.where(distance_points > tolerance) - indexes_too_far_filtered = _indexes_filter(indexes_too_far,mymat,xlim,ylim) - - print('length of indexes too_far_filtered') - print(len(indexes_too_far_filtered)) - print(indexes_too_far_filtered) - - if len(indexes_too_far_filtered) > 0: - print('points are added') - # if indexes_too_far_filtered[0] != 0: - # point_before_start = indexes_too_far_filtered[0] -1 - # indexes_too_far_filtered.insert(0,point_before_start) - # - # if indexes_too_far_filtered[-1] != len(mymat): - # point_at_end = indexes_too_far_filtered[-1]+1 - # indexes_too_far_filtered.append(point_at_end) - - print('before while loop') - print(len(indexes_too_far_filtered)) - print(indexes_too_far_filtered) - while (len(indexes_too_far_filtered) > 0) and (kvect.size < 5000): - counter = 0 - for list_index, index in enumerate(indexes_too_far_filtered): - index+= counter*3 - new_gains = np.linspace(kvect[index], kvect[index+1], 5) - new_points = _RLFindRoots(num, den, new_gains[1:4]) - - print('inside while loop') - print(index) - print(kvect[index],kvect[index+1]) - print(new_gains[1:4]) - - kvect = np.insert(kvect, index+1, new_gains[1:4]) - - mymat = np.insert(mymat, index+1, new_points, axis=0) - counter +=1 - - - - mymat = _RLSortRoots(mymat) - distance_points = np.abs(np.diff(mymat, axis=0)) > tolerance - indexes_too_far = np.where(distance_points) - indexes_too_far_filtered = _indexes_filter(indexes_too_far, mymat, xlim, ylim) - print('new indexes too far') - print(indexes_too_far_filtered) - print('K SORTED?') - print(all(kvect[i] <= kvect[i+1] for i in range(len(kvect)-1))) - - #print(kvect) - - new_gains = kvect[-1] * np.hstack((np.logspace(0, 3, 4))) - new_points = _RLFindRoots(num, den, new_gains[1:4]) - kvect = np.append(kvect, new_gains[1:4]) - mymat = np.concatenate((mymat, new_points), axis=0) - - - return mymat,kvect - -def _indexes_filter(indexes_too_far,mymat,xlim,ylim): - """Filter the indexes so only the resolution of points within the xlim and ylim is improved""" - if xlim== None and ylim == None: - indexes_too_far_filtered = list(np.unique(indexes_too_far[0])) - else: + indexes_too_far = list(np.unique(np.where(distance_points > tolerance)[0])) + + if zoom_xlim != None and zoom_ylim != None: + x_tolerance_zoom = 0.05 * (zoom_xlim[1] - zoom_xlim[0]) + y_tolerance_zoom = 0.05 * (zoom_ylim[1] - zoom_ylim[0]) + tolerance_zoom = np.min([x_tolerance_zoom, y_tolerance_zoom]) + distance_points_zoom_ = np.abs(np.diff(mymat, axis=0)) + indexes_too_far_zoom = list(np.unique(np.where(distance_points_zoom_ > tolerance_zoom)[0])) indexes_too_far_filtered = [] - for index in np.unique(indexes_too_far[0]): + + for index in indexes_too_far_zoom: for point in mymat[index]: - if (xlim[0] <= point.real <= xlim[1]) and (ylim[0] <= point.imag <=ylim[1]): + if (zoom_xlim[0] <= point.real <= zoom_xlim[1]) and (zoom_ylim[0] <= point.imag <= zoom_ylim[1]): indexes_too_far_filtered.append(index) break - return indexes_too_far_filtered + # Check if the zoom box is not overshot and insert points where neccessary + if len(indexes_too_far_filtered) == 0 and len(mymat) <300: + limits = [zoom_xlim[0],zoom_xlim[1],zoom_ylim[0],zoom_ylim[1]] + for index,limit in enumerate(limits): + if index <= 1: + asign = np.sign(real(mymat)-limit) + else: + asign = np.sign(imag(mymat) - limit) + signchange = ((np.roll(asign, 1, axis=0) - asign) != 0).astype(int) + signchange[0] = np.zeros((len(mymat[0]))) + if len(np.where(signchange ==1)) > 0: + indexes_too_far_filtered.append(np.where(signchange == 1)[0][0]) + + if len(indexes_too_far_filtered) > 0 : + if indexes_too_far_filtered[0] != 0: + indexes_too_far_filtered.insert(0,indexes_too_far_filtered[0]-1) + if not indexes_too_far_filtered[-1] +1 >= len(mymat)-2: + indexes_too_far_filtered.append(indexes_too_far_filtered[-1]+1) + + indexes_too_far.extend(indexes_too_far_filtered) + + indexes_too_far = list(np.unique(indexes_too_far)) + indexes_too_far.sort() + return indexes_too_far def _break_points(num, den): """Extract break points over real axis and the gains give these location""" @@ -405,7 +372,6 @@ def _systopoly1d(sys): def _RLFindRoots(nump, denp, kvect): """Find the roots for the root locus.""" - # Convert numerator and denominator to polynomials if they aren't roots = [] for k in kvect: @@ -421,7 +387,6 @@ def _RLFindRoots(nump, denp, kvect): mymat = row_stack(roots) return mymat - def _RLSortRoots(mymat): """Sort the roots from sys._RLFindRoots, so that the root locus doesn't show weird pseudo-branches as roots jump from @@ -446,20 +411,17 @@ def _RLSortRoots(mymat): def _RLClickDispatcher(event,sys,fig,ax_rlocus,plotstr,sisotool=False,bode_plot_params=None,tvect=None): """Rootlocus plot click dispatcher""" - # If zoom is used on the rootlocus plot smooth and update it if plt.get_current_fig_manager().toolbar.mode == 'zoom rect' and event.inaxes == ax_rlocus.axes: (nump, denp) = _systopoly1d(sys) - xlim = ax_rlocus.get_xlim() - ylim = ax_rlocus.get_ylim() + xlim,ylim = ax_rlocus.get_xlim(),ax_rlocus.get_ylim() + kvect,mymat, xlim,ylim = _default_gains(nump, denp,xlim=None,ylim=None, zoom_xlim=xlim,zoom_ylim=ylim) _removeLine('rootlocus', ax_rlocus) - for index,col in enumerate(mymat.T): - ax_rlocus.plot(real(col), imag(col), plotstr,label=index) - - fig.canvas.draw() + for i,col in enumerate(mymat.T): + ax_rlocus.plot(real(col), imag(col), plotstr,label='rootlocus') # if a point is clicked on the rootlocus plot visually emphasize it else: From 2ba827764aa14534bc2a698e1473b63652115f35 Mon Sep 17 00:00:00 2001 From: David de Jong Date: Sun, 3 Jun 2018 14:07:19 +0200 Subject: [PATCH 11/21] Fixed paning on plot and clicking on the root locus plot while zoomed --- control/rlocus.py | 20 +++++++++++++------- 1 file changed, 13 insertions(+), 7 deletions(-) diff --git a/control/rlocus.py b/control/rlocus.py index 65b95641f..2c68f716e 100644 --- a/control/rlocus.py +++ b/control/rlocus.py @@ -174,6 +174,7 @@ def _default_gains(num, den, xlim, ylim,zoom_xlim=None,zoom_ylim=None): k_break, real_break = _break_points(num, den) kmax = _k_max(num, den, real_break, k_break) + print(kmax) kvect = np.hstack((np.linspace(0, kmax, 50), np.real(k_break))) kvect.sort() @@ -242,7 +243,6 @@ def _default_gains(num, den, xlim, ylim,zoom_xlim=None,zoom_ylim=None): def _indexes_filt(mymat,tolerance,zoom_xlim=None,zoom_ylim=None): """Calculate the distance between points and return the indexes. Filter the indexes so only the resolution of points within the xlim and ylim is improved when zoom is used""" - distance_points = np.abs(np.diff(mymat, axis=0)) indexes_too_far = list(np.unique(np.where(distance_points > tolerance)[0])) @@ -250,8 +250,7 @@ def _indexes_filt(mymat,tolerance,zoom_xlim=None,zoom_ylim=None): x_tolerance_zoom = 0.05 * (zoom_xlim[1] - zoom_xlim[0]) y_tolerance_zoom = 0.05 * (zoom_ylim[1] - zoom_ylim[0]) tolerance_zoom = np.min([x_tolerance_zoom, y_tolerance_zoom]) - distance_points_zoom_ = np.abs(np.diff(mymat, axis=0)) - indexes_too_far_zoom = list(np.unique(np.where(distance_points_zoom_ > tolerance_zoom)[0])) + indexes_too_far_zoom = list(np.unique(np.where(distance_points > tolerance_zoom)[0])) indexes_too_far_filtered = [] for index in indexes_too_far_zoom: @@ -270,8 +269,8 @@ def _indexes_filt(mymat,tolerance,zoom_xlim=None,zoom_ylim=None): asign = np.sign(imag(mymat) - limit) signchange = ((np.roll(asign, 1, axis=0) - asign) != 0).astype(int) signchange[0] = np.zeros((len(mymat[0]))) - if len(np.where(signchange ==1)) > 0: - indexes_too_far_filtered.append(np.where(signchange == 1)[0][0]) + if len(np.where(signchange ==1)[0]) > 0: + indexes_too_far_filtered.append(np.where(signchange == 1)[0][0]-1) if len(indexes_too_far_filtered) > 0 : if indexes_too_far_filtered[0] != 0: @@ -411,8 +410,9 @@ def _RLSortRoots(mymat): def _RLClickDispatcher(event,sys,fig,ax_rlocus,plotstr,sisotool=False,bode_plot_params=None,tvect=None): """Rootlocus plot click dispatcher""" + # If zoom is used on the rootlocus plot smooth and update it - if plt.get_current_fig_manager().toolbar.mode == 'zoom rect' and event.inaxes == ax_rlocus.axes: + if plt.get_current_fig_manager().toolbar.mode in ['zoom rect','pan/zoom'] and event.inaxes == ax_rlocus.axes: (nump, denp) = _systopoly1d(sys) xlim,ylim = ax_rlocus.get_xlim(),ax_rlocus.get_ylim() @@ -446,7 +446,13 @@ def _RLFeedbackClicksPoint(event,sys,fig,ax_rlocus,sisotool=False): except TypeError: K = float('inf') - if abs(K.real) > 1e-8 and abs(K.imag / K.real) < 0.04 and event.inaxes == ax_rlocus.axes: + xlim = ax_rlocus.get_xlim() + ylim = ax_rlocus.get_ylim() + x_tolerance = 0.05 * (xlim[1] - xlim[0]) + y_tolerance = 0.05 * (ylim[1] - ylim[0]) + gain_tolerance = np.min([x_tolerance, y_tolerance])*1e-1 + + if abs(K.real) > 1e-8 and abs(K.imag / K.real) < gain_tolerance and event.inaxes == ax_rlocus.axes: # Display the parameters in the output window and figure print("Clicked at %10.4g%+10.4gj gain %10.4g damp %10.4g" % From c26187f9d46f1d80353b16d74bd033be260aabf7 Mon Sep 17 00:00:00 2001 From: David de Jong Date: Sun, 3 Jun 2018 15:15:12 +0200 Subject: [PATCH 12/21] Added unit test for rootlocus zoom functionality --- control/rlocus.py | 1 - control/tests/rlocus_test.py | 30 +++++++++++++++++++++++++++++- 2 files changed, 29 insertions(+), 2 deletions(-) diff --git a/control/rlocus.py b/control/rlocus.py index 2c68f716e..e36a1afe2 100644 --- a/control/rlocus.py +++ b/control/rlocus.py @@ -174,7 +174,6 @@ def _default_gains(num, den, xlim, ylim,zoom_xlim=None,zoom_ylim=None): k_break, real_break = _break_points(num, den) kmax = _k_max(num, den, real_break, k_break) - print(kmax) kvect = np.hstack((np.linspace(0, kmax, 50), np.real(k_break))) kvect.sort() diff --git a/control/tests/rlocus_test.py b/control/tests/rlocus_test.py index d2522c881..6b143ae7c 100644 --- a/control/tests/rlocus_test.py +++ b/control/tests/rlocus_test.py @@ -5,10 +5,13 @@ import unittest import numpy as np -from control.rlocus import root_locus +from control.rlocus import root_locus, _RLClickDispatcher from control.xferfcn import TransferFunction from control.statesp import StateSpace from control.bdalg import feedback +import matplotlib.pyplot as plt +from control.tests.margin_test import assert_array_almost_equal + class TestRootLocus(unittest.TestCase): """These are tests for the feedback function in rlocus.py.""" @@ -42,6 +45,31 @@ def test_without_gains(self): roots, kvect = root_locus(sys, Plot=False) self.check_cl_poles(sys, roots, kvect) + def testRootLocusZoom(self): + """Check the zooming functionality of the Root locus plot""" + system = TransferFunction([1000], [1, 25, 100, 0]) + root_locus(system) + fig = plt.gcf() + ax_rlocus = fig.axes[0] + + event = type('test', (object,), {'xdata': 0.9957380594313321, 'ydata': 1.7825491928580846, 'inaxes': ax_rlocus.axes})() + ax_rlocus.set_xlim((-4.420022219849855, 0.9957380594313321)) + ax_rlocus.set_ylim((1.7825491928580846, 10.695295157148486)) + plt.get_current_fig_manager().toolbar.mode = 'zoom rect' + _RLClickDispatcher(event, system, fig, ax_rlocus, '-') + + zoom_x = ax_rlocus.lines[-2].get_data()[0][65:75] + zoom_y = ax_rlocus.lines[-2].get_data()[1][65:75] + zoom_y = [abs(y) for y in zoom_y] + + zoom_x_valid = [-2.23659192, -2.23659121, -2.23659103, -2.23659086, -2.23659068, -2.23659063, + -2.23659059, -2.23659055, -2.2365905, -2.23658766] + zoom_y_valid = [1.78253589, 1.78254318, 1.782545, 1.78254682, 1.78254864, 1.7825491, + 1.78254955, 1.78255001, 1.78255046, 1.7825796] + + assert_array_almost_equal(zoom_x,zoom_x_valid) + assert_array_almost_equal(zoom_y,zoom_y_valid) + def test_suite(): return unittest.TestLoader().loadTestsFromTestCase(TestRootLocus) From 9962ceafdcf592f1147287f7dd68a5e15eb9dac0 Mon Sep 17 00:00:00 2001 From: David de Jong Date: Sun, 3 Jun 2018 15:29:06 +0200 Subject: [PATCH 13/21] Fixed sisotool unittest --- control/tests/sisotool_test.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/control/tests/sisotool_test.py b/control/tests/sisotool_test.py index 8db6fe9a6..84d0fe8d6 100644 --- a/control/tests/sisotool_test.py +++ b/control/tests/sisotool_test.py @@ -2,7 +2,7 @@ import numpy as np from control.sisotool import sisotool from control.tests.margin_test import assert_array_almost_equal -from control.rlocus import _RLFeedbackClicksSisotool +from control.rlocus import _RLClickDispatcher from control.xferfcn import TransferFunction import matplotlib.pyplot as plt @@ -49,12 +49,14 @@ def test_sisotool(self): # Move the rootlocus to another point event = type('test', (object,), {'xdata': 2.31206868287,'ydata':15.5983051046, 'inaxes':ax_rlocus.axes})() - _RLFeedbackClicksSisotool(event=event, sys=self.system, fig=fig, bode_plot_params=bode_plot_params, tvect=None) + _RLClickDispatcher(event=event, sys=self.system, fig=fig,ax_rlocus=ax_rlocus,sisotool=True, plotstr='-' ,bode_plot_params=bode_plot_params, tvect=None) # Check the moved root locus plot points moved_point_0 = (np.array([-29.91742755]), np.array([0.])) moved_point_1 = (np.array([2.45871378]), np.array([-15.52647768])) moved_point_2 = (np.array([2.45871378]), np.array([15.52647768])) + for line in ax_rlocus.lines: + print(line.get_data()) assert_array_almost_equal(ax_rlocus.lines[-3].get_data(),moved_point_0) assert_array_almost_equal(ax_rlocus.lines[-2].get_data(),moved_point_1) assert_array_almost_equal(ax_rlocus.lines[-1].get_data(),moved_point_2) From 938ae886c8d3fbcdf4e8721c6a81944488cad1b0 Mon Sep 17 00:00:00 2001 From: David de Jong Date: Sun, 3 Jun 2018 15:41:57 +0200 Subject: [PATCH 14/21] Fixed frequency response unittest --- control/tests/freqresp_test.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/control/tests/freqresp_test.py b/control/tests/freqresp_test.py index 781c40496..149db9fa9 100644 --- a/control/tests/freqresp_test.py +++ b/control/tests/freqresp_test.py @@ -112,7 +112,7 @@ def test_bode_margin(self): num = [1000] den = [1, 25, 100, 0] sys = ctrl.tf(num, den) - ctrl.bode_plot(sys, margins=True,dB=False,deg = True) + ctrl.bode_plot(sys, margins=True,dB=False,deg = True, Hz=False) fig = plt.gcf() allaxes = fig.get_axes() From 9f13a778cf81c1a42c20b58d609921f5c90dd902 Mon Sep 17 00:00:00 2001 From: David de Jong Date: Sun, 3 Jun 2018 15:54:57 +0200 Subject: [PATCH 15/21] small change to sisotool test --- control/tests/sisotool_test.py | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/control/tests/sisotool_test.py b/control/tests/sisotool_test.py index 84d0fe8d6..68a2173ed 100644 --- a/control/tests/sisotool_test.py +++ b/control/tests/sisotool_test.py @@ -14,7 +14,7 @@ def setUp(self): self.system = TransferFunction([1000],[1,25,100,0]) def test_sisotool(self): - sisotool(self.system) + sisotool(self.system,Hz=False) fig = plt.gcf() ax_mag,ax_rlocus,ax_phase,ax_step = fig.axes[0],fig.axes[1],fig.axes[2],fig.axes[3] @@ -31,7 +31,6 @@ def test_sisotool(self): assert_array_almost_equal(ax_mag.lines[0].get_data()[1][10:20],bode_mag_original) # Check the step response before moving the point - print(ax_step.lines[0].get_data()[1][:10]) step_response_original = np.array([ 0., 0.02233651, 0.13118374, 0.33078542, 0.5907113, 0.87041549, 1.13038536, 1.33851053, 1.47374666, 1.52757114]) assert_array_almost_equal(ax_step.lines[0].get_data()[1][:10],step_response_original) @@ -55,8 +54,6 @@ def test_sisotool(self): moved_point_0 = (np.array([-29.91742755]), np.array([0.])) moved_point_1 = (np.array([2.45871378]), np.array([-15.52647768])) moved_point_2 = (np.array([2.45871378]), np.array([15.52647768])) - for line in ax_rlocus.lines: - print(line.get_data()) assert_array_almost_equal(ax_rlocus.lines[-3].get_data(),moved_point_0) assert_array_almost_equal(ax_rlocus.lines[-2].get_data(),moved_point_1) assert_array_almost_equal(ax_rlocus.lines[-1].get_data(),moved_point_2) From 6e06a9cafee486108a167e0b67452ede05ea5d2d Mon Sep 17 00:00:00 2001 From: David de Jong Date: Sun, 3 Jun 2018 16:23:21 +0200 Subject: [PATCH 16/21] small change to unittests --- control/rlocus.py | 2 +- control/tests/rlocus_test.py | 2 +- control/tests/sisotool_test.py | 4 ---- 3 files changed, 2 insertions(+), 6 deletions(-) diff --git a/control/rlocus.py b/control/rlocus.py index e36a1afe2..d551a42f7 100644 --- a/control/rlocus.py +++ b/control/rlocus.py @@ -102,7 +102,7 @@ def root_locus(sys, kvect=None, xlim=None, ylim=None, plotstr='b' if int(matplot start_mat = _RLFindRoots(nump, denp, [1]) kvect, mymat, xlim, ylim = _default_gains(nump, denp, xlim, ylim) else: - start_mat = _RLFindRoots(nump, denp, kvect[0]) + start_mat = _RLFindRoots(nump, denp, [kvect[0]]) mymat = _RLFindRoots(nump, denp, kvect) mymat = _RLSortRoots(mymat) diff --git a/control/tests/rlocus_test.py b/control/tests/rlocus_test.py index 6b143ae7c..56ac39206 100644 --- a/control/tests/rlocus_test.py +++ b/control/tests/rlocus_test.py @@ -45,7 +45,7 @@ def test_without_gains(self): roots, kvect = root_locus(sys, Plot=False) self.check_cl_poles(sys, roots, kvect) - def testRootLocusZoom(self): + def test_root_locus_zoom(self): """Check the zooming functionality of the Root locus plot""" system = TransferFunction([1000], [1, 25, 100, 0]) root_locus(system) diff --git a/control/tests/sisotool_test.py b/control/tests/sisotool_test.py index 68a2173ed..40ef0f966 100644 --- a/control/tests/sisotool_test.py +++ b/control/tests/sisotool_test.py @@ -26,10 +26,6 @@ def test_sisotool(self): assert_array_almost_equal(ax_rlocus.lines[1].get_data(),initial_point_1) assert_array_almost_equal(ax_rlocus.lines[2].get_data(),initial_point_2) - # Check the bode plot magnitude line - bode_mag_original = np.array([ 15.12670407, 12.48358602, 10.28367521, 8.44961576, 6.91740918, 5.63436301, 4.55750665, 3.65237575, 2.89200358, 2.25587634]) - assert_array_almost_equal(ax_mag.lines[0].get_data()[1][10:20],bode_mag_original) - # Check the step response before moving the point step_response_original = np.array([ 0., 0.02233651, 0.13118374, 0.33078542, 0.5907113, 0.87041549, 1.13038536, 1.33851053, 1.47374666, 1.52757114]) assert_array_almost_equal(ax_step.lines[0].get_data()[1][:10],step_response_original) From cdb9f42d0e28a3590d1a0e917ebb6a052f2028d9 Mon Sep 17 00:00:00 2001 From: David de Jong Date: Sun, 3 Jun 2018 17:21:35 +0200 Subject: [PATCH 17/21] fixed python 2 integer division error in false gain of rootlocus --- control/rlocus.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/control/rlocus.py b/control/rlocus.py index d551a42f7..28c381bc1 100644 --- a/control/rlocus.py +++ b/control/rlocus.py @@ -192,7 +192,8 @@ def _default_gains(num, den, xlim, ylim,zoom_xlim=None,zoom_ylim=None): important_points = np.concatenate((singular_points, real_break), axis=0) important_points = np.concatenate((important_points, np.zeros(2)), axis=0) mymat_xl = np.append(mymat_xl, important_points) - false_gain = den.coeffs[0] / num.coeffs[0] + + false_gain = float(den.coeffs[0]) / float(num.coeffs[0]) if false_gain < 0 and not den.order > num.order: raise ValueError("Not implemented support for 0 degrees root " "locus with equal order of numerator and denominator.") From 90959bdc57dc2cc43bf5dc27c87ccae79d766815 Mon Sep 17 00:00:00 2001 From: David de Jong Date: Sun, 3 Jun 2018 17:52:55 +0200 Subject: [PATCH 18/21] changed root locus test to a numerically more significant one. --- control/rlocus.py | 1 - control/tests/rlocus_test.py | 14 +++++++------- 2 files changed, 7 insertions(+), 8 deletions(-) diff --git a/control/rlocus.py b/control/rlocus.py index 28c381bc1..0dcc4cbee 100644 --- a/control/rlocus.py +++ b/control/rlocus.py @@ -413,7 +413,6 @@ def _RLClickDispatcher(event,sys,fig,ax_rlocus,plotstr,sisotool=False,bode_plot_ # If zoom is used on the rootlocus plot smooth and update it if plt.get_current_fig_manager().toolbar.mode in ['zoom rect','pan/zoom'] and event.inaxes == ax_rlocus.axes: - (nump, denp) = _systopoly1d(sys) xlim,ylim = ax_rlocus.get_xlim(),ax_rlocus.get_ylim() diff --git a/control/tests/rlocus_test.py b/control/tests/rlocus_test.py index 56ac39206..58550b1d2 100644 --- a/control/tests/rlocus_test.py +++ b/control/tests/rlocus_test.py @@ -52,9 +52,9 @@ def test_root_locus_zoom(self): fig = plt.gcf() ax_rlocus = fig.axes[0] - event = type('test', (object,), {'xdata': 0.9957380594313321, 'ydata': 1.7825491928580846, 'inaxes': ax_rlocus.axes})() - ax_rlocus.set_xlim((-4.420022219849855, 0.9957380594313321)) - ax_rlocus.set_ylim((1.7825491928580846, 10.695295157148486)) + event = type('test', (object,), {'xdata': 14.7607954359, 'ydata': -35.6171379864, 'inaxes': ax_rlocus.axes})() + ax_rlocus.set_xlim((-10.813628105112421, 14.760795435937652)) + ax_rlocus.set_ylim((-35.61713798641108, 33.879716621220311)) plt.get_current_fig_manager().toolbar.mode = 'zoom rect' _RLClickDispatcher(event, system, fig, ax_rlocus, '-') @@ -62,10 +62,10 @@ def test_root_locus_zoom(self): zoom_y = ax_rlocus.lines[-2].get_data()[1][65:75] zoom_y = [abs(y) for y in zoom_y] - zoom_x_valid = [-2.23659192, -2.23659121, -2.23659103, -2.23659086, -2.23659068, -2.23659063, - -2.23659059, -2.23659055, -2.2365905, -2.23658766] - zoom_y_valid = [1.78253589, 1.78254318, 1.782545, 1.78254682, 1.78254864, 1.7825491, - 1.78254955, 1.78255001, 1.78255046, 1.7825796] + zoom_x_valid = [4.35145783, 4.49519318, 4.63559911 , 4.7728639 , 5.2937363 , 5.77581159, + 6.22564674 , 6.64818196 , 7.0472438 , 7.4258642] + zoom_y_valid = [19.34886165, 19.63109635, 19.90593613 , 20.17384159 , 21.18390302, + 22.11041787, 22.96863881, 23.76981422 , 24.52250243, 25.23338241] assert_array_almost_equal(zoom_x,zoom_x_valid) assert_array_almost_equal(zoom_y,zoom_y_valid) From 70ec8be87d5a9077ee5148adb83acf021186d71b Mon Sep 17 00:00:00 2001 From: David de Jong Date: Sun, 3 Jun 2018 18:40:14 +0200 Subject: [PATCH 19/21] changed root locus test indices. --- control/tests/rlocus_test.py | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/control/tests/rlocus_test.py b/control/tests/rlocus_test.py index 58550b1d2..464f04066 100644 --- a/control/tests/rlocus_test.py +++ b/control/tests/rlocus_test.py @@ -58,14 +58,12 @@ def test_root_locus_zoom(self): plt.get_current_fig_manager().toolbar.mode = 'zoom rect' _RLClickDispatcher(event, system, fig, ax_rlocus, '-') - zoom_x = ax_rlocus.lines[-2].get_data()[0][65:75] - zoom_y = ax_rlocus.lines[-2].get_data()[1][65:75] + zoom_x = ax_rlocus.lines[-2].get_data()[0][0:5] + zoom_y = ax_rlocus.lines[-2].get_data()[1][0:5] zoom_y = [abs(y) for y in zoom_y] - zoom_x_valid = [4.35145783, 4.49519318, 4.63559911 , 4.7728639 , 5.2937363 , 5.77581159, - 6.22564674 , 6.64818196 , 7.0472438 , 7.4258642] - zoom_y_valid = [19.34886165, 19.63109635, 19.90593613 , 20.17384159 , 21.18390302, - 22.11041787, 22.96863881, 23.76981422 , 24.52250243, 25.23338241] + zoom_x_valid = [-5. ,- 4.61281263, - 4.16689986, - 4.04122642, - 3.90736502] + zoom_y_valid = [0. ,0., 0., 0., 0.] assert_array_almost_equal(zoom_x,zoom_x_valid) assert_array_almost_equal(zoom_y,zoom_y_valid) From 03f87a35b683d38a4951119be332ea67a94a31c5 Mon Sep 17 00:00:00 2001 From: Richard Murray Date: Thu, 11 Apr 2019 21:33:56 -0700 Subject: [PATCH 20/21] Fix small typo from manual merge --- control/freqplot.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/control/freqplot.py b/control/freqplot.py index 746962baa..d264966ab 100644 --- a/control/freqplot.py +++ b/control/freqplot.py @@ -280,7 +280,7 @@ def bode_plot(syslist, omega=None, dB=None, Hz=None, deg=None, ax_phase.semilogx([Wcp, Wcp], [1e5, math.radians(phase_limit) + math.radians(pm)], color='k', linestyle=':', zorder=-20) - ax_phase.semilogx([Wcp, Wcp], [math.radians(phase_limit) + math.radians(pm) + ax_phase.semilogx([Wcp, Wcp], [math.radians(phase_limit) + math.radians(pm), math.radians(phase_limit)], color='k', zorder=-20) From 17c4a95c447f76ae2a6f72e594da3aa3cf4d5741 Mon Sep 17 00:00:00 2001 From: Richard Murray Date: Thu, 11 Apr 2019 22:22:25 -0700 Subject: [PATCH 21/21] fixed issues identified in review --- control/freqplot.py | 110 +++++++++++++++++++++++++++++--------------- control/sisotool.py | 12 +++-- control/timeresp.py | 2 +- 3 files changed, 81 insertions(+), 43 deletions(-) diff --git a/control/freqplot.py b/control/freqplot.py index d264966ab..6600a5b4a 100644 --- a/control/freqplot.py +++ b/control/freqplot.py @@ -217,8 +217,10 @@ def bode_plot(syslist, omega=None, dB=None, Hz=None, deg=None, # 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', + ax_mag = plt.subplot(211, + label='control-bode-magnitude') + ax_phase = plt.subplot(212, + label='control-bode-phase', sharex=ax_mag) # Magnitude plot @@ -226,10 +228,8 @@ def bode_plot(syslist, omega=None, dB=None, Hz=None, deg=None, 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()) @@ -259,7 +259,8 @@ def bode_plot(syslist, omega=None, dB=None, Hz=None, deg=None, if Hz: Wcg, Wcp = Wcg/(2*math.pi),Wcp/(2*math.pi) - ax_mag.axhline(y=0 if dB else 1, color='k', linestyle=':',zorder=-20) + ax_mag.axhline(y=0 if dB else 1, color='k', linestyle=':', + zorder=-20) ax_phase.axhline(y=phase_limit if deg else math.radians(phase_limit), color='k', linestyle=':', zorder=-20) mag_ylim = ax_mag.get_ylim() @@ -267,61 +268,87 @@ def bode_plot(syslist, omega=None, dB=None, Hz=None, deg=None, if pm != float('inf') and Wcp != float('nan'): if dB: - ax_mag.semilogx([Wcp, Wcp], [0.,-1e5],color='k', linestyle=':',zorder=-20) + ax_mag.semilogx([Wcp, Wcp], [0.,-1e5], + color='k', linestyle=':', + zorder=-20) else: - ax_mag.loglog([Wcp,Wcp], [1.,1e-8],color='k',linestyle=':',zorder=-20) + ax_mag.loglog([Wcp,Wcp], [1.,1e-8],color='k', + linestyle=':', zorder=-20) if deg: - ax_phase.semilogx([Wcp, Wcp], [1e5, phase_limit+pm], - color='k', linestyle=':', zorder=-20) - ax_phase.semilogx([Wcp, Wcp], [phase_limit + pm, phase_limit], + ax_phase.semilogx([Wcp, Wcp], + [1e5, phase_limit+pm], + color='k', linestyle=':', + zorder=-20) + ax_phase.semilogx([Wcp, Wcp], + [phase_limit + pm, phase_limit], color='k', zorder=-20) else: - ax_phase.semilogx([Wcp, Wcp], [1e5, math.radians(phase_limit) + - math.radians(pm)], - color='k', linestyle=':', zorder=-20) - ax_phase.semilogx([Wcp, Wcp], [math.radians(phase_limit) + math.radians(pm), - math.radians(phase_limit)], + ax_phase.semilogx([Wcp, Wcp], + [1e5, math.radians(phase_limit) + + math.radians(pm)], + color='k', linestyle=':', + zorder=-20) + ax_phase.semilogx([Wcp, Wcp], + [math.radians(phase_limit) + + math.radians(pm), + math.radians(phase_limit)], color='k', zorder=-20) if gm != float('inf') and Wcg != float('nan'): if dB: - ax_mag.semilogx([Wcg, Wcg], [-20.*np.log10(gm), -1e5], - color='k', linestyle=':',zorder=-20) + ax_mag.semilogx([Wcg, Wcg], + [-20.*np.log10(gm), -1e5], + color='k', linestyle=':', + zorder=-20) ax_mag.semilogx([Wcg, Wcg], [0,-20*np.log10(gm)], color='k', zorder=-20) else: - ax_mag.loglog([Wcg, Wcg], [1./gm,1e-8],color='k', linestyle=':', zorder=-20) - ax_mag.loglog([Wcg, Wcg], [1.,1./gm],color='k', zorder=-20) + ax_mag.loglog([Wcg, Wcg], + [1./gm,1e-8],color='k', + linestyle=':', zorder=-20) + ax_mag.loglog([Wcg, Wcg], + [1.,1./gm],color='k', zorder=-20) if deg: ax_phase.semilogx([Wcg, Wcg], [1e-8, phase_limit], - color='k', linestyle=':', zorder=-20) + color='k', linestyle=':', + zorder=-20) else: - ax_phase.semilogx([Wcg, Wcg], [1e-8, math.radians(phase_limit)], - color='k', linestyle=':', zorder=-20) + ax_phase.semilogx([Wcg, Wcg], + [1e-8, math.radians(phase_limit)], + color='k', linestyle=':', + zorder=-20) ax_mag.set_ylim(mag_ylim) ax_phase.set_ylim(phase_ylim) if sisotool: - ax_mag.text(0.04, 0.06, 'G.M.: %.2f %s\nFreq: %.2f %s' % + ax_mag.text(0.04, 0.06, + 'G.M.: %.2f %s\nFreq: %.2f %s' % (20*np.log10(gm) if dB else gm, - 'dB ' if dB else '', Wcg, 'Hz' if Hz else 'rad/s'), - horizontalalignment='left', verticalalignment='bottom', + 'dB ' if dB else '', + Wcg, 'Hz' if Hz else 'rad/s'), + horizontalalignment='left', + verticalalignment='bottom', transform=ax_mag.transAxes, fontsize=8 if int(matplotlib.__version__[0]) == 1 else 6) - ax_phase.text(0.04, 0.06, 'P.M.: %.2f %s\nFreq: %.2f %s' % + ax_phase.text(0.04, 0.06, + 'P.M.: %.2f %s\nFreq: %.2f %s' % (pm if deg else math.radians(pm), - 'deg' if deg else 'rad', Wcp, 'Hz' if Hz else 'rad/s'), - horizontalalignment='left', verticalalignment='bottom', - transform=ax_phase.transAxes, + 'deg' if deg else 'rad', + Wcp, 'Hz' if Hz else 'rad/s'), + horizontalalignment='left', + verticalalignment='bottom', + transform=ax_phase.transAxes, fontsize=8 if int(matplotlib.__version__[0]) == 1 else 6) else: plt.suptitle('Gm = %.2f %s(at %.2f %s), Pm = %.2f %s (at %.2f %s)' % (20*np.log10(gm) if dB else gm, - 'dB ' if dB else '\b', Wcg, 'Hz' if Hz else 'rad/s', - pm if deg else math.radians(pm),'deg' if deg else 'rad', + 'dB ' if dB else '\b', + Wcg, 'Hz' if Hz else 'rad/s', + pm if deg else math.radians(pm), + 'deg' if deg else 'rad', Wcp, 'Hz' if Hz else 'rad/s')) if nyquistfrq_plot: @@ -336,12 +363,19 @@ def gen_zero_centered_series(val_min, val_max, period): return np.arange(v1, v2 + 1) * period if deg: ylim = ax_phase.get_ylim() - ax_phase.set_yticks(gen_zero_centered_series(ylim[0], ylim[1], 45.)) - ax_phase.set_yticks(gen_zero_centered_series(ylim[0], ylim[1], 15.), minor=True) + ax_phase.set_yticks(gen_zero_centered_series(ylim[0], + ylim[1], 45.)) + ax_phase.set_yticks(gen_zero_centered_series(ylim[0], + ylim[1], 15.), + minor=True) else: ylim = ax_phase.get_ylim() - ax_phase.set_yticks(gen_zero_centered_series(ylim[0], ylim[1], math.pi / 4.)) - ax_phase.set_yticks(gen_zero_centered_series(ylim[0], ylim[1], math.pi / 12.), + ax_phase.set_yticks(gen_zero_centered_series(ylim[0], + ylim[1], + math.pi / 4.)) + ax_phase.set_yticks(gen_zero_centered_series(ylim[0], + ylim[1], + math.pi / 12.), minor=True) ax_phase.grid(False if margins else True, which='both') # ax_mag.grid(which='minor', alpha=0.3) @@ -350,7 +384,8 @@ def gen_zero_centered_series(val_min, val_max, period): # ax_phase.grid(which='major', alpha=0.9) # Label the frequency axis - ax_phase.set_xlabel("Frequency (Hz)" if Hz else "Frequency (rad/sec)") + ax_phase.set_xlabel("Frequency (Hz)" if Hz + else "Frequency (rad/sec)") if len(syslist) == 1: return mags[0], phases[0], omegas[0] @@ -358,7 +393,8 @@ def gen_zero_centered_series(val_min, val_max, period): return mags, phases, omegas -def nyquist_plot(syslist, omega=None, Plot=True, color=None, labelFreq=0, *args, **kwargs): +def nyquist_plot(syslist, omega=None, Plot=True, color=None, + labelFreq=0, *args, **kwargs): """ Nyquist plot for a system diff --git a/control/sisotool.py b/control/sisotool.py index 916eef1cd..7a312cf5c 100644 --- a/control/sisotool.py +++ b/control/sisotool.py @@ -7,11 +7,13 @@ import matplotlib.pyplot as plt import warnings -def sisotool(sys, kvect = None, xlim_rlocus = None, ylim_rlocus = None, plotstr_rlocus = 'b' if int(matplotlib.__version__[0]) == 1 else 'C0',rlocus_grid = False, omega = None, dB = None, Hz = None, deg = None, omega_limits = None, omega_num = None,margins_bode = True, tvect=None): - - """Sisotool - - Sisotool style collection of plots inspired by the matlab sisotool. +def sisotool(sys, kvect = None, xlim_rlocus = None, ylim_rlocus = None, + plotstr_rlocus = 'b' if int(matplotlib.__version__[0]) == 1 else 'C0', + rlocus_grid = False, omega = None, dB = None, Hz = None, + deg = None, omega_limits = None, omega_num = None, + margins_bode = True, tvect=None): + """ + Sisotool style collection of plots inspired by MATLAB's sisotool. The left two plots contain the bode magnitude and phase diagrams. The top right plot is a clickable root locus plot, clicking on the root locus will change the gain of the system. The bottom left plot diff --git a/control/timeresp.py b/control/timeresp.py index e9ca866fd..14e7072d2 100644 --- a/control/timeresp.py +++ b/control/timeresp.py @@ -375,7 +375,7 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False, def _get_ss_simo(sys, input=None, output=None): """Return a SISO or SIMO state-space version of sys - If input is not sfpecified, select first input and issue warning + If input is not specified, select first input and issue warning """ sys_ss = _convertToStateSpace(sys) if sys_ss.issiso():