混合整數非線性規劃的自適應變異差分進化算法

車林仙， 何兵， 劉永波

混合整數非線性規劃的自適應變異差分進化算法

車林仙1,2， 何兵2,3a， 劉永波3b

(1.重慶工程職業技術學院機械工程學院, 重慶402260;2.人工智能四川省重點實驗室, 四川自貢643000;3.瀘州職業技術學院a.機械工程系;b.信息工程系, 四川瀘州646005)

提出了一種求解混合整數非線性規劃(Mixed integer nonlinear programming，MINLP)的混合差分進化(Differential evolution，DE)算法。為提高DE算法的優化性能，設計了混沌初始化種群、可平衡全局探索與精細開采能力的混合變異版本、基于種群進化停滯代數記錄的自適應二次變異算子等新型策略。將前述策略融入DE算法，形成面向MINLP的自適應變異差分進化(Adaptive mutation differential evolution，AMDE)算法。6個MINLP數值實例的對比實驗表明了新算法的有效性和可靠性。最后，應用AMDE算法求解齒輪傳動體積最小化工程優化設計實例，顯示了該算法的工程應用價值。

混合整數非線性規劃；差分進化算法；人工蜂群算法；自適應變異；齒輪傳動優化設計

引言

工程設計、系統工程和經濟管理等領域存在許多混合整數非線性規劃(Mixed integer nonlinear programming，MINLP)問題，它同時包含連續變量和整數變量，其目標函數和約束條件通常呈現強烈非線性。MINLP的求解方法主要有兩類。一類為確定性方法[1]，如分支定界算法、外逼近算法、割平面算法和廣義Bender分解方法等混合變量規劃方法。但這類方法的求解過程通常較復雜，并且對數學理論基礎有較高要求。另一類為智能優化算法，如遺傳算法(Genetic algorithm，GA)、進化策略(Evolution strategy，ES)、粒子群優化(Particle swarm optimization，PSO)算法和差分進化(Differential evolution，DE)算法等。

很多學者針對MINLP，提出了多種智能優化方法。張甲江等[2]、譚躍等[3]分別提出適于MINLP的量子粒子群優化算法和基于交叉與多混沌策略的粒子群優化算法。張莉等[4]、鄧長壽等[5]、譚躍等[6]和擺亮等[7]先后設計多種面向MINLP的混合差分進化算法，如嵌入正交雜交算子的差分進化算法[4]、交替變異版本的差分進化算法[5]、混沌局部搜索差分進化(DE with chaotic local search strategy，CLSDE)算法[6]和差分進化分布估計混合算法[7]。文獻[2-6]僅進行了經典函數優化問題測試，并未給出工程應用實例。因此，面向MINLP的智能算法，尚未在工程優化領域顯示出應用價值。

為提高DE算法求解MINLP的精度和可靠性，本文提出一種自適應變異差分進化(Adaptive mutation differential evolution，AMDE)算法。該算法采用混沌序列初始化種群，融合人工蜂群(Artificial bee colony，ABC)算法[8-10]搜索方式和動態調整變異概率的二次變異，并采用可行性規則[11]處理約束條件。最后，應用AMDE算法求解6個MINLP數值實例和1個二級斜齒圓柱齒輪傳動優化設計案例，得到了滿意結果。

1 MINLP數學模型描述

不失一般性，將最小化MINLP問題表示為

(1)

式中，x為連續決策變量，x=(x1,x2,…,xD1)T；y為整數決策變量，y=(y1,y2,…,yD2)T；xU,d,xL,d為各維連續決策變量的上、下界；yU,d,yL,d為各維整數決策變量的上、下界；D1,D2為連續、整數決策變量的維數；f(x,y)是目標函數；gi(x,y)為不等式約束；hi(x,y)為等式約束；I1,I2為不等式、等式約束數。

為便于敘述，將決策矢量(x,y)統一記為混合整數決策矢量：

z=(z1,z2,…,zD1+D2)T= (x1,x2,…,xD1,y1,y2,…,yD2)T

