Repository navigation
Expand file tree
/
Copy pathlti.py
More file actions
748 lines (597 loc) · 24.4 KB
/
Copy pathlti.py
File metadata and controls
748 lines (597 loc) · 24.4 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
# lti.py - LTI class and functions for linear systems
"""LTI class and functions for linear systems.
This module contains the LTI parent class to the child classes
StateSpace and TransferFunction.
"""
import math
from warnings import warn
import numpy as np
from numpy import abs, real
from . import config
from .iosys import InputOutputSystem
__all__ = ['poles', 'zeros', 'damp', 'evalfr', 'frequency_response',
'freqresp', 'dcgain', 'bandwidth', 'LTI']
class LTI(InputOutputSystem):
"""Parent class for linear time-invariant system objects.
LTI is the parent to the `FrequencyResponseData`, `StateSpace`, and
`TransferFunction` child classes. It contains the number of inputs and
outputs, and the timebase (dt) for the system. This class is not
generally accessed directly by the user.
See Also
--------
InputOutputSystem, StateSpace, TransferFunction, FrequencyResponseData
"""
def __init__(self, inputs=1, outputs=1, states=None, name=None, **kwargs):
"""Assign the LTI object's numbers of inputs and outputs."""
super().__init__(
name=name, inputs=inputs, outputs=outputs, states=states, **kwargs)
def __call__(self, x, squeeze=None, warn_infinite=True):
"""Evaluate system transfer function at point in complex plane.
Returns the value of the system's transfer function at a point `x`
in the complex plane, where `x` is `s` for continuous-time systems
and `z` for discrete-time systems.
By default, a (complex) scalar will be returned for SISO systems
and a p x m array will be return for MIMO systems with m inputs and
p outputs. This can be changed using the `squeeze` keyword.
To evaluate at a frequency `omega` in radians per second,
enter ``x = omega * 1j`` for continuous-time systems,
``x = exp(1j * omega * dt)`` for discrete-time systems, or
use the `~LTI.frequency_response` method.
Parameters
----------
x : complex or complex 1D array_like
Complex value(s) at which transfer function will be evaluated.
squeeze : bool, optional
Squeeze output, as described below. Default value can be set
using `config.defaults['control.squeeze_frequency_response']`.
warn_infinite : bool, optional
If set to False, turn off divide by zero warning.
Returns
-------
fresp : complex ndarray
The value of the system transfer function at `x`. If the system
is SISO and `squeeze` is not True, the shape of the array matches
the shape of `x`. If the system is not SISO or `squeeze` is
False, the first two dimensions of the array are indices for the
output and input and the remaining dimensions match `x`. If
`squeeze` is True then single-dimensional axes are removed.
Notes
-----
See `FrequencyResponseData.__call__`, `StateSpace.__call__`,
`TransferFunction.__call__` for class-specific details.
"""
raise NotImplementedError("not implemented in subclass")
def damp(self):
"""Natural frequency, damping ratio of system poles.
Returns
-------
wn : array
Natural frequency for each system pole.
zeta : array
Damping ratio for each system pole.
poles : array
System pole locations.
"""
poles = self.poles()
if self.isdtime(strict=True):
splane_poles = np.log(poles.astype(complex))/self.dt
else:
splane_poles = poles
wn = abs(splane_poles)
zeta = -real(splane_poles)/wn
return wn, zeta, poles
def feedback(self, other=1, sign=-1):
"""Feedback interconnection between two input/output systems.
Parameters
----------
other : `InputOutputSystem`
System in the feedback path.
sign : float, optional
Gain to use in feedback path. Defaults to -1.
"""
raise NotImplementedError("feedback not implemented in subclass")
def frequency_response(self, omega=None, squeeze=None):
"""Evaluate LTI system response at an array of frequencies.
See `frequency_response` for more detailed information.
"""
from .frdata import FrequencyResponseData
if omega is None:
# Use default frequency range
from .freqplot import _default_frequency_range
omega = _default_frequency_range(self)
omega = np.sort(np.array(omega, ndmin=1))
if self.isdtime(strict=True):
# Convert the frequency to discrete time
if np.any(omega * self.dt > np.pi):
warn("__call__: evaluation above Nyquist frequency")
s = np.exp(1j * omega * self.dt)
else:
s = 1j * omega
# Return the data as a frequency response data object
response = self(s)
return FrequencyResponseData(
response, omega, return_magphase=True, squeeze=squeeze,
dt=self.dt, sysname=self.name, inputs=self.input_labels,
outputs=self.output_labels, plot_type='bode')
def dcgain(self):
"""Return the zero-frequency (DC) gain."""
raise NotImplementedError("dcgain not defined for subclass")
def _dcgain(self, warn_infinite):
zeroresp = self(0 if self.isctime() else 1,
warn_infinite=warn_infinite)
if np.all(np.logical_or(np.isreal(zeroresp), np.isnan(zeroresp.imag))):
return zeroresp.real
else:
return zeroresp
def bandwidth(self, dbdrop=-3):
"""Evaluate bandwidth of an LTI system for a given dB drop.
Evaluate the first frequency that the response magnitude is lower than
DC gain by `dbdrop` dB.
Parameters
----------
dbdrop : float, optional
A strictly negative scalar in dB (default = -3) defines the
amount of gain drop for deciding bandwidth.
Returns
-------
bandwidth : ndarray
The first frequency (rad/time-unit) where the gain drops below
`dbdrop` of the dc gain of the system, or nan if the system has
infinite dc gain, inf if the gain does not drop for all frequency.
Raises
------
TypeError
If `sys` is not an SISO LTI instance.
ValueError
If `dbdrop` is not a negative scalar.
"""
# check if system is SISO and dbdrop is a negative scalar
if not self.issiso():
raise TypeError("system should be a SISO system")
if (not np.isscalar(dbdrop)) or dbdrop >= 0:
raise ValueError("expecting dbdrop be a negative scalar in dB")
dcgain = self.dcgain()
if np.isinf(dcgain):
# infinite dcgain, return np.nan
return np.nan
# use frequency range to identify the 0-crossing (dbdrop) bracket
from control.freqplot import _default_frequency_range
omega = _default_frequency_range(self)
mag, phase, omega = self.frequency_response(omega)
idx_dropped = np.nonzero(mag - dcgain*10**(dbdrop/20) < 0)[0]
if idx_dropped.shape[0] == 0:
# no frequency response is dbdrop below the dc gain, return np.inf
return np.inf
else:
# solve for the bandwidth, use scipy.optimize.root_scalar() to
# solve using bisection
import scipy
result = scipy.optimize.root_scalar(
lambda w: np.abs(self(w*1j)) - np.abs(dcgain)*10**(dbdrop/20),
bracket=[omega[idx_dropped[0] - 1], omega[idx_dropped[0]]],
method='bisect')
# check solution
if result.converged:
return np.abs(result.root)
else:
raise Exception(result.message)
def ispassive(self):
r"""Indicate if a linear time invariant (LTI) system is passive.
See `ispassive` for details.
"""
# importing here prevents circular dependency
from control.passivity import ispassive
return ispassive(self)
#
# Convenience aliases for conversion functions
#
# Allow conversion between state space and transfer function types
# as methods. These are just pass throughs to factory functions.
#
# Note: in order for docstrings to created, these have to set these up
# as independent methods, not just assigned to ss() and tf().
#
# Imports are done within the function to avoid circular imports.
#
def to_ss(self, *args, **kwargs):
"""Convert to state space representation.
See `ss` for details.
"""
from .statesp import ss
return ss(self, *args, **kwargs)
def to_tf(self, *args, **kwargs):
"""Convert to transfer function representation.
See `tf` for details.
"""
from .xferfcn import tf
return tf(self, *args, **kwargs)
#
# Convenience aliases for plotting and response functions
#
# Allow standard plots to be generated directly from the system object
# in addition to standalone plotting and response functions.
#
# Note: in order for docstrings to created, these have to set these up as
# independent methods, not just assigned to plotting/response functions.
#
# Imports are done within the function to avoid circular imports.
#
def bode_plot(self, *args, **kwargs):
"""Generate a Bode plot for the system.
See `bode_plot` for more information.
"""
from .freqplot import bode_plot
return bode_plot(self, *args, **kwargs)
def nichols_plot(self, *args, **kwargs):
"""Generate a Nichols plot for the system.
See `nichols_plot` for more information.
"""
from .nichols import nichols_plot
return nichols_plot(self, *args, **kwargs)
def nyquist_plot(self, *args, **kwargs):
"""Generate a Nyquist plot for the system.
See `nyquist_plot` for more information.
"""
from .freqplot import nyquist_plot
return nyquist_plot(self, *args, **kwargs)
def forced_response(self, *args, **kwargs):
"""Generate the forced response for the system.
See `forced_response` for more information.
"""
from .timeresp import forced_response
return forced_response(self, *args, **kwargs)
def impulse_response(self, *args, **kwargs):
"""Generate the impulse response for the system.
See `impulse_response` for more information.
"""
from .timeresp import impulse_response
return impulse_response(self, *args, **kwargs)
def initial_response(self, *args, **kwargs):
"""Generate the initial response for the system.
See `initial_response` for more information.
"""
from .timeresp import initial_response
return initial_response(self, *args, **kwargs)
def step_response(self, *args, **kwargs):
"""Generate the step response for the system.
See `step_response` for more information.
"""
from .timeresp import step_response
return step_response(self, *args, **kwargs)
def poles(sys):
"""
Compute system poles.
Parameters
----------
sys : `StateSpace` or `TransferFunction`
Linear system.
Returns
-------
poles : ndarray
Array that contains the system's poles.
See Also
--------
zeros, StateSpace.poles, TransferFunction.poles
"""
return sys.poles()
def zeros(sys):
"""
Compute system zeros.
Parameters
----------
sys : `StateSpace` or `TransferFunction`
Linear system.
Returns
-------
zeros : ndarray
Array that contains the system's zeros.
See Also
--------
poles, StateSpace.zeros, TransferFunction.zeros
"""
return sys.zeros()
def damp(sys, doprint=True):
"""Compute system's natural frequencies, damping ratios, and poles.
Parameters
----------
sys : `StateSpace` or `TransferFunction`
A linear system object.
doprint : bool (optional)
If True, print table with values.
Returns
-------
wn : array
Natural frequency for each system pole.
zeta : array
Damping ratio for each system pole.
poles : array
System pole locations.
See Also
--------
poles
Notes
-----
If the system is continuous
| ``wn = abs(poles)``
| ``zeta = -real(poles)/poles``
If the system is discrete, the discrete poles are mapped to their
equivalent location in the s-plane via
| ``s = log(poles)/dt``
and
| ``wn = abs(s)``
| ``zeta = -real(s)/wn``
Examples
--------
>>> G = ct.tf([1], [1, 4])
>>> wn, zeta, poles = ct.damp(G)
Eigenvalue (pole) Damping Frequency
-4 1 4
"""
wn, zeta, poles = sys.damp()
if doprint:
print(' Eigenvalue (pole) Damping Frequency')
for p, z, w in zip(poles, zeta, wn):
if abs(p.imag) < 1e-12:
print(" %10.4g %10.4g %10.4g" %
(p.real, 1.0, w))
else:
print("%10.4g%+10.4gj %10.4g %10.4g" %
(p.real, p.imag, z, w))
return wn, zeta, poles
# TODO: deprecate this function
def evalfr(sys, x, squeeze=None):
"""Evaluate transfer function of LTI system at complex frequency.
Returns the complex frequency response ``sys(x)`` where `x` is `s` for
continuous-time systems and `z` for discrete-time systems, with
``m = sys.ninputs`` number of inputs and ``p = sys.noutputs`` number of
outputs.
To evaluate at a frequency omega in radians per second, enter
``x = omega * 1j`` for continuous-time systems, or
``x = exp(1j * omega * dt)`` for discrete-time systems, or use
``freqresp(sys, omega)``.
Parameters
----------
sys : `StateSpace` or `TransferFunction`
Linear system.
x : complex scalar or 1D array_like
Complex frequency(s).
squeeze : bool, optional (default=True)
If `squeeze` = True, remove single-dimensional entries from the
shape of the output even if the system is not SISO. If
`squeeze` = False, keep all indices (output, input and, if omega is
array_like, frequency) even if the system is SISO. The default
value can be set using
`config.defaults['control.squeeze_frequency_response']`.
Returns
-------
fresp : complex ndarray
The frequency response of the system. If the system is SISO and
squeeze is not True, the shape of the array matches the shape of
omega. If the system is not SISO or squeeze is False, the first two
dimensions of the array are indices for the output and input and the
remaining dimensions match omega. If `squeeze` is True then
single-dimensional axes are removed.
See Also
--------
LTI.__call__, frequency_response, bode_plot
Notes
-----
This function is a wrapper for `StateSpace.__call__` and
`TransferFunction.__call__`.
Examples
--------
>>> G = ct.ss([[-1, -2], [3, -4]], [[5], [7]], [[6, 8]], [[9]])
>>> fresp = ct.evalfr(G, 1j) # evaluate at s = 1j
"""
return sys(x, squeeze=squeeze)
def frequency_response(
sysdata, omega=None, omega_limits=None, omega_num=None,
Hz=None, squeeze=None):
"""Frequency response of an LTI system.
For continuous-time systems with transfer function G, computes the
frequency response as
G(j*omega) = mag * exp(j*phase)
For discrete-time systems, the response is evaluated around the unit
circle such that
G(exp(j*omega*dt)) = mag * exp(j*phase).
In general the system may be multiple input, multiple output (MIMO),
where ``m = self.ninputs`` number of inputs and ``p = self.noutputs``
number of outputs.
Parameters
----------
sysdata : LTI system or list of LTI systems
Linear system(s) for which frequency response is computed.
omega : float or 1D array_like, optional
A list, tuple, array, or scalar value of frequencies in radians/sec
at which the system will be evaluated. Can be a single frequency
or array of frequencies, which will be sorted before evaluation.
If None (default), a common set of frequencies that works across
all given systems is computed.
omega_limits : array_like of two values, optional
Limits to the range of frequencies, in rad/sec. Specifying
`omega` as a list of two elements is equivalent to providing
`omega_limits`. Ignored if omega is provided.
omega_num : int, optional
Number of frequency samples at which to compute the response.
Defaults to `config.defaults['freqplot.number_of_samples']`. Ignored
if omega is provided.
Returns
-------
response : `FrequencyResponseData`
Frequency response data object representing the frequency
response. When accessed as a tuple, returns ``(magnitude,
phase, omega)``. If `sysdata` is a list of systems, returns a
`FrequencyResponseList` object. Results can be plotted using
the `~FrequencyResponseData.plot` method. See
`FrequencyResponseData` for more detailed information.
response.magnitude : array
Magnitude of the frequency response (absolute value, not dB or
log10). If the system is SISO and squeeze is not True, the
array is 1D, indexed by frequency. If the system is not SISO
or squeeze is False, the array is 3D, indexed by the output,
input, and, if omega is array_like, frequency. If `squeeze` is
True then single-dimensional axes are removed.
response.phase : array
Wrapped phase, in radians, with same shape as `magnitude`.
response.omega : array
Sorted list of frequencies at which response was evaluated.
Other Parameters
----------------
Hz : bool, optional
If True, when computing frequency limits automatically set
limits to full decades in Hz instead of rad/s. Omega is always
returned in rad/sec.
squeeze : bool, optional
If `squeeze` = True, remove single-dimensional entries from the
shape of the output even if the system is not SISO. If
`squeeze` = False, keep all indices (output, input and, if omega is
array_like, frequency) even if the system is SISO. The default
value can be set using
`config.defaults['control.squeeze_frequency_response']`.
See Also
--------
LTI.__call__, bode_plot
Notes
-----
This function is a wrapper for `StateSpace.frequency_response` and
`TransferFunction.frequency_response`. You can also use the
lower-level methods ``sys(s)`` or ``sys(z)`` to generate the frequency
response for a single system.
All frequency data should be given in rad/sec. If frequency limits are
computed automatically, the `Hz` keyword can be used to ensure that
limits are in factors of decades in Hz, so that Bode plots with
`Hz` = True look better.
The frequency response data can be plotted by calling the `bode_plot`
function or using the `plot` method of the `FrequencyResponseData`
class.
Examples
--------
>>> G = ct.ss([[-1, -2], [3, -4]], [[5], [7]], [[6, 8]], [[9]])
>>> mag, phase, omega = ct.frequency_response(G, [0.1, 1., 10.])
>>> sys = ct.rss(3, 2, 2)
>>> mag, phase, omega = ct.frequency_response(sys, [0.1, 1., 10.])
>>> mag[0, 1, :] # Magnitude of second input to first output
array([..., ..., ...])
>>> phase[1, 0, :] # Phase of first input to second output
array([..., ..., ...])
"""
from .frdata import FrequencyResponseData
from .freqplot import _determine_omega_vector
# Process keyword arguments
omega_num = config._get_param('freqplot', 'number_of_samples', omega_num)
# Convert the first argument to a list
syslist = sysdata if isinstance(sysdata, (list, tuple)) else [sysdata]
# Get the common set of frequencies to use
omega_syslist, omega_range_given = _determine_omega_vector(
syslist, omega, omega_limits, omega_num, Hz=Hz)
responses = []
for sys_ in syslist:
if isinstance(sys_, FrequencyResponseData) and sys_._ifunc is None \
and not omega_range_given:
omega_sys = sys_.omega # use system properties
else:
omega_sys = omega_syslist.copy() # use common omega vector
# Add the Nyquist frequency for discrete-time systems
if sys_.isdtime(strict=True):
nyquistfrq = math.pi / sys_.dt
if not omega_range_given:
# Limit up to the Nyquist frequency
omega_sys = omega_sys[omega_sys < nyquistfrq]
# Compute the frequency response
responses.append(sys_.frequency_response(omega_sys, squeeze=squeeze))
if isinstance(sysdata, (list, tuple)):
from .freqplot import FrequencyResponseList
return FrequencyResponseList(responses)
else:
return responses[0]
# Alternative name (legacy)
def freqresp(sys, omega):
"""Legacy version of frequency_response.
.. deprecated:: 0.9.0
This function will be removed in a future version of python-control.
Use `frequency_response` instead.
"""
warn("freqresp() is deprecated; use frequency_response()", FutureWarning)
return frequency_response(sys, omega)
def dcgain(sys):
"""Return the zero-frequency (or DC) gain of the given system.
Parameters
----------
sys : LTI
System for which the zero-frequency gain is computed.
Returns
-------
gain : ndarray
The zero-frequency gain, or (inf + nanj) if the system has a pole at
the origin, (nan + nanj) if there is a pole/zero cancellation at the
origin.
Examples
--------
>>> G = ct.tf([1], [1, 2])
>>> ct.dcgain(G) # doctest: +SKIP
np.float(0.5)
"""
return sys.dcgain()
def bandwidth(sys, dbdrop=-3):
"""Find first frequency where gain drops by 3 dB.
Parameters
----------
sys : `StateSpace` or `TransferFunction`
Linear system for which the bandwidth should be computed.
dbdrop : float, optional
By how much the gain drop in dB (default = -3) that defines the
bandwidth. Should be a negative scalar.
Returns
-------
bandwidth : ndarray
The first frequency where the gain drops below `dbdrop` of the zero
frequency (DC) gain of the system, or nan if the system has infinite
zero frequency gain, inf if the gain does not drop for any frequency.
Raises
------
TypeError
If `sys` is not an SISO LTI instance.
ValueError
If `dbdrop` is not a negative scalar.
Examples
--------
>>> G = ct.tf([1], [1, 1])
>>> ct.bandwidth(G)
np.float64(0.9976283451102316)
>>> G1 = ct.tf(0.1, [1, 0.1])
>>> wn2 = 1
>>> zeta2 = 0.001
>>> G2 = ct.tf(wn2**2, [1, 2*zeta2*wn2, wn2**2])
>>> ct.bandwidth(G1*G2)
np.float64(0.10184838823897456)
"""
if not isinstance(sys, LTI):
raise TypeError("sys must be a LTI instance.")
return sys.bandwidth(dbdrop)
# Process frequency responses in a uniform way
def _process_frequency_response(sys, omega, out, squeeze=None):
# Set value of squeeze argument if not set
if squeeze is None:
squeeze = config.defaults['control.squeeze_frequency_response']
if np.asarray(omega).ndim < 1:
# received a scalar x, squeeze down the array along last dim
out = np.squeeze(out, axis=2)
#
# Get rid of unneeded dimensions
#
# There are three possible values for the squeeze keyword at this point:
#
# squeeze=None: squeeze input/output axes iff SISO
# squeeze=True: squeeze all single dimensional axes (ala numpy)
# squeeze-False: don't squeeze any axes
#
if squeeze is True:
# Squeeze everything that we can if that's what the user wants
return np.squeeze(out)
elif squeeze is None and sys.issiso():
# SISO system output squeezed unless explicitly specified otherwise
return out[0][0]
elif squeeze is False or squeeze is None:
return out
else:
raise ValueError("unknown squeeze value")