How to find symmetry transformations? For a given Lagrangian 
$$ {\cal L} = - \frac{1}{4} F_{\mu \nu} F^{\mu\nu} + |D_{\mu} \phi|^2 -V (\phi) $$
with $\phi = \frac{1}{\sqrt{2}} (\phi^1 + i \phi^2)$, there are the infinitesimal local symmetry transformations 
$$
\begin{align}
\delta \phi^1 &= - \alpha (x) \phi^2  \\
\delta \phi^2 &=  \alpha (x) \phi^1  \\
\delta A_\mu &= - \frac{1}{e} \partial_\mu \alpha
\end{align}
$$
I can plug them in and show that the Lagrangian is invariant. But, where do they come from? Is it just by experience or does the structure of the transformations results from the Lie algebra or other constraints? 
In general the transformations look like they are god given for me even if they reappear in similar way for different types of fields. In other words, why do $\delta \phi^1 $ and $\delta \phi^2 $ mix?
How do I know the structure of the transformation for a field? Is it just trial and error?
 A: The transformations can be deduced essentially from representations of Lie algebras. For the Lagrangian you have written, we are describing a $U(1)$ gauge theory. Irreducible representations of $U(1)$ are one complex dimensional and are always of the form $e^{i q \alpha}$ where $q\in {\mathbb R}$. By definition, this is a representation and therefore acts on a one complex dimensional vector space as
$$
\phi \to e^{i q \alpha} \phi
$$
We can now apply this principle in physics where everything is a function of space-time, $\phi \to \phi(x)$, $\alpha \to \alpha(x)$. You can check that writing $\phi = \frac{1}{\sqrt{2}} ( \phi^1 + i \phi^2)$ reproduces precisely the transformations you have written. 
The transformation of $A_\mu$ can then be deduced from geometric considerations. To do this, let us try and construct a Lagrangian that is invariant under the local gauge transformations
$$
\phi (x) \to e^{i q \alpha(x)} \phi(x)
$$
The first natural thing to do is to try and write down a kinetic term for the scalar. However, the usual one, namely $\partial_\mu \phi \partial^\mu \phi$ no longer works since it transforms weirdly under gauge transformations. The reason of course is that the derivative of a field is defined as
$$
n^\mu \partial_\mu \phi(x) = \lim_{\epsilon \to 0} \frac{ \phi(x+ \epsilon n ) - \phi(x) }{ \epsilon} 
$$
Clearly, the issue is that in the derivative, we are taking a difference of fields at different space-time points! Since the gauge transformation is local, it acts differently on the fields at different space-time points. This is basically the issue. To remedy this, we introduce a different field $W(x,y)$ such that $W(x,x) = 1$ and under gauge transformations transforms as
$$
W(x,y) \to e^{i q \alpha(x) } W(x,y) e^{- i q \alpha(y) } 
$$
We can define a "new derivative" as
$$
n^\mu D_\mu \phi(x) = \lim_{\epsilon \to 0} \frac{W(x,x+\epsilon n) \phi(x+ \epsilon n ) - \phi(x) }{ \epsilon} 
$$
Take the limit described above and you will reproduce precisely the covariant derivative you have in your Lagrangian where
$$
W(x,x+\epsilon n) = 1 - i \epsilon n^\mu A_\mu + {\cal O}(\epsilon^2)
$$
You can also further check that the gauge transformation of $W$ implies the specific gauge transformation of $A_\mu$ that you wrote down. 
Thus, all the gauge transformations can be "derived" from arguments like above. Note also that this discussion can easily be generalized to non-abelian gauge theories.
DISCLAIMER: I HAVE NOT KEPT TRACK OF FACTORS OF $e$ AND SIGNS. PLEASE FIX THEM YOURSELF. 
