控制圖經濟設計的隨機梯度優(yōu)化方法

張公伯 彭一杰 楊書淮

摘 要：控制圖作為一種實施統(tǒng)計過程控制的工具，在現(xiàn)代制造業(yè)的質量管理中具有廣泛應用?？刂茍D的經濟設計可以提高過程控制效率，降低過程控制成本。本文提出一種基于隨機梯度估計的控制圖經濟設計優(yōu)化方法，首先應用廣義似然比法估計平均運行成本的梯度，并分別結合逐步增加仿真樣本量和隨機化仿真樣本量的隨機近似方法得到控制圖的最優(yōu)控制上界。最后，重點展開介紹了本文提出的控制圖經濟設計方法在休哈特控制圖和指數(shù)加權移動平均控制圖中的應用。

關鍵詞：控制圖;隨機梯度估計;廣義似然比法;隨機近似

中圖分類號：F273.2文獻標識碼：A文章編號：2097-0145（2022）03-0045-08doi：10.11847/fj.41.3.45

Stochastic Gradient Method for Economic Design of Control Charts

ZHANG Gong-bo1， PENG Yi-jie1， YANG Shu-huai2

（1.Guanghua School of Management， Peking University， Beijing 100871， China; 2.College of Engineering， Peking University， Beijing 100871， China）

Abstract：As a tool for statistical process control， control chart has been widely used in quality management in modern manufacturing. The economic design of a control chart can improve the efficiency and lower the cost of process control. This paper proposes a gradient-based simulation optimization method for economic design of control charts. In the first part， the generalized likelihood ratio method is used to obtain the gradient estimate of the average running cost. Then two stochastic approximation methods with increasing sample size in iterations and randomized sample size are utilized respectively to obtain the optimal upper control limit. In the end， the paper applies the proposed method to Shewhart and exponentially weighted moving average control charts.

Key words：control chart; stochastic gradient estimation; generalized likelihood ratio method; stochastic approximation

1 引言

當今世界正高速邁入以數(shù)字化、智能化為特點的新階段，人們對產品和服務質量的重視程度不斷提高。統(tǒng)計過程控制作為一種應用統(tǒng)計方法對生產和服務過程進行評估和監(jiān)控，從而保證產品和服務符合規(guī)定要求的質量管理技術，已廣泛應用于制造業(yè)和服務業(yè)。控制圖是一種實施統(tǒng)計過程控制的工具，它是基于系統(tǒng)運行中測定的關鍵特性值，通過繪制相應的系統(tǒng)特性值統(tǒng)計圖，并根據(jù)假設檢驗原理判斷系統(tǒng)是否處于控制狀態(tài)的統(tǒng)計圖形方法?？刂茍D的概念最早由休哈特[1]提出，在統(tǒng)計過程控制領域中備受學界和業(yè)界關注。近年來，控制圖廣泛應用于集成電路制造[2]、醫(yī)學與流行病學[3]、畜牧業(yè)養(yǎng)殖管理[4]、鐵路接觸網(wǎng)檢測[5]和臨床化學質量控制[6]等領域。

系統(tǒng)平均運行長度（Average Run Length， ARL）是普遍使用的控制圖統(tǒng)計性評價指標，它是系統(tǒng)偏移或宣告失控前所需要的平均觀測周期數(shù)。蒙特卡洛仿真、積分方程和馬爾可夫鏈方法可以估計控制圖的平均運行長度[7]。對控制圖進行優(yōu)化設計和敏感性分析是統(tǒng)計質量控制中具有挑戰(zhàn)性的問題。受生產人員、機械設備、生產環(huán)境和原材料等因素的影響，生產過程會突發(fā)故障處于失控狀態(tài)，將影響制造產品的質量，導致產品異常損耗及故障維修成本?？刂茍D的參數(shù)設計往往會影響樣本的抽樣和檢測、生產過程異常排查及維修等控制措施的成本。在控制圖統(tǒng)計性質的基礎上，需要考慮控制圖的經濟設計，在提高控制圖監(jiān)控效率的同時，降低監(jiān)控過程產生的成本，維持生產過程處于受控狀態(tài)，減少制造階段因產品異常帶來的損失，提升控制圖的經濟效益。Lorenzen和Vance[8]提出了一種靜態(tài)控制圖經濟設計的成本函數(shù)模型。有關可變抽樣區(qū)間和可變樣本容量的動態(tài)控制圖經濟設計的研究，可參考文獻[9～11]。這些研究往往通過馬爾科夫鏈方法估計控制圖的平均運行長度，采用遺傳算法等啟發(fā)式算法優(yōu)化成本目標函數(shù)求解最優(yōu)設計參數(shù)。

