序列二次规划法在多水源管网优化调度中的应用研究(1.浙江大学 水工结构与水环境研究所，浙江 杭州 310027；2.杭州市自来水总公司，浙江 杭州 310016)
1 问题的提出
城市供水管网是城市的生命线之一。这一复杂的网络系统，主要通过几个供水泵站为城市血液提供能量，送至城市的各个角落。通过对供水泵站的优化调度，可以降低企业的制水成本、使管网的供水压力分布更合理，根据初步估计，对于一个日供水量为10万吨的自来水公司，如果供水扬程降低1m，每年可以节电15万kW·h；由于管道系统的渗漏与水头有关，降低供水水头也可以在一定程度上减少管网的渗漏；供水水头的降低还可以减少爆管的风险，这对于管网的管理有更深刻的意义。因此管网合理调度研究一直是供水企业一个重要课题，同时也是一个难题。
供水管网运行的合理调度可以用一个最优化问题来描述，管网的水力方程组是一组非线性方程，各水源水泵的开启状态作为离散变量，因此这是一个混合变量的非线性最优化问题。由于离散变量与连续变量的同时存在，求解极为不方便，最常用的方法是将该问题分作两级进行优化：一级优化是针对管网而言，目的在于求各水源的最佳供水量或最佳供水扬程；二级优化是在一级最优化的基础上，根据水源的具体情况，确定满意的水泵开启方案和水泵的调速比。采用以上方法可以在一定程度上降低求解的困难，但一级优化也是一个非线性的优化问题，求解起来相当麻烦，目前国内外最常用的方法是广义简约梯度法。广义简约梯度法虽属较优秀的约束非线性规划算法，根据作者在以往其他优化应用方面的研究，其重分析次数相当多。在本优化问题中水力计算是计算量的主体部分，由数值试验的经验知，在目前中等配置的微机上完成一个2000个左右节点的供水管网，一次水力计算需要10s左右，如果采用广义简约梯度法，需要反复迭代计算，花费的时间是相当可观的。由此可见采用广义简约梯度法实现管网的在线优化调度存在较大难度。
本文针对一级优化问题，采用序列二次规划法进行求解。M.J.D.Powell所给出的序列二次规划法实质上是运用KuhnTucker最优化条件所形成的非线性方程进行迭代计算，而这一迭代过程恰好可以用求解一相应的二次规划问题替代，故原问题的求解过程转化为求解一个二次规划的序列。其中二次规划问题的二次目标函数是原问题Lagrange函数的二次展开式，包含了目标与约束函数的二次信息。通常其二阶导数矩阵由变尺度的思想通过先前迭代点的梯度信息逐步生成。序列二次规划法综合利用了KＴ条件、变尺度、线性及二次近似等有效手段，在理论上是一种比广义简约梯度法优秀的算法[1]，但是它的迭代序列通常从不可行域逐步逼近可行域，需要在极限情况下才能完全达到约束要求，这显然不利于尽快获得可行的较优解，故约束条件的妥善处理非常重要，本文将结构优化中齿行法的思想和水力学的基本概念相结合，提出了一种新的算法，可以方便地将迭代中的非可行点拉回到约束界面上，获得了较高的计算效率，有助于实现管网的在线优化调度。
2 供水优化调度一级优化的数学模型
管网的运行调度一般以经济性作为目标函数，与水源的供水量、供水水头有关，据此可以建立供水管网的目标函数：
式中：FG为各水源的制水成本和供水的动力费用；Qs、Hs为各水源的供水量和供水水头。
供水调度的主要约束条件有：管网的水力关系，各水源的水量和水压的约束，管网中各节点的最小服务水头。这些约束条件分别表示如下：
管网水力关系
各水源的供水水头约束
各水源的供水量约束
Qsmin(Hs)≤Qs≤Qsmax(Hs) | (4) |
管网各节点服务水头约束
| (5) |
其中：Hsmax、Hsmin分为水源的最大、最小供水水头；Qsmin(Hs)、Qsmax(Hs)分为水源的最大、最小供水能力，通常水源的供水量的能力与供水水头有关。HN为管网中各节点的服务水头；HNmax、HNmin分为管网中各节点的最大、最小服务水头；QN为管网中各节点的节点流量。
3 模型的求解
模型求解主要有2个难点：(1)约束条件太多，一个中等复杂的城市管网可能会有上千个约束；(2)目标函数中各变量隐式相关——水源的供水水头Hs和供水水量Qs隐式相关。如果能对以上两个方面进行适当的处理，可以大大的降低难度，提高求解效率。针对以上两点，本文从管网的水力条件出发，提出了一套求解方法：
在一定负荷ＱN下，将管网的水力计算公式(2)在H0处作一阶泰勒展开有：
| (6) |
, 
称为敏度矩阵。如果用哈真-威廉公式表示管道的能量损失，用矩阵A、B可以分别表示为 , ，管网的水力学公式可以用式(7)表 达。A和B仅与管网中管道的水力坡度有关。当任一水源的供水水头发生变化，由于管网自身的调节作用，每根管道的水力坡度的变化幅度要比节点水头变化小得多，A、B的变化都比较小。管网的水力计算公式(2)在H0附近可以线性近似为式(7)，且方程有足够的精度(算例的数值计算结果参见附录)。
| 
| (7) |
在文献[2]中已证明B是正定对称矩阵，其逆矩阵存在。 | |
令 | 
| |
则 | 
| (8) |
矩阵C的分量ci,j反映了第j个水源对节点i的影响，矩阵C也称为影响矩阵。如果管网中所有水源的供水水头同步上升Δh，即ΔHN＝[Δh1,Δh2,…,Δhs]T，相当于管网的参考水位提高了Δh。由式(8)知管网中任一点的水头上升的水位 
，因此其中矩阵C的行向量的各分量之和必等于1，各管段的水力坡度不变。如果各水源的供水水头和节点流量已知，可求得管网中的各节点的水头，同样可以求出各水源的供水量。水源泵站供水的动力能耗可以表示为Qs·Hs/γ(γ为水源效率)，供水的动力费用与耗能成正比。水源供水量在H0附近可以线性近似为：Qs=LsHs。水源的制水0费用(除动力费用)可以表示成RsLsHs(Rs表示各水源的单位制水成本)，因此目标函数在H0处可以近似用水源水头的二次函数表示如下：
| (9) |
各水源的供水量约束在H0处可以线性近似表示：
| (10) |
由于在管网中往往只是一部分的最不利节点违反约束，只要最不利的节点满足了服务水头的要求，其他节点也满足了要求，因此可以将最不利的一些节点与水源节点的水头关系从式(9)中的影响矩阵C中抽取出来，表示成矩阵G，管网节点水头的约束方程(5)可以简化表示如下：
| (11) |
其中Ω为最不利节点的集合。
通过上述方法，一级优化模型在H0附近可以近似表示为线性约束的二次规划问题：
| (12a) |
s.t . Hsmin≤Hs≤Hsmax | (12b) |
Qsmin(H0s) +KsminHs≤LsHs≤Qsmax(H0s) +KsmaxHs | (12c) |
| (12d) |
在原优化问题中，各水源的供水量、管网中节点的水头是水源供水水头的函数，是隐式关系，求解起来非常不方便。通过把管网水力关系式(2)线性化，消去原目标函数(1)中的变量——水源供水量Qs，可以把目标函数表示仅含水源水头变量的形式，将管网中各节点的水头HN表示成水源的供水水头Hs的线性函数，只取其中最不利一部分作为每次优化计算的约束条件，这样大大地减少了约束条件。如果管网有上千个节点，只要保证最不利的10%左右节点满足服务水头约束，就能基本上保证每次优化计算结果不会离约束边界太远，同时优化计算的计算量成倍的减少。
由于采用了线性近似的方法简化约束条件和目标函数，采用二次规划法(QP法)优化之后会导致结果越过实际约束边界，其中主要是最不利点不满足管网最小服务水头的要求。在结构优化设计中经常采用齿行法进行优化迭代，其基本思想是在每次优化迭代后，通过射线步(即将所有设计变量以同一倍数放大或缩小)将结果拉到最严格的约束边界上。根据管网水力学，所有的水源的供水水头同时都提高或降低相同的水位，使管网的最不利点的水位恰好处于约束边界上，不会改变各个水源的供水关系。利用这一特性，可以构造一修正的射线步，能够方便的将中间优化迭代点拉回到约束界面上(见图1)。由于在这一修正的射线步中，每个水源提高的水位相同，因此变化后的值与变化前的值与坐标轴正好成45°。上述优化方法的流程图如图2所示。
4 算 例
为了检验上述模型的可行性和可靠性，本文采用Fortran语言编制了优化计算程序，算例中管网基本形状如图3所示。管网有两个水源，12个节点，19根管道，管网节点的基本数据见表1。
在计算中，不失一般性，以泵站供水的能耗作为目标函数，即  (γi是水源泵站的工作效率)，如果考虑取水、净化等制水成本的话，目标增加一个线性项或一个二次项，对本文方法没有影响。收敛准则采用目标值｜objk-objk+1｜/objk≤0.0001和前后两次迭代的步长‖Hk+1-Hk‖≤0.1。为了测试计算方法的有效性和计算效率，本文对多种情况下进行了测试，试验1和试验2是比较在不同效率下，优化获得的目标值是否能够充分反映目标函数对计算结果的影响。由于最大供水量和最大供水扬程往往是制约管网供水能力的约束条件，试验3、试验4是检验最大供水量约束与最大扬程是作用约束时计算结果是否可靠。从表2中的计算结果可以看出计算结果在以上情况下都能获得较满意的结果。
优化迭代计算量的分析：在计算过程中，计算量主要包括以下部分：1)水力计算；2)约束关系的线性化；3)二次规划问题的求解。水力计算是求解非线性方程组，一般是通过线性化，把非线性方程转化为求解线性方程组。在文献[2]中已经证明该线性后方程组的系数矩阵主对角占优，对称正定，可以采用迭代法、高斯消去法、LU分解法或乔斯基分解法进行计算。数值试验表明对于较大规模的管网采用迭代法具有一定的优势。关于约束关系的线性化，如果是采用LU分解法进行的水力计算，只需要进行少量的矩阵变换就可获得约束关系的线性表示。 | 表1管网基本数据 | 编号 | 节点流量/(L/S) | 地面标高/m | 最小服务水头/m | 1
2
3
4
5
6
7
8
9
10
11
12
水源1
水源2 | 36.2
36.8
82.5
36.4
48.7
81.5
198.7
66.1
50.6
43.2
105.8
35.5
—
— | 34.0
35.0
36.0
35.0
37.0
41.0
40.0
43.0
42.0
45.0
46.0
48.0
33
30 | 24.0
24.0
24.0
24.0
24.0
24.0
24.0
24.0
24.0
24.0
24.0
24.0
—
— | |
如果是采用迭代法求解线性方程组，那么还需要额外求解S+1个线性方程组(S表示水源数)，由于水源数一般比较少，约束线性化的计算量相比进行一次水力计算来讲要少得多。再有是二次规划问题的求解，影响二次规划问题的计算量主要是约束量，通过线性化之后，剔除了大部分非作用约束，即使对于上千个节点的管网，也可以把约束量控制在100个左右。在本文算例的迭代计算过程中，服务水头约束条件的最大的违反值是个别最不利点比最小服务水头约束小0.45m，这一误差能够满足实际应用的要求。二次规划问题采用Lemke算法，该算法是在单纯形法基础上加以适当修改来求解二次规划的K-T点，其计算量与单纯行法计算量相当，所以二次规划计算中的计算量也比较小。综合以上的分析，可以估计每次寻优迭代的计算一般不会超过2～3次水力计算，比常用的广义简约梯度法计算量要小得多。本文中的算例特意在配置较低的计算机(PI-120)上进行测试，在以上所有算例的计算时间均少于1s，由此可见本文方法的效率相当高。
表2 优化约束条件与优化结果 |
| 试验1 | 试验2 | 试验3 | 试验4 |
水源1 | 水源2 | 水源1 | 水源2 | 水源1 | 水源2 | 水源1 | 水源2 |
工作效率
最大扬程约束
最小扬程约束
最大供水量约束
最小供水量约束
供水水位
供水量 | 0.7
80
60
0.8
0.02
98.48
0.586 | 0.7
80
60
0.5
0.02
100.95
0.235 | 0.6
80
60
0.8
0.02
96.90
0.573 | 0.7
80
60
0.5
0.02
104.87
0.249 | 0.6
80
60
0.8
0.02
100.43
0.602 | 0.7
80
60
0.22
0.02
96.83
0.220 | 0.6
80
60
0.8
0.02
100.87
0.605 | 0.7
66
60
0.23
0.02
95.99
0.217 |
目标值
迭代次数 | 78.83
11 | 87.76
7 | 88.66
4 | 88.66
4 |
注：试验1和试验2水源水头初值采用附表1中的状态1，试验1和试验2水源水头初值采用附表1中的状态2。 |
5 小 结
供水管网优化调度的一级优化问题是一个复杂的非线性优化问题，针对求解过程中的难点，本文通过分析管网的水力关系，对管网水力关系进行了合理的线性化，使目标函数和约束条件显式化，将问题转化为序列二次规划问题。在求解二次规划问题中，同时考虑到大部分节点水头的约束是非作用约束，利用线性化结果，将非作用约束从约束集中剔除，大大地减少了求解的工作量，并借鉴结构优化中齿行法的思想，结合水力学上的一些基本概念，提出了适合于本文的修正齿行法，将二次规划结果拉回到约束面，保证了解的可行性。最后本文还简要分析了优化迭代的计算量，数值试验表明采用本文的方法计算量小、效率高，结果可靠，为大规模管网的实时优化调度提供了良好的基础。
参 考 文 献：
[1] 汪树玉，杨德铨，刘国华，张科锋.优化原理、优化方法与工程应用[M].杭州：浙江大学出版社，1991，10：375-381，221-228.
[2] 程伟平，包志仁.管网恒定流水力计算的最优化方法[J].水利学报，2000，7：86-90.
[3] 杨钦，严煦世.给水工程[M].北京：中国建筑工业出版社，1987.
附录：
本文在进行线性化推导过程中，认为矩阵C和D在计算过程中变化比较小，本文在计算中选取了供水水头差异较大的情况比较了矩阵C和向量E的值(见附表1)，计算结果表明即使在2种供水状态差异较大的情况下C和向量E的变化比较小。表中 
，与理论分析有一点差别，主要是由计算中截断误差引起。
同时本文还比较了计算了在不同的比例负荷作用下，(即各节点的流量同时乘以一个相同的比例系数)，此时敏度矩阵C矩阵仍然保持基本不变，矩阵假设基本合理。向量E的各系数与比例系数的1852次方成正比，数字计算中间结果见附表2。
附表１ 不同状态下C、E矩阵的变化情况 |
| 状态1
水源1供水水头：100.227m
水源2供水水头： 97.227m | 状态2
水源1供水水头： 98.402m
水源2供水水头：101.131m |
1
2
3
4
5
6
7
8
9
10
11
12
| Ci,1
0.90448
0.90326
0.66842
0.04095
0.94500
0.89536
0.66775
0.37115
0.91821
0.86529
0.73265
0.66456 | Ci,2
0.03544
0.09666
0.33150
0.95897
0.05491
0.10456
0.33217
0.62877
0.08171
0.13463
0.26726
0.33534 | Ei
-4.3374
-9.9901
-22.4775
-2.0184
-6.8753
-11.7309
-24.0387
-17.4159
-11.4742
-16.0033
-26.5912
-27.2211 | Ci,1
0.96663
0.90943
0.69408
0.04145
0.94833
0.90198
0.69326
0.37706
0.92307
0.87392
0.07514
0.69799 | Ci,2
0.03329
0.09048
0.37706
0.95837
0.05158
0.09793
0.30663
0.62293
0.07684
0.12598
0.24847
0.30189 | Ei
-4.3379
-9.9913
-22.4826
-2.0183
-6.8760
-11.7323
-24.0439
-17.4166
-11.4752
-16.0051
-26.5951
-27.2264 |
| | | | | | |
附表2 不同比例负荷状态下C、E矩阵的变化情况 |
| 状态1
水源1供水水头：100.227m，水源2供水水头：97.227m，
负荷比Ｑ1/Q0=1.2，且(Q1/Q0)1.852＝1.4017 | 状态2
水源1供水水头：100.227m，水源2供水水头：97.227m，
负荷比Ｑ2/Q0=1.5，且(Q2/Q0)1.852＝2.1189 |
1
2
3
4
5
6
7
8
9
10
11
12
| Ci,1
0.96482
0.90425
0.67251
0.04103
0.94554
0.89642
0.67181
0.37203
0.91899
0.86668
0.73567
0.66970 | Ci,2
0.03510
0.09567
0.32741
0.95889
0.05439
0.10350
0.32811
0.62789
0.08093
0.13324
0.26425
0.33022 | Ei
-6.0807
-14.006
-31.5198
-2.8293
-9.6385
-16.4462
-33.7080
-24.4138
-16.0853
-22.4356
-37.2818
-38.1722 | Ei/Ei0*
1.4019
1.4020
1.4018
1.4019
1.4019
1.4020
1.4022
1.4018
1.4018
1.4019
1.4020
1.4048 | Ci,1
0.96511
0.90506
0.67588
0.04109
0.94600
0.89727
0.67517
0.37278
0.91963
0.86782
0.73815
0.67400 | Ci,2
0.034815
0.094856
0.32404
0.95883
0.053944
0.10263
0.32475
0.62714
0.080289
0.13210
0.26177
0.32592 | Ei
-9.1932
-21.1757
-47.66.9
4.2773
-14.5721
-24.8649
-50.9687
-36.9088
-24.3185
-33.9200
-56.3678
-57.7025 | Ei/Ei0*
2.1195
2.1197
2.1204
2.1192
2.1151
2.1196
2.1203
2.1193
2.1194
2.1195
2.1198
2.1204 |
注：其中E0i采用附表1中状态1的计算值 |
| | | | | | | | |
综上所述，管网中各节点在比例负荷的情况下，管网中任一节点的水头可以线性近似如下：其中：S为所有水源的集合，Ki为与管网总流量Q0相对应的系数，其中。
收稿日期：2002-01-04
作者简介：刘国华(1963-)，男，福建莆田人，教授，博导，主要从事水工结构优化与系统分析方面的研究。
