混料D-最優近似設計的改進模擬退火算法

許秀鈿,張崇岐

(廣州大學 經濟與統計學院,廣東廣州 510006)

§1 引言

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

通過近年來混料的應用研究[6-7]可以發現,如何科學的設計試驗方案,以較少的工作量獲得最優的試驗結果是當前混料設計的研究熱點.由于傳統的計算方法對于求解高階混料模型或是約束域較為復雜的混料模型難度較大,利用優化算法求解混料最優設計[8]逐漸成為該領域關注的對象.

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

§2 預備知識

最基本的混料模型為

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

假設設計ξ由設計點zi和對應的測度ri組成,

其中r1+r2+···+rN=1,Ξ是所有設計構成的集合.

設計ξ的信息矩陣[15]定義為

稱信息矩陣的行列式最大化即max detM(ξ)的設計ξ為D-最優設計.

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

即在無附加約束的混料問題上,試驗的可行域是一個q-1維正規單純形

在實際的問題中,有時出于經濟與技術方面的考慮,混料試驗往往還需滿足其他約束條件,常見的附加約束類型[16]有

(1) 同時具有上,下界約束條件(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都是已知的常數,φi(x1,x2,···,xq)是各變量的非線性函數.

1958年,Scheffé首次提出單純形格子設計,奠定了混料試驗設計的基礎.目前許多設計方法是由該設計衍生出來的,其中混料規范多項式應用較為廣泛.在尋找最優設計的過程中,常使用的q分量m階混料規范多項式有

(1)m=1時,

(2)m=2時,

(3)m=3時,

顯然當q或m增加時,模型的參數個數也會隨之增加.

如果試驗設計的試驗點數與相應的混料規范多項式模型中未知參數的個數相等,則設計為飽和設計.易知飽和設計是試驗點數最少情況下的最優設計.

§3 改進SA算法

在實際的混料最優設計的問題求解中,由于模型參數個數較多或是約束域較為復雜,模型的D-最優設計的解析解往往較難求出.此時借助計算機編程算法求解是一種很好的解決方法.

由于D-最優設計具有等測度的性質,故求解D-最優近似最優設計的過程也是求解測度為參數個數的倒數的最優支撐點的過程.

模擬退火算法是借鑒固體退火過程原理設計而成的算法.固體溫度很高時,內能很大,固體內部粒子的無規則運動很快,此時整個試驗域的點都有一定概率得以取到.隨著溫度參數慢慢降低,內能逐漸減小,固體內部粒子逐漸趨于有序,即算法中的解趨于穩定.但是,這個穩定的解可能只是局部最優解.所以在解逐漸收斂的過程中,模擬退火算法會以一定的概率跳出這個局部最優解.即固體內能增加,取到的解遠離目標函數的時候,算法以一定概率的接納,而不是一味拒絕.這也是模擬退火算法的核心思想.

3.1 改進SA算法的實現過程

基于改進的SA算法求解混料模型最優設計的基本步驟如下.

1.設定初始值.設置迭代的初始設計點,為確保迭代能遍歷整個可行域,初始設計點一般隨機生成并均勻分布在整個可行域上.

2.產生隨機擾動項.由一個產生函數在初始設計點附近產生一個位于試驗可行域上的新設計點.本文提出的改進的SA算法相比傳統的SA算法,考慮了隨機擾動項波動范圍的處理,將步長設置為溫度的平方根.溫度越低,得到的解約趨于穩定時,步長越短,波動范圍越小,能更好的提高搜索設計點的效率.

3.判斷新解是否能被接受.最常用的判斷依據是Metropolis準則,即當新設計點的目標函數小于當前設計點的目標函數時,接受新設計點作為新的當前設計點;否則則以一定的概率接受新設計點為新的當前設計點.

4.生成最優設計點集.當“溫度”降到“冷卻溫度”時,迭代停止,最終保留下來的目標函數最小的設計點即為所得最優支撐點.本文提出的改進SA算法將馬爾科夫鏈長度設置為溫度的函數,使得溫度越低時,馬爾科夫鏈越長.相較于一般的SA算法,極大的提高了搜索和迭代的效率.

將最優支撐點與其對應的測度組合即可得到該問題的D-最優近似設計.此外,對于求解出的混料模型,也可應用該算法求出最優解.

3.2 改進SA算法求解D-最優近似設計的工作流程

第1步 初始化

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

k=L/t為每一溫度的迭代次數,顯然溫度降低時,迭代的次數將會增加.

第2步 產生隨機擾動項,隨機生成一個領域N(i)

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

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

第5步 若達到熱平衡,即某一溫度的迭代結束時,則繼續降溫,直至降為最終冷卻溫度,否則返回第2步

第6步 生成最優設計點集

將最終冷卻溫度時保留的最優支撐點集與取值為參數個數倒數的測度組合即可得到所求D-最優近似設計.

將最優支撐點集進行試驗,得到相應的響應值.代入混料模型中求解出各參數估計值,即可得到經驗回歸方程.通過改進的SA算法,將上述步驟中目標函數改為經驗回歸方程,初始化時只隨機抽取1個設計點,便可類似求解出試驗的最優設計點.

3.3 改進SA算法求解的偽代碼

§4 混料試驗的應用:紗線伸長度[18]

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

此次試驗采用三分量二階混料規范多項式模型,其表達式為

其可行域為S3-1,分別利用改進SA算法和傳統的SA算法求出其D-最優近似設計,結果如表1和圖1所示,圖1中改進的SA算法得出的設計點標記為三角形,傳統的SA算法得出的設計點標記為正方形,理論上的設計點標記為圓形.

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

表1 改進的SA算法與傳統的SA算法得出的D-最優近似設計解析解比較

易知此時參數個數為6,當給定初始溫度T0=100,最終冷卻溫度Tf=1,溫度衰減速率α=0.99,馬爾科夫鏈長度L=1000000時,2種算法的三分量二階混料規范多項式模型D-最優近似飽和設計的結果如表1所示.對于三分量二階混料規范多項式,D-最優設計的理論解[19-20]是單純形上3個頂點與3個邊界中點的組合,即

因此從圖1中的結果可知,對于D-近似最優設計而言,用改進的SA算法比傳統的SA算法求出的D-最優近似設計數值解跟理論解更為接近.

從圖2中可以看到,隨著迭代次數的增加,溫度逐漸降低.改進的模擬退火算法在目標函數值降為-4.7162×10-9之后達到平衡,此時代碼運行總時間為3233.792秒;傳統的模擬退火算法在目標函數值降為-3.9007×10-11之后達到平衡，此時代碼運行總時間為13465.094秒.可以發現改進SA算法的目標函數值更低,即由改進的SA算法得到的設計較傳統SA算法更優.同時,根據這兩種算法的代碼運行時間可知,改進的SA算法相比傳統的SA 算法求解速度更快.

圖2 改進(左圖)和傳統(右圖)的SA算法的D-最優近似設計目標函數值與迭代情況

本次試驗的部分試驗數據和設計點的響應平均值分別如表2和圖3所示.

圖3 設計點的響應平均值

表2 紗線伸長度部分試驗數據

通過最小二乘法即可算出本例所擬合的經驗回歸方程為

該模型的響應等值線圖如圖4所示.

圖4 紗線伸長度的等值線圖

從算法運行結果可以看出,通過改進的SA算法得到最優設計點為(0.3014,0.0082,0.6904),其對應的響應為17.2905,如圖5所示.

圖5 響應目標函數值與迭代情況

§5 結束語

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