基于遺傳算法優(yōu)化粒子群算法的支斗兩級(jí)渠系優(yōu)化配水研究

高 建，張運(yùn)鑫，2

（1.河北工程大學(xué)水利水電學(xué)院，河北 邯鄲 056038；2.河北省水資源高效利用工程技術(shù)研究中心，河北 邯鄲 056038）

0 引 言

中國是世界上13 個(gè)貧水國家之一。2021 年我國用水總量為5 920.2 億m3，其中農(nóng)業(yè)用水量3 644.3 億m3，占總用水量的61.5%[1]。水資源短缺問題日益成為制約我國農(nóng)業(yè)和經(jīng)濟(jì)社會(huì)發(fā)展的重要因素。因此，優(yōu)化灌區(qū)輸配水渠系運(yùn)行管理對(duì)提升灌區(qū)渠道輸配水效率，改善灌區(qū)灌溉供用水狀況具有重要意義。

明渠是被廣泛應(yīng)用于灌區(qū)的輸配水方式[2]。灌區(qū)輸配水渠系一般由干、支、斗、農(nóng)、毛五級(jí)渠道構(gòu)成，其中支、斗兩級(jí)渠道是灌區(qū)渠系配水由續(xù)灌轉(zhuǎn)變?yōu)檩喒嗟年P(guān)鍵銜接部分，對(duì)實(shí)現(xiàn)灌區(qū)渠系優(yōu)化配水和提高渠系水利用系數(shù)方面具有重要作用。

有學(xué)者對(duì)支、斗渠合理配水調(diào)度進(jìn)行了研究，Suryavanshi和Reddy提出將輪灌分組問題轉(zhuǎn)化為數(shù)學(xué)問題，引入了流管概念，將0-1規(guī)劃模型中的二元決策變量Xij看作第i個(gè)流管的第j個(gè)分水口，它的值為0或1，0表示該分水口關(guān)閉，1表示該分水口打開[3]，通過這種數(shù)學(xué)規(guī)劃的方式獲得最優(yōu)配水調(diào)度方案。呂宏興等改進(jìn)了Suryavanshi 的模型，通過進(jìn)行引水時(shí)間的均一化處理，解決了渠道配水時(shí)需要多次調(diào)節(jié)進(jìn)水閘的問題[4]。

目前，大多數(shù)灌區(qū)仍憑借以往管理經(jīng)驗(yàn)來確定灌溉渠系的水量和流量分配，即采用經(jīng)驗(yàn)配水法。盡管采用經(jīng)驗(yàn)配水方法可以進(jìn)行灌區(qū)配水，但經(jīng)常出現(xiàn)渠系輸配水時(shí)間長、小流量配水以及水資源浪費(fèi)等問題。經(jīng)驗(yàn)配水方法無法實(shí)現(xiàn)精準(zhǔn)配水，也難以對(duì)供需變化做出及時(shí)響應(yīng)[5]。隨著水資源供需矛盾日益加劇，對(duì)提升灌區(qū)供需水管理水平的要求也越來越高[6]。為了實(shí)現(xiàn)渠系的高效優(yōu)化配水，學(xué)者們開始研究將智能優(yōu)化算法引入到渠系配水中，如遺傳算法[6-12]、粒子群算法[12-15]、模擬退火算法[16]、灰狼優(yōu)化算法[17]等，通過對(duì)這些算法應(yīng)用獲得渠系最優(yōu)調(diào)度方案。馬孝義等建立了下級(jí)渠道引水流量不相等情況下的渠道優(yōu)化配水模型，通過遺傳算法獲得優(yōu)化配水方案，輸水損失更低[9]；劉照等以渠道輸水損失最小與輪灌組間引水持續(xù)時(shí)間差異值最小為目標(biāo)建立多目標(biāo)優(yōu)化配水模型，通過雙層粒子群算法求解優(yōu)化方案，不但可減少渠道棄水量，同時(shí)統(tǒng)一了輪灌組閘門關(guān)閉時(shí)間，減少了閘門啟閉次數(shù)[13]；劉葉等以渠系輸水損失最小為目標(biāo)構(gòu)建兩級(jí)灌溉渠系優(yōu)化配水模型，通過遺傳算法確定的方案，不但能滿足灌溉要求，同時(shí)減少渠系配水時(shí)間以及閘門啟閉次數(shù)[10]；徐淑琴等以灌區(qū)渠道輸水損失最小及干渠水位變化均勻?yàn)閮?yōu)化目標(biāo)，運(yùn)用帶精英策略非支配排序遺傳算法確定運(yùn)行方案，在滿足優(yōu)化目標(biāo)的前提下，可以縮短各生育期的配水時(shí)間，且流量變化幅度較小[11]。

