一周前，我的遗传算法在数值优化上再次失败。在此之前我尝试过遗传算法结合局部搜索(随机爬山、SQP、模式搜索)，由于INTJ型人格认为结合局部信息的进化不是好进化，转而尝试了其他的方法：小生境技术(niching)，以适应度共享的手段来维持种群的多样性——在某些函数上取得了极好的效果，而参数的自适应控制却非常困难；多种群技术(MPGA)——由于我的迁徙算法错误使得我的程序本质上就是用不同的控制参数进行多次计算；正交技术(OGA)，通过构造正交矩阵，使得在初始化种群和产生子代时能够均匀搜索整个解空间——在几个测试函数上都取得了不错的成果，但是后来我意识到这明明就是带启发式策略的网格搜索(grid

search)，以及在一次离散事件仿真模拟的优化上得到了很糟糕的结果，加上计算开销大、运算速度慢，很快被我打入冷宫。再加上遗传算法在集合覆盖问题、TSP、Packing问题上令我蛋疼的表现，一度我对遗传算法的态度是“鸡肋，还不如用分支定界”，转而研究传统优化算法——序列二次规划(sequential

quadratic programming,SQP)，在组合优化上也开始关注遗传算法结合回溯搜索。

由于在五月份的对仿真模拟的优化问题上，我使用的NPMEA算法在处理带约束的问题上让我印象深刻(之前我只用著名的G8函数做过测试)，上周我怀着极大的信心采用它的单纯形交叉算子和自适应的高斯变异来做无约束优化(或者说是bound

constraint)，但令人失望的是在某个五元函数上几乎连全局最优点的附近都搜索不到，增加了多种群和小生境运算速度又太慢，实在令人垂头丧气。

最后我抱着死马当活马医的态度打开gatool——我对这东西从来不屑一顾的，花里胡哨，不如自己coding来得快——结果太可怕了，一切参数采用默认设置，Schaffer函数以1的概率收敛到全局最优，颠覆三观！！我当时又是失落又是激动，心里只有两个字：“GADST，你为何这么叼？”一怒之下，把GADST的ga.m拆了，恍然大悟。

二、拆ga.m

step1：ga.m设置默认defaultopt

%-----------------------------e.g.---------------------------------

%代码取自ga.m 181th~202th rows

%默认编码类型为双精度编码

%默认生成初始种群向量下界为0上界为1

%默认种群大小为20

%默认使用gacreationuniform函数来生成均匀随机分布的初始种群

defaultopt = struct('PopulationType', 'doubleVector', ...

'PopInitRange', [0;1],

...

'PopulationSize', 20,

...

'CreationFcn',@gacreationuniform);

step2：ga.m作各种检查和判断

%-----------------------------e.g.---------------------------------

%代码取自ga.m 251th~277th rows

%检查第十个输入变量，若存在且是结构体，则为整数约束规划

if nargin == 10 &&

isstruct(intcon)

options = intcon;

intcon = [];

end

%检查输入的FitnessFcn是不是function handles或者inlines，若不是则返回error

if isempty(FitnessFcn) ||

~(isa(FitnessFcn,'inline') ||

isa(FitnessFcn,'function_handle'))

error(message('globaloptim:ga:needFunctionHandle'));

end

step3：ga.m对空的options结构域使用默认defaultopt

%-----------------------------e.g.---------------------------------

%代码取自ga.m 304th~308th rows

%若options为空，则使用默认的defaultopt

if ~isempty(options) &&

~isa(options,'struct')

error(message('globaloptim:ga:optionsNotAStruct'));

elseif isempty(options)

options =

defaultopt;

end

step4：ga.m调用子函数gacommon判断问题类型和种群编码类型

%-----------------------------e.g.---------------------------------

%代码取自gacommon.m

14th~33th rows

%若结构体nonlcon非空则为非线性约束规划，若线性等式约束矩阵非空则为线性约束规划，若上下界矩阵非空

%则为边界约束规划，若都为空则为无约束规划

if ~isempty(nonlcon)

type =

'nonlinearconstr';

elseif ~isempty(Aeq) || ~isempty(beq) || ~isempty(Aineq) ||

~isempty(bineq)

type =

'linearconstraints';

elseif ~isempty(lb) || ~isempty(ub)

