基于数值仿真的单目标多变量渠道控制参数寻优算法研究

张雨萌，管光华

（武汉大学 水资源与水电工程科学国家重点实验室，武汉 430072）

0 引 言

【研究进展】自20 世纪30年代水力自动闸门的应用开创了渠道系统自动化控制的先河以来，国内外学者不断研究创新，使得渠系自动化控制理论日趋成熟，相关的技术得到了一定范围的应用。渠系自动化控制的核心是控制算法，即描述渠道输入（水位或流量误差）与输出（闸门动作）的逻辑关系的方法。渠系控制算法先后出现了基于经典控制理论的渠道控制逻辑（如小人算法、科文算法以及后来的PID 类算法等[1]）；基于现代控制理论的渠道线性二次最优控制算法、最优下游水位控制器等；以及基于智能控制算法的鲁棒算法等[2]。

近年来，渠系控制理论的各个方面又有了一些新的进展。姚雄等[3]对大型渠道控制蓄量的实现方法进行了研究，建立了包含倒虹吸、渡槽、取水口等建筑物在内的大型输水渠道自动运行控制仿真模型，并对京石段进行了仿真验证。Clemmens[4]提出了一种解决用水管理部门与用户的供需不匹配问题的基于相邻渠池的水位差异的控制器，其本质是将渠道系统看作一个大型水库。Belaud 等[5]针对渠系响应时间展开研究，运用IDZ模型寻找基于用水需求改变的开环系统的恰当时间，并将其在一个渠系系统的不同的下游边界条件下与SV 方程做比较，提出了灌溉渠系的响应时间概念。崔巍等[6]提出了一种具有较好实用性的闸前常水位运行控制算法，综合前馈反馈算法提出动态调节渠道蓄量实现闸前常水位运行的控制算法，在南水北调的仿真测试中取得了良好效果。管光华等[7]对蓄量阶跃补偿、动力波原理、水量平衡模型3种时滞参数算法在ASCE的测试渠系和漳河灌区上进行仿真模拟，比较了它们的控制效果，从而推荐蓄量阶跃补偿算法为主要的前馈控制算法，并为简化前馈规则提出了时滞参数显式算法。

【拟解决的关键问题】上述控制理论及算法尚处于试验和研究阶段，在实际工程中，经典控制应用较为广泛，特别是原理较为简单、鲁棒性较强的PID算法[2,8-12]。该算法在确定控制器参数时，主要有2种方法：一是传统的经典设计方法，主要包括构建ID模型、IDZ模型等，并将其结合最优控制理论推算得到PID 参数；二是进行程序仿真试算，用圣维南方程结合寻优程序确定较优PID 参数。但由于方法一过于理想，没有考虑真实渠道的局限性、非线性以及系统的通信延迟，故基于经典控制理论的设计方法、模型所得结果往往不够理想，甚至出现控制失稳的情况。【切入点】目前，对于渠道控制系统PI 控制参数的寻优大多采用方法二中的网格法，本文采用仿真模型来率定和优化控制器参数，旨在将一些在其他领域运用较为成熟的算法（如梯度算法、PSO算法）改进运用于渠道控制系统中。【研究意义】以期能够在渠系运行过程中较快地取得较好的PI 控制参数，使渠道达到控制目标，安全可靠及时地保障渠道系统高效运行。

1 PID 控制逻辑及其渠道控制性能指标

PID算法又叫作比例积分微分算法。此处的应用原理为：传感器实时监控下游水位，将该水位与预定的常水位做比较，将二者之间的差记为误差值e(t)。然后将该误差值进行比例积分微分的线性组合计算，找到一个输出的闸门开度u(t)作为控制器输出变量[2]，计算式为：

式中：Kp为比例项因子；Ti为积分项因子；Td为微分项因子；e(t)为水位误差值（m）；u(t)为闸门开度（m）。比例项是为了减少水位的误差值；积分项是为了减少水位的稳态误差；微分项是为了加快闸门的动作量，本文只考虑比例项和积分项[13]。

2 单渠池控制参数寻优算法

选取美国加利福尼亚州ASCE的小型测试渠道的单个渠池作为本文初步的寻优仿真对象，渠道参数见表1。流量变化工况设置为设计流量的15%～30%的阶跃流量变化。

表1 单个渠池工况仿真渠道参数Table1 simulation parameters of single channel pool

2.1 网格法寻优及优化目标

网格法是PI 参数寻优最为经典的算法，其优点在于只要计算时间足够长，取足够密的网格进行寻优，就一定可以得到性能指标较好的一组控制参数。本文先以经典网格法进行寻优，将其作为一个比较对象与其他PI 参数寻优算法在运算效率、控制效果上做出比较，来评判其他算法是否有效或更优。

