☰
基于MATLAB的Stewart平台逆解与动力学仿真实践
2026/10/15 8:04:39 网站建设 项目流程

简介:围绕MATLAB环境下的Stewart并联平台运动学逆解与动力学仿真,这份PDF资料面向机器人控制及并联机构研究人员,系统解决从坐标变换、支腿长度计算到雅可比矩阵应用、动力学建模与三维仿真的完整问题。文中基于牛顿-欧拉法推导动力学方程,并结合Simscape Multibody完成三维仿真,针对奇异位形检测、轨迹规划不合理等实际问题给出排错思路;内容预览中的多段MATLAB代码涵盖逆解函数、实时可视化、螺旋轨迹生成与S-Function动力学模型,并提示了eul2rotm旋转顺序、drawnow limitrate性能优化、ode15s求解器选择等实战细节。资料为单份PDF文档,压缩包大小116KB,便于离线阅读。目前已有189人学习,可帮助科研人员与技术爱好者掌握从理论建模、算法实现到仿真验证的完整链路,避免常见陷阱。

1. 为什么说 Stewart 平台的逆解是整个仿真链条里最难绕过去的坎

把“基于MATLAB的Stewart并联平台运动学逆解与动力学仿真”做成能落地的东西,我首先要劝你收起“逆向不就是算个距离吗”的心态。Stewart 平台有六个自由度,动平台位姿由位置加姿态决定,逆解要从给定的六个自由度到位计算六根驱动杆长度;公式确实只有两步——旋转坐标、求范数——但真正决定精度的是坐标约定、铰链布局和旋转顺序。坐标错一处,所有杆长都会连带出错;逆解错,后面的动力学仿真就是在错误输入上做无意义计算。

这篇记录适合手里已有平台机械尺寸、想在仿真阶段就把控制逻辑和极限位姿验证清楚的人。这里不会谈太多抽象理论,而是先把模型立住,再给复现代码、参数设定和连续踩坑记录。读完你能判断这个方向值不值得投入,也知道杆长跳变、仿真发散时该改哪里。

2. 先把数学模型立住:坐标系、铰链点与逆解推导

2.1 铰链布局怎么定:半径、半角与初始高度

动笔写任何一行 MATLAB 代码之前,我习惯先做一张表,把六对铰链点的空间坐标固定下来。Stewart 平台的常见布局是:静平台和动平台各有一个铰链圆,两个圆上各有六个铰链点,按“两两一组、组间夹角 120°”的方式分布。两个铰链点绕组对称轴的夹角叫半角,通过半径和半角能唯一确定每个铰链点的位置。

某模拟项目X使用的参数如下:

参数符号取值说明
静平台半径Rb0.6 m六个静铰链所在圆半径
动平台半径Rt0.35 m六个动铰链所在圆半径
静铰链半角bet_b10°静平台每对铰链的半夹角
动铰链半角bet_t5°动平台每对铰链的半夹角
初始高度h00.75 m动坐标系原点在静系中的 z 坐标

这里有个容易被忽略的细节:动平台的六个铰链点定义在动坐标系里,当平台处于初始位姿(不旋转、原点在 h0)时,这六个点应当落在静坐标系里半径为 Rt、高度为 h0 的圆上。半角选得是否合理,直接决定后续杆长范围。若两平台半角相同且组间相位完全对齐,在某个姿态下会出现连杆互相平行、机构失去刚度的奇异位形;实际中常把动平台组间相位相对静平台偏转一个角度,降低奇异风险。

生成铰链点坐标的 MATLAB 函数如下:

% 生成静平台与动平台的铰链点坐标 % 输出: Ab(3x6) 静铰链, At0(3x6) 动铰链(动系坐标) function [Ab, At0] = hull_hinge(Rb, Rt, bet_b, bet_t) % 静平台六点: 分三组, 每组绕0, 120, 240度对称 ang_b = [bet_b, -bet_b, ... 2*pi/3+bet_b, 2*pi/3-bet_b, ... 4*pi/3+bet_b, 4*pi/3-bet_b]; % 动平台六点: 组间相位整体偏转60度, 降低奇异风险 ang_t = [bet_t+pi/3, -bet_t+pi/3, ... 2*pi/3+bet_t+pi/3, 2*pi/3-bet_t+pi/3, ... 4*pi/3+bet_t+pi/3, 4*pi/3-bet_t+pi/3]; Ab = [Rb*cos(ang_b); Rb*sin(ang_b); zeros(1,6)]; At0 = [Rt*cos(ang_t); Rt*sin(ang_t); zeros(1,6)]; end

