function [Fbest,BestChart,BestFitness,farmlayout,farmlayout_NA]=IWO(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,:);
% C = S = 1.5,a1 = a2  1, P [[0.8,1], FL[[0.5,0.9], FQ ¼ 3

Smin = 2;       % Minimum Number of Seeds
Smax = 5;       % Maximum Number of Seeds
Exponent = 2;           % Variance Reduction Exponent
sigma_initial = 0.5;    % Initial Value of Standard Deviation
sigma_final = 0.001;




for it = 1:iterations
    sigma = ((iterations - it)/(iterations - 1))^Exponent * (sigma_initial - sigma_final) + sigma_final;
    Bestfitness = max(pBestFit);
    Worstfitness = min(pBestFit);
    offspringpopu = [];
    offfitness = [];
    for i = 1:popsize
        ratio = (pBestFit(i) - Worstfitness) / (Bestfitness - Worstfitness);
        S = floor(Smin + (Smax - Smin) * ratio);
        for j = 1:S
            % Generate Random Location
            offsPbest = popu(i) + sigma * randn([1,D]);
            % Apply Lower/Upper Bounds
            offsPbest = windfarm_constraint(offsPbest, wf.NA_loc, D,down,up);
            % Evaluate Offsring
            [offsPbestFitness,~,~]= wf_fitness(wf,offsPbest);
            offspringpopu  = [offspringpopu; offsPbest];
            offfitness = [offfitness; offsPbestFitness];
        end

    end
    popu = [popu; offspringpopu];
    pBestFit = [pBestFit; offfitness];
    [pBestFit,sort_index] = sort(pBestFit,'descend');
    popu =  popu(sort_index,:);
    if size(popu,1) > popsize
        popu = popu(1:popsize,:);
        pBestFit = pBestFit(1:popsize);
    end

    %% Update pBest
    [gBestFit,index]=max(pBestFit);
    gBest = popu(index,:);
    %% Survival

    Fbest = gBestFit;
    BestChart(it) = gBestFit / wf.power_total;
    BestFitness(it) = gBestFit;
    [best_farmlayout,best_farmlayout_NA]  = gene_layout_by_indices_one(wf,gBest);
    farmlayout(it,:) = best_farmlayout;
    farmlayout_NA(it,:) = 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,it,(Fbest / wf.power_total),Fbest)
end
end




