目前, 我国的供水部门仅凭经验决策调度, 导致管网耗能巨大。同时, 由于很多地区的管网没有维护, 导致了严重的漏损, 无疑造成了严重的经济损失及能源的损耗。为了解决这些问题, 我国大规模地开展节能降耗技术的研究, 而水泵及阀门的优化调度是节能降耗必不可少的环节之一。许多研究人员提出多种不同方法来解决水泵与阀门联合调度的问题, 例如陈晨^[1]以水泵运行与管网漏失损耗费用之和作为目标函数, 以水泵转速及减压阀的开启度为自变量因素, 运用NSGA-Ⅱ模型求解, 成功地降低了管网的供水成本;Skworcow等^[2]提出了一种水泵-阀门两级调度的模型, 运用GAMS算法进行一、二级调度, 降低了水泵的能耗及管网的漏损。

　　 水泵及阀门一体化调度是一个大计算量的目标优化问题, 通常使用优化算法进行求解。梯度下降法 (GD) 是一种易于实现的算法, 其克服了传统优化算法的缺点, 但对于多自变量求解, 若自变量因素个数过多, 将导致结果陷入局部最优。Yang等^[3]在2009年提出杜鹃算法 (CS) , 采用莱维飞行并通过概率偏移, 可快速搜索到全局最优解, 同时也很少陷入局部最优。综上所述, 为达到均衡管网压力, 降低管网漏失量的目的, 本文通过两级目标函数对水泵、阀门实时一体化调度, 并选取CS进行求解。

　　 水泵及阀门实时一体化调度的具体流程如图1所示。首先将各节点用户需水量参数导入MAT-LAB-Open Water Analytics (OWA) 二次开发工具箱中^[4]。以水厂中各水泵的转速比为自变量因素, 以各时刻水泵运行供电成本最小为目标, 通过杜鹃算法 (CS) 进行优化求解, 从而降低水泵的运行费用。其次, 将调度后的水泵转速比导入管网模型中, 进行一次水力计算, 求出各时刻管段的流量。接着, 每次改变一个节流阀的开启度, 进行水力计算, 求出各时刻各个管段的压力差, 建立灵敏度矩阵, 并采用K-means算法, 对管网进行实时分区。最后, 根据同分区同开启度原则, 以各分区阀门开启度为自变量因素, 以各时刻管段漏损费用最小为目标, 通过CS进行优化求解。

　　 一级调度是以水泵转速比为自变量, 最小水泵供电成本为目标函数, 应用EPANET工具箱对管网进行水力计算, 并通过CS进行求解的模型。

　　 最小水泵能耗费用是指以各水泵的供电成本最小值为目标的函数, 如式 (1) 所示:

　　 式中W_t———时刻t各水泵的供电成本, 元;

　　 w_t———时刻t的单位电度电费, 元/ (kW·h) ;

　　 N_p———水泵的总数, 台;

　　 Q_jt、h_jt———分别为时刻t的水泵j的出流量, m³/h;出水压力, m;

　　 K———常数, 取0.010 19;

　　 η_jt———水泵j的效率, %;

　　 X_jt———水泵j的转速比, 水泵关闭为0, 全负荷运行为1。

　　 目标函数的计算除满足连续性方程及能量守恒方程外, 还应满足以下4个约束条件。

　　 水泵出水流量约束见式 (2) :

　　 式中q_{min, j}———水泵j的最小出流量, m³/h;

　　 q_{max, j}———水泵j的最大出流量, m³/h;

　　 q_j———水泵j的出流量, m³/h。

　　 水泵转速比约束见式 (3) :

　　 式中X_{min, j}———水泵j的最小转速比;

　　 X_{max, j}———水泵j的最大转速比;

　　 X_j———水泵j的转速比。

　　 水泵效率约束见式 (4) :

　　 式中η_{min, j}———水泵j的最小效率, %;

　　 η_{max, j}———水泵j的最大效率, %;

　　 η_j———水泵j的效率, %。

　　 各压力检测点压力约束见式 (5) :

　　 式中H_{min, i}———监测点i的最小压力, m;

　　 H_{max, i}———监测点i的最大压力, m;

