Estimation of Performance Parameters of Turbine Engine Components Using Experimental Data in Parametric Uncertainty Conditions

: Zero-dimensional models based on the description of the thermo-gas-dynamic process are widely used in the design of engines and their control and diagnostic systems. The models are subjected to an identiﬁcation procedure to bring their outputs as close as possible to experimental data and assess engine health. This paper aims to improve the stability of engine model identiﬁcation when the number of measured parameters is small, and their measurement error is not negligible. The proposed method for the estimation of engine components’ parameters, based on multi-criteria identiﬁcation, provides stable estimations and their conﬁdence intervals within known measurement errors. A priori information about the engine, its parameters and performance is used directly in the regularized identiﬁcation procedure. The mathematical basis for this approach is the fuzzy sets theory. Synthesis of objective functions and subsequent scalar convolutions of these functions are used to estimate gas-path components’ parameters. A comparison with traditional methods showed that the main advantage of the proposed approach is the high stability of estimation in the parametric uncertainty conditions. Regularization reduces scattering, excludes incorrect solutions that do not correspond to a priori assumptions and also helps to implement the gas path analysis with a limited number of measured parameters. The method can be used for matching thermodynamic models to experimental data, gas path analysis and adapting dynamic models to the needs of the engine control system.


Introduction
Aircraft gas turbine engines are the most expensive and vital components of the aircraft. Their maintenance and operating costs depend on the engine performance, which is mainly determined by a condition of the engine gas path components such as the compressor, turbine, combustor, ducts, etc. For this reason, engine maintenance is usually supported by a prognostics and health management system (PHM), which usually implements a diagnostic procedure known as gas path analysis (GPA).
Nowadays, different approaches are used for processing flight and test-cell data to estimate component performance and detect faults [1]. These methods can be divided into two groups [2]: • Model-based methods which use thermodynamic models to relate failures to measured parameters of the gas path [3]

Basic Identification Procedure
The model calculates parameters of the engine gas path Y  at steady-state modes depending on an operating point, external conditions U  and component performance   . Hence, in the general case, it is presented as: The linear model may be formulated as which relates small deviations of the gas path parameters Y   and parameters of components' performances   at a single operating point. H is an influence coefficient matrix (ICM). The linear model (2) where is the Euclidean norm of the vector, Φ-least squares functional to be minimized.
The precision of estimation will improve if the amount of a priori information is higher than the number of unknown parameters. Therefore, the identification can be provided by a set of measurements that are done at N different operating points. For this purpose, the generalized vector of residuals is minimized: The proposed approach is a significant modification of the standard model-based procedure for estimation of the performance parameters of engine components. The ability to implement prior knowledge in different forms to stabilize the solution is the main contribution of the paper. Regularized multi-criteria identification enables engine users to implement Non-Linear GPA [1,2] in ill-conditioned configurations with a limited number of measured parameters.

Basic Identification Procedure
The model calculates parameters of the engine gas path → Y at steady-state modes depending on an operating point, external conditions → U and component performance → Θ. Hence, in the general case, it is presented as: The linear model may be formulated as which relates small deviations of the gas path parameters δ where is the Euclidean norm of the vector, Φ-least squares functional to be minimized. The precision of estimation will improve if the amount of a priori information is higher than the number of unknown parameters. Therefore, the identification can be provided by a set of measurements that are done at N different operating points. For this purpose, the generalized vector of residuals is minimized: Model (1) which is included in (4) is numerical. Therefore, the minimization of residuals (4) is a numerical iterative procedure, which in each iteration determines the solution as a sum of the previous solution and current correction: (5) and iterations continue until corrections become negligible.
The correction ∆ → Θ n+1 is determined as a solution of the overdetermined system of linear algebraic equations: is a generalized matrix of influence coefficients, which is composed of elementary matrixes H i ( → Θ n ) determined for each operating point. In accordance with Equation (6), for each iteration, such correction ∆ → Θ n+1 of required parameters is found, which minimizes residuals The known solution of the linear system (6) by the least squares method takes the form: where A = C T GC-Fisher information matrix, G-weight diagonal matrix, whose elements are inverse to variations of measuring errors σ 2 y i . Unfortunately, the least squares method is sensitive to outliers in the right part of the system (6), which can be related to faults in experimental data. Therefore, practical applications of the classic LSM suffer from the instability of estimations ∆ → Θ n+1 and poor convergence of the identification algorithm as a whole. In some cases, the estimations→ Θ are far from the expected values. The reasons for such results may be the lack of empirical information, an excessive number of estimated parameters or correlation between two or more state parameters. This causes ill-conditioning of the Fisher Matrix and excessive deviations of estimations. These estimations lose a physical sense (are out of range of possible components' performances parameters variation-for example, efficiency is more than one) and can cause the model calculation program to crash.
Hence, the LSM identification procedure needs modification to provide its stability and physically adequate estimations.

Regularized Identification Procedure
In the known regularization procedure named for Tikhonov, the generalized functional is minimized, which besides the residuals of measured parameters includes a norm of the finding parameter vector with a weighting factor α. The identification task takes the following form: where g i is the weight coefficient that describes measurement error (g i = 1 The regularized identification is based on the least squares method described above. As it is shown in Equation (9), the minimized functional is extended by elements which are related to parameters → Θ to be found. Thus, the main changes in the identification procedure are the generalized vector of residuals → Z( → Θ) and generalized matrix of influence coefficients C( → Θ). The modified terms have the following structure: where I is an identity matrix. This identification procedure will deliver the non-regularized solution if the regularization coefficient is zero. Otherwise, the influence of coefficient α depends on the proportion between components: and r k=1 θ 2 k of the generalized functional Φ (Θ). Therefore, the effect of regularization is determined by the conditions of identification such as the number of measured parameters m, number of operating points N, measuring errors G and number of estimated parameters r.

Numerical Simulation
To assess the impact of the regularization coefficient α, the identification procedure was tested with the zero-dimensional thermodynamic model of a turboshaft engine [17] (Figure 2), based on experimental component maps. Performed numerical experiments involved calculation of N = 12 operating points for two scenarios: (1) Without measuring noise, for a wide range of α values; (2) With measuring noise, for selected α values.
where gi is the weight coefficient that describes measurement error ( i The regularized identification is based on the least squares method described above. As it is shown in Equation (9), the minimized functional is extended by elements which are related to parameters   to be found. Thus, the main changes in the identification procedure are the generalized vector of residuals Z ( ) and generalized matrix of influence coefficients C ( )   . The modified terms have the following structure: where I is an identity matrix.
This identification procedure will deliver the non-regularized solution if the regularization coefficient is zero. Otherwise, the influence of coefficient α depends on the proportion between components: the effect of regularization is determined by the conditions of identification such as the number of measured parameters m, number of operating points N, measuring errors G and number of estimated parameters r.

Numerical Simulation
To assess the impact of the regularization coefficient α, the identification procedure was tested with the zero-dimensional thermodynamic model of a turboshaft engine [17] (Figure 2), based on experimental component maps. Performed numerical experiments involved calculation of N = 12 operating points for two scenarios: (1) Without measuring noise, for a wide range of α values; (2) With measuring noise, for selected α values. Five gas path parameters i Y were used for model identification: (1) Compressor discharge pressure p2; (2) Turbine inlet temperature TIT; (3) Power turbine discharge temperature; (4) Rotational speeds of rotors N2 and N1.
Initial deviations of two compressor performance parameters were assumed as follows:  Five gas path parameters Y i were used for model identification: (1) Compressor discharge pressure p 2 ; (2) Turbine inlet temperature TIT; (3) Power turbine discharge temperature; (4) Rotational speeds of rotors N 2 and N 1 .
Four estimated parameters included:  for different values of α, which varied in a range from 0 to 1. ∆Y av is the average deviation between initial and estimated models: and ∆Y * av is the average deviation between measurements and estimated model: which varied in a range from 0 to 1. ΔYav is the average deviation between initial and estimated models: and * av Y  is the average deviation between measurements and estimated model: Figure 3. Influence of the regularization coefficient on estimations when noise is absent. Figure 3 shows that despite the initial growth of 3  and 4  absolute values, the average value a v  decreases monotonically, whereas the residual between measurements and the corrected model increases. This corresponds to theoretical predictions of the influence of regularization on the identification process. The simulation results help to choose the value of the regularization coefficient α. If the deviation of parameters Y  is to be lower than 5 % and the deviation of parameters  lower than 1%, then  Figure 3 shows that despite the initial growth of Θ 3 and Θ 4 absolute values, the average value Θ av decreases monotonically, whereas the residual between measurements and the corrected model increases. This corresponds to theoretical predictions of the influence of regularization on the identification process.
The simulation results help to choose the value of the regularization coefficient α. If the deviation of parameters → Y is to be lower than 5 % and the deviation of parameters Θ lower than 1%, then the value of α cannot exceed 0.03.
To obtain precise values of variance and to analyse the distributions of estimations→ Θ, a sequence of data sets with random measuring errors with variance σ 2 Y i was generated and identified 1000 times. This provided the average precision of about 1% of initial residuals for parameters Y and about 0.5 % of simulated deviation 0.03 for the estimated parameters→ Θ. Figure 4 presents histograms for calculations with two estimated parameters at α = 0 and α = 0.04. Both non-regularized and regularized estimations have distributions close to normal ones. Total scatter of estimations is also similar, but centres of regularized estimations moved from their true values by 0.008 and 0.005, respectively.    At low values of the regularization coefficient α, centres of distributions are close to their true values, but scatter is high, and distributions are strongly asymmetrical. When α increases, then mean values move, but scatter decreases, and its distribution approaches the normal one.
As shown above, the value of the regularization coefficient α has to be set appropriately to provide meaningful results of regularized least squares. This procedure is useful in ill-conditioned    At low values of the regularization coefficient α, centres of distributions are close to their true values, but scatter is high, and distributions are strongly asymmetrical. When α increases, then mean values move, but scatter decreases, and its distribution approaches the normal one.
As shown above, the value of the regularization coefficient α has to be set appropriately to provide meaningful results of regularized least squares. This procedure is useful in ill-conditioned engine model fitting but needs preliminary adjustment and sensitivity analysis. At low values of the regularization coefficient α, centres of distributions are close to their true values, but scatter is high, and distributions are strongly asymmetrical. When α increases, then mean values move, but scatter decreases, and its distribution approaches the normal one.
As shown above, the value of the regularization coefficient α has to be set appropriately to provide meaningful results of regularized least squares. This procedure is useful in ill-conditioned engine model fitting but needs preliminary adjustment and sensitivity analysis.

Regularized Identification Procedure Development Using a Priori Information
This paper demonstrates an idea to implement engine heuristics involving component parameters and performances to improve the precision and stability of model identification. The main difficulty is in the diversity of this information, which is presented in one of the following forms: • Exact statement (for example, a part-load performance in a determined area is smooth); • Statement in the form of limitations of the area of acceptable solutions (for example, the efficiency of the individual compressor cannot differ more than 3% from the efficiency of an "average" compressor, the performance of which is used in the initial model); • Statement in the form of fuzzy information (for example, the gas temperature in the turbine will grow with the engine's life); • Statistical form (for example, probability density functions of parameters).
The next difficulty is in the formalization of parameters that characterize the model quality. These parameters are set on the basis of subjective preferences of decision makers (DM). The same difficulties appear in the ranking of partial criteria and limitations according to their significance for estimation of the model quality. The analysis showed that the main problem is that theoretical-probabilistic methods are difficult to apply in the presence of uncertainties, which are related to subjective preferences whose nature is not statistical. Actually, the choice of the model's structure is the DM procedure, which in multi-criteria case is inevitably highly subjective. If the complexity of the task increases, the role of quality factors will grow. Therefore, it is possible to take into account all criteria using the proper mathematical tool.
We propose to use the fuzzy set theory [33] as this tool. This theory makes the uniform base to describe the information given in all the forms listed above, thus providing the correct mathematical definition of the identification task.
An example of applying fuzzy sets in the GPA and engine model matching is given by M. Zwingenberg et al. [34]. They used fuzzy logic for the evaluation of sensor failures. In contrast, we introduce a priori information about engine performance and experimental data directly into the stabilizing functional through the fuzzy logic approach.
Generalization of the functional (9) gives: where Φ e ( → Θ) is the least squares functional (3), which minimizes measurement error, so it is called the empirical risk functional and Φ a ( → Θ) is the stabilizing functional, which is considered as the functional of a priori risk.
The determined a priori information is set as a limitation, which may take the form of an equality or inequality. The more general form of setting a priori information is its representation as a fuzzy set. The fuzzy set of parameter x is represented as a definition domain and a membership function µ(x) in this domain. For the considered task, the parameters x may be the model parameters → Θ or the output (calculated) variables → X. Next, we will use limited normal fuzzy sets with limited definition domain (x min , x max ) and 0 ≤ µ(x) ≤ 1. We will express each a priori information as a particular functional of a priori risk Φ a q ( → θ ) and will determine the general functional of a priori risk as a linear composition of particular functionals: where α q are weight coefficients. Let us consider some types of a priori information and its expression in view of fuzzy sets: Case 1. Limitations of some parameters are known. For example, it is known that 0 < η < 1 and 0 < σ < 1, etc. Using experience and calculation results, these limits can be significantly reduced; for example, 0.5 < η < 0.9 and 0.9 < σ < 0.99 ( Figure 6).

Case 2.
The a priori mathematical model is known. This can be the model with design maps of components or the model of the average engine, which is matched with previous testing results. This information may be expressed as a relationship between μa or Φa and the difference between parameters that correspond to the matched and a priori models. These parameters may be the model parameters Θ  as measured (for example fuel flow) or non-measured calculated parameters (for example thrust). The membership function, in this case, can be of symmetrical triangular shape and the functional of a priori risk a q ( )    that characterizes the similarity between values of the parameter q x of matched and a priori models can be formed as N 2 a q q q j q j j 1 By analogy, the stabilizing functional can be formed as a second power of a distance from estimated parameter k  (or function x( Θ  ) that is determined using the estimated parameters).
The main problem in this analogy is that the fuzzy set that contains a priori information is infinite, so the sum in the above equation must be replaced with integral. For example, the second power of the distance from the engine map parameter θ = x to the fuzzy set with a given membership This information may be expressed as a relationship between µ a or Φ a and the difference between parameters that correspond to the matched and a priori models. These parameters may be the model parameters → Θ as measured (for example fuel flow) or non-measured calculated parameters (for example thrust). The membership function, in this case, can be of symmetrical triangular shape and the functional of a priori risk Φ a q ( → θ ) that characterizes the similarity between values of the parameter x q of matched and a priori models can be formed as By analogy, the stabilizing functional can be formed as a second power of a distance from estimated parameter θ k (or function x( → Θ) that is determined using the estimated parameters).
The main problem in this analogy is that the fuzzy set that contains a priori information is infinite, so the sum in the above equation must be replaced with integral. For example, the second power of the distance from the engine map parameter θ = x to the fuzzy set with a given membership function µ(θ) equals: In the iteration process of the task solution, each integral must be calculated numerically. This requires a lot of time and high computational capacity. A numerical solution can be found for certain types of the membership function. We considered the trapezoidal function, the particular cases of which are the rectangle, triangle and isosceles trapezium (Figure 7).
In the iteration process of the task solution, each integral must be calculated numerically. This requires a lot of time and high computational capacity. A numerical solution can be found for certain types of the membership function. We considered the trapezoidal function, the particular cases of which are the rectangle, triangle and isosceles trapezium (Figure 7). The functional derived for the trapezium: 12a 12b 4c x 6a 6b 12ab 4ac 4bc c } / ( b ). 12 2 2 (18) The proposed modifications of the objective functional cause the non-linearity of the estimation problem, so traditional solution methods are not applicable. This problem is overcome by adapting genetic optimization algorithm to the specifics of the engine model matching [35]. Genetic algorithms are increasingly used in gas turbines to solve complex optimization problems [4,10,36,37].

Results
The designed methods were implemented using experimental data obtained during the test-cell testing of the helicopter turboshaft engine. Table 1 contains the values of measured parameters. The last row presents mean squared measurement errors σi, which were used to form the diagonal weight matrix G (Equation (8)). The model parameters to be corrected:  The functional derived for the trapezium: The proposed modifications of the objective functional cause the non-linearity of the estimation problem, so traditional solution methods are not applicable. This problem is overcome by adapting genetic optimization algorithm to the specifics of the engine model matching [35]. Genetic algorithms are increasingly used in gas turbines to solve complex optimization problems [4,10,36,37].

Results
The designed methods were implemented using experimental data obtained during the test-cell testing of the helicopter turboshaft engine. Table 1 contains the values of measured parameters. The last row presents mean squared measurement errors σ i , which were used to form the diagonal weight matrix G (Equation (8)). The model parameters to be corrected: π C * (c)-Scaling factor of CPR (Compressor Pressure Ratio); π C * (a)-Factor of rotation in CPR-W Ccor plane of the compressor pressure map; π C * (b)-Factor of rotation in CPR-N cor plane; W HPT (c)-HPT flow coefficient; W PT (c)-Power turbine flow coefficient; η C * (a)-Factor of rotation in the η C * -W Ccor plane of the compressor efficiency map; η C * (c)-Scaling factor of the compressor efficiency map; η CC -Scaling factor of combustion efficiency;

Least Squares Identification
The initial values of model parameters → Θ 0 that correspond to the average engine (a priori information) are shown in Table 2. Table 2. Initial values of the model parameters.
Such large corrections make it impossible to calculate part-load performance and influence coefficient matrices for the next iteration. Therefore, the calculation procedure was modified: Corrections to the component maps were limited. At corrections as high as 7%, the procedure slowly becomes convergent (by 7-9 iterations). The solution corresponds to small (less than 1%) corrections to the component maps.
The results of the standard least squares (the Newton-Raphson method) show that this method has low stability due to a weak conditionality of the engine model. Therefore, practical application requires improvement of the method's stability.

Regularized Identification
The same experimental data were processed by the regularized method based on Equation (9).
Parameter α was set at 0.01. Table 3 shows corrections δ → Θ s for the first three iterations while Table 4 presents the parameters calculated by the corrected model.  Table 3 shows the high stability of the procedure: Estimations change slowly, without significant oscillations. The changes become smaller from iteration to iteration (this is a sign of stability). Tables 1 and 4 shows that parameters determined by the corrected model are well correlated with the measured parameters. Thus, the developed procedure of component performance estimation and the engine model matching can be implemented in practice.

Regularized Multi-Criteria Identification
Finally, the above task of multi-mode diagnostics was considered with a priori information about the model of the average engine, the parameters of which were estimated by previous testing data ( Table 2). It is also known that a scatter of measured parameters of different engines at P net = const, N 1 = const follows the normal distribution with the following mean squared error (Table 5). Table 5. Mean squared error of performance parameters in the engine population. Parameter Mean squared error 0.85% 0.7% 1.1% 0.6% This information was then transformed into the functional Φ a (Equation (13)), and taking into account that in this case, the functionals Φ e and Φ a are homogeneous as they contain residuals by parameters of the same names. This facilitated the setting of weight coefficients in Equation (14). They had to relate to a scatter of measurements and a scatter of the engine parameters by series. This provided the composition of the functionals: In the considered example  Table 6 shows the estimated parameters of the model in subsequent iterations. The final results of the corrected model are presented in Table 7. The comparison of Tables 1, 4 and 7 confirms that the proposed procedure is stable. For the corrected model, deviations of output parameters from experimental data are within the measurement error. The introduction of a priori information results in a smaller deviation of results compared to the conventional regularization method.  Figure 8 presents the experimental data and modelling results for part-load performance. As the initial model is far enough from the field data, the diagnostic algorithm needs a preliminary adjustment to make the model closer to the average engine and to improve the precision of GPA. The visible lines correspond to the initial model before the correction, the average engine model and corrected model of the analyzed engine. The deviations of performance parameters from the reference indicate the engine's health status.

Discussion
Matching turbine engine models to experimental data is an inverse problem of mathematical modelling which is often characterized by parametric uncertainty. This results from the fact that the number of measured parameters is significantly lower than the number of components' performance parameters needed to describe the real engine. It was shown that in such conditions, even small measurement errors resulted in a high variation of results, hence the obtained efficiency, loss factors, etc., exceeded physically feasible ranges.
It was confirmed that extending the least squares functional by regularizing term improved stability of estimation but caused significant bias, which could lead to incorrect diagnoses. Stable operation of the algorithm and tolerant bias of the estimates in wide range of engine operation conditions are difficult to achieve as only selected values of the regularization coefficient α can ensure that. Finding its optimal solution is complicated because the regularizing term is not related to the physical properties of the object.  Figure 9 illustrates the correction of compressor maps and confirms that the procedure was effective even though the initial model was inaccurate.

Discussion
Matching turbine engine models to experimental data is an inverse problem of mathematical modelling which is often characterized by parametric uncertainty. This results from the fact that the number of measured parameters is significantly lower than the number of components' performance parameters needed to describe the real engine. It was shown that in such conditions, even small measurement errors resulted in a high variation of results, hence the obtained efficiency, loss factors, etc., exceeded physically feasible ranges.
It was confirmed that extending the least squares functional by regularizing term improved stability of estimation but caused significant bias, which could lead to incorrect diagnoses. Stable operation of the algorithm and tolerant bias of the estimates in wide range of engine operation conditions are difficult to achieve as only selected values of the regularization coefficient α can ensure that. Finding its optimal solution is complicated because the regularizing term is not related to the physical properties of the object.

Discussion
Matching turbine engine models to experimental data is an inverse problem of mathematical modelling which is often characterized by parametric uncertainty. This results from the fact that the number of measured parameters is significantly lower than the number of components' performance parameters needed to describe the real engine. It was shown that in such conditions, even small measurement errors resulted in a high variation of results, hence the obtained efficiency, loss factors, etc., exceeded physically feasible ranges.
It was confirmed that extending the least squares functional by regularizing term improved stability of estimation but caused significant bias, which could lead to incorrect diagnoses. Stable operation of the algorithm and tolerant bias of the estimates in wide range of engine operation conditions are difficult to achieve as only selected values of the regularization coefficient α can ensure that. Finding its optimal solution is complicated because the regularizing term is not related to the physical properties of the object.
The new regularization procedure of multi-criteria identification, introduced in Section 2.4, uses a priori information about the engine, its measurement system, mathematical model and expected performance (Figure 1). This information is heterogeneous: Partially, it is presented in the form of statistic parameters and partially in the form of heuristic expert representations. Some other methods, e.g., Extended Kalman Filter can also use a priori information [22]. Undoubtedly, the quality of implemented expert knowledge should be carefully checked in advance.
The framework that supported the formalization of the regularization problem used fuzzy sets theory [33,37]. Similarly to the solution of Tanaka et al. [38], membership functions were introduced directly into the minimized functional, which contained a priori information about an acceptable range of estimated parameters. However, fuzzy logic is not used here for the non-linear mapping of parameters like in AI-based diagnostics [39]. In this work, inputs and outputs are measured variables with Gaussian noise.
The integral form of the stabilization functional was proposed which operates with continuous membership functions. To accelerate numerical calculations, the trapezoidal membership function was recommended. The analytical equation was derived for the stabilization functional with this function, which was used in place of numerical integration, thus decreasing the computing capacity required for the estimation procedure.
In the example presented in Section 3.3, available component maps were significantly different from the actual maps of the tested engine. Nevertheless, the estimation algorithm worked stable and provided a perfect matching of the model to the experimental data. Thanks to regularizing terms implemented in the estimation procedure, intermediate estimations did not exceed their limits when the computational algorithm was unstable.
Further studies and more field data are needed to quantify the accuracy and the stability of the procedure in relation to the number of estimated parameters. However, when a large number of engine parameters is measured, measurement uncertainty is low, and the number of coefficients of the model to be found is significantly less than the number of measured parameters, the estimates will not exceed the acceptable range, stability is not a problem, and the regularization does not have a positive effect and is not required. Hence, the proposed method is efficient when the number of measured parameters is small, and errors of measurement are significant. This is often the case in real applications.

Conclusions
A new method for the stable estimation of engine performance parameters is proposed. It uses a priori information about the engine and data acquisition system, i.e., the expected performance and acceptable variations of parameters. The method improves the accuracy of gas-turbine performance models in off-design conditions. A priori information is introduced in the model identification procedure in the form of membership functions, which are transformed into additional stabilization terms of the functional to be minimized.
The proposed approach was used to process test-cell data from a helicopter engine. The comparison with the conventional least squares and Tikhonov regularization proved the following: (1) LSM has low stability that can cause the engine simulation software to crash.
(2) Conventional regularization improves the stability, but the intensive bias of the estimates needs to be excluded in advance.
(3) A priori information has a physical sense, so it improves the stability and precision of the estimation. The introduced identification method is effective in ill-conditioned configurations and provides small bias of estimates. Acknowledgments: The authors would like to acknowledge the constructive comments, and suggestions provided by the anonymous reviewers that greatly improved the quality of the article. We are also grateful to Michail Ugryumov and Anatolyi Mazurkov for their help in the application of the genetic optimization procedure and Kacper Grzędziński for his comments on an earlier version of the manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.

Nomenclature
The