企业官网建设流程全解析

PharmaPy 干法制粒二次开发教程（17）：参数估计与模型标定实战——用 curve_fit/least_squares 反演 S_i、B_ij、破碎核，AIC/R² 选模

版本声明块

• 工具/软件：scipy.optimize.curve_fit/scipy.optimize.least_squares（SciPy）＋ 本项目二次开发的破碎 PBM（示例）
• 语言/环境：Python 3.9 + numpy/scipy
• 本文目标：把前面所有"示意参数"换成"从实验数据反演出来的真参数"，并学会用统计指标守住"别过拟合"。

一句话结论：把实验测得的全时间 PSD 作为目标，用scipy.optimize.least_squares最小化"模型 PSD − 实测 PSD"的残差，一次反演出选择函数S_i（或其中S_ref/q）与破碎函数B_ij（或其中 GGS 的k/α）与 Vogel–Peukert 参数；用AIC 与 R²在多候选模型中选模（R² 高但 AIC 更苛刻），并把拟合残差并入cov得到参数误差带——这才是把"示意"变成"标定"的最后一跃。

〇、本篇要解决的认知问题

1. "参数标定"和"参数设置"差在哪？为什么非要一个最优化过程而手调不行？
2. curve_fit和least_squares什么场景用哪个？破碎 PBM 这种多输出/离散的怎么设目标？
3. 反演S_i、B_ij时，怎么做到"一次全都估"而不是逐个手试？待估量的自由度怎么控制？
4. AIC 和 R² 到底各自惩罚什么？为什么 R² 高还不够、还要看 AIC 与残差？
5. 拟合出来的参数"准不准"怎么量化？误差带（cov）怎么来的？

一、机制解析

为什么必须"优化"而非"手调"

篇 10-13 用的S_ref/q/k/α/β全是示意常数。但真实参数是从实验来的：你在几个转速/筛网下测了一组"进料 PSD → 出料 PSD"的时间序列（或停留时间分布），要找一个参数集使模型预测最贴近实测。手调几万个组合不现实，而且人调会偏好"眼睛觉得像"而非"统计上最贴"。所以参数标定 = 一个最小化残差的最优化问题：min_θ Σ_i (PSD_model(t_i;θ) − PSD_exp(t_i))²。这就是"标定"相对"设置"的本质差异——有可复现的目标函数与统计指标。

curve_fitvsleast_squares分工（底账 G）

• curve_fit：适合"稀疏参数、标量输出、可写成f(x, *θ)"的标准非线性回归，自动算cov（参数协方差）→ 误差带。
• least_squares：更一般——支持多输出残差向量（可以逐 bin、逐时刻展平）、逐残基加权、边界约束。破碎 PBM 的观测是"整个离散 PSD"而非单个标量，天然是多输出，所以首选least_squares：把(时刻 t_i × bin)的偏差展平成一个大残差向量去最小化。

反演策略：层叠逼近，控制自由度

直接反演每一个B_ij（n×n 个）在数据有限时会严重欠定（参数比数据点还多→过拟合）。所以工程上分两步"层叠"（底账 G，back-calculation + 直接测量）：

1. 先估S_i：低细度批次（几乎无二次破碎）下，S_i主导减质量速率，可用t_n/t_10类或单粒级法先估；再纳入least_squares一起优化。
2. 再估B_ij的紧凑参数：不给每个B_ij一个自由参数，而是用低维参数化（GGS 的k/α、分支比β）约束B_ij的形状——这是避免过拟合的关键。B_ij相对不敏感（底账 C2），几个低维参数就够好。

这样总自由度降到O(几个参数)，least_squares才稳健、可解释、能过 AIC 检验。

AIC 与 R²：一个爱复杂、一个讨厌水晶鞋

• R²=1−SSres/SStot，衡量"模型解释了观测方差的比例"，越高拟合越贴——但它不惩罚参数多少，模型加参数 R² 只会涨不会降。
• AIC=n·ln(RSS/n) + 2k（k=参数个数，底账 G）：显式惩罚参数个数，RSS低但k大 → AIC 升高。所以 AIC 是"在拟合好与参数省之间做权衡"的得分，越低越好。

它俩互补：R² 看"贴不贴"，AIC 看"值不值得"。选模准则：AIC 更低者优先（兼顾简约），R² 高而 AIC 也高说明你在用过多参数堆出来的"虚好"。这正是"示意→标定"里防过拟合的统计那道闸。

误差带（cov）怎么来

curve_fit/least_squares都基于残差在最优解附近的一阶近似给出参数协方差cov≈(JᵀJ)⁻¹·s²（J 是雅可比）。参数标准差就是sqrt(diag(cov))。若残留误差的尺度s²估计合理、模型-数据近似线性，这套误差带就可靠；若残差很大或非线性强，误差带会低估——这也是为什么要把"残差形状"（是否随机）一起检查，别只信误差带数字。

二、完整代码与逐行剖析

