Python实现单形重心法混料设计

介绍及概念

单纯形重心设计法是一种效率较高的混料试验统计模型，可用来设计试验，它能根据试验点和响应值给出响应曲面（模型）。该法具有较高的精度和可靠性，以及较少的试验量，很适合混料试验设计。

概念

单形（单纯形）：为空间中顶点比维数多1的凸多边形，比如〇维的点、一维的线段，二维的三角形、三维的四面体等。
响应值：某一试验点试验后的评价。
三角坐标系：由三角形构成的坐标系，顶点代表对应成分的最大值。特点是三角形内（含边）一点 $P(ax,by,cz)$，过P作三条对三角形边的平行线，可以确定x,y,z的比例。如图所示。

$单形重心法示意图$

等值线：将某一指标中数值相同的点连线。等值线实际上是将空间（Z轴）投影到平面（XY轴）的表示。如下图。

$等值线示意图$

单形格子设计

设试验中考察的指标为 $y$ ，那么 $y$ 与 $p$ 个因子 $x_1,x_2,...,x_p$ 的关系可以表示为：$y=f(x_1,x_2,\dot,x_p )+\varepsilon$

其中， $\varepsilon$ 是随机误差，且要保证 $0\leq x\leq 1 ,\sum_{i=1}^p x_i=1$

称 $y=f(x_1,x_2,\dots,x_p)$ 为响应函数，当响应函数中的未知参数用估计值代替后便得到回归方程，也称响应曲面方程。

由于 $f(x_1,x_2,...,x_p)$ 形式往往是未知的，通常用 $x_1,x_2,\dots,x_p$ 的一个 $d$ 次多项式表示，此时一个混料试验由因子数 $p$ 与响应多项式的次数 $d$ 来确定，以后用 $M\{p,d\}$ 表示一个混料试验。即：（公式1）

\begin{align*} y&=f(x_1,x_2,\dots,x_p) \\ &=\sum_{i=1}^p\eta_ix_i+\sum_{i\leq j}\eta_{ij}x_ix_j+\sum_{i\leq j\leq k}\eta_{ijk}x_ix_jx_k+\cdots \end{align*}

单形格子设计是 Scheffé 提出的一种混料设计。 $M\{p,d\}$ 的单形格子设计，为 $d$ 阶格子设计，它将单形的边划分成 $d$ 等份，在等分点做与其它边平行的直线，形成许多格子，故名单形格子设计。

$p=1$ ，为点； $p=2$ ，为线段； $p=3$ ，为三角形； $p=4$ ，为四面体； $p=5$ ，为超四面体……

举例： $p=3$ ，一阶、二阶和三阶单形格子设计的点分布图。

$格子设计示意图$

单形重心设计

单形格子设计的缺陷是在 $M\{p,d\}$ 单形格子设计中，当 $d>2$ 时某些混料设计中格子点的非零坐标（如 $M\{3,3\}$ 中的8、9及其对称点）并不相等，这种非对称性会使某些点对回归系数的估计产生较大的影响，为改进这一点， Scheffé 提出了一种只考虑有相等非零坐标的单形重心设计，既消除了以上缺陷，又不至于试验点数太多。

单形重心设计 [1] 的试验点为 $1$ 到 $p$ 个顶点的重心，顶点本身就是重心，两个顶点的重心是它们连线的中点，三个顶点的重心是它们组成正三角形的中心，……， $p$ 个顶点的重心就是该单形的中心。

此时给出的多项式模型称为 Scheffé单形重心设计的多项式模型 。具体如下，注意最后一项。（公式2）

\begin{align*} y&=f(x_1,x_2,\cdots,x_p) \\ &=\sum_{i=1}^p\eta_ix_i+\sum_{i\leq j}\eta_{ij}x_ix_j+\sum_{i\leq j\leq k}\eta_{ijk}x_ix_jx_k+\cdots+\eta_{12\cdots p}x_1x_2\cdots x_p \end{align*}

其中，若 $p=3$ ，由于只取重心试验，故有：

\begin{align*} y&=f(x_1,x_2,x_3) \\ &=\sum_{i=1}^3\eta_ix_i+\sum_{i\leq j}\eta_{ij}x_ix_j+\eta_{123}x_1x_2x_3 \end{align*}