控制圖的優(yōu)化問題屬于復雜隨機系統(tǒng)優(yōu)化問題，其目標函數(shù)的解析表達通常未知且含有大量系統(tǒng)設計變量。蒙特卡洛仿真得到的系統(tǒng)樣本表現(xiàn)可以估計目標函數(shù)。梯度包含了目標函數(shù)隨參數(shù)變化的信息，在迭代搜索優(yōu)化算法和敏感性分析中有著重要的應用。通過蒙特卡洛仿真估計的隨機系統(tǒng)中目標函數(shù)關于參數(shù)的梯度稱為隨機梯度估計。常用的無偏隨機梯度估計方法有無窮小擾動分析法[12]，似然比法[13]和弱導數(shù)法[14]。無窮小擾動分析法無法解決樣本表現(xiàn)不連續(xù)的問題，而似然比法和弱導數(shù)法無法解決樣本表現(xiàn)帶有結構參數(shù)（出現(xiàn)在樣本路徑與系統(tǒng)表現(xiàn)度量的函數(shù)關系中的參數(shù)）的問題?？刂茍D的控制邊界屬于結構參數(shù)，且蒙特卡洛仿真得到的系統(tǒng)樣本表現(xiàn)關于控制邊界不連續(xù)，無法保障常用的隨機梯度估計方法的無偏性。Peng等[15，16]提出的廣義似然比法為無偏隨機梯度估計，該方法既適用于帶有結構參數(shù)和不連續(xù)的樣本表現(xiàn)的問題，同時具有無偏性和單樣本軌道性。A67D4071-5767-49FE-868D-67CB7CFF2701

采用隨機梯度估計方法優(yōu)化設計控制圖的現(xiàn)有研究具有局限性。Lele[17]采用有限差分法分別估計指數(shù)加權移動平均和貝葉斯控制圖的梯度，但有限差分法得到的梯度是有偏的，且需要選取攝動值的大小，在處理高維參數(shù)向量時需要很大的計算量。Fu等[18，19]采用平滑擾動分析法得到指數(shù)加權移動平均控制圖的無偏梯度估計和控制圖經濟設計的最優(yōu)控制邊界，但平滑擾動分析法需要解析地計算一個條件期望，在實際應用中需要視控制圖的設定而選取適當?shù)臈l件。Peng等[15，16]采用廣義似然比法得到指數(shù)加權移動平均控制圖的無偏梯度估計，但沒有考慮控制圖的經濟設計問題。本文提出一種考慮過程平均運行長度和控制成本的控制圖經濟性評價指標，并提出一種基于隨機梯度估計的控制圖經濟設計方法。本文首先采用廣義似然比法估計控制圖平均運行成本的梯度，解決了常用的隨機梯度估計方法不適用及需要取條件期望的問題，提供了成本目標函數(shù)隨控制圖設計參數(shù)變化的梯度信息，豐富了對控制圖進行經濟優(yōu)化設計的方法。進而考慮服從正態(tài)分布和均勻分布的質量特性值，通過逐步增加仿真樣本量和隨機化仿真樣本量的隨機近似方法得到控制圖經濟設計的最優(yōu)控制邊界，同時討論了兩種隨機近似方法的理論性質。并以休哈特（Shewart）和指數(shù)加權移動平均（Exponential Weighted Moving Average， EWMA）控制圖的應用為例，說明提出方法的有效性。

2 問題構建

考慮包含單個可測的系統(tǒng)過程變量的控制圖，其中系統(tǒng)過程變量具有穩(wěn)定的受控狀態(tài)和不穩(wěn)定的失控狀態(tài)。系統(tǒng)失控時間τ不可觀測，其累積分布函數(shù)為F0?？刂茍D的采樣間隔時間為Δ。在每個采樣點，可以觀測到系統(tǒng)的檢驗統(tǒng)計量：Y1=X1;Yi=（Xi，Yi-1;θ），i>1，其中θ∈Θ為控制圖設計參數(shù)，Xi為第i個仿真輸出的質量特性值，是與系統(tǒng)參數(shù)無關的馬爾可夫鏈的轉移函數(shù)。系統(tǒng)處于穩(wěn)定的受控狀態(tài)時，Xi的累積分布函數(shù)為F1;系統(tǒng)處于不穩(wěn)定的失控狀態(tài)時，Xi的累積分布函數(shù)為F2。

例1 休哈特控制圖：Yi=Xi，i1。

例2 指數(shù)加權移動平均控制圖：Y1=X1;Yi=αXi+（1-α）Yi-1，i2，0<α1。

令θ1，i和θ2，i分別為第i個檢驗統(tǒng)計量的控制下界和控制上界?？刂茍D的設計參數(shù)為控制界限：θ=θ1，i或θ=θ2，i。當檢驗統(tǒng)計量Yi偏移至控制界限之外時，判定系統(tǒng)處于不穩(wěn)定的失控狀態(tài)。判定系統(tǒng)處于失控狀態(tài)的時間為：

N=min{i：Yi[θ1，i，θ2，i]}，即檢驗統(tǒng)計量Yi不在控制界限內的停時，是一個離散型隨機變量?？刂茍D的平均運行長度是判定系統(tǒng)處于失控狀態(tài)的時間的期望，即

ARL（θ）=E[N]

=∑∞n=1nPr（N=n）=∑∞n=1nE[1（N=n）]（1）

