趙亞威，方洪斌

（1.黃河勘測規(guī)劃設(shè)計研究院有限公司，河南 鄭州 450003；2.水利部黃河流域水治理與水安全重點實驗室（籌），河南 鄭州 450003）

1 引 言

中長期水量調(diào)度是在中長期時間尺度下（年計劃月調(diào)度或月計劃日調(diào)度），協(xié)調(diào)當(dāng)?shù)氐乇硭?、地下水、非常?guī)水、外調(diào)水等多個水源供水量和用戶需水量的一項水資源優(yōu)化配置工作，是實現(xiàn)水資源系統(tǒng)水量供需平衡和水資源系統(tǒng)綜合效益最大化的重要手段。隨著用水結(jié)構(gòu)、供水來源的多樣化發(fā)展，區(qū)域水資源系統(tǒng)的復(fù)雜性不斷提升，使得中長期水量調(diào)度模型變量和約束增多、規(guī)模變大，求解難度隨之增加，尋求快速、精確的模型求解方法成為目前中長期水量調(diào)度的重要研究方向之一。

國內(nèi)外許多學(xué)者對中長期水量調(diào)度模型求解算法進行了研究，主要分為模擬算法和優(yōu)化算法。模擬算法［1－4］是一種經(jīng)驗方法，主要依靠調(diào)度經(jīng)驗將模型的一些復(fù)雜約束條件概化為對應(yīng)的調(diào)度規(guī)則，進而降低模型的求解難度，實際應(yīng)用較為廣泛，然而這種方法缺乏最優(yōu)化理論支撐，獲取的解往往不是最優(yōu)解；優(yōu)化算法種類較多，主要分為以線性規(guī)劃、動態(tài)規(guī)劃等為代表的傳統(tǒng)優(yōu)化算法［5－8］和以遺傳算法、粒子群算法等為代表的智能優(yōu)化算法［9－10］。傳統(tǒng)優(yōu)化算法比智能優(yōu)化算法求解效率低，例如線性規(guī)劃在求解大規(guī)模復(fù)雜模型時計算量較大、計算時間較長，但線性規(guī)劃全局尋優(yōu)的優(yōu)良特性驅(qū)使國內(nèi)外學(xué)者對提升其求解效率的研究不斷深入。隨著計算機硬件和軟件技術(shù)的不斷突破，計算機整體性能逐步提升，并行計算成為加速線性規(guī)劃計算效率的重要方式，例如目前比較流行的CPU并行、GPU并行可以實現(xiàn)多核、多線程并行計算，能夠大幅提升計算效率。針對線性規(guī)劃求解效率低的問題，許多學(xué)者和軟件廠商基于并行技術(shù)開發(fā)了線性規(guī)劃軟件。線性求解器是專門用于求解線性規(guī)劃問題的軟件或程序集（大部分求解器還可求解非線性規(guī)劃問題，如二次規(guī)劃），目前國際上流行的線性求解器主要有Cplex、Gurobi、Xpress、Lingo／Lindo等，免費開源求解器主要有Lpsolver、Scip、Glpk、Matlab等；國內(nèi)線性求解器發(fā)展起步較晚，目前應(yīng)用較為廣泛的主要有Cmip、Leaves等。經(jīng)過大量實例驗證，線性求解器在求解大規(guī)模復(fù)雜線性模型時顯現(xiàn)出優(yōu)秀的計算性能。

本文將線性求解器運用于中長期水量調(diào)度模型求解中，并將其與遺傳算法進行對比，以期探索模型的快速求解技術(shù)，為中長期水量調(diào)度方案編制提供參考。

2 中長期水量調(diào)度模型

2.1 水資源系統(tǒng)網(wǎng)絡(luò)拓?fù)?/h3>

