A solution to your exercise is obtained by using the method of lagrange multipliers.
The constraints we have to satisfy are three:
1) $max_{P_j} S[p]$
2) $\sum_{j} P_j E_j = E$
3) $\sum_{j} P_j =1$
Choosing as variation parameters $\lambda$ and $\gamma$, calculations follow as:
$S_{\lambda,\gamma}=-\sum_{j} P_j ln{P_j} + \lambda(\sum_{j} P_j E_j - E)+\gamma(\sum_{j} P_j -1)$
$\frac{ \partial S_{\lambda,\gamma} }{ \partial{P_j} }=0\Longrightarrow -ln{P_j}-1+\lambda E_j+\gamma=0 \Longrightarrow P_j=\frac{ e^{- \lambda E_j} }{ e^{1-\gamma} }$
Calling $Z$ the normalization constant we have:
$P_j= \frac{ e^{-\lambda E_j} }{Z}$
$\sum_{j} P_j=1 \rightarrow e^{1-\gamma}=\sum_{j}e^{-\lambda E_j}:=Z $
The values of the $P_j$ maximize entropy since:
$\frac{ \partial^2 S_{\lambda , \gamma}}{ \partial P^2_j}=\frac{-1}{P_j}$
The parameter $\lambda$ can be determined from:
$\frac{ \sum_{j} E_j e^{-\lambda E_j} }{ Z } = E$
Now i will try to show how the parameter $\lambda$ can be interpreted as the inverse of temperature. Let's introduce the function $F(\lambda , E_j ):=ln{Z}$. We have:
$dF=\frac{\partial F}{\partial \lambda} d\lambda + \sum_{j} \frac{\partial F}{\partial E_j} dE_j=-Ed\lambda-\lambda \sum_{j} \frac{N_j}{N} dE_j $ where $P_j=\frac{N_j}{N}$
The last formula can be rewritten as:
$d(F+E\lambda)=\lambda(dE-\sum_{j} \frac{N_j}{N} dE_j):=\lambda dQ$
The last association follows from the physical fact that (looking from the point of view of quantum mechanics if you like) $-\sum_{j} \frac{N_j}{N} dE_j$ can be interpreted as the work done on the system to vary the energy levels from $E_j$ to $E_j+dE_j$ and $dE$ is the variation of internal energy. So $dE-\sum_{j} \frac{N_j}{N} dE_j$ is the quantity of heat $dQ$ exchanged by the ensemble with the outer.
Since we only have an exact differential involving heat from thermodynamics we can conclude that $\lambda=1/T$
Now we can define:
$f:=\frac{-1}{\beta} ln{Z}=E-TS$
So the maximum values of the entropy are embedded in the definition of free energy and they determine the minimum value of it.
The same procedure applies if we start from the free energy:
Constraints: $E=\sum_{j} P_j E_j$ and $\sum_{j} P_j = 1$
Parameter as lagrange multiplier: $\mu$
$f_\mu=\sum_{j}P_jE_j+T\sum_{j}P_j ln(P_j)+\mu(\sum_{j} P_j -1)$
$\frac{\partial f}{\partial P_j}=E_j + T ln{P_j} + T + \mu =0 \Longrightarrow P_j=e^{\frac{-E_j}{T}} e^{\frac{\mu}{T}-1} $
$\sum_{i}P_i=1 \rightarrow e^{\frac{-\mu}{T}-1}=\frac{1}{\sum_{i} e^{\frac{-E_i}{T}} }:=\frac{1}{Z} \Longrightarrow P_i= \frac{e^{\frac{-E_i}{T}}}{Z}$