企业官网建设流程全解析

简介：针对二自由度机械臂的逆动力学分析需求，一套MATLAB源码包给出了欧拉-拉格朗日动力学模型的完整实现，适合机器人控制、机械工程及自动化相关专业学生与研究人员学习和二次开发。代码以关节位置、速度、加速度为输入，利用拉格朗日欧拉公式求解关节力矩或力，能够反映关节间的耦合与非线性特性，可依据期望运动反解驱动力矩，也为后续控制律设计提供了可验证的模型基础。包内共7个文件，6个.m脚本分别实现动能、势能、惯性矩阵、离心力与科氏力及整体动力学计算，另含1张PNG机械臂示意图，便于对照结构理解；压缩包整体仅53KB，轻量易读，适合快速移植与课程实验。目前已有305人学习下载，对需要掌握机器人建模与MATLAB实现的读者，尤其适合作为课程设计与毕业设计的参考资料。

1. 从Simulink调PID到手抖：为什么二自由度机器人要先写好欧拉-拉格朗日模型

在MATLAB里调二自由度机械臂的控制，很多人第一步就栽在动力学上：给关节1一个阶跃力矩，关节2会先反着跑一下；或者两个连杆同时运动时，Simulink里的PID怎么调都压不住振荡。追到底，你会发现问题往往不在控制器，而是手里那个用于仿真的二自由度机器人欧拉-拉格朗日动力学模型本身就是错的。欧拉-拉格朗日建模法用系统的动能、势能推导出关节力矩与角度、角速度、角加速度之间的完整约束关系，把惯性耦合、科氏力、离心力和重力四类效应显式写出来，是轨迹跟踪、计算力矩控制、阻抗控制以及Simscape虚拟仿真验证的公共底座。这篇文章配套MATLAB源码，带你从动能势能推导、符号计算、ode45数值求解一路做到Simulink闭环，并给出三个能立刻判断模型对错的验证手法。适合刚接触机械臂仿真的研究生，也适合需要快速把自有动力学模型跑起来的工程师。

2. 建立二自由度机械臂的欧拉-拉格朗日方程：从动能势能到M、C、G矩阵

推导平面二自由度机械臂的动力学历来是固定套路：定义坐标、写动能势能、代入拉格朗日方程、整理成标准形式。这里不急着背公式，先把物理图像建立起来，后面写MATLAB源码时才知道每一步在算什么。

2.1 二自由度旋转机械臂的坐标约定与参数表

考虑平面2R机械臂，两个旋转关节的转轴都垂直于纸面。关节1连接杆1与基座，关节2连接杆1与杆2。角度定义上，q1表示杆1与水平线的夹角，q2表示杆2相对杆1的关节角；质心位置用lc1和lc2表示，分别是从各自关节中心到连杆质心的距离。小角度假设在这里不成立，后续推导全部保留cos(q2)、sin(q2)等非线性项。

这些值本身不是关键，重要的是量级关系：关节2负载明显小于关节1，惯量比大概是2:1。实际项目中如果从URDF或CAD拿到的是质心坐标和惯性张量，先折算成这里的lc和I再往下走，否则后面Simulink里调PID的规律会失真。

2.2 动能与势能：交叉项就是耦合的来源

写出两个连杆质心的速度平方。杆1质心速度只有关节1的转动贡献：v₁² = lc1²·q̇₁²。杆2质心速度要复杂一些，是"关节2随杆1转动的牵连速度"与"杆2相对杆1转动的相对速度"的向量合成：

v₂² = l1²·q̇₁² + lc2²·(q̇₁+q̇₂)² + 2·l1·lc2·q̇₁·(q̇₁+q̇₂)·cos(q2)

第三项是交叉项，它出现的原因是牵连速度与相对速度之间的夹角刚好由q2决定，于是内积里带出cos(q2)。这个交叉项是整篇模型的关键：没有它，方程退化成两个独立的单摆；有它，才会出现后面M矩阵里的cos(q2)和C矩阵里的sin(q2)，也就是两个关节之间的惯性耦合和科氏力。

动能由平动加转动组成：

K = 0.5·m1·v₁² + 0.5·I1·q̇₁² + 0.5·m2·v₂² + 0.5·I2·(q̇₁+q̇₂)²

势能以关节1轴心为零点，按高度计算：

V = m1·g·lc1·sin(q1) + m2·g·(l1·sin(q1) + lc2·sin(q1+q2))

写出拉格朗日函数L = K - V，对每个关节代入欧拉-拉格朗日方程：