中長期水量調(diào)度模型一般基于網(wǎng)絡(luò)拓?fù)錁?gòu)建，網(wǎng)絡(luò)拓?fù)溆晒?jié)點、連線構(gòu)成，節(jié)點表示輸入水量、輸出水量、蓄水量等，連線表示水流方向。對于水資源系統(tǒng)，節(jié)點可對應(yīng)水庫、配水樞紐、控制斷面、水廠、用水戶等，其中水庫有蓄水庫容、配水樞紐有調(diào)節(jié)池，均具有一定調(diào)蓄能力，故蓄水量可變化；連線可對應(yīng)河流、渠道、隧洞、管線等，在進行中長期水量調(diào)度時，一般將連線蓄水量設(shè)為恒定值或0。節(jié)點和連線遵循水量平衡原則，構(gòu)成供需水網(wǎng)絡(luò)。節(jié)點、連線水量平衡示意圖分別見圖1、圖2。圖1中，R上游、R區(qū)間、R退水、R地表、R地下、R其他分別為節(jié)點i的上游相鄰節(jié)點出流水量、區(qū)間水量、退水量、當(dāng)?shù)氐乇硭?、地下水量、其他水量（再生水、收集雨水、淡化水等），均為?jié)點i的輸入水量；C下游、C退水、C損失、C生活、C工業(yè)、C農(nóng)業(yè)、C生態(tài)分別為節(jié)點i的下游相鄰節(jié)點入流水量、退水量、損失水量、生活供水量、工業(yè)供水量、農(nóng)業(yè)供水量、生態(tài)供水量，均為節(jié)點i的輸出水量；S時段初、S時段末分別為節(jié)點i時段初、末蓄水量。圖2中，R節(jié)點i為連線ij上游節(jié)點i的出流水量（輸入水量）；C損失、C節(jié)點j分別為連線ij的損失水量、下游節(jié)點j的入流水量（輸出水量）；連線ij的箭頭表示水流方向。

圖1 節(jié)點水量平衡示意

圖2 連線水量平衡示意

2.2 基于網(wǎng)絡(luò)拓?fù)涞乃空{(diào)度模型

網(wǎng)絡(luò)拓?fù)涫敲枋鏊Y源系統(tǒng)各單元水量供需關(guān)系的重要工具，依據(jù)網(wǎng)絡(luò)拓?fù)涔?jié)點和連線的上下游關(guān)系、水量平衡約束、其他約束條件以及調(diào)度目標(biāo)可構(gòu)建中長期水量調(diào)度模型。

2.2.1 目標(biāo)函數(shù)

以水資源系統(tǒng)總?cè)彼孔钚閮?yōu)化目標(biāo)，目標(biāo)函數(shù)表達式為

式中：I為節(jié)點總數(shù)；J、j為節(jié)點上用水戶總數(shù)、序號，j取值為1、2、3、4，分別表示生活、工業(yè)、農(nóng)業(yè)、生態(tài)用水戶；T為調(diào)度時段總數(shù)；Qijt、Cijt分別為節(jié)點i用戶j時段t的需水量和實際供水量。

2.2.2 約束條件

（1）節(jié)點水量平衡。表達式為

其中

式中：Ni、n分別為節(jié)點i上游相鄰節(jié)點總數(shù)、序號；Mi、m分別為節(jié)點i下游相鄰節(jié)點總數(shù)、序號；Rsyint、Rqint、Rtsint及ain、bin、fin分別為節(jié)點i時段t的上游第n個相鄰節(jié)點出流水量、區(qū)間入流水量、退水量及對應(yīng)的沿程損耗系數(shù)；Cxyimt、Ctsimt、Cgsit、Cdcit分別為節(jié)點i時段t的下游第m個相鄰節(jié)點出流水量、退水量、總供水量、調(diào)出到外系統(tǒng)水量；Vi（t－1）、Vit分別為節(jié)點i時段t初、末蓄水量；cij為節(jié)點i用戶j的退水系數(shù)；dim為節(jié)點i與下游相鄰節(jié)點m之間的退水量占節(jié)點i總退水量的比例；ei為節(jié)點i上供水管網(wǎng)損耗系數(shù)；K、k分別為節(jié)點上水源數(shù)量、序號，其中k取值為1、2、3、4，分別表示當(dāng)?shù)氐乇硭?、地下水、非常?guī)水和過境水（上游相鄰節(jié)點出流水量、區(qū)間水量、退水量之和）；Dijkt為時段t水源k向節(jié)點i用戶j的實際供水量。

（2）連線水量平衡。若網(wǎng)絡(luò)拓?fù)渲泄?jié)點i與節(jié)點p直接相連，水流由節(jié)點i流向節(jié)點p，則連線ip滿足如下水量平衡方程：

