Change API to work with TimeResponseData and ndarray as input, add pytest, fix small things

KybernetikJo · KybernetikJo · commit 614a0808a452 · 2024-07-07T19:31:31.000+02:00
diff --git a/control/modelsimp.py b/control/modelsimp.py
@@ -48,6 +48,7 @@
 from .iosys import isdtime, isctime
 from .statesp import StateSpace
 from .statefbk import gram
+from .timeresp import TimeResponseData
 
 __all__ = ['hsvd', 'balred', 'modred', 'era', 'markov', 'minreal']
 
@@ -368,8 +369,10 @@ def minreal(sys, tol=None, verbose=True):
     return sysr
 
 
-def era(data, r, m=None, n=None, dt=True):
-    """Calculate an ERA model of order `r` based on the impulse-response data.
+def era(arg, r, m=None, n=None, dt=True, transpose=False):
+    r"""era(YY, r)
+
+    Calculate an ERA model of order `r` based on the impulse-response data.
 
     This function computes a discrete time system
 
@@ -380,8 +383,19 @@ def era(data, r, m=None, n=None, dt=True):
 
     for a given impulse-response data (see [1]_).
 
+    The function can be called with 2 arguments:
+
+    * ``sysd, S = era(data, r)``
+    * ``sysd, S = era(YY, r)``
+
+    where `response` is an `TimeResponseData` object, and `YY` is 1D or 3D
+    array and r is an integer.
+
     Parameters
     ----------
+    YY : array_like
+        impulse-response data from which the StateSpace model is estimated,
+        1D or 3D array.
     data : TimeResponseData
         impulse-response data from which the StateSpace model is estimated.
     r : integer
@@ -398,11 +412,16 @@ def era(data, r, m=None, n=None, dt=True):
         It can be used to scale the StateSpace model in order to match
         the impulse response of this library.
         Default values is True.
+    transpose : bool, optional
+        Assume that input data is transposed relative to the standard
+        :ref:`time-series-convention`. For TimeResponseData this parameter 
+        is ignored.
+        Default value is False.
 
     Returns
     -------
     sys : StateSpace
-        A reduced order model sys=StateSpace(Ar,Br,Cr,Dr,dt)
+        A reduced order model sys=StateSpace(Ar,Br,Cr,Dr,dt).
     S : array
         Singular values of Hankel matrix.
         Can be used to choose a good r value.
@@ -416,6 +435,10 @@ def era(data, r, m=None, n=None, dt=True):
 
     Examples
     --------
+    >>> T = np.linspace(0, 10, 100)
+    >>> _, YY = ct.impulse_response(ct.tf([1], [1, 0.5], True), T)
+    >>> sysd, _ = ct.era(YY, r=1)
+
     >>> T = np.linspace(0, 10, 100)
     >>> response = ct.impulse_response(ct.tf([1], [1, 0.5], True), T)
     >>> sysd, _ = ct.era(response, r=1)
@@ -434,10 +457,16 @@ def block_hankel_matrix(Y, m, n):
                 
         return H
     
-    Y = np.array(data.outputs, ndmin=3)
-    if data.transpose:
-        Y = np.transpose(Y)
-    q, p, l = Y.shape
+    if isinstance(arg, TimeResponseData):
+        YY = np.array(arg.outputs, ndmin=3)
+        if arg.transpose:
+            YY = np.transpose(YY)
+    else:
+        YY = np.array(arg, ndmin=3)
+        if transpose:
+            YY = np.transpose(YY)
+    
+    q, p, l = YY.shape
 
     if m is None:
         m = 2*r
@@ -450,7 +479,7 @@ def block_hankel_matrix(Y, m, n):
     if (l-1) < m+n:
         raise ValueError("Not enough data for requested number of parameters")
     
-    H = block_hankel_matrix(Y[:,:,1:], m, n+1) # Hankel matrix (q*m, p*(n+1))
+    H = block_hankel_matrix(YY[:,:,1:], m, n+1) # Hankel matrix (q*m, p*(n+1))
     Hf = H[:,:-p] # first p*n columns of H
     Hl = H[:,p:] # last p*n columns of H
     
@@ -463,7 +492,7 @@ def block_hankel_matrix(Y, m, n):
     Ar = Sigma_inv @ Ur.T @ Hl @ Vhr.T @ Sigma_inv
     Br = Sigma_inv @ Ur.T @ Hf[:,0:p]*dt
     Cr = Hf[0:q,:] @ Vhr.T @ Sigma_inv
