In the previous lessons, the controllability Gramian was used to answer a precise energy question: how difficult is it to reach a particular state? For a continuous-time LTI system

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

assume in this lesson that \( \mathbf{A} \) is Hurwitz, meaning \( \operatorname{Re}(\lambda_i(\mathbf{A}))<0 \) for every eigenvalue. Then the infinite-horizon controllability Gramian is

\[ \mathbf{W}_c=\int_0^\infty e^{\mathbf{A}t}\mathbf{B}\mathbf{B}^T e^{\mathbf{A}^Tt}\,dt, \]

and it solves the continuous Lyapunov equation

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

If \( \mathbf{W}_c \) is positive definite, the minimum input energy required to drive the state from \( \mathbf{0} \) at \( -\infty \) to \( \mathbf{x}_f \) at \( 0 \) is

\[ \mathcal{E}_{u,\min}(\mathbf{x}_f)= \mathbf{x}_f^T\mathbf{W}_c^{-1}\mathbf{x}_f. \]

Balanced realizations start from this exact energy idea but ask a second question: after a state is created, how strongly does it appear in the output? The answer uses an output-energy matrix introduced only as a preview in this lesson.

A Gramian is not merely a property of the transfer function; it is also expressed in a chosen state coordinate basis. Let \( \mathbf{x}=\mathbf{T}\mathbf{z} \), where \( \mathbf{T} \) is nonsingular. Then

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

The controllability Gramian in the new coordinates is

\[ \bar{\mathbf{W} }_c= \mathbf{T}^{-1}\mathbf{W}_c\mathbf{T}^{-T}. \]

Therefore, changing coordinates can make the energy ellipsoid look rounder, more elongated, diagonal, or coupled, even though the physical input-output behavior has not changed.

\[ \mathbf{x}^T\mathbf{W}_c^{-1}\mathbf{x} = \mathbf{z}^T\bar{\mathbf{W} }_c^{-1}\mathbf{z}, \quad \mathbf{x}=\mathbf{T}\mathbf{z}. \]

This identity shows that the minimum energy is invariant, but the matrix representation of that energy is coordinate-dependent.

flowchart TD
  A["Original state x"] --> B["Controllability Gramian Wc"]
  B --> C["Energy ellipsoid: xT inv(Wc) x = constant"]
  A --> D["Change coordinates: x = T z"]
  D --> E["New Gramian: inv(T) Wc inv(T)^T"]
  E --> F["Same energy, different coordinate picture"]
  F --> G["Balanced coordinates seek a useful picture"]
        

Balanced coordinates require not only an input-energy measure but also an output-energy measure. Full observability theory will be developed in Chapters 13 through 15. For now, consider the zero-input response from an initial state \( \mathbf{x}(0)=\mathbf{x}_0 \). With \( \mathbf{u}(t)=\mathbf{0} \) and \( \mathbf{D} \) irrelevant to the zero-input output,

\[ \mathbf{y}(t)=\mathbf{C}e^{\mathbf{A}t}\mathbf{x}_0. \]

The total output energy generated by the initial condition is

\[ \mathcal{E}_y(\mathbf{x}_0)= \int_0^\infty \mathbf{y}(t)^T\mathbf{y}(t)\,dt = \mathbf{x}_0^T\mathbf{W}_o\mathbf{x}_0, \]

where

\[ \mathbf{W}_o= \int_0^\infty e^{\mathbf{A}^Tt}\mathbf{C}^T\mathbf{C} e^{\mathbf{A}t}\,dt. \]

The matrix \( \mathbf{W}_o \) solves

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

Under the same coordinate change \( \mathbf{x}=\mathbf{T}\mathbf{z} \), the output-energy Gramian transforms oppositely to the controllability Gramian:

\[ \bar{\mathbf{W} }_o=\mathbf{T}^T\mathbf{W}_o\mathbf{T}. \]

This opposite transformation law is the reason it is nontrivial to make both Gramians simultaneously simple.

