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

技术文档

平衡截断
- balreal 算法原理
- 平衡截断过程
- - 求解 HSV 为什么不使用定义而是使用 Cholesy 和SVD 分解？
- MATLAB 实践
- Python 实现
先验知识：可控性 Gramian $W_c$ 、可观性 Gramian $W_o$ 以及 Hankel 奇异值（HSV） $\\sigma_i$
- 可控性 Gramian $W_c$
- 可观性 Gramian $W_o$
- Hankel 奇异值 $\\sigma_i$

平衡截断

平衡截断（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$ 范数下满足经典误差界：
$\\|G - G_r\\|_\\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，以 提示截断时勿删除不稳定状态。

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

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

稳定性保证：Gramian 只有在系统稳定时才存在有限、正定的解；对不稳定部分直接调用 Lyapunov 会发散或得到非正定结果。

误差界仅对稳定部分：平衡截断的 $∞\\mathcal H_\\infty$ 误差界 $i\\|G - G_r\\|_\\infty \\le 2\\sum_{i>r} \\sigma_i$ 只对稳定子系统有意义。

保留不稳定行为：不稳定极点通常对系统行为至关重要，不能被截断。

平衡截断过程

下面用坐标变换
$\\bar x = T\\,x$

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

对于原系统，考虑连续时间 LTI 系统
$\\dot x = A\\,x + B\\,u, \\\\[5pt] y = C\\,x + D\\,u$
其中 $x\\in\\mathbb R^n$ 、 $u\\in\\mathbb R^m$ 、 $y\\in\\mathbb R^p$ 。

求解 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.$
平衡变换的构造
- 对 $W_c$ 和 $W_o$ 做 Cholesky 分解：
  $W_c = R\\,R^T,\\quad W_o = S\\,S^T.$
- 对 $R^T S$ 做 SVD：
  $R^T S = U\\,\\Sigma\\,V^T,\\\\[5pt] \\Sigma = \\mathrm{diag}(\\sigma_1,\\dots,\\sigma_n),$
  $\\{\\sigma_i\\}$ 即 Hankel 奇异值，按降序排列。
- 于是取
  $\\Sigma^{-\\tfrac12} U^T S^T,\\\\ \\qquad\\\\ T^{-1} = R\\,V\\,\\Sigma^{-\\tfrac12},$
  
  能使得在新坐标下， $\\tilde W_c = T\\,W_c\\,T^T = \\Sigma \\\\[5pt] \\tilde W_o = T^{-T}W_o\\,T^{-1} = \\Sigma$
坐标变换
定义
$\\bar x = T\\,x \\quad\\Longrightarrow\\quad x = T^{-1}\\,\\bar x.$
将原系统变换得
$\\dot{\\bar x} = T\\,\\dot x = T\\bigl(A\\,x + B\\,u\\bigr) = \\bigl(T\\,A\\,T^{-1}\\bigr)\\,\\bar x \\;+\\; \\bigl(T\\,B\\bigr)\\,u, \\\\[5pt] y = C\\,x + D\\,u = \\bigl(C\\,T^{-1}\\bigr)\\,\\bar x + D\\,u.$
设
$A_b = T\\,A\\,T^{-1},\\quad B_b = T\\,B,\\quad C_b = C\\,T^{-1},$
则平衡化系统为
$\\dot{\\bar x} = A_b\\,\\bar x + B_b\\,u,\\\\[5pt] y = C_b\\,\\bar x + D\\,u.$

此时可控性和可观性 Gramian 均为 $\\Sigma$ ，对角元素 $i\\sigma_i$ 刚好是各状态的能量度量。
截断降阶
根据 $\\Sigma=\\mathrm{diag}(\\sigma_1,\\dots,\\sigma_n)$ ，保留能量大的前 $r$ 个分量，舍弃后 $n - r$ 个分量。将
$\\bar x = \\begin{bmatrix}\\bar x_1 \\\\[4pt]\\bar x_2\\end{bmatrix},\\quad \\bar x_1\\in\\mathbb R^r,\\; \\bar x_2\\in\\mathbb R^{n-r} \\\\[5pt] A_b = \\begin{bmatrix}A_{11}&A_{12}\\\\[3pt]A_{21}&A_{22}\\end{bmatrix},\\;\\\\[5pt] B_b = \\begin{bmatrix}B_1\\\\[3pt]B_2\\end{bmatrix},\\;\\\\[5pt] C_b = \\begin{bmatrix}C_1 & C_2\\end{bmatrix},$