其中1（·）為示性函數(shù)，Pr（·）為括號中事件發(fā)生的概率。ARL（θ）作為控制圖的一種統(tǒng)計性評價指標，沒有考慮成本等經濟因素。如下提出一種考慮系統(tǒng)在失控狀態(tài)下的產品損耗及判定系統(tǒng)處于失控狀態(tài)時的維修成本的控制圖經濟性評價指標。假設系統(tǒng)處于不穩(wěn)定的失控狀態(tài)時，由于制造流程出現(xiàn)異常，生產產品的質量無法得到保障，每單位時間支付損耗成本c;當判定系統(tǒng)處于失控狀態(tài)時，需一次性支付系統(tǒng)維修成本C。當系統(tǒng)完成維修時，系統(tǒng)會恢復正常運行狀態(tài)，因此，該過程可以視為一個再生過程。每進行一次循環(huán)，該過程的成本期望值為：cE[（N-「τ/Δ）+]+C，其中「x為向上取整函數(shù);若x0，（x）+=x，反之，（x）+=0。單位時間內的平均運行成本可以表示為

C（θ）=cE[（N-「τ/Δ）+]+CE[N]

C（θ）是一種控制圖的經濟性評價指標?？梢杂^察到，增大控制圖的控制上界與控制下界之間的距離會導致判定系統(tǒng)處于失控狀態(tài)的頻率降低，從而減少一次性維修成本;但同時會增加系統(tǒng)處于不穩(wěn)定的失控狀態(tài)的時間，導致更大的單位時間產品異常成本。因此，優(yōu)化設計最優(yōu)的控制界限可以實現(xiàn)產品損耗成本與系統(tǒng)維修成本的平衡。在給定產品損耗成本c、系統(tǒng)維修成本C的條件下，通過優(yōu)化設計控制圖的控制界限θ，達到極小化單位時間的平均運行成本C（θ）的目標，即θ*=arg minθC（θ）。隨機近似方法是一種基于梯度的參數(shù)優(yōu)化方法，適用于優(yōu)化求解高維參數(shù)。本文采用隨機近似方法搜尋控制圖的最優(yōu)控制界限θ*。參數(shù)的第k次迭代可以表示為：

θ（k+1）=θ（k）+αkC^k（θ（k）），其中αk為第k次迭代的步長，C^k（θ（k））為C（θ）/θ|θ=θ（k）的估計量

C（θ）θ=cE[（N-「τ/Δ）+]θE[N]-（cE[（N-「τ/Δ）+]+C）E[N]θ（E[N]）2

注意到C（θ）/θ的計算中包含E[N]/θ和E[（N-「τ/Δ）+]/θ。記E[N]/θ=E[G1（θ）]，E[（N-「τ/Δ）+]/θ=E[G2（θ）]。由于通過仿真得到的N關于控制圖的控制界限θ不連續(xù)，無窮小擾動分析法、似然比法和弱導數(shù)法均不適用。本文采用的廣義似然比估計法可以在仿真N的同時，分別得到隨機梯度的單樣本無偏估計量G1（θ）和G2（θ）。將仿真得到的估計量代入C（θ）/θ中，得到梯度估計C^k（·）

C^k（θ）=c∑ki=1G2，i∑ki=1ni-（c∑ki=1（ni-「τ/Δ）++kC）∑ki=1G1，i（∑ki=1ni）2

C^k（·）是由若干估計量構成的非線性函數(shù)，通常不是C（θ）/θ的無偏梯度估計量。梯度估計

C^k（·）與真實梯度C（θ）/θ之間的偏差可以通過增加仿真樣本數(shù)量的方法進行消除，即滿足漸近無偏性：limk→∞ E[C^k（θ）]=C（θ）/θ。

3 基于廣義似然比法的無偏梯度估計

本節(jié)給出當仿真輸出的質量特性值Xi分別服從正態(tài)分布和均勻分布時，采用廣義似然比法得到的E[N]/θ和E[（N-「τ/Δ）+]/θ的單樣本無偏梯度估計量。A67D4071-5767-49FE-868D-67CB7CFF2701

3.1 質量特性值服從正態(tài)分布的廣義似然比估計量

假設質量特性值Xi處于穩(wěn)定的受控狀態(tài)時，服從均值μ1，方差σ21的正態(tài)分布;質量特性值Xi處于不穩(wěn)定的失控狀態(tài)時，服從均值μ2，方差σ22的正態(tài)分布。（1）式中的N可以改寫為

N=min{i：Yi[θ1，θ2]}

=min{i：gi（X1，…，Xi;θ）[0，1]}

其中gi（X1，…，Xi;θ）=（Yi-θ1）/（θ2-θ1），i1。則1（N=n）=（∏n-1i=11（0gi（X1，…，Xi;θ）1））×（1-1（0gn（X1，…，Xn;θ）1）），∏0i=11。取條件于τ/Δ=z，Xi的條件概率密度為：

fi（xi|z）=1（i

記Jg（x（n）;θ）g（x（n）;θ）為g（g1（x1;θ），g2（x1，x2;θ），…，gn（x1，x2，…，xn;θ））′的雅克比矩陣。如果Jg（x（n）;θ）是可逆矩陣且g（x（n）;θ）二階連續(xù)可微，采用廣義似然比法得到