分别取 $x_i=1, x_j=x_k=0;x_i=x_j=1/2,x_k=0;x_i=x_j=x_k=1/3$ ,分别代入单形重心设计的多项式，可得各系数与各响应值的关系。

各系数与各响应值的关系如下：

\begin{equation*} \begin{cases} \eta_i&=y_i \\ \eta_{ij}&=2\big(2^1y_{ij}-1^1(y_i+y_j)\big) \\ \eta_{ijk}&=3\big(3^2y_{ijk}-2^2(y_{ij}+y_{ik}+y_{jk})+1^2(y_i+y_j+y_k)\big) \end{cases} \end{equation*}

其实经过推导，可以获得一般公式的，若 $S_r=\{i_1,i_2,\cdots,i_r\}$ 为关于 $(1,2,\cdots,p)$ 的某个 $r$ 元素的子集，则:（公式3)

\begin{equation*} \eta_{S_r}=r\Big(\sum_{t=1}^r(-1)^{r-t}t^{r-1}L_t(S_r)\Big) \end{equation*}

其中，从 $S_r$ 对应的 $r$ 个分量中取 $t$ 个， $L_t(S_r)$ 表示所有 $C_r^t$ 个分量等比例混料响应的和（简单的说就是 $1$ 到 $r$ 取 $t$ 个的不重复组合）。

下面举例计算 $\eta_{1234}$ ：

在 $\eta_{1234}$ 中，

\begin{equation*} \begin{cases} r& = 4 \\ S_4& = \{i_1,i_2,\cdots,i_r\}=\{1,2,3,4\} \\ L_1(S_4)& = y_1+y_2+y_3+y_4 \\ L_2(S_4)& = y_{12}+y_{13}+y_{14}+y_{23}+y_{24}+y_{34} \\ L_3(S_4)& = y_{123}+y_{124}+y_{134}+y_{234} \\ L_4(S_4)& = y_{1234} \end{cases} \end{equation*}

所以根据公式3有：

\begin{equation*} \begin{split} \eta_{1234} & = r\Big(\sum_{t=1}^r(-1)^{r-t}t^{r-1}L_t(S_r)\Big) \\ & = 4\Big(\sum_{t=1}^4(-1)^{4-t}t^{4-1}L_t(S_4)\Big) \\ & = 4\Big((-1)^{4-1}1^{4-1}L_1(S_4)+(-1)^{4-2}2^{4-1}L_2(S_4)+(-1)^{4-3}3^{4-1}L_3(S_4)+(-1)^{4-4}4^{4-1}L_4(S_4)\Big) \\ & =-4(y_1+y_2+y_3+y_4)+32(y_{12}+y_{13}+y_{14}+y_{23}+y_{24}+y_{34})-108(y_{123}+y_{124}+y_{134}+y_{234})+256y_{1234} \end{split} \end{equation*}

再举例计算 $\eta_{23}$ ：

在 $\eta_{23}$ 中，

\begin{equation*} \begin{cases} r& = 2 \\ S_2& = \{i_1,i_2,\cdots,i_r\}=\{2,3\} \\ L_1(S_2)& = y_2+y_3\\ L_2(S_2)& = y_{23} \\ \end{cases} \end{equation*}

所以根据公式3有：

\begin{equation*} \begin{split} \eta_{23} & = r\Big(\sum_{t=1}^r(-1)^{r-t}t^{r-1}L_t(S_r)\Big) \\ & = 2\Big(\sum_{t=1}^2(-1)^{2-t}t^{2-1}L_t(S_2)\Big) \\ & = 2\Big((-1)^{2-1}1^{2-1}L_1(S_2)+(-1)^{2-2}2^{2-1}L_2(S_2)\Big) \\ & =2\Big(-(y_2+y_3)+2y_{23}\Big) \end{split} \end{equation*}

单形重心法代码实现

先导入相应的包。

import numpy as np
from itertools import chain, combinations

由公式3可知，我们需要生成 $S_r$ ，它是一个不含空集的幂集，简单地，生成一个不含空集的 Power Set ，这里演示 n_point 为对应的 $p$ 。