　　 H_i———监测点i的压力, m。

　　 模型采用CS进行优化调度, 具体方法参考文献[3], 步骤如下: (1) 设置时刻为t, 通过MATLAB编写程序, 随机生成M组鸟群, 即M组水泵转速比X, 带入OWA工具箱进行水力计算, 并通过式 (1) 求得水泵的电耗成本; (2) 通过莱维飞行对自变量进行更新, 再次进行水力计算, 并与式 (2) ~式 (5) 限定条件进行比较, 若满足则通过式 (1) 求得水泵的电费, 不满足舍弃; (3) 多次迭代计算, 至电费趋于最优解为止; (4) 设置时刻t=t+1, 返回 (1) 重复运算, 至时刻t大于最大时刻为止。

　　 一级水泵调度后, 接下来对管网进行实时分区。K-means模型是通过各个质心的距离大小, 将m组数据分为c个类别的方法^[5]。若管网拓扑模型中有管段m个, 节流阀 (TCV) n个, 则K-means模型的计算过程分为6步:

　　 (1) 将时刻t的水泵调速比导入管网模型中, 阀门默认全开启的状况下, 通过OWA进行一次水力计算, 得出各管段的流量q_t⁰={q⁰_1t, q⁰_2t, …, q⁰_mt}, 其中q⁰_kt为在时刻t默认状态下管段k的流量, 单位为m³/h。

　　 (2) 调整阀门i的开启度β_i, 再次进行水力计算, 得出各管段的新流量q_tⁱ={qⁱ_{1 t}, qⁱ_{2 t}, …, qⁱ_mt}, 其中qⁱ_kt为在时刻t, 阀门i的开启度改变后的管段k流量, 单位为m³/h。

　　 (3) 求流量的变化值, 即流量灵敏度Δq_tⁱ={qⁱ_{1 t}q⁰_1t, qⁱ_{2 t}-q⁰_2t, …, qⁱ_mt-q⁰_mt}, 后还原阀门i的开启度。同理改变其他阀门的开启度, 对管段的灵敏度值进行计算, 最终形成流量灵敏度矩阵Q_t, 见式 (6) :

　　 式中Δqⁱ_kt———在时刻t, 阀门i开启度改变所导致的管段k流量变化, m³/h。

　　 (4) 在式 (6) 中选取前c列作为质心x_it, 其中i为聚类编号 (1≤i≤c) , 并通过式 (7) 计算出式 (6) 中各列到每个质心的距离。

　　 式中d (x_tⁱ, q_t^j) ———矩阵第j行到质心i的距离;

　　 xⁱ_kt———质心i中第k个元素。

　　 找出q_jt到距离最小的质心, 并将q_jt划分至其所在类中。

　　 (5) 分类完成后, 通过式 (8) 求出第4步中每一类中元素的平均值, 作为质心新位置, 并通过式 (7) 并重新聚类, 至类别中的元素不再改变为止。

　　 式中———第i类中元素的个数。

　　 (6) 令时刻t=t+1, 重复上述步骤对管网每一时刻的管段都进行一次分区。最终, m个阀门在每一时刻都被划分为c个区域。

　　 二级调度是动态实时分区后, 根据同分区同开启度的原则, 以各分区阀门开启度为自变量, 以最小漏失费用为目标函数的模型。

　　 动态分区后, 以同分区同开启度的原则, 对阀门进行调节, 从而最大程度地降低分区内由于压力过大所造成的漏失, 最小漏失费用为目标函数, 如式 (9) 所示。

　　 式中F_t———时刻t的管网漏失费用, 元;

　　 f_t———时刻t的水费, 元;

　　 N———管段总数;

　　 L_k———管段k的管长, m;

　　 c_k———管段k的漏失系数, c_k=2.5×10^-6;

　　 H_i、z_i、H_j、z_j———连接管段l的节点i、j绝对水头及地面标高, m;

　　 α———管段的漏失指数, α=1.18^[6]。

　　 4个约束条件应被考虑:

　　 (1) 各压力检测点条件约束, 同式 (5) 。

　　 (2) 各节点需水量与漏失水量之和与该时段水泵总出流量相等, 见式 (10) 。