針對(duì)優(yōu)化離散空間約束問題，Kennedy 和Eberhart 在傳統(tǒng)粒子群算法(PSO)基礎(chǔ)上改進(jìn)并提出離散二進(jìn)制粒子群算法(BPSO)[18,19]。本文在已有研究基礎(chǔ)上，綜合考慮離散二進(jìn)制粒子群算法(BPSO)和遺傳算法(GA)的優(yōu)缺點(diǎn)[20]，研究提出了混合二進(jìn)制粒子群算法（GA-BPSO），用于灌區(qū)支、斗兩級(jí)渠系優(yōu)化配水研究。用二元決策變量Xij描述每一個(gè)支渠向斗渠配水的分水口配水運(yùn)行狀態(tài)，其中賦值0表示該分水口關(guān)閉；賦值為1 則表示該分水口開啟。應(yīng)用MATLAB 對(duì)BPSO 算法和GABPSO算法進(jìn)行模擬計(jì)算，并通過應(yīng)用案例進(jìn)行檢驗(yàn)分析。

1 支斗渠系優(yōu)化配水模型

當(dāng)支渠向斗渠配水時(shí)，往往將斗渠分為幾個(gè)組，按照“組間續(xù)灌，組內(nèi)輪灌”的模式運(yùn)行，支渠向斗渠配水灌溉如圖1所示。

圖1 支渠向斗渠配水示意圖Fig.1 Irrigation schematic from branch canal to lateral canal

設(shè)支渠的凈流量為Q凈，各斗渠的流量為q，斗渠渠道水利用系數(shù)為η，并用向上取整公式ceil 和向下取整公式floor 進(jìn)行取整，灌溉編組數(shù)M計(jì)算公式如下：

本研究提出的模型以支渠向斗渠輸配水運(yùn)行時(shí)間最短作為目標(biāo)，模型的決策變量主要包括灌溉組數(shù)、分水口數(shù)和每個(gè)分水口的連續(xù)引水時(shí)間。約束條件包括引水流量和任意分水口開啟次數(shù)等，優(yōu)化模型如下：

式中：Tb表示支渠運(yùn)行時(shí)間，s；Tl表示M 個(gè)灌溉斗渠編組中運(yùn)行時(shí)間最長的灌溉斗渠編組運(yùn)行時(shí)間，s；M 為灌溉編組個(gè)數(shù)；N 為被配水斗渠條數(shù)；ti為第i 組灌溉斗渠編組運(yùn)行時(shí)間，s，其中i=1,2,…,M；tj表示第j條斗渠配水運(yùn)行時(shí)間，s，其中j=1,2,…,N；Sj表示第j 條斗渠所需的灌水量，m3；qj表示第j 條斗渠配水流量，m3/s；ηj表示第j條斗渠渠道水利效率；Xij表示第i個(gè)灌溉編組中第j個(gè)分水口的狀態(tài)，Xij=0 為i個(gè)灌溉編組中第j個(gè)分水口是關(guān)閉狀態(tài)，即第i個(gè)灌溉編組內(nèi)的第j個(gè)分水口沒有配水計(jì)劃；Xij=1 為i 個(gè)灌溉編組中第j 個(gè)分水口是打開狀態(tài)，即第i個(gè)灌溉編組內(nèi)的第j個(gè)分水口有配水計(jì)劃。

