汪建均， 屠雅楠， 馬義中

(南京理工大學(xué)經(jīng)濟(jì)管理學(xué)院， 南京 210094)

0 引 言

產(chǎn)品質(zhì)量不僅是企業(yè)的生命線，更是在全球市場(chǎng)上贏得顧客的關(guān)鍵，因此連續(xù)性質(zhì)量改進(jìn)已經(jīng)成為世界各大企業(yè)永恒追求的目標(biāo).現(xiàn)代質(zhì)量工程的主流觀點(diǎn)認(rèn)為：波動(dòng)是導(dǎo)致產(chǎn)品產(chǎn)生質(zhì)量問(wèn)題的根本原因.因此，為了改進(jìn)和提高產(chǎn)品質(zhì)量，必須最大限度地減小和控制圍繞設(shè)計(jì)目標(biāo)值的波動(dòng).如何減小或控制產(chǎn)品實(shí)現(xiàn)過(guò)程中的波動(dòng)，已經(jīng)成為連續(xù)質(zhì)量改進(jìn)活動(dòng)的核心內(nèi)容[1].質(zhì)量設(shè)計(jì)應(yīng)用于產(chǎn)品形成的早期階段，能夠從源頭上查找導(dǎo)致產(chǎn)品產(chǎn)生缺陷的原因，因此通常認(rèn)為在產(chǎn)品設(shè)計(jì)階段實(shí)施質(zhì)量改進(jìn)活動(dòng)是最具有生產(chǎn)效率和成本效益的生產(chǎn)管理方法[2].近年來(lái)，產(chǎn)品設(shè)計(jì)趨于復(fù)雜化，顧客需求層次也呈現(xiàn)多樣化，在產(chǎn)品的質(zhì)量設(shè)計(jì)過(guò)程中往往需要同時(shí)考慮多個(gè)質(zhì)量特性[3]，因此考慮多質(zhì)量特性的優(yōu)化設(shè)計(jì)在連續(xù)質(zhì)量改進(jìn)活動(dòng)中顯示出越來(lái)越重要的作用與地位[4].

響應(yīng)曲面方法涉及統(tǒng)計(jì)試驗(yàn)與數(shù)據(jù)分析、模型構(gòu)建方法以及數(shù)值優(yōu)化技術(shù)等一系列的連續(xù)質(zhì)量改進(jìn)方法或技術(shù)[5].響應(yīng)曲面方法(response surface methodology, RSM)已經(jīng)廣泛地應(yīng)用到產(chǎn)品或過(guò)程的連續(xù)質(zhì)量改進(jìn)活動(dòng)中，成為了工業(yè)試驗(yàn)的核心[6].因此，基于RSM的多質(zhì)量特性優(yōu)化設(shè)計(jì)方法已引起眾多研究者的關(guān)注，并取得了一系列新的研究成果[4].通常來(lái)說(shuō)，以質(zhì)量設(shè)計(jì)為基礎(chǔ)的連續(xù)質(zhì)量改進(jìn)活動(dòng)主要由變量篩選、響應(yīng)建模和參數(shù)優(yōu)化三個(gè)階段所構(gòu)成[7].在變量篩選階段，研究者旨在通過(guò)科學(xué)的試驗(yàn)設(shè)計(jì)，篩選出響應(yīng)曲面模型的顯著性變量，從而為后續(xù)的響應(yīng)曲面模型構(gòu)建奠定堅(jiān)實(shí)的基礎(chǔ).著名統(tǒng)計(jì)學(xué)家Wu和Hamada[8]在其經(jīng)典著作《試驗(yàn)計(jì)劃、分析與參數(shù)設(shè)計(jì)優(yōu)化》中指出“因子效應(yīng)原則(效應(yīng)稀疏原則、效應(yīng)排序原則以及效應(yīng)遺傳原則)能夠評(píng)判因子設(shè)計(jì)結(jié)構(gòu)的合理性，是否考慮了因子效應(yīng)原則是區(qū)分一般回歸問(wèn)題與試驗(yàn)設(shè)計(jì)問(wèn)題的關(guān)鍵，因此考慮因子效應(yīng)原則的變量篩選對(duì)提高產(chǎn)品或過(guò)程的設(shè)計(jì)質(zhì)量具有重要的意義”.在響應(yīng)建模階段，研究者旨在通過(guò)響應(yīng)曲面模型反映試驗(yàn)因子與多質(zhì)量特性之間的函數(shù)關(guān)系.然而，正如著名統(tǒng)計(jì)學(xué)家Box等[9]所言 “所有的模型都是錯(cuò)誤的，但是有些模型是有用的”.傳統(tǒng)的響應(yīng)曲面建模方法往往忽略模型參數(shù)[10]或模型結(jié)構(gòu)的不確定性[11]以及響應(yīng)預(yù)測(cè)值的波動(dòng)問(wèn)題[12]，不僅可能導(dǎo)致得到不科學(xué)的研究結(jié)果，而且將極大地降低了試驗(yàn)結(jié)果的可重復(fù)性.在參數(shù)優(yōu)化階段，研究者旨在通過(guò)構(gòu)建穩(wěn)健可靠的目標(biāo)函數(shù)和高效的優(yōu)化算法，尋求一組最佳的參數(shù)設(shè)計(jì)值，使得多個(gè)質(zhì)量特性同時(shí)滿足設(shè)定的目標(biāo)值.多優(yōu)化目標(biāo)之間的沖突、多響應(yīng)過(guò)程的穩(wěn)健性度量以及優(yōu)化結(jié)果的可靠性已成為該階段需要重點(diǎn)考慮的研究問(wèn)題.上述三個(gè)階段的質(zhì)量設(shè)計(jì)問(wèn)題往往相互疊加、相互影響，通常會(huì)造成響應(yīng)預(yù)測(cè)值與設(shè)計(jì)目標(biāo)值之間存在較大的波動(dòng).然而，以往的研究往往關(guān)注于其中某個(gè)階段的問(wèn)題，未能從連續(xù)質(zhì)量改進(jìn)的視角全面地系統(tǒng)考慮上述多質(zhì)量特性的質(zhì)量設(shè)計(jì)問(wèn)題.為此，本文將在SUR響應(yīng)曲面的貝葉斯建模與優(yōu)化框架下，首先，結(jié)合因子效應(yīng)原則篩選出顯著性的模型變量；然后，在此基礎(chǔ)上運(yùn)用貝葉斯抽樣技術(shù)構(gòu)建多變量過(guò)程能力指數(shù)；最后運(yùn)用混合遺傳算法對(duì)所構(gòu)建的貝葉斯多變量過(guò)程能力指數(shù)進(jìn)行參數(shù)優(yōu)化，從而獲得產(chǎn)品或過(guò)程的最佳參數(shù)設(shè)計(jì)值.

1 國(guó)內(nèi)外研究現(xiàn)狀及分析

