diff --git a/control/bdalg.py b/control/bdalg.py index b92496477..cd6f59068 100644 --- a/control/bdalg.py +++ b/control/bdalg.py @@ -247,7 +247,8 @@ def feedback(sys1, sys2=1, sign=-1): return sys1.feedback(sys2, sign) def append(*sys): - ''' + '''append(sys1, sys2, ... sysn) + Group models by appending their inputs and outputs Forms an augmented system model, and appends the inputs and diff --git a/control/ctrlutil.py b/control/ctrlutil.py index 273257c74..35035c7e9 100644 --- a/control/ctrlutil.py +++ b/control/ctrlutil.py @@ -43,12 +43,12 @@ # Packages that we need access to from . import lti import numpy as np -from numpy import pi +import math __all__ = ['unwrap', 'issys', 'db2mag', 'mag2db'] # Utility function to unwrap an angle measurement -def unwrap(angle, period=2*pi): +def unwrap(angle, period=2*math.pi): """Unwrap a phase angle to give a continuous curve Parameters diff --git a/control/frdata.py b/control/frdata.py index 313da12b2..32ba8f11b 100644 --- a/control/frdata.py +++ b/control/frdata.py @@ -58,7 +58,9 @@ __all__ = ['FRD', 'frd'] class FRD(LTI): - """A class for models defined by Frequency Response Data (FRD) + """FRD(d, w) + + A class for models defined by frequency response data (FRD) The FRD class is used to represent systems in frequency response data form. @@ -81,7 +83,9 @@ class FRD(LTI): epsw = 1e-8 def __init__(self, *args, **kwargs): - """Construct an FRD object + """FRD(d, w) + + Construct an FRD object The default constructor is FRD(d, w), where w is an iterable of frequency points, and d is the matching frequency data. @@ -470,8 +474,9 @@ def _convertToFRD(sys, omega, inputs=1, outputs=1): sys.__class__) def frd(*args): - """ - Construct a Frequency Response Data model, or convert a system + """frd(d, w) + + Construct a frequency response data model frd models store the (measured) frequency response of a system. @@ -501,6 +506,6 @@ def frd(*args): See Also -------- - ss, tf + FRD, ss, tf """ return FRD(*args) diff --git a/control/freqplot.py b/control/freqplot.py index 315a9635b..c53e83f31 100644 --- a/control/freqplot.py +++ b/control/freqplot.py @@ -44,6 +44,7 @@ import matplotlib.pyplot as plt import scipy as sp import numpy as np +import math from .ctrlutil import unwrap from .bdalg import feedback @@ -128,7 +129,7 @@ def bode_plot(syslist, omega=None, dB=None, Hz=None, deg=None, else: omega_limits = np.array(omega_limits) if Hz: - omega_limits *= 2.*np.pi + omega_limits *= 2.*math.pi if omega_num: omega = sp.logspace(np.log10(omega_limits[0]), np.log10(omega_limits[1]), num=omega_num, endpoint=True) else: @@ -142,7 +143,7 @@ def bode_plot(syslist, omega=None, dB=None, Hz=None, deg=None, else: omega_sys = np.array(omega) if sys.isdtime(True): - nyquistfrq = 2. * np.pi * 1. / sys.dt / 2. + nyquistfrq = 2. * math.pi * 1. / sys.dt / 2. omega_sys = omega_sys[omega_sys < nyquistfrq] # TODO: What distance to the Nyquist frequency is appropriate? else: @@ -154,9 +155,9 @@ def bode_plot(syslist, omega=None, dB=None, Hz=None, deg=None, phase = unwrap(phase) nyquistfrq_plot = None if Hz: - omega_plot = omega_sys / (2. * np.pi) + omega_plot = omega_sys / (2. * math.pi) if nyquistfrq: - nyquistfrq_plot = nyquistfrq / (2. * np.pi) + nyquistfrq_plot = nyquistfrq / (2. * math.pi) else: omega_plot = omega_sys if nyquistfrq: @@ -187,7 +188,7 @@ def bode_plot(syslist, omega=None, dB=None, Hz=None, deg=None, # Phase plot ax_phase = plt.subplot(212, sharex=ax_mag); if deg: - phase_plot = phase * 180. / np.pi + phase_plot = phase * 180. / math.pi else: phase_plot = phase ax_phase.semilogx(omega_plot, phase_plot, *args, **kwargs) @@ -208,8 +209,8 @@ def genZeroCenteredSeries(val_min, val_max, period): ax_phase.set_yticks(genZeroCenteredSeries(ylim[0], ylim[1], 15.), minor=True) else: ylim = ax_phase.get_ylim() - ax_phase.set_yticks(genZeroCenteredSeries(ylim[0], ylim[1], np.pi / 4.)) - ax_phase.set_yticks(genZeroCenteredSeries(ylim[0], ylim[1], np.pi / 12.), minor=True) + ax_phase.set_yticks(genZeroCenteredSeries(ylim[0], ylim[1], math.pi / 4.)) + ax_phase.set_yticks(genZeroCenteredSeries(ylim[0], ylim[1], math.pi / 12.), minor=True) ax_phase.grid(True, which='both') # ax_mag.grid(which='minor', alpha=0.3) # ax_mag.grid(which='major', alpha=0.9) @@ -449,7 +450,7 @@ def default_frequency_range(syslist, Hz=None, number_of_samples=None, feature_pe features_ = features_[features_ != 0.0]; features = np.concatenate((features, features_)) elif sys.isdtime(strict=True): - fn = np.pi * 1. / sys.dt + fn = math.pi * 1. / sys.dt # TODO: What distance to the Nyquist frequency is appropriate? freq_interesting.append(fn * 0.9) @@ -475,12 +476,12 @@ def default_frequency_range(syslist, Hz=None, number_of_samples=None, feature_pe features = np.array([1.]); if Hz: - features /= 2.*np.pi + features /= 2.*math.pi features = np.log10(features) lsp_min = np.floor(np.min(features) - feature_periphery_decade) lsp_max = np.ceil(np.max(features) + feature_periphery_decade) - lsp_min += np.log10(2.*np.pi) - lsp_max += np.log10(2.*np.pi) + lsp_min += np.log10(2.*math.pi) + lsp_max += np.log10(2.*math.pi) else: features = np.log10(features) lsp_min = np.floor(np.min(features) - feature_periphery_decade) diff --git a/control/margins.py b/control/margins.py index 26ec47325..7f4f06c23 100644 --- a/control/margins.py +++ b/control/margins.py @@ -48,21 +48,21 @@ Date: 14 July 2011 $Id$ - """ +import math import numpy as np +import scipy as sp from . import xferfcn from .lti import issiso from . import frdata -import scipy as sp __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 - imaginary part with w applied""" + imaginary part with w applied""" rpencil = np.zeros_like(pol) ipencil = np.zeros_like(pol) rpencil[-1::-4] = 1. @@ -141,7 +141,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 * np.pi/180), + sys = frdata.FRD(mag * np.exp(1j * phase * math.pi/180), omega, smooth=True) else: sys = xferfcn._convertToTransferFunction(sysdata) @@ -294,8 +294,7 @@ def dstab(w): # Contributed by Steffen Waldherr #! TODO - need to add test functions def phase_crossover_frequencies(sys): - """ - Compute frequencies and gains at intersections with real axis + """Compute frequencies and gains at intersections with real axis in Nyquist plot. Call as: @@ -338,11 +337,13 @@ def phase_crossover_frequencies(sys): def margin(*args): - """Calculate gain and phase margins and associated crossover frequencies + """margin(sysdata) + + Calculate gain and phase margins and associated crossover frequencies Parameters ---------- - sysdata: LTI system or (mag, phase, omega) sequence + sysdata : LTI system or (mag, phase, omega) sequence sys : StateSpace or TransferFunction Linear SISO system mag, phase, omega : sequence of array_like @@ -360,8 +361,8 @@ def margin(*args): Wcp : float Phase crossover frequency (corresponding to gain margin) (in rad/sec) - Margins are of SISO open-loop. If more than one crossover frequency is - detected, returns the lowest corresponding margin. + Margins are of SISO open-loop. If more than one crossover frequency is + detected, returns the lowest corresponding margin. Examples -------- diff --git a/control/matlab/wrappers.py b/control/matlab/wrappers.py index 8c1cfb005..d83890a33 100644 --- a/control/matlab/wrappers.py +++ b/control/matlab/wrappers.py @@ -10,7 +10,9 @@ __all__ = ['bode', 'ngrid', 'dcgain'] def bode(*args, **keywords): - """Bode plot of the frequency response + """bode(syslist[, omega, dB, Hz, deg, ...]) + + Bode plot of the frequency response Plots a bode gain and phase diagram diff --git a/control/robust.py b/control/robust.py index 1b33e93d0..0b9b98c41 100644 --- a/control/robust.py +++ b/control/robust.py @@ -72,7 +72,6 @@ def h2syn(P,nmeas,ncon): >>> K = h2syn(P,nmeas,ncon) """ - #Check for ss system object, need a utility for this? #TODO: Check for continous or discrete, only continuous supported right now @@ -116,11 +115,11 @@ def hinfsyn(P,nmeas,ncon): CL: closed loop system (State-space sys) gam: infinity norm of closed loop system rcond: 4-vector, reciprocal condition estimates of: - 1: control transformation matrix - 2: measurement transformation matrix - 3: X-Ricatti equation - 4: Y-Ricatti equation - TODO: document significance of rcond + 1: control transformation matrix + 2: measurement transformation matrix + 3: X-Ricatti equation + 4: Y-Ricatti equation + TODO: document significance of rcond Raises ------ diff --git a/control/statefbk.py b/control/statefbk.py index 977f83ae4..634922131 100644 --- a/control/statefbk.py +++ b/control/statefbk.py @@ -144,7 +144,9 @@ def acker(A, B, poles): return K def lqr(*args, **keywords): - """Linear quadratic regulator design + """lqr(A, B, Q, R[, N]) + + Linear quadratic regulator design The lqr() function computes the optimal state feedback controller that minimizes the quadratic cost diff --git a/control/statesp.py b/control/statesp.py index 816485f77..4abc4175f 100644 --- a/control/statesp.py +++ b/control/statesp.py @@ -51,9 +51,10 @@ $Id$ """ +import math import numpy as np from numpy import all, angle, any, array, asarray, concatenate, cos, delete, \ - dot, empty, exp, eye, matrix, ones, pi, poly, poly1d, roots, shape, sin, \ + dot, empty, exp, eye, matrix, ones, poly, poly1d, roots, shape, sin, \ zeros, squeeze from numpy.random import rand, randn from numpy.linalg import solve, eigvals, matrix_rank @@ -69,7 +70,9 @@ __all__ = ['StateSpace', 'ss', 'rss', 'drss', 'tf2ss', 'ssdata'] class StateSpace(LTI): - """A class for representing state-space models + """StateSpace(A, B, C, D[, dt]) + + A class for representing state-space models The StateSpace class is used to represent state-space realizations of linear time-invariant (LTI) systems: @@ -89,15 +92,19 @@ class StateSpace(LTI): means the system timebase is not specified. If 'dt' is set to True, the system will be treated as a discrete time system with unspecified sampling time. - """ def __init__(self, *args): - """Construct a state space object. + """ + StateSpace(A, B, C, D[, dt]) + + Construct a state space object. - The default constructor is StateSpace(A, B, C, D), where A, B, C, D are - matrices or equivalent objects. To call the copy constructor, call - StateSpace(sys), where sys is a StateSpace object. + The default constructor is StateSpace(A, B, C, D), where A, B, C, D + are matrices or equivalent objects. To create a discrete time system, + use StateSpace(A, B, C, D, dt) where 'dt' is the sampling time (or + True for unspecified sampling time). To call the copy constructor, + call StateSpace(sys), where sys is a StateSpace object. """ @@ -111,8 +118,7 @@ def __init__(self, *args): elif len(args) == 1: # Use the copy constructor. if not isinstance(args[0], StateSpace): - raise TypeError("The one-argument constructor can only take in \ -a StateSpace object. Recived %s." % type(args[0])) + raise TypeError("The one-argument constructor can only take in a StateSpace object. Received %s." % type(args[0])) A = args[0].A B = args[0].B C = args[0].C @@ -363,7 +369,7 @@ def evalfr(self, omega): if isdtime(self, strict=True): dt = timebase(self) s = exp(1.j * omega * dt) - if (omega * dt > pi): + if (omega * dt > math.pi): warnings.warn("evalfr: frequency evaluation above Nyquist frequency") else: s = omega * 1.j @@ -793,7 +799,7 @@ def _rss_generate(states, inputs, outputs, type): poles[i] = complex(-exp(randn()), 3. * exp(randn())) elif type == 'd': mag = rand() - phase = 2. * pi * 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) @@ -956,7 +962,8 @@ def _mimo2simo(sys, input, warn_conversion=False): return sys def ss(*args): - """ + """ss(A, B, C, D[, dt]) + Create a state space system. The function accepts either 1, 4 or 5 parameters: @@ -1014,6 +1021,7 @@ def ss(*args): See Also -------- + StateSpace tf ss2tf tf2ss @@ -1045,7 +1053,8 @@ def ss(*args): raise ValueError("Needs 1 or 4 arguments; received %i." % len(args)) def tf2ss(*args): - """ + """tf2ss(sys) + Transform a transfer function to a state space system. The function accepts either 1 or 2 parameters: diff --git a/control/statesp.py+ b/control/statesp.py+ new file mode 100644 index 000000000..1007feb0d --- /dev/null +++ b/control/statesp.py+ @@ -0,0 +1,1211 @@ +"""statesp.py + +State space representation and functions. + +This file contains the StateSpace class, which is used to represent linear +systems in state space. This is the primary representation for the +python-control library. + +""" + +# Python 3 compatibility (needs to go here) +from __future__ import print_function +from __future__ import division # for _convertToStateSpace + +"""Copyright (c) 2010 by California Institute of Technology +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions +are met: + +1. Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + +2. Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in the + documentation and/or other materials provided with the distribution. + +3. Neither the name of the California Institute of Technology nor + the names of its contributors may be used to endorse or promote + products derived from this software without specific prior + written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS +FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL CALTECH +OR THE CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF +USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND +ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, +OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT +OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF +SUCH DAMAGE. + +Author: Richard M. Murray +Date: 24 May 09 +Revised: Kevin K. Chen, Dec 10 + +$Id$ +""" + +import numpy as np +from numpy import all, angle, any, array, asarray, concatenate, cos, delete, \ + dot, empty, exp, eye, matrix, ones, pi, poly, poly1d, roots, shape, sin, \ + zeros, squeeze +from numpy.random import rand, randn +from numpy.linalg import solve, eigvals, matrix_rank +from numpy.linalg.linalg import LinAlgError +from scipy.signal import lti, cont2discrete +# from exceptions import Exception +import warnings +from .lti import LTI, timebase, timebaseEqual, isdtime +from .xferfcn import _convertToTransferFunction +from copy import deepcopy + +__all__ = ['StateSpace', 'ss', 'rss', 'drss', 'tf2ss', 'ssdata'] + +class StateSpace(LTI): + """StateSpace(A, B, C, D[, dt]) + + A class for representing state-space models + + The StateSpace class is used to represent state-space realizations of linear + time-invariant (LTI) systems: + + dx/dt = A x + B u + y = C x + D u + + where u is the input, y is the output, and x is the state. + + The main data members are the A, B, C, and D matrices. The class also + keeps track of the number of states (i.e., the size of A). + + Discrete-time state space system are implemented by using the 'dt' instance + variable and setting it to the sampling period. If 'dt' is not None, + then it must match whenever two state space systems are combined. + Setting dt = 0 specifies a continuous system, while leaving dt = None + means the system timebase is not specified. If 'dt' is set to True, the + system will be treated as a discrete time system with unspecified + sampling time. + """ + + def __init__(self, *args): + """ + StateSpace(A, B, C, D[, dt]) + + Construct a state space object. + + The default constructor is StateSpace(A, B, C, D), where A, B, C, D + are matrices or equivalent objects. To create a discrete time system, + use StateSpace(A, B, C, D, dt) where 'dt' is the sampling time (or + True for unspecified sampling time). To call the copy constructor, + call StateSpace(sys), where sys is a StateSpace object. + + """ + + if len(args) == 4: + # The user provided A, B, C, and D matrices. + (A, B, C, D) = args + dt = None; + elif len(args) == 5: + # Discrete time system + (A, B, C, D, dt) = args + elif len(args) == 1: + # Use the copy constructor. + if not isinstance(args[0], StateSpace): + raise TypeError("The one-argument constructor can only take in a StateSpace object. Received %s." % type(args[0])) + A = args[0].A + B = args[0].B + C = args[0].C + D = args[0].D + try: + dt = args[0].dt + except NameError: + dt = None; + else: + raise ValueError("Needs 1 or 4 arguments; received %i." % len(args)) + + A, B, C, D = map(matrix, [A, B, C, D]) + + # TODO: use super here? + LTI.__init__(self, inputs=D.shape[1], outputs=D.shape[0], dt=dt) + self.A = A + self.B = B + self.C = C + self.D = D + + self.states = A.shape[1] + + 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) + + # Check that the matrix sizes are consistent. + if self.states != A.shape[0]: + raise ValueError("A must be square.") + if self.states != B.shape[0]: + raise ValueError("A and B must have the same number of rows.") + if self.states != C.shape[1]: + raise ValueError("A and C must have the same number of columns.") + if self.inputs != B.shape[1]: + raise ValueError("B and D must have the same number of columns.") + if self.outputs != C.shape[0]: + raise ValueError("C and D must have the same number of rows.") + + # Check for states that don't do anything, and remove them. + self._remove_useless_states() + + def _remove_useless_states(self): + """Check for states that don't do anything, and remove them. + + Scan the A, B, and C matrices for rows or columns of zeros. If the + zeros are such that a particular state has no effect on the input-output + dynamics, then remove that state from the A, B, and C matrices. + + """ + + # Indices of useless states. + useless = [] + + # 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)))): + useless.append(i) + # To avoid duplucate 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)))): + useless.append(i) + + # Remove the useless states. + self.A = delete(self.A, useless, 0) + self.A = delete(self.A, useless, 1) + self.B = delete(self.B, useless, 0) + self.C = delete(self.C, useless, 1) + + self.states = self.A.shape[0] + self.inputs = self.B.shape[1] + self.outputs = self.C.shape[0] + + 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" + #! TODO: replace with standard calls to lti functions + if (type(self.dt) == bool and self.dt == 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__ + + # Negation of a system + def __neg__(self): + """Negate a state space system.""" + + return StateSpace(self.A, self.B, -self.C, -self.D, self.dt) + + # Addition of two state space systems (parallel interconnection) + def __add__(self, other): + """Add two LTI systems (parallel connection).""" + + # Check for a couple of special cases + if (isinstance(other, (int, float, complex))): + # Just adding a scalar; put it in the D matrix + A, B, C = self.A, self.B, self.C; + D = self.D + other; + dt = self.dt + else: + other = _convertToStateSpace(other) + + # Check to make sure the dimensions are OK + 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 == None and other.dt != None): + dt = other.dt # use dt from second argument + elif (other.dt == None and self.dt != 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) + B = concatenate((self.B, other.B), axis=0) + C = concatenate((self.C, other.C), axis=1) + D = self.D + other.D + + return StateSpace(A, B, C, D, dt) + + # Right addition - just switch the arguments + def __radd__(self, other): + """Right add two LTI systems (parallel connection).""" + + return self + other + + # Subtraction of two state space systems (parallel interconnection) + def __sub__(self, other): + """Subtract two LTI systems.""" + + return self + (-other) + + def __rsub__(self, other): + """Right subtract two LTI systems.""" + + return other + (-self) + + # Multiplication of two state space systems (series interconnection) + def __mul__(self, other): + """Multiply two LTI objects (serial connection).""" + + # Check for a couple of special cases + if isinstance(other, (int, float, complex)): + # Just multiplying by a scalar; change the output + A, B = self.A, self.B + C = self.C * other + D = self.D * other + dt = self.dt + else: + other = _convertToStateSpace(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)) + + # Figure out the sampling time to use + if (self.dt == None and other.dt != None): + dt = other.dt # use dt from second argument + elif (other.dt == None and self.dt != 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) + B = concatenate((other.B, self.B * other.D), axis=0) + C = concatenate((self.D * other.C, self.C),axis=1) + D = self.D * other.D + + return StateSpace(A, B, C, D, dt) + + # Right multiplication of two state space systems (series interconnection) + # Just need to convert LH argument to a state space object + # TODO: __rmul__ only works for special cases (??) + def __rmul__(self, other): + """Right multiply two LTI objects (serial connection).""" + + # Check for a couple of special cases + if isinstance(other, (int, float, complex)): + # Just multiplying by a scalar; change the input + A, C = self.A, self.C; + B = self.B * other; + D = self.D * other; + return StateSpace(A, B, C, D, self.dt) + + # is lti, and convertible? + if isinstance(other, LTI): + return _convertToStateSpace(other) * self + + # try to treat this as a matrix + try: + X = matrix(other) + C = X * self.C + D = X * self.D + return StateSpace(self.A, self.B, C, D, self.dt) + + except Exception as e: + print(e) + pass + raise TypeError("can't interconnect systems") + + # TODO: __div__ and __rdiv__ are not written yet. + def __div__(self, other): + """Divide two LTI systems.""" + + raise NotImplementedError("StateSpace.__div__ is not implemented yet.") + + def __rdiv__(self, other): + """Right divide two LTI systems.""" + + raise NotImplementedError("StateSpace.__rdiv__ is not implemented yet.") + + # TODO: add discrete time check + def evalfr(self, omega): + """Evaluate a SS system's transfer function at a single frequency. + + self.evalfr(omega) returns the value of the transfer function matrix with + input value s = i * omega. + + """ + # Figure out the point to evaluate the transfer function + if isdtime(self, strict=True): + dt = timebase(self) + s = exp(1.j * omega * dt) + if (omega * dt > pi): + warnings.warn("evalfr: frequency evaluation above Nyquist frequency") + else: + s = omega * 1.j + + return self.horner(s) + + def horner(self, s): + '''Evaluate the systems's transfer function for a complex variable + + Returns a matrix of values evaluated at complex variable s. + ''' + 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 + def freqresp(self, omega): + """Evaluate the system's transfer func. at a list of ang. frequencies. + + mag, phase, omega = self.freqresp(omega) + + reports the value of the magnitude, phase, and angular frequency of the + system's transfer function matrix evaluated at s = i * omega, where + omega is a list of angular frequencies, and is a sorted version of the + input omega. + + """ + # when evaluating at many frequencies, much faster to convert to + # transfer function first and then evaluate, than to solve an + # n-dimensional linear system at each frequency + tf = _convertToTransferFunction(self) + return tf.freqresp(omega) + + # Compute poles and zeros + def pole(self): + """Compute the poles of a state space system.""" + + return eigvals(self.A) if self.states else np.array([]) + + def zero(self): + """Compute the zeros of a state space system.""" + + if self.inputs > 1 or self.outputs > 1: + raise NotImplementedError("StateSpace.zeros is currently \ +implemented only for SISO systems.") + + den = poly1d(poly(self.A)) + # Compute the numerator based on zeros + #! TODO: This is currently limited to SISO systems + num = poly1d(poly(self.A - dot(self.B, self.C)) + ((self.D[0, 0] - 1) * + den)) + + return roots(num) + + # Feedback around a state space system + def feedback(self, other=1, sign=-1): + """Feedback interconnection between two LTI systems.""" + + other = _convertToStateSpace(other) + + # 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.") + + # Figure out the sampling time to use + if (self.dt == None and other.dt != None): + dt = other.dt # use dt from second argument + elif (other.dt == None and self.dt != None) or \ + timebaseEqual(self, other): + dt = self.dt # use dt from first argument + else: + raise ValueError("Systems have different sampling times") + + A1 = self.A + B1 = self.B + C1 = self.C + D1 = self.D + A2 = other.A + B2 = other.B + C2 = other.C + D2 = other.D + + F = eye(self.inputs) - sign * D2 * D1 + if matrix_rank(F) != self.inputs: + raise ValueError("I - sign * D2 * D1 is singular to working precision.") + + # Precompute F\D2 and F\C2 (E = inv(F)) + # We can solve two linear systems in one pass, since the + # coefficients matrix F is the same. Thus, we perform the LU + # decomposition (cubic runtime complexity) of F only once! + # The remaining back substitutions are only quadratic in runtime. + E_D2_C2 = solve(F, concatenate((D2, C2), axis=1)) + E_D2 = E_D2_C2[:, :other.inputs] + E_C2 = E_D2_C2[:, other.inputs:] + + T1 = eye(self.outputs) + sign * D1 * E_D2 + T2 = eye(self.inputs) + sign * E_D2 * D1 + + A = concatenate( + (concatenate( + (A1 + sign * B1 * E_D2 * C1, sign * B1 * E_C2), axis=1), + concatenate( + (B2 * T1 * C1, A2 + sign * B2 * D1 * E_C2), axis=1)), + axis=0) + B = concatenate((B1 * T2, B2 * D1 * T2), axis=0) + C = concatenate((T1 * C1, sign * D1 * E_C2), axis=1) + D = D1 * T2 + + return StateSpace(A, B, C, D, dt) + + def minreal(self, tol=0.0): + """Calculate a minimal realization, removes unobservable and + uncontrollable states""" + if self.states: + try: + from slycot import tb01pd + B = empty((self.states, max(self.inputs, self.outputs))) + B[:,:self.inputs] = self.B + C = empty((max(self.outputs, self.inputs), self.states)) + 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) + 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. + + For instance, + + >>> out = ssobject.returnScipySignalLTI() + >>> out[3][5] + + is a signal.scipy.lti object corresponding to the transfer function from + the 6th input to the 4th output.""" + + # Preallocate the output. + out = [[[] for j in range(self.inputs)] for i in range(self.outputs)] + + for i in range(self.outputs): + for j in range(self.inputs): + out[i][j] = lti(asarray(self.A), asarray(self.B[:, j]), + asarray(self.C[i, :]), asarray(self.D[i, j])) + + return out + + def append(self, other): + """Append a second model to the present model. The second + model is converted to state-space if necessary, inputs and + outputs are appended and their order is preserved""" + if not isinstance(other, StateSpace): + other = _convertToStateSpace(other) + + if self.dt != other.dt: + raise ValueError("Systems must have the same time step") + + n = self.states + other.states + m = self.inputs + other.inputs + p = self.outputs + other.outputs + A = zeros( (n, n) ) + B = zeros( (n, m) ) + C = zeros( (p, n) ) + D = zeros( (p, m) ) + A[:self.states,:self.states] = self.A + A[self.states:,self.states:] = other.A + B[:self.states,:self.inputs] = self.B + B[self.states:,self.inputs:] = other.B + C[:self.outputs,:self.states] = self.C + C[self.outputs:,self.states:] = other.C + D[:self.outputs,:self.inputs] = self.D + D[self.outputs:,self.inputs:] = other.D + return StateSpace(A, B, C, D, self.dt) + + def __getitem__(self, indices): + """Array style acces""" + if len(indices) != 2: + raise IOError('must provide indices of length 2 for state space') + i = indices[0] + j = indices[1] + return StateSpace(self.A, + self.B[:,j], + self.C[i,:], + self.D[i,j], self.dt) + + def sample(self, Ts, method='zoh', alpha=None): + """Convert a continuous time system to discrete time + + Creates a discrete-time system from a continuous-time system by + sampling. Multiple methods of conversion are supported. + + Parameters + ---------- + Ts : float + Sampling period + method : {"gbt", "bilinear", "euler", "backward_diff", "zoh"} + Which method to use: + + * gbt: generalized bilinear transformation + * bilinear: Tustin's approximation ("gbt" with alpha=0.5) + * euler: Euler (or forward differencing) method ("gbt" with alpha=0) + * backward_diff: Backwards differencing ("gbt" with alpha=1.0) + * zoh: zero-order hold (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 + ------- + sysd : StateSpace system + Discrete time system, with sampling rate Ts + + Notes + ----- + Uses the command 'cont2discrete' from scipy.signal + + Examples + -------- + >>> sys = StateSpace(0, 1, 1, 0) + >>> sysd = sys.sample(0.5, method='bilinear') + + """ + if not self.isctime(): + raise ValueError("System must be continuous time system") + + sys = (self.A, self.B, self.C, self.D) + Ad, Bd, C, D, dt = cont2discrete(sys, Ts, method, alpha) + return StateSpace(Ad, Bd, C, D, dt) + + def dcgain(self): + """Return the zero-frequency gain + + The zero-frequency gain of a continuous-time state-space + system is given by: + + .. math: G(0) = - C A^{-1} B + D + + and of a discrete-time state-space system by: + + .. math: G(1) = C (I - A)^{-1} B + D + + Returns + ------- + gain : ndarray + An array of shape (outputs,inputs); the array will either + be the zero-frequency (or DC) gain, or, if the frequency + response is singular, the array will be filled with np.nan. + """ + try: + if self.isctime(): + gain = np.asarray(self.D - + self.C.dot(np.linalg.solve(self.A, self.B))) + else: + gain = self.horner(1) + except LinAlgError: + # eigenvalue at DC + gain = np.tile(np.nan,(self.outputs,self.inputs)) + return np.squeeze(gain) + + +# TODO: add discrete time check +def _convertToStateSpace(sys, **kw): + """Convert a system to state space form (if needed). + + If sys is already a state space, then it is returned. If sys is a transfer + function object, then it is converted to a state space and returned. If sys + is a scalar, then the number of inputs and outputs can be specified + manually, as in: + + >>> sys = _convertToStateSpace(3.) # Assumes inputs = outputs = 1 + >>> sys = _convertToStateSpace(1., inputs=3, outputs=2) + + In the latter example, A = B = C = 0 and D = [[1., 1., 1.] + [1., 1., 1.]]. + + """ + + from .xferfcn import TransferFunction + import itertools + if isinstance(sys, StateSpace): + if len(kw): + raise TypeError("If sys is a StateSpace, _convertToStateSpace \ +cannot take keywords.") + + # Already a state space system; just return it + return sys + elif isinstance(sys, TransferFunction): + try: + from slycot import td04ad + if len(kw): + 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. + num, den = sys._common_den() + # Make a list of the orders of the denominator polynomials. + index = [len(den) - 1 for i in range(sys.outputs)] + # Repeat the common denominator along the rows. + den = array([den for i in range(sys.outputs)]) + #! TODO: transfer function to state space conversion is still buggy! + #print num + #print shape(num) + ssout = td04ad('R',sys.inputs, sys.outputs, index, den, num,tol=0.0) + + states = ssout[0] + return StateSpace(ssout[1][:states, :states], + ssout[2][:states, :sys.inputs], + ssout[3][:sys.outputs, :states], + ssout[4], sys.dt) + except ImportError: + # No Slycot. Scipy tf->ss can't handle MIMO, but static + # MIMO is an easy special case we can check for here + maxn = max(max(len(n) for n in nrow) + for nrow in sys.num) + maxd = max(max(len(d) for d in drow) + for drow in sys.den) + if 1==maxn and 1==maxd: + D = empty((sys.outputs,sys.inputs),dtype=float) + for i,j in itertools.product(range(sys.outputs),range(sys.inputs)): + D[i,j] = sys.num[i][j][0] / sys.den[i][j][0] + return StateSpace([], [], [], D, sys.dt) + else: + if (sys.inputs != 1 or sys.outputs != 1): + raise TypeError("No support for MIMO without slycot") + + # TODO: do we want to squeeze first and check dimenations? + # I think this will fail if num and den aren't 1-D after + # the squeeze + lti_sys = lti(squeeze(sys.num), squeeze(sys.den)) + return StateSpace(lti_sys.A, lti_sys.B, lti_sys.C, lti_sys.D, + sys.dt) + + elif isinstance(sys, (int, float, complex)): + if "inputs" in kw: + inputs = kw["inputs"] + else: + inputs = 1 + if "outputs" in kw: + outputs = kw["outputs"] + else: + outputs = 1 + + # Generate a simple state space system of the desired dimension + # 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))) + + # If this is a matrix, try to create a constant feedthrough + try: + D = matrix(sys) + outputs, inputs = D.shape + + return StateSpace(0., zeros((1, inputs)), zeros((outputs, 1)), D) + except Exception(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. + + This does the actual random state space generation expected from rss and + drss. type is 'c' for continuous systems and 'd' for discrete systems. + + """ + + # Probability of repeating a previous root. + pRepeat = 0.05 + # Probability of choosing a real root. Note that when choosing a complex + # root, the conjugate gets chosen as well. So the expected proportion of + # real roots is pReal / (pReal + 2 * (1 - pReal)). + pReal = 0.6 + # Probability that an element in B or C will not be masked out. + pBCmask = 0.8 + # Probability that an element in D will not be masked out. + pDmask = 0.3 + # Probability that D = 0. + pDzero = 0.5 + + # Check for valid input arguments. + if states < 1 or states % 1: + 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) + if outputs < 1 or outputs % 1: + 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 + i = 0 + + while i < states: + 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: + # Copy previous real pole. + poles[i] = poles[i-1] + i += 1 + else: + # Copy previous complex conjugate pair of poles. + poles[i:i+2] = poles[i-2:i] + i += 2 + elif rand() < pReal or i == states - 1: + # No-oscillation pole. + if type == 'c': + poles[i] = -exp(randn()) + 0.j + elif type == 'd': + poles[i] = 2. * rand() - 1. + i += 1 + else: + # Complex conjugate pair of oscillating poles. + if type == 'c': + poles[i] = complex(-exp(randn()), 3. * exp(randn())) + elif type == 'd': + mag = rand() + phase = 2. * pi * rand() + poles[i] = complex(mag * cos(phase), + mag * sin(phase)) + poles[i+1] = complex(poles[i].real, -poles[i].imag) + i += 2 + + # Now put the poles in A as real blocks on the diagonal. + A = zeros((states, states)) + i = 0 + while i < states: + if poles[i].imag == 0: + 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 + i += 2 + # Finally, apply a transformation so that A is not block-diagonal. + while True: + T = randn(states, states) + try: + A = dot(solve(T, A), T) # A = T \ A * T + break + except LinAlgError: + # In the unlikely event that T is rank-deficient, iterate again. + pass + + # Make the remaining matrices. + B = randn(states, inputs) + C = randn(outputs, states) + D = randn(outputs, inputs) + + # Make masks to zero out some of the elements. + while True: + Bmask = rand(states, inputs) < pBCmask + if any(Bmask): # Retry if we get all zeros. + break + while True: + Cmask = rand(outputs, states) < pBCmask + if any(Cmask): # Retry if we get all zeros. + break + if rand() < pDzero: + Dmask = zeros((outputs, inputs)) + else: + Dmask = rand(outputs, inputs) < pDmask + + # Apply masks. + B = B * Bmask + C = C * Cmask + D = D * Dmask + + return StateSpace(A, B, C, D) + +# Convert a MIMO system to a SISO system +# TODO: add discrete time check +def _mimo2siso(sys, input, output, warn_conversion=False): + #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.) + + The input and output that are used in the SISO system can be selected + with the parameters ``input`` and ``output``. All other inputs are set + to 0, all other outputs are ignored. + + If ``sys`` is already a SISO system, it will be returned unaltered. + + Parameters + ---------- + sys: StateSpace + Linear (MIMO) system that should be converted. + input: int + Index of the input that will become the SISO system's only input. + output: int + Index of the output that will become the SISO system's only output. + warn_conversion: bool + If True: print a warning message when sys is a MIMO system. + 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.") + if not (0 <= input < sys.inputs): + raise ValueError("Selected input does not exist. " + "Selected input: {sel}, " + "number of system inputs: {ext}." + .format(sel=input, ext=sys.inputs)) + if not (0 <= output < sys.outputs): + raise ValueError("Selected output does not exist. " + "Selected output: {sel}, " + "number of system outputs: {ext}." + .format(sel=output, ext=sys.outputs)) + #Convert sys to SISO if necessary + if sys.inputs > 1 or sys.outputs > 1: + if warn_conversion: + warnings.warn("Converting MIMO system to SISO system. " + "Only input {i} and output {o} are used." + .format(i=input, o=output)) + # $X = A*X + B*U + # Y = C*X + D*U + new_B = sys.B[:, input] + new_C = sys.C[output, :] + new_D = sys.D[output, input] + sys = StateSpace(sys.A, new_B, new_C, new_D, sys.dt) + + return sys + +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 + multiple outputs.) + + The input that is used in the SIMO system can be selected with the + parameter ``input``. All other inputs are set to 0, all other + outputs are ignored. + + If ``sys`` is already a SIMO system, it will be returned unaltered. + + Parameters + ---------- + sys: StateSpace + Linear (MIMO) system that should be converted. + input: int + Index of the input that will become the SIMO system's only input. + warn_conversion: bool + If True: print a warning message when sys is a MIMO system. + Warn that a conversion will take place. + + Returns: + -------- + sys: StateSpace + The converted (SIMO) system. + """ + if not (isinstance(input, int)): + raise TypeError("Parameter ``input`` be an integer number.") + if not (0 <= input < sys.inputs): + raise ValueError("Selected input does not exist. " + "Selected input: {sel}, " + "number of system inputs: {ext}." + .format(sel=input, ext=sys.inputs)) + #Convert sys to SISO if necessary + if sys.inputs > 1: + if warn_conversion: + warnings.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] + new_D = sys.D[:, input] + sys = StateSpace(sys.A, new_B, sys.C, new_D, sys.dt) + + return sys + +def ss(*args): + """ + Create a state space system. + + The function accepts either 1, 4 or 5 parameters: + + ``ss(sys)`` + Convert a linear system into space system form. Always creates a + new system, even if sys is already a StateSpace object. + + ``ss(A, B, C, D)`` + Create a state space system from the matrices of its state and + output equations: + + .. math:: + \dot x = A \cdot x + B \cdot u + + y = C \cdot x + D \cdot u + + ``ss(A, B, C, D, dt)`` + Create a discrete-time state space system from the matrices of + its state and output equations: + + .. math:: + x[k+1] = A \cdot x[k] + B \cdot u[k] + + y[k] = C \cdot x[k] + D \cdot u[ki] + + The matrices can be given as *array like* data types or strings. + Everything that the constructor of :class:`numpy.matrix` accepts is + permissible here too. + + Parameters + ---------- + sys: StateSpace or TransferFunction + A linear system + A: array_like or string + System matrix + B: array_like or string + Control matrix + C: array_like or string + Output matrix + D: array_like or string + Feed forward matrix + dt: If present, specifies the sampling period and a discrete time + system is created + + Returns + ------- + out: :class:`StateSpace` + The new linear system + + Raises + ------ + ValueError + if matrix sizes are not self-consistent + + See Also + -------- + tf + ss2tf + tf2ss + + Examples + -------- + >>> # Create a StateSpace object from four "matrices". + >>> sys1 = ss("1. -2; 3. -4", "5.; 7", "6. 8", "9.") + + >>> # Convert a TransferFunction to a StateSpace object. + >>> sys_tf = tf([2.], [1., 3]) + >>> sys2 = ss(sys_tf) + + """ + + if len(args) == 4 or len(args) == 5: + return StateSpace(*args) + elif len(args) == 1: + from .xferfcn import TransferFunction + sys = args[0] + if isinstance(sys, StateSpace): + return deepcopy(sys) + 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)) + else: + raise ValueError("Needs 1 or 4 arguments; received %i." % len(args)) + +def tf2ss(*args): + """ + Transform a transfer function to a state space system. + + The function accepts either 1 or 2 parameters: + + ``tf2ss(sys)`` + Convert a linear system into transfer function form. Always creates + a new system, even if sys is already a TransferFunction object. + + ``tf2ss(num, den)`` + Create a transfer function system from its numerator and denominator + polynomial coefficients. + + For details see: :func:`tf` + + Parameters + ---------- + sys: LTI (StateSpace or TransferFunction) + A linear system + num: array_like, or list of list of array_like + Polynomial coefficients of the numerator + den: array_like, or list of list of array_like + Polynomial coefficients of the denominator + + Returns + ------- + out: StateSpace + New linear system in state space form + + Raises + ------ + ValueError + if `num` and `den` have invalid or unequal dimensions, or if an + invalid number of arguments is passed in + TypeError + if `num` or `den` are of incorrect type, or if sys is not a + TransferFunction object + + See Also + -------- + ss + tf + ss2tf + + Examples + -------- + >>> num = [[[1., 2.], [3., 4.]], [[5., 6.], [7., 8.]]] + >>> den = [[[9., 8., 7.], [6., 5., 4.]], [[3., 2., 1.], [-1., -2., -3.]]] + >>> sys1 = tf2ss(num, den) + + >>> sys_tf = tf(num, den) + >>> sys2 = tf2ss(sys_tf) + + """ + + from .xferfcn import TransferFunction + if len(args) == 2 or len(args) == 3: + # Assume we were given the num, den + return _convertToStateSpace(TransferFunction(*args)) + + elif len(args) == 1: + sys = args[0] + if not isinstance(sys, TransferFunction): + 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)) + +def rss(states=1, outputs=1, inputs=1): + """ + Create a stable **continuous** random state space object. + + Parameters + ---------- + states: integer + Number of state variables + inputs: integer + Number of system inputs + outputs: integer + Number of system outputs + + Returns + ------- + sys: StateSpace + The randomly created linear system + + Raises + ------ + ValueError + if any input is not a positive integer + + See Also + -------- + drss + + Notes + ----- + If the number of states, inputs, or outputs is not specified, then the + missing numbers are assumed to be 1. The poles of the returned system + will always have a negative real part. + + """ + + return _rss_generate(states, inputs, outputs, 'c') + +def drss(states=1, outputs=1, inputs=1): + """ + Create a stable **discrete** random state space object. + + Parameters + ---------- + states: integer + Number of state variables + inputs: integer + Number of system inputs + outputs: integer + Number of system outputs + + Returns + ------- + sys: StateSpace + The randomly created linear system + + Raises + ------ + ValueError + if any input is not a positive integer + + See Also + -------- + rss + + Notes + ----- + If the number of states, inputs, or outputs is not specified, then the + missing numbers are assumed to be 1. The poles of the returned system + will always have a magnitude less than 1. + + """ + + return _rss_generate(states, inputs, outputs, 'd') + +def ssdata(sys): + ''' + Return state space data objects for a system + + Parameters + ---------- + sys: LTI (StateSpace, or TransferFunction) + LTI system whose data will be returned + + Returns + ------- + (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) diff --git a/control/xferfcn.py b/control/xferfcn.py index 982b8c5ba..82cbf4a9b 100644 --- a/control/xferfcn.py +++ b/control/xferfcn.py @@ -66,7 +66,9 @@ class TransferFunction(LTI): - """A class for representing transfer functions + """TransferFunction(num, den[, dt]) + + A class for representing transfer functions The TransferFunction class is used to represent systems in transfer function form. @@ -84,17 +86,20 @@ class TransferFunction(LTI): non-zero value, then it must match whenever two transfer functions are combined. If 'dt' is set to True, the system will be treated as a discrete time system with unspecified sampling time. - """ def __init__(self, *args): - """Construct a transfer function. + """TransferFunction(num, den[, dt]) + + Construct a transfer function. The default constructor is TransferFunction(num, den), where num and den are lists of lists of arrays containing polynomial coefficients. To create a discrete time transfer funtion, use TransferFunction(num, - den, dt). To call the copy constructor, call TransferFunction(sys), - where sys is a TransferFunction object (continuous or discrete). + den, dt) where 'dt' is the sampling time (or True for unspecified + sampling time). To call the copy constructor, call + TransferFunction(sys), where sys is a TransferFunction object + (continuous or discrete). """ args = deepcopy(args) @@ -1134,10 +1139,11 @@ def _convertToTransferFunction(sys, **kw): def tf(*args): - """ + """tf(num, den[, dt]) + Create a transfer function system. Can create MIMO systems. - The function accepts either 1 or 2 parameters: + The function accepts either 1, 2, or 3 parameters: ``tf(sys)`` Convert a linear system into transfer function form. Always creates @@ -1182,6 +1188,7 @@ def tf(*args): See Also -------- + TransferFunction ss ss2tf tf2ss @@ -1226,7 +1233,8 @@ def tf(*args): raise ValueError("Needs 1 or 2 arguments; received %i." % len(args)) def ss2tf(*args): - """ + """ss2tf(sys) + Transform a state space system to a transfer function. The function accepts either 1 or 4 parameters: diff --git a/doc/control.rst b/doc/control.rst index f3d124ee5..294701f8f 100644 --- a/doc/control.rst +++ b/doc/control.rst @@ -24,14 +24,14 @@ System creation System interconnections ======================= .. autosummary:: - :toctree: generated/ + :toctree: generated/ - append - connect - feedback - negate - parallel - series + append + connect + feedback + negate + parallel + series Frequency domain plotting ========================= @@ -118,6 +118,8 @@ Model simplification tools era markov +.. _utility-and-conversions: + Utility functions and conversions ================================= .. autosummary:: diff --git a/doc/conventions.rst b/doc/conventions.rst index 881a0a6bb..9c3b46ce5 100644 --- a/doc/conventions.rst +++ b/doc/conventions.rst @@ -9,6 +9,98 @@ Library conventions The python-control library uses a set of standard conventions for the way that different types of standard information used by the library. +LTI system representation +========================= + +Linear time invariant (LTI) systems are represented in python-control in +state space, transfer function, or frequency response data (FRD) form. Most +functions in the toolbox will operate on any of these data types and +functions for converting between between compatible types is provided. + +State space systems +------------------- +The :class:`StateSpace` class is used to represent state-space realizations +of linear time-invariant (LTI) systems: + +.. math:: + + \frac{dx}{dt} &= A x + B u \\ + y &= C x + D u + +where u is the input, y is the output, and x is the state. + +To create a state space system, use the :class:`StateSpace` constructor: + + sys = StateSpace(A, B, C, D) + +State space systems can be manipulated using standard arithmetic operations +as well as the :func:`feedback`, :func:`parallel`, and :func:`series` +function. A full list of functions can be found in :ref:`function-ref`. + +Transfer functions +------------------ +The :class:`TransferFunction` class is used to represent input/output +transfer functions + +.. math:: + + G(s) = \frac{\text{num}(s)}{\text{den}(s)} + = \frac{a_0 s^n + a_1 s^{n-1} + \cdots + a_n} + {b_0 s^m + b_1 s^{m-1} + \cdots + b_m}, + +where n is generally greater than or equal to m (for a proper transfer +function). + +To create a transfer function, use the :class:`TransferFunction` +constructor: + + sys = TransferFunction(num, den) + +Transfer functions can be manipulated using standard arithmetic operations +as well as the :func:`feedback`, :func:`parallel`, and :func:`series` +function. A full list of functions can be found in :ref:`function-ref`. + +FRD (frequency response data) systems +------------------------------------- +The :class:`FRD` class is used to represent systems in frequency response +data form. + +The main data members are `omega` and `fresp`, where `omega` is a 1D array +with the frequency points of the response, and `fresp` is a 3D array, with +the first dimension corresponding to the output index of the FRD, the second +dimension corresponding to the input index, and the 3rd dimension +corresponding to the frequency points in omega. + +FRD systems have a somewhat more limited set of functions that are +available, although all of the standard algebraic manipulations can be +performed. + +Discrete time systems +--------------------- +By default, all systems are considered to be continuous time systems. A +discrete time system is created by specifying the 'time base' dt. The time +base argument can be given when a system is constructed: + +* dt = None: continuous time +* dt = number: discrete time system with sampling period 'dt' +* dt = True: discrete time with unspecified sampling period + +Only the :class:`StateSpace` and :class:`TransferFunction` classes allow +explicit representation of discrete time systems. + +Systems must have the same time base in order to be combined. For +continuous time systems, the :func:`sample_system` function or the +:meth:`StateSpace.sample` and :meth:`TransferFunction.sample` methods can be +used to create a discrete time system from a continuous time system. See +:ref:`utility-and-conversions`. + +Conversion between representations +---------------------------------- +LTI systems can be converted between representations either by calling the +constructor for the desired data type using the original system as the sole +argument or using the explicit conversion functions :func:`ss2tf` and +:func:`tf2ss`. + .. _time-series-convention: Time series data