ga.m

用遗传算法求解数学模型中参数的应用实例

%Generic Algorithm for function f(x) optimum
clear all;
close all;

K1=1;K2=197;mdn=97;    %设定炉温、钢温数据数组大小
%装载数据并存入数组中%
TE101=load('d:\GAprogram&data\data\data1\PDE#HD#HistStoveTemp#TE_101.txt');
F1=TE101(K1:K2);
TE102=load('d:\GAprogram&data\data\data1\PDE#HD#HistStoveTemp#TE_102.txt');
F2=TE102(K1:K2);
TE103=load('d:\GAprogram&data\data\data1\PDE#HD#HistStoveTemp#TE_103.txt');F3=TE103(K1:K2);
TE104=load('d:\GAprogram&data\data\data1\PDE#HD#HistStoveTemp#TE_104.txt');F4=TE104(K1:K2);
TE105=load('d:\GAprogram&data\data\data1\PDE#HD#HistStoveTemp#TE_105.txt');F5=TE105(K1:K2);
TE106=load('d:\GAprogram&data\data\data1\PDE#HD#HistStoveTemp#TE_106.txt');F6=TE106(K1:K2);
TE107=load('d:\GAprogram&data\data\data1\PDE#HD#HistStoveTemp#TE_107.txt');F7=TE107(K1:K2);
TE108=load('d:\GAprogram&data\data\data1\PDE#HD#HistStoveTemp#TE_108.txt');F8=TE108(K1:K2);
TE109=load('d:\GAprogram&data\data\data1\PDE#HD#HistStoveTemp#TE_109.txt');F9=TE109(K1:K2);
TE198=load('d:\GAprogram&data\data\data1\PDE#HD#HistStoveTemp#TE_198.txt');F98=TE198(K1:K2);
%处理钢温数据滤除不佳点%
for i=1:1:K2,
    if i<=(K2-1) 
        if F98(i)-F98(i+1)>60
            F98(i+1)=F98(i);
        end
        if F98(i+1)-F98(i)>60
            F98(i)=F98(i+1);
        end
    end
    Tjun(i)=(F1(i)+F2(i))/2;
    Tjia(i)=(F4(i)+F5(i))/2;
    TY(i)=(F7(i)+F8(i))/2;
    
    if(i>=2)
    F98(i)=F98(i)*0.3+F98(i-1)*0.7;
    end
    F98(i)=1090+(F98(i)-1090)*0.85;
end;
%入炉钢坯各点温度%
for  m=1:1:mdn,
    T1(1,m)=450;
    T2(1,m)=450;
    T3(1,m)=450;
    T4(1,m)=450;
    T5(1,m)=450;
    T6(1,m)=450;
end;
%遗传算法寻优找到最佳的总括热吸收率
%具体过程参考书目：玄光男,程润伟.遗传算法与工程优化.北京：清华大学出版社，2004.
%Parameters
Size=80;   %种群大小
G=50;      %最大繁衍代数     
CodeL=10;  %二进制编码长度 
 
umax=1.0;  %总括热吸收率最大值 
umin=0;    %总括热吸收率最小值

E=round(rand(Size,CodeL));    %Initial Code

%Main Program
for gek=1:1:G
time(gek)=gek;