代码中把动平台角度整体加了 60°,是为了在初始位姿就把六条腿错开。你可以在命令行调用一次:[Ab, At0] = hull_hinge(0.6, 0.35, 10*pi/180, 5*pi/180);,然后用scatter3(Ab(1,:), Ab(2,:), Ab(3,:))标出两组点。能看到六个静铰链点和六个动铰链点不重合、连线自然展开,说明布局对路子。

提示:如果不想让铰链点都在圆周上,也可以按矩形或圆加偏心的方式生成,但对纯运动学逆解来说,圆周布局计算量最小,也最容易对照手算结果。

2.2 姿态表示与旋转顺序:为什么 ZYX 约定要统一

逆解中最常翻车的不是公式,而是姿态表示。MATLAB Robotics System Toolbox 提供eul2rotm(angles, sequence),默认顺序是'ZYX',也就是先绕 Z 轴、再绕新 Y 轴、最后绕新 X 轴的内旋。机械臂控制里'ZYX'很常见,但很多从视觉端切入的开发者会习惯'XYZ'。两者算出来的旋转矩阵完全不同,而且不是简单差一个符号。

我一般会在项目最开始就写一个验证脚本,把同一组角度分别转成两种矩阵并做差。代码很简单:

% 验证不同旋转顺序造成的差异 angles = [20, 15, 10] * pi/180; % 横滚, 俯仰, 偏航(弧度) R_zyx = eul2rotm(angles, 'ZYX'); R_xyz = eul2rotm(angles, 'XYZ'); disp('ZYX - XYZ 矩阵差:'); disp(R_zyx - R_xyz);

这段代码不会报错,但会输出一个非零矩阵,说明同样一组欧拉角在不同 sequence 参数下描述的是完全不同的姿态。在有视觉或惯性传感器参与的系统中,传感器输出频率高、姿态更新快,一旦顺序不一致,逆解输出的六根杆长在轨迹上就会表现为突然跳变。

