Repository: arminwitte/llsi

Issue:
SISO -> MIMO

Discussion:
@arminwitte: Resolved by #48 @arminwitte: Resolved by #48 

Solution patch:
diff --git a/src/llsi/ltimodel.py b/src/llsi/ltimodel.py
index 644e209..b42f164 100644
--- a/src/llsi/ltimodel.py
+++ b/src/llsi/ltimodel.py
@@ -12,14 +12,16 @@
 
 
 class LTIModel(ABC):
-    def __init__(self, Ts=1.0):
+    def __init__(self, Ts=1.0, nu=1, ny=1):
         self.Ts = Ts
         self.info = {}
+        self.nu = nu
+        self.ny = ny
 
     def impulse_response(self, N=100):
         t = np.linspace(0, (N - 1) * self.Ts, N)
-        u = np.zeros((N,))
-        u[0] = 1 / self.Ts
+        u = np.zeros((N, self.nu))
+        u[0, :] = 1 / self.Ts
 
         y = self.simulate(u)
 
@@ -27,7 +29,7 @@ def impulse_response(self, N=100):
 
     def step_response(self, N=100):
         t = np.linspace(0, (N - 1) * self.Ts, N)
-        u = np.ones((N,))
+        u = np.ones((N, self.nu))
 
         y = self.simulate(u)
 
diff --git a/src/llsi/statespacemodel.py b/src/llsi/statespacemodel.py
index 0b42745..abda121 100644
--- a/src/llsi/statespacemodel.py
+++ b/src/llsi/statespacemodel.py
@@ -20,7 +20,7 @@ class StateSpaceModel(LTIModel):
     https://en.wikipedia.org/wiki/State-space_representation
     """
 
-    def __init__(self, A=None, B=None, C=None, D=None, Ts=1.0, Nx=0, Nu=1, Ny=1):
+    def __init__(self, A=None, B=None, C=None, D=None, Ts=1.0, nx=0, nu=1, ny=1):
         """
 
         Parameters
@@ -48,31 +48,31 @@ def __init__(self, A=None, B=None, C=None, D=None, Ts=1.0, Nx=0, Nu=1, Ny=1):
         # set A matrix and number of states
         if A is not None:
             self.A = np.array(A)
-            self.Nx = self.A.shape[0]
+            self.nx = self.A.shape[0]
         else:
-            self.Nx = Nx
-            self.A = np.zeros((self.Nx, self.Nx))
+            self.nx = nx
+            self.A = np.zeros((self.nx, self.nx))
 
         # set B matrix and number of inputs
         if B is not None:
-            self.B = np.array(B).reshape(self.Nx, -1)
-            self.Nu = self.B.shape[1]
+            self.B = np.array(B).reshape(self.nx, -1)
+            self.nu = self.B.shape[1]
         else:
-            self.Nu = Nu
-            self.B = np.zeros((self.Nx, self.Nu))
+            self.nu = nu
+            self.B = np.zeros((self.nx, self.nu))
 
         # set C matrix and number of outputs
         if C is not None:
-            self.C = np.array(C).reshape(-1, self.Nx)
-            self.Ny = self.C.shape[0]
+            self.C = np.array(C).reshape(-1, self.nx)
+            self.ny = self.C.shape[0]
         else:
-            self.Ny = Ny
-            self.C = np.zeros((self.Ny, self.Nx))
+            self.ny = ny
+            self.C = np.zeros((self.ny, self.nx))
 
         if D is not None:
-            self.D = np.array(D).reshape(self.Ny, self.Nu)
+            self.D = np.array(D).reshape(self.ny, self.nu)
         else:
-            self.D = np.zeros((self.Ny, self.Nu))
+            self.D = np.zeros((self.ny, self.nu))
 
         self.cov = None
 
@@ -89,9 +89,9 @@ def vectorize(self):
         return np.array(theta).ravel()
 
     def reshape(self, theta: np.ndarray):
