function [Fbest,BestChart,BestFitness,farmlayout,farmlayout_NA]=GWO(popsize,wf,iterations,t,NA_type,tn,wt,algorithmDir)

BestChart=zeros(iterations,1);
BestFitness=zeros(iterations,1);
farmlayout = zeros(iterations,wf.rows*wf.cols);
farmlayout_NA = zeros(iterations,wf.rows*wf.cols);
pop_power = zeros(popsize,wf.cols * wf.rows);

D = wf.turbine_num; down = 1;up = wf.rows *wf.cols;

% Population initialization
[popu,lu] = windfarm_init(popsize, wf.turbine_num,wf);
% Population evaluation
[popuFitness,power_order,lp_power_accum]= wf_fitness(wf,popu);
optimal = max(popuFitness);
for i = 1:popsize
    pop_power(i,power_order(i,:))  = pop_power(i,power_order(i,:)) + (lp_power_accum(i,:) / sum(lp_power_accum(i,:)));
end





% Local parameters setting
vel = zeros(popsize, D);
pBest = popu;
pBestFit = popuFitness;
pBest_power = power_order;
[~, gBestId] = max(pBestFit);
gBest = pBest(gBestId,:);
Alpha_score = 0;
Alpha_pos = zeros(1, D);
Beta_pos = zeros(1,D);
Beta_score = 0;
Delta_pos = zeros(1,D);
Delta_score = 0;

for i = 1 : popsize
    if popuFitness(i) > Alpha_score
        Alpha_score = popuFitness(i);
        Alpha_pos = popu(i, :);
    end
    if popuFitness(i)<Alpha_score && popuFitness(i)>Beta_score
        Beta_score = popuFitness(i); % Update beta
        Beta_pos = popu(i,:);
    end
    if popuFitness(i)<Alpha_score && popuFitness(i)<Beta_score && popuFitness(i)>Delta_score
        Delta_score = popuFitness(i); % Update delta
        Delta_pos = popu(i,:);
    end
end

for iter = 1:iterations
    %% Search machine
    [popu_offspring]=GWO_operator(popsize,popu,1,144,iterations,Alpha_pos,Beta_pos,Delta_pos,iter);
    popu_offspring = windfarm_constraint(popu_offspring, wf.NA_loc, D,down,up);
    [offsPbestFitness,offs_power_order,offs_lp_power_accum]= wf_fitness(wf,popu_offspring);

    for i = 1 : popsize
        if popuFitness(i) > Alpha_score
            Alpha_score = offsPbestFitness(i);
            Alpha_pos = popu_offspring(i, :);
        end
        if offsPbestFitness(i) > Alpha_score
            Alpha_score = offsPbestFitness(i);
            Alpha_pos = popu_offspring(i, :);
        end
        if offsPbestFitness(i)<Alpha_score && offsPbestFitness(i)>Beta_score
            Beta_score = offsPbestFitness(i); % Update beta
            Beta_pos = popu_offspring(i,:);
        end
        if offsPbestFitness(i)<Alpha_score && offsPbestFitness(i)<Beta_score && offsPbestFitness(i)>Delta_score
            Delta_score = offsPbestFitness(i); % Update delta
            Delta_pos = popu_offspring(i,:);
        end
    end
    
    
    

    %% Update pBest
    temp_fitness = zeros(size(offsPbestFitness));
    popu1 = zeros(size(popu));

    pos = offsPbestFitness > popuFitness;
    temp_fitness(pos) = offsPbestFitness(pos);
    temp_fitness(~pos) = popuFitness(~pos);

    popu1(pos,:) = popu_offspring(pos,:);
    popu1(~pos,:) = popu(~pos,:);
    [gBestFit, gBestId] = max(temp_fitness);
    gBest = popu1(gBestId,:);
    [popuFitness, ~] = max([popuFitness, offsPbestFitness], [], 2);
    
    popu = popu_offspring;
    %% Survival

    Fbest = gBestFit;
    BestChart(iter) = gBestFit / wf.power_total;
    BestFitness(iter) = gBestFit;
    [best_farmlayout,best_farmlayout_NA]  = gene_layout_by_indices_one(wf,gBest);
    farmlayout(iter,:) = best_farmlayout;
    farmlayout_NA(iter,:) = best_farmlayout_NA;
    fprintf('%s NA %d Turbine Num:%d Wind %s run: %d iteration: %d  eta %f fitness %f\n',algorithmDir,NA_type,tn,wt,t,iter,(Fbest / wf.power_total),Fbest)
end

end