How to apply two or more gauge transformations? Let $f=f(t)$ be a physical quantity of a system, $t$ is a variable e.g time. For infinitesimal variation $\delta t$:
$$f(t_{0}+\delta t)=f(t_{0})+\frac{\mathrm{d}f}{\mathrm{d}t}\delta t=f_{0}+\{f,H_{T}\}$$
where $H_{T}=H+u_{m}\phi_{m}$ with $\phi_{m}\approx 0$ is primary constraint of the system. Afterward, $f(t_{0}+\delta t)$ can be written as:
$$f(t+\delta t)=f_{0}+\{f,H\}\delta t + u_{m}\{f,\phi_{m}\}\delta t$$
This is the way we apply gauge transformation with the infinitesimal change of $t$.
I think, if we apply the gauge transformation one more time, it will look like:
$$f'=f_{0}+\{f(t_{0}+\delta t),H\}\delta t + u_{m}\{f(t_{0}+\delta t),\phi_{m}\}\delta t$$
Then the procedure will be continued with expansion of $f(t_{0}+\delta t)$. Am I correct? Assume that I am correct at that point, the calculation continues, we obtain:
$$f'=f_{0}+\{f_{0},H\}\delta t+\{\{f,H\},H\}\delta^2 t+\{u_{n}\{f,\phi_{n}\},H\}\delta^2 t+u_{m}\{f_{0},\phi_{m}\}\delta t+u_{m}\{\{f,H\},\phi_{m}\}\delta^2 t+u_{m}\{u_{n}\{f_,\phi_{n}\},\phi_{m}\}\delta^2 t=f_{0}+\{f_{0},H\}\delta t+u_{m}\{f_{0},\phi_{m}\}\delta t $$
I smell something wrong here because after applying second gauge, I can not see where its effect is. There is no appearance of $\phi_{n}$ and $u_{n}$.
Finally, how does it look if we apply $n$ gauge transformations?
For related documents:
http://homepages.uni-r.de/~bon39708/lectures/2017_ws/skript_ws_2017.pdf page 18
http://crypto.fmf.ktu.lt/lt/telekonf/archyvas/KriptoTeorija/2001%20-%20Paul%20Adrien%20Maurice%20Dirac%20-%20Lectures%20on%20Quantum%20Mechanics.pdf page 20-21 of Dirac's lectures.
Update with answer of my question
I had just figured it out how to apply two successive gauge transformation. 
At first, calculate $\Delta f=\epsilon_{m}\{f,\phi_{m}\}$, then first gauge is $f=f_{0}+\Delta f=f_{0}+\epsilon_{m}\{f,\phi_{m}\}$. 
Next, treat $f+\epsilon_{m}\{f,\phi_{m}\}$ as new $f$, we have 
$$f'=f_{0}+\epsilon_{m}\{f_{0},\phi_{m}\}+\delta t\{f+\epsilon_{m}\{f,\phi_{m}\},H\}+u_{n}\{f+\epsilon_{m}\{f,\phi_{m}\},\phi_{n}\}\delta t$$compute $\Delta f'$, we get $\epsilon_{n}\{f+\epsilon_{m}\{f,\phi_{m}\},\phi_{n}\}$. Thus, after appying two successive gauge, what we get is
$$f=f_{0}+\epsilon_{m}\{f,\phi_{m}\}+\epsilon_{n}\{f+\epsilon_{m}\{f,\phi_{m}\},\phi_{n}\}$$
That's how it is done. Still, I wonder why $\epsilon_{m}\phi_{m}$ is generator of gauge.
 A: Ok, let's make a quick review. As you can see from the Dirac's lecture, Hamiltonian is not uniquely determined, but one can add a linear combination of primary constrains on it definition:
$$H_T=H+u_m\phi_m$$
without changing the physics of the system.
This is same as the case in Electrodynamics, where the four-potential $A^\mu$ is not determined uniquely, but:
$$A^\mu\to A'^\mu = A^\mu - \partial^\mu \Lambda$$
This transformation, which leaves physical results unchanged ,is called gauge transformation. Here the additional term gives the gauge freedom.
Returning to our case, the gauge transformation, initially, is the transformation on Hamiltonian as above, with the gauge freedom characterised by the coefficient $u_m$. The realization of this transformation on a phase space function $f$ is:
$$\Delta f(t_0+\delta t)=f_{u_1}(t_0+\delta t)-f_{u_2}(t_0+\delta t) = \epsilon_a[g,\phi_a]$$
This definitely determines the generators of the gauge transformation, which are to be $(\epsilon_a\phi_a)$ characterised by gauge-free parameters $\epsilon_a$. Equivalently, we can present the gauge transformation operator as:
$$U=e^{-\epsilon_a\phi_a}$$
(you can check that by Taylor expanding the gauge transformation $g_2 = Ug_1U^{-1}$ to the first order and realizing that the result is the above equation).
Now, suppose that we start with $g_0$ at some gauge. Now, we apply two successive gauge transformations with generators $(\epsilon_a\phi_a)$ and $(\gamma_b\phi_b)$ to get $g$ and, then, $g'$, respectively:
$$g' = e^{-\gamma_b\phi_b}e^{-\epsilon_a\phi_a}g_0e^{\epsilon_a\phi_a}e^{\gamma_b\phi_b}$$
Now, Taylor expand all the exponents to the first order and expand everything up to the mixed term between $\gamma$ and $\epsilon$ (neglect higher-than-second-oder terms and the second order terms containing only $\gamma$ or only $\epsilon$, also notice that $g_0$, $\phi_a$ and $\phi_b$ do not commute (for $a\neq b$)). And finally, you will get the Dirac's answer.
The multi-gauge transformation is done in the same manner.
I hope that this answers your question.