A stable realization is called balanced when the controllability and output-energy Gramians are equal and diagonal in the same coordinates:

\[ \bar{\mathbf{W} }_c=\bar{\mathbf{W} }_o= \boldsymbol{\Sigma}= \operatorname{diag}(\sigma_1,\sigma_2,\ldots,\sigma_n), \quad \sigma_1\ge \sigma_2\ge \cdots \ge \sigma_n >0. \]

The positive diagonal entries \( \sigma_i \) are called the Hankel singular values. In balanced coordinates, the same scalar \( \sigma_i \) measures two things for the \( i \)-th coordinate direction:

\[ \text{large }\sigma_i \quad \Longleftrightarrow \quad \text{easy to reach and strongly visible at the output}. \]

This lesson stops here conceptually. A later model-reduction lesson would use small \( \sigma_i \) to decide which balanced coordinates can be truncated with controlled error. Here, no coordinate is removed.

The Hankel singular values can be computed from the eigenvalues of \( \mathbf{W}_c\mathbf{W}_o \):

\[ \sigma_i=\sqrt{\lambda_i(\mathbf{W}_c\mathbf{W}_o)}, \quad i=1,\ldots,n. \]

Although \( \mathbf{W}_c \) and \( \mathbf{W}_o \) individually change with the coordinate basis, the eigenvalues of their product are invariant under similarity transformations.

Proof. From Section 2 and Section 3,

\[ \bar{\mathbf{W} }_c\bar{\mathbf{W} }_o = (\mathbf{T}^{-1}\mathbf{W}_c\mathbf{T}^{-T}) (\mathbf{T}^T\mathbf{W}_o\mathbf{T}) = \mathbf{T}^{-1}\mathbf{W}_c\mathbf{W}_o\mathbf{T}. \]

Thus \( \bar{\mathbf{W} }_c\bar{\mathbf{W} }_o \) is similar to \( \mathbf{W}_c\mathbf{W}_o \), so they have the same eigenvalues. Therefore the \( \sigma_i \) are coordinate-invariant.

A constructive balancing transformation can be obtained as follows. Let \( \mathbf{W}_c=\mathbf{R}\mathbf{R}^T \) be a Cholesky factorization and diagonalize the symmetric positive matrix

\[ \mathbf{R}^T\mathbf{W}_o\mathbf{R} = \mathbf{U}\boldsymbol{\Sigma}^2\mathbf{U}^T. \]

Then one valid balancing transformation for \( \mathbf{x}=\mathbf{T}\mathbf{z} \) is

\[ \mathbf{T}=\mathbf{R}\mathbf{U}\boldsymbol{\Sigma}^{-1/2}. \]

Substitution gives

\[ \mathbf{T}^{-1}\mathbf{W}_c\mathbf{T}^{-T} = \boldsymbol{\Sigma}, \qquad \mathbf{T}^T\mathbf{W}_o\mathbf{T} = \boldsymbol{\Sigma}. \]

flowchart TD
  A["Stable minimal state-space realization"] --> B["Solve input-energy Lyapunov equation"]
  A --> C["Solve output-energy Lyapunov equation"]
  B --> D["Wc"]
  C --> E["Wo"]
  D --> F["Eigenvalues of Wc Wo"]
  E --> F
  F --> G["Hankel singular values"]
  G --> H["Build coordinate transform T"]
  H --> I["Balanced coordinates: Wc and Wo become same diagonal matrix"]
  I --> J["Stop here: no state removal in this lesson"]
        

The controllability Gramian alone identifies directions that are easy or difficult to reach. The output-energy Gramian alone identifies directions that are highly visible or weakly visible at the output. Balanced coordinates align these two viewpoints.

In balanced coordinates, the minimum energy needed to reach a coordinate vector \( \mathbf{z}=\alpha\mathbf{e}_i \) is

\[ \mathcal{E}_{u,\min}(\alpha\mathbf{e}_i) = \frac{\alpha^2}{\sigma_i}. \]

