## 3 Proposed Solutions
### Robust Iterative Solution  
The solution above is not optimal, due to the optimization problem (16)-(18) not considering the inner relationship between the last two elements of \(\mathbf{\phi}^{o}\); thus, the result cannot attain the CRLB, as shown in Section 5. Ref. [11] formulated the same optimization problem and solved it using a constrained WLS (CWLS), but the CWLS does not reach the CRLB either. Ref. [11] applied Taylor expansion to formulate a new problem that refines the accuracy. Although the final solution has an appropriate mean square error (MSE) for the CRLB level when the noise is small, it is not sufficiently robust to noise and deviates from the CRLB if the noise is relatively large. Factor graphing (FG) has been applied to the complex localization problem [31; 32; 33]. Once the problem is formulated through factor graphs, finding an effective solution using the sum-product algorithm (SPA) operated by message passing is trivial. FG can provide an iteration solution resilient to noise and initialization [34], which is feasible to the sequential TOA JLAS problem. FG can converge to the optimal, notwithstanding the estimator in Section 3.2 being suboptimal. Factor graphs of sequential TOA JLAS are shown in Figure 2. The rounded rectangles represent the factor nodes, and the circles are variable nodes. Edges connect factor nodes, and variable nodes are the arguments of a local function. Let us start from the maximum likelihood (ML) cost function:  
\[\mathcal{J}(\mathbf{\mu})=\exp\Big{\{}(\mathbf{r}-\mathbf{\bar{r}}(\mathbf{\mu}))^{T} \mathbf{\mathcal{W}}(\mathbf{r}-\mathbf{\bar{r}}(\mathbf{\mu}))\Big{\}}, \tag{40}\]  
where  
\[\mathbf{\mathcal{W}}^{-1}=\mathbf{Q}_{\varepsilon}+\left(\frac{\partial\mathbf{r} ^{o}}{\partial\mathbf{s}^{oT}}\mathbf{Q}_{s}\frac{\partial\mathbf{r}^{T}}{ \partial\mathbf{s}^{o}}\right). \tag{41}\]  
In assuming that \(\mathbf{\mathcal{W}}\) is diagonal, although it is commonly not, the conclusion does not affect this assumption. ML cost (40) can be rewritten as  
\[\mathcal{J}(\mathbf{\mu})=\prod_{i=1}^{M}\exp\Big{\{}w_{i}(r_{i}-\mathbf{\bar{r}} _{i}(\mathbf{\mu}))^{2}\Big{\}}, \tag{42}\]  
where \(w_{i}\) is the \(i\)-_th_ diagonal element of \(\mathbf{\mathcal{W}}\). Denote  
\[h_{i}(\varsigma_{i},\mathbf{\mu}) =\delta(\varsigma_{i}-\mathbf{\bar{r}}_{i}(\mathbf{\mu})), \tag{43}\] \[q_{i}(\varsigma_{i}) =\exp\Big{\{}w_{i}(r_{i}-\varsigma_{i})^{2}\Big{\}}, \tag{44}\]  
where \(\delta(\cdot)\) is the Dirac delta function. Figure 2 indicates that the marginal or belief of \(\mathbf{\mu}\) is the product of all \(M\) incoming messages:  
\[\mathcal{B}(\mathbf{\mu})=\prod_{i=1}^{M}\varpi_{h_{i}\to\mathbf{\mu}}(\mathbf{\mu}). \tag{45}\]  
Since the noise and errors are assumed to be Gaussian, the message from node \(q_{i}\) to \(\varsigma_{i}\) is  
\[\varpi_{q_{i}\to\varsigma_{i}}(\varsigma_{i})\sim\mathcal{N}(\varsigma_{i};r _{i},w_{i}). \tag{46}\]  
The message \(\varpi_{q_{i}\to\varsigma_{i}}(\varsigma_{i})\) is forwarded to the function node \(h_{i}\) by  
\[\varpi_{\varsigma_{i}\to h_{i}}(\varsigma_{i})=\varpi_{q_{i}\to\varsigma_{i}} (\varsigma_{i}), \tag{47}\]which follows a Gaussian distribution.