近年來(lái)，隨著以質(zhì)量設(shè)計(jì)為基礎(chǔ)的連續(xù)質(zhì)量改進(jìn)活動(dòng)在企業(yè)中廣泛應(yīng)用和推廣，以響應(yīng)曲面建模為主的質(zhì)量設(shè)計(jì)問(wèn)題引起一些研究者的廣泛關(guān)注和重視.在響應(yīng)曲面模型變量的篩選階段，一些學(xué)者提出了一系列新穎的研究方法.例如，Beattie[13]針對(duì)超飽和試驗(yàn)設(shè)計(jì)，在隨機(jī)搜索的基礎(chǔ)上，結(jié)合貝葉斯因子方法提出了一種兩階段的因子篩選方法.Song等[14]針對(duì)因子數(shù)遠(yuǎn)大于試驗(yàn)次數(shù)的超高維試驗(yàn)設(shè)計(jì)，提出一種基于隨機(jī)近似蒙特卡洛(stochastic approximate Monte Carlo，SAMC)算法的兩階段貝葉斯變量選擇方法.該方法首先將高維數(shù)據(jù)切割成許多子集，從每個(gè)子集中篩選因子，然后將篩選出的因子進(jìn)行匯總，從中識(shí)別出顯著性因子.針對(duì)非正態(tài)響應(yīng)的部分因子試驗(yàn)，當(dāng)篩選試驗(yàn)所涉及的因子數(shù)目較大時(shí)，汪建均和馬義中[15]提出基于廣義線性模型(generalized linear models，GLM)的貝葉斯變量選擇與模型選擇方法.該方法通過(guò)經(jīng)驗(yàn)貝葉斯先驗(yàn)考慮了模型參數(shù)不確定性的影響，然后在廣義線性模型的線性預(yù)測(cè)器中對(duì)每個(gè)變量設(shè)置一個(gè)二元變量指示器，在此基礎(chǔ)上建立起變量指示器與模型指示器之間的轉(zhuǎn)換關(guān)系，然后利用變量指示器與模型指示器的后驗(yàn)概率來(lái)識(shí)別顯著性因子與選擇最佳模型.上述方法雖然都篩選出了模型的顯著性變量，但是在因子篩選或變量選擇過(guò)程中忽略了因子效應(yīng)原則的影響.因此，上述研究結(jié)果往往只是統(tǒng)計(jì)意義上的顯著，所獲得研究結(jié)果可能并不符合試驗(yàn)設(shè)計(jì)的因子效應(yīng)原則.為此，一些學(xué)者在進(jìn)行響應(yīng)曲面建模時(shí)考慮因子效應(yīng)原則以便篩選出更加符合實(shí)際的顯著性變量.例如，針對(duì)含有復(fù)雜別名的試驗(yàn)設(shè)計(jì)問(wèn)題，Chipman等[16]在一般回歸模型基礎(chǔ)上考慮因子效應(yīng)原則，并運(yùn)用貝葉斯變量隨機(jī)搜索方法來(lái)識(shí)別顯著性因子，提出了一種新的貝葉斯變量選擇方法.針對(duì)兩水平的部分因子試驗(yàn)設(shè)計(jì)問(wèn)題，Bergquist等[17]在響應(yīng)曲面建模過(guò)程中依次考慮效應(yīng)稀疏原則、效應(yīng)層次原則和效應(yīng)遺傳原則，提出了一種新的貝葉斯變量選擇方法.針對(duì)含有復(fù)雜別名結(jié)構(gòu)的非正態(tài)響應(yīng)部分因子試驗(yàn)，汪建均等[18]在廣義線性模型的框架下將因子效應(yīng)原則融入到了模型參數(shù)的先驗(yàn)信息之中，提出了一種多階段的貝葉斯變量選擇方法.Noguchi等[19]提出了一種改進(jìn)的貝葉斯LASSO (least absolute shrinkage and selection operator, LASSO)方法.該方法將效應(yīng)遺傳原則和效應(yīng)層次原則加入LASSO模型，并利用馬爾科夫蒙特卡羅(Markov chain Monte Carlo, MCMC)方法計(jì)算了各因子和模型結(jié)構(gòu)的后驗(yàn)概率.Taylorrodriguez等[20]提出考慮因子效應(yīng)遺傳原則的貝葉斯變量選擇方法.該方法為模型空間設(shè)定不同的先驗(yàn)，并利用MH算法搜索符合遺傳效應(yīng)原則的模型.針對(duì)單個(gè)質(zhì)量特性的穩(wěn)健參數(shù)設(shè)計(jì)問(wèn)題，歐陽(yáng)林寒等[21]考慮了模型參數(shù)以及模型結(jié)構(gòu)不確定性的影響，結(jié)合貝葉斯模型平均方法考慮了因子效應(yīng)原則，從而獲得更加符合實(shí)際、更為精確的響應(yīng)曲面模型.上述研究往往只涉及到單響應(yīng)情形下的變量篩選問(wèn)題，很少涉及到多響應(yīng)的變量篩選以及后續(xù)的響應(yīng)曲面建模與參數(shù)優(yōu)化問(wèn)題.

在響應(yīng)曲面模型構(gòu)建過(guò)程中，模型參數(shù)和模型結(jié)構(gòu)的不確定性以及響應(yīng)預(yù)測(cè)值的波動(dòng)問(wèn)題引起一些研究者廣泛的興趣和普遍的關(guān)注.針對(duì)模型參數(shù)和模型結(jié)構(gòu)不確定性，研究人員通常運(yùn)用置信區(qū)間或貝葉斯方法開展相關(guān)的研究.例如，針對(duì)多響應(yīng)優(yōu)化中模型的不確定性問(wèn)題，He等[11]通過(guò)擬合響應(yīng)曲面模型的置信水平考慮了模型不確定性對(duì)優(yōu)化結(jié)果的影響，從而構(gòu)建出了一種解決多響應(yīng)曲面優(yōu)化問(wèn)題的穩(wěn)健理想函數(shù).該方法由于考慮模型預(yù)測(cè)結(jié)果的不確定性，因此能夠在實(shí)踐中獲得比傳統(tǒng)理想函數(shù)更為可靠的優(yōu)化結(jié)果.針對(duì)多響應(yīng)優(yōu)化設(shè)計(jì)的問(wèn)題，Ouyang等[22]運(yùn)用多響應(yīng)的預(yù)測(cè)區(qū)間考慮了模型不確定性的影響，提出了一種新的質(zhì)量損失函數(shù).該方法不僅克服了傳統(tǒng)質(zhì)量損失函數(shù)忽視模型不確定性對(duì)研究結(jié)果的影響，而且所提的優(yōu)化策略能夠適應(yīng)不同的特殊情形(如不考慮模型不確定性的最佳策略以及考慮模型不確定性的穩(wěn)健性保守策略).He等[23]通過(guò)置信區(qū)間考慮模型的不確定性影響，并運(yùn)用穩(wěn)健的模糊規(guī)劃方法同時(shí)考慮了多響應(yīng)曲面優(yōu)化過(guò)程所涉及到的位置效應(yīng)、散度效應(yīng)以及模型不確定性等問(wèn)題，從而確保了研究結(jié)果的穩(wěn)健性.針對(duì)多響應(yīng)優(yōu)化問(wèn)題，Ouyang等[24]考慮模型參數(shù)的不確定性以及實(shí)施誤差(即參數(shù)設(shè)計(jì)值難以在實(shí)際操作中準(zhǔn)確無(wú)誤地實(shí)施)，提出一種新的質(zhì)量損失函數(shù).針對(duì)多響應(yīng)的優(yōu)化問(wèn)題，張旭濤和何楨[25]考慮考慮響應(yīng)間相關(guān)性和可控因子波動(dòng)的影響，提出了一種基于似無(wú)關(guān)回歸的多元穩(wěn)健損失函數(shù)方法.該方法不僅利用似不相關(guān)回歸模型考慮響應(yīng)間的相關(guān)性，而且利用給定點(diǎn)的梯度信息估計(jì)了可控因子波動(dòng)對(duì)過(guò)程穩(wěn)健性的影響.針對(duì)作用關(guān)系復(fù)雜且多極值的制造過(guò)程，崔慶安[26]考慮了模型結(jié)構(gòu)的不確定性對(duì)優(yōu)化結(jié)果的影響，運(yùn)用支持向量機(jī)作為復(fù)雜作用關(guān)系過(guò)程的近似模型，提出了基于支持向量聚類與序列二次規(guī)劃的參數(shù)全局性優(yōu)化方法.在多響應(yīng)優(yōu)化設(shè)計(jì)中，模型參數(shù)的不確定性以及生產(chǎn)過(guò)程的噪聲因子不可避免地會(huì)導(dǎo)致預(yù)測(cè)響應(yīng)值出現(xiàn)較大的波動(dòng).針對(duì)上述的問(wèn)題，汪建均和屠雅楠[27]結(jié)合貝葉斯抽樣技術(shù)、帕累托優(yōu)化策略以及灰色關(guān)聯(lián)分析方法提出了一種多響應(yīng)優(yōu)化設(shè)計(jì)方法.Chapman等[28]提出一種利用置信區(qū)間度量預(yù)測(cè)響應(yīng)值波動(dòng)的方法，該方法利用隨機(jī)抽樣法得到響應(yīng)預(yù)測(cè)值的置信區(qū)間，將該區(qū)間端點(diǎn)的最優(yōu)值作為響應(yīng)預(yù)測(cè)值，同時(shí)利用圖形顯示響應(yīng)預(yù)測(cè)值波動(dòng)對(duì)優(yōu)化結(jié)果的影響.由于貝葉斯方法能夠在建模過(guò)程中有效地利用相關(guān)的先驗(yàn)信息和專家知識(shí)，而且能夠在后續(xù)的建模過(guò)程中對(duì)以往的先驗(yàn)信息進(jìn)行不斷地更新和修正，因此貝葉斯方法在不確定性質(zhì)量設(shè)計(jì)(模型參數(shù)、模型結(jié)構(gòu)以及響應(yīng)分布等不確定性)的建模過(guò)程中也得到了廣泛地應(yīng)用.例如，針對(duì)多響應(yīng)優(yōu)化問(wèn)題，Peterson[29]提出了一種貝葉斯的后驗(yàn)預(yù)測(cè)方法.該方法在貝葉斯建模與優(yōu)化過(guò)程中考慮了試驗(yàn)數(shù)據(jù)之間的相關(guān)結(jié)構(gòu)、過(guò)程分布的變化以及模型參數(shù)的不確定性，然后利用多響應(yīng)的后驗(yàn)概率分布計(jì)算響應(yīng)預(yù)測(cè)值滿足所設(shè)定質(zhì)量條件的概率.Miro-Quesada等[30]擴(kuò)展了Peterson[29]的研究工作，在多響應(yīng)的后驗(yàn)概率密度函數(shù)中考慮了噪聲因子的影響.此外，Peterson等[31]考察了多響應(yīng)之間存在不同協(xié)方差的問(wèn)題，結(jié)合多變量似不相關(guān)回歸模型進(jìn)一步改進(jìn)以往的研究工作.上述方法過(guò)分關(guān)注響應(yīng)預(yù)測(cè)值落在規(guī)格限內(nèi)的概率，卻忽視了多元過(guò)程的穩(wěn)健性[32].事實(shí)上，在很多情況下僅僅考慮優(yōu)化結(jié)果的可靠性所獲得的參數(shù)設(shè)計(jì)值，其可靠性結(jié)果往往令人滿意，然而其質(zhì)量損失卻相當(dāng)大[7].為此，汪建均等[4]在貝葉斯多元回歸模型的統(tǒng)一框架下，結(jié)合多元質(zhì)量損失函數(shù)與貝葉斯后驗(yàn)概率方法提出一種新的多響應(yīng)優(yōu)化方法.該方法能夠在優(yōu)化過(guò)程中較好地兼顧多元過(guò)程的穩(wěn)健性和優(yōu)化結(jié)果的可靠性，從而為實(shí)現(xiàn)多響應(yīng)穩(wěn)健參數(shù)設(shè)計(jì)提供了各方面(如多元過(guò)程的穩(wěn)健性、優(yōu)化結(jié)果的可靠性)均較滿意的優(yōu)化結(jié)果.此外，Vanli等[33]在考慮模型參數(shù)不確定性的穩(wěn)健參數(shù)設(shè)計(jì)時(shí)，融入了在線可觀測(cè)的噪聲因子，提出了一種自適應(yīng)的貝葉斯穩(wěn)健參數(shù)設(shè)計(jì)方法.該方法旨在通過(guò)在線觀察噪聲因子的水平，從而不斷調(diào)節(jié)最佳的輸入?yún)?shù)水平.

