MATLAB实现LHS采样与二阶响应面：从代理建模到多目标优化实战

1. 整体方案设计与思路拆解

先说说我为什么觉得这套组合值得专门写一篇。做结构优化或者参数优化的人应该都有这种体验：设计变量十几个，目标函数来自有限元仿真或者 CFD 计算，一次仿真少则几分钟、多则几个小时。想直接在里面套遗传算法，理论上要迭代上千次，算下来几个月都跑不完。更麻烦的是，很多时候仿真软件根本不给你解析梯度，传统梯度优化算法都没法直接用。

这种情况下最实用的思路就是把问题拆成三段式流水线。第一步用拉丁超立方采样在变量空间里布点，代替盲目试算，用尽量少的点覆盖尽量大的范围。第二步用这些采样点去跑仿真，得到一组输入输出数据，然后用二阶多项式回归去拟合它们，建一个响应面模型。第三步真正想要的寻优过程，不管是找工程上的最优解还是做多目标权衡，都在这棵响应面上进行，而不是回到昂贵的仿真里去搜。这也就是标题里“LHS 采样 + 二阶响应面 + 优化求解”要做的事，MATLAB 恰好能把这三个环节全部串起来。

1.1 为什么把优化拆成“抽样—建模—寻优”三步

传统的直接优化方法，看起来逻辑很简单，就是优化算法反复调用仿真程序，用当前设计点算目标函数值，然后迭代更新设计变量。问题是仿真程序本身是黑盒，没有梯度表达式，算法只能靠有限差分去近似梯度，或者靠启发式搜索。有限差分意味着每个变量至少要额外跑一两次仿真，启发式搜索意味着动辄几百上千次评估。两条路落在仿真代价高的场景里都是死局。

三步式优化真正聪明的地方在于，它把一个需要几千次仿真的搜索过程，压缩成一次几十次仿真的抽样过程加一个毫秒级的计算过程。抽样阶段只调用几十次仿真，建立响应面，之后不管是用非线性规划还是遗传算法，都是在解析表达式上计算，速度完全不可同日而语。对于工程实际问题，这个思路特别适合那种设计空间不大、变量个数适中、目标函数虽然计算耗时但整体趋势相对光滑的场景。

1.2 为什么选 LHS 而不是随机抽样或全因子设计

实验设计方法里，全因子设计最直接，但 k 个变量每个取 m 个水平就需要跑 m 的 k 次方次仿真。变量稍微一多，组合就爆炸了。纯随机抽样虽然次数可以自己定，但小样本下很容易出现变量空间里的点扎堆，或者在某个角落完全空白，导致后面的回归模型存在严重偏差。

拉丁超立方采样是一种分层随机抽样方式，它先把每个变量的取值范围均匀地分成 n 个区间，然后在每个区间里随机取一个样本值，最后再把各维度的值随机配对组合成 n 个样本点。这样做的好处是每个变量的一维投影都完全均匀覆盖，不会出现某个区间一个点都没有的情况。同时样本量 n 可以由我们根据需要灵活设定，不会像正交表那样被因子和水平数卡死。工程上做代理模型实验设计时，LHS 基本属于默认选择，原因就是它在样本数量可控和空间覆盖率这两个指标之间取得了一个很实用的平衡。

1.3 为什么是二阶多项式而不是克里金或神经网络

代理模型的种类非常多，除了二阶多项式响应面，还有克里金插值、径向基函数、神经网络、支持向量回归等。我在实际项目里见过有人一上来就上神经网络，其实很多时候没有必要。选代理模型要看数据量、非线性程度、以及你是否需要解释性。

二阶多项式最明显的优势是形式透明，回归系数一眼就能看出哪个变量影响大、哪个变量几乎可以忽略。配合方差分析还能做显著性检验，这在工程报告和前期方案论证阶段是非常实用的。它的短板也很明显，如果目标函数在局部有剧烈的峰谷变化，或者高度非线性，二阶多项式的拟合精度会迅速下降。但在变量个数十个以内、响应较为光滑的常见工程问题上，它往往是性价比最高的起步方案。

提示：建议把二阶多项式响应面当成流程的第一版模型。用它跑通整个优化链路后，再根据精度需求决定是否换用克里金等更复杂的代理模型。一上来就用复杂模型，样本量不够、调参繁琐，反而容易翻车。

