PharmaPy 干法制粒二次开发教程（14）：PSD 报告层——从离散 bin 算 d10/d50/d90、span 与细粉占比，lognormal/GGS/RRB 拟合

PharmaPy 干法制粒二次开发教程14PSD 报告层——从离散 bin 算 d10/d50/d90、span 与细粉占比lognormal/GGS/RRB 拟合版本声明块工具/软件本项目二次开发的 PSD 报告函数本系列二次开发示例基于 lognormal/GGS/RRB 标准分布与分位数定义语言/环境Python 3.9 numpy/scipyscipy.stats.lognorm、scipy.optimize.curve_fit本文目标把篇 13 跑出的离散 PSD 提炼成工程师真的会写进报告的那几个数字并学会用三种分布拟合还原连续曲线。一句话结论用PSD 质量份额的累积插值在跨过 10% / 50% / 90% 的地方算出d10/d50/d90再得span(d90−d10)/d50把离散分布在累积坐标上用scipy.optimize.curve_fit拟合lognormald50exp(μ)、σln GSD、GGSY100(x/k)α**、**RRBF(d)1−exp(−(d/λ)k)、d50λ(ln2)^(1/k)三种连续曲线用 R² 选更贴者——这套函数可直接对接篇 13 的输出与篇 16 的细粉 CQA。〇、本篇要解决的认知问题从一堆离散 bind10/d50/d90 到底怎么插出来而不是拍脑袋span 为什么要用(d90−d10)/d50而不是诸如 d90−d10它告诉工程师什么细粉占比本项目 125µm质量份额在离散网格上怎么精确算含 bin 边界lognormal / GGS / RRB 三种拟合参数各是什么、各自什么物理/工程语义什么时候用哪个拟合怎么打分curve_fit R² 会不会过拟合水分一、机制解析为什么要报告层而不是直接读 bin离散 PBM 的终点是N_f每 bin 数量密度一长串数字。工程师不关心bin 17 的 N...关心的是D50 是多少、分布宽不宽、细粉多不多。报告层的职责就是把高维的 PSD 压成几个互不冗余、可跨批次对比的标量。它的难点不在算而在一致性与可对比性——同一数据用不同网格、不同口径会得到略有差异的 d 值所以必须固定算法与口径比如质量份额累积 线性插值 10/50/90才能让篇 16 的工艺窗口、篇 17 的标定跨工况可比。分位数怎么从累积分布取给一份每 bin 质量份额w_i就能得到累积质量通过率Q(x)Σ_{x_i≤x} w_i。d10/d50/d90 就是使Q10/50/90%的粒径。离散网格上两格之间用线性插值在 log–liner 粒径轴上先在 bin 边界上把Q累出来找到Q跨过目标百分位的两相邻 bin再线性找粒径。这与动态图像分析 Q3 分布分位数的定义一致底账 F。span 与细粉占比两个互补的展宽/细端读数span (d90−d10)/d50无量纲的相对展宽扣除了d50的绝对尺度影响——所以跨不同批次、不同规模even 不同物料可比。它回答分布收得紧不紧。所以是比值而非差若用d90−d10大粒径批次天然更宽就没有可比性了。细粉占比本系列固定125µm的质量份额底账 F。网格上125未必落在 bin 边界上需在Q(x)上插出恰好到 125µm 已通过多少质量那才是125的质量份额。这比数 bin 个数精确得多bin 往往上百 µm 宽、跨 125 界。三种拟合log-normal / GGS / RRB同一份离散 PSD 可以用多种连续模型刻画各有侧重底账 C4/G模型形式关键参数工程语义适用lognormalln d ~ N(μ,σ)d50exp(μ)、几何标准差GSDexp(σ)中位径 几何展宽对称对数轴大多数粉碎/整粒产物最常用GGSY100(x/k)^α粒度模数k≈最大粒径、陡度α对数轴直线用α描述斜率贫粒度范围、早期累积SchuhmannRRB/WeibullF(d)1−exp(−(d/λ)^k)特征粒径λ、模数k中位径d50λ(ln2)^(1/k)宽分布、尾端约束好拟合用curve_fit对累积通过率拟合residual 取Q_model−Q_dataR² 打分。注意过拟合红线R² 高不代表物理对还看残差形状是否系统偏与参数合理性如 GGS 的k是否≈观察到的最大粒径。三模型各给一行结论报告里写用 RMSE 最小者即可不强求唯一。二、完整代码与逐行剖析下面一份可运行报告层从离散x_gr/distrib质量份额算 d10/d50/d90、span、细粉占比再拟合成 lognormal/GGS/RRB 并打分。importnumpyasnpfromscipy.optimizeimportcurve_fitfromscipy.statsimportlognormfromscipy.integrateimportcumulative_trapezoid# ---- 输入一份离散 PSD质量份额来自篇13输出----x_lo,x_hi,nb10.0,2000.0,40ratio(x_hi/x_lo)**(1.0/(nb-1))x_grx_lo*ratio**np.arange(nb)x_cnp.sqrt(x_gr[:-1]*x_gr[1:])volx_c**3.0wnp.where(x_c600,1.0,0.0)**3np.exp(-((np.log(x_c)-np.log(300.0))**2)/(2*0.6**2))ww/w.sum()# 归一质量份额# ---- 1) 累积质量通过率 Q(x) ----Qnp.concatenate(([0.0],np.cumsum(w)))# 在节点上累到 0..1x_nodenp.concatenate(([x_gr[0]],x_c))# 与 Q 对齐的粒径轴defd_pct(p):线性(对数轴)插值回 Qp 处的粒径。pp/100.0ifp0:returnx_node[0]ifp1:returnx_node[-1]inp.searchsorted(Q,p)-1imax(0,min(len(Q)-2,i))# log 轴上线性插值粒径取 log 求 dz(p-Q[i])/(Q[i1]-Q[i])returnnp.exp(np.log(x_node[i])z*(np.log(x_node[i1])-np.log(x_node[i])))d10,d50,d90d_pct(10),d_pct(50),d_pct(90)span(d90-d10)/d50 f125float(d_pct(125)/1.0)*1.0# Q(125)125µm质量份额print(fd10{d10:.1f}d50{d50:.1f}d90{d90:.1f}µm)print(fspan{(d90-d10)/d50:.3f}细粉(125µm)份额{f125:.3f})# ---- 2) 三种拟合对累积通过率 Q ----deflogn_cdf(x__,d50,gsd):returnlognorm.cdf(x__,snp.log(gsd),scaled50)defggs(x__,k,al):returnnp.minimum(1.0,(x__/k)**al)defrrb(x__,lam,kk):return1.0-np.exp(-((x__/lam)**kk))x_fitx_c Q_fitnp.cumsum(w)/w.sum()p_logncurve_fit(lambdax,p,q:logn_cdf(x,p,q),x_fit,Q_fit,p0[300.0,2.0],maxfev20000)p_ggscurve_fit(ggs,x_fit,Q_fit,p0[2000.0,0.8],maxfev20000)p_rrbcurve_fit(rrb,x_fit,Q_fit,p0[300.0,1.2],maxfev20000)defr2(y,yhat):sstnp.sum((y-y.mean())**2);ssrnp.sum((y-yhat)**2);return1-ssr/sstprint(flognormal: d50{p_logn[0][0]:.1f}, GSD{p_logn[0][1]:.3f}, R²{r2(Q_fit,logn_cdf(x_fit,*p_logn[0])):.4f})print(fGGS : k {p_ggs[0][0]:.1f}, alpha{p_ggs[0][1]:.3f}, R²{r2(Q_fit,ggs(x_fit,*p_ggs[0])):.4f})print(fRRB : lam{p_rrb[0][0]:.1f}, k {p_rrb[0][1]:.3f}, R²{r2(Q_fit,rrb(x_fit,*p_rrb[0])):.4f})逐行要点d_pct在对数粒径轴上线性插值Qp处的粒径z是相邻节点间的线性插值权重。这是分位数累积插值的标准实现p50即 d50。span(d90-d10)/d50是无量纲相对展宽底账 F这是它能跨批次/规模可比的原因。f125直接用d_pct(125)得到恰好 125µm 已通过的质量份额——因为d_pct本质是Q的反函数Q(125)正是125质量份额网格 125 边界精确插值胜过数 bin。拟合把curve_fit的目标设在累积通过率Q_fitlognormal 用lognorm.cdf(s,lscale)sln GSD、scaled50GGS 用(x/k)^α截断 1RRB 用指数形式。p0给合理初值辅助收敛。R²用ssr/sst打分这里的更好只作选型参考还要看残差形状与参数合理性见报错。三、常见报错与排查d_pct 在 90% 处插值失败 / 返回末端粒径。现象d90恒定等于x_gr[-1]。根因Q在尾部没到 1cumsum后最后一格权重小或网格太粗np.searchsorted(Q,0.9)落在最后且z插不到位。解法确认w归一w.sum()必要时把最右 bin 也纳入插值用min(1, ...)夹住 p。拟合报OptimizeWarning: Covariance of parameters...。现象curve_fit警告但给出结果参数可能是inf/噪声。根因初值p0差、模型对某参数不敏感如 GGS 的k在尾部数据少时无法识别。解法换更合理的p0d50直觉值、α≈0.5~1.5并增加maxfev或对 GGS 用对数‑对数在直线段先手算k/α当初值。lognormal 用lognorm但gsd与直觉不符。现象拟出GSD偏大/偏小d50与d_pct(50)差多。根因lognorm.cdf的形参辨析——s是σln GSDscaled50若误把s当GSD填2.0分布会被压窄/拉宽。解法参照定义ln d~N(μ,σ)d50exp(μ)、GSDexp(σ)据此正确传s。细粉占比 num bin 计数与插值法不符。现象f125与统计x_c125的 bin 质量差很多。根因bin 中心制 vs 边界制——125常落在某 bin 中心附近按 bin 中心数会把该档整档归入或不归。解法坚持累积插值Q(125)本篇口径并把口径写进注释避免与数 bin混用报告要一致可比。三种拟合 R² 都高却拿不准选哪个。现象R²都 0.98 但结论不同。根因R²只看整体吻合掩盖残差系统性。解法看残差符号分布与参数合理性——若 lognormal 中段系统偏 、GGS 尾端超调选残差最随机的那个工程上多数产物用 lognormal 就够GGS/RRB 留作敏感位置尾端/细端的补救。四、动手练习对账基准造一个标准对数正态x~lognorm(scale300,s0.5)把它离散到网格再跑d10/d50/d90。判定d50≈300、d90≈300*exp(1.2816*0.5)exp-上浮与解析值偏差2%证明插值分位数算法正确。口径一致性对篇 13 的同一份d_fine合格品 PSD用插值Q(125)与bin 中心计数两种口径都算细粉占比。判定两者大幅不一致时你能否解释为什么、并把口径统一——训练跨批次可比性。三模型选型分别对窄分布σ0.3宽分布σ1.0“拟三种模型。判定窄分布 lognormal 最贴R² 最高、宽分布 RRB 常胜出——你能用数据说出什么时候用哪种”为篇 16-20 报告定型。五、小结与下一篇预告本篇补齐了从离散 PSD 到报告数字的最后一层d_pct对数轴累积插值算 d10/d50/d90span无量纲展宽Q(125)精确细粉占比并用curve_fit把与离散数据拟合成 lognormal/GGS/RRB 三种连续曲线、R² 选型。你手里现在有了统一的 PSD 指标体系——它会让篇 16 的工艺窗口、篇 17 的标定、篇 20 的报告全都站在同一把尺子上。下一篇篇 15《颗粒强度模块》——离开粒径这一维度加强度从 ribbon 抗张强度/孔隙率角度预测颗粒破坏载荷理解预压缩硬化与粒度增大 → 结合面积↓为何让干法颗粒压片性损失。本篇认知问题回显FAQQ1离散 bin 怎么把 d10/d50/d90 插出来A1先按每档质量份额累出累积通过率Q(x)d10/d50/d90 就是Q10/50/90%处的粒径。离散网格上在对数粒径轴做线性插值找到Q跨过目标百分位的两相邻 bind_pct(p)用searchsortedz插值回粒径。这与动态图像分析 Q3 分位数口径一致。Q2span 为什么用(d90−d10)/d50A2因为 span 要以无量纲表示相对展宽扣掉d50的绝对尺度才能跨不同批次、规模、物料可比。若用d90−d10绝对差大粒径批次天然更宽失去可比性。它告诉工程师分布收得紧不紧。Q3细粉占比125µm离散网格上怎么精确算A3不数 binbin 常上百 µm 宽、125 往往落在 bin 中间而是用累积插值Q(125)——直接在累积通过率曲线上插出恰好 125µm 已通过多少质量那就是125的质量份额。口径固定成累积插值才能在篇 16 跨工况可比。Q4log-normal / GGS / RRB 参数与语义A4lognormalln d~N(μ,σ)d50exp(μ)、几何标准差GSDexp(σ)对称对数轴多数产物最常用GGSY100(x/k)^αk≈最大粒径、α陡度对数轴直线描述斜率适合贫范围/早期累积RRBF1−exp(−(d/λ)^k)特征粒径λ、模数k、d50λ(ln2)^(1/k)适合宽分布/尾端约束好。Q5拟合怎么打分、会不会过拟合数轴多数产物最常用GGSY100(x/k)^αk≈最大粒径、α陡度对数轴直线描述斜率适合贫范围/早期累积RRBF1−exp(−(d/λ)^k)特征粒径λ、模数k、d50λ(ln2)^(1/k)适合宽分布/尾端约束好。Q5拟合怎么打分、会不会过拟合A5用curve_fit对累积通过率拟合残差取Q_model−Q_dataR²1−ssr/sst打分。但 R² 高不是全部——还要看残差形状是否系统偏与参数合理性如 GGS 的k应≈观察到的最大粒径。工程上多数产物 lognormal 够用GGS/RRB 用作尾端/细端敏感位置的补救。