利用門(mén)限接受法生成均勻設(shè)計(jì)表

王浩宇

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

利用門(mén)限接受法生成均勻設(shè)計(jì)表

王浩宇

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

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

均勻設(shè)計(jì)；門(mén)限接受法；局部鄰表；目標(biāo)函數(shù)

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

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

1 門(mén)限接受法

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

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

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

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

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

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

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

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

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

（2）X中相互對(duì)稱(chēng)的元素各水平取值互換，f（X）保持不變；

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

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

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

當(dāng)q為奇數(shù)時(shí)，理論下限為

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

1.3 局部鄰表的生成

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

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

（2）生成過(guò)程不宜太復(fù)雜，否則會(huì)浪費(fèi)大量計(jì)算資源.

本文采取的方法：首先隨機(jī)選取當(dāng)前設(shè)計(jì)表的1列，再?gòu)拇肆兄须S機(jī)選取2個(gè)元素，若不同，則進(jìn)行互換；若相同，則再次隨機(jī)選取2個(gè)元素，以此類(lèi)推.此方法生成的新表與原表差別很小，適合較為細(xì)致的搜索，不容易產(chǎn)生“過(guò)度跳躍”.

1.4 門(mén)限序列的生成

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

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

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

第三步：生成門(mén)限序列T＝h-range*a.

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

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

1.5 設(shè)計(jì)表的評(píng)估

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

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

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

2 結(jié) 論

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

［1］ 方開(kāi)泰.均勻試驗(yàn)設(shè)計(jì)的理論、方法和應(yīng)用——?dú)v史回顧［J］.數(shù)理統(tǒng)計(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，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

【責(zé)任編輯：周 全】

2015-09-30；

2015-10-23

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