# 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 numpy as np import matplotlib as mpl import matplotlib.pyplot as plt import math import warnings import itertools from os.path import commonprefix from .ctrlutil import unwrap from .bdalg import feedback from .margins import stability_margins from .exception import ControlMIMONotImplemented from .statesp import StateSpace from .lti import LTI, frequency_response, _process_frequency_response from .xferfcn import TransferFunction from .frdata import FrequencyResponseData from .timeplot import _make_legend_labels from . import config __all__ = ['bode_plot', 'NyquistResponseData', 'nyquist_response', 'nyquist_plot', 'singular_values_response', 'singular_values_plot', 'gangof4_plot', 'gangof4_response', 'bode', 'nyquist', 'gangof4'] # Default font dictionary _freqplot_rcParams = mpl.rcParams.copy() _freqplot_rcParams.update({ 'axes.labelsize': 'small', 'axes.titlesize': 'small', 'figure.titlesize': 'medium', 'legend.fontsize': 'x-small', 'xtick.labelsize': 'small', 'ytick.labelsize': 'small', }) # Default values for module parameter variables _freqplot_defaults = { 'freqplot.rcParams': _freqplot_rcParams, '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 [%s]", 'freqplot.share_magnitude': 'row', 'freqplot.share_phase': 'row', 'freqplot.share_frequency': 'col', } # # 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, display_margins=None, margins_method='best', legend_map=None, legend_loc=None, sharex=None, sharey=None, title=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 List of frequencies in rad/sec over 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. margins_method : str, optional Method to use in computing margins (see :func:`stability_margins`). **kwargs : :func:`matplotlib.pyplot.plot` keyword properties, optional Additional keywords passed to `matplotlib` to specify line properties. Returns ------- lines : array of Line2D Array of 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. Other Parameters ---------------- grid : bool If True, plot grid lines on gain and phase plots. Default is set by `config.defaults['freqplot.grid']`. initial_phase : float 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. 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. omega_num : int Number of samples to use for the frequeny range. Defaults to config.defaults['freqplot.number_of_samples']. Ignore if data is not a list of systems. 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. rcParams : dict Override the default parameters used for generating plots. Default is set by config.default['freqplot.rcParams']. 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']. The default values for Bode plot configuration parameters can be reset using the `config.defaults` dictionary, with module name 'bode'. Notes ----- 1. Starting with python-control version 0.10, `bode_plot`returns an array of lines 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. 2. 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. 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) freqplot_rcParams = config._get_param( 'freqplot', 'rcParams', 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 all([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 values of mag, phase, omega is deprecated; " "use frequency_response()", DeprecationWarning) 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 # See if we can use the current figure axes fig = plt.gcf() # get current figure (or create new one) if ax is None and plt.get_fignums(): ax = fig.get_axes() if len(ax) == nrows * ncols: # Assume that the shape is right (no easy way to infer this) ax = np.array(ax).reshape(nrows, ncols) # Clear out any old text from the current figure for text in fig.texts: text.set_visible(False) # turn off the text del text # get rid of it completely elif len(ax) != 0: # Need to generate a new figure fig, ax = plt.figure(), None else: # Blank figure, just need to recreate axes ax = None # Create new axes, if needed, and customize them if ax is None: with plt.rc_context(_freqplot_rcParams): ax_array = fig.subplots(nrows, ncols, squeeze=False) fig.set_layout_engine('tight') fig.align_labels() # 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] else: # Make sure the axes are the right shape if ax.shape != (nrows, ncols): raise ValueError( "specified axes are not the right shape; " f"got {ax.shape} but expecting ({nrows}, {ncols})") ax_array = ax fig = ax_array[0, 0].figure # just in case this is not gcf() # 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 off 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 # 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 label = _make_line_label(response, 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) mag_ylim = ax_mag.get_ylim() if plot_phase: ax_phase.axhline(y=phase_limit if deg else math.radians(phase_limit), color='k', linestyle=':', zorder=-20) phase_ylim = ax_phase.get_ylim() # 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(_freqplot_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 phase 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 limts 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) # # 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 title is None: if data[0].title is None: title = "Bode plot for " + ", ".join(sysnames) else: title = data[0].title if fig is not None and isinstance(title, str): # Get the current title, if it exists old_title = None if fig._suptitle is None else fig._suptitle._text new_title = title if old_title is not None: # Find the common part of the titles common_prefix = commonprefix([old_title, new_title]) # Back up to the last space last_space = common_prefix.rfind(' ') if last_space > 0: common_prefix = common_prefix[:last_space] common_len = len(common_prefix) # Add the new part of the title (usually the system name) if old_title[common_len:] != new_title[common_len:]: separator = ',' if len(common_prefix) > 0 else ';' new_title = old_title + separator + new_title[common_len:] # Add the title with plt.rc_context(freqplot_rcParams): fig.suptitle(new_title) # # 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(_freqplot_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 % ("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()) # TODO: remove? # Redraw the figure to get the proper locations for everything # fig.tight_layout() # 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())) phase_bbox = inv_transform.transform( ax_phase.get_tightbbox(fig.canvas.get_renderer())) # Get the axes limits without labels for use in the y position mag_bot = inv_transform.transform( ax_mag.transAxes.transform((0, 0)))[1] phase_top = inv_transform.transform( ax_phase.transAxes.transform((0, 1)))[1] # Figure out location for the text (center left in figure frame) xpos = mag_bbox[0, 0] # left edge ypos = (mag_bot + phase_top) / 2 # centered between axes # 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=_freqplot_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()) # # Create legends # # Legends can be placed manually by passing a legend_map array that # matches the shape of the suplots, with each item being a string # indicating the location of the legend for that axes (or None for no # legend). # # If no legend spec is passed, a minimal number of legends are used so # that each line in each axis can be uniquely identified. The details # depends on the various plotting parameters, but the general rule is # to place legends in the top row and right column. # # Because plots can be built up by multiple calls to plot(), the legend # strings are created from the line labels manually. Thus an initial # call to plot() may not generate any legends (eg, if no signals are # overlaid), but subsequent calls to plot() will need a legend for each # different response (system). # # Figure out where to put legends if legend_map is None: legend_map = np.full(ax_array.shape, None, dtype=object) if legend_loc == None: legend_loc = 'center right' # TODO: add in additional processing later # Put legend in the upper right legend_map[0, -1] = legend_loc # Create axis legends for i in range(nrows): for j in range(ncols): ax = ax_array[i, j] # Get the labels to use, removing common strings lines = [line for line in ax.get_lines() if line.get_label()[0] != '_'] labels = _make_legend_labels([line.get_label() for line in lines]) # Generate the label, if needed if len(labels) > 1 and legend_map[i, j] != None: with plt.rc_context(freqplot_rcParams): ax.legend(lines, labels, loc=legend_map[i, j]) # # Legacy return pocessing # if plot is True: # legacy usage; remove in future release # 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 return out # # Nyquist plot # # Default values for module parameter variables _nyquist_defaults = { 'nyquist.primary_style': ['-', '-.'], # style for primary curve 'nyquist.mirror_style': ['--', ':'], # style for mirror curve 'nyquist.arrows': 2, # number of arrows around curve 'nyquist.arrow_size': 8, # pixel size for arrows 'nyquist.encirclement_threshold': 0.05, # warning threshold 'nyquist.indent_radius': 1e-4, # indentation radius 'nyquist.indent_direction': 'right', # indentation direction 'nyquist.indent_points': 50, # number of points to insert 'nyquist.max_curve_magnitude': 20, # clip large values 'nyquist.max_curve_offset': 0.02, # offset of primary/mirror 'nyquist.start_marker': 'o', # marker at start of curve 'nyquist.start_marker_size': 4, # size of the maker } class NyquistResponseData: """Nyquist response data object. Nyquist contour analysis allows the stability and robustness of a closed loop linear system to be evaluated using the open loop response of the loop transfer function. The NyquistResponseData class is used by the :func:`~control.nyquist_response` function to return the response of a linear system along the Nyquist 'D' contour. The response object can be used to obtain information about the Nyquist response or to generate a Nyquist plot. Attributes ---------- count : integer Number of encirclements of the -1 point by the Nyquist curve for a system evaluated along the Nyquist contour. contour : complex array The Nyquist 'D' contour, with appropriate indendtations to avoid open loop poles and zeros near/on the imaginary axis. response : complex array The value of the linear system under study along the Nyquist contour. dt : None or float The system timebase. sysname : str The name of the system being analyzed. return_contour: bool If true, when the object is accessed as an iterable return two elements": `count` (number of encirlements) and `contour`. If false (default), then return only `count`. """ def __init__( self, count, contour, response, dt, sysname=None, return_contour=False): self.count = count self.contour = contour self.response = response self.dt = dt self.sysname = sysname self.return_contour = return_contour # Implement iter to allow assigning to a tuple def __iter__(self): if self.return_contour: return iter((self.count, self.contour)) else: return iter((self.count, )) # Implement (thin) getitem to allow access via legacy indexing def __getitem__(self, index): return list(self.__iter__())[index] # Implement (thin) len to emulate legacy testing interface def __len__(self): return 2 if self.return_contour else 1 def plot(self, *args, **kwargs): return nyquist_plot(self, *args, **kwargs) class NyquistResponseList(list): def plot(self, *args, **kwargs): return nyquist_plot(self, *args, **kwargs) def nyquist_response( sysdata, omega=None, plot=None, omega_limits=None, omega_num=None, return_contour=False, warn_encirclements=True, warn_nyquist=True, check_kwargs=True, **kwargs): """Nyquist response for a system. Computes a Nyquist contour for the system over a (optional) frequency range and evaluates the number of net encirclements. The curve is computed by evaluating the Nyqist segment along the positive imaginary axis, with a mirror image generated to reflect the negative imaginary axis. Poles on or near the imaginary axis are avoided using a small indentation. The portion of the Nyquist contour at infinity is not explicitly computed (since it maps to a constant value for any system with a proper transfer function). Parameters ---------- sysdata : LTI or list of LTI List of linear input/output systems (single system is OK). Nyquist curves for each system are plotted on the same graph. omega : array_like, optional Set of frequencies to be evaluated, in rad/sec. omega_limits : array_like of two values, optional Limits to the range of frequencies. Ignored if omega is provided, and auto-generated if omitted. omega_num : int, optional Number of frequency samples to plot. Defaults to config.defaults['freqplot.number_of_samples']. Returns ------- responses : list of :class:`~control.NyquistResponseData` For each system, a Nyquist response data object is returned. If `sysdata` is a single system, a single elemeent is returned (not a list). For each response, the following information is available: response.count : int Number of encirclements of the point -1 by the Nyquist curve. If multiple systems are given, an array of counts is returned. response.contour : ndarray The contour used to create the primary Nyquist curve segment. To obtain the Nyquist curve values, evaluate system(s) along contour. Other Parameters ---------------- encirclement_threshold : float, optional Define the threshold for generating a warning if the number of net encirclements is a non-integer value. Default value is 0.05 and can be set using config.defaults['nyquist.encirclement_threshold']. indent_direction : str, optional For poles on the imaginary axis, set the direction of indentation to be 'right' (default), 'left', or 'none'. indent_points : int, optional Number of points to insert in the Nyquist contour around poles that are at or near the imaginary axis. indent_radius : float, optional Amount to indent the Nyquist contour around poles on or near the imaginary axis. Portions of the Nyquist plot corresponding to indented portions of the contour are plotted using a different line style. warn_nyquist : bool, optional If set to 'False', turn off warnings about frequencies above Nyquist. warn_encirclements : bool, optional If set to 'False', turn off warnings about number of encirclements not meeting the Nyquist criterion. Notes ----- 1. If a discrete time model is given, the frequency response is computed 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. 2. If a continuous-time system contains poles on or near the imaginary axis, a small indentation will be used to avoid the pole. The radius of the indentation is given by `indent_radius` and it is taken to the right of stable poles and the left of unstable poles. If a pole is exactly on the imaginary axis, the `indent_direction` parameter can be used to set the direction of indentation. Setting `indent_direction` to `none` will turn off indentation. If `return_contour` is True, the exact contour used for evaluation is returned. 3. For those portions of the Nyquist plot in which the contour is indented to avoid poles, resuling in a scaling of the Nyquist plot, the line styles are according to the settings of the `primary_style` and `mirror_style` keywords. By default the scaled portions of the primary curve use a dotted line style and the scaled portion of the mirror image use a dashdot line style. 4. If the legacy keyword `return_contour` is specified as True, the response object can be iterated over to return `count, contour`. This behavior is deprecated and will be removed in a future release. Examples -------- >>> G = ct.zpk([], [-1, -2, -3], gain=100) >>> response = ct.nyquist_response(G) >>> count = response.count >>> lines = response.plot() """ # Get values for params omega_num_given = omega_num is not None omega_num = config._get_param('freqplot', 'number_of_samples', omega_num) indent_radius = config._get_param( 'nyquist', 'indent_radius', kwargs, _nyquist_defaults, pop=True) encirclement_threshold = config._get_param( 'nyquist', 'encirclement_threshold', kwargs, _nyquist_defaults, pop=True) indent_direction = config._get_param( 'nyquist', 'indent_direction', kwargs, _nyquist_defaults, pop=True) indent_points = config._get_param( 'nyquist', 'indent_points', kwargs, _nyquist_defaults, pop=True) if check_kwargs and kwargs: raise TypeError("unrecognized keywords: ", str(kwargs)) # Convert the first argument to a list syslist = sysdata if isinstance(sysdata, (list, tuple)) else [sysdata] # Determine the range of frequencies to use, based on args/features omega, omega_range_given = _determine_omega_vector( syslist, omega, omega_limits, omega_num, feature_periphery_decades=2) # If omega was not specified explicitly, start at omega = 0 if not omega_range_given: if omega_num_given: # Just reset the starting point omega[0] = 0.0 else: # Insert points between the origin and the first frequency point omega = np.concatenate(( np.linspace(0, omega[0], indent_points), omega[1:])) # Go through each system and keep track of the results responses = [] for idx, sys in enumerate(syslist): if not sys.issiso(): # TODO: Add MIMO nyquist plots. raise ControlMIMONotImplemented( "Nyquist plot currently only supports SISO systems.") # Figure out the frequency range omega_sys = np.asarray(omega) # Determine the contour used to evaluate the Nyquist curve if sys.isdtime(strict=True): # Restrict frequencies for discrete-time systems nyq_freq = math.pi / sys.dt if not omega_range_given: # limit up to and including Nyquist frequency omega_sys = np.hstack(( omega_sys[omega_sys < nyq_freq], nyq_freq)) # Issue a warning if we are sampling above Nyquist if np.any(omega_sys * sys.dt > np.pi) and warn_nyquist: warnings.warn("evaluation above Nyquist frequency") # do indentations in s-plane where it is more convenient splane_contour = 1j * omega_sys # Bend the contour around any poles on/near the imaginary axis if isinstance(sys, (StateSpace, TransferFunction)) \ and indent_direction != 'none': if sys.isctime(): splane_poles = sys.poles() splane_cl_poles = sys.feedback().poles() else: # map z-plane poles to s-plane. We ignore any at the origin # to avoid numerical warnings because we know we # don't need to indent for them zplane_poles = sys.poles() zplane_poles = zplane_poles[~np.isclose(abs(zplane_poles), 0.)] splane_poles = np.log(zplane_poles) / sys.dt zplane_cl_poles = sys.feedback().poles() # eliminate z-plane poles at the origin to avoid warnings zplane_cl_poles = zplane_cl_poles[ ~np.isclose(abs(zplane_cl_poles), 0.)] splane_cl_poles = np.log(zplane_cl_poles) / sys.dt # # Check to make sure indent radius is small enough # # If there is a closed loop pole that is near the imaginary axis # at a point that is near an open loop pole, it is possible that # indentation might skip or create an extraneous encirclement. # We check for that situation here and generate a warning if that # could happen. # for p_cl in splane_cl_poles: # See if any closed loop poles are near the imaginary axis if abs(p_cl.real) <= indent_radius: # See if any open loop poles are close to closed loop poles if len(splane_poles) > 0: p_ol = splane_poles[ (np.abs(splane_poles - p_cl)).argmin()] if abs(p_ol - p_cl) <= indent_radius and \ warn_encirclements: warnings.warn( "indented contour may miss closed loop pole; " "consider reducing indent_radius to below " f"{abs(p_ol - p_cl):5.2g}", stacklevel=2) # # See if we should add some frequency points near imaginary poles # for p in splane_poles: # See if we need to process this pole (skip if on the negative # imaginary axis or not near imaginary axis + user override) if p.imag < 0 or abs(p.real) > indent_radius or \ omega_range_given: continue # Find the frequencies before the pole frequency below_points = np.argwhere( splane_contour.imag - abs(p.imag) < -indent_radius) if below_points.size > 0: first_point = below_points[-1].item() start_freq = p.imag - indent_radius else: # Add the points starting at the beginning of the contour assert splane_contour[0] == 0 first_point = 0 start_freq = 0 # Find the frequencies after the pole frequency above_points = np.argwhere( splane_contour.imag - abs(p.imag) > indent_radius) last_point = above_points[0].item() # Add points for half/quarter circle around pole frequency # (these will get indented left or right below) splane_contour = np.concatenate(( splane_contour[0:first_point+1], (1j * np.linspace( start_freq, p.imag + indent_radius, indent_points)), splane_contour[last_point:])) # Indent points that are too close to a pole if len(splane_poles) > 0: # accomodate no splane poles if dtime sys for i, s in enumerate(splane_contour): # Find the nearest pole p = splane_poles[(np.abs(splane_poles - s)).argmin()] # See if we need to indent around it if abs(s - p) < indent_radius: # Figure out how much to offset (simple trigonometry) offset = np.sqrt( indent_radius ** 2 - (s - p).imag ** 2) \ - (s - p).real # Figure out which way to offset the contour point if p.real < 0 or (p.real == 0 and indent_direction == 'right'): # Indent to the right splane_contour[i] += offset elif p.real > 0 or (p.real == 0 and indent_direction == 'left'): # Indent to the left splane_contour[i] -= offset else: raise ValueError( "unknown value for indent_direction") # change contour to z-plane if necessary if sys.isctime(): contour = splane_contour else: contour = np.exp(splane_contour * sys.dt) # Compute the primary curve resp = sys(contour) # Compute CW encirclements of -1 by integrating the (unwrapped) angle phase = -unwrap(np.angle(resp + 1)) encirclements = np.sum(np.diff(phase)) / np.pi count = int(np.round(encirclements, 0)) # Let the user know if the count might not make sense if abs(encirclements - count) > encirclement_threshold and \ warn_encirclements: warnings.warn( "number of encirclements was a non-integer value; this can" " happen is contour is not closed, possibly based on a" " frequency range that does not include zero.") # # Make sure that the enciriclements match the Nyquist criterion # # If the user specifies the frequency points to use, it is possible # to miss enciriclements, so we check here to make sure that the # Nyquist criterion is actually satisfied. # if isinstance(sys, (StateSpace, TransferFunction)): # Count the number of open/closed loop RHP poles if sys.isctime(): if indent_direction == 'right': P = (sys.poles().real > 0).sum() else: P = (sys.poles().real >= 0).sum() Z = (sys.feedback().poles().real >= 0).sum() else: if indent_direction == 'right': P = (np.abs(sys.poles()) > 1).sum() else: P = (np.abs(sys.poles()) >= 1).sum() Z = (np.abs(sys.feedback().poles()) >= 1).sum() # Check to make sure the results make sense; warn if not if Z != count + P and warn_encirclements: warnings.warn( "number of encirclements does not match Nyquist criterion;" " check frequency range and indent radius/direction", UserWarning, stacklevel=2) elif indent_direction == 'none' and any(sys.poles().real == 0) and \ warn_encirclements: warnings.warn( "system has pure imaginary poles but indentation is" " turned off; results may be meaningless", RuntimeWarning, stacklevel=2) # Decide on system name sysname = sys.name if sys.name is not None else f"Unknown-{idx}" responses.append(NyquistResponseData( count, contour, resp, sys.dt, sysname=sysname, return_contour=return_contour)) if isinstance(sysdata, (list, tuple)): return NyquistResponseList(responses) else: return responses[0] def nyquist_plot( data, omega=None, plot=None, label_freq=0, color=None, return_contour=None, title=None, legend_loc='upper right', **kwargs): """Nyquist plot for a system. Generates a Nyquist plot for the system over a (optional) frequency range. The curve is computed by evaluating the Nyqist segment along the positive imaginary axis, with a mirror image generated to reflect the negative imaginary axis. Poles on or near the imaginary axis are avoided using a small indentation. The portion of the Nyquist contour at infinity is not explicitly computed (since it maps to a constant value for any system with a proper transfer function). Parameters ---------- data : list of LTI or NyquistResponseData List of linear input/output systems (single system is OK) or Nyquist ersponses (computed using :func:`~control.nyquist_response`). Nyquist curves for each system are plotted on the same graph. omega : array_like, optional Set of frequencies to be evaluated, in rad/sec. omega_limits : array_like of two values, optional Limits to the range of frequencies. Ignored if omega is provided, and auto-generated if omitted. omega_num : int, optional Number of frequency samples to plot. Defaults to config.defaults['freqplot.number_of_samples']. color : string, optional Used to specify the color of the line and arrowhead. return_contour : bool, optional If 'True', return the contour used to evaluate the Nyquist plot. **kwargs : :func:`matplotlib.pyplot.plot` keyword properties, optional Additional keywords (passed to `matplotlib`) Returns ------- lines : array of Line2D 2D array of Line2D objects for each line in the plot. The shape of the array is given by (nsys, 4) where nsys is the number of systems or Nyquist responses passed to the function. The second index specifies the segment type: * lines[idx, 0]: unscaled portion of the primary curve * lines[idx, 1]: scaled portion of the primary curve * lines[idx, 2]: unscaled portion of the mirror curve * lines[idx, 3]: scaled portion of the mirror curve Other Parameters ---------------- arrows : int or 1D/2D array of floats, optional Specify the number of arrows to plot on the Nyquist curve. If an integer is passed. that number of equally spaced arrows will be plotted on each of the primary segment and the mirror image. If a 1D array is passed, it should consist of a sorted list of floats between 0 and 1, indicating the location along the curve to plot an arrow. If a 2D array is passed, the first row will be used to specify arrow locations for the primary curve and the second row will be used for the mirror image. arrow_size : float, optional Arrowhead width and length (in display coordinates). Default value is 8 and can be set using config.defaults['nyquist.arrow_size']. arrow_style : matplotlib.patches.ArrowStyle, optional Define style used for Nyquist curve arrows (overrides `arrow_size`). encirclement_threshold : float, optional Define the threshold for generating a warning if the number of net encirclements is a non-integer value. Default value is 0.05 and can be set using config.defaults['nyquist.encirclement_threshold']. indent_direction : str, optional For poles on the imaginary axis, set the direction of indentation to be 'right' (default), 'left', or 'none'. indent_points : int, optional Number of points to insert in the Nyquist contour around poles that are at or near the imaginary axis. indent_radius : float, optional Amount to indent the Nyquist contour around poles on or near the imaginary axis. Portions of the Nyquist plot corresponding to indented portions of the contour are plotted using a different line style. label_freq : int, optiona Label every nth frequency on the plot. If not specified, no labels are generated. max_curve_magnitude : float, optional Restrict the maximum magnitude of the Nyquist plot to this value. Portions of the Nyquist plot whose magnitude is restricted are plotted using a different line style. max_curve_offset : float, optional When plotting scaled portion of the Nyquist plot, increase/decrease the magnitude by this fraction of the max_curve_magnitude to allow any overlaps between the primary and mirror curves to be avoided. mirror_style : [str, str] or False Linestyles for mirror image of the Nyquist curve. The first element is used for unscaled portions of the Nyquist curve, the second element is used for portions that are scaled (using max_curve_magnitude). If `False` then omit completely. Default linestyle (['--', ':']) is determined by config.defaults['nyquist.mirror_style']. 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. primary_style : [str, str], optional Linestyles for primary image of the Nyquist curve. The first element is used for unscaled portions of the Nyquist curve, the second element is used for portions that are scaled (using max_curve_magnitude). Default linestyle (['-', '-.']) is determined by config.defaults['nyquist.mirror_style']. start_marker : str, optional Matplotlib marker to use to mark the starting point of the Nyquist plot. Defaults value is 'o' and can be set using config.defaults['nyquist.start_marker']. start_marker_size : float, optional Start marker size (in display coordinates). Default value is 4 and can be set using config.defaults['nyquist.start_marker_size']. warn_nyquist : bool, optional If set to 'False', turn off warnings about frequencies above Nyquist. warn_encirclements : bool, optional If set to 'False', turn off warnings about number of encirclements not meeting the Nyquist criterion. Notes ----- 1. If a discrete time model is given, the frequency response is computed 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. 2. If a continuous-time system contains poles on or near the imaginary axis, a small indentation will be used to avoid the pole. The radius of the indentation is given by `indent_radius` and it is taken to the right of stable poles and the left of unstable poles. If a pole is exactly on the imaginary axis, the `indent_direction` parameter can be used to set the direction of indentation. Setting `indent_direction` to `none` will turn off indentation. If `return_contour` is True, the exact contour used for evaluation is returned. 3. For those portions of the Nyquist plot in which the contour is indented to avoid poles, resuling in a scaling of the Nyquist plot, the line styles are according to the settings of the `primary_style` and `mirror_style` keywords. By default the scaled portions of the primary curve use a dotted line style and the scaled portion of the mirror image use a dashdot line style. Examples -------- >>> G = ct.zpk([], [-1, -2, -3], gain=100) >>> out = ct.nyquist_plot(G) """ # # Keyword processing # # Keywords for the nyquist_plot function can either be keywords that # are unique to this function, keywords that are intended for use by # nyquist_response (if data is a list of systems), or keywords that # are intended for the plotting commands. # # We first pop off all keywords that are used directly by this # function. If data is a list of systems, when then pop off keywords # that correspond to nyquist_response() keywords. The remaining # keywords are passed to matplotlib (and will generate an error if # unrecognized). # # Get values for params (and pop from list to allow keyword use in plot) arrows = config._get_param( 'nyquist', 'arrows', kwargs, _nyquist_defaults, pop=True) arrow_size = config._get_param( 'nyquist', 'arrow_size', kwargs, _nyquist_defaults, pop=True) arrow_style = config._get_param('nyquist', 'arrow_style', kwargs, None) max_curve_magnitude = config._get_param( 'nyquist', 'max_curve_magnitude', kwargs, _nyquist_defaults, pop=True) max_curve_offset = config._get_param( 'nyquist', 'max_curve_offset', kwargs, _nyquist_defaults, pop=True) start_marker = config._get_param( 'nyquist', 'start_marker', kwargs, _nyquist_defaults, pop=True) start_marker_size = config._get_param( 'nyquist', 'start_marker_size', kwargs, _nyquist_defaults, pop=True) # Set line styles for the curves def _parse_linestyle(style_name, allow_false=False): style = config._get_param( 'nyquist', style_name, kwargs, _nyquist_defaults, pop=True) if isinstance(style, str): # Only one style provided, use the default for the other style = [style, _nyquist_defaults['nyquist.' + style_name][1]] warnings.warn( "use of a single string for linestyle will be deprecated " " in a future release", PendingDeprecationWarning) if (allow_false and style is False) or \ (isinstance(style, list) and len(style) == 2): return style else: raise ValueError(f"invalid '{style_name}': {style}") primary_style = _parse_linestyle('primary_style') mirror_style = _parse_linestyle('mirror_style', allow_false=True) # Parse the arrows keyword if not arrows: arrow_pos = [] elif isinstance(arrows, int): N = arrows # Space arrows out, starting midway along each "region" arrow_pos = np.linspace(0.5/N, 1 + 0.5/N, N, endpoint=False) elif isinstance(arrows, (list, np.ndarray)): arrow_pos = np.sort(np.atleast_1d(arrows)) else: raise ValueError("unknown or unsupported arrow location") # Set the arrow style if arrow_style is None: arrow_style = mpl.patches.ArrowStyle( 'simple', head_width=arrow_size, head_length=arrow_size) # If argument was a singleton, turn it into a tuple if not isinstance(data, (list, tuple)): data = [data] # If we are passed a list of systems, compute response first if all([isinstance( sys, (StateSpace, TransferFunction, FrequencyResponseData)) for sys in data]): # Get the response, popping off keywords used there nyquist_responses = nyquist_response( data, omega=omega, return_contour=return_contour, omega_limits=kwargs.pop('omega_limits', None), omega_num=kwargs.pop('omega_num', None), warn_encirclements=kwargs.pop('warn_encirclements', True), warn_nyquist=kwargs.pop('warn_nyquist', True), check_kwargs=False, **kwargs) else: nyquist_responses = data # Legacy return value processing if plot is not None or return_contour is not None: warnings.warn( "`nyquist_plot` return values of count[, contour] is deprecated; " "use nyquist_response()", DeprecationWarning) # Extract out the values that we will eventually return counts = [response.count for response in nyquist_responses] contours = [response.contour for response in nyquist_responses] if plot is False: # Make sure we used all of the keywrods if kwargs: raise TypeError("unrecognized keywords: ", str(kwargs)) if len(data) == 1: counts, contours = counts[0], contours[0] # Return counts and (optionally) the contour we used return (counts, contours) if return_contour else counts # Create a list of lines for the output out = np.empty(len(nyquist_responses), dtype=object) for i in range(out.shape[0]): out[i] = [] # unique list in each element for idx, response in enumerate(nyquist_responses): resp = response.response if response.dt in [0, None]: splane_contour = response.contour else: splane_contour = np.log(response.contour) / response.dt # Find the different portions of the curve (with scaled pts marked) reg_mask = np.logical_or( np.abs(resp) > max_curve_magnitude, splane_contour.real != 0) # reg_mask = np.logical_or( # np.abs(resp.real) > max_curve_magnitude, # np.abs(resp.imag) > max_curve_magnitude) scale_mask = ~reg_mask \ & np.concatenate((~reg_mask[1:], ~reg_mask[-1:])) \ & np.concatenate((~reg_mask[0:1], ~reg_mask[:-1])) # Rescale the points with large magnitude rescale = np.logical_and( reg_mask, abs(resp) > max_curve_magnitude) resp[rescale] *= max_curve_magnitude / abs(resp[rescale]) # Plot the regular portions of the curve (and grab the color) x_reg = np.ma.masked_where(reg_mask, resp.real) y_reg = np.ma.masked_where(reg_mask, resp.imag) p = plt.plot( x_reg, y_reg, primary_style[0], color=color, label=response.sysname, **kwargs) c = p[0].get_color() out[idx] += p # Figure out how much to offset the curve: the offset goes from # zero at the start of the scaled section to max_curve_offset as # we move along the curve curve_offset = _compute_curve_offset( resp, scale_mask, max_curve_offset) # Plot the scaled sections of the curve (changing linestyle) x_scl = np.ma.masked_where(scale_mask, resp.real) y_scl = np.ma.masked_where(scale_mask, resp.imag) if x_scl.count() >= 1 and y_scl.count() >= 1: out[idx] += plt.plot( x_scl * (1 + curve_offset), y_scl * (1 + curve_offset), primary_style[1], color=c, **kwargs) else: out[idx] += [None] # Plot the primary curve (invisible) for setting arrows x, y = resp.real.copy(), resp.imag.copy() x[reg_mask] *= (1 + curve_offset[reg_mask]) y[reg_mask] *= (1 + curve_offset[reg_mask]) p = plt.plot(x, y, linestyle='None', color=c) # Add arrows ax = plt.gca() _add_arrows_to_line2D( ax, p[0], arrow_pos, arrowstyle=arrow_style, dir=1) # Plot the mirror image if mirror_style is not False: # Plot the regular and scaled segments out[idx] += plt.plot( x_reg, -y_reg, mirror_style[0], color=c, **kwargs) if x_scl.count() >= 1 and y_scl.count() >= 1: out[idx] += plt.plot( x_scl * (1 - curve_offset), -y_scl * (1 - curve_offset), mirror_style[1], color=c, **kwargs) else: out[idx] += [None] # Add the arrows (on top of an invisible contour) x, y = resp.real.copy(), resp.imag.copy() x[reg_mask] *= (1 - curve_offset[reg_mask]) y[reg_mask] *= (1 - curve_offset[reg_mask]) p = plt.plot(x, -y, linestyle='None', color=c, **kwargs) _add_arrows_to_line2D( ax, p[0], arrow_pos, arrowstyle=arrow_style, dir=-1) else: out[idx] += [None, None] # Mark the start of the curve if start_marker: plt.plot(resp[0].real, resp[0].imag, start_marker, color=c, markersize=start_marker_size) # Mark the -1 point plt.plot([-1], [0], 'r+') # Label the frequencies of the points if label_freq: ind = slice(None, None, label_freq) omega_sys = np.imag(splane_contour[np.real(splane_contour) == 0]) for xpt, ypt, omegapt in zip(x[ind], y[ind], omega_sys[ind]): # Convert to Hz f = omegapt / (2 * np.pi) # Factor out multiples of 1000 and limit the # result to the range [-8, 8]. pow1000 = max(min(get_pow1000(f), 8), -8) # Get the SI prefix. prefix = gen_prefix(pow1000) # Apply the text. (Use a space before the text to # prevent overlap with the data.) # # np.round() is used because 0.99... appears # instead of 1.0, and this would otherwise be # truncated to 0. plt.text(xpt, ypt, ' ' + str(int(np.round(f / 1000 ** pow1000, 0))) + ' ' + prefix + 'Hz') # Label the axes fig, ax = plt.gcf(), plt.gca() ax.set_xlabel("Real axis") ax.set_ylabel("Imaginary axis") ax.grid(color="lightgray") # List of systems that are included in this plot lines, labels = _get_line_labels(ax) # Add legend if there is more than one system plotted if len(labels) > 1: ax.legend(lines, labels, loc=legend_loc) # Add the title if title is None: title = "Nyquist plot for " + ", ".join(labels) fig.suptitle(title) # Legacy return pocessing if plot is True or return_contour is not None: if len(data) == 1: counts, contours = counts[0], contours[0] # Return counts and (optionally) the contour we used return (counts, contours) if return_contour else counts return out # Internal function to add arrows to a curve def _add_arrows_to_line2D( axes, line, arrow_locs=[0.2, 0.4, 0.6, 0.8], arrowstyle='-|>', arrowsize=1, dir=1, transform=None): """ Add arrows to a matplotlib.lines.Line2D at selected locations. Parameters: ----------- axes: Axes object as returned by axes command (or gca) line: Line2D object as returned by plot command arrow_locs: list of locations where to insert arrows, % of total length arrowstyle: style of the arrow arrowsize: size of the arrow transform: a matplotlib transform instance, default to data coordinates Returns: -------- arrows: list of arrows Based on https://stackoverflow.com/questions/26911898/ """ if not isinstance(line, mpl.lines.Line2D): raise ValueError("expected a matplotlib.lines.Line2D object") x, y = line.get_xdata(), line.get_ydata() arrow_kw = { "arrowstyle": arrowstyle, } color = line.get_color() use_multicolor_lines = isinstance(color, np.ndarray) if use_multicolor_lines: raise NotImplementedError("multicolor lines not supported") else: arrow_kw['color'] = color linewidth = line.get_linewidth() if isinstance(linewidth, np.ndarray): raise NotImplementedError("multiwidth lines not supported") else: arrow_kw['linewidth'] = linewidth if transform is None: transform = axes.transData # Compute the arc length along the curve s = np.cumsum(np.sqrt(np.diff(x) ** 2 + np.diff(y) ** 2)) arrows = [] for loc in arrow_locs: n = np.searchsorted(s, s[-1] * loc) # Figure out what direction to paint the arrow if dir == 1: arrow_tail = (x[n], y[n]) arrow_head = (np.mean(x[n:n + 2]), np.mean(y[n:n + 2])) elif dir == -1: # Orient the arrow in the other direction on the segment arrow_tail = (x[n + 1], y[n + 1]) arrow_head = (np.mean(x[n:n + 2]), np.mean(y[n:n + 2])) else: raise ValueError("unknown value for keyword 'dir'") p = mpl.patches.FancyArrowPatch( arrow_tail, arrow_head, transform=transform, lw=0, **arrow_kw) axes.add_patch(p) arrows.append(p) return arrows # # Function to compute Nyquist curve offsets # # This function computes a smoothly varying offset that starts and ends at # zero at the ends of a scaled segment. # def _compute_curve_offset(resp, mask, max_offset): # Compute the arc length along the curve s_curve = np.cumsum( np.sqrt(np.diff(resp.real) ** 2 + np.diff(resp.imag) ** 2)) # Initialize the offset offset = np.zeros(resp.size) arclen = np.zeros(resp.size) # Walk through the response and keep track of each continous component i, nsegs = 0, 0 while i < resp.size: # Skip the regular segment while i < resp.size and mask[i]: i += 1 # Increment the counter if i == resp.size: break # Keep track of the arclength arclen[i] = arclen[i-1] + np.abs(resp[i] - resp[i-1]) nsegs += 0.5 if i == resp.size: break # Save the starting offset of this segment seg_start = i # Walk through the scaled segment while i < resp.size and not mask[i]: i += 1 if i == resp.size: # See if we are done with this segment break # Keep track of the arclength arclen[i] = arclen[i-1] + np.abs(resp[i] - resp[i-1]) nsegs += 0.5 if i == resp.size: break # Save the ending offset of this segment seg_end = i # Now compute the scaling for this segment s_segment = arclen[seg_end-1] - arclen[seg_start] offset[seg_start:seg_end] = max_offset * s_segment/s_curve[-1] * \ np.sin(np.pi * (arclen[seg_start:seg_end] - arclen[seg_start])/s_segment) return offset # # Gang of Four plot # def gangof4_response(P, C, omega=None, Hz=False): """Compute the response of the "Gang of 4" transfer functions for a system. Generates a 2x2 frequency response for the "Gang of 4" sensitivity functions [T, PS; CS, S]. Parameters ---------- P, C : LTI Linear input/output systems (process and control) omega : array Range of frequencies (list or bounds) in rad/sec Returns ------- response : :class:`~control.FrequencyResponseData` Frequency response with inputs 'r' and 'd' and outputs 'y', and 'u' representing the 2x2 matrix of transfer functions in the Gang of 4. Examples -------- >>> P = ct.tf([1], [1, 1]) >>> C = ct.tf([2], [1]) >>> response = ct.gangof4_response(P, C) >>> lines = response.plot() """ if not P.issiso() or not C.issiso(): # TODO: Add MIMO go4 plots. raise ControlMIMONotImplemented( "Gang of four is currently only implemented for SISO systems.") # Compute the senstivity functions L = P * C S = feedback(1, L) T = L * S # Select a default range if none is provided # TODO: This needs to be made more intelligent if omega is None: omega = _default_frequency_range((P, C, S), Hz=Hz) # # bode_plot based implementation # # Compute the response of the Gang of 4 resp_T = T(1j * omega) resp_PS = (P * S)(1j * omega) resp_CS = (C * S)(1j * omega) resp_S = S(1j * omega) # Create a single frequency response data object with the underlying data data = np.empty((2, 2, omega.size), dtype=complex) data[0, 0, :] = resp_T data[0, 1, :] = resp_PS data[1, 0, :] = resp_CS data[1, 1, :] = resp_S return FrequencyResponseData( data, omega, outputs=['y', 'u'], inputs=['r', 'd'], title=f"Gang of Four for P={P.name}, C={C.name}", plot_phase=False) def gangof4_plot(P, C, omega=None, **kwargs): """Legacy Gang of 4 plot; use gangof4_response().plot() instead.""" return gangof4_response(P, C).plot(**kwargs) # # Singular values plot # def singular_values_response( sysdata, omega=None, omega_limits=None, omega_num=None, Hz=False): """Singular value response for a system. Computes the singular values for a system or list of systems over a (optional) frequency range. Parameters ---------- sysdata : LTI or list of LTI List of linear input/output systems (single system is OK). omega : array_like List of frequencies in rad/sec to be used for frequency response. omega_limits : array_like of two values Limits of the frequency vector to generate, in rad/s. omega_num : int Number of samples to plot. Default value (1000) set by config.defaults['freqplot.number_of_samples']. 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. Returns ------- response : FrequencyResponseData Frequency response with the number of outputs equal to the number of singular values in the response, and a single input. Examples -------- >>> omegas = np.logspace(-4, 1, 1000) >>> den = [75, 1] >>> G = ct.tf([[[87.8], [-86.4]], [[108.2], [-109.6]]], ... [[den, den], [den, den]]) >>> response = ct.singular_values_response(G, omega=omegas) """ # Convert the first argument to a list syslist = sysdata if isinstance(sysdata, (list, tuple)) else [sysdata] if any([not isinstance(sys, LTI) for sys in syslist]): ValueError("singular values can only be computed for LTI systems") # Compute the frequency responses for the systems responses = frequency_response( syslist, omega=omega, omega_limits=omega_limits, omega_num=omega_num, Hz=Hz, squeeze=False) # Calculate the singular values for each system in the list svd_responses = [] for response in responses: # Compute the singular values (permute indices to make things work) fresp_permuted = response.fresp.transpose((2, 0, 1)) sigma = np.linalg.svd(fresp_permuted, compute_uv=False).transpose() sigma_fresp = sigma.reshape(sigma.shape[0], 1, sigma.shape[1]) # Save the singular values as an FRD object svd_responses.append( FrequencyResponseData( sigma_fresp, response.omega, _return_singvals=True, outputs=[f'$\\sigma_{{{k+1}}}$' for k in range(sigma.shape[0])], inputs='inputs', dt=response.dt, plot_phase=False, sysname=response.sysname, plot_type='svplot', title=f"Singular values for {response.sysname}")) if isinstance(sysdata, (list, tuple)): return FrequencyResponseList(svd_responses) else: return svd_responses[0] def singular_values_plot( data, omega=None, *fmt, plot=None, omega_limits=None, omega_num=None, title=None, legend_loc='center right', **kwargs): """Plot the singular values for a system. Plot the singular values as a function of frequency for a system or list of systems. If multiple systems are plotted, each system in the list is plotted in a different color. Parameters ---------- data : list of `FrequencyResponseData` List of :class:`FrequencyResponseData` objects. For backward compatibility, a list of LTI systems can also be given. omega : array_like List of frequencies in rad/sec over to plot over. *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']. legend_loc : str, optional For plots with multiple lines, a legend will be included in the given location. Default is 'center right'. Use False to supress. **kwargs : :func:`matplotlib.pyplot.plot` keyword properties, optional Additional keywords passed to `matplotlib` to specify line properties. Returns ------- lines : array of Line2D 1-D array of Line2D objects. The size of the array matches the number of systems and the value of the array is a list of Line2D objects for that system. mag : ndarray (or list of ndarray if len(data) > 1)) If plot=False, magnitude of the response (deprecated). phase : ndarray (or list of ndarray if len(data) > 1)) If plot=False, phase in radians of the response (deprecated). omega : ndarray (or list of ndarray if len(data) > 1)) If plot=False, frequency in rad/sec (deprecated). Other Parameters ---------------- grid : bool If True, plot grid lines on gain and phase plots. Default is set by `config.defaults['freqplot.grid']`. 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. omega_num : int Number of samples to use for the frequeny range. Defaults to config.defaults['freqplot.number_of_samples']. Ignore if data is not a list of systems. plot : bool, optional (legacy) If given, `singular_values_plot` returns the legacy return values of magnitude, phase, and frequency. If False, just return the values with no plot. rcParams : dict Override the default parameters used for generating plots. Default is set up config.default['freqplot.rcParams']. """ # Keyword processing dB = config._get_param( 'freqplot', 'dB', 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) freqplot_rcParams = config._get_param( 'freqplot', 'rcParams', kwargs, _freqplot_defaults, pop=True) # If argument was a singleton, turn it into a tuple data = data if isinstance(data, (list, tuple)) else (data,) # Convert systems into frequency responses if any([isinstance(response, (StateSpace, TransferFunction)) for response in data]): responses = singular_values_response( data, omega=omega, omega_limits=omega_limits, omega_num=omega_num) 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=}") responses = data # Process (legacy) plot keyword if plot is not None: warnings.warn( "`singular_values_plot` return values of sigma, omega is " "deprecated; use singular_values_response()", DeprecationWarning) # Warn the user if we got past something that is not real-valued if any([not np.allclose(np.imag(response.fresp[:, 0, :]), 0) for response in responses]): warnings.warn("data has non-zero imaginary component") # Extract the data we need for plotting sigmas = [np.real(response.fresp[:, 0, :]) for response in responses] omegas = [response.omega for response in responses] # Legacy processing for no plotting case if plot is False: if len(data) == 1: return sigmas[0], omegas[0] else: return sigmas, omegas fig = plt.gcf() # get current figure (or create new one) ax_sigma = None # axes for plotting singular values # Get the current axes if they already exist for ax in fig.axes: if ax.get_label() == 'control-sigma': ax_sigma = ax # If no axes present, create them from scratch if ax_sigma is None: if len(fig.axes) > 0: # Create a new figure to avoid overwriting in the old one fig = plt.figure() with plt.rc_context(_freqplot_rcParams): ax_sigma = plt.subplot(111, label='control-sigma') # Handle color cycle manually as all singular values # of the same systems are expected to be of the same color color_cycle = plt.rcParams['axes.prop_cycle'].by_key()['color'] color_offset = 0 if len(ax_sigma.lines) > 0: last_color = ax_sigma.lines[-1].get_color() if last_color in color_cycle: color_offset = color_cycle.index(last_color) + 1 # Create a list of lines for the output out = np.empty(len(data), dtype=object) # Plot the singular values for each response for idx_sys, response in enumerate(responses): sigma = sigmas[idx_sys].transpose() # frequency first for plotting omega = omegas[idx_sys] / (2 * math.pi) if Hz else omegas[idx_sys] if response.isdtime(strict=True): nyq_freq = (0.5/response.dt) if Hz else (math.pi/response.dt) else: nyq_freq = None # See if the color was specified, otherwise rotate if kwargs.get('color', None) or any( [isinstance(arg, str) and any([c in arg for c in "bgrcmykw#"]) for arg in fmt]): color_arg = {} # color set by *fmt, **kwargs else: color_arg = {'color': color_cycle[ (idx_sys + color_offset) % len(color_cycle)]} # Decide on the system name sysname = response.sysname if response.sysname is not None \ else f"Unknown-{idx_sys}" # Plot the data if dB: with plt.rc_context(freqplot_rcParams): out[idx_sys] = ax_sigma.semilogx( omega, 20 * np.log10(sigma), *fmt, label=sysname, **color_arg, **kwargs) else: with plt.rc_context(freqplot_rcParams): out[idx_sys] = ax_sigma.loglog( omega, sigma, label=sysname, *fmt, **color_arg, **kwargs) # Plot the Nyquist frequency if nyq_freq is not None: ax_sigma.axvline( nyq_freq, linestyle='--', label='_nyq_freq_' + sysname, **color_arg) # If specific omega_limits were given, use them if omega_limits is not None: ax_sigma.set_xlim(omega_limits) # Add a grid to the plot + labeling if grid: ax_sigma.grid(grid, which='both') with plt.rc_context(freqplot_rcParams): ax_sigma.set_ylabel( "Singular Values [dB]" if dB else "Singular Values") ax_sigma.set_xlabel("Frequency [Hz]" if Hz else "Frequency [rad/sec]") # List of systems that are included in this plot lines, labels = _get_line_labels(ax_sigma) # Add legend if there is more than one system plotted if len(labels) > 1 and legend_loc is not False: with plt.rc_context(freqplot_rcParams): ax_sigma.legend(lines, labels, loc=legend_loc) # Add the title if title is None: title = "Singular values for " + ", ".join(labels) with plt.rc_context(freqplot_rcParams): fig.suptitle(title) # Legacy return processing if plot is not None: if len(responses) == 1: return sigmas[0], omegas[0] else: return sigmas, omegas return out # # Utility functions # # This section of the code contains some utility functions for # generating frequency domain plots # # Determine the frequency range to be used def _determine_omega_vector(syslist, omega_in, omega_limits, omega_num, Hz=None, feature_periphery_decades=None): """Determine the frequency range for a frequency-domain plot according to a standard logic. If omega_in and omega_limits are both None, then omega_out is computed on omega_num points according to a default logic defined by _default_frequency_range and tailored for the list of systems syslist, and omega_range_given is set to False. If omega_in is None but omega_limits is an array-like of 2 elements, then omega_out is computed with the function np.logspace on omega_num points within the interval [min, max] = [omega_limits[0], omega_limits[1]], and omega_range_given is set to True. If omega_in is not None, then omega_out is set to omega_in, and omega_range_given is set to True Parameters ---------- syslist : list of LTI List of linear input/output systems (single system is OK) omega_in : 1D array_like or None Frequency range specified by the user omega_limits : 1D array_like or None Frequency limits specified by the user omega_num : int Number of points to be used for the frequency range (if the frequency range is not user-specified) Hz : bool, optional If True, the limits (first and last value) of the frequencies are set to full decades in Hz so it fits plotting with logarithmic scale in Hz otherwise in rad/s. Omega is always returned in rad/sec. Returns ------- omega_out : 1D array Frequency range to be used omega_range_given : bool True if the frequency range was specified by the user, either through omega_in or through omega_limits. False if both omega_in and omega_limits are None. """ omega_range_given = True if omega_in is None: if omega_limits is None: omega_range_given = False # Select a default range if none is provided omega_out = _default_frequency_range( syslist, number_of_samples=omega_num, Hz=Hz, feature_periphery_decades=feature_periphery_decades) else: omega_limits = np.asarray(omega_limits) if len(omega_limits) != 2: raise ValueError("len(omega_limits) must be 2") omega_out = np.logspace(np.log10(omega_limits[0]), np.log10(omega_limits[1]), num=omega_num, endpoint=True) else: omega_out = np.copy(omega_in) return omega_out, omega_range_given # Compute reasonable defaults for axes def _default_frequency_range(syslist, Hz=None, number_of_samples=None, feature_periphery_decades=None): """Compute a default frequency range for frequency domain plots. This code looks at the poles and zeros of all of the systems that we are plotting and sets the frequency range to be one decade above and below the min and max feature frequencies, rounded to the nearest integer. If no features are found, it returns logspace(-1, 1) Parameters ---------- syslist : list of LTI List of linear input/output systems (single system is OK) Hz : bool, optional If True, the limits (first and last value) of the frequencies are set to full decades in Hz so it fits plotting with logarithmic scale in Hz otherwise in rad/s. Omega is always returned in rad/sec. number_of_samples : int, optional Number of samples to generate. The default value is read from ``config.defaults['freqplot.number_of_samples']. If None, then the default from `numpy.logspace` is used. feature_periphery_decades : float, optional Defines how many decades shall be included in the frequency range on both sides of features (poles, zeros). The default value is read from ``config.defaults['freqplot.feature_periphery_decades']``. Returns ------- omega : array Range of frequencies in rad/sec Examples -------- >>> G = ct.ss([[-1, -2], [3, -4]], [[5], [7]], [[6, 8]], [[9]]) >>> omega = ct._default_frequency_range(G) >>> omega.min(), omega.max() (0.1, 100.0) """ # Set default values for options number_of_samples = config._get_param( 'freqplot', 'number_of_samples', number_of_samples) feature_periphery_decades = config._get_param( 'freqplot', 'feature_periphery_decades', feature_periphery_decades, 1) # Find the list of all poles and zeros in the systems features = np.array(()) freq_interesting = [] # detect if single sys passed by checking if it is sequence-like if not hasattr(syslist, '__iter__'): syslist = (syslist,) for sys in syslist: try: # Add new features to the list if sys.isctime(): features_ = np.concatenate( (np.abs(sys.poles()), np.abs(sys.zeros()))) # Get rid of poles and zeros at the origin toreplace = np.isclose(features_, 0.0) if np.any(toreplace): features_ = features_[~toreplace] elif sys.isdtime(strict=True): fn = math.pi / sys.dt # TODO: What distance to the Nyquist frequency is appropriate? freq_interesting.append(fn * 0.9) features_ = np.concatenate((sys.poles(), sys.zeros())) # Get rid of poles and zeros on the real axis (imag==0) # * origin and real < 0 # * at 1.: would result in omega=0. (logaritmic plot!) toreplace = np.isclose(features_.imag, 0.0) & ( (features_.real <= 0.) | (np.abs(features_.real - 1.0) < 1.e-10)) if np.any(toreplace): features_ = features_[~toreplace] # TODO: improve (mapping pack to continuous time) features_ = np.abs(np.log(features_) / (1.j * sys.dt)) else: # TODO raise NotImplementedError( "type of system in not implemented now") features = np.concatenate((features, features_)) except NotImplementedError: pass # Make sure there is at least one point in the range if features.shape[0] == 0: features = np.array([1.]) if Hz: features /= 2. * math.pi features = np.log10(features) lsp_min = np.rint(np.min(features) - feature_periphery_decades) lsp_max = np.rint(np.max(features) + feature_periphery_decades) if Hz: lsp_min += np.log10(2. * math.pi) lsp_max += np.log10(2. * math.pi) if freq_interesting: lsp_min = min(lsp_min, np.log10(min(freq_interesting))) lsp_max = max(lsp_max, np.log10(max(freq_interesting))) # TODO: Add a check in discrete case to make sure we don't get aliasing # (Attention: there is a list of system but only one omega vector) # Set the range to be an order of magnitude beyond any features if number_of_samples: omega = np.logspace( lsp_min, lsp_max, num=number_of_samples, endpoint=True) else: omega = np.logspace(lsp_min, lsp_max, endpoint=True) return omega # Get labels for all lines in an axes def _get_line_labels(ax, use_color=True): labels, lines = [], [] last_color, counter = None, 0 # label unknown systems for i, line in enumerate(ax.get_lines()): label = line.get_label() if use_color and label.startswith("Unknown"): label = f"Unknown-{counter}" if last_color is None: last_color = line.get_color() elif last_color != line.get_color(): counter += 1 last_color = line.get_color() elif label[0] == '_': continue if label not in labels: lines.append(line) labels.append(label) return lines, labels # # Utility functions to create nice looking labels (KLD 5/23/11) # def get_pow1000(num): """Determine exponent for which significand of a number is within the range [1, 1000). """ # Based on algorithm from http://www.mail-archive.com/ # matplotlib-users@lists.sourceforge.net/msg14433.html, accessed 2010/11/7 # by Jason Heeris 2009/11/18 from decimal import Decimal from math import floor dnum = Decimal(str(num)) if dnum == 0: return 0 elif dnum < 0: dnum = -dnum return int(floor(dnum.log10() / 3)) def gen_prefix(pow1000): """Return the SI prefix for a power of 1000. """ # Prefixes according to Table 5 of [BIPM 2006] (excluding hecto, # deca, deci, and centi). if pow1000 < -8 or pow1000 > 8: raise ValueError( "Value is out of the range covered by the SI prefixes.") return ['Y', # yotta (10^24) 'Z', # zetta (10^21) 'E', # exa (10^18) 'P', # peta (10^15) 'T', # tera (10^12) 'G', # giga (10^9) 'M', # mega (10^6) 'k', # kilo (10^3) '', # (10^0) 'm', # milli (10^-3) r'$\mu$', # micro (10^-6) 'n', # nano (10^-9) 'p', # pico (10^-12) 'f', # femto (10^-15) 'a', # atto (10^-18) 'z', # zepto (10^-21) 'y'][8 - pow1000] # yocto (10^-24) # Function aliases bode = bode_plot nyquist = nyquist_plot gangof4 = gangof4_plot