In previous chapters, a state-space realization of a SISO LTI system was written as \( \dot{\mathbf{x}}=\mathbf{A}\mathbf{x}+\mathbf{B}u \), \( y=\mathbf{C}\mathbf{x}+Du \). A realization is not unique: if \( \mathbf{x}=\mathbf{P}\mathbf{z} \) with nonsingular \( \mathbf{P} \), then

\[ \dot{\mathbf{z}}=\mathbf{P}^{-1}\mathbf{A}\mathbf{P}\mathbf{z} +\mathbf{P}^{-1}\mathbf{B}u,\qquad y=\mathbf{C}\mathbf{P}\mathbf{z}+Du . \]

The observable canonical form is a special coordinate representation whose output equation is extremely simple and whose observability is explicit from its matrix pattern. It is the natural dual of the controllable canonical form studied in Chapter 16.

flowchart TD
  A["Start with G(s)=N(s)/D(s)"] --> B["Normalize D(s) so leading coefficient is 1"]
  B --> C["Separate direct term D_o if numerator degree equals denominator degree"]
  C --> D["Build observer companion matrix A_o"]
  D --> E["Place numerator coefficients into B_o"]
  E --> F["Set C_o = [0 ... 0 1]"]
  F --> G["Verify rank of observability matrix is n"]
  G --> H["Use realization: xdot=A_o x+B_o u, y=C_o x+D_o u"]
        

Consider a proper SISO transfer function \( G(s)=N(s)/D(s) \) with monic denominator of degree \( n \):

\[ D(s)=s^n+a_1s^{n-1}+a_2s^{n-2}+\cdots+a_{n-1}s+a_n . \]

For the strictly proper part, write the numerator as

\[ N_{sp}(s)=b_0s^{n-1}+b_1s^{n-2}+\cdots+b_{n-2}s+b_{n-1}. \]

If the original transfer function is proper but not strictly proper, first divide:

\[ G(s)=d+\frac{N_{sp}(s)}{D(s)},\qquad \deg N_{sp} < n . \]

The scalar \( d \) becomes the feedthrough term \( D_o \). The canonical state matrices below realize only the strictly proper part; the feedthrough is then added to the output equation.

With the polynomial convention above, the observable canonical realization is

\[ \mathbf{A}_o= \begin{bmatrix} 0 & 0 & \cdots & 0 & -a_n\\ 1 & 0 & \cdots & 0 & -a_{n-1}\\ 0 & 1 & \cdots & 0 & -a_{n-2}\\ \vdots & \vdots & \ddots & \vdots & \vdots\\ 0 & 0 & \cdots & 1 & -a_1 \end{bmatrix}, \quad \mathbf{B}_o= \begin{bmatrix} b_{n-1}\\ b_{n-2}\\ \vdots\\ b_1\\ b_0 \end{bmatrix}, \quad \mathbf{C}_o= \begin{bmatrix} 0 & 0 & \cdots & 0 & 1 \end{bmatrix}, \quad D_o=d . \]

Notice that \( \mathbf{C}_o \) simply reads the last canonical state. The denominator coefficients appear in the final column of \( \mathbf{A}_o \), while the numerator coefficients appear in \( \mathbf{B}_o \) in reverse order from constant term to highest strictly proper term.

This convention is the transpose-dual of the controllable companion realization. Some texts use the opposite state ordering; the transfer function is unchanged, but the displayed positions of the coefficients differ.

The observability matrix for a pair \( (\mathbf{C},\mathbf{A}) \) is

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

For \( \mathbf{C}_o=[0\;0\;\cdots\;0\;1] \), the first row selects the last state. Multiplication by \( \mathbf{A}_o \) shifts this selection leftward while adding combinations of already observed columns through the final-column coefficients. Consequently, the rows of \( \mathcal{O}(\mathbf{A}_o,\mathbf{C}_o) \) have the anti-triangular structure

\[ \mathcal{O}(\mathbf{A}_o,\mathbf{C}_o)= \begin{bmatrix} 0 & 0 & \cdots & 0 & 1\\ 0 & 0 & \cdots & 1 & *\\ \vdots & \vdots & \cdots & \vdots & \vdots\\ 0 & 1 & \cdots & * & *\\ 1 & * & \cdots & * & * \end{bmatrix}. \]

The anti-diagonal entries are all equal to one. Hence