約束條件如下：

設(shè)第i 灌溉編組供水時(shí)間為ti，則其總供水時(shí)間不應(yīng)大于灌溉周期T，即：

任意分水口在灌溉周期T內(nèi)只能開啟一次，即：

決策變量取值約束：

各灌溉斗渠編組凈流量之和不大于對(duì)應(yīng)支渠來水凈流量，即：

式中：Qi表示第i個(gè)灌溉斗渠編組的凈流量，m3/s。

2 混合二進(jìn)制粒子群算法

2.1 粒子群算法原理（PSO）

傳統(tǒng)粒子群算法（Particle Swarm Optimization，PSO）[18]是在連續(xù)區(qū)域內(nèi)對(duì)所求問題進(jìn)行優(yōu)化求解的一種算法，粒子群中的每一個(gè)粒子個(gè)體均表示一種可行的解決方案，其算法原理是通過粒子之間的簡單行為，經(jīng)過交互得出最優(yōu)解，具有較快的收斂效率和局部搜索能力但容易陷入局部最優(yōu)解的特點(diǎn)。

粒子群算法的基本求解計(jì)算表達(dá)式如下：

式中：v(i)(t)表示第i 個(gè)粒子在第t 次的移動(dòng)速度；x(i)(t)表示第i個(gè)粒子在第t次粒子所在的位置；x(i)(t+1)表示第i個(gè)粒子在第t+1次所在的位置；c1、c2為粒子的學(xué)習(xí)因子；rand1、rand2為分布在[0,1]的隨機(jī)數(shù)；p(i)best為第i個(gè)粒子歷史最優(yōu)解；gbest為當(dāng)前歷史最優(yōu)解。

進(jìn)一步分析可將式(10)中第i個(gè)粒子在第t+1次的移動(dòng)的速度可分解為三部分：w v(i)(t)為第一部分，稱為“慣性”部分，這一部分反映了粒子原有的運(yùn)動(dòng)趨勢；c1rand1[p(i)best-x(i)(t)]為第二部分，稱為“個(gè)人認(rèn)知”部分，這一部分就反映了粒子向自身粒子最優(yōu)解移動(dòng)的趨勢；c2rand2[gbest- x(i)(t)]為第三部分，稱為“社會(huì)認(rèn)知”部分，這一部分反映了粒子向全部粒子中的最優(yōu)解移動(dòng)的趨勢。粒子移動(dòng)過程圖和粒子群算法基本流程圖分別如圖2和圖3所示。

圖2 PSO流程圖Fig.2 Flow chart of PSO

圖3 PSO粒子移動(dòng)過程圖Fig.3 PSO particle movement process diagram

典型粒子群算法步驟如下：

（1）按照需求對(duì)粒子進(jìn)行初始化，包括粒子速度和粒子位置；

（2）構(gòu)建適應(yīng)度函數(shù)，對(duì)每個(gè)已經(jīng)生成的粒子按照適應(yīng)度函數(shù)計(jì)算相應(yīng)的適應(yīng)度；

（3）根據(jù)已經(jīng)計(jì)算出的適應(yīng)度評(píng)價(jià)出個(gè)體最佳粒子以及全局最佳粒子；

（4）根據(jù)公式（10）和公式（11）更新每個(gè)粒子的速度和位置；

（5）重新計(jì)算更新過后的每個(gè)粒子的適應(yīng)度；

（6）重新計(jì)算過后的適應(yīng)度與個(gè)體最佳粒子和全局最佳粒子進(jìn)行比較，重新選出個(gè)體最佳粒子和全局最佳粒子；

（7）依照結(jié)束條件進(jìn)行判斷，滿足結(jié)束條件結(jié)束算法；若不滿足則重復(fù)執(zhí)行上述步驟4）、5）、6）直到滿足條件。

2.2 離散二進(jìn)制粒子群算法原理（BPSO）