所以项目里应当建立一个约定:所有位姿都用位置向量加'ZYX'欧拉角表达,所有接触外部接口的地方都强制校验旋转矩阵为正交矩阵。校验方法是norm(R'*R - eye(3)),如果超过 1e-6 就需要停下检查,这比肉眼盯着杆长曲线判断更快。

2.3 逆解公式:位姿到六根杆长的两步推导

逆解本身只有两步。设动平台在静系中的位置为 p,姿态旋转矩阵为 R,静平台第 i 个铰链坐标 Ab_i,动平台第 i 个铰链在动系中的坐标为 At0_i,则该铰链在静系中的实际坐标是:

At_i = p + R * At0_i

第 i 根驱动杆的杆长就是两个点之间的欧几里得距离:

L_i = || At_i - Ab_i ||_2

这里不需要解方程,也不需要做矩阵求逆,逐个点计算即可。写成 MATLAB 函数非常短:

% Stewart 平台位置逆解 % 输入: pos(3x1) 动系原点位置; R(3x3) 姿态矩阵; % Ab(3x6) 静铰链; At0(3x6) 动铰链(动系) % 输出: L(6x1) 六根驱动杆长度 function L = stewart_inv(pos, R, Ab, At0) L = zeros(6,1); for i = 1:6 At_w = pos(:) + R * At0(:,i); % 动铰链转到静系 L(i) = norm(At_w - Ab(:,i)); % 杆长即两点距离 end end

参数说明:pos必须是一个 3 行 1 列的列向量,否则pos(:)虽然能修正形状,但会掩盖外部传参的问题,所以调用前统一pos = reshape(pos,3,1)。R必须是姿态旋转矩阵,不能是欧拉角;如果使用者直接传角度向量,R * At0就会因为维度不匹配直接报错。这个报错其实是好事,它迫使接口保持清晰。

逆解不需要考虑六根杆的安装角度,因为杆长由两个球铰固定点之间的距离决定。相比正解的迭代求解,逆解开销极小,一万个位姿点在一台普通机器上也就是几十毫秒量级,完全可以做实时控制前的预计算。

3. 在 MATLAB 里把逆解跑起来:轨迹代码与参数设定

3.1 生成一段连续位姿轨迹并计算杆长变化

逆解函数写好后,下一步不是直接上 Simulink,而是在 MATLAB 脚本里给动平台规划一段连续轨迹,观察六根杆长随时间的变化是否平滑。这是成本最低、也最能暴露问题的验证手段。

常见做法是让动平台沿正弦形的轨迹运动:高度做低频正弦波动,横滚和偏航也做小幅度摆动。这样既能覆盖位置和姿态两个自由度,又不会在起步阶段就撞上奇异位形。

%% 连续轨迹逆解 Ab = ...; % 由 hull_hinge 生成, 此处省略运行结果 At0 = ...; % 由 hull_hinge 生成, 此处省略运行结果 t = linspace(0, 10, 1000); L = zeros(6, length(t)); for k = 1:length(t) tk = t(k); % 位置: z 方向波动, 水平方向小幅偏移 pos = [0.02*sin(0.5*tk); 0.01*cos(0.3*tk); 0.75 + 0.05*sin(2*tk)]; % 姿态: 横滚, 俯仰, 偏航均是小角度(单位弧度) eul = [3*pi/180*sin(tk), ... 2*pi/180*cos(0.7*tk), ... 5*pi/180*sin(0.3*tk)]; R = eul2rotm(eul, 'ZYX'); L(:,k) = stewart_inv(pos, R, Ab, At0); end plot(t, L'); grid on; xlabel('时间 (s)'); ylabel('杆长 (m)');

轨迹参数需要说明:z 方向振幅 0.05 m、频率 2 rad/s,位置平移振幅 0.02 m 和 0.01 m,姿态角振幅从 2° 到 5° 不等。这些数值对应一个重 10 kg 级别、工作空间半径约 0.3 m 的典型研究平台,不是标准值,只是一个安全起点。如果手头平台尺寸不同,可以把幅度缩放到动平台半径的 10% 到 20%,这个比例一般能保证六根杆的杆长曲线始终在极限长度内平滑变化。

如果 plot 出来的曲线有几条出现台阶或不连续跳变,往往是欧拉角溢出或铰链点坐标顺序不对,先不要往控制方向排查。杆长连续只是必要条件,还不能说明逆解正确,因为姿态矩阵可能描绘了一个错误的姿态,但对同一段轨迹仍保持连续。

3.2 铰链配对顺序:最容易让逆解“看起来对”的隐藏问题

逆解代码中Ab和At0的第 i 列,必须按“第 i 根杆连接静铰链 i 与动铰链 i”的配对关系排列。如果平台 CAD 模型里的编号顺序和代码不一致,逆解不会报错,但杆长曲线会表现出奇怪的高频振荡,因为每一根杆都在试图连接错误的两个端点。

我通常把检查步骤写成一个独立函数,将六对铰链一次画出连线:

% 可视化六杆连接关系: 检查铰链配对是否合理 function plot_platform(Ab, At0, pos, R) At_w = pos(:) + R * At0; % 所有动铰链转到静系 figure; hold on; axis equal; scatter3(Ab(1,:), Ab(2,:), Ab(3,:), 'ro', 'filled'); scatter3(At_w(1,:), At_w(2,:), At_w(3,:), 'bo', 'filled'); for i = 1:6 plot3([Ab(1,i) At_w(1,i)], [Ab(2,i) At_w(2,i)], ... [Ab(3,i) At_w(3,i)], 'k-'); end end

这个函数把静铰链画成红点、动铰链画成蓝点、六根杆画成黑线。正确配对下,六条线在两圆之间均匀展开,不会出现交叉或某根杆横穿平台中心的情况。这里有一条血泪经验:某次在动平台上顺时针编号,静平台上逆时针编号,逆解算出的杆长曲线看起来完全正常,直到伺服仿真阶段发现平台无法归零,排查半天才发现是配对顺序写反。

造成问题的原因是逆解只用范数,不感知方向。范数对两个端点的顺序对称,静铰链和动铰链的编号顺序却没有对称性,所以配对检查必须依赖可视化而非数值。

3.3 逆解结果如何作为 Simulink 的输入准备

轨迹预计算完成后,可以把L保存成时间序列,作为 Simulink 模型驱动的参考信号。常见做法有两种:一是直接把离散时间点写入timeseries对象;二是在 Simulink 的 MATLAB Function 块里实时调用逆解。前者适合先验证机械结构是否能承受轨迹,后者适合把逆解和控制闭环放在一起仿真。

我一般先走第一种,因为timeseries可以直接连接 Simulink 的信号线,省去函数块内部的变量声明。代码用L_ts = timeseries(L', t);一行即可生成。注意L的每一列对应一个时刻的六杆长度,所以先转置再传给 timeseries。

到这里,运动学部分的验证闭环已经闭合。是否还需要做动力学,取决于目标:如果只是算杆长极限和干涉检查,逆解已经足够;如果要评估电机选型、验证重力补偿和抗扰动能力,就必须进入动力学仿真。

4. 动力学仿真:从速度雅可比矩阵到闭环等效质量模型

4.1 雅可比矩阵:把平台速度映射到六根杆的伸缩速度

动力学的第一步,是要先建立连杆伸缩速度和动平台速度之间的关系。逆解给了长度,但没有给速度映射;动力学仿真里速度映射由速度雅可比矩阵承担。对每条腿 i,单位方向向量 u_i 和动铰链相对质心的位置矢量 r_i 定义如下:

u_i = (At_i - Ab_i) / L_i

r_i = At_i - pos

腿长变化率等于平台质心平移速度 v_p 在 u_i 上的投影,加上角速度 ω 引起的分量,合起来就是:

Ldot_i = u_i^T v_p + (r_i × u_i)^T ω

把六条腿的表达式堆叠起来,就得到一个 6 行 6 列的雅可比矩阵 J,其中第 i 行为[u_i', (r_i × u_i)']。代码实现如下:

% 速度雅可比矩阵: 平台速度 <-> 六杆伸缩速度 % 输入: pos(3x1), R(3x3), Ab(3x6), At0(3x6) % 输出: J(6x6), L(6x1) function [J, L] = stewart_jacobian(pos, R, Ab, At0) J = zeros(6,6); L = zeros(6,1); for i = 1:6 At_w = pos(:) + R * At0(:,i); % 动铰链在静系坐标 leg = At_w - Ab(:,i); % 杆长矢量 L(i) = norm(leg); u = leg / L(i); % 单位方向向量 r = At_w - pos(:); % 动铰链相对质心矢量 J(i, 1:3) = u'; J(i, 4:6) = cross(r, u)'; end end

需要注意,这是一个与位形相关的矩阵,平台每换一个姿态,J 就变一次。在导入动力学模型之前,先检查 J 的条件数cond(J)。条件数过大说明当前姿态接近奇异,此时逆解可能仍得到合法杆长,但速度层面的放大关系失真,动力学结果不可信。

提示:奇异位形附近最容易出现“逆解正常但动力学发散”的现象。条件数曲线应当和杆长曲线一起绘制,作为第一道筛选防线。

4.2 动力学方程:工程上够用的等效质量模型

对一个以研究为主要目的的项目,我通常不直接展开六维拉格朗日方程,而是用一个“平台质量主导、六条腿质量折算”的动力学模型。假设动平台质量 m、惯性张量 I,平台广义坐标取位置和 ZYX 欧拉角,那么动平台产生的广义力与平台加速度的关系可以写成:

tau_cart = M_cart * qddot + g(q)

其中 M_cart 是 6x6 质量矩阵,由平动质量和转动惯量组成,g(q) 是重力项。这个公式的前提是各条腿的质量远小于平台质量,或在控制中已做前馈补偿。对比完整拉格朗日方程,该模型不包含科氏力,只在低速小姿态角范围内足够准;但作为仿真研究已经能回答“平台在这个轨迹下需要多大驱动力”这类问题。

MATLAB 函数写出来大致是这样:

% 平台等效动力学模型 % 输入: pos(3x1), eul(3x1), qddot(6x1), param(结构体) % 输出: tau_cart(6x1) 广义力 function tau_cart = platform_dyn(pos, eul, qddot, param) R = eul2rotm(eul, 'ZYX'); % 质量矩阵: 前三行为平动质量, 后三行为转动惯量 M_cart = diag([param.m, param.m, param.m, ... param.Ix, param.Iy, param.Iz]); % 重力广义力(简化: 纯力, 不含姿态耦合) grav = [0; 0; -param.m * param.g; 0; 0; 0]; tau_cart = M_cart * qddot + grav; end

参数结构体param中可以包含m=12、g=9.81、Ix=0.35、Iy=0.35、Iz=0.6。这里重力项是近似处理,把重力当作作用于质心的纯力;更高精度模型应当把重力在质心上的分量通过旋转矩阵映射到广义坐标,但姿态角在 10° 以内时误差很小。

真正的广义驱动力还要通过雅可比矩阵变换成六根杆上的力:

tau_leg = J^T * tau_cart

tau_leg 才是电机需要提供的驱动杆轴向力。这条映射是对称的,它由逆运动学和速度雅可比推导,不需要额外的力学分析,前提是忽略杆件惯量、摩擦和铰链柔性。

4.3 在 Simulink 里搭一个闭环仿真骨架

动力学仿真我不建议把整个模型都塞进一个脚本里迭代,那样排错效率低。常见做法是用 Simulink 搭三层结构:

  • 第一层是参考轨迹源,用 Signal Builder 或前面的 timeseries 输出位姿和角速度;
  • 第二层是逆解与雅可比计算,用 MATLAB Function 块承载;
  • 第三层是动力学模型,把 tau_leg 作为输入,经过积分得到平台位姿和速度,再反馈回第二层。

最简洁的闭环是:期望杆长 L_des 由轨迹逆解得到,将实际杆长反馈做差,经过一个比例增益 Kp 后换算为驱动器力,同时把动力学模型输出的加速度积分成速度、再积分成位姿。这个仿真不涉及电机换向和摩擦,属于理想闭环,但能把运动学、动力学、控制三者的数据流通顺验证一遍。

具体搭建步骤是:新建空白模型,添加两个 MATLAB Function 块。第一个块输入时间 t、输出 L_des 和 J,内部调用stewart_inv与stewart_jacobian;第二个块输入当前位姿与 qddot,输出 tau_cart。两个积分器按“加速度—速度—位置”级联。运行前在模型配置里把求解器设为 ode45,最大步长设为 0.001 秒。

最大步长是动力学仿真的关键参数。步长过大,关节力和加速度会出现数值振荡;步长过小,仿真速度被拖慢。0.001 秒对 10 秒仿真通常是一个平衡点。仿真结束后,把六根杆的杆长和平台位置信号都保存下来,与期望曲线对比。如果跟踪误差在姿态角不超过 10° 的轨迹段上小于初始杆长的 1%,基本可以认为运动学和动力学模型没有原则性错误。

5. 避坑常见问题与排查: Stewart 平台仿真中的 5 个真实踩坑记录

5.1 现象:六根杆长曲线出现突兀跳变,位置却看起来正常

原因:旋转顺序不一致或角度单位混用。最常见的错误是直接把角度传感器传来的值当弧度送进eul2rotm。比如横滚 30° 被当成 30 弧度,旋转矩阵在几个采样点后翻转,杆长自然剧烈变化。

解决:在数据入口处统一成弧度制,并在调用eul2rotm前加一行断言:assert(all(abs(eul) < pi), '欧拉角超出合理范围')。角度超过 ±180° 立即停止,可以拦截大部分单位错误,却几乎不会误杀正常轨迹。

5.2 现象:逆解杆长为零或负值

杆长出现零,说明动平台和静平台的某个铰链点重合了;负数不可能出现在范数计算里,如果有人算出了负数,一定是直接算了矢量分量而没有取模。

真正需要排查的是重合:当某根杆的铰链配对写错,例如把第 3 个静铰链对应到第 3 个动铰链,而这两个点恰好属于同一物理位置时,norm会得到 0。解决方式是先检查两组点是否沿圆周均匀分布,再画连线图查看配对关系。杆长接近零意味着该腿在压缩极限位置,并联平台在这一点附近易发生奇异,仿真应当跳过或标记为不可达。

5.3 现象:Simulink 中加速度和驱动力高频振荡,杆长仍平滑

原因在动力学侧:闭环里用了过大的比例增益,或者积分步长太大。杆长平滑说明逆解正确,振荡说明动力学模型与控制回路之间的数值刚性很强。把最大步长从 0.01 改成 0.001 后若振荡消失,问题在步长;若仍然存在,则要检查 Kp 是否超过平台的等效刚度极限。

解决方式:先调步长,后调增益。我常用的做法是先把 Kp 设为 0,观察模型在重力下是否稳定;再以 100、200、500 逐步递增,找到振荡临界点。这个临界点可以作为后续控制器设计的参考,而不是当作参数搜索的终点。

5.4 现象:静态位姿下平台缓慢漂移,重力补偿看似失效

原因:重力项处理过于粗糙。4.2 节把重力当纯力处理,一旦姿态角增大,重力会产生与姿态耦合的力矩,只补偿力不补偿力矩就会造成漂移。

解决:把重力项改为R * [0;0;-m*g]在广义坐标中的映射。做法是先计算质心在惯性系下的重力矢量,再通过雅可比矩阵的力映射投影到六条腿上。做这一步时坐标系要保持一致,否则会出现“补偿了一半”的玄学现象,程序看起来没错,平台就是不肯停住。

5.5 现象:仿真在某个轨迹点突然算不下去,cond(J) 急剧增大

原因:轨迹规划穿过了奇异位形。Stewart 平台的工作空间不是球体,而是六条腿长度和关节角共同约束下的不规则体。轨迹规划如果只用位置和姿态的正弦叠加,很容易扫过奇异区域。

解决:动力学仿真前对整段轨迹逐点计算cond(J),并绘制条件数曲线。条件数超过临界值(我一般取 50)的区间直接从参考轨迹中剔除,或调整该时间段的运动幅度。之后再跑动力学就不会出现矩阵接近奇异导致的数值发散。这个检查放在逆解之后、动力学之前,是成本最低的一道防线。

6. 进阶:用闭环跟踪误差验证逆解和动力学模型到底准不准

模型建完、参数调完,最后还需要一个能拿得出手的验证指标。我推荐做“闭环轨迹跟踪误差”实验:给平台一个同时包含位置和姿态变化的参考轨迹,用上文的理想闭环跑完,比较实际位姿和参考位姿。

实现上可以先跑一次开环逆解,把六杆长度和雅可比矩阵存成时间序列。然后在 Simulink 里让实际平台从初始位姿出发,通过杆长反馈跟踪这条序列,结束后在 MATLAB 里计算误差:

% 计算闭环跟踪误差 pos_err = sqrt(sum((pos_actual - pos_ref).^2, 2)); eul_err = sqrt(sum((eul_actual - eul_ref).^2, 2)); disp(['位置最大误差: ', num2str(max(pos_err))]); disp(['姿态最大误差: ', num2str(max(eul_err))]);

pos_err的 2 范数应始终小于参考轨迹振幅的 2%,姿态角误差应小于初始姿态幅度的 5%。如果误差超限,优先检查重力补偿方向,其次是增益参数。这里有一个我自己的习惯:跑验证前先把所有角度都减去初始偏移,否则零位偏差会被误判为跟踪误差。

如果位置误差很小、姿态误差却很大,问题通常出在动平台惯性张量上。转动惯量如果是估算的,误差会导致姿态通道刚度不足。修正办法是用一个简单实验:只给偏航一个正弦激励,记录偏航响应,反推惯性张量中的对应项。这样标定出来的惯量比直接查表可靠。

整套流程走到这里,你已经有一个从逆解、雅可比、动力学到闭环验证的完整链路。这个方向值不值得投入,答案很清楚:如果目标是姿态级控制和驱动力评估,这套基于 MATLAB 的方法足够;如果目标要做到复杂力控制,还需要补充摩擦模型和驱动建模。我的建议是先把这条链路完整跑通,再做扩展。希望这套验证习惯能帮你少走弯路,也希望这些参数与坑位记录对你有实际帮助。

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

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询