并將z的各維分量上下界記為zU,d和zL,d(d=1, 2, …,D1+D。

應用DE算法求解MINLP式(1)時，涉及候選解(個體)比較與選擇，本文應用Deb提出的可行性規則[11]比較兩個候選解。該規則為：(1)任何可行解都優于不可行解；(2)對于兩個可行解，目標函數值較小者為優；(3)對于兩個不可行解，約束違反度較小者為優。為此，將任意候選解z的約束違反度函數fvio(z)定義為：

2 混合變量AMDE算法

DE算法的4個控制參數為：種群規模NP，縮放因子F，交叉因子CR和最大進化代數tmax。

應用DE算法求解式(1)時，將進化至第t代的第n個個體編碼記為：

t=0, 1, …,tmax；n=1, 2, …,NP

2.1 混沌初始化

Kent映射生成的混沌序列分布均勻，且具有良好遍歷性，本文應用該混沌序列生成初始種群。一維Kent映射可表示為：

d=1, 2,…,D1+D2；k=0, 1, …，NI

(2)

式中：μ0∈(0, 1)，為常數，且μ0≠0.5，本文取μ0=0.4；q0,d為(0, 1)內服從均勻分布的隨機量，且滿足q0,d≠μ0；NI為混沌序列個數。

由式(2)生成NI個(D1+D維混沌序列點之后，再按下式載波為NI個初始個體

在NI個初始個體中，應用Deb規則選擇較優的NP個體作為初始種群。

2.2 自適應變異概率

基本DE算法的種群多樣性將隨進化代數的增大而降低，求解某些多極值MINLP問題時，可能出現進化停滯和早熟收斂現象。為增強算法的全局優化能力，克服該缺陷，有必要對新生成的試驗個體執行二次變異操作。當連續停滯代數較小時，宜設置較小的變異概率，避免過多無用操作而降低搜索效率；當連續停滯代數較大時，宜增大變異概率，利于保持種群多樣性，增大獲得全局最優解的概率。

定義1對于足夠小的精度閾值δ(t)(隨進化代數改變)，若滿足

(3)

或

(4)

則認為種群在第t代進化停滯。根據實踐經驗，取δ(0)=10-5，δ(t)=δ(t-1)/1.035時，效果較理想。

若連續兩代的最優個體均為可行解，則按式(3)判斷；若連續兩代的最優個體均為不可行解，則按式(4)判斷；若連續兩代的最優個體分別為不可行解和可行解，則認為種群在第t代進化未停滯。

算法設置計數器tst記錄種群連續停滯代數，隨著tst的變化，自適應調整變異概率pm，其計算式為：

(5)

由式(5)可知，當tst較小時，pm也較小；當tst增大時，pm也隨之增大，但最大不超過0.1。因pm過大，將加大隨機搜索成分，反而降低搜索效率。

2.3 計算流程

計算流程各步驟及生成新解時，參照2.1節。

(1) 初始化

(2) 變異操作

DE算法的關鍵在于先應用差分算子生成變異個體，再執行交叉操作生成試驗個體。應用DE/rand/1策略生成變異個體時，算法具有較強全局探索能力，但局部精細開采能力偏弱。為提高收斂速度和精度，本文借鑒文獻[10]中改進ABC算法的引領蜂搜索方式，引入一種基于最優個體的變異策略，記為DE/ rand/best-rand。該策略仍以隨機個體為基向量，但差分向量指向最優個體，有利于精細開采且不易陷入局部最優區域。為了在全局探索和精細開采之間取得平衡，本文算法每次生成變異個體時，均在上述兩種策略中隨機選擇一種。現將兩種策略的變異算子描述如下：

(i) DE/rand/1策略

(6)

(ii) DE/rand/best-rand策略

式中，R(D1+D為 (-1, 1)內服從均勻分布的(D1+D維隨機數；?表示點對點乘法運算符。

若u之某維分量越界，則在其取值區間內隨機生成該維分量。

(3) 交叉操作

式中，rand( )表示(0,1)內服從均勻分布的隨機數；dr為隨機自然數，且dr∈[1,D]∩Z；CR為交叉因子，且CR∈[0,1]。

(4) 二次變異操作

式中，M(·)為二次變異算子。

為增強算法的適應性，文中對連續變量和整數變量，分別執行非均勻變異和均勻變異操作。即，當vd(d=1, 2,…,D1)為連續變量時，二次變異算子為：

式中，r表示在(0,1)內服從均勻分布的隨機數。

當vd(d=D1+1,D1+2,…,D1+D為整數變量時，又分為兩種情形：若為0-1整數規劃問題，二次變異算子為：

M(vd)=1-vd

若為一般整數規劃問題，二次變異算子為：

M(vd)=zL,d+round[rand( )·(zU,d-zL,d)]

(5) 選擇操作

對n=1,2,…,NP，重復執行NP次步驟(2)～(5)。

(6) 更新最優個體

(7) 終止判斷

2.4 計算復雜度

可據第2.3節的計算流程分析AMDE算法的計算復雜度。步驟(1)的計算復雜度為O(NID)，步驟(2)～(4)的計算復雜度均為O(NPD)；步驟(5)的計算復雜度為O(NP)；步驟(6)的計算復雜度：搜索NP個個體中的最優個體為O(NP)，更新連續停滯代數計數器為O(1)。求和并忽略低階項，可求得AMDE算法的計算復雜度為O(NPDtmax+NID)。算法獨立運行1次的函數評價次數為NI+NPtmax。

3MINLP算例測試

3.1 測試實例

以下6個測試實例均引自文獻[6]。

P1：minf(x,y)=2x+y

s.t.g1(x,y)=1.25-x2-y≤0

g2(x,y)=x+y-1.6≤0

0≤x≤1.6，y∈{0,1}

已知最優解(x,y;f*)=(0.5, 1; 2)。

P2：minf(x,y)=2x-ln(x/2)-y

s.t.g(x,y)=-x-ln(x/2)+y≤0

0.5≤x≤1.4，y∈{0,1}

已知最優解(x,y;f*)=(1.375, 1; 2.124)。

P3：minf(x,y)=5(x1-0.5)2-0.7y+0.8

s.t.g1(x,y)=-exp(x1-0.2)-x2≤0

g2(x,y)=x2+1.1y+1≤0

g3(x,y)=x1-1.2y-0.2≤0

0.2≤x1≤1，-2.22554≤x2≤-1，y∈{0,1}

已知最優解(x,y;f*)=(0.94194, -2.1, 1; 1.07654)。

P4：minf(x,y)=

s.t.g(x,y)=

0≤x≤10，y∈{0,1}

已知最優解(x,y;f*)=(3.51424, 1; 99.23964)。

P5：minf(x,y)=(x1-1)2+(x2-2)2+(x3-3)2+

(y1-1)2+(y2-1)2+(y3-1)2-ln(y4+1)

s.t.g1(x,y)=x1+x2+x3+y1+y2+y3-5≤0

g3(x,y)=x1+y1-1.2≤0

g4(x,y)=x2+y2-1.8≤0

g5(x,y)=x3+y3-2.5≤0

g6(x,y)=x1+y4-1.2≤0

0≤x1≤1.2，0≤x2≤1.8，0≤x3≤2.5

y1,y2,y3,y4∈{0,1}

已知最優解(x,y;f*)=(0.2,1.280624,1.954483,1,0,0,1;3.557463)。

37.29329y1+40792.141

s.t.g1(x,y)=a1+a2x3y2+a3x2y1-a4x1x3- 92≤0

g3(x,y)=a9+a10x1x3+a11x1y1+a12x1x2-25≤0

27≤x1,x2,x3≤45，y1∈[78,102]∩Z，

y2∈[33,45]∩Z

已知最優解(x,y;f*)=(27,*,27,78,*;32217.42778) (*表示最優解不確定)。對應系數取值見表1 。

表1P6中的系數取值

3.2 實驗設置

AMDE算法控制參數：NP和tmax取值見表2 ，NI= 100NP，F=0.5+ 0.05rand( )，CR=0.9。

表2 AMDE算法的控制參數

為考察AMDE算法求解測試實例的優化性能，文中還以PSO算法，人工蜂群(Artificial bee colony，ABC)算法[11-12]和基本DE算法(第2.3節的算法流程中，步驟(2)僅用式(6)生成變異個體，且不執行步驟(4))等作為對比算法，進行比較分析。

為公平起見，3種對比算法均采用與AMDE算法相同的初始化、約束函數、取整和越界處理方法；對比算法的NP和tmax均與AMDE算法一致。DE算法的F,CR與AMDE算法一致。

PSO算法的慣性權重w線性遞減，取w=1.2～0.2，加速度系數為c1=c2=2，最大速度取vmax,d= 0.3(zU,d-zL,d) (d=1, 2, …,D1+D。

ABC算法的引領蜂和跟隨蜂各占種群規模的50%；跟隨蜂采用2元錦標賽方法選擇被跟隨的引領蜂；引領蜂和跟隨蜂執行鄰域搜索生成新蜜源時，均以概率CR隨機改變多維分量。

3.3 實驗結果

4種算法分別獨立運行50次的統計結果見表3 。其中，fb,fav,fw和σf分別表示最優目標函數值的最好、平均、最差值及標準差。同時，將CLSDE算法[9]的性能指標也列入表5 ，以便比較。

表3 5種算法求解實例的性能指標

為了便于應用對數坐標清晰顯示各種算法的收斂精度，首先給出關于算法種群最小誤差的概念。

定義2設待求MINLP問題的理論最優目標函數值為f*，定義第t代種群的最小誤差為：

4種算法的平均最小誤差eb,av(圖中取對數)進化曲線對比如圖1 所示。

由表3 和圖1 可知，求解P1、P4、P5和P6時，PSO算法的eb,av最大，求解P2和P3時，PSO算法的eb,av與DE算法相當，但大于ABC和AMDE算法；求解6個實例時，PSO算法的σf均最大，穩健性最差。在5種算法中，PSO算法的優化性能最差。AMDE與ABC算法比較，求解P1、P2、P4和P6時，2種算法相當；而求解P3和P5時，AMDE明顯優于ABC算法。AMDE與DE算法比較，求解P1、P4和P6時，2種算法相當；求解P2和P3時，AMDE明顯優于DE算法；求解P5時，AMDE差于DE算法。求解6個實例時，AMDE算法的eb,av和σf均很小，表明其收斂精度高、穩健性好。對于其余4種對比算法，無一算法在求解6個實例時同時具有良好優化性能。由表2 可知，AMDE算法求解實例時的最大函數評價次數為15 000，而CLSDE算法為52 000[6]，因此前者的計算效率更高。綜上可知，AMDE算法的性能指標優于4種對比算法。

4 齒輪傳動優化設計應用

4.1 齒輪傳動優化設計模型

已知某企業使用的二級斜齒圓柱齒輪傳動原始數據為：高速軸輸入功率P1=8 kW；高速級主動軸轉速n1=1450 r/min；傳動比itr=9.7，允許傳動比誤差±5%；動力機為電動機，工作機載荷均勻平穩。大、小齒輪材料為45鋼，大齒輪正火處理，硬度185～205 HBS；小齒輪調質處理，硬度為230～255 HBS。許用應力[σH]I=[σH]II= 520 MPa，[σF]I=155 MPa，[σF]II=145 MPa。載荷系數K=1.25。要求以齒輪總體積最小為目標，設計該齒輪傳動機構。

圖1 平均最小誤差進化曲線

該問題的混合設計變量為：第I級傳動，法面模數mnI、小齒輪齒數z1、螺旋角βI、齒寬系數ψdI和傳動比iI；第II級傳動，法面模數mnII、小齒輪齒數z3、螺旋角βII和齒寬系數ψdII。其中，mnI和mnII為非等間隔離散變量，z1和z3為整數變量；βI、βII、ψdI、ψdII和iI為連續變量。將決策變量記為x=(x1,x2,…,x5)T=(βI,βII,ψdI,ψdII,iI)T，y=(y1,y2,y3,y4)T=(mnI,mnII,z1,z3)T?？筛鶕覙藴?，建立模數的離散變量取值列表，再映射為整數變量[12-14]。故，該問題可轉化為等效MINLP，再應用第2節的方法求解。

若以齒輪分度圓柱體積近似表示齒輪體積，則總體積最小化設計問題的目標函數為：

minf(x,y)=

式中，z2表示第I級傳動大齒輪齒數，取z2=round(x5y3)；z4表示第II級傳動大齒輪齒數，取z4=round(itry3y4/z。

參照文獻[12]，列出以下約束函數。其中，g1～g4為強度條件，g5為避免第II級大齒輪與第I級傳動軸發生干涉的限制條件，g6～g7為最小齒數約束，g8～g9為軸向重合度約束，g10為傳動比精度約束，g11～g14為齒寬約束。

式中，T1,T3為第I,II級傳動計算轉矩，N·mm，且T1=9.55×106P1/n1，T3=9.55×106x5P1/n1；ZEI,ZEII為第I,II級傳動彈性系數；ZHI,ZHII為第I, II級傳動節點區域系數；ZεI,ZεII為第I, II級傳動重合度系數；ZβI,ZβII為第I, II級傳動螺旋角系數；YεI,YεII為第I, II級傳動重合度系數；YβI,YβII為第I, II級傳動螺旋角系數；YFSI,YFSII為第I, II級傳動復合齒形系數。

以上系數的確定方法見文獻[15]，其中，將YFSI,YFSII表達為當量齒數的擬合函數。

g5(x,y)=y1z2cosx2+2(y1+50)cosx1cosx2-y2(y4+z4)cosx1≤0

g6(x,y)=17cos3x1-y3≤0

g7(x,y)=17cos3x2-y4≤0

g8(x,y)=π-x3y3tanx1≤0

g9(x,y)=π-x4y4tanx2≤0

g11(x,y)=40-x3y1y3secx1≤0

g12(x,y)=x3y1y3secx1-120≤0

g13(x,y)=40-x4y2y4secx2≤0

g14(x,y)=x4y2y4secx2-120≤0

設計變量的取值范圍設置為：0.139 6 rad≤x1,x2≤0.261 8 rad，0.6≤x3,x4≤1.2，2≤x5≤4，y1∈{2, 2.5, 2.75, 3, 3.25, 3.5, 3.75, 4, 4.5, 5} mm，y2∈{3.5, 3.75, 4, 4.5, 5, 5.5, 6} mm，y3,y4∈[15, 30]∩Z。

4.2 優化設計結果

應用AMDE算法求解該設計實例，除取NP=100及tmax=1500外，其余參數設置同第3.2節。算法獨立運行10次，最好、最差和平均最優目標函數值分別為2.7623 dm3, 2.7758 dm3和2.7690 dm3，標準差為0.003 949 dm3。該標準差很小，表明算法具有良好穩健性。平均最優值進化曲線如圖2 ，由圖可知，算法大約進化1200代即收斂于近優解。

圖2 AMDE算法求解案例的平均最優值進化曲線

AMDE算法的最優解見表4 ，將原設計方案也列入表中，以便比較。由表4 可知，應用AMDE算法求得的結果明顯優于原設計方案，齒輪體積和較原設計下降約25.5%。

表4 優化設計與原方案結果比較

5 結論

(1)為提高DE算法求解MINLP的優化性能，設計混沌初始化種群、可平衡全局探索與精細開采能力的混合變異版本、基于種群進化停滯代數記錄的自適應二次變異算子等新型策略。將諸策略融入DE算法，形成求解MINLP的混合智能算法AMDE。

(2)應用6個MINLP數值實例測試AMDE算法，結果表明：AMDE算法可靠、有效，穩健性優于4種對比算法。

(3)以二級斜齒圓柱齒輪傳動體積最小化設計問題，驗證AMDE算法求解混合變量工程優化問題的可行性。根據映射關系，將混合離散變量約束優化問題轉化為等效MINLP問題，再應用AMDE算法求解該問題。結果顯示：AMDE算法可行、有效，具有良好穩健性。與原設計方案相比，AMDE算法的最優解下降了約25.5%。

[1]劉明明,崔春風,童小嬌,等.混合整數非線性規劃的算法軟件及最新進展.中國科學:數學,2016,46(1):1-20.

[2]張甲江,高岳林,高晨陽.非線性混合整數規劃問題的改進量子粒子群算法.太原理工大學學報,2015,46(2):196-200.

[3]譚躍,譚冠政,鄧曙光.基于遺傳交叉和多混沌策略改進的粒子群優化算法.計算機應用研究,2016,33(12):3643-3647.

[4]張莉,李宏,馮大政.求解混合整數規劃的嵌入正交雜交的差分進化算法.系統工程與電子技術,2011,33(9):2126-2132.

[5]鄧長壽,任紅衛,彭虎.混合整數非線性規劃問題的改進差分進化算法.計算機應用研究,2012,29(2):445-448.

[6]譚躍,譚冠政,楊冰,等.解決混合整數非線性規劃問題的混沌局部搜索差分進化算法.小型微型計算機系統,2012,33(6):1306-1309.

[7]BAI L,WANG J Y,JIANG Y H,et al.Improved hybrid differential evolution-estimation of distribution algorithm with feasibility rules for NLP/MINLP engineering optimization problems.Chinese Journal of Chemical Engineering,2012,20(6):1074-1080.

[8]AKAY B,KARABOGA D.A modified artificial bee colony algorithm for real-parameter optimization.Information Sciences,2012,192:120-142.

[9]GAO W F,LIU S Y,HUANG L L.A global best artificial bee colony algorithm for global optimization.Journal of Computational and Applied Mathematics,2012,236(11):2741-2753.

[10]羅鈞,林于晴,劉學明,等.改進蜂群算法及其在圓度誤差評定中的應用.機械工程學報,2016,52(16):27-32.

[11]MALLIPEDDI R,SUGANTHAN P N.Ensemble of constraint handling techniques.IEEE Transactions on Evolutionary Computation,2010,14(4):561-579.

[12]車林仙.面向機構分析與設計的差分進化算法研究.徐州:中國礦業大學,2012.

[13]車林仙,程志紅.工程約束優化的自適應罰函數混合離散差分進化算法.機械工程學報,2011,47(3):141-151.

[14]車林仙.多樣性保持離散差分進化算法及齒輪傳動優化應用.機械工程學報,2016,52(21):44-55.

DOI:10.11863/j.suse.2017.01.05

Adaptive Mutation Differential Evolution Algorithm for Mixed Integer Nonlinear Programming

CHELinxian1,2,HEBing2,3a,LIUYongbo3b

(1.School of Mechanical Engineering, Chongqing Vocational Institute of Engineering, Chongqing 402260, China;2.Artificial Intelligence Key Laboratory of Sichuan Province, Zigong 643000, China; 3a.Department of Mechanical Engineering; 3b.Department of Information Engineering, Luzhou Vocational and Technical College,Luzhou 646005, China)

A hybrid differential evolution (DE) algorithm is presented to solve the mixed integer nonlinear programming (MINLP). To enhance the optimization performance of DE algorithm, several new strategies are designed, such as chaotic initializing population, hybrid mutation scheme which can balance global exploration and fine mining abilities, and adaptive second mutation operator based on the recording generations of evolutionary stagnation for a population. The aforementioned strategies are embedded in DE algorithm and an adaptive mutation DE (AMDE) algorithm is formed for solving MINLP. Six numerical examples of MINLP are tested comparatively to show that this new approach is valid and reliable. Finally the proposed AMDE algorithm is used to solve a real case of engineering optimization for minimizing volumes of gear transmission, so the practical applicability of the new algorithm is indicated.

mixed integer nonlinear programming; differential evolution algorithm; artificial bee colony algorithm; adaptive mutation operator; optimal design of gear transmission

2016-10-24

重慶市基礎科學與前沿技術研究專項資助項目(cstc2015jcyjA70006);重慶市教育委員會科學技術研究項目(KJ1403201);人工智能四川省重點實驗室開放基金項目(2013RYJ02)

車林仙(1971-),男,四川瀘州人,教授,博士,主要從事機械優化設計、機構學及智能信息處理方面的研究,(E-mail)lx.che@163.com

1673-1549(2017)01-0019-08

10.11863/j.suse.2017.01.04

TP18;O221.4

A