式中：Cxyivt為時段t從節(jié)點i流向節(jié)點p的出流水量；Rsyput為時段t由節(jié)點i流入節(jié)點p的水量；v為節(jié)點p在節(jié)點i下游相鄰節(jié)點中的序號，1≤v≤Mi；u為節(jié)點i在節(jié)點p上游相鄰節(jié)點中的序號，1≤u≤Ni。

（3）節(jié)點出流上限、下限約束。約束條件為

式中：、分別為節(jié)點i下游相鄰的第m個節(jié)點時段t的出流量下限、上限，為節(jié)點的過流能力，一般取0，若當(dāng)前節(jié)點存在生態(tài)流量約束，則取生態(tài)流量。

（4）節(jié)點蓄水量上限、下限約束。約束條件為

（5）節(jié)點各水源向各用戶的實際供水量上限、下限約束。約束條件為

式中：、分別為時段t節(jié)點i上水源k向用戶j的實際供水量下限、上限。

3 基于線性求解器的模型求解

本文所構(gòu)建的中長期水量調(diào)度模型中不含非線性表達式，為標(biāo)準(zhǔn)的線性規(guī)劃模型，可采用線性求解器進行求解。線性求解器的種類較多，但主要的數(shù)學(xué)原理、求解方法、操作流程有相似之處。數(shù)學(xué)原理方面，大多數(shù)求解器是基于單純形法、分支定界法等傳統(tǒng)的線性規(guī)劃方法；求解方法方面，大多數(shù)求解器是輸入標(biāo)準(zhǔn)化線性規(guī)劃模型的目標(biāo)函數(shù)系數(shù)向量、約束條件系數(shù)矩陣和常數(shù)向量、目標(biāo)優(yōu)化方式Max／Min、決策變量個數(shù)、整數(shù)變量個數(shù)等信息，然后進行求解，最后得到最優(yōu)決策向量和目標(biāo)值；操作流程方面，使用線性求解器求解模型一般需要編寫程序，主要有集成開發(fā)環(huán)境（IDE）和基于應(yīng)用程序編程接口（API）兩種方式。IDE的優(yōu)點是程序編寫方式與數(shù)學(xué)模型的書寫習(xí)慣較為統(tǒng)一，程序結(jié)構(gòu)簡單；缺點是不同求解器的程序語法往往不同，在使用多種求解器時需要熟悉多種IDE程序語法，另外IDE一般不能與其他語言進行程序混編，無法集成到C／C＋＋、Java、Python等其他應(yīng)用程序開發(fā)系統(tǒng)中。API優(yōu)點是支持C／C＋＋、Java、Python等當(dāng)前主流編程語言，程序語法易于學(xué)習(xí)，調(diào)用方便，便于進行多語言混編，可用于大型應(yīng)用系統(tǒng)的優(yōu)化計算和模塊開發(fā)；缺點是很多求解器不提供API，有的僅支持部分主流編程語言，且API功能相較于IDE相對單薄。

本文采用C＃與Cplex、Lpsolver線性求解器API混合編程實現(xiàn)模型求解，步驟如下。

（1）模型前處理。將模型進行標(biāo)準(zhǔn)化處理，包括將目標(biāo)函數(shù)轉(zhuǎn)化為“Min”形式以及將所有不等式約束條件轉(zhuǎn)化為“≤”形式。

（2）模型關(guān)鍵參數(shù)提取。提取中長期水量調(diào)度數(shù)學(xué)模型中的關(guān)鍵信息，包括目標(biāo)函數(shù)中決策變量系數(shù)向量C、約束條件中決策變量系數(shù)矩陣A和常數(shù)向量b、決策變量上下限、決策變量為整數(shù)或非整數(shù)等信息，形成矩陣數(shù)據(jù)或向量數(shù)據(jù)。

（3）模型求解。將步驟（2）提取的所有數(shù)學(xué)模型矩陣數(shù)據(jù)和向量數(shù)據(jù)，整理成Cplex或Lpsolver的API標(biāo)準(zhǔn)化輸入形式，調(diào)用求解器中的線性規(guī)劃方法對模型進行求解。

