何 楨，王 晶,2 ， 李湧 范

(1. 天津大學(xué)管理學(xué)院，天津 300072；2. 南開大學(xué)泰達(dá)學(xué)院，天津 300457；3. 中國質(zhì)量認(rèn)證中心，北京 100070)

響應(yīng)曲面方法(response surface methodology，RSM)最早由Box和Wilson提出，由于其設(shè)計(jì)思想先于田口穩(wěn)健性設(shè)計(jì)方法，所以大多數(shù) RSM 都沒有考慮穩(wěn)健參數(shù)設(shè)計(jì)問題．Vining和 Myers[1]提出了一種將響應(yīng)曲面方法與穩(wěn)健性思想結(jié)合起來的雙響應(yīng)方法，通過對一個(gè)質(zhì)量特性的均值和方差進(jìn)行擬合，根據(jù)擬合的響應(yīng)曲面模型來優(yōu)化均值并減少方差，從而得到穩(wěn)健的過程參數(shù)．然而，該方法主要解決的是單響應(yīng)穩(wěn)健參數(shù)設(shè)計(jì)問題．在現(xiàn)實(shí)生產(chǎn)過程中，經(jīng)常會遇到多響應(yīng)參數(shù)設(shè)計(jì)問題，如手機(jī)、半導(dǎo)體、鋰電池生產(chǎn)、冶金、化工等產(chǎn)品設(shè)計(jì)和過程設(shè)計(jì)中，往往存在多個(gè)質(zhì)量特性值．多響應(yīng)問題通常不存在一組特定的輸入變量，可使所有響應(yīng)變量同時(shí)達(dá)到最優(yōu)．解決多響應(yīng)問題的一個(gè)簡單有效的方法就是重疊等值線圖法，該方法在一個(gè)圖形中為多個(gè)響應(yīng)繪制等值線圖，然后通過觀察找出最優(yōu)解．但這種方法通常只適于輸入變量和響應(yīng)變量很少的情況．另一類解決多響應(yīng)問題的思路是運(yùn)用數(shù)學(xué)的方法將多個(gè)響應(yīng)轉(zhuǎn)化為單個(gè)響應(yīng)，然后再進(jìn)行優(yōu)化，如滿意度函數(shù)法、多元損失函數(shù)法和廣義距離函數(shù)法等．Derringer與Suich[2]提出的滿意度函數(shù)法由于未考慮響應(yīng)之間的相關(guān)性，該方法難以應(yīng)對在諸如化工、電子等行業(yè)中經(jīng)常遇到的多響應(yīng)相關(guān)的情況．Pignatiello[3]、Vining[4]和Ko等[5]提出的質(zhì)量損失函數(shù)法雖然考慮了過程的經(jīng)濟(jì)性，但是其成本矩陣的確定過多依賴于定性分析．廣義距離函數(shù)法是由Khuri和Conlon等[6]提出的，使用的是能夠消除方差和交互作用負(fù)面影響的馬氏距離，因此也稱為馬氏距離函數(shù)法．該方法使用響應(yīng)間的方差-協(xié)方差矩陣，考慮了響應(yīng)間的相關(guān)性，是一種應(yīng)用廣泛的多響應(yīng)優(yōu)化方法，但仍未考慮產(chǎn)品或過程的穩(wěn)健性問題．近年來，雖有很多學(xué)者對該問題進(jìn)行研究[7-8]，但是，有關(guān)如何充分利用 RSM 的響應(yīng)曲面信息，對多響應(yīng)穩(wěn)健性參數(shù)設(shè)計(jì)進(jìn)行的研究還非常有限[9]．

如何解決多響應(yīng)問題中的穩(wěn)健參數(shù)設(shè)計(jì)問題，是當(dāng)今國內(nèi)外質(zhì)量工程領(lǐng)域研究的熱點(diǎn)問題．多響應(yīng)穩(wěn)健參數(shù)設(shè)計(jì)方法發(fā)展至今，大多忽略了可控因子的波動對響應(yīng)的影響．然而在實(shí)際生產(chǎn)過程中，可控因子也存在一定變異，當(dāng)可控因子發(fā)生微小變動時(shí)，若響應(yīng)輸出值的變動較大，則認(rèn)為該產(chǎn)品或過程不穩(wěn)?。P者針對目前研究中存在的不足，考慮響應(yīng)對可控因子波動的穩(wěn)健性，對廣義距離函數(shù)進(jìn)行了改進(jìn)，考慮了多響應(yīng)問題的最優(yōu)性及穩(wěn)健性，提出了穩(wěn)健優(yōu)化總體廣義距離函數(shù)方法，并用該方法對文獻(xiàn)[10]中的實(shí)例進(jìn)行了應(yīng)用驗(yàn)證．

