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. 化工仿真基准搬进 PythonpyTEP 能干什么第一次打开田纳西伊士曼过程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.mdplant.py是整个项目的重点里面一个Plant类装着过程的全部状态反应器温度、压力、液位、各组分浓度以及 12 个阀位。control.py里的Controls类不持有工艺状态它只读测量值、算阀位属于“贴在对象外面的壳”。这种把被控对象和控制器分开的做法是我在仿真项目里最愿意看到的写法——换成自己的控制器时不用去动工艺模型那一坨。2.2 依赖与安装纯 Python 包依赖很轻主力就是 NumPy画图和做分析时才需要 Matplotlib。两种装法# 方式一直接从源码目录装成可编辑模式改完代码立即生效 git clone 项目地址 pytep cd pytep pip install -e . # 方式二只想跑一下示例装依赖即可 pip install numpy matplotlibpip 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 Plant41 个测量值和 12 个阀位从哪来Plant对外暴露的接口基本就是三个数组和几个方法属性/方法维度含义xmeas41测量变量流量、温度、压力、液位、组分分析xmv12操纵变量各进料阀、排放阀、冷却水阀、搅拌转速idv21故障开关向量对应 TEP 的 21 种预设扰动set_time(dt)—设定采样周期set_integration_time(dt)—设定内部积分步长set_steady_state()—把状态置为标准稳态工作点step()—按采样周期推进一次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_stepNone): 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, steps300, on_stepinject_fault)注入时机是有讲究的。前 50 拍留给系统在稳态附近跑一段“正常数据”这段就是你后面做对比的基准第 50 拍注入故障开始发酵第 150 拍撤销能拿到一段恢复过程对做故障隔离特别有用。如果一上来就注入模型还在往稳态收敛你根本分不清数据里的波动是故障引起的还是初始过渡引起的。plant.idv这个字段名也要以源码为准有的版本把它包装成plant.set_fault(index)这样的方法。打开源码搜一下idv这个字符串一秒就能确认。4.2 采样周期和积分步长的配比这两个参数是新手最容易配错的一对。经验值是积分步长比采样周期小一到两个数量级采样周期积分步长适用场景0.10.001常规故障诊断数据集数值稳定0.050.0005需要更细时间分辨率的小幅故障0.50.005长时间慢漂移省算力配错的后果很直观积分步长太接近采样周期模型在一个采样周期里只积了几步数值积分误差被放大压力、温度这些快变量会出现明显的高频抖动积分步长太小、采样周期又大计算时间会成倍上涨但数据精度提升很有限。一个实用的检查方法把同一组故障跑两遍只改积分步长如果两次结果的均值差超过 1%说明步长还不够小。4.3 批量生成带标签的数据集单个实验跑通之后批量造数据就是套循环的事。核心是把“工况描述”和“故障标签”一起存下来import json import numpy as np def make_dataset(fault_ids, steps300, onset50, seed_base0): 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, fidfid): if k onset and fid 0: p.idv[fid - 1] 1 log run_closed_loop(plant, ctrl, steps, on_stephook) samples.append(log) labels.append(fid) np.savez_compressed( tep_dataset.npz, Xnp.array(samples), # 形状 (样本数, 时间步, 53) ynp.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_temp122.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的实例化和时间参数设置抽成工厂函数批量实验的时候每个样本都新建一个干净对象别为了省内存复用。仿真对象的内存开销相比一次跑错的代价完全不值一提。从那以后我每次改控制器都强制先跑一遍“复现标准方案”的对照实验两条曲线不重合就不往下走。这套笨办法救过我至少三次。希望帮到你。本文还有配套的精品资源点击获取