function [pop,pop_NA,pop_indices] = mutation(pop,pop_NA, pop_indices, pop_size, mutation_rate,wf)

for i = 1:pop_size
    if rand > mutation_rate
        while 1
            turbine_pos = randi(wf.cols * wf.rows);
            if pop_NA(i, turbine_pos) == 1
                break
            end
        end
        while 1
            null_turbine_pos = randi(wf.cols * wf.rows);
            if pop_NA(i, null_turbine_pos) == 0
                break
            end
        end

        pop(i, turbine_pos) = 0;
        pop(i, null_turbine_pos) = 1;

        pop_NA(i, turbine_pos) = 0;
        pop_NA(i, null_turbine_pos) = 1;

        for j =1:wf.turbine_num
            if pop_indices(i, j) == turbine_pos
                pop_indices(i, j) = null_turbine_pos;
                break
            end
        end
        pop_indices(i, :) = sort(pop_indices(i, :));
    end
end
end