n_point = ...
nums = range(n_point)
# generate a powerset except zeroset
test_points = tuple(chain.from_iterable(
    map(lambda num: combinations(nums, num + 1), nums)))

有了 $S_r$ ，我们就可以逐个生成 $L_t(S_r)$ ，及 $\eta_{S_r}$ ，然后传入各点实测的 $y$ ，即可计算出响应曲面的系数： _response_surface_coef 。

y = np.array([...])
_response_surface_coef = []
for i, test_point in enumerate(test_points):
    r = len(test_point)
    temp = 0
    for j in range(1, r + 1):
        for test_point_pos in combinations(test_point, j):
            t = len(test_point_pos)
            temp += y[test_points.index(test_point_pos)] * \
                r * (-1)**(r - t) * t**(r - 1)
    _response_surface_coef.append(temp)

有了响应曲面的系数，即可算出每个 $\begin{pmatrix}x_1,\dots,x_n\end{pmatrix}$ 的响应值。 此处输入的是x'，也就是编码矩阵。
```
prediction = _response_surface_coef.dot(
  [X.take(test_point_pos, axis=1).prod(axis=1)
  for test_point_pos in test_points]
)
```

带下界的设计

由于单形重心要求 $0\leq x_i\leq 1$ ，但在实际试验中，很难做到。如混凝土掺合料设计试验中，取煤灰或石粉等为100%时无法获取有效数据。

$带下界的单形重心法示意图1$

如上图， $a_i$ 为各成分最小含量，我们将 $\{x'_i\}$ 称为自然空间， $\{x_i\}$ 称为编码空间（实际值），要在编码空间上进行单形重心设计，必须将其映射到自然空间上（编码值）。

\begin{equation*} \begin{cases} 0\leq x'_i\leq 1 \\ \sum_{i=1}^pa_i\leq 1 \\ 0\leq x'_i\leq x_i\leq 1 \end{cases} \end{equation*}

$带下界的单形重心法示意图2$

编码值与实际值可以进行线性变换 $\{x'_i\}\leftrightarrow\{x_i\}$ ，由上图可得变换矩阵 $M$ ：

\begin{equation*} M = \begin{bmatrix} 1+a_1-\sum_{i=1}^pa_i & a_1 & a_1 & \dots & a_1 \\ a_2 & 1+a_2-\sum_{i=1}^pa_i & a_1 & \dots & a_1 \\ \vdots & \vdots & \ddots & \dots &\vdots \\ a_p & a_p & a_p & \dots & 1+a_p-\sum_{i=1}^pa_i \\ \end{bmatrix} \end{equation*}

再通过坐标与自然空间相乘得 $x_i$ （公式4)

\begin{equation*} \begin{pmatrix} x_1 \\ x_2 \\ \vdots \\ x_p \end{pmatrix} = M\begin{pmatrix} x'_1 \\ x'_2 \\ \vdots \\ x'_p \end{pmatrix} \end{equation*}

实际值与编码值也可以按以下公式进行转换（线性变换化简的公式，虽然文献中用的多，但个人觉得不够直观，不推荐使用）：

\begin{equation*} \begin{cases} x'_i&=(1-\sum_{i=1}^pa_i)x_i+a_i \\ x_i&=\frac{x'_i-a_i}{1-\sum_{i=1}^pa_i} \end{cases} \end{equation*}

编码空间上的点 $X$ 为配料的真实比例，而自然空间的点 $X'$ （编码值）则为变换后的符合单纯形设计的比例。

编码空间映射到的自然空间并不是成分均为100%的单纯形，仅仅是概念上的自然空间。故应直接按单纯形重点设计，得到数据之后通过等值线找到自然空间中的点，再通过编码转换得到真实的配比。

带下界的单形重心法代码实现

假设下界为变量 lower_bounds ， $lower\_bounds = [a_1,a_2,\dots,a_p]$ ，当 $lower\_bounds = [0,\dots,0]$ 时，即是普通的单形重心法。

由公式3生成 $S_r$

n_point = ...
nums = range(n_point)
# generate a powerset except zeroset
test_points = tuple(chain.from_iterable(
    map(lambda num: combinations(nums, num + 1), nums)))

