Why is chemical potential constant for the system in Gibbs distribution? Consider the grand canonical Gibbs distribution where a system is in thermal and diffusive equilibrium with a large reservoir. 
The system can take on different energies $E_i$ and contain different number of particles $N_i$ when it is in different accessible microstates. 
The definition of chemical potential for the system is 
$$\mu = -T\left(\frac{\partial S}{\partial N}\right) $$
where $T$, $S$, $N$ refers to the temperature, entropy and number of particles of the system.
If the system in the Gibbs distribution take on different $E_i$ and $N_i$, why is $\mu$ the same? For example if system is in $E_1$, $N_1$ state it will have a corresponding $S_1$ and a calculated $\mu_1$ using the formula above. If system is now in $E_2$, $N_2$ state it will now have entropy $S_2$ and $\mu_2$. 
 A: The answer is essentially the same as why a system in contact with a thermal bath is at the temperature of the bath, or a system subject to an external pressure will be at that pressure.  With a thermal bath, the system is allowed to exchange energy with the environment; with external pressure, the system is allowed to exchange volume with the environment.  In the case of diffusive equilibrium with a reservoir of particles, the system is allowed to exchange particles with the external environment.
In all cases, the equilibrium will occur when the whole system+reservoir is in the maximum state of entropy.  Taking particles out the reservoir decreases the entropy of the reservoir, while increasing the entropy of the system of interest.  At equilibrium, the two are exactly in balance.  This works exactly the same way as insisting that changing the volume of system leads to no change in $S$ for the combined whole.  The math showing that the pressure of the system is the reservoir pressure is exactly the same as for showing that the chemical potential of the system is the reservoir value.  This kind of analysis should be found in almost any statistical mechanics text—for example, Fundamentals of Statistical and Thermal Physics by F. Rief (section 3.9, pages 114–116).
A: In Thermodynamics and in Statistical Mechanics, many routes may reach the same result. However, mixing them, even if the final result is a correct formula, may leave some doubts about the consistency of the arguments with the starting hypothesis and definition.
In this specific case, I see a potential conceptual  problem in mixing the first statement of the question:

Consider the grand canonical Gibbs distribution...

and the formula for the chemical potential

$$ \mu = -T \frac{\partial{S}}{\partial{N}}. $$

Such a formula is correct from the thermodynamic point of view, but, to be meaningful in the present context, one should specify which state variables are kept constant in the partial derivative, i.e., which variables $S$ is a function of. Here, some ambiguity of the description enters.
Indeed, it is true that one traditional way of deriving and interpreting the grand canonical ensemble is through the analysis of a "big"  system with a finite number of particles, partitioning it into a "system of interest" with a fixed subvolume $V$ of the total volume, considering all the possible constraints corresponding to fix the values of the number of particles of the subvolume, and looking for the equilibrium condition with respect to the exchange of particles at a fixed temperature.
Following this way of deriving the grand canonical ensemble, or even through a direct thermodynamic analysis, it turns out that the above formula for chemical potential should involve a partial derivative of $S$ with respect to $N$ at fixed $T,V$ (or with respect to $E,V$; $E$ being the internal energy of the system). Thus, one should know $S(T,V,N)$ (or $S(E,V,N)$) for the system of interest.
However, once the grand canonical formulae have been obtained, the grand canonical partition function $\Xi $ and the resulting Grand Potential $\Omega = -k_BT \log \Xi$ are functions of $T,V,\mu$. The number of particles is an ill-defined concept in the grand-canonical ensemble. It cannot be fixed since it fluctuates. The average number of particles $<N>$ should replace it, which can be obtained as $-1/k_BT$ times the partial derivative of $\Omega$ with respect to $\mu$ at fixed $T,V$.
Said in another way, a formula like  $$ \mu = -T \left. \frac{\partial{S}}{\partial{N}}\right|_{T,V} = -T \left. \frac{\partial{S}}{\partial{N}}\right|_{E,V} $$ looks like a foreigner in the grand canonical ensemble because $N$ is not a well defined quantity (neither $E$, referring to the second form). On the other side, $\mu$, exactly like $T$ and $V$ is an independent physical parameter characterizing the ensemble and the thermodynamic macrostate one is interested in.
Put this way, the question looks like a similar question which could be asked for the canonical ensemble: since we know that $\frac{1}{T} = \frac{\partial{S}}{\partial{E}}$, if the system may take different energies, and then different entropies, why the temperature is unique?
Once one has already gone through the derivation of the canonical ensemble, the question, in this form, is not well-posed because energy is not an independent variable anymore. Missing this point may drive towards even worse ill-posed questions (sometimes appeared even in the scientific literature) like wondering about temperature fluctuations in the canonical ensemble or fluctuations of chemical potential in the grand-canonical ensemble. By construction, both quantities are fixed, and speaking about their fluctuations is meaningless.
