通俗理解线性与非线性、时变与时不变系统，和数值不稳定性机制

在工程中，系统 可看作“把输入信号变成输出信号的黑盒子”。

• 时不变：无论你什么时候给它输入，表现都一样；
• 时变：它随着“时间”本身而改变；
• 线性：对输入做加权叠加，输出也是相应的加权叠加；
• 非线性：输出和输入之间有“怪异”或“曲折”的关系，不能直接拆分成简单之和。

系统	输入–输出关系	线性？	时变？
$y(t)=2x(t)$	纯比例放大	线性	时不变
$y(t)=x^2(t)$	平方运算	非线性	时不变
$y(t)=tx(t)$	随时间线性调整增益	线性	时变
$y(t)=\sin [x(t)]$	正弦非线性转换	非线性	时不变

控制系统

一、时不变 vs. 时变

把输入信号整体向右（或左）平移，再送入系统；如果输出也严格同样平移 → 时不变，否则 → 时变。

给定状态空间模型， $$\begin{gathered} \dot{x}(t)=A(t)x(t)+B(t)u(t) \\ y(t)=C(t)x(t)+D(t)u(t) \end{gathered}$$

真正决定 时变/时不变 的是 $A(t)$、$B(t)$、$C(t)$、$D(t)$ 是否随 $t$ 变化，

• 若 $A(t)$ 随时间变化，系统在不同 $t$ 的“响应法则”就不同，输入平移后输出并不能简单平移，则为时变系统；
• 若 $A$ 恒为常数，即与时间无关，系统对任意时刻的“对状态增益”一致，不会因时刻不同而改变，则为时不变系统。

时不变系统（Time-Invariant）

时不变（Time-Invariant）系统的定义是：对任意输入 $u(t)$ 和任意时移 $\tau$，若系统满足

S { u ( t − τ ) } = y ( t − τ ) , S\{u(t-\tau)\} \;=\; y(t-\tau), S{u(t−τ)}=y(t−τ),

则该系统为时不变系统。

• “你在早上给它一杯水，它做 A；下午给它同样一杯水，它还是做 A。”
• 换句话说，系统的“规则”不会因为时间流逝而改变。

时变系统（Time-Varying）

• “你早上给它一杯水，它做 A；下午给它同样一杯水，它却做 C。”
• 系统的“规则”会随时刻不同而发生改变。

二、线性 vs. 非线性

取任意两个输入 $x_1$、$x_2$，验证

$$S(a{x}_1+b{x}_2)\stackrel{?}{=}aS(x_1)+bS(x_2).$$
若总是相等 → 线性（满足“叠加性 + 齐次性”），否则 → 非线性（对输入做平方、绝对值、饱和裁剪等操作后，不再是简单倍数或叠加）。

线性叠加原理

一个系统 S { ⋅ } S\{\cdot\} S{⋅} 若对任意两组输入 $(x_1,u_1)$、$(x_2,u_2)$ 及任意常数 $a,b$ 都满足

S { a ( x 1 , u 1 ) + b ( x 2 , u 2 ) } = a S { ( x 1 , u 1 ) } + b S { ( x 2 , u 2 ) } , S\{\,a\,(x_1,u_1) + b\,(x_2,u_2)\} = a\,S\{(x_1,u_1)\} \;+\; b\,S\{(x_2,u_2)\}, S{a(x1​,u1​)+b(x2​,u2​)}=aS{(x1​,u1​)}+bS{(x2​,u2​)},

则称为线性系统。

状态方程可视为映射，

$$f(x,u)=Ax+Bu.$$

若对任意 $x_1,x_2$ 与 $u_1,u_2$、常数 $a,b$ 有

$$\begin{gathered} f(a{x}_1+b{x}_2,a{u}_1+b{u}_2)=A(a{x}_1+b{x}_2)+B(a{u}_1+b{u}_2) \\ =a(A{x}_1+B{u}_1)+b(A{x}_2+B{u}_2)=af(x_1,u_1)+bf(x_2,u_2), \end{gathered}$$

则 $f$ 为线性映射，系统即为线性。

数值不稳定性机制

Moler & Van Loan 指出，算术抵消尤其会在乘方次数较多的区域（the “hump” region）发生，其后续的误差可能永远无法补偿，造成 最终结果严重偏离真实值。

• Moler & Van Loan 在其“十九种矩阵指数算法”一文中强调，只有当中间矩阵 条件数较小 或进行充分平衡（balancing）后，重复平方才可能稳定，否则极易触发上述抵消问题。
• Higham 在后续研究中也给出了精确的舍入误差边界，证明若不进行 缩放或平衡处理，重复平方往往会超出可接受的数值误差范围。

