PharmaPy 干法制粒二次开发教程（5）：自定义单元操作基座——仿照 Crystallizers.py 写第一个可运行的自定义单元并接入 flowsheet

PharmaPy 干法制粒二次开发教程（5）：自定义单元操作基座——仿照 Crystallizers.py 写第一个可运行的自定义单元并接入 flowsheet

版本声明块

• 工具/软件：CryPTSys/PharmaPy（参考实现 BatchCryst / SolidPhase）
• 语言/环境：Python 3.9 + numpy<2 + scipy<1.14，积分用 scipy.integrate.solve_ivp
• 本文目标：用框架里"长得最像固相 PBM"的参考实现（Crystallizers.BatchCryst + Phases.SolidPhase）做模板，写出本系列第一个真正能跑的最小可行自定义单元，并把它挂进 SimulationExec 的 flowsheet 解出结果——为今后所有干法制粒单元立起"最小契约"样板。

一句话结论：在 PharmaPy 上新增干法制粒单元的最小做法是仿照Crystallizers.BatchCryst写一个类（本系列二次开发示例，如RollerCompactor）：声明Phases/Kinetics属性、在求解方法内实现破碎 PBM（dN_i/dt = −S_i·N_i + Σ_j B_ij·S_j·N_j）、再挂到SimulationExec的拓扑串IN-->X-->OUT上求解。框架"期待"你的类长什么样、怎么装配、怎么取结果，是这篇要讲透的。

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

• 官方框架里没有辊压机/粉碎机/分级机，那我要新增一个单元操作，从哪里"抄"数据结构最省力气而不踩官方 API 的雷？
• 一个能被sim.xxx = ...挂上去、能被SolveFlowsheet(kwargs_run=...)正常求值的 PharmaPy 单元类，到底必须满足哪几条"硬"契约？
• SolidPhase把粒径分布（PSD）承载成x_distrib和distrib两个量，我自己的破碎单元该如何读写它们，才能保证接口不跟框架断掉？
• 为什么只写一个"动力学退化成简单破碎 PBM"的骨架单元就够？最小可行单元的边界在哪里？
• 挂上拓扑串之后我要怎么把结果拿回来——result.time、result.字段、GetStreamTable是怎么一步步接通的？

一、机制解析

我们先直面一个事实，这也是本系列反复立论的核心：CryPTSys/PharmaPy 的公开仓库（master/develop）只覆盖原料药侧（drug substance）单元操作——反应器、结晶器、蒸发、蒸馏、萃取、干燥、固液分离、三相沉降；没有任何固体制剂（drug product）单元，feeder、blender、roller compactor、mill、tablet press 全缺。干法制粒三件套在框架里是空白，必须二次开发。这恰恰说明：本系列的难点不在"调模型"，而在"学会在框架的地基上把一个新的单元操作类立起来"。而本篇要做的，就是把这条扩展链路第一次完整打通。

