雍龍泉,賈 偉,黎延海

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

線性互補(bǔ)問題是應(yīng)用數(shù)學(xué)、計(jì)算數(shù)學(xué)與運(yùn)籌學(xué)相互交叉的一個(gè)研究方向，在經(jīng)濟(jì)均衡理論、雙矩陣對(duì)策、接觸力學(xué)、交通流平衡等領(lǐng)域具有廣泛的應(yīng)用[1-3].線性互補(bǔ)問題就是求一個(gè)向量x∈n，使得

x≥0,Mx+q≥0,xT(Mx+q)=0,

其中：M∈n×n,q∈n.

線性互補(bǔ)問題簡(jiǎn)記為L(zhǎng)CP(M,q)或LCP(linear complementary problem).針對(duì)存在唯一解的LCP，研究的重點(diǎn)是探尋有效的求解算法.目前已有的方法有：各類內(nèi)點(diǎn)法[4-9]、矩陣分裂迭代法[10-13]、混合整數(shù)線性規(guī)劃方法[14-15]等，其中智能優(yōu)化算法、神經(jīng)網(wǎng)絡(luò)方法[16-17]被用于求解各類互補(bǔ)問題.近年來，隨機(jī)線性互補(bǔ)、對(duì)稱錐互補(bǔ)、參數(shù)線性互補(bǔ)、張量互補(bǔ)、垂直線性互補(bǔ)[18-24]等成為研究熱點(diǎn).

針對(duì)存在多個(gè)解的LCP，在眾多解中尋找稀疏解、最小一范數(shù)解已成為當(dāng)前的一個(gè)研究熱點(diǎn)[25-27].線性互補(bǔ)問題的稀疏解即求解如下零范數(shù)優(yōu)化問題

(1)

其中：‖x‖0表示向量x中非零元素個(gè)數(shù).

零范數(shù)優(yōu)化通常的做法是采用凸松弛技巧[28-33]，把稀疏解問題L0(x)轉(zhuǎn)化為如下L1(x)范數(shù)優(yōu)化

(2)

L1(x)問題是一個(gè)凸優(yōu)化，在一定條件下，L1(x)的解恰好就是L0(x)的稀疏解.

問題(2)可以采用約束優(yōu)化算法、交替方向法[29,34]等.論文采用一種仿物理現(xiàn)象(正弦波與余弦波)的優(yōu)化算法——正弦余弦算法(sine cosine algorithm,簡(jiǎn)稱SCA)求解線性互補(bǔ)問題盡可能多的解，進(jìn)而在其中選取范數(shù)最小的解.

1 解的存在性與問題的轉(zhuǎn)化

定義[1]矩陣M∈n×n，對(duì)?x∈n，都有xTMx≥0，則稱矩陣M∈n×n為半正定矩陣；對(duì)?x∈n,x≠0，都有xTMx>0，則稱M為正定矩陣.

在線性互補(bǔ)問題研究中，矩陣的正定性起著重要的作用，這里的半正定矩陣與正定矩陣無對(duì)稱性要求.當(dāng)矩陣M是半正定矩陣時(shí)，稱LCP為單調(diào)線性互補(bǔ)問題.矩陣M對(duì)稱時(shí)，特征值全大于零，則矩陣M是正定矩陣.需要注意的是，非對(duì)稱實(shí)矩陣的特征值全大于零時(shí)不一定是正定矩陣[35].下面給出線性互補(bǔ)問題解存在的2個(gè)充分條件.

定理[1](i)設(shè)矩陣M∈n×n是一正定矩陣，則對(duì)于任意q∈n，LCP有唯一解.

(ii)設(shè)矩陣M∈n×n是一半正定矩陣，對(duì)于任意q∈n，若LCP是可行的，則其必有解，且解集為凸集.

利用文獻(xiàn)[36]中的結(jié)論，任意的LCP都等價(jià)于如下絕對(duì)值方程

(M+I)x+q=|(M-I)x+q|，

于是，論文定義智能算法的適應(yīng)值函數(shù)

(3)

作為優(yōu)化目標(biāo)函數(shù)，當(dāng)且僅當(dāng)f(x)=0時(shí)，則找到LCP的解.通過多次執(zhí)行算法，就有可能找到LCP盡可能多的解，再?gòu)闹羞x取范數(shù)最小的解.所以，問題的核心還是求解(3).

2 正弦余弦算法

SCA是一種新型仿物理現(xiàn)象(正弦波與余弦波)的優(yōu)化算法[37]，具有結(jié)構(gòu)清晰、容易實(shí)現(xiàn)等優(yōu)點(diǎn).下面給出SCA的步驟:

SCA算法參數(shù)初始化: 設(shè)定種群規(guī)模N，優(yōu)化問題變量個(gè)數(shù)D，轉(zhuǎn)換參數(shù)a及Tmax；