丢弃 $\\bar x_2$ 及其关联块，令 $x_r\\equiv\\bar x_1$ ，得到约简模型：

$\\dot x_r = A_{11}\\,x_r + B_1\\,u,\\\\ \\qquad\\\\ y = C_1\\,x_r + D\\,u.$

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

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

降阶模型的初始状态：
- 给定原始系统的初始状态 $x (0)$ ，降阶模型的初始状态应取
  $x_r(0) \\;=\\;\\bar x_1(0) \\;=\\;\\Bigl[T\\,x(0)\\Bigr]_{1:r},$
  
  即先用 $T$ 映射到平衡坐标，再截取前 $r$ 分量。
- 这样，降阶系统从正确的低维初始条件开始演化，才能近似再现原系统的初始响应。
原系统的初始输出：
- 原始系统在 $t = 0$ 的输出是 $C\\,x(0) + D\\,u(0)$ ，若已知 $D = 0$ ，且若假设初始输入 $u (0) = 0$ ，则初始输出简化为 $C\\,x(0)$ ；
- 降阶系统的输出方程 截断后，用 $A_r,B_r,C_r)$ 表示低维模型，则
  $\\approx C_r\\,x_r(t) + D\\,u(t).$
- 将 降阶初始状态 按上面方式设置为 $x_r(0)=\\bigl[T\\,x(0)\\bigr]_{1:r}$ ，则
  $C_r\\,x_r(0) = C_1\\,\\bar x_1(0) = C_1\\,\\bigl[T\\,x(0)\\bigr]_{1:r} \\approx C\\,x(0),$
  其中 $C_1$ 是对 $C\\,T$ 的前 $r$ 列截取。
- 这样，降阶系统的初始输出就能尽量贴合原系统。

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

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

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

避免显式构造乘积矩阵，降低计算成本与内存开销
- 如果直接计算 $W_cW_o$ ，首先要把两个 $n\\times n$ 的 Gramian 显式地存储并相乘，生成一个新的 $n\\times n$ 矩阵，所需的存储量和乘法运算均为 $O(n^3)$ 量级，且中间结果通常比原来更“密集”甚至更难存储。
- 而通过 Cholesky 分解：
  
  $W_c = R\\,R^T,\\quad W_o = S\\,S^T,$
  
  只需分别存储和操作 $R$ 与 $S$ （同样是三角矩阵），然后对 $S^T R$ 做一次 SVD，就能得到奇异值 $\\Sigma$ ，即 HSV，无需生成 $W_cW_o$ 本身，从而节省了额外的存储和矩阵乘法开销。
提高数值稳定性
- Gramian 通常是病态的：它的特征值会很快衰减（尤其是在高阶系统中），导致 $W_c$ 和 $W_o$ 的条件数都非常大。若再相乘得到 $W_cW_o$ ，其条件数大约是原来条件数的平方，数值误差成倍放大，可能导致特征值计算崩溃。
- 相反，Cholesky 分解对正定矩阵是数值稳定且无需列主元（pivoting）的操作，而对 $S^T R$ 做 SVD（而非对不对称的 $W_cW_o$ 做特征分解）可以直接获取奇异值，且 SVD 本身对误差更不敏感，能够在高病态情况下仍然给出可靠的奇异值排序和数值结果。
天然支持低秩/截断优化
- 在大规模系统中，Gramian 往往是低秩或数值低秩的。Cholesky 分解（或其它“平方根”算法）可以直接求出低秩因子 $R, S$ ，保留有效维度 $r\\ll n$ ，随后只对 $S^T R\\in\\mathbb R^{r\\times r}$ 做 SVD，大幅降低运算量。
- 这在“平方根算法”（square-root method）中被广泛采纳，也是 MATLAB balreal 在大规模场景中常用的实现策略。