参考ASCE 提出的基于水位误差变化幅度、渠系达到稳定所需时间、闸门动作量等关键性控制目标的性能指标，包括水位误差平方积分（ISE）、最大绝对误差（MAE）、稳态误差（STE）、绝对值误差积分（IAE）、绝对流量变化积分（IAQ）、绝对闸门开度积分（IAW等）[14]。其中，ISE是用于衡量水位波动平稳性能的指标；IAQ是描述系统在过渡过程中对流量的控制能力的指标；IAW是衡量闸门的启闭幅度和频率的指标。ISE越小，水位波动越小。初步探究时，本文认为若能控制水位波动在一个合理的范围内，就能够使系统达到一个较好的被控状态。故采用较为传统的ISE（水位误差平方积分）作为寻优指标，计算式为：

式中：Δt为时间间隔；T为总的仿真时间；yt为当前时刻水位；ytarget为目标水位；ISE为水位误差平方积分（此处为去量纲化，且多渠池工况下取平均值）。

综合参考文献[15]，本文中Kp、Ti的取值范围为（0,100）。采取参数计算步长为2，对Kp、Ti进行（0,100）内的单渠池寻优，需计算50×50 共2 500个性能指标，共耗时约20 h（用于运算的电脑配置如下：Intel(R)Xeon(R)Gold5118 处理器，32GB 内存）。若再从中选取区间进行加密，需要约160 h 才能得到理想的控制参数，而本文2.4 节的计算结果显示，运用本文提出的2种单渠池寻优方法，最多只需要8 h 就可以找到同样控制效果的控制参数。可见网格法所付出的时间代价十分大，亟须改进。

2.2 复合梯度算法寻优

2.2.1 梯度法

梯度法是数学中较为经典的算法，广泛应用于军事、工业等方面的寻优。它利用负梯度方向是下降最快的方向的原理进行寻优，其迭代公式为：

式中：a0为步长；p(k)为目标值y对于变量的负导数；k为寻优次数[16]。

2.2.2 复合梯度法及其寻优效果

寻优分为2 部分，即随机二分寻找梯度法初值以及梯度法寻优，称为复合梯度法。由于不合适的控制器参数会使得系统失稳或性能指标显著劣化，且梯度法寻优需要一个较好的初值。故认为梯度法控制参数的用于寻优的初始值对应的性能指标应优于没有反馈控制（只有前馈控制）时的性能指标。

首先为了适应梯度法的初值敏感性，先使得PI参数（Kp、Ti）在一定范围内（0～100）随机生成n×n个值，并计算对应的性能指标值。若得到的最小的性能指标仍大于没有反馈控制时系统对应的性能指标，则先数量级二分再区间二分，不断缩小逼近最优范围直到找到合适的梯度法初值。然后进行梯度法运算，以负梯度方向为最速下降方向为基础，使得PI 参数对应的性能指标越来越好，直至达到预期性能指标（目前经仿真实验验证，预期性能指标取无反馈控制对应的性能指标的1/50～1/100，渠系控制效果良好），不断逼近最优的PI 参数。

程序设计为当型循环，采取MPIEC（第i次的最小性能指标和i-1次最小性能指标的差）和MPIEL（第i-1次的最小性能指标和i-2次最小性能指标的差）作为循环条件，当二者都达到较小的程度（取预期性能指标精度的1/10），认为寻优进入了一个较小值的平缓区域，则停止循环，并将性能指标达到预期值作为结束循环的另一个条件。当遇到控制出错的情况，例如闸门上游水位低于下游水位、上游闸门全关等情况。自出现控制出错的时刻起，将性能指标扩充到整个寻优时间段并将对应性能指标扩大1 000倍（经试验，每次遇到控制出错的情况都会加长寻优时间，且此组性能指标必然不能用）。

使用梯度法原理进行PI 参数的迭代，计算式为：

式中：i为迭代次数；a0、b0为步长；IK、IT为ISE对Kp、Ti的偏导数的负值。运用式（4）和式（5）时，将其中a0、b0采用控制数量级的方法，使得加数的第2 项的数量级近似比加数第1 项的数量级少1个。经过多次验算，该方法较其他方法更加有效稳定。

图1 过闸流量过程线 Fig.1 Flow process line

图2 闸门开度过程线Fig.2 Gate opening process line

图3 下游水位误差过程线Fig.3 Error line of downstream water level

表2 梯度法时间及循环次数统计Table2 Statistics of gradient method time and cycle times