最近遇到的问题：我在求解 $\dot{x}=Ax+Bu$ 时使用了数值解，
$$\mathbf{x}(t)=e^{\mathbf{A}t}{\mathbf{x}}_0+{\mathbf{A}}^{-1}\left(e^{\mathbf{A}t}-\mathbf{I}\right)\mathbf{B}\mathbf{u}$$
但是遇到了 同样的求解代码 在使用 float32 和 float64 下得到的 结果相差两个数量级，但 float64 的结果看上去更合理，于是，了解到下面的内容。

• 条件数放大效应：对于高条件数矩阵 $R_k$，小的特征值分量在相对误差放大下，易被舍入至零 或产生巨大扭曲。
• 重复累积：每一步平方不仅有新的舍入误差，还要处理上一步的误差后结果，若不加控制，误差将以指数形式累积。

上述原理共同导致了在 float32 精度下，对条件数大或“hump”区域的矩阵进行多次平方时，输出结果可能与理论值产生戏剧性偏离，从而影响整个解析解求解过程。

直观理解：

1. 几何视角：矩阵作用可视为把单位球拉伸成椭球，高条件数意味着椭球的主轴长短极端不一。误差项相当于在“拉伸后的空间”中加上一层厚度，当椭球某一最短轴（对应极小特征值）的半径小于这层厚度时，就会被“覆盖”看不到。
2. 频域比喻：将矩阵分解到特征模态后，某些高频或低频分量信号极弱（幅度很小），而 舍入误差像是一片噪声地毯，当信号振幅低于噪声地毯高度时，信号就彻底淹没，计算结果只剩下噪声部分。
3. 算法应对：关键在于降低条件数或限制每一步运算的数值跨度，如对矩阵平衡（balancing）、缩放-平方、分段迭代等技术，都是在每步计算中尽量避免极端模态与中间误差相互抵消。

条件数（condition number）

What is the Condition Number of a Matrix?

条件数（condition number）衡量数值问题——包括标量函数、矩阵运算等——对输入微小扰动的敏感程度。

• 若条件数很大，表示输出会被输入误差显著放大，是病态（ill‐conditioned）问题；
• 若条件数接近 1，则是良态（well‐conditioned）问题。

对矩阵 $A$ 而言，其条件数定义为

$${\kappa }_p(A)=\Vert A{\Vert }_p\Vert A^{-1}{\Vert }_p,$$

等价于其奇异值比值 ${\kappa }_2(A)={\sigma }_{\max }(A)/{\sigma }_{\min }(A)$。

1. $\Vert A{\Vert }_p$——矩阵的 $p$ 范数（诱导范数）
   $$\Vert A{\Vert }_p=\underset{x\ne 0}{\sup }\frac{\Vert Ax{\Vert }_p}{\Vert x{\Vert }_p}=\underset{\Vert x{\Vert }_p=1}{\max }\Vert Ax{\Vert }_p$$
   这里 $\Vert \cdot {\Vert }_p$ 对向量是常见的 $p$-范数（如 $p=2$ 时为欧几里得范数），而对矩阵则取该范数在输入和输出空间上的算子（诱导）范数。

• 直观：想象单位 $p$ 范数球 V p , n = { x : ∥ x ∥ p ≤ 1 } V_{p,n}=\{x:\|x\|_p\le1\} Vp,n​={x:∥x∥p​≤1} 被线性映射 $A$ 后变形为凸集 $A(V_{p,n})$，$\Vert A{\Vert }_p$ 就是这个凸集在某方向上的最大“伸长”距离。

2. $\Vert A^{-1}{\Vert }_p$——矩阵逆的 $p$ 范数
   同样是对逆映射 $A^{-1}$ 的诱导范数
   $$\Vert A^{-1}{\Vert }_p=\underset{y\ne 0}{\sup }\frac{\Vert A^{-1}y{\Vert }_p}{\Vert y{\Vert }_p}=\underset{\Vert y{\Vert }_p=1}{\max }\Vert A^{-1}y{\Vert }_p$$

• 直观：度量如果输出发生单位量级扰动 $(\Vert y{\Vert }_p=1)$，对应的输入（原始 $x=A^{-1}y$）会被放大到多大。

3. ${\kappa }_p(A)$——条件数
   $${\kappa }_p(A)=\Vert A{\Vert }_p\Vert A^{-1}{\Vert }_p=\underset{x\ne 0}{\max }\frac{\Vert Ax{\Vert }_p}{\Vert x{\Vert }_p}\times \underset{y\ne 0}{\max }\frac{\Vert A^{-1}y{\Vert }_p}{\Vert y{\Vert }_p}.$$

