先说个场景你手里有个仿真模型算一次结果可能要几十秒甚至几分钟而领导的要求是“把几个关键参数调到最优”。如果你直接用遗传算法去套这个仿真器跑300代、每代100个个体意味着要调用仿真器三万次按每次1分钟算得跑20天项目早就黄了。就算你狠下心等得起仿真器偶尔数值不收敛、中途报错整个优化进程也容易直接崩掉。我最早碰这套东西也是被逼的。当时要在一个设计空间里找一组参数让某个结构件的综合性能最优仿真模型又贵又慢直接套遗传算法根本不现实。后来用的就是标题里这条路先用拉丁超立方采样在设计空间里均匀抽取几十个点把这几十个点跑完仿真拿结果去拟合一个二阶多项式响应面模型把这个模型当“替身”接下来不管是用非线性规划还是遗传算法都是在替身上做文章速度快到可以忽略不计。等优化算法给出一个推荐解最后再调一次原始仿真或者实验做复核。整套流程在MATLAB里实现非常顺手而且对新手也足够友好。下面的内容我不打算讲太多抽象理论重点是把代码怎么写、每一步为什么这么干、哪些坑我踩过原原本本讲清楚。如果你也在做代理模型辅助优化、试验设计、参数寻优这一类工作这篇应该能帮你省不少时间。1. 这套技术组合到底是干什么的1.1 问题源头一次仿真太贵很多工程优化问题的难点不在优化算法本身而在于目标函数的计算成本。结构有限元分析、流体仿真、电磁场仿真单次求解耗时可长可短但普遍不便宜做实验就更不用提了一组试件做下来可能就得好几天。你手上有一堆设计变量想在约束范围内找最优理论上可以用优化算法直接驱动仿真器每一次“算一版参数好不好”就要完整跑一遍仿真这个代价通常承受不起。这时候工程上最主流的思路就是“用响应面当替身”。说白了既然原模型太贵我先在设计空间里选一批有代表性的试验点跑有限次数仿真拿到“输入-输出”的对应关系然后用一个数学表达式拟合这批数据让表达式尽量逼近真实关系后面反复调用这个表达式做优化只在最后阶段用真实仿真器对推荐的候选参数做复核。这样真实仿真的调用次数被压缩到几十次整个优化过程在几分钟甚至几秒内完成。1.2 为什么是LHS、二阶响应面和混合优化这么凑在一起采样、建模、寻优这三个阶段各自有大量候选方法但它们不是随便拼的。拉丁超立方采样LHS抽出的点数量可控且覆盖均匀正适合给二阶多项式回归提供训练数据二阶多项式模型结构简单、有明确的系数表达式既可以解释参数影响方向又能非常廉价地求梯度非线性规划和遗传算法则分别承担“局部精搜”和“全局粗搜”的角色两者组合能把代理模型的价值榨干。这套组合最吸引人的地方在于它很“皮实”。哪怕你不打算深究数学细节只要照流程把采样点做了、把fitlm跑起来、再把优化器调用起来大概率能得到一个可用结果。对大多数工程场景来说我们不追求数学意义上的绝对最优而是追求“用尽量少的昂贵仿真换一个足够好的可行方案”这正好是这套方法的强项。1.3 适合谁来参考如果你做的是参数优化相关的工作手里有一个可以输入设计变量、输出性能指标的黑盒仿真或实验系统并且这个系统单次计算成本不低那么这篇文章的内容可以直接落到你的项目里。如果你刚接触代理模型和优化算法连LHS样本怎么生成、响应面代码怎么下笔都不太清楚这篇文章同样能作为一条完整可跑的入门路径代码拿过去改改边界和函数就能用。2. 用LHS合理取点采样不是越多越好2.1 拉丁超立方采样到底在做什么设计空间里取点这个事最朴素的做法是随机抽样和全因子网格。但全因子网格在变量多的时候组合数爆炸比如5个变量、每个变量取10个水平就是10万次组合实际工程中根本跑不完。随机抽样又容易出现扎堆现象点与点挤在一起样本信息冗余有些设计区域却完全没有覆盖到。LHS的核心思想是分层。把每个设计变量的取值范围分成n个等概率小区间然后保证每个变量在每个小区间里只被取到一个点。这样说比较抽象我用二维的例子解释假设两个变量都在0到1之间取10个样本点如果画一个10乘10的网格LHS会确保在这10行里每行恰好有一个点在这10列里每列也恰好有一个点。这些点的分布像棋盘上互不攻击的车整体铺得比较均匀同时样本量可以由你自由控制。与简单随机抽样相比LHS能以小得多的样本量覆盖整个设计空间与全因子相比它又把实验次数压缩到可以接受的范围。所以LHS特别适合给响应面建模当训练数据来源。2.2 MATLAB里一分钟生成漂亮样本点MATLAB里生成LHS样本非常简单核心函数就是lhsdesign。最基础的用法是% 生成 n 行 d 列的样本取值在 [0,1] 区间 n 40; % 样本点数 d 2; % 设计变量个数 X01 lhsdesign(n, d);默认情况下lhsdesign会返回n行d列的矩阵每一列是一个设计变量在[0,1]区间的取值。如果直接拿这个矩阵去建模物理含义不对因为实际变量的上下界往往不是0和1。因此要把样本映射到真实设计空间。我实际用的时候通常会再加一个参数criterion设为maximin也就是让样本点之间的最小距离尽可能大避免出现低质量的拥挤样本。lhsdesign默认其实已经做了一定优化但显式指定这个准则会更稳X01 lhsdesign(n, d, criterion, maximin, iterations, 100);iterations控制优化过程的迭代次数越大越好但耗时也相应增加。对常规工程问题100次迭代已经很足够了。2.3 从[0,1]区间映射到设计空间有了标准化样本X01之后通过简单线性映射就能得到真实变量lb [-5, -5]; % 各变量下界 ub [5, 5]; % 各变量上界 X lb X01 .* (ub - lb);这么做本质上是把[0,1]区间按比例平移到实际的上下界范围里。如果各个变量的数量级差异很大这个映射方式本身也自带了归一化的逆过程训练回归模型时高阶项和交叉项的数值稳定性会更好一些。这也是LHS样本相比直接在实际坐标系里随机抽样的又一个优势——至少能让你下意识地先想清楚设计空间的边界。2.4 样本数量怎么定才不浪费样本数量少了二阶多项式回归系数都估计不准样本数量多了意味着大量昂贵的仿真调用违背了代理模型的初衷。我个人的经验公式是和待估系数个数挂钩。假设有d个设计变量二阶多项式回归模型包含常数项、一次项、二次项和所有两两交叉项总系数个数是m (d1)(d2)/2当d2时m6也就是 y b0 b1x1 b2x2 b3x1x2 b4x1^2 b5x2^2。当d6时m28意味着你至少要准备28个点回归方程才可解但为了拟合稳定和有一定余量采样数量最好在系数个数的1.5到2倍以上。实际项目中我一般这样把握变量2到4个时采样数量取30到60个变量5到8个时采样数量取60到120个变量超过10个后二阶多项式模型的系数数量会涨得很快这时候要考虑先用筛选试验剔除次要变量或者换更强的代理模型Kriging、径向基函数等。2.5 固定随机种子保证结果可复现LHS本身带随机性每次运行生成的样本可能不同。这个特性本身不是问题但如果你在调试代码、对比不同建模方案的效果结果一会一个样会非常痛苦。所以我每次写采样代码第一件事就是固定随机种子rng(2024);这样无论你运行多少次只要样本量和变量边界不变得到的LHS样本就完全一致。工程交流中这一点也很重要你的同事拿到你的脚本后跑出来的结果和你文章里写的对得上才能顺利复现。3. 二阶多项式回归响应面把黑盒模型“翻译”成公式3.1 二阶多项式模型长什么样响应面建模就是用一个相对简单的数学函数去逼近仿真系统输入和输出之间的关系。二阶多项式回归是其中最常用的一种原因是它能表达三类信息每个变量单独对响应的影响线性项和平方项变量之间的交互作用交叉项响应在空间中的弯曲变化趋势。以两个变量为例模型形式是这样的y b0 b1x1 b2x2 b3x1x2 b4x1^2 b5x2^2b0是常数项b1、b2反映主效应b3反映x1和x2的交互作用b4和b5反映各自的曲率。用生活化的方式理解一阶模型相当于在一块斜板上拟合平面二阶模型则在斜板的基础上加上了“碗形”“马鞍形”等弯曲形态表达能力明显更强。工程中很多性能指标在局部范围内趋势连续平滑二阶模型已经可以拟合得足够好。3.2 用fitlm一步完成回归拟合与显著性检验MATLAB里做多项式回归最顺手的函数是fitlm。它对常规最小二乘回归做了完整封装返回一个LinearModel对象里面包含系数估计、p值、R方、均方根误差等一堆现成结果。对单个响应变量拟合二阶模型的代码% X 是采样点矩阵每行一个样本每列一个变量 % y 是对应每个样本的仿真响应值列向量 mdl fitlm(X, y, quadratic);这里的关键是第三个参数quadratic它告诉fitlm要生成所有一次项、二次项和两两交叉项对应前面说的完整二阶多项式结构。相比自己手写一大串项名直接用关键字最省事。拟合完成后我最先看的永远是这几个结果mdl.Rsquared.Ordinary % 普通R方 mdl.Rsquared.Adjusted % 调整R方 mdl.RMSE % 均方根误差 mdl.Coefficients % 系数表包含估计值、标准误、t统计量、p值如果一个项的p值大于0.05说明该项对响应的解释能力不显著可以考虑手动去掉。不过在代理模型用于优化时我通常不会太激进地删项因为优化器最终要在整个设计空间里横冲直撞少一个交互项可能影响的是边界区域的预测精度。3.3 模型质量好坏不能只看R方很多新手拿到模型先看R方R方接近0.99就觉得万事大吉。这是个非常危险的误区。R方衡量的是“训练样本点上模型能解释多少响应波动”可响应面建模真正关心的是“没参与训练的位置预测得准不准”。我常用的验证方式有三种第一看调整R方Adjusted R-squared与普通R方的差距。如果一个模型拼命加项普通R方必然上升但调整R方会在加入无效项后下降两者的差值如果过大基本说明模型存在过拟合风险。第二看残差图。fitlm自带绘图函数plotResiduals(mdl, fitted);如果残差随机散布在零线附近没有喇叭状扩张或明显弯曲说明模型假设基本成立如果残差出现系统性规律说明二阶模型可能不够用了或者有重要变量没进来。第三留出验证集。我在完整案例里会演示这种做法把LHS采样点分成训练集和测试集模型只用训练集拟合然后对测试集做预测算预测值和真实仿真结果的误差。这个误差才是代理模型真实预测能力的体现。3.4 模型精度不够时的方向性思路如果测试集误差偏大不要急着调优化参数问题多半出在代理模型这一环。我的排查顺序大致是这样第一步检查采样数量是不是太少太少就补点重新建模。第二步检查设计变量的范围和响应量级是否差异过大如果差异大对变量和响应做标准化处理后再拟合。第三步考虑设计空间里是否存在剧烈突变。二阶多项式天然擅长拟合光滑响应如果真实模型存在阶跃、谐振或强非线性多项式响应面会非常吃力。这时候要么把设计空间切分成子区域分别建模要么换Kriging或径向基函数这类更灵活的代理模型。我的建议是不要在多项式模型上死磕该换方法就换方法。多项式模型的优势是简单、好解释、好求梯度但它的表达上限也在那里。4. 非线性规划与遗传算法组合起来做优化4.1 两类优化算法的定位完全不同响应面模型建好之后优化搜索可以在代理模型上飞速执行。此时面临的问题是用哪个优化器来搜非线性规划在MATLAB里对应fmincon和遗传算法ga是两类思路完全不同的方法。fmincon属于基于梯度的局部优化算法速度极快对光滑问题收敛精度高但它极度依赖初始点——初始点选得不好很容易收敛到局部最优解。遗传算法属于元启发式全局优化算法靠种群进化和随机搜索在大范围内寻找最优不要求目标函数光滑可导但收敛慢最终精度通常不如梯度法。既然代理模型的计算成本极低最优策略就显而易见了让遗传算法先做全局粗搜找到有希望的区域然后把这个结果作为fmincon的初始点让fmincon在局部精修。4.2 用fmincon做非线性规划初始点多试几个如果在代理模型上直接用fmincon目标函数就是predict。假设我们用单个响应目标objF代码如下objF (x) predict(mdl, x(:)); x0 [0; 0]; % 初始点 opts optimoptions(fmincon, Algorithm, interior-point, Display, final); [x_opt, f_opt] fmincon(objF, x0, [], [], [], [], lb, ub, [], opts);这里有一个细节fmincon传递给目标函数的决策变量默认是列向量而predict函数期望输入是行向量或者矩阵格式所以匿名函数内部用x(:)把向量改成行向量。这个小坑很多人第一次跑都会碰到症状是预测结果维度对不上。因为fmincon对初值敏感我的习惯是不要只从一个初始点出发。比较实用的做法是生成5到10个随机初始点分别跑fmincon然后综合比较结果。当然这里说的随机也可以用LHS来生成初始点效果会更均匀。4.3 用ga做全局搜索拿到初值再精修MATLAB的遗传算法调用格式同样干净opts_ga optimoptions(ga, PopulationSize, 100, MaxGenerations, 300, Display, final); [x_ga, f_ga] ga(objF, d, [], [], [], [], lb, ub, [], opts_ga);ga的第二个参数d是设计变量个数这里就是2。种群大小和最大代数是两个最需要关注的参数。种群太小搜索容易过早收敛种群太大虽然搜索充分但对于代理模型来说反正计算便宜大一点无所谓。既然是在代理模型上跑通常我会把种群和代数都设置得偏大一些让搜索更充分。跑完ga之后再用x_ga作为fmincon的初值精修一圈[x_final, f_final] fmincon(objF, x_ga, [], [], [], [], lb, ub, [], opts);这是整个代理优化流程里最值得养成的习惯全局粗搜负责定位有希望的山头局部精搜负责爬上山头最高点。两个步骤分工明确缺一不可。4.4 多目标优化gamultiobj和Pareto前沿工程优化里经常不只有一个目标。比如既要性能好又要重量轻既要精度高又要成本低。这种多目标问题通常不存在一个让所有目标同时达到最优的解我们需要的是Pareto前沿——在这条前沿上想改善任何一个目标都会牺牲另一个目标。MATLAB的多目标遗传算法函数是gamultiobj用法和ga类似objM (x) [predict(mdl1, x(:)), predict(mdl2, x(:))]; opts_multi optimoptions(gamultiobj, PopulationSize, 100, MaxGenerations, 300, Display, final); [x_pareto, f_pareto] gamultiobj(objM, d, [], [], [], [], lb, ub, [], opts_multi);mdl1和mdl2分别是对两个响应建立的二阶多项式模型。x_pareto是Pareto前沿上的解集矩阵f_pareto是这些解对应的两个目标预测值。注意这里目标函数返回的是一个包含两个预测值的行向量两个目标分别来自两个独立的响应面模型。实际项目中如果两个响应的数量级差异很大建议在展示和决策时对目标值做归一化处理否则Pareto前沿的几何形状会被量纲大的目标带偏。4.5 Pareto前沿算出来之后怎么选点很多刚接触多目标优化的同学跑出Pareto前沿就不知道该干什么了其实真正需要人拍板的部分才刚刚开始。Pareto前沿给的是“一票候选方案”最终工程方案要从里面挑一个。常用的选点方法有几种。工程偏好明确时直接根据权重确定各目标优先级没有明确偏好时可以用“到理想点最近”的策略在Pareto前沿上找到离“理想点”最近的一个点。理想点就是把每个目标单独优化时的最优值拼成的合成向量实际中虽然达不到但它提供了一个参照基准。这个选点策略的MATLAB代码大概是这样% 在Pareto前沿中找一个折中解到理想点的归一化欧氏距离最近 ideal min(f_pareto, [], 1); f_norm (f_pareto - ideal) ./ (max(f_pareto) - min(f_pareto) 1e-12); dist sqrt(sum(f_norm.^2, 2)); [~, idx_best] min(dist); x_compromise x_pareto(idx_best, :);需要特别提醒的是这里的f_pareto是代理模型的预测。用选出来的x_compromise再去调用一次真实仿真器看看真实目标值是否和代理预测一致这是整个流程中绝对不能省略的步骤。5. 一个可直接复现的完整案例从采样到寻优写通5.1 算例定义一个返回两个响应的“仿真器”为了演示完整流程我构造一个“仿真器”。现实中它可能是有限元模型、流体仿真模型或者实验系统这里我用一个带噪声和周期项的双响应数学函数来模拟它返回两个响应y1和y2分别代表两个相互冲突的工程性能目标function y mySimulator(X) x1 X(:,1); x2 X(:,2); % 响应1可以理解为“性能指标”越小越好 y1 (x1 - 2).^2 2*(x2 1).^2 0.5*sin(2*x1).*cos(0.8*x2); % 响应2可以理解为“成本或重量”越小越好 y2 (x1 1).^2 1.5*(x2 - 1.5).^2 0.8*cos(1.2*x1).*sin(0.6*x2); y [y1, y2]; end请你记住这个函数在实际项目中要被替换成你自己的仿真模型。它的特点是两个目标的最优解不重合——一个想让x1靠近2、x2靠近-1另一个想让x1靠近-1、x2靠近1.5所以天然存在折中问题非常适合演示多目标优化。5.2 第1步LHS采样并生成训练数据我先写主脚本生成LHS样本调用仿真器得到训练输出%% 1. 清理环境并固定随机种子 clear; clc; close all; rng(2024); %% 2. 问题定义 lb [-5, -5]; ub [5, 5]; d length(lb); n 40; % LHS样本总量 %% 3. 拉丁超立方采样 X01 lhsdesign(n, d, criterion, maximin, iterations, 100); X lb X01 .* (ub - lb); %% 4. 调用“真实仿真器”计算响应 Y mySimulator(X); % 分别取出两个响应 y1 Y(:,1); y2 Y(:,2);可以看到采样这段代码非常短真正贵的是第4步的mySimulator调用但整个流程只需要调用它40次。如果换成真实仿真可能需要在程序里等待较长时间但这是无法避免的核心成本。5.3 第2步训练/测试集划分并进行二阶多项式回归为了客观验证模型预测能力我不会把所有样本都用于训练而是随机抽一部分作为测试集最后用测试集上的预测误差来评估建模质量%% 5. 划分训练集和测试集 idx_perm randperm(n); n_train 32; idx_train idx_perm(1:n_train); idx_test idx_perm(n_train1:end); %% 6. 对两个响应分别建立二阶多项式响应面 mdl1 fitlm(X(idx_train,:), y1(idx_train), quadratic); mdl2 fitlm(X(idx_train,:), y2(idx_train), quadratic); %% 7. 在测试集上评估模型 y1_pred predict(mdl1, X(idx_test,:)); y2_pred predict(mdl2, X(idx_test,:)); err1 y1_pred - y1(idx_test); err2 y2_pred - y2(idx_test); fprintf(响应1测试集RMSE: %.4f\n, sqrt(mean(err1.^2))); fprintf(响应2测试集RMSE: %.4f\n, sqrt(mean(err2.^2))); % 也可以画一个预测值-真实值的散点图来直观判断 figure; plot(y1(idx_test), y1_pred, o); hold on; plot(y2(idx_test), y2_pred, s); plot([min(Y(:)), max(Y(:))], [min(Y(:)), max(Y(:))], k--); xlabel(真实响应值); ylabel(模型预测值); legend(响应1, 响应2, yx参考线, Location, best);训练/测试集划分是很多人写代理模型时会忽略的一步但它在工程上的价值非常大。测试集上的RMSE代表了模型“外推”到未采样区域的实际表现如果这个误差大后面所有优化结果都要打问号。5.4 第3步用非线性规划和遗传算法做单目标优化工程上如果只有一个综合目标可以把两个响应加权合并成一个。这里我用权重w10.65和w20.35代表项目上更看重y1一点%% 8. 单目标加权优化用代理模型代替真实仿真 w1 0.65; w2 0.35; objS (x) w1 * predict(mdl1, x(:)) w2 * predict(mdl2, x(:)); % 先跑遗传算法做全局粗搜 opts_ga optimoptions(ga, PopulationSize, 100, MaxGenerations, 300, Display, final); [x_ga, f_ga] ga(objS, d, [], [], [], [], lb, ub, [], opts_ga); % 再用遗传算法结果作为初值fmincon做局部精修 opts_fmin optimoptions(fmincon, Algorithm, interior-point, Display, final, MaxFunctionEvaluations, 5000); [x_final, f_final] fmincon(objS, x_ga, [], [], [], [], lb, ub, [], opts_fmin); %% 9. 把优化结果放回真实仿真器验证 y_true mySimulator(x_final(:)); y_true_weighted w1*y_true(1) w2*y_true(2); fprintf(ga得到最优解: (%.4f, %.4f), 代理模型目标: %.4f\n, x_ga(1), x_ga(2), f_ga); fprintf(fmincon精修后最优解: (%.4f, %.4f), 代理模型目标: %.4f\n, x_final(1), x_final(2), f_final); fprintf(真实仿真验证加权目标: %.4f\n, y_true_weighted);如果你跑完会发现代理模型的预测值和真实仿真器算出的值有差异但只要差异在可接受范围内这个代理优化流程就成立。如果差异超出预期就要回头补样本、优化模型然后再跑优化。5.5 第4步多目标优化与Pareto前沿决策接下来用gamultiobj在同一个代理模型上做多目标优化%% 10. 多目标优化直接用两个响应面模型作为目标 objM (x) [predict(mdl1, x(:)), predict(mdl2, x(:))]; opts_multi optimoptions(gamultiobj, PopulationSize, 120, MaxGenerations, 400, Display, final); [x_pareto, f_pareto] gamultiobj(objM, d, [], [], [], [], lb, ub, [], opts_multi); %% 11. 可视化Pareto前沿 figure; plot(f_pareto(:,1), f_pareto(:,2), o); xlabel(响应1预测值性能指标); ylabel(响应2预测值成本指标); title(Pareto前沿代理模型预测); grid on; %% 12. 折中解到理想点最近的解 ideal min(f_pareto, [],