matlab遗传算法

function [xv,fv] = myGA(fitness, a, b, NP, NG, Pc, Pm, eps)
%       用遗传算法求解一维无约束优化问题
%
%       待优化的目标函数 fitness
%       自变量下界 a 
%       自变量上界 b 
%       种群个体数 NP
%       最大进化代数 NG
%       杂交概率 Pc
%       变异概率 Pm
%       自变量离散精度 eps
%       目标变量取最大值时自变量的值: xm
%       目标函数的最大值 fv
%
%       Example:
%               function F = fitness(x)
%               F = x^3-60*x^2+900*x+100;   
%   -------------------------------
%               [xv,fv] = myGA(@fitness,0, 30, 50, 100, 0.9, 0.04, 0.01);
%       --------------------------------------------------
%               xv = 10
%               fv = 4100
%
%       本程序在《精通MATLAB最优化计算》页315程序的基础上修改

L = ceil(log2((b-a) / eps + 1));                                %编码长度
x = zeros(NP, L);                                                               %种群
nx = zeros(size(x));                                                    %滚动数组
fx = zeros(NP, 1);                                                              %适应度
for i = 1:NP
        x(i,:) = Initial(L);
end

fv = -inf;

for k = 1 : NG
        for i = 1 : NP
                fx(i) = fitness(Dec(a, b, x(i, :), L));
                if (fx(i) > fv)
                        xv = Dec(a, b, x(i, :), L);
                        fv = fx(i);
                end
        end
        
        sumfx = sum(fx);
        Px = fx / sumfx;
        
        PPx = zeros(NP, 1);
        PPx(1) = Px(1);                                                                 %概率叠加
        for i = 2 : NP
                PPx(i) = PPx(i - 1) + Px(i);
        end
        
        selFather = 0;
        for i = 1 : NP
                sita = rand();
                for j = 1 : NP
                        if (sita <= PPx(j))
                                selFather = j;                                          %使用轮盘赌法进行选择父亲
                                break;
                        end
                end
                
                selMother = floor(rand() * NP) + 1;                     %母亲随机选择
                posCut = floor(rand() * (L - 1)) + 1;           %交叉点
                
                r1 = rand();
                if (r1 <= Pc)
                        nx(i, 1 : posCut) = x(selFather, 1:posCut);
                        nx(i, (posCut + 1) : L) = x(selMother, (posCut + 1) : L);
                        r2 = rand();
                        if (r2 <= Pm)
                                posMut = floor(rand() * L) + 1;
                                nx(i, posMut) = ~nx(i, posMut);
                        end
                else
                        nx(i, :) = x(selFather, :);
                end
        end
        
        x = nx;
end

%--------------------------------------------------------
%       初始化种群
function  result = Initial(length)
result = zeros(size(length()));
for i = 1 : length
        r = rand();
        result(i) = round(r);
end

%----------------------------------------------------------
%       编码转换
function y = Dec(a, b, x, L)
base = 2 .^ ((L - 1) : -1: 0);
y = dot(base, x);
y = a + y * (b - a) / (2 ^ L - 1);