coga_v.m

遗传算法例程：通过MATLAB

violthres = 50;     % Threshold of Violation
%%%%%%%%%%%%%%%%%%%%%%% Generating Test Data
FS = [];

F = 100;
T = 40;

% 只考虑单一机场的容量
FSu = [1,2,2,1,2,2,2,1,2,1,2,2,1,2,2,2,1,2,1,2,2,2,2,2,1,1,2,1,2,2,1,1,1,1,2,1,2,2,2,1,1,1,1,2,1,2,2,2,2,2,2,1,2,2,1,1,2,2,2,2,2,2,1,1,1,1,1,2,1,1,2,1,2,2,2,1,2,1,2,1,2,1,2,2,2,1,1,2,2,1,2,1,2,1,2,1,2,1,2,2;];
FSd = [2,1,1,2,1,1,1,2,1,2,1,1,2,1,1,1,2,1,2,1,1,1,1,1,2,2,1,2,1,1,2,2,2,2,1,2,1,1,1,2,2,2,2,1,2,1,1,1,1,1,1,2,1,1,2,2,1,1,1,1,1,1,2,2,2,2,2,1,2,2,1,2,1,1,1,2,1,2,1,2,1,2,1,1,1,2,2,1,1,2,1,2,1,2,1,2,1,2,1,1;];
% 起降机场表，FSu表示起飞机场序列，FSu中的1表示考虑容量限制的机场，其他不考虑者一律用2表示；
% FSd表示降落的机场序列，和上面类似；

FSt = [1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,1,2,1,2,2,1,1,1,1,1,1,1,1,1,1,1,2,1,2,1,1,1,2,2,1,1,1,2,1,1,1,1,2,1,1,1,2,1,1,2,1,2,1,2,1,1,1,2,1,1,1,2,1,1,1,1,1,1,1,1,1,1,2,1,1,2,1,2,1,1,1,2,1,1,2,1,1,1,1,2,1,1,1,1,1;];
% 飞机的类型序列，1表示中型飞机，2表示大型飞机

% Generate zs according to TdfsM & TrfsM,  

Cgf = 7;
Caf = 10;

Tpdfs = [11,10,11,12,18,4,2,14,27,29,8,5,27,8,20,30,20,27,1,5,25,13,27,23,21,11,5,5,6,13,26,15,25,14,14,14,13,28,1,9,2,21,20,30,17,13,6,19,23,12,1,13,23,24,28,26,12,19,22,6,28,18,19,8,17,28,11,20,12,19,21,12,13,20,26,12,13,18,17,22,16,24,15,6,22,30,25,22,15,4,20,11,5,18,25,21,30,29,2,11;];
Tprfs = [27,18,29,14,31,30,39,34,37,34,29,34,29,39,39,37,26,39,22,31,32,39,28,33,27,36,17,29,8,23,33,26,34,31,17,38,34,38,33,30,10,27,33,36,19,16,16,28,32,38,25,22,32,34,30,38,39,32,23,34,36,34,21,22,26,30,36,37,25,40,24,37,34,29,35,39,31,24,26,26,28,31,27,27,24,34,35,26,31,13,32,26,22,30,40,28,35,33,18,20;];
% 以上是标准起降时间表

Mgfs = 4.*ones(1,F);    % 地面最多等待4个时间段
Mafs = Mgfs./2;         % 空中最多等待2个时间段

Qfs = FSt.*20;         % 取消费用，大概相当于单位时间损失的10倍

StdTime = Tprfs - Tpdfs;   % 标准飞行时间，实际飞行时间不能比这个短

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Generation Process Complete

%%%%%%%%%%%%%%%%%%%%%%% Generating Chrom
NIND1=100;
Chrom11 = crtbp(NIND1,T.*ones(1,F));
for ii = 1:NIND1
    for jj = 1:F
        if Chrom11(ii,jj) <= Tpdfs(jj)
            Chrom11(ii,jj) = Tpdfs(jj);
        end
    end
end               % Repair chrom so that all the requirments are satisfied

Chrom12 = ceil(10.*rand(NIND1,F)) + Chrom11;
for ii = 1:NIND1
    Chrom12(ii,:) = Chrom12(ii,:) + StdTime;
end
Chrom = [Chrom11,Chrom12];

Chromzs = zeros(NIND1,F);    % 0 -- 不取消；1 -- 取消；

for ii = 1:NIND1
    for jj = 1:F
        if (Chrom11(ii,jj) > Tpdfs(jj) + Mgfs(jj))|(Chrom12(ii,jj) - Chrom11(ii,jj) > StdTime(jj) + Mafs(jj))
            Chromzs(ii,jj) = 1;                    % 取消航班的情况（延时过长）
        end
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Value of mincost
for i=1:NIND1
    TdfsM = zeros(F,T);
    TrfsM = zeros(F,T);
    for mm = 1:F
        for nn = 1:T
            if Chrom(i,mm) == nn
                TdfsM(mm,nn) = 1;
            end
            if Chrom(i,mm+F) == nn
                TrfsM(mm,nn) = 1;
            end
        end
    end
    ObjV(i,:) = mincost(FS,FSt,Chromzs(i,:),Cgf,Caf,Tpdfs,Tprfs,TdfsM,TrfsM,Qfs);