\[ \det\mathcal{O}(\mathbf{A}_o,\mathbf{C}_o)=(-1)^{n(n-1)/2}\neq 0, \qquad \operatorname{rank}\mathcal{O}(\mathbf{A}_o,\mathbf{C}_o)=n . \]

Therefore the OCF realization is observable for every monic denominator polynomial, independent of the numerator coefficients.

The transfer function of any state-space realization is

\[ G_o(s)=\mathbf{C}_o(s\mathbf{I}-\mathbf{A}_o)^{-1}\mathbf{B}_o+D_o . \]

Because \( \mathbf{A}_o \) is the transpose of the controllable companion matrix, it has the same characteristic polynomial:

\[ \det(s\mathbf{I}-\mathbf{A}_o)=s^n+a_1s^{n-1}+\cdots+a_n=D(s). \]

The vector \( \mathbf{B}_o \) was chosen so that the numerator generated by \( \mathbf{C}_o\operatorname{adj}(s\mathbf{I}-\mathbf{A}_o)\mathbf{B}_o \) is exactly \( N_{sp}(s) \). Thus

\[ \mathbf{C}_o(s\mathbf{I}-\mathbf{A}_o)^{-1}\mathbf{B}_o = \frac{b_0s^{n-1}+b_1s^{n-2}+\cdots+b_{n-1}} {s^n+a_1s^{n-1}+\cdots+a_n}. \]

Adding \( D_o=d \) gives the original proper transfer function. Therefore the OCF is not merely a convenient matrix pattern; it is a realization of the same input-output map.

Let \( (\mathbf{A}_c,\mathbf{B}_c,\mathbf{C}_c,D) \) be the controllable canonical realization from the previous chapter. The observable canonical realization can be viewed as the dual system:

\[ \mathbf{A}_o=\mathbf{A}_c^{T},\qquad \mathbf{B}_o=\mathbf{C}_c^{T},\qquad \mathbf{C}_o=\mathbf{B}_c^{T},\qquad D_o=D . \]

The transfer function is invariant under scalar transpose:

\[ \mathbf{C}_c(s\mathbf{I}-\mathbf{A}_c)^{-1}\mathbf{B}_c = \left[ \mathbf{B}_c^{T}(s\mathbf{I}-\mathbf{A}_c^{T})^{-1}\mathbf{C}_c^{T} \right]^{T}. \]

Since the transfer function is scalar, the transpose does not alter its value. This is the algebraic reason why CCF and OCF realize the same transfer function while emphasizing different structural properties.

flowchart TD
  CCF["Controllable canonical form"] -->|transpose A and swap B/C| OCF["Observable canonical form"]
  CCF --> C1["Rank controllability matrix is obvious"]
  OCF --> O1["Rank observability matrix is obvious"]
  C1 --> SAME["Same scalar transfer function"]
  O1 --> SAME
        

Suppose an observable realization \( (\mathbf{A},\mathbf{B},\mathbf{C},D) \) is already available. We seek a nonsingular matrix \( \mathbf{P} \) such that \( \mathbf{x}=\mathbf{P}\mathbf{z} \) gives

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

Define the two observability matrices

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

Since \( \mathbf{C}_o\mathbf{A}_o^k=\mathbf{C}\mathbf{A}^k\mathbf{P} \), we have \( \mathcal{O}_o=\mathcal{O}\mathbf{P} \). Therefore

\[ \boxed{\mathbf{P}=\mathcal{O}^{-1}\mathcal{O}_o} \]

This formula is valid only when the original system is observable, so that \( \mathcal{O}^{-1} \) exists. Numerically, avoid explicitly forming the inverse; solve the linear system \( \mathcal{O}\mathbf{P}=\mathcal{O}_o \).

Consider

\[ G(s)=\frac{2s^2+5s+3}{s^3+4s^2+6s+8}. \]

Here \( a_1=4,\;a_2=6,\;a_3=8 \) and \( b_0=2,\;b_1=5,\;b_2=3 \). Therefore

\[ \mathbf{A}_o= \begin{bmatrix} 0 & 0 & -8\\ 1 & 0 & -6\\ 0 & 1 & -4 \end{bmatrix}, \quad \mathbf{B}_o= \begin{bmatrix} 3\\5\\2 \end{bmatrix}, \quad \mathbf{C}_o= \begin{bmatrix} 0 & 0 & 1 \end{bmatrix}, \quad D_o=0 . \]