1 最優(yōu)點(diǎn)及穩(wěn)健點(diǎn)的確定

試驗(yàn)區(qū)域中各個(gè)響應(yīng)的最優(yōu)值通常利用Derringer與Suich[2]的滿意度函數(shù)法通過對各個(gè)響應(yīng)單獨(dú)優(yōu)化求得．令 Ti表示目標(biāo)值，Li表示規(guī)格下限，Ui表示規(guī)格上限，則對于望大特性的響應(yīng)變量 yi，其滿意度函數(shù)為

對于望小特性的響應(yīng)變量yi，其滿意度函數(shù)為

對于望目特性的響應(yīng)變量yi，其滿意度函數(shù)為

通過最大化各個(gè)響應(yīng)的滿意度函數(shù)，可以分別得到各響應(yīng)的最優(yōu)點(diǎn)和最優(yōu)值．

在多響應(yīng)穩(wěn)健參數(shù)設(shè)計(jì)中，結(jié)果的最優(yōu)性和穩(wěn)健性同等重要．本文利用穩(wěn)健滿意度函數(shù)法[11]尋找試驗(yàn)區(qū)域中的穩(wěn)健點(diǎn)．對于試驗(yàn)區(qū)域中的任意一點(diǎn)，以該點(diǎn)為中心定義一個(gè)半徑為 r的圓形區(qū)域，用以描述該點(diǎn)處可控因子的波動范圍．在此區(qū)域中可分別計(jì)算各響應(yīng)的極差

則對于響應(yīng)變量yi，其穩(wěn)健滿意度函數(shù)為

對于各個(gè)響應(yīng)變量，使其穩(wěn)健滿意度函數(shù)最大的點(diǎn)可認(rèn)為是試驗(yàn)區(qū)域中該響應(yīng)變量的穩(wěn)健點(diǎn)．在各響應(yīng)的穩(wěn)健點(diǎn)處，響應(yīng)受可控因子波動的影響最小．

2 改進(jìn)的穩(wěn)健優(yōu)化總體廣義距離函數(shù)方法

傳統(tǒng)的廣義距離函數(shù)法使用廣義距離來衡量各響應(yīng)與其最優(yōu)值之間的偏差，用響應(yīng)的方差-協(xié)方差結(jié)構(gòu)來估計(jì)各響應(yīng)距其最優(yōu)值距離的權(quán)重，將多響應(yīng)優(yōu)化問題轉(zhuǎn)化為使廣義距離最小化問題，來確定最佳的操作條件．但是，該方法僅考慮了響應(yīng)的最優(yōu)性問題．考慮到可控因子的波動對響應(yīng)穩(wěn)健性的影響，利用穩(wěn)健滿意度函數(shù)法計(jì)算出各響應(yīng)的穩(wěn)健解，定義一個(gè)穩(wěn)健廣義距離函數(shù)來衡量各個(gè)響應(yīng)與其穩(wěn)健值之間的偏差，然后使其最小化，從而使得多個(gè)響應(yīng)同時(shí)達(dá)到或接近其穩(wěn)健值，考慮了響應(yīng)的穩(wěn)健性要求．然而從統(tǒng)計(jì)分析及質(zhì)量管理的角度出發(fā)，響應(yīng)的最優(yōu)性以及穩(wěn)健性同等重要．穩(wěn)健優(yōu)化總體廣義距離函數(shù)法通過對優(yōu)化廣義距離和穩(wěn)健廣義距離進(jìn)行權(quán)衡，求二者的幾何均值，然后使其最小化，得到穩(wěn)健最優(yōu)的操作條件，使優(yōu)化的響應(yīng)對可控因子的波動具有穩(wěn)健性．

