In previous lessons we constructed observable canonical form (OCF), related it to controllable canonical form (CCF), introduced diagonal modal form for distinct eigenvalues, and transformed between physical and modal coordinates. In this lesson we decide when each form should be used. By the end, students should be able to:

• identify which system properties are invariant under similarity transformation;
• explain why OCF is useful for observer-side algebra but often poor numerically;
• explain why modal form is useful for mode interpretation but can be ill-conditioned;
• compute transformation matrices and compare their conditioning;
• choose a coordinate representation for analysis, simulation, and implementation.

We consider a continuous-time SISO or MIMO LTI realization \( (\mathbf{A},\mathbf{B},\mathbf{C},\mathbf{D}) \) with

\[ \dot{\mathbf{x} } = \mathbf{A}\mathbf{x} + \mathbf{B}\mathbf{u}, \qquad \mathbf{y} = \mathbf{C}\mathbf{x} + \mathbf{D}\mathbf{u}. \]

A nonsingular coordinate transformation \( \mathbf{x} = \mathbf{T}\mathbf{z} \) gives

\[ \dot{\mathbf{z} } = \bar{\mathbf{A} }\mathbf{z} + \bar{\mathbf{B} }\mathbf{u}, \qquad \mathbf{y} = \bar{\mathbf{C} }\mathbf{z} + \mathbf{D}\mathbf{u}, \]

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

Choosing a canonical form is therefore choosing a coordinate system. The transfer function matrix is preserved, but internal numerical behavior may change dramatically.

The first practical rule is simple: do not confuse the physics of the system with the appearance of a particular matrix. Similarity changes coordinates, not the underlying finite-dimensional linear system.

Characteristic polynomial invariance. Since

\[ s\mathbf{I}-\bar{\mathbf{A} } = s\mathbf{I}-\mathbf{T}^{-1}\mathbf{A}\mathbf{T} = \mathbf{T}^{-1}(s\mathbf{I}-\mathbf{A})\mathbf{T}, \]

we have

\[ \det(s\mathbf{I}-\bar{\mathbf{A} })= \det(\mathbf{T}^{-1})\det(s\mathbf{I}-\mathbf{A}) \det(\mathbf{T})=\det(s\mathbf{I}-\mathbf{A}). \]

Hence poles/eigenvalues are invariant. For the transfer function, assuming zero initial condition,

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

Using the identity above,

\[ (s\mathbf{I}-\bar{\mathbf{A} })^{-1} = \mathbf{T}^{-1}(s\mathbf{I}-\mathbf{A})^{-1}\mathbf{T}, \]

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

Therefore poles, transfer function, controllability rank, observability rank, minimality, and stability are invariant. However, the entries of \( \bar{\mathbf{A} },\bar{\mathbf{B} }, \bar{\mathbf{C} } \), the numerical condition of the transformation, and the ease of interpretation are not invariant.

The same minimal SISO system may be represented in physical coordinates, observable canonical coordinates, or modal coordinates:

\[ \text{physical coordinates} \quad \longleftrightarrow \quad \text{observable canonical coordinates} \quad \longleftrightarrow \quad \text{modal coordinates}. \]

For a SISO observable system with characteristic polynomial

\[ p(s)=s^n+a_{n-1}s^{n-1}+\cdots+a_1s+a_0, \]

a common OCF convention is

\[ \mathbf{A}_o=\begin{bmatrix} 0 & 0 & \cdots & 0 & -a_0\\ 1 & 0 & \cdots & 0 & -a_1\\ 0 & 1 & \cdots & 0 & -a_2\\ \vdots & \vdots & \ddots & \vdots & \vdots\\ 0 & 0 & \cdots & 1 & -a_{n-1} \end{bmatrix}, \qquad \mathbf{C}_o=\begin{bmatrix}0 & 0 & \cdots & 1\end{bmatrix}. \]

This form makes observability transparent because \( \mathbf{C}_o \) selects the final canonical state, and the rows generated by \( \mathbf{C}_o\mathbf{A}_o^k \) recover the companion structure. By contrast, if \( \mathbf{A} \) has distinct eigenvalues \( \lambda_1,\dots,\lambda_n \) and eigenvector matrix \( \mathbf{V}=[\mathbf{v}_1\;\cdots\;\mathbf{v}_n] \), modal coordinates use

\[ \mathbf{A}\mathbf{V}=\mathbf{V}\boldsymbol{\Lambda}, \qquad \boldsymbol{\Lambda}=\operatorname{diag}(\lambda_1, \dots,\lambda_n), \qquad \mathbf{x}=\mathbf{V}\mathbf{z}. \]

Then \( \dot{z}_i=\lambda_i z_i+\tilde{b}_i u \) in the diagonal case, so each coordinate is directly associated with a mode. This is excellent for interpretation but can be fragile when the eigenvectors are nearly linearly dependent.

