从MATLAB仿真到理论证明：实战验证能控性与格拉姆矩阵的关系-编程实验室

从MATLAB仿真到理论证明：实战验证能控性与格拉姆矩阵的关系

在控制系统的设计与分析中，能控性是一个基础而重要的概念。它决定了我们能否通过合适的控制输入，将系统从任意初始状态驱动到期望的终态。格拉姆矩阵判据作为判断能控性的有力工具，其理论证明往往让初学者感到抽象难懂。本文将采用一种独特的"代码复现理论"方法，通过MATLAB仿真直观展示格拉姆矩阵与能控性的关系，再回溯到严格的数学证明，为工程师和学生提供一条从实践到理论的理解路径。

1. 实验环境搭建与系统建模

1.1 MATLAB基础配置

在开始之前，确保你的MATLAB环境已经安装了Control System Toolbox。这个工具箱提供了我们需要的所有控制相关函数。可以通过以下命令检查：

ver('control')

如果未安装，可以通过MATLAB的Add-Ons菜单进行安装。接下来，我们创建一个简单的脚本文件gramian_controllability.m作为我们的实验基础。

1.2 构建线性时不变系统

我们选择一个典型的二阶系统作为研究对象，其状态空间表示为：

A = [0 1; -2 -3]; % 状态矩阵 B = [0; 1]; % 输入矩阵 C = [1 0]; % 输出矩阵（本实验暂不使用） D = 0; % 直接传输矩阵 sys = ss(A, B, C, D); % 创建状态空间模型

这个系统代表了一个简单的质量-弹簧-阻尼系统，其中：

状态x₁表示位置
状态x₂表示速度
输入u表示外力

2. 格拉姆矩阵的计算与能控性分析

2.1 格拉姆矩阵的理论基础

对于线性时不变系统，能控性格拉姆矩阵定义为：

W_c(t) = ∫₀ᵗ e^{Aτ}BBᵀe^{Aᵀτ} dτ

其中e^{Aτ}是状态转移矩阵。MATLAB提供了直接计算格拉姆矩阵的函数gram：

t = 5; % 选择的时间区间 Wc = gram(sys, 'c', t); % 计算能控性格拉姆矩阵

2.2 格拉姆矩阵的数值计算验证

为了深入理解计算过程，我们可以手动实现格拉姆矩阵的计算：

% 定义时间步长 dt = 0.01; time = 0:dt:t; % 初始化格拉姆矩阵 Wc_manual = zeros(size(A)); % 数值积分计算 for tau = time expA = expm(A*tau); Wc_manual = Wc_manual + expA*B*B'*expA'*dt; end

比较两种方法的结果差异：

error = norm(Wc - Wc_manual); disp(['计算误差：', num2str(error)]);

2.3 能控性判据的数值验证

格拉姆矩阵判据指出：系统完全能控当且仅当存在t>0使得W_c(t)非奇异。我们可以通过计算行列式来验证：

det_Wc = det(Wc); if abs(det_Wc) > 1e-6 % 考虑数值误差 disp('系统是完全能控的'); else disp('系统不是完全能控的'); end

同时，MATLAB还提供了直接的能控性检查函数：

Co = ctrb(sys); % 计算能控性矩阵 rank_Co = rank(Co); if rank_Co == size(A,1) disp('能控性矩阵满秩，系统能控'); else disp('系统不能控'); end

3. 参数变化对能控性的影响

3.1 改变系统参数

让我们修改系统矩阵A，观察能控性的变化：

A_new = [0 1; -1 -1]; % 修改后的系统矩阵 sys_new = ss(A_new, B, C, D);

重新计算格拉姆矩阵和能控性矩阵：

Wc_new = gram(sys_new, 'c', t); Co_new = ctrb(sys_new);

3.2 能控性丢失的情况

考虑一个明显不能控的系统：

A_unctrl = [1 0; 0 2]; B_unctrl = [1; 0]; % 只能控制第一个状态 sys_unctrl = ss(A_unctrl, B_unctrl, C, D);

计算其能控性格拉姆矩阵：

Wc_unctrl = gram(sys_unctrl, 'c', t); disp(['行列式：', num2str(det(Wc_unctrl))]);

3.3 可视化分析

我们可以绘制格拉姆矩阵行列式随时间变化的曲线：

time_points = linspace(0.1, 10, 50); det_values = zeros(size(time_points)); for i = 1:length(time_points) Wc_temp = gram(sys, 'c', time_points(i)); det_values(i) = det(Wc_temp); end plot(time_points, det_values); xlabel('时间 t'); ylabel('det(Wc(t))'); title('格拉姆矩阵行列式随时间变化'); grid on;

4. 从仿真回溯到理论证明

4.1 充分性证明的数值解释

仿真结果显示，当格拉姆矩阵非奇异时，我们确实能够找到控制输入将系统从任意状态驱动到原点。这对应着理论证明中的构造性部分：

u(τ) = -Bᵀe^{Aᵀ(T-τ)}W_c(T)^{-1}x₀

在MATLAB中，我们可以实现这个控制律：

x0 = [1; -1]; % 任意初始状态 T = 5; % 计算控制输入 Wc_T = gram(sys, 'c', T); tau = linspace(0, T, 100); u = zeros(size(tau)); for i = 1:length(tau) u(i) = -B'*expm(A'*(T-tau(i)))*inv(Wc_T)*x0; end % 仿真系统响应 lsim(sys, u, tau, x0);

4.2 必要性证明的反例验证

对于不能控的系统，如前面定义的sys_unctrl，其格拉姆矩阵始终是奇异的，这与理论预测一致。我们可以通过特征值分析来理解：

[V, D] = eig(A_unctrl); disp('特征向量：'); disp(V); disp('能控性矩阵的秩：'); disp(rank(ctrb(A_unctrl, B_unctrl)));

结果显示存在与输入不耦合的模式，验证了理论中的反证法逻辑。

4.3 格拉姆矩阵与能控性子空间

格拉姆矩阵的秩揭示了能控性子空间的维度：

rank_Wc = rank(Wc); disp(['格拉姆矩阵的秩：', num2str(rank_Wc)]);

这与能控性矩阵的秩理论结果一致，提供了另一种验证途径。

5. 高级应用与扩展思考

5.1 时变系统的扩展

虽然本文聚焦线性时不变系统，但格拉姆矩阵概念可以推广到时变系统：

W_c(t) = ∫₀ᵗ Φ(t,τ)B(τ)B(τ)ᵀΦ(t,τ)ᵀ dτ

其中Φ(t,τ)是状态转移矩阵。MATLAB中可以通过数值积分实现。

5.2 数值稳定性考虑

在实际计算中，需要注意数值稳定性问题：

% 更稳健的行列式计算 log_det = 2*sum(log(diag(chol(Wc)))); disp(['对数行列式：', num2str(log_det)]);

5.3 能控性格拉姆矩阵的其他应用

格拉姆矩阵不仅能判断能控性，还可用于：

最优控制中的最小能量控制计算
模型降阶中的能控性度量
系统辨识中的输入设计

例如，最小能量控制可以通过：

x0 = [1; -1]; xf = [0; 0]; T = 5; Wc_T = gram(sys, 'c', T); u_opt = @(t) -B'*expm(A'*(T-t))*inv(Wc_T)*(expm(A*T)*x0 - xf); t_opt = linspace(0, T, 100); u_opt_values = arrayfun(u_opt, t_opt);