假設(shè)某一生產(chǎn)過程中所有m個(gè)響應(yīng)都依賴于同一組 k個(gè)輸入變量{x1，x2，…，xk}，共計(jì) n次試驗(yàn)，且均可以用同階次的多項(xiàng)式回歸模型表示．則第i個(gè)響應(yīng)的二階模型可以用向量形式表示為

式中：yi為一個(gè)n×1響應(yīng)矩陣，表示第i個(gè)響應(yīng)的n個(gè)觀測值；X為n×p設(shè)計(jì)矩陣，p =(k + 1)(k +2)/2；βi為表示回歸系數(shù)的p×1列向量；εi為一個(gè)n×1列向量，表示第 i個(gè)響應(yīng)的隨機(jī)誤差，服從均值為 0、方差為σii的正態(tài)分布．進(jìn)一步假設(shè)第i個(gè)響應(yīng)的隨機(jī)誤差與第 j個(gè)響應(yīng)的隨機(jī)誤差的協(xié)方差cov(εi,εj)=σijIn，i ≠ j = 1,2,… ,m，其中In為n×n單位矩陣．βi的最小二乘估計(jì)值可推導(dǎo)為=(X ′X )?1，則第 i個(gè)響應(yīng)的二階擬合模型為

var[y?(x ) ] = z ′(x ) (X ′X )?1z ( x)σ i = 1,2,…,m

i ii y?i(x)和y?j(x)的協(xié)方差為

其中ijσ為方差-協(xié)方差矩陣Σ的第(i，j)個(gè)元素．Σ的無偏估計(jì)可表示為

用?Σ替代式(4)中的Σ可以得到? var[()]yx 的無偏估計(jì)值

用 θi表示響應(yīng)yi在試驗(yàn)區(qū)域內(nèi)單獨(dú)優(yōu)化時(shí)得到的最優(yōu)值，θ ＝ [θ1,θ2,… ,θm]′ 表示由這些最優(yōu)值構(gòu)造的m×1向量，則衡量各個(gè)響應(yīng)與其最優(yōu)值之間偏差的優(yōu)化廣義距離函數(shù)為

因?yàn)槭?6)中 v ar[y?(x ) ]是正定矩陣，故而 { var[y?(x)]}?1也是正定矩陣，因此當(dāng)y?(x)≠ θ 時(shí)，距離ρ[y?(x),θ]嚴(yán)格為正．將響應(yīng)的方差-協(xié)方差矩陣的估計(jì)形式(5)代入式(6)中，可得

通過最小化廣義距離 ρ?[y?(x),θ]可確定最佳的操作條件，使各響應(yīng)同時(shí)達(dá)到或接近最優(yōu)值．

用φi表示響應(yīng)yi在試驗(yàn)區(qū)域內(nèi)得到的各自穩(wěn)健點(diǎn)的響應(yīng)值，φ ＝ [φ1,φ2,… ,φm]′ 表示由這些穩(wěn)健值構(gòu)造的 m×1向量，則衡量各個(gè)響應(yīng)與其穩(wěn)健值之間偏差的穩(wěn)健廣義距離函數(shù)為

將式(5)代入上式，可得其估計(jì)值為

通過最小化 ρ?′[y?(x),φ]可確定最穩(wěn)健的操作條件．

考慮多響應(yīng)的最優(yōu)性和穩(wěn)健性的權(quán)衡問題，穩(wěn)健優(yōu)化總體廣義距離表達(dá)式為

式中 w1和 w2分別為優(yōu)化廣義距離和穩(wěn)健廣義距離的權(quán)重，均為正實(shí)數(shù)且滿足 w1+ w2= 1．使穩(wěn)健優(yōu)化總體廣義距離最小的解可認(rèn)為是過程的穩(wěn)健優(yōu)化解．根據(jù)多響應(yīng)問題的最優(yōu)性及穩(wěn)健性的重要程度，權(quán)重w1和w2可選擇不同的取值．當(dāng)w1＝1、w2=0時(shí)，ρoverall=ρ[y?(x),θ]，只考慮多響應(yīng)的最優(yōu)性問題；當(dāng)w1＝0、w2=1 時(shí)，ρoverall= ρ ′[y?(x),φ]，只考慮多響應(yīng)的穩(wěn)健性問題；當(dāng)響應(yīng)的最優(yōu)性及穩(wěn)健性具有相同的權(quán)重，即 w1＝w2=0.5時(shí)，ρoverall= ( ρ [ y? (x ) ,θ ] )12(ρ′[y?(x ),φ] )12，在實(shí)際生產(chǎn)中這種情況應(yīng)用最為廣泛．

