Add unit tests, change okid output for siso

KybernetikJo · KybernetikJo · commit 860acb7da77f · 2024-07-16T11:24:15.000+02:00
diff --git a/control/modelsimp.py b/control/modelsimp.py
@@ -764,5 +764,9 @@ def okid(*args, m=None, transpose=False, dt=True, truncate=False):
     H[:,:,1:] = Y[:,:,:]
     H = H/dt # scaling
 
+    # for siso return a 1D array instead of a 3D array
+    if q == 1 and p == 1:
+        H = np.squeeze(H)
+
     # Return the first m Markov parameters
     return H if not transpose else np.transpose(H)
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, 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, okid, modred
 
 
 class TestModelsimp:
@@ -111,6 +111,35 @@ 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 testOKIDSignature(self):
+
+        # Example 6.1, Applied System Identification
+        m1, k1, c1 = 1., 1., 0.01
+        A = np.array([
+            [0., 1.],
+            [-k1/m1, -c1/m1],
+        ])
+        B = np.array([[0.],[1./m1]])
+        C = np.array([[-k1/m1, -c1/m1]])
+        D = np.array([[1.]])
+        sys = StateSpace(A, B, C, D)
+        dt = 0.1
+        sysd = sys.sample(dt, method='zoh')
+
+        T = np.arange(0,200,dt)
+        U = np.random.randn(sysd.B.shape[-1], len(T))
+        response = forced_response(sysd, U=U)
+        Y = response.outputs
+
+        m = 5
+        ir_true = impulse_response(sysd, T=T)
+        Htrue = ir_true.outputs[:m+1]*dt
+        H = okid(Y, U, m, dt=True)
+
+        np.testing.assert_allclose(Htrue, H, atol=1e-1)
+
+
+
     def testModredMatchDC(self):
         #balanced realization computed in matlab for the transfer function:
         # num = [1 11 45 32], den = [1 15 60 200 60]