下面用一个合成实验演示标定闭环：先用"真参数"跑出合成实测 PSD，再用least_squares反演回参数、验 AIC/R²、给出误差带。合成数据让"真值已知"，才能验收反演是否真的还原。

importnumpyasnpfromscipy.optimizeimportleast_squaresfromscipy.integrateimportsolve_ivp# ---- 网格 ----x_lo,x_hi,nb=10.0,2000.0,40ratio=(x_hi/x_lo)**(1.0/(nb-1));x_gr=x_lo*ratio**np.arange(nb)x_c=np.sqrt(x_gr[:-1]*x_gr[1:]);vol=x_c**3.0# ---- PBM 求解（复用篇10）----defsolve_pbm(S,B,N0,t_end=30.0):defrhs(t,N):returnnp.dot(B,S*N)-S*Nreturnsolve_ivp(rhs,(0,t_end),N0,method='BDF',rtol=1e-8,atol=1e-10)# ---- "真参数"造合成数据 ----k_ref,q_ref=0.8,0.9S_true=k_ref*(x_c/500.0)**q_ref B=np.zeros((nb,nb))forjinrange(nb):i=max(1,min(nb-1,int(np.searchsorted(vol,0.5*vol[j]))))w=(vol[i]-0.5*vol[j])/(vol[i]-vol[i-1])B[i-1,j]+=w;B[i,j]+=(1-w)N0=np.where(x_c<=500,1.0,0.0)sol=solve_pbm(S_true,B,N0)t_exp=np.linspace(10,30,6)# 观测时刻PSD_exp=np.array([(solve_pbm(S_true,B,N0,te).y[:,-1])forteint_exp])# (ntime x nb)# ---- least_squares 反演 S_i 的紧凑参数 (k,q)；B_ij 已知 ----defresid_inv(theta):k_,q_=theta S_=k_*(x_c/500.0)**q_ out=np.zeros(0)forte,p_expinzip(t_exp,PSD_exp):p_mod=solve_pbm(S_,B,N0,te).y[:,-1]out=np.concatenate([out,(p_mod-p_exp)/(np.max(p_exp)+1e-9)])returnout res=least_squares(resid_inv,x0=[0.5,0.8],bounds=([0,0],[5,3]))k_est,q_est=res.xprint(f"反演: k={k_est:.3f}(真{k_ref}) q={q_est:.3f}(真{q_ref})")# ---- R² 与 AIC ----SSres=np.sum(res.fun**2);SSTot=np.sum((PSD_exp-PSD_exp.mean())**2)R2=1-SSres/SSTot n=PSD_exp.size;kk=len(res.x)AIC=n*np.log(SSres/n)+2*kkprint(f"R²={R2:.4f}AIC={AIC:.2f}(RSS={SSres:.2e}, k={kk}个参数)")# ---- 参数误差带 (cov) ----J=res.jac cov=np.linalg.inv(J.T@J)*(np.sum(res.fun**2)/(n-kk))perr=np.sqrt(np.abs(np.diag(cov)))print(f"k ±{perr[0]:.3f}q ±{perr[1]:.3f}")assertabs(k_est-k_ref)<0.1andabs(q_est-q_ref)<0.1,"反演应还原真参数"print("OK：least_squares 从合成数据还原真参数，误差带合理")

逐行要点：

• 合成实验：先用真参数k_ref/q_ref造出PSD_exp（solve_pbm在每个观测时刻取y[:,-1]），作为"标定目标"。真值已知，才能验收反演是否还原。
• least_squares的目标resid_inv返回一个展平的残差向量：把(时刻×bin)的偏差全部拼起来，用np.max(p_exp)归一（多输出的合并）。这正是破碎 PBM 这类"非标量"目标用least_squares而非curve_fit的原因。
• bounds=([0,0],[5,3])约束参数非负且q合理范围，防止最优化找负值/荒谬值。
• R²/AIC：res.fun是最优残差，R²=1−SSres/SSTot，AIC=n·ln(RSS/n)+2k——k=len(res.x)=2只惩罚 2 个自由参数，这就是"紧凑参数化"让 AIC 有意义的关键（自由度高会劣化 AIC）。
• 误差带：cov≈(JᵀJ)⁻¹·s²，s²=RSS/(n−k)，参数标准差sqrt(diag(cov))——J=res.jac。assert把"反演还原真参数"变成硬验收（合成数据里应该准）。
• 这也呼应篇 10/13 铁律：标定后的参数要回填到 PBM，再走一遍守恒+行为断言，确认标定没把模型"标歪"（本轮最优化只改参数、不改结构）。

三、常见报错与排查

least_squares不收敛 / 参数乱跑。
现象：res.success=False或k_est/q_est远超真值。根因：目标残差尺度不当（不同时刻/bin 量级差太多，被大项主导）、初值远离最优、或无边界。解法：残差按每档最大值归一、给合理bounds、x0用经验初值（如k≈0.3~1/ q≈0.6~1）；加分max_nfev。

