Menu
▾
▴
[python-control-discuss] SF.net python-control commit:[274] trunk
|
From: <mur...@us...> - 2013-06-09 14:58:04
|
Revision: 274
http://sourceforge.net/p/python-control/code/274
Author: murrayrm
Date: 2013-06-09 14:57:59 +0000 (Sun, 09 Jun 2013)
Log Message:
-----------
merging changes from Rane van Paassen - see 2013-05-23 ChangeLog entry for details
Modified Paths:
--------------
trunk/ChangeLog
trunk/src/bdalg.py
trunk/src/frdata.py
trunk/src/lti.py
trunk/src/margins.py
trunk/src/matlab.py
trunk/src/modelsimp.py
trunk/src/rlocus.py
trunk/src/statesp.py
trunk/src/timeresp.py
trunk/src/xferfcn.py
trunk/tests/convert_test.py
trunk/tests/matlab_test.py
trunk/tests/slycot_convert_test.py
trunk/tests/statesp_test.py
trunk/tests/timeresp_test.py
trunk/tests/xferfcn_test.py
Removed Paths:
-------------
branches/add-frd/
branches/mix-matrices-and-sys/
Property Changed:
----------------
trunk/
trunk/src/ctrlutil.py
Index: trunk
===================================================================
--- trunk 2013-05-31 11:53:29 UTC (rev 273)
+++ trunk 2013-06-09 14:57:59 UTC (rev 274)
Property changes on: trunk
___________________________________________________________________
Modified: svn:mergeinfo
## -1 +1,2 ##
/branches/add-frd:188-219
+/branches/mix-matrices-and-sys:249-273
\ No newline at end of property
Modified: trunk/ChangeLog
===================================================================
--- trunk/ChangeLog 2013-05-31 11:53:29 UTC (rev 273)
+++ trunk/ChangeLog 2013-06-09 14:57:59 UTC (rev 274)
@@ -1,3 +1,24 @@
+2013-05-23 Rene van Paassen <ren...@gm...>
+
+ * src/margin.py: re-vamped stability_margins function. Now
+ uses analytical method for TF and state-space, works on FRD
+ objects too, and can optionally return all margins found for
+ the analytical method.
+
+ * src/xferfcn.py and src/frdata.py: Now return array result
+ for fresp when called with an iterable frequency vector
+
+ * src/xferfcn.py and src/statesp.py: added minreal functions
+
+ * src/bdalg.py: implemented append and connect functions (old
+ version, with indexes, not with names)
+
+ * src/matlab.py: added connect, append, minreal
+
+ * src/xferfcn.py and src/statesp.py: added/corrected
+ __rmul__. In most cases a gain matrix (np.matrix) may now be
+ used in combination with ss or tf systems.
+
2012-11-10 Richard Murray <murray@altura.local>
* src/canonical.py: new module implementing conversions to selected
Modified: trunk/src/bdalg.py
===================================================================
--- trunk/src/bdalg.py 2013-05-31 11:53:29 UTC (rev 273)
+++ trunk/src/bdalg.py 2013-06-09 14:57:59 UTC (rev 274)
@@ -4,10 +4,12 @@
Routines in this module:
+append
series
parallel
negate
feedback
+connect
"""
@@ -236,3 +238,97 @@
sys2 = tf._convertToTransferFunction(sys2)
return sys1.feedback(sys2, sign)
+
+def append(*sys):
+ '''
+ Group models by appending their inputs and outputs
+
+ Forms an augmented system model, and appends the inputs and
+ outputs together. The system type will be the type of the first
+ system given; if you mix state-space systems and gain matrices,
+ make sure the gain matrices are not first.
+
+ Parameters.
+ -----------
+ sys1, sys2, ... sysn: StateSpace or Transferfunction
+ LTI systems to combine
+
+
+ Returns
+ -------
+ sys: LTI system
+ Combined LTI system, with input/output vectors consisting of all
+ input/output vectors appended
+
+ Examples
+ --------
+ >>> sys1 = ss("1. -2; 3. -4", "5.; 7", "6. 8", "9.")
+ >>> sys2 = ss("-1.", "1.", "1.", "0.")
+ >>> sys = append(sys1, sys2)
+
+ .. 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
+ multiple inputs and outputs. The inputs and outputs are connected
+ according to the interconnection matrix Q, and then the final
+ inputs and outputs are trimmed according to the inputs and outputs
+ listed in inputv and outputv.
+
+ Note: to have this work, inputs and outputs start counting at 1!!!!
+
+ Parameters.
+ -----------
+ sys: StateSpace Transferfunction
+ System to be connected
+ Q: 2d array
+ Interconnection matrix. First column gives the input to be connected
+ second column gives the output to be fed into this input. Negative
+ values for the second column mean the feedback is negative, 0 means
+ no connection is made
+ inputv: 1d array
+ list of final external inputs
+ outputv: 1d array
+ list of final external outputs
+
+ Returns
+ -------
+ sys: LTI system
+ Connected and trimmed LTI system
+
+ Examples
+ --------
+ >>> sys1 = ss("1. -2; 3. -4", "5.; 7", "6, 8", "9.")
+ >>> sys2 = ss("-1.", "1.", "1.", "0.")
+ >>> 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) )
+ for r in sp.array(Q).astype(int):
+ inp = r[0]-1
+ for outp in r[1:]:
+ if outp > 0 and outp <= sys.outputs:
+ K[inp,outp-1] = 1.
+ elif outp < 0 and -outp >= -sys.outputs:
+ 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)
Index: trunk/src/ctrlutil.py
===================================================================
--- trunk/src/ctrlutil.py 2013-05-31 11:53:29 UTC (rev 273)
+++ trunk/src/ctrlutil.py 2013-06-09 14:57:59 UTC (rev 274)
Property changes on: trunk/src/ctrlutil.py
___________________________________________________________________
Modified: svn:mergeinfo
## -1,2 +1,3 ##
/branches/add-frd/src/ctrlutil.py:188-219
+/branches/mix-matrices-and-sys/src/ctrlutil.py:249-273
/trunk/src/ctrlutil.py:140-142
\ No newline at end of property
Modified: trunk/src/frdata.py
===================================================================
--- trunk/src/frdata.py 2013-05-31 11:53:29 UTC (rev 273)
+++ trunk/src/frdata.py 2013-06-09 14:57:59 UTC (rev 274)
@@ -107,7 +107,7 @@
def __init__(self, *args, **kwargs):
"""Construct a transfer function.
- The default constructor is FRD(w, d), where w is an iterable of
+ The default constructor is FRD(d, w), where w is an iterable of
frequency points, and d is the matching frequency data.
If d is a single list, 1d array, or tuple, a SISO system description
@@ -142,8 +142,8 @@
# The user provided a response and a freq vector
self.fresp = array(args[0], dtype=complex)
if len(self.fresp.shape) == 1:
- self.fresp.reshape(1, 1, len(args[0]))
- self.omega = array(args[1])
+ self.fresp = self.fresp.reshape(1, 1, len(args[0]))
+ self.omega = array(args[1], dtype=float)
if len(self.fresp.shape) != 3 or \
self.fresp.shape[-1] != self.omega.shape[-1] or \
len(self.omega.shape) != 1:
@@ -189,9 +189,9 @@
if mimo:
outstr.append("Input %i to output %i:" % (i + 1, j + 1))
outstr.append('Freq [rad/s] Response ')
- outstr.append('------------ ------------------------')
+ outstr.append('------------ ---------------------')
outstr.extend(
- [ '%12.3f %10.4g + %10.4g' % (w, m, p)
+ [ '%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) ])
@@ -340,7 +340,10 @@
"""
# Preallocate the output.
- out = empty((self.outputs, self.inputs), dtype=complex)
+ if getattr(omega, '__iter__', False):
+ out = empty((self.outputs, self.inputs, len(omega)), dtype=complex)
+ else:
+ out = empty((self.outputs, self.inputs), dtype=complex)
if self.ifunc is None:
try:
@@ -350,11 +353,18 @@
"Frequency %f not in frequency list, try an interpolating"
" FRD if you want additional points")
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]
-
+ 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]
+ 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]
+
return out
# Method for generating the frequency response of the system
Modified: trunk/src/lti.py
===================================================================
--- trunk/src/lti.py 2013-05-31 11:53:29 UTC (rev 273)
+++ trunk/src/lti.py 2013-06-09 14:57:59 UTC (rev 274)
@@ -12,6 +12,8 @@
timebaseEqual()
"""
+from numpy import absolute, real
+
class Lti:
"""Lti is a parent class to linear time invariant control (LTI) objects.
@@ -64,6 +66,12 @@
self.outputs = outputs
self.dt = dt
+ def damp(self):
+ poles = self.pole()
+ wn = absolute(poles)
+ Z = -real(poles)/wn
+ return wn, Z, poles
+
# Test to see if a system is SISO
def issiso(sys, strict=False):
if isinstance(sys, (int, float, complex)) and not strict:
Modified: trunk/src/margins.py
===================================================================
--- trunk/src/margins.py 2013-05-31 11:53:29 UTC (rev 273)
+++ trunk/src/margins.py 2013-06-09 14:57:59 UTC (rev 274)
@@ -53,15 +53,37 @@
import numpy as np
import control.xferfcn as xferfcn
from control.freqplot import bode
-from control.lti import isdtime
+from control.lti import isdtime, issiso
+import control.frdata as frdata
+import scipy as sp
-# gain and phase margins
-# contributed by Sawyer B. Fuller <mi...@ca...>
-#! TODO - need to add unit test functions
-def stability_margins(sysdata, deg=True):
+# helper functions for stability_margins
+def _polyimsplit(pol):
+ """split a polynomial with (iw) applied into a real and an
+ imaginary part with w applied"""
+ rpencil = np.zeros_like(pol)
+ ipencil = np.zeros_like(pol)
+ rpencil[-1::-4] = 1.
+ rpencil[-3::-4] = -1.
+ ipencil[-2::-4] = 1.
+ ipencil[-4::-4] = -1.
+ 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 <mi...@ca...>, removed a lot of the innards
+# and replaced with analytical polynomial functions for Lti systems.
+#
+# idea for the frequency data solution copied/adapted from
+# https://github.com/alchemyst/Skogestad-Python/blob/master/BODE.py
+# Rene van Paassen <ren...@gm...>
+def stability_margins(sysdata, deg=True, returnall=False):
"""Calculate gain, phase and stability margins and associated
crossover frequencies.
-
+
Usage
-----
gm, pm, sm, wg, wp, ws = stability_margins(sysdata, deg=True)
@@ -76,86 +98,117 @@
bode frequency response data
deg=True: boolean
If true, all input and output phases in degrees, else in radians
-
+ returnall=False: boolean
+ If true, return all margins found. Note that for frequency data or
+ FRD systems, only one margin is found and returned.
+
Returns
-------
- gm, pm, sm, wg, wp, ws: float
+ gm, pm, sm, wg, wp, ws: float or array_like
Gain margin gm, phase margin pm, stability margin sm, and
associated crossover
frequencies wg, wp, and ws of SISO open-loop. If more than
one crossover frequency is detected, returns the lowest corresponding
margin.
- """
- #TODO do this precisely without the effects of discretization of frequencies?
- #TODO assumes SISO
- #TODO unit tests, margin plot
+ When requesting all margins, the return values are array_like,
+ and all margins are returns for linear systems not equal to FRD
+ """
+
+ try:
+ if isinstance(sysdata, frdata.FRD):
+ sys = frdata.FRD(sysdata, smooth=True)
+ elif isinstance(sysdata, xferfcn.TransferFunction):
+ sys = sysdata
+ elif getattr(sysdata, '__iter__', False) and len(sysdata) == 3:
+ mag, phase, omega = sysdata
+ sys = frdata.FRD(mag*np.exp((1j/360.)*phase), omega, smooth=True)
+ else:
+ sys = xferfcn._convertToTransferFunction(sysdata)
+ except Exception, e:
+ print (e)
+ raise ValueError("Margin sysdata must be either a linear system or "
+ "a 3-sequence of mag, phase, omega.")
- if (not getattr(sysdata, '__iter__', False)):
- sys = sysdata
+ # calculate gain of system
+ if isinstance(sys, xferfcn.TransferFunction):
+
+ # check for siso
+ if not issiso(sys):
+ raise ValueError("Can only do margins for SISO system")
+
+ # real and imaginary part polynomials in omega:
+ rnum, inum = _polyimsplit(sys.num[0][0])
+ rden, iden = _polyimsplit(sys.den[0][0])
- # TODO: implement for discrete time systems
- if (isdtime(sys, strict=True)):
- raise NotImplementedError("Function not implemented in discrete time")
+ # test imaginary part of tf == 0, for phase crossover/gain margins
+ test_w_180 = np.polyadd(np.polymul(inum, rden), np.polymul(rnum, -iden))
+ w_180 = np.roots(test_w_180)
+ w_180 = np.real(w_180[(np.imag(w_180) == 0) * (w_180 > 0.)])
+ w_180.sort()
- mag, phase, omega = bode(sys, deg=deg, Plot=False)
- elif len(sysdata) == 3:
- # TODO: replace with FRD object type?
- mag, phase, omega = sysdata
- else:
- raise ValueError("Margin sysdata must be either a linear system or a 3-sequence of mag, phase, omega.")
-
- if deg:
- cycle = 360.
- crossover = 180.
+ # test magnitude is 1 for gain crossover/phase margins
+ test_wc = np.polysub(np.polyadd(_polysqr(rnum), _polysqr(inum)),
+ np.polyadd(_polysqr(rden), _polysqr(iden)))
+ wc = np.roots(test_wc)
+ wc = np.real(wc[(np.imag(wc) == 0) * (wc > 0.)])
+ wc.sort()
+
+ # stability margin was a bitch to elaborate, relies on magnitude to
+ # point -1, then take the derivative. Second derivative needs to be >0
+ # to have a minimum
+ test_wstabn = np.polyadd(_polysqr(rnum), _polysqr(inum))
+ test_wstabd = 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))
+
+ # find the solutions
+ wstab = np.roots(test_wstab)
+
+ # 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 > 0.) *
+ (np.abs(wstabplus) > 0.)])
+ wstab.sort()
+
else:
- cycle = 2 * np.pi
- crossover = np.pi
-
- wrapped_phase = -np.mod(phase, cycle)
-
- # phase margin from minimum phase among all gain crossovers
- neg_mag_crossings_i = np.nonzero(np.diff(mag < 1) > 0)[0]
- mag_crossings_p = wrapped_phase[neg_mag_crossings_i]
- if len(neg_mag_crossings_i) == 0:
- if mag[0] < 1: # gain always less than one
- wp = np.nan
- pm = np.inf
- else: # gain always greater than one
- print("margin: no magnitude crossings found")
- wp = np.nan
- pm = np.nan
- else:
- min_mag_crossing_i = neg_mag_crossings_i[np.argmin(mag_crossings_p)]
- wp = omega[min_mag_crossing_i]
- pm = crossover + phase[min_mag_crossing_i]
- if pm < 0:
- print("warning: system unstable: negative phase margin")
-
- # gain margin from minimum gain margin among all phase crossovers
- neg_phase_crossings_i = np.nonzero(np.diff(wrapped_phase < -crossover) > 0)[0]
- neg_phase_crossings_g = mag[neg_phase_crossings_i]
- if len(neg_phase_crossings_i) == 0:
- wg = np.nan
- gm = np.inf
- else:
- min_phase_crossing_i = neg_phase_crossings_i[
- np.argmax(neg_phase_crossings_g)]
- wg = omega[min_phase_crossing_i]
- gm = abs(1/mag[min_phase_crossing_i])
- if gm < 1:
- print("warning: system unstable: gain margin < 1")
+ # a bit coarse, have the interpolated frd evaluated again
+ def mod(w):
+ """to give the function to calculate |G(jw)| = 1"""
+ return [np.abs(sys.evalfr(w[0])[0][0]) - 1]
- # stability margin from minimum abs distance from -1 point
- if deg:
- phase_rad = phase * np.pi / 180.
+ def arg(w):
+ """function to calculate the phase angle at -180 deg"""
+ return [np.angle(sys.evalfr(w[0])[0][0]) + np.pi]
+
+ def dstab(w):
+ """function to calculate the distance from -1 point"""
+ return np.abs(sys.evalfr(w[0])[0][0] + 1.)
+
+ # how to calculate the frequency at which |G(jw)| = 1
+ wc = np.array([sp.optimize.fsolve(mod, sys.omega[0])])[0]
+ w_180 = np.array([sp.optimize.fsolve(arg, sys.omega[0])])[0]
+ wstab = np.real(
+ np.array([sp.optimize.fmin(dstab, sys.omega[0], disp=0)])[0])
+
+ # margins, as iterables, converted frdata and xferfcn calculations to
+ # vector for this
+ PM = np.angle(sys.evalfr(wc)[0][0], deg=True) + 180
+ GM = 1/(np.abs(sys.evalfr(w_180)[0][0]))
+ SM = np.abs(sys.evalfr(wstab)[0][0]+1)
+
+ if returnall:
+ return GM, PM, SM, wc, w_180, wstab
else:
- phase_rad = phase
- L = mag * np.exp(1j * phase_rad) # complex loop response to -1 pt
- min_Lplus1_i = np.argmin(np.abs(L + 1))
- sm = np.abs(L[min_Lplus1_i] + 1)
- ws = phase[min_Lplus1_i]
+ return (
+ (GM.shape[0] or None) and GM[0],
+ (PM.shape[0] or None) and PM[0],
+ (SM.shape[0] or None) and SM[0],
+ (wc.shape[0] or None) and wc[0],
+ (w_180.shape[0] or None) and w_180[0],
+ (wstab.shape[0] or None) and wstab[0])
- return gm, pm, sm, wg, wp, ws
# Contributed by Steffen Waldherr <wal...@is...>
#! TODO - need to add test functions
Modified: trunk/src/matlab.py
===================================================================
--- trunk/src/matlab.py 2013-05-31 11:53:29 UTC (rev 273)
+++ trunk/src/matlab.py 2013-06-09 14:57:59 UTC (rev 274)
@@ -89,6 +89,7 @@
from control.statesp import StateSpace, _rss_generate, _convertToStateSpace
from control.xferfcn import TransferFunction, _convertToTransferFunction
from control.lti import Lti #base class of StateSpace, TransferFunction
+from control.frdata import FRD
from control.dtime import sample_system
from control.exception import ControlArgument
@@ -96,11 +97,11 @@
from control.ctrlutil import unwrap
from control.freqplot import nyquist, gangof4
from control.nichols import nichols
-from control.bdalg import series, parallel, negate, feedback
+from control.bdalg import series, parallel, negate, feedback, append, connect
from control.pzmap import pzmap
from control.statefbk import ctrb, obsv, gram, place, lqr
from control.delay import pade
-from control.modelsimp import hsvd, balred, modred
+from control.modelsimp import hsvd, balred, modred, minreal
from control.mateqn import lyap, dlyap, dare, care
__doc__ += r"""
@@ -125,7 +126,7 @@
\* :func:`ss` create state-space (SS) models
\ dss create descriptor state-space models
\ delayss create state-space models with delayed terms
-\ frd create frequency response data (FRD) models
+\* :func:`frd` create frequency response data (FRD) models
\ lti/exp create pure continuous-time delays (TF and
ZPK only)
\ filt specify digital filters
@@ -156,7 +157,7 @@
\* :func:`tf` conversion to transfer function
\ zpk conversion to zero/pole/gain
\* :func:`ss` conversion to state space
-\ frd conversion to frequency data
+\* :func:`frd` conversion to frequency data
\ c2d continuous to discrete conversion
\ d2c discrete to continuous conversion
\ d2d resample discrete-time model
@@ -174,7 +175,7 @@
----------------------------------------------------------------------------
== ========================== ============================================
-\ append group LTI models by appending inputs/outputs
+\* :func:`~bdalg.append` group LTI models by appending inputs/outputs
\* :func:`~bdalg.parallel` connect LTI models in parallel
(see also overloaded ``+``)
\* :func:`~bdalg.series` connect LTI models in series
@@ -200,7 +201,7 @@
\ lti/order model order (number of states)
\* :func:`~pzmap.pzmap` pole-zero map (TF only)
\ lti/iopzmap input/output pole-zero map
-\ damp natural frequency, damping of system poles
+\* :func:`damp` natural frequency, damping of system poles
\ esort sort continuous poles by real part
\ dsort sort discrete poles by magnitude
\ lti/stabsep stable/unstable decomposition
@@ -243,7 +244,7 @@
----------------------------------------------------------------------------
== ========================== ============================================
-\ minreal minimal realization; pole/zero cancellation
+\* :func:`~modelsimp.minreal` minimal realization; pole/zero cancellation
\ ss/sminreal structurally minimal realization
\* :func:`~modelsimp.hsvd` hankel singular values (state contributions)
\* :func:`~modelsimp.balred` reduced-order approximations of LTI models
@@ -569,7 +570,43 @@
else:
raise ValueError("Needs 1 or 2 arguments; received %i." % len(args))
+def frd(*args):
+ '''
+ Construct a Frequency Response Data model, or convert a system
+ frd models store the (measured) frequency response of a system.
+
+ This function can be called in different ways:
+
+ ``frd(response, freqs)``
+ Create an frd model with the given response data, in the form of
+ complex response vector, at matching frequency freqs [in rad/s]
+
+ ``frd(sys, freqs)``
+ Convert an Lti system into an frd model with data at frequencies
+ freqs.
+
+ Parameters
+ ----------
+ response: array_like, or list
+ complex vector with the system response
+ freq: array_lik or lis
+ vector with frequencies
+ sys: Lti (StateSpace or TransferFunction)
+ A linear system
+
+ Returns
+ -------
+ sys: FRD
+ New frequency response system
+
+ See Also
+ --------
+ ss, tf
+ '''
+ return FRD(*args)
+
+
def ss2tf(*args):
"""
Transform a state space system to a transfer function.
@@ -714,7 +751,7 @@
else:
raise ValueError("Needs 1 or 2 arguments; received %i." % len(args))
-def rss(states=1, inputs=1, outputs=1):
+def rss(states=1, outputs=1, inputs=1):
"""
Create a stable **continuous** random state space object.
@@ -751,7 +788,7 @@
return _rss_generate(states, inputs, outputs, 'c')
-def drss(states=1, inputs=1, outputs=1):
+def drss(states=1, outputs=1, inputs=1):
"""
Create a stable **discrete** random state space object.
@@ -852,16 +889,20 @@
return sys.zero()
-def evalfr(sys, omega):
+def evalfr(sys, x):
"""
- Evaluate the transfer function of an LTI system at an angular frequency.
+ 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
- omega: scalar
- Frequency
+ x: scalar
+ Complex number
Returns
-------
@@ -880,15 +921,17 @@
Examples
--------
>>> sys = ss("1. -2; 3. -4", "5.; 7", "6. 8", "9.")
- >>> evalfr(sys, 1.)
+ >>> 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)
+
- return sys.evalfr(omega)
-
def freqresp(sys, omega):
"""
Frequency response of an LTI system at multiple angular frequencies.
@@ -1057,10 +1100,11 @@
Returns
-------
- gm, pm, wg, wp : float
- Gain margin gm, phase margin pm (in deg), and associated crossover
- frequencies wg and wp (in rad/sec) of SISO open-loop. If more than
- one crossover frequency is detected, returns the lowest
+ gm, pm, Wcg, Wcp : float
+ Gain margin gm, phase margin pm (in deg), gain crossover frequency
+ (corresponding to phase margin) and phase crossover frequency
+ (corresponding to gain margin), in rad/sec of SISO open-loop.
+ If more than one crossover frequency is detected, returns the lowest
corresponding margin.
Examples
@@ -1083,7 +1127,7 @@
raise ValueError("Margin needs 1 or 3 arguments; received %i."
% len(args))
- return margin[0], margin[1], margin[3], margin[4]
+ return margin[0], margin[1], margin[4], margin[3]
def dcgain(*args):
'''
@@ -1139,10 +1183,48 @@
gain = sys.D - sys.C * sys.A.I * sys.B
return gain
+def damp(sys, doprint=True):
+ '''
+ Compute natural frequency, damping 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
+
+ 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
+
# Simulation routines
# Call corresponding functions in timeresp, with arguments transposed
-def step(sys, T=None, X0=0., input=0, output=0, **keywords):
+def step(sys, T=None, X0=0., input=0, output=None, **keywords):
'''
Step response of a linear system
@@ -1192,7 +1274,7 @@
Examples
--------
- >>> T, yout = step(sys, T, X0)
+ >>> yout, T = step(sys, T, X0)
'''
T, yout = timeresp.step_response(sys, T, X0, input, output,
transpose = True, **keywords)
@@ -1405,3 +1487,4 @@
# TODO: add docstring
# Call the sample_system() function to do the work
return sample_system(sysc, Ts, method)
+
Modified: trunk/src/modelsimp.py
===================================================================
--- trunk/src/modelsimp.py 2013-05-31 11:53:29 UTC (rev 273)
+++ trunk/src/modelsimp.py 2013-06-09 14:57:59 UTC (rev 274)
@@ -269,6 +269,33 @@
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
+ characteristics as the original model sys.
+
+ Parameters
+ ----------
+ sys: StateSpace or TransferFunction
+ Original system
+ tol: real
+ Tolerance
+ verbose: bool
+ Print results if True
+
+ Returns
+ -------
+ 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())))
+ return sysr
+
def era(YY, m, n, nin, nout, r):
"""
Calculate an ERA model of order `r` based on the impulse-response data `YY`.
Modified: trunk/src/rlocus.py
===================================================================
--- trunk/src/rlocus.py 2013-05-31 11:53:29 UTC (rev 273)
+++ trunk/src/rlocus.py 2013-06-09 14:57:59 UTC (rev 274)
@@ -50,9 +50,11 @@
import scipy.signal # signal processing toolbox
import pylab # plotting routines
import control.xferfcn as xferfcn
+from functools import partial
# Main function: compute a root locus diagram
-def root_locus(sys, kvect, xlim=None, ylim=None, plotstr='-', Plot=True):
+def root_locus(sys, kvect, xlim=None, ylim=None, plotstr='-', Plot=True,
+ PrintGain=True):
"""Calculate the root locus by finding the roots of 1+k*TF(s)
where TF is self.num(s)/self.den(s) and each k is an element
of kvect.
@@ -65,7 +67,9 @@
List of gains to use in computing diagram
Plot : boolean (default = True)
If True, plot magnitude and phase
-
+ PrintGain: boolean (default = True)
+ If True, report mouse clicks when close to the root-locus branches,
+ calculate gain, damping and print
Return values
-------------
rlist : list of computed root locations
@@ -80,6 +84,10 @@
# Create the plot
if (Plot):
+ f = pylab.figure()
+ if PrintGain:
+ cid = f.canvas.mpl_connect(
+ 'button_release_event', partial(_RLFeedbackClicks, sys=sys))
ax = pylab.axes();
# plot open loop poles
@@ -165,3 +173,12 @@
sorted[n,ind] = elem
prevrow = sorted[n,:]
return sorted
+
+def _RLFeedbackClicks(event, sys):
+ """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:
+ print("Clicked at %10.4g%+10.4gj gain %10.4g damp %10.4g" %
+ (s.real, s.imag, K.real, -1*s.real/abs(s)))
Modified: trunk/src/statesp.py
===================================================================
--- trunk/src/statesp.py 2013-05-31 11:53:29 UTC (rev 273)
+++ trunk/src/statesp.py 2013-06-09 14:57:59 UTC (rev 274)
@@ -27,6 +27,7 @@
StateSpace.zero
StateSpace.feedback
StateSpace.returnScipySignalLti
+StateSpace.append
_convertToStateSpace
_rss_generate
@@ -338,10 +339,23 @@
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
- else:
- raise TypeError("can't interconnect systems")
+ # 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, 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."""
@@ -370,10 +384,16 @@
else:
s = omega * 1.j
- fresp = self.C * solve(s * eye(self.states) - self.A,
- self.B) + self.D
+ return self.horner(s)
- return array(fresp)
+ 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
# TODO: add discrete time check
@@ -475,6 +495,17 @@
return StateSpace(A, B, C, D, dt)
+ def minreal(self, tol=None):
+ """Calculate a minimal realization, removes unobservable and
+ uncontrollable states"""
+ try:
+ from slycot import tb01pd
+ A, B, C, nr = tb01pd(self.states, self.inputs, self.outputs,
+ self.A, self.B, self.C, tol=tol)
+ return StateSpace(A[:nr,:nr], B[:nr,:], C[:,:nr], self.D)
+ except ImportError:
+ raise TypeError("minreal requires slycot tb01pd")
+
# TODO: add discrete time check
def returnScipySignalLti(self):
"""Return a list of a list of scipy.signal.lti objects.
@@ -497,6 +528,33 @@
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)
+
# TODO: add discrete time check
def _convertToStateSpace(sys, **kw):
"""Convert a system to state space form (if needed).
@@ -757,3 +815,51 @@
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
Modified: trunk/src/timeresp.py
===================================================================
--- trunk/src/timeresp.py 2013-05-31 11:53:29 UTC (rev 273)
+++ trunk/src/timeresp.py 2013-06-09 14:57:59 UTC (rev 274)
@@ -121,7 +121,7 @@
from copy import deepcopy
import warnings
from control.lti import Lti # base class of StateSpace, TransferFunction
-from control. statesp import StateSpace, _rss_generate, _convertToStateSpace, _mimo2siso
+from control. statesp import StateSpace, _rss_generate, _convertToStateSpace, _mimo2simo, _mimo2siso
from control.lti import isdtime, isctime
# Helper function for checking array-like parameters
@@ -386,15 +386,15 @@
return T, yout, xout
-def step_response(sys, T=None, X0=0., input=0, output=0, \
- transpose = False, **keywords):
+def step_response(sys, T=None, X0=0., input=0, output=None,
+ transpose = False, **keywords):
#pylint: disable=W0622
"""Step response of a linear system
- If the system has multiple inputs or outputs (MIMO), one input and one
- output have to be selected for the simulation. The parameters `input`
- and `output` do this. All other inputs are set to 0, all other outputs
- are ignored.
+ If the system has multiple inputs or outputs (MIMO), one input has
+ to be selected for the simulation. Optionally, one output may be
+ selected. The parameters `input` and `output` do this. All other
+ inputs are set to 0, all other outputs are ignored.
For information on the **shape** of parameters `T`, `X0` and
return values `T`, `yout` see: :ref:`time-series-convention`
@@ -416,7 +416,8 @@
Index of the input that will be used in this simulation.
output: int
- Index of the output that will be used in this simulation.
+ Index of the output that will be used in this simulation. Set to None
+ to not trim outputs
transpose: bool
If True, transpose all input and output arrays (for backward
@@ -447,7 +448,10 @@
>>> T, yout = step_response(sys, T, X0)
"""
sys = _convertToStateSpace(sys)
- sys = _mimo2siso(sys, input, output, warn_conversion=True)
+ if output == None:
+ sys = _mimo2simo(sys, input, warn_conversion=True)
+ else:
+ sys = _mimo2siso(sys, input, output, warn_conversion=True)
if T is None:
if isctime(sys):
T = _default_response_times(sys.A, 100)
@@ -459,13 +463,13 @@
U = np.ones_like(T)
T, yout, _xout = forced_response(sys, T, U, X0,
- transpose=transpose, **keywords)
+ transpose=transpose, **keywords)
return T, yout
-def initial_response(sys, T=None, X0=0., input=0, output=0, transpose=False,
- **keywords):
+def initial_response(sys, T=None, X0=0., input=0, output=None, transpose=False,
+ **keywords):
#pylint: disable=W0622
"""Initial condition response of a linear system
@@ -494,7 +498,8 @@
Index of the input that will be used in this simulation.
output: int
- Index of the output that will be used in this simulation.
+ Index of the output that will be used in this simulation. Set to None
+ to not trim outputs
transpose: bool
If True, transpose all input and output arrays (for backward
@@ -525,7 +530,10 @@
>>> T, yout = initial_response(sys, T, X0)
"""
sys = _convertToStateSpace(sys)
- sys = _mimo2siso(sys, input, output, warn_conversion=True)
+ if output == None:
+ sys = _mimo2simo(sys, input, warn_conversion=False)
+ else:
+ sys = _mimo2siso(sys, input, output, warn_conversion=False)
# Create time and input vectors; checking is done in forced_response(...)
# The initial vector X0 is created in forced_response(...) if necessary
@@ -534,11 +542,11 @@
U = np.zeros_like(T)
T, yout, _xout = forced_response(sys, T, U, X0, transpose=transpose,
- **keywords)
+ **keywords)
return T, yout
-def impulse_response(sys, T=None, X0=0., input=0, output=0,
+def impulse_response(sys, T=None, X0=0., input=0, output=None,
transpose=False, **keywords):
#pylint: disable=W0622
"""Impulse response of a linear system
@@ -568,7 +576,8 @@
Index of the input that will be used in this simulation.
output: int
- Index of the output that will be used in this simulation.
+ Index of the output that will be used in this simulation. Set to None
+ to not trim outputs
transpose: bool
If True, transpose all input and output arrays (for backward
@@ -599,7 +608,10 @@
>>> T, yout = impulse_response(sys, T, X0)
"""
sys = _convertToStateSpace(sys)
- sys = _mimo2siso(sys, input, output, warn_conversion=True)
+ if output == None:
+ sys = _mimo2simo(sys, input, warn_conversion=True)
+ else:
+ sys = _mimo2siso(sys, input, output, warn_conversion=True)
# System has direct feedthrough, can't simulate impulse response numerically
if np.any(sys.D != 0) and isctime(sys):
Modified: trunk/src/xferfcn.py
===================================================================
--- trunk/src/xferfcn.py 2013-05-31 11:53:29 UTC (rev 273)
+++ trunk/src/xferfcn.py 2013-06-09 14:57:59 UTC (rev 274)
@@ -404,7 +404,8 @@
for k in range(self.inputs): # Multiply & add.
num_summand[k] = polymul(self.num[i][k], other.num[k][j])
den_summand[k] = polymul(self.den[i][k], other.den[k][j])
- num[i][j], den[i][j] = _addSISO(num[i][j], den[i][j],
+ num[i][j], den[i][j] = _addSISO(
+ num[i][j], den[i][j],
num_summand[k], den_summand[k])
return TransferFunction(num, den, dt)
@@ -412,8 +413,48 @@
def __rmul__(self, other):
"""Right multiply two LTI objects (serial connection)."""
- return self * other
+ # Convert the second argument to a transfer function.
+ if isinstance(other, (int, float, complex)):
+ other = _convertToTransferFunction(other, inputs=self.inputs,
+ outputs=self.inputs)
+ else:
+ other = _convertToTransferFunction(other)
+
+ # Check that the input-output sizes are consistent.
+ if other.inputs != self.outputs:
+ raise ValueError("C = A * B: A has %i column(s) (input(s)), but B \
+has %i row(s)\n(output(s))." % (other.inputs, self.outputs))
+ inputs = self.inputs
+ outputs = 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 (self.dt == other.dt):
+ dt = self.dt # use dt from first argument
+ else:
+ 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)]
+
+ # Temporary storage for the summands needed to find the (i, j)th element
+ # of the product.
+ num_summand = [[] for k in range(other.inputs)]
+ den_summand = [[] for k in range(other.inputs)]
+
+ for i in range(outputs): # Iterate through rows of product.
+ for j in range(inputs): # Iterate through columns of product.
+ for k in range(other.inputs): # Multiply & add.
+ num_summand[k] = polymul(other.num[i][k], self.num[k][j])
+ den_summand[k] = polymul(other.den[i][k], self.den[k][j])
+ num[i][j], den[i][j] = _addSISO(num[i][j], den[i][j],
+ num_summand[k], den_summand[k])
+
+ return TransferFunction(num, den, dt)
+
# TODO: Division of MIMO transfer function objects is not written yet.
def __div__(self, other):
"""Divide two LTI objects."""
@@ -486,9 +527,20 @@
warn("evalfr: frequency evaluation above Nyquist frequency")
else:
s = 1.j * omega
+
+ 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.
+ '''
+
# Preallocate the output.
- out = empty((self.outputs, self.inputs), dtype=complex)
+ if getattr(s, '__iter__', False):
+ out = empty((self.outputs, self.inputs, len(s)), dtype=complex)
+ else:
+ out = empty((self.outputs, self.inputs), dtype=complex)
for i in range(self.outputs):
for j in range(self.inputs):
@@ -620,8 +672,12 @@
newzeros.append(z)
# keep result
- num[i][j] = gain * real(poly(newzeros))
+ if len(newzeros):
+ num[i][j] = gain * real(poly(newzeros))
+ else:
+ num[i][j] = array([gain])
den[i][j] = real(poly(poles))
+
# end result
return TransferFunction(num, den)
@@ -906,6 +962,15 @@
In the latter example, sys's matrix transfer function is [[1., 1., 1.]
[1., 1., 1.]].
+
+ If sys is an array-like type, then it is converted to a constant-gain
+ transfer function.
+
+ >>> sys = _convertToTransferFunction([[1. 0.], [2. 3.]])
+
+ In this example, the numerator matrix will be
+ [[[1.0], [0.0]], [[2.0], [3.0]]]
+ and the denominator matrix [[[1.0], [1.0]], [[1.0], [1.0]]]
"""
from control.statesp import StateSpace
@@ -966,5 +1031,16 @@
den = [[[1] for j in range(inputs)] for i in range(outputs)]
return TransferFunction(num, den)
- else:
- raise TypeError("Can't convert given type to TransferFunction system.")
+
+ # If this is array-like, try to create a constant feedthrough
+ try:
+ D = array(sys)
+ outputs, inputs = D.shape
+ num = [[[D[i,j]] for j in range(inputs)] for i in range(outputs)]
+ den = [[[1] for j in range(inputs)] for i in range(outputs)]
+ return TransferFunction(num, den)
+ except Exception, e:
+ print("Failure to assume argument is matrix-like in"
+ " _convertToTransferFunction, result %s" % e)
+
+ raise TypeError("Can't convert given type to TransferFunction system.")
Modified: trunk/tests/convert_test.py
===================================================================
--- trunk/tests/convert_test.py 2013-05-31 11:53:29 UTC (rev 273)
+++ trunk/tests/convert_test.py 2013-06-09 14:57:59 UTC (rev 274)
@@ -57,7 +57,7 @@
for outputs in range(1, self.maxIO):
# start with a random SS system and transform to TF then
# back to SS, check that the matrices are the same.
- ssOriginal = matlab.rss(states, inputs, outputs)
+ ssOriginal = matlab.rss(states, outputs, inputs)
if (verbose):
self.printSys(ssOriginal, 1)
Modified: trunk/tests/matlab_test.py
===================================================================
--- trunk/tests/matlab_test.py 2013-05-31 11:53:29 UTC (rev 273)
+++ trunk/tests/matlab_test.py 2013-06-09 14:57:59 UTC (rev 274)
@@ -11,9 +11,48 @@
from __future__ import print_function
import unittest
import numpy as np
+from scipy.linalg import eigvals
import scipy as sp
from control.matlab import *
+from control.frdata import FRD
+# for running these through Matlab or Octave
+'''
+siso_ss1 = ss([1. -2.; 3. -4.], [5.; 7.], [6. 8.], [0])
+
+siso_tf1 = tf([1], [1, 2, 1])
+siso_tf2 = tf([1, 1], [1, 2, 3, 1])
+
+siso_tf3 = tf(siso_ss1)
+siso_ss2 = ss(siso_tf2)
+siso_ss3 = ss(siso_tf3)
+siso_tf4 = tf(siso_ss2)
+
+A =[ 1. -2. 0. 0.;
+ 3. -4. 0. 0.;
+ 0. 0. 1. -2.;
+ 0. 0. 3. -4. ]
+B = [ 5. 0.;
+ 7. 0.;
+ 0. 5.;
+ 0. 7. ]
+C = [ 6. 8. 0. 0.;
+ 0. 0. 6. 8. ]
+D = [ 9. 0.;
+ 0. 9. ]
+mimo_ss1 = ss(A, B, C, D)
+
+% all boring, since no cross-over
+margin(siso_tf1)
+margin(siso_tf2)
+margin(siso_ss1)
+margin(siso_ss2)
+
+% make a bit better
+[gm, pm, gmc, pmc] = margin(siso_ss2*siso_ss2*2)
+
+'''
+
class TestMatlab(unittest.TestCase):
def setUp(self):
"""Set up some systems for testing out MATLAB functions"""
@@ -194,6 +233,9 @@
gm, pm, wg, wp = margin(self.siso_tf2);
gm, pm, wg, wp = margin(self.siso_ss1);
gm, pm, wg, wp = margin(self.siso_ss2);
+ gm, pm, wg, wp = margin(self.siso_ss2*self.siso_ss2*2);
+ np.testing.assert_array_almost_equal(
+ [gm, pm, wg, wp], [1.5451, 75.9933, 1.2720, 0.6559], decimal=3)
def testDcgain(self):
#Create different forms of a SISO system
@@ -274,13 +316,16 @@
freqresp(self.siso_tf3, w)
def testEvalfr(self):
- w = 1
- evalfr(self.siso_ss1, w)
+ w = 1j
+ self.assertEqual(evalfr(self.siso_ss1, w), 44.8-21.4j)
evalfr(self.siso_ss2, w)
evalfr(self.siso_ss3, w)
evalfr(self.siso_tf1, w)
evalfr(self.siso_tf2, w)
evalfr(self.siso_tf3, w)
+ np.testing.assert_array_almost_equal(
+ evalfr(self.mimo_ss1, w),
+ np.array( [[44.8-21.4j, 0.], [0., 44.8-21.4j]]))
def testHsvd(self):
hsvd(self.siso_ss1)
@@ -309,12 +354,12 @@
def testRss(self):
rss(1)
rss(2)
- rss(2, 3, 1)
+ rss(2, 1, 3)
def testDrss(self):
drss(1)
drss(2)
- drss(2, 3, 1)
+ drss(2, 1, 3)
def testCtrb(self):
ctrb(self.siso_ss1.A, self.siso_ss1.B)
@@ -372,6 +417,102 @@
for i in range(len(tfdata_1)):
np.testing.assert_array_almost_equal(tfdata_1[i], tfdata_2[i])
+ def testDamp(self):
+ A = np.mat('''-0.2 0.06 0 -1;
+ 0 0 1 0;
+ -17 0 -3.8 1;
+ 9.4 0 -0.4 -0.6''')
+ B = np.mat('''-0.01 0.06;
+ 0 0;
+ -32 5.4;
+ 2.6 -7''')
+ C = np.eye(4)
+ D = np.zeros((4,2))
+ sys = ss(A, B, C, D)
+ wn, Z, p = damp(sys, False)
+ print (wn)
+ np.testing.assert_array_almost_equal(
+ wn, np.array([4.07381994, 3.28874827, 3.28874827,
+ 1.08937685e-03]))
+ np.testing.assert_array_almost_equal(
+ Z, np.array([1.0, 0.07983139, 0.07983139, 1.0]))
+
+ def testConnect(self):
+ sys1 = ss("1. -2; 3. -4", "5.; 7", "6, 8", "9.")
+ sys2 = ss("-1.", "1.", "1.", "0.")
+ sys = append(sys1, sys2)
+ Q= np.mat([ [ 1, 2], [2, -1] ]) # basically feedback, output 2 in 1
+ sysc = connect(sys, Q, [2], [1, 2])
+ print(sysc)
+ np.testing.assert_array_almost_equal(
+ sysc.A, np.mat('1 -2 5; 3 -4 7; -6 -8 -10'))
+ np.testing.assert_array_almost_equal(
+ sysc.B, np.mat('0; 0; 1'))
+ np.testing.assert_array_almost_equal(
+ sysc.C, np.mat('6 8 9; 0 0 1'))
+ np.testing.assert_array_almost_equal(
+ sysc.D, np.mat('0; 0'))
+
+ def testConnect2(self):
+ sys = append(ss([[-5, -2.25], [4, 0]], [[2], [0]],
+ [[0, 1.125]], [[0]]),
+ ss([[-1.6667, 0], [1, 0]], [[2], [0]],
+ [[0, 3.3333]], [[0]]),
+ 1)
+ Q = [ [ 1, 3], [2, 1], [3, -2]]
+ sysc = connect(sys, Q, [3], [3, 1, 2])
+ np.testing.assert_array_almost_equal(
+ sysc.A, np.mat([[-5, -2.25, 0, -6.6666],
+ [4, 0, 0, 0],
+ [0, 2.25, -1.6667, 0],
+ [0, 0, 1, 0]]))
+ np.testing.assert_array_almost_equal(
+ sysc.B, np.mat([[2], [0], [0], [0]]))
+ np.testing.assert_array_almost_equal(
+ sysc.C, np.mat([[0, 0, 0, -3.3333],
+ [0, 1.125, 0, 0],
+ [0, 0, 0, 3.3333]]))
+ 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)
+ assert isinstance(frd2, FRD)
+
+ def testMinreal(self, verbose=False):
+ """Test a minreal model reduction"""
+ #A = [-2, 0.5, 0; 0.5, -0.3, 0; 0, 0, -0.1]
+ A = [[-2, 0.5, 0], [0.5, -0.3, 0], [0, 0, -0.1]]
+ #B = [0.3, -1.3; 0.1, 0; 1, 0]
+ B = [[0.3, -1.3], [0.1, 0.], [1.0, 0.0]]
+ #C = [0, 0.1, 0; -0.3, -0.2, 0]
+ C = [[0., 0.1, 0.0], [-0.3, -0.2, 0.0]]
+ #D = [0 -0.8; -0.3 0]
+ D = [[0., -0.8], [-0.3, 0.]]
+ # sys = ss(A, B, C, D)
+
+ sys = ss(A, B, C, D)
+ sysr = minreal(sys)
+ self.assertEqual(sysr.states, 2)
+ self.assertEqual(sysr.inputs, sys.inputs)
+ self.assertEqual(sysr.outputs, sys.outputs)
+ np.testing.assert_array_almost_equal(
+ eigvals(sysr.A), [-2.136154, -0.1638459])
+
+ s = tf([1, 0], [1])
+ h = (s+1)*(s+2.00000000001)/(s+2)/(s**2+s+1)
+ hm = minreal(h)
+ hr = (s+1)/(s**2+s+1)
+ np.testing.assert_array_almost_equal(hm.num[0][0], hr.num[0][0])
+ np.testing.assert_array_almost_equal(hm.den[0][0], hr.den[0][0])
+
+
#! TODO: not yet implemented
# def testMIMOtfdata(self):
# sisotf = ss2tf(self.siso_ss1)
Modified: trunk/tests/slycot_convert_test.py
===================================================================
--- trunk/tests/slycot_convert_test.py 2013-05-31 11:53:29 UTC (rev 273)
+++ trunk/tests/slycot_convert_test.py 2013-06-09 14:57:59 UTC (rev 274)
@@ -37,7 +37,7 @@
for inputs in range(1, self.maxI+1):
for outputs in range(1, self.maxO+1):
for testNum in range(self.numTests):
- ssOriginal = matlab.rss(states, inputs, outputs)
+ ssOriginal = matlab.rss(states, outputs, inputs)
if (verbose):
print('====== Original SS ==========')
print(ssOriginal)
@@ -82,7 +82,7 @@
for testNum in range(self.numTests):
for inputs in range(1,1):
for outputs in range(1,1):
- ssOriginal = matlab.rss(states, inputs, outputs)
+ ssOriginal = matlab.rss(states, outputs, inputs)
tfOriginal_Actrb, tfOriginal_Bctrb, tfOriginal_Cctrb, tfOrigingal_nctrb, tfOriginal_index,\
tfOriginal_dcoeff, tfOriginal_ucoeff = tb04ad(states,inputs,outputs,\
Modified: trunk/tests/statesp_test.py
===================================================================
--- trunk/tests/statesp_test.py 2013-05-31 11:53:29 UTC (rev 273)
+++ trunk/tests/statesp_test.py 2013-06-09 14:57:59 UTC (rev 274)
@@ -5,8 +5,10 @@
import unittest
import numpy as np
+from scipy.linalg import eigvals
import control.matlab as matlab
-from control.statesp import StateSpace
+from control.statesp import StateSpace, _convertToStateSpace
+from control.xferfcn import TransferFunction
class TestStateSpace(unittest.TestCase):
"""Tests for the StateSpace class."""
@@ -134,7 +136,74 @@
np.testing.assert_almost_equal(mag, truemag)
np.testing.assert_almost_equal(phase, truephase)
np.testing.assert_equal(omega, trueomega)
+
+ def testMinreal(self):
+ """Test a minreal model reduction"""
+ #A = [-2, 0.5, 0; 0.5, -0.3, 0; 0, 0, -0.1]
+ A = [[-2, 0.5, 0], [0.5, -0.3, 0], [0, 0, -0.1]]
+ #B = [0.3, -1.3; 0.1, 0; 1, 0]
+ B = [[0.3, -1.3], [0.1, 0.], [1.0, 0.0]]
+ #C = [0, 0.1, 0; -0.3, -0.2, 0]
+ C = [[0., 0.1, 0.0], [-0.3, -0.2, 0.0]]
+ #D = [0 -0.8; -0.3 0]
+ D = [[0., -0.8], [-0.3, 0.]]
+ # sys = ss(A, B, C, D)
+
+ sys = StateSpace(A, B, C, D)
+ sysr = sys.minreal()
+ self.assertEqual(sysr.states, 2)
+ self.assertEqual(sysr.inputs, sys.inputs)
+ self.assertEqual(sysr.outputs, sys.outputs)
+ np.testing.assert_array_almost_equal(
+ eigvals(sysr.A), [-2.136154, -0.1638459])
+
+ def testAppendSS(self):
+ """Test appending two state-space systems"""
+ A1 = [[-2, 0.5, 0], [0.5, -0.3, 0], [0, 0, -0.1]]
+ B1 = [[0.3, -1.3], [0.1, 0.], [1.0, 0.0]]
+ C1 = [[0., 0.1, 0.0], [-0.3, -0.2, 0.0]]
+ D1 = [[0., -0.8], [-0.3, 0.]]
+ A2 = [[-1.]]
+ B2 = [[1.2]]
+ C2 = [[0.5]]
+ D2 = [[0.4]]
+ A3 = [[-2, 0.5, 0, 0], [0.5, -0.3, 0, 0], [0, 0, -0.1, 0],
+ [0, 0, 0., -1.]]
+ B3 = [[0.3, -1.3, 0], [0.1, 0., 0], [1.0, 0.0, 0], [0., 0, 1.2]]
+ C3 = [[0., 0.1, 0.0, 0.0], [-0.3, -0.2, 0.0, 0.0], [0., 0., 0., 0.5]]
+ D3 = [[0., -0.8, 0.], [-0.3, 0., 0.], [0., 0., 0.4]]
+ sys1 = StateSpace(A1, B1, C1, D1)
+ sys2 = StateSpace(A2, B2, C2, D2)
+ sys3 = StateSpace(A3, B3, C3, D3)
+ sys3c = sys1.append(sys2)
+ np.testing.assert_array_almost_equal(sys3.A, sys3c.A)
+ np.testing.assert_array_almost_equal(sys3.B, sys3c.B)
+ np.testing.assert_array_almost_equal(sys3.C, sys3c.C)
+ np.testing.assert_array_almost_equal(sys3.D, sys3c.D)
+
+ def testAppendTF(self):
+ """Test appending a state-space system with a tf"""
+ A1 = [[-2, 0.5, 0], [0.5, -0.3, 0], [0, 0, -0.1]]
+ B1 = [[0.3...
[truncated message content] |