The transformation matrix \( \mathbf{T} \) is not a harmless object numerically. Let the computed coordinate vector be \( \mathbf{z}+\delta\mathbf{z} \). Since \( \mathbf{x}=\mathbf{T}\mathbf{z} \), a relative perturbation satisfies the standard bound

\[ \frac{\|\delta\mathbf{z}\|}{\|\mathbf{z}\|} \lesssim \kappa(\mathbf{T}) \frac{\|\delta\mathbf{x}\|}{\|\mathbf{x}\|}, \qquad \kappa(\mathbf{T})=\|\mathbf{T}\|\,\|\mathbf{T}^{-1}\|. \]

Thus a coordinate basis with large \( \kappa(\mathbf{T}) \) magnifies sensor noise, rounding error, and parameter uncertainty in the transformed state.

First-order perturbation of a transformed state matrix.

Suppose \( \mathbf{T} \) is perturbed to \( \mathbf{T}+\Delta\mathbf{T} \). To first order,

\[ (\mathbf{T}+\Delta\mathbf{T})^{-1} \approx \mathbf{T}^{-1}-\mathbf{T}^{-1} \Delta\mathbf{T}\mathbf{T}^{-1}. \]

Therefore

\[ \Delta\bar{\mathbf{A} } \approx -\mathbf{T}^{-1}\Delta\mathbf{T}\bar{\mathbf{A} }+ \mathbf{T}^{-1}\Delta\mathbf{A}\mathbf{T}+ \bar{\mathbf{A} }\mathbf{T}^{-1}\Delta\mathbf{T}. \]

This expression explains why canonical transformations may be poor in finite precision even though they are exact algebraically. If \( \mathbf{T} \) or \( \mathbf{T}^{-1} \) has large norm, the same physical uncertainty produces a much larger uncertainty in the canonical matrices.

OCF is usually the right coordinate system for symbolic observer-side derivations, transfer-function-to-state-space conversion, and proofs based on duality with controllable canonical form. It is usually not the right coordinate system for high-order numerical simulation or embedded implementation unless the realization is well scaled.

Given an observable pair \( (\mathbf{A},\mathbf{C}) \), define the observability matrix

\[ \mathcal{O}(\mathbf{A},\mathbf{C})= \begin{bmatrix} \mathbf{C}\\ \mathbf{C}\mathbf{A}\\ \vdots\\ \mathbf{C}\mathbf{A}^{n-1} \end{bmatrix}. \]

If the desired canonical pair is \( (\mathbf{A}_o,\mathbf{C}_o) \), then the transformation \( \mathbf{x}=\mathbf{T}\mathbf{z} \) satisfies

\[ \mathcal{O}(\mathbf{A}_o,\mathbf{C}_o)= \mathcal{O}(\mathbf{A},\mathbf{C})\mathbf{T}, \qquad \mathbf{T}=\mathcal{O}(\mathbf{A},\mathbf{C})^{-1} \mathcal{O}(\mathbf{A}_o,\mathbf{C}_o). \]

This formula is exact, but it also reveals the danger: if the observability matrix is ill-conditioned, OCF construction amplifies error. In practice, check

\[ \operatorname{rank}\mathcal{O}=n, \qquad \kappa(\mathcal{O})\ \text{moderate}, \qquad \kappa(\mathbf{T})\ \text{moderate}. \]

OCF is especially sensitive to polynomial coefficients. For large \( n \), the coefficients \( a_0,\dots,a_{n-1} \) may vary over many orders of magnitude even when the eigenvalues are physically reasonable. This is why high-order transfer functions are frequently realized by cascaded low-order sections rather than one large companion matrix.

Modal form is most useful when the objective is interpretation: natural modes, dominant poles, slow/fast separation, oscillatory components, and qualitative stability. If \( \mathbf{A} \) is diagonalizable,

\[ \mathbf{A}=\mathbf{V}\boldsymbol{\Lambda}\mathbf{V}^{-1}, \qquad e^{\mathbf{A}t}= \mathbf{V}e^{\boldsymbol{\Lambda}t}\mathbf{V}^{-1}. \]

Thus the homogeneous response is

\[ \mathbf{x}(t)=\sum_{i=1}^n \alpha_i e^{\lambda_i t}\mathbf{v}_i, \qquad \boldsymbol{\alpha}=\mathbf{V}^{-1}\mathbf{x}(0). \]

For a real system with a complex conjugate pair \( \lambda=\sigma+j\omega \) and \( \bar{\lambda}=\sigma-j\omega \), one may use a real modal block instead of complex coordinates:

\[ \mathbf{A}_r=\begin{bmatrix} \sigma & \omega\\ -\omega & \sigma \end{bmatrix}, \qquad e^{\mathbf{A}_r t}=e^{\sigma t} \begin{bmatrix} \cos\omega t & \sin\omega t\\ -\sin\omega t & \cos\omega t \end{bmatrix}. \]

