changes so that a default dt=0 passes unit tests. fixed code everywhere that combines systems with different timebases

sawyerbfuller · bnavigator · commit 1ed9ee272d63 · 2020-12-31T01:14:00.000+01:00
diff --git a/control/config.py b/control/config.py
@@ -59,6 +59,9 @@ def reset_defaults():
     from .statesp import _statesp_defaults
     defaults.update(_statesp_defaults)
 
+    from .iosys import _iosys_defaults
+    defaults.update(_iosys_defaults)
+
 
 def _get_param(module, param, argval=None, defval=None, pop=False):
     """Return the default value for a configuration option.
@@ -208,8 +211,13 @@ def use_legacy_defaults(version):
     # Go backwards through releases and reset defaults
     #
 
-    # Version 0.9.0: switched to 'array' as default for state space objects
+    # Version 0.9.0:
     if major == 0 and minor < 9:
+        # switched to 'array' as default for state space objects
         set_defaults('statesp', use_numpy_matrix=True)
+        # switched to 0 (=continuous) as default timestep
+        set_defaults('statesp', default_dt=None)
+        set_defaults('xferfcn', default_dt=None)
+        set_defaults('iosys', default_dt=None)
 
     return (major, minor, patch)
diff --git a/control/iosys.py b/control/iosys.py
@@ -38,12 +38,16 @@
 
 from .statesp import StateSpace, tf2ss
 from .timeresp import _check_convert_array
-from .lti import isctime, isdtime, _find_timebase
+from .lti import isctime, isdtime, common_timebase
+from . import config
 
 __all__ = ['InputOutputSystem', 'LinearIOSystem', 'NonlinearIOSystem',
            'InterconnectedSystem', 'input_output_response', 'find_eqpt',
            'linearize', 'ss2io', 'tf2io']
 
+# Define module default parameter values
+_iosys_defaults = {
+    'iosys.default_dt': 0}
 
 class InputOutputSystem(object):
     """A class for representing input/output systems.
@@ -118,7 +122,7 @@ def name_or_default(self, name=None):
         return name
 
     def __init__(self, inputs=None, outputs=None, states=None, params={},
-                 dt=None, name=None):
+                 name=None, **kwargs):
         """Create an input/output system.
 
         The InputOutputSystem contructor is used to create an input/output
@@ -163,7 +167,7 @@ def __init__(self, inputs=None, outputs=None, states=None, params={},
         """
         # Store the input arguments
         self.params = params.copy()             # default parameters
-        self.dt = dt                            # timebase
+        self.dt = kwargs.get('dt', config.defaults['iosys.default_dt'])                    # timebase
         self.name = self.name_or_default(name)  # system name
 
         # Parse and store the number of inputs, outputs, and states
@@ -210,9 +214,7 @@ def __mul__(sys2, sys1):
                              "inputs and outputs")
 
         # Make sure timebase are compatible
-        dt = _find_timebase(sys1, sys2)
-        if dt is False:
-            raise ValueError("System timebases are not compabile")
+        dt = common_timebase(sys1.dt, sys2.dt)
 
         inplist = [(0,i) for i in range(sys1.ninputs)]
         outlist = [(1,i) for i in range(sys2.noutputs)]
@@ -464,12 +466,11 @@ def feedback(self, other=1, sign=-1, params={}):
                              "inputs and outputs")
 
         # Make sure timebases are compatible
-        dt = _find_timebase(self, other)
-        if dt is False:
-            raise ValueError("System timebases are not compabile")
+        dt = common_timebase(self.dt, other.dt)
 
         inplist = [(0,i) for i in range(self.ninputs)]
         outlist = [(0,i) for i in range(self.noutputs)]
+
         # Return the series interconnection between the systems
         newsys = InterconnectedSystem((self, other), inplist=inplist, outlist=outlist,
                                       params=params, dt=dt)
@@ -650,7 +651,8 @@ class NonlinearIOSystem(InputOutputSystem):
 
     """
     def __init__(self, updfcn, outfcn=None, inputs=None, outputs=None,
-                 states=None, params={}, dt=None, name=None):
+                 states=None, params={},
+                 name=None, **kwargs):
         """Create a nonlinear I/O system given update and output functions.
 
         Creates an `InputOutputSystem` for a nonlinear system by specifying a