Kennedy 和Eberhart 為解決離散型或二進(jìn)制類型的問題從而改進(jìn)了常規(guī)的粒子群算法，提出了離散二進(jìn)制粒子群算法（Discrete Binary Particle Swarm Optimization Algorithm，BPSO）[19]。該算法是采用二進(jìn)制編碼0或1對(duì)粒子中某一維度的狀態(tài)進(jìn)行表示，速度更新公式同基本粒子群算法的速度更新公式相同，但沒有使用粒子群算法位置更新公式，取而代之使用sigmoid函數(shù)將速度值映射到[0,1]區(qū)間，以此表示取0或1 的概率。通過rand 隨機(jī)數(shù)與sigmoid 函數(shù)映射之后的值進(jìn)行比較，得到位置的具體值[21]。

下式分別為sigmoid函數(shù)以及位置更新公式：

2.3 遺傳算法（GA）

遺傳算法（Genetic Algorithm，GA）借鑒了達(dá)爾文的進(jìn)化論和孟德爾的遺傳學(xué)說，是通過模擬生物在自然環(huán)境中的遺傳和進(jìn)化的過程而形成的一種自適應(yīng)全局優(yōu)化搜索算法，具有全局搜索能力強(qiáng)，局部搜索能力較弱的特點(diǎn)[20]。圖4 為遺傳算法流程圖。

圖4 GA流程圖Fig.4 Flow chart of GA

其算法流程大致如下：

（1）按照編碼需求（二進(jìn)制編碼，實(shí)數(shù)編碼等）對(duì)種群進(jìn)行初始化；

（2）構(gòu)建適應(yīng)度函數(shù)，對(duì)每個(gè)已經(jīng)生成的個(gè)體按照適應(yīng)度函數(shù)計(jì)算相應(yīng)的適應(yīng)度；

（3）根據(jù)需求對(duì)已經(jīng)生成的個(gè)體進(jìn)行選擇，一般選擇適應(yīng)度較高的個(gè)體為父代；

（4）按照設(shè)定好的交叉概率及交叉規(guī)則，將選擇出來的父代進(jìn)行隨機(jī)配對(duì)，進(jìn)行交叉運(yùn)算，生成新的個(gè)體；

（5）按照設(shè)定好的變異概率及變異規(guī)則，將選擇出來的父代進(jìn)行變異操作，生成新的個(gè)體；

（6）重新計(jì)算適應(yīng)度并更新種群；

（7）依照結(jié)束條件進(jìn)行判斷，滿足結(jié)束條件結(jié)束算法；若不滿足重復(fù)執(zhí)行上述步驟（3）、（4）、（5）、（6）直到滿足條件。

2.4 混合二進(jìn)制粒子群算法（GA-BPSO）

根據(jù)GA 和BPSO 的特點(diǎn)，在BPSO 的基礎(chǔ)上融入GA 的方法，融合兩種算法的優(yōu)點(diǎn)并彌補(bǔ)相應(yīng)不足，得到一種整體性能更優(yōu)的混合二進(jìn)制粒子群算法(GA-BPSO)。圖5為GA-BPSO基本流程圖。

圖5 GA-BPSO流程圖Fig.5 Flow chart of GA-BPSO

其算法流程大致如下：

（1）按照需求對(duì)粒子進(jìn)行初始化，包括粒子速度和粒子位置；

（2）對(duì)每個(gè)已經(jīng)生成的粒子按照適應(yīng)度函數(shù)計(jì)算相應(yīng)的適應(yīng)度；

（3）根據(jù)已經(jīng)計(jì)算出的適應(yīng)度評(píng)價(jià)出個(gè)體最佳粒子以及全局最佳粒子；

（4）根據(jù)公式（10）更粒子新速度，根據(jù)公式（12）、（13）更新粒子的位置；

（5）計(jì)算更新后的粒子適應(yīng)度；

（6）重新選出個(gè)體最佳粒子和全局最佳粒子；

（7）按照需求選擇出一定數(shù)量的個(gè)體作為父代，并進(jìn)行交叉的相關(guān)操作；

