一种基于D-最优内表设计的田口试验设计方法

摘要

本发明提供一种基于D‑最优内表设计的田口试验设计方法，该方法的具体步骤是：1.根据试验目的、条件及工程经验，给出试验响应、因子、因子间约束等基本试验信息；2.根据噪声因子的个数以及拟采用的水平个数，选择合适的均匀设计表；3.确定响应与因子之间的关系模型以及内表试验次数；4.在设计区域内，采用D‑最优设计方法给出给定试验次数的内表设计方案，完成内表设计；5.采用基于经验分布函数等分位点方法，完成外表设计；6.根据试验方案进行试验，根据试验得到性能指标计算质量特性信噪比，计算得到使信噪比最大化的可控因子最佳水平；7.基于灵敏度分析法确定稳定因子最佳水平；8.变动调整因子，确定最佳因子水平组合。

权利要求

1.一种基于D-最优内表设计的田口试验设计方法，其特征在于：该方法具体步骤如下：步骤一：根据试验目的、条件以及工程经验，给出试验基本试验信息，包括试验响应即所关注的质量特性、影响因子包括可控因子和噪声因子及其取值范围，各因子之间的约束关系、以及试验资源允许的试验总次数ntotal；

步骤二：根据噪声因子的个数以及拟采用的水平个数m，选择合适的均匀设计表，并以此确定外表试验次数nouter，选择方法简介如下：假设有噪声因子l个，拟采用水平数m，则从均匀设计表库中选择均匀表同时外表试验次数确定为nouter＝ω；

步骤三：确定响应与因子之间的关系模型以及内表试验次数ninner；

该步骤中响应与因子之间的关系模型具体确定方法如下：

基于步骤一中的试验基本信息，采用数学模型的形式表达出响应与因子之间的关系模型，记为y＝Xβ+ε；其中, 为已知的ninner×p维参数模型矩阵；f(xi)为关于xi的已知函数，反映全部可控因子及其之间的约束关系；β＝(β0,β1,β2,...,βp)T为p+

1个待估参数，反映各个因子与响应之间的影响关系；

内表试验次数ninner的确定方法如下：

由一般性等价定理，对于有p个参数的线性模型，存在由n*个试验点组成的D-最优设计ξ*,其中,p≤n*≤p(p+1)/2；所以，根据总试验次数ntotal与外表试验次数nouter，不难给出D-最优内表试验次数范围[p,min(p(p+1)/2,ntotal/nouter)]；

因此，根据实际试验条件，选择任意ninner∈[p,min(p(p+1)/2,ntotal/nouter)]即可；

步骤四：在设计区域内，采用D-最优设计方法给出给定试验次数ninner的内表设计方案，完成内表设计，其具体设计方法如下：首先，将各因子取值范围规范至[-1，1]区间范围内，根据实际因子间的约束关系，确定设计区域；然后，在该设计区域内，选择满足给定次数的设计ξ，使得信息矩阵行列式最大，至此完成D-最优内表设计；进一步，使用G-最优效率衡量该设计的优良性，G-最优效率的定义为：

其中，d(x,ξ)＝fT(x)M-1(ξ)f(x),G-最优效率值Geff越接近于1，则表示设计ξ越好；

步骤五：采用基于经验分布函数等分位点方法，针对内表中每个外表试验点，确定噪声因子水平，完成外表设计；

该步骤中噪声因子水平的确定方法如下：

记x[i]为变量x的基于经验分布函数第i等分位点，其定义如下：其中，F(x)为随机变量x的分布函数，m为噪声因子水平数；argF(·)表示函数F(·)的反函数；换言之，上式表示，当 时，所对应的随机变量x的值记为x[i]，x[i]即为因子x的第i个水平值；

步骤六：根据试验方案进行具体试验，根据试验得到的产品性能指标计算质量特性信噪比，拟合得到信噪比与因子之间的回归函数η(xi)，计算得到使信噪比最大化的可控因子最佳水平；

步骤七:基于灵敏度分析方法确定稳定因子及其最佳水平；

