ON THE NUMERICAL MINIMISATION OF THE OBJECTIVE FUNCTION APPLIED TO SPHERICAL HARMONICS FITTING

The paper presents some considerations on the performance of various objective function minimization methods in the process of GNSS antenna PCV determination. It is particulary important in the case of structural health monitoring and diagnostics. PCV are used as an additional feature to improve the GNSS positioning accuracy. The process of PCV derivation is complex and involves fitting spherical harmonics into a set of observables. The paper compares computing performance and accuracy of few methods used in the fitting process.


INTRODUCTION
GNSS antenna calibration is a process in which phase centre variations (PCV) are estimated for a given GNSS receiver antenna.It is important in high-accuracy applications of GNSS like coordinate systems definition or deformation measurements.This aspect is particularly important in the diagnostics of the structural health of constructions where the deformations are often at the level of single millimetres.The PCVs obtained from the calibration process are disseminated as "individual" or "type" ANTEX files.These files are then used in GNSS data processing software.To derive the PCVs, a few steps must be completed: (1) measurements, (2) data preprocessing, and (3) deriving a model.This article is focused on the last part of the process.Antenna PCV models are delivered to end users as tabular data containing PCV values in the regular grid on the hemisphere.These values are used to interpolate corrections to observed phase measurements in receiver-satellite line of sight (LOS) in GNSS data postprocessing.The observables used to estimate the PCV are uniformly distributed on the hemisphere.To estimate the correction values in ANTEX files, spherical harmonics are used.
Spherical harmonics are special functions defined on the surface of a sphere.Their main advantage is that each function defined on the surface of a sphere can be written as a sum of spherical harmonics.
If we assume that the square-integrable 2-D spatial function f(x) = f(λ, ϕ) is specified on a sphere ΩR, it can be defined as a series expansion in terms of spherical harmonics Yn,m (λ, ϕ), which is computable using relations from the associated Legendre function Pn,m(sin ϕ) of order m and degree n.Values cn,m are the coefficients of the spherical harmonics.In numerical applications, truncating the expansion at finite degree n = N is needed, such that the double sum would be composed of all N′ = (N + 1) 2 words (11).
In the PCV model estimation, spherical harmonics are fitted into a set of observation data with the least squares (LS) method.The minimum of the LS objective function, which is: (3) must be calculated in the process (where v is a residuum).The minimum can be found using various approaches to objective function minimization.In this article we have tested following methods: SLSQP, CG, Powell, BFGS, Newton-CG, L-BFGS-B, TNC, Cobyla, trust-constr.LS fitting can be done by solving a set of normal equations defined as: This approach is often not optimal in the scope of computational efficiency.Therefore in this paper, we compare other, numerical methods of finding the objective function minimum.

METHODS
There are many methods for finding the minimum of the objective function.In order to find a more optimal algorithm, a total of nine methods were selected and then tested.Below is a brief description of each algorithm: 1. SLSQP: is a gradient-based method for solving nonlinear optimization problems with constraints.The problem is solved iteratively for a local minimum (6). 2. CG: it is based on selecting orthogonal search directions n and then, in each direction, performing one step, the size of which should correspond to the proposed solution x in that direction.The CG method is useful for large and sparse problems since the number of iterations in this method is not greater than the number n (13).3. Powell: the algorithm uses a two-way search.It first moves in one direction until it finds a local minimum.Then it performs the same in the other direction.This process is repeated for subsequent directions until the fit statistic for a given iteration is minimized.The algorithm terminates when no significant improvement is achieved in subsequent iterations (4).4. BFGS: The descent direction is determined by preconditioning the gradient using curvature information.This is done by successively improving the approximation of the loss function's hessian matrix, which is extracted from the gradient evaluations or their approximations using the generalized secant method (1). 5. Newton-CG: uses conjugate gradient (CG) method to the second-order Taylor-series approximation of f around the current iterate xk (10).6. L-BFGS-B: this method is an extension of the L-BFGS algorithm, which aims to minimize the differentiable scalar function f(x) over unconstrained values of the real vector x.The L-BFGS-B method, unlike the L-BFGS method, can handle bounds imposed on the variables (14).7. TNC: this method to determine an update to the function's parameters approximately solves Newton's equations by repeatedly applying the optimization algorithm.The inner solver runs for only a limited number of iterations.This implies that the inner solver must return a good approximation in a finite number of iterations for the method to work (7) (8).8. Cobyla: is based on approximating an actual constrained optimization problem using linear programming problems.A candidate for the optimal solution is obtained by solving the problem in a single iteration.The acquired candidate is then evaluated using objective and constraint functions, resulting in a new data point in the optimization space.Then, the approximation of the linear programming problem used in the next iteration is improved using the acquired information.If the solution cannot be improved, the step size is reduced, thus improving the search.The cycle repeats until the step size is small enough (9).9. trust-constr: depending on the problem definition, the method uses one of two implementations: trust-region sequential quadratic programming solver or trust-region interior point method.The first is used for problems with inequality constraints.The second is when inequality constraints are imposed (15).

