利用門限接受法生成均勻設計表

王浩宇

（北京師范大學珠海分校應用數學學院，廣東 珠海 519000）

利用門限接受法生成均勻設計表

王浩宇

（北京師范大學珠海分校應用數學學院，廣東 珠海 519000）

在試驗設計中，均勻設計表的生成通常需要大量的計算并伴有陷入局部最小值的危險.而門限接受法（threshold-accepting algorithm，簡稱TA）的使用可以有效的避免這種情況，從而得到更優解.文章目標在MATLAB上實現門限接受法對均勻設計表的生成，具體包括初始表的選取，局部鄰表的生成，目標函數的確定，以及接受準則的確立等.

均勻設計；門限接受法；局部鄰表；目標函數

科學試驗通常涉及若干因素，且每個因素有若干的水平.一般來說，采用多種參數水平組合進行多次實驗可為問題的解決提供足夠的信息.但是，基于現實方面的諸多考慮，試驗次數不能無限制的增加.于是，如何在有限的試驗次數中獲取足夠信息的問題就擺在面前.均勻設計［1］（Uniform Design）的提出即是為了應對此類問題.其做法是從整個設計空間“均勻”地抽取有限的試驗點，使試驗點具有均勻分布的統計特征，比傳統的試驗設計方法具有更好的穩健性.

試驗設計表的“均勻性”（uniformity）的度量通常用“偏差”（Discrepancy）來進行.在同因素同水平的設計表中，均勻設計表具有最低的偏差，以及最高的均勻性.于是，在已知因素和水平的前提下，均勻設計表的構造或計算問題即是在整個設計空間中，找到偏差值最低表格的優化問題.但是，隨著試驗因素個數和水平個數的增加，設計空間會迅速增大，最終成為一個NP難問題.同時，傳統的優化方法容易使得在整個空間的搜索陷入“局部最優解”而無法達到“全局最優解”.針對這個問題，一些新的優化算法如模擬退火法［2-3］等被提出，其中，門限接受法［4］（Threshold-accepting method，簡稱TA）作為對模擬退火法的一種改進，被證明是一種行之有效的算法.

1 門限接受法

利用門限接受法解決均勻設計表生成的問題，首先需要確定以下幾個要素：①目標函數的選取；②門限序列的確定；③初始設計表格的生成；④局部鄰表的生成規則；⑤接受準則和停止準則等.

圖1表示的是門限接受法的一般流程［5］，其中2表示整個設計空間，xc表示當前設計，J表示單個門限的使用次數，I表示使用門限的個數，T表示使用中的門限，f（x）即目標函數.可見，與一般的優化算法不同，門限接受法的接受準則為Δf≤T，其中在最后一次循環之前T＞0.這也就意味著，即使新表的均勻性低于舊表，舊表仍有可能被替代.這樣做的好處是可以避免陷入“局部最優”，從而在更大的范圍內，尋找若干“局部最優”中的最優解.門限T來自于一個門限序列，通常是一個從大到小的數列，最后一個值為0.在整個搜索過程中，門限被從大到小選取，最后才變成0.如果一開始就設為0，那么在找到第一個局部最優解的時候搜索就停止了.需要注意的是，即使是門限接受法也無法保證最后找到的解是全局最優解，因此最終輸出的表格嚴格意義上并不能成為均勻設計表，但其均勻性已經足夠高，有些甚至可以達到相應偏差的下限，對實際應用不會造成太大影響.

1.1 初始設計表格的生成

正如上節所述，通過門限接受法得到的設計表并不是嚴格意義上的均勻設計表，而是具有良好均勻性的“近似均勻設計表”.因此，為了進一步減少計算次數，可以將搜索區域從整個試驗設計空間縮小到部分均勻性較好的設計表構成的集合，如U型設計表［6］.

圖1 門限接受法生成均勻設計表流程圖Fig.1 Flowchart for the generation of uniform design table using TA method

初始設計表可從轉化后的U型設計中選取.由于沒有有效的證據證明初始表格的均勻性會對最終表格造成影響，在選取過程中可以盡量采取簡化原則，或者利用計算機隨機生成.

1.2 目標函數的選取

一般來說，若利用函數f（X）來度量X的均勻性（或偏差），則其必須滿足以下3個條件：

（1）X中發生行互換或者列互換時，f（X）保持不變；

（2）X中相互對稱的元素各水平取值互換，f（X）保持不變；

（3）將X投射到更低維度，f（X）仍可以計算其均勻性.

在偽蒙特卡羅方法中經常用星LP-偏差（star Lp-discrepancy）來計算偏差［7］，但是星LP-偏差并不滿足上述3個條件.因此，可以利用其變形環繞L2-偏差［8-9］（wrap-around L2-discrepancy，以下簡稱WD2）.若單位超立方Cs表示整個試驗設計空間，其中的n個試驗點構成P＝（xk1，…，xks），k＝1，2，…，n.那么P的WD2-值為

WD2不僅滿足上述3個條件，更重要的是，可以比較容易的得到WD2的下限值［10］，這對之后評估近似均勻設計表有重要的意義.當q為偶數時，設計U（n；qs）的WD2-值的理論下限為

當q為奇數時，理論下限為

需要指出的是，其WD2-值能達到理論下限的設計表并不一定存在.如果搜索過程中發現設計表的WD2-值達到了對應下限，那么搜索過程可以立即停止，此設計表即為嚴格意義上的均勻設計表.但是，更多情況中，WD2-值的理論下限是作為評估最終表格均勻性的一個手段.

1.3 局部鄰表的生成

