許秀鈿,張崇岐

(廣州大學(xué) 經(jīng)濟(jì)與統(tǒng)計(jì)學(xué)院,廣東廣州 510006)

§1 引言

混料試驗(yàn)設(shè)計(jì)自1958年Scheffé[1]提出以來發(fā)展十分迅速,當(dāng)前已被廣泛應(yīng)用于工業(yè),食品,醫(yī)藥,金融等領(lǐng)域.例如,工業(yè)領(lǐng)域,孫世昌等[2]在2021年利用混料試驗(yàn)設(shè)計(jì)優(yōu)化煙用膠囊壁材配方,提高了工業(yè)生產(chǎn)的煙用膠囊質(zhì)量;食品領(lǐng)域,張麗霞等[3]在2021年通過最優(yōu)混料試驗(yàn)設(shè)計(jì)優(yōu)化風(fēng)味蛋白酶,堿性蛋白酶和復(fù)合蛋白酶的組合,為花生餅粕制備花生呈味基料提供一種高效復(fù)合酶;醫(yī)藥領(lǐng)域,張雪等[4]在2022年通過混料D-最優(yōu)設(shè)計(jì)優(yōu)選出適宜的擠出搓圓法制備丸劑的復(fù)合輔料;金融領(lǐng)域,燕飛和張崇岐[5]用混料模型G-最優(yōu)設(shè)計(jì)來研究證券投資組合.

通過近年來混料的應(yīng)用研究[6-7]可以發(fā)現(xiàn),如何科學(xué)的設(shè)計(jì)試驗(yàn)方案,以較少的工作量獲得最優(yōu)的試驗(yàn)結(jié)果是當(dāng)前混料設(shè)計(jì)的研究熱點(diǎn).由于傳統(tǒng)的計(jì)算方法對(duì)于求解高階混料模型或是約束域較為復(fù)雜的混料模型難度較大,利用優(yōu)化算法求解混料最優(yōu)設(shè)計(jì)[8]逐漸成為該領(lǐng)域關(guān)注的對(duì)象.

模擬退火算法[9](Simulated Annealing)又稱SA算法,是在1953年由Metropolis提出的一種理論上的全局優(yōu)化算法.在1983年Kirkpatrick[10]成功將其應(yīng)用到求解組合優(yōu)化問題上.當(dāng)前應(yīng)用較為廣泛,比較具有代表性的智能優(yōu)化算法[11]是遺傳算法[12],禁忌搜索算法,模擬退火算法[13-14].與其他種類的優(yōu)化算法相比,模擬退火算法理論較為完善.本文介紹一種改進(jìn)的SA算法,該算法保留了SA算法易于求出全局最優(yōu)解的優(yōu)點(diǎn),同時(shí)在求解精度方面相比一般的SA算法有一定程度的提升.

§2 預(yù)備知識(shí)

最基本的混料模型為

其中y是響應(yīng)變量,β是未知參數(shù)組成的d維列向量,f(x)是在x ∈χ上給定的d維函數(shù)列向量,χ是一個(gè)封閉緊集,ε為隨機(jī)誤差,設(shè)其服從均值為0,方差為常數(shù)的正態(tài)分布.

假設(shè)設(shè)計(jì)ξ由設(shè)計(jì)點(diǎn)zi和對(duì)應(yīng)的測(cè)度ri組成,

其中r1+r2+···+rN=1,Ξ是所有設(shè)計(jì)構(gòu)成的集合.

設(shè)計(jì)ξ的信息矩陣[15]定義為

稱信息矩陣的行列式最大化即max detM(ξ)的設(shè)計(jì)ξ為D-最優(yōu)設(shè)計(jì).

混料試驗(yàn)是研究各個(gè)分量比例x1,x2,···,xq的設(shè)計(jì)問題,q分量混料系統(tǒng)各種成份的占比xi(i=1,2,···,q)必須滿足約束條件

即在無附加約束的混料問題上,試驗(yàn)的可行域是一個(gè)q-1維正規(guī)單純形

在實(shí)際的問題中,有時(shí)出于經(jīng)濟(jì)與技術(shù)方面的考慮,混料試驗(yàn)往往還需滿足其他約束條件,常見的附加約束類型[16]有

(1) 同時(shí)具有上,下界約束條件(SCCs):C′s:0≤ai ≤xi ≤bi ≤1,i=1,2,···,q;

