"""lti.py
The lti module contains the LTI parent class to the child classes StateSpace
and TransferFunction. It is designed for use in the python-control library.
Routines in this module:
LTI.__init__
isdtime()
isctime()
timebase()
timebaseEqual()
"""
import numpy as np
from numpy import absolute, real
__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.
LTI is the parent to the StateSpace and TransferFunction child
classes. It contains the number of inputs and outputs, and the
timebase (dt) for the system.
The timebase for the system, dt, is used to specify whether the
system is operating in continuous or discrete time. It can have
the following values:
* dt = None No timebase specified
* dt = 0 Continuous time system
* dt > 0 Discrete time system with sampling time dt
* dt = True Discrete time system with unspecified sampling time
When two LTI systems are combined, their timebases much match. A system
with timebase None can be combined with a system having a specified
timebase, and the result will have the timebase of the latter system.
"""
def __init__(self, inputs=1, outputs=1, dt=None):
"""Assign the LTI object's numbers of inputs and ouputs."""
# Data members common to StateSpace and TransferFunction.
self.inputs = inputs
self.outputs = outputs
self.dt = dt
def isdtime(self, strict=False):
"""
Check to see if a system is a discrete-time system
Parameters
----------
strict: bool, optional
If strict is True, make sure that timebase is not None. Default
is False.
"""
# If no timebase is given, answer depends on strict flag
if self.dt == None:
return True if not strict else False
# Look for dt > 0 (also works if dt = True)
return self.dt > 0
def isctime(self, strict=False):
"""
Check to see if a system is a continuous-time system
Parameters
----------
sys : LTI system
System to be checked
strict: bool, optional
If strict is True, make sure that timebase is not None. Default
is False.
"""
# If no timebase is given, answer depends on strict flag
if self.dt is None:
return True if not strict else False
return self.dt == 0
def issiso(self):
'''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
Returns
-------
wn : array
Natural frequencies for each system pole
zeta : array
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
else:
splane_poles = poles
wn = absolute(splane_poles)
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__))
# Test to see if a system is SISO
def issiso(sys, strict=False):
"""
Check to see if a system is single input, single output
Parameters
----------
sys : LTI system
System to be checked
strict: bool (default = False)
If strict is True, do not treat scalars as SISO
"""
if isinstance(sys, (int, float, complex, np.number)) and not strict:
return True
elif not isinstance(sys, LTI):
raise ValueError("Object is not an LTI system")
# 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
dt = timebase(sys)
returns the timebase for a system 'sys'. If the strict option is
set to False, dt = True will be returned as 1.
"""
# System needs to be either a constant or an LTI system
if isinstance(sys, (int, float, complex, np.number)):
return None
elif not isinstance(sys, LTI):
raise ValueError("Timebase not defined")
# Return the sample time, with converstion to float if strict is false
if (sys.dt == None):
return None
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
timebaseEqual(sys1, sys2)
returns True if the timebases for the two systems are compatible. By
default, systems with timebase 'None' are compatible with either
discrete or continuous timebase systems. If two systems have a discrete
timebase (dt > 0) then their timebases must be equal.
"""
if (type(sys1.dt) == bool or type(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):
# One or the other is unspecified => the other can be anything
return True
else:
return sys1.dt == sys2.dt
# Find a common timebase between two or more systems
def _find_timebase(sys1, *sysn):
"""Find the common timebase between systems, otherwise return False"""
# Create a list of systems to check
syslist = [sys1]
syslist.append(*sysn)
# Look for a common timebase
dt = None
for sys in syslist:
# Make sure time bases are consistent
if (dt is None and sys.dt is not None) or \
(dt is True and isdiscrete(sys)):
# Timebase was not specified; set to match this system
dt = sys.dt
elif dt != sys.dt:
return False
return dt
# Check to see if a system is a discrete time system
def isdtime(sys, strict=False):
"""
Check to see if a system is a discrete time system
Parameters
----------
sys : LTI system
System to be checked
strict: bool (default = False)
If strict is True, make sure that timebase is not None
"""
# Check to see if this is a constant
if isinstance(sys, (int, float, complex, np.number)):
# OK as long as strict checking is off
return True if not strict else False
# Check for a transfer function or state-space object
if isinstance(sys, LTI):
return sys.isdtime(strict)
# Check to see if object has a dt object
if hasattr(sys, 'dt'):
# If no timebase is given, answer depends on strict flag
if sys.dt == None:
return True if not strict else False
# Look for dt > 0 (also works if dt = True)
return sys.dt > 0
# 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):
"""
Check to see if a system is a continuous-time system
Parameters
----------
sys : LTI system
System to be checked
strict: bool (default = False)
If strict is True, make sure that timebase is not None
"""
# Check to see if this is a constant
if isinstance(sys, (int, float, complex, np.number)):
# OK as long as strict checking is off
return True if not strict else False
# Check for a transfer function or state space object
if isinstance(sys, LTI):
return sys.isctime(strict)
# Check to see if object has a dt object
if hasattr(sys, 'dt'):
# If no timebase is given, answer depends on strict flag
if sys.dt is None:
return True if not strict else False
return sys.dt == 0
# Got passed something we don't recognize
return False
def pole(sys):
"""
Compute system poles.
Parameters
----------
sys: StateSpace or TransferFunction
Linear system
Returns
-------
poles: ndarray
Array that contains the system's poles.
Raises
------
NotImplementedError
when called on a TransferFunction object
See Also
--------
zero
TransferFunction.pole
StateSpace.pole
"""
return sys.pole()
def zero(sys):
"""
Compute system zeros.
Parameters
----------
sys: StateSpace or TransferFunction
Linear system
Returns
-------
zeros: ndarray
Array that contains the system's zeros.
Raises
------
NotImplementedError
when called on a MIMO system
See Also
--------
pole
StateSpace.zero
TransferFunction.zero
"""
return sys.zero()
def damp(sys, doprint=True):
"""
Compute natural frequency, damping ratio, and poles of a system
The function takes 1 or 2 parameters
Parameters
----------
sys: LTI (StateSpace or TransferFunction)
A linear system object
doprint:
if true, print table with values
Returns
-------
wn: array
Natural frequencies of the poles
damping: array
Damping values
poles: array
Pole locations
Algorithm
---------
If the system is continuous,
wn = abs(poles)
Z = -real(poles)/poles.
If the system is discrete, the discrete poles are mapped to their
equivalent location in the s-plane via
s = log10(poles)/dt
and
wn = abs(s)
Z = -real(s)/wn.
See Also
--------
pole
"""
wn, damping, poles = sys.damp()
if doprint:
print('_____Eigenvalue______ Damping___ Frequency_')
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))
else:
print("%10.4g%+10.4gj %10.4g %10.4g" %
(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
number x.
To evaluate at a frequency, enter x = omega*j, where omega is the
frequency in radians
Parameters
----------
sys: StateSpace or TransferFunction
Linear system
x: scalar
Complex number
Returns
-------
fresp: ndarray
See Also
--------
freqresp
bode
Notes
-----
This function is a wrapper for StateSpace.evalfr and
TransferFunction.evalfr.
Examples
--------
>>> sys = ss("1. -2; 3. -4", "5.; 7", "6. 8", "9.")
>>> evalfr(sys, 1j)
array([[ 44.8-21.4j]])
>>> # This is the transfer function matrix evaluated at s = i.
.. todo:: Add example with MIMO system
"""
if issiso(sys):
return sys.horner(x)[0][0]
return sys.horner(x)
def freqresp(sys, omega):
"""
Frequency response of an LTI system at multiple angular frequencies.
Parameters
----------
sys: StateSpace or TransferFunction
Linear system
omega: array_like
List of frequencies
Returns
-------
mag: ndarray
phase: ndarray
omega: list, tuple, or ndarray
See Also
--------
evalfr
bode
Notes
-----
This function is a wrapper for StateSpace.freqresp and
TransferFunction.freqresp. The output omega is a sorted version of the
input omega.
Examples
--------
>>> sys = ss("1. -2; 3. -4", "5.; 7", "6. 8", "9.")
>>> mag, phase, omega = freqresp(sys, [0.1, 1., 10.])
>>> mag
array([[[ 58.8576682 , 49.64876635, 13.40825927]]])
>>> phase
array([[[-0.05408304, -0.44563154, -0.66837155]]])
.. todo::
Add example with MIMO system
#>>> sys = rss(3, 2, 2)
#>>> mag, phase, omega = freqresp(sys, [0.1, 1., 10.])
#>>> mag[0, 1, :]
#array([ 55.43747231, 42.47766549, 1.97225895])
#>>> phase[1, 0, :]
#array([-0.12611087, -1.14294316, 2.5764547 ])
#>>> # This is the magnitude of the frequency response from the 2nd
#>>> # input to the 1st output, and the phase (in radians) of the
#>>> # frequency response from the 1st input to the 2nd output, for
#>>> # s = 0.1i, i, 10i.
"""
return sys.freqresp(omega)
def dcgain(sys):
"""Return the zero-frequency (or DC) gain of the given system
Returns
-------
gain : ndarray
The zero-frequency gain, or np.nan if the system has a pole
at the origin
"""
return sys.dcgain()
