企业官网建设流程全解析

悬臂梁的振动抑制是结构动力学与振动控制里最经典的问题。把一根悬臂梁用有限元建模成高维状态空间方程，再设计 LQR 控制器，是主动振动控制里最完整、也最容易被低估的一套流程。不论你是研究生阶段接触结构主动控制，还是工作中要处理传感器支架、翼面或机械臂末端的振动问题，这套方法都能直接复用。这篇文章我会把悬臂梁有限元模型的建立过程、LQR 控制器设计的核心逻辑、以及可运行的 Matlab 代码全部拆开来讲，重点不是堆公式，而是讲明白为什么要这么做、代码每句话在干什么、调参时你会遇到哪些问题。

1. 悬臂梁 LQR 控制器：从力学问题到控制问题

1.1 为什么主动振动控制非要先做有限元建模

悬臂梁从物理本质上看是一个连续体，理论上有无限多个振动模态。真要在控制器里处置一个“无限维系统”，任何算法都会卡死，工程上必然要做离散化。有限元建模是其中最通用的一种手段，它把梁切成一连串单元，每个单元的位移场用形函数近似，最终得到一个有限自由度的二阶常微分方程组。这个过程的核心价值有两个：一是把偏微分方程变成计算机能算的代数问题，二是只要合理加密网格，模型的精度就能切近真实结构，这个收敛性比集中质量法和假设模态法都更可控。

我见过很多同学直接跳过模型细节，拿一本教材上的解析模态公式就开整 LQR。解析法不是不行，但碰到变截面梁、带集中质量、局部开孔这些情况，解析解基本写不出来。而有限元模型在同样代码框架下只要改几个单元参数，结构变化就能被吸收进局部矩阵，在工程场景里的通用性要高得多。

1.2 LQR 在结构控制里的定位

LQR 是线性二次型调节器，它在结构振动控制里的地位，有点像“最优控制入门第一课”。控制目标可以浓缩成一个性能指标：

J = ∫(xᵀQx + uᵀRu)dt

其中 x 是系统状态，u 是控制输入，Q 和 R 是加权矩阵。这个式子的物理含义很直白：Q 大，意味着你更在意状态偏差快速收敛；R 大，意味着你更舍不得用力去控制。LQR 的输出是一个状态反馈增益矩阵 K，控制律 u = -Kx。它的魅力在于，只要系统能控、权矩阵设置合理，闭环系统就必然稳定，而且这个稳定是有严格理论保证的，不像 PID 要靠经验试凑。

和极点配置相比，LQR 的另一个好处是不需要预先指定闭环极点位置。极点配置看着精确，但对模型误差很敏感，一旦能控性矩阵病态，配出来的极点根本无法落地。LQR 本质上是在控制代价和收敛速度之间做权衡，理念更符合工程直觉。

1.3 有限元离散到状态空间再到 LQR 的三步走

整个项目可以简化成三步。第一步，用有限元方法把悬臂梁离散成 M、K 矩阵，加上阻尼 C，得到二阶振动方程；第二步，把二阶方程改写成状态空间形式 ẋ = Ax + Bu，让 LQR 能直接使用；第三步，选取 Q、R 并求解 Riccati 方程得到反馈增益 K，做闭环仿真验证。

这三步里最容易出问题的是第二步。很多人把 A、B 矩阵直接交给 lqr 函数，却根本不清楚状态变量是什么、控制力作用在哪个自由度上。后面一旦仿真结果不对，也无从下手排查。所以我下面会先花足够篇幅把有限元建模和状态空间转换讲透，再放出完整代码，这样你拿到代码后才知道每一行到底在算什么东西。

提醒一句：控制设计里最常见的“仿真漂亮、实验翻车”情形，八成的原因出在模型环节。固定端约束没删干净、自由度顺序错乱、阻尼矩阵随手给一个数，这些都会让 LQR 设计失去实际意义。

2. 悬臂梁有限元模型建立的完整推导

2.1 单元选择与矩阵组装

对于细长的悬臂梁，欧拉-伯努利梁理论已经足够精确，它忽略了剪切变形和转动惯量，每个节点只保留两个自由度：横向位移 w 和转角 θ。于是一个两节点梁单元就有四个自由度，分别是 w₁、θ₁、w₂、θ₂。位移场用 Hermite 插值函数来构造，既能保证节点处位移连续，也能保证转角连续。

单元的一致质量阵和刚度阵可以直接写成：

