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}^{n}(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}\tag{8}$$

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