(2) 具有線性約束條件(MCCs):C′s:0≤ai ≤≤bi ≤1,i=1,2,···,q;

(3) 具有非線性約束條件(NCCs):C′s:ai ≤φi(x1,x2,···,xq)≤bi ≤1,i=1,2,···,q;其中ai,bi都是已知的常數(shù),φi(x1,x2,···,xq)是各變量的非線性函數(shù).

1958年,Scheffé首次提出單純形格子設(shè)計(jì),奠定了混料試驗(yàn)設(shè)計(jì)的基礎(chǔ).目前許多設(shè)計(jì)方法是由該設(shè)計(jì)衍生出來的,其中混料規(guī)范多項(xiàng)式應(yīng)用較為廣泛.在尋找最優(yōu)設(shè)計(jì)的過程中,常使用的q分量m階混料規(guī)范多項(xiàng)式有

(1)m=1時(shí),

(2)m=2時(shí),

(3)m=3時(shí),

顯然當(dāng)q或m增加時(shí),模型的參數(shù)個(gè)數(shù)也會(huì)隨之增加.

如果試驗(yàn)設(shè)計(jì)的試驗(yàn)點(diǎn)數(shù)與相應(yīng)的混料規(guī)范多項(xiàng)式模型中未知參數(shù)的個(gè)數(shù)相等,則設(shè)計(jì)為飽和設(shè)計(jì).易知飽和設(shè)計(jì)是試驗(yàn)點(diǎn)數(shù)最少情況下的最優(yōu)設(shè)計(jì).

§3 改進(jìn)SA算法

在實(shí)際的混料最優(yōu)設(shè)計(jì)的問題求解中,由于模型參數(shù)個(gè)數(shù)較多或是約束域較為復(fù)雜,模型的D-最優(yōu)設(shè)計(jì)的解析解往往較難求出.此時(shí)借助計(jì)算機(jī)編程算法求解是一種很好的解決方法.

由于D-最優(yōu)設(shè)計(jì)具有等測(cè)度的性質(zhì),故求解D-最優(yōu)近似最優(yōu)設(shè)計(jì)的過程也是求解測(cè)度為參數(shù)個(gè)數(shù)的倒數(shù)的最優(yōu)支撐點(diǎn)的過程.

模擬退火算法是借鑒固體退火過程原理設(shè)計(jì)而成的算法.固體溫度很高時(shí),內(nèi)能很大,固體內(nèi)部粒子的無規(guī)則運(yùn)動(dòng)很快,此時(shí)整個(gè)試驗(yàn)域的點(diǎn)都有一定概率得以取到.隨著溫度參數(shù)慢慢降低,內(nèi)能逐漸減小,固體內(nèi)部粒子逐漸趨于有序,即算法中的解趨于穩(wěn)定.但是,這個(gè)穩(wěn)定的解可能只是局部最優(yōu)解.所以在解逐漸收斂的過程中,模擬退火算法會(huì)以一定的概率跳出這個(gè)局部最優(yōu)解.即固體內(nèi)能增加,取到的解遠(yuǎn)離目標(biāo)函數(shù)的時(shí)候,算法以一定概率的接納,而不是一味拒絕.這也是模擬退火算法的核心思想.

3.1 改進(jìn)SA算法的實(shí)現(xiàn)過程

基于改進(jìn)的SA算法求解混料模型最優(yōu)設(shè)計(jì)的基本步驟如下.

1.設(shè)定初始值.設(shè)置迭代的初始設(shè)計(jì)點(diǎn),為確保迭代能遍歷整個(gè)可行域,初始設(shè)計(jì)點(diǎn)一般隨機(jī)生成并均勻分布在整個(gè)可行域上.

2.產(chǎn)生隨機(jī)擾動(dòng)項(xiàng).由一個(gè)產(chǎn)生函數(shù)在初始設(shè)計(jì)點(diǎn)附近產(chǎn)生一個(gè)位于試驗(yàn)可行域上的新設(shè)計(jì)點(diǎn).本文提出的改進(jìn)的SA算法相比傳統(tǒng)的SA算法,考慮了隨機(jī)擾動(dòng)項(xiàng)波動(dòng)范圍的處理,將步長(zhǎng)設(shè)置為溫度的平方根.溫度越低,得到的解約趨于穩(wěn)定時(shí),步長(zhǎng)越短,波動(dòng)范圍越小,能更好的提高搜索設(shè)計(jì)點(diǎn)的效率.