信噪比对影响因子xi,i＝1,2,...,n的灵敏度定义如下：Sij的数值表示 对变量xi的敏感程度，Sij的符号表示 对变量xi的单调性，数值较大的Sij对应的因子为稳定因子，取其最优水平组合以达到系统的稳定性，使系统输出的波动减少到最低程度；其余因子为调整因子，调整系统的输出值使之达到或接近于目标值，从而在尽量减少波动的情况下，得到使得质量特性最接近其目标值的最佳参数组合；

步骤八:回归拟合求得关系模型y中未知参数，并依此对调整因子进行调整，确定最佳因子水平组合，使得在该试验方案下质量特性估计值最接近目标值。

说明书

一种基于D-最优内表设计的田口试验设计方法

技术领域

[0001] 本发明涉及一种基于D-最优内表设计(即最以大化试验点行列式(Determinant)为准则的设计方法)的田口试验设计方法，不仅可以有效解决传统田口设计难以处理的设计区域不规则的情况，而且可有效缓解田口设计采用直积表形式导致试验设计次数剧增的问题，适用于产品设计、生产制造、质量控制以及工艺优化等相关技术领域。

背景技术

[0002] 一个合理的试验设计方案可以在有限的时间、成本等资源条件下，充分揭示各个影响因子对目标响应的影响，在有限的试验次数下获得最优的目标响应结果，为提高产品质量寻求最优的工艺参数组合等。

[0003] 为了提高产品性能的稳定性，田口玄一博士提出了田口设计方法，采用内表、外表直积的形式，在内表中安排可控因子的各个水平，在外表中安排噪声因子进行试验，而所谓直积法，是指对内表中可控因子的每个搭配，使用外表的误差因子模拟各种干扰，计算该搭配的抗干扰能力，即信噪比(SN比)，从而通过对内表设计各种方案的比较，利用SN比寻找最佳的可控因子搭配。

[0004] 广义上说，所有的试验设计本身都是资源受限的设计，试验设计的目的就是在尽量节约资源的情况下，获取尽量多的系统或者过程的信息，优化设计和工艺。在传统的田口设计中，内表、外表均使用正交设计进行，通过方差分析进行数据分析，其设计区域要求是规则的超立方体，一般不考虑因子之间存在约束而导致设计区域不规则的情况。在实际工程中，由于影响因子之间经常是不独立的、存在一定的相互作用，从而产生设计区域不规则的情况，这时，使用正交设计进行内表设计，需切割或填补试验区域进行反复的折衷、对比，不仅设计过程复杂，而且由于降维过多，导致设计效率较低。同时，田口设计的直积表形式，使得在噪声因子较多的情况下，试验次数剧增，经常超出试验资源允许的范围，导致企业难以承受；另外，由于传统的田口内外表设计使用正交设计进行，而正交设计是基于设计表的设计，导致田口设计试验次数不能根据实际需求灵活调整，容易造成资源不能充分利用或资源不足的情况。因此，研究在影响因子存在约束的情况下，合理地选择试验点，从而科学地收集和分析数据，获得更优的参数组合，优化设计和工艺，具有重要的理论意义和迫切的现实需求。

[0005] 为此，本发明给出了一种基于D-最优内表设计的田口试验设计方法。

发明内容

[0006] (1)本发明的目的：本发明针对传统的田口设计在有限的试验次数内，难以灵活高效解决试验区域不规则时的试验设计问题，提供一种基于D-最优内表设计的田口试验设计方法，以筛选出更具有代表性的试验点，从而进行精确、高效、灵活的试验方案设计。

[0007] D-最优设计是一种基于数学模型的设计方法，可以根据模型参数个数，在给定的设计区域内，给出满足一定条件范围内任意试验次数的最优设计；本发明在田口设计的框架下，利用D-最优设计进行内表设计，可以在因子间存在约束的情况下，给出更加灵活、高效的试验设计方案，建立更加贴近产品加工或生产过程实际的数学模型，进行更准确的数据分析和预测，从而确定最佳的因子水平组合。

[0008] (2)技术方案：

[0009] 本发明给出了一种基于D-最优内表设计的田口试验设计方法。

[0010] 田口设计运用内表和外表相结合的方法进行稳健设计，并使用信噪比(SN)作为质量评价指标，以寻求最稳定的参数组合。

[0011] 传统的田口设计，内表采用正交设计安排可控因子的试验，外表采用正交设计安排噪声因子的试验，使用信噪比(SN)度量指标波动的大小。由于正交设计的局限性，当设计区域不规则时，传统的田口设计方法无法给出高效的设计方案。另外，正交试验次数是由设计方法确定的，而非试验精度，在田口设计直积表的作用下，正交设计“整齐可比”的性质使得试验次数剧增，并且调整非常不灵活。然而，在信噪比计算过程中，并不要求试验方案中各误差因素“整齐可比”，因此，可以放松“整齐可比”的要求，探索更为有效的试验设计方法。

[0012] 本发明采用D-最优设计进行内表设计，采用基于经验分布函数等分位点的均匀设* *计进行外表设计。由一般性等价定理，对于有p个参数的线性模型，存在由n (p≤n≤p(p+
1)/2)个试验点组成的D-最优设计ξ*。即D-最优设计可用在一定范围内，给出任意设计区域、任意试验次数的内表试验设计方案，且随着试验因子个数的增加，D-最优设计所需要的最小试验次数将远远小于其他设计。在满足试验精度的前提下，使用D-最优设计进行内表设计，并结合均匀设计进行外表设计，不仅可以有效解决传统田口设计难以处理的因子间存在约束的情况，而且可以在较少的试验次数下，给出高效的试验设计方案。

[0013] 基于上述理论和思路，本发明一种基于D-最优内表设计的田口试验设计方法，具体实施步骤如下：

[0014] 步骤一：根据试验目的、条件以及工程经验，给出试验基本试验信息，包括试验响应(即所关注的质量特性)、影响因子(包括可控因子和噪声因子)及其取值范围，各因子之间的约束关系、以及试验资源允许的试验总次数ntotal。

[0015] 步骤二：根据噪声因子的个数以及拟采用的水平个数m，选择合适的均匀设计表，并以此确定外表试验次数nouter。选择方法简介如下：

[0016] 假设有噪声因子l个，拟采用水平数m，则从均匀设计表库中选择均匀表同时外表试验次数确定为nouter＝ω。

[0017] 步骤三：确定响应与因子之间的关系模型以及内表试验次数ninner。

[0018] 该步骤中响应与因子之间的关系模型具体确定方法如下：

[0019] 基于步骤一中的试验基本信息，采用数学模型的形式表达出响应与因子之间的关系模型，记为y＝Xβ+ε。其中, 为已知的ninner×p维参数模型矩阵；f(xi)为关于xi的已知函数，反应全部可控因子及其之间的约束关系；β＝(β0,β1,β2,...,βp)T为p个待估参数，反应各个因子与响应之间的影响关系。

[0020] 内表试验次数ninner的确定方法如下：

[0021] 由一般性等价定理，对于有p个参数的线性模型，存在由n*(p≤n*≤p(p+1)/2)个试验点组成的D-最优设计ξ*。所以，根据总试验次数ntotal与外表试验次数nouter，不难给出D-最优内表试验次数范围[p,min(p(p+1)/2,ntotal/nouter)]。

[0022] 因此，根据实际试验条件，选择任意ninner∈[p,min(p(p+1)/2,ntotal/nouter)]即可。

[0023] 步骤四：在设计区域内，采用D-最优设计方法给出给定试验次数ninner的内表设计方案，完成内表设计，其具体设计方法如下：

[0024] 首先，将各因子取值范围规范至[-1，1]区间范围内，根据实际因子间的约束关系，确定设计区域；然后，在该设计区域内，选择满足给定次数的设计ξ，使得信息矩阵行列式最大，至此完成D-最优内表设计。进一步，使用G-最优效率衡量该设计的优良性。G-最优效率的定义为：

[0025]

[0026] 这里，G-最优效率值Geff越接近于1，表示设计ξ越好。

[0027] 步骤五：采用基于经验分布函数等分位点方法，针对内表中每个外表试验点，确定噪声因子水平，完成外表设计。

[0028] 该步骤中噪声因子水平的确定方法如下：

[0029] 记x[i]为变量x的基于经验分布函数第i等分位点，其定义如下：

[0030]

[0031] 其中，F(x)为随机变量x的分布函数，m为噪声因子水平数；argF(·)表示函数F(·)的反函数。换言之，上式表示，当 时，所对应的随机变量x的值记为x[i]，x[i]即为因子x的第i个水平值。

