基于高斯過程模型的在線穩健參數設計

周曉劍，高云龍

（南京郵電大學 管理學院，南京 210023）

0 引言

穩健參數設計（Robust Parameter Design，RPD）是基于試驗設計，結合統計知識和工程技術發展起來的一種質量改進方法[1]。其基本原理為選擇最佳的組合參數水平以使產品及過程對環境噪聲等不確定因素帶來的變化不敏感，從而降低產品質量波動，提高產品及過程設計的穩健性，使產品性能更加穩定可靠[2]。計算機試驗通過構建元模型近似模擬物理試驗，在保證仿真優化精度的同時大大降低了試驗成本。對于元模型，人們進行了大量的研究，其中，多項式模型由于其簡單高效的優勢在模型擬合的使用中最為廣泛，但模型本身的局限性也很明顯，僅適用于低階問題仿真場景[3]。機器學習的蓬勃發展也推動了對元模型的探索，如支持向量機（SVM）[4]、徑向基模型（RBF）[5]、神經網絡模型（NN）[6]等。然而在實際工業生產中，相關數據的輸入輸出越趨復雜，越來越表現出高維性、非線性等特點，以上提到的元模型很容易因數據的高維性造成計算困難，甚至會陷入維數災難。高斯過程（Gaussian Process,GP）模型是一種基于貝葉斯推斷的機器學習方法，適用于處理高維性、非線性等復雜回歸問題[7]。基于以上優良特性，GP模型被廣泛應用于復雜工業過程建模領域，且在穩健參數設計研究中具有優良表現[8,9]。

傳統的穩健參數設計是離線設計，即在生產之前就確定好可控因子的最優水平，且該水平在整個生產過程中并不改變。但在實際的工業生產中，環境因素（噪聲因子）實時變化會給參數設計帶來影響，導致可控因子最優設置水平的準確性往往得不到保證。為了提高生產設計中參數設計的可靠性，希望新的觀測數據能夠更新響應面模型，從而實現穩健參數設計的在線調整。對此，Zhou 等（2021）[10]提出了一種序貫SVR的在線穩健參數設計方法，他將先前得到的可控因子最優設置水平以及現階段觀測到的噪聲因子水平作為新的數據樣本更新響應面模型，進而在下一階段得到更可靠的參數設置水平。然而其研究仍存在一些不足，主要表現在精度和效率二者不能兼顧。

基于上述問題，本文將采用在線學習策略改進的OGP模型作為元模型構建響應面，在實現穩健參數設計實時調整的同時進一步提高了建模效率，進而提出了一種在線穩健參數設計方法（OGP-RPD）。

1 高斯過程（GP）模型

GP 模型是一種基于貝葉斯框架的非參數建模方法，它可以用未知的先驗信息來探索輸入變量和輸出響應之間的復雜關系[11]。

給定數據樣本集合D={(xi，yi)|i=1，…，n}，xi∈Rm，其中，xi為輸入，yi為輸出。假設輸入x與輸出y之間的函數關系為：

其中，ε是與f(x)獨立同分布的隨機噪聲，即高斯噪聲（ε～N(0，)）。f(x)是預測函數，其值為預測輸出，服從高斯過程分布：

