8000 Implement c2d without slycot (fixes #23) · MorS25/python-control@03c6637 · GitHub
[go: up one dir, main page]
Skip to content

Commit 03c6637

Browse files
cwrowleycwrowley
committed
Implement c2d without slycot (fixes python-control#23)
Conversion from continuous-time to discrete-time are now implemented as methods of the TransferFunction or StateSpace classes. The function sample_system() is still provided, and takes either a state-space system or transfer function as an argument. The implementations simply call the routine `cont2discrete` from `scipy.signal` This new organization into methods avoids some messy type checking in sample_system(), and also avoids some unnecessary conversions between state-space systems and transfer functions, which fixes python-control#23.
1 parent 40c9ccd commit 03c6637

File tree

4 files changed

+166
-108
lines changed

4 files changed

+166
-108
lines changed

‎control/dtime.py

Lines changed: 4 additions & 89 deletions
Original file line numberDiff line numberDiff line change
@@ -5,7 +5,6 @@
55
Routines in this module:
66
77
sample_system()
8-
_c2dmatched()
98
"""
109

1110
"""Copyright (c) 2012 by California Institute of Technology
@@ -47,16 +46,10 @@
4746
4847
"""
4948

50-
from scipy.signal import zpk2tf, tf2zpk
51-
import numpy as np
52-
from cmath import exp
53-
from warnings import warn
5449
from .lti import isctime
55-
from .statesp import StateSpace, _convertToStateSpace
56-
from .xferfcn import TransferFunction, _convertToTransferFunction
5750

5851
# Sample a continuous time system
59-
def sample_system(sysc, Ts, method='matched'):
52+
def sample_system(sysc, Ts, method='zoh', alpha=None):
6053
"""Convert a continuous time system to discrete time
6154
6255
Creates a discrete time system from a continuous time system by
@@ -78,11 +71,8 @@ def sample_system(sysc, Ts, method='matched'):
7871
7972
Notes
8073
-----
81-
1. The conversion methods 'tustin' and 'zoh' require the
82-
cont2discrete() function, including in SciPy 0.10.0 and above.
83-
84-
2. Additional methods 'foh' and 'impulse' are planned for future
85-
implementation.
74+
See `TransferFunction.sample` and `StateSpace.sample` for
75+
further details.
8676
8777
Examples
8878
--------
@@ -94,79 +84,4 @@ def sample_system(sysc, Ts, method='matched'):
9484
if not isctime(sysc):
9585
raise ValueError("First argument must be continuous time system")
9686

97-
# If we are passed a state space system, convert to transfer function first
98-
if isinstance(sysc, StateSpace) and method == 'zoh':
99-
try:
100-
# try with slycot routine
101-
from slycot import mb05nd
102-
F, H = mb05nd(sysc.A, Ts)
103-
return StateSpace(F, H*sysc.B, sysc.C, sysc.D, Ts)
104-
except ImportError:
105-
if sysc.inputs != 1 or sysc.outputs != 1:
106-
raise TypeError(
107-
"mb05nd not found in slycot, or slycot not installed")
108-
109-
# TODO: implement MIMO version for other than ZOH state-space
110-
if (sysc.inputs != 1 or sysc.outputs != 1):
111-
raise NotImplementedError("MIMO implementation not available")
112-
113-
# SISO state-space, with other than ZOH, or failing slycot import,
114-
# is handled by conversion to TF
115-
if isinstance(sysc, StateSpace):
116-
warn("sample_system: converting to transfer function")
117-
sysc = _convertToTransferFunction(sysc)
118-
119-
# Decide what to do based on the methods available
120-
if method == 'matched':
121-
sysd = _c2dmatched(sysc, Ts)
122-
123-
elif method == 'tustin':
124-
try:
125-
from scipy.signal import cont2discrete
126-
sys = [sysc.num[0][0], sysc.den[0][0]]
127-
scipySysD = cont2discrete(sys, Ts, method='bilinear')
128-
sysd = TransferFunction(scipySysD[0][0], scipySysD[1], Ts)
129-
except ImportError:
130-
raise TypeError("cont2discrete not found in scipy.signal; upgrade to v0.10.0+")
131-
132-
elif method == 'zoh':
133-
try:
134-
from scipy.signal import cont2discrete
135-
sys = [sysc.num[0][0], sysc.den[0][0]]
136-
scipySysD = cont2discrete(sys, Ts, method='zoh')
137-
sysd = TransferFunction(scipySysD[0][0],scipySysD[1], Ts)
138-
except ImportError:
139-
raise TypeError("cont2discrete not found in scipy.signal; upgrade to v0.10.0+")
140-
141-
elif method == 'foh' or method == 'impulse':
142-
raise ValueError("Method not developed yet")
143-
144-
else:
145-
raise ValueError("Invalid discretization method: %s" % method)
146-
147-
# TODO: Convert back into the input form
148-
# Set sampling time
149-
return sysd
150-
151-
# c2d function contributed by Benjamin White, Oct 2012
152-
def _c2dmatched(sysC, Ts):
153-
# Pole-zero match method of continuous to discrete time conversion
154-
szeros, spoles, sgain = tf2zpk(sysC.num[0][0], sysC.den[0][0])
155-
zzeros = [0] * len(szeros)
156-
zpoles = [0] * len(spoles)
157-
pregainnum = [0] * len(szeros)
158-
pregainden = [0] * len(spoles)
159-
for idx, s in enumerate(szeros):
160-
sTs = s*Ts
161-
z = exp(sTs)
162-
zzeros[idx] = z
163-
pregainnum[idx] = 1-z
164-
for idx, s in enumerate(spoles):
165-
sTs = s*Ts
166-
z = exp(sTs)
167-
zpoles[idx] = z
168-
pregainden[idx] = 1-z
169-
zgain = np.multiply.reduce(pregainnum)/np.multiply.reduce(pregainden)
170-
gain = sgain/zgain
171-
sysDnum, sysDden = zpk2tf(zzeros, zpoles, gain)
172-
return TransferFunction(sysDnum, sysDden, Ts)
87+
return sysc.sample(Ts, method, alpha)

‎control/statesp.py

Lines changed: 47 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -83,7 +83,7 @@
8383
from numpy.random import rand, randn
8484
from numpy.linalg import inv, det, solve
8585
from numpy.linalg.linalg import LinAlgError
86-
from scipy.signal import lti
86+
from scipy.signal import lti, cont2discrete
8787
# from exceptions import Exception
8888
import warnings
8989
from .lti import Lti, timebase, timebaseEqual, isdtime
@@ -576,6 +576,52 @@ def __getitem__(self, indices):
576576
self.C[i,:],
577577
self.D[i,j], self.dt)
578578

579+
def sample(self, Ts, method='zoh', alpha=None):
580+
"""Convert a continuous time system to discrete time
581+
582+
Creates a discrete-time system from a continuous-time system by
583+
sampling. Multiple methods of conversion are supported.
584+
585+
Parameters
586+
----------
587+
Ts : float
588+
Sampling period
589+
method : {"gbt", "bilinear", "euler", "backward_diff", "zoh"}
590+
Which method to use:
591+
592+
* gbt: generalized bilinear transformation
593+
* bilinear: Tustin's approximation ("gbt" with alpha=0.5)
594+
* euler: Euler (or forward differencing) method ("gbt" with alpha=0)
595+
* backward_diff: Backwards differencing ("gbt" with alpha=1.0)
596+
* zoh: zero-order hold (default)
597+
598+
alpha : float within [0, 1]
599+
The generalized bilinear transformation weighting parameter, which
600+
should only be specified with method="gbt", and is ignored otherwise
601+
602+
Returns
603+
-------
604+
sysd : StateSpace system
605+
Discrete time system, with sampling rate Ts
606+
607+
Notes
608+
-----
609+
Uses the command 'cont2discrete' from scipy.signal
610+
611+
Examples
612+
--------
613+
>>> sys = StateSpace(0, 1, 1, 0)
614+
>>> sysd = sys.sample(0.5, method='bilinear')
615+
616+
"""
617+
if not self.isctime():
618+
raise ValueError("System must be continuous time system")
619+
620+
sys = (self.A, self.B, self.C, self.D)
621+
Ad, Bd, C, D, dt = cont2discrete(sys, Ts, method, alpha)
622+
return StateSpace(Ad, Bd, C, D, dt)
623+
624+
579625
# TODO: add discrete time check
580626
def _convertToStateSpace(sys, **kw):
581627
"""Convert a system to state space form (if needed).

‎control/tests/discrete_test.py

Lines changed: 40 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -259,29 +259,53 @@ def testSimulation(self):
259259
tout, yout, xout = forced_response(self.siso_ss2d, T, U, 0)
260260
tout, yout, xout = forced_response(self.siso_ss3d, T, U, 0)
261261

262-
@unittest.skip("skipping test_sample_system: not implemented for MIMO")
263262
def test_sample_system(self):
264263
# Make sure we can convert various types of systems
265-
for sysc in (self.siso_ss1, self.siso_ss1c, self.siso_tf1c):
266-
sysd = sample_system(sysc, 1, method='matched')
264+
for sysc in (self.siso_tf1, self.siso_tf1c,
265+
self.siso_ss1, self.siso_ss1c,
266+
self.mimo_ss1, self.mimo_ss1c):
267+
for method in ("zoh", "bilinear", "euler", "backward_diff"):
268+
sysd = sample_system(sysc, 1, method=method)
269+
self.assertEqual(sysd.dt, 1)
270+
271+
# Check "matched", defined only for SISO transfer functions
272+
for sysc in (self.siso_tf1, self.siso_tf1c):
273+
sysd = sample_system(sysc, 1, method="matched")
267274
self.assertEqual(sysd.dt, 1)
268275

269-
sysd = sample_system(sysc, 1, method='tustin')
270-
self.assertEqual(sysd.dt, 1)
271-
272-
sysd = sample_system(sysc, 1, method='zoh')
273-
self.assertEqual(sysd.dt, 1)
274-
# TODO: put in other generic checks
275-
276-
for sysc in (self.mimo_ss1, self.mimo_ss1c):
277-
sysd = sample_system(sysc, 1, method='zoh')
278-
self.assertEqual(sysd.dt, 1)
279-
280-
# TODO: check results of converstion
281-
282276
# Check errors
283277
self.assertRaises(ValueError, sample_system, self.siso_ss1d, 1)
284278
self.assertRaises(ValueError, sample_system, self.siso_ss1, 1, 'unknown')
279+
280+
def test_sample_ss(self):
281+
# double integrators, two different ways
282+
sys1 = StateSpace([[0.,1.],[0.,0.]], [[0.],[1.]], [[1.,0.]], 0.)
283+
sys2 = StateSpace([[0.,0.],[1.,0.]], [[1.],[0.]], [[0.,1.]], 0.)
284+
I = np.eye(2)
285+
for sys in (sys1, sys2):
286+
for h in (0.1, 0.5, 1, 2):
287+
Ad = I + h * sys.A
288+
Bd = h * sys.B + 0.5 * h**2 * (sys.A * sys.B)
289+
sysd = sample_system(sys, h, method='zoh')
290+
np.testing.assert_array_almost_equal(sysd.A, Ad)
291+
np.testing.assert_array_almost_equal(sysd.B, Bd)
292+
np.testing.assert_array_almost_equal(sysd.C, sys.C)
293+
np.testing.assert_array_almost_equal(sysd.D, sys.D)
294+
self.assertEqual(sysd.dt, h)
295+
296+
def test_sample_tf(self):
297+
# double integrator
298+
sys = TransferFunction(1, [1,0,0])
299+
for h in (0.1, 0.5, 1, 2):
300+
numd_expected = 0.5 * h**2 * np.array([1.,1.])
301+
dend_expected = np.array([1.,-2.,1.])
302+
sysd = sample_system(sys, h, method='zoh')
303+
self.assertEqual(sysd.dt, h)
304+
numd = sysd.num[0][0]
305+
dend = sysd.den[0][0]
306+
np.testing.assert_array_almost_equal(numd, numd_expected)
307+
np.testing.assert_array_almost_equal(dend, dend_expected)
308+
285309
def test_discrete_bode(self):
286310
# Create a simple discrete time system and check the calculation
287311
sys = TransferFunction([1], [1, 0.5], 1)

‎control/xferfcn.py

Lines changed: 75 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -76,7 +76,7 @@
7676
7777
Author: Richard M. Murray
7878
Date: 24 May 09
79-
Revised: Kevin K. Chewn, Dec 10
79+
Revised: Kevin K. Chen, Dec 10
8080
8181
$Id$
8282
@@ -86,7 +86,8 @@
8686
from numpy import angle, any, array, empty, finfo, insert, ndarray, ones, \
8787
polyadd, polymul, polyval, roots, sort, sqrt, zeros, squeeze, exp, pi, \
8888
where, delete, real, poly, poly1d
89-
from scipy.signal import lti
89+
import numpy as np
90+
from scipy.signal import lti, tf2zpk, zpk2tf, cont2discrete
9091
from copy import deepcopy
9192
from warnings import warn
9293
from .lti import Lti, timebaseEqual, timebase, isdtime
@@ -928,6 +929,78 @@ def _common_den(self, imag_tol=None):
928929

929930
return num, den
930931

932+
def sample(self, Ts, method='zoh', alpha=None):
933+
"""Convert a continuous-time system to discrete time
934+
935+
Creates a discrete-time system from a continuous-time system by
936+
sampling. Multiple methods of conversion are supported.
937+
938+
Parameters
939+
----------
940+
Ts : float
941+
Sampling period
942+
method : {"gbt", "bilinear", "euler", "backward_diff", "zoh", "matched"}
943+
Which method to use:
944+
945+
* gbt: generalized bilinear transformation
946+
* bilinear: Tustin's approximation ("gbt" with alpha=0.5)
947+
* euler: Euler (or forward differencing) method ("gbt" with alpha=0)
948+
* backward_diff: Backwards differencing ("gbt" with alpha=1.0)
949+
* zoh: zero-order hold (default)
950+
951+
alpha : float within [0, 1]
952+
The generalized bilinear transformation weighting parameter, which
953+
should only be specified with method="gbt", and is ignored otherwise
954+
955+
Returns
956+
-------
957+
sysd : StateSpace system
958+
Discrete time system, with sampling rate Ts
959+
960+
Notes
961+
-----
962+
1. Available only for SISO systems
963+
964+
2. Uses the command `cont2discrete` from `scipy.signal`
965+
966+
Examples
967+
--------
968+
>>> sys = TransferFunction(1, [1,1])
969+
>>> sysd = sys.sample(0.5, method='bilinear')
970+
971+
"""
972+
if not self.isctime():
973+
raise ValueError("System must be continuous time system")
974+
if not self.issiso():
975+
raise NotImplementedError("MIMO implementation not available")
976+
if method == "matched":
977+
return _c2dmatched(self, Ts)
978+
sys = (self.num[0][0], self.den[0][0])
979+
numd, dend, dt = cont2discrete(sys, Ts, method, alpha)
980+
return TransferFunction(numd[0,:], dend, dt)
981+
982+
# c2d function contributed by Benjamin White, Oct 2012
983+
def _c2dmatched(sysC, Ts):
984+
# Pole-zero match method of continuous to discrete time conversion
985+
szeros, spoles, sgain = tf2zpk(sysC.num[0][0], sysC.den[0][0])
986+
zzeros = [0] * len(szeros)
987+
zpoles = [0] * len(spoles)
988+
pregainnum = [0] * len(szeros)
989+
pregainden = [0] * len(spoles)
990+
for idx, s in enumerate(szeros):
991+
sTs = s*Ts
992+
z = exp(sTs)
993+
zzeros[idx] = z
994+
pregainnum[idx] = 1-z
995+
for idx, s in enumerate(spoles):
996+
sTs = s*Ts
997+
z = exp(sTs)
998+
zpoles[idx] = z
999+
pregainden[idx] = 1-z
1000+
zgain = np.multiply.reduce(pregainnum)/np.multiply.reduce(pregainden)
1001+
gain = sgain/zgain
1002+
sysDnum, sysDden = zpk2tf(zzeros, zpoles, gain)
1003+
return TransferFunction(sysDnum, sysDden, Ts)
9311004

9321005
# Utility function to convert a transfer function polynomial to a string
9331006
# Borrowed from poly1d library

0 commit comments

Comments
 (0)
0