3 實(shí)例分析

本文以 Myers與 Montgomery[10]介紹過的一個(gè)經(jīng)典的多響應(yīng)試驗(yàn)設(shè)計(jì)問題為例，利用本文提出的改進(jìn)的穩(wěn)健優(yōu)化總體廣義距離函數(shù)方法對其進(jìn)行計(jì)算，并與傳統(tǒng)方法進(jìn)行比較分析．該試驗(yàn)通過調(diào)整一個(gè)化學(xué)過程的反應(yīng)時(shí)間(x1)和反應(yīng)溫度(x2)，從而使其產(chǎn)量(y1)最大化，黏性(y2)達(dá)到目標(biāo)值 65，同時(shí)使相對分子質(zhì)量(y3)最小化，假定各響應(yīng)變量可接受的范圍分別是 y1≥78.5，62≤y2≤68，y3≤3,300．該試驗(yàn)是一個(gè)含有 5個(gè)中心點(diǎn)的中心復(fù)合設(shè)計(jì)，試驗(yàn)數(shù)據(jù)見表 1．

首先分別擬合出各響應(yīng)的二階響應(yīng)曲面模型為

表1 試驗(yàn)數(shù)據(jù)Tab.1 Test data

方差-協(xié)方差矩陣的無偏估計(jì)值可通過式(3)計(jì)算得到，即

因?yàn)樵撛囼?yàn)設(shè)計(jì)是一個(gè)球形中心復(fù)合設(shè)計(jì)，因此試驗(yàn)區(qū)域是一個(gè)以原點(diǎn)為圓心、半徑為 2的圓，根據(jù)擬合得到的二階響應(yīng)曲面模型，對各響應(yīng)單獨(dú)優(yōu)化得到的最優(yōu)值θi以及取得最優(yōu)值的可控因子的參數(shù)組合如表 2所示．各響應(yīng)在試驗(yàn)區(qū)域內(nèi)得到的穩(wěn)健值iφ及其對應(yīng)的穩(wěn)健點(diǎn)坐標(biāo)見表3．

表2 各響應(yīng)的最優(yōu)值θi與其對應(yīng)點(diǎn)坐標(biāo)Tab.2 Optimal value θi and its corresponding location for each response

表3 各響應(yīng)的穩(wěn)健值iφ及其對應(yīng)點(diǎn)坐標(biāo)Tab.3 Robust value iφ and its corresponding location for each response

首先利用Minitab軟件中多響應(yīng)優(yōu)化的滿意度函數(shù)方法解決該問題．選擇T1＝80作為響應(yīng)產(chǎn)量y1的目標(biāo)值，L1=78.5作為規(guī)格下限，為望大特性．令 T2＝65作為響應(yīng)黏性 y2的目標(biāo)值，L2=62作為規(guī)格下限，U2＝68作為規(guī)格上限，為望目特性．選擇響應(yīng)相對分子質(zhì)量y3的目標(biāo)值T3=3.100，規(guī)格上限為U3＝3.300，為望小特性．各響應(yīng)單獨(dú)滿意度函數(shù)中所有權(quán)重均取 1．若考慮對3個(gè)響應(yīng)同時(shí)優(yōu)化，Minitab給出一組可控變量的最佳參數(shù)組合：x1=-0.157.1，x2=-0.785,7．在該水平下，產(chǎn)量的預(yù)測值為78.76，黏性的預(yù)測值為 66.47，相對分子質(zhì)量的預(yù)測值為3.229.5．但是滿意度函數(shù)法沒有考慮最優(yōu)點(diǎn)處是否穩(wěn)健．在求得的最優(yōu)點(diǎn)(－0.157.1，－0.785.7)處，定義可控因子的變動幅度為 0.01，則在此范圍內(nèi)各響應(yīng)的極差分別為：R1＝0.065.1，R2＝0.215.2，R3＝7.526.2．

