3. 双目标DOE的因子水平编码与试验设计表生成
四个因子的水平编码
本例的配方场景沿用前两例已确定的因子范围：增韧剂含量（A）、稀释剂含量（B）、固化温度（C）、促进剂用量（D）。响应曲面设计需要把物理量换算成无量纲编码值，中心点记为 0，低水平记为 −1，高水平记为 +1，星点（轴向点）记为 ±α。编码公式为：
x = (X − X₀) / (ΔX)，其中 X₀ 为水平中心值，ΔX 为半区间（高水平与低水平之差的一半）。
| 因子 | 物理量 | 低水平 (−1) | 中心 (0) | 高水平 (+1) | ΔX | 单位 |
|---|---|---|---|---|---|---|
| A | 增韧剂含量 | 8 | 12 | 16 | 4 | phr |
| B | 稀释剂含量 | 3 | 6 | 9 | 3 | phr |
| C | 固化温度 | 60 | 80 | 100 | 20 | °C |
| D | 促进剂用量 | 0.5 | 1.0 | 1.5 | 0.5 | phr |
取 α = 2 的可旋转中心复合设计（CCD），每个因子有 5 个水平：−2、−1、0、+1、+2。星点对应的物理量为 X₀ ± 2ΔX，例如增韧剂含量星点为 4 phr 与 20 phr，稀释剂为 0 phr 与 12 phr，固化温度为 40 °C 与 120 °C，促进剂为 0 phr 与 2.0 phr。
24 run 中心复合设计的结构
标准 CCD 由三部分组成：2⁴ = 16 个析因点（立方点，编码值取 ±1）、8 个星点（轴向点，每个因子各取 ±α，其余因子为 0）、以及若干中心点。16 + 8 = 24，正好满足 24 run 的要求，因此本例采用零中心点的 CCD 配置。若需要估计纯误差，可把部分析因点替换为中心点，但会损失部分高阶交互信息，本例优先保证设计表规模与可执行性。
生成 24 run 设计表
下表按“析因点 → 星点”顺序列出编码值，物理量由编码公式反算：X = X₀ + x·ΔX。表中同时给出物理配方，便于直接称量配料。
| Run | 类型 | A 编码 | B 编码 | C 编码 | D 编码 | A (phr) | B (phr) | C (°C) | D (phr) |
|---|---|---|---|---|---|---|---|---|---|
| 1 | 析因 | −1 | −1 | −1 | −1 | 8 | 3 | 60 | 0.5 |
| 2 | 析因 | +1 | −1 | −1 | −1 | 16 | 3 | 60 | 0.5 |
| 3 | 析因 | −1 | +1 | −1 | −1 | 8 | 9 | 60 | 0.5 |
| 4 | 析因 | +1 | +1 | −1 | −1 | 16 | 9 | 60 | 0.5 |
| 5 | 析因 | −1 | −1 | +1 | −1 | 8 | 3 | 100 | 0.5 |
| 6 | 析因 | +1 | −1 | +1 | −1 | 16 | 3 | 100 | 0.5 |
| 7 | 析因 | −1 | +1 | +1 | −1 | 8 | 9 | 100 | 0.5 |
| 8 | 析因 | +1 | +1 | +1 | −1 | 16 | 9 | 100 | 0.5 |
| 9 | 析因 | −1 | −1 | −1 | +1 | 8 | 3 | 60 | 1.5 |
| 10 | 析因 | +1 | −1 | −1 | +1 | 16 | 3 | 60 | 1.5 |
| 11 | 析因 | −1 | +1 | −1 | +1 | 8 | 9 | 60 | 1.5 |
| 12 | 析因 | +1 | +1 | −1 | +1 | 16 | 9 | 60 | 1.5 |
| 13 | 析因 | −1 | −1 | +1 | +1 | 8 | 3 | 100 | 1.5 |
| 14 | 析因 | +1 | −1 | +1 | +1 | 16 | 3 | 100 | 1.5 |
| 15 | 析因 | −1 | +1 | +1 | +1 | 8 | 9 | 100 | 1.5 |
| 16 | 析因 | +1 | +1 | +1 | +1 | 16 | 9 | 100 | 1.5 |
| 17 | 星点 | −2 | 0 | 0 | 0 | 4 | 6 | 80 | 1.0 |
| 18 | 星点 | +2 | 0 | 0 | 0 | 20 | 6 | 80 | 1.0 |
| 19 | 星点 | 0 | −2 | 0 | 0 | 12 | 0 | 80 | 1.0 |
| 20 | 星点 | 0 | +2 | 0 | 0 | 12 | 12 | 80 | 1.0 |
| 21 | 星点 | 0 | 0 | −2 | 0 | 12 | 6 | 40 | 1.0 |
| 22 | 星点 | 0 | 0 | +2 | 0 | 12 | 6 | 120 | 1.0 |
| 23 | 星点 | 0 | 0 | 0 | −2 | 12 | 6 | 80 | 0 |
| 24 | 星点 | 0 | 0 | 0 | +2 | 12 | 6 | 80 | 2.0 |
设计点可执行性检查
编码表生成后必须逐点核对物理可行性，否则会出现“设计上合法、实验上做不出来”的点。本例需要重点确认以下四项：
- 稀释剂 0 phr（Run 19）：无稀释剂时体系黏度可能过高，需确认能否在 80 °C 下搅拌并浇注。若不可行，应把稀释剂下限从 3 phr 上调，或改用 α 较小的设计。
- 促进剂 0 phr（Run 23）：无促进剂时 80 °C 是否能在合理时间内凝胶，需做一次预试验确认，否则该点响应可能无法测量。
- 固化温度 120 °C（Run 22）：确认烘箱控温上限与模具/基材耐温，避免超温导致试样报废。
- 增韧剂 20 phr（Run 18）：高增韧剂含量下混合均匀性与黏度是否仍在可操作范围，需确认搅拌设备扭矩余量。
用代码生成并校验设计表
下面用 Python 生成编码矩阵、反算物理量，并做基本校验（析因点数量、星点数量、编码值范围、物理量是否落在设定区间内）。代码自包含，可直接运行。
import itertools
import numpy as np
# 因子定义：名称, 中心值, 半区间, 单位
factors = [
("A_增韧剂", 12.0, 4.0, "phr"),
("B_稀释剂", 6.0, 3.0, "phr"),
("C_固化温度", 80.0, 20.0, "°C"),
("D_促进剂", 1.0, 0.5, "phr"),
]
alpha = 2.0
# 1) 析因点：2^4 全因子，编码 ±1
factorial = list(itertools.product([-1, 1], repeat=4))
# 2) 星点：每个因子 ±alpha，其余为 0
star = []
for i in range(4):
for s in (-alpha, alpha):
row = [0.0] * 4
row[i] = s
star.append(tuple(row))
# 3) 合并
coded = factorial + star
assert len(coded) == 24, f"run 数应为 24，实际 {len(coded)}"
# 4) 反算物理量
def decode(row):
return [c + x * d for (_, c, d, _), x in zip(factors, row)]
physical = np.array([decode(r) for r in coded])
# 5) 校验：析因点 16 个、星点 8 个
n_fact = sum(1 for r in coded if all(abs(v) == 1 for v in r))
n_star = sum(1 for r in coded if sum(abs(v) for v in r) == alpha)
print(f"析因点 {n_fact} 个，星点 {n_star} 个，合计 {len(coded)} run")
# 6) 校验：物理量是否落在 [中心-2*半区间, 中心+2*半区间]
for j, (name, c, d, unit) in enumerate(factors):
lo, hi = physical[:, j].min(), physical[:, j].max()
print(f"{name}: {lo:.2f} ~ {hi:.2f} {unit}")
# 7) 输出设计表
header = ["Run", "A", "B", "C", "D", "A_phr", "B_phr", "C_C", "D_phr"]
print(",".join(header))
for i, (cr, pr) in enumerate(zip(coded, physical), start=1):
print(f"{i}," + ",".join(f"{v:+.0f}" for v in cr) + "," +
",".join(f"{v:.1f}" for v in pr))
运行后应看到：析因点 16 个、星点 8 个、合计 24 run；四个因子的物理量范围分别为 4–20 phr、0–12 phr、40–120 °C、0–2.0 phr，与编码表一致。若某因子范围与预期不符，说明中心值或半区间录入有误，应先修正再进入试验。