平衡截断（Balanced Truncation）—— MTALAB 和 Python 实现

平衡截断

平衡截断（Balanced Truncation）是一种经典的 模型降阶方法，它通过衡量 各状态对系统输入–输出响应的贡献（Hankel 奇异值）来丢弃“能量”较小的状态，从而得到低阶近似模型。平衡截断既能 保证降阶后模型与原模型在频域响应上的接近性，也给出了严格的误差界。

在做平衡截断（或其他降阶）之后，不会直接再用原来高维的 状态向量 $x$——而是引入一个新的、维度已降到 $r$ 的降阶状态向量 $x_r$。

balreal 算法原理

balreal 函数 是 MATLAB Control Toolbox 中用于 平衡截断（balanced truncation）模型约简的核心函数。

[sysb, g, T, Ti] = balreal(sys, opts);

其中，sysb 是平衡化后的系统（ss 对象）；g 是 Hankel 奇异值向量，对应 Gramian 对角线；T、Ti：相似变换矩阵及其逆，用于在原始和平衡化坐标间转换。

算法起源与扩展

• 历史：最早由 Moore 提出“平衡化（balancing）”概念，后由 Glover 完善，成为经典的 SVD 基模型约简方法。
• 数学基础：本质上是 对可控性 Gramian 和可观性 Gramian 的同时对角化问题，等价于 对两正定矩阵进行相似对角化。

其基本思路是：

• 将原始系统通过相似变换（similarity transformation）转化到可控性和可观性 Gramian 相等且对角化的坐标下，所得对角线即为 Hankel 奇异值（Hankel singular values）。
• 奇异值较小的状态对系统输入-输出行为贡献较小，可 予以截断以简化模型，同时可利用理论上严格的误差界评估截断误差。

  截断阶段可选择多种策略，包括简单删除（Truncate）和匹配直流增益（MatchDC），并提供基于奇异值的 $\infty$-范数误差上界估计。

  • Truncate：直接删除最后 $n-r$ 个状态，频域近似效果较好，但不保证直流增益匹配。
  • MatchDC（默认）：在截断同时重新计算系统矩阵，保证截断系统与原系统的直流增益一致。

  平衡截断在 ${\mathcal{H}}_{\infty }$ 范数下满足经典误差界：
  $$\Vert G-G_r{\Vert }_{\infty }\le 2\sum _{i=r+1}^n{\sigma }_i$$
  其中 ${\sigma }_i$ 为被截断的 Hankel 奇异值，提供了保留前 $r$ 个状态时最坏情况的频域误差估计。

对于含有 不稳定极点 的系统，balreal（或 balred）会首先调用 stabsep 对系统进行“稳定/不稳定”子系统分离：
$$G(s)=G_s(s)+G_u(s),$$
其中 $G_s$ 是 所有极点实部<0 的稳定子系统，$G_u$ 是不稳定子系统。

• 针对稳定部分 $G_s$ 求解 Lyapunov 方程（李雅普诺夫方程）得到 Gramian，再 通过 Cholesky 分解及奇异值分解（SVD）构造平衡变换，最后输出平衡化系统及相应的 Hankel 奇异值和变换矩阵。
• 不稳定部分 $G_u$ 则原样保留，并在输出中 将对应奇异值设置为无穷大 Inf，以 提示截断时勿删除不稳定状态。

最后将两部分按并联结构拼回。

为什么要先分离 稳定/不稳定 子系统？

1. 稳定性保证：Gramian 只有在系统稳定时才存在有限、正定的解；对不稳定部分直接调用 Lyapunov 会发散或得到非正定结果。
2. 误差界仅对稳定部分：平衡截断的 ${\mathcal{H}}_{\infty }$ 误差界$\Vert G-G_r{\Vert }_{\infty }\le 2{\sum }_{i>r}{\sigma }_i$只对稳定子系统有意义。
3. 保留不稳定行为：不稳定极点通常对系统行为至关重要，不能被截断。

平衡截断过程

