diff --git a/control/config.py b/control/config.py index 188a4efd7..98198aae7 100644 --- a/control/config.py +++ b/control/config.py @@ -11,7 +11,7 @@ bode_dB = False # Bode plot magnitude units bode_deg = True # Bode Plot phase units bode_Hz = False # Bode plot frequency units -bode_number_of_samples = None # Bode plot number of samples +bode_number_of_samples = None # Bode plot number of samples bode_feature_periphery_decade = 1.0 # Bode plot feature periphery in decades # Set defaults to match MATLAB @@ -29,7 +29,7 @@ def use_matlab_defaults(): def use_fbs_defaults(): """ Use `Astrom and Murray `_ compatible settings - * Bode plots plot gain in powers of ten, phase in degrees, + * Bode plots plot gain in powers of ten, phase in degrees, frequency in Hertz """ # Bode plot defaults diff --git a/control/exception.py b/control/exception.py index 7579df4f6..2c4f29704 100644 --- a/control/exception.py +++ b/control/exception.py @@ -2,7 +2,7 @@ # # Author: Richard M. Murray # Date: 31 May 2010 -# +# # This file contains definitions of standard exceptions for the control package # # Copyright (c) 2010 by California Institute of Technology @@ -14,16 +14,16 @@ # # 1. Redistributions of source code must retain the above copyright # notice, this list of conditions and the following disclaimer. -# +# # 2. Redistributions in binary form must reproduce the above copyright # notice, this list of conditions and the following disclaimer in the # documentation and/or other materials provided with the distribution. -# +# # 3. Neither the name of the California Institute of Technology nor # the names of its contributors may be used to endorse or promote # products derived from this software without specific prior # written permission. -# +# # THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS # "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT # LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS @@ -36,19 +36,19 @@ # OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT # OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF # SUCH DAMAGE. -# +# # $Id$ -class ControlSlycot(Exception): +class ControlSlycot(Exception): """Exception for Slycot import. Used when we can't import a function from the slycot package""" pass -class ControlDimension(Exception): +class ControlDimension(Exception): """Raised when dimensions of system objects are not correct""" pass -class ControlArgument(Exception): +class ControlArgument(Exception): """Raised when arguments to a function are not correct""" pass diff --git a/control/freqplot.py b/control/freqplot.py index 3e05d8937..9bb713216 100644 --- a/control/freqplot.py +++ b/control/freqplot.py @@ -78,11 +78,11 @@ def bode_plot(syslist, omega=None, dB=None, Hz=None, deg=None, If True, plot phase in degrees (else radians) Plot : boolean If True, plot magnitude and phase - omega_limits: tuple, list, ... of two values + 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 + number of samples *args, **kwargs: Additional options to matplotlib (color, linestyle, etc) @@ -120,7 +120,7 @@ def bode_plot(syslist, omega=None, dB=None, Hz=None, deg=None, # If argument was a singleton, turn it into a list if (not getattr(syslist, '__iter__', False)): syslist = (syslist,) - + if omega is None: if omega_limits is None: # Select a default range if none is provided @@ -133,7 +133,7 @@ def bode_plot(syslist, omega=None, dB=None, Hz=None, deg=None, omega = sp.logspace(np.log10(omega_limits[0]), np.log10(omega_limits[1]), num=omega_num, endpoint=True) else: omega = sp.logspace(np.log10(omega_limits[0]), np.log10(omega_limits[1]), endpoint=True) - + mags, phases, omegas, nyquistfrqs = [], [], [], [] for sys in syslist: if (sys.inputs > 1 or sys.outputs > 1): @@ -142,8 +142,8 @@ def bode_plot(syslist, omega=None, dB=None, Hz=None, deg=None, else: omega_sys = np.array(omega) if sys.isdtime(True): - nyquistfrq = 2. * np.pi * 1. / sys.dt / 2. - omega_sys = omega_sys[omega_sys < nyquistfrq] + nyquistfrq = 2. * np.pi * 1. / sys.dt / 2. + omega_sys = omega_sys[omega_sys < nyquistfrq] # TODO: What distance to the Nyquist frequency is appropriate? else: nyquistfrq = None @@ -179,7 +179,7 @@ def bode_plot(syslist, omega=None, dB=None, Hz=None, deg=None, plt.hold(True); if nyquistfrq_plot: ax_mag.axvline(nyquistfrq_plot, color=pltline[0].get_color()) - + # Add a grid to the plot + labeling ax_mag.grid(True, which='both') ax_mag.set_ylabel("Magnitude (dB)" if dB else "Magnitude") @@ -194,30 +194,31 @@ def bode_plot(syslist, omega=None, dB=None, Hz=None, deg=None, ax_phase.hold(True); if nyquistfrq_plot: ax_phase.axvline(nyquistfrq_plot, color=pltline[0].get_color()) - + # Add a grid to the plot + labeling - ax_phase.set_ylabel("Phase (deg)" if deg else "Phase (rad)") + ax_phase.set_ylabel("Phase (deg)" if deg else "Phase (rad)") + def genZeroCenteredSeries(val_min, val_max, period): v1 = np.ceil(val_min / period - 0.2) v2 = np.floor(val_max / period + 0.2) - return np.arange(v1, v2 + 1) * period + return np.arange(v1, v2 + 1) * period if deg: ylim = ax_phase.get_ylim() - ax_phase.set_yticks(genZeroCenteredSeries(ylim[0], ylim[1], 45.)) - ax_phase.set_yticks(genZeroCenteredSeries(ylim[0], ylim[1], 15.), minor=True) + ax_phase.set_yticks(genZeroCenteredSeries(ylim[0], ylim[1], 45.)) + ax_phase.set_yticks(genZeroCenteredSeries(ylim[0], ylim[1], 15.), minor=True) else: ylim = ax_phase.get_ylim() - ax_phase.set_yticks(genZeroCenteredSeries(ylim[0], ylim[1], np.pi / 4.)) + ax_phase.set_yticks(genZeroCenteredSeries(ylim[0], ylim[1], np.pi / 4.)) ax_phase.set_yticks(genZeroCenteredSeries(ylim[0], ylim[1], np.pi / 12.), minor=True) ax_phase.grid(True, which='both') - # ax_mag.grid(which='minor', alpha=0.3) + # ax_mag.grid(which='minor', alpha=0.3) # ax_mag.grid(which='major', alpha=0.9) - # ax_phase.grid(which='minor', alpha=0.3) - # ax_phase.grid(which='major', alpha=0.9) - + # ax_phase.grid(which='minor', alpha=0.3) + # 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] else: @@ -314,7 +315,7 @@ def nyquist_plot(syslist, omega=None, Plot=True, color='b', # instead of 1.0, and this would otherwise be # truncated to 0. plt.text(xpt, ypt, - ' ' + str(int(np.round(f / 1000 ** pow1000, 0))) + + ' ' + str(int(np.round(f / 1000 ** pow1000, 0))) + ' ' + prefix + 'Hz') return x, y, omega @@ -394,18 +395,18 @@ def default_frequency_range(syslist, Hz=None, number_of_samples=None, feature_pe syslist : list of LTI List of linear input/output systems (single system is OK) Hz: boolean - If True, the limits (first and last value) of the frequencies - are set to full decades in Hz so it fits plotting with logarithmic + If True, the limits (first and last value) of the frequencies + are set to full decades in Hz so it fits plotting with logarithmic scale in Hz otherwise in rad/s. Omega is always returned in rad/sec. number_of_samples: int Number of samples to generate feature_periphery_decade: float - Defines how many decades shall be included in the frequency range on - both sides of features (poles, zeros). + Defines how many decades shall be included in the frequency range on + both sides of features (poles, zeros). Example: If there is a feature, e.g. a pole, at 1Hz and feature_periphery_decade=1. - then the range of frequencies shall span 0.1 .. 10 Hz. + then the range of frequencies shall span 0.1 .. 10 Hz. The default value is read from config.bode_feature_periphery_decade. - + Returns ------- omega : array @@ -425,14 +426,14 @@ def default_frequency_range(syslist, Hz=None, number_of_samples=None, feature_pe # Set default values for options from . import config - if (number_of_samples is None): + if (number_of_samples is None): number_of_samples = config.bode_number_of_samples - if (feature_periphery_decade is None): - feature_periphery_decade = config.bode_feature_periphery_decade + if (feature_periphery_decade is None): + feature_periphery_decade = config.bode_feature_periphery_decade # Find the list of all poles and zeros in the systems features = np.array(()) - freq_interesting = [] + freq_interesting = [] # detect if single sys passed by checking if it is sequence-like if (not getattr(syslist, '__iter__', False)): @@ -443,34 +444,34 @@ def default_frequency_range(syslist, Hz=None, number_of_samples=None, feature_pe # Add new features to the list if sys.isctime(): features_ = np.concatenate((np.abs(sys.pole()), - np.abs(sys.zero()))) + np.abs(sys.zero()))) # Get rid of poles and zeros at the origin features_ = features_[features_ != 0.0]; features = np.concatenate((features, features_)) elif sys.isdtime(strict=True): fn = np.pi * 1. / sys.dt - # TODO: What distance to the Nyquist frequency is appropriate? + # TODO: What distance to the Nyquist frequency is appropriate? freq_interesting.append(fn * 0.9) features_ = np.concatenate((sys.pole(), - sys.zero())) - # Get rid of poles and zeros + sys.zero())) + # Get rid of poles and zeros # * at the origin and real <= 0 & imag==0: log! # * at 1.: would result in omega=0. (logaritmic plot!) features_ = features_[(features_.imag != 0.0) | (features_.real > 0.)] features_ = features_[np.bitwise_not((features_.imag == 0.0) & (np.abs(features_.real - 1.0) < 1.e-10))] # TODO: improve features__ = np.abs(np.log(features_) / (1.j * sys.dt)) - features = np.concatenate((features, features__)) + features = np.concatenate((features, features__)) else: # TODO - raise NotImplementedError('type of system in not implemented now') + raise NotImplementedError('type of system in not implemented now') except: pass # Make sure there is at least one point in the range - if (features.shape[0] == 0): + if (features.shape[0] == 0): features = np.array([1.]); if Hz: diff --git a/control/margins.py b/control/margins.py index 5fb5a8b40..655a38619 100644 --- a/control/margins.py +++ b/control/margins.py @@ -130,7 +130,7 @@ def stability_margins(sysdata, returnall=False, epsw=1e-8): sys = sysdata elif getattr(sysdata, '__iter__', False) and len(sysdata) == 3: mag, phase, omega = sysdata - sys = frdata.FRD(mag * np.exp(1j * phase * np.pi/180), + sys = frdata.FRD(mag * np.exp(1j * phase * np.pi/180), omega, smooth=True) else: sys = xferfcn._convertToTransferFunction(sysdata) @@ -163,7 +163,7 @@ def stability_margins(sysdata, returnall=False, epsw=1e-8): # evaluate response at remaining frequencies, to test for phase 180 vs 0 resp_w_180 = np.real(np.polyval(sys.num[0][0], 1.j*w_180) / np.polyval(sys.den[0][0], 1.j*w_180)) - #print ('resp_w_180', resp_w_180) + #print ('resp_w_180', resp_w_180) # only keep frequencies where the negative real axis is crossed w_180 = w_180[np.real(resp_w_180) < 0.0] @@ -192,7 +192,7 @@ def stability_margins(sysdata, returnall=False, epsw=1e-8): # find the solutions, for positive omega, and only real ones wstab = np.roots(test_wstab) #print('wstabr', wstab) - wstab = np.real(wstab[(np.imag(wstab) == 0) * + wstab = np.real(wstab[(np.imag(wstab) == 0) * (np.real(wstab) >= 0)]) #print('wstab', wstab) @@ -224,7 +224,7 @@ def dstab(w): wc = np.array( [ sp.optimize.brentq(mod, sys.omega[i], sys.omega[i+1]) for i in widx if i+1 < len(sys.omega)]) - + # find the phase crossings ang(H(jw) == -180 widx = np.where(np.diff(np.sign(arg(sys.omega))))[0] #print('widx (180)', widx, sys.omega[widx]) @@ -246,17 +246,17 @@ def dstab(w): for i in widx if i+1 < len(sys.omega) and np.diff(np.diff(dstab(sys.omega[i-1:i+2])))[0] > 0 ]) #print('wstabf0', wstab) - wstab = wstab[(wstab >= sys.omega[0]) * + wstab = wstab[(wstab >= sys.omega[0]) * (wstab <= sys.omega[-1])] #print ('wstabf', wstab) - + # margins, as iterables, converted frdata and xferfcn calculations to # vector for this GM = 1/np.abs(sys.evalfr(w_180)[0][0]) SM = np.abs(sys.evalfr(wstab)[0][0]+1) PM = np.angle(sys.evalfr(wc)[0][0], deg=True) + 180 - + if returnall: return GM, PM, SM, w_180, wc, wstab else: @@ -337,7 +337,7 @@ def margin(*args): Gain crossover frequency (corresponding to phase margin) Wcp : float Phase crossover frequency (corresponding to gain margin) (in rad/sec) - + Margins are of SISO open-loop. If more than one crossover frequency is detected, returns the lowest corresponding margin. diff --git a/control/tests/bdalg_test.py b/control/tests/bdalg_test.py index e04c96b4e..7db39d54f 100644 --- a/control/tests/bdalg_test.py +++ b/control/tests/bdalg_test.py @@ -40,7 +40,7 @@ def testScalarSS(self): ans1 = feedback(self.x1, self.sys2) ans2 = feedback(self.x1, self.sys2, 1.) - + np.testing.assert_array_almost_equal(ans1.A, [[-1.5, 4.], [13., 2.]]) np.testing.assert_array_almost_equal(ans1.B, [[2.5], [-10.]]) np.testing.assert_array_almost_equal(ans1.C, [[-2.5, 0.]]) @@ -77,7 +77,7 @@ def testScalarTF(self): def testSSScalar(self): """State space system with scalar feedback block.""" - + ans1 = feedback(self.sys2, self.x1) ans2 = feedback(self.sys2, self.x1, 1.) @@ -98,7 +98,7 @@ def testSSSS1(self): np.testing.assert_array_almost_equal(ans1.A, [[1., 4., -1., 0.], [3., 2., 4., 0.], [1., 0., 1., 4.], [-4., 0., 3., 2]]) - np.testing.assert_array_almost_equal(ans1.B, [[1.], [-4.], [0.], [0.]]) + np.testing.assert_array_almost_equal(ans1.B, [[1.], [-4.], [0.], [0.]]) np.testing.assert_array_almost_equal(ans1.C, [[1., 0., 0., 0.]]) np.testing.assert_array_almost_equal(ans1.D, [[0.]]) np.testing.assert_array_almost_equal(ans2.A, [[1., 4., 1., 0.], @@ -115,7 +115,7 @@ def testSSSS2(self): [[-2.]]) sys4 = StateSpace([[-3., -2.], [1., 4.]], [[-2.], [-6.]], [[2., -3.]], [[3.]]) - + ans1 = feedback(sys3, sys4) ans2 = feedback(sys3, sys4, 1.) @@ -144,7 +144,7 @@ def testSSSS2(self): def testSSTF(self): """State space system with transfer function feedback block.""" - + # This functionality is not implemented yet. pass diff --git a/control/tests/freqresp_test.py b/control/tests/freqresp_test.py index cc25bc37c..c7660e919 100644 --- a/control/tests/freqresp_test.py +++ b/control/tests/freqresp_test.py @@ -3,8 +3,8 @@ # freqresp_test.py - test frequency response functions # RMM, 30 May 2016 (based on timeresp_test.py) # -# This is a rudimentary set of tests for frequency response functions, -# including bode plots. +# This is a rudimentary set of tests for frequency response functions, +# including bode plots. import unittest import numpy as np @@ -53,7 +53,7 @@ def test_mimo(self): #bode(sysMIMO) # - should throw not implemented exception #bode(tfMIMO) # - should throw not implemented exception - + #plt.figure(3) #plt.semilogx(self.omega,20*np.log10(np.squeeze(frq[0]))) diff --git a/control/tests/margin_test.py b/control/tests/margin_test.py index 3ef74e7d4..c8b52e2f1 100644 --- a/control/tests/margin_test.py +++ b/control/tests/margin_test.py @@ -25,12 +25,12 @@ def setUp(self): (StateSpace([[1., 4.], [3., 2.]], [[1.], [-4.]], [[1., 0.]], [[0.]]), [], [], [147.0743], [2.5483]), - ((8.75*(4*s**2+0.4*s+1))/((100*s+1)*(s**2+0.22*s+1)) * - 1./(s**2/(10.**2)+2*0.04*s/10.+1), + ((8.75*(4*s**2+0.4*s+1))/((100*s+1)*(s**2+0.22*s+1)) * + 1./(s**2/(10.**2)+2*0.04*s/10.+1), [2.2716], [10.0053], [97.5941, 360-157.7904, 134.7359], [0.0850, 0.9373, 1.0919])) - - + + self.sys1 = TransferFunction([1, 2], [1, 2, 3]) # alternative # sys1 = tf([1, 2], [1, 2, 3]) @@ -55,7 +55,7 @@ def test_stability_margins(self): np.testing.assert_array_almost_equal( out[out != np.array(None)], outf[outf != np.array(None)], 2) - + # final one with fixed values np.testing.assert_array_almost_equal( [gm, pm, sm, wg, wp, ws], @@ -77,7 +77,7 @@ def test_stability_margins_all(self): print(res, '\n', comp) np.testing.assert_array_almost_equal( res, comp, 2) - + def test_phase_crossover_frequencies(self): omega, gain = phase_crossover_frequencies(self.sys2) np.testing.assert_array_almost_equal(omega, [1.73205, 0.]) @@ -163,8 +163,8 @@ def test_nocross(self): out1b = stability_margins(FRD(h1, omega)) out2b = stability_margins(FRD(h2, omega)) out3b = stability_margins(FRD(h3, omega)) - - + + def test_suite(): return unittest.TestLoader().loadTestsFromTestCase(TestMargin) diff --git a/control/tests/nichols_test.py b/control/tests/nichols_test.py index 9aca2c125..297c63f2d 100644 --- a/control/tests/nichols_test.py +++ b/control/tests/nichols_test.py @@ -17,8 +17,8 @@ def setUp(self): B = [[1.], [-3.], [-2.]] C = [[4., 2., -3.]] D = [[0.]] - - self.sys = StateSpace(A, B, C, D) + + self.sys = StateSpace(A, B, C, D) def testNicholsPlain(self): """Generate a Nichols plot.""" @@ -32,6 +32,6 @@ def testNgrid(self): def suite(): return unittest.TestLoader().loadTestsFromTestCase(TestStateSpace) - + if __name__ == "__main__": unittest.main() diff --git a/control/tests/phaseplot_test.py b/control/tests/phaseplot_test.py index d584858d1..e1cc25a35 100644 --- a/control/tests/phaseplot_test.py +++ b/control/tests/phaseplot_test.py @@ -24,19 +24,19 @@ def testInvPendNoSims(self): phase_plot(self.invpend_ode, (-6,6,10), (-6,6,10)); def testInvPendSims(self): - phase_plot(self.invpend_ode, (-6,6,10), (-6,6,10), + phase_plot(self.invpend_ode, (-6,6,10), (-6,6,10), X0 = ([1,1], [-1,1])); def testInvPendTimePoints(self): - phase_plot(self.invpend_ode, (-6,6,10), (-6,6,10), + phase_plot(self.invpend_ode, (-6,6,10), (-6,6,10), X0 = ([1,1], [-1,1]), T=np.linspace(0,5,100)); def testInvPendLogtime(self): - phase_plot(self.invpend_ode, X0 = + phase_plot(self.invpend_ode, X0 = [ [-2*pi, 1.6], [-2*pi, 0.5], [-1.8, 2.1], [-1, 2.1], [4.2, 2.1], [5, 2.1], [2*pi, -1.6], [2*pi, -0.5], [1.8, -2.1], - [1, -2.1], [-4.2, -2.1], [-5, -2.1] ], + [1, -2.1], [-4.2, -2.1], [-5, -2.1] ], T = np.linspace(0, 40, 200), logtime=(3, 0.7), verbose=False) @@ -47,7 +47,7 @@ def testInvPendAuto(self): def testOscillatorParams(self): m = 1; b = 1; k = 1; # default values - phase_plot(self.oscillator_ode, timepts = [0.3, 1, 2, 3], X0 = + phase_plot(self.oscillator_ode, timepts = [0.3, 1, 2, 3], X0 = [[-1,1], [-0.3,1], [0,1], [0.25,1], [0.5,1], [0.7,1], [1,1], [1.3,1], [1,-1], [0.3,-1], [0,-1], [-0.25,-1], [-0.5,-1], [-0.7,-1], [-1,-1], [-1.3,-1]], diff --git a/control/tests/run_all.py b/control/tests/run_all.py index 0b18c0827..b21248432 100755 --- a/control/tests/run_all.py +++ b/control/tests/run_all.py @@ -26,7 +26,7 @@ def test_all(verbosity=0): t.run(tests) print('Completed tests in', mod) - except: + except: testModules = findTests('./tests/') # Now go through each module and run all of its tests. @@ -46,13 +46,13 @@ def test_all(verbosity=0): 'a proper suite() function') t=unittest.TextTestRunner(verbosity=verbosity) - t.run(unittest.TestSuite(unittest.TestSuite(suiteList))) + t.run(unittest.TestSuite(unittest.TestSuite(suiteList))) print('Completed tests in', mod) def findTests(testdir = './', pattern = "[^.#]*_test.py$"): - """Since python <2.7 doesn't have test discovery, this finds tests in the + """Since python <2.7 doesn't have test discovery, this finds tests in the provided directory. The default is to check the current directory. Any files - that match test* or Test* are considered unittest modules and checked for + that match test* or Test* are considered unittest modules and checked for a module.suite() function (in tests()).""" # Get list of files in test directory diff --git a/control/timeresp.py b/control/timeresp.py index ed95e319a..30da52462 100644 --- a/control/timeresp.py +++ b/control/timeresp.py @@ -2,7 +2,7 @@ """Time domain simulation. This file contains a collection of functions that calculate time -responses for linear systems. +responses for linear systems. See doc/conventions.rst#time-series-conventions_ for more information on how time series data are represented. diff --git a/doc/intro.rst b/doc/intro.rst index a55913291..5310be55b 100644 --- a/doc/intro.rst +++ b/doc/intro.rst @@ -58,10 +58,11 @@ and verifying that no error message appears. It may be necessary to install `lapack` library. More information on the slycot package can be obtained from the `slycot project page `_. -For users with a working the Anaconda distribution of Python, the following -command can be used:: +For users with the Anaconda distribution of Python, the following +commands can be used:: - conda install -c python-control -c cyclus slycot python-control + conda install numpy scipy matplotlib # if not yet installed + conda install -c python-control -c cyclus slycot control This installs `slycot` and `python-control` from the `python-control` channel and uses the `cyclus` channel to obtain the required `lapack`