Previous lessons introduced the companion matrix, construction from transfer-function coefficients, and structural properties of controllable canonical form (CCF). Here the goal is operational: given a concrete rational transfer function, form \( \mathbf{A}_c,\mathbf{B}_c,\mathbf{C}_c,D \), verify the Kalman rank condition, and prove that the realization reproduces the original input-output map.

We restrict attention to monic strictly proper SISO transfer functions, because this is the cleanest setting for the standard CCF introduced in the previous lessons:

\[ G(s)=\frac{N(s)}{d(s)},\qquad d(s)=s^n+a_{n-1}s^{n-1}+\cdots+a_1s+a_0,\qquad \deg N(s) < n. \]

If the numerator is \( N(s)=b_{n-1}s^{n-1}+\cdots+b_1s+b_0 \), then missing high-order numerator powers are filled by zero coefficients.

The standard CCF construction maps denominator coefficients into the last row of the companion matrix and numerator coefficients into the output row:

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

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

flowchart TD
  A["Transfer function G(s) = N(s)/d(s)"] --> B["Normalize denominator so leading coefficient is 1"]
  B --> C["Read denominator: a0, a1, ..., a(n-1)"]
  C --> D["Build companion A: \nlast row is negative denominator coefficients"]
  B --> E["Read numerator: \nb0, b1, ..., b(n-1)"]
  E --> F["Build C row from numerator coefficients"]
  D --> G["Use B = last basis vector"]
  F --> H["State-space realization"]
  G --> H
  H --> I["Verify: rank controllability matrix equals n"]
  H --> J["Verify: C (sI - A)^(-1) B + D equals G(s)"]
        

Let \( \mathbf{z}(s)=(s\mathbf{I}-\mathbf{A}_c)^{-1}\mathbf{B}_c \). Solving the companion linear system gives the structured vector

\[ \mathbf{z}(s)=\frac{1}{d(s)} \begin{bmatrix} 1\\s\\s^2\\\vdots\\s^{n-1} \end{bmatrix}. \]

Therefore,

\[ \mathbf{C}_c(s\mathbf{I}-\mathbf{A}_c)^{-1}\mathbf{B}_c =\frac{b_0+b_1s+\cdots+b_{n-1}s^{n-1}}{d(s)} =\frac{N(s)}{d(s)}. \]

The controllability matrix is \( \mathcal{C}=[\mathbf{B}_c\;\mathbf{A}_c\mathbf{B}_c\;\cdots\;\mathbf{A}_c^{n-1}\mathbf{B}_c] \). For the companion pair \( (\mathbf{A}_c,\mathbf{B}_c) \), this matrix has determinant \( (-1)^{n(n-1)/2} \) up to the ordering convention used for the state basis. Hence \( \operatorname{rank}\mathcal{C}=n \) for every choice of finite denominator coefficients.

Consider the second-order transfer function

\[ G_1(s)=\frac{2}{s^2+3s+2}. \]

Here \( a_1=3 \), \( a_0=2 \), \( b_1=0 \), and \( b_0=2 \). Thus

\[ \mathbf{A}_1= \begin{bmatrix}0&1\\-2&-3\end{bmatrix},\qquad \mathbf{B}_1=\begin{bmatrix}0\\1\end{bmatrix},\qquad \mathbf{C}_1=\begin{bmatrix}2&0\end{bmatrix},\qquad D_1=0. \]

The controllability matrix is

\[ \mathcal{C}_1= \begin{bmatrix}\mathbf{B}_1&\mathbf{A}_1\mathbf{B}_1\end{bmatrix} = \begin{bmatrix}0&1\\1&-3\end{bmatrix}, \qquad \det(\mathcal{C}_1)=-1\neq 0. \]

The input-output verification is direct:

\[ (s\mathbf{I}-\mathbf{A}_1)^{-1}\mathbf{B}_1 =\frac{1}{s^2+3s+2}\begin{bmatrix}1\\s\end{bmatrix}, \qquad \mathbf{C}_1(s\mathbf{I}-\mathbf{A}_1)^{-1}\mathbf{B}_1 =\frac{2}{s^2+3s+2}. \]

Now consider

\[ G_2(s)=\frac{s+2}{s^3+6s^2+11s+6}. \]

The denominator coefficients are \( a_2=6 \), \( a_1=11 \), and \( a_0=6 \). The numerator is \( b_2=0 \), \( b_1=1 \), \( b_0=2 \). Hence

\[ \mathbf{A}_2= \begin{bmatrix} 0&1&0\\ 0&0&1\\ -6&-11&-6 \end{bmatrix},\qquad \mathbf{B}_2= \begin{bmatrix}0\\0\\1\end{bmatrix},\qquad \mathbf{C}_2= \begin{bmatrix}2&1&0\end{bmatrix}. \]

The controllability matrix becomes

\[ \mathcal{C}_2= \begin{bmatrix} 0&0&1\\ 0&1&-6\\ 1&-6&25 \end{bmatrix}, \qquad \det(\mathcal{C}_2)=-1. \]

Its rank is three, so all three state directions are reachable from the input. The transfer-function recovery gives

\[ \mathbf{C}_2(s\mathbf{I}-\mathbf{A}_2)^{-1}\mathbf{B}_2 =\begin{bmatrix}2&1&0\end{bmatrix} \frac{1}{d_2(s)}\begin{bmatrix}1\\s\\s^2\end{bmatrix} =\frac{s+2}{s^3+6s^2+11s+6}. \]

A third example keeps the denominator third order but uses a numerator of degree two:

\[ G_3(s)=\frac{0.5s^2+1.5s+1}{s^3+4s^2+5s+2}. \]

\[ \mathbf{A}_3= \begin{bmatrix} 0&1&0\\ 0&0&1\\ -2&-5&-4 \end{bmatrix},\qquad \mathbf{B}_3= \begin{bmatrix}0\\0\\1\end{bmatrix},\qquad \mathbf{C}_3= \begin{bmatrix}1&1.5&0.5\end{bmatrix}. \]

The output row is not arbitrary: it stores the numerator coefficients in ascending powers of \( s \). The high-order numerator term \( 0.5s^2 \) appears in the third component of \( \mathbf{C}_3 \).

flowchart LR
  E1["G1: denominator order 2, \nconstant numerator"] --> C1["C = [2, 0]"]
  E2["G2: denominator order 3, \nfirst-order numerator"] --> C2["C = [2, 1, 0]"]
  E3["G3: denominator order 3, \nsecond-order numerator"] --> C3["C = [1, 1.5, 0.5]"]
  C1 --> V["Same verification: \nC (sI - A)^(-1) B = N(s)/d(s)"]
  C2 --> V
  C3 --> V
        

In CCF, the state coordinates are usually not the physical variables of a plant. They are algebraic coordinates chosen so that the input enters through the last state equation and the denominator coefficients appear in the last row of the state matrix. The chain structure

\[ \dot{x}_1=x_2,\qquad \dot{x}_2=x_3,\qquad \ldots,\qquad \dot{x}_n=-a_0x_1-a_1x_2-\cdots-a_{n-1}x_n+u \]

exposes the differential equation represented by the transfer function. For the second-order example, this gives

\[ \dot{x}_1=x_2,\qquad \dot{x}_2=-2x_1-3x_2+u,\qquad y=2x_1. \]

Eliminating the states recovers \( \ddot{y}+3\dot{y}+2y=2u \). This confirms that CCF is a realization of the same input-output differential equation, but in a coordinate basis chosen for controllability analysis.

The Python implementation uses numpy only. It constructs the CCF matrices from coefficient arrays and numerically verifies both controllability and equality of the state-space and transfer-function evaluations.

Chapter16_Lesson4.py

#!/usr/bin/env python3
"""
Chapter16_Lesson4.py

Example systems in controllable canonical form (CCF) for continuous-time SISO
state-space systems.

The convention used here is

    d(s) = s^n + a_{n-1}s^{n-1} + ... + a_1 s + a_0
    n(s) = b_{n-1}s^{n-1} + ... + b_1 s + b_0

with

    A = [[0, 1, 0, ..., 0],
         [0, 0, 1, ..., 0],
         ...
         [-a0, -a1, ..., -a_{n-1}]]

    B = [0, 0, ..., 1]^T
    C = [b0, b1, ..., b_{n-1}]
    D = 0

so that G(s) = C(sI-A)^(-1)B = n(s)/d(s).
"""

from __future__ import annotations

import numpy as np


def controllable_canonical_form(den_desc: list[float], num_desc: list[float]):
    """Return A, B, C, D for a strictly proper SISO transfer function."""
    den = np.asarray(den_desc, dtype=float)
    num = np.asarray(num_desc, dtype=float)

    if den.ndim != 1 or num.ndim != 1:
        raise ValueError("den_desc and num_desc must be one-dimensional lists.")
    if abs(den[0]) < 1e-14:
        raise ValueError("Leading denominator coefficient must be nonzero.")

    num = num / den[0]
    den = den / den[0]
    n = len(den) - 1

    if len(num) > n:
        raise ValueError("This script handles strictly proper transfer functions only.")

    num_padded_desc = np.pad(num, (n - len(num), 0), mode="constant")

    A = np.zeros((n, n))
    if n > 1:
        A[:-1, 1:] = np.eye(n - 1)
    A[-1, :] = -den[1:][::-1]

    B = np.zeros((n, 1))
    B[-1, 0] = 1.0

    C = num_padded_desc[::-1].reshape(1, n)
    D = np.array([[0.0]])
    return A, B, C, D


def controllability_matrix(A: np.ndarray, B: np.ndarray) -> np.ndarray:
    """Compute [B AB ... A^(n-1)B]."""
    n = A.shape[0]
    cols = []
    Ak = np.eye(n)
    for _ in range(n):
        cols.append(Ak @ B)
        Ak = A @ Ak
    return np.hstack(cols)


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


def print_system(name: str, den_desc: list[float], num_desc: list[float]):
    A, B, C, D = controllable_canonical_form(den_desc, num_desc)
    Wc = controllability_matrix(A, B)

    print("\n" + "=" * 72)
    print(name)
    print("A =\n", A)
    print("B =\n", B)
    print("C =\n", C)
    print("D =\n", D)
    print("Controllability matrix Wc =\n", Wc)
    print("rank(Wc) =", np.linalg.matrix_rank(Wc))
    print("det(Wc)  =", np.linalg.det(Wc))

    test_points = [0.0, -0.5, 1.0 + 0.2j]
    for z in test_points:
        g_ss = eval_transfer(A, B, C, D, z)
        den_val = np.polyval(den_desc, z)
        num_val = np.polyval(num_desc, z)
        g_tf = num_val / den_val
        print(f"s={z!s:>10}: state-space={g_ss: .8f}, transfer-function={g_tf: .8f}")


if __name__ == "__main__":
    print_system("Example 1: G1(s) = 2/(s^2 + 3s + 2)",
                 den_desc=[1, 3, 2],
                 num_desc=[2])

    print_system("Example 2: G2(s) = (s + 2)/(s^3 + 6s^2 + 11s + 6)",
                 den_desc=[1, 6, 11, 6],
                 num_desc=[1, 2])

    print_system("Example 3: G3(s) = (0.5s^2 + 1.5s + 1)/(s^3 + 4s^2 + 5s + 2)",
                 den_desc=[1, 4, 5, 2],
                 num_desc=[0.5, 1.5, 1.0])

      

The C++ version uses only the standard library. The code includes elementary matrix multiplication, Gaussian-elimination rank computation, and determinant calculation.

Chapter16_Lesson4.cpp

/*
Chapter16_Lesson4.cpp

Construct example continuous-time SISO systems in controllable canonical form
and verify the rank of the controllability matrix.

Build:
    g++ -std=c++17 Chapter16_Lesson4.cpp -O2 -o Chapter16_Lesson4
*/

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

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

Matrix zeros(int r, int c) {
    return Matrix(r, std::vector<double>(c, 0.0));
}

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

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

double determinant(Matrix A) {
    const int n = static_cast<int>(A.size());
    double det = 1.0;

    for (int i = 0; i < n; ++i) {
        int pivot = i;
        for (int r = i + 1; r < n; ++r) {
            if (std::fabs(A[r][i]) > std::fabs(A[pivot][i])) pivot = r;
        }
        if (std::fabs(A[pivot][i]) < 1e-12) return 0.0;
        if (pivot != i) {
            std::swap(A[pivot], A[i]);
            det = -det;
        }
        det *= A[i][i];
        for (int r = i + 1; r < n; ++r) {
            double factor = A[r][i] / A[i][i];
            for (int c = i; c < n; ++c) {
                A[r][c] -= factor * A[i][c];
            }
        }
    }
    return det;
}

int rank(Matrix A) {
    const int rows = static_cast<int>(A.size());
    const int cols = static_cast<int>(A[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(A[i][c]) > std::fabs(A[pivot][c])) pivot = i;
        }
        if (std::fabs(A[pivot][c]) < 1e-10) continue;
        std::swap(A[pivot], A[r]);

        double piv = A[r][c];
        for (int j = c; j < cols; ++j) A[r][j] /= piv;

        for (int i = 0; i < rows; ++i) {
            if (i == r) continue;
            double factor = A[i][c];
            for (int j = c; j < cols; ++j) A[i][j] -= factor * A[r][j];
        }
        ++r;
    }
    return r;
}

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

StateSpace ccf(const Vector& den_desc, const Vector& num_desc) {
    if (den_desc.empty() || std::fabs(den_desc[0]) < 1e-14) {
        throw std::invalid_argument("Leading denominator coefficient must be nonzero.");
    }
    const int n = static_cast<int>(den_desc.size()) - 1;
    if (static_cast<int>(num_desc.size()) > n) {
        throw std::invalid_argument("Only strictly proper systems are handled.");
    }

    Vector den(den_desc.size());
    Vector num(num_desc.size());
    for (std::size_t i = 0; i < den.size(); ++i) den[i] = den_desc[i] / den_desc[0];
    for (std::size_t i = 0; i < num.size(); ++i) num[i] = num_desc[i] / den_desc[0];

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

    Matrix B = zeros(n, 1);
    B[n - 1][0] = 1.0;

    Matrix C = zeros(1, n);
    int offset = n - static_cast<int>(num.size());
    for (int k = 0; k < static_cast<int>(num.size()); ++k) {
        int descending_index = offset + k;
        int ascending_index = n - 1 - descending_index;
        C[0][ascending_index] = num[k];
    }

    Matrix D = zeros(1, 1);
    return {A, B, C, D};
}

Matrix controllabilityMatrix(const Matrix& A, const Matrix& B) {
    const int n = static_cast<int>(A.size());
    Matrix W = zeros(n, n);
    Matrix Ak = identity(n);

    for (int k = 0; k < n; ++k) {
        Matrix col = multiply(Ak, B);
        for (int i = 0; i < n; ++i) W[i][k] = col[i][0];
        Ak = multiply(A, Ak);
    }
    return W;
}

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

void runExample(const std::string& name, const Vector& den, const Vector& num) {
    StateSpace sys = ccf(den, num);
    Matrix W = controllabilityMatrix(sys.A, sys.B);

    std::cout << "\n" << std::string(72, '=') << "\n";
    std::cout << name << "\n";
    printMatrix("A", sys.A);
    printMatrix("B", sys.B);
    printMatrix("C", sys.C);
    printMatrix("Wc", W);
    std::cout << "rank(Wc) = " << rank(W) << "\n";
    std::cout << "det(Wc)  = " << determinant(W) << "\n";
}

int main() {
    runExample("Example 1: G1(s) = 2/(s^2 + 3s + 2)", {1, 3, 2}, {2});
    runExample("Example 2: G2(s) = (s + 2)/(s^3 + 6s^2 + 11s + 6)", {1, 6, 11, 6}, {1, 2});
    runExample("Example 3: G3(s) = (0.5s^2 + 1.5s + 1)/(s^3 + 4s^2 + 5s + 2)", {1, 4, 5, 2}, {0.5, 1.5, 1.0});
    return 0;
}

      