The modal form should not be selected merely because it diagonalizes \( \mathbf{A} \). Repeated or clustered eigenvalues can make \( \mathbf{V} \) nearly singular. If \( \kappa(\mathbf{V}) \) is large, a small perturbation in \( \mathbf{A} \) can cause a large change in the computed eigenvectors. In such cases, a real Schur form, which uses an orthogonal basis, is often better for numerical work, although it is not a canonical form in the strict algebraic sense.

The following flow summarizes a practical decision rule. It is not a theorem; it is an engineering procedure based on exact invariance results and numerical conditioning diagnostics.

flowchart TD
  A["Start with physical realization"] --> B["Need proof or symbolic observer algebra?"]
  B -->|yes| C["Check observability rank and condition"]
  C --> D{"Observable and \nwell conditioned?"}
  D -->|yes| E["Use OCF for derivation; \navoid high-order \nimplementation"]
  D -->|no| F["Do not force OCF; \nrescale or \ndecompose first"]
  B -->|no| G["Need mode interpretation \nor dominant modes?"]
  G -->|yes| H["Compute eigenvectors \nand condition of V"]
  H --> I{"V well conditioned?"}
  I -->|yes| J["Use modal or \nreal modal form"]
  I -->|no| K["Use physical \ncoordinates or orthogonal \nSchur form"]
  G -->|no| L["Keep physical / \nscaled coordinates"]
  E --> M["Document transformation T"]
  F --> M
  J --> M
  K --> M
  L --> M
        

In implementation, a robust workflow is: first scale physical states, then test controllability/observability, then compute the candidate transforms, then compare \( \kappa(\mathbf{T}) \), \( \kappa(\mathcal{O}) \), and \( \kappa(\mathbf{V}) \). The final representation should serve the task, not the other way around.

For a candidate coordinate basis \( \mathbf{T}_q \), define a simple representation score

\[ J(q)=w_1\log\kappa(\mathbf{T}_q)+w_2\rho_q+w_3\eta_q-w_4\iota_q. \]

Here \( \rho_q \) penalizes poor scaling of matrix entries, \( \eta_q \) penalizes loss of physical interpretability, and \( \iota_q \) rewards the information exposed by the representation: observability structure for OCF, modal separation for modal form, or physical meaning for original coordinates. Smaller \( J(q) \) is better. A simple scaling penalty is

\[ \rho_q=\log\left( \frac{\max_{i,j} |\bar{a}_{ij}|}{\min_{i,j}:\bar{a}_{ij}\ne 0 |\bar{a}_{ij}|}\right). \]

This formula is not unique, but it forces students to ask the right questions: Is the form mathematically valid? Is the transformation well-conditioned? Does the form reveal the property needed for the design task?

Canonical transformations are supported differently across languages. For teaching, it is useful to implement the transformations from observability matrices and eigenvectors. For production computation, prefer mature numerical libraries.

• Python: numpy for matrix algebra, scipy.linalg for eigenvalue and Schur routines, python-control for state-space utilities.
• C++: Eigen, Armadillo, Blaze, or LAPACK bindings for reliable inverses, decompositions, and condition estimates.
• Java: EJML, JBLAS, or Apache Commons Math for matrix operations and eigensolvers.
• MATLAB/Simulink: Control System Toolbox functions such as ss, canon, obsv, ctrb, eig, and Simulink State-Space blocks.
• Wolfram Mathematica: Eigensystem, MatrixExp, JordanDecomposition, and symbolic polynomial operations.

The examples below use a third-order system with \( p(s)=s^3+6s^2+11s+6=(s+1)(s+2)(s+3) \). The system is deliberately scaled so that students can see why algebraic canonical transformations may be unattractive numerically even for a stable, observable system.

Chapter17_Lesson5.py

# Chapter17_Lesson5.py
# Practical Considerations in Choosing Canonical Forms
# Requires: numpy, scipy (optional for matrix exponential demo)

import numpy as np

np.set_printoptions(precision=6, suppress=True)


def inv(M: np.ndarray) -> np.ndarray:
    """Numerically invert a square matrix."""
    return np.linalg.inv(M)


def observability_matrix(A: np.ndarray, C: np.ndarray) -> np.ndarray:
    """Build O = [C; C A; ...; C A^(n-1)]."""
    n = A.shape[0]
    rows = []
    Ak = np.eye(n)
    for _ in range(n):
        rows.append(C @ Ak)
        Ak = Ak @ A
    return np.vstack(rows)


def frobenius_condition(M: np.ndarray) -> float:
    """A simple conditioning proxy: ||M||_F ||M^{-1}||_F."""
    return np.linalg.norm(M, ord="fro") * np.linalg.norm(inv(M), ord="fro")


