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?).