此外，在參數(shù)優(yōu)化過(guò)程中，優(yōu)化指標(biāo)的選擇對(duì)研究結(jié)果的影響至關(guān)重要.針對(duì)多響應(yīng)曲面優(yōu)化設(shè)計(jì)，研究者通常采用維度縮減技術(shù)構(gòu)建一個(gè)簡(jiǎn)化的綜合性能指標(biāo)，從而將多響應(yīng)曲面優(yōu)化問(wèn)題轉(zhuǎn)換為單目標(biāo)優(yōu)化問(wèn)題.在多響應(yīng)曲面優(yōu)化設(shè)計(jì)中，常見的綜合性能指標(biāo)主要有：理想函數(shù)[11, 34, 35]、質(zhì)量損失函數(shù)[10, 36, 37]、馬氏距離函數(shù)、后驗(yàn)概率函數(shù)[7, 29, 31]以及多變量過(guò)程能力指數(shù)[38-40].多變量過(guò)程能力指數(shù)不僅考慮過(guò)程是否滿足規(guī)格要求，而且還考慮過(guò)程與目標(biāo)值之間的偏差與波動(dòng)，因此在多響應(yīng)優(yōu)化設(shè)計(jì)中選擇多變量過(guò)程能力作為綜合性能度量指標(biāo)能夠較好地權(quán)衡產(chǎn)品的符合性與質(zhì)量損失二者之間的關(guān)系.因此，一些研究開始嘗試運(yùn)用多變量過(guò)程能力指數(shù)來(lái)優(yōu)化多響應(yīng)優(yōu)化設(shè)計(jì)問(wèn)題.例如，Amir等[38]利用多變量過(guò)程能力指數(shù)解決非正態(tài)響應(yīng)的優(yōu)化問(wèn)題.Awad等[39]將多變量過(guò)程能力指數(shù)作為目標(biāo)函數(shù)，輸入變量的限制條件作為約束，提出一種基于多變量過(guò)程能力指數(shù)的多響應(yīng)優(yōu)化設(shè)計(jì)方法.Ouyang等[24]利用多變量過(guò)程能力指數(shù)來(lái)評(píng)價(jià)不同方法的優(yōu)化效果.

本文擬在SUR模型的框架下，結(jié)合貝葉斯建模與優(yōu)化方法，從變量篩選、響應(yīng)曲面建模以及參數(shù)優(yōu)化三個(gè)階段開展連續(xù)質(zhì)量改進(jìn).首先，為SUR模型中的各個(gè)變量設(shè)置二元變量指示器，在此基礎(chǔ)上構(gòu)建符合因子效應(yīng)原則的混合二元變量指示器；然后利用MCMC方法計(jì)算混合二元變量指示器和模型結(jié)構(gòu)的后驗(yàn)概率，從而確定模型中的顯著性變量；最后，優(yōu)化多變量過(guò)程能力指數(shù)獲得最佳的參數(shù)設(shè)計(jì)值.

2 似不相關(guān)回歸模型的貝葉斯分析

在多響應(yīng)優(yōu)化設(shè)計(jì)中，假設(shè)有m個(gè)響應(yīng)，各個(gè)響應(yīng)包含n個(gè)觀測(cè)值和q個(gè)因子效應(yīng)，則一般多響應(yīng)回歸模型(SMR)為

(1)

其中Yj為n×1的響應(yīng)矩陣；βj為q×1的回歸系數(shù)矩陣；Xj為n×q的因子效應(yīng)矩陣；ej為n×1的隨機(jī)誤差矩陣，服從均值向量為0、方差-協(xié)方差矩陣為H的正態(tài)分布,而且H為對(duì)角陣.式(1)中隨機(jī)誤差項(xiàng)的方差-協(xié)方差矩陣為對(duì)角陣，其非對(duì)角線元素均為0，表明該式中的m個(gè)響應(yīng)是相互獨(dú)立的.但是在實(shí)際生產(chǎn)中，需要解決大量相關(guān)多響應(yīng)的優(yōu)化問(wèn)題.Zellner[41]在1962年首次提出SUR模型，該模型與SMR模型最大的區(qū)別在于：該模型中的每個(gè)響應(yīng)包含的隨機(jī)誤差項(xiàng)是相關(guān)的，即響應(yīng)之間是相互聯(lián)系的.假設(shè)需要考慮m個(gè)響應(yīng)模型，每個(gè)響應(yīng)有n個(gè)觀測(cè)值和mj個(gè)解釋變量，則SUR模型的結(jié)構(gòu)如下

(2)

其中Yj為n×1的響應(yīng)矩陣；Xj為n×mj的因子效應(yīng)矩陣；βj為mj×1的回歸系數(shù)矩陣；εj為n×1的隨機(jī)誤差矩陣.由于式(2)中的隨機(jī)誤差項(xiàng)εj不獨(dú)立，因此滿足以下性質(zhì)

(3)

還可將式(2)改寫為矩陣形式

Y=Xβ+ε,ε～N(0,Σ?I)

(4)

根據(jù)式(2)～式(4)，可計(jì)算出關(guān)于Y的似然函數(shù)為

(5)

其中|Σ|表示矩陣Σ的行列式；tr表示矩陣的跡(主對(duì)角線元素之和)；R為m×m的矩陣，其元素為rij=(Yi-Xiβi)T(Yj-Xjβj).在式(5)的基礎(chǔ)上，可求得參數(shù)β和Σ的極大似然估計(jì)，比如Zellner[41]采用迭代的SUR方法獲得了β和Σ的參數(shù)估計(jì)值.