E[N]θ=∑∞n=1nE[1（N=n）]θ=∑∞n=1nE[1（N=n）w（X（n）;θ）]

其中w（X（n）;θ）=lnf（x（n）;θ）/θ+d（x（n）;θ）。d（x（n）;θ）=∑ni=1（J-1g（x（n）;θ）xiJg（x（n）;θ）J-1g（x（n）;θ）ei）′θg（x（n）;θ）-tr（J-1g（x（n）;θ）θJg（x（n）;θ））-（θg（x（n）;θ））′J-1g（x（n）;θ）x（n）ln f（x（n）;θ）。

qJg表示對Jg矩陣中與q有關的元素求導數(shù)，ei為第i個元素為1的單位向量，tr表示矩陣的跡。

x（n）ln f（x（n）;θ）=（ln f（x（n）;θ）/x1，…，ln f（x（n）;θ）/xn）′

=-（xi-μ1σ211（i

采用廣義似然比法分別得到E[N]/θ和E[（N-「τ/Δ）+]/θ的單樣本無偏梯度估計量為：Gnormal1（θ）=N×w（X（n）;θ），Gnormal2（θ）=（N-「τ/Δ）+×w（X（n）;θ）。

例3 考慮常數(shù)的控制圖界限，即：θ1，i=θ1，θ2，i=θ2。在給定控制圖的控制下界θ1的條件下，優(yōu)化設計控制圖的控制上界θ=θ2。當仿真輸出的質量特性值Xi服從正態(tài)分布時，對于指數(shù)加權移動平均控制圖，

w（X（n）;θ）=Nθ2-θ1-（θg（x（n）;θ））′J-1g（x（n）;θ）x（n）ln f（n）（x（n）;z）|n=N，x（n）=X（n），其中θg（x（n）;θ）=-1θ2-θ1（g1（x1;θ），g2（x1，x2;θ），…，gn（x1，x2，…，xn;θ））′，

Jg=1θ2-θ1×11-α… （1-α）n-2（1-α）n-1

0α…α（1-α）n-3α（1-α）n-2

00… αα（1-α）

00… 0αn×n

當α=1時，Jg為單位矩陣，可以相應地得到休哈特控制圖的單樣本無偏梯度估計量。

3.2 質量特性值由均勻分布驅動的廣義似然比估計量

均勻分布隨機數(shù)可以產生服從任意分布的隨機數(shù)，考慮由均勻分布驅動的質量特性值具有更廣泛的應用場景。假設質量特性值Xi處于穩(wěn)定的受控狀態(tài)時的累積分布函數(shù)為F1（·）;質量特性值Xi處于不穩(wěn)定的失控狀態(tài)時的累積分布函數(shù)為F2（·）。仿真輸出的質量特性值由獨立同分布的服從均勻分布U（0，1）的隨機變量構成的向量U（n）=（U1，U2，…，Un）驅動。（1）式中的N可以改寫為：N=min{i：Yi[θ1，θ2]}=min{i：Zi[0，1]}，其中Zi=（Yi-θ1）/（θ2-θ1），i1。{Yi，i1}是一條馬爾科夫鏈，則{Zi，i1}也是一條馬爾可夫鏈：Z1=h1（U1;θ），Zi=（Zi-1，hi（Ui;θ））。（1）式中：1（N=n）=（∏n-1i=11（0Zi1））×（1-1（0Zn1）），∏0i=11。記Jh=Th為h（h1（u;θ），…，hn（u;θ））′的雅克比矩陣。如果Jh是可逆矩陣且h二階連續(xù)可微，采用廣義似然比法得到

E[N]θ=∑∞n=1nE[1（N=n）]θ

=∑∞n=1n∑ni=1E[（∏n-1j=11（0Z+ij1））×（1-1（0Z+in1））×ri（1-;θ）]-

∑∞n=1n∑ni=1E[（∏n-1j=11（0Z-ij1））×（1-1（0Z-in1））×ri（0+;θ）]+

∑∞n=1nE[（∏n-1i=11（0Zi1））×

（1-1（0Zn1））×dn（U（n）;θ）]

其中{Z+ij，j1}和{Z-ij，j1}是分別由U（n）i（U1，…，

1-第i個元素，

…，Un）和U（n）i（U1，…，0-第i個元素，…，Un）通過{Zi，i1}相應生成的馬爾可夫鏈。ri（u;θ）=（J-1h（u;θ）θh（u;θ））′ei，i=1，…，n。dn（u;θ）=∑ni=1ei′J-1h（u;θ）（uiJh（u;θ））J-1h（u;θ）θh（u;θ）-tr（J-1h（u;θ）θJh（u;θ））。qJh表示對矩陣Jh中與q有關的相應元素求導數(shù)，x-和x+分別為x的左極限和右極限。A67D4071-5767-49FE-868D-67CB7CFF2701

對任意函數(shù)h～（·），記h～（x-）limx→x-h～（x），h～（x+）limx→x+h～（x）。

結合文獻[20]的隨機化方法，令N′1為概率質量函數(shù)f～的離散隨機變量，得到

