Why do we insert the Lagrangian with some Lagrange multipliers into the Euler-Lagrange-equation? Let's assume that we have two coordinates $q_1$ and $q_2$, a Lagrangian function $\mathcal{L}(q,\dot{q},t):=\mathcal{L}(q_1, q_2, \dot{q_1}, \dot{q_2},t)$ and a constraint $f(q,t) = 0$. Then I understand the procedure to obtain the equations of motions as constructing a new Lagrangian function $\mathcal{L}'(q,\dot{q},t) = \mathcal{L}(q,\dot q, t) + \lambda(t) f(q,t)$ where $\lambda$ is the Lagrange multiplier. Then you insert this in the Euler-Lagrange-equations (where you treat $\lambda$ like a coordinate $q_3$, so you get three equations) and these equations of motions then describe your constrained system.
I don't really understand why this is working. In regular non-physical optimization problems I already constructed Lagrangian functions with the multipliers analogous to the above one, but then I always took the gradient of the Lagrangian function and set it equal to zero. I also understand why this works, because this procedure results in a calculation of points where the gradient of the function of interest and the gradient of the constraint function are parallel, and this is a necessary condition for a maximum which also satisfies the constraint. But in physics we don't do the gradient of the Lagrangian - we just insert it in the Euler-Lagrange-equation (or is this somehow equivalent here? If yes, why?).
 A: So we are given a lagrangian of the form $\mathcal{L}(q_1,q_2,\dot{q_1},\dot{q_2},t)$ with constraint equation $f(q_1,q_2,t) = 0$.
Well, according to Hamilton's principle we have
$$\delta\int_{t_1}^{t_2}\mathcal{L}dt=\int_{t_1}^{t_2}\left( \left[ \frac{\partial\mathcal{L}}{\partial q_1}-\frac{d}{dt}\left(\frac{\partial\mathcal{L}}{\partial \dot{q_1}}\right )\right]\delta q_1 +\left[ \frac{\partial\mathcal{L}}{\partial q_2}-\frac{d}{dt}\left(\frac{\partial\mathcal{L}}{\partial \dot{q_2}}\right )\right]\delta q_2 \right) dt=0 $$
At given instant $\delta q_1$ and $\delta q_2$ are not independent but
$$\delta f= \frac{\partial f}{\partial q_1}\delta q_1+\frac{\partial f}{\partial q_2}\delta q_2=0$$
Thus we have
$$\int_{t_1}^{t_2}\left( \left[ \frac{\partial\mathcal{L}}{\partial q_1}-\frac{d}{dt}\left(\frac{\partial\mathcal{L}}{\partial \dot{q_1}}\right )\right] \frac{1}{\frac{\partial f}{\partial q_1}} -\left[ \frac{\partial\mathcal{L}}{\partial q_2}-\frac{d}{dt}\left(\frac{\partial\mathcal{L}}{\partial \dot{q_2}}\right )\right]\frac{1}{\frac{\partial f}{\partial q_2}} \right)\delta q_1 dt=0 $$
which is true for all choices of $\delta q_1$
$$\frac{\frac{\partial\mathcal{L}}{\partial q_1}-\frac{d}{dt}\left(\frac{\partial\mathcal{L}}{\partial \dot{q_1}}\right )}{\frac{\partial f}{\partial q_1}}=\frac{\frac{\partial\mathcal{L}}{\partial q_2}-\frac{d}{dt}\left(\frac{\partial\mathcal{L}}{\partial \dot{q_2}}\right )}{\frac{\partial f}{\partial q_2}}$$
One obvious way in which we can solve this equation is to separately set both sides equal to the same function of time, which we shall denote $-\lambda(t)$. It follows that the Lagrangian equations of motion of the system can be written
$\displaystyle \frac{d}{dt}\!\left(\frac{\partial L}{\partial \dot{q}_1}\right)-\frac{\partial L}{\partial q_1} - \lambda(t)\,\frac{\partial f}{\partial q_1}$ $\textstyle =$ $\displaystyle 0,$
$\displaystyle \frac{d}{dt}\!\left(\frac{\partial L}{\partial \dot{q}_2}\right)-\frac{\partial L}{\partial q_2} - \lambda(t)\,\frac{\partial f}{\partial q_2}$ $\textstyle =$ $\displaystyle 0.$
which is same as to solve
$$\mathcal{L}′(q,\dot{q},t)=\mathcal{L}(q,\dot{q},t)+\lambda(t)f(q,t)$$
