基于线性预测控制算法的明渠实时控制研究

孔令仲，王 浩，雷晓辉，权 锦，杨 迁

（1.扬州大学 水利科学与工程学院，江苏 扬州 225009；2.中国水利水电科学研究院 水资源研究所，北京 100038；3.河海大学 水利水电学院，江苏 南京 210000）

1 研究背景

修建大型长距离输水工程是解决水资源时空分布不均和缓解水资源供需矛盾的主要途径之一［1］。明渠输水作为长距离输水的重要方式，是指由一条长输水干渠和若干节制闸、分水闸等构成的复杂的输水方式，其存在的输水效率低、输水损失大是明渠输水系统中亟待解决的关键问题［2］。过去的明渠输水调控多依靠传统的人工经验，这种方法落后且不能达到精确调控。采用先进的明渠控制算法代替人工经验来进行明渠精确调控，是提高明渠输水效率的一个重要手段。长距离输水的水力学关键科学问题是水力控制［3］。

将水量适时、适量的供给沿线用水户是调水工程的首要运行目标之一［4］，而实现上述目标的关键在于将渠道运行控制点的水位控制在其目标水位。在每个渠池内，渠池下游端的水深最大导致边坡高于水位的值最小，且分水口一般靠近下游，因此每个渠池的水位控制点通常情况设为明渠下游端［5］。目前，关于明渠水位控制算法的研究主要集中在两个方面：明渠前馈控制算法和明渠反馈控制算法［6］。前馈控制算法主要包括蓄量补偿算法［7］和闸门步进法［8］，其基本思路为根据明渠仿真模型的模拟结果进行控制策略优化。这种方法的使用条件较为严格，一是只能用于仿真模型的边界已知情况，比如分水流量变化信息可知的情况；二是对仿真模型的精度要求较高。因此其实用性受到很大的限制，前馈算法一般不单独用于渠道控制。反馈控制的基本逻辑是根据实时的观测水位信息来生成调控方案，反馈控制算法既不需要预知分水扰动信息也不依靠于仿真模型，因此可以考虑渠系的随机扰动特性，满足用户取水灵活性和可靠性的要求。在反馈控制算法中，分布式比例积分（PI）类反馈控制算法［9-10］研究最为广泛，但PI 控制器属于单输入单输出（SISO）形式，在应用于多级串联输水工程时由于控制逻辑没有考虑渠池间的耦合作用，会出现水位变幅和流量变幅逐级放大的现象［11］。因此在多级串联输水工程中，为了在控制算法中考虑渠池耦合作用，国内外学者也开展了多输入多输出（MI⁃MO）控制模式的控制算法——集中控制算法的研究［12］。明渠集中控制算法主要包括线性二次型控制算法（LQR）和预测控制算法（MPC）。针对LQR 控制算法，Maleterre［13］基于有限线性圣维南方程构造了状态空间组来实现线性二次型最优控制；崔巍等［14］则进一步讨论了控制参数的取值问题，利用二次性能指标得到最优控制方案；王忠静等［15］在4 级串联渠池上设计了LQR 控制算法并与PI 控制算法的控制结果进行了对比，展示了LQR 控制算法在降低水位偏差和缩短稳定时间上的优势。针对MPC 控制算法的研究，Wahlin［16］讨论了基于简化渠道控制模型的MPC 算法设计并将算法应用于美国土木工程师协会推荐的模拟测试渠道Canal 1，结果表明基于简化模型设计的控制算法同样具有较好的鲁棒性；Overloop 等［17］首次将MPC 控制算法应用于实际的工程渠道中，结果显示尽管实际工程中存在着噪声和干扰问题，MPC 算法仍能控制水位保持在约束范围以内。综合而言，基于状态空间方程构造的集中控制算法能考虑渠系固有的时滞、耦合及随机扰动特性，是目前最适用于多级串联渠池中的控制算法。

基于此，本文以南水北调中线工程的最后6 节渠池为案例，研究在这种纵坡较缓、渠池的槽蓄量大且完全处于回水区的大型串联渠池中的水位控制问题，设计了MPC 控制算法并将控制算法应用于仿真渠道以检验算法的控制性能。