type =

'boundconstraints';

else

type =

'unconstrained';

end

step5：ga.m对不同类型的问题调用不同子求解器

%-----------------------------e.g.---------------------------------

%代码取自ga.m 349th~359th rows

%'unconstrained'调用无约束规划求解器gaunc,'boundconstraints'和'linearconstraints'调用线性规划求解

%器galincon,'nonlinearconstr'调用非线性规划求解器gacon。

switch (output.problemtype)

case

'unconstrained'

[x,fval,exitFlag,output,population,scores] =

gaunc(FitnessFcn,nvars, ...

options,output,Iterate);

case

{'boundconstraints', 'linearconstraints'}

[x,fval,exitFlag,output,population,scores] =

galincon(FitnessFcn,nvars, ...

Aineq,bineq,Aeq,beq,lb,ub,options,output,Iterate);

case

'nonlinearconstr'

[x,fval,exitFlag,output,population,scores] =

gacon(FitnessFcn,nvars, ...

Aineq,bineq,Aeq,beq,lb,ub,NonconFcn,options,output,Iterate,type);

end

(下面以boundconstraints为例,调用galincon)

step6：在galincon.m中调用子函数makeState来初始化结构数组state

%-----------------------------e.g.---------------------------------

%代码取自makeState.m 10th~30th rows

%初始化当前进化代数和当前停止代数

state.Generation = 0;

state.LastImprovement = 1;

%调用options.CreationFcn里的函数来生成初始种群，见step7

state.Population =

feval_r(options.CreationFcn,GenomeLength,FitnessFcn,options,

options.CreationFcnArgs{:});

(下面以默认的gacreationuniform函数为例，生成初始种群)

step7：在gacreationuniform.m里通过用户提供的初始种群和初始种群范围来随机生产种群

%-----------------------------e.g.---------------------------------

%代码取自gacreationuniform.m 42th~46th rows

%在PopInitRange范围内生成均匀分布的随机初始种群

%等价于Population=

unifrnd(lb,ub,initPopProvided+1,PopSize);

range = options.PopInitRange;

lowerBound = range(1,:);

span = range(2,:) - lowerBound;

Population(initPopProvided+1:end,:) =

repmat(lowerBound,individualsToCreate,1) + ...

repmat(span,individualsToCreate,1) .*

rand(individualsToCreate,GenomeLength);

step8：在galincon.m中调用子函数gadsplot继续初始化结构数组state，为作图做准备

略..

step9：在galincon.m中判断是否满足进化结束条件，若不满足则进化一次

%-----------------------------e.g.---------------------------------

%代码取自galincon.m 48th~91th rows

while isempty(exitFlag)

%进化代数加一

state.Generation =

state.Generation + 1;

%函数stepGA内完成排序\选择\交叉\变异和种群的进化，见step10

[score,population,state]

=

stepGA(score,population,options,state,GenomeLength,FitnessFcn);

%记录最佳个体

best =

min(state.Score);

generation =

state.Generation;

state.Best(generation) =

best;

%判断和进行迁徙操作(仅针对多种群)

state =

migrate(FitnessFcn,GenomeLength,options,state);

%更新图像输出

state =

gadsplot(options,state,currentState,'Genetic Algorithm');

end

step10：stepGA.m产生精英后代、交叉后代和变异后代

%-----------------------------e.g.---------------------------------

%代码取自stepGA.m 9th~36th rows

%精英后代数目

nEliteKids = options.EliteCount;

%交叉后代数目

nXoverKids = round(options.CrossoverFraction *

(size(thisPopulation,1) - nEliteKids));

%变异后代数目

nMutateKids = size(thisPopulation,1) - nEliteKids -

nXoverKids;

%用于产生交叉和变异所需的父代数目

nParents = 2 * nXoverKids + nMutateKids;

%适应度排序操作，见step11

state.Expectation =

feval_r(options.FitnessScalingFcn,thisScore,nParents,...

options.FitnessScalingFcnArgs{:});

%选择交叉和变异所需的后代，见step12

parents =

feval_r(options.SelectionFcn,state.Expectation,nParents,options,...

options.SelectionFcnArgs{:});

[unused,k] = sort(thisScore);

%产生精英后代

eliteKids  =