由于貝葉斯方法能夠在建模過(guò)程中有效地考慮先驗(yàn)信息或?qū)＜抑R(shí)，因此被廣泛地應(yīng)用到復(fù)雜模型的參數(shù)估計(jì)與統(tǒng)計(jì)推斷之中.近年來(lái)，一些學(xué)者構(gòu)建了SUR模型的貝葉斯推斷方法.貝葉斯的先驗(yàn)分布大體上可分為兩類：無(wú)信息先驗(yàn)和共軛先驗(yàn)，因此SUR模型的貝葉斯參數(shù)估計(jì)與統(tǒng)計(jì)推斷通常也是從上述兩個(gè)角度展開分析.

第一類是在對(duì)參數(shù)信息缺乏了解的情況下，利用無(wú)信息先驗(yàn)作為參數(shù)的先驗(yàn)信息，此處傾向于選擇Jeffrey無(wú)信息先驗(yàn)[42].假設(shè)參數(shù)β和Σ是相互獨(dú)立的，則參數(shù)β和Σ的聯(lián)合先驗(yàn)分布為

(6)

根據(jù)式(5)和式(6)，可以求得參數(shù)的聯(lián)合后驗(yàn)密度

(7)

在參數(shù)β和Σ聯(lián)合后驗(yàn)密度的基礎(chǔ)上，可計(jì)算參數(shù)β和Σ各自的條件后驗(yàn)分布為

(8)

h1(Σ|Y,X,β)∝IW(R,n)

(9)

(10)

(11)

第二類是在根據(jù)過(guò)去的經(jīng)驗(yàn)知識(shí)，假設(shè)參數(shù)的先驗(yàn)分布與后驗(yàn)分布具有相同的形式，則參數(shù)β和Σ的聯(lián)合先驗(yàn)分布為

π2(β,Σ)=π2(β)π2(Σ)

(12)

其中

π2(β)=N(β0,A-1)

(13)

π2(Σ)=IW(Λ0,ν0)

(14)

根據(jù)式(5)和參數(shù)β和Σ的先驗(yàn)分布，可推斷得參數(shù)β和Σ的條件后驗(yàn)分布為

(15)

h2(Σ|Y,X,β)∝IW(Λ0+R,n+ν0)

(16)

其中

(17)

(18)

3 結(jié)合SUR模型與因子效應(yīng)原則的連續(xù)質(zhì)量改進(jìn)

3.1 考慮因子效應(yīng)原則的SUR模型

Wu和Hamada[8]指出：因子效應(yīng)原則通常包括效應(yīng)稀疏原則、效應(yīng)排序原則以及效應(yīng)遺傳原則，主要用于選擇模型的顯著性變量.效應(yīng)稀疏原則指模型包含的重要性因子是少數(shù)的，可稱之為試驗(yàn)設(shè)計(jì)的帕累托原則.效應(yīng)排序原則認(rèn)為：同階效應(yīng)具有相同的重要性；低階效應(yīng)與高階效應(yīng)相比，低階效應(yīng)更可能是重要性因子.因此，當(dāng)試驗(yàn)次數(shù)有限時(shí)，應(yīng)當(dāng)優(yōu)先估計(jì)低階效應(yīng).效應(yīng)遺傳原則指：在考慮交互項(xiàng)時(shí)，如果該交互項(xiàng)是顯著的，那么它的親本因子中應(yīng)該至少有一個(gè)是顯著的，即不可能出現(xiàn)沒(méi)有親本因子的交互項(xiàng).若一個(gè)模型僅有交互項(xiàng)而未包含其親本因子，則該模型違背了效應(yīng)遺傳原則，應(yīng)該將其剔除.在篩選試驗(yàn)中，需要篩選出能夠解釋模型的顯著性變量，同時(shí)需要剔除那些擬合效果很好但違反因子效應(yīng)原則的模型.如果將這部分模型用于后續(xù)的質(zhì)量改進(jìn)階段，可能會(huì)導(dǎo)致預(yù)測(cè)模型更加偏離真實(shí)模型，從而影響試驗(yàn)最終的優(yōu)化結(jié)果.

在以往關(guān)于變量選擇的研究中，大都沒(méi)有考慮因子效應(yīng)原則，通常都是計(jì)算各個(gè)變量指示器γyj={γ1,γ2,…}j=1,2,…,m的后驗(yàn)概率來(lái)判斷該因子是否顯著.本文借鑒文獻(xiàn)的經(jīng)驗(yàn)[19]，通過(guò)構(gòu)建混合二元變量指示器來(lái)進(jìn)行因子篩選.假設(shè)共有m個(gè)響應(yīng)，p個(gè)主效應(yīng)，并考慮其二階交互項(xiàng)和平方項(xiàng)，則該混合二元變量指示器如下

(19)

其中γyj表示第j個(gè)響應(yīng)的一組混合二元變量指示器向量；γp表示第p個(gè)主效應(yīng)的變量指示器；γ(p-1)p表示第p-1個(gè)主效應(yīng)與第p個(gè)主效應(yīng)之間的交互項(xiàng)對(duì)應(yīng)的變量指示器；γp2表示第p個(gè)主效應(yīng)的平方項(xiàng)對(duì)應(yīng)的變量指示器.

式(19)中的模型結(jié)構(gòu)表明該設(shè)定符合效應(yīng)排序原則和效應(yīng)遺傳原則.以交互項(xiàng)為例，該模型利用三個(gè)變量指示器的乘積γ(p-1)p×γp-1×γp來(lái)判斷該交互效應(yīng)是否顯著.如果該交互項(xiàng)的兩個(gè)親本因子中至少有一項(xiàng)為非顯著性主效應(yīng)，即γp-1和γp中至少有一個(gè)值等于0.那么該交互項(xiàng)的變量指示器γ(p-1)p等于1，否則γ(p-1)p×γp-1×γp的值仍然等于0，因此該交互效應(yīng)將判定為非顯著性因子，符合效應(yīng)遺傳原則的假設(shè).此外，增加了γp-1和γp的約束后，γ(p-1)p×γp-1×γp=1的可能性要遠(yuǎn)遠(yuǎn)小于兩個(gè)親本因子的變量指示器γp-1=1和γp=1的可能性，因此低階主效應(yīng)為顯著性因子的概率將遠(yuǎn)遠(yuǎn)高于高階交互效應(yīng)，這與效應(yīng)排序原則的思想保持一致.

為了在SUR模型中更好地融入因子效應(yīng)原則，在一般SUR模型式(4)的基礎(chǔ)上，引入式(19)中的混合二元變量指示器，從而構(gòu)建出一種結(jié)合SUR模型與因子效應(yīng)原則的新模型

(20)

即

(21)

本文擬在上述新模型的框架下構(gòu)建考慮因子效應(yīng)原則的貝葉斯變量選擇方法.借鑒一般SUR模型的第二類貝葉斯分析方法，擬對(duì)模型參數(shù)采用一種經(jīng)驗(yàn)貝葉斯先驗(yàn).首先，式(20)中參數(shù)β的先驗(yàn)形式采用多變量正態(tài)分布，即運(yùn)用式(13)獲得；然后，式(20)中參數(shù)Σ的先驗(yàn)形式采用逆Wishart分布，即運(yùn)用式(14)獲得.最后，針對(duì)混合二元變量指示器γym采用貝努利(Bernoulli)分布，即假設(shè)其包含的各個(gè)因子效應(yīng)的變量指示器滿足概率pj=P(γj=1)的0-1分布.

由于新模型中加入了混合二元變量指示器，增加了參數(shù)計(jì)算的復(fù)雜度，從而無(wú)法直接得出模型各參數(shù)的封閉后驗(yàn)分布.因此，本文采用MCMC方法動(dòng)態(tài)模擬各參數(shù)后驗(yàn)分布的馬爾科夫鏈，利用各參數(shù)的后驗(yàn)樣本進(jìn)行參數(shù)估計(jì).結(jié)合以上先驗(yàn)假設(shè)，利用R和JAGS軟件實(shí)現(xiàn)相應(yīng)的MCMC過(guò)程，以兩個(gè)響應(yīng)為例，其先驗(yàn)形式如下.

