Let us consider a grand canonical ensemble of $M$ identical systems such that the total number of particles in the ensemble is $M\bar{N}$ and a total energy of $M\bar{E}$. Let $n_{r,s}$ be the number of systems having energy $E_{s}$ and number of particles $N_{r}$. Putting these statements mathematically, we have
$\sum_{r,s} n_{r,s}= M$
$\sum_{r,s} n_{r,s} N_{r}=M\bar{N}$
$\sum_{r,s} n_{r,s} E_{s}=M \bar{E}$
The above part lays the groundwork. The solutions begins here.
Let $P_{r,s}$ be the probability that a given member of the ensemble has an energy $E_{s}$ and number of particles $n_{r}$. From basic considerations of the grand canonical ensemble, we have
$P_{r,s}=\frac{exp(-\alpha N_{r}-\beta E_{s})}{\sum_{r,s}exp(-\alpha N_{r}-\beta E_{s})}$ .............. (1)
Now, the denominator is the grand canonical partition function $Z_{G}$
$Z_{G}=\sum_{r,s}exp(-\alpha N_{r}-\beta E_{s})$ .............. (2)
Therefore,
$P_{r,s}=\frac{exp(-\alpha N_{r}-\beta E_{s})}{Z_{G}}$ .............. (3)
Now, from the definition of entropy, we have
$S=-k_{B}<ln(P_{r,s})>$ .............. (4)
$=> S=-k_{B}\sum_{r,s} P_{r,s} ln(P_{r,s}) $ .............. (5)
Now, we substitute the expression of $P_{r,s}$ obtained in equation (3) inside the $ln$. We leave the $P_{r,s}$ outside the $ln$ unchanged. (This makes the calculations easier as we will see shortly).
$S=-k_{B} \sum_{r,s} P_{r,s} ln[\frac{exp(-\alpha N_{r}-\beta E_{s})}{Z_G}]$ .............. (6)
$=>S=-k_{B} \sum_{r,s} P_{r,s}[-\alpha N_{r}-\beta E_{s}-ln(Z_{G})]$ .............. (7)
$=>S=k_{B} \alpha \sum_{r,s} P_{r,s}N_{r}+k_{B} \beta \sum_{r,s} P_{r,s} E_{s}+k_{B}ln(Z_{G}) \sum_{r,s} P_{r,s} $ .............. (8)
From normalization of probabilities,
$\sum_{r,s} P_{r,s}=1$.
From the basic definition of mean, we have
$\sum_{r,s} P_{r,s}N_{r}=\bar{N}$
and
$\sum_{r,s} P_{r,s}E_{s}=\bar{E}$
Therefore, equation (8) can be written as
$S=k_{B} \alpha \bar{N}+k_{B} \beta \bar{E}+k_{B}ln(Z_{G}) $ .............. (9)
Now, using the first law and second law of thermodynamics,
$d\bar{E}=TdS-PdV+\mu d\bar{N}$ .............. (10)
Since $S$, $V$, and $N$ are extensive quantities, using Euler's homogeneous function theorem, we get
$\bar{E}=TS-PV+\mu\bar{N}$ .............. (11)
Rearranging the terms of (11) and using the fact that $\beta=\frac{1}{k_{B}T}$, we get
$S=k_{B}\beta\bar{E}+k_{B}\beta PV - k_{B}\beta \mu \bar{N} $ .............. (12)
Comparing (9) and (12), we get
$k_{B}\beta PV =k_{B}ln(Z_{G})$
$=>\beta PV=ln(Z_{G})$
$=>\frac{PV}{k_{B}T}=ln(Z_{G})$ .............. (13)
[Hence proved]
Note : The key here is to bridge statistical mechanics, and thermodynamics. Entropy is often a useful tool for this purpose as we saw here.