在搜索過程中，需要不斷地計算當前設計表的局部鄰表，用兩者WD2-值之差與當前門限進行比較，從而決定是否將當前設計替換為鄰表.局部鄰表的生成需要滿足以下2個條件：

（1）新生成的表格仍為U型設計表；

（2）生成過程不宜太復雜，否則會浪費大量計算資源.

本文采取的方法：首先隨機選取當前設計表的1列，再從此列中隨機選取2個元素，若不同，則進行互換；若相同，則再次隨機選取2個元素，以此類推.此方法生成的新表與原表差別很小，適合較為細致的搜索，不容易產生“過度跳躍”.

1.4 門限序列的生成

門限序列中的數值從左到右依次減小，最終變為0.在每一次循環中，當前表格與其鄰表的WD2-值之差與當前門限進行比較，從而決定是否替換，在經過J次循環后，當前門限在門限序列中向右取下一個值，以此類推.門限序列的生成沒有固定的規則，但一般來說，所有門限值的選取必須大小適中，并且門限序列的長度應隨著試驗因素和水平個數的增加而有所增長.本文采取的方法如下所示［11］：

第一步：選取一個初始設計表N，并以此為基礎生成足夠多（由試驗因素和水平個數決定）的鄰表.依次計算這些鄰表的WD2-值，并記最小值為WDmin，最大值為WDmax，設h-range＝WDmax-WDmin；

第二步：生成數列a，第i個元素為ai＝0.1ln（i），i＝1，2，…，n×s；

第三步：生成門限序列T＝h-range*a.

圖2表示對一次優化過程中鄰表WD2-值的追蹤，其中試驗次數n＝25，因素個數s＝10，水平數q＝5.可見收斂過程迅速而平穩.

圖2 WD2-值追蹤圖Fig.2 The tracing of the WD2-value during one convergence process

1.5 設計表的評估

至此，可以利用算法程序對特定試驗次數n，因素個數s以及水平數q的試驗生成相應的近似均勻設計表.由于可以比較容易地得到WD2-值的理論下限，可以對已經生成的表格進行一定程度的評估，結果見表1.

表1 不同元素、水平數和試驗次數下生成表格WD2-值與相應理論下限的對比Table 1 The contrast between WD2-value of generated table and corresponding lower limit of WD2-value

可見，由本文方法生成的近似均勻設計表的偏差值已經十分接近WD2-值的下限，基本上可以作為實際試驗設計用表來使用.

2 結 論

隨著計算機技術的不斷發展，一些NP難問題逐漸可以利用“全局搜索”的方式找到全局最優解，但是在計算資源仍舊緊缺的情況下，利用門限接受法等方法避免落入局部最優從而得到更好的局部最優解仍不失為一種好的替代方式.

［1］ 方開泰.均勻試驗設計的理論、方法和應用——歷史回顧［J］.數理統計與管理，2004，23（3）：69-80.FANG K T.The theory，method and application of uniform experimental design：A historical review［J］.Math Statist Manag，2004，23（3）：69-80.

［2］ KIRKPATRICK S，GELATT C D，VECCHI M P.Optimization by simulated annealing［J］.Science，1983，220（4598）：671-680.

［3］ MORRIS M D，MITCHELL T J.Exploratory design for computational experiments［J］.Statist Plann Infer，1995，43（3）：381-402.

［4］ DUECK G，SCHEUER T.Threshold accepting：A general purpose optimization algorithm appearing superior to simulated annealing［J］.Comput Phys，1990，90（1）：161-175.

［5］ WINKER P，FANG K T.Application of threshold-accepting to the evaluation of the discrepancy of a set of points［J］.Siam J Numer Anal，1997，34（5）：2028-2042.

［6］ GEORGIOU S D，KOUKOUVINOS C，LIU M Q.U-type and column-orthogonal designs for computer experiments［J］.Metrika，2014，77（8）：1057-1073.

［7］ GROSSWALD E，HUA L K，WANG Y.Review：Applications of number theory to numerical analysis［J］.Bull Amer Math Soc，1983，8（3）：489-496.

［8］ HICKERNELL F J.Lattice rules：How well do they measure up？［M］∥HELLEKALEK P，LARCHER G.Random and Quasi-Random Point Sets.New York：Springer-Verlag，1998：109-166.

［9］ FANG K T，MA C X.Wrap-around L2-discrepancy of random sampling，Latin hypercube and uniform designs［J］.Complexity，2001，17（4）：608-624.

［10］FANG K T，TANG Y，YIN J X.Lower bounds for wrap-around L2-discrepancy and constructions of symmetrical uniform designs［J］.Complexity，2005，21（5）：757-771.

［11］FANG K T，LI R，SUDJIANTO A.Design and modeling for computer experiments［M］.London：Chapman＆Hall/CRC，2006：113-117.

Generate uniform design table using threshold accepting algorithm

WANG Hao-yu
（School of Applied Mathematics，Zhuhai Campus of Beijing Normal University，Zhuhai 519000，China）

Using traditional optimization methods，the generation of uniform design table is usually“trapped”in a local minimizer.To avoid this，the threshold-accepting algorithm can be used.This article focuses on realizing this idea in software MATLAB，which includes the choosing of initial table，the generation of a local neighbor design，the object function and the law of accepting in the iterative process.

uniform design；threshold-accepting algorithm；local neighbor；object function

O 212

A

1671-4229（2016）01-0032-04

【責任編輯：周 全】

2015-09-30；

2015-10-23

王浩宇（1988-），男，助教，碩士.E-mail：tkzz0909＠163.com