d/dt(∂L/∂q̇ᵢ) - ∂L/∂qᵢ = τᵢ

这里的τᵢ是作用在关节i上的广义力矩。展开后方程是一堆二阶非线性微分方程，直接看很乱，工程上习惯把它整理成下面这种标准形式。

2.3 标准形式 M(q)·q̈ + C(q,q̇)·q̇ + G(q) = τ

把上面的展开式重新聚合，得到机械臂动力学通用结构。惯性矩阵M(q)是对称的：

M₁₁ = m1·lc1² + I1 + m2·(l1² + lc2² + 2·l1·lc2·cos(q2)) + I2
M₁₂ = M₂₁ = m2·(lc2² + l1·lc2·cos(q2)) + I2
M₂₂ = m2·lc2² + I2

重力项G(q) = ∂V/∂q：

G₁ = m1·g·lc1·cos(q1) + m2·g·(l1·cos(q1) + lc2·cos(q1+q2))
G₂ = m2·g·lc2·cos(q1+q2)

科氏与离心项可以写成C(q,q̇)·q̇。令B = m2·l1·lc2，则：

C = [ -B·sin(q2)·q̇₂, -B·sin(q2)·(q̇₁+q̇₂);
B·sin(q2)·q̇₁, 0 ]

C矩阵第一行描述的是杆2运动对关节1产生的反作用力矩，第二行描述杆1运动对关节2产生的科氏力矩。注意C不是对称的，这符合科氏力的物理特性：两个关节互相施加的"扰动"方向相反、强度不同，后面验证模型时会用到这个性质。

关于符号约定，这里特意把G(q)写成加号形式，对应方程M·q̈ + C·q̇ + G = τ。有些教材写成M·q̈ + C·q̇ + G = τ + τ_ext，或者把重力项挪到左侧，只要最终仿真代码里按同一个约定写，结果一致。最容易出错的恰恰是符号不统一。

2.4 用MATLAB符号工具箱把推导做成源码

手推M和G容易漏项，尤其M₁₂那个交叉项，少了它整个仿真的耦合特性就丢了。用MATLAB Symbolic Toolbox可以一步到位，常见的做法是直接用hessian和jacobian，不需要手工对L求导。

% 二自由度机器人欧拉-拉格朗日模型：符号推导 M、G 矩阵 syms q1 q2 q1d q2d real syms m1 m2 l1 lc1 lc2 I1 I2 g real % 质心速度平方：q1 为杆1与水平夹角，q2 为杆2相对转角 v1_sq = lc1^2 * q1d^2; v2_sq = l1^2*q1d^2 + lc2^2*(q1d+q2d)^2 ... + 2*l1*lc2*q1d*(q1d+q2d)*cos(q2); % 动能 = 平动动能 + 绕质心转动动能 K = 0.5*m1*v1_sq + 0.5*I1*q1d^2 ... + 0.5*m2*v2_sq + 0.5*I2*(q1d+q2d)^2; % 势能：取关节1轴心为势能零点 V = m1*g*lc1*sin(q1) + m2*g*(l1*sin(q1) + lc2*sin(q1+q2)); % 动能是广义速度的二次型: K = 0.5*qd'*M*qd，hessian 直接给出 M M = simplify(hessian(K, [q1d, q2d])); G = simplify(jacobian(V, [q1, q2])'); % 列向量，与 M*qdd + C*qd + G 一致 fprintf('M11 = %s\n', M(1,1)); fprintf('M12 = %s\n', M(1,2)); fprintf('M22 = %s\n', M(2,2)); fprintf('G1 = %s\n', G(1)); fprintf('G2 = %s\n', G(2));

hessian(K, [q1d, q2d])返回K对广义速度的二阶偏导矩阵，因为动能对q̇是二次齐次的，这个二阶导矩阵就等于2倍的系数矩阵，也就是M(q)。jacobian(V, [q1, q2])'得到重力列向量，这里转置是为了保持列方向与方程左侧一致。运行后打印的结果应该和2.3节手推的公式逐项相同。

C矩阵没有用符号推导，原因是Christoffel符号展开比较繁琐，而2.3给出的解析式足够简单。如果你坚持要符号化推导C，可以借助dM/dt和偏导组合来提取，但实际项目中很少这么做。符号推导的价值在于验证M和G，C矩阵的正确性用第5章的斜对称校验来保证。

3. 把欧拉-拉格朗日模型变成可求解的MATLAB源码：m函数与ode45数值求解