2.3 复合PSO算法寻优

2.3.1 PSO算法 PSO算法，即粒子群算法，是基于鸟群觅食原理设计的[17]算法。先一次性撒下许多粒子，随机生成N行D列的位置和速度数组。通过计算每个粒子的评价指标，可以筛选出此时的全局最优点。再进行M次迭代，对每个粒子的位置和速度进行一定规律的运算。该运算公式有许多不同的公式，综合考虑寻优时间与寻优效果[18]。本文选取了较为基础经典的一种计算方法，计算式为：

式中：v为粒子的速度；w为惯性权重；x为粒子所在的位置；i为粒子序号；c1、c2为比例因子；y为存储x值的临时存储变量；pg为当前最优控制参数的存储变量。

关键参数选取如下：N为初次撒下的点，本文认为初次撒下的点的效率较低，迭代计算的效率较高，故初选N=20；D为计算的维度，本例的研究变量仅有Kp、Ti共2个值，故D=2；w为惯性权重，认为粒子下一时刻的运动速度受上一时刻的运动速度所影响。此处取常规值0.6；M为迭代次数，是根据PSO算法的核心原理进行迭代计算，故取M=80；c1为粒子受自身最优值影响的比例因子；c2为粒子受群体最优值影响的比例因子。二者之和一般为4，目前取二者均为2[19]。

按此公式对每个粒子进行速度和位置的改变，并计算对应的性能指标来进行评判。不断循环，直至找到满足条件的性能指标，并筛选出对应的PI 参数，完成粒子群算法寻优[20]。

2.3.2 复合PSO算法及寻优效果

为提高寻优效率，参考2.2 节先使用随机二分法生成N个寻优粒子（粒子含有位置：即Kp、Ti值；速度：下一步寻优的方向和步长），再利用PSO算法原理，使得N个粒子以不同的权重沿着上一时刻前进的方向、自身最优的方向、种群最优的方法进行前进循环，直到找到符合预期性能指标的PI 参数。

随机二分法设计同梯度法，粒子群算法设计为：定义种群数（N=20）、粒子维数（D=2）、自身最优值学习因子（c1=2）、群体最优值学习因子（c2=2）、迭代次数（M=50）。先随机生成N个二位数组x(i，n)（i=1～n），其物理意义为：随机生成20 组Kp、Ti，并用x(i,1)储存Kp的值；x(i,2)储存Ti的值。再随机生成N个二维数组速度v，以后期改变对应的第i组Kp、Ti值。

若该区域仅从空调风机盘管的角度考虑，可以采用以上公共走道和电梯厅采用的联网型风机盘管的做法。本案例采用另外一种做法，借鉴客房客控RCU控制的方式，既可以实现远程对风机盘管的控制，实现效果等同于联网型风机盘管，又可以通过设置红外探测器联动公共卫生间内的照明开关，以实现节能效果。裙楼公共卫生间灯光和空调控制布置如图6所示。

计算每组PI 参数对应的性能指标。并将当前20组Kp、Ti值参数存入到y矩阵里，并从中选择对应性能指标最小的那组性能参数对应的PI 参数存入全局最佳性能指标对应PI 参数值pg里。对应的最佳性能指标值存入到c变量里。

再执行当型循环，循环条件为：迭代进行50次以内、性能指标达到预期值或者寻优进入平滑区域。多次根据粒子群算法核心公式更新鸟群的速度和位置。并不断以更小的性能指标赋给变量c，不断以对应的PI 参数赋给变量pg。若跳出当前循环后发现性能指标没有达到预期值，则重新寻优。以此重复，直到鸟群的最佳性能指标小于预期值为止。最终的位置为pg，pg的第1 列为Kp，第2 列为Ti。c为最后的最佳性能指标[21-22]。

同样，罗列出15 组计算结果如表3所示，表中TP表示粒子群法所用时间，NP表示粒子群法循环次数。水位误差、闸门开度、过闸流量过程线显示与复合梯度法得到的图1—图3控制效果相当：控制很快达到了稳定，水位误差较小、控制效果也较为良好。

表3 PSO算法时间及循环次数统计Table3 Statistics of PSO algorithm time and cycle times

2.4 复合梯度法与复合PSO算法比较

本文以梯度算法与PSO算法寻优的各15 组算例进行相同寻优效果下寻优时间的对比，结果见表4。在控制效果相当的情况下，就其用时来说，2种算法较均匀网格法的效率均得到很大程度的提高，最多用时仅约为网格法的5%。但总体来说，复合梯度法所用平均时间仅仅为复合粒子群算法的35.76%，复合梯度法较复合粒子群算法更高效一些。二者对于初值的选取都较敏感，都需要先使用随机二分法初步寻找主体算法初始值。二者的收敛性都较好，寻优一段时间均能找到良好的控制参数。