The Java version mirrors the C++ implementation with two-dimensional arrays. It is deliberately written without external dependencies so that students can inspect every matrix operation.

Chapter16_Lesson4.java

/*
Chapter16_Lesson4.java

Construct example continuous-time SISO systems in controllable canonical form
and verify the rank of the controllability matrix.

Build:
    javac Chapter16_Lesson4.java
    java Chapter16_Lesson4
*/

import java.util.Arrays;

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

    static double[][] zeros(int r, int c) {
        return new double[r][c];
    }

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

    static double[][] multiply(double[][] A, double[][] B) {
        int r = A.length;
        int m = A[0].length;
        int c = B[0].length;
        double[][] C = zeros(r, c);
        for (int i = 0; i < r; i++)
            for (int k = 0; k < m; k++)
                for (int j = 0; j < c; j++)
                    C[i][j] += A[i][k] * B[k][j];
        return C;
    }

    static double determinant(double[][] input) {
        int n = input.length;
        double[][] A = copy(input);
        double det = 1.0;

        for (int i = 0; i < n; i++) {
            int pivot = i;
            for (int r = i + 1; r < n; r++) {
                if (Math.abs(A[r][i]) > Math.abs(A[pivot][i])) pivot = r;
            }
            if (Math.abs(A[pivot][i]) < 1e-12) return 0.0;
            if (pivot != i) {
                double[] tmp = A[pivot];
                A[pivot] = A[i];
                A[i] = tmp;
                det = -det;
            }
            det *= A[i][i];

            for (int r = i + 1; r < n; r++) {
                double factor = A[r][i] / A[i][i];
                for (int c = i; c < n; c++) A[r][c] -= factor * A[i][c];
            }
        }
        return det;
    }

    static int rank(double[][] input) {
        double[][] A = copy(input);
        int rows = A.length;
        int cols = A[0].length;
        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(A[i][c]) > Math.abs(A[pivot][c])) pivot = i;
            }
            if (Math.abs(A[pivot][c]) < 1e-10) continue;

            double[] tmp = A[pivot];
            A[pivot] = A[r];
            A[r] = tmp;

            double piv = A[r][c];
            for (int j = c; j < cols; j++) A[r][j] /= piv;

            for (int i = 0; i < rows; i++) {
                if (i == r) continue;
                double factor = A[i][c];
                for (int j = c; j < cols; j++) A[i][j] -= factor * A[r][j];
            }
            r++;
        }
        return r;
    }

    static double[][] copy(double[][] M) {
        double[][] out = new double[M.length][M[0].length];
        for (int i = 0; i < M.length; i++) out[i] = Arrays.copyOf(M[i], M[i].length);
        return out;
    }

    static StateSpace ccf(double[] denDesc, double[] numDesc) {
        if (denDesc.length == 0 || Math.abs(denDesc[0]) < 1e-14) {
            throw new IllegalArgumentException("Leading denominator coefficient must be nonzero.");
        }
        int n = denDesc.length - 1;
        if (numDesc.length > n) {
            throw new IllegalArgumentException("Only strictly proper systems are handled.");
        }

        double[] den = new double[denDesc.length];
        double[] num = new double[numDesc.length];
        for (int i = 0; i < den.length; i++) den[i] = denDesc[i] / denDesc[0];
        for (int i = 0; i < num.length; i++) num[i] = numDesc[i] / denDesc[0];

        double[][] A = zeros(n, n);
        for (int i = 0; i < n - 1; i++) A[i][i + 1] = 1.0;
        for (int j = 0; j < n; j++) A[n - 1][j] = -den[n - j];

        double[][] B = zeros(n, 1);
        B[n - 1][0] = 1.0;

        double[][] C = zeros(1, n);
        int offset = n - num.length;
        for (int k = 0; k < num.length; k++) {
            int descendingIndex = offset + k;
            int ascendingIndex = n - 1 - descendingIndex;
            C[0][ascendingIndex] = num[k];
        }

        return new StateSpace(A, B, C, zeros(1, 1));
    }

    static double[][] controllabilityMatrix(double[][] A, double[][] B) {
        int n = A.length;
        double[][] W = zeros(n, n);
        double[][] Ak = identity(n);

        for (int k = 0; k < n; k++) {
            double[][] col = multiply(Ak, B);
            for (int i = 0; i < n; i++) W[i][k] = col[i][0];
            Ak = multiply(A, Ak);
        }
        return W;
    }

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

    static void runExample(String name, double[] den, double[] num) {
        StateSpace sys = ccf(den, num);
        double[][] W = controllabilityMatrix(sys.A, sys.B);

        System.out.println("\n" + "=".repeat(72));
        System.out.println(name);
        printMatrix("A", sys.A);
        printMatrix("B", sys.B);
        printMatrix("C", sys.C);
        printMatrix("Wc", W);
        System.out.println("rank(Wc) = " + rank(W));
        System.out.println("det(Wc)  = " + determinant(W));
    }

    public static void main(String[] args) {
        runExample("Example 1: G1(s) = 2/(s^2 + 3s + 2)", new double[]{1, 3, 2}, new double[]{2});
        runExample("Example 2: G2(s) = (s + 2)/(s^3 + 6s^2 + 11s + 6)", new double[]{1, 6, 11, 6}, new double[]{1, 2});
        runExample("Example 3: G3(s) = (0.5s^2 + 1.5s + 1)/(s^3 + 4s^2 + 5s + 2)", new double[]{1, 4, 5, 2}, new double[]{0.5, 1.5, 1.0});
    }
}

      

