8000 Merge pull request #286 from murrayrm/timeresp-timevec · python-control/python-control@1254ccb · GitHub
[go: up one dir, main page]
Skip to content

Commit 1254ccb

Browse files
repaghrepagh
authored
Merge pull request #286 from murrayrm/timeresp-timevec
fix discrete-time simulation time output
2 parents 2ffb9c8 + 14261cc commit 1254ccb

File tree

2 files changed

+178
-7
lines changed

2 files changed

+178
-7
lines changed

‎control/tests/timeresp_test.py‎

Lines changed: 125 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -45,6 +45,14 @@ def setUp(self):
4545
"0. 9. ")
4646
self.mimo_ss1 = StateSpace(A, B, C, D)
4747

48+
# Create discrete time systems
49+
self.siso_dtf1 = TransferFunction([1], [1, 1, 0.25], True)
50+
self.siso_dtf2 = TransferFunction([1], [1, 1, 0.25], 0.2)
51+
self.siso_dss1 = tf2ss(self.siso_dtf1)
52+
self.siso_dss2 = tf2ss(self.siso_dtf2)
53+
self.mimo_dss1 = StateSpace(A, B, C, D, True)
54+
self.mimo_dss2 = c2d(self.mimo_ss1, 0.2)
55+
4856
def test_step_response(self):
4957
# Test SISO system
5058
sys = self.siso_ss1
@@ -320,6 +328,123 @@ def test_step_robustness(self):
320328
t2, y2 = step_response(sys2, input=0)
321329
np.testing.assert_array_almost_equal(y1, y2)
322330