[0032] 步骤六：根据试验方案进行具体试验，根据试验得到的产品性能指标计算质量特性信噪比，拟合得到信噪比与因子之间的回归函数η(xi)，计算得到使信噪比最大化的可控因子最佳水平。

[0033] 步骤七:基于灵敏度分析方法确定稳定因子及其最佳水平。

[0034] 信噪比对影响因子xi,i＝1,2,...,n的灵敏度定义如下：

[0035]

[0036] Sij的数值表示 对变量xi的敏感程度，Sij的符号表示 对变量xi的单调性。数值较大的Sij对应的因子为稳定因子，取其最优水平组合以达到系统的稳定性，使系统输出的波动减少到最低程度。其余因子为调整因子，调整系统的输出值使之达到或接近于目标值，从而在尽量减少波动的情况下，得到使得质量特性最接近其目标值的最佳参数组合。

[0037] 步骤八:回归拟合求得关系模型y中未知参数，并依此对调整因子进行调整，确定最佳因子水平组合，使得在该试验方案下质量特性估计值最接近目标值。

[0038] (3)优点和功效：

[0039] 本发明提供一种基于D-最优内表设计的田口试验设计方法，其优点为：

[0040] ①本发利用D-最优设计进行内表设计，可以在试验因子之间存在约束条件下，给出高效合理的田口试验设计方案。

[0041] ②本发明有效缓解了田口设计直积表形式导致的试验次数剧增问题，并且内表可以给出满足一定条件下任意次数的试验，使得试验方案次数可以灵活调整。

附图说明

[0042] 图1是本发明方法流程图。

具体实施方式

[0043] 下面以小球氮化过程为例，结合附图，对本发明做进一步详细说明。

[0044] 反馈弹簧是电液伺服阀的关键部件，其反馈杆上的小球与阀芯上的反馈槽要求处于零间隙配合状态。在反馈弹簧的使用过程中，若小球硬度表面不够，容易引发磨损现象，造成空载流量曲线突变，从而导致错误伺服信号，使电液伺服阀产生致命性故障。为了提高小球表面硬度，需要对小球进行离子软氮化处理。氮化层深度越深，则小球的耐磨性越好，但同时氮化层深度受到小球直径变化的制约。为了保证小球的直径在允许范围内的同时，提高小球表面硬度，设定氮化层深度目标值为0.05mm。见图1，基于本发明具体实施步骤如下：

[0045] 步骤一：确定试验响应(即所关注的质量特性)、影响因子(包括可控因子和噪声因子)及其取值范围，各因子之间的约束关系、以及试验资源允许的试验总次数ntotal。这里，考虑试验资源和进度要求，总试验次数要求不超过150次。

[0046] 将氮化层的深度作为响应，在离子软氮化的过程中，每次都要经历一个溅射阶段和两个氮化阶段，表1汇总了各个影响因子及其取值范围：

[0047] 表1 影响因子及其取值范围

[0048]

[0049] 根据离子软氮化理论，为保证氮化效果，设定约束条件：x5-x2≥20℃，x9-x5≥10℃。另外，Ar,H2,N2通气量(L/min)难以精确控制，存在波动，为噪声因子，其容差分别为[0050] 步骤二：根据噪声因子的个数以及拟采用的水平个数m，选择合适的均匀设计表，并以此确定外表试验次数nouter。选择方法简介如下：

[0051] 假设有噪声因子l个，拟采用水平数m，则从均匀设计表库中，依据偏差选择均匀表或 且外表试验次数确定为nouter＝ω。

[0052] 假设有噪声因子l个，拟采用水平数m，则从均匀设计表库中选择均匀表Uω(ml)，且外表试验次数确定为nouter＝ω。

[0053] 在本例中，共有3个噪声因子，为在保证精度的前提下减少试验次数，对每个因子选取6个水平进行试验，选用均匀设计表 进行外表设计，外表试验次数nouter＝6。

[0054] 步骤三：确定响应与因子之间的关系模型以及内表试验次数ninner。

[0055] 该步骤中响应与因子之间的关系模型具体确定方法如下：