下面用坐标变换
$$\bar{x}=Tx$$

（即 $T$ 将原始状态 $x$ 映到“平衡坐标” $\bar{x}$）来重新推导平衡截断的全过程。

对于 原系统，考虑连续时间 LTI 系统
$$\begin{gathered} \dot{x}=Ax+Bu, \\ y=Cx+Du \end{gathered}$$
其中 $x\in {\mathbb{R}}^n$、$u\in {\mathbb{R}}^m$、$y\in {\mathbb{R}}^p$。

1. 求解 Gramian 并构造变换矩阵

   • 可控性 Gramian
     $$W_c={\int }_0^{\infty }{e}^{A\tau }B{B}^T{e}^{A^T\tau }d\tau ,$$
     解 Lyapunov 方程
     $$A{W}_c+W_c{A}^T+B{B}^T=0.$$
   • 可观性 Gramian
     $$W_o={\int }_0^{\infty }{e}^{A^T\tau }{C}^{T}C{e}^{A\tau }d\tau ,$$
     解 Lyapunov 方程
     $$A^T{W}_o+W_{o}A+C^{T}C=0.$$

2. 平衡变换的构造

   • 对 $W_c$ 和 $W_o$ 做 Cholesky 分解：
     $$W_c=R{R}^T,W_o=S{S}^T.$$

   • 对 $R^{T}S$ 做 SVD：
     $$\begin{gathered} R^{T}S=U\Sigma V^T, \\ \Sigma =\mathrm{diag}({\sigma }_1,\dots ,{\sigma }_n), \end{gathered}$$
     { σ i } \{\sigma_i\} {σi​} 即 Hankel 奇异值，按降序排列。

   • 于是取
     $$\begin{gathered} T={\Sigma }^{-\frac{1}{2}}{U}^T{S}^T, \\ {T}^{-1}=RV{\Sigma }^{-\frac{1}{2}}, \end{gathered}$$

     能使得在新坐标下，$$\begin{gathered} {\tilde{W}}_c=T{W}_c{T}^T=\Sigma \\ {\tilde{W}}_o=T^{-T}{W}_o{T}^{-1}=\Sigma \end{gathered}$$

3. 坐标变换
   定义
   $$\bar{x}=Tx⟹x=T^{-1}\bar{x}.$$
   将原系统变换得
   $$\begin{gathered} \dot{\bar{x}}=T\dot{x}=T(Ax+Bu)=\left(TA{T}^{-1}\right)\bar{x}+\left(TB\right)u, \\ y=Cx+Du=\left(C{T}^{-1}\right)\bar{x}+Du. \end{gathered}$$
   设
   $$A_b=TA{T}^{-1},B_b=TB,C_b=C{T}^{-1},$$
   则平衡化系统为
   $$\begin{gathered} \dot{\bar{x}}=A_b\bar{x}+B_{b}u, \\ y=C_b\bar{x}+Du. \end{gathered}$$

   此时可控性和可观性 Gramian 均为 $\Sigma$，对角元素 ${\sigma }_i$ 刚好是各状态的能量度量。

4. 截断降阶
   根据 $\Sigma =\mathrm{diag}({\sigma }_1,\dots ,{\sigma }_n)$，保留能量大的前 $r$ 个分量，舍弃后 $n-r$ 个分量。将
   $$\begin{gathered} \bar{x}=\left[\begin{array}{c}{\bar{x}}_1\\ {\bar{x}}_2\end{array}\right],{\bar{x}}_1\in {\mathbb{R}}^r,{\bar{x}}_2\in {\mathbb{R}}^{n-r} \\ {A}_b=\left[\begin{array}{cc}{A}_{11}& A_{12}\\ A_{21}& A_{22}\end{array}\right], \\ {B}_b=\left[\begin{array}{c}{B}_1\\ B_2\end{array}\right], \\ {C}_b=\left[\begin{array}{cc}{C}_1& C_2\end{array}\right], \end{gathered}$$

