pyTEP：Python化工仿真与TEP故障注入闭环实践

简介：pyTEP 是面向过程控制、仿真建模与故障诊断研究的开源 Python 工具包，封装田纳西伊士曼（Tennessee Eastman）流程的交互式仿真接口。该工具包让使用者无需理解模拟器底层实现，即可通过高级 API 完成仿真初始化、操作条件调整与运行数据存储，适合化工过程教学、控制算法验证及机器学习实验等场景，对具备 Python 基础的中高级学习者尤为友好。资源包共 76 个文件，约 754KB，以 16 个 py 源码文件为主体，配合 23 个 m 文件及 mexw64、mexa64 二进制实现 MATLAB 桥接，另含 pkl、mat 数据文件、md 与 rst 文档及构建配置，目录覆盖模拟器、仿真接口、工具函数和测试模块。已有 546 人学习下载。借助其中的 API 源码、示例脚本与桥接模块，读者可快速搭建可复现的 TEP 仿真环境，理解仿真数据组织方式，并在此基础上扩展控制器、采集运行日志或对比不同工况，减少从零配置复杂模拟场景的成本。

1. 化工仿真基准搬进 Python：pyTEP 能干什么

第一次打开田纳西伊士曼过程（TEP）的原始仿真代码时，我花了半天才在密密麻麻的 COMMON 块里找到那个控制反应器温度的阀位变量；改一个故障参数要重新编译，想看中间状态还得自己往文件里写。pyTEP 干的事就是把这一整套搬到 Python 里：它用 Python 重写了 TEP 的工艺对象、控制回路和故障开关，几十行代码就能拉起一次带故障注入的闭环仿真。

它面向的不是跑分的人，而是做过程监控、故障诊断、软测量和数据集生成的开发者。你要 41 个测量变量和 12 个操纵变量的时序数据，它按你给的采样节奏吐出来；你要复现某个故障工况，它给你一个开关。对没接触过化工流程的 Python 开发者来说，它也是一份结构相当工整的仿真引擎样本，读源码比读论文快。

2. 项目结构与运行环境：从克隆到跑出第一帧数据

拿到包先别急着写循环，把目录翻开看一眼。pyTEP 的结构不大，核心就两三个文件，剩下的是例子和测试。搞清哪一层负责过程模型、哪一层负责控制器，后面改参数才知道该动哪个文件。

2.1 源码目录里的分工

常见的组织方式是这样，不同版本可能在文件名上有出入，但层次基本一致：

pytep/ ├── pytep/ │ ├── __init__.py # 对外导出 Plant、Controls 等入口 │ ├── plant.py # 工艺对象：状态向量、物料衡算、故障开关 │ ├── control.py # 控制器：标准控制方案的 PI 回路集合 │ └── utils.py # 单位换算、数据记录的辅助函数 ├── examples/ # 开环、闭环、故障注入的示例脚本 ├── tests/ # 回归测试，改模型后先跑这里 ├── setup.py / pyproject.toml └── README.md

plant.py是整个项目的重点，里面一个Plant类装着过程的全部状态：反应器温度、压力、液位、各组分浓度，以及 12 个阀位。control.py里的Controls类不持有工艺状态，它只读测量值、算阀位，属于“贴在对象外面的壳”。这种把被控对象和控制器分开的做法，是我在仿真项目里最愿意看到的写法——换成自己的控制器时，不用去动工艺模型那一坨。

2.2 依赖与安装

纯 Python 包，依赖很轻，主力就是 NumPy，画图和做分析时才需要 Matplotlib。两种装法：

# 方式一：直接从源码目录装成可编辑模式，改完代码立即生效 git clone <项目地址> pytep cd pytep pip install -e . # 方式二：只想跑一下示例，装依赖即可 pip install numpy matplotlib

pip install -e .这个-e是关键，做二次开发时一定要用可编辑安装，否则你改完plant.py发现行为没变，会怀疑人生——其实就是装到了 site-packages 里的那份旧拷贝。

提示：项目里的依赖版本尽量别自己升到最新。仿真代码对 NumPy 的数组广播行为比较敏感，老项目在新版 NumPy 上偶尔会冒出形状不匹配的警告，先把版本锁在 README 写的那一档。

2.3 最小可运行脚本：从稳态出发走一步

先不急着接控制器，跑一个开环步进，确认模型本身能转起来：

import numpy as np from pytep.plant import Plant from pytep.control import Controls # 1. 实例化工艺对象和控制器 plant = Plant() ctrl = Controls(plant) # 2. 设置时间参数 plant.set_time(0.1) # 采样周期，单位与模型一致（TEP 里是小时） plant.set_integration_time(0.001) # 内部积分步长，必须明显小于采样周期 # 3. 把对象置到稳态初值，否则一堆变量是 0，第一步就会算飞 plant.set_steady_state() # 4. 手动推进 10 个采样周期，观察测量值 for k in range(10): plant.step() print(k, plant.xmeas[6], plant.xmeas[8], plant.xmv[9])