[0056] 基于步骤一中的试验基本信息，采用数学模型的形式表达出响应与因子之间的关系模型，记为y＝Xβ+ε。其中, 为已知的ninner×p维参数模型矩阵；fT
(xi)为关于xi的已知函数，反应全部可控因子及其之间的约束关系；β＝(β0,β1,β2,...,βp)为p个待估参数，反应各个因子与响应之间的影响关系。

[0057] 内表试验次数ninner的确定方法如下：

[0058] 由一般性等价定理，对于有p个参数的线性模型，存在由n*(p≤n*≤p(p+1)/2)个试* total outer验点组成的D-最优设计ξ。所以，根据总试验次数n 与外表试验次数n ，不难给出D-最优内表试验次数范围[p,min(p(p+1)/2,ntotal/nouter)]。因此，根据实际试验条件，选择任意ninner∈[p,min(p(p+1)/2,ntotal/nouter)]即可。

[0059] 在本实施例中，根据工程经验，选取如下模型：

[0060] y＝β0+x1β1+x2β2+xβ3+x4β4+x5β5+x6β6+x7β7+x8β8+x9β9+x2x3β23+[0061] x2x5β25+x2x8β28+x3x5β35+x3x6β36+x3x7β37+x3x8β38+x3x9β39+x6x9β69+e[0062] 其中，误差e～N(0,σ2)。

[0063] 该模型中共有19个未知参数，即p＝19。结合内表所需次数，可确定外表试验次数的取值范围19≤ninner≤min(19(19+1)/2,ntotal/6)。考虑总试验资源，选用内表试验次数ninner＝19。

[0064] 步骤四：在设计区域内，采用D-最优设计方法给出给定试验次数ninner的内表设计方案，完成内表设计。具体设计方法如下：

[0065] 首先，将各因子取值范围规范至[-1，1]区间范围内，根据实际因子间的约束关系，确定设计区域；然后，在该设计区域内，选择满足给定次数的设计ξ，使得信息矩阵行列式最大，其中，f(xi)为关于xi的已知函数，至此完成D-最优内表设计。进一步，使用G-最优效率衡量该设计的优良性。G-最优效率的定义为：

[0066]

[0067] 其中，d(x,ξ)＝fT(x)M-1(ξ)f(x)。G-最优效率值Geff越接近于1，则表示设计ξ越好。

[0068] 在本实施例中，D-最优内表设计表如下所示：

[0069] 表2 基于D-最优设计的内表试验方案(编码)

[0070]序号 x1 x2 x3 x4 x5 x6 x7 x8 x9
1 -1 -1 -1 -1 -1 -1 1 -0.38418 1
2 -1 -0.12675 1 -1 0.66667 1 1 -1 1
3 0.02211 -0.19135 0.001343 0.021786 -0.1083 -0.01131 0.050151 -0.03109 1
4 1 1 0.022368 1 0.66667 -0.00126 1 -0.01406 1
5 -1 1 1 -1 -0.33333 -1 -1 1 1
6 1 0.31204 0.036341 -1 -1 1 1 -1 -0.66667
7 -1 1 0.073784 1 0.66667 1 -1 1 1
8 1 -1 -1 -0.15154 -1 1 1 0.025892 1
9 0.072806 -1 -1 1 -1 -0.03687 -1 -1 -0.66667
10 1 0.18165 1 1 -1 -1 1 -1 -0.66667
11 -1 1 -1 -1 0.66667 -1 -0.10588 -0.045 1
12 -1 -1 0.020447 1 -0.09342 -1 -1 -1 0.23992
13 -1 -0.03265 1 1 -1 1 1 1 1
14 1 -1 1 -1 -1 -1 1 -0.00697 0.27705
15 1 -1 -0.26793 -1 -1 -1 -1 -1 1
16 -1 -1 1 1 0.66667 -1 -1 1 1

[0071]17 1 0.020626 0.074386 1 -1 -1 -1 0.055004 1
18 -1 0.33333 1 -0.04891 -1 1 1 1 -0.66667
19 -0.08383 1 1 0.048094 0.66667 -1 -1 -1 1

[0072] 该设计的G效率 可见其设计效率较高。

[0073] 步骤五：采用基于经验分布函数等分位点方法，针对内表中每个外表试验点，确定噪声因子水平，完成外表设计。

[0074] 该步骤中噪声因子水平的确定方法如下：