The observability matrix is

\[ \mathcal{O}_o= \begin{bmatrix} 0 & 0 & 1\\ 0 & 1 & -4\\ 1 & -4 & 10 \end{bmatrix}, \qquad \det(\mathcal{O}_o)=-1 . \]

Hence the canonical realization is observable. Direct substitution into \( \mathbf{C}_o(s\mathbf{I}-\mathbf{A}_o)^{-1}\mathbf{B}_o \) recovers the stated rational function.

Canonical forms are excellent for theory, hand calculation, and introductory algorithms. They are not always ideal for high-order numerical computation because polynomial coefficients and observability matrices may be poorly conditioned. A useful warning indicator is

\[ \kappa(\mathcal{O})=\|\mathcal{O}\|\,\|\mathcal{O}^{-1}\|. \]

Large \( \kappa(\mathcal{O}) \) means that the coordinate transformation to OCF may amplify roundoff errors. In practical software, prefer numerically stable state-space algorithms for high-order models, and use canonical forms mainly when their structure is analytically useful.

Useful libraries: Python (numpy, scipy.signal, python-control), C++ (Eigen, Armadillo), Java (EJML, Apache Commons Math), MATLAB/Simulink (Control System Toolbox, State-Space block), and Wolfram Mathematica (StateSpaceModel, symbolic matrix tools).

Chapter17_Lesson1.py

"""
Chapter17_Lesson1.py

Observable Canonical Form (OCF) construction for a SISO continuous-time
transfer function

        G(s) = N(s) / D(s)

with monic denominator

        D(s) = s^n + a1 s^(n-1) + ... + an.

This script constructs the observer/observable companion realization

        x_dot = A_o x + B_o u
        y     = C_o x + D_o u

using the convention

        A_o = A_c^T,   C_o = [0 ... 0 1],

where A_c is the controllable companion matrix.
"""

from __future__ import annotations

import numpy as np


def trim_leading_zeros(coeffs, tol=1e-12):
    """Remove leading coefficients that are numerically zero."""
    coeffs = list(map(float, coeffs))
    while len(coeffs) > 1 and abs(coeffs[0]) < tol:
        coeffs.pop(0)
    return np.array(coeffs, dtype=float)


def observable_canonical_form(den, num):
    """
    Build OCF from denominator and numerator coefficients in descending powers.

    Parameters
    ----------
    den : sequence of float
        [1, a1, a2, ..., an] after normalization; non-monic input is allowed.
    num : sequence of float
        Numerator coefficients in descending powers. It may be strictly proper
        or proper. If degree(num) == degree(den), the direct term D_o is
        separated by polynomial division.

    Returns
    -------
    A_o, B_o, C_o, D_o : numpy arrays
        Observable canonical realization.
    """
    den = trim_leading_zeros(den)
    num = trim_leading_zeros(num)

    if den[0] == 0:
        raise ValueError("The denominator leading coefficient must be nonzero.")

    # Normalize denominator to be monic.
    num = num / den[0]
    den = den / den[0]

    n = len(den) - 1
    if n <= 0:
        raise ValueError("The denominator must have positive degree.")

    if len(num) > n + 1:
        raise ValueError("This implementation expects a proper transfer function.")

    # Pad numerator to length n+1, aligned with powers s^n, ..., s, 1.
    num_pad = np.zeros(n + 1)
    num_pad[-len(num):] = num

    # Separate direct feedthrough if numerator degree is n.
    D_o = float(num_pad[0])
    remainder = num_pad - D_o * den
    # Now remainder[0] is zero, and remainder[1:] corresponds to
    # coefficients of s^(n-1), ..., s, 1.
    numerator_without_direct = remainder[1:]

    a = den[1:]  # [a1, ..., an]

    A_o = np.zeros((n, n))
    for i in range(1, n):
        A_o[i, i - 1] = 1.0
    A_o[:, -1] = -a[::-1]  # [-an, ..., -a1]^T

    B_o = numerator_without_direct[::-1].reshape((n, 1))  # [constant, ..., s^(n-1)]^T
    C_o = np.zeros((1, n))
    C_o[0, -1] = 1.0

    return A_o, B_o, C_o, np.array([[D_o]])


