From 8bb8e2219081a705e5376cd3d295f2620f353b01 Mon Sep 17 00:00:00 2001
From: Richard Murray <murray@cds.caltech.edu>
Date: Sat, 16 Jan 2021 15:50:00 -0800
Subject: [PATCH 1/3] standardize time response return values, return_x/squeeze
 keywords

---
 control/config.py               |   8 +-
 control/iosys.py                |  28 ++--
 control/matlab/timeresp.py      |  38 +++---
 control/tests/config_test.py    |  11 +-
 control/tests/discrete_test.py  |   8 +-
 control/tests/flatsys_test.py   |   2 +-
 control/tests/iosys_test.py     |  55 ++++----
 control/tests/modelsimp_test.py |   4 +-
 control/tests/timeresp_test.py  | 150 ++++++++++++++++++++-
 control/timeresp.py             | 231 +++++++++++++++++++-------------
 10 files changed, 357 insertions(+), 178 deletions(-)

diff --git a/control/config.py b/control/config.py
index 9ac953e11..6a8fb8db6 100644
--- a/control/config.py
+++ b/control/config.py
@@ -16,7 +16,9 @@
 # Package level default values
 _control_defaults = {
     'control.default_dt': 0,
-    'control.squeeze_frequency_response': None
+    'control.squeeze_frequency_response': None,
+    'control.squeeze_time_response': True,
+    'forced_response.return_x': False,
 }
 defaults = dict(_control_defaults)
 
@@ -211,6 +213,7 @@ def use_legacy_defaults(version):
     #
     # Go backwards through releases and reset defaults
     #
+    reset_defaults()            # start from a clean slate
 
     # Version 0.9.0:
     if major == 0 and minor < 9:
@@ -230,4 +233,7 @@ def use_legacy_defaults(version):
         # turned off _remove_useless_states
         set_defaults('statesp', remove_useless_states=True)
 
+        # forced_response no longer returns x by default
+        set_defaults('forced_response', return_x=True)
+
     return (major, minor, patch)
diff --git a/control/iosys.py b/control/iosys.py
index 913e8d471..8fc17c016 100644
--- a/control/iosys.py
+++ b/control/iosys.py
@@ -37,7 +37,7 @@
 from warnings import warn
 
 from .statesp import StateSpace, tf2ss
-from .timeresp import _check_convert_array
+from .timeresp import _check_convert_array, _process_time_response
 from .lti import isctime, isdtime, common_timebase
 from . import config
 