生成变换矩阵 $M$ ，

_M = lower_bounds.repeat(n_point).reshape(
    (n_point, n_point)) \
    - np.eye(n_point) \
    - (1 - lower_bounds.sum())

如前所述，有了 $S_r$ ，然后传入各点实测的 $y$ ，即可计算出响应曲面的系数： _response_surface_coef ，注意，此处不需要转换矩阵 M 。

y = np.array([...])
_response_surface_coef = []
for i, test_point in enumerate(test_points):
    r = len(test_point)
    temp = 0
    for j in range(1, r + 1):
        for test_point_pos in combinations(test_point, j):
            t = len(test_point_pos)
            temp += y[test_points.index(test_point_pos)] * \
                r * (-1)**(r - t) * t**(r - 1)
    _response_surface_coef.append(temp)

有了响应曲面的系数，即可通过转换矩阵 _M 和真实比例 X 算出每个 $\begin{pmatrix}x'_1,\dots,x'_n\end{pmatrix}$ 的响应值。 此处输入的是X，也就是编码矩阵。
```
XX = X.dot(np.linalg.inv(_M.T))
prediction = _response_surface_coef.dot(
  [XX.take(test_point_pos, axis=1).prod(axis=1)
  for test_point_pos in test_points]
)
```

应用举例说明

这里对一种调料和一种混凝土的带下界约束单形重心试验设计进行了应用的举例。

两个例子均为3成分混料试验，由于公式2、公式3、公式4均是一般情形的公式，可以推广到任意数量成分的混料试验，不赘述。

例1，调料配制

一种调料由三种成分 $A、B、C$ 混合制成 $A、B、C$ 各为味精、盐、五香粉。 $a\geq 0.2,b\geq 0.4,c\geq0.2$ 。求设计方案。

此处可以按公式1采用 $M\{3,2\}$ 单形格子设计，本文主要讲单形重心设计，就不按格子设计来了。采用单形重心设计，取点如下图：

$示例1设计图$

根据 $a\geq 0.2,b\geq 0.4,c\geq0.2$ 画出小单形，即黑色小三角，然后根据单形重心设计标出7个重心（红色数字表示）

由公式4列出编码矩阵：

\begin{align*} M &= \begin{bmatrix} 1+a-(a+b+c) & a & a \\ b & 1+b-(a+b+c) & b \\ c & c & 1+c-(a+b+c) \\ \end{bmatrix} \\ &= \begin{bmatrix} 0.4 & 0.2 & 0.2 \\ 0.4 & 0.6 & 0.4 \\ 0.2 & 0.2 & 0.4 \\ \end{bmatrix} \end{align*}

根据公式2、公式3和 $X'*(Z.T)$ 可轻松列出试验表和试验结果以及口感得分（A*、B*、C*为编码值，味精_ 、盐_ 、五香粉_ 则为实际成分）

试验号	A*	B*	C*	味精_	盐_	五香粉_	口感
	$x'_1$	$x'_2$	$x'_3$	$x_1$	$x_2$	$x_3$
$y_1$	1	0	0	0.4	0.4	0.2	5
$y_2$	0	1	0	0.2	0.6	0.2	11
$y_3$	0	0	1	0.2	0.4	0.4	8
$y_{12}$	1/2	1/2	0	0.3	0.5	0.2	10
$y_{13}$	1/2	0	1/2	0.3	0.4	0.3	2
$y_{23}$	0	1/2	1/2	0.2	0.5	0.3	10
$y_{123}$	1/3	1/3	1/3	4/15	7/15	4/15	13

将结果和编码值代入公式2得：

\begin{equation*} y=159x'_1x'_2x'_3+8x'_1x'_2-18x'_1x'_3+5x'_1+2x'_2x'_3+11x'_2+8x'_3 \end{equation*}

$示例1结果图$