这段代码干了四件事：建对象、定时间尺度、置初值、推进。set_time是采样周期，也就是你“多久看一次数据”；set_integration_time是模型内部数值积分的步长，必须比采样周期小一个量级以上，否则一个采样周期内只积一两步，结果会明显失真。set_steady_state把全部状态量拉到一个已经收敛的工作点，跳过从冷态开车那段又慢又容易发散的过渡过程。

plant.xmeas是长度 41 的测量值数组，索引从 0 开始，所以xmeas[6]对应文档里的xmeas_7（反应器压力），xmeas[8]是反应器温度，plant.xmv[9]是反应器冷却水阀位。这个“文档编号从 1 起、数组下标从 0 起”的错位是新手第一个绊脚石，我在下面还会专门说。

如果这十行能跑出稳定的数值，说明环境和模型都没问题，可以往上叠控制器了。如果第一步就出现nan或者数值炸到 1e10，多半是漏了set_steady_state，或者积分步长设得比采样周期还大。

3. 核心模块拆解：一个闭环步进里谁先动

把开环跑通之后，真正要用的是闭环。TEP 的经典配置是 20 来个 PI 回路，pyTEP 把这套标准控制方案打包进了Controls。理解它的关键不是背回路参数，而是搞清楚一次步进里“测量—控制—积分”谁先谁后。

3.1 Plant：41 个测量值和 12 个阀位从哪来

Plant对外暴露的接口基本就是三个数组和几个方法：

xmeas里前 22 个是连续测量，后面是成分分析，采样间隔更长、带延迟，这是 TEP 本身的设计，不是 pyTEP 的 bug。你在做软测量或者故障诊断时如果发现某几个通道看起来“卡住不动”，先确认是不是成分分析通道还没到更新时刻。

xmv的 12 个通道可以直接赋值。开环实验里最常见的做法是把某个阀位钉死，看系统怎么反应：

# 把反应器冷却水阀固定在一个偏离稳态的位置 plant.xmv[9] = plant.xmv[9] * 0.9 for _ in range(100): plant.step()

注意这里改的是“阀位指令”，模型内部还有阀门动态，所以第一拍的反应不会立刻体现，通常要滞后几个采样点才看得到温度变化。这也是验证模型有没有正常工作的土办法：动阀位、看温度，方向对不上就是索引搞错了。

3.2 Controls：标准控制方案与手动接管

Controls在构造时绑定了Plant实例，内部维护一组 PI 控制器的设定值、增益和积分时间。它的默认行为是维持 TEP 的标准控制方案：反应器温度靠冷却水阀调，压力靠排放阀调，各液位有自己的回路。

最常见的用法是：大部分回路让它自动跑，只接管其中一两个。

# 把反应器温度设定值抬高一点，看控制器怎么把冷却水阀拉大 ctrl.set_xmeas_sp(8, 122.0) # 第 9 个测量量（反应器温度）的设定值 for k in range(200): ctrl.set_external_measurement(plant.xmeas) # 把当前测量喂给控制器 ctrl.step() # 控制器算一拍，写 xmv 的设定 plant.step() # 对象向前积分

上面两个方法名以你手上源码为准，set_xmeas_sp这类接口在不同版本里可能叫set_setpoint或在ctrl.sp数组上直接赋值。但顺序是不变的：先把测量值送进控制器，控制器更新阀位，对象再积分。顺序反了就会出现“用上一拍的阀位算这一拍的测量”，闭环延迟凭空多一拍，调参时怎么调都觉得控制器慢半拍。

ctrl.step()之后，plant.xmv会被控制器的输出覆盖。所以如果你想做手动接管，比如在某个故障工况下强制关小某个阀，必须放在ctrl.step()之后、plant.step()之前，否则会被下一拍的控制输出冲掉。这个顺序问题我在两个项目里都踩过，第一次排查了整整一下午。

3.3 一次闭环步进的三段链路

把上面的顺序固定成模板，写进自己的脚本里：

def run_closed_loop(plant, ctrl, steps, on_step=None): log = [] for k in range(steps): # 第一段：测量，把对象状态交给控制器 ctrl.set_external_measurement(plant.xmeas) # 第二段：控制，控制器更新阀位 ctrl.step() # 第三段：可选的外部干预，比如故障注入、手动阀位覆盖 if on_step is not None: on_step(k, plant, ctrl) # 第四段：积分，对象推进一个采样周期 plant.step() # 记录：把这一拍的测量和阀位存下来 log.append(np.concatenate([plant.xmeas, plant.xmv])) return np.array(log)