• 含义：它衡量了解线性系统 $Ax=b$ 时，右端 $b$ 上的相对微小扰动会在解 $x$ 上放大多少倍。
  例如，当 ${\kappa }_p(A)=10^6$ 时，$b$ 上 1e-6 的相对误差可能被放大到解上 1 的相对误差。

1. 原理剖析：条件数如何放大抵消

在 高条件数 矩阵运算中，“算术抵消”（catastrophic cancellation）是由于 小模态（极小特征值方向）的幅度远小于大模态，导致 舍入误差完全掩盖了真正的信号部分。

当一个矩阵 $R_k$ 的 条件数 $\kappa (R_k)=\Vert R_k\Vert \Vert R_k^{-1}\Vert$ 很大 时，它在某些方向上会 “极度伸缩”——有的特征向量被放大到极大模态，有的被压缩到极小模态。

1.1 舍入误差与矩阵模

1. 乘积误差：单次乘法 $r_{ik}\times r_{kj}$ 在浮点下的相对误差为 ${\epsilon }_{\mathrm{mach}}$, 绝对误差约为 ${\epsilon }_{\mathrm{mach}}|r_{ik}{r}_{kj}|$。
2. 范数级放大：取最大元素或算子范数时，所有 $|r_{ik}{r}_{kj}|$ 均被上界为 $\Vert R_k{\Vert }^2$，故绝对误差上界为 ${\epsilon }_{\mathrm{mach}}\Vert R_k{\Vert }^2$。
3. 累加累积：再经过 $n-1$ 次加法的舍入，误差进一步被 ${\gamma }_n\approx n{\epsilon }_{\mathrm{mach}}$ 放大，但仍与 $\Vert R_k{\Vert }^2$ 成正比。

因此，矩阵乘法的总舍入误差自然表现为 ${\epsilon }_{\mathrm{mach}}\Vert R_k{\Vert }^2$。

在浮点运算中，矩阵乘法 运算中产生的 舍入误差 规模大约是 ${ϵ}_{\mathrm{mach}}\Vert R_k{\Vert }^2$，

$$\widehat{R_k^2}=R_k^2+E_k,\Vert E_k\Vert \lesssim c{ϵ}_{\mathrm{mach}}\Vert R_k{\Vert }^2,$$

其中 ${ϵ}_{\mathrm{mach}}$ 是机器精度（float32 约 $10^{-7}$，）。当 $\Vert R_k\Vert$ 很大时，$\Vert E_k\Vert$ 也相应放大。

对于 IEEE 754 标准下的浮点数表示：

• float32（单精度） 的机器精度（machine epsilon，$ϵ\_\mathrm{mach}$）约为$${ϵ}_{\mathrm{mach}}^{\text{float32}}\approx 1.19\times 10^{-7}$$
• float64（双精度） 的机器精度约为
  $${ϵ}_{\mathrm{mach}}^{\text{float64}}\approx 2.22\times 10^{-16}$$
  $ϵ\_\mathrm{mach}$ 定义为满足：
  $$1+{ϵ}_{\mathrm{mach}}>1$$
  的最小浮点数，即浮点数系统中能与 1 区分开的最小量，反映了浮点数系统的相对精度。

如果你用 Python 的 numpy 验证可以这样写：

import numpy as np
print(np.finfo(np.float32).eps)  # 输出: 1.1920929e-07
print(np.finfo(np.float64).eps)  # 输出: 2.220446049250313e-16

1.2 极小模态与抵消

设 $R_k$ 在某个方向上的真实平方值分量为 $\alpha \ll \Vert R_k{\Vert }^2$。若舍入误差与极小模态同量级甚至更大，

$$\Vert E_k\Vert \approx {ϵ}_{\mathrm{mach}}\Vert R_k{\Vert }^2\gg \alpha ,$$

则计算结果中该分量变为

$$\hat{\alpha }=\alpha +\delta ,|\delta |\approx \Vert E_k\Vert \gg \alpha ,$$

因此小模态 $\hat{\alpha }$ 会被舍入成与 $\delta$ 同级或直接被舍入误差吞没为零——典型的算术抵消，也就是我们看不到原本就微弱的那部分信号（相对误差剧增）。

2. 数值例子：2×2 矩阵的抵消演示

考虑对角矩阵

$$R=\left(\begin{array}{cc}10^4& 0\\ 0& 10^{-4}\end{array}\right),$$
其条件数 $\kappa (R)=10^4/10^{-4}=10^8$ 非常大，属于高度病态矩阵。用 float32 模拟一次“自乘”运算——对应缩放-平方中的一次平方步骤。

import numpy as np# 原始矩阵对角元
r = np.array([1e4, 1e-4], dtype=np.float32)# 模拟 R^2 的对角运算
r_sq = r * r
print("r_sq =", r_sq)  # 期望 [1e8, 1e-8]

