利用門限接受法生成均勻設(shè)計表

王浩宇

（北京師范大學珠海分校應(yīng)用數(shù)學學院，廣東 珠海 519000）

利用門限接受法生成均勻設(shè)計表

王浩宇

（北京師范大學珠海分校應(yīng)用數(shù)學學院，廣東 珠海 519000）

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

均勻設(shè)計；門限接受法；局部鄰表；目標函數(shù)

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

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

1 門限接受法

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

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

1.1 初始設(shè)計表格的生成

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

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

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

1.2 目標函數(shù)的選取

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

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

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

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

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

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

當q為奇數(shù)時，理論下限為

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

1.3 局部鄰表的生成

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

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

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

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

1.4 門限序列的生成

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

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

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

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

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

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

1.5 設(shè)計表的評估

至此，可以利用算法程序?qū)μ囟ㄔ囼灤螖?shù)n，因素個數(shù)s以及水平數(shù)q的試驗生成相應(yīng)的近似均勻設(shè)計表.由于可以比較容易地得到WD2-值的理論下限，可以對已經(jīng)生成的表格進行一定程度的評估，結(jié)果見表1.

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

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

2 結(jié) 論

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

［1］ 方開泰.均勻試驗設(shè)計的理論、方法和應(yīng)用——歷史回顧［J］.數(shù)理統(tǒng)計與管理，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，F(xiàn)ANG 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