I am interested to write down a derivation of Lagrange equations from Newton's second law for a non-holonomic system of particles. Here, I mention my derivation where I am stuck write at the last step.

Consider a system of $N$ particles where their position vectors are written as

$$\mathbf{r}_i=\mathscr{R}_i(q_1(t),\dots,q_M(t),t),\quad i=1,\dots,N\,,\tag{1}$$

where the functions $q_i$ are called the generalized coordinates which are subjected to holonomic and non-holonomic constraints as below

\begin{align*} f_i(q_1(t),\dots,q_M(t),t)&=0,\quad i=1,\dots,C_h\,, \\ g_i(q_1(t),\dots,q_(t),\dot q_1(t),\dots,\dot q_M(t),t)&=0,\quad i=1,\dots,C_n\,, \tag{2} \end{align*}

where $C_h$ and $C_n$ are the number of holonomic and non-holonomic constraints, respectively. Also, if the degree of freedom of the system is $n$ then $n=M-C\ge1$ where $C=C_n+C_h$ is the total number of constraints. Using the chain rule of differentiation we have

\begin{align*} \mathscr{\dot R}_i := \mathbf{v}_i &= \mathbf{v}^*_i+\frac{\partial\mathscr{R}_i}{\partial t},\quad i=1,\dots,N\,, \\ \mathbf{v}^*_i&:=\sum_{j=1}^{M}\frac{\partial \mathscr{R}_i}{\partial q_j}\dot q_j\,,\tag{3} \end{align*}

where we defined the virtual velocity of a particle by $\mathbf{v}^*_i$. Also, from Newton's second law we have

$$\mathbf{F}_i=m \mathbf{a}_i\tag{4},\quad i=1,\dots,N\,.$$

Multiplying both sides of $(4)$ by $\mathbf{v}^*_i$, summing over the number of particles $N$ and interchanging the the order of summations we get

$$\sum_{j=1}^{M}\sum_{i=1}^{N}(\mathbf{F}_i-m\mathbf{a}_i)\cdot\frac{\partial \mathscr{R}_i}{\partial q_j}\dot q_j=0\,.\tag{5}$$

Then using the following definitions and identities

\begin{align*} Q_j&:=\sum_{i=1}^{N}\mathbf{F}_i\cdot\frac{\partial \mathscr{R}_i}{\partial q_j},\quad j=1,\dots,M\,, \\ S_j&:=\sum_{i=1}^{N}m\mathbf{a}_i\cdot\frac{\partial \mathscr{R}_i}{\partial q_j}=\frac{d}{dt}\frac{\partial T}{\partial \dot q_j}-\frac{\partial T}{\partial q_j},\quad j=1,\dots,M\,, \\ T&:=\sum_{i=1}^{N}\frac{1}{2}m\mathbf{v}_i\cdot\mathbf{v}_i\,, \tag{6} \end{align*}

Eq. $(5)$ reduces to

$$\sum_{j=1}^{M}(Q_j-S_j)\dot q_j=0.\tag{7}$$

If there were no constraint equations at all, either holonomic or non-holonomic as mentioned in Eq.$(2)$, then the functions $q_i$ were linearly independent and from this we could conclude that the functions $\dot q_i$ are also linearly independent. Then Eq.$(7)$ would result in the well known form of Lagrange equations $S_j=Q_j$. But here is my question, what if there are constraint equations like Eq.$(2)$. Note that sometimes we are inclined not to eliminate the holomonic constraints by using a transformation. So I am insisting to have both holonomic and non-holonomic constraints at the same time.

As the functions $\dot q_i$ are not (linearly) independent in this case, I am wondering that how the last step works here?

I know that the final result should be

$$\frac{d}{dt}\frac{\partial T}{\partial \dot q_j}-\frac{\partial T}{\partial q_j}=Q_j+\sum_{i=1}^{C_h}\lambda_i\frac{\partial f_i}{\partial q_j}+\sum_{i=1}^{C_n}\mu_i\frac{\partial g_i}{\partial \dot q_j},\quad j=1,\dots,M\,,\tag{8}$$

where $\lambda_i$ and $\mu_i$ are some real numbers which are called Lagrange Multipliers.