2. LHS 拉丁超立方采样原理与 MATLAB 实现

2.1 分层逻辑与均匀覆盖的了解

LHS 的核心思想可以理解成棋盘上放棋子。假设有两个变量，每个变量都想抽 5 个样本点，传统随机抽样就像往棋盘上随手撒豆子，有的格子可能挤了好几颗，有的格子一颗没有。LHS 的做法是先在 x 方向把棋盘切成 5 条等概率的带子，再在 y 方向也切成 5 条带子，要求最终的 5 个点必须分布在不同的行和不同的列里。当然，具体配对可以随机，但每一行每一列都必须恰好有一个点。这样就保证了无论从 x 方向看还是从 y 方向看，样本点都是均匀铺开的。

对于三个变量以上的高维空间，LHS 的推广方式是每个维度划分成等概率区间，然后在保证每个维度每区间恰好取一个样本的前提下，随机排列组合出 n 个样本点。这种设计会让样本点的边缘分布非常均匀，同时对变量之间的交互也有合理的覆盖。

这里有个细节值得注意：LHS 保证的是一维投影均匀，并不保证任意两个变量之间的联合分布也均匀。样本量不大时，某两个变量之间仍可能出现伪相关。所以抽样完成后做相关性检查是有必要的，后面我会专门说这个问题。

2.2 MATLAB 中的 LHS 命令与样本量选择

MATLAB 里做 LHS 最直接的工具是lhsdesign。基本用法是X = lhsdesign(n,k)，生成 n 行 k 列的矩阵，每一列对应一个变量，所有取值都在区间 [0,1] 内。如果你不想让点边集中在边界上，或者想控制迭代次数来提高均匀性，可以带上参数，比如：

n = 30; % 样本点数 k = 3; % 变量个数 X_norm = lhsdesign(n, k, 'Smooth', 'off', 'Criterion', 'maximin');

代码里'Criterion','maximin'表示让样本点之间的最小距离尽量大，这样点与点之间不会太拥挤，空间填充性更好。'Smooth','off'的作用是让点在每个区间内落在区间中心附近，而不是完全随机，这样能减少极端聚簇。这里要提醒一句，不同 MATLAB 版本的lhsdesign参数选项略有差异，如果你的版本提示参数不支持，直接在命令行里输入doc lhsdesign查一下当前版本的帮助文档为准。

抽样完成之后，因为lhsdesign默认输出单位区间里的值，需要把每一列映射回真实的物理取值范围：

lb = [0.5, 100, 0.001]; % 变量下界 ub = [2.0, 500, 0.01]; % 变量上界 X_phys = X_norm .* (ub - lb) + lb;

样本量 n 怎么定？这个问题新手最容易问。二阶多项式响应面的待定系数总数是 (k+1)(k+2)/2，其中 k 是变量个数。为了能拟合出全部系数，样本量 n 必须大于等于这个数。但实际做工程时，我一般不会卡着这个下限来取。原因很简单，样本点如果正好等于待定系数个数，回归模型表面上看能解出来，但没有任何冗余，一旦数据带点噪声，拟合结果就很脆弱。我的经验是先按 n = 4k 到 8k 来取，复杂问题或者预算充足时取 10k 也不过分。比如 3 个变量的二阶模型需要 10 个系数，取 30 到 40 个样本点已经比较稳。

2.3 样本空间均匀性检查

抽样完成后不要直接丢进仿真里，先花两分钟检查一下样本分布是否正常。一个很简单有效的做法是画散点图矩阵：

plotmatrix(X_phys);

看各个变量两两之间的散点分布，确认不存在明显的聚簇或空洞。再算一下相关性矩阵：

R = corr(X_phys); disp(R);

理想情况下，任意两个变量之间的相关系数应该接近 0。如果某个变量对之间的相关系数绝对值超过 0.5，说明这次抽样的随机配对出现了伪相关，最保险的做法是换随机种子重新抽一次，或者增大 n 重新抽样。

我在实际项目里遇到过一个小坑：有个同事用 LHS 做了一个 5 变量的采样，结果第 2 列和第 5 列之间的相关系数高达 0.72。他自己没注意，后面响应面模型里这两个变量的交互项系数特别奇怪，怎么看都解释不通，最后排查半天发现是样本相关性出了问题。所以这个检查环节真的不能省。