從分布IW(Λ0,ν0)獲取Σ的先驗(yàn)值從分布N(β0,A-1)獲取β的先驗(yàn)值gamma1=(γ1,γ2,…,γ12,…,γ21,γ22,…)for j=1 to p1 do 從分布Bernoulli(gamma1[j])獲取gamma1[j]的先驗(yàn)值end forgamma2=(γ1,γ2,…,γ12,…,γ21,γ22,…)for j=1 to p2 do 從分布Bernoulli(pgamma2[j])獲取gamma2[j]的先驗(yàn)值end for

其中p1和p2分別代表響應(yīng)y1和響應(yīng)y2的候選變量個(gè)數(shù)，gamma1[j](j=1,2,…,p1)和gamma2[j](j=1,2,…,p2)分別代表兩個(gè)響應(yīng)候選變量所對(duì)應(yīng)的變量指示器的先驗(yàn)值，pgamma1[j](j=1,2,…,p1)和pgamma2[j](j=1,2,…,p2) 代表相應(yīng)0-1分布的先驗(yàn)概率值.

該方法主要通過(guò)混合變量指示器和模型的后驗(yàn)概率來(lái)識(shí)別顯著性因子，利用每個(gè)混合變量指示器γ和模型結(jié)構(gòu)M的后驗(yàn)樣本，可以估計(jì)每個(gè)混合變量指示器γ和模型結(jié)構(gòu)M的后驗(yàn)概率，計(jì)算公式如下

(22)

(23)

其中T為迭代總次數(shù)，B為舍棄的燃燒期次數(shù)，γ(t)為該混合變量指示器在第t次迭代的后驗(yàn)樣本值，M(t)為第t次迭代后的模型結(jié)構(gòu).

3.2 多變量過(guò)程能力指數(shù)

工業(yè)生產(chǎn)中常用過(guò)程能力指數(shù)評(píng)估產(chǎn)品處于受控狀態(tài)下的質(zhì)量水平，以體現(xiàn)生產(chǎn)過(guò)程的工藝水平和產(chǎn)品符合性的程度.所謂過(guò)程能力指數(shù)[43]，就是判斷過(guò)程是否滿足規(guī)格要求的一種度量方法，即評(píng)估過(guò)程能力滿足產(chǎn)品規(guī)格要求程度的數(shù)量值.在連續(xù)質(zhì)量改進(jìn)過(guò)程中，計(jì)算與分析過(guò)程能力指數(shù)是至關(guān)重要的一項(xiàng)工作.

質(zhì)量管理學(xué)家朱蘭(Juran)[44]提出了第一代過(guò)程能力指數(shù)Cp，該指數(shù)假設(shè)過(guò)程輸出的均值與目標(biāo)值相重合，具體的定義如下

(24)

其中USL表示規(guī)格上限，LSL表示規(guī)格下限，σ表示過(guò)程特性的標(biāo)準(zhǔn)差.自從朱蘭(Juran)提出第一代過(guò)程能力指數(shù)之后，相繼又有學(xué)者提出了第二代和第三代過(guò)程能力指數(shù)，具體定義見文獻(xiàn)[44].

然而，上述過(guò)程能力指數(shù)只能用來(lái)評(píng)價(jià)單個(gè)質(zhì)量特性過(guò)程的質(zhì)量水平.在實(shí)際生產(chǎn)過(guò)程中通常包含多個(gè)質(zhì)量特性，并且需要考慮各響應(yīng)間可能存在的相關(guān)性，衡量各個(gè)響應(yīng)對(duì)產(chǎn)品質(zhì)量影響的差異，從而全面地評(píng)價(jià)該產(chǎn)品的質(zhì)量和可靠性水平.因此，一些研究人員將單變量過(guò)程能力指數(shù)擴(kuò)展到了多變量情形.Wang等[45]將多變量過(guò)程能力指數(shù)定義為兩塊區(qū)域的比值，具體的表達(dá)式如下

(25)

其中v表示響應(yīng)個(gè)數(shù)；R1表示工藝的規(guī)范區(qū)域面積或體積；R2表示過(guò)程區(qū)域的面積或體積.假設(shè)R1是一個(gè)橢圓形，則R1區(qū)域內(nèi)(包括邊界上)的點(diǎn)都滿足以下不等式[47]

(26)

結(jié)合式(25)和式(26), Chan等[46]提出了以下多變量過(guò)程能力指數(shù)

(27)

其中n表示試驗(yàn)樣本數(shù)，其余參數(shù)與式(26)中的一致.

3.3 結(jié)合SUR模型和因子效應(yīng)原則的連續(xù)質(zhì)量改進(jìn)

本文在SUR模型的框架下結(jié)合因子效應(yīng)原則，運(yùn)用貝葉斯統(tǒng)計(jì)、隨機(jī)搜索技術(shù)以及混合遺傳算法等方法和技術(shù)，從變量篩選、模型構(gòu)建以及參數(shù)優(yōu)化實(shí)現(xiàn)多響應(yīng)的質(zhì)量設(shè)計(jì).通過(guò)結(jié)合因子效應(yīng)原則篩選出符合試驗(yàn)設(shè)計(jì)原則的因子效應(yīng)，不僅能夠?yàn)楹罄m(xù)的響應(yīng)曲面設(shè)計(jì)節(jié)約試驗(yàn)成本，也能夠構(gòu)建更加精確的響應(yīng)曲面模型.響應(yīng)曲面模型構(gòu)建和參數(shù)優(yōu)化則可以視為對(duì)篩選試驗(yàn)分析(即顯著性變量識(shí)別)工作的進(jìn)一步改進(jìn)與優(yōu)化.通過(guò)響應(yīng)曲面建模與優(yōu)化工作，獲得產(chǎn)品或過(guò)程最佳參數(shù)設(shè)計(jì)值，從而能夠從產(chǎn)品形成的源頭減少或控制設(shè)計(jì)值與目標(biāo)值之間的波動(dòng)，提升產(chǎn)品的質(zhì)量與可靠性水平.首先，在SUR模型中針對(duì)每個(gè)因子設(shè)置一個(gè)二元變量指示器，然后運(yùn)用貝葉斯抽樣方法計(jì)算出變量指示器與不同模型結(jié)構(gòu)的后驗(yàn)概率，從而確定顯著性變量和模型結(jié)構(gòu)；其次，根據(jù)所選擇的模型結(jié)構(gòu)，運(yùn)用MCMC方法估計(jì)模型參數(shù)，獲得各響應(yīng)的抽樣值，在此基礎(chǔ)上考慮模型參數(shù)不確定性與預(yù)測(cè)響應(yīng)值的波動(dòng)，構(gòu)建出新的貝葉斯多變量過(guò)程能力指數(shù)；最后，結(jié)合混合遺傳算法對(duì)所構(gòu)建的多變量過(guò)程能力指數(shù)進(jìn)行優(yōu)化，獲得最佳的參數(shù)設(shè)計(jì)值.所提方法的基本流程如圖1所示，具體實(shí)施步驟如下：

步驟1確定試驗(yàn)因子和試驗(yàn)響應(yīng)，選擇合適的試驗(yàn)設(shè)計(jì)開展試驗(yàn)，并收集試驗(yàn)數(shù)據(jù)；

步驟2假設(shè)各個(gè)響應(yīng)的模型結(jié)構(gòu)一致，確定初始的因子組合；

步驟3針對(duì)各個(gè)響應(yīng)的每個(gè)因子設(shè)置一個(gè)二元變量指示器，然后根據(jù)式(19)，構(gòu)建混合二元變量指示器；

步驟4根據(jù)步驟3的研究結(jié)果，構(gòu)建初始的SUR模型.針對(duì)模型參數(shù)設(shè)置相應(yīng)的先驗(yàn)分布，然后利用貝葉斯方法計(jì)算出變量指示器的后驗(yàn)概率，確定模型的顯著性變量，并選擇最佳的模型結(jié)構(gòu)；

步驟5根據(jù)步驟4的結(jié)果，判斷各個(gè)響應(yīng)的模型結(jié)構(gòu)是否一致，如果模型結(jié)構(gòu)一致，則采用SMR模型建模，轉(zhuǎn)入步驟6，否則，采用SUR模型建模，轉(zhuǎn)入步驟6；

步驟6利用MCMC方法對(duì)模型參數(shù)進(jìn)行估計(jì)，獲得各響應(yīng)的后驗(yàn)抽樣值；

步驟7利用步驟6所獲得各響應(yīng)后驗(yàn)抽樣值，構(gòu)建多變量過(guò)程能力指數(shù)；

步驟8針對(duì)所構(gòu)建的多變量過(guò)程能力指數(shù)，通過(guò)混合遺傳算法優(yōu)化獲得最佳的參數(shù)設(shè)計(jì)值.

圖1 本文所提方法的基本流程圖

4 實(shí)例研究

4.1 實(shí)例1

4.1.1 實(shí)例背景

實(shí)例1來(lái)源于文獻(xiàn)[48]，該試驗(yàn)具有兩個(gè)相關(guān)的響應(yīng)，即聚合物的轉(zhuǎn)化效率y1和熱活動(dòng)放射性y2，其中轉(zhuǎn)化效率為望大質(zhì)量特性，取值范圍為[80,100]，目標(biāo)值設(shè)為100；而熱活動(dòng)放射性為望目質(zhì)量特性值，取值范圍為[55,60]，目標(biāo)值設(shè)為57.5.影響上述響應(yīng)的可控因子主要包括：反應(yīng)時(shí)間(x1)、反應(yīng)溫度(x2)和催化劑的使用量(x3).試驗(yàn)者選擇中心復(fù)合設(shè)計(jì)(central composite design，CDD)開展了相關(guān)的試驗(yàn)，其試驗(yàn)計(jì)劃與試驗(yàn)結(jié)果如表1所示.

表1 實(shí)例1的試驗(yàn)計(jì)劃與試驗(yàn)結(jié)果

4.1.2 變量篩選

首先，假設(shè)各響應(yīng)的模型結(jié)構(gòu)相同，則式(2)中的因子效應(yīng)滿足

其次，針對(duì)兩個(gè)響應(yīng)的各個(gè)因子設(shè)置變量指示器γi,i=1,2,…,10 ，在此基礎(chǔ)上結(jié)合因子效應(yīng)原則構(gòu)建混合二元變量指示器，其結(jié)構(gòu)如下

γy1={γ0,γ1,γ2,γ3,γ12×γ1×γ2,γ13×γ1×γ3,

γ23×γ2×γ3,γ12×γ1,γ22×γ2,γ32×γ3}

γy2={γ0,γ1,γ2,γ3,γ12×γ1×γ2,γ13×γ1×γ3,

γ23×γ2×γ3,γ12×γ1,γ22×γ2,γ32×γ3}

本例中因子效應(yīng)數(shù)量較少，此處不考慮效應(yīng)稀疏原則，所以假設(shè)先驗(yàn)信息中的向量pgamma1和pgamma2中的元素均為0.5.然后，運(yùn)用MCMC方法對(duì)每個(gè)參數(shù)迭代了100 000次，舍棄了前20 000次的燃燒期樣本.為了弱化各參數(shù)抽樣值之間自相關(guān)性對(duì)研究結(jié)果的影響，在此對(duì)所獲取的各參數(shù)后驗(yàn)樣本每間隔4步抽樣一次，共20 000個(gè)有效的參數(shù)抽樣值.在此基礎(chǔ)上對(duì)所獲取的樣本進(jìn)行收斂性檢驗(yàn)，并提供一些可視化的收斂性(如蹤跡圖等)來(lái)幫助判斷參數(shù)抽樣值的收斂性.考慮到篇幅限制，在此僅僅提供響應(yīng)y1和響應(yīng)y2的交互效應(yīng)x1x2所對(duì)應(yīng)的參數(shù)后驗(yàn)抽樣值的蹤跡圖如圖2所示.

(a)響應(yīng)y1的交互效應(yīng)x1x2參數(shù)蹤跡圖

(b) 響應(yīng)y2的交互效應(yīng)x1x2參數(shù)蹤跡圖

根據(jù)圖2的蹤跡圖觀察可知：在抽樣過(guò)程中，由于參數(shù)的不確定性在其抽樣值會(huì)圍繞著某個(gè)確定的均值呈現(xiàn)上下波動(dòng)的趨勢(shì)，且其波動(dòng)幅度基本保持一致，呈現(xiàn)出穩(wěn)態(tài)分布的特征，因此可以利用其參數(shù)的抽樣值進(jìn)行后續(xù)的統(tǒng)計(jì)推斷與數(shù)據(jù)分析.在對(duì)模型參數(shù)進(jìn)行收斂性診斷之后，可以計(jì)算變量指示器的后驗(yàn)概率，其研究結(jié)果如表2所示.

表2 混合二元變量指示器的后驗(yàn)概率

4.1.3 響應(yīng)建模與參數(shù)優(yōu)化

變量篩選的研究結(jié)果表明：在本例中，兩個(gè)響應(yīng)的最佳模型為不一致的模型結(jié)構(gòu)，因此后續(xù)將運(yùn)用SUR模型進(jìn)行響應(yīng)曲面建模與優(yōu)化工作.因此式(2)中的因子效應(yīng)向量確定為

對(duì)所構(gòu)建的SUR模型運(yùn)用MCMC方法對(duì)模型參數(shù)進(jìn)行10 000次迭代，并舍棄了前2 000次的燃燒期樣本，獲得了8 000個(gè)模型參數(shù)與響應(yīng)的后驗(yàn)抽樣值.在對(duì)模型參數(shù)和響應(yīng)的抽樣值進(jìn)行收斂性診斷后，運(yùn)用所獲得的響應(yīng)抽樣值構(gòu)建式(27)中的多變量過(guò)程能力指數(shù).其中，樣本數(shù)n設(shè)為8 000，響應(yīng)個(gè)數(shù)v設(shè)為2，目標(biāo)值向量T為[100,57.5],方差-協(xié)方差矩陣Ω利用其響應(yīng)后驗(yàn)抽樣值的樣本方差計(jì)算獲得.事實(shí)上，由于模型參數(shù)的不確定性以及隨機(jī)誤差的影響，其響應(yīng)的后驗(yàn)抽樣值也將呈現(xiàn)出一定的波動(dòng)，因此根據(jù)其響應(yīng)的后驗(yàn)抽樣值及其方差-協(xié)方差矩陣所構(gòu)建的多變量過(guò)程能力指數(shù)將呈現(xiàn)出高度復(fù)雜的非線性、多峰的特征.鑒于上述復(fù)雜的情形，傳統(tǒng)的優(yōu)化算法(如線性規(guī)劃方法等)將難以獲得理想的優(yōu)化結(jié)果.為此，本文擬充分利用遺傳算法的全局搜索與模式搜索的局部?jī)?yōu)化的優(yōu)勢(shì)，構(gòu)建結(jié)合遺傳算法與模式搜索技術(shù)的混合優(yōu)化算法，從而有效地解決上述復(fù)雜曲面函數(shù)的優(yōu)化求解問(wèn)題.在本試驗(yàn)中，將上述多變量過(guò)程能力指數(shù)作為目標(biāo)函數(shù)，運(yùn)用MATLAB優(yōu)化工具箱中的混合遺傳算法對(duì)其進(jìn)行最優(yōu)化求解，其可控因子的參數(shù)設(shè)計(jì)值為(-0.552 0,1.679 9,-0.407 0)，此時(shí)多變量過(guò)程能力指數(shù)等于0.937 0，優(yōu)化過(guò)程見圖3.在上述優(yōu)化結(jié)果中，若固定可控因子x3=-0.407 0，則可以繪制出關(guān)于參數(shù)x1和x2以多變量過(guò)程能力指數(shù)為優(yōu)化目標(biāo)的響應(yīng)曲面，其結(jié)果如圖4所示.從圖4觀察可知：考慮響應(yīng)預(yù)測(cè)值波動(dòng)所構(gòu)建的響應(yīng)曲面往往表現(xiàn)出復(fù)雜的非線性、多峰等特征.

圖3 混合遺傳算法的優(yōu)化過(guò)程

圖4 多變量過(guò)程能力指數(shù)的響應(yīng)曲面圖

4.1.4 不同研究方法的比較與分析

實(shí)例1的試驗(yàn)曾被眾多文獻(xiàn)廣泛應(yīng)用，在此選擇汪建均等[4]，Ko等[37]，Park等[47]，Vining[48]文獻(xiàn)的研究結(jié)果開展比較分析.將上述方法所獲得參數(shù)設(shè)計(jì)值代入到本文所構(gòu)建的SUR模型之中，從而能夠獲得不同方法響應(yīng)預(yù)測(cè)值波動(dòng)的情況，其結(jié)果如圖5所示.

(a)響應(yīng)y1的箱形圖

(b)響應(yīng)y2的箱形圖

