Matlab中多目标优化遗传算法求解问题方法

gamultiobj 跑了一晚上，Pareto 前沿出来两百多个点，但 matlab 画出来一团糊，两个目标函数的值域差了一个量级，横轴从0到1，纵轴从0到10000，点全挤在底部。回去分别归一化之后再画，前沿的弧线就出来了。D:\Matlab_work\multiobj\pareto_plot.m 里存着这个画图脚本，横轴用归一化后的 f1，纵轴用归一化后的 f2，plot(x(:,1), x(:,2), 'o')，grid on，一眼能看出拐点在哪。

调用格式先记清楚，x = gamultiobj(fun, nvars) 是最简单的形式，fun 是目标函数句柄，返回两个目标值，nvars 是变量个数。带线性约束的话往后排 A、b、Aeq、beq、lb、ub，非线性约束用 nonlcon 函数句柄。[x, fval] = gamultiobj(...) 多拿一个返回值，fval 是 Pareto 前沿上每个点对应的两个目标函数值，画图的时候直接用 fval 不用重新算目标函数，省一遍计算。D:\Matlab_work\multiobj\basic_call.m 里写了三种调用形式的对比。

选项用 optimoptions('gamultiobj', ...) 建，常用的几个：PopulationSize 默认50，变量多或者前沿复杂的时候给200到500，有人做到500才把前沿跑完整；MaxGenerations 默认100，目标函数是仿真的话一代几分钟，先给30代看趋势，曲线平了就停；ParetoFraction 控制前沿上保留的解的比例，默认0.35，想要更密的前沿把它加到0.5；PlotFcn 设成 @gaplotpareto，每代自动画一次前沿，跑的时候肉眼看着收不收敛。FunctionTolerance 默认1e-4，前沿的 spread 变化小于这个值求解器就停。MaxStallGenerations 默认100，连续100代没有改善才停，跑仿真的时候这个值给30到50省时间。

并行计算对 gamultiobj 是刚需，每代要评估整个种群，串行跑几十代根本受不了。optimoptions('gamultiobj', 'UseParallel', true, 'UseVectorized', false)，UseParallel 设成 true 或者 'auto'，UseVectorized 必须 false，gamultiobj 靠 parfor 把种群拆到多个 worker 上算，UseVectorized 设成 true 的话求解器用一次函数调用评估整个种群，不走 parfor，并行就失效了。启动之前 parpool 开好并行池，或者让 UseParallel 自动开。delete(gcp) 先清掉现有的池子，避免残留的 worker 干扰。并行有个副作用记一下，并行种群的随机数序列不可控，每次跑出来的 Pareto 前沿都不完全一样，做对比实验的时候要么固定串行跑，要么接受这种随机性。UseParallel 设成 'auto' 跟 true 的区别不大，auto 是求解器自己决定用不用。

非线性约束用 nonlcon 函数句柄传进去，函数返回两个行向量：ineqnonlin(x) 和 eqnonlin(x)，前者是不等式约束 c(x)≤0，后者是等式约束 ceq(x)=0。x = gamultiobj(@myfun, nvars, [], [], [], [], lb, ub, @mycon, options)，@mycon 放在第八个参数的位置。gamultiobj 的非线性约束算法固定是 'penalty'，改不了，用惩罚函数把约束违反量加到目标函数上，罚因子迭代调整。等式约束尽量少用，gamultiobj 对等式约束的处理不如 ga 好，能转成不等式就转。整数变量 gamultiobj 不支持，官方文档里明确写了“gamultiobj does not accept integer constraints”，有整数变量的话要么先把连续解算出来再取整，要么换 paretosearch 或者自己写 NSGA-II。

从 Pareto 前沿上选点的时候别只看两端的极端解，两端的解通常一个目标极好另一个目标极差，工程上没法用。中间那段曲率最大的地方就是拐点，在拐点附近选，两个目标的平衡最好。选点的方法是在前沿上算每个点到原点的距离，归一化之后找距离最小的那个点，或者算相邻点之间的斜率变化，斜率变化最大的位置就是拐点。D:\Matlab_work\multiobj\knee_point.m 里写了一个简单的拐点选取脚本，算的是归一化后前沿上每个点的曲率，返回曲率最大的那个点的索引。选出来的点回仿真验证一遍，约束违反量和目标函数值都对上再确认。

下次把那个支架轻量化的案例跑一遍，目标函数是质量和一阶频率，变量是三个板的厚度，lb 给 1mm ub 给 5mm，非线性约束是最大应力不超过 200MPa，那个 stress_con 函数还没写，应力从 Abaqus 里读的路径是……

免责声明：本文系网络转载或改编，未找到原创作者，版权归原作者所有。如涉及版权，请联系删