[0075] 记x[i]为变量x的基于经验分布函数第i等分位点，其定义如下：

[0076]

[0077] 其中，F(x)为随机变量x的分布函数，m为噪声因子水平数；argF(·)表示函数F(·)的反函数。换言之，上式表示，当 时，所对应的随机变量x的值记为x[i]，x[i]即为因子x的第i个水平值。

[0078] 本例中，假定噪声因子均服从正态分布N(0,σ2)，容差Δ与标准差σ之间存在关系Δ＝3σ，可得三个噪声因子的方差分别为：均值μ则根据相应的内表确定。以表2中2号内表为例，该条件下确定的三个噪声因子水平数见表3。

[0079] 表3 2号内表对应噪声因子水平数设定

[0080]

[0081] 则对应的外表设计为：

[0082] 表4 2号内表对应的外表设计

[0083]

[0084] 步骤六：根据试验方案进行具体试验，根据试验得到的产品性能指标计算质量特性信噪比，拟合得到信噪比与因子之间的回归函数η(xi)，计算得到使信噪比最大化的可控因子最佳水平。

[0085] 该例中，所关注的响应氮化层深度呈现望目特性，根据所得试验数据，计算各个条件下的信噪比见表5

[0086] 表5 试验设计方案(解码)及信噪比

[0087]

[0088]

[0089] 利用上述数据，对信噪比η与各因子进行模型参数估计，采用逐步回归法可得到以下模型

[0090] η＝3465.311+0.033862x1-1.34615x2+16.60486x3-0.1968x4-6.99205x5-384.59x6[0091] +84.20613x7+0.3755x8-5.12495x9+0.148254x1x3+0.002458x2x5+0.328684x4x6[0092] +0.255114x4x7+2.575394x6x8+0.88726x7x8+0.010374x5x9

[0093] 利用非线性求解法，求出在满足试验约束的情况下使得信噪比最大的因子水平组合为：

[0094] x*＝[30,500,0.6,120,520,0.07,0.475,150,530]

[0095] 步骤七:基于灵敏度分析方法确定稳定因子及其最佳水平。

[0096] 信噪比对影响因子xi,i＝1,2,...,n的灵敏度定义如下：

[0097]

[0098] Sij的数值表示 对变量xi的敏感程度，Sij的符号表示 对变量xi的单调性。数值较大的Sij对应的因子为稳定因子，取其最优水平组合以达到系统的稳定性，使系统输出的波动减少到最低程度。其余因子为调整因子，调整系统的输出值使之达到或接近于目标值，从而在尽量减少波动的情况下，得到使得质量特性最接近其目标值的最佳参数组合。

[0099] 本实施例中，目标函数η对各个影响因子的灵敏度计算如下：

[0100]

[0101]

[0102]

[0103]

[0104]

[0105]

[0106]

[0107]

[0108]

[0109] 目标函数η对所有影响因子的灵敏度自大至小汇总为：

[0110] 表6 η对各影响因子灵敏度汇总表

[0111]

[0112] 显然，目标函数η对N2通气量最为敏感，其次是H2通气量和Ar通气量。因此，确定N2、H2Ar通气量三种因子为稳定因子，选取其对信噪比的最佳水平；确定三个阶段通气时间和通气温度为调整因子，利用调整因子将设计方案的质量特性调整到目标值。

[0113] 步骤八:回归拟合求得关系模型y中未知参数，并依此对调整因子进行调整，确定最佳因子水平组合，使得在该试验方案下质量特性估计值最接近目标值。

[0114] 本实施例中，利用最小二乘法，求得关系模型y中未知参数值为：

[0115] β＝(β0,β1,...,β69)

[0116] ＝(1.9666,0.0003,-0.086,-0.7368,-0.0002,-0.001,-0.0019,0.1384,0.78,[0117] -0.0083,0,0572,-0.0006,0.0028,-0.0006,0.0011,-0.1771,0.001,0.003,0.231)

[0118] 根据关系模型，调动调整因子，确定最终的最优方案为：

[0119] x*＝[30,440,0.6,120,350,0.07,0.475,120,500]

[0120] 该参数组合下，小球氮化层深度预测值为5.01mm，达到目标值要求，圆满完成了小球表面氮化试验设计的目标。