Skip to content

Commit e1e250c

Browse files
murrayrmmurrayrm
committed
Initial implementation of sample_system (c2d) functionality based on code contributed by Benjamin White
 * Recreated dtime.py to hold new functions * New function sample_system to convert continuous time to discrete time * Use signal.cont2discrete for tustin and zohn * First cut at unit tests (need to check results, not just calling syntax) * Most functions are missing docstrings (will insert later)
1 parent e0e76be commit e1e250c

File tree

7 files changed

+164
-11
lines changed

7 files changed

+164
-11
lines changed

‎ChangeLog

Lines changed: 12 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,17 @@
11
2012-10-07 Richard Murray <murray@altura.local>
22

3+
* src/matlab.py (c2d): MATLAB compatible function (just calls
4+
sample_system)
5+
6+
* tests/discrete_test.py (TestDiscrete.test_sample_system): unit
7+
test for sampling
8+
9+
* src/dtime.py (sample_system, _c2dmatched): new functions, based on
10+
contributiosn from Benjamin White
11+
12+
* src/dtime.py: added back this file, this time with sampling routes
13+
to convert from continuous to discrete time
14+
315
* src/timeresp.py (forced_response): added discrete time simulator,
416
using dlsim from scipy.signal
517

‎src/__init__.py

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -59,6 +59,7 @@
5959
#! Should probably only import the exact functions we use...
6060
from bdalg import series, parallel, negate, feedback
6161
from delay import pade
62+
from dtime import sample_system
6263
from freqplot import bode_plot, nyquist_plot, gangof4_plot
6364
from freqplot import bode, nyquist, gangof4
6465
from lti import timebase, timebaseEqual, isdtime, isctime

‎src/dtime.py

Lines changed: 120 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,120 @@
1+
"""dtime.py
2+
3+
Functions for manipulating discrete time systems.
4+
5+
Routines in this module:
6+
7+
sample_system()
8+
_c2dmatched()
9+
"""
10+
11+
"""Copyright (c) 2012 by California Institute of Technology
12+
All rights reserved.
13+
14+
Redistribution and use in source and binary forms, with or without
15+
modification, are permitted provided that the following conditions
16+
are met:
17+
18+
1. Redistributions of source code must retain the above copyright
19+
notice, this list of conditions and the following disclaimer.
20+
21+
2. Redistributions in binary form must reproduce the above copyright
22+
notice, this list of conditions and the following disclaimer in the
23+
documentation and/or other materials provided with the distribution.
24+
25+
3. Neither the name of the California Institute of Technology nor
26+
the names of its contributors may be used to endorse or promote
27+
products derived from this software without specific prior
28+
written permission.
29+
30+
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
31+
"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
32+
LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
33+
FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL CALTECH
34+
OR THE CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
35+
SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
36+
LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
37+
USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
38+
ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
39+
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
40+
OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
41+
SUCH DAMAGE.
42+
43+
Author: Richard M. Murray
44+
Date: 6 October 2012
45+
46+
$Id: dtime.py 185 2012-08-30 05:44:32Z murrayrm $
47+
48+
"""
49+
50+
from scipy.signal import zpk2tf, tf2zpk, cont2discrete
51+
import numpy as np
52+
from cmath import exp
53+
from warnings import warn
54+
from lti import isctime
55+
from statesp import StateSpace, _convertToStateSpace
56+
from xferfcn import TransferFunction, _convertToTransferFunction
57+
58+
# Sample a continuous time system
59+
def sample_system(sysc, Ts, method='matched'):
60+
# TODO: add docstring
61+
62+
# Make sure we have a continuous time system
63+
if not isctime(sysc):
64+
raise ValueError("First argument must be continuous time system")
65+
66+
# TODO: impelement MIMO version
67+
if (sysc.inputs != 1 or sysc.outputs != 1):
68+
raise NotImplementedError("MIMO implementation not available")
69+
70+
# If we are passed a state space system, convert to transfer function first
71+
if isinstance(sysc, StateSpace):
72+
warn("sample_system: converting to transfer function")
73+
sysc = _convertToTransferFunction(sysc)
74+
75+
# Decide what to do based on the methods available
76+
if method == 'matched':
77+
sysd = _c2dmatched(sysc, Ts)
78+
79+
elif method == 'tustin':
80+
sys = [sysc.num[0][0], sysc.den[0][0]]
81+
scipySysD = cont2discrete(sys, Ts, method='bilinear')
82+
sysd = TransferFunction(scipySysD[0][0], scipySysD[1], dt)
83+
84+
elif method == 'zoh':
85+
sys = [sysc.num[0][0], sysc.den[0][0]]
86+
scipySysD = cont2discrete(sys, Ts, method='zoh')
87+
sysd = TransferFunction(scipySysD[0][0],scipySysD[1], dt)
88+
89+
elif method == 'foh' or method == 'impulse':
90+
raise ValueError("Method not developed yet")
91+
92+
else:
93+
raise ValueError("Invalid discretization method: %s" % method)
94+
95+
# TODO: Convert back into the input form
96+
# Set sampling time
97+
return sysd
98+
99+
# c2d function contributed by Benjamin White, Oct 2012
100+
def _c2dmatched(sysC, Ts):
101+
# Pole-zero match method of continuous to discrete time conversion
102+
szeros, spoles, sgain = tf2zpk(sysC.num[0][0], sysC.den[0][0])
103+
zzeros = [0] * len(szeros)
104+
zpoles = [0] * len(spoles)
105+
pregainnum = [0] * len(szeros)
106+
pregainden = [0] * len(spoles)
107+
for idx, s in enumerate(szeros):
108+
sTs = s*Ts
109+
z = exp(sTs)
110+
zzeros[idx] = z
111+
pregainnum[idx] = 1-z
112+
for idx, s in enumerate(spoles):
113+
sTs = s*Ts
114+
z = exp(sTs)
115+
zpoles[idx] = z
116+
pregainden[idx] = 1-z
117+
zgain = np.multiply.reduce(pregainnum)/np.multiply.reduce(pregainden)
118+
gain = sgain/zgain
119+
sysDnum, sysDden = zpk2tf(zzeros, zpoles, gain)
120+
return TransferFunction(sysDnum, sysDden, Ts)