def observability_matrix(A, C):
    """Construct O = [C; C A; ...; C A^(n-1)]."""
    A = np.array(A, dtype=float)
    C = np.array(C, dtype=float)
    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 transfer_value(A, B, C, D, s):
    """Evaluate G(s) = C (sI - A)^(-1) B + D at a scalar s."""
    A = np.array(A, dtype=float)
    n = A.shape[0]
    return C @ np.linalg.solve(s * np.eye(n) - A, B) + D


def similarity_to_ocf(A, B, C):
    """
    For an observable SISO realization, compute the similarity map P satisfying

        A_ocf = P^(-1) A P, B_ocf = P^(-1) B, C_ocf = C P.

    The formula uses observability matrices:

        O_ocf = O_original P  =>  P = O_original^(-1) O_ocf.
    """
    A = np.array(A, dtype=float)
    B = np.array(B, dtype=float).reshape((-1, 1))
    C = np.array(C, dtype=float).reshape((1, -1))

    # Get characteristic polynomial and use a zero numerator only to build the
    # target OCF A and C. B is not needed for the transformation matrix.
    den = np.poly(A)  # [1, a1, ..., an]
    A_target, _, C_target, _ = observable_canonical_form(den, [0.0])

    O_original = observability_matrix(A, C)
    O_target = observability_matrix(A_target, C_target)

    if np.linalg.matrix_rank(O_original) < A.shape[0]:
        raise ValueError("The original realization is not observable.")

    P = np.linalg.solve(O_original, O_target)
    A_ocf = np.linalg.solve(P, A @ P)
    B_ocf = np.linalg.solve(P, B)
    C_ocf = C @ P
    return P, A_ocf, B_ocf, C_ocf


def main():
    # Example:
    # G(s) = (2 s^2 + 5 s + 3) / (s^3 + 4 s^2 + 6 s + 8)
    den = [1.0, 4.0, 6.0, 8.0]
    num = [2.0, 5.0, 3.0]

    A_o, B_o, C_o, D_o = observable_canonical_form(den, num)

    print("A_o =\n", A_o)
    print("B_o =\n", B_o)
    print("C_o =\n", C_o)
    print("D_o =\n", D_o)

    # Check at selected complex frequencies.
    for s in [1.0, 2.0, 1.0 + 1.0j]:
        g_state = transfer_value(A_o, B_o, C_o, D_o, s)[0, 0]
        g_poly = np.polyval(num, s) / np.polyval(den, s)
        print(f"s={s:>8}: state={g_state}, polynomial={g_poly}, error={abs(g_state-g_poly):.2e}")

    O_o = observability_matrix(A_o, C_o)
    print("rank(O_o) =", np.linalg.matrix_rank(O_o))


if __name__ == "__main__":
    main()

      

Chapter17_Lesson1.cpp

/*
Chapter17_Lesson1.cpp

Observable Canonical Form (OCF) construction for a SISO transfer function.

Compile:
    g++ -std=c++17 Chapter17_Lesson1.cpp -o Chapter17_Lesson1

This version uses only the C++ standard library. For larger systems, prefer
Eigen or Armadillo for numerical linear algebra.
*/

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

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

struct StateSpace {
    Matrix A;
    Vector B;
    Vector C;
    double D;
};

void print_matrix(const std::string& name, const Matrix& M) {
    std::cout << name << " =\n";
    for (const auto& row : M) {
        for (double v : row) {
            std::cout << std::setw(12) << v << " ";
        }
        std::cout << "\n";
    }
}

void print_vector(const std::string& name, const Vector& v) {
    std::cout << name << " = [ ";
    for (double x : v) {
        std::cout << x << " ";
    }
    std::cout << "]\n";
}

Matrix multiply(const Matrix& A, const Matrix& B) {
    int n = static_cast<int>(A.size());
    int m = static_cast<int>(B[0].size());
    int p = static_cast<int>(B.size());
    Matrix C(n, Vector(m, 0.0));
    for (int i = 0; i < n; ++i) {
        for (int k = 0; k < p; ++k) {
            for (int j = 0; j < m; ++j) {
                C[i][j] += A[i][k] * B[k][j];
            }
        }
    }
    return C;
}

Matrix identity(int n) {
    Matrix I(n, Vector(n, 0.0));
    for (int i = 0; i < n; ++i) {
        I[i][i] = 1.0;
    }
    return I;
}

