## 3 Identification Algorithm  
If we have \(m\) observations from \(t=1\) to \(t=m\), Equation (4) can be written in a stacked form  
\[\mathbf{Y}=\Phi\mathbf{\theta}+\mathbf{V}, \tag{8}\]  
where  
\[Y :=[y(1),y(2),\cdots,y(m)]^{T}\in\mathbb{R}^{m},\] (9) \[\Phi :=[\mathbf{\varphi}(1),\mathbf{\varphi}(2),\cdots,\mathbf{\varphi}(m)]^{T} \in\mathbb{R}^{m\times n},\] (10 From Equations (2), (5) and (10), we can see that the information matrix \(\Phi\) contains many unknown intermediate terms. Therefore, it is difficult to perform the identification directly. According to the auxiliary model identification idea [24; 27], the information matrix \(\Phi\) can be replaced with its estimate  
\[\hat{\mathbf{\Phi}}_{k}:=[\hat{\mathbf{\varphi}}_{k}(1),\hat{\mathbf{\varphi}}_{k}(2), \cdots,\hat{\mathbf{\varphi}}_{k}(m)]^{T}, \tag{11}\]  
where  
\[\hat{\mathbf{\varphi}}_{k}(t) :=[\hat{\mathbf{\varphi}}_{1,k}^{T}(t),\hat{\mathbf{\varphi}}_{2,k}^{T}(t ),\cdots,\hat{\mathbf{\varphi}}_{r,k}^{T}(t)]^{T}, \tag{12}\] \[\hat{\mathbf{\varphi}}_{i,k}(t) :=[-\hat{x}_{i,k-1}(t-1),\cdots,-\hat{x}_{i,k-1}(t-n_{ai}),u_{i}(t -1),\cdots,u_{i}(t-l)]^{T}. \tag{13}\]  
Note that the unmeasurable terms \(x_{i}(t-j)\) are replaced with their auxiliary model output estimates \(\hat{x}_{i}(t-j)\). Then, the parameter vector \(\theta\) can be estimated by the auxiliary model based least squares iterative (AM-LSI) algorithm [28],  
\[\hat{\mathbf{\theta}}_{k} =[\hat{\mathbf{\Phi}}_{k}^{T}\hat{\mathbf{\Phi}}_{k}]^{-1}\hat{\mathbf{\Phi} }_{k}^{T}\mathbf{Y},\ k=1,2,3,\cdots \tag{14}\] \[\hat{x}_{1,k}(t) =\hat{\mathbf{\varphi}}_{1,k}^{T}(t)\hat{\mathbf{\theta}}_{k}(1:n_{a1}+l),\] (15) \[\hat{x}_{i,k}(t) =\hat{\mathbf{\varphi}}_{i,k}^{T}(t)\hat{\mathbf{\theta}}_{k}\left(1+ \sum_{j=1}^{i-1}(n_{aj}+l):\sum_{j=1}^{i}(n_{aj}+l)\right),\ i=2,3,\cdots,r, \tag{16}\]  
where \(\hat{\mathbf{\theta}}_{k}\) denotes the parameter vector estimate at the \(k\)th iteration. According to the least squares (LS) theory, the AM-LSI algorithm is efficient if it is satisfied that \(m\gg n\). However, from Equation (7), we can see that the dimension of the system in Equation (8) is high. Therefore, it would take a lot of time and efforts to obtain enough observations to meet the identification requirement. Moreover, Equation (14) shows that the AM-LSI algorithm requires computing the inverse matrix \([\hat{\mathbf{\Phi}}_{k}^{T}\hat{\mathbf{\Phi}}_{k}]^{-1}\) at each iteration, which leads to a heavy computational burden. Furthermore, the sparse solution cannot be obtained [23], and the time-delays cannot be effectively estimated. Thus, the AM-LSI algorithm is infeasible for the high dimensional and sparse system identification. Inspired by the compressive sensing theory, the identification of the sparse system in Equation (8) can be further expressed as an optimization problem  
\[\hat{\mathbf{\theta}}=\arg\min\|\mathbf{\theta}\|_{0},\ \ \text{s.t.}\ \ \|\mathbf{Y}-\Phi\mathbf{ \theta}\|_{2}^{2}\leq\varepsilon,\]  
where \(\|\cdot\|_{0}\) represents the \(l_{0}\) norm, \(\|\cdot\|_{2}\) the \(l_{2}\) norm, and \(\varepsilon\) the error tolerance, which is a priori chosen. However, it is a non-convex problem and is difficult to solve in practice. A commonly used alternative is to replace the \(l_{0}\) norm with a relaxed convex \(l_{1}\) norm [9],  
\[Y :=[y(1),y(2),\cdots,y(m)]^{T}\in\mathbb{R}^{m},\]0  
where \(\|\cdot\|_{1}\) is the \(l_{1}\) norm. It is proved that the minimization of \(l_{1}\) norm is equivalent to the minimization of \(l_{0}\) norm when the Restricted Isometry Property (RIP) is satisfied [13]. Then, the convex optimization problem in Equation (17) is solvable with the BPDN criterion [16]. The BPDN approach can effectively reduce the interference of noise and has good robustness. Inspired by the auxiliary model and the iterative identification idea in the AM-LSI algorithm, we apply the BPDN criterion to modify the AM-LSI algorithm by replacing the LS step with the BPDN criterion.