（4）結(jié)果輸出。按照Cplex或Lpsolver輸出的規(guī)則，將計算得到的最優(yōu)結(jié)果解析為方便進行方案分析的數(shù)據(jù)格式。

小兒柴桂退熱顆粒的UPLC指紋圖譜及聚類、主成分分析…………………………………………………… 林 源等（4）：474

4 實例分析

以北方某城市2025年規(guī)劃水平年中長期水量調(diào)度為例進行分析。該市境內(nèi)主要河流有5條，分別命名為河流A、B、C、D、E。該市現(xiàn)狀供水水源主要由當(dāng)?shù)氐乇硭?、地下水、非常?guī)水和從A河引水組成。2019年該市從A河引水量8.7億m3，占該市總用水量的38.3%。從A河引水進入該市水資源系統(tǒng)后，主要由位于河流B的水庫1和位于河流D的水庫2進行調(diào)節(jié)，供給該市生活、工業(yè)、農(nóng)業(yè)、生態(tài)環(huán)境需水。該市水資源系統(tǒng)概化見圖3，網(wǎng)絡(luò)拓?fù)涓呕妶D4。采用Cplex、Lpsolver、GA（遺傳算法）對模型進行求解。

圖3 某北方城市水資源系統(tǒng)概化

圖4 某北方城市水資源系統(tǒng)網(wǎng)絡(luò)拓?fù)涓呕?/p>

4.1 不同方案優(yōu)化結(jié)果分析

該市每年各月按照一定比例（見表1）從河流A的控制水庫引水，按照河流A來水量豐枯變化，分別設(shè)置從A河年引水量4億、5億、6億、7億、8億m35種方案，優(yōu)化結(jié)果見表2。

表1 某北方城市各月從河流A的引水量比例

表2 從A河引水不同方案優(yōu)化結(jié)果

從優(yōu)化結(jié)果來看，各方案Cplex、Lpsolver優(yōu)化得到的缺水量和缺水率均小于GA的優(yōu)化結(jié)果，原因是Cplex、Lpsolver是基于線性規(guī)劃的全局優(yōu)化算法，而GA則存在易陷入局部最優(yōu)的缺點，無法保證收斂到全局最優(yōu)解。從計算時間來看，采用Cplex、Lpsolver的各方案計算時間均小于0.1 s，而采用GA的計算時間最短需978.89 s，最長則需1 188.83 s，Cplex、Lpsolver計算時間遠遠少于GA，原因是Cplex、Lpsolver采用了并行計算，充分發(fā)揮了計算機性能，極大提升了計算效率，而GA不僅存在繁瑣的選擇、交叉、變異和迭代等過程，且未采用并行計算，計算效率較低。

4.2 不同規(guī)模模型優(yōu)化結(jié)果分析

為了進一步檢驗不同算法的性能，改變模型規(guī)模進行對比分析。以圖4模型19個節(jié)點（大規(guī)模）為基礎(chǔ)，將節(jié)點數(shù)精簡為8個（小規(guī)模），分別采用Cplex、Lpsolver、GA進行優(yōu)化求解，結(jié)果見表3。

表3 小規(guī)模模型從A河引水不同方案優(yōu)化結(jié)果

通過對比表2、表3可知，由于Cplex、Lpsolver求解效率較高，因此模型規(guī)模大小對其計算時間影響不大；小規(guī)模模型和大規(guī)模模型不同A河引水方案下，GA的平均計算時長分別為111.70、1 107.84 s，GA限于自身算法的復(fù)雜性和未采用并行機制，隨著模型規(guī)模的增大，計算時間增加。

4.3 從A河引水時機分析

由表2、表3可知，從A河引水8億m3時，仍不能滿足總水量需求，原因是各月從A河引水比例與該市需水比例變化趨勢不完全一致，而該市境內(nèi)主要依靠水庫1、水庫2來調(diào)節(jié)從A河引水分配過程，這兩座水庫興利庫容分別為0.651億m3和0.684億m3，相較于從A河引水量4億～8億m3來說均很小，調(diào)蓄能力非常有限，無法進行時空調(diào)節(jié)。為解決這一問題，可以從增加水庫興利庫容或者改變引水時機兩方面入手，但該市近期沒有新建水庫或增加現(xiàn)有水庫興利庫容的規(guī)劃，因此只能從改變引水時機角度考慮改善缺水狀況。