end

[R,S]=size(ObjV);
if R==1
    ObjV=ObjV';
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% feasible rate
for i=1:NIND1
    c=0;
    feasible(i)=0;
    yviol(i,:)=Avi_Voil(FS,FSt,Chromzs(i,:),Chrom(i,1:F),Chrom(i,F+1:2*F),Cgf,Caf,Tpdfs,Tprfs,Qfs,Mgfs,Mafs,FSu,FSd,F,T);
    if yviol(i,:) <= voilthres;
       c=c+1;
       feasible(i)=1;
    end
end

if c>0
    ppfirst=(NIND1-c)./NIND1;
    t=1;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

tic;

for gen=1:MAXGEN
    index=unidrnd(NIND1,2.*NIND1,1);        % 开始协同进化计算
    for n=1:NIND1
        if feasible(index(n))==1 & feasible(index(NIND+n))==1
            if ObjV(index(n))>ObjV(index(NIND+n))
                SelChmin(n,:)=Chrom(index(NIND+n),:);
            else
                SelChmin(n,:)=Chrom(index(n),:);
            end
        elseif feasible(index(n))==0 & feasible(index(NIND+n))==0
                if yviol(index(n))>yviol(index(NIND+n))
                    SelChmin(n,:)=Chrom(index(NIND+n),:);
                else
                    SelChmin(n,:)=Chrom(index(n),:);
                end
        elseif feasible(index(n))==1 & feasible(index(NIND+n))==0
            if yviol(index(NIND+n))<=esp
                if ObjV(index(n))>ObjV(index(NIND+n))
                    SelChmin(n,:)=Chrom(index(NIND+n),:);
                else
                    SelChmin(n,:)=Chrom(index(n),:);
                end
            else
                SelChmin(n,:)=Chrom(index(n),:);
            end
        elseif feasible(index(n))==0 & feasible(index(NIND+n))==1
            if yviol(index(n))<esp
                if ObjV(index(n))>ObjV(index(NIND+n))
                    SelChmin(n,:)=Chrom(index(NIND+n),:);
                else
                    SelChmin(n,:)=Chrom(index(n),:);
                end
            else
                SelChmin(n,:)=Chrom(index(NIND+n),:);
            end
        end
    end                         
    
    %for i=1:NIND
    %    ObjVSelmin(i,:)=CND_Obj(Chrom(i,:),nodes,cost);
    %end
    %[R,S]=size(ObjVSelmin);
    %if R==1
    %   ObjVSelmin=ObjVSelmin';
    %end
    SelChmin = recombin('xovsp',SelChmin,0.25);
    SelChmin = mut(SelChmin);
    
    for mm = 1:NIND1                       % Executing Simple Repair Algorithm
        for nn = 1:F
            if SelChmin(mm,nn) < Tpdfs(nn)
                SelChmin(mm,nn) = Tpdfs(nn);
            end
            if SelChmin(mm,nn+F) - SelChmin(mm,nn) < StdTime(nn)
                SelChmin(mm,nn+F) = SelChmin(mm,nn) + StdTime(nn);
            end
        end
    end        
    
    Chrom = reins(Chrom,SelChmin);
 
    for i=1:NIND                  % 计算每一个染色体对应的Objv和Viol
        c=0;
        feasible(i)=0;
        ObjV(i,:)=CND_Obj(Chrom(i,:),nodes,cost);
        yviol(i,:)=CND_Voil(Chrom(i,:),nodes,eff,rthres,nthres);
        if yviol(i)<=voilthres
           c=c+1;
           feasible(i)=1;
        end
    end
    
    [R,S]=size(ObjV);
    if R==1
       ObjV=ObjV';
    end
    if c>0
       pp(gen)=(NIND-c)./NIND;
       t=1;
    else
        pp(gen)=1;
    end
    
    if t==1 & mod(gen,5)==0
       if pp(gen-4)<=P & pp(gen-3)<=P & pp(gen-2)<=P & pp(gen-1)<=P & pp(gen)<=P
           esp=1.1.*esp;
       elseif pp(gen-4)>P & pp(gen-3)>P & pp(gen-2)>P & pp(gen-1)>P & pp(gen)>P
              esp=0.3.*esp;
       else
              esp=esp;
       end
    end
    
    if gen>500       % 选择本次的Chrom是什么
        a=1;
        ObjVNew=[];
       for i=1:NIND
           if (feasible(i)==1)
               ObjVNew(a,:)=CND_Obj(Chrom(i,:),nodes,cost);
               Chrom1(a,:)=Chrom(i,:);
               a=a+1;
           end 
       end
    else
        ObjVNew=ObjV;
        Chrom1=Chrom;
    end
    [Y,I]=min(ObjVNew);
    Out(gen,:)=Y;
    Chromoutp(gen,:)=Chrom1(I,:);
    disp(gen);
    toc
end

[j,k]=size(Chromoutp);
for i=1:j
    VoilNew(i,:)=CND_Voil(Chromoutp(i,:),nodes,eff,rthres,nthres);
end
ChromOut=[Chromoutp,Out,VoilNew];
x=1:MAXGEN;
plot(x,Out);