(转载)基于遗传算法的多目标优化算法(matlab实现)

1.3 函数gamultiobj中的一些基本概念

        (2)在命令行调用GUI：

gamultiobj是使用遗传算法解决多目标优化问题的函数。它试图解决以下形式的多目标问题：
min F（X） subject to: AX <= b, AeqX = beq (线性约束)
X lb <= X <= ub (边界约束)

X = gamultiobj（FITNESSFCN，NVARS）找到在FITNESSFCN中定义的目标函数的局部帕累托集X。NVARS是优化问题的维数（决策变量的数量）。FITNESSFCN接受大小为1-by-NVARS的向量X，并返回在决策变量处评估的大小为1-by-numberOfObjectives的向量。X是具有NVARS列的矩阵。X中的行数与帕累托解的数量相同。帕累托集中的所有解同样优，取决于应用程序的设计者选择帕累托集中的解。

X = gamultiobj（FITNESSFCN，NVARS，A，b）找到在FITNESSFCN中定义的目标函数的局部帕累托集X，受制于线性不等式A*X <= B。仅支持默认的PopulationType选项（'doubleVector'）的线性约束；不支持其他种群类型，例如'bitString'和'custom'。

X = gamultiobj（FITNESSFCN，NVARS，A，b，Aeq，beq）找到在FITNESSFCN中定义的目标函数的局部帕累托集X，受制于线性等式AeqX = beq以及线性不等式AX <= b。（如果不存在不等式，则设置A=[]和b=[]。）仅支持默认的PopulationType选项（'doubleVector'）的线性约束；不支持其他种群类型，例如'bitString'和'custom'。

X = gamultiobj（FITNESSFCN，NVARS，A，b，Aeq，beq，lb，ub）定义设计变量X的下限和上限，以便在范围lb <= X <= ub中找到局部帕累托集。如果不存在边界，则使用空矩阵。如果X（i）无下限，则设置lb（i）= -Inf；如果X（i）无上限，则设置ub（i）= Inf。仅支持默认的PopulationType选项（'doubleVector'）的边界约束；不支持其他种群类型，例如'bitString'和'custom'。

X = gamultiobj（FITNESSFCN，NVARS，A，b，Aeq，beq，lb，ub，NONLCON）找到满足NONLCON中定义的约束的局部帕累托集X。函数NONLCON接受X并返回向量C和Ceq，分别表示非线性不等式和等式。帕累托集X必须满足C（X）<=0和Ceq（X）=0。如果不存在边界，则设置lb=[]和/或ub=[]。当PopulationType选项设置为'bitString'或'custom'时，不满足非线性约束。有关详细信息，请参见文档。

X = gamultiobj（FITNESSFCN，NVARS，A，b，Aeq，beq，lb，ub，NONLCON，options）使用OPTIONS中的值替换默认优化参数，找到帕累托集X。OPTIONS可以通过OPTIMOPTIONS函数创建。有关详细信息，请参见OPTIMOPTIONS。有关gamultiobj接受的选项列表，请参见文档。

X = gamultiobj（PROBLEM）查找PROBLEM的最小值。PROBLEM是具有以下字段的结构：
fitnessfcn：
nvars：
Aineq：<不等式约束的A矩阵>
bineq：<不等式约束的b向量>
Aeq：<等式约束的Aeq矩阵>
beq：<等式约束的beq向量>
lb：
ub：
nonlcon：<非线性约束函数>
options：<使用optimoptions（'gamultiobj'，...）创建的选项>
solver：
rngstate：<随机数生成器的状态>

[X，FVAL] = gamultiobj（FITNESSFCN，NVARS，...）还返回矩阵FVAL，FITNESSFCN中定义的所有目标函数在X的所有解中的值。FVAL具有numberOfObjectives列，行数与X相同。

[X，FVAL，EXITFLAG] = gamultiobj（FITNESSFCN，NVARS，...）还返回描述gamultiobj的退出条件的EXITFLAG。EXITFLAG的可能值以及相应的退出条件为