‎src/matlab.py

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -77,6 +77,7 @@
7777
from statesp import StateSpace, _rss_generate, _convertToStateSpace
7878
from xferfcn import TransferFunction, _convertToTransferFunction
7979
from lti import Lti #base class of StateSpace, TransferFunction
80+
from dtime import sample_system
8081
from exception import ControlArgument
8182

8283
# Import MATLAB-like functions that can be used as-is
@@ -1386,3 +1387,9 @@ def tfdata(sys, **kw):
13861387
tf = _convertToTransferFunction(sys, **kw)
13871388

13881389
return (tf.num, tf.den)
1390+
1391+
# Convert a continuous time system to a discrete time system
1392+
def c2d(sysc, Ts, method):
1393+
# TODO: add docstring
1394+
# Call the sample_system() function to do the work
1395+
return sample_system(sysc, Ts, method)

‎src/timeresp.py

Lines changed: 10 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -356,19 +356,19 @@ def f_dot(x, _t):
356356
if len(U.shape) == 1:
357357
U = U.reshape(1,-1) #pylint: disable=E1103
358358

359-
# Create a callable that uses linear interpolation to
360-
# calculate the input at any time.
361-
compute_u = sp.interpolate.interp1d(T, U, kind='linear',
362-
copy=False,
363-
axis=-1, bounds_error=False,
364-
fill_value=0)
359+
# Create a callable that uses linear interpolation to
360+
# calculate the input at any time.
361+
compute_u = \
362+
sp.interpolate.interp1d(T, U, kind='linear', copy=False,
363+
axis=-1, bounds_error=False,
364+
fill_value=0)
365365

366366
# Function that computes the time derivative of the linear system
367-
def f_dot(x, t):
368-
return dot(A,x) + squeeze(dot(B,compute_u([t])))
367+
def f_dot(x, t):
368+
return dot(A,x) + squeeze(dot(B,compute_u([t])))
369369

370-
xout = sp.integrate.odeint(f_dot, X0, T, **keywords)
371-
yout = dot(C, xout.T) + dot(D, U)
370+
xout = sp.integrate.odeint(f_dot, X0, T, **keywords)
371+
yout = dot(C, xout.T) + dot(D, U)
372372

373373
yout = squeeze(yout)
374374
xout = xout.T

‎tests/discrete_test.py

Lines changed: 13 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -256,6 +256,19 @@ def testSimulation(self):
256256
tout, yout, xout = forced_response(self.siso_ss2d, T, U, 0)
257257
tout, yout, xout = forced_response(self.siso_ss3d, T, U, 0)
258258

259+
def test_sample_system(self):
260+
# Make sure we can convert various types of systems
261+
for sysc in (self.siso_ss1, self.siso_ss1c, self.siso_tf1c):
262+
sysd = sample_system(sysc, 1, method='matched')
263+
self.assertEqual(sysd.dt, 1)
264+
# TODO: put in other generic checks
265+
266+
# TODO: check results of converstion
267+
268+
# Check errors
269+
self.assertRaises(ValueError, sample_system, self.siso_ss1d, 1)
270+
self.assertRaises(ValueError, sample_system, self.siso_ss1, 1, 'unknown')
271+
259272
def suite():
260273
return unittest.TestLoader().loadTestsFromTestCase(TestDiscrete)
261274

‎tests/xferfcn_test.py

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -7,7 +7,7 @@
77
import numpy as np
88
from control.statesp import StateSpace
99
from control.xferfcn import TransferFunction, _convertToTransferFunction
10-
from control.dtime import isdtime
10+
from control.lti import isdtime
1111

1212
class TestXferFcn(unittest.TestCase):
1313
"""These are tests for functionality and correct reporting of the transfer

0 commit comments

Comments
 (0)