连续粉体进料与混合建模教程4第一个进料混合串联仿真——阶跃扰动下的最小可跑闭环版本声明块工具/软件numpy卷积、scipyscipy.signal.fftconvolve快速卷积、matplotlib阶跃响应图实际产线建模可再接入 PharmaPyPurdue CryPTSys作为 flowsheet 骨架。语言/环境Python 3.9numpy2、scipy1.14中文字体使用SimHei或Microsoft YaHei。本文全部数值为教学示意参数非真实产线实测数据。本文目标从零写出第一个进料 → 混合 → 出口浓度的最小串联仿真把卷积这条扰动传播的唯一法则亲手跑通一遍。一句话结论一个最小可跑的进料-混合串联模型只需要三个部件——进料器生成质量流量序列、混合器用np.convolve(feed_profile, rtd) * dt做 RTD 卷积、出口浓度输出在入口阶跃扰动下出口会呈现先迟滞由 PFR 延迟 θ 造成后平滑由返混造成的 S 形响应且只有当sum(E)*dt 1即∫E(t)dt 1时稳态出口浓度才等于入口浓度。〇、本篇要解决的认知问题前 3 篇都在讲清楚这一篇开始动手。五个必须先想明白的问题一个进料-混合串联系统的最小可跑模型最少需要哪三个部件各自负责什么卷积到底在算什么为什么串联 卷积而不是串联 相加为什么在阶跃扰动下出口会同时出现迟滞与平滑两个特征它们分别来自哪里离散时间卷积和连续卷积差在哪里采样步长dt该怎么选∫E(t)dt 1为什么是数值陷阱长序列卷积为什么要用scipy.signal.fftconvolve一、机制解析1.1 最小模型只需三个部件连续上游段的动态最核心的一条链是进料器把物料送出去混合器把扰动抹平出口浓度随之变化。把它写成最小模型┌───────────────┐ 质量流量序列 ┌───────────────┐ 出口浓度 ┌──────────┐ │ 进料器 │ ────────────────► │ 混合器 │ ──────────► │ 出口输出 │ │ feed_profile │ c_in(t) / F(t) │ RTD 卷积 │ c_out(t) │ 记录/判据 │ │ 生成质量流量 │ │ E(t) │ │ │ └───────────────┘ └───────────────┘ └──────────┘ ① ② ③部件职责输入输出① 进料器生成入口质量流量/浓度序列设定值、扰动feed_profile随时间变化② 混合器用 RTD 把入口序列变换为出口序列feed_profile、E(t)c_out(t)③ 出口输出记录与判定c_out(t)曲线 / 是否超规格为什么这对你重要把系统拆成生成 → 变换 → 输出三段是通用套路。进料器负责造扰动混合器负责变扰动输出负责判扰动。你后面写PowderFeeder、ContinuousBlender类时本质就是给①和②各写一个更真实的实现接口顺序不变。1.2 卷积串联系统的组合法则先给定义。连续卷积是c_out(t) ∫₀ᵗ c_in(t − x) · E(x) dx E 为 RTD量纲 1/时间两个 RTD 的串联组合也是卷积E₁ ∗ E₂ ∫₀ᵗ E₁(x) · E₂(t − x) dx这就是扰动沿产线传播全部计算的基石。它的直观含义是加权求和出口在t时刻的浓度是入口在过去所有时刻t−x的浓度各自乘上这段物料在设备里刚好待了x秒的概率E(x)再加起来。手算一下最粗糙的版本把时间切成 3 个格子入口 c_in [0, 1, 1] RTD 采样 E [0.5, 0.5, 0] t0: c_out 0*0.5 0 t1: c_out 0*0 1*0.5 0.5 t2: c_out 0*0 1*0 1*0.5 0.5 ← 注意每格还要乘 dt 才是积分为什么串联 卷积而不是串联 相加因为串联不是把浓度叠加而是把停留时间叠加。物料穿过第一台设备的耗时是一个随机变量穿过第二台又是一个随机变量两个独立的随机变量之和的分布正是二者分布的卷积。相加是两份物料混在一起卷积是同一份物料先后穿过两台设备——两者物理含义完全不同。1.3 迟滞与平滑同一个响应里的两个特征在出口曲线上你会同时看到两个特征它们来自两种不同的物理机制特征表现来源数学对应迟滞延迟出口比入口晚了一段时间才开始变化设备的活塞流段物料整体被推着走先到先出PFR 的延迟时间θ平滑返混出口不是阶跃而是缓慢的 S 形爬升设备的混合段颗粒停留时间有分布CSTR 的返混τ_tank一个常用的教学 RTD 把两者叠在一起PFR 延迟 两个等体积 CSTRE(t) (t − θ) · exp[ −(t − θ)/τ_tank ] / τ_tank² t θt ≤ θ 时为 0 总 MRT θ 2·τ_tank第一项(t − θ)保证了t θ之前出口完全为 0——这就是迟滞而指数项带来的拖尾就是平滑。入口阶跃 出口响应 │ ┌────── │ ,------ ← 平滑返混决定爬升形状 │ │ │ ,- │ │ │ ,- │──────┘ │,- ← 迟滞 θ 决定何时开始动 └──────────► t └──────────────► t t₀ t₀θ为什么这对你重要这两个特征对应两种完全不同的工程对策。迟滞主要由输送路径长度与流量决定改不了多少而平滑程度由混合器持料量与转速决定可调。分不清这两者你调参时就会按错旋钮。1.4 离散时间卷积 vs 连续卷积×dt不能漏计算机只能算离散序列。连续卷积∫ c_in(t−x)E(x)dx用离散近似是c_out[n] Σₖ c_in[n−k] · E[k] · dt dt 是采样步长三件事必须记住×dt是积分变离散的必备因子。漏掉它输出会整体放大1/dt倍——这是最隐蔽的量纲错误。E必须按dt采样。你构造的是E(t)在t 0, dt, 2dt, …上的取值数组而不是若干个概率值。采样步长dt的选择有硬约束dt必须显著小于被解析的最短时间尺度。经验法则是dt ≤ θ/10且dt ≤ τ_tank/10同时序列总长要覆盖 RTD 尾部一般取≥ 5×MRT否则尾部被截断ΣE·dt会小于 1。dt 太大 ──► 迟滞被抹平、爬升失真把 S 形算成折线 dt 太小 ──► 计算量暴涨且数值噪声被放大 dt 合适 ──► 能分辨 θ 与 τ_tank 的最小时间常数为什么这对你重要第 3 篇讲的时间尺度错配在这里第一次变成具体的dt选择。你现在选的dt直接决定后面所有卷积结果的可信度。1.5 归一化陷阱∫E(t)dt 1RTD 是一个概率密度必须满足∫₀^∞ E(t) dt 1数值上等价于Σ E[k] · dt 1。三个最常见的坑截断尾部时间轴太短RTD 的长尾被切掉积分 1出口稳态浓度会偏低漏乘dt归一化检查时忘了×dt看到Σ E远大于 1 就误以为自己构造错了离散化偏差用稀疏采样构造E(t)时梯形积分的误差在大dt下不可忽略。为什么要较真归一化不成立意味着物料不守恒。在连续制造里物料守恒是硬约束——你把 1 kg/h 的示踪剂打进产线出口累计就不可能超过 1 kg/h。归一化错误会直接让你在第 18 篇的物料追溯与分流计算里得出荒谬结论。1.6 用fftconvolve加速长序列直接卷积的复杂度是O(N·M)。当你为了解析毫秒级压片循环、或为了跑几百秒的产线动态而把序列拉到十万点以上时直接卷积会明显变慢。scipy.signal.fftconvolve用快速傅里叶变换把复杂度降到约O(N log N)序列越长加速比越大。一个反直觉的默认值必须提前说np.convolve的默认modefull输出长度是N M − 1比输入长出一截。如果你不做截断向下游传的序列会越来越长几级串联后时间轴就错位了。正确做法是modefull之后截断回输入长度或用modesame并核对对齐方式。二、完整代码与逐行剖析下面第一段脚本就是本系列的第一口活气端到端最小串联 阶跃响应。全部数值为教学示意参数。# -*- coding: utf-8 -*- 第 4 篇配套脚本 A进料-混合串联最小可跑闭环 阶跃响应。 依赖numpy, matplotlib。教学示意参数非真实产线数据。 importnumpyasnpimportmatplotlib matplotlib.use(Agg)importmatplotlib.pyplotasplt plt.rcParams[font.sans-serif][SimHei,Microsoft YaHei]plt.rcParams[axes.unicode_minus]Falsedefbuild_rtd(t,theta,tau_tank): 构造“PFR 延迟 theta 两个等体积 CSTR”的教学示意 RTD。 E(t) (t - theta) * exp(-(t - theta)/tau_tank) / tau_tank^2t theta否则 0。 总 MRT theta 2 * tau_tank。 Enp.zeros_like(t)maskttheta shiftt[mask]-theta E[mask]shift*np.exp(-shift/tau_tank)/tau_tank**2returnEdefsimulate_line(feed_profile,rtd,dt): 端到端最小串联出口 入口 * RTD离散卷积近似连续卷积须乘 dt。 feed_profile : 入口浓度/流量序列与 rtd 同一时间基准 rtd : E(t) 序列1/s按 dt 采样 dt : 采样步长s 返回与 feed_profile 等长的出口序列保证时间轴对齐。 fullnp.convolve(feed_profile,rtd,modefull)*dtreturnfull[:len(feed_profile)]# 截断回输入长度避免时间轴错位defmain():dt1.0# 采样步长 1 stnp.arange(0.0,400.0,dt)theta,tau_tank20.0,15.0# 教学示意延迟 20 s单釜 15 srtdbuild_rtd(t,theta,tau_tank)areartd.sum()*dtprint(fRTD 归一化检查: sum(E)*dt {area:.6f}应 ≈ 1)ifabs(area-1.0)1e-3:print(警告归一化偏差偏大请延长时间轴或减小 dt。)c_in(t100.0).astype(float)# 入口阶跃t100 s 起由 0 变 1c_outsimulate_line(c_in,rtd,dt)idx_halfint(np.argmax(c_out0.5))# 首次达到终值 50% 的时刻print(f入口阶跃时刻:{100.0:.0f}s)print(f出口达到终值 50% 的时刻:{t[idx_half]:.0f}s迟滞 爬升的合成)print(f出口稳态值末点:{c_out[-1]:.4f}应 ≈ 1前提是 RTD 归一化)fig,axplt.subplots(figsize(9,4))ax.plot(t,c_in,--,label入口浓度 c_in阶跃)ax.plot(t,c_out,-,linewidth2,label出口浓度 c_outS 形响应)ax.axvline(100.0,colorgray,linestyle:,linewidth1)ax.set_xlabel(时间 ts)ax.set_ylabel(归一化浓度-)ax.set_title(进料-混合串联的阶跃响应迟滞 平滑教学示意)ax.legend()fig.tight_layout()fig.savefig(step_response.png,dpi150)print(已生成 step_response.png)if__name____main__:main()逐段剖析build_rtd的第一行E np.zeros_like(t)是必须的。先生成全零数组再只给t θ的部分赋值自然保证了t ≤ θ时E 0——这就是迟滞的数值实现。若用向量化公式直接算t θ时(t − θ)为负、会产生负的概率物理上荒谬。simulate_line里的* dt是全篇最关键的一个乘号。漏掉它出口浓度会放大1/dt 1倍这里 dt1 恰好掩盖了错误——所以你必须用dt 0.5或2.0各跑一次来验证。这是本篇最反直觉的地方dt 1时漏乘dt和不漏乘结果一样会让你误以为代码是对的。full[: len(feed_profile)]的截断。np.convolve(modefull)输出长度是N M − 1多出M − 1个点。截断回N保证c_in和c_out在同一时间轴上一一对应。c_in (t 100.0).astype(float)用布尔数组生成阶跃。比循环赋值快得多而且不会漏点.astype(float)把True/False转成1.0/0.0。idx_half用np.argmax(c_out 0.5)。返回第一个为真的下标用来量化迟滞 爬升合成后的半程时刻。教学参数θ 20 s、τ_tank 15 s、阶跃在t 100 s下它约等于145 s 100 θ 1.68·τ_tank必然大于100 θ 120 s——出口还要花约 25 s 爬升这个数字比单纯的θ更能说明迟滞与平滑是叠加的。第二段脚本演示归一化陷阱与fftconvolve加速并顺手给出第 10 篇会用到的 n-CSTR 形状 RTD。# -*- coding: utf-8 -*- 第 4 篇配套脚本 B归一化陷阱与 fftconvolve 加速。 依赖numpy, scipy。教学示意参数。 importtimefrommathimportfactorialimportnumpyasnpfromscipy.signalimportfftconvolvedefrtd_gamma(t,tau,n): n 个等体积 CSTR 串联的无量纲形状 RTD教学示意深入讨论见第 10 篇。 E(t) n^n / (n-1)! * (t/tau)^(n-1) * exp(-n * t/tau) / tau Enp.zeros_like(t)maskt0xt[mask]/tau E[mask](n**n)/factorial(n-1)*x**(n-1)*np.exp(-n*x)/taureturnEdefnormalize_rtd(E,dt):把离散 RTD 归一化使 sum(E)*dt 1。返回 (归一化后序列, 原始积分)。areaE.sum()*dtreturnE/area,areadefmain():dt0.5tnp.arange(0.0,300.0,dt)Ertd_gamma(t,tau20.0,n3)E_norm,areanormalize_rtd(E,dt)print(f归一化前 sum(E)*dt {area:.6f})print(f归一化后 sum(E)*dt {E_norm.sum()*dt:.6f}应 ≈ 1)# 陷阱演示np.convolve 默认 modefull输出长度 N M - 1N,M200000,2000rngnp.random.default_rng(0)arng.random(N)brng.random(M)print(f\n输入长度 N{N}RTD 长度 M{M})t0time.time()conv_directnp.convolve(a,b,modefull)t_directtime.time()-t0 t0time.time()conv_fftfftconvolve(a,b,modefull)t_ffttime.time()-t0print(fnp.convolve 耗时{t_direct*1000:8.1f}ms输出长度{len(conv_direct)})print(ffftconvolve 耗时{t_fft*1000:8.1f}ms输出长度{len(conv_fft)})print(f两者最大绝对差 {np.max(np.abs(conv_direct-conv_fft)):.3e})print(f加速比 ≈{t_direct/max(t_fft,1e-9):.1f}x序列越长越明显)if__name____main__:main()逐段剖析rtd_gamma的分母是factorial(n - 1)而不是factorial(n)。n-CSTR 的标准形状是E(θ) n^n/(n−1)! · θ^(n−1) · exp(−nθ)再除以 τ 换成有量纲的E(t)。写错分母会让积分不等于 1——正好被第二段代码的归一化检查抓出来。normalize_rtd返回原始积分是为了让你看到截断或采样不足到底差了多少。工程上习惯先归一化再使用但先看原始积分是纪律如果原始积分只有 0.7那说明你的时间轴或dt根本不够归一化只是掩盖了问题。np.random.default_rng(0)用固定种子保证你我的加速比测试跑在同样的输入上结果可复现。这是性能对比脚本的通用写法。fftconvolve与np.convolve结果应几乎相同差异在1e-12量级浮点误差。如果你看到显著差异几乎一定是边界处理或长度不匹配的问题。time.time()做基准测试是粗粒度手段适合这种差一个数量级的对比要精确测时请用time.perf_counter。三、常见报错与排查1现象出口稳态浓度不等于入口浓度比如入口阶跃到 1出口只到 0.83。根因归一化不成立——Σ E · dt ≠ 1通常是时间轴太短、RTD 尾部被截断。解法打印rtd.sum() * dt把时间轴延长到≥ 5×MRT或减小dt。归一化前先看原始积分别直接用归一化掩盖截断问题。2现象出口曲线整体比预期高一到几个数量级。根因离散卷积漏乘dt。连续卷积∫c_in·E·dx近似为离散和时必须×dt。解法在np.convolve(...)后乘dt并用dt 0.5与dt 2.0各跑一次验证结果的dt无关性只用dt 1是验证不出这个错误的。3现象多级串联后出口序列越来越长和输入的t轴对不上。根因np.convolve默认modefull输出长度是N M − 1。解法卷积后截断回输入长度full[:len(x)]或改用modesame并统一对齐规则每级串联都把长度钉死在N。4现象阶跃响应曲线看起来是折线没有明显的 S 形平滑段。根因采样步长dt太大把返混造成的爬升抹平成了直线段。解法按dt ≤ θ/10且dt ≤ τ_tank/10收紧步长教学示例里θ 20 s、τ_tank 15 s故dt应取1 s或更小。5现象长序列卷积跑得很慢几百秒的产线动态要等很久。根因直接用O(N·M)的np.convolve序列上十万点时耗时明显。解法改用scipy.signal.fftconvolve约O(N log N)序列越长加速比越大。注意两种实现的边界与长度保持一致。四、动手练习练习 1基础跑通最小闭环并核对三个数。把脚本 A 存为first_lineup.py运行。成功标准量化打印的sum(E)*dt落在1.000 ± 0.005出口稳态末值落在0.98–1.02出口达到终值 50% 的时刻落在145 ± 5 s。练习 2进阶验证×dt的必要性本练习是本篇核心。把dt改成0.5与2.0各跑一次先故意删掉simulate_line里的* dt观察结果。成功标准量化删掉* dt后dt0.5时出口稳态约为2.0、dt2.0时约为0.5即偏离真值1/dt倍把* dt加回来后两种步长下出口稳态都落在0.98–1.02。练习 3挑战给串联加一段并核对整线 MRT。把simulate_line改成两段串联先用build_rtd生成两段 RTD参数可自定如θ₁5, τ₁8与θ₂10, τ₂12先分别卷积、再把两段 RTD 卷积后再卷积入口。成功标准量化两种算法得到的出口曲线最大绝对差 1e-9且整线 MRT数值积分Σ t·E·dt / Σ E·dt落在(θ₁2τ₁)(θ₂2τ₂)的±3%以内。五、小结与下一篇预告本篇把卷积这条扰动传播的唯一法则亲手跑通了第一最小可跑模型只有三个部件——进料器生成质量流量、混合器做 RTD 卷积、出口输出记录判据第二c_out c_in ∗ E因为串联叠加的是停留时间这一随机变量而不是浓度所以是卷积不是相加第三阶跃响应上迟滞来自 PFR 延迟 θ、平滑来自 CSTR 返混二者叠加成 S 形总 MRT θ 2τ_tank第四离散卷积必须×dt且Σ E·dt 1步长要满足dt ≤ θ/10且dt ≤ τ_tank/10时间轴要覆盖≥5×MRT第五长序列用scipy.signal.fftconvolve加速且np.convolve默认modefull会改变输出长度。术语上请守住两条写卷积convolution“不写褶积”凡涉及扰动传播一律走入口 ∗ RTD绝不允许直接乘一个衰减系数——那是物理错误。下一篇第 5 篇从函数升级到类给出仿照 PharmaPy 单元操作组织方式的自研示例类骨架PowderFeeder与ContinuousBlender讲清构造入参螺杆直径/螺距、料斗几何、转速、设定流量、混合器体积/转速、step()接口与结果对象以及如何把它们接入SimulationExec的 flowsheet 拓扑串形如F01 -- B01 -- TP01。第 6~7 篇将用 feed factor 半经验模型与料斗质量平衡把进料器做真第 9 篇详述 RTD 数学E(t)、F(t)、矩第 12 篇讲 PFR 2CSTR 组合模型第 16 篇讲三类扰动如何经卷积衰减、以及分流diversion的触发逻辑。本篇认知问题回显FAQQ1一个进料-混合串联系统的最小可跑模型需要哪三个部件A进料器生成入口质量流量/浓度序列、混合器用 RTD 做卷积变换、出口输出记录 c_out 并判定是否超规格接口顺序为生成 → 变换 → 输出后续自研类只是把前两个部件写得更真实。Q2卷积到底在算什么为什么串联等于卷积而不是相加A卷积是加权求和出口 t 时刻浓度等于入口各历史时刻浓度按物料在设备内停留 x 秒的概率 E(x) 加权后累加串联叠加的是停留时间这一随机变量两个独立随机变量之和的分布就是二者分布的卷积。Q3为什么阶跃扰动下出口会同时出现迟滞与平滑两个特征A迟滞来自设备中的活塞流段物料整体被推着先到先出由 PFR 延迟 θ 决定何时开始变化平滑来自混合段的返混颗粒停留时间呈分布由 τ_tank 决定 S 形爬升形状二者叠加即总 MRT θ 2τ_tank。Q4离散时间卷积和连续卷积差在哪采样步长怎么选A连续卷积是积分离散近似必须乘 dt 且要求 sum(E)*dt 1采样步长要显著小于被解析的最短时间尺度经验法则是 dt ≤ θ/10 且 dt ≤ τ_tank/10时间轴还需覆盖至少 5 倍 MRT 以免尾部截断。Q5为什么 ∫E(t)dt 1 是数值陷阱长序列卷积怎么加速ARTD 是概率密度积分不为 1 意味着物料不守恒常见原因是时间轴太短截断了长尾、归一化检查漏乘 dt长序列用 scipy.signal.fftconvolve 把复杂度从 O(N·M) 降到约 O(N log N)且 np.convolve 默认 mode‘full’ 会改变输出长度。 
长河君 · 2026/10/14 3:08:44 · 阅读 1.2k 连续粉体进料与混合建模教程（4）：第一个进料混合串联仿真——阶跃扰动下的最小可跑闭环
连续粉体进料与混合建模教程4第一个进料混合串联仿真——阶跃扰动下的最小可跑闭环版本声明块工具/软件numpy卷积、scipyscipy.signal.fftconvolve快速卷积、matplotlib阶跃响应图实际产线建模可再接入 PharmaPyPurdue CryPTSys作为 flowsheet 骨架。语言/环境Python 3.9numpy2、scipy1.14中文字体使用SimHei或Microsoft YaHei。本文全部数值为教学示意参数非真实产线实测数据。本文目标从零写出第一个进料 → 混合 → 出口浓度的最小串联仿真把卷积这条扰动传播的唯一法则亲手跑通一遍。一句话结论一个最小可跑的进料-混合串联模型只需要三个部件——进料器生成质量流量序列、混合器用np.convolve(feed_profile, rtd) * dt做 RTD 卷积、出口浓度输出在入口阶跃扰动下出口会呈现先迟滞由 PFR 延迟 θ 造成后平滑由返混造成的 S 形响应且只有当sum(E)*dt 1即∫E(t)dt 1时稳态出口浓度才等于入口浓度。〇、本篇要解决的认知问题前 3 篇都在讲清楚这一篇开始动手。五个必须先想明白的问题一个进料-混合串联系统的最小可跑模型最少需要哪三个部件各自负责什么卷积到底在算什么为什么串联 卷积而不是串联 相加为什么在阶跃扰动下出口会同时出现迟滞与平滑两个特征它们分别来自哪里离散时间卷积和连续卷积差在哪里采样步长dt该怎么选∫E(t)dt 1为什么是数值陷阱长序列卷积为什么要用scipy.signal.fftconvolve一、机制解析1.1 最小模型只需三个部件连续上游段的动态最核心的一条链是进料器把物料送出去混合器把扰动抹平出口浓度随之变化。把它写成最小模型┌───────────────┐ 质量流量序列 ┌───────────────┐ 出口浓度 ┌──────────┐ │ 进料器 │ ────────────────► │ 混合器 │ ──────────► │ 出口输出 │ │ feed_profile │ c_in(t) / F(t) │ RTD 卷积 │ c_out(t) │ 记录/判据 │ │ 生成质量流量 │ │ E(t) │ │ │ └───────────────┘ └───────────────┘ └──────────┘ ① ② ③部件职责输入输出① 进料器生成入口质量流量/浓度序列设定值、扰动feed_profile随时间变化② 混合器用 RTD 把入口序列变换为出口序列feed_profile、E(t)c_out(t)③ 出口输出记录与判定c_out(t)曲线 / 是否超规格为什么这对你重要把系统拆成生成 → 变换 → 输出三段是通用套路。进料器负责造扰动混合器负责变扰动输出负责判扰动。你后面写PowderFeeder、ContinuousBlender类时本质就是给①和②各写一个更真实的实现接口顺序不变。1.2 卷积串联系统的组合法则先给定义。连续卷积是c_out(t) ∫₀ᵗ c_in(t − x) · E(x) dx E 为 RTD量纲 1/时间两个 RTD 的串联组合也是卷积E₁ ∗ E₂ ∫₀ᵗ E₁(x) · E₂(t − x) dx这就是扰动沿产线传播全部计算的基石。它的直观含义是加权求和出口在t时刻的浓度是入口在过去所有时刻t−x的浓度各自乘上这段物料在设备里刚好待了x秒的概率E(x)再加起来。手算一下最粗糙的版本把时间切成 3 个格子入口 c_in [0, 1, 1] RTD 采样 E [0.5, 0.5, 0] t0: c_out 0*0.5 0 t1: c_out 0*0 1*0.5 0.5 t2: c_out 0*0 1*0 1*0.5 0.5 ← 注意每格还要乘 dt 才是积分为什么串联 卷积而不是串联 相加因为串联不是把浓度叠加而是把停留时间叠加。物料穿过第一台设备的耗时是一个随机变量穿过第二台又是一个随机变量两个独立的随机变量之和的分布正是二者分布的卷积。相加是两份物料混在一起卷积是同一份物料先后穿过两台设备——两者物理含义完全不同。1.3 迟滞与平滑同一个响应里的两个特征在出口曲线上你会同时看到两个特征它们来自两种不同的物理机制特征表现来源数学对应迟滞延迟出口比入口晚了一段时间才开始变化设备的活塞流段物料整体被推着走先到先出PFR 的延迟时间θ平滑返混出口不是阶跃而是缓慢的 S 形爬升设备的混合段颗粒停留时间有分布CSTR 的返混τ_tank一个常用的教学 RTD 把两者叠在一起PFR 延迟 两个等体积 CSTRE(t) (t − θ) · exp[ −(t − θ)/τ_tank ] / τ_tank² t θt ≤ θ 时为 0 总 MRT θ 2·τ_tank第一项(t − θ)保证了t θ之前出口完全为 0——这就是迟滞而指数项带来的拖尾就是平滑。入口阶跃 出口响应 │ ┌────── │ ,------ ← 平滑返混决定爬升形状 │ │ │ ,- │ │ │ ,- │──────┘ │,- ← 迟滞 θ 决定何时开始动 └──────────► t └──────────────► t t₀ t₀θ为什么这对你重要这两个特征对应两种完全不同的工程对策。迟滞主要由输送路径长度与流量决定改不了多少而平滑程度由混合器持料量与转速决定可调。分不清这两者你调参时就会按错旋钮。1.4 离散时间卷积 vs 连续卷积×dt不能漏计算机只能算离散序列。连续卷积∫ c_in(t−x)E(x)dx用离散近似是c_out[n] Σₖ c_in[n−k] · E[k] · dt dt 是采样步长三件事必须记住×dt是积分变离散的必备因子。漏掉它输出会整体放大1/dt倍——这是最隐蔽的量纲错误。E必须按dt采样。你构造的是E(t)在t 0, dt, 2dt, …上的取值数组而不是若干个概率值。采样步长dt的选择有硬约束dt必须显著小于被解析的最短时间尺度。经验法则是dt ≤ θ/10且dt ≤ τ_tank/10同时序列总长要覆盖 RTD 尾部一般取≥ 5×MRT否则尾部被截断ΣE·dt会小于 1。dt 太大 ──► 迟滞被抹平、爬升失真把 S 形算成折线 dt 太小 ──► 计算量暴涨且数值噪声被放大 dt 合适 ──► 能分辨 θ 与 τ_tank 的最小时间常数为什么这对你重要第 3 篇讲的时间尺度错配在这里第一次变成具体的dt选择。你现在选的dt直接决定后面所有卷积结果的可信度。1.5 归一化陷阱∫E(t)dt 1RTD 是一个概率密度必须满足∫₀^∞ E(t) dt 1数值上等价于Σ E[k] · dt 1。三个最常见的坑截断尾部时间轴太短RTD 的长尾被切掉积分 1出口稳态浓度会偏低漏乘dt归一化检查时忘了×dt看到Σ E远大于 1 就误以为自己构造错了离散化偏差用稀疏采样构造E(t)时梯形积分的误差在大dt下不可忽略。为什么要较真归一化不成立意味着物料不守恒。在连续制造里物料守恒是硬约束——你把 1 kg/h 的示踪剂打进产线出口累计就不可能超过 1 kg/h。归一化错误会直接让你在第 18 篇的物料追溯与分流计算里得出荒谬结论。1.6 用fftconvolve加速长序列直接卷积的复杂度是O(N·M)。当你为了解析毫秒级压片循环、或为了跑几百秒的产线动态而把序列拉到十万点以上时直接卷积会明显变慢。scipy.signal.fftconvolve用快速傅里叶变换把复杂度降到约O(N log N)序列越长加速比越大。一个反直觉的默认值必须提前说np.convolve的默认modefull输出长度是N M − 1比输入长出一截。如果你不做截断向下游传的序列会越来越长几级串联后时间轴就错位了。正确做法是modefull之后截断回输入长度或用modesame并核对对齐方式。二、完整代码与逐行剖析下面第一段脚本就是本系列的第一口活气端到端最小串联 阶跃响应。全部数值为教学示意参数。# -*- coding: utf-8 -*- 第 4 篇配套脚本 A进料-混合串联最小可跑闭环 阶跃响应。 依赖numpy, matplotlib。教学示意参数非真实产线数据。 importnumpyasnpimportmatplotlib matplotlib.use(Agg)importmatplotlib.pyplotasplt plt.rcParams[font.sans-serif][SimHei,Microsoft YaHei]plt.rcParams[axes.unicode_minus]Falsedefbuild_rtd(t,theta,tau_tank): 构造“PFR 延迟 theta 两个等体积 CSTR”的教学示意 RTD。 E(t) (t - theta) * exp(-(t - theta)/tau_tank) / tau_tank^2t theta否则 0。 总 MRT theta 2 * tau_tank。 Enp.zeros_like(t)maskttheta shiftt[mask]-theta E[mask]shift*np.exp(-shift/tau_tank)/tau_tank**2returnEdefsimulate_line(feed_profile,rtd,dt): 端到端最小串联出口 入口 * RTD离散卷积近似连续卷积须乘 dt。 feed_profile : 入口浓度/流量序列与 rtd 同一时间基准 rtd : E(t) 序列1/s按 dt 采样 dt : 采样步长s 返回与 feed_profile 等长的出口序列保证时间轴对齐。 fullnp.convolve(feed_profile,rtd,modefull)*dtreturnfull[:len(feed_profile)]# 截断回输入长度避免时间轴错位defmain():dt1.0# 采样步长 1 stnp.arange(0.0,400.0,dt)theta,tau_tank20.0,15.0# 教学示意延迟 20 s单釜 15 srtdbuild_rtd(t,theta,tau_tank)areartd.sum()*dtprint(fRTD 归一化检查: sum(E)*dt {area:.6f}应 ≈ 1)ifabs(area-1.0)1e-3:print(警告归一化偏差偏大请延长时间轴或减小 dt。)c_in(t100.0).astype(float)# 入口阶跃t100 s 起由 0 变 1c_outsimulate_line(c_in,rtd,dt)idx_halfint(np.argmax(c_out0.5))# 首次达到终值 50% 的时刻print(f入口阶跃时刻:{100.0:.0f}s)print(f出口达到终值 50% 的时刻:{t[idx_half]:.0f}s迟滞 爬升的合成)print(f出口稳态值末点:{c_out[-1]:.4f}应 ≈ 1前提是 RTD 归一化)fig,axplt.subplots(figsize(9,4))ax.plot(t,c_in,--,label入口浓度 c_in阶跃)ax.plot(t,c_out,-,linewidth2,label出口浓度 c_outS 形响应)ax.axvline(100.0,colorgray,linestyle:,linewidth1)ax.set_xlabel(时间 ts)ax.set_ylabel(归一化浓度-)ax.set_title(进料-混合串联的阶跃响应迟滞 平滑教学示意)ax.legend()fig.tight_layout()fig.savefig(step_response.png,dpi150)print(已生成 step_response.png)if__name____main__:main()逐段剖析build_rtd的第一行E np.zeros_like(t)是必须的。先生成全零数组再只给t θ的部分赋值自然保证了t ≤ θ时E 0——这就是迟滞的数值实现。若用向量化公式直接算t θ时(t − θ)为负、会产生负的概率物理上荒谬。simulate_line里的* dt是全篇最关键的一个乘号。漏掉它出口浓度会放大1/dt 1倍这里 dt1 恰好掩盖了错误——所以你必须用dt 0.5或2.0各跑一次来验证。这是本篇最反直觉的地方dt 1时漏乘dt和不漏乘结果一样会让你误以为代码是对的。full[: len(feed_profile)]的截断。np.convolve(modefull)输出长度是N M − 1多出M − 1个点。截断回N保证c_in和c_out在同一时间轴上一一对应。c_in (t 100.0).astype(float)用布尔数组生成阶跃。比循环赋值快得多而且不会漏点.astype(float)把True/False转成1.0/0.0。idx_half用np.argmax(c_out 0.5)。返回第一个为真的下标用来量化迟滞 爬升合成后的半程时刻。教学参数θ 20 s、τ_tank 15 s、阶跃在t 100 s下它约等于145 s 100 θ 1.68·τ_tank必然大于100 θ 120 s——出口还要花约 25 s 爬升这个数字比单纯的θ更能说明迟滞与平滑是叠加的。第二段脚本演示归一化陷阱与fftconvolve加速并顺手给出第 10 篇会用到的 n-CSTR 形状 RTD。# -*- coding: utf-8 -*- 第 4 篇配套脚本 B归一化陷阱与 fftconvolve 加速。 依赖numpy, scipy。教学示意参数。 importtimefrommathimportfactorialimportnumpyasnpfromscipy.signalimportfftconvolvedefrtd_gamma(t,tau,n): n 个等体积 CSTR 串联的无量纲形状 RTD教学示意深入讨论见第 10 篇。 E(t) n^n / (n-1)! * (t/tau)^(n-1) * exp(-n * t/tau) / tau Enp.zeros_like(t)maskt0xt[mask]/tau E[mask](n**n)/factorial(n-1)*x**(n-1)*np.exp(-n*x)/taureturnEdefnormalize_rtd(E,dt):把离散 RTD 归一化使 sum(E)*dt 1。返回 (归一化后序列, 原始积分)。areaE.sum()*dtreturnE/area,areadefmain():dt0.5tnp.arange(0.0,300.0,dt)Ertd_gamma(t,tau20.0,n3)E_norm,areanormalize_rtd(E,dt)print(f归一化前 sum(E)*dt {area:.6f})print(f归一化后 sum(E)*dt {E_norm.sum()*dt:.6f}应 ≈ 1)# 陷阱演示np.convolve 默认 modefull输出长度 N M - 1N,M200000,2000rngnp.random.default_rng(0)arng.random(N)brng.random(M)print(f\n输入长度 N{N}RTD 长度 M{M})t0time.time()conv_directnp.convolve(a,b,modefull)t_directtime.time()-t0 t0time.time()conv_fftfftconvolve(a,b,modefull)t_ffttime.time()-t0print(fnp.convolve 耗时{t_direct*1000:8.1f}ms输出长度{len(conv_direct)})print(ffftconvolve 耗时{t_fft*1000:8.1f}ms输出长度{len(conv_fft)})print(f两者最大绝对差 {np.max(np.abs(conv_direct-conv_fft)):.3e})print(f加速比 ≈{t_direct/max(t_fft,1e-9):.1f}x序列越长越明显)if__name____main__:main()逐段剖析rtd_gamma的分母是factorial(n - 1)而不是factorial(n)。n-CSTR 的标准形状是E(θ) n^n/(n−1)! · θ^(n−1) · exp(−nθ)再除以 τ 换成有量纲的E(t)。写错分母会让积分不等于 1——正好被第二段代码的归一化检查抓出来。normalize_rtd返回原始积分是为了让你看到截断或采样不足到底差了多少。工程上习惯先归一化再使用但先看原始积分是纪律如果原始积分只有 0.7那说明你的时间轴或dt根本不够归一化只是掩盖了问题。np.random.default_rng(0)用固定种子保证你我的加速比测试跑在同样的输入上结果可复现。这是性能对比脚本的通用写法。fftconvolve与np.convolve结果应几乎相同差异在1e-12量级浮点误差。如果你看到显著差异几乎一定是边界处理或长度不匹配的问题。time.time()做基准测试是粗粒度手段适合这种差一个数量级的对比要精确测时请用time.perf_counter。三、常见报错与排查1现象出口稳态浓度不等于入口浓度比如入口阶跃到 1出口只到 0.83。根因归一化不成立——Σ E · dt ≠ 1通常是时间轴太短、RTD 尾部被截断。解法打印rtd.sum() * dt把时间轴延长到≥ 5×MRT或减小dt。归一化前先看原始积分别直接用归一化掩盖截断问题。2现象出口曲线整体比预期高一到几个数量级。根因离散卷积漏乘dt。连续卷积∫c_in·E·dx近似为离散和时必须×dt。解法在np.convolve(...)后乘dt并用dt 0.5与dt 2.0各跑一次验证结果的dt无关性只用dt 1是验证不出这个错误的。3现象多级串联后出口序列越来越长和输入的t轴对不上。根因np.convolve默认modefull输出长度是N M − 1。解法卷积后截断回输入长度full[:len(x)]或改用modesame并统一对齐规则每级串联都把长度钉死在N。4现象阶跃响应曲线看起来是折线没有明显的 S 形平滑段。根因采样步长dt太大把返混造成的爬升抹平成了直线段。解法按dt ≤ θ/10且dt ≤ τ_tank/10收紧步长教学示例里θ 20 s、τ_tank 15 s故dt应取1 s或更小。5现象长序列卷积跑得很慢几百秒的产线动态要等很久。根因直接用O(N·M)的np.convolve序列上十万点时耗时明显。解法改用scipy.signal.fftconvolve约O(N log N)序列越长加速比越大。注意两种实现的边界与长度保持一致。四、动手练习练习 1基础跑通最小闭环并核对三个数。把脚本 A 存为first_lineup.py运行。成功标准量化打印的sum(E)*dt落在1.000 ± 0.005出口稳态末值落在0.98–1.02出口达到终值 50% 的时刻落在145 ± 5 s。练习 2进阶验证×dt的必要性本练习是本篇核心。把dt改成0.5与2.0各跑一次先故意删掉simulate_line里的* dt观察结果。成功标准量化删掉* dt后dt0.5时出口稳态约为2.0、dt2.0时约为0.5即偏离真值1/dt倍把* dt加回来后两种步长下出口稳态都落在0.98–1.02。练习 3挑战给串联加一段并核对整线 MRT。把simulate_line改成两段串联先用build_rtd生成两段 RTD参数可自定如θ₁5, τ₁8与θ₂10, τ₂12先分别卷积、再把两段 RTD 卷积后再卷积入口。成功标准量化两种算法得到的出口曲线最大绝对差 1e-9且整线 MRT数值积分Σ t·E·dt / Σ E·dt落在(θ₁2τ₁)(θ₂2τ₂)的±3%以内。五、小结与下一篇预告本篇把卷积这条扰动传播的唯一法则亲手跑通了第一最小可跑模型只有三个部件——进料器生成质量流量、混合器做 RTD 卷积、出口输出记录判据第二c_out c_in ∗ E因为串联叠加的是停留时间这一随机变量而不是浓度所以是卷积不是相加第三阶跃响应上迟滞来自 PFR 延迟 θ、平滑来自 CSTR 返混二者叠加成 S 形总 MRT θ 2τ_tank第四离散卷积必须×dt且Σ E·dt 1步长要满足dt ≤ θ/10且dt ≤ τ_tank/10时间轴要覆盖≥5×MRT第五长序列用scipy.signal.fftconvolve加速且np.convolve默认modefull会改变输出长度。术语上请守住两条写卷积convolution“不写褶积”凡涉及扰动传播一律走入口 ∗ RTD绝不允许直接乘一个衰减系数——那是物理错误。下一篇第 5 篇从函数升级到类给出仿照 PharmaPy 单元操作组织方式的自研示例类骨架PowderFeeder与ContinuousBlender讲清构造入参螺杆直径/螺距、料斗几何、转速、设定流量、混合器体积/转速、step()接口与结果对象以及如何把它们接入SimulationExec的 flowsheet 拓扑串形如F01 -- B01 -- TP01。第 6~7 篇将用 feed factor 半经验模型与料斗质量平衡把进料器做真第 9 篇详述 RTD 数学E(t)、F(t)、矩第 12 篇讲 PFR 2CSTR 组合模型第 16 篇讲三类扰动如何经卷积衰减、以及分流diversion的触发逻辑。本篇认知问题回显FAQQ1一个进料-混合串联系统的最小可跑模型需要哪三个部件A进料器生成入口质量流量/浓度序列、混合器用 RTD 做卷积变换、出口输出记录 c_out 并判定是否超规格接口顺序为生成 → 变换 → 输出后续自研类只是把前两个部件写得更真实。Q2卷积到底在算什么为什么串联等于卷积而不是相加A卷积是加权求和出口 t 时刻浓度等于入口各历史时刻浓度按物料在设备内停留 x 秒的概率 E(x) 加权后累加串联叠加的是停留时间这一随机变量两个独立随机变量之和的分布就是二者分布的卷积。Q3为什么阶跃扰动下出口会同时出现迟滞与平滑两个特征A迟滞来自设备中的活塞流段物料整体被推着先到先出由 PFR 延迟 θ 决定何时开始变化平滑来自混合段的返混颗粒停留时间呈分布由 τ_tank 决定 S 形爬升形状二者叠加即总 MRT θ 2τ_tank。Q4离散时间卷积和连续卷积差在哪采样步长怎么选A连续卷积是积分离散近似必须乘 dt 且要求 sum(E)*dt 1采样步长要显著小于被解析的最短时间尺度经验法则是 dt ≤ θ/10 且 dt ≤ τ_tank/10时间轴还需覆盖至少 5 倍 MRT 以免尾部截断。Q5为什么 ∫E(t)dt 1 是数值陷阱长序列卷积怎么加速ARTD 是概率密度积分不为 1 意味着物料不守恒常见原因是时间轴太短截断了长尾、归一化检查漏乘 dt长序列用 scipy.signal.fftconvolve 把复杂度从 O(N·M) 降到约 O(N log N)且 np.convolve 默认 mode‘full’ 会改变输出长度。