Change default value of dt to None, add equally spaced check for response time data

KybernetikJo · KybernetikJo · commit 8d62775303e1 · 2024-07-15T12:08:30.000+02:00
diff --git a/control/modelsimp.py b/control/modelsimp.py
@@ -402,7 +402,7 @@ def era(YY, m, n, nin, nout, r):
     raise NotImplementedError('This function is not implemented yet.')
 
 
-def markov(*args, m=None, transpose=False, dt=True, truncate=False):
+def markov(*args, m=None, transpose=False, dt=None, truncate=False):
     """markov(Y, U, [, m])
     
     Calculate the first `m` Markov parameters [D CB CAB ...] from data.
@@ -448,7 +448,9 @@ def markov(*args, m=None, transpose=False, dt=True, truncate=False):
         True indicates discrete time with unspecified sampling time and a
         positive float is discrete time with the specified sampling time.
         It can be used to scale the Markov parameters in order to match
-        the unit-area impulse response of python-control. Default is True.
+        the unit-area impulse response of python-control. Default is True
+        for array_like and dt=data.time[1]-data.time[0] for
+        TimeResponseData as input.
     truncate : bool, optional
         Do not use first m equation for least squares. Default is False.
     transpose : bool, optional
@@ -483,10 +485,16 @@ def markov(*args, m=None, transpose=False, dt=True, truncate=False):
         raise ControlArgument("not enough input arguments")
     
     if isinstance(args[0], TimeResponseData):
-        Umat = np.array(args[0].inputs, ndmin=2)
-        Ymat = np.array(args[0].outputs, ndmin=2)
-        transpose = args[0].transpose
-        if args[0].transpose and not args[0].issiso:
+        data = args[0]
+        Umat = np.array(data.inputs, ndmin=2)
+        Ymat = np.array(data.outputs, ndmin=2)
+        if dt is None:
+            dt = data.time[1] - data.time[0]
+            if not np.allclose(np.diff(data.time), dt):
+                raise ValueError("response time values must be equally "
+                                 "spaced.")
+        transpose = data.transpose
+        if data.transpose and not data.issiso:
             Umat, Ymat = np.transpose(Umat), np.transpose(Ymat)
         if len(args) == 2:
             m = args[1]
@@ -497,6 +505,8 @@ def markov(*args, m=None, transpose=False, dt=True, truncate=False):
             raise ControlArgument("not enough input arguments")
         Umat = np.array(args[1], ndmin=2)
         Ymat = np.array(args[0], ndmin=2)
+        if dt is None:
+            dt = True
         if transpose:
             Umat, Ymat = np.transpose(Umat), np.transpose(Ymat)
         if len(args) == 3:
diff --git a/control/tests/modelsimp_test.py b/control/tests/modelsimp_test.py
@@ -81,11 +81,11 @@ def testMarkovSignature(self):
         T = np.linspace(0, 10, 100)
         U = np.ones((1, 100))
         T, Y = forced_response(tf([1], [1, 0.5], True), T, U)
-        H = markov(Y, U, 4, transpose=False)
+        H = markov(Y, U, 4, dt=True)
         np.testing.assert_array_almost_equal(H[:3], Htrue[:3])
 
         response = forced_response(tf([1], [1, 0.5], True), T, U)
-        H = markov(response, 4)
+        H = markov(response, 4, dt=True)
         np.testing.assert_array_almost_equal(H[:3], Htrue[:3])
 
         # Test example from issue #395
@@ -164,14 +164,28 @@ 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)
-        response = forced_response(Hd, T, U, squeeze=True)
-        Mcomp = markov(response, m)
+
+        ir_true = impulse_response(Hd,T)
+        Mtrue_scaled = ir_true[1][:m]
+
+        T, Y = forced_response(Hd, T, U, squeeze=True)
+        Mcomp = markov(Y, U, m, dt=True)
+        Mcomp_scaled = markov(Y, U, m, dt=Ts)
 
         # Compare to results from markov()
         # experimentally determined probability to get non matching results
         # with rtot=1e-6 and atol=1e-8 due to numerical errors
         # for k=5, m=n=10: 0.015 %
         np.testing.assert_allclose(Mtrue, Mcomp, rtol=1e-6, atol=1e-8)
+        np.testing.assert_allclose(Mtrue_scaled, Mcomp_scaled, rtol=1e-6, atol=1e-8)
+
+        response = forced_response(Hd, T, U, squeeze=True)
+        Mcomp = markov(response, m, dt=True)
+        Mcomp_scaled = markov(response, m, dt=Ts)
+
+        np.testing.assert_allclose(Mtrue, Mcomp, rtol=1e-6, atol=1e-8)
+        np.testing.assert_allclose(Mtrue_scaled, Mcomp_scaled, rtol=1e-6, atol=1e-8)
+        
 
     def testModredMatchDC(self):
         #balanced realization computed in matlab for the transfer function:
