雍龍泉

(陜西理工大學(xué) 數(shù)學(xué)與計算機(jī)科學(xué)學(xué)院, 陜西 漢中 723000)

考慮非線性方程組

(1)

記向量x=(x1,x2,…,xn)，向量函數(shù)F(x)=(f1(x),f2(x),…,fn(x))T，則式(1)等價于非線性方程

F(x)=0。

(2)

很多物理、化工、工業(yè)里復(fù)雜的問題都可以轉(zhuǎn)化為非線性方程組，求解非線性方程組最基本的方法是牛頓法以及各類修正牛頓法等[1-8]。多步迭代的高階牛頓迭代法(諸如5階、7階、8階、9階牛頓迭代法等[9-19])成為當(dāng)前研究的熱點。當(dāng)然，收斂階越高，計算代價也就越大。

實際計算過程中，若其雅克比矩陣奇異，即行列式det(F′(x)) =0，可以采用阻尼牛頓法或延拓法進(jìn)行處理。目前多數(shù)牛頓迭代法都依賴于初始點的選取和F(x)的性態(tài)；此外針對存在多個根的非線性方程組，牛頓迭代法無法同時找到多個根。本文提出了一種求非線性方程組的新方法，具體做法是把非線性方程組轉(zhuǎn)化為一個無約束優(yōu)化，采用正弦余弦算法求解與之等價的無約束優(yōu)化問題而獲得原問題盡可能多的根。該方法的優(yōu)點是無需計算F(x)的雅克比矩陣，且能找到多個根。

1 建立優(yōu)化模型

為了求解非線性方程組，建立如下無約束優(yōu)化模型

(3)

(3)式是一個標(biāo)準(zhǔn)的無約束優(yōu)化問題，為了避免計算F(x)的雅克比矩陣，下面采用正弦余弦算法求解無約束優(yōu)化(3)式。

2 正弦余弦算法

澳大利亞學(xué)者M(jìn)irjalili[20]在2016年提出的一種采用正弦函數(shù)和余弦函數(shù)對種群進(jìn)行更新的智能優(yōu)化算法，之后便形成了正弦余弦算法(Sine Cosine Algorithm，SCA)。

SCA 算法結(jié)構(gòu)簡單，容易實現(xiàn)，其最顯著的特點是基于正弦函數(shù)和余弦函數(shù)值的變化來達(dá)到尋優(yōu)目的。在SCA算法中，主要有r1、r2、r3、r44個參數(shù)。其中，最關(guān)鍵的參數(shù)是r1，控制算法從全局搜索到局部開發(fā)的轉(zhuǎn)換。當(dāng)r1的值較大時，算法傾向于全局探索；當(dāng)r1的值較小時，算法偏向于局部開發(fā)。

下面給出SCA算法的主要步驟。

初始化算法參數(shù): 種群規(guī)模N，空間維數(shù)D，控制參數(shù)a，最大迭代次數(shù)Tmax；

在可行域空間中隨機(jī)初始化N個個體組成初始種群；t=1；

計算當(dāng)前每個個體的適應(yīng)值，并記錄最優(yōu)個體位置Pit；