def transform_system(A: np.ndarray, B: np.ndarray, C: np.ndarray, T: np.ndarray):
    """For x = T z, return z-dot = Abar z + Bbar u, y = Cbar z."""
    Ti = inv(T)
    Abar = Ti @ A @ T
    Bbar = Ti @ B
    Cbar = C @ T
    return Abar, Bbar, Cbar


def choose_form(kappa_ocf: float, kappa_modal: float, observable_rank: int, n: int) -> str:
    """Rule-based recommendation for instructional use."""
    if observable_rank < n:
        return "Do not use OCF before separating the unobservable subspace."
    if kappa_modal < 50 and kappa_modal <= kappa_ocf:
        return "Modal form is preferred for mode interpretation and decoupled simulation."
    if kappa_ocf < 50 and kappa_ocf < kappa_modal:
        return "OCF is acceptable for observer-side algebra and transfer-function realization."
    return "Use the physical/scaled coordinates or a numerically orthogonal form instead."


def main() -> None:
    # Companion model with p(s) = s^3 + 6 s^2 + 11 s + 6 = (s+1)(s+2)(s+3)
    A0 = np.array([[0.0, 1.0, 0.0],
                   [0.0, 0.0, 1.0],
                   [-6.0, -11.0, -6.0]])
    B0 = np.array([[0.0], [0.0], [1.0]])
    C0 = np.array([[1.0, 0.0, 0.0]])

    # A deliberately scaled coordinate system: x = S z0.
    S = np.diag([1.0, 0.05, 20.0])
    A = S @ A0 @ inv(S)
    B = S @ B0
    C = C0 @ inv(S)

    print("Physical/scaled A:")
    print(A)
    print("B =", B.T)
    print("C =", C)

    # Observable canonical target: A_o = A0^T, C_o = [0 0 1].
    Ao = A0.T
    Co = np.array([[0.0, 0.0, 1.0]])
    O = observability_matrix(A, C)
    Oo = observability_matrix(Ao, Co)
    T_ocf = inv(O) @ Oo
    A_ocf, B_ocf, C_ocf = transform_system(A, B, C, T_ocf)

    print("\nRank of observability matrix:", np.linalg.matrix_rank(O))
    print("Frobenius condition proxy of OCF transform:", frobenius_condition(T_ocf))
    print("A in observable canonical coordinates:")
    print(A_ocf)
    print("C in observable canonical coordinates:", C_ocf)

    # Modal transform: right eigenvectors of A.
    eigvals, V = np.linalg.eig(A)
    A_modal, B_modal, C_modal = transform_system(A, B, C, V)

    print("\nEigenvalues:", eigvals)
    print("Frobenius condition proxy of modal eigenvector matrix:", frobenius_condition(V))
    print("A in modal coordinates:")
    print(A_modal)

    print("\nRecommendation:")
    print(choose_form(frobenius_condition(T_ocf), frobenius_condition(V), np.linalg.matrix_rank(O), A.shape[0]))


if __name__ == "__main__":
    main()

Chapter17_Lesson5.cpp

// Chapter17_Lesson5.cpp
// Practical Considerations in Choosing Canonical Forms
// Compile: g++ -std=c++17 Chapter17_Lesson5.cpp -o Chapter17_Lesson5

#include <cmath>
#include <iomanip>
#include <iostream>
#include <stdexcept>
#include <string>
#include <vector>

using Matrix = std::vector<std::vector<double>>;

Matrix eye(int n) {
    Matrix I(n, std::vector<double>(n, 0.0));
    for (int i = 0; i < n; ++i) I[i][i] = 1.0;
    return I;
}

Matrix matmul(const Matrix& A, const Matrix& B) {
    int r = (int)A.size(), k = (int)B.size(), c = (int)B[0].size();
    Matrix M(r, std::vector<double>(c, 0.0));
    for (int i = 0; i < r; ++i)
        for (int j = 0; j < c; ++j)
            for (int t = 0; t < k; ++t)
                M[i][j] += A[i][t] * B[t][j];
    return M;
}

Matrix transpose(const Matrix& A) {
    int r = (int)A.size(), c = (int)A[0].size();
    Matrix T(c, std::vector<double>(r));
    for (int i = 0; i < r; ++i)
        for (int j = 0; j < c; ++j)
            T[j][i] = A[i][j];
    return T;
}

Matrix inverse(Matrix A) {
    int n = (int)A.size();
    Matrix Aug(n, std::vector<double>(2 * n, 0.0));
    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 i = col + 1; i < n; ++i)
            if (std::abs(Aug[i][col]) > std::abs(Aug[pivot][col])) pivot = i;
        if (std::abs(Aug[pivot][col]) < 1e-12) throw std::runtime_error("Singular matrix");
        std::swap(Aug[pivot], Aug[col]);
        double div = Aug[col][col];
        for (int j = 0; j < 2 * n; ++j) Aug[col][j] /= div;
        for (int i = 0; i < n; ++i) {
            if (i == col) continue;
            double factor = Aug[i][col];
            for (int j = 0; j < 2 * n; ++j) Aug[i][j] -= factor * Aug[col][j];
        }
    }
    Matrix Inv(n, std::vector<double>(n));
    for (int i = 0; i < n; ++i)
        for (int j = 0; j < n; ++j)
            Inv[i][j] = Aug[i][n + j];
    return Inv;
}

