I'm studying from Goldstein's Classical Mechanics, 3rd edition. In section 2.4, he discusses non-holonomic systems. We assume that the constraints can be put in the form
$$f_\alpha(q, \dot{q}, t) =0, \qquad \alpha = 1 \dots m.\tag{2.20}$$
Then it also holds that $$\sum \lambda_\alpha f_\alpha = 0.\tag{2.21}$$ Using Hamilton's principle (i.e. that the action must be stationary), we get that
$$\delta \int_1^2 L\ dt =
\int_1^2 dt\ \sum_{k=1}^n \left(\frac{\partial L}{\partial q_k} – \frac{d}{dt}\frac{\partial L}{\partial \dot{q_k}}\right)\delta q_k = 0\tag{2.22}.
$$
But we can't get Lagrange's equations from this because the $\delta q_k$ aren't independent. However, if we add this with $\sum \lambda_\alpha f_\alpha = 0$, it follows that
$$\delta \int_{t_1}^{t_2} \left(L +\sum_{\alpha=1}^m \lambda_\alpha f_\alpha\right)\ dt = 0.\tag{2.23}$$
And then Goldstein says that
The variation can now be performed with the $n\, \delta q_i$ and $m\, \lambda_\alpha$ for $m+n$ independent variables.
Why have the variables suddenly become independent? First we had $n$ dependent variables, why do we now have $m+n$ independent ones?
Best Answer
Let there be $n$ coordinates $q^j$. The treatment of non-holonomic constraints in Ref. 1 is subpar for various reasons, see e.g. this & this related Phys.SE posts. However we interpret OP's question (v2) as mostly being about counting independent degrees of freedom in constrained systems, and not so much about non-holonomic constraints per se. Therefore, to gain intuition, let us for simplicity just consider $m$ holonomic constraints
$$\tag{A}f_{\alpha}(q)~=~0,$$
where $m\leq n$ (and where we have suppress possible explicit time dependence in the notation). Granted some regularity assumptions, we may in principle solve the $m$ constraints (A) locally so that the coordinates
$$\tag{B}q^j~=~g^j(\xi, \varphi)$$
become functions of $n-m$ independent physical coordinates $\xi^a$, and $m$ coordinates $\varphi^{\alpha}$, in such a way that locally the $n-m$ dimensional constraint surface
$$\tag{C}\{q\in\mathbb{R}^n|f(q)=0\}$$
is parametrized as
$$\tag{D}\{g(\xi, \varphi=0)\in\mathbb{R}^n| \xi\in \mathbb{R}^{n-m}\}.$$
Thus we have at least two equivalent variational formulations:
Reduced formalism: Replace $q^j$ with $g^j(\xi, \varphi=0)$ in the action $S[q]$. Vary the corresponding action $S[\xi]$ wrt. the $n-m$ independent variables $\xi^a$.
Extended formalism: Replace the action $S[q]$ with $$\tag{E} S[q,\lambda]=S[q]+\int\!dt~\lambda^{\alpha}f_{\alpha}(q).$$ Vary the corresponding action $S[\xi,\lambda]$ wrt. the $n+m$ independent variables $q^j$ and $\lambda^{\alpha}$.
The role of the $m$ Lagrange multipliers $\lambda^{\alpha}$ can be view as putting the $m$ variables $\varphi^{\alpha}=0$, so that only the $n-m$ physical variables $\xi^a$ remains, and formulation (2) reduces to (1).
References: