Change API to work with ndarray and TimeResponseData as input

KybernetikJo · KybernetikJo · commit 2b1cd219431a · 2024-07-07T16:07:01.000+02:00
diff --git a/control/modelsimp.py b/control/modelsimp.py
@@ -43,8 +43,7 @@
 # External packages and modules
 import numpy as np
 import warnings
-from .exception import ControlSlycot, ControlMIMONotImplemented, \
-    ControlDimension
+from .exception import ControlSlycot, ControlArgument, ControlDimension
 from .iosys import isdtime, isctime
 from .statesp import StateSpace
 from .statefbk import gram
@@ -403,8 +402,10 @@ def era(YY, m, n, nin, nout, r):
     raise NotImplementedError('This function is not implemented yet.')
 
 
-def markov(data, m=None, dt=True, truncate=False):
-    """Calculate the first `m` Markov parameters [D CB CAB ...]
+def markov(*args, **kwargs):
+    """markov(Y, U, [, m])
+    
+    Calculate the first `m` Markov parameters [D CB CAB ...]
     from data
 
     This function computes the Markov parameters for a discrete time system
@@ -419,14 +420,31 @@ def markov(data, m=None, dt=True, truncate=False):
     the input data is less than the desired number of Markov parameters (a
     warning message is generated in this case).
 
+    The function can be called with either 1, 2, or 3 arguments:
+
+    * ``K, S, E = lqr(response)``
+    * ``K, S, E = lqr(respnose, m)``
+    * ``K, S, E = lqr(Y, U)``
+    * ``K, S, E = lqr(Y, U, m)``
+
+    where `response` is an `TimeResponseData` object, and `Y`, `U`, are 1D or 2D
+    array and m is an integer.
+
     Parameters
     ----------
+    Y : array_like
+        Output data.  If the array is 1D, the system is assumed to be single
+        input.  If the array is 2D and transpose=False, the columns of `Y`
+        are taken as time points, otherwise the rows of `Y` are taken as
+        time points.
+    U : array_like
+        Input data, arranged in the same way as `Y`.
     data : TimeResponseData
         Response data from which the Markov parameters where estimated.
         Input and output data must be 1D or 2D array.
     m : int, optional
         Number of Markov parameters to output.  Defaults to len(U).
-    dt : (True of float, optional)
+    dt : True of float, optional
         True indicates discrete time with unspecified sampling time,
         positive number is discrete time with specified sampling time.
         It can be used to scale the markov parameters in order to match
@@ -460,17 +478,41 @@ def markov(data, m=None, dt=True, truncate=False):
     --------
     >>> T = np.linspace(0, 10, 100)
     >>> U = np.ones((1, 100))