The zero-input output energy generated by the same coordinate direction is

\[ \mathcal{E}_y(\alpha\mathbf{e}_i) = \alpha^2\sigma_i. \]

Thus small \( \sigma_i \) has a double meaning: the coordinate is expensive to create and weakly visible at the output. This dual interpretation is the conceptual foundation of balanced model reduction, but reduction is deliberately postponed.

The following script computes both Gramians, Hankel singular values, a balancing transformation, and verifies that the transformed Gramians are approximately equal and diagonal.

Chapter12_Lesson5.py

# Chapter12_Lesson5.py
# Balanced-realization preview for a continuous-time stable LTI system.
# Requires: numpy, scipy

import numpy as np
from scipy.linalg import solve_continuous_lyapunov, cholesky, eigh

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

# x_dot = A x + B u, y = C x
# A is Hurwitz, so infinite-horizon Gramians are finite.
A = np.array([[-1.0, 0.3],
              [ 0.0,-2.0]])
B = np.array([[1.0],
              [0.5]])
C = np.array([[1.0, -0.2]])

# Controllability Gramian:
# A Wc + Wc A^T + B B^T = 0
Wc = solve_continuous_lyapunov(A, -(B @ B.T))

# Output-energy / observability Gramian preview:
# A^T Wo + Wo A + C^T C = 0
Wo = solve_continuous_lyapunov(A.T, -(C.T @ C))

# Symmetrize against floating-point roundoff.
Wc = 0.5 * (Wc + Wc.T)
Wo = 0.5 * (Wo + Wo.T)

# Hankel singular values are sqrt(eig(Wc Wo)).
# A numerically stable computation uses Cholesky(Wc) and eig(R^T Wo R).
R = cholesky(Wc, lower=True)
M = R.T @ Wo @ R
eigvals, U = eigh(M)
idx = np.argsort(eigvals)[::-1]
eigvals = eigvals[idx]
U = U[:, idx]
sigma = np.sqrt(np.maximum(eigvals, 0.0))
Sigma = np.diag(sigma)

# Balancing transformation for x = T z.
# In z-coordinates:
# Wc_z = T^{-1} Wc T^{-T}, Wo_z = T^T Wo T.
T = R @ U @ np.diag(1.0 / np.sqrt(sigma))
Tinv = np.linalg.inv(T)

Ab = Tinv @ A @ T
Bb = Tinv @ B
Cb = C @ T
Wc_bal = Tinv @ Wc @ Tinv.T
Wo_bal = T.T @ Wo @ T

print("A =\n", A)
print("B =\n", B)
print("C =\n", C)
print("\nControllability Gramian Wc =\n", Wc)
print("\nOutput-energy Gramian Wo =\n", Wo)
print("\nHankel singular values sigma =\n", sigma)
print("\nBalancing transformation T for x = T z =\n", T)
print("\nBalanced A =\n", Ab)
print("Balanced B =\n", Bb)
print("Balanced C =\n", Cb)
print("\nWc in balanced coordinates =\n", Wc_bal)
print("\nWo in balanced coordinates =\n", Wo_bal)
print("\nTarget diagonal Sigma =\n", Sigma)

# Energy check for one final state.
x_final = np.array([[1.0], [0.2]])
z_final = Tinv @ x_final
energy_original = float(x_final.T @ np.linalg.inv(Wc) @ x_final)
energy_balanced = float(z_final.T @ np.linalg.inv(Wc_bal) @ z_final)
print("\nMinimum input energy to reach x_final:", energy_original)
print("Same energy computed in balanced coordinates:", energy_balanced)

      

This C++ example solves the two-dimensional Lyapunov equations from scratch by vectorizing the matrix equation into a four-variable linear system. It then computes the Hankel singular values from \( \mathbf{W}_c\mathbf{W}_o \).

Chapter12_Lesson5.cpp

// Chapter12_Lesson5.cpp
// Gramian and Hankel-singular-value preview without external libraries.
// Compile: g++ -std=c++17 Chapter12_Lesson5.cpp -o Chapter12_Lesson5