2 渠道仿真模型与控制模型

2.1 渠池仿真模型渠池仿真模型是进行控制算法研究的基础。本文中明渠仿真模型主要起到代替实际工程的作用，此外在控制模型的参数确定中也需要用到仿真模型。由于控制模型参数的确定并不严格要求仿真模型具有较高的精度，因此模型的精度问题本文不进行展开讨论。

渠池的主要水力建筑物为输水渠道、渠道两端的节制闸以及中间的分水口。其中，输水渠道采用Saint-Venant 方程组来进行仿真模拟：

式中：x 和t 为空间和时间坐标；A 为过流面积，m2；Q 为流量，m3/s；h 为水深，m；S0为渠道底坡；g 为重力加速度，m/s2； ql为单位长度渠道上的侧向出流量，m2/s，在这里不考虑渗漏的情况下设置为0；Sf为水力坡度，定义为：

式中：n 为糙率系数，s/m1/3；R 为水力半径，m，其中R = A/P，P 为湿周，m。节制闸采用过闸流量计算方程［18-19］模拟：

式中：Q 为过闸流量，m3/s；Cd为过闸流量系数；G 为闸门开度，m；b 为闸门底宽，m；hu为节制闸闸前水深，m； hd为闸后水深，m，在自由出流情况下设hd为0。其他建筑物可参考方神光等［20］的研究成果进行概化处理。

分水口是渠道向外界供水的通道，也正是分水口的分水流量变化造成渠池的水位变化，因此分水流量变化一般被认为是外界的扰动。分水口模拟主要是依据流量平衡方法：

式中：Qe为分水口前渠道断面的流量，m3/s；Qf为分水口后渠道断面的流量，m3/s；Qi为分水口分水流量，m3/s，在进行仿真模拟的时候需要给出。当分水口间距很小时，可以认为其瞬时水位相等，即有：

式中：He为分水口前渠道断面的水位，m； Hf为分水口后渠道断面的水位，m。

2.2 渠池线性模型圣维南方程组在进行渠池控制算法设计时，由于方程的微分特性，无法直接建立起控制点水位与调控点流量之间的关系，而离散后的圣维南方程组又存在着节点过多、难以实施优化算法的问题。鉴于此，Schuurmans［21］对圣维南方程在初始稳定态附近进行线性化假设和拉氏变换，将单个渠池概化为时间滞后特性为主的均匀流区和积分特性为主的非均匀流区。这样单个渠池的水力特性可由积分时滞模型来进行描述：

式中：e 为渠池下游端控制点的水位相对于初始设定水位的变化量，m；t 为时间，s； qin为渠池入口流量相对于初始稳定状态的变化量，m3/s； qout为渠池出口流量相对于初始稳定状态的变化量，m3/s；qd为渠池分水量相对于初始稳定状态的变化量，m3/s；As 为回水区面积，m2；τu为渠池的滞后时间，s。其中，参数As 和τu可通过前述仿真模型进行参数识别方法得到［22-23］。方程（6）可写成离散形式

式中：Ts为控制时间间隔，s；kd为滞后步长，kd=τuTs。将方程（7）的上一步展开并与方程（7）相减，可得：

式中：Δe( k )=e( k )-e( k-1) 为e( k )相对于上一步的增量；同理，Δqin( k )、Δqout( k )和Δqd( k )分别为qin( k )、qout( k )和qd( k )相对于上一步的增量。若将e( k )、Δe( k )以及Δqin( k-i )（i=1，2，…，kd）记为状态变量x( k )，Δqin( k )记为控制输入变量u( k )，Δqd( k )记为扰动变量d ( k )，则方程（7）和方程（8）可联合写为以下形式：

式中：A为系统矩阵；Bu为控制矩阵；Bd为扰动矩阵。而状态变量x( k )中包含的参数较多，还需要对状态变量x( k )和控制输出变量y( k )进行一定的转化：

