Repository navigation
Expand file tree
/
Copy pathtimeresp.py
More file actions
1348 lines (1098 loc) · 52.2 KB
/
Copy pathtimeresp.py
File metadata and controls
1348 lines (1098 loc) · 52.2 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
# timeresp.py - time-domain simulation routines
#
# This file contains a collection of functions that calculate time
# responses for linear systems.
"""The :mod:`~control.timeresp` module contains a collection of
functions that are used to compute time-domain simulations of LTI
systems.
Arguments to time-domain simulations include a time vector, an input
vector (when needed), and an initial condition vector. The most
general function for simulating LTI systems the
:func:`forced_response` function, which has the form::
t, y = forced_response(sys, T, U, X0)
where `T` is a vector of times at which the response should be
evaluated, `U` is a vector of inputs (one for each time point) and
`X0` is the initial condition for the system.
See :ref:`time-series-convention` for more information on how time
series data are represented.
"""
"""Copyright (c) 2011 by California Institute of Technology
All rights reserved.
Copyright (c) 2011 by Eike Welk
Copyright (c) 2010 by SciPy Developers
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions
are met:
1. Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.
2. Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in the
documentation and/or other materials provided with the distribution.
3. Neither the name of the California Institute of Technology nor
the names of its contributors may be used to endorse or promote
products derived from this software without specific prior
written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL CALTECH
OR THE CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
SUCH DAMAGE.
Initial Author: Eike Welk
Date: 12 May 2011
Modified: Sawyer B. Fuller (minster@uw.edu) to add discrete-time
capability and better automatic time vector creation
Date: June 2020
Modified by Ilhan Polat to improve automatic time vector creation
Date: August 17, 2020
$Id$
"""
# Libraries that we make use of
import scipy as sp # SciPy library (used all over)
import numpy as np # NumPy library
from scipy.linalg import eig, eigvals, matrix_balance, norm
from numpy import (einsum, maximum, minimum,
atleast_1d)
import warnings
from .lti import LTI # base class of StateSpace, TransferFunction
from .xferfcn import TransferFunction
from .statesp import _convert_to_statespace, _mimo2simo, _mimo2siso, ssdata
from .lti import isdtime, isctime
from . import config
__all__ = ['forced_response', 'step_response', 'step_info', 'initial_response',
'impulse_response']
# Helper function for checking array-like parameters
def _check_convert_array(in_obj, legal_shapes, err_msg_start, squeeze=False,
transpose=False):
"""Helper function for checking array_like parameters.
* Check type and shape of ``in_obj``.
* Convert ``in_obj`` to an array if necessary.
* Change shape of ``in_obj`` according to parameter ``squeeze``.
* If ``in_obj`` is a scalar (number) it is converted to an array with
a legal shape, that is filled with the scalar value.
The function raises an exception when it detects an error.
Parameters
----------
in_obj : array like object
The array or matrix which is checked.
legal_shapes : list of tuple
A list of shapes that in_obj can legally have.
The special value "any" means that there can be any
number of elements in a certain dimension.
* ``(2, 3)`` describes an array with 2 rows and 3 columns
* ``(2, "any")`` describes an array with 2 rows and any number of
columns
err_msg_start : str
String that is prepended to the error messages, when this function
raises an exception. It should be used to identify the argument which
is currently checked.
squeeze : bool
If True, all dimensions with only one element are removed from the
array. If False the array's shape is unmodified.
For example:
``array([[1,2,3]])`` is converted to ``array([1, 2, 3])``
transpose : bool, optional
If True, assume that 2D input arrays are transposed from the standard
format. Used to convert MATLAB-style inputs to our format.
Returns
-------
out_array : array
The checked and converted contents of ``in_obj``.
"""
# convert nearly everything to an array.
out_array = np.asarray(in_obj)
if (transpose):
out_array = np.transpose(out_array)
# Test element data type, elements must be numbers
legal_kinds = set(("i", "f", "c")) # integer, float, complex
if out_array.dtype.kind not in legal_kinds:
err_msg = "Wrong element data type: '{d}'. Array elements " \
"must be numbers.".format(d=str(out_array.dtype))
raise TypeError(err_msg_start + err_msg)
# If array is zero dimensional (in_obj is scalar):
# create array with legal shape filled with the original value.
if out_array.ndim == 0:
for s_legal in legal_shapes:
# search for shape that does not contain the special symbol any.
if "any" in s_legal:
continue
the_val = out_array[()]
out_array = np.empty(s_legal, 'd')
out_array.fill(the_val)
break
# Test shape
def shape_matches(s_legal, s_actual):
"""Test if two shape tuples match"""
# Array must have required number of dimensions
if len(s_legal) != len(s_actual):
return False
# All dimensions must contain required number of elements. Joker: "all"
for n_legal, n_actual in zip(s_legal, s_actual):
if n_legal == "any":
continue
if n_legal != n_actual:
return False
return True
# Iterate over legal shapes, and see if any matches out_array's shape.
for s_legal in legal_shapes:
if shape_matches(s_legal, out_array.shape):
break
else:
legal_shape_str = " or ".join([str(s) for s in legal_shapes])
err_msg = "Wrong shape (rows, columns): {a}. Expected: {e}." \
.format(e=legal_shape_str, a=str(out_array.shape))
raise ValueError(err_msg_start + err_msg)
# Convert shape
if squeeze:
out_array = np.squeeze(out_array)
# We don't want zero dimensional arrays
if out_array.shape == tuple():
out_array = out_array.reshape((1,))
return out_array
# Forced response of a linear system
def forced_response(sys, T=None, U=0., X0=0., transpose=False,
interpolate=False, return_x=None, squeeze=None):
"""Simulate the output of a linear system.
As a convenience for parameters `U`, `X0`:
Numbers (scalars) are converted to constant arrays with the correct shape.
The correct shape is inferred from arguments `sys` and `T`.
For information on the **shape** of parameters `U`, `T`, `X0` and
return values `T`, `yout`, `xout`, see :ref:`time-series-convention`.
Parameters
----------
sys : LTI (StateSpace or TransferFunction)
LTI system to simulate
T : array_like, optional for discrete LTI `sys`
Time steps at which the input is defined; values must be evenly spaced.
U : array_like or float, optional
Input array giving input at each time `T` (default = 0).
If `U` is ``None`` or ``0``, a special algorithm is used. This special
algorithm is faster than the general algorithm, which is used
otherwise.
X0 : array_like or float, optional
Initial condition (default = 0).
transpose : bool, optional
If True, transpose all input and output arrays (for backward
compatibility with MATLAB and :func:`scipy.signal.lsim`). Default
value is False.
interpolate : bool, optional (default=False)
If True and system is a discrete time system, the input will
be interpolated between the given time steps and the output
will be given at system sampling rate. Otherwise, only return
the output at the times given in `T`. No effect on continuous
time simulations (default = False).
return_x : bool, optional
If True (default), return the the state vector. Set to False to
return only the time and output vectors.
squeeze : bool, optional
By default, if a system is single-input, single-output (SISO) then
the output response is returned as a 1D array (indexed by time). If
squeeze=True, remove single-dimensional entries from the shape of
the output even if the system is not SISO. If squeeze=False, keep
the output as a 2D array (indexed by the output number and time)
even if the system is SISO. The default value can be set using
config.defaults['control.squeeze_time_response'].
Returns
-------
T : array
Time values of the output.
yout : array
Response of the system. If the system is SISO and squeeze is not
True, the array is 1D (indexed by time). If the system is not SISO or
squeeze is False, the array is 2D (indexed by the output number and
time).
xout : array
Time evolution of the state vector. Not affected by squeeze.
See Also
--------
step_response, initial_response, impulse_response
Notes
-----
For discrete time systems, the input/output response is computed using the
:func:`scipy.signal.dlsim` function.
For continuous time systems, the output is computed using the matrix
exponential `exp(A t)` and assuming linear interpolation of the inputs
between time points.
Examples
--------
>>> T, yout, xout = forced_response(sys, T, u, X0)
See :ref:`time-series-convention`.
"""
if not isinstance(sys, LTI):
raise TypeError('Parameter ``sys``: must be a ``LTI`` object. '
'(For example ``StateSpace`` or ``TransferFunction``)')
# If return_x was not specified, figure out the default
if return_x is None:
return_x = config.defaults['forced_response.return_x']
# If return_x is used for TransferFunction, issue a warning
if return_x and isinstance(sys, TransferFunction):
warnings.warn(
"return_x specified for a transfer function system. Internal "
"conversion to state space used; results may meaningless.")
sys = _convert_to_statespace(sys)
A, B, C, D = np.asarray(sys.A), np.asarray(sys.B), np.asarray(sys.C), \
np.asarray(sys.D)
# d_type = A.dtype
n_states = A.shape[0]
n_inputs = B.shape[1]
n_outputs = C.shape[0]
# Convert inputs to numpy arrays for easier shape checking
if U is not None:
U = np.asarray(U)
if T is not None:
T = np.asarray(T)
# Set and/or check time vector in discrete time case
if isdtime(sys, strict=True):
if T is None:
if U is None:
raise ValueError('Parameters ``T`` and ``U`` can\'t both be'
'zero for discrete-time simulation')
# Set T to equally spaced samples with same length as U
if U.ndim == 1:
n_steps = U.shape[0]
else:
n_steps = U.shape[1]
T = np.array(range(n_steps)) * (1 if sys.dt is True else sys.dt)
else:
# Make sure the input vector and time vector have same length
# TODO: allow interpolation of the input vector
if (U.ndim == 1 and U.shape[0] != T.shape[0]) or \
(U.ndim > 1 and U.shape[1] != T.shape[0]):
ValueError('Pamameter ``T`` must have same elements as'
' the number of columns in input array ``U``')
# Test if T has shape (n,) or (1, n);
# T must be array-like and values must be increasing.
# The length of T determines the length of the input vector.
if T is None:
raise ValueError('Parameter ``T``: must be array-like, and contain '
'(strictly monotonic) increasing numbers.')
T = _check_convert_array(T, [('any',), (1, 'any')],
'Parameter ``T``: ', squeeze=True,
transpose=transpose)
dt = T[1] - T[0]
if not np.allclose(T[1:] - T[:-1], dt):
raise ValueError("Parameter ``T``: time values must be "
"equally spaced.")
n_steps = T.shape[0] # number of simulation steps
# create X0 if not given, test if X0 has correct shape
X0 = _check_convert_array(X0, [(n_states,), (n_states, 1)],
'Parameter ``X0``: ', squeeze=True)
# If we are passed a transfer function and X0 is non-zero, warn the user
if isinstance(sys, TransferFunction) and np.any(X0 != 0):
warnings.warn(
"Non-zero initial condition given for transfer function system. "
"Internal conversion to state space used; may not be consistent "
"with given X0.")
xout = np.zeros((n_states, n_steps))
xout[:, 0] = X0
yout = np.zeros((n_outputs, n_steps))
# Separate out the discrete and continuous time cases
if isctime(sys):
# Solve the differential equation, copied from scipy.signal.ltisys.
dot = np.dot # Faster and shorter code
# Faster algorithm if U is zero
if U is None or (isinstance(U, (int, float)) and U == 0):
# Solve using matrix exponential
expAdt = sp.linalg.expm(A * dt)
for i in range(1, n_steps):
xout[:, i] = dot(expAdt, xout[:, i-1])
yout = dot(C, xout)
# General algorithm that interpolates U in between output points
else:
# Test if U has correct shape and type
legal_shapes = [(n_steps,), (1, n_steps)] if n_inputs == 1 else \
[(n_inputs, n_steps)]
U = _check_convert_array(U, legal_shapes,
'Parameter ``U``: ', squeeze=False,
transpose=transpose)
# convert 1D array to 2D array with only one row
if len(U.shape) == 1:
U = U.reshape(1, -1) # pylint: disable=E1103
# Algorithm: to integrate from time 0 to time dt, with linear
# interpolation between inputs u(0) = u0 and u(dt) = u1, we solve
# xdot = A x + B u, x(0) = x0
# udot = (u1 - u0) / dt, u(0) = u0.
#
# Solution is
# [ x(dt) ] [ A*dt B*dt 0 ] [ x0 ]
# [ u(dt) ] = exp [ 0 0 I ] [ u0 ]
# [u1 - u0] [ 0 0 0 ] [u1 - u0]
M = np.block([[A * dt, B * dt, np.zeros((n_states, n_inputs))],
[np.zeros((n_inputs, n_states + n_inputs)),
np.identity(n_inputs)],
[np.zeros((n_inputs, n_states + 2 * n_inputs))]])
expM = sp.linalg.expm(M)
Ad = expM[:n_states, :n_states]
Bd1 = expM[:n_states, n_states+n_inputs:]
Bd0 = expM[:n_states, n_states:n_states + n_inputs] - Bd1
for i in range(1, n_steps):
xout[:, i] = (dot(Ad, xout[:, i-1]) + dot(Bd0, U[:, i-1]) +
dot(Bd1, U[:, i]))
yout = dot(C, xout) + dot(D, U)
tout = T
else:
# Discrete type system => use SciPy signal processing toolbox
if sys.dt is not True:
# Make sure that the time increment is a multiple of sampling time
# First make sure that time increment is bigger than sampling time
# (with allowance for small precision errors)
if dt < sys.dt and not np.isclose(dt, sys.dt):
raise ValueError("Time steps ``T`` must match sampling time")
# Now check to make sure it is a multiple (with check against
# sys.dt because floating point mod can have small errors
elif not (np.isclose(dt % sys.dt, 0) or
np.isclose(dt % sys.dt, sys.dt)):
raise ValueError("Time steps ``T`` must be multiples of "
"sampling time")
sys_dt = sys.dt
else:
sys_dt = dt # For unspecified sampling time, use time incr
# Discrete time simulation using signal processing toolbox
dsys = (A, B, C, D, sys_dt)
# Use signal processing toolbox for the discrete time simulation
# Transpose the input to match toolbox convention
tout, yout, xout = sp.signal.dlsim(dsys, np.transpose(U), T, X0)
if not interpolate:
# If dt is different from sys.dt, resample the output
inc = int(round(dt / sys_dt))
tout = T # Return exact list of time steps
yout = yout[::inc, :]
xout = xout[::inc, :]
# Transpose the output and state vectors to match local convention
xout = np.transpose(xout)
yout = np.transpose(yout)
return _process_time_response(sys, tout, yout, xout, transpose=transpose,
return_x=return_x, squeeze=squeeze)
# Process time responses in a uniform way
def _process_time_response(
sys, tout, yout, xout, transpose=None, return_x=False,
squeeze=None, input=None, output=None):
"""Process time response signals.
This function processes the outputs of the time response functions and
processes the transpose and squeeze keywords.
Parameters
----------
sys : LTI or InputOutputSystem
System that generated the data (used to check if SISO/MIMO).
T : 1D array
Time values of the output. Ignored if None.
yout : ndarray
Response of the system. This can either be a 1D array indexed by time
(for SISO systems), a 2D array indexed by output and time (for MIMO
systems with no input indexing, such as initial_response or forced
response) or a 3D array indexed by output, input, and time.
xout : array, optional
Individual response of each x variable (if return_x is True). For a
SISO system (or if a single input is specified), this should be a 2D
array indexed by the state index and time (for single input systems)
or a 3D array indexed by state, input, and time. Ignored if None.
transpose : bool, optional
If True, transpose all input and output arrays (for backward
compatibility with MATLAB and :func:`scipy.signal.lsim`). Default
value is False.
return_x : bool, optional
If True, return the state vector (default = False).
squeeze : bool, optional
By default, if a system is single-input, single-output (SISO) then the
output response is returned as a 1D array (indexed by time). If
squeeze=True, remove single-dimensional entries from the shape of the
output even if the system is not SISO. If squeeze=False, keep the
output as a 3D array (indexed by the output, input, and time) even if
the system is SISO. The default value can be set using
config.defaults['control.squeeze_time_response'].
input : int, optional
If present, the response represents only the listed input.
output : int, optional
If present, the response represents only the listed output.
Returns
-------
T : 1D array
Time values of the output
yout : ndarray
Response of the system. If the system is SISO and squeeze is not
True, the array is 1D (indexed by time). If the system is not SISO or
squeeze is False, the array is either 2D (indexed by output and time)
or 3D (indexed by input, output, and time).
xout : array, optional
Individual response of each x variable (if return_x is True). For a
SISO system (or if a single input is specified), xout is a 2D array
indexed by the state index and time. For a non-SISO system, xout is a
3D array indexed by the state, the input, and time. The shape of xout
is not affected by the ``squeeze`` keyword.
"""
# If squeeze was not specified, figure out the default (might remain None)
if squeeze is None:
squeeze = config.defaults['control.squeeze_time_response']
# Determine if the system is SISO
issiso = sys.issiso() or (input is not None and output is not None)
# Figure out whether and how to squeeze output data
if squeeze is True: # squeeze all dimensions
yout = np.squeeze(yout)
elif squeeze is False: # squeeze no dimensions
pass
elif squeeze is None: # squeeze signals if SISO
if issiso:
if len(yout.shape) == 3:
yout = yout[0][0] # remove input and output
else:
yout = yout[0] # remove input
else:
raise ValueError("unknown squeeze value")
# Figure out whether and how to squeeze the state data
if issiso and xout is not None and len(xout.shape) > 2:
xout = xout[:, 0, :] # remove input
# See if we need to transpose the data back into MATLAB form
if transpose:
# Transpose time vector in case we are using np.matrix
tout = np.transpose(tout)
# For signals, put the last index (time) into the first slot
yout = np.transpose(yout, np.roll(range(yout.ndim), 1))
if xout is not None:
xout = np.transpose(xout, np.roll(range(xout.ndim), 1))
# Return time, output, and (optionally) state
return (tout, yout, xout) if return_x else (tout, yout)
def _get_ss_simo(sys, input=None, output=None, squeeze=None):
"""Return a SISO or SIMO state-space version of sys.
This function converts the given system to a state space system in
preparation for simulation and sets the system matrixes to match the
desired input and output.
If input is not specified, select first input and issue warning (legacy
behavior that should eventually not be used).
If the output is not specified, report on all outputs.
"""
# If squeeze was not specified, figure out the default
if squeeze is None:
squeeze = config.defaults['control.squeeze_time_response']
sys_ss = _convert_to_statespace(sys)
if sys_ss.issiso():
return squeeze, sys_ss
elif squeeze == None and (input is None or output is None):
# Don't squeeze outputs if resulting system turns out to be siso
# Note: if we expand input to allow a tuple, need to update this check
squeeze = False
warn = False
if input is None:
# issue warning if input is not given
warn = True
input = 0
if output is None:
return squeeze, _mimo2simo(sys_ss, input, warn_conversion=warn)
else:
return squeeze, _mimo2siso(sys_ss, input, output, warn_conversion=warn)
def step_response(sys, T=None, X0=0., input=None, output=None, T_num=None,
transpose=False, return_x=False, squeeze=None):
# pylint: disable=W0622
"""Compute the step response for a linear system.
If the system has multiple inputs and/or multiple outputs, the step
response is computed for each input/output pair, with all other inputs set
to zero. Optionally, a single input and/or single output can be selected,
in which case all other inputs are set to 0 and all other outputs are
ignored.
For information on the **shape** of parameters `T`, `X0` and
return values `T`, `yout`, see :ref:`time-series-convention`.
Parameters
----------
sys : StateSpace or TransferFunction
LTI system to simulate
T : array_like or float, optional
Time vector, or simulation time duration if a number. If T is not
provided, an attempt is made to create it automatically from the
dynamics of sys. If sys is continuous-time, the time increment dt
is chosen small enough to show the fastest mode, and the simulation
time period tfinal long enough to show the slowest mode, excluding
poles at the origin and pole-zero cancellations. If this results in
too many time steps (>5000), dt is reduced. If sys is discrete-time,
only tfinal is computed, and final is reduced if it requires too
many simulation steps.
X0 : array_like or float, optional
Initial condition (default = 0). Numbers are converted to constant
arrays with the correct shape.
input : int, optional
Only compute the step response for the listed input. If not
specified, the step responses for each independent input are computed.
output : int, optional
Only report the step response for the listed output. If not
specified, all outputs are reported.
T_num : int, optional
Number of time steps to use in simulation if T is not provided as an
array (autocomputed if not given); ignored if sys is discrete-time.
transpose : bool, optional
If True, transpose all input and output arrays (for backward
compatibility with MATLAB and :func:`scipy.signal.lsim`). Default
value is False.
return_x : bool, optional
If True, return the state vector (default = False).
squeeze : bool, optional
By default, if a system is single-input, single-output (SISO) then the
output response is returned as a 1D array (indexed by time). If
squeeze=True, remove single-dimensional entries from the shape of the
output even if the system is not SISO. If squeeze=False, keep the
output as a 3D array (indexed by the output, input, and time) even if
the system is SISO. The default value can be set using
config.defaults['control.squeeze_time_response'].
Returns
-------
T : 1D array
Time values of the output
yout : ndarray
Response of the system. If the system is SISO and squeeze is not
True, the array is 1D (indexed by time). If the system is not SISO or
squeeze is False, the array is 3D (indexed by the input, output, and
time).
xout : array, optional
Individual response of each x variable (if return_x is True). For a
SISO system (or if a single input is specified), xout is a 2D array
indexed by the state index and time. For a non-SISO system, xout is a
3D array indexed by the state, the input, and time. The shape of xout
is not affected by the ``squeeze`` keyword.
See Also
--------
forced_response, initial_response, impulse_response
Notes
-----
This function uses the `forced_response` function with the input set to a
unit step.
Examples
--------
>>> T, yout = step_response(sys, T, X0)
"""
# Create the time and input vectors
if T is None or np.asarray(T).size == 1:
T = _default_time_vector(sys, N=T_num, tfinal=T, is_step=True)
U = np.ones_like(T)
# If we are passed a transfer function and X0 is non-zero, warn the user
if isinstance(sys, TransferFunction) and np.any(X0 != 0):
warnings.warn(
"Non-zero initial condition given for transfer function system. "
"Internal conversion to state space used; may not be consistent "
"with given X0.")
# Convert to state space so that we can simulate
sys = _convert_to_statespace(sys)
# Set up arrays to handle the output
ninputs = sys.ninputs if input is None else 1
noutputs = sys.noutputs if output is None else 1
yout = np.empty((noutputs, ninputs, np.asarray(T).size))
xout = np.empty((sys.nstates, ninputs, np.asarray(T).size))
# Simulate the response for each input
for i in range(sys.ninputs):
# If input keyword was specified, only simulate for that input
if isinstance(input, int) and i != input:
continue
# Create a set of single inputs system for simulation
squeeze, simo = _get_ss_simo(sys, i, output, squeeze=squeeze)
out = forced_response(simo, T, U, X0, transpose=False,
return_x=return_x, squeeze=True)
inpidx = i if input is None else 0
yout[:, inpidx, :] = out[1]
if return_x:
xout[:, i, :] = out[2]
return _process_time_response(
sys, out[0], yout, xout, transpose=transpose, return_x=return_x,
squeeze=squeeze, input=input, output=output)
def step_info(sys, T=None, T_num=None, SettlingTimeThreshold=0.02,
RiseTimeLimits=(0.1, 0.9)):
'''
Step response characteristics (Rise time, Settling Time, Peak and others).
Parameters
----------
sys : SISO or MIMO dynamic system model. Dynamic systems that you can use include:
StateSpace or TransferFunction
LTI system to simulate
T : array_like or float, optional
Time vector, or simulation time duration if a number (time vector is
autocomputed if not given, see :func:`step_response` for more detail)
T_num : int, optional
Number of time steps to use in simulation if T is not provided as an
array (autocomputed if not given); ignored if sys is discrete-time.
SettlingTimeThreshold : float value, optional
Defines the error to compute settling time (default = 0.02)
RiseTimeLimits : tuple (lower_threshold, upper_theshold)
Defines the lower and upper threshold for RiseTime computation
Returns
-------
S: a dictionary for SISO systems or a list of Ny-by-Nu of dictionaries
for MIMO systems, each dictionary containing:
RiseTime: Rise time (by default time from 10% to 90% of the steady-state value).
SettlingTime: Time to enter inside a default error of 2%
SettlingMin: Minimum value of `y(t)` after RiseTime
SettlingMax: Maximum value of `y(t)` after RiseTime
Overshoot: Percentage of the Peak relative to steady value
Undershoot: Percentage of undershoot
Peak: Peak absolute value of `y(t)`
PeakTime: Time of the peak value
SteadyStateValue: Steady-state value
See Also
--------
step_response, initial_response, impulse_response, forced_response
Examples
--------
>>> from control import tf, step_info
>>> sys = tf([3.32, 0, -162.8],[1, 24.56, 186.5, 457.8, 116.2])
>>> step_info(sys)
>>>
>>> {'RiseTime': 7.70245002940221,
>>> 'SettlingTime': 14.13151114265325,
>>> 'SettlingMin': -1.3994264225116284,
>>> 'SettlingMax': -1.2612640212267976,
>>> 'Overshoot': 0,
>>> 'Undershoot': 0.6864433321178778,
>>> 'Peak': 1.3994264225116284,
>>> 'PeakTime': 24.13227287437709,
>>> 'SteadyStateValue': -1.4010327022375215}
'''
if T is None:
T, y = step_response(sys)
elif np.asarray(T).size == 1:
T = _default_time_vector(sys, N=T_num, tfinal=T, is_step=True)
T, y = step_response(sys, T)
else:
T, y = step_response(sys, T)
# Steady state value
yfinal = sys.dcgain().real
# TODO: this could be a function step_info4data(t,y,yfinal)
# y: Step-response data, specified as:
# For SISO response data, a vector of length Ns, where Ns is the number of samples in the response data.
# For MIMO response data, an Ny-by-Nu-by-Ns array, where Ny is the number of system outputs, and Nu is the number of system inputs.
# T: Time vector corresponding to the response data in y, specified as a vector of length Ns.
# TODO: check dimentions for time vs step output and final output
if len(y.shape) == 3: # MIMO
Ny,Nu,Ns = y.shape
elif len(y.shape) == 1: # SISO
Ns = y.shape
Ny = 1
Nu = 1
else:
raise TypeError("Poner un mensaje")
S = [[[] for x in range(Ny)] for x in range(Nu)]
# Struct for each couple of inputs
for i in range(Ny): # Ny
for j in range(Nu): # Nu
rise_time: float = np.NaN
settling_time: float = np.NaN
settling_min: float = np.NaN
settling_max: float = np.NaN
peak_value: float = np.Inf
peak_time: float = np.Inf
undershoot: float = np.NaN
overshoot: float = np.NaN
steady_state_value: float = np.NaN
if len(y.shape) > 1:
yout = y[i,j,:]
InfValue = yfinal[i,j]
else:
yout = y
InfValue = yfinal
InfValue_sign = np.sign(InfValue)
if not np.isnan(InfValue) and not np.isinf(InfValue):
# SteadyStateValue
steady_state_value = InfValue
# Peak
peak_index = np.abs(yout).argmax()
peak_value = np.abs(yout[peak_index])
peak_time = T[peak_index]
sup_margin = (1. + SettlingTimeThreshold) * InfValue
inf_margin = (1. - SettlingTimeThreshold) * InfValue
# RiseTime
tr_lower_index = (np.where(InfValue_sign * (yout- RiseTimeLimits[0] * InfValue) >= 0 )[0])[0]
tr_upper_index = (np.where(InfValue_sign * yout >= InfValue_sign * RiseTimeLimits[1] * InfValue)[0])[0]
# SettlingTime
settling_time = T[np.where(np.abs(yout-InfValue) >= np.abs(SettlingTimeThreshold*InfValue))[0][-1]+1]
# Overshoot and Undershoot
y_os = (InfValue_sign*yout).max()
dy_os = np.abs(y_os) - np.abs(InfValue)
if dy_os > 0:
overshoot = np.abs(100. * dy_os / InfValue)
else:
overshoot = 0
y_us = (InfValue_sign*yout).min()
y_us_index = (InfValue_sign*yout).argmin()
if (InfValue_sign * yout[y_us_index]) < 0: # must have oposite sign
undershoot = np.abs(100. * np.abs(y_us) / InfValue)
else:
undershoot = 0
# RiseTime
rise_time = T[tr_upper_index] - T[tr_lower_index]
settling_max = (yout[tr_upper_index:]).max()
settling_min = (yout[tr_upper_index:]).min()
S[i][j] = {
'RiseTime': rise_time,
'SettlingTime': settling_time,
'SettlingMin': settling_min,
'SettlingMax': settling_max,
'Overshoot': overshoot,
'Undershoot': undershoot,
'Peak': peak_value,
'PeakTime': peak_time,
'SteadyStateValue': steady_state_value
}
if Ny > 1 or Nu > 1:
return S
else:
return S[0][0]
def initial_response(sys, T=None, X0=0., input=0, output=None, T_num=None,
transpose=False, return_x=False, squeeze=None):
# pylint: disable=W0622
"""Initial condition response of a linear system
If the system has multiple outputs (MIMO), optionally, one output
may be selected. If no selection is made for the output, all
outputs are given.
For information on the **shape** of parameters `T`, `X0` and
return values `T`, `yout`, see :ref:`time-series-convention`.
Parameters
----------
sys : StateSpace or TransferFunction
LTI system to simulate
T : array_like or float, optional
Time vector, or simulation time duration if a number (time vector is
autocomputed if not given; see :func:`step_response` for more detail)
X0 : array_like or float, optional
Initial condition (default = 0). Numbers are converted to constant
arrays with the correct shape.
input : int
Ignored, has no meaning in initial condition calculation. Parameter
ensures compatibility with step_response and impulse_response.
output : int
Index of the output that will be used in this simulation. Set to None
to not trim outputs.
T_num : int, optional
Number of time steps to use in simulation if T is not provided as an
array (autocomputed if not given); ignored if sys is discrete-time.
transpose : bool, optional
If True, transpose all input and output arrays (for backward
compatibility with MATLAB and :func:`scipy.signal.lsim`). Default
value is False.
return_x : bool, optional
If True, return the state vector (default = False).
squeeze : bool, optional
By default, if a system is single-input, single-output (SISO) then the
output response is returned as a 1D array (indexed by time). If
squeeze=True, remove single-dimensional entries from the shape of the
output even if the system is not SISO. If squeeze=False, keep the
output as a 2D array (indexed by the output number and time) even if
the system is SISO. The default value can be set using
config.defaults['control.squeeze_time_response'].
Returns
-------
T : array
Time values of the output
yout : array
Response of the system. If the system is SISO and squeeze is not
True, the array is 1D (indexed by time). If the system is not SISO or
squeeze is False, the array is 2D (indexed by the output number and
time).
xout : array, optional
Individual response of each x variable (if return_x is True).
See Also
--------
forced_response, impulse_response, step_response
Notes
-----
This function uses the `forced_response` function with the input set to
zero.
Examples
--------
>>> T, yout = initial_response(sys, T, X0)
"""
squeeze, sys = _get_ss_simo(sys, input, output, squeeze=squeeze)
# Create time and input vectors; checking is done in forced_response(...)
# The initial vector X0 is created in forced_response(...) if necessary
if T is None or np.asarray(T).size == 1:
T = _default_time_vector(sys, N=T_num, tfinal=T, is_step=False)
U = np.zeros_like(T)