A state vector is not unique. If \( \mathbf{x}\in\mathbb{R}^n \) is a valid state vector and \( \mathbf{T}\in\mathbb{R}^{n\times n} \) is nonsingular, then \( \mathbf{z}=\mathbf{T}^{-1}\mathbf{x} \) is an equally valid coordinate vector. The transformed realization is

\[ \dot{\mathbf{z}} = \underbrace{\mathbf{T}^{-1}\mathbf{A}\mathbf{T}}_{\mathbf{A}_z}\mathbf{z} + \underbrace{\mathbf{T}^{-1}\mathbf{B}}_{\mathbf{B}_z}\mathbf{u},\qquad \mathbf{y} = \underbrace{\mathbf{C}\mathbf{T}}_{\mathbf{C}_z}\mathbf{z}+ \mathbf{D}\mathbf{u}. \]

The transfer function is invariant under this transformation:

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

Therefore, coordinates do not change input-output behavior, eigenvalues, controllability, or observability in exact arithmetic. They do, however, strongly affect floating-point algorithms. A poor coordinate choice can make a theoretically controllable system look nearly uncontrollable numerically.

flowchart TD
  A["Physical model with mixed units"] --> B["Choose state coordinates"]
  B --> C["Check state magnitudes and matrix norms"]
  C --> D["Apply diagonal scaling or structured transformation"]
  D --> E["Compute ranks, Gramians, poles, feedback gains"]
  E --> F["Validate in original physical units"]
  F --> G["Implement controller with documented coordinate map"]
        

Similarity preserves characteristic polynomials:

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

The controllability matrix transforms as

\[ \mathcal{C}_z = \left[\mathbf{B}_z,\mathbf{A}_z\mathbf{B}_z,\ldots, \mathbf{A}_z^{n-1}\mathbf{B}_z\right] =\mathbf{T}^{-1}\mathcal{C}. \]

The observability matrix transforms as

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

Hence exact ranks are invariant: \( \operatorname{rank}(\mathcal{C}_z)=\operatorname{rank}(\mathcal{C}) \) and \( \operatorname{rank}(\mathcal{O}_z)=\operatorname{rank}(\mathcal{O}) \). But finite-precision algorithms see singular values, not symbolic ranks. For any nonsingular \( \mathbf{T} \),

\[ \kappa_2(\mathcal{C}_z) =\kappa_2(\mathbf{T}^{-1}\mathcal{C}) \leq \kappa_2(\mathbf{T})\kappa_2(\mathcal{C}),\qquad \kappa_2(\mathcal{O}_z) \leq \kappa_2(\mathcal{O})\kappa_2(\mathbf{T}). \]

This inequality explains why arbitrary transformations should not be used casually. A badly conditioned transformation can amplify roundoff error even when it looks algebraically harmless.

In physical models, state components often have incompatible units and magnitudes: meters, radians, amperes, pressures, angular velocities, and temperatures may appear in the same vector. A common first correction is diagonal scaling. Let

\[ \mathbf{S}=\operatorname{diag}(x_{1,\mathrm{nom}},\ldots, x_{n,\mathrm{nom}}),\qquad \mathbf{z}=\mathbf{S}^{-1}\mathbf{x}. \]

The numbers \( x_{i,\mathrm{nom}} \) should be representative state magnitudes: rated operating points, expected maximum deviations, or values from a preliminary simulation. The scaled variables \( z_i \) are dimensionless and ideally of order one:

\[ |z_i|=\left|\frac{x_i}{x_{i,\mathrm{nom}}}\right| \approx 1 \quad \text{during typical operation}. \]

If inputs and outputs also have widely different magnitudes, introduce \( \mathbf{u}=\mathbf{S}_u\mathbf{v} \) and \( \mathbf{y}=\mathbf{S}_y\boldsymbol{\eta} \). Then

\[ \dot{\mathbf{z}}=\mathbf{S}^{-1}\mathbf{A}\mathbf{S}\mathbf{z} +\mathbf{S}^{-1}\mathbf{B}\mathbf{S}_u\mathbf{v},\qquad \boldsymbol{\eta}=\mathbf{S}_y^{-1}\mathbf{C}\mathbf{S}\mathbf{z} +\mathbf{S}_y^{-1}\mathbf{D}\mathbf{S}_u\mathbf{v}. \]

For controller implementation, gains must be mapped back. If the scaled controller is \( \mathbf{v}=-\mathbf{K}_z\mathbf{z} \), then the physical control law is