运行结果（float32）：

r_sq = [1.000000e+08 9.999999e-09]

• 第一项正确，但第二项本该是 $10^{-8}$，却被舍入为 9.999999e-09（注意虽未完全归零，但相对误差已高达 100%），有时甚至会被舍入成 0，这取决于具体硬件与中间误差累积程度。

3. 对角缩放（平衡）恢复小模态

为了均衡数值量级，我们 做相似变换：

$$\begin{gathered} D=\mathrm{diag}(10^2,10^{-2}), \\ \tilde{R}=D^{-1}RD=\mathrm{diag}(10^2,10^{-2}). \end{gathered}$$

此时 $\kappa (\tilde{R})=10^4$ 大幅降低，重复上述运算：

r_bal = np.array([1e2, 1e-2], dtype=np.float32)
r_bal_sq = r_bal * r_bal
print("r_bal_sq =", r_bal_sq)  # 期望 [1e4, 1e-4]

运行结果（float32）：

r_bal_sq = [1.e+04 1.e-04]

两项均能准确表示，避免了算术抵消。最后通过逆变换恢复到原尺度。

缩放策略降低条件数 原理及公式推导

• 对矩阵或数值同时左、右乘（或除以）对角因子 $D$，可在 不改变本征结构（相似变换） 的前提下，将数值量级压缩到精度允许的安全区间。

• 缩放后，显著降低矩阵条件数，从而避免重复平方过程中的舍入误差累积与抵消。

总结：大条件数矩阵的极端伸缩作用，使得极小模态在浮点舍入误差面前显得微不足道，从而在“减法”或“取差”时被完全抵消。对策则是让矩阵在计算时“穿戴护甲”——即 缩放，使得 每一步运算都在安全的数值范围内进行。

回到我的问题，我需要对矩阵 $A$ 进行缩放得到变换矩阵 $T$，然后对系统从 $\dot{x}=Ax+Bu,y=Cx+Du+Const$ 进行坐标变换到 $z$-空间，计算完结果后再变换回原坐标 $x$-空间。
关键点：相似变换不改变输入–输出响应

• 虽然对 $A$ 做了“缩放”得到 $A_b=T^{-1}AT$，但只要同时对状态 $x$、输入矩阵 $B$ 及输出矩阵 $C$ 做配套变换，系统的 传递函数 和 时域响应 完全不变。
• 对 $z$-空间使用解析解 或数值积分后，转换回 $x$-空间 即可得到与原系统完全一致的结果。

对状态空间模型

$$\dot{x}=Ax+Bu,y=Cx+Du+\mathrm{Const}$$

施行 相似变换

$$x=Tz$$

并 不会改变系统的输入–输出特性（物理响应），仅是将状态表征从原坐标系 $x$ 换到平衡后坐标系 $z$。为保持等价性，所有与状态相关的矩阵和向量——包括 $A$、$B$、$C$、初始状态 $x_0$、输出常数 $\mathrm{Const}$——必须同步变换，而输入矩阵 $D$ 不受影响。

这保证了新坐标系下的系统与原系统在输入 $u$ 相同情况下，输出 $y$ 完全相同——只不过是 在不同的“状态基”下描述同一个动力学系统。

下面给出完整的变换原理与公式，演示从平衡变换到仿真再到结果还原的全流程。

相似变换原理与公式：

令

$$x=Tz,T\in {\mathbb{R}}^{n\times n}\text{可逆},$$

则

$$\dot{z}=T^{-1}\dot{x}=T^{-1}(Ax+Bu)=(T^{-1}AT)z+(T^{-1}B)u.$$

输出方程

$$y=Cx+Du+\mathrm{Const}=CTz+Du+\mathrm{Const}.$$

因此 缩放后模型为

$$\begin{gathered} \dot{z}=A_{b}z+B_{b}u, \\ y=C_{b}z+Du+{\mathrm{Const}}_b, \end{gathered}$$

其中

$$\begin{gathered} A_b=T^{-1}AT, \\ {B}_b=T^{-1}B, \\ {C}_b=CT, \\ {\mathrm{Const}}_b=\mathrm{Const}. \end{gathered}$$

注意：$D$ 和常数项 $\mathrm{Const}$ 保持不变，因为它们不与状态乘以 $T$ 相关。

初始状态变换：

若原初始状态为 $x(0)=x_0$，则对应新坐标为

$$z(0)=T^{-1}{x}_0.$$

仿真或解析解后得 $z(t)$，再通过

$$x(t)=Tz(t)$$

还原至原坐标。

带常数 offset 的处理：