E[N]/θ的單樣本無偏梯度估計量為

Guniform1（θ）=N′×（SN′（U（N′））+DN′（U（N′）））/f～（N′）

E[（N-「τ/Δ）+]/θ的單樣本無偏梯度估計量為

Guniform2（θ）=（N′-「τ/Δ）+×（SN′（U（N′））+

DN′（U（N′）））/f～（N′）

其中Dn（U（n））=（∏n-1i=11（0Zi1））×

（1-1（0Zn1））×dn（U（n）;θ）

Sn（U（n））=∑ni=1（∏n-1j=11（0Z+ij1））×

（1-1（0Z+in1））×ri（1-;θ）-

∑ni=1（∏n-1j=11（0Z-ij1））×

（1-1（0Z-in1））×ri（0+;θ）

例4 考慮例3中的控制圖優(yōu)化設計問題。當仿真輸出的質量特性值Xi的分布由均勻分布驅動時，對于指數(shù)加權移動平均控制圖，馬爾可夫鏈{Zi，i1}可以表示為：Z0=0，Zi=（1-α）Zi-1+hi（Ui;θ），i1。其中hi（Ui;θ）=αθ2-θ1[1（i<τ/Δ）F-11（Ui）+1（iτ/Δ）×F-12（Ui）-θ1]。進而得到

ri（u;θ）=-1θ2-θ1[1（i<τ/Δ）F-11（ui）f1（F-11（ui））+1（iτ/Δ）F-12（ui）f2（F-12（ui））-θ1]

dn（u;θ）=1θ2-θ1∑ni=1[1（i<τ/Δ）F-11（ui）f1（F-11（ui））+

1（iτ/Δ）F-12（ui）f2（F-12（ui））+1]

當α=1時，Zi=hi（Ui;θ），i1，可以相應地得到休哈特控制圖的梯度估計量。

4 基于隨機近似方法的最優(yōu)控制上界求解

本節(jié)通過隨機近似方法求解控制圖經濟設計的最優(yōu)控制上界，在消除梯度估計量與真實梯度之間的偏差的同時，達到極小化控制圖的單位時間內的平均運行成本的目標。考慮逐步增加仿真樣本量和隨機化仿真樣本量的兩種隨機近似方法，并討論相應的理論性質。

4.1 逐步增加仿真樣本數(shù)量

該方法在參數(shù)優(yōu)化迭代時，使用逐步增加仿真樣本數(shù)量的方法計算C^k（θ）。為了證明逐步增加仿真樣本數(shù)量的隨機近似方法的全局最優(yōu)性，提出如下假設：①C（θ）是θ∈R的嚴格凹函數(shù);②C^k（θ）的二階矩一致有界，即：supk∈NE[supθ∈RC^2k（θ）]<∞;且C^k（θ）是漸近無偏的，即：limk→∞E[C^k（θ）]=C（θ）/θ+limk→∞ Γk（θ）=C（θ）/θ，其中Γk（θ）為殘差項;③步長序列與偏差序列滿足：∑∞k=1αk=∞，∑∞k=1α2k<∞，∑∞k=1αkβk<∞和βk=supθ∈R|Γk（θ）|。

定理1 如果條件①～③成立，則limk→∞θ（k）=θ*，a.s.，其中a.s.表示幾乎處處。

證明 令Πk（θ（k）-θ*）2，則E[Πk+1|θ（k），…，θ（1）]=Πk+2αkζk+α2kE[C^2（θ（k））]=Πk+2αk（ζ+k+ζ-k）+α2kE[C^2（θ（k））|θ（k））]，其中ζk（C（θ）/（θ|θ=θ（k）+Γk（θ（k）））（θ（k）-θ*），ζ+k12（ζk+|ζk|），ζ-k12（ζk-|ζk|）。令{kt}是滿足ζ+kt≠0的子序列，則

|C（θ）θ|θ=θ（kt）<|Γkt（θ（kt））|

由條件①和③

limt→∞（θ（kt）-θ*）=0， limk→∞∑∞t=kαtζ+t=0

令ΠkΠk+2∑∞t=kαtζ+t+∑∞t=kα2tE[C^2（θ（t））|θ（k），…，θ（1）]0，{Πk}是正的上鞅，因為E[Πk+1|θ（k），…，θ（1）]=Πk+2αkζ-kΠk。由條件②和③

limk→∞∑∞t=kα2tE[C^2（θ（t））|θ（k），…，θ（1）]=0

根據(jù)鞅收斂定理[21]，存在隨機變量Π，使得

limk→∞Πk=limk→∞Πk=Π， a.s.

注意到

E[Πk]=E[Π1]+2E[∑kt=1αtζt]+∑kt=1α2tE[C^2（θ（t））]

由條件①，對任意δ>0，存在ε>0，使得

sup|θ-θ*|>δ|C（θ）/θ|>ε

若Pr（Π=0）≠1，則