1：在options.MaxStallGenerations代内帕累托集的扩展的平均值的值的变化小于options.FunctionTolerance。
0：超过最大代数。
-1：由输出或绘图函数终止优化。
-2：未找到可行点。
-5：超时。

[X，FVAL，EXITFLAG，OUTPUT] = gamultiobj（FITNESSFCN，NVARS，...）还返回带有以下字段的结构OUTPUT：
rngstate：
generations：<总代数，不包括HybridFcn迭代>
funccount：<函数评估的总数>
maxconstraint：<最大约束违规（如果有）>
message：

[X，FVAL，EXITFLAG，OUTPUT，POPULATION] = gamultiobj（FITNESSFCN，...）在终止时还返回最终POPULATION。

[X，FVAL，EXITFLAG，OUTPUT，POPULATION，SCORE] = gamultiobj（FITNESSFCN，...）在终止时还返回最终POPULATION的SCORE。

3 MATLAB程序实现

3.1 gamultiobj组织结构

3.2 函数stepgamultiobj分析

1.函数stepgamultiobj结构及图形描述

        函数stepgamultiobj的代码位于：MATLAB安装目录\toolbox\gads\gads\private,其结构如图3所示。

function parents = selectiontournament(expectation,nParents,options,tournamentSize)
 
if nargin < 4 || isempty(tournamentSize)
    tournamentSize = 4;
end
 
% Choose the players
playerlist = ceil(size(expectation,1) * rand(nParents,tournamentSize));
% Play tournament
parents = tournament(playerlist,expectation);
 
function champions = tournament(playerlist,expectation)
%tournament between players based on their expectation
 
playerSize = size(playerlist,1);
champions = zeros(1,playerSize);
% For each set of players
for i = 1:playerSize
    players = playerlist(i,:);
    % For each tournament
    winner = players(1); % Assume that the first player is the winner
    for j = 2:length(players) % Winner plays against each other consecutively
        score1 = expectation(winner,:);
        score2 = expectation(players(j),:);
        if score2(1) > score1(1) 
            winner = players(j);
        elseif score2(1) == score1(1)
            try % socre(2) may not be present for single objective problems
                if score2(2) > score1(2)
                    winner = players(j);
                end
            catch
            end
        end
    end
    champions(i) = winner;
end

3.交叉、变异、产生子种群和父子种群合并

        函数gamultiobj默认的交叉函数和变异函数分别为函数crossoverintermediate和函数 mutationadaptfeasible,实现的功能与函数ga是一致的，这里不再展开。需要说明的是，由于函数gamultiobj中使用了支配和排序的思想，其精英是自动保留的，因此，不再具有函数ga中的精英保留操作，父子种群的合并通过函数stepgamultiobj中的以下语句实现：

population=[population;nextPopulation];

4.计算序值和拥挤距离

        函数rankAndDistance依次调用函数nonDominatedRank、函数DistanceMeasureFcn和函数trimPopulation,下面分别介绍。

(1)函数nonDominatedRank

        非支配排序函数nonDominatedRank的作用是对父、子种群合并后的种群中的个体进行排序。非支配排序函数nonDominatedRank的基本思想是：序值rankCounter从1开始，依次加1,在每一轮rankCounter中，依次将种群中未被排序(判断个体是否已被排序的矩阵每轮更新一次)的个体p与其余所有未被排序的个体q进行比较，检查个体q是否支配个体p,若否，则个体p被赋予当前序值rankCounter;反之，因为个体p受个体q支配，故个体p的序值高于当前rankCounter,应参与下一轮的排序。通过排序，种群中的所有个体被分到了不同的前端。接着进行的是前端中的拥挤距离计算。

