coga.m

遗传算法例程：通过MATLAB

% FS--计划航班编号序列，大小1*F（函数中没有出现，可以不用）,FSt--计划航班类型序列，大小1*F(2表示大型机,1表示中型机,没有小型机);
% Fdashs--连续飞行的航班序列（为F的一个子集,表示为一个2进制串（大小1*F'),第i个位置为1表示连续航班序列包括第i个航班); --本函数没有用到
% FSu--计划航班序列的各个起飞机场,;FSd--计划航班序列的各个降落机场; --本函数没有用到，二者大小都是1*F；
% zs--各个航班是否取消的向量,为一个二进制串,大小1*F，第i个位置为1表示第i个航班取消,为0表示第i个航班不取消(这是一个决策向量);
% Cgf--地面延误的单位时间成本(对大型机); Caf--空中延误的单位时间成本(对大型机)，这两个都是单一数字
% Tpdfs/Tprfs--这是计划航班序列的计划起飞和降落的时间段(都是1*F向量);
% Mgfs/Mafs--计划航班序列的依次最大地面/空中的等待时间数(都是1*F向量);--本函数没有用到
% TdfsM/TrfsM--计划航班序列的实际起飞/降落时间(为两个矩阵，大小都为F行T列，每一行代表一个航班的起（降）情况，且每行只有一个1，代表航班在该时刻位置起（降）);
% Qfs--计划航班序列的取消成本，为1*F向量

% Only zs plays the role of j_v in this function 
% F -- Num of total flights
% T -- Num of total time periods
% LagY & r -- Lagrangian Multiplexer

%%%%%%%%%%%%%%%%%%%%%%% Generating Test Data
FS = [];

T = 40;

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

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

% Generate zs according to TdfsM & TrfsM,  

Cgf = 7;
Caf = 10;

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

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

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

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Generating Process Complete

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
Chrom1 = [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

NIND2=20;
PRECI=10;
FieldD2=[rep([PRECI],[1,1]); rep([0;10],[1,1]);rep([1;0;1;1],[1,1])];
Chrom2=crtbp(NIND2,PRECI);
Chrom1u=Chrom1(5,:);
Chrom2u=Chrom2(1,:);
r=1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
MAXGEN=80;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for i=1:NIND1
%    reqmin(i)=T_P2Lag(Chrom1(i,:),bs2rv(Chrom2u,FieldD2),r);
    reqmin(i) = Avi_Lag(FS,FSt,Chromzs(i,:),Chrom11(i,:),Chrom12(i,:),Cgf,Caf,Tpdfs,Tprfs,Qfs,Mgfs,Mafs,FSu,FSd,F,T,bs2rv(Chrom2u,FieldD2),r);
end
[R,S]=size(reqmin);
if R==1
    reqmin=reqmin';
end
%disp(reqmin);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%disp('part1 done');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for j=1:NIND2
    reqmax(j)=Avi_Lag(FS,FSt,Chromzs(5,:),Chrom11(5,:),Chrom12(5,:),Cgf,Caf,Tpdfs,Tprfs,Qfs,Mgfs,Mafs,FSu,FSd,F,T,bs2rv(Chrom2(j,:),FieldD2),r);
end
[M,N]=size(reqmax);
if M==1
    reqmax=reqmax';
end

%disp(Qfs);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%disp('part2 done');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Field = rep([1;T],[1,F]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% CoGA Starts
for gen=1:MAXGEN
    %tic;
    for gen1=1:20
        % GA begin
        FitnVmin=ranking(reqmin);
        SelChmin=select('rws',Chrom1,FitnVmin);
        SelChmin=recombin('xovsp',SelChmin,0.8);
        SelChmin=mutate('mutbga',SelChmin,[Field,Field],[1/4]); 
        SelChmin = round(SelChmin);
        
        for mm = 1:NIND1                       % 修补染色体
            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
        
        Chrom1=reins(Chrom1,SelChmin);
        Chrom11 = Chrom1(:,1:F);
        Chrom12 = Chrom1(:,1+F:2*F);

        for ii = 1:NIND1            % 生成新的zs
            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;                    % 取消航班的情况（延时过长）
                else
                    Chromzs(ii,jj) = 0;
                end
            end
        end
        
        for i=1:NIND1
            reqmin(i)=Avi_Lag(FS,FSt,Chromzs(i,:),Chrom11(i,:),Chrom12(i,:),Cgf,Caf,Tpdfs,Tprfs,Qfs,Mgfs,Mafs,FSu,FSd,F,T,bs2rv(Chrom2u,FieldD2),r);
        end
        [R,S]=size(reqmin);
        if R==1
           reqmin=reqmin';
        end
        %t1=t1.*0.85; % can adjust 0.85
    end
    [dmin,ii]= min(reqmin);
    
    Chrom1u=Chrom1(ii,:);
    Chrom11u = Chrom1u(:,1:F);
    Chrom12u = Chrom1u(:,F+1:2*F);
    
    for jj = 1:F
        if (Chrom11u(jj) > Tpdfs(jj) + Mgfs(jj))|(Chrom12u(jj) - Chrom11u(jj) > StdTime(jj) + Mafs(jj))
            Chromzsu(jj) = 1;                    % 取消航班的情况（延时过长）
        else
            Chromzsu(jj) = 0;
        end
    end
    
    for gen2=1:10
        FitnVmax=ranking(-reqmax);
        SelChmax=select('rws',Chrom2,FitnVmax);
        SelChmax=recombin('xovsp',SelChmax,0.8);
        SelChmax=mut(SelChmax);

        Chrom2=reins(Chrom2,SelChmax);
        
        for j=1:20
            reqmax(j)=Avi_Lag(FS,FSt,Chromzsu,Chrom11u,Chrom12u,Cgf,Caf,Tpdfs,Tprfs,Qfs,Mgfs,Mafs,FSu,FSd,F,T,bs2rv(Chrom2(j,:),FieldD2),r);
        end
        
        [M,N]=size(reqmax);
        if M==1
        reqmax=reqmax';
        end
        %t2=t2.*0.85;
    end
    [dmax,jj]=max(reqmax);
    Chrom2u=Chrom2(jj,:);
    
    r=r.*1.1;
    if r>10000
        r=10000;
    end
    
    x=Chrom1u;

    disp(gen);
    
    TdfsMr = zeros(F,T);
    TrfsMr = TdfsMr;
    for mm = 1:F
        for nn = 1:T
            if Chrom11u(mm) == nn
                TdfsMr(mm,nn) = 1;
            end
            if Chrom12u(mm) == nn
                TrfsMr(mm,nn) = 1;
            end
        end
    end
    
    Obj(gen) = mincost(FS,FSt,Chromzsu,Cgf,Caf,Tpdfs,Tprfs,TdfsMr,TrfsMr,Qfs);
    %toc;
end

mm=1:MAXGEN;
plot(mm,Obj);