for s=1:1:Size
    m=E(s,:);
    y1=0;

    %Uncoding
    m1=m(1:1:CodeL);
    for i=1:1:CodeL
        y1=y1+m1(i)*2^(i-1);
    end
    w=(umax-umin)*y1/1023+umin;%总括热吸收率
    %钢温模型
    for  j=1:1:mdn,   
        for k=2:1:34,
            Qb=w*(5.67E-8)*((TY(k+j-2)+273)^4-(T1((k-1),j)+273)^4);
            Qd=w*(5.67E-8)*((TY(k+j-2)+273)^4-(T6((k-1),j)+273)^4);
        
            T1(k,j)=0.814815*T1((k-1),j)+0.185185*T2((k-1),j)+0.000193*Qb;
            T2(k,j)=0.185185*T1((k-1),j)+0.629630*T2((k-1),j)+0.185185*T3((k-1),j);    
            T3(k,j)=0.185185*T2((k-1),j)+0.629630*T3((k-1),j)+0.185185*T4((k-1),j);
            T4(k,j)=0.185185*T3((k-1),j)+0.629630*T4((k-1),j)+0.185185*T5((k-1),j);
            T5(k,j)=0.185185*T4((k-1),j)+0.629630*T5((k-1),j)+0.185185*T6((k-1),j);
            T6(k,j)=0.814815*T6((k-1),j)+0.185185*T5((k-1),j)+0.000193*Qd;
        end;   
        for k=35:1:67,
            Qb=w*(5.67E-8)*((Tjia(k+j-2)+273)^4-(T1((k-1),j)+273)^4);
            Qd=w*(5.67E-8)*((Tjia(k+j-2)+273)^4-(T6((k-1),j)+273)^4);
        
            T1(k,j)=0.814815*T1((k-1),j)+0.185185*T2((k-1),j)+0.000193*Qb;
            T2(k,j)=0.185185*T1((k-1),j)+0.629630*T2((k-1),j)+0.185185*T3((k-1),j);    
            T3(k,j)=0.185185*T2((k-1),j)+0.629630*T3((k-1),j)+0.185185*T4((k-1),j);
            T4(k,j)=0.185185*T3((k-1),j)+0.629630*T4((k-1),j)+0.185185*T5((k-1),j);
            T5(k,j)=0.185185*T4((k-1),j)+0.629630*T5((k-1),j)+0.185185*T6((k-1),j);
            T6(k,j)=0.814815*T6((k-1),j)+0.185185*T5((k-1),j)+0.000193*Qd;
        end;       
        for k=68:1:102,
            Qb=w*(5.67E-8)*((Tjun(k+j-2)+273)^4-(T1((k-1),j)+273)^4);
            Qd=w*(5.67E-8)*((Tjun(k+j-2)+273)^4-(T6((k-1),j)+273)^4);
        
            T1(k,j)=0.814815*T1((k-1),j)+0.185185*T2((k-1),j)+0.000193*Qb;
            T2(k,j)=0.185185*T1((k-1),j)+0.629630*T2((k-1),j)+0.185185*T3((k-1),j);    
            T3(k,j)=0.185185*T2((k-1),j)+0.629630*T3((k-1),j)+0.185185*T4((k-1),j);
            T4(k,j)=0.185185*T3((k-1),j)+0.629630*T4((k-1),j)+0.185185*T5((k-1),j);
            T5(k,j)=0.185185*T4((k-1),j)+0.629630*T5((k-1),j)+0.185185*T6((k-1),j);
            T6(k,j)=0.814815*T6((k-1),j)+0.185185*T5((k-1),j)+0.000193*Qd;
        end;
    end;
    %计算目标函数SUM
    SUM=0;
    for i=1:1:mdn,
        A(i,1)=T1(102,i); A(i,2)=F98(100+i);
    end;
    for i=1:1:mdn,
        B(i,1)=A(i,1)-A(i,2);
    end;
    for i=1:1:mdn,
         SUM=SUM+B(i,1)^2;
    end;       
  F(s)=1.0/SUM;         
end

Ji=1./F;
%****** Step 1 : Evaluate BestJ ******
BestJ(gek)=min(Ji);
%适应度函数
fi=F;                          %Fitness Function
[Oderfi,Indexfi]=sort(fi);     %Arranging fi small to bigger
Bestfi=Oderfi(Size);           %Let Bestfi=max(fi)
BestS=E(Indexfi(Size),:);      %Let BestS=E(m), m is the Indexfi belong to max(fi)
bfi(gek)=Bestfi;

%****** Step 2 : Select and Reproduct Operation******
   fi_sum=sum(fi);
   fi_Size=(Oderfi/fi_sum)*Size;
   
   fi_S=floor(fi_Size);        %Selecting Bigger fi value
   
   kk=1;
   for i=1:1:Size
      for j=1:1:fi_S(i)        %Select and Reproduce 
         TempE(kk,:)=E(Indexfi(i),:);  
         kk=kk+1;              %kk is used to reproduce
      end
   end
   
%************ Step 3 : Crossover Operation ************
pc=0.60;
n=ceil(10*rand);
for i=1:2:(Size-1)
    temp=rand;
    if pc>temp                  %Crossover Condition
    for j=n:1:10
        TempE(i,j)=E(i+1,j);
        TempE(i+1,j)=E(i,j);
    end
    end
end
TempE(Size,:)=BestS;
E=TempE;
   
%************ Step 4: Mutation Operation **************
%pm=0.001;
%pm=0.001-[1:1:Size]*(0.001)/Size; %Bigger fi, smaller Pm
%pm=0.0;    %No mutation
pm=0.1;     %Big mutation

   for i=1:1:Size
      for j=1:1:CodeL
         temp=rand;
         if pm>temp                %Mutation Condition
            if TempE(i,j)==0
               TempE(i,j)=1;
            else
               TempE(i,j)=0;
            end
        end
      end
   end
   
%Guarantee TempPop(30,:) is the code belong to the best individual(max(fi))
TempE(Size,:)=BestS;
E=TempE;
end
 
Max_Value=Bestfi;
Min_Value=1/Bestfi;
BestS;
w   %显示最佳总括热吸收率

figure(1);
plot(time,BestJ); 
xlabel('Times');ylabel('Best J');
figure(2);
plot(time,bfi);
xlabel('times');ylabel('Best F');