Me = (ρAle / 420) * [156 22le 54 -13le; 22le 4le² 13le -3le²; 54 13le 156 -22le; -13le -3le² -22le 4le²]

Ke = (EI / le³) * [12 6le -12 6le; 6le 4le² -6le 2le²; -12 -6le 12 -6le; 6le 2le² -6le 4le²]

这里 ρ 是密度，A 是截面积，l e 是单元长度，E 是弹性模量，I 是截面惯性矩。我第一次做这个模型的时候想偷懒用集中质量阵，因为 M 变成对角阵，组装和计算都省事。但后来对比发现，集中质量阵会让高频模态频率偏高，尤其是做控制设计需要前几阶模态时，误差会影响 LQR 的权矩阵设计。所以这里的建议很明确：做 LQR 控制建模，老老实实用一致质量阵，不要省这点计算量。

组装全局矩阵是有限元里最容易出错的地方。全局自由度编号按“节点 1 的两个自由度、节点 2 的两个自由度……”依次排列，那么第 e 个单元涉及的全局自由度索引就是 [2e-1, 2e, 2e+1, 2e+2]。比如 8 个单元、9 个节点，总自由度就是 18，第 3 个单元的四个局部自由度对应全局的第 5、6、7、8 个自由度。代码里用一个 for 循环把这个映射写出来，就可以完成 M、K 的组装。

2.2 固定端约束的正确处理方式

悬臂梁的固定端意味着 w = 0、θ = 0，这两个自由度在振动过程中始终为零。处理方式是在全局矩阵里把这两行两列直接删掉。有人会觉得把对应行和列置零、对角置 1 更“统一”，但对 LQR 设计来说，残余的自由度会让状态空间里混入毫无物理意义的零特征值，有时候还会干扰控制器设计。我建议干脆删除，让系统矩阵更干净。

代码里的实现是：

fixed_dofs = [1, 2]; free_dofs = 3:ndof; Mff = Mg(free_dofs, free_dofs); Kff = Kg(free_dofs, free_dofs);

需要注意的是，如果固定端不在第一个节点，而是在梁中间某处支撑，那固定自由度的编号也要相应调整。只要记住“先组装、后删自由度”这个原则，就不会乱。

2.3 阻尼矩阵与模态验证

阻尼矩阵 C 不像 M、K 那样能从材料参数直接推导，它取决于结构内部的能量耗散机制。工程上最常用的模型是瑞利阻尼，也叫比例阻尼：

C = αM + βK

其中 α、β 是待定系数。任意一阶模态的阻尼比可以写成：

ξ = α / (2ω) + βω / 2

如果知道前两阶固有频率 ω₁、ω₂，并且想要它们对应的阻尼比都等于某个目标值（比如 0.01），就可以解一个二元一次方程组把 α、β 求出来。代码里用模态分析先算固有频率，再反解阻尼系数，比直接拍一个 α、β 更靠谱，也让你对模型本身的频率分布有直觉。

这里强烈建议做一次模态验证：用 eig(Kff, Mff) 算出前几阶频率，和欧拉-伯努利梁的解析解对比。第一阶固有频率的解析公式是：

f₁ = (β₁² / (2πL²)) * sqrt(EI / (ρA))

其中 β₁ ≈ 1.875。如果你的有限元网格加密之后，第一阶频率越来越接近这个值，说明建模链路没有大问题。这一步十几分钟就能做完，但能给你后续控制设计吃一颗定心丸。

2.4 把二阶振动方程写进状态空间

离散后的振动方程是二阶常微分方程组：

M q̈ + C q̇ + K q = B_f u

要让 LQR 能用，必须把它改写成状态空间标准形式。定义状态变量：

x = [q; q̇]

那么一阶状态方程为：

ẋ = [0, I; -M⁻¹K, -M⁻¹C] x + [0; M⁻¹B_f] u

这里 0 是零矩阵，I 是单位阵。状态维度是自由度数 N 的两倍。比如 8 单元悬臂梁，总自由度 18，删掉固定端两个自由度之后 N = 16，状态空间维度就是 32。这个规模在 Matlab 里毫无压力，哪怕加满 20 个单元也没问题。

输入矩阵 B_f 的物理意义很关键。它描述控制力作用在哪些自由度上。如果压电作动器贴在自由端，等效外力作用在自由端节点的平动自由度上，那么 B_f 里对应位置取 1，其余取 0。这种建模方式让你可以灵活模拟不同位置的作动器，也是在仿真里对比“作动器安装在哪更有效”的根基。

3. Matlab 代码实现：从零搭建仿真链路