#include <array>
#include <cmath>
#include <iomanip>
#include <iostream>
#include <stdexcept>

using Mat2 = std::array<std::array<double, 2>, 2>;
using Vec4 = std::array<double, 4>;
using Mat4 = std::array<std::array<double, 4>, 4>;

Mat2 add(const Mat2& X, const Mat2& Y) {
    Mat2 Z{};
    for (int i = 0; i < 2; ++i)
        for (int j = 0; j < 2; ++j)
            Z[i][j] = X[i][j] + Y[i][j];
    return Z;
}

Mat2 mul(const Mat2& X, const Mat2& Y) {
    Mat2 Z{};
    for (int i = 0; i < 2; ++i)
        for (int j = 0; j < 2; ++j)
            for (int k = 0; k < 2; ++k)
                Z[i][j] += X[i][k] * Y[k][j];
    return Z;
}

Mat2 transpose(const Mat2& X) {
    return { { {X[0][0], X[1][0]}, {X[0][1], X[1][1]} } };
}

void printMat(const char* name, const Mat2& X) {
    std::cout << name << "\n";
    for (int i = 0; i < 2; ++i) {
        std::cout << "  ";
        for (int j = 0; j < 2; ++j)
            std::cout << std::setw(12) << X[i][j] << " ";
        std::cout << "\n";
    }
}

Vec4 solve4(Mat4 A, Vec4 b) {
    for (int k = 0; k < 4; ++k) {
        int pivot = k;
        for (int i = k + 1; i < 4; ++i)
            if (std::abs(A[i][k]) > std::abs(A[pivot][k])) pivot = i;
        if (std::abs(A[pivot][k]) < 1e-12) throw std::runtime_error("Singular system");
        std::swap(A[k], A[pivot]);
        std::swap(b[k], b[pivot]);

        double diag = A[k][k];
        for (int j = k; j < 4; ++j) A[k][j] /= diag;
        b[k] /= diag;

        for (int i = 0; i < 4; ++i) {
            if (i == k) continue;
            double factor = A[i][k];
            for (int j = k; j < 4; ++j) A[i][j] -= factor * A[k][j];
            b[i] -= factor * b[k];
        }
    }
    return b;
}

Mat2 basis(int index) {
    Mat2 E{ { {0.0, 0.0}, {0.0, 0.0} } };
    E[index / 2][index % 2] = 1.0;
    return E;
}

Vec4 flatten(const Mat2& X) {
    return {X[0][0], X[0][1], X[1][0], X[1][1]};
}

// Solve A W + W A^T + Q = 0.
Mat2 lyapunovSolve(const Mat2& A, const Mat2& Q) {
    Mat4 K{};
    for (int col = 0; col < 4; ++col) {
        Mat2 E = basis(col);
        Mat2 L = add(mul(A, E), mul(E, transpose(A)));
        Vec4 f = flatten(L);
        for (int row = 0; row < 4; ++row) K[row][col] = f[row];
    }

    Vec4 rhs = flatten(Q);
    for (double& v : rhs) v = -v;
    Vec4 sol = solve4(K, rhs);

    return { { {sol[0], sol[1]}, {sol[2], sol[3]} } };
}

std::array<double, 2> eigenvalues2(const Mat2& X) {
    double tr = X[0][0] + X[1][1];
    double det = X[0][0] * X[1][1] - X[0][1] * X[1][0];
    double disc = std::max(0.0, tr * tr - 4.0 * det);
    double r1 = 0.5 * (tr + std::sqrt(disc));
    double r2 = 0.5 * (tr - std::sqrt(disc));
    return (r1 >= r2) ? std::array<double,2>{r1, r2} : std::array<double,2>{r2, r1};
}

