Use of generating function in canonical transformation In the theory of Canonical transformations, initially we use the fact that the new and the old system of $(q_i, p_i)$ with the Hamiltonian $H$ satisfy the modified Hamilton's principle. Now here, the use of the theory of generating functions makes sense. But after the development of Poisson bracket formalism of canonical transformations, what exactly is the use of generating function formalism? As such even, if it was tough to guess the transformation equation $$P_i=P_i(q_i,p_i,t) ; Q_i=Q_i(q_i,p_i,t)$$ as such, how would it be any easy to guess the generating function?
 A: In canonical perturbation theory, the generating function (of the $F_2$-type) is determined order by order from the perturbative term, i.e. if 
$H(q,p,t)=H_0(q,p)+\epsilon H_1(q,p,t)$ then one can write a series expansion
$$
S= qP+ \epsilon S^{(1)}(q,P,t)+\epsilon^2 S^{(2)}(q,P,t)+\ldots
$$
and determine $S^{(1)}$ given $H_1(q,p,t)$ and $H_0(q,p)$.  Note that $S^{(1)}$ usually involves some averaging process.  From $S^{(1)}$ one can then determine the corrections to the unperturbed motion (to first order in this case) and so forth to higher order.  It gets messy but it's doable.
Possibly the cleanest example is seen in the action-angle formalism.  If $K(I)$ is the transformed Hamiltonian in the new action-angle coordinates $(I,\theta)$, we have
\begin{align}
K(I)&=K_0(I)+\epsilon K^{(1)}(I)+\ldots\, ,\\
H(\phi,J)&= H_0(J)+\epsilon H^{(1)}(\phi,J)
\end{align}
with $(\phi,J)$ the "old" action angle variables.  Then 
$$
S^{(1)}(\phi,I)=\displaystyle\int_0^{2\pi} 
\frac{K^{(1)}(I)-H^{(1)}(I,\phi)}{\omega_0(I)}d\phi
$$
with $K^{(1)}(I)$ the average of the perturbation
$$
K^{(1)}(I)=\frac{1}{2\pi}\int_0^{2\pi} d\theta H^{(1)}(\theta,I)\, .
$$
expressed in terms of the new variables, with $J=I, \phi=\theta$ to leading order in $\epsilon$.
A good discussion of this can be found in the book by Jose and Saletan, Classical Dynamics.  Also some examples with one position and one momentum can be found in the book by Percival and Richards Introduction to dynamics.
A: *

*For various notions of canonical transformations (CTs), see this Phys.SE post.


*Presumably what OP calls "Poisson bracket formalism of CTs" refers to symplectomorphisms (at least if there is no explicit time dependence).
If we furthermore assume that the $2n$ new coordinates $(Q^1,\ldots, Q^n, P_1, \ldots, P_n)$ are canonical/Darboux coordinates, then we get $n (2n-1)$ differential conditions$^1$
$$ \{Q^i,P_j\}~=~\delta^i_j , \qquad\{Q^i,Q^j\} ~=~0~=~\{P_i,P_j\}  ,\qquad i,j~\in~\{1,\dots,n\},\tag{A}$$
to satisfy in terms of the old canonical/Darboux coordinates $(q^1,\ldots, q^n, p_1, \ldots, p_n)$.


*On the other hand, if we instead use a generating function $F$ to define a CT, we are automatically guaranteed that the new coordinates $(Q^1,\ldots, Q^n, P_1, \ldots, P_n)$ are canonical/Darboux coordinates. Also it is often much simpler to solve for a generating function $F$ than to solve the full equation system (A). This is one possible answer to OP's question.
--
$^1$ The number of differential conditions corresponds to the number of independent entries in an antisymmetric $2n\times 2n$ matrix.