-        nx = self.Nx
-        nu = self.Nu
-        ny = self.Ny
+        nx = self.nx
+        nu = self.nu
+        ny = self.ny
 
         na = nx * nx
         nb = nx * nu
@@ -104,13 +104,17 @@ def reshape(self, theta: np.ndarray):
         self.D = theta[na + nb + nc :].reshape(ny, nu)
 
     def simulate(self, u: np.ndarray):
-        u = u.reshape(-1, self.Nu)
-        N = u.shape[0]
+        u = np.atleast_2d(u)
+        u = u.reshape(self.nu, -1)
+        N = u.shape[1]
         # TODO: initialize x properly
-        x1 = np.zeros((self.Nx, 1))
-        y = np.empty((N, self.Ny))
-        for i, u_ in enumerate(u):
-            u_ = u_.reshape(1, self.Nu)
+        x1 = np.zeros((self.nx, 1))
+        y = np.empty((N, self.ny))
+        # assert u.shape[1] == self.nu
+        for i, u_ in enumerate(u.T):
+            u_ = u_.T
+
+            u_ = u_.reshape(self.nu, 1)
             x = x1
             with np.errstate(over="ignore", invalid="ignore"):
                 x1 = self.A @ x + self.B @ u_
@@ -121,7 +125,7 @@ def simulate(self, u: np.ndarray):
             # print(f"x1:{x1}")
             # print(f"y:{y_}")
 
-            y[i, :] = y_
+            y[i, :] = y_.ravel()
 
         return y
 
@@ -138,7 +142,7 @@ def from_PT1(cls, K: float, tauC: float, Ts=1.0):
         A = [[-a]]
         C = [1]
 
-        mod = cls(A=A, B=B, C=C, D=D, Ts=Ts, Nx=1)
+        mod = cls(A=A, B=B, C=C, D=D, Ts=Ts, nx=1)
 
         return mod
 
@@ -302,9 +306,9 @@ def to_json(self, filename=None):
             data["info"] = self.info.__repr__()
         except:
             data["info"] = {}
-        data["Nx"] = self.Nx
-        data["Nu"] = self.Nu
-        data["Ny"] = self.Ny
+        data["nx"] = self.ny
+        data["nu"] = self.ny
+        data["ny"] = self.ny
 
         if filename is not None:
             with open(filename, "w") as f:
@@ -323,9 +327,9 @@ def from_json(cls, filename):
             C=data["C"],
             D=data["D"],
             Ts=data["Ts"],
-            Nx=data["Nx"],
-            Nu=data["Nu"],
-            Ny=data["Ny"],
+            nx=data["nx"],
+            nu=data["nu"],
+            ny=data["ny"],
         )
         mod.info = data["info"]
         return mod
diff --git a/src/llsi/subspace.py b/src/llsi/subspace.py
index 3676597..7ef9669 100644
--- a/src/llsi/subspace.py
+++ b/src/llsi/subspace.py
@@ -23,9 +23,9 @@ def __init__(self, data, y_name, u_name, settings):
     def hankel(x, n):
         A = []
         n = n // x.shape[1]
-        for x_ in x.T:
-            x_ = x_.ravel().T
-            for i in range(n):
+        for i in range(n):
+            for x_ in x.T:
+                x_ = x_.ravel().T
                 A.append(x_[i : -n + i])
 
         return np.array(A)
@@ -37,12 +37,12 @@ def _abcd_state(self, Xf, s, n, r):
         # Theta = Y_ @ np.linalg.pinv(X_)
 
         ###########################
-        l = self.settings.get("lambda", 0.0)
+        lmbd = self.settings.get("lambda", 0.0)
         # print(l)
         U, s, Vh = scipy.linalg.svd(X_, full_matrices=False)
         Sigma = np.diag(1 / s)
 
-        rho = np.diag(s**2 / (s**2 + l))
+        rho = np.diag(s**2 / (s**2 + lmbd))
         Theta = Y_ @ (Vh.T @ rho @ Sigma @ U.T)
         ###########################
 