int rank_gaussian(Matrix M, double tol = 1e-10) {
    int rows = static_cast<int>(M.size());
    int cols = static_cast<int>(M[0].size());
    int r = 0;
    for (int c = 0; c < cols && r < rows; ++c) {
        int pivot = r;
        for (int i = r + 1; i < rows; ++i) {
            if (std::fabs(M[i][c]) > std::fabs(M[pivot][c])) {
                pivot = i;
            }
        }
        if (std::fabs(M[pivot][c]) < tol) {
            continue;
        }
        std::swap(M[pivot], M[r]);
        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;
}

StateSpace observable_canonical_form(Vector den, Vector num) {
    if (den.empty() || std::fabs(den[0]) < 1e-14) {
        throw std::runtime_error("Denominator leading coefficient must be nonzero.");
    }

    double leading = den[0];
    for (double& v : den) {
        v /= leading;
    }
    for (double& v : num) {
        v /= leading;
    }

    int n = static_cast<int>(den.size()) - 1;
    if (n <= 0) {
        throw std::runtime_error("Denominator degree must be positive.");
    }
    if (static_cast<int>(num.size()) > n + 1) {
        throw std::runtime_error("Expected a proper transfer function.");
    }

    Vector num_pad(n + 1, 0.0);
    int offset = n + 1 - static_cast<int>(num.size());
    for (int i = 0; i < static_cast<int>(num.size()); ++i) {
        num_pad[offset + i] = num[i];
    }

    double D = num_pad[0];
    Vector remainder(n + 1, 0.0);
    for (int i = 0; i <= n; ++i) {
        remainder[i] = num_pad[i] - D * den[i];
    }

    Matrix A(n, Vector(n, 0.0));
    for (int i = 1; i < n; ++i) {
        A[i][i - 1] = 1.0;
    }
    for (int i = 0; i < n; ++i) {
        A[i][n - 1] = -den[n - i];
    }

    Vector B(n, 0.0);
    for (int i = 0; i < n; ++i) {
        B[i] = remainder[n - i];
    }

    Vector C(n, 0.0);
    C[n - 1] = 1.0;

    return {A, B, C, D};
}

Matrix observability_matrix(const Matrix& A, const Vector& C) {
    int n = static_cast<int>(A.size());
    Matrix O(n, Vector(n, 0.0));
    Matrix Ak = identity(n);

    for (int row = 0; row < n; ++row) {
        for (int j = 0; j < n; ++j) {
            double value = 0.0;
            for (int k = 0; k < n; ++k) {
                value += C[k] * Ak[k][j];
            }
            O[row][j] = value;
        }
        Ak = multiply(Ak, A);
    }
    return O;
}

int main() {
    // G(s) = (2 s^2 + 5 s + 3) / (s^3 + 4 s^2 + 6 s + 8)
    Vector den = {1.0, 4.0, 6.0, 8.0};
    Vector num = {2.0, 5.0, 3.0};

    StateSpace sys = observable_canonical_form(den, num);

    print_matrix("A_o", sys.A);
    print_vector("B_o", sys.B);
    print_vector("C_o", sys.C);
    std::cout << "D_o = " << sys.D << "\n";

    Matrix O = observability_matrix(sys.A, sys.C);
    print_matrix("O_o", O);
    std::cout << "rank(O_o) = " << rank_gaussian(O) << "\n";

    return 0;
}

      

Chapter17_Lesson1.java

/*
Chapter17_Lesson1.java

Observable Canonical Form (OCF) construction for a SISO transfer function.

Compile:
    javac Chapter17_Lesson1.java
Run:
    java Chapter17_Lesson1

This version uses plain Java arrays. For larger numerical systems, use EJML,
Apache Commons Math, or ND4J.
*/

public class Chapter17_Lesson1 {
    static class StateSpace {
        double[][] A;
        double[] B;
        double[] C;
        double D;

        StateSpace(double[][] A, double[] B, double[] C, double D) {
            this.A = A;
            this.B = B;
            this.C = C;
            this.D = D;
        }
    }

    static StateSpace observableCanonicalForm(double[] denInput, double[] numInput) {
        double[] den = denInput.clone();
        double[] num = numInput.clone();

        if (den.length < 2 || Math.abs(den[0]) < 1e-14) {
            throw new IllegalArgumentException("Denominator leading coefficient must be nonzero.");
        }

        double leading = den[0];
        for (int i = 0; i < den.length; i++) {
            den[i] /= leading;
        }
        for (int i = 0; i < num.length; i++) {
            num[i] /= leading;
        }

        int n = den.length - 1;
        if (num.length > n + 1) {
            throw new IllegalArgumentException("Expected a proper transfer function.");
        }

        double[] numPad = new double[n + 1];
        int offset = n + 1 - num.length;
        for (int i = 0; i < num.length; i++) {
            numPad[offset + i] = num[i];
        }

        double D = numPad[0];
        double[] remainder = new double[n + 1];
        for (int i = 0; i <= n; i++) {
            remainder[i] = numPad[i] - D * den[i];
        }

        double[][] A = new double[n][n];
        for (int i = 1; i < n; i++) {
            A[i][i - 1] = 1.0;
        }
        for (int i = 0; i < n; i++) {
            A[i][n - 1] = -den[n - i];
        }

        double[] B = new double[n];
        for (int i = 0; i < n; i++) {
            B[i] = remainder[n - i];
        }

        double[] C = new double[n];
        C[n - 1] = 1.0;

        return new StateSpace(A, B, C, D);
    }

    static double[][] multiply(double[][] A, double[][] B) {
        int n = A.length;
        int p = B.length;
        int m = B[0].length;
        double[][] C = new double[n][m];

        for (int i = 0; i < n; i++) {
            for (int k = 0; k < p; k++) {
                for (int j = 0; j < m; j++) {
                    C[i][j] += A[i][k] * B[k][j];
                }
            }
        }
        return C;
    }

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

    static double[][] observabilityMatrix(double[][] A, double[] C) {
        int n = A.length;
        double[][] O = new double[n][n];
        double[][] Ak = identity(n);

        for (int row = 0; row < n; row++) {
            for (int j = 0; j < n; j++) {
                double value = 0.0;
                for (int k = 0; k < n; k++) {
                    value += C[k] * Ak[k][j];
                }
                O[row][j] = value;
            }
            Ak = multiply(Ak, A);
        }
        return O;
    }

    static int rankGaussian(double[][] input, double tol) {
        int rows = input.length;
        int cols = input[0].length;
        double[][] M = new double[rows][cols];

        for (int i = 0; i < rows; i++) {
            System.arraycopy(input[i], 0, M[i], 0, cols);
        }

        int 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]) < tol) {
                continue;
            }
            double[] tmp = M[pivot];
            M[pivot] = M[r];
            M[r] = 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 void printMatrix(String name, double[][] M) {
        System.out.println(name + " =");
        for (double[] row : M) {
            for (double v : row) {
                System.out.printf("%12.6f ", v);
            }
            System.out.println();
        }
    }

    static void printVector(String name, double[] v) {
        System.out.print(name + " = [ ");
        for (double x : v) {
            System.out.printf("%10.6f ", x);
        }
        System.out.println("]");
    }

    public static void main(String[] args) {
        // G(s) = (2 s^2 + 5 s + 3) / (s^3 + 4 s^2 + 6 s + 8)
        double[] den = {1.0, 4.0, 6.0, 8.0};
        double[] num = {2.0, 5.0, 3.0};

        StateSpace sys = observableCanonicalForm(den, num);

        printMatrix("A_o", sys.A);
        printVector("B_o", sys.B);
        printVector("C_o", sys.C);
        System.out.println("D_o = " + sys.D);

        double[][] O = observabilityMatrix(sys.A, sys.C);
        printMatrix("O_o", O);
        System.out.println("rank(O_o) = " + rankGaussian(O, 1e-10));
    }
}

      

