搞污水处理的控制、节能优化，或者工艺动态仿真的朋友，十有八九都会撞上一个叫BSM1的模型。它的全称是Benchmark Simulation Model No.1，中文一般叫“一号基准仿真模型”或者BSM1基准模型，由国际水协会（IWA）主导开发，专门用来给活性污泥法污水处理工艺的动态仿真和自动控制研究提供一个相对标准、可重复、能横向比较的实验平台。
我第一次真正动手搭BSM1，是在做曝气系统节能优化的时候。当时我满怀信心地把ASM1的微分方程抄进Simulink，结果第一天仿真就发散，满屏NaN。折腾了一周才意识到：BSM1表面上是一个“模型”，实际上是一整套规范，包括进水数据、工艺拓扑、模型参数、评价指标，甚至仿真流程都有明确规定，少一个环节，结果就不可比。这篇文章把我从零搭起BSM1的经验完整写下来，内容涉及模型结构、ASM1方程细节、Simulink搭建步骤、求解器配置、性能指标计算，以及我实际踩过的坑，希望能给准备入门或者正在被BSM1折磨的同行一点帮助。
1. BSM1到底解决了什么问题
1.1 为什么研究污水控制需要“统一考题”
活性污泥法是个典型的大时滞、强耦合、多变量非线性系统。进水流量和水质时刻在变，微生物种群之间相互影响，溶解氧、氨氮、硝酸盐之间又存在明显的耦合关系。早些年大家做控制策略研究，基本都是“各玩各的”：有人用自己的厂站数据，有人用自建模拟装置，有人拿某个实际项目的负荷曲线做输入，结果同样是“提升了控制效果”，却根本没法放在一起比较，审稿人也很难判断你的改进是算法真的有效，还是单纯因为你的场景更简单。
所以国际水协会的专家们在20世纪90年代开始推动一件事：设计一套统一的污水处理仿真基准模型，规定好工艺结构、进水数据、模型参数和评价规则，让全世界的研究者在同一份“考卷”上做题。于是就有了BSM1（Benchmark Simulation Model No.1）。它是目前国际上研究活性污泥工艺控制、能耗优化、联表参数辨识、数据驱动控制等方向使用最广的参考模型之一。对于做算法研究的人来说，它是标准测试台；对于工艺工程师来说，它是理解污水厂动态响应的极好学习工具。
如果你刚接触这个领域，我建议把BSM1理解为“污水处理控制领域的MNIST数据集”或者“自动驾驶里的CARLA仿真平台”。它不是某个真实污水厂的精确数字化映射，而是一个足够复杂、又足够规范的虚拟被控对象，目的是让控制算法可以在去除“数据隐私、设备差异、天气不确定性”等干扰因素的前提下，公平、稳定地比较优劣。
1.2 BSM1的工艺结构
BSM1的工艺拓扑可以概括为“5个生物反应器+1个二沉池”，生物反应器部分分前段缺氧区和后段好氧区。官方定义里，反应器1和反应器2是缺氧池，体积各为1000 m³；反应器3、4、5是好氧池，体积各为1333 m³，生物反应器总体积6000 m³。二沉池是标准的10层分层沉降模型，表面积1500 m²，深度4 m，总体积6000 m³。
水流关系上，原水进入反应器1，依次流经5个反应器后进入二沉池。二沉池底部有两条重要的回流管线：一条是外回流（污泥回流），把底部浓缩后的活性污泥送回反应器1，维持生物系统内的污泥浓度；另一条是剩余污泥排放，用于控制泥龄。除了外回流，BSM1还设置了一路内回流，把反应器5（好氧末端）的混合液部分送回反应器1，以提供缺氧池反硝化所需的硝酸盐，这一路相对于外回流，流量通常更大。
整套系统只有当所有回流、所有反应器体积、所有参数都对齐后，才能复现出BSM1标准工况下“出水氨氮、总氮、COD在一定范围内波动”的典型响应。很多新手最容易犯的错就是少连了一根回流管线，或者回流流量方向接反，结果仿真也能“跑起来”，但出水的动态曲线跟官方基准结果完全对不上。
1.3 为什么用ASM1而不是其他模型
BSM1的生化反应部分采用活性污泥1号模型，也就是ASM1（Activated Sludge Model No.1）。这个选择有它的历史必然性。ASM1是IWA在1987年发布的第一个国际公认的活性污泥模型，核心是用“Monod动力学+开关函数”描述有机底物的降解、硝化和反硝化过程。相比后来的ASM2（增加了生物除磷）、ASM3（修正了储存代谢机制），ASM1结构相对简单，参数数量适中，但已经足以体现脱氮除碳过程中的主要动态特性。
从实用角度讲，ASM1包含13个状态变量和多个动力学反应过程，既能表现出进水冲击下的动态响应，又不至于复杂到让控制算法仿真慢到无法迭代。BSM1把它作为默认生化模型，正是看中了这种“复杂度与可解释性之间的平衡”。你自己做研究时，也可以基于ASM1改扩展模型，但前提是要先跑通原始的BSM1基线，否则后续策略改进的对比没有说服力。
2. 在Simulink里动手前必须搞懂的隐藏逻辑
2.1 ASM1的13个状态变量怎么看
很多人拿到BSM1的文档，第一反应是被ASM1那一大张化学计量矩阵吓住。其实不用慌，你只需要把13个状态变量当成13个“浓度状态”即可。官方文档和绝大多数仿真代码里，状态变量的顺序是固定的，常见顺序如下：
| 序号 | 符号 | 含义 | 单位 |
|---|---|---|---|
| 1 | S_I | 惰性可溶性有机物 | mg COD/L |
| 2 | S_S | 易生物降解底物 | mg COD/L |
| 3 | X_I | 惰性颗粒态有机物 | mg COD/L |
| 4 | X_S | 慢速可生物降解底物 | mg COD/L |
| 5 | X_BH | 异养菌生物量 | mg COD/L |
| 6 | X_BA | 自养菌生物量 | mg COD/L |
| 7 | X_P | 微生物内源呼吸产生的惰性颗粒物 | mg COD/L |
| 8 | S_O | 溶解氧 | mg COD/L |
| 9 | S_NO | 硝酸盐氮 | mg N/L |
| 10 | S_NH | 氨氮 | mg N/L |
| 11 | S_ND | 可生物降解可溶性有机氮 | mg N/L |
| 12 | X_ND | 可生物降解颗粒态有机氮 | mg N/L |
| 13 | S_ALK | 碱度 | mol HCO3/m³ |
注意，第8项溶解氧的单位用的是“负COD当量”，这是因为ASM1里氧的化学计量系数按COD平衡推导，这容易让新手困惑：反应器里溶解氧明明是mg/L级别的数值，怎么代进模型后符号却要参与COD衡算。实际写代码时，最后输出给控制器的溶解氧浓度仍按mg/L处理，但状态更新过程中必须遵守ASM1的矩阵规则。
2.2 二沉池不是简单分离器
二沉池在BSM1里绝不是简单地把悬浮固体“按比例一分为二”，而是一个分层沉降模型。官方采用Takács双指数沉降速率方程来模拟活性污泥的沉降特性，并根据固体通量理论计算每一层中的颗粒物浓度变化。
Takács方程的核心公式可以写成：
vt = v0 * (exp(-rh * X) - exp(-rp * X))
其中v0是最大沉降速度，rh和rp是沉降参数，X是该层悬浮固体浓度。这个方程描述的是：低浓度时沉降速度会因“澄清区”的阻滞而下降，高浓度时又会因“压缩区”的拥挤而下降，中间某个浓度时沉降速度最大。二沉池被垂直划分为10层，每一层按固体通量平衡进行更新，进水从中间某一层进入，上层清水从顶部流出，底层污泥从底部排出并回流。
实际搭建时，二沉池模块的输出会直接影响外回流污泥浓度和内回流硝酸盐浓度。如果简化成固定比例分配，你会发现系统很难稳定，或者出水SS曲线跟官方结果差得很远。我的建议是：初学阶段不要重新发明轮子，直接用官方固化的沉降模型公式，至少先保住结果可对比。
2.3 回流、曝气和初始条件的工程设定
BSM1的标准工况有一组默认的工艺参数。内回流量Qa通常设为进水流量Q0的倍数，外回流量Qr根据污泥回流比设定，剩余污泥排放量Qw按泥龄控制。这里每个回流量的单位都是m³/d，时间基准是天，Simulink里所有积分器的时间单位也统一用“天”，否则流量和浓度的时间导数很容易差出几十倍。
溶解氧控制方面，BSM1标准结构里好氧池的曝气通过氧气传质系数KLa来调节。你可以把KLa设成固定常数模拟“有阀门的恒定曝气”，也可以加一个PI控制器，根据好氧池溶解氧测量值与设定值的偏差实时调整KLa。后者是很多控制文章的起点。
初始条件方面，BSM1有一个约定俗成的做法：先在稳态进水条件下运行足够长时间，让系统收敛到“动态周期的起点”，再切换成14天动态进水开始正式仿真。如果没有先跑稳态，直接拿零浓度初值做动态仿真，大概率前几个时刻就出现负浓度或振荡发散。
3. 在Matlab/Simulink中逐步搭建BSM1
3.1 准备阶段：数据文件和参数表
在打开Simulink画框图之前，先把三样东西准备好：ASM1模型方程、进水数据文件、动力学参数表。官方BSM1的14天进水数据包括干天、雨天、暴风雨三种典型工况，默认最常用的就是干天数据。文件通常是CSV或MAT格式，每列分别对应流量和各进水组分浓度，采样间隔为15分钟，总共1345个左右的数据点。
数据导入Simulink时，我建议用readtable读入Excel或CSV，再用timeseries封装成连续信号，或者直接用From Workspace模块按时间序列读取。这里有个细节：如果进水CSV里的时间单位是小时，而你的Simulink求解器用的是天，那流量数值不会错，但“时间轴”一定要做换算。进水流量的峰值、负荷冲击到达反应器的时间都会受此影响。
动力学参数部分，BSM1官方参数表定义了大量化学计量系数和动力学常数，包括异养菌最大比增长速率、自养菌最大比增长速率、半饱和常数、衰减系数等。推荐的做法是把这些参数统一放进一个MATLAB结构体里，例如param.mu_H = 4.0，然后在每个S-Function里直接引用同一个结构体。不要分散写在各个模块里，不然后面做灵敏度分析或者参数校准时会改到怀疑人生。
3.2 用Level-2 S-Function写反应器
BSM1的5个生物反应器可以用一个通用的Level-2 MATLAB S-Function实现。每个反应器本质上是基于质量平衡的常微分方程组：
V * dC/dt = Q_in * C_in - Q_out * C + V * r(C)
其中C是13维状态向量，Q_in和C_in来自上一级反应器或进水与回流混合液，V是该反应器体积，r(C)是ASM1生化反应速率向量。
下面给一个S-Function的骨架，方便理解结构：
function BSM1Reactor(block) setup(block); block.RegistrationMethod('PostProcess', @PostProcess); end function setup(block) block.NumInputPorts = 1; % 输入: 上游浓度向量(13维) block.NumOutputPorts = 1; % 输出: 本反应器浓度(13维) block.SetPreCompInpPortInfoToDynamic; block.SetPreCompOutPortInfoToDynamic; block.InputPort(1).Dimensions = 13; block.OutputPort(1).Dimensions = 13; block.NumContStates = 13; block.SampleTimes = [0 0]; % 连续时间 block.SetAllowConstantSampleTime(true); end function derivatives(block) x = block.ContStates.Data; u = block.InputPort(1).Data; % 根据工艺参数计算Qin, Qout以及进水混合后的Cin dxdt = ASM1_reactor_ode(x, u, Q_in, Q_out, V, param); block.Derivatives.Data = dxdt; end实际使用中，流量Q_in不能简单等于上一级出流，还要叠加外回流和内回流混入的部分。所以更合适的做法是把每个反应器的输入分成“来自上级的跨流”和“来自回流的混合流”，在连接时用混合模块先算好混合后的浓度和流量，再送给反应器。
ASM1的动力学矩阵比较复杂，建议单独写一个函数文件，输入当前状态和参数，输出反应速率向量。注意所有组分单位要保持一致，比如有机物的单位是mg COD/L，氮组分的单位是mg N/L，碱度是mol/m³，你不能把它们混在一个向量里不做区分。
3.3 搭建二沉池与回流管线
二沉池模块是整个BSM1里最容易出错的部分。标准做法是把它写成一个单独的S-Function，输入为二沉池入流的流量和浓度（实际上是反应器5出流），输出为顶层出水浓度和底层回流污泥浓度，中间通过分层计算得到每层的固体浓度分布。
由于二沉池内部对溶解性组分和颗粒性组分的处理方式不同，你需要把13个状态变量拆成两类：溶解性组分随水流走，颗粒性组分额外受沉降过程影响。这也是为什么很多简化模型在二沉池上“翻车”的原因——只按总SS沉降，却忽略了不同组分在澄清液和回流污泥中的再分配。
搭建回流管线时，可以从二沉池底部输出直接通过一个增益模块或流量拆分模块分出一路作为外回流，另一路作为剩余污泥排放。内回流则是从反应器5取混合液送回反应器1，流量大、浓度接近好氧末端状态。接线时要特别注意代数环问题：如果回流模块直接读当前时刻二沉池输出再反馈到反应器1，但二沉池输入又来自反应器5，Simulink可能生成代数环，导致求解器报错或仿真极慢。可以在回流管线上加Memory或Unit Delay来打破代数环，但要注意引入的人为延迟足够小，不影响动态特性。
3.4 求解器怎么设置不让仿真发散
BSM1属于典型的刚性系统：溶解氧和碱度等组分的动态响应快，而微生物生物量的动态响应慢，不同状态变量的时间尺度差了好几个数量级。如果采用普通变步长ode45，仿真往往需要极长的时间，还可能在某个时刻突然不稳定。所以一定要选择刚性求解器，Simulink里推荐ode15s。
我的经验参数是：相对容差RelTol设为1e-4到1e-6，绝对容差AbsTol按各组分数量级分别设置，最大步长MaxStep设为0.01天（约14.4分钟），这样既能捕捉进水15分钟采样间隔引起的动态冲击，又不会让求解器在快变量上浪费时间。如果你发现仿真初期出现NaN或负浓度，优先检查初始化条件和限幅设置，而不是盲目调小步长。更稳妥的做法是在S-Function的导数函数里对所有状态做max(x, 0)限幅，但这只能作为兜底，真正的问题可能出在回流浓度或者动力参数单位上。
4. 跑完仿真怎么算指标：EQI、OCI与合规判定
4.1 出水水质指数EQI的计算逻辑
BSM1的价值不在于“仿真曲线好看”，而在于它可以量化比较不同控制策略的好坏。官方定义了出水水质指数EQI（Effluent Quality Index）来统一评价出水质量。EQI不是简单的某个污染物浓度，而是把多个关键出水指标按权重加权后，在整个仿真周期内积分：
EQI = (1 / (T * 1000)) * ∫(2SS_e + COD_e + 30TKN_e + 10NO_e + 2BOD5_e) * Q_e * dt
这里的SS_e是出水悬浮固体，COD_e是出水化学需氧量，TKN_e是出水总凯氏氮，NO_e是出水硝态氮，BOD5_e是出水五日生化需氧量，Q_e是出水流量。权重系数反映了不同污染物对环境影响的相对重要程度，其中氮类物质权重明显更高，因为脱氮是BSM1这种脱氮除碳工艺的核心目标。
计算EQI时，建议不要用Simulink里的积分器直接读，因为步长变化可能导致积分误差。更稳妥的做法是把仿真结果保存到工作区，再用trapz函数做数值积分，脚本如下：
EQI = trapz(t, (2*SS + COD + 30*TKN + 10*NO + 2*BOD5) .* Qe) / (t(end) - t(0)) / 1000;注意单位统一，流量用m³/d，浓度用mg/L，时间用d，计算出来的EQI单位才是kg/d级别。
4.2 运行成本指数OCI的构成
控制策略不能只看水质，还要看经济性，于是BSM1又定义了运行成本指数OCI（Overall Cost Index）。OCI综合考虑了曝气能耗、泵送能耗、污泥处理成本等。曝气能耗通常与KLa和好氧池体积、溶解氧水平相关；泵送能耗主要来自内回流和外回流的流量；污泥处理成本跟剩余污泥排放量直接相关。
一个典型公式大致是：
OCI = AE + PE + SP + CE - ME
其中AE是曝气能耗，PE是泵送能耗，SP是污泥处理成本，CE是外加碳源成本，ME是混合能耗。各部分都有基于物理量的折算系数。比如曝气能耗往往写成KLa的函数，泵送能耗跟Qa和Qr成正比，剩余污泥处理量则跟Qw和污泥浓度有关。如果你用了外加碳源投加控制策略，CE也要算进去。BSM1官方文档给出了完整的系数表，直接按公式套即可，不需要自己推导。
很多初学者只关注水质指标，忽略OCI，结果论文里说“控制效果变好”，但一看能耗也涨了一大截。在控制策略研究里，EQI和OCI往往是一对矛盾目标，类似投资里的收益和风险，很多文章的核心贡献就是通过优化算法在两者之间找Pareto前沿。
4.3 出水合规判定
除了连续积分指标，BSM1还设置了出水超标判定。官方限值大概包括：出水COD不超过100 mg/L，总氮不超过18 mg/L，氨氮不超过4 mg/L，SS不超过30 mg/L。仿真结束后统计每个指标的超标次数或超标持续时间，整理成表格，就能直观看出某个控制策略在“达标率”上的表现。
这个步骤看似简单，但有很多细节：超标判定要基于出水浓度，而BSM1中出水浓度是二沉池顶层出口的浓度，不是反应器5的浓度；统计周期通常是第7天到第14天（排除前几天的初始效应）；不同论文对“超标”的定义略有差异，有的看瞬时值，有的看日平均值。写论文时务必交代清楚，否则别人无法复现你的合规统计结果。
5. 常见问题排查与扩展玩法
5.1 仿真NaN、负浓度、发散怎么查
我把这些年遇到过的BSM1仿真问题整理成一张速查表，按出现频率排序：
| 现象 | 可能原因 | 处理建议 |
|---|---|---|
| 一开始就NaN | 初始浓度字段没初始化或初值离谱 | 使用官方初始状态文件，不要从零开始 |
| 前几秒就负浓度 | 回流接线方向不对，形成正反馈 | 检查内回流、外回流是否指向正确反应器 |
| 特定组分长期为0 | 进水数据列顺序与状态变量顺序不对应 | 核对ASM1组分顺序，确认S_S、X_BH等列没串位 |
| 仿真步长极小、进展缓慢 | 求解器不适合刚性系统 | 换ode15s，并适当放宽RelTol到1e-4 |
| 结果与官方曲线差异大 | 二沉池模型简化为比例分离 | 改用Takács分层沉降模型 |
| 控制器不稳定 | 溶解氧回路采样时间不合理 | PI控制器采样时间不要远大于0.01天 |
我的习惯是：遇到发散先不猜，直接在几个关键时间点打印反应器1和反应器5的状态向量，观察是哪一个组分先出问题。如果是S_NH或S_O先变成负值，多半是初始化或限幅问题；如果是X_BH暴涨，那大概率是回流流量或底物浓度单位错了。
5.2 初始稳态怎么建立
BSM1的正式仿真流程要求在动态仿真前先建立初始稳态。最简单的做法是让模型在“恒定进水工况”下运行足够长的时间，直到所有状态变量不再明显变化，再把最终状态作为动态仿真的初值。实际操作中，官方提供了初始状态的参考值文件，你可以直接把这些值写进Integrator模块的Initial Condition，省去跑稳态的过程。
但这里有个坑：官方初始状态是基于特定参数的，如果你自己改了动力学常数或工艺条件，直接套用官方初始值可能导致初值不匹配，动态仿真前几个周期出现漂移。所以更严谨的做法是：用你修改后的参数重新跑一段较长稳态仿真，确认状态稳定后再取初值。动态仿真正式评价时，一般从第7天开始统计，让前7天作为“预热期”。
5.3 在BSM1上做控制策略的常见玩法
BSM1最经典的控制场景是两条回路：一是好氧池溶解氧控制，通常通过调节曝气量或KLa实现；二是缺氧池硝酸盐控制，通常通过调节内回流量或外加碳源实现。在此基础上，你可以尝试PID参数整定、前馈补偿、模型预测控制（MPC）、模糊控制，甚至强化学习。
我在实际训练强化学习智能体时发现一个重要的技巧：如果直接拿BSM1的Simulink模型做每一步的训练交互，仿真速度会非常慢，因为每个控制步长都要调用一次刚性求解器。更经济的做法是先用BSM1生成数据集，训练一个简化RNN或神经状态空间模型作为代理模型，再让强化学习算法在代理模型上试错，最后回到BSM1里做验证。这本质上是现在很主流的“仿真到现实”思路，可以大幅压缩研究周期。
另外，如果你的目标是做嵌入式在线控制而不是纯仿真研究，可以考虑用Simulink Coder把BSM1里的控制算法模型生成C代码，部署到PLC或实时控制器上。在硬件在环测试中，BSM1可以当作被控对象的实时仿真器，控制算法跑在真实控制器里，这个组合在工业界和学术界都有不少应用。
5.4 从BSM1到BSM2的平滑过渡
如果你觉得BSM1不够过瘾，或者研究的系统不止于“水线”，BSM2是自然的下一步。BSM2在BSM1基础上加入了污泥处理单元，包括厌氧消化池、贮泥池、脱水单元等，形成更完整的污水处理厂全流程仿真平台。BSM2的进水数据更长、系统时间常数更大，模型刚性也更明显，求解器设置需要额外注意。
从学习方法看，BSM1是把活性污泥工艺动态仿真搞清楚的最佳起点。BSM2则更接近“整个厂级系统”的研究，适合做全厂能量优化、污泥资源化等方向。我的建议是：先把BSM1的每个模块、每条回流、每个指标彻底吃透，再扩展BSM2，你会发现两者在建模思想和评价框架上一脉相承，迁移成本不高。
最后说点个人体会。BSM1之所以能在污水处理仿真领域“活”这么多年，不是因为它的模型有多精确，而是因为它的基准属性：同一份进水、同一套参数、同样的评价公式，让所有人的结果都可以被严格复现和比较。这个特性在研究里太宝贵了。对新手来说，与其直接找一个大而全的污水厂模型，不如从BSM1入手，亲手搭一遍、跑一遍、算一遍指标，你会发现以前那些“工艺参数对出水影响”的模糊认知，会突然变得非常具体。我做仿真这些年，一直保留着BSM1的原始Simulink模型，偶尔还会拿它测试新的控制思想。数据可以换、参数可以调，但那份“统一考题”的框架，始终是入行最值得先学会的东西。