作三角坐标图，根据三角坐标系在最大值作图，可得自然空间中最大值坐标 ${x'_i}$ 为：

\begin{equation*} \begin{pmatrix} x'_1 \\ x'_2 \\ x'_3\end{pmatrix}=\begin{pmatrix}0.26 \\ 0.48 \\ 0.26 \end{pmatrix} \end{equation*}

用编码矩阵转换得真实比例为：

\begin{equation*} \begin{pmatrix}0.252,0.496,0.252\end{pmatrix} \end{equation*}

即，按这个比例配制的调料味道最好。

例2，混凝土强度预测

$示例2设计图$

如上图，混凝土用胶凝材料为水泥，矿粉，煤灰，其中水泥用量在25%以上，求单形重心试验方案。

这是一个约束设计问题，水泥用量25%以上，即 $a_1\geq 0.25, a_2=a_3=0$ 。编码矩阵 $M$ 为：

\begin{equation*} \begin{bmatrix} 1 & 0.25 & 0.25 \\ 0 & 0.75 & 0 \\ 0 & 0 & 0.75 \end{bmatrix} \end{equation*}

标上7个实验点，列出试验表并根据表来做实验得到结果 [2] ：

试验号	A*	B*	C*	水泥_	矿粉_	煤灰_		强度(MPa)
	$x'_1$	$x'_2$	$x'_3$	$x_1$	$x_2$	$x_3$	3d	28d	180d
$y_1$	1	0	0	1	0	0	63.1	88.3	96
$y_2$	0	1	0	0.25	0.75	0	29.0	56.2	77
$y_3$	0	0	1	0.25	0	0.75	22.2	53.5	75.4
$y_{12}$	1/2	1/2	0	0.625	0.375	0	50.6	84.5	90.1
$y_{13}$	1/2	0	1/2	0.625	0	0.375	44.5	92.3	102
$y_{23}$	0	1/2	1/2	0.25	0.375	0.375	26.5	62.8	86
$y_{123}$	1/3	1/3	1/3	0.5	0.25	0.25	40.3	80.5	96.5

将结果和编码值代入公式2得：

\begin{align*} y_{3d}&=63.1x'_1+29.0x'_2+22.2x'_3+18.2x'_1 x'_2+7.4x'_1 x'_3+3.6x'_2 x'_3-28.2x'_1 x'_2 x'_3 \\ y_{28d}&=88.3x'_1+56.2x'_2+53.5x'_3+49x'_1 x'_2+85.6x'_1 x'_3+31.8x'_2 x'_3-107.7x'_1 x'_2 x'_3 \\ y_{180d}&=96x'_1+77x'_2+75.4x'_3+14.9x'_1 x'_2+65.2x'_1 x'_3+39.2x'_2 x'_3+13.5x'_1 x'_2 x'_3 \end{align*}

对上面三式作三角坐标图，可以清晰地看出各组分对强度的贡献。从图中求出 $x'_1$ , $x'_2$ , $x'_3$ 坐标，再利用编码矩阵即可换算出实际各组分比例。另外，从3d与28d、180d对比，说明前期是水泥、矿粉对强度贡献大，后期煤灰贡献逐渐超过矿粉。

作图

对组分的单形重心设计可以作图。

$示例2结果图$

[1]	关颖男. 混料试验设计. 上海科学技术出版社, 1990.

[2]	孙伟, 严捍东. 复合胶凝材料组成与混凝土抗压强度定量关系研究[J]. 東南大學學報 (自然科學版), 2003, 33(4): 450-453.

Read more:

≪ dnspod DDNS 自动更新用 Docker 编译 Openwrt ≫

Python数据分析-安装与配置

试验号	A*	B*	C*	味精_	盐_	五香粉_	口感
	\(x'_1\)	\(x'_2\)	\(x'_3\)	\(x_1\)	\(x_2\)	\(x_3\)
\(y_1\)	1	0	0	0.4	0.4	0.2	5
\(y_2\)	0	1	0	0.2	0.6	0.2	11
\(y_3\)	0	0	1	0.2	0.4	0.4	8
\(y_{12}\)	1/2	1/2	0	0.3	0.5	0.2	10
\(y_{13}\)	1/2	0	1/2	0.3	0.4	0.3	2
\(y_{23}\)	0	1/2	1/2	0.2	0.5	0.3	10
\(y_{123}\)	1/3	1/3	1/3	4/15	7/15	4/15	13