符号表达式不能直接喂给ode45，必须转换成数值函数。常见做法有两种：用matlabFunction把符号表达式转成函数句柄，或者手写一个独立的m函数。手写m函数的优势在于可读性好、便于加断点、也方便后续封装成Simulink模块，下面采用这种方案。

3.1 从符号表达到可复用m函数

动力学模型输入是当前时刻的状态x = [q1; q2; q̇₁; q̇₂]和外力矩tau = [τ₁; τ₂]，输出是状态导数xdot = [q̇₁; q̇₂; q̈₁; q̈₂]。把M、C、G全部显式算出来之后，加速度由q̈ = M⁻¹·(tau - C·q̇ - G)解出。以下是完整实现：

function xdot = twoR_dynamics(t, x, tau, p) % 二自由度机器人欧拉-拉格朗日动力学 % 状态 x = [q1; q2; q1d; q2d]，输入 tau = [tau1; tau2] % 返回 xdot = [q1d; q2d; q1dd; q2dd] q1 = x(1); q2 = x(2); q1d = x(3); q2d = x(4); % 惯性矩阵 M(q) M11 = p.m1*p.lc1^2 + p.I1 ... + p.m2*(p.l1^2 + p.lc2^2 + 2*p.l1*p.lc2*cos(q2)) + p.I2; M12 = p.m2*(p.lc2^2 + p.l1*p.lc2*cos(q2)) + p.I2; M22 = p.m2*p.lc2^2 + p.I2; M = [M11 M12; M12 M22]; % 科氏/离心项，B = m2*l1*lc2 B = p.m2*p.l1*p.lc2; C = [-B*sin(q2)*q2d, -B*sin(q2)*(q1d+q2d); B*sin(q2)*q1d, 0]; Cqd = C*[q1d; q2d]; % 重力项 G(q) G = [p.m1*p.g*p.lc1*cos(q1) ... + p.m2*p.g*(p.l1*cos(q1) + p.lc2*cos(q1+q2)); p.m2*p.g*p.lc2*cos(q1+q2)]; % 解出角加速度，M qdd = tau - Cqd - G qdd = M \ (tau - Cqd - G); xdot = [q1d; q2d; qdd(1); qdd(2)]; end

函数签名里的t参数在本次实现中没有用到，但必须保留。ode45要求被积函数的签名是f(t, x)，即使方程不显含时间也要占住这个位置。参数用结构体p传递而不是写成全局变量，是为了方便之后做参数扫描：改p的值不需要动函数内部代码。qdd用M \ 而不是inv(M) *，前者走数值线性求解，对2×2矩阵差异不大，但这个习惯在更高自由度模型里能明显提升精度和速度。

3.2 用ode45计算自由运动和阶跃力矩响应

模型函数就绪后，剩下的就是调用MATLAB自带求解器。先跑两个典型工况：零力矩自由摆动和常值力矩持续加速。

% 参数与初始条件 p.l1 = 1.0; p.lc1 = 0.5; p.lc2 = 0.5; p.m1 = 2.0; p.m2 = 1.5; p.I1 = 0.15; p.I2 = 0.08; p.g = 9.81; x0 = [pi/4; pi/6; 0; 0]; % 初始角度 q1=45°, q2=30° tspan = [0 5]; % 工况1：零力矩，机械臂在重力作用下自由摆动 u_zero = @(t) [0; 0]; opts = odeset('RelTol', 1e-6, 'AbsTol', 1e-8); [t, x] = ode45(@(t, x) twoR_dynamics(t, x, u_zero(t), p), ... tspan, x0, opts); subplot(2,1,1); plot(t, x(:,1)*180/pi, t, x(:,2)*180/pi, 'LineWidth', 1.2); legend('q1 (deg)', 'q2 (deg)'); grid on; ylabel('角度'); xlabel('t (s)'); % 工况2：关节1、2分别施加常值力矩 5N·m 和 3N·m u_const = @(t) [5; 3]; [t2, x2] = ode45(@(t, x) twoR_dynamics(t, x, u_const(t), p), ... tspan, x0, opts); figure; plot(t2, x2(:,3), t2, x2(:,4), 'LineWidth', 1.2); grid on; legend('q1d', 'q2d'); xlabel('t (s)'); ylabel('角速度(rad/s)');

