heixianghanshuyouhua.txt

黑箱优化算法

function [Opt_prmt,Opt_obj,Opt_obj2,LC]=RODGSA(M,N,Pm,K,t0,alpha,net1,net2,net3)
%输入参数列表
%  M            遗传进化迭代次数
%  N            种群规模
%  Pm           变异概率
%  K            同一温度下的状态跳转次数
%  t0           初始温度
%  alpha        降温系数
%
%输出参数列表
%  Opt_prmt     最优结构参数
%  Opt_obj      目标函数的最优值
%  LC           算法的收敛曲线

%第一步：变量初始化
load LHC%导入原始数据
%[net1,net2,net3]=BPNET(P1,T1,P2,T2,P3,T3)%构造并训练BP神经网络
[g_D,g_L,g_hm,g_d,g_px,g_n,g_pin,g_kl]=Gene_Construct;%构造基因库
Opt_prmt=zeros(M,8);
Opt_obj=zeros(1,M);
LC=zeros(1,M);

%第二步：产生初始种群
farm=zeros(N,8);
farm(1,:)=[120,75,0.06,0.6,0.3,8,0.4,0.6];
Len_g_D=length(g_D);
Len_g_L=length(g_L);
Len_g_hm=length(g_hm);
Len_g_d=length(g_d);
Len_g_px=length(g_px);
Len_g_n=length(g_n);
Len_g_pin=length(g_pin);
Len_g_kl=length(g_kl);
for i=2:N
    farm(i,1)=g_D(unidrnd(Len_g_D));
    farm(i,2)=g_L(unidrnd(Len_g_L));
    farm(i,3)=g_hm(unidrnd(Len_g_hm));
    farm(i,4)=g_d(unidrnd(Len_g_d));
    farm(i,5)=g_px(unidrnd(Len_g_px));
    farm(i,6)=g_n(unidrnd(Len_g_n));
    farm(i,7)=g_pin(unidrnd(Len_g_pin));
    farm(i,8)=g_kl(unidrnd(Len_g_kl));
end

counter=0;%设置迭代计数器
T=t0;
while counter<M%停止条件为达到最大迭代次数
       
    %第三步：双亲双子交叉
    newfarm=zeros(2*N,8);
    Ser=randperm(N);
    AA=farm(Ser(1),:);
    BB=farm(Ser(2),:);
    cp=unidrnd(7);
    aa=[AA(1:cp),BB((cp+1):end)];
    bb=[BB(1:cp),AA((cp+1):end)];
    newfarm(2*N-1,:)=aa;
    newfarm(2*N,:)=bb;
    for i=1:(N-1)
        AA=farm(Ser(i),:);
        BB=farm(Ser(i+1),:);
        cp=unidrnd(7);
        aa=[AA(1:cp),BB((cp+1):end)];
        bb=[BB(1:cp),AA((cp+1):end)];
        newfarm(2*i-1,:)=aa;
        newfarm(2*i,:)=bb;
    end
    FARM=[farm;newfarm];
   
    %第四步：选择复制
    Qer=randperm(3*N);
    FITNESS=zeros(3*N,1);
    fitness=zeros(N,1);
    for i=1:(3*N)
        Ind=FARM(i,:);
        FITNESS(i)=Fit(Ind,net1,net2,net3,P1,T1,P2,T2,P3,T3);
    end
    for i=1:N
        f1=FITNESS(Qer(3*i-2));
        f2=FITNESS(Qer(3*i-1));
        f3=FITNESS(Qer(3*i));
        if f1>=f2&&f1>=f3
            farm(i)=FARM(Qer(3*i-2));
            fitness(i)=FITNESS(Qer(3*i-2));
        elseif f2>=f1&&f2>=f3
            farm(i)=FARM(Qer(3*i-1));
            fitness(i)=FITNESS(Qer(3*i-1));
        elseif f3>=f1&&f3>=f2
            farm(i)=FARM(Qer(3*i));
            fitness(i)=FITNESS(Qer(3*i));
        else
        end
    end
    %记载最佳个体和平均适应值
    maxfit=max(fitness);
    Opt_obj(counter+1)=maxfit;
    pos=find(fitness==maxfit);
    Opt_prmt(counter+1,:)=farm(pos(1),:);
    LC(counter+1)=mean(fitness);
   
    %第五步：变异退火
    for i=1:N
        kk=1;
        if Pm>rand
            AAA=farm(i,:);
            while kk<K
                BBB=AAA;
                Bpos=unidrnd(8);
                switch Bpos
                    case 1,
                        BBB(1)=g_D(unidrnd(Len_g_D));
                    case 2,
                        BBB(2)=g_L(unidrnd(Len_g_L));
                    case 3,
                        BBB(3)=g_hm(unidrnd(Len_g_hm));
                    case 4,
                        BBB(4)=g_d(unidrnd(Len_g_d));
                    case 5,
                        BBB(5)=g_px(unidrnd(Len_g_px));
                    case 6,
                        BBB(6)=g_n(unidrnd(Len_g_n));
                    case 7,
                        BBB(7)=g_pin(unidrnd(Len_g_pin));
                    case 8,
                        BBB(8)=g_kl(unidrnd(Len_g_kl));
                    otherwise,
                end
                kk=kk+1;
                fitAAA=Fit(AAA,net1,net2,net3,P1,T1,P2,T2,P3,T3);
                fitBBB=Fit(BBB,net1,net2,net3,P1,T1,P2,T2,P3,T3);
                if fitBBB>=fitAAA
                    AAA=BBB;
                    farm(i,:)=AAA;
                elseif exp((fitBBB-fitAAA)/T)>rand
                    AAA=BBB;
                    farm(i,:)=AAA;
                else
                end
            end
        end
    end
   
    %第六步：进化退火
    counter=counter+1
    T=T*alpha;
