薛 丹，姚若俠

(陜西師范大學(xué) 計(jì)算機(jī)科學(xué)學(xué)院，陜西 西安 710119)

0 引 言

在科學(xué)和工程領(lǐng)域中，許多重要的問題常?？梢詺w結(jié)為對大型稀疏線性代數(shù)系統(tǒng)的求解問題。這類線性代數(shù)系統(tǒng)的系數(shù)矩陣階數(shù)較高，采用直接法進(jìn)行Gauss消去或者是矩陣的三角分解，計(jì)算量大并且復(fù)雜度高，對于這種情況，通常采用迭代法求解[1]。迭代法是指從一個(gè)初始向量出發(fā)，按照給定的迭代公式，逐次迭代來逼近方程組的解，因此，迭代法是否收斂，收斂速度以及迭代次數(shù)成為衡量迭代算法非常重要的指標(biāo)。逐次超松弛迭代法(successive overrelaxation，SOR)是對Gauss-Seidel(G-S)迭代法的改進(jìn)，由Frankel和Young提出，隨后被廣泛應(yīng)用于求解核反應(yīng)擴(kuò)散、石油天然氣存儲(chǔ)、天氣預(yù)測等大型實(shí)際問題的數(shù)學(xué)模型[2]。求解方程組Ax=b的SOR迭代公式如下[1]：

(1)

其中，k=1,2,…表示每一步迭代計(jì)算；i=1,2,…n，xi表示解向量x的各個(gè)分量。松弛因子w是影響SOR迭代法的關(guān)鍵，w選取恰當(dāng)，可以使SOR迭代法快速收斂，而不恰當(dāng)?shù)膚可能減緩SOR迭代法的收斂速度甚至導(dǎo)致發(fā)散。因此，研究如何選取最優(yōu)的松弛因子w在數(shù)值代數(shù)領(lǐng)域具有重要意義。

定理1[1]：若SOR迭代收斂，則松弛因子w需滿足w(0,2)。

文獻(xiàn)[3]指出，通過數(shù)學(xué)公式確定最優(yōu)松弛因子w具有一定難度，因而另辟蹊徑提出通過計(jì)算機(jī)算法選取最優(yōu)松弛因子，即在區(qū)間(0,2)上分別采用二分比較法、黃金分割法、逐步搜索法選取最優(yōu)松弛因子w。通過對這三種算法的研究和分析，發(fā)現(xiàn)傳統(tǒng)的分割算法等雖然容易實(shí)現(xiàn)，但卻對參數(shù)設(shè)置有較強(qiáng)的依賴，且無法保證在(0,2)區(qū)間上可以成功選取全局最優(yōu)松弛因子。因此，該文根據(jù)群智能優(yōu)化算法可以在解空間根據(jù)迭代規(guī)則智能靠近最優(yōu)解的思路，擬通過粒子群優(yōu)化算法(particle swarm optimization，PSO)來實(shí)現(xiàn)在(0,2)區(qū)間上選取SOR迭代法最優(yōu)松弛因子w，以進(jìn)一步推動(dòng)通過計(jì)算機(jī)算法選取SOR最優(yōu)松弛因子w的研究發(fā)展。

1 粒子群優(yōu)化算法

粒子群優(yōu)化算法是由Eberhart和Kennedy于1995年提出的一類基于群智能的隨機(jī)優(yōu)化算法[4]。PSO算法的基本思想來源于對鳥群捕食行為的研究：一群鳥在隨機(jī)搜尋食物，如果這個(gè)區(qū)域里只有一塊食物，那么找到食物的最簡單有效的策略是，搜尋目前離食物最近的鳥的周圍區(qū)域。粒子群優(yōu)化算法自提出以來，由于形式簡單，實(shí)現(xiàn)容易，需要調(diào)整的參數(shù)較少，已被廣泛應(yīng)用于多個(gè)學(xué)科和工程領(lǐng)域[5-6]。

1.1 PSO算法

PSO算法首先初始化一組種群數(shù)為n的隨機(jī)粒子，每個(gè)隨機(jī)粒子代表D維解空間的一個(gè)解，粒子i在d維的位置用xid表示，粒子i在d維的速度用vid表示，粒子通過改變vid來改變自己的位置，不斷靠近潛在最優(yōu)解。vid的更改受到兩個(gè)極值的影響：一個(gè)極值是粒子本身找到的最優(yōu)解，稱為個(gè)體極值，記為pbestid；另一個(gè)極值是整個(gè)種群找到的最優(yōu)解，稱為全局極值，記為gbestd[4]。vid體現(xiàn)了粒子在群體中追隨潛在最優(yōu)解流動(dòng)的速度和方向，文獻(xiàn)[4]最早提出的粒子群追蹤潛在最優(yōu)解的進(jìn)化方程如下式所示[4]：

