简介pyTEP 是面向过程控制、化工仿真与机器学习研究者的开源 Python 仿真 API用于田纳西伊士曼Tennessee Eastman流程的交互式模拟。它通过高级接口简化复杂工况的设置让用户无需深入理解模拟器底层原理即可调整操作条件、运行仿真并存储数据适合控制算法验证、故障诊断与数据驱动建模等场景。资源包共 76 个文件约 754KB包含 16 个 py 源码、23 个 m 文件、8 个 pkl 数据、7 个 mdl 模型及 mat、mexw64、mexa64 等混合文件覆盖 Python 接口、MATLAB 桥接、仿真引擎与测试模块并附 README、许可证与构建配置。目前已有 532 人学习下载。借助该包读者可快速搭建 TEP 仿真环境复现标准工况数据理解 siminterface、engineutils 等模块的调用方式为控制策略对比与异常检测实验提供可复用的工具基础。1. 用 pyTEP 把 Tennessee Eastman 流程装进 Python为什么它值得你花一个下午跑通如果你做过程控制、故障诊断或者强化学习在工业场景的落地大概率绕不开 Tennessee Eastman 流程。这个 1993 年就被提出的化工仿真基准至今仍是论文里的常客原因很简单它足够复杂又足够标准。但真正动手时你会发现很多现成的仿真环境要么绑死在商业软件里要么接口老旧、改一个参数要翻半天文档。pyTEP 这个开源 API 的出现把这件事拉回到了 Python 生态里——用几行代码就能建起一个可交互的闭环仿真改工况、加扰动、接控制器都在同一个脚本里完成。它适合三类人做数据驱动故障诊断、需要批量生成带标签工况数据的研究者做 MPC、RL 控制器验证、想要一个可控被控对象的工程师以及教学场景里需要让学生直观看到多变量耦合效应的讲师。这篇文章不讲空泛概念直接按「装什么、怎么建、参数怎么调、哪里会翻车」的顺序走一遍。2. 环境搭建与最小可运行闭环从零到第一条趋势曲线2.1 安装路径选择与依赖确认pyTEP 本质是对底层 Fortran 仿真核心的 Python 封装所以安装时最容易出问题的不是 Python 包本身而是编译工具链。常见做法有两种一是直接用 pip 从源码构建二是先拿到预编译的仿真核心再装 Python 层。我一般推荐后者因为省去 Fortran 编译器版本匹配的麻烦。# 建议先建独立环境避免和现有科学计算栈冲突 python -m venv pytep_env source pytep_env/bin/activate # Windows 下用 pytep_env\Scripts\activate # 升级构建工具老版本 setuptools 经常在编译扩展时翻车 pip install --upgrade pip setuptools wheel # 安装 pyTEP若本地已有仿真核心可加 --no-build-isolation pip install pytep装完后别急着写业务代码先做一次导入自检。这一步能过滤掉 80% 的「装上了但跑不起来」问题。# check_env.py import pytep import numpy as np # 打印版本和可用接口确认核心库被正确加载 print(pyTEP version:, getattr(pytep, __version__, unknown)) print(numpy version:, np.__version__) # 尝试实例化一个最简仿真对象不设任何扰动 sim pytep.TEP() print(simulation object created:, type(sim).__name__)逻辑说明这段代码只做三件事——确认包能导入、确认数值库版本兼容、确认仿真对象能实例化。参数上不需要改任何东西如果这一步报ImportError或OSError基本是编译产物没找到回去检查 Fortran 运行库路径。如果实例化时报参数错误说明 API 签名和文档有出入用help(pytep.TEP)看实际入参。2.2 构建第一个闭环仿真设定工况与采样pyTEP 的核心用法是「先定义仿真配置再逐步推进时间」。它不像某些环境那样一次性跑完返回数组而是给你一个可以逐拍交互的对象这对做在线控制非常关键。import pytep import numpy as np import matplotlib.pyplot as plt # 创建仿真实例指定采样周期为 0.01 小时约 36 秒 # 这个值对应 TEP 标准文献里的控制周期改大改小都会影响动态响应 sim pytep.TEP(sampling_time0.01) # 设定初始工况这里用标准工况 1即基本进料比例 sim.set_mode(1) # 初始化让系统进入稳态 sim.initialize() # 准备记录容器 time_log [] reactor_pressure_log [] # 推进 200 个采样点每步读取反应器压力 for step in range(200): sim.step() # 推进一个采样周期 time_log.append(step * 0.01) # 变量索引需查 pyTEP 变量表反应器压力通常是第 7 个测量量 reactor_pressure_log.append(sim.get_variable(reactor_pressure)) # 画趋势图确认仿真在正常波动而非发散 plt.plot(time_log, reactor_pressure_log) plt.xlabel(Time (h)) plt.ylabel(Reactor Pressure) plt.title(pyTEP Startup Trend) plt.grid(True) plt.show()逻辑说明set_mode(1)决定了进料和产物的基础配比TEP 一共定义了 6 种标准工况模式 1 是最常用的基准。initialize()会跑一段内部稳态计算不调用它直接step()可能拿到未收敛的初值。get_variable的字符串名需要和 pyTEP 的变量字典对齐不同版本命名可能有差异稳妥做法是先打印sim.list_variables()看一遍。参数说明sampling_time不建议随意改因为 TEP 的控制器参数和扰动设计都是基于 0.01 小时这个基准。如果你做强化学习需要更细的时间粒度正确做法是在step()之间做插值而不是改采样周期。step()每次推进一个周期返回的是当前状态字典还是 None 取决于版本别假设它有返回值。2.3 变量索引与测量量映射TEP 有 41 个测量变量和 12 个操纵变量新手最容易在这里迷路。pyTEP 一般提供名称到索引的映射但名称规范不一定和你的论文一致。# 打印所有可读变量名建立自己的映射表 variables sim.list_variables() for idx, name in enumerate(variables): print(f{idx:3d} {name}) # 常见做法把关键变量名存成常量避免代码里到处写魔法字符串 REACTOR_PRESSURE reactor_pressure REACTOR_LEVEL reactor_level SEPARATOR_LEVEL separator_level PRODUCT_FLOW product_flow # 批量读取多个变量减少函数调用开销 state sim.get_variables([REACTOR_PRESSURE, REACTOR_LEVEL, PRODUCT_FLOW]) print(current state:, state)逻辑说明list_variables()返回的顺序在不同版本里可能不同所以不要硬编码索引用名称查。get_variables批量读取比循环单次读取快做长时间仿真时这个差异会累积。如果你发现某个变量名不存在先确认仿真模式是否支持该变量——某些工况下部分测量量是冻结的。3. 扰动注入与故障数据生成把 TEP 变成你的数据工厂3.1 内置扰动类型与触发时机TEP 的价值一半在正常工况一半在它预定义的 20 多种故障。pyTEP 把这些故障做成了可编程接口你可以在任意时刻注入也可以设定持续时长。import pytep import numpy as np sim pytep.TEP(sampling_time0.01) sim.set_mode(1) sim.initialize() # 先跑 100 步让系统稳定 for _ in range(100): sim.step() # 在第 100 步注入故障 ID 4反应器冷却水入口温度阶跃 # 这个故障在文献里属于较难检测的类型适合测试诊断算法 sim.set_fault(4) # 再跑 300 步记录故障后的动态 fault_data [] for _ in range(300): sim.step() fault_data.append(sim.get_variables([ reactor_temperature, reactor_cooling_water_outlet_temp, product_flow ])) fault_data np.array(fault_data) print(fault data shape:, fault_data.shape)逻辑说明set_fault(4)是即时生效的不需要额外调用step()来激活。故障 ID 和 TEP 标准表对应但 pyTEP 可能对部分故障做了简化建议先用一个已知容易检测的故障比如 ID 1验证流程再换难检测的。fault_data的 shape 是(300, 3)可以直接喂给 sklearn 或 pytorch 做分类。参数说明故障注入的强度在 pyTEP 里通常是固定的对应原始 TEP 的设计值。如果你想做变强度故障常见做法是叠加一个缩放因子到测量量上而不是改故障本身——因为改故障定义会破坏和文献的可比性。3.2 批量生成带标签数据集做故障诊断模型训练你需要的是成千上万条带标签的轨迹而不是一条。pyTEP 支持重置和重新初始化但要注意重置的彻底性。import pytep import numpy as np def generate_trajectory(fault_id, n_steps500, seedNone): 生成一条带故障标签的轨迹 if seed is not None: np.random.seed(seed) sim pytep.TEP(sampling_time0.01) sim.set_mode(1) sim.initialize() # 稳态段 for _ in range(50): sim.step() # 注入故障fault_id 为 0 表示正常工况 if fault_id 0: sim.set_fault(fault_id) data [] for _ in range(n_steps): sim.step() data.append(sim.get_variables(sim.list_variables())) return np.array(data), fault_id # 批量生成正常 5 种故障每种 20 条 dataset [] labels [] for fid in [0, 1, 2, 4, 8, 13]: for run in range(20): traj, label generate_trajectory(fid, n_steps500, seedrun) dataset.append(traj) labels.append(label) dataset np.array(dataset) labels np.array(labels) print(dataset shape:, dataset.shape) # (120, 500, 41) print(labels shape:, labels.shape)逻辑说明每条轨迹都重新实例化TEP()这是最干净的做法避免状态残留。seed控制的是 numpy 的随机数pyTEP 内部如果有随机扰动需要确认它是否也接受种子——如果没有轨迹之间的差异只来自故障本身多样性会受限。list_variables()返回全部 41 个测量量实际建模时通常要做特征选择因为很多变量高度相关。参数说明n_steps500对应 5 小时仿真时长对慢动态故障足够对快动态故障可能太长。稳态段 50 步是经验值太短系统没稳太长浪费计算。如果你要做早期故障检测故障后的前 50 步才是关键窗口后面的数据反而会稀释信号。3.3 与 scikit-learn 对接的预处理原始 TEP 数据不能直接扔给模型标准化和滑窗是标配。from sklearn.preprocessing import StandardScaler from sklearn.model_selection import train_test_split # dataset shape: (n_samples, n_timesteps, n_features) n_samples, n_timesteps, n_features dataset.shape # 展平后做标准化再还原回三维 flat dataset.reshape(-1, n_features) scaler StandardScaler() flat_scaled scaler.fit_transform(flat) dataset_scaled flat_scaled.reshape(n_samples, n_timesteps, n_features) # 用滑窗构造样本窗口长度 20步长 5 window 20 stride 5 X, y [], [] for i in range(n_samples): for t in range(0, n_timesteps - window, stride): X.append(dataset_scaled[i, t:twindow, :]) y.append(labels[i]) X np.array(X) y np.array(y) print(X shape:, X.shape) # (n_windows, 20, 41) print(y shape:, y.shape) # 划分训练测试集注意按轨迹划分避免泄漏 X_train, X_test, y_train, y_test train_test_split( X, y, test_size0.3, random_state42, stratifyy )逻辑说明标准化必须用训练集的均值和方差但这里为了简化先在全量上 fit实际项目里要严格分开。滑窗的window20对应 0.2 小时对 TEP 的多数故障来说足够捕捉动态。stride5是为了减少窗口间的相关性如果做在线检测可以设 stride1。按轨迹划分比随机划分更严格因为同一条轨迹的窗口高度相似随机划分会导致测试集泄漏。4. 控制器接入与交互式仿真让 pyTEP 跑在你的算法里4.1 读取操纵变量并写入控制量pyTEP 的交互性体现在你可以每步读取状态、计算控制量、再写回去。这是做 MPC 或 RL 的标准循环。import pytep import numpy as np sim pytep.TEP(sampling_time0.01) sim.set_mode(1) sim.initialize() # 简单比例控制器控制反应器温度到设定值 setpoint 120.0 Kp 2.0 for step in range(500): # 读取当前反应器温度 current_temp sim.get_variable(reactor_temperature) # 计算控制量这里操纵冷却水流量 error setpoint - current_temp control Kp * error # 写入操纵变量注意要遵守物理上下限 control np.clip(control, 0.0, 100.0) sim.set_variable(reactor_cooling_water_flow, control) # 推进仿真 sim.step() if step % 50 0: print(fstep {step}: temp{current_temp:.2f}, control{control:.2f})逻辑说明set_variable必须在step()之前调用否则控制量会在下一步才生效造成一步延迟。np.clip是必须的因为 pyTEP 对超出物理范围的操纵变量可能直接报错或产生非物理行为。比例控制器的Kp2.0是随手设的实际整定需要看阶跃响应——TEP 的反应器温度回路增益不高Kp 太大容易振荡。参数说明setpoint120.0是反应器温度的典型设定值但不同工况下最优值不同。操纵变量名reactor_cooling_water_flow需要和list_variables()里的名称一致。如果你发现写入后温度没反应先确认这个变量在当前模式下是否可操纵——有些变量在特定工况下是固定的。4.2 与 Gym 接口对接做强化学习做 RL 的人希望 pyTEP 长得像 Gym 环境。自己包一层是最稳妥的因为现成的包装器往往版本滞后。import gym from gym import spaces import pytep import numpy as np class TEPEnv(gym.Env): 把 pyTEP 包装成 Gym 环境 def __init__(self, fault_id0, max_steps500): super().__init__() self.fault_id fault_id self.max_steps max_steps self.sim None self.step_count 0 # 动作空间3 个操纵变量归一化到 [-1, 1] self.action_space spaces.Box( low-1.0, high1.0, shape(3,), dtypenp.float32 ) # 观测空间选 10 个关键测量量 self.obs_vars [ reactor_temperature, reactor_pressure, reactor_level, separator_temperature, separator_level, stripper_temperature, stripper_level, product_flow, feed_flow, cooling_water_outlet_temp ] self.observation_space spaces.Box( low-np.inf, highnp.inf, shape(10,), dtypenp.float32 ) def reset(self): self.sim pytep.TEP(sampling_time0.01) self.sim.set_mode(1) self.sim.initialize() for _ in range(50): self.sim.step() if self.fault_id 0: self.sim.set_fault(self.fault_id) self.step_count 0 return self._get_obs() def _get_obs(self): obs self.sim.get_variables(self.obs_vars) return np.array(obs, dtypenp.float32) def step(self, action): # 把 [-1,1] 的动作映射到操纵变量的物理范围 # 这里简化处理实际要按变量逐个缩放 control (action 1.0) / 2.0 * 100.0 self.sim.set_variable(reactor_cooling_water_flow, control[0]) self.sim.set_variable(separator_cooling_water_flow, control[1]) self.sim.set_variable(stripper_steam_flow, control[2]) self.sim.step() self.step_count 1 obs self._get_obs() # 奖励设计负的温度偏差 惩罚项 temp_error abs(obs[0] - 120.0) reward -temp_error done self.step_count self.max_steps info {} return obs, reward, done, info逻辑说明reset()里重新实例化仿真对象保证每个 episode 独立。动作映射这里做了简化把三个动作统一缩放到 0-100实际项目里每个操纵变量的范围不同要分别处理。奖励函数只用了温度偏差做多目标控制时需要加权组合。done只按步数判断实际还要加安全约束——比如压力超限直接终止。参数说明obs_vars选了 10 个变量是为了控制观测维度全 41 维对 RL 来说太大。max_steps500对应 5 小时对 RL 训练来说 episode 偏长可以降到 200。动作空间的 shape 要和实际操纵变量数量一致改的时候别忘了同步改step()里的写入逻辑。4.3 交互式调试在 Jupyter 里逐步观察pyTEP 的交互性在调试时特别有用。你可以在 Jupyter 里一个 cell 推进一步实时看变量变化。# Jupyter cell 1: 初始化 import pytep sim pytep.TEP(sampling_time0.01) sim.set_mode(1) sim.initialize() print(initialized) # Jupyter cell 2: 推进一步并查看关键变量 sim.step() state sim.get_variables([ reactor_temperature, reactor_pressure, product_flow ]) print(temp{:.2f}, pressure{:.2f}, flow{:.2f}.format(*state)) # Jupyter cell 3: 注入故障后再看 sim.set_fault(1) for _ in range(10): sim.step() state sim.get_variables([ reactor_temperature, reactor_pressure, product_flow ]) print(after fault: temp{:.2f}, pressure{:.2f}, flow{:.2f}.format(*state))逻辑说明这种逐步交互适合验证你对故障动态的直觉。比如故障 1 是进料比例阶跃你应该看到产品流量和反应器温度在几十步内明显偏移。如果没变化要么故障没生效要么你观察的变量和该故障无关。Jupyter 里每次step()都会改变仿真状态所以 cell 的执行顺序很重要别跳着跑。5. 避坑与排查pyTEP 落地时最容易翻车的五个地方5.1 仿真发散但没有任何报错现象跑了几百步后变量变成 NaN 或极大值但step()不抛异常。原因TEP 是开环不稳定的如果你在初始化后没有及时接入控制器或者控制器参数太激进系统会自然发散。解决在initialize()后先跑一段开环稳态确认变量在合理范围接入控制器后监控关键变量一旦超限就重置。我一般会在循环里加一个检查发现异常立即break并打印当前状态。5.2 故障注入后动态和文献对不上现象你注入了故障 ID 4但反应器温度的响应曲线和论文里的图差很多。原因pyTEP 的故障定义可能和原始 TEP 有细微差异尤其是故障强度和作用位置。另外你的初始工况可能不是标准工况 1。解决先用故障 ID 1 做基准测试它的动态最明确确认set_mode的值和文献一致如果还是对不上检查采样周期是否被改过。5.3 批量生成数据时内存爆掉现象生成几千条轨迹后内存占用超过 32G。原因每条轨迹存了全部 41 个变量、500 个时间点float64 下一条就是 164KB一万条就是 1.6GB加上中间变量很容易爆。解决生成时只存你需要的变量用float32而不是默认的float64或者边生成边写入磁盘np.savez 或 hdf5。如果做在线训练用生成器逐条产出不要一次性全加载。5.4 控制器写入被忽略现象你调用了set_variable但下一步读回来发现操纵变量没变。原因变量名拼写错误或者该变量在当前模式下不可操纵。pyTEP 对未知变量名可能静默忽略而不是报错。解决先用list_variables()确认名称再检查该变量是否在操纵变量列表里。另一个常见原因是你在step()之后才调用set_variable那它要等下一步才生效看起来像被忽略了。5.5 多线程/多进程下仿真结果不可复现现象同样的种子单线程跑结果一致多进程跑每条轨迹都不一样。原因pyTEP 底层可能用了全局状态或非线程安全的随机数生成器。解决每个进程独立实例化TEP()不要共享如果 pyTEP 内部有随机性且不接受种子那只能接受轨迹间的随机差异或者改用单进程串行生成。做论文实验时我习惯把生成脚本和种子一起存档确保别人能复现。6. 进阶技巧用 pyTEP 做闭环性能基准测试前面讲的都是怎么把 pyTEP 跑起来、怎么生成数据、怎么接控制器。最后一章说一个更具体的用法把 pyTEP 当成控制器性能的基准平台用一套标准化的测试流程来对比不同算法的表现。这件事的价值在于你不需要自己搭化工模型也不用担心被审稿人质疑被控对象的合理性——TEP 本身就是公认基准。具体做法是定义一个测试套件包含若干标准场景工况 1 下的设定值跟踪、工况 1 下的故障 4 抑制、工况 3 下的进料扰动拒绝。每个场景跑固定时长记录关键变量的积分绝对误差和超调量。下面是一个简化的基准测试框架。import pytep import numpy as np def run_benchmark(controller_factory, scenario, n_steps500): 运行一个基准场景 controller_factory: 返回控制器对象的可调用对象 scenario: dict包含 mode, fault_id, setpoint sim pytep.TEP(sampling_time0.01) sim.set_mode(scenario[mode]) sim.initialize() # 稳态段 for _ in range(50): sim.step() # 注入故障 if scenario.get(fault_id, 0) 0: sim.set_fault(scenario[fault_id]) controller controller_factory() setpoint scenario[setpoint] errors [] for step in range(n_steps): current sim.get_variable(reactor_temperature) control controller.compute(setpoint, current) control np.clip(control, 0.0, 100.0) sim.set_variable(reactor_cooling_water_flow, control) sim.step() errors.append(abs(setpoint - current)) # 返回积分绝对误差和最大偏差 iae np.sum(errors) * 0.01 max_dev np.max(errors) return {IAE: iae, max_deviation: max_dev} # 定义两个控制器做对比 class PIDController: def __init__(self, Kp, Ki, Kd): self.Kp, self.Ki, self.Kd Kp, Ki, Kd self.integral 0.0 self.prev_error 0.0 def compute(self, setpoint, current): error setpoint - current self.integral error * 0.01 derivative (error - self.prev_error) / 0.01 self.prev_error error return self.Kp * error self.Ki * self.integral self.Kd * derivative # 场景定义 scenarios [ {name: setpoint_tracking, mode: 1, fault_id: 0, setpoint: 122.0}, {name: fault_rejection, mode: 1, fault_id: 4, setpoint: 120.0}, {name: feed_disturbance, mode: 3, fault_id: 0, setpoint: 120.0}, ] # 跑对比 for sc in scenarios: result_a run_benchmark(lambda: PIDController(2.0, 0.1, 0.5), sc) result_b run_benchmark(lambda: PIDController(1.0, 0.05, 0.2), sc) print(f{sc[name]}: PID-A IAE{result_a[IAE]:.2f}, PID-B IAE{result_b[IAE]:.2f})逻辑说明run_benchmark把场景配置和控制器解耦方便你换算法。controller_factory用可调用对象而不是实例是为了每个场景拿到干净的控制器状态。IAE 的计算乘了采样周期 0.01得到的是物理时间上的积分方便跨采样周期比较。max_deviation反映的是最坏情况对安全评估比 IAE 更敏感。参数说明PID 的两组参数是随手设的实际对比时应该先各自整定到最优否则比的是整定水平而不是算法本身。场景里的setpoint改成了 122.0 做跟踪测试因为 120.0 是稳态值跟踪测试需要给一个阶跃。mode3对应另一种进料工况它的稳态操纵变量值和模式 1 不同控制器需要重新整定。这个框架的扩展方向很直接把 PID 换成 MPC或者换成你正在做的 RL 策略跑同一套场景看 IAE 和 max_deviation 的对比。我自己的习惯是任何新控制器先在 setpoint_tracking 上跑通再上 fault_rejection最后才碰 feed_disturbance——因为前两个场景的失败模式最容易诊断第三个涉及工况切换变量耦合更复杂新手容易在这里迷失。一个血泪教训别在没跑通开环稳态的情况下直接上闭环。我有一次为了赶实验跳过稳态段直接注入故障加控制器结果系统在 50 步内就发散了排查了半天才发现是初始条件没收敛。后来我养成了一个习惯任何 pyTEP 脚本开头都先跑 50 步开环并打印关键变量确认在合理范围再往下走。这个检查花不了几秒钟但能省掉大量后悔药。希望帮到你。本文还有配套的精品资源点击获取
简介pyTEP 是面向过程控制、化工仿真与机器学习研究者的开源 Python 仿真 API用于田纳西伊士曼Tennessee Eastman流程的交互式模拟。它通过高级接口简化复杂工况的设置让用户无需深入理解模拟器底层原理即可调整操作条件、运行仿真并存储数据适合控制算法验证、故障诊断与数据驱动建模等场景。资源包共 76 个文件约 754KB包含 16 个 py 源码、23 个 m 文件、8 个 pkl 数据、7 个 mdl 模型及 mat、mexw64、mexa64 等混合文件覆盖 Python 接口、MATLAB 桥接、仿真引擎与测试模块并附 README、许可证与构建配置。目前已有 532 人学习下载。借助该包读者可快速搭建 TEP 仿真环境复现标准工况数据理解 siminterface、engineutils 等模块的调用方式为控制策略对比与异常检测实验提供可复用的工具基础。1. 用 pyTEP 把 Tennessee Eastman 流程装进 Python为什么它值得你花一个下午跑通如果你做过程控制、故障诊断或者强化学习在工业场景的落地大概率绕不开 Tennessee Eastman 流程。这个 1993 年就被提出的化工仿真基准至今仍是论文里的常客原因很简单它足够复杂又足够标准。但真正动手时你会发现很多现成的仿真环境要么绑死在商业软件里要么接口老旧、改一个参数要翻半天文档。pyTEP 这个开源 API 的出现把这件事拉回到了 Python 生态里——用几行代码就能建起一个可交互的闭环仿真改工况、加扰动、接控制器都在同一个脚本里完成。它适合三类人做数据驱动故障诊断、需要批量生成带标签工况数据的研究者做 MPC、RL 控制器验证、想要一个可控被控对象的工程师以及教学场景里需要让学生直观看到多变量耦合效应的讲师。这篇文章不讲空泛概念直接按「装什么、怎么建、参数怎么调、哪里会翻车」的顺序走一遍。2. 环境搭建与最小可运行闭环从零到第一条趋势曲线2.1 安装路径选择与依赖确认pyTEP 本质是对底层 Fortran 仿真核心的 Python 封装所以安装时最容易出问题的不是 Python 包本身而是编译工具链。常见做法有两种一是直接用 pip 从源码构建二是先拿到预编译的仿真核心再装 Python 层。我一般推荐后者因为省去 Fortran 编译器版本匹配的麻烦。# 建议先建独立环境避免和现有科学计算栈冲突 python -m venv pytep_env source pytep_env/bin/activate # Windows 下用 pytep_env\Scripts\activate # 升级构建工具老版本 setuptools 经常在编译扩展时翻车 pip install --upgrade pip setuptools wheel # 安装 pyTEP若本地已有仿真核心可加 --no-build-isolation pip install pytep装完后别急着写业务代码先做一次导入自检。这一步能过滤掉 80% 的「装上了但跑不起来」问题。# check_env.py import pytep import numpy as np # 打印版本和可用接口确认核心库被正确加载 print(pyTEP version:, getattr(pytep, __version__, unknown)) print(numpy version:, np.__version__) # 尝试实例化一个最简仿真对象不设任何扰动 sim pytep.TEP() print(simulation object created:, type(sim).__name__)逻辑说明这段代码只做三件事——确认包能导入、确认数值库版本兼容、确认仿真对象能实例化。参数上不需要改任何东西如果这一步报ImportError或OSError基本是编译产物没找到回去检查 Fortran 运行库路径。如果实例化时报参数错误说明 API 签名和文档有出入用help(pytep.TEP)看实际入参。2.2 构建第一个闭环仿真设定工况与采样pyTEP 的核心用法是「先定义仿真配置再逐步推进时间」。它不像某些环境那样一次性跑完返回数组而是给你一个可以逐拍交互的对象这对做在线控制非常关键。import pytep import numpy as np import matplotlib.pyplot as plt # 创建仿真实例指定采样周期为 0.01 小时约 36 秒 # 这个值对应 TEP 标准文献里的控制周期改大改小都会影响动态响应 sim pytep.TEP(sampling_time0.01) # 设定初始工况这里用标准工况 1即基本进料比例 sim.set_mode(1) # 初始化让系统进入稳态 sim.initialize() # 准备记录容器 time_log [] reactor_pressure_log [] # 推进 200 个采样点每步读取反应器压力 for step in range(200): sim.step() # 推进一个采样周期 time_log.append(step * 0.01) # 变量索引需查 pyTEP 变量表反应器压力通常是第 7 个测量量 reactor_pressure_log.append(sim.get_variable(reactor_pressure)) # 画趋势图确认仿真在正常波动而非发散 plt.plot(time_log, reactor_pressure_log) plt.xlabel(Time (h)) plt.ylabel(Reactor Pressure) plt.title(pyTEP Startup Trend) plt.grid(True) plt.show()逻辑说明set_mode(1)决定了进料和产物的基础配比TEP 一共定义了 6 种标准工况模式 1 是最常用的基准。initialize()会跑一段内部稳态计算不调用它直接step()可能拿到未收敛的初值。get_variable的字符串名需要和 pyTEP 的变量字典对齐不同版本命名可能有差异稳妥做法是先打印sim.list_variables()看一遍。参数说明sampling_time不建议随意改因为 TEP 的控制器参数和扰动设计都是基于 0.01 小时这个基准。如果你做强化学习需要更细的时间粒度正确做法是在step()之间做插值而不是改采样周期。step()每次推进一个周期返回的是当前状态字典还是 None 取决于版本别假设它有返回值。2.3 变量索引与测量量映射TEP 有 41 个测量变量和 12 个操纵变量新手最容易在这里迷路。pyTEP 一般提供名称到索引的映射但名称规范不一定和你的论文一致。# 打印所有可读变量名建立自己的映射表 variables sim.list_variables() for idx, name in enumerate(variables): print(f{idx:3d} {name}) # 常见做法把关键变量名存成常量避免代码里到处写魔法字符串 REACTOR_PRESSURE reactor_pressure REACTOR_LEVEL reactor_level SEPARATOR_LEVEL separator_level PRODUCT_FLOW product_flow # 批量读取多个变量减少函数调用开销 state sim.get_variables([REACTOR_PRESSURE, REACTOR_LEVEL, PRODUCT_FLOW]) print(current state:, state)逻辑说明list_variables()返回的顺序在不同版本里可能不同所以不要硬编码索引用名称查。get_variables批量读取比循环单次读取快做长时间仿真时这个差异会累积。如果你发现某个变量名不存在先确认仿真模式是否支持该变量——某些工况下部分测量量是冻结的。3. 扰动注入与故障数据生成把 TEP 变成你的数据工厂3.1 内置扰动类型与触发时机TEP 的价值一半在正常工况一半在它预定义的 20 多种故障。pyTEP 把这些故障做成了可编程接口你可以在任意时刻注入也可以设定持续时长。import pytep import numpy as np sim pytep.TEP(sampling_time0.01) sim.set_mode(1) sim.initialize() # 先跑 100 步让系统稳定 for _ in range(100): sim.step() # 在第 100 步注入故障 ID 4反应器冷却水入口温度阶跃 # 这个故障在文献里属于较难检测的类型适合测试诊断算法 sim.set_fault(4) # 再跑 300 步记录故障后的动态 fault_data [] for _ in range(300): sim.step() fault_data.append(sim.get_variables([ reactor_temperature, reactor_cooling_water_outlet_temp, product_flow ])) fault_data np.array(fault_data) print(fault data shape:, fault_data.shape)逻辑说明set_fault(4)是即时生效的不需要额外调用step()来激活。故障 ID 和 TEP 标准表对应但 pyTEP 可能对部分故障做了简化建议先用一个已知容易检测的故障比如 ID 1验证流程再换难检测的。fault_data的 shape 是(300, 3)可以直接喂给 sklearn 或 pytorch 做分类。参数说明故障注入的强度在 pyTEP 里通常是固定的对应原始 TEP 的设计值。如果你想做变强度故障常见做法是叠加一个缩放因子到测量量上而不是改故障本身——因为改故障定义会破坏和文献的可比性。3.2 批量生成带标签数据集做故障诊断模型训练你需要的是成千上万条带标签的轨迹而不是一条。pyTEP 支持重置和重新初始化但要注意重置的彻底性。import pytep import numpy as np def generate_trajectory(fault_id, n_steps500, seedNone): 生成一条带故障标签的轨迹 if seed is not None: np.random.seed(seed) sim pytep.TEP(sampling_time0.01) sim.set_mode(1) sim.initialize() # 稳态段 for _ in range(50): sim.step() # 注入故障fault_id 为 0 表示正常工况 if fault_id 0: sim.set_fault(fault_id) data [] for _ in range(n_steps): sim.step() data.append(sim.get_variables(sim.list_variables())) return np.array(data), fault_id # 批量生成正常 5 种故障每种 20 条 dataset [] labels [] for fid in [0, 1, 2, 4, 8, 13]: for run in range(20): traj, label generate_trajectory(fid, n_steps500, seedrun) dataset.append(traj) labels.append(label) dataset np.array(dataset) labels np.array(labels) print(dataset shape:, dataset.shape) # (120, 500, 41) print(labels shape:, labels.shape)逻辑说明每条轨迹都重新实例化TEP()这是最干净的做法避免状态残留。seed控制的是 numpy 的随机数pyTEP 内部如果有随机扰动需要确认它是否也接受种子——如果没有轨迹之间的差异只来自故障本身多样性会受限。list_variables()返回全部 41 个测量量实际建模时通常要做特征选择因为很多变量高度相关。参数说明n_steps500对应 5 小时仿真时长对慢动态故障足够对快动态故障可能太长。稳态段 50 步是经验值太短系统没稳太长浪费计算。如果你要做早期故障检测故障后的前 50 步才是关键窗口后面的数据反而会稀释信号。3.3 与 scikit-learn 对接的预处理原始 TEP 数据不能直接扔给模型标准化和滑窗是标配。from sklearn.preprocessing import StandardScaler from sklearn.model_selection import train_test_split # dataset shape: (n_samples, n_timesteps, n_features) n_samples, n_timesteps, n_features dataset.shape # 展平后做标准化再还原回三维 flat dataset.reshape(-1, n_features) scaler StandardScaler() flat_scaled scaler.fit_transform(flat) dataset_scaled flat_scaled.reshape(n_samples, n_timesteps, n_features) # 用滑窗构造样本窗口长度 20步长 5 window 20 stride 5 X, y [], [] for i in range(n_samples): for t in range(0, n_timesteps - window, stride): X.append(dataset_scaled[i, t:twindow, :]) y.append(labels[i]) X np.array(X) y np.array(y) print(X shape:, X.shape) # (n_windows, 20, 41) print(y shape:, y.shape) # 划分训练测试集注意按轨迹划分避免泄漏 X_train, X_test, y_train, y_test train_test_split( X, y, test_size0.3, random_state42, stratifyy )逻辑说明标准化必须用训练集的均值和方差但这里为了简化先在全量上 fit实际项目里要严格分开。滑窗的window20对应 0.2 小时对 TEP 的多数故障来说足够捕捉动态。stride5是为了减少窗口间的相关性如果做在线检测可以设 stride1。按轨迹划分比随机划分更严格因为同一条轨迹的窗口高度相似随机划分会导致测试集泄漏。4. 控制器接入与交互式仿真让 pyTEP 跑在你的算法里4.1 读取操纵变量并写入控制量pyTEP 的交互性体现在你可以每步读取状态、计算控制量、再写回去。这是做 MPC 或 RL 的标准循环。import pytep import numpy as np sim pytep.TEP(sampling_time0.01) sim.set_mode(1) sim.initialize() # 简单比例控制器控制反应器温度到设定值 setpoint 120.0 Kp 2.0 for step in range(500): # 读取当前反应器温度 current_temp sim.get_variable(reactor_temperature) # 计算控制量这里操纵冷却水流量 error setpoint - current_temp control Kp * error # 写入操纵变量注意要遵守物理上下限 control np.clip(control, 0.0, 100.0) sim.set_variable(reactor_cooling_water_flow, control) # 推进仿真 sim.step() if step % 50 0: print(fstep {step}: temp{current_temp:.2f}, control{control:.2f})逻辑说明set_variable必须在step()之前调用否则控制量会在下一步才生效造成一步延迟。np.clip是必须的因为 pyTEP 对超出物理范围的操纵变量可能直接报错或产生非物理行为。比例控制器的Kp2.0是随手设的实际整定需要看阶跃响应——TEP 的反应器温度回路增益不高Kp 太大容易振荡。参数说明setpoint120.0是反应器温度的典型设定值但不同工况下最优值不同。操纵变量名reactor_cooling_water_flow需要和list_variables()里的名称一致。如果你发现写入后温度没反应先确认这个变量在当前模式下是否可操纵——有些变量在特定工况下是固定的。4.2 与 Gym 接口对接做强化学习做 RL 的人希望 pyTEP 长得像 Gym 环境。自己包一层是最稳妥的因为现成的包装器往往版本滞后。import gym from gym import spaces import pytep import numpy as np class TEPEnv(gym.Env): 把 pyTEP 包装成 Gym 环境 def __init__(self, fault_id0, max_steps500): super().__init__() self.fault_id fault_id self.max_steps max_steps self.sim None self.step_count 0 # 动作空间3 个操纵变量归一化到 [-1, 1] self.action_space spaces.Box( low-1.0, high1.0, shape(3,), dtypenp.float32 ) # 观测空间选 10 个关键测量量 self.obs_vars [ reactor_temperature, reactor_pressure, reactor_level, separator_temperature, separator_level, stripper_temperature, stripper_level, product_flow, feed_flow, cooling_water_outlet_temp ] self.observation_space spaces.Box( low-np.inf, highnp.inf, shape(10,), dtypenp.float32 ) def reset(self): self.sim pytep.TEP(sampling_time0.01) self.sim.set_mode(1) self.sim.initialize() for _ in range(50): self.sim.step() if self.fault_id 0: self.sim.set_fault(self.fault_id) self.step_count 0 return self._get_obs() def _get_obs(self): obs self.sim.get_variables(self.obs_vars) return np.array(obs, dtypenp.float32) def step(self, action): # 把 [-1,1] 的动作映射到操纵变量的物理范围 # 这里简化处理实际要按变量逐个缩放 control (action 1.0) / 2.0 * 100.0 self.sim.set_variable(reactor_cooling_water_flow, control[0]) self.sim.set_variable(separator_cooling_water_flow, control[1]) self.sim.set_variable(stripper_steam_flow, control[2]) self.sim.step() self.step_count 1 obs self._get_obs() # 奖励设计负的温度偏差 惩罚项 temp_error abs(obs[0] - 120.0) reward -temp_error done self.step_count self.max_steps info {} return obs, reward, done, info逻辑说明reset()里重新实例化仿真对象保证每个 episode 独立。动作映射这里做了简化把三个动作统一缩放到 0-100实际项目里每个操纵变量的范围不同要分别处理。奖励函数只用了温度偏差做多目标控制时需要加权组合。done只按步数判断实际还要加安全约束——比如压力超限直接终止。参数说明obs_vars选了 10 个变量是为了控制观测维度全 41 维对 RL 来说太大。max_steps500对应 5 小时对 RL 训练来说 episode 偏长可以降到 200。动作空间的 shape 要和实际操纵变量数量一致改的时候别忘了同步改step()里的写入逻辑。4.3 交互式调试在 Jupyter 里逐步观察pyTEP 的交互性在调试时特别有用。你可以在 Jupyter 里一个 cell 推进一步实时看变量变化。# Jupyter cell 1: 初始化 import pytep sim pytep.TEP(sampling_time0.01) sim.set_mode(1) sim.initialize() print(initialized) # Jupyter cell 2: 推进一步并查看关键变量 sim.step() state sim.get_variables([ reactor_temperature, reactor_pressure, product_flow ]) print(temp{:.2f}, pressure{:.2f}, flow{:.2f}.format(*state)) # Jupyter cell 3: 注入故障后再看 sim.set_fault(1) for _ in range(10): sim.step() state sim.get_variables([ reactor_temperature, reactor_pressure, product_flow ]) print(after fault: temp{:.2f}, pressure{:.2f}, flow{:.2f}.format(*state))逻辑说明这种逐步交互适合验证你对故障动态的直觉。比如故障 1 是进料比例阶跃你应该看到产品流量和反应器温度在几十步内明显偏移。如果没变化要么故障没生效要么你观察的变量和该故障无关。Jupyter 里每次step()都会改变仿真状态所以 cell 的执行顺序很重要别跳着跑。5. 避坑与排查pyTEP 落地时最容易翻车的五个地方5.1 仿真发散但没有任何报错现象跑了几百步后变量变成 NaN 或极大值但step()不抛异常。原因TEP 是开环不稳定的如果你在初始化后没有及时接入控制器或者控制器参数太激进系统会自然发散。解决在initialize()后先跑一段开环稳态确认变量在合理范围接入控制器后监控关键变量一旦超限就重置。我一般会在循环里加一个检查发现异常立即break并打印当前状态。5.2 故障注入后动态和文献对不上现象你注入了故障 ID 4但反应器温度的响应曲线和论文里的图差很多。原因pyTEP 的故障定义可能和原始 TEP 有细微差异尤其是故障强度和作用位置。另外你的初始工况可能不是标准工况 1。解决先用故障 ID 1 做基准测试它的动态最明确确认set_mode的值和文献一致如果还是对不上检查采样周期是否被改过。5.3 批量生成数据时内存爆掉现象生成几千条轨迹后内存占用超过 32G。原因每条轨迹存了全部 41 个变量、500 个时间点float64 下一条就是 164KB一万条就是 1.6GB加上中间变量很容易爆。解决生成时只存你需要的变量用float32而不是默认的float64或者边生成边写入磁盘np.savez 或 hdf5。如果做在线训练用生成器逐条产出不要一次性全加载。5.4 控制器写入被忽略现象你调用了set_variable但下一步读回来发现操纵变量没变。原因变量名拼写错误或者该变量在当前模式下不可操纵。pyTEP 对未知变量名可能静默忽略而不是报错。解决先用list_variables()确认名称再检查该变量是否在操纵变量列表里。另一个常见原因是你在step()之后才调用set_variable那它要等下一步才生效看起来像被忽略了。5.5 多线程/多进程下仿真结果不可复现现象同样的种子单线程跑结果一致多进程跑每条轨迹都不一样。原因pyTEP 底层可能用了全局状态或非线程安全的随机数生成器。解决每个进程独立实例化TEP()不要共享如果 pyTEP 内部有随机性且不接受种子那只能接受轨迹间的随机差异或者改用单进程串行生成。做论文实验时我习惯把生成脚本和种子一起存档确保别人能复现。6. 进阶技巧用 pyTEP 做闭环性能基准测试前面讲的都是怎么把 pyTEP 跑起来、怎么生成数据、怎么接控制器。最后一章说一个更具体的用法把 pyTEP 当成控制器性能的基准平台用一套标准化的测试流程来对比不同算法的表现。这件事的价值在于你不需要自己搭化工模型也不用担心被审稿人质疑被控对象的合理性——TEP 本身就是公认基准。具体做法是定义一个测试套件包含若干标准场景工况 1 下的设定值跟踪、工况 1 下的故障 4 抑制、工况 3 下的进料扰动拒绝。每个场景跑固定时长记录关键变量的积分绝对误差和超调量。下面是一个简化的基准测试框架。import pytep import numpy as np def run_benchmark(controller_factory, scenario, n_steps500): 运行一个基准场景 controller_factory: 返回控制器对象的可调用对象 scenario: dict包含 mode, fault_id, setpoint sim pytep.TEP(sampling_time0.01) sim.set_mode(scenario[mode]) sim.initialize() # 稳态段 for _ in range(50): sim.step() # 注入故障 if scenario.get(fault_id, 0) 0: sim.set_fault(scenario[fault_id]) controller controller_factory() setpoint scenario[setpoint] errors [] for step in range(n_steps): current sim.get_variable(reactor_temperature) control controller.compute(setpoint, current) control np.clip(control, 0.0, 100.0) sim.set_variable(reactor_cooling_water_flow, control) sim.step() errors.append(abs(setpoint - current)) # 返回积分绝对误差和最大偏差 iae np.sum(errors) * 0.01 max_dev np.max(errors) return {IAE: iae, max_deviation: max_dev} # 定义两个控制器做对比 class PIDController: def __init__(self, Kp, Ki, Kd): self.Kp, self.Ki, self.Kd Kp, Ki, Kd self.integral 0.0 self.prev_error 0.0 def compute(self, setpoint, current): error setpoint - current self.integral error * 0.01 derivative (error - self.prev_error) / 0.01 self.prev_error error return self.Kp * error self.Ki * self.integral self.Kd * derivative # 场景定义 scenarios [ {name: setpoint_tracking, mode: 1, fault_id: 0, setpoint: 122.0}, {name: fault_rejection, mode: 1, fault_id: 4, setpoint: 120.0}, {name: feed_disturbance, mode: 3, fault_id: 0, setpoint: 120.0}, ] # 跑对比 for sc in scenarios: result_a run_benchmark(lambda: PIDController(2.0, 0.1, 0.5), sc) result_b run_benchmark(lambda: PIDController(1.0, 0.05, 0.2), sc) print(f{sc[name]}: PID-A IAE{result_a[IAE]:.2f}, PID-B IAE{result_b[IAE]:.2f})逻辑说明run_benchmark把场景配置和控制器解耦方便你换算法。controller_factory用可调用对象而不是实例是为了每个场景拿到干净的控制器状态。IAE 的计算乘了采样周期 0.01得到的是物理时间上的积分方便跨采样周期比较。max_deviation反映的是最坏情况对安全评估比 IAE 更敏感。参数说明PID 的两组参数是随手设的实际对比时应该先各自整定到最优否则比的是整定水平而不是算法本身。场景里的setpoint改成了 122.0 做跟踪测试因为 120.0 是稳态值跟踪测试需要给一个阶跃。mode3对应另一种进料工况它的稳态操纵变量值和模式 1 不同控制器需要重新整定。这个框架的扩展方向很直接把 PID 换成 MPC或者换成你正在做的 RL 策略跑同一套场景看 IAE 和 max_deviation 的对比。我自己的习惯是任何新控制器先在 setpoint_tracking 上跑通再上 fault_rejection最后才碰 feed_disturbance——因为前两个场景的失败模式最容易诊断第三个涉及工况切换变量耦合更复杂新手容易在这里迷失。一个血泪教训别在没跑通开环稳态的情况下直接上闭环。我有一次为了赶实验跳过稳态段直接注入故障加控制器结果系统在 50 步内就发散了排查了半天才发现是初始条件没收敛。后来我养成了一个习惯任何 pyTEP 脚本开头都先跑 50 步开环并打印关键变量确认在合理范围再往下走。这个检查花不了几秒钟但能省掉大量后悔药。希望帮到你。本文还有配套的精品资源点击获取
