pyTEP：TEP故障诊断仿真数据生成与二次开发

简介pyTEP 是面向 Tennessee Eastman田纳西伊士曼化工过程的 Python 开源仿真 API适合过程控制、故障诊断与优化算法的学习者与研究者。它把复杂的仿真场景搭建封装为高层接口用户无需理解模拟器底层原理即可设置仿真、调整操作条件并存储仿真数据同时保留交互式仿真的灵活性。压缩包共 76 个文件、约 754KB以 23 个 .m 和 7 个 .mdl 的 MATLAB/Simulink 模型为仿真引擎16 个 .py 脚本承担 Python 侧的接口与工具函数另含 .pkl、.mat 数据文件与 .rst 文档、测试用例形成引擎、接口、文档齐全的完整结构。读者可据此了解 pyTEP 的工程组织方式包括 simulator、utils、matlab_bridge 等模块的职责划分与仿真实例的调用路径便于在本地复现并扩展自己的实验。目前已有 546 人学习下载。1. 用 Python 造假数据不难难的是造出物理上说得通的数据做过程监控和故障诊断的人应该都遇到过这个场面算法写完了评估指标也定义好了回头一看——没有数据。真实产线的传感器数据不可能落到个人手上公开的那几份表格数据集又固定死了想验证自己调的检测器在“进料组成缓慢漂移”这种场景下到底灵不灵根本没有对应的样本。pyTEP 就是冲这个缺口来的它把 TEP 这个流程控制领域用了三十年的标准仿真基准直接封装成一个 Python 包tep.TEP()一行实例化simulate()一行出数41 个测量变量、12 个操纵变量的闭环时序数据就到手了还能顺手把 21 种预设故障打开。适合谁做故障诊断、过程监控、控制算法验证的开发者和研究生以及需要离线批量造带标签异常数据的人。它不做界面不画图就是一个干净的 API剩下的活自己接。2. 架构拆解从 TEP 基准到 pyTEP 的三层映射2.1 三层结构数学模型、状态推进、API 封装理解 pyTEP 的关键是把它看成三层而不是一个“仿真函数”。最底下是 TEP 基准本身的数学模型一个带再循环的放热反应流程包含反应器、冷凝器、汽提塔、压缩机和气液分离器。这套模型在上世纪九十年代初被公开之后成了过程控制界的公共考卷——几十个微分方程加代数方程描述物料平衡、能量平衡和相平衡。这一层决定了“物理上说得通”也是为什么它造出来的异常曲线比手工拼的随机数可信。中间一层是状态推进内核。每一步要干三件事把当前状态按微分方程积分一小步用控制回路把 12 个操纵变量算出来并施加回流程再把测量变量加噪后输出。这一层是纯数值计算不涉及任何 Python 对象设计。最上面一层才是 pyTEP一个薄薄的 Python 封装把前两层包成“构造对象 → 调 simulate → 拿 numpy 数组”的接口。它只做三件事参数收口、状态保存、噪声注入。所以整个包不大你不需要啃几千行代码就能改。提示正因为是三层你改东西之前先想清楚自己要改哪一层。只想换噪声水平动第一层的参数就够了想换控制策略得动中间层想改接口风格才碰最上面这层。2.2 目录结构与模块职责这类科学仿真包的结构通常很规整我拿到手习惯先看目录基本能判断改动的成本落在哪里。常见做法是这样组织的pytep/ ├── tep/ │ ├── __init__.py # 导出 TEP 类决定 import tep 之后能拿到什么 │ ├── tep.py # 对外 API 层TEP 类的构造、simulate、状态保存 │ ├── plant.py # 流程模型层微分方程、状态向量、积分推进 │ └── controller.py # 控制回路层12 个操纵变量的闭环计算 ├── examples/ # 跑通示例通常就是 minimal 那几行 ├── tests/ # 回归测试验证给定故障下的输出是否漂移 ├── setup.py # 打包脚本pip install -e . 走这里 └── README.md # API 说明和示例tep.py是入口plant.py是心脏。二次开发里 90% 的改动要么发生在tep.py加接口、加字段、加包装要么发生在controller.py换控制器真正动plant.py的很少——改了流程模型就不叫 TEP 了别人也没法拿你的结果跟公开基线对比。这个分层还有个现实好处plant.py和controller.py通常不依赖框架纯 numpy 加 Python 数学你可以单独把它们抽出来放进自己的项目里不装整个包也能跑。我在做嵌入式侧的离线验证时就这么干过只搬了状态推进那两个文件省掉一堆依赖。2.3 41 和 12每个数字都对应一个物理量simulate()返回两张表一张 41 列一张 12 列。这个列数不能靠猜得知道每一列是什么否则后面做特征工程全错。XMEAS 系列按功能大致分成四组组别编号范围内容数量流量类1-6, 10, 14, 17各路进料、循环、反应器进料、排放、分离器与汽提塔底部流量9温度与压力7-9, 11, 13, 16, 18, 21, 22反应器温度压力液位、分离器温度压力、汽提塔温度压力、两路冷却水出口温度9阀门开度19, 20汽提塔蒸汽阀、压缩机循环阀2液位8, 12, 15反应器、分离器、汽提塔液位3成分分析23-41各路物料中 A-H 组分的摩尔分数19XMV 系列是 12 个操纵变量也就是你能“拧”的旋钮编号变量编号变量XMV(1)D 进料流量XMV(7)分离器液体产品阀XMV(2)E 进料流量XMV(8)汽提塔产品阀XMV(3)A 进料流量XMV(9)汽提塔蒸汽阀XMV(4)A 和 C 进料流量XMV(10)反应器冷却水阀XMV(5)循环阀XMV(11)冷凝器冷却水阀XMV(6)排放阀XMV(12)搅拌速度有个细节值得注意XMV(12) 搅拌速度在标准 TEP 里其实是个近似常数公开数据集中常常不带这一列。所以如果你拿别人整理的 52 列数据集41 11来跟 pyTEP 的输出对拍会发现多了一列这不是 bug是取舍不同。这一点后面避坑章节还会再说。3. 跑通流程安装、首跑与 21 种故障数据集的批量生成3.1 依赖与安装依赖很轻正常情况就是 numpy绘图可选 matplotlib。装法两种用 pip 从包索引装或者把源码包下下来走开发模式装——我强烈建议后者因为你在二次开发时几乎一定要改源码装成 site-packages 里的正式包会让你改得很难受。# 方式一从包索引安装快速试水 pip install tep # 方式二源码包内开发模式安装准备改代码时用这个 cd pytep python -m venv .venv source .venv/bin/activate # Windows 用 .venv\Scripts\activate pip install -e . # -e 表示可编辑安装改源码立即生效 pip install numpy matplotlib # 补上常用依赖-e这个参数很关键。它把包目录软链到环境的 site-packages 里你改tep/tep.py保存后下一次 import 直接生效不用重装。血泪经验第一次接触这类包的时候我用的是普通安装改了半天代码发现行为没变排查了快一个小时才发现改的是源码目录、跑的是装好的副本。验证装好了没有别急着跑仿真先做个最小 importimport tep print(tep.__file__) # 确认加载路径指向你改的那个目录 print(dir(tep)) # 看看顶层导出了哪些名字第二行尤其有用。如果dir(tep)里只有TEP说明包的出口很窄所有能力都收在一个类上如果还能看到plant、controller之类的子模块说明作者把内部结构也暴露了你可以直接调用里面的函数做更细粒度的实验。3.2 十行代码跑出第一条曲线别一上来就想着造数据集先跑通一次最小仿真确认环境没问题。最小可运行代码如下import numpy as np import tep # fault0 表示无故障运行xmeas_noise 是测量噪声强度 model tep.TEP(fault0, xmeas_noise0.1) # n_samples 是采样点个数timestep 单位是小时0.05 小时 3 分钟 xmeas, xmv model.simulate(n_samples500, timestep0.05) print(测量变量形状:, xmeas.shape) # 期望 (500, 41) print(操纵变量形状:, xmv.shape) # 期望 (500, 12) print(反应器温度前 5 个点:, xmeas[:5, 8]) print(D 进料前 5 个点:, xmv[:5, 0])这段代码里有三个参数决定了你拿到的数据长什么样fault是故障编号0 表示正常运行1 到 21 对应 21 种预设异常进料阶跃、冷却水温度漂移、阀门粘滞等。xmeas_noise是测量噪声幅值TEP 基准在正常情况下本身就是有噪声的这个值给 0 会得到一条光滑得不像真的曲线拿去做检测算法会严重高估性能。timestep单位是小时这是最容易翻车的地方——很多人以为单位是秒写了个timestep1想表示 1 秒实际含义是一步走 1 小时等于把时间轴压缩了 20 倍。500 个采样点、每点 3 分钟对应 25 小时的流程时间。标准 TEP 实验一般跑 48 小时也就是 960 个点你可以按自己的需求缩放。跑完先别急着保存画一眼看看形状对不对import matplotlib.pyplot as plt fig, axes plt.subplots(2, 1, figsize(10, 6), sharexTrue) axes[0].plot(xmeas[:, 8]) # XMEAS(9) 反应器温度 axes[0].set_ylabel(Reactor Temp) axes[1].plot(xmv[:, 9]) # XMV(10) 反应器冷却水阀 axes[1].set_ylabel(CW Valve) axes[1].set_xlabel(Sample) plt.tight_layout() plt.savefig(first_run.png, dpi120)正常情况下你会看到反应器温度围绕某个稳态小幅波动冷却水阀开度也在稳态附近微调——这就是那个闭环控制回路在工作。如果温度一路单调漂走或者剧烈震荡八成是构造参数给得太离谱比如噪声给了个 10 这种量级。3.3 21 种故障的批量生成单次仿真只是热身pyTEP 真正的价值在于批量造带标签的数据集。下面这个脚本是我常用的一版直接抄改成自己的路径就能用import os import numpy as np import tep OUT_DIR tep_dataset N_SAMPLES 960 # 960 * 3 分钟 48 小时对齐标准实验时长 TIMESTEP 0.05 # 单位小时 SEED_BASE 20240501 os.makedirs(OUT_DIR, exist_okTrue) # fault0 是正常工况1..21 是 21 种预设故障 for fault in range(0, 22): # 每个故障单独构造一个模型实例避免状态从上一次仿真里带过来 model tep.TEP(faultfault, xmeas_noise0.1, xmv_noise0.0) xmeas, xmv model.simulate(n_samplesN_SAMPLES, timestepTIMESTEP) # 标签单独存一份方便后面拼成分类任务的 y label 0 if fault 0 else 1 np.savez_compressed( os.path.join(OUT_DIR, ffault_{fault:02d}.npz), xmeasxmeas.astype(np.float32), xmvxmv.astype(np.float32), faultfault, labellabel, ) print(ffault{fault:02d} done, xmeas{xmeas.shape}, xmv{xmv.shape})几个决策点说明一下。astype(np.float32)是为了省空间——41 列 × 960 行 × 22 个文件float64 存下来是几十兆float32 直接减半做机器学习训练时精度完全够。np.savez_compressed而不是 CSV因为列多行多CSV 读写慢且容易在浮点精度上出现往返误差。range(0, 22)这个区间要特别注意0 到 21 一共 22 个取值不是 21 个。很多人写range(21)结果是跑了 0 到 20漏掉最后一个故障还顺手把无故障工况也包含进去了标签对不上的时候极难发现。如果要做训练集/测试集划分常规做法是每个故障跑多个随机种子前几个当训练、后几个当测试这样测试集里是同一种故障但不同噪声实现能测出算法对噪声的鲁棒性而不只是记住了某条曲线。3.4 数据落盘格式与目录组织文件组织方式决定了你后面做实验的效率。22 个文件平铺在一层目录里写到第三个算法脚本时就会开始难受。我一般按工况分两级tep_dataset/ ├── normal/ │ ├── seed_00.npz │ └── seed_01.npz ├── fault_01/ │ ├── seed_00.npz │ └── seed_01.npz ├── ... └── fault_21/好处是加载时用glob一把扫目录名本身就是标签不需要再读文件里的字段做映射。写一个 dataset 加载函数import glob import os import numpy as np def load_dataset(root): 扫描两级目录返回 (样本列表, 标签列表) X, y [], [] for case_dir in sorted(glob.glob(os.path.join(root, *))): if not os.path.isdir(case_dir): continue name os.path.basename(case_dir) label 0 if name normal else 1 # 二分类正常 / 异常 for fp in sorted(glob.glob(os.path.join(case_dir, *.npz))): data np.load(fp) X.append(data[xmeas]) y.append(label) return X, np.array(y)注意这里只取了xmeas。要不要把xmv也拼进去是个策略问题操纵变量反映的是控制器的动作对诊断“控制阀粘滞”这类故障是有信息量的但对“进料组成漂移”这种故障反而可能引入噪声。我的习惯是先只用xmeas建基线等基线稳定了再把xmv加进来做对照实验看提升是否显著。4. 核心链路simulate() 一次调用内部的采样、控制与故障注入4.1 对象初始化一次构造决定了后面的一切tep.TEP()被调用的那一刻做了三件事把初始稳态装进状态向量、把故障编号记下来、把噪声参数存好。这决定了后面所有采样的起点。# 构造时发生的事按常见实现的顺序大致是 # 1. 载入流程的初始稳态各容器液位、温度、压力、组分浓度 # 2. 记录 fault 编号决定后续在积分时叠加哪一路扰动 # 3. 保存噪声参数供输出阶段加噪使用 # 4. 初始化控制器的积分项避免第一步就出现巨大偏差 model tep.TEP(fault3, xmeas_noise0.1, xmv_noise0.0)第 4 步容易被忽略但很关键。控制器里通常有积分环节如果积分初值给 0 而实际偏差非 0前几十个采样点会出现一段明显的“爬升”暂态。这段暂态在真实数据里是存在的开车阶段但如果你要做稳态故障检测它就是个干扰项。这也是为什么标准实验流程要求先跑一段无故障预热再引入故障。对象一旦构造完内部就保存了完整的状态向量。simulate()不重新初始化它从当前状态往下接着算——这一点后面避坑章节会专门讲。4.2 单步推进一个采样周期内的时序simulate(n_samples500)内部是一个 500 次的循环每次循环对应一个采样周期。理解这个循环里各步骤的先后顺序是排查“数据对不上”的第一把钥匙。按常见实现一个周期大致长这样def _step(self, timestep): # 1) 采样读取当前时刻的测量值叠加测量噪声 xmeas_raw self._read_measurements() # 2) 控制用带噪声的测量值算控制量这才是真实的闭环 xmv self._controller.compute(xmeas_raw, self._setpoint) # 3) 扰动按故障编号注入对应异常可能是阶跃、漂移或随机游走 xmv_disturbed, state_disturbance self._apply_fault(xmv) # 4) 积分把控制量和扰动一起喂给流程模型推进一个时间步 self._state self._plant.integrate(self._state, xmv_disturbed, timestep) # 5) 输出记录本周期结果xmv 叠加执行机构噪声 return xmeas_raw, xmv self._xmv_noise * np.random.randn(len(xmv))顺序上有个细节值得注意控制器用的是加了噪声的测量值而不是真值。这符合真实工业场景——控制器只能看到传感器读数。如果你在二次开发时把这里改成用真值闭环性能会明显变好但数据就不真实了拿去做算法验证会高估效果。第 4 步的积分是计算量的大头。timestep越小数值积分越准但越慢0.05 小时是个折中正好等于标准 TEP 的 3 分钟采样间隔。不建议往更小调除非你在改控制器的响应速度。4.3 返回值的形状、顺序与内存布局返回值就是两个二维数组行是时间列是变量。没有列名没有表头位置即含义——这正是需要把 2.3 节那张表存在手边的原因。xmeas, xmv model.simulate(n_samples1000, timestep0.05) # 取单个变量的整条时间序列 reactor_temp xmeas[:, 8] # XMEAS(9)索引 8 对应第 9 个变量 cw_valve xmv[:, 9] # XMV(10)索引 9 对应第 10 个变量 # 取某个时间窗的切片第 20 到第 200 个采样点 window xmeas[20:200, :]这里最常见的翻车是索引偏一位。XMEAS 编号从 1 开始Python 索引从 0 开始所以xmeas[:, 8]才是 XMEAS(9)。我给自己的规矩是代码里凡是用到列索引的地方都写一行注释标明对应的物理量名三个月后回头看的时候能省下大量时间。另一个容易忽略的点是返回的数组是不是视图。如果simulate()内部用的是预分配缓冲区、逐步填充、最后返回切片那多次调用可能共享同一块内存你保存下来的数组会被后一次调用覆盖。稳妥做法是在拿到数据后立即.copy()或者立刻落盘。4.4 随机性与复现把种子管起来pyTEP 输出里至少有三种随机源测量噪声、执行机构噪声、故障本身的随机分量比如第 8、12 号这类的随机波动型故障。这意味着同一个模型跑两次结果不会完全相同。这在做算法对比时是个麻烦事。你改了检测器的阈值指标涨了两个点到底是改进有效还是这次运气好解决方法是固定种子import numpy as np def run_case(fault, seed, n_samples960, timestep0.05): np.random.seed(seed) # 固定全局随机源 model tep.TEP(faultfault, xmeas_noise0.1) return model.simulate(n_samplesn_samples, timesteptimestep)如果包内部用的是np.random.randn这类全局函数np.random.seed就够了如果它自己持有一个RandomState或者Generator实例那就需要在构造时把种子传进去。两种写法都试一遍跑两次看输出是否逐位相等就能判断出来。验证方法很简单同一段代码跑两次np.array_equal(a, b)返回 True 才算真复现。稳态工况下曲线看着差不多不代表数值相同——噪声幅值 0.1 的时候肉眼看不出差别但对检测算法来说那点差异可能正好跨过阈值。所以别靠眼睛判断复现性。5. 避坑与排查数据看着不对劲时先查这五处5.1 现象曲线变化快得离谱故障几秒就“发作”完了原因时间单位搞错了。timestep的单位是小时不是秒也不是分钟。写timestep1想表达 1 秒实际一步推进了 1 小时等于标准采样间隔的 20 倍。积分步长变大后流程的动态响应被压缩到少数几个采样点里看起来就是“瞬间爆发”。解决默认用0.05也就是 3 分钟跟标准 TEP 实验对齐。要算总时长就用n_samples * timestep小时960 个点对应 48 小时。改步长之前先问自己一句我要的是更细的时间分辨率还是更长的仿真时长这两个诉求不能靠调一个参数同时满足。5.2 现象明明设了fault5前 20 个采样点跟无故障跑出来一模一样原因故障不是从第 0 个采样点开始的。标准 TEP 实验规程是先在正常工况下跑 1 小时20 个采样点建立稳态然后再引入故障。pyTEP 通常遵循这个约定所以前 20 个点是干净的正常数据。解决造数据集时保留这段预热区训练时可以选择性丢弃前 20 个点评估检测延迟时把故障起点定义在第 20 个采样点而不是第 0 个。如果你只跑 20 个点就结束那整个样本里根本没有故障——这个坑我踩过排查了半天以为是故障注入失效其实就是样本太短。5.3 现象连续调两次simulate()第二次的数据开头跟第一次结尾接不上或者数据分布明显偏移原因simulate()是增量式的它从对象当前保存的状态接着往下算不做重置。第二次调用的起点是第一次的终点而不是初始稳态。如果你的控制器在这期间累积了积分偏差第二次的起点可能已经偏离稳态很远了。解决需要独立样本时每个样本重新tep.TEP(...)构造一个新对象别复用。这个成本不高构造本身不怎么耗时真正贵的是积分循环。要写循环造数据的话把model tep.TEP(...)放进循环体内这一行别省。5.4 现象拿 pyTEP 的输出跟公开的 TEP 数据集对拍列数对不上4112 和 4111 差了 1 列原因标准 TEP 的操纵变量有 12 个但其中 XMV(12) 搅拌速度在很多公开数据集里被视为常数、不单独列出所以流传的表格数据集常常是 52 列41 个 XMEAS 11 个 XMV。pyTEP 输出完整的 12 列自然对不上。解决对拍之前先对齐列定义把xmv[:, 11]搅拌速度去掉再比。另外注意 XMV 的顺序在不同实现里可能不完全一致别只核列数核一下具体变量的数值范围——比如进料流量和阀门开度的量级差好几个数量级错位了立刻能看出来。5.5 现象把xmeas直接喂给分类器准确率 99% 以上一到真实数据就崩原因量纲不统一。41 个变量里流量类可能在几十到几百的量级成分分析类XMEAS 23-41是摩尔分数量级在 0.01 到 0.5 之间。不做标准化直接送进基于距离的模型KNN、SVM、PCA 重构误差结果完全被大数值变量主导成分漂移这类故障根本检测不出来。999% 的准确率是在拟合噪声。解决按变量维度做 z-score 标准化均值方差只用训练集的正常样本算别用全体数据——用全体算等于把测试集信息漏进了预处理。另外成分分析变量本身有较强的自相关做差分或者滑窗统计时要注意直接对原始值求一阶差分会被噪声淹没。6. 二次开发把 pyTEP 接进算法流水线的三种改法改 pyTEP 之前先明确你要改哪一层。大多数需求其实不需要动流程模型本身包一层就够了。第一种改法包成 Gym 风格的交互环境。如果你要做强化学习或者在线自适应控制需要的是reset()和step(action)接口而不是一次性吐出 960 个点。这个改动量很小核心就三件事reset里重建 TEP 对象并返回初始观测step里推进一个采样周期并返回(obs, reward, done, info)done按最大步数判断。难点在 reward 设计而不是代码通常用“偏离设定点的加权平方和”取负值权重按变量量纲归一化别让温度项压过一切。第二种改法替换控制器。controller.compute()是控制回路的入口把里面的 PID 换成你的 MPC 或者数据驱动策略就能对比同一个流程在不同控制律下的表现。改的时候务必保留“控制器只能看到加噪测量值”这个约束去掉它你得到的性能提升是假的。验证方法是先用原控制器跑一遍记下反应器温度的标准差作为基线换控制器后再跑一遍只有明显低于基线才算真改进。第三种改法加自定义故障。21 种预设故障覆盖不了所有场景比如“传感器缓慢偏置”就没有标准编号。加自定义故障的常规位置是_apply_fault()在里面按fault编号返回不同的扰动。这里有个约定要坚持新故障用 100 以上的编号别占用 1-21否则后面有人拿你的代码跟标准数据集对比时会疯掉。最后说一个我自己的习惯。每次改完 pyTEP 相关代码我都会跑一个固定的回归脚本对故障 1、4、8 各跑 200 个采样点检查反应器温度的均值和方差是否落在预设区间内。这个脚本不做任何算法评估只回答一个问题——流程还是原来那个流程吗。有一次我把积分步长调小了想提升精度结果忘了同步调整故障注入的时序回归脚本立刻报警省了我至少半天时间。从那以后我每次动plant.py和controller.py之前都强制先跑一遍这个基线快照把改动前后的数字并排打出来看。这个习惯看着笨但它是我在这个项目上唯一没后悔过的规矩。希望帮到你。本文还有配套的精品资源点击获取