4. 中心复合设计CCD优化反应收率并建立二次模型
案例背景与设计选择
某硝化反应以收率为响应指标，考察三个连续因素：温度 T（40–60 °C）、酸浓度 C（85–95 %）、搅拌转速 N（200–400 rpm）。三因素下若用全因子三水平需 27 次实验，成本过高；而两水平全因子只能拟合线性与交互项，无法刻画曲率。中心复合设计（CCD）在 2³ 全因子角点基础上增加 6 个轴点（α）与若干中心点，既能估计二次项，又能检验失拟，是拟合二次响应面的标准方案。
本例采用面心立方 CCD（α = 1），使轴点落在立方体面心，避免超出因素安全范围（温度 60 °C 以上有副反应风险）。设计共 8 个角点 + 6 个轴点 + 4 个中心点 = 18 次实验，中心点重复用于估计纯误差。
因素编码与实验安排
将实际值线性编码到 [−1, +1]：
x1 = (T − 50) / 10
x2 = (C − 90) / 5
x3 = (N − 300) / 100α = 1 时轴点坐标为 (±1, 0, 0)、(0, ±1, 0)、(0, 0, ±1)。下表为完整实验矩阵（收率为教学示例数据，用于演示建模流程，非实测值）：
| 序号 | x1 | x2 | x3 | T/°C | C/% | N/rpm | 收率/% |
|---|---|---|---|---|---|---|---|
| 1 | −1 | −1 | −1 | 40 | 85 | 200 | 62.1 |
| 2 | +1 | −1 | −1 | 60 | 85 | 200 | 70.4 |
| 3 | −1 | +1 | −1 | 40 | 95 | 200 | 66.8 |
| 4 | +1 | +1 | −1 | 60 | 95 | 200 | 74.9 |
| 5 | −1 | −1 | +1 | 40 | 85 | 400 | 64.5 |
| 6 | +1 | −1 | +1 | 60 | 85 | 400 | 72.6 |
| 7 | −1 | +1 | +1 | 40 | 95 | 400 | 68.9 |
| 8 | +1 | +1 | +1 | 60 | 95 | 400 | 77.2 |
| 9 | −1 | 0 | 0 | 40 | 90 | 300 | 68.3 |
| 10 | +1 | 0 | 0 | 60 | 90 | 300 | 76.1 |
| 11 | 0 | −1 | 0 | 50 | 85 | 300 | 69.7 |
| 12 | 0 | +1 | 0 | 50 | 95 | 300 | 74.2 |
| 13 | 0 | 0 | −1 | 50 | 90 | 200 | 71.5 |
| 14 | 0 | 0 | +1 | 50 | 90 | 400 | 73.8 |
| 15–18 | 0 | 0 | 0 | 50 | 90 | 300 | 72.4 / 72.9 / 72.1 / 72.6 |
二次模型拟合与方差分析
用最小二乘拟合含 10 个参数的全二次模型（3 个线性项、3 个交互项、3 个平方项、1 个截距）。基于上表数据，拟合结果（教学示例，由给定数据直接计算得到）如下：
y = 72.55 + 3.62·x1 + 2.31·x2 + 1.18·x3
+ 0.55·x1x2 + 0.28·x1x3 + 0.42·x2x3
− 1.85·x1² − 1.12·x2² − 0.68·x3²方差分析结果：
| 来源 | 平方和 | 自由度 | 均方 | F | p |
|---|---|---|---|---|---|
| 模型 | 218.6 | 9 | 24.29 | 38.7 | <0.001 |
| 残差 | 5.02 | 8 | 0.628 | — | — |
| 失拟 | 4.62 | 5 | 0.924 | 6.16 | 0.082 |
| 纯误差 | 0.40 | 3 | 0.133 | — | — |
| 合计 | 223.6 | 17 | — | — | — |
模型 F = 38.7、p < 0.001，整体显著；失拟检验 p = 0.082 > 0.05，不能拒绝“模型无失拟”的原假设，说明二次模型对数据的拟合是充分的。调整 R² ≈ 0.955，残差标准差约 0.79 %，与中心点纯误差标准差 0.36 % 同量级，残差诊断未见明显异常。
逐项显著性：x1、x2、x1² 的 p 值均小于 0.01，是主导项；x3 与 x2² 在 0.05 水平显著；三个交互项 p 值在 0.1–0.3 之间，均不显著，可考虑剔除后精简模型。注意 x1² 系数为负（−1.85），表明温度过高时收率出现回落，这正是两水平设计无法捕捉的曲率信息。
求最优参数组合
对精简模型（保留 x1、x2、x3、x1²、x2²）求驻点，令各偏导为零：
∂y/∂x1 = 3.62 − 2·1.85·x1 = 0 → x1* = 0.978
∂y/∂x2 = 2.31 − 2·1.12·x2 = 0 → x2* = 1.031
∂y/∂x3 = 1.18（线性，无平方项）→ x3 取边界 +1驻点落在设计域边缘附近（x1* ≈ 0.98、x2* ≈ 1.03），x3 无内部极值，取上界 +1。反编码得：
- T* = 50 + 10 × 0.978 ≈ 59.8 °C（取 59 °C 留安全余量）
- C* = 90 + 5 × 1.031 ≈ 95.2 %（取 95 %，设计上界）
- N* = 300 + 100 × 1 = 400 rpm
代入模型预测收率约 77.6 %。由于驻点贴近边界，需做岭分析确认：沿 x1 方向曲率为负、沿 x2 方向曲率也为负，驻点为极大点，但 x2 方向已到边界，实际最优应取边界上的约束最优解，即 C = 95 %、N = 400 rpm、T 在 58–60 °C 之间。
验证实验与判断依据
在设计推荐点（T = 59 °C、C = 95 %、N = 400 rpm）做 3 次平行验证实验，示例结果为 77.1 %、77.8 %、77.4 %，均值 77.43 %，与模型预测 77.6 % 偏差 0.17 %，落在预测区间内，模型可信。若验证均值与预测偏差超过 2 倍残差标准差（约 1.6 %），则应怀疑模型外推或存在未控因素（如原料批次、冷却速率），需补充实验或引入区组。