总结：

定义法（直接特征值分解 $W_cW_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 的好处带来的弊端吧，开源的同时带来了很多不规范。

注意， $V,W\\in\\mathbb R^{n\\times r}$ 即为降维映射和重构映射基底，

把原 $n$ 阶 $x$ 映射到降阶空间 $x_r$ ： $\\displaystyle x_{r} = W^T\\,x, \\quad x_r \\in\\mathbb R^r$
重构到原空间： $\\displaystyle \\hat x = V\\,x_{r}, \\quad \\hat x \\in\\mathbb R^n$ ，计算原空间输出 $\\displaystyle \\hat y = C\\,\\hat x + D\\,u$ 。

注意：不同版本的 pyMOR 中，基底可能命名不同，看源码。

在 pyMOR 中，bt.V 和 bt.W 并不是以一个 NumPy “矩阵”直接存储为形状 $(n, r)$ 的二维数组，而是以长度为 $r$ 、每个向量长度为 $n$ 的 VectorArray 存储。

当调用 V = bt.V.to_numpy()，实际上是在把一个长度为 $r$ 的 VectorArray（其中每个“向量”都是维度 $n$ ）展平成一个常规的 NumPy 二维数组。默认行为是将 “向量数组” 的每一行对应一个 VectorArray 中的向量，因此得到的形状就是 $(r, n)$ 。

在平衡截断的理论推导里，习惯把基矩阵写成 $V,W\\in\\mathbb R^{n\\times r}$ ，所以再进行转置得到常规的 shape。

import scipy.sparse as spsfrom pymor.models.iosys import LTIModelfrom 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) V = bt.V.to_numpy() # (r, n)W = bt.W.to_numpy() # (r, n)V, W = V.T, W.T # (n, r)# Project initial statex0_r = W.T @ x0# One‐step simulate in reduced modelx1_r = A_r @ x0_r + B_r @ u0y1_r = C_r @ x1_r + D_r @ u0# Reconstruct back to full orderx1_full = V @ x1_ry1_full = C @ x1_full + D @ u0plot_bode_and_error(syso, sysb)

def plot_bode_and_error(fom, rom, w=None, save_prefix=\'\'): # 1) Compute H-infinity error bounds for all orders bt = 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 order r = rom.order err_bound = bounds[r-1] if r <= len(bounds) else None # 2) Prepare default frequency range if not specified if w is None: w = (1e-5, 1e6) # 3) Plot Bode magnitude & phase for full and reduced models fig, 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 axes for 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 bound if 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_cW_o$ 的特征值的平方根，即 $\\sigma_i = \\sqrt{\\lambda_i W_cW_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]=A\\,x[k]+B\\,u[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)=C\\,x(t)+D\\,u(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]=C\\,x[k]+D\\,u[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 奇异值 $\\{\\sigma_i\\}_{i=1}^n$ 定义为矩阵 $W_cW_o$ 的特征值 $\\{\\lambda_i\\}$ 的平方根：

$\\sigma_i \\;=\\;\\sqrt{\\lambda_i\\!\\bigl(W_c\\,W_o\\bigr)},\\quad i=1,\\dots,n.$

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

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

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

平衡截断

balreal 算法原理

平衡截断过程

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

MATLAB 实践

Python 实现

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

可控性 Gramian $W_c$

可观性 Gramian $W_o$

Hankel 奇异值 $\\sigma_i$

公告

DeepSeek全套部署资料免费下载

免费可商用字体批量下载

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

平衡截断

balreal 算法原理

平衡截断过程

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

MATLAB 实践

Python 实现

先验知识：可控性 Gramian W c W_c Wc​、可观性 Gramian W o W_o Wo​ 以及 Hankel 奇异值（HSV） σ i \\sigma_i σi​

可控性 Gramian W c W_c Wc​

可观性 Gramian W o W_o Wo​

Hankel 奇异值 σ i \\sigma_i σi​

相关问题

公告

DeepSeek全套部署资料免费下载

免费可商用字体批量下载

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

可控性 Gramian $W_c$

可观性 Gramian $W_o$

Hankel 奇异值 $\\sigma_i$