Change output to ndarray

murrayrm · Aug 6, 2024 · Jul 3, 2024 · Jul 4, 2024 · Jul 4, 2024 · Jul 4, 2024
commit 60234806736889bf19adfb51fa77a153fd189120
diff --git a/control/modelsimp.py b/control/modelsimp.py
@@ -438,21 +438,9 @@ def markov(data, m=None, dt=True, truncate=False):
 
     Returns
     -------
-    H : TimeResponseData
-        Markov parameters / impulse response [D CB CAB ...] represented as
-        a :class:`TimeResponseData` object containing the following properties:
-
-        * time (array): Time values of the output.
-
-        * outputs (array): Response of the system.  If the system is SISO,
-          the array is 1D (indexed by time).  If the
-          system is not SISO, the array is 3D (indexed
-          by the output, trace, and time).
-
-        * inputs (array): Inputs of the system.  If the system is SISO,
-          the array is 1D (indexed by time).  If the
-          system is not SISO, the array is 3D (indexed
-          by the output, trace, and time).
+    H : ndarray
+        First m Markov parameters, [D CB CAB ...]
+
 
     Notes
     -----
@@ -557,23 +545,8 @@ def markov(data, m=None, dt=True, truncate=False):
     H = H.reshape(q,m,p) # output, time*input -> output, time, input
     H = H.transpose(0,2,1) # output, input, time
 
-    # Create unit area impulse inputs
-    inputs = np.zeros((p,p,m))
-    trace_labels, trace_types = [], []
-    for i in range(p):
-        inputs[i,i,0] = 1/dt # unit area impulse
-        trace_labels.append(f"From {data.input_labels[i]}")
-        trace_types.append('impulse')
+    if q == 1 and p == 1:
+        H = np.squeeze(H)
 
-    # Markov parameters as TimeResponseData with unit area impulse inputs
     # Return the first m Markov parameters
-    return TimeResponseData(time=data.time[:m],
-                            outputs=H,
-                            output_labels=data.output_labels,
-                            inputs=inputs,
-                            input_labels=data.input_labels,
-                            trace_labels=trace_labels,
-                            trace_types=trace_types,
-                            transpose=data.transpose,
-                            plot_inputs=False,
-                            issiso=data.issiso)
+    return H if not data.transpose else np.transpose(H)
diff --git a/control/tests/modelsimp_test.py b/control/tests/modelsimp_test.py
@@ -44,17 +44,17 @@ def testMarkovSignature(self):
         m = 3
         H = markov(response, m)
         Htrue = np.array([1., 0., 0.])
-        np.testing.assert_array_almost_equal(H.outputs, Htrue)
+        np.testing.assert_array_almost_equal(H, Htrue)
 
         # Make sure that transposed data also works
         response.transpose=True
         H = markov(response, m)
-        np.testing.assert_array_almost_equal(H.outputs, np.transpose(Htrue))
+        np.testing.assert_array_almost_equal(H, np.transpose(Htrue))
 
         # Generate Markov parameters without any arguments
         response.transpose=False
         H = markov(response, m)
-        np.testing.assert_array_almost_equal(H.outputs, Htrue)
+        np.testing.assert_array_almost_equal(H, Htrue)
 
         # Test example from docstring
         T = np.linspace(0, 10, 100)
@@ -119,7 +119,7 @@ def testMarkovResults(self, k, m, n):
         # 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.outputs, rtol=1e-6, atol=1e-8)
+        np.testing.assert_allclose(Mtrue, Mcomp, rtol=1e-6, atol=1e-8)
 
     def testModredMatchDC(self):
         #balanced realization computed in matlab for the transfer function:

diff --git a/examples/markov.py b/examples/markov.py
@@ -10,6 +10,43 @@
 
 import control as ct
 
+def create_impulse_response(H, time, transpose, dt):
+    """Helper function to use TimeResponseData type for plotting"""
+
+    H = np.array(H, ndmin=3)
+
+    if transpose:
+        H = np.transpose(H)
+
+    q, p, m = H.shape
+    inputs = np.zeros((p,p,m))
+
+    issiso = True if (q == 1 and p == 1) else False
+
+    input_labels = []
+    trace_labels, trace_types = [], []
+    for i in range(p):
+        inputs[i,i,0] = 1/dt # unit area impulse
+        input_labels.append(f"u{[i]}")
+        trace_labels.append(f"From u{[i]}")
+        trace_types.append('impulse')
+
+    output_labels = []
+    for i in range(q):
+        output_labels.append(f"y{[i]}")
+
+    return ct.TimeResponseData(time=time[:m],
+                            outputs=H,
+                            output_labels=output_labels,
+                            inputs=inputs,
+                            input_labels=input_labels,
+                            trace_labels=trace_labels,
+                            trace_types=trace_types,
+                            sysname="H_est",
+                            transpose=transpose,
+                            plot_inputs=False,
+                            issiso=issiso)
+
 # set up a mass spring damper system (2dof, MIMO case)
 # m q_dd + c q_d + k q = u
 m1, k1, c1 = 1., 1., .1
@@ -41,19 +78,29 @@
 
 dt = 0.5
 sysd = sys.sample(dt, method='zoh')
+sysd.name = "H_true"
 
-t = np.arange(0,5000,dt)
-u = np.random.randn(sysd.B.shape[-1], len(t)) # random forcing input
+ # random forcing input
+t = np.arange(0,500,dt)
+u = np.random.randn(sysd.B.shape[-1], len(t))
 
 response = ct.forced_response(sysd, U=u)
 response.plot()
 plt.show()
 
-markov_true = ct.impulse_response(sysd,T=dt*100)
-markov_est = ct.markov(response,m=100,dt=dt)
+m = 100
+ir_true = ct.impulse_response(sysd,T=dt*m)
+ir_true.tranpose = True
+
+H_est = ct.markov(response,m=m,dt=dt)
+# Helper function for plotting only
+ir_est = create_impulse_response(H_est,
+                                 ir_true.time,
+                                 ir_true.transpose,
+                                 dt)
 
-markov_true.plot(title=xixo)
-markov_est.plot(color='orange',linestyle='dashed')
+ir_true.plot(title=xixo)
+ir_est.plot(color='orange',linestyle='dashed')
 plt.show()
 
 if 'PYCONTROL_TEST_EXAMPLES' not in os.environ: