In previous lessons, a minimal realization was defined as a realization that is both reachable and observable. We now ask a subtler question: when do two different state-space descriptions represent the same system? For continuous-time LTI systems,

\[ \Sigma: \quad \dot{\mathbf{x} }(t)=\mathbf{A}\mathbf{x}(t)+\mathbf{B}\mathbf{u}(t), \qquad \mathbf{y}(t)=\mathbf{C}\mathbf{x}(t)+\mathbf{D}\mathbf{u}(t), \]

the internal state dimension is \( n \), the input dimension is \( m \), and the output dimension is \( p \). The corresponding transfer matrix is

\[ \mathbf{G}(s)=\mathbf{C}(s\mathbf{I}-\mathbf{A})^{-1}\mathbf{B}+\mathbf{D}. \]

The state variables are not directly unique. A coordinate change can alter \( \mathbf{A},\mathbf{B},\mathbf{C} \) without altering physical behavior. But an even weaker phenomenon is possible: two realizations may have exactly the same transfer matrix while one contains hidden unreachable or unobservable dynamics. This is the difference between internal and external equivalence.

flowchart TD
  A["State-space realization Sigma = A,B,C,D"] --> B["Internal equivalence"]
  A --> C["External equivalence"]
  B --> D["Same state dimension and similarity transform"]
  C --> E["Same transfer matrix and impulse response"]
  D --> F["Always implies external equivalence"]
  E --> G["May hide unreachable or unobservable modes"]
  F --> H["If both realizations are minimal, external equivalence implies internal equivalence"]
  G --> H
        

Two realizations \( \Sigma_1=(\mathbf{A}_1,\mathbf{B}_1,\mathbf{C}_1,\mathbf{D}_1) \) and \( \Sigma_2=(\mathbf{A}_2,\mathbf{B}_2,\mathbf{C}_2,\mathbf{D}_2) \) of the same input-output dimensions are internally equivalent if there exists a nonsingular matrix \( \mathbf{T}\in\mathbb{R}^{n\times n} \) such that, for the coordinate relation \( \mathbf{x}_2=\mathbf{T}\mathbf{x}_1 \),

\[ \mathbf{A}_2=\mathbf{T}\mathbf{A}_1\mathbf{T}^{-1},\qquad \mathbf{B}_2=\mathbf{T}\mathbf{B}_1, \qquad \mathbf{C}_2=\mathbf{C}_1\mathbf{T}^{-1}, \qquad \mathbf{D}_2=\mathbf{D}_1. \]

This is precisely the similarity transformation introduced earlier in the course. It is an equivalence relation because it is reflexive (\( \mathbf{T}=\mathbf{I} \)), symmetric (\( \mathbf{T}^{-1} \) reverses the transformation), and transitive (products of nonsingular matrices are nonsingular).

Theorem 1: Internal equivalence preserves the transfer matrix.

Proof: Substitute the similarity relations into the transfer matrix of \( \Sigma_2 \):

\[ \begin{aligned} \mathbf{G}_2(s) &=\mathbf{C}_2(s\mathbf{I}-\mathbf{A}_2)^{-1}\mathbf{B}_2+\mathbf{D}_2 \\ &=\mathbf{C}_1\mathbf{T}^{-1} \left(s\mathbf{I}-\mathbf{T}\mathbf{A}_1\mathbf{T}^{-1}\right)^{-1} \mathbf{T}\mathbf{B}_1+\mathbf{D}_1 \\ &=\mathbf{C}_1\mathbf{T}^{-1} \left(\mathbf{T}(s\mathbf{I}-\mathbf{A}_1)\mathbf{T}^{-1}\right)^{-1} \mathbf{T}\mathbf{B}_1+\mathbf{D}_1 \\ &=\mathbf{C}_1\mathbf{T}^{-1} \mathbf{T}(s\mathbf{I}-\mathbf{A}_1)^{-1}\mathbf{T}^{-1} \mathbf{T}\mathbf{B}_1+\mathbf{D}_1 \\ &=\mathbf{C}_1(s\mathbf{I}-\mathbf{A}_1)^{-1}\mathbf{B}_1+\mathbf{D}_1 =\mathbf{G}_1(s). \quad \square \end{aligned} \]

Therefore, internal equivalence is a stronger relation than equality of transfer matrices. It preserves all coordinate-free internal invariants: eigenvalues of \( \mathbf{A} \), Jordan block sizes, reachability, observability, controllability indices, and the minimal state dimension.

Two realizations are externally equivalent if they have the same input-output behavior for every admissible input under zero initial condition. For finite-dimensional proper rational LTI systems, this is equivalent to equality of transfer matrices:

\[ \Sigma_1 \sim_{ext} \Sigma_2 \quad \Longleftrightarrow \quad \mathbf{G}_1(s)=\mathbf{G}_2(s). \]

In the time domain, the impulse response is

\[ \mathbf{g}(t)=\mathbf{C}e^{\mathbf{A}t}\mathbf{B}+\mathbf{D}\delta(t), \qquad t\ge 0. \]

Expanding the resolvent at large \( s \) gives the Markov-parameter characterization:

\[ \mathbf{G}(s)=\mathbf{D}+\sum_{k=0}^{\infty} \mathbf{C}\mathbf{A}^k\mathbf{B}\,s^{-(k+1)}. \]

Hence, external equivalence implies

\[ \mathbf{D}_1=\mathbf{D}_2, \qquad \mathbf{C}_1\mathbf{A}_1^k\mathbf{B}_1 =\mathbf{C}_2\mathbf{A}_2^k\mathbf{B}_2, \qquad k=0,1,2,\ldots. \]

The important limitation is that the transfer matrix only sees modes that are both reachable from the input and observable at the output. A hidden mode may exist in the state equation and still be absent from \( \mathbf{G}(s) \). For example,

\[ \mathbf{A}_3=\begin{bmatrix} -1&0&0\\0&-2&0\\0&0&5 \end{bmatrix},\quad \mathbf{B}_3=\begin{bmatrix}1\\1\\0\end{bmatrix},\quad \mathbf{C}_3=\begin{bmatrix}1&1&0\end{bmatrix},\quad \mathbf{D}_3=0. \]

The third state has eigenvalue \( 5 \), but it is not reached by \( u \) and is not measured by \( y \). Therefore it is externally invisible. The transfer matrix is the same as the two-state realization

\[ \mathbf{A}_1=\begin{bmatrix} -1&0\\0&-2 \end{bmatrix},\quad \mathbf{B}_1=\begin{bmatrix}1\\1\end{bmatrix},\quad \mathbf{C}_1=\begin{bmatrix}1&1\end{bmatrix},\quad \mathbf{D}_1=0, \]

\[ G_1(s)=G_3(s)=\frac{1}{s+1}+\frac{1}{s+2}. \]

Nevertheless, these two realizations are not internally equivalent, because internal equivalence requires equal state dimension and a nonsingular similarity map.

Minimal realization theory says that external behavior determines only the reachable-and-observable part of a system. If no hidden part is present, equality of input-output behavior forces equality of internal dynamics up to a coordinate change.

Theorem 2: If \( \Sigma_1 \) and \( \Sigma_2 \) are minimal realizations of the same transfer matrix, then they are internally equivalent. That is, there exists a nonsingular matrix \( \mathbf{T} \) such that

\[ \mathbf{A}_2\mathbf{T}=\mathbf{T}\mathbf{A}_1, \qquad \mathbf{B}_2=\mathbf{T}\mathbf{B}_1, \qquad \mathbf{C}_2\mathbf{T}=\mathbf{C}_1, \qquad \mathbf{D}_2=\mathbf{D}_1. \]

The displayed form is equivalent to the previous similarity convention \( \mathbf{A}_2=\mathbf{T}\mathbf{A}_1\mathbf{T}^{-1} \).

Proof sketch: Since \( \Sigma_1 \) is reachable, every state can be written as a finite linear combination of columns of the form \( \mathbf{A}_1^k\mathbf{B}_1\mathbf{v} \). Define a linear map on reachable generators by

\[ \mathbf{T}\left(\sum_k \mathbf{A}_1^k\mathbf{B}_1\mathbf{v}_k\right) =\sum_k \mathbf{A}_2^k\mathbf{B}_2\mathbf{v}_k. \]

Equality of all Markov parameters gives \( \mathbf{C}_1\mathbf{A}_1^j\mathbf{A}_1^k\mathbf{B}_1 =\mathbf{C}_2\mathbf{A}_2^j\mathbf{A}_2^k\mathbf{B}_2 \) for all nonnegative \( j,k \). Therefore, if a linear combination of reachable generators is zero in the first realization, its image produces zero output derivatives in the second realization. Observability of \( \Sigma_2 \) forces that image to be zero, so \( \mathbf{T} \) is well-defined. Symmetry of the argument gives nonsingularity. The same generator definition gives \( \mathbf{B}_2=\mathbf{T}\mathbf{B}_1 \) and \( \mathbf{A}_2\mathbf{T}=\mathbf{T}\mathbf{A}_1 \). Equality of the output maps on reachable states gives \( \mathbf{C}_2\mathbf{T}=\mathbf{C}_1 \). Thus the two minimal realizations differ only by a similarity transformation. \( \square \)

flowchart TD
  S["Compare two realizations"] --> Q1["Same transfer matrix?"]
  Q1 -->|"no"| N1["Not externally equivalent"]
  Q1 -->|"yes"| Q2["Both reachable and observable?"]
  Q2 -->|"no"| E1["Externally equivalent; \ninternal equivalence not guaranteed"]
  Q2 -->|"yes"| E2["Minimal realizations of \nsame transfer matrix"]
  E2 --> I1["There exists a nonsingular state map T"]
  I1 --> I2["Internally equivalent"]
        

The distinction has several important consequences for analysis and design.

1. Transfer functions do not reveal all internal modes. The poles of \( \mathbf{G}(s) \) are only the modes that survive after unreachable/unobservable cancellations. The spectrum of \( \mathbf{A} \) may contain extra modes that do not appear as transfer-function poles.

2. External stability is weaker than internal stability. A transfer matrix can be stable while a nonminimal realization contains an unstable hidden mode. This is why internal stability in state-space design must be checked using the state matrix of the implemented realization, not only the transfer matrix.

3. Minimal realizations are unique up to similarity. Thus the minimal state dimension is a property of the transfer matrix, not of a particular coordinate choice. This dimension is also called the McMillan degree.

4. Numerical realization algorithms should remove hidden modes. If the controllability matrix \( \mathcal{C}=[\mathbf{B}\;\mathbf{A}\mathbf{B}\;\cdots\;\mathbf{A}^{n-1}\mathbf{B}] \) or the observability matrix \( \mathcal{O}=[\mathbf{C}^\top\;(\mathbf{C}\mathbf{A})^\top\;\cdots\;(\mathbf{C}\mathbf{A}^{n-1})^\top]^\top \) is rank deficient, then the realization is not minimal:

\[ \operatorname{rank}(\mathcal{C})=n, \qquad \operatorname{rank}(\mathcal{O})=n \quad \Longleftrightarrow \quad \Sigma \text{ is minimal}. \]

The following programs compare three systems: a minimal two-state realization, an internally equivalent realization obtained by a similarity transformation, and a nonminimal externally equivalent realization containing a hidden unstable mode. The same conceptual checks are implemented in Python, C++, Java, MATLAB/Simulink, and Wolfram Mathematica.

6.1 Python Implementation

File: Chapter20_Lesson3.py

# Chapter20_Lesson3.py
# Internal vs external equivalence of continuous-time LTI systems.
# Required library: numpy

import numpy as np


def ctrb(A, B):
    """Kalman controllability matrix [B AB ... A^(n-1)B]."""
    n = A.shape[0]
    blocks = []
    Ak = np.eye(n)
    for _ in range(n):
        blocks.append(Ak @ B)
        Ak = A @ Ak
    return np.hstack(blocks)


def obsv(A, C):
    """Kalman observability matrix [C; CA; ...; CA^(n-1)]."""
    n = A.shape[0]
    blocks = []
    Ak = np.eye(n)
    for _ in range(n):
        blocks.append(C @ Ak)
        Ak = Ak @ A
    return np.vstack(blocks)


def rank(M, tol=1e-10):
    return np.linalg.matrix_rank(M, tol=tol)


def transfer_value(A, B, C, D, s):
    """Evaluate G(s) = C(sI-A)^(-1)B + D."""
    n = A.shape[0]
    return C @ np.linalg.solve(s * np.eye(n) - A, B) + D


def markov_parameters(A, B, C, D, count=6):
    """Return D and the first count continuous-time Markov matrices C A^k B."""
    n = A.shape[0]
    params = [D.copy()]
    Ak = np.eye(n)
    for _ in range(count):
        params.append(C @ Ak @ B)
        Ak = A @ Ak
    return params


def is_internally_equivalent(A1, B1, C1, D1, A2, B2, C2, D2, T, tol=1e-9):
    """Check A2=T A1 T^(-1), B2=T B1, C2=C1 T^(-1), D2=D1."""
    Ti = np.linalg.inv(T)
    tests = {
        "A2 = T A1 T^(-1)": np.linalg.norm(A2 - T @ A1 @ Ti),
        "B2 = T B1": np.linalg.norm(B2 - T @ B1),
        "C2 = C1 T^(-1)": np.linalg.norm(C2 - C1 @ Ti),
        "D2 = D1": np.linalg.norm(D2 - D1),
    }
    return tests, all(v < tol for v in tests.values())


def compare_transfers(systems, sample_points):
    """Compare scalar transfer functions at real sample points."""
    base = systems[0]
    rows = []
    for s in sample_points:
        g0 = transfer_value(*base, s=s)
        row = [float(g0[0, 0])]
        for sys in systems[1:]:
            gi = transfer_value(*sys, s=s)
            row.append(float(gi[0, 0]))
        rows.append((s, row))
    return rows


def print_system_report(name, A, B, C, D):
    n = A.shape[0]
    print(f"\n{name}")
    print("A =\n", A)
    print("B =\n", B)
    print("C =\n", C)
    print("D =\n", D)
    print("rank controllability =", rank(ctrb(A, B)), "of", n)
    print("rank observability   =", rank(obsv(A, C)), "of", n)
    print("minimal?", rank(ctrb(A, B)) == n and rank(obsv(A, C)) == n)
    print("eigenvalues(A) =", np.linalg.eigvals(A))


if __name__ == "__main__":
    # A minimal two-state SISO realization.
    A1 = np.array([[-1.0, 0.0],
                   [ 0.0,-2.0]])
    B1 = np.array([[1.0],
                   [1.0]])
    C1 = np.array([[1.0, 1.0]])
    D1 = np.array([[0.0]])

    # Internal equivalence: choose a nonsingular state-coordinate map x2 = T x1.
    T = np.array([[1.0, 2.0],
                  [0.5, 1.5]])
    A2 = T @ A1 @ np.linalg.inv(T)
    B2 = T @ B1
    C2 = C1 @ np.linalg.inv(T)
    D2 = D1.copy()

    # External equivalence only: append a hidden mode. This mode is neither
    # reachable from u nor visible in y, so it does not change G(s).
    A3 = np.array([[-1.0, 0.0, 0.0],
                   [ 0.0,-2.0, 0.0],
                   [ 0.0, 0.0, 5.0]])
    B3 = np.array([[1.0],
                   [1.0],
                   [0.0]])
    C3 = np.array([[1.0, 1.0, 0.0]])
    D3 = D1.copy()

    # Another externally equivalent nonminimal realization with a different hidden mode.
    A4 = np.array([[-1.0, 0.0,  0.0],
                   [ 0.0,-2.0,  0.0],
                   [ 0.0, 0.0,-10.0]])
    B4 = B3.copy()
    C4 = C3.copy()
    D4 = D1.copy()

    for item in [("Sigma_1 minimal", A1, B1, C1, D1),
                 ("Sigma_2 internally equivalent to Sigma_1", A2, B2, C2, D2),
                 ("Sigma_3 externally equivalent but nonminimal", A3, B3, C3, D3),
                 ("Sigma_4 externally equivalent but different hidden eigenvalue", A4, B4, C4, D4)]:
        print_system_report(*item)

    tests, ok = is_internally_equivalent(A1, B1, C1, D1, A2, B2, C2, D2, T)
    print("\nInternal-equivalence residuals for Sigma_1 and Sigma_2:")
    for k, v in tests.items():
        print(f"  {k:22s}: {v:.3e}")
    print("Internally equivalent?", ok)

    print("\nTransfer-function samples G_i(s):")
    systems = [(A1, B1, C1, D1), (A2, B2, C2, D2), (A3, B3, C3, D3), (A4, B4, C4, D4)]
    for s, values in compare_transfers(systems, [0.1, 1.0, 3.0, 10.0]):
        print(f"s={s:4.1f} ->", values)

    print("\nFirst Markov parameters D, C B, C A B, ... for Sigma_1 and Sigma_3:")
    for k, (m1, m3) in enumerate(zip(markov_parameters(A1, B1, C1, D1),
                                     markov_parameters(A3, B3, C3, D3))):
        print(f"k={k}: Sigma_1 {m1.ravel()}   Sigma_3 {m3.ravel()}")

6.2 C++ Implementation

File: Chapter20_Lesson3.cpp

// Chapter20_Lesson3.cpp
// Internal vs external equivalence of continuous-time LTI systems.
// Dependency: Eigen 3 (header-only). Example compile command:
// g++ -std=c++17 Chapter20_Lesson3.cpp -I /path/to/eigen -O2 -o Chapter20_Lesson3

#include <Eigen/Dense>
#include <iostream>
#include <iomanip>
#include <vector>

using Eigen::MatrixXd;

MatrixXd ctrb(const MatrixXd& A, const MatrixXd& B) {
    int n = A.rows();
    MatrixXd Ctrb(n, n * B.cols());
    MatrixXd Ak = MatrixXd::Identity(n, n);
    for (int k = 0; k < n; ++k) {
        Ctrb.block(0, k * B.cols(), n, B.cols()) = Ak * B;
        Ak = A * Ak;
    }
    return Ctrb;
}

MatrixXd obsv(const MatrixXd& A, const MatrixXd& C) {
    int n = A.rows();
    MatrixXd Obsv(n * C.rows(), n);
    MatrixXd Ak = MatrixXd::Identity(n, n);
    for (int k = 0; k < n; ++k) {
        Obsv.block(k * C.rows(), 0, C.rows(), n) = C * Ak;
        Ak = Ak * A;
    }
    return Obsv;
}

int rankOf(const MatrixXd& M, double tol = 1e-10) {
    Eigen::FullPivLU<MatrixXd> lu(M);
    lu.setThreshold(tol);
    return lu.rank();
}

MatrixXd transferValue(const MatrixXd& A, const MatrixXd& B,
                       const MatrixXd& C, const MatrixXd& D, double s) {
    int n = A.rows();
    MatrixXd resolvent = (s * MatrixXd::Identity(n, n) - A).inverse();
    return C * resolvent * B + D;
}

std::vector<MatrixXd> markovParameters(const MatrixXd& A, const MatrixXd& B,
                                       const MatrixXd& C, const MatrixXd& D,
                                       int count = 6) {
    int n = A.rows();
    std::vector<MatrixXd> params;
    params.push_back(D);
    MatrixXd Ak = MatrixXd::Identity(n, n);
    for (int k = 0; k < count; ++k) {
        params.push_back(C * Ak * B);
        Ak = A * Ak;
    }
    return params;
}

void report(const std::string& name, const MatrixXd& A, const MatrixXd& B,
            const MatrixXd& C, const MatrixXd& D) {
    int n = A.rows();
    std::cout << "\n" << name << "\n";
    std::cout << "A =\n" << A << "\nB =\n" << B << "\nC =\n" << C << "\nD =\n" << D << "\n";
    std::cout << "rank controllability = " << rankOf(ctrb(A, B)) << " of " << n << "\n";
    std::cout << "rank observability   = " << rankOf(obsv(A, C)) << " of " << n << "\n";
    std::cout << "minimal? " << ((rankOf(ctrb(A, B)) == n && rankOf(obsv(A, C)) == n) ? "yes" : "no") << "\n";
}

int main() {
    MatrixXd A1(2, 2), B1(2, 1), C1(1, 2), D1(1, 1);
    A1 << -1.0, 0.0,
           0.0,-2.0;
    B1 << 1.0,
          1.0;
    C1 << 1.0, 1.0;
    D1 << 0.0;

    MatrixXd T(2, 2);
    T << 1.0, 2.0,
         0.5, 1.5;
    MatrixXd A2 = T * A1 * T.inverse();
    MatrixXd B2 = T * B1;
    MatrixXd C2 = C1 * T.inverse();
    MatrixXd D2 = D1;

    MatrixXd A3(3, 3), B3(3, 1), C3(1, 3), D3(1, 1);
    A3 << -1.0, 0.0, 0.0,
           0.0,-2.0, 0.0,
           0.0, 0.0, 5.0;
    B3 << 1.0, 1.0, 0.0;
    C3 << 1.0, 1.0, 0.0;
    D3 << 0.0;

    report("Sigma_1 minimal", A1, B1, C1, D1);
    report("Sigma_2 internally equivalent to Sigma_1", A2, B2, C2, D2);
    report("Sigma_3 externally equivalent but nonminimal", A3, B3, C3, D3);

    std::cout << "\nInternal-equivalence residuals for Sigma_1 and Sigma_2:\n";
    std::cout << "A residual = " << (A2 - T * A1 * T.inverse()).norm() << "\n";
    std::cout << "B residual = " << (B2 - T * B1).norm() << "\n";
    std::cout << "C residual = " << (C2 - C1 * T.inverse()).norm() << "\n";
    std::cout << "D residual = " << (D2 - D1).norm() << "\n";

    std::cout << "\nTransfer-function samples G_i(s):\n";
    for (double s : {0.1, 1.0, 3.0, 10.0}) {
        std::cout << std::fixed << std::setprecision(6)
                  << "s=" << s
                  << "  G1=" << transferValue(A1, B1, C1, D1, s)(0, 0)
                  << "  G2=" << transferValue(A2, B2, C2, D2, s)(0, 0)
                  << "  G3=" << transferValue(A3, B3, C3, D3, s)(0, 0)
                  << "\n";
    }

    std::cout << "\nFirst Markov parameters for Sigma_1 and Sigma_3:\n";
    auto m1 = markovParameters(A1, B1, C1, D1);
    auto m3 = markovParameters(A3, B3, C3, D3);
    for (size_t k = 0; k < m1.size(); ++k) {
        std::cout << "k=" << k << "  Sigma_1=" << m1[k](0, 0)
                  << "  Sigma_3=" << m3[k](0, 0) << "\n";
    }
    return 0;
}

6.3 Java Implementation

File: Chapter20_Lesson3.java

// Chapter20_Lesson3.java
// Internal vs external equivalence of continuous-time LTI systems.
// No external library is required; this file implements small matrix utilities from scratch.

import java.util.Arrays;

public class Chapter20_Lesson3 {
    static double[][] eye(int n) {
        double[][] I = new double[n][n];
        for (int i = 0; i < n; i++) I[i][i] = 1.0;
        return I;
    }

    static double[][] add(double[][] A, double[][] B) {
        double[][] C = new double[A.length][A[0].length];
        for (int i = 0; i < A.length; i++)
            for (int j = 0; j < A[0].length; j++) C[i][j] = A[i][j] + B[i][j];
        return C;
    }

    static double[][] sub(double[][] A, double[][] B) {
        double[][] C = new double[A.length][A[0].length];
        for (int i = 0; i < A.length; i++)
            for (int j = 0; j < A[0].length; j++) C[i][j] = A[i][j] - B[i][j];
        return C;
    }

    static double[][] scale(double a, double[][] A) {
        double[][] C = new double[A.length][A[0].length];
        for (int i = 0; i < A.length; i++)
            for (int j = 0; j < A[0].length; j++) C[i][j] = a * A[i][j];
        return C;
    }

    static double[][] mul(double[][] A, double[][] B) {
        int m = A.length, n = B[0].length, p = B.length;
        double[][] C = new double[m][n];
        for (int i = 0; i < m; i++)
            for (int k = 0; k < p; k++)
                for (int j = 0; j < n; j++) C[i][j] += A[i][k] * B[k][j];
        return C;
    }

    static double[][] hcat(double[][][] blocks) {
        int rows = blocks[0].length;
        int cols = 0;
        for (double[][] B : blocks) cols += B[0].length;
        double[][] out = new double[rows][cols];
        int offset = 0;
        for (double[][] B : blocks) {
            for (int i = 0; i < rows; i++)
                for (int j = 0; j < B[0].length; j++) out[i][offset + j] = B[i][j];
            offset += B[0].length;
        }
        return out;
    }

    static double[][] vcat(double[][][] blocks) {
        int cols = blocks[0][0].length;
        int rows = 0;
        for (double[][] B : blocks) rows += B.length;
        double[][] out = new double[rows][cols];
        int offset = 0;
        for (double[][] B : blocks) {
            for (int i = 0; i < B.length; i++)
                for (int j = 0; j < cols; j++) out[offset + i][j] = B[i][j];
            offset += B.length;
        }
        return out;
    }

    static double[][] inverse(double[][] A) {
        int n = A.length;
        double[][] aug = new double[n][2 * n];
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < n; j++) aug[i][j] = A[i][j];
            aug[i][n + i] = 1.0;
        }
        for (int col = 0; col < n; col++) {
            int pivot = col;
            for (int r = col + 1; r < n; r++)
                if (Math.abs(aug[r][col]) > Math.abs(aug[pivot][col])) pivot = r;
            double[] tmp = aug[col]; aug[col] = aug[pivot]; aug[pivot] = tmp;
            double div = aug[col][col];
            if (Math.abs(div) < 1e-12) throw new RuntimeException("singular matrix");
            for (int j = 0; j < 2 * n; j++) aug[col][j] /= div;
            for (int r = 0; r < n; r++) {
                if (r == col) continue;
                double factor = aug[r][col];
                for (int j = 0; j < 2 * n; j++) aug[r][j] -= factor * aug[col][j];
            }
        }
        double[][] inv = new double[n][n];
        for (int i = 0; i < n; i++)
            for (int j = 0; j < n; j++) inv[i][j] = aug[i][n + j];
        return inv;
    }

    static int rank(double[][] A) {
        double[][] M = new double[A.length][A[0].length];
        for (int i = 0; i < A.length; i++) M[i] = Arrays.copyOf(A[i], A[i].length);
        int rows = M.length, cols = M[0].length, r = 0;
        for (int c = 0; c < cols && r < rows; c++) {
            int pivot = r;
            for (int i = r + 1; i < rows; i++)
                if (Math.abs(M[i][c]) > Math.abs(M[pivot][c])) pivot = i;
            if (Math.abs(M[pivot][c]) < 1e-10) continue;
            double[] tmp = M[r]; M[r] = M[pivot]; M[pivot] = tmp;
            double div = M[r][c];
            for (int j = c; j < cols; j++) M[r][j] /= div;
            for (int i = 0; i < rows; i++) {
                if (i == r) continue;
                double factor = M[i][c];
                for (int j = c; j < cols; j++) M[i][j] -= factor * M[r][j];
            }
            r++;
        }
        return r;
    }

    static double norm(double[][] A) {
        double s = 0.0;
        for (double[] row : A) for (double v : row) s += v * v;
        return Math.sqrt(s);
    }

    static double[][] ctrb(double[][] A, double[][] B) {
        int n = A.length;
        double[][][] blocks = new double[n][][];
        double[][] Ak = eye(n);
        for (int k = 0; k < n; k++) {
            blocks[k] = mul(Ak, B);
            Ak = mul(A, Ak);
        }
        return hcat(blocks);
    }

    static double[][] obsv(double[][] A, double[][] C) {
        int n = A.length;
        double[][][] blocks = new double[n][][];
        double[][] Ak = eye(n);
        for (int k = 0; k < n; k++) {
            blocks[k] = mul(C, Ak);
            Ak = mul(Ak, A);
        }
        return vcat(blocks);
    }

    static double transfer(double[][] A, double[][] B, double[][] C, double[][] D, double s) {
        double[][] resolvent = inverse(sub(scale(s, eye(A.length)), A));
        return add(mul(mul(C, resolvent), B), D)[0][0];
    }

    static void report(String name, double[][] A, double[][] B, double[][] C, double[][] D) {
        int n = A.length;
        System.out.println("\n" + name);
        System.out.println("rank controllability = " + rank(ctrb(A, B)) + " of " + n);
        System.out.println("rank observability   = " + rank(obsv(A, C)) + " of " + n);
        System.out.println("minimal? " + (rank(ctrb(A, B)) == n && rank(obsv(A, C)) == n));
    }

    public static void main(String[] args) {
        double[][] A1 = { {-1.0, 0.0}, {0.0, -2.0} };
        double[][] B1 = { {1.0}, {1.0} };
        double[][] C1 = { {1.0, 1.0} };
        double[][] D1 = { {0.0} };

        double[][] T = { {1.0, 2.0}, {0.5, 1.5} };
        double[][] Ti = inverse(T);
        double[][] A2 = mul(mul(T, A1), Ti);
        double[][] B2 = mul(T, B1);
        double[][] C2 = mul(C1, Ti);
        double[][] D2 = D1;

        double[][] A3 = { {-1.0, 0.0, 0.0}, {0.0, -2.0, 0.0}, {0.0, 0.0, 5.0} };
        double[][] B3 = { {1.0}, {1.0}, {0.0} };
        double[][] C3 = { {1.0, 1.0, 0.0} };
        double[][] D3 = D1;

        report("Sigma_1 minimal", A1, B1, C1, D1);
        report("Sigma_2 internally equivalent to Sigma_1", A2, B2, C2, D2);
        report("Sigma_3 externally equivalent but nonminimal", A3, B3, C3, D3);

        System.out.println("\nInternal-equivalence residuals for Sigma_1 and Sigma_2:");
        System.out.println("A residual = " + norm(sub(A2, mul(mul(T, A1), Ti))));
        System.out.println("B residual = " + norm(sub(B2, mul(T, B1))));
        System.out.println("C residual = " + norm(sub(C2, mul(C1, Ti))));
        System.out.println("D residual = " + norm(sub(D2, D1)));

        System.out.println("\nTransfer-function samples G_i(s):");
        for (double s : new double[]{0.1, 1.0, 3.0, 10.0}) {
            System.out.printf("s=%4.1f  G1=% .8f  G2=% .8f  G3=% .8f%n",
                    s,
                    transfer(A1, B1, C1, D1, s),
                    transfer(A2, B2, C2, D2, s),
                    transfer(A3, B3, C3, D3, s));
        }
    }
}

6.4 MATLAB/Simulink Implementation

File: Chapter20_Lesson3.m

% Chapter20_Lesson3.m
% Internal vs external equivalence of continuous-time LTI systems.
% Uses only base MATLAB for matrices. The optional Simulink section creates
% a State-Space block if Simulink is installed.

clear; clc;

A1 = [-1 0; 0 -2];
B1 = [1; 1];
C1 = [1 1];
D1 = 0;

% Internal equivalence: x2 = T x1.
T  = [1 2; 0.5 1.5];
A2 = T*A1/T;
B2 = T*B1;
C2 = C1/T;
D2 = D1;

% External equivalence only: add an unreachable and unobservable hidden mode.
A3 = [-1 0 0; 0 -2 0; 0 0 5];
B3 = [1; 1; 0];
C3 = [1 1 0];
D3 = 0;

reportSystem('Sigma_1 minimal', A1, B1, C1, D1);
reportSystem('Sigma_2 internally equivalent to Sigma_1', A2, B2, C2, D2);
reportSystem('Sigma_3 externally equivalent but nonminimal', A3, B3, C3, D3);

fprintf('\nInternal-equivalence residuals for Sigma_1 and Sigma_2:\n');
fprintf('A residual = %.3e\n', norm(A2 - T*A1/T));
fprintf('B residual = %.3e\n', norm(B2 - T*B1));
fprintf('C residual = %.3e\n', norm(C2 - C1/T));
fprintf('D residual = %.3e\n', norm(D2 - D1));

fprintf('\nTransfer-function samples G_i(s):\n');
for s = [0.1 1.0 3.0 10.0]
    g1 = transferValue(A1, B1, C1, D1, s);
    g2 = transferValue(A2, B2, C2, D2, s);
    g3 = transferValue(A3, B3, C3, D3, s);
    fprintf('s=%4.1f  G1=% .8f  G2=% .8f  G3=% .8f\n', s, g1, g2, g3);
end

fprintf('\nFirst Markov parameters D, C B, C A B, ... for Sigma_1 and Sigma_3:\n');
M1 = markovParameters(A1, B1, C1, D1, 6);
M3 = markovParameters(A3, B3, C3, D3, 6);
for k = 1:numel(M1)
    fprintf('k=%d  Sigma_1=% .8f  Sigma_3=% .8f\n', k-1, M1{k}, M3{k});
end

% Optional Simulink generation: create a State-Space block for Sigma_1.
% This demonstrates how the same state-space matrices can be used in Simulink.
try
    modelName = 'Chapter20_Lesson3_Simulink';
    if bdIsLoaded(modelName)
        close_system(modelName, 0);
    end
    new_system(modelName);
    open_system(modelName);
    assignin('base', 'A1', A1);
    assignin('base', 'B1', B1);
    assignin('base', 'C1', C1);
    assignin('base', 'D1', D1);
    add_block('simulink/Sources/Step', [modelName '/StepInput'], 'Position', [40 80 80 110]);
    add_block('simulink/Continuous/State-Space', [modelName '/MinimalRealization'], 'Position', [140 70 280 120]);
    set_param([modelName '/MinimalRealization'], 'A', 'A1', 'B', 'B1', 'C', 'C1', 'D', 'D1');
    add_block('simulink/Sinks/Scope', [modelName '/Scope'], 'Position', [340 75 380 115]);
    add_line(modelName, 'StepInput/1', 'MinimalRealization/1');
    add_line(modelName, 'MinimalRealization/1', 'Scope/1');
    save_system(modelName);
    fprintf('\nSimulink model Chapter20_Lesson3_Simulink.slx was created.\n');
catch ME
    fprintf('\nSimulink section skipped: %s\n', ME.message);
end

function R = ctrbLocal(A, B)
    n = size(A, 1);
    R = [];
    Ak = eye(n);
    for k = 1:n
        R = [R Ak*B]; %#ok<AGROW>
        Ak = A*Ak;
    end
end

function O = obsvLocal(A, C)
    n = size(A, 1);
    O = [];
    Ak = eye(n);
    for k = 1:n
        O = [O; C*Ak]; %#ok<AGROW>
        Ak = Ak*A;
    end
end

function g = transferValue(A, B, C, D, s)
    n = size(A, 1);
    g = C*((s*eye(n) - A)\B) + D;
end

function params = markovParameters(A, B, C, D, count)
    n = size(A, 1);
    params = cell(count + 1, 1);
    params{1} = D;
    Ak = eye(n);
    for k = 1:count
        params{k + 1} = C*Ak*B;
        Ak = A*Ak;
    end
end

