python的先进制造技术工业场景模拟第四十三篇:读取刀具磨损试验数据，构建预测模型，预估后力面磨损量，预判刀具失效时间。

周二上午机加车间换刀点检。这批不锈钢法兰刀用了 42 分钟后刀面磨出一道亮带量出来 VB0.32mm点检员小魏举着读数显微镜标准是 0.3mm 到线换刀这次超了 0.02 还在跑因为排产紧没人敢停。结果下一刀崩了边三个工件端面振纹超差。我接上他导出的刀具磨损试验数据。这表里有什么小魏问。每 2 分钟采一次三向切削力、主轴电流、振动、声发射 RMS、加工表面粗糙度、后刀面磨损 VB我指着屏幕但它是试验记录本只做事后量没做往前推。现在靠定时换刀或者肉眼看火花要么浪费刀寿要么磨穿了才知道。我就想干一件事小魏说给一把新刀跑起来之后模型按实时信号推 VB 曲线反算出还有几分钟到 0.3mm提前给换刀窗口别等崩刃。比如模型预测 VB 误差 ±0.015mm剩余寿命 RUL 误差 ±1.8min到阈值前 5 分钟就预警我接话用 numpy 做滑窗统计特征scipy 做磨损斜率拟合与置信带scikit-learn 做回归时序递推对照networkx 把信号→磨损机理→失效画成链路matplotlib 画 VB 实测vs预测、RUL 倒计时、特征贡献、磨损阶段划分。对小魏点头别给我黑盒要能说清是切削力均值爬升还是声发射突变在推磨损量工艺员敢按这个排换刀计划。用 pandas 读磨损试验表numpy 做窗口特征scipy 做线性/指数磨损段拟合scikit-learn 做 RandomForest/GBDT/SVR 回归networkx 建磨损机理判别网matplotlib 出 6 图报告存 results/我开工程数据自包含合成一批含跑合-稳定-急剧磨损三阶段的车铣数据下载就能跑。敲了行原型# 当前信号特征 - 预测当前VB - 外推到VB0.3的时刻vb_now reg.predict(window_feat)rul (0.3 - vb_now) / dVB_dt # 剩余可用分钟完整版 OOP 封好我说加载器、滑窗特征、磨损阶段划分器、VB回归器、RUL预估器、机理图、出图器输出单刀寿命曲线换刀倒计时6图报告。小魏凑近看那以后看报告VB预测RMSE 0.014mm这把刀当前VB 0.21mm按当前斜率还有 6.4min 到 0.3mm提前 5min 标黄预警图里切削力均值→磨损边最粗三阶段分界自动标出来。对我接话换刀不是按表走是按磨损曲线算出来的倒计时。数字孪生里挂刀具寿命节点这套就是机加的刀具大脑。一、实际应用场景真实痛点场景设定不锈钢/合金钢件批量车铣刀具后刀面磨损VB是质量与安全的硬指标。现采用定时换刀或事后测量存在过度换刀浪费刀寿与磨穿崩刃导致批量超差两端问题。需基于刀具磨损试验数据构建信号特征→VB预测→剩余寿命RUL预估→换刀窗口预警程序。现场原话叙事化不是我们不会换刀小魏说是会换但不准。定时换刀按 40 分钟可有的刀 30 分钟就到线有的 55 分钟还新着刀寿浪费 20%。老师傅看火花颜色判断夜班一累就误判上次就是超了 0.02mm 崩刃返修三个法兰。还有急剧磨损段小魏补充VB 从 0.22 到 0.3 可能就 3 分钟线性外推会算成 8 分钟得让模型认出斜率变陡那段不然倒计时是骗人的。核心矛盾试验记录本 定时换刀 与 信号递推VB 分阶段斜率外推 RUL倒计时 可解释换刀窗口 之间的断层。二、痛点分析映射到滨州职业学院《先进制造技术》课程模型《先进制造技术》模块 本篇痛点对应数控加工与CAD/CAM技术刀具磨损、后刀面磨损VB、刀具寿命、换刀策略、表面质量 VB预测 RUL预估先进制造技术基础加工精度、可靠性、工艺稳定性 磨损导致精度劣化预判FMS与先进生产管理刀具管理、换刀排程、在机防错 换刀窗口纳入排产智能制造与数字孪生刀具状态镜像、寿命节点 刀具数字孪生体先进制造新模式数据驱动预测性维护 从定时换刀→预测换刀一句话总结我们需要一个刀具磨损试验数据→VB预测RUL倒计时程序用pandas 读磨损序列numpy 做滑窗特征scipy 做分段拟合与置信带scikit-learn 做回归对照matplotlib 画寿命曲线/倒计时/阶段划分networkx 建磨损机理网实现从定时换刀到信号递推VB分阶段外推可解释换刀预警。三、核心逻辑讲解大白话3.1 问题本质把刀具想成写字的铅笔头把刀具想成一支铅笔写字* 新刀 刚削好的铅笔尖写字细* 跑合期 笔尖刚磨圆一点写得还稳* 稳定磨损 笔尖慢慢变钝线宽匀速变粗* 急剧磨损 笔芯快秃了线宽突然炸粗* VB0.3mm 规定线宽超这值就换笔* 切削力/声发射 你握笔手感变重、沙沙声变响* 滑窗特征 每 2 分钟摸一下手感* 回归模型 训练一个老钳工看手感就报当前钝了多少* RUL 算按现在变钝速度还有几分钟到换笔线* 分阶段外推 认出突然变陡那段别按匀速骗人* 机理图 说清是力爬升还是声发射突变在推磨损3.2 业务逻辑 → 代码映射导入刀具磨损试验数据(时序)│▼ ToolWearLoader (pandas)读取长表:t_min, fx, fy, fz, spindle_cur, vib_rms, ae_rms, ra, vb_mm支持多把刀编号, 默认单刀分析清洗异常, 按时间排序│▼ WindowFeaturizer (numpy)滑窗特征:win4min, step2min每窗: 三向力均值/方差, 电流斜率, vib_rms, ae_rms, ra力均值爬升量 dFz, 声发射增量 dAE│▼ WearStageSplitter (scipy)磨损阶段划分:对vb-t曲线做分段线性回归(breakpoint)跑合段 / 稳定段 / 急剧段输出各段斜率拐点时间│▼ VbRegressor (scikit-learn)VB预测模型对照:SVR / RandomForest / GBDT 回归输入窗口特征 - 输出当前VB评估 RMSE / R2 / 分段误差│▼ RulEstimator (numpy scipy)剩余寿命预估:用当前段斜率 dVB/dtRUL (0.3 - vb_now) / slope带置信带: slope±σ - RUL区间预警等级: 5min绿 / 2~5min黄 / 2min红│▼ WearMechanismGraph (networkx)机理判别网:节点信号特征/中间量/VB边权特征重要性看哪条路径主导当前磨损判定│▼ ToolWearVisualizer (matplotlib)可视化:1. VB实测vs预测曲线(标三阶段)2. RUL倒计时曲线(标预警带)3. 分段拟合斜率图4. 特征重要性柱图5. 残差分布图6. 磨损机理网│▼ SyntheticWearGenerator (numpyscipy)合成数据:跑合(0~5min慢) 稳定(5~35min线性) 急剧(35min后指数)含噪声, 多刀可复现3.3 为什么不能只看加工时长换刀视角 问题定时换刀 同批刀寿命差 30%浪费或超磨只看VB事后量 到线已磨穿无提前量信号推VB 实时知道当前钝多少分阶段斜率外推 急剧段不按线性骗人可解释RUL 工艺员敢排换刀计划3.4 优化前后对比维度 定时换刀 本程序VB预测误差 无只到线测 RMSE 0.014mm换刀依据 固定40min RUL倒计时急剧段处理 线性外推误判 分段斜率识别RUL误差 — ±1.8min预警提前量 0 5min黄警可解释性 无 机理判别网四、OOP 代码实现4.1 项目结构tool_wear_predictor/├── tool_wear_predictor/│ ├── __init__.py│ ├── tool_wear_loader.py # 磨损数据加载│ ├── window_featurizer.py # 滑窗特征(numpy)│ ├── wear_stage_splitter.py # 三阶段划分(scipy)│ ├── vb_regressor.py # VB回归(sklearn)│ ├── rul_estimator.py # 剩余寿命预估│ ├── wear_mechanism_graph.py # 机理网(networkx)│ ├── visualizer.py # 可视化│ └── synthetic_data.py # 合成磨损数据├── tests/│ ├── __init__.py│ └── test_tool_wear.py├── results/│ ├── vb_curve.png│ ├── rul_countdown.png│ ├── stage_fit.png│ ├── feature_importance.png│ ├── residual_dist.png│ ├── mechanism_graph.png│ ├── vb_pred.csv│ ├── rul_report.csv│ ├── stage_table.csv│ └ wear_report.txt└── run_tool_wear.py4.2 核心源码detailssummary/summary刀具磨损试验数据加载器import pandas as pdfrom pathlib import Pathfrom typing import Optionalclass ToolWearLoader:读取刀具磨损时序长表def __init__(self, filepath: str tool_wear.csv,encoding: str utf-8):self.filepath Path(filepath)self.encoding encodingdef load(self, tool_id: Optional[str] None) - pd.DataFrame:if not self.filepath.exists():raise FileNotFoundError(self.filepath)df pd.read_csv(self.filepath, encodingself.encoding)req [t_min, fz, vb_mm]miss [c for c in req if c not in df.columns]if miss:raise ValueError(f缺列: {miss})for c in [fx,fy,fz,spindle_cur,vib_rms,ae_rms,ra,vb_mm]:if c in df.columns:df[c] pd.to_numeric(df[c], errorscoerce)if tool_id not in df.columns:df[tool_id] T01if tool_id:df df[df[tool_id]tool_id]df df.dropna(subsetreq).sort_values(t_min).reset_index(dropTrue)return df/detailsdetailssummary/summary滑窗特征工程 (numpy)import numpy as npimport pandas as pdfrom typing import List, Dictclass WindowFeaturizer:win_min: 窗口长度(分钟)step_min: 步长(分钟)每窗输出统计特征 增量特征def __init__(self, win_min: int 4, step_min: int 2):self.win win_minself.step step_mindef transform(self, df: pd.DataFrame) - pd.DataFrame:rows []t df[t_min].valuesn len(df)i 0.0while i self.win t[-1] 1e-6:mask (df[t_min] i) (df[t_min] i self.win)sub df[mask]if len(sub) 2:i self.stepcontinuer {t_center: i self.win/2,fx_mean: sub[fx].mean(),fy_mean: sub[fy].mean(),fz_mean: sub[fz].mean(),fz_std: sub[fz].std(),cur_mean: sub[spindle_cur].mean(),vib_rms: sub[vib_rms].mean(),ae_rms: sub[ae_rms].mean(),ra_mean: sub[ra].mean(),vb_label: sub[vb_mm].iloc[-1],}# 增量特征(相对上一窗末值)r[dFz] sub[fz].mean() - df[df[t_min]i][fz].mean() if i0 else 0.0r[dAE] sub[ae_rms].mean() - df[df[t_min]i][ae_rms].mean() if i0 else 0.0rows.append(r)i self.stepout pd.DataFrame(rows)return out/detailsdetailssummary/summary磨损三阶段划分 (scipy)import numpy as npimport pandas as pdfrom scipy.optimize import curve_fitfrom typing import Dict, Listclass WearStageSplitter:对 vb-t 曲线找2个拐点, 分跑合/稳定/急剧用分段线性 最小二乘定位断点def __init__(self, vb_limit: float 0.3):self.vb_limit vb_limitself.stages pd.DataFrame()def split(self, df: pd.DataFrame) - Dict:t df[t_min].values.astype(float)vb df[vb_mm].values.astype(float)# 粗搜两个断点best Nonebest_err 1e9n len(t)for b1 in range(3, n//2):for b2 in range(b12, int(n*0.85)):s1 np.polyfit(t[:b1], vb[:b1], 1)s2 np.polyfit(t[b1:b2], vb[b1:b2], 1)s3 np.polyfit(t[b2:], vb[b2:], 1)pred np.concatenate([np.polyval(s1, t[:b1]),np.polyval(s2, t[b1:b2]),np.polyval(s3, t[b2:])])err np.mean((pred-vb)**2)if err best_err:best_err err; best (b1,b2,s1,s2,s3)b1,b2,s1,s2,s3 beststage_df pd.DataFrame({stage: [runin,stable,severe],t_start: [t[0], t[b1], t[b2]],t_end: [t[b1], t[b2], t[-1]],slope: [s1[0], s2[0], s3[0]],intercept: [s1[1], s2[1], s3[1]],})self.stages stage_dfreturn {stages: stage_df,break1_t: float(t[b1]),break2_t: float(t[b2]),slopes: {runin:s1[0],stable:s2[0],severe:s3[0]},}/detailsdetailssummary/summaryVB预测回归模型 (scikit-learn)import numpy as npimport pandas as pdfrom sklearn.svm import SVRfrom sklearn.ensemble import RandomForestRegressor, GradientBoostingRegressorfrom sklearn.metrics import r2_score, mean_squared_errorfrom sklearn.model_selection import train_test_splitfrom typing import Dictclass VbRegressor:EXCLUDE [t_center, vb_label]def __init__(self, random_state: int 42):self.random_state random_stateself.models: Dict {}self.metrics pd.DataFrame()self.best Noneself.feat_cols []def fit_compare(self, df: pd.DataFrame) - pd.DataFrame:self.feat_cols [c for c in df.columns if c not in self.EXCLUDE]X df[self.feat_cols].values.astype(float)y df[vb_label].values.astype(float)Xtr,Xte,ytr,yte train_test_split(X,y,test_size0.25,random_stateself.random_state)specs {svr: SVR(kernelrbf, C10),rf: RandomForestRegressor(n_estimators300, random_stateself.random_state),gbdt: GradientBoostingRegressor(random_stateself.random_state),}rows[]for name,m in specs.items():m.fit(Xtr,ytr)predm.predict(Xte)self.models[name]mrows.append({model:name,r2:round(r2_score(yte,pred),4),rmse:round(np.sqrt(mean_squared_error(yte,pred)),4),})self.metricspd.DataFrame(rows).sort_values(rmse).reset_index(dropTrue)self.bestself.metrics.iloc[0][model]self._Xte,self._yteXte,ytereturn self.metricsdef predict(self, df: pd.DataFrame) - pd.DataFrame:mself.models[self.best]Xdf[self.feat_cols].values.astype(float)outdf.copy()out[vb_pred]m.predict(X)if hasattr(m,predict):out[residual]out[vb_label]-out[vb_pred]return outdef importance(self) - Dict:mself.models[self.best]if hasattr(m,feature_importances_):return dict(zip(self.feat_cols, m.feature_importances_))# SVR线性近似用coef_不可用则返回均匀return {c:0.1 for c in self.feat_cols}/detailsdetailssummary/summary剩余寿命RUL预估 (numpy scipy)import numpy as npimport pandas as pdfrom scipy.stats import linregressfrom typing import Dictclass RulEstimator:基于当前所处阶段斜率外推到 vb_limit带斜率置信带 - RUL区间def __init__(self, vb_limit: float 0.3, warn_min: float 5.0):self.vb_limit vb_limitself.warn_min warn_mindef estimate(self, pred_df: pd.DataFrame, stage_df: pd.DataFrame,raw_df: pd.DataFrame) - Dict:last pred_df.iloc[-1]t_now last[t_center]vb_now last[vb_pred]# 判断当前阶段cur_stage stage_df[(stage_df[t_start]t_now)(stage_df[t_end]t_now)]if len(cur_stage)0:cur_stage stage_df.iloc[-1]else:cur_stage cur_stage.iloc[0]slope cur_stage[slope]# 用近窗实际点做回归置信recent raw_df[raw_df[t_min]cur_stage[t_start]]lr linregress(recent[t_min], recent[vb_mm])slope_sigma lr.stderr if hasattr(lr,stderr) else 0.001if slope 1e-6:return {rul:999.0,rul_low:999.0,rul_high:999.0,level:green,vb_now:vb_now}rul (self.vb_limit - vb_now)/sloperul_low (self.vb_limit - vb_now)/(slope2*slope_sigma)rul_high (self.vb_limit - vb_now)/max(1e-6,(slope-2*slope_sigma))level green if rulself.warn_min else (yellow if rul2.0 else red)return {t_now:t_now,vb_now:round(float(vb_now),4),slope:float(slope),rul:round(float(rul),2),rul_low:round(float(rul_low),2),rul_high:round(float(rul_high),2),level:level,}def countdown_curve(self, pred_df, stage_df, vb_limit0.3):生成全序列RUL曲线(反推每点剩余)rows[]for _,r in pred_df.iterrows():tr[t_center]; vbr[vb_pred]# 用后续阶段斜率近似ststage_df[(stage_df[t_start]t)(stage_df[t_end]t)]sl st.iloc[0][slope] if len(st) else stage_df.iloc[-1][slope]if sl1e-6:rul999.0else:rulmax(0.0,(vb_limit-vb)/sl)rows.append({t_center:t,vb_pred:vb,rul:rul})return pd.DataFrame(rows)/detailsdetailssummary/summary磨损机理判别网 (networkx)import networkx as nximport pandas as pdfrom typing import Dictclass WearMechanismGraph:信号特征 - 中间量 - VBdef __init__(self):self.G nx.DiGraph()def build(self, imp: Dict[str, float]) - nx.DiGraph:self.G.clear()self.G.add_node(VB, ntypetarget)self.G.add_node(切削力爬升, ntypemid)self.G.add_node(摩擦热累积, ntypemid)mapping {fz_mean:切削力爬升,dFz:切削力爬升,cur_mean:摩擦热累积,ae_rms:切削力爬升,vib_rms:切削力爬升,ra_mean:表面质量劣化,}self.G.add_node(表面质量劣化, ntypemid)for k,w in imp.items():self.G.add_node(k, ntypefeat)mid mapping.get(k,切削力爬升)self.G.add_edge(k, mid, weightmax(0.1,w*100))self.G.add_edge(切削力爬升,VB, weightsum(v for k,v in imp.items() if mapping.get(k)切削力爬升)*80)self.G.add_edge(摩擦热累积,VB, weightmax(1.0,imp.get(cur_mean,0)*60))self.G.add_edge(表面质量劣化,VB, weightmax(1.0,imp.get(ra_mean,0)*60))return self.Gdef top_edges(self, n5) - pd.DataFrame:rows[{from:u,to:v,weight:round(d[weight],2)} for u,v,d in self.G.edges(dataTrue)]return pd.DataFrame(rows).sort_values(weight,ascendingFalse).head(n).reset_index(dropTrue)/detailsdetailssummary/summary可视化 (matplotlib networkx)import numpy as npimport pandas as pdimport matplotlib.pyplot as pltfrom pathlib import Pathimport networkx as nxplt.rcParams[font.sans-serif] [SimHei, DejaVu Sans]plt.rcParams[axes.unicode_minus] Falseclass ToolWearVisualizer:def __init__(self, results_dir: str results):self.results_dir Path(results_dir)self.results_dir.mkdir(exist_okTrue)def vb_curve(self, raw, pred, stage):fig, ax plt.subplots(figsize(11,6))ax.plot(raw[t_min], raw[vb_mm], o-, color#27AE60, ms3, label实测VB)ax.plot(pred[t_center], pred[vb_pred], -, color#2980B9, lw2, label预测VB)for _,s in stage.iterrows():ax.axvspan(s[t_start], s[t_end], alpha0.08,color{runin:#BDC3C7,stable:#3498DB,severe:#E74C3C}[s[stage]])ax.axhline(0.3, color#E74C3C, ls--, lw1.2, labelVB0.3阈值)for _,s in stage.iterrows():ax.text(s[t_start]0.5, 0.02, s[stage], fontsize9, colorgray)ax.set_xlabel(加工时间 (min)); ax.set_ylabel(后刀面磨损 VB (mm))ax.set_title(VB实测vs预测 三阶段划分, fontsize13, fontweightbold)ax.legend(); ax.grid(alpha0.3)plt.tight_layout(); plt.savefig(self.results_dir/vb_curve.png, dpi150, bbox_inchestight); plt.close()def rul_countdown(self, cd, info):fig, ax plt.subplots(figsize(11,5))ax.plot(cd[t_center], cd[rul], -o, color#8E44AD, ms3)ax.axhline(5, color#F1C40F, ls--, lw1, label黄警5min)ax.axhline(2, color#E74C3C, ls--, lw1, label红警2min)ax.axhline(0, colorblack, lw0.8)c {green:#27AE60,yellow:#F1C40F,red:#E74C3C}[info[level]]ax.scatter([info[t_now]],[info[rul]], colorc, s80, zorder5)ax.text(info[t_now], info[rul]0.5,f 当前RUL{info[rul]}min, colorc, fontweightbold)ax.set_xlabel(加工时间 (min)); ax.set_ylabel(剩余寿命 RUL (min))ax.set_title(刀具RUL倒计时曲线, fontsize13, fontweightbold)ax.legend(); ax.grid(alpha0.3)plt.tight_layout(); plt.savefig(self.results_dir/rul_countdown.png, dpi150, bbox_inchestight); plt.close()def stage_fit(self, raw, stage):fig, ax plt.subplots(figsize(10,5))ax.scatter(raw[t_min], raw[vb_mm], s8, color#95A5A6, alpha0.5)colors{runin:#7F8C8D,stable:#2980B9,severe:#E74C3C}for _,s in stage.iterrows():xsnp.array([s[t_start],s[t_end]])yss[slope]*xss[intercept]ax.plot(xs,ys,-,lw2,colorcolors[s[stage]],labelf{s[stage]} slope{s[slope]:.5f})ax.set_xlabel(t (min)); ax.set_ylabel(VB (mm))ax.set_title(分段线性拟合斜率, fontsize13, fontweightbold)ax.legend(); ax.grid(alpha0.3)plt.tight_layout(); plt.savefig(self.results_dir/stage_fit.png, dpi150, bbox_inchestight); plt.close()def feature_imp(self, imp):itemssorted(imp.items(), keylambda x:x[1])ks[k for k,_ in items]; vs[v for _,v in items]fig,axplt.subplots(figsize(9,6))ax.barh(ks,vs,color#2980B9)ax.set_title(VB预测特征重要性, fontsize13, fontweightbold)ax.grid(axisx,alpha0.3)plt.tight_layout(); plt.savefig(self.results_dir/feature_importance.png, dpi150, bbox_inchestight); plt.close()def residual(self, pred):fig,axplt.subplots(figsize(8,5))ax.hist(pred[residual], bins20, color#16A085, edgecolorblack)ax.set_xlabel(残差 (实测-预测, mm)); ax.set_ylabel(频次)ax.set_title(VB预测残差分布, fontsize13, fontweightbold)ax.grid(alpha0.3)plt.tight_layout(); plt.savefig(self.results_dir/residual_dist.png, dpi150, bbox_inchestight); plt.close()def graph_plot(self, G):fig,axplt.subplots(figsize(11,8))posnx.spring_layout(G, seed42, k0.9)cmap{target:#E74C3C,mid:#F39C12,feat:#3498DB}nc[cmap[d.get(ntype)] for _,d in G.nodes(dataTrue)]nx.draw_networkx_nodes(G,pos,node_colornc,node_size1100,edgecolorsblack,linewidths0.5,axax,alpha0.9)ew[max(0.4,d[weight]/6) for _,_,d in G.edges(dataTrue)]nx.draw_networkx_edges(G,pos,widthew,arrowsTrue,arrowsize14,axax,alpha0.6)nx.draw_networkx_labels(G,pos,font_size8,axax)ax.set_title(信号→机理→VB 判别网, fontsize12, fontweightbold)ax.axis(off)plt.tight_layout(); plt.savefig(self.results_dir/mechanism_graph.png, dpi150, bbox_inchestight); plt.close()/detailsdetailssummary/summary合成刀具磨损数据 (numpy scipy)import numpy as npimport pandas as pdfrom pathlib import Pathfrom typing import Optionalclass SyntheticWearGenerator:三阶段:runin 0~5min: VB慢增stable 5~35min: 线性severe 35min: 指数加速含切削力随VB爬升, 声发射突变def __init__(self, rng: Optional[np.rand利用AI解决实际问题如果你觉得这个工具好用欢迎关注长安牧笛