然后利用廣義距離函數(shù)法對該問題進(jìn)行分析．前面已經(jīng)求得方差-協(xié)方差矩陣的無偏估計(jì)值和各響應(yīng)單獨(dú)優(yōu)化時(shí)得到的最優(yōu)值θi．利用廣義距離函數(shù)法求得一組可控因子的最優(yōu)參數(shù)組合：x1＝0.048,2，x2=-0.626,2．在該水平下，產(chǎn)量的預(yù)測值為79.26，黏性的預(yù)測值為 68，相對分子質(zhì)量的預(yù)測值為3,300．同樣，若考慮響應(yīng)的穩(wěn)健性，在求得的最優(yōu)點(diǎn)(0.048,2，-0.626,2)處，定義可控因子的變動幅度為0.01，在該區(qū)域內(nèi)各個(gè)響應(yīng)的極差分別為：R1＝0.049,7，R2＝0.158,5，R3＝7.034,3．與滿意度函數(shù)法相比，利用廣義距離函數(shù)法計(jì)算得到的最優(yōu)點(diǎn)處各響應(yīng)的極差均較小，因此可以認(rèn)為各響應(yīng)在(0.048,2，-0.626,2)處，較之在(-0.157,1，-0.785,7)處更穩(wěn)?。?/p>

在實(shí)際生產(chǎn)過程中，找到既穩(wěn)健又可使響應(yīng)達(dá)到最優(yōu)的參數(shù)組合十分重要．下面利用本文提出的穩(wěn)健優(yōu)化總體廣義距離函數(shù)法求解．當(dāng) w1=1、w2＝0時(shí)，得到的結(jié)果與利用傳統(tǒng)廣義距離函數(shù)法的結(jié)果相同．當(dāng) w1＝0、w2＝1時(shí)，求得的可控因子的穩(wěn)健參數(shù)組合為 x1＝0.058,3，x2＝-0.641,6．在該水平下，產(chǎn)量的預(yù)測值為 79.24，黏性的預(yù)測值為 67.89，相對分子質(zhì)量的預(yù)測值為 3,301，各響應(yīng)的極差值為 R1＝0.048,9，R2＝0.161,0，R3＝6.886,3．雖然在該參數(shù)組合下響應(yīng)的極差值最小，但是由于相對分子質(zhì)量的預(yù)測值超出其規(guī)格上限，因此是不可接受的．通常情況下，需要同時(shí)考慮響應(yīng)的最優(yōu)性和穩(wěn)健性問題，因此賦予二者相同的權(quán)重 w1＝w2＝0.5，利用該方法求得一組可控因子的穩(wěn)健最優(yōu)參數(shù)組合：x1＝0.051,2，x2＝-0.649．在該水平下，產(chǎn)量的預(yù)測值為79.22，黏性的預(yù)測值為 67.83，相對分子質(zhì)量的預(yù)測值為3,298.5．當(dāng)可控因子的變動幅度為0.01時(shí)，各響應(yīng)的極差值為 R1＝0.049,5，R2＝0.163,5，R3＝6.904,0．比較利用3種方法計(jì)算得到的結(jié)果，如表4所示．

表4 改進(jìn)的廣義距離函數(shù)方法與傳統(tǒng)方法的結(jié)果對比Tab.4 Comparison of results between improved generalized distance function method and traditional methods