The MATLAB script computes the same matrices from scratch. If the Control System Toolbox is installed, it also forms ss and tf objects. If Simulink is available, it creates a simple state-space block diagram for the first example.

Chapter16_Lesson4.m

% Chapter16_Lesson4.m
% Example systems in controllable canonical form (CCF).
%
% Requires only base MATLAB for the from-scratch computations.
% The optional ss/tf/ctrb section requires Control System Toolbox.

clear; clc;

examples = {
    'Example 1: G1(s) = 2/(s^2 + 3s + 2)', [1 3 2], [2];
    'Example 2: G2(s) = (s + 2)/(s^3 + 6s^2 + 11s + 6)', [1 6 11 6], [1 2];
    'Example 3: G3(s) = (0.5s^2 + 1.5s + 1)/(s^3 + 4s^2 + 5s + 2)', [1 4 5 2], [0.5 1.5 1.0]
};

for k = 1:size(examples, 1)
    name = examples{k, 1};
    den = examples{k, 2};
    num = examples{k, 3};

    [A, B, C, D] = ccf_from_tf(den, num);
    Wc = controllability_matrix(A, B);

    fprintf('\n%s\n', repmat('=', 1, 72));
    disp(name);
    disp('A ='); disp(A);
    disp('B ='); disp(B);
    disp('C ='); disp(C);
    disp('D ='); disp(D);
    disp('Wc ='); disp(Wc);
    fprintf('rank(Wc) = %d\n', rank(Wc));
    fprintf('det(Wc)  = %.6g\n', det(Wc));

    if exist('ss', 'file') == 2 && exist('tf', 'file') == 2
        sys = ss(A, B, C, D);
        fprintf('Transfer function from ss(A,B,C,D):\n');
        tf(sys)
    end
end

% Optional Simulink generation: creates a simple continuous-time state-space
% block model for Example 1 if Simulink is available.
if exist('new_system', 'file') == 2
    [A1, B1, C1, D1] = ccf_from_tf([1 3 2], [2]);
    modelName = 'Chapter16_Lesson4_Simulink_CCF';
    if bdIsLoaded(modelName)
        close_system(modelName, 0);
    end
    new_system(modelName);
    open_system(modelName);
    add_block('simulink/Sources/Step', [modelName '/Step']);
    add_block('simulink/Continuous/State-Space', [modelName '/CCF_StateSpace']);
    add_block('simulink/Sinks/Scope', [modelName '/Scope']);
    set_param([modelName '/CCF_StateSpace'], ...
        'A', mat2str(A1), 'B', mat2str(B1), 'C', mat2str(C1), 'D', mat2str(D1));
    add_line(modelName, 'Step/1', 'CCF_StateSpace/1');
    add_line(modelName, 'CCF_StateSpace/1', 'Scope/1');
    set_param(modelName, 'StopTime', '10');
    save_system(modelName);
    fprintf('\nSimulink model saved as %s.slx\n', modelName);