331+
def test_time_vector(self):
332+
"Unit test: https://github.com/python-control/python-control/issues/239"
333+
# Discrete time simulations with specified time vectors
334+
Tin1 = np.arange(0, 5, 1) # matches dtf1, dss1; multiple of 0.2
335+
Tin2 = np.arange(0, 5, 0.2) # matches dtf2, dss2
336+
Tin3 = np.arange(0, 5, 0.5) # incompatible with 0.2
337+
338+
# Initial conditions to use for the different systems
339+
siso_x0 = [1, 2]
340+
mimo_x0 = [1, 2, 3, 4]
341+
342+
#
343+
# Easy cases: make sure that output sample time matches input
344+
#
345+
# No timebase in system => output should match input
346+
#
347+
# Initial response
348+
tout, yout = initial_response(self.siso_dtf1, Tin2, siso_x0)
349+
self.assertEqual(np.shape(tout), np.shape(yout[0,:]))
350+
np.testing.assert_array_equal(tout, Tin2)
351+
352+
# Impulse response
353+
tout, yout = impulse_response(self.siso_dtf1, Tin2)
354+
self.assertEqual(np.shape(tout), np.shape(yout[0,:]))
355+
np.testing.assert_array_equal(tout, Tin2)
356+
357+
# Step response
358+
tout, yout = step_response(self.siso_dtf1, Tin2)
359+
self.assertEqual(np.shape(tout), np.shape(yout[0,:]))
360+
np.testing.assert_array_equal(tout, Tin2)
361+
362+
# Forced response with specified time vector
363+
tout, yout, xout = forced_response(self.siso_dtf1, Tin2, np.sin(Tin2))
364+
self.assertEqual(np.shape(tout), np.shape(yout[0,:]))
365+
np.testing.assert_array_equal(tout, Tin2)
366+
367+
# Forced response with no time vector, no sample time (should use 1)
368+
tout, yout, xout = forced_response(self.siso_dtf1, None, np.sin(Tin1))
369+
self.assertEqual(np.shape(tout), np.shape(yout[0,:]))
370+
np.testing.assert_array_equal(tout, Tin1)
371+
372+
# MIMO forced response
373+
tout, yout, xout = forced_response(self.mimo_dss1, Tin1,
374+
(np.sin(Tin1), np.cos(Tin1)),
375+
mimo_x0)
376+
self.assertEqual(np.shape(tout), np.shape(yout[0,:]))
377+
self.assertEqual(np.shape(tout), np.shape(yout[1,:]))
378+
np.testing.assert_array_equal(tout, Tin1)
379+
380+
# Matching timebase in system => output should match input
381+
#
382+
# Initial response
383+
tout, yout = initial_response(self.siso_dtf2, Tin2, siso_x0)
384+
self.assertEqual(np.shape(tout), np.shape(yout[0,:]))
385+
np.testing.assert_array_equal(tout, Tin2)
386+
387+
# Impulse response
388+
tout, yout = impulse_response(self.siso_dtf2, Tin2)
389+
self.assertEqual(np.shape(tout), np.shape(yout[0,:]))
390+
np.testing.assert_array_equal(tout, Tin2)
391+
392+
# Step response
393+
tout, yout = step_response(self.siso_dtf2, Tin2)
394+
self.assertEqual(np.shape(tout), np.shape(yout[0,:]))
395+
np.testing.assert_array_equal(tout, Tin2)
396+
397+
# Forced response
398+
tout, yout, xout = forced_response(self.siso_dtf2, Tin2, np.sin(Tin2))
399+
self.assertEqual(np.shape(tout), np.shape(yout[0,:]))
400+
np.testing.assert_array_equal(tout, Tin2)
401+
402+
# Forced response with no time vector, use sample time
403+
tout, yout, xout = forced_response(self.siso_dtf2, None, np.sin(Tin2))
404+
self.assertEqual(np.shape(tout), np.shape(yout[0,:]))
405+
np.testing.assert_array_equal(tout, Tin2)
406+
407+
# Compatible timebase in system => output should match input
408+
#
409+
# Initial response
410+
tout, yout = initial_response(self.siso_dtf2, Tin1, siso_x0)
411+
self.assertEqual(np.shape(tout), np.shape(yout[0,:]))
412+
np.testing.assert_array_equal(tout, Tin1)
413+
414+
# Impulse response
415+
tout, yout = impulse_response(self.siso_dtf2, Tin1)
416+
self.assertEqual(np.shape(tout), np.shape(yout[0,:]))
417+
np.testing.assert_array_equal(tout, Tin1)
418+
419+
# Step response
420+
tout, yout = step_response(self.siso_dtf2, Tin1)
421+
self.assertEqual(np.shape(tout), np.shape(yout[0,:]))
422+
np.testing.assert_array_equal(tout, Tin1)
423+
424+
# Forced response
425+
tout, yout, xout = forced_response(self.siso_dtf2, Tin1, np.sin(Tin1))
426+
self.assertEqual(np.shape(tout), np.shape(yout[0,:]))
427+
np.testing.assert_array_equal(tout, Tin1)
428+
429+
#
430+
# Interpolation of the input (to match scipy.signal.dlsim)
431+
#
432+
# Initial response
433+
tout, yout, xout = forced_response(self.siso_dtf2, Tin1,
434+
np.sin(Tin1), interpolate=True)
435+
self.assertEqual(np.shape(tout), np.shape(yout[0,:]))
436+
self.assertTrue(np.allclose(tout[1:] - tout[:-1], self.siso_dtf2.dt))
437+
438+
#
439+
# Incompatible cases: make sure an error is thrown
440+
#
441+
# System timebase and given time vector are incompatible
442+
#
443+
# Initial response
444+
with self.assertRaises(Exception) as context:
445+
tout, yout = initial_response(self.siso_dtf2, Tin3, siso_x0)
446+
self.assertTrue(isinstance(context.exception, ValueError))
447+
323448
def suite():
324449
return unittest.TestLoader().loadTestsFromTestCase(TestTimeresp)
325450

‎control/timeresp.py‎

Lines changed: 53 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -169,7 +169,8 @@ def shape_matches(s_legal, s_actual):
169169

170170

