x = linspace(0, 12.1, 200);
y = linspace(4.1, 5.8, 200);
[x, y] = meshgrid(x,y);
f = 21.5+x.*sin(4*pi*x)+y.*sin(20*pi*y);
figure;
surf(x,y,f);
colormap(gca, 'jet')
xlabel('x');
ylabel('y');
zlabel('f(x,y)');
title('函数图像')

function [Chrom,ObjV]=immigrant(Chrom,ObjV)
%% 移民算子
MP=length(Chrom);
for i=1:MP
    [MaxO,maxI]=max(ObjV{i});  % 找出第i种群中最优的个体
    next_i=i+1;                % 目标种群（移民操作中）
    if next_i>MP;next_i=mod(next_i,MP);end
    [MinO,minI]=min(ObjV{next_i});          %  找出目标种群中最劣的个体
    %% 目标种群最劣个体替换为源种群最优个体
    Chrom{next_i}(minI,:)=Chrom{i}(maxI,:);
    ObjV{next_i}(minI)=ObjV{i}(maxI);
end

function [MaxObjV,MaxChrom]=EliteInduvidual(Chrom,ObjV,MaxObjV,MaxChrom)
%% 人工选择算子
MP=length(Chrom);  %种群数
for i=1:MP
    [MaxO,maxI]=max(ObjV{i});   %找出第i种群中最优个体
    if MaxO>MaxObjV(i)
        MaxObjV(i)=MaxO;         %记录各种群的精华个体
        MaxChrom(i,:)=Chrom{i}(maxI,:);  %记录各种群精华个体的编码
    end
end

function obj=ObjectFunction(X)
%% 待优化的目标函数
col=size(X,1);
for i=1:col
    obj(i,1)=21.5+X(i,1)*sin(4*pi*X(i,1))+X(i,2)*sin(20*pi*X(i,2));
%     obj(i,1)=exp(((X(i,1)-0.1)/0.8)^2)*(sin(5*pi*X(i,1)))^6;
end

%% 标准遗传算法SGA
clear;
%clc
pc=0.7;     % 交叉概率
pm=0.05; % 变异概率
%定义遗传算法参数
NIND=40;        %个体数目
MAXGEN=500;     %最大遗传代数
NVAR=2;               %变量的维数
PRECI=20;             %变量的二进制位数
GGAP=0.9;             %代沟
trace=zeros(MAXGEN,1); %记录优化轨迹
FieldD=[rep(PRECI,[1,NVAR]);[-3,4.1;12.1,5.8];rep([1;0;1;1],[1,NVAR])]; %建立区域描述器
Chrom=crtbp(NIND, NVAR*PRECI);                       %创建初始种群
gen=0;                                               %代计数器   
ObjV=ObjectFunction(bs2rv(Chrom, FieldD));%计算初始种群个体的目标函数值
[maxY,I]=max(ObjV); %最优值
X=bs2rv(Chrom, FieldD);
maxX=X(I,:);
while gen<MAXGEN                                     %迭代
    FitnV=ranking(-ObjV);                            %分配适应度值(Assign fitness values)
    SelCh=select('sus', Chrom, FitnV, GGAP);         %选择
    SelCh=recombin('xovsp', SelCh, pc);              %重组
    SelCh=mut(SelCh,pm);                             %变异
    ObjVSel=ObjectFunction(bs2rv(SelCh, FieldD));           %计算子代目标函数值
    [Chrom ObjV]=reins(Chrom, SelCh, 1, 1, ObjV, ObjVSel);  %重插入
    gen=gen+1;           %代计数器增加
    if maxY<max(ObjV)
        [maxY,I]=max(ObjV);
        X=bs2rv(Chrom, FieldD);
        maxX=X(I,:);
    end
    trace(gen,1)=maxY;
end
%% 进化过程图
plot(1:gen,trace(:,1));
hold on
xlabel('进化代数');
ylabel('最优解变化');
title('SGA进化过程');
%% 输出最优解
disp(['最优值为：',num2str(maxY)]);
disp(['对应的自变量取值：',num2str(maxX)]);

%%
%% 2、多种群遗传算法
clear;
clc
NIND=40;               %个体数目
NVAR=2;                %变量的维数
PRECI=20;              %变量的二进制位数
GGAP=0.9;              %代沟
MP=10;                 %种群数目
FieldD=[rep(PRECI,[1,NVAR]);[-3,4.1;12.1,5.8];rep([1;0;1;1],[1,NVAR])];  %译码矩阵
for i=1:MP
    Chrom{i}=crtbp(NIND, NVAR*PRECI);                       %创建初始种群
end
pc=0.7+(0.9-0.7)*rand(MP,1);    %在【0.7,0.9】范围i内随机产生交叉概率
pm=0.001+(0.05-0.001)*rand(MP,1);  %在【0.001,0.05】范围内随机产生变异概率
gen=0;  %初始遗传代数
gen0=0; %初始保持代数
MAXGEN=10;  %最优个体最少保持代数
maxY=0; %最优值
for i=1:MP
    ObjV{i}=ObjectFunction(bs2rv(Chrom{i}, FieldD));%计算各初始种群个体的目标函数值
end
MaxObjV=zeros(MP,1);           %记录精华种群
MaxChrom=zeros(MP,PRECI*NVAR); %记录精华种群的编码
while gen0<=MAXGEN
    gen=gen+1;       %遗传代数加1
    for i=1:MP
        FitnV{i}=ranking(-ObjV{i});                      % 各种群的适应度
        SelCh{i}=select('sus', Chrom{i}, FitnV{i},GGAP); % 选择操作
        SelCh{i}=recombin('xovsp',SelCh{i}, pc(i));      % 交叉操作
        SelCh{i}=mut(SelCh{i},pm(i));                    % 变异操作
        ObjVSel=ObjectFunction(bs2rv(SelCh{i}, FieldD)); % 计算子代目标函数值
        [Chrom{i},ObjV{i}]=reins(Chrom{i},SelCh{i},1,1,ObjV{i},ObjVSel);    %重插入操作
    end
    [Chrom,ObjV]=immigrant(Chrom,ObjV);     % 移民操作
    [MaxObjV,MaxChrom]=EliteInduvidual(Chrom,ObjV,MaxObjV,MaxChrom);     % 人工选择精华种群
    YY(gen)=max(MaxObjV);    %找出精华种群中最优的个体
    if YY(gen)>maxY   %判断当前优化值是否与前一次优化值相同
        maxY=YY(gen); %更新最优值
        gen0=0;
    else
        gen0=gen0+1; %最优值保持次数加1
    end
end
%% 进化过程图
plot(1:gen,YY);
xlabel('进化代数')
ylabel('最优解变化')
title('MPGA进化过程')
%% 输出最优解
[Y,I]=max(MaxObjV);    %找出精华种群中最优的个体
X=(bs2rv(MaxChrom(I,:), FieldD));   %最优个体的解码解
disp(['最优值为：',num2str(Y)])
disp(['对应的自变量取值：',num2str(X)])