自由运动工况下，两个关节会做类似摆的运动，q2由于没有力矩输入会在重力作用下往复振荡，振荡过程中还能看到q1和q2之间的相位差，这正是耦合项的体现。常值力矩工况下系统持续加速、角速度单调增加是正常现象，因为模型里没有任何人为阻尼。如果看到角速度发散或者曲线出现突变，优先回头查M矩阵是否对称、G项符号是否反了。

3.3 ode45求解器参数怎么设：误差容差与步长控制

求解这类刚柔混合的机械系统，最常调整的是误差容差。默认的RelTol=1e-3对控制初筛够用，但轨迹要再喂给下游控制器或做参数辨识时，建议收紧到1e-6。

注意MaxStep默认并非固定，MATLAB会按积分区间自动决定一个上限。如果你在结果里看到明显的高频锯齿，先不要怀疑算法，把MaxStep显式设成1e-3再跑一次，往往能定位是数值问题还是模型问题。AbsTol和RelTol的配合原则是：如果角度量级在0.1 rad以下，AbsTol至少要比真实幅值小一个数量级，否则小角度区间被相对误差放大器放大，轨迹会显得很毛糙。

4. 把二自由度欧拉-拉格朗日模型装进Simulink做PID轨迹跟踪仿真

纯脚本仿真适合验证模型本身，但做控制器设计就费劲了：反馈回路、力矩饱和、示波器观察、参数扫描每样都得手写。把动力学模型封装进Simulink是更常见的工程做法，后面还能对接代码生成。

4.1 为什么建议用Simulink而不是纯脚本

Simulink的价值在于把"被控对象"和"控制器"分成两个清晰的可替换模块。调试PD参数时，只需要改两个增益模块的数值，不需要重新运行整个脚本；观察力矩曲线、跟踪误差、关节加速度也远比命令行快捷。另外，如果项目后续要部署到实时硬件，Simulink模型里的控制器可以直接生成C代码，而手写脚本则要多做一层翻译。

封装方式选择上，不建议一上来写Level-2 MATLAB S-Function，样板代码多、排错麻烦。常用做法是用MATLAB Function块加Integrator积分环，内部直接调用第3章的twoR_dynamics函数。这样既复用了已验证的源码，又保留了Simulink的可视化优势。

4.2 用MATLAB Function块封装动力学源码并接线

在Simulink模型里新建一个MATLAB Function模块，输入为力矩tau（2×1）和状态x（4×1），输出为状态导数xdot（4×1）。积分器设置初始条件为x0 = [pi/4; pi/6; 0; 0]，积分器输出x反馈回MATLAB Function输入，形成一个完整的连续积分闭环。

MATLAB Function块内部代码，直接用参数结构体初始化，代码清晰且不依赖工作区变量：

function xdot = fcn(tau, x) % Simulink中复用二自由度欧拉-拉格朗日动力学源码 p = struct('l1',1.0,'lc1',0.5,'lc2',0.5, ... 'm1',2.0,'m2',1.5, ... 'I1',0.15,'I2',0.08, 'g',9.81); xdot = twoR_dynamics(0, x, tau, p); end

这段代码里t直接传0，因为动力学本身不显含时间，传任意值都不影响结果。如果你在多个模型里复用，更规范的做法是把p定义成模型工作区变量，然后用参数对象传入，避免每个MATLAB Function块都重复构造结构体。

4.3 PID控制器参数与重力补偿调试顺序

动力学模型配上最简单的PD控制器就能做轨迹跟踪，但直接给两个关节各加一组增益往往是调不通的。问题出在重力项：低速运动时，重力矩占主导，纯PD要靠静差来抵抗重力，表现为稳态误差大、响应迟缓。工程上第一步先把重力前馈加上：

τ = Kp·(q_des − q) + Kd·(q̇_des − q̇) + G(q)

其中G(q)就是第2章推导出的重力项。

% 参考控制器：重力补偿 + PD，Kp Kd 为2x2对角阵 e = q_des - [q1; q2]; ed = qd_des - [q1d; q2d]; tau = Kp*e + Kd*ed + G; % G 来自 twoR_dynamics 中的表达式

调试顺序我一般这样安排：先把Kp、Kd全部置0，只保留重力前馈，看机械臂能否静止在目标位置；能静止，说明G项符号和数值正确。然后逐个关节加Kp，先关节1后关节2，每加一次观察超调量；出现等幅振荡时，再按Kp的0.2到0.5倍加Kd。

Kp2明显小于Kp1不是因为关节2不重要，而是因为关节2的等效惯量随q2变化，在q2接近0度时等效惯量最大，Kp给大了容易在伸展位形激励出高频振荡。