Chapter17_Lesson1.m

% Chapter17_Lesson1.m
%
% Observable Canonical Form (OCF) construction for a SISO transfer function.
% This script uses a from-scratch construction and then validates it with
% Control System Toolbox functions if they are available.

clear; clc;

% G(s) = (2 s^2 + 5 s + 3)/(s^3 + 4 s^2 + 6 s + 8)
den = [1 4 6 8];
num = [2 5 3];

[Ao, Bo, Co, Do] = observableCanonicalForm(den, num);

disp('Ao ='); disp(Ao);
disp('Bo ='); disp(Bo);
disp('Co ='); disp(Co);
disp('Do ='); disp(Do);

Oo = obsv(Ao, Co);
disp('rank(Oo) ='); disp(rank(Oo));

% Validation with Control System Toolbox.
if exist('ss', 'file') == 2 && exist('tf', 'file') == 2
    sys_ocf = ss(Ao, Bo, Co, Do);
    disp('Transfer function recovered from OCF:');
    tf(sys_ocf)
end

% Simulink use:
% 1. Insert a State-Space block.
% 2. Set A = Ao, B = Bo, C = Co, D = Do in the block parameters.
% 3. Connect a Step block to the input and a Scope block to the output.
% 4. Compare with a Transfer Fcn block using numerator num and denominator den.