double frobeniusNorm(const Matrix& A) {
    double s = 0.0;
    for (const auto& row : A)
        for (double v : row) s += v * v;
    return std::sqrt(s);
}

double frobeniusCondition(const Matrix& A) {
    return frobeniusNorm(A) * frobeniusNorm(inverse(A));
}

Matrix observabilityMatrix(const Matrix& A, const Matrix& C) {
    int n = (int)A.size();
    Matrix O;
    Matrix Ak = eye(n);
    for (int k = 0; k < n; ++k) {
        Matrix row = matmul(C, Ak);
        O.push_back(row[0]);
        Ak = matmul(Ak, A);
    }
    return O;
}

void printMatrix(const std::string& name, const Matrix& A) {
    std::cout << name << "\n";
    for (const auto& row : A) {
        for (double v : row) std::cout << std::setw(12) << std::setprecision(6) << std::fixed << v << " ";
        std::cout << "\n";
    }
}

Matrix diag(const std::vector<double>& d) {
    int n = (int)d.size();
    Matrix D(n, std::vector<double>(n, 0.0));
    for (int i = 0; i < n; ++i) D[i][i] = d[i];
    return D;
}

int main() {
    Matrix A0 = { {0, 1, 0}, {0, 0, 1}, {-6, -11, -6} };
    Matrix B0 = { {0}, {0}, {1} };
    Matrix C0 = { {1, 0, 0} };
    Matrix S = diag({1.0, 0.05, 20.0});
    Matrix A = matmul(matmul(S, A0), inverse(S));
    Matrix B = matmul(S, B0);
    Matrix C = matmul(C0, inverse(S));

    Matrix Ao = transpose(A0);
    Matrix Co = { {0, 0, 1} };
    Matrix O = observabilityMatrix(A, C);
    Matrix Oo = observabilityMatrix(Ao, Co);
    Matrix T_ocf = matmul(inverse(O), Oo);
    Matrix Ti_ocf = inverse(T_ocf);
    Matrix A_ocf = matmul(matmul(Ti_ocf, A), T_ocf);
    Matrix C_ocf = matmul(C, T_ocf);

    // Modal eigenvectors are known for the companion model: [1, lambda, lambda^2]^T.
    // For A = S A0 S^{-1}, the corresponding eigenvectors are S times those vectors.
    std::vector<double> lam = {-1.0, -2.0, -3.0};
    Matrix V(3, std::vector<double>(3));
    for (int j = 0; j < 3; ++j) {
        double l = lam[j];
        Matrix v0 = { {1.0}, {l}, {l * l} };
        Matrix v = matmul(S, v0);
        for (int i = 0; i < 3; ++i) V[i][j] = v[i][0];
    }
    Matrix A_modal = matmul(matmul(inverse(V), A), V);

    printMatrix("Physical/scaled A:", A);
    printMatrix("Observable canonical A:", A_ocf);
    printMatrix("Observable canonical C:", C_ocf);
    std::cout << "OCF transform condition proxy = " << frobeniusCondition(T_ocf) << "\n\n";
    printMatrix("Modal A:", A_modal);
    std::cout << "Modal eigenvector condition proxy = " << frobeniusCondition(V) << "\n";

    double kModal = frobeniusCondition(V);
    double kOcf = frobeniusCondition(T_ocf);
    if (std::min(kModal, kOcf) > 5000.0)
        std::cout << "Recommendation: both canonical transforms are poorly conditioned; keep physical/scaled coordinates or use an orthogonal Schur form.\n";
    else if (kModal < kOcf)
        std::cout << "Recommendation: modal coordinates are cleaner for this example.\n";
    else
        std::cout << "Recommendation: OCF is acceptable for observer algebra in this example.\n";
    return 0;
}

Chapter17_Lesson5.java

// Chapter17_Lesson5.java
// Practical Considerations in Choosing Canonical Forms
// Compile: javac Chapter17_Lesson5.java
// Run:     java Chapter17_Lesson5

public class Chapter17_Lesson5 {
    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[][] matmul(double[][] A, double[][] B) {
        int r = A.length, k = B.length, c = B[0].length;
        double[][] M = new double[r][c];
        for (int i = 0; i < r; i++)
            for (int j = 0; j < c; j++)
                for (int t = 0; t < k; t++)
                    M[i][j] += A[i][t] * B[t][j];
        return M;
    }