thisPopulation(k(1:options.EliteCount),:);

%产生交叉后代，见step13

xoverKids  = feval_r(options.CrossoverFcn,

parents(1:(2 * nXoverKids)),options,GenomeLength,...

FitnessFcn,thisScore,thisPopulation,options.CrossoverFcnArgs{:});

%产生交叉后代，见step14

mutateKids = feval_r(options.MutationFcn,

parents((1 + 2 * nXoverKids):end),

options,GenomeLength,

FitnessFcn,state,thisScore,thisPopulation,options.MutationFcnArgs{:});

step11：在stepGA.m中调用options.FitnessScalingFcn中的适应度排序函数进行适应度排序操作

这里用默认的fitscalingrank函数

%-----------------------------e.g.---------------------------------

%代码取自fitscalingrank.m 16th~21th rows

%expectation将score映射到从1开始的整数的平方根的倒数，并作比例处理

[~,i] = sort(scores);

expectation = zeros(size(scores));

expectation(i) = 1 ./ ((1:length(scores))  .^

0.5);

expectation = nParents * expectation ./

sum(expectation);

step12：在stepGA.m中调用options.SelectionFcn中的选择函数，根据各父代的适应度选出参与交叉变异的父代。

这里用默认的selectionstochunif函数

%-----------------------------e.g.---------------------------------

%代码取自selectionstochunif.m 22th~44th rows

%对expectation累计求和及对[0,1]等分为stepSize后，position从第一等分的任意位置开始，找到第一个大于

%position的wheel的索引作为选出的一个子代，position前进一个stepSize

wheel = cumsum(expectation) / nParents;

stepSize = 1/nParents;

position = rand * stepSize;

lowest = 1;

for i = 1:nParents % for each parent needed,

for j =

lowest:length(wheel)

if(position <

wheel(j))

parents(i)

= j;

lowest =

j;

break;

end

end

position = position +

stepSize;

end

step13：在stepGA.m中调用options.CrossoverFcn中的交叉函数，产生交叉后代。

这里用默认的selectionstochunif函数

%-----------------------------e.g.---------------------------------

%代码取自crossoverscattered.m 42th~48th rows

%对子代的每个位置，取[0,1]随机数，若随机数大于0.5则取父代A对应位置的数，反之取父代B

for j = 1:GenomeLength

if(rand >

0.5)

xoverKids(i,j) = thisPopulation(r1,j);

else

xoverKids(i,j) = thisPopulation(r2,j);

end

end

step14：在stepGA.m中调用options.MutationFcn中的变异函数，让选出的父代变异产生后代。

MATLAB默认的mutationadaptfeasible函数采用了复杂的算法来确定变异的搜索方向(search

directions)

%-----------------------------e.g.---------------------------------

%代码简化自mutationadaptfeasible.m 58th~126th rows

%MeshSize是自适应网格步长，和模式搜索类似，当遗传算法进化一代搜索到更优解时步长扩大为四倍，反之缩

%小为四分之一

%注意到MeshSize = 2^(-4^k)，可知pollParam一定为整数

pollParam = 1/sqrt(MeshSize);

%等价于lowerT = tril(unidrnd(pollParam, nGenome, nGenome));

lowerT = tril((round((pollParam+1)*rand(nGenome) -

0.5)),-1);

%等价于diagtemp = unidrnd(2*pollParam, nGenome)

- pollParam;

diagtemp = pollParam*sign(rand(nGenome,1) - 0.5);

diagT  = diag(diagtemp);

Basis = lowerT + diagT;

order = randperm(nGenome);

%方阵Basis每一列决定一个搜索方向向量

Basis = Basis(order,order);

%随机取搜索方向向量和方向符号来生成后代

indexVec = [1:nGenome 1:nGenome];

dirSign = [ones(1,nGenome) -ones(1,nGenome)];

OrderVec = randperm(2*nGenome);

for jj = 1:2*nGenome

direction =

dirSign(jj).*Basis(:,indexVec(OrderVec(jj)));

mutant = Parent(ii,:) +

MeshSize*direction';

%若生成后代满足约束则不再搜索其他方向

if all(mutant

>= lb & mutant <=

ub)

Kid(ii,:) = mutant;

break;

else

Kid(ii,:) = Parent(ii,:);

end

end