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 系列按功能大致分成四组：

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 直接生效，不用重装。血泪经验：第一次接触这类包的时候我用的是普通安装，改了半天代码发现行为没变，排查了快一个小时才发现改的是源码目录、跑的是装好的副本。

验证装好了没有，别急着跑仿真，先做个最小 import：

import tep print(tep.__file__) # 确认加载路径指向你改的那个目录 print(dir(tep)) # 看看顶层导出了哪些名字

第二行尤其有用。如果dir(tep)里只有TEP，说明包的出口很窄，所有能力都收在一个类上；如果还能看到plant、controller之类的子模块，说明作者把内部结构也暴露了，你可以直接调用里面的函数做更细粒度的实验。

3.2 十行代码跑出第一条曲线

别一上来就想着造数据集，先跑通一次最小仿真，确认环境没问题。最小可运行代码如下：

import numpy as np import tep # fault=0 表示无故障运行；xmeas_noise 是测量噪声强度 model = tep.TEP(fault=0, xmeas_noise=0.1) # n_samples 是采样点个数；timestep 单位是小时，0.05 小时 = 3 分钟 xmeas, xmv = model.simulate(n_samples=500, timestep=0.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单位是小时，这是最容易翻车的地方——很多人以为单位是秒，写了个timestep=1想表示 1 秒，实际含义是一步走 1 小时，等于把时间轴压缩了 20 倍。

500 个采样点、每点 3 分钟，对应 25 小时的流程时间。标准 TEP 实验一般跑 48 小时，也就是 960 个点，你可以按自己的需求缩放。

跑完先别急着保存，画一眼看看形状对不对：

import matplotlib.pyplot as plt fig, axes = plt.subplots(2, 1, figsize=(10, 6), sharex=True) 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", dpi=120)

正常情况下你会看到反应器温度围绕某个稳态小幅波动，冷却水阀开度也在稳态附近微调——这就是那个闭环控制回路在工作。如果温度一路单调漂走或者剧烈震荡，八成是构造参数给得太离谱，比如噪声给了个 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_ok=True) # fault=0 是正常工况，1..21 是 21 种预设故障 for fault in range(0, 22): # 每个故障单独构造一个模型实例，避免状态从上一次仿真里带过来 model = tep.TEP(fault=fault, xmeas_noise=0.1, xmv_noise=0.0) xmeas, xmv = model.simulate(n_samples=N_SAMPLES, timestep=TIMESTEP) # 标签单独存一份，方便后面拼成分类任务的 y label = 0 if fault == 0 else 1 np.savez_compressed( os.path.join(OUT_DIR, f"fault_{fault:02d}.npz"), xmeas=xmeas.astype(np.float32), xmv=xmv.astype(np.float32), fault=fault, label=label, ) print(f"fault={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(fault=3, xmeas_noise=0.1, xmv_noise=0.0)

第 4 步容易被忽略但很关键。控制器里通常有积分环节，如果积分初值给 0 而实际偏差非 0，前几十个采样点会出现一段明显的“爬升”暂态。这段暂态在真实数据里是存在的（开车阶段），但如果你要做稳态故障检测，它就是个干扰项。这也是为什么标准实验流程要求先跑一段无故障预热再引入故障。

对象一旦构造完，内部就保存了完整的状态向量。simulate()不重新初始化，它从当前状态往下接着算——这一点后面避坑章节会专门讲。

4.2 单步推进：一个采样周期内的时序

simulate(n_samples=500)内部是一个 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_samples=1000, timestep=0.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_samples=960, timestep=0.05): np.random.seed(seed) # 固定全局随机源 model = tep.TEP(fault=fault, xmeas_noise=0.1) return model.simulate(n_samples=n_samples, timestep=timestep)

如果包内部用的是np.random.randn这类全局函数，np.random.seed就够了；如果它自己持有一个RandomState或者Generator实例，那就需要在构造时把种子传进去。两种写法都试一遍，跑两次看输出是否逐位相等，就能判断出来。验证方法很简单：同一段代码跑两次，np.array_equal(a, b)返回 True 才算真复现。

稳态工况下曲线看着差不多，不代表数值相同——噪声幅值 0.1 的时候肉眼看不出差别，但对检测算法来说，那点差异可能正好跨过阈值。所以别靠眼睛判断复现性。

5. 避坑与排查：数据看着不对劲时先查这五处

5.1 现象：曲线变化快得离谱，故障几秒就“发作”完了

原因：时间单位搞错了。timestep的单位是小时，不是秒也不是分钟。写timestep=1想表达 1 秒，实际一步推进了 1 小时，等于标准采样间隔的 20 倍。积分步长变大后，流程的动态响应被压缩到少数几个采样点里，看起来就是“瞬间爆发”。

解决：默认用0.05，也就是 3 分钟，跟标准 TEP 实验对齐。要算总时长就用n_samples * timestep小时，960 个点对应 48 小时。改步长之前先问自己一句：我要的是更细的时间分辨率，还是更长的仿真时长？这两个诉求不能靠调一个参数同时满足。

5.2 现象：明明设了fault=5，前 20 个采样点跟无故障跑出来一模一样

原因：故障不是从第 0 个采样点开始的。标准 TEP 实验规程是先在正常工况下跑 1 小时（20 个采样点）建立稳态，然后再引入故障。pyTEP 通常遵循这个约定，所以前 20 个点是干净的正常数据。

解决：造数据集时保留这段预热区，训练时可以选择性丢弃前 20 个点；评估检测延迟时，把故障起点定义在第 20 个采样点而不是第 0 个。如果你只跑 20 个点就结束，那整个样本里根本没有故障——这个坑我踩过，排查了半天以为是故障注入失效，其实就是样本太短。

5.3 现象：连续调两次simulate()，第二次的数据开头跟第一次结尾接不上，或者数据分布明显偏移

原因：simulate()是增量式的，它从对象当前保存的状态接着往下算，不做重置。第二次调用的起点是第一次的终点，而不是初始稳态。如果你的控制器在这期间累积了积分偏差，第二次的起点可能已经偏离稳态很远了。

解决：需要独立样本时，每个样本重新tep.TEP(...)构造一个新对象，别复用。这个成本不高，构造本身不怎么耗时，真正贵的是积分循环。要写循环造数据的话，把model = tep.TEP(...)放进循环体内，这一行别省。

5.4 现象：拿 pyTEP 的输出跟公开的 TEP 数据集对拍，列数对不上，41+12 和 41+11 差了 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之前，都强制先跑一遍这个基线快照，把改动前后的数字并排打出来看。这个习惯看着笨，但它是我在这个项目上唯一没后悔过的规矩。希望帮到你。

本文还有配套的精品资源，点击获取