MATLAB 柔性机械臂抑振
- 刚柔耦合建模（Euler-Bernoulli + 集中质量）
- LQR 状态反馈抑振仿真
- 实验数据回放/对比接口（含双目视觉采集格式）
一、列表
main_flexibleLQR.m% 主脚本（一键跑通）flexibleDynamics.m% 刚柔耦合动力学（含模态阻尼）lqrDesign.m% LQR 权重自动调节（可实验迭代）experimentInterface.m% 实验数据回放/对比（双目视觉 CSV 格式）plotFlexibleResult.m% 仿真-实验对比图（含残差振动能量）
二、刚柔耦合模型（单杆弹性简化）
- 杆长 L = 0.4 m，线密度 ρA = 0.12 kg/m，EI = 60 N·m²
- 集中质量 m_tip = 0.05 kg（代表末端小球）
- 前 2 阶模态即可覆盖 95% 振动能量
- 状态向量：x = [θ, θ̇, q₁, q̇₁, q₂, q̇₂]ᵀ
（θ：关节角，qᵢ：模态坐标）
动力学方程（已线性化）
M ẍ + C ẋ + K x = B u
其中 u 为关节力矩（N·m）
三、核心代码（精简版）
- 主脚本
%% 0 参数
clear; clc; close all
global M C K B % 供动力学使用
L = 0.4; rhoA = 0.12; EI = 60; m_tip = 0.05;
nMod = 2; % 模态阶数
dt = 1e-3; T = 5; % 5 s 仿真
t = 0:dt:T;%% 1 建立状态空间
[M,C,K,B] = flexibleModel(L,rhoA,EI,m_tip,nMod);
A = [zeros(3+nMod) eye(3+nMod);-M\K -M\C];
B = [zeros(3+nMod,1); M\B];
C = eye(3+nMod); D = 0;
sys = ss(A,B,C,D);%% 2 LQR 设计
Q = diag([1e4 1e2 1e1 1e1 1e1 1e1]); % 权重可实验迭代
R = 0.1;
K_lqr = lqr(A,B,Q,R);%% 3 期望轨迹（软启动圆轨迹，抑振关键）
theta_des = @(t) 0.5*sin(2*pi*0.5*t); % 0.5 Hz 圆
dtheta_des = @(t) 0.5*2*pi*0.5*cos(2*pi*0.5*t);
q_des = zeros(nMod,1); % 模态目标=0
x_des = [theta_des(0); q_des; dtheta_des(0); zeros(nMod,1)];%% 4 仿真（闭环 LQR）
x0 = [0; zeros(3+nMod-1,1)]; % 初始静止
[t_sim, x_sim] = lsim(sys, ...@(t,x) -K_lqr*(x - x_des_func(t)), t, x0);%% 5 实验数据对比（可选）
expData = experimentInterface('flexible_exp.csv'); % 双目视觉 CSV
plotFlexibleResult(t_sim, x_sim, expData);
- 动力学模型函数
function [M,C,K,B] = flexibleModel(L,rhoA,EI,m_tip,nMod)
% 模态频率与振型（ clamped-free ）
for i = 1:nModbeta(i) = (2*i-1)*pi/2; % 近似omega(i) = beta(i)^2 * sqrt(EI/(rhoA*L^4));
end
M = diag([1, rhoA*L/2, rhoA*L/2] + [0, m_tip, m_tip]); % 广义质量
C = diag([0.02, 0.02, 0.02]); % 模态阻尼 2%
K = diag([0, omega.^2]); % 刚度
B = [1; zeros(nMod,1)]; % 输入矩阵
end
- 轨迹函数（软启动，抑振关键）
function x_ref = x_des_func(t)
theta = 0.5*sin(2*pi*0.5*t);
dtheta = 0.5*2*pi*0.5*cos(2*pi*0.5*t);
x_ref = [theta; 0; 0; dtheta; 0; 0]; % 模态目标=0
end
参考代码 柔性机械臂的抑振实验研究 www.3dddown.com/cna/64105.html
四、实验接口（双目视觉 CSV 格式）
CSV 列：time, theta_exp, q1_exp, q2_exp, theta_dot_exp, q1_dot_exp, q2_dot_exp
函数自动对齐时间基线，计算残差振动能量：
function exp = experimentInterface(csvFile)
exp = readtable(csvFile);
exp.time = exp.time - exp.time(1); % 对齐 t=0
exp.Evib = exp.q1_exp.^2 + exp.q2_exp.^2; % 模态振动能量
end