-    >>> response = ct.forced_response(ct.tf([1], [1, 0.5], True), T, U)
-    >>> H = ct.markov(response, 3)
+    >>> T, Y = ct.forced_response(ct.tf([1], [1, 0.5], True), T, U)
+    >>> H = ct.markov(Y, U, 3, transpose=False)
 
     """
     # Convert input parameters to 2D arrays (if they aren't already)
-    Umat = np.array(data.inputs, ndmin=2)
-    Ymat = np.array(data.outputs, ndmin=2)
 
-    # If data is in transposed format, switch it around
-    if data.transpose and not data.issiso:
-        Umat, Ymat = np.transpose(Umat), np.transpose(Ymat)
+    # Get the system description
+    if (len(args) < 1):
+        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:
+            Umat, Ymat = np.transpose(Umat), np.transpose(Ymat)
+        index = 1
+    else:
+        if (len(args) < 2):
+            raise ControlArgument("not enough input arguments")
+        Umat = np.array(args[0], ndmin=2)
+        Ymat = np.array(args[1], ndmin=2)
+        transpose = kwargs.pop('transpose', False)
+        if transpose:
+            Umat, Ymat = np.transpose(Umat), np.transpose(Ymat)
+        index = 2
+
+
+    if (len(args) > index):
+        m = args[index]
+    else:
+        m = None
+
+    dt = kwargs.pop('dt', True)
+    truncate = kwargs.pop('truncate', False)
 
     # Make sure the number of time points match
     if Umat.shape[1] != Ymat.shape[1]:
@@ -549,4 +591,4 @@ def markov(data, m=None, dt=True, truncate=False):
         H = np.squeeze(H)
 
     # Return the first m Markov parameters
-    return H if not data.transpose else np.transpose(H)
+    return H if not transpose else np.transpose(H)
diff --git a/control/tests/kwargs_test.py b/control/tests/kwargs_test.py
@@ -252,6 +252,7 @@ def test_response_plot_kwargs(data_fcn, plot_fcn, mimo):
     'linearize': test_unrecognized_kwargs,
     'lqe': test_unrecognized_kwargs,
     'lqr': test_unrecognized_kwargs,
+    'markov': test_unrecognized_kwargs,
     'nichols_plot': test_matplotlib_kwargs,
     'nichols': test_matplotlib_kwargs,
     'nlsys': test_unrecognized_kwargs,
diff --git a/control/tests/modelsimp_test.py b/control/tests/modelsimp_test.py
@@ -7,8 +7,7 @@
 import pytest
 
 
-from control import StateSpace, forced_response, tf, rss, c2d, TimeResponseData
-from control.exception import ControlMIMONotImplemented
+from control import StateSpace, forced_response, impulse_response, tf, rss, c2d, TimeResponseData
 from control.tests.conftest import slycotonly
 from control.modelsimp import balred, hsvd, markov, modred
 
@@ -33,44 +32,88 @@ def testHSVD(self):
         assert not isinstance(hsv, np.matrix)
 
     def testMarkovSignature(self):
-        U = np.array([1., 1., 1., 1., 1.])
+        U = np.array([[1., 1., 1., 1., 1.]])
         Y = U
         response = TimeResponseData(time=np.arange(U.shape[-1]),
                                     outputs=Y,
                                     output_labels='y',
                                     inputs=U,
                                     input_labels='u',
                                     )
+        
+        # Basic usage
         m = 3
+        H = markov(Y, U, m, transpose=False)
+        Htrue = np.array([1., 0., 0.])
+        np.testing.assert_array_almost_equal(H, Htrue)
+
+        response.transpose=False
         H = markov(response, m)
         Htrue = np.array([1., 0., 0.])
         np.testing.assert_array_almost_equal(H, Htrue)
 
         # Make sure that transposed data also works
+        H = markov(Y.T, U.T, m, transpose=True)
+        np.testing.assert_array_almost_equal(H, np.transpose(Htrue))
+
         response.transpose=True
         H = markov(response, m)
         np.testing.assert_array_almost_equal(H, np.transpose(Htrue))
+        response.transpose=False
 
         # Generate Markov parameters without any arguments
-        response.transpose=False
+        H = markov(Y, U, m)
+        np.testing.assert_array_almost_equal(H, Htrue)
+
         H = markov(response, m)
         np.testing.assert_array_almost_equal(H, Htrue)
 
         # Test example from docstring
+        T = np.linspace(0, 10, 100)
+        U = np.ones((1, 100))
+        _, Y = forced_response(tf([1], [1, 0.5], True), T, U)
+        H = markov(Y, U, 3)
+
         T = np.linspace(0, 10, 100)
         U = np.ones((1, 100))
         response = forced_response(tf([1], [1, 0.5], True), T, U)
         H = markov(response, 3)
 
         # Test example from issue #395
-        #inp = np.array([1, 2])
-        #outp = np.array([2, 4])
-        #mrk = markov(outp, inp, 1, transpose=False)
-
-        # Make sure MIMO generates an error
-        #U = np.ones((2, 100))   # 2 inputs (Y unchanged, with 1 output)
-        #with pytest.raises(ControlMIMONotImplemented):
-        #    markov(Y, U, m)
+        inp = np.array([1, 2])
+        outp = np.array([2, 4])
+        mrk = markov(outp, inp, 1, transpose=False)
+
+        # Test mimo example
+        # Mechanical Vibrations: Theory and Application, SI Edition, 1st ed.
+        # Figure 6.5 / Example 6.7
+        m1, k1, c1 = 1., 4., 1.
+        m2, k2, c2 = 2., 2., 1.
+        k3, c3 = 6., 2.
+
+        A = np.array([
+            [0., 0., 1., 0.],
+            [0., 0., 0., 1.],
+            [-(k1+k2)/m1, (k2)/m1, -(c1+c2)/m1, c2/m1],
+            [(k2)/m2, -(k2+k3)/m2, c2/m2, -(c2+c3)/m2]
+        ])
+        B = np.array([[0.,0.],[0.,0.],[1/m1,0.],[0.,1/m2]])
+        C = np.array([[1.0, 0.0, 0.0, 0.0],[0.0, 1.0, 0.0, 0.0]])
+        D = np.zeros((2,2))
+
+        sys = StateSpace(A, B, C, D)
+        dt = 0.25
+        sysd = sys.sample(dt, method='zoh')
+
+        t = np.arange(0,100,dt)
+        u = np.random.randn(sysd.B.shape[-1], len(t))
+        response = forced_response(sysd, U=u)
+
+        m = 100
+        H = markov(response, m, dt=dt)
+        _, Htrue = impulse_response(sysd, T=dt*(m-1))
+
+        np.testing.assert_array_almost_equal(H, Htrue)
 
     # Make sure markov() returns the right answer
     @pytest.mark.parametrize("k, m, n",
diff --git a/examples/markov.py b/examples/markov.py
@@ -48,7 +48,7 @@ def create_impulse_response(H, time, transpose, dt):
                             issiso=issiso)
 
 # set up a mass spring damper system (2dof, MIMO case)
-# Mechanical Vibartions: Theory and Application, SI Edition, 1st ed.
+# Mechanical Vibrations: Theory and Application, SI Edition, 1st ed.
 # Figure 6.5 / Example 6.7
 # m q_dd + c q_d + k q = f
 m1, k1, c1 = 1., 4., 1.
@@ -91,10 +91,9 @@ def create_impulse_response(H, time, transpose, dt):
 plt.show()
 
 m = 50
-ir_true = ct.impulse_response(sysd,T=dt*m)
-ir_true.tranpose = True
+ir_true = ct.impulse_response(sysd, T=dt*m)
 
-H_est = ct.markov(response,m=m,dt=dt)
+H_est = ct.markov(response, m, dt=dt)
 # Helper function for plotting only
 ir_est = create_impulse_response(H_est,
                                  ir_true.time,
