## 2 Materials and Methods
### The Model of Geodetic Network Accuracy  
The method proposed in this paper is based on the \"greedy algorithm\" [35]. Considering the principle of this algorithm (i.e., maximizing immediate gain; in our case, the accuracy increment), its use could lead to the identification of a local rather than global optimum. Several strategies to overcome this problem (or minimize its effects) have been proposed and will be tested in this paper on four generated 3D geodetic networks. The inputs for the proposed method of optimization include (a) the approximate coordinates of the points that will form the geodetic network, (b) the information on what measurements can be made, and (c) the accuracy characteristics of the device (total station) that will be used for measurements. The algorithm itself then consists, in principle, of four steps:  
1. The measurement (or, rather, a group of measurements) is added into the network if it causes the greatest improvement in the used accuracy characteristic of the (at the time) least accurate point in the network. 2. This is repeated until the required accuracy is achieved for all points. 3. Subsequently, the backward analysis follows, in which a group of measurements is removed from the network if its deletion has the least negative effect on the respective accuracy characteristic of all possible measurement deletions. 4. This is repeated until the removal of the next group of measurements would worsen the accuracy of any point in the network below the required limit. This approach uses the geodetic network adjustment by the least squares method. The accuracy of the network is characterized by a covariance matrix \(M_{\mathbf{x}}\). In the case of a geodetic network with two or more fixed points, the covariance matrix can be calculated as  
\[M_{\mathbf{x}}=c\;(\mathbf{A}^{\mathrm{T}}\mathbf{P}\mathbf{A})^{-1}, \tag{1}\]  
where \(c\) is an arbitrary constant, which has to be identical for all elements in the matrix, \(\mathbf{A}\) is the Jacobi's matrix, i.e., a matrix of partial derivations of the observation equations by network point coordinates (see, e.g., [36] for more details, dimensions are \(m\times u\), where \(m\) is number of measurement in the network and \(u\) is the number of unknowns), and \(\mathbf{P}\) is the weight matrix (dimensions \(m\times m\)). Covariance matrix dimension is \(u\times u\). To be able to use integers as the numbers of measurement repetitions instead of measurement weights, we have modified the definition of the weight matrix that takes them into account as follows:  
\[\mathbf{P}=\mathbf{N}\cdot\mathbf{P}_{0}= \begin{pmatrix}n_{1}&\cdots&0\\ \vdots&\ddots&\vdots\\ 0&\cdots&n_{k}\end{pmatrix} \cdot \begin{pmatrix}p_{1}&\cdots&0\\ \vdots&\ddots&\vdots\\ 0&\cdots&p_{k}\end{pmatrix} , \tag{2}\]  
where \(\mathbf{N}\) is a diagonal matrix containing the numbers of measurement repetitions, and \(\mathbf{P_{0}}\) is a weight matrix with weights for a single repetition of each measurement, both with dimension (\(m\times m\)). Each \(p_{i}\) (\(i=1\) to \(k\)) is then defined as  
\[p_{i}=\frac{c}{\sigma_{i}^{2}}, \tag{3}\]where \(p_{i}\) is the weight of the \(i\)th measurement and \(\sigma_{i}\) is the standard deviation of the \(i\)th measurement. However, if the network is solved as a free network (without a sufficient number of fixed points) additional conditions must be added into the calculation. In this work, conditions for Helmert transformation to all points of the network (\(\mathbf{B}\), with dimensions \(4\times 3l\), where \(l\) is number of points in the network) were added to transform the network so that the sum of squares of displacements (for all points) is minimized as follows:  
\[\mathbf{B}= \begin{pmatrix}Y_{1,0}&X_{1,0}&0&Y_{1,0}\cdots&Y_{l,0}&X_{l,0}\\ 1&0&0&1&\cdots&1&0\\ 0&1&0&0&\cdots&0&1\\ 0&0&1&0&\cdots&0&0\end{pmatrix} , \tag{4}\]  
where \(Y_{j,0}\) and \(X_{j,0}\) are the approximate coordinates of the \(j\)-th point of the network with \(l\) points in total.