简介面向从事摩擦建模、伺服控制与机械系统仿真的Matlab学习者这套代码围绕Stribeck摩擦模型展开重点复现低速阶段的负斜率摩擦现象克服静摩擦后摩擦力随速度先降后升并以指数形式描述该过程。内容适合需要理解非线性摩擦特性、开展控制补偿或完成课程实验的初、中级读者使用。包内共3个.m文件压缩包约1KB包含模型实现、演示脚本及辅助函数便于直接运行观察参数变化及δ、vs等对曲线形态的影响也方便嵌入已有仿真流程。已有1138人学习。读者可据此理解静摩擦、库仑摩擦与Stribeck速度尺度等参数的含义掌握从公式到代码的映射并通过示例快速验证负斜率段与速度拐点为控制算法补偿和论文复现提供可改可扩的基础脚本。1. 低速爬行不是机械问题Stribeck 摩擦力模型在 Matlab 里的一次完整落地精密进给台做 0.5 mm/s 匀速跟随速度曲线却走成台阶位置环一抖一抖地呼吸换成更贵的导轨还是原样。遇到这种低速爬行很多人的第一反应是查机械装配和伺服增益但真正的原因常常在摩擦力本身它在零速附近不是常数而是一条先陡降、再回升的曲线。Stribeck 摩擦力模型描述的就是这段曲线而这组 Matlab 实现代码与使用示例想干的事很直接——把这条曲线变成一个能算、能拟合、能接进 Simulink 的解析函数让你的仿真里出现和台架上一样的爬行也让前馈补偿有参数可调。它适合做运动控制仿真的工程师、做摩擦补偿算法验证的研究生以及手上有一批速度-摩擦力实测数据但不知道怎么反推参数的人。2. 拆开 Stribeck 曲线五个参数各自的物理含义与选型理由动手写函数之前先把模型本身的骨架立住。不然后面拟合出来的数字看着能过一放进闭环就会翻车。2.1 从边界润滑到流体润滑负斜率段是哪来的两个接触面在相对运动时润滑油膜的状态随速度变化分三段。速度极低时处于边界润滑微凸体直接接触摩擦力最大这就是最大静摩擦力Fs速度稍微上来一点油膜开始建立、微凸体接触面积下降摩擦力反而随速度升高而下降这一段就是 Stribeck 负斜率段也是低速爬行的直接来源速度继续升高流体动压效应占据主导粘性阻力随速度线性上升斜率变正。稳态下把这三段写成一个关于速度v的解析式最常用的形式是F(v) Fc·sgn(v) (Fs - Fc)·exp(-|v/vs|^d)·sgn(v) Fv·v五个参数含义如下Fc是库仑摩擦力也就是负斜率段走完之后剩下的那部分近似恒定的摩擦力Fs是最大静摩擦力必须满足Fs ≥ FcFv是粘性阻尼系数单位 N·s/m描述高速段的线性上升vs是 Stribeck 特征速度量纲上就是速度它决定负斜率段在哪个速度量级上塌下去d是形状指数控制从Fs衰减到Fc的曲线陡峭程度工程上常见取值在 0.5 到 2 之间取 1 的指数衰减和取 2 的高斯型衰减在低速段差别很大。这里有个容易忽略的细节vs和d是强相关的。同一批数据用vs0.01, d1和vs0.014, d1.6可能拟合出几乎一样的残差平方和但把速度外推到 0.05 mm/s 时两条曲线会分叉。这也是后面参数辨识章节要反复强调固定d或者做多初值搜索的原因。2.2 五种常见写法怎么选指数型、tanh 平滑型、LuGre 与 Karnopp市面上叫Stribeck 模型的实现其实不止一种选错了会在特定场景下白干。经典指数型就是上面那个式子sgn函数在零速处跳变。优点是形式简单、参数物理意义清楚、拟合稳定缺点是零速不连续直接丢进 ode45 会在零速附近反复穿越、步长被压到极小仿真时间成倍增加。适合稳态摩擦特性绘制、参数辨识、不关心零速细节的场合。tanh 平滑型把sgn(v)换成tanh(v/v_eps)v_eps是平滑带宽。好处是全域连续可导数值求解器不会在零速处卡死代价是丢掉了真正的静摩擦——零速时输出为 0 而不是Fs。做带位置保持的伺服仿真时这一点必须注意否则你的系统会自己滑走。Karnopp 模型走的是另一条路设定一个零速判定带|v| v_0带内直接按静摩擦处理用一个力的平衡判断是否起步。它能把静摩擦和粘滑切换都表达出来但要处理状态判断写起来啰嗦且v_0选得不好会引入非物理的极限环。LuGre 动态模型在 Stribeck 稳态曲线的基础上加了一个内部状态z描述鬃毛平均变形能复现预滑移阶段的滞回和微位移下的非线性刚度。如果研究对象是微米级定位、需要关注换向时的位移滞后只靠稳态 Stribeck 是不够的必须上 LuGre 或 GMS。LuGre 的代价是多了两个参数鬃毛刚度sigma0、微观阻尼sigma1辨识难度和对初值的敏感性都上一个台阶。多项式型在零速附近用低阶多项式拼接好处是求导方便、便于做控制器设计时的线性化坏处是拼接点不光滑、参数没有物理含义。一份实用的选型对照使用场景推荐写法关键理由绘制摩擦-速度特性曲线、教学演示经典指数型参数直观曲线形状可控闭环伺服仿真需要长时间积分tanh 平滑型零速连续求解器不易卡步位置保持、粘滑振动复现Karnopp 或 tanh 静摩擦补偿项需要明确的静摩擦区间微米级定位、换向滞后分析LuGre / GMS需要预滑移滞回摩擦前馈补偿的工程实现经典指数型离线算表在线计算量小可查表插值多数工程仿真里我一般会准备两个版本辨识和画图用经典指数型保证参数干净闭环仿真用 tanh 平滑型保证数值稳定两者的Fs/Fc/Fv/vs/d完全共用一套。2.3 参数从哪来三种获取路径和它们的可信度Fs、Fc、Fv、vs、d五个参数来源无非三条路。第一条是查样本和估算。导轨、丝杠的厂商样本一般会给动摩擦系数和静摩擦系数乘上法向载荷就得到Fc和Fs的量级Fv按润滑脂粘性粗估vs按经验取 0.01~0.05 m/sd先取 1。这条路的好处是十分钟就能让仿真跑起来坏处是误差可能到 50% 以上只能用来验证控制结构是否合理不能用来下补偿后能到多少精度的结论。第二条是准静态扫频实验。让执行机构以一系列恒定速度单向运动每个速度点稳定后采集驱动力减去惯性力和已知外载剩下的就是摩擦力。速度点建议在低速段密、高速段疏比如 0.01、0.02、0.05、0.1、0.2、0.5、1、2、5、10 mm/s 各采 2~3 秒取均值。正反两个方向分别采一遍能顺带看出方向依赖的非对称性。第三条是闭环辨识。给定控制器让系统做正弦或扫频跟随用输入输出数据通过最小二乘或智能优化算法反推参数。这条路的难点在于摩擦力和控制器输出、惯性力完全耦合在一起必须先有一个足够准的动力学模型否则辨识出来的摩擦参数会把模型误差一起吸收进去。三条路的可信度是递增的成本也是递增的。实际项目里比较稳妥的做法是先用第一条路给出一组初值用第二条路的数据做最小二乘拟合把拟合残差控制在满量程摩擦力的 5% 以内再往下走。3. Matlab 实现一个可复用函数、一组默认参数、两种调用方式代码部分按能直接抄的标准写。整个实现拆成三个文件参数默认值函数、摩擦计算函数、被控对象导数函数。这样拆的好处是参数只在一个地方维护仿真脚本、辨识脚本、Simulink 模块都能共用。3.1 stribeck_friction.m核心计算函数function F stribeck_friction(v, p) %STRIBECK_FRICTION 稳态 Stribeck 摩擦模型 % F stribeck_friction(v, p) % v : 相对速度标量或任意维数组单位 m/s % p : 参数结构体字段 Fs/Fc/Fv/vs/d/v_eps % F : 与 v 同尺寸的摩擦力单位 N % % F(v) Fc*s(v) (Fs-Fc)*exp(-|v/vs|^d)*s(v) Fv*v % 其中 s(v) tanh(v/v_eps)v_eps 为 0 时退化为 sign(v) if nargin 2 || isempty(p) p stribeck_default_params(); end v double(v); av abs(v); % 零速邻域平滑v_eps0 用 tanh 过渡否则用严格 sign if isfield(p, v_eps) p.v_eps 0 s tanh(v ./ p.v_eps); else s sign(v); end % 逐元素运算保证 v 为矩阵时也能正确广播 stribeck (p.Fs - p.Fc) .* exp(-(av ./ p.vs) .^ p.d) .* s; F p.Fc .* s stribeck p.Fv .* v; end逻辑说明先做参数缺省判断再算|v|和符号项。符号项是唯一有分支的地方v_eps 0走tanh、否则退回sign这样同一份代码既能画教科书上的理想曲线v_eps 0也能直接喂给 ode 求解器v_eps 1e-4。三项分别对应库仑项、Stribeck 衰减项和粘性项全部用.运算符保证向量化v传 2001 点的等距向量时一次算完画图不用循环。参数说明v必须是 SI 单位 m/s如果实验数据是 mm/s调用前统一除以 1000。p.vs的量纲和v一致p.d无量纲p.Fv是 N·s/m。v_eps的取值建议是仿真中最小关注速度的 1/10 到 1/100太小起不到平滑作用太大会把真实的低速非线性抹平一般 1e-4 到 1e-3 之间。3.2 stribeck_default_params.m把参数收在一处function p stribeck_default_params() %STRIBECK_DEFAULT_PARAMS 一组中等负载滑动导轨的量级参数 % 数值仅作仿真起点实际项目必须用实测数据替换 p.Fs 1.20; % 最大静摩擦力N p.Fc 0.85; % 库仑摩擦力N p.Fv 0.060; % 粘性阻尼系数N·s/m p.vs 0.012; % Stribeck 特征速度m/s p.d 1.0; % 形状指数建议 0.5~2 p.v_eps 1e-4; % 零速平滑带宽m/s置 0 退化为 sign end这里的数值是给中等负载滑动导轨的一个量级参考Fs/Fc比值约 1.4属于偏干的接触状态润滑充分时这个比值会更接近 1.05~1.15。每次改参数只改这个文件其他脚本用p stribeck_default_params(); p.vs 0.02;这种覆盖方式避免同一批数字散落在十几个文件里——这一点吃过亏后面在避坑章节还会提。3.3 画曲线验证形状三分钟确认参数没填反p stribeck_default_params(); % 低速段加密才能看清负斜率 v [-flip(linspace(0, 0.05, 1000)), linspace(0, 0.05, 1000)]; F stribeck_friction(v, p); figure(Color, w); plot(v * 1000, F, LineWidth, 1.6); hold on; yline(p.Fc, --, Fc); yline(-p.Fc, --, -Fc); xline(0, :, Color, [0.5 0.5 0.5]); grid on; box off; xlabel(速度 v (mm/s)); ylabel(摩擦力 F (N)); title(Stribeck 摩擦特性曲线);画出来应该看到原点附近跳变到Fs随速度绝对值增加迅速下凹到接近Fc然后被粘性项缓慢抬起整体关于原点奇对称。如果曲线是单调上升或者原点附近没有尖峰多半是Fs Fc或者v_eps给大了。这一步花三分钟能挡掉后面几小时的困惑。3.4 接进动力学ode15s 求解质量-弹簧-摩擦系统function dx plant(t, x, p, c) %PLANT 单自由度质量-弹簧-阻尼 Stribeck 摩擦 % x(1)位移 mx(2)速度 m/s % c.m 质量 kgc.c 阻尼 N·s/mc.k 刚度 N/mc.F 驱动力 N Ff stribeck_friction(x(2), p); dx [ x(2); (c.F - c.k * x(1) - c.c * x(2) - Ff) / c.m ]; endp stribeck_default_params(); c struct(m, 0.5, c, 2.0, k, 120, F, 1.0); opts odeset(RelTol, 1e-8, AbsTol, 1e-10, MaxStep, 1e-3); [t, x] ode15s((t, x) plant(t, x, p, c), [0 2.0], [0; 0], opts); figure(Color, w); subplot(2,1,1); plot(t, x(:,1)*1000); ylabel(位移 (mm)); grid on; subplot(2,1,2); plot(t, x(:,2)*1000); ylabel(速度 (mm/s)); xlabel(时间 (s)); grid on;逻辑说明plant把摩擦力直接作为状态导数的一部分没有引入额外的内部状态所以这是纯稳态模型。求解器选ode15s而不是ode45因为 tanh 平滑后系统在零速附近仍然接近刚性——摩擦力对速度的等效刚度可达Fs/v_eps量级ode45会把步长压到 1e-9 以下仿真几秒钟要跑几分钟。参数上MaxStep限制在 1e-3 是为了不错过粘滑切换的细节RelTol放到 1e-8 是因为位移量级只有毫米、速度量级只有毫米每秒默认的 1e-3 相对误差会吃掉有效数字。驱动力取 1.0 N 略大于Fc而小于Fs是为了让系统从静止起步后进入低速滑移从而把爬行现象逼出来。3.5 Simulink 里的 MATLAB Function 模块写法模型结构不复杂常值力源 → Sum → Gain(1/m) → Integrator → Integrator → 速度反馈回摩擦模块。摩擦模块用 MATLAB Function 模块内容如下function F fcn(v, Fs, Fc, Fv, vs, d) % 代码生成友好的写法全部用标量输入不用结构体和 isfield av abs(v); s tanh(v / 1e-4); F Fc * s (Fs - Fc) * exp(-(av / vs) ^ d) * s Fv * v; end三个必须注意的点。一是不要在这个模块里用结构体MATLAB Function 模块在代码生成时对结构体字段的动态访问支持有限isfield这类判断在部分配置下会直接报错全部展开成标量输入最稳。二是参数不要写死在函数体里在基础工作区建Simulink.Parameter对象再连到模块输入端改参数不用重新编译模型。三是零速平滑带宽1e-4是硬编码的如果模型用定步长求解器比如ode1步长 1e-3这个带宽比步长还小等于没平滑此时要么把带宽放大到 1e-3 量级要么改用变步长求解器。4. 参数辨识实战用 lsqcurvefit 从实测数据反推五个参数有了模型和函数下一步是让五个参数从数据里长出来而不是拍脑袋。4.1 数据预处理把台架数据整理成速度和摩擦力两张表前提是有准静态匀速实验数据。每个速度点采一段稳定区间取均值和标准差剔除标准差超过均值 10% 的坏点。正反行程分开存先用正向数据拟合再用反向数据验证如果两向的Fc差超过 15%说明机构存在明显的方向依赖这时要么做成非对称模型Fc、Fc-分开要么回去查装配。数据处理脚本raw readmatrix(friction_sweep.csv); % 第1列速度 mm/s第2列力 N v_mm raw(:,1); F_N raw(:,2); % 剔除反向行程与异常点 mask v_mm 0 abs(F_N) 10; v v_mm(mask) / 1000; % 统一到 m/s F F_N(mask); % 按速度排序便于后面画残差 [v, idx] sort(v); F F(idx);单位统一是这一步唯一但最致命的动作。速度用 mm/s、力用 N 混着进拟合最后出来的Fv会差三个数量级模型在闭环里表现完全不对。4.2 lsqcurvefit 拟合初值、边界与多起点% 待拟合参数向量 x [Fs, Fc, Fv, vs, d] model (x, v) x(2) .* tanh(v/1e-4) ... (x(1) - x(2)) .* exp(-(abs(v)./x(4)).^x(5)) .* tanh(v/1e-4) ... x(3) .* v; x0 [max(F)*1.05, min(F)*0.95, 0.05, 0.010, 1.0]; lb [0, 0, 0, 1e-4, 0.5]; ub [inf, inf, inf, 1.0, 2.0]; opts optimoptions(lsqcurvefit, ... Display, iter, ... Algorithm, trust-region-reflective, ... MaxFunctionEvaluations, 3000, ... FunctionTolerance, 1e-12, ... StepTolerance, 1e-10); [xopt, resnorm, residual] lsqcurvefit(model, x0, v, F, lb, ub, opts); fprintf(Fs%.4f N Fc%.4f N Fv%.5f Ns/m vs%.5f m/s d%.3f\n, xopt); fprintf(RMSE %.4f N\n, sqrt(mean(residual.^2)));逻辑说明model是匿名函数形式的模型参数顺序和x0严格对应这里最容易错位的是把Fv和vs的位置搞反——错位后拟合照样收敛但参数毫无意义一定要打印出来和量级核对。x0的构造有讲究Fs取实测力的最大值略放大Fc取最小值略缩小这两个初值基本靠谱vs的初值用低速段中摩擦力衰减到(FsFc)/2对应的速度估计懒得估就取 0.01d固定从 1.0 起步。边界条件的作用Fs ≥ Fc这个物理约束没法直接写成lb/ub要靠初值和后面的校验保证d限制在 [0.5, 2] 是防止拟合出 5 以上的形状指数那种曲线在拟合区间内贴合、外推时完全失真。trust-region-reflective算法要求残差个数不少于参数个数这里十几个速度点对五个参数满足条件。如果数据点少于 8 个改用levenberg-marquardt但要接受它更容易陷局部极小。单次拟合很可能停在局部极小尤其是vs和d的耦合方向。稳妥做法是撒一组初值做多起点best struct(rmse, inf); for vs0 [0.003 0.008 0.015 0.030] for d0 [0.6 1.0 1.5] try [x, ~, r] lsqcurvefit(model, [max(F)*1.05, min(F)*0.95, 0.05, vs0, d0], ... v, F, lb, ub, opts); rmse sqrt(mean(r.^2)); if rmse best.rmse best struct(rmse, rmse, x, x); end catch continue; % 某些初值组合会发散直接跳过 end end end xopt best.x;十二次拟合通常十几秒跑完比手工调初值省时间。try-catch不能省边界附近偶尔会抛Local minimum possible之外的错误一次异常不该中断整个搜索。4.3 结果验证残差图比 RMSE 更能说明问题Ffit model(xopt, v); figure(Color,w); subplot(2,1,1); semilogx(v*1000, F, o, v*1000, Ffit, -, LineWidth, 1.4); grid on; xlabel(速度 (mm/s)); ylabel(摩擦力 (N)); legend(实测,拟合,Location,southeast); subplot(2,1,2); semilogx(v*1000, F - Ffit, o-); grid on; yline(0, --); xlabel(速度 (mm/s)); ylabel(残差 (N));看三件事。第一残差是否随速度出现系统性趋势如果低速段残差普遍为正、高速段普遍为负说明形状指数d没选对曲线弯的程度不匹配。第二残差绝对值是否小于满量程摩擦力的 5%超了就别急着往下做补偿。第三拟合出的vs是否落在实验速度点的覆盖范围内如果vs落在了最小测试速度以下说明数据对负斜率段根本没有约束这个参数是编出来的闭环里会出问题。还有一步常被跳过拿反向行程数据当验证集不参与拟合只看预测误差。两边误差都在可接受范围内参数才算可信。5. 避坑与排查五个让仿真结果失真的细节这一节全是踩出来的按现象 → 原因 → 解决组织。5.1 现象仿真跑不动oder45 步长被压到 1e-9原因用了严格sign版本的摩擦函数系统在零速附近反复穿越摩擦力在Fc和-Fc之间跳变求解器为了满足误差要求不断缩步。这在闭环位置保持仿真里几乎必然发生。解决把v_eps从 0 改成 1e-4 量级或直接换ode15s并把MaxStep限制在 1e-3。如果必须保留理想静摩擦用 Karnopp 的零速判定带替代sign带宽取求解器步长对应的速度变化量的 2~3 倍。5.2 现象辨识残差很小但换一组数据参数全变原因vs和d在拟合区间内高度相关数据点又集中在高速段低速负斜率段只有两三个点参数不可辨识。残差小是因为模型在高速段被粘性项主导低速段的形状怎么变残差都差不多。解决补低速实验点在vs附近至少布 5 个点拟合时固定d 1.0只放Fs/Fc/Fv/vs四个参数稳定后再考虑放开d或者对vs加一个弱的正则项把解往物理合理区间拉。5.3 现象Simulink 仿真结果和 m 脚本对不上原因多半是单位或者求解器设置不一致。脚本里v是 m/sSimulink 里某个 Gain 忘了换算脚本用ode15s变步长Simulink 用定步长ode1且步长 1e-3而v_eps只有 1e-4平滑完全失效。解决在 Simulink 模型里显式标注每个信号的单位速度信号进摩擦模块之前加一个单位换算的注释块把v_eps设成不小于求解器步长对应速度变化的量级定步长 1e-3 时v_eps至少给 1e-3。5.4 现象位置保持工况下系统缓慢漂移原因用了 tanh 平滑版本零速时tanh(0) 0摩擦力输出为 0静摩擦被彻底抹掉了。系统在零速附近没有任何阻力矩积分器上的微小残差会让位置慢慢滑走。解决位置保持场景不能只用平滑版。要么改成 Karnopp 结构在|v| v_0时按静摩擦做力平衡判断要么在平滑版基础上叠加一个零速保持项当|v| v_0且外力小于Fs时输出等于外力反向。这一点没得绕。5.5 现象模型在台架上验证通过换一台同型号设备就失效原因摩擦参数是装配预紧、润滑状态、表面粗糙度的函数同型号设备的Fs差异可以到 30%。把一台设备辨识出的参数当型号参数用必然翻车。解决参数按单机管理每台设备建一份参数文件标定流程固化成脚本控制算法里给Fs/Fc留 20%~30% 的在线自适应余量或者用低速段的跟踪误差做慢速在线修正。别指望一套参数通吃。6. 进阶用法把摩擦模型接进前馈补偿附一份验证清单模型算得准只是前提真正产生价值的是拿它去做补偿。思路很直接既然F(v)能算那就在控制量里预先加上同样大小的力把摩擦顶掉。function u control_with_friction_ff(v_ref, u_pid, p, alpha) %CONTROL_WITH_FRICTION_FF 带摩擦前馈补偿的控制量 % v_ref : 参考速度m/s % u_pid : 反馈控制器输出N 或 A % alpha : 补偿系数0~1工程上常取 0.6~0.8 F_ff stribeck_friction(v_ref, p); u u_pid alpha * F_ff; endalpha不取 1 是有原因的。前馈力完全等于模型输出时一旦参数偏大系统会在换向瞬间被推过头产生反向的超调取 0.6~0.8 相当于留了一条退路让反馈环还有活干。另外v_ref在零速附近符号会抖直接用会让前馈力来回跳常见做法是给v_ref加一个一阶低通截止频率 20~50 Hz或者在前馈通道里对F_ff做速率限幅。验证环节比补偿本身更容易被糊弄过去。我习惯按下面这张表逐项走每项都有明确的判据不看感觉好多了。验证项做法通过判据常见异常曲线形状画v ∈ [-50, 50]mm/s 的 F-v 曲线原点有尖峰低速有负斜率高速线性上升Fs Fc或曲线单调参数物理性打印五个参数并与样本量级比对Fs/Fc在 1.05~1.5vs在 1~50 mm/svs小于最小测试速度拟合残差残差-速度图无系统趋势绝对值 满量程 5%低速正、高速负 →d偏小开环爬行复现恒力驱动力介于Fc与Fs之间速度曲线呈台阶状台阶周期与台架同量级完全平滑 → 静摩擦被抹掉补偿前后对比0.5 mm/s 匀速跟随比较位置误差峰峰值补偿后峰峰值下降 40% 以上下降不足 →alpha太小或vs偏大换向死区正弦 0.1 Hz 跟随看过零附近误差死区宽度缩小无新增超调出现反向超调 →alpha偏大表里最容易出问题的是开环爬行复现这一项。很多人的模型参数是从样本上抄的曲线画出来漂亮但一跑恒力驱动就是平滑的指数趋近根本没有台阶——这说明Fs/Fc比值给得太小或者v_eps太大把负斜率段抹平了。爬行复现不出来后面的补偿效果就无从对比整个验证链条断在第一环。还有一个我后来固定加上的动作把辨识出的参数、拟合残差、验证曲线一起写进一个friction_calib_YYYYMMDD.mat和当次的实验数据放在同一个目录。原因是三个月后回头查当时为什么把vs定成 0.018光看代码里的数字根本想不起来有这份存档才能复盘。同一个项目里换过三次参数、每次都没留档最后只能重新做一遍扫频实验那个下午的教训记到现在。从那以后我每次改摩擦参数都强制走一遍拟合 → 残差图 → 爬行复现 → 存档这四步一步不跳。希望帮到你。本文还有配套的精品资源点击获取