@@ -722,6 +724,7 @@ def __init__(self, updfcn, outfcn=None, inputs=None, outputs=None,
         self.outfcn = outfcn
 
         # Initialize the rest of the structure
+        dt = kwargs.get('dt', config.defaults['iosys.default_dt'])
         super(NonlinearIOSystem, self).__init__(
             inputs=inputs, outputs=outputs, states=states,
             params=params, dt=dt, name=name
@@ -888,20 +891,14 @@ def __init__(self, syslist, connections=[], inplist=[], outlist=[],
         # Check to make sure all systems are consistent
         self.syslist = syslist
         self.syslist_index = {}
-        dt = None
         nstates = 0; self.state_offset = []
         ninputs = 0; self.input_offset = []
         noutputs = 0; self.output_offset = []
         sysobj_name_dct = {}
         sysname_count_dct = {}
         for sysidx, sys in enumerate(syslist):
             # Make sure time bases are consistent
-            # TODO: Use lti._find_timebase() instead?
-            if dt is None and sys.dt is not None:
-                # Timebase was not specified; set to match this system
-                dt = sys.dt
-            elif dt != sys.dt:
-                raise TypeError("System timebases are not compatible")
+            dt = common_timebase(dt, sys.dt)
 
             # Make sure number of inputs, outputs, states is given
             if sys.ninputs is None or sys.noutputs is None or \
diff --git a/control/lti.py b/control/lti.py
@@ -9,13 +9,13 @@
 isdtime()
 isctime()
 timebase()
-timebaseEqual()
+common_timebase()
 """
 
 import numpy as np
 from numpy import absolute, real
 
-__all__ = ['issiso', 'timebase', 'timebaseEqual', 'isdtime', 'isctime',
+__all__ = ['issiso', 'timebase', 'common_timebase', 'isdtime', 'isctime',
            'pole', 'zero', 'damp', 'evalfr', 'freqresp', 'dcgain']
 
 class LTI:
@@ -157,48 +157,31 @@ def timebase(sys, strict=True):
 
     return sys.dt
 
-# 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.
-    """
-
-    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
-
-# Find a common timebase between two or more systems
-def _find_timebase(sys1, *sysn):
-    """Find the common timebase between systems, otherwise return False"""
-
-    # Create a list of systems to check
-    syslist = [sys1]
-    syslist.append(*sysn)
-
-    # Look for a common timebase
-    dt = None
-
-    for sys in syslist:
-        # Make sure time bases are consistent
-        if (dt is None and sys.dt is not None) or \
-           (dt is True and isdiscrete(sys)):
-            # Timebase was not specified; set to match this system
-            dt = sys.dt
-        elif dt != sys.dt:
-            return False
-    return dt
-
+def common_timebase(dt1, dt2):
+    """Find the common timebase when interconnecting systems."""
+    # cases: 
+    # 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 they must be equal (holds for both cont and discrete systems)
+    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 a system is a discrete time system
 def isdtime(sys, strict=False):
diff --git a/control/statesp.py b/control/statesp.py
@@ -62,7 +62,7 @@
 from scipy.signal import cont2discrete
 from scipy.signal import StateSpace as signalStateSpace
 from warnings import warn
-from .lti import LTI, timebase, timebaseEqual, isdtime
+from .lti import LTI, common_timebase, isdtime
 from . import config
 from copy import deepcopy
 
@@ -180,7 +180,10 @@ def __init__(self, *args, **kw):
         if len(args) == 4:
             # The user provided A, B, C, and D matrices.
             (A, B, C, D) = args
-            dt = config.defaults['statesp.default_dt']
+            if _isstaticgain(A, B, C, D): 
+                dt = None
+            else:
+                dt = config.defaults['statesp.default_dt']
         elif len(args) == 5:
             # Discrete time system
             (A, B, C, D, dt) = args
@@ -196,9 +199,12 @@ def __init__(self, *args, **kw):
             try:
                 dt = args[0].dt
             except NameError:
-                dt = config.defaults['statesp.default_dt']
+                if _isstaticgain(A, B, C, D): 
+                    dt = None
+                else:
+                    dt = config.defaults['statesp.default_dt']
         else:
-            raise ValueError("Needs 1 or 4 arguments; received %i." % len(args))
+            raise ValueError("Expected 1, 4, or 5 arguments; received %i." % len(args))
 
         # Process keyword arguments
         remove_useless = kw.get('remove_useless',
@@ -330,14 +336,7 @@ def __add__(self, other):
                     (self.outputs != other.outputs)):
                 raise ValueError("Systems have different shapes.")
 
-            # Figure out the sampling time to use
-            if self.dt is None and other.dt is not None:
-                dt = other.dt       # use dt from second argument
-            elif (other.dt is None and self.dt is not None) or \
-                    (timebaseEqual(self, other)):
-                dt = self.dt        # use dt from first argument
-            else:
-                raise ValueError("Systems have different sampling times")
+            dt = common_timebase(self.dt, other.dt)
 
             # Concatenate the various arrays
             A = concatenate((
@@ -386,16 +385,8 @@ def __mul__(self, other):
             # Check to make sure the dimensions are OK
             if self.inputs != other.outputs:
                 raise ValueError("C = A * B: A has %i column(s) (input(s)), \
-but B has %i row(s)\n(output(s))." % (self.inputs, other.outputs))
-
-            # Figure out the sampling time to use
-            if (self.dt == None and other.dt != None):
-                dt = other.dt       # use dt from second argument
-            elif (other.dt == None and self.dt != None) or \
-                    (timebaseEqual(self, other)):
-                dt = self.dt        # use dt from first argument
-            else:
-                raise ValueError("Systems have different sampling times")
+                    but B has %i row(s)\n(output(s))." % (self.inputs, other.outputs))
+            dt = common_timebase(self.dt, other.dt)
 
             # Concatenate the various arrays
             A = concatenate(
@@ -467,9 +458,8 @@ def _evalfr(self, omega):
         """Evaluate a SS system's transfer function at a single frequency"""
         # Figure out the point to evaluate the transfer function
         if isdtime(self, strict=True):
-            dt = timebase(self)
-            s = exp(1.j * omega * dt)
-            if omega * dt > math.pi:
+            s = exp(1.j * omega * self.dt)
+            if omega * self.dt > math.pi:
                 warn("_evalfr: frequency evaluation above Nyquist frequency")
         else:
             s = omega * 1.j
@@ -526,9 +516,8 @@ def freqresp(self, omega):
         # axis (continuous time) or unit circle (discrete time).
         omega.sort()
         if isdtime(self, strict=True):
-            dt = timebase(self)
-            cmplx_freqs = exp(1.j * omega * dt)
-            if max(np.abs(omega)) * dt > math.pi:
+            cmplx_freqs = exp(1.j * omega * self.dt)
+            if max(np.abs(omega)) * self.dt > math.pi:
                 warn("freqresp: frequency evaluation above Nyquist frequency")
         else:
             cmplx_freqs = omega * 1.j
@@ -631,14 +620,7 @@ def feedback(self, other=1, sign=-1):
         if (self.inputs != other.outputs) or (self.outputs != other.inputs):
                 raise ValueError("State space systems don't have compatible inputs/outputs for "
                                  "feedback.")
-
-        # Figure out the sampling time to use
-        if self.dt is None and other.dt is not None:
-            dt = other.dt       # use dt from second argument
-        elif other.dt is None and self.dt is not None or timebaseEqual(self, other):
-            dt = self.dt        # use dt from first argument
-        else:
-            raise ValueError("Systems have different sampling times")
+        dt = common_timebase(self.dt, other.dt)
 
         A1 = self.A
         B1 = self.B
@@ -708,14 +690,7 @@ def lft(self, other, nu=-1, ny=-1):
         # dimension check
         # TODO
 
-        # Figure out the sampling time to use
-        if (self.dt == None and other.dt != None):
-            dt = other.dt       # use dt from second argument
-        elif (other.dt == None and self.dt != None) or \
-                timebaseEqual(self, other):
-            dt = self.dt        # use dt from first argument
-        else:
-            raise ValueError("Systems have different time bases")
+        dt = common_timebase(self.dt, other.dt)
 
         # submatrices
         A = self.A
@@ -858,8 +833,7 @@ def append(self, other):
         if not isinstance(other, StateSpace):
             other = _convertToStateSpace(other)
 
-        if self.dt != other.dt:
-            raise ValueError("Systems must have the same time step")
+        self.dt = common_timebase(self.dt, other.dt)
 
         n = self.states + other.states
         m = self.inputs + other.inputs
@@ -1293,6 +1267,11 @@ def _mimo2simo(sys, input, warn_conversion=False):
 
     return sys
 
+def _isstaticgain(A, B, C, D):
+    """returns True if and only if the system has no dynamics, that is, 
+    if A and B are zero. """
+    return not np.any(np.matrix(A, dtype=float)) \
+        and not np.any(np.matrix(B, dtype=float))
 
 def ss(*args):
     """ss(A, B, C, D[, dt])
diff --git a/control/tests/config_test.py b/control/tests/config_test.py
@@ -246,7 +246,8 @@ def test_change_default_dt(self, dt):
 
     def test_change_default_dt_static(self):
         """Test that static gain systems always have dt=None"""
-        ct.set_defaults('control', default_dt=0)
+        ct.set_defaults('xferfcn', default_dt=0)
         assert ct.tf(1, 1).dt is None
+        ct.set_defaults('statesp', default_dt=0)
         assert ct.ss(0, 0, 0, 1).dt is None
         # TODO: add in test for static gain iosys
diff --git a/control/tests/matlab_test.py b/control/tests/matlab_test.py
diff --git a/control/xferfcn.py b/control/xferfcn.py