丢弃 ${\bar{x}}_2$ 及其关联块，令 $x_r\equiv {\bar{x}}_1$，得到约简模型：

$$\begin{gathered} {\dot{x}}_r=A_{11}{x}_r+B_{1}u, \\ y=C_1{x}_r+Du. \end{gathered}$$

这样，整个过程严格按照“$\bar{x}=Tx$”的坐标变换来推导，变换后直接在平衡坐标下截断，就得到了降阶模型。

当对原系统做平衡截断降阶后，原来的初始状态 $x(0)$ 以及 对应的初始输出 $y(0)=Cx(0)$ 都必须 映射到降阶系统的坐标空间，否则就无法直接在低维系统中使用它们来启动仿真或分析。

• 降阶模型的初始状态：

  • 给定原始系统的初始状态 $x(0)$，降阶模型的初始状态应取
    $$x_r(0)={\bar{x}}_1(0)=[Tx(0){]}_{1:r},$$

    即先用 $T$ 映射到平衡坐标，再截取前 $r$ 分量。

  • 这样，降阶系统从正确的低维初始条件开始演化，才能近似再现原系统的初始响应。

• 原系统的初始输出：
  • 原始系统在 $t=0$ 的输出是 $y(0)=Cx(0)+Du(0)$，若已知 $D=0$，且若假设初始输入 $u(0)=0$，则初始输出简化为 $y(0)=Cx(0)$；
  • 降阶系统的输出方程 截断后，用 $(A_r,B_r,C_r)$ 表示低维模型，则
    $$y(t)\approx C_r{x}_r(t)+Du(t).$$
  • 将 降阶初始状态 按上面方式设置为 $x_r(0)=[Tx(0){]}_{1:r}$，则
    $$C_r{x}_r(0)=C_1{\bar{x}}_1(0)=C_1[Tx(0){]}_{1:r}\approx Cx(0),$$
    其中 $C_1$ 是对 $CT$ 的前 $r$ 列截取 。
  • 这样，降阶系统的初始输出就能尽量贴合原系统。

这样就解决了“维度不匹配”的问题，确保降阶模型能够从正确的低维初始条件开始，重现原系统的初始响应。

求解 HSV 为什么不使用定义而是使用 Cholesy 和SVD 分解？

在数值实现中，很少把 Hankel 奇异值（HSV）直接定义为可控 Gramian $W_c$ 与可观 Gramian $W_o$ 乘积的特征值平方根来计算，而是先做 Cholesky 分解再做 SVD，其主要原因有以下几点：

1. 避免显式构造乘积矩阵，降低计算成本与内存开销

   • 如果直接计算 $W_c{W}_o$，首先要把两个 $n\times n$ 的 Gramian 显式地存储并相乘，生成一个新的 $n\times n$ 矩阵，所需的存储量和乘法运算均为 $O(n^3)$ 量级，且中间结果通常比原来更“密集”甚至更难存储。

   • 而通过 Cholesky 分解：

     $$W_c=R{R}^T,W_o=S{S}^T,$$

     只需分别存储和操作 $R$ 与 $S$（同样是三角矩阵），然后对 $S^{T}R$ 做一次 SVD，就能得到奇异值 $\Sigma$，即 HSV，无需生成 $W_c{W}_o$ 本身，从而节省了额外的存储和矩阵乘法开销。

2. 提高数值稳定性

   • Gramian 通常是病态的：它的特征值会很快衰减（尤其是在高阶系统中），导致 $W_c$ 和 $W_o$ 的条件数都非常大。若再相乘得到 $W_c{W}_o$，其条件数大约是原来条件数的平方，数值误差成倍放大，可能导致特征值计算崩溃。
   • 相反，Cholesky 分解对正定矩阵是数值稳定且无需列主元（pivoting）的操作，而对 $S^{T}R$ 做 SVD（而非对不对称的 $W_c{W}_o$ 做特征分解）可以直接获取奇异值，且 SVD 本身对误差更不敏感，能够在高病态情况下仍然给出可靠的奇异值排序和数值结果。