式中C 为对角矩阵形式的输出矩阵。方程（9）和方程（10）构成了单渠池的状态空间方程。对于由串联多渠池组成的输水工程，可将各个渠池的状态空间方程联合写为具有方程（9）和方程（10）的形式的方程组，这样即可得到整个输水工程的状态空间方程。

3 预测控制算法

3.1 基本方程预测控制算法（MPC）使用真实系统的简化控制模型来进行未来预测，并通过最小化目标函数以获得最佳控制动作。本文采用线性状态空间方程（9）和（10）作为其控制模型，并利用控制模型对未来的输出进行预测，其基本原理是基于当前步长k 的实际状态变量x( k )、未来的控制输入变量u( k+i|k )以及未来的扰动变量d ( k+i|k )，对未来的输出( k+i|k )进行预测。其中，预测输出变量( k+i|k )中的i 的变化范围为从1 到p，p 为预测时域；而控制输入变量u( k+i|k )中的i 的变化范围为从0 到控制时域m-1，m 为控制时域，小于等于p。可将预测时域内各个时间步长的预测输出写成以下形式，则基于前述线性状态空间方程组，Y ( k+1|k )可采用以下公式计算：

式中：U ( k )=[u ( k )，u( k+1)，…，u( k+m-1) ]，D( k )=[ d ( k )，d ( k+1)，…，d ( k+p-1) ]，

预测控制的目标通常可表述为通过调整控制输入，使得预测时域内的控制输出变量与目标轨迹的偏差以及控制输入值最小。这里将目标轨迹定为水位偏差e( k )为零，则目标函数可表述为：

式中：Q 为输出偏差差权重；R 为控制输入变量权重。公式（12）为典型的二次规划问题，通过求解此二次规划问题可得到最优控制指令u∗( k )=[u∗( k )，u∗( k+1)，…，u∗( k+m-1)]T。但是只有第一个控制动作u∗( k )将被传输给控制结构用于执行，此时控制系统的状态变量滚动更新到下一个时段，并再次进行模型输出预测和目标函数优化求解等过程。

3.2 实测水位信息处理若渠池的实时水位波动明显，由于预测控制包含反馈控制环节，会根据实时的波动水位信息生成波动状态的流量动作，这样最终可能导致渠池发生共振而使得反馈环节失稳。因此，需要先对水位信号进行处理后再代入到预测控制算法中。这里，采用低通滤波器将波动水位信号进行平滑处理：

式中：eFJ( k )为k 时刻滤波后的水位差；Fc为滤波常数，，其中Tf为滤波时间常数，对于共振系统而言，其滤波时间常数可用以下公式计算：

式中：Rp为共振峰高；ωr为共振频率，，其中Tr为共振周期，这里用水波往返渠池上下游的时间估算：

式中：L 为渠池长度，v 为流速，c 为波速。通过式（14）和（15）即可计算每个渠池的滤波常数Fc，再利用每个渠池的滤波常数Fc即可对渠池的水位信号进行高频噪声过滤，采用滤波后的水位偏差eFJ( k )来代替实际的水位偏差e( k )。

4 计算实例

4.1 研究渠池研究案例渠池系统这里选择为南水北调中线工程的最后6 段渠池，整个6 段渠池被当成是一个独立的渠池系统。研究区域全长112 km，明渠断面为梯形断面，每个渠池的特征参数及初始流量值见表1所示。渠池系统上游假设为固定常水位边界，水深为7 m；下游也假设为一个固定常水位边界，水深为3 m。每个渠池末端有一个分水口向用水户供水。研究渠池系统的概化结构见图1所示。假设初始时刻水位都稳定在目标水位，输水系统为稳定状态。

图1 研究渠池系统概化图

表1 研究渠池基本参数及初始流量设置