@@ -1353,7 +1353,7 @@ def __init__(self, io_sys, ss_sys=None):
 
 
 def input_output_response(sys, T, U=0., X0=0, params={}, method='RK45',
-                          return_x=False, squeeze=True):
+                          transpose=False, return_x=False, squeeze=None):
 
     """Compute the output response of a system to a given input.
 
@@ -1373,9 +1373,9 @@ def input_output_response(sys, T, U=0., X0=0, params={}, method='RK45',
     return_x : bool, optional
         If True, return the values of the state at each time (default = False).
     squeeze : bool, optional
-        If True (default), squeeze unused dimensions out of the output
-        response.  In particular, for a single output system, return a
-        vector of shape (nsteps) instead of (nsteps, 1).
+        If True and if the system has a single output, return the
+        system output as a 1D array rather than a 2D array.  Default
+        value (True) set by config.defaults['control.squeeze_time_response'].
 
     Returns
     -------
@@ -1420,12 +1420,8 @@ 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)
-        if squeeze:
-            y = np.squeeze(y)
-        if return_x:
-            return T, y, []
-        else:
-            return T, y
+        return _process_time_response(sys, T, y, [], 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)],
@@ -1500,14 +1496,8 @@ def ivp_rhs(t, x): return sys._rhs(t, x, u(t))
     else:                       # Neither ctime or dtime??
         raise TypeError("Can't determine system type")
 
-    # Get rid of extra dimensions in the output, of desired
-    if squeeze:
-        y = np.squeeze(y)
-
-    if return_x:
-        return soln.t, y, soln.y
-    else:
-        return soln.t, y
+    return _process_time_response(sys, soln.t, y, soln.y, transpose=transpose,
+                                  return_x=return_x, squeeze=squeeze)
 
 
 def find_eqpt(sys, x0, u0=[], y0=None, t=0, params={},
diff --git a/control/matlab/timeresp.py b/control/matlab/timeresp.py
index 1ba7b2a0a..31b761bcd 100644
--- a/control/matlab/timeresp.py
+++ b/control/matlab/timeresp.py
@@ -59,15 +59,13 @@ def step(sys, T=None, X0=0., input=0, output=None, return_x=False):
     '''
     from ..timeresp import step_response
 
-    T, yout, xout = step_response(sys, T, X0, input, output,
-                                  transpose=True, return_x=True)
+    # Switch output argument order and transpose outputs
+    out = step_response(sys, T, X0, input, output,
+                        transpose=True, return_x=return_x)
+    return (out[1], out[0], out[2]) if return_x else (out[1], out[0])
 
-    if return_x:
-        return yout, T, xout
-
-    return yout, T
-
-def stepinfo(sys, T=None, SettlingTimeThreshold=0.02, RiseTimeLimits=(0.1, 0.9)):
+def stepinfo(sys, T=None, SettlingTimeThreshold=0.02,
+             RiseTimeLimits=(0.1, 0.9)):
     '''
     Step response characteristics (Rise time, Settling Time, Peak and others).
 
@@ -110,6 +108,7 @@ def stepinfo(sys, T=None, SettlingTimeThreshold=0.02, RiseTimeLimits=(0.1, 0.9))
     '''
     from ..timeresp import step_info
 
+    # Call step_info with MATLAB defaults
     S = step_info(sys, T, None, SettlingTimeThreshold, RiseTimeLimits)
 
     return S
@@ -164,13 +163,11 @@ def impulse(sys, T=None, X0=0., input=0, output=None, return_x=False):
     >>> yout, T = impulse(sys, T)
     '''
     from ..timeresp import impulse_response
-    T, yout, xout = impulse_response(sys, T, X0, input, output,
-                                     transpose = True, return_x=True)
-
-    if return_x:
-        return yout, T, xout
 
-    return yout, T
+    # Switch output argument order and transpose outputs
+    out = impulse_response(sys, T, X0, input, output,
+                           transpose = True, return_x=return_x)
+    return (out[1], out[0], out[2]) if return_x else (out[1], out[0])
 
 def initial(sys, T=None, X0=0., input=None, output=None, return_x=False):
     '''
@@ -222,13 +219,12 @@ def initial(sys, T=None, X0=0., input=None, output=None, return_x=False):
 
     '''
     from ..timeresp import initial_response
+
+    # Switch output argument order and transpose outputs
     T, yout, xout = initial_response(sys, T, X0, output=output,
                                      transpose=True, return_x=True)
+    return (yout, T, xout) if return_x else (yout, T)
 
-    if return_x:
-        return yout, T, xout
-
-    return yout, T
 
 def lsim(sys, U=0., T=None, X0=0.):
     '''
@@ -273,5 +269,7 @@ def lsim(sys, U=0., T=None, X0=0.):
     >>> yout, T, xout = lsim(sys, U, T, X0)
     '''
     from ..timeresp import forced_response
-    T, yout, xout = forced_response(sys, T, U, X0, transpose = True)
-    return yout, T, xout
+
+    # Switch output argument order and transpose outputs (and always return x)
+    out = forced_response(sys, T, U, X0, return_x=True, transpose=True)
+    return out[1], out[0], out[2]
diff --git a/control/tests/config_test.py b/control/tests/config_test.py
index 0e68ec8a7..02d0ad51c 100644
--- a/control/tests/config_test.py
+++ b/control/tests/config_test.py
@@ -211,12 +211,15 @@ def test_legacy_defaults(self):
             ct.use_legacy_defaults('0.8.3')
             assert(isinstance(ct.ss(0, 0, 0, 1).D, np.matrix))
         ct.reset_defaults()
-        assert(isinstance(ct.ss(0, 0, 0, 1).D, np.ndarray))
-        assert(not isinstance(ct.ss(0, 0, 0, 1).D, np.matrix))
+        assert isinstance(ct.ss(0, 0, 0, 1).D, np.ndarray)
+        assert not isinstance(ct.ss(0, 0, 0, 1).D, np.matrix)
+
+        ct.use_legacy_defaults('0.8.4')
+        assert ct.config.defaults['forced_response.return_x'] is True
 
         ct.use_legacy_defaults('0.9.0')
-        assert(isinstance(ct.ss(0, 0, 0, 1).D, np.ndarray))
-        assert(not isinstance(ct.ss(0, 0, 0, 1).D, np.matrix))
+        assert isinstance(ct.ss(0, 0, 0, 1).D, np.ndarray)
+        assert not isinstance(ct.ss(0, 0, 0, 1).D, np.matrix)
 
         # test that old versions don't raise a problem
         ct.use_legacy_defaults('REL-0.1')
diff --git a/control/tests/discrete_test.py b/control/tests/discrete_test.py
index 3dcbb7f3b..379098ff2 100644
--- a/control/tests/discrete_test.py
+++ b/control/tests/discrete_test.py
@@ -353,9 +353,11 @@ def testSimulation(self, tsys):
         tout, yout = step_response(tsys.siso_ss1d, T)
         tout, yout = impulse_response(tsys.siso_ss1d)
         tout, yout = impulse_response(tsys.siso_ss1d, T)
-        tout, yout, xout = forced_response(tsys.siso_ss1d, T, U, 0)
-        tout, yout, xout = forced_response(tsys.siso_ss2d, T, U, 0)
-        tout, yout, xout = forced_response(tsys.siso_ss3d, T, U, 0)
+        tout, yout = forced_response(tsys.siso_ss1d, T, U, 0)
+        tout, yout = forced_response(tsys.siso_ss2d, T, U, 0)
+        tout, yout = forced_response(tsys.siso_ss3d, T, U, 0)
+        tout, yout, xout = forced_response(tsys.siso_ss1d, T, U, 0,
+                                           return_x=True)
 
     def test_sample_system(self, tsys):
         # Make sure we can convert various types of systems
diff --git a/control/tests/flatsys_test.py b/control/tests/flatsys_test.py
index 1281c519a..0239d9455 100644
--- a/control/tests/flatsys_test.py
+++ b/control/tests/flatsys_test.py
@@ -48,7 +48,7 @@ def test_double_integrator(self, xf, uf, Tf):
         T = np.linspace(0, Tf, 100)
         xd, ud = traj.eval(T)
 
-        t, y, x = ct.forced_response(sys, T, ud, x1)
+        t, y, x = ct.forced_response(sys, T, ud, x1, return_x=True)
         np.testing.assert_array_almost_equal(x, xd, decimal=3)
 
     def test_kinematic_car(self):
diff --git a/control/tests/iosys_test.py b/control/tests/iosys_test.py
index 0dcbf3325..32e5b102f 100644
--- a/control/tests/iosys_test.py
+++ b/control/tests/iosys_test.py
@@ -61,7 +61,7 @@ def test_linear_iosys(self, tsys):
 
         # Make sure that simulations also line up
         T, U, X0 = tsys.T, tsys.U, tsys.X0
-        lti_t, lti_y, lti_x = ct.forced_response(linsys, T, U, X0)
+        lti_t, lti_y = ct.forced_response(linsys, T, U, X0)
         ios_t, ios_y = ios.input_output_response(iosys, T, U, X0)
         np.testing.assert_array_almost_equal(lti_t, ios_t)
         np.testing.assert_allclose(lti_y, ios_y, atol=0.002, rtol=0.)
@@ -75,7 +75,7 @@ def test_tf2io(self, tsys):
 
         # Verify correctness via simulation
         T, U, X0 = tsys.T, tsys.U, tsys.X0
-        lti_t, lti_y, lti_x = ct.forced_response(linsys, T, U, X0)
+        lti_t, lti_y = ct.forced_response(linsys, T, U, X0)
         ios_t, ios_y = ios.input_output_response(iosys, T, U, X0)
         np.testing.assert_array_almost_equal(lti_t, ios_t)
         np.testing.assert_allclose(lti_y, ios_y, atol=0.002, rtol=0.)
@@ -84,7 +84,7 @@ def test_tf2io(self, tsys):
         tfsys = ct.tf('s')
         with pytest.raises(ValueError):
             iosys=ct.tf2io(tfsys)
-            
+
     def test_ss2io(self, tsys):
         # Create an input/output system from the linear system
         linsys = tsys.siso_linsys
@@ -162,7 +162,7 @@ def test_nonlinear_iosys(self, tsys):
 
         # Make sure that simulations also line up
         T, U, X0 = tsys.T, tsys.U, tsys.X0
-        lti_t, lti_y, lti_x = ct.forced_response(linsys, T, U, X0)
+        lti_t, lti_y = ct.forced_response(linsys, T, U, X0)
         ios_t, ios_y = ios.input_output_response(nlsys, T, U, X0)
         np.testing.assert_array_almost_equal(lti_t, ios_t)
         np.testing.assert_allclose(lti_y, ios_y,atol=0.002,rtol=0.)
@@ -256,7 +256,7 @@ def test_connect(self, tsys):
         X0 = np.concatenate((tsys.X0, tsys.X0))
         ios_t, ios_y, ios_x = ios.input_output_response(
             iosys_series, T, U, X0, return_x=True)
-        lti_t, lti_y, lti_x = ct.forced_response(linsys_series, T, U, X0)
+        lti_t, lti_y = ct.forced_response(linsys_series, T, U, X0)
         np.testing.assert_array_almost_equal(lti_t, ios_t)
         np.testing.assert_allclose(lti_y, ios_y,atol=0.002,rtol=0.)
 
@@ -273,7 +273,7 @@ def test_connect(self, tsys):
         assert ct.isctime(iosys_series, strict=True)
         ios_t, ios_y, ios_x = ios.input_output_response(
             iosys_series, T, U, X0, return_x=True)
-        lti_t, lti_y, lti_x = ct.forced_response(linsys_series, T, U, X0)
+        lti_t, lti_y = ct.forced_response(linsys_series, T, U, X0)
         np.testing.assert_array_almost_equal(lti_t, ios_t)
         np.testing.assert_allclose(lti_y, ios_y,atol=0.002,rtol=0.)
 
@@ -288,7 +288,7 @@ def test_connect(self, tsys):
         )
         ios_t, ios_y, ios_x = ios.input_output_response(
             iosys_feedback, T, U, X0, return_x=True)
-        lti_t, lti_y, lti_x = ct.forced_response(linsys_feedback, T, U, X0)
+        lti_t, lti_y = ct.forced_response(linsys_feedback, T, U, X0)
         np.testing.assert_array_almost_equal(lti_t, ios_t)
         np.testing.assert_allclose(lti_y, ios_y,atol=0.002,rtol=0.)
 
@@ -325,7 +325,8 @@ def test_connect_spec_variants(self, tsys, connections, inplist, outlist):
         # Create a simulation run to compare against
         T, U = tsys.T, tsys.U
         X0 = np.concatenate((tsys.X0, tsys.X0))
-        lti_t, lti_y, lti_x = ct.forced_response(linsys_series, T, U, X0)
+        lti_t, lti_y, lti_x = ct.forced_response(
+            linsys_series, T, U, X0, return_x=True)
 
         # Create the input/output system with different parameter variations
         iosys_series = ios.InterconnectedSystem(
@@ -360,7 +361,8 @@ def test_connect_spec_warnings(self, tsys, connections, inplist, outlist):
         # Create a simulation run to compare against
         T, U = tsys.T, tsys.U
         X0 = np.concatenate((tsys.X0, tsys.X0))
-        lti_t, lti_y, lti_x = ct.forced_response(linsys_series, T, U, X0)
+        lti_t, lti_y, lti_x = ct.forced_response(
+            linsys_series, T, U, X0, return_x=True)
 
         # Set up multiple gainst and make sure a warning is generated
         with pytest.warns(UserWarning, match="multiple.*Combining"):
@@ -388,7 +390,8 @@ def test_static_nonlinearity(self, tsys):
 
         # Make sure saturation works properly by comparing linear system with
         # saturated input to nonlinear system with saturation composition
-        lti_t, lti_y, lti_x = ct.forced_response(linsys, T, Usat, X0)
+        lti_t, lti_y, lti_x = ct.forced_response(
+            linsys, T, Usat, X0, return_x=True)
         ios_t, ios_y, ios_x = ios.input_output_response(
             ioslin * nlsat, T, U, X0, return_x=True)
         np.testing.assert_array_almost_equal(lti_t, ios_t)
@@ -424,7 +427,7 @@ def test_algebraic_loop(self, tsys):
         # Nonlinear system composed with LTI system (series) -- with states
         ios_t, ios_y = ios.input_output_response(
             nlios * lnios * nlios, T, U, X0)
-        lti_t, lti_y, lti_x = ct.forced_response(linsys, T, U*U, X0)
+        lti_t, lti_y = ct.forced_response(linsys, T, U*U, X0)
         np.testing.assert_array_almost_equal(ios_y, lti_y*lti_y, decimal=3)
 
         # Nonlinear system in feeback loop with LTI system
@@ -480,7 +483,7 @@ def test_summer(self, tsys):
         U = [np.sin(T), np.cos(T)]
         X0 = 0
 
-        lin_t, lin_y, lin_x = ct.forced_response(linsys_parallel, T, U, X0)
+        lin_t, lin_y = ct.forced_response(linsys_parallel, T, U, X0)
         ios_t, ios_y = ios.input_output_response(iosys_parallel, T, U, X0)
         np.testing.assert_allclose(ios_y, lin_y,atol=0.002,rtol=0.)
 
@@ -502,7 +505,7 @@ def test_rmul(self, tsys):
 
         # Make sure we got the right thing (via simulation comparison)
         ios_t, ios_y = ios.input_output_response(sys2, T, U, X0)
-        lti_t, lti_y, lti_x = ct.forced_response(ioslin, T, U*U, X0)
+        lti_t, lti_y = ct.forced_response(ioslin, T, U*U, X0)
         np.testing.assert_array_almost_equal(ios_y, lti_y*lti_y, decimal=3)
 
     @noscipy0
@@ -525,7 +528,7 @@ def test_neg(self, tsys):
 
         # Make sure we got the right thing (via simulation comparison)
         ios_t, ios_y = ios.input_output_response(sys, T, U, X0)
-        lti_t, lti_y, lti_x = ct.forced_response(ioslin, T, U*U, X0)
+        lti_t, lti_y = ct.forced_response(ioslin, T, U*U, X0)
         np.testing.assert_array_almost_equal(ios_y, -lti_y, decimal=3)
 
     @noscipy0
@@ -541,7 +544,7 @@ def test_feedback(self, tsys):
         linsys = ct.feedback(tsys.siso_linsys, 1)
 
         ios_t, ios_y = ios.input_output_response(iosys, T, U, X0)
-        lti_t, lti_y, lti_x = ct.forced_response(linsys, T, U, X0)
+        lti_t, lti_y = ct.forced_response(linsys, T, U, X0)
         np.testing.assert_allclose(ios_y, lti_y,atol=0.002,rtol=0.)
 
     @noscipy0
@@ -561,33 +564,33 @@ def test_bdalg_functions(self, tsys):
         # Series interconnection
         linsys_series = ct.series(linsys1, linsys2)
         iosys_series = ct.series(linio1, linio2)
-        lin_t, lin_y, lin_x = ct.forced_response(linsys_series, T, U, X0)
+        lin_t, lin_y = ct.forced_response(linsys_series, T, U, X0)
         ios_t, ios_y = ios.input_output_response(iosys_series, T, U, X0)
         np.testing.assert_allclose(ios_y, lin_y,atol=0.002,rtol=0.)
 
         # Make sure that systems don't commute
         linsys_series = ct.series(linsys2, linsys1)
-        lin_t, lin_y, lin_x = ct.forced_response(linsys_series, T, U, X0)
+        lin_t, lin_y = ct.forced_response(linsys_series, T, U, X0)
         assert not (np.abs(lin_y - ios_y) < 1e-3).all()
 
         # Parallel interconnection
         linsys_parallel = ct.parallel(linsys1, linsys2)
         iosys_parallel = ct.parallel(linio1, linio2)
-        lin_t, lin_y, lin_x = ct.forced_response(linsys_parallel, T, U, X0)
+        lin_t, lin_y = ct.forced_response(linsys_parallel, T, U, X0)
         ios_t, ios_y = ios.input_output_response(iosys_parallel, T, U, X0)
         np.testing.assert_allclose(ios_y, lin_y,atol=0.002,rtol=0.)
 
         # Negation
         linsys_negate = ct.negate(linsys1)
         iosys_negate = ct.negate(linio1)
-        lin_t, lin_y, lin_x = ct.forced_response(linsys_negate, T, U, X0)
+        lin_t, lin_y = ct.forced_response(linsys_negate, T, U, X0)
         ios_t, ios_y = ios.input_output_response(iosys_negate, T, U, X0)
         np.testing.assert_allclose(ios_y, lin_y,atol=0.002,rtol=0.)
 
         # Feedback interconnection
         linsys_feedback = ct.feedback(linsys1, linsys2)
         iosys_feedback = ct.feedback(linio1, linio2)
-        lin_t, lin_y, lin_x = ct.forced_response(linsys_feedback, T, U, X0)
+        lin_t, lin_y = ct.forced_response(linsys_feedback, T, U, X0)
         ios_t, ios_y = ios.input_output_response(iosys_feedback, T, U, X0)
         np.testing.assert_allclose(ios_y, lin_y,atol=0.002,rtol=0.)
 
@@ -614,13 +617,13 @@ def test_nonsquare_bdalg(self, tsys):
         # Multiplication
         linsys_multiply = linsys_3i2o * linsys_2i3o
         iosys_multiply = iosys_3i2o * iosys_2i3o
-        lin_t, lin_y, lin_x = ct.forced_response(linsys_multiply, T, U2, X0)
+        lin_t, lin_y = ct.forced_response(linsys_multiply, T, U2, X0)
         ios_t, ios_y = ios.input_output_response(iosys_multiply, T, U2, X0)
         np.testing.assert_allclose(ios_y, lin_y,atol=0.002,rtol=0.)
 
         linsys_multiply = linsys_2i3o * linsys_3i2o
         iosys_multiply = iosys_2i3o * iosys_3i2o
-        lin_t, lin_y, lin_x = ct.forced_response(linsys_multiply, T, U3, X0)
+        lin_t, lin_y = ct.forced_response(linsys_multiply, T, U3, X0)
         ios_t, ios_y = ios.input_output_response(iosys_multiply, T, U3, X0)
         np.testing.assert_allclose(ios_y, lin_y,atol=0.002,rtol=0.)
 
@@ -633,7 +636,7 @@ def test_nonsquare_bdalg(self, tsys):
         # Feedback
         linsys_multiply = ct.feedback(linsys_3i2o, linsys_2i3o)
         iosys_multiply = iosys_3i2o.feedback(iosys_2i3o)
-        lin_t, lin_y, lin_x = ct.forced_response(linsys_multiply, T, U3, X0)
+        lin_t, lin_y = ct.forced_response(linsys_multiply, T, U3, X0)
         ios_t, ios_y = ios.input_output_response(iosys_multiply, T, U3, X0)
         np.testing.assert_allclose(ios_y, lin_y,atol=0.002,rtol=0.)
 
@@ -655,7 +658,7 @@ def test_discrete(self, tsys):
 
         # Simulate and compare to LTI output
         ios_t, ios_y = ios.input_output_response(lnios, T, U, X0)
-        lin_t, lin_y, lin_x = ct.forced_response(linsys, T, U, X0)
+        lin_t, lin_y = ct.forced_response(linsys, T, U, X0)
         np.testing.assert_allclose(ios_t, lin_t,atol=0.002,rtol=0.)
         np.testing.assert_allclose(ios_y, lin_y,atol=0.002,rtol=0.)
 
@@ -671,7 +674,7 @@ def test_discrete(self, tsys):
 
         # Simulate and compare to LTI output
         ios_t, ios_y = ios.input_output_response(lnios, T, U, X0)
-        lin_t, lin_y, lin_x = ct.forced_response(linsys, T, U, X0)
+        lin_t, lin_y = ct.forced_response(linsys, T, U, X0)
         np.testing.assert_allclose(ios_t, lin_t,atol=0.002,rtol=0.)
         np.testing.assert_allclose(ios_y, lin_y,atol=0.002,rtol=0.)
 
@@ -839,7 +842,7 @@ def test_params(self, tsys):
         linsys = tsys.siso_linsys
         iosys = ios.LinearIOSystem(linsys)
         T, U, X0 = tsys.T, tsys.U, tsys.X0
-        lti_t, lti_y, lti_x = ct.forced_response(linsys, T, U, X0)
+        lti_t, lti_y = ct.forced_response(linsys, T, U, X0)
         with pytest.warns(UserWarning, match="LinearIOSystem.*ignored"):
             ios_t, ios_y = ios.input_output_response(
                 iosys, T, U, X0, params={'something':0})
diff --git a/control/tests/modelsimp_test.py b/control/tests/modelsimp_test.py
index 4def0b4d7..df656e1fc 100644
--- a/control/tests/modelsimp_test.py
+++ b/control/tests/modelsimp_test.py
@@ -51,7 +51,7 @@ def testMarkovSignature(self, matarrayout, matarrayin):
         # Test example from docstring
         T = np.linspace(0, 10, 100)
         U = np.ones((1, 100))
-        T, Y, _ = forced_response(tf([1], [1, 0.5], True), T, U)
+        T, Y = forced_response(tf([1], [1, 0.5], True), T, U)
         H = markov(Y, U, 3, transpose=False)
 
         # Test example from issue #395
@@ -102,7 +102,7 @@ def testMarkovResults(self, k, m, n):
         # Generate input/output data
         T = np.array(range(n)) * Ts
         U = np.cos(T) + np.sin(T/np.pi)
-        _, Y, _ = forced_response(Hd, T, U, squeeze=True)
+        _, Y = forced_response(Hd, T, U, squeeze=True)
         Mcomp = markov(Y, U, m)
 
         # Compare to results from markov()
diff --git a/control/tests/timeresp_test.py b/control/tests/timeresp_test.py
index 6977973ff..41bf05c4b 100644
--- a/control/tests/timeresp_test.py
+++ b/control/tests/timeresp_test.py
@@ -16,12 +16,14 @@
 import scipy as sp
 
 
+import control as ct
 from control import (StateSpace, TransferFunction, c2d, isctime, isdtime,
                      ss2tf, tf2ss)
 from control.timeresp import (_default_time_vector, _ideal_tfinal_and_dt,
                               forced_response, impulse_response,
                               initial_response, step_info, step_response)
 from control.tests.conftest import slycotonly
+from control.exception import slycot_check
 
 
 class TSys:
@@ -409,7 +411,7 @@ def test_forced_response_step(self, tsystem):
         u = np.ones_like(t, dtype=np.float)
         yref = tsystem.ystep
 
-        tout, yout, _xout = forced_response(sys, t, u)
+        tout, yout = forced_response(sys, t, u)
         np.testing.assert_array_almost_equal(tout, t)
         np.testing.assert_array_almost_equal(yout, yref, decimal=4)
 
@@ -424,7 +426,7 @@ def test_forced_response_initial(self, siso_ss1, u):
         x0 = np.array([[.5], [1.]])
         yref = siso_ss1.yinitial
 
-        tout, yout, _xout = forced_response(sys, t, u, X0=x0)
+        tout, yout = forced_response(sys, t, u, X0=x0)
         np.testing.assert_array_almost_equal(tout, t)
         np.testing.assert_array_almost_equal(yout, yref, decimal=4)
 
@@ -444,11 +446,31 @@ def test_forced_response_mimo(self, tsystem, useT):
         yref = np.vstack([tsystem.yinitial, tsystem.ystep])
 
         if useT:
-            _t, yout, _xout = forced_response(sys, t, u, x0)
+            _t, yout = forced_response(sys, t, u, x0)
         else:
-            _t, yout, _xout = forced_response(sys, U=u, X0=x0)
+            _t, yout = forced_response(sys, U=u, X0=x0)
         np.testing.assert_array_almost_equal(yout, yref, decimal=4)
 
+    @pytest.mark.usefixtures("editsdefaults")
+    def test_forced_response_legacy(self):
+        # Define a system for testing
+        sys = ct.rss(2, 1, 1)
+        T = np.linspace(0, 10, 10)
+        U = np.sin(T)
+
+        """Make sure that legacy version of forced_response works"""
+        ct.config.use_legacy_defaults("0.8.4")
+        # forced_response returns x by default
+        t, y = ct.step_response(sys, T)
+        t, y, x = ct.forced_response(sys, T, U)
+
+        ct.config.use_legacy_defaults("0.9.0")
+        # forced_response returns input/output by default
+        t, y = ct.step_response(sys, T)
+        t, y = ct.forced_response(sys, T, U)
+        t, y, x = ct.forced_response(sys, T, U, return_x=True)
+
+
     @pytest.mark.parametrize("u, x0, xtrue",
                              [(np.zeros((10,)),
                                np.array([2., 3.]),
@@ -475,7 +497,7 @@ def test_lsim_double_integrator(self, u, x0, xtrue):
         sys = StateSpace(A, B, C, D)
         t = np.linspace(0, 1, 10)
 
-        _t, yout, xout = forced_response(sys, t, u, x0)
+        _t, yout, xout = forced_response(sys, t, u, x0, return_x=True)
         np.testing.assert_array_almost_equal(xout, xtrue, decimal=6)
         ytrue = np.squeeze(np.asarray(C.dot(xtrue)))
         np.testing.assert_array_almost_equal(yout, ytrue, decimal=6)
@@ -643,8 +665,8 @@ def test_time_vector_interpolation(self, siso_dtf2, squeeze):
 
         squeezekw = {} if squeeze is None else {"squeeze": squeeze}
 
-        tout, yout, xout = forced_response(sys, t, u, x0,
-                                           interpolate=True, **squeezekw)
+        tout, yout = forced_response(sys, t, u, x0,
+                                     interpolate=True, **squeezekw)
         if squeeze is False or sys.outputs > 1:
             assert yout.shape[0] == sys.outputs
             assert yout.shape[1] == tout.shape[0]
@@ -700,3 +722,117 @@ def test_time_series_data_convention_2D(self, siso_ss1):
         assert t.ndim == 1
         assert y.ndim == 1  # SISO returns "scalar" output
         assert t.shape == y.shape  # Allows direct plotting of output
+
+    @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)],
+#       [4, 1, 1, 'siso', (8,)],        # Use for later 'siso' implementation
+        [1, 2, 1, None, (2, 8)],
+        [2, 2, 1, True, (2, 8)],
+        [3, 2, 1, False, (2, 8)],
+#       [4, 2, 1, 'siso', (2, 8)],      # Use for later 'siso' implementation
+        [1, 1, 2, None, (8,)],
+        [2, 1, 2, True, (8,)],
+        [3, 1, 2, False, (1, 8)],
+#       [4, 1, 2, 'siso', (1, 8)],      # Use for later 'siso' implementation
+        [1, 2, 2, None, (2, 8)],
+        [2, 2, 2, True, (2, 8)],
+        [3, 2, 2, False, (2, 8)],
+#       [4, 2, 2, 'siso', (2, 8)],      # Use for later 'siso' implementation
+    ])
+    def test_squeeze(self, fcn, nstate, nout, ninp, squeeze, shape):
+        # Figure out if we have SciPy 1+
+        scipy0 = StrictVersion(sp.__version__) < '1.0'
+
+        # Define the system
+        if fcn == ct.tf and (nout > 1 or ninp > 1) and not slycot_check():
+            pytest.skip("Conversion of MIMO systems to transfer functions "
+                        "requires slycot.")
+        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
+
+        _, yvec = ct.initial_response(sys, tvec, 1, squeeze=squeeze)
+        assert yvec.shape == shape
+
+        with pytest.warns(warntype, match="Converting MIMO system"):
+            _, yvec = ct.step_response(sys, tvec, squeeze=squeeze)
+        assert yvec.shape == shape
+
+        if isinstance(sys, StateSpace):
+            # Check the states as well
+            _, yvec, xvec = ct.forced_response(
+                sys, tvec, uvec, 0, return_x=True, squeeze=squeeze)
+            assert xvec.shape == (sys.states, 8)
+        else:
+            # Just check the input/output response
+            _, yvec = ct.forced_response(sys, tvec, uvec, 0, squeeze=squeeze)
+        assert yvec.shape == shape
+
+        # 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 or (squeeze == 'siso' and not sys.issiso()):
+        if squeeze is False:
+            # Shape should be unsqueezed
+            assert yvec.shape == (1, 8)
+        else:
+            # Shape should be squeezed
+            assert yvec.shape == (8, )
+
+        # 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
+
+        #
+        # 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.initial_response(sys, tvec, 1)
+        assert yvec.shape == (sys.outputs, 8)
+
+        with pytest.warns(warntype, match="Converting MIMO system"):
+            _, yvec = ct.step_response(sys, tvec)
+        assert yvec.shape == (sys.outputs, 8)
+
+        _, yvec, xvec = ct.forced_response(
+            sys, tvec, uvec, 0, return_x=True)
+        assert yvec.shape == (sys.outputs, 8)
+        if isinstance(sys, ct.StateSpace):
+            assert xvec.shape == (sys.states, 8)
+        else:
+            assert xvec.shape[1] == 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)
+
+    @pytest.mark.parametrize("fcn", [ct.ss, ct.tf, ct.ss2io])
+    def test_squeeze_exception(self, fcn):
+        sys = fcn(ct.rss(2, 1, 1))
+        with pytest.raises(ValueError, match="unknown squeeze value"):
+            step_response(sys, squeeze=1)
diff --git a/control/timeresp.py b/control/timeresp.py
index a5cc245bf..7a48a5c5f 100644
--- a/control/timeresp.py
+++ b/control/timeresp.py
@@ -79,8 +79,10 @@
                    atleast_1d)
 import warnings
 from .lti import LTI     # base class of StateSpace, TransferFunction
+from .xferfcn import TransferFunction
 from .statesp import _convert_to_statespace, _mimo2simo, _mimo2siso, ssdata
 from .lti import isdtime, isctime
+from . import config
 
 __all__ = ['forced_response', 'step_response', 'step_info', 'initial_response',
            'impulse_response']
@@ -102,10 +104,10 @@ def _check_convert_array(in_obj, legal_shapes, err_msg_start, squeeze=False,
 
     Parameters
     ----------
-    in_obj: array like object
+    in_obj : array like object
         The array or matrix which is checked.
 
-    legal_shapes: list of tuple
+    legal_shapes : list of tuple
         A list of shapes that in_obj can legally have.
         The special value "any" means that there can be any
         number of elements in a certain dimension.
@@ -114,25 +116,26 @@ def _check_convert_array(in_obj, legal_shapes, err_msg_start, squeeze=False,
         * ``(2, "any")`` describes an array with 2 rows and any number of
           columns
 
-    err_msg_start: str
+    err_msg_start : str
         String that is prepended to the error messages, when this function
         raises an exception. It should be used to identify the argument which
         is currently checked.
 
-    squeeze: bool
+    squeeze : bool
         If True, all dimensions with only one element are removed from the
         array. If False the array's shape is unmodified.
 
         For example:
         ``array([[1,2,3]])`` is converted to ``array([1, 2, 3])``
 
-   transpose: bool
+   transpose : bool
         If True, assume that input arrays are transposed for the standard
         format.  Used to convert MATLAB-style inputs to our format.
 
-    Returns:
+    Returns
+    -------
 
-    out_array: array
+    out_array : array
         The checked and converted contents of ``in_obj``.
     """
     # convert nearly everything to an array.
