Matlab布谷鸟算法成本、时间、质量为目标的多目标优化求解代码既要省钱、又要赶工、还要交活不被挑毛病——这怕不是所有做工程优化、产品参数调优、甚至是毕设卡壳多目标建模的朋友都听过的老板/导师版“灵魂三连问”。之前试过用NSGA-II调遗传参数调到头秃今天换个更佛系更有生物感的算法多目标布谷鸟算法MOCS就拿最经典的三目标——成本、时间、质量——当小白鼠揉一揉。先别着急看代码先唠五毛钱布谷鸟的核心“操作”毕竟算法都是原理顺了改起来才敢改。单目标的布谷鸟就靠俩事 Lévy飞行找新窝随机但有长有短的跳能跳出局部最优这点比高斯散步强、帕累托淘汰劣质窝蛋随机选主人窝把自己生的新解和原窝蛋旧解比主人看不惯就扔劣质的还有1-pa的概率保留优质蛋pa一般设成0.25左右就行别太挑剔也别太随便。Matlab布谷鸟算法成本、时间、质量为目标的多目标优化求解代码那多目标怎么办无非是把“谁更好”的标准从“单数值大/小”换成“帕累托支配关系”解A支配解B就是所有目标都不比B差至少有一个比B强。剩下的就是把布谷鸟的淘汰逻辑改成先建个外部存档Archive存当前找到的所有非支配解然后每生成新解要么进Archive替换掉被它支配的旧解要么如果没人能支配它就加进去同时控制Archive的大小不然存太多内存炸了。废话够多了直接上可直接跑的简化版成本-时间-质量三目标MOCS代码我就拿某车间生产线的参数调优当背景编个目标函数参数比如是机器数量、工人加班时长、原材料检验率这些编得真实一点但别太复杂免得看代码分心。%% 参数预设 clear; clc; close all; % 布谷鸟自身参数 n_pop 30; % 种群大小30-50差不多够用 max_iter 200; % 最大迭代次数看收敛情况三目标可能要多跑点 p_a 0.25; % 蛋被主人发现的概率 dim 5; % 决策变量维度机器数(1-10整数)、加班时长(0-20h)、检验率(50%-100%)、零件损耗率阈值(1%-5%)、换线时间(5-30min) lb [1, 0, 0.5, 0.01, 5]; % 决策变量下界 ub [10, 20, 1, 0.05, 30]; % 决策变量上界 archive_size 50; % 外部存档最大容量 % 目标函数参数随便编的大家可以换成自己的问题 fixed_cost 10000; % 固定成本 machine_unit_cost 5000; % 每台机器的成本 overtime_unit_cost 50; % 每小时加班工资 inspect_unit_cost 0.1; % 每单位原材料检验费假设一批10w件这里系数调小 defect_rework_cost 20; % 每件缺陷品返工费 setup_unit_time_cost 100;% 每分钟换线时间对应的机会成本 base_time 240; % 无加班、无换线假设的基础生产时间 quality_coef1 0.2; % 检验率对质量的正影响系数 quality_coef2 0.3; % 加班时长对质量的负影响系数太累容易出活糙 quality_coef3 0.1; % 换线时间对质量的负影响系数频繁换线参数不稳 %% 初始化种群和外部存档 pop zeros(n_pop, dim); for i 1:n_pop pop(i,1) randi([lb(1), ub(1)],1); % 机器数是整数单独处理 pop(i,2:end) lb(2:end) (ub(2:end)-lb(2:end)).*rand(1, dim-1); end % 计算初始种群的三个目标值目标1最小化成本目标2最小化时间目标3最大化质量 fitness zeros(n_pop, 3); for i 1:n_pop x pop(i,:); % 目标1总成本 固定 机器 加班 检验 返工返工率和质量负相关 cost fixed_cost machine_unit_cost*x(1) overtime_unit_cost*x(2) ... inspect_unit_cost*100000*x(3) defect_rework_cost*100000*(1 - (quality_coef1*x(3) - quality_coef2*x(2)/20 - quality_coef3*x(30)/30 0.5)/1.1); % 哦对质量得先算个0-1之间的数不然返工率要炸刚才的返工率里是临时瞎调的不如单独拎出来 quality quality_coef1*x(3) - quality_coef2*(x(2)/ub(2)) - quality_coef3*(x(5)/ub(5)) 0.5; quality max(0, min(1, quality)); % 卡到0-1 cost fixed_cost machine_unit_cost*x(1) overtime_unit_cost*x(2) ... inspect_unit_cost*100000*x(3) defect_rework_cost*100000*(1 - quality); % 目标2总时间 基础 加班 换线时间假设有机器数*2次换线随便编个小公式 time base_time x(2) x(5)*(x(1)/2); % 目标3最大化质量后面帕累托支配如果统一成最小化会方便所以转成负的 fitness(i,:) [cost, time, -quality]; end % 初始化外部存档只存非支配解 archive []; archive_fitness []; for i 1:n_pop current_fit fitness(i,:); dominated false; % 先看当前存档里有没有支配这个解的 for j 1:size(archive_fitness,1) if all(archive_fitness(j,:) current_fit) any(archive_fitness(j,:) current_fit) dominated true; break; end end if ~dominated % 如果没有就删存档里被当前解支配的 to_keep true(size(archive_fitness,1),1); for j 1:size(archive_fitness,1) if all(current_fit archive_fitness(j,:)) any(current_fit archive_fitness(j,:)) to_keep(j) false; end end archive archive(to_keep,:); archive_fitness archive_fitness(to_keep,:); % 加进去 archive [archive; pop(i,:)]; archive_fitness [archive_fitness; current_fit]; end end %% 迭代优化 for iter 1:max_iter % Lévy飞行生成新解 new_pop pop; for i 1:n_pop % Lévy飞行的经典公式s u/|v|^(1/β)β一般取1.5 beta 1.5; sigma_u (gamma(1beta)*sin(pi*beta/2)/(gamma((1beta)/2)*beta*2^((beta-1)/2)))^(1/beta); u randn(1, dim)*sigma_u; v randn(1, dim); step u./abs(v).^(1/beta); % 步长缩放不然跳出决策变量范围太狠 step_size 0.01*(step).*(new_pop(i,:)-archive(randi(size(archive,1)),:)); % 加个存档解的引导让搜索更有方向 new_x new_pop(i,:) step_size; % 边界处理机器数取最近整数其他夹逼 new_x(1) round(new_x(1)); new_x max(lb, min(ub, new_x)); % 计算新解的适应度 new_quality quality_coef1*new_x(3) - quality_coef2*(new_x(2)/ub(2)) - quality_coef3*(new_x(5)/ub(5)) 0.5; new_quality max(0, min(1, new_quality)); new_cost fixed_cost machine_unit_cost*new_x(1) overtime_unit_cost*new_x(2) ... inspect_unit_cost*100000*new_x(3) defect_rework_cost*100000*(1 - new_quality); new_time base_time new_x(2) new_x(5)*(new_x(1)/2); new_fit [new_cost, new_time, -new_quality]; % 第一步淘汰和旧解比新解支配旧解就替换或者都进候选或者随机这里简化成支配就换不然保留旧解 if all(new_fit fitness(i,:)) any(new_fit fitness(i,:)) new_pop(i,:) new_x; fitness(i,:) new_fit; end end % 第二步淘汰蛋被主人发现随机替换p_a比例的解这里用高斯散步随机换 n_replace round(p_a*n_pop); replace_idx randperm(n_pop, n_replace); for i replace_idx % 随机选两个存档里的解或者随机换 new_x pop(i,:) 0.1*(ub - lb).*randn(1, dim); new_x(1) round(new_x(1)); new_x max(lb, min(ub, new_x)); new_quality quality_coef1*new_x(3) - quality_coef2*(new_x(2)/ub(2)) - quality_coef3*(new_x(5)/ub(5)) 0.5; new_quality max(0, min(1, new_quality)); new_cost fixed_cost machine_unit_cost*new_x(1) overtime_unit_cost*new_x(2) ... inspect_unit_cost*100000*new_x(3) defect_rework_cost*100000*(1 - new_quality); new_time base_time new_x(2) new_x(5)*(new_x(1)/2); new_fit [new_cost, new_time, -new_quality]; % 还是支配就换 if all(new_fit fitness(i,:)) any(new_fit fitness(i,:)) pop(i,:) new_x; fitness(i,:) new_fit; end end pop new_pop; % 更新完两步替换同步主种群 % 更新外部存档 for i 1:n_pop current_fit fitness(i,:); dominated false; for j 1:size(archive_fitness,1) if all(archive_fitness(j,:) current_fit) any(archive_fitness(j,:) current_fit) dominated true; break; end end if ~dominated to_keep true(size(archive_fitness,1),1); for j 1:size(archive_fitness,1) if all(current_fit archive_fitness(j,:)) any(current_fit archive_fitness(j,:)) to_keep(j) false; end end archive archive(to_keep,:); archive_fitness archive_fitness(to_keep,:); % 加之前检查一下存档大小如果超了删掉拥挤度最小的这个很重要不然存档会全是集中在某几个角落的解帕累托前沿不好看 if size(archive_fitness,1) archive_size archive [archive; pop(i,:)]; archive_fitness [archive_fitness; current_fit]; else % 计算拥挤度 crowding zeros(size(archive_fitness,1),1); for obj 1:3 [sorted_fit, sorted_idx] sort(archive_fitness(:,obj)); crowding(sorted_idx(1)) inf; crowding(sorted_idx(end)) inf; obj_range sorted_fit(end) - sorted_fit(1); if obj_range 0 for j 2:size(archive_fitness,1)-1 crowding(sorted_idx(j)) crowding(sorted_idx(j)) (sorted_fit(j1) - sorted_fit(j-1))/obj_range; end end end % 删掉最小的那个 [~, min_crowd_idx] min(crowding); archive(min_crowd_idx,:) []; archive_fitness(min_crowd_idx,:) []; archive [archive; pop(i,:)]; archive_fitness [archive_fitness; current_fit]; end end end % 每迭代10次打印一下进度和当前存档大小 if mod(iter,10) 0 fprintf(迭代次数%d当前非支配解数量%d\n, iter, size(archive_fitness,1)); end end %% 结果可视化三目标帕累托前沿用3D散点图 figure(Color,w,Position,[100,100,800,600]); scatter3(archive_fitness(:,1), archive_fitness(:,2), -archive_fitness(:,3), 100, filled, MarkerFaceAlpha,0.8); xlabel(总成本元,FontSize,12,FontWeight,bold); ylabel(总时间小时,FontSize,12,FontWeight,bold); zlabel(产品质量0-1,FontSize,12,FontWeight,bold); title(三目标MOCS帕累托前沿,FontSize,14,FontWeight,bold); grid on; view(30,30); % 调整一下视角方便看 colorbar; % 颜色可以代表迭代最后加进去的顺序或者其他这里随便放这段代码里有几个不是最“标准”但新手友好的修改说一下我为啥这么改Lévy飞行加了存档引导原版是直接用当前解和全局解不对原版单目标有时候用的是全局最优但多目标没有全局最优所以拉个随机的存档解当“指南针”能让搜索更快往帕累托前沿靠比瞎跳强。边界处理单独拎了机器数很多新手把所有变量都当连续结果优化出来机器数是2.7台老板看了要打人这里直接round卡上下界再取整完美。质量目标转成负的统一最小化帕累托支配要么全最小要么全最大全最小写代码最顺不用来回换if里的符号省得出错。外部存档加了拥挤度删除这个必须加不然你跑200次迭代存档里全是成本最低但时间最长质量最差的或者反过来的几个极端解中间的过渡解全被挤没了没法做决策。拥挤度的逻辑很简单每个目标方向上把解排个序最左边最右边的拥挤度设成无穷大必须留边界解重要中间的解拥挤度等于左右两个解在这个目标上的差除以总范围最后把所有目标的拥挤度加起来越小说明这个解周围越挤删了不可惜。跑出来的结果大概是什么样3D散点图上会有一条或者一片连续的“曲面”曲面左边的点成本低、质量差、时间可能长可能短曲面右边的点成本高、质量好、时间可能短机器多、检验率高、换线时间短曲面中间的点就是老板/导师想要的“折中方案”比如成本在XX万时间在XX小时质量在95%左右你可以挑几个离原点近统一最小化嘛原点虽然可能不存在但方向是对的的解给他们选。最后补一句这段代码是简化版的通用框架大家可以直接把目标函数那块换成自己的三目标或者多目标模型决策变量如果全是连续的就把单独处理机器数的部分删掉如果全是整数的话Lévy飞行之后可以全round高斯散步也可以改成随机整数替换灵活得很。