從圖5(a)可知：針對(duì)響應(yīng)y1，本文方法的抽樣值波動(dòng)情況與Ko方法、Park方法、Vining方法基本一致，明顯優(yōu)于汪方法的結(jié)果.其主要原因是汪方法采用與本文的初始模型相同的結(jié)構(gòu)，未能進(jìn)一步篩選顯著性變量識(shí)別出最佳的模型結(jié)構(gòu)，從而導(dǎo)致其抽樣結(jié)果明顯偏離目標(biāo)值.從圖5(b)可知：針對(duì)響應(yīng)y2，除了Ko方法以外，本文方法的抽樣值波動(dòng)情況與其他方法比較接近，且都接近目標(biāo)值.需要特別指出的是，雖然Ko方法也采用結(jié)構(gòu)不一致的模型，但是由于其在試驗(yàn)過(guò)程中采用重復(fù)的響應(yīng)數(shù)據(jù)，且需要同時(shí)考慮響應(yīng)與目標(biāo)值的偏差、響應(yīng)的預(yù)測(cè)性能和穩(wěn)健性這三方面的質(zhì)量損失，無(wú)法最小化每一部分的質(zhì)量損失值，所以導(dǎo)致其響應(yīng)y2明顯偏離目標(biāo)值.

考慮到模型參數(shù)不確定性以及響應(yīng)預(yù)測(cè)值的波動(dòng)，將上述方法的最優(yōu)參數(shù)設(shè)計(jì)值代入到所提的多變量過(guò)程能力指數(shù)目標(biāo)函數(shù)之中，并重復(fù)上述抽樣過(guò)程1 000次，其結(jié)果如表3所示.

表3 不同研究方法結(jié)

比較表3的結(jié)果可知：與其他方法比較而言，本文所提方法的多變量過(guò)程能力指數(shù)最大，且兩個(gè)響應(yīng)預(yù)測(cè)值都比較接近相應(yīng)的目標(biāo)值.盡管汪的方法考慮了模型參數(shù)不確定性和預(yù)測(cè)響應(yīng)值的波動(dòng)，但是一方面該方法未實(shí)施變量篩選，而是選擇全因子模型來(lái)進(jìn)行建模.另外該方法將質(zhì)量損失函數(shù)作為目標(biāo)函數(shù)，將響應(yīng)位于規(guī)格限內(nèi)的概率作為約束，得到的優(yōu)化結(jié)果需要最大限度地滿足這兩個(gè)函數(shù)的優(yōu)化目標(biāo).上述原因?qū)е缕淠Ｐ皖A(yù)測(cè)結(jié)果與目標(biāo)值之間存在較大的偏差.因此導(dǎo)致其所獲得的MCpm較低.雖然Ko方法、Vining方法和Park方法針對(duì)兩個(gè)響應(yīng)采用了完全不同的模型結(jié)構(gòu)，但是并未考慮因子效應(yīng)原則，且優(yōu)化函數(shù)與本文也存在差異.此外，上述三種方法在建模過(guò)程中未能考慮模型參數(shù)的不確定性和預(yù)測(cè)響應(yīng)值波動(dòng)的影響，因此導(dǎo)致上述方法的MCpm值仍然低于本文結(jié)果.本文方法不僅結(jié)合因子效應(yīng)原則篩選出顯著性變量，識(shí)別出更加合理的模型，從而提高了模型的預(yù)測(cè)性能，而且利用貝葉斯抽樣方法，考慮了模型參數(shù)不確定性和預(yù)測(cè)響應(yīng)值波動(dòng)對(duì)優(yōu)化結(jié)果的影響，從而獲得較高的MCpm值.

4.2 實(shí)例2

4.2.1 實(shí)例背景

該實(shí)例來(lái)自文獻(xiàn)[49]，主要研究某試驗(yàn)中表面活性劑和乳化變量的設(shè)計(jì)問(wèn)題.在該試驗(yàn)中，重點(diǎn)考察的質(zhì)量特性為粒度(particle size)y1，玻璃轉(zhuǎn)變溫度(glass transition temperature)y2.影響上述兩個(gè)關(guān)鍵質(zhì)量特性的可控因子為普朗尼克F68的含量(x1)、聚氧乙烯40單硬脂酸酯的含量(x2)、聚氧乙烯脫水山梨糖醇脂肪酸酯NF的含量(x3).兩個(gè)響應(yīng)均為望目質(zhì)量特性，其目標(biāo)值分別為150、8；三個(gè)因子的取值范圍均為0≤x1≤1.此試驗(yàn)的目的是希望通過(guò)試驗(yàn)設(shè)計(jì)確定獲得可控因子的最佳參數(shù)設(shè)計(jì).試驗(yàn)者采用中心復(fù)合設(shè)計(jì)開展相關(guān)試驗(yàn)，試驗(yàn)數(shù)據(jù)見表4.

表4 實(shí)例2的試驗(yàn)計(jì)劃及其試驗(yàn)結(jié)果

4.2.2 篩選試驗(yàn)

由于上述試驗(yàn)的次數(shù)較少，因此不適合選擇包含全因子效應(yīng)的初始模型.根據(jù)Wang等[7]，Peterson等[31]的相關(guān)研究，本文選擇的初始模型所包含的因子效應(yīng)向量為

Xj={x1,x2,x3,x1x2,x1x3,x2x3},j=1,2

表5 混合二元變量指示器的后驗(yàn)概率

根據(jù)表5，比較響應(yīng)y1與響應(yīng)y2的顯著性變量可知：響應(yīng)y2模型的因子效應(yīng)x1x2的后驗(yàn)概率較低，是否納入到最終模型有待進(jìn)一步分析.在上述試驗(yàn)中兩個(gè)響應(yīng)模型之間的差異較小，若針對(duì)兩個(gè)響應(yīng)采用相同模型結(jié)構(gòu)對(duì)后續(xù)統(tǒng)計(jì)建模與抽樣更為容易.為此，本文進(jìn)一步比較了兩種不同模型結(jié)構(gòu)之間的后驗(yàn)概率.模型1中兩個(gè)響應(yīng)包含相同的因子效應(yīng){x1,x2,x3,x1x3,x2x3}，其模型后驗(yàn)概率為0.210 2.模型2中響應(yīng)y1包含的因子效應(yīng)為{x1,x2,x3,x1x3,x2x3}，響應(yīng)y2包含的因子效應(yīng)為{x1,x2,x3,x1x2,x1x3,x2x3}.模型2所包含的因子效應(yīng)與運(yùn)用本文所提方法所識(shí)別出模型結(jié)構(gòu)一致，其模型后驗(yàn)概率為0.358 8.盡管模型1與模型2的后驗(yàn)概率存在一定差異，但是為了快速實(shí)現(xiàn)貝葉斯統(tǒng)計(jì)推斷與數(shù)據(jù)分析，同時(shí)比較SUR模型和SMR模型的抽樣效果，在后續(xù)的建模與優(yōu)化中本文傾向選擇相對(duì)簡(jiǎn)單的模型1，從而能夠建立其更為簡(jiǎn)單的SMR模型.

4.2.3 響應(yīng)建模和目標(biāo)優(yōu)化

首先，運(yùn)用SMR模型開展后續(xù)的貝葉斯統(tǒng)計(jì)建模與優(yōu)化，其中式(1)中的因子效應(yīng)向量為

Xj={x1,x2,x3,x1x3,x2x3},j=1,2

然后，運(yùn)用MCMC方法對(duì)SMR模型進(jìn)行抽樣.利用所獲得的響應(yīng)抽樣值構(gòu)建式(27)中的多變量過(guò)程能力指數(shù).其中，樣本數(shù)n設(shè)為10 000，響應(yīng)個(gè)數(shù)v設(shè)為2，目標(biāo)值向量T為[150, 8]，方差-協(xié)方差矩陣Ω利用響應(yīng)抽樣值之間的樣本方差代替.最后，將上述多變量過(guò)程能力指數(shù)作為目標(biāo)函數(shù)，運(yùn)用混合遺傳算法進(jìn)行最優(yōu)化求解，獲得最佳的參數(shù)設(shè)計(jì)值為(0.028 4, 0.499 6, 0.013 7)，目標(biāo)函數(shù)值等于0.999 9.將其他研究方法的優(yōu)化結(jié)果代入到本文所選擇的SMR模型中，則其獲得的抽樣結(jié)果如圖6所示.

