fix discrete time X0, U processing

murrayrm · murrayrm · commit 994fedc63730 · 2022-02-28T21:32:03.000-08:00
diff --git a/control/iosys.py b/control/iosys.py
@@ -1820,6 +1820,7 @@ def input_output_response(
         legal_shapes = [(sys.ninputs, n_steps)]
     U = _check_convert_array(U, legal_shapes,
                              'Parameter ``U``: ', squeeze=False)
+    U = U.reshape(-1, n_steps)
 
     # Check to make sure this is not a static function
     nstates = _find_size(sys.nstates, X0)
@@ -1908,7 +1909,7 @@ def ivp_rhs(t, x):
         # Compute the solution
         soln = sp.optimize.OptimizeResult()
         soln.t = T                      # Store the time vector directly
-        x = [float(x0) for x0 in X0]    # State vector (store as floats)
+        x = X0                          # Initilize state
         soln.y = []                     # Solution, following scipy convention
         y = []                          # System output
         for i in range(len(T)):
diff --git a/control/tests/iosys_test.py b/control/tests/iosys_test.py
@@ -731,6 +731,32 @@ def test_discrete(self, tsys):
         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.)
 
+    def test_discrete_iosys(self, tsys):
+        """Create a discrete time system from scratch"""
+        linsys = ct.StateSpace(
+            [[-1, 1], [0, -2]], [[0], [1]], [[1, 0]], [[0]], True)
+
+        # Create nonlinear version of the same system
+        def nlsys_update(t, x, u, params):
+            A, B = params['A'], params['B']
+            return A @ x + B @ u
+        def nlsys_output(t, x, u, params):
+            C = params['C']
+            return C @ x
+        nlsys = ct.NonlinearIOSystem(
+            nlsys_update, nlsys_output, inputs=1, outputs=1, states=2, dt=True)
+        
+        # Set up parameters for simulation
+        T, U, X0 = tsys.T, tsys.U, tsys.X0
+
+        # Simulate and compare to LTI output
+        ios_t, ios_y = ios.input_output_response(
+            nlsys, T, U, X0,
+            params={'A': linsys.A, 'B': linsys.B, 'C': linsys.C})
+        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.)
+
     def test_find_eqpts(self, tsys):
         """Test find_eqpt function"""
         # Simple equilibrium point with no inputs
@@ -1526,7 +1552,6 @@ def secord_update(t, x, u, params={}):
     """Second order system dynamics"""
     omega0 = params.get('omega0', 1.)
     zeta = params.get('zeta', 0.5)
-    u = np.array(u, ndmin=1)
     return np.array([
         x[1],
         -2 * zeta * omega0 * x[1] - omega0*omega0 * x[0] + u[0]