3.1 参数设置与网格划分

代码最前面几行是材料和几何参数。我给的铝梁尺寸是长 1 米、宽 2 厘米、高 5 毫米，这个尺寸在实验室里很常见，算出来的前几阶固有频率也在很舒服的频段。你可以换成钢、复合材料，或者换成更长的细梁，只需要改动前六行即可。

单元数量建议从 8 开始。单元太少，模型过于刚硬，频率偏差大；单元太多，状态维度暴涨，LQR 增益容易出现数值上的高频放大。实际用 8 到 20 个单元之间就可以，本文代码默认 8 个。

3.2 完整代码与逐步说明

下面是一个可以直接运行的完整脚本，从有限元建模到 LQR 设计再到闭环仿真一条龙。代码里的注释我写得很密，建议拿到手之后一段一段对照执行。

clear; clc; close all; %% 1. 结构与材料参数 L = 1.0; % 悬臂梁长度 (m) ne = 8; % 单元数量 b = 0.02; % 梁宽 (m) h = 0.005; % 梁高 (m) A = b * h; % 截面积 (m^2) I = b * h^3 / 12; % 截面惯性矩 (m^4) E = 70e9; % 弹性模量 (Pa) rho = 2700; % 密度 (kg/m^3) le = L / ne; % 单元长度 nnode = ne + 1; % 节点数量 ndof = 2 * nnode; % 总自由度，每节点含位移和转角 %% 2. 单元一致质量阵与刚度阵 Me = (rho*A*le/420) * [156 22*le 54 -13*le; 22*le 4*le^2 13*le -3*le^2; 54 13*le 156 -22*le; -13*le -3*le^2 -22*le 4*le^2]; Ke = (E*I/le^3) * [12 6*le -12 6*le; 6*le 4*le^2 -6*le 2*le^2; -12 -6*le 12 -6*le; 6*le 2*le^2 -6*le 4*le^2]; %% 3. 组装全局 M、K 矩阵 Mg = zeros(ndof); Kg = zeros(ndof); for e = 1:ne dof_idx = [2*e-1, 2*e, 2*e+1, 2*e+2]; Mg(dof_idx, dof_idx) = Mg(dof_idx, dof_idx) + Me; Kg(dof_idx, dof_idx) = Kg(dof_idx, dof_idx) + Ke; end %% 4. 固定端边界条件：删除前两个自由度 fixed_dofs = [1, 2]; free_dofs = 3:ndof; Mff = Mg(free_dofs, free_dofs); Kff = Kg(free_dofs, free_dofs); %% 5. 模态分析：用于验证模型并计算比例阻尼 [V, D] = eig(Kff, Mff); omega = sqrt(diag(D)); [omega, idx] = sort(omega); % 固有频率升序排列 freq_hz = omega / (2*pi); fprintf('有限元前三阶固有频率: %.2f Hz, %.2f Hz, %.2f Hz\n', ... freq_hz(1), freq_hz(2), freq_hz(3)); beta1 = 1.875104; f_analytic = beta1^2 * sqrt(E*I/(rho*A)) / (2*pi*L^2); fprintf('解析第一阶固有频率: %.2f Hz\n', f_analytic); %% 6. 比例阻尼，目标阻尼比 1% xi1 = 0.01; xi2 = 0.01; damp_mat = [1/(2*omega(1)), omega(1)/2; 1/(2*omega(2)), omega(2)/2]; ab = damp_mat \ [xi1; xi2]; alpha = ab(1); beta_coef = ab(2); Cff = alpha * Mff + beta_coef * Kff; %% 7. 状态空间模型 n = size(Mff, 1); Bf = zeros(n, 1); Bf(n-1) = 1; % 控制力作用在自由端平动自由度 A_s = [zeros(n), eye(n); -Mff\Kff, -Mff\Cff]; B_s = [zeros(n, 1); Mff\Bf]; C_s = eye(2*n); D_s = zeros(2*n, 1); %% 8. LQR 控制器设计 Q = blkdiag(1e4*eye(n), 50*eye(n)); % 位置权重、速度权重 R = 1e-4; % 控制力权重 [K_lqr, S, e_cl] = lqr(A_s, B_s, Q, R); Acl = A_s - B_s * K_lqr; %% 9. 闭环仿真：自由端初始位移 5mm t = 0:1e-4:2; x0 = zeros(2*n, 1); x0(n-1) = 0.005; sys_open = ss(A_s, B_s, C_s, D_s); sys_cl = ss(Acl, B_s, C_s, D_s); y_open = initial(sys_open, x0, t); y_cl = initial(sys_cl, x0, t); u_ctrl = -K_lqr * y_cl'; %% 10. 绘图 figure('Color', 'w'); subplot(2,1,1); plot(t, y_open(:, n-1)*1000, 'k--', 'LineWidth', 1.0); hold on; plot(t, y_cl(:, n-1)*1000, 'b-', 'LineWidth', 1.4); grid on; xlabel('时间 (s)'); ylabel('自由端位移 (mm)'); legend('开环响应', '闭环响应', 'Location', 'northeast'); title('悬臂梁自由端位移响应'); subplot(2,1,2); plot(t, u_ctrl, 'r-', 'LineWidth', 1.2); grid on; xlabel('时间 (s)'); ylabel('控制力 (N)'); title('LQR 控制输入');