171171
# Forced response of a linear system
172-
def forced_response(sys, T=None, U=0., X0=0., transpose=False):
172+
def forced_response(sys, T=None, U=0., X0=0., transpose=False,
173+
interpolate=False):
173174
"""Simulate the output of a linear system.
174175
175176
As a convenience for parameters `U`, `X0`:
@@ -200,6 +201,13 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False):
200201
If True, transpose all input and output arrays (for backward
201202
compatibility with MATLAB and scipy.signal.lsim)
202203
204+
interpolate:bool
205+
If True and system is a discrete time system, the input will
206+
be interpolated between the given time steps and the output
207+
will be given at system sampling rate. Otherwise, only return
208+
the output at the times given in `T`. No effect on continuous
209+
time simulations (default = False).
210+
203211
Returns
204212
-------
205213
T: array
@@ -236,8 +244,8 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False):
236244
if U is None:
237245
raise ValueError('Parameters ``T`` and ``U`` can\'t both be'
238246
'zero for discrete-time simulation')
239-
# Set T to integers with same length as U
240-
T = range(len(U))
247+
# Set T to equally spaced samples with same length as U
248+
T = np.array(range(len(U))) * (1 if sys.dt == True else sys.dt)
241249
else:
242250
# Make sure the input vector and time vector have same length
243251
# TODO: allow interpolation of the input vector
@@ -316,26 +324,53 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False):
316324
dot(Bd1, U[:, i]))
317325
yout = dot(C, xout) + dot(D, U)
318326

327+
tout = T
319328
yout = squeeze(yout)
320329
xout = squeeze(xout)
321330

322331
else:
332+
# Discrete type system => use SciPy signal processing toolbox
333+
if (sys.dt != True):
334+
# Make sure that the time increment is a multiple of sampling time
335+
336+
# First make sure that time increment is bigger than sampling time
337+
if dt < sys.dt:
338+
raise ValueError("Time steps ``T`` must match sampling time")
339+
340+
# Now check to make sure it is a multiple (with check against
341+
# sys.dt because floating point mod can have small errors
342+
elif not (np.isclose(dt % sys.dt, 0) or
343+
np.isclose(dt % sys.dt, sys.dt)):
344+
raise ValueError("Time steps ``T`` must be multiples of " \
345+
"sampling time")
346+
else:
347+
sys.dt = dt # For unspecified sampling time, use time incr
348+
323349
# Discrete time simulation using signal processing toolbox
324350
dsys = (A, B, C, D, sys.dt)
351+
352+
# Use signal processing toolbox for the discrete time simulation
325353
# Transpose the input to match toolbox convention
326354
tout, yout, xout = sp.signal.dlsim(dsys, np.transpose(U), T, X0)
327355

356+
if not interpolate:
357+
# If dt is different from sys.dt, resample the output
358+
inc = int(round(dt / sys.dt))
359+
tout = T # Return exact list of time steps
360+
yout = yout[::inc,:]
361+
xout = xout[::inc,:]
362+
328363
# Transpose the output and state vectors to match local convention
329364
xout = sp.transpose(xout)
330365
yout = sp.transpose(yout)
331366

332367
# See if we need to transpose the data back into MATLAB form
333368
if (transpose):
334-
T = np.transpose(T)
369+
tout = np.transpose(tout)
335370
yout = np.transpose(yout)
336371
xout = np.transpose(xout)
337372

338-
return T, yout, xout
373+
return tout, yout, xout
339374

340375
def _get_ss_simo(sys, input=None, output=None):
341376
"""Return a SISO or SIMO state-space version of sys
@@ -588,7 +623,12 @@ def initial_response(sys, T=None, X0=0., input=0, output=None,
588623
# Create time and input vectors; checking is done in forced_response(...)
589624
# The initial vector X0 is created in forced_response(...) if necessary
590625
if T is None:
591-
T = _default_response_times(sys.A, 100)
626+
if isctime(sys):
627+
T = _default_response_times(sys.A, 1000)
628+
else:
629+
# For discrete time, use integers
630+
tvec = _default_response_times(sys.A, 1000)
631 A593 +
T = range(int(np.ceil(max(tvec))))
592632
U = np.zeros_like(T)
593633

594634
T, yout, _xout = forced_response(sys, T, U, X0, transpose=transpose)
@@ -672,7 +712,13 @@ def impulse_response(sys, T=None, X0=0., input=0, output=None,
672712

673713
# Compute T and U, no checks necessary, they will be checked in lsim
674714
if T is None:
675-
T = _default_response_times(sys.A, 100)
715+
if isctime(sys):
716+
T = _default_response_times(sys.A, 100)
717+
else:
718+
# For discrete time, use integers
719+
tvec = _default_response_times(sys.A, 100)
720+
T = range(int(np.ceil(max(tvec))))
721+
676722
U = np.zeros_like(T)
677723

678724
# Compute new X0 that contains the impulse

0 commit comments

Comments
 (0)
0