    static double[][] transpose(double[][] A) {
        int r = A.length, c = A[0].length;
        double[][] T = new double[c][r];
        for (int i = 0; i < r; i++)
            for (int j = 0; j < c; j++)
                T[j][i] = A[i][j];
        return T;
    }

    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 i = col + 1; i < n; i++)
                if (Math.abs(Aug[i][col]) > Math.abs(Aug[pivot][col])) pivot = i;
            if (Math.abs(Aug[pivot][col]) < 1e-12) throw new RuntimeException("Singular matrix");
            double[] temp = Aug[pivot]; Aug[pivot] = Aug[col]; Aug[col] = temp;
            double div = Aug[col][col];
            for (int j = 0; j < 2 * n; j++) Aug[col][j] /= div;
            for (int i = 0; i < n; i++) {
                if (i == col) continue;
                double factor = Aug[i][col];
                for (int j = 0; j < 2 * n; j++) Aug[i][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 double frobeniusNorm(double[][] A) {
        double s = 0.0;
        for (double[] row : A)
            for (double v : row) s += v * v;
        return Math.sqrt(s);
    }

    static double frobeniusCondition(double[][] A) {
        return frobeniusNorm(A) * frobeniusNorm(inverse(A));
    }

    static double[][] observabilityMatrix(double[][] A, double[][] C) {
        int n = A.length;
        double[][] O = new double[n][n];
        double[][] Ak = eye(n);
        for (int k = 0; k < n; k++) {
            double[][] row = matmul(C, Ak);
            for (int j = 0; j < n; j++) O[k][j] = row[0][j];
            Ak = matmul(Ak, A);
        }
        return O;
    }

    static double[][] diag(double[] d) {
        int n = d.length;
        double[][] D = new double[n][n];
        for (int i = 0; i < n; i++) D[i][i] = d[i];
        return D;
    }

    static void printMatrix(String name, double[][] A) {
        System.out.println(name);
        for (double[] row : A) {
            for (double v : row) System.out.printf("%12.6f ", v);
            System.out.println();
        }
    }

    public static void main(String[] args) {
        double[][] A0 = { {0, 1, 0}, {0, 0, 1}, {-6, -11, -6} };
        double[][] C0 = { {1, 0, 0} };
        double[][] S = diag(new double[]{1.0, 0.05, 20.0});
        double[][] A = matmul(matmul(S, A0), inverse(S));
        double[][] C = matmul(C0, inverse(S));

        double[][] Ao = transpose(A0);
        double[][] Co = { {0, 0, 1} };
        double[][] O = observabilityMatrix(A, C);
        double[][] Oo = observabilityMatrix(Ao, Co);
        double[][] T_ocf = matmul(inverse(O), Oo);
        double[][] A_ocf = matmul(matmul(inverse(T_ocf), A), T_ocf);
        double[][] C_ocf = matmul(C, T_ocf);

        double[] lam = {-1.0, -2.0, -3.0};
        double[][] V = new double[3][3];
        for (int j = 0; j < 3; j++) {
            double l = lam[j];
            double[][] v0 = { {1.0}, {l}, {l * l} };
            double[][] v = matmul(S, v0);
            for (int i = 0; i < 3; i++) V[i][j] = v[i][0];
        }
        double[][] A_modal = matmul(matmul(inverse(V), A), V);

        printMatrix("Physical/scaled A:", A);
        printMatrix("Observable canonical A:", A_ocf);
        printMatrix("Observable canonical C:", C_ocf);
        System.out.printf("OCF transform condition proxy = %.6f%n%n", frobeniusCondition(T_ocf));
        printMatrix("Modal A:", A_modal);
        System.out.printf("Modal eigenvector condition proxy = %.6f%n", frobeniusCondition(V));

        double kModal = frobeniusCondition(V);
        double kOcf = frobeniusCondition(T_ocf);
        if (Math.min(kModal, kOcf) > 5000.0)
            System.out.println("Recommendation: both canonical transforms are poorly conditioned; keep physical/scaled coordinates or use an orthogonal Schur form.");
        else if (kModal < kOcf)
            System.out.println("Recommendation: modal coordinates are cleaner for this example.");
        else
            System.out.println("Recommendation: OCF is acceptable for observer algebra in this example.");
    }
}

The script below uses base MATLAB operations and optionally calls Control System Toolbox commands if they are available. In Simulink, the matrices \( \mathbf{A},\mathbf{B},\mathbf{C},\mathbf{D} \) can be inserted into a State-Space block; the transformed matrices \( \bar{\mathbf{A} },\bar{\mathbf{B} },\bar{\mathbf{C} },\mathbf{D} \) give the same zero-initial-condition input-output response.

Chapter17_Lesson5.m

% Chapter17_Lesson5.m
% Practical Considerations in Choosing Canonical Forms
% Requires: base MATLAB. Control System Toolbox is optional for ss/canon.

clear; clc;

A0 = [0 1 0; 0 0 1; -6 -11 -6];
B0 = [0; 0; 1];
C0 = [1 0 0];
D0 = 0;

% Scaled coordinates x = S*z0
S = diag([1, 0.05, 20]);
A = S*A0/S;
B = S*B0;
C = C0/S;
D = D0;

% Observable canonical target: transpose of controllable companion form
Ao = A0.';
Co = [0 0 1];
O = obsv_local(A, C);
Oo = obsv_local(Ao, Co);
T_ocf = O \ Oo;
A_ocf = T_ocf \ A * T_ocf;
B_ocf = T_ocf \ B;
C_ocf = C * T_ocf;

fprintf('Rank of observability matrix: %d\n', rank(O));
fprintf('OCF transform condition number: %.6g\n', cond(T_ocf));
disp('A in observable canonical coordinates:'); disp(A_ocf);
disp('C in observable canonical coordinates:'); disp(C_ocf);

% Modal transformation
[V,Lambda] = eig(A);
A_modal = V \ A * V;
B_modal = V \ B;
C_modal = C * V;

fprintf('Modal eigenvector condition number: %.6g\n', cond(V));
disp('A in modal coordinates:'); disp(A_modal);

% Optional Control System Toolbox checks
if exist('ss','file') == 2
    sys = ss(A,B,C,D);
    disp('State-space object created. Use canon(sys,''modal'') if Control System Toolbox supports it.');
    try
        sys_modal = canon(sys, 'modal'); %#ok<NASGU>
        disp('MATLAB canon(sys,''modal'') succeeded.');
    catch ME
        fprintf('canon(sys,''modal'') not available or failed: %s\n', ME.message);
    end
end

if min(cond(V), cond(T_ocf)) > 5000
    disp('Recommendation: both canonical transforms are poorly conditioned; keep physical/scaled coordinates or use an orthogonal Schur form.');
elseif cond(V) < cond(T_ocf)
    disp('Recommendation: modal form is numerically preferable for this example.');
else
    disp('Recommendation: OCF is acceptable for observer algebra in this example.');
end

function O = obsv_local(A,C)
    n = size(A,1);
    O = zeros(n,n);
    Ak = eye(n);
    for k = 1:n
        O(k,:) = C*Ak;
        Ak = Ak*A;
    end
end

Chapter17_Lesson5.nb

Notebook[{Cell["Chapter17_Lesson5.nb", "Title"],
Cell["Practical Considerations in Choosing Canonical Forms", "Subtitle"],
Cell[BoxData[ToBoxes[
  ClearAll[obs, froCond, transformSystem];
  obs[A_, C_] := Module[{n = Length[A], Ak = IdentityMatrix[Length[A]]},
    Table[With[{row = C.Ak}, Ak = Ak.A; row[[1]]], {Length[A]}]
  ];
  froCond[M_] := Norm[M, "Frobenius"] Norm[Inverse[M], "Frobenius"];
  transformSystem[A_, B_, C_, T_] := {Inverse[T].A.T, Inverse[T].B, C.T};
]], "Input"],
Cell[BoxData[ToBoxes[
  A0 = { {0, 1, 0}, {0, 0, 1}, {-6, -11, -6} };
  B0 = { {0}, {0}, {1} };
  C0 = { {1, 0, 0} };
  S = DiagonalMatrix[{1, 0.05, 20}];
  A = S.A0.Inverse[S];
  B = S.B0;
  Cmat = C0.Inverse[S];
  Ao = Transpose[A0];
  Co = { {0, 0, 1} };
  Omat = obs[A, Cmat];
  Oocf = obs[Ao, Co];
  Tocf = Inverse[Omat].Oocf;
  {Aocf, Bocf, Cocf} = transformSystem[A, B, Cmat, Tocf];
  {vals, V} = Eigensystem[A];
  Vmat = Transpose[V];
  {Amodal, Bmodal, Cmodal} = transformSystem[A, B, Cmat, Vmat];
]], "Input"],
Cell[BoxData[ToBoxes[
  Column[{
    "Physical/scaled A:", MatrixForm[A],
    "Observable canonical A:", MatrixForm[Chop[Aocf]],
    "Observable canonical C:", MatrixForm[Chop[Cocf]],
    "OCF transform Frobenius condition proxy:", froCond[Tocf],
    "Modal A:", MatrixForm[Chop[Amodal]],
    "Modal eigenvector Frobenius condition proxy:", froCond[Vmat],
    Which[
      Min[froCond[Vmat], froCond[Tocf]] > 5000,
        "Recommendation: both canonical transforms are poorly conditioned; keep physical/scaled coordinates or use an orthogonal Schur form.",
      froCond[Vmat] < froCond[Tocf],
        "Recommendation: modal coordinates are cleaner for this example.",
      True,
        "Recommendation: OCF is acceptable for observer algebra in this example."]
  }]
]], "Input"]
}]

Problem 1 (Transfer Function Invariance): Let \( \mathbf{x}=\mathbf{T}\mathbf{z} \) with nonsingular \( \mathbf{T} \). Prove that the transfer function is unchanged by the transformation.

Solution: The transformed realization is

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

Then

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

After cancellation, this equals

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

Problem 2 (OCF Transformation from Observability Matrices): Suppose \( (\mathbf{A},\mathbf{C}) \) and \( (\mathbf{A}_o,\mathbf{C}_o) \) are observable pairs with the same characteristic polynomial. Derive the formula for \( \mathbf{T} \) such that \( \mathbf{T}^{-1}\mathbf{A}\mathbf{T}=\mathbf{A}_o \) and \( \mathbf{C}\mathbf{T}=\mathbf{C}_o \).

Solution: Under \( \mathbf{x}=\mathbf{T}\mathbf{z} \),

\[ \mathbf{C}_o\mathbf{A}_o^k = \mathbf{C}\mathbf{T}(\mathbf{T}^{-1}\mathbf{A}\mathbf{T})^k = \mathbf{C}\mathbf{A}^k\mathbf{T}, \qquad k=0,1,\dots,n-1. \]

Stacking these rows gives

\[ \mathcal{O}(\mathbf{A}_o,\mathbf{C}_o)= \mathcal{O}(\mathbf{A},\mathbf{C})\mathbf{T}. \]

Since the pair is observable, \( \mathcal{O}(\mathbf{A},\mathbf{C}) \) is nonsingular, so

\[ \mathbf{T}=\mathcal{O}(\mathbf{A},\mathbf{C})^{-1} \mathcal{O}(\mathbf{A}_o,\mathbf{C}_o). \]

Problem 3 (Modal Coordinates): Let \( \mathbf{A}\mathbf{V}=\mathbf{V}\boldsymbol{\Lambda} \) with nonsingular \( \mathbf{V} \). Show that \( \mathbf{x}=\mathbf{V}\mathbf{z} \) yields \( \dot{\mathbf{z} }=\boldsymbol{\Lambda}\mathbf{z}+ \mathbf{V}^{-1}\mathbf{B}\mathbf{u} \).

Solution: Substitute \( \mathbf{x}=\mathbf{V}\mathbf{z} \):

\[ \mathbf{V}\dot{\mathbf{z} }=\mathbf{A}\mathbf{V}\mathbf{z}+ \mathbf{B}\mathbf{u}=\mathbf{V}\boldsymbol{\Lambda}\mathbf{z}+ \mathbf{B}\mathbf{u}. \]

Multiplying by \( \mathbf{V}^{-1} \) gives

\[ \dot{\mathbf{z} }=\boldsymbol{\Lambda}\mathbf{z}+ \mathbf{V}^{-1}\mathbf{B}\mathbf{u}. \]

Problem 4 (Conditioning Interpretation): A modal transformation has \( \kappa(\mathbf{V})=10^7 \). A sensor reconstruction error in physical coordinates has relative size \( 10^{-6} \). Estimate the possible relative error in modal coordinates.

Solution: The first-order bound gives

\[ \frac{\|\delta\mathbf{z}\|}{\|\mathbf{z}\|} \lesssim \kappa(\mathbf{V}) \frac{\|\delta\mathbf{x}\|}{\|\mathbf{x}\|} =10^7\cdot 10^{-6}=10. \]

A relative error of order ten means the modal state can be essentially unreliable. The eigenvalues may still be accurate, but the modal coordinates should not be used for implementation.

Problem 5 (Decision Flow): A third-order observable system has distinct real eigenvalues, but both the observability matrix and eigenvector matrix are badly conditioned. Which form should be used for numerical simulation?

Solution: The correct choice is usually the original physical/scaled coordinates or an orthogonal Schur representation, not OCF or modal form. The reasoning is summarized below.

flowchart TD
  A["Observable system with distinct eigenvalues"] --> B["Check condition of observability matrix"]
  A --> C["Check condition of eigenvector matrix"]
  B --> D{"OCF transform reliable?"}
  C --> E{"Modal transform reliable?"}
  D -->|no| F["Reject OCF for simulation"]
  E -->|no| G["Reject modal coordinates for implementation"]
  F --> H["Use physical/scaled coordinates or Schur form"]
  G --> H
        

Similarity transformations preserve the system's input-output behavior, eigenvalues, controllability, observability, and minimality, but they do not preserve numerical conditioning or physical interpretability. OCF is valuable for observer-side symbolic derivations and transfer-function realization, while modal form is valuable for interpreting natural modes. For high-order or poorly scaled systems, both can be numerically fragile. In practical modern control, the best coordinate system is the one that exposes the needed structure while keeping the transformation well-conditioned and the state interpretation defensible.