　　 式中q_i———节点i的需水量, m³;

　　 q_ij———连接节点i的管段j流量, m³, 其流入为正, 流出为负;

　　 q_j^L———管段j的漏损量, m³;

　　 N_L———连接节点i的管段总数。

　　 (3) 阀门通过流量约束见式 (11) :

　　 式中T_{min, j}———分区j的阀门最小通过流量, m³;

　　 T_j———分区j的阀门通过流量, m³;

　　 T_{max, j}———分区j的阀门最大通过流量, m³。

　　 (4) 阀门开启度约束见式 (12) :

　　 式中β_{min, j}———分区j的阀门最小开启度;

　　 β_{max, j}———分区j的阀门最大开启度;

　　 β_j———分区j的阀门开启度。

　　 (1) 同水泵优化调度, 设置时刻为t, 生成M组阀门开启度β, 进行水力计算, 通过式 (9) 求得管网漏失费用。

　　 (2) 通过莱维飞行对自变量进行更新, 进行水力计算, 并与约束条件比较, 若满足则求漏失费用, 否则舍弃。

　　 (3) 多次迭代计算, 至漏失费用趋于最优解。

　　 (4) 设置时刻t=t+1, 返回 (1) 重复运算, 至时刻t大于最大时刻为止。

　　 通过M市的一实际供水管网作为验证对象, 管网中有15 587个管段, 1个定速泵 (800S-32) , 4个变速泵 (800S-47) , 2 142个TCV, 21个压力监测点, 平均供水量为5 468m³/h, 管网拓扑结构如图2所示。

　　 以式 (2) ~式 (5) 、式 (10) ~式 (12) 作为约束条件, 水泵的出流量q_j取值范围在0至出厂流量之间, 且各水泵出流量之和小于该时刻总供水量;水泵转速比X_j在0.3~1;水泵效率η_j在0~100%;监测点压力H_i在30~40m;阀门通过流量T_j在0~Q_t;阀门开启度β_j在5%~90%。在动态分区过程中, 设置开启度β=60%, 对管网进行分区。

　　 以式 (9) 的漏损最小值为目标, 分区数目c从2~10进行试验, 结果如图3所示。由图3可知, 将管网分为4个区域时, 漏损成本最小, 因此设定分区数目c=4。

　　 设M市的工业高、中、低峰供电成本w_t为1.29元/ (kW·h) 、0.84元/ (kW·h) 、0.42元/ (kW·h) , 供水成本3.19元/m³。选取管径DN100~400的TCV进行开启度优化。在EPANET软件中, 通常以局部阻力来模拟开启度的大小, 阀门阻力系数与开启度关系公式见文献[7]。根据上述已知条件, M市管网的水泵及节流阀实时一体化调度过程如下: (1) 导入2016年1月1日12h各节点需水量数据至管网模型; (2) 运用CS算法对各时段的供电成本最小值进行求解, 求得最优水泵转速比, 如表1所示, 调度后可将定速泵设置为备用泵; (3) 计算各时段的灵敏度矩阵Q_t, 运用K-means算法对各时刻阀门分为4个区域, 并运用CS算法对各时段的管网漏失费用最小值进行求解, 得最优阀门开启度, 如表2所示。

　　 计算每一时刻水泵能耗及各分区的漏损费用之和, 并与单独运用式 (1) 或式 (9) 对单水泵、单阀门优化调度以及优化前的费用结果相比较, 如图4所示。

　　 通过水泵与阀门两级调度后, 总费用共16 175.49元, 比单水泵调度费用16 216.18元, 节约了40.7元;比单阀门调度费用16 504.74元, 节约了329.25元;比传统方法费用16 543.91元, 节约了368.42元。

　　 运用水泵-阀门两级优化的方法, 解决了动态分区、水泵-阀门实时一体化的调度的问题, 通过M市供水管网系统的实例验证可知, 案例模型优化后, 水泵能耗及管网漏失费用共16 175.49元, 比单水泵调度花费节约0.25%, 比单阀门调度花费节约1.99%, 比传统方法花费减少2.23%。因此, 该方法效果良好, 可应用于实际管网。