混料试验设计.doc

王中来：混料试验设计混料试验设计The Design of Mixture Experiments主要参考文献：1、 栾军. 现代试验设计优化方法. 上海：上海交通大学出版社，19952、 茆诗松等. 回归分析及其试验设计. 上海：华东师范大学出版社, 1981一、 混料问题与混料试验 (栾军, 1995；茆诗松 等, 1981)日常生活中和工业生产上经常遇到配方配比一类的问题，即所谓混料问题。这里所说的混料是指由若干不同成分的元素混合形成一种新的物品。由不同成分组成的钢、铁、铝、药方、饲料以及燃料等都是混料，某些分配问题，如企业的材料、资金、设备和人员等的分配也可看着混料问题。混料试验就是通过实物试验或非实物试验，考察各种混料成分与试验指标之间的关系。例如，人们吃的糕点是将面粉、水、油、糖发酵及某些香料混合后经烘烤制成的，考察这些成分对糕点的柔软性、口味等试验指标的影响所进行的试验就是混料试验。应该指出，混料试验中的混料成分至少应有三种，并且混料成分中的不变成分不应作为混料成分。混料试验设计，不同于以前所介绍的各种试验设计。混料试验设计的试验指标只与每种成分的含量有关，而与混料的总量无关，且每种成分的比例必须是非负的，且在01之间变化，各种成分的含量之和必须等于1（即100%）。也就是说，各种成分不能完全自由地变化，受到一定条件的约束。设：y为试验指标，x是第i种成分的含量，则混料问题的约束条件，即混料条件为： （1）其中xi称为混料成分或混料分量，即混料试验中的试验因素。混料试验设计是一种受特殊条件约束的回归设计，它是通过合理地安排混料试验，以求得各种线性或非线性回归方程的技术方法。它具有试验点数少、计算简便、容易分析、迅速得到最佳混料条件等优点。混料条件（1）决定了混料试验设计不能采用一般多项式作为回归模型，否则会由于混料条件的约束而引起信息矩阵的退化。混料试验设计常采用Scheff 多项式回归模型。例如，一般的三元二次回归方程为(2)而混料试验设计中，三分量二次回归方程应为(3)比较式（2）和式（3）可知，Scheff多项式没有常数项和平方项。这是因为，将约束条件代入式（2），即可推导得到式（3）。通常，混料试验设计的p分量d次多项式回归方程，其Scheff多项式（或称为规范多项式）为一次式（d=1）：(4)二次式（d=2）：(5)三次式（d=3）：(6)式中为三次项的回归系数。由此看来，混料试验设计的 (p, d) Scheff 多项式回归方程中，待估计的回归系数的个数，比一般的p因素d次多项式回归方程要少。例如，对于混料试验设计（p, d）的回归方程式（5），无常数项和二次项。于是，减少了 p+1 个回归系数，所以至少可以少做 p+1 次试验。混料试验设计由 H. Scheff 于1958年首先提出，至今已有40多年。由于这种试验设计方法与工农业生产及科学试验有密切的关系，所以无论在理论研究中，还是实际应用中都有了很大的发展。在工业试验方面，合金、混凝土、陶瓷、油漆、混纺纤维、医药、食品等的配方和生产制造都广泛地应用混料试验设计方法。二、单（纯）形格子设计 （茆诗松等，1981；栾军，1995）1. 引言 （1）单（纯）形在混料试验设计方法中，单纯形格子设计是最早出现的，是Scheff 于1958年提出的。它是混料试验设计中最基本的方法，其它一些方法都要用到单纯形格子设计。在混料问题中，各分量 xi (i=1, 2, , p) 的变化范围受混料条件式（1）的制约。在几何上，称为p维平面，而（x1, x2, , xp）为p维平面上点的坐标。在p维平面上满足 的区域构成一个图形称为单形（或单纯形）。单形上的点，若其p个坐标中有一个坐标 xi=1 , 而其余的 p-1 个坐标 xj=0 (ji), 则这种点称为单形的顶点。因此，在p因子混料试验中，单形的顶点有p个。例如，p=3时，单形的三个顶点为（1，0，0）、（0，1，0）和（0，0，1）。所以单型的图形为一等边三角形，如图1（a）所示。（2）单形上点的坐标下面，以p=3为例讨论单形上点的坐标问题。对于三因子混料试验，这个试验的单形是一个等边三角形，其三个顶点分别为A（1，0，0）、B（0，1，0）和C（0，0，1）。x1设P（x1、x2、x3）为这单形的内点，定义x1表示P点到边BC的距离，x2为P点到边AC的距离，x3为P点到边AB的距离。为简单起见，使用时不再画出三个坐标轴，只画出一个等边（正）三角形，如图1（b）所示。A(1, 0, 0)A(1, 0, 0)x3x2P(x1, x2, x3)ox1x3x2C(0, 0, 1)B(0, 1, 0)B(0, 1, 0)C(0, 0, 1)(a) (b)图1 p=3时的单形（x1+x2+x3 = 正三角形的高 = 1）取此等边（正）三角形的高为1，则由初等几何学可知，ABC内任一点P到三个边的距离之和为1，即 或 所以，三因子混料试验可以用等边三角形这样一个单形上的点表示。一般情况下，对p因子混料试验，其p个顶点分别为A1（1，0，0，0）、A2（0，1，0，0）、 Ap（0，0，0，1）。设P（x1, x2, , xp）为单形的内点，定义x1，x2，xp分别表示P 点到A2Ap面的距离、A1A3Ap面的距离、A1Ap-1面的距离，并取p-1维空间内正多边形的高为1正多边形又称为正规单（纯）形。于是，我们就建立了p因子混料试验的单形坐标系。2. 单（纯）形格子点的概念对于由混料条件式（1）构成的正规单（纯）形因素空间，当采用式（5）、式（6）等完全型规范多项式回归模型时，试验点可以取在正规单（纯）形格子点上，构成单（纯）形格子设计。对于三因素（p=3时）的格子点集，其单（纯）形是一个高为1的等边三角形，它的三个顶点的全体称为一阶格子点集，记为3，1，如图2（a）所示。 31x1 = 1x1 = 1x1 = 1x3 = 1x3 = 1x3 = 1x2 = 1x2 = 1x2 = 1 (a) 3, 1(b) 3, 2 (c) 3, 3图2 单（纯）形格子点分布图 将高为1的等边三角形的三条边各二等分，则该三角形的三个顶点与三条边的中点的全体称为二阶格子点集，记为3，2，如图2（b）所示。其中共有6个点，各点坐标如表1所示。将等边三角形的各边进行三等分，对应分点连成与各边平行的直线，在等边三角形上形成许多格子，则这些格子顶点的全体称为三阶格子点集，记为3，3，如图2（c）所示，其中共有10个点。各点坐标如表2所示。表1 3，2各点坐标（三因素二阶格子点集坐标）点号坐 标x1x2x3 110020103001405060表2 3，3各点坐标（三因素三阶格子点集坐标）点号坐 标x1x2x311002010300140506070809010四因素二阶格子点集，记为4，2，其中共有10个点，各点坐标见表3。表3 4，2各点坐标（四因素二阶格子点集坐标）点号坐 标x1x2x3x4110002010030010400015006007008009001000一般情况下，p因素d阶格子点集记为p，d，其中p表示单（纯）形顶点的个数，d表示单（纯）形边长被等分的段数，并且总的点数为。3、单（纯）形格子设计法Scheff提出的单（纯）形格子设计，具有以下两个特点：（1）、每个p，d设计所要做的试验次数为，恰好等于完全型规范多项式回归方程如式（5）、式（6）所示中的回归系数的个数。因而单（纯）形格子设计是饱和设计，是一种优化设计。代表试验的点对称地排列在单形上，构成单形的一个格子，称为p，d格子。每一点的p个坐标代表p个因素的成分值，它们加起来的和等于1；（2）、试验点的成分与模型的次数（或阶数）d有关，我们约定每一成分xi取值为的倍数，即。并且在设计中因素成分量的各种配合都使用到。单（纯）形格子设计的试验点数，与相应的完全型规范多项式回归方程的阶数（或次数）d的关系，见表4。表4 单（纯）形格子设计的试验点数（）Pd234361015410203551535706215612683612033010552207154、回归系数的计算在单（纯）形格子设计中，每个回归系数的值只取决于所对应的一些格子上的观测值，而与其它设计点上的观测值无关，故使得用最小二乘法计算回归系数变得很简单。各回归系数均可表达为相应设计点上观测值的简单线性组合。例1 3，2单形格子设计 是三因素二次单形格子设计，即 p=3，d=2， 响应方程（又称Scheff多项式或规范多项式或正则多项式）为： 成分的取值0，1。此时单形格子设计及试验结果如表5所示。表5 单形格子3, 2设计及试验结果试验号x1x2x3观测值1100y12010y23001y340y450y560y6为了求出模型中系数的估计，可以通过把同一号的试验成分值分别代入上述模型中的x1、x2和x3，并把观测值yi代入对应的，这样便得到一组方程，解这组方程便可获得回归系数的估计值。将第1号试验（1，0，0）代入（x1，x2，x3），便得到 同理可得 ， 例2 4，2单形格子设计 p=4，d=2 未知参数的个数单形格子设计的试验点个数也是10。成分的取值xi=0，1。此时的单形格子设计如表3所示，并将10个试验观测值放在最后一列。4，2单形格子设计的规范多项式（即响应方程）为： 与前面一样，可得估计参数如下：， ， ， ， ， ， 例3 3，3单形格子设计 p=3，d=3， 未知参数的个数=试验次数也是10 ，成分的取值0，1。响应方程为 此时的单形格子设计，如表2所示。将10个试验观测值填入表中最后一列。同理可得待估系数如下：， ， ， ， ，单形格子设计的很大优点是：响应函数（即规范多项式）各项的意义很容易理解。在上面的系数估计公式中可以看出，一次项的系数bi只受单形的第i个顶点上观测值的影响。一次项bixi反映了纯组分（只有一个坐标xi为1，其余均为零）的响应。若响应曲面是一个平面，则方程 （7）表示成分与观测值（响应值）的线性组合（即xi之间的线性组合）。二次项系数bij仅受单形中连接第i个和第j个顶点的棱上试验点观测值的影响。全体二次项表示响应曲面与式（7）所表示的平面之间的离差。应当指出的是，二次项不能单纯理解为和的交互效应，这是因为它们受约束条件式（1）的限制，不能独立地变动，所以它们只表示一种非线性混合的关系。当0 时，Scheff称这种非线性混合关系为协调的；而当2时，某些混料试验中格子点的非零坐标不相等，这种非对称性反映到估计响应函数（即回归多项式或规范多项式）的系数时，就会出现某些观测值对回归方程影响大，而某些观测值对回归方程影响小。此外，单形格子设计的试验次数还是比较多（用计算）。为了改进上述两个缺点，Scheff（1958）提出单（纯）形重心设计，对单（纯）形格子设计进行改进，使混料试验中格子点的非零坐标相等。在一个p因子单形重心设计中，试验点数目是，包括：p个顶点（1，0，0），（0，0，1），共有个点；两个顶点的重心点（），（），共有个点；三个顶点的重心点（），（），共有个点；p个顶点的重心点（），共有个点。显然，单形重心设计的全部点的坐标不依赖于d。当p=3时，单形重心设计的试验点数目=，包括：以（1，0，0）为代表的个点；以（）为代表的个点；以（）为代表的个点。所以，总的试验点数目为.试验点分布和试验方案，见图3和表7。x1=1 x2=1x3=1图3 试验点分布图表 7 试验方案试验点x1x2x3y1100y12010y23001y341/21/20y1251/201/2y13601/21/2y2371/31/31/3y123显然，当p=3时，有(10)例5 3，3单形重心设计三元三次回归方程（p=3和d=3时），即单形重心设计的回归方程为： （11）上式中回归系数的计算公式为（回归系数的推导与单形格子设计相同） （12）一般情况下，p因素单形重心设计的总试验次数为 （13）通常，对单形重心设计，各个回归系数的计算通式为： （14）式中 Dr p个成分中某r个的集合； yt(Dr) r个成分中取出t个的全部C个组合的试验指标值的总和。（注：如果对上式通式不理解，那么也可象单形格子设计一样，针对具体情况，自己推导具体的计算公式）例6 4，3单形重心设计当p=4和d=3时，即4, 3单形重心设计的回归方程为 （15）上式中各回归系数，按式（14）计算如下：r=1时，则t=1, r=2时，则t=2和1,r=3时，则t=3、2和1， ，()回归系数的总数，也就是要做15次试验，而4，3单形格子设计的试验次数为20次（）。显然，单形重心设计的试验次数，比单形格子设计的要少。4, 3单形重心设计的设计表，见表8。表8 4，3单形重心试验设计表试验号X1X2X3X4测试值y11000201003001040001500600700800900100011012013014015表8中的第15号试验可以省去不做，因为4，3单形重心设计是四元三次模型，即式（15）中已设有一项。一般情况下，当我们考虑p，d单形重心试验设计（p元d次）时，若d p，截尾模型时可以省去单形重心设计中后面的一些设计点。通常，p，d单形重心设计的响应曲面方程为： （16） 2、单（纯）形重心设计实例例7 用单形重心设计检验两种评分法无差异的假设是否成立RM评分法是在试验室里对燃料抗震性能使用的一种评分方法，而实际中采用的是道路行驶性能评分法。今欲设计一组试验，系统地改变燃料的特性来检验这两种评分方法无差异的假设是否成立。本试验考核的指标yi 为两种评分法之差值。试验燃料的三个成分为：x1 石蜡环烷；x2 二芳香烃；x3 二烯烃。试验的约束条件为：（17）本试验选用3，3单形重心设计，试验次数为：试验方案及结果如表9所示。表9 3，3单形重心设计试验方案及结果试验号x1x2x3y1100y1 = 4.62010y2 = 4.93001y3 = 0.841/21/20y12 = 4.851/201/2y13 = 3.8601/21/2y23 = 3.071/31/31/3y123 = 3.7根据式（12）或式（14），可求出回归系数：将回归系数代入方程式（11），得 （18）所以（19）讨论：（1） 如果根据实际情况认为响应值相差半个单位（即y=0.5）是显著的，那么只含x1 的石蜡环烷的燃料和只含x2 的二芳香烃的燃料有几乎相等的响应值，因测试值中y1 与y2 的差小于0.5（y1 = 4.6x1 = 4.61 = 4.6, y2 = 4.9x2 = 4.91 = 4.9, ）。但它们都大大地超过只含x3 的二烯烃燃料的响应值y3 (y3 = 0.8 x3 = 0.81 =0.8 0, 所以存在着一种协调的非线性混合关系，b13 x1 x3 项的极大值 = ，所以这一项对响应值的影响很大。（3）由于 b123 = -8.4 3000 的等高线全部落在坐标系内部。为满足试验要求，还应求出在下列条件约束下的适宜配料比。24502650255027502850300029503050图5 式（28）的等高线图（示意图）利用式（28）画出指标估计值的等高线（见图5中的红色点），估计出较好的混料点为：于是，实际成分的混料点为：即较好的喷气剂配方为：粘合剂21.0%，氧化剂48.2%，燃料30.8%，用这一配方的喷气剂，实测的弹性模量为3010，并且粘合剂用量接近下界。 五、极端顶点设计（栾军，1995；茆诗松等，1981）有些混料问题，常常会同时兼有上、下界约束条件的限制，其混料条件为： （29）这种问题有几种设计方法，这里仅介绍一种简便的极端顶点设计法。我们称在极限平面和的等交线上满足的点为极端顶点。利用极端顶点子集所构成的混料试验设计称为极端顶点设计。如果响应函数的未知参数个数多于极端顶点个数的话，那么需要补充一些由极端顶点构成的棱、面、体的中心作为试验点，这样得到得设计也称极端顶点设计。换言之，对于单（纯）形坐标系，满足式（29）约束条件得点的总体，就是（p-1）维正单（纯）形上的一个（p-1）维凸面体。该多面体的（p-2）维边界面分别与相应坐标面（）平行。极端顶点设计就是把试验点取在该凸多面体的顶点及各个（p-2）维边界面的重心上，或者再加上各顶点的重心，构成兼有上、下界约束的混料问题的极端顶点设计。有关极端顶点设计的具体方法，一般都是通过一个实例进行讲解，极端顶点的构造分两步进行：第一步是找出以上界或下界为成分值的极端顶点，第二步是找出全部边界面重心试验点和整个多面体的重心坐标。然后，在极端顶点上进行试验，测得指标值，并用最小二乘法求得多项式回归方程。最后，通过求极值得到最佳配方。请参见茆诗松等人编著的回归分析及其试验设计一书（1981）中的例7（p.320）。录入：96食工 肖瑶煌、张映斌、翁小建、郭晓玲校对：王中来2000年5月19日星期五2004年3月23日星期二