From 573b75f635964839831619c214e8d8af45593e55 Mon Sep 17 00:00:00 2001 From: Joris Geysens Date: Wed, 27 Mar 2019 12:43:44 +0100 Subject: [PATCH 1/6] Configure flake8. Fix flake8 errors where possible. --- .travis.yml | 3 +- control/__init__.py | 46 ++++---- control/bdalg.py | 69 +++++------ control/bench/time_freqresp.py | 4 +- control/canonical.py | 22 ++-- control/config.py | 20 +++- control/ctrlutil.py | 8 +- control/delay.py | 30 ++--- control/dtime.py | 8 +- control/exception.py | 10 +- control/frdata.py | 119 +++++++++---------- control/freqplot.py | 64 +++++----- control/grid.py | 93 ++++++++------- control/lti.py | 38 +++--- control/margins.py | 82 ++++++------- control/mateqn.py | 207 ++++++++++++++++----------------- control/matlab/__init__.py | 44 +++---- control/matlab/timeresp.py | 26 +++-- control/matlab/wrappers.py | 45 ++++--- control/modelsimp.py | 142 +++++++++++----------- control/nichols.py | 27 ++--- control/phaseplot.py | 163 +++++++++++++------------- control/pzmap.py | 12 +- control/rlocus.py | 67 ++++++----- control/robust.py | 137 +++++++++++----------- control/statefbk.py | 154 ++++++++++++------------ control/statesp.py | 161 +++++++++++-------------- control/tests/matlab_test.py | 4 +- control/timeresp.py | 55 ++++----- control/xferfcn.py | 150 ++++++++++-------------- setup.cfg | 7 ++ setup.py | 3 +- 32 files changed, 1003 insertions(+), 1017 deletions(-) diff --git a/.travis.yml b/.travis.yml index a28f57419..b7eb18ba1 100644 --- a/.travis.yml +++ b/.travis.yml @@ -15,7 +15,7 @@ python: # Test against multiple version of SciPy, with and without slycot # -# Because there were significant changes in SciPy between v0 and v1, we +# Because there were significant changes in SciPy between v0 and v1, we # test against both of these using the Travis CI environment capability # # We also want to test with and without slycot @@ -73,6 +73,7 @@ install: # command to run tests script: + - flake8 control/ - 'if [ $SLYCOT != "" ]; then python -c "import slycot"; fi' - coverage run setup.py test diff --git a/control/__init__.py b/control/__init__.py index 99cf9a73d..0f94a151c 100644 --- a/control/__init__.py +++ b/control/__init__.py @@ -46,34 +46,34 @@ # Import functions from within the control system library # Note: the functions we use are specified as __all__ variables in the modules -from .bdalg import * -from .delay import * -from .dtime import * -from .freqplot import * -from .lti import * -from .margins import * -from .mateqn import * -from .modelsimp import * -from .nichols import * -from .phaseplot import * -from .pzmap import * -from .rlocus import * -from .statefbk import * -from .statesp import * -from .timeresp import * -from .xferfcn import * -from .ctrlutil import * -from .frdata import * -from .canonical import * -from .robust import * -from .config import * +from .bdalg import * # noqa F403, F401 +from .delay import * # noqa F403, F401 +from .dtime import * # noqa F403, F401 +from .freqplot import * # noqa F403, F401 +from .lti import * # noqa F403, F401 +from .margins import * # noqa F403, F401 +from .mateqn import * # noqa F403, F401 +from .modelsimp import * # noqa F403, F401 +from .nichols import * # noqa F403, F401 +from .phaseplot import * # noqa F403, F401 +from .pzmap import * # noqa F403, F401 +from .rlocus import * # noqa F403, F401 +from .statefbk import * # noqa F403, F401 +from .statesp import * # noqa F403, F401 +from .timeresp import * # noqa F403, F401 +from .xferfcn import * # noqa F403, F401 +from .ctrlutil import * # noqa F403, F401 +from .frdata import * # noqa F403, F401 +from .canonical import * # noqa F403, F401 +from .robust import * # noqa F403, F401 +from .config import * # noqa F403, F401 # Exceptions -from .exception import * +from .exception import * # noqa F403, F401 # Version information try: - from ._version import __version__, __commit__ + from ._version import __version__, __commit__ # noqa F403, F401 except ImportError: __version__ = "dev" diff --git a/control/bdalg.py b/control/bdalg.py index 0f4a14c1a..a7e865343 100644 --- a/control/bdalg.py +++ b/control/bdalg.py @@ -61,8 +61,10 @@ __all__ = ['series', 'parallel', 'negate', 'feedback', 'append', 'connect'] + def series(sys1, *sysn): - """Return the series connection (... \* sys3 \*) sys2 \* sys1 + """ + Return the series connection (... \* sys3 \*) sys2 \* sys1 Parameters ---------- @@ -103,7 +105,8 @@ def series(sys1, *sysn): """ from functools import reduce - return reduce(lambda x, y:y*x, sysn, sys1) + return reduce(lambda x, y: y * x, sysn, sys1) + def parallel(sys1, *sysn): """ @@ -148,7 +151,8 @@ def parallel(sys1, *sysn): """ from functools import reduce - return reduce(lambda x, y:x+y, sysn, sys1) + return reduce(lambda x, y: x + y, sysn, sys1) + def negate(sys): """ @@ -178,9 +182,10 @@ def negate(sys): """ - return -sys; + return -sys -#! TODO: expand to allow sys2 default to work in MIMO case? + +# TODO: expand to allow sys2 default to work in MIMO case? def feedback(sys1, sys2=1, sign=-1): """ Feedback interconnection between two I/O systems. @@ -226,14 +231,10 @@ def feedback(sys1, sys2=1, sign=-1): """ # Check for correct input types. - if not isinstance(sys1, (int, float, complex, np.number, - tf.TransferFunction, ss.StateSpace, frd.FRD)): - raise TypeError("sys1 must be a TransferFunction, StateSpace " + - "or FRD object, or a scalar.") - if not isinstance(sys2, (int, float, complex, np.number, - tf.TransferFunction, ss.StateSpace, frd.FRD)): - raise TypeError("sys2 must be a TransferFunction, StateSpace " + - "or FRD object, or a scalar.") + if not isinstance(sys1, (int, float, complex, np.number, tf.TransferFunction, ss.StateSpace, frd.FRD)): + raise TypeError("sys1 must be a TransferFunction, StateSpace or FRD object, or a scalar.") + if not isinstance(sys2, (int, float, complex, np.number, tf.TransferFunction, ss.StateSpace, frd.FRD)): + raise TypeError("sys2 must be a TransferFunction, StateSpace or FRD object, or a scalar.") # If sys1 is a scalar, convert it to the appropriate LTI type so that we can # its feedback member function. @@ -244,15 +245,15 @@ def feedback(sys1, sys2=1, sign=-1): sys1 = ss._convertToStateSpace(sys1) elif isinstance(sys2, frd.FRD): sys1 = ss._convertToFRD(sys1) - else: # sys2 is a scalar. + else: # sys2 is a scalar. sys1 = tf._convert_to_transfer_function(sys1) sys2 = tf._convert_to_transfer_function(sys2) return sys1.feedback(sys2, sign) -def append(*sys): - '''append(sys1, sys2, ..., sysn) +def append(*sys): + """ Group models by appending their inputs and outputs Forms an augmented system model, and appends the inputs and @@ -278,16 +279,17 @@ def append(*sys): >>> sys2 = ss("-1.", "1.", "1.", "0.") >>> sys = append(sys1, sys2) - .. todo:: - also implement for transfer function, zpk, etc. - ''' + # TODO : also implement for transfer function, zpk, etc. + """ + s1 = sys[0] for s in sys[1:]: s1 = s1.append(s) return s1 + def connect(sys, Q, inputv, outputv): - ''' + """ Index-base interconnection of system The system sys is a system typically constructed with append, with @@ -324,23 +326,24 @@ def connect(sys, Q, inputv, outputv): >>> sys = append(sys1, sys2) >>> Q = sp.mat([ [ 1, 2], [2, -1] ]) # basically feedback, output 2 in 1 >>> sysc = connect(sys, Q, [2], [1, 2]) - ''' + """ + # first connect - K = sp.zeros( (sys.inputs, sys.outputs) ) + K = sp.zeros((sys.inputs, sys.outputs)) for r in sp.array(Q).astype(int): - inp = r[0]-1 + inp = r[0] - 1 for outp in r[1:]: - if outp > 0 and outp <= sys.outputs: - K[inp,outp-1] = 1. + if 0 < outp <= sys.outputs: + K[inp, outp - 1] = 1. elif outp < 0 and -outp >= -sys.outputs: - K[inp,-outp-1] = -1. + K[inp, -outp - 1] = -1. sys = sys.feedback(sp.matrix(K), sign=1) # now trim - Ytrim = sp.zeros( (len(outputv), sys.outputs) ) - Utrim = sp.zeros( (sys.inputs, len(inputv)) ) - for i,u in enumerate(inputv): - Utrim[u-1,i] = 1. - for i,y in enumerate(outputv): - Ytrim[i,y-1] = 1. - return sp.matrix(Ytrim)*sys*sp.matrix(Utrim) + Ytrim = sp.zeros((len(outputv), sys.outputs)) + Utrim = sp.zeros((sys.inputs, len(inputv))) + for i, u in enumerate(inputv): + Utrim[u - 1, i] = 1. + for i, y in enumerate(outputv): + Ytrim[i, y - 1] = 1. + return sp.matrix(Ytrim) * sys * sp.matrix(Utrim) diff --git a/control/bench/time_freqresp.py b/control/bench/time_freqresp.py index 1945cbc24..7462e21c4 100644 --- a/control/bench/time_freqresp.py +++ b/control/bench/time_freqresp.py @@ -6,9 +6,11 @@ nstates = 10 sys = rss(nstates) sys_tf = tf(sys) -w = logspace(-1,1,50) +w = logspace(-1, 1, 50) ntimes = 1000 + time_ss = timeit("sys.freqresp(w)", setup="from __main__ import sys, w", number=ntimes) time_tf = timeit("sys_tf.freqresp(w)", setup="from __main__ import sys_tf, w", number=ntimes) + print("State-space model on %d states: %f" % (nstates, time_ss)) print("Transfer-function model on %d states: %f" % (nstates, time_tf)) diff --git a/control/canonical.py b/control/canonical.py index c0244d75f..3f1e30247 100644 --- a/control/canonical.py +++ b/control/canonical.py @@ -11,6 +11,7 @@ __all__ = ['canonical_form', 'reachable_form', 'observable_form'] + def canonical_form(xsys, form='reachable'): """Convert a system into canonical form @@ -40,8 +41,7 @@ def canonical_form(xsys, form='reachable'): elif form == 'modal': return modal_form(xsys) else: - raise ControlNotImplemented( - "Canonical form '%s' not yet implemented" % form) + raise ControlNotImplemented("Canonical form '%s' not yet implemented" % form) # Reachable canonical form @@ -62,8 +62,7 @@ def reachable_form(xsys): """ # Check to make sure we have a SISO system if not issiso(xsys): - raise ControlNotImplemented( - "Canonical forms for MIMO systems not yet supported") + raise ControlNotImplemented("Canonical forms for MIMO systems not yet supported") # Create a new system, starting with a copy of the old one zsys = StateSpace(xsys) @@ -74,9 +73,9 @@ def reachable_form(xsys): zsys.A = zeros(shape(xsys.A)) Apoly = poly(xsys.A) # characteristic polynomial for i in range(0, xsys.states): - zsys.A[0, i] = -Apoly[i+1] / Apoly[0] - if (i+1 < xsys.states): - zsys.A[i+1, i] = 1.0 + zsys.A[0, i] = -Apoly[i + 1] / Apoly[0] + if i + 1 < xsys.states: + zsys.A[i + 1, i] = 1.0 # Compute the reachability matrices for each set of states Wrx = ctrb(xsys.A, xsys.B) @@ -126,9 +125,9 @@ def observable_form(xsys): zsys.A = zeros(shape(xsys.A)) Apoly = poly(xsys.A) # characteristic polynomial for i in range(0, xsys.states): - zsys.A[i, 0] = -Apoly[i+1] / Apoly[0] - if (i+1 < xsys.states): - zsys.A[i, i+1] = 1 + zsys.A[i, 0] = -Apoly[i + 1] / Apoly[0] + if i + 1 < xsys.states: + zsys.A[i, i + 1] = 1 # Compute the observability matrices for each set of states Wrx = obsv(xsys.A, xsys.C) @@ -145,6 +144,7 @@ def observable_form(xsys): return zsys, Tzx + def modal_form(xsys): """Convert a system into modal canonical form @@ -176,7 +176,7 @@ def modal_form(xsys): # Sorting eigenvalues and respective vectors by largest to smallest eigenvalue idx = eigval.argsort()[::-1] eigval = eigval[idx] - eigvec = eigvec[:,idx] + eigvec = eigvec[:, idx] # If all eigenvalues are real, the matrix of eigenvectors is Tzx directly if not iscomplex(eigval).any(): diff --git a/control/config.py b/control/config.py index 10ab8a1ed..62e3fcf35 100644 --- a/control/config.py +++ b/control/config.py @@ -14,6 +14,7 @@ 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 def use_matlab_defaults(): """ @@ -23,9 +24,13 @@ def use_matlab_defaults(): * Bode plots plot gain in dB, phase in degrees, frequency in Hertz """ # Bode plot defaults - global bode_dB; bode_dB = True - global bode_deg; bode_deg = True - global bode_Hz; bode_Hz = True + global bode_dB + bode_dB = True + global bode_deg + bode_deg = True + global bode_Hz + bode_Hz = True + # Set defaults to match FBS (Astrom and Murray) def use_fbs_defaults(): @@ -37,6 +42,9 @@ def use_fbs_defaults(): frequency in Hertz """ # Bode plot defaults - global bode_dB; bode_dB = False - global bode_deg; bode_deg = True - global bode_Hz; bode_Hz = True + global bode_dB + bode_dB = False + global bode_deg + bode_deg = True + global bode_Hz + bode_Hz = True diff --git a/control/ctrlutil.py b/control/ctrlutil.py index 35035c7e9..b0af051e5 100644 --- a/control/ctrlutil.py +++ b/control/ctrlutil.py @@ -47,8 +47,9 @@ __all__ = ['unwrap', 'issys', 'db2mag', 'mag2db'] + # Utility function to unwrap an angle measurement -def unwrap(angle, period=2*math.pi): +def unwrap(angle, period=2 * math.pi): """Unwrap a phase angle to give a continuous curve Parameters @@ -72,15 +73,17 @@ def unwrap(angle, period=2*math.pi): """ dangle = np.diff(angle) - dangle_desired = (dangle + period/2.) % period - period/2. + dangle_desired = (dangle + period / 2.) % period - period / 2. correction = np.cumsum(dangle_desired - dangle) angle[1:] += correction return angle + def issys(obj): """Return True if an object is a system, otherwise False""" return isinstance(obj, lti.LTI) + def db2mag(db): """Convert a gain in decibels (dB) to a magnitude @@ -101,6 +104,7 @@ def db2mag(db): """ return 10. ** (db / 20.) + def mag2db(mag): """Convert a magnitude to decibels (dB) diff --git a/control/delay.py b/control/delay.py index d6350d45b..328f97bee 100644 --- a/control/delay.py +++ b/control/delay.py @@ -1,5 +1,5 @@ # -*-coding: utf-8-*- -#! TODO: add module docstring + # delay.py - functions involving time delays # # Author: Sawyer Fuller @@ -46,6 +46,7 @@ __all__ = ['pade'] + def pade(T, n=1, numdeg=None): """ Create a linear system that approximates a delay. @@ -65,7 +66,7 @@ def pade(T, n=1, numdeg=None): Returns ------- - num, den : array + num, den : list Polynomial coefficients of the delay model, in descending powers of s. Notes @@ -76,6 +77,7 @@ def pade(T, n=1, numdeg=None): 2. M. Vajta, "Some remarks on Padé-approximations", 3rd TEMPUS-INTCOM Symposium """ + if numdeg is None: numdeg = n elif numdeg < 0: @@ -89,26 +91,26 @@ def pade(T, n=1, numdeg=None): raise ValueError("require 0 <= numdeg <= n") if T == 0: - num = [1,] - den = [1,] + num = [1, ] + den = [1, ] else: - num = [0. for i in range(numdeg+1)] + num = [0. for _ in range(numdeg + 1)] num[-1] = 1. cn = 1. - for k in range(1, numdeg+1): + for k in range(1, numdeg + 1): # derived from Gloub and van Loan eq. for Dpq(z) on p. 572 # this accumulative style follows Alg 11.3.1 - cn *= -T * (numdeg - k + 1)/(numdeg + n - k + 1)/k - num[numdeg-k] = cn + cn *= -T * (numdeg - k + 1) / (numdeg + n - k + 1) / k + num[numdeg - k] = cn - den = [0. for i in range(n+1)] + den = [0. for _ in range(n + 1)] den[-1] = 1. cd = 1. - for k in range(1, n+1): + for k in range(1, n + 1): # see cn above - cd *= T * (n - k + 1)/(numdeg + n - k + 1)/k - den[n-k] = cd + cd *= T * (n - k + 1) / (numdeg + n - k + 1) / k + den[n - k] = cd - num = [coeff/den[0] for coeff in num] - den = [coeff/den[0] for coeff in den] + num = [coeff / den[0] for coeff in num] + den = [coeff / den[0] for coeff in den] return num, den diff --git a/control/dtime.py b/control/dtime.py index 36053da33..6bd1fef19 100644 --- a/control/dtime.py +++ b/control/dtime.py @@ -51,6 +51,7 @@ __all__ = ['sample_system', 'c2d'] + # Sample a continuous time system def sample_system(sysc, Ts, method='zoh', alpha=None): """Convert a continuous time system to discrete time @@ -66,6 +67,9 @@ def sample_system(sysc, Ts, method='zoh', alpha=None): Sampling period method : string Method to use for conversion: 'matched', 'tustin', 'zoh' (default) + alpha : float within [0, 1] + The generalized bilinear transformation weighting parameter, which + should only be specified with method="gbt", and is ignored otherwise Returns ------- @@ -91,7 +95,7 @@ def sample_system(sysc, Ts, method='zoh', alpha=None): def c2d(sysc, Ts, method='zoh'): - ''' + """ Return a discrete-time system Parameters @@ -109,7 +113,7 @@ def c2d(sysc, Ts, method='zoh'): 'impulse' Impulse-invariant discretization, currently not implemented 'tustin' Bilinear (Tustin) approximation, only SISO 'matched' Matched pole-zero method, only SISO - ''' + """ # Call the sample_system() function to do the work sysd = sample_system(sysc, Ts, method) if isinstance(sysc, StateSpace) and not isinstance(sysd, StateSpace): diff --git a/control/exception.py b/control/exception.py index 2c4f29704..b6a5f9fe5 100644 --- a/control/exception.py +++ b/control/exception.py @@ -39,32 +39,38 @@ # # $Id$ + class ControlSlycot(Exception): """Exception for Slycot import. Used when we can't import a function from the slycot package""" pass + class ControlDimension(Exception): """Raised when dimensions of system objects are not correct""" pass + class ControlArgument(Exception): """Raised when arguments to a function are not correct""" pass + class ControlMIMONotImplemented(Exception): """Function is not currently implemented for MIMO systems""" pass + class ControlNotImplemented(Exception): """Functionality is not yet implemented""" pass + # Utility function to see if slycot is installed def slycot_check(): try: - import slycot - except: + import slycot # noqa F401 + except ImportError: return False else: return True diff --git a/control/frdata.py b/control/frdata.py index d34200455..24ad2f51e 100644 --- a/control/frdata.py +++ b/control/frdata.py @@ -1,4 +1,6 @@ from __future__ import division + + """ Frequency response data representation and functions. @@ -49,6 +51,7 @@ """ # External function declarations +from warnings import warn import numpy as np from numpy import angle, array, empty, ones, \ real, imag, matrix, absolute, eye, linalg, where, dot @@ -57,6 +60,7 @@ __all__ = ['FRD', 'frd'] + class FRD(LTI): """FRD(d, w) @@ -149,10 +153,10 @@ def __init__(self, *args, **kwargs): dtype=tuple) for i in range(self.fresp.shape[0]): for j in range(self.fresp.shape[1]): - self.ifunc[i,j],u = splprep( + self.ifunc[i, j], u = splprep( u=self.omega, x=[real(self.fresp[i, j, :]), imag(self.fresp[i, j, :])], - w=1.0/(absolute(self.fresp[i, j, :])+0.001), s=0.0) + w=1.0 / (absolute(self.fresp[i, j, :]) + 0.001), s=0.0) else: self.ifunc = None LTI.__init__(self, self.fresp.shape[1], self.fresp.shape[0]) @@ -161,7 +165,7 @@ def __str__(self): """String representation of the transfer function.""" mimo = self.inputs > 1 or self.outputs > 1 - outstr = [ 'frequency response data ' ] + outstr = ['frequency response data '] mt, pt, wt = self.freqresp(self.omega) for i in range(self.inputs): @@ -171,9 +175,8 @@ def __str__(self): outstr.append('Freq [rad/s] Response ') outstr.append('------------ ---------------------') outstr.extend( - [ '%12.3f %10.4g%+10.4gj' % (w, m, p) - for m, p, w in zip(real(self.fresp[j,i,:]), imag(self.fresp[j,i,:]), wt) ]) - + ['%12.3f %10.4g%+10.4gj' % (w, m, p) + for m, p, w in zip(real(self.fresp[j, i, :]), imag(self.fresp[j, i, :]), wt)]) return '\n'.join(outstr) @@ -208,7 +211,7 @@ def __add__(self, other): def __radd__(self, other): """Right add two LTI objects (parallel connection).""" - return self + other; + return self + other def __sub__(self, other): """Subtract two LTI objects.""" @@ -233,18 +236,16 @@ def __mul__(self, other): # Check that the input-output sizes are consistent. if self.inputs != other.outputs: raise ValueError("H = G1*G2: input-output size mismatch" - " G1 has %i input(s), G2 has %i output(s)." % - (self.inputs, other.outputs)) + " G1 has %i input(s), G2 has %i output(s)." % + (self.inputs, other.outputs)) inputs = other.inputs outputs = self.outputs fresp = empty((outputs, inputs, len(self.omega)), dtype=self.fresp.dtype) for i in range(len(self.omega)): - fresp[:,:,i] = dot(self.fresp[:,:,i], other.fresp[:,:,i]) - return FRD(fresp, self.omega, - smooth=(self.ifunc is not None) and - (other.ifunc is not None)) + fresp[:, :, i] = dot(self.fresp[:, :, i], other.fresp[:, :, i]) + return FRD(fresp, self.omega, smooth=(self.ifunc is not None and other.ifunc is not None)) def __rmul__(self, other): """Right Multiply two LTI objects (serial connection).""" @@ -259,39 +260,34 @@ def __rmul__(self, other): # Check that the input-output sizes are consistent. if self.outputs != other.inputs: raise ValueError("H = G1*G2: input-output size mismatch" - " G1 has %i input(s), G2 has %i output(s)." % - (other.inputs, self.outputs)) + " G1 has %i input(s), G2 has %i output(s)." % + (other.inputs, self.outputs)) inputs = self.inputs outputs = other.outputs - fresp = empty((outputs, inputs, len(self.omega)), - dtype=self.fresp.dtype) + fresp = empty((outputs, inputs, len(self.omega)), dtype=self.fresp.dtype) for i in range(len(self.omega)): - fresp[:,:,i] = dot(other.fresp[:,:,i], self.fresp[:,:,i]) - return FRD(fresp, self.omega, - smooth=(self.ifunc is not None) and - (other.ifunc is not None)) + fresp[:, :, i] = dot(other.fresp[:, :, i], self.fresp[:, :, i]) + + return FRD(fresp, self.omega, smooth=(self.ifunc is not None) and (other.ifunc is not None)) # TODO: Division of MIMO transfer function objects is not written yet. def __truediv__(self, other): """Divide two LTI objects.""" if isinstance(other, (int, float, complex, np.number)): - return FRD(self.fresp * (1/other), self.omega, + return FRD(self.fresp * (1 / other), self.omega, smooth=(self.ifunc is not None)) else: other = _convertToFRD(other, omega=self.omega) - - if (self.inputs > 1 or self.outputs > 1 or - other.inputs > 1 or other.outputs > 1): + if self.inputs > 1 or self.outputs > 1 or other.inputs > 1 or other.outputs > 1: raise NotImplementedError( "FRD.__truediv__ is currently implemented only for SISO systems.") - return FRD(self.fresp/other.fresp, self.omega, - smooth=(self.ifunc is not None) and - (other.ifunc is not None)) + return FRD(self.fresp / other.fresp, self.omega, + smooth=(self.ifunc is not None) and (other.ifunc is not None)) # TODO: Remove when transition to python3 complete def __div__(self, other): @@ -306,8 +302,7 @@ def __rtruediv__(self, other): else: other = _convertToFRD(other, omega=self.omega) - if (self.inputs > 1 or self.outputs > 1 or - other.inputs > 1 or other.outputs > 1): + if self.inputs > 1 or self.outputs > 1 or other.inputs > 1 or other.outputs > 1: raise NotImplementedError( "FRD.__rtruediv__ is currently implemented only for SISO systems.") @@ -317,17 +312,17 @@ def __rtruediv__(self, other): def __rdiv__(self, other): return self.__rtruediv__(other) - def __pow__(self,other): + def __pow__(self, other): if not type(other) == int: raise ValueError("Exponent must be an integer") if other == 0: - return FRD(ones(self.fresp.shape),self.omega, - smooth=(self.ifunc is not None)) #unity + return FRD(ones(self.fresp.shape), self.omega, + smooth=(self.ifunc is not None)) # unity if other > 0: - return self * (self**(other-1)) + return self * (self**(other - 1)) if other < 0: return (FRD(ones(self.fresp.shape), self.omega) / self) * \ - (self**(other+1)) + (self**(other + 1)) def evalfr(self, omega): """Evaluate a transfer function at a single angular frequency. @@ -340,8 +335,7 @@ def evalfr(self, omega): intermediate values. """ - warn("FRD.evalfr(omega) will be deprecated in a future release of python-control; use sys.eval(omega) instead", - PendingDeprecationWarning) + warn("FRD.evalfr(omega) will be deprecated in a future release of python-control; use sys.eval(omega) instead", PendingDeprecationWarning) return self._evalfr(omega) # Define the `eval` function to evaluate an FRD at a given (real) @@ -375,21 +369,19 @@ def _evalfr(self, omega): try: out = self.fresp[:, :, where(self.omega == omega)[0][0]] except: - raise ValueError( - "Frequency %f not in frequency list, try an interpolating" - " FRD if you want additional points" % omega) + raise ValueError("Frequency %f not in frequency list, try an interpolating FRD if you want additional points" % omega) else: if getattr(omega, '__iter__', False): for i in range(self.outputs): for j in range(self.inputs): - for k,w in enumerate(omega): - frraw = splev(w, self.ifunc[i,j], der=0) - out[i,j,k] = frraw[0] + 1.0j*frraw[1] + for k, w in enumerate(omega): + frraw = splev(w, self.ifunc[i, j], der=0) + out[i, j, k] = frraw[0] + 1.0j * frraw[1] else: for i in range(self.outputs): for j in range(self.inputs): - frraw = splev(omega, self.ifunc[i,j], der=0) - out[i,j] = frraw[0] + 1.0j*frraw[1] + frraw = splev(omega, self.ifunc[i, j], der=0) + out[i, j] = frraw[0] + 1.0j * frraw[1] return out @@ -424,24 +416,20 @@ def feedback(self, other=1, sign=-1): other = _convertToFRD(other, omega=self.omega) - if (self.outputs != other.inputs or - self.inputs != other.outputs): - raise ValueError( - "FRD.feedback, inputs/outputs mismatch") - fresp = empty((self.outputs, self.inputs, len(other.omega)), - dtype=complex) + if self.outputs != other.inputs or self.inputs != other.outputs: + raise ValueError("FRD.feedback, inputs/outputs mismatch") + + fresp = empty((self.outputs, self.inputs, len(other.omega)), dtype=complex) + # TODO: vectorize this # TODO: handle omega re-mapping for k, w in enumerate(other.omega): - fresp[:, :, k] = self.fresp[:, :, k].view(type=matrix)* \ - linalg.solve( - eye(self.inputs) + - other.fresp[:, :, k].view(type=matrix) * - self.fresp[:, :, k].view(type=matrix), - eye(self.inputs)) + fresp[:, :, k] = self.fresp[:, :, k].view(type=matrix) * \ + linalg.solve(eye(self.inputs) + other.fresp[:, :, k].view(type=matrix) * self.fresp[:, :, k].view(type=matrix), eye(self.inputs)) return FRD(fresp, other.omega, smooth=(self.ifunc is not None)) + def _convertToFRD(sys, omega, inputs=1, outputs=1): """Convert a system to frequency response data form (if needed). @@ -466,8 +454,7 @@ def _convertToFRD(sys, omega, inputs=1, outputs=1): # frequencies match, and system was already frd; simply use return sys - raise NotImplementedError( - "Frequency ranges of FRD do not match, conversion not implemented") + raise NotImplementedError("Frequency ranges of FRD do not match, conversion not implemented") elif isinstance(sys, LTI): omega.sort() @@ -478,25 +465,25 @@ def _convertToFRD(sys, omega, inputs=1, outputs=1): return FRD(fresp, omega, smooth=True) elif isinstance(sys, (int, float, complex, np.number)): - fresp = ones((outputs, inputs, len(omega)), dtype=float)*sys + fresp = ones((outputs, inputs, len(omega)), dtype=float) * sys return FRD(fresp, omega, smooth=True) # try converting constant matrices try: sys = array(sys) - outputs,inputs = sys.shape + outputs, inputs = sys.shape fresp = empty((outputs, inputs, len(omega)), dtype=float) for i in range(outputs): for j in range(inputs): - fresp[i,j,:] = sys[i,j] + fresp[i, j, :] = sys[i, j] return FRD(fresp, omega, smooth=True) except: pass - raise TypeError('''Can't convert given type "%s" to FRD system.''' % - sys.__class__) + raise TypeError('''Can't convert given type "%s" to FRD system.''' % sys.__class__) + -def frd(*args): +def frd(*args, **kwargs): """frd(d, w) Construct a frequency response data model @@ -531,4 +518,4 @@ def frd(*args): -------- FRD, ss, tf """ - return FRD(*args) + return FRD(*args, **kwargs) diff --git a/control/freqplot.py b/control/freqplot.py index 26783d794..15d7ff763 100644 --- a/control/freqplot.py +++ b/control/freqplot.py @@ -90,7 +90,7 @@ def bode_plot(syslist, omega=None, dB=None, Hz=None, deg=None, number of samples margins : boolean if True, plot gain and phase margin - \*args, \**kwargs: + *args, **kwargs: Additional options to matplotlib (color, linestyle, etc) Returns @@ -109,8 +109,8 @@ def bode_plot(syslist, omega=None, dB=None, Hz=None, deg=None, but it returns a MIMO response. 2. If a discrete time model is given, the frequency response is plotted - along the upper branch of the unit circle, using the mapping z = exp(j - \omega dt) where omega ranges from 0 to pi/dt and dt is the discrete + along the upper branch of the unit circle, using the mapping z = exp(j/ omega dt) + where omega ranges from 0 to pi/dt and dt is the discrete timebase. If not timebase is specified (dt = True), dt is set to 1. Examples @@ -140,13 +140,13 @@ def bode_plot(syslist, omega=None, dB=None, Hz=None, deg=None, if Hz: omega_limits *= 2. * math.pi if omega_num: - omega = sp.logspace(np.log10(omega_limits[0]), - np.log10(omega_limits[1]), - num=omega_num, + 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]), + omega = sp.logspace(np.log10(omega_limits[0]), + np.log10(omega_limits[1]), endpoint=True) mags, phases, omegas, nyquistfrqs = [], [], [], [] @@ -243,7 +243,7 @@ def bode_plot(syslist, omega=None, dB=None, Hz=None, deg=None, phase_limit = 180. 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), + ax_phase.axhline(y=phase_limit if deg else math.radians(phase_limit), color='k', linestyle=':') mag_ylim = ax_mag.get_ylim() phase_ylim = ax_phase.get_ylim() @@ -255,39 +255,30 @@ def bode_plot(syslist, omega=None, dB=None, Hz=None, deg=None, ax_mag.loglog([Wcp, Wcp], [1., 1e-8], color='k', linestyle=':') 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=':') + ax_phase.semilogx([Wcp, Wcp], [phase_limit + pm, phase_limit], color='k') 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=':') + ax_phase.semilogx([Wcp, Wcp], [math.radians(phase_limit) + math.radians(pm), math.radians(phase_limit)], color='k') 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], [-20. * np.log10(gm), -1e5], color='k', linestyle=':') ax_mag.semilogx([Wcg, Wcg], [0, -20 * np.log10(gm)], color='k') else: ax_mag.loglog([Wcg, Wcg], [1. / gm, 1e-8], color='k', linestyle=':') ax_mag.loglog([Wcg, Wcg], [1., 1. / gm], color='k') 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=':') 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=':') 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)' % + 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), + 'dB ' if dB else '\b', Wcg, pm if deg else math.radians(pm), 'deg' if deg else 'rad', Wcp)) if nyquistfrq_plot: @@ -342,7 +333,7 @@ def nyquist_plot(syslist, omega=None, Plot=True, color=None, labelFreq=0, *args, Used to specify the color of the plot labelFreq : int Label every nth frequency on the plot - \*args, \**kwargs: + *args, **kwargs: Additional options to matplotlib (color, linestyle, etc) Returns @@ -362,7 +353,10 @@ def nyquist_plot(syslist, omega=None, Plot=True, color=None, labelFreq=0, *args, """ # If argument was a singleton, turn it into a list if not getattr(syslist, '__iter__', False): - syslist = (syslist,) + syslist = [syslist] + + if len(syslist) <= 0: + raise ValueError("syslist should contain at least one element") # Select a default range if none is provided if omega is None: @@ -397,11 +391,11 @@ def nyquist_plot(syslist, omega=None, Plot=True, color=None, labelFreq=0, *args, c = p[0].get_color() ax = plt.gca() # Plot arrow to indicate Nyquist encirclement orientation - ax.arrow(x[0], y[0], (x[1]-x[0])/2, (y[1]-y[0])/2, fc=c, ec=c, + ax.arrow(x[0], y[0], (x[1] - x[0]) / 2, (y[1] - y[0]) / 2, fc=c, ec=c, head_width=0.2, head_length=0.2) plt.plot(x, -y, '-', color=c, *args, **kwargs) - ax.arrow(x[-1], -y[-1], (x[-1]-x[-2])/2, (y[-1]-y[-2])/2, fc=c, ec=c, + ax.arrow(x[-1], -y[-1], (x[-1] - x[-2]) / 2, (y[-1] - y[-2]) / 2, fc=c, ec=c, head_width=0.2, head_length=0.2) # Mark the -1 point @@ -427,8 +421,7 @@ def nyquist_plot(syslist, omega=None, Plot=True, color=None, labelFreq=0, *args, # np.round() is used because 0.99... appears # instead of 1.0, and this would otherwise be # truncated to 0. - plt.text(xpt, ypt, ' ' + str(int(np.round(f / 1000 ** pow1000, 0))) + ' ' + - prefix + 'Hz') + plt.text(xpt, ypt, ' ' + str(int(np.round(f / 1000 ** pow1000, 0))) + ' ' + prefix + 'Hz') if Plot: ax = plt.gca() @@ -521,7 +514,7 @@ def gangof4_plot(P, C, omega=None): # Compute reasonable defaults for axes -def default_frequency_range(syslist, Hz=None, number_of_samples=None, +def default_frequency_range(syslist, Hz=None, number_of_samples=None, feature_periphery_decade=None): """Compute a reasonable default frequency range for frequency domain plots. @@ -598,8 +591,7 @@ def default_frequency_range(syslist, Hz=None, number_of_samples=None, # * 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))] + 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__)) diff --git a/control/grid.py b/control/grid.py index 33fc9e975..ae1d3dede 100644 --- a/control/grid.py +++ b/control/grid.py @@ -7,16 +7,20 @@ from matplotlib.projections import PolarAxes from matplotlib.transforms import Affine2D + class FormatterDMS(object): - '''Transforms angle ticks to damping ratios''' - def __call__(self,direction,factor,values): - angles_deg = values/factor - damping_ratios = np.cos((180-angles_deg)*np.pi/180) - ret = ["%.2f"%val for val in damping_ratios] + """Transforms angle ticks to damping ratios""" + + def __call__(self, direction, factor, values): + angles_deg = values / factor + damping_ratios = np.cos((180 - angles_deg) * np.pi / 180) + ret = ["%.2f" % val for val in damping_ratios] return ret + class ModifiedExtremeFinderCycle(angle_helper.ExtremeFinderCycle): - '''Changed to allow only left hand-side polar grid''' + """Changed to allow only left hand-side polar grid""" + def __call__(self, transform_xy, x1, y1, x2, y2): x_, y_ = np.linspace(x1, x2, self.nx), np.linspace(y1, y2, self.ny) x, y = np.meshgrid(x_, y_) @@ -25,10 +29,10 @@ def __call__(self, transform_xy, x1, y1, x2, y2): with np.errstate(invalid='ignore'): if self.lon_cycle is not None: lon0 = np.nanmin(lon) - lon -= 360. * ((lon - lon0) > 360.) # Changed from 180 to 360 to be able to span only 90-270 (left hand side) + lon -= 360. * ((lon - lon0) > 360.) # Changed from 180 to 360 to be able to span only 90-270 (left hand side) if self.lat_cycle is not None: lat0 = np.nanmin(lat) - lat -= 360. * ((lat - lat0) > 360.) # Changed from 180 to 360 to be able to span only 90-270 (left hand side) + lat -= 360. * ((lat - lat0) > 360.) # Changed from 180 to 360 to be able to span only 90-270 (left hand side) lon_min, lon_max = np.nanmin(lon), np.nanmax(lon) lat_min, lat_max = np.nanmin(lat), np.nanmax(lat) @@ -38,6 +42,7 @@ def __call__(self, transform_xy, x1, y1, x2, y2): return lon_min, lon_max, lat_min, lat_max + def sgrid(): # From matplotlib demos: # https://matplotlib.org/gallery/axisartist/demo_curvelinear_grid.html @@ -45,7 +50,7 @@ def sgrid(): # PolarAxes.PolarTransform takes radian. However, we want our coordinate # system in degree - tr = Affine2D().scale(np.pi/180., 1.) + PolarAxes.PolarTransform() + tr = Affine2D().scale(np.pi / 180., 1.) + PolarAxes.PolarTransform() # polar projection, which involves cycle, and also has limits in # its coordinates, needs a special method to find the extremes # (min, max of the coordinate within the view). @@ -53,10 +58,10 @@ def sgrid(): # 20, 20 : number of sampling points along x, y direction sampling_points = 20 extreme_finder = ModifiedExtremeFinderCycle(sampling_points, sampling_points, - lon_cycle=360, - lat_cycle=None, - lon_minmax=(90,270), - lat_minmax=(0, np.inf),) + lon_cycle=360, + lat_cycle=None, + lon_minmax=(90, 270), + lat_minmax=(0, np.inf),) grid_locator1 = angle_helper.LocatorDMS(15) tick_formatter1 = FormatterDMS() @@ -97,24 +102,25 @@ def sgrid(): fig.add_subplot(ax) - ### RECTANGULAR X Y AXES WITH SCALE - #par2 = ax.twiny() - #par2.axis["top"].toggle(all=False) - #par2.axis["right"].toggle(all=False) - #new_fixed_axis = par2.get_grid_helper().new_fixed_axis - #par2.axis["left"] = new_fixed_axis(loc="left", + # RECTANGULAR X Y AXES WITH SCALE + # par2 = ax.twiny() + # par2.axis["top"].toggle(all=False) + # par2.axis["right"].toggle(all=False) + # new_fixed_axis = par2.get_grid_helper().new_fixed_axis + # par2.axis["left"] = new_fixed_axis(loc="left", # axes=par2, # offset=(0, 0)) - #par2.axis["bottom"] = new_fixed_axis(loc="bottom", + # par2.axis["bottom"] = new_fixed_axis(loc="bottom", # axes=par2, # offset=(0, 0)) - ### FINISH RECTANGULAR + # FINISH RECTANGULAR - ax.grid(True, zorder=0,linestyle='dotted') + ax.grid(True, zorder=0, linestyle='dotted') _final_setup(ax) return ax, fig + def _final_setup(ax): ax.set_xlabel('Real') ax.set_ylabel('Imaginary') @@ -122,6 +128,7 @@ def _final_setup(ax): ax.axvline(x=0, color='black', lw=1) plt.axis('equal') + def nogrid(): f = plt.figure() ax = plt.axes() @@ -129,8 +136,9 @@ def nogrid(): _final_setup(ax) return ax, f + def zgrid(zetas=None, wns=None): - '''Draws discrete damping and frequency grid''' + """Draws discrete damping and frequency grid""" fig = plt.figure() ax = fig.gca() @@ -140,43 +148,42 @@ def zgrid(zetas=None, wns=None): zetas = linspace(0, 0.9, 10) for zeta in zetas: # Calculate in polar coordinates - factor = zeta/sqrt(1-zeta**2) - x = linspace(0, sqrt(1-zeta**2),200) - ang = pi*x - mag = exp(-pi*factor*x) + factor = zeta / sqrt(1 - zeta**2) + x = linspace(0, sqrt(1 - zeta**2), 200) + ang = pi * x + mag = exp(-pi * factor * x) # Draw upper part in retangular coordinates - xret = mag*cos(ang) - yret = mag*sin(ang) - ax.plot(xret,yret, 'k:', lw=1) + xret = mag * cos(ang) + yret = mag * sin(ang) + ax.plot(xret, yret, 'k:', lw=1) # Draw lower part in retangular coordinates - xret = mag*cos(-ang) - yret = mag*sin(-ang) - ax.plot(xret,yret,'k:', lw=1) + xret = mag * cos(-ang) + yret = mag * sin(-ang) + ax.plot(xret, yret, 'k:', lw=1) # Annotation - an_i = int(len(xret)/2.5) + an_i = int(len(xret) / 2.5) an_x = xret[an_i] an_y = yret[an_i] - ax.annotate(str(round(zeta,2)), xy=(an_x, an_y), xytext=(an_x, an_y), size=7) + ax.annotate(str(round(zeta, 2)), xy=(an_x, an_y), xytext=(an_x, an_y), size=7) # Constant natural frequency lines if wns is None: wns = linspace(0, 1, 10) for a in wns: # Calculate in polar coordinates - x = linspace(-pi/2,pi/2,200) - ang = pi*a*sin(x) - mag = exp(-pi*a*cos(x)) + x = linspace(-pi / 2, pi / 2, 200) + ang = pi * a * sin(x) + mag = exp(-pi * a * cos(x)) # Draw in retangular coordinates - xret = mag*cos(ang) - yret = mag*sin(ang) - ax.plot(xret,yret,'k:', lw=1) + xret = mag * cos(ang) + yret = mag * sin(ang) + ax.plot(xret, yret, 'k:', lw=1) # Annotation an_i = -1 an_x = xret[an_i] an_y = yret[an_i] num = '{:1.1f}'.format(a) - ax.annotate("$\\frac{"+num+"\pi}{T}$", xy=(an_x, an_y), xytext=(an_x, an_y), size=9) + ax.annotate("$\\frac{" + num + "\pi}{T}$", xy=(an_x, an_y), xytext=(an_x, an_y), size=9) _final_setup(ax) return ax, fig - diff --git a/control/lti.py b/control/lti.py index 5950d9d58..9f6eb5da3 100644 --- a/control/lti.py +++ b/control/lti.py @@ -18,6 +18,7 @@ __all__ = ['issiso', 'timebase', 'timebaseEqual', 'isdtime', 'isctime', 'pole', 'zero', 'damp', 'evalfr', 'freqresp', 'dcgain'] + class LTI: """LTI is a parent class to linear time-invariant (LTI) system objects. @@ -59,7 +60,7 @@ def isdtime(self, strict=False): """ # If no timebase is given, answer depends on strict flag - if self.dt == None: + if self.dt is None: return True if not strict else False # Look for dt > 0 (also works if dt = True) @@ -71,8 +72,6 @@ def isctime(self, strict=False): Parameters ---------- - sys : LTI system - System to be checked strict: bool (default = False) If strict is True, make sure that timebase is not None """ @@ -82,11 +81,11 @@ def isctime(self, strict=False): return self.dt == 0 def issiso(self): - '''Check to see if a system is single input, single output''' + """Check to see if a system is single input, single output""" return self.inputs == 1 and self.outputs == 1 def damp(self): - '''Natural frequency, damping ratio of system poles + """Natural frequency, damping ratio of system poles Returns ------- @@ -96,21 +95,21 @@ def damp(self): Damping ratio for each system pole poles : array Array of system poles - ''' + """ poles = self.pole() if isdtime(self, strict=True): - splane_poles = np.log(poles)/self.dt + splane_poles = np.log(poles) / self.dt else: splane_poles = poles wn = absolute(splane_poles) - Z = -real(splane_poles)/wn + Z = -real(splane_poles) / wn return wn, Z, poles def dcgain(self): """Return the zero-frequency gain""" - raise NotImplementedError("dcgain not implemented for %s objects" % - str(self.__class__)) + raise NotImplementedError("dcgain not implemented for %s objects" % str(self.__class__)) + # Test to see if a system is SISO def issiso(sys, strict=False): @@ -132,6 +131,7 @@ def issiso(sys, strict=False): # Done with the tricky stuff... return sys.issiso() + # Return the timebase (with conversion if unspecified) def timebase(sys, strict=True): """Return the timebase for an LTI system @@ -148,13 +148,14 @@ def timebase(sys, strict=True): raise ValueError("Timebase not defined") # Return the sample time, with converstion to float if strict is false - if (sys.dt == None): + if sys.dt is None: return None - elif (strict): + elif strict: return float(sys.dt) return sys.dt + # Check to see if two timebases are equal def timebaseEqual(sys1, sys2): """Check to see if two systems have the same timebase @@ -167,15 +168,16 @@ def timebaseEqual(sys1, sys2): timebase (dt > 0) then their timebases must be equal. """ - if (type(sys1.dt) == bool or type(sys2.dt) == bool): + if isinstance(sys1.dt, bool) or isinstance(sys2.dt, bool): # Make sure both are unspecified discrete timebases return type(sys1.dt) == type(sys2.dt) and sys1.dt == sys2.dt - elif (sys1.dt is None or sys2.dt is None): + elif sys1.dt is None or sys2.dt is None: # One or the other is unspecified => the other can be anything return True else: return sys1.dt == sys2.dt + # Check to see if a system is a discrete time system def isdtime(sys, strict=False): """ @@ -201,6 +203,7 @@ def isdtime(sys, strict=False): # Got passed something we don't recognize return False + # Check to see if a system is a continuous time system def isctime(sys, strict=False): """ @@ -226,6 +229,7 @@ def isctime(sys, strict=False): # Got passed something we don't recognize return False + def pole(sys): """ Compute system poles. @@ -285,6 +289,7 @@ def zero(sys): return sys.zero() + def damp(sys, doprint=True): """ Compute natural frequency, damping ratio, and poles of a system @@ -330,7 +335,7 @@ def damp(sys, doprint=True): wn, damping, poles = sys.damp() if doprint: print('_____Eigenvalue______ Damping___ Frequency_') - for p, d, w in zip(poles, damping, wn) : + for p, d, w in zip(poles, damping, wn): if abs(p.imag) < 1e-12: print("%10.4g %10.4g %10.4g" % (p.real, 1.0, -p.real)) @@ -339,6 +344,7 @@ def damp(sys, doprint=True): (p.real, p.imag, d, w)) return wn, damping, poles + def evalfr(sys, x): """ Evaluate the transfer function of an LTI system for a single complex @@ -381,6 +387,7 @@ def evalfr(sys, x): return sys.horner(x)[0][0] return sys.horner(x) + def freqresp(sys, omega): """ Frequency response of an LTI system at multiple angular frequencies. @@ -435,6 +442,7 @@ def freqresp(sys, omega): return sys.freqresp(omega) + def dcgain(sys): """Return the zero-frequency (or DC) gain of the given system diff --git a/control/margins.py b/control/margins.py index 193b6c599..5f756d273 100644 --- a/control/margins.py +++ b/control/margins.py @@ -59,6 +59,7 @@ __all__ = ['stability_margins', 'phase_crossover_frequencies', 'margin'] + # helper functions for stability_margins def _polyimsplit(pol): """split a polynomial with (iw) applied into a real and an @@ -69,12 +70,14 @@ def _polyimsplit(pol): rpencil[-3::-4] = -1. ipencil[-2::-4] = 1. ipencil[-4::-4] = -1. - return pol * rpencil, pol*ipencil + return pol * rpencil, pol * ipencil + def _polysqr(pol): """return a polynomial squared""" return np.polymul(pol, pol) + # Took the framework for the old function by # Sawyer B. Fuller , removed a lot of the innards # and replaced with analytical polynomial functions for LTI systems. @@ -141,8 +144,7 @@ def stability_margins(sysdata, returnall=False, epsw=0.0): sys = sysdata elif getattr(sysdata, '__iter__', False) and len(sysdata) == 3: mag, phase, omega = sysdata - sys = frdata.FRD(mag * np.exp(1j * phase * math.pi/180), - omega, smooth=True) + sys = frdata.FRD(mag * np.exp(1j * phase * math.pi / 180), omega, smooth=True) else: sys = xferfcn._convert_to_transfer_function(sysdata) except Exception as e: @@ -171,9 +173,7 @@ def stability_margins(sysdata, returnall=False, epsw=0.0): # evaluate response at remaining frequencies, to test for phase 180 vs 0 with np.errstate(all='ignore'): - 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)) + 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)) # only keep frequencies where the negative real axis is crossed w_180 = w_180[np.real(resp_w_180) <= 0.0] @@ -192,21 +192,19 @@ def stability_margins(sysdata, returnall=False, epsw=0.0): # point -1, then take the derivative. Second derivative needs to be >0 # to have a minimum test_wstabd = np.polyadd(_polysqr(rden), _polysqr(iden)) - test_wstabn = np.polyadd(_polysqr(np.polyadd(rnum,rden)), - _polysqr(np.polyadd(inum,iden))) + test_wstabn = np.polyadd(_polysqr(np.polyadd(rnum, rden)), + _polysqr(np.polyadd(inum, iden))) test_wstab = np.polysub( - np.polymul(np.polyder(test_wstabn),test_wstabd), - np.polymul(np.polyder(test_wstabd),test_wstabn)) + np.polymul(np.polyder(test_wstabn), test_wstabd), + np.polymul(np.polyder(test_wstabd), test_wstabn)) # find the solutions, for positive omega, and only real ones wstab = np.roots(test_wstab) - wstab = np.real(wstab[(np.imag(wstab) == 0) * - (np.real(wstab) >= 0)]) + wstab = np.real(wstab[(np.imag(wstab) == 0) * (np.real(wstab) >= 0)]) # and find the value of the 2nd derivative there, needs to be positive wstabplus = np.polyval(np.polyder(test_wstab), wstab) - wstab = np.real(wstab[(np.imag(wstab) == 0) * (wstab > epsw) * - (wstabplus > 0.)]) + wstab = np.real(wstab[(np.imag(wstab) == 0) * (wstab > epsw) * (wstabplus > 0.)]) wstab.sort() else: @@ -226,42 +224,37 @@ def dstab(w): # Find all crossings, note that this depends on omega having # a correct range widx = np.where(np.diff(np.sign(mod(sys.omega))))[0] - wc = np.array( - [ sp.optimize.brentq(mod, sys.omega[i], sys.omega[i+1]) - for i in widx if i+1 < len(sys.omega)]) + 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] widx = widx[np.real(sys._evalfr(sys.omega[widx])[0][0]) <= 0] w_180 = np.array( - [ sp.optimize.brentq(arg, sys.omega[i], sys.omega[i+1]) - for i in widx if i+1 < len(sys.omega) ]) + [sp.optimize.brentq(arg, sys.omega[i], sys.omega[i + 1]) + for i in widx if i + 1 < len(sys.omega)]) # find all stab margins? widx = np.where(np.diff(np.sign(np.diff(dstab(sys.omega)))))[0] - wstab = np.array([ sp.optimize.minimize_scalar( - dstab, bracket=(sys.omega[i], sys.omega[i+1])).x - for i in widx if i+1 < len(sys.omega) and - np.diff(np.diff(dstab(sys.omega[i-1:i+2])))[0] > 0 ]) - wstab = wstab[(wstab >= sys.omega[0]) * - (wstab <= sys.omega[-1])] + wstab = np.array([sp.optimize.minimize_scalar(dstab, + bracket=(sys.omega[i], + sys.omega[i + 1])).x + for i in widx if i + 1 < len(sys.omega) and np.diff(np.diff(dstab(sys.omega[i - 1:i + 2])))[0] > 0]) + wstab = wstab[(wstab >= sys.omega[0]) * (wstab <= sys.omega[-1])] - - # margins, as iterables, converted frdata and xferfcn calculations to - # vector for this + # margins, as iterables, converted frdata and xferfcn calculations to vector for this with np.errstate(all='ignore'): gain_w_180 = np.abs(sys._evalfr(w_180)[0][0]) - GM = 1.0/gain_w_180 - SM = np.abs(sys._evalfr(wstab)[0][0]+1) + GM = 1.0 / gain_w_180 + + SM = np.abs(sys._evalfr(wstab)[0][0] + 1) PM = np.remainder(np.angle(sys._evalfr(wc)[0][0], deg=True), 360.0) - 180.0 - + if returnall: return GM, PM, SM, w_180, wc, wstab else: if GM.shape[0] and not np.isinf(GM).all(): with np.errstate(all='ignore'): - gmidx = np.where(np.abs(np.log(GM)) == - np.min(np.abs(np.log(GM)))) + gmidx = np.where(np.abs(np.log(GM)) == np.min(np.abs(np.log(GM)))) else: gmidx = -1 if PM.shape[0]: @@ -272,11 +265,11 @@ def dstab(w): (not SM.shape[0] and float('inf')) or np.amin(SM), (not gmidx != -1 and float('nan')) or w_180[gmidx][0], (not wc.shape[0] and float('nan')) or wc[pmidx][0], - (not wstab.shape[0] and float('nan')) or wstab[SM==np.amin(SM)][0]) + (not wstab.shape[0] and float('nan')) or wstab[SM == np.amin(SM)][0]) # Contributed by Steffen Waldherr -#! TODO - need to add test functions +# TODO - need to add test functions def phase_crossover_frequencies(sys): """Compute frequencies and gains at intersections with real axis in Nyquist plot. @@ -302,14 +295,14 @@ def phase_crossover_frequencies(sys): tf = xferfcn._convert_to_transfer_function(sys) # if not siso, fall back to (0,0) element - #! TODO: should add a check and warning here + # TODO: should add a check and warning here num = tf.num[0][0] den = tf.den[0][0] # Compute frequencies that we cross over the real axis - numj = (1.j)**np.arange(len(num)-1,-1,-1)*num - denj = (-1.j)**np.arange(len(den)-1,-1,-1)*den - allfreq = np.roots(np.imag(np.polymul(numj,denj))) + numj = (1.j)**np.arange(len(num) - 1, -1, -1) * num + denj = (-1.j)**np.arange(len(den) - 1, -1, -1) * den + allfreq = np.roots(np.imag(np.polymul(numj, denj))) realfreq = np.real(allfreq[np.isreal(allfreq)]) realposfreq = realfreq[realfreq >= 0.] @@ -355,16 +348,15 @@ def margin(*args): Examples -------- >>> sys = tf(1, [1, 2, 1, 0]) - >>> gm, pm, wg, wp = margin(sys) + >>> gm, pm, wg, wp = margins(sys) """ if len(args) == 1: sys = args[0] - margin = stability_margins(sys) + margins = stability_margins(sys) elif len(args) == 3: - margin = stability_margins(args) + margins = stability_margins(args) else: - raise ValueError("Margin needs 1 or 3 arguments; received %i." - % len(args)) + raise ValueError("Margin needs 1 or 3 arguments; received %i." % len(args)) - return margin[0], margin[1], margin[3], margin[4] + return margins[0], margins[1], margins[3], margins[4] diff --git a/control/mateqn.py b/control/mateqn.py index e5f7ba4f6..9163e63f9 100644 --- a/control/mateqn.py +++ b/control/mateqn.py @@ -47,9 +47,8 @@ __all__ = ['lyap', 'dlyap', 'dare', 'care'] -#### Lyapunov equation solvers lyap and dlyap -def lyap(A,Q,C=None,E=None): +def lyap(A, Q, C=None, E=None): """ X = lyap(A,Q) solves the continuous-time Lyapunov equation :math:`A X + X A^T + Q = 0` @@ -84,34 +83,33 @@ def lyap(A,Q,C=None,E=None): # Reshape 1-d arrays if len(shape(A)) == 1: - A = A.reshape(1,A.size) + A = A.reshape(1, A.size) if len(shape(Q)) == 1: - Q = Q.reshape(1,Q.size) + Q = Q.reshape(1, Q.size) if C is not None and len(shape(C)) == 1: - C = C.reshape(1,C.size) + C = C.reshape(1, C.size) if E is not None and len(shape(E)) == 1: - E = E.reshape(1,E.size) + E = E.reshape(1, E.size) # Determine main dimensions if size(A) == 1: n = 1 else: - n = size(A,0) + n = size(A, 0) if size(Q) == 1: m = 1 else: - m = size(Q,0) + m = size(Q, 0) # Solve standard Lyapunov equation if C is None and E is None: # Check input data for consistency if shape(A) != shape(Q): - raise ControlArgument("A and Q must be matrices of identical \ - sizes.") + raise ControlArgument("A and Q must be matrices of identical sizes.") if size(A) > 1 and shape(A)[0] != shape(A)[1]: raise ControlArgument("A must be a quadratic matrix.") @@ -124,18 +122,16 @@ def lyap(A,Q,C=None,E=None): # Solve the Lyapunov equation by calling Slycot function sb03md try: - X,scale,sep,ferr,w = sb03md(n,-Q,A,eye(n,n),'C',trana='T') + X, scale, sep, ferr, w = sb03md(n, -Q, A, eye(n, n), 'C', trana='T') except ValueError as ve: if ve.info < 0: e = ValueError(ve.message) e.info = ve.info - elif ve.info == n+1: - e = ValueError("The matrix A and -A have common or very \ - close eigenvalues.") + elif ve.info == n + 1: + e = ValueError("The matrix A and -A have common or very close eigenvalues.") e.info = ve.info else: - e = ValueError("The QR algorithm failed to compute all \ - the eigenvalues (see LAPACK Library routine DGEES).") + e = ValueError("The QR algorithm failed to compute all the eigenvalues (see LAPACK Library routine DGEES).") e.info = ve.info raise e @@ -150,19 +146,19 @@ def lyap(A,Q,C=None,E=None): if (size(C) > 1 and shape(C)[0] != n) or \ (size(C) > 1 and shape(C)[1] != m) or \ - (size(C) == 1 and size(A) != 1) or (size(C) == 1 and size(Q) != 1): + (size(C) == 1 and size(A) != 1) or (size(C) == 1 and size(Q) != 1): raise ControlArgument("C matrix has incompatible dimensions.") # Solve the Sylvester equation by calling the Slycot function sb04md try: - X = sb04md(n,m,A,Q,-C) + X = sb04md(n, m, A, Q, -C) except ValueError as ve: if ve.info < 0: e = ValueError(ve.message) e.info = ve.info elif ve.info > m: e = ValueError("A singular matrix was encountered whilst \ - solving for the %i-th column of matrix X." % ve.info-m) + solving for the %i-th column of matrix X." % ve.info - m) e.info = ve.info else: e = ValueError("The QR algorithm failed to compute all the \ @@ -175,15 +171,13 @@ def lyap(A,Q,C=None,E=None): # Check input data for consistency if (size(Q) > 1 and shape(Q)[0] != shape(Q)[1]) or \ (size(Q) > 1 and shape(Q)[0] != n) or \ - (size(Q) == 1 and n > 1): - raise ControlArgument("Q must be a square matrix with the same \ - dimension as A.") + (size(Q) == 1 and n > 1): + raise ControlArgument("Q must be a square matrix with the same dimension as A.") if (size(E) > 1 and shape(E)[0] != shape(E)[1]) or \ (size(E) > 1 and shape(E)[0] != n) or \ - (size(E) == 1 and n > 1): - raise ControlArgument("E must be a square matrix with the same \ - dimension as A.") + (size(E) == 1 and n > 1): + raise ControlArgument("E must be a square matrix with the same dimension as A.") if not (asarray(Q) == asarray(Q).T).all(): raise ControlArgument("Q must be a symmetric matrix.") @@ -197,8 +191,8 @@ def lyap(A,Q,C=None,E=None): # Solve the generalized Lyapunov equation by calling Slycot # function sg03ad try: - A,E,Q,Z,X,scale,sep,ferr,alphar,alphai,beta = \ - sg03ad('C','B','N','T','L',n,A,E,eye(n,n),eye(n,n),-Q) + A, E, Q, Z, X, scale, sep, ferr, alphar, alphai, beta = \ + sg03ad('C', 'B', 'N', 'T', 'L', n, A, E, eye(n, n), eye(n, n), -Q) except ValueError as ve: if ve.info < 0 or ve.info > 4: e = ValueError(ve.message) @@ -231,7 +225,7 @@ def lyap(A,Q,C=None,E=None): return X -def dlyap(A,Q,C=None,E=None): +def dlyap(A, Q, C=None, E=None): """ dlyap(A,Q) solves the discrete-time Lyapunov equation :math:`A X A^T - X + Q = 0` @@ -271,27 +265,27 @@ def dlyap(A,Q,C=None,E=None): # Reshape 1-d arrays if len(shape(A)) == 1: - A = A.reshape(1,A.size) + A = A.reshape(1, A.size) if len(shape(Q)) == 1: - Q = Q.reshape(1,Q.size) + Q = Q.reshape(1, Q.size) if C is not None and len(shape(C)) == 1: - C = C.reshape(1,C.size) + C = C.reshape(1, C.size) if E is not None and len(shape(E)) == 1: - E = E.reshape(1,E.size) + E = E.reshape(1, E.size) # Determine main dimensions if size(A) == 1: n = 1 else: - n = size(A,0) + n = size(A, 0) if size(Q) == 1: m = 1 else: - m = size(Q,0) + m = size(Q, 0) # Solve standard Lyapunov equation if C is None and E is None: @@ -311,7 +305,7 @@ def dlyap(A,Q,C=None,E=None): # Solve the Lyapunov equation by calling the Slycot function sb03md try: - X,scale,sep,ferr,w = sb03md(n,-Q,A,eye(n,n),'D',trana='T') + X, scale, sep, ferr, w = sb03md(n, -Q, A, eye(n, n), 'D', trana='T') except ValueError as ve: if ve.info < 0: e = ValueError(ve.message) @@ -333,19 +327,19 @@ def dlyap(A,Q,C=None,E=None): if (size(C) > 1 and shape(C)[0] != n) or \ (size(C) > 1 and shape(C)[1] != m) or \ - (size(C) == 1 and size(A) != 1) or (size(C) == 1 and size(Q) != 1): + (size(C) == 1 and size(A) != 1) or (size(C) == 1 and size(Q) != 1): raise ControlArgument("C matrix has incompatible dimensions") # Solve the Sylvester equation by calling Slycot function sb04qd try: - X = sb04qd(n,m,-A,asarray(Q).T,C) + X = sb04qd(n, m, -A, asarray(Q).T, C) except ValueError as ve: if ve.info < 0: e = ValueError(ve.message) e.info = ve.info elif ve.info > m: e = ValueError("A singular matrix was encountered whilst \ - solving for the %i-th column of matrix X." % ve.info-m) + solving for the %i-th column of matrix X." % ve.info - m) e.info = ve.info else: e = ValueError("The QR algorithm failed to compute all the \ @@ -358,13 +352,13 @@ def dlyap(A,Q,C=None,E=None): # Check input data for consistency if (size(Q) > 1 and shape(Q)[0] != shape(Q)[1]) or \ (size(Q) > 1 and shape(Q)[0] != n) or \ - (size(Q) == 1 and n > 1): + (size(Q) == 1 and n > 1): raise ControlArgument("Q must be a square matrix with the same \ dimension as A.") if (size(E) > 1 and shape(E)[0] != shape(E)[1]) or \ (size(E) > 1 and shape(E)[0] != n) or \ - (size(E) == 1 and n > 1): + (size(E) == 1 and n > 1): raise ControlArgument("E must be a square matrix with the same \ dimension as A.") @@ -374,8 +368,8 @@ def dlyap(A,Q,C=None,E=None): # Solve the generalized Lyapunov equation by calling Slycot # function sg03ad try: - A,E,Q,Z,X,scale,sep,ferr,alphar,alphai,beta = \ - sg03ad('D','B','N','T','L',n,A,E,eye(n,n),eye(n,n),-Q) + A, E, Q, Z, X, scale, sep, ferr, alphar, alphai, beta = \ + sg03ad('D', 'B', 'N', 'T', 'L', n, A, E, eye(n, n), eye(n, n), -Q) except ValueError as ve: if ve.info < 0 or ve.info > 4: e = ValueError(ve.message) @@ -408,10 +402,9 @@ def dlyap(A,Q,C=None,E=None): return X +# Riccati equation solvers care and dare -#### Riccati equation solvers care and dare - -def care(A,B,Q,R=None,S=None,E=None): +def care(A, B, Q, R=None, S=None, E=None): """ (X,L,G) = care(A,B,Q,R=None) solves the continuous-time algebraic Riccati equation @@ -453,35 +446,35 @@ def care(A,B,Q,R=None,S=None,E=None): # Reshape 1-d arrays if len(shape(A)) == 1: - A = A.reshape(1,A.size) + A = A.reshape(1, A.size) if len(shape(B)) == 1: - B = B.reshape(1,B.size) + B = B.reshape(1, B.size) if len(shape(Q)) == 1: - Q = Q.reshape(1,Q.size) + Q = Q.reshape(1, Q.size) if R is not None and len(shape(R)) == 1: - R = R.reshape(1,R.size) + R = R.reshape(1, R.size) if S is not None and len(shape(S)) == 1: - S = S.reshape(1,S.size) + S = S.reshape(1, S.size) if E is not None and len(shape(E)) == 1: - E = E.reshape(1,E.size) + E = E.reshape(1, E.size) # Determine main dimensions if size(A) == 1: n = 1 else: - n = size(A,0) + n = size(A, 0) if size(B) == 1: m = 1 else: - m = size(B,1) + m = size(B, 1) if R is None: - R = eye(m,m) + R = eye(m, m) # Solve the standard algebraic Riccati equation if S is None and E is None: @@ -491,12 +484,12 @@ def care(A,B,Q,R=None,S=None,E=None): if (size(Q) > 1 and shape(Q)[0] != shape(Q)[1]) or \ (size(Q) > 1 and shape(Q)[0] != n) or \ - size(Q) == 1 and n > 1: + size(Q) == 1 and n > 1: raise ControlArgument("Q must be a quadratic matrix of the same \ dimension as A.") if (size(B) > 1 and shape(B)[0] != n) or \ - size(B) == 1 and n > 1: + size(B) == 1 and n > 1: raise ControlArgument("Incompatible dimensions of B matrix.") if not (asarray(Q) == asarray(Q).T).all(): @@ -512,12 +505,12 @@ def care(A,B,Q,R=None,S=None,E=None): # Solve the standard algebraic Riccati equation by calling Slycot # functions sb02mt and sb02md try: - A_b,B_b,Q_b,R_b,L_b,ipiv,oufact,G = sb02mt(n,m,B,R) + A_b, B_b, Q_b, R_b, L_b, ipiv, oufact, G = sb02mt(n, m, B, R) except ValueError as ve: if ve.info < 0: e = ValueError(ve.message) e.info = ve.info - elif ve.info == m+1: + elif ve.info == m + 1: e = ValueError("The matrix R is numerically singular.") e.info = ve.info else: @@ -527,7 +520,7 @@ def care(A,B,Q,R=None,S=None,E=None): raise e try: - X,rcond,w,S_o,U,A_inv = sb02md(n,A,G,Q,'C') + X, rcond, w, S_o, U, A_inv = sb02md(n, A, G, Q, 'C') except ValueError as ve: if ve.info < 0 or ve.info > 5: e = ValueError(ve.message) @@ -556,13 +549,13 @@ def care(A,B,Q,R=None,S=None,E=None): # Calculate the gain matrix G if size(R_b) == 1: - G = dot(dot(1/(R_ba), asarray(B_ba).T), X) + G = dot(dot(1 / R_ba, asarray(B_ba).T), X) else: G = dot(solve(R_ba, asarray(B_ba).T), X) # Return the solution X, the closed-loop eigenvalues L and # the gain matrix G - return (X , w[:n] , G ) + return X, w[:n], G # Solve the generalized algebraic Riccati equation elif S is not None and E is not None: @@ -572,30 +565,30 @@ def care(A,B,Q,R=None,S=None,E=None): if (size(Q) > 1 and shape(Q)[0] != shape(Q)[1]) or \ (size(Q) > 1 and shape(Q)[0] != n) or \ - size(Q) == 1 and n > 1: + size(Q) == 1 and n > 1: raise ControlArgument("Q must be a quadratic matrix of the same \ dimension as A.") if (size(B) > 1 and shape(B)[0] != n) or \ - size(B) == 1 and n > 1: + size(B) == 1 and n > 1: raise ControlArgument("Incompatible dimensions of B matrix.") if (size(E) > 1 and shape(E)[0] != shape(E)[1]) or \ (size(E) > 1 and shape(E)[0] != n) or \ - size(E) == 1 and n > 1: + size(E) == 1 and n > 1: raise ControlArgument("E must be a quadratic matrix of the same \ dimension as A.") if (size(R) > 1 and shape(R)[0] != shape(R)[1]) or \ (size(R) > 1 and shape(R)[0] != m) or \ - size(R) == 1 and m > 1: + size(R) == 1 and m > 1: raise ControlArgument("R must be a quadratic matrix of the same \ dimension as the number of columns in the B matrix.") if (size(S) > 1 and shape(S)[0] != n) or \ (size(S) > 1 and shape(S)[1] != m) or \ size(S) == 1 and n > 1 or \ - size(S) == 1 and m > 1: + size(S) == 1 and m > 1: raise ControlArgument("Incompatible dimensions of S matrix.") if not (asarray(Q) == asarray(Q).T).all(): @@ -613,8 +606,8 @@ def care(A,B,Q,R=None,S=None,E=None): # Solve the generalized algebraic Riccati equation by calling the # Slycot function sg02ad try: - rcondu,X,alfar,alfai,beta,S_o,T,U,iwarn = \ - sg02ad('C','B','N','U','N','N','S','R',n,m,0,A,E,B,Q,R,S) + rcondu, X, alfar, alfai, beta, S_o, T, U, iwarn = \ + sg02ad('C', 'B', 'N', 'U', 'N', 'N', 'S', 'R', n, m, 0, A, E, B, Q, R, S) except ValueError as ve: if ve.info < 0 or ve.info > 7: e = ValueError(ve.message) @@ -652,26 +645,27 @@ def care(A,B,Q,R=None,S=None,E=None): raise e # Calculate the closed-loop eigenvalues L - L = zeros((n,1)) + L = zeros((n, 1)) L.dtype = 'complex64' for i in range(n): - L[i] = (alfar[i] + alfai[i]*1j)/beta[i] + L[i] = (alfar[i] + alfai[i] * 1j) / beta[i] # Calculate the gain matrix G if size(R_b) == 1: - G = dot(1/(R_b), dot(asarray(B_b).T, dot(X,E_b)) + asarray(S_b).T) + G = dot(1 / R_b, dot(asarray(B_b).T, dot(X, E_b)) + asarray(S_b).T) else: G = solve(R_b, dot(asarray(B_b).T, dot(X, E_b)) + asarray(S_b).T) # Return the solution X, the closed-loop eigenvalues L and # the gain matrix G - return (X , L , G) + return X, L, G # Invalid set of input parameters else: raise ControlArgument("Invalid set of input parameters.") -def dare(A,B,Q,R,S=None,E=None): + +def dare(A, B, Q, R, S=None, E=None): """ (X,L,G) = dare(A,B,Q,R) solves the discrete-time algebraic Riccati equation @@ -702,7 +696,8 @@ def dare(A,B,Q,R,S=None,E=None): L = eigvals(A - B.dot(G)) return X, L, G -def dare_old(A,B,Q,R,S=None,E=None): + +def dare_old(A, B, Q, R, S=None, E=None): # Make sure we can import required slycot routine try: from slycot import sb02md @@ -722,33 +717,33 @@ def dare_old(A,B,Q,R,S=None,E=None): # Reshape 1-d arrays if len(shape(A)) == 1: - A = A.reshape(1,A.size) + A = A.reshape(1, A.size) if len(shape(B)) == 1: - B = B.reshape(1,B.size) + B = B.reshape(1, B.size) if len(shape(Q)) == 1: - Q = Q.reshape(1,Q.size) + Q = Q.reshape(1, Q.size) if R is not None and len(shape(R)) == 1: - R = R.reshape(1,R.size) + R = R.reshape(1, R.size) if S is not None and len(shape(S)) == 1: - S = S.reshape(1,S.size) + S = S.reshape(1, S.size) if E is not None and len(shape(E)) == 1: - E = E.reshape(1,E.size) + E = E.reshape(1, E.size) # Determine main dimensions if size(A) == 1: n = 1 else: - n = size(A,0) + n = size(A, 0) if size(B) == 1: m = 1 else: - m = size(B,1) + m = size(B, 1) # Solve the standard algebraic Riccati equation if S is None and E is None: @@ -758,12 +753,12 @@ def dare_old(A,B,Q,R,S=None,E=None): if (size(Q) > 1 and shape(Q)[0] != shape(Q)[1]) or \ (size(Q) > 1 and shape(Q)[0] != n) or \ - size(Q) == 1 and n > 1: + size(Q) == 1 and n > 1: raise ControlArgument("Q must be a quadratic matrix of the same \ dimension as A.") if (size(B) > 1 and shape(B)[0] != n) or \ - size(B) == 1 and n > 1: + size(B) == 1 and n > 1: raise ControlArgument("Incompatible dimensions of B matrix.") if not (asarray(Q) == asarray(Q).T).all(): @@ -780,12 +775,12 @@ def dare_old(A,B,Q,R,S=None,E=None): # Solve the standard algebraic Riccati equation by calling Slycot # functions sb02mt and sb02md try: - A_b,B_b,Q_b,R_b,L_b,ipiv,oufact,G = sb02mt(n,m,B,R) + A_b, B_b, Q_b, R_b, L_b, ipiv, oufact, G = sb02mt(n, m, B, R) except ValueError as ve: if ve.info < 0: e = ValueError(ve.message) e.info = ve.info - elif ve.info == m+1: + elif ve.info == m + 1: e = ValueError("The matrix R is numerically singular.") e.info = ve.info else: @@ -795,7 +790,7 @@ def dare_old(A,B,Q,R,S=None,E=None): raise e try: - X,rcond,w,S,U,A_inv = sb02md(n,A,G,Q,'D') + X, rcond, w, S, U, A_inv = sb02md(n, A, G, Q, 'D') except ValueError as ve: if ve.info < 0 or ve.info > 5: e = ValueError(ve.message) @@ -824,15 +819,15 @@ def dare_old(A,B,Q,R,S=None,E=None): # Calculate the gain matrix G if size(R_b) == 1: - G = dot(1/(dot(asarray(B_ba).T, dot(X, B_ba)) + R_ba), \ - dot(asarray(B_ba).T, dot(X, A_ba))) + G = dot(1 / (dot(asarray(B_ba).T, dot(X, B_ba)) + R_ba), + dot(asarray(B_ba).T, dot(X, A_ba))) else: - G = solve(dot(asarray(B_ba).T, dot(X, B_ba)) + R_ba, \ - dot(asarray(B_ba).T, dot(X, A_ba))) + G = solve(dot(asarray(B_ba).T, dot(X, B_ba)) + R_ba, + dot(asarray(B_ba).T, dot(X, A_ba))) # Return the solution X, the closed-loop eigenvalues L and # the gain matrix G - return (X , w[:n] , G) + return X, w[:n], G # Solve the generalized algebraic Riccati equation elif S is not None and E is not None: @@ -842,30 +837,30 @@ def dare_old(A,B,Q,R,S=None,E=None): if (size(Q) > 1 and shape(Q)[0] != shape(Q)[1]) or \ (size(Q) > 1 and shape(Q)[0] != n) or \ - size(Q) == 1 and n > 1: + size(Q) == 1 and n > 1: raise ControlArgument("Q must be a quadratic matrix of the same \ dimension as A.") if (size(B) > 1 and shape(B)[0] != n) or \ - size(B) == 1 and n > 1: + size(B) == 1 and n > 1: raise ControlArgument("Incompatible dimensions of B matrix.") if (size(E) > 1 and shape(E)[0] != shape(E)[1]) or \ (size(E) > 1 and shape(E)[0] != n) or \ - size(E) == 1 and n > 1: + size(E) == 1 and n > 1: raise ControlArgument("E must be a quadratic matrix of the same \ dimension as A.") if (size(R) > 1 and shape(R)[0] != shape(R)[1]) or \ (size(R) > 1 and shape(R)[0] != m) or \ - size(R) == 1 and m > 1: + size(R) == 1 and m > 1: raise ControlArgument("R must be a quadratic matrix of the same \ dimension as the number of columns in the B matrix.") if (size(S) > 1 and shape(S)[0] != n) or \ (size(S) > 1 and shape(S)[1] != m) or \ size(S) == 1 and n > 1 or \ - size(S) == 1 and m > 1: + size(S) == 1 and m > 1: raise ControlArgument("Incompatible dimensions of S matrix.") if not (asarray(Q) == asarray(Q).T).all(): @@ -884,8 +879,8 @@ def dare_old(A,B,Q,R,S=None,E=None): # Solve the generalized algebraic Riccati equation by calling the # Slycot function sg02ad try: - rcondu,X,alfar,alfai,beta,S_o,T,U,iwarn = \ - sg02ad('D','B','N','U','N','N','S','R',n,m,0,A,E,B,Q,R,S) + rcondu, X, alfar, alfai, beta, S_o, T, U, iwarn = \ + sg02ad('D', 'B', 'N', 'U', 'N', 'N', 'S', 'R', n, m, 0, A, E, B, Q, R, S) except ValueError as ve: if ve.info < 0 or ve.info > 7: e = ValueError(ve.message) @@ -922,22 +917,22 @@ def dare_old(A,B,Q,R,S=None,E=None): e.info = ve.info raise e - L = zeros((n,1)) + L = zeros((n, 1)) L.dtype = 'complex64' for i in range(n): - L[i] = (alfar[i] + alfai[i]*1j)/beta[i] + L[i] = (alfar[i] + alfai[i] * 1j) / beta[i] # Calculate the gain matrix G if size(R_b) == 1: - G = dot(1/(dot(asarray(B_b).T, dot(X,B_b)) + R_b), \ - dot(asarray(B_b).T, dot(X,A_b)) + asarray(S_b).T) + G = dot(1 / (dot(asarray(B_b).T, dot(X, B_b)) + R_b), + dot(asarray(B_b).T, dot(X, A_b)) + asarray(S_b).T) else: - G = solve(dot(asarray(B_b).T, dot(X,B_b)) + R_b, \ - dot(asarray(B_b).T, dot(X,A_b)) + asarray(S_b).T) + G = solve(dot(asarray(B_b).T, dot(X, B_b)) + R_b, + dot(asarray(B_b).T, dot(X, A_b)) + asarray(S_b).T) # Return the solution X, the closed-loop eigenvalues L and # the gain matrix G - return (X , L , G) + return X, L, G # Invalid set of input parameters else: diff --git a/control/matlab/__init__.py b/control/matlab/__init__.py index 3c9111b3b..2244188a7 100644 --- a/control/matlab/__init__.py +++ b/control/matlab/__init__.py @@ -51,8 +51,8 @@ """ # Import MATLAB-like functions that are defined in other packages -from scipy.signal import zpk2ss, ss2zpk, tf2zpk, zpk2tf -from numpy import linspace, logspace +from scipy.signal import zpk2ss, ss2zpk, tf2zpk, zpk2tf # noqa F403, F401 +from numpy import linspace, logspace # noqa F403, F401 # If configuration is not yet set, import and use MATLAB defaults import sys @@ -61,30 +61,30 @@ config.use_matlab_defaults() # Control system library -from ..statesp import * -from ..xferfcn import * -from ..lti import * -from ..frdata import * -from ..dtime import * -from ..exception import ControlArgument +from ..statesp import * # noqa F403, F401 +from ..xferfcn import * # noqa F403, F401 +from ..lti import * # noqa F403, F401 +from ..frdata import * # noqa F403, F401 +from ..dtime import * # noqa F403, F401 +from ..exception import ControlArgument # noqa F403, F401 # Import MATLAB-like functions that can be used as-is -from ..ctrlutil import * -from ..freqplot import nyquist, gangof4 -from ..nichols import nichols -from ..bdalg import * -from ..pzmap import * -from ..statefbk import * -from ..delay import * -from ..modelsimp import * -from ..mateqn import * -from ..margins import margin -from ..rlocus import rlocus -from ..dtime import c2d +from ..ctrlutil import * # noqa F403, F401 +from ..freqplot import nyquist, gangof4 # noqa F403, F401 +from ..nichols import nichols # noqa F403, F401 +from ..bdalg import * # noqa F403, F401 +from ..pzmap import * # noqa F403, F401 +from ..statefbk import * # noqa F403, F401 +from ..delay import * # noqa F403, F401 +from ..modelsimp import * # noqa F403, F401 +from ..mateqn import * # noqa F403, F401 +from ..margins import margin # noqa F403, F401 +from ..rlocus import rlocus # noqa F403, F401 +from ..dtime import c2d # noqa F403, F401 # Import functions specific to Matlab compatibility package -from .timeresp import * -from .wrappers import * +from .timeresp import * # noqa F403, F401 +from .wrappers import * # noqa F403, F401 r""" The following tables give an overview of the module ``control.matlab``. diff --git a/control/matlab/timeresp.py b/control/matlab/timeresp.py index ffdf9dd6f..3ee3ba2aa 100644 --- a/control/matlab/timeresp.py +++ b/control/matlab/timeresp.py @@ -6,8 +6,9 @@ __all__ = ['step', 'impulse', 'initial', 'lsim'] + def step(sys, T=None, X0=0., input=0, output=None, return_x=False): - ''' + """ Step response of a linear system If the system has multiple inputs or outputs (MIMO), one input has @@ -55,19 +56,20 @@ def step(sys, T=None, X0=0., input=0, output=None, return_x=False): Examples -------- >>> yout, T = step(sys, T, X0) - ''' + """ from ..timeresp import step_response T, yout, xout = step_response(sys, T, X0, input, output, - transpose = True, return_x=True) + transpose=True, return_x=True) if return_x: return yout, T, xout return yout, T + def impulse(sys, T=None, X0=0., input=0, output=None, return_x=False): - ''' + """ Impulse response of a linear system If the system has multiple inputs or outputs (MIMO), one input has @@ -113,18 +115,19 @@ def impulse(sys, T=None, X0=0., input=0, output=None, return_x=False): Examples -------- >>> yout, T = impulse(sys, T) - ''' + """ from ..timeresp import impulse_response T, yout, xout = impulse_response(sys, T, X0, input, output, - transpose = True, return_x=True) + transpose=True, return_x=True) if return_x: return yout, T, xout return yout, T + def initial(sys, T=None, X0=0., input=None, output=None, return_x=False): - ''' + """ Initial condition response of a linear system If the system has multiple outputs (?IMO), optionally, one output @@ -170,7 +173,7 @@ def initial(sys, T=None, X0=0., input=None, output=None, return_x=False): -------- >>> yout, T = initial(sys, T, X0) - ''' + """ from ..timeresp import initial_response T, yout, xout = initial_response(sys, T, X0, output=output, transpose=True, return_x=True) @@ -180,8 +183,9 @@ def initial(sys, T=None, X0=0., input=None, output=None, return_x=False): return yout, T + def lsim(sys, U=0., T=None, X0=0.): - ''' + """ Simulate the output of a linear system. As a convenience for parameters `U`, `X0`: @@ -222,7 +226,7 @@ def lsim(sys, U=0., T=None, X0=0.): Examples -------- >>> yout, T, xout = lsim(sys, U, T, X0) - ''' + """ from ..timeresp import forced_response - T, yout, xout = forced_response(sys, T, U, X0, transpose = True) + T, yout, xout = forced_response(sys, T, U, X0, transpose=True) return yout, T, xout diff --git a/control/matlab/wrappers.py b/control/matlab/wrappers.py index d83890a33..0271240d7 100644 --- a/control/matlab/wrappers.py +++ b/control/matlab/wrappers.py @@ -3,12 +3,15 @@ """ import numpy as np + +from control import ControlArgument from ..statesp import ss from ..xferfcn import tf from scipy.signal import zpk2tf __all__ = ['bode', 'ngrid', 'dcgain'] + def bode(*args, **keywords): """bode(syslist[, omega, dB, Hz, deg, ...]) @@ -56,22 +59,25 @@ def bode(*args, **keywords): # If the first argument is a list, then assume python-control calling format from ..freqplot import bode as bode_orig - if (getattr(args[0], '__iter__', False)): + if getattr(args[0], '__iter__', False): return bode_orig(*args, **keywords) # Otherwise, run through the arguments and collect up arguments - syslist = []; plotstyle=[]; omega=None; - i = 0; + syslist = [] + plotstyle = [] + omega = None + i = 0 + while i < len(args): # Check to see if this is a system of some sort from ..ctrlutil import issys - if (issys(args[i])): + if issys(args[i]): # Append the system to our list of systems syslist.append(args[i]) i += 1 # See if the next object is a plotsytle (string) - if (i < len(args) and isinstance(args[i], str)): + if i < len(args) and isinstance(args[i], str): plotstyle.append(args[i]) i += 1 @@ -79,7 +85,7 @@ def bode(*args, **keywords): continue # See if this is a frequency list - elif (isinstance(args[i], (list, np.ndarray))): + elif isinstance(args[i], (list, np.ndarray)): omega = args[i] i += 1 break @@ -88,29 +94,35 @@ def bode(*args, **keywords): raise ControlArgument("unrecognized argument type") # Check to make sure that we processed all arguments - if (i < len(args)): + if i < len(args): raise ControlArgument("not all arguments processed") # Check to make sure we got the same number of plotstyles as systems - if (len(plotstyle) != 0 and len(syslist) != len(plotstyle)): + if len(plotstyle) != 0 and len(syslist) != len(plotstyle): raise ControlArgument("number of systems and plotstyles should be equal") # Warn about unimplemented plotstyles - #! TODO: remove this when plot styles are implemented in bode() - #! TODO: uncomment unit test code that tests this out - if (len(plotstyle) != 0): - print("Warning (matlab.bode): plot styles not implemented"); + # TODO: remove this when plot styles are implemented in bode() + # TODO: uncomment unit test code that tests this out + if len(plotstyle) != 0: + print("Warning (matlab.bode): plot styles not implemented") # Call the bode command return bode_orig(syslist, omega, **keywords) + from ..nichols import nichols_grid + + def ngrid(): return nichols_grid() + + ngrid.__doc__ = nichols_grid.__doc__ + def dcgain(*args): - ''' + """ Compute the gain of the system in steady state. The function takes either 1, 2, 3, or 4 parameters: @@ -141,8 +153,8 @@ def dcgain(*args): computes: .. math:: gain = - C \cdot A^{-1} \cdot B + D - ''' - #Convert the parameters to state space form + """ + # Convert the parameters to state space form if len(args) == 4: A, B, C, D = args return ss(A, B, C, D).dcgain() @@ -157,5 +169,4 @@ def dcgain(*args): sys, = args return sys.dcgain() else: - raise ValueError("Function ``dcgain`` needs either 1, 2, 3 or 4 " - "arguments.") + raise ValueError("Function ``dcgain`` needs either 1, 2, 3 or 4 arguments.") diff --git a/control/modelsimp.py b/control/modelsimp.py index 46f468685..d163e89ad 100644 --- a/control/modelsimp.py +++ b/control/modelsimp.py @@ -52,10 +52,12 @@ __all__ = ['hsvd', 'balred', 'modred', 'era', 'markov', 'minreal'] + # Hankel Singular Value Decomposition -# The following returns the Hankel singular values, which are singular values -#of the matrix formed by multiplying the controllability and observability -#grammians +# +# The following returns the Hankel singular values, which are singular values +# of the matrix formed by multiplying the controllability and observability +# grammians def hsvd(sys): """Calculate the Hankel singular values. @@ -87,11 +89,11 @@ def hsvd(sys): """ # TODO: implement for discrete time systems - if (isdtime(sys, strict=True)): + if isdtime(sys, strict=True): raise NotImplementedError("Function not implemented in discrete time") - Wc = gram(sys,'c') - Wo = gram(sys,'o') + Wc = gram(sys, 'c') + Wo = gram(sys, 'o') WoWc = np.dot(Wo, Wc) w, v = np.linalg.eig(WoWc) @@ -102,6 +104,7 @@ def hsvd(sys): # Return the Hankel singular values return hsv + def modred(sys, ELIM, method='matchdc'): """ Model reduction of `sys` by eliminating the states in `ELIM` using a given @@ -134,21 +137,17 @@ def modred(sys, ELIM, method='matchdc'): >>> rsys = modred(sys, ELIM, method='truncate') """ - #Check for ss system object, need a utility for this? - - #TODO: Check for continous or discrete, only continuous supported right now - # if isCont(): - # dico = 'C' - # elif isDisc(): - # dico = 'D' - # else: - if (isctime(sys)): - dico = 'C' - else: + # TODO: Check for ss system object, need a utility for this? + # TODO: Check for continous or discrete, only continuous supported right now + # if isCont(): + # dico = 'C' + # elif isDisc(): + # dico = 'D' + # else: + if not isctime(sys): raise NotImplementedError("Function not implemented in discrete time") - - #Check system is stable + # Check system is stable if np.any(np.linalg.eigvals(sys.A).real >= 0.0): raise ValueError("Oops, the system is unstable!") @@ -156,24 +155,24 @@ def modred(sys, ELIM, method='matchdc'): # Create list of elements not to eliminate (NELIM) NELIM = [i for i in range(len(sys.A)) if i not in ELIM] # A1 is a matrix of all columns of sys.A not to eliminate - A1 = sys.A[:,NELIM[0]] + A1 = sys.A[:, NELIM[0]] for i in NELIM[1:]: - A1 = np.hstack((A1, sys.A[:,i])) - A11 = A1[NELIM,:] - A21 = A1[ELIM,:] + A1 = np.hstack((A1, sys.A[:, i])) + A11 = A1[NELIM, :] + A21 = A1[ELIM, :] # A2 is a matrix of all columns of sys.A to eliminate - A2 = sys.A[:,ELIM[0]] + A2 = sys.A[:, ELIM[0]] for i in ELIM[1:]: - A2 = np.hstack((A2, sys.A[:,i])) - A12 = A2[NELIM,:] - A22 = A2[ELIM,:] + A2 = np.hstack((A2, sys.A[:, i])) + A12 = A2[NELIM, :] + A22 = A2[ELIM, :] - C1 = sys.C[:,NELIM] - C2 = sys.C[:,ELIM] - B1 = sys.B[NELIM,:] - B2 = sys.B[ELIM,:] + C1 = sys.C[:, NELIM] + C2 = sys.C[:, ELIM] + B1 = sys.B[NELIM, :] + B2 = sys.B[ELIM, :] - if method=='matchdc': + if method == 'matchdc': # if matchdc, residualize # Check if the matrix A22 is invertible @@ -189,11 +188,11 @@ def modred(sys, ELIM, method='matchdc'): A22I_A21 = A22I_A21_B2[:, :A21.shape[1]] A22I_B2 = A22I_A21_B2[:, A21.shape[1]:] - Ar = A11 - A12*A22I_A21 - Br = B1 - A12*A22I_B2 - Cr = C1 - C2*A22I_A21 - Dr = sys.D - C2*A22I_B2 - elif method=='truncate': + Ar = A11 - A12 * A22I_A21 + Br = B1 - A12 * A22I_B2 + Cr = C1 - C2 * A22I_A21 + Dr = sys.D - C2 * A22I_B2 + elif method == 'truncate': # if truncate, simply discard state x2 Ar = A11 Br = B1 @@ -202,9 +201,10 @@ def modred(sys, ELIM, method='matchdc'): else: raise ValueError("Oops, method is not supported!") - rsys = StateSpace(Ar,Br,Cr,Dr) + rsys = StateSpace(Ar, Br, Cr, Dr) return rsys + def balred(sys, orders, method='truncate', alpha=None): """ Balanced reduced order model of sys of a given order. @@ -221,7 +221,7 @@ def balred(sys, orders, method='truncate', alpha=None): ---------- sys: StateSpace Original system to reduce - orders: integer or array of integer + orders: integer or iterable of integer Desired order of reduced order model (if a vector, returns a vector of systems) method: string @@ -254,22 +254,23 @@ def balred(sys, orders, method='truncate', alpha=None): >>> rsys = balred(sys, orders, method='truncate') """ - if method!='truncate' and method!='matchdc': - raise ValueError("supported methods are 'truncate' or 'matchdc'") - elif method=='truncate': + + if method == 'truncate': try: from slycot import ab09md, ab09ad except ImportError: raise ControlSlycot("can't find slycot subroutine ab09md or ab09ad") - elif method=='matchdc': + elif method == 'matchdc': try: from slycot import ab09nd except ImportError: raise ControlSlycot("can't find slycot subroutine ab09nd") + else: + raise ValueError("supported methods are 'truncate' or 'matchdc'") - #Check for ss system object, need a utility for this? + # Check for ss system object, need a utility for this? - #TODO: Check for continous or discrete, only continuous supported right now + # TODO: Check for continous or discrete, only continuous supported right now # if isCont(): # dico = 'C' # elif isDisc(): @@ -277,7 +278,7 @@ def balred(sys, orders, method='truncate', alpha=None): # else: dico = 'C' - job = 'B' # balanced (B) or not (N) + job = 'B' # balanced (B) or not (N) equil = 'N' # scale (S) or not (N) if alpha is None: if dico == 'C': @@ -285,41 +286,42 @@ def balred(sys, orders, method='truncate', alpha=None): elif dico == 'D': alpha = 1. - rsys = [] #empty list for reduced systems + rsys = [] # empty list for reduced systems - #check if orders is a list or a scalar + # check if orders is a list or a scalar try: - order = iter(orders) - except TypeError: #if orders is a scalar + _ = iter(orders) + except TypeError: # if orders is a scalar orders = [orders] for i in orders: - n = np.size(sys.A,0) - m = np.size(sys.B,1) - p = np.size(sys.C,0) + n = np.size(sys.A, 0) + m = np.size(sys.B, 1) + p = np.size(sys.C, 0) if method == 'truncate': - #check system stability + # check system stability if np.any(np.linalg.eigvals(sys.A).real >= 0.0): - #unstable branch - Nr, Ar, Br, Cr, Ns, hsv = ab09md(dico,job,equil,n,m,p,sys.A,sys.B,sys.C,alpha=alpha,nr=i,tol=0.0) + # unstable branch + Nr, Ar, Br, Cr, Ns, hsv = ab09md(dico, job, equil, n, m, p, sys.A, sys.B, sys.C, alpha=alpha, nr=i, tol=0.0) else: - #stable branch - Nr, Ar, Br, Cr, hsv = ab09ad(dico,job,equil,n,m,p,sys.A,sys.B,sys.C,nr=i,tol=0.0) + # stable branch + Nr, Ar, Br, Cr, hsv = ab09ad(dico, job, equil, n, m, p, sys.A, sys.B, sys.C, nr=i, tol=0.0) rsys.append(StateSpace(Ar, Br, Cr, sys.D)) elif method == 'matchdc': - Nr, Ar, Br, Cr, Dr, Ns, hsv = ab09nd(dico,job,equil,n,m,p,sys.A,sys.B,sys.C,sys.D,alpha=alpha,nr=i,tol1=0.0,tol2=0.0) + Nr, Ar, Br, Cr, Dr, Ns, hsv = ab09nd(dico, job, equil, n, m, p, sys.A, sys.B, sys.C, sys.D, alpha=alpha, nr=i, tol1=0.0, tol2=0.0) rsys.append(StateSpace(Ar, Br, Cr, Dr)) - #if orders was a scalar, just return the single reduced model, not a list + # if orders was a scalar, just return the single reduced model, not a list if len(orders) == 1: return rsys[0] - #if orders was a list/vector, return a list/vector of systems + # if orders was a list/vector, return a list/vector of systems else: return rsys + def minreal(sys, tol=None, verbose=True): - ''' + """ Eliminates uncontrollable or unobservable states in state-space models or cancelling pole-zero pairs in transfer functions. The output sysr has minimal order and the same response @@ -338,13 +340,14 @@ def minreal(sys, tol=None, verbose=True): ------- rsys: StateSpace or TransferFunction Cleaned model - ''' + """ sysr = sys.minreal(tol) if verbose: print("{nstates} states have been removed from the model".format( - nstates=len(sys.pole()) - len(sysr.pole()))) + nstates=len(sys.pole()) - len(sysr.pole()))) return sysr + def era(YY, m, n, nin, nout, r): """ Calculate an ERA model of order `r` based on the impulse-response data `YY`. @@ -377,6 +380,7 @@ def era(YY, m, n, nin, nout, r): """ raise NotImplementedError('This function is not implemented yet.') + def markov(Y, U, M): """ Calculate the first `M` Markov parameters [D CB CAB ...] @@ -412,12 +416,12 @@ def markov(Y, U, M): # Construct a matrix of control inputs to invert UU = Umat - for i in range(1, M-1): - newCol = np.vstack((0, UU[0:n-1,i-2])) + for i in range(1, M - 1): + newCol = np.vstack((0, UU[0:n - 1, i - 2])) UU = np.hstack((UU, newCol)) - Ulast = np.vstack((0, UU[0:n-1,M-2])) - for i in range(n-1,0,-1): - Ulast[i] = np.sum(Ulast[0:i-1]) + Ulast = np.vstack((0, UU[0:n - 1, M - 2])) + for i in range(n - 1, 0, -1): + Ulast[i] = np.sum(Ulast[0:i - 1]) UU = np.hstack((UU, Ulast)) # Invert and solve for Markov parameters diff --git a/control/nichols.py b/control/nichols.py index a4c5795a4..b0c6484be 100644 --- a/control/nichols.py +++ b/control/nichols.py @@ -94,7 +94,7 @@ def nichols_plot(sys_list, omega=None, grid=True): # Convert to Nichols-plot format (phase in degrees, # and magnitude in dB) x = unwrap(sp.degrees(phase), 360) - y = 20*sp.log10(mag) + y = 20 * sp.log10(mag) # Generate the plot plt.plot(x, y) @@ -154,8 +154,7 @@ def nichols_grid(cl_mags=None, cl_phases=None, line_style='dotted'): # magnitudes are approximately aligned with open-loop magnitudes beyond # the value of np.min(key_cl_mags) cl_mag_step = -20.0 # dB - extended_cl_mags = np.arange(np.min(key_cl_mags), - ol_mag_min + cl_mag_step, cl_mag_step) + extended_cl_mags = np.arange(np.min(key_cl_mags), ol_mag_min + cl_mag_step, cl_mag_step) cl_mags = np.concatenate((extended_cl_mags, key_cl_mags)) # N-circle phases (should be in the range -360 to 0) @@ -173,12 +172,12 @@ def nichols_grid(cl_mags=None, cl_phases=None, line_style='dotted'): # Find the M-contours m = m_circles(cl_mags, phase_min=np.min(cl_phases), phase_max=np.max(cl_phases)) - m_mag = 20*sp.log10(np.abs(m)) + m_mag = 20 * sp.log10(np.abs(m)) m_phase = sp.mod(sp.degrees(sp.angle(m)), -360.0) # Unwrap # Find the N-contours n = n_circles(cl_phases, mag_min=np.min(cl_mags), mag_max=np.max(cl_mags)) - n_mag = 20*sp.log10(np.abs(n)) + n_mag = 20 * sp.log10(np.abs(n)) n_phase = sp.mod(sp.degrees(sp.angle(n)), -360.0) # Unwrap # Plot the contours behind other plot elements. @@ -187,16 +186,14 @@ def nichols_grid(cl_mags=None, cl_phases=None, line_style='dotted'): # over the range -360 < phase < 0. Given the range # the base chart is computed over, the phase offset should be 0 # for -360 < ol_phase_min < 0. - phase_offset_min = 360.0*np.ceil(ol_phase_min/360.0) - phase_offset_max = 360.0*np.ceil(ol_phase_max/360.0) + 360.0 + phase_offset_min = 360.0 * np.ceil(ol_phase_min / 360.0) + phase_offset_max = 360.0 * np.ceil(ol_phase_max / 360.0) + 360.0 phase_offsets = np.arange(phase_offset_min, phase_offset_max, 360.0) for phase_offset in phase_offsets: # Draw M and N contours - plt.plot(m_phase + phase_offset, m_mag, color='lightgray', - linestyle=line_style, zorder=0) - plt.plot(n_phase + phase_offset, n_mag, color='lightgray', - linestyle=line_style, zorder=0) + plt.plot(m_phase + phase_offset, m_mag, color='lightgray', linestyle=line_style, zorder=0) + plt.plot(n_phase + phase_offset, n_mag, color='lightgray', linestyle=line_style, zorder=0) # Add magnitude labels for x, y, m in zip(m_phase[:][-1] + phase_offset, m_mag[:][-1], cl_mags): @@ -235,10 +232,10 @@ def closed_loop_contours(Gcl_mags, Gcl_phases): # Compute the contours in Gcl-space. Since we're given closed-loop # magnitudes and phases, this is just a case of converting them into # a complex number. - Gcl = Gcl_mags*sp.exp(1.j*Gcl_phases) + Gcl = Gcl_mags * sp.exp(1.j * Gcl_phases) # Invert Gcl = Gol/(1+Gol) to map the contours into the open-loop space - return Gcl/(1.0 - Gcl) + return Gcl / (1.0 - Gcl) def m_circles(mags, phase_min=-359.75, phase_max=-0.25): @@ -263,7 +260,7 @@ def m_circles(mags, phase_min=-359.75, phase_max=-0.25): # Convert magnitudes and phase range into a grid suitable for # building contours phases = sp.radians(sp.linspace(phase_min, phase_max, 2000)) - Gcl_mags, Gcl_phases = sp.meshgrid(10.0**(mags/20.0), phases) + Gcl_mags, Gcl_phases = sp.meshgrid(10.0**(mags / 20.0), phases) return closed_loop_contours(Gcl_mags, Gcl_phases) @@ -288,7 +285,7 @@ def n_circles(phases, mag_min=-40.0, mag_max=12.0): """ # Convert phases and magnitude range into a grid suitable for # building contours - mags = sp.linspace(10**(mag_min/20.0), 10**(mag_max/20.0), 2000) + mags = sp.linspace(10**(mag_min / 20.0), 10**(mag_max / 20.0), 2000) Gcl_phases, Gcl_mags = sp.meshgrid(sp.radians(phases), mags) return closed_loop_contours(Gcl_mags, Gcl_phases) diff --git a/control/phaseplot.py b/control/phaseplot.py index 10fbec640..e68a98e80 100644 --- a/control/phaseplot.py +++ b/control/phaseplot.py @@ -41,7 +41,6 @@ import matplotlib.pyplot as mpl from scipy.integrate import odeint -from .exception import ControlNotImplemented __all__ = ['phase_plot', 'box_grid'] @@ -54,8 +53,8 @@ def _find(condition): def phase_plot(odefun, X=None, Y=None, scale=1, X0=None, T=None, - lingrid=None, lintime=None, logtime=None, timepts=None, - parms=(), verbose=True): + lingrid=None, lintime=None, logtime=None, timepts=None, + parms=(), verbose=True): """ Phase plot for 2D dynamical systems @@ -127,31 +126,33 @@ def phase_plot(odefun, X=None, Y=None, scale=1, X0=None, T=None, # # Figure out ranges for phase plot (argument processing) # - #! TODO: need to add error checking to arguments - #! TODO: think through proper action if multiple options are given + # TODO: need to add error checking to arguments + # TODO: think through proper action if multiple options are given # - autoFlag = False; logtimeFlag = False; timeptsFlag = False; Narrows = 0; + autoFlag = False + logtimeFlag = False + timeptsFlag = False if lingrid is not None: - autoFlag = True; - Narrows = lingrid; - if (verbose): + autoFlag = True + Narrows = lingrid + if verbose: print('Using auto arrows\n') elif logtime is not None: - logtimeFlag = True; - Narrows = logtime[0]; - timefactor = logtime[1]; - if (verbose): + logtimeFlag = True + Narrows = logtime[0] + timefactor = logtime[1] + if verbose: print('Using logtime arrows\n') elif timepts is not None: - timeptsFlag = True; - Narrows = len(timepts); + timeptsFlag = True + Narrows = len(timepts) # Figure out the set of points for the quiver plot - #! TODO: Add sanity checks - elif (X is not None and Y is not None): + # TODO: Add sanity checks + elif X is not None and Y is not None: (x1, x2) = np.meshgrid( np.linspace(X[0], X[1], X[2]), np.linspace(Y[0], Y[1], Y[2])) @@ -159,43 +160,43 @@ def phase_plot(odefun, X=None, Y=None, scale=1, X0=None, T=None, else: # If we weren't given any grid points, don't plot arrows - Narrows = 0; + Narrows = 0 - if ((not autoFlag) and (not logtimeFlag) and (not timeptsFlag) - and (Narrows > 0)): + if not autoFlag and not logtimeFlag and not timeptsFlag and Narrows > 0: # Now calculate the vector field at those points - (nr,nc) = x1.shape; + (nr, nc) = x1.shape dx = np.empty((nr, nc, 2)) for i in range(nr): for j in range(nc): - dx[i, j, :] = np.squeeze(odefun((x1[i,j], x2[i,j]), 0, *parms)) + dx[i, j, :] = np.squeeze(odefun((x1[i, j], x2[i, j]), 0, *parms)) # Plot the quiver plot - #! TODO: figure out arguments to make arrows show up correctly + # TODO: figure out arguments to make arrows show up correctly if scale is None: - mpl.quiver(x1, x2, dx[:,:,1], dx[:,:,2], angles='xy') - elif (scale != 0): - #! TODO: optimize parameters for arrows - #! TODO: figure out arguments to make arrows show up correctly - xy = mpl.quiver(x1, x2, dx[:,:,0]*np.abs(scale), - dx[:,:,1]*np.abs(scale), angles='xy') - # set(xy, 'LineWidth', PP_arrow_linewidth, 'Color', 'b'); - - #! TODO: Tweak the shape of the plot + mpl.quiver(x1, x2, dx[:, :, 1], dx[:, :, 2], angles='xy') + elif scale != 0: + # TODO: optimize parameters for arrows + # TODO: figure out arguments to make arrows show up correctly + _ = mpl.quiver(x1, x2, dx[:, :, 0] * np.abs(scale), dx[:, :, 1] * np.abs(scale), angles='xy') + # set(xy, 'LineWidth', PP_arrow_linewidth, 'Color', 'b') + + # TODO: Tweak the shape of the plot # a=gca; set(a,'DataAspectRatio',[1,1,1]); # set(a,'XLim',X(1:2)); set(a,'YLim',Y(1:2)); - mpl.xlabel('x1'); mpl.ylabel('x2'); + mpl.xlabel('x1') + mpl.ylabel('x2') # See if we should also generate the streamlines if X0 is None or len(X0) == 0: return # Convert initial conditions to a numpy array - X0 = np.array(X0); - (nr, nc) = np.shape(X0); + X0 = np.array(X0) + (nr, nc) = np.shape(X0) # Generate some empty matrices to keep arrow information - x1 = np.empty((nr, Narrows)); x2 = np.empty((nr, Narrows)); + x1 = np.empty((nr, Narrows)) + x2 = np.empty((nr, Narrows)) dx = np.empty((nr, Narrows, 2)) # See if we were passed a simulation time @@ -203,71 +204,75 @@ def phase_plot(odefun, X=None, Y=None, scale=1, X0=None, T=None, T = 50 # Parse the time we were passed - TSPAN = T; - if (isinstance(T, (int, float))): - TSPAN = np.linspace(0, T, 100); + TSPAN = T + if isinstance(T, (int, float)): + TSPAN = np.linspace(0, T, 100) # Figure out the limits for the plot if scale is None: # Assume that the current axis are set as we want them - alim = mpl.axis(); - xmin = alim[0]; xmax = alim[1]; - ymin = alim[2]; ymax = alim[3]; + alim = mpl.axis() + xmin = alim[0] + xmax = alim[1] + ymin = alim[2] + ymax = alim[3] else: # Use the maximum extent of all trajectories - xmin = np.min(X0[:,0]); xmax = np.max(X0[:,0]); - ymin = np.min(X0[:,1]); ymax = np.max(X0[:,1]); + xmin = np.min(X0[:, 0]) + xmax = np.max(X0[:, 0]) + ymin = np.min(X0[:, 1]) + ymax = np.max(X0[:, 1]) # Generate the streamlines for each initial condition for i in range(nr): - state = odeint(odefun, X0[i], TSPAN, args=parms); + state = odeint(odefun, X0[i], TSPAN, args=parms) time = TSPAN - mpl.plot(state[:,0], state[:,1]) - #! TODO: add back in colors for stream lines + mpl.plot(state[:, 0], state[:, 1]) + # TODO: add back in colors for stream lines # PP_stream_color(np.mod(i-1, len(PP_stream_color))+1)); # set(h[i], 'LineWidth', PP_stream_linewidth); # Plot arrows if quiver parameters were 'auto' - if (autoFlag or logtimeFlag or timeptsFlag): + if autoFlag or logtimeFlag or timeptsFlag: # Compute the locations of the arrows - #! TODO: check this logic to make sure it works in python + # TODO: check this logic to make sure it works in python for j in range(Narrows): # Figure out starting index; headless arrows start at 0 - k = -1 if scale is None else 0; + k = -1 if scale is None else 0 # Figure out what time index to use for the next point - if (autoFlag): + if autoFlag: # Use a linear scaling based on ODE time vector - tind = np.floor((len(time)/Narrows) * (j-k)) + k; - elif (logtimeFlag): + tind = np.floor((len(time) / Narrows) * (j - k)) + k + elif logtimeFlag: # Use an exponential time vector # MATLAB: tind = find(time < (j-k) / lambda, 1, 'last'); - tarr = _find(time < (j-k) / timefactor); - tind = tarr[-1] if len(tarr) else 0; - elif (timeptsFlag): + tarr = _find(time < (j - k) / timefactor) + tind = tarr[-1] if len(tarr) else 0 + elif timeptsFlag: # Use specified time points # MATLAB: tind = find(time < Y[j], 1, 'last'); - tarr = _find(time < timepts[j]); - tind = tarr[-1] if len(tarr) else 0; + tarr = _find(time < timepts[j]) + tind = tarr[-1] if len(tarr) else 0 # For tailless arrows, skip the first point if tind == 0 and scale is None: - continue; + continue # Figure out the arrow at this point on the curve - x1[i,j] = state[tind, 0]; - x2[i,j] = state[tind, 1]; + x1[i, j] = state[tind, 0] + x2[i, j] = state[tind, 1] # Skip arrows outside of initial condition box - if (scale is not None or - (x1[i,j] <= xmax and x1[i,j] >= xmin and - x2[i,j] <= ymax and x2[i,j] >= ymin)): - v = odefun((x1[i,j], x2[i,j]), 0, *parms) - dx[i, j, 0] = v[0]; dx[i, j, 1] = v[1]; + if scale is not None or (xmax >= x1[i, j] >= xmin and ymax >= x2[i, j] >= ymin): + v = odefun((x1[i, j], x2[i, j]), 0, *parms) + dx[i, j, 0] = v[0] + dx[i, j, 1] = v[1] else: - dx[i, j, 0] = 0; dx[i, j, 1] = 0; + dx[i, j, 0] = 0 + dx[i, j, 1] = 0 # Set the plot shape before plotting arrows to avoid warping # a=gca; @@ -280,21 +285,21 @@ def phase_plot(odefun, X=None, Y=None, scale=1, X0=None, T=None, # Plot arrows on the streamlines if scale is None and Narrows > 0: # Use a tailless arrow - #! TODO: figure out arguments to make arrows show up correctly - mpl.quiver(x1, x2, dx[:,:,0], dx[:,:,1], angles='xy') - elif (scale != 0 and Narrows > 0): - #! TODO: figure out arguments to make arrows show up correctly - xy = mpl.quiver(x1, x2, dx[:,:,0]*abs(scale), dx[:,:,1]*abs(scale), - angles='xy') + # TODO: figure out arguments to make arrows show up correctly + mpl.quiver(x1, x2, dx[:, :, 0], dx[:, :, 1], angles='xy') + elif scale != 0 and Narrows > 0: + # TODO: figure out arguments to make arrows show up correctly + _ = mpl.quiver(x1, x2, dx[:, :, 0] * abs(scale), dx[:, :, 1] * abs(scale), angles='xy') # set(xy, 'LineWidth', PP_arrow_linewidth); # set(xy, 'AutoScale', 'off'); # set(xy, 'AutoScaleFactor', 0); - if (scale < 0): - bp = mpl.plot(x1, x2, 'b.'); # add dots at base + if scale < 0: + _ = mpl.plot(x1, x2, 'b.') # add dots at base # set(bp, 'MarkerSize', PP_arrow_markersize); - return; + return + # Utility function for generating initial conditions around a box def box_grid(xlimp, ylimp): @@ -308,7 +313,7 @@ def box_grid(xlimp, ylimp): sx10 = np.linspace(xlimp[0], xlimp[1], xlimp[2]) sy10 = np.linspace(ylimp[0], ylimp[1], ylimp[2]) - sx1 = np.hstack((0, sx10, 0*sy10+sx10[0], sx10, 0*sy10+sx10[-1])) - sx2 = np.hstack((0, 0*sx10+sy10[0], sy10, 0*sx10+sy10[-1], sy10)) + sx1 = np.hstack((0, sx10, 0 * sy10 + sx10[0], sx10, 0 * sy10 + sx10[-1])) + sx2 = np.hstack((0, 0 * sx10 + sy10[0], sy10, 0 * sx10 + sy10[-1], sy10)) - return np.transpose( np.vstack((sx1, sx2)) ) + return np.transpose(np.vstack((sx1, sx2))) diff --git a/control/pzmap.py b/control/pzmap.py index 252d10011..6a98a508f 100644 --- a/control/pzmap.py +++ b/control/pzmap.py @@ -40,13 +40,13 @@ # # $Id:pzmap.py 819 2009-05-29 21:28:07Z murray $ -from numpy import real, imag, linspace, exp, cos, sin, sqrt -from math import pi -from .lti import LTI, isdtime, isctime +from numpy import real, imag +from .lti import LTI, isdtime from .grid import sgrid, zgrid, nogrid __all__ = ['pzmap'] + # TODO: Implement more elegant cross-style axes. See: # http://matplotlib.sourceforge.net/examples/axes_grid/demo_axisline_style.html # http://matplotlib.sourceforge.net/examples/axes_grid/demo_curvelinear_grid.html @@ -77,7 +77,7 @@ def pzmap(sys, Plot=True, grid=False, title='Pole Zero Map'): poles = sys.pole() zeros = sys.zero() - if (Plot): + if Plot: import matplotlib.pyplot as plt if grid: @@ -92,9 +92,7 @@ def pzmap(sys, Plot=True, grid=False, title='Pole Zero Map'): if len(poles) > 0: ax.scatter(real(poles), imag(poles), s=50, marker='x', facecolors='k') if len(zeros) > 0: - ax.scatter(real(zeros), imag(zeros), s=50, marker='o', - facecolors='none', edgecolors='k') - + ax.scatter(real(zeros), imag(zeros), s=50, marker='o', facecolors='none', edgecolors='k') plt.title(title) diff --git a/control/rlocus.py b/control/rlocus.py index 61322624c..039fed168 100644 --- a/control/rlocus.py +++ b/control/rlocus.py @@ -47,8 +47,8 @@ # Packages used by this module import numpy as np -from scipy import array, poly1d, row_stack, zeros_like, real, imag, exp, sin, cos, linspace, sqrt -from math import pi +from scipy import array, poly1d, row_stack, zeros_like, real, imag + import scipy.signal # signal processing toolbox import pylab # plotting routines from .xferfcn import _convert_to_transfer_function @@ -86,6 +86,7 @@ def root_locus(sys, kvect=None, xlim=None, ylim=None, plotstr='-', Plot=True, calculate gain, damping and print grid: boolean (default = False) If True plot omega-damping grid. + plotstr: ? Returns ------- @@ -93,6 +94,7 @@ def root_locus(sys, kvect=None, xlim=None, ylim=None, plotstr='-', Plot=True, Computed root locations, given as a 2d array klist : ndarray or list Gains used. Same as klist keyword argument if provided. + """ # Convert numerator and denominator to polynomials if they aren't @@ -125,9 +127,8 @@ def root_locus(sys, kvect=None, xlim=None, ylim=None, plotstr='-', Plot=True, 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)) + 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)) # plot open loop poles poles = array(denp.r) @@ -151,8 +152,8 @@ def root_locus(sys, kvect=None, xlim=None, ylim=None, plotstr='-', Plot=True, def _default_gains(num, den, xlim, ylim): - """Unsupervised gains calculation for root locus plot. - + """Unsupervised gains calculation for root locus plot. + References: Ogata, K. (2002). Modern control engineering (4th ed.). Upper Saddle River, NJ : New Delhi: Prentice Hall..""" @@ -178,16 +179,15 @@ def _default_gains(num, den, xlim, ylim): mymat_xl = np.append(mymat_xl, important_points) false_gain = den.coeffs[0] / 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.") + raise ValueError("Not implemented support for 0 degrees root 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))) 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 = 0.05 * (axmax - axmin) @@ -206,10 +206,10 @@ def _default_gains(num, den, xlim, ylim): 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_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) + kvect = np.insert(kvect, index + 1, new_gains[1:4]) + mymat = np.insert(mymat, index + 1, new_points, axis=0) mymat = _RLSortRoots(mymat) distance_points = np.abs(np.diff(mymat, axis=0)) > tolerance # distance between points @@ -261,9 +261,9 @@ def _k_max(num, den, real_break_points, k_break_points): false_gain = den.coeffs[0] / num.coeffs[0] if asymp_number > 0: - asymp_center = (np.sum(den.roots) - np.sum(num.roots))/asymp_number + asymp_center = (np.sum(den.roots) - np.sum(num.roots)) / asymp_number distance_max = 4 * np.max(np.abs(important_points - asymp_center)) - asymp_angles = (2 * np.arange(0, asymp_number)-1) * np.pi / asymp_number + asymp_angles = (2 * np.arange(0, asymp_number) - 1) * np.pi / asymp_number if false_gain > 0: farthest_points = asymp_center + distance_max * np.exp(asymp_angles * 1j) # farthest points over asymptotes else: @@ -282,7 +282,7 @@ def _k_max(num, den, real_break_points, k_break_points): def _systopoly1d(sys): """Extract numerator and denominator polynomails for a system""" # Allow inputs from the signal processing toolbox - if (isinstance(sys, scipy.signal.lti)): + if isinstance(sys, scipy.signal.lti): nump = sys.num denp = sys.den @@ -291,7 +291,7 @@ def _systopoly1d(sys): sys = _convert_to_transfer_function(sys) # Make sure we have a SISO system - if (sys.inputs > 1 or sys.outputs > 1): + if sys.inputs > 1 or sys.outputs > 1: raise ControlMIMONotImplemented() # Start by extracting the numerator and denominator from system object @@ -299,13 +299,13 @@ def _systopoly1d(sys): denp = sys.den[0][0] # Check to see if num, den are already polynomials; otherwise convert - if (not isinstance(nump, poly1d)): + if not isinstance(nump, poly1d): nump = poly1d(nump) - if (not isinstance(denp, poly1d)): + if not isinstance(denp, poly1d): denp = poly1d(denp) - return (nump, denp) + return nump, denp def _RLFindRoots(nump, denp, kvect): @@ -332,37 +332,42 @@ def _RLSortRoots(mymat): locus doesn't show weird pseudo-branches as roots jump from one branch to another.""" - sorted = zeros_like(mymat) + sorted_mat = zeros_like(mymat) + prevrow = None + for n, row in enumerate(mymat): if n == 0: - sorted[n, :] = row + sorted_mat[n, :] = row + prevrow = row else: # sort the current row by finding the element with the # smallest absolute distance to each root in the # previous row available = list(range(len(prevrow))) for elem in row: - evect = elem-prevrow[available] + evect = elem - prevrow[available] ind1 = abs(evect).argmin() ind = available.pop(ind1) - sorted[n, ind] = elem - prevrow = sorted[n, :] - return sorted + sorted_mat[n, ind] = elem + prevrow = sorted_mat[n, :] + return sorted_mat -def _RLFeedbackClicks(event, sys,fig,point): + +def _RLFeedbackClicks(event, sys, fig, point): """Print root-locus gain feedback for clicks on the root-locus plot """ 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))) 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))) + (s.real, s.imag, K.real, -1 * s.real / abs(s))) fig.canvas.draw() + rlocus = root_locus diff --git a/control/robust.py b/control/robust.py index 54e970a5d..2b5c2b1b9 100644 --- a/control/robust.py +++ b/control/robust.py @@ -41,11 +41,11 @@ # External packages and modules import numpy as np -from .exception import * +from .exception import ControlSlycot from .statesp import StateSpace -from .statefbk import * -def h2syn(P,nmeas,ncon): + +def h2syn(P, nmeas, ncon): """H_2 control synthesis for plant P. Parameters @@ -72,25 +72,25 @@ def h2syn(P,nmeas,ncon): >>> K = h2syn(P,nmeas,ncon) """ - #Check for ss system object, need a utility for this? + # Check for ss system object, need a utility for this? - #TODO: Check for continous or discrete, only continuous supported right now - # if isCont(): - # dico = 'C' - # elif isDisc(): - # dico = 'D' - # else: - dico = 'C' + # TODO: Check for continous or discrete, only continuous supported right now + # if isCont(): + # dico = 'C' + # elif isDisc(): + # dico = 'D' + # else: + # dico = 'C' try: from slycot import sb10hd except ImportError: raise ControlSlycot("can't find slycot subroutine sb10hd") - n = np.size(P.A,0) - m = np.size(P.B,1) - np_ = np.size(P.C,0) - out = sb10hd(n,m,np_,ncon,nmeas,P.A,P.B,P.C,P.D) + n = np.size(P.A, 0) + m = np.size(P.B, 1) + np_ = np.size(P.C, 0) + out = sb10hd(n, m, np_, ncon, nmeas, P.A, P.B, P.C, P.D) Ak = out[0] Bk = out[1] Ck = out[2] @@ -100,7 +100,8 @@ def h2syn(P,nmeas,ncon): return K -def hinfsyn(P,nmeas,ncon): + +def hinfsyn(P, nmeas, ncon): """H_{inf} control synthesis for plant P. Parameters @@ -136,26 +137,26 @@ def hinfsyn(P,nmeas,ncon): """ - #Check for ss system object, need a utility for this? + # Check for ss system object, need a utility for this? - #TODO: Check for continous or discrete, only continuous supported right now - # if isCont(): - # dico = 'C' - # elif isDisc(): - # dico = 'D' - # else: - dico = 'C' + # TODO: Check for continous or discrete, only continuous supported right now + # if isCont(): + # dico = 'C' + # elif isDisc(): + # dico = 'D' + # else: + # dico = 'C' try: from slycot import sb10ad except ImportError: raise ControlSlycot("can't find slycot subroutine sb10ad") - n = np.size(P.A,0) - m = np.size(P.B,1) - np_ = np.size(P.C,0) + n = np.size(P.A, 0) + m = np.size(P.B, 1) + np_ = np.size(P.C, 0) gamma = 1.e100 - out = sb10ad(n,m,np_,ncon,nmeas,gamma,P.A,P.B,P.C,P.D) + out = sb10ad(n, m, np_, ncon, nmeas, gamma, P.A, P.B, P.C, P.D) gam = out[0] Ak = out[1] Bk = out[2] @@ -202,21 +203,20 @@ def _size_as_needed(w, wname, n): """ from . import append, ss if w is not None: - if not isinstance(w,StateSpace): + if not isinstance(w, StateSpace): w = ss(w) - if 1==w.inputs and 1==w.outputs: - w = append(*(w,)*n) + if 1 == w.inputs and 1 == w.outputs: + w = append(*(w,) * n) else: if w.inputs != n: - msg=("{}: weighting function has {} inputs, expected {}". - format(wname,w.inputs,n)) + msg = ("{}: weighting function has {} inputs, expected {}".format(wname, w.inputs, n)) raise ValueError(msg) else: - w = ss([],[],[],[]) + w = ss([], [], [], []) return w -def augw(g,w1=None,w2=None,w3=None): +def augw(g, w1=None, w2=None, w3=None): """Augment plant for mixed sensitivity problem. Parameters @@ -254,12 +254,12 @@ def augw(g,w1=None,w2=None,w3=None): ny = g.outputs nu = g.inputs - w1,w2,w3 = [_size_as_needed(w,wname,n) - for w,wname,n in zip((w1,w2,w3), - ('w1','w2','w3'), - (ny,nu,ny))] + w1, w2, w3 = [_size_as_needed(w, wname, n) + for w, wname, n in zip((w1, w2, w3), + ('w1', 'w2', 'w3'), + (ny, nu, ny))] - if not isinstance(g,StateSpace): + if not isinstance(g, StateSpace): g = ss(g) # w u @@ -270,11 +270,11 @@ def augw(g,w1=None,w2=None,w3=None): # v [ I | -g ] # error summer: inputs are -y and r=w - Ie = ss([],[],[],np.eye(ny)) + Ie = ss([], [], [], np.eye(ny)) # control: needed to "distribute" control input - Iu = ss([],[],[],np.eye(nu)) + Iu = ss([], [], [], np.eye(nu)) - sysall = append(w1,w2,w3,Ie,g,Iu) + sysall = append(w1, w2, w3, Ie, g, Iu) niw1 = w1.inputs niw2 = w2.inputs @@ -284,43 +284,44 @@ def augw(g,w1=None,w2=None,w3=None): now2 = w2.outputs now3 = w3.outputs - q = np.zeros((niw1+niw2+niw3+ny+nu,2)) - q[:,0] = np.arange(1,q.shape[0]+1) + q = np.zeros((niw1 + niw2 + niw3 + ny + nu, 2)) + q[:, 0] = np.arange(1, q.shape[0] + 1) # Ie -> w1 - q[:niw1,1] = np.arange(1+now1+now2+now3, - 1+now1+now2+now3+niw1) + q[:niw1, 1] = np.arange(1 + now1 + now2 + now3, + 1 + now1 + now2 + now3 + niw1) # Iu -> w2 - q[niw1:niw1+niw2,1] = np.arange(1+now1+now2+now3+2*ny, - 1+now1+now2+now3+2*ny+niw2) + q[niw1:niw1 + niw2, 1] = np.arange(1 + now1 + now2 + now3 + 2 * ny, + 1 + now1 + now2 + now3 + 2 * ny + niw2) # y -> w3 - q[niw1+niw2:niw1+niw2+niw3,1] = np.arange(1+now1+now2+now3+ny, - 1+now1+now2+now3+ny+niw3) + q[niw1 + niw2:niw1 + niw2 + niw3, 1] = np.arange(1 + now1 + now2 + now3 + ny, + 1 + now1 + now2 + now3 + ny + niw3) # -y -> Iy; note the leading - - q[niw1+niw2+niw3:niw1+niw2+niw3+ny,1] = -np.arange(1+now1+now2+now3+ny, - 1+now1+now2+now3+2*ny) + q[niw1 + niw2 + niw3:niw1 + niw2 + niw3 + ny, 1] = -np.arange(1 + now1 + now2 + now3 + ny, + 1 + now1 + now2 + now3 + 2 * ny) # Iu -> G - q[niw1+niw2+niw3+ny:niw1+niw2+niw3+ny+nu,1] = np.arange(1+now1+now2+now3+2*ny, - 1+now1+now2+now3+2*ny+nu) + q[niw1 + niw2 + niw3 + ny:niw1 + niw2 + niw3 + ny + nu, 1] = np.arange(1 + now1 + now2 + now3 + 2 * ny, + 1 + now1 + now2 + now3 + 2 * ny + nu) # input indices: to Ie and Iu - ii = np.hstack((np.arange(1+now1+now2+now3, - 1+now1+now2+now3+ny), - np.arange(1+now1+now2+now3+ny+nu, - 1+now1+now2+now3+ny+nu+nu))) + ii = np.hstack((np.arange(1 + now1 + now2 + now3, + 1 + now1 + now2 + now3 + ny), + np.arange(1 + now1 + now2 + now3 + ny + nu, + 1 + now1 + now2 + now3 + ny + nu + nu))) # output indices - oi = np.arange(1,1+now1+now2+now3+ny) + oi = np.arange(1, 1 + now1 + now2 + now3 + ny) - p = connect(sysall,q,ii,oi) + p = connect(sysall, q, ii, oi) return p -def mixsyn(g,w1=None,w2=None,w3=None): + +def mixsyn(g, w1=None, w2=None, w3=None): """Mixed-sensitivity H-infinity synthesis. mixsyn(g,w1,w2,w3) -> k,cl,info @@ -336,9 +337,9 @@ def mixsyn(g,w1=None,w2=None,w3=None): Returns ------- k: synthesized controller; StateSpace object - cl: closed system mapping evaluation inputs to evaluation outputs; if + cl: closed system mapping evaluation inputs to evaluation outputs; if p is the augmented plant, with - [z] = [p11 p12] [w], + [z] = [p11 p12] [w], [y] [p21 g] [u] then cl is the system from w->z with u=-k*y. StateSpace object. @@ -356,8 +357,8 @@ def mixsyn(g,w1=None,w2=None,w3=None): """ nmeas = g.outputs ncon = g.inputs - p = augw(g,w1,w2,w3) + p = augw(g, w1, w2, w3) - k,cl,gamma,rcond=hinfsyn(p,nmeas,ncon) - info = gamma,rcond - return k,cl,info + k, cl, gamma, rcond = hinfsyn(p, nmeas, ncon) + info = gamma, rcond + return k, cl, info diff --git a/control/statefbk.py b/control/statefbk.py index 0fb377a47..0d3e2549c 100644 --- a/control/statefbk.py +++ b/control/statefbk.py @@ -41,7 +41,6 @@ # External packages and modules import numpy as np -import scipy as sp from . import statesp from .exception import ControlSlycot, ControlArgument, ControlDimension @@ -98,10 +97,10 @@ def place(A, B, p): # Convert the system inputs to NumPy arrays A_mat = np.array(A) B_mat = np.array(B) - if (A_mat.shape[0] != A_mat.shape[1]): + if A_mat.shape[0] != A_mat.shape[1]: raise ControlDimension("A must be a square matrix") - if (A_mat.shape[0] != B_mat.shape[0]): + if A_mat.shape[0] != B_mat.shape[0]: err_str = "The number of rows of A must equal the number of rows in B" raise ControlDimension(err_str) @@ -176,8 +175,7 @@ def place_varga(A, B, p, dtime=False, alpha=None): # Convert the system inputs to NumPy arrays A_mat = np.array(A) B_mat = np.array(B) - if (A_mat.shape[0] != A_mat.shape[1] or - A_mat.shape[0] != B_mat.shape[0]): + if (A_mat.shape[0] != A_mat.shape[1] or A_mat.shape[0] != B_mat.shape[0]): raise ControlDimension("matrix dimensions are incorrect") # Compute the system eigenvalues and convert poles to numpy array @@ -205,19 +203,19 @@ def place_varga(A, B, p, dtime=False, alpha=None): # the same as what slicot thinks are the eigs. So we need some # numerical breathing room. The following is pretty heuristic, # but does the trick - alpha = -2*abs(min(system_eigs.real)) + alpha = -2 * abs(min(system_eigs.real)) elif dtime and alpha < 0.0: raise ValueError("Need alpha > 0 when DICO='D'") - # Call SLICOT routine to place the eigenvalues - A_z,w,nfp,nap,nup,F,Z = \ + A_z, w, nfp, nap, nup, F, Z = \ sb01bd(B_mat.shape[0], B_mat.shape[1], len(placed_eigs), alpha, A_mat, B_mat, placed_eigs, DICO) # Return the gain matrix, with MATLAB gain convention return -F + # Contributed by Roberto Bucher def acker(A, B, poles): """Pole placement using Ackermann method @@ -239,12 +237,12 @@ def acker(A, B, poles): """ # Convert the inputs to matrices - a = np.mat(A) - b = np.mat(B) + a_matrix = np.mat(A) + b_matrix = np.mat(B) # Make sure the system is controllable - ct = ctrb(A, B) - if np.linalg.matrix_rank(ct) != a.shape[0]: + ct = ctrb(a_matrix, b_matrix) + if np.linalg.matrix_rank(ct) != a_matrix.shape[0]: raise ValueError("System not reachable; pole placement invalid") # Compute the desired characteristic polynomial @@ -252,15 +250,16 @@ def acker(A, B, poles): # Place the poles using Ackermann's method n = np.size(p) - pmat = p[n-1]*a**0 - for i in np.arange(1,n): - pmat = pmat + p[n-i-1]*a**i + pmat = p[n - 1] * a_matrix**0 + for i in np.arange(1, n): + pmat = pmat + p[n - i - 1] * a_matrix**i K = np.linalg.solve(ct, pmat) K = K[-1][:] # Extract the last row return K -def lqr(*args, **keywords): + +def lqr(*args): """lqr(A, B, Q, R[, N]) Linear quadratic regulator design @@ -319,56 +318,56 @@ def lqr(*args, **keywords): # # Get the system description - if (len(args) < 3): + if len(args) < 3: raise ControlArgument("not enough input arguments") try: # If this works, we were (probably) passed a system as the # first argument; extract A and B - A = np.array(args[0].A, ndmin=2, dtype=float); - B = np.array(args[0].B, ndmin=2, dtype=float); - index = 1; + A = np.array(args[0].A, ndmin=2, dtype=float) + B = np.array(args[0].B, ndmin=2, dtype=float) + index = 1 except AttributeError: # Arguments should be A and B matrices - A = np.array(args[0], ndmin=2, dtype=float); - B = np.array(args[1], ndmin=2, dtype=float); - index = 2; + A = np.array(args[0], ndmin=2, dtype=float) + B = np.array(args[1], ndmin=2, dtype=float) + index = 2 # Get the weighting matrices (converting to matrices, if needed) - Q = np.array(args[index], ndmin=2, dtype=float); - R = np.array(args[index+1], ndmin=2, dtype=float); - if (len(args) > index + 2): - N = np.array(args[index+2], ndmin=2, dtype=float); + Q = np.array(args[index], ndmin=2, dtype=float) + R = np.array(args[index + 1], ndmin=2, dtype=float) + if len(args) > index + 2: + N = np.array(args[index + 2], ndmin=2, dtype=float) else: - N = np.zeros((Q.shape[0], R.shape[1])); + N = np.zeros((Q.shape[0], R.shape[1])) # Check dimensions for consistency - nstates = B.shape[0]; - ninputs = B.shape[1]; - if (A.shape[0] != nstates or A.shape[1] != nstates): + nstates = B.shape[0] + ninputs = B.shape[1] + if A.shape[0] != nstates or A.shape[1] != nstates: raise ControlDimension("inconsistent system dimensions") - elif (Q.shape[0] != nstates or Q.shape[1] != nstates or - R.shape[0] != ninputs or R.shape[1] != ninputs or - N.shape[0] != nstates or N.shape[1] != ninputs): + elif Q.shape[0] != nstates or Q.shape[1] != nstates or\ + R.shape[0] != ninputs or R.shape[1] != ninputs or \ + N.shape[0] != nstates or N.shape[1] != ninputs: raise ControlDimension("incorrect weighting matrix dimensions") # Compute the G matrix required by SB02MD - A_b,B_b,Q_b,R_b,L_b,ipiv,oufact,G = \ - sb02mt(nstates, ninputs, B, R, A, Q, N, jobl='N'); + A_b, B_b, Q_b, R_b, L_b, ipiv, oufact, G = sb02mt(nstates, ninputs, B, R, A, Q, N, jobl='N') # Call the SLICOT function - X,rcond,w,S,U,A_inv = sb02md(nstates, A_b, G, Q_b, 'C') + X, rcond, w, S, U, A_inv = sb02md(nstates, A_b, G, Q_b, 'C') # Now compute the return value # We assume that R is positive definite and, hence, invertible - K = np.linalg.solve(R, np.dot(B.T, X) + N.T); - S = X; - E = w[0:nstates]; + K = np.linalg.solve(R, np.dot(B.T, X) + N.T) + S = X + E = w[0:nstates] return K, S, E -def ctrb(A,B): + +def ctrb(A, B): """Controllabilty matrix Parameters @@ -394,9 +393,10 @@ def ctrb(A,B): # Construct the controllability matrix ctrb = bmat for i in range(1, n): - ctrb = np.hstack((ctrb, amat**i*bmat)) + ctrb = np.hstack((ctrb, amat**i * bmat)) return ctrb + def obsv(A, C): """Observability matrix @@ -422,12 +422,13 @@ def obsv(A, C): n = np.shape(amat)[0] # Construct the controllability matrix - obsv = cmat + return_value = cmat for i in range(1, n): - obsv = np.vstack((obsv, cmat*amat**i)) - return obsv + return_value = np.vstack((return_value, cmat * amat**i)) + return return_value + -def gram(sys,type): +def gram(sys, type): """Gramian (controllability or observability) Parameters @@ -461,67 +462,68 @@ def gram(sys,type): >>> Wo = gram(sys,'o') >>> Rc = gram(sys,'cf'), where Wc=Rc'*Rc >>> Ro = gram(sys,'of'), where Wo=Ro'*Ro - """ - #Check for ss system object - if not isinstance(sys,statesp.StateSpace): + # Check for ss system object + if not isinstance(sys, statesp.StateSpace): raise ValueError("System must be StateSpace!") if type not in ['c', 'o', 'cf', 'of']: raise ValueError("That type is not supported!") - #TODO: Check for continous or discrete, only continuous supported right now - # if isCont(): - # dico = 'C' - # elif isDisc(): - # dico = 'D' - # else: + # TODO: Check for continous or discrete, only continuous supported right now + # if isCont(): + # dico = 'C' + # elif isDisc(): + # dico = 'D' + dico = 'C' - #TODO: Check system is stable, perhaps a utility in ctrlutil.py - # or a method of the StateSpace class? + # TODO: Check system is stable, perhaps a utility in ctrlutil.py or a method of the StateSpace class? if np.any(np.linalg.eigvals(sys.A).real >= 0.0): raise ValueError("Oops, the system is unstable!") - if type=='c' or type=='o': - #Compute Gramian by the Slycot routine sb03md - #make sure Slycot is installed + if type == 'c' or type == 'o': + # Compute Gramian by the Slycot routine sb03md + # make sure Slycot is installed try: from slycot import sb03md except ImportError: raise ControlSlycot("can't find slycot module 'sb03md'") - if type=='c': + + if type == 'c': tra = 'T' - C = -np.dot(sys.B,sys.B.transpose()) - elif type=='o': + C = -np.dot(sys.B, sys.B.transpose()) + elif type == 'o': tra = 'N' - C = -np.dot(sys.C.transpose(),sys.C) + C = -np.dot(sys.C.transpose(), sys.C) + n = sys.states - U = np.zeros((n,n)) + U = np.zeros((n, n)) A = np.array(sys.A) # convert to NumPy array for slycot - X,scale,sep,ferr,w = sb03md(n, C, A, U, dico, job='X', fact='N', trana=tra) + X, scale, sep, ferr, w = sb03md(n, C, A, U, dico, job='X', fact='N', trana=tra) gram = X return gram - elif type=='cf' or type=='of': - #Compute cholesky factored gramian from slycot routine sb03od + elif type == 'cf' or type == 'of': + # Compute cholesky factored gramian from slycot routine sb03od try: from slycot import sb03od except ImportError: raise ControlSlycot("can't find slycot module 'sb03od'") - tra='N' + + tra = 'N' n = sys.states - Q = np.zeros((n,n)) + Q = np.zeros((n, n)) A = np.array(sys.A) # convert to NumPy array for slycot - if type=='cf': + if type == 'cf': m = sys.B.shape[1] B = np.zeros_like(A) - B[0:m,0:n] = sys.B.transpose() - X,scale,w = sb03od(n, m, A.transpose(), Q, B, dico, fact='N', trans=tra) - elif type=='of': + B[0:m, 0:n] = sys.B.transpose() + X, scale, w = sb03od(n, m, A.transpose(), Q, B, dico, fact='N', trans=tra) + elif type == 'of': m = sys.C.shape[0] C = np.zeros_like(A) - C[0:n,0:m] = sys.C.transpose() - X,scale,w = sb03od(n, m, A, Q, C.transpose(), dico, fact='N', trans=tra) + C[0:n, 0:m] = sys.C.transpose() + X, scale, w = sb03od(n, m, A, Q, C.transpose(), dico, fact='N', trans=tra) gram = X return gram diff --git a/control/statesp.py b/control/statesp.py index 172b392c0..516990799 100644 --- a/control/statesp.py +++ b/control/statesp.py @@ -9,8 +9,8 @@ """ # Python 3 compatibility (needs to go here) +from __future__ import division # for _convertToStateSpace from __future__ import print_function -from __future__ import division # for _convertToStateSpace """Copyright (c) 2010 by California Institute of Technology All rights reserved. @@ -53,8 +53,8 @@ import math import numpy as np -from numpy import all, angle, any, array, asarray, concatenate, cos, delete, \ - dot, empty, exp, eye, isinf, matrix, ones, pad, shape, sin, zeros, squeeze +from numpy import all, any, array, asarray, concatenate, cos, delete, \ + dot, empty, exp, eye, isinf, ones, pad, sin, zeros, squeeze from numpy.random import rand, randn from numpy.linalg import solve, eigvals, matrix_rank from numpy.linalg.linalg import LinAlgError @@ -62,7 +62,6 @@ from scipy.signal import lti, cont2discrete from warnings import warn from .lti import LTI, timebase, timebaseEqual, isdtime -from .xferfcn import _convert_to_transfer_function from copy import deepcopy __all__ = ['StateSpace', 'ss', 'rss', 'drss', 'tf2ss', 'ssdata'] @@ -166,9 +165,9 @@ def __init__(self, *args): if 0 == self.states: # static gain # matrix's default "empty" shape is 1x0 - A.shape = (0,0) - B.shape = (0,self.inputs) - C.shape = (self.outputs,0) + A.shape = (0, 0) + B.shape = (0, self.inputs) + C.shape = (self.outputs, 0) # Check that the matrix sizes are consistent. if self.states != A.shape[0]: @@ -199,14 +198,14 @@ def _remove_useless_states(self): # Search for useless states. for i in range(self.states): - if (all(self.A[i, :] == zeros((1, self.states))) and - all(self.B[i, :] == zeros((1, self.inputs)))): + if all(self.A[i, :] == zeros((1, self.states))) and \ + all(self.B[i, :] == zeros((1, self.inputs))): useless.append(i) # To avoid duplicate indices in useless, jump to the next # iteration. continue - if (all(self.A[:, i] == zeros((self.states, 1))) and - all(self.C[:, i] == zeros((self.outputs, 1)))): + if all(self.A[:, i] == zeros((self.states, 1))) and \ + all(self.C[:, i] == zeros((self.outputs, 1))): useless.append(i) # Remove the useless states. @@ -222,16 +221,16 @@ def _remove_useless_states(self): def __str__(self): """String representation of the state space.""" - str = "A = " + self.A.__str__() + "\n\n" - str += "B = " + self.B.__str__() + "\n\n" - str += "C = " + self.C.__str__() + "\n\n" - str += "D = " + self.D.__str__() + "\n" + _str = "A = " + self.A.__str__() + "\n\n" + _str += "B = " + self.B.__str__() + "\n\n" + _str += "C = " + self.C.__str__() + "\n\n" + _str += "D = " + self.D.__str__() + "\n" # TODO: replace with standard calls to lti functions - if (type(self.dt) == bool and self.dt is True): - str += "\ndt unspecified\n" - elif (not (self.dt is None) and type(self.dt) != bool and self.dt > 0): - str += "\ndt = " + self.dt.__str__() + "\n" - return str + if type(self.dt) == bool and self.dt is True: + _str += "\ndt unspecified\n" + elif not (self.dt is None) and type(self.dt) != bool and self.dt > 0: + _str += "\ndt = " + self.dt.__str__() + "\n" + return _str # represent as string, makes display work for IPython __repr__ = __str__ @@ -256,26 +255,20 @@ def __add__(self, other): other = _convertToStateSpace(other) # Check to make sure the dimensions are OK - if ((self.inputs != other.inputs) or - (self.outputs != other.outputs)): + if self.inputs != other.inputs or self.outputs != other.outputs: raise ValueError("Systems have different shapes.") # Figure out the sampling time to use if self.dt is None and other.dt is not None: dt = other.dt # use dt from second argument - elif (other.dt is None and self.dt is not None) or \ - (timebaseEqual(self, other)): + elif (other.dt is None and self.dt is not None) or timebaseEqual(self, other): dt = self.dt # use dt from first argument else: raise ValueError("Systems have different sampling times") # Concatenate the various arrays - A = concatenate(( - concatenate((self.A, zeros((self.A.shape[0], - other.A.shape[-1]))),axis=1), - concatenate((zeros((other.A.shape[0], self.A.shape[-1])), - other.A),axis=1) - ),axis=0) + A = concatenate((concatenate((self.A, zeros((self.A.shape[0], other.A.shape[-1]))), axis=1), + concatenate((zeros((other.A.shape[0], self.A.shape[-1])), other.A), axis=1)), axis=0) B = concatenate((self.B, other.B), axis=0) C = concatenate((self.C, other.C), axis=1) D = self.D + other.D @@ -315,25 +308,21 @@ def __mul__(self, other): # Check to make sure the dimensions are OK if self.inputs != other.outputs: - raise ValueError("C = A * B: A has %i column(s) (input(s)), \ -but B has %i row(s)\n(output(s))." % (self.inputs, other.outputs)) + raise ValueError("C = A * B: A has %i column(s) (input(s)), but B has %i row(s)\n(output(s))." % (self.inputs, other.outputs)) # Figure out the sampling time to use - if (self.dt == None and other.dt != None): + if self.dt is None and other.dt is not None: dt = other.dt # use dt from second argument - elif (other.dt == None and self.dt != None) or \ - (timebaseEqual(self, other)): + elif (other.dt is None and self.dt is not None) or timebaseEqual(self, other): dt = self.dt # use dt from first argument else: raise ValueError("Systems have different sampling times") # Concatenate the various arrays - A = concatenate( - (concatenate((other.A, zeros((other.A.shape[0], self.A.shape[1]))), - axis=1), - concatenate((self.B * other.C, self.A), axis=1)), axis=0) + A = concatenate((concatenate((other.A, zeros((other.A.shape[0], self.A.shape[1]))), axis=1), + concatenate((self.B * other.C, self.A), axis=1)), axis=0) B = concatenate((other.B, self.B * other.D), axis=0) - C = concatenate((self.D * other.C, self.C),axis=1) + C = concatenate((self.D * other.C, self.C), axis=1) D = self.D * other.D return StateSpace(A, B, C, D, dt) @@ -409,8 +398,7 @@ def horner(self, s): Returns a matrix of values evaluated at complex variable s. """ - resp = self.C * solve(s * eye(self.states) - self.A, - self.B) + self.D + resp = self.C * solve(s * eye(self.states) - self.A, self.B) + self.D return array(resp) # Method for generating the frequency response of the system @@ -452,8 +440,7 @@ def freqresp(self, omega): omega = np.asarray(omega) numFreqs = len(omega) - Gfrf = np.empty((self.outputs, self.inputs, numFreqs), - dtype=np.complex128) + Gfrf = np.empty((self.outputs, self.inputs, numFreqs), dtype=np.complex128) # Sort frequency and calculate complex frequencies on either imaginary # axis (continuous time) or unit circle (discrete time). @@ -476,8 +463,7 @@ def freqresp(self, omega): p = self.outputs # The first call both evaluates C(sI-A)^-1 B and also returns # Hessenberg transformed matrices at, bt, ct. - result = tb05ad(n, m, p, cmplx_freqs[0], self.A, - self.B, self.C, job='NG') + result = tb05ad(n, m, p, cmplx_freqs[0], self.A, self.B, self.C, job='NG') # When job='NG', result = (at, bt, ct, g_i, hinvb, info) at = result[0] bt = result[1] @@ -491,13 +477,12 @@ def freqresp(self, omega): # Start at the second frequency, already have the first. for kk, cmplx_freqs_kk in enumerate(cmplx_freqs[1:numFreqs]): - result = tb05ad(n, m, p, cmplx_freqs_kk, at, - bt, ct, job='NH') + result = tb05ad(n, m, p, cmplx_freqs_kk, at, bt, ct, job='NH') # When job='NH', result = (g_i, hinvb, info) # kk+1 because enumerate starts at kk = 0. # but zero-th spot is already filled. - Gfrf[:, :, kk+1] = result[0] + self.D + Gfrf[:, :, kk + 1] = result[0] + self.D except ImportError: # Slycot unavailable. Fall back to horner. for kk, cmplx_freqs_kk in enumerate(cmplx_freqs): @@ -562,8 +547,7 @@ def feedback(self, other=1, sign=-1): # Check to make sure the dimensions are OK if (self.inputs != other.outputs) or (self.outputs != other.inputs): - raise ValueError("State space systems don't have compatible inputs/outputs for " - "feedback.") + raise ValueError("State space systems don't have compatible inputs/outputs for feedback.") # Figure out the sampling time to use if self.dt is None and other.dt is not None: @@ -614,20 +598,20 @@ def minreal(self, tol=0.0): try: from slycot import tb01pd B = empty((self.states, max(self.inputs, self.outputs))) - B[:,:self.inputs] = self.B + B[:, :self.inputs] = self.B C = empty((max(self.outputs, self.inputs), self.states)) - C[:self.outputs,:] = self.C + C[:self.outputs, :] = self.C A, B, C, nr = tb01pd(self.states, self.inputs, self.outputs, self.A, B, C, tol=tol) - return StateSpace(A[:nr,:nr], B[:nr,:self.inputs], - C[:self.outputs,:nr], self.D) + return StateSpace(A[:nr, :nr], B[:nr, :self.inputs], + C[:self.outputs, :nr], self.D) except ImportError: raise TypeError("minreal requires slycot tb01pd") else: return StateSpace(self) - # TODO: add discrete time check + def returnScipySignalLTI(self): """Return a list of a list of scipy.signal.lti objects. @@ -750,7 +734,7 @@ def dcgain(self): """ try: if self.isctime(): - gain = np.asarray(self.D-self.C.dot(np.linalg.solve(self.A, self.B))) + gain = np.asarray(self.D - self.C.dot(np.linalg.solve(self.A, self.B))) else: gain = self.horner(1) except LinAlgError: @@ -779,8 +763,7 @@ def _convertToStateSpace(sys, **kw): import itertools if isinstance(sys, StateSpace): if len(kw): - raise TypeError("If sys is a StateSpace, _convertToStateSpace \ -cannot take keywords.") + raise TypeError("If sys is a StateSpace, _convertToStateSpace cannot take keywords.") # Already a state space system; just return it return sys @@ -788,8 +771,7 @@ def _convertToStateSpace(sys, **kw): try: from slycot import td04ad if len(kw): - raise TypeError("If sys is a TransferFunction, " - "_convertToStateSpace cannot take keywords.") + raise TypeError("If sys is a TransferFunction, _convertToStateSpace cannot take keywords.") # Change the numerator and denominator arrays so that the transfer # function matrix has a common denominator. @@ -797,8 +779,7 @@ def _convertToStateSpace(sys, **kw): num, den, denorder = sys.minreal()._common_den() # transfer function to state space conversion now should work! - ssout = td04ad('C', sys.inputs, sys.outputs, - denorder, den, num, tol=0) + ssout = td04ad('C', sys.inputs, sys.outputs, denorder, den, num, tol=0) states = ssout[0] return StateSpace(ssout[1][:states, :states], ssout[2][:states, :sys.inputs], @@ -839,19 +820,20 @@ def _convertToStateSpace(sys, **kw): # The following Doesn't work due to inconsistencies in ltisys: # return StateSpace([[]], [[]], [[]], eye(outputs, inputs)) return StateSpace(0., zeros((1, inputs)), zeros((outputs, 1)), - sys * ones((outputs, inputs))) + sys * ones((outputs, inputs))) # If this is a matrix, try to create a constant feedthrough try: D = _matrix(sys) return StateSpace([], [], [], D) except Exception as e: - print("Failure to assume argument is matrix-like in" \ - " _convertToStateSpace, result %s" % e) + print("Failure to assume argument is matrix-like in _convertToStateSpace, result %s" % e) raise TypeError("Can't convert given type to StateSpace system.") # TODO: add discrete time option + + def _rss_generate(states, inputs, outputs, type): """Generate a random state space. @@ -875,14 +857,13 @@ def _rss_generate(states, inputs, outputs, type): # Check for valid input arguments. if states < 1 or states % 1: - raise ValueError("states must be a positive integer. states = %g." % - states) + raise ValueError("states must be a positive integer. states = %g." % states) + if inputs < 1 or inputs % 1: - raise ValueError("inputs must be a positive integer. inputs = %g." % - inputs) + raise ValueError("inputs must be a positive integer. inputs = %g." % inputs) + if outputs < 1 or outputs % 1: - raise ValueError("outputs must be a positive integer. outputs = %g." % - outputs) + raise ValueError("outputs must be a positive integer. outputs = %g." % outputs) # Make some poles for A. Preallocate a complex array. poles = zeros(states) + zeros(states) * 0.j @@ -892,13 +873,13 @@ def _rss_generate(states, inputs, outputs, type): if rand() < pRepeat and i != 0 and i != states - 1: # Small chance of copying poles, if we're not at the first or last # element. - if poles[i-1].imag == 0: + if poles[i - 1].imag == 0: # Copy previous real pole. - poles[i] = poles[i-1] + poles[i] = poles[i - 1] i += 1 else: # Copy previous complex conjugate pair of poles. - poles[i:i+2] = poles[i-2:i] + poles[i:i + 2] = poles[i - 2:i] i += 2 elif rand() < pReal or i == states - 1: # No-oscillation pole. @@ -915,7 +896,7 @@ def _rss_generate(states, inputs, outputs, type): mag = rand() phase = 2. * math.pi * rand() poles[i] = complex(mag * cos(phase), mag * sin(phase)) - poles[i+1] = complex(poles[i].real, -poles[i].imag) + poles[i + 1] = complex(poles[i].real, -poles[i].imag) i += 2 # Now put the poles in A as real blocks on the diagonal. @@ -926,9 +907,9 @@ def _rss_generate(states, inputs, outputs, type): A[i, i] = poles[i].real i += 1 else: - A[i, i] = A[i+1, i+1] = poles[i].real - A[i, i+1] = poles[i].imag - A[i+1, i] = -poles[i].imag + A[i, i] = A[i + 1, i + 1] = poles[i].real + A[i, i + 1] = poles[i].imag + A[i + 1, i] = -poles[i].imag i += 2 # Finally, apply a transformation so that A is not block-diagonal. while True: @@ -970,7 +951,7 @@ def _rss_generate(states, inputs, outputs, type): # Convert a MIMO system to a SISO system # TODO: add discrete time check def _mimo2siso(sys, input, output, warn_conversion=False): - #pylint: disable=W0622 + # pylint: disable=W0622 """ Convert a MIMO system to a SISO system. (Convert a system with multiple inputs and/or outputs, to a system with a single input and output.) @@ -994,12 +975,12 @@ def _mimo2siso(sys, input, output, warn_conversion=False): Warn that a conversion will take place. Returns + ------- sys: StateSpace The converted (SISO) system. """ if not (isinstance(input, int) and isinstance(output, int)): - raise TypeError("Parameters ``input`` and ``output`` must both " - "be integer numbers.") + raise TypeError("Parameters ``input`` and ``output`` must both be integer numbers.") if not (0 <= input < sys.inputs): raise ValueError("Selected input does not exist. " "Selected input: {sel}, " @@ -1010,7 +991,7 @@ def _mimo2siso(sys, input, output, warn_conversion=False): "Selected output: {sel}, " "number of system outputs: {ext}." .format(sel=output, ext=sys.outputs)) - #Convert sys to SISO if necessary + # Convert sys to SISO if necessary if sys.inputs > 1 or sys.outputs > 1: if warn_conversion: warn("Converting MIMO system to SISO system. " @@ -1027,7 +1008,6 @@ def _mimo2siso(sys, input, output, warn_conversion=False): def _mimo2simo(sys, input, warn_conversion=False): - # pylint: disable=W0622 """ Convert a MIMO system to a SIMO system. (Convert a system with multiple inputs and/or outputs, to a system with a single input but possibly @@ -1064,8 +1044,7 @@ def _mimo2simo(sys, input, warn_conversion=False): # Convert sys to SISO if necessary if sys.inputs > 1: if warn_conversion: - warn("Converting MIMO system to SIMO system. " - "Only input {i} is used." .format(i=input)) + warn("Converting MIMO system to SIMO system. Only input {i} is used." .format(i=input)) # $X = A*X + B*U # Y = C*X + D*U new_B = sys.B[:, input] @@ -1161,8 +1140,7 @@ def ss(*args): elif isinstance(sys, TransferFunction): return tf2ss(sys) else: - raise TypeError("ss(sys): sys must be a StateSpace or \ -TransferFunction object. It is %s." % type(sys)) + raise TypeError("ss(sys): sys must be a StateSpace or TransferFunction object. It is %s." % type(sys)) else: raise ValueError("Needs 1 or 4 arguments; received %i." % len(args)) @@ -1232,8 +1210,7 @@ def tf2ss(*args): elif len(args) == 1: sys = args[0] if not isinstance(sys, TransferFunction): - raise TypeError("tf2ss(sys): sys must be a TransferFunction \ -object.") + raise TypeError("tf2ss(sys): sys must be a TransferFunction object.") return _convertToStateSpace(sys) else: raise ValueError("Needs 1 or 2 arguments; received %i." % len(args)) @@ -1329,5 +1306,5 @@ def ssdata(sys): (A, B, C, D): list of matrices State space data for the system """ - ss = _convertToStateSpace(sys) - return ss.A, ss.B, ss.C, ss.D + _ss = _convertToStateSpace(sys) + return _ss.A, _ss.B, _ss.C, _ss.D diff --git a/control/tests/matlab_test.py b/control/tests/matlab_test.py index d187f6125..d1c3bf7d8 100644 --- a/control/tests/matlab_test.py +++ b/control/tests/matlab_test.py @@ -551,14 +551,12 @@ def testConnect2(self): np.testing.assert_array_almost_equal( sysc.D, np.mat([[1], [0], [0]])) - - def testFRD(self): h = tf([1], [1, 2, 2]) omega = np.logspace(-1, 2, 10) frd1 = frd(h, omega) assert isinstance(frd1, FRD) - frd2 = frd(frd1.fresp[0,0,:], omega) + frd2 = frd(frd1.fresp[0, 0, :], omega) assert isinstance(frd2, FRD) @unittest.skipIf(not slycot_check(), "slycot not installed") diff --git a/control/timeresp.py b/control/timeresp.py index a0ecf9d93..923166145 100644 --- a/control/timeresp.py +++ b/control/timeresp.py @@ -61,6 +61,7 @@ __all__ = ['forced_response', 'step_response', 'initial_response', 'impulse_response'] + # Helper function for checking array-like parameters def _check_convert_array(in_obj, legal_shapes, err_msg_start, squeeze=False, transpose=False): @@ -112,11 +113,11 @@ def _check_convert_array(in_obj, legal_shapes, err_msg_start, squeeze=False, """ # convert nearly everything to an array. out_array = np.asarray(in_obj) - if (transpose): + if transpose: out_array = np.transpose(out_array) # Test element data type, elements must be numbers - legal_kinds = set(("i", "f", "c")) # integer, float, complex + legal_kinds = {"i", "f", "c"} # integer, float, complex if out_array.dtype.kind not in legal_kinds: err_msg = "Wrong element data type: '{d}'. Array elements " \ "must be numbers.".format(d=str(out_array.dtype)) @@ -223,9 +224,8 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False): raise TypeError('Parameter ``sys``: must be a ``LTI`` object. ' '(For example ``StateSpace`` or ``TransferFunction``)') sys = _convertToStateSpace(sys) - A, B, C, D = np.asarray(sys.A), np.asarray(sys.B), np.asarray(sys.C), \ - np.asarray(sys.D) -# d_type = A.dtype + A, B, C, D = np.asarray(sys.A), np.asarray(sys.B), np.asarray(sys.C), np.asarray(sys.D) + n_states = A.shape[0] n_inputs = B.shape[1] n_outputs = C.shape[0] @@ -234,23 +234,20 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False): if isdtime(sys, strict=True): if T is None: if U is None: - raise ValueError('Parameters ``T`` and ``U`` can\'t both be' - 'zero for discrete-time simulation') + raise ValueError('Parameters ``T`` and ``U`` can\'t both be zero for discrete-time simulation') # Set T to integers with same length as U T = range(len(U)) else: # Make sure the input vector and time vector have same length # TODO: allow interpolation of the input vector if len(U) != len(T): - ValueError('Pamameter ``T`` must have same length as' - 'input vector ``U``') + ValueError('Pamameter ``T`` must have same length as input vector ``U``') # Test if T has shape (n,) or (1, n); # T must be array-like and values must be increasing. # The length of T determines the length of the input vector. if T is None: - raise ValueError('Parameter ``T``: must be array-like, and contain ' - '(strictly monotonic) increasing numbers.') + raise ValueError('Parameter ``T``: must be array-like, and contain (strictly monotonic) increasing numbers.') T = _check_convert_array(T, [('any',), (1, 'any')], 'Parameter ``T``: ', squeeze=True, transpose=transpose) @@ -260,8 +257,7 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False): n_steps = len(T) # number of simulation steps # create X0 if not given, test if X0 has correct shape - X0 = _check_convert_array(X0, [(n_states,), (n_states, 1)], - 'Parameter ``X0``: ', squeeze=True) + X0 = _check_convert_array(X0, [(n_states,), (n_states, 1)], 'Parameter ``X0``: ', squeeze=True) xout = np.zeros((n_states, n_steps)) xout[:, 0] = X0 @@ -277,20 +273,17 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False): # Solve using matrix exponential expAdt = sp.linalg.expm(A * dt) for i in range(1, n_steps): - xout[:, i] = dot(expAdt, xout[:, i-1]) + xout[:, i] = dot(expAdt, xout[:, i - 1]) yout = dot(C, xout) # General algorithm that interpolates U in between output points else: # Test if U has correct shape and type - legal_shapes = [(n_steps,), (1, n_steps)] if n_inputs == 1 else \ - [(n_inputs, n_steps)] - U = _check_convert_array(U, legal_shapes, - 'Parameter ``U``: ', squeeze=False, - transpose=transpose) + legal_shapes = [(n_steps,), (1, n_steps)] if n_inputs == 1 else [(n_inputs, n_steps)] + U = _check_convert_array(U, legal_shapes, 'Parameter ``U``: ', squeeze=False, transpose=transpose) # convert 1D array to 2D array with only one row if len(U.shape) == 1: - U = U.reshape(1, -1) # pylint: disable=E1103 + U = U.reshape(1, -1) # Algorithm: to integrate from time 0 to time dt, with linear # interpolation between inputs u(0) = u0 and u(dt) = u1, we solve @@ -308,12 +301,12 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False): [np.zeros((n_inputs, n_states + 2 * n_inputs))]]) expM = sp.linalg.expm(M) Ad = expM[:n_states, :n_states] - Bd1 = expM[:n_states, n_states+n_inputs:] + Bd1 = expM[:n_states, n_states + n_inputs:] Bd0 = expM[:n_states, n_states:n_states + n_inputs] - Bd1 for i in range(1, n_steps): - xout[:, i] = (dot(Ad, xout[:, i-1]) + dot(Bd0, U[:, i-1]) + - dot(Bd1, U[:, i])) + xout[:, i] = (dot(Ad, xout[:, i - 1]) + dot(Bd0, U[:, i - 1]) + dot(Bd1, U[:, i])) + yout = dot(C, xout) + dot(D, U) yout = squeeze(yout) @@ -329,13 +322,14 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False): yout = sp.transpose(yout) # See if we need to transpose the data back into MATLAB form - if (transpose): + if transpose: T = np.transpose(T) yout = np.transpose(yout) xout = np.transpose(xout) return T, yout, xout + def _get_ss_simo(sys, input=None, output=None): """Return a SISO or SIMO state-space version of sys @@ -354,6 +348,7 @@ def _get_ss_simo(sys, input=None, output=None): else: return _mimo2siso(sys_ss, input, output, warn_conversion=warn) + def step_response(sys, T=None, X0=0., input=None, output=None, transpose=False, return_x=False): # pylint: disable=W0622 @@ -424,8 +419,7 @@ def step_response(sys, T=None, X0=0., input=None, output=None, U = np.ones_like(T) - T, yout, xout = forced_response(sys, T, U, X0, - transpose=transpose) + T, yout, xout = forced_response(sys, T, U, X0, transpose=transpose) if return_x: return T, yout, xout @@ -433,8 +427,7 @@ def step_response(sys, T=None, X0=0., input=None, output=None, return T, yout -def initial_response(sys, T=None, X0=0., input=0, output=None, - transpose=False, return_x=False): +def initial_response(sys, T=None, X0=0., input=0, output=None, transpose=False, return_x=False): # pylint: disable=W0622 """Initial condition response of a linear system @@ -506,8 +499,7 @@ def initial_response(sys, T=None, X0=0., input=0, output=None, return T, yout -def impulse_response(sys, T=None, X0=0., input=0, output=None, - transpose=False, return_x=False): +def impulse_response(sys, T=None, X0=0., input=0, output=None, transpose=False, return_x=False): # pylint: disable=W0622 """Impulse response of a linear system @@ -574,8 +566,7 @@ def impulse_response(sys, T=None, X0=0., input=0, output=None, # create X0 if not given, test if X0 has correct shape. # Must be done here because it is used for computations here. n_states = sys.A.shape[0] - X0 = _check_convert_array(X0, [(n_states,), (n_states, 1)], - 'Parameter ``X0``: \n', squeeze=True) + X0 = _check_convert_array(X0, [(n_states,), (n_states, 1)], 'Parameter ``X0``: \n', squeeze=True) # Compute T and U, no checks necessary, they will be checked in lsim if T is None: diff --git a/control/xferfcn.py b/control/xferfcn.py index 44a038773..d8e6c3695 100644 --- a/control/xferfcn.py +++ b/control/xferfcn.py @@ -231,7 +231,7 @@ def __str__(self, var=None): mimo = self.inputs > 1 or self.outputs > 1 if var is None: - #! TODO: replace with standard calls to lti functions + # TODO: replace with standard calls to lti functions var = 's' if self.dt is None or self.dt == 0 else 'z' outstr = "" @@ -250,11 +250,9 @@ def __str__(self, var=None): # Center the numerator or denominator if len(numstr) < dashcount: - numstr = (' ' * int(round((dashcount - len(numstr)) / 2)) + - numstr) + numstr = (' ' * int(round((dashcount - len(numstr)) / 2)) + numstr) if len(denstr) < dashcount: - denstr = (' ' * int(round((dashcount - len(denstr)) / 2)) + - denstr) + denstr = (' ' * int(round((dashcount - len(denstr)) / 2)) + denstr) outstr += "\n" + numstr + "\n" + dashes + "\n" + denstr + "\n" @@ -291,11 +289,10 @@ def __add__(self, other): # Check that the input-output sizes are consistent. if self.inputs != other.inputs: - raise ValueError("The first summand has %i input(s), but the second has %i." - % (self.inputs, other.inputs)) + raise ValueError("The first summand has %i input(s), but the second has %i." % (self.inputs, other.inputs)) + if self.outputs != other.outputs: - raise ValueError("The first summand has %i output(s), but the second has %i." - % (self.outputs, other.outputs)) + raise ValueError("The first summand has %i output(s), but the second has %i." % (self.outputs, other.outputs)) # Figure out the sampling time to use if self.dt is None and other.dt is not None: @@ -311,9 +308,7 @@ def __add__(self, other): for i in range(self.outputs): for j in range(self.inputs): - num[i][j], den[i][j] = _add_siso(self.num[i][j], self.den[i][j], - other.num[i][j], - other.den[i][j]) + num[i][j], den[i][j] = _add_siso(self.num[i][j], self.den[i][j], other.num[i][j], other.den[i][j]) return TransferFunction(num, den, dt) @@ -333,8 +328,7 @@ def __mul__(self, other): """Multiply two LTI objects (serial connection).""" # Convert the second argument to a transfer function. if isinstance(other, (int, float, complex, np.number)): - other = _convert_to_transfer_function(other, inputs=self.inputs, - outputs=self.inputs) + other = _convert_to_transfer_function(other, inputs=self.inputs, outputs=self.inputs) else: other = _convert_to_transfer_function(other) @@ -355,12 +349,12 @@ def __mul__(self, other): raise ValueError("Systems have different sampling times") # Preallocate the numerator and denominator of the sum. - num = [[[0] for j in range(inputs)] for i in range(outputs)] - den = [[[1] for j in range(inputs)] for i in range(outputs)] + num = [[[0] for _ in range(inputs)] for _ in range(outputs)] + den = [[[1] for _ in range(inputs)] for _ in range(outputs)] # Temporary storage for the summands needed to find the (i, j)th element of the product. - num_summand = [[] for k in range(self.inputs)] - den_summand = [[] for k in range(self.inputs)] + num_summand = [[] for _ in range(self.inputs)] + den_summand = [[] for _ in range(self.inputs)] # Multiply & add. for row in range(outputs): @@ -368,9 +362,7 @@ def __mul__(self, other): for k in range(self.inputs): num_summand[k] = polymul(self.num[row][k], other.num[k][col]) den_summand[k] = polymul(self.den[row][k], other.den[k][col]) - num[row][col], den[row][col] = _add_siso( - num[row][col], den[row][col], - num_summand[k], den_summand[k]) + num[row][col], den[row][col] = _add_siso(num[row][col], den[row][col], num_summand[k], den_summand[k]) return TransferFunction(num, den, dt) @@ -427,17 +419,12 @@ def __truediv__(self, other): """Divide two LTI objects.""" if isinstance(other, (int, float, complex, np.number)): - other = _convert_to_transfer_function( - other, inputs=self.inputs, - outputs=self.inputs) + other = _convert_to_transfer_function(other, inputs=self.inputs, outputs=self.inputs) else: other = _convert_to_transfer_function(other) - if (self.inputs > 1 or self.outputs > 1 or - other.inputs > 1 or other.outputs > 1): - raise NotImplementedError( - "TransferFunction.__truediv__ is currently \ - implemented only for SISO systems.") + if self.inputs > 1 or self.outputs > 1 or other.inputs > 1 or other.outputs > 1: + raise NotImplementedError("TransferFunction.__truediv__ is currently implemented only for SISO systems.") # Figure out the sampling time to use if self.dt is None and other.dt is not None: @@ -466,10 +453,8 @@ def __rtruediv__(self, other): else: other = _convert_to_transfer_function(other) - if (self.inputs > 1 or self.outputs > 1 or - other.inputs > 1 or other.outputs > 1): - raise NotImplementedError( - "TransferFunction.__rtruediv__ is currently implemented only for SISO systems.") + if self.inputs > 1 or self.outputs > 1 or other.inputs > 1 or other.outputs > 1: + raise NotImplementedError("TransferFunction.__rtruediv__ is currently implemented only for SISO systems.") return other / self @@ -567,8 +552,7 @@ def horner(self, s): for i in range(self.outputs): for j in range(self.inputs): - out[i][j] = (polyval(self.num[i][j], s) / - polyval(self.den[i][j], s)) + out[i][j] = (polyval(self.num[i][j], s) / polyval(self.den[i][j], s)) return out @@ -602,8 +586,7 @@ def freqresp(self, omega): # Compute frequency response for each input/output pair for i in range(self.outputs): for j in range(self.inputs): - fresp = (polyval(self.num[i][j], slist) / - polyval(self.den[i][j], slist)) + fresp = (polyval(self.num[i][j], slist) / polyval(self.den[i][j], slist)) mag[i, j, :] = abs(fresp) phase[i, j, :] = angle(fresp) @@ -620,8 +603,7 @@ def pole(self): def zero(self): """Compute the zeros of a transfer function.""" if self.inputs > 1 or self.outputs > 1: - raise NotImplementedError("TransferFunction.zero is currently only implemented " - "for SISO systems.") + raise NotImplementedError("TransferFunction.zero is currently only implemented for SISO systems.") else: # for now, just give zeros of a SISO tf return roots(self.num[0][0]) @@ -630,11 +612,9 @@ def feedback(self, other=1, sign=-1): """Feedback interconnection between two LTI objects.""" other = _convert_to_transfer_function(other) - if (self.inputs > 1 or self.outputs > 1 or - other.inputs > 1 or other.outputs > 1): + if self.inputs > 1 or self.outputs > 1 or other.inputs > 1 or other.outputs > 1: # TODO: MIMO feedback - raise NotImplementedError("TransferFunction.feedback is currently only implemented " - "for SISO functions.") + raise NotImplementedError("TransferFunction.feedback is currently only implemented for SISO functions.") # Figure out the sampling time to use if self.dt is None and other.dt is not None: @@ -668,8 +648,8 @@ def minreal(self, tol=None): sqrt_eps = sqrt(float_info.epsilon) # pre-allocate arrays - num = [[[] for j in range(self.inputs)] for i in range(self.outputs)] - den = [[[] for j in range(self.inputs)] for i in range(self.outputs)] + num = [[[] for _ in range(self.inputs)] for _ in range(self.outputs)] + den = [[[] for _ in range(self.inputs)] for _ in range(self.outputs)] for i in range(self.outputs): for j in range(self.inputs): @@ -682,8 +662,7 @@ def minreal(self, tol=None): # check all zeros for z in zeros: - t = tol or \ - 1000 * max(float_info.epsilon, abs(z) * sqrt_eps) + t = tol or 1000 * max(float_info.epsilon, abs(z) * sqrt_eps) idx = where(abs(z - poles) < t)[0] if len(idx): # cancel this zero against one of the poles @@ -714,11 +693,10 @@ def returnScipySignalLTI(self): # TODO: implement for discrete time systems if self.dt != 0 and self.dt is not None: - raise NotImplementedError("Function not \ - implemented in discrete time") + raise NotImplementedError("Function not implemented in discrete time") # Preallocate the output. - out = [[[] for j in range(self.inputs)] for i in range(self.outputs)] + out = [[[] for _ in range(self.inputs)] for _ in range(self.outputs)] for i in range(self.outputs): for j in range(self.inputs): @@ -742,7 +720,7 @@ def _common_den(self, imag_tol=None): ---------- imag_tol: float Threshold for the imaginary part of a root to use in detecting - complex poles + complex poles. NOT USED!! Returns ------- @@ -773,12 +751,11 @@ def _common_den(self, imag_tol=None): eps = finfo(float).eps # Decide on the tolerance to use in deciding of a pole is complex - if (imag_tol is None): - imag_tol = 1e-8 # TODO: figure out the right number to use + # TODO: figure out the right imag_tol to use. # A list to keep track of cumulative poles found as we scan # self.den[..][..] - poles = [[] for j in range(self.inputs)] + poles = [[] for _ in range(self.inputs)] # RvP, new implementation 180526, issue #194 @@ -792,8 +769,7 @@ def _common_den(self, imag_tol=None): poleset.append([]) for j in range(self.inputs): if abs(self.num[i][j]).max() <= eps: - poleset[-1].append([array([], dtype=float), - roots(self.den[i][j]), 0.0, [], 0]) + poleset[-1].append([array([], dtype=float), roots(self.den[i][j]), 0.0, [], 0]) else: z, p, k = tf2zpk(self.num[i][j], self.den[i][j]) poleset[-1].append([z, p, k, [], 0]) @@ -805,8 +781,7 @@ def _common_den(self, imag_tol=None): currentpoles = poleset[i][j][1] nothave = ones(currentpoles.shape, dtype=bool) for ip, p in enumerate(poles[j]): - idx, = nonzero( - (abs(currentpoles - p) < epsnm) * nothave) + idx, = nonzero((abs(currentpoles - p) < epsnm) * nothave) if len(idx): nothave[idx[0]] = False else: @@ -821,8 +796,7 @@ def _common_den(self, imag_tol=None): # figure out maximum number of poles, for sizing the den npmax = max([len(p) for p in poles]) den = zeros((self.inputs, npmax + 1), dtype=float) - num = zeros((max(1, self.outputs, self.inputs), - max(1, self.outputs, self.inputs), npmax + 1), dtype=float) + num = zeros((max(1, self.outputs, self.inputs), max(1, self.outputs, self.inputs), npmax + 1), dtype=float) denorder = zeros((self.inputs,), dtype=int) for j in range(self.inputs): @@ -845,8 +819,7 @@ def _common_den(self, imag_tol=None): nwzeros = list(poleset[i][j][0]) # add all poles not found in the original denominator, # and the ones later added from other denominators - for ip in chain(poleset[i][j][3], - range(poleset[i][j][4], np)): + for ip in chain(poleset[i][j][3], range(poleset[i][j][4], np)): nwzeros.append(poles[j][ip]) numpoly = poleset[i][j][2] * polyfromroots(nwzeros).real @@ -858,8 +831,8 @@ def _common_den(self, imag_tol=None): # print(num[i, j]) if (abs(den.imag) > epsnm).any(): - print("Warning: The denominator has a nontrivial imaginary part: %f" - % abs(den.imag).max()) + print("Warning: The denominator has a nontrivial imaginary part: %f" % abs(den.imag).max()) + den = den.real return num, den, denorder @@ -906,10 +879,13 @@ def sample(self, Ts, method='zoh', alpha=None): """ if not self.isctime(): raise ValueError("System must be continuous time system") + if not self.issiso(): raise NotImplementedError("MIMO implementation not available") + if method == "matched": return _c2d_matched(self, Ts) + sys = (self.num[0][0], self.den[0][0]) numd, dend, dt = cont2discrete(sys, Ts, method, alpha) return TransferFunction(numd[0, :], dend, dt) @@ -948,6 +924,7 @@ def _dcgain_cont(self): else: # numerator is zero too: give up gain[i][j] = np.nan + return np.squeeze(gain) # c2d function contributed by Benjamin White, Oct 2012 @@ -960,19 +937,23 @@ def _c2d_matched(sysC, Ts): zpoles = [0] * len(spoles) pregainnum = [0] * len(szeros) pregainden = [0] * len(spoles) + for idx, s in enumerate(szeros): sTs = s * Ts z = exp(sTs) zzeros[idx] = z pregainnum[idx] = 1 - z + for idx, s in enumerate(spoles): sTs = s * Ts z = exp(sTs) zpoles[idx] = z pregainden[idx] = 1 - z + zgain = np.multiply.reduce(pregainnum) / np.multiply.reduce(pregainden) gain = sgain / zgain sysDnum, sysDden = zpk2tf(zzeros, zpoles, gain) + return TransferFunction(sysDnum, sysDden, Ts) # Utility function to convert a transfer function polynomial to a string @@ -1025,6 +1006,7 @@ def _tf_polynomial_to_string(coeffs, var='s'): thestr = "-%s" % (newstr,) else: thestr = newstr + return thestr @@ -1069,8 +1051,7 @@ def _convert_to_transfer_function(sys, **kw): if isinstance(sys, TransferFunction): if len(kw): - raise TypeError("If sys is a TransferFunction, " + - "_convertToTransferFunction cannot take keywords.") + raise TypeError("If sys is a TransferFunction, _convertToTransferFunction cannot take keywords.") return sys elif isinstance(sys, StateSpace): @@ -1085,9 +1066,7 @@ def _convert_to_transfer_function(sys, **kw): try: from slycot import tb04ad if len(kw): - raise TypeError( - "If sys is a StateSpace, " + - "_convertToTransferFunction cannot take keywords.") + raise TypeError("If sys is a StateSpace, _convertToTransferFunction cannot take keywords.") # Use Slycot to make the transformation # Make sure to convert system matrices to numpy arrays @@ -1095,8 +1074,8 @@ def _convert_to_transfer_function(sys, **kw): array(sys.B), array(sys.C), array(sys.D), tol1=0.0) # Preallocate outputs. - num = [[[] for j in range(sys.inputs)] for i in range(sys.outputs)] - den = [[[] for j in range(sys.inputs)] for i in range(sys.outputs)] + num = [[[] for _ in range(sys.inputs)] for _ in range(sys.outputs)] + den = [[[] for _ in range(sys.inputs)] for _ in range(sys.outputs)] for i in range(sys.outputs): for j in range(sys.inputs): @@ -1141,8 +1120,7 @@ def _convert_to_transfer_function(sys, **kw): den = [[[1] for j in range(inputs)] for i in range(outputs)] return TransferFunction(num, den) except Exception as e: - print("Failure to assume argument is matrix-like in" - " _convertToTransferFunction, result %s" % e) + print("Failure to assume argument is matrix-like in _convertToTransferFunction, result %s" % e) raise TypeError("Can't convert given type to TransferFunction system.") @@ -1235,8 +1213,7 @@ def tf(*args): elif isinstance(sys, TransferFunction): return deepcopy(sys) else: - raise TypeError("tf(sys): sys must be a StateSpace or TransferFunction object. " - "It is %s." % type(sys)) + raise TypeError("tf(sys): sys must be a StateSpace or TransferFunction object. It is %s." % type(sys)) else: raise ValueError("Needs 1 or 2 arguments; received %i." % len(args)) @@ -1353,22 +1330,19 @@ def _clean_part(data): valid_types = (int, float, complex, np.number) valid_collection = (list, tuple, ndarray) - if (isinstance(data, valid_types) or - (isinstance(data, ndarray) and data.ndim == 0)): + if isinstance(data, valid_types) or (isinstance(data, ndarray) and data.ndim == 0): # Data is a scalar (including 0d ndarray) data = [[array([data])]] - elif (isinstance(data, ndarray) and data.ndim == 3 and - isinstance(data[0, 0, 0], valid_types)): - data = [[array(data[i, j]) - for j in range(data.shape[1])] - for i in range(data.shape[0])] - elif (isinstance(data, valid_collection) and - all([isinstance(d, valid_types) for d in data])): + elif isinstance(data, ndarray) and \ + data.ndim == 3 and \ + isinstance(data[0, 0, 0], valid_types): + data = [[array(data[i, j]) for j in range(data.shape[1])] for i in range(data.shape[0])] + elif isinstance(data, valid_collection) and \ + all([isinstance(d, valid_types) for d in data]): data = [[array(data)]] - elif (isinstance(data, (list, tuple)) and - isinstance(data[0], (list, tuple)) and - (isinstance(data[0][0], valid_collection) and - all([isinstance(d, valid_types) for d in data[0][0]]))): + elif isinstance(data, (list, tuple)) and \ + isinstance(data[0], (list, tuple)) and \ + (isinstance(data[0][0], valid_collection) and all([isinstance(d, valid_types) for d in data[0][0]])): data = list(data) for j in range(len(data)): data[j] = list(data[j]) diff --git a/setup.cfg b/setup.cfg index 3c6e79cf3..73dbcef01 100644 --- a/setup.cfg +++ b/setup.cfg @@ -1,2 +1,9 @@ [bdist_wheel] universal=1 + + +[flake8] +max-line-length = 160 +# ignore E402 (import not at top of file) and W605 (invalid escape sequence) +ignore = E402, W605 +exclude = .git,control/tests/*,build,dist,__pycache__ \ No newline at end of file diff --git a/setup.py b/setup.py index cd4bcbf9f..dd311fb3c 100644 --- a/setup.py +++ b/setup.py @@ -44,6 +44,7 @@ 'matplotlib'], tests_require=['scipy', 'matplotlib', - 'nose'], + 'nose', + 'flake8'], test_suite = 'nose.collector', ) From c6d21f2d9bb229e57dd4c8ce05d5d4e0e6dc3d89 Mon Sep 17 00:00:00 2001 From: Joris Geysens Date: Wed, 27 Mar 2019 13:38:03 +0100 Subject: [PATCH 2/6] Install flake8 manually for travis CI. Run flake8 in before_script. --- .travis.yml | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/.travis.yml b/.travis.yml index b7eb18ba1..fa2e84606 100644 --- a/.travis.yml +++ b/.travis.yml @@ -62,7 +62,7 @@ before_install: # Install packages install: # Install packages needed by python-control - - conda install $SCIPY matplotlib + - conda install $SCIPY matplotlib flake8 # Build slycot from source # For python 3, need to provide pointer to python library #! git clone https://github.com/repagh/Slycot.git slycot; @@ -71,9 +71,11 @@ install: cd slycot; python setup.py install; cd ..; fi +before_script: + - flake8 control/ + # command to run tests script: - - flake8 control/ - 'if [ $SLYCOT != "" ]; then python -c "import slycot"; fi' - coverage run setup.py test From 38327007b853801a2fc1cb4ba5174da5a1e0be71 Mon Sep 17 00:00:00 2001 From: Joris Geysens Date: Wed, 27 Mar 2019 17:30:03 +0100 Subject: [PATCH 3/6] Fix remaining flake8 issues. --- control/frdata.py | 22 +++++++++++----------- control/freqplot.py | 4 ++-- control/lti.py | 6 +++++- control/mateqn.py | 1 - control/modelsimp.py | 1 - 5 files changed, 18 insertions(+), 16 deletions(-) diff --git a/control/frdata.py b/control/frdata.py index 24ad2f51e..d4c428c24 100644 --- a/control/frdata.py +++ b/control/frdata.py @@ -368,7 +368,7 @@ def _evalfr(self, omega): if self.ifunc is None: try: out = self.fresp[:, :, where(self.omega == omega)[0][0]] - except: + except IndexError: raise ValueError("Frequency %f not in frequency list, try an interpolating FRD if you want additional points" % omega) else: if getattr(omega, '__iter__', False): @@ -469,18 +469,18 @@ def _convertToFRD(sys, omega, inputs=1, outputs=1): return FRD(fresp, omega, smooth=True) # try converting constant matrices + + sys = array(sys) + outputs, inputs = sys.shape + fresp = empty((outputs, inputs, len(omega)), dtype=float) + for i in range(outputs): + for j in range(inputs): + fresp[i, j, :] = sys[i, j] + try: - sys = array(sys) - outputs, inputs = sys.shape - fresp = empty((outputs, inputs, len(omega)), dtype=float) - for i in range(outputs): - for j in range(inputs): - fresp[i, j, :] = sys[i, j] return FRD(fresp, omega, smooth=True) - except: - pass - - raise TypeError('''Can't convert given type "%s" to FRD system.''' % sys.__class__) + except(TypeError, ValueError): + raise TypeError('''Can't convert given type "%s" to FRD system.''' % sys.__class__) def frd(*args, **kwargs): diff --git a/control/freqplot.py b/control/freqplot.py index 15d7ff763..07e743833 100644 --- a/control/freqplot.py +++ b/control/freqplot.py @@ -524,7 +524,7 @@ def default_frequency_range(syslist, Hz=None, number_of_samples=None, Parameters ---------- - syslist : list of LTI + syslist : tuple 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 @@ -598,7 +598,7 @@ def default_frequency_range(syslist, Hz=None, number_of_samples=None, else: # TODO raise NotImplementedError('type of system in not implemented now') - except: + except NotImplementedError: pass # Make sure there is at least one point in the range diff --git a/control/lti.py b/control/lti.py index 9f6eb5da3..9037fa4da 100644 --- a/control/lti.py +++ b/control/lti.py @@ -170,7 +170,11 @@ def timebaseEqual(sys1, sys2): if isinstance(sys1.dt, bool) or isinstance(sys2.dt, bool): # Make sure both are unspecified discrete timebases - return type(sys1.dt) == type(sys2.dt) and sys1.dt == sys2.dt + if (isinstance(sys1.dt, bool) and not isinstance(sys2.dt, bool)) or \ + (not isinstance(sys1.dt, bool) and isinstance(sys2.dt, bool)): + return False + + return sys1.dt == sys2.dt elif sys1.dt is None or sys2.dt is None: # One or the other is unspecified => the other can be anything return True diff --git a/control/mateqn.py b/control/mateqn.py index 9163e63f9..2bf818521 100644 --- a/control/mateqn.py +++ b/control/mateqn.py @@ -873,7 +873,6 @@ def dare_old(A, B, Q, R, S=None, E=None): A_b = copy(A) R_b = copy(R) B_b = copy(B) - E_b = copy(E) S_b = copy(S) # Solve the generalized algebraic Riccati equation by calling the diff --git a/control/modelsimp.py b/control/modelsimp.py index d163e89ad..ff1714387 100644 --- a/control/modelsimp.py +++ b/control/modelsimp.py @@ -410,7 +410,6 @@ def markov(Y, U, M): """ # Convert input parameters to matrices (if they aren't already) - Ymat = np.mat(Y) Umat = np.mat(U) n = np.size(U) From b769d3d52ad9edd11bf738c2f90e5584c1972893 Mon Sep 17 00:00:00 2001 From: Joris Geysens Date: Wed, 27 Mar 2019 17:50:27 +0100 Subject: [PATCH 4/6] Fix python 3 flake8 errors. --- control/margins.py | 5 ++--- control/modelsimp.py | 2 +- control/phaseplot.py | 2 +- 3 files changed, 4 insertions(+), 5 deletions(-) diff --git a/control/margins.py b/control/margins.py index 5f756d273..b606f8257 100644 --- a/control/margins.py +++ b/control/margins.py @@ -148,9 +148,8 @@ def stability_margins(sysdata, returnall=False, epsw=0.0): else: sys = xferfcn._convert_to_transfer_function(sysdata) except Exception as e: - print (e) - raise ValueError("Margin sysdata must be either a linear system or " - "a 3-sequence of mag, phase, omega.") + print(e) + raise ValueError("Margin sysdata must be either a linear system or a 3-sequence of mag, phase, omega.") # calculate gain of system if isinstance(sys, xferfcn.TransferFunction): diff --git a/control/modelsimp.py b/control/modelsimp.py index ff1714387..2c05f50e6 100644 --- a/control/modelsimp.py +++ b/control/modelsimp.py @@ -290,7 +290,7 @@ def balred(sys, orders, method='truncate', alpha=None): # check if orders is a list or a scalar try: - _ = iter(orders) + iter(orders) except TypeError: # if orders is a scalar orders = [orders] diff --git a/control/phaseplot.py b/control/phaseplot.py index e68a98e80..c36405a11 100644 --- a/control/phaseplot.py +++ b/control/phaseplot.py @@ -295,7 +295,7 @@ def phase_plot(odefun, X=None, Y=None, scale=1, X0=None, T=None, # set(xy, 'AutoScaleFactor', 0); if scale < 0: - _ = mpl.plot(x1, x2, 'b.') # add dots at base + mpl.plot(x1, x2, 'b.') # add dots at base # set(bp, 'MarkerSize', PP_arrow_markersize); return From c0491b7ea5897aa214b1275a334f3b7eb4bfe1b9 Mon Sep 17 00:00:00 2001 From: Joris Geysens Date: Thu, 28 Mar 2019 10:31:18 +0100 Subject: [PATCH 5/6] Python 3.5 flake8 error fix --- control/phaseplot.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/control/phaseplot.py b/control/phaseplot.py index c36405a11..23607dbbf 100644 --- a/control/phaseplot.py +++ b/control/phaseplot.py @@ -289,7 +289,7 @@ def phase_plot(odefun, X=None, Y=None, scale=1, X0=None, T=None, mpl.quiver(x1, x2, dx[:, :, 0], dx[:, :, 1], angles='xy') elif scale != 0 and Narrows > 0: # TODO: figure out arguments to make arrows show up correctly - _ = mpl.quiver(x1, x2, dx[:, :, 0] * abs(scale), dx[:, :, 1] * abs(scale), angles='xy') + mpl.quiver(x1, x2, dx[:, :, 0] * abs(scale), dx[:, :, 1] * abs(scale), angles='xy') # set(xy, 'LineWidth', PP_arrow_linewidth); # set(xy, 'AutoScale', 'off'); # set(xy, 'AutoScaleFactor', 0); From abb517178a795e4d73021fac5bf41efab90b9b08 Mon Sep 17 00:00:00 2001 From: Joris Geysens Date: Thu, 28 Mar 2019 10:31:18 +0100 Subject: [PATCH 6/6] Python 3.5 flake8 error fix --- control/phaseplot.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/control/phaseplot.py b/control/phaseplot.py index c36405a11..b4036022d 100644 --- a/control/phaseplot.py +++ b/control/phaseplot.py @@ -177,7 +177,7 @@ def phase_plot(odefun, X=None, Y=None, scale=1, X0=None, T=None, elif scale != 0: # TODO: optimize parameters for arrows # TODO: figure out arguments to make arrows show up correctly - _ = mpl.quiver(x1, x2, dx[:, :, 0] * np.abs(scale), dx[:, :, 1] * np.abs(scale), angles='xy') + mpl.quiver(x1, x2, dx[:, :, 0] * np.abs(scale), dx[:, :, 1] * np.abs(scale), angles='xy') # set(xy, 'LineWidth', PP_arrow_linewidth, 'Color', 'b') # TODO: Tweak the shape of the plot @@ -289,7 +289,7 @@ def phase_plot(odefun, X=None, Y=None, scale=1, X0=None, T=None, mpl.quiver(x1, x2, dx[:, :, 0], dx[:, :, 1], angles='xy') elif scale != 0 and Narrows > 0: # TODO: figure out arguments to make arrows show up correctly - _ = mpl.quiver(x1, x2, dx[:, :, 0] * abs(scale), dx[:, :, 1] * abs(scale), angles='xy') + mpl.quiver(x1, x2, dx[:, :, 0] * abs(scale), dx[:, :, 1] * abs(scale), angles='xy') # set(xy, 'LineWidth', PP_arrow_linewidth); # set(xy, 'AutoScale', 'off'); # set(xy, 'AutoScaleFactor', 0);