@@ -195,7 +198,7 @@ def shape_matches(s_legal, s_actual):
 
 # Forced response of a linear system
 def forced_response(sys, T=None, U=0., X0=0., transpose=False,
-                    interpolate=False, squeeze=True):
+                    interpolate=False, return_x=None, squeeze=None):
     """Simulate the output of a linear system.
 
     As a convenience for parameters `U`, `X0`:
@@ -207,45 +210,50 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False,
 
     Parameters
     ----------
-    sys: LTI (StateSpace or TransferFunction)
+    sys : LTI (StateSpace or TransferFunction)
         LTI system to simulate
 
-    T: array_like, optional for discrete LTI `sys`
+    T : array_like, optional for discrete LTI `sys`
         Time steps at which the input is defined; values must be evenly spaced.
 
-    U: array_like or float, optional
+    U : array_like or float, optional
         Input array giving input at each time `T` (default = 0).
 
         If `U` is ``None`` or ``0``, a special algorithm is used. This special
         algorithm is faster than the general algorithm, which is used
         otherwise.
 
-    X0: array_like or float, optional
+    X0 : array_like or float, optional
         Initial condition (default = 0).
 
-    transpose: bool, optional (default=False)
+    transpose : bool, optional (default=False)
         If True, transpose all input and output arrays (for backward
         compatibility with MATLAB and :func:`scipy.signal.lsim`)
 
-    interpolate: bool, optional (default=False)
+    interpolate : bool, optional (default=False)
         If True and system is a discrete time system, the input will
         be interpolated between the given time steps and the output
         will be given at system sampling rate.  Otherwise, only return
         the output at the times given in `T`.  No effect on continuous
         time simulations (default = False).
 
-    squeeze: bool, optional (default=True)
-        If True, remove single-dimensional entries from the shape of
-        the output.  For single output systems, this converts the
-        output response to a 1D array.
+    return_x : bool, optional
+        If True (default), return the the state vector.  Set to False to
+        return only the time and output vectors.
+
+    squeeze : bool, optional
+        If True (default), remove single-dimensional entries from the shape
+        of the output.  For single output systems, this converts the output
+        response to a 1D array.  The default value can be set using
+        config.defaults['control.squeeze_time_response'].
 
     Returns
     -------
-    T: array
+    T : array
         Time values of the output.
-    yout: array
+    yout : array
         Response of the system.
-    xout: array
+    xout : array
         Time evolution of the state vector.
 
     See Also
@@ -271,6 +279,17 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False,
     if not isinstance(sys, LTI):
         raise TypeError('Parameter ``sys``: must be a ``LTI`` object. '
                         '(For example ``StateSpace`` or ``TransferFunction``)')
+
+    # If return_x was not specified, figure out the default
+    if return_x is None:
+        return_x = config.defaults['forced_response.return_x']
+
+    # If return_x is used for TransferFunction, issue a warning
+    if return_x and not isinstance(sys, TransferFunction):
+        warnings.warn(
+            "return_x specified for a transfer function system. Internal "
+            "conversation to state space used; results may meaningless.")
+
     sys = _convert_to_statespace(sys)
     A, B, C, D = np.asarray(sys.A), np.asarray(sys.B), np.asarray(sys.C), \
         np.asarray(sys.D)
@@ -324,6 +343,13 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False,
     X0 = _check_convert_array(X0, [(n_states,), (n_states, 1)],
                               'Parameter ``X0``: ', squeeze=True)
 
+    # If we are passed a transfer function and X0 is non-zero, warn the user
+    if isinstance(sys, TransferFunction) and np.any(X0 != 0):
+        warnings.warn(
+            "Non-zero initial condition given for transfer function system. "
+            "Internal conversation to state space used; may not be consistent "
+            "with given X0.")
+
     xout = np.zeros((n_states, n_steps))
     xout[:, 0] = X0
     yout = np.zeros((n_outputs, n_steps))
@@ -417,10 +443,28 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False,
         xout = np.transpose(xout)
         yout = np.transpose(yout)
 
-    # Get rid of unneeded dimensions
-    if squeeze:
+    return _process_time_response(sys, tout, yout, xout, transpose=transpose,
+                                  return_x=return_x, squeeze=squeeze)
+
+
+# Process time responses in a uniform way
+def _process_time_response(sys, tout, yout, xout, transpose=None,
+                           return_x=False, squeeze=None):
+    # If squeeze was not specified, figure out the default
+    if squeeze is None:
+        squeeze = config.defaults['control.squeeze_time_response']
+
+    # Figure out whether and now to squeeze output data
+    if squeeze is True:         # squeeze all dimensions
         yout = np.squeeze(yout)
-        xout = np.squeeze(xout)
+    elif squeeze is False:      # squeeze no dimensions
+        pass
+    # Possible code if we implement a squeeze='siso' option
+    # elif squeeze == 'siso':     # squeeze signals if SISO
+    #     yout = yout[0] if sys.issiso() else yout
+    else:
+        # In preparation for more complicated values for squeeze
+        raise ValueError("unknown squeeze value")
 
     # See if we need to transpose the data back into MATLAB form
     if transpose:
@@ -428,30 +472,37 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False,
         yout = np.transpose(yout)
         xout = np.transpose(xout)
 
-    return tout, yout, xout
+    # Return time, output, and (optionally) state
+    return (tout, yout, xout) if return_x else (tout, yout)
 
 
-def _get_ss_simo(sys, input=None, output=None):
+def _get_ss_simo(sys, input=None, output=None, squeeze=None):
     """Return a SISO or SIMO state-space version of sys
 
     If input is not specified, select first input and issue warning
     """
     sys_ss = _convert_to_statespace(sys)
     if sys_ss.issiso():
-        return sys_ss
+        return squeeze, sys_ss
+    # Possible code if we implement a squeeze='siso' option
+    # elif squeeze == 'siso':
+    #     # Don't squeeze outputs if resulting system turns out to be siso
+    #     squeeze = False
+
     warn = False
     if input is None:
         # issue warning if input is not given
         warn = True
         input = 0
+
     if output is None:
-        return _mimo2simo(sys_ss, input, warn_conversion=warn)
+        return squeeze, _mimo2simo(sys_ss, input, warn_conversion=warn)
     else:
-        return _mimo2siso(sys_ss, input, output, warn_conversion=warn)
+        return squeeze, _mimo2siso(sys_ss, input, output, warn_conversion=warn)
 
 
 def step_response(sys, T=None, X0=0., input=None, output=None, T_num=None,
-                  transpose=False, return_x=False, squeeze=True):
+                  transpose=False, return_x=False, squeeze=None):
     # pylint: disable=W0622
     """Step response of a linear system
 