（8）按照需求選擇出一定數(shù)量的個(gè)體作為父代，并進(jìn)行變異的相關(guān)操作；

（9）對(duì)新生成粒子進(jìn)行適應(yīng)度計(jì)算，根據(jù)結(jié)果更新粒子并重新選出個(gè)體最佳粒子和全局最佳粒子；

（10）依照結(jié)束條件進(jìn)行判斷，滿足結(jié)束條件結(jié)束算法；若不滿不足重復(fù)執(zhí)行步驟（4）、（5）、（6）、（7）、（8）、（9）直到滿足條件。

3 實(shí)例應(yīng)用

某灌區(qū)一條支渠下有7 條斗渠，支渠設(shè)計(jì)流量2.5 m3/s，渠道水利用系數(shù)為0.92，各斗渠設(shè)計(jì)流量、控制面積以及渠道水利用系數(shù)如表1 所示。支渠灌溉周期為7 d，本次灌水定額確定為900 m3/hm2，優(yōu)化配水目標(biāo)為渠系輸配水時(shí)間最短。

表1 各斗渠參數(shù)Tab.1 Parameters of each lateral canal

3.1 案例分析

根據(jù)該灌區(qū)支渠、斗渠相關(guān)信息以及公式(1)確定灌溉編組數(shù)M=4 最佳。根據(jù)相關(guān)約束可知每條斗渠只能在一個(gè)灌溉編組中運(yùn)行，結(jié)合算法可將每一條斗渠看作一個(gè)1×4 的列向量，均由0 或1 組成。其中0 代表該斗渠在某灌溉編組內(nèi)的配水口為關(guān)閉狀態(tài)，即該灌溉編組內(nèi)該支渠沒有配水計(jì)劃；1代表該斗渠在某灌溉編組內(nèi)的配水口為打開狀態(tài)，即該灌溉編組內(nèi)該斗渠有配水計(jì)劃。如[1,0,0,0]T表示該斗渠在第一個(gè)灌溉組中運(yùn)行，而不在其余灌溉組中運(yùn)行，示意如圖6(a)所示。

圖6 算法分析圖Fig.6 Algorithm analysis graph

按照上述設(shè)想，如果將灌溉編組和渠道成生一個(gè)只包含0或1 的二維數(shù)組，初始化各灌溉斗渠編組分水口狀態(tài)如表2 所示，則當(dāng)前配水狀態(tài)為7條斗渠均在第一個(gè)灌溉編組中，其余灌溉編組無斗渠配水計(jì)劃，按此設(shè)定方式會(huì)降低算法的運(yùn)行效率。所以對(duì)BPSO 的粒子初始化做出改變，在生成粒子時(shí)，生成一個(gè)7 維的向量，每一個(gè)維度可看成每一條斗渠變化狀態(tài)，其中：假設(shè)這個(gè)粒子初始化時(shí)為[0,0,0,0,0,0,0]，當(dāng)粒子某一維度發(fā)生變化時(shí)則對(duì)應(yīng)斗渠分水口在灌溉編組中工作狀態(tài)就會(huì)發(fā)生變，從第i個(gè)灌溉編組變換到第i+1個(gè)，具體如圖6(b)所示。當(dāng)粒子由[0,0,1,0,0,0,0]變?yōu)閇0,0,0,0,0,0,0]則第三條斗渠會(huì)從第三個(gè)灌溉編中有配水計(jì)劃變?yōu)樵诘谒膫€(gè)灌溉編組中有配水計(jì)劃。當(dāng)某斗渠在最后一個(gè)灌溉編組中有配水計(jì)劃時(shí)下一次變換就會(huì)變換到在第一個(gè)灌溉編組中有配水計(jì)劃；如果粒子某一維度未發(fā)生改變則對(duì)應(yīng)的斗渠分水口在灌溉編組中的狀態(tài)就保持不變。

表2 初始化各灌溉斗渠編組分水口狀態(tài)Tab.2 Initialize the status of each irrigation lateral grouping manifold

3.2 算法應(yīng)用

