## 4 MATLAB Simulation Results  
We verify the performance of Algorithm 1 with MATLAB simulations. The sampling interval is \(dt=1\) second. We use \(Q=36*36\), and consider a given workspace with a known size \([0,200]*[0,200]*[0,200]\) in meters. This implies that the workspace has \([0,200]\) as its x-coordinate range, y-coordinate range, or z-coordinate range. Note that obstacles are unknown a priori. In Algorithm 1, a robot \(n_{i}^{s}\) moves along a directed path to reach its target position. Designing path following controls is not within the scope of this paper. In simulations, the dynamics of \(n_{i}^{s}\) are given by  
\[\mathbf{n}_{i}^{s}(k+1)=\mathbf{n}_{i}^{s}(k)+\mathbf{u}_{i}^{s}(k)dt, \tag{10}\]  
where \(dt\) presents the sampling interval in discrete-time systems. Also, \(\mathbf{u}_{i}^{s}(k)\) presents the velocity of \(n_{i}^{s}\) at sample-step \(k\). The dynamic model in (10) is commonly used in multi-robot systems [10; 33; 34]. Let \(U_{i}^{s}=max_{k}\|\mathbf{u}_{i}^{s}(k)\|\) denote the maximum speed of \(n_{i}^{s}\). Considering heterogeneous robots, \(U_{i}^{s}\
eq U_{j}^{s}\) is possible. The controllers for \(n_{i}^{s}\) are designed in discrete-time systems. Let \(W\) indicate the next waypoint that \(n_{i}^{s}\) will encounter as \(n_{i}^{s}\) travels along the directed path to its target position \(f_{n_{i}^{s}}\). Let \(\mathbf{W}\) indicate the coordinates of \(W\). Let \(\mathbf{g}=\mathbf{W}-\mathbf{n}_{i}^{s}(k)\). \(n_{i}^{s}\) is controlled to move towards \(\mathbf{W}\). If \(\|\mathbf{g}\|>U_{i}^{s}dt\), then (10) is set as  
\[\mathbf{n}_{i}^{s}(k+1)=\mathbf{n}_{i}^{s}(k)+U_{i}^{s}\frac{\mathbf{g}}{\| \mathbf{g}\|}dt. \tag{11}\]  
If \(\|\mathbf{g}\|\leq U_{i}^{s}dt\) and \(W\
eq f_{n_{i}^{s}}\), then \(n_{i}^{s}\) heads towards the next waypoint after \(\mathbf{W}\). If \(W==f_{n_{i}^{s}}\) and \(\|\mathbf{g}\|\leq U_{i}^{s}dt\), then (10) is set as  
\[\mathbf{n}_{i}^{s}(k+1)=\mathbf{n}_{i}^{s}(k)+(\frac{\|\mathbf{g}\|}{dt}) \frac{\mathbf{g}}{\|\mathbf{g}\|}dt. \tag{12}\]  
This leads to  
\[\mathbf{n}_{i}^{s}(k+1)=\mathbf{W}. \tag{13}\]  
We consider two ships in total, i.e., \(S=2\). Using the simulation in the workspace without obstacles (Section 4.1), we can estimate the number of required robots as \(32\). Hence, we conjecture that \(5,10,20\) robots are not sufficient to cover the workspace. We hence use \(50\) robots in the Simulation section. Each ship contains \(25\) robots in total. The position of one ship carrying \(N^{1}=25\) robots is \((30,30,0)\) in meters. In other words, \(n_{1}^{1},n_{2}^{1},\ldots,n_{25}^{1}\) are located at \((30,30,0)\) initially. The position of another ship carrying \(N^{2}=25\) robots is \((150,150,0)\) in meters. In other words, \(n_{1}^{2},n_{2}^{2},\ldots,n_{25}^{2}\) are located at \((150,150,0)\) initially. We consider heterogeneous robots as follows. Each robot may move with distinct speed, while having distinct USBL sensors. In the case where \(i\) mod 3 is zero, \(r_{i}^{dc}\) is 100 m, and the maximum speed is \(U_{i}^{s}=5\) m/s. In the case where \(i\) mod 3 is one, \(r_{i}^{dc}\) is 50 m, and the maximum speed is \(U_{i}^{s}=3\) m/s. In the case where \(i\) mod 3 is two, \(r_{i}^{dc}\) is 75 m, and the maximum speed is \(U_{i}^{s}=2\) m/s. In practice, there exists sonar sensing noise and external disturbance. Considering these practical aspects, a robot may not be located at its target position accurately. As we localize \(n_{i}^{s}\) using the 3D coordinates of \(f_{n_{i}^{s}}\) (Algorithm 1), we added a Gaussian noise with mean 0 and standard deviation 1 m to each element in the 3D coordinates of \(f_{n_{i}^{s}}\). In this way, \(n_{i}^{s}\) is not accurately located at \(f_{n_{i}^{s}}\). Figure 2 shows the final sensor configuration after the robots are deployed to cover the 3D workspace with a known size \(200*200*200\) in meters. 966 s are spent to cover the workspace completely. Using MATLAB simulations, 7 s are spent to build the complete network without sensing holes.