3. 二阶多项式响应面模型的构建与回归校验

3.1 模型形式、待定系数数量与样本下限

二阶响应面的标准形式是：

y = b0 + Σ(bi·xi) + Σ(bii·xi²) + ΣΣ(bij·xi·xj)

也就是常数项加一次项加平方项再加交互项。k 个变量时，待定系数的个数为：

n_coefs = (k+1)(k+2) / 2

如果样本数小于待定系数数，最小二乘问题就变成欠定问题，解不唯一，模型完全不可信。哪怕刚好等于待定系数数，拟合出来的响应面也会把所有样本点全部完美穿过，但样本点以外区域的预测往往一塌糊涂，这其实就是过拟合。

3.2 用 MATLAB 构造二阶设计矩阵并完成回归

MATLAB 里做多项式回归，我推荐直接用fitlm，因为它的输出里自带 R 方、F 检验和回归系数显著性 p 值，省去很多手动统计的麻烦。关键一步是构造二阶设计矩阵。假设原始样本存储在矩阵X_raw中，行是样本，列是变量，响应值存储在列向量Y中：

X_quad = [X_raw, X_raw(:,1).^2, X_raw(:,2).^2, X_raw(:,3).^2, ... X_raw(:,1).*X_raw(:,2), X_raw(:,1).*X_raw(:,3), ... X_raw(:,2).*X_raw(:,3)]; mdl = fitlm(X_quad, Y);

变量不多的时候，手工写这些平方项和交互项完全没问题。如果变量比较多，想自动生成多项式基函数，可以用x2fx：

X_quad = x2fx(X_raw, 'quadratic'); mdl = fitlm(X_quad(:,2:end), Y);

x2fx的第二列开始就是所有一次项、平方项和交互项。用fitlm的好处是，输出里每个回归系数都带一列pValue。一般来说，p 值小于 0.05 的项可以认为是显著项，p 值大于 0.2 的项可以考虑剔除后重新拟合，让模型更精简。如果显著项太少，或者整体模型 R 方很低，那就说明要么样本量不足，要么目标函数和输入参数的关系并不是简单的二阶关系。

3.3 拟合质量评估：R 方、调整 R 方与残差检查

拟合完之后，先看两个指标。第一个是决定系数 R²，它刻画的是模型解释了总变异有多大比例。第二个是调整 R²，它对变量个数做了惩罚，可以避免盲目增加项数带来的虚高 R²。工程上我给自己定的及格线是调整 R² 大于 0.9，如果这个值不到 0.9，说明模型很大程度上没有捕捉到真实响应趋势，这时候拿着它去做优化，最后找出来的“最优解”很可能是仿真评价下很差的点。

除了 R 方，残差分析也很关键。直接在 MATLAB 里画：

plotResiduals(mdl, 'fitted');

如果残差随机分布在零线上下，没有明显的喇叭口形状，说明模型假设基本合理。如果残差随着拟合值的变大而明显变大或变小，说明方差不是常数，经常是模型缺少高次项或者存在变量交互结构导致的，这时考虑换克里金模型或者补样本点，比单纯在二阶模型里加项更有效。

3.4 变量归一化与系数解释

这里想单独强调一下变量归一化，因为它直接关系到回归的数值稳定性和系数解释。物理量纲不同的变量，直接放进回归模型时，量级大的变量可能在数值计算中直接影响矩阵条件数，导致求解精度下降甚至失败。标准做法是在建模前将每个变量映射到 [-1, 1] 或 [0, 1] 区间：

X_mean = mean(X_phys, 1); X_std = std(X_phys, 1); X_norm = (X_phys - X_mean) ./ X_std;

用归一化后的变量建模，得到的回归系数之间就可以直接比较大小。系数绝对值越大，说明该变量或者该交互效应对响应的影响越强。我在做项目报告时经常把这个归一化系数表直接放在结论页，效果很好，评审一眼就能看出来哪些参数是主导因素。