end

function [A, B, C, D] = ccf_from_tf(den_desc, num_desc)
    if abs(den_desc(1)) < 1e-14
        error('Leading denominator coefficient must be nonzero.');
    end

    lead = den_desc(1);
    den_desc = den_desc / lead;
    num_desc = num_desc / lead;
    n = length(den_desc) - 1;

    if length(num_desc) > n
        error('Only strictly proper systems are handled.');
    end

    A = zeros(n, n);
    if n > 1
        A(1:n-1, 2:n) = eye(n - 1);
    end
    A(n, :) = -fliplr(den_desc(2:end));

    B = zeros(n, 1);
    B(n) = 1;

    padded = [zeros(1, n - length(num_desc)), num_desc];
    C = fliplr(padded);
    D = 0;
end

function Wc = controllability_matrix(A, B)
    n = size(A, 1);
    Wc = zeros(n, n);
    Ak = eye(n);
    for k = 1:n
        Wc(:, k) = Ak * B;
        Ak = A * Ak;
    end
end

      

Mathematica is useful for symbolic verification because it can simplify \( \mathbf{C}(s\mathbf{I}-\mathbf{A})^{-1}\mathbf{B} \) exactly.

Chapter16_Lesson4.nb

(* Chapter16_Lesson4.nb *)

ClearAll[CCF, ControllabilityMatrixManual, PrintSystem];

CCF[denDesc_List, numDesc_List] := Module[
  {den = N[denDesc], num = N[numDesc], n, A, B, C, D, padded},
  If[Abs[den[[1]]] < 10^-14,
    Print["Leading denominator coefficient must be nonzero."]; Abort[]
  ];
  num = num/den[[1]];
  den = den/den[[1]];
  n = Length[den] - 1;
  If[Length[num] > n,
    Print["Only strictly proper systems are handled."]; Abort[]
  ];

  A = ConstantArray[0, {n, n}];
  Do[A[[i, i + 1]] = 1, {i, 1, n - 1}];
  A[[n, All]] = -Reverse[Rest[den]];

  B = ConstantArray[0, {n, 1}];
  B[[n, 1]] = 1;

  padded = Join[ConstantArray[0, n - Length[num]], num];
  C = {Reverse[padded]};
  D = {{0}};
  {A, B, C, D}
];

ControllabilityMatrixManual[A_, B_] := Module[
  {n = Length[A]},
  Transpose[Table[Flatten[MatrixPower[A, k].B], {k, 0, n - 1}]]
];

PrintSystem[name_String, denDesc_List, numDesc_List] := Module[
  {A, B, C, D, Wc, s, Gss, Gtf},
  {A, B, C, D} = CCF[denDesc, numDesc];
  Wc = ControllabilityMatrixManual[A, B];

  Print[""];
  Print[StringRepeat["=", 72]];
  Print[name];
  Print["A = ", MatrixForm[A]];
  Print["B = ", MatrixForm[B]];
  Print["C = ", MatrixForm[C]];
  Print["Wc = ", MatrixForm[Wc]];
  Print["rank(Wc) = ", MatrixRank[Wc]];
  Print["det(Wc) = ", Det[Wc]];

  Gss = Simplify[(C.Inverse[s IdentityMatrix[Length[A]] - A].B + D)[[1, 1]]];
  Gtf = Apart[FromDigits[Reverse[numDesc], s]/FromDigits[Reverse[denDesc], s]];
  Print["G(s) from state space = ", Gss];
  Print["G(s) from coefficients  = ", Gtf];
];

PrintSystem[
  "Example 1: G1(s) = 2/(s^2 + 3s + 2)",
  {1, 3, 2},
  {2}
];