(2)函数DistanceMeasureFcn

        拥挤距离计算函数DistanceMeasureFcn的作用是计算某一前端内每个个体与其相邻个体的距离，与函数nonDominatedRank计算出的序值一起，为函数selectiontournament的选择提供依据，同时，也为接下来的函数trimPopulation做好准备。默认的拥挤距离计算函数是distancecrowding,

        拥挤距离计算函数DistanceMeasureFcn的基本思想是：对多目标中的每一个目标，分别计算一次相应的拥挤距离，再将这些拥挤距离相加得到最后的拥挤距离。在每个目标对应的拥挤距离的计算中，前端两头的两个个体的拥挤距离为inf,其余个体的拥挤距离为与该个体相邻的前后两个个体在(一1,1)区间映射后的目标函数值之差，这里的相邻是指同一前端中个体间的目标函数值大小接近。显然，某个体的拥挤距离越大，表示该个体与相邻个体的目标函数值差别越大，也就是多样性越好，故在序值相同的条件下，在函数selectiontournament中就越应该被选中，在接下来的函数rimPopulation中就越不会被裁减掉。

(3)函数trimPopulation

        由于父子种群的合并，使得函数rankAndDistance中的popSize为两倍的种群大小(该种群大小在options中设置，稍后述及),而函数rankAndDistance中的nParents即为种群大小，故其中的条件“nParents==popSize”不成立，需要调用函数trimPopulation以修剪种群。种群修剪函数trimPopulation的作用是：在两倍于种群大小的个体中修剪出个数等于种群大小的个体

种群修剪函数trimPopulation的基本思想是：

        ①根据设定的系数ParetoFraction计算第一前端中允许保留的个体数目，按照一定的公式计算其余前端中允许保留的个体数目，则某前端中保留的个体数目为min{允许保留的个体数目，现存的个体数目},也就是说，对于第一前端，所设定的系数ParetoFraction直接决定了该前端中允许保留的个体数目，当允许保留的个体数目小于前端中现存的个体数目时，系数ParetoFraction所决定的允许保留的个体数目对该前端中保留的个体数目有限制作用，对于其他前端，也有类似思想。

        ②某前端中保留的个体数目计算出来以后，剩下的就是执行了，也就是说，将该前端中的个体数目修剪至保留的个体数目，这是通过锦标赛选择实现的。前已述及，对于同一前端，个体的拥挤距离越小，则多样性越差，因此，在锦标赛选择中，就越应该成为失败者而被选中淘汰，这种思想与函数selectiontournament中将expectation赋值为[-rank,Distance]是异曲同工的。另外，个体的去除通过将其赋值为空矩阵巧妙地实现。

        ③通过以上操作后，若保留下来的各前端的个体之和比Options中设置的种群大小多，那么需要进一步修剪，将最终的种群大小精减为Options中设置的种群大小。其中所用的思想与函数selectiontournament的选择思想也是一致的。

        5.函数distanceAndSpread

        该函数的作用是计算种群的平均距离和spread,后者参与图2中函数gamultiobjConverged对终止条件的判断。

3.3使用函数gamultiobj求解多目标优化问题

        (1)使用函数gamultiobj求解多目标优化问题的第一步是编写目标函数的M文件。对于以上问题，函数名为my_first_multi,目标函数代码如下：

function f = my_first_multi(x)
 
f(1) = x(1)^4 - 10*x(1)^2+x(1)*x(2) + x(2)^4 - (x(1)^2)*(x(2)^2);
f(2) = x(2)^4 - (x(1)^2)*(x(2)^2) + x(1)^4 + x(1)*x(2);

        (2)使用命令行方式调用gamultiobj函数，代码如下：

clear
clc
 
 
fitnessfcn = @my_first_multi;   % Function handle to the fitness function
nvars = 2;                      % Number of decision variables
lb = [-5,-5];                   % Lower bound
ub = [5,5];                     % Upper bound
A = []; b = [];                 % No linear inequality constraints
Aeq = []; beq = [];             % No linear equality constraints
options = gaoptimset('ParetoFraction',0.3,'PopulationSize',100,'Generations',200,'StallGenLimit',200,'TolFun',1e-100,'PlotFcns',@gaplotpareto);
 
