用免疫克隆算法求二元函数的最优值

求下面函数的最小值（最优解）

​ Goldstein-Price函数：

​ f₁(x)=[1+(x₁+x₂+1)²(19-14x₁+3x₁²-14x₂+6x₁x₂+3x₂²)]*[30+(2x₁-3x₂)²(18-32x₁+12x₁²+48x₂-36x₁x₂+27x₂²)]

​ -2≤x₁,x₂≤2

​ 要求：（1）用免疫克隆算法编程求解问题

​ 编程语言用MATLAB 或C 或Java

​ 输出问题的最优解、并绘图显示

​ 代码：bestvalue.m

clear all;          
clc;       

figure
[x,y]=meshgrid(-2:0.01:2,-2:0.01:2);
z=(1+((x+y+1).^2).*(19-14.*x+3.*x.^2-14.*y+6.*x.*y+3.*y.^2)).*(30+(2.*x-3.*y).^2.*(18-32.*x+12.*x.^2+48.*y-36.*x.*y+27.*y.^2));
mesh(x,y,z)
xlabel('x')
ylabel('y')
hold on
%初始化
         
D=2;                                  %免疫个体维数
NP=50;                                %免疫个体数目
Xs=4;                                 %上限
Xx=-4;                                %下限
G=50;                                %最大免疫代数
pm=0.7;                               %变异概率
alfa=2;                               %激励度参数
belta=1;                              %激励度参数
detas=0.2;                            %相似度阈值
gen=0;                                %免疫代数
Nc1=5;                                %克隆个数
deta0=0.5*Xs;                         %邻域范围初值
%初始种群
f=rand(D,NP)*(Xs-Xx)+Xx;
for np=1:NP
    MSLL(np)=fun(f(:,np));
end
%计算个体浓度和激励度
for np=1:NP
    for j=1:NP
        nd(j)=sum(sqrt((f(:,np)-f(:,j)).^2));
        if nd(j)<detas
            nd(j)=1;
        else
            nd(j)=0;
        end
    end
    ND(np)=sum(nd)/NP;
end
MSLL=alfa*MSLL-belta*ND;
%激励度按升序排列
[SortMSLL,Index]=sort(MSLL);
Sortf=f(:,Index);
%免疫循环
while gen<G
    for i=1:NP/2
        a=Sortf(:,i);
        Na=repmat(a,1,Nc1);
        deta=deta0/gen;
        for j=1:Nc1
            for ii=1:D
                if rand<pm
                    Na(ii,j)=Na(ii,j)+(rand-0.5)*deta;
                end
                %边界条件处理
                if (Na(ii,j)>Xs)||(Na(ii,j)<Xx)
                    Na(ii,j)=rand*(Xs-Xx)+Xx;
                end
            end
        end
        Na(:,1)=Sortf(:,i);%保留克隆源个体
        %克隆抑制，保留亲和度最高的个体
        for j=1:Nc1
            NaMSLL=fun(Na(:,j));
        end
        [NaSortMSLL,Index]=sort(NaMSLL);
        aMSLL(i)=NaSortMSLL(1);
        NaSortf=Na(:,Index);
        af(:,i)=NaSortf(:,1);
    end
    %免疫种群激励度
    for np=1:NP/2
        for j=1:NP/2
            nda(j)=sum(sqrt((af(:,np)-af(:,j)).^2));
            if nda(j)<detas
                nda(j)=1;
            else
                nda(j)=0;
            end
        end
        aND(np)=sum(nda)/NP/2;
    end
    aMSLL=alfa*aMSLL-belta*aND;
    %总群刷新
    bf=rand(D,NP/2)*(Xs-Xx)+Xx;
    for np=1:NP/2
        bMSLL(np)=fun(bf(:,np));
    end
    %生成新的种群激励度
    for np=1:NP/2
        for j=1:NP/2
            ndc(j)=sum(sqrt((bf(:,np)-bf(:,j)).^2));
            if ndc(j)<detas
                ndc(j)=1;
            else 
                ndc(j)=0;
            end
        end
        bND=sum(ndc)/NP/2;
    end
    bMSLL=alfa*bMSLL-belta*bND;
    %免疫种群与新生种群合并
    f1=[af,bf];
    MSLL=[aMSLL,bMSLL];
    [SortfMSLL,Index]=sort(MSLL);
    Sortf=f1(:,Index);
    gen=gen+1;
    trace(gen)=fun(Sortf(:,1));
end


%输出优化结果
Bestf=Sortf(:,1);
trace(end);
disp('最优值为：')
disp(trace(end))
disp('最优个体为：')
disp(Bestf)
plot3(Bestf(1), Bestf(2), trace(end),'bo','linewidth',1.5)
figure
plot(trace)
xlabel('迭代次数')
ylabel('目标函数值')
title('亲和度进化曲线')
%亲和度函数

%目标函数
function v=fun(x)
v=((1+(x(1)+x(2)+1).^2).*(19-14.*x(1)+3.*x(1).^2-14.*x(2)+6.*x(1).*x(2)+3.*x(2).^2)).*(30+(2.*x(1)-3.*x(2)).^2.*(18-32.*x(1)+12.*x(1).^2+48.*x(2)-36.*x(1).*x(2)+27.*x(2).^2));
end

结果：