若僅比較預(yù)測的響應(yīng)值，3種方法中利用滿意度函數(shù)法可以得到最好的結(jié)果，但是預(yù)測的響應(yīng)極差較大，說明最優(yōu)點(diǎn)處的穩(wěn)健性較差．雖然另 2種方法得到的響應(yīng)的預(yù)測值稍差，但是響應(yīng)的穩(wěn)健性都要好于滿意度函數(shù)法計(jì)算的結(jié)果．另一方面，利用廣義距離函數(shù)法預(yù)測的響應(yīng) y2和 y3都在響應(yīng)的規(guī)格限上，可控因子的微小波動就可能會導(dǎo)致響應(yīng)變量超出規(guī)格限，造成大量廢品的產(chǎn)生．利用穩(wěn)健優(yōu)化廣義距離函數(shù)進(jìn)行分析：當(dāng) w1=1、w2=0時(shí)，可簡化為傳統(tǒng)廣義距離函數(shù)法；當(dāng)w1=w2=0.5時(shí)，預(yù)測的響應(yīng)y1和y3的極差都要小于傳統(tǒng)方法，表明點(diǎn)(0.051.2，－0.649)處的響應(yīng)y1和y3更穩(wěn)?。m然響應(yīng)y2的極差R2相對廣義距離函數(shù)法更大，但是預(yù)測的響應(yīng)值 67.83與廣義距離函數(shù)法相比更接近目標(biāo)值，不再是響應(yīng)的規(guī)格限，可大大減少廢品的產(chǎn)生．這說明利用改進(jìn)的距離函數(shù)法對多響應(yīng)問題進(jìn)行優(yōu)化時(shí)，權(quán)衡了響應(yīng)的穩(wěn)健性和最優(yōu)性．隨著穩(wěn)健性的權(quán)重w2變大，預(yù)測的響應(yīng)值的極差相應(yīng)變小，但是同時(shí)響應(yīng)值也逐漸遠(yuǎn)離目標(biāo)值，當(dāng) w1=0、w2=1時(shí)，響應(yīng) y3的預(yù)測值 3,301已超出其規(guī)格限 3,300，不能滿足顧客要求．因此在實(shí)際生產(chǎn)中通常選取 w1=w2=0.5，即式(10)取幾何均值作為穩(wěn)健優(yōu)化廣義距離函數(shù)．該方法與傳統(tǒng)方法相比，考慮了響應(yīng)的方差-協(xié)方差結(jié)構(gòu)和可控因子的波動，使優(yōu)化的響應(yīng)對可控因子的波動具有穩(wěn)健性．

4 結(jié) 語

多響應(yīng)穩(wěn)健參數(shù)設(shè)計(jì)是改進(jìn)產(chǎn)品或過程質(zhì)量的一個(gè)重要方法．本文基于廣義距離函數(shù)法對多響應(yīng)穩(wěn)健參數(shù)設(shè)計(jì)進(jìn)行了研究，考慮到響應(yīng)的方差-協(xié)方差結(jié)構(gòu)和可控因子的波動，提出了改進(jìn)的穩(wěn)健優(yōu)化總體廣義距離函數(shù)法，通過對優(yōu)化廣義距離和穩(wěn)健廣義距離進(jìn)行權(quán)衡，將多響應(yīng)穩(wěn)健優(yōu)化問題轉(zhuǎn)化為使穩(wěn)健優(yōu)化總體廣義距離最小化問題，從而得到穩(wěn)健最優(yōu)的操作條件，使優(yōu)化的響應(yīng)對可控因子的波動具有穩(wěn)健性．

［1］ Vining G G，Myers R H. Combining Taguchi and response surface philosophies：A dual response approach[J]. Journal of Quality Technology，1990，22(1)：38-45.

［2］ Derringer G C，Suich R. Simultaneous optimization of several response variables[J]. Journal of Quality Technology，1980，12(4)：214-219.

［3］ Pignatiello J J. Strategies for robust multiresponse quality engineering [J]. IIE Transactions，1993，25(3)：5-15.

［4］ Vining G G. A compromise approach to multiresponse optimization[J]. Journal of Quality Technology，1998，30(4)：309-313.

［5］ Ko Y H，Kim K J，Jun C H. A new loss function-based method for multiresponse optimization[J]. Journal of Quality Technology，2005，37(1)：50-59.

［6］ Khuri A I，Conlon M. Simultaneous optimization of multiple responses represented by polynomial regression functions[J]. Technometrics，1981，23(4)：363-375.

［7］ Park K S，Kim K J. Optimizing multi-response surface problems：How to use multi-objective optimization techniques[J]. IIE Transactions，2005，37(6)：523-532.

［8］ Onur K?ksoy. A nonlinear programming solution to robust multi-response quality problem[J]. Applied Mathematics and Computation，2008，196(2)：603-612.

［9］ Beyer H G，Sendhoff B. Robust optimization—A comprehensive survey[J]. Computer Methods in Applied Mechanics and Engineering，2007，196(33/34)：3190-3218.［10］ Myers R H，Montgomery D C. Response Surface Methodology[M]. 2nd ed. New York：John Wiley & Sons，2002.

［11］ Wang Jing，He Zhen，Oh J H，et al. Multi-response robust optimization using desirability function[C] // IEEE Symposium on Advanced Management of Information for Globalized Enterprises. Tianjin，China，2008：313-315.