A3：科技智慧
引 言
城市排水系统的实时控制（Real-time control, RTC）方法是指通过对泵、闸等城市排水系统中的流量控制设备进行实时操作，对降雨及管网内部的水力学特征变化做出及时的响应，以达到“蓄排结合”的目的，从而尽可能避免合流制溢流及城市内涝的发生。模型预测控制（Model Predictive Control, MPC）是一种常用的预测式实时控制方法,其基本结构与步骤如图1所示。与响应式控制（Reactive control）不同，MPC在每一步控制中都需要通过模型模拟的方法对未来影响进行预测和分析，以得到效果更好的控制方案，并逐时间步滚动进行控制的实时操作。
图1 模型预测控制基本结构及步骤
排水系统的实时控制需要对水力、水文学状态进行实时的数据分析和预测，因此在模拟环境中可以采用pyswmm的逐时间步模拟方法（Simulation模块）实现MPC的实时操作。然而，pyswmm不是一个线程安全的模拟方式，因此，无法在逐时间步实时模拟的同时，进行某一时间步控制方案的预测评估（Sadler et al., 2019）。因此需要用EPA-SWMM的计算引擎，加载此时pyswmm提供的系统状态信息（水文、水力变量）作为初始边界条件进行热启动，对未来一段时间的控制效果进行预测模拟，实现系统的实时预测。MPC的实时模拟和实时预测主要流程如下图所示。
图2 基于pyswmm和EPA-SWMM的排水系统MPC实时模拟和实时预测
一、实时模拟
尽管EPA-SWMM提供了控制模块，可以通过条件动作语句来编写控制规则，但其只能实现预定义规则的启发式控制，无法实现复杂的控制逻辑以及实时的寻优搜索。MPC的实施需要在模型模拟之外，逐时间步开发和设定控制寻优算法，并将得到的控制方案作用于模型的实时模拟过程。pyswmm提供了一个逐时间步的SWMM模型模拟接口，可以实时追踪模拟变量，并且让基于python编程的控制逻辑作用于系统模拟过程。
代码片段：基于pyswmm的实时模拟控制接口
from pyswmm import Simulation,Links,Nodes
from mpc import run_mpc
sim = Simulation('model.inp') # 初始化Simulation对
links = Links(sim) # 初始化Links对象
nodes = Nodes(sim) # 初始化Nodes对象
tanks = ['T1','T2','T3','T4'] # 4个调蓄池
orifices = ['V1','V2','V3','V4']# 受控对象：4个调蓄池的出流孔口
for st in sim:# 逐时间步遍历模拟
sim.step_advance(15*60) # 15 min为一个控制步长
ct = sim.current_time # 获取当前模拟时间
depths = [nodes[tank].depth fortankintanks]# 获取当前调蓄池水深
ctrls =run_mpc(sim,depths) # 运行MPC得到控制方案
for orf in orifices:
links[orf].target_setting = ctrls[orf]# 设置各孔口的设定值
sim.close() # 实时模拟结束
通过对pyswmm的Simulation模块对象的遍历，实现逐时间步的模拟，在每一次遍历中可以通过模型外部的python代码进行数据交换，获取实时模拟信息，设置控制方案。在pyswmm中，simulation、nodes、links、subcatchments等模块提供了对象的参数和变量的接口，可以在各模拟时间步中调用，查询当前的状态信息，具体的接口形式见pyswmm参考文档Reference — pyswmm 1.1.1 documentation。基于查询得到的状态信息，调用实时预测模块进行模拟评估和算法寻优，并将优化得到的当前步控制方案通过Links模块的target_setting接口传输给管网模拟过程，实现在各时间步的实时决策控制。
二、实时预测
在每一个时间步，MPC系统需要实时地对未来一段时间内的控制效果进行预测评估，结合优化算法实现控制方案的搜索和寻优。因此，实时预测过程需要在管网系统当前状态的基础上设置未来的控制方案。EPA-SWMM的热启动文件（hotstart file, hsf）可以记录模型模拟结束时的水文水力学状态，用于下一次模拟的起始边界条件。但由于实时控制在每一个时间步都需要预测评估，故存在需要实时导出热启动文件的技术问题。在加载了实时预测的起始条件之后，需要将优化算法生成的控制方案设置在SWMM模型中，对未来一段预测时长进行模拟，以评估其控制效果。
1.技术问题：SWMM模型热启动文件的导出
Sadler et al. (2019)提供了一个排水系统的MPC框架swmm_mpc，但框架建立在python 2的编程环境中，且2019年之后已不再更新维护，环境的兼容性和适用性较差。通过对该框架源码的解析，发现框架对OpenWaterAnalytics (OWA)开发的pyswmm进行了更改，建立了单独的pyswmm-feature_save_hotstart，增加了保存热启动文件的功能，并且通过对SWMM源码的重新封装将保存热启动文件的函数封装在swmm5.dll动态链接库中，在python环境中通过ctypes调用，实现将模拟中的实时水力学状态保存为热启动文件，作为EPA-SWMM模拟评估的起始边界条件。但不巧的是，该框架提供的swmm5.dll是32位的，在目前常用的64位python环境中无法调用，而重新配置32位环境过于复杂，因此只能另想办法。
代码片段：pyswmm-feature_save_hotstart中的swmm5.py部分片段
def save_hotstart(self, hs_file=None):
"""
saves a hotstart file
"""
self.SWMMlibobj.save_hotstart(ctypes.c_char_p(six.b(hs_file)))
既然无法直接调用SWMM程序导出hsf格式热启动文件，那么考虑用python语言将其C语言源码复现，实现在python环境中hsf文件的写出。hsf格式文件与SWMM的输出结果out格式文件本质上类似，均为按照固定顺序和格式将相应信息编码的二进制文件，因此只要能从pyswmm的模拟接口中获取到所需信息，完全可以在python环境中完成hsf文件的导出。问题的重点就在于解读hsf文件的编码顺序和格式。根据EPA-SWMM提供的源代码中的hotstart.c文件，可以看到hsf文件主要由以下三个部分组成：文件戳及对象数量基本信息、径流（runoff）和路由计算（routing）相关信息。以下内容基于EPA-SWMM 5.1.015版本的源代码总结整理，未经过准确求证，仅供参考。
（1）文件戳及对象数量
表1 hsf热启动文件头部基本信息
其中，流量单位flowUnits的对应关系为:
表2 流量单位代码对应关系
（2）径流Runoff
径流部分将当前时刻子汇水区的所有模拟状态信息保存至热启动文件中。对于每一个子汇水区（subcatchment），首先写入积水深度（depth）及总径流量（runoff），再写入渗透量（infiltration），若存在地下水（Ground Water）和积雪（Snowpack）的模拟，也会依次写入相应模块的参数。若存在水质模拟，首先写入径流中各个污染物（POLLUTANT）的浓度，再写入积水中各污染物浓度；最后对于每一个用地类型（LANDUSE），先写出各个污染物的积累值（Buildup），然后写出上一次清扫（Last Swept）的时间。当前状态的所有径流信息格式及顺序如表3所示。
表3 热启动文件中子汇水区径流变量编译顺序及格式
注意：
若子汇水区（subcatchment）只有一个分区（subarea），那么其参数则会分配给分区2；
渗透、地下水和积雪模型的基本参数详见相应的源码文件和SWMM水文学模块参考文档；
上表径流信息的所有变量数值均为双精度浮点数Double，占8个字节；
对于系统中每一个子汇水区，依次写入上表的所有信息。
（3）路由计算Routing
路由计算部分将当前时刻管道及节点的所有实时信息保存至热启动文件中。对于每一个节点node，首先写入水深（Depth）及横向入流量（Lateral Inflow）；若节点为水池（STORAGE）类型，则写出该节点的水力停留时间（HRT）；若存在水质模拟，则写出该节点的各污染物浓度。对于每一个管道link，首先写出当前的流量（Flow）、水深（Depth）和设定值（Setting）；若存在水质模拟，则写出该管道的各污染物浓度。
表4 热启动文件中管道节点计算路由变量编译顺序及格式
注意：
对于系统中每一个节点，依次写入上表中节点模块的所有信息；
对于系统中每一个管道，依次写入上表中管道模块的所有信息；
上表节点和管道信息的所有变量数值均为单精度浮点数Float，占4个字节；
在pyswmm模拟的某一时间步，创建hsf格式的二进制文件，利用bytes方法编译字符串，写入模拟的基本信息；通过调用subcatchments、nodes和links模块读取子汇水区、节点及管道的全部状态数据，采用python内置的struct.pack方法将以上读取的int、float及double格式的数据编译为二进制字符，并写出到hsf格式文件中，即可得到当前模拟时间步的热启动文件。
代码片段：基于pyswmm导出模拟过程当前状态的hsf热启动文件
from pyswmm import Nodes,Links,Subcatchments
from struct import pack
def save_hotstart(sim,config):
filestamp = 'SWMM5-HOTSTART4'
ct = sim.current_time
links = [link for link in Links(sim)]
nodes = [node for node in Nodes(sim)]
subs = [sub for sub in Subcatchments(sim)]
hsf_file = '%s.hsf'%ct.strftime('%Y-%m-%d-%H-%M')
hsf_file = path.join(path.dirname(config['inp_file']),
config['hsf_dir'],hsf_file)
with open(hsf_file,'wb') as f:
f.write(bytes(filestamp,encoding='utf-8'))
f.write(pack('i',len(subs)))
f.write(pack('i',0)) #nLandUses
f.write(pack('i',len(nodes)))
f.write(pack('i',len(links)))
f.write(pack('i',0)) #nPollutants
f.write(pack('i',['CFS','GPM','MGD','CMS','LPS','MLD'].index(sim.flow_units))) #FlowUnits
for sub in subs:
x = (0.0,0.0,0.0,sub.runoff) # ponded depths in 3 subareas, runoff
f.write(pack('dddd',*x))
x = (0.0,sub.infiltration_loss,0.0,0.0,0.0,0.0)
f.write(pack('dddddd',*x))
for node in nodes:
x = (node.depth,node.lateral_inflow)
f.write(pack('ff',*x))
if node.is_storage():
f.write(pack('f',0)) # HRT for storage
for link in links:
x = (link.flow,link.depth,link.current_setting)
f.write(pack('fff',*x))
return hsf_file
注：上述代码中的模型不包括水质模块、用地类型模块、地下水模块和积雪模块，仅包含基础的管网水力学状态信息和子汇水区径流信息。
2.预测评估文件的构建
在某一时刻（2009年11月3日15:30），构建并输出该时刻的hsf格式热启动文件之后，可用于EPA-SWMM的实时预测模型对未来控制方案的预测评估。具体而言，实时预测评估模型需要对原SWMM模型进行如下几点修改：
（1）在原来的inp格式SWMM输入文件的基础上插入如下的FILES模块，即可采用该时刻的hsf文件作为系统模拟热启动的初始状态值；
[FILES]
USE HOTSTART 2009-11-03-15-30.hsf
（2）将inp文件的模拟起始时间修改为当前时刻，根据预测时长设定模拟终止时间；
START_DATE 11/03/2009
START_TIME 15:30:00
REPORT_START_DATE 11/03/2009
REPORT_START_TIME 15:30:00
END_DATE 11/03/2009
END_TIME 16:30:00
（3）创建控制规则（CONTROLS），根据控制步的时间顺序（时间单位为小时）将未来一个小时的控制方案（15分钟一个控制步）依次写入不同优先级（PRIORITY）的规则中。
[CONTROLS]
RULE P1
IF SIMULATION TIME < 0.25
THEN ORIFICE V1 SETTING = 1.0
AND ORIFICE V2 SETTING = 1.0
AND ORIFICE V3 SETTING = 1.0
AND ORIFICE V4 SETTING = 0.4
PRIORITY 5
RULE P2
IF SIMULATION TIME < 0.5
THEN ORIFICE V1 SETTING = 1.0
AND ORIFICE V2 SETTING = 0.6
AND ORIFICE V3 SETTING = 1.0
AND ORIFICE V4 SETTING = 0.2
PRIORITY 4
RULE P3
IF SIMULATION TIME < 0.75
THEN ORIFICE V1 SETTING = 0.1
AND ORIFICE V2 SETTING = 1.0
AND ORIFICE V3 SETTING = 0.2
AND ORIFICE V4 SETTING = 0.2
PRIORITY 3
RULE P4
IF SIMULATION TIME < 1.0
THEN ORIFICE V1 SETTING = 0.1
AND ORIFICE V2 SETTING = 1.0
AND ORIFICE V3 SETTING = 0.5
AND ORIFICE V4 SETTING = 0.3
PRIORITY 2
在完成以上三步的修改之后，即可调用SWMM命令行工具进行模拟，预测该控制方案未来一段时间内的控制效果，并根据读取rpt文件得到的模拟结果不断优化控制方案。在python环境中可以采用swmm_api工具箱编写代码对inp文件进行操作和编辑，调用SWMM计算引擎模拟，并从rpt格式报告文件和out格式输出文件中读取结果（Picher, 2022），有关swmm_api工具的介绍之后有机会再整理。
结语
此次技术探索通过热启动文件的导出和预测评估文件的构建，可以解决排水系统MPC的实时预测部分存在的技术问题，构建出了城市排水系统的预测控制环境。在此基础上，将虚拟环境耦合进化算法、邻域搜索等求解组合优化问题的启发式优化算法，结合基于pyswmm的实时模拟模块，即可实现城市排水系统的实时优化控制，构建起一套完整的城市排水系统的MPC系统。
近期我在探索SWMM系统中MPC的技术实现，看到Sadler et al. (2019)文章之后受到较大的启发，无奈无法直接配置环境使用该框架，便有了这篇记录。此次探索是一个简易的技术尝试，若有错误内容欢迎交流讨论。
（城市智慧水务）