int main() {
    Mat2 A{ { {-1.0, 0.3}, {0.0, -2.0} } };

    // B = [1.0; 0.5], so B B^T:
    Mat2 BBt{ { {1.0, 0.5}, {0.5, 0.25} } };

    // C = [1.0, -0.2], so C^T C:
    Mat2 CtC{ { {1.0, -0.2}, {-0.2, 0.04} } };

    Mat2 Wc = lyapunovSolve(A, BBt);
    Mat2 Wo = lyapunovSolve(transpose(A), CtC);
    Mat2 P = mul(Wc, Wo);
    auto lambda = eigenvalues2(P);

    std::cout << std::fixed << std::setprecision(6);
    printMat("Controllability Gramian Wc:", Wc);
    printMat("\nOutput-energy Gramian Wo:", Wo);
    printMat("\nProduct Wc*Wo:", P);

    std::cout << "\nHankel singular values:\n";
    std::cout << "  sigma1 = " << std::sqrt(std::max(0.0, lambda[0])) << "\n";
    std::cout << "  sigma2 = " << std::sqrt(std::max(0.0, lambda[1])) << "\n";

    std::cout << "\nInterpretation: large sigma means a state direction is both reachable "
              << "with small input energy and visible at the output.\n";
    return 0;
}

      

The Java version mirrors the C++ implementation: no external linear algebra package is required for this two-state demonstration.

Chapter12_Lesson5.java

// Chapter12_Lesson5.java
// Gramian and Hankel-singular-value preview without external libraries.
// Compile: javac Chapter12_Lesson5.java
// Run:     java Chapter12_Lesson5

public class Chapter12_Lesson5 {
    static double[][] transpose(double[][] X) {
        return new double[][]{ {X[0][0], X[1][0]}, {X[0][1], X[1][1]} };
    }

    static double[][] add(double[][] X, double[][] Y) {
        double[][] Z = new double[2][2];
        for (int i = 0; i < 2; i++)
            for (int j = 0; j < 2; j++)
                Z[i][j] = X[i][j] + Y[i][j];
        return Z;
    }

    static double[][] mul(double[][] X, double[][] Y) {
        double[][] Z = new double[2][2];
        for (int i = 0; i < 2; i++)
            for (int j = 0; j < 2; j++)
                for (int k = 0; k < 2; k++)
                    Z[i][j] += X[i][k] * Y[k][j];
        return Z;
    }

    static double[] flatten(double[][] X) {
        return new double[]{X[0][0], X[0][1], X[1][0], X[1][1]};
    }

    static double[][] basis(int index) {
        double[][] E = new double[2][2];
        E[index / 2][index % 2] = 1.0;
        return E;
    }

    static double[] solve4(double[][] A, double[] b) {
        int n = 4;
        for (int k = 0; k < n; k++) {
            int pivot = k;
            for (int i = k + 1; i < n; i++)
                if (Math.abs(A[i][k]) > Math.abs(A[pivot][k])) pivot = i;

            if (Math.abs(A[pivot][k]) < 1e-12)
                throw new RuntimeException("Singular system");

            double[] tempRow = A[k]; A[k] = A[pivot]; A[pivot] = tempRow;
            double tempVal = b[k]; b[k] = b[pivot]; b[pivot] = tempVal;

            double diag = A[k][k];
            for (int j = k; j < n; j++) A[k][j] /= diag;
            b[k] /= diag;

            for (int i = 0; i < n; i++) {
                if (i == k) continue;
                double factor = A[i][k];
                for (int j = k; j < n; j++) A[i][j] -= factor * A[k][j];
                b[i] -= factor * b[k];
            }
        }
        return b;
    }

    // Solve A W + W A^T + Q = 0.
    static double[][] lyapunovSolve(double[][] A, double[][] Q) {
        double[][] K = new double[4][4];
        for (int col = 0; col < 4; col++) {
            double[][] E = basis(col);
            double[][] L = add(mul(A, E), mul(E, transpose(A)));
            double[] f = flatten(L);
            for (int row = 0; row < 4; row++) K[row][col] = f[row];
        }

        double[] rhs = flatten(Q);
        for (int i = 0; i < rhs.length; i++) rhs[i] = -rhs[i];

        double[] sol = solve4(K, rhs);
        return new double[][]{ {sol[0], sol[1]}, {sol[2], sol[3]} };
    }