这个函数是整个二次开发的地基。on_step这个钩子看着不起眼，但它是故障注入、设定值调度、控制器切换的统一入口——所有“在第 k 步做点什么”的需求都能塞进去，不用去改循环体。返回值是一个steps × 53的矩阵，前 41 列是测量，后 12 列是操纵变量，直接喂给 pandas 或者 sklearn 都行。

steps这个参数的取值要配合采样周期来算。TEP 的时间单位是小时，set_time(0.1)意味着每拍 0.1 小时，跑 200 拍就是 20 小时的过程时间。故障诊断任务里，阶跃型故障通常要观察 8 到 10 小时才看得出明显偏离，所以一次实验至少 100 拍起步，别指望跑 20 拍就能看出差别。

4. 故障注入与数据采样：把 21 种工况真正用起来

pyTEP 对我来说最大的价值不是能跑仿真，而是能批量造带标签的数据。TEP 自带 21 种预设扰动，覆盖进料浓度阶跃、进料温度变化、冷却水温度变化、反应动力学漂移、阀门粘滞、阀门固定等几大类，正好对应故障诊断里最常见的几种场景。

4.1 故障开关的注入时机

故障通过idv向量控制，长度 21，把对应位置置 1 就打开：

def inject_fault(k, plant, ctrl): # 在第 50 拍打开 4 号故障（反应器冷却水入口温度阶跃） if k == 50: plant.idv[3] = 1 # 在第 150 拍撤销，观察系统恢复过程 if k == 150: plant.idv[3] = 0 log = run_closed_loop(plant, ctrl, steps=300, on_step=inject_fault)

注入时机是有讲究的。前 50 拍留给系统在稳态附近跑一段“正常数据”，这段就是你后面做对比的基准；第 50 拍注入，故障开始发酵；第 150 拍撤销，能拿到一段恢复过程，对做故障隔离特别有用。如果一上来就注入，模型还在往稳态收敛，你根本分不清数据里的波动是故障引起的还是初始过渡引起的。

plant.idv这个字段名也要以源码为准，有的版本把它包装成plant.set_fault(index)这样的方法。打开源码搜一下idv这个字符串，一秒就能确认。

4.2 采样周期和积分步长的配比

这两个参数是新手最容易配错的一对。经验值是积分步长比采样周期小一到两个数量级：

配错的后果很直观：积分步长太接近采样周期，模型在一个采样周期里只积了几步，数值积分误差被放大，压力、温度这些快变量会出现明显的高频抖动；积分步长太小、采样周期又大，计算时间会成倍上涨，但数据精度提升很有限。

一个实用的检查方法：把同一组故障跑两遍，只改积分步长，如果两次结果的均值差超过 1%，说明步长还不够小。

4.3 批量生成带标签的数据集

单个实验跑通之后，批量造数据就是套循环的事。核心是把“工况描述”和“故障标签”一起存下来：

