NSGA-II算法求解多目标优化问题
多目标优化问题介绍
1.对于多余一个的目标函数在给定区域内的最优化问题称为多目标优化问题。
例如：在给定条件下，设计一款汽车，既要满足安全（重量大），又要满足经济（耗油量小）即为多目标最优化问题。
该模型通常可总述为：
其中x=(x1,x2,x3…xn)所在的空间Ω称为决策空间(可行解空间)，向量F(x)所在的空间为目标空间
不同于单目标优化，在多目标优化中，各目标函数之间是相互冲突的。导致不一定存在在所有目标函数上都是最优解，某个解可能在一个目标函数中是最优的，但在另一个目标函数中是最差(判断多目标优化问题的关键)
因此，多目标最优化问题得到的通常是一个解集，他们之间不能简单地比较优劣，这样的被称为非支配解（有效解），Pareto最优解。
2.多目标优化相关概念：
定义一
向量支配对于最小化问题，一个向量u=(u1,u2…un)支配另一个向量v=(v1,v2…vn)，当且仅当ui≤vi
定义二
如果Ω中没有支配(优于)x的解，则称x是该问题的一个Pareto最优解(非劣解，有效解，非支配解)
Pareto最优解的全体组成的集合称为Pareto最优解集，Pareto最优解集在目标函数空间的像集称为Pareto Front(Pareto 前沿)
如上图所示，A，B两点中是一个最优解，他们相应的自变量向量没有被其他任何的自变量向量所支配，也就意味着在任何一个目标上他们都不能被其他个体所支配。这样的解就是Pareto最优解。
由定义一和定义二可知C点优于D点，C点优于F点，C点优于E点
NSGA-II
NSGA-II（Non-dominated Sorting Genetic Algorithm II）是一种流行的多目标优化遗传算法，由Srinivas和Deb在2000年提出。它在NSGA的基础上进行了改进，具有运行速度快和解集收敛性好的优点，常被用作其他多目标优化算法性能的基准。
NSGA-II的核心算法主要由快速非支配排序和拥挤比较算子两部分组成。快速非支配排序算法用于将种群中的个体根据非支配关系进行分层排序，而拥挤比较算子则用于在同一层内保持种群的多样性。此外，NSGA-II还采用了锦标赛选择法作为选择策略，并引入了精英策略，即在父代和子代中选择最好的个体进入下一代，以保证优良种群个体在进化过程中不会被丢弃。
NSGA-II算法的主要特点包括：
- 非支配排序：算法首先将种群中的个体根据非支配关系进行排序，非支配关系意味着一个解在所有目标上都不比另一个解差。
- 拥挤距离：为了维持种群多样性，NSGA-II引入了拥挤距离的概念，这是一个度量，用于衡量一个解周围有多少空间可以用于生成新的解，从而避免解的聚集。
- 精英策略：在生成下一代种群时，NSGA-II使用精英策略，即在父代和子代中选择最好的个体进入下一代。
- 快速非支配排序：NSGA-II改进了非支配排序的算法，使其更加高效。
NSGA-II算法在多个领域有广泛的应用，包括电力系统、制造系统调度、热交换器设计、材料设计等。例如，在材料开发领域，NSGA-II可以用于优化材料的设计参数，以满足多个性能指标，如寻找具有最佳强度、韧性和成本效益的合金配方，或确定具有特定电导率、热导率和机械性能的复合材料的最佳制造过程。
总的来说，NSGA-II算法因其高效性和良好的解集收敛性，在多目标优化问题中得到了广泛的应用和研究。
使用MATLAB的gamultiobj工具箱求解NSGA-II算法
gamultiobj的使用范式
编写程序
清除所有变量（非必须，但注意变量不能和下面所用的冲突）
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
RUN
求解时间受求解器设置影响，可能会较长，请耐心等待
求解过程中会实时显示当前种群的情况：
如果已经达到满意，也可点击stop按钮提前结束求解
最后的求解结果，即Pareto最优解集储存在[x,fval]中，fval是x对应的目标函数值。fval大致构成了一条空间曲线——Pareto前沿。若各个解较为均匀分布，说明该图包含了大部分最优解情况，全局性优，适用性强。在满足Pareto最优的条件下，是没有办法在不让某一优化目标受损的情况下，令另一方目标获得更优的。所以这些解均为最优，对最优解的具体选择可以根据实际情况。
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
————————————————
例子2
- 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.
扫一扫
1185
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