\[ \mathbf{u}=-\underbrace{\mathbf{S}_u\mathbf{K}_z\mathbf{S}^{-1}}_{\mathbf{K}_x} \mathbf{x}. \]

The 2-norm condition number of a full-rank matrix \( \mathbf{M} \) is

\[ \kappa_2(\mathbf{M})= \frac{\sigma_{\max}(\mathbf{M})}{\sigma_{\min}(\mathbf{M})}. \]

A linear solve \( \mathbf{M}\mathbf{w}=\mathbf{b} \) is sensitive when \( \kappa_2(\mathbf{M}) \) is large. A standard first-order perturbation estimate is

\[ \frac{\|\delta \mathbf{w}\|_2}{\|\mathbf{w}\|_2} \lesssim \kappa_2(\mathbf{M}) \left( \frac{\|\delta \mathbf{M}\|_2}{\|\mathbf{M}\|_2} +\frac{\|\delta \mathbf{b}\|_2}{\|\mathbf{b}\|_2} \right). \]

Numerical rank should be decided from singular values, not from exact determinant tests. For tolerance \( \varepsilon_{\mathrm{rank}} \), declare

\[ \operatorname{rank}_\varepsilon(\mathbf{M}) =\#\left\{i:\sigma_i(\mathbf{M}) > \varepsilon_{\mathrm{rank}}\sigma_1(\mathbf{M})\right\}. \]

In controllability testing, a very small singular value means one state direction is reachable only with very large input energy or is indistinguishable from numerical roundoff.

flowchart TD
  A["Build A, B, C in physical units"] --> B["Compute singular values of controllability and observability matrices"]
  B --> C["Smallest singular value near machine precision?"]
  C -->|yes| D["Rescale states and inputs"]
  C -->|no| E["Proceed with design"]
  D --> F["Repeat rank and condition checks"]
  F --> G["Design in scaled coordinates"]
  G --> H["Map gains and initial conditions back to physical variables"]
        

For a stable continuous-time LTI system, the controllability and observability Gramians satisfy

\[ \mathbf{A}\mathbf{W}_c+\mathbf{W}_c\mathbf{A}^T+ \mathbf{B}\mathbf{B}^T=\mathbf{0},\qquad \mathbf{A}^T\mathbf{W}_o+\mathbf{W}_o\mathbf{A}+ \mathbf{C}^T\mathbf{C}=\mathbf{0}. \]

Under \( \mathbf{x}=\mathbf{T}\mathbf{z} \), the Gramians transform as

\[ \mathbf{W}_{c,z}=\mathbf{T}^{-1}\mathbf{W}_c\mathbf{T}^{-T}, \qquad \mathbf{W}_{o,z}=\mathbf{T}^T\mathbf{W}_o\mathbf{T}. \]

Balanced coordinates choose \( \mathbf{T} \) so that

\[ \mathbf{W}_{c,z}=\mathbf{W}_{o,z} =\boldsymbol{\Sigma} =\operatorname{diag}(\sigma_1,\ldots,\sigma_n),\qquad \sigma_1\geq\sigma_2\geq\cdots\geq\sigma_n>0. \]

These \( \sigma_i \) are Hankel singular values. Large values indicate states that are simultaneously easy to control and easy to observe. Small values indicate dynamically weak input-output directions. Balanced coordinates are central in model reduction, but in this lesson their immediate value is diagnostic: they expose hidden state directions that physical coordinates may obscure.

Proposition 1: Similarity transformations preserve controllability.

Proof: Since \( \mathcal{C}_z=\mathbf{T}^{-1}\mathcal{C} \) and \( \mathbf{T}^{-1} \) is nonsingular, left multiplication cannot change column rank. Therefore \( (\mathbf{A},\mathbf{B}) \) is controllable if and only if \( (\mathbf{T}^{-1}\mathbf{A}\mathbf{T}, \mathbf{T}^{-1}\mathbf{B}) \) is controllable.

Proposition 2: Similarity transformations preserve the transfer function.

Proof: Use \( s\mathbf{I}-\mathbf{T}^{-1}\mathbf{A}\mathbf{T} =\mathbf{T}^{-1}(s\mathbf{I}-\mathbf{A})\mathbf{T} \). Inverting gives \( (s\mathbf{I}-\mathbf{A}_z)^{-1} =\mathbf{T}^{-1}(s\mathbf{I}-\mathbf{A})^{-1}\mathbf{T} \). Substitution into \( \mathbf{C}_z(s\mathbf{I}-\mathbf{A}_z)^{-1}\mathbf{B}_z+\mathbf{D} \) cancels \( \mathbf{T} \) and \( \mathbf{T}^{-1} \).

Proposition 3: If \( \mathbf{K}_z \) assigns poles for \( \mathbf{A}_z-\mathbf{B}_z\mathbf{K}_z \), then \( \mathbf{K}_x=\mathbf{K}_z\mathbf{T}^{-1} \) assigns the same poles in physical coordinates.

\[ \mathbf{A}-\mathbf{B}\mathbf{K}_x =\mathbf{T}(\mathbf{A}_z-\mathbf{B}_z\mathbf{K}_z)\mathbf{T}^{-1}. \]

The closed-loop matrices are similar; hence they have the same characteristic polynomial and eigenvalues.

Good coordinate choices are not only mathematically valid; they are interpretable, documented, and numerically moderate. Use the following rules in engineering projects:

1. Prefer physical coordinates when their magnitudes are comparable and the implementation needs direct sensor-state interpretation.
2. Use diagonal scaling when state components differ by several orders of magnitude.
3. Avoid transformations with large \( \kappa_2(\mathbf{T}) \), because they amplify roundoff in both simulation and controller gains.
4. Use SVD-based rank checks for controllability and observability. Avoid determinant-based rank decisions.
5. For pole placement, design in scaled coordinates and transform the feedback gain back to physical variables before implementation.
6. For state estimators, scale process and measurement noise consistently; otherwise Kalman-filter covariance matrices become numerically misleading.

Scaling is not a substitute for modeling. If a system has a physically weak actuator, scaling can reveal the weakness more clearly, but it cannot create true controllability.

Libraries used here include numpy, scipy.linalg, scipy.signal, and optionally python-control. The script compares physical and scaled realizations and maps the feedback gain back to physical coordinates.

Chapter30_Lesson2.py

"""
Chapter30_Lesson2.py
Coordinate selection, diagonal scaling, and numerical conditioning
for state-space systems.

Dependencies:
    pip install numpy scipy control

The script compares a badly scaled physical-coordinate realization with a
diagonally scaled realization and shows how controllability/observability
matrices, Gramians, and pole placement become numerically better behaved.
"""

import numpy as np
from numpy.linalg import cond, eigvals
from scipy.linalg import solve_continuous_lyapunov as lyap
from scipy.signal import place_poles

try:
    import control
    HAS_CONTROL = True
except Exception:
    HAS_CONTROL = False


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


def observability_matrix(A: np.ndarray, C: np.ndarray) -> np.ndarray:
    """Return stacked [C; CA; ...; CA^(n-1)]."""
    n = A.shape[0]
    blocks = [C]
    Apow = np.eye(n)
    for _ in range(1, n):
        Apow = Apow @ A
        blocks.append(C @ Apow)
    return np.vstack(blocks)


def diagonal_state_scaling(A, B, C, x_nom):
    """
    Use z = S^{-1} x with S = diag(x_nom).
    Then A_z = S^{-1} A S, B_z = S^{-1} B, C_z = C S.
    """
    S = np.diag(x_nom)
    Sinv = np.diag(1.0 / x_nom)
    return Sinv @ A @ S, Sinv @ B, C @ S, S, Sinv


def gramian_condition_numbers(A, B, C):
    """
    For stable A:
        A Wc + Wc A^T + B B^T = 0
        A^T Wo + Wo A + C^T C = 0
    """
    Wc = lyap(A, -(B @ B.T))
    Wo = lyap(A.T, -(C.T @ C))
    return cond(Wc), cond(Wo), Wc, Wo


def main():
    # A simple electromechanical model with states of very different magnitudes:
    # x1 = position [m], x2 = velocity [m/s], x3 = current [A].
    A = np.array([
        [0.0,       1.0,       0.0],
        [-2.0e3,   -5.0e1,    8.0e4],
        [0.0,      -2.0e-2,  -4.0e3]
    ], dtype=float)
    B = np.array([[0.0], [0.0], [2.0e3]], dtype=float)
    C = np.array([[1.0, 0.0, 0.0]], dtype=float)

    # Nominal magnitudes chosen from engineering units or preliminary simulations.
    x_nom = np.array([1.0e-3, 1.0e-1, 1.0e1], dtype=float)
    Az, Bz, Cz, S, Sinv = diagonal_state_scaling(A, B, C, x_nom)

    print("Eigenvalues are invariant under similarity:")
    print("eig(A)  =", np.sort_complex(eigvals(A)))
    print("eig(Az) =", np.sort_complex(eigvals(Az)))

    Mc = controllability_matrix(A, B)
    Mcz = controllability_matrix(Az, Bz)
    Mo = observability_matrix(A, C)
    Moz = observability_matrix(Az, Cz)

    print("\nCondition numbers:")
    print(f"cond(A)                  = {cond(A):.3e}")
    print(f"cond(Az)                 = {cond(Az):.3e}")
    print(f"cond(controllability)    = {cond(Mc):.3e}")
    print(f"cond(scaled controll.)   = {cond(Mcz):.3e}")
    print(f"cond(observability)      = {cond(Mo):.3e}")
    print(f"cond(scaled observ.)     = {cond(Moz):.3e}")

    # Stable-system Gramian conditioning.
    kWc, kWo, Wc, Wo = gramian_condition_numbers(A, B, C)
    kWcz, kWoz, Wcz, Woz = gramian_condition_numbers(Az, Bz, Cz)
    print(f"cond(Wc), cond(Wo)       = {kWc:.3e}, {kWo:.3e}")
    print(f"cond(Wc_z), cond(Wo_z)   = {kWcz:.3e}, {kWoz:.3e}")

    # Pole placement in scaled coordinates. Convert Kz back to physical coordinates:
    # u = -Kz z = -Kz S^{-1} x, so Kx = Kz S^{-1}.
    desired_poles = np.array([-20.0, -35.0, -1200.0])
    Kz = place_poles(Az, Bz, desired_poles, method="YT").gain_matrix
    Kx_from_scaled = Kz @ Sinv
    closed_loop_poles = eigvals(A - B @ Kx_from_scaled)

    print("\nPole placement using scaled coordinates:")
    print("Kz =", Kz)
    print("Kx =", Kx_from_scaled)
    print("eig(A - B Kx) =", np.sort_complex(closed_loop_poles))

    if HAS_CONTROL:
        sys = control.ss(A, B, C, 0.0)
        sys_z = control.ss(Az, Bz, Cz, 0.0)
        print("\npython-control DC gains:")
        print(control.dcgain(sys), control.dcgain(sys_z))
        try:
            # balred/balanced realizations require optional slycot in many installs.
            hsv = control.hsvd(sys)
            print("Hankel singular values:", hsv)
        except Exception as exc:
            print("Hankel singular values unavailable:", exc)


if __name__ == "__main__":
    main()

C++ control projects often use Eigen for dense linear algebra. The code below implements controllability and observability matrices and computes condition numbers using singular values.

Chapter30_Lesson2.cpp

/*
Chapter30_Lesson2.cpp
Coordinate scaling and conditioning for a state-space realization.

Dependencies:
    Eigen 3.4 or later
Compile example:
    g++ -std=c++17 Chapter30_Lesson2.cpp -I /path/to/eigen -O2 -o Chapter30_Lesson2
*/

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

using Eigen::MatrixXd;
using Eigen::VectorXd;

MatrixXd controllabilityMatrix(const MatrixXd& A, const MatrixXd& B) {
    const int n = static_cast<int>(A.rows());
    MatrixXd Mc(n, n * B.cols());
    MatrixXd Apow = MatrixXd::Identity(n, n);
    for (int k = 0; k < n; ++k) {
        Mc.block(0, k * B.cols(), n, B.cols()) = Apow * B;
        Apow = Apow * A;
    }
    return Mc;
}

MatrixXd observabilityMatrix(const MatrixXd& A, const MatrixXd& C) {
    const int n = static_cast<int>(A.rows());
    MatrixXd Mo(n * C.rows(), n);
    MatrixXd Apow = MatrixXd::Identity(n, n);
    for (int k = 0; k < n; ++k) {
        Mo.block(k * C.rows(), 0, C.rows(), n) = C * Apow;
        Apow = Apow * A;
    }
    return Mo;
}

double conditionNumber2(const MatrixXd& M) {
    Eigen::JacobiSVD<MatrixXd> svd(M);
    auto s = svd.singularValues();
    return s(0) / s(s.size() - 1);
}

int main() {
    MatrixXd A(3, 3);
    A << 0.0,      1.0,      0.0,
        -2.0e3,  -5.0e1,   8.0e4,
         0.0,    -2.0e-2, -4.0e3;

    MatrixXd B(3, 1);
    B << 0.0, 0.0, 2.0e3;

    MatrixXd C(1, 3);
    C << 1.0, 0.0, 0.0;

    VectorXd xNom(3);
    xNom << 1.0e-3, 1.0e-1, 1.0e1;

    MatrixXd S = xNom.asDiagonal();
    MatrixXd Sinv = xNom.cwiseInverse().asDiagonal();

    MatrixXd Az = Sinv * A * S;
    MatrixXd Bz = Sinv * B;
    MatrixXd Cz = C * S;

    Eigen::EigenSolver<MatrixXd> eigA(A);
    Eigen::EigenSolver<MatrixXd> eigAz(Az);

    std::cout << std::scientific << std::setprecision(6);
    std::cout << "Eigenvalues of A:\n" << eigA.eigenvalues() << "\n\n";
    std::cout << "Eigenvalues of Az:\n" << eigAz.eigenvalues() << "\n\n";

    MatrixXd Mc = controllabilityMatrix(A, B);
    MatrixXd Mcz = controllabilityMatrix(Az, Bz);
    MatrixXd Mo = observabilityMatrix(A, C);
    MatrixXd Moz = observabilityMatrix(Az, Cz);

    std::cout << "cond(A)                = " << conditionNumber2(A) << "\n";
    std::cout << "cond(Az)               = " << conditionNumber2(Az) << "\n";
    std::cout << "cond(Mc)               = " << conditionNumber2(Mc) << "\n";
    std::cout << "cond(Mc scaled)        = " << conditionNumber2(Mcz) << "\n";
    std::cout << "cond(Mo)               = " << conditionNumber2(Mo) << "\n";
    std::cout << "cond(Mo scaled)        = " << conditionNumber2(Moz) << "\n";

    std::cout << "\nScaled realization matrices:\n";
    std::cout << "Az =\n" << Az << "\n";
    std::cout << "Bz =\n" << Bz << "\n";
    std::cout << "Cz =\n" << Cz << "\n";

    return 0;
}

EJML provides practical dense-matrix tools for Java. The implementation mirrors the C++ workflow and is suitable for educational desktop or server-side numerical experiments.

Chapter30_Lesson2.java

/*
Chapter30_Lesson2.java
Coordinate scaling and numerical conditioning for state-space systems.

Dependencies:
    EJML dense row module, for example:
    ejml-simple, ejml-ddense, ejml-core

Compile example:
    javac -cp ".;ejml-simple.jar;ejml-ddense.jar;ejml-core.jar" Chapter30_Lesson2.java
Run example:
    java  -cp ".;ejml-simple.jar;ejml-ddense.jar;ejml-core.jar" Chapter30_Lesson2
*/

import org.ejml.simple.SimpleEVD;
import org.ejml.simple.SimpleMatrix;
import org.ejml.simple.SimpleSVD;

public class Chapter30_Lesson2 {
    static SimpleMatrix controllabilityMatrix(SimpleMatrix A, SimpleMatrix B) {
        int n = A.numRows();
        int m = B.numCols();
        SimpleMatrix Mc = new SimpleMatrix(n, n * m);
        SimpleMatrix Apow = SimpleMatrix.identity(n);
        for (int k = 0; k < n; k++) {
            Mc.insertIntoThis(0, k * m, Apow.mult(B));
            Apow = Apow.mult(A);
        }
        return Mc;
    }

    static SimpleMatrix observabilityMatrix(SimpleMatrix A, SimpleMatrix C) {
        int n = A.numRows();
        int p = C.numRows();
        SimpleMatrix Mo = new SimpleMatrix(n * p, n);
        SimpleMatrix Apow = SimpleMatrix.identity(n);
        for (int k = 0; k < n; k++) {
            Mo.insertIntoThis(k * p, 0, C.mult(Apow));
            Apow = Apow.mult(A);
        }
        return Mo;
    }

    static double conditionNumber2(SimpleMatrix M) {
        SimpleSVD<SimpleMatrix> svd = M.svd();
        double[] s = svd.getSingularValues();
        double smax = 0.0;
        double smin = Double.POSITIVE_INFINITY;
        for (double value : s) {
            smax = Math.max(smax, value);
            smin = Math.min(smin, value);
        }
        return smax / smin;
    }

    static SimpleMatrix diagonal(double[] values) {
        SimpleMatrix D = new SimpleMatrix(values.length, values.length);
        for (int i = 0; i < values.length; i++) {
            D.set(i, i, values[i]);
        }
        return D;
    }

    public static void main(String[] args) {
        SimpleMatrix A = new SimpleMatrix(new double[][] {
            {0.0,       1.0,       0.0},
            {-2.0e3,   -5.0e1,    8.0e4},
            {0.0,      -2.0e-2,  -4.0e3}
        });

        SimpleMatrix B = new SimpleMatrix(new double[][] {
            {0.0},
            {0.0},
            {2.0e3}
        });

        SimpleMatrix C = new SimpleMatrix(new double[][] {
            {1.0, 0.0, 0.0}
        });

        double[] xNom = {1.0e-3, 1.0e-1, 1.0e1};
        double[] xNomInv = {1.0 / xNom[0], 1.0 / xNom[1], 1.0 / xNom[2]};

        SimpleMatrix S = diagonal(xNom);
        SimpleMatrix Sinv = diagonal(xNomInv);

        SimpleMatrix Az = Sinv.mult(A).mult(S);
        SimpleMatrix Bz = Sinv.mult(B);
        SimpleMatrix Cz = C.mult(S);

        SimpleEVD<SimpleMatrix> evdA = A.eig();
        SimpleEVD<SimpleMatrix> evdAz = Az.eig();

        System.out.println("Eigenvalues of A:");
        for (int i = 0; i < evdA.getNumberOfEigenvalues(); i++) {
            System.out.println(evdA.getEigenvalue(i));
        }

        System.out.println("\nEigenvalues of Az:");
        for (int i = 0; i < evdAz.getNumberOfEigenvalues(); i++) {
            System.out.println(evdAz.getEigenvalue(i));
        }

        SimpleMatrix Mc = controllabilityMatrix(A, B);
        SimpleMatrix Mcz = controllabilityMatrix(Az, Bz);
        SimpleMatrix Mo = observabilityMatrix(A, C);
        SimpleMatrix Moz = observabilityMatrix(Az, Cz);

        System.out.printf("%ncond(A)             = %.6e%n", conditionNumber2(A));
        System.out.printf("cond(Az)            = %.6e%n", conditionNumber2(Az));
        System.out.printf("cond(Mc)            = %.6e%n", conditionNumber2(Mc));
        System.out.printf("cond(Mc scaled)     = %.6e%n", conditionNumber2(Mcz));
        System.out.printf("cond(Mo)            = %.6e%n", conditionNumber2(Mo));
        System.out.printf("cond(Mo scaled)     = %.6e%n", conditionNumber2(Moz));

        System.out.println("\nAz:");
        Az.print();
        System.out.println("Bz:");
        Bz.print();
        System.out.println("Cz:");
        Cz.print();
    }
}

MATLAB Control System Toolbox directly supports state-space models, controllability matrices, observability matrices, Gramians, pole-placement, matrix balancing, and balanced realizations. In Simulink, use the scaled matrices in a State-Space block if the simulated state is \( \mathbf{z} \), and add gain blocks \( \mathbf{S} \) or \( \mathbf{S}^{-1} \) at interfaces where physical variables are required.

Chapter30_Lesson2.m

% Chapter30_Lesson2.m
% Coordinate selection, diagonal scaling, conditioning, and Control System
% Toolbox workflows for state-space models.

clear; clc;

A = [0        1        0;
    -2.0e3   -5.0e1   8.0e4;
     0       -2.0e-2 -4.0e3];

B = [0; 0; 2.0e3];
C = [1 0 0];
D = 0;

% Engineering nominal magnitudes for x = [position; velocity; current].
x_nom = [1.0e-3; 1.0e-1; 1.0e1];

% z = inv(S) x, with S = diag(x_nom).
S = diag(x_nom);
Sinv = inv(S);

Az = Sinv*A*S;
Bz = Sinv*B;
Cz = C*S;

fprintf('Eigenvalues of A:\n');
disp(eig(A));
fprintf('Eigenvalues of Az:\n');
disp(eig(Az));

Mc  = ctrb(A, B);
Mcz = ctrb(Az, Bz);
Mo  = obsv(A, C);
Moz = obsv(Az, Cz);

fprintf('cond(A)             = %.6e\n', cond(A));
fprintf('cond(Az)            = %.6e\n', cond(Az));
fprintf('cond(ctrb(A,B))     = %.6e\n', cond(Mc));
fprintf('cond(ctrb(Az,Bz))   = %.6e\n', cond(Mcz));
fprintf('cond(obsv(A,C))     = %.6e\n', cond(Mo));
fprintf('cond(obsv(Az,Cz))   = %.6e\n', cond(Moz));

% Stable Gramians.
Wc  = gram(A, B, 'c');
Wo  = gram(A, C, 'o');
Wcz = gram(Az, Bz, 'c');
Woz = gram(Az, Cz, 'o');

fprintf('cond(Wc),  cond(Wo)  = %.6e, %.6e\n', cond(Wc),  cond(Wo));
fprintf('cond(Wcz), cond(Woz) = %.6e, %.6e\n', cond(Wcz), cond(Woz));

% Pole placement in scaled coordinates, then transform the gain back.
p = [-20 -35 -1200];
Kz = place(Az, Bz, p);
Kx = Kz*Sinv;

fprintf('Kz =\n');
disp(Kz);
fprintf('Kx = Kz*inv(S) =\n');
disp(Kx);
fprintf('eig(A - B*Kx) =\n');
disp(eig(A - B*Kx));

% MATLAB balance acts on a matrix, not directly on a state-space realization.
% Use it as a diagnostic for eigenvalue computation, not blindly as a physical
% coordinate transformation for controller implementation.
[Abal, Tbal] = balance(A);
fprintf('cond(Abal) from balance(A) = %.6e\n', cond(Abal));

% Optional balanced-realization workflow if available:
sys = ss(A, B, C, D);
try
    [sysb, g] = balreal(sys);
    fprintf('Hankel singular values from balreal:\n');
    disp(g);
catch ME
    fprintf('balreal unavailable or not applicable: %s\n', ME.message);
end

% Simulink note:
% Use Az, Bz, Cz, D in a State-Space block when the simulated states are z.
% To display physical states, add a Gain block S after the state output.

Mathematica is useful for symbolic verification, exact eigenvalue experiments, and rapid exploration of state-space transformations.

Chapter30_Lesson2.nb

(* Chapter30_Lesson2.nb
   Wolfram Mathematica code cells for coordinate scaling and conditioning.
   Paste these cells into a notebook if a .nb front end is desired. *)

ClearAll["Global`*"];

A = {{0., 1., 0.},
     {-2.0*^3, -5.0*^1, 8.0*^4},
     {0., -2.0*^-2, -4.0*^3}};

B = {{0.}, {0.}, {2.0*^3}};
Cmat = {{1., 0., 0.}};
Dmat = {{0.}};

xNom = {1.0*^-3, 1.0*^-1, 1.0*^1};
S = DiagonalMatrix[xNom];
Sinv = Inverse[S];

Az = Sinv.A.S;
Bz = Sinv.B;
Cz = Cmat.S;

controllabilityMatrix[AA_, BB_] := Module[
  {n = Length[AA]},
  ArrayFlatten[{Table[MatrixPower[AA, k].BB, {k, 0, n - 1}]}]
];

observabilityMatrix[AA_, CC_] := Module[
  {n = Length[AA]},
  Join @@ Table[CC.MatrixPower[AA, k], {k, 0, n - 1}]
];

conditionNumber2[M_] := Max[SingularValueList[M]]/Min[SingularValueList[M]];

Print["Eigenvalues of A: ", Eigenvalues[A]];
Print["Eigenvalues of Az: ", Eigenvalues[Az]];

Mc = controllabilityMatrix[A, B];
Mcz = controllabilityMatrix[Az, Bz];
Mo = observabilityMatrix[A, Cmat];
Moz = observabilityMatrix[Az, Cz];

Print["cond(A) = ", conditionNumber2[A]];
Print["cond(Az) = ", conditionNumber2[Az]];
Print["cond(Mc) = ", conditionNumber2[Mc]];
Print["cond(Mc scaled) = ", conditionNumber2[Mcz]];
Print["cond(Mo) = ", conditionNumber2[Mo]];
Print["cond(Mo scaled) = ", conditionNumber2[Moz]];

sys = StateSpaceModel[{A, B, Cmat, Dmat}];
sysz = StateSpaceModel[{Az, Bz, Cz, Dmat}];

Print["Transfer functions are equivalent:"];
Print[TransferFunctionModel[sys]];
Print[TransferFunctionModel[sysz]];

(* Pole placement through scaled coordinates. *)
desiredPoles = {-20, -35, -1200};
Kz = StateFeedbackGains[StateSpaceModel[{Az, Bz, IdentityMatrix[3], 0}], desiredPoles];
Kx = Kz.Sinv;

Print["Kz = ", Kz];
Print["Kx = ", Kx];
Print["Closed-loop eigenvalues = ", Eigenvalues[A - B.Kx]];

Problem 1 (Similarity and Transfer Functions): Let \( \mathbf{z}=\mathbf{T}^{-1}\mathbf{x} \). Prove that the transfer function of the transformed realization equals the original transfer function.

Solution: The transformed matrices are \( \mathbf{A}_z=\mathbf{T}^{-1}\mathbf{A}\mathbf{T} \), \( \mathbf{B}_z=\mathbf{T}^{-1}\mathbf{B} \), and \( \mathbf{C}_z=\mathbf{C}\mathbf{T} \). Since \( s\mathbf{I}-\mathbf{A}_z =\mathbf{T}^{-1}(s\mathbf{I}-\mathbf{A})\mathbf{T} \), its inverse is \( \mathbf{T}^{-1}(s\mathbf{I}-\mathbf{A})^{-1}\mathbf{T} \). Therefore all transformation matrices cancel in \( \mathbf{C}_z(s\mathbf{I}-\mathbf{A}_z)^{-1}\mathbf{B}_z+\mathbf{D} \).

Problem 2 (Gain Mapping): Suppose \( \mathbf{z}=\mathbf{S}^{-1}\mathbf{x} \) and the scaled controller is \( u=-\mathbf{K}_z\mathbf{z} \). Find the controller gain in physical coordinates.

Solution: Substitute \( \mathbf{z}=\mathbf{S}^{-1}\mathbf{x} \):

\[ u=-\mathbf{K}_z\mathbf{S}^{-1}\mathbf{x} =-\mathbf{K}_x\mathbf{x},\qquad \mathbf{K}_x=\mathbf{K}_z\mathbf{S}^{-1}. \]

Problem 3 (Condition Number Bound): Show that \( \kappa_2(\mathbf{T}^{-1}\mathcal{C}) \leq \kappa_2(\mathbf{T})\kappa_2(\mathcal{C}) \) when the matrices have full column rank and compatible dimensions.

Solution: Using submultiplicativity of the 2-norm, \( \|\mathbf{T}^{-1}\mathcal{C}\|_2 \leq \|\mathbf{T}^{-1}\|_2\|\mathcal{C}\|_2 \). Also, \( \sigma_{\min}(\mathbf{T}^{-1}\mathcal{C}) \geq \sigma_{\min}(\mathbf{T}^{-1})\sigma_{\min}(\mathcal{C}) \). Combining these inequalities gives the stated bound.

Problem 4 (Numerical Rank): A controllability matrix has singular values \( 10^4, 10^1, 10^{-9} \). With tolerance \( \varepsilon_{\mathrm{rank}}=10^{-10} \), what numerical rank is reported?

Solution: The threshold is \( 10^{-10}\sigma_1=10^{-6} \). Since \( 10^{-9}<10^{-6} \), only the first two singular values are counted. The numerical rank is 2, although exact symbolic rank might be 3.

Problem 5 (Scaled Input and Output): Let \( \mathbf{x}=\mathbf{S}_x\mathbf{z} \), \( \mathbf{u}=\mathbf{S}_u\mathbf{v} \), and \( \mathbf{y}=\mathbf{S}_y\boldsymbol{\eta} \). Derive the scaled realization.

Solution: Substitute into \( \dot{\mathbf{x}}=\mathbf{A}\mathbf{x}+\mathbf{B}\mathbf{u} \). Since \( \dot{\mathbf{x}}=\mathbf{S}_x\dot{\mathbf{z}} \),

\[ \dot{\mathbf{z}}=\mathbf{S}_x^{-1}\mathbf{A}\mathbf{S}_x\mathbf{z} +\mathbf{S}_x^{-1}\mathbf{B}\mathbf{S}_u\mathbf{v},\qquad \boldsymbol{\eta}=\mathbf{S}_y^{-1}\mathbf{C}\mathbf{S}_x\mathbf{z} +\mathbf{S}_y^{-1}\mathbf{D}\mathbf{S}_u\mathbf{v}. \]

Coordinate transformations preserve exact state-space structure but can radically change numerical reliability. Diagonal scaling based on engineering units is often the most defensible first step. SVD-based rank tests, Gramian conditioning, and careful gain mapping are essential in practical state-feedback, observer, and simulation workflows.