简介:面向需要在 Matlab 中复现低速负斜率摩擦特性的学习者与工程人员,这份代码包围绕 Stribeck 摩擦力模型展开,重点解决静摩擦向动摩擦过渡阶段摩擦力随速度先降后升的建模问题,并帮助理解静摩擦、库仑摩擦与粘滞摩擦之间的衔接关系。压缩包共 3 个文件,均为 .m 脚本,整体约 1KB,体量轻巧,包含模型实现与调用示例,便于直接放入现有仿真流程中调试和二次修改。已有 1138 人学习,说明该模型在摩擦补偿、伺服控制与机械传动仿真等场景中具有较稳定的参考需求。读者可获得指数型 Stribeck 曲线的参数化写法、负斜率段的速度分区处理思路,以及从示例脚本反推参数含义与曲线形态的完整线索,适合作为课程实验、算法验证和论文复现的基础素材。
1. 低速段总拟合歪:Stribeck 摩擦力模型到底卡在哪
低速运动控制里,Stribeck 摩擦力模型最容易让人产生错觉:高速段看起来一条直线,到了零速附近却怎么调都拟合不上去。这份 Matlab 实现代码和使用示例,核心就是把 Fc、Fs、vs、δ 和粘性项拆开,让你能在脚本里复现一条可解释的 Stribeck 曲线,而不是只会画 sign(v)。它适合做伺服进给、机器人关节、液压执行器低速摩擦补偿的开发者,也适合正在写论文实验但需要一份能跑通的参考实现的人。先记住一个边界:Stribeck 是稳态速度-摩擦力模型,不是万能的预滑动模型,零速附近的处理方式决定后面仿真会不会抖。
2. 从公式到 m 文件:Stribeck 摩擦力的 Matlab 函数落地
Stribeck 曲线真正难写的部分不是公式,而是把“符号函数在零速处跳变”这件事处理得既能算、又不至于把求解器逼疯。常见写法是把摩擦力拆成三块:库仑项、Stribeck 衰减项和粘性项。写成 Matlab 函数时,最重要的不是行数少,而是参数含义清楚、输入输出维度稳定、零速附近有可控的平滑尺度。
2.1 符号平滑与三个速度区间
稳态 Stribeck 模型一般写成:
Ff(v) = Fc * sgn(v) + (Fs - Fc) * exp(-(|v|/vs)^delta) * sgn(v) + sigma * v其中 Fc 是库仑摩擦力,Fs 是静摩擦附近的峰值,vs 是 Stribeck 速度,delta 控制曲线从 Fs 衰减到 Fc 的形状,sigma 是粘性阻尼系数。v 为正时摩擦力反向阻碍运动,v 为负时符号反过来;真正麻烦的是 v 接近 0 的时候,sgn(v) 在数值上直接跳变,ode45 这类求解器会在零速附近反复缩小步长,表现出来就是速度曲线毛刺、仿真变慢,甚至直接报步长过小。
我一般不会在代码里直接写sign(v),而是用tanh(v/vEps)做平滑符号。vEps 不是物理参数,它只是一个数值正则化尺度,单位跟速度一致。vEps 取太大,低速段的摩擦力会被抹平,拟合出来的 vs 和 delta 会偏;vEps 取太小,零速附近又回到接近跳变的状态,求解器照样难受。经验上先让 vEps 落在 vs 的 1/100 到 1/20 之间,再根据仿真步长和曲线光滑程度微调。如果 vs 是 0.04 m/s,vEps 可以先从 1e-3 或 5e-4 试起,不要一上来就设 1e-8。
三个速度区间要分开看。高速段由 Fc 和 sigma 主导,曲线近似一条斜线;中低速段由(Fs-Fc)*exp(...)主导,会出现明显的负斜率;零速附近则由平滑符号和数值步长共同决定。很多人拟合时只盯高速段,结果 Fc 和 sigma 看起来不错,低速段却完全不对,原因就是 vs 和 delta 没有被低速数据约束住。
| 参数 | 含义 | 单位 | 使用时的注意点 |
|---|---|---|---|
| Fc | 库仑摩擦力 | N | 不能为负,通常由高速段渐近线估计 |
| Fs | 静摩擦峰值 | N | 一般应大于等于 Fc,否则曲线形状会反过来 |
| vs | Stribeck 速度 | m/s | 必须大于 0,决定负斜率出现的位置 |
| delta | 形状指数 | 无量纲 | 常用 1 到 3,越大过渡越陡 |
| sigma | 粘性系数 | N·s/m | 由高速段斜率辅助估计 |
| vEps | 符号平滑速度 | m/s | 不是物理量,只用于数值正则化 |
2.2 函数签名、参数校验和向量化
下面这个函数可以直接保存为stribeckForce.m。它只做一件事:输入速度向量和参数结构体,输出同样大小的摩擦力向量。把参数放进结构体而不是一长串位置参数,是为了后面做辨识和扫参时不容易把顺序搞错。
function Ff = stribeckForce(v, p) %STRIBECKFORCE 计算 Stribeck 稳态摩擦力 % v : 速度向量,单位 m/s % p : 参数结构体,字段为 Fc、Fs、vs、delta、sigma、vEps if nargin < 2 error('需要输入速度 v 和参数结构体 p'); end required = {'Fc','Fs','vs','delta','sigma','vEps'}; for k = 1:numel(required) if ~isfield(p, required{k}) error('参数结构体缺少字段: %s', required{k}); end end if p.vs <= 0 || p.vEps <= 0 || p.delta <= 0 error('vs、vEps、delta 必须大于 0'); end if p.Fc < 0 || p.Fs < 0 || p.sigma < 0 error('Fc、Fs、sigma 不能为负'); end v = v(:); % 统一成列向量,避免维度混用 sv = tanh(v ./ p.vEps); % 平滑符号,替代 sign(v) av = abs(v); shape = exp(-(av ./ p.vs).^p.delta); Ff = p.Fc .* sv ... + (p.Fs - p.Fc) .* shape .* sv ... + p.sigma .* v; end这段代码里有两个地方值得留意。第一,v = v(:)不是多余操作,Matlab 里行向量和列向量混用很容易在.*处报维度错误,统一成列向量后,函数内部和外部数据都好对齐。第二,所有乘法、除法、幂运算都用了点运算,原因是 v 既可能是标量,也可能是几千点的扫速向量;用*和^只适合标量,换到向量就会翻车。
参数校验也不要省。Fs 小于 Fc 时,公式仍然能算,但曲线形状会变成先下后上,和常见 Stribeck 曲线不一致;这种错误在拟合时表现为残差怎么都下不去。代码不一定要强行禁止 Fs < Fc,但至少要在示例脚本里打印参数,看到 Fs 比 Fc 小就先检查数据或初值。
2.3 示例脚本:正弦扫速、阶跃和结果检查
拿到函数后,先别急着接动力学模型,先画一条静态曲线确认形状。下面的示例用一组典型参数生成速度-摩擦力曲线,并检查零速附近有没有异常跳变。这里的数值只是示例,实际使用时必须换成你台架或仿真对象的数据。
% demo_stribeck_curve.m clear; clc; p = struct(); p.Fc = 1.8; % N p.Fs = 3.2; % N p.vs = 0.04; % m/s p.delta = 2.0; % 无量纲 p.sigma = 0.12; % N*s/m p.vEps = 1e-3; % m/s,先取 vs 的 1/40 左右 v = linspace(-0.4, 0.4, 4001); Ff = stribeckForce(v, p); figure; plot(v, Ff, 'LineWidth', 1.2); grid on; xlabel('速度 v (m/s)'); ylabel('摩擦力 F_f (N)'); title('Stribeck 摩擦力曲线'); % 检查零速附近相邻点的最大跳变 idx = abs(v) < 5e-3; jumpMax = max(abs(diff(Ff(idx)))); fprintf('零速附近最大相邻跳变: %.6f N\n', jumpMax); % 检查关键点 fprintf('v=0 时 Ff=%.6f N\n', stribeckForce(0, p)); fprintf('v=vs 时 Ff=%.6f N\n', stribeckForce(p.vs, p)); fprintf('v=0.4 时 Ff=%.6f N\n', stribeckForce(0.4, p));运行后先看三件事。曲线在正负速度上应该近似奇对称;随着 |v| 增大,摩擦力应该从峰值逐步过渡到库仑项加粘性项;零速附近不应该出现肉眼可见的竖直跳变。v=0时输出接近 0 是正常的,因为平滑符号tanh(0)=0,这并不代表真实静摩擦力为零,只代表这个稳态公式在零速点没有给出静摩擦约束。真实系统在零速附近可能处于粘着状态,摩擦力可以在一定范围内取任意值,这部分要靠状态机、事件检测或 LuGre 类模型补,不要指望一条 Stribeck 曲线全包。
如果你要做时间域示例,可以用一个低频正弦速度指令去扫过低速段:
t = linspace(0, 10, 2001); vCmd = 0.08 * sin(2*pi*0.2*t); Fcmd = stribeckForce(vCmd, p); figure; plot(t, vCmd, 'LineWidth', 1.0); hold on; plot(t, Fcmd, 'LineWidth', 1.0); grid on; xlabel('时间 t (s)'); ylabel('速度 / 摩擦力'); legend('v (m/s)', 'F_f (N)');这段代码的意义不是模拟真实运动,而是让你在接 ODE 之前确认函数能接受时变速度输入。真正做参数辨识时,速度指令不能变化太快,否则测到的力包含加速度项和结构振动,拟合出来的 vs 会明显偏小。
3. 参数辨识不靠手调:从扫速数据反推五个 Stribeck 参数
手调参数最大的问题是“看起来像”和“物理上对”之间没有判据。你今天把 vs 调到 0.03,曲线在屏幕上顺眼了,明天换一组数据又得重来。更稳的做法是把 Stribeck 曲线当成一个带约束的非线性最小二乘问题:用一组稳态扫速数据,让 Matlab 去最小化模型输出和实测力的残差,再用残差分布反查模型结构有没有问题。
3.1 稳态扫速数据的预处理
辨识 Stribeck 参数的数据应该是“稳态速度-摩擦力”对,不是随便一段加减速过程。常见做法是让执行器在一组离散速度点上运行,每个速度点保持足够长时间,等速度波动和力信号稳定后取平均。速度点要覆盖低速段和高速段,尤其是 vs 附近要多放几个点,否则 vs 和 delta 的可辨识性很差。
数据进 Matlab 后先做四件事:去掉启动和停止的瞬态段;把速度统一成 m/s,力统一成 N;按速度排序并检查有没有重复点;把正反两个方向的数据分开看。如果正反方向在同一速度下的摩擦力差异很大,先别急着加补偿项,可能是导轨预紧、重力分量或传感器零偏没有处理。基础 Stribeck 模型假设正反方向对称,强行拟合不对称数据会让参数失去物理意义。
% 假设 raw 是导入后的表格,包含 speed 和 force 两列 vData = raw.speed(:); % m/s FData = raw.force(:); % N % 去掉零速点和明显异常点 mask = abs(vData) > 1e-4 & isfinite(vData) & isfinite(FData); vData = vData(mask); FData = FData(mask); % 按速度排序,便于画残差 [vData, idx] = sort(vData); FData = FData(idx); % 简单的三点滑动平均,只用于查看趋势,不替代原始数据拟合 F smooth = movmean(FData, 3);最后一行只是示例,实际代码里不要把变量名写成带空格的形式。更重要的是,滤波要在辨识之前决定好,并且仿真回代时使用同样的滤波设置。很多人拟合时对实测力做了低通,回代时却拿未滤波的仿真去比,残差里混进了滤波器相位差,最后怪模型不对。
3.2 lsqcurvefit 的目标函数与初值边界
为了避免 Fs 小于 Fc,我不直接把[Fc, Fs, vs, delta, sigma]作为优化变量,而是用dF = Fs - Fc参数化,约束dF >= 0。这样优化器不会跑到物理上奇怪的区域,初值也更好给。模型函数单独保存为stribeckModel.m:
function Fhat = stribeckModel(x, v, vEps) % x = [Fc, dF, vs, delta, sigma] Fc = x(1); dF = x(2); vs = x(3); delta = x(4); sigma = x(5); v = v(:); sv = tanh(v ./ vEps); Fhat = Fc .* sv ... + dF .* exp(-(abs(v) ./ vs).^delta) .* sv ... + sigma .* v; end拟合脚本如下。初值不要随便填 1,先从数据里估一个粗糙范围:Fc 可以用高速段力的中位数近似,dF 可以用低速峰值和高速段之差近似,vs 先给 0.05 m/s 左右,delta 从 2 开始,sigma 从高速段斜率估。边界的作用是挡住负的 vs、过大的 delta 和负的粘性项。
vEps = 1e-3; Fc0 = median(abs(FData(abs(vData) > 0.2))); if isempty(Fc0) || ~isfinite(Fc0) Fc0 = min(abs(FData)); end dF0 = max(abs(FData)) - Fc0; if dF0 <= 0 dF0 = 0.5 * max(abs(FData)); end x0 = [Fc0, dF0, 0.05, 2.0, 0.1]; lb = [0, 0, 1e-4, 1.0, 0]; ub = [max(abs(FData))*2, max(abs(FData))*2, 1.0, 4.0, 10]; opts = optimoptions('lsqcurvefit', ... 'Display', 'iter', ... 'MaxFunctionEvaluations', 2e4, ... 'MaxIterations', 1e4, ... 'StepTolerance', 1e-10); [xFit, resnorm, residual, exitflag, output] = lsqcurvefit( ... @(x,v) stribeckModel(x, v, vEps), ... x0, vData, FData, lb, ub, opts); Fc = xFit(1); dF = xFit(2); Fs = Fc + dF; vs = xFit(3); delta = xFit(4); sigma = xFit(5); rmse = sqrt(mean(residual.^2)); fprintf('Fc=%.4f N, Fs=%.4f N, vs=%.4f m/s, delta=%.3f, sigma=%.4f\n', ... Fc, Fs, vs, delta, sigma); fprintf('RMSE=%.4f N, exitflag=%d\n', rmse, exitflag);这段代码里,vEps必须和最终仿真使用的平滑尺度一致。辨识时用 1e-3,仿真时改成 1e-6,参数看起来拟合得很好,接进 ODE 后零速附近却会抖,这是很常见的“离线拟合”和“在线仿真”不一致问题。另一个细节是ub里的 vs 上限。如果你的速度范围只有 0.1 m/s,vs 上限给 1 m/s 也可以,但给得太大容易让优化器在低灵敏度方向上乱走;一般让上限略大于数据最大速度即可。
3.3 回代验证和残差怎么看
拟合完不要只看 RMSE 一个数。把模型曲线和实测点画在一起,再把残差对速度画出来,很多问题会直接暴露。残差在低速段一边倒,通常是 vs 或 delta 不对;残差在高速段呈线性增长,通常是 sigma 或 Fc 不对;残差在正负速度上符号相反,说明数据本身有零偏或不对称,基础模型需要扩展。
Fhat = stribeckModel(xFit, vData, vEps); res = FData - Fhat; figure; subplot(2,1,1); plot(vData, FData, 'o', 'MarkerSize', 3); hold on; plot(vData, Fhat, 'LineWidth', 1.2); grid on; xlabel('v (m/s)'); ylabel('F (N)'); legend('实测', '拟合'); subplot(2,1,2); plot(vData, res, '.'); grid on; xlabel('v (m/s)'); ylabel('残差 (N)'); r2 = 1 - sum(res.^2) / sum((FData - mean(FData)).^2); fprintf('R^2 = %.4f\n', r2);如果 R² 很高但低速残差仍然有台阶,不要急着加高阶项。先检查每个速度点是不是真的达到了稳态,再检查力传感器零偏和夹具柔性。Stribeck 参数辨识最怕拿动态数据当稳态数据用,尤其是低速爬行阶段,速度本身就在波动,拟合出来的 vs 会把爬行频率也吸收进去。
4. 接进动力学仿真:ode45/ode15s 里摩擦项的三处硬约束
把 Stribeck 函数接进 ODE 后,问题会从“曲线像不像”变成“数值稳不稳”。摩擦项本质上是速度的函数,放在状态方程里没有错,但符号方向、零速平滑尺度和求解器选择会直接决定仿真能不能跑完。尤其是低速往复运动,摩擦曲线存在负斜率,系统可能自激振荡,这不是代码写错了,而是模型本身允许的现象。
4.1 状态方程里摩擦项放在哪
以最简单的单质量-弹簧-阻尼系统为例,运动方程写成:
m*x'' + c*x' + k*x = u - Ff(x')其中 x 是位移,v=x' 是速度,u 是驱动力。摩擦力方向始终阻碍相对运动,所以方程右边是减去 Ff。写成状态空间时,状态取[x; v],导数就是[v; (u - c*v - k*x - Ff)/m]。这里最容易错的是把 Ff 加到右边,或者在外面又乘了一次 sign,导致摩擦力方向反了,仿真出来的速度会越来越大。
function dx = plantStribeck(t, x, p, m, k, c, uFcn) % x(1): 位移 m % x(2): 速度 m/s v = x(2); Ff = stribeckForce(v, p); u = uFcn(t); dx = [v; ... (u - c*v - k*x(1) - Ff) / m]; end调用时用ode15s比ode45更稳,因为平滑符号在 vEps 较小时会让方程变刚。下面这段是完整调用示例,驱动力的设置要让系统能够越过静摩擦峰值,否则质量块可能停在原地,速度曲线看起来像“没有响应”。
p = struct('Fc',1.8,'Fs',3.2,'vs',0.04, ... 'delta',2.0,'sigma',0.12,'vEps',1e-3); m = 1.0; % kg k = 20; % N/m c = 0.4; % N*s/m uFcn = @(t) 4.0 * (t >= 0.2); % 0.2 s 后施加 4 N 驱动力 odefun = @(t,x) plantStribeck(t, x, p, m, k, c, uFcn); opts = odeset('RelTol', 1e-7, 'AbsTol', 1e-9, 'MaxStep', 1e-3); [t, x] = ode15s(odefun, [0 5], [0; 0], opts); figure; plot(t, x(:,2), 'LineWidth', 1.1); grid on; xlabel('时间 t (s)'); ylabel('速度 v (m/s)'); title('Stribeck 摩擦下的速度响应');MaxStep不要设得太大。低速段的速度变化很小,如果步长跨过了 vs 附近的负斜率区,求解器可能直接错过摩擦力峰值。RelTol和AbsTol也不用盲目追小,先保证模型方向正确,再逐步收紧容差。
4.2 零速附近的正则化与刚性
零速附近是 Stribeck 仿真最容易出问题的地方。tanh(v/vEps)把符号函数抹平后,摩擦力在零速附近变成连续函数,求解器能算,但代价是系统在零速附近存在一个很窄的高梯度区。vEps 越小,梯度越陡,ode45 越容易把步长压得很小;vEps 越大,低速摩擦力越“软”,粘着阶段的保持力会被低估。
我一般按两步调:先让 vEps 等于 vs 的 1/50 左右,跑通仿真;再看零速附近的加速度有没有高频抖振。如果抖振明显,先把 vEps 放大到 vs 的 1/20,确认是不是数值刚性问题;如果放大后抖振消失,但低速保持力明显不对,就不能继续放大,应该考虑事件检测或更完整的预滑动模型。不要用“把 vEps 调到很小”来追求物理真实,数值上付出的代价可能远大于收益。
还有一个常见误解:平滑后的 Stribeck 模型能自动处理静摩擦。它不能。tanh只是把符号函数变连续,零速时摩擦力仍然接近 0,系统在驱动力小于 Fs 时可能会缓慢爬行,而不是完全锁死。要做真正的粘着-滑移切换,需要在 ODE 外面加状态机,检测速度过零和驱动力是否超过阈值,或者换 LuGre 这类带内部状态的模型。
4.3 采样噪声、滤波和求解器顺序
仿真里没有传感器噪声,但如果你要拿仿真和实测对比,就必须考虑实测信号的采样和滤波。常见错误是:实测力做了 50 Hz 低通,仿真输出直接和原始数据比,结果残差里全是相位差;或者辨识时用了滤波后的数据,仿真时却把摩擦函数接在一个未滤波的速度信号上,零速附近被噪声反复触发。
顺序应该是:先确定物理单位和采样率,再对实测速度、力做同样的零相位滤波,然后用滤波后的稳态数据辨识参数,最后在仿真输出上做同样的滤波再比较。采样率至少覆盖你关心的最高动态,低速爬行阶段看起来频率不高,但摩擦力切换的等效频率可能不低。若采样率只有 100 Hz,却去拟合 vs=0.005 m/s 的曲线,低速段的数据点可能根本不够。
fs = 1000; % 采样频率 Hz fc = 50; % 低通截止频率 Hz [b,a] = butter(2, fc/(fs/2)); Ffilt = filtfilt(b, a, FData); % 零相位滤波,避免相位偏移filtfilt需要信号长度足够,数据太短时会报错。遇到这种情况,不要改成filter后继续做相位敏感的比较,先把数据段截长,或者只做趋势检查。
5. 常见问题排查:Stribeck 代码跑不通先看这五处
Stribeck 代码的报错通常不复杂,复杂的是“不报错但结果不对”。下面五类问题我在调试时遇到得最多,基本覆盖了从函数调用到 ODE 仿真的主要翻车点。每一条都按现象、原因、解决来拆,看到类似症状可以先对号入座。
5.1 曲线出现负斜率后的抖动
现象:低速仿真速度曲线出现等幅振荡,或者 ode45 跑到某一段后步长越来越小,最后提示无法满足积分容差。
原因:Fs 大于 Fc 时,Stribeck 曲线在 vs 附近存在负斜率,系统本身可能进入自激振动;如果 vEps 又设得很小,摩擦力在零速附近梯度很陡,数值求解器会把这种负阻尼放大。很多人以为是代码写错了,其实是模型和求解器在一起“打架”。
解决:先把vEps放大到 vs 的 1/20 左右,换ode15s,并限制MaxStep,确认抖动是否来自数值刚性。如果放大后曲线仍然抖,而且抖振频率随 k 或 m 变化,那更可能是真实的 stick-slip 现象,需要事件检测、状态机或 LuGre 模型,不要继续靠调容差硬压。
5.2 参数单位不统一导致拟合飘移
现象:拟合结果看起来不错,但 vs 只有 0.0005,或者 sigma 大得离谱,换一组数据参数就完全变样。
原因:速度列一边用 mm/s,一边用 m/s;力传感器输出是 mV,却直接当成 N 送进lsqcurvefit;时间轴是毫秒,导致速度差分后量纲不对。单位不统一时,优化器仍然能下降,但它找到的是数值上的局部最优,不是物理参数。
解决:在数据入口处统一成 SI 单位,变量名里带单位后缀,例如v_mps、F_N。画图前打印max(abs(vData))和max(abs(FData)),凭数量级就能发现异常。vs 的初值和边界也要跟速度单位一致,不能数据是 mm/s,边界却按 m/s 给。
5.3 初值选错收敛到局部极小
现象:exitflag是正的,RMSE 也不大,但拟合曲线只贴住了高速段,低速峰值完全没出来;或者 Fs 贴着下界,delta 贴着上界。
原因:Stribeck 模型对 vs 和 delta 比较敏感,初值离真实低速区太远时,优化器会先牺牲低速段去拟合高速直线。再加上如果数据里低速点很少,残差对 vs 的梯度很弱,局部极小很常见。
解决:用参数化dF = Fs - Fc保证 Fs 不小于 Fc;初值从数据估,不要全填 1;把 vs 的边界限制在数据速度范围内;实在不稳就做几组不同初值的lsqcurvefit,取 RMSE 最小且参数不贴边的那组。贴边参数不是不能要,但必须能解释。
5.4 零速判断写成 v==0
现象:代码里写了if v == 0的分支,但仿真中这个分支几乎从不执行;或者零速附近摩擦力突然变成 0,速度出现一个小台阶。
原因:浮点数很少精确等于 0,ode15s 给出的速度可能是 1e-12、-3e-15 这种量级。用==判断零速,等于把控制逻辑交给舍入误差。
解决:改成abs(v) < vTol,vTol跟 vEps 或采样分辨率挂钩。如果要处理粘着状态,不要只靠一个阈值判断,应该结合速度符号变化、驱动力是否超过 Fs 和停留时间做状态机。阈值本身要写在参数里,不要在函数内部写魔数。
5.5 向量化后维度对不上
现象:单独传标量没问题,一传linspace生成的速度向量就报Matrix dimensions must agree,或者 plot 时 x 和 y 长度不一致。
原因:函数里用了*、/、^而不是.*、./、.^;或者输入 v 是行向量,内部又用zeros(n,1)生成了列向量,拼在一起维度错位。Matlab 不会自动帮你把“元素对元素”的意图翻译成点运算。
解决:函数入口统一v = v(:);所有跟 v 有关的乘除幂都用点运算;输出前用assert(isequal(size(Ff), size(v)))做一次自检。示例脚本里故意传一个行向量和一个列向量各跑一遍,能提前发现大部分维度问题。
6. 进阶用法:把 Stribeck 模型做成参数扫描和验收脚本
拟合出一条曲线只是起点,真正要交付的是“这组参数在什么范围内可信”。我习惯在拟合之后做一次 vs 和 delta 的二维扫描,把 RMSE 画成热图,看最优点是不是落在一条很扁的谷底上。如果谷底很宽,说明数据对这两个参数约束不足,换一组工况参数就会漂;如果最优点紧贴边界,说明模型结构或数据范围有问题。
vGrid = linspace(0.005, 0.30, 60); dGrid = linspace(1.0, 3.0, 41); [D, V] = ndgrid(dGrid, vGrid); % 行对应 delta,列对应 vs rmseMap = zeros(size(V)); for ii = 1:numel(V) xTmp = xFit; xTmp(3) = V(ii); xTmp(4) = D(ii); Ftmp = stribeckModel(xTmp, vData, vEps); rmseMap(ii) = sqrt(mean((Ftmp - FData).^2)); end figure; imagesc(vGrid, dGrid, rmseMap); axis xy; colorbar; xlabel('v_s (m/s)'); ylabel('\delta'); hold on; plot(xFit(3), xFit(4), 'r+', 'MarkerSize', 12, 'LineWidth', 1.5); [~, idx] = min(rmseMap(:)); [i1, i2] = ind2sub(size(rmseMap), idx); fprintf('网格最优: vs=%.4f, delta=%.3f, RMSE=%.4f\n', ... vGrid(i2), dGrid(i1), rmseMap(idx));扫参时不要只盯着最小 RMSE。把参数代回后检查三件事:低速段残差有没有系统偏置;Fs 是否大于 Fc;参数有没有贴着你设定的上下界。验收表可以简单写成下面这样,够用就行。
| 验收项 | 通过标准 | 不通过先看 |
|---|---|---|
| RMSE | 小于传感器噪声的 2 倍左右 | vs、delta 和稳态数据质量 |
| 低速残差 | 正负速度上无明显单边偏置 | 零偏、夹具柔性、vEps |
| 参数边界 | 最优参数不贴上下界 | 初值、数据速度覆盖范围 |
| 回代仿真 | 速度响应无异常高频抖振 | 求解器、MaxStep、平滑尺度 |
从那次低速爬行仿真连续抖了三天之后,我每次把 Stribeck 参数交给别人之前,都会强制走一遍“曲线形状检查、残差对速度检查、ode15s 回代、参数边界检查”这四步,不再只看一条拟合曲线顺不顺眼。希望帮到你。
本文还有配套的精品资源,点击获取