c2r2(gbestd-xid)

(2)

(3)

其中，i=1,2,…,n；d=1,2,…,D；r1和r2是[0,1]之間的隨機(jī)數(shù)；學(xué)習(xí)因子c1、c2為非負(fù)常數(shù)，通常取c1=c2=2；vid∈[-vmax,vmax]，vmax是用戶設(shè)定的常數(shù)，表示速度更改的最大值。進(jìn)化終止條件為預(yù)設(shè)的最大迭代次數(shù)或預(yù)定的最小適應(yīng)度閾值，該算法也被稱為全局版PSO算法。

1.2 PSO算法的發(fā)展

全局版PSO算法雖然形式簡單、易于理解、容易實(shí)現(xiàn)且能快速收斂至潛在最優(yōu)解，但存在容易陷入局部極值、進(jìn)化后期收斂速度慢等問題，許多學(xué)者也一直致力于研究如何提高PSO算法的尋優(yōu)能力。1998年，Yuhui Shi和Russell Eberhart通過給式(3)添加慣性權(quán)重(weight)來優(yōu)化PSO算法的能力，并通過實(shí)驗(yàn)得出結(jié)論：當(dāng)慣性權(quán)重在區(qū)間[0.9,1.2]時(shí)，PSO算法有較高的概率可以搜索到全局最優(yōu)解[7]，修改后的進(jìn)化方程如式(4)和式(5)[7]式所示：

c2r2(gbestd-xid)

(4)

(5)

很多學(xué)者將式(4)和式(5)稱為基本粒子群優(yōu)化算法(basic particle swarm optimization，bPSO)。

2006年，胡旺等提出簡化粒子群優(yōu)化算法(simple particle swarm optimization，sPSO)[8]，sPSO算法僅由粒子位置控制進(jìn)化過程。實(shí)驗(yàn)證明，該算法能夠有效避免粒子群進(jìn)化后期收斂速度變慢和精度低等問題[8]。sPSO算法的進(jìn)化方程為[8]：

(6)

同時(shí)，為了克服bPSO、sPSO算法在進(jìn)化后期容易陷入局部極值的缺點(diǎn)，提出了帶極值擾動(dòng)的粒子群優(yōu)化算法(extremum disturbed particle swarm optimization，tPSO)[8]和帶極值擾動(dòng)的簡化粒子群優(yōu)化算法(extremum disturbed and simple particle swarm optimization，tsPSO)[8]。tPSO、tsPSO算法的基本思路是：在粒子群進(jìn)化方程中，通過擾動(dòng)因子r3、r4分別對個(gè)體極值和全局極值添加隨機(jī)擾動(dòng)，以提升算法跳出局部極值的能力。具體擾動(dòng)策略為：將粒子群的進(jìn)化停滯步數(shù)作為隨機(jī)擾動(dòng)的觸發(fā)條件，用t0和tg分別表示個(gè)體極值和全局極值進(jìn)化停滯步數(shù)，用T0和Tg分別表示個(gè)體極值和全局極值需要擾動(dòng)的停滯步數(shù)閾值。T0和Tg的具體取值由用戶設(shè)定，當(dāng)t0≤T0時(shí)，個(gè)體極值不需要擾動(dòng)，r3=1；否則對個(gè)體極值進(jìn)行隨機(jī)擾動(dòng)，r3=U(0,1)。同理，當(dāng)tg≤Tg時(shí)，全局極值不需要擾動(dòng)，r4=1；否則對全局極值進(jìn)行隨機(jī)擾動(dòng)，r4=U(0,1)。tPSO算法的進(jìn)化方程為：

c2r2(r4*gbestd-xid)

(7)

(8)

tsPSO算法的進(jìn)化方程為：

(9)

擾動(dòng)因子r3、r4的取值分別為：

(10)

(11)

通過對粒子群優(yōu)化算法的研究及應(yīng)用現(xiàn)狀的調(diào)查和分析[9-14]，選取bPSO、sPSO、tPSO、tsPSO等四個(gè)算法分別應(yīng)用在(0,2)區(qū)間上針對不同方程組Ax=b選取SOR迭代法的最優(yōu)松弛因子w，實(shí)驗(yàn)具體設(shè)計(jì)、部署和結(jié)果分析將在下文詳細(xì)介紹。

2 實(shí)驗(yàn)與分析

2.1 實(shí)驗(yàn)環(huán)境及數(shù)據(jù)準(zhǔn)備

實(shí)驗(yàn)環(huán)境信息如下：Win 10操作系統(tǒng)，Intel?Core(TM)i7-8550U CPU處理器，8 GB內(nèi)存，Python開發(fā)語言，Python3.5解釋器版本，PyCharm2018.1.4開發(fā)環(huán)境。