function [Ao, Bo, Co, Do] = observableCanonicalForm(den, num)
    if abs(den(1)) < 1e-14
        error('Denominator leading coefficient must be nonzero.');
    end

    denLeading = den(1);
    den = den / denLeading;
    num = num / denLeading;

    n = length(den) - 1;
    if length(num) > n + 1
        error('Expected a proper transfer function.');
    end

    numPad = zeros(1, n + 1);
    numPad(end - length(num) + 1:end) = num;

    Do = numPad(1);
    remainder = numPad - Do * den;

    Ao = zeros(n, n);
    for i = 2:n
        Ao(i, i - 1) = 1;
    end
    Ao(:, n) = -flipud(den(2:end).');

    Bo = flipud(remainder(2:end).');
    Co = zeros(1, n);
    Co(n) = 1;
end

      

Chapter17_Lesson1.nb

(* Chapter17_Lesson1.nb *)

ClearAll[observableCanonicalForm, observabilityMatrix, s];

observableCanonicalForm[denInput_List, numInput_List] := Module[
  {den = N[denInput], num = N[numInput], n, numPad, d, remainder, ao, bo, co},

  If[Abs[den[[1]]] < 10^-14,
    Print["Denominator leading coefficient must be nonzero."];
    Return[$Failed];
  ];

  num = num/den[[1]];
  den = den/den[[1]];
  n = Length[den] - 1;

  If[Length[num] > n + 1,
    Print["Expected a proper transfer function."];
    Return[$Failed];
  ];

  numPad = Join[ConstantArray[0, n + 1 - Length[num]], num];
  d = numPad[[1]];
  remainder = numPad - d den;

  ao = ConstantArray[0, {n, n}];
  Do[ao[[i, i - 1]] = 1, {i, 2, n}];
  Do[ao[[i, n]] = -den[[n - i + 2]], {i, 1, n}];

  bo = Reverse[Rest[remainder]];
  co = ConstantArray[0, n];
  co[[n]] = 1;

  {ao, bo, co, d}
];

observabilityMatrix[a_, c_] := Module[{n = Length[a]},
  Table[c . MatrixPower[a, k], {k, 0, n - 1}]
];

den = {1, 4, 6, 8};
num = {2, 5, 3};

{Ao, Bo, Co, DoFeed} = observableCanonicalForm[den, num];

Print["Ao = "]; MatrixForm[Ao]
Print["Bo = "]; MatrixForm[Bo]
Print["Co = "]; MatrixForm[Co]
Print["Do = ", DoFeed]

Oo = observabilityMatrix[Ao, Co];
Print["Observability matrix rank = ", MatrixRank[Oo]];

gState = Co . Inverse[s IdentityMatrix[Length[Ao]] - Ao] . Bo + DoFeed // Simplify;
gPoly = (2 s^2 + 5 s + 3)/(s^3 + 4 s^2 + 6 s + 8);

Print["G_state(s) = ", gState // Simplify];
Print["G_poly(s)  = ", gPoly];
Print["Difference = ", Simplify[gState - gPoly]];

      

Problem 1 (Constructing OCF from a Second-Order Transfer Function): Construct an observable canonical realization for \( G(s)=\frac{3s+7}{s^2+5s+6} \).

Solution: Here \( a_1=5,\;a_2=6 \), \( b_0=3 \), and \( b_1=7 \). Hence

\[ \mathbf{A}_o= \begin{bmatrix} 0 & -6\\ 1 & -5 \end{bmatrix}, \quad \mathbf{B}_o= \begin{bmatrix} 7\\3 \end{bmatrix}, \quad \mathbf{C}_o= \begin{bmatrix} 0 & 1 \end{bmatrix}, \quad D_o=0 . \]

The observability matrix is \( \mathcal{O}_o=\begin{bmatrix}0 & 1\\1 & -5\end{bmatrix} \), whose determinant is \( -1 \). Therefore the realization is observable.

Problem 2 (Proper but Non-Strictly Proper Case): Construct the OCF data for \( G(s)=\frac{2s^2+3s+4}{s^2+5s+6} \).

Solution: Since numerator and denominator have the same degree, separate the direct term:

\[ \frac{2s^2+3s+4}{s^2+5s+6} = 2+\frac{-7s-8}{s^2+5s+6}. \]

Thus \( D_o=2 \), \( b_0=-7 \), and \( b_1=-8 \). The OCF matrices are

\[ \mathbf{A}_o= \begin{bmatrix} 0 & -6\\ 1 & -5 \end{bmatrix}, \quad \mathbf{B}_o= \begin{bmatrix} -8\\-7 \end{bmatrix}, \quad \mathbf{C}_o= \begin{bmatrix} 0 & 1 \end{bmatrix}, \quad D_o=2 . \]

Problem 3 (Observability Proof for \( n=3 \)): For \( \mathbf{A}_o=\begin{bmatrix}0&0&-a_3\\1&0&-a_2\\0&1&-a_1\end{bmatrix} \) and \( \mathbf{C}_o=[0\;0\;1] \), compute \( \mathcal{O}_o \) and prove observability.

Solution:

\[ \mathbf{C}_o= \begin{bmatrix}0 & 0 & 1\end{bmatrix},\qquad \mathbf{C}_o\mathbf{A}_o= \begin{bmatrix}0 & 1 & -a_1\end{bmatrix}, \]

\[ \mathbf{C}_o\mathbf{A}_o^2= \begin{bmatrix} 1 & -a_1 & a_1^2-a_2 \end{bmatrix}. \]

Therefore

\[ \mathcal{O}_o= \begin{bmatrix} 0 & 0 & 1\\ 0 & 1 & -a_1\\ 1 & -a_1 & a_1^2-a_2 \end{bmatrix}, \qquad \det(\mathcal{O}_o)=-1\neq 0 . \]

Hence the pair \( (\mathbf{C}_o,\mathbf{A}_o) \) is observable for every coefficient choice.

Problem 4 (Coordinate Transformation to OCF): Suppose \( (\mathbf{A},\mathbf{B},\mathbf{C},D) \) is observable and has the same transfer function as a known OCF realization. Derive the transformation matrix.

Solution: Let \( \mathcal{O}=[\mathbf{C};\mathbf{C}\mathbf{A};\cdots;\mathbf{C}\mathbf{A}^{n-1}] \) and \( \mathcal{O}_o=[\mathbf{C}_o;\mathbf{C}_o\mathbf{A}_o;\cdots;\mathbf{C}_o\mathbf{A}_o^{n-1}] \). If \( \mathbf{x}=\mathbf{P}\mathbf{z} \), then \( \mathcal{O}_o=\mathcal{O}\mathbf{P} \). Since the original realization is observable, \( \mathcal{O} \) is nonsingular, so

\[ \mathbf{P}=\mathcal{O}^{-1}\mathcal{O}_o . \]

In numerical computation, solve \( \mathcal{O}\mathbf{P}=\mathcal{O}_o \) rather than explicitly computing \( \mathcal{O}^{-1} \).

Problem 5 (Duality Check): Let \( \mathbf{A}_o=\mathbf{A}_c^T \), \( \mathbf{B}_o=\mathbf{C}_c^T \), and \( \mathbf{C}_o=\mathbf{B}_c^T \). Prove that the OCF and CCF have the same scalar transfer function.

Solution: Starting from the OCF transfer function,

\[ \mathbf{C}_o(s\mathbf{I}-\mathbf{A}_o)^{-1}\mathbf{B}_o = \mathbf{B}_c^T(s\mathbf{I}-\mathbf{A}_c^T)^{-1}\mathbf{C}_c^T . \]

Since \( (s\mathbf{I}-\mathbf{A}_c^T)^{-1}=[(s\mathbf{I}-\mathbf{A}_c)^{-1}]^T \), the expression equals the transpose of the scalar \( \mathbf{C}_c(s\mathbf{I}-\mathbf{A}_c)^{-1}\mathbf{B}_c \). A scalar equals its transpose, so both realizations produce the same transfer function.

Observable canonical form places the denominator coefficients in an observer companion matrix, the numerator coefficients in the input vector, and the output map as a simple selector of the final state. Its observability follows from an anti-triangular observability matrix with nonzero determinant. OCF is the transpose-dual of CCF and is especially useful for theoretical proofs, realization exercises, and understanding how output measurements encode internal state coordinates.