@@ -53,38 +53,43 @@ def _abcd_state(self, Xf, s, n, r):
         return A, B, C, D
 
     def _abcd_observability_matrix(self, U1, U2, L11, L31, Sigma_sqrt, n, r):
+        ny = self.ny
+        nu = self.nu
+
         Or = U1 @ Sigma_sqrt  # extended observability matrix
 
-        C = Or[0, :]  # TODO: might be wrong!!!
-        A = scipy.linalg.pinv(Or[0:-1, :]) @ Or[1:, :]
+        C = Or[:ny, :]  # TODO: might be wrong!!!
+        A = scipy.linalg.pinv(Or[0:-ny, :]) @ Or[ny:, :]
         # A = scipy.linalg.lstsq(Or[0:-1,:],Or[1:,:])
         # print(A)
 
-        P = U2.T
-        # print(P)
-        A1_ = P.ravel(order="F")
+        # ====================================================================
+
+        P = np.split(U2.T, U2.T.shape[1] // ny, axis=1)
 
-        nn = A1_.shape[0]
-        ny = 1
+        nn = len(P) * P[0].shape[0]
+        rny = r
+        r = r // ny
         A_ = np.zeros((nn, ny + n))
-        A_[:, 0:ny] = A1_.reshape(-1, ny)
+        P1 = np.vstack(P)
+        rr = P[0].shape[0]
+        A_[:, :ny] = P1
 
         for i in range(1, r):
-            Pi = P[:, i:r]
-            Oi = Or[0 : r - i, :]
-            Ni = Pi @ Oi
-            j = (i - 1) * (r - n)
-            A_[j : j + Ni.shape[0], ny : ny + Ni.shape[1]] = Ni
+            Pi_tilda = np.hstack(P[i:])
+            Ori = Or[: rny - (ny * i)]
+            N = Pi_tilda @ Ori
+            A_[(i - 1) * rr : i * rr, ny:] = N
 
-        # print(A_)
-
-        M = (U2.T @ L31 @ np.linalg.inv(L11)).ravel(order="F")
+        M = U2.T @ L31 @ np.linalg.inv(L11)
+        Mi = np.split(M, M.shape[1] // nu, axis=1)
+        M = np.vstack(Mi)
+        print(f"M shape: {M.shape}")
 
         x_, *_ = scipy.linalg.lstsq(A_, M)
-        # x_ = scipy.linalg.pinv(A_) @ M
 
-        D = x_[0:ny]
-        B = x_[ny : ny + n]
+        D = x_[:ny, :]
+        B = x_[ny:, :]
 
         return A, B, C, D
 
@@ -121,45 +126,34 @@ def ident(self, order):
         Uf = U[r : 2 * r, :]
 
         Wp = np.vstack((Up, Yp))
-        # Wf = np.vstack((Uf,Yf))
         Psi = np.vstack((np.vstack((Uf, Wp)), Yf))
 
         L, Q = self.lq(Psi)
-
-        L11 = L[0:r, 0:r]
-        # L12 = L[  0:r  ,r:2*r]
-        # L21 = L[  r:3*r,0:r]
         L22 = L[r : 3 * r, r : 3 * r]
-        L31 = L[3 * r : 4 * r, 0:r]
         L32 = L[3 * r : 4 * r, r : 3 * r]
 
-        # print(L32.shape, np.linalg.pinv(L22).shape, Wp.shape)
         Gamma_r = L32 @ np.linalg.pinv(L22) @ Wp  # oblique projection
-        # print(Gamma_r.shape)
+
         U_, s_, V_ = scipy.linalg.svd(Gamma_r, full_matrices=False)
-        # print(U_.shape, s_.shape, V_.shape)
 
         self.singular_values = s_
+        Sigma_sqrt = np.diag(np.sqrt(s_[:n]))
 
+        # ====================================================================
+        # V1 = V_[0:n, :]
+        # Xf = (
+        #     Sigma_sqrt @ V1
+        # )  # state matrix # TANGIRALA SAYS IT SHOULD BE TRANSPOSED !?!
+
+        # A, B, C, D = self._abcd_state(Xf, s, n, r)
+        # ====================================================================
         U1 = U_[:, 0:n]
         U2 = U_[:, n:r]
-        Sigma_sqrt = np.diag(np.sqrt(s_[:n]))
-        # Sigma_sqrt = np.zeros(Gamma_r.shape, s_.dtype)
-        # np.fill_diagonal(Sigma_sqrt, np.sqrt(s_))
-        # Sigma_sqrt = Sigma_sqrt[:n, :]
-        # Sigma_sqrt = scipy.linalg.diagsvd(s_, *Gamma_r.shape)
-        V1 = V_[0:n, :]
-        # V1 = V_[:, :n]
-        # V2 = V_[:,n:r]
-
-        Or = U1 @ Sigma_sqrt  # extended observability matrix
-        Xf = (
-            Sigma_sqrt @ V1
-        )  # state matrix # TANGIRALA SAYS IT SHOULD BE TRANSPOSED !?!
+        L11 = L[0:r, 0:r]
+        L31 = L[3 * r : 4 * r, 0:r]
 
-        # print(s,n,r)
-        A, B, C, D = self._abcd_state(Xf, s, n, r)
-        # A, B, C, D = self._abcd_observability_matrix(self, U1, U2, L11, L31, Sigma_sqrt, n, r):
+        A, B, C, D = self._abcd_observability_matrix(U1, U2, L11, L31, Sigma_sqrt, n, r)
+        # ====================================================================
 
         mod = StateSpaceModel(A=A, B=B, C=C, D=D, Ts=self.Ts)
         mod.info["Hankel singular values"] = s_
@@ -171,91 +165,6 @@ def name():
         return "n4sid"
 
 
-# class N4SID2(SubspaceIdent):
-#     def __init__(self, data, y_name, u_name, settings={}):
-#         super().__init__(data, y_name, u_name, settings=settings)
-
-#         # estimate extended observability matrix and states. Then estimate A, B, C, and D in one go.
-#         # (Tangirala 2014)
-
-#     def ident(self, order):
-#         if isinstance(order, (tuple, list)):
-#             n = order[0]
-#         else:
-#             n = order
-
-#         r = 2 * n + 1  # window length
-
-#         Y = self.hankel(self.y, 2 * r)
-#         U = self.hankel(self.u, 2 * r)
-
-#         s = Y.shape[1]
-
-#         Yp = Y[0:r, :]
-#         YpPlus = Y[0:r+1, :]
-#         Up = U[0:r, :]
-#         UpPlus = U[0:r+1, :]
-
-#         Yf = Y[r : 2 * r, :]
-#         YfMinus = Y[r+1 : 2 * r, :]
-#         Uf = U[r : 2 * r, :]
-#         UfMinus = U[r+1 : 2 * r, :]
-
-#         Wp = np.vstack((Up, Yp))
-#         WpPlus = np.vstack((UpPlus, YpPlus))
-#         # Wf = np.vstack((Uf,Yf))
-
-
-#         # Projection 1
-#         Psi = np.vstack((np.vstack((Uf, Wp)), Yf))
-
-#         L, _ = self.lq(Psi)
-
-#         L22 = L[r : 3 * r, r : 3 * r]
-#         L32 = L[3 * r : 4 * r, r : 3 * r]
-
-
-#         # Projection 2
-#         PsiAst = np.vstack((np.vstack((Uf, Wp)), Yf))
-
-#         L, _ = self.lq(PsiAst)
-
-#         LAst22 = L[r-1 : 3 * r+1, r-1 : 3 * r+1]
-#         LAst32 = L[3 * r+1 : 4 * r, r-1 : 3 * r+1]
-
-
-#         # print(L32.shape, np.linalg.pinv(L22).shape, Wp.shape)
-#         Gamma_r = L32 @ np.linalg.pinv(L22) @ Wp  # oblique projection
-#         GammaAst_r = LAst32 @ np.linalg.pinv(LAst22) @ WpPlus  # oblique projection
-#         U_, s_, _ = scipy.linalg.svd(Gamma_r, full_matrices=False)
-
-#         self.singular_values = s_
-
-#         U1 = U_[:,0:n]
-#         Sigma_sqrt = np.diag(np.sqrt(s_[:n]))
-
-#         Or = U1 @ Sigma_sqrt
-#         Xf = np.linalg.pinv(Or) @ Gamma_r
-#         XfPlus =  np.linalg.pinv(Or[:-1,:]) @ GammaAst_r
-
-#         Y_ = np.vstack((XfPlus[:, : s - 2], self.y[r : r + s - 2].T))
-#         X_ = np.vstack((Xf[:, : s - 2], self.u[r : r + s - 2].T))
-#         Theta = Y_ @ np.linalg.pinv(X_)
-#         A = Theta[:n, :n]
-#         B = Theta[:n, n].ravel()
-#         C = Theta[n, :n].ravel()
-#         D = Theta[n, n]
-
-#         mod = StateSpaceModel(A=A, B=B, C=C, D=D, Ts=self.Ts)
-#         mod.info["Hankel singular values"] = s_
-
-#         return mod
-
-#     @staticmethod
-#     def name():
-#         return "n4sid2"
-
-
 class PO_MOESP(SubspaceIdent):
     def __init__(self, data, y_name, u_name, settings={}):
         super().__init__(data, y_name, u_name, settings=settings)
@@ -271,15 +180,11 @@ def ident(self, order):
         else:
             n = order
 
-        N = self.y.shape[0]
-
-        r = 2 * n + 1  # window length
+        r = (2 * n + 1) * self.nu * self.ny  # window length
 
         Y = self.hankel(self.y, 2 * r)
         U = self.hankel(self.u, 2 * r)
 
-        s = Y.shape[1]
-
         Yp = Y[0:r, :]
         Up = U[0:r, :]
 
@@ -287,22 +192,16 @@ def ident(self, order):
         Uf = U[r : 2 * r, :]
 
         Wp = np.vstack((Up, Yp))
-        # Wf = np.vstack((Uf,Yf))
-        # Psi = 1.0 / N * np.vstack((np.vstack((Uf, Wp)), Yf))
         Psi = np.vstack((np.vstack((Uf, Wp)), Yf))
 
         L, Q = self.lq(Psi)
 
         L11 = L[0:r, 0:r]
-        # L12 = L[  0:r  ,r:3*r]
-        # L21 = L[  r:3*r,0:r]
-        # L22 = L[  r:3*r,r:3*r]
         L31 = L[3 * r : 4 * r, 0:r]
         L32 = L[3 * r : 4 * r, r : 3 * r]
 
-        # Gamma_r = 1./np.sqrt(N) * L32
         Gamma_r = L32
-        # print(Gamma_r.shape)
+
         U_, s_, V_ = scipy.linalg.svd(Gamma_r, full_matrices=False)
 
         self.singular_values = s_
@@ -312,48 +211,8 @@ def ident(self, order):
         U1 = U_[:, 0:n]
         U2 = U_[:, n:r]
         Sigma_sqrt = np.diag(np.sqrt(s_[:n]))
-        V1 = V_[:n, :]
-        # V1 = V_[:,0:n]
-        # V2 = V_[:,n:r]
-
-        Or = U1 @ Sigma_sqrt  # extended observability matrix
-        Xf = np.linalg.pinv(Or) @ Gamma_r
-        # print(Xf.shape)
-
-        # C = Or[0, :]# TODO: might be wrong!!!
-        # A = scipy.linalg.pinv(Or[0:-1, :]) @ Or[1:, :]
-        # # A = scipy.linalg.lstsq(Or[0:-1,:],Or[1:,:])
-        # # print(A)
-
-        # P = U2.T
-        # # print(P)
-        # A1_ = P.ravel(order="F")
 
-        # nn = A1_.shape[0]
-        # ny = 1
-        # A_ = np.zeros((nn, ny + n))
-        # A_[:, 0:ny] = A1_.reshape(-1, ny)
-
-        # for i in range(1, r):
-        #     Pi = P[:, i:r]
-        #     Oi = Or[0 : r - i, :]
-        #     Ni = Pi @ Oi
-        #     j = (i - 1) * (r - n)
-        #     A_[j : j + Ni.shape[0], ny : ny + Ni.shape[1]] = Ni
-
-        # # print(A_)
-
-        # M = (U2.T @ L31 @ np.linalg.inv(L11)).ravel(order="F")
-
-        # x_, *_ = scipy.linalg.lstsq(A_, M)
-        # # x_ = scipy.linalg.pinv(A_) @ M
-
-        # D = x_[0:ny]
-        # B = x_[ny : ny + n]
-
-        # print(s,n,r)
         A, B, C, D = self._abcd_observability_matrix(U1, U2, L11, L31, Sigma_sqrt, n, r)
-        # A, B, C, D = self._abcd_state(Xf, s, n, r)
 
         mod = StateSpaceModel(A=A, B=B, C=C, D=D, Ts=self.Ts)
         mod.info["Hankel singular values"] = s_
diff --git a/src/llsi/utils.py b/src/llsi/utils.py
index ab2b8c1..4f12546 100644
--- a/src/llsi/utils.py
+++ b/src/llsi/utils.py
@@ -1,7 +1,8 @@
-from .sysidalg import sysid
 import numpy as np
 import scipy.optimize
 
+from .sysidalg import sysid
+
 
 def cv(
     training_data,
@@ -13,7 +14,6 @@ def cv(
     settings={},
     bounds=(0, 100),
 ):
-
     def fun(l):
         settings["lambda"] = 10**l
         mod = sysid(


Test patch:
diff --git a/tests/conftest.py b/tests/conftest.py
index c87ba06..5b7bae2 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -73,8 +73,8 @@ def data_siso_deterministic_stochastic(ss_mod):
 
 @pytest.fixture
 def data_mimo_deterministic():
-    filt0 = StateSpaceModel.from_PT1(1.0, 0.8, Ts=1.0)
-    filt1 = StateSpaceModel.from_PT1(2.0, 0.6, Ts=1.0)
+    filt0 = StateSpaceModel.from_PT1(1.0, 7.0, Ts=1.0)
+    filt1 = StateSpaceModel.from_PT1(2.0, 5.0, Ts=1.0)
     data0 = generate_data(filt0, noise=1e-9, seed=42)
     data0.center()
     data1 = generate_data(filt1, noise=1e-9, seed=43)
diff --git a/tests/test_subspace.py b/tests/test_subspace.py
index ed3df77..ffe642b 100644
--- a/tests/test_subspace.py
+++ b/tests/test_subspace.py
@@ -47,34 +47,35 @@ def test_po_moesp(data_siso_deterministic, ss_mod):
     )
 
 
-def test_n4sid_mimo(data_mimo_deterministic, ss_mod):
-    # data_mimo_deterministic.equidistant()
+def test_n4sid_mimo(data_mimo_deterministic):
     identifyer = N4SID(data_mimo_deterministic, ["y0", "y1"], ["u0", "u1"])
+    sys0 = N4SID(data_mimo_deterministic, "y0", "u0").ident(1)
+    sys1 = N4SID(data_mimo_deterministic, "y1", "u1").ident(1)
     mod = identifyer.ident(2)
     print(mod.info["Hankel singular values"])
     np.testing.assert_allclose(
         mod.info["Hankel singular values"],
         [
-            6.40829670e00,
-            3.94039974e-04,
-            2.47748028e-08,
-            1.01130860e-08,
-            9.63664129e-09,
-            7.82223423e-09,
-            7.34439653e-09,
-            7.30785370e-09,
-            6.41191499e-09,
-            5.96869347e-09,
-            5.92269632e-09,
-            5.38759994e-09,
-            5.29767778e-09,
-            5.02836434e-09,
-            4.37323782e-09,
-            4.03590185e-09,
-            3.79029097e-09,
-            3.32988710e-09,
-            2.74121403e-09,
-            2.40889178e-09,
+            5.23238711e01,
+            2.46007773e01,
+            9.15587701e-02,
+            9.90325106e-09,
+            9.65156298e-09,
+            8.85547497e-09,
+            8.51780379e-09,
+            8.32136576e-09,
+            7.44491698e-09,
+            7.03837731e-09,
+            5.99893485e-09,
+            5.85910395e-09,
+            5.52005230e-09,
+            5.26083503e-09,
+            5.10508405e-09,
+            4.67083640e-09,
+            3.59163096e-09,
+            3.53484111e-09,
+            3.01876150e-09,
+            2.64576961e-09,
         ],
         rtol=1e-6,
         atol=1e-6,
@@ -85,14 +86,68 @@ def test_n4sid_mimo(data_mimo_deterministic, ss_mod):
     assert mod.C.shape == (2, 2)
     assert mod.D.shape == (2, 2)
 
-    # np.testing.assert_allclose(
-    #    mod.to_controllable_form().A, ss_mod.A, rtol=1e-1, atol=1e-1
-    # )
+    print(mod)
+
+    t, i = mod.step_response(10)
+    assert i.shape == (10, 2)
+    np.testing.assert_allclose(
+        i,
+        np.array(
+            [sys0.step_response(10)[1].ravel(), sys1.step_response(10)[1].ravel()]
+        ).T,
+        rtol=1e-5,
+        atol=1e-5,
+    )
 
-    t, i = mod.impulse_response()
-    assert i.shape == (100, 2)
-    np.testing.assert_allclose(i[:10, :], [[0], [0]])
 
-    d = data_mimo_deterministic
-    y = mod.simulate([d["u0"], d["u1"]])
-    np.testing.assert_allclose(y[:10, :], [[0], [0]])
+def test_po_moesp_mimo(data_mimo_deterministic):
+    identifyer = PO_MOESP(data_mimo_deterministic, ["y0", "y1"], ["u0", "u1"])
+    sys0 = PO_MOESP(data_mimo_deterministic, "y0", "u0").ident(1)
+    sys1 = PO_MOESP(data_mimo_deterministic, "y1", "u1").ident(1)
+    mod = identifyer.ident(2)
+    print(mod.info["Hankel singular values"])
+    np.testing.assert_allclose(
+        mod.info["Hankel singular values"],
+        [
+            5.13375445e01,
+            2.45203750e01,
+            9.15569690e-02,
+            9.88027549e-09,
+            9.63367326e-09,
+            8.84678201e-09,
+            8.50474100e-09,
+            8.31302996e-09,
+            7.42880630e-09,
+            7.02472812e-09,
+            5.99054895e-09,
+            5.85392065e-09,
+            5.51099680e-09,
+            5.25372488e-09,
+            5.10043257e-09,
+            4.66624062e-09,
+            3.58670073e-09,
+            3.52980561e-09,
+            3.01474135e-09,
+            2.64218261e-09,
+        ],
+        rtol=1e-6,
+        atol=1e-6,
+    )
+
+    assert mod.A.shape == (2, 2)
+    assert mod.B.shape == (2, 2)
+    assert mod.C.shape == (2, 2)
+    assert mod.D.shape == (2, 2)
+
+    print(mod)
+
+    t, i = mod.step_response(10)
+    assert i.shape == (10, 2)
+    np.testing.assert_allclose(
+        i,
+        np.array(
+            [sys0.step_response(10)[1].ravel(), sys1.step_response(10)[1].ravel()]
+        ).T,
+        rtol=1e-5,
+        atol=1e-5,
+    )