t=1；在搜索空間中隨機(jī)生成N個(gè)個(gè)體作為初始種群，計(jì)算個(gè)體的適應(yīng)值并記錄當(dāng)前最優(yōu)個(gè)體位置P(t)；

while(t

fori= 1 toNdo

forj= 1 toDdo

r1=a-at/Tmax；r2∈U[0,2π],r3∈U[0,2],r4∈U[0,1]；

ifr4<0.5

(4)

else

(5)

end

end

end

進(jìn)行越界處理；

更新種群的最優(yōu)位置P(t)；

t=t+1；

end.

SCA算法利用正弦波與余弦波的周期震蕩實(shí)現(xiàn)搜索過程[37-38].圖1給出了a=2、Tmax=100與a=2、Tmax=200時(shí)，r1sin(r2)與r1cos(r2)的圖像.

圖1 r1sin(r2)與r1cos(r2)的圖像

當(dāng)r1>1時(shí)，r1sin(r2)與r1cos(r2)才有可能大于1或者小于-1；當(dāng)r1≤1時(shí)，r1sin(r2)與r1cos(r2)值必介于-1和1之間.因此當(dāng)|r1sin(r2)|>1或者|r1cos(r2)|>1時(shí)，算法進(jìn)行全局搜索；當(dāng)|r1sin(r2)|≤1或者|r1cos(r2)|≤1時(shí)，算法進(jìn)行局部開發(fā).如圖2所示.

圖2 r1=2、r3=1時(shí)，r1sin(r2)和r1cos(r2)的波動(dòng)性及算法搜索策略

在SCA算法中，參數(shù)r1較為關(guān)鍵，控制算法從全局搜索到局部開發(fā)的轉(zhuǎn)換.更多的SCA算法見文獻(xiàn)[39-41].下面應(yīng)用SCA算法求解(3).

3 數(shù)值算例

該節(jié)給出幾個(gè)LCP，算例(i)～(iii)的共同特點(diǎn)是矩陣M半正定，(iv)的特點(diǎn)是M非半正定,4個(gè)LCP解的集合均為非空凸集.通過將LCP轉(zhuǎn)化為問題(3)，采用SCA算法求解.設(shè)置參數(shù)N=30,a=2,Tmax=1 000.算例(i)搜索空間為0≤xi≤2，(ii)搜索空間為0≤xi≤5，(iii)，(iv)搜索空間為0≤xi≤1,i=1,2,…,n.為消除隨機(jī)數(shù)對(duì)算法的影響，算法運(yùn)行10次.

(i)考慮線性互補(bǔ)問題

其中：對(duì)稱矩陣M的特征值為λ1=0,λ2=2,λ3=2，因此M是半正定矩陣.

文獻(xiàn)[25]中指出該問題具有無窮解x*=(λ+1,0,λ)T,λ≥0.采用SCA算法求解，得到最小范數(shù)解為x*=(1,0,0)T.事實(shí)上，x*=(1,0,0)T也是其稀疏解.

(ii)考慮線性互補(bǔ)問題

其中：對(duì)稱矩陣M的特征值為λ1=0,λ2=0,λ3=0,λ3=4，因此M是半正定矩陣.

文獻(xiàn)[25]指出該問題具有無窮解x*=(λ,0,0,5-λ,0)T,0≤λ≤5.采用SCA算法求解，可以獲得2個(gè)最小范數(shù)解同時(shí)也是稀疏解為x*=(0,0,0,5,0)T,(5,0,0,0,0)T.

(iii)考慮線性互補(bǔ)問題

其中：對(duì)稱矩陣M的特征值為λ1=0,λ2=0,λ3=3，因此M是半正定矩陣.

下面用定義來判別對(duì)稱矩陣M的正定性.對(duì)于?x, 這里x=(x1,x2,x3)T，由xTMx=(x1+x1+x3)2≥0，知矩陣M是半正定矩陣.文獻(xiàn)[26]中指出該問題具有無窮解

x*=(λ1,λ2,λ3)T,λ1+λ2+λ3=1,λ1≥0,λ2≥0,λ3≥0.

采用SCA算法求解，可以獲得3個(gè)最小范數(shù)解同時(shí)也是稀疏解為x*=(1,0,0)T,(0,1,0)T,(0,0,1)T.

(iv)考慮線性互補(bǔ)問題

對(duì)于?x=(x1,x2,x3,x4,x5)T，由于

xTMx=3x1x3+5x1x4+3x1x5+3x2x3+3x2x4+5x2x5，

4 結(jié)束語

最小范數(shù)解與稀疏解具有一定的內(nèi)在關(guān)系，但是二者是否等價(jià)難以判斷.目前研究結(jié)果證明了在有限等距條件下最小范數(shù)解和稀疏解的等價(jià)性.因此通過研究最小范數(shù)解，可以為研究稀疏解提供一些近似的結(jié)果.