## 3 Formation Tracking Control of Nonholonomic Mobile Robots
### Control Design  
Note that (6) constitutes the new system to be controlled. In (6), \(u_{i1}\), \(\dot{\zeta}_{i}\) and \(u_{i2}\) are the control inputs while \(s_{xi}\), \(s_{yi}\) and \(\dot{\theta}_{i}\) act as the outputs, which has difference with the traditional under-actuated kinematic model in [37]. The transformed system can be treated as three separate single-input single-output systems by applying transverse function technique. Define local error variables as  
\[z_{ix}=\sum_{j=1}^{N}a_{ij}\Big{(}s_{xi}+\varrho_{ix}-\sum_{i=1}^{m_{i}}E_{ix}- s_{xj}-\varrho_{yx}+\sum_{j=1}^{m_{j}}E_{jx}\Big{)}+\mu_{i}\Big{(}s_{xi}+ \varrho_{ix}-s_{sr}-\sum_{i=1}^{m_{i}}E_{ix}\Big{)}\]  
\[z_{iy}=\sum_{j=1}^{N}a_{ij}\Big{(}s_{yi}+\varrho_{iy}-\sum_{i=1}^{m_{i}}E_{iy}- s_{yi}-\varrho_{iy}+\sum_{j=1}^{m_{j}}E_{iy}\Big{)}+\mu_{i}\Big{(}s_{yi}+ \varrho_{iy}-s_{yr}-\sum_{i=1}^{m_{i}}E_{iy}\Big{)}\]  
\[e_{ix}=\tilde{s}_{xi}-\mu_{i}s_{sr}-(1-\mu_{i})\big{(}\rho_{rx}\hat{a}_{rx,i} -\hat{d}_{rx,i}\big{)}+\varrho_{ix}-\sum_{i=1}^{m_{i}}E_{ix}\]  
\[e_{iy}=\tilde{s}_{yi}-\mu_{i}s_{yr}-(1-\mu_{i})\big{(}\rho_{ry}\hat{a}_{ry,i} -\hat{d}_{ry,i}\big{)}+\varrho_{iy}-\sum_{i=1}^{m_{i}}E_{iy}\]  
\[\xi_{i\theta}=\theta_{i}-\theta_{r}\]  
where \(\hat{a}_{rx,i}\), \(\hat{a}_{ry,i}\), \(\hat{d}_{rx,i}\) and \(\hat{d}_{ry,i}\) are used by the \(i\)th robot to estimate unknown trajectory parameters when \(\mu_{i}=0\). The control laws \(u_{i1}\), \(u_{i2}\) and \(\dot{\zeta}\) in transverse function technique are selected as  
\[\left[ \begin{array}{c}u_{i1}\\ \dot{\zeta}_{i}\end{array} \right]=\left[ \begin{array}{cc}1/\tau&0\\ 0&1\end{array} \right]Q_{i}^{-1}\Omega_{i} \tag{13}\] \[u_{i2}=b_{i}/r_{i}\Big{(}-k_{2}\xi_{i\theta}+\frac{\partial h_{ i3}(\zeta_{i})}{\partial\zeta_{i}}\zeta_{i}+\dot{\theta}_{r}\Big{)}\]  
where  
\[\Omega_{i}= -k_{1}p_{i}\bigg{[} \begin{array}{c}z_{ix}\\ z_{iy}\end{array} \bigg{]}-\left[ \begin{array}{c}\dot{\varrho}_{ix}\\ \dot{\varrho}_{iy}\end{array} \right]+\mu_{i}\left[ \begin{array}{c}\dot{\rho}_ {rx}\hat{a}_{rx}\\ \dot{\rho}_{ry}\hat{a}_{ry}\end{array} \right]+(1-\mu_{i})\times\left[ \begin{array} []{c}\dot{\rho}_{rx}\hat{a}_{rx,i}+\hat{d}_{rx,i}\\ \rho_{ry}\hat{a}_{ry,i}+\rho_{ry}\hat{a}_{ry,i}+\hat{d}_{ry,i}\end{array} \right]\]where \(k_{1}\) and \(k_{2}\) are positive constants to be selected. \(P_{i}\) is defined in Lemma 1. Based on the above design, we obtain the following formula:  
\[ \begin{split}\left[\begin{array}{c}\dot{e}_{ix}\\ \dot{e}_{iy}\end{array}\right]&=-k_{1}P_{i}\left[\begin{array}{c}z_{ix}\\ z_{iy}\end{array}\right]\\ \dot{\xi}_{i\theta}&=-k_{2}\xi_{i\theta}\end{split}  \tag{14}\]  
The structure of the parameter estimators in this step is as follows:  
\[ \begin{split}\dot{a}_{rx,i}&=-\gamma_{ri}\rho_{rx}z_{ix}\\ \dot{a}_{ry,i}&=-\gamma_{ri}\rho_{ry}z_{iy}\\ \dot{d}_{rx,i}&=-\gamma_{ri}z_{ix}\\ \dot{d}_{ry,i}&=-\gamma_{ri}z_{iy}\end{split}  \tag{15}\]  
In this step, the Lyapunov candidate function is selected as  
\[ \begin{split} V&=\frac{1}{2}\sum\limits_{i=1}^{N}\Big{(}e _{ix}^{2}+e_{iy}^{2}+\xi_{i\theta}^{2}\Big{)}+\frac{k_{1}}{2}\sum\limits_{i=1 }^{N}\left(1-\mu_{i}\right)\\ &\times\frac{P_{i}}{\gamma_{ri}}(\tilde{a}_{rx,i}^{2}+\tilde{a}_ {ry,i}^{2}+\tilde{d}_{rx,i}^{2}+\tilde{d}_{ry,i}^{2})\end{split}  \tag{16}\]  
where \(\tilde{a}_{rx,i}=a_{rx,i}-\hat{a}_{rx,i}\), \(\tilde{a}_{ry,i}=a_{ry,i}-\hat{a}_{ry,i}\), \(\tilde{d}_{rx,i}=d_{rx,i}-\hat{d}_{rx,i}\), \(\tilde{d}_{ry,i}=d_{ry,i}-\hat{d}_{ry,i}\) are parameter estimation errors. Then, combining Equation (14), the derivative of \(V\) is obtained as follows:  
\[ \begin{split}\dot{V}&=-k_{1}\sum\limits_{i=1}^{N}P_ {i}\big{(}e_{ix}z_{ix}+e_{iy}z_{iy}\big{)}-k_{2}\sum\limits_{i=1}^{N}\xi_{i \theta}^{2}+\sum\limits_{i=1}^{N}\left(1-\mu_{i}\right)\\ &\times\frac{P_{i}}{\gamma_{ri}}(-\tilde{a}_{rx,i}\hat{a}_{rx,i} -\tilde{a}_{ry,i}\hat{a}_{ry,i}-\tilde{d}_{rx,i}\hat{d}_{rx,i}-\tilde{d}_{ry,i}^{2}\hat{d}_{ry,i}^{2})\end{split}