Balanced truncation of a large thermal system: a certified 5-state surrogate for a 100-state rod¶
A finite-difference model of a heated rod hands you a state vector with one temperature per node. Refine the mesh for accuracy and the state dimension balloons — a hundred nodes here, ten thousand for a 3-D casting — while the input/output behaviour you actually care about (heater in, thermocouple out) stays low-dimensional. Balanced truncation collapses that gap: it produces a small state-space model that reproduces the transfer function of the large one, with an a priori bound on the worst-case error you incur. This notebook builds a 100-state rod, ranks its internal modes by Hankel singular value, truncates to a 5-state surrogate with jaxonomy.library.rom, certifies the error against the theoretical bound, and simulates both models side by side in jaxonomy.simulate.
What you will be able to do afterwards. Discretize a 1-D PDE into a state-space LTI; read a Hankel-singular-value spectrum and choose a reduction order from it; run balanced truncation (and its modal / residualization cousins) on a LinearizedSystem; verify the $H_\infty$ error bound holds; and drop the reduced LTISystem straight back into a simulation.
Estimated reading time: ~30 min. Runtime on CPU: ~30 s from a fresh kernel.
Prerequisites¶
- State-space LTI models — $\dot{\mathbf x}=A\mathbf x+B\mathbf u,\ \mathbf y=C\mathbf x+D\mathbf u$, transfer functions, poles, DC gain.
LinearizedSystem/LTISystem— Jaxonomy's frozen-linear-model container and its simulatable block. If these are new, skim the linearization workflow notebook first; we reuse itsLinearizedSystem.to_scipy_lti()bridge for Bode and step plots.- Lyapunov equations — helpful but not required; we state everything we use.
The reduction routines live in jaxonomy.library.rom and are re-exported from jaxonomy.library.
The problem: a model that is too big for its own dynamics¶
Picture a thin metal rod — a busbar in a power module, say — heated by a resistive element clamped a quarter of the way along it, with a thermocouple reading temperature three quarters of the way along. Both ends are bolted to a cold chassis that holds them at ambient. We want a dynamic model of heater power in $\to$ probe temperature out to design a controller and run it in real time.
First-principles heat conduction is a partial differential equation. To simulate it we discretize space into $n$ nodes, turning the PDE into $n$ coupled ODEs — one temperature state per node. Accuracy wants $n$ large; a real-time controller and an observer want $n$ small. At $n=100$ the model is already too heavy to embed comfortably, yet the map from a single heater to a single probe cannot possibly need 100 states to describe. That mismatch — a high-dimensional realization of a low-dimensional input/output map — is exactly what model-order reduction exploits.
import warnings
import numpy as np
import scipy.signal as sig
import matplotlib.pyplot as plt
import jaxonomy
from jaxonomy import DiagramBuilder, library, SimulatorOptions
from jaxonomy.library import (
LinearizedSystem,
reduce,
balred,
modal_truncation,
residualize,
hankel_singular_values,
)
from jaxonomy import logging as jaxlog
jaxlog.set_log_level(jaxlog.ERROR) # keep the simulator quiet in notebook output
warnings.filterwarnings("ignore", category=sig.BadCoefficients)
np.set_printoptions(precision=4, suppress=True)
Mathematical model¶
From the heat equation to a state-space LTI¶
Let $\theta(z,t)$ be the temperature rise above ambient at position $z\in[0,L]$ along a rod of length $L$. One-dimensional heat conduction with a distributed source $q(z,t)$ is
$$\frac{\partial \theta}{\partial t} \;=\; \alpha\,\frac{\partial^2 \theta}{\partial z^2} \;+\; \frac{1}{\rho c_p}\,q(z,t), \tag{1}$$
with thermal diffusivity $\alpha=k/(\rho c_p)$. The ends are clamped to the chassis (ambient), giving homogeneous Dirichlet boundary conditions $\theta(0,t)=\theta(L,t)=0$. Dimensionally $[\alpha]=\mathrm{m^2/s}$ and $[\alpha\,\partial^2\theta/\partial z^2]=\mathrm{K/s}$ — a rate of temperature change, as a transient balance must be.
We place $n$ interior nodes at spacing $\Delta z = L/(n+1)$ and approximate the second derivative by the standard three-point stencil (method of lines):
$$\left.\frac{\partial^2\theta}{\partial z^2}\right|_{i} \approx \frac{\theta_{i-1}-2\theta_i+\theta_{i+1}}{\Delta z^2}. \tag{2}$$
Collecting the node temperatures into $\mathbf x=(\theta_1,\dots,\theta_n)^\top$ turns (1) into a linear time-invariant system
$$\dot{\mathbf x} = A\,\mathbf x + B\,u,\qquad y = C\,\mathbf x, \tag{3}$$
where $A=\dfrac{\alpha}{\Delta z^2}\,\mathrm{tridiag}(1,-2,1)$ is the discrete Laplacian, the input $u$ is the heater flux injected at node $i_\text{in}$ (so $B$ is a scaled unit column), and the output $y$ is the temperature at the probe node $i_\text{out}$ (so $C$ is a unit row). The clamped ends never enter $\mathbf x$, which is what makes $A$ Hurwitz (all eigenvalues negative) rather than merely marginally stable — a fully insulated rod would instead conserve total heat, leaving a zero eigenvalue and unbounded Gramians.
| symbol | meaning | units |
|---|---|---|
| $\theta_i$ | temperature rise at node $i$ | K |
| $\alpha$ | thermal diffusivity | m²/s |
| $L,\ \Delta z$ | rod length, node spacing | m |
| $n$ | number of interior nodes (= state dimension) | – |
| $A$ | discrete Laplacian $\tfrac{\alpha}{\Delta z^2}\mathrm{tridiag}(1,-2,1)$ | 1/s |
| $B$ | heater-injection column at node $i_\text{in}$ | K/s per input unit |
| $C$ | probe-selection row at node $i_\text{out}$ | – |
| $u,\ y$ | heater input, probe temperature | (input unit), K |
| $\sigma_i$ | Hankel singular value $i$ | K |
| $r$ | reduced order | – |
def build_rod(n=100, alpha=1.0e-2, L=1.0, i_in=None, i_out=None):
'''1-D heat-conduction rod, Dirichlet ends, as a state-space (A, B, C, D). n interior
nodes -> an n-state stable LTI mapping heater flux at i_in to temperature at i_out.'''
i_in = n // 4 if i_in is None else i_in # heater a quarter along
i_out = 3 * n // 4 if i_out is None else i_out # probe three quarters along
dz = L / (n + 1)
lap = np.diag(-2.0 * np.ones(n)) + np.diag(np.ones(n - 1), 1) + np.diag(np.ones(n - 1), -1)
A = lap * alpha / dz**2
B = np.zeros((n, 1)); B[i_in, 0] = 1.0 / dz # flux -> temperature-rate
C = np.zeros((1, n)); C[0, i_out] = 1.0 # read one node
D = np.zeros((1, 1))
return A, B, C, D
N = 100
A, B, C, D = build_rod(N)
poles = np.linalg.eigvals(A)
print(f"state dimension n = {A.shape[0]}")
print(f"slowest pole (dominant time constant): {poles.real.max():.4f} 1/s "
f"-> tau ~= {-1.0/poles.real.max():.2f} s")
print(f"fastest pole: {poles.real.min():.2f} 1/s")
print(f"stable (all Re(lambda) < 0)? {np.all(poles.real < 0)}")
state dimension n = 100 slowest pole (dominant time constant): -0.0987 1/s -> tau ~= 10.13 s fastest pole: -407.94 1/s stable (all Re(lambda) < 0)? True
The spectrum spans two decades of time constants: a slow ~10 s mode governs the settling we care about, while nodes near the stencil resolution relax in milliseconds. That stiffness is the physical reason the model is compressible — the fast modes are barely excited by a single heater and barely seen by a single probe.
We wrap the matrices in a LinearizedSystem, the frozen-linear-model container the reduction routines consume. Reduction is an analysis-time operation on a fixed realization, so the rom numerics run host-side in NumPy/SciPy; the result comes back as a simulatable Jaxonomy block.
full = LinearizedSystem(A, B, C, D, {}) # operating point is empty (already at ambient)
dc_gain = float((-C @ np.linalg.solve(A, B) + D).ravel()[0]) # steady probe temp for a unit step, y = -C A^{-1} B u
print(f"full model: {full.A.shape[0]} states, stable={full.is_stable()}")
print(f"DC gain (steady probe temperature per unit heater input): {dc_gain:.4f} K")
full model: 100 states, stable=True DC gain (steady probe temperature per unit heater input): 6.3719 K
Which states matter? Gramians and Hankel singular values¶
Not all internal states pull equal weight in the input/output map. Two energies quantify a state's relevance:
- the controllability Gramian $W_c$ measures how strongly the input can move each direction of state space,
- the observability Gramian $W_o$ measures how strongly each direction shows up at the output.
For a stable continuous-time system they solve the Lyapunov equations
$$A W_c + W_c A^\top + BB^\top = 0,\qquad A^\top W_o + W_o A + C^\top C = 0. \tag{4}$$
A state direction only matters to the transfer function if it is both reachable and observable. The invariant that fuses the two is the set of Hankel singular values
$$\sigma_i \;=\; \sqrt{\lambda_i(W_c W_o)}, \qquad \sigma_1 \ge \sigma_2 \ge \dots \ge \sigma_n \ge 0. \tag{5}$$
Each $\sigma_i$ is the "input-to-output energy gain" of one balanced mode; they are invariant under any change of state coordinates, so they are an intrinsic ranking of what the system actually does. A sharp drop in the spectrum is a green light for truncation.
hsv = hankel_singular_values(full)
print("first 10 Hankel singular values:")
for i, s in enumerate(hsv[:10], start=1):
print(f" sigma_{i:<2d} = {s:.3e}")
energy = np.cumsum(hsv**2) / np.sum(hsv**2)
print(f"\nenergy captured by the first 5 modes: {100*energy[4]:.4f}%")
first 10 Hankel singular values: sigma_1 = 3.932e+00 sigma_2 = 8.270e-01 sigma_3 = 8.590e-02 sigma_4 = 6.087e-03 sigma_5 = 9.392e-04 sigma_6 = 1.221e-04 sigma_7 = 5.217e-05 sigma_8 = 5.690e-06 sigma_9 = 9.160e-07 sigma_10 = 1.259e-07 energy captured by the first 5 modes: 100.0000%
fig, ax = plt.subplots(figsize=(7, 4.2))
idx = np.arange(1, len(hsv) + 1)
ax.semilogy(idx, hsv, "o-", ms=4, color="tab:blue")
r = 5
ax.axvline(r + 0.5, color="tab:red", ls="--", lw=1.2, label=f"truncation at r = {r}")
ax.set_xlim(0, 25)
ax.set_xlabel("mode index $i$")
ax.set_ylabel(r"Hankel singular value $\sigma_i$ [K]")
ax.set_title("Hankel singular-value spectrum of the 100-state rod")
ax.grid(True, which="both", alpha=0.3)
ax.legend()
plt.tight_layout(); plt.show()
Figure 1: Hankel singular values on a log scale (blue), truncation cut at $r=5$ (red dashed). The spectrum falls roughly an order of magnitude per mode: $\sigma_1\approx 3.9$, $\sigma_5\approx 9\times10^{-4}$, and everything past the tenth mode is below $10^{-6}$ K. The first five modes already carry essentially all of the input/output energy — a textbook signal that a 5-state surrogate can stand in for the 100-state model.
The super-exponential decay is characteristic of diffusion: successive modes are both harder to reach and harder to see, so their product collapses fast. We take $r=5$.
Balanced truncation¶
Balanced truncation, introduced by Moore (1981), works in two steps. First find the state coordinates in which the two Gramians are equal and diagonal, $W_c=W_o=\mathrm{diag}(\sigma_1,\dots,\sigma_n)$ — the balanced realization, computed stably by the square-root algorithm of Laub, Heath, Paige & Ward (1987). In these coordinates every state is exactly as controllable as it is observable, and its importance is its Hankel singular value. Then simply delete the states with the smallest $\sigma_i$.
The reason this is the workhorse of linear MOR is the accompanying a priori guarantee (Glover 1984, Enns 1984): the reduced model $G_r$ satisfies
$$\lVert G - G_r\rVert_\infty \;\le\; 2\sum_{i=r+1}^{n}\sigma_i. \tag{6}$$
The worst-case gain error, across all frequencies, is bounded by twice the tail of the discarded singular values — a number you can read off before building the reduced model. Jaxonomy exposes this as balred (a MATLAB-style alias), or through the unified reduce(...) front door, which returns a ReducedOrderModel carrying the reduced system plus provenance in .info.
The one-line API, on a small rod¶
The whole workflow is a single call. Here is a coarse 6-node rod reduced to 3 states, so you can see the shape of the result — the reduced LTISystem, the error bound, and the retained spectrum:
demo = reduce(build_rod(n=6), method="balred", order=3)
print(demo)
print(f" reduced system : {type(demo.system).__name__} ({demo.reduced_order} states)")
print(f" H-inf bound : {demo.info['error_bound']:.3e}")
print(f" retained HSV : {np.array2string(demo.info['hsv'][:3], precision=3)}")
ReducedOrderModel(method='balred', full_order=6, reduced_order=3) reduced system : LTISystem (3 states) H-inf bound : 4.303e-03 retained HSV : [4.927 0.905 0.061]
The full 100-state reduction¶
Now the real model. The same one-line call takes the 100-state rod down to five states. Balanced truncation here forms the Gramian square roots from a symmetric eigendecomposition, so it stays valid even where the Hankel spectrum underflows to zero (Figure 1) and the Gramians are numerically rank-deficient — exactly the regime a stiff single-input/single-output diffusion model lands in.
reduced = balred(full, order=5)
error_bound = float(2.0 * hsv[5:].sum()) # a priori H-inf bound, Eq. (6)
print(f"reduced order : {reduced.A.shape[0]}")
print(f"reduced poles : {np.sort(np.linalg.eigvals(np.asarray(reduced.A)).real)}")
print(f"a priori H-inf error bound 2*sum_(i>5) sigma_i = {error_bound:.3e} K")
reduced order : 5 reduced poles : [-2.2098 -2.2098 -0.9555 -0.3922 -0.0987] a priori H-inf error bound 2*sum_(i>5) sigma_i = 3.632e-04 K
The five retained poles span the same slow-to-fast range as the dominant modes of the full model — balancing kept the input/output-relevant dynamics and discarded the rest.
Verifying the a priori bound¶
The bound in Eq. (6) is a promise: the peak of $|G(j\omega)-G_r(j\omega)|$ over all frequencies must not exceed $2\sum_{i>r}\sigma_i$. We measure the left-hand side directly by sweeping the frequency response of both models — Jaxonomy's LinearizedSystem.to_scipy_lti() hands us a SciPy LTI object for exactly this kind of analysis.
w = np.logspace(-3, 3, 800) # rad/s
_, H_full = full.to_scipy_lti().freqresp(w)
_, H_red = reduced.to_scipy_lti().freqresp(w)
measured_hinf = float(np.max(np.abs(H_full - H_red)))
print(f"measured ||G - G_r||_inf = {measured_hinf:.3e} K")
print(f"a priori bound = {error_bound:.3e} K")
assert measured_hinf <= error_bound, "a priori H-inf bound violated!"
print(f"bound holds: measured is {error_bound/measured_hinf:.1f}x below the certificate")
measured ||G - G_r||_inf = 1.500e-04 K a priori bound = 3.632e-04 K bound holds: measured is 2.4x below the certificate
The measured worst-case error comes in comfortably under the certificate — as it must, since the bound is guaranteed, and in practice conservative. Now let us see the agreement, on a Bode plot.
mag_full, ph_full = np.abs(H_full), np.unwrap(np.angle(H_full))
mag_red, ph_red = np.abs(H_red), np.unwrap(np.angle(H_red))
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(7.2, 6), sharex=True)
ax1.loglog(w, mag_full, color="tab:blue", lw=2.4, alpha=0.55, label="full (100 states)")
ax1.loglog(w, mag_red, color="tab:red", lw=1.4, ls="--", label="reduced (5 states)")
ax1.set_ylabel("magnitude [K / input]"); ax1.grid(True, which="both", alpha=0.3); ax1.legend()
ax1.set_title("Bode plot: full vs. balanced-truncation surrogate")
ax2.semilogx(w, np.degrees(ph_full), color="tab:blue", lw=2.4, alpha=0.55)
ax2.semilogx(w, np.degrees(ph_red), color="tab:red", lw=1.4, ls="--")
ax2.set_ylabel("phase [deg]"); ax2.set_xlabel(r"frequency $\omega$ [rad/s]")
ax2.grid(True, which="both", alpha=0.3)
plt.tight_layout(); plt.show()
Figure 2: Frequency response of the full 100-state rod (blue, thick) and the 5-state balanced surrogate (red, dashed). Magnitude (top) and phase (bottom) overlie across six decades; the reduced model matches the DC gain, the dominant roll-off, and the high-frequency slope. The two curves are visually indistinguishable wherever the response has appreciable gain — the $10^{-4}$-scale mismatch certified above lives where $|G|$ has already fallen far below its DC value.
T_end = 50.0
t = np.linspace(0.0, T_end, 500)
_, y_full_step = sig.step(full.to_scipy_lti(), T=t)
_, y_red_step = sig.step(reduced.to_scipy_lti(), T=t)
fig, ax = plt.subplots(figsize=(7, 4.2))
ax.plot(t, y_full_step, color="tab:blue", lw=2.4, alpha=0.55, label="full (100 states)")
ax.plot(t, y_red_step, color="tab:red", lw=1.4, ls="--", label="reduced (5 states)")
ax.axhline(dc_gain, color="k", lw=0.8, ls=":", label=f"DC gain = {dc_gain:.3f} K")
ax.set_xlabel("time [s]"); ax.set_ylabel("probe temperature [K]")
ax.set_title("Step response: unit heater step at $t=0$")
ax.grid(True, alpha=0.3); ax.legend()
plt.tight_layout(); plt.show()
print(f"peak step-response error: {np.max(np.abs(y_full_step - y_red_step)):.3e} K")
peak step-response error: 7.327e-05 K
Figure 3: Probe temperature after a unit heater step. The blue full-order response and the red 5-state surrogate track through the diffusive delay (temperature takes seconds to travel from heater to probe), the rise, and the settle to the DC gain (dotted). Peak deviation is a few millikelvin over a ~6 K swing.
Same input, both models, in jaxonomy.simulate¶
The reduced LinearizedSystem becomes a first-class simulation block via .to_lti(). We now leave the frequency domain entirely: build a tiny diagram — a constant heater driving an LTISystem — for the full model and for the reduced model, and integrate both with the same solver under the same input. This is the payoff: the 5-state block is a drop-in replacement that simulates at a fraction of the cost.
def simulate_lti(sys_matrices, u_level=1.0, t_end=50.0):
'''Drive an (A,B,C,D) LTISystem with a constant input; return (t, y, wall_time).'''
import time
A, B, C, D = sys_matrices
n = A.shape[0]
builder = DiagramBuilder()
heater = builder.add(library.Constant(np.full(B.shape[1], u_level))) # input port
plant = builder.add(library.LTISystem(A, B, C, D, np.zeros(n))) # x0 = ambient
builder.connect(heater.output_ports[0], plant.input_ports[0])
diagram = builder.build()
ctx = diagram.create_context()
opts = SimulatorOptions(buffer_length=8192) # hold every sample of the stiff full-model trace
t0 = time.perf_counter()
res = jaxonomy.simulate(diagram, ctx, (0.0, t_end), options=opts,
recorded_signals={"y": plant.output_ports[0]})
return np.asarray(res.time), np.asarray(res.outputs["y"]).ravel(), time.perf_counter() - t0
full_mats = (A, B, C, D)
red_mats = tuple(np.asarray(M) for M in (reduced.to_lti().A, reduced.to_lti().B,
reduced.to_lti().C, reduced.to_lti().D))
t_full, y_full, wt_full = simulate_lti(full_mats)
t_red, y_red, wt_red = simulate_lti(red_mats)
print(f"full (100 states): {wt_full:.2f} s wall")
print(f"reduced (5 states): {wt_red:.2f} s wall ({wt_full/max(wt_red,1e-9):.1f}x faster)")
full (100 states): 0.28 s wall reduced (5 states): 0.16 s wall (1.7x faster)
tt = np.linspace(0.0, 50.0, 500)
yf = np.interp(tt, t_full, y_full)
yr = np.interp(tt, t_red, y_red)
sim_err = np.abs(yf - yr)
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(7.2, 6), sharex=True,
gridspec_kw={"height_ratios": [3, 1]})
ax1.plot(tt, yf, color="tab:blue", lw=2.4, alpha=0.55, label="full (100 states)")
ax1.plot(tt, yr, color="tab:red", lw=1.4, ls="--", label="reduced (5 states)")
ax1.set_ylabel("probe temperature [K]")
ax1.set_title("jaxonomy.simulate: full vs. reduced under the same heater input")
ax1.grid(True, alpha=0.3); ax1.legend()
ax2.semilogy(tt, np.maximum(sim_err, 1e-12), color="tab:purple", lw=1.4)
ax2.set_ylabel("|error| [K]"); ax2.set_xlabel("time [s]"); ax2.grid(True, which="both", alpha=0.3)
plt.tight_layout(); plt.show()
l2_ratio = np.linalg.norm(yf - yr) / np.linalg.norm(yf)
print(f"peak simulate error: {sim_err.max():.3e} K ({100*sim_err.max()/dc_gain:.3f}% of DC gain)")
print(f"relative L2 error: {l2_ratio:.3e}")
peak simulate error: 4.885e-03 K (0.077% of DC gain) relative L2 error: 2.684e-04
Figure 4: Time-domain simulation in jaxonomy.simulate. Top: the two output traces overlie (blue full, red dashed reduced). Bottom: the pointwise error stays at the millikelvin level throughout, a fraction of a percent of the ~6 K signal — the transient error is largest during the fast initial rise and decays as the slow dominant mode takes over. The 5-state block integrates several times faster than the 100-state one, and the gap widens with mesh refinement.
The time-domain error is consistent with the frequency-domain certificate: $\lVert y_\text{full}-y_\text{red}\rVert$ is bounded by $\lVert G-G_r\rVert_\infty\,\lVert u\rVert$, and here it lands well inside that envelope.
Two cousins: modal truncation and residualization¶
Balanced truncation is not the only reducer in jaxonomy.library.rom. Two others reduce by keeping the slowest physical modes rather than the most input/output-relevant ones, and they differ in what they do with the fast modes they drop:
modal_truncationtransforms to a block-diagonal modal form and simply deletes the fast modes. Cheap and interpretable — the retained poles are exact eigenvalues — but because it throws away the discarded modes' static contribution, the DC gain generally shifts.residualize(singular-perturbation reduction) instead sets the fast modes' derivatives to zero and solves for their quasi-steady value, folding that back into the retained model. This preserves the DC gain exactly, at the cost of the high-frequency fidelity that balanced truncation protects.
Both run directly on the 100-state model. Let us reduce to $r=5$ with each and compare DC gains.
modal5 = modal_truncation(full, order=5)
resid5 = residualize(full, order=5)
def dc_of(sys):
Am, Bm, Cm, Dm = (np.asarray(sys.A), np.asarray(sys.B), np.asarray(sys.C), np.asarray(sys.D))
return float((-Cm @ np.linalg.solve(Am, Bm) + Dm).ravel()[0])
print(f"DC gain, full model : {dc_gain:.5f} K")
print(f"DC gain, balanced truncation (r=5): {dc_of(reduced):.5f} K (error {abs(dc_of(reduced)-dc_gain):.2e})")
print(f"DC gain, modal truncation (r=5): {dc_of(modal5):.5f} K (error {abs(dc_of(modal5)-dc_gain):.2e})")
print(f"DC gain, residualization (r=5): {dc_of(resid5):.5f} K (error {abs(dc_of(resid5)-dc_gain):.2e})")
DC gain, full model : 6.37192 K DC gain, balanced truncation (r=5): 6.37177 K (error 1.50e-04) DC gain, modal truncation (r=5): 6.72964 K (error 3.58e-01) DC gain, residualization (r=5): 6.37192 K (error 2.13e-14)
_, y_modal = sig.step(modal5.to_scipy_lti(), T=t)
_, y_resid = sig.step(resid5.to_scipy_lti(), T=t)
fig, ax = plt.subplots(figsize=(7.2, 4.4))
ax.plot(t, y_full_step, color="tab:blue", lw=2.6, alpha=0.5, label="full (100 states)")
ax.plot(t, y_red_step, color="tab:red", lw=1.4, ls="--", label="balanced truncation")
ax.plot(t, y_modal, color="tab:green", lw=1.4, ls="-.", label="modal truncation")
ax.plot(t, y_resid, color="tab:orange", lw=1.4, ls=":", label="residualization")
ax.axhline(dc_gain, color="k", lw=0.8, ls=":", alpha=0.6)
ax.set_xlabel("time [s]"); ax.set_ylabel("probe temperature [K]")
ax.set_title("Three reducers, r = 5: step response")
ax.grid(True, alpha=0.3); ax.legend()
plt.tight_layout(); plt.show()
Figure 5: Step responses of the three order-5 reducers against the full model (blue). Balanced truncation (red) tracks best overall. Residualization (orange) nails the final value — it is built to match DC — but lags slightly through the transient. Modal truncation (green) captures the transient shape yet settles to a visibly wrong steady state, because it discards the DC contribution of the fast modes it deletes. The right reducer depends on what you need: overall $H_\infty$ fidelity (balanced), exact steady state (residualize), or interpretable physical modes (modal).
Validation summary¶
Two independent checks confirm the surrogate:
- Certified $H_\infty$ error. The measured worst-case gain error is below the a priori bound $2\sum_{i>r}\sigma_i$ — the promise of balanced truncation, verified rather than assumed.
- DC-gain preservation. Residualization reproduces the full model's steady-state gain to machine precision, while modal truncation does not — exactly as the singular-perturbation theory predicts.
Failure modes¶
Reduction is not magic; here is where it breaks, concretely.
- Slowly-decaying Hankel spectrum. If the $\sigma_i$ do not drop off — e.g. a lightly-damped structure or an advection-dominated (high-Péclet) transport line, where many modes are comparably controllable and observable — then $2\sum_{i>r}\sigma_i$ stays large and no small $r$ gives a tight bound. Read the spectrum (Figure 1) before committing to an order.
- Extreme ill-conditioning. When the $\sigma_i$ underflow to zero the Gramians are numerically rank-deficient. Jaxonomy's
balredhandles this by forming its square-root factors from a symmetric eigendecomposition rather than a Cholesky factorization, so balanced truncation stays valid at scale; for very large models ($n$ in the thousands), project onto a modal or POD subspace first to obtain a numerically minimal core, then balance. - Marginal stability. A perfectly insulated rod (Neumann ends) has a conserved-heat mode with eigenvalue zero. The Lyapunov equations (4) then have no bounded solution and the Gramians blow up — balanced truncation assumes a Hurwitz $A$. Reduce the stable part separately, or add the smallest physical loss (a Robin/convective boundary) to regularize.
- Wrong reducer for the spec. Using
modal_truncationwhere steady-state accuracy matters gives the offset in Figure 5. Match the method to the requirement.
Exercises¶
Move the sensor (code). In
build_rod, seti_out = n // 4so the probe sits on the heater node. Re-run the HSV spectrum: how does co-locating input and output change the decay rate, the order-5 error bound, and the reduced poles? Explain physically why a co-located sensor makes the system "more observable" in the modes the heater excites.Sweep the order (code). Loop $r$ from 1 to 10, and for each compute the a priori bound $2\sum_{i>r}\sigma_i$ and the measured $\lVert G-G_r\rVert_\infty$ with
balred. Plot both against $r$ on a log axis. Where does the bound stop improving, and how conservative is it (the ratio) at each order?Beat the diffusion (conceptual). Replace pure conduction with advection–diffusion by adding a first-derivative term $-v\,\partial\theta/\partial z$ (a moving rod or a coolant flow) — an upwind stencil adds a sub-/super-diagonal to $A$. Argue, before computing, whether the Hankel spectrum should decay faster or slower than pure diffusion, and what that does to the achievable reduced order. (Hint: advection makes modes travel rather than dissipate, so they stay observable downstream.)
Open-ended. The
rompackage also offers data-driven reduction —reduce(snapshots, method="dmd")and Koopman/edmd— that build a reduced model from simulated trajectories instead of the matrices. Generate step- and chirp-response snapshots of the full rod, fit a DMD model, and compare its accuracy and its identified eigenvalues against the balanced-truncation surrogate. When would you reach for data-driven reduction over balanced truncation, and what does each assume about access to the model?
Key takeaways¶
- Hankel singular values are the ranking that matters. They fuse controllability and observability into a coordinate-invariant measure of each mode's input/output relevance; their decay (Figure 1) tells you whether and how far a model can be reduced.
- Balanced truncation comes with a certificate. The $H_\infty$ error is bounded a priori by twice the tail of the discarded $\sigma_i$ — computable before you build the reduced model, and verified here against the measured error.
- A 5-state surrogate reproduced a 100-state rod to millikelvin accuracy across six frequency decades and in direct
jaxonomy.simulatetime-domain runs, at a fraction of the cost. - Pick the reducer for the spec: balanced truncation for overall fidelity,
residualizefor exact steady state,modal_truncationfor interpretable physical modes. - Watch the numerics. Stiff, ill-conditioned Gramians are the norm for discretized PDEs; Jaxonomy's
balredforms its square-root factors from a symmetric eigendecomposition so it stays robust there, and a modal or POD pre-projection helps at very large scale.
Next steps¶
- Feed the reduced
LTISysteminto an LQR or Kalman design (see the linearization and control examples) — a 5-state model makes real-time estimation and control tractable. - Explore the data-driven reducers in
jaxonomy.library.rom(dmd,dmdc,edmd) when you have trajectories but not matrices (Exercise 4). - Scale to 2-D/3-D conduction, where $n$ reaches the thousands and reduction stops being optional.
References¶
- B. C. Moore, "Principal component analysis in linear systems: controllability, observability, and model reduction," IEEE Transactions on Automatic Control, 26(1):17–32, 1981.
- A. J. Laub, M. T. Heath, C. C. Paige, R. C. Ward, "Computation of system balancing transformations and other applications of simultaneous diagonalization algorithms," IEEE Transactions on Automatic Control, 32(2):115–122, 1987.
- K. Glover, "All optimal Hankel-norm approximations of linear multivariable systems and their $L^\infty$-error bounds," International Journal of Control, 39(6):1115–1193, 1984.
- D. F. Enns, "Model reduction with balanced realizations: an error bound and a frequency weighted generalization," Proc. 23rd IEEE Conference on Decision and Control, 1984.
- P. Kokotović, H. K. Khalil, J. O'Reilly, Singular Perturbation Methods in Control: Analysis and Design, Academic Press, 1986.
- A. C. Antoulas, Approximation of Large-Scale Dynamical Systems, SIAM, 2005 (comprehensive modern reference).