当前位置：首页 > news >正文

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

news 来源：原创 2025/7/3 4:55:15

平衡截断
- balreal 算法原理
- 平衡截断过程
- - 求解 HSV 为什么不使用定义而是使用 Cholesy 和SVD 分解？
- MATLAB 实践
- Python 实现
先验知识：可控性 Gramian $W_c$ 、可观性 Gramian $W_o$ 以及 Hankel 奇异值（HSV） $\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$ 误差界 $\|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$ ，对角元素 $\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 的好处带来的弊端吧，开源的同时带来了很多不规范。

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_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 向量。

平衡截断

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​

相关文章：

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

可控性 Gramian $W_c$

可观性 Gramian $W_o$

Hankel 奇异值 $\sigma_i$