diff --git a/.github/workflows/python-package-conda.yml b/.github/workflows/python-package-conda.yml index 10cf2d1a9..3f1910697 100644 --- a/.github/workflows/python-package-conda.yml +++ b/.github/workflows/python-package-conda.yml @@ -4,7 +4,7 @@ on: [push, pull_request] jobs: test-linux: - name: Python ${{ matrix.python-version }}${{ matrix.slycot && format(' with Slycot from {0}', matrix.slycot) || ' without Slycot' }}${{ matrix.array-and-matrix == 1 && ', array and matrix' || '' }} + name: Python ${{ matrix.python-version }}${{ matrix.slycot && format(' with Slycot from {0}', matrix.slycot) || ' without Slycot' }}${{ matrix.pandas && ', with pandas' || '' }}${{ matrix.array-and-matrix == 1 && ', array and matrix' || '' }} runs-on: ubuntu-latest strategy: @@ -12,10 +12,12 @@ jobs: matrix: python-version: [3.7, 3.9] slycot: ["", "conda"] + pandas: [""] array-and-matrix: [0] include: - python-version: 3.9 slycot: conda + pandas: conda array-and-matrix: 1 steps: @@ -41,6 +43,9 @@ jobs: if [[ '${{matrix.slycot}}' == 'conda' ]]; then conda install -c conda-forge slycot fi + if [[ '${{matrix.pandas}}' == 'conda' ]]; then + conda install -c conda-forge pandas + fi - name: Test with pytest env: diff --git a/control/__init__.py b/control/__init__.py index 386fa91c1..ad2685273 100644 --- a/control/__init__.py +++ b/control/__init__.py @@ -55,6 +55,7 @@ from .margins import * from .mateqn import * from .modelsimp import * +from .namedio import * from .nichols import * from .phaseplot import * from .pzmap import * diff --git a/control/canonical.py b/control/canonical.py index 7b2b58ef7..e714e5b8d 100644 --- a/control/canonical.py +++ b/control/canonical.py @@ -2,7 +2,7 @@ # RMM, 10 Nov 2012 from .exception import ControlNotImplemented, ControlSlycot -from .lti import issiso +from .namedio import issiso from .statesp import StateSpace, _convert_to_statespace from .statefbk import ctrb, obsv diff --git a/control/dtime.py b/control/dtime.py index c60778d00..b05d22b96 100644 --- a/control/dtime.py +++ b/control/dtime.py @@ -47,7 +47,7 @@ """ -from .lti import isctime +from .namedio import isctime from .statesp import StateSpace __all__ = ['sample_system', 'c2d'] diff --git a/control/exception.py b/control/exception.py index e28ba8609..f66eb7f30 100644 --- a/control/exception.py +++ b/control/exception.py @@ -71,3 +71,16 @@ def slycot_check(): except: slycot_installed = False return slycot_installed + + +# Utility function to see if pandas is installed +pandas_installed = None +def pandas_check(): + global pandas_installed + if pandas_installed is None: + try: + import pandas + pandas_installed = True + except: + pandas_installed = False + return pandas_installed diff --git a/control/flatsys/flatsys.py b/control/flatsys/flatsys.py index 2f20aa1e9..c01eb9127 100644 --- a/control/flatsys/flatsys.py +++ b/control/flatsys/flatsys.py @@ -156,6 +156,11 @@ def __init__(self, # Save the length of the flat flag + def __str__(self): + return f"{NonlinearIOSystem.__str__(self)}\n\n" \ + + f"Forward: {self.forward}\n" \ + + f"Reverse: {self.reverse}" + def forward(self, x, u, params={}): """Compute the flat flag given the states and input. diff --git a/control/flatsys/linflat.py b/control/flatsys/linflat.py index 94523cc0b..e4a31c6de 100644 --- a/control/flatsys/linflat.py +++ b/control/flatsys/linflat.py @@ -140,3 +140,13 @@ def reverse(self, zflag, params): x = self.Tinv @ z u = zflag[0][-1] - self.F @ z return np.reshape(x, self.nstates), np.reshape(u, self.ninputs) + + # Update function + def _rhs(self, t, x, u, params={}): + # Use LinearIOSystem._rhs instead of default (MRO) NonlinearIOSystem + return LinearIOSystem._rhs(self, t, x, u) + + # output function + def _out(self, t, x, u, params={}): + # Use LinearIOSystem._out instead of default (MRO) NonlinearIOSystem + return LinearIOSystem._out(self, t, x, u) diff --git a/control/flatsys/systraj.py b/control/flatsys/systraj.py index 5e390a7b5..9d425295b 100644 --- a/control/flatsys/systraj.py +++ b/control/flatsys/systraj.py @@ -37,6 +37,7 @@ # SUCH DAMAGE. import numpy as np +from ..timeresp import TimeResponseData class SystemTrajectory: """Class representing a system trajectory. @@ -117,3 +118,75 @@ def eval(self, tlist): self.system.reverse(zflag, self.params) return xd, ud + + # Return the system trajectory as a TimeResponseData object + def response(self, tlist, transpose=False, return_x=False, squeeze=None): + """Return the trajectory of a system as a TimeResponseData object + + Evaluate the trajectory at a list of time points, returning the state + and input vectors for the trajectory: + + response = traj.response(tlist) + time, yd, ud = response.time, response.outputs, response.inputs + + Parameters + ---------- + tlist : 1D array + List of times to evaluate the trajectory. + + transpose : bool, optional + If True, transpose all input and output arrays (for backward + compatibility with MATLAB and :func:`scipy.signal.lsim`). + Default value is False. + + return_x : bool, optional + If True, return the state vector when assigning to a tuple + (default = False). See :func:`forced_response` for more details. + + squeeze : bool, optional + By default, if a system is single-input, single-output (SISO) then + the output response is returned as a 1D array (indexed by time). + If squeeze=True, remove single-dimensional entries from the shape + of the output even if the system is not SISO. If squeeze=False, + keep the output as a 3D array (indexed by the output, input, and + time) even if the system is SISO. The default value can be set + using config.defaults['control.squeeze_time_response']. + + Returns + ------- + results : TimeResponseData + Time response represented as a :class:`TimeResponseData` object + containing the following properties: + + * time (array): Time values of the output. + + * outputs (array): Response of the system. If the system is SISO + and squeeze is not True, the array is 1D (indexed by time). If + the system is not SISO or ``squeeze`` is False, the array is 3D + (indexed by the output, trace, and time). + + * states (array): Time evolution of the state vector, represented + as either a 2D array indexed by state and time (if SISO) or a 3D + array indexed by state, trace, and time. Not affected by + ``squeeze``. + + * inputs (array): Input(s) to the system, indexed in the same + manner as ``outputs``. + + The return value of the system can also be accessed by assigning + the function to a tuple of length 2 (time, output) or of length 3 + (time, output, state) if ``return_x`` is ``True``. + + """ + # Compute the state and input response using the eval function + sys = self.system + xout, uout = self.eval(tlist) + yout = np.array([ + sys.output(tlist[i], xout[:, i], uout[:, i]) + for i in range(len(tlist))]).transpose() + + return TimeResponseData( + tlist, yout, xout, uout, issiso=sys.issiso(), + input_labels=sys.input_labels, output_labels=sys.output_labels, + state_labels=sys.state_labels, + transpose=transpose, return_x=return_x, squeeze=squeeze) diff --git a/control/frdata.py b/control/frdata.py index 19e865821..a33775afb 100644 --- a/control/frdata.py +++ b/control/frdata.py @@ -44,12 +44,17 @@ """ # External function declarations +from copy import copy from warnings import warn + import numpy as np from numpy import angle, array, empty, ones, \ real, imag, absolute, eye, linalg, where, sort from scipy.interpolate import splprep, splev + from .lti import LTI, _process_frequency_response +from .exception import pandas_check +from .namedio import NamedIOSystem, _process_namedio_keywords from . import config __all__ = ['FrequencyResponseData', 'FRD', 'frd'] @@ -112,7 +117,7 @@ class FrequencyResponseData(LTI): # Allow NDarray * StateSpace to give StateSpace._rmul_() priority # https://docs.scipy.org/doc/numpy/reference/arrays.classes.html - __array_priority__ = 11 # override ndarray and matrix types + __array_priority__ = 13 # override ndarray, StateSpace, I/O sys # # Class attributes @@ -139,23 +144,22 @@ def __init__(self, *args, **kwargs): The default constructor is FRD(d, w), where w is an iterable of frequency points, and d is the matching frequency data. - If d is a single list, 1d array, or tuple, a SISO system description + If d is a single list, 1D array, or tuple, a SISO system description is assumed. d can also be To call the copy constructor, call FRD(sys), where sys is a FRD object. To construct frequency response data for an existing LTI - object, other than an FRD, call FRD(sys, omega) + object, other than an FRD, call FRD(sys, omega). """ # TODO: discrete-time FRD systems? smooth = kwargs.pop('smooth', False) - # Make sure there were no extraneous keywords - if kwargs: - raise TypeError("unrecognized keywords: ", str(kwargs)) - + # + # Process positional arguments + # if len(args) == 2: if not isinstance(args[0], FRD) and isinstance(args[0], LTI): # not an FRD, but still a system, second argument should be @@ -172,13 +176,12 @@ def __init__(self, *args, **kwargs): else: # The user provided a response and a freq vector - self.fresp = array(args[0], dtype=complex) - if len(self.fresp.shape) == 1: - self.fresp = self.fresp.reshape(1, 1, len(args[0])) - self.omega = array(args[1], dtype=float) - if len(self.fresp.shape) != 3 or \ - self.fresp.shape[-1] != self.omega.shape[-1] or \ - len(self.omega.shape) != 1: + self.fresp = array(args[0], dtype=complex, ndmin=1) + if self.fresp.ndim == 1: + self.fresp = self.fresp.reshape(1, 1, -1) + self.omega = array(args[1], dtype=float, ndmin=1) + if self.fresp.ndim != 3 or self.omega.ndim != 1 or \ + self.fresp.shape[-1] != self.omega.shape[-1]: raise TypeError( "The frequency data constructor needs a 1-d or 3-d" " response data array and a matching frequency vector" @@ -196,6 +199,29 @@ def __init__(self, *args, **kwargs): raise ValueError( "Needs 1 or 2 arguments; received %i." % len(args)) + # + # Process key word arguments + # + # Keep track of return type + self.return_magphase=kwargs.pop('return_magphase', False) + if self.return_magphase not in (True, False): + raise ValueError("unknown return_magphase value") + + # Determine whether to squeeze the output + self.squeeze=kwargs.pop('squeeze', None) + if self.squeeze not in (None, True, False): + raise ValueError("unknown squeeze value") + + # Process namedio keywords + defaults = { + 'inputs': self.fresp.shape[1], 'outputs': self.fresp.shape[0]} + name, inputs, outputs, states, dt = _process_namedio_keywords( + kwargs, defaults, end=True) + + # Process signal names + NamedIOSystem.__init__( + self, name=name, inputs=inputs, outputs=outputs, dt=dt) + # create interpolation functions if smooth: self.ifunc = empty((self.fresp.shape[0], self.fresp.shape[1]), @@ -208,7 +234,29 @@ def __init__(self, *args, **kwargs): w=1.0/(absolute(self.fresp[i, j, :]) + 0.001), s=0.0) else: self.ifunc = None - super().__init__(self.fresp.shape[1], self.fresp.shape[0]) + + # + # Frequency response properties + # + # Different properties of the frequency response that can be used for + # analysis and characterization. + # + + @property + def magnitude(self): + return np.abs(self.fresp) + + @property + def phase(self): + return np.angle(self.fresp) + + @property + def frequency(self): + return self.omega + + @property + def response(self): + return self.fresp def __str__(self): """String representation of the transfer function.""" @@ -260,11 +308,13 @@ def __add__(self, other): # Check that the input-output sizes are consistent. if self.ninputs != other.ninputs: - raise ValueError("The first summand has %i input(s), but the \ -second has %i." % (self.ninputs, other.ninputs)) + raise ValueError( + "The first summand has %i input(s), but the " \ + "second has %i." % (self.ninputs, other.ninputs)) if self.noutputs != other.noutputs: - raise ValueError("The first summand has %i output(s), but the \ -second has %i." % (self.noutputs, other.noutputs)) + raise ValueError( + "The first summand has %i output(s), but the " \ + "second has %i." % (self.noutputs, other.noutputs)) return FRD(self.fresp + other.fresp, other.omega) @@ -460,7 +510,7 @@ def eval(self, omega, squeeze=None): return _process_frequency_response(self, omega, out, squeeze=squeeze) - def __call__(self, s, squeeze=None): + def __call__(self, s=None, squeeze=None, return_magphase=None): """Evaluate system's transfer function at complex frequencies. Returns the complex frequency response `sys(s)` of system `sys` with @@ -473,17 +523,31 @@ def __call__(self, s, squeeze=None): For a frequency response data object, the argument must be an imaginary number (since only the frequency response is defined). + If ``s`` is not given, this function creates a copy of a frequency + response data object with a different set of output settings. + Parameters ---------- s : complex scalar or 1D array_like - Complex frequencies - squeeze : bool, optional (default=True) + Complex frequencies. If not specified, return a copy of the + frequency response data object with updated settings for output + processing (``squeeze``, ``return_magphase``). + + squeeze : bool, optional If squeeze=True, remove single-dimensional entries from the shape of the output even if the system is not SISO. If squeeze=False, keep all indices (output, input and, if omega is array_like, frequency) even if the system is SISO. The default value can be set using config.defaults['control.squeeze_frequency_response']. + return_magphase : bool, optional + If True, then a frequency response data object will enumerate as a + tuple of the form (mag, phase, omega) where where ``mag`` is the + magnitude (absolute value, not dB or log10) of the system + frequency response, ``phase`` is the wrapped phase in radians of + the system frequency response, and ``omega`` is the (sorted) + frequencies at which the response was evaluated. + Returns ------- fresp : complex ndarray @@ -502,6 +566,17 @@ def __call__(self, s, squeeze=None): frequency values. """ + if s is None: + # Create a copy of the response with new keywords + response = copy(self) + + # Update any keywords that we were passed + response.squeeze = self.squeeze if squeeze is None else squeeze + response.return_magphase = self.return_magphase \ + if return_magphase is None else return_magphase + + return response + # Make sure that we are operating on a simple list if len(np.atleast_1d(s).shape) > 1: raise ValueError("input list must be 1D") @@ -516,6 +591,22 @@ def __call__(self, s, squeeze=None): else: return self.eval(complex(s).imag, squeeze=squeeze) + # Implement iter to allow assigning to a tuple + def __iter__(self): + fresp = _process_frequency_response( + self, self.omega, self.fresp, squeeze=self.squeeze) + if not self.return_magphase: + return iter((self.omega, fresp)) + return iter((np.abs(fresp), np.angle(fresp), self.omega)) + + # Implement (thin) getitem to allow access via legacy indexing + def __getitem__(self, index): + return list(self.__iter__())[index] + + # Implement (thin) len to emulate legacy testing interface + def __len__(self): + return 3 if self.return_magphase else 2 + def freqresp(self, omega): """(deprecated) Evaluate transfer function at complex frequencies. @@ -551,6 +642,22 @@ def feedback(self, other=1, sign=-1): return FRD(fresp, other.omega, smooth=(self.ifunc is not None)) + # Convert to pandas + def to_pandas(self): + if not pandas_check(): + ImportError('pandas not installed') + import pandas + + # Create a dict for setting up the data frame + data = {'omega': self.omega} + data.update( + {'H_{%s, %s}' % (out, inp): self.fresp[i, j] \ + for i, out in enumerate(self.output_labels) \ + for j, inp in enumerate(self.input_labels)}) + + return pandas.DataFrame(data) + + # # Allow FRD as an alias for the FrequencyResponseData class # @@ -561,8 +668,6 @@ def feedback(self, other=1, sign=-1): # FrequenceResponseData and then assigning FRD to point to the same object # fixes this problem. # - - FRD = FrequencyResponseData diff --git a/control/freqplot.py b/control/freqplot.py index 7a1243c6c..7f29dce36 100644 --- a/control/freqplot.py +++ b/control/freqplot.py @@ -204,8 +204,9 @@ def bode_plot(syslist, omega=None, initial_phase = config._get_param( 'freqplot', 'initial_phase', kwargs, None, pop=True) omega_num = config._get_param('freqplot', 'number_of_samples', omega_num) + # If argument was a singleton, turn it into a tuple - if not hasattr(syslist, '__iter__'): + if not isinstance(syslist, (list, tuple)): syslist = (syslist,) omega, omega_range_given = _determine_omega_vector( @@ -678,8 +679,8 @@ def nyquist_plot(syslist, omega=None, plot=True, omega_limits=None, indent_direction = config._get_param( 'nyquist', 'indent_direction', kwargs, _nyquist_defaults, pop=True) - # If argument was a singleton, turn it into a list - if not hasattr(syslist, '__iter__'): + # If argument was a singleton, turn it into a tuple + if not isinstance(syslist, (list, tuple)): syslist = (syslist,) omega, omega_range_given = _determine_omega_vector( @@ -721,11 +722,11 @@ def nyquist_plot(syslist, omega=None, plot=True, omega_limits=None, if isinstance(sys, (StateSpace, TransferFunction)) \ and indent_direction != 'none': if sys.isctime(): - splane_poles = sys.pole() + splane_poles = sys.poles() else: # map z-plane poles to s-plane, ignoring any at the origin # because we don't need to indent for them - zplane_poles = sys.pole() + zplane_poles = sys.poles() zplane_poles = zplane_poles[~np.isclose(abs(zplane_poles), 0.)] splane_poles = np.log(zplane_poles)/sys.dt @@ -1109,7 +1110,7 @@ def singular_values_plot(syslist, omega=None, omega_num = config._get_param('freqplot', 'number_of_samples', omega_num) # If argument was a singleton, turn it into a tuple - if not hasattr(syslist, '__iter__'): + if not isinstance(syslist, (list, tuple)): syslist = (syslist,) omega, omega_range_given = _determine_omega_vector( @@ -1327,8 +1328,8 @@ def _default_frequency_range(syslist, Hz=None, number_of_samples=None, try: # Add new features to the list if sys.isctime(): - features_ = np.concatenate((np.abs(sys.pole()), - np.abs(sys.zero()))) + features_ = np.concatenate( + (np.abs(sys.poles()), np.abs(sys.zeros()))) # Get rid of poles and zeros at the origin toreplace = np.isclose(features_, 0.0) if np.any(toreplace): @@ -1338,8 +1339,7 @@ def _default_frequency_range(syslist, Hz=None, number_of_samples=None, # TODO: What distance to the Nyquist frequency is appropriate? freq_interesting.append(fn * 0.9) - features_ = np.concatenate((sys.pole(), - sys.zero())) + features_ = np.concatenate((sys.poles(), sys.zeros())) # Get rid of poles and zeros on the real axis (imag==0) # * origin and real < 0 # * at 1.: would result in omega=0. (logaritmic plot!) diff --git a/control/iosys.py b/control/iosys.py index 161f06510..e3719614b 100644 --- a/control/iosys.py +++ b/control/iosys.py @@ -31,13 +31,14 @@ import copy from warnings import warn -from .namedio import _NamedIOStateSystem, _process_signal_list +from .lti import LTI +from .namedio import NamedIOSystem, _process_signal_list, \ + _process_namedio_keywords, isctime, isdtime, common_timebase from .statesp import StateSpace, tf2ss, _convert_to_statespace -from .statesp import _ss, _rss_generate +from .statesp import _rss_generate from .xferfcn import TransferFunction from .timeresp import _check_convert_array, _process_time_response, \ TimeResponseData -from .lti import isctime, isdtime, common_timebase from . import config __all__ = ['InputOutputSystem', 'LinearIOSystem', 'NonlinearIOSystem', @@ -55,7 +56,7 @@ } -class InputOutputSystem(_NamedIOStateSystem): +class InputOutputSystem(NamedIOSystem): """A class for representing input/output systems. The InputOutputSystem class allows (possibly nonlinear) input/output @@ -68,7 +69,7 @@ class for a set of subclasses that are used to implement specific ---------- inputs : int, list of str, or None Description of the system inputs. This can be given as an integer - count or as a list of strings that name the individual signals. If an + count or a list of strings that name the individual signals. If an integer count is specified, the names of the signal will be of the form `s[i]` (where `s` is one of `u`, `y`, or `x`). If this parameter is not given or given as `None`, the relevant quantity will be @@ -79,17 +80,16 @@ class for a set of subclasses that are used to implement specific states : int, list of str, or None Description of the system states. Same format as `inputs`. dt : None, True or float, optional - System timebase. 0 (default) indicates continuous - time, True indicates discrete time with unspecified sampling - time, positive number is discrete time with specified - sampling time, None indicates unspecified timebase (either - continuous or discrete time). - params : dict, optional - Parameter values for the systems. Passed to the evaluation functions - for the system as default values, overriding internal defaults. + System timebase. 0 (default) indicates continuous time, True + indicates discrete time with unspecified sampling time, positive + number is discrete time with specified sampling time, None indicates + unspecified timebase (either continuous or discrete time). name : string, optional System name (used for specifying signals). If unspecified, a generic name is generated with a unique integer id. + params : dict, optional + Parameter values for the systems. Passed to the evaluation functions + for the system as default values, overriding internal defaults. Attributes ---------- @@ -126,8 +126,7 @@ class for a set of subclasses that are used to implement specific # Allow ndarray * InputOutputSystem to give IOSystem._rmul_() priority __array_priority__ = 12 # override ndarray, matrix, SS types - def __init__(self, inputs=None, outputs=None, states=None, params={}, - name=None, **kwargs): + def __init__(self, params={}, **kwargs): """Create an input/output system. The InputOutputSystem constructor is used to create an input/output @@ -139,19 +138,18 @@ def __init__(self, inputs=None, outputs=None, states=None, params={}, """ # Store the system name, inputs, outputs, and states - _NamedIOStateSystem.__init__( - self, inputs=inputs, outputs=outputs, states=states, name=name) + name, inputs, outputs, states, dt = _process_namedio_keywords( + kwargs, end=True) + + # Initialize the data structure + # Note: don't use super() to override LinearIOSystem/StateSpace MRO + NamedIOSystem.__init__( + self, inputs=inputs, outputs=outputs, + states=states, name=name, dt=dt) # default parameters self.params = params.copy() - # timebase - self.dt = kwargs.pop('dt', config.defaults['control.default_dt']) - - # Make sure there were no extraneous keywords - if kwargs: - raise TypeError("unrecognized keywords: ", str(kwargs)) - def __mul__(sys2, sys1): """Multiply two input/output systems (series interconnection)""" # Note: order of arguments is flipped so that self = sys2, @@ -165,7 +163,8 @@ def __mul__(sys2, sys1): elif isinstance(sys1, np.ndarray): sys1 = LinearIOSystem(StateSpace([], [], [], sys1)) - elif isinstance(sys1, (StateSpace, TransferFunction)): + elif isinstance(sys1, (StateSpace, TransferFunction)) and \ + not isinstance(sys1, LinearIOSystem): sys1 = LinearIOSystem(sys1) elif not isinstance(sys1, InputOutputSystem): @@ -211,7 +210,8 @@ def __rmul__(sys1, sys2): elif isinstance(sys2, np.ndarray): sys2 = LinearIOSystem(StateSpace([], [], [], sys2)) - elif isinstance(sys2, (StateSpace, TransferFunction)): + elif isinstance(sys2, (StateSpace, TransferFunction)) and \ + not isinstance(sys2, LinearIOSystem): sys2 = LinearIOSystem(sys2) elif not isinstance(sys2, InputOutputSystem): @@ -229,7 +229,8 @@ def __add__(sys1, sys2): elif isinstance(sys2, np.ndarray): sys2 = LinearIOSystem(StateSpace([], [], [], sys2)) - elif isinstance(sys2, (StateSpace, TransferFunction)): + elif isinstance(sys2, (StateSpace, TransferFunction)) and \ + not isinstance(sys2, LinearIOSystem): sys2 = LinearIOSystem(sys2) elif not isinstance(sys2, InputOutputSystem): @@ -266,7 +267,8 @@ def __radd__(sys1, sys2): elif isinstance(sys2, np.ndarray): sys2 = LinearIOSystem(StateSpace([], [], [], sys2)) - elif isinstance(sys2, (StateSpace, TransferFunction)): + elif isinstance(sys2, (StateSpace, TransferFunction)) and \ + not isinstance(sys2, LinearIOSystem): sys2 = LinearIOSystem(sys2) elif not isinstance(sys2, InputOutputSystem): @@ -284,7 +286,8 @@ def __sub__(sys1, sys2): elif isinstance(sys2, np.ndarray): sys2 = LinearIOSystem(StateSpace([], [], [], sys2)) - elif isinstance(sys2, (StateSpace, TransferFunction)): + elif isinstance(sys2, (StateSpace, TransferFunction)) and \ + not isinstance(sys2, LinearIOSystem): sys2 = LinearIOSystem(sys2) elif not isinstance(sys2, InputOutputSystem): @@ -321,7 +324,8 @@ def __rsub__(sys1, sys2): elif isinstance(sys2, np.ndarray): sys2 = LinearIOSystem(StateSpace([], [], [], sys2)) - elif isinstance(sys2, (StateSpace, TransferFunction)): + elif isinstance(sys2, (StateSpace, TransferFunction)) and \ + not isinstance(sys2, LinearIOSystem): sys2 = LinearIOSystem(sys2) elif not isinstance(sys2, InputOutputSystem): @@ -579,15 +583,6 @@ def linearize(self, x0, u0, t=0, params={}, eps=1e-6, return linsys - def copy(self, newname=None): - """Make a copy of an input/output system.""" - dup_prefix = config.defaults['iosys.duplicate_system_name_prefix'] - dup_suffix = config.defaults['iosys.duplicate_system_name_suffix'] - newsys = copy.copy(self) - newsys.name = self._name_or_default( - dup_prefix + self.name + dup_suffix if not newname else newname) - return newsys - class LinearIOSystem(InputOutputSystem, StateSpace): """Input/output representation of a linear (state space) system. @@ -616,12 +611,12 @@ class LinearIOSystem(InputOutputSystem, StateSpace): discrete time with unspecified sampling time, positive number is discrete time with specified sampling time, None indicates unspecified timebase (either continuous or discrete time). - params : dict, optional - Parameter values for the systems. Passed to the evaluation functions - for the system as default values, overriding internal defaults. name : string, optional System name (used for specifying signals). If unspecified, a generic name is generated with a unique integer id. + params : dict, optional + Parameter values for the systems. Passed to the evaluation functions + for the system as default values, overriding internal defaults. Attributes ---------- @@ -632,8 +627,7 @@ class LinearIOSystem(InputOutputSystem, StateSpace): See :class:`~control.StateSpace` for inherited attributes. """ - def __init__(self, linsys, inputs=None, outputs=None, states=None, - name=None, **kwargs): + def __init__(self, linsys, **kwargs): """Create an I/O system from a state space linear system. Converts a :class:`~control.StateSpace` system into an @@ -649,33 +643,19 @@ def __init__(self, linsys, inputs=None, outputs=None, states=None, raise TypeError("Linear I/O system must be a state space " "or transfer function object") - # Look for 'input' and 'output' parameter name variants - states = _parse_signal_parameter(states, 'state', kwargs) - inputs = _parse_signal_parameter(inputs, 'input', kwargs) - outputs = _parse_signal_parameter(outputs, 'output', kwargs, end=True) + # Process keyword arguments + name, inputs, outputs, states, dt = _process_namedio_keywords( + kwargs, linsys, end=True) # Create the I/O system object - super(LinearIOSystem, self).__init__( - inputs=linsys.ninputs, outputs=linsys.noutputs, - states=linsys.nstates, params={}, dt=linsys.dt, name=name) + # Note: don't use super() to override StateSpace MRO + InputOutputSystem.__init__( + self, inputs=inputs, outputs=outputs, states=states, + params={}, dt=dt, name=name) # Initalize additional state space variables - StateSpace.__init__(self, linsys, remove_useless_states=False) - - # Process input, output, state lists, if given - # Make sure they match the size of the linear system - ninputs, self.input_index = _process_signal_list( - inputs if inputs is not None else linsys.ninputs, prefix='u') - if ninputs is not None and linsys.ninputs != ninputs: - raise ValueError("Wrong number/type of inputs given.") - noutputs, self.output_index = _process_signal_list( - outputs if outputs is not None else linsys.noutputs, prefix='y') - if noutputs is not None and linsys.noutputs != noutputs: - raise ValueError("Wrong number/type of outputs given.") - nstates, self.state_index = _process_signal_list( - states if states is not None else linsys.nstates, prefix='x') - if nstates is not None and linsys.nstates != nstates: - raise ValueError("Wrong number/type of states given.") + StateSpace.__init__( + self, linsys, remove_useless_states=False, init_namedio=False) # The following text needs to be replicated from StateSpace in order for # this entry to show up properly in sphinx doccumentation (not sure why, @@ -705,6 +685,10 @@ def _out(self, t, x, u): + self.D @ np.reshape(u, (-1, 1)) return np.array(y).reshape((-1,)) + def __repr__(self): + # Need to define so that I/O system gets used instead of StateSpace + return InputOutputSystem.__repr__(self) + def __str__(self): return InputOutputSystem.__str__(self) + "\n\n" \ + StateSpace.__str__(self) @@ -752,11 +736,6 @@ class NonlinearIOSystem(InputOutputSystem): states : int, list of str, or None, optional Description of the system states. Same format as `inputs`. - params : dict, optional - Parameter values for the systems. Passed to the evaluation - functions for the system as default values, overriding internal - defaults. - dt : timebase, optional The timebase for the system, used to specify whether the system is operating in continuous or discrete time. It can have the @@ -771,28 +750,27 @@ class NonlinearIOSystem(InputOutputSystem): System name (used for specifying signals). If unspecified, a generic name is generated with a unique integer id. + params : dict, optional + Parameter values for the systems. Passed to the evaluation + functions for the system as default values, overriding internal + defaults. + """ - def __init__(self, updfcn, outfcn=None, inputs=None, outputs=None, - states=None, params={}, name=None, **kwargs): + def __init__(self, updfcn, outfcn=None, params={}, **kwargs): """Create a nonlinear I/O system given update and output functions.""" - # Look for 'input' and 'output' parameter name variants - inputs = _parse_signal_parameter(inputs, 'input', kwargs) - outputs = _parse_signal_parameter(outputs, 'output', kwargs) - - # Store the update and output functions - self.updfcn = updfcn - self.outfcn = outfcn + # Process keyword arguments + name, inputs, outputs, states, dt = _process_namedio_keywords( + kwargs, end=True) # Initialize the rest of the structure - dt = kwargs.pop('dt', config.defaults['control.default_dt']) - super(NonlinearIOSystem, self).__init__( + super().__init__( inputs=inputs, outputs=outputs, states=states, params=params, dt=dt, name=name ) - # Make sure there were no extraneous keywords - if kwargs: - raise TypeError("unrecognized keywords: ", str(kwargs)) + # Store the update and output functions + self.updfcn = updfcn + self.outfcn = outfcn # Check to make sure arguments are consistent if updfcn is None: @@ -815,6 +793,11 @@ def __init__(self, updfcn, outfcn=None, inputs=None, outputs=None, # Initialize current parameters to default parameters self._current_params = params.copy() + def __str__(self): + return f"{InputOutputSystem.__str__(self)}\n\n" + \ + f"Update: {self.updfcn}\n" + \ + f"Output: {self.outfcn}" + # Return the value of a static nonlinear system def __call__(sys, u, params=None, squeeze=None): """Evaluate a (static) nonlinearity at a given input value @@ -880,31 +863,33 @@ class InterconnectedSystem(InputOutputSystem): """ def __init__(self, syslist, connections=[], inplist=[], outlist=[], - inputs=None, outputs=None, states=None, - params={}, dt=None, name=None, **kwargs): + params={}, warn_duplicate=None, **kwargs): """Create an I/O system from a list of systems + connection info.""" - - # Look for 'input' and 'output' parameter name variants - inputs = _parse_signal_parameter(inputs, 'input', kwargs) - outputs = _parse_signal_parameter(outputs, 'output', kwargs, end=True) - # Convert input and output names to lists if they aren't already if not isinstance(inplist, (list, tuple)): inplist = [inplist] if not isinstance(outlist, (list, tuple)): outlist = [outlist] - # Check to make sure all systems are consistent + # Process keyword arguments + defaults = {'inputs': len(inplist), 'outputs': len(outlist)} + name, inputs, outputs, states, dt = _process_namedio_keywords( + kwargs, defaults, end=True) + + # Initialize the system list and index self.syslist = syslist self.syslist_index = {} - nstates = 0 - self.state_offset = [] - ninputs = 0 - self.input_offset = [] - noutputs = 0 - self.output_offset = [] + + # Initialize the input, output, and state counts, indices + nstates, self.state_offset = 0, [] + ninputs, self.input_offset = 0, [] + noutputs, self.output_offset = 0, [] + + # Keep track of system objects and names we have already seen sysobj_name_dct = {} sysname_count_dct = {} + + # Go through the system list and keep track of counts, offsets for sysidx, sys in enumerate(syslist): # Make sure time bases are consistent dt = common_timebase(dt, sys.dt) @@ -929,17 +914,32 @@ def __init__(self, syslist, connections=[], inplist=[], outlist=[], # Check for duplicate systems or duplicate names # Duplicates are renamed sysname_1, sysname_2, etc. if sys in sysobj_name_dct: - sys = sys.copy() - warn("Duplicate object found in system list: %s. " - "Making a copy" % str(sys.name)) + # Make a copy of the object using a new name + if warn_duplicate is None and sys._generic_name_check(): + # Make a copy w/out warning, using generic format + sys = sys.copy(use_prefix_suffix=False) + warn_flag = False + else: + sys = sys.copy() + warn_flag = warn_duplicate + + # Warn the user about the new object + if warn_flag is not False: + warn("duplicate object found in system list; " + "created copy: %s" % str(sys.name), stacklevel=2) + + # Check to see if the system name shows up more than once if sys.name is not None and sys.name in sysname_count_dct: count = sysname_count_dct[sys.name] sysname_count_dct[sys.name] += 1 sysname = sys.name + "_" + str(count) sysobj_name_dct[sys] = sysname self.syslist_index[sysname] = sysidx - warn("Duplicate name found in system list. " - "Renamed to {}".format(sysname)) + + if warn_duplicate is not False: + warn("duplicate name found in system list; " + "renamed to {}".format(sysname), stacklevel=2) + else: sysname_count_dct[sys.name] = 1 sysobj_name_dct[sys] = sys.name @@ -952,23 +952,18 @@ def __init__(self, syslist, connections=[], inplist=[], outlist=[], states += [sysname + state_name_delim + statename for statename in sys.state_index.keys()] + # Make sure we the state list is the right length (internal check) + if isinstance(states, list) and len(states) != nstates: + raise RuntimeError( + f"construction of state labels failed; found: " + f"{len(states)} labels; expecting {nstates}") + # Create the I/O system - super(InterconnectedSystem, self).__init__( - inputs=len(inplist), outputs=len(outlist), + # Note: don't use super() to override LinearICSystem/StateSpace MRO + InputOutputSystem.__init__( + self, inputs=inputs, outputs=outputs, states=states, params=params, dt=dt, name=name) - # If input or output list was specified, update it - if inputs is not None: - nsignals, self.input_index = \ - _process_signal_list(inputs, prefix='u') - if nsignals is not None and len(inplist) != nsignals: - raise ValueError("Wrong number/type of inputs given.") - if outputs is not None: - nsignals, self.output_index = \ - _process_signal_list(outputs, prefix='y') - if nsignals is not None and len(outlist) != nsignals: - raise ValueError("Wrong number/type of outputs given.") - # Convert the list of interconnections to a connection map (matrix) self.connect_map = np.zeros((ninputs, noutputs)) for connection in connections: @@ -1505,8 +1500,8 @@ class LinearICSystem(InterconnectedSystem, LinearIOSystem): :class:`StateSpace` class structure, allowing it to be passed to functions that expect a :class:`StateSpace` system. - This class is usually generated using :func:`~control.interconnect` and - not called directly + This class is generated using :func:`~control.interconnect` and + not called directly. """ @@ -1514,18 +1509,15 @@ def __init__(self, io_sys, ss_sys=None): if not isinstance(io_sys, InterconnectedSystem): raise TypeError("First argument must be an interconnected system.") - # Create the I/O system object + # Create the (essentially empty) I/O system object InputOutputSystem.__init__( self, name=io_sys.name, params=io_sys.params) - # Copy over the I/O systems attributes + # Copy over the named I/O system attributes self.syslist = io_sys.syslist - self.ninputs = io_sys.ninputs - self.noutputs = io_sys.noutputs - self.nstates = io_sys.nstates - self.input_index = io_sys.input_index - self.output_index = io_sys.output_index - self.state_index = io_sys.state_index + self.ninputs, self.input_index = io_sys.ninputs, io_sys.input_index + self.noutputs, self.output_index = io_sys.noutputs, io_sys.output_index + self.nstates, self.state_index = io_sys.nstates, io_sys.state_index self.dt = io_sys.dt # Copy over the attributes from the interconnected system @@ -1545,13 +1537,14 @@ def __init__(self, io_sys, ss_sys=None): # Initialize the state space attributes if isinstance(ss_sys, StateSpace): - # Make sure the dimension match + # Make sure the dimensions match if io_sys.ninputs != ss_sys.ninputs or \ io_sys.noutputs != ss_sys.noutputs or \ io_sys.nstates != ss_sys.nstates: raise ValueError("System dimensions for first and second " "arguments must match.") - StateSpace.__init__(self, ss_sys, remove_useless_states=False) + StateSpace.__init__( + self, ss_sys, remove_useless_states=False, init_namedio=False) else: raise TypeError("Second argument must be a state space system.") @@ -1662,15 +1655,14 @@ def input_output_response( # Process keyword arguments # - # Allow method as an alternative to solve_ivp_method - if kwargs.get('method', None): - solve_ivp_kwargs['method'] = kwargs.pop('method') - # Figure out the method to be used if kwargs.get('solve_ivp_method', None): if kwargs.get('method', None): raise ValueError("ivp_method specified more than once") solve_ivp_kwargs['method'] = kwargs.pop('solve_ivp_method') + elif kwargs.get('method', None): + # Allow method as an alternative to solve_ivp_method + solve_ivp_kwargs['method'] = kwargs.pop('method') # Set the default method to 'RK45' if solve_ivp_kwargs.get('method', None) is None: @@ -1736,6 +1728,23 @@ def input_output_response( U = U.reshape(-1, ntimepts) ninputs = U.shape[0] + # If we were passed a list of initial states, concatenate them + if isinstance(X0, (tuple, list)): + X0_list = [] + for i, x0 in enumerate(X0): + x0 = np.array(x0).reshape(-1) # convert everyting to 1D array + X0_list += x0.tolist() # add elements to initial state + + # Save the newly created input vector + X0 = np.array(X0_list) + + # If the initial state is too short, make it longer (NB: sys.nstates + # could be None if nstates comes from size of initial condition) + if sys.nstates and isinstance(X0, np.ndarray) and X0.size < sys.nstates: + if X0[-1] != 0: + warn("initial state too short; padding with zeros") + X0 = np.hstack([X0, np.zeros(sys.nstates - X0.size)]) + # If we were passed a list of initial states, concatenate them if isinstance(X0, (tuple, list)): X0_list = [] @@ -2193,24 +2202,23 @@ def linearize(sys, xeq, ueq=[], t=0, params={}, **kw): The linearization of the system, as a :class:`~control.LinearIOSystem` object (which is also a :class:`~control.StateSpace` object. + Additional Parameters + --------------------- + inputs : int, list of str or None, optional + Description of the system inputs. If not specified, the origional + system inputs are used. See :class:`InputOutputSystem` for more + information. + outputs : int, list of str or None, optional + Description of the system outputs. Same format as `inputs`. + states : int, list of str, or None, optional + Description of the system states. Same format as `inputs`. + """ if not isinstance(sys, InputOutputSystem): raise TypeError("Can only linearize InputOutputSystem types") return sys.linearize(xeq, ueq, t=t, params=params, **kw) -# Utility function to parse a signal parameter -def _parse_signal_parameter(value, name, kwargs, end=False): - # Check kwargs for a variant of the parameter name - if value is None and name in kwargs: - value = kwargs.pop(name) - - if end and kwargs: - raise TypeError("unrecognized keywords: ", str(kwargs)) - - return value - - def _find_size(sysval, vecval): """Utility function to find the size of a system parameter @@ -2236,12 +2244,17 @@ def ss(*args, **kwargs): Create a state space system. - The function accepts either 1, 4 or 5 parameters: + The function accepts either 1, 2, 4 or 5 parameters: ``ss(sys)`` Convert a linear system into space system form. Always creates a new system, even if sys is already a state space system. + ``ss(updfcn, outfucn)``` + Create a nonlinear input/output system with update function ``updfcn`` + and output function ``outfcn``. See :class:`NonlinearIOSystem` for + more information. + ``ss(A, B, C, D)`` Create a state space system from the matrices of its state and output equations: @@ -2264,6 +2277,10 @@ def ss(*args, **kwargs): Everything that the constructor of :class:`numpy.matrix` accepts is permissible here too. + ``ss(args, inputs=['u1', ..., 'up'], outputs=['y1', ..., 'yq'], + states=['x1', ..., 'xn']) + Create a system with named input, output, and state signals. + Parameters ---------- sys : StateSpace or TransferFunction @@ -2279,7 +2296,8 @@ def ss(*args, **kwargs): inputs, outputs, states : str, or list of str, optional List of strings that name the individual signals. If this parameter is not given or given as `None`, the signal names will be of the - form `s[i]` (where `s` is one of `u`, `y`, or `x`). + form `s[i]` (where `s` is one of `u`, `y`, or `x`). See + :class:`InputOutputSystem` for more information. name : string, optional System name (used for specifying signals). If unspecified, a generic name is generated with a unique integer id. @@ -2310,32 +2328,61 @@ def ss(*args, **kwargs): >>> sys2 = ss(sys_tf) """ - # Extract the keyword arguments needed for StateSpace (via _ss) - ss_kwlist = ('dt', 'remove_useless_states') - ss_kwargs = {} - for kw in ss_kwlist: - if kw in kwargs: - ss_kwargs[kw] = kwargs.pop(kw) + # See if this is a nonlinear I/O system + if len(args) > 0 and (hasattr(args[0], '__call__') or args[0] is None) \ + and not isinstance(args[0], (InputOutputSystem, LTI)): + # Function as first (or second) argument => assume nonlinear IO system + return NonlinearIOSystem(*args, **kwargs) + + elif len(args) == 4 or len(args) == 5: + # Create a state space function from A, B, C, D[, dt] + sys = LinearIOSystem(StateSpace(*args, **kwargs)) + + elif len(args) == 1: + sys = args[0] + if isinstance(sys, LTI): + # Check for system with no states and specified state names + if sys.nstates is None and 'states' in kwargs: + warn("state labels specified for " + "non-unique state space realization") + + # Create a state space system from an LTI system + sys = LinearIOSystem(_convert_to_statespace(sys), **kwargs) + else: + raise TypeError("ss(sys): sys must be a StateSpace or " + "TransferFunction object. It is %s." % type(sys)) + else: + raise TypeError( + "Needs 1, 4, or 5 arguments; received %i." % len(args)) - # Create the statespace system and then convert to I/O system - sys = _ss(*args, keywords=ss_kwargs) - return LinearIOSystem(sys, **kwargs) + return sys def rss(states=1, outputs=1, inputs=1, strictly_proper=False, **kwargs): - """ - Create a stable *continuous* random state space object. + """Create a stable random state space object. Parameters ---------- - states : int - Number of state variables - outputs : int - Number of system outputs - inputs : int - Number of system inputs + inputs : int, list of str, or None + Description of the system inputs. This can be given as an integer + count or as a list of strings that name the individual signals. If an + integer count is specified, the names of the signal will be of the + form `s[i]` (where `s` is one of `u`, `y`, or `x`). + outputs : int, list of str, or None + Description of the system outputs. Same format as `inputs`. + states : int, list of str, or None + Description of the system states. Same format as `inputs`. strictly_proper : bool, optional If set to 'True', returns a proper system (no direct term). + dt : None, True or float, optional + System timebase. 0 (default) indicates continuous + time, True indicates discrete time with unspecified sampling + time, positive number is discrete time with specified + sampling time, None indicates unspecified timebase (either + continuous or discrete time). + name : string, optional + System name (used for specifying signals). If unspecified, a generic + name is generated with a unique integer id. Returns ------- @@ -2347,85 +2394,136 @@ def rss(states=1, outputs=1, inputs=1, strictly_proper=False, **kwargs): ValueError if any input is not a positive integer - See Also - -------- - drss - Notes ----- If the number of states, inputs, or outputs is not specified, then the - missing numbers are assumed to be 1. The poles of the returned system - will always have a negative real part. + missing numbers are assumed to be 1. If dt is not specified or is given + as 0 or None, the poles of the returned system will always have a + negative real part. If dt is True or a postive float, the poles of the + returned system will have magnitude less than 1. """ - # Process states, inputs, outputs (ignoring names) + # Process keyword arguments + kwargs.update({'states': states, 'outputs': outputs, 'inputs': inputs}) + name, inputs, outputs, states, dt = _process_namedio_keywords( + kwargs, end=True) + + # Figure out the size of the sytem nstates, _ = _process_signal_list(states) ninputs, _ = _process_signal_list(inputs) noutputs, _ = _process_signal_list(outputs) sys = _rss_generate( - nstates, ninputs, noutputs, 'c', strictly_proper=strictly_proper) + nstates, ninputs, noutputs, 'c' if not dt else 'd', name=name, + strictly_proper=strictly_proper) + return LinearIOSystem( - sys, states=states, inputs=inputs, outputs=outputs, **kwargs) + sys, name=name, states=states, inputs=inputs, outputs=outputs, dt=dt) + + +def drss(*args, **kwargs): + """Create a stable, discrete-time, random state space system + Create a stable *discrete time* random state space object. This + function calls :func:`rss` using either the `dt` keyword provided by + the user or `dt=True` if not specified. -def drss(states=1, outputs=1, inputs=1, strictly_proper=False, **kwargs): """ - Create a stable *discrete* random state space object. + # Make sure the timebase makes sense + if 'dt' in kwargs: + dt = kwargs['dt'] + + if dt == 0: + raise ValueError("drss called with continuous timebase") + elif dt is None: + warn("drss called with unspecified timebase; " + "system may be interpreted as continuous time") + kwargs['dt'] = True # force rss to generate discrete time sys + else: + dt = True + kwargs['dt'] = True + + # Create the system + sys = rss(*args, **kwargs) + + # Reset the timebase (in case it was specified as None) + sys.dt = dt + + return sys + + +# Convert a state space system into an input/output system (wrapper) +def ss2io(*args, **kwargs): + return LinearIOSystem(*args, **kwargs) +ss2io.__doc__ = LinearIOSystem.__init__.__doc__ + + +# Convert a transfer function into an input/output system (wrapper) +def tf2io(*args, **kwargs): + """tf2io(sys) + + Convert a transfer function into an I/O system + + The function accepts either 1 or 2 parameters: + + ``tf2io(sys)`` + Convert a linear system into space space form. Always creates + a new system, even if sys is already a StateSpace object. + + ``tf2io(num, den)`` + Create a linear I/O system from its numerator and denominator + polynomial coefficients. + + For details see: :func:`tf` Parameters ---------- - states : int - Number of state variables - inputs : integer - Number of system inputs - outputs : int - Number of system outputs - strictly_proper: bool, optional - If set to 'True', returns a proper system (no direct term). + sys : LTI (StateSpace or TransferFunction) + A linear system. + num : array_like, or list of list of array_like + Polynomial coefficients of the numerator. + den : array_like, or list of list of array_like + Polynomial coefficients of the denominator. Returns ------- - sys : StateSpace - The randomly created linear system + out : LinearIOSystem + New I/O system (in state space form). + + Other Parameters + ---------------- + inputs, outputs : str, or list of str, optional + List of strings that name the individual signals of the transformed + system. If not given, the inputs and outputs are the same as the + original system. + name : string, optional + System name. If unspecified, a generic name is generated + with a unique integer id. Raises ------ ValueError - if any input is not a positive integer + if `num` and `den` have invalid or unequal dimensions, or if an + invalid number of arguments is passed in. + TypeError + if `num` or `den` are of incorrect type, or if sys is not a + TransferFunction object. See Also -------- - rss - - Notes - ----- - If the number of states, inputs, or outputs is not specified, then the - missing numbers are assumed to be 1. The poles of the returned system - will always have a magnitude less than 1. - - """ - # Process states, inputs, outputs (ignoring names) - nstates, _ = _process_signal_list(states) - ninputs, _ = _process_signal_list(inputs) - noutputs, _ = _process_signal_list(outputs) - - sys = _rss_generate( - nstates, ninputs, noutputs, 'd', strictly_proper=strictly_proper) - return LinearIOSystem( - sys, states=states, inputs=inputs, outputs=outputs, **kwargs) - + ss2io + tf2ss -# Convert a state space system into an input/output system (wrapper) -def ss2io(*args, **kwargs): - return LinearIOSystem(*args, **kwargs) -ss2io.__doc__ = LinearIOSystem.__init__.__doc__ + Examples + -------- + >>> num = [[[1., 2.], [3., 4.]], [[5., 6.], [7., 8.]]] + >>> den = [[[9., 8., 7.], [6., 5., 4.]], [[3., 2., 1.], [-1., -2., -3.]]] + >>> sys1 = tf2ss(num, den) + >>> sys_tf = tf(num, den) + >>> sys2 = tf2ss(sys_tf) -# Convert a transfer function into an input/output system (wrapper) -def tf2io(*args, **kwargs): - """Convert a transfer function into an I/O system""" - # TODO: add remaining documentation + """ # Convert the system to a state space system linsys = tf2ss(*args) @@ -2434,11 +2532,9 @@ def tf2io(*args, **kwargs): # Function to create an interconnected system -def interconnect(syslist, connections=None, inplist=[], outlist=[], - inputs=None, outputs=None, states=None, - params={}, dt=None, name=None, +def interconnect(syslist, connections=None, inplist=[], outlist=[], params={}, check_unused=True, ignore_inputs=None, ignore_outputs=None, - **kwargs): + warn_duplicate=None, **kwargs): """Interconnect a set of input/output systems. This function creates a new system that is an interconnection of a set of @@ -2588,6 +2684,12 @@ def interconnect(syslist, connections=None, inplist=[], outlist=[], outputs from all sub-systems with that base name are considered ignored. + warn_duplicate : None, True, or False + Control how warnings are generated if duplicate objects or names are + detected. In `None` (default), then warnings are generated for + systems that have non-generic names. If `False`, warnings are not + generated and if `True` then warnings are always generated. + Example ------- @@ -2610,7 +2712,7 @@ def interconnect(syslist, connections=None, inplist=[], outlist=[], >>> P = control.tf2io(control.tf(1, [1, 0]), inputs='u', outputs='y') >>> C = control.tf2io(control.tf(10, [1, 1]), inputs='e', outputs='u') >>> sumblk = control.summing_junction(inputs=['r', '-y'], output='e') - >>> T = control.interconnect([P, C, sumblk], input='r', output='y') + >>> T = control.interconnect([P, C, sumblk], inputs='r', outputs='y') Notes ----- @@ -2639,9 +2741,9 @@ def interconnect(syslist, connections=None, inplist=[], outlist=[], `outputs`, for more natural naming of SISO systems. """ - # Look for 'input' and 'output' parameter name variants - inputs = _parse_signal_parameter(inputs, 'input', kwargs) - outputs = _parse_signal_parameter(outputs, 'output', kwargs, end=True) + dt = kwargs.pop('dt', None) # by pass normal 'dt' processing + name, inputs, outputs, states, _ = _process_namedio_keywords( + kwargs, end=True) if not check_unused and (ignore_inputs or ignore_outputs): raise ValueError('check_unused is False, but either ' @@ -2658,10 +2760,10 @@ def interconnect(syslist, connections=None, inplist=[], outlist=[], # For each system input, look for outputs with the same name connections = [] for input_sys in syslist: - for input_name in input_sys.input_index.keys(): + for input_name in input_sys.input_labels: connect = [input_sys.name + "." + input_name] for output_sys in syslist: - if input_name in output_sys.output_index.keys(): + if input_name in output_sys.output_labels: connect.append(output_sys.name + "." + input_name) if len(connect) > 1: connections.append(connect) @@ -2673,11 +2775,6 @@ def interconnect(syslist, connections=None, inplist=[], outlist=[], # Use an empty connections list connections = [] - if isinstance(inputs, str): - inputs = [inputs] - if isinstance(outputs, str): - outputs = [outputs] - # If inplist/outlist is not present, try using inputs/outputs instead if not inplist and inputs is not None: inplist = list(inputs) @@ -2743,11 +2840,12 @@ def interconnect(syslist, connections=None, inplist=[], outlist=[], newsys = InterconnectedSystem( syslist, connections=connections, inplist=inplist, outlist=outlist, inputs=inputs, outputs=outputs, states=states, - params=params, dt=dt, name=name) + params=params, dt=dt, name=name, warn_duplicate=warn_duplicate) # check for implicity dropped signals if check_unused: newsys.check_unused_signals(ignore_inputs, ignore_outputs) + # If all subsystems are linear systems, maintain linear structure if all([isinstance(sys, LinearIOSystem) for sys in syslist]): return LinearICSystem(newsys, None) @@ -2757,8 +2855,7 @@ def interconnect(syslist, connections=None, inplist=[], outlist=[], # Summing junction def summing_junction( - inputs=None, output=None, dimension=None, name=None, - prefix='u', **kwargs): + inputs=None, output=None, dimension=None, prefix='u', **kwargs): """Create a summing junction as an input/output system. This function creates a static input/output system that outputs the sum of @@ -2797,10 +2894,10 @@ def summing_junction( Example ------- - >>> P = control.tf2io(ct.tf(1, [1, 0]), input='u', output='y') - >>> C = control.tf2io(ct.tf(10, [1, 1]), input='e', output='u') + >>> P = control.tf2io(ct.tf(1, [1, 0]), inputs='u', outputs='y') + >>> C = control.tf2io(ct.tf(10, [1, 1]), inputs='e', outputs='u') >>> sumblk = control.summing_junction(inputs=['r', '-y'], output='e') - >>> T = control.interconnect((P, C, sumblk), input='r', output='y') + >>> T = control.interconnect((P, C, sumblk), inputs='r', outputs='y') """ # Utility function to parse input and output signal lists @@ -2833,15 +2930,15 @@ def _parse_list(signals, signame='input', prefix='u'): # Return the parsed list return nsignals, names, gains - # Look for 'input' and 'output' parameter name variants - inputs = _parse_signal_parameter(inputs, 'input', kwargs) - output = _parse_signal_parameter(output, 'outputs', kwargs, end=True) - - # Default values for inputs and output + # Parse system and signal names (with some minor pre-processing) + if input is not None: + kwargs['inputs'] = inputs # positional/keyword -> keyword + if output is not None: + kwargs['output'] = output # positional/keyword -> keyword + name, inputs, output, states, dt = _process_namedio_keywords( + kwargs, {'inputs': None, 'outputs': 'y'}, end=True) if inputs is None: raise TypeError("input specification is required") - if output is None: - output = 'y' # Read the input list ninputs, input_names, input_gains = _parse_list( diff --git a/control/lti.py b/control/lti.py index b56c2bb44..fdb4946cd 100644 --- a/control/lti.py +++ b/control/lti.py @@ -16,12 +16,12 @@ from numpy import absolute, real, angle, abs from warnings import warn from . import config +from .namedio import NamedIOSystem, isdtime -__all__ = ['issiso', 'timebase', 'common_timebase', 'timebaseEqual', - 'isdtime', 'isctime', 'pole', 'zero', 'damp', 'evalfr', - 'freqresp', 'dcgain'] +__all__ = ['poles', 'zeros', 'damp', 'evalfr', 'frequency_response', + 'freqresp', 'dcgain', 'pole', 'zero'] -class LTI: +class LTI(NamedIOSystem): """LTI is a parent class to linear time-invariant (LTI) system objects. LTI is the parent to the StateSpace and TransferFunction child classes. It @@ -41,15 +41,13 @@ class LTI: with timebase None can be combined with a system having a specified timebase, and the result will have the timebase of the latter system. - """ + Note: dt processing has been moved to the NamedIOSystem class. - def __init__(self, inputs=1, outputs=1, dt=None): + """ + def __init__(self, inputs=1, outputs=1, states=None, name=None, **kwargs): """Assign the LTI object's numbers of inputs and ouputs.""" - - # Data members common to StateSpace and TransferFunction. - self.ninputs = inputs - self.noutputs = outputs - self.dt = dt + super().__init__( + name=name, inputs=inputs, outputs=outputs, states=states, **kwargs) # # Getter and setter functions for legacy state attributes @@ -77,9 +75,9 @@ def _set_inputs(self, value): """ Deprecated attribute; use :attr:`ninputs` instead. - The ``input`` attribute was used to store the number of system inputs. - It is no longer used. If you need access to the number of inputs for - an LTI system, use :attr:`ninputs`. + The ``inputs`` attribute was used to store the number of system + inputs. It is no longer used. If you need access to the number + of inputs for an LTI system, use :attr:`ninputs`. """) def _get_outputs(self): @@ -100,50 +98,11 @@ def _set_outputs(self, value): """ Deprecated attribute; use :attr:`noutputs` instead. - The ``output`` attribute was used to store the number of system + The ``outputs`` attribute was used to store the number of system outputs. It is no longer used. If you need access to the number of outputs for an LTI system, use :attr:`noutputs`. """) - def isdtime(self, strict=False): - """ - Check to see if a system is a discrete-time system - - Parameters - ---------- - strict: bool, optional - If strict is True, make sure that timebase is not None. Default - is False. - """ - - # If no timebase is given, answer depends on strict flag - if self.dt == None: - return True if not strict else False - - # Look for dt > 0 (also works if dt = True) - return self.dt > 0 - - def isctime(self, strict=False): - """ - Check to see if a system is a continuous-time system - - Parameters - ---------- - sys : LTI system - System to be checked - strict: bool, optional - If strict is True, make sure that timebase is not None. Default - is False. - """ - # If no timebase is given, answer depends on strict flag - if self.dt is None: - return True if not strict else False - return self.dt == 0 - - def issiso(self): - '''Check to see if a system is single input, single output''' - return self.ninputs == 1 and self.noutputs == 1 - def damp(self): '''Natural frequency, damping ratio of system poles @@ -156,9 +115,9 @@ def damp(self): poles : array Array of system poles ''' - poles = self.pole() + poles = self.poles() - if isdtime(self, strict=True): + if self.isdtime(strict=True): splane_poles = np.log(poles.astype(complex))/self.dt else: splane_poles = poles @@ -172,16 +131,16 @@ def frequency_response(self, omega, squeeze=None): Reports the frequency response of the system, - G(j*omega) = mag*exp(j*phase) + G(j*omega) = mag * exp(j*phase) - for continuous time systems. For discrete time systems, the response is - evaluated around the unit circle such that + for continuous time systems. For discrete time systems, the response + is evaluated around the unit circle such that - G(exp(j*omega*dt)) = mag*exp(j*phase). + G(exp(j*omega*dt)) = mag * exp(j*phase). In general the system may be multiple input, multiple output (MIMO), - where `m = self.ninputs` number of inputs and `p = self.noutputs` number - of outputs. + where `m = self.ninputs` number of inputs and `p = self.noutputs` + number of outputs. Parameters ---------- @@ -197,29 +156,37 @@ def frequency_response(self, omega, squeeze=None): Returns ------- - mag : ndarray - The magnitude (absolute value, not dB or log10) of the system - frequency response. If the system is SISO and squeeze is not - True, the array is 1D, indexed by frequency. If the system is not - SISO or squeeze is False, the array is 3D, indexed by the output, + response : :class:`FrequencyReponseData` + Frequency response data object representing the frequency + response. This object can be assigned to a tuple using + + mag, phase, omega = response + + where ``mag`` is the magnitude (absolute value, not dB or + log10) of the system frequency response, ``phase`` is the wrapped + phase in radians of the system frequency response, and ``omega`` + is the (sorted) frequencies at which the response was evaluated. + If the system is SISO and squeeze is not True, ``magnitude`` and + ``phase`` are 1D, indexed by frequency. If the system is not SISO + or squeeze is False, the array is 3D, indexed by the output, input, and frequency. If ``squeeze`` is True then single-dimensional axes are removed. - phase : ndarray - The wrapped phase in radians of the system frequency response. - omega : ndarray - The (sorted) frequencies at which the response was evaluated. """ omega = np.sort(np.array(omega, ndmin=1)) - if isdtime(self, strict=True): + if self.isdtime(strict=True): # Convert the frequency to discrete time if np.any(omega * self.dt > np.pi): warn("__call__: evaluation above Nyquist frequency") s = np.exp(1j * omega * self.dt) else: s = 1j * omega - response = self.__call__(s, squeeze=squeeze) - return abs(response), angle(response), omega + + # Return the data as a frequency response data object + from .frdata import FrequencyResponseData + response = self.__call__(s) + return FrequencyResponseData( + response, omega, return_magphase=True, squeeze=squeeze) def dcgain(self): """Return the zero-frequency gain""" @@ -234,191 +201,22 @@ def _dcgain(self, warn_infinite): else: return zeroresp -# Test to see if a system is SISO -def issiso(sys, strict=False): - """ - Check to see if a system is single input, single output - - Parameters - ---------- - sys : LTI system - System to be checked - strict: bool (default = False) - If strict is True, do not treat scalars as SISO - """ - if isinstance(sys, (int, float, complex, np.number)) and not strict: - return True - elif not isinstance(sys, LTI): - raise ValueError("Object is not an LTI system") - - # Done with the tricky stuff... - return sys.issiso() - -# Return the timebase (with conversion if unspecified) -def timebase(sys, strict=True): - """Return the timebase for an LTI system - - dt = timebase(sys) - - returns the timebase for a system 'sys'. If the strict option is - set to False, dt = True will be returned as 1. - """ - # System needs to be either a constant or an LTI system - if isinstance(sys, (int, float, complex, np.number)): - return None - elif not isinstance(sys, LTI): - raise ValueError("Timebase not defined") - - # Return the sample time, with converstion to float if strict is false - if (sys.dt == None): - return None - elif (strict): - return float(sys.dt) - - return sys.dt - -def common_timebase(dt1, dt2): - """ - Find the common timebase when interconnecting systems - - Parameters - ---------- - dt1, dt2: number or system with a 'dt' attribute (e.g. TransferFunction - or StateSpace system) - - Returns - ------- - dt: number - The common timebase of dt1 and dt2, as specified in - :ref:`conventions-ref`. - - Raises - ------ - ValueError - when no compatible time base can be found - """ - # explanation: - # if either dt is None, they are compatible with anything - # if either dt is True (discrete with unspecified time base), - # use the timebase of the other, if it is also discrete - # otherwise both dts must be equal - if hasattr(dt1, 'dt'): - dt1 = dt1.dt - if hasattr(dt2, 'dt'): - dt2 = dt2.dt - - if dt1 is None: - return dt2 - elif dt2 is None: - return dt1 - elif dt1 is True: - if dt2 > 0: - return dt2 - else: - raise ValueError("Systems have incompatible timebases") - elif dt2 is True: - if dt1 > 0: - return dt1 - else: - raise ValueError("Systems have incompatible timebases") - elif np.isclose(dt1, dt2): - return dt1 - else: - raise ValueError("Systems have incompatible timebases") - -# Check to see if two timebases are equal -def timebaseEqual(sys1, sys2): - """ - Check to see if two systems have the same timebase - - timebaseEqual(sys1, sys2) - - returns True if the timebases for the two systems are compatible. By - default, systems with timebase 'None' are compatible with either - discrete or continuous timebase systems. If two systems have a discrete - timebase (dt > 0) then their timebases must be equal. - """ - warn("timebaseEqual will be deprecated in a future release of " - "python-control; use :func:`common_timebase` instead", - PendingDeprecationWarning) - - if (type(sys1.dt) == bool or type(sys2.dt) == bool): - # Make sure both are unspecified discrete timebases - return type(sys1.dt) == type(sys2.dt) and sys1.dt == sys2.dt - elif (sys1.dt is None or sys2.dt is None): - # One or the other is unspecified => the other can be anything - return True - else: - return sys1.dt == sys2.dt - - -# Check to see if a system is a discrete time system -def isdtime(sys, strict=False): - """ - Check to see if a system is a discrete time system - - Parameters - ---------- - sys : LTI system - System to be checked - strict: bool (default = False) - If strict is True, make sure that timebase is not None - """ - - # Check to see if this is a constant - if isinstance(sys, (int, float, complex, np.number)): - # OK as long as strict checking is off - return True if not strict else False - - # Check for a transfer function or state-space object - if isinstance(sys, LTI): - return sys.isdtime(strict) - - # Check to see if object has a dt object - if hasattr(sys, 'dt'): - # If no timebase is given, answer depends on strict flag - if sys.dt == None: - return True if not strict else False - - # Look for dt > 0 (also works if dt = True) - return sys.dt > 0 - - # Got passed something we don't recognize - return False - -# Check to see if a system is a continuous time system -def isctime(sys, strict=False): - """ - Check to see if a system is a continuous-time system - - Parameters - ---------- - sys : LTI system - System to be checked - strict: bool (default = False) - If strict is True, make sure that timebase is not None - """ - - # Check to see if this is a constant - if isinstance(sys, (int, float, complex, np.number)): - # OK as long as strict checking is off - return True if not strict else False + # + # Deprecated functions + # - # Check for a transfer function or state space object - if isinstance(sys, LTI): - return sys.isctime(strict) + def pole(self): + warn("pole() will be deprecated; use poles()", + PendingDeprecationWarning) + return self.poles() - # Check to see if object has a dt object - if hasattr(sys, 'dt'): - # If no timebase is given, answer depends on strict flag - if sys.dt is None: - return True if not strict else False - return sys.dt == 0 + def zero(self): + warn("zero() will be deprecated; use zeros()", + PendingDeprecationWarning) + return self.zeros() - # Got passed something we don't recognize - return False -def pole(sys): +def poles(sys): """ Compute system poles. @@ -432,23 +230,23 @@ def pole(sys): poles: ndarray Array that contains the system's poles. - Raises - ------ - NotImplementedError - when called on a TransferFunction object - See Also -------- - zero - TransferFunction.pole - StateSpace.pole + zeros + TransferFunction.poles + StateSpace.poles """ - return sys.pole() + return sys.poles() -def zero(sys): +def pole(sys): + warn("pole() will be deprecated; use poles()", PendingDeprecationWarning) + return poles(sys) + + +def zeros(sys): """ Compute system zeros. @@ -462,20 +260,21 @@ def zero(sys): zeros: ndarray Array that contains the system's zeros. - Raises - ------ - NotImplementedError - when called on a MIMO system - See Also -------- - pole - StateSpace.zero - TransferFunction.zero + poles + StateSpace.zeros + TransferFunction.zeros """ - return sys.zero() + return sys.zeros() + + +def zero(sys): + warn("zero() will be deprecated; use zeros()", PendingDeprecationWarning) + return zeros(sys) + def damp(sys, doprint=True): """ @@ -589,7 +388,7 @@ def evalfr(sys, x, squeeze=None): """ return sys.__call__(x, squeeze=squeeze) -def freqresp(sys, omega, squeeze=None): +def frequency_response(sys, omega, squeeze=None): """Frequency response of an LTI system at multiple angular frequencies. In general the system may be multiple input, multiple output (MIMO), where @@ -613,18 +412,21 @@ def freqresp(sys, omega, squeeze=None): Returns ------- - mag : ndarray - The magnitude (absolute value, not dB or log10) of the system - frequency response. If the system is SISO and squeeze is not True, - the array is 1D, indexed by frequency. If the system is not SISO or - squeeze is False, the array is 3D, indexed by the output, input, and + response : FrequencyResponseData + Frequency response data object representing the frequency response. + This object can be assigned to a tuple using + + mag, phase, omega = response + + where ``mag`` is the magnitude (absolute value, not dB or log10) of + the system frequency response, ``phase`` is the wrapped phase in + radians of the system frequency response, and ``omega`` is the + (sorted) frequencies at which the response was evaluated. If the + system is SISO and squeeze is not True, ``magnitude`` and ``phase`` + are 1D, indexed by frequency. If the system is not SISO or squeeze + is False, the array is 3D, indexed by the output, input, and frequency. If ``squeeze`` is True then single-dimensional axes are removed. - phase : ndarray - The wrapped phase in radians of the system frequency response. - omega : ndarray - The list of sorted frequencies at which the response was - evaluated. See Also -------- @@ -663,6 +465,10 @@ def freqresp(sys, omega, squeeze=None): return sys.frequency_response(omega, squeeze=squeeze) +# Alternative name (legacy) +freqresp = frequency_response + + def dcgain(sys): """Return the zero-frequency (or DC) gain of the given system diff --git a/control/margins.py b/control/margins.py index 41739704e..662634086 100644 --- a/control/margins.py +++ b/control/margins.py @@ -52,7 +52,8 @@ import numpy as np import scipy as sp from . import xferfcn -from .lti import issiso, evalfr +from .lti import evalfr +from .namedio import issiso from . import frdata from . import freqplot from .exception import ControlMIMONotImplemented diff --git a/control/matlab/__init__.py b/control/matlab/__init__.py index f10a76c54..80f2a0a65 100644 --- a/control/matlab/__init__.py +++ b/control/matlab/__init__.py @@ -65,6 +65,7 @@ from ..iosys import ss, rss, drss # moved from .statesp from ..xferfcn import * from ..lti import * +from ..namedio import * from ..frdata import * from ..dtime import * from ..exception import ControlArgument @@ -84,6 +85,9 @@ from ..dtime import c2d from ..sisotool import sisotool +# Functions that are renamed in MATLAB +pole, zero = poles, zeros + # Import functions specific to Matlab compatibility package from .timeresp import * from .wrappers import * diff --git a/control/modelsimp.py b/control/modelsimp.py index f43acc2fd..432b76b96 100644 --- a/control/modelsimp.py +++ b/control/modelsimp.py @@ -45,7 +45,7 @@ import warnings from .exception import ControlSlycot, ControlMIMONotImplemented, \ ControlDimension -from .lti import isdtime, isctime +from .namedio import isdtime, isctime from .statesp import StateSpace from .statefbk import gram @@ -354,7 +354,7 @@ def minreal(sys, tol=None, verbose=True): sysr = sys.minreal(tol) if verbose: print("{nstates} states have been removed from the model".format( - nstates=len(sys.pole()) - len(sysr.pole()))) + nstates=len(sys.poles()) - len(sysr.poles()))) return sysr diff --git a/control/namedio.py b/control/namedio.py index 4ea82d819..254f310ff 100644 --- a/control/namedio.py +++ b/control/namedio.py @@ -1,25 +1,21 @@ -# namedio.py - internal named I/O object class +# namedio.py - named I/O system class and helper functions # RMM, 13 Mar 2022 # -# This file implements the _NamedIOSystem and _NamedIOStateSystem classes, -# which are used as a parent classes for FrequencyResponseData, -# InputOutputSystem, LTI, TimeResponseData, and other similar classes to -# allow naming of signals. +# This file implements the NamedIOSystem class, which is used as a parent +# class for FrequencyResponseData, InputOutputSystem, LTI, TimeResponseData, +# and other similar classes to allow naming of signals. import numpy as np +from copy import copy +from warnings import warn +from . import config +__all__ = ['issiso', 'timebase', 'common_timebase', 'timebaseEqual', + 'isdtime', 'isctime'] -class _NamedIOSystem(object): - _idCounter = 0 - - def _name_or_default(self, name=None): - if name is None: - name = "sys[{}]".format(_NamedIOSystem._idCounter) - _NamedIOSystem._idCounter += 1 - return name - +class NamedIOSystem(object): def __init__( - self, inputs=None, outputs=None, name=None): + self, name=None, inputs=None, outputs=None, states=None, **kwargs): # system name self.name = self._name_or_default(name) @@ -27,6 +23,31 @@ def __init__( # Parse and store the number of inputs and outputs self.set_inputs(inputs) self.set_outputs(outputs) + self.set_states(states) + + # Process timebase: if not given use default, but allow None as value + self.dt = _process_dt_keyword(kwargs) + + # Make sure there were no other keywords + if kwargs: + raise TypeError("unrecognized keywords: ", str(kwargs)) + + # + # Functions to manipulate the system name + # + _idCounter = 0 # Counter for creating generic system name + + # Return system name + def _name_or_default(self, name=None): + if name is None: + name = "sys[{}]".format(NamedIOSystem._idCounter) + NamedIOSystem._idCounter += 1 + return name + + # Check if system name is generic + def _generic_name_check(self): + import re + return re.match(r'^sys\[\d*\]$', self.name) is not None # # Class attributes @@ -38,33 +59,64 @@ def __init__( #: Number of system inputs. #: #: :meta hide-value: - ninputs = 0 + ninputs = None #: Number of system outputs. #: #: :meta hide-value: - noutputs = 0 + noutputs = None + + #: Number of system states. + #: + #: :meta hide-value: + nstates = None def __repr__(self): - return str(type(self)) + ": " + self.name if self.name is not None \ - else str(type(self)) + return f'<{self.__class__.__name__}:{self.name}:' + \ + f'{list(self.input_labels)}->{list(self.output_labels)}>' def __str__(self): """String representation of an input/output object""" - str = "Object: " + (self.name if self.name else "(None)") + "\n" - str += "Inputs (%s): " % self.ninputs - for key in self.input_index: - str += key + ", " - str += "\nOutputs (%s): " % self.noutputs - for key in self.output_index: - str += key + ", " + str = f"<{self.__class__.__name__}>: {self.name}\n" + str += f"Inputs ({self.ninputs}): {self.input_labels}\n" + str += f"Outputs ({self.noutputs}): {self.output_labels}\n" + if self.nstates is not None: + str += f"States ({self.nstates}): {self.state_labels}" return str # Find a signal by name def _find_signal(self, name, sigdict): return sigdict.get(name, None) + def copy(self, name=None, use_prefix_suffix=True): + """Make a copy of an input/output system + + A copy of the system is made, with a new name. The `name` keyword + can be used to specify a specific name for the system. If no name + is given and `use_prefix_suffix` is True, the name is constructed + by prepending config.defaults['iosys.duplicate_system_name_prefix'] + and appending config.defaults['iosys.duplicate_system_name_suffix']. + Otherwise, a generic system name of the form `sys[]` is used, + where `` is based on an internal counter. + + """ + # Create a copy of the system + newsys = copy(self) + + # Update the system name + if name is None and use_prefix_suffix: + # Get the default prefix and suffix to use + dup_prefix = config.defaults['iosys.duplicate_system_name_prefix'] + dup_suffix = config.defaults['iosys.duplicate_system_name_suffix'] + newsys.name = self._name_or_default( + dup_prefix + self.name + dup_suffix) + else: + newsys.name = self._name_or_default(name) + + return newsys + def set_inputs(self, inputs, prefix='u'): + """Set the number/names of the system inputs. Parameters @@ -122,44 +174,6 @@ def find_output(self, name): lambda self: list(self.output_index.keys()), # getter set_outputs) # setter - def issiso(self): - """Check to see if a system is single input, single output""" - return self.ninputs == 1 and self.noutputs == 1 - - -class _NamedIOStateSystem(_NamedIOSystem): - def __init__( - self, inputs=None, outputs=None, states=None, name=None): - # Parse and store the system name, inputs, and outputs - super().__init__(inputs=inputs, outputs=outputs, name=name) - - # Parse and store the number of states - self.set_states(states) - - # - # Class attributes - # - # These attributes are defined as class attributes so that they are - # documented properly. They are "overwritten" in __init__. - # - - #: Number of system states. - #: - #: :meta hide-value: - nstates = 0 - - def __str__(self): - """String representation of an input/output system""" - str = _NamedIOSystem.__str__(self) - str += "\nStates (%s): " % self.nstates - for key in self.state_index: - str += key + ", " - return str - - def _isstatic(self): - """Check to see if a system is a static system (no states)""" - return self.nstates == 0 - def set_states(self, states, prefix='x'): """Set the number/names of the system states. @@ -189,6 +203,344 @@ def find_state(self, name): lambda self: list(self.state_index.keys()), # getter set_states) # setter + def isctime(self, strict=False): + """ + Check to see if a system is a continuous-time system + + Parameters + ---------- + sys : Named I/O system + System to be checked + strict: bool, optional + If strict is True, make sure that timebase is not None. Default + is False. + """ + # If no timebase is given, answer depends on strict flag + if self.dt is None: + return True if not strict else False + return self.dt == 0 + + def isdtime(self, strict=False): + """ + Check to see if a system is a discrete-time system + + Parameters + ---------- + strict: bool, optional + If strict is True, make sure that timebase is not None. Default + is False. + """ + + # If no timebase is given, answer depends on strict flag + if self.dt == None: + return True if not strict else False + + # Look for dt > 0 (also works if dt = True) + return self.dt > 0 + + def issiso(self): + """Check to see if a system is single input, single output""" + return self.ninputs == 1 and self.noutputs == 1 + + def _isstatic(self): + """Check to see if a system is a static system (no states)""" + return self.nstates == 0 + + +# Test to see if a system is SISO +def issiso(sys, strict=False): + """ + Check to see if a system is single input, single output + + Parameters + ---------- + sys : I/O or LTI system + System to be checked + strict: bool (default = False) + If strict is True, do not treat scalars as SISO + """ + if isinstance(sys, (int, float, complex, np.number)) and not strict: + return True + elif not isinstance(sys, NamedIOSystem): + raise ValueError("Object is not an I/O or LTI system") + + # Done with the tricky stuff... + return sys.issiso() + +# Return the timebase (with conversion if unspecified) +def timebase(sys, strict=True): + """Return the timebase for a system + + dt = timebase(sys) + + returns the timebase for a system 'sys'. If the strict option is + set to False, dt = True will be returned as 1. + """ + # System needs to be either a constant or an I/O or LTI system + if isinstance(sys, (int, float, complex, np.number)): + return None + elif not isinstance(sys, NamedIOSystem): + raise ValueError("Timebase not defined") + + # Return the sample time, with converstion to float if strict is false + if (sys.dt == None): + return None + elif (strict): + return float(sys.dt) + + return sys.dt + +def common_timebase(dt1, dt2): + """ + Find the common timebase when interconnecting systems + + Parameters + ---------- + dt1, dt2: number or system with a 'dt' attribute (e.g. TransferFunction + or StateSpace system) + + Returns + ------- + dt: number + The common timebase of dt1 and dt2, as specified in + :ref:`conventions-ref`. + + Raises + ------ + ValueError + when no compatible time base can be found + """ + # explanation: + # if either dt is None, they are compatible with anything + # if either dt is True (discrete with unspecified time base), + # use the timebase of the other, if it is also discrete + # otherwise both dts must be equal + if hasattr(dt1, 'dt'): + dt1 = dt1.dt + if hasattr(dt2, 'dt'): + dt2 = dt2.dt + + if dt1 is None: + return dt2 + elif dt2 is None: + return dt1 + elif dt1 is True: + if dt2 > 0: + return dt2 + else: + raise ValueError("Systems have incompatible timebases") + elif dt2 is True: + if dt1 > 0: + return dt1 + else: + raise ValueError("Systems have incompatible timebases") + elif np.isclose(dt1, dt2): + return dt1 + else: + raise ValueError("Systems have incompatible timebases") + +# Check to see if two timebases are equal +def timebaseEqual(sys1, sys2): + """ + Check to see if two systems have the same timebase + + timebaseEqual(sys1, sys2) + + returns True if the timebases for the two systems are compatible. By + default, systems with timebase 'None' are compatible with either + discrete or continuous timebase systems. If two systems have a discrete + timebase (dt > 0) then their timebases must be equal. + """ + warn("timebaseEqual will be deprecated in a future release of " + "python-control; use :func:`common_timebase` instead", + PendingDeprecationWarning) + + if (type(sys1.dt) == bool or type(sys2.dt) == bool): + # Make sure both are unspecified discrete timebases + return type(sys1.dt) == type(sys2.dt) and sys1.dt == sys2.dt + elif (sys1.dt is None or sys2.dt is None): + # One or the other is unspecified => the other can be anything + return True + else: + return sys1.dt == sys2.dt + + +# Check to see if a system is a discrete time system +def isdtime(sys, strict=False): + """ + Check to see if a system is a discrete time system + + Parameters + ---------- + sys : I/O or LTI system + System to be checked + strict: bool (default = False) + If strict is True, make sure that timebase is not None + """ + + # Check to see if this is a constant + if isinstance(sys, (int, float, complex, np.number)): + # OK as long as strict checking is off + return True if not strict else False + + # Check for a transfer function or state-space object + if isinstance(sys, NamedIOSystem): + return sys.isdtime(strict) + + # Check to see if object has a dt object + if hasattr(sys, 'dt'): + # If no timebase is given, answer depends on strict flag + if sys.dt == None: + return True if not strict else False + + # Look for dt > 0 (also works if dt = True) + return sys.dt > 0 + + # Got passed something we don't recognize + return False + +# Check to see if a system is a continuous time system +def isctime(sys, strict=False): + """ + Check to see if a system is a continuous-time system + + Parameters + ---------- + sys : I/O or LTI system + System to be checked + strict: bool (default = False) + If strict is True, make sure that timebase is not None + """ + + # Check to see if this is a constant + if isinstance(sys, (int, float, complex, np.number)): + # OK as long as strict checking is off + return True if not strict else False + + # Check for a transfer function or state space object + if isinstance(sys, NamedIOSystem): + return sys.isctime(strict) + + # Check to see if object has a dt object + if hasattr(sys, 'dt'): + # If no timebase is given, answer depends on strict flag + if sys.dt is None: + return True if not strict else False + return sys.dt == 0 + + # Got passed something we don't recognize + return False + + +# Utility function to parse nameio keywords +def _process_namedio_keywords( + keywords={}, defaults={}, static=False, end=False): + """Process namedio specification + + This function processes the standard keywords used in initializing a named + I/O system. It first looks in the `keyword` dictionary to see if a value + is specified. If not, the `default` dictionary is used. The `default` + dictionary can also be set to a NamedIOSystem object, which is useful for + copy constructors that change system and signal names. + + If `end` is True, then generate an error if there are any remaining + keywords. + + """ + # If default is a system, redefine as a dictionary + if isinstance(defaults, NamedIOSystem): + sys = defaults + defaults = { + 'name': sys.name, 'inputs': sys.input_labels, + 'outputs': sys.output_labels, 'dt': sys.dt} + + if sys.nstates is not None: + defaults['states'] = sys.state_labels + + elif not isinstance(defaults, dict): + raise TypeError("default must be dict or sys") + + else: + sys = None + + # Sort out singular versus plural signal names + for singular in ['input', 'output', 'state']: + kw = singular + 's' + if singular in keywords and kw in keywords: + raise TypeError(f"conflicting keywords '{singular}' and '{kw}'") + + if singular in keywords: + keywords[kw] = keywords.pop(singular) + + # Utility function to get keyword with defaults, processing + def pop_with_default(kw, defval=None, return_list=True): + val = keywords.pop(kw, None) + if val is None: + val = defaults.get(kw, defval) + if return_list and isinstance(val, str): + val = [val] # make sure to return a list + return val + + # Process system and signal names + name = pop_with_default('name', return_list=False) + inputs = pop_with_default('inputs') + outputs = pop_with_default('outputs') + states = pop_with_default('states') + + # If we were given a system, make sure sizes match list lengths + if sys: + if isinstance(inputs, list) and sys.ninputs != len(inputs): + raise ValueError("Wrong number of input labels given.") + if isinstance(outputs, list) and sys.noutputs != len(outputs): + raise ValueError("Wrong number of output labels given.") + if sys.nstates is not None and \ + isinstance(states, list) and sys.nstates != len(states): + raise ValueError("Wrong number of state labels given.") + + # Process timebase: if not given use default, but allow None as value + dt = _process_dt_keyword(keywords, defaults, static=static) + + # If desired, make sure we processed all keywords + if end and keywords: + raise TypeError("unrecognized keywords: ", str(keywords)) + + # Return the processed keywords + return name, inputs, outputs, states, dt + +# +# Parse 'dt' in for named I/O system +# +# The 'dt' keyword is used to set the timebase for a system. Its +# processing is a bit unusual: if it is not specified at all, then the +# value is pulled from config.defaults['control.default_dt']. But +# since 'None' is an allowed value, we can't just use the default if +# dt is None. Instead, we have to look to see if it was listed as a +# variable keyword. +# +# In addition, if a system is static and dt is not specified, we set dt = +# None to allow static systems to be combined with either discrete-time or +# continuous-time systems. +# +# TODO: update all 'dt' processing to call this function, so that +# everything is done consistently. +# +def _process_dt_keyword(keywords, defaults={}, static=False): + if static and 'dt' not in keywords and 'dt' not in defaults: + dt = None + elif 'dt' in keywords: + dt = keywords.pop('dt') + elif 'dt' in defaults: + dt = defaults.pop('dt') + else: + dt = config.defaults['control.default_dt'] + + # Make sure that the value for dt is valid + if dt is not None and not isinstance(dt, (bool, int, float)) or \ + isinstance(dt, (bool, int, float)) and dt < 0: + raise ValueError(f"invalid timebase, dt = {dt}") + + return dt + # Utility function to parse a list of signals def _process_signal_list(signals, prefix='s'): diff --git a/control/pzmap.py b/control/pzmap.py index 7d3836d7f..09f58b79c 100644 --- a/control/pzmap.py +++ b/control/pzmap.py @@ -41,7 +41,8 @@ from numpy import real, imag, linspace, exp, cos, sin, sqrt from math import pi -from .lti import LTI, isdtime, isctime +from .lti import LTI +from .namedio import isdtime, isctime from .grid import sgrid, zgrid, nogrid from . import config @@ -104,8 +105,8 @@ def pzmap(sys, plot=None, grid=None, title='Pole Zero Map', **kwargs): if not isinstance(sys, LTI): raise TypeError('Argument ``sys``: must be a linear system.') - poles = sys.pole() - zeros = sys.zero() + poles = sys.poles() + zeros = sys.zeros() if (plot): import matplotlib.pyplot as plt diff --git a/control/rlocus.py b/control/rlocus.py index 5cf7983a3..9d531de94 100644 --- a/control/rlocus.py +++ b/control/rlocus.py @@ -55,7 +55,7 @@ import matplotlib.pyplot as plt from numpy import array, poly1d, row_stack, zeros_like, real, imag import scipy.signal # signal processing toolbox -from .lti import isdtime +from .namedio import isdtime from .xferfcn import _convert_to_transfer_function from .exception import ControlMIMONotImplemented from .sisotool import _SisotoolUpdate diff --git a/control/sisotool.py b/control/sisotool.py index 41f21ecbe..52c061249 100644 --- a/control/sisotool.py +++ b/control/sisotool.py @@ -3,14 +3,12 @@ from control.exception import ControlMIMONotImplemented from .freqplot import bode_plot from .timeresp import step_response -from .lti import issiso, isdtime +from .namedio import issiso, common_timebase, isctime, isdtime from .xferfcn import tf from .iosys import ss from .bdalg import append, connect from .iosys import tf2io, ss2io, summing_junction, interconnect from control.statesp import _convert_to_statespace, StateSpace -from control.lti import common_timebase, isctime -import matplotlib import matplotlib.pyplot as plt import warnings diff --git a/control/statefbk.py b/control/statefbk.py index 0aaf49f61..97f314da5 100644 --- a/control/statefbk.py +++ b/control/statefbk.py @@ -45,7 +45,8 @@ from . import statesp from .mateqn import care, dare, _check_shape from .statesp import StateSpace, _ssmatrix, _convert_to_statespace -from .lti import LTI, isdtime, isctime +from .lti import LTI +from .namedio import isdtime, isctime from .iosys import InputOutputSystem, NonlinearIOSystem, LinearIOSystem, \ interconnect, ss from .exception import ControlSlycot, ControlArgument, ControlDimension, \ diff --git a/control/statesp.py b/control/statesp.py index 435ff702f..58412e57a 100644 --- a/control/statesp.py +++ b/control/statesp.py @@ -58,8 +58,10 @@ from scipy.signal import cont2discrete from scipy.signal import StateSpace as signalStateSpace from warnings import warn -from .lti import LTI, common_timebase, isdtime, _process_frequency_response -from .namedio import _NamedIOStateSystem, _process_signal_list +from .frdata import FrequencyResponseData +from .lti import LTI, _process_frequency_response +from .namedio import common_timebase, isdtime +from .namedio import _process_namedio_keywords from . import config from copy import deepcopy @@ -153,7 +155,7 @@ def _f2s(f): return s -class StateSpace(LTI, _NamedIOStateSystem): +class StateSpace(LTI): """StateSpace(A, B, C, D[, dt]) A class for representing state-space models. @@ -161,7 +163,9 @@ class StateSpace(LTI, _NamedIOStateSystem): The StateSpace class is used to represent state-space realizations of linear time-invariant (LTI) systems: + .. math:: dx/dt = A x + B u + y = C x + D u where u is the input, y is the output, and x is the state. @@ -217,6 +221,8 @@ class StateSpace(LTI, _NamedIOStateSystem): The default value of dt can be changed by changing the value of ``control.config.defaults['control.default_dt']``. + Note: timebase processing has moved to namedio. + A state space system is callable and returns the value of the transfer function evaluated at a point in the complex plane. See :meth:`~control.StateSpace.__call__` for a more detailed description. @@ -244,7 +250,7 @@ class StateSpace(LTI, _NamedIOStateSystem): # Allow ndarray * StateSpace to give StateSpace._rmul_() priority __array_priority__ = 11 # override ndarray and matrix types - def __init__(self, *args, keywords=None, **kwargs): + def __init__(self, *args, init_namedio=True, **kwargs): """StateSpace(A, B, C, D[, dt]) Construct a state space object. @@ -262,18 +268,27 @@ def __init__(self, *args, keywords=None, **kwargs): value is read from `config.defaults['statesp.remove_useless_states']` (default = False). - """ - # Use keywords object if we received one (and pop keywords we use) - if keywords is None: - keywords = kwargs + The `init_namedio` keyword can be used to turn off initialization of + system and signal names. This is used internally by the + :class:`LinearIOSystem` class to avoid renaming. - # first get A, B, C, D matrices + """ + # + # Process positional arguments + # if len(args) == 4: # The user provided A, B, C, and D matrices. (A, B, C, D) = args + elif len(args) == 5: # Discrete time system - (A, B, C, D, _) = args + (A, B, C, D, dt) = args + if 'dt' in kwargs: + warn("received multiple dt arguments, " + "using positional arg dt = %s" % dt) + kwargs['dt'] = dt + args = args[:-1] + elif len(args) == 1: # Use the copy constructor. if not isinstance(args[0], StateSpace): @@ -284,15 +299,11 @@ def __init__(self, *args, keywords=None, **kwargs): B = args[0].B C = args[0].C D = args[0].D + else: - raise ValueError( + raise TypeError( "Expected 1, 4, or 5 arguments; received %i." % len(args)) - # Process keyword arguments - remove_useless_states = keywords.pop( - 'remove_useless_states', - config.defaults['statesp.remove_useless_states']) - # Convert all matrices to standard form A = _ssmatrix(A) # if B is a 1D array, turn it into a column vector if it fits @@ -309,41 +320,38 @@ def __init__(self, *args, keywords=None, **kwargs): D = np.zeros((C.shape[0], B.shape[1])) D = _ssmatrix(D) - super().__init__(inputs=D.shape[1], outputs=D.shape[0]) + # Matrices definining the linear system self.A = A self.B = B self.C = C self.D = D - # now set dt - if len(args) == 4: - if 'dt' in keywords: - dt = keywords.pop('dt') - elif self._isstatic(): - dt = None - else: - dt = config.defaults['control.default_dt'] - elif len(args) == 5: - dt = args[4] - if 'dt' in keywords: - warn("received multiple dt arguments, " - "using positional arg dt = %s" % dt) - keywords.pop('dt') - elif len(args) == 1: - try: - dt = args[0].dt - except AttributeError: - if self._isstatic(): - dt = None - else: - dt = config.defaults['control.default_dt'] - self.dt = dt - self.nstates = A.shape[1] - - # Make sure there were no extraneous keywords - if keywords: - raise TypeError("unrecognized keywords: ", str(keywords)) + # + # Process keyword arguments + # + remove_useless_states = kwargs.pop( + 'remove_useless_states', + config.defaults['statesp.remove_useless_states']) + + # Initialize the instance variables + if init_namedio: + # Process namedio keywords + defaults = args[0] if len(args) == 1 else \ + {'inputs': D.shape[1], 'outputs': D.shape[0], + 'states': A.shape[0]} + static = (A.size == 0) + name, inputs, outputs, states, dt = _process_namedio_keywords( + kwargs, defaults, static=static, end=True) + + # Initialize LTI (NamedIOSystem) object + super().__init__( + name=name, inputs=inputs, outputs=outputs, + states=states, dt=dt) + elif kwargs: + raise TypeError("unrecognized keyword(s): ", str(kwargs)) + + # Reset shapes (may not be needed once np.matrix support is removed) if 0 == self.nstates: # static gain # matrix's default "empty" shape is 1x0 @@ -351,8 +359,11 @@ def __init__(self, *args, keywords=None, **kwargs): B.shape = (0, self.ninputs) C.shape = (self.noutputs, 0) - # Check that the matrix sizes are consistent. - if self.nstates != A.shape[0]: + # + # Check to make sure everything is consistent + # + # Check that the matrix sizes are consistent + if A.shape[0] != A.shape[1] or self.nstates != A.shape[0]: raise ValueError("A must be square.") if self.nstates != B.shape[0]: raise ValueError("A and B must have the same number of rows.") @@ -363,7 +374,10 @@ def __init__(self, *args, keywords=None, **kwargs): if self.noutputs != C.shape[0]: raise ValueError("C and D must have the same number of rows.") - # Check for states that don't do anything, and remove them. + # + # Final processing + # + # Check for states that don't do anything, and remove them if remove_useless_states: self._remove_useless_states() @@ -464,9 +478,10 @@ def _remove_useless_states(self): self.B = delete(self.B, useless, 0) self.C = delete(self.C, useless, 1) - self.nstates = self.A.shape[0] - self.ninputs = self.B.shape[1] - self.noutputs = self.C.shape[0] + # Remove any state names that we don't need + self.set_states( + [self.state_labels[i] for i in range(self.nstates) + if i not in useless]) def __str__(self): """Return string representation of the state space system.""" @@ -655,6 +670,13 @@ def __add__(self, other): D = self.D + other dt = self.dt else: + # Check to see if the right operator has priority + if getattr(other, '__array_priority__', None) and \ + getattr(self, '__array_priority__', None) and \ + other.__array_priority__ > self.__array_priority__: + return other.__radd__(self) + + # Convert the other argument to state space other = _convert_to_statespace(other) # Check to make sure the dimensions are OK @@ -705,6 +727,13 @@ def __mul__(self, other): D = self.D * other dt = self.dt else: + # Check to see if the right operator has priority + if getattr(other, '__array_priority__', None) and \ + getattr(self, '__array_priority__', None) and \ + other.__array_priority__ > self.__array_priority__: + return other.__rmul__(self) + + # Convert the other argument to state space other = _convert_to_statespace(other) # Check to make sure the dimensions are OK @@ -933,7 +962,7 @@ def horner(self, x, warn_infinite=True): # Evaluating at a pole. Return value depends if there # is a zero at the same point or not. - if x_idx in self.zero(): + if x_idx in self.zeros(): out[:, :, idx] = complex(np.nan, np.nan) else: out[:, :, idx] = complex(np.inf, np.nan) @@ -955,12 +984,12 @@ def freqresp(self, omega): return self.frequency_response(omega) # Compute poles and zeros - def pole(self): + def poles(self): """Compute the poles of a state space system.""" return eigvals(self.A) if self.nstates else np.array([]) - def zero(self): + def zeros(self): """Compute the zeros of a state space system.""" if not self.nstates: @@ -982,9 +1011,9 @@ def zero(self): except ImportError: # Slycot unavailable. Fall back to scipy. if self.C.shape[0] != self.D.shape[1]: - raise NotImplementedError("StateSpace.zero only supports " - "systems with the same number of " - "inputs as outputs.") + raise NotImplementedError( + "StateSpace.zero only supports systems with the same " + "number of inputs as outputs.") # This implements the QZ algorithm for finding transmission zeros # from @@ -1369,7 +1398,7 @@ def dynamics(self, t, x, u=None): The first argument `t` is ignored because :class:`StateSpace` systems are time-invariant. It is included so that the dynamics can be passed - to most numerical integrators, such as :func:`scipy.integrate.solve_ivp` + to numerical integrators, such as :func:`scipy.integrate.solve_ivp` and for consistency with :class:`IOSystem` systems. Parameters @@ -1384,6 +1413,7 @@ def dynamics(self, t, x, u=None): Returns ------- dx/dt or x[t+dt] : ndarray + """ x = np.reshape(x, (-1, 1)) # force to a column in case matrix if np.size(x) != self.nstates: @@ -1447,7 +1477,7 @@ def _isstatic(self): # TODO: add discrete time check -def _convert_to_statespace(sys, **kw): +def _convert_to_statespace(sys): """Convert a system to state space form (if needed). If sys is already a state space, then it is returned. If sys is a @@ -1460,16 +1490,15 @@ def _convert_to_statespace(sys, **kw): In the latter example, A = B = C = 0 and D = [[1., 1., 1.] [1., 1., 1.]]. + + Note: no renaming of inputs and outputs is performed; this should be done + by the calling function. + """ from .xferfcn import TransferFunction import itertools if isinstance(sys, StateSpace): - if len(kw): - raise TypeError("If sys is a StateSpace, _convert_to_statespace " - "cannot take keywords.") - - # Already a state space system; just return it return sys elif isinstance(sys, TransferFunction): @@ -1478,11 +1507,9 @@ def _convert_to_statespace(sys, **kw): [[len(num) for num in col] for col in sys.den]): raise ValueError("Transfer function is non-proper; can't " "convert to StateSpace system.") + try: from slycot import td04ad - if len(kw): - raise TypeError("If sys is a TransferFunction, " - "_convert_to_statespace cannot take keywords.") # Change the numerator and denominator arrays so that the transfer # function matrix has a common denominator. @@ -1494,10 +1521,10 @@ def _convert_to_statespace(sys, **kw): denorder, den, num, tol=0) states = ssout[0] - return StateSpace(ssout[1][:states, :states], - ssout[2][:states, :sys.ninputs], - ssout[3][:sys.noutputs, :states], ssout[4], - sys.dt) + return StateSpace( + ssout[1][:states, :states], ssout[2][:states, :sys.ninputs], + ssout[3][:sys.noutputs, :states], ssout[4], sys.dt, + inputs=sys.input_labels, outputs=sys.output_labels) except ImportError: # No Slycot. Scipy tf->ss can't handle MIMO, but static # MIMO is an easy special case we can check for here @@ -1520,34 +1547,25 @@ def _convert_to_statespace(sys, **kw): # the squeeze A, B, C, D = \ sp.signal.tf2ss(squeeze(sys.num), squeeze(sys.den)) - return StateSpace(A, B, C, D, sys.dt) + return StateSpace( + A, B, C, D, sys.dt, inputs=sys.input_labels, + outputs=sys.output_labels) - elif isinstance(sys, (int, float, complex, np.number)): - if "inputs" in kw: - inputs = kw["inputs"] - else: - inputs = 1 - if "outputs" in kw: - outputs = kw["outputs"] - else: - outputs = 1 - - # Generate a simple state space system of the desired dimension - # The following Doesn't work due to inconsistencies in ltisys: - # return StateSpace([[]], [[]], [[]], eye(outputs, inputs)) - return StateSpace([], zeros((0, inputs)), zeros((outputs, 0)), - sys * ones((outputs, inputs))) + elif isinstance(sys, FrequencyResponseData): + raise TypeError("Can't convert FRD to StateSpace system.") # If this is a matrix, try to create a constant feedthrough try: - D = _ssmatrix(sys) - return StateSpace([], [], [], D) + D = _ssmatrix(np.atleast_2d(sys)) + return StateSpace([], [], [], D, dt=None) + except Exception: raise TypeError("Can't convert given type to StateSpace system.") # TODO: add discrete time option -def _rss_generate(states, inputs, outputs, cdtype, strictly_proper=False): +def _rss_generate( + states, inputs, outputs, cdtype, strictly_proper=False, name=None): """Generate a random state space. This does the actual random state space generation expected from rss and @@ -1665,7 +1683,7 @@ def _rss_generate(states, inputs, outputs, cdtype, strictly_proper=False): ss_args = (A, B, C, D) else: ss_args = (A, B, C, D, True) - return StateSpace(*ss_args) + return StateSpace(*ss_args, name=name) # Convert a MIMO system to a SISO system @@ -1776,31 +1794,7 @@ def _mimo2simo(sys, input, warn_conversion=False): return sys -def _ss(*args, keywords=None, **kwargs): - """Internal function to create StateSpace system""" - if len(args) == 4 or len(args) == 5: - return StateSpace(*args, keywords=keywords, **kwargs) - - elif len(args) == 1: - # Make sure there were no extraneous keywords - if kwargs: - raise TypeError("unrecognized keywords: ", str(kwargs)) - - from .xferfcn import TransferFunction - sys = args[0] - if isinstance(sys, StateSpace): - return deepcopy(sys) - elif isinstance(sys, TransferFunction): - return tf2ss(sys) - else: - raise TypeError("ss(sys): sys must be a StateSpace or " - "TransferFunction object. It is %s." % type(sys)) - else: - raise ValueError( - "Needs 1, 4, or 5 arguments; received %i." % len(args)) - - -def tf2ss(*args): +def tf2ss(*args, **kwargs): """tf2ss(sys) Transform a transfer function to a state space system. @@ -1808,11 +1802,11 @@ def tf2ss(*args): The function accepts either 1 or 2 parameters: ``tf2ss(sys)`` - Convert a linear system into transfer function form. Always creates - a new system, even if sys is already a TransferFunction object. + Convert a linear system into space space form. Always creates + a new system, even if sys is already a StateSpace object. ``tf2ss(num, den)`` - Create a transfer function system from its numerator and denominator + Create a state space system from its numerator and denominator polynomial coefficients. For details see: :func:`tf` @@ -1831,6 +1825,16 @@ def tf2ss(*args): out : StateSpace New linear system in state space form + Other Parameters + ---------------- + inputs, outputs : str, or list of str, optional + List of strings that name the individual signals of the transformed + system. If not given, the inputs and outputs are the same as the + original system. + name : string, optional + System name. If unspecified, a generic name is generated + with a unique integer id. + Raises ------ ValueError @@ -1860,14 +1864,15 @@ def tf2ss(*args): from .xferfcn import TransferFunction if len(args) == 2 or len(args) == 3: # Assume we were given the num, den - return _convert_to_statespace(TransferFunction(*args)) + return StateSpace( + _convert_to_statespace(TransferFunction(*args)), **kwargs) elif len(args) == 1: sys = args[0] if not isinstance(sys, TransferFunction): raise TypeError("tf2ss(sys): sys must be a TransferFunction " "object.") - return _convert_to_statespace(sys) + return StateSpace(_convert_to_statespace(sys), **kwargs) else: raise ValueError("Needs 1 or 2 arguments; received %i." % len(args)) diff --git a/control/stochsys.py b/control/stochsys.py index fd276b92c..2b8233070 100644 --- a/control/stochsys.py +++ b/control/stochsys.py @@ -21,7 +21,8 @@ from math import sqrt from .iosys import InputOutputSystem, NonlinearIOSystem -from .lti import LTI, isctime, isdtime +from .lti import LTI +from .namedio import isctime, isdtime from .mateqn import care, dare, _check_shape from .statesp import StateSpace, _ssmatrix from .exception import ControlArgument, ControlNotImplemented diff --git a/control/tests/bdalg_test.py b/control/tests/bdalg_test.py index 433a584cc..2f6b5523f 100644 --- a/control/tests/bdalg_test.py +++ b/control/tests/bdalg_test.py @@ -11,7 +11,7 @@ from control.xferfcn import TransferFunction from control.statesp import StateSpace from control.bdalg import feedback, append, connect -from control.lti import zero, pole +from control.lti import zeros, poles class TestFeedback: @@ -188,52 +188,52 @@ def testLists(self, tsys): # Series sys1_2 = ctrl.series(sys1, sys2) - np.testing.assert_array_almost_equal(sort(pole(sys1_2)), [-4., -2.]) - np.testing.assert_array_almost_equal(sort(zero(sys1_2)), [-3., -1.]) + np.testing.assert_array_almost_equal(sort(poles(sys1_2)), [-4., -2.]) + np.testing.assert_array_almost_equal(sort(zeros(sys1_2)), [-3., -1.]) sys1_3 = ctrl.series(sys1, sys2, sys3) - np.testing.assert_array_almost_equal(sort(pole(sys1_3)), + np.testing.assert_array_almost_equal(sort(poles(sys1_3)), [-6., -4., -2.]) - np.testing.assert_array_almost_equal(sort(zero(sys1_3)), + np.testing.assert_array_almost_equal(sort(zeros(sys1_3)), [-5., -3., -1.]) sys1_4 = ctrl.series(sys1, sys2, sys3, sys4) - np.testing.assert_array_almost_equal(sort(pole(sys1_4)), + np.testing.assert_array_almost_equal(sort(poles(sys1_4)), [-8., -6., -4., -2.]) - np.testing.assert_array_almost_equal(sort(zero(sys1_4)), + np.testing.assert_array_almost_equal(sort(zeros(sys1_4)), [-7., -5., -3., -1.]) sys1_5 = ctrl.series(sys1, sys2, sys3, sys4, sys5) - np.testing.assert_array_almost_equal(sort(pole(sys1_5)), + np.testing.assert_array_almost_equal(sort(poles(sys1_5)), [-8., -6., -4., -2., -0.]) - np.testing.assert_array_almost_equal(sort(zero(sys1_5)), + np.testing.assert_array_almost_equal(sort(zeros(sys1_5)), [-9., -7., -5., -3., -1.]) # Parallel sys1_2 = ctrl.parallel(sys1, sys2) - np.testing.assert_array_almost_equal(sort(pole(sys1_2)), [-4., -2.]) - np.testing.assert_array_almost_equal(sort(zero(sys1_2)), - sort(zero(sys1 + sys2))) + np.testing.assert_array_almost_equal(sort(poles(sys1_2)), [-4., -2.]) + np.testing.assert_array_almost_equal(sort(zeros(sys1_2)), + sort(zeros(sys1 + sys2))) sys1_3 = ctrl.parallel(sys1, sys2, sys3) - np.testing.assert_array_almost_equal(sort(pole(sys1_3)), + np.testing.assert_array_almost_equal(sort(poles(sys1_3)), [-6., -4., -2.]) - np.testing.assert_array_almost_equal(sort(zero(sys1_3)), - sort(zero(sys1 + sys2 + sys3))) + np.testing.assert_array_almost_equal(sort(zeros(sys1_3)), + sort(zeros(sys1 + sys2 + sys3))) sys1_4 = ctrl.parallel(sys1, sys2, sys3, sys4) - np.testing.assert_array_almost_equal(sort(pole(sys1_4)), + np.testing.assert_array_almost_equal(sort(poles(sys1_4)), [-8., -6., -4., -2.]) np.testing.assert_array_almost_equal( - sort(zero(sys1_4)), - sort(zero(sys1 + sys2 + sys3 + sys4))) + sort(zeros(sys1_4)), + sort(zeros(sys1 + sys2 + sys3 + sys4))) sys1_5 = ctrl.parallel(sys1, sys2, sys3, sys4, sys5) - np.testing.assert_array_almost_equal(sort(pole(sys1_5)), + np.testing.assert_array_almost_equal(sort(poles(sys1_5)), [-8., -6., -4., -2., -0.]) np.testing.assert_array_almost_equal( - sort(zero(sys1_5)), - sort(zero(sys1 + sys2 + sys3 + sys4 + sys5))) + sort(zeros(sys1_5)), + sort(zeros(sys1 + sys2 + sys3 + sys4 + sys5))) def testMimoSeries(self, tsys): """regression: bdalg.series reverses order of arguments""" diff --git a/control/tests/config_test.py b/control/tests/config_test.py index b495a0f6f..295c68bdd 100644 --- a/control/tests/config_test.py +++ b/control/tests/config_test.py @@ -292,8 +292,12 @@ def test_change_default_dt_static(self): """Test that static gain systems always have dt=None""" ct.set_defaults('control', default_dt=0) assert ct.tf(1, 1).dt is None - assert ct.ss(0, 0, 0, 1).dt is None - # TODO: add in test for static gain iosys + assert ct.ss([], [], [], 1).dt is None + + # Make sure static gain is preserved for the I/O system + sys = ct.ss([], [], [], 1) + sys_io = ct.ss2io(sys) + assert sys_io.dt is None def test_get_param_last(self): """Test _get_param last keyword""" diff --git a/control/tests/convert_test.py b/control/tests/convert_test.py index 36eac223c..6c4586471 100644 --- a/control/tests/convert_test.py +++ b/control/tests/convert_test.py @@ -225,7 +225,7 @@ def testTf2SsDuplicatePoles(self): [[1], [1, 0]]] g = tf(num, den) s = ss(g) - np.testing.assert_allclose(g.pole(), s.pole()) + np.testing.assert_allclose(g.poles(), s.poles()) @slycotonly def test_tf2ss_robustness(self): @@ -241,10 +241,10 @@ def test_tf2ss_robustness(self): sys2ss = tf2ss(sys2tf) # Make sure that the poles match for StateSpace and TransferFunction - np.testing.assert_array_almost_equal(np.sort(sys1tf.pole()), - np.sort(sys1ss.pole())) - np.testing.assert_array_almost_equal(np.sort(sys2tf.pole()), - np.sort(sys2ss.pole())) + np.testing.assert_array_almost_equal(np.sort(sys1tf.poles()), + np.sort(sys1ss.poles())) + np.testing.assert_array_almost_equal(np.sort(sys2tf.poles()), + np.sort(sys2ss.poles())) def test_tf2ss_nonproper(self): """Unit tests for non-proper transfer functions""" diff --git a/control/tests/flatsys_test.py b/control/tests/flatsys_test.py index 8b182a17a..a12852759 100644 --- a/control/tests/flatsys_test.py +++ b/control/tests/flatsys_test.py @@ -378,3 +378,29 @@ def test_point_to_point_errors(self): with pytest.raises(TypeError, match="unrecognized keyword"): traj_method = fs.point_to_point( flat_sys, timepts, x0, u0, xf, uf, solve_ivp_method=None) + + @pytest.mark.parametrize( + "xf, uf, Tf", + [([1, 0], [0], 2), + ([0, 1], [0], 3), + ([1, 1], [1], 4)]) + def test_response(self, xf, uf, Tf): + # Define a second order integrator + sys = ct.StateSpace([[-1, 1], [0, -2]], [[0], [1]], [[1, 0]], 0) + flatsys = fs.LinearFlatSystem(sys) + + # Define the basis set + poly = fs.PolyFamily(6) + + x1, u1, = [0, 0], [0] + traj = fs.point_to_point(flatsys, Tf, x1, u1, xf, uf, basis=poly) + + # Compute the response the regular way + T = np.linspace(0, Tf, 10) + x, u = traj.eval(T) + + # Recompute using response() + response = traj.response(T, squeeze=False) + np.testing.assert_equal(T, response.time) + np.testing.assert_equal(u, response.inputs) + np.testing.assert_equal(x, response.states) diff --git a/control/tests/frd_test.py b/control/tests/frd_test.py index af7d18bc1..ff88c3dea 100644 --- a/control/tests/frd_test.py +++ b/control/tests/frd_test.py @@ -15,6 +15,7 @@ from control.frdata import FRD, _convert_to_FRD, FrequencyResponseData from control import bdalg, evalfr, freqplot from control.tests.conftest import slycotonly +from control.exception import pandas_check class TestFRD: @@ -478,3 +479,87 @@ def test_unrecognized_keyword(self): omega = np.logspace(-1, 2, 10) with pytest.raises(TypeError, match="unrecognized keyword"): frd = FRD(h, omega, unknown=None) + + +def test_named_signals(): + ct.namedio.NamedIOSystem._idCounter = 0 + h1 = TransferFunction([1], [1, 2, 2]) + h2 = TransferFunction([1], [0.1, 1]) + omega = np.logspace(-1, 2, 10) + f1 = FRD(h1, omega) + f2 = FRD(h2, omega) + + # Make sure that systems were properly named + assert f1.name == 'sys[2]' + assert f2.name == 'sys[3]' + assert f1.ninputs == 1 + assert f1.input_labels == ['u[0]'] + assert f1.noutputs == 1 + assert f1.output_labels == ['y[0]'] + + # Change names + f1 = FRD(h1, omega, name='mysys', inputs='u0', outputs='y0') + assert f1.name == 'mysys' + assert f1.ninputs == 1 + assert f1.input_labels == ['u0'] + assert f1.noutputs == 1 + assert f1.output_labels == ['y0'] + + +@pytest.mark.skipif(not pandas_check(), reason="pandas not installed") +def test_to_pandas(): + # Create a SISO frequency response + h1 = TransferFunction([1], [1, 2, 2]) + omega = np.logspace(-1, 2, 10) + resp = FRD(h1, omega) + + # Convert to pandas + df = resp.to_pandas() + + # Check to make sure the data make senses + np.testing.assert_equal(df['omega'], resp.omega) + np.testing.assert_equal(df['H_{y[0], u[0]}'], resp.fresp[0, 0]) + + +def test_frequency_response(): + # Create an SISO frequence response + sys = ct.rss(2, 2, 2) + omega = np.logspace(-2, 2, 20) + resp = ct.frequency_response(sys, omega) + eval = sys(omega*1j) + + # Make sure we get the right answers in various ways + np.testing.assert_equal(resp.magnitude, np.abs(eval)) + np.testing.assert_equal(resp.phase, np.angle(eval)) + np.testing.assert_equal(resp.omega, omega) + + # Make sure that we can change the properties of the response + sys = ct.rss(2, 1, 1) + resp_default = ct.frequency_response(sys, omega) + mag_default, phase_default, omega_default = resp_default + assert mag_default.ndim == 1 + assert phase_default.ndim == 1 + assert omega_default.ndim == 1 + assert mag_default.shape[0] == omega_default.shape[0] + assert phase_default.shape[0] == omega_default.shape[0] + + resp_nosqueeze = ct.frequency_response(sys, omega, squeeze=False) + mag_nosqueeze, phase_nosqueeze, omega_nosqueeze = resp_nosqueeze + assert mag_nosqueeze.ndim == 3 + assert phase_nosqueeze.ndim == 3 + assert omega_nosqueeze.ndim == 1 + assert mag_nosqueeze.shape[2] == omega_nosqueeze.shape[0] + assert phase_nosqueeze.shape[2] == omega_nosqueeze.shape[0] + + # Try changing the response + resp_def_nosq = resp_default(squeeze=False) + mag_def_nosq, phase_def_nosq, omega_def_nosq = resp_def_nosq + assert mag_def_nosq.shape == mag_nosqueeze.shape + assert phase_def_nosq.shape == phase_nosqueeze.shape + assert omega_def_nosq.shape == omega_nosqueeze.shape + + resp_nosq_sq = resp_nosqueeze(squeeze=True) + mag_nosq_sq, phase_nosq_sq, omega_nosq_sq = resp_nosq_sq + assert mag_nosq_sq.shape == mag_default.shape + assert phase_nosq_sq.shape == phase_default.shape + assert omega_nosq_sq.shape == omega_default.shape diff --git a/control/tests/freqresp_test.py b/control/tests/freqresp_test.py index 4d1ac55e0..0e35a38ea 100644 --- a/control/tests/freqresp_test.py +++ b/control/tests/freqresp_test.py @@ -418,11 +418,11 @@ def test_dcgain_consistency(): """Test to make sure that DC gain is consistently evaluated""" # Set up transfer function with pole at the origin sys_tf = ctrl.tf([1], [1, 0]) - assert 0 in sys_tf.pole() + assert 0 in sys_tf.poles() # Set up state space system with pole at the origin sys_ss = ctrl.tf2ss(sys_tf) - assert 0 in sys_ss.pole() + assert 0 in sys_ss.poles() # Finite (real) numerator over 0 denominator => inf + nanj np.testing.assert_equal( @@ -440,8 +440,8 @@ def test_dcgain_consistency(): # Set up transfer function with pole, zero at the origin sys_tf = ctrl.tf([1, 0], [1, 0]) - assert 0 in sys_tf.pole() - assert 0 in sys_tf.zero() + assert 0 in sys_tf.poles() + assert 0 in sys_tf.zeros() # Pole and zero at the origin should give nan + nanj for the response np.testing.assert_equal( @@ -456,7 +456,7 @@ def test_dcgain_consistency(): ctrl.tf2ss(ctrl.tf([1], [1, 0])) # Different systems give different representations => test accordingly - if 0 in sys_ss.pole() and 0 in sys_ss.zero(): + if 0 in sys_ss.poles() and 0 in sys_ss.zeros(): # Pole and zero at the origin => should get (nan + nanj) np.testing.assert_equal( sys_ss(0, warn_infinite=False), complex(np.nan, np.nan)) @@ -464,7 +464,7 @@ def test_dcgain_consistency(): sys_ss(0j, warn_infinite=False), complex(np.nan, np.nan)) np.testing.assert_equal( sys_ss.dcgain(), np.nan) - elif 0 in sys_ss.pole(): + elif 0 in sys_ss.poles(): # Pole at the origin, but zero elsewhere => should get (inf + nanj) np.testing.assert_equal( sys_ss(0, warn_infinite=False), complex(np.inf, np.nan)) @@ -479,11 +479,11 @@ def test_dcgain_consistency(): # Pole with non-zero, complex numerator => inf + infj s = ctrl.tf('s') sys_tf = (s + 1) / (s**2 + 1) - assert 1j in sys_tf.pole() + assert 1j in sys_tf.poles() # Set up state space system with pole on imaginary axis sys_ss = ctrl.tf2ss(sys_tf) - assert 1j in sys_tf.pole() + assert 1j in sys_tf.poles() # Make sure we get correct response if evaluated at the pole np.testing.assert_equal( diff --git a/control/tests/iosys_test.py b/control/tests/iosys_test.py index 4e0adfa03..ecb30c316 100644 --- a/control/tests/iosys_test.py +++ b/control/tests/iosys_test.py @@ -130,7 +130,8 @@ def test_iosys_print(self, tsys, capsys): print(ios_linearized) @noscipy0 - def test_nonlinear_iosys(self, tsys): + @pytest.mark.parametrize("ss", [ios.NonlinearIOSystem, ct.ss]) + def test_nonlinear_iosys(self, tsys, ss): # Create a simple nonlinear I/O system nlsys = ios.NonlinearIOSystem(predprey) T = tsys.T @@ -239,9 +240,9 @@ def test_linearize_named_signals(self, kincar): def test_connect(self, tsys): # Define a couple of (linear) systems to interconnection linsys1 = tsys.siso_linsys - iosys1 = ios.LinearIOSystem(linsys1) + iosys1 = ios.LinearIOSystem(linsys1, name='iosys1') linsys2 = tsys.siso_linsys - iosys2 = ios.LinearIOSystem(linsys2) + iosys2 = ios.LinearIOSystem(linsys2, name='iosys2') # Connect systems in different ways and compare to StateSpace linsys_series = linsys2 * linsys1 @@ -407,8 +408,8 @@ def test_algebraic_loop(self, tsys): lnios = ios.LinearIOSystem(linsys) nlios = ios.NonlinearIOSystem(None, \ lambda t, x, u, params: u*u, inputs=1, outputs=1) - nlios1 = nlios.copy() - nlios2 = nlios.copy() + nlios1 = nlios.copy(name='nlios1') + nlios2 = nlios.copy(name='nlios2') # Set up parameters for simulation T, U, X0 = tsys.T, tsys.U, tsys.X0 @@ -473,8 +474,8 @@ def test_algebraic_loop(self, tsys): def test_summer(self, tsys): # Construct a MIMO system for testing linsys = tsys.mimo_linsys1 - linio1 = ios.LinearIOSystem(linsys) - linio2 = ios.LinearIOSystem(linsys) + linio1 = ios.LinearIOSystem(linsys, name='linio1') + linio2 = ios.LinearIOSystem(linsys, name='linio2') linsys_parallel = linsys + linsys iosys_parallel = linio1 + linio2 @@ -1042,8 +1043,12 @@ def test_sys_naming_convention(self, tsys): ct.config.use_legacy_defaults('0.8.4') # changed delims in 0.9.0 ct.config.use_numpy_matrix(False) # np.matrix deprecated - ct.namedio._NamedIOSystem._idCounter = 0 - sys = ct.LinearIOSystem(tsys.mimo_linsys1) + + # Create a system with a known ID + ct.namedio.NamedIOSystem._idCounter = 0 + sys = ct.ss( + tsys.mimo_linsys1.A, tsys.mimo_linsys1.B, + tsys.mimo_linsys1.C, tsys.mimo_linsys1.D) assert sys.name == "sys[0]" assert sys.copy().name == "copy of sys[0]" @@ -1093,7 +1098,7 @@ def test_sys_naming_convention(self, tsys): # Same system conflict with pytest.warns(UserWarning): - unnamedsys1 * unnamedsys1 + namedsys * namedsys @pytest.mark.usefixtures("editsdefaults") def test_signals_naming_convention_0_8_4(self, tsys): @@ -1106,8 +1111,13 @@ def test_signals_naming_convention_0_8_4(self, tsys): ct.config.use_legacy_defaults('0.8.4') # changed delims in 0.9.0 ct.config.use_numpy_matrix(False) # np.matrix deprecated - ct.namedio._NamedIOSystem._idCounter = 0 - sys = ct.LinearIOSystem(tsys.mimo_linsys1) + + # Create a system with a known ID + ct.namedio.NamedIOSystem._idCounter = 0 + sys = ct.ss( + tsys.mimo_linsys1.A, tsys.mimo_linsys1.B, + tsys.mimo_linsys1.C, tsys.mimo_linsys1.D) + for statename in ["x[0]", "x[1]"]: assert statename in sys.state_index for inputname in ["u[0]", "u[1]"]: @@ -1154,9 +1164,9 @@ def test_signals_naming_convention_0_8_4(self, tsys): # Same system conflict with pytest.warns(UserWarning): - same_name_series = unnamedsys * unnamedsys - assert "sys[1].x[0]" in same_name_series.state_index - assert "copy of sys[1].x[0]" in same_name_series.state_index + same_name_series = namedsys * namedsys + assert "namedsys.x0" in same_name_series.state_index + assert "copy of namedsys.x0" in same_name_series.state_index def test_named_signals_linearize_inconsistent(self, tsys): """Mare sure that providing inputs or outputs not consistent with @@ -1206,13 +1216,13 @@ def outfcn(t, x, u, params): def test_lineariosys_statespace(self, tsys): """Make sure that a LinearIOSystem is also a StateSpace object""" - iosys_siso = ct.LinearIOSystem(tsys.siso_linsys) - iosys_siso2 = ct.LinearIOSystem(tsys.siso_linsys) + iosys_siso = ct.LinearIOSystem(tsys.siso_linsys, name='siso') + iosys_siso2 = ct.LinearIOSystem(tsys.siso_linsys, name='siso2') assert isinstance(iosys_siso, ct.StateSpace) # Make sure that state space functions work for LinearIOSystems np.testing.assert_allclose( - iosys_siso.pole(), tsys.siso_linsys.pole()) + iosys_siso.poles(), tsys.siso_linsys.poles()) omega = np.logspace(.1, 10, 100) mag_io, phase_io, omega_io = iosys_siso.frequency_response(omega) mag_ss, phase_ss, omega_ss = tsys.siso_linsys.frequency_response(omega) @@ -1390,7 +1400,7 @@ def test_duplicates(self, tsys): name="sys") # Duplicate objects - with pytest.warns(UserWarning, match="Duplicate object"): + with pytest.warns(UserWarning, match="duplicate object"): ios_series = nlios * nlios # Nonduplicate objects @@ -1398,7 +1408,7 @@ def test_duplicates(self, tsys): ct.config.use_numpy_matrix(False) # np.matrix deprecated nlios1 = nlios.copy() nlios2 = nlios.copy() - with pytest.warns(UserWarning, match="Duplicate name"): + with pytest.warns(UserWarning, match="duplicate name"): ios_series = nlios1 * nlios2 assert "copy of sys_1.x[0]" in ios_series.state_index.keys() assert "copy of sys.x[0]" in ios_series.state_index.keys() @@ -1412,7 +1422,7 @@ def test_duplicates(self, tsys): lambda t, x, u, params: u * u, inputs=1, outputs=1, name="sys") - with pytest.warns(UserWarning, match="Duplicate name"): + with pytest.warns(UserWarning, match="duplicate name"): ct.InterconnectedSystem([nlios1, iosys_siso, nlios2], inputs=0, outputs=0, states=0) @@ -1711,6 +1721,7 @@ def test_interconnect_unused_output(): outputs=['y'], ignore_outputs=['v']) + def test_input_output_broadcasting(): # Create a system, time vector, and noisy input sys = ct.rss(6, 2, 3) @@ -1752,6 +1763,7 @@ def test_input_output_broadcasting(): resp_bad = ct.input_output_response( sys, T, (U[0, :], U[:2, :-1]), [X0, P0]) + def test_nonuniform_timepts(): """Test non-uniform time points for simulations""" sys = ct.LinearIOSystem(ct.rss(2, 1, 1)) @@ -1773,3 +1785,77 @@ def test_nonuniform_timepts(): t_even, y_even = ct.input_output_response( sys, noufpts, nonunif, t_eval=unifpts) np.testing.assert_almost_equal(y_unif, y_even, decimal=6) + + +def test_ss_nonlinear(): + """Test ss() for creating nonlinear systems""" + secord = ct.ss(secord_update, secord_output, inputs='u', outputs='y', + states = ['x1', 'x2'], name='secord') + assert secord.name == 'secord' + assert secord.input_labels == ['u'] + assert secord.output_labels == ['y'] + assert secord.state_labels == ['x1', 'x2'] + + # Make sure we get the same answer for simulations + T = np.linspace(0, 10, 100) + U = np.sin(T) + X0 = np.array([1, -1]) + secord_nlio = ct.NonlinearIOSystem( + secord_update, secord_output, inputs=1, outputs=1, states=2) + ss_response = ct.input_output_response(secord, T, U, X0) + io_response = ct.input_output_response(secord_nlio, T, U, X0) + np.testing.assert_almost_equal(ss_response.time, io_response.time) + np.testing.assert_almost_equal(ss_response.inputs, io_response.inputs) + np.testing.assert_almost_equal(ss_response.outputs, io_response.outputs) + + # Make sure that optional keywords are allowed + secord = ct.ss(secord_update, secord_output, dt=True) + assert ct.isdtime(secord) + + # Make sure that state space keywords are flagged + with pytest.raises(TypeError, match="unrecognized keyword"): + ct.ss(secord_update, remove_useless_states=True) + + +def test_rss(): + # Basic call, with no arguments + sys = ct.rss() + assert sys.ninputs == 1 + assert sys.noutputs == 1 + assert sys.nstates == 1 + assert sys.dt == 0 + assert np.all(np.real(sys.poles()) < 0) + + # Set the timebase explicitly + sys = ct.rss(inputs=2, outputs=3, states=4, dt=None, name='sys') + assert sys.name == 'sys' + assert sys.ninputs == 2 + assert sys.noutputs == 3 + assert sys.nstates == 4 + assert sys.dt == None + assert np.all(np.real(sys.poles()) < 0) + + # Discrete time + sys = ct.rss(inputs=['a', 'b'], outputs=1, states=1, dt=True) + assert sys.ninputs == 2 + assert sys.input_labels == ['a', 'b'] + assert sys.noutputs == 1 + assert sys.nstates == 1 + assert sys.dt == True + assert np.all(np.abs(sys.poles()) < 1) + + # Call drss directly + sys = ct.drss(inputs=['a', 'b'], outputs=1, states=1, dt=True) + assert sys.ninputs == 2 + assert sys.input_labels == ['a', 'b'] + assert sys.noutputs == 1 + assert sys.nstates == 1 + assert sys.dt == True + assert np.all(np.abs(sys.poles()) < 1) + + with pytest.raises(ValueError, match="continuous timebase"): + sys = ct.drss(2, 1, 1, dt=0) + + with pytest.warns(UserWarning, match="may be interpreted as continuous"): + sys = ct.drss(2, 1, 1, dt=None) + assert np.all(np.abs(sys.poles()) < 1) diff --git a/control/tests/kwargs_test.py b/control/tests/kwargs_test.py index 818a906a5..ada16a46a 100644 --- a/control/tests/kwargs_test.py +++ b/control/tests/kwargs_test.py @@ -99,7 +99,9 @@ def test_unrecognized_kwargs(): [control.summing_junction, (2,), {}], [control.tf, ([1], [1, 1]), {}], [control.tf2io, (control.tf([1], [1, 1]),), {}], - [control.InputOutputSystem, (1, 1, 1), {}], + [control.tf2ss, (control.tf([1], [1, 1]),), {}], + [control.InputOutputSystem, (), + {'inputs': 1, 'outputs': 1, 'states': 1}], [control.InputOutputSystem.linearize, (sys, 0, 0), {}], [control.StateSpace, ([[-1, 0], [0, -1]], [[1], [1]], [[1, 1]], 0), {}], [control.TransferFunction, ([1], [1, 1]), {}], @@ -113,18 +115,23 @@ def test_unrecognized_kwargs(): with pytest.raises(TypeError, match="unrecognized keyword"): function(*args, **kwargs, unknown=None) + # If we opened any figures, close them to avoid matplotlib warnings + if plt.gca(): + plt.close('all') + def test_matplotlib_kwargs(): # Create a SISO system for use in parameterized tests sys = control.ss([[-1, 1], [0, -1]], [[0], [1]], [[1, 0]], 0, dt=None) + ctl = control.ss([[-1, 1], [0, -1]], [[0], [1]], [[1, 0]], 0, dt=None) table = [ [control.bode, (sys, ), {}], [control.bode_plot, (sys, ), {}], [control.describing_function_plot, (sys, control.descfcn.saturation_nonlinearity(1), [1, 2, 3, 4]), {}], - [control.gangof4, (sys, sys), {}], - [control.gangof4_plot, (sys, sys), {}], + [control.gangof4, (sys, ctl), {}], + [control.gangof4_plot, (sys, ctl), {}], [control.nyquist, (sys, ), {}], [control.nyquist_plot, (sys, ), {}], [control.singular_values_plot, (sys, ), {}], @@ -138,7 +145,7 @@ def test_matplotlib_kwargs(): with pytest.raises(AttributeError, match="has no property"): function(*args, **kwargs, unknown=None) - # If we opened any figures, close them + # If we opened any figures, close them to avoid matplotlib warnings if plt.gca(): plt.close('all') @@ -157,7 +164,7 @@ def test_matplotlib_kwargs(): 'bode_plot': test_matplotlib_kwargs, 'describing_function_plot': test_matplotlib_kwargs, 'dlqe': test_unrecognized_kwargs, - 'dlqr': statefbk_test.TestStatefbk.test_lqr_errors, + 'dlqr': test_unrecognized_kwargs, 'drss': test_unrecognized_kwargs, 'gangof4': test_matplotlib_kwargs, 'gangof4_plot': test_matplotlib_kwargs, @@ -165,10 +172,10 @@ def test_matplotlib_kwargs(): 'interconnect': interconnect_test.test_interconnect_exceptions, 'linearize': test_unrecognized_kwargs, 'lqe': test_unrecognized_kwargs, - 'lqr': statefbk_test.TestStatefbk.test_lqr_errors, + 'lqr': test_unrecognized_kwargs, 'nyquist': test_matplotlib_kwargs, 'nyquist_plot': test_matplotlib_kwargs, - 'pzmap': test_matplotlib_kwargs, + 'pzmap': test_unrecognized_kwargs, 'rlocus': test_unrecognized_kwargs, 'root_locus': test_unrecognized_kwargs, 'rss': test_unrecognized_kwargs, @@ -180,6 +187,7 @@ def test_matplotlib_kwargs(): 'summing_junction': interconnect_test.test_interconnect_exceptions, 'tf': test_unrecognized_kwargs, 'tf2io' : test_unrecognized_kwargs, + 'tf2ss' : test_unrecognized_kwargs, 'flatsys.point_to_point': flatsys_test.TestFlatSys.test_point_to_point_errors, 'FrequencyResponseData.__init__': diff --git a/control/tests/lti_test.py b/control/tests/lti_test.py index e2f7f2e03..8e45ea482 100644 --- a/control/tests/lti_test.py +++ b/control/tests/lti_test.py @@ -5,23 +5,44 @@ from .conftest import editsdefaults import control as ct -from control import c2d, tf, tf2ss, NonlinearIOSystem -from control.lti import (LTI, common_timebase, evalfr, damp, dcgain, isctime, - isdtime, issiso, pole, timebaseEqual, zero) +from control import c2d, tf, ss, tf2ss, NonlinearIOSystem +from control.lti import LTI, evalfr, damp, dcgain, zeros, poles +from control import common_timebase, isctime, isdtime, issiso, timebaseEqual from control.tests.conftest import slycotonly from control.exception import slycot_check class TestLTI: + @pytest.mark.parametrize("fun, args", [ + [tf, (126, [-1, 42])], + [ss, ([[42]], [[1]], [[1]], 0)] + ]) + def test_poles(self, fun, args): + sys = fun(*args) + np.testing.assert_allclose(sys.poles(), 42) + np.testing.assert_allclose(poles(sys), 42) + + with pytest.warns(PendingDeprecationWarning): + pole_list = sys.pole() + assert pole_list == sys.poles() + + with pytest.warns(PendingDeprecationWarning): + pole_list = ct.pole(sys) + assert pole_list == sys.poles() + + @pytest.mark.parametrize("fun, args", [ + [tf, (126, [-1, 42])], + [ss, ([[42]], [[1]], [[1]], 0)] + ]) + def test_zero(self, fun, args): + sys = fun(*args) + np.testing.assert_allclose(sys.zeros(), 42) + np.testing.assert_allclose(zeros(sys), 42) - def test_pole(self): - sys = tf(126, [-1, 42]) - np.testing.assert_allclose(sys.pole(), 42) - np.testing.assert_allclose(pole(sys), 42) + with pytest.warns(PendingDeprecationWarning): + sys.zero() - def test_zero(self): - sys = tf([-1, 42], [1, 10]) - np.testing.assert_allclose(sys.zero(), 42) - np.testing.assert_allclose(zero(sys), 42) + with pytest.warns(PendingDeprecationWarning): + ct.zero(sys) def test_issiso(self): assert issiso(1) @@ -267,7 +288,7 @@ def test_squeeze_exceptions(self, fcn): sys = fcn(ct.rss(2, 1, 1)) with pytest.raises(ValueError, match="unknown squeeze value"): - sys.frequency_response([1], squeeze=1) + resp = sys.frequency_response([1], squeeze='siso') with pytest.raises(ValueError, match="unknown squeeze value"): sys([1j], squeeze='siso') with pytest.raises(ValueError, match="unknown squeeze value"): diff --git a/control/tests/minreal_test.py b/control/tests/minreal_test.py index 466f9384d..10c56d4ca 100644 --- a/control/tests/minreal_test.py +++ b/control/tests/minreal_test.py @@ -7,7 +7,7 @@ from scipy.linalg import eigvals import pytest -from control import rss, ss, zero +from control import rss, ss, zeros from control.statesp import StateSpace from control.xferfcn import TransferFunction from itertools import permutations @@ -64,8 +64,8 @@ def testMinrealBrute(self): # Check that the zeros match # Note: sorting doesn't work => have to do the hard way - z1 = zero(s1) - z2 = zero(s2) + z1 = zeros(s1) + z2 = zeros(s2) # Start by making sure we have the same # of zeros assert len(z1) == len(z2) diff --git a/control/tests/namedio_test.py b/control/tests/namedio_test.py index 9278136b5..3a96203a8 100644 --- a/control/tests/namedio_test.py +++ b/control/tests/namedio_test.py @@ -9,11 +9,13 @@ """ import re +from copy import copy import numpy as np import control as ct import pytest + def test_named_ss(): # Create a system to play with sys = ct.rss(2, 2, 2) @@ -23,29 +25,229 @@ def test_named_ss(): # Get the state matrices for later use A, B, C, D = sys.A, sys.B, sys.C, sys.D - + # Set up a named state space systems with default names - ct.namedio._NamedIOSystem._idCounter = 0 + ct.namedio.NamedIOSystem._idCounter = 0 sys = ct.ss(A, B, C, D) assert sys.name == 'sys[0]' assert sys.input_labels == ['u[0]', 'u[1]'] assert sys.output_labels == ['y[0]', 'y[1]'] assert sys.state_labels == ['x[0]', 'x[1]'] + assert repr(sys) == \ + "['y[0]', 'y[1]']>" # Pass the names as arguments sys = ct.ss( A, B, C, D, name='system', inputs=['u1', 'u2'], outputs=['y1', 'y2'], states=['x1', 'x2']) assert sys.name == 'system' - assert ct.namedio._NamedIOSystem._idCounter == 1 + assert ct.namedio.NamedIOSystem._idCounter == 1 assert sys.input_labels == ['u1', 'u2'] assert sys.output_labels == ['y1', 'y2'] assert sys.state_labels == ['x1', 'x2'] + assert repr(sys) == \ + "['y1', 'y2']>" # Do the same with rss sys = ct.rss(['x1', 'x2', 'x3'], ['y1', 'y2'], 'u1', name='random') assert sys.name == 'random' - assert ct.namedio._NamedIOSystem._idCounter == 1 + assert ct.namedio.NamedIOSystem._idCounter == 1 assert sys.input_labels == ['u1'] assert sys.output_labels == ['y1', 'y2'] assert sys.state_labels == ['x1', 'x2', 'x3'] + assert repr(sys) == \ + "['y1', 'y2']>" + + +# List of classes that are expected +fun_instance = { + ct.rss: (ct.InputOutputSystem, ct.LinearIOSystem, ct.StateSpace), + ct.drss: (ct.InputOutputSystem, ct.LinearIOSystem, ct.StateSpace), + ct.FRD: (ct.lti.LTI), + ct.NonlinearIOSystem: (ct.InputOutputSystem), + ct.ss: (ct.InputOutputSystem, ct.LinearIOSystem, ct.StateSpace), + ct.StateSpace: (ct.StateSpace), + ct.tf: (ct.TransferFunction), + ct.TransferFunction: (ct.TransferFunction), +} + +# List of classes that are not expected +fun_notinstance = { + ct.FRD: (ct.InputOutputSystem, ct.LinearIOSystem, ct.StateSpace), + ct.StateSpace: (ct.InputOutputSystem, ct.TransferFunction), + ct.TransferFunction: (ct.InputOutputSystem, ct.StateSpace), +} + + +@pytest.mark.parametrize("fun, args, kwargs", [ + [ct.rss, (4, 1, 1), {}], + [ct.rss, (3, 2, 1), {}], + [ct.drss, (4, 1, 1), {}], + [ct.drss, (3, 2, 1), {}], + [ct.FRD, ([1, 2, 3,], [1, 2, 3]), {}], + [ct.NonlinearIOSystem, + (lambda t, x, u, params: -x, None), + {'inputs': 2, 'outputs':2, 'states':2}], + [ct.ss, ([[1, 2], [3, 4]], [[0], [1]], [[1, 0]], 0), {}], + [ct.StateSpace, ([[1, 2], [3, 4]], [[0], [1]], [[1, 0]], 0), {}], + [ct.tf, ([1, 2], [3, 4, 5]), {}], + [ct.TransferFunction, ([1, 2], [3, 4, 5]), {}], +]) +def test_io_naming(fun, args, kwargs): + # Reset the ID counter to get uniform generic names + ct.namedio.NamedIOSystem._idCounter = 0 + + # Create the system w/out any names + sys_g = fun(*args, **kwargs) + + # Make sure the class are what we expect + if fun in fun_instance: + assert isinstance(sys_g, fun_instance[fun]) + + if fun in fun_notinstance: + assert not isinstance(sys_g, fun_notinstance[fun]) + + # Make sure the names make sense + assert sys_g.name == 'sys[0]' + assert sys_g.input_labels == [f'u[{i}]' for i in range(sys_g.ninputs)] + assert sys_g.output_labels == [f'y[{i}]' for i in range(sys_g.noutputs)] + if sys_g.nstates: + assert sys_g.state_labels == [f'x[{i}]' for i in range(sys_g.nstates)] + + # + # Reset the names to something else and make sure they stick + # + sys_r = copy(sys_g) + + input_labels = [f'u{i}' for i in range(sys_g.ninputs)] + sys_r.set_inputs(input_labels) + assert sys_r.input_labels == input_labels + + output_labels = [f'y{i}' for i in range(sys_g.noutputs)] + sys_r.set_outputs(output_labels) + assert sys_r.output_labels == output_labels + + if sys_g.nstates: + state_labels = [f'x{i}' for i in range(sys_g.nstates)] + sys_r.set_states(state_labels) + assert sys_r.state_labels == state_labels + + # + # Set names using keywords and make sure they stick + # + + # How the keywords are used depends on the type of system + if fun in (ct.rss, ct.drss): + # Pass the labels instead of the numbers + sys_k = fun(state_labels, output_labels, input_labels, name='mysys') + + elif sys_g.nstates is None: + # Don't pass state labels + sys_k = fun( + *args, inputs=input_labels, outputs=output_labels, name='mysys') + + else: + sys_k = fun( + *args, inputs=input_labels, outputs=output_labels, + states=state_labels, name='mysys') + + assert sys_k.name == 'mysys' + assert sys_k.input_labels == input_labels + assert sys_k.output_labels == output_labels + if sys_g.nstates: + assert sys_k.state_labels == state_labels + + # + # Convert the system to state space and make sure labels transfer + # + if ct.slycot_check() and not isinstance( + sys_r, (ct.FrequencyResponseData, ct.NonlinearIOSystem)): + sys_ss = ct.ss(sys_r) + assert sys_ss != sys_r + assert sys_ss.input_labels == input_labels + assert sys_ss.output_labels == output_labels + + # Reassign system and signal names + sys_ss = ct.ss( + sys_g, inputs=input_labels, outputs=output_labels, name='new') + assert sys_ss.name == 'new' + assert sys_ss.input_labels == input_labels + assert sys_ss.output_labels == output_labels + + # + # Convert the system to a transfer function and make sure labels transfer + # + if not isinstance( + sys_r, (ct.FrequencyResponseData, ct.NonlinearIOSystem)) and \ + ct.slycot_check(): + sys_tf = ct.tf(sys_r) + assert sys_tf != sys_r + assert sys_tf.input_labels == input_labels + assert sys_tf.output_labels == output_labels + + # Reassign system and signal names + sys_tf = ct.tf( + sys_g, inputs=input_labels, outputs=output_labels, name='new') + assert sys_tf.name == 'new' + assert sys_tf.input_labels == input_labels + assert sys_tf.output_labels == output_labels + + +# Internal testing of StateSpace initialization +def test_init_namedif(): + # Set up the initial system + sys = ct.rss(2, 1, 1) + + # Rename the system, inputs, and outouts + sys_new = sys.copy() + ct.StateSpace.__init__( + sys_new, sys, inputs='u', outputs='y', name='new') + assert sys_new.name == 'new' + assert sys_new.input_labels == ['u'] + assert sys_new.output_labels == ['y'] + + # Call constructor without re-initialization + sys_keep = sys.copy() + ct.StateSpace.__init__(sys_keep, sys, init_namedio=False) + assert sys_keep.name == sys_keep.name + assert sys_keep.input_labels == sys_keep.input_labels + assert sys_keep.output_labels == sys_keep.output_labels + + # Make sure that passing an unrecognized keyword generates an error + with pytest.raises(TypeError, match="unrecognized keyword"): + ct.StateSpace.__init__( + sys_keep, sys, inputs='u', outputs='y', init_namedio=False) + +# Test state space conversion +def test_convert_to_statespace(): + # Set up the initial system + sys = ct.tf(ct.rss(2, 1, 1)) + + # Make sure we can rename system name, inputs, outputs + sys_new = ct.ss(sys, inputs='u', outputs='y', name='new') + assert sys_new.name == 'new' + assert sys_new.input_labels == ['u'] + assert sys_new.output_labels == ['y'] + + # Try specifying the state names (via low level test) + with pytest.warns(UserWarning, match="non-unique state space realization"): + sys_new = ct.ss(sys, inputs='u', outputs='y', states=['x1', 'x2']) + assert sys_new.input_labels == ['u'] + assert sys_new.output_labels == ['y'] + assert sys_new.state_labels == ['x1', 'x2'] + + +# Duplicate name warnings +def test_duplicate_sysname(): + # Start with an unnamed system + sys = ct.rss(4, 1, 1) + + # No warnings should be generated if we reuse an an unnamed system + with pytest.warns(None) as record: + res = sys * sys + assert not any([type(msg) == UserWarning for msg in record]) + + # Generate a warning if the system is named + sys = ct.rss(4, 1, 1, name='sys') + with pytest.warns(UserWarning, match="duplicate object found"): + res = sys * sys diff --git a/control/tests/nyquist_test.py b/control/tests/nyquist_test.py index c77d94c86..a001598a6 100644 --- a/control/tests/nyquist_test.py +++ b/control/tests/nyquist_test.py @@ -19,11 +19,11 @@ # Utility function for counting unstable poles of open loop (P in FBS) def _P(sys, indent='right'): if indent == 'right': - return (sys.pole().real > 0).sum() + return (sys.poles().real > 0).sum() elif indent == 'left': - return (sys.pole().real >= 0).sum() + return (sys.poles().real >= 0).sum() elif indent == 'none': - if any(sys.pole().real == 0): + if any(sys.poles().real == 0): raise ValueError("indent must be left or right for imaginary pole") else: raise TypeError("unknown indent value") @@ -31,7 +31,7 @@ def _P(sys, indent='right'): # Utility function for counting unstable poles of closed loop (Z in FBS) def _Z(sys): - return (sys.feedback().pole().real >= 0).sum() + return (sys.feedback().poles().real >= 0).sum() # Basic tests @@ -308,6 +308,6 @@ def test_nyquist_exceptions(): print("Unusual Nyquist plot") sys = ct.tf([1], [1, 3, 2]) * ct.tf([1], [1, 0, 1]) plt.figure() - plt.title("Poles: %s" % np.array2string(sys.pole(), precision=2, separator=',')) + plt.title("Poles: %s" % np.array2string(sys.poles(), precision=2, separator=',')) count = ct.nyquist_plot(sys) assert _Z(sys) == count + _P(sys) diff --git a/control/tests/rlocus_test.py b/control/tests/rlocus_test.py index ef9bd7ecb..a0ecebb15 100644 --- a/control/tests/rlocus_test.py +++ b/control/tests/rlocus_test.py @@ -41,7 +41,7 @@ def sys(self, request): def check_cl_poles(self, sys, pole_list, k_list): for k, poles in zip(k_list, pole_list): - poles_expected = np.sort(feedback(sys, k).pole()) + poles_expected = np.sort(feedback(sys, k).poles()) poles = np.sort(poles) np.testing.assert_array_almost_equal(poles, poles_expected) diff --git a/control/tests/statefbk_test.py b/control/tests/statefbk_test.py index 9f04b3723..13f164e1f 100644 --- a/control/tests/statefbk_test.py +++ b/control/tests/statefbk_test.py @@ -7,7 +7,7 @@ import pytest import control as ct -from control import lqe, dlqe, pole, rss, ss, tf +from control import lqe, dlqe, poles, rss, ss, tf from control.exception import ControlDimension, ControlSlycot, \ ControlArgument, slycot_check from control.mateqn import care, dare @@ -167,12 +167,12 @@ def testAcker(self, fixedseed): # Place the poles at random locations des = rss(states, 1, 1) - poles = pole(des) + desired = poles(des) # Now place the poles using acker - K = acker(sys.A, sys.B, poles) + K = acker(sys.A, sys.B, desired) new = ss(sys.A - sys.B * K, sys.B, sys.C, sys.D) - placed = pole(new) + placed = poles(new) # Debugging code # diff = np.sort(poles) - np.sort(placed) @@ -181,8 +181,8 @@ def testAcker(self, fixedseed): # print(sys) # print("desired = ", poles) - np.testing.assert_array_almost_equal(np.sort(poles), - np.sort(placed), decimal=4) + np.testing.assert_array_almost_equal( + np.sort(desired), np.sort(placed), decimal=4) def checkPlaced(self, P_expected, P_placed): """Check that placed poles are correct""" @@ -679,7 +679,7 @@ def test_lqr_integral_continuous(self): np.testing.assert_array_almost_equal(clsys.D, D_clsys) # Check the poles of the closed loop system - assert all(np.real(clsys.pole()) < 0) + assert all(np.real(clsys.poles()) < 0) # Make sure controller infinite zero frequency gain if slycot_check(): diff --git a/control/tests/statesp_test.py b/control/tests/statesp_test.py index be6cd9a6b..f7757f2e9 100644 --- a/control/tests/statesp_test.py +++ b/control/tests/statesp_test.py @@ -125,7 +125,7 @@ def test_constructor(self, sys322ABCD, dt, argfun): @pytest.mark.parametrize("args, exc, errmsg", [((True, ), TypeError, "(can only take in|sys must be) a StateSpace"), - ((1, 2), ValueError, "1, 4, or 5 arguments"), + ((1, 2), TypeError, "1, 4, or 5 arguments"), ((np.ones((3, 2)), np.ones((3, 2)), np.ones((2, 2)), np.ones((2, 2))), ValueError, "A must be square"), @@ -180,16 +180,17 @@ def test_copy_constructor(self): linsys.A[0, 0] = -3 np.testing.assert_allclose(cpysys.A, [[-1]]) # original value + @pytest.mark.skip("obsolete test") def test_copy_constructor_nodt(self, sys322): """Test the copy constructor when an object without dt is passed""" sysin = sample_system(sys322, 1.) - del sysin.dt + del sysin.dt # this is a nonsensical thing to do sys = StateSpace(sysin) assert sys.dt == defaults['control.default_dt'] # test for static gain sysin = StateSpace([], [], [], [[1, 2], [3, 4]], 1.) - del sysin.dt + del sysin.dt # this is a nonsensical thing to do sys = StateSpace(sysin) assert sys.dt is None @@ -230,7 +231,7 @@ def test_D_broadcast(self, sys623): def test_pole(self, sys322): """Evaluate the poles of a MIMO system.""" - p = np.sort(sys322.pole()) + p = np.sort(sys322.poles()) true_p = np.sort([3.34747678408874, -3.17373839204437 + 1.47492908003839j, -3.17373839204437 - 1.47492908003839j]) @@ -240,7 +241,7 @@ def test_pole(self, sys322): def test_zero_empty(self): """Test to make sure zero() works with no zeros in system.""" sys = _convert_to_statespace(TransferFunction([1], [1, 2, 1])) - np.testing.assert_array_equal(sys.zero(), np.array([])) + np.testing.assert_array_equal(sys.zeros(), np.array([])) @slycotonly def test_zero_siso(self, sys222): @@ -252,9 +253,9 @@ def test_zero_siso(self, sys222): # compute zeros as root of the characteristic polynomial at the numerator of tf111 # this method is simple and assumed as valid in this test - true_z = np.sort(tf111[0, 0].zero()) + true_z = np.sort(tf111[0, 0].zeros()) # Compute the zeros through ab08nd, which is tested here - z = np.sort(sys111.zero()) + z = np.sort(sys111.zeros()) np.testing.assert_almost_equal(true_z, z) @@ -262,7 +263,7 @@ def test_zero_siso(self, sys222): def test_zero_mimo_sys322_square(self, sys322): """Evaluate the zeros of a square MIMO system.""" - z = np.sort(sys322.zero()) + z = np.sort(sys322.zeros()) true_z = np.sort([44.41465, -0.490252, -5.924398]) np.testing.assert_array_almost_equal(z, true_z) @@ -270,7 +271,7 @@ def test_zero_mimo_sys322_square(self, sys322): def test_zero_mimo_sys222_square(self, sys222): """Evaluate the zeros of a square MIMO system.""" - z = np.sort(sys222.zero()) + z = np.sort(sys222.zeros()) true_z = np.sort([-10.568501, 3.368501]) np.testing.assert_array_almost_equal(z, true_z) @@ -278,7 +279,7 @@ def test_zero_mimo_sys222_square(self, sys222): def test_zero_mimo_sys623_non_square(self, sys623): """Evaluate the zeros of a non square MIMO system.""" - z = np.sort(sys623.zero()) + z = np.sort(sys623.zeros()) true_z = np.sort([2., -1.]) np.testing.assert_array_almost_equal(z, true_z) @@ -570,15 +571,24 @@ def test_scalar_static_gain(self): """ g1 = StateSpace([], [], [], [2]) g2 = StateSpace([], [], [], [3]) + assert g1.dt == None + assert g2.dt == None g3 = g1 * g2 assert 6 == g3.D[0, 0] + assert g3.dt == None + g4 = g1 + g2 assert 5 == g4.D[0, 0] + assert g4.dt == None + g5 = g1.feedback(g2) np.testing.assert_allclose(2. / 7, g5.D[0, 0]) + assert g5.dt == None + g6 = g1.append(g2) np.testing.assert_allclose(np.diag([2, 3]), g6.D) + assert g6.dt == None def test_matrix_static_gain(self): """Regression: can we create matrix static gains?""" @@ -749,9 +759,9 @@ def test_str(self, sys322): assert str(sysdt1) == tref + "\ndt = {}\n".format(1.) def test_pole_static(self): - """Regression: pole() of static gain is empty array.""" + """Regression: poles() of static gain is empty array.""" np.testing.assert_array_equal(np.array([]), - StateSpace([], [], [], [[1]]).pole()) + StateSpace([], [], [], [[1]]).poles()) def test_horner(self, sys322): """Test horner() function""" @@ -853,7 +863,7 @@ def test_shape(self, states, outputs, inputs): def test_pole(self, states, outputs, inputs): """Test that the poles of rss outputs have a negative real part.""" sys = rss(states, outputs, inputs) - p = sys.pole() + p = sys.poles() for z in p: assert z.real < 0 @@ -905,7 +915,7 @@ def test_shape(self, states, outputs, inputs): def test_pole(self, states, outputs, inputs): """Test that the poles of drss outputs have less than unit magnitude.""" sys = drss(states, outputs, inputs) - p = sys.pole() + p = sys.poles() for z in p: assert abs(z) < 1 diff --git a/control/tests/timeresp_test.py b/control/tests/timeresp_test.py index 61c0cae38..fe73ab4a9 100644 --- a/control/tests/timeresp_test.py +++ b/control/tests/timeresp_test.py @@ -9,7 +9,7 @@ import control as ct from control import StateSpace, TransferFunction, c2d, isctime, ss2tf, tf2ss -from control.exception import slycot_check +from control.exception import slycot_check, pandas_check from control.tests.conftest import slycotonly from control.timeresp import (_default_time_vector, _ideal_tfinal_and_dt, forced_response, impulse_response, @@ -875,7 +875,7 @@ def test_time_vector(self, tsystem, fun, squeeze, matarrayout): kw['U'] = np.vstack([np.sin(t) for i in range(sys.ninputs)]) elif fun == forced_response and isctime(sys, strict=True): pytest.skip("No continuous forced_response without time vector.") - if hasattr(sys, "nstates"): + if hasattr(sys, "nstates") and sys.nstates is not None: kw['X0'] = np.arange(sys.nstates) + 1 if sys.ninputs > 1 and fun in [step_response, impulse_response]: kw['input'] = 1 @@ -1180,3 +1180,55 @@ def test_response_transpose( assert t.shape == (T.size, ) assert y.shape == ysh_no assert x.shape == (T.size, sys.nstates) + + +@pytest.mark.skipif(not pandas_check(), reason="pandas not installed") +def test_to_pandas(): + # Create a SISO time response + sys = ct.rss(2, 1, 1) + timepts = np.linspace(0, 10, 10) + resp = ct.input_output_response(sys, timepts, 1) + + # Convert to pandas + df = resp.to_pandas() + + # Check to make sure the data make senses + np.testing.assert_equal(df['time'], resp.time) + np.testing.assert_equal(df['u[0]'], resp.inputs) + np.testing.assert_equal(df['y[0]'], resp.outputs) + np.testing.assert_equal(df['x[0]'], resp.states[0]) + np.testing.assert_equal(df['x[1]'], resp.states[1]) + + # Create a MIMO time response + sys = ct.rss(2, 2, 1) + resp = ct.input_output_response(sys, timepts, np.sin(timepts)) + df = resp.to_pandas() + np.testing.assert_equal(df['time'], resp.time) + np.testing.assert_equal(df['u[0]'], resp.inputs[0]) + np.testing.assert_equal(df['y[0]'], resp.outputs[0]) + np.testing.assert_equal(df['y[1]'], resp.outputs[1]) + np.testing.assert_equal(df['x[0]'], resp.states[0]) + np.testing.assert_equal(df['x[1]'], resp.states[1]) + + # Change the time points + sys = ct.rss(2, 1, 1) + T = np.linspace(0, timepts[-1]/2, timepts.size * 2) + resp = ct.input_output_response(sys, timepts, np.sin(timepts), t_eval=T) + df = resp.to_pandas() + np.testing.assert_equal(df['time'], resp.time) + np.testing.assert_equal(df['u[0]'], resp.inputs) + np.testing.assert_equal(df['y[0]'], resp.outputs) + np.testing.assert_equal(df['x[0]'], resp.states[0]) + np.testing.assert_equal(df['x[1]'], resp.states[1]) + + +@pytest.mark.skipif(pandas_check(), reason="pandas installed") +def test_no_pandas(): + # Create a SISO time response + sys = ct.rss(2, 1, 1) + timepts = np.linspace(0, 10, 10) + resp = ct.input_output_response(sys, timepts, 1) + + # Convert to pandas + with pytest.raises(ImportError, match="pandas"): + df = resp.to_pandas() diff --git a/control/tests/type_conversion_test.py b/control/tests/type_conversion_test.py index d8c2d2b71..cdf302015 100644 --- a/control/tests/type_conversion_test.py +++ b/control/tests/type_conversion_test.py @@ -13,13 +13,13 @@ @pytest.fixture() def sys_dict(): sdict = {} - sdict['ss'] = ct.ss([[-1]], [[1]], [[1]], [[0]]) - sdict['tf'] = ct.tf([1],[0.5, 1]) - sdict['tfx'] = ct.tf([1, 1],[1]) # non-proper transfer function - sdict['frd'] = ct.frd([10+0j, 9 + 1j, 8 + 2j], [1,2,3]) + sdict['ss'] = ct.StateSpace([[-1]], [[1]], [[1]], [[0]]) + sdict['tf'] = ct.TransferFunction([1],[0.5, 1]) + sdict['tfx'] = ct.TransferFunction([1, 1], [1]) # non-proper TF + sdict['frd'] = ct.frd([10+0j, 9 + 1j, 8 + 2j, 7 + 3j], [1, 2, 3, 4]) sdict['lio'] = ct.LinearIOSystem(ct.ss([[-1]], [[5]], [[5]], [[0]])) sdict['ios'] = ct.NonlinearIOSystem( - sdict['lio']._rhs, sdict['lio']._out, 1, 1, 1) + sdict['lio']._rhs, sdict['lio']._out, inputs=1, outputs=1, states=1) sdict['arr'] = np.array([[2.0]]) sdict['flt'] = 3. return sdict @@ -59,39 +59,39 @@ def sys_dict(): rtype_list = ['ss', 'tf', 'frd', 'lio', 'ios', 'arr', 'flt'] conversion_table = [ # op left ss tf frd lio ios arr flt - ('add', 'ss', ['ss', 'ss', 'xrd', 'ss', 'ios', 'ss', 'ss' ]), - ('add', 'tf', ['tf', 'tf', 'xrd', 'tf', 'xos', 'tf', 'tf' ]), - ('add', 'frd', ['xrd', 'xrd', 'frd', 'xrd', 'E', 'xrd', 'xrd']), + ('add', 'ss', ['ss', 'ss', 'frd', 'ss', 'ios', 'ss', 'ss' ]), + ('add', 'tf', ['tf', 'tf', 'frd', 'lio', 'ios', 'tf', 'tf' ]), + ('add', 'frd', ['frd', 'frd', 'frd', 'frd', 'E', 'frd', 'frd']), ('add', 'lio', ['lio', 'lio', 'xrd', 'lio', 'ios', 'lio', 'lio']), ('add', 'ios', ['ios', 'ios', 'E', 'ios', 'ios', 'ios', 'ios']), - ('add', 'arr', ['ss', 'tf', 'xrd', 'lio', 'ios', 'arr', 'arr']), - ('add', 'flt', ['ss', 'tf', 'xrd', 'lio', 'ios', 'arr', 'flt']), + ('add', 'arr', ['ss', 'tf', 'frd', 'lio', 'ios', 'arr', 'arr']), + ('add', 'flt', ['ss', 'tf', 'frd', 'lio', 'ios', 'arr', 'flt']), # op left ss tf frd lio ios arr flt - ('sub', 'ss', ['ss', 'ss', 'xrd', 'ss', 'ios', 'ss', 'ss' ]), - ('sub', 'tf', ['tf', 'tf', 'xrd', 'tf', 'xos', 'tf', 'tf' ]), - ('sub', 'frd', ['xrd', 'xrd', 'frd', 'xrd', 'E', 'xrd', 'xrd']), + ('sub', 'ss', ['ss', 'ss', 'frd', 'ss', 'ios', 'ss', 'ss' ]), + ('sub', 'tf', ['tf', 'tf', 'frd', 'lio', 'ios', 'tf', 'tf' ]), + ('sub', 'frd', ['frd', 'frd', 'frd', 'frd', 'E', 'frd', 'frd']), ('sub', 'lio', ['lio', 'lio', 'xrd', 'lio', 'ios', 'lio', 'lio']), ('sub', 'ios', ['ios', 'ios', 'E', 'ios', 'ios', 'ios', 'ios']), - ('sub', 'arr', ['ss', 'tf', 'xrd', 'lio', 'ios', 'arr', 'arr']), - ('sub', 'flt', ['ss', 'tf', 'xrd', 'lio', 'ios', 'arr', 'flt']), + ('sub', 'arr', ['ss', 'tf', 'frd', 'lio', 'ios', 'arr', 'arr']), + ('sub', 'flt', ['ss', 'tf', 'frd', 'lio', 'ios', 'arr', 'flt']), # op left ss tf frd lio ios arr flt - ('mul', 'ss', ['ss', 'ss', 'xrd', 'ss', 'ios', 'ss', 'ss' ]), - ('mul', 'tf', ['tf', 'tf', 'xrd', 'tf', 'xos', 'tf', 'tf' ]), - ('mul', 'frd', ['xrd', 'xrd', 'frd', 'xrd', 'E', 'xrd', 'frd']), + ('mul', 'ss', ['ss', 'ss', 'frd', 'ss', 'ios', 'ss', 'ss' ]), + ('mul', 'tf', ['tf', 'tf', 'frd', 'lio', 'ios', 'tf', 'tf' ]), + ('mul', 'frd', ['frd', 'frd', 'frd', 'frd', 'E', 'frd', 'frd']), ('mul', 'lio', ['lio', 'lio', 'xrd', 'lio', 'ios', 'lio', 'lio']), ('mul', 'ios', ['ios', 'ios', 'E', 'ios', 'ios', 'ios', 'ios']), - ('mul', 'arr', ['ss', 'tf', 'xrd', 'lio', 'ios', 'arr', 'arr']), + ('mul', 'arr', ['ss', 'tf', 'frd', 'lio', 'ios', 'arr', 'arr']), ('mul', 'flt', ['ss', 'tf', 'frd', 'lio', 'ios', 'arr', 'flt']), # op left ss tf frd lio ios arr flt - ('truediv', 'ss', ['xs', 'tf', 'xrd', 'xio', 'xos', 'xs', 'xs' ]), + ('truediv', 'ss', ['xs', 'tf', 'frd', 'xio', 'xos', 'xs', 'xs' ]), ('truediv', 'tf', ['tf', 'tf', 'xrd', 'tf', 'xos', 'tf', 'tf' ]), - ('truediv', 'frd', ['xrd', 'xrd', 'frd', 'xrd', 'E', 'xrd', 'frd']), - ('truediv', 'lio', ['xio', 'tf', 'xrd', 'xio', 'xio', 'xio', 'xio']), + ('truediv', 'frd', ['frd', 'frd', 'frd', 'frd', 'E', 'frd', 'frd']), + ('truediv', 'lio', ['xio', 'tf', 'frd', 'xio', 'xio', 'xio', 'xio']), ('truediv', 'ios', ['xos', 'xos', 'E', 'xos', 'xos' 'xos', 'xos']), - ('truediv', 'arr', ['xs', 'tf', 'xrd', 'xio', 'xos', 'arr', 'arr']), + ('truediv', 'arr', ['xs', 'tf', 'frd', 'xio', 'xos', 'arr', 'arr']), ('truediv', 'flt', ['xs', 'tf', 'frd', 'xio', 'xos', 'arr', 'flt'])] # Now create list of the tests we actually want to run @@ -147,9 +147,8 @@ def test_operator_type_conversion(opname, ltype, rtype, expected, sys_dict): # Note: tfx = non-proper transfer function, order(num) > order(den) # -type_list = ['ss', 'tf', 'tfx', 'frd', 'lio', 'ios', 'arr', 'flt'] +type_list = ['ss', 'tf', 'tfx', 'frd', 'lio', 'ios', 'arr', 'flt'] conversion_table = [ - # L \ R ['ss', 'tf', 'tfx', 'frd', 'lio', 'ios', 'arr', 'flt'] ('ss', ['ss', 'ss', 'tf' 'frd', 'lio', 'ios', 'ss', 'ss' ]), ('tf', ['tf', 'tf', 'tf' 'frd', 'lio', 'ios', 'tf', 'tf' ]), ('tfx', ['tf', 'tf', 'tf', 'frd', 'E', 'E', 'tf', 'tf' ]), @@ -161,6 +160,7 @@ def test_operator_type_conversion(opname, ltype, rtype, expected, sys_dict): @pytest.mark.skip(reason="future test; conversions not yet fully implemented") # @pytest.mark.parametrize("opname", ['add', 'sub', 'mul', 'truediv']) +# @pytest.mark.parametrize("opname", ['add', 'sub', 'mul']) # @pytest.mark.parametrize("ltype", type_list) # @pytest.mark.parametrize("rtype", type_list) def test_binary_op_type_conversions(opname, ltype, rtype, sys_dict): @@ -185,3 +185,9 @@ def test_binary_op_type_conversions(opname, ltype, rtype, sys_dict): # Print out what we are testing in case something goes wrong assert isinstance(result, type_dict[expected]) + + # Make sure that input, output, and state names make sense + assert len(result.input_labels) == result.ninputs + assert len(result.output_labels) == result.noutputs + if result.nstates is not None: + assert len(result.state_labels) == result.nstates diff --git a/control/tests/xferfcn_test.py b/control/tests/xferfcn_test.py index 7821ce54d..79273f31b 100644 --- a/control/tests/xferfcn_test.py +++ b/control/tests/xferfcn_test.py @@ -39,7 +39,7 @@ def test_constructor_bad_input_type(self): TransferFunction([1]) # Too many arguments - with pytest.raises(ValueError): + with pytest.raises(TypeError): TransferFunction(1, 2, 3, 4) # Different numbers of elements in numerator rows @@ -85,18 +85,19 @@ def test_constructor_zero_denominator(self): TransferFunction([[[1.], [2., 3.]], [[-1., 4.], [3., 2.]]], [[[1., 0.], [0.]], [[0., 0.], [2.]]]) + @pytest.mark.skip("outdated test") def test_constructor_nodt(self): """Test the constructor when an object without dt is passed""" sysin = TransferFunction([[[0., 1.], [2., 3.]]], [[[5., 2.], [3., 0.]]]) - del sysin.dt + del sysin.dt # this doesn't make sense and now breaks sys = TransferFunction(sysin) assert sys.dt == defaults['control.default_dt'] # test for static gain sysin = TransferFunction([[[2.], [3.]]], [[[1.], [.1]]]) - del sysin.dt + del sysin.dt # this doesn't make sense and now breaks sys = TransferFunction(sysin) assert sys.dt is None @@ -596,7 +597,7 @@ def test_pole_mimo(self): sys = TransferFunction( [[[1.], [1.]], [[1.], [1.]]], [[[1., 2.], [1., 3.]], [[1., 4., 4.], [1., 9., 14.]]]) - p = sys.pole() + p = sys.poles() np.testing.assert_array_almost_equal(p, [-2., -2., -7., -3., -2.]) @@ -604,14 +605,14 @@ def test_pole_mimo(self): sys2 = TransferFunction( [[[1., 2., 3., 4.], [1.]], [[1.], [1.]]], [[[1., 2.], [1., 3.]], [[1., 4., 4.], [1., 9., 14.]]]) - p2 = sys2.pole() + p2 = sys2.poles() np.testing.assert_array_almost_equal(p2, [-2., -2., -7., -3., -2.]) def test_double_cancelling_poles_siso(self): H = TransferFunction([1, 1], [1, 2, 1]) - p = H.pole() + p = H.poles() np.testing.assert_array_almost_equal(p, [-1, -1]) # Tests for TransferFunction.feedback diff --git a/control/timeresp.py b/control/timeresp.py index bf826b539..aa1261ccd 100644 --- a/control/timeresp.py +++ b/control/timeresp.py @@ -79,7 +79,8 @@ from copy import copy from . import config -from .lti import isctime, isdtime +from .exception import pandas_check +from .namedio import isctime, isdtime from .statesp import StateSpace, _convert_to_statespace, _mimo2simo, _mimo2siso from .xferfcn import TransferFunction @@ -462,7 +463,7 @@ def __call__(self, **kwargs): # Update any keywords that we were passed response.transpose = kwargs.pop('transpose', self.transpose) response.squeeze = kwargs.pop('squeeze', self.squeeze) - response.return_x = kwargs.pop('return_x', self.squeeze) + response.return_x = kwargs.pop('return_x', self.return_x) # Check for new labels input_labels = kwargs.pop('input_labels', None) @@ -638,6 +639,23 @@ def __getitem__(self, index): def __len__(self): return 3 if self.return_x else 2 + # Convert to pandas + def to_pandas(self): + if not pandas_check(): + ImportError('pandas not installed') + import pandas + + # Create a dict for setting up the data frame + data = {'time': self.time} + data.update( + {name: self.u[i] for i, name in enumerate(self.input_labels)}) + data.update( + {name: self.y[i] for i, name in enumerate(self.output_labels)}) + data.update( + {name: self.x[i] for i, name in enumerate(self.state_labels)}) + + return pandas.DataFrame(data) + # Process signal labels def _process_labels(labels, signal, length): diff --git a/control/xferfcn.py b/control/xferfcn.py index a171a1143..93a66ce9d 100644 --- a/control/xferfcn.py +++ b/control/xferfcn.py @@ -59,8 +59,10 @@ from warnings import warn from itertools import chain from re import sub -from .lti import LTI, common_timebase, isdtime, _process_frequency_response +from .lti import LTI, _process_frequency_response +from .namedio import common_timebase, isdtime, _process_namedio_keywords from .exception import ControlMIMONotImplemented +from .frdata import FrequencyResponseData from . import config __all__ = ['TransferFunction', 'tf', 'ss2tf', 'tfdata'] @@ -161,13 +163,22 @@ def __init__(self, *args, **kwargs): (continuous or discrete). """ - args = deepcopy(args) + # + # Process positional arguments + # if len(args) == 2: # The user provided a numerator and a denominator. - (num, den) = args + num, den = args + elif len(args) == 3: # Discrete time transfer function - (num, den, dt) = args + num, den, dt = args + if 'dt' in kwargs: + warn("received multiple dt arguments, " + "using positional arg dt = %s" % dt) + kwargs['dt'] = dt + args = args[:-1] + elif len(args) == 1: # Use the copy constructor. if not isinstance(args[0], TransferFunction): @@ -176,43 +187,68 @@ def __init__(self, *args, **kwargs): % type(args[0])) num = args[0].num den = args[0].den + else: - raise ValueError("Needs 1, 2 or 3 arguments; received %i." + raise TypeError("Needs 1, 2 or 3 arguments; received %i." % len(args)) num = _clean_part(num) den = _clean_part(den) - inputs = len(num[0]) - outputs = len(num) + # + # Process keyword arguments + # + + # Determine if the transfer function is static (needed for dt) + static = True + for col in num + den: + for poly in col: + if len(poly) > 1: + static = False + defaults = args[0] if len(args) == 1 else \ + {'inputs': len(num[0]), 'outputs': len(num)} + + name, inputs, outputs, states, dt = _process_namedio_keywords( + kwargs, defaults, static=static, end=True) + if states: + raise TypeError( + "states keyword not allowed for transfer functions") + + # Initialize LTI (NamedIOSystem) object + super().__init__( + name=name, inputs=inputs, outputs=outputs, dt=dt) + + # + # Check to make sure everything is consistent + # # Make sure numerator and denominator matrices have consistent sizes - if inputs != len(den[0]): + if self.ninputs != len(den[0]): raise ValueError( "The numerator has %i input(s), but the denominator has " - "%i input(s)." % (inputs, len(den[0]))) - if outputs != len(den): + "%i input(s)." % (self.ninputs, len(den[0]))) + if self.noutputs != len(den): raise ValueError( "The numerator has %i output(s), but the denominator has " - "%i output(s)." % (outputs, len(den))) + "%i output(s)." % (self.noutputs, len(den))) # Additional checks/updates on structure of the transfer function - for i in range(outputs): + for i in range(self.noutputs): # Make sure that each row has the same number of columns - if len(num[i]) != inputs: + if len(num[i]) != self.ninputs: raise ValueError( "Row 0 of the numerator matrix has %i elements, but row " - "%i has %i." % (inputs, i, len(num[i]))) - if len(den[i]) != inputs: + "%i has %i." % (self.ninputs, i, len(num[i]))) + if len(den[i]) != self.ninputs: raise ValueError( "Row 0 of the denominator matrix has %i elements, but row " - "%i has %i." % (inputs, i, len(den[i]))) + "%i has %i." % (self.ninputs, i, len(den[i]))) # Check for zeros in numerator or denominator # TODO: Right now these checks are only done during construction. # It might be worthwhile to think of a way to perform checks if the # user modifies the transfer function after construction. - for j in range(inputs): + for j in range(self.ninputs): # Check that we don't have any zero denominators. zeroden = True for k in den[i][j]: @@ -233,42 +269,16 @@ def __init__(self, *args, **kwargs): if zeronum: den[i][j] = ones(1) - super().__init__(inputs, outputs) + # Store the numerator and denominator self.num = num self.den = den + # + # Final processing + # + # Truncate leading zeros self._truncatecoeff() - # get dt - if len(args) == 2: - # no dt given in positional arguments - if 'dt' in kwargs: - dt = kwargs.pop('dt') - elif self._isstatic(): - dt = None - else: - dt = config.defaults['control.default_dt'] - elif len(args) == 3: - # Discrete time transfer function - if 'dt' in kwargs: - warn('received multiple dt arguments, ' - 'using positional arg dt=%s' % dt) - kwargs.pop('dt') - elif len(args) == 1: - # TODO: not sure this can ever happen since dt is always present - try: - dt = args[0].dt - except AttributeError: - if self._isstatic(): - dt = None - else: - dt = config.defaults['control.default_dt'] - self.dt = dt - - # Make sure there were no extraneous keywords - if kwargs: - raise TypeError("unrecognized keywords: ", str(kwargs)) - # # Class attributes # @@ -528,6 +538,12 @@ def __add__(self, other): """Add two LTI objects (parallel connection).""" from .statesp import StateSpace + # Check to see if the right operator has priority + if getattr(other, '__array_priority__', None) and \ + getattr(self, '__array_priority__', None) and \ + other.__array_priority__ > self.__array_priority__: + return other.__radd__(self) + # Convert the second argument to a transfer function. if isinstance(other, StateSpace): other = _convert_to_transfer_function(other) @@ -573,6 +589,12 @@ def __rsub__(self, other): def __mul__(self, other): """Multiply two LTI objects (serial connection).""" + # Check to see if the right operator has priority + if getattr(other, '__array_priority__', None) and \ + getattr(self, '__array_priority__', None) and \ + other.__array_priority__ > self.__array_priority__: + return other.__rmul__(self) + # Convert the second argument to a transfer function. if isinstance(other, (int, float, complex, np.number)): other = _convert_to_transfer_function(other, inputs=self.ninputs, @@ -770,7 +792,7 @@ def freqresp(self, omega): "MATLAB compatibility module instead", DeprecationWarning) return self.frequency_response(omega) - def pole(self): + def poles(self): """Compute the poles of a transfer function.""" _, den, denorder = self._common_den(allow_nonproper=True) rts = [] @@ -778,7 +800,7 @@ def pole(self): rts.extend(roots(d[:o + 1])) return np.array(rts) - def zero(self): + def zeros(self): """Compute the zeros of a transfer function.""" if self.ninputs > 1 or self.noutputs > 1: raise NotImplementedError( @@ -1225,10 +1247,9 @@ def _c2d_matched(sysC, Ts): sysDnum, sysDden = zpk2tf(zzeros, zpoles, gain) return TransferFunction(sysDnum, sysDden, Ts) + # Utility function to convert a transfer function polynomial to a string # Borrowed from poly1d library - - def _tf_polynomial_to_string(coeffs, var='s'): """Convert a transfer function polynomial to a string""" @@ -1318,6 +1339,9 @@ def _convert_to_transfer_function(sys, inputs=1, outputs=1): If sys is an array-like type, then it is converted to a constant-gain transfer function. + Note: no renaming of inputs and outputs is performed; this should be done + by the calling function. + >>> sys = _convert_to_transfer_function([[1., 0.], [2., 3.]]) In this example, the numerator matrix will be @@ -1326,6 +1350,7 @@ def _convert_to_transfer_function(sys, inputs=1, outputs=1): """ from .statesp import StateSpace + kwargs = {} if isinstance(sys, TransferFunction): return sys @@ -1373,13 +1398,19 @@ def _convert_to_transfer_function(sys, inputs=1, outputs=1): num = squeeze(num) # Convert to 1D array den = squeeze(den) # Probably not needed - return TransferFunction(num, den, sys.dt) + return TransferFunction( + num, den, sys.dt, inputs=sys.input_labels, + outputs=sys.output_labels) elif isinstance(sys, (int, float, complex, np.number)): num = [[[sys] for j in range(inputs)] for i in range(outputs)] den = [[[1] for j in range(inputs)] for i in range(outputs)] - return TransferFunction(num, den) + return TransferFunction( + num, den, inputs=inputs, outputs=outputs) + + elif isinstance(sys, FrequencyResponseData): + raise TypeError("Can't convert given FRD to TransferFunction system.") # If this is array-like, try to create a constant feedthrough try: @@ -1388,6 +1419,7 @@ def _convert_to_transfer_function(sys, inputs=1, outputs=1): num = [[[D[i, j]] for j in range(inputs)] for i in range(outputs)] den = [[[1] for j in range(inputs)] for i in range(outputs)] return TransferFunction(num, den) + except Exception: raise TypeError("Can't convert given type to TransferFunction system.") @@ -1437,6 +1469,16 @@ def tf(*args, **kwargs): out: :class:`TransferFunction` The new linear system + Other Parameters + ---------------- + inputs, outputs : str, or list of str, optional + List of strings that name the individual signals of the transformed + system. If not given, the inputs and outputs are the same as the + original system. + name : string, optional + System name. If unspecified, a generic name is generated + with a unique integer id. + Raises ------ ValueError @@ -1483,7 +1525,8 @@ def tf(*args, **kwargs): if len(args) == 2 or len(args) == 3: return TransferFunction(*args, **kwargs) - elif len(args) == 1: + + elif len(args) == 1 and isinstance(args[0], str): # Make sure there were no extraneous keywords if kwargs: raise TypeError("unrecognized keywords: ", str(kwargs)) @@ -1494,12 +1537,14 @@ def tf(*args, **kwargs): elif args[0] == 'z': return TransferFunction.z + elif len(args) == 1: from .statesp import StateSpace sys = args[0] if isinstance(sys, StateSpace): - return ss2tf(sys) + return ss2tf(sys, **kwargs) elif isinstance(sys, TransferFunction): - return deepcopy(sys) + # Use copy constructor + return TransferFunction(sys, **kwargs) else: raise TypeError("tf(sys): sys must be a StateSpace or " "TransferFunction object. It is %s." % type(sys)) @@ -1542,6 +1587,16 @@ def ss2tf(*args, **kwargs): out: TransferFunction New linear system in transfer function form + Other Parameters + ---------------- + inputs, outputs : str, or list of str, optional + List of strings that name the individual signals of the transformed + system. If not given, the inputs and outputs are the same as the + original system. + name : string, optional + System name. If unspecified, a generic name is generated + with a unique integer id. + Raises ------ ValueError @@ -1574,14 +1629,11 @@ def ss2tf(*args, **kwargs): # Assume we were given the A, B, C, D matrix and (optional) dt return _convert_to_transfer_function(StateSpace(*args, **kwargs)) - # Make sure there were no extraneous keywords - if kwargs: - raise TypeError("unrecognized keywords: ", str(kwargs)) - if len(args) == 1: sys = args[0] if isinstance(sys, StateSpace): - return _convert_to_transfer_function(sys) + return TransferFunction( + _convert_to_transfer_function(sys), **kwargs) else: raise TypeError( "ss2tf(sys): sys must be a StateSpace object. It is %s." diff --git a/doc/.gitignore b/doc/.gitignore new file mode 100644 index 000000000..d948f64d2 --- /dev/null +++ b/doc/.gitignore @@ -0,0 +1 @@ +*.fig.bak diff --git a/doc/Makefile b/doc/Makefile index b72312be4..3f372684c 100644 --- a/doc/Makefile +++ b/doc/Makefile @@ -14,7 +14,11 @@ help: .PHONY: help Makefile +# Rules to create figures +FIGS = classes.pdf +classes.pdf: classes.fig; fig2dev -Lpdf $< $@ + # Catch-all target: route all unknown targets to Sphinx using the new # "make mode" option. $(O) is meant as a shortcut for $(SPHINXOPTS). -%: Makefile - @$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) \ No newline at end of file +html pdf: Makefile $(FIGS) + @$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) diff --git a/doc/classes.fig b/doc/classes.fig new file mode 100644 index 000000000..950510c01 --- /dev/null +++ b/doc/classes.fig @@ -0,0 +1,149 @@ +#FIG 3.2 Produced by xfig version 3.2.8b +Landscape +Center +Inches +Letter +100.00 +Single +-2 +1200 2 +2 2 1 1 0 7 50 -1 -1 4.000 0 0 -1 0 0 5 + 9750 3375 12075 3375 12075 4725 9750 4725 9750 3375 +2 2 1 1 0 7 50 -1 -1 4.000 0 0 -1 0 0 5 + 9750 6000 12075 6000 12075 7350 9750 7350 9750 6000 +2 1 0 2 4 7 50 -1 -1 0.000 0 0 7 1 0 2 + 1 1 1.00 60.00 120.00 + 8925 3600 9750 3600 +2 1 0 2 4 7 50 -1 -1 0.000 0 0 7 1 0 2 + 1 1 1.00 60.00 120.00 + 10875 3750 10875 4350 +2 1 0 2 4 7 50 -1 -1 0.000 0 0 7 1 0 2 + 1 1 1.00 60.00 120.00 + 6375 3750 9975 6150 +2 1 0 2 4 7 50 -1 -1 0.000 0 0 7 1 0 2 + 1 1 1.00 60.00 120.00 + 10875 6375 10875 6975 +2 1 0 2 4 7 50 -1 -1 0.000 0 0 7 1 0 2 + 1 1 1.00 60.00 120.00 + 6750 6225 9975 6225 +2 1 0 2 4 7 50 -1 -1 0.000 0 0 7 1 0 2 + 1 1 1.00 60.00 120.00 + 6000 6075 6000 6975 +2 1 0 2 4 7 50 -1 -1 0.000 0 0 7 1 0 2 + 1 1 1.00 60.00 120.00 + 2700 5400 3075 5850 +2 2 1 1 0 7 50 -1 -1 4.000 0 0 -1 0 0 5 + 1650 4500 6750 4500 6750 7425 1650 7425 1650 4500 +2 2 1 1 0 7 50 -1 -1 4.000 0 0 -1 0 0 5 + 1650 7950 6150 7950 6150 8550 1650 8550 1650 7950 +2 1 0 2 4 7 50 -1 -1 0.000 0 0 -1 1 0 2 + 1 1 1.00 60.00 120.00 + 2775 8175 4200 8175 +2 1 0 2 4 7 50 -1 -1 0.000 0 0 -1 1 0 2 + 1 1 1.00 60.00 120.00 + 9075 8100 9675 8100 +2 2 1 1 0 7 50 -1 -1 4.000 0 0 -1 0 0 5 + 9075 8250 9675 8250 9675 8550 9075 8550 9075 8250 +2 1 0 2 4 7 50 -1 -1 0.000 0 0 7 0 1 2 + 1 1 1.00 60.00 120.00 + 4725 5925 5175 5925 +2 1 0 2 4 7 50 -1 -1 0.000 0 0 7 1 1 2 + 1 1 1.00 60.00 120.00 + 1 1 1.00 60.00 120.00 + 6525 3600 7275 3600 +2 1 0 2 4 7 50 -1 -1 0.000 0 0 7 1 0 2 + 1 1 1.00 60.00 120.00 + 5775 8175 9975 6300 +2 2 1 1 0 7 50 -1 -1 4.000 0 0 -1 0 0 5 + 5400 3375 6600 3375 6600 3900 5400 3900 5400 3375 +2 2 1 1 0 7 50 -1 -1 4.000 0 0 -1 0 0 5 + 7050 2175 8100 2175 8100 2700 7050 2700 7050 2175 +2 2 1 1 1 7 50 -1 -1 4.000 0 0 -1 0 0 5 + 4500 975 6525 975 6525 1500 4500 1500 4500 975 +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 + 5250 1350 3825 4575 +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 + 5775 1350 7575 2250 +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 + 7875 2550 10875 3450 +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 + 7575 2550 8025 3450 +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 + 7350 2550 6225 3450 +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 + 3300 4875 3000 5100 +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 + 3825 4875 3825 5775 +2 1 0 2 4 7 50 -1 -1 0.000 0 0 7 1 0 2 + 1 1 1.00 60.00 120.00 + 4350 4875 5625 5775 +2 2 1 1 0 7 50 -1 -1 4.000 0 0 -1 0 0 5 + 7350 3375 8925 3375 8925 3900 7350 3900 7350 3375 +2 1 0 2 1 7 50 -1 -1 0.000 0 0 -1 0 1 2 + 1 0 1.00 60.00 90.00 + 9075 7800 9675 7800 +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 + 4350 6075 5625 6975 +2 1 0 2 1 7 50 -1 -1 0.000 0 0 -1 0 1 2 + 1 0 1.00 60.00 90.00 + 2400 5400 2400 8025 +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 + 5850 6075 5850 6975 +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 + 4125 4875 5400 5775 +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 + 5925 3750 5925 5775 +4 0 0 50 -1 0 12 0.0000 4 165 885 5400 3300 statesp.py\001 +4 0 0 50 -1 0 12 0.0000 4 195 420 8175 2325 lti.py\001 +4 2 0 50 -1 0 12 0.0000 4 195 885 8925 3300 xferfcn.py\001 +4 2 0 50 -1 0 12 0.0000 4 195 780 12075 3300 frdata.py\001 +4 2 0 50 -1 0 12 0.0000 4 195 780 12075 5925 trdata.py\001 +4 1 1 50 -1 0 12 0.0000 4 150 345 7575 2475 LTI\001 +4 1 1 50 -1 0 12 0.0000 4 195 1440 5925 6000 LinearIOSystem\001 +4 0 0 50 -1 0 12 0.0000 4 195 615 1650 7875 flatsys/\001 +4 0 0 50 -1 0 12 0.0000 4 195 705 1650 4425 iosys.py\001 +4 0 0 50 -1 0 12 0.0000 4 195 720 8700 7575 Legend:\001 +4 1 1 50 -1 16 12 0.0000 4 210 1590 5475 1275 NamedIOSystem\001 +4 1 1 50 -1 16 12 0.0000 4 210 1770 3975 4800 InputOutputSystem\001 +4 1 1 50 -1 16 12 0.0000 4 210 1830 2625 5325 NonlinearIOSystem\001 +4 0 0 50 -1 0 12 0.0000 4 195 1005 6600 1125 namedio.py\001 +4 0 4 50 -1 16 12 0.0000 4 210 945 4800 5100 linearize()\001 +4 1 1 50 -1 16 12 0.0000 4 210 2115 3750 6000 InterconnectedSystem\001 +4 0 4 50 -1 16 12 0.0000 4 210 1875 3000 6750 ic() = interconnect()\001 +4 1 1 50 -1 16 12 0.0000 4 210 1500 5925 7200 LinearICSystem\001 +4 1 1 50 -1 16 12 0.0000 4 210 1035 2250 8250 FlatSystem\001 +4 1 4 50 -1 16 12 0.0000 4 210 1500 3525 8400 point_to_point()\001 +4 1 1 50 -1 16 12 0.0000 4 210 1095 6000 3675 StateSpace\001 +4 1 1 50 -1 16 12 0.0000 4 165 1605 8100 3675 TransferFunction\001 +4 1 1 50 -1 16 12 0.0000 4 210 2400 10875 3675 FrequencyResponseData\001 +4 0 4 50 -1 16 12 0.0000 4 210 1155 10950 4050 to_pandas()\001 +4 1 1 50 -1 16 12 0.0000 4 210 1800 10875 4575 pandas.DataFrame\001 +4 0 4 50 -1 16 12 0.0000 4 210 1560 7950 4725 step_response()\001 +4 0 4 50 -1 16 12 0.0000 4 210 1635 8400 5025 initial_response()\001 +4 0 4 50 -1 16 12 0.0000 4 210 1755 8850 5325 forced_response()\001 +4 1 1 50 -1 16 12 0.0000 4 210 1875 10875 6300 TimeResponseData\001 +4 0 4 50 -1 16 12 0.0000 4 210 1155 10950 6675 to_pandas()\001 +4 1 1 50 -1 16 12 0.0000 4 210 1800 10875 7200 pandas.DataFrame\001 +4 0 1 50 -1 16 12 0.0000 4 210 1755 9750 7875 Class dependency\001 +4 0 4 50 -1 16 12 0.0000 4 210 2475 9750 8175 Conversion [via function()]\001 +4 0 0 50 -1 0 12 0.0000 4 150 1380 9750 8475 Source code file\001 +4 1 4 50 -1 16 12 0.0000 4 210 300 3150 5625 ic()\001 +4 0 4 50 -1 16 12 0.0000 4 210 300 6075 6600 ic()\001 +4 1 1 50 -1 16 12 0.0000 4 210 1650 4950 8250 SystemTrajectory\001 +4 1 4 50 -1 16 12 0.0000 4 210 945 9375 3825 freqresp()\001 +4 1 4 50 -1 16 12 0.0000 4 210 600 6975 3825 tf2ss()\001 +4 1 4 50 -1 16 12 0.0000 4 210 600 6975 3450 ss2tf()\001 +4 1 4 50 -1 16 12 0.0000 4 210 300 5025 6150 ic()\001 +4 1 4 50 -1 16 12 0.0000 4 210 2295 8325 6075 input_output_response()\001 +4 2 4 50 -1 16 12 0.0000 4 210 1035 8175 6975 response()\001 diff --git a/doc/classes.pdf b/doc/classes.pdf new file mode 100644 index 000000000..66ef25e10 Binary files /dev/null and b/doc/classes.pdf differ diff --git a/doc/classes.rst b/doc/classes.rst index 0753271c4..87ce457de 100644 --- a/doc/classes.rst +++ b/doc/classes.rst @@ -5,20 +5,30 @@ Control system classes ********************** -The classes listed below are used to represent models of linear time-invariant -(LTI) systems. 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 classes listed below are used to represent models of input/output +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. .. autosummary:: :toctree: generated/ :template: custom-class-template.rst - TransferFunction StateSpace + TransferFunction + InputOutputSystem FrequencyResponseData TimeResponseData +The following figure illustrates the relationship between the classes and +some of the functions that can be used to convert objects from one class to +another: + +.. image:: classes.pdf + :width: 800 + +| + Input/output system subclasses ============================== Input/output systems are accessed primarily via a set of subclasses diff --git a/doc/control.rst b/doc/control.rst index 20f363a1e..fc6618d24 100644 --- a/doc/control.rst +++ b/doc/control.rst @@ -83,8 +83,8 @@ Control system analysis margin stability_margins phase_crossover_frequencies - pole - zero + poles + zeros pzmap root_locus sisotool diff --git a/doc/conventions.rst b/doc/conventions.rst index 462a71408..1832b9525 100644 --- a/doc/conventions.rst +++ b/doc/conventions.rst @@ -29,9 +29,9 @@ of linear time-invariant (LTI) systems: where u is the input, y is the output, and x is the state. -To create a state space system, use the :class:`StateSpace` constructor: +To create a state space system, use the :fun:`ss` function: - sys = StateSpace(A, B, C, D) + sys = ct.ss(A, B, C, D) State space systems can be manipulated using standard arithmetic operations as well as the :func:`feedback`, :func:`parallel`, and :func:`series` @@ -51,10 +51,9 @@ transfer functions where n is generally greater than or equal to m (for a proper transfer function). -To create a transfer function, use the :class:`TransferFunction` -constructor: +To create a transfer function, use the :func:`tf` function: - sys = TransferFunction(num, den) + sys = ct.tf(num, den) Transfer functions can be manipulated using standard arithmetic operations as well as the :func:`feedback`, :func:`parallel`, and :func:`series` @@ -75,6 +74,28 @@ FRD systems have a somewhat more limited set of functions that are available, although all of the standard algebraic manipulations can be performed. +The FRD class is also used as the return type for the +:func:`frequency_response` function (and the equivalent method for the +:class:`StateSpace` and :class:`TransferFunction` classes). This +object can be assigned to a tuple using + + mag, phase, omega = response + +where `mag` is the magnitude (absolute value, not dB or log10) of the +system frequency response, `phase` is the wrapped phase in radians of +the system frequency response, and `omega` is the (sorted) frequencies +at which the response was evaluated. If the system is SISO and the +`squeeze` argument to :func:`frequency_response` is not True, +`magnitude` and `phase` are 1D, indexed by frequency. If the system +is not SISO or `squeeze` is False, the array is 3D, indexed by the +output, input, and frequency. If `squeeze` is True then +single-dimensional axes are removed. The processing of the `squeeze` +keyword can be changed by calling the response function with a new +argument: + + mag, phase, omega = response(squeeze=False) + + Discrete time systems --------------------- A discrete time system is created by specifying a nonzero 'timebase', dt. @@ -89,14 +110,16 @@ Only the :class:`StateSpace`, :class:`TransferFunction`, and :class:`InputOutputSystem` classes allow explicit representation of discrete time systems. -Systems must have compatible timebases in order to be combined. A discrete time -system with unspecified sampling time (`dt = True`) can be combined with a system -having a specified sampling time; the result will be a discrete time system with the sample time of the latter -system. Similarly, a system with timebase `None` can be combined with a system having a specified -timebase; the result will have the timebase of the latter system. For continuous -time systems, the :func:`sample_system` function or the :meth:`StateSpace.sample` and :meth:`TransferFunction.sample` methods -can be used to create a discrete time system from a continuous time system. -See :ref:`utility-and-conversions`. The default value of 'dt' can be changed by +Systems must have compatible timebases in order to be combined. A discrete +time system with unspecified sampling time (`dt = True`) can be combined with +a system having a specified sampling time; the result will be a discrete time +system with the sample time of the latter system. Similarly, a system with +timebase `None` can be combined with a system having a specified timebase; the +result will have the timebase of the latter system. For continuous time +systems, the :func:`sample_system` function or the :meth:`StateSpace.sample` +and :meth:`TransferFunction.sample` methods can be used to create a discrete +time system from a continuous time system. See +:ref:`utility-and-conversions`. The default value of 'dt' can be changed by changing the value of ``control.config.defaults['control.default_dt']``. Conversion between representations @@ -129,11 +152,6 @@ and :func:`initial_response`. Thus, all 2D values must be transposed when they are used with functions from `scipy.signal`_. -Types: - - * **Arguments** can be **arrays**, **matrices**, or **nested lists**. - * **Return values** are **arrays** (not matrices). - The time vector is a 1D array with shape (n, ):: T = [t1, t2, t3, ..., tn ] @@ -170,8 +188,8 @@ Functions that return time responses (e.g., :func:`forced_response`, response. These data can be accessed via the ``time``, ``outputs``, ``states`` and ``inputs`` properties:: - sys = rss(4, 1, 1) - response = step_response(sys) + sys = ct.rss(4, 1, 1) + response = ct.step_response(sys) plot(response.time, response.outputs) The dimensions of the response properties depend on the function being @@ -185,12 +203,12 @@ The time response functions can also be assigned to a tuple, which extracts the time and output (and optionally the state, if the `return_x` keyword is used). This allows simple commands for plotting:: - t, y = step_response(sys) + t, y = ct.step_response(sys) plot(t, y) -The output of a MIMO system can be plotted like this:: +The output of a MIMO LTI system can be plotted like this:: - t, y = forced_response(sys, t, u) + t, y = ct.forced_response(sys, t, u) plot(t, y[0], label='y_0') plot(t, y[1], label='y_1') @@ -201,6 +219,16 @@ response, can be computed like this:: ft = D @ U +Finally, the `to_pandas()` function can be used to create a pandas dataframe: + + df = response.to_pandas() + +The column labels for the data frame are `time` and the labels for the input, +output, and state signals (`u[i]`, `y[i]`, and `x[i]` by default, but these +can be changed using the `inputs`, `outputs`, and `states` keywords when +constructing the system, as described in :func:`ss`, :func:`tf`, and other +system creation function. Note that when exporting to pandas, "rows" in the +data frame correspond to time and "cols" (DataSeries) correspond to signals. .. currentmodule:: control .. _package-configuration-parameters: @@ -218,14 +246,14 @@ element of the `control.config.defaults` dictionary: .. code-block:: python - control.config.defaults['module.parameter'] = value + ct.config.defaults['module.parameter'] = value The `~control.config.set_defaults` function can also be used to set multiple configuration parameters at the same time: .. code-block:: python - control.config.set_defaults('module', param1=val1, param2=val2, ...] + ct.config.set_defaults('module', param1=val1, param2=val2, ...] Finally, there are also functions available set collections of variables based on standard configurations. diff --git a/doc/iosys.rst b/doc/iosys.rst index 41e37cfec..1da7f5884 100644 --- a/doc/iosys.rst +++ b/doc/iosys.rst @@ -13,25 +13,22 @@ The dynamics of the system can be in continuous or discrete time. To simulate an input/output system, use the :func:`~control.input_output_response` function:: - t, y = input_output_response(io_sys, T, U, X0, params) + t, y = ct.input_output_response(io_sys, T, U, X0, params) An input/output system can be linearized around an equilibrium point to obtain a :class:`~control.StateSpace` linear system. Use the :func:`~control.find_eqpt` function to obtain an equilibrium point and the :func:`~control.linearize` function to linearize about that equilibrium point:: - xeq, ueq = find_eqpt(io_sys, X0, U0) - ss_sys = linearize(io_sys, xeq, ueq) + xeq, ueq = ct.find_eqpt(io_sys, X0, U0) + ss_sys = ct.linearize(io_sys, xeq, ueq) -Input/output systems can be created from state space LTI systems by using the -:class:`~control.LinearIOSystem` class`:: - - io_sys = LinearIOSystem(ss_sys) - -Nonlinear input/output systems can be created using the -:class:`~control.NonlinearIOSystem` class, which requires the definition of an -update function (for the right hand side of the differential or different -equation) and and output function (computes the outputs from the state):: +Input/output systems are automatically created for state space LTI systems +when using the :func:`ss` function. Nonlinear input/output systems can be +created using the :class:`~control.NonlinearIOSystem` class, which requires +the definition of an update function (for the right hand side of the +differential or different equation) and an output function (computes the +outputs from the state):: io_sys = NonlinearIOSystem(updfcn, outfcn, inputs=M, outputs=P, states=N) @@ -64,7 +61,7 @@ We begin by defining the dynamics of the system .. code-block:: python - import control + import control as ct import numpy as np import matplotlib.pyplot as plt @@ -94,7 +91,7 @@ We now create an input/output system using these dynamics: .. code-block:: python - io_predprey = control.NonlinearIOSystem( + io_predprey = ct.NonlinearIOSystem( predprey_rhs, None, inputs=('u'), outputs=('H', 'L'), states=('H', 'L'), name='predprey') @@ -110,7 +107,7 @@ of the system: T = np.linspace(0, 70, 500) # Simulation 70 years of time # Simulate the system - t, y = control.input_output_response(io_predprey, T, 0, X0) + t, y = ct.input_output_response(io_predprey, T, 0, X0) # Plot the response plt.figure(1) @@ -125,9 +122,9 @@ system and computing the linearization about that point. .. code-block:: python - eqpt = control.find_eqpt(io_predprey, X0, 0) + eqpt = ct.find_eqpt(io_predprey, X0, 0) xeq = eqpt[0] # choose the nonzero equilibrium point - lin_predprey = control.linearize(io_predprey, xeq, 0) + lin_predprey = ct.linearize(io_predprey, xeq, 0) We next compute a controller that stabilizes the equilibrium point using eigenvalue placement and computing the feedforward gain using the number of @@ -135,7 +132,7 @@ lynxes as the desired output (following FBS2e, Example 7.5): .. code-block:: python - K = control.place(lin_predprey.A, lin_predprey.B, [-0.1, -0.2]) + K = ct.place(lin_predprey.A, lin_predprey.B, [-0.1, -0.2]) A, B = lin_predprey.A, lin_predprey.B C = np.array([[0, 1]]) # regulated output = number of lynxes kf = -1/(C @ np.linalg.inv(A - B @ K) @ B) @@ -147,7 +144,7 @@ constructed using the `~control.ios.NonlinearIOSystem` class: .. code-block:: python - io_controller = control.NonlinearIOSystem( + io_controller = ct.NonlinearIOSystem( None, lambda t, x, u, params: -K @ (u[1:] - xeq) + kf * (u[0] - xeq[1]), inputs=('Ld', 'u1', 'u2'), outputs=1, name='control') @@ -161,7 +158,7 @@ function: .. code-block:: python - io_closed = control.interconnect( + io_closed = ct.interconnect( [io_predprey, io_controller], # systems connections=[ ['predprey.u', 'control.y[0]'], @@ -177,7 +174,7 @@ Finally, we simulate the closed loop system: .. code-block:: python # Simulate the system - t, y = control.input_output_response(io_closed, T, 30, [15, 20]) + t, y = ct.input_output_response(io_closed, T, 30, [15, 20]) # Plot the response plt.figure(2) @@ -245,10 +242,10 @@ interconnecting systems, especially when combined with the :func:`~control.summing_junction` function. For example, the following code will create a unity gain, negative feedback system:: - P = control.tf2io(control.tf(1, [1, 0]), inputs='u', outputs='y') - C = control.tf2io(control.tf(10, [1, 1]), inputs='e', outputs='u') - sumblk = control.summing_junction(inputs=['r', '-y'], output='e') - T = control.interconnect([P, C, sumblk], inplist='r', outlist='y') + P = ct.tf2io([1], [1, 0], inputs='u', outputs='y') + C = ct.tf2io([10], [1, 1], inputs='e', outputs='u') + sumblk = ct.summing_junction(inputs=['r', '-y'], output='e') + T = ct.interconnect([P, C, sumblk], inplist='r', outlist='y') If a signal name appears in multiple outputs then that signal will be summed when it is interconnected. Similarly, if a signal name appears in multiple diff --git a/doc/optimal.rst b/doc/optimal.rst index 8da08e7af..bb952e9cc 100644 --- a/doc/optimal.rst +++ b/doc/optimal.rst @@ -99,7 +99,10 @@ The optimal control module provides a means of computing optimal trajectories for nonlinear systems and implementing optimization-based controllers, including model predictive control. It follows the basic problem setup described above, but carries out all computations in *discrete -time* (so that integrals become sums) and over a *finite horizon*. +time* (so that integrals become sums) and over a *finite horizon*. To local +the optimal control modules, import `control.optimal`: + + import control.optimal as obc To describe an optimal control problem we need an input/output system, a time horizon, a cost function, and (optionally) a set of constraints on the diff --git a/examples/pvtol-lqr.py b/examples/pvtol-lqr.py index 8654c77ad..8a9ff55d9 100644 --- a/examples/pvtol-lqr.py +++ b/examples/pvtol-lqr.py @@ -9,8 +9,8 @@ import os import numpy as np -import matplotlib.pyplot as plt # MATLAB plotting functions -from control.matlab import * # MATLAB-like functions +import matplotlib.pyplot as plt # MATLAB-like plotting functions +import control as ct # # System dynamics @@ -28,14 +28,13 @@ # State space dynamics xe = [0, 0, 0, 0, 0, 0] # equilibrium point of interest -ue = [0, m*g] # (note these are lists, not matrices) +ue = [0, m * g] # (note these are lists, not matrices) # TODO: The following objects need converting from np.matrix to np.array # This will involve re-working the subsequent equations as the shapes # See below. -# Dynamics matrix (use matrix type so that * works for multiplication) -A = np.matrix( +A = np.array( [[0, 0, 0, 1, 0, 0], [0, 0, 0, 0, 1, 0], [0, 0, 0, 0, 0, 1], @@ -45,7 +44,7 @@ ) # Input matrix -B = np.matrix( +B = np.array( [[0, 0], [0, 0], [0, 0], [np.cos(xe[2])/m, -np.sin(xe[2])/m], [np.sin(xe[2])/m, np.cos(xe[2])/m], @@ -53,8 +52,8 @@ ) # Output matrix -C = np.matrix([[1, 0, 0, 0, 0, 0], [0, 1, 0, 0, 0, 0]]) -D = np.matrix([[0, 0], [0, 0]]) +C = np.array([[1, 0, 0, 0, 0, 0], [0, 1, 0, 0, 0, 0]]) +D = np.array([[0, 0], [0, 0]]) # # Construct inputs and outputs corresponding to steps in xy position @@ -74,8 +73,8 @@ # so that xd corresponds to the desired steady state. # -xd = np.matrix([[1], [0], [0], [0], [0], [0]]) -yd = np.matrix([[0], [1], [0], [0], [0], [0]]) +xd = np.array([[1], [0], [0], [0], [0], [0]]) +yd = np.array([[0], [1], [0], [0], [0], [0]]) # # Extract the relevant dynamics for use with SISO library @@ -93,14 +92,14 @@ # Decoupled dynamics Ax = A[np.ix_(lat, lat)] -Bx = B[lat, 0] -Cx = C[0, lat] -Dx = D[0, 0] +Bx = B[np.ix_(lat, [0])] +Cx = C[np.ix_([0], lat)] +Dx = D[np.ix_([0], [0])] Ay = A[np.ix_(alt, alt)] -By = B[alt, 1] -Cy = C[1, alt] -Dy = D[1, 1] +By = B[np.ix_(alt, [1])] +Cy = C[np.ix_([1], alt)] +Dy = D[np.ix_([1], [1])] # Label the plot plt.clf() @@ -113,44 +112,24 @@ # Start with a diagonal weighting Qx1 = np.diag([1, 1, 1, 1, 1, 1]) Qu1a = np.diag([1, 1]) -K, X, E = lqr(A, B, Qx1, Qu1a) -K1a = np.matrix(K) +K1a, X, E = ct.lqr(A, B, Qx1, Qu1a) # Close the loop: xdot = Ax - B K (x-xd) +# # Note: python-control requires we do this 1 input at a time # H1a = ss(A-B*K1a, B*K1a*concatenate((xd, yd), axis=1), C, D); -# (T, Y) = step(H1a, T=np.linspace(0,10,100)); - -# TODO: The following equations will need modifying when converting from np.matrix to np.array -# because the results and even intermediate calculations will be different with numpy arrays -# For example: -# Bx = B[lat, 0] -# Will need to be changed to: -# Bx = B[lat, 0].reshape(-1, 1) -# (if we want it to have the same shape as before) - -# For reference, here is a list of the correct shapes of these objects: -# A: (6, 6) -# B: (6, 2) -# C: (2, 6) -# D: (2, 2) -# xd: (6, 1) -# yd: (6, 1) -# Ax: (4, 4) -# Bx: (4, 1) -# Cx: (1, 4) -# Dx: () -# Ay: (2, 2) -# By: (2, 1) -# Cy: (1, 2) +# (T, Y) = step_response(H1a, T=np.linspace(0,10,100)); +# # Step response for the first input -H1ax = ss(Ax - Bx*K1a[0, lat], Bx*K1a[0, lat]*xd[lat, :], Cx, Dx) -Yx, Tx = step(H1ax, T=np.linspace(0, 10, 100)) +H1ax = ct.ss(Ax - Bx @ K1a[np.ix_([0], lat)], + Bx @ K1a[np.ix_([0], lat)] @ xd[lat, :], Cx, Dx) +Tx, Yx = ct.step_response(H1ax, T=np.linspace(0, 10, 100)) # Step response for the second input -H1ay = ss(Ay - By*K1a[1, alt], By*K1a[1, alt]*yd[alt, :], Cy, Dy) -Yy, Ty = step(H1ay, T=np.linspace(0, 10, 100)) +H1ay = ct.ss(Ay - By @ K1a[np.ix_([1], alt)], + By @ K1a[np.ix_([1], alt)] @ yd[alt, :], Cy, Dy) +Ty, Yy = ct.step_response(H1ay, T=np.linspace(0, 10, 100)) plt.subplot(221) plt.title("Identity weights") @@ -164,20 +143,23 @@ # Look at different input weightings Qu1a = np.diag([1, 1]) -K1a, X, E = lqr(A, B, Qx1, Qu1a) -H1ax = ss(Ax - Bx*K1a[0, lat], Bx*K1a[0, lat]*xd[lat, :], Cx, Dx) +K1a, X, E = ct.lqr(A, B, Qx1, Qu1a) +H1ax = ct.ss(Ax - Bx @ K1a[np.ix_([0], lat)], + Bx @ K1a[np.ix_([0], lat)] @ xd[lat, :], Cx, Dx) Qu1b = (40 ** 2)*np.diag([1, 1]) -K1b, X, E = lqr(A, B, Qx1, Qu1b) -H1bx = ss(Ax - Bx*K1b[0, lat], Bx*K1b[0, lat]*xd[lat, :], Cx, Dx) +K1b, X, E = ct.lqr(A, B, Qx1, Qu1b) +H1bx = ct.ss(Ax - Bx @ K1b[np.ix_([0], lat)], + Bx @ K1b[np.ix_([0], lat)] @ xd[lat, :], Cx, Dx) Qu1c = (200 ** 2)*np.diag([1, 1]) -K1c, X, E = lqr(A, B, Qx1, Qu1c) -H1cx = ss(Ax - Bx*K1c[0, lat], Bx*K1c[0, lat]*xd[lat, :], Cx, Dx) +K1c, X, E = ct.lqr(A, B, Qx1, Qu1c) +H1cx = ct.ss(Ax - Bx @ K1c[np.ix_([0], lat)], + Bx @ K1c[np.ix_([0], lat)] @ xd[lat, :], Cx, Dx) -[Y1, T1] = step(H1ax, T=np.linspace(0, 10, 100)) -[Y2, T2] = step(H1bx, T=np.linspace(0, 10, 100)) -[Y3, T3] = step(H1cx, T=np.linspace(0, 10, 100)) +T1, Y1 = ct.step_response(H1ax, T=np.linspace(0, 10, 100)) +T2, Y2 = ct.step_response(H1bx, T=np.linspace(0, 10, 100)) +T3, Y3 = ct.step_response(H1cx, T=np.linspace(0, 10, 100)) plt.subplot(222) plt.title("Effect of input weights") @@ -189,21 +171,22 @@ plt.axis([0, 10, -0.1, 1.4]) # arcarrow([1.3, 0.8], [5, 0.45], -6) -plt.text(5.3, 0.4, 'rho') +plt.text(5.3, 0.4, r'$\rho$') # Output weighting - change Qx to use outputs -Qx2 = C.T*C -Qu2 = 0.1*np.diag([1, 1]) -K, X, E = lqr(A, B, Qx2, Qu2) -K2 = np.matrix(K) +Qx2 = C.T @ C +Qu2 = 0.1 * np.diag([1, 1]) +K2, X, E = ct.lqr(A, B, Qx2, Qu2) -H2x = ss(Ax - Bx*K2[0, lat], Bx*K2[0, lat]*xd[lat, :], Cx, Dx) -H2y = ss(Ay - By*K2[1, alt], By*K2[1, alt]*yd[alt, :], Cy, Dy) +H2x = ct.ss(Ax - Bx @ K2[np.ix_([0], lat)], + Bx @ K2[np.ix_([0], lat)] @ xd[lat, :], Cx, Dx) +H2y = ct.ss(Ay - By @ K2[np.ix_([1], alt)], + By @ K2[np.ix_([1], alt)] @ yd[alt, :], Cy, Dy) plt.subplot(223) plt.title("Output weighting") -[Y2x, T2x] = step(H2x, T=np.linspace(0, 10, 100)) -[Y2y, T2y] = step(H2y, T=np.linspace(0, 10, 100)) +T2x, Y2x = ct.step_response(H2x, T=np.linspace(0, 10, 100)) +T2y, Y2y = ct.step_response(H2y, T=np.linspace(0, 10, 100)) plt.plot(T2x.T, Y2x.T, T2y.T, Y2y.T) plt.ylabel('position') plt.xlabel('time') @@ -220,19 +203,21 @@ Qx3 = np.diag([100, 10, 2*np.pi/5, 0, 0, 0]) Qu3 = 0.1*np.diag([1, 10]) -(K, X, E) = lqr(A, B, Qx3, Qu3) -K3 = np.matrix(K) +K3, X, E = ct.lqr(A, B, Qx3, Qu3) -H3x = ss(Ax - Bx*K3[0, lat], Bx*K3[0, lat]*xd[lat, :], Cx, Dx) -H3y = ss(Ay - By*K3[1, alt], By*K3[1, alt]*yd[alt, :], Cy, Dy) +H3x = ct.ss(Ax - Bx @ K3[np.ix_([0], lat)], + Bx @ K3[np.ix_([0], lat)] @ xd[lat, :], Cx, Dx) +H3y = ct.ss(Ay - By @ K3[np.ix_([1], alt)], + By @ K3[np.ix_([1], alt)] @ yd[alt, :], Cy, Dy) plt.subplot(224) -# step(H3x, H3y, 10) -[Y3x, T3x] = step(H3x, T=np.linspace(0, 10, 100)) -[Y3y, T3y] = step(H3y, T=np.linspace(0, 10, 100)) +# step_response(H3x, H3y, 10) +T3x, Y3x = ct.step_response(H3x, T=np.linspace(0, 10, 100)) +T3y, Y3y = ct.step_response(H3y, T=np.linspace(0, 10, 100)) plt.plot(T3x.T, Y3x.T, T3y.T, Y3y.T) plt.title("Physically motivated weights") plt.xlabel('time') plt.legend(('x', 'y'), loc='lower right') +plt.tight_layout() if 'PYCONTROL_TEST_EXAMPLES' not in os.environ: plt.show() diff --git a/examples/pvtol-nested.py b/examples/pvtol-nested.py index 24cd7d1c5..040b4a1f4 100644 --- a/examples/pvtol-nested.py +++ b/examples/pvtol-nested.py @@ -9,8 +9,8 @@ # import os -import matplotlib.pyplot as plt # MATLAB plotting functions -from control.matlab import * # MATLAB-like functions +import matplotlib.pyplot as plt # MATLAB-like plotting functions +import control as ct import numpy as np # System parameters @@ -21,8 +21,8 @@ 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) # # Inner loop control design @@ -34,59 +34,58 @@ # Design a simple lead controller for the system k, a, b = 200, 2, 50 -Ci = k*tf([1, a], [1, b]) # lead compensator -Li = Pi*Ci +Ci = k * ct.tf([1, a], [1, b]) # lead compensator +Li = Pi * Ci # Bode plot for the open loop process plt.figure(1) -bode(Pi) +ct.bode_plot(Pi) # Bode plot for the loop transfer function, with margins plt.figure(2) -bode(Li) +ct.bode_plot(Li) # Compute out the gain and phase margins -#! Not implemented -# gm, pm, wcg, wcp = margin(Li) +gm, pm, wcg, wcp = ct.margin(Li) # Compute the sensitivity and complementary sensitivity functions -Si = feedback(1, Li) -Ti = Li*Si +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() plt.subplot(221) -bode(Hi) +ct.bode_plot(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 Lo = -m*g*Po*Co plt.figure(5) -bode(Lo) # margin(Lo) +ct.bode_plot(Lo) # margin(Lo) # Finally compute the real outer-loop loop gain + responses -L = Co*Hi*Po -S = feedback(1, L) -T = feedback(L, 1) +L = Co * Hi * Po +S = ct.feedback(1, L) +T = ct.feedback(L, 1) # Compute stability margins -gm, pm, wgc, wpc = margin(L) +gm, pm, wgc, wpc = ct.margin(L) print("Gain margin: %g at %g" % (gm, wgc)) print("Phase margin: %g at %g" % (pm, wpc)) plt.figure(6) plt.clf() -bode(L, np.logspace(-4, 3)) +ct.bode_plot(L, np.logspace(-4, 3)) # Add crossover line to the magnitude plot # @@ -113,7 +112,7 @@ break # Recreate the frequency response and shift the phase -mag, phase, w = freqresp(L, np.logspace(-4, 3)) +mag, phase, w = ct.freqresp(L, np.logspace(-4, 3)) phase = phase - 360 # Replot the phase by hand @@ -130,7 +129,7 @@ # plt.figure(7) plt.clf() -nyquist(L, (0.0001, 1000)) +ct.nyquist_plot(L, (0.0001, 1000)) # Add a box in the region we are going to expand plt.plot([-2, -2, 1, 1, -2], [-4, 4, 4, -4, -4], 'r-') @@ -138,7 +137,7 @@ # Expanded region plt.figure(8) plt.clf() -nyquist(L) +ct.nyquist_plot(L) plt.axis([-2, 1, -4, 4]) # set up the color @@ -154,21 +153,21 @@ # 'EdgeColor', color, 'FaceColor', color); plt.figure(9) -Yvec, Tvec = step(T, np.linspace(0, 20)) +Tvec, Yvec = ct.step_response(T, np.linspace(0, 20)) plt.plot(Tvec.T, Yvec.T) -Yvec, Tvec = step(Co*S, np.linspace(0, 20)) +Tvec, Yvec = ct.step_response(Co*S, np.linspace(0, 20)) plt.plot(Tvec.T, Yvec.T) plt.figure(10) plt.clf() -P, Z = pzmap(T, plot=True, grid=True) +P, Z = ct.pzmap(T, plot=True, grid=True) print("Closed loop poles and zeros: ", P, Z) # Gang of Four plt.figure(11) plt.clf() -gangof4(Hi*Po, Co) +ct.gangof4_plot(Hi * Po, Co) if 'PYCONTROL_TEST_EXAMPLES' not in os.environ: plt.show() diff --git a/examples/steering-gainsched.py b/examples/steering-gainsched.py index 7db2d9a73..7ddc6b5b8 100644 --- a/examples/steering-gainsched.py +++ b/examples/steering-gainsched.py @@ -10,7 +10,7 @@ import numpy as np import control as ct from cmath import sqrt -import matplotlib.pyplot as mpl +import matplotlib.pyplot as plt # # Vehicle steering dynamics @@ -137,7 +137,7 @@ def trajgen_output(t, x, u, params): # We construct the system using the InterconnectedSystem constructor and using # signal labels to keep track of everything. -steering = ct.InterconnectedSystem( +steering = ct.interconnect( # List of subsystems (trajgen, controller, vehicle), name='steering', @@ -167,10 +167,10 @@ def trajgen_output(t, x, u, params): T = np.linspace(0, 5, 100) # Set up a figure for plotting the results -mpl.figure(); +plt.figure(); # Plot the reference trajectory for the y position -mpl.plot([0, 5], [yref, yref], 'k--') +plt.plot([0, 5], [yref, yref], 'k--') # Find the signals we want to plot y_index = steering.find_output('y') @@ -183,13 +183,13 @@ def trajgen_output(t, x, u, params): steering, T, [vref * np.ones(len(T)), yref * np.ones(len(T))]) # Plot the reference speed - mpl.plot([0, 5], [vref, vref], 'k--') + plt.plot([0, 5], [vref, vref], 'k--') # Plot the system output - y_line, = mpl.plot(tout, yout[y_index, :], 'r') # lateral position - v_line, = mpl.plot(tout, yout[v_index, :], 'b') # vehicle velocity + y_line, = plt.plot(tout, yout[y_index, :], 'r') # lateral position + v_line, = plt.plot(tout, yout[v_index, :], 'b') # vehicle velocity # Add axis labels -mpl.xlabel('Time (s)') -mpl.ylabel('x vel (m/s), y pos (m)') -mpl.legend((v_line, y_line), ('v', 'y'), loc='center right', frameon=False) +plt.xlabel('Time (s)') +plt.ylabel('x vel (m/s), y pos (m)') +plt.legend((v_line, y_line), ('v', 'y'), loc='center right', frameon=False) diff --git a/examples/tfvis.py b/examples/tfvis.py index 30a084ffb..0cb789db4 100644 --- a/examples/tfvis.py +++ b/examples/tfvis.py @@ -270,8 +270,8 @@ def button_release(self, event): tfcn = self.tfi.get_tf() if (tfcn): - self.zeros = tfcn.zero() - self.poles = tfcn.pole() + self.zeros = tfcn.zeros() + self.poles = tfcn.poles() self.sys = tfcn self.redraw() @@ -314,8 +314,8 @@ def apply(self): tfcn = self.tfi.get_tf() if (tfcn): - self.zeros = tfcn.zero() - self.poles = tfcn.pole() + self.zeros = tfcn.zeros() + self.poles = tfcn.poles() self.sys = tfcn self.redraw()