-    Dr = Y[:,:,0]
+    Dr = YY[:,:,0]
 
     return StateSpace(Ar,Br,Cr,Dr,dt), S
 
diff --git a/control/tests/modelsimp_test.py b/control/tests/modelsimp_test.py
@@ -7,10 +7,10 @@
 import pytest
 
 
-from control import StateSpace, forced_response, tf, rss, c2d
+from control import StateSpace, impulse_response, step_response, forced_response, tf, rss, c2d
 from control.exception import ControlMIMONotImplemented
 from control.tests.conftest import slycotonly
-from control.modelsimp import balred, hsvd, markov, modred
+from control.modelsimp import balred, hsvd, markov, modred, era
 
 
 class TestModelsimp:
@@ -111,6 +111,85 @@ def testMarkovResults(self, k, m, n):
         # for k=5, m=n=10: 0.015 %
         np.testing.assert_allclose(Mtrue, Mcomp, rtol=1e-6, atol=1e-8)
 
+    def testERASignature(self):
+
+        # test siso
+        # Katayama, Subspace Methods for System Identification
+        # Example 6.1, Fibonacci sequence
+        H_true = np.array([0.,1.,1.,2.,3.,5.,8.,13.,21.,34.])
+
+        # A realization of fibonacci impulse response
+        A = np.array([[0., 1.],[1., 1.,]])
+        B = np.array([[1.],[1.,]])
+        C = np.array([[1., 0.,]])
+        D = np.array([[0.,]])
+        
+        T = np.arange(0,10,1)
+        sysd_true = StateSpace(A,B,C,D,True)
+        ir_true = impulse_response(sysd_true,T=T)
+
+        # test TimeResponseData
+        sysd_est, _  = era(ir_true,r=2)
+        ir_est = impulse_response(sysd_est, T=T)
+        _, H_est = ir_est
+    
+        np.testing.assert_allclose(H_true, H_est, rtol=1e-6, atol=1e-8)
+
+        # test ndarray
+        _, YY_true = ir_true
+        sysd_est, _  = era(YY_true,r=2)
+        ir_est = impulse_response(sysd_est, T=T)
+        _, H_est = ir_est
+    
+        np.testing.assert_allclose(H_true, H_est, rtol=1e-6, atol=1e-8)
+
+        # test mimo
+        # 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.
+        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.1
+        T = np.arange(0,10,dt)
+        sysd_true = sys.sample(dt, method='zoh')
+        ir_true = impulse_response(sysd_true, T=T)
+
+        # test TimeResponseData
+        sysd_est, _ = era(ir_true,r=4,dt=dt)
+
+        step_true = step_response(sysd_true)
+        step_est = step_response(sysd_est)
+
+        np.testing.assert_allclose(step_true.outputs, 
+                                   step_est.outputs,
+                                   rtol=1e-6, atol=1e-8)
+        
+        # test ndarray
+        _, YY_true = ir_true
+        sysd_est, _  = era(YY_true,r=4,dt=dt)
+
+        step_true = step_response(sysd_true, T=T)
+        step_est = step_response(sysd_est, T=T)
+
+        np.testing.assert_allclose(step_true.outputs, 
+                                   step_est.outputs,
+                                   rtol=1e-6, atol=1e-8)
+
+
     def testModredMatchDC(self):
         #balanced realization computed in matlab for the transfer function:
         # num = [1 11 45 32], den = [1 15 60 200 60]
diff --git a/examples/era_msd.py b/examples/era_msd.py
@@ -11,10 +11,12 @@
 import control as ct
 
 # 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
-m2, k2, c2 = 2., .5, .1
-k3, c3 = .5, .1
+# 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.
+m2, k2, c2 = 2., 2., 1.
+k3, c3 = 6., 2.
 
 A = np.array([
     [0., 0., 1., 0.],
@@ -39,7 +41,7 @@
         sys = ct.StateSpace(A, B, C, D)
 
 
-dt = 0.5
+dt = 0.1
 sysd = sys.sample(dt, method='zoh')
 response = ct.impulse_response(sysd)
 response.plot()
@@ -48,7 +50,9 @@
 sysd_est, _ = ct.era(response,r=4,dt=dt)
 
 step_true = ct.step_response(sysd)
+step_true.sysname="H_true"
 step_est = ct.step_response(sysd_est)
+step_est.sysname="H_est"
 
 step_true.plot(title=xixo)
 step_est.plot(color='orange',linestyle='dashed')
