Why we develop formalism of Bose Einstein Condensation in framework of grand canonical ensemble ?
When dealing with a collection of particles, there are 2 ways we can represent the state of the entire system. We can take each particle and specify what state it is in or we can take each state and specify how many particles are in that state. This second representation is known as the occupation number representation.
Now when are dealing with identical bosons (or respectively fermions) we have an additional requirement that the state of the system must be (anti)symmetric under particle exchange. In the standard representation this means that after specifying the state of each particle we must perform an appropriate (anti)symmeturization procedure. For a large number of particles this calculation can become very involved and cumbersome. In the occupation number representation, on the other hand, we can simply limit our attention to the occupation numbers of appropriately (anti)symmetric states and so essentially take care of the symmeturization automatically. This dramatically simplifies many calculations.
The only complication with the occupation number representation is what happens when we try to impose a constraint on the total number of particles. If I have ten particles and I know that the occupation number of state $A$ is 1, then I know that the occupation number of state $B$ cannot be 10, but if state $A$ is unoccupied then all the particles can be in state $B$. In other words the constraint means that different occupation numbers are not independent of each other, since they must sum to a fixed total. In the extreme case where I have only one particle, then as soon as I know that some state has an occupation number of 1, I immediately know all the other occupation numbers. This correlation between occupation numbers makes calculation much more complicated.
Therefore, when dealing with identical particles, it is far simpler to work in the grand canonical ensemble, where there is no hard constraint on particle number, than in the canonical ensemble, where there is. We can simply set our chemical potential so that the mean particle number comes out the way we want it too and significant deviations from this value are exponentially suppressed and so do not significantly contribute. In the thermodynamic limit the fluctuations in particle number should go to zero and so the calculations in either ensemble should agree.