【MATLAB】多目标优化算法 NSGA-II(gamultiobj) 的使用精解
原始博文因为写的比较潦草,评论中有疑问的网友较多,所以重新写了一下 2021-4-24
增加了一些说明与参考文献,修改了几处笔误 2021-5-20
对于多目标优化(multiobjective optimization)算法 NSGA-II 实现的细节与原理不在此说明。感兴趣的读者可另行查阅
清除所有变量(非必须,但注意变量不能和下面所用的冲突)
clear
- 需求解模型的参数设置部分:(模型导入)
%% 模型设置
% 适应度函数的函数句柄
fitnessfcn=@Fun;
% 变量个数
nvars=4;
% 约束条件形式1:下限与上限(若无取空数组[])
% lb<=X <=ub
lb=[0,0,0,0];
ub=[];
% 约束条件形式2:线性不等式约束(若无取空数组[])
% A*X <=b
A = [0 0 1 1
-1/3 0 0 0
0 -1/2 0 0
0 0 0 0];
b = [48 ; 30 ; 30 ; 0];
% 约束条件形式3:线性等式约束(若无取空数组[])
% Aeq*X==beq
Aeq=[1 1 0 0;0 0 0 0; 0 0 0 0; 0 0 0 0];
beq=[120;0;0;0];
目标函数:(这一段需放在脚本最后或单独放在一个文件里)
function y=Fun(x)
% y是目标函数向量。有几个目标函数y就有多少个维度(数组y的长度)
% 因为gamultiobj是以目标函数分量取极小值为目标,
% 因此有些取极大值的目标函数注意取相反数
y(1)=-(x(1)*100/3 + x(3)*90/3 + x(2)*80/2+x(4)*70/2);
y(2)=x(3)+x(4);
end
gamultiobj
求解器设置部分:
%% 求解器设置
% 最优个体系数paretoFraction
% 种群大小populationsize
% 最大进化代数generations
% 停止代数stallGenLimit
% 适应度函数偏差TolFun
% 函数gaplotpareto:绘制Pareto前沿
options=gaoptimset('paretoFraction',0.3,'populationsize',200,'generations',300,'stallGenLimit',200,'TolFun',1e-10,'PlotFcns',@gaplotpareto);
gamultiobj
求解与结果输出部分:
%% 主求解
[x,fval]=gamultiobj(fitnessfcn,nvars,A,b,Aeq,beq,lb,ub,options)
%% 结果提取
% 因为gamultiobj是以目标函数分量取极小值为目标,
% 因此在y=Fun(x)里取相反数的目标函数再取相反数画出原始情况
plot(-fval(:,1),fval(:,2),'pr')
xlabel('f_1(x)')
ylabel('f_2(x)')
title('Pareto front')
grid on
求解时间受求解器设置影响,可能会较长,请耐心等待
求解过程中会实时显示当前种群的情况:
如果已经达到满意,也可点击stop
按钮提前结束求解
最后的求解结果,即 Pareto 最优解集储存在[x,fval]
中,fval
是x
对应的目标函数值。fval
大致构成了一条空间曲线——Pareto 前沿。若各个解较为均匀分布,说明该图包含了大部分最优解情况,全局性优,适用性强。在满足 Pareto 最优的条件下,是没有办法在不让某一优化目标受损的情况下,令另一方目标获得更优的。所以这些解均为最优,对最优解的具体选择可以根据实际情况。
表 1 工厂产品生产规格表
| 产品 | 生产时间(h/公斤) | 利润(元/公斤) | 加班时利润(元/公斤) |
| ---- | ------------------ | --------------- | --------------------- |
| A | 3 | 100 | 90 |
| B | 2 | 80 | 70 |
设工厂每周生产产品 A、B 的常规生产时长为、(h),加班生产时长为、 (h)。令 。设每周的利润函数为,加班时长函数为。
则目标函数为:
约束条件为:
因此,数学模型可以归纳为:
MATLAB 求解如下:
clear
clc
fitnessfcn=@Fun;
% 变量个数
nvars=4;
% lb<=X <=ub
lb=[0,0,0,0];
ub=[270 240 460 130];
% A*X <=b
A = [[13 13.5 14 11.5]
-[13 13.5 14 11.5]
0 0 -1 0
[0.015 0.02 0.018 0.011]];
b = [300*48 ; -300*40 ; -150 ; 20];
% Aeq*X=beq
Aeq=[];beq=[];
%最优个体系数paretoFraction
%种群大小populationsize
%最大进化代数generations
%停止代数stallGenLimit
%适应度函数偏差TolFun
%函数gaplotpareto:绘制Pareto前沿
options=gaoptimset('paretoFraction',0.3,'populationsize',200,'generations',300,'stallGenLimit',200,'TolFun',1e-10,'PlotFcns',@gaplotpareto);
[x,fval]=gamultiobj(fitnessfcn,nvars,A,b,Aeq,beq,lb,ub,options)
plot(-fval(:,1),fval(:,2),'pr')
xlabel('f_1(x)')
ylabel('f_2(x)')
title('Pareto front')
grid on
function y=Fun(x)
b = [270 240 460 130];
c = [300 300 600 200];
t = [190 210 148 100];
s = [200 230 160 114];
a = [0.015 0.02 0.018 0.011];
d = [13 13.5 14 11.5];
y(1)=-sum(x.*(s-t));
y(2)=sum(a.*x);
end
得到结果:
x1 x2 x3 x4 如下:
119.231391258967 0.769608488165712 0 0
119.231391258967 0.769608488165712 0 0
0.000499510291359209 120.000482867813 13.1757491344817 34.8252467588411
71.3391090591218 48.6614868846125 3.11344686170493 9.36001224382937
…… …… …… ……
27.0549008917871 92.9459248170927 9.47711506311976 25.3969178809142
8.65187243477257 111.349067997015 13.3680683558073 31.7052095195761
用分别表示 A、B、C、D 四种产品,表示第 i 种产品的产量(kg)。设最大产量为,销售量为,成本为,售价为,能耗为,生产时间为。设该问题的利润函数为 ,能耗函数为。 则利润函数为:
能耗函数为:
- 建立模型如下
- MATLAB 求解如下
% 清除所有变量(非必须)
clear
%% 模型设置
% 获取目标函数的函数句柄
fitnessfcn=@Fun;
% 变量个数
nvars=4;
% 约束条件形式1:(若无取空数组[])
% lb<=X <=ub
lb=[0,0,0,0];
ub=[];
% 约束条件形式2:(若无取空数组[])
% A*X <=b
A = [0 0 1 1
-1/3 0 0 0
0 -1/2 0 0
0 0 0 0];
b = [48 ; 30 ; 30 ; 0];
% 约束条件形式3:(若无取空数组[])
% Aeq*X=beq
Aeq=[1 1 0 0;0 0 0 0; 0 0 0 0; 0 0 0 0];
beq=[120;0;0;0];
%% 求解器设置
% 最优个体系数paretoFraction
% 种群大小populationsize
% 最大进化代数generations
% 停止代数stallGenLimit
% 适应度函数偏差TolFun
% 函数gaplotpareto:绘制Pareto前沿
options=gaoptimset('paretoFraction',0.3,'populationsize',200,'generations',300,'stallGenLimit',200,'TolFun',1e-10,'PlotFcns',@gaplotpareto);
%% 主求解
[x,fval]=gamultiobj(fitnessfcn,nvars,A,b,Aeq,beq,lb,ub,options)
%% 结果提取
% 因为gamultiobj是以目标函数分量取极小值为目标,
% 因此在y=Fun(x)里取相反数的目标函数再取相反数画出原始情况
plot(-fval(:,1),fval(:,2),'pr')
xlabel('f_1(x)')
ylabel('f_2(x)')
title('Pareto front')
grid on
function y=Fun(x)
% y是目标函数向量。有几个目标函数y就有多少个维度(数组y的长度)
% 因为gamultiobj是以目标函数分量取极小值为目标,
% 因此有些取极大值的目标函数注意取相反数
y(1)=-(x(1)*100/3 + x(3)*90/3 + x(2)*80/2+x(4)*70/2);
y(2)=x(3)+x(4);
end
求解结果为:
x1 x2 x3 x4 如下:
257.911499184609 147.920309053797 368.392357989384 129.803238959239
204.527452370415 215.831670926692 376.135110383218 129.541339137201
251.942563886570 239.988410149935 456.713154231118 129.635721179650
…… …… …… ……
245.261897051381 238.784443203755 429.044675830294 129.548264747046
216.531507460989 201.737471951873 368.856283121162 129.726418090506
参考文献:
K. Deb, S. Agrawal, A. Pratap, T. Meyarivan. A fast and elitist multiobjective genetic algorithm: NSGA-II. IEEE Trans. Evol. Comput. 2002 6(2): 182-197.
本文使用 Zhihu On VSCode 创作并发布