Repository navigation
Expand file tree
/
Copy pathfreqplot.py
More file actions
2848 lines (2476 loc) · 118 KB
/
Copy pathfreqplot.py
File metadata and controls
2848 lines (2476 loc) · 118 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
# freqplot.py - frequency domain plots for control systems
#
# Initial author: Richard M. Murray
# Date: 24 May 09
#
# This file contains some standard control system plots: Bode plots,
# Nyquist plots and other frequency response plots. The code for Nichols
# charts is in nichols.py. The code for pole-zero diagrams is in pzmap.py
# and rlocus.py.
import itertools
import math
import warnings
import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
from . import config
from .bdalg import feedback
from .ctrlplot import ControlPlot, _add_arrows_to_line2D, _find_axes_center, \
_get_color, _get_color_offset, _get_line_labels, _make_legend_labels, \
_process_ax_keyword, _process_legend_keywords, _process_line_labels, \
_update_plot_title
from .ctrlutil import unwrap
from .exception import ControlMIMONotImplemented
from .frdata import FrequencyResponseData
from .lti import LTI, _process_frequency_response, frequency_response
from .margins import stability_margins
from .statesp import StateSpace
from .xferfcn import TransferFunction
__all__ = ['bode_plot', 'NyquistResponseData', 'nyquist_response',
'nyquist_plot', 'singular_values_response',
'singular_values_plot', 'gangof4_plot', 'gangof4_response',
'bode', 'nyquist', 'gangof4', 'FrequencyResponseList']
# Default values for module parameter variables
_freqplot_defaults = {
'freqplot.feature_periphery_decades': 1,
'freqplot.number_of_samples': 1000,
'freqplot.dB': False, # Plot gain in dB
'freqplot.deg': True, # Plot phase in degrees
'freqplot.Hz': False, # Plot frequency in Hertz
'freqplot.grid': True, # Turn on grid for gain and phase
'freqplot.wrap_phase': False, # Wrap the phase plot at a given value
'freqplot.freq_label': "Frequency [{units}]",
'freqplot.share_magnitude': 'row',
'freqplot.share_phase': 'row',
'freqplot.share_frequency': 'col',
'freqplot.title_frame': 'axes',
}
#
# Frequency response data list class
#
# This class is a subclass of list that adds a plot() method, enabling
# direct plotting from routines returning a list of FrequencyResponseData
# objects.
#
class FrequencyResponseList(list):
def plot(self, *args, plot_type=None, **kwargs):
if plot_type == None:
for response in self:
if plot_type is not None and response.plot_type != plot_type:
raise TypeError(
"inconsistent plot_types in data; set plot_type "
"to 'bode', 'nichols', or 'svplot'")
plot_type = response.plot_type
# Use FRD plot method, which can handle lists via plot functions
return FrequencyResponseData.plot(
self, plot_type=plot_type, *args, **kwargs)
#
# Bode plot
#
# This is the default method for plotting frequency responses. There are
# lots of options available for tuning the format of the plot, (hopefully)
# covering most of the common use cases.
#
def bode_plot(
data, omega=None, *fmt, ax=None, omega_limits=None, omega_num=None,
plot=None, plot_magnitude=True, plot_phase=None,
overlay_outputs=None, overlay_inputs=None, phase_label=None,
magnitude_label=None, label=None, display_margins=None,
margins_method='best', title=None, sharex=None, sharey=None, **kwargs):
"""Bode plot for a system.
Plot the magnitude and phase of the frequency response over a
(optional) frequency range.
Parameters
----------
data : list of `FrequencyResponseData` or `LTI`
List of LTI systems or :class:`FrequencyResponseData` objects. A
single system or frequency response can also be passed.
omega : array_like, optoinal
Set of frequencies in rad/sec to plot over. If not specified, this
will be determined from the proporties of the systems. Ignored if
`data` is not a list of systems.
*fmt : :func:`matplotlib.pyplot.plot` format string, optional
Passed to `matplotlib` as the format string for all lines in the plot.
The `omega` parameter must be present (use omega=None if needed).
dB : bool
If True, plot result in dB. Default is False.
Hz : bool
If True, plot frequency in Hz (omega must be provided in rad/sec).
Default value (False) set by config.defaults['freqplot.Hz'].
deg : bool
If True, plot phase in degrees (else radians). Default value (True)
set by config.defaults['freqplot.deg'].
display_margins : bool or str
If True, draw gain and phase margin lines on the magnitude and phase
graphs and display the margins at the top of the graph. If set to
'overlay', the values for the gain and phase margin are placed on
the graph. Setting display_margins turns off the axes grid.
**kwargs : :func:`matplotlib.pyplot.plot` keyword properties, optional
Additional keywords passed to `matplotlib` to specify line properties.
Returns
-------
cplt : :class:`ControlPlot` object
Object containing the data that were plotted:
* cplt.lines: Array of :class:`matplotlib.lines.Line2D` objects
for each line in the plot. The shape of the array matches the
subplots shape and the value of the array is a list of Line2D
objects in that subplot.
* cplt.axes: 2D array of :class:`matplotlib.axes.Axes` for the plot.
* cplt.figure: :class:`matplotlib.figure.Figure` containing the plot.
* cplt.legend: legend object(s) contained in the plot
See :class:`ControlPlot` for more detailed information.
Other Parameters
----------------
ax : array of matplotlib.axes.Axes, optional
The matplotlib axes to draw the figure on. If not specified, the
axes for the current figure are used or, if there is no current
figure with the correct number and shape of axes, a new figure is
created. The shape of the array must match the shape of the
plotted data.
freq_label, magnitude_label, phase_label : str, optional
Labels to use for the frequency, magnitude, and phase axes.
Defaults are set by `config.defaults['freqplot.<keyword>']`.
grid : bool, optional
If True, plot grid lines on gain and phase plots. Default is set by
`config.defaults['freqplot.grid']`.
initial_phase : float, optional
Set the reference phase to use for the lowest frequency. If set, the
initial phase of the Bode plot will be set to the value closest to the
value specified. Units are in either degrees or radians, depending on
the `deg` parameter. Default is -180 if wrap_phase is False, 0 if
wrap_phase is True.
label : str or array_like of str, optional
If present, replace automatically generated label(s) with the given
label(s). If sysdata is a list, strings should be specified for each
system. If MIMO, strings required for each system, output, and input.
legend_map : array of str, optional
Location of the legend for multi-axes plots. Specifies an array
of legend location strings matching the shape of the subplots, with
each entry being either None (for no legend) or a legend location
string (see :func:`~matplotlib.pyplot.legend`).
legend_loc : int or str, optional
Include a legend in the given location. Default is 'center right',
with no legend for a single response. Use False to suppress legend.
margins_method : str, optional
Method to use in computing margins (see :func:`stability_margins`).
omega_limits : array_like of two values
Set limits for plotted frequency range. If Hz=True the limits are
in Hz otherwise in rad/s. Specifying ``omega`` as a list of two
elements is equivalent to providing ``omega_limits``. Ignored if
data is not a list of systems.
omega_num : int
Number of samples to use for the frequeny range. Defaults to
config.defaults['freqplot.number_of_samples']. Ignored if data is
not a list of systems.
overlay_inputs, overlay_outputs : bool, optional
If set to True, combine input and/or output signals onto a single
plot and use line colors, labels, and a legend to distinguish them.
plot : bool, optional
(legacy) If given, `bode_plot` returns the legacy return values
of magnitude, phase, and frequency. If False, just return the
values with no plot.
plot_magnitude, plot_phase : bool, optional
If set to `False`, don't plot the magnitude or phase, respectively.
rcParams : dict
Override the default parameters used for generating plots.
Default is set by config.default['ctrlplot.rcParams'].
share_frequency, share_magnitude, share_phase : str or bool, optional
Determine whether and how axis limits are shared between the
indicated variables. Can be set set to 'row' to share across all
subplots in a row, 'col' to set across all subplots in a column, or
`False` to allow independent limits. Note: if `sharex` is given,
it sets the value of `share_frequency`; if `sharey` is given, it
sets the value of both `share_magnitude` and `share_phase`.
Default values are 'row' for `share_magnitude` and `share_phase',
'col', for `share_frequency`, and can be set using
config.defaults['freqplot.share_<axis>'].
show_legend : bool, optional
Force legend to be shown if ``True`` or hidden if ``False``. If
``None``, then show legend when there is more than one line on an
axis or ``legend_loc`` or ``legend_map`` has been specified.
title : str, optional
Set the title of the plot. Defaults to plot type and system name(s).
title_frame : str, optional
Set the frame of reference used to center the plot title. If set to
'axes' (default), the horizontal position of the title will be
centered relative to the axes. If set to 'figure', it will be
centered with respect to the figure (faster execution).
wrap_phase : bool or float
If wrap_phase is `False` (default), then the phase will be unwrapped
so that it is continuously increasing or decreasing. If wrap_phase is
`True` the phase will be restricted to the range [-180, 180) (or
[:math:`-\\pi`, :math:`\\pi`) radians). If `wrap_phase` is specified
as a float, the phase will be offset by 360 degrees if it falls below
the specified value. Default value is `False` and can be set using
config.defaults['freqplot.wrap_phase'].
See Also
--------
frequency_response
Notes
-----
Starting with python-control version 0.10, `bode_plot` returns a
:class:`ControlPlot` object instead of magnitude, phase, and
frequency. To recover the old behavior, call `bode_plot` with
`plot=True`, which will force the legacy values (mag, phase, omega) to
be returned (with a warning). To obtain just the frequency response of
a system (or list of systems) without plotting, use the
:func:`~control.frequency_response` command.
If a discrete time model is given, the frequency response is plotted
along the upper branch of the unit circle, using the mapping ``z =
exp(1j * omega * dt)`` where `omega` ranges from 0 to `pi/dt` and `dt`
is the discrete timebase. If timebase not specified (``dt=True``),
`dt` is set to 1.
The default values for Bode plot configuration parameters can be reset
using the `config.defaults` dictionary, with module name 'bode'.
Examples
--------
>>> G = ct.ss([[-1, -2], [3, -4]], [[5], [7]], [[6, 8]], [[9]])
>>> out = ct.bode_plot(G)
"""
#
# Process keywords and set defaults
#
# Make a copy of the kwargs dictionary since we will modify it
kwargs = dict(kwargs)
# Get values for params (and pop from list to allow keyword use in plot)
dB = config._get_param(
'freqplot', 'dB', kwargs, _freqplot_defaults, pop=True)
deg = config._get_param(
'freqplot', 'deg', kwargs, _freqplot_defaults, pop=True)
Hz = config._get_param(
'freqplot', 'Hz', kwargs, _freqplot_defaults, pop=True)
grid = config._get_param(
'freqplot', 'grid', kwargs, _freqplot_defaults, pop=True)
wrap_phase = config._get_param(
'freqplot', 'wrap_phase', kwargs, _freqplot_defaults, pop=True)
initial_phase = config._get_param(
'freqplot', 'initial_phase', kwargs, None, pop=True)
rcParams = config._get_param('ctrlplot', 'rcParams', kwargs, pop=True)
title_frame = config._get_param(
'freqplot', 'title_frame', kwargs, _freqplot_defaults, pop=True)
# Set the default labels
freq_label = config._get_param(
'freqplot', 'freq_label', kwargs, _freqplot_defaults, pop=True)
if magnitude_label is None:
magnitude_label = "Magnitude [dB]" if dB else "Magnitude"
if phase_label is None:
phase_label = "Phase [deg]" if deg else "Phase [rad]"
# Use sharex and sharey as proxies for share_{magnitude, phase, frequency}
if sharey is not None:
if 'share_magnitude' in kwargs or 'share_phase' in kwargs:
ValueError(
"sharey cannot be present with share_magnitude/share_phase")
kwargs['share_magnitude'] = sharey
kwargs['share_phase'] = sharey
if sharex is not None:
if 'share_frequency' in kwargs:
ValueError(
"sharex cannot be present with share_frequency")
kwargs['share_frequency'] = sharex
# Legacy keywords for margins
display_margins = config._process_legacy_keyword(
kwargs, 'margins', 'display_margins', display_margins)
if kwargs.pop('margin_info', False):
warnings.warn(
"keyword 'margin_info' is deprecated; "
"use 'display_margins='overlay'")
if display_margins is False:
raise ValueError(
"conflicting_keywords: `display_margins` and `margin_info`")
margins_method = config._process_legacy_keyword(
kwargs, 'method', 'margins_method', margins_method)
if not isinstance(data, (list, tuple)):
data = [data]
#
# Pre-process the data to be plotted (unwrap phase, limit frequencies)
#
# To maintain compatibility with legacy uses of bode_plot(), we do some
# initial processing on the data, specifically phase unwrapping and
# setting the initial value of the phase. If bode_plot is called with
# plot == False, then these values are returned to the user (instead of
# the list of lines created, which is the new output for _plot functions.
#
# If we were passed a list of systems, convert to data
if any([isinstance(
sys, (StateSpace, TransferFunction)) for sys in data]):
data = frequency_response(
data, omega=omega, omega_limits=omega_limits,
omega_num=omega_num, Hz=Hz)
else:
# Generate warnings if frequency keywords were given
if omega_num is not None:
warnings.warn("`omega_num` ignored when passed response data")
elif omega is not None:
warnings.warn("`omega` ignored when passed response data")
# Check to make sure omega_limits is sensible
if omega_limits is not None and \
(len(omega_limits) != 2 or omega_limits[1] <= omega_limits[0]):
raise ValueError(f"invalid limits: {omega_limits=}")
# If plot_phase is not specified, check the data first, otherwise true
if plot_phase is None:
plot_phase = True if data[0].plot_phase is None else data[0].plot_phase
if not plot_magnitude and not plot_phase:
raise ValueError(
"plot_magnitude and plot_phase both False; no data to plot")
mag_data, phase_data, omega_data = [], [], []
for response in data:
noutputs, ninputs = response.noutputs, response.ninputs
if initial_phase is None:
# Start phase in the range 0 to -360 w/ initial phase = 0
# TODO: change this to 0 to 270 (?)
# If wrap_phase is true, use 0 instead (phase \in (-pi, pi])
initial_phase_value = -math.pi if wrap_phase is not True else 0
elif isinstance(initial_phase, (int, float)):
# Allow the user to override the default calculation
if deg:
initial_phase_value = initial_phase/180. * math.pi
else:
initial_phase_value = initial_phase
else:
raise ValueError("initial_phase must be a number.")
# Reshape the phase to allow standard indexing
phase = response.phase.copy().reshape((noutputs, ninputs, -1))
# Shift and wrap the phase
for i, j in itertools.product(range(noutputs), range(ninputs)):
# Shift the phase if needed
if abs(phase[i, j, 0] - initial_phase_value) > math.pi:
phase[i, j] -= 2*math.pi * round(
(phase[i, j, 0] - initial_phase_value) / (2*math.pi))
# Phase wrapping
if wrap_phase is False:
phase[i, j] = unwrap(phase[i, j]) # unwrap the phase
elif wrap_phase is True:
pass # default calc OK
elif isinstance(wrap_phase, (int, float)):
phase[i, j] = unwrap(phase[i, j]) # unwrap phase first
if deg:
wrap_phase *= math.pi/180.
# Shift the phase if it is below the wrap_phase
phase[i, j] += 2*math.pi * np.maximum(
0, np.ceil((wrap_phase - phase[i, j])/(2*math.pi)))
else:
raise ValueError("wrap_phase must be bool or float.")
# Put the phase back into the original shape
phase = phase.reshape(response.magnitude.shape)
# Save the data for later use (legacy return values)
mag_data.append(response.magnitude)
phase_data.append(phase)
omega_data.append(response.omega)
#
# Process `plot` keyword
#
# We use the `plot` keyword to track legacy usage of `bode_plot`.
# Prior to v0.10, the `bode_plot` command returned mag, phase, and
# omega. Post v0.10, we return an array with the same shape as the
# axes we use for plotting, with each array element containing a list
# of lines drawn on that axes.
#
# There are three possibilities at this stage in the code:
#
# * plot == True: set explicitly by the user. Return mag, phase, omega,
# with a warning.
#
# * plot == False: set explicitly by the user. Return mag, phase,
# omega, with a warning.
#
# * plot == None: this is the new default setting. Return an array of
# lines that were drawn.
#
# If `bode_plot` was called with no `plot` argument and the return
# values were used, the new code will cause problems (you get an array
# of lines instead of magnitude, phase, and frequency). To recover the
# old behavior, call `bode_plot` with `plot=True`.
#
# All of this should be removed in v0.11+ when we get rid of deprecated
# code.
#
if plot is not None:
warnings.warn(
"bode_plot() return value of mag, phase, omega is deprecated; "
"use frequency_response()", FutureWarning)
if plot is False:
# Process the data to match what we were sent
for i in range(len(mag_data)):
mag_data[i] = _process_frequency_response(
data[i], omega_data[i], mag_data[i], squeeze=data[i].squeeze)
phase_data[i] = _process_frequency_response(
data[i], omega_data[i], phase_data[i], squeeze=data[i].squeeze)
if len(data) == 1:
return mag_data[0], phase_data[0], omega_data[0]
else:
return mag_data, phase_data, omega_data
#
# Find/create axes
#
# Data are plotted in a standard subplots array, whose size depends on
# which signals are being plotted and how they are combined. The
# baseline layout for data is to plot everything separately, with
# the magnitude and phase for each output making up the rows and the
# columns corresponding to the different inputs.
#
# Input 0 Input m
# +---------------+ +---------------+
# | mag H_y0,u0 | ... | mag H_y0,um |
# +---------------+ +---------------+
# +---------------+ +---------------+
# | phase H_y0,u0 | ... | phase H_y0,um |
# +---------------+ +---------------+
# : :
# +---------------+ +---------------+
# | mag H_yp,u0 | ... | mag H_yp,um |
# +---------------+ +---------------+
# +---------------+ +---------------+
# | phase H_yp,u0 | ... | phase H_yp,um |
# +---------------+ +---------------+
#
# Several operations are available that change this layout.
#
# * Omitting: either the magnitude or the phase plots can be omitted
# using the plot_magnitude and plot_phase keywords.
#
# * Overlay: inputs and/or outputs can be combined onto a single set of
# axes using the overlay_inputs and overlay_outputs keywords. This
# basically collapses data along either the rows or columns, and a
# legend is generated.
#
# Decide on the maximum number of inputs and outputs
ninputs, noutputs = 0, 0
for response in data: # TODO: make more pythonic/numpic
ninputs = max(ninputs, response.ninputs)
noutputs = max(noutputs, response.noutputs)
# Figure how how many rows and columns to use + offsets for inputs/outputs
if overlay_outputs and overlay_inputs:
nrows = plot_magnitude + plot_phase
ncols = 1
elif overlay_outputs:
nrows = plot_magnitude + plot_phase
ncols = ninputs
elif overlay_inputs:
nrows = (noutputs if plot_magnitude else 0) + \
(noutputs if plot_phase else 0)
ncols = 1
else:
nrows = (noutputs if plot_magnitude else 0) + \
(noutputs if plot_phase else 0)
ncols = ninputs
if ax is None:
# Set up default sharing of axis limits if not specified
for kw in ['share_magnitude', 'share_phase', 'share_frequency']:
if kw not in kwargs or kwargs[kw] is None:
kwargs[kw] = config.defaults['freqplot.' + kw]
fig, ax_array = _process_ax_keyword(
ax, (nrows, ncols), squeeze=False, rcParams=rcParams, clear_text=True)
legend_loc, legend_map, show_legend = _process_legend_keywords(
kwargs, (nrows,ncols), 'center right')
# Get the values for sharing axes limits
share_magnitude = kwargs.pop('share_magnitude', None)
share_phase = kwargs.pop('share_phase', None)
share_frequency = kwargs.pop('share_frequency', None)
# Set up axes variables for easier access below
if plot_magnitude and not plot_phase:
mag_map = np.empty((noutputs, ninputs), dtype=tuple)
for i in range(noutputs):
for j in range(ninputs):
if overlay_outputs and overlay_inputs:
mag_map[i, j] = (0, 0)
elif overlay_outputs:
mag_map[i, j] = (0, j)
elif overlay_inputs:
mag_map[i, j] = (i, 0)
else:
mag_map[i, j] = (i, j)
phase_map = np.full((noutputs, ninputs), None)
share_phase = False
elif plot_phase and not plot_magnitude:
phase_map = np.empty((noutputs, ninputs), dtype=tuple)
for i in range(noutputs):
for j in range(ninputs):
if overlay_outputs and overlay_inputs:
phase_map[i, j] = (0, 0)
elif overlay_outputs:
phase_map[i, j] = (0, j)
elif overlay_inputs:
phase_map[i, j] = (i, 0)
else:
phase_map[i, j] = (i, j)
mag_map = np.full((noutputs, ninputs), None)
share_magnitude = False
else:
mag_map = np.empty((noutputs, ninputs), dtype=tuple)
phase_map = np.empty((noutputs, ninputs), dtype=tuple)
for i in range(noutputs):
for j in range(ninputs):
if overlay_outputs and overlay_inputs:
mag_map[i, j] = (0, 0)
phase_map[i, j] = (1, 0)
elif overlay_outputs:
mag_map[i, j] = (0, j)
phase_map[i, j] = (1, j)
elif overlay_inputs:
mag_map[i, j] = (i*2, 0)
phase_map[i, j] = (i*2 + 1, 0)
else:
mag_map[i, j] = (i*2, j)
phase_map[i, j] = (i*2 + 1, j)
# Identity map needed for setting up shared axes
ax_map = np.empty((nrows, ncols), dtype=tuple)
for i, j in itertools.product(range(nrows), range(ncols)):
ax_map[i, j] = (i, j)
#
# Set up axes limit sharing
#
# This code uses the share_magnitude, share_phase, and share_frequency
# keywords to decide which axes have shared limits and what ticklabels
# to include. The sharing code needs to come before the plots are
# generated, but additional code for removing tick labels needs to come
# *during* and *after* the plots are generated (see below).
#
# Note: if the various share_* keywords are None then a previous set of
# axes are available and no updates should be made.
#
# Utility function to turn on sharing
def _share_axes(ref, share_map, axis):
ref_ax = ax_array[ref]
for index in np.nditer(share_map, flags=["refs_ok"]):
if index.item() == ref:
continue
if axis == 'x':
ax_array[index.item()].sharex(ref_ax)
elif axis == 'y':
ax_array[index.item()].sharey(ref_ax)
else:
raise ValueError("axis must be 'x' or 'y'")
# Process magnitude, phase, and frequency axes
for name, value, map, axis in zip(
['share_magnitude', 'share_phase', 'share_frequency'],
[ share_magnitude, share_phase, share_frequency],
[ mag_map, phase_map, ax_map],
[ 'y', 'y', 'x']):
if value in [True, 'all']:
_share_axes(map[0 if axis == 'y' else -1, 0], map, axis)
elif axis == 'y' and value in ['row']:
for i in range(noutputs if not overlay_outputs else 1):
_share_axes(map[i, 0], map[i], 'y')
elif axis == 'x' and value in ['col']:
for j in range(ncols):
_share_axes(map[-1, j], map[:, j], 'x')
elif value in [False, 'none']:
# TODO: turn off any sharing that is on
pass
elif value is not None:
raise ValueError(
f"unknown value for `{name}`: '{value}'")
#
# Plot the data
#
# The mag_map and phase_map arrays have the indices axes needed for
# making the plots. Labels are used on each axes for later creation of
# legends. The generic labels if of the form:
#
# To output label, From input label, system name
#
# The input and output labels are omitted if overlay_inputs or
# overlay_outputs is False, respectively. The system name is always
# included, since multiple calls to plot() will require a legend that
# distinguishes which system signals are plotted. The system name is
# stripped off later (in the legend-handling code) if it is not needed.
#
# Note: if we are building on top of an existing plot, tick labels
# should be preserved from the existing axes. For log scale axes the
# tick labels seem to appear no matter what => we have to detect if
# they are present at the start and, it not, remove them after calling
# loglog or semilogx.
#
# Create a list of lines for the output
out = np.empty((nrows, ncols), dtype=object)
for i in range(nrows):
for j in range(ncols):
out[i, j] = [] # unique list in each element
# Process label keyword
line_labels = _process_line_labels(label, len(data), ninputs, noutputs)
# Utility function for creating line label
def _make_line_label(response, output_index, input_index):
label = "" # start with an empty label
# Add the output name if it won't appear as an axes label
if noutputs > 1 and overlay_outputs:
label += response.output_labels[output_index]
# Add the input name if it won't appear as a column label
if ninputs > 1 and overlay_inputs:
label += ", " if label != "" else ""
label += response.input_labels[input_index]
# Add the system name (will strip off later if redundant)
label += ", " if label != "" else ""
label += f"{response.sysname}"
return label
for index, response in enumerate(data):
# Get the (pre-processed) data in fully indexed form
mag = mag_data[index].reshape((noutputs, ninputs, -1))
phase = phase_data[index].reshape((noutputs, ninputs, -1))
omega_sys, sysname = omega_data[index], response.sysname
for i, j in itertools.product(range(noutputs), range(ninputs)):
# Get the axes to use for magnitude and phase
ax_mag = ax_array[mag_map[i, j]]
ax_phase = ax_array[phase_map[i, j]]
# Get the frequencies and convert to Hz, if needed
omega_plot = omega_sys / (2 * math.pi) if Hz else omega_sys
if response.isdtime(strict=True):
nyq_freq = (0.5/response.dt) if Hz else (math.pi/response.dt)
# Save the magnitude and phase to plot
mag_plot = 20 * np.log10(mag[i, j]) if dB else mag[i, j]
phase_plot = phase[i, j] * 180. / math.pi if deg else phase[i, j]
# Generate a label
if line_labels is None:
label = _make_line_label(response, i, j)
else:
label = line_labels[index, i, j]
# Magnitude
if plot_magnitude:
pltfcn = ax_mag.semilogx if dB else ax_mag.loglog
# Plot the main data
lines = pltfcn(
omega_plot, mag_plot, *fmt, label=label, **kwargs)
out[mag_map[i, j]] += lines
# Save the information needed for the Nyquist line
if response.isdtime(strict=True):
ax_mag.axvline(
nyq_freq, color=lines[0].get_color(), linestyle='--',
label='_nyq_mag_' + sysname)
# Add a grid to the plot
ax_mag.grid(grid and not display_margins, which='both')
# Phase
if plot_phase:
lines = ax_phase.semilogx(
omega_plot, phase_plot, *fmt, label=label, **kwargs)
out[phase_map[i, j]] += lines
# Save the information needed for the Nyquist line
if response.isdtime(strict=True):
ax_phase.axvline(
nyq_freq, color=lines[0].get_color(), linestyle='--',
label='_nyq_phase_' + sysname)
# Add a grid to the plot
ax_phase.grid(grid and not display_margins, which='both')
#
# Display gain and phase margins (SISO only)
#
if display_margins:
if ninputs > 1 or noutputs > 1:
raise NotImplementedError(
"margins are not available for MIMO systems")
# Compute stability margins for the system
margins = stability_margins(response, method=margins_method)
gm, pm, Wcg, Wcp = (margins[i] for i in [0, 1, 3, 4])
# Figure out sign of the phase at the first gain crossing
# (needed if phase_wrap is True)
phase_at_cp = phase[
0, 0, (np.abs(omega_data[0] - Wcp)).argmin()]
if phase_at_cp >= 0.:
phase_limit = 180.
else:
phase_limit = -180.
if Hz:
Wcg, Wcp = Wcg/(2*math.pi), Wcp/(2*math.pi)
# Draw lines at gain and phase limits
if plot_magnitude:
ax_mag.axhline(y=0 if dB else 1, color='k', linestyle=':',
zorder=-20)
if plot_phase:
ax_phase.axhline(y=phase_limit if deg else
math.radians(phase_limit),
color='k', linestyle=':', zorder=-20)
# Annotate the phase margin (if it exists)
if plot_phase and pm != float('inf') and Wcp != float('nan'):
# Draw dotted lines marking the gain crossover frequencies
if plot_magnitude:
ax_mag.axvline(Wcp, color='k', linestyle=':', zorder=-30)
ax_phase.axvline(Wcp, color='k', linestyle=':', zorder=-30)
# Draw solid segments indicating the margins
if deg:
ax_phase.semilogx(
[Wcp, Wcp], [phase_limit + pm, phase_limit],
color='k', zorder=-20)
else:
ax_phase.semilogx(
[Wcp, Wcp], [math.radians(phase_limit) +
math.radians(pm),
math.radians(phase_limit)],
color='k', zorder=-20)
# Annotate the gain margin (if it exists)
if plot_magnitude and gm != float('inf') and \
Wcg != float('nan'):
# Draw dotted lines marking the phase crossover frequencies
ax_mag.axvline(Wcg, color='k', linestyle=':', zorder=-30)
if plot_phase:
ax_phase.axvline(Wcg, color='k', linestyle=':', zorder=-30)
# Draw solid segments indicating the margins
if dB:
ax_mag.semilogx(
[Wcg, Wcg], [0, -20*np.log10(gm)],
color='k', zorder=-20)
else:
ax_mag.loglog(
[Wcg, Wcg], [1., 1./gm], color='k', zorder=-20)
if display_margins == 'overlay':
# TODO: figure out how to handle case of multiple lines
# Put the margin information in the lower left corner
if plot_magnitude:
ax_mag.text(
0.04, 0.06,
'G.M.: %.2f %s\nFreq: %.2f %s' %
(20*np.log10(gm) if dB else gm,
'dB ' if dB else '',
Wcg, 'Hz' if Hz else 'rad/s'),
horizontalalignment='left',
verticalalignment='bottom',
transform=ax_mag.transAxes,
fontsize=8 if int(mpl.__version__[0]) == 1 else 6)
if plot_phase:
ax_phase.text(
0.04, 0.06,
'P.M.: %.2f %s\nFreq: %.2f %s' %
(pm if deg else math.radians(pm),
'deg' if deg else 'rad',
Wcp, 'Hz' if Hz else 'rad/s'),
horizontalalignment='left',
verticalalignment='bottom',
transform=ax_phase.transAxes,
fontsize=8 if int(mpl.__version__[0]) == 1 else 6)
else:
# Put the title underneath the suptitle (one line per system)
ax = ax_mag if ax_mag else ax_phase
axes_title = ax.get_title()
if axes_title is not None and axes_title != "":
axes_title += "\n"
with plt.rc_context(rcParams):
ax.set_title(
axes_title + f"{sysname}: "
"Gm = %.2f %s(at %.2f %s), "
"Pm = %.2f %s (at %.2f %s)" %
(20*np.log10(gm) if dB else gm,
'dB ' if dB else '',
Wcg, 'Hz' if Hz else 'rad/s',
pm if deg else math.radians(pm),
'deg' if deg else 'rad',
Wcp, 'Hz' if Hz else 'rad/s'))
#
# Finishing handling axes limit sharing
#
# This code handles labels on Bode plots and also removes tick labels
# on shared axes. It needs to come *after* the plots are generated,
# in order to handle two things:
#
# * manually generated labels and grids need to reflect the limits for
# shared axes, which we don't know until we have plotted everything;
#
# * the loglog and semilog functions regenerate the labels (not quite
# sure why, since using sharex and sharey in subplots does not have
# this behavior).
#
# Note: as before, if the various share_* keywords are None then a
# previous set of axes are available and no updates are made. (TODO: true?)
#
for i in range(noutputs):
for j in range(ninputs):
# Utility function to generate phase labels
def gen_zero_centered_series(val_min, val_max, period):
v1 = np.ceil(val_min / period - 0.2)
v2 = np.floor(val_max / period + 0.2)
return np.arange(v1, v2 + 1) * period
# Label the phase axes using multiples of 45 degrees
if plot_phase:
ax_phase = ax_array[phase_map[i, j]]
# Set the labels
if deg:
ylim = ax_phase.get_ylim()
num = np.floor((ylim[1] - ylim[0]) / 45)
factor = max(1, np.round(num / (32 / nrows)) * 2)
ax_phase.set_yticks(gen_zero_centered_series(
ylim[0], ylim[1], 45 * factor))
ax_phase.set_yticks(gen_zero_centered_series(
ylim[0], ylim[1], 15 * factor), minor=True)
else:
ylim = ax_phase.get_ylim()
num = np.ceil((ylim[1] - ylim[0]) / (math.pi/4))
factor = max(1, np.round(num / (36 / nrows)) * 2)
ax_phase.set_yticks(gen_zero_centered_series(
ylim[0], ylim[1], math.pi / 4. * factor))
ax_phase.set_yticks(gen_zero_centered_series(
ylim[0], ylim[1], math.pi / 12. * factor), minor=True)
# Turn off y tick labels for shared axes
for i in range(0, noutputs):
for j in range(1, ncols):
if share_magnitude in [True, 'all', 'row']:
ax_array[mag_map[i, j]].tick_params(labelleft=False)
if share_phase in [True, 'all', 'row']:
ax_array[phase_map[i, j]].tick_params(labelleft=False)
# Turn off x tick labels for shared axes
for i in range(0, nrows-1):
for j in range(0, ncols):
if share_frequency in [True, 'all', 'col']:
ax_array[i, j].tick_params(labelbottom=False)
# If specific omega_limits were given, use them
if omega_limits is not None:
for i, j in itertools.product(range(nrows), range(ncols)):
ax_array[i, j].set_xlim(omega_limits)
#
# Label the axes (including header labels)
#
# Once the data are plotted, we label the axes. The horizontal axes is
# always frequency and this is labeled only on the bottom most row. The
# vertical axes can consist either of a single signal or a combination
# of signals (when overlay_inputs or overlay_outputs is True)
#
# Input/output signals are give at the top of columns and left of rows
# when these are individually plotted.
#
# Label the columns (do this first to get row labels in the right spot)
for j in range(ncols):
# If we have more than one column, label the individual responses
if (noutputs > 1 and not overlay_outputs or ninputs > 1) \
and not overlay_inputs:
with plt.rc_context(rcParams):
ax_array[0, j].set_title(f"From {data[0].input_labels[j]}")
# Label the frequency axis
ax_array[-1, j].set_xlabel(
freq_label.format(units="Hz" if Hz else "rad/s"))
# Label the rows
for i in range(noutputs if not overlay_outputs else 1):
if plot_magnitude:
ax_mag = ax_array[mag_map[i, 0]]
ax_mag.set_ylabel(magnitude_label)
if plot_phase:
ax_phase = ax_array[phase_map[i, 0]]
ax_phase.set_ylabel(phase_label)
if (noutputs > 1 or ninputs > 1) and not overlay_outputs:
if plot_magnitude and plot_phase:
# Get existing ylabel for left column and add a blank line
ax_mag.set_ylabel("\n" + ax_mag.get_ylabel())
ax_phase.set_ylabel("\n" + ax_phase.get_ylabel())
# Find the midpoint between the row axes (+ tight_layout)
_, ypos = _find_axes_center(fig, [ax_mag, ax_phase])
# Get the bounding box including the labels
inv_transform = fig.transFigure.inverted()
mag_bbox = inv_transform.transform(
ax_mag.get_tightbbox(fig.canvas.get_renderer()))
# Figure out location for the text (center left in figure frame)
xpos = mag_bbox[0, 0] # left edge
# Put a centered label as text outside the box
fig.text(
0.8 * xpos, ypos, f"To {data[0].output_labels[i]}\n",
rotation=90, ha='left', va='center',
fontsize=rcParams['axes.titlesize'])
else:
# Only a single axes => add label to the left
ax_array[i, 0].set_ylabel(
f"To {data[0].output_labels[i]}\n" +
ax_array[i, 0].get_ylabel())
#
# Update the plot title (= figure suptitle)
#
# If plots are built up by multiple calls to plot() and the title is
# not given, then the title is updated to provide a list of unique text
# items in each successive title. For data generated by the frequency
# response function this will generate a common prefix followed by a
# list of systems (e.g., "Step response for sys[1], sys[2]").
#
# Set the initial title for the data (unique system names, preserving order)
seen = set()
sysnames = [response.sysname for response in data \
if not (response.sysname in seen or seen.add(response.sysname))]
if ax is None and title is None:
if data[0].title is None:
title = "Bode plot for " + ", ".join(sysnames)
else:
# Allow data to set the title (used by gangof4)
title = data[0].title
_update_plot_title(title, fig, rcParams=rcParams, frame=title_frame)
elif ax is None:
_update_plot_title(