3.判斷新解是否能被接受.最常用的判斷依據(jù)是Metropolis準(zhǔn)則,即當(dāng)新設(shè)計(jì)點(diǎn)的目標(biāo)函數(shù)小于當(dāng)前設(shè)計(jì)點(diǎn)的目標(biāo)函數(shù)時(shí),接受新設(shè)計(jì)點(diǎn)作為新的當(dāng)前設(shè)計(jì)點(diǎn);否則則以一定的概率接受新設(shè)計(jì)點(diǎn)為新的當(dāng)前設(shè)計(jì)點(diǎn).

4.生成最優(yōu)設(shè)計(jì)點(diǎn)集.當(dāng)“溫度”降到“冷卻溫度”時(shí),迭代停止,最終保留下來的目標(biāo)函數(shù)最小的設(shè)計(jì)點(diǎn)即為所得最優(yōu)支撐點(diǎn).本文提出的改進(jìn)SA算法將馬爾科夫鏈長(zhǎng)度設(shè)置為溫度的函數(shù),使得溫度越低時(shí),馬爾科夫鏈越長(zhǎng).相較于一般的SA算法,極大的提高了搜索和迭代的效率.

將最優(yōu)支撐點(diǎn)與其對(duì)應(yīng)的測(cè)度組合即可得到該問題的D-最優(yōu)近似設(shè)計(jì).此外,對(duì)于求解出的混料模型,也可應(yīng)用該算法求出最優(yōu)解.

3.2 改進(jìn)SA算法求解D-最優(yōu)近似設(shè)計(jì)的工作流程

第1步 初始化

設(shè)置初始溫度T0,最終冷卻溫度Tf,溫度衰減速率α ∈(0,1),馬爾科夫鏈長(zhǎng)度L,步長(zhǎng)Step.

k=L/t為每一溫度的迭代次數(shù),顯然溫度降低時(shí),迭代的次數(shù)將會(huì)增加.

第2步 產(chǎn)生隨機(jī)擾動(dòng)項(xiàng),隨機(jī)生成一個(gè)領(lǐng)域N(i)

第3步 計(jì)算目標(biāo)值增量?f=f(j)-f(i),其中j ∈N(i)

第4步 判斷新解是否被接受,若?f<0或?f>0且exp(-)>U(0,1),則令i=j

第5步 若達(dá)到熱平衡,即某一溫度的迭代結(jié)束時(shí),則繼續(xù)降溫,直至降為最終冷卻溫度,否則返回第2步

第6步 生成最優(yōu)設(shè)計(jì)點(diǎn)集

將最終冷卻溫度時(shí)保留的最優(yōu)支撐點(diǎn)集與取值為參數(shù)個(gè)數(shù)倒數(shù)的測(cè)度組合即可得到所求D-最優(yōu)近似設(shè)計(jì).

將最優(yōu)支撐點(diǎn)集進(jìn)行試驗(yàn),得到相應(yīng)的響應(yīng)值.代入混料模型中求解出各參數(shù)估計(jì)值,即可得到經(jīng)驗(yàn)回歸方程.通過改進(jìn)的SA算法,將上述步驟中目標(biāo)函數(shù)改為經(jīng)驗(yàn)回歸方程,初始化時(shí)只隨機(jī)抽取1個(gè)設(shè)計(jì)點(diǎn),便可類似求解出試驗(yàn)的最優(yōu)設(shè)計(jì)點(diǎn).

3.3 改進(jìn)SA算法求解的偽代碼

§4 混料試驗(yàn)的應(yīng)用:紗線伸長(zhǎng)度[18]

紡紗線前需要將三種原料:聚乙烯(x1),聚苯乙烯(x2)和聚丙烯(x3)混合制成纖維.固定在作用力下紗線的伸長(zhǎng)度是試驗(yàn)者感興趣的響應(yīng).

此次試驗(yàn)采用三分量二階混料規(guī)范多項(xiàng)式模型,其表達(dá)式為