表4 2种算法寻优时间比较Table4 Comparison of optimization time between the two algorithms

3 多渠池工况下高效寻优

本文选取美国加利福尼亚州ASCE的小型测试渠道中连续3个渠池作为本文进一步的寻优仿真对象，渠道糙率为0.02，底坡为0.000 1，底宽为7 m，边坡为1.5，渠道长度见表5。工况设置为30%～40%的阶跃流量变化。将复合梯度法运用于多渠池工况，进行多变量寻优。并考虑到其复杂以及耗时性，先采用固定权重的多渠池多变量寻优，再采取随机多渠池多变量寻优，由浅及深逐步探索。

3.1 固定权重的多渠池多变量寻优

将单渠池问题扩充到多渠池问题后，应该是6个变量的寻优，寻优难度会大大增加，效率也会迅速降低，故先采用简化算法进行探索。即运算中第1个渠池的控制参数Kp、Ti仍然沿着负梯度方向进行寻优，而其他2个渠池的控制参数利用渠池之间的关系进行直接赋值（寻优指标为多渠池对应的性能指标）。并在其过程中，以单渠池梯度法为基础，充分利用单渠池寻优所发现的结果和所利用的技巧，如解决梯度法初值敏感性的初值的随机二分，选取步长前进所采用的控制数量级法，将肯定会出错的性能指标根据时间占比扩充到整个时间段内等。

根据不同渠池的水面面积之间的关系决定不同渠池的控制参数的关系做法较为合理。根据相关文献研究可知，渠池长度越长（水面面积越大），敏感性越小，对应的控制参数越大[23]。故本研究采用各个渠池的控制参数与其对应水面面积成正比的方式来选取不同渠池的Kp、Ti的值（Ti为反比关系）。

表5 多渠池各渠道长度 Table5 Length of channel

为了解决梯度法的初值敏感性问题，多渠池仍然采用先随机二分，来取得大于无反馈控制时的性能指标所对应的寻优参数值，再将取得的初始值代入多渠池梯度法的核心算法中的方法，即随机给定第1个渠池的Kp、Ti，再根据水面面积之比给定其他2个渠池的寻优参数。将这6个有一定关系的寻优参数先进行随机二分法寻优，找到对应性能指标大于无反馈控制的一组值，作为梯度法的合理初始值代入多渠池的梯度算法，以期得到目标性能指标（此处的目标性能指标的数量级选取无反馈控制时性能指标的约1/100。若真正用于工程中，选取真实工况的预期性能指标精度即可）。

取得3个渠池的初始寻优参数后，沿着负梯度方向寻找第1个渠池的性能指标，而其他渠池的性能指标则采取每一步都根据水面面积之比给定。同单渠池算法，比较本次寻优参数所求得的性能指标和之前的性能指标，将更小值存在变量minJUG里。并求得MPISC、MPIEL不断循环，直到MPISC和MPIEL都较小，说明寻优进入了一个较为平坦的区域停止循环。或者minJUG达到预期精度，则停止循环。最终输出对应的性能指标和对应的寻优参数Kp1、Kp2、Kp3、Ti1、Ti2、Ti3。进入梯度法循环80次，若80次后仍未找到预期的性能指标，则进行重新寻优，重复以上步骤，不断循环直到找到预期值。（根据实验结果，此种情况循环2次左右即可找到预期值）完成多渠池梯度法运算。

得到的结果任取10 组取最大值、最小值、平均值列于表6。由表6可知，固定权重的多渠池多变量寻优法的平均寻优时间为4 h，较网格法大大提高了寻优效率。各项性能指标均较小、稳定时间较短、闸门运动幅度较小。从图4—图6所示过程线可得，水位误差、流量过程变化均较为合理，控制效果良好。

3.2 随机多渠池多变量寻优

为取得更好的控制效果，本文进行了6个变量的复杂多渠池寻优。仍然采用将寻优分成2个部分。第一部分为随机二分粗选梯度法初值，第2 部分为梯度法寻优。与3.1 节的区别在于，第1 部分采取随机生成6个毫无关系的变量，以及第2 部分梯度法寻优也是多渠池性能指标分别对6个变量求差分，6个变量都沿着自身的负梯度方向进行寻优。剩余算法处理同上一节大致相同，只是运算量由于求了3次差分而大大增加。

得到的结果任取10 组，其最大值、最小值、平均值列于表7。可见随机多渠池多变量寻优性能指标较小、控制效果良好。只是寻优时间远大于固定权重的多渠池多变量寻优。

图4 过闸流量过程线 Fig.4 Flow process line

图5 闸门开度过程线Fig.5 Process line of gate opening

图6 水位误差过程线Fig.6 Downstream water level error line

表6 固定权重梯度法寻优结果Table6 Optimization results of fixed weight gradient method

表7 多渠池多变量寻优结果Table7 Multivariate optimization results of multi-channel pools

表8 多渠池工况下2种算法结果比较Table8 Comparison of two algorithms under multi-channel pool conditions

3.3 2种算法比较

2种算法的各项性能指标结合表6—表8 进行比较，随机多渠池多变量寻优的控制效果较固定权重的多渠池多变量寻优有较大的改善，随机法较固定权重法的IAQ小90.48%、IAW小77.57%，Tstable（系统稳定时间）小17.05%，ISE小50.42%。根据2.1 节各性能指标所示物理意义可知，随机法所得的控制器对流量的控制能力更强，闸门的起闭幅度和频率更小，稳定时间更短，水位波动更小。但随机法的所付出的时间代价也较大，较固定权重法多63.19%。用户应综合考虑利弊，根据自身需要来选取更为适用的方法。

4 讨 论

目前，在水利工程领域研究PI 参数试错法寻优的学者较少，其中，Stringam 等[12]提到的迭代法寻优就是实际应用中PI 参数寻优的经典方法。该法高度依赖经验，需要渠道运行人员能够依靠经验给定PI参数较为合理的比例项因子初值。首先只进行比例控制，选取较为合理的比例项值；再固定比例项因子的值，以一定的规律改变积分项因子，直至渠系达到较好的控制效果。该法较为耗时，仅进行单渠池寻优也需要近8 h。对于缺乏良好管理经验的管理人员的渠系，给定不恰当的PI 参数容易使渠系控制失调。且在较为紧急的渠系水量需求改变的工况下，也不能及时满足用户需求。而本文提出的几种算法具有较好的普适性，对渠道运行管理人员经验的要求不高，且控制的稳定性和安全性更有保障。

但本文的研究仍存在一些有待改进的方面：多渠池工况由于其存在不同渠池之间的耦合性，可能会造成误差累计或放大。且随着寻优渠道数目的增加，寻优变量也成倍增加，造成寻优时间大幅增加。故本文先采用固定权重的多渠池多变量寻优算法。之后，为寻找更好的控制指标，采用了随机多渠池多变量寻优法，将所有变量均沿梯度法进行寻优，虽然得到的控制效果更好，但所付出的时间也多了约60%。在接下来的研究中，笔者将尝试寻找控制参数与渠道特性及对应工况的更深层次的关系，以求能够在尽量减少寻优时间的情况下，寻找控制效果更好的控制参数。

本文仅针对ISE进行寻优，只考虑了控制水位误差较小这一优化目标。采用该性能指标的缺点是会造成闸门的过度调节，实际输配水工程的运行调度需综合考虑流量误差、闸门运动、系统稳定时间等方面的要求。故可采用更加综合的寻优指标，以获得更为均衡的控制效果。

实际应用中，寻找仿真性能最优的控制器往往不是解决控制问题的良策，因糙率、闸门流量关系等水力要素的不确定性及时变性，“最优”的控制器在应用中并不一定有最佳的性能，而调度者更为关注的是控制性能在较大的范围内能保持较好的水准，及控制器应具有较好的鲁棒性。故可考虑进一步将控制性能的梯度引入到寻优目标中。

5 结 论

1）本文提出的复合梯度法和复合PSO算法都能较好地实现渠道控制系统PI 控制参数高效寻优。在控制效果相当的情况下，复合梯度法的寻优时间缩短为原来网格法的0.5%，PSO算法的寻优时间缩短为网格法的1.4%。

2）2种复合寻优算法都可以较好地解决寻优初值不稳定及寻优效率低下等问题。

3）多渠池工况下，随机多渠池多变量寻优算法与固定权重的多渠池多变量寻优算法所取得的PI 控制参数都能较好地对渠道实现控制。随机多渠池多变量寻优算法的控制效果明显优于固定权重的多渠池多变量寻优算法。前者水位误差平方积分（ISE）是后者的50%，稳定时间（Tstable）是后者的17%。但固定权重的多渠池多变量寻优算法的时间成本远低于随机多渠池多变量寻优，是随机多渠池多变量寻优时间的37%。在进行控制器设计时，应视需求选择合适的算法。