function reportSystem(name, A, B, C, D)
    n = size(A, 1);
    rc = rank(ctrbLocal(A, B));
    ro = rank(obsvLocal(A, C));
    fprintf('\n%s\n', name);
    disp('A ='); disp(A);
    disp('B ='); disp(B);
    disp('C ='); disp(C);
    disp('D ='); disp(D);
    fprintf('rank controllability = %d of %d\n', rc, n);
    fprintf('rank observability   = %d of %d\n', ro, n);
    fprintf('minimal? %d\n', rc == n && ro == n);
    fprintf('eigenvalues(A) = '); disp(eig(A).');
end

6.5 Wolfram Mathematica Implementation

File: Chapter20_Lesson3.nb

(* Chapter20_Lesson3.nb *)
(* Internal vs external equivalence of continuous-time LTI systems. *)

ClearAll["Global`*"];

ctrb[A_, B_] := Module[{n = Length[A]}, ArrayFlatten[{Table[MatrixPower[A, k].B, {k, 0, n - 1}]}]];
obsv[A_, C_] := Module[{n = Length[A]}, Join @@ Table[{C.MatrixPower[A, k]}, {k, 0, n - 1}]];
tf[A_, B_, C_, D_, s_] := Simplify[D + C.Inverse[s IdentityMatrix[Length[A]] - A].B];
markov[A_, B_, C_, D_, count_] := Join[{D}, Table[C.MatrixPower[A, k].B, {k, 0, count - 1}]];

A1 = { {-1, 0}, {0, -2} };
B1 = { {1}, {1} };
C1 = { {1, 1} };
D1 = { {0} };

T = { {1, 2}, {1/2, 3/2} };
A2 = Simplify[T.A1.Inverse[T]];
B2 = Simplify[T.B1];
C2 = Simplify[C1.Inverse[T]];
D2 = D1;

A3 = { {-1, 0, 0}, {0, -2, 0}, {0, 0, 5} };
B3 = { {1}, {1}, {0} };
C3 = { {1, 1, 0} };
D3 = D1;

Print["Transfer matrix of Sigma_1:"];
Print[tf[A1, B1, C1, D1, s] // MatrixForm];

Print["Transfer matrix of Sigma_2:"];
Print[tf[A2, B2, C2, D2, s] // MatrixForm];

Print["Transfer matrix of Sigma_3:"];
Print[tf[A3, B3, C3, D3, s] // MatrixForm];

Print["Sigma_1 and Sigma_2 internally equivalent residuals:"];
Print[Simplify[A2 - T.A1.Inverse[T]] // MatrixForm];
Print[Simplify[B2 - T.B1] // MatrixForm];
Print[Simplify[C2 - C1.Inverse[T]] // MatrixForm];
Print[Simplify[D2 - D1] // MatrixForm];

Print["External equivalence checks:"];
Print[FullSimplify[tf[A1, B1, C1, D1, s] == tf[A2, B2, C2, D2, s]]];
Print[FullSimplify[tf[A1, B1, C1, D1, s] == tf[A3, B3, C3, D3, s]]];

Print["Controllability and observability ranks:"];
Print[{MatrixRank[ctrb[A1, B1]], MatrixRank[obsv[A1, C1]], Length[A1]}];
Print[{MatrixRank[ctrb[A3, B3]], MatrixRank[obsv[A3, C3]], Length[A3]}];

Print["Markov parameters D, C B, C A B, ...:"];
Print[markov[A1, B1, C1, D1, 6]];
Print[markov[A3, B3, C3, D3, 6]];

Problem 1 (Similarity and Transfer Matrices): Let \( \Sigma_2 \) be obtained from \( \Sigma_1 \) by the transformation \( \mathbf{x}_2=\mathbf{T}\mathbf{x}_1 \), where \( \mathbf{T} \) is nonsingular. Prove that \( \mathbf{G}_1(s)=\mathbf{G}_2(s) \).

Solution: By definition, \( \mathbf{A}_2=\mathbf{T}\mathbf{A}_1\mathbf{T}^{-1} \), \( \mathbf{B}_2=\mathbf{T}\mathbf{B}_1 \), \( \mathbf{C}_2=\mathbf{C}_1\mathbf{T}^{-1} \), and \( \mathbf{D}_2=\mathbf{D}_1 \). Therefore,

\[ \begin{aligned} \mathbf{G}_2(s) &=\mathbf{C}_1\mathbf{T}^{-1} \left(\mathbf{T}(s\mathbf{I}-\mathbf{A}_1)\mathbf{T}^{-1}\right)^{-1} \mathbf{T}\mathbf{B}_1+\mathbf{D}_1 \\ &=\mathbf{C}_1(s\mathbf{I}-\mathbf{A}_1)^{-1}\mathbf{B}_1+\mathbf{D}_1 =\mathbf{G}_1(s). \end{aligned} \]

Problem 2 (External Equivalence without Internal Equivalence): Construct two externally equivalent systems that cannot be internally equivalent.

Solution: Consider the two-state system

\[ \mathbf{A}_1=\begin{bmatrix}-1&0\\0&-2\end{bmatrix},\quad \mathbf{B}_1=\begin{bmatrix}1\\1\end{bmatrix},\quad \mathbf{C}_1=\begin{bmatrix}1&1\end{bmatrix},\quad \mathbf{D}_1=0, \]

and the three-state system

\[ \mathbf{A}_3=\begin{bmatrix}-1&0&0\\0&-2&0\\0&0&5\end{bmatrix},\quad \mathbf{B}_3=\begin{bmatrix}1\\1\\0\end{bmatrix},\quad \mathbf{C}_3=\begin{bmatrix}1&1&0\end{bmatrix},\quad \mathbf{D}_3=0. \]

Both have transfer function \( G(s)=1/(s+1)+1/(s+2) \), but their state dimensions are different. Since internal equivalence requires a nonsingular map between state spaces of the same dimension, they are not internally equivalent.

Problem 3 (Markov Parameter Test): Suppose two strictly proper SISO realizations satisfy \( C_1A_1^kB_1=C_2A_2^kB_2 \) for \( k=0,1,2,\ldots \). Show that their transfer functions are equal.

Solution: For sufficiently large \( |s| \), the Neumann expansion gives

\[ (s\mathbf{I}-\mathbf{A})^{-1} =s^{-1}\left(\mathbf{I}-s^{-1}\mathbf{A}\right)^{-1} =\sum_{k=0}^{\infty}\mathbf{A}^k s^{-(k+1)}. \]

Thus \( G_i(s)=\sum_{k=0}^{\infty}C_iA_i^kB_i s^{-(k+1)} \). Equality of all Markov parameters gives equality of the Laurent series at infinity. Since both transfer functions are rational functions, this equality on a nonempty analytic region implies equality as rational functions.

Problem 4 (Minimality and Uniqueness): Let \( \Sigma_1 \) and \( \Sigma_2 \) be minimal realizations of the same transfer matrix. Explain why they must have the same state dimension.

Solution: The minimal state dimension is the smallest possible dimension among all realizations of the transfer matrix. If \( \Sigma_1 \) is minimal with dimension \( n_1 \) and \( \Sigma_2 \) is minimal with dimension \( n_2 \), then \( n_1\le n_2 \) because \( \Sigma_1 \) is a realization of the same transfer matrix, and \( n_2\le n_1 \) by the same argument using \( \Sigma_2 \). Hence \( n_1=n_2 \). The uniqueness theorem then strengthens this equality of dimensions to similarity equivalence.

Problem 5 (Hidden Unstable Mode): The nonminimal realization in Problem 2 has hidden eigenvalue \( 5 \). Is the transfer function externally stable? Is the realization internally asymptotically stable?

Solution: The transfer function is \( G(s)=1/(s+1)+1/(s+2) \), so its transfer poles are \( -1 \) and \( -2 \). Hence the zero-state input-output map is externally stable. However, the state matrix of the nonminimal realization has eigenvalues \( -1,-2,5 \). Since \( 5 > 0 \), the realization is not internally asymptotically stable. This is the core reason minimality matters when relating transfer-function behavior to internal state behavior.

Internal equivalence means two realizations differ only by a nonsingular state-coordinate transformation. External equivalence means two realizations have the same transfer matrix or impulse response. Internal equivalence always implies external equivalence, but external equivalence does not imply internal equivalence unless both realizations are minimal. Minimal realization theory therefore provides the precise bridge between input-output descriptions and coordinate-free internal dynamics.