diff --git a/.github/conda-env/doctest-env.yml b/.github/conda-env/doctest-env.yml index a03b64ae7..4c0d36728 100644 --- a/.github/conda-env/doctest-env.yml +++ b/.github/conda-env/doctest-env.yml @@ -7,7 +7,7 @@ dependencies: - numpy - matplotlib - scipy - - sphinx + - sphinx<8.2 - sphinx_rtd_theme - ipykernel - nbsphinx diff --git a/.github/scripts/set-conda-test-matrix.py b/.github/scripts/set-conda-test-matrix.py index 6bcd0fa6f..49714e7eb 100644 --- a/.github/scripts/set-conda-test-matrix.py +++ b/.github/scripts/set-conda-test-matrix.py @@ -27,5 +27,19 @@ 'blas_lib': cbl} conda_jobs.append(cjob) +# Make sure Windows jobs are included even if we didn't build any +windows_pythons = ['3.11'] # Whatever you want to test + +for py in windows_pythons: + for blas in combinations['windows']: + cjob = { + 'packagekey': f'windows-{py}', + 'os': 'windows', + 'python': py, + 'blas_lib': blas, + 'package_source': 'conda-forge' + } + conda_jobs.append(cjob) + matrix = { 'include': conda_jobs } print(json.dumps(matrix)) diff --git a/.github/workflows/install_examples.yml b/.github/workflows/install_examples.yml index cfbf40fe7..6893a99fb 100644 --- a/.github/workflows/install_examples.yml +++ b/.github/workflows/install_examples.yml @@ -20,7 +20,8 @@ jobs: --quiet --yes \ python=3.12 pip \ numpy matplotlib scipy \ - slycot pmw jupyter + slycot pmw jupyter \ + ipython!=9.0 - name: Install from source run: | diff --git a/.github/workflows/os-blas-test-matrix.yml b/.github/workflows/os-blas-test-matrix.yml index 263afb7a4..3661c0499 100644 --- a/.github/workflows/os-blas-test-matrix.yml +++ b/.github/workflows/os-blas-test-matrix.yml @@ -107,7 +107,6 @@ jobs: os: - 'ubuntu' - 'macos' - - 'windows' python: # build on one, expand matrix in conda-build from the Sylcot/conda-recipe/conda_build_config.yaml - '3.11' @@ -332,7 +331,13 @@ jobs: echo "libblas * *mkl" >> $CONDA_PREFIX/conda-meta/pinned ;; esac - conda install -c ./slycot-conda-pkgs slycot + if [ "${{ matrix.os }}" = "windows" ]; then + echo "Installing slycot from conda-forge on Windows" + conda install slycot + else + echo "Installing built conda package from local channel" + conda install -c ./slycot-conda-pkgs slycot + fi conda list - name: Test with pytest run: JOBNAME="$JOBNAME" pytest control/tests diff --git a/benchmarks/flatsys_bench.py b/benchmarks/flatsys_bench.py index 05a2e7066..a2f8ae1d2 100644 --- a/benchmarks/flatsys_bench.py +++ b/benchmarks/flatsys_bench.py @@ -7,7 +7,6 @@ import numpy as np import math -import control as ct import control.flatsys as flat import control.optimal as opt diff --git a/benchmarks/optestim_bench.py b/benchmarks/optestim_bench.py index 612ee6bb3..534d1024d 100644 --- a/benchmarks/optestim_bench.py +++ b/benchmarks/optestim_bench.py @@ -6,7 +6,6 @@ # used for optimization-based estimation. import numpy as np -import math import control as ct import control.optimal as opt diff --git a/benchmarks/optimal_bench.py b/benchmarks/optimal_bench.py index 997b5a241..bd0c0cd6b 100644 --- a/benchmarks/optimal_bench.py +++ b/benchmarks/optimal_bench.py @@ -6,7 +6,6 @@ # performance of the functions used for optimization-base control. import numpy as np -import math import control as ct import control.flatsys as fs import control.optimal as opt @@ -21,7 +20,6 @@ 'RK23': ('RK23', {}), 'RK23_sloppy': ('RK23', {'atol': 1e-4, 'rtol': 1e-2}), 'RK45': ('RK45', {}), - 'RK45': ('RK45', {}), 'RK45_sloppy': ('RK45', {'atol': 1e-4, 'rtol': 1e-2}), 'LSODA': ('LSODA', {}), } @@ -129,9 +127,6 @@ def time_optimal_lq_methods(integrator_name, minimizer_name, method): Tf = 10 timepts = np.linspace(0, Tf, 20) - # Create the basis function to use - basis = get_basis('poly', 12, Tf) - res = opt.solve_ocp( sys, timepts, x0, traj_cost, constraints, terminal_cost=term_cost, solve_ivp_method=integrator[0], solve_ivp_kwargs=integrator[1], @@ -223,8 +218,6 @@ def time_discrete_aircraft_mpc(minimizer_name): # compute the steady state values for a particular value of the input ud = np.array([0.8, -0.3]) xd = np.linalg.inv(np.eye(5) - A) @ B @ ud - yd = C @ xd - # provide constraints on the system signals constraints = [opt.input_range_constraint(sys, [-5, -6], [5, 6])] @@ -234,7 +227,6 @@ def time_discrete_aircraft_mpc(minimizer_name): cost = opt.quadratic_cost(model, Q, R, x0=xd, u0=ud) # Set the time horizon and time points - Tf = 3 timepts = np.arange(0, 6) * 0.2 # Get the minimizer parameters to use diff --git a/control/config.py b/control/config.py index 8da7e2fc2..12d7b0052 100644 --- a/control/config.py +++ b/control/config.py @@ -297,7 +297,7 @@ def use_legacy_defaults(version): Parameters ---------- version : string - Version number of the defaults desired. Ranges from '0.1' to '0.10.1'. + Version number of `python-control` to use for setting defaults. Examples -------- @@ -342,6 +342,14 @@ def use_legacy_defaults(version): # reset_defaults() # start from a clean slate + # Version 0.10.2: + if major == 0 and minor < 10 or (minor == 10 and patch < 2): + from math import inf + + # Reset Nyquist defaults + set_defaults('nyquist', arrows=2, max_curve_magnitude=20, + blend_fraction=0, indent_points=50) + # Version 0.9.2: if major == 0 and minor < 9 or (minor == 9 and patch < 2): from math import inf diff --git a/control/ctrlplot.py b/control/ctrlplot.py index dbdb4e1ec..50b2751d8 100644 --- a/control/ctrlplot.py +++ b/control/ctrlplot.py @@ -41,7 +41,10 @@ # # Customize axes (curvilinear grids, shared axes, etc) # # # Plot the data -# lines = np.full(ax_array.shape, []) +# lines = np.empty(ax_array.shape, dtype=object) +# for i in range(ax_array.shape[0]): +# for j in range(ax_array.shape[1]): +# lines[i, j] = [] # line_labels = _process_line_labels(label, ntraces, nrows, ncols) # color_offset, color_cycle = _get_color_offset(ax) # for i, j in itertools.product(range(nrows), range(ncols)): @@ -355,7 +358,7 @@ def _process_ax_keyword( the calling function to do the actual axis creation (needed for curvilinear grids that use the AxisArtist module). - Legacy behavior: some of the older plotting commands use a axes label + Legacy behavior: some of the older plotting commands use an axes label to identify the proper axes for plotting. This behavior is supported through the use of the label keyword, but will only work if shape == (1, 1) and squeeze == True. diff --git a/control/descfcn.py b/control/descfcn.py index bfe2d1a7e..6f3f5169d 100644 --- a/control/descfcn.py +++ b/control/descfcn.py @@ -9,7 +9,6 @@ import math from warnings import warn -import matplotlib.pyplot as plt import numpy as np import scipy @@ -521,16 +520,18 @@ def describing_function_plot( # Plot the Nyquist response cplt = dfresp.response.plot(**kwargs) - lines[0] = cplt.lines[0] # Return Nyquist lines for first system + ax = cplt.axes[0, 0] # Get the axes where the plot was made + lines[0] = np.concatenate( # Return Nyquist lines for first system + cplt.lines.flatten()).tolist() # Add the describing function curve to the plot - lines[1] = plt.plot(dfresp.N_vals.real, dfresp.N_vals.imag) + lines[1] = ax.plot(dfresp.N_vals.real, dfresp.N_vals.imag) # Label the intersection points if point_label: for pos, (a, omega) in zip(dfresp.positions, dfresp.intersections): # Add labels to the intersection points - plt.text(pos.real, pos.imag, point_label % (a, omega)) + ax.text(pos.real, pos.imag, point_label % (a, omega)) return ControlPlot(lines, cplt.axes, cplt.figure) diff --git a/control/freqplot.py b/control/freqplot.py index fe258d636..475467147 100644 --- a/control/freqplot.py +++ b/control/freqplot.py @@ -1100,13 +1100,14 @@ def gen_zero_centered_series(val_min, val_max, period): _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.arrows': 3, # 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.indent_points': 200, # number of points to insert + 'nyquist.max_curve_magnitude': 15, # rescale large values + 'nyquist.blend_fraction': 0.15, # when to start scaling '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 marker @@ -1638,6 +1639,10 @@ def nyquist_plot( The matplotlib axes to draw the figure on. If not specified and the current figure has a single axes, that axes is used. Otherwise, a new figure is created. + blend_fraction : float, optional + For portions of the Nyquist curve that are scaled to have a maximum + magnitude of `max_curve_magnitude`, begin a smooth rescaling at + this fraction of `max_curve_magnitude`. Default value is 0.15. 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 @@ -1654,7 +1659,7 @@ def nyquist_plot( portions of the contour are plotted using a different line style. 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 + label(s). If `data` is a list, strings should be specified for each system. label_freq : int, optional Label every nth frequency on the plot. If not specified, no labels @@ -1685,8 +1690,8 @@ def nyquist_plot( elements is equivalent to providing `omega_limits`. omega_num : int, optional Number of samples to use for the frequency range. Defaults to - `config.defaults['freqplot.number_of_samples']`. Ignored if data is - not a list of systems. + `config.defaults['freqplot.number_of_samples']`. Ignored if `data` + is not a system or list of systems. plot : bool, optional (legacy) If given, `nyquist_plot` returns the legacy return values of (counts, contours). If False, return the values with no plot. @@ -1751,8 +1756,8 @@ def nyquist_plot( to avoid poles, resulting 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. + curve use a dashdot line style and the scaled portions of the mirror + image use a dotted line style. Examples -------- @@ -1784,6 +1789,8 @@ def nyquist_plot( ax_user = ax max_curve_magnitude = config._get_param( 'nyquist', 'max_curve_magnitude', kwargs, _nyquist_defaults, pop=True) + blend_fraction = config._get_param( + 'nyquist', 'blend_fraction', kwargs, _nyquist_defaults, pop=True) max_curve_offset = config._get_param( 'nyquist', 'max_curve_offset', kwargs, _nyquist_defaults, pop=True) rcParams = config._get_param('ctrlplot', 'rcParams', kwargs, pop=True) @@ -1878,10 +1885,16 @@ def _parse_linestyle(style_name, allow_false=False): legend_loc, _, show_legend = _process_legend_keywords( kwargs, None, 'upper right') + # Figure out where the blended curve should start + if blend_fraction < 0 or blend_fraction > 1: + raise ValueError("blend_fraction must be between 0 and 1") + blend_curve_start = (1 - blend_fraction) * max_curve_magnitude + # 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 + out = np.empty((len(nyquist_responses), 4), dtype=object) + for i in range(len(nyquist_responses)): + for j in range(4): + out[i, j] = [] # unique list in each element for idx, response in enumerate(nyquist_responses): resp = response.response @@ -1892,20 +1905,31 @@ def _parse_linestyle(style_name, allow_false=False): # 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) + np.abs(resp) > blend_curve_start, + np.logical_not(np.isclose(splane_contour.real, 0))) 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]) + rescale_idx = (np.abs(resp) > blend_curve_start) + + if np.any(rescale_idx): # Only process if rescaling is needed + subset = resp[rescale_idx] + abs_subset = np.abs(subset) + unit_vectors = subset / abs_subset # Preserve phase/direction + + if blend_curve_start == max_curve_magnitude: + # Clip at max_curve_magnitude + resp[rescale_idx] = max_curve_magnitude * unit_vectors + else: + # Logistic scaling + newmag = blend_curve_start + \ + (max_curve_magnitude - blend_curve_start) * \ + (abs_subset - blend_curve_start) / \ + (abs_subset + max_curve_magnitude - 2 * blend_curve_start) + resp[rescale_idx] = newmag * unit_vectors # Get the label to use for the line label = response.sysname if line_labels is None else line_labels[idx] @@ -1913,10 +1937,10 @@ def _parse_linestyle(style_name, allow_false=False): # 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( + p = ax.plot( x_reg, y_reg, primary_style[0], color=color, label=label, **kwargs) c = p[0].get_color() - out[idx] += p + out[idx, 0] += 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 @@ -1928,54 +1952,56 @@ def _parse_linestyle(style_name, allow_false=False): 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( + out[idx, 1] += ax.plot( x_scl * (1 + curve_offset), y_scl * (1 + curve_offset), primary_style[1], color=c, **kwargs) else: - out[idx] += [None] + out[idx, 1] += [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) + p = ax.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( + out[idx, 2] += ax.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( + out[idx, 3] += ax.plot( x_scl * (1 - curve_offset), -y_scl * (1 - curve_offset), mirror_style[1], color=c, **kwargs) else: - out[idx] += [None] + out[idx, 3] += [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) + p = ax.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] + out[idx, 2] += [None] + out[idx, 3] += [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) + segment = 0 if 0 in rescale_idx else 1 # regular vs scaled + out[idx, segment] += ax.plot( + resp[0].real, resp[0].imag, start_marker, + color=c, markersize=start_marker_size) # Mark the -1 point - plt.plot([-1], [0], 'r+') + ax.plot([-1], [0], 'r+') # # Draw circles for gain crossover and sensitivity functions @@ -1987,16 +2013,16 @@ def _parse_linestyle(style_name, allow_false=False): # Display the unit circle, to read gain crossover frequency if unit_circle: - plt.plot(cos, sin, **config.defaults['nyquist.circle_style']) + ax.plot(cos, sin, **config.defaults['nyquist.circle_style']) # Draw circles for given magnitudes of sensitivity if ms_circles is not None: for ms in ms_circles: pos_x = -1 + (1/ms)*cos pos_y = (1/ms)*sin - plt.plot( + ax.plot( pos_x, pos_y, **config.defaults['nyquist.circle_style']) - plt.text(pos_x[label_pos], pos_y[label_pos], ms) + ax.text(pos_x[label_pos], pos_y[label_pos], ms) # Draw circles for given magnitudes of complementary sensitivity if mt_circles is not None: @@ -2006,17 +2032,17 @@ def _parse_linestyle(style_name, allow_false=False): rt = mt/(mt**2-1) # Mt radius pos_x = ct+rt*cos pos_y = rt*sin - plt.plot( + ax.plot( pos_x, pos_y, **config.defaults['nyquist.circle_style']) - plt.text(pos_x[label_pos], pos_y[label_pos], mt) + ax.text(pos_x[label_pos], pos_y[label_pos], mt) else: - _, _, ymin, ymax = plt.axis() + _, _, ymin, ymax = ax.axis() pos_y = np.linspace(ymin, ymax, 100) - plt.vlines( + ax.vlines( -0.5, ymin=ymin, ymax=ymax, **config.defaults['nyquist.circle_style']) - plt.text(-0.5, pos_y[label_pos], 1) + ax.text(-0.5, pos_y[label_pos], 1) # Label the frequencies of the points on the Nyquist curve if label_freq: @@ -2039,7 +2065,7 @@ def _parse_linestyle(style_name, allow_false=False): # np.round() is used because 0.99... appears # instead of 1.0, and this would otherwise be # truncated to 0. - plt.text(xpt, ypt, ' ' + + ax.text(xpt, ypt, ' ' + str(int(np.round(f / 1000 ** pow1000, 0))) + ' ' + prefix + 'Hz') @@ -2556,12 +2582,12 @@ def singular_values_plot( nyq_freq = None # Determine the color to use for this response - color = _get_color( + current_color = _get_color( color, fmt=fmt, offset=color_offset + idx_sys, color_cycle=color_cycle) # To avoid conflict with *fmt, only pass color kw if non-None - color_arg = {} if color is None else {'color': color} + color_arg = {} if current_color is None else {'color': current_color} # Decide on the system name sysname = response.sysname if response.sysname is not None \ diff --git a/control/iosys.py b/control/iosys.py index 110552138..42cf4094d 100644 --- a/control/iosys.py +++ b/control/iosys.py @@ -46,6 +46,24 @@ class NamedSignal(np.ndarray): This class modifies the `numpy.ndarray` class and allows signals to be accessed using the signal name in addition to indices and slices. + Signals can be either a 2D array, index by signal and time, or a 3D + array, indexed by signal, trace, and time. + + Attributes + ---------- + signal_labels : list of str + Label names for each of the signal elements in the signal. + trace_labels : list of str, optional + Label names for each of the traces in the signal (if multi-trace). + + Examples + -------- + >>> sys = ct.rss( + ... states=['p1', 'p2', 'p3'], inputs=['u1', 'u2'], outputs=['y']) + >>> resp = ct.step_response(sys) + >>> resp.states['p1', 'u1'] # Step response from u1 to p1 + NamedSignal(...) + """ def __new__(cls, input_array, signal_labels=None, trace_labels=None): # See https://numpy.org/doc/stable/user/basics.subclassing.html @@ -315,6 +333,9 @@ def _repr_html_(self): # Defaults to using __repr__; override in subclasses return None + def _repr_markdown_(self): + return self._repr_html_() + @property def repr_format(self): """String representation format. diff --git a/control/margins.py b/control/margins.py index d7c7992be..57e825c65 100644 --- a/control/margins.py +++ b/control/margins.py @@ -11,12 +11,17 @@ import numpy as np import scipy as sp -from . import frdata, freqplot, xferfcn +from . import frdata, freqplot, xferfcn, statesp from .exception import ControlMIMONotImplemented from .iosys import issiso +from .ctrlutil import mag2db +try: + from slycot import ab13md +except ImportError: + ab13md = None -__all__ = ['stability_margins', 'phase_crossover_frequencies', 'margin'] - +__all__ = ['stability_margins', 'phase_crossover_frequencies', 'margin', + 'disk_margins'] # private helper functions def _poly_iw(sys): @@ -168,6 +173,7 @@ def fun(wdt): return z, w + def _likely_numerical_inaccuracy(sys): # crude, conservative check for if # num(z)*num(1/z) << den(z)*den(1/z) for DT systems @@ -517,3 +523,137 @@ def margin(*args): % len(args)) return margin[0], margin[1], margin[3], margin[4] + + +def disk_margins(L, omega, skew=0.0, returnall=False): + """Disk-based stability margins of loop transfer function. + + Parameters + ---------- + L : `StateSpace` or `TransferFunction` + Linear SISO or MIMO loop transfer function. + omega : sequence of array_like + 1D array of (non-negative) frequencies (rad/s) at which + to evaluate the disk-based stability margins. + skew : float or array_like, optional + Skew parameter(s) for disk margin (default = 0.0): + + * skew = 0.0 "balanced" sensitivity function 0.5*(S - T) + * skew = 1.0 sensitivity function S + * skew = -1.0 complementary sensitivity function T + + returnall : bool, optional + If True, return frequency-dependent margins. + If False (default), return worst-case (minimum) margins. + + Returns + ------- + DM : float or array_like + Disk margin. + DGM : float or array_like + Disk-based gain margin. + DPM : float or array_like + Disk-based phase margin. + + Examples + -------- + >> omega = np.logspace(-1, 3, 1001) + >> P = control.ss([[0, 10], [-10, 0]], np.eye(2), [[1, 10], + [-10, 1]], 0) + >> K = control.ss([], [], [], [[1, -2], [0, 1]]) + >> L = P * K + >> DM, DGM, DPM = control.disk_margins(L, omega, skew=0.0) + """ + + # First argument must be a system + if not isinstance(L, (statesp.StateSpace, xferfcn.TransferFunction)): + raise ValueError( + "Loop gain must be state-space or transfer function object") + + # Loop transfer function must be square + if statesp.ss(L).B.shape[1] != statesp.ss(L).C.shape[0]: + raise ValueError("Loop gain must be square (n_inputs = n_outputs)") + + # Need slycot if L is MIMO, for mu calculation + if not L.issiso() and ab13md == None: + raise ControlMIMONotImplemented( + "Need slycot to compute MIMO disk_margins") + + # Get dimensions of feedback system + num_loops = statesp.ss(L).C.shape[0] + I = statesp.ss([], [], [], np.eye(num_loops)) + + # Loop sensitivity function + S = I.feedback(L) + + # Compute frequency response of the "balanced" (according + # to the skew parameter "sigma") sensitivity function [1-2] + ST = S + 0.5 * (skew - 1) * I + ST_mag, ST_phase, _ = ST.frequency_response(omega) + ST_jw = (ST_mag * np.exp(1j * ST_phase)) + if not L.issiso(): + ST_jw = ST_jw.transpose(2, 0, 1) + + # Frequency-dependent complex disk margin, computed using + # upper bound of the structured singular value, a.k.a. "mu", + # of (S + (skew - I)/2). + DM = np.zeros(omega.shape) + DGM = np.zeros(omega.shape) + DPM = np.zeros(omega.shape) + for ii in range(0, len(omega)): + # Disk margin (a.k.a. "alpha") vs. frequency + if L.issiso() and ab13md == None: + # For the SISO case, the norm on (S + (skew - I)/2) is + # unstructured, and can be computed as the magnitude + # of the frequency response. + DM[ii] = 1.0 / ST_mag[ii] + else: + # For the MIMO case, the norm on (S + (skew - I)/2) + # assumes a single complex uncertainty block diagonal + # uncertainty structure. AB13MD provides an upper bound + # on this norm at the given frequency omega[ii]. + DM[ii] = 1.0 / ab13md(ST_jw[ii], np.array(num_loops * [1]), + np.array(num_loops * [2]))[0] + + # Disk-based gain margin (dB) and phase margin (deg) + with np.errstate(divide='ignore', invalid='ignore'): + # Real-axis intercepts with the disk + gamma_min = (1 - 0.5 * DM[ii] * (1 - skew)) \ + / (1 + 0.5 * DM[ii] * (1 + skew)) + gamma_max = (1 + 0.5 * DM[ii] * (1 - skew)) \ + / (1 - 0.5 * DM[ii] * (1 + skew)) + + # Gain margin (dB) + DGM[ii] = mag2db(np.minimum(1 / gamma_min, gamma_max)) + if np.isnan(DGM[ii]): + DGM[ii] = float('inf') + + # Phase margin (deg) + if np.isinf(gamma_max): + DPM[ii] = 90.0 + else: + DPM[ii] = (1 + gamma_min * gamma_max) \ + / (gamma_min + gamma_max) + if abs(DPM[ii]) >= 1.0: + DPM[ii] = float('Inf') + else: + DPM[ii] = np.rad2deg(np.arccos(DPM[ii])) + + if returnall: + # Frequency-dependent disk margin, gain margin and phase margin + return DM, DGM, DPM + else: + # Worst-case disk margin, gain margin and phase margin + if DGM.shape[0] and not np.isinf(DGM).all(): + with np.errstate(all='ignore'): + gmidx = np.where(DGM == np.min(DGM)) + else: + gmidx = -1 + + if DPM.shape[0]: + pmidx = np.where(DPM == np.min(DPM)) + + return ( + float('inf') if DM.shape[0] == 0 else np.amin(DM), + float('inf') if gmidx == -1 else DGM[gmidx][0], + float('inf') if DPM.shape[0] == 0 else DPM[pmidx][0]) diff --git a/control/matlab/__init__.py b/control/matlab/__init__.py index 6414c9131..facaf28bb 100644 --- a/control/matlab/__init__.py +++ b/control/matlab/__init__.py @@ -6,7 +6,7 @@ This subpackage contains a number of functions that emulate some of the functionality of MATLAB. The intent of these functions is to -provide a simple interface to the python control systems library +provide a simple interface to the Python Control Systems Library (python-control) for people who are familiar with the MATLAB Control Systems Toolbox (tm). diff --git a/control/nichols.py b/control/nichols.py index 3c4edcdbd..98775ddaf 100644 --- a/control/nichols.py +++ b/control/nichols.py @@ -132,15 +132,15 @@ def nichols_plot( out[idx] = ax_nichols.plot(x, y, *fmt, label=label_, **kwargs) # Label the plot axes - plt.xlabel('Phase [deg]') - plt.ylabel('Magnitude [dB]') + ax_nichols.set_xlabel('Phase [deg]') + ax_nichols.set_ylabel('Magnitude [dB]') # Mark the -180 point - plt.plot([-180], [0], 'r+') + ax_nichols.plot([-180], [0], 'r+') # Add grid if grid: - nichols_grid() + nichols_grid(ax=ax_nichols) # List of systems that are included in this plot lines, labels = _get_line_labels(ax_nichols) diff --git a/control/nlsys.py b/control/nlsys.py index 30f06f819..04fda2e40 100644 --- a/control/nlsys.py +++ b/control/nlsys.py @@ -1572,12 +1572,13 @@ def input_output_response( If discontinuous inputs are given, the underlying SciPy numerical integration algorithms can sometimes produce erroneous results due to - the default tolerances that are used. The `ivp_method` and - `ivp_keywords` parameters can be used to tune the ODE solver and - produce better results. In particular, using 'LSODA' as the - `ivp_method` or setting the `rtol` parameter to a smaller value - (e.g. using ``ivp_kwargs={'rtol': 1e-4}``) can provide more accurate - results. + the default tolerances that are used. The `solve_ivp_method` and + `solve_ivp_keywords` parameters can be used to tune the ODE solver and + produce better results. In particular, using 'LSODA' as the + `solve_ivp_method`, setting the `rtol` parameter to a smaller value + (e.g. using ``solve_ivp_kwargs={'rtol': 1e-4}``), or setting the + maximum step size to a smaller value (e.g. ``solve_ivp_kwargs= + {'max_step': 0.01}``) can provide more accurate results. """ # diff --git a/control/optimal.py b/control/optimal.py index 3242ac3fb..6b60c5b25 100644 --- a/control/optimal.py +++ b/control/optimal.py @@ -746,9 +746,9 @@ def _compute_states_inputs(self, coeffs): states = self.last_states else: states = self._simulate_states(self.x, inputs) - self.last_x = self.x - self.last_states = states - self.last_coeffs = coeffs + self.last_x = self.x.copy() # save initial state + self.last_states = states # always a new object + self.last_coeffs = coeffs.copy() # save coefficients return states, inputs diff --git a/control/pzmap.py b/control/pzmap.py index 42ba8e087..6bf928f56 100644 --- a/control/pzmap.py +++ b/control/pzmap.py @@ -124,6 +124,17 @@ def plot(self, *args, **kwargs): """ return pole_zero_plot(self, *args, **kwargs) + def replot(self, cplt: ControlPlot): + """Update the pole/zero loci of an existing plot. + + Parameters + ---------- + cplt : `ControlPlot` + Graphics handles of the existing plot. + """ + pole_zero_replot(self, cplt) + + # Pole/zero map def pole_zero_map(sysdata): @@ -513,6 +524,35 @@ def _click_dispatcher(event): return ControlPlot(out, ax, fig, legend=legend) +def pole_zero_replot(pzmap_responses, cplt): + """Update the loci of a plot after zooming/panning. + + Parameters + ---------- + pzmap_responses : PoleZeroMap list + Responses to update. + cplt : ControlPlot + Collection of plot handles. + """ + + for idx, response in enumerate(pzmap_responses): + + # remove the old data + for l in cplt.lines[idx, 2]: + l.set_data([], []) + + # update the line data + if response.loci is not None: + + for il, locus in enumerate(response.loci.transpose()): + try: + cplt.lines[idx,2][il].set_data(real(locus), imag(locus)) + except IndexError: + # not expected, but more lines apparently needed + cplt.lines[idx,2].append(cplt.ax[0,0].plot( + real(locus), imag(locus))) + + # Utility function to find gain corresponding to a click event def _find_root_locus_gain(event, sys, ax): # Get the current axis limits to set various thresholds diff --git a/control/rlocus.py b/control/rlocus.py index c4ef8b40e..b7344e31d 100644 --- a/control/rlocus.py +++ b/control/rlocus.py @@ -33,7 +33,7 @@ # Root locus map -def root_locus_map(sysdata, gains=None): +def root_locus_map(sysdata, gains=None, xlim=None, ylim=None): """Compute the root locus map for an LTI system. Calculate the root locus by finding the roots of 1 + k * G(s) where G @@ -46,6 +46,10 @@ def root_locus_map(sysdata, gains=None): gains : array_like, optional Gains to use in computing plot of closed-loop poles. If not given, gains are chosen to include the main features of the root locus map. + xlim : tuple or list, optional + Set limits of x axis (see `matplotlib.axes.Axes.set_xlim`). + ylim : tuple or list, optional + Set limits of y axis (see `matplotlib.axes.Axes.set_ylim`). Returns ------- @@ -75,7 +79,7 @@ def root_locus_map(sysdata, gains=None): nump, denp = _systopoly1d(sys[0, 0]) if gains is None: - kvect, root_array, _, _ = _default_gains(nump, denp, None, None) + kvect, root_array, _, _ = _default_gains(nump, denp, xlim, ylim) else: kvect = np.atleast_1d(gains) root_array = _RLFindRoots(nump, denp, kvect) @@ -205,6 +209,11 @@ def root_locus_plot( # Plot the root loci cplt = responses.plot(grid=grid, **kwargs) + # Add a reaction to axis scale changes, if given LTI systems, and + # there is no set of pre-defined gains + if gains is None: + add_loci_recalculate(sysdata, cplt, cplt.axes[0,0]) + # Legacy processing: return locations of poles and zeros as a tuple if plot is True: return responses.loci, responses.gains @@ -212,6 +221,40 @@ def root_locus_plot( return ControlPlot(cplt.lines, cplt.axes, cplt.figure) +def add_loci_recalculate(sysdata, cplt, axis): + """Add a callback to re-calculate the loci data fitting a zoom action. + + Parameters + ---------- + sysdata: LTI object or list + Linear input/output systems (SISO only, for now). + cplt: ControlPlot + Collection of plot handles. + axis: matplotlib.axes.Axis + Axis on which callbacks are installed. + """ + + # if LTI, treat everything as a list of lti + if isinstance(sysdata, LTI): + sysdata = [sysdata] + + # check that we can actually recalculate the loci + if isinstance(sysdata, list) and all( + [isinstance(sys, LTI) for sys in sysdata]): + + # callback function for axis change (zoom, pan) events + # captures the sysdata object and cplt + def _zoom_adapter(_ax): + newresp = root_locus_map(sysdata, None, + _ax.get_xlim(), + _ax.get_ylim()) + newresp.replot(cplt) + + # connect the callback to axis changes + axis.callbacks.connect('xlim_changed', _zoom_adapter) + axis.callbacks.connect('ylim_changed', _zoom_adapter) + + def _default_gains(num, den, xlim, ylim): """Unsupervised gains calculation for root locus plot. @@ -288,7 +331,7 @@ def _default_gains(num, den, xlim, ylim): # Root locus is on imaginary axis (rare), use just y distance tolerance = y_tolerance elif y_tolerance == 0: - # Root locus is on imaginary axis (common), use just x distance + # Root locus is on real axis (common), use just x distance tolerance = x_tolerance else: tolerance = np.min([x_tolerance, y_tolerance]) diff --git a/control/sisotool.py b/control/sisotool.py index 78be86b16..2d4506f52 100644 --- a/control/sisotool.py +++ b/control/sisotool.py @@ -22,6 +22,7 @@ from .statesp import ss, summing_junction from .timeresp import step_response from .xferfcn import tf +from .rlocus import add_loci_recalculate _sisotool_defaults = { 'sisotool.initial_gain': 1 @@ -105,7 +106,7 @@ def sisotool(sys, initial_gain=None, xlim_rlocus=None, ylim_rlocus=None, fig = plt.gcf() if fig.canvas.manager.get_window_title() != 'Sisotool': plt.close(fig) - fig,axes = plt.subplots(2, 2) + fig, axes = plt.subplots(2, 2) fig.canvas.manager.set_window_title('Sisotool') else: axes = np.array(fig.get_axes()).reshape(2, 2) @@ -137,8 +138,8 @@ def sisotool(sys, initial_gain=None, xlim_rlocus=None, ylim_rlocus=None, # sys[0, 0], initial_gain=initial_gain, xlim=xlim_rlocus, # ylim=ylim_rlocus, plotstr=plotstr_rlocus, grid=rlocus_grid, # ax=fig.axes[1]) - ax_rlocus = fig.axes[1] - root_locus_map(sys[0, 0]).plot( + ax_rlocus = axes[0,1] # fig.axes[1] + cplt = root_locus_map(sys[0, 0]).plot( xlim=xlim_rlocus, ylim=ylim_rlocus, initial_gain=initial_gain, ax=ax_rlocus) if rlocus_grid is False: @@ -146,6 +147,9 @@ def sisotool(sys, initial_gain=None, xlim_rlocus=None, ylim_rlocus=None, from .grid import nogrid nogrid(sys.dt, ax=ax_rlocus) + # install a zoom callback on the root-locus axis + add_loci_recalculate(sys, cplt, ax_rlocus) + # Reset the button release callback so that we can update all plots fig.canvas.mpl_connect( 'button_release_event', partial( @@ -155,9 +159,8 @@ def sisotool(sys, initial_gain=None, xlim_rlocus=None, ylim_rlocus=None, def _click_dispatcher(event, sys, ax, bode_plot_params, tvect): # Zoom handled by specialized callback in rlocus, only handle gain plot - if event.inaxes == ax.axes and \ - plt.get_current_fig_manager().toolbar.mode not in \ - {'zoom rect', 'pan/zoom'}: + if event.inaxes == ax.axes: + fig = ax.figure # if a point is clicked on the rootlocus plot visually emphasize it diff --git a/control/tests/ctrlplot_test.py b/control/tests/ctrlplot_test.py index 958d855b2..bf8a075ae 100644 --- a/control/tests/ctrlplot_test.py +++ b/control/tests/ctrlplot_test.py @@ -243,6 +243,15 @@ def test_plot_ax_processing(resp_fcn, plot_fcn): # No response function available; just plot the data plot_fcn(*args, **kwargs, **plot_kwargs, ax=ax) + # Make sure the plot ended up in the right place + assert len(axs[0, 0].get_lines()) == 0 # upper left + assert len(axs[0, 1].get_lines()) != 0 # top middle + assert len(axs[1, 0].get_lines()) == 0 # lower left + if resp_fcn != ct.gangof4_response: + assert len(axs[1, 2].get_lines()) == 0 # lower right (normally empty) + else: + assert len(axs[1, 2].get_lines()) != 0 # gangof4 uses this axes + # Check to make sure original settings did not change assert fig._suptitle.get_text() == title assert fig._suptitle.get_fontsize() == titlesize diff --git a/control/tests/descfcn_test.py b/control/tests/descfcn_test.py index e91738e82..4fbb42913 100644 --- a/control/tests/descfcn_test.py +++ b/control/tests/descfcn_test.py @@ -190,11 +190,11 @@ def test_describing_function_plot(): cplt = response.plot() assert len(plt.gcf().get_axes()) == 1 # make sure there is a plot - assert len(cplt.lines[0]) == 4 and len(cplt.lines[1]) == 1 + assert len(cplt.lines[0]) == 5 and len(cplt.lines[1]) == 1 # Call plot directly cplt = ct.describing_function_plot(H_larger, F_saturation, amp, omega) - assert len(cplt.lines[0]) == 4 and len(cplt.lines[1]) == 1 + assert len(cplt.lines[0]) == 5 and len(cplt.lines[1]) == 1 def test_describing_function_exceptions(): diff --git a/control/tests/discrete_test.py b/control/tests/discrete_test.py index 9b87bd61b..7296c0f31 100644 --- a/control/tests/discrete_test.py +++ b/control/tests/discrete_test.py @@ -26,9 +26,11 @@ class Tsys: sys = rss(3, 1, 1) T.siso_ss1 = StateSpace(sys.A, sys.B, sys.C, sys.D, None) T.siso_ss1c = StateSpace(sys.A, sys.B, sys.C, sys.D, 0.0) - T.siso_ss1d = StateSpace(sys.A, sys.B, sys.C, sys.D, 0.1) - T.siso_ss2d = StateSpace(sys.A, sys.B, sys.C, sys.D, 0.2) - T.siso_ss3d = StateSpace(sys.A, sys.B, sys.C, sys.D, True) + + dsys = ct.sample_system(sys, 1) + T.siso_ss1d = StateSpace(dsys.A, dsys.B, dsys.C, dsys.D, 0.1) + T.siso_ss2d = StateSpace(dsys.A, dsys.B, dsys.C, dsys.D, 0.2) + T.siso_ss3d = StateSpace(dsys.A, dsys.B, dsys.C, dsys.D, True) # Two input, two output continuous-time system A = [[-3., 4., 2.], [-1., -3., 0.], [2., 5., 3.]] @@ -39,17 +41,18 @@ class Tsys: T.mimo_ss1c = StateSpace(A, B, C, D, 0) # Two input, two output discrete-time system - T.mimo_ss1d = StateSpace(A, B, C, D, 0.1) + T.mimo_ss1d = ct.sample_system(T.mimo_ss1c, 0.1) # Same system, but with a different sampling time - T.mimo_ss2d = StateSpace(A, B, C, D, 0.2) + T.mimo_ss2d = StateSpace( + T.mimo_ss1d.A, T.mimo_ss1d.B, T.mimo_ss1d.C, T.mimo_ss1d.D, 0.2) # Single input, single output continuus and discrete transfer function T.siso_tf1 = TransferFunction([1, 1], [1, 2, 1], None) - T.siso_tf1c = TransferFunction([1, 1], [1, 2, 1], 0) - T.siso_tf1d = TransferFunction([1, 1], [1, 2, 1], 0.1) - T.siso_tf2d = TransferFunction([1, 1], [1, 2, 1], 0.2) - T.siso_tf3d = TransferFunction([1, 1], [1, 2, 1], True) + T.siso_tf1c = TransferFunction([1, 1], [1, 0.2, 1], 0) + T.siso_tf1d = TransferFunction([1, 1], [1, 0.2, 0.1], 0.1) + T.siso_tf2d = TransferFunction([1, 1], [1, 0.2, 0.1], 0.2) + T.siso_tf3d = TransferFunction([1, 1], [1, 0.2, 0.1], True) return T diff --git a/control/tests/freqplot_test.py b/control/tests/freqplot_test.py index c0dfa4030..b3770486c 100644 --- a/control/tests/freqplot_test.py +++ b/control/tests/freqplot_test.py @@ -680,6 +680,39 @@ def test_display_margins(nsys, display_margins, gridkw, match): assert cplt.axes[0, 0].get_title() == '' +def test_singular_values_plot_colors(): + # Define some systems for testing + sys1 = ct.rss(4, 2, 2, strictly_proper=True) + sys2 = ct.rss(4, 2, 2, strictly_proper=True) + + # Get the default color cycle + color_cycle = plt.rcParams['axes.prop_cycle'].by_key()['color'] + + # Plot the systems individually and make sure line colors are OK + cplt = ct.singular_values_plot(sys1) + assert cplt.lines.size == 1 + assert len(cplt.lines[0]) == 2 + assert cplt.lines[0][0].get_color() == color_cycle[0] + assert cplt.lines[0][1].get_color() == color_cycle[0] + + cplt = ct.singular_values_plot(sys2) + assert cplt.lines.size == 1 + assert len(cplt.lines[0]) == 2 + assert cplt.lines[0][0].get_color() == color_cycle[1] + assert cplt.lines[0][1].get_color() == color_cycle[1] + plt.close('all') + + # Plot the systems as a list and make sure colors are OK + cplt = ct.singular_values_plot([sys1, sys2]) + assert cplt.lines.size == 2 + assert len(cplt.lines[0]) == 2 + assert len(cplt.lines[1]) == 2 + assert cplt.lines[0][0].get_color() == color_cycle[0] + assert cplt.lines[0][1].get_color() == color_cycle[0] + assert cplt.lines[1][0].get_color() == color_cycle[1] + assert cplt.lines[1][1].get_color() == color_cycle[1] + + if __name__ == "__main__": # # Interactive mode: generate plots for manual viewing diff --git a/control/tests/margin_test.py b/control/tests/margin_test.py index 43cd68ae3..679c1c685 100644 --- a/control/tests/margin_test.py +++ b/control/tests/margin_test.py @@ -14,7 +14,8 @@ from control import ControlMIMONotImplemented, FrequencyResponseData, \ StateSpace, TransferFunction, margin, phase_crossover_frequencies, \ - stability_margins + stability_margins, disk_margins, tf, ss +from control.exception import slycot_check s = TransferFunction.s @@ -372,3 +373,106 @@ def test_stability_margins_discrete(cnum, cden, dt, else: out = stability_margins(tf) assert_allclose(out, ref, rtol=rtol) + +def test_siso_disk_margin(): + # Frequencies of interest + omega = np.logspace(-1, 2, 1001) + + # Loop transfer function + L = tf(25, [1, 10, 10, 10]) + + # Balanced (S - T) disk-based stability margins + DM, DGM, DPM = disk_margins(L, omega, skew=0.0) + assert_allclose([DM], [0.46], atol=0.1) # disk margin of 0.46 + assert_allclose([DGM], [4.05], atol=0.1) # disk-based gain margin of 4.05 dB + assert_allclose([DPM], [25.8], atol=0.1) # disk-based phase margin of 25.8 deg + + # For SISO systems, the S-based (S) disk margin should match the third output + # of existing library "stability_margins", i.e., minimum distance from the + # Nyquist plot to -1. + _, _, SM = stability_margins(L)[:3] + DM = disk_margins(L, omega, skew=1.0)[0] + assert_allclose([DM], [SM], atol=0.01) + +def test_mimo_disk_margin(): + # Frequencies of interest + omega = np.logspace(-1, 3, 1001) + + # Loop transfer gain + P = ss([[0, 10], [-10, 0]], np.eye(2), [[1, 10], [-10, 1]], 0) # plant + K = ss([], [], [], [[1, -2], [0, 1]]) # controller + Lo = P * K # loop transfer function, broken at plant output + Li = K * P # loop transfer function, broken at plant input + + if slycot_check(): + # Balanced (S - T) disk-based stability margins at plant output + DMo, DGMo, DPMo = disk_margins(Lo, omega, skew=0.0) + assert_allclose([DMo], [0.3754], atol=0.1) # disk margin of 0.3754 + assert_allclose([DGMo], [3.3], atol=0.1) # disk-based gain margin of 3.3 dB + assert_allclose([DPMo], [21.26], atol=0.1) # disk-based phase margin of 21.26 deg + + # Balanced (S - T) disk-based stability margins at plant input + DMi, DGMi, DPMi = disk_margins(Li, omega, skew=0.0) + assert_allclose([DMi], [0.3754], atol=0.1) # disk margin of 0.3754 + assert_allclose([DGMi], [3.3], atol=0.1) # disk-based gain margin of 3.3 dB + assert_allclose([DPMi], [21.26], atol=0.1) # disk-based phase margin of 21.26 deg + else: + # Slycot not installed. Should throw exception. + with pytest.raises(ControlMIMONotImplemented,\ + match="Need slycot to compute MIMO disk_margins"): + DMo, DGMo, DPMo = disk_margins(Lo, omega, skew=0.0) + +def test_siso_disk_margin_return_all(): + # Frequencies of interest + omega = np.logspace(-1, 2, 1001) + + # Loop transfer function + L = tf(25, [1, 10, 10, 10]) + + # Balanced (S - T) disk-based stability margins + DM, DGM, DPM = disk_margins(L, omega, skew=0.0, returnall=True) + assert_allclose([omega[np.argmin(DM)]], [1.94],\ + atol=0.01) # sensitivity peak at 1.94 rad/s + assert_allclose([min(DM)], [0.46], atol=0.1) # disk margin of 0.46 + assert_allclose([DGM[np.argmin(DM)]], [4.05],\ + atol=0.1) # disk-based gain margin of 4.05 dB + assert_allclose([DPM[np.argmin(DM)]], [25.8],\ + atol=0.1) # disk-based phase margin of 25.8 deg + +def test_mimo_disk_margin_return_all(): + # Frequencies of interest + omega = np.logspace(-1, 3, 1001) + + # Loop transfer gain + P = ss([[0, 10], [-10, 0]], np.eye(2),\ + [[1, 10], [-10, 1]], 0) # plant + K = ss([], [], [], [[1, -2], [0, 1]]) # controller + Lo = P * K # loop transfer function, broken at plant output + Li = K * P # loop transfer function, broken at plant input + + if slycot_check(): + # Balanced (S - T) disk-based stability margins at plant output + DMo, DGMo, DPMo = disk_margins(Lo, omega, skew=0.0, returnall=True) + assert_allclose([omega[np.argmin(DMo)]], [omega[0]],\ + atol=0.01) # sensitivity peak at 0 rad/s (or smallest provided) + assert_allclose([min(DMo)], [0.3754], atol=0.1) # disk margin of 0.3754 + assert_allclose([DGMo[np.argmin(DMo)]], [3.3],\ + atol=0.1) # disk-based gain margin of 3.3 dB + assert_allclose([DPMo[np.argmin(DMo)]], [21.26],\ + atol=0.1) # disk-based phase margin of 21.26 deg + + # Balanced (S - T) disk-based stability margins at plant input + DMi, DGMi, DPMi = disk_margins(Li, omega, skew=0.0, returnall=True) + assert_allclose([omega[np.argmin(DMi)]], [omega[0]],\ + atol=0.01) # sensitivity peak at 0 rad/s (or smallest provided) + assert_allclose([min(DMi)], [0.3754],\ + atol=0.1) # disk margin of 0.3754 + assert_allclose([DGMi[np.argmin(DMi)]], [3.3],\ + atol=0.1) # disk-based gain margin of 3.3 dB + assert_allclose([DPMi[np.argmin(DMi)]], [21.26],\ + atol=0.1) # disk-based phase margin of 21.26 deg + else: + # Slycot not installed. Should throw exception. + with pytest.raises(ControlMIMONotImplemented,\ + match="Need slycot to compute MIMO disk_margins"): + DMo, DGMo, DPMo = disk_margins(Lo, omega, skew=0.0, returnall=True) diff --git a/control/tests/nyquist_test.py b/control/tests/nyquist_test.py index 42bb210c4..243a291d2 100644 --- a/control/tests/nyquist_test.py +++ b/control/tests/nyquist_test.py @@ -291,7 +291,7 @@ def test_nyquist_indent_default(indentsys): def test_nyquist_indent_dont(indentsys): # first value of default omega vector was 0.1, replaced by 0. for contour - # indent_radius is larger than 0.1 -> no extra quater circle around origin + # indent_radius is larger than 0.1 -> no extra quarter circle around origin with pytest.warns() as record: count, contour = ct.nyquist_response( indentsys, omega=[0, 0.2, 0.3, 0.4], indent_radius=.1007, @@ -406,15 +406,16 @@ def test_linestyle_checks(): # Set the line styles cplt = ct.nyquist_plot( sys, primary_style=[':', ':'], mirror_style=[':', ':']) - assert all([line.get_linestyle() == ':' for line in cplt.lines[0]]) + assert all([lines[0].get_linestyle() == ':' for lines in cplt.lines[0, :]]) # Set the line colors cplt = ct.nyquist_plot(sys, color='g') - assert all([line.get_color() == 'g' for line in cplt.lines[0]]) + assert all([line.get_color() == 'g' for line in cplt.lines[0, 0]]) # Turn off the mirror image cplt = ct.nyquist_plot(sys, mirror_style=False) - assert cplt.lines[0][2:] == [None, None] + assert cplt.lines[0, 2] == [None] + assert cplt.lines[0, 3] == [None] with pytest.raises(ValueError, match="invalid 'primary_style'"): ct.nyquist_plot(sys, primary_style=False) @@ -428,6 +429,7 @@ def test_linestyle_checks(): ct.nyquist_plot(sys, primary_style=':', mirror_style='-.') @pytest.mark.usefixtures("editsdefaults") +@pytest.mark.xfail(reason="updated code avoids warning") def test_nyquist_legacy(): ct.use_legacy_defaults('0.9.1') @@ -526,6 +528,42 @@ def test_no_indent_pole(): sys, warn_encirclements=False, indent_direction='none') +def test_nyquist_rescale(): + sys = 2 * ct.tf([1], [1, 1]) * ct.tf([1], [1, 0])**2 + sys.name = 'How example' + + # Default case + resp = ct.nyquist_response(sys, indent_direction='left') + cplt = resp.plot(label='default [0.15]') + assert len(cplt.lines[0, 0]) == 2 + assert all([len(cplt.lines[0, i]) == 1 for i in range(1, 4)]) + + # Sharper corner + cplt = ct.nyquist_plot( + sys*4, indent_direction='left', + max_curve_magnitude=17, blend_fraction=0.05, label='fraction=0.05') + assert len(cplt.lines[0, 0]) == 2 + assert all([len(cplt.lines[0, i]) == 1 for i in range(1, 4)]) + + # More gradual corner + cplt = ct.nyquist_plot( + sys*0.25, indent_direction='left', + max_curve_magnitude=13, blend_fraction=0.25, label='fraction=0.25') + assert len(cplt.lines[0, 0]) == 2 + assert all([len(cplt.lines[0, i]) == 1 for i in range(1, 4)]) + + # No corner + cplt = ct.nyquist_plot( + sys*12, indent_direction='left', + max_curve_magnitude=19, blend_fraction=0, label='fraction=0') + assert len(cplt.lines[0, 0]) == 2 + assert all([len(cplt.lines[0, i]) == 1 for i in range(1, 4)]) + + # Bad value + with pytest.raises(ValueError, match="blend_fraction must be between"): + ct.nyquist_plot(sys, indent_direction='left', blend_fraction=1.2) + + if __name__ == "__main__": # # Interactive mode: generate plots for manual viewing @@ -566,8 +604,8 @@ def test_no_indent_pole(): sys = 3 * (s+6)**2 / (s * (s**2 + 1e-4 * s + 1)) plt.figure() ct.nyquist_plot(sys) - ct.nyquist_plot(sys, max_curve_magnitude=15) - ct.nyquist_plot(sys, indent_radius=1e-6, max_curve_magnitude=25) + ct.nyquist_plot(sys, max_curve_magnitude=10) + ct.nyquist_plot(sys, indent_radius=1e-6, max_curve_magnitude=20) print("Unusual Nyquist plot") sys = ct.tf([1], [1, 3, 2]) * ct.tf([1], [1, 0, 1]) @@ -595,3 +633,29 @@ def test_no_indent_pole(): plt.figure() cplt = ct.nyquist_plot([sys, sys1, sys2]) cplt.set_plot_title("Mixed FRD, tf data") + + plt.figure() + print("Jon How example") + test_nyquist_rescale() + + # + # Save the figures in a PDF file for later comparisons + # + import subprocess + from matplotlib.backends.backend_pdf import PdfPages + from datetime import date + + # Create the file to store figures + try: + git_info = subprocess.check_output( + ['git', 'describe'], text=True).strip() + except subprocess.CalledProcessError: + git_info = 'UNKNOWN-REPO-INFO' + pdf = PdfPages( + f'nyquist_gallery-{git_info}-{date.today().isoformat()}.pdf') + + # Go through each figure and save it + for fignum in plt.get_fignums(): + pdf.savefig(plt.figure(fignum)) + + pdf.close() diff --git a/control/tests/phaseplot_test.py b/control/tests/phaseplot_test.py index fc4edcbea..ac5249948 100644 --- a/control/tests/phaseplot_test.py +++ b/control/tests/phaseplot_test.py @@ -167,7 +167,7 @@ def invpend_ode(t, x, m=0, l=0, b=0, g=0): ct.phase_plane_plot( invpend_ode, [-5, 5, 2, 2], params={'stuff': (1, 1, 0.2, 1)}, plot_streamlines=True) - + with pytest.raises(ValueError, match="gridtype must be 'meshgrid' when using streamplot"): ct.phase_plane_plot(ct.rss(2, 1, 1), plot_streamlines=False, plot_streamplot=True, gridtype='boxgrid') @@ -190,7 +190,7 @@ def invpend_ode(t, x, m=0, l=0, b=0, g=0): sys, [-12, 12, -10, 10], 15, gridspec=[2, 9], plot_streamlines=True, plot_separatrices=False, suppress_warnings=True) - + @pytest.mark.usefixtures('mplcleanup') def test_phase_plot_zorder(): # some of these tests are a bit akward since the streamlines and separatrices @@ -211,7 +211,7 @@ def get_zorders(cplt): assert cplt.lines[3] == None or isinstance(cplt.lines[3], mpl.streamplot.StreamplotSet) streamplot = max(cplt.lines[3].lines.get_zorder(), cplt.lines[3].arrows.get_zorder()) if cplt.lines[3] else None return streamlines, quiver, streamplot, separatrices, equilpoints - + def assert_orders(streamlines, quiver, streamplot, separatrices, equilpoints): print(streamlines, quiver, streamplot, separatrices, equilpoints) if streamlines is not None: @@ -261,8 +261,6 @@ def sys(t, x): # make sure changing the norm at least doesn't throw an error ct.phase_plane_plot(sys, plot_streamplot=dict(vary_color=True, norm=mpl.colors.LogNorm())) - - @pytest.mark.usefixtures('mplcleanup') def test_basic_phase_plots(savefigs=False): diff --git a/control/tests/statesp_test.py b/control/tests/statesp_test.py index 0806696f0..3c1411f04 100644 --- a/control/tests/statesp_test.py +++ b/control/tests/statesp_test.py @@ -136,6 +136,7 @@ def test_constructor(self, sys322ABCD, dt, argfun): ((np.ones((3, 3)), np.ones((3, 2)), np.ones((2, 3)), np.ones((2, 3))), ValueError, r"Incompatible dimensions of D matrix; expected \(2, 2\)"), + (([1j], 2, 3, 0), TypeError, "real number, not 'complex'"), ]) def test_constructor_invalid(self, args, exc, errmsg): """Test invalid input to StateSpace() constructor""" @@ -1440,6 +1441,7 @@ def test_html_repr(gmats, ref, dt, dtref, repr_type, num_format, editsdefaults): dt_html = dtref.format(dt=dt, fmt=defaults['statesp.latex_num_format']) ref_html = ref[refkey].format(dt=dt_html) assert g._repr_html_() == ref_html + assert g._repr_html_() == g._repr_markdown_() @pytest.mark.parametrize( diff --git a/control/tests/xferfcn_test.py b/control/tests/xferfcn_test.py index 87c852395..d3db08ef6 100644 --- a/control/tests/xferfcn_test.py +++ b/control/tests/xferfcn_test.py @@ -49,6 +49,10 @@ def test_constructor_bad_input_type(self): [[4, 5], [6, 7]]], [[[6, 7], [4, 5]], [[2, 3]]]) + + with pytest.raises(TypeError, match="unsupported data type"): + ct.tf([1j], [1, 2, 3]) + # good input TransferFunction([[[0, 1], [2, 3]], [[4, 5], [6, 7]]], diff --git a/control/timeresp.py b/control/timeresp.py index bd549589a..3c49d213e 100644 --- a/control/timeresp.py +++ b/control/timeresp.py @@ -551,10 +551,10 @@ def outputs(self): def states(self): """Time response state vector. - Time evolution of the state vector, indexed indexed by either the - state and time (if only a single trace is given) or the state, trace, - and time (for multiple traces). See `TimeResponseData.squeeze` - for a description of how this can be modified using the `squeeze` + Time evolution of the state vector, indexed by either the state and + time (if only a single trace is given) or the state, trace, and + time (for multiple traces). See `TimeResponseData.squeeze` for a + description of how this can be modified using the `squeeze` keyword. Input and output signal names can be used to index the data in @@ -616,9 +616,9 @@ def inputs(self): def _legacy_states(self): """Time response state vector (legacy version). - Time evolution of the state vector, indexed indexed by either the - state and time (if only a single trace is given) or the state, - trace, and time (for multiple traces). + Time evolution of the state vector, indexed by either the state and + time (if only a single trace is given) or the state, trace, and + time (for multiple traces). The `legacy_states` property is not affected by the `squeeze` keyword and hence it will always have these dimensions. diff --git a/control/xferfcn.py b/control/xferfcn.py index 16d7c5054..8e51534d7 100644 --- a/control/xferfcn.py +++ b/control/xferfcn.py @@ -1350,7 +1350,7 @@ def _c2d_matched(sysC, Ts, **kwargs): zpoles[idx] = z pregainden[idx] = 1 - z zgain = np.multiply.reduce(pregainnum) / np.multiply.reduce(pregainden) - gain = sysC.dcgain() / zgain + gain = sysC.dcgain() / zgain.real sysDnum, sysDden = zpk2tf(zzeros, zpoles, gain) return TransferFunction(sysDnum, sysDden, Ts, **kwargs) @@ -1954,6 +1954,7 @@ def _clean_part(data, name=""): """ valid_types = (int, float, complex, np.number) + unsupported_types = (complex, np.complexfloating) valid_collection = (list, tuple, ndarray) if isinstance(data, np.ndarray) and data.ndim == 2 and \ @@ -1998,8 +1999,11 @@ def _clean_part(data, name=""): for i in range(out.shape[0]): for j in range(out.shape[1]): for k in range(len(out[i, j])): - if isinstance(out[i, j][k], (int, np.int32, np.int64)): + if isinstance(out[i, j][k], (int, np.integer)): out[i, j][k] = float(out[i, j][k]) + elif isinstance(out[i, j][k], unsupported_types): + raise TypeError( + f"unsupported data type: {type(out[i, j][k])}") return out diff --git a/doc/Makefile b/doc/Makefile index 493fd7da5..4029dd70f 100644 --- a/doc/Makefile +++ b/doc/Makefile @@ -1,7 +1,7 @@ # Makefile for python-control Sphinx documentation # RMM, 15 Jan 2025 -FIGS = figures/classes.pdf +FIGS = figures/classes.svg RST_FIGS = figures/flatsys-steering-compare.png \ figures/iosys-predprey-open.png \ figures/timeplot-servomech-combined.png \ diff --git a/doc/_templates/extended-class-template.rst b/doc/_templates/extended-class-template.rst new file mode 100644 index 000000000..6e1e4ccd7 --- /dev/null +++ b/doc/_templates/extended-class-template.rst @@ -0,0 +1,9 @@ +{{ fullname | escape | underline}} + +.. currentmodule:: {{ module }} + +.. autoclass:: {{ objname }} + :no-members: + :show-inheritance: + :no-inherited-members: + :no-special-members: diff --git a/doc/classes.rst b/doc/classes.rst index 0ab508a3a..000761724 100644 --- a/doc/classes.rst +++ b/doc/classes.rst @@ -14,10 +14,10 @@ systems (both linear time-invariant and nonlinear). They are usually created from factory functions such as :func:`tf` and :func:`ss`, so the user should normally not need to instantiate these directly. -The following figure illustrates the relationship between the classes. +The following figure illustrates the relationship between the classes: -.. image:: figures/classes.pdf - :width: 800 +.. figure:: figures/classes.svg + :width: 640 :align: center .. autosummary:: @@ -34,6 +34,17 @@ The following figure illustrates the relationship between the classes. InterconnectedSystem LinearICSystem +The time response of an input/output system is represented using a +special :class:`NamedSignal` class that allows the individual signal +elements to be access using signal names in place of integer offsets: + +.. autosummary:: + :toctree: generated/ + :template: extended-class-template.rst + :nosignatures: + + NamedSignal + Response and Plotting Classes ============================= @@ -95,5 +106,5 @@ operations: optimal.OptimalEstimationProblem optimal.OptimalEstimationResult -More informaton on the functions used to create these classes can be +More information on the functions used to create these classes can be found in the :ref:`nonlinear-systems` chapter. diff --git a/doc/develop.rst b/doc/develop.rst index c9b6738a8..3ab4f8a94 100644 --- a/doc/develop.rst +++ b/doc/develop.rst @@ -110,8 +110,8 @@ Filenames * Source files are lower case, usually less than 10 characters (and 8 or less is better). -* Unit tests (in `control/tests/`) are of the form `module_test.py` or - `module_function.py`. +* Unit tests (in `control/tests/`) are of the form `module_test.py`, + `module_functionality_test.py`, or `functionality_test.py`. Class names diff --git a/doc/examples.rst b/doc/examples.rst index 2937fecab..b8d71807a 100644 --- a/doc/examples.rst +++ b/doc/examples.rst @@ -39,6 +39,7 @@ other sources. examples/mrac_siso_lyapunov examples/markov examples/era_msd + examples/disk_margins Jupyter Notebooks ================= diff --git a/doc/examples/disk_margins.py b/doc/examples/disk_margins.py new file mode 120000 index 000000000..a1dbcb7b1 --- /dev/null +++ b/doc/examples/disk_margins.py @@ -0,0 +1 @@ +../../examples/disk_margins.py \ No newline at end of file diff --git a/doc/examples/disk_margins.rst b/doc/examples/disk_margins.rst new file mode 100644 index 000000000..e7938f4ac --- /dev/null +++ b/doc/examples/disk_margins.rst @@ -0,0 +1,19 @@ +Disk margin example +------------------------------------------ + +This example demonstrates the use of the `disk_margins` routine +to compute robust stability margins for a feedback system, i.e., +variation in gain and phase one or more loops. The SISO examples +are drawn from the published paper and the MIMO example is the +"spinning satellite" example from the MathWorks documentation. + +Code +.... +.. literalinclude:: disk_margins.py + :language: python + :linenos: + +Notes +..... +1. The environment variable `PYCONTROL_TEST_EXAMPLES` is used for +testing to turn off plotting of the outputs. diff --git a/doc/figures/Makefile b/doc/figures/Makefile index 1ca54b372..26bdf22c2 100644 --- a/doc/figures/Makefile +++ b/doc/figures/Makefile @@ -2,7 +2,7 @@ # RMM, 26 Dec 2024 # List of figures that need to be created (first figure generated is OK) -FIGS = classes.pdf +FIGS = classes.svg # Location of the control package SRCDIR = ../.. @@ -12,5 +12,5 @@ all: $(FIGS) clean: /bin/rm -f $(FIGS) -classes.pdf: classes.fig - fig2dev -Lpdf $< $@ +classes.svg: classes.fig + fig2dev -Lsvg $< $@ diff --git a/doc/figures/classes.fig b/doc/figures/classes.fig index 4e63b8bff..17c112cc7 100644 --- a/doc/figures/classes.fig +++ b/doc/figures/classes.fig @@ -1,4 +1,5 @@ -#FIG 3.2 Produced by xfig version 3.2.8b +#FIG 3.2 Produced by xfig version 3.2.9 +#encoding: UTF-8 Landscape Center Inches @@ -37,12 +38,13 @@ Single 2 1 0 2 1 7 50 -1 -1 0.000 0 0 7 0 1 2 1 0 1.00 60.00 90.00 6525 1950 7050 2550 -4 1 1 50 -1 16 12 0.0000 4 210 2115 4050 3675 InterconnectedSystem\001 -4 1 1 50 -1 16 12 0.0000 4 165 1605 7950 3675 TransferFunction\001 -4 1 1 50 -1 0 12 0.0000 4 150 345 7050 2775 LTI\001 -4 1 1 50 -1 16 12 0.0000 4 210 1830 5175 2775 NonlinearIOSystem\001 -4 1 1 50 -1 16 12 0.0000 4 210 1095 6150 3675 StateSpace\001 -4 1 1 50 -1 16 12 0.0000 4 210 1500 5175 4575 LinearICSystem\001 -4 2 1 50 -1 16 12 0.0000 4 210 1035 3375 3225 FlatSystem\001 -4 0 1 50 -1 16 12 0.0000 4 165 420 8400 3225 FRD\001 -4 1 1 50 -1 16 12 0.0000 4 210 1770 6300 1875 InputOutputSystem\001 +4 1 1 50 -1 16 12 0.0000 4 191 1958 4050 3675 InterconnectedSystem\001 +4 1 1 50 -1 16 12 0.0000 4 151 1496 7950 3675 TransferFunction\001 +4 1 1 50 -1 0 12 0.0000 4 133 305 7050 2775 LTI\001 +4 1 1 50 -1 16 12 0.0000 4 191 1705 5175 2775 NonlinearIOSystem\001 +4 1 1 50 -1 16 12 0.0000 4 190 1016 6150 3675 StateSpace\001 +4 1 1 50 -1 16 12 0.0000 4 191 1394 5175 4575 LinearICSystem\001 +4 2 1 50 -1 16 12 0.0000 4 191 970 3375 3225 FlatSystem\001 +4 0 1 50 -1 16 12 0.0000 4 145 384 8400 3225 FRD\001 +4 1 1 50 -1 16 12 0.0000 4 191 1681 6300 1875 InputOutputSystem\001 +4 1 7 50 -1 16 12 0.0000 4 22 21 5175 4800 .\001 diff --git a/doc/figures/classes.pdf b/doc/figures/classes.pdf deleted file mode 100644 index 2c51b0193..000000000 Binary files a/doc/figures/classes.pdf and /dev/null differ diff --git a/doc/figures/classes.svg b/doc/figures/classes.svg new file mode 100644 index 000000000..98fedb596 --- /dev/null +++ b/doc/figures/classes.svg @@ -0,0 +1,151 @@ + + + + + + + +. + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +InterconnectedSystem + +TransferFunction + +LTI + +NonlinearIOSystem + +StateSpace + +LinearICSystem + +FlatSystem + +FRD + +InputOutputSystem + + + + + + + + + + + diff --git a/doc/figures/freqplot-nyquist-custom.png b/doc/figures/freqplot-nyquist-custom.png index 5cd2c19d0..29e7e476d 100644 Binary files a/doc/figures/freqplot-nyquist-custom.png and b/doc/figures/freqplot-nyquist-custom.png differ diff --git a/doc/figures/freqplot-nyquist-default.png b/doc/figures/freqplot-nyquist-default.png index c511509fa..968987d2c 100644 Binary files a/doc/figures/freqplot-nyquist-default.png and b/doc/figures/freqplot-nyquist-default.png differ diff --git a/doc/functions.rst b/doc/functions.rst index d657fd431..8432f7fcf 100644 --- a/doc/functions.rst +++ b/doc/functions.rst @@ -142,6 +142,7 @@ Frequency domain analysis: bandwidth dcgain + disk_margins linfnorm margin stability_margins diff --git a/doc/intro.rst b/doc/intro.rst index e1e5fb8e6..0054bb668 100644 --- a/doc/intro.rst +++ b/doc/intro.rst @@ -126,7 +126,7 @@ some things to keep in mind: * Vectors and matrices used as arguments to functions can be written using lists, with commas required between elements and column - vectors implemented as nested list . So [1 2 3] must be written as + vectors implemented as nested lists. So [1 2 3] must be written as [1, 2, 3] and matrices are written using 2D nested lists, e.g., [[1, 2], [3, 4]]. * Functions that in MATLAB would return variable numbers of values @@ -150,12 +150,12 @@ This documentation has a number of notional conventions and functionality: Manual, which contains documentation for all functions, classes, configurable default parameters, and other detailed information. -* Class, functions, and methods with additional documentation appear +* Classes, functions, and methods with additional documentation appear in a bold, code font that link to the Reference Manual. Example: `ss`. * Links to other sections appear in blue. Example: :ref:`nonlinear-systems`. -* Parameters appear in a (non-bode) code font, as do code fragments. +* Parameters appear in a (non-bold) code font, as do code fragments. Example: `omega`. * Example code is contained in code blocks that can be copied using diff --git a/doc/linear.rst b/doc/linear.rst index b7b8f7137..200260303 100644 --- a/doc/linear.rst +++ b/doc/linear.rst @@ -39,7 +39,7 @@ of linear time-invariant (LTI) systems: y &= C x + D u where :math:`u` is the input, :math:`y` is the output, and :math:`x` -is the state. +is the state. All vectors and matrices must be real-valued. To create a state space system, use the :func:`ss` function: @@ -94,7 +94,8 @@ transfer functions {b_0 s^n + b_1 s^{n-1} + \cdots + b_n}, where :math:`n` is greater than or equal to :math:`m` for a proper -transfer function. Improper transfer functions are also allowed. +transfer function. Improper transfer functions are also allowed. All +coefficients must be real-valued. To create a transfer function, use the :func:`tf` function:: @@ -530,6 +531,8 @@ equilibrium values (thereby keeping the input/output gain unchanged at zero frequency ["DC"]). +.. _displaying-lti-system-information: + Displaying LTI System Information ================================= diff --git a/doc/phaseplot.rst b/doc/phaseplot.rst index d2a3e6353..8c014415f 100644 --- a/doc/phaseplot.rst +++ b/doc/phaseplot.rst @@ -43,7 +43,7 @@ on a grid, equilibrium points and separatrices if they exist. A variety of options are available to modify the information that is plotted, including plotting a grid of vectors instead of streamlines, plotting streamlines from arbitrary starting points and turning on and off -various features of the plot. +various features of the plot. To illustrate some of these possibilities, consider a phase plane plot for an inverted pendulum system, which is created using a mesh grid: diff --git a/doc/releases.rst b/doc/releases.rst index 88a76775a..7bc5c0f46 100644 --- a/doc/releases.rst +++ b/doc/releases.rst @@ -27,6 +27,7 @@ the ability to index systems and signal using signal labels. .. toctree:: :maxdepth: 1 + releases/0.10.2-notes releases/0.10.1-notes releases/0.10.0-notes diff --git a/doc/releases/0.10.1-notes.rst b/doc/releases/0.10.1-notes.rst index dd0939021..8b99100f2 100644 --- a/doc/releases/0.10.1-notes.rst +++ b/doc/releases/0.10.1-notes.rst @@ -2,8 +2,8 @@ .. _version-0.10.1: -Version 0.10.1 Release Notes (current) --------------------------------------- +Version 0.10.1 Release Notes +---------------------------- * Released: 17 Aug 2024 * `GitHub release page diff --git a/doc/releases/0.10.2-notes.rst b/doc/releases/0.10.2-notes.rst new file mode 100644 index 000000000..175fdaff2 --- /dev/null +++ b/doc/releases/0.10.2-notes.rst @@ -0,0 +1,240 @@ +.. currentmodule:: control + +.. _version-0.10.2: + +Version 0.10.2 Release Notes (current) +-------------------------------------- + +* Released: date of release +* `GitHub release page + `_ + +This release includes numerous bug fixes and improvements, with major +changes such as a substantial reorganization of the documentation into +a User Guide and Reference Manual, more consistent and complete +docstrings, and support for referencing signals and subsystems by name +as well as by index. Phase plane plots now use matplotlib’s +`streamplot` for better visuals. New functions include `combine_tf` +and `split_tf` for MIMO/SISO conversion and `disk_margins` for +stability analysis. Additional improvements include consistent keyword +usage, expanded LTI system methods for plotting and responses, better +error messages, and legacy aliases to maintain backward compatibility. + +This version of `python-control` requires Python 3.10 or higher, NumPy +1.23 or higher (2.x recommended), and SciPy 1.8 or higher. + + +New classes, functions, and methods +................................... + +The following new classes, functions, and methods have been added in +this release: + +* `find_operating_point`: this function replaces (with a legacy alias) + the `find_eqpt` function and now returns an `OperatingPoint` object + containing the information about the operating point. + +* `combine_tf` and `split_tf`: these two new functions allow you to + create an MIMO transfer function from SISO transfer functions and + vice versa. + +* `create_statefbk_iosystem` now allows the creation of state feedback + controllers using a "reference gain" pattern (:math:`u = k_\text{f}\, + r - K x`) in addition to the default "trajectory generation" pattern + (:math:`u = u_\text{d} - K(x - x_\text{d})`). + +* `disk_margins`: compute disk-based stability margins for SISO and + MIMO systems. + +* `model_reduction`: allow specific states, inputs, or outputs to be + either eliminated or retained. + +* `place_acker`: renamed version of `acker` (which is still accessible + via an alias). + + +Bug fixes +......... + +The following bugs have been fixed in this release: + +* `phase_plane_plot`: fixed a bug in which the return value was + returning a sublist of lines rather than just a list of lines in + `cplt.lines`. + +* Processing of the timebase parameter (`dt`) for I/O systems is now + handled uniformly across all I/O system factory functions. This + affected the `zpk` function, which was defaulting to a discrete time + system to have timebase None instead of 0. + +* Multiplying (*), adding (+), or subtracting (-) a constant from any + (MIMO) LTI object now acts element-wise (same as ndarray's). This + fixes a bug where multiplying a MIMO LTI system by a constant was + multiplying by a matrix filled with the constant rather than a + diagonal matrix (scaled identity). + +* Fixed a bug where specifying an FRD system with fewer than 4 + frequency points was generating an error because the default + settings try to set up a smooth (interpolating) response and the + default degree of the fit was 3. + +* Fixed some bugs where computing poles and zeros of transfer + functions could generate spurious error messages about unsafe + casting of complex numbers to real numbers. + +* `TimeResponseData.to_pandas`: multi-trace data (e.g., the output + from a MIMO step response) was not being processed correctly. A new + column 'trace' is now generated for multi-trace responses. + +* Fixed a bug where where some arguments to `nyquist_plot` were not + being processed correctly and generated errors about unrecognized + keywords. + +* Updated `ctrb` and `obsv` to handle 1D `B` or `C` matrix correctly. + +* `bode_plot`: Fixed missing plot title when `display_margin` keyword + was used. + +* `singular_values_plot`: color cycling was not working correctly when + a list of systems or responses was provided. + +* `nyquist_plot`: The `lines` parameter of the `ControlPlot` object + now matches the documentation. A 2D array is returned with the + first index corresponding to the response (system) index and the + second index corresponding to the segment type (primary, mirror x + unscaled, scaled). + +* Fix some internal bugs that cropped up when using NumPy 2.3.1 but + were latent prior to that. + + +Improvements +............ + +The following additional improvements and changes in functionality +were implemented in this release: + +* User documentation is now divided into a User Guide that provides a + description of the main functionality of the python-control package, + along with a Reference Manual describing classes, functions, and + parameters in more detail. + +* Signal responses and I/O subsystem specifications can now use signal + names in addition to indices to get the desired inputs, outputs, and + states (e.g., `response.outputs['y0', 'y1']`). This is implemented + via a new `NamedSignal` object, which generalizes `numpy.ndarray`. + +* `find_operating_point` (legacy `find_eqpt`): accepts new parameters + `root_method` and `root_kwargs` to set the root finding algorithm + that is used. + +* `root_locus_map` now correctly handles the degenerate case of being + passed a single gain. + +* The `PoleZeroData` object now takes a `sort_loci` parameter when it + is created, with a default value of True. This is useful if you + create a `PoleZeroData` object by hand (e.g., for producing stability + diagrams). + +* Factory functions for I/O system creation are now consistent in + terms of copying signal/system names, overriding system/signal + names, and converting between classes. + +* The `tf` factory function to allow a 2D list of SISO transfer + functions to be given as a means of creating a MIMO transfer + function (use the new `combine_tf` function). + +* The `nlsys` factory function can now create a `NonlinearIOSystem` + representation of a `StateSpace` system (passed as the first + argument to `nlsys`). + +* LTI systems now have member functions for computing the various time + responses and generating frequency domain plots. See `LTI.to_ss`, + `LTI.to_tf`, `LTI.bode_plot`, `LTI.nyquist_plot`, `LTI.nichols_plot` + and `LTI.forced_response`, `LTI.impulse_response`, + `LTI.initial_response`, `LTI.step_response`. + +* String representations of I/O systems (accessed via `repr`, `print`, + and `str`) have been updated to create a more consistent form and + provide more useful information. See + :ref:`displaying-lti-system-information` for more information. + +* Binary operations between MIMO and SISO functions are now supported, + with the SISO system being converted to a MIMO system as if it were + a scalar. + +* `nyquist_response`: generates an error if you force the system to + evaluate the dynamics at a pole. + +* `phase_crossover_frequencies`: turned off spurious warning messages. + +* `ss2tf`: added new `method=scipy` capability, allowing `ss2tf` to + work on MIMO systems even if Slycot is not present. + +* `flatsys.solve_flat_optimal` (legacy `flatsys.solve_flat_ocp`): + allows scalar time vector. + +* Improved checking of matrix shapes and better error messages in + state space factory functions and other operations where matrices + are passed as arguments. + +* `FrequencyResponseData`: use `~FrequencyResponseData.complex` to + access the (squeeze processed) complex frequency response (instead + of the legacy `response` property) and + `~FrequencyResponseData.frdata ` to access + the 3D frequency response data array (instead of the legacy `fresp` + attribute). + +* Time response and optimization function keywords have been + regularized to allow consistent use of keywords across related + functions: + + - Parameters specifying the inputs, outputs, and states are referred + to as `inputs`, `outputs`, and `states` consistently throughout the + functions. + + - Variables associated with inputs, outputs, states and time use + those words plus an appropriate modifier: `initial_state`, + `final_output`, `input_indices`, etc. + + - Aliases are used both to maintain backward compatibility and to + allow shorthand descriptions: e.g., `U`, `Y`, `X0`. Short form + aliases are documented in docstrings by listing the parameter as + ``long_form (or sf) : type``. + + - Existing legacy keywords are allowed and generate a + `PendingDeprecationWarning`. Specifying a parameter value in two + different ways (e.g., via long form and an alias) generates a + `TypeError`. + +* `phase_plane_plot`: makes use of the matplotlib + `~matplotlib.pyplot.streamplot` function to provide better default + phase plane diagrams. + +* `root_locus_plot`: added by the ability to recompute the root locus + when zooming in on portions of the root locus diagram. + +* `nyquist_plot`: updated the rescaling algorithm to use a more + gradual change in the magnitude of the Nyquist curve. The + `blend_fraction` parameter can be used to start the rescaling prior + to reaching `max_curve_magnitude`, giving less confusing plots. Some + default parameter values have been adjusted to improve Nyquist + plots. + + +Deprecations +............ + +The following functions have been newly deprecated in this release and +generate a warning message when used: + +* `FrequencyResponseData.response`: use + `FrequencyResponseData.complex` to return the complex value of the + frequency response. + +* `FrequencyResponseData.fresp`: use `FrequencyResponseData.frdata + ` to access the raw 3D frequency response + data array. + +The listed items are slated to be removed in future releases (usually +the next major or minor version update). diff --git a/doc/requirements.txt b/doc/requirements.txt index 5fdf9113d..ded0c7087 100644 --- a/doc/requirements.txt +++ b/doc/requirements.txt @@ -7,4 +7,4 @@ sphinx-copybutton numpydoc ipykernel nbsphinx -docutils==0.16 # pin until sphinx_rtd_theme is compatible with 0.17 or later +docutils diff --git a/doc/response.rst b/doc/response.rst index 0058a500d..8ccdccba8 100644 --- a/doc/response.rst +++ b/doc/response.rst @@ -547,7 +547,7 @@ the computation of the Nyquist curve and the way the data are plotted: sys = ct.tf([1, 0.2], [1, 0, 1]) * ct.tf([1], [1, 0]) nyqresp = ct.nyquist_response(sys) nyqresp.plot( - max_curve_magnitude=6, max_curve_offset=1, + max_curve_magnitude=6, max_curve_offset=1, blend_fraction=0.05, arrows=[0, 0.15, 0.3, 0.6, 0.7, 0.925], title='Custom Nyquist plot') print("Encirclements =", nyqresp.count) diff --git a/doc/test_sphinxdocs.py b/doc/test_sphinxdocs.py index 1a49f357c..80024c23d 100644 --- a/doc/test_sphinxdocs.py +++ b/doc/test_sphinxdocs.py @@ -48,7 +48,6 @@ # Functons that we can skip object_skiplist = [ - control.NamedSignal, # np.ndarray members cause errors control.FrequencyResponseList, # Use FrequencyResponseData control.TimeResponseList, # Use TimeResponseData control.common_timebase, # mainly internal use diff --git a/examples/bdalg-matlab.py b/examples/bdalg-matlab.py index 8911d6579..eaafaa59a 100644 --- a/examples/bdalg-matlab.py +++ b/examples/bdalg-matlab.py @@ -1,7 +1,7 @@ -# bdalg-matlab.py - demonstrate some MATLAB commands for block diagram altebra +# bdalg-matlab.py - demonstrate some MATLAB commands for block diagram algebra # RMM, 29 May 09 -from control.matlab import * # MATLAB-like functions +from control.matlab import ss, ss2tf, tf, tf2ss # MATLAB-like functions # System matrices A1 = [[0, 1.], [-4, -1]] diff --git a/examples/cds110-L8b_pvtol-complete-limits.ipynb b/examples/cds110-L8b_pvtol-complete-limits.ipynb index 0b482c865..b0ae80d76 100644 --- a/examples/cds110-L8b_pvtol-complete-limits.ipynb +++ b/examples/cds110-L8b_pvtol-complete-limits.ipynb @@ -830,13 +830,13 @@ "\n", "# Gain margin for Lx\n", "neg1overgm_x = -0.67 # vary this manually to find intersection with curve\n", - "color = cplt.lines[0][0].get_color()\n", + "color = cplt.lines[0, 0][0].get_color()\n", "plt.plot(neg1overgm_x, 0, color=color, marker='o', fillstyle='none')\n", "gm_x = -1/neg1overgm_x\n", "\n", "# Gain margin for Ly\n", "neg1overgm_y = -0.32 # vary this manually to find intersection with curve\n", - "color = cplt.lines[1][0].get_color()\n", + "color = cplt.lines[1, 0][0].get_color()\n", "plt.plot(neg1overgm_y, 0, color=color, marker='o', fillstyle='none')\n", "gm_y = -1/neg1overgm_y\n", "\n", @@ -885,13 +885,13 @@ "# Phase margin of Lx:\n", "th_pm_x = 0.14*np.pi\n", "th_plt_x = np.pi + th_pm_x\n", - "color = cplt.lines[0][0].get_color()\n", + "color = cplt.lines[0, 0][0].get_color()\n", "plt.plot(np.cos(th_plt_x), np.sin(th_plt_x), color=color, marker='o')\n", "\n", "# Phase margin of Ly\n", "th_pm_y = 0.19*np.pi\n", "th_plt_y = np.pi + th_pm_y\n", - "color = cplt.lines[1][0].get_color()\n", + "color = cplt.lines[1, 0][0].get_color()\n", "plt.plot(np.cos(th_plt_y), np.sin(th_plt_y), color=color, marker='o')\n", "\n", "print('Margins obtained visually:')\n", @@ -936,12 +936,12 @@ "\n", "# Stability margin:\n", "sm_x = 0.3 # vary this manually to find min which intersects\n", - "color = cplt.lines[0][0].get_color()\n", + "color = cplt.lines[0, 0][0].get_color()\n", "sm_circle = plt.Circle((-1, 0), sm_x, color=color, fill=False, ls=':')\n", "cplt.axes[0, 0].add_patch(sm_circle)\n", "\n", "sm_y = 0.5 # vary this manually to find min which intersects\n", - "color = cplt.lines[1][0].get_color()\n", + "color = cplt.lines[1, 0][0].get_color()\n", "sm_circle = plt.Circle((-1, 0), sm_y, color=color, fill=False, ls=':')\n", "cplt.axes[0, 0].add_patch(sm_circle)\n", "\n", @@ -1023,8 +1023,22 @@ } ], "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, "language_info": { - "name": "python" + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.13.1" } }, "nbformat": 4, diff --git a/examples/check-controllability-and-observability.py b/examples/check-controllability-and-observability.py index 67ecdf26c..a8fc5c6ad 100644 --- a/examples/check-controllability-and-observability.py +++ b/examples/check-controllability-and-observability.py @@ -4,8 +4,8 @@ RMM, 6 Sep 2010 """ -import numpy as np # Load the scipy functions -from control.matlab import * # Load the controls systems library +import numpy as np # Load the numpy functions +from control.matlab import ss, ctrb, obsv # Load the controls systems library # Parameters defining the system diff --git a/examples/cruise-control.py b/examples/cruise-control.py index 5bb263830..77768aa86 100644 --- a/examples/cruise-control.py +++ b/examples/cruise-control.py @@ -247,7 +247,6 @@ def pi_update(t, x, u, params={}): # Assign variables for inputs and states (for readability) v = u[0] # current velocity vref = u[1] # reference velocity - z = x[0] # integrated error # Compute the nominal controller output (needed for anti-windup) u_a = pi_output(t, x, u, params) @@ -394,7 +393,7 @@ def sf_output(t, z, u, params={}): ud = params.get('ud', 0) # Get the system state and reference input - x, y, r = u[0], u[1], u[2] + x, r = u[0], u[2] return ud - K * (x - xd) - ki * z + kf * (r - yd) @@ -440,13 +439,13 @@ def sf_output(t, z, u, params={}): 4./180. * pi for t in T] t, y = ct.input_output_response( cruise_sf, T, [vref, gear, theta_hill], [X0[0], 0], - params={'K': K, 'kf': kf, 'ki': 0.0, 'kf': kf, 'xd': xd, 'ud': ud, 'yd': yd}) + params={'K': K, 'kf': kf, 'ki': 0.0, 'xd': xd, 'ud': ud, 'yd': yd}) subplots = cruise_plot(cruise_sf, t, y, label='Proportional', linetype='b--') # Response of the system with state feedback + integral action t, y = ct.input_output_response( cruise_sf, T, [vref, gear, theta_hill], [X0[0], 0], - params={'K': K, 'kf': kf, 'ki': 0.1, 'kf': kf, 'xd': xd, 'ud': ud, 'yd': yd}) + params={'K': K, 'kf': kf, 'ki': 0.1, 'xd': xd, 'ud': ud, 'yd': yd}) cruise_plot(cruise_sf, t, y, label='PI control', t_hill=8, linetype='b-', subplots=subplots, legend=True) diff --git a/examples/cruise.ipynb b/examples/cruise.ipynb index 16935b15e..08a1583ac 100644 --- a/examples/cruise.ipynb +++ b/examples/cruise.ipynb @@ -420,8 +420,8 @@ "name": "stdout", "output_type": "stream", "text": [ - "system: a = (0.010124405669387215-0j) , b = (1.3203061238159202+0j)\n", - "pzcancel: kp = 0.5 , ki = (0.005062202834693608+0j) , 1/(kp b) = (1.5148002148317266+0j)\n", + "system: a = 0.010124405669387215 , b = 1.3203061238159202\n", + "pzcancel: kp = 0.5 , ki = 0.005062202834693608 , 1/(kp b) = 1.5148002148317266\n", "sfb_int: K = 0.5 , ki = 0.1\n" ] }, @@ -442,7 +442,7 @@ "\n", "# Construction a controller that cancels the pole\n", "kp = 0.5\n", - "a = -P.poles()[0]\n", + "a = -P.poles()[0].real\n", "b = np.real(P(0)) * a\n", "ki = a * kp\n", "control_pz = ct.TransferFunction(\n", diff --git a/examples/disk_margins.py b/examples/disk_margins.py new file mode 100644 index 000000000..1b9934156 --- /dev/null +++ b/examples/disk_margins.py @@ -0,0 +1,548 @@ +"""disk_margins.py + +Demonstrate disk-based stability margin calculations. + +References: +[1] Blight, James D., R. Lane Dailey, and Dagfinn Gangsaas. “Practical + Control Law Design for Aircraft Using Multivariable Techniques.” + International Journal of Control 59, no. 1 (January 1994): 93-137. + https://doi.org/10.1080/00207179408923071. + +[2] Seiler, Peter, Andrew Packard, and Pascal Gahinet. “An Introduction + to Disk Margins [Lecture Notes].” IEEE Control Systems Magazine 40, + no. 5 (October 2020): 78-95. + +[3] P. Benner, V. Mehrmann, V. Sima, S. Van Huffel, and A. Varga, "SLICOT + - A Subroutine Library in Systems and Control Theory", Applied and + Computational Control, Signals, and Circuits (Birkhauser), Vol. 1, Ch. + 10, pp. 505-546, 1999. + +[4] S. Van Huffel, V. Sima, A. Varga, S. Hammarling, and F. Delebecque, + "Development of High Performance Numerical Software for Control", IEEE + Control Systems Magazine, Vol. 24, Nr. 1, Feb., pp. 60-76, 2004. + +[5] Deodhare, G., & Patel, V. (1998, August). A "Modern" Look at Gain + and Phase Margins: An H-Infinity/mu Approach. In Guidance, Navigation, + and Control Conference and Exhibit (p. 4134). +""" + +import os +import control +import matplotlib.pyplot as plt +import numpy as np + +def plot_allowable_region(alpha_max, skew, ax=None): + """Plot region of allowable gain/phase variation, given worst-case disk margin. + + Parameters + ---------- + alpha_max : float (scalar or list) + worst-case disk margin(s) across all frequencies. May be a scalar or list. + skew : float (scalar or list) + skew parameter(s) for disk margin calculation. + skew=0 uses the "balanced" sensitivity function 0.5*(S - T) + skew=1 uses the sensitivity function S + skew=-1 uses the complementary sensitivity function T + ax : axes to plot bounding curve(s) onto + + Returns + ------- + DM : ndarray + 1D array of frequency-dependent disk margins. DM is the same + size as "omega" parameter. + GM : ndarray + 1D array of frequency-dependent disk-based gain margins, in dB. + GM is the same size as "omega" parameter. + PM : ndarray + 1D array of frequency-dependent disk-based phase margins, in deg. + PM is the same size as "omega" parameter. + """ + + # Create axis if needed + if ax is None: + ax = plt.gca() + + # Allow scalar or vector arguments (to overlay plots) + if np.isscalar(alpha_max): + alpha_max = np.asarray([alpha_max]) + else: + alpha_max = np.asarray(alpha_max) + + if np.isscalar(skew): + skew=np.asarray([skew]) + else: + skew=np.asarray(skew) + + # Add a plot for each (alpha, skew) pair present + theta = np.linspace(0, np.pi, 500) + legend_list = [] + for ii in range(0, skew.shape[0]): + legend_str = "$\\sigma$ = %.1f, $\\alpha_{max}$ = %.2f" %(\ + skew[ii], alpha_max[ii]) + legend_list.append(legend_str) + + # Complex bounding curve of stable gain/phase variations + f = (2 + alpha_max[ii] * (1 - skew[ii]) * np.exp(1j * theta))\ + /(2 - alpha_max[ii] * (1 + skew[ii]) * np.exp(1j * theta)) + + # Allowable combined gain/phase variations + gamma_dB = control.ctrlutil.mag2db(np.abs(f)) # gain margin (dB) + phi_deg = np.rad2deg(np.angle(f)) # phase margin (deg) + + # Plot the allowable combined gain/phase variations + out = ax.plot(gamma_dB, phi_deg, alpha=0.25, label='_nolegend_') + ax.fill_between(ax.lines[ii].get_xydata()[:,0],\ + ax.lines[ii].get_xydata()[:,1], alpha=0.25) + + plt.ylabel('Phase Variation (deg)') + plt.xlabel('Gain Variation (dB)') + plt.title('Range of Gain and Phase Variations') + plt.legend(legend_list) + plt.grid() + plt.tight_layout() + + return out + +def test_siso1(): + # + # Disk-based stability margins for example + # SISO loop transfer function(s) + # + + # Frequencies of interest + omega = np.logspace(-1, 2, 1001) + + # Loop transfer gain + L = control.tf(25, [1, 10, 10, 10]) + + print("------------- Python control built-in (S) -------------") + GM_, PM_, SM_ = control.stability_margins(L)[:3] # python-control default (S-based...?) + print(f"SM_ = {SM_}") + print(f"GM_ = {GM_} dB") + print(f"PM_ = {PM_} deg\n") + + print("------------- Sensitivity function (S) -------------") + DM, GM, PM = control.disk_margins(L, omega, skew=1.0, returnall=True) # S-based (S) + print(f"min(DM) = {min(DM)} (omega = {omega[np.argmin(DM)]})") + print(f"GM = {GM[np.argmin(DM)]} dB") + print(f"PM = {PM[np.argmin(DM)]} deg") + print(f"min(GM) = {min(GM)} dB") + print(f"min(PM) = {min(PM)} deg\n") + + plt.figure(1) + plt.subplot(3, 3, 1) + plt.semilogx(omega, DM, label='$\\alpha$') + plt.ylabel('Disk Margin (abs)') + plt.legend() + plt.title('S-Based Margins') + plt.grid() + plt.xlim([omega[0], omega[-1]]) + plt.ylim([0, 2]) + + plt.figure(1) + plt.subplot(3, 3, 4) + plt.semilogx(omega, GM, label='$\\gamma_{m}$') + plt.ylabel('Gain Margin (dB)') + plt.legend() + #plt.title('Gain-Only Margin') + plt.grid() + plt.xlim([omega[0], omega[-1]]) + plt.ylim([0, 40]) + + plt.figure(1) + plt.subplot(3, 3, 7) + plt.semilogx(omega, PM, label='$\\phi_{m}$') + plt.ylabel('Phase Margin (deg)') + plt.legend() + #plt.title('Phase-Only Margin') + plt.grid() + plt.xlim([omega[0], omega[-1]]) + plt.ylim([0, 90]) + plt.xlabel('Frequency (rad/s)') + + print("------------- Complementary sensitivity function (T) -------------") + DM, GM, PM = control.disk_margins(L, omega, skew=-1.0, returnall=True) # T-based (T) + print(f"min(DM) = {min(DM)} (omega = {omega[np.argmin(DM)]})") + print(f"GM = {GM[np.argmin(DM)]} dB") + print(f"PM = {PM[np.argmin(DM)]} deg") + print(f"min(GM) = {min(GM)} dB") + print(f"min(PM) = {min(PM)} deg\n") + + plt.figure(1) + plt.subplot(3, 3, 2) + plt.semilogx(omega, DM, label='$\\alpha$') + plt.ylabel('Disk Margin (abs)') + plt.legend() + plt.title('T_Based Margins') + plt.grid() + plt.xlim([omega[0], omega[-1]]) + plt.ylim([0, 2]) + + plt.figure(1) + plt.subplot(3, 3, 5) + plt.semilogx(omega, GM, label='$\\gamma_{m}$') + plt.ylabel('Gain Margin (dB)') + plt.legend() + #plt.title('Gain-Only Margin') + plt.grid() + plt.xlim([omega[0], omega[-1]]) + plt.ylim([0, 40]) + + plt.figure(1) + plt.subplot(3, 3, 8) + plt.semilogx(omega, PM, label='$\\phi_{m}$') + plt.ylabel('Phase Margin (deg)') + plt.legend() + #plt.title('Phase-Only Margin') + plt.grid() + plt.xlim([omega[0], omega[-1]]) + plt.ylim([0, 90]) + plt.xlabel('Frequency (rad/s)') + + print("------------- Balanced sensitivity function (S - T) -------------") + DM, GM, PM = control.disk_margins(L, omega, skew=0.0, returnall=True) # balanced (S - T) + print(f"min(DM) = {min(DM)} (omega = {omega[np.argmin(DM)]})") + print(f"GM = {GM[np.argmin(DM)]} dB") + print(f"PM = {PM[np.argmin(DM)]} deg") + print(f"min(GM) = {min(GM)} dB") + print(f"min(PM) = {min(PM)} deg\n") + + plt.figure(1) + plt.subplot(3, 3, 3) + plt.semilogx(omega, DM, label='$\\alpha$') + plt.ylabel('Disk Margin (abs)') + plt.legend() + plt.title('Balanced Margins') + plt.grid() + plt.xlim([omega[0], omega[-1]]) + plt.ylim([0, 2]) + + plt.figure(1) + plt.subplot(3, 3, 6) + plt.semilogx(omega, GM, label='$\\gamma_{m}$') + plt.ylabel('Gain Margin (dB)') + plt.legend() + #plt.title('Gain-Only Margin') + plt.grid() + plt.xlim([omega[0], omega[-1]]) + plt.ylim([0, 40]) + + plt.figure(1) + plt.subplot(3, 3, 9) + plt.semilogx(omega, PM, label='$\\phi_{m}$') + plt.ylabel('Phase Margin (deg)') + plt.legend() + #plt.title('Phase-Only Margin') + plt.grid() + plt.xlim([omega[0], omega[-1]]) + plt.ylim([0, 90]) + plt.xlabel('Frequency (rad/s)') + + # Disk margin plot of admissible gain/phase variations for which + DM_plot = [] + DM_plot.append(control.disk_margins(L, omega, skew=-2.0)[0]) + DM_plot.append(control.disk_margins(L, omega, skew=0.0)[0]) + DM_plot.append(control.disk_margins(L, omega, skew=2.0)[0]) + plt.figure(10); plt.clf() + plot_allowable_region(DM_plot, skew=[-2.0, 0.0, 2.0]) + + return + +def test_siso2(): + # + # Disk-based stability margins for example + # SISO loop transfer function(s) + # + + # Frequencies of interest + omega = np.logspace(-1, 2, 1001) + + # Laplace variable + s = control.tf('s') + + # Loop transfer gain + L = (6.25 * (s + 3) * (s + 5)) / (s * (s + 1)**2 * (s**2 + 0.18 * s + 100)) + + print("------------- Python control built-in (S) -------------") + GM_, PM_, SM_ = control.stability_margins(L)[:3] # python-control default (S-based...?) + print(f"SM_ = {SM_}") + print(f"GM_ = {GM_} dB") + print(f"PM_ = {PM_} deg\n") + + print("------------- Sensitivity function (S) -------------") + DM, GM, PM = control.disk_margins(L, omega, skew=1.0, returnall=True) # S-based (S) + print(f"min(DM) = {min(DM)} (omega = {omega[np.argmin(DM)]})") + print(f"GM = {GM[np.argmin(DM)]} dB") + print(f"PM = {PM[np.argmin(DM)]} deg") + print(f"min(GM) = {min(GM)} dB") + print(f"min(PM) = {min(PM)} deg\n") + + plt.figure(2) + plt.subplot(3, 3, 1) + plt.semilogx(omega, DM, label='$\\alpha$') + plt.ylabel('Disk Margin (abs)') + plt.legend() + plt.title('S-Based Margins') + plt.grid() + plt.xlim([omega[0], omega[-1]]) + plt.ylim([0, 2]) + + plt.figure(2) + plt.subplot(3, 3, 4) + plt.semilogx(omega, GM, label='$\\gamma_{m}$') + plt.ylabel('Gain Margin (dB)') + plt.legend() + #plt.title('Gain-Only Margin') + plt.grid() + plt.xlim([omega[0], omega[-1]]) + plt.ylim([0, 40]) + + plt.figure(2) + plt.subplot(3, 3, 7) + plt.semilogx(omega, PM, label='$\\phi_{m}$') + plt.ylabel('Phase Margin (deg)') + plt.legend() + #plt.title('Phase-Only Margin') + plt.grid() + plt.xlim([omega[0], omega[-1]]) + plt.ylim([0, 90]) + plt.xlabel('Frequency (rad/s)') + + print("------------- Complementary sensitivity function (T) -------------") + DM, GM, PM = control.disk_margins(L, omega, skew=-1.0, returnall=True) # T-based (T) + print(f"min(DM) = {min(DM)} (omega = {omega[np.argmin(DM)]})") + print(f"GM = {GM[np.argmin(DM)]} dB") + print(f"PM = {PM[np.argmin(DM)]} deg") + print(f"min(GM) = {min(GM)} dB") + print(f"min(PM) = {min(PM)} deg\n") + + plt.figure(2) + plt.subplot(3, 3, 2) + plt.semilogx(omega, DM, label='$\\alpha$') + plt.ylabel('Disk Margin (abs)') + plt.legend() + plt.title('T-Based Margins') + plt.grid() + plt.xlim([omega[0], omega[-1]]) + plt.ylim([0, 2]) + + plt.figure(2) + plt.subplot(3, 3, 5) + plt.semilogx(omega, GM, label='$\\gamma_{m}$') + plt.ylabel('Gain Margin (dB)') + plt.legend() + #plt.title('Gain-Only Margin') + plt.grid() + plt.xlim([omega[0], omega[-1]]) + plt.ylim([0, 40]) + + plt.figure(2) + plt.subplot(3, 3, 8) + plt.semilogx(omega, PM, label='$\\phi_{m}$') + plt.ylabel('Phase Margin (deg)') + plt.legend() + #plt.title('Phase-Only Margin') + plt.grid() + plt.xlim([omega[0], omega[-1]]) + plt.ylim([0, 90]) + plt.xlabel('Frequency (rad/s)') + + print("------------- Balanced sensitivity function (S - T) -------------") + DM, GM, PM = control.disk_margins(L, omega, skew=0.0, returnall=True) # balanced (S - T) + print(f"min(DM) = {min(DM)} (omega = {omega[np.argmin(DM)]})") + print(f"GM = {GM[np.argmin(DM)]} dB") + print(f"PM = {PM[np.argmin(DM)]} deg") + print(f"min(GM) = {min(GM)} dB") + print(f"min(PM) = {min(PM)} deg\n") + + plt.figure(2) + plt.subplot(3, 3, 3) + plt.semilogx(omega, DM, label='$\\alpha$') + plt.ylabel('Disk Margin (abs)') + plt.legend() + plt.title('Balanced Margins') + plt.grid() + plt.xlim([omega[0], omega[-1]]) + plt.ylim([0, 2]) + + plt.figure(2) + plt.subplot(3, 3, 6) + plt.semilogx(omega, GM, label='$\\gamma_{m}$') + plt.ylabel('Gain Margin (dB)') + plt.legend() + #plt.title('Gain-Only Margin') + plt.grid() + plt.xlim([omega[0], omega[-1]]) + plt.ylim([0, 40]) + + plt.figure(2) + plt.subplot(3, 3, 9) + plt.semilogx(omega, PM, label='$\\phi_{m}$') + plt.ylabel('Phase Margin (deg)') + plt.legend() + #plt.title('Phase-Only Margin') + plt.grid() + plt.xlim([omega[0], omega[-1]]) + plt.ylim([0, 90]) + plt.xlabel('Frequency (rad/s)') + + # Disk margin plot of admissible gain/phase variations for which + # the feedback loop still remains stable, for each skew parameter + DM_plot = [] + DM_plot.append(control.disk_margins(L, omega, skew=-1.0)[0]) # T-based (T) + DM_plot.append(control.disk_margins(L, omega, skew=0.0)[0]) # balanced (S - T) + DM_plot.append(control.disk_margins(L, omega, skew=1.0)[0]) # S-based (S) + plt.figure(20) + plot_allowable_region(DM_plot, skew=[-1.0, 0.0, 1.0]) + + return + +def test_mimo(): + # + # Disk-based stability margins for example + # MIMO loop transfer function(s) + # + + # Frequencies of interest + omega = np.logspace(-1, 3, 1001) + + # Loop transfer gain + P = control.ss([[0, 10],[-10, 0]], np.eye(2), [[1, 10], [-10, 1]], 0) # plant + K = control.ss([], [], [], [[1, -2], [0, 1]]) # controller + L = P * K # loop gain + + print("------------- Sensitivity function (S) -------------") + DM, GM, PM = control.disk_margins(L, omega, skew=1.0, returnall=True) # S-based (S) + print(f"min(DM) = {min(DM)} (omega = {omega[np.argmin(DM)]})") + print(f"GM = {GM[np.argmin(DM)]} dB") + print(f"PM = {PM[np.argmin(DM)]} deg") + print(f"min(GM) = {min(GM)} dB") + print(f"min(PM) = {min(PM)} deg\n") + + plt.figure(3) + plt.subplot(3, 3, 1) + plt.semilogx(omega, DM, label='$\\alpha$') + plt.ylabel('Disk Margin (abs)') + plt.legend() + plt.title('S-Based Margins') + plt.grid() + plt.xlim([omega[0], omega[-1]]) + plt.ylim([0, 2]) + + plt.figure(3) + plt.subplot(3, 3, 4) + plt.semilogx(omega, GM, label='$\\gamma_{m}$') + plt.ylabel('Gain Margin (dB)') + plt.legend() + #plt.title('Gain-Only Margin') + plt.grid() + plt.xlim([omega[0], omega[-1]]) + plt.ylim([0, 40]) + + plt.figure(3) + plt.subplot(3, 3, 7) + plt.semilogx(omega, PM, label='$\\phi_{m}$') + plt.ylabel('Phase Margin (deg)') + plt.legend() + #plt.title('Phase-Only Margin') + plt.grid() + plt.xlim([omega[0], omega[-1]]) + plt.ylim([0, 90]) + plt.xlabel('Frequency (rad/s)') + + print("------------- Complementary sensitivity function (T) -------------") + DM, GM, PM = control.disk_margins(L, omega, skew=-1.0, returnall=True) # T-based (T) + print(f"min(DM) = {min(DM)} (omega = {omega[np.argmin(DM)]})") + print(f"GM = {GM[np.argmin(DM)]} dB") + print(f"PM = {PM[np.argmin(DM)]} deg") + print(f"min(GM) = {min(GM)} dB") + print(f"min(PM) = {min(PM)} deg\n") + + plt.figure(3) + plt.subplot(3, 3, 2) + plt.semilogx(omega, DM, label='$\\alpha$') + plt.ylabel('Disk Margin (abs)') + plt.legend() + plt.title('T-Based Margins') + plt.grid() + plt.xlim([omega[0], omega[-1]]) + plt.ylim([0, 2]) + + plt.figure(3) + plt.subplot(3, 3, 5) + plt.semilogx(omega, GM, label='$\\gamma_{m}$') + plt.ylabel('Gain Margin (dB)') + plt.legend() + #plt.title('Gain-Only Margin') + plt.grid() + plt.xlim([omega[0], omega[-1]]) + plt.ylim([0, 40]) + + plt.figure(3) + plt.subplot(3, 3, 8) + plt.semilogx(omega, PM, label='$\\phi_{m}$') + plt.ylabel('Phase Margin (deg)') + plt.legend() + #plt.title('Phase-Only Margin') + plt.grid() + plt.xlim([omega[0], omega[-1]]) + plt.ylim([0, 90]) + plt.xlabel('Frequency (rad/s)') + + print("------------- Balanced sensitivity function (S - T) -------------") + DM, GM, PM = control.disk_margins(L, omega, skew=0.0, returnall=True) # balanced (S - T) + print(f"min(DM) = {min(DM)} (omega = {omega[np.argmin(DM)]})") + print(f"GM = {GM[np.argmin(DM)]} dB") + print(f"PM = {PM[np.argmin(DM)]} deg") + print(f"min(GM) = {min(GM)} dB") + print(f"min(PM) = {min(PM)} deg\n") + + plt.figure(3) + plt.subplot(3, 3, 3) + plt.semilogx(omega, DM, label='$\\alpha$') + plt.ylabel('Disk Margin (abs)') + plt.legend() + plt.title('Balanced Margins') + plt.grid() + plt.xlim([omega[0], omega[-1]]) + plt.ylim([0, 2]) + + plt.figure(3) + plt.subplot(3, 3, 6) + plt.semilogx(omega, GM, label='$\\gamma_{m}$') + plt.ylabel('Gain Margin (dB)') + plt.legend() + #plt.title('Gain-Only Margin') + plt.grid() + plt.xlim([omega[0], omega[-1]]) + plt.ylim([0, 40]) + + plt.figure(3) + plt.subplot(3, 3, 9) + plt.semilogx(omega, PM, label='$\\phi_{m}$') + plt.ylabel('Phase Margin (deg)') + plt.legend() + #plt.title('Phase-Only Margin') + plt.grid() + plt.xlim([omega[0], omega[-1]]) + plt.ylim([0, 90]) + plt.xlabel('Frequency (rad/s)') + + # Disk margin plot of admissible gain/phase variations for which + # the feedback loop still remains stable, for each skew parameter + DM_plot = [] + DM_plot.append(control.disk_margins(L, omega, skew=-1.0)[0]) # T-based (T) + DM_plot.append(control.disk_margins(L, omega, skew=0.0)[0]) # balanced (S - T) + DM_plot.append(control.disk_margins(L, omega, skew=1.0)[0]) # S-based (S) + plt.figure(30) + plot_allowable_region(DM_plot, skew=[-1.0, 0.0, 1.0]) + + return + +if __name__ == '__main__': + #test_siso1() + #test_siso2() + test_mimo() + if 'PYCONTROL_TEST_EXAMPLES' not in os.environ: + #plt.tight_layout() + plt.show() diff --git a/examples/kincar.py b/examples/kincar.py index a12cdc774..ab026cba6 100644 --- a/examples/kincar.py +++ b/examples/kincar.py @@ -3,7 +3,6 @@ import numpy as np import matplotlib.pyplot as plt -import control as ct import control.flatsys as fs # diff --git a/examples/mrac_siso_mit.py b/examples/mrac_siso_mit.py index f8940e694..a821b65d0 100644 --- a/examples/mrac_siso_mit.py +++ b/examples/mrac_siso_mit.py @@ -46,7 +46,6 @@ def adaptive_controller_state(t, xc, uc, params): # Parameters gam = params["gam"] Am = params["Am"] - Bm = params["Bm"] signB = params["signB"] # Controller inputs diff --git a/examples/phase_plane_plots.py b/examples/phase_plane_plots.py index db989d5d9..44a47a29c 100644 --- a/examples/phase_plane_plots.py +++ b/examples/phase_plane_plots.py @@ -5,9 +5,8 @@ # using the phaseplot module. Most of these figures line up with examples # in FBS2e, with different display options shown as different subplots. -import time import warnings -from math import pi, sqrt +from math import pi import matplotlib.pyplot as plt import numpy as np diff --git a/examples/pvtol-nested-ss.py b/examples/pvtol-nested-ss.py index f53ac70f1..e8542a828 100644 --- a/examples/pvtol-nested-ss.py +++ b/examples/pvtol-nested-ss.py @@ -10,7 +10,6 @@ import os import matplotlib.pyplot as plt # MATLAB plotting functions -from control.matlab import * # MATLAB-like functions import numpy as np import math import control as ct @@ -23,12 +22,12 @@ c = 0.05 # damping factor (estimated) # Transfer functions for dynamics -Pi = tf([r], [J, 0, 0]) # inner loop (roll) -Po = tf([1], [m, c, 0]) # outer loop (position) +Pi = ct.tf([r], [J, 0, 0]) # inner loop (roll) +Po = ct.tf([1], [m, c, 0]) # outer loop (position) # Use state space versions -Pi = tf2ss(Pi) -Po = tf2ss(Po) +Pi = ct.tf2ss(Pi) +Po = ct.tf2ss(Po) # # Inner loop control design @@ -40,10 +39,10 @@ # Design a simple lead controller for the system k, a, b = 200, 2, 50 -Ci = k*tf([1, a], [1, b]) # lead compensator +Ci = k*ct.tf([1, a], [1, b]) # lead compensator # Convert to statespace -Ci = tf2ss(Ci) +Ci = ct.tf2ss(Ci) # Compute the loop transfer function for the inner loop Li = Pi*Ci @@ -51,49 +50,49 @@ # Bode plot for the open loop process plt.figure(1) -bode(Pi) +ct.bode(Pi) # Bode plot for the loop transfer function, with margins plt.figure(2) -bode(Li) +ct.bode(Li) # Compute out the gain and phase margins #! Not implemented # (gm, pm, wcg, wcp) = margin(Li); # Compute the sensitivity and complementary sensitivity functions -Si = feedback(1, Li) +Si = ct.feedback(1, Li) Ti = Li*Si # Check to make sure that the specification is met plt.figure(3) -gangof4(Pi, Ci) +ct.gangof4(Pi, Ci) # Compute out the actual transfer function from u1 to v1 (see L8.2 notes) # Hi = Ci*(1-m*g*Pi)/(1+Ci*Pi); -Hi = parallel(feedback(Ci, Pi), -m*g*feedback(Ci*Pi, 1)) +Hi = ct.parallel(ct.feedback(Ci, Pi), -m*g*ct.feedback(Ci*Pi, 1)) plt.figure(4) plt.clf() -bode(Hi) +ct.bode(Hi) # Now design the lateral control system a, b, K = 0.02, 5, 2 -Co = -K*tf([1, 0.3], [1, 10]) # another lead compensator +Co = -K*ct.tf([1, 0.3], [1, 10]) # another lead compensator # Convert to statespace -Co = tf2ss(Co) +Co = ct.tf2ss(Co) # Compute the loop transfer function for the outer loop Lo = -m*g*Po*Co plt.figure(5) -bode(Lo, display_margins=True) # margin(Lo) +ct.bode(Lo, display_margins=True) # margin(Lo) # Finally compute the real outer-loop loop gain + responses L = Co*Hi*Po -S = feedback(1, L) -T = feedback(L, 1) +S = ct.feedback(1, L) +T = ct.feedback(L, 1) # Compute stability margins #! Not yet implemented @@ -101,7 +100,7 @@ plt.figure(6) plt.clf() -out = ct.bode(L, logspace(-4, 3), initial_phase=-math.pi/2) +out = ct.bode(L, np.logspace(-4, 3), initial_phase=-math.pi/2) axs = ct.get_plot_axes(out) # Add crossover line to magnitude plot @@ -111,7 +110,7 @@ # Nyquist plot for complete design # plt.figure(7) -nyquist(L) +ct.nyquist(L) # set up the color color = 'b' @@ -126,10 +125,10 @@ # 'EdgeColor', color, 'FaceColor', color); plt.figure(9) -Yvec, Tvec = step(T, linspace(1, 20)) +Yvec, Tvec = ct.step_response(T, np.linspace(1, 20)) plt.plot(Tvec.T, Yvec.T) -Yvec, Tvec = step(Co*S, linspace(1, 20)) +Yvec, Tvec = ct.step_response(Co*S, np.linspace(1, 20)) plt.plot(Tvec.T, Yvec.T) #TODO: PZmap for statespace systems has not yet been implemented. @@ -142,7 +141,7 @@ # Gang of Four plt.figure(11) plt.clf() -gangof4(Hi*Po, Co, linspace(-2, 3)) +ct.gangof4(Hi*Po, Co, np.linspace(-2, 3)) if 'PYCONTROL_TEST_EXAMPLES' not in os.environ: plt.show() diff --git a/examples/pvtol.py b/examples/pvtol.py index 4f92f12fa..bc826a564 100644 --- a/examples/pvtol.py +++ b/examples/pvtol.py @@ -64,8 +64,6 @@ def _pvtol_flat_forward(states, inputs, params={}): F1, F2 = inputs # Use equations of motion for higher derivates - x1ddot = (F1 * cos(theta) - F2 * sin(theta)) / m - x2ddot = (F1 * sin(theta) + F2 * cos(theta) - m * g) / m thddot = (r * F1) / J # Flat output is a point above the vertical axis @@ -110,7 +108,6 @@ def _pvtol_flat_reverse(zflag, params={}): J = params.get('J', 0.0475) # inertia around pitch axis r = params.get('r', 0.25) # distance to center of force g = params.get('g', 9.8) # gravitational constant - c = params.get('c', 0.05) # damping factor (estimated) # Given the flat variables, solve for the state theta = np.arctan2(-zflag[0][2], zflag[1][2] + g) @@ -185,10 +182,6 @@ def _windy_update(t, x, u, params): def _noisy_update(t, x, u, params): # Get the inputs F1, F2, Dx, Dy = u[:4] - if u.shape[0] > 4: - Nx, Ny, Nth = u[4:] - else: - Nx, Ny, Nth = 0, 0, 0 # Get the system response from the original dynamics xdot, ydot, thetadot, xddot, yddot, thddot = \ @@ -196,7 +189,6 @@ def _noisy_update(t, x, u, params): # Get the parameter values we need m = params.get('m', 4.) # mass of aircraft - J = params.get('J', 0.0475) # inertia around pitch axis # Now add the disturbances xddot += Dx / m @@ -219,7 +211,6 @@ def _noisy_output(t, x, u, params): def pvtol_noisy_A(x, u, params={}): # Get the parameter values we need m = params.get('m', 4.) # mass of aircraft - J = params.get('J', 0.0475) # inertia around pitch axis c = params.get('c', 0.05) # damping factor (estimated) # Get the angle and compute sine and cosine diff --git a/examples/python-control_tutorial.ipynb b/examples/python-control_tutorial.ipynb index 4d718b050..6ac127758 100644 --- a/examples/python-control_tutorial.ipynb +++ b/examples/python-control_tutorial.ipynb @@ -94,7 +94,7 @@ "id": "qMVGK15gNQw2" }, "source": [ - "## Example 1: Open loop analysis of a coupled mass spring system\n", + "## Example 1: Open Loop Analysis of a Coupled Mass Spring System\n", "\n", "Consider the spring mass system below:\n", "\n", @@ -781,7 +781,7 @@ "id": "2f27f767-e012-45f9-8b76-cc040cfc89e2", "metadata": {}, "source": [ - "## Example 2: Trajectory tracking for a kinematic vehicle model\n", + "## Example 2: Trajectory Tracking for a Kinematic Vehicle Model\n", "\n", "This example illustrates the use of python-control to model, analyze, and design nonlinear control systems.\n", "\n", @@ -1213,7 +1213,7 @@ "id": "03b1fd75-579c-47da-805d-68f155957084", "metadata": {}, "source": [ - "## Computing environment" + "## Computing Environment" ] }, { diff --git a/examples/secord-matlab.py b/examples/secord-matlab.py index 6cef881c1..53fe69e6f 100644 --- a/examples/secord-matlab.py +++ b/examples/secord-matlab.py @@ -3,7 +3,8 @@ import os import matplotlib.pyplot as plt # MATLAB plotting functions -from control.matlab import * # MATLAB-like functions +import numpy as np +from control.matlab import ss, step, bode, nyquist, rlocus # MATLAB-like functions # Parameters defining the system m = 250.0 # system mass @@ -24,7 +25,7 @@ # Bode plot for the system plt.figure(2) -mag, phase, om = bode(sys, logspace(-2, 2), plot=True) +mag, phase, om = bode(sys, np.logspace(-2, 2), plot=True) plt.show(block=False) # Nyquist plot for the system @@ -32,7 +33,7 @@ nyquist(sys) plt.show(block=False) -# Root lcous plot for the system +# Root locus plot for the system rlocus(sys) if 'PYCONTROL_TEST_EXAMPLES' not in os.environ: diff --git a/examples/sisotool_example.py b/examples/sisotool_example.py index 6453bec74..44d7c0443 100644 --- a/examples/sisotool_example.py +++ b/examples/sisotool_example.py @@ -10,24 +10,24 @@ #%% import matplotlib.pyplot as plt -from control.matlab import * +import control as ct # first example, aircraft attitude equation -s = tf([1,0],[1]) +s = ct.tf([1,0],[1]) Kq = -24 T2 = 1.4 damping = 2/(13**.5) omega = 13**.5 H = (Kq*(1+T2*s))/(s*(s**2+2*damping*omega*s+omega**2)) plt.close('all') -sisotool(-H) +ct.sisotool(-H) #%% # a simple RL, with multiple poles in the origin plt.close('all') H = (s+0.3)/(s**4 + 4*s**3 + 6.25*s**2) -sisotool(H) +ct.sisotool(H) #%% @@ -43,4 +43,4 @@ plt.close('all') H = (b0 + b1*s + b2*s**2) / (a0 + a1*s + a2*s**2 + a3*s**3) -sisotool(H) +ct.sisotool(H) diff --git a/examples/slycot-import-test.py b/examples/slycot-import-test.py index 2df9b5b23..9c92fd2dc 100644 --- a/examples/slycot-import-test.py +++ b/examples/slycot-import-test.py @@ -5,7 +5,7 @@ """ import numpy as np -from control.matlab import * +import control as ct from control.exception import slycot_check # Parameters defining the system @@ -17,12 +17,12 @@ A = np.array([[1, -1, 1.], [1, -k/m, -b/m], [1, 1, 1]]) B = np.array([[0], [1/m], [1]]) C = np.array([[1., 0, 1.]]) -sys = ss(A, B, C, 0) +sys = ct.ss(A, B, C, 0) # Python control may be used without slycot, for example for a pole placement. # Eigenvalue placement w = [-3, -2, -1] -K = place(A, B, w) +K = ct.place(A, B, w) print("[python-control (from scipy)] K = ", K) print("[python-control (from scipy)] eigs = ", np.linalg.eig(A - B*K)[0]) diff --git a/examples/type2_type3.py b/examples/type2_type3.py index 52e0645e2..f0d79dc51 100644 --- a/examples/type2_type3.py +++ b/examples/type2_type3.py @@ -4,9 +4,10 @@ import os import matplotlib.pyplot as plt # Grab MATLAB plotting functions -from control.matlab import * # MATLAB-like functions -from numpy import pi -integrator = tf([0, 1], [1, 0]) # 1/s +import control as ct +import numpy as np + +integrator = ct.tf([0, 1], [1, 0]) # 1/s # Parameters defining the system J = 1.0 @@ -29,20 +30,20 @@ # System Transfer Functions # tricky because the disturbance (base motion) is coupled in by friction -closed_loop_type2 = feedback(C_type2*feedback(P, friction), gyro) +closed_loop_type2 = ct.feedback(C_type2*ct.feedback(P, friction), gyro) disturbance_rejection_type2 = P*friction/(1. + P*friction+P*C_type2) -closed_loop_type3 = feedback(C_type3*feedback(P, friction), gyro) +closed_loop_type3 = ct.feedback(C_type3*ct.feedback(P, friction), gyro) disturbance_rejection_type3 = P*friction/(1. + P*friction + P*C_type3) # Bode plot for the system plt.figure(1) -bode(closed_loop_type2, logspace(0, 2)*2*pi, dB=True, Hz=True) # blue -bode(closed_loop_type3, logspace(0, 2)*2*pi, dB=True, Hz=True) # green +ct.bode(closed_loop_type2, np.logspace(0, 2)*2*np.pi, dB=True, Hz=True) # blue +ct.bode(closed_loop_type3, np.logspace(0, 2)*2*np.pi, dB=True, Hz=True) # green plt.show(block=False) plt.figure(2) -bode(disturbance_rejection_type2, logspace(0, 2)*2*pi, Hz=True) # blue -bode(disturbance_rejection_type3, logspace(0, 2)*2*pi, Hz=True) # green +ct.bode(disturbance_rejection_type2, np.logspace(0, 2)*2*np.pi, Hz=True) # blue +ct.bode(disturbance_rejection_type3, np.logspace(0, 2)*2*np.pi, Hz=True) # green if 'PYCONTROL_TEST_EXAMPLES' not in os.environ: plt.show() diff --git a/examples/vehicle.py b/examples/vehicle.py index 07af35c9f..f89702d4e 100644 --- a/examples/vehicle.py +++ b/examples/vehicle.py @@ -3,7 +3,6 @@ import numpy as np import matplotlib.pyplot as plt -import control as ct import control.flatsys as fs # diff --git a/pyproject.toml b/pyproject.toml index d3754dc52..b47f7462c 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -10,13 +10,12 @@ build-backend = "setuptools.build_meta" name = "control" description = "Python Control Systems Library" authors = [{name = "Python Control Developers", email = "python-control-developers@lists.sourceforge.net"}] -license = {text = "BSD-3-Clause"} +license = "BSD-3-Clause" readme = "README.rst" classifiers = [ "Development Status :: 4 - Beta", "Intended Audience :: Science/Research", "Intended Audience :: Developers", - "License :: OSI Approved :: BSD License", "Programming Language :: Python :: 3", "Programming Language :: Python :: 3.10", "Programming Language :: Python :: 3.11", @@ -60,7 +59,7 @@ filterwarnings = [ [tool.ruff] # TODO: expand to cover all code -include = ['control/**.py'] +include = ['control/**.py', 'benchmarks/*.py', 'examples/*.py'] [tool.ruff.lint] select = [