    static double[] eigenvalues2(double[][] X) {
        double tr = X[0][0] + X[1][1];
        double det = X[0][0] * X[1][1] - X[0][1] * X[1][0];
        double disc = Math.max(0.0, tr * tr - 4.0 * det);
        double r1 = 0.5 * (tr + Math.sqrt(disc));
        double r2 = 0.5 * (tr - Math.sqrt(disc));
        return (r1 >= r2) ? new double[]{r1, r2} : new double[]{r2, r1};
    }

    static void printMat(String name, double[][] X) {
        System.out.println(name);
        for (int i = 0; i < 2; i++)
            System.out.printf("  %12.6f %12.6f%n", X[i][0], X[i][1]);
    }

    public static void main(String[] args) {
        double[][] A = { {-1.0, 0.3}, {0.0, -2.0} };

        // B = [1.0; 0.5], so B B^T:
        double[][] BBt = { {1.0, 0.5}, {0.5, 0.25} };

        // C = [1.0, -0.2], so C^T C:
        double[][] CtC = { {1.0, -0.2}, {-0.2, 0.04} };

        double[][] Wc = lyapunovSolve(A, BBt);
        double[][] Wo = lyapunovSolve(transpose(A), CtC);
        double[][] product = mul(Wc, Wo);
        double[] lambda = eigenvalues2(product);

        printMat("Controllability Gramian Wc:", Wc);
        printMat("\nOutput-energy Gramian Wo:", Wo);
        printMat("\nProduct Wc*Wo:", product);

        System.out.println("\nHankel singular values:");
        System.out.printf("  sigma1 = %.6f%n", Math.sqrt(Math.max(0.0, lambda[0])));
        System.out.printf("  sigma2 = %.6f%n", Math.sqrt(Math.max(0.0, lambda[1])));

        System.out.println("\nInterpretation: large sigma means a state direction is both reachable with small input energy and visible at the output.");
    }
}

      

MATLAB provides direct Lyapunov solvers and state-space objects. The optional Simulink section creates a State-Space block using the balanced matrices.

Chapter12_Lesson5.m

% Chapter12_Lesson5.m
% Balanced-realization preview for a continuous-time stable LTI system.
% Requires Control System Toolbox for lyap, ss, and optional Simulink model setup.

clear; clc; close all;

A = [-1.0  0.3;
      0.0 -2.0];
B = [1.0; 0.5];
C = [1.0 -0.2];
D = 0;