R² 很高但 AIC 也高，拿不准。
现象：R²>0.99但加了参数 AIC 反而升。根因：这是典型"过拟合"信号——参数多了把噪声也拟合进去，R² 涨但 AIC 惩罚高。解法：换更紧凑参数化（减少自由k），比较不同模型的 AIC；选 AIC 最低者，而不是 R² 最高者。

cov出现负值或inf。
现象：sqrt(diag(cov))报错或给 NaN。根因：J.T@J奇异（参数完全相关/数据欠定），或n−k为 0/负。解法：检查秩（np.linalg.matrix_rank(J)），欠定则减少待估参数（用 GGSk/α而非每个 B_ij）；n>k是 AIC/cov 成立的前提。

把每个B_ij都当自由参数，反演爆炸。
现象：自由度比数据点还多，least_squares过拟合、误差带巨大。根因：B_ij是 n×n=1600 个，数据只有几十个点，完全欠定（底账 G 反算").解法：绝不逐格估B_ij——用 GGSk/α+ 分支比β的低维参数化约束它的形状，让自由度降到个位数。

标定后的参数回填 PBM，守恒断了。
现象：把估算参数塞回篇 10 PBM 后drift>1e-6。根因：最优化只改参数不改结构，理论上不应破坏守恒——但若标定时改了B的归一或网格，就会失衡。解法：标定函数里也复用同一套B构造（保持Σ_i B_ij=1）与同一网格，回填后再跑守恒断言确认没被"标歪"。

四、动手练习

1. 还原验收（必做）：把k_ref/q_ref改成别组真值（如(1.2,0.7)、(0.5,1.1)）各跑一遍合成-反演。判定：每次k_est/q_est都还原到真值 ±0.1，误差带覆盖真值，证明标定工具可信。
2. 噪声鲁棒：给PSD_exp加 5% 高斯噪声再反演。判定：k/q仍能大致还原（偏差放宽），误差带随噪声变大而变宽——理解"数据噪声会直接进误差带"。
3. 模型选型用 AIC：对同一合成数据，比较"2 参数幂律 S_i"与"3 参数（含偏移）S_i"的 AIC。判定：真参数是 2 参数的，3 参数模型 R² 略升但AIC 更高（过拟合）——用数据让"选 AIC 最低者"落地。

五、小结与下一篇预告

本篇完成了最后的"示意→标定"一跃：用least_squares从（合成）实验 PSD 反演S_i紧凑参数，用 R² 看拟合、AIC 防过拟合、cov给误差带。你掌握了参数标定的完整闭环——目标函数＋紧凑参数化＋统计选模＋误差带，它让此前的所有动力学变得可溯源、可复现、可审计。

下一篇：篇 18《灵敏度分析实战》——标定完参数，要问"哪些 CPP 最要命"：用 OAT/Morris 与 Sobol（SALib）定位对 CQA（细粉/D50）最敏感的关键工艺参数，排序出 Top3。

本篇认知问题回显（FAQ）

Q1：为什么必须"优化"而非"手调"？
A1：真实参数来自实验，需要一个可复现目标函数把"模型预测 vs 实测"最小化，min_θ Σ(model(θ)−exp)²。手调几万组合不现实且受主观偏好影响，最优化才有可复现的目标函数与统计指标，让标定可验证、可审计。

Q2：curve_fit和least_squares怎么分工？
A2：curve_fit适合"稀疏参数、标量输出、f(x,*θ)“的标准回归；破碎 PBM 的观测是"整个离散 PSD”（多输出、逐时刻×bin），所以首选least_squares——把(时刻×bin)偏差展平成一个大残差向量，并支持逐残基加权与边界约束。

Q3：怎么一次反演出S_i、B_ij，自由度怎么控制？
A3：层叠逼近：先低细度批次估S_i（主控减质量），再用least_squares一起优化；B_ij绝不给每个自由参数（欠定），用GGSk/α+ 分支比β的低维参数化约束形状。总自由度降到个位数，最优化才稳健、不过拟合、能过 AIC。

Q4：AIC 和 R² 各自惩罚什么，为什么 R² 高不够？
A4：R²=1−SSres/SStot只看拟合贴不贴、不惩罚参数数（加参数必涨）；AIC=n·ln(RSS/n)+2k显式惩罚参数个数，权衡"拟合好 vs 参数省"，越低越好。所以 R² 高而 AIC 也高 = 用过多参数堆出的"虚好"。选模看 AIC 更低者。

Q5：拟合参数"准不准"怎么量化（误差带）？
A5：least_squares/curve_fit在一阶近似下给cov≈(JᵀJ)⁻¹·s²（J 雅可比、s²=RSS/(n−k)），参数标准差sqrt(diag(cov))。可套性前提是残留误差尺度合理、模型近线性；若残差大或非线性强误差带会低估——所以还要看残差形状是否随机，别只信数字。