• 常数项 $\mathrm{Const}$：由于输出方程是仿射关系，$x$ 变换不影响 $\mathrm{Const}$，故 ${\mathrm{Const}}_b=\mathrm{Const}$.
• 若想把常数项也纳入状态，可用 增广法（augmented state）把常数视为额外输入通道，但在纯相似变换下无需改变常数维度。

总结：

1. 相似变换 $x=Tz$ 只改变状态坐标，不改变系统动态与输出。
2. 所有与状态相关 的矩阵 $A,B,C$ 及初始状态 $x_0$ 必须同步变换。
3. 输出常数项 $\mathrm{Const}$ 与 $D$ 都保持不变，除非采用增广状态方法才需要另行处理。
4. 仿真后通过 $x=Tz$ 将结果还原，即可与原模型输出严格匹配。

Is variable scaling essential when solving some PDE problems numerically?

在半导体模拟中，通常会将方程缩放以使其具有归一化值。例如，在极端情况下，半导体中的电子密度可以变化超过 18 个数量级，电场可以变化很大，超过 6 个（或更多）数量级。然而，论文中从未真正解释为什么要这样做。

• 如果考虑双精度中保留的位数，你会看到如果“双精度下是否有足够的位数来应对这些波动”确实成立…真正的问题在于当将这些数字输入数值算法时：取平方，突然就有 36 个数量级的差异…
• 解决（线性）偏微分方程包括将方程离散化以得到线性系统，该系统由线性求解器求解，其收敛性（速率）取决于矩阵的条件数。变量缩放通常可以降低条件数，从而提高收敛性。（这基本上相当于应用对角预处理器，参见 Nicholas Higham 的《数值算法的精度与稳定性》。）

刚性（Stiffness）与快慢模态

系统“刚性”（stiffness）常对应于雅可比矩阵的特征值跨度极大，对于线性常系数系统 $\dot{x}=Ax$，若 $A$ 的特征值 { λ i } \{\lambda_i\} {λi​} 满足
$$\frac{\max |\Re (\lambda )|}{\min |\Re (\lambda )|}\gg 1$$，则该比值称为 stiffness ratio，常与矩阵的条件数同义（在正规矩阵情形下，两者同量级）。

对称或正规矩阵 时，$\kappa (A)=\frac{|{\lambda }_{\text{max}}|}{|{\lambda }_{\text{min}}|}$ 正好等于 stiffness ratio。

这导致显式积分器对步长 $dt$ 有严格的稳定性限制，必须取非常小的 $dt$ 才能避免数值发散。将 $dt$ 取得很小，实际上相当于将 $Adt$ 的谱半径压缩到显式方法的绝对稳定区域内，从而保证稳定性。

在刚性（stiff）系统中，状态包含 “快模态”（快速衰减分量）和“慢模态”（缓慢演化分量），二者时间尺度相差很大。

• 显式积分法（如显式 Euler 或 Runge‑Kutta）只有在步长 $h$ 足够小以将所有特征值乘积 ${\lambda }_{i}h$ 落在其绝对稳定域内时才能收敛，否则就会发生 “数值发散”——也就是解的振荡或爆炸，看似 “失踪”了快速分量的正确衰减。
• 换言之，步长过大时，显式方法 无法“捕捉”快模态的衰减行为，从而在迭代中不断放大对应模态的误差，最终导致解发散。

总结：

• 刚性大 → 特征值跨距大 → 快模态时间常数极短 → 必需 $h\ll {\tau }_{\mathrm{fast}}$ 才将所有 ${\lambda }_{i}h$ 保证在稳定域。
• 步长过大 → 快模态 ${\lambda }_{\mathrm{fast}}h$ 落在显式方法的发散区间 → 模态误差被指数放大 → 数值解发散。
• 理解核心：并非“捕捉不到”仅是丢失精度，而是该模态对应的稳定函数 $\varphi (z)$ 在该 $z$ 区域本身就是发散的多项式运算，必然导致解的爆炸性错误。

扩展思考

平衡截断 与 条件数

平衡截断（Balanced Truncation, BT）是一种 模型降阶 方法，其 目标是最小化输入–输出误差的 ${\mathcal{H}}_{\infty }$ 或 ${\mathcal{H}}_2$ 范数，而不是专门为降低矩阵 $A$ 的算子范数条件数 $\kappa (A)$ 而设计。

Balanced Truncation 的核心思想是对可控性 Gramian $P$ 与可观测性 Gramian $Q$ 做相似变换使之相等并对角化，然后截断那些对应 Hankel 奇异值（${\sigma }_i$）较小的态，得到低阶模型。此过程保证了
$$\Vert G-G_r{\Vert }_{\infty }\le 2\sum _{i=r+1}^n{\sigma }_i$$
的误差上界，却不保证或专门优化原始矩阵 $A$ 的条件数 $\kappa (A)$。