简单过一遍关键段落。第 5 节用 eig(Kff, Mff) 解广义特征值问题，得到系统固有频率，一方面验证模型，另一方面为之后的瑞利阻尼提供频率依据。第 6 节把阻尼比设成 1%，通过前两阶频率反凑 α、β，这个操作比凭感觉拍一个 C 矩阵严谨得多。

第 7 节是整段代码的灵魂。A_s 的左上角是零矩阵，右上角是单位阵，左下角是 -Mff\Kff，右下角是 -Mff\Cff。这个结构对应的是 q 和 q̇ 的一阶状态方程。B_s 里 Mff\Bf 这一步，是把节点的集中力换算成加速度影响。注意 Bf(n-1) = 1，因为自由端节点平动自由度在删掉固定端后的局部索引是 n-1。

第 8 节的 lqr 函数会返回三个输出：增益 K、Riccati 方程的解 S 和闭环极点 e_cl。实际调试时可以看 e_cl 是否都在左半平面，以及虚部大小是否和系统模态对应。

3.3 Q、R 权矩阵怎么选才合理

Q、R 的选择没有唯一答案，但有很实用的经验规律。常见做法是让 Q 为对角矩阵，对每一位移状态取一个权值 q_d、每一位移速度状态取另一权值 q_v。这里的位置权重 1e4、速度权重 50，意味着控制目标里“位移收敛”比“速度归零”更优先，所以控制器会主动消耗能量去压振幅。

R 代表控制力的代价。R 越大，控制器越“抠门”，控制力小、振动衰减慢；R 越小，控制器越“敢下狠手”，衰减快但控制力峰值大。一个更系统的办法是 Bryson 规则：把期望的位移范围平方的倒数放进 Q 对应位置，把期望最大控制力的平方倒数放进 R。比如期望位移限幅 1mm，那么 q ≈ 1/(0.001)² = 10⁶；期望控制力不超过 1N，R ≈ 1/1² = 1。这个规则能给你一个合理的起步值，再根据仿真微调。

注意：LQR 的“最优”是在给定 Q、R 意义下的最优，不是绝对最优。所以调整权矩阵是控制设计里完全正当且必要的工作，不要觉得改参数就是“不严谨”。

4. 仿真结果分析与参数调优

4.1 开环对比闭环响应

代码里同时跑了开环和闭环两种响应。开环情况下，悬臂梁只靠 1% 的结构阻尼慢慢衰减，自由端位移从 5mm 缓缓回到零，这个过程可能要持续十几秒。闭环情况下，LQR 控制力会持续作用于自由端，相当于把能量“主动抽走”，位移曲线会明显短时间衰减到接近零。

这里特别值得看两个指标。一是调节时间，即响应进入稳态误差带（比如 ±5% 振幅）所需时间；二是控制峰值，即 u 的最大值。两个指标往往相互矛盾：想让调节时间更短，就得加大控制力峰值。你调 R 的时候，本质上就是在两者之间找一个自己能接受的平衡点。

4.2 控制力过冲与饱和问题

如果仿真里控制力峰值达到很大的值，比如几十牛甚至上百牛，而现实中的压电作动器最大只能输出几牛，那这个 LQR 设计就无法落地。解决办法有两个方向。第一，增大 R，直接压低控制力；第二，在 Q 里降低对位移的高权重，别再那么“着急”地收敛。

另一种情况是控制力曲线出现明显的高频抖振，这说明控制器在高频段增益过大，激发了模型中那些你以为不存在的高阶模态。处理手段通常是降低模型阶次，只保留前几阶模态，或者对控制输入做带宽限制。这部分问题我放到下一节的“控制溢出”里细讲。