while (t

fori=1 toNdo

forj=1 toDdo

根據(jù)式r1=a-at/Tmax計算r1的值；

隨機(jī)產(chǎn)生r2∈U[0,2π],r3∈U[0,2],r4∈U[0,1]；

ifr4<0.5

else

end if

end for

end for

越界處理；

計算每個個體的適應(yīng)值并更新整個種群的最優(yōu)值；t=t+1；

end while

圖1分別給出了a=2、Tmax=100與a=2、Tmax=500時r1sinr2與r1cosr2的圖像。且只有當(dāng)r1>1時，函數(shù)r1sinr2與r1cosr2值才有可能大于1或者小于-1；當(dāng)r1≤1時，函數(shù)r1sinr2與r1cosr2值必在-1和1之間。

各級人影部門均按需設(shè)崗，市級人影業(yè)務(wù)人員分為：作業(yè)天氣過程預(yù)報崗、作業(yè)條件潛力預(yù)報崗、作業(yè)條件監(jiān)測崗、聯(lián)合作業(yè)方案制定崗、作業(yè)跟蹤指揮崗、登機(jī)作業(yè)指揮崗、登機(jī)技術(shù)保障崗、外場地面保障崗、作業(yè)信息管理崗、作業(yè)效益評估崗、云水資源評估崗、指揮平臺運(yùn)行保障崗、探測設(shè)備保障崗、網(wǎng)絡(luò)運(yùn)行管理崗、信息上報及效果評估崗；旗縣級人影業(yè)務(wù)人員分為：綜合管理崗、火箭作業(yè)崗、應(yīng)急保障崗、高炮管理崗、信息上報及效果評估崗等；作業(yè)點人員分為現(xiàn)場指揮崗、作業(yè)實施崗、空域請示及信息上報崗。這樣調(diào)整以后，做到了每個崗位均由專人負(fù)責(zé)，專職專崗、崗崗有人、崗崗有責(zé)。

根據(jù)SCA算法設(shè)計原理，算法先進(jìn)行全局搜索再進(jìn)行局部開發(fā)。如圖2所示,當(dāng)|r1sinr2|>1或者|r1cosr2|>1時，算法進(jìn)行全局搜索；當(dāng)|r1sinr2|<1或者|r1cosr2|<1時，算法進(jìn)行局部開發(fā)。

劉勇等[21]對其參數(shù)r1展開研究，提出轉(zhuǎn)換參數(shù)r1非線性遞減的正弦余弦算法;曲良東等[22]提出了正弦算法。目前SCA算法已成功應(yīng)用到工農(nóng)業(yè)中，并獲得了較好的結(jié)果。更多的SCA算法見文獻(xiàn)[23-26]。下面應(yīng)用SCA算法求解(3)式。

3 數(shù)值算例

給出幾個非線性方程組，算例1—3的共同特點是其雅克比矩陣存在奇異點，故牛頓法難以進(jìn)行。通過將其轉(zhuǎn)化為無約束優(yōu)化，采用SCA算法求解，程序用MATLABR2009a編寫，參數(shù)設(shè)置N=30,a=2,Tmax=1000，算例1—3搜索空間為-100≤xi≤100,i=1,2；算例4搜索空間為-0.5≤xi≤1.5,i=1,2,…,n。為消除隨機(jī)數(shù)對算法的影響，算法運(yùn)行10次。

(a) a=2,Tmax=100 (b) a=2,Tmax=500圖1 r1sinr2與r1cosr2的圖像

圖2 SCA算法搜索原理(r1=2,r3=1)

算例1 考慮非線性方程組

(a) 適應(yīng)值收斂曲線 (b) 平均適應(yīng)值收斂曲線圖3 算例1運(yùn)行10次的適應(yīng)值和平均適應(yīng)值收斂曲線

算例2 考慮非線性方程組

其根為x*=(3,0.5)T，x*=(81/32，-1/3)T。該方程的雅克比矩陣存在奇異點(說明：奇異點分布在x1=0,x2=1,x2=0,x2=-2 四條直線上)，因此牛頓法難以進(jìn)行。圖4給出了SCA算法運(yùn)行10次的適應(yīng)值收斂曲線和平均適應(yīng)值收斂曲線圖；表1中給出了SCA算法運(yùn)行10次后得到的結(jié)果。

(a) 適應(yīng)值收斂曲線 (b) 平均適應(yīng)值收斂曲線圖4 算例2運(yùn)行10次的適應(yīng)值和平均適應(yīng)值收斂曲線

算例3 考慮非線性方程組

(a) 適應(yīng)值收斂曲線 (b) 平均適應(yīng)值收斂曲線圖5 算例3運(yùn)行10次的適應(yīng)值和平均適應(yīng)值收斂曲線

算例4 考慮非線性方程組

其根為x*=(1,1,…,1)T。圖6分別給出了n=5時SCA算法運(yùn)行10次的適應(yīng)值收斂曲線和平均適應(yīng)值收斂曲線圖。

(a) 適應(yīng)值收斂曲線 (b) 平均適應(yīng)值收斂曲線圖6 算例4運(yùn)行10次的適應(yīng)值和平均適應(yīng)值收斂曲線

算例5 考慮非線性方程組

其解為x*=(1,1)T,x*=(0,0)T。該方程非光滑，因此牛頓法難以進(jìn)行。圖7分別給出了SCA算法運(yùn)行10次的適應(yīng)值收斂曲線和平均適應(yīng)值收斂曲線圖；表1中給出了SCA算法運(yùn)行10次后得到的結(jié)果。

(a) 適應(yīng)值收斂曲線 (b) 平均適應(yīng)值收斂曲線圖7 算例5運(yùn)行10次的適應(yīng)值和平均適應(yīng)值收斂曲線

表1 算例1、2、3、5計算10次的結(jié)果

4 結(jié)束語

通過把非線性方程組轉(zhuǎn)化為無約束優(yōu)化，采用SCA算法求解，針對唯一根的非線性方程組，該方法能夠收斂到其唯一根；針對具有多個根的非線性方程組，該方法能夠找到盡可能多的根；且該算法無需計算雅克比矩陣，因此對雅克比矩陣是否存在奇異點無限制，適用于多數(shù)非線性方程組。計算結(jié)果表明, 該方法計算效率高, 對非線性方程組求解較為有效。如果要得到較高精度的數(shù)值解，可以通過增加迭代次數(shù)實現(xiàn)(如需源代碼，請與作者聯(lián)系)。