• BT 并非为降低 $\kappa (A)$ 而生，BT 关注的是 系统输入–输出传递函数的逼近误差，与矩阵指数或线性求逆算子的数值病态性没有直接关联。
• BT 可能在转换后产生一个具有更差$\kappa (A_b)$ 的平衡矩阵 $A_b$，尤其当系统含不稳定模态或 Gramian 求解本身数值不稳定时。

矩阵条件数 $\kappa (A)=\Vert A\Vert \Vert A^{-1}\Vert$ 衡量的是前向算子对输入扰动的放大倍数，与模型的阶数无必然对应关系：高阶系统也可以有低条件数，低阶系统也可能高条件数。

对于求解线性方程或矩阵运算，条件数

$$\kappa (A)=\Vert A\Vert \Vert A^{-1}\Vert$$

直接决定前向误差最多会被放大多少倍—— 若 $\kappa (A)=10^k$，则 最终误差最多可损失 $k$ 位精度。高条件数会导致小的舍入误差被放大，使得双精度（float64）也不能完全消除数值偏差。

舍入误差源于病态，而非高阶：数值舍入误差 主要取决于所用算法的 数值稳定性 与矩阵的 病态程度（condition number），而非模型阶数本身。要提高数值精度，应针对 降低 $\kappa (A)$ 或选用更稳定的求解方法，而非单纯降低阶数。

• 阶数 $n$ 与 条件数 $\kappa (A)$ 并不一一对应。一个 $n=172$ 的系统如果是对角矩阵且对角元素幅度齐一，依然可以有 $\kappa (A)=1$，精度问题几乎不存在。
• 相反，一个小阶系统（例如 $n=10$）如果极度病态（如 Hilbert 矩阵），$\kappa (A)\approx 10^{12}$，同样会遭遇严重舍入抵消问题。
• 因此，数值精度的关键在于条件数小，而非阶数低。

• 如果目标是消除数值误差或更好地求解刚性系统，应采用矩阵等比缩放（equilibration）或刚性求解器，而模型降阶只是对 I/O 行为逼近的一种手段，不能保证消除数值病态。

结论：

• 降阶可减小模型的维度和在线计算成本，但不必然降低条件数。
• 数值精度依赖于 $\kappa (A)$ 与所用算法的稳定性，而非模型的阶数。
• 若需要提高数值可靠性，请重点：
  1. 评估并压缩 $\kappa (A)$（equilibration、预条件）；
  2. 选用 对刚性友好的数值方法（隐式积分、Krylov expm）；
  3. 在极端情况下，提升浮点精度。
• 仅当 BT 后显著减少 Hankel 奇异值总和且关心的是I/O逼近误差时，才以降阶为首选；否则，应直接面向条件数与数值算法下手。

矩阵等比缩放（Equilibration / Ruiz 算法）

矩阵等比缩放（Equilibration），也称Ruiz 算法，是一种 迭代预处理方法，旨在通过左右对角缩放将矩阵 $A\in {\mathbb{R}}^{n\times n}$ 的每一行和每一列在某一范数（通常是 ${\ell }_2$ 或 ${\ell }_{\infty }$ 范数）下的长度逼近相同，从而显著降低条件数 $\kappa (A)$，提升后续数值计算的稳定性和精度。该方法源自 Daniel Ruiz 等人提出的交替缩放思想，并在大规模稀疏系统、Krylov 子空间方法等场景中被广泛应用。下文先给出原理与公式推导，再提供完整的 Python 实现示例。

目标：最小化条件数

给定矩阵 $A\in {\mathbb{R}}^{n\times n}$，希望找到行缩放矩阵

$$D_r=\mathrm{diag}(d_1,\dots ,d_n)$$

和列缩放矩阵

$$D_c=\mathrm{diag}(e_1,\dots ,e_n)$$

使得缩放后矩阵

$$A^{\prime }=D_{r}A{D}_c$$

的条件数 $\kappa (A^{\prime })=\Vert A^{\prime }\Vert \Vert A^{\prime -1}\Vert$ 较原矩阵显著降低。直接优化此目标是非凸且计算昂贵的，因此采用 迭代启发式方案：通过 轮流缩放行与列，使得 每行／列的 $p$-范数逐步接近 。

迭代等比缩放（Ruiz 算法）

对于范数参数 $p\ge 1$（常用 $p=2$ 或 $p=\infty$），定义第 $k$ 次迭代的矩阵为 $A^{(k)}$，初始化

$$D_r^{(0)}=I,D_c^{(0)}=I,A^{(0)}=A.$$

然后重复以下两步，直到收敛或达到最大迭代次数：