3. 天然支持低秩/截断优化

   • 在大规模系统中，Gramian 往往是低秩或数值低秩的。Cholesky 分解（或其它“平方根”算法）可以直接求出低秩因子 $R,S$，保留有效维度 $r\ll n$，随后只对 $S^{T}R\in {\mathbb{R}}^{r\times r}$ 做 SVD，大幅降低运算量。
   • 这在“平方根算法”（square-root method）中被广泛采纳，也是 MATLAB balreal 在大规模场景中常用的实现策略 。

总结：

• 定义法（直接特征值分解 $W_c{W}_o$）在实现上既费内存，又数值极不稳定；
• Cholesky+SVD 法（先分解再奇异值分解）既避免了构造病态乘积矩阵，也利用了数值线性代数中对称正定矩阵和奇异值分解的稳定性优势，因而成为计算 HSV 的标准做法。

MATLAB 实践

MATLAB 平衡截断 步骤：MATLAB -> 模型降阶器 -> 导入模型（如工作区中的状态空间模型 sys_siso = ss(A, B_siso, C_siso, D_siso) 变量）-> 选中模型后平衡截断 -> 输入目标阶数。

Python 实现

Python-Control 0.9.4 文档总览 — 无 balreal 定义

不过有现成的轮子 pyMor，使用 pip install pymor 安装。

注意，之前由于安装了一个 pip install slycot，导致一运行就 报错 AttributeError: module 'slycot' has no attribute 'sb03md57'，通常是因为 Python 中安装的 Slycot 包缺少底层 Fortran/C 扩展，导致无法导入 SLICOT 的新接口 sb03md57。

• 解决方法：卸载掉即可，pip uninstall slycot 和 pip uninstall pymor，然后再重新安装 pip install pymor。
• 可能这就是 Python 的好处带来的弊端吧，开源的同时带来了很多不规范。

import scipy.sparse as sps
from pymor.models.iosys import LTIModel
from pymor.reductors.bt import BTReductorA, B, C, D, x0, y0 = ...  # 自定义syso = LTIModel.from_matrices(A, B, C, D)
bt = BTReductor(fom)
sysb = bt.reduce(20)A_r = sysb.A.matrix   # Reduced A matrix, shape (r, r)
B_r = sysb.B.matrix   # Reduced B matrix, shape (r, m)
C_r = sysb.C.matrix   # Reduced C matrix, shape (p, r)
D_r = sysb.D.matrix   # Reduced D matrix, shape (p, m)Ti = bt.W   # shape (n, r)
Ti_mat = Ti.to_numpy()  # shape (r, n)
x0_r = Ti_mat.dot(x0)  # shape (r,)
y0_r = C_r.dot(x0_r)plot_bode_and_error(syso, sysb)

def plot_bode_and_error(fom, rom, w=None, save_prefix=''):# 1) Compute H-infinity error bounds for all ordersbt = BTReductor(fom)bounds = bt.error_bounds()  # Returns list of error bounds for orders 1,2,… :contentReference[oaicite:0]{index=0}# Select the bound corresponding to the reduced orderr = rom.ordererr_bound = bounds[r-1] if r <= len(bounds) else None# 2) Prepare default frequency range if not specifiedif w is None:w = (1e-5, 1e6)# 3) Plot Bode magnitude & phase for full and reduced modelsfig, axs = plt.subplots(2, 1, figsize=(8, 6), sharex=True)fig.suptitle('Bode Plot Comparison', fontsize=14)# Use pyMOR's built-in transfer_function bode_plot methods :contentReference[oaicite:1]{index=1}fom.transfer_function.bode_plot(w, ax=axs, label='Full-order')rom.transfer_function.bode_plot(w, ax=axs, linestyle='--', label='Reduced-order')# Finalize axesfor ax in axs:ax.grid(True)ax.legend(loc='best')axs[1].set_xlabel('Frequency (rad/s)')fig.tight_layout(rect=[0, 0, 1, 0.95])fig.savefig(f'{save_prefix}bode_comparison.png', dpi=300)plt.close(fig)# 4) Report the H-infinity error boundif err_bound is not None:print(f"Estimated H∞ error bound for reduced order {r}: {err_bound:.2e}")  # :contentReference[oaicite:2]{index=2}else:print("No error bound available for order", r)return err_bound