4.4 Simulink代数环与求解器设置

注意：MATLAB Function块里如果直接从输出端取状态再参与同一时刻的运算，Simulink会报代数环错误。解决方法是把Integrator的输出作为唯一状态来源，控制器与动力学模块之间只走前向信号，不在同一个Function块内部做瞬时闭环。

代数环最常见的触发场景是：想把加速度q̈直接反馈给控制器做前馈，但q̈又是当前力矩的函数，于是Simulink陷入"先有鸡还是先有蛋"的迭代。正确处理办法是让控制器只使用Integrator输出的q和q̇，前馈量用期望轨迹的二阶导数计算，而不是用实际加速度。仿真采样时间也值得确认：如果只是连续仿真，变步长ode45即可；若模型里有零阶保持器或离散控制器，把求解器改成ode23t更稳妥，固定步长仿真则用ode4，步长设为1ms，与真实控制周期的量级匹配。

5. 验证二自由度欧拉-拉格朗日模型正确性的三个实用技巧

模型写完不等于模型正确，下面三个验证手法能快速暴露多数推导和实现错误，建议每一步都做一遍再开始调控制器。

5.1 能量守恒检查：零输入下总能量应当恒定

无外力矩时，系统总能量E = K + V应该保持不变。数值积分会有缓慢漂移，但量级应当远小于单次摆动周期内的能量变化。把能量曲线画出来，如果有明显上升或下降趋势，说明M矩阵或C矩阵有错。

E = zeros(size(t)); for k = 1:numel(t) q2k = x(k,2); M11 = p.m1*p.lc1^2 + p.I1 + p.m2*(p.l1^2 + p.lc2^2 + 2*p.l1*p.lc2*cos(q2k)) + p.I2; M12 = p.m2*(p.lc2^2 + p.l1*p.lc2*cos(q2k)) + p.I2; M22 = p.m2*p.lc2^2 + p.I2; qd = x(k,3:4)'; K = 0.5*qd'*[M11 M12; M12 M22]*qd; V = p.m1*p.g*p.lc1*sin(x(k,1)) + p.m2*p.g*(p.l1*sin(x(k,1)) + p.lc2*sin(x(k,1)+x(k,2))); E(k) = K + V; end plot(t, E);

能量曲线在5秒内漂移不超过初始值的1%就说明积分精度达标。如果初始段能量就跳变，多半是M矩阵漏了交叉项，比如M₁₂少了m2·l1·lc2·cos(q2)这一项。

5.2 检查M矩阵对称正定与Ṁ − 2C斜对称

M矩阵对称正定是机械臂模型的物理底线，任何一处符号写错都会破坏对称性。更强的是Ṁ − 2C应满足斜对称关系，即S + Sᵀ = 0。这个性质在Lyapunov稳定性证明中是核心，也最适合用来检验C矩阵的实现。

% 抽取某个时刻做校验 q1 = 0.3; q2 = 0.5; q1d = 1.2; q2d = -0.8; [M, C] = get_MC(q1, q2, q1d, q2d, p); % 从twoR_dynamics拆分出M和C Mdot = [-2*B*sin(q2)*q2d, -B*sin(q2)*q2d; -B*sin(q2)*q2d, 0]; S = Mdot - 2*C; assert(norm(S + S') < 1e-12);

其中B = m2·l1·lc2。对于教科书上给出的C矩阵，斜对称条件严格成立；如果你用的是Christoffel符号推导出来的C，这个条件也应当满足。误差大于1e-10说明C矩阵实现和M矩阵的时间导数不是同一套动力学，至少有一处错了。

5.3 退化单摆对照：锁死关节2比较周期

让q2 ≡ 0、q̇₂ ≡ 0，模型就退化成一根由两个连杆拼成的单摆。此时关于关节1的等效转动惯量J = m1·lc1² + I1 + m2·(l1+lc2)² + I2，等效重心矩L = m1·g·lc1 + m2·g·(l1+lc2)，小角度振荡频率应当满足ω² = L/J。在模型里把q2初始值设成0，给q1一个1°的小初始偏角，测振荡周期，与2π·sqrt(J/L)对照，误差应在1%以内。这一步能同时暴露重力项G的符号问题和M₁₂的交叉项错误，把模型质量一次性锁死。三个验证都通过后，这个二自由度机器人欧拉-拉格朗日模型就可以放心交给下游的轨迹规划、力控制或参数辨识使用了。

本文还有配套的精品资源，点击获取