MATLAB roots函数实战：5分钟搞定高阶系统稳定性判断（附完整代码）

MATLAB roots函数实战：高阶系统稳定性分析的黄金法则

在控制工程和自动化领域，系统稳定性分析是每个工程师的必修课。面对复杂的高阶系统特征方程，传统的手工计算方法不仅耗时耗力，还容易出错。而MATLAB的roots函数配合简单的可视化技巧，能在几分钟内完成从特征根求解到稳定性判断的全过程。本文将带你掌握这套高效工作流的核心要点，即使你是MATLAB新手也能快速上手。

1. 理解特征根与系统稳定性的关系

任何动态系统的稳定性都可以通过其特征方程的根（即特征根）来判断。对于连续时间系统，稳定性判据非常简单：所有特征根的实部均为负数时系统稳定；只要有一个特征根的实部为正数，系统就不稳定。

举个例子，考虑一个六阶系统的特征方程：

s⁶ + 3s⁵ + 16s⁴ + 2s³ - 4s² - 25s - 60 = 0

手工求解这个方程的根几乎是不可能的任务，但MATLAB可以轻松解决。

特征根在复平面上的分布直观反映了系统动态特性：左半平面代表衰减，右半平面代表发散，虚轴代表等幅振荡。

2. 使用roots函数求解特征根的完整流程

2.1 准备系数向量

MATLAB中多项式用系数向量表示，从最高次项到常数项依次排列。对于上面的六阶方程，构建系数向量的代码如下：

p = [1, 3, 16, 2, -4, -25, -60]; % s⁶系数为1，s⁵系数为3，...，常数项-60

如果多项式有缺项，必须用0补位。例如，对于方程s⁴ + 2s² + 1 = 0，对应的系数向量应为：

p = [1, 0, 2, 0, 1]; % 补上s³和s的0系数

2.2 调用roots函数求解

有了系数向量后，求解特征根只需一行代码：

r = roots(p)

执行后会返回一个包含所有特征根的复数向量。例如，对于我们的六阶系统，结果可能类似：

r = -1.4204 + 3.6565i -1.4204 - 3.6565i 1.4475 + 0.0000i -0.1563 + 1.4343i -0.1563 - 1.4343i -1.2941 + 0.0000i

2.3 快速稳定性判断

检查特征根实部是否有正值的最直接方法：

if any(real(r) > 0) disp('系统不稳定！存在右半平面极点') else disp('系统稳定') end

3. 特征根分布可视化技巧

图形化展示能让稳定性判断更加直观。下面这段代码可以生成专业的特征根分布图：

% 绘制特征根分布 figure h = plot(real(r), imag(r), 'rx', 'MarkerSize', 12, 'LineWidth', 2); hold on % 绘制坐标轴 plot([0 0], ylim, 'k-', 'LineWidth', 1) % 虚轴 plot(xlim, [0 0], 'k-', 'LineWidth', 1) % 实轴 % 美化图形 grid on xlabel('实部') ylabel('虚部') title('系统特征根分布图') set(gca, 'FontSize', 12)

这段代码会生成一个带坐标轴和网格的散点图，其中：

• 每个"x"标记代表一个特征根
• 右半平面（实部>0）的根表示系统不稳定
• 虚轴附近的根可能表示系统有较差的动态性能

4. 高级应用与实用技巧

4.1 批量处理多个特征方程

当需要分析多个系统的稳定性时，可以编写函数自动化处理：

function stability = checkStability(p) r = roots(p); stability = all(real(r) < 0); % 可视化 figure plot(real(r), imag(r), 'rx', 'MarkerSize', 12) hold on plot([0 0], ylim, 'k-') plot(xlim, [0 0], 'k-') grid on title(['系统稳定性: ' num2str(stability)]) end

调用示例：

p1 = [1, 3, 16, 2, -4, -25, -60]; checkStability(p1);

4.2 性能优化与数值稳定性

对于极高阶系统（如100阶以上），roots函数可能出现数值不稳定的情况。此时可以考虑：

1. 使用balance函数预处理系数矩阵
2. 转换为状态空间形式后使用eig函数
3. 对于稀疏系统，考虑专门的稀疏矩阵求解器

% 高阶系统处理示例 p = rand(1,101); % 随机生成100阶多项式 A = compan(p); % 构建伴随矩阵 r = eig(A); % 使用特征值求解

4.3 实际工程中的注意事项

• 系数精度问题：工程中获得的系数往往有测量误差，可以通过蒙特卡洛分析评估稳定性鲁棒性
• 临界稳定判断：理论上实部刚好为零的根在实际中极少出现，应设置合理的阈值（如abs(real(r)) < 1e-6）
• 复根共轭对验证：物理系统的特征根总是成共轭对出现，可以添加验证代码：

[~,idx] = sort(imag(r)); % 按虚部排序 r_sorted = r(idx); for i = 1:2:length(r_sorted) if abs(r_sorted(i) - conj(r_sorted(i+1))) > 1e-6 warning('发现非共轭复根，请检查系统模型有效性') end end

5. 完整案例：从理论到实践

让我们通过一个完整的案例来巩固所学内容。假设有一个四阶控制系统，其特征方程为：

s⁴ + 5s³ + 10s² + 10s + 4 = 0

完整的MATLAB分析代码如下：

% 案例：四阶系统稳定性分析 p = [1, 5, 10, 10, 4]; % 系数向量 % 求解特征根 r = roots(p); disp('特征根为：') disp(r) % 稳定性判断 if any(real(r) > 0) disp('系统不稳定') else disp('系统稳定') end % 可视化 figure plot(real(r), imag(r), 'rx', 'MarkerSize', 15, 'LineWidth', 2) hold on plot([0 0], ylim, 'k-', 'LineWidth', 1.5) plot(xlim, [0 0], 'k-', 'LineWidth', 1.5) grid on xlabel('实部', 'FontSize', 12) ylabel('虚部', 'FontSize', 12) title('四阶系统特征根分布', 'FontSize', 14) set(gca, 'FontSize', 11) % 添加根值标签 for i = 1:length(r) text(real(r(i)), imag(r(i)), sprintf(' %.2f%+.2fi', real(r(i)), imag(r(i))), ... 'VerticalAlignment','bottom', 'FontSize', 10) end

执行这段代码后，MATLAB会输出特征根并显示分布图。从图中可以直观看到所有根都位于左半平面，因此系统是稳定的。