零极点到状态空间的直接高数值鲁棒性转换方法问询

问题背景 #

给定零极点增益模型$(z,p,k)$，需将其转换为具备高数值精度与稳定性的状态空间表示$(A,B,C,D)$，且完全规避传递函数作为中间环节。常规方法通过poly(z)、poly(p)生成分子分母多项式再转换为状态空间，但在极点分布特殊场景（如Wilkinson示例中$p=(1,2,3,…,19,20)$）下，双精度计算无法准确表示多项式系数，导致数值精度损失。

已知可构造$A=\text{diag}(p)$（极点为$A$的特征值），严格真系统下$D=0$，核心问题是直接由零点确定矩阵$B$和$C$。

直接构造方法 #

利用部分分式展开的数值稳定特性，直接基于极点、零点的两两差值计算$B$和$C$，无需生成多项式：

1. 确定矩阵$A$和$D$

   • $A$为极点构成的对角矩阵：$A = \text{diag}(p_1, p_2, ..., p_n)$，其中$n$为极点数
   • 严格真系统（零点数≤极点数）下，$D=0$

2. 构造向量$B$
   选择任意非零$n$维列向量即可，为简化计算通常取全1向量：

   B = ones(n, 1);
   

3. 计算向量$C$
   $C$的每个元素$c_i$对应极点$p_i$的部分分式留数，直接通过极点与零点的差值计算：
   $$c_i = k \cdot \frac{\prod_{j=1}^m (p_i - z_j)}{\prod_{\substack{m=1 \ m \neq i}}^n (p_i - p_m)}$$
   其中$m$为零点数。

   在代码中可通过逐元素计算实现（以Matlab为例）：

   n = length(p);
   m = length(z);
   C = zeros(1, n);
   for i = 1:n
       % 计算分子：k乘以(p_i与所有零点的差值乘积)
       numerator = k;
       for j = 1:m
           numerator = numerator * (p(i) - z(j));
       end
       % 计算分母：p_i与所有其他极点的差值乘积
       denominator = 1;
       for idx = 1:n
           if idx ~= i
               denominator = denominator * (p(i) - p(idx));
           end
       end
       C(i) = numerator / denominator;
   end
   

方法优势 #

• 完全规避多项式展开步骤，避免了poly函数带来的数值精度损失
• 所有计算基于极点、零点的两两差值，在极点分布密集或存在相近极点的场景下，数值稳定性远优于传递函数中转法
• 直接对应部分分式展开形式，物理意义清晰，且易于验证传递函数匹配性

内容的提问来源于stack exchange，提问作者DaveC