function [popu] = BoundaryDetection(popu,lu)
[NP, D] = size(popu);
xl = repmat(lu(1,:), NP, 1);
xu = repmat(lu(2,:), NP, 1);
ranX=rand(NP, D);
%check the lower bound
pos = popu < xl;
popu(pos) = xl(pos) + ranX(pos) .* (xu(pos) - xl(pos));

%check the upper bound
pos = popu > xu;
popu(pos) = xl(pos) + ranX(pos) .* (xu(pos) - xl(pos));
end