算法的參數(shù)設(shè)置如下：種群粒子數(shù)=3，進(jìn)化世代數(shù)閾值=50，c1=c2=2，慣性系數(shù)：weight=0.9，bPSO、tPSO算法中速度更改的最大值vmax=0.3；sPSO、tsPSO算法中粒子位置更改的最小值=0.001，最大值=1.999；tPSO、tsPSO算法中個(gè)體極值需要擾動(dòng)的停滯步數(shù)閾值To=3，全局極值需要擾動(dòng)的停滯步數(shù)閾值Tg=5。實(shí)驗(yàn)中用到的兩個(gè)形為Ax=b的線性方程組數(shù)據(jù)[1, 3,16]如下：

b[1]=[-2,-6, 6,12]T

S[1]=[ 1,-2,-1, 3]T

b[2]=[0,5,0,6,-2,6]T

S[2]=[1,2,1,2,1,2]T

其中，S[i],i=1,2是方程組的精確解。

2.2 算法設(shè)計(jì)

應(yīng)用bPSO、sPSO、tPSO、tsPSO算法在(0,2)區(qū)間上對上述兩個(gè)線性方程組分別求最優(yōu)松弛因子w，算法的程序流程如圖1～圖4所示。

圖1 bPSO算法程序流程

圖2 sPSO算法程序流程

圖3 tPSO算法程序流程

圖4 tsPSO算法程序流程

2.3 實(shí)驗(yàn)結(jié)果及分析

運(yùn)行以上四個(gè)算法，分別對兩個(gè)線性方程組在(0,2)區(qū)間各執(zhí)行一次搜索SOR迭代法最優(yōu)松弛因子w，程序執(zhí)行結(jié)果見表1。由表1數(shù)據(jù)可知：

(2)相對于文獻(xiàn)[3]中應(yīng)用二分查找法、黃金分割法等搜索方程組1的SOR最優(yōu)松弛因子w，該文應(yīng)用bPSO、sPSO、tPSO、tsPSO算法搜索的w對應(yīng)的SOR迭代次數(shù)為8次，優(yōu)于文獻(xiàn)[3]中的12次。

表1 四種算法選取最優(yōu)w的實(shí)驗(yàn)結(jié)果

該文在實(shí)驗(yàn)編碼過程中調(diào)用Python的matplotlib庫，對粒子的執(zhí)行過程進(jìn)行了可視化展示。僅列出以上四個(gè)算法在搜索方程組2的SOR最優(yōu)松弛因子時(shí)的執(zhí)行過程，見圖5。

(a)bPSO算法

(b)sPSO算法

(c)tPSO算法

(d)tsPSO算法

由圖5(a)可以看出，bPSO算法搜索方程組2最優(yōu)松弛因子w時(shí)，主要搜索區(qū)間集中在(0.5,1.5)，在進(jìn)化第8次時(shí)找到最優(yōu)松弛因子1.122，對應(yīng)的SOR迭代次數(shù)為8次。圖5(b)反映了sPSO算法實(shí)際搜索區(qū)間在(0,2)，覆蓋的解區(qū)間的范圍更廣，粒子在進(jìn)化第4次時(shí)找到最優(yōu)松弛因子1.122，對應(yīng)的SOR迭代次數(shù)為8次。圖5(c)反映了tPSO算法主要搜索區(qū)間在(0,1.5)，在進(jìn)化第26次時(shí)找到最優(yōu)松弛因子1.119，對應(yīng)的SOR迭代次數(shù)為8次。圖5(d)反映了tsPSO算法主要搜索區(qū)間在(0,2)，幾乎覆蓋全部解區(qū)間，其結(jié)果最具全局最優(yōu)性，且粒子在進(jìn)化第35次時(shí)找到最優(yōu)松弛因子1.122，對應(yīng)的SOR迭代次數(shù)為8次。

3 結(jié)束語

應(yīng)用bPSO、sPSO、tPSO、tsPSO等四個(gè)算法分別對兩個(gè)不同的線性方程組搜索SOR最優(yōu)松弛因子w，獲得以下結(jié)論：

(1)bPSO、sPSO、tPSO、tsPSO算法均成功實(shí)現(xiàn)在(0,2)區(qū)間上選取SOR全局最優(yōu)松弛因子w，且實(shí)驗(yàn)結(jié)果優(yōu)于文獻(xiàn)[3]中的二分比較法等算法。

(2)對于同一個(gè)方程組的應(yīng)用中，bPSO算法的實(shí)際搜索區(qū)間較小，容易陷入局部極值。sPSO算法、tPSO算法、tsPSO算法，實(shí)際搜索范圍相較bPSO算法有所擴(kuò)大，其中tsPSO算法的搜索范圍幾乎覆蓋了全部解區(qū)間，其搜索結(jié)果也更具全局最優(yōu)性。