1. 行等比缩放

   • 计算每行范数

     $$r_i^{(k)}=\Vert (A^{(k)}{)}_{i,:}{\Vert }_p,i=1,\dots ,n.$$

   • 构造行缩放因子

     $$s_i^{(k)}=\frac{1}{\sqrt{r_i^{(k)}}}(\text{若 }{r}_i^{(k)}=0\text{ 则 }{s}_i^{(k)}\leftarrow 1).$$

   • 更新

     $$\begin{gathered} D_r^{(k+1)}=\mathrm{diag}\left(s_i^{(k)}\right)D_r^{(k)}, \\ {A}^{(k+\frac{1}{2})}=\mathrm{diag}\left(s_i^{(k)}\right)A^{(k)}. \end{gathered}$$

2. 列等比缩放

   • 计算每列范数

     $$c_j^{(k)}=\Vert (A^{(k+\frac{1}{2})}{)}_{:,j}{\Vert }_p,j=1,\dots ,n.$$

   • 构造列缩放因子

     $$t_j^{(k)}=\frac{1}{\sqrt{c_j^{(k)}}}(\text{若 }{c}_j^{(k)}=0\text{ 则 }{t}_j^{(k)}\leftarrow 1).$$

   • 更新

     $$\begin{gathered} D_c^{(k+1)}=D_c^{(k)}\mathrm{diag}\left(t_j^{(k)}\right), \\ {A}^{(k+1)}=A^{(k+\frac{1}{2})}\mathrm{diag}\left(t_j^{(k)}\right). \end{gathered}$$

迭代结束后，得到

$$\begin{gathered} A_{\mathrm{eq}}=D_r^{(K)}A{D}_c^{(K)}, \\ {D}_r=D_r^{(K)},D_c=D_c^{(K)}, \end{gathered}$$

并可证明行列范数均接近 1，$\kappa (A_{\mathrm{eq}})$ 通常比 $\kappa (A)$ 小数个数量级。

通过 Ruiz 算法得到了两组对角缩放因子

$$D_r=\mathrm{diag}(d_i),D_c=\mathrm{diag}(e_j)$$

缩放后的矩阵 $A^{\prime }=D_{r}A{D}_c$

坐标变换：原状态 $x$ 与等比缩放后的状态 $z$ 的关系 为

$$z=D_c^{-1}x,x=D_{c}z$$
这里 $D_c^{-1}=\mathrm{diag}(1/e_j)$。

在 Ruiz 平衡（equilibration）中，我们对矩阵 $A$ 同时左右两端做对角缩放
$$A_{\mathrm{eq}}=D_{r}A{D}_c,$$
其中

• $D_c=\mathrm{diag}(c_1,\dots ,c_n)$ 对应“列缩放”，
• $D_r=\mathrm{diag}(r_1,\dots ,r_n)$ 对应“行缩放”。

但在动力系统的状态空间变换里，我们只做“列缩放”——也就是只用 $D_c$——其原理如下：

1. 列缩放 $D_c$
   • 它相当于对每个状态分量 $x_j$ 做缩放：
     $$x=D_{c}z⟺z=D_c^{-1}x.$$
   • 变换到新坐标 $z$ 后，系统矩阵变为
     $$\dot{z}=D_c^{-1}A{D}_{c}z+D_c^{-1}Bu.$$
   • 这正是我们在求解、仿真或优化时经常使用的“状态变量坐标变换”。
2. 行缩放 $D_r$
   • 它相当于对方程的“输出残差”或“度量”做缩放：
     $$y=Cx⟹\tilde{y}=D_{r}y=(D_{r}C)x.$$
   • 在平衡 $A$ 的条件数时，左乘 $D_r$ 能让每一行（即每个状态方程）的系数范数也趋于统一、平衡；它并不改变状态的坐标，只是等价地改变了“我们看向方程”的度量。

• 状态坐标变换 的本质是：
  $$x=D_{c}z⟹\dot{x}=D_c\dot{z}.$$
  把它带入原系统
  $$D_c\dot{z}=A(D_{c}z)+Bu⟹\dot{z}=D_c^{-1}A{D}_{c}z+D_c^{-1}Bu.$$
  这里根本没有出现 $D_r$，因为我们并没有去重新定义“输出”或“方程的度量”，而只是换了一套“状态坐标”。
• 行缩放 $D_r$ 虽然也出现在平衡算法里，但它对应的是“把方程的左侧（$\dot{x}$）或者输出 $y$ 做度量变换”，不属于“状态坐标变化”的一部分。它对状态轨迹本身并不做变换，只是等价地把不同方程的量纲／大小平衡了。

python 实现