Data acquisition
The research aims to determine the most optimal method for calculating the parameters of a spherical harmonic in terms of time and correctness of the result.The research scope includes 9 minimization algorithms: SLSQP, CG, Powell, BFGS, Newton-CG, L-BFGS-B, TNC, Cobyla, trust-constr.The most optimal algorithm is selected by testing the methods on a single set of measurement data acquired with a GNSS antenna during actual field measurements and then comparing the tested quantities, i.e., calculation time and root mean square (RMS) value.The field measurements were collected by authors using the GNSS antenna calibration facility developed at the University of Warmia and Mazury in Olsztyn, Poland.Czyrzniak M, Rapiński J.: On the numerical minimisation of the objective function applied to spherical …

Compared values
To compare the efficacy of objective function minimization methods, two quantities were used: the method processing time understood as a time required by a computer to obtain a correct solution, and the root mean square (RMS) calculated on the basis of residuals.
Due to processes taking place in computer processors, individual time measurements of a given method may vary slightly.In order to obtain a reliable result, the execution time of each method was measured 10 times and then the average value was calculated.Time measurement was started exactly before the method was performed and ended immediately after.Below is a fragment of code for time measurement written in Python.start_time = datetime.now()x=minimize(residuals_differential, x0,args=(m,l,dane), method=method, tol=0.001,options={'maxiter':10000}) end_time = datetime.now() The first and last lines of the code recorded the time before the algorithm was run and after it was completed.Based on the difference of these values, the execution time of the method was determined.Between the time measurements, the algorithm of the method is run.The parameters of the function vary depending on the method used.

Input data
In order to obtain a correct PCV model, the acquired observations should be spread over the entire hemisphere of the antenna.In order to acquire such data, a precision robot was used, which has the ability to rotate and tilt the antenna.The antenna took measurements at 1-second intervals.The resulting observations include zenith angle and azimuth for the antenna in two epochs t and t+1, along with time differences of a single difference (TDSD), which can be defined as the difference of two consecutive single differences.Given a phase observation:  (5).Since the distance between two antennas is very small (up to few meters) the ionospheric and tropospheric effects are cancelled.To remove the receiver clock error, both receivers are connected to the common frequency source (the rubidium oscillator in this case).Hence TDSD is defined as: The geometric term    can be calculated on the basis of known receiver positions and precise orbits.On this basis the set of normal equations are created and the PVC model can be estimated.

Processing
The compared methods were tested in the same environment to avoid the impact of changes in technology on calculation time measurements and for the same parameters n and m (degree and order expansion), which were set to the value of 5.However, it is difficult to determine the correctness of the methods without knowing the result they should return.For the comparison of methods to make sense, it is necessary to assume some valid values for the polynomial parameters.For this purpose, the polynomial parameters obtained by the method of least squares were taken as correct values.These values, burdened by randomly generated numbers from 0 to 1, were then used as initial parameters for estimation using tested methods.
Each method takes a tolerance parameter, which is responsible for the precision of the results and affects the number of iterations performed by the algorithm.This parameter was defined individually for each method in order to obtain enough iterations to observe the stabilization of the value of the sum of residuals difference calculated at the end of each iteration using the callback function.The process of adjusting the tolerance parameter was done through Convergence charts analysis.Czyrzniak M, Rapiński J.: On the numerical minimisation of the objective function applied to spherical …

DISCUSSION
The results of the analyses were the average execution time obtained from 10 samples and RMS values for each method (Table 2).A Convergence graph was also plotted for each method.

Time and RMS values
The fastest method was the Newton-CG, which performed calculations in just 0.01 s.The second fastest algorithm was the Cobyla method, for which the average calculation time was 5.71 s.In contrast, the longest calculation was performed by the Powell method 113.37 s.The RMS values are very similar for all methods and vary by 1 of ten-thousandths of a millimeter for most of them.Only the RMS value for the Newton-CG method deviates from the rest and has the highest value.Therefore, it can be observed that the Newton-CG method is the fastest but also the least accurate.

Convergence graphs
Convergence graphs for each method except the Newton-CG method are included below.This method, with the assumed parameters, performed in one iteration therefore the graph was abandoned.

Fig. 1. Convergence of SLSQP method
The algorithm completed calculations in 16 iterations.In the first iteration, the value of sum of relative residuals reached 26.2672 mm.From the 10th iteration onward, the value began to circulate close to the value of 14 mm.In the last iteration, it amounted to 14.13141 mm.

Fig. 2. Convergence of CG method
The algorithm completed calculations in 53 iterations.In the first iteration, the value of sum of relative residuals reached 14.1254 mm.After that it dropped drastically and from around the 25th iteration onward, the value began to circulate close to the value of 14.04 mm.In the last iteration, it amounted to 14.0332 mm.

Fig. 3. Convergence of Powell method
The algorithm completed calculations in 13 iterations.In the first iteration, the value of sum of relative residuals reached 13.9232 mm.From 7th iteration onward, the value began to circulate close to the value of 14.04 mm.In the last iteration, it amounted to 14.0496 mm.

Table 2 .
Time and RMS values