4.3 从连续仿真到离散控制

本文的仿真把控制律当连续信号处理，但真实控制器一定是离散的。Matlab 里做离散化通常用 c2d 函数，采样时间要根据闭环系统带宽来定。一个经验法则是采样频率要高于系统闭环带宽的 10 倍。如果悬臂梁第一阶频率是 5Hz，但 LQR 把闭环极点拉高到 30Hz 对应频率，那采样频率至少 300Hz，也就是说采样周期要小于 3.3ms。

你可以在代码里加这几行：

Ts = 0.001; sys_d = c2d(ss(Acl, B_s, C_s, D_s), Ts);

再用离散模型做仿真，对比连续结果的差距。如果闭环特性发生明显变化，多半是采样率不够，需要进一步缩小 Ts。

5. 实际工程中常踩的坑与排查速查

5.1 五个典型坑及根因

第一个坑是控制溢出。模型阶次太高、R 太小，控制器在高频段增益过大，激发了有限元模型里物理意义不大但数值上存在的高阶模态。表现为仿真响应在高频上抖得厉害。解决办法是模型降阶，只保留前几阶主模态，或者增大 R、对控制输出做低通处理。

第二个坑是传感器位置。LQR 默认全状态可测，但实验里往往只有少数传感器。如果传感器恰好安装在某一阶模态的节点位置，这一阶模态在测量信号里就“隐身”了，状态估计会丢模态。做实验布置传感器之前，先用模态振型图查一下目标模态的节点位置。

第三个坑是自由度编号不一致。很多人在组装矩阵时用“先位移后转角”的排列，到施加边界条件或提取自由端响应时又搞混编号，结果仿真曲线完全不对。建议一开始就在注释里写清楚地步，并且用一个小段代码打印几个关键位置的自由度索引。

第四个坑是矩阵量纲问题。M 的数量级可能到 10⁻²，K 可能到 10⁴，两者相差六个数量级。状态空间矩阵里 -M⁻¹K 的数值可能很大，后续矩阵求逆和 Riccati 方程求解都会积累数值误差。可以把单位换成 mm、N、ton 重新统一量纲，或者对系统做平衡化处理。

第五个坑是激励方式不对。很多作业只做初始条件响应，但工程上更需要看正弦扫频或随机扰动下的抑制效果。LQR 设计本身处理的是状态反馈，输入 B_d 的建模方式和控制力 B_f 是分开的，扰动要从另一个通道进模型。可以用 lsim 来做这类仿真，而不是只用 initial。

5.2 常见问题速查表

5.3 扩展方向：LQG、降阶与实验台验证

如果实际系统只有位移传感器、没有速度传感器，可以直接在 LQR 基础上扩展成 LQG。做法是先做 Kalman 滤波来估计状态，再配合 LQR 增益输出控制力。Matlab 里有 lqg 命令可以一步完成。

另一个常见扩展是做模态截断。先对 M、K 做模态分解，只保留前 r 阶模态，再把控制输入和测量投影到模态坐标上，得到一个 r 维降阶模型。这样 LQR 控制器阶数更小、更稳，也更容易在嵌入式控制器上实现。

我建议你在跑通全状态反馈之后，试着把执行器位置从自由端改到梁中间、或者改到靠近固定端的位置，重新运行仿真。你会发现同一个 LQR 参数在不同位置上的控制效果差很多，这就是“作动器布局影响能控性”的直观体验，比看书上的能控性格拉姆矩阵有意思得多。

6. 写在最后：实验验证里的几条心得

这套代码我在不同项目上改过很多次，最深刻的体会是：有限元模型的自由度编号和边界条件是所有问题的根源，比 LQR 的公式难搞多了。每次仿真结果不对，我都先回头检查 M、K 的组装和固定端删行，而不是急着调 Q、R。

另一个经验是，仿真里“最优”的控制器到了实验台上往往需要再降增益。因为有限元模型永远比真实结构“干净”，真实结构的模型误差、传感器噪声和执行器饱和都会让高增益控制器变得不安定。所以我在实验前一定会把仿真里的 R 上调两个数量级再试试，然后再慢慢往下压。

如果你想把模型变得更真实，可以考虑加上压电片本身的机电耦合刚度、传感器噪声模型，以及执行器的饱和特性。这些扩展会让控制器不再“纸上谈兵”，也是从课程作业走向工程应用最重要的一步。先把这里的悬臂梁 LQR 跑通、跑透，后面学 LQG、H∞、自适应控制都会轻松很多。