Merge pull request #514 from murrayrm/mimo_impulse_step_response

MIMO impulse and step response

control/config.py

@@ -17,7 +17,6 @@
 _control_defaults = {
     'control.default_dt': 0,
     'control.squeeze_frequency_response': None,
-    'control.squeeze_time_response': True,
     'control.squeeze_time_response': None,
     'forced_response.return_x': False,
 }

control/iosys.py

@@ -1428,8 +1428,9 @@ def input_output_response(sys, T, U=0., X0=0, params={}, method='RK45',
         for i in range(len(T)):
             u = U[i] if len(U.shape) == 1 else U[:, i]
             y[:, i] = sys._out(T[i], [], u)
-        return _process_time_response(sys, T, y, [], transpose=transpose,
-                                      return_x=return_x, squeeze=squeeze)
+        return _process_time_response(
+            sys, T, y, np.array((0, 0, np.asarray(T).size)),
+            transpose=transpose, return_x=return_x, squeeze=squeeze)
 
     # create X0 if not given, test if X0 has correct shape
     X0 = _check_convert_array(X0, [(nstates,), (nstates, 1)],

control/tests/matlab2_test.py

@@ -121,25 +121,25 @@ def test_step(self, SISO_mats, MIMO_mats, mplcleanup):
         #print("gain:", dcgain(sys))
 
         subplot2grid(plot_shape, (0, 0))
-        t, y = step(sys)
+        y, t = step(sys)
         plot(t, y)
 
         subplot2grid(plot_shape, (0, 1))
         T = linspace(0, 2, 100)
         X0 = array([1, 1])
-        t, y = step(sys, T, X0)
+        y, t = step(sys, T, X0)
         plot(t, y)
 
         # Test output of state vector
-        t, y, x = step(sys, return_x=True)
+        y, t, x = step(sys, return_x=True)
 
         #Test MIMO system
         A, B, C, D = MIMO_mats
         sys = ss(A, B, C, D)
 
         subplot2grid(plot_shape, (0, 2))
-        t, y = step(sys)
-        plot(t, y)
+        y, t = step(sys)
+        plot(t, y[:, 0, 0])
 
     def test_impulse(self, SISO_mats, mplcleanup):
         A, B, C, D = SISO_mats
@@ -168,8 +168,8 @@ def test_impulse_mimo(self, MIMO_mats, mplcleanup):
         #Test MIMO system
         A, B, C, D = MIMO_mats
         sys = ss(A, B, C, D)
-        t, y = impulse(sys)
-        plot(t, y, label='MIMO System')
+        y, t = impulse(sys)
+        plot(t, y[:, :, 0], label='MIMO System')
 
         legend(loc='best')
         #show()

control/tests/timeresp_test.py

@@ -347,7 +347,7 @@ def test_impulse_response_mimo(self, mimo_ss2):
         yref_notrim = np.zeros((2, len(t)))
         yref_notrim[:1, :] = yref
         _t, yy = impulse_response(sys, T=t, input=0)
-        np.testing.assert_array_almost_equal(yy, yref_notrim, decimal=4)
+        np.testing.assert_array_almost_equal(yy[:,0,:], yref_notrim, decimal=4)
 
     @pytest.mark.skipif(StrictVersion(sp.__version__) < "1.3",
                         reason="requires SciPy 1.3 or greater")
@@ -639,9 +639,10 @@ def test_time_vector(self, tsystem, fun, squeeze, matarrayout):
         if hasattr(tsystem, 't'):
             # tout should always match t, which has shape (n, )
             np.testing.assert_allclose(tout, tsystem.t)
-        if squeeze is False or sys.outputs > 1:
+
+        if squeeze is False or not sys.issiso():
             assert yout.shape[0] == sys.outputs
-            assert yout.shape[1] == tout.shape[0]
+            assert yout.shape[-1] == tout.shape[0]
         else:
             assert yout.shape == tout.shape
 
@@ -725,21 +726,22 @@ def test_time_series_data_convention_2D(self, siso_ss1):
 
     @pytest.mark.usefixtures("editsdefaults")
     @pytest.mark.parametrize("fcn", [ct.ss, ct.tf, ct.ss2io])
-    @pytest.mark.parametrize("nstate, nout, ninp, squeeze, shape", [
-        [1, 1, 1, None, (8,)],
-        [2, 1, 1, True, (8,)],
-        [3, 1, 1, False, (1, 8)],
-        [3, 2, 1, None, (2, 8)],
-        [4, 2, 1, True, (2, 8)],
-        [5, 2, 1, False, (2, 8)],
-        [3, 1, 2, None, (1, 8)],
-        [4, 1, 2, True, (8,)],
-        [5, 1, 2, False, (1, 8)],
-        [4, 2, 2, None, (2, 8)],
-        [5, 2, 2, True, (2, 8)],
-        [6, 2, 2, False, (2, 8)],
+    @pytest.mark.parametrize("nstate, nout, ninp, squeeze, shape1, shape2", [
+    #  state  out   in   squeeze  in/out      out-only
+        [1,    1,    1,  None,   (8,),       (8,)],
+        [2,    1,    1,  True,   (8,),       (8,)],
+        [3,    1,    1,  False,  (1, 1, 8),  (1, 8)],
+        [3,    2,    1,  None,   (2, 1, 8),  (2, 8)],
+        [4,    2,    1,  True,   (2, 8),     (2, 8)],
+        [5,    2,    1,  False,  (2, 1, 8),  (2, 8)],
+        [3,    1,    2,  None,   (1, 2, 8),  (1, 8)],
+        [4,    1,    2,  True,   (2, 8),     (8,)],
+        [5,    1,    2,  False,  (1, 2, 8),  (1, 8)],
+        [4,    2,    2,  None,   (2, 2, 8),  (2, 8)],
+        [5,    2,    2,  True,   (2, 2, 8),  (2, 8)],
+        [6,    2,    2,  False,  (2, 2, 8),  (2, 8)],
     ])
-    def test_squeeze(self, fcn, nstate, nout, ninp, squeeze, shape):
+    def test_squeeze(self, fcn, nstate, nout, ninp, squeeze, shape1, shape2):
         # Figure out if we have SciPy 1+
         scipy0 = StrictVersion(sp.__version__) < '1.0'
 
@@ -750,27 +752,56 @@ def test_squeeze(self, fcn, nstate, nout, ninp, squeeze, shape):
         else:
             sys = fcn(ct.rss(nstate, nout, ninp, strictly_proper=True))
 
-        # Keep track of expect users warnings
-        warntype = UserWarning if sys.inputs > 1 else None
-
         # Generate the time and input vectors
         tvec = np.linspace(0, 1, 8)
         uvec = np.dot(
             np.ones((sys.inputs, 1)),
             np.reshape(np.sin(tvec), (1, 8)))
 
+        #
         # Pass squeeze argument and make sure the shape is correct
-        with pytest.warns(warntype, match="Converting MIMO system"):
-            _, yvec = ct.impulse_response(sys, tvec, squeeze=squeeze)
-        assert yvec.shape == shape
+        #
+        # For responses that are indexed by the input, check against shape1
+        # For responses that have no/fixed input, check against shape2
+        #
 
-        _, yvec = ct.initial_response(sys, tvec, 1, squeeze=squeeze)
-        assert yvec.shape == shape
+        # Impulse response
+        if isinstance(sys, StateSpace):
+            # Check the states as well
+            _, yvec, xvec = ct.impulse_response(
+                sys, tvec, squeeze=squeeze, return_x=True)
+            if sys.issiso():
+                assert xvec.shape == (sys.states, 8)
+            else:
+                assert xvec.shape == (sys.states, sys.inputs, 8)
+        else:
+            _, yvec = ct.impulse_response(sys, tvec, squeeze=squeeze)
+        assert yvec.shape == shape1
 
-        with pytest.warns(warntype, match="Converting MIMO system"):
+        # Step response
+        if isinstance(sys, StateSpace):
+            # Check the states as well
+            _, yvec, xvec = ct.step_response(
+                sys, tvec, squeeze=squeeze, return_x=True)
+            if sys.issiso():
+                assert xvec.shape == (sys.states, 8)
+            else:
+                assert xvec.shape == (sys.states, sys.inputs, 8)
+        else:
             _, yvec = ct.step_response(sys, tvec, squeeze=squeeze)
-        assert yvec.shape == shape
+        assert yvec.shape == shape1
+
+        # Initial response (only indexed by output)
+        if isinstance(sys, StateSpace):
+            # Check the states as well
+            _, yvec, xvec = ct.initial_response(
+                sys, tvec, 1, squeeze=squeeze, return_x=True)
+            assert xvec.shape == (sys.states, 8)
+        else:
+            _, yvec = ct.initial_response(sys, tvec, 1, squeeze=squeeze)
+        assert yvec.shape == shape2
 
+        # Forced response (only indexed by output)
         if isinstance(sys, StateSpace):
             # Check the states as well
             _, yvec, xvec = ct.forced_response(
@@ -779,52 +810,54 @@ def test_squeeze(self, fcn, nstate, nout, ninp, squeeze, shape):
         else:
             # Just check the input/output response
             _, yvec = ct.forced_response(sys, tvec, uvec, 0, squeeze=squeeze)
-        assert yvec.shape == shape
+        assert yvec.shape == shape2
 
         # Test cases where we choose a subset of inputs and outputs
         _, yvec = ct.step_response(
             sys, tvec, input=ninp-1, output=nout-1, squeeze=squeeze)
-        # Possible code if we implemenet a squeeze='siso' option
         if squeeze is False:
             # Shape should be unsqueezed
-            assert yvec.shape == (1, 8)
+            assert yvec.shape == (1, 1, 8)
         else:
             # Shape should be squeezed
             assert yvec.shape == (8, )
 
-        # For InputOutputSystems, also test input_output_response
+        # For InputOutputSystems, also test input/output response
         if isinstance(sys, ct.InputOutputSystem) and not scipy0:
             _, yvec = ct.input_output_response(sys, tvec, uvec, squeeze=squeeze)
-            assert yvec.shape == shape
+            assert yvec.shape == shape2
 
         #
         # Changing config.default to False should return 3D frequency response
         #
         ct.config.set_defaults('control', squeeze_time_response=False)
 
-        with pytest.warns(warntype, match="Converting MIMO system"):
-            _, yvec = ct.impulse_response(sys, tvec)
-        assert yvec.shape == (sys.outputs, 8)
+        _, yvec = ct.impulse_response(sys, tvec)
+        if squeeze is not True or sys.inputs > 1 or sys.outputs > 1:
+            assert yvec.shape == (sys.outputs, sys.inputs, 8)
 
-        _, yvec = ct.initial_response(sys, tvec, 1)
-        assert yvec.shape == (sys.outputs, 8)
+        _, yvec = ct.step_response(sys, tvec)
+        if squeeze is not True or sys.inputs > 1 or sys.outputs > 1:
+            assert yvec.shape == (sys.outputs, sys.inputs, 8)
 
-        with pytest.warns(warntype, match="Converting MIMO system"):
-            _, yvec = ct.step_response(sys, tvec)
-        assert yvec.shape == (sys.outputs, 8)
+        _, yvec = ct.initial_response(sys, tvec, 1)
+        if squeeze is not True or sys.outputs > 1:
+            assert yvec.shape == (sys.outputs, 8)
 
         if isinstance(sys, ct.StateSpace):
             _, yvec, xvec = ct.forced_response(
                 sys, tvec, uvec, 0, return_x=True)
             assert xvec.shape == (sys.states, 8)
         else:
             _, yvec = ct.forced_response(sys, tvec, uvec, 0)
-        assert yvec.shape == (sys.outputs, 8)
+        if squeeze is not True or sys.outputs > 1:
+            assert yvec.shape == (sys.outputs, 8)
 
         # For InputOutputSystems, also test input_output_response
         if isinstance(sys, ct.InputOutputSystem) and not scipy0:
             _, yvec = ct.input_output_response(sys, tvec, uvec)
-            assert yvec.shape == (sys.noutputs, 8)
+            if squeeze is not True or sys.outputs > 1:
+                assert yvec.shape == (sys.outputs, 8)
 
     @pytest.mark.parametrize("fcn", [ct.ss, ct.tf, ct.ss2io])
     def test_squeeze_exception(self, fcn):
@@ -861,3 +894,37 @@ def test_squeeze_0_8_4(self, nstate, nout, ninp, squeeze, shape):
 
         _, yvec = ct.initial_response(sys, tvec, 1, squeeze=squeeze)
         assert yvec.shape == shape
+
+    @pytest.mark.parametrize(
+        "nstate, nout, ninp, squeeze, ysh_in, ysh_no, xsh_in", [
+        [4,    1,    1,  None,   (8,),       (8,),    (8, 4)],
+        [4,    1,    1,  True,   (8,),       (8,),    (8, 4)],
+        [4,    1,    1,  False,  (8, 1, 1),  (8, 1),  (8, 4)],
+        [4,    2,    1,  None,   (8, 2, 1),  (8, 2),  (8, 4, 1)],
+        [4,    2,    1,  True,   (8, 2),     (8, 2),  (8, 4, 1)],
+        [4,    2,    1,  False,  (8, 2, 1),  (8, 2),  (8, 4, 1)],
+        [4,    1,    2,  None,   (8, 1, 2),  (8, 1),  (8, 4, 2)],
+        [4,    1,    2,  True,   (8, 2),     (8,),    (8, 4, 2)],
+        [4,    1,    2,  False,  (8, 1, 2),  (8, 1),  (8, 4, 2)],
+        [4,    2,    2,  None,   (8, 2, 2),  (8, 2),  (8, 4, 2)],
+        [4,    2,    2,  True,   (8, 2, 2),  (8, 2),  (8, 4, 2)],
+        [4,    2,    2,  False,  (8, 2, 2),  (8, 2),  (8, 4, 2)],
+    ])
+    def test_response_transpose(
+            self, nstate, nout, ninp, squeeze, ysh_in, ysh_no, xsh_in):
+        sys = ct.rss(nstate, nout, ninp)
+        T = np.linspace(0, 1, 8)
+
+        # Step response - input indexed
+        t, y, x = ct.step_response(
+            sys, T, transpose=True, return_x=True, squeeze=squeeze)
+        assert t.shape == (T.size, )
+        assert y.shape == ysh_in
+        assert x.shape == xsh_in
+
+        # Initial response - no input indexing
+        t, y, x = ct.initial_response(
+            sys, T, 1, transpose=True, return_x=True, squeeze=squeeze)
+        assert t.shape == (T.size, )
+        assert y.shape == ysh_no
+        assert x.shape == (T.size, sys.states)