import numpy as np
from scipy.linalg import normdef ruiz_equilibrate(A, n_iter=10, p=2):"""Ruiz 等比缩放：迭代对 A 做行/列等比缩放，使每行/列 p-范数趋于1。Returns: A_eq, Dr, Dc"""n = A.shape[0]Dr = np.ones(n)Dc = np.ones(n)A_eq = A.copy().astype(float)for _ in range(n_iter):# 行缩放row_norms = np.linalg.norm(A_eq, ord=p, axis=1)row_norms[row_norms == 0] = 1s = 1.0 / np.sqrt(row_norms)Dr *= sA_eq = (s[:, None] * A_eq)# 列缩放col_norms = np.linalg.norm(A_eq, ord=p, axis=0)col_norms[col_norms == 0] = 1t = 1.0 / np.sqrt(col_norms)Dc *= tA_eq = (A_eq * t[None, :])return A_eq, Dr, Dc# 示例：构造病态矩阵
n = 172
diag = np.logspace(0, 8, n)
A = np.diag(diag) + 1e-3 * np.ones((n,n))# 预处理
A_eq, Dr, Dc = ruiz_equilibrate(A, n_iter=20, p=2)print("cond before:", np.linalg.cond(A, p=2))
print("cond after: ", np.linalg.cond(A_eq, p=2))# 状态映射示例
x0 = np.random.randn(n,1)
Dc_inv = 1.0 / Dc# 压缩到平衡坐标
z0 = Dc_inv[:, None] * x0# 在平衡坐标下仿真（示例：x'=A_eq x）
z1 = A_eq @ z0# 重构回原始坐标
x1 = Dc[:, None] * z1

系统状态矩阵条件数降低技术研究

模型预处理阶段

• 状态变量等比缩放：对系统状态向量的各分量按合适比例重新标度（即对矩阵 $A$ 做 $A^{\prime }=D^{-1}AD$ 的相似变换）。这种对角缩放相当于施加对角预处理器，常能显著降低矩阵的条件数。
  • 例如，MATLAB 的 equilibrate 函数即对 $A$ 的行列进行缩放，使对角线元素幅值为1，从而改善条件并提高数值稳定性。
• 矩阵平衡 (Matrix Balancing)：类似于对角缩放的一种相似变换，通过 调整矩阵的行范数和列范数使其更加均衡，进而改善条件数。矩阵平衡技术可显著降低非对称矩阵的条件数，从而加速迭代求解收敛并提升数值稳定性。
• 平衡实现 / 平衡状态空间：对线性系统进行平衡变换，使 可控性 Gramian 和可观测性 Gramian 相等并对角化。这种变换平衡了系统的输入-输出能量分布，可用于后续的平衡截断或模型降维，通常可减少状态幅值的不平衡，从而有助于改善数值特性。

数值求解阶段

• 矩阵指数和逆运算：直接计算矩阵指数 $\exp (At)$ 时，推荐使用高质量的数值算法（如 SciPy/ MATLAB 的 expm），其基于 Pade 近似和缩放-平方技术。注意，当 $A$ 条件数极高时，

  • $\exp (At)$ 的计算仍会受敏感误差影响；
  • 尤其是显式求逆 $A^{-1}$ 的操作几乎是病态问题，误差难以控制。因此应 避免显式求逆，必要时通过解线性方程组代替或直接计算积分项。例如，可 利用扩展矩阵指数方法或块矩阵方法隐式求解，无需孤立求 $A^{-1}$。

• 数值积分替代：对于高阶 LTI 系统，使用 数值积分器（如 Runge-Kutta 或隐式步进法）求解 $\dot{x}=Ax+Bu$ 可以 避免直接计算矩阵指数，从而提高稳定性。

  • 对于非刚性系统，显式 Runge-Kutta（如 RK4/RK45）通常足够；
  • 但若系统存在 刚性（特征值实部范围很宽），应采用 隐式方法（如 Radau IIA、BDF 等）。推荐使用 MATLAB 的 ode15s、SciPy 的 solve_ivp（可选用 BDF/Radau 模式）等，并合理调整时间步长和误差公差，以兼顾精度和效率。

• 矩阵指数作用计算：若 只需计算 $\exp (A)t$ 的作用，可使用 Krylov 子空间方法（如 SciPy 的 expm_multiply）直接对向量作用，这种方法在计算大矩阵指数时往往更高效且更稳定。该函数允许预处理矩阵（例如通过提供迹 traceA 来做预处理），从而改善迭代性能。

• 多精度和预处理：若双精度仍不足以保证精度，可以尝试 多精度算术（如 Python 的 mpmath、MATLAB 的变量精度计算 vpa 等），以 减少舍入误差。但是多精度计算成本极高，仅建议在非常关键的精度验证时使用。对于 线性方程求解，可使用预条件器（如对称对角预条件、Incomplete LU 等）改善迭代收敛，但对直接计算矩阵指数的帮助有限。