function reproductedInds = reproduction(indsDec, mutatedInds, elitismRate,y1,N)
% worst individuals from mutatedInds are replaced by best individuals
% from indsDec. The number of replaced individuals depends on elitismRate
    mutatedIndsDec=bin2dec(mutatedInds);
    % number of elite solutions to replace
    numEliteInds=round(N*elitismRate);
    
    % previous fitness
    fitnessOld=y1(indsDec);
    [fitnessOld i]=sort(fitnessOld,'descend');
    % avgFitnessOld=sum(fitnessOld)/N;
    eliteInds=indsDec(i(1:numEliteInds));
    
    % fitness of individuals after selection, recombination and mutation
    fitnessNew=y1(mutatedIndsDec);
    % avgFitnessNew=sum(fitnessNew)/N;
    
    [fitnessNew i]=sort(fitnessNew); 
%     fitnessNew=fitnessNew';
    i=i';
    % worst ind in new population
%     worstInds=mutatedIndsDec(i(1:numEliteInds));
%     mutatedIndsDec(i(1:numEliteInds))=[];
    for j=1:numEliteInds
        mutatedIndsDec(i(j))=eliteInds(j);
    end
    reproductedInds=mutatedIndsDec;
end