趙亞威,方洪斌
(1.黃河勘測(cè)規(guī)劃設(shè)計(jì)研究院有限公司,河南 鄭州 450003;2.水利部黃河流域水治理與水安全重點(diǎn)實(shí)驗(yàn)室(籌),河南 鄭州 450003)
中長(zhǎng)期水量調(diào)度是在中長(zhǎng)期時(shí)間尺度下(年計(jì)劃月調(diào)度或月計(jì)劃日調(diào)度),協(xié)調(diào)當(dāng)?shù)氐乇硭?、地下水、非常?guī)水、外調(diào)水等多個(gè)水源供水量和用戶(hù)需水量的一項(xiàng)水資源優(yōu)化配置工作,是實(shí)現(xiàn)水資源系統(tǒng)水量供需平衡和水資源系統(tǒng)綜合效益最大化的重要手段。隨著用水結(jié)構(gòu)、供水來(lái)源的多樣化發(fā)展,區(qū)域水資源系統(tǒng)的復(fù)雜性不斷提升,使得中長(zhǎng)期水量調(diào)度模型變量和約束增多、規(guī)模變大,求解難度隨之增加,尋求快速、精確的模型求解方法成為目前中長(zhǎng)期水量調(diào)度的重要研究方向之一。
國(guó)內(nèi)外許多學(xué)者對(duì)中長(zhǎng)期水量調(diào)度模型求解算法進(jìn)行了研究,主要分為模擬算法和優(yōu)化算法。模擬算法[1-4]是一種經(jīng)驗(yàn)方法,主要依靠調(diào)度經(jīng)驗(yàn)將模型的一些復(fù)雜約束條件概化為對(duì)應(yīng)的調(diào)度規(guī)則,進(jìn)而降低模型的求解難度,實(shí)際應(yīng)用較為廣泛,然而這種方法缺乏最優(yōu)化理論支撐,獲取的解往往不是最優(yōu)解;優(yōu)化算法種類(lèi)較多,主要分為以線性規(guī)劃、動(dòng)態(tài)規(guī)劃等為代表的傳統(tǒng)優(yōu)化算法[5-8]和以遺傳算法、粒子群算法等為代表的智能優(yōu)化算法[9-10]。傳統(tǒng)優(yōu)化算法比智能優(yōu)化算法求解效率低,例如線性規(guī)劃在求解大規(guī)模復(fù)雜模型時(shí)計(jì)算量較大、計(jì)算時(shí)間較長(zhǎng),但線性規(guī)劃全局尋優(yōu)的優(yōu)良特性驅(qū)使國(guó)內(nèi)外學(xué)者對(duì)提升其求解效率的研究不斷深入。隨著計(jì)算機(jī)硬件和軟件技術(shù)的不斷突破,計(jì)算機(jī)整體性能逐步提升,并行計(jì)算成為加速線性規(guī)劃計(jì)算效率的重要方式,例如目前比較流行的CPU并行、GPU并行可以實(shí)現(xiàn)多核、多線程并行計(jì)算,能夠大幅提升計(jì)算效率。針對(duì)線性規(guī)劃求解效率低的問(wèn)題,許多學(xué)者和軟件廠商基于并行技術(shù)開(kāi)發(fā)了線性規(guī)劃軟件。線性求解器是專(zhuān)門(mén)用于求解線性規(guī)劃問(wèn)題的軟件或程序集(大部分求解器還可求解非線性規(guī)劃問(wèn)題,如二次規(guī)劃),目前國(guó)際上流行的線性求解器主要有Cplex、Gurobi、Xpress、Lingo/Lindo等,免費(fèi)開(kāi)源求解器主要有Lpsolver、Scip、Glpk、Matlab等;國(guó)內(nèi)線性求解器發(fā)展起步較晚,目前應(yīng)用較為廣泛的主要有Cmip、Leaves等。經(jīng)過(guò)大量實(shí)例驗(yàn)證,線性求解器在求解大規(guī)模復(fù)雜線性模型時(shí)顯現(xiàn)出優(yōu)秀的計(jì)算性能。
本文將線性求解器運(yùn)用于中長(zhǎng)期水量調(diào)度模型求解中,并將其與遺傳算法進(jìn)行對(duì)比,以期探索模型的快速求解技術(shù),為中長(zhǎng)期水量調(diào)度方案編制提供參考。
中長(zhǎng)期水量調(diào)度模型一般基于網(wǎng)絡(luò)拓?fù)錁?gòu)建,網(wǎng)絡(luò)拓?fù)溆晒?jié)點(diǎn)、連線構(gòu)成,節(jié)點(diǎn)表示輸入水量、輸出水量、蓄水量等,連線表示水流方向。對(duì)于水資源系統(tǒng),節(jié)點(diǎn)可對(duì)應(yīng)水庫(kù)、配水樞紐、控制斷面、水廠、用水戶(hù)等,其中水庫(kù)有蓄水庫(kù)容、配水樞紐有調(diào)節(jié)池,均具有一定調(diào)蓄能力,故蓄水量可變化;連線可對(duì)應(yīng)河流、渠道、隧洞、管線等,在進(jìn)行中長(zhǎng)期水量調(diào)度時(shí),一般將連線蓄水量設(shè)為恒定值或0。節(jié)點(diǎn)和連線遵循水量平衡原則,構(gòu)成供需水網(wǎng)絡(luò)。節(jié)點(diǎn)、連線水量平衡示意圖分別見(jiàn)圖1、圖2。圖1中,R上游、R區(qū)間、R退水、R地表、R地下、R其他分別為節(jié)點(diǎn)i的上游相鄰節(jié)點(diǎn)出流水量、區(qū)間水量、退水量、當(dāng)?shù)氐乇硭?、地下水量、其他水量(再生水、收集雨水、淡化水等),均為?jié)點(diǎn)i的輸入水量;C下游、C退水、C損失、C生活、C工業(yè)、C農(nóng)業(yè)、C生態(tài)分別為節(jié)點(diǎn)i的下游相鄰節(jié)點(diǎn)入流水量、退水量、損失水量、生活供水量、工業(yè)供水量、農(nóng)業(yè)供水量、生態(tài)供水量,均為節(jié)點(diǎn)i的輸出水量;S時(shí)段初、S時(shí)段末分別為節(jié)點(diǎn)i時(shí)段初、末蓄水量。圖2中,R節(jié)點(diǎn)i為連線ij上游節(jié)點(diǎn)i的出流水量(輸入水量);C損失、C節(jié)點(diǎn)j分別為連線ij的損失水量、下游節(jié)點(diǎn)j的入流水量(輸出水量);連線ij的箭頭表示水流方向。
圖1 節(jié)點(diǎn)水量平衡示意
圖2 連線水量平衡示意
網(wǎng)絡(luò)拓?fù)涫敲枋鏊Y源系統(tǒng)各單元水量供需關(guān)系的重要工具,依據(jù)網(wǎng)絡(luò)拓?fù)涔?jié)點(diǎn)和連線的上下游關(guān)系、水量平衡約束、其他約束條件以及調(diào)度目標(biāo)可構(gòu)建中長(zhǎng)期水量調(diào)度模型。
2.2.1 目標(biāo)函數(shù)
以水資源系統(tǒng)總?cè)彼孔钚閮?yōu)化目標(biāo),目標(biāo)函數(shù)表達(dá)式為
式中:I為節(jié)點(diǎn)總數(shù);J、j為節(jié)點(diǎn)上用水戶(hù)總數(shù)、序號(hào),j取值為1、2、3、4,分別表示生活、工業(yè)、農(nóng)業(yè)、生態(tài)用水戶(hù);T為調(diào)度時(shí)段總數(shù);Qijt、Cijt分別為節(jié)點(diǎn)i用戶(hù)j時(shí)段t的需水量和實(shí)際供水量。
2.2.2 約束條件
(1)節(jié)點(diǎn)水量平衡。表達(dá)式為
其中
式中:Ni、n分別為節(jié)點(diǎn)i上游相鄰節(jié)點(diǎn)總數(shù)、序號(hào);Mi、m分別為節(jié)點(diǎn)i下游相鄰節(jié)點(diǎn)總數(shù)、序號(hào);Rsyint、Rqint、Rtsint及ain、bin、fin分別為節(jié)點(diǎn)i時(shí)段t的上游第n個(gè)相鄰節(jié)點(diǎn)出流水量、區(qū)間入流水量、退水量及對(duì)應(yīng)的沿程損耗系數(shù);Cxyimt、Ctsimt、Cgsit、Cdcit分別為節(jié)點(diǎn)i時(shí)段t的下游第m個(gè)相鄰節(jié)點(diǎn)出流水量、退水量、總供水量、調(diào)出到外系統(tǒng)水量;Vi(t-1)、Vit分別為節(jié)點(diǎn)i時(shí)段t初、末蓄水量;cij為節(jié)點(diǎn)i用戶(hù)j的退水系數(shù);dim為節(jié)點(diǎn)i與下游相鄰節(jié)點(diǎn)m之間的退水量占節(jié)點(diǎn)i總退水量的比例;ei為節(jié)點(diǎn)i上供水管網(wǎng)損耗系數(shù);K、k分別為節(jié)點(diǎn)上水源數(shù)量、序號(hào),其中k取值為1、2、3、4,分別表示當(dāng)?shù)氐乇硭?、地下水、非常?guī)水和過(guò)境水(上游相鄰節(jié)點(diǎn)出流水量、區(qū)間水量、退水量之和);Dijkt為時(shí)段t水源k向節(jié)點(diǎn)i用戶(hù)j的實(shí)際供水量。
(2)連線水量平衡。若網(wǎng)絡(luò)拓?fù)渲泄?jié)點(diǎn)i與節(jié)點(diǎn)p直接相連,水流由節(jié)點(diǎn)i流向節(jié)點(diǎn)p,則連線ip滿(mǎn)足如下水量平衡方程:
式中:Cxyivt為時(shí)段t從節(jié)點(diǎn)i流向節(jié)點(diǎn)p的出流水量;Rsyput為時(shí)段t由節(jié)點(diǎn)i流入節(jié)點(diǎn)p的水量;v為節(jié)點(diǎn)p在節(jié)點(diǎn)i下游相鄰節(jié)點(diǎn)中的序號(hào),1≤v≤Mi;u為節(jié)點(diǎn)i在節(jié)點(diǎn)p上游相鄰節(jié)點(diǎn)中的序號(hào),1≤u≤Ni。
(3)節(jié)點(diǎn)出流上限、下限約束。約束條件為
式中:、分別為節(jié)點(diǎn)i下游相鄰的第m個(gè)節(jié)點(diǎn)時(shí)段t的出流量下限、上限,為節(jié)點(diǎn)的過(guò)流能力,一般取0,若當(dāng)前節(jié)點(diǎn)存在生態(tài)流量約束,則取生態(tài)流量。
(4)節(jié)點(diǎn)蓄水量上限、下限約束。約束條件為
(5)節(jié)點(diǎn)各水源向各用戶(hù)的實(shí)際供水量上限、下限約束。約束條件為
式中:、分別為時(shí)段t節(jié)點(diǎn)i上水源k向用戶(hù)j的實(shí)際供水量下限、上限。
本文所構(gòu)建的中長(zhǎng)期水量調(diào)度模型中不含非線性表達(dá)式,為標(biāo)準(zhǔn)的線性規(guī)劃模型,可采用線性求解器進(jìn)行求解。線性求解器的種類(lèi)較多,但主要的數(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、決策變量個(gè)數(shù)、整數(shù)變量個(gè)數(shù)等信息,然后進(jìn)行求解,最后得到最優(yōu)決策向量和目標(biāo)值;操作流程方面,使用線性求解器求解模型一般需要編寫(xiě)程序,主要有集成開(kāi)發(fā)環(huán)境(IDE)和基于應(yīng)用程序編程接口(API)兩種方式。IDE的優(yōu)點(diǎn)是程序編寫(xiě)方式與數(shù)學(xué)模型的書(shū)寫(xiě)習(xí)慣較為統(tǒng)一,程序結(jié)構(gòu)簡(jiǎn)單;缺點(diǎn)是不同求解器的程序語(yǔ)法往往不同,在使用多種求解器時(shí)需要熟悉多種IDE程序語(yǔ)法,另外IDE一般不能與其他語(yǔ)言進(jìn)行程序混編,無(wú)法集成到C/C++、Java、Python等其他應(yīng)用程序開(kāi)發(fā)系統(tǒng)中。API優(yōu)點(diǎn)是支持C/C++、Java、Python等當(dāng)前主流編程語(yǔ)言,程序語(yǔ)法易于學(xué)習(xí),調(diào)用方便,便于進(jìn)行多語(yǔ)言混編,可用于大型應(yīng)用系統(tǒng)的優(yōu)化計(jì)算和模塊開(kāi)發(fā);缺點(diǎn)是很多求解器不提供API,有的僅支持部分主流編程語(yǔ)言,且API功能相較于IDE相對(duì)單薄。
本文采用C#與Cplex、Lpsolver線性求解器API混合編程實(shí)現(xiàn)模型求解,步驟如下。
(1)模型前處理。將模型進(jìn)行標(biāo)準(zhǔn)化處理,包括將目標(biāo)函數(shù)轉(zhuǎn)化為“Min”形式以及將所有不等式約束條件轉(zhuǎn)化為“≤”形式。
(2)模型關(guān)鍵參數(shù)提取。提取中長(zhǎng)期水量調(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ī)劃方法對(duì)模型進(jìn)行求解。
(4)結(jié)果輸出。按照Cplex或Lpsolver輸出的規(guī)則,將計(jì)算得到的最優(yōu)結(jié)果解析為方便進(jìn)行方案分析的數(shù)據(jù)格式。
小兒柴桂退熱顆粒的UPLC指紋圖譜及聚類(lèi)、主成分分析…………………………………………………… 林 源等(4):474
以北方某城市2025年規(guī)劃水平年中長(zhǎng)期水量調(diào)度為例進(jìn)行分析。該市境內(nèi)主要河流有5條,分別命名為河流A、B、C、D、E。該市現(xiàn)狀供水水源主要由當(dāng)?shù)氐乇硭?、地下水、非常?guī)水和從A河引水組成。2019年該市從A河引水量8.7億m3,占該市總用水量的38.3%。從A河引水進(jìn)入該市水資源系統(tǒng)后,主要由位于河流B的水庫(kù)1和位于河流D的水庫(kù)2進(jìn)行調(diào)節(jié),供給該市生活、工業(yè)、農(nóng)業(yè)、生態(tài)環(huán)境需水。該市水資源系統(tǒng)概化見(jiàn)圖3,網(wǎng)絡(luò)拓?fù)涓呕?jiàn)圖4。采用Cplex、Lpsolver、GA(遺傳算法)對(duì)模型進(jìn)行求解。
圖3 某北方城市水資源系統(tǒng)概化
圖4 某北方城市水資源系統(tǒng)網(wǎng)絡(luò)拓?fù)涓呕?/p>
該市每年各月按照一定比例(見(jiàn)表1)從河流A的控制水庫(kù)引水,按照河流A來(lái)水量豐枯變化,分別設(shè)置從A河年引水量4億、5億、6億、7億、8億m35種方案,優(yōu)化結(jié)果見(jiàn)表2。
表1 某北方城市各月從河流A的引水量比例
表2 從A河引水不同方案優(yōu)化結(jié)果
從優(yōu)化結(jié)果來(lái)看,各方案Cplex、Lpsolver優(yōu)化得到的缺水量和缺水率均小于GA的優(yōu)化結(jié)果,原因是Cplex、Lpsolver是基于線性規(guī)劃的全局優(yōu)化算法,而GA則存在易陷入局部最優(yōu)的缺點(diǎn),無(wú)法保證收斂到全局最優(yōu)解。從計(jì)算時(shí)間來(lái)看,采用Cplex、Lpsolver的各方案計(jì)算時(shí)間均小于0.1 s,而采用GA的計(jì)算時(shí)間最短需978.89 s,最長(zhǎng)則需1 188.83 s,Cplex、Lpsolver計(jì)算時(shí)間遠(yuǎn)遠(yuǎn)少于GA,原因是Cplex、Lpsolver采用了并行計(jì)算,充分發(fā)揮了計(jì)算機(jī)性能,極大提升了計(jì)算效率,而GA不僅存在繁瑣的選擇、交叉、變異和迭代等過(guò)程,且未采用并行計(jì)算,計(jì)算效率較低。
為了進(jìn)一步檢驗(yàn)不同算法的性能,改變模型規(guī)模進(jìn)行對(duì)比分析。以圖4模型19個(gè)節(jié)點(diǎn)(大規(guī)模)為基礎(chǔ),將節(jié)點(diǎn)數(shù)精簡(jiǎn)為8個(gè)(小規(guī)模),分別采用Cplex、Lpsolver、GA進(jìn)行優(yōu)化求解,結(jié)果見(jiàn)表3。
表3 小規(guī)模模型從A河引水不同方案優(yōu)化結(jié)果
通過(guò)對(duì)比表2、表3可知,由于Cplex、Lpsolver求解效率較高,因此模型規(guī)模大小對(duì)其計(jì)算時(shí)間影響不大;小規(guī)模模型和大規(guī)模模型不同A河引水方案下,GA的平均計(jì)算時(shí)長(zhǎng)分別為111.70、1 107.84 s,GA限于自身算法的復(fù)雜性和未采用并行機(jī)制,隨著模型規(guī)模的增大,計(jì)算時(shí)間增加。
由表2、表3可知,從A河引水8億m3時(shí),仍不能滿(mǎn)足總水量需求,原因是各月從A河引水比例與該市需水比例變化趨勢(shì)不完全一致,而該市境內(nèi)主要依靠水庫(kù)1、水庫(kù)2來(lái)調(diào)節(jié)從A河引水分配過(guò)程,這兩座水庫(kù)興利庫(kù)容分別為0.651億m3和0.684億m3,相較于從A河引水量4億~8億m3來(lái)說(shuō)均很小,調(diào)蓄能力非常有限,無(wú)法進(jìn)行時(shí)空調(diào)節(jié)。為解決這一問(wèn)題,可以從增加水庫(kù)興利庫(kù)容或者改變引水時(shí)機(jī)兩方面入手,但該市近期沒(méi)有新建水庫(kù)或增加現(xiàn)有水庫(kù)興利庫(kù)容的規(guī)劃,因此只能從改變引水時(shí)機(jī)角度考慮改善缺水狀況。
在實(shí)際調(diào)度中,該市年內(nèi)各月從A河的引水量為固定值,可表示為
式中:qt、q0t分別為第t月該市從A河的引水量、固定引水量。
對(duì)原模型進(jìn)行改進(jìn),將“各月從A河引水量為固定值”等式約束修改為“各月從A河引水量之和為固定值”等式約束,計(jì)算公式為
式中:Q為從A河年引水總量。Q與q0t的關(guān)系為
以16個(gè)節(jié)點(diǎn)模型為例,以從A河引水量4億m3為起點(diǎn)、8億m3為終點(diǎn),等值離散為若干個(gè)方案,采用Cplex算法對(duì)不同引水量下不考慮引水時(shí)機(jī)(原模型)和考慮引水時(shí)機(jī)(改進(jìn)模型)兩種模型進(jìn)行求解,得到從A河引水量與缺水量的關(guān)系,見(jiàn)圖5。
圖5 從A河引水量與缺水量的關(guān)系
(1)不同引水方案下,改進(jìn)模型的缺水量均小于原模型,且隨著從A河引水量的增加,改進(jìn)模型較原模型缺水量逐漸減小,說(shuō)明改變引水時(shí)機(jī)能夠減小從A河引水量,節(jié)約系統(tǒng)供水成本。
(2)當(dāng)從A河引水量為7.65億m3時(shí),改進(jìn)模型的缺水量為0,該引水量下水庫(kù)1、水庫(kù)2的水位變化過(guò)程見(jiàn)圖6和圖7,可以看出,原模型水庫(kù)1僅在2月、9月發(fā)揮了調(diào)蓄作用,水庫(kù)2僅在5月、8月、9月發(fā)揮了調(diào)蓄作用,而改進(jìn)模型的水庫(kù)1和水庫(kù)2的蓄泄變化過(guò)程更加明顯,說(shuō)明改進(jìn)模型充分發(fā)揮了水庫(kù)的調(diào)蓄作用。原模型、改進(jìn)模型從A河引水過(guò)程見(jiàn)圖8,可知:7月、9—12月當(dāng)?shù)厮枯^富裕,改進(jìn)模型比原模型減小了引水量,而在1—4月、6月加大了從A河引水量,彌補(bǔ)需水缺口;5月需水量較大,而改進(jìn)模型從A河引水量小于原模型,原因是改進(jìn)模型1—4月從A河引水量較大,增加了水庫(kù)1、水庫(kù)2的蓄水量(5月初水庫(kù)1、水庫(kù)2蓄水量較大),5月從A河引水,水庫(kù)1和水庫(kù)2同時(shí)向用戶(hù)供水,故改進(jìn)模型5月從A河引水量小于原模型;8月當(dāng)?shù)厮糠浅8辉?,原模型從A河引水量較小,而改進(jìn)模型引水量較大,原因是8月當(dāng)?shù)厮饕獊?lái)自節(jié)點(diǎn)N1002、N1004、N1006、N2001、N2002、N3002、N4001、N5002、N5003,除N2001(水 庫(kù)1)、N4001(水庫(kù)2)外,其他節(jié)點(diǎn)當(dāng)?shù)厮浚?.35億m3,占8月當(dāng)?shù)厮偭康?2.76%)均無(wú)法存入水庫(kù)1和水庫(kù)2,不能被調(diào)蓄,因此需要加大從A河引水量,來(lái)滿(mǎn)足9—12月用戶(hù)需求。從圖6、圖7可以看出,8月末水庫(kù)1和水庫(kù)2蓄水量明顯增加,且9—12月為開(kāi)閘供水狀態(tài),水庫(kù)水位不斷降低,直至水庫(kù)放空。綜上所述,與原模型相比,改進(jìn)模型充分發(fā)揮了當(dāng)?shù)厮畮?kù)的調(diào)蓄作用。
圖6 水庫(kù)1水位變化過(guò)程
圖7 水庫(kù)2水位變化過(guò)程
圖8 原模型、改進(jìn)模型從A河引水過(guò)程
本文以北方某城市水資源系統(tǒng)為例,建立中長(zhǎng)期水量調(diào)度模型,采用線性求解器對(duì)模型進(jìn)行求解,并與遺傳算法進(jìn)行了對(duì)比。主要結(jié)論如下:①線性求解器Cplex、Lpsolver編程簡(jiǎn)單,能夠獲得全局最優(yōu)解,求解精度高、速度快;②當(dāng)模型規(guī)模增大時(shí),遺傳算法計(jì)算效率明顯下降,而Cplex、Lpsolver計(jì)算效率變化不明顯;③通過(guò)改變引水時(shí)機(jī)能夠充分發(fā)揮現(xiàn)有水庫(kù)的調(diào)蓄能力,增加系統(tǒng)總供水量。
線性求解器應(yīng)用的前提是模型必須為線性模型,因此如何實(shí)現(xiàn)含非線性約束的中長(zhǎng)期水量調(diào)度模型的線性化處理,是下一步的研究重點(diǎn)。