@@ -465,10 +516,10 @@ def step_response(sys, T=None, X0=0., input=None, output=None, T_num=None,
 
     Parameters
     ----------
-    sys: StateSpace or TransferFunction
+    sys : StateSpace or TransferFunction
         LTI system to simulate
 
-    T: array_like or float, optional
+    T : array_like or float, optional
         Time vector, or simulation time duration if a number. If T is not
         provided, an attempt is made to create it automatically from the
         dynamics of sys. If sys is continuous-time, the time increment dt
@@ -479,44 +530,45 @@ def step_response(sys, T=None, X0=0., input=None, output=None, T_num=None,
         only tfinal is computed, and final is reduced if it requires too
         many simulation steps.
 
-    X0: array_like or float, optional
+    X0 : array_like or float, optional
         Initial condition (default = 0)
 
         Numbers are converted to constant arrays with the correct shape.
 
-    input: int
-        Index of the input that will be used in this simulation.
+    input : int
+        Index of the input that will be used in this simulation.  Default = 0.
 
-    output: int
+    output : int
         Index of the output that will be used in this simulation. Set to None
         to not trim outputs
 
-    T_num: int, optional
+    T_num : int, optional
         Number of time steps to use in simulation if T is not provided as an
         array (autocomputed if not given); ignored if sys is discrete-time.
 
-    transpose: bool
+    transpose : bool
         If True, transpose all input and output arrays (for backward
         compatibility with MATLAB and :func:`scipy.signal.lsim`)
 
-    return_x: bool
+    return_x : bool
         If True, return the state vector (default = False).
 
-    squeeze: bool, optional (default=True)
-        If True, remove single-dimensional entries from the shape of
-        the output.  For single output systems, this converts the
-        output response to a 1D array.
+    squeeze : bool, optional
+        If True (default), remove single-dimensional entries from the shape
+        of the output.  For single output systems, this converts the output
+        response to a 1D array.  The default value can be set using
+        config.defaults['control.squeeze_time_response'].
 
     Returns
     -------
-    T: array
+    T : array
         Time values of the output
 
-    yout: array
+    yout : array
         Response of the system
 
-    xout: array
-        Individual response of each x variable
+    xout : array, optional
+        Individual response of each x variable (if return_x is True)
 
     See Also
     --------
@@ -532,18 +584,13 @@ def step_response(sys, T=None, X0=0., input=None, output=None, T_num=None,
     >>> T, yout = step_response(sys, T, X0)
 
     """
-    sys = _get_ss_simo(sys, input, output)
+    squeeze, sys = _get_ss_simo(sys, input, output, squeeze=squeeze)
     if T is None or np.asarray(T).size == 1:
         T = _default_time_vector(sys, N=T_num, tfinal=T, is_step=True)
     U = np.ones_like(T)
 
-    T, yout, xout = forced_response(sys, T, U, X0, transpose=transpose,
-                                    squeeze=squeeze)
-
-    if return_x:
-        return T, yout, xout
-
-    return T, yout
+    return forced_response(sys, T, U, X0, transpose=transpose,
+                           return_x=return_x, squeeze=squeeze)
 
 
 def step_info(sys, T=None, T_num=None, SettlingTimeThreshold=0.02,
@@ -592,7 +639,7 @@ def step_info(sys, T=None, T_num=None, SettlingTimeThreshold=0.02,
     --------
     >>> info = step_info(sys, T)
     '''
-    sys = _get_ss_simo(sys)
+    _, sys = _get_ss_simo(sys)
     if T is None or np.asarray(T).size == 1:
         T = _default_time_vector(sys, N=T_num, tfinal=T, is_step=True)
 
@@ -630,7 +677,7 @@ def step_info(sys, T=None, T_num=None, SettlingTimeThreshold=0.02,
 
 
 def initial_response(sys, T=None, X0=0., input=0, output=None, T_num=None,
-                     transpose=False, return_x=False, squeeze=True):
+                     transpose=False, return_x=False, squeeze=None):
     # pylint: disable=W0622
     """Initial condition response of a linear system
 
@@ -651,33 +698,34 @@ def initial_response(sys, T=None, X0=0., input=0, output=None, T_num=None,
         autocomputed if not given; see  :func:`step_response` for more detail)
 
     X0 : array_like or float, optional
-        Initial condition (default = 0)
-
-        Numbers are converted to constant arrays with the correct shape.
+        Initial condition (default = 0).  Numbers are converted to constant
+        arrays with the correct shape.
 
     input : int
         Ignored, has no meaning in initial condition calculation. Parameter
-        ensures compatibility with step_response and impulse_response
+        ensures compatibility with step_response and impulse_response.
 
     output : int
         Index of the output that will be used in this simulation. Set to None
-        to not trim outputs
+        to not trim outputs.
 
     T_num : int, optional
         Number of time steps to use in simulation if T is not provided as an
         array (autocomputed if not given); ignored if sys is discrete-time.
 
-    transpose : bool
+    transpose : bool, optional
         If True, transpose all input and output arrays (for backward
-        compatibility with MATLAB and :func:`scipy.signal.lsim`)
+        compatibility with MATLAB and :func:`scipy.signal.lsim`).  Default
+        value is False.
 
-    return_x : bool
+    return_x : bool, optional
         If True, return the state vector (default = False).
 
-    squeeze : bool, optional (default=True)
-        If True, remove single-dimensional entries from the shape of
-        the output.  For single output systems, this converts the
-        output response to a 1D array.
+    squeeze : bool, optional
+        If True (default), remove single-dimensional entries from the shape
+        of the output.  For single output systems, this converts the output
+        response to a 1D array.  The default value can be set using
+        config.defaults['control.squeeze_time_response'].
 
     Returns
     -------
@@ -685,8 +733,8 @@ def initial_response(sys, T=None, X0=0., input=0, output=None, T_num=None,
         Time values of the output
     yout : array
         Response of the system
-    xout : array
-        Individual response of each x variable
+    xout : array, optional
+        Individual response of each x variable (if return_x is True)
 
     See Also
     --------
@@ -700,8 +748,9 @@ def initial_response(sys, T=None, X0=0., input=0, output=None, T_num=None,
     Examples
     --------
     >>> T, yout = initial_response(sys, T, X0)
+
     """
-    sys = _get_ss_simo(sys, input, output)
+    squeeze, sys = _get_ss_simo(sys, input, output, squeeze=squeeze)
 
     # Create time and input vectors; checking is done in forced_response(...)
     # The initial vector X0 is created in forced_response(...) if necessary
@@ -709,17 +758,12 @@ def initial_response(sys, T=None, X0=0., input=0, output=None, T_num=None,
         T = _default_time_vector(sys, N=T_num, tfinal=T, is_step=False)
     U = np.zeros_like(T)
 
-    T, yout, _xout = forced_response(sys, T, U, X0, transpose=transpose,
-                                     squeeze=squeeze)
-
-    if return_x:
-        return T, yout, _xout
-
-    return T, yout
+    return forced_response(sys, T, U, X0, transpose=transpose,
+                           return_x=return_x, squeeze=squeeze)
 
 
-def impulse_response(sys, T=None, X0=0., input=0, output=None, T_num=None,
-                     transpose=False, return_x=False, squeeze=True):
+def impulse_response(sys, T=None, X0=0., input=None, output=None, T_num=None,
+                     transpose=False, return_x=False, squeeze=None):
     # pylint: disable=W0622
     """Impulse response of a linear system
 
@@ -746,7 +790,7 @@ def impulse_response(sys, T=None, X0=0., input=0, output=None, T_num=None,
         Numbers are converted to constant arrays with the correct shape.
 
     input : int
-        Index of the input that will be used in this simulation.
+        Index of the input that will be used in this simulation.  Default = 0.
 
     output : int
         Index of the output that will be used in this simulation. Set to None
@@ -763,10 +807,11 @@ def impulse_response(sys, T=None, X0=0., input=0, output=None, T_num=None,
     return_x : bool
         If True, return the state vector (default = False).
 
-    squeeze : bool, optional (default=True)
-        If True, remove single-dimensional entries from the shape of
-        the output.  For single output systems, this converts the
-        output response to a 1D array.
+    squeeze : bool, optional
+        If True (default), remove single-dimensional entries from the shape
+        of the output.  For single output systems, this converts the output
+        response to a 1D array.  The default value can be set using
+        config.defaults['control.squeeze_time_response'].
 
     Returns
     -------
@@ -774,8 +819,8 @@ def impulse_response(sys, T=None, X0=0., input=0, output=None, T_num=None,
         Time values of the output
     yout : array
         Response of the system
-    xout : array
-        Individual response of each x variable
+    xout : array, optional
+        Individual response of each x variable (if return_x is True)
 
     See Also
     --------
@@ -792,7 +837,7 @@ def impulse_response(sys, T=None, X0=0., input=0, output=None, T_num=None,
     >>> T, yout = impulse_response(sys, T, X0)
 
     """
-    sys = _get_ss_simo(sys, input, output)
+    squeeze, sys = _get_ss_simo(sys, input, output, squeeze=squeeze)
 
     # if system has direct feedthrough, can't simulate impulse response
     # numerically
@@ -824,13 +869,9 @@ def impulse_response(sys, T=None, X0=0., input=0, output=None, T_num=None,
         new_X0 = X0
         U[0] = 1./sys.dt # unit area impulse
 
-    T, yout, _xout = forced_response(sys, T, U, new_X0, transpose=transpose,
-                                     squeeze=squeeze)
-
-    if return_x:
-        return T, yout, _xout
+    return forced_response(sys, T, U, new_X0, transpose=transpose,
+                           return_x=return_x, squeeze=squeeze)
 
-    return T, yout
 
 # utility function to find time period and time increment using pole locations
 def _ideal_tfinal_and_dt(sys, is_step=True):

From 05e6fe6797af7ae699ddaed7b9c174f91d779ef8 Mon Sep 17 00:00:00 2001
From: Richard Murray <murray@cds.caltech.edu>
Date: Sat, 16 Jan 2021 21:48:16 -0800
Subject: [PATCH 2/3] get rid of warnings in unit tests

---
 control/tests/matlab_test.py   |  3 ++-
 control/tests/timeresp_test.py | 34 +++++++++++++++++-----------------
 control/timeresp.py            |  4 ++--
 3 files changed, 21 insertions(+), 20 deletions(-)

diff --git a/control/tests/matlab_test.py b/control/tests/matlab_test.py
index 3a15a5aff..c9ba818cb 100644
--- a/control/tests/matlab_test.py
+++ b/control/tests/matlab_test.py
@@ -321,7 +321,8 @@ def testLsim(self, siso):
         yout, tout, _xout = lsim(siso.ss1, u, t)
         np.testing.assert_array_almost_equal(yout, youttrue, decimal=4)
         np.testing.assert_array_almost_equal(tout, t)
-        yout, _t, _xout = lsim(siso.tf3, u, t)
+        with pytest.warns(UserWarning, match="Internal conversion"):
+            yout, _t, _xout = lsim(siso.tf3, u, t)
         np.testing.assert_array_almost_equal(yout, youttrue, decimal=4)
 
         # test with initial value and special algorithm for ``U=0``
diff --git a/control/tests/timeresp_test.py b/control/tests/timeresp_test.py
index 41bf05c4b..b52c9dad7 100644
--- a/control/tests/timeresp_test.py
+++ b/control/tests/timeresp_test.py
@@ -723,26 +723,26 @@ def test_time_series_data_convention_2D(self, siso_ss1):
         assert y.ndim == 1  # SISO returns "scalar" output
         assert t.shape == y.shape  # Allows direct plotting of output
 
-    @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)],
 #       [4, 1, 1, 'siso', (8,)],        # Use for later 'siso' implementation
-        [1, 2, 1, None, (2, 8)],
-        [2, 2, 1, True, (2, 8)],
-        [3, 2, 1, False, (2, 8)],
-#       [4, 2, 1, 'siso', (2, 8)],      # Use for later 'siso' implementation
-        [1, 1, 2, None, (8,)],
-        [2, 1, 2, True, (8,)],
-        [3, 1, 2, False, (1, 8)],
-#       [4, 1, 2, 'siso', (1, 8)],      # Use for later 'siso' implementation
-        [1, 2, 2, None, (2, 8)],
-        [2, 2, 2, True, (2, 8)],
-        [3, 2, 2, False, (2, 8)],
-#       [4, 2, 2, 'siso', (2, 8)],      # Use for later 'siso' implementation
+        [3, 2, 1, None, (2, 8)],
+        [4, 2, 1, True, (2, 8)],
+        [5, 2, 1, False, (2, 8)],
+#       [6, 2, 1, 'siso', (2, 8)],      # Use for later 'siso' implementation
+        [3, 1, 2, None, (8,)],
+        [4, 1, 2, True, (8,)],
+        [5, 1, 2, False, (1, 8)],
+#       [6, 1, 2, 'siso', (1, 8)],      # Use for later 'siso' implementation
+        [4, 2, 2, None, (2, 8)],
+        [5, 2, 2, True, (2, 8)],
+        [6, 2, 2, False, (2, 8)],
+#       [7, 2, 2, 'siso', (2, 8)],      # Use for later 'siso' implementation
     ])
+    @pytest.mark.usefixtures("editsdefaults")
     def test_squeeze(self, fcn, nstate, nout, ninp, squeeze, shape):
         # Figure out if we have SciPy 1+
         scipy0 = StrictVersion(sp.__version__) < '1.0'
@@ -818,13 +818,13 @@ def test_squeeze(self, fcn, nstate, nout, ninp, squeeze, shape):
             _, yvec = ct.step_response(sys, tvec)
         assert yvec.shape == (sys.outputs, 8)
 
-        _, yvec, xvec = ct.forced_response(
-            sys, tvec, uvec, 0, return_x=True)
-        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:
-            assert xvec.shape[1] == 8
+            _, yvec = ct.forced_response(sys, tvec, uvec, 0)
+        assert yvec.shape == (sys.outputs, 8)
 
         # For InputOutputSystems, also test input_output_response
         if isinstance(sys, ct.InputOutputSystem) and not scipy0:
diff --git a/control/timeresp.py b/control/timeresp.py
index 7a48a5c5f..72013e41f 100644
--- a/control/timeresp.py
+++ b/control/timeresp.py
@@ -285,10 +285,10 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False,
         return_x = config.defaults['forced_response.return_x']
 
     # If return_x is used for TransferFunction, issue a warning
-    if return_x and not isinstance(sys, TransferFunction):
+    if return_x and isinstance(sys, TransferFunction):
         warnings.warn(
             "return_x specified for a transfer function system. Internal "
-            "conversation to state space used; results may meaningless.")
+            "conversion to state space used; results may meaningless.")
 
     sys = _convert_to_statespace(sys)
     A, B, C, D = np.asarray(sys.A), np.asarray(sys.B), np.asarray(sys.C), \

From 47bbf169b3c291b7fb4c589bd6af0b8dc2ec99df Mon Sep 17 00:00:00 2001
From: Richard Murray <murray@cds.caltech.edu>
Date: Sat, 16 Jan 2021 23:17:44 -0800
Subject: [PATCH 3/3] update squeeze processing for time responses + doc
 updates

---
 control/config.py              |   4 ++
 control/iosys.py               |  18 ++++--
 control/tests/iosys_test.py    |   2 +-
 control/tests/timeresp_test.py |  39 +++++++++---
 control/timeresp.py            | 107 +++++++++++++++++++++------------
 5 files changed, 119 insertions(+), 51 deletions(-)

diff --git a/control/config.py b/control/config.py
index 6a8fb8db6..ef3184a5e 100644
--- a/control/config.py
+++ b/control/config.py
@@ -18,6 +18,7 @@
     'control.default_dt': 0,
     'control.squeeze_frequency_response': None,
     'control.squeeze_time_response': True,
+    'control.squeeze_time_response': None,
     'forced_response.return_x': False,
 }
 defaults = dict(_control_defaults)
@@ -236,4 +237,7 @@ def use_legacy_defaults(version):
         # forced_response no longer returns x by default
         set_defaults('forced_response', return_x=True)
 
+        # time responses are only squeezed if SISO
+        set_defaults('control', squeeze_time_response=True)
+
     return (major, minor, patch)
diff --git a/control/iosys.py b/control/iosys.py
index 8fc17c016..d16cbce93 100644
--- a/control/iosys.py
+++ b/control/iosys.py
@@ -450,6 +450,10 @@ def find_state(self, name):
         """Find the index for a state given its name (`None` if not found)"""
         return self.state_index.get(name, None)
 
+    def issiso(self):
+        """Check to see if a system is single input, single output"""
+        return self.ninputs == 1 and self.noutputs == 1
+
     def feedback(self, other=1, sign=-1, params={}):
         """Feedback interconnection between two input/output systems
 
@@ -1373,18 +1377,22 @@ def input_output_response(sys, T, U=0., X0=0, params={}, method='RK45',
     return_x : bool, optional
         If True, return the values of the state at each time (default = False).
     squeeze : bool, optional
-        If True and if the system has a single output, return the
-        system output as a 1D array rather than a 2D array.  Default
-        value (True) set by config.defaults['control.squeeze_time_response'].
+        If True and if the system has a single output, return the system
+        output as a 1D array rather than a 2D array.  If False, return the
+        system output as a 2D array even if the system is SISO.  Default value
+        set by config.defaults['control.squeeze_time_response'].
 
     Returns
     -------
     T : array
         Time values of the output.
     yout : array
-        Response of the system.
+        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 2D (indexed by the output number and
+        time).
     xout : array
-        Time evolution of the state vector (if return_x=True)
+        Time evolution of the state vector (if return_x=True).
 
     Raises
     ------
diff --git a/control/tests/iosys_test.py b/control/tests/iosys_test.py
index 32e5b102f..acdf43422 100644
--- a/control/tests/iosys_test.py
+++ b/control/tests/iosys_test.py
@@ -158,7 +158,7 @@ def test_nonlinear_iosys(self, tsys):
         nlout = lambda t, x, u, params: \
             np.reshape(np.dot(linsys.C, np.reshape(x, (-1, 1)))
                        + np.dot(linsys.D, u), (-1,))
-        nlsys = ios.NonlinearIOSystem(nlupd, nlout)
+        nlsys = ios.NonlinearIOSystem(nlupd, nlout, inputs=1, outputs=1)
 
         # Make sure that simulations also line up
         T, U, X0 = tsys.T, tsys.U, tsys.X0
diff --git a/control/tests/timeresp_test.py b/control/tests/timeresp_test.py
index b52c9dad7..e5c7471b6 100644
--- a/control/tests/timeresp_test.py
+++ b/control/tests/timeresp_test.py
@@ -723,26 +723,22 @@ def test_time_series_data_convention_2D(self, siso_ss1):
         assert y.ndim == 1  # SISO returns "scalar" output
         assert t.shape == y.shape  # Allows direct plotting of output
 
+    @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)],
-#       [4, 1, 1, 'siso', (8,)],        # Use for later 'siso' implementation
         [3, 2, 1, None, (2, 8)],
         [4, 2, 1, True, (2, 8)],
         [5, 2, 1, False, (2, 8)],
-#       [6, 2, 1, 'siso', (2, 8)],      # Use for later 'siso' implementation
-        [3, 1, 2, None, (8,)],
+        [3, 1, 2, None, (1, 8)],
         [4, 1, 2, True, (8,)],
         [5, 1, 2, False, (1, 8)],
-#       [6, 1, 2, 'siso', (1, 8)],      # Use for later 'siso' implementation
         [4, 2, 2, None, (2, 8)],
         [5, 2, 2, True, (2, 8)],
         [6, 2, 2, False, (2, 8)],
-#       [7, 2, 2, 'siso', (2, 8)],      # Use for later 'siso' implementation
     ])
-    @pytest.mark.usefixtures("editsdefaults")
     def test_squeeze(self, fcn, nstate, nout, ninp, squeeze, shape):
         # Figure out if we have SciPy 1+
         scipy0 = StrictVersion(sp.__version__) < '1.0'
@@ -789,7 +785,6 @@ def test_squeeze(self, fcn, nstate, nout, ninp, squeeze, shape):
         _, 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 or (squeeze == 'siso' and not sys.issiso()):
         if squeeze is False:
             # Shape should be unsqueezed
             assert yvec.shape == (1, 8)
@@ -836,3 +831,33 @@ def test_squeeze_exception(self, fcn):
         sys = fcn(ct.rss(2, 1, 1))
         with pytest.raises(ValueError, match="unknown squeeze value"):
             step_response(sys, squeeze=1)
+
+    @pytest.mark.usefixtures("editsdefaults")
+    @pytest.mark.parametrize("nstate, nout, ninp, squeeze, shape", [
+        [1, 1, 1, None, (8,)],
+        [2, 1, 1, True, (8,)],
+        [3, 1, 1, False, (1, 8)],
+        [1, 2, 1, None, (2, 8)],
+        [2, 2, 1, True, (2, 8)],
+        [3, 2, 1, False, (2, 8)],
+        [1, 1, 2, None, (8,)],
+        [2, 1, 2, True, (8,)],
+        [3, 1, 2, False, (1, 8)],
+        [1, 2, 2, None, (2, 8)],
+        [2, 2, 2, True, (2, 8)],
+        [3, 2, 2, False, (2, 8)],
+    ])
+    def test_squeeze_0_8_4(self, nstate, nout, ninp, squeeze, shape):
+        # Set defaults to match release 0.8.4
+        ct.config.use_legacy_defaults('0.8.4')
+        ct.config.use_numpy_matrix(False)
+
+        # Generate system, time, and input vectors
+        sys = ct.rss(nstate, nout, ninp, strictly_proper=True)
+        tvec = np.linspace(0, 1, 8)
+        uvec = np.dot(
+            np.ones((sys.inputs, 1)),
+            np.reshape(np.sin(tvec), (1, 8)))
+
+        _, yvec = ct.initial_response(sys, tvec, 1, squeeze=squeeze)
+        assert yvec.shape == shape
diff --git a/control/timeresp.py b/control/timeresp.py
index 72013e41f..b8dc6631a 100644
--- a/control/timeresp.py
+++ b/control/timeresp.py
@@ -242,19 +242,27 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False,
         return only the time and output vectors.
 
     squeeze : bool, optional
-        If True (default), remove single-dimensional entries from the shape
-        of the output.  For single output systems, this converts the output
-        response to a 1D array.  The default value can be set using
+        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 2D array (indexed by the output number and time)
+        even if the system is SISO. The default value can be set using
         config.defaults['control.squeeze_time_response'].
 
     Returns
     -------
     T : array
         Time values of the output.
+
     yout : array
-        Response of the system.
+        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 2D (indexed by the output number and
+        time).
+
     xout : array
-        Time evolution of the state vector.
+        Time evolution of the state vector. Not affected by squeeze.
 
     See Also
     --------
@@ -459,11 +467,9 @@ def _process_time_response(sys, tout, yout, xout, transpose=None,
         yout = np.squeeze(yout)
     elif squeeze is False:      # squeeze no dimensions
         pass
-    # Possible code if we implement a squeeze='siso' option
-    # elif squeeze == 'siso':     # squeeze signals if SISO
-    #     yout = yout[0] if sys.issiso() else yout
+    elif squeeze is None:       # squeeze signals if SISO
+        yout = yout[0] if sys.issiso() else yout
     else:
-        # In preparation for more complicated values for squeeze
         raise ValueError("unknown squeeze value")
 
     # See if we need to transpose the data back into MATLAB form
@@ -481,13 +487,17 @@ def _get_ss_simo(sys, input=None, output=None, squeeze=None):
 
     If input is not specified, select first input and issue warning
     """
+    # If squeeze was not specified, figure out the default
+    if squeeze is None:
+        squeeze = config.defaults['control.squeeze_time_response']
+
     sys_ss = _convert_to_statespace(sys)
     if sys_ss.issiso():
         return squeeze, sys_ss
-    # Possible code if we implement a squeeze='siso' option
-    # elif squeeze == 'siso':
-    #     # Don't squeeze outputs if resulting system turns out to be siso
-    #     squeeze = False
+    elif squeeze == None and (input is None or output is None):
+        # Don't squeeze outputs if resulting system turns out to be siso
+        # Note: if we expand input to allow a tuple, need to update this check
+        squeeze = False
 
     warn = False
     if input is None:
@@ -531,14 +541,13 @@ def step_response(sys, T=None, X0=0., input=None, output=None, T_num=None,
         many simulation steps.
 
     X0 : array_like or float, optional
-        Initial condition (default = 0)
-
-        Numbers are converted to constant arrays with the correct shape.
+        Initial condition (default = 0). Numbers are converted to constant
+        arrays with the correct shape.
 
-    input : int
+    input : int, optional
         Index of the input that will be used in this simulation.  Default = 0.
 
-    output : int
+    output : int, optional
         Index of the output that will be used in this simulation. Set to None
         to not trim outputs
 
@@ -546,17 +555,20 @@ def step_response(sys, T=None, X0=0., input=None, output=None, T_num=None,
         Number of time steps to use in simulation if T is not provided as an
         array (autocomputed if not given); ignored if sys is discrete-time.
 
-    transpose : bool
+    transpose : bool, optional
         If True, transpose all input and output arrays (for backward
         compatibility with MATLAB and :func:`scipy.signal.lsim`)
 
-    return_x : bool
+    return_x : bool, optional
         If True, return the state vector (default = False).
 
     squeeze : bool, optional
-        If True (default), remove single-dimensional entries from the shape
-        of the output.  For single output systems, this converts the output
-        response to a 1D array.  The default value can be set using
+        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 2D array (indexed by the output number and time) even if
+        the system is SISO. The default value can be set using
         config.defaults['control.squeeze_time_response'].
 
     Returns
@@ -565,10 +577,13 @@ def step_response(sys, T=None, X0=0., input=None, output=None, T_num=None,
         Time values of the output
 
     yout : array
-        Response of the system
+        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 2D (indexed by the output number and
+        time).
 
     xout : array, optional
-        Individual response of each x variable (if return_x is True)
+        Individual response of each x variable (if return_x is True).
 
     See Also
     --------
@@ -722,19 +737,27 @@ def initial_response(sys, T=None, X0=0., input=0, output=None, T_num=None,
         If True, return the state vector (default = False).
 
     squeeze : bool, optional
-        If True (default), remove single-dimensional entries from the shape
-        of the output.  For single output systems, this converts the output
-        response to a 1D array.  The default value can be set using
+        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 2D array (indexed by the output number and time) even if
+        the system is SISO. The default value can be set using
         config.defaults['control.squeeze_time_response'].
 
     Returns
     -------
     T : array
         Time values of the output
+
     yout : array
-        Response of the system
+        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 2D (indexed by the output number and
+        time).
+
     xout : array, optional
-        Individual response of each x variable (if return_x is True)
+        Individual response of each x variable (if return_x is True).
 
     See Also
     --------
@@ -789,10 +812,10 @@ def impulse_response(sys, T=None, X0=0., input=None, output=None, T_num=None,
 
         Numbers are converted to constant arrays with the correct shape.
 
-    input : int
+    input : int, optional
         Index of the input that will be used in this simulation.  Default = 0.
 
-    output : int
+    output : int, optional
         Index of the output that will be used in this simulation. Set to None
         to not trim outputs
 
@@ -804,23 +827,31 @@ def impulse_response(sys, T=None, X0=0., input=None, output=None, T_num=None,
         If True, transpose all input and output arrays (for backward
         compatibility with MATLAB and :func:`scipy.signal.lsim`)
 
-    return_x : bool
+    return_x : bool, optional
         If True, return the state vector (default = False).
 
     squeeze : bool, optional
-        If True (default), remove single-dimensional entries from the shape
-        of the output.  For single output systems, this converts the output
-        response to a 1D array.  The default value can be set using
+        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 2D array (indexed by the output number and time) even if
+        the system is SISO. The default value can be set using
         config.defaults['control.squeeze_time_response'].
 
     Returns
     -------
     T : array
         Time values of the output
+
     yout : array
-        Response of the system
+        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 2D (indexed by the output number and
+        time).
+
     xout : array, optional
-        Individual response of each x variable (if return_x is True)
+        Individual response of each x variable (if return_x is True).
 
     See Also
     --------