Why do we use grand canonical ensemble for BEC description?

The system we consider has constant $N$, $V$ and $T$ (the number of particles, volume and temperature) This is just the thermodynamic variables for the canonical ensemble, why we use fugacity $z$ or chemical potential $\mu$ there?

In the thermodynamic limit fluctuations in particle number should tend to $$0$$ and so we should get the same result whether we use the cannonical or grand canonical ensamble. If we are ultimately interested in a system with a fixed particle number, we can, if we so choose, set up the problem in the grand cannonical ensamble and then invert the equation for $$\langle N \rangle$$ to find what chemical potential will give us the desired mean particle number.
In the grand cannonical ensamble we do not encounter this problem, since all values of $$N$$ are allowed. We simply consider all possible sets of occupation numbers, but since states with very large or very small particle number are exponentially weighted against, they do not contirbute. In effect we allow the statistics of the problem to "automatically" find the states with the correct particle number for us.
If you choose to set $N$, $V$ and $T$, then you are choosing the Canonical Ensemble, so $N$ is fixed. In this approach, the chemical potential $\mu$ is just a parameter that you've to determine by imposing that the total number of particles is $N$.
Viceversa, if you choose the grandcanonical ensemble, you fix $\mu$ and so the number of particles $N$, fluctuates around an average value $\langle N \rangle$. I hope that I've well understood your problem and answered properly.