先验知识：可控性 Gramian $W_c$、可观性 Gramian $W_o$ 以及 Hankel 奇异值（HSV）${\sigma }_i$

• 可控性 Gramian $W_c$ 衡量 系统从零初始状态经输入到达某一状态所需的能量，反映了状态可控程度；
• 可观性 Gramian $W_o$ 衡量 系统从某一状态到零输出所需的能量，反映了状态可观测程度。
• Hankel 奇异值 ${\sigma }_i$ 定义为 可控性 Gramian 与可观性 Gramian 乘积 $W_c{W}_o$ 的特征值的平方根，即 ${\sigma }_i=\sqrt{{\lambda }_i{W}_c{W}_o},i=1,2,\dots ,n$，用以 量化每个状态在输入–输出能量传递中的贡献大小。

可控性 Gramian $W_c$

在 MATLAB Control Toolbox 中，可使用命令

Wc = gram(sys,'c');

分别计算可控性 Gramian 或其 Cholesky 因子 (gram - MathWorks)。

对于 连续时间 LTI 系统，

$$\dot{x}(t)=Ax(t)+Bu(t),$$

其无限时域可控性 Gramian 定义为

$$W_c={\int }_0^{\infty }{e}^{A\tau }B{B}^T{e}^{A^T\tau }d\tau ,$$

并且它是 Lyapunov 方程

$$A{W}_c+W_c{A}^T+B{B}^T=0$$

的唯一正定解 (Controllability Gramian)。

对于 离散时间 LTI 系统，

$$x[k+1]=Ax[k]+Bu[k],$$

可控性 Gramian 定义为

$$W_{c,d}=\sum _{m=0}^{\infty }{A}^{m}B{B}^T(A^T{)}^m,$$

它满足离散 Lyapunov 方程

$$W_{c,d}-A{W}_{c,d}{A}^T=B{B}^T,$$

且当 $A$ 稳定（谱半径 <1）时，$W_{c,d}$ 为正定矩阵。

可观性 Gramian $W_o$

对于 同一连续系统附带输出方程，

$$y(t)=Cx(t)+Du(t),$$

无限时域可观性 Gramian 定义为

$$W_o={\int }_0^{\infty }{e}^{A^T\tau }{C}^{T}C{e}^{A\tau }d\tau ,$$

它是 Lyapunov 方程

$$A^T{W}_o+W_{o}A+C^{T}C=0$$

的唯一正定解 (Observability Gramian)。

对于 离散系统，

$$y[k]=Cx[k]+Du[k],$$

可观性 Gramian 定义为

$$W_{o,d}=\sum _{m=0}^{\infty }(A^T{)}^m{C}^{T}C{A}^m,$$

满足离散 Lyapunov 方程

$$W_{o,d}-A^T{W}_{o,d}A=C^{T}C,$$

当 $A$ 稳定时，$W_{o,d}$ 为正定。

Hankel 奇异值 ${\sigma }_i$

Hankel 奇异值 { σ i } i = 1 n \{\sigma_i\}_{i=1}^n {σi​}i=1n​ 定义为矩阵 $W_c{W}_o$ 的特征值 { λ i } \{\lambda_i\} {λi​} 的平方根：

$${\sigma }_i=\sqrt{{\lambda }_i\left(W_c{W}_o\right)},i=1,\dots ,n.$$

它们可视为 系统在输入–输出能量传递路径上的“能量谱”，按降序排列时，较大的 ${\sigma }_i$ 对模型行为贡献更大。

在 MATLAB 中，通过平衡变换函数 balreal 可直接得到 $W_c$ 与 $W_o$ 对角化后的奇异值序列，即 HSV 向量。