在實際調(diào)度中，該市年內(nèi)各月從A河的引水量為固定值，可表示為

式中：qt、q0t分別為第t月該市從A河的引水量、固定引水量。

對原模型進行改進，將“各月從A河引水量為固定值”等式約束修改為“各月從A河引水量之和為固定值”等式約束，計算公式為

式中：Q為從A河年引水總量。Q與q0t的關(guān)系為

以16個節(jié)點模型為例，以從A河引水量4億m3為起點、8億m3為終點，等值離散為若干個方案，采用Cplex算法對不同引水量下不考慮引水時機（原模型）和考慮引水時機（改進模型）兩種模型進行求解，得到從A河引水量與缺水量的關(guān)系，見圖5。

圖5 從A河引水量與缺水量的關(guān)系

（1）不同引水方案下，改進模型的缺水量均小于原模型，且隨著從A河引水量的增加，改進模型較原模型缺水量逐漸減小，說明改變引水時機能夠減小從A河引水量，節(jié)約系統(tǒng)供水成本。

（2）當(dāng)從A河引水量為7.65億m3時，改進模型的缺水量為0，該引水量下水庫1、水庫2的水位變化過程見圖6和圖7，可以看出，原模型水庫1僅在2月、9月發(fā)揮了調(diào)蓄作用，水庫2僅在5月、8月、9月發(fā)揮了調(diào)蓄作用，而改進模型的水庫1和水庫2的蓄泄變化過程更加明顯，說明改進模型充分發(fā)揮了水庫的調(diào)蓄作用。原模型、改進模型從A河引水過程見圖8，可知：7月、9—12月當(dāng)?shù)厮枯^富裕，改進模型比原模型減小了引水量，而在1—4月、6月加大了從A河引水量，彌補需水缺口；5月需水量較大，而改進模型從A河引水量小于原模型，原因是改進模型1—4月從A河引水量較大，增加了水庫1、水庫2的蓄水量（5月初水庫1、水庫2蓄水量較大），5月從A河引水，水庫1和水庫2同時向用戶供水，故改進模型5月從A河引水量小于原模型；8月當(dāng)?shù)厮糠浅８辉＃Ｐ蛷腁河引水量較小，而改進模型引水量較大，原因是8月當(dāng)?shù)厮饕獊碜怨?jié)點N1002、N1004、N1006、N2001、N2002、N3002、N4001、N5002、N5003，除N2001（水 庫1）、N4001（水庫2）外，其他節(jié)點當(dāng)?shù)厮浚?.35億m3，占8月當(dāng)?shù)厮偭康?2.76%）均無法存入水庫1和水庫2，不能被調(diào)蓄，因此需要加大從A河引水量，來滿足9—12月用戶需求。從圖6、圖7可以看出，8月末水庫1和水庫2蓄水量明顯增加，且9—12月為開閘供水狀態(tài)，水庫水位不斷降低，直至水庫放空。綜上所述，與原模型相比，改進模型充分發(fā)揮了當(dāng)?shù)厮畮斓恼{(diào)蓄作用。

圖6 水庫1水位變化過程

圖7 水庫2水位變化過程

圖8 原模型、改進模型從A河引水過程

5 結(jié) 語

本文以北方某城市水資源系統(tǒng)為例，建立中長期水量調(diào)度模型，采用線性求解器對模型進行求解，并與遺傳算法進行了對比。主要結(jié)論如下：①線性求解器Cplex、Lpsolver編程簡單，能夠獲得全局最優(yōu)解，求解精度高、速度快；②當(dāng)模型規(guī)模增大時，遺傳算法計算效率明顯下降，而Cplex、Lpsolver計算效率變化不明顯；③通過改變引水時機能夠充分發(fā)揮現(xiàn)有水庫的調(diào)蓄能力，增加系統(tǒng)總供水量。

線性求解器應(yīng)用的前提是模型必須為線性模型，因此如何實現(xiàn)含非線性約束的中長期水量調(diào)度模型的線性化處理，是下一步的研究重點。