import json import numpy as np def make_dataset(fault_ids, steps=300, onset=50, seed_base=0): samples, labels = [], [] for fid in fault_ids: for rep in range(5): np.random.seed(seed_base + fid * 100 + rep) plant, ctrl = Plant(), Controls(plant) plant.set_time(0.1) plant.set_integration_time(0.001) plant.set_steady_state() def hook(k, p, c, fid=fid): if k == onset and fid > 0: p.idv[fid - 1] = 1 log = run_closed_loop(plant, ctrl, steps, on_step=hook) samples.append(log) labels.append(fid) np.savez_compressed( "tep_dataset.npz", X=np.array(samples), # 形状 (样本数, 时间步, 53) y=np.array(labels), # 每个样本对应的故障编号，0 表示正常 )

fault_ids传[0, 1, 4, 8, 13, 21]这样的列表，0代表无故障样本，其余对应具体故障。每个工况跑 5 次是为了引入随机初值、避免模型只看过一条轨迹，虽然set_steady_state之后初值差不多，但加上微小扰动还是有帮助的。

np.savez_compressed比直接np.save省不少空间，三维数组压缩率很高。存成(样本, 时间步, 通道)这种形状，后面切训练集和验证集的时候一行X[:n]就完事，不用来回转置。

5. 避坑排查：跑不起来、数值发散、结果对不上

这个项目的坑不算多，但集中在那几个地方，踩一次记一辈子。下面几条是我实际遇到过的，按现象、原因、解决写。

5.1 环境和调用层面的坑

现象：ModuleNotFoundError: No module named 'pytep'，但明明pip install过了。原因基本只有一个——你装的是 PyPI 上的旧版本，而你在项目目录里改的是另一份源码。解决方式是用pip install -e .从源码目录装，或者干脆把示例脚本放在项目根目录下跑，让 Python 优先加载当前目录的包。装完用python -c "import pytep; print(pytep.__file__)"看一眼路径对不对，这一步能省掉大量“我改了代码怎么没生效”的困惑。

现象：闭环跑起来之后，阀位纹丝不动，测量值也不变。原因是漏了把测量值送进控制器，或者顺序写反了。Controls需要显式接收当前的xmeas才能算输出，它不是自动绑定的。检查你的循环里有没有ctrl.set_external_measurement(plant.xmeas)这一行，以及它是不是在ctrl.step()之前。

现象：IndexError或某个通道一直在 0 附近不动。典型的下标错位。文档里说“第 7 个测量量是反应器压力”，对应xmeas[6]而不是xmeas[7]。另外xmv只有 12 个，别拿xmeas的索引去索引它。最快的核对方式是把plant.xmeas打印出来，对照文档里的量程看哪几个值明显不合理——比如压力应该是 2700 到 3000 之间，如果某个通道是 0.0 或者 1e-3 这种量级，基本可以确定索引偏了。

5.2 数值和结果层面的坑

现象：跑几十拍之后出现nan，或者数值直接冲到 1e8。原因是没置稳态初值，或者积分步长设得太大。冷态启动时，反应器里没有物料，各组分浓度都是 0，物料衡算里的除法立刻出问题。解决方式是先plant.set_steady_state()，再把set_integration_time调到采样周期的十分之一以下。如果已经出了nan，把积分步长减半再跑一遍，同时打印每拍的xmeas前几项，定位是哪一步开始崩的。

现象：故障注入了，但系统响应几乎看不出来。原因可能是故障设得太晚、观察窗口太短，或者故障本身在你的采样周期下变化太慢。TEP 里的阀门粘滞、动力学漂移这类故障本身就是慢过程，跑 30 拍确实看不出差异。解决方式是延长仿真时长、把故障注入点提前，或者换成阶跃型故障（进料浓度、冷却水温度这类）先验证注入链路通不通，确认通了再换慢故障。

现象：两次跑的同一组参数，结果差得离谱。原因是模型内部有随机初始扰动，或者控制器的积分项带着上一轮的残余。解决方式是在每次实验前显式np.random.seed()，并且每次都重新Plant()和Controls()实例化，不要复用对象。复用对象看起来省事，实际上前一次的状态还挂在里面，第二次实验的起点就变了。

6. 二次开发：把控制器换掉，让仿真为你的算法服务

用熟了之后，我最常做的改造是把Controls换掉。标准控制方案适合生成基准数据，但如果你要验证自己写的模型预测控制或者强化学习策略，就得让控制器完全听你的。

最省事的做法是保留Plant，绕开Controls，自己写一个只含step()的轻量对象，接口对齐set_external_measurement和step两个方法就行：

class MyController: """只保留 pyTEP 控制器的最小接口，内部换成自己的算法""" def __init__(self, plant, target_temp=122.0): self.plant = plant self.target = target_temp self.integral = 0.0 self.kp, self.ki = 0.8, 0.02 def set_external_measurement(self, xmeas): # 控制器拿到的就是当前这一拍的 41 个测量值 self.xmeas = xmeas.copy() def step(self): # 单回路温度控制：误差驱动冷却水阀 err = self.target - self.xmeas[8] self.integral += err * self.plant.get_time() u = self.kp * err + self.ki * self.integral # 阀位有物理上下限，必须夹住，否则会送出非法值让模型崩 self.plant.xmv[9] = min(max(u, 0.0), 100.0)

替换进去只需要改一行ctrl = MyController(plant)，run_closed_loop完全不用动——这就是提前把循环写成模板的好处。几个我踩过的细节：self.xmeas = xmeas.copy()一定要 copy，否则你持有的是对象内部数组的引用，下一拍被覆盖后你的控制器看到的“历史”全是最后一次的值；阀位输出必须夹在物理范围内，TEP 的阀位超出范围不会报错，而是悄悄算出离谱的流量，几拍之后整个模型就崩了；积分项要做抗饱和，不然设定值大幅变化时积分会一路累积，撤掉扰动后阀位很久回不来。

验证方法很简单：先让MyController复现标准控制方案的行为——把kp、ki调到和Controls里反应器温度回路差不多，跑同一组故障，两条曲线应该基本重合。重合了就说明接口接对了，再逐步改参数、加算法。如果一开始就上复杂算法、结果不理想，你会分不清是算法不行还是接口接错了，这时候后悔药可没有。

另外一个小技巧：把Plant的实例化和时间参数设置抽成工厂函数，批量实验的时候每个样本都新建一个干净对象，别为了省内存复用。仿真对象的内存开销相比一次跑错的代价，完全不值一提。

从那以后我每次改控制器，都强制先跑一遍“复现标准方案”的对照实验，两条曲线不重合就不往下走。这套笨办法救过我至少三次。希望帮到你。

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