%function y = Aviation_Lag(FS,FSt,zs,Cgf,Caf,Tpdfs,Tprfs,Qfs,Mgfs,Mafs,F,T,LagY,r)

%LagY = 10;
%r = 1.1;

FS = [];

F = 100;
T = 40;

FSu = ones(1,F);
FSd = 2*FSu;

for i = 1:F
    z = rand(1);
    if z > 0.7
        FSt(i) = 2;
    else
        FSt(i) = 1;
    end
end

Fdashs = [];

for i = 1:F
    z = rand(1);
    if z > 0.8
        zs(i) = 1;
    else
        zs(i) = 0;
    end
end

Cgf = 7;
Caf = 10;

for i = 1:F
Tpdfs(i) = ceil(rand(1).*10);
Tprfs(i) = ceil(T - 10.*rand(1));
end

Mgfs = T.*ones(1,F);
Mafs = Mgfs;

Qfs = FSt.*100;

% 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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create TdfsM/TrfsM Matrix Using standard GA
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
alpha = 1/3;   % Set a threshold for GA process
Cthres = sum(Qfs).*alpha;   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

StdTime = Tprfs - Tpdfs;
% Inner GA starts

NIND = 20;            % num of individuals
GGAP=1;               % generation gap
XOVR=0.7;               % rate of recombine
MUTR=1/4;          % rate of mutation
MAXGEN=10;           % Max number of generation
INSR=0.8;             % probabily of insertion
SUBPOP=1;             % number of subpop
MIGR=0.2;             % rate of migration
MIGGEN=20; % every 20 gen miggration
kk=1;

BaseVT = (T-StdTime-1).*ones(1,F);   % 計時時間段數從1開始，最大到T - StdTime,因航班而異，而這是起飛的時間段向量范圍；
ChromT = crtbp(NIND,BaseVT); % 隨機構造標準的起飛/降落時間表，嚴格滿足要求：先起飛后降落，且都在預定的整個時間段內完成飛行任務；
ChromT1 = ChromT + 1;
BaseVL = StdTime + Mgfs + Mafs;       % 這是降落的時間段向量范圍，要疊加在起飛時間段的基礎上
BX = max(ChromT1);
TempB = T.*ones(1,F) - BX;
for i = 1:F
    if BaseVL(i) > TempB(i)
        BaseVL(i) = TempB(i);
    end
end
ChromL = crtbp(NIND,BaseVL);
ChromL = ChromL + ChromT1;
ChromStd = [ChromT1,ChromL];
%zox = ChromStd(:,1+F:2*F) - ChromStd(:,1:F)

gen = 0;

TdfsM = zeros(F,T);
TrfsM = zeros(F,T);

for i = 1:NIND
    
% TdfsM/TrfsM--計劃航班序列的實際起飛/降落時間(為兩個矩陣，大小都為F行T列，每一行代表一個航班的起（降）情況，且每行只有一個1，代表航班在該時刻位置起（降）);
    for mm = 1:F
        for nn = 1:T
            if ChromT(i,mm) == nn
                TdfsM(mm,nn) = 1;
            end
            if ChromL(i,mm) == nn
                TrfsM(mm,nn) = 1;
            end
        end
    end
    ObjV(i,:) = mincost(FS,FSt,zs,Cgf,Caf,Tpdfs,Tprfs,TdfsM,TrfsM,Qfs);
end

cnt = 1;

%disp(min(ObjV));

while gen < MAXGEN
    [YYY,MN] = min(ObjV);
    rec(cnt) = YYY;
    cnt = cnt + 1;
    
    FitnV = ranking(ObjV,[2,1],SUBPOP);
    SelCh=select('sus',ChromStd,FitnV,GGAP,SUBPOP); % select
    SelCh=recombin('recdis',SelCh,XOVR,SUBPOP);  % recombin
    SelCh=mut(SelCh); % mutate, next generation is exist
    
    for ii = 1:NIND
        REC = SelCh(ii,1:F);
        for jj = 1:F
            if REC(jj) > T - StdTime(jj) - 1
                REC(jj) = T - StdTime(jj) - 1;
            end
        end
        SelCh(ii,1:F) = REC;
        REC = SelCh(ii,F+1:2*F);
        for jj = 1:F
            if REC(jj) > T
                REC(jj) = T;
            end
        end
        
        for mm = 1:F
            for nn = 1:T
                if SelCh(i,mm) == nn
                    TdfsM(mm,nn) = 1;
                end
                if SelCh(i,mm+F) == nn
                    TrfsM(mm,nn) = 1;
                end
            end
        end
        
        ObjVOff(i,:) = feval('mincost',FS,FSt,zs,Cgf,Caf,Tpdfs,Tprfs,TdfsM,TrfsM,Qfs);
    end
    
    [ChromStd, ObjV]=reins(ChromStd, SelCh, SUBPOP, [1,INSR],ObjV,ObjVOff); 
    
    gen = gen + 1;
end
        
%    Chromm1(jj,:)=Chrom(MN,:);
%          R0(jj,:)=YYY;
%          Chrommin=[Chromm1,R0];
%          jj=jj+1;
%          FitnV=ranking(ObjV,[2,1],SUBPOP);% assign fitness values
%          SelCh=select('sus',Chrom,FitnV,GGAP,SUBPOP); % select
%          SelCh=recombin('recdis',SelCh,XOVR,SUBPOP);  % recombin
%          SelCh=mutate('mutbga',SelCh, FieldDD,[MUTR],SUBPOP); % mutate, next generation is exist
%          ObjVOff=feval('findcal2forGA_Accurate_withAR_CRTRP_examine',SelCh(:,1),SelCh(:,2),SelCh(:,3),SelCh(:,4),SelCh(:,5),gen);   % object function value of next generation
%          [Chrom, ObjV]=reins(Chrom, SelCh, SUBPOP, [1,INSR],ObjV,ObjVOff); % replace
%          if (rem(gen, MIGGEN)==0)
%              [Chrom,ObjV]=migrate(Chrom, SUBPOP, [MIGR,1,1],ObjV);    % every 20 generations, migrate
%          end

    

% 任意航班降落的時間總是落后于起飛的時間
% 實際飛行時間不得小于預定飛行時間

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

yobj = min(rec);

ylagrangian = 0;

% The Taking-off-Capacity Constraint

% serial_to = 10;

for f = 1:F
    for t = 1:T
        recc = OffCon(FS,zs,TdfsM,FSu,f,t);
        yc = TdfsM(f,t) - recc;                       % 拉格朗日乘子 -- 起飛容量約束，任意機場任意時刻都要加上，yc >= 0;
        ylagrangian = ylagrangian + r.*(min(0,yc+LagY./r).^2-(LagY./r).^2);
    end
end

for f = 1:F
    for t = 1:T
        recc = LandCon(FS,zs,TrfsM,FSd,f,t);
        yc = TrfsM(f,t) - recc;                       % 拉格朗日乘子 -- 降落容量約束，任意機場任意時刻都要加上，yc >= 0;
        ylagrangian = ylagrangian + r.*(min(0,yc+LagY./r).^2 - (LagY./r).^2); 
    end
end

for f = 1:F
    rec = TimeCon(zs,TdfsM,TrfsM,Tpdfs,Tprfs,f);     % 拉格朗日乘子 -- 飛行時間約束，所有航班都要加上，
    ylagrangian = ylagrangian + r.*(min(0,rec+LagY./r).^2 - (LagY./r).^2);
end