摘要
拉丁超立方设计是最常用的计算机试验设计方法之一，针对现有拉丁超立方设计方法采样一次性且难以兼顾设计的空间均匀性和计算效率的问题，提出了一种演化排列拉丁超立方试验设计方法。通过对小样本设计的演化、排列信息继承和扩充等操作，以较小的计算量实现了样本的扩充与优化。此外，所提方法可以兼顾现有样本和新增采样点之间的关系，实现样本的序列扩充，这在实际近似建模过程中十分方便。通过多组数值试验，验证了本文方法在空间均匀性和计算效率等方面的优越性。
Abstract
Latin hypercube design is one of the most commonly used computer experimental design methods, in response to the problem of one-time sampling and difficulty in balancing spatial uniformity and computational efficiency in existing Latin hypercube design methods, an experimental design method of permutation evolution Latin hypercube was proposed. By evolving small-sample designs, inheriting and expanding permutation information, the expansion and optimization of samples were achieved with a relatively small computational effort. In addition, this method gives consideration to the relationship between existing samples and new added samples and realizes the sequence expansion of samples, which is very convenient in the actual approximate modeling process. Through several groups of numerical experiments, the advantages of this method in space filling quality and calculation efficiency were verified.
工程实践中，计算机数值模拟通常用于表征复杂的现实过程[1-3]。但随着设计变量和仿真精度的增加，系统的单次仿真将变得非常耗时。为了利用有限的计算资源，在最大限度上探索系统的全局特性，研究者提出了多种试验设计（design of experiments，DOE）方法[4-5]。DOE可以在设计空间内生成均匀分布的样本，对设计空间的有效探索[6-8]。其中，拉丁超立方设计（Latin hypercube design，LHD）以其优秀的空间均匀性能被广泛应用于试验设计领域。
为提高LHD的性能，研究人员将优化准则与排列算法相结合，提出了一系列适用于不同条件下的优化拉丁超立方设计（optimal LHD，OLHD）方法。试验设计的均匀性可以采用优化准则进行表征，常用的准则主要有最大熵准则、最大距离准则和φp准则等。例如：Morris等利用模拟退火算法构造拉丁超立方设计[9]；Jin等使用增强随机进化（enhanced stochastic evolutionary，ESE）算法来优化设计的空间均匀性能[10]；Viana等通过平移传播（translational propagation algorithm，TPA）算法直接构建LHD[11]；Zhu等提出连续局部枚举（successive local enumeration，SLE）算法[12]；张甜甜等基于数独分组提出一种新型LHD构造方法[13]；陈浩等改进了分片拉丁超立方设计[14]。现有方法在低维小样本条件下效果较好，但随着样本维度和数量的增加，其空间均匀性和计算效率将急剧下降。此外，该方法通常是一次性生成所有样本的单阶采样方法，难以应用于随机性较大的实际工程中。
针对上述问题，本文采用OLHD方法生成初始样本，提出了演化排列拉丁超立方设计（evolution permutation LHD，ELHD）方法。应用序列扩充思想，保证了新增样本和现有样本之间的分布均匀性，以较小的计算代价实现设计的更新与优化，大幅提升实际建模工作效率。
1 优化拉丁超立方设计
针对维度为d、样本数量为N的试验设计问题，采用LHD方法获得分布均匀的样本。LHD的结果可以表示为一个N×d大小的设计矩阵，该矩阵的每一行数据表示一个样本，每一列对应所属变量的不同取值，并且可以被视作是整数1~N的排列。LHD的设计矩阵存在（N!）d种排列方式，这使得随机生成的LHD通常难以保持良好的空间均匀性能。
φp准则被广泛应用于评价样本空间均匀性能，其核心思想为确保样本点间最小距离最大，引入φp准则改进LHD可大幅提高设计空间均匀性。将设计中第i个和第j个采样点之间的距离表示为dij（1≤i＜j≤n），则φp定义为：
(1)
式中，n为采样点个数，p为正整数（本文中p=50），s为不同距离值的数量。
考虑一个N×d的试验设计问题，样本集A={S1，S2，···，SN}的矩阵形式如式（2）所示，Si表示第i个样本点。
(2)
采用φp准则的OLHD方法的一般步骤可总结为：
步骤1：随机生成初始化设计L；
步骤2：在L上执行不同的变换操作从而得到一系列新的设计；
步骤3：计算设计的φp准则值，选择性能最优的结果作为当前最优解；
步骤4：重复前3个步骤，直到满足终止条件。
在处理低维小样本问题时，OLHD方法具有良好的表现。但随着变量维数d和样本数量N的增加，OLHD的计算复杂度（N!）d将急剧增加，其效率和空间均匀性能也会大幅降低。针对OLHD在高维大样本条件下性能不理想的问题，提出基于演化排列方法更新和优化现有设计的ELHD方法，该方法在继承初始设计排列信息的同时可有效降低设计所需计算量。
2 演化排列改进方法
针对样本数量较大的试验设计，采用序列扩充方法可以有效提高采样效率。在样本扩充过程中，性能良好的初始样本包含的排列信息可为后续设计提供参考[15]。基于此，本节提出了一种继承初始样本排列信息，并在该信息指导下完成序列扩充的演化排列拉丁超立方设计方法。
2.1 算法流程
在ELHD中，通过对初始设计执行演化操作以获得扩充所需的种子样本，并采用填充操作和规模调整步骤以满足设计的数量需求。此外，使用排列信息继承算法以较小的计算量更新和优化当前设计。ELHD的流程图如图1所示，其具体步骤如下，并在第2.2节和第2.3节中给出相应步骤的详细说明。
图1ELHD方法流程图
Fig.1Flow chart of the ELHD method
步骤1：利用OLHD方法得到p+1个均匀分布的样本作为初始样本集A1={S1，S2，···，Sp，Sp+1}，其中Si表示第i个样本。
步骤2：通过对A1进行演化操作得到演化矩阵A*1best。
步骤3：A1中的p个间隔将形成一个p×d的试验空间B1。将演化矩阵A*1best叠加到试验空间B1中，得到扩充矩阵A1B1。
步骤4：对A1B1进行2.3节中的基于排列信息继承的更新操作，保持A*1best各列间的变换关系不变，以A*1best的第一维排列为设计变量，以扩充矩阵A1B1的φp值为目标函数进行优化，得到整体样本均匀性优良的设计。
步骤5：判断样本数量2p+1与需求的样本数量n的大小关系，若2p+1＜n则令p=2p+1，重复步骤2到步骤5直至2p+1≥n。
步骤6：进行样本规模调整操作，除去冗余点，完成试验设计。
2.2 演化、扩充和数量调整算法
2.2.1 初始设计
(3)
A1中的p+1个样本点将每个维度的设计域划分为p个区间，这些区间将形成一个新的设计试验空间B1。样本的空间一维投影如图2所示，需要p个种子实现对试验空间B1的填充。
图2样本一维投影示意图
Fig.2One dimensional projection diagram of samples
2.2.2 演化操作
对初始设计执行演化操作以获得填充所需的p个种子。由于初始样本和种子的数量之差仅为1，通过移除初始样本中的某个特定点来生成种子是较为便捷的方法。演化是指遍历初始样本点，使得去除该点后剩余样本的空间均匀性能最佳。通过该操作，实现初始样本排列信息的继承，从而降低设计所需计算量，具体步骤如下：
步骤1：删除初始设计样本点，以获得演化矩阵A*1i。
(4)
步骤2：遍历所有样本点，以min φp为目标，从所有演化矩阵A*1i中选出最优演化矩阵A*1best。
(5)
2.2.3 叠加操作
将A*1best中的样本种子叠加到试验空间B1中，得到扩展矩阵A1B1：
(6)
2.2.4 样本规模调整
若在最后一个循环中所获样本数量大于所需数量，则根据文献[11]中的方法删除距离当前设计空间中心最远的冗余点。
2.3 排列信息继承算法
多种OLHD方法采用列元素交换更新设计。列元素交换是指在设计矩阵的同一列中交换两个不同的元素[11，16]。然而，随着样本规模和维度的增加，列元素交换效率制约着试验设计速度。针对该问题，提出一种基于排列信息继承操作（permutation information inheritance operation，PIO）的列操作方法用于更新和优化设计矩阵。
LHD矩阵的任意两列Pi和Pj均是满足XL≤x≤XU的整数排列，XL和XU分别为设计区间的上下限。
(7)
结合线性代数的相关知识可知，任意Pj均可由Pi经式（8）变换得到。
(8)
式中，Tij是大小为N×N的变换矩阵。
通过线性变换，LHD的各列都可与其第i列相关联。若设置i=1，则整个设计与矩阵第一列相关联。在保持变换矩阵Tij不变的情况下，通过对LHD所得矩阵的第一列执行列元素交换，可以实现整个设计的更新，该操作即为排列信息继承操作。以大小为5×4的设计为例，对其进行排列信息继承操作的结果如图3所示。
图3PIO示意图
Fig.3PIO diagram
显然，在PIO之后，LHD的每一列都发生了变化，而设计仍为LHD。因此，在试验空间B1中执行PIO，并保持扩展矩阵A1B1中初始样本的位置不变。通过继承初始设计的排列属性，以较低的计算成本实现整个设计的更新。在所有结果中，选择空间均匀质量最好的一个作为优化后的样本A1A2：
(9)
式中，xoii表示空间B1中种子的优化排列。
通过该操作，遍历更新LHD的计算量由原先的（N!）d下降至（N!），后续章节将对PIO的效果进行验证。二维ELHD构建过程如下：
步骤1：由ESE算法生成样本数量为11的初始设计A1，如图4（a）所示。
步骤2：对初始设计执行演化操作，得到样本数量为10的最优演化设计A*1best，如图4（b）所示。
步骤3：通过叠加操作得到扩展矩阵A1B1，如图4（c）所示，其中蓝色三角表示叠加的种子，黑色圆点表示已有的采样点。
步骤4：执行PIO以构建优化后的设计A1A2，如图4（d）所示。重复上述过程以获得第二和第三个循环的结果，如图4（e）和图4（f）所示，其中红点表示在该循环中获得的新采样点，黑点表示现有采样点。
步骤5：删除冗余点以获得所需的设计。
图4ELHD的二维示例
Fig.42D example of ELHD
与采用直接优化方法相比，ELHD在扩充样本时实现了新样本与原有样本整体的均匀性。通过排列信息继承算法，将整个LHD的排列与其第一列关联，极大减少了更新样本的计算量。
合理选择初始样本的大小，给定参数p的取值范围，仍需要做进一步讨论。当p过小时，初始样本所包含的排列信息不够全面，容易出现明显的偏差趋势；当p过大时，会给设计的优化求解带来困难。所以，下面将结合数值试验对参数p的影响进行讨论。
2.4 参数p的影响
为了分析参数p的影响，进行维度从2到10、样本数量从180到500的一系列试验。
为避免2.2节中样本规模的调整过程对结果产生影响，在给定设计需求后，通过倒推的方式选择初始样本的大小。例如，要完成一个大小为100的设计，则可选取的初始样本大小为：51、26和14等，与之对应的扩充次数分别为：1次、2次和3次。重复进行20次试验以降低随机因素所造成的影响。此外，为了便于比较，利用式（10）对φp值进行归一化。
(10)
其中，max φp和min φp分别表示同组试验所获取φp的最大值和最小值。
所得结果如图5所示，箱型图的值越小越集中，对应结果的性能越均匀越稳健。在每个系列中均有一个明显的转折点，同一系列中，扩充次数的增加对应着参数p数值的减小。可以看出，随着p的减小，ELHD的稳健性和空间均匀性能呈现出先优后劣的趋势，与2.3节中的分析一致。表1总结了图5中每个系列表现最优的p值。从该表可以看出，当p∈[14，28]时，可获得最佳结果。
图5φp归一化后的箱型图
Fig.5Box diagram of the normalized φp
表1同一系列中性能最佳的参数p
Tab.1 The best performance p in the same series
3 数值试验
本节通过数值算例，比较ESE、TPA和lhsdesign三种算法与ELHD算法在样本空间均匀质量和计算效率上的差异。以上所有算法均是在Windows 10下使用配备Intel Core i7-9750H 2.6 GHz CPU、16 GB RAM和MATLAB R2019b的个人计算机上运行。为降低随机影响，以每种算法独立运行50次所获得的平均值作为结果。
3.1 空间均匀性能比较
为探索ELHD的空间均匀性能，进行了一系列试验。LHD的采样设置如表2所示，本节考虑了Viana等[11]提出的采样设置问题。ELHD以小样本OLHD所得结果为初始设计，结合序列扩充思想后该方法更适用于大中型样本设计。由于高维条件下（维度达到10维及以上）ESE的计算时间过长，本节主要研究维度范围为2到8的大中型采样问题。
表2LHD的采样设置
Tab.2 Sampling settings of LHD
图6给出了大中型采样问题空间均匀性能的比较结果。可以看出，lhdesign的性能最差且随机性较强。TPA在维度较低且样本数量较大时表现良好，随着维度的增加，其竞争力下降。在大多数情况下，ESE与ELHD有着相似的空间均匀性，且明显优于其他两种方法。ELHD在绝大多数问题上有着较优的表现，其优势随着样本数量和维度的增加而更加明显。
图6大中型采样问题对比
Fig.6Comparison of large and medium sampling problems
3.2 计算效率比较
由于ELHD和ESE的空间均匀质量在大多情况下较为接近，且ELHD方法的初始样本通过ESE算法生成，对ELHD和ESE的计算效率进行对比，以验证ELHD性能。
图7显示了不同维度下（d=2，5，10，20），样本数量为100和200时计算效率的比较结果。从图7可以看出，得益于LHD的直接构造，TPA和lhsdesign相较于ELHD耗时更少。然而，TPA和lhsdesign的空间均匀性往往并不令人满意。在满足良好均匀性的前提下，ELHD的计算效率要明显优于ESE。并且，随着维数的增加，ELHD的优势将愈发明显。
图7不同维度下算法耗时的比较
Fig.7Comparison of algorithm time consumption in different dimensions
为了进一步体现ELHD的高效性，下面进行给定初始样本条件下的扩充试验。计算效率体现为同等条件下完成设计所需的工作时间，并通过式（11）计算上述两种方法的耗时比。
(11)
式中，TESE和TELHD分别表示两种算法的计算时间。
设置参数p分别为20和25，在不同维度 （d=2，5，10）上进行序列扩充，同等条件下所获得的耗时比如图8所示。
结果表明，相较于ESE，ELHD的计算效果与样本数量和维度正相关。随着样本数量的增加，ELHD的效率显著提高。此外，在样本数量相同的情况下，维数越高，ELHD的优势越明显。尤其针对800×10的设计，ELHD的计算效率较ESE提升约96%。ELHD兼具空间均匀性能及计算效率，在大中型采样问题中具有显著优势。这主要得益于ELHD通过演化、扩充等操作可以直接实现样本数量的序列扩增。并且，基于排列信息继承操作，整个设计的更新优化问题转变为初始排列信息指导下的单列优化问题，使得优化整个设计的计算复杂度显著降低，从而进一步提高了设计效率。
图8耗时比对比图
Fig.8Comparison of time-consuming ratio
4 结论
本文提出了一种基于演化排列算法的拉丁超立方设计方法。该方法采用排列信息继承算法，以较小的计算代价实现了对整个设计的更新和优化。采用数值算例对ELHD与其他3种常见方法的空间均匀质量进行对比，再选择空间均匀性较好的两种方法进行计算效率的研究。在8个大中型采样问题中进行了空间均匀性能的比较，结果显示，ELHD的空间均匀性能明显优于lhdesign、TPA两种算法，在大部分算例中与ESE性能相近。其后，针对同一条件下ELHD与ESE的计算效率进行对比。结果表明，ELHD的效率明显优于ESE，其优势随着样本维度和数量的增加而提升，在设计为800×10时，ELHD耗时仅需ESE的4%左右。综上所述，ELHD是进行大中型采样的有效方法。并且，其所具有的动态扩充特性使其在实际工程应用中可以更灵活地选择样本数量，从而避免了过拟合以及计算资源浪费等问题，相较于其他一次性采样方法具有明显优势。
