classical mechanics - Why can we assume independent variables when using Lagrange multipliers in non-holonomic systems?

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?

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:

1. 
   

   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$.

   
   


2. 
   

   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:

1. H. Goldstein, Classical Mechanics, 3rd ed.; Section 2.4.