其中：μ(x)表示高斯過程的均值函數，為簡化計算，通常取值為0；k(x，x')表示高斯過程的協方差函數，本文選取靜態平穩高斯過程中最常用的平方指數協方差函數。樣本xi和樣本xj之間的協方差函數定義如下：

其中：wd為模型的測度函數，用于衡量不同輸入維度的權重；為幅度參數，表示局部相關性的程度；為服從高斯分布的噪聲方差。以上參數組合成為協方差函數的超參數集合，表示為，Dim是輸入數據樣本的維數，δij為Kronecker 算子（δij=1 ifi=jand 0 otherwise）。

基于以上模型假設，對于給定訓練集X，可以得到測試樣本x*處的預測輸出和預測方差：

其中：K是K(X，X)的縮寫，表示訓練數據集的協方差矩陣；K*是K(X，x*)的縮寫，表示訓練集X與測試樣本x*之間協方差的向量；k*是k(x*，x*)的縮寫，表示測試樣本本身的協方差。

在確定協方差函數類型后，需要對其超參數集合Θ進行優化調整。本文根據ARD原理[12]，選擇極大似然法來優化超參數。模型超參數的對數邊際似然表示如下：

其中，Ky=K+I。將求解上述對數邊際似然的最大值轉換為其負對數邊際似然的最小值，即可得到：

接下來利用梯度信息求解，得到超參數Θ的導數：

其中，α=K-1y。一般情況下采用梯度下降的共軛梯度法搜索超參數的最優解。

2 在線高斯過程（OGP）模型

本文采用在線學習策略改進現有的GP 模型，使其具備處理流式數據的能力。隨著樣本數據的不斷增加，增量算法會面臨存儲空間飽和、歷史樣本價值降低等問題。尤其是當協方差矩陣K的規模不斷增大時，矩陣處理的計算成本呈指數級增長，甚至會陷入維數災難，因此，設置數據集的容量上限是有必要的。當樣本數量增加至數據集容量上限時，需要在新增樣本的同時刪去一個對模型影響程度最小的舊數據，也就是減量算法。這樣一來，在增量算法和減量算法的共同作用下，可以始終保持訓練數據的規模不變，能夠在保證合理的計算復雜度的基礎上實現GP模型的在線更新以及對新樣本的預測。在線高斯過程（Online Gaussian Process,OGP）模型的基本原理如圖1 所示。

圖1 OGP模型原理

從圖1可以看出，OGP模型主要需要解決兩個關鍵問題：（1）篩選出舊樣本；（2）更新模型與超參數。

針對第一個問題，由于需要從訓練數據集中刪除一個對模型影響程度最小的樣本，因此引入信息論中信息增益的概念（數據樣本間的信息增益反映影響模型的程度）。利用式（6）中的對數邊際似然來表示數據樣本間的信息增益。當從含有n個樣本的訓練數據集中刪除一個數據樣本時，需要計算n-1 個數據樣本的對數邊際似然。對數邊際似然越高意味著n-1 個數據樣本的信息增益越高，同時也就意味著剩余的那個數據對模型的影響越小。這樣一來，數據樣本的信息增益可以通過式（7）來計算。此外，本文引入一種遺忘機制來調整減量規則，從時間序列的角度分析，越晚加入訓練數據集的數據對模型的影響程度越大，所以刪除較早加入的數據是合理的。基于上述思想，本文給出減量算法的準則：

其中，λ為遺忘因子，n為當前數據樣本的序列號，c為當前加入訓練數據集中的數據樣本的序列號。利用式（9）計算得到舊數據的信息增益，并從中找到信息增益最小的舊數據樣本，將其從訓練數據集中刪除。

針對第二個問題，在線算法基于前一階段的計算結果進行更新迭代，以最小的計算代價實現模型的更新。所以無論是增加樣本還是刪除樣本，更新GP 模型實質上體現在協方差矩陣的更新。本文的在線迭代策略是一種基于貝葉斯方法的近似技術[13,14]。

根據GP 模型的定義，當加入新樣本xn+1時，均值μn與協方差函數kn(x，x')增量更新的近似值可以由式（10）迭代計算：

其中，系數qn+1和rn+1分別是對數似然的一、二階導數：

其中，Z=E{P(yn+1|fn+1(x))}n是已知的完全數據似然函數。

展開式（11）中的更新遞歸步驟得到GP模型后驗概率的近似值：

其中：kn+1=kxn+1；en+1是n+1 階單位向量；Tn+1和Un+1是增加樣本后擴展矩陣的n+1 維算子，是通過在向量的末尾和矩陣的最后一行和一列添0得到的。

利用上述遞歸方式可以實現增量過程的迭代更新，系數規模也將隨著新數據樣本的到來不斷增加。當數據樣本增加至訓練數據集容量上限時，需要進行減量操作，系數更新如下：

對于模型的超參數來說，本文沿用前文提到的GP 模型的超參數尋優方法，利用式（7）求取最大似然，之后采用共軛梯度法尋優超參數。

以上闡述了OGP模型的兩個關鍵性問題。由此根據圖1所示的OGP模型原理圖總結OGP模型具體步驟如下：

步驟1：加入新數據樣本。

步驟2：判斷訓練數據集是否達到規模上限。若未達到則轉步驟4；若已達到則轉步驟3。

步驟3：計算數據樣本iGain，刪除信息增益最小的舊數據樣本。

步驟4：更新OGP模型。

步驟5：更新OGP模型超參數。

重復步驟1至步驟5，直至沒有新樣本用于更新模型。

3 基于OGP模型的在線RPD(OGP-RPD)方法

在響應面方法（Response Surface Methodology,RSM）的理論指導下，發展出兩種RPD建模策略：一種是雙響應模型方法，分別對響應值的均值和方差建模；另一種是單響應模型方法，控制變量和噪聲變量共同作用于同一模型。其中方差響應面是由均值模型來確定的[15]。本文擬采用單響應模型方法構建優化模型。

構建優化模型的方法同樣有很多，大致思想為在保持均值的同時最小化方差或者在保持方差小于一定閾值的同時最小化均值和目標值的偏差。

不同于以上兩種思想，Lin 和Tu（1995）[16]提出了一個更簡單的解決方案：

該方法不是圍繞特定響應（均值或方差）進行優化，而是圍繞目標值進行優化。本文將式（15）作為RPD優化模型的目標函數。其中，y^ 表示響應面模型的預測輸出，T表示相應質量特征的目標值。

在OGP-RPD 方法中，將OGP 模型用于構建響應面，其輸入是可控因子和噪聲因子的組合(xc，xe)，輸出則是與生產過程/產品質量相關的屬性值y(xc，xe)，其理想目標值為T。在線RPD的思想為利用上一階段獲得的可控因子設置水平，將由當前階段觀察到的噪聲因子以及對應的響應值組成的數據樣本{(xc，xe)，y(xc，xe)}作為新樣本更新RPD 過程，即通過新樣本更新OGP 模型進而得到下一階段的可控因子設置水平。按照上述規則不斷對生產過程/產品質量進行在線調整，直至獲得可控因子的最優設置水平。下面詳細闡述OGP-RPD方法的具體實現步驟：

步驟1：選擇初始數據樣本。利用拉丁超立方抽樣（Latin Hypercube Sampling，LHS），同時基于maxmin 準則抽取初始樣本點Sn={ }s1，s2，…，sn，每個樣本輸入均包含可控因子和噪聲因子兩個部分，即表示為si=(xci，xei)，樣本輸入對應的響應輸出為y(·) 。

步驟2：構建響應面模型。利用OGP模型擬合抽取的初始數據樣本，構建響應面模型，探索樣本輸入與輸出之間的函數關系

步驟3：構建噪聲樣本分布簇。采用蒙特卡洛方法模擬高斯噪聲因子的分布，隨機生成ne個噪聲樣本。

步驟4：求解可控因子最優設置水平。利用優化策略式（15）中最小化偏差的思想構建RPD 優化模型，求解可獲得當前階段的可控因子最優設置水平。優化模型的目標函數如下：

其中，預測方差參見式（5），可得：

步驟6：更新OGP 模型，重復RPD 流程。用步驟5 中獲得的新樣本更新響應面模型（OGP），重復步驟3至步驟5，在線更新RPD 流程，直至得到滿足優化終止條件的可控因子最優設置水平，算法終止，在線RPD結束。

在步驟6 中，主要從以下兩個方面衡量RPD 是否有效：一方面，當獲得的參數最優設置水平與真實值（產品或過程潛在的真實參數最優設置，在生產實踐中往往可以根據大量實驗獲得）之間的偏差很小時，可認為已經達到了好的RPD 效果，算法可以終止；另一方面，當獲得的參數最優設置水平維持在一定水平且相對波動較小時，算法也可以終止。

在生產實踐中，當工程師們獲得推薦的參數最優設置水平時，并非直接將其應用于產品生產或過程控制中，而是需要在該參數水平下進行數次驗證試驗。由于噪聲因子的影響，可能會得出多個不同的輸出響應值，這時可以計算這些響應的均值和方差，進而計算當前參數水平下的均方誤差（MSE）。可以將MSE作為一種有效的評估指標，如果算法推薦的參數水平不夠理想，那么可以選擇在整個在線RPD過程中MSE值最小的參數設計水平作為參數最優設置水平。

4 仿真與實證算例

本文選取了一個仿真算例和一個實證算例，均屬于研究復雜產品質量設計問題的經典案例。通過與Zhou 等（2021）[10]提出的在線穩健參數設計方法（SSVR-RPD）進行對比，從而驗證本文所提方法（OGP-RPD）在穩健參數設計問題上的有效性和優越性。

4.1 仿真算例

本實驗采用Branin函數[17]作為仿真測試函數，Branin函數被廣泛應用于計算機實驗。其函數定義在χ=[-5，10]×[0，15]范圍內，表示如下：

假設上述定義式中變量x1是可控因子，x2是噪聲因子且服從高斯分布（x2～N(0.5，1)），y(·) 表示在可控因子和噪聲因子共同作用下的真實系統的響應值。

為了驗證所提方法（OGP-RPD）是否能夠有效地找到可控因子的最優設置水平，需要先知道仿真系統可控因子潛在的真實最優水平。本文將式（15）中的響應面替換為本實驗中的Branin函數，然后通過求解優化問題即可得到仿真系統潛在的真實可控因子最優設置水平。式（15）中T表示函數系統的響應目標值，在本次實驗中是根據Branin函數的真實分布情況來設置的。然而從數學的角度精準計算是比較困難的，尤其在不可積的情況下更是難以完成。因此本文統一采用蒙特卡洛近似方法進行積分近似求解，為此隨機產生1×105個高斯白噪聲樣本。在得到的近似解后，利用MATLAB工具箱中的fmincon函數求解式（15）表示的最優化問題，得出Branin函數仿真系統的可控因子最優設置水平。可控因子最優設置水平被確定后，即可以根據生成的1×105個噪聲樣本計算對應響應值y(·) 及其波動情況，將其記為MeanY和VarY。由于本文采用的是近似求解，而且每次生成的噪聲因子樣本是不同的，因此為了降低不確定性，重復上述優化求解過程10次并求取平均值，最終得到Branin函數仿真系統的數據結果（如表1所示）。從表1中可以看出，Branin函數仿真系統的可控因子潛在真實最優設置水平為5.1808，其MeanY的值穩定在57.8268左右，且波動較小。

表1 仿真算例可控因子真實最優設置水平

為了驗證OGP-RPD是否能夠準確有效地找到可控因子的真實最優設置水平，基于高斯過程的在線穩健參數設計（OGP-RPD）的過程如下：

第一步：確定初始數據集。在實際生產中，工程師們根據生產過程中生成的樣本構建模型。而在本文的仿真實驗中，生產過程由式（19）表示。本文采用LHS方法，遵循maxmin準則隨機從系統中抽取20個樣本點，并獲得其對應響應值。每個數據樣本由影響因子（可控因子與噪聲因子的組合）及其對應響應值組成，記為{(x1，x2)，y(·) }。

第二步：構建響應面模型。將上述抽取的20個數據樣本作為初始訓練數據集，采用OGP模型構建響應面。

第三步：求解優化策略。用得到的OGP模型替換優化策略（見式（15））中的響應面，當前狀態下的優化策略變成式（16），求解式（16）中的優化問題需要計算期望值，考慮到所建模型的函數可能不可積，因此依舊采用蒙特卡洛方法來近似求解，為此隨機生成了ne個噪聲樣本（ne=1×104），其服從的分布為N(0，0.5)。然后可以通過MATLAB工具箱中的fmincon 函數進行求解，進而可以得到當前的可控因子最優設置水平，如表2所示。

表2 OGP-RPD仿真算例實驗結果

第四步：判斷當前可控因子的最優設置水平是否達到優化終止條件。通過當前可控因子的最優設置水平處ne個噪聲樣本可以計算MeanY和VarY，結果如表2 所示。可以看到，第一次參數設計得到的可控因子最優設置水平為-3.0343，對應的MeanY值為24.6524。很顯然與表1中的真實數據相比，二者都相差甚遠。得到的可控因子設置水平不理想，其中一個重要原因是噪聲因子的影響導致構建的OGP 模型不夠準確，這時就需要增加新樣本更新OGP模型。

第五步：確定新樣本并更新響應面模型。根據OGP-RPD 方法原理，在當前階段獲得的可控因子最優設置水平下，利用式（17）找出影響波動最大的噪聲因子并將其作為新樣本噪聲因子的觀測值。同時可以根據式（19）求得當前的輸出響應值。由此得到新樣本，將其加入訓練數據集更新響應面模型y^OGP。

重復第三步到第五步的實驗流程，直至找到滿足優化終止條件的可控因子最優設置水平。從表2 展示的實驗結果來看，當更新迭代到第9 次時，可控因子的最優設置水平為5.1768，其對應的MeanY值為57.8564，這與表1中Branin 函數仿真系統的可控因子真實最優設置水平極為接近，且MSE值也相對較低。由此可以證明OGP-RPD方法能夠有效找到理想的可控因子最優設置水平。

為了證明本文所提OGP-RPD 方法的先進性，采用目前在線RPD 研究中表現較為優異的SSVR-RPD 方法[7]進行對比分析，實驗結果如表3所示。可以看到，SSVR-RPD方法更新迭代至第18到20次時才得到了較為理想的可控因子最優設置水平。由此可以得出結論：OGP-RPD 方法可以更精準地得到可控因子的最優設置水平，且效率遠高于SSVR-RPD方法。

表3 SSVR-RPD仿真算例實驗結果

4.2 實證算例

某個半導體零件加工過程實驗[18]共包含5 個影響因子，其中包括2 個可控因子(x1，x2)和3 個噪聲因子(z1，z2，z3) 。此實驗以5 個影響因子的標準CCD 設計為基礎，刪去了與3個噪聲因子有關的軸向試驗點的中心復合設計，共進行23次實驗。最后得出在可控因子、三個噪聲因子與響應輸出之間具有二階模型，擬合的響應模型為：

在生產過程中，質量特性的目標值為T=30。噪聲因子服從均值為0、方差為1 的正態分布，即z1～N(0，1),z2～N(0，1) 和z3～N(0，1) 。

與之前仿真算例的流程類似，利用式（20）所表示的真實響應模型確定半導體零件加工過程中可控因子的真實最優設置水平，相關結果見表4。

表4 實證研究算例可控因子真實最優設置水平

就OGP-RPD方法的有效性而言，整個在線RPD過程與仿真實驗程序一致。實驗的最終結果如表5所示。

表5 OGP-RPD實證研究算例實驗結果

當更新迭代至第10次時，可以看到MeanY已趨于平穩且非常接近表4中真實的響應輸出值，此外MSE值也達到了相對較低的水平，因此能夠判定達到了優化終止條件。下頁表6展示了SSVR-RPD方法的實驗結果，可以看到更新迭代至第20 至25 次時，MeanY值趨于穩定，且MSE 值已達到全局最低水平。通過對比兩個方法的實驗結果，很顯然OGP-RPD 方法尋找到的可控因子最優水平更理想，輸出響應更接近真實響應目標值，且效率更高。

表6 SSVR-RPD實證研究算例實驗結果

綜合上述兩組實驗結果可知，OGP-RPD 方法能夠有效找到可控因子的最優設置水平，且較好地兼顧了輸出響應的最優性和穩健性。

5 結論

傳統的RPD方法是離線設計，受環境因素的影響，通過一次性建模獲得的可控因子設置水平可能并不理想。對此本文提出了一種基于GP 模型的在線RPD 方法（OGP-RPD）。本文采用滑動窗口的在線策略改進傳統GP模型，使其具備處理流式數據的能力，將其稱為OGP模型。利用OGP 模型構建響應面，并基于響應面理論中的單響應優化策略提出了一套完整的在線RPD流程。通過仿真算例與實證研究，將OGP-RPD 方法與現有的SSVR-RPD 方法進行比較。實驗結果表明，OGP-RPD 方法在尋找可控因子的最優設置水平方面準確度更高，能夠更好地兼顧輸出響應的最優性和穩健性，且效率更高。由此可以得出結論，OGP-RPD 方法可以有效解決RPD 離線設計帶來的問題，能夠更好地應用于產品或過程的在線質量設計，對于提高產品質量具有重要理論和實踐意義。