其可行域?yàn)镾3-1,分別利用改進(jìn)SA算法和傳統(tǒng)的SA算法求出其D-最優(yōu)近似設(shè)計(jì),結(jié)果如表1和圖1所示,圖1中改進(jìn)的SA算法得出的設(shè)計(jì)點(diǎn)標(biāo)記為三角形,傳統(tǒng)的SA算法得出的設(shè)計(jì)點(diǎn)標(biāo)記為正方形,理論上的設(shè)計(jì)點(diǎn)標(biāo)記為圓形.

圖1 兩種算法得出的解與理論解的三元圖

表1 改進(jìn)的SA算法與傳統(tǒng)的SA算法得出的D-最優(yōu)近似設(shè)計(jì)解析解比較

易知此時(shí)參數(shù)個(gè)數(shù)為6,當(dāng)給定初始溫度T0=100,最終冷卻溫度Tf=1,溫度衰減速率α=0.99,馬爾科夫鏈長(zhǎng)度L=1000000時(shí),2種算法的三分量二階混料規(guī)范多項(xiàng)式模型D-最優(yōu)近似飽和設(shè)計(jì)的結(jié)果如表1所示.對(duì)于三分量二階混料規(guī)范多項(xiàng)式,D-最優(yōu)設(shè)計(jì)的理論解[19-20]是單純形上3個(gè)頂點(diǎn)與3個(gè)邊界中點(diǎn)的組合,即

因此從圖1中的結(jié)果可知,對(duì)于D-近似最優(yōu)設(shè)計(jì)而言,用改進(jìn)的SA算法比傳統(tǒng)的SA算法求出的D-最優(yōu)近似設(shè)計(jì)數(shù)值解跟理論解更為接近.

從圖2中可以看到,隨著迭代次數(shù)的增加,溫度逐漸降低.改進(jìn)的模擬退火算法在目標(biāo)函數(shù)值降為-4.7162×10-9之后達(dá)到平衡,此時(shí)代碼運(yùn)行總時(shí)間為3233.792秒;傳統(tǒng)的模擬退火算法在目標(biāo)函數(shù)值降為-3.9007×10-11之后達(dá)到平衡，此時(shí)代碼運(yùn)行總時(shí)間為13465.094秒.可以發(fā)現(xiàn)改進(jìn)SA算法的目標(biāo)函數(shù)值更低,即由改進(jìn)的SA算法得到的設(shè)計(jì)較傳統(tǒng)SA算法更優(yōu).同時(shí),根據(jù)這兩種算法的代碼運(yùn)行時(shí)間可知,改進(jìn)的SA算法相比傳統(tǒng)的SA 算法求解速度更快.

圖2 改進(jìn)(左圖)和傳統(tǒng)(右圖)的SA算法的D-最優(yōu)近似設(shè)計(jì)目標(biāo)函數(shù)值與迭代情況

本次試驗(yàn)的部分試驗(yàn)數(shù)據(jù)和設(shè)計(jì)點(diǎn)的響應(yīng)平均值分別如表2和圖3所示.

圖3 設(shè)計(jì)點(diǎn)的響應(yīng)平均值

表2 紗線伸長(zhǎng)度部分試驗(yàn)數(shù)據(jù)

通過最小二乘法即可算出本例所擬合的經(jīng)驗(yàn)回歸方程為

該模型的響應(yīng)等值線圖如圖4所示.

圖4 紗線伸長(zhǎng)度的等值線圖

從算法運(yùn)行結(jié)果可以看出,通過改進(jìn)的SA算法得到最優(yōu)設(shè)計(jì)點(diǎn)為(0.3014,0.0082,0.6904),其對(duì)應(yīng)的響應(yīng)為17.2905,如圖5所示.

圖5 響應(yīng)目標(biāo)函數(shù)值與迭代情況

§5 結(jié)束語(yǔ)

在給定混料模型和可行域的情況下,通過改進(jìn)的SA算法來求解混料D-最優(yōu)近似設(shè)計(jì)不失為一個(gè)好方法.本文給出了一般情況下的混料規(guī)范多項(xiàng)式求解D-最優(yōu)近似設(shè)計(jì)解析解的改進(jìn)SA算法,通過比較可知,該算法相比傳統(tǒng)的SA算法,求解的速度更快,精度更高.此外,本文提出的改進(jìn)的SA算法也可用于求解其他最優(yōu)設(shè)計(jì)的設(shè)計(jì)點(diǎn),由于篇幅限制,本文僅以D-最優(yōu)設(shè)計(jì)為例進(jìn)行展開.