∑∞k=1αkζk=∑∞k=1αt（C（θ）/θ|θ=θ（k）+Γk（θ（k）））（θ（k）-θ*）=-∞以正的概率成立，這與E[Πk]0和∑∞k=1α2kE[C^2（θ（k））]<∞的結論矛盾。因此，Pr（Π=0）=1。綜上所述，定理的結論成立。證畢。

4.2 隨機化仿真樣本數(shù)量

該方法在參數(shù)優(yōu)化迭代時，使用隨機化仿真樣本數(shù)量的方法計算C^k（θ）。具體地，令C^k（θ）∑kt=1Δt（θ），其中Δt（θ）C^t（θ）-C^t-1（θ），C^00。定義C（θ）∑t=1Δt（θ）/Pr（t），其中是與{C^k（θ）}獨立的離散隨機變量。根據(jù)文獻[20]的定理1，C（θ）是C（θ）/θ的無偏梯度估計量，該證明可由limk→∞ E[C^k（θ）]=C（θ）/θ，∑∞t=1E[|Δt（θ）|]<∞，和下述結論得到：如果∑∞t=1E[|C^t-1（θ）-C（θ）/θ|2]/Pr（t）<∞，則C（θ）∈L2，E[C（θ）]=C（θ）/θ，其中L2為平方可積的隨機變量的希爾伯特空間。隨機化仿真樣本數(shù)量的方法可能會導致較大的方差[22]，因此需要推導樣本數(shù)量的最優(yōu)隨機分布。C（θ）的方差可以表示為A67D4071-5767-49FE-868D-67CB7CFF2701

E[C2（θ）]=∑∞t=1vt（θ）/Pr（t）

其中vt（θ）E[|C^t-1（θ）-C（θ）/θ|2]-E[|C^t（θ）-C（θ）/θ|2]。令P（N）為支撐集N的分布族，考慮如下有約束的優(yōu)化問題

infPr∈P（N）：E[]=T{∑∞t=1vt（θ）/Pr（t）}

（2）

優(yōu)化問題（2）本質上是一個最優(yōu)控制問題，該問題的最優(yōu)解即為最優(yōu)隨機分布。

定理2 如果vt（θ）0，vt（θ）vt+1（θ），優(yōu)化問題（2）的解為

Pr*（t）=1，1t

（T+1-t*）vt/2/∑∞t=t*vt，tt*

其中t*min{t∈N：T+1-t∑∞l=tvl/vt}。

證明 考慮如下離散時間的最優(yōu)控制問題

inf{μ（t）∈（0，1）：t∈N}{∑∞t=1vt/μ（t）}

s.t. （t+1）-（t）=μ（t）， t∈N

（1）=0， limt→∞（t）=T（3）