本文在PC 端使用MATLAB 軟件求解算法。將渠系配水時(shí)間作為適應(yīng)度函數(shù)值，最佳配水方案將具有最小的適應(yīng)度函數(shù)值。算法參數(shù)確定：種群Z=20，迭代次數(shù)N=100，慣性權(quán)重W=0.9，個(gè)體學(xué)習(xí)因子C1=2，社會(huì)學(xué)習(xí)因子C2=2，交叉概率Pc=0.7，變異概率Pm=0.05。按照公式（1）確定灌溉斗渠分組數(shù)M=4，根據(jù)各灌溉斗渠編組分水口狀態(tài)，按照公式（4）計(jì)算出每個(gè)灌溉斗渠分組運(yùn)行的時(shí)間ti，再根據(jù)公式（3）得出所有ti運(yùn)行時(shí)間最長的灌溉斗渠分組運(yùn)行的時(shí)間Tl，將Tl作為算法的適應(yīng)度函數(shù)值，使其盡可能小以達(dá)到算法求解要求。

3.3 結(jié)果分析

針對(duì)上述兩級(jí)渠系輸配水問題案例，可得到如下分析：

（1）兩種算法得到的配水編組的方式都可以將渠系配水時(shí)間降到最低，具體方案如表3和表4所示。

表3 BPSO和GA-BPSO計(jì)算最優(yōu)配水方案結(jié)果（1）Tab.3 Results of BPSO and GA-BPSO calculations of optimal water distribution schemes (1)

表4 BPSO和GA-BPSO計(jì)算最優(yōu)配水方案結(jié)果（2）Tab.4 Results of BPSO and GA-BPSO calculations of optimal water distribution schemes (2)

（2）通過表5 以及圖7 可知，在參數(shù)相同的情況下BPSO在大概21 代左右能得到最優(yōu)解，而GA-BPSO 則在大概12 代左右得到最優(yōu)解，且兩者最優(yōu)解相同?？梢?，針對(duì)BPSO 的不足之處結(jié)合GA 之后的新算法GA-BPSO 的整體性能優(yōu)于BPSO。

表5 BPSO和GA-BPSO迭代終止時(shí)的性能對(duì)比Tab.5 Performance comparison between BPSO and GA-BPSO at iterative termination

圖7 算法結(jié)果對(duì)比Fig.7 Comparison of algorithm results

（3）由于在BPSO 中包含速度和位置兩個(gè)變量，在與GA結(jié)合時(shí)交叉、變異操作會(huì)導(dǎo)致速度、位置均發(fā)生相應(yīng)改變，而速度變化會(huì)影響下一次迭代時(shí)位置發(fā)生改變的概率，所以在模擬過程中會(huì)出現(xiàn)個(gè)別GA-BPSO效果不如BPSO的情況。

4 結(jié) 論

本文針對(duì)灌區(qū)支斗兩級(jí)渠系輸配水優(yōu)化問題，根據(jù)BPSO和GA的優(yōu)缺點(diǎn)，在BPSO基礎(chǔ)上融合GA的思想，構(gòu)建一種更為高效的GA-BPSO。通過GA-BPSO 對(duì)案例問題進(jìn)行優(yōu)化，模擬結(jié)果表明GA-BPSO 達(dá)到最優(yōu)解的平均收斂代數(shù)要比BPSO達(dá)到最優(yōu)解的平均收斂代數(shù)要低，說明結(jié)合后的GA-BPSO 的整體性能相比BPSO 要更好。證實(shí)了支斗兩級(jí)渠系在“組間續(xù)灌，組內(nèi)輪灌”的輸配水模式下GA-BPSO 具有更好的可行性和高效性，可為灌區(qū)用水調(diào)度提供一種較優(yōu)的解決方案。由于遺傳操作對(duì)速度的改變可能會(huì)對(duì)下一次迭代時(shí)位置的改變產(chǎn)生影響，在未來算法優(yōu)化的方向上可以進(jìn)一步研究GABPSO 運(yùn)行過程中的速度更新方式及更高效的交叉、變異方式。