[x,fval] = gamultiobj(fitnessfcn,nvars, A,b,Aeq,beq,lb,ub,options);

       其中，fitnessfcn即(1)中定义的目标函数M文件，设置的最优前端个体系数ParetoFraction为0.3,种群大小PopulationSize为100,最大进化代数Generations为200,停止代数StallGenLimit也为200,适应度函数值偏差TolFun为le-100,绘制Pareto前端。当然，也可以通过GUI方式调用函数gamultiobj,其方法与函数ga的调用相同，这里不再赘述。

3.4 结果分析

        可以看到，在基于遗传算法的多目标优化算法的运行过程中，自动绘制了第一前端中个体的分布情况，且分布随着算法进化一代而更新一次。当迭代停止后，得到如图5所示的第一前端个体分布图。同时，Worksapce中返回了函数gamultiobj得到的Pareto解集x及与x对应的目标函数值，如表1所列。需要说明的是，由于算法的初始种群是随机产生的，因此每次运行的结果不一样。

X1	X2	F1	F2
-0.702604978119867	0.702952585998864	-5.18650007950382	-0.249962526715632
-0.702604978119867	0.702952585998864	-5.18650007950382	-0.249962526715632
-2.67029628249862	1.97094984179847	-38.3329919666111	32.9718303966485
-2.63333082043443	1.96439479904977	-38.2990862211675	31.0452258773312
-2.08276662583749	1.31511944032696	-31.8120240265547	11.5671441504702
-1.15140918299647	0.929595004170188	-12.9690674486452	0.288363618240700
-2.24644109783116	1.74843483243066	-35.0074620102752	15.4575140499734
-2.41514714936595	1.79906399485147	-37.0545025942865	21.2748549366181
-2.38150536923245	1.86154289865626	-36.6275981206039	20.0880801162262
-1.82290570989828	1.23859537968056	-27.1897975855276	6.04005468627002
-1.30683826246950	1.00106276130680	-16.1770194512696	0.901242991273376
-2.51924521816231	1.97762256309168	-37.6944252011462	25.7715394911906
-1.52418626573825	1.14877588640926	-20.9096383290698	2.32179939758130
-1.06426303292729	0.940627721164028	-11.2640394037852	0.0625186287707547
-2.03691502287105	1.46637748311609	-31.5605834251883	9.92964467878934
-1.25902473471341	0.812532300017947	-14.9724017318298	0.879031094371823
-2.54579676153300	1.91445377012727	-38.0010331247212	26.8097783855977
-1.20396912943971	1.11637502706332	-13.9916338566617	0.503782789776543
-2.48983479640961	1.85584414778789	-37.6715458295670	24.3212273045537
-2.20516184528621	1.75854178206263	-34.3335326972675	14.2938549417933
-2.01646480190617	1.57180853250172	-31.2393168080679	9.42198616519718
-2.29454764988585	1.76104970025566	-35.6807426169950	16.9687465589719
-1.45000820746889	1.40337975638453	-18.9015806104910	2.12365740678036
-1.97425293081786	1.31617729318077	-30.1344339638333	8.84231238459566
-1.74540792332427	1.19107078666463	-25.5718188853058	4.89266930272549
-2.67029628249862	1.97094984179847	-38.3329919666111	32.9718303966485
-0.903690328921961	0.874168873899340	-8.32972021598921	-0.163158110118390
-1.00010614649929	0.757833629162209	-10.0042121434185	-0.00208910076199753
-2.16416470915156	1.47815692188379	-33.5583157841867	13.2777730991840
-2.57929181596638	1.94522918966795	-38.1411557075827	28.3863070115286

        从图5可以看到，第一前端的Pareto最优解分布均匀。从表1可以看到，返回的Pareto最优解个数为30个，而种群大小为100,可见，ParetoFraction为0.3的设置发挥了作用。另外，个体被限制在了[-5,5]的上下限范围内。