简介面向多目标优化研究者与MATLAB开发者该资源提供基于超体积期望改进HVEI算法的MATLAB实现适用于在成本、效率、安全性等互相冲突的目标中获取Pareto最优解集的场景特别适合有一定优化基础、希望获得可运行代码的读者。压缩包共7个文件6个.m代码文件涵盖核心HVEI计算、二维超体积计算以及高斯分布函数等辅助模块并提供示例脚本演示调用流程另含1个txt说明文件讲解使用要点整体大小仅4KB轻量精简便于快速部署和调试。资源当前已有379人学习下载。读者借助代码可完整复现HVEI模型理解期望改进如何引导搜索最大化超体积覆盖同时可根据自身问题重新定义目标函数与约束调整种群初始化与迭代逻辑从而迁移到工程优化、能源调度等实际任务中。1. 从文件清单看 HVEI 优化模型的结构设计拿到HV_based_expected_improvement.zip解压以后就是六个.m文件和一份README.txt。Untitled.m是入口脚本exipsi.m、exi2d.m负责期望改进的计算hvolume2d.m算二维超体积gausspdf.m和gausscdf.m提供高斯分布的概率密度和累积分布函数。这个包解决的不是普通单目标优化而是多目标优化里最难的一类问题如何在不知道全局 Pareto 前沿的前提下用尽可能少的真实评估找到一组覆盖度好、分布均匀的非支配解。HVEI 的核心思想是把单目标贝叶斯优化里的 Expected Improvement 推广到超体积指标上每一步都选择对当前非支配集超体积增量期望最大的候选点在探索未知区域和利用已有优势之间取得平衡。如果你正在做多目标贝叶斯优化、代理模型辅助设计或者想在一个比较干净的 MATLAB 代码基础上改造出适合自己的采集函数这个包值得读一遍。它没有依赖额外工具箱核心逻辑都写在几个函数里适合边读边改。2. 超体积与期望改进HVEI 的数学原理和 MATLAB 实现基础2.1 超体积HV的定义与 hvolume2d.m 计算超体积是指由非支配解集与参考点围成的多维区域大小在多目标优化中它是唯一一个同时满足严格单调性、不受目标尺度影响只要参考点一致的纯量化指标。二维问题里超体积就是若干矩形面积并集。hvolume2d.m正是针对二维情况做了快速实现输入是前沿点集和参考点输出面积。常见的实现逻辑是function hv hvolume2d(points, ref) % POINTS: 每个点为一行目标值越小越好 % REF: 参考点例如 [1.1, 1.1] pts sortrows(points, 1); % 按第一维排序 hv 0; prev_y ref(2); for i 1:size(pts, 1) if pts(i,2) prev_y hv hv (pts(i,1) - ref(1)) * (prev_y - pts(i,2)); prev_y pts(i,2); end end hv abs(hv); end这里按第一维升序排列prev_y从参考点第二维开始逐步下移累计每个点的矩形条面积。注意如果目标是最小化那么参考点应当大于所有目标值而排序后第一个点贡献面积从参考点往左abs用于防止负面积。实际使用时要先做非支配排序把被支配的点过滤掉否则面积会偏大。对于三维以上二维扫描线就失效需要改用deldir或蒙特卡洛采样这也是这个包里只提供二维hvolume的原因。2.2 高斯密度与分布函数gausspdf.m 和 gausscdf.mHVEI 的期望运算离不开高斯分布。代理模型给出的预测不只是一个均值还有一个方差可以用正态分布刻画目标值的不确定性。gausspdf.m和gausscdf.m提供了标准高斯 pdf 和 cdf防止依赖工具箱版本。其内部实现等价于function y gausspdf(x, mu, sigma) if nargin 1, mu 0; sigma 1; end y exp(-0.5 * ((x - mu) ./ sigma).^2) ./ (sigma * sqrt(2 * pi)); end function p gausscdf(x, mu, sigma) if nargin 1, mu 0; sigma 1; end p 0.5 * (1 erf((x - mu) ./ (sigma * sqrt(2)))); end两个函数都支持向量输入参数mu和sigma与x同维度或标量。和 MATLAB 内置的normpdf、normcdf相比自定义版本少一次函数寻址开销在循环里调用更快。但数值上要注意sigma为 0 时会出现除零后续章节会给出处理建议。2.3 从单目标 EI 到多目标 HVEI单目标期望改进定义是EI(x) E[ max(0, f_best - f(x)) ]当预测分布为高斯时有解析表达式EI(x) (f_best - μ) * Φ((f_best - μ)/σ) σ * φ((f_best - μ)/σ)其中Φ是 cdfφ是 pdf。多目标场景下f_best变成当前 Pareto 前沿改进量变成“新点加入后超体积的提升量”。HVEI 的自然定义HVEI(x) E[ HV(P ∪ {f(x)}) - HV(P) ]但f(x)是随机向量各目标间可能有相关性通常先假设独立再利用联合分布做数值积分。这个包里的exipsi.m就是把连续积分转化为有限求和配合二维超体积的解析式得到较精确的 HVEI 近似值。函数文件作用依赖hvolume2d.m计算二维超体积排序gausspdf.m标准正态密度无gausscdf.m标准正态分布函数无exi2d.m二维HVEI积分gausspdf, gausscdfexipsi.mHVEI及改进概率近似exi2dUntitled.m主优化循环以上全部3. 核心函数拆解exipsi.m 和 exi2d.m 的 HVEI 计算流程3.1 exipsi.m 的签名与结构exipsi.m的名称来自 Expected Improvementexi和 Probability of Improvementpsi。在实际实现中它接收预测均值、标准差、当前 Pareto 集和参考点返回期望超体积增量和改进概率。框架可以这样写function [exi, psi] exipsi(mu, sigma, pareto, ref) if length(mu) ~ 2 % 高维问题退化为采样近似本包未实现时提示 error(当前版本只支持二维目标); end [exi, psi] exi2d(mu, sigma, pareto, ref); end从这个结构可以看出作者刻意保留了扩展接口二维走解析路线高维可以在这里替换成蒙特卡洛采样或数值积分。psi的实用价值在于给采集函数提供另一个参考当 HVEI 接近 0 但改进概率很高时说明该点性能波动大值得赌一把反之 HVEI 高而psi低说明期望收益来自极端情况实际风险大。3.2 exi2d.m 的二维积分实现二维 HVEI 计算的核心是对于任意候选点预测分布 (μ1,σ1) 和 (μ2,σ2)求出其目标向量落在当前前沿“超体积缺口”内的概率与贡献量的乘积积分。常见的做法是先由当前 Pareto 前沿生成一组不重叠矩形区域再对每个矩形区域做二维正态期望计算。伪代码级别的实现如下function [exi, psi] exi2d(mu, sigma, pareto, ref) exi 0; psi 0; np size(pareto, 1); for i 1:np % 取出当前前沿上两个相邻点围成的矩形区间 if i 1 lower [ref(1), pareto(1,2)]; else lower [pareto(i-1,1), pareto(i,2)]; end upper pareto(i, :); if lower(1) upper(1) || lower(2) upper(2) continue; end % 计算二维矩形内均值的累积概率 p (gausscdf(upper(1), mu(1), sigma(1)) - gausscdf(lower(1), mu(1), sigma(1))) ... * (gausscdf(upper(2), mu(2), sigma(2)) - gausscdf(lower(2), mu(2), sigma(2))); psi psi p; % 用中心点乘以概率近似贡献精细实现需用条件期望 mid 0.5 * (lower upper); exi exi p * (ref - mid); % 向量差值二维时求范数或矩形面积 end % 归一化避免坐标系影响 exi abs(sum(exi)); end这段代码在严格推导上并不是精确积分但它抓住了二维 HVEI 的几何本质把超体积增量拆成前沿点附近的若干个矩形条候选点落入这些矩形条的概率乘以面积增量就是期望贡献。gausscdf的两次差值给出了落入矩形区域的概率。实际使用时pareto必须预先按第一目标排序同时过滤掉被支配点否则矩形条会重叠概率被重复计算。3.3 HVEI 与单目标 EI 的权衡差异单目标 EI 在均值附近和标准差较大时都会偏高而 HVEI 还依赖候选点相对参考点的位置。远离参考点的区域即使均值不高但一旦能大幅提升超体积也会呈现较高采集值。因此 HVEI 更容易在探索早期选中极端目标区域后期则需要前沿点分布更密。这跟多目标算法里常用的拥挤距离淘汰策略正好互补。对比项单目标 EIHVEI采集函数输入最优单目标值整个Pareto前沿参考点不确定性处理高斯分布多维高斯联合分布计算复杂度O(1)O(npoints) 及以上探索倾向中等高尤其超体积缺口区域4. 运行与调试利用 Untitled.m 搭建多目标优化迭代框架4.1 主脚本 Untitled.m 的典型流程Untitled.m一看就是作者临时起名的主脚本。通常一个完整的 HVEI 多目标优化循环会这样组织先随机生成初始样本评估真实目标函数得到初始 Pareto 集然后循环(1) 用高斯过程回归拟合每个目标(2) 调用exipsi.m计算候选点的 HVEI(3) 用优化器找使 HVEI 最大的位置(4) 在该位置评估真实目标并更新 Pareto 集。下面是一段可运行的骨架% Untitled.m 骨架 clear; clc; rng(1); % 初始采样lhsdesign 来自 Statistics Toolbox X lhsdesign(10, 2); % 示例目标: 两个冲突函数 Y [X(:,1).^2 X(:,2), X(:,1) X(:,2).^2]; ref [1.2, 1.2]; for iter 1:30 % 非支配筛选non_dominated 需自行实现 nd non_dominated(Y); pareto_Y Y(nd, :); % 为每个目标独立拟合高斯过程模型 gp1 fitrgp(X, Y(:,1), KernelFunction, ardsquaredexponential); gp2 fitrgp(X, Y(:,2), KernelFunction, ardsquaredexponential); % 最大化HVEI: 通过最小化 -HVEI options optimoptions(fmincon, Display, off); lb [0,0]; ub [1,1]; x0 rand(1,2); [xbest, ~] fmincon((x)hvei_opt(x, gp1, gp2, pareto_Y, ref), ... x0, [], [], [], [], lb, ub, [], options); % 真实评估新候选点 ybest [xbest(1)^2 xbest(2), xbest(1) xbest(2)^2]; X [X; xbest]; Y [Y; ybest]; end function f hvei_opt(x, gp1, gp2, pareto_Y, ref) [mu1, s1] predict(gp1, x); [mu2, s2] predict(gp2, x); [exi, ~] exipsi([mu1, mu2], [s1, s2], pareto_Y, ref); f -exi; end这段代码里fitrgp和lhsdesign依赖 Statistics and Machine Learning Toolbox而包内核心函数并不依赖任何附加工具箱。predict的第二个输出是标准差不是方差s1、s2会随着x变化这正是采集函数需要的不确定性信息。fmincon来自 MATLAB 优化工具箱用来最大化hvei的负值由于 HVEI 往往非凸建议用MultiStart或GlobalSearch多起点搜索避免陷入局部最优。4.2 参数设置与调优建议HVEI 对几个参数非常敏感参考点ref、代理模型核函数、候选搜索初始点数量。下面是我习惯的初始值。参数建议值影响初始样本数2dim5dim太少则GP拟合不可靠参考点每个目标上界*1.1参考点过大HV差异不敏感GP核函数ARD Squared Exponential适应性最好fmincon起点数10防止只找到局部最优停止条件HV增量1e-4 或最大迭代控制过拟合提示exipsi的第二个输出psi可以用作停止条件当它连续 5 次低于 0.1% 时说明新点改进概率很小可以提前终止。4.3 调试中常见的维度与数值错误这个包最容易出错的是把pareto_Y与pareto混淆。exipsi.m里要求传入的是目标值空间的 Pareto 前沿不是决策变量。如果是做三维以上的问题exi2d.m会直接报错所以主脚本里要先判断目标数。另一个高频错误是标准差为零比如初始样本太少或者重复采样导致gausscdf出现 NaN。可以在调用前加一行sigma max(sigma, 1e-6)兜底。5. 收敛性评估与超体积计算实践5.1 记录每次迭代的超体积主循环里每更新一次 Pareto 集就应该记录一次超体积用来画收敛曲线。利用hvolume2d.m非常简单hv_hist zeros(maxIter, 1); for iter 1:maxIter % ... 迭代逻辑 ... hv_hist(iter) hvolume2d(pareto_Y, ref); end semilogy(1:maxIter, hv_hist); xlabel(迭代次数); ylabel(Hypervolume);如果hv_hist出现震荡先检查pareto_Y中是否混入了被支配点。hvolume2d不负责非支配检查它假定输入已经是非支配点。在多目标优化库里常见的做法是调用paretoset或自己写非支配排序再传给hvolume2d。5.2 对比 HVEI 与 NSGA-II、MOEA/D为了验证这个包的实际效果可以把 HVEI 作为采集函数嵌入贝叶斯优化与直接使用 NSGA-II、MOEA/D 做固定评估次数预算对比。下表展示了典型结果以 ZDT1 问题为例算法最终超体积均值需要真实评估次数运行时间秒NSGA-II0.7832006.2MOEA/D0.7912005.8HVEIGP0.836603.1HVEIGP0.8421005.0这里的数字只是示意实际效果取决于问题和代理模型精度。HVEI 的优势是真实评估次数少代价是每次选点都要优化一个复杂的采集函数耗时集中在 GP 拟合和exipsi数值积分上。在机器学习模型优化方案里超参搜索经常被建模成多目标问题用 HVEI 替代网格搜索可以明显减少训练次数。5.3 参考点选择对 HV 收敛的影响参考点决定了超体积度量的基准。参考点太远所有解的覆盖面积普遍较大难以区分优劣参考点太近部分前沿点可能被参考点截断。一个实用的动态策略是每一轮用当前Y的最大值乘以 1.1 作为参考点并让参考点在迭代过程中保持稳定避免曲线跳变。在带噪声的鲁棒优化模型里参考点必须位于可行域的可达区域之外否则超体积值会失真。6. 踩坑与边界情况高维扩展与数值稳定性调优6.1 数值稳定性处理gausscdf在参数极端时会返回 0 或 1导致exi2d中出现零概率区间。可以在exi2d里对概率加一个下限例如p max(p, eps); sigma max(sigma, 1e-6);sigma太小时gausspdf计算出的密度值可能超过realmax这时最好对sigma做下限裁剪或改用log形式计算后再取指数。这些看似小的改动在迭代几百次之后能显著减少 NaN 和 Inf 的传播。6.2 高维目标的粗粒度扩展原包里只有hvolume2d但你需要扩展三维超体积时不必从头实现。可以使用蒙特卡洛采样近似 HV然后在exipsi.m中把解析积分替换成样本均值每采样一个候选点预测值就计算一次HV(P ∪ {y}) - HV(P)最后平均。样本数建议取 1000 到 2000维度再高就配合拉丁超立方采样。注意采样方差会拖慢收敛所以高维场景下应优先考虑参考点裁剪和前沿精简。6.3 带约束和离散变量的处理对于带约束问题可以给不可行解的 HVEI 乘以惩罚系数或者在exipsi.m中把不可行域的目标值推到参考点之外。离散变量则需要对候选解做整数约束如果用fmincon就把整数变量四舍五入后再评估但这样会破坏 GP 的平滑性更稳妥的是用遗传算法做采集函数优化。本文还有配套的精品资源点击获取 HVEI多目标优化算法：超体积期望改进原理与MATLAB实现
简介面向多目标优化研究者与MATLAB开发者该资源提供基于超体积期望改进HVEI算法的MATLAB实现适用于在成本、效率、安全性等互相冲突的目标中获取Pareto最优解集的场景特别适合有一定优化基础、希望获得可运行代码的读者。压缩包共7个文件6个.m代码文件涵盖核心HVEI计算、二维超体积计算以及高斯分布函数等辅助模块并提供示例脚本演示调用流程另含1个txt说明文件讲解使用要点整体大小仅4KB轻量精简便于快速部署和调试。资源当前已有379人学习下载。读者借助代码可完整复现HVEI模型理解期望改进如何引导搜索最大化超体积覆盖同时可根据自身问题重新定义目标函数与约束调整种群初始化与迭代逻辑从而迁移到工程优化、能源调度等实际任务中。1. 从文件清单看 HVEI 优化模型的结构设计拿到HV_based_expected_improvement.zip解压以后就是六个.m文件和一份README.txt。Untitled.m是入口脚本exipsi.m、exi2d.m负责期望改进的计算hvolume2d.m算二维超体积gausspdf.m和gausscdf.m提供高斯分布的概率密度和累积分布函数。这个包解决的不是普通单目标优化而是多目标优化里最难的一类问题如何在不知道全局 Pareto 前沿的前提下用尽可能少的真实评估找到一组覆盖度好、分布均匀的非支配解。HVEI 的核心思想是把单目标贝叶斯优化里的 Expected Improvement 推广到超体积指标上每一步都选择对当前非支配集超体积增量期望最大的候选点在探索未知区域和利用已有优势之间取得平衡。如果你正在做多目标贝叶斯优化、代理模型辅助设计或者想在一个比较干净的 MATLAB 代码基础上改造出适合自己的采集函数这个包值得读一遍。它没有依赖额外工具箱核心逻辑都写在几个函数里适合边读边改。2. 超体积与期望改进HVEI 的数学原理和 MATLAB 实现基础2.1 超体积HV的定义与 hvolume2d.m 计算超体积是指由非支配解集与参考点围成的多维区域大小在多目标优化中它是唯一一个同时满足严格单调性、不受目标尺度影响只要参考点一致的纯量化指标。二维问题里超体积就是若干矩形面积并集。hvolume2d.m正是针对二维情况做了快速实现输入是前沿点集和参考点输出面积。常见的实现逻辑是function hv hvolume2d(points, ref) % POINTS: 每个点为一行目标值越小越好 % REF: 参考点例如 [1.1, 1.1] pts sortrows(points, 1); % 按第一维排序 hv 0; prev_y ref(2); for i 1:size(pts, 1) if pts(i,2) prev_y hv hv (pts(i,1) - ref(1)) * (prev_y - pts(i,2)); prev_y pts(i,2); end end hv abs(hv); end这里按第一维升序排列prev_y从参考点第二维开始逐步下移累计每个点的矩形条面积。注意如果目标是最小化那么参考点应当大于所有目标值而排序后第一个点贡献面积从参考点往左abs用于防止负面积。实际使用时要先做非支配排序把被支配的点过滤掉否则面积会偏大。对于三维以上二维扫描线就失效需要改用deldir或蒙特卡洛采样这也是这个包里只提供二维hvolume的原因。2.2 高斯密度与分布函数gausspdf.m 和 gausscdf.mHVEI 的期望运算离不开高斯分布。代理模型给出的预测不只是一个均值还有一个方差可以用正态分布刻画目标值的不确定性。gausspdf.m和gausscdf.m提供了标准高斯 pdf 和 cdf防止依赖工具箱版本。其内部实现等价于function y gausspdf(x, mu, sigma) if nargin 1, mu 0; sigma 1; end y exp(-0.5 * ((x - mu) ./ sigma).^2) ./ (sigma * sqrt(2 * pi)); end function p gausscdf(x, mu, sigma) if nargin 1, mu 0; sigma 1; end p 0.5 * (1 erf((x - mu) ./ (sigma * sqrt(2)))); end两个函数都支持向量输入参数mu和sigma与x同维度或标量。和 MATLAB 内置的normpdf、normcdf相比自定义版本少一次函数寻址开销在循环里调用更快。但数值上要注意sigma为 0 时会出现除零后续章节会给出处理建议。2.3 从单目标 EI 到多目标 HVEI单目标期望改进定义是EI(x) E[ max(0, f_best - f(x)) ]当预测分布为高斯时有解析表达式EI(x) (f_best - μ) * Φ((f_best - μ)/σ) σ * φ((f_best - μ)/σ)其中Φ是 cdfφ是 pdf。多目标场景下f_best变成当前 Pareto 前沿改进量变成“新点加入后超体积的提升量”。HVEI 的自然定义HVEI(x) E[ HV(P ∪ {f(x)}) - HV(P) ]但f(x)是随机向量各目标间可能有相关性通常先假设独立再利用联合分布做数值积分。这个包里的exipsi.m就是把连续积分转化为有限求和配合二维超体积的解析式得到较精确的 HVEI 近似值。函数文件作用依赖hvolume2d.m计算二维超体积排序gausspdf.m标准正态密度无gausscdf.m标准正态分布函数无exi2d.m二维HVEI积分gausspdf, gausscdfexipsi.mHVEI及改进概率近似exi2dUntitled.m主优化循环以上全部3. 核心函数拆解exipsi.m 和 exi2d.m 的 HVEI 计算流程3.1 exipsi.m 的签名与结构exipsi.m的名称来自 Expected Improvementexi和 Probability of Improvementpsi。在实际实现中它接收预测均值、标准差、当前 Pareto 集和参考点返回期望超体积增量和改进概率。框架可以这样写function [exi, psi] exipsi(mu, sigma, pareto, ref) if length(mu) ~ 2 % 高维问题退化为采样近似本包未实现时提示 error(当前版本只支持二维目标); end [exi, psi] exi2d(mu, sigma, pareto, ref); end从这个结构可以看出作者刻意保留了扩展接口二维走解析路线高维可以在这里替换成蒙特卡洛采样或数值积分。psi的实用价值在于给采集函数提供另一个参考当 HVEI 接近 0 但改进概率很高时说明该点性能波动大值得赌一把反之 HVEI 高而psi低说明期望收益来自极端情况实际风险大。3.2 exi2d.m 的二维积分实现二维 HVEI 计算的核心是对于任意候选点预测分布 (μ1,σ1) 和 (μ2,σ2)求出其目标向量落在当前前沿“超体积缺口”内的概率与贡献量的乘积积分。常见的做法是先由当前 Pareto 前沿生成一组不重叠矩形区域再对每个矩形区域做二维正态期望计算。伪代码级别的实现如下function [exi, psi] exi2d(mu, sigma, pareto, ref) exi 0; psi 0; np size(pareto, 1); for i 1:np % 取出当前前沿上两个相邻点围成的矩形区间 if i 1 lower [ref(1), pareto(1,2)]; else lower [pareto(i-1,1), pareto(i,2)]; end upper pareto(i, :); if lower(1) upper(1) || lower(2) upper(2) continue; end % 计算二维矩形内均值的累积概率 p (gausscdf(upper(1), mu(1), sigma(1)) - gausscdf(lower(1), mu(1), sigma(1))) ... * (gausscdf(upper(2), mu(2), sigma(2)) - gausscdf(lower(2), mu(2), sigma(2))); psi psi p; % 用中心点乘以概率近似贡献精细实现需用条件期望 mid 0.5 * (lower upper); exi exi p * (ref - mid); % 向量差值二维时求范数或矩形面积 end % 归一化避免坐标系影响 exi abs(sum(exi)); end这段代码在严格推导上并不是精确积分但它抓住了二维 HVEI 的几何本质把超体积增量拆成前沿点附近的若干个矩形条候选点落入这些矩形条的概率乘以面积增量就是期望贡献。gausscdf的两次差值给出了落入矩形区域的概率。实际使用时pareto必须预先按第一目标排序同时过滤掉被支配点否则矩形条会重叠概率被重复计算。3.3 HVEI 与单目标 EI 的权衡差异单目标 EI 在均值附近和标准差较大时都会偏高而 HVEI 还依赖候选点相对参考点的位置。远离参考点的区域即使均值不高但一旦能大幅提升超体积也会呈现较高采集值。因此 HVEI 更容易在探索早期选中极端目标区域后期则需要前沿点分布更密。这跟多目标算法里常用的拥挤距离淘汰策略正好互补。对比项单目标 EIHVEI采集函数输入最优单目标值整个Pareto前沿参考点不确定性处理高斯分布多维高斯联合分布计算复杂度O(1)O(npoints) 及以上探索倾向中等高尤其超体积缺口区域4. 运行与调试利用 Untitled.m 搭建多目标优化迭代框架4.1 主脚本 Untitled.m 的典型流程Untitled.m一看就是作者临时起名的主脚本。通常一个完整的 HVEI 多目标优化循环会这样组织先随机生成初始样本评估真实目标函数得到初始 Pareto 集然后循环(1) 用高斯过程回归拟合每个目标(2) 调用exipsi.m计算候选点的 HVEI(3) 用优化器找使 HVEI 最大的位置(4) 在该位置评估真实目标并更新 Pareto 集。下面是一段可运行的骨架% Untitled.m 骨架 clear; clc; rng(1); % 初始采样lhsdesign 来自 Statistics Toolbox X lhsdesign(10, 2); % 示例目标: 两个冲突函数 Y [X(:,1).^2 X(:,2), X(:,1) X(:,2).^2]; ref [1.2, 1.2]; for iter 1:30 % 非支配筛选non_dominated 需自行实现 nd non_dominated(Y); pareto_Y Y(nd, :); % 为每个目标独立拟合高斯过程模型 gp1 fitrgp(X, Y(:,1), KernelFunction, ardsquaredexponential); gp2 fitrgp(X, Y(:,2), KernelFunction, ardsquaredexponential); % 最大化HVEI: 通过最小化 -HVEI options optimoptions(fmincon, Display, off); lb [0,0]; ub [1,1]; x0 rand(1,2); [xbest, ~] fmincon((x)hvei_opt(x, gp1, gp2, pareto_Y, ref), ... x0, [], [], [], [], lb, ub, [], options); % 真实评估新候选点 ybest [xbest(1)^2 xbest(2), xbest(1) xbest(2)^2]; X [X; xbest]; Y [Y; ybest]; end function f hvei_opt(x, gp1, gp2, pareto_Y, ref) [mu1, s1] predict(gp1, x); [mu2, s2] predict(gp2, x); [exi, ~] exipsi([mu1, mu2], [s1, s2], pareto_Y, ref); f -exi; end这段代码里fitrgp和lhsdesign依赖 Statistics and Machine Learning Toolbox而包内核心函数并不依赖任何附加工具箱。predict的第二个输出是标准差不是方差s1、s2会随着x变化这正是采集函数需要的不确定性信息。fmincon来自 MATLAB 优化工具箱用来最大化hvei的负值由于 HVEI 往往非凸建议用MultiStart或GlobalSearch多起点搜索避免陷入局部最优。4.2 参数设置与调优建议HVEI 对几个参数非常敏感参考点ref、代理模型核函数、候选搜索初始点数量。下面是我习惯的初始值。参数建议值影响初始样本数2dim5dim太少则GP拟合不可靠参考点每个目标上界*1.1参考点过大HV差异不敏感GP核函数ARD Squared Exponential适应性最好fmincon起点数10防止只找到局部最优停止条件HV增量1e-4 或最大迭代控制过拟合提示exipsi的第二个输出psi可以用作停止条件当它连续 5 次低于 0.1% 时说明新点改进概率很小可以提前终止。4.3 调试中常见的维度与数值错误这个包最容易出错的是把pareto_Y与pareto混淆。exipsi.m里要求传入的是目标值空间的 Pareto 前沿不是决策变量。如果是做三维以上的问题exi2d.m会直接报错所以主脚本里要先判断目标数。另一个高频错误是标准差为零比如初始样本太少或者重复采样导致gausscdf出现 NaN。可以在调用前加一行sigma max(sigma, 1e-6)兜底。5. 收敛性评估与超体积计算实践5.1 记录每次迭代的超体积主循环里每更新一次 Pareto 集就应该记录一次超体积用来画收敛曲线。利用hvolume2d.m非常简单hv_hist zeros(maxIter, 1); for iter 1:maxIter % ... 迭代逻辑 ... hv_hist(iter) hvolume2d(pareto_Y, ref); end semilogy(1:maxIter, hv_hist); xlabel(迭代次数); ylabel(Hypervolume);如果hv_hist出现震荡先检查pareto_Y中是否混入了被支配点。hvolume2d不负责非支配检查它假定输入已经是非支配点。在多目标优化库里常见的做法是调用paretoset或自己写非支配排序再传给hvolume2d。5.2 对比 HVEI 与 NSGA-II、MOEA/D为了验证这个包的实际效果可以把 HVEI 作为采集函数嵌入贝叶斯优化与直接使用 NSGA-II、MOEA/D 做固定评估次数预算对比。下表展示了典型结果以 ZDT1 问题为例算法最终超体积均值需要真实评估次数运行时间秒NSGA-II0.7832006.2MOEA/D0.7912005.8HVEIGP0.836603.1HVEIGP0.8421005.0这里的数字只是示意实际效果取决于问题和代理模型精度。HVEI 的优势是真实评估次数少代价是每次选点都要优化一个复杂的采集函数耗时集中在 GP 拟合和exipsi数值积分上。在机器学习模型优化方案里超参搜索经常被建模成多目标问题用 HVEI 替代网格搜索可以明显减少训练次数。5.3 参考点选择对 HV 收敛的影响参考点决定了超体积度量的基准。参考点太远所有解的覆盖面积普遍较大难以区分优劣参考点太近部分前沿点可能被参考点截断。一个实用的动态策略是每一轮用当前Y的最大值乘以 1.1 作为参考点并让参考点在迭代过程中保持稳定避免曲线跳变。在带噪声的鲁棒优化模型里参考点必须位于可行域的可达区域之外否则超体积值会失真。6. 踩坑与边界情况高维扩展与数值稳定性调优6.1 数值稳定性处理gausscdf在参数极端时会返回 0 或 1导致exi2d中出现零概率区间。可以在exi2d里对概率加一个下限例如p max(p, eps); sigma max(sigma, 1e-6);sigma太小时gausspdf计算出的密度值可能超过realmax这时最好对sigma做下限裁剪或改用log形式计算后再取指数。这些看似小的改动在迭代几百次之后能显著减少 NaN 和 Inf 的传播。6.2 高维目标的粗粒度扩展原包里只有hvolume2d但你需要扩展三维超体积时不必从头实现。可以使用蒙特卡洛采样近似 HV然后在exipsi.m中把解析积分替换成样本均值每采样一个候选点预测值就计算一次HV(P ∪ {y}) - HV(P)最后平均。样本数建议取 1000 到 2000维度再高就配合拉丁超立方采样。注意采样方差会拖慢收敛所以高维场景下应优先考虑参考点裁剪和前沿精简。6.3 带约束和离散变量的处理对于带约束问题可以给不可行解的 HVEI 乘以惩罚系数或者在exipsi.m中把不可行域的目标值推到参考点之外。离散变量则需要对候选解做整数约束如果用fmincon就把整数变量四舍五入后再评估但这样会破坏 GP 的平滑性更稳妥的是用遗传算法做采集函数优化。本文还有配套的精品资源点击获取