end

%第七步：绘收敛曲线图
Opt_obj2=Opt_obj;
for i=1:M
    leijia=Opt_obj2(1:i);
    Opt_obj2(i)=max(leijia);
end
plot(Opt_obj2)
hold on
plot(LC)
Best_obj=max(Opt_obj)
Pos=find(Opt_obj==Best_obj);
Best_prmt=Opt_prmt(Pos(end),:)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 构造并训练BP神经网络
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [net1,net2,net3]=BPNET(P1,T1,P2,T2,P3,T3)
[PN1,minP1,maxP1]=premnmx(P1);
[TN1,minT1,maxT1]=premnmx(T1);
net1=newff(minmax(PN1),[6,20,3],{'tansig','tansig','purelin'},'trainlm');
net1.trainParam.show=5;
net1.trainParam.lr=0.1;
net1.trainParam.epochs=2000;
net1.trainParam.goal=0.000001;
[net1,tr1]=train(net1,PN1,TN1)
AN1=sim(net1,PN1)
[A1]=postmnmx(AN1,minT1,maxT1)
error1=(A1-T1)./T1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
P2=[P2(:,1:16),P2(:,18:end)];%剔除奇异样本点
T2=[T2(:,1:16),T2(:,18:end)];
[PN2,minP2,maxP2]=premnmx(P2);
[TN2,minT2,maxT2]=premnmx(T2);
net2=newff(minmax(PN2),[6,20,2],{'tansig','tansig','purelin'},'trainlm');
net2.trainParam.show=5;
net2.trainParam.lr=0.1;
net2.trainParam.epochs=2000;
net2.trainParam.goal=0.0000001;
[net2,tr2]=train(net2,PN2,TN2)
AN2=sim(net2,PN2)
[A2]=postmnmx(AN2,minT2,maxT2)
error2=(A2-T2)./T2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[PN3,minP3,maxP3]=premnmx(P3);
[TN3,minT3,maxT3]=premnmx(T3);
net3=newff(minmax(PN3),[4,10,1],{'tansig','tansig','purelin'},'trainlm');
net3.trainParam.show=5;
net3.trainParam.lr=0.1;
net3.trainParam.epochs=2000;
net3.trainParam.goal=0;
[net3,tr3]=train(net3,PN3,TN3)
AN3=sim(net3,PN3)
[A3]=postmnmx(AN3,minT3,maxT3)
error3=(A3-T3)