PrintSystem[
  "Example 2: G2(s) = (s + 2)/(s^3 + 6s^2 + 11s + 6)",
  {1, 6, 11, 6},
  {1, 2}
];

PrintSystem[
  "Example 3: G3(s) = (0.5s^2 + 1.5s + 1)/(s^3 + 4s^2 + 5s + 2)",
  {1, 4, 5, 2},
  {0.5, 1.5, 1.0}
];

      

Problem 1: Construct the CCF realization for \( G(s)=\frac{5}{s^2+4s+5} \), and verify the controllability rank condition.

Solution: Here \( a_1=4 \), \( a_0=5 \), \( b_1=0 \), and \( b_0=5 \).

\[ \mathbf{A}= \begin{bmatrix}0&1\\-5&-4\end{bmatrix},\qquad \mathbf{B}=\begin{bmatrix}0\\1\end{bmatrix},\qquad \mathbf{C}=\begin{bmatrix}5&0\end{bmatrix}. \]

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

Therefore the realization is controllable.

Problem 2: For \( G(s)=\frac{3s+1}{s^3+2s^2+4s+8} \), write the CCF matrices.

Solution: The denominator gives \( a_2=2,a_1=4,a_0=8 \). The numerator has \( b_2=0,b_1=3,b_0=1 \).

\[ \mathbf{A}= \begin{bmatrix} 0&1&0\\ 0&0&1\\ -8&-4&-2 \end{bmatrix},\qquad \mathbf{B}= \begin{bmatrix}0\\0\\1\end{bmatrix},\qquad \mathbf{C}= \begin{bmatrix}1&3&0\end{bmatrix}. \]

Problem 3: Prove that the output row \( \mathbf{C}=[b_0\;b_1\;\cdots\;b_{n-1}] \) reproduces the numerator polynomial.

Solution: From the companion linear-system identity,

\[ (s\mathbf{I}-\mathbf{A}_c)^{-1}\mathbf{B}_c = \frac{1}{d(s)} \begin{bmatrix}1&s&s^2&\cdots&s^{n-1}\end{bmatrix}^T. \]

Multiplication by the output row yields

\[ \mathbf{C}_c(s\mathbf{I}-\mathbf{A}_c)^{-1}\mathbf{B}_c = \frac{b_0+b_1s+\cdots+b_{n-1}s^{n-1}}{d(s)}. \]

This is exactly \( N(s)/d(s) \).

Problem 4: For the second-order realization in Example 1, eliminate the state variables and derive the input-output differential equation.

Solution: The CCF equations are

\[ \dot{x}_1=x_2,\qquad \dot{x}_2=-2x_1-3x_2+u,\qquad y=2x_1. \]

Since \( x_1=y/2 \) and \( x_2=\dot{y}/2 \), differentiating gives

\[ \frac{\ddot{y}}{2}=-2\frac{y}{2}-3\frac{\dot{y}}{2}+u. \]

Multiplying by two:

\[ \ddot{y}+3\dot{y}+2y=2u. \]

Problem 5: Explain why the CCF pair \( (\mathbf{A}_c,\mathbf{B}_c) \) is controllable even if the output row \( \mathbf{C}_c \) is zero.

Solution: Controllability depends only on \( \mathbf{A} \) and \( \mathbf{B} \), not on the output matrix. The Kalman matrix \( [\mathbf{B}_c\;\mathbf{A}_c\mathbf{B}_c\;\cdots\;\mathbf{A}_c^{n-1}\mathbf{B}_c] \) has full rank for the companion pair. Setting \( \mathbf{C}_c=\mathbf{0} \) would make the output uninformative, but it would not remove actuator authority over the internal state.

CCF gives a systematic realization for SISO transfer functions by storing denominator coefficients in the companion matrix and numerator coefficients in the output row. The companion input pair is always controllable under the standard construction, and the identity \( (s\mathbf{I}-\mathbf{A}_c)^{-1}\mathbf{B}_c =d(s)^{-1}[1\;s\;\cdots\;s^{n-1}]^T \) explains why the realization reproduces the transfer function exactly.