In Gibbs ensemble, we assume that along with energy exchanges (as heat), there is also an extensive parameter $x$ that changes in a field $f$ (so that $dE = dQ + fdx$) such that $E_{\text{sys}}+E_{\text{env}}=E_0$ and $x_{\text{sys}}+x_{\text{env}}=x_0$ where $E_0$ and $x_0$ are both constants.
Then we observe that the probability $P(\mu)$ of finding the system in a microstate $\mu$ is proportional to $\mathcal N_{\text{env}}(E_0-E_\mu, x_0-x_\mu)$.
Then using that the number of particles in the system is much smaller than that of those in the environment (to ignore the higher derivatives), and using that $\partial S/\partial E = 1/T$, we get that $P(\mu)$ is proportional to $$ \exp\left(-\frac{E_\mu}{k_BT} -\frac{1}{k_B}\frac{\partial S}{\partial x}x_\mu\right). $$
Now, all that remains to show is that $$ \frac{\partial S}{\partial x}= -\frac{f}{T} $$, and we'll get the celebrated $e^{-\beta(E_\mu-fx)}$. But I've no clue how to prove this.
I tried considering two systems labeled 1 and 2, isolated from the rest of the world, exchanging energy and $x$ such that \begin{align} E_1 + E_2 &= \text{const}, \quad\text{and}\\ x_1+x_2 &=\text{const}. \end{align}
Then we need to extremize $$ S_1(E_1, x_1) + S_2(E_2, x_2)-\lambda(E_1+E_2)-\mu(x_1+x_2), $$ which yields \begin{alignat}{2} \frac{\partial S_1}{\partial E_1} &= \frac{\partial S_2}{\partial E_2} &= \lambda &\stackrel{\text{def}}{=}\frac{1}{T},\quad\text{and}\\ \frac{\partial S_1}{\partial x_1} &= \frac{\partial S_2}{\partial x_2} &= \mu. \end{alignat}
But this was a dead end as I could not show $\mu=-f/T$ here.