% Controllability Gramian:
% A*Wc + Wc*A' + B*B' = 0
Wc = lyap(A, B*B');

% Output-energy / observability Gramian preview:
% A'*Wo + Wo*A + C'*C = 0
Wo = lyap(A', C'*C);

Wc = 0.5*(Wc + Wc');
Wo = 0.5*(Wo + Wo');

% Stable computation of Hankel singular values and balancing transform.
R = chol(Wc, 'lower');
[V, Lambda] = eig(R' * Wo * R);
lambda = diag(Lambda);
[lambda, idx] = sort(lambda, 'descend');
V = V(:, idx);
sigma = sqrt(max(lambda, 0));
Sigma = diag(sigma);

% Balancing transformation for x = T*z.
T = R * V * diag(1 ./ sqrt(sigma));
Tinv = inv(T);

Ab = Tinv * A * T;
Bb = Tinv * B;
Cb = C * T;

Wc_bal = Tinv * Wc * Tinv';
Wo_bal = T' * Wo * T;

disp('Controllability Gramian Wc ='); disp(Wc);
disp('Output-energy Gramian Wo ='); disp(Wo);
disp('Hankel singular values sigma ='); disp(sigma);
disp('Balanced-coordinate Wc ='); disp(Wc_bal);
disp('Balanced-coordinate Wo ='); disp(Wo_bal);
disp('Target Sigma ='); disp(Sigma);

% State-space objects before and after coordinate transformation.
sys_original = ss(A, B, C, D);
sys_balanced = ss(Ab, Bb, Cb, D);

% Simulink setup preview:
% This creates a simple model with one State-Space block using balanced matrices.
make_simulink_model = false;
if make_simulink_model
    mdl = 'Chapter12_Lesson5_Simulink';
    if bdIsLoaded(mdl), close_system(mdl, 0); end
    new_system(mdl);
    add_block('simulink/Sources/Step', [mdl '/Step Input']);
    add_block('simulink/Continuous/State-Space', [mdl '/Balanced State Space']);
    add_block('simulink/Sinks/Scope', [mdl '/Scope']);

    set_param([mdl '/Balanced State Space'], ...
        'A', mat2str(Ab), ...
        'B', mat2str(Bb), ...
        'C', mat2str(Cb), ...
        'D', mat2str(D));

    add_line(mdl, 'Step Input/1', 'Balanced State Space/1');
    add_line(mdl, 'Balanced State Space/1', 'Scope/1');
    save_system(mdl);
    open_system(mdl);
end

      

The Mathematica code solves the Lyapunov equations by vectorization, computes Hankel singular values, and constructs the balancing transformation.

Chapter12_Lesson5.nb

(* Chapter12_Lesson5.nb *)
(* Balanced-realization preview for a continuous-time stable LTI system. *)

ClearAll["Global`*"];

A = { {-1.0, 0.3}, {0.0, -2.0} };
B = { {1.0}, {0.5} };
Cmat = { {1.0, -0.2} };

(* Solve A.W + W.Transpose[A] + Q == 0 by vectorization. *)
lyapunovSolve[Amat_, Qmat_] := Module[
  {n, K, rhs, x},
  n = Length[Amat];
  K = KroneckerProduct[IdentityMatrix[n], Amat] +
      KroneckerProduct[Amat, IdentityMatrix[n]];
  rhs = -Flatten[Qmat];
  x = LinearSolve[K, rhs];
  ArrayReshape[x, {n, n}]
];

Wc = lyapunovSolve[A, B.Transpose[B]];
Wo = lyapunovSolve[Transpose[A], Transpose[Cmat].Cmat];

Wc = (Wc + Transpose[Wc])/2;
Wo = (Wo + Transpose[Wo])/2;

R = Transpose[CholeskyDecomposition[Wc]];
M = Transpose[R].Wo.R;

{vals, vecs} = Eigensystem[M];
ord = Reverse[Ordering[vals]];
vals = vals[[ord]];
vecs = vecs[[ord]];
U = Transpose[vecs];

sigma = Sqrt[vals];
Sigma = DiagonalMatrix[sigma];

(* Balancing transformation for x = T.z. *)
Tmat = R.U.DiagonalMatrix[1/Sqrt[sigma]];
Tinv = Inverse[Tmat];

Abal = Tinv.A.Tmat;
Bbal = Tinv.B;
Cbal = Cmat.Tmat;

WcBal = Tinv.Wc.Transpose[Tinv];
WoBal = Transpose[Tmat].Wo.Tmat;

Print["Controllability Gramian Wc = "]; MatrixForm[Wc]
Print["Output-energy Gramian Wo = "]; MatrixForm[Wo]
Print["Hankel singular values = ", sigma]
Print["Balanced Wc = "]; MatrixForm[WcBal]
Print["Balanced Wo = "]; MatrixForm[WoBal]
Print["Target Sigma = "]; MatrixForm[Sigma]

      

Problem 1 (Coordinate Invariance of Minimum Energy): Let \( \mathbf{x}=\mathbf{T}\mathbf{z} \) and suppose \( \bar{\mathbf{W} }_c= \mathbf{T}^{-1}\mathbf{W}_c\mathbf{T}^{-T} \). Show that \( \mathbf{x}^T\mathbf{W}_c^{-1}\mathbf{x}= \mathbf{z}^T\bar{\mathbf{W} }_c^{-1}\mathbf{z} \).

Solution: Since

\[ \bar{\mathbf{W} }_c^{-1} = \mathbf{T}^T\mathbf{W}_c^{-1}\mathbf{T}, \]

we obtain

\[ \mathbf{z}^T\bar{\mathbf{W} }_c^{-1}\mathbf{z} = \mathbf{z}^T\mathbf{T}^T\mathbf{W}_c^{-1}\mathbf{T}\mathbf{z} = (\mathbf{T}\mathbf{z})^T\mathbf{W}_c^{-1} (\mathbf{T}\mathbf{z}) = \mathbf{x}^T\mathbf{W}_c^{-1}\mathbf{x}. \]

Therefore, the numerical value of the minimum energy is coordinate invariant even though the Gramian matrix changes.

Problem 2 (Diagonal Stable System): Consider \( \mathbf{A}=\operatorname{diag}(-a_1,-a_2) \), \( \mathbf{B}=\mathbf{I} \), with \( a_1,a_2 \) positive. Compute \( \mathbf{W}_c \).

Solution: Because the modes are decoupled, write

\[ \mathbf{W}_c= \operatorname{diag}(w_1,w_2). \]

Substituting into \( \mathbf{A}\mathbf{W}_c+ \mathbf{W}_c\mathbf{A}^T+\mathbf{I}=\mathbf{0} \) gives

\[ -2a_iw_i+1=0,\quad w_i=\frac{1}{2a_i},\quad i=1,2. \]

Hence

\[ \mathbf{W}_c= \operatorname{diag}\left(\frac{1}{2a_1}, \frac{1}{2a_2}\right). \]

Faster modes with larger \( a_i \) have smaller controllability Gramian entries under the same identity input scaling.

Problem 3 (Balanced Scaling for a Diagonal SISO Pair): Suppose \( A_i=-a_i \), \( B_i=b_i \), and \( C_i=c_i \) for one decoupled scalar mode. Find the scalar balanced coordinate scaling.

Solution: The scalar Gramians are

\[ W_{c,i}=\frac{b_i^2}{2a_i},\qquad W_{o,i}=\frac{c_i^2}{2a_i}. \]

Let \( x_i=t_i z_i \). Then

\[ \bar{W}_{c,i}=t_i^{-2}W_{c,i},\qquad \bar{W}_{o,i}=t_i^2W_{o,i}. \]

Balance requires \( t_i^{-2}W_{c,i}=t_i^2W_{o,i} \), so

\[ t_i=\left(\frac{W_{c,i} }{W_{o,i} }\right)^{1/4}, \qquad \bar{W}_{c,i}=\bar{W}_{o,i}= \sqrt{W_{c,i}W_{o,i} } = \frac{|b_ic_i|}{2a_i}. \]

The common balanced value is exactly the scalar Hankel singular value.

Problem 4 (Meaning of Small Hankel Singular Values): In balanced coordinates, explain why a small \( \sigma_i \) indicates a state coordinate that is both difficult to excite and weakly visible.

Solution: For a coordinate displacement \( \alpha\mathbf{e}_i \), the input energy and output energy are

\[ \mathcal{E}_{u,\min}= \frac{\alpha^2}{\sigma_i},\qquad \mathcal{E}_y=\alpha^2\sigma_i. \]

If \( \sigma_i \) is small, the required input energy is large, while the produced output energy is small. This is why the same scalar captures both poor reachability and poor output visibility.

The controllability Gramian gives a coordinate-dependent matrix representation of a coordinate-invariant input-energy quantity. Balanced realizations combine this input-energy geometry with a previewed output-energy geometry. In balanced coordinates, both Gramians are equal to the same diagonal matrix of Hankel singular values. Large values identify directions that are both easy to reach and strongly visible, while small values identify directions that are energetically difficult and weakly visible. No model reduction was performed in this lesson; the purpose was to prepare the conceptual bridge from Gramians to balanced coordinates.