注意：归一化用的均值和标准差要保存下来。后续任何新点想用这个模型做预测，必须先做相同的归一化变换，否则预测结果没有意义。建议通过mdl.Coefficients.Estimate导出回归系数后，把归一化参数和系数一起存成字典或者结构体，写成一个独立的预测函数。

4. 单目标快速寻优：非线性规划求解

4.1 用响应面替换仿真作为目标函数

响应面建好之后，优化目标就从黑盒仿真变成了一个显式的二次多项式，这时候很多经典优化算法都可以用了。最常用的是 MATLAB 的fmincon，它用于求解带约束的非线性规划问题。实践里我通常是直接从mdl对象中取出回归系数，然后构造一个匿名函数：

b = mdl.Coefficients.Estimate; fun = @(x) b(1) + b(2)*x(1) + b(3)*x(2) + b(4)*x(3) + ... b(5)*x(1)^2 + b(6)*x(2)^2 + b(7)*x(3)^2 + ... b(8)*x(1)*x(2) + b(9)*x(1)*x(3) + b(10)*x(2)*x(3);

注意这里的x必须使用和建模时完全相同的变量顺序和归一化尺度，所以我一般会把归一化参数合并在一个函数里写，避免每次手动转换时混乱：

fun = @(x) predictResponse(x, mdl, X_mean, X_std);

predictResponse是自己写的一个小函数：先对输入做归一化，再调用mdl.predict就行了。这种做法虽然增加了一次函数调用的开销，但直观、不易出错。

fmincon需要一个初始点x0，并且支持同时添加线性约束、边界约束和非线性约束：

lb = [0.5, 100, 0.001]; ub = [2.0, 500, 0.01]; A = []; b = []; Aeq = []; beq = []; options = optimoptions('fmincon', 'Display', 'iter', 'Algorithm', 'sqp'); [x_opt, f_opt] = fmincon(fun, x0, A, b, Aeq, beq, lb, ub, @nonlcon, options);

非线性约束函数nonlcon要写成[c, ceq] = nonlcon(x)的格式，其中c是不等式约束，需要满足 c(x) ≤ 0；ceq是等式约束，需要满足 ceq(x) = 0。没有等式约束时直接返回空矩阵。

4.2 多起点策略：绕开局部最优

用fmincon求解响应面模型虽然在单次计算上很快，但存在一个绕不过去的问题：对于一个二次多项式来说，目标函数如果是非凸的，那么局部最优解并不等于全局最优解。很多人一开始直接拿一个随机初值去跑，跑出来一个解就当结果了，这是有风险的。

我自己的习惯是从 LHS 样本点中挑预测值最优的几个点作为多个初始值，再在它们附近加一点随机扰动，分别跑fmincon，最后从多个结果里选目标值最优点作为最终解：

num_trials = 8; x_best = []; f_best = inf; for i = 1:num_trials x0_try = X_phys(randi(n), :) + 0.05 * randn(1, k); x0_try = max(lb, min(ub, x0_try)); [x_tmp, f_tmp] = fmincon(fun, x0_try, A, b, Aeq, beq, lb, ub, @nonlcon, options); if f_tmp < f_best f_best = f_tmp; x_best = x_tmp; end end

初值点的选择有一定讲究，直接覆盖设计空间几个不同区域，比全部挤在一个角落要好得多。

4.3 单目标方法的适用边界

fmincon虽然快，但本质上只能解决单目标问题。如果我们要同时优化两个互相冲突的目标，比如既想让结构重量最小，又想让刚度最大，这时简单加权合并成单目标的策略容易出现主观性强和权重敏感性高的问题。碰到这种需求，下一步就要换多目标遗传算法来跑帕累托前沿。

5. 多目标遗传算法求解帕累托前沿

5.1 gamultiobj 的用法与配置

MATLAB 的多目标优化工具是gamultiobj，基于 NSGA-II 改进算法，可以直接输出一组帕累托最优解集。使用它的第一步是定义多目标函数：

fitness = @(x) [response_one(x), response_two(x)];

第一个目标可以是重量，第二个可以是最大应力。然后设置求解器参数：

options_ga = optimoptions('gamultiobj', ... 'PopulationSize', 100, ... 'MaxGenerations', 200, ... 'ParetoFraction', 0.35, ... 'Display', 'iter'); [x_pareto, f_pareto] = gamultiobj(fitness, k, A, b, Aeq, beq, lb, ub, @nonlcon, options_ga);

ParetoFraction控制在最终非支配解集里的比例，默认 0.35 通常够用。种群大小和迭代代数需要根据问题复杂度调整，变量 3 到 5 个时种群 100、代数 200 是一个比较保守的起步设置。

这里注意一个容易踩坑的地方：由于目标是响应面模型，计算本身非常快，所以遗传代数可以放心加大。如果你在真实仿真上跑遗传算法，200 代根本不敢想，但在响应面上跑 500 代也就几秒钟的事。

5.2 帕累托前沿的解读与收敛判断

跑完之后，最常见的操作是画出两个目标的散点图：

scatter(f_pareto(:,1), f_pareto(:,2), 20, 'filled'); xlabel('Objective 1'); ylabel('Objective 2');

画出来的这条前沿上，每个点都是一个非支配解，意思是任何一个目标要继续变好，另一个目标必然会变差。判断算法是否收敛，可以分两批跑一次，比如第一批跑 150 代，第二批跑 300 代，看前沿是否变化不大。如果两批结果差得很远，说明代数不够或者种群太小。

我自己的做法是直接设置一个相对较高代数和较大种群，因为反正是在响应面上跑，代价小。完全没必要省这个计算量。

5.3 从帕累托前沿到最终决策

多目标优化跑完之后，很多人不知道下一步怎么选点。帕累托前沿提供的是候选集，不是唯一答案。实际做决策时，通常要看几个位置：端点代表某个目标最优但另一个目标很差的极端方案，拐点代表两个目标折中比较好的方案。如果没有明确偏好，可以用一些客观权重方法比如熵权法或者 TOPSIS 从前沿中挑一个综合排名最高的点，然后把这个点丢回真实仿真里做验证。

我在地下项目里试过多次，前沿拐点附近的设计点往往是最平衡的，拐点判断可以直接从图上目视，也可以计算相邻点构成的折线变化来辅助判断。最终选定设计后，必须回到仿真中完整计算一次真实目标值，并和响应面预测值做误差比较，这一步余量和闭环验证是绝对不建议跳过的。

6. 常见问题与排查技巧实录

6.1 问题速查表

下面这张表是我在多个项目里反复遇到并整理出来的高频问题与排查路线，直接对着现象查就行。

6.2 关于代理模型的闭环验证

前面反复提到仿真验证，但很多初次接触这套流程的读者可能低估了它的重要性。响应面模型本质上是对真实仿真的一种近似，近似必然有误差。优化算法找到的那个最优点，在响应面上可能是一座漂亮的山峰，但在真实仿真里可能只是一个普通坡地。

所以一套严谨的操作流程应当是这样：在 LHS 采样点基础上建响应面，跑优化，得到候选解，回到原仿真里计算这个候选解的真实目标值。如果预测值和真实值误差小于 5%，说明模型在这个区域可信；如果误差大于 5%，就把这个点作为新增样本补进原来的数据集合里，重新拟合响应面，再跑优化。这种“加点迭代”的策略非常实用，每一次补点都会让模型在最优区域附近变得更准。

我第一次做这类项目时没有做闭环验证，直接拿着响应面优化结果交差，后来实测数据出来差了将近 12%，虽然方案方向是对的，但具体参数几乎不能用。从那以后，我就养成了“优化解必验证、验证完必补点”的习惯。

6.3 不同阶段适配不同复杂度的模型

最后想再强调一点我自己的取舍逻辑：如果工程精度要求不高，或者只是做初步设计阶段的方案筛选，二阶多项式响应面搭配 LHS 和遗传算法是性价比极高的选择。如果变量个数大于 12，或者响应面精度一直提不上去，可考虑克里金或径向基函数模型。但要注意，换复杂代理模型时样本量必须相应提高，否则复杂模型会因为数据不足而发生过拟合，反而比简单模型更差。

从我的经验来看，LHS、二阶多项式响应面和遗传算法的搭配组合，特别适合把单个仿真时长较长但设计变量不算太多的问题快速跑通。它也许不是理论上最精确的方法，却是在项目周期和计算资源限制下最稳健、最能落地的解法。这也是我一直在各种场景下推荐先把它作为第一版优化流程的原因。