通过施加分水扰动来让仿真模型中的渠池水位发生变化，并基于控制算法的计算结果来调控节制闸，以维持闸前的水位稳定在目标水位附近。已有研究表明分水扰动所在的渠池越靠近下游，渠池系统中需要参与调控的节制闸就越多，这就会使得渠池的水位控制越困难［24］，因此这里将分水扰动工况设置为分水口6 的分水流量增加5 m3/s。由公式（10）可知，MPC 控制算法可将未来的分水扰动信息代入到预测模型中去，并以此进行水位预测并进行调控指令优化，这种情况下的控制就包含了前馈控制环节。而多数情况下的小分水扰动是不可知的，只有在分水变化导致的渠池水位变化后，控制算法通过感知水位变化再进行渠池的调控，这种情况下控制算法就属于纯反馈控制。因此，对于MPC 控制中的分水变化工况，又可分为可预知分水流量变化和不可预知分水流量变化工况。为了对比MPC 控制算法的控制效果，设置了LQR 算法［14］调控下的结果作为对比。LQR 算法为纯反馈算法，其算法中不包含分水流量信息，因此不需要对分水是否可预知情况进行区分。考虑到研究渠池的底坡较缓且控制水深较大的特征会使得渠池中的水位波动坦化不明显，这里采用低通滤波器对算法中收集的水位信号进行处理，记这种加了低通滤波器后的MPC 和LQR 分别为MPC-F 和LQR-F 控制。同时，将不采用水位信号低通滤波器的控制算法作为对比。这样可设置5 种计算工况，见表2所示。

表2 计算工况

表3 渠池简化模型基本参数以及各渠池滤波参数

本研究中渠道仿真模型的模拟步长选择为1 min，控制步长选择为10 min。MPC 中需要确定的控制参数主要是权重矩阵Q、权重矩阵R、预测时域p 和控制时域m。对于权重Q 和R 的取值，经过试算，将R 设置为单位矩阵，将矩阵Q 中的对应于水位偏差e( k )部分的权重值设为5。预测时域p 应当尽可能包含所有渠池的滞后步长数，所以预测时域p 取值为50，而控制时域m 设为20。LQR 算法中只需要确定控制输出偏差权重Q 和控制输入权重R，采用和MPC 算法中相同的权重值。6 个渠池的简化模型基本参数以及滤波常数Fc分别见表3所示。

4.2 结果5 种计算工况下的结果分别见图2—图6所示。图2是MPC 算法应用于下游流量变化5 m3/s情况下的渠池1 和渠池6 中调控结果。图2（a）为分水口6 的分水增加后，渠池1 和渠池6 的进口节制闸的闸前水位偏差变化情况，图2（b）为渠池1 和渠池6 进口节制闸施加控制作用后的流量相对于初始流量的变化情况。在图2中，渠池1 和6 的水位以及进口节制闸的流量发生剧烈波动现象，导致渠池的水位变幅较大且水位变化一直处于共振状态，这种情况下的控制算法失效。图3则为对观测到的水位信息进行滤波处理后的MPC-F 控制算法的控制效果。从图3中可以看出相比于图2中的水位和流量变化过程，图3中的水位和流量变化过程较为平滑且最终水位稳定能稳定在目标水位、分水口6 上游的所有节制闸都增加了5 m3/s 以满足最终流量平衡。从图3（b）的节制闸流量变化过程可以看出在MPC 控制的作用下，尽管这里假设最下游分水变化情况不可知，但是分水扰动发生后扰动上游的所有节制闸能够快速进行流量响应，使得扰动上游节制闸流量动作幅度并没有发生明显的放大现象。这说明在MPC 控制算法下，节制闸能根据下游取水扰动引起的水位变化快速进行响应，使得扰动造成的水位变幅较小并使水位快速稳定在目标水位。

图2 MPC 算法应用于不可预知分水变化下的渠池1 和6 的控制结果

图3 MPC-F 算法应用于不可预知分水变化下的控制结果