第一步是找准"模板源"。扩展框架切忌凭空发明接口。我们翻一翻框架自带包里谁最接近固相 PBM：PharmaPy.Crystallizers.BatchCryst就是用SolidPhase承载粒径分布的结晶器，它内部同时支持method='1D-FVM'（离散 bin 的有限体积）和'moments'（各自矩），这两条正是离散 PBM 的两种落法；而PharmaPy.Phases.SolidPhase/MixedPhases.Slurry就是把它那个"双表示"——x_distrib（粒径网格 [µm]）加distrib（数目密度分布 [#/m³/µm]）——和一组矩moments同时存起来的数据容器。所以，我们新单元的RollerCompactor/Comminution/Classifier（本系列二次开发示例类名）就要以BatchCryst+SolidPhase为参照系对接口，而不是另起炉灶。

再说类的最小契约。综合官方 guide（SimulationExec、单元挂载、SolveFlowsheet）的写法，一个能进 flowsheet 的 PharmaPy 单元，至少得满足这么几条"硬"约定：

1. 构造签名是"数据/参数不是容器"。像sim.R01 = PlugFlowReactor(diam_in, num_discr=50)这样，构造时传入数值/数组参数，决不直接往里塞Phases。这也意味着初始 PSD 不是构造参数，而是构造之后再注入。
2. 暴露属性名要能被框架识别。框架靠sim.单元名.Phases/.Kinetics/.Inlet/.Utility这些固定属性名来做装配。所以我们的新类必须定义self.Phases = None这样的属性占位，并在(kw)init之外提供赋值入口。
3. 要给self.Phases放一个真的承载 PSD 的对象，也就是SolidPhase(path_thermo, x_distrib=..., distrib=...)。它既是状态容器，也是我们读取初始粒径分布、回写演化后分布的通道。
4. 提供可供求解器调用的核心方法。框架的SolveFlowsheet(kwargs_run=...)会给每个单元准备一组运行参数（如runtime），然后调用单元内部的积分求解例程。一个对齐到BatchCryst风格的做法是：写一个类内方法（我们姑且称solve_unit，本系列二次开发示例方法名，非官方 API），内部用scipy.integrate.solve_ivp把这个单元自己的 ODE 系统（这里是破碎 PBM）积分完，返回一个结果对象。
5. 返回结果对象并让它与官方结果结构兼容。官方把结果封装成DynamicResult：DynamicResult.__dict__.update(**results)，于是result.time是时间向量、result.<字段>直接取各状态；而sim.xxx挂载的各单元状态会汇总进SimulationResult，用GetStreamTable(basis='mass')输出 pandas DataFrame。我们不硬造一套新的结果抽象，而是复用这套对象，把单元的状态塞进去。

为什么最小可行就够？你可能会问：这个骨架单元的动力学很简陋，真的有意义吗？有。这一篇的交付目标不是动力学精度，而是把"类结构 → 装配 → 接入 flowsheet → 取结果"这条链路第一次跑通。类写对了、挂得上、解得出、取得到，后面（篇 6-12）我们只是把solve_unit里的动力学内核换得更物理：线条从"占位 PBM"换成真正由篇 6-7 的 Johanson/Midoux、篇 9-10 的 Vogel–Peukert + 离散破碎 PBM、篇 12 的 Tromp 分级驱动的内核。骨架搭错，后面全是白忙；骨架搭对，剩下的就是填内核。这套"最小可行单元"就是我们整个干法制粒模块的原型。

把这套扩展逻辑画成图，就是figs/fig12_arch.png：底层是官方框架的可复用基座（SolidPhase的 PSD 承载、SFlowsheet拓扑执行器、DynamicResult/GetStreamTable结果层），我们在基座上新增一组"干法制粒单元"（辊压—破碎—分级，本系列二次开发示例），它们复用基座的 PSD 数据结构与结果通道，只自创动力学内核。

一句话把机制钉住：框架不给你现成单元，但给你离固相 PBM 最近的参考实现（BatchCryst+SolidPhase）以及一套固定的装配/求解/取结果协议。你要做的，就是让新类"长得像"参考实现、满足协议，然后把动力学内核换成干法制粒的。

二、完整代码与逐行剖析

下面给出一个完整可运行的最小骨架。请注意：RollerCompactor、solve_unit、MYU01这些都是本系列二次开发示例的类名/方法名/拓扑单元名，不是 PharmaPy 官方 API；只有SimulationExec、SolidPhase、BatchCryst、SolveFlowsheet等是官方逐字签名。

importnumpyasnpfromscipy.integrateimportsolve_ivp# ============ 唯一的外部依赖：官方 PSD 承载类 ============fromPharmaPy.PhasesimportSolidPhase# ============ 对策 1：构造签名收"数据/参数"，不收 Phases ============classRollerCompactor:# ★ 本系列二次开发示例类（非官方 API）"""最小可行自定义单元骨架，模拟固相"辊压-破碎"。 method 固定选 'discrete'（离散 bin 破碎 PBM）。"""def__init__(self,path_phys,x_distrib,distrib,method='discrete',num_mom=4):# 记住工作粒径网格 [µm]，本系列示例骨架用"常量选择率"拟合后续 Vogel-Peukertself.x_distrib=np.asarray(x_distrib,dtype=float)self.num_bins=len(self.x_distrib)self.method=method# ★ 契约 2：暴露框架能识别的属性名（占位，构造时置空）self.Phases=None# 将由外部注入一个 SolidPhaseself.Kinetics=None# 本骨架留空，后篇接 Johanson/Vogel-Peukert# 把初始 PSD 也暂存一份（供 PBM 出生项用）self.distrib0=np.asarray(distrib,dtype=float)# 本骨架使用"常量选择率 S_i"，物理上由篇6-12 补齐self.Selection=self._default_selection()def_default_selection(self):# ★ 自造方法"""占位选择函数：大粒更易破（一次方按粒径归一），本系列示例近似。"""x_ref=float(np.max(self.x_distrib))return0.05*(self.x_distrib/x_ref)# 单位 1/s，量级示意def_breakage_rhs(self,t,N,S,B):# ★ 自造方法（破碎 PBM 右端）"""离散破碎 PBM 的 ODE 右端： dN_i/dt = Σ_{j>=i} S_j*B_ij*N_j − S_i*N_i （birth−death） N 为各 bin 的数目密度向量；同时累计总质量用于守恒校验。"""birth=np.einsum('ij,j,j->i',B,S,N)# 出生项：来自更大/同级粒death=S*N# 死亡项returnbirth-deathdef_build_B(self):# ★ 自造方法（破碎函数矩阵，质量守恒近似）"""构造下三角归一破碎矩阵 B_ij（j 破碎到 i 的分数）。 用最简"等比均分到至少两个子bin"，本系列骨架示例，后篇换 GGS/RRB。"""nb,B=self.num_bins,np.zeros((self.num_bins,self.num_bins))forjinrange(1,nb):# 母粒 bin ji_lo=max(0,j//2)# 子粒去向：j 与 j//2 两个相邻 binw1=(self.x_distrib[j]-self.x_distrib[i_lo])/self.x_distrib[j]B[i_lo,j]=w1# 分配到 i_lo 的质量分数B[j,j]=1.0-w1# 剩余留在 j（粗粒近似保持守恒）returnBdefsolve_unit(self,runtime,**opt):# ★ 自造方法（非官方 API）"""按 BatchCryst 风格：读 Phases 的 PSD → 积分解 → 回写并返回结果。"""ifself.PhasesisNone:raiseRuntimeError("RollerCompactor.solve_unit 需要先注入 .Phases")N0=self.Phases.distrib.copy()# 取注入的 SolidPhase 的初始分布B=self._build_B()S=self.Selection sol=solve_ivp(self._breakage_rhs,(0.0,runtime),N0,method='BDF',rtol=1e-5,atol=1e-9,args=(S,B),dense_output=False)# 质量守恒自检（相对漂移应 <1e-3，作为骨架硬门槛）mass0=np.sum(N0*self.x_distrib**3)mass1=np.sum(sol.y[:,-1]*self.x_distrib**3)drift=abs(mass1-mass0)/mass0assertdrift<1e-3,f"破碎 PBM 质量漂移{drift:.2e}超容差"# ★ 对策：回写演化后的分布，让 Phases 携带"本地状态"self.Phases.distrib=sol.y[:,-1].copy()# ★ 结果结构对齐：仿 DynamicResult.__dict__.update(**results)classFakeResult:# ★ 自造简化结果壳，仅演示字段访问passres=FakeResult()res.time=sol.t.copy()# → result.timeres.N=sol.y.copy()# → result.Nres.defects_num_of_phases=0res.defects={'num_of_phases':0,'mass_drift':drift}returnres

装配并求解整个 flowsheet：

# ============ 对策 3：装配 + 接入 flowsheet ============importnumpyasnpfromPharmaPy.SimExecimportSimulationExecfromPharmaPy.PhasesimportSolidPhase# 粒径网格：与官方结晶器一致的对数等距，35 bin（本教程亦可取 30-50）x_gr=np.geomspace(1,1500,num=35)# [µm]sim=SimulationExec(path_phys,flowsheet='IN --> MYU01 --> OUT')# 拓扑串# 构造单元（样例 stat）：构造签名收数据不收 Phasessim.MYU01=RollerCompactor(path_phys,x_gr,distrib0,# ★ 示例单元名 MYU01method='discrete')# 与示例类名一致，非官方# 注入初始 PSD：这才把分布交到单元手上sim.MYU01.Phases=SolidPhase(path_phys,x_distrib=x_gr,distrib=distrib0,mass_frac=[0,0,1,0,0])# （骨架单元暂无 Kinetics/Inlet/Utility 注入项，仅挂 Phases 即可求解）

跑通并把结果取回来：

run_kwargs={'MYU01':{'runtime':30.0}}# 求解时间窗（秒）sim.SolveFlowsheet(kwargs_run=run_kwargs)# ★ 官方流程：按单元给运行参数# 取结果三种通道（与官方结果结构对齐）res=sim.MYU01.resultsifhasattr(sim.MYU01,'results')elseNonet=res.timeifresisnotNoneelseNone# → 时间向量N_last=res.N[:,-1]ifresisnotNoneelseNone# → 末时刻各 bin 数目密度table=sim.GetStreamTable(basis='mass')# → 全流程流表（pandas MultiIndex）

逐行要点：SimulationExec是官方执行器，flowsheet='IN --> MYU01 --> OUT'是官方拓扑字符串语法（箭头用-->），我们要把示例单元MYU01嵌在中间；SolidPhase的构造吃x_distrib+distrib，正好对应我们单元内部使用的网格；sim.单元名.方法(...)属性注入是官方装配范式，我们让RollerCompactor提供Phases占位以对号入座；solve_unit内部用solve_ivp(method='BDF')积分离散破碎 ODE，出生项用np.einsum('ij,j->i', B, S, N)双重求和、死亡项S*N，完全对应离散破碎 PBM 的 birth-death 结构；末了做质量守恒断言，再把演化分布回写self.Phases.distrib——这样你连"单元内部状态"都能在 flowsheet 语境里查得到。

三、常见报错与排查

报错一：AttributeError: 'RollerCompactor' object has no attribute 'Phases'（或挂到sim.MYU01前报未注入）。
现象：一给单元赋值sim.MYU01 = RollerCompactor(...)或在SolveFlowsheet前访问就报缺属性。根因：类里没在__init__占位self.Phases = None，框架装配时找不到官方约定的属性名。解法：在构造方法里把所有框架约定属性（Phases/Kinetics/Inlet/Utility）先占位成None，再通过属性注入赋值。

报错二：拓扑串里名字对不上，SolveFlowsheet说找不到MYU01。
现象：flowsheet='IN --> MYU -> OUT'写错一个字母，或kwargs_run={'MYU01':...}与挂载名不一致，结果求解器直接抛 KeyError/找不到单元。根因：拓扑串、sim.单元名、kwargs_run键三者必须严格同名同写。解法：三步统一——flowsheet里的节点名、sim.那个名 = ...、kwargs_run的字典键，三处用同一个字符串常量，别手打。

报错三：solve_ivp返回失败或distrib维度不匹配。
现象：BDF 一步都走不动、返回RuntimeError: solver，或 einsum 报维度错误。根因：初始N0=self.Phases.distrib与self.x_distrib的 bin 数不一致（len(distrib) != len(x_distrib)），或distrib里含负值/零破坏数值。解法：构造时校验len(distrib)==len(x_distrib)，初始分布强制非负；积分泛型容差收到rtol=1e-5, atol=1e-9；对骨架 PBM 足够。

报错四：没实现求解方法，SolveFlowsheet求值时抛TypeError或"object is not callable"。
现象：单元能创建、能挂载，但一进SolveFlowsheet(kwargs_run=...)就报"类型错误，无法调用"。根因：框架在 flowsheet 求值时需要一个可调用的求解入口（对应BatchCryst内部的求解例程），而你的类只写了solve_unit或连它也没写，名字/签名对不上框架的调用约定。解法：确认类里实现的是名字风格对齐BatchCryst的求解方法（本骨架叫solve_unit，自造示例名），并在方法内真正调用solve_ivp并返回结果对象；别只留 stub。

报错五：结果对象是空的，res.time/GetStreamTable取不到值。
现象：仿真"跑完"但sim.单元名.results为空，或GetStreamTable的 DataFrame 缺该单元的行。根因：单元没把状态写进结果结构——要么没返回结果对象，要么没在result.time/字段上赋值，要么该单元与官方结果聚合逻辑没有接上。解法：对照本篇结果壳做法，保证方法返回的对象有time和状态字段（result.字段直接可取值）；需要进流表时，把单元的出口流写成官方SolidPhase/Stream兼容结构并由GetStreamTable汇总。

四、动手练习

练习一（契约落地）：把上面的骨架类抄下来，改一个你自己的名字（如MyBreakageUnit，本系列二次开发示例），按篇 2 的环境在官方SimulationExec(flowsheet='IN --> X --> OUT')里挂载一个SolidPhase初始 PSD 并SolveFlowsheet。客观判定标准：能正常求值到runtime=30，且从返回结果拿到res.time与末时刻res.N[:, -1]，无缺属性/缺方法报错。

练习二（守恒门槛）：给骨架加上"质量随 bin 的 x³ 加权求和全程守恒"的自检，跑runtime拉长到 300 s。客观判定标准：每一步积分末的drift<1e-3；若超容差，回到_build_B检查母粒破碎后子粒质量是否漏分配（这是最常出错的地方）。

练习三（结果通道）：把单元出口连成官方可读流，调用sim.GetStreamTable(basis='max mass')打印该单元所在行的质量流与各 bin 分布。客观判定标准：表里出现你的单元名，且其质量流与sum(N * x³)的对账误差 <1%（提示：骨架用basis='mass'并检查 MultiIndex 各级索引是否齐全）。

五、小结与下一篇预告

本篇我们第一次在 PharmaPy 上立起了一个属于干法制粒的最小可行单元：以官方Crystallizers.BatchCryst+Phases.SolidPhase为模板源，兑现了"构造签名收数据、属性占位可注入、solve_unit内积分解破碎 PBM、结果进DynamicResult/GetStreamTable"整套最小契约的演示（含本系列二次开发示例的RollerCompactor骨架与my_breakage内核）。动力学现在是简化的，但链路是通的——这正是接下来所有干法制粒单元的落点。

下一篇（篇 6）进入辊压模块①：我们把这条链上的第一块物理内核——Johanson 滚压理论——落进来，讲咬入角（nip angle）与滑移区/咬入区的应力分析、压力-密度关系，以及"直接套用会高估 ribbon 密度、未计入弹性回弹"这个经典局限，为后篇的 Midoux 无量纲数与 Reynolds 模型做铺垫。

本篇认知问题回显（FAQ）

Q1：官方没有干法制粒单元，我新增单元该从哪里"抄"数据结构最合适？
A：从框架里"离固相 PBM 最近"的实现抄——Crystallizers.BatchCryst（支持method='1D-FVM'离散 bin 与'moments'两种求解法）和Phases.SolidPhase/MixedPhases.Slurry的"双表示"数据容器（x_distrib粒径网格 [µm] +distrib数目密度分布 [#/m³/µm] + 一组矩moments）。扩展单元复用这套 PSD 承载结构，只自创动力学内核，是最省力也最不踩官方 API 雷的路径。

Q2：能被sim.xxx挂上、能被SolveFlowsheet求值的单元类必须满足哪几条硬契约？
A：① 构造签名收"数据/参数"而非 Phase 容器；② 暴露官方约定属性名（Phases/Kinetics/Inlet/Utility），在__init__先占位为None；③ 提供一个对齐BatchCryst风格的求解方法（本骨架叫solve_unit，自造示例名），内部用scipy.integrate.solve_ivp积分真正返回结果对象；④ 结果对象结构对齐DynamicResult（result.time+ 状态字段），出口流兼容官方SolidPhase/流对象以便GetStreamTable汇总。

Q3：SolidPhase把 PSD 承载成x_distrib和distrib两个量，破碎单元怎么读写才能不断接口？
A：两量必须len相等、分布非负；构造时从x_distrib建立自己的 bin 网格，从注入的self.Phases.distrib复制初始N0，积分完把末态回写self.Phases.distrib，这样单元内部状态在 flowsheet 语境也可查；任何质量换算都用Σ N·x³加权（单位对齐 µm）。

Q4：为什么只写一个动力学退化成简单破碎 PBM 的骨架单元就够？最小可行单元的边界在哪？
A：本篇目标是打通"类结构→装配→接入 flowsheet→取结果"链路，不是求动力学精度。骨架把dN_i/dt = −S_i·N_i + Σ_j B_ij·S_j·N_j的 birth-death 结构、单位约定、守恒自检、结果通道全部跑通；动力学内核是唯一"最简/占位"部分（常量选择率 + 均分子粒），后篇逐步替换为 Johanson/Midoux（辊压）、Vogel–Peukert+离散破碎 PBM（粉碎）、Tromp（分级）。骨架搭对，填充内核即可；这是本项目所有干法制粒单元的原型。

Q5：挂上拓扑串后如何把结果拿回来？
A：三条通道：① 单元返回的结果对象 →res.time时间向量、res.字段状态；② 把单元出口连成官方可读流 →sim.GetStreamTable(basis='mass')输出包含该单元的 pandas MultiIndex 流表（temp/pres/mass/vol_flow 等列）；③ 单元内部状态可查self.Phases.distrib回写值。三者都与官方DynamicResult.__dict__.update(**results)的结果结构对齐。