比較圖6中不同方法的抽樣結(jié)果可知：Peterson方法的抽樣結(jié)果明顯偏離了目標(biāo)值.其主要原因在于Peterson方法將各個(gè)響應(yīng)位于規(guī)格限內(nèi)的概率作為優(yōu)化函數(shù)，忽略了響應(yīng)與目標(biāo)值之間的偏差，從而導(dǎo)致該方法的后驗(yàn)抽樣值嚴(yán)重偏離目標(biāo)值.Wang方法采用一種有約束的目標(biāo)函數(shù)進(jìn)行優(yōu)化.該方法將質(zhì)量損失作為目標(biāo)函數(shù)，將響應(yīng)位于規(guī)格限內(nèi)的概率作為約束.盡管該方法利用質(zhì)量損失函數(shù)考慮了響應(yīng)預(yù)測(cè)值與目標(biāo)值之間的偏差，但是由于其決策者在選擇成本矩陣以及期望的符合性概率時(shí)往往具有一定的主觀性，從而導(dǎo)致其響應(yīng)抽樣值與目標(biāo)值之間的偏差比本文所提方法更大.

(a)y1的箱形

(b)y2的箱形

另外，為了減少響應(yīng)預(yù)測(cè)值波動(dòng)的影響，將上述方法的最優(yōu)參數(shù)設(shè)計(jì)值代入到所提的多變量過(guò)程指數(shù)目標(biāo)函數(shù)之中，并重復(fù)上述抽樣過(guò)程1 000次，其結(jié)果如表6所示.

表6 不同研究方法結(jié)果

比較表6的研究結(jié)果可知：與其他兩種方法比較而言，本文方法得到的MCpm值最高.上述研究結(jié)果進(jìn)一步驗(yàn)證本文所提方法的有效性與可靠性.

5 討論

5.1 所提方法的優(yōu)勢(shì)

比較而言，本文所提方在處理多響應(yīng)的質(zhì)量設(shè)計(jì)問(wèn)題時(shí)具有一些明顯的優(yōu)勢(shì).首先，本文所提方法在SUR模型的框架下結(jié)合因子效應(yīng)原則(即效應(yīng)系數(shù)原則、 效應(yīng)排序原則和效應(yīng)遺傳原則)篩選出顯著性效應(yīng)，能夠快速排除不符合試驗(yàn)設(shè)計(jì)基本原則的因子效應(yīng)組合，從而大大地縮減模型選擇的空間.因此本文所提方法能夠處理試驗(yàn)因子數(shù)目較大的質(zhì)量設(shè)計(jì)問(wèn)題.其次，在建模過(guò)程中充分結(jié)合貝葉斯統(tǒng)計(jì)的優(yōu)勢(shì)，因此對(duì)樣本量大小沒(méi)有特殊要求，即使是小樣本的試驗(yàn)數(shù)據(jù)也能夠順利開展相關(guān)研究.例如，本文所選擇案例均為三個(gè)因子、兩個(gè)響應(yīng)的質(zhì)量設(shè)計(jì)問(wèn)題，但是其試驗(yàn)數(shù)據(jù)的樣本量差異較大.第一個(gè)案例試驗(yàn)運(yùn)行次數(shù)為20，第二個(gè)案例試驗(yàn)運(yùn)行次數(shù)僅為11，但是從研究結(jié)果來(lái)看，所提方法對(duì)小樣本量的質(zhì)量設(shè)計(jì)問(wèn)題同樣是有效的.最后，本文利用貝葉斯似不相關(guān)(seemingly unrelated regression, SUR)響應(yīng)曲面模型，并在此基礎(chǔ)上運(yùn)用Gibbs抽樣方法獲得多響應(yīng)的抽樣值.通過(guò)上述多響應(yīng)的抽樣值來(lái)估計(jì)響應(yīng)之間的方差-協(xié)方差矩陣能夠充分地考慮由于模型參數(shù)的不確定性以及其他隨機(jī)誤差所引起響應(yīng)波動(dòng)，從而能夠更好地構(gòu)建符合實(shí)際生產(chǎn)過(guò)程的多響應(yīng)優(yōu)化指標(biāo)，即考慮響應(yīng)波動(dòng)的多變量過(guò)程能力指數(shù).這種改進(jìn)的貝葉斯多變量過(guò)程能力指數(shù)，由于結(jié)合貝葉斯響應(yīng)曲面建模與抽樣技術(shù)，因此能夠有效考慮模型參數(shù)不確定性以及其他隨機(jī)誤差對(duì)優(yōu)化結(jié)果的影響.

5.2 所提方法的不足

本文所提方法除了具有上述的優(yōu)點(diǎn)之外，也存在一些明顯的不足之處.首先，所提方法若涉及到大量的響應(yīng)優(yōu)化問(wèn)題(比如10個(gè)或10個(gè)以上)，在運(yùn)行速度方面可能會(huì)受到比較大的影響.由于本文所提方法會(huì)涉及到SUR模型的貝葉斯抽樣技術(shù)，因此當(dāng)響應(yīng)變量為高維數(shù)據(jù)時(shí)在運(yùn)行速度上將會(huì)有很大的影響.其次，當(dāng)響應(yīng)變量個(gè)數(shù)較多時(shí)，在優(yōu)化過(guò)程中可能會(huì)出現(xiàn)響應(yīng)之間的方差-協(xié)方差矩陣行列式接近零情況，從而無(wú)法獲得理想的優(yōu)化結(jié)果.最后，本文所提方法在響應(yīng)曲面建模過(guò)程中通常是假設(shè)響應(yīng)變量滿足為正態(tài)分布.此外，本文利用貝葉斯抽樣方法所構(gòu)建的多變量過(guò)程指數(shù)也僅適用于正態(tài)響應(yīng)的情況.因此本文所提方法無(wú)法有效地處理非正態(tài)、多類別響應(yīng)(如同時(shí)包含離散的屬性變量和連續(xù)性變量)的質(zhì)量設(shè)計(jì)問(wèn)題.

6 結(jié)束語(yǔ)

在連續(xù)性質(zhì)量改進(jìn)中，往往需要同時(shí)考慮變量篩選、響應(yīng)建模和參數(shù)優(yōu)化過(guò)程所涉及到的一系列問(wèn)題.針對(duì)多響應(yīng)的質(zhì)量設(shè)計(jì)問(wèn)題，本文結(jié)合SUR模型和因子效應(yīng)原則，提出了一種連續(xù)性質(zhì)量改進(jìn)方法.該方法在SUR模型融入了因子效應(yīng)原則，提出了一種新的貝葉斯變量選擇方法，不僅保證了模型結(jié)構(gòu)的合理性，而且提高了模型的預(yù)測(cè)性能.另外，本文還采用多變量過(guò)程能力指數(shù)來(lái)衡量過(guò)程能力滿足產(chǎn)品規(guī)格要求的質(zhì)量水平，同時(shí)在優(yōu)化過(guò)程中，利用貝葉斯抽樣技術(shù)考慮了模型參數(shù)的不確定性和預(yù)測(cè)響應(yīng)值的波動(dòng)對(duì)優(yōu)化結(jié)果的影響，從而提高了優(yōu)化結(jié)果的可靠性.最后，通過(guò)兩個(gè)實(shí)例分別展示了SUR模型和SMR模型的抽樣效果.

需要特別指出的是，本文在整個(gè)分析過(guò)程中均假設(shè)響應(yīng)服從多元正態(tài)分布，如何在貝葉斯廣義線性模型的框架下結(jié)合帕累托優(yōu)化策略與灰色關(guān)聯(lián)度分析開展相關(guān)研究解決非正態(tài)多響應(yīng)的質(zhì)量設(shè)計(jì)問(wèn)題，將是下一步重點(diǎn)研究的課題之一.另外，本文所提出的方法僅研究了待測(cè)因子數(shù)小于試驗(yàn)次數(shù)的試驗(yàn)設(shè)計(jì).然而，由于經(jīng)費(fèi)、試驗(yàn)條件等方面的限制，試驗(yàn)者經(jīng)常會(huì)選擇大幅減少試驗(yàn)次數(shù)，采用超飽和試驗(yàn)設(shè)計(jì)(待篩選因子數(shù)目大于試驗(yàn)次數(shù)的情形)，如何選擇合適的方法解決超飽和試驗(yàn)的質(zhì)量設(shè)計(jì)問(wèn)題，也是未來(lái)有待深入研究的課題之一.