图4和图5是LQR 算法以及LQR-F 算法应用于下游流量变化5 m3/s 情况下的调控结果。在图4中，尽管渠池的水位最终稳定在目标水位，但是在调控过程中也发生了轻微的波动现象。图5则为对观测到的水位信号进行滤波处理后的LQR-F 控制算法的控制效果，从图中可以看出最终水位稳定在初始水位，而分水口6 上游的节制闸也都增加了5 m3/s 以满足最终流量平衡，且调控过程水位变化相对平滑，这说明加入低通滤波器后的调控算法能一定程度上的消除明渠中的波动现象。

对比图5中的LQR-F 算法以及图3中的MPC-F 算法，可以看出在两种控制算法控制下的水位的变化过程和流量的变化过程比较相似。不同的是在图5（a）中渠池的最大水位偏差为-0.11 m，而图3（a）中渠池的最大水位偏差为-0.07 m，MPC-F 算法控制下的水位偏差相对较小；而流量的结果则相反，MPC-F 算法控制下的流量最大变幅相对更大，在图3（b）中节制闸的最大流量偏差为11.5 m3/s 而图5（b）中节制闸的最大流量偏差为7.5 m3/s。而且在LQR-F 控制结果的图5（a）中，当时间大于20 h 后渠池的水位是从低于目标水位恢复到目标水位；而在MPC-F 的控制结果图3（a）中，当时间大于20 h后渠池的水位有一个升高后回落的过程。这说明在LQR-F 控制中，渠池水位的上下回落现象不明显且节制闸的动作的回调次数较少。考虑到MPC 控制算法的目标为预测时间内最优，而LQR 算法为无限时间内的最优，可理解为MPC 的控制在初始阶段会更集中于水位差控制，故其水位差控制效果好于LQR 算法；但是从整个控制过程中的水位回落情况和流量回调情况可看出，LQR 控制算法显然更加合理。

图4 LQR 算法应用于不可预知分水下的控制结果

图5 LQR-F 算法应用于不可预知分水下的控制结果

图6 MPC-F 算法应用于可预知分水下的控制结果

图6为MPC-F 应用于分水可预知情况下的控制结果。从图6（b）的流量变化结果可以看出，在分水发生变化前，分水口上游渠池的进口节制闸的流量就开始发生变化，并且在10 h 的时候上游各节制闸的流量就增加了大约4 m3/s。相对应地，在图6（a）中，10 h 前节制闸的闸前水位也因为流量调整而发生变化。其中最明显的是渠池6 进口节制闸的闸前水位，在10 h 时已经升高0.04 m；而在10 h 分水口6 增加分水后，其水位立刻开始下降，在之后降为-0.02 m 后开始回升到目标水位。而其他渠池进口节制闸的闸前水位几乎没发生大的变化。这说明在预知分水情况下，MPC 控制算法能够预测未来分水变化可能造成的影响而采取提前动作，从而降低由于分水变化造成的水位变幅。

5 结论

本研究开展了基于线性控制模型的明渠MPC 控制算法研究。针对线性控制模型中无法处理的水位波动问题采用低通滤波对渠池的实时观测水位信息进行处理，并将上述方法以及LQR 算法分别应用于6 级串联仿真渠池。主要结论如下：

（1）MPC 控制算法和LQR 控制算法都属于多输入多输出的算法结构，在分水变化不可预知情况下，都充当纯反馈控制算法且控制结果类似。

（2）在完全处于回水区的渠池中，由于渠池下游端的水位在扰动下的波动较为明显，若直接采用实测水位值进行反馈控制，则会引起节制闸的流量变化过程也呈现波动状态，使得水位波动更为明显甚至发生共振导致控制失效；而在采用低通滤波器对实测水位值进行噪声过滤处理后再代入到算法中进行计算，能在一定程度上消除波动，实现渠池水位的安全稳定控制。

（3）MPC 控制算法的控制目标为有限时域内的目标最优，可将在有限时域内的优化目标与控制输入以及扰动的关系进行公式表达从而可通过算法实现考虑扰动下的目标最优。在可预知分水扰动信息的情况下算法能根据预测到的分水扰动进行前馈控制，在分水变化引起水位变化前就进行闸门调控，预先从渠首往下游补充对应的流量，从而起到降低分水扰动影响的效果。