建立漢密爾頓函數(shù)：H（（t），μ（t），y（t+1），t）vt/μ（t）+y（t+1）μ（t），其中（t）為狀態(tài)變量，μ（t）為控制變量。對于t∈N，最優(yōu)控制μ*（n）滿足極大值原理，其中最優(yōu)條件：μ*（t）=infμ∈[0，1）H（（t），μ（t），y（t+1），t），伴隨方程：y（t+1）-y（t）=-Hz（（t），μ（t），y（t+1），t）=0。

根據(jù)伴隨方程，存在v∈R，y*（t）=v，t∈N。對于t0，雖然控制變量μ*（t）=1滿足最優(yōu)條件，但不滿足最優(yōu)控制問題（3）中的約束。因此，v∈R+，且μ*（t）=1，0t

若v1v，則t*=1，limt→∞ *（t）=∑∞t=1vt/v=T，即v=（∑∞t=1vt）2/T2，μ*（t）=Tvt/∑∞l=0vl。

若v1>v，則t*>1，limt→∞ *（t）=t*-1+∑∞t=t*vt/v=T，即v=（∑∞t=t*vt/（T+1-t*））2。因此，t*=min{t∈N：T+1-t∑∞l=tvl/vt}。相似地，可以證明t*>1當且僅當v1>（∑∞l=1vl）2/T2。如果vt非負且非增，則最優(yōu)控制μ*（t）非增，且limt→∞ μ*（t）=0。因此，Pr*（t）=μ*（t）。證畢。

5 數(shù)值實驗

本節(jié)通過數(shù)值實驗介紹本文提出的控制圖經濟設計方法在休哈特和指數(shù)加權移動平均控制圖中的應用。為了方便進行數(shù)值實驗，考慮常數(shù)的控制圖界限，即：θ1，i=θ1，θ2，i=θ2。在給定控制圖的控制下界θ1的條件下，優(yōu)化設計控制圖的控制上界θ=θ2。數(shù)值實驗在搭載英特爾酷睿i5-11300H處理器的Windows系統(tǒng)計算機上通過Matlab R2021a軟件進行測算。

5.1 質量特性值服從正態(tài)分布

假設質量特性值Xi處于穩(wěn)定的受控狀態(tài)時，服從均值為1，方差為1的正態(tài)分布;質量特性值Xi處于不穩(wěn)定的失控狀態(tài)時，服從均值為2，方差為4的正態(tài)分布。設定系統(tǒng)失控時間τ服從均值為20的指數(shù)分布，采樣間隔時間Δ=1，控制圖的控制下界θ1=-1。令損耗成本c=1，系統(tǒng)維修成本C=3?？刂茍D的初始控制上界θ2=1。設定隨機近似方法迭代步長αk=1/k。設定逐步增加仿真樣本量的隨機近似方法中，第k次迭代使用1000+k個仿真樣本;隨機化仿真樣本量的隨機近似方法中，第k次迭代使用的仿真樣本數(shù)量服從成功概率為10-4的幾何分布。對于休哈特控制圖，兩種隨機近似方法分別迭代K=104和K=5×103次。對于不同參數(shù)設定的指數(shù)加權移動平均控制圖：α=0.2，0.4，0.6，0.8，兩種隨機近似方法分別迭代K=2×103和K=103次。表1展示了兩種隨機近似方法迭代估計的最優(yōu)參數(shù)值θ（K）和分別通過108次獨立仿真實驗估計的初始平均運行成本C（θ（0））和最優(yōu)平均運行成本C（θ（K））。圖1展示了對于α=0.2的指數(shù)加權移動平均控制圖，兩種隨機近似方法迭代得到的樣本軌道。

從表1中可以觀察到，兩種隨機近似方法得到相近的最優(yōu)控制上界的估計值。通過對控制圖的控制上界進行優(yōu)化設計，過程的平均運行成本有顯著的下降。雖然逐步增加仿真樣本量的隨機近似方法消耗了更多的迭代次數(shù)，但隨機化仿真樣本量的隨機近似方法在每步迭代時隨機消耗均值為104個仿真樣本，導致運算時間更長。例如，在不同參數(shù)設定的指數(shù)加權移動平均控制圖中，雖然逐步增加仿真樣本量的隨機近似方法的迭代次數(shù)是隨機化仿真樣本量的隨機近似方法的迭代次數(shù)的2倍，但后者的平均運算時間約為前者的3.5倍（前者：約638.47秒;后者：約2220.00秒）。

（a）逐步增加仿真樣本量（b）隨機化仿真樣本量圖1 α=0.2的指數(shù)加權移動平均控制圖控制上界θ2的樣本軌道

5.2 質量特性值由均勻分布驅動

假設質量特性值Xi處于穩(wěn)定的受控狀態(tài)時，服從[-1，1]區(qū)間上的均勻分布，即Xi=2Ui-1;質量特性值Xi處于不穩(wěn)定的失控狀態(tài)時，服從[1，2]區(qū)間上的均勻分布，即Xi=Ui+1。此時，F(xiàn)-11（Ui）=2Ui-1，F(xiàn)-12（Ui）=Ui+1。相應地，

ri（u;θ）=-1θ2-θ1[1（i<τ/Δ）（ui-12）+

1（iτ/Δ）（ui+1）-θ1]

dn（u;θ）=1θ2-θ1∑ni=1[1（i<τ/Δ）（4ui-2）+1（iτ/Δ）（ui+1）+1]

設定系統(tǒng)失控時間τ服從均值為5的指數(shù)分布，其余系統(tǒng)參數(shù)設置與5.1節(jié)的實驗相同。隨機化方法中的離散隨機變量N′服從成功概率為0.1的幾何分布。設定隨機近似方法迭代步長αk=1/k。設定逐步增加仿真樣本量的隨機近似方法中，第k次迭代使用1000+k個仿真樣本;隨機化仿真樣本量的隨機近似方法中，第k次迭代使用的仿真樣本數(shù)量服從成功概率為10-5的幾何分布。對于休哈特控制圖和不同參數(shù)設定的指數(shù)加權移動平均控制圖：α=0.1，0.3，0.5，0.7，兩種隨機近似方法分別迭代K=5×104次和K=104次。表2展示了兩種隨機近似方法迭代估計的最優(yōu)參數(shù)值和分別通過108次獨立仿真實驗估計的初始平均運行成本C（θ（0））和最優(yōu)平均運行成本C（θ（K））。圖2展示了對于休哈特控制圖，兩種隨機近似方法迭代得到的樣本軌道。A67D4071-5767-49FE-868D-67CB7CFF2701

從表2中可以觀察到與表1相同的結論，即：兩種隨機近似方法得到相近的最優(yōu)控制上界的估計值、過程的平均運行成本有顯著的下降。隨機化仿真樣本量的隨機近似方法在每步迭代時隨機消耗均值為105個仿真樣本，導致運算時間更長。在不同參數(shù)設定的指數(shù)加權移動平均控制圖中，雖然逐步增加仿真樣本量的隨機近似方法的迭代次數(shù)是隨機化仿真樣本量的隨機近似方法的迭代次數(shù)的5倍，但后者的平均運算時間約為前者的1.2倍（前者：約2.09×104秒;后者：約2.59×104秒）。

6 結論與啟示

本文基于隨機梯度估計方法提出一種控制圖的經濟設計方法。管理者期望在制造生產流程運作的同時實現(xiàn)控制圖界限優(yōu)化的經濟設計目標。控制圖的經濟性評價指標同時考慮系統(tǒng)失控后的損耗成本和系統(tǒng)判定失控后的維修成本。由于控制圖邊界屬于結構參數(shù)且關于控制圖的系統(tǒng)表現(xiàn)不連續(xù)，本文采用隨機梯度估計方法中的廣義似然比法估計控制圖的經濟性評價指標關于控制圖界限的梯度。結合逐步增加仿真樣本量和隨機化仿真樣本量的隨機近似方法優(yōu)化控制圖界限，并討論了隨機近似方法的理論性質。特別地，本文考慮對控制圖控制上界的優(yōu)化，給出了當質量特性值服從正態(tài)分布和其分布由均勻分布驅動時的梯度估計量。在數(shù)值實驗中，介紹了本文提出的控制圖經濟設計方法在休哈特控制圖和指數(shù)加權移動平均控制圖中的應用，表明了方法的有效性。研究成果為控制圖的經濟優(yōu)化設計提供了一種新的求解思路。未來的研究可以考慮更多影響系統(tǒng)運行成本的因素。采用強化學習方法對控制圖進行經濟優(yōu)化設計也是可以進一步研究的方向。

參 考 文 獻：

[1]Shewhart W A. Economic control of quality of manufactured product[M]. Macmillan And Co Ltd， London， 1931.

[2]Hsieh K L， Tong L I， Wang M C. The application of control chart for defects and defect clustering in IC manufacturing based on fuzzy theory[J]. Expert Systems with Applications， 2007， 32（3）： 765-776.

[3]Sego L H. Applications of control charts in medicine and epidemiology[M]. Virginia： Virginia Polytechnic Institute and State University， 2006.

[4]Stewart B， James R， Knowlton K， et al.. An example of application of process control charts to feed management on dairy farms[J]. The Professional Animal Scientist， 2011， 27（6）： 571-573.

[5]程宏波，王佳鑫，劉杰，等.基于改進多元統(tǒng)計控制圖的接觸網(wǎng)狀態(tài)綜合評價方法[J].鐵道科學與工程學報，2021，18（11）：3048-3056.

[6]Westgard J O， Barry P L， Hunt M R， et al.. A multi-rule shewhart chart for quality control in clinical chemistry[J]. Clin Chem， 1981， 27（3）： 493-501.

[7]Huang W， Shu L， Jiang W. A gradient approach to efficient design and analysis of multivariate EWMA control charts[J]. Journal of Statistical Computation and Simulation， 2018， 88（14）： 2707-2725.

[8]Lorenzen T J， Vance L C. The economic design of control charts： a unified approach[J]. Technometrics， 1986， 28（1）： 3-10.

[9]Bai D S， Lee K T. An economic design of variable sampling interval X control charts[J]. International Journal of Production Economics， 1998， 54（1）： 57-64.

[10]Chen Y K. Economic design of X control charts for non-normal data using variable sampling policy[J]. International Journal of Production Economics， 2004， 92（1）： 61-74.

[11]Naderkhani F， Makis V. Economic design of multivariate Bayesian control chart with two sampling intervals[J]. International Journal of Production Economics， 2016， 174： 29-42.

[12]Ho Y C， Cao X R. Perturbation analysis of discrete-event dynamic systems[M]. Boston： Kluwer Academic Publisher， 1991.A67D4071-5767-49FE-868D-67CB7CFF2701

[13]Rubinstein R Y. The score function approach for sensitivity analysis of computer simulation models[J]. Mathematics and Computers in Simulation， 1986， 28（5）： 351-379.

[14]Pflug G C. Sampling derivatives of probabilities[J]. Computing， 1989， 42（4）： 315-328.

[15]Peng Y J， Fu M C， Hu J Q， et al.. A new unbiased stochastic derivative estimator for discontinuous sample performances with structural parameters[J]. Operations Research， 2018， 66（2）： 487-499.

[16]Peng Y J， Fu M C， Hu J Q， et al.. Generalized likelihood ratio method for stochastic models with uniform random numbers as inputs[J/OL]. Submitted to Operations Research. Preprint in https：//hal.inria.fr/hal02652068/document， 2020.

[17]Lele S. Steady state analysis of three process monitoring procedures in quality control[D]. Michigan： University of Michigan， 1996.

[18]Fu M C， Hu J Q. Efficient design and sensitivity analysis of control charts using monte carlo simulation[J]. Management Science， 1999， 45（3）： 395-413.

[19]Fu M C， Lele S， Vossen T W. Conditional monte carlo gradient estimation in economic design of control limits[J]. Production and Operations Management， 2009， 18（1）： 60-77.

[20]Rhee C H， Glynn P W. Unbiased estimation with square root convergence for SDE models[J]. Operations Research， 2015， 63（5）： 1026-1043.

[21]Karatzas I， Shreve S. Brownian motion and stochastic calculus[M]. Germaoy： Springer Science & Business Media， 2012.

[22]Cui Z， Fu M C， Peng Y， et al.. Optimal unbiased estimation for expected cumulative discounted cost[J]. European Journal of Operational Research， 2020， 286（2）： 604-618.A67D4071-5767-49FE-868D-67CB7CFF2701