A New Network for Particle Filtering of Multivariable Nonlinear Objects †

: In this paper, a new object in the form of a theoretical network is presented, which is useful as a benchmark for particle ﬁltering algorithms designed for multivariable nonlinear systems (potentially linear, nonlinear, and even semi-Markovian jump system). The main goal of the paper is to propose an object that potentially can have similar to the power system grid properties


Introduction
The particle filter (PF) is potentially a very good estimation method because it is based on the optimal solution resulting from Bayesian filtering rules. On the contrary, the biggest disadvantage of PFs is their need for computational power as the number of calculations grows exponentially with the number of system variables [1]. This is the reason why PF methods are usually used only for state estimation with plants of small order.
Some solutions to this problem are hybrid filters, e.g., Rao-Blackwellized PF (RBPF) [2,3] or marginalized PF [4], in which state variables are divided into two groups (one group is estimated using the PF method and the second group using Kalman Filter (KF) methods (linear, extended or unscented)).
Another solution was proposed in [5], where all state variables were divided into groups in which estimation processes were performed; however, the disadvantage of this method was a loss of information contained in measurements, which were based on state variables from two or more groups.
The Dispersed Particle Filter (DPF) was proposed by the authors in one of the previous works [6]; this method assumed that dependencies between state variables and measurements were relatively sparse (transition and measurement functions depended on a relatively small number of different state variables). These research works were performed with the power system grid as an object. Unfortunately, the number of state variables was twice higher than the number of its nodes. Therefore, the performed calculations lasted relatively long, even for small systems.
The Power System State Estimation (PSSE) issue dates back to the 1970s, when Fred Schweppe in his three papers [7][8][9] proposed the PSSE approach using the Weighted Least Squares (WLS) method. This static estimation algorithm is still used in many applications, even after fifty years.
The literature about the PSSE topic is very rich. One can find studies relying not only on the WLS algorithm, but also Kalman filters (extended EKF [10] and unscented UKF [11]), as well as artificial neural networks [12,13]. In [14], the Mixed Kalman Particle Filter (MKPF) method was proposed, which combined elements of the EKF, UKF, and PF algorithms.
Regardless of the estimation method used, a frequent solution is a division of the whole power system grid into subareas [15][16][17][18]. A similar solution was proposed by the authors for the DPF algorithm [6], where the number of subareas was equal to the node number. However, the studies were problematic due to the calculation time. Based on studies from [19], the reasonable number of particles should be about: where N p is the particle number and N x is the number of state variables. The main purpose of the current article is to propose a new network-type system with the number of state variables equal to the number of nodes. Thanks to this, it will be possible to conduct studies even on quite large plants in a relatively short time. The final task will be to propose a particle filter method that will operate satisfactorily (in terms of calculation time and estimation quality) even with large networks. Based on the proposed plant, the subsequent algorithm propositions will be able to be checked quickly, and calculations in power system grids will be performed for results' confirmation only.
However, the proposed network is much more general. It enables fast preparation of even very complex objects. Hence, it can also be used by researchers who are mainly interested in nonlinear estimation issues. Therefore, we provide a very elaborate description of the proposed system in the next section. There are many specification possibilities for both transition and measurement functions. It is also planned to share the source code, which may be used to perform calculations for a given network, without the need for code modification.
In the second section, one can find the description of the particle filter algorithm and the presentation of the proposed network. The simulations performed to evaluate the approach and the results are shown in the third section. The last section contains the conclusions and future plans in this research direction. This paper is a greatly extended version of the article presented during the 13th Federated Conference on Computer Science and Information Systems (FedCSIS 2018) [20].
For a better understanding of the content presented herein, the authors describe the most important symbols in Table 1.
branch measurements between the ith and jth nodes, placed near the ith node q i,(k) normalized weight of the ith particle at the kth time step

Particle Filter
A plant can be expressed by the following state-space equations: where x (k) is a state vector, u (k) is an input vector, z (k) is an output vector, and vectors v (k) and n (k) are internal and measurement noises, respectively; all at the kth time step. The main task of the considered particle filter is to estimate the state vector x (2) based on the measurements and input signals. The particle filters' operation principle is based on recursive Bayesian filtering [21]: where Z (k) is a set of measurement vectors from the first to the kth time step, p(x (k) |Z (k) ) is a posterior Probability Density Function (PDF), p(x (k) |Z (k−1) ) is a prior PDF, p(z (k) |x (k) ) is a likelihood, and p(z (k) |Z (k−1) ) is evidence (normalization factor).
The key idea in the considered PF is to implement the posterior PDF as a set of particles. The ith particle is represented by a pair x i,(k) , q i,(k) composed of the state value and the weight. The higher the weight, the higher the probability that the value x i,(k) is close to the real state vector. If the number of particles, N p , is high enough, the information about the posterior PDF contained in the particle set is the same as in the continuous PDF function. The first particle filter was proposed in 1993 by Gordon, Salmond, and Smith [22] and was called the bootstrap filter. The operation principle of this PF is presented in Algorithm 1.
Algorithm 1: bootstrap filter 1. Initialization: draw N p initial values x i,(0) from initial PDF p(x (0) ); set time step k = 1; 2. Prediction: draw N p new particles from the transition model x i,(k) ∼ p(x (k) |x i,(k−1) ); 3. Update: compute the particles weights based on the measurement modelq i,(k) = p(z (k) |x i,(k) ); 4. Normalization: normalize weights so that their sum is equal to 1: , whereq (k) are unnormalized weights, calculated in Step 3. 5. Resampling: draw N p new particles using the posterior PDF obtained in previous steps (the chance that the particle will be drawn is equal to the normalized weight); 6. End of iteration: calculate the estimated state vector; increase the time step k = k + 1; go to Step 2.
The considered PF can be used for both non-Gaussian distributions and complex transition models. An example of such a complex system model can be a semi-Markovian jump system [23], in which the whole system model can switch to a different structure (with some probability) described by other equations and also can go back to the previous "system equations".

Resampling
Resampling is a very important step of the algorithm. Its main task is to prevent particles from degeneration. The Sequential Importance Sampling (SIS) is an example of the algorithm without resampling and is a direct predecessor of the PF. The above-mentioned problem relates to the effect when all but one particles have their weights equal to or near zero (and the weight of that one particle tends to 1).
The particle filter is inherently connected to resampling; hence, the first usage of resampling was presented in the same article in which the particle filter was proposed [22]. During the resampling stage, N p particles must be drawn from the PDF obtained in previous steps. At every draw stage, the probability that the specific particle will be drawn is equal to its normalized weight.
The first resampling algorithm that was proposed in the literature is an algorithm called later "multinomial resampling", which requires N p random numbers u i from uniform PDF U (0, 1). The jth particle x j,(k) is duplicated (one can also say that it is reproduced) if the condition: holds. However, this approach has a computational complexity of order at least O(N p log N p ). The linear complexity provided later resulted in resampling being proposed, e.g., "systematic resampling", which was used by the authors in all experiments presented here. In this resampling, only one random value is needed, and each of the other "drawn" values is higher by 1 N p from the previous one. The systematic resampling is presented in Algorithm 2 (it is assumed that rand() means a random value drawn from uniform distribution U (0, 1)).
u=u+(1/N_p); 11. end This particular resampling, as well as many other schemes can be found in [24], in which over 20 different algorithms, which can be used in PFs, were presented and compared.
However, even a linearly-executed resampling is the "bottleneck" of the whole PF algorithm, because all other steps can be performed in parallel. This is the reason why in some solutions, the resampling is not executed in every iteration, but only if the Effective Sample Size (ESS) (marked by N e f f ) is lower than a specified threshold (usually taken as a half of the particle number [25,26]). This parameter can be calculated by: between Steps 4 and 5 of Algorithm 1 (after weight normalization). However, based on [27], (5) is not the only estimator of ESS, and quite good properties are also given by the measure 1 max{q (k) } , because the range of the practical threshold influence is wider. The ESS parameter can be treated as the information of how many particles give significant information about the current PDF. It can take any value in the range between 1 and N p . For example, if all the particles have weights equal to 0. 25, then: For the case with particle weights q 1,(k) = 0.1, q 2,(k) = 0.2, q 3,(k) = 0.3, q 4,(k) = 0.4, the ESS is equal to 3.33, and if only one particle has a non-zero weight, then N e f f = 1.

Proposed Network
The authors proposed a grid, which on the one hand could be easily adapted for particular requirements, and on the other hand provided wide possibilities in the new plants creation. This explained why there are many complex and maybe even incomprehensible options presented below; however, the most common networks will be presented in a very simple way. Objects created in this way will not make any physical sense, but will be an appropriate and demanding benchmark for any proposed algorithm. Moreover, if one needs to add any dependence, which is not specified, there still is a possibility to present this on the scheme.
The proposed network is composed of nodes and lines (branches). Two different nodes can be connected by a line. Exactly one state variable x i is associated with the ith node (i = 1, . . . , N x , where N x is the number of state variables). Nodes can be represented in two ways: by circles or by squares. Transition function f i depends on the shape of the ith node.
In the first transition function type (circles on the scheme), a state variable depends only on itself (value from the previous time step). In the second type (squares), a state value depends on all state variables, which are connected (by branches) to the calculated state variable. Expressions for transition functions of the first type (circles) are presented below, and their connections with scheme designations are described in Table 2.
It should be noted that more than one of these expressions can be used for a specific node. For example, one can obtain transition function (11) (wide equation) by the graphical notation presented in Figure 1. In such cases, an internal noise v (k) i is of course added only once. For the second type of transition functions, connections between nodes play a vital role. The branch between the ith and jth nodes has a value µ i,j = µ j,i and order m i,j = m j,i ; whereas, if there is no connection between the ith and jth nodes, the branch has a value µ i,j = µ j,i = 0 (line order m i,j does not matter in such a case). It is also assumed that µ i,i = 1 and m i,i = 1.

Equation
Designation of Node Explanation (7) α and n parameters should be given in the scheme.
However, if one of these parameters is omitted, it is assumed that this value is equal to 1.
Furthermore, both values can be omitted; in such a case, α = 1 and n = 1.
(8) and (9) When all three values are given, one must take into account that α should be written before β (above or on the left side of β).
If any parameter (α, β, or n) is not presented, it is assumed that it is equal to 1; however, one should keep in mind that β can be omitted only if α is also omitted.
The designations for Equation (9) are the same, but triple lines (through circles) should be used.
Parameter ϕ should be always presented in degrees ( • ) to be able to distinguish it from other values. Omitted ϕ means that it is equal to zero. (10) p J should be written outside of the circle as J's sub-or super-script (also from the left side).
If p J is omitted, it is assumed that p J = 0.5.
-To use another function (which must be explained in the text), a double circle should be used.
The types of lines are specified, which differ in line functions f l i,j (where i is the number of nodes, from which line function "was called"). These functions are defined in the remaining part of the paper, and their connections with scheme designations are described in Table 3. Line functions f l i,j will be used in both transition and measurement functions. If anywhere, f l i,i occurs (e.g., in Equation (17)), it is assumed that the first type of line function should be used, and all line parameters are the default. Table 3. Explanation of designations: lines.

Equation Line Designation Explanation
Values m i,j and µ i,j should be written on the scheme near the line center. Furthermore, information about the branch type, in the form of diagonal lines (one for (12), two for (13), and three for (14)) should be presented there. If one (or both) of the parameters µ i,j and m i,j = 1 is omitted, it is assumed that its value is equal to 1.
If branches in a whole grid are only of the first type, diagonal lines can be omitted.
-To use another line function, one should draw a double line between nodes; however, this line function should be explained in the text.
The second type of transition functions (squares), which depend also on other state variables, are presented below, i.e., in Equations (15)- (18), and their designations are presented and described in Table 4.

Equation Designation of Node Explanation
The n value should be written inside a square bracket. If the α or n parameter is omitted, it is assumed that the omitted parameter is equal to 1. (16) and (17) When all three values are given, one must take into account that α should be written before β (above or on the left side of β).
If any parameter (α, β, or n) is not presented on the scheme, it is assumed that it is equal to 1; however, one should keep in mind that β can be omitted only if α is also omitted.
Parameter ϕ should always be presented in degrees ( • ) to be able to distinguish it from other values. Omitted ϕ means that it is equal to zero. (18) p J should be written outside of the square as J's sub-or super-script (can be also from the left side). Omitted p J means that this value is equal to 0.5.
-To use another transition function (which must be explained in the text), a double square should be used.
It is also possible to combine both transition function types (circles and squares). An example is presented in Figure 2, and the appropriate transition function of the first state variable is written in Equation (19). To obtain the jump function (Equation (18)) for both types of transition function, the "J" mark and p J value should be given in the place where the square and circle are connected.  (19).
Inputs should be presented by the element shown in Figure 3. The specific input signal should be presented inside the element, and as long as it is possible, all needed parameters for simulation should be described inside the input figure (see Figure 3). 3. Exemplary input signal notation; the input signal should be presented as an arrow-shaped pentagon. As long as it is possible, the input signal should be integrally described inside the shape. If it is not possible, the pentagon should be doubled, and an explanation should be given in the text.
One state variable can have more than one input signal. If it is not obvious that the signal is periodic, one should insert information about the period length (in time steps), e.g., "T = 20". If all information about the input signal cannot be presented on the scheme, a double "input figure" should be used (examples are presented in Figure 3 on the right side).
The type of internal noise also should be presented on the scheme, or optionally, it should be described in the text. The PDF type with parameters should be connected to a specific node by a dashed line. Examples are presented in Figure 4. Measurements are presented on the scheme by filled marks. There are two types of measurements: nodal (associated with a specific node) and branch (associated with a specific line). The nodal measurement can depend on a state variable of this node and also on all state variables of the nodes that are connected with this node by lines. The branch measurement depends only on the two state variables of the two nodes that are connected by the line where the measurement is placed. In the case of branch measurements, one should keep in mind that the first index informs about near to which node a specific branch measurement is placed (branch measurements are not symmetric in general).
All nodal measurements proposed in the new network (P i , Q i , R i , S i , T i ) are described by (20)- (24), and examples of the designations on the scheme are presented in Table 5; whereas all proposed branch measurements are written by (25)- (29), and the exemplary designations are presented in Table 6.

Measurement designation Explanation
Value λ should be always below or on the right side of the α parameter.
If α or λ are skipped, one can assume that they are equal to 1. However, if γ is omitted, one should assume that γ = 0. Parameter α can be skipped only if γ is also omitted.
Nodal measurements should be placed near the node or should be connected with the node by the dashed line. In both cases, one can add information about measurement noise or the measurement scaling parameter (α). Table 6. Explanation of designations: measurements, Part 2; for branch measurements between the ith and jth nodes, but at the ith node.

Equation Measurement Designation Explanation
Branch measurements should be placed on the specific line. One should keep in mind that measurements are different at both ends of the branch, so the position of the mark should be unambiguous. For information about measurement noise or its scale parameter, a dashed line should be used.
The measurement noise should be marked in the same way as the internal noise: by a description in the text or using a dashed line on the scheme. Examples of measurement noise designations are presented in Tables 5 and 6.
To use another measurement function, one can simply use a new filled mark on the scheme (such a function should be explained in the text).

Application for Power Systems
As was mentioned in the Introduction, the network described above can be created like the power system grid. Mainly, the S-type measurements were proposed based on the power flows and power injections, which occur in power systems. Below, 4 measurements from power grids with a description are presented; however, one should take into account that the notation presented in these equations absolutely cannot be associated with any markings of the rest of the paper (the notation is typical for the power system subject). The first two are power injections (active and reactive, respectively): and the two subsequent are power flows (active and reactive, respectively; in power flows, the first index informs about near to which bus (node) the measurement is placed): In the above equations, B is the number of buses in a system. Y i,j , µ i,j and y i,j are system parameters.
The state vector in power system grids is composed of voltage magnitudes U and phase angles δ. One phasor (voltage magnitude and phase angle) is associated with every bus (node). Hence, it is visible that the length of the state vector is two times higher than the number of nodes (actually, the number of state variables is given by the 2B − 1 expression, because one node must be a reference node, in which the phase angle is constant, usually taken as 0).
The system composed only of S-type measurements (nodal and branch) and P-type measurements (nodal only) should successfully simulate the power system grid properties. Therefore, one can use a system created in such a way to speed up studies (due to the lower number of state variables), and after finding a promising method or algorithm, it can be checked on power systems.
The reader is directed to [32,33] for more information about power systems.

Exemplary Objects from the Literature
There are a few plants, which were used in literature, presented below. The first object is usually used in research associated with particle filters. It was proposed in 1978 by Netto, Gimeno, and Mendes (according to [34]). The object is a one-variable plant, the state space model of which is described by Equations (34)- (36). This model can be presented using the designations proposed in this paper, which one can see in Figure 5. The second plant, described by Equations (37)- (39), was proposed in [35]. It can also be presented using the proposed designations, which is shown in Figure 6. This object was created based on the state Equation (34), but both (transition and measurement) functions were simplified. Figure 6. Scheme of the plant, which was presented in [35] and described by Equations (37)-(39).
The subsequent object was used in [36] and is described by Equations (40)- (42). The description, using the proposed designations, is presented in Figure 7. As one can see, parameter ϕ = −90 • was used on the scheme to obtain "sin" instead of "cos" from Equation (17). The plant described by Equations (43)-(45) differs from other objects in such a way that its measurement function is switched after 30 time steps. It is marked in Figure 8 by two measurements, but with additional information, for the range of time steps for which it should be used. Moreover, the gamma PDF was used as internal noise. The last object, which is mentioned below, was used in [38]. It is a tracking system, and its state-space model is given by Equation (46). The first two variables (x 1 , x 2 ) are coordinates; the other two (x 3 , x 4 ) are velocities; and the system scheme is presented in Figure 9. In the referenced article, the authors did not provide information about internal and measurement noises (it is only known that they were Gaussian). State variables x 3 and x 4 do not have internal noise, and this is why on the scheme, one can find "PDF" U (0; 0). Figure 9. Scheme of the object, which was presented in [38] and described by Equation (46). As was presented above, it is possible to model systems presented in the literature using the network proposed in this paper; however, one should keep in mind that the main goal was the easy creation of complex, multivariable, nonlinear systems, and not modeling the existing ones. Thanks to this, research on the plants with a network-type structure is possible with a relatively short state vector (number of state variables equal to the number of nodes).

Plants Used in the Experiments
The experiments were performed for the networks presented in Figures 10-12. On all these schemes, there are no measurements; this is because they were considered in 5 different cases: P, Q, R, S, and T. In each case, the measurements were located in every possible place, but without repetitions (10 branch measurements and 4 nodal measurements). One can see an example in Figure 13: it is a S-type case for the Ob402 object.    Finally, one should obtain 15 different systems; however, objects Ob402 and Ob403 had the same transition functions, and the difference was only in line functions f l i,j and in line parameters, which impacted the measurements; however, P-type measurements were identical for Ob402 and Ob403. Hence, the case with P-type measurements for Ob402 was omitted.
The transition function of these three systems is presented below: Ob401 is described by (47) and (48), and both objects, Ob402 and Ob403 are described by (49) and (50). All the measurement functions (for any case) one may find in the appendices: for system Ob401 in Appendix A, for system Ob402 in Appendix B (but the equations are the same as for Ob401) and for Ob403 plant in Appendix C (all measurements had default parameters:

Results
The experiments were performed to check the dependence of the estimation quality on ESS and also on the complexity level of the transition and measurement functions.
The estimation quality was evaluated using the average root mean squared error, which is defined as: where N x is the number of state variables, and: where x with a hat is the estimated state variable, x with a plus is the real value of state variable, and M is the simulation length (set on the 1000 time steps in all experiments). The aRMSE index can be used only for comparison; it does not make any physical sense. The lower the value of this index, the better the estimation quality. Every experiment was different, values of the initial state vector (at zero time step) were drawn from uniform distribution U (−1; 1). It was also assumed that initial values of particles were the same as the initial state vector.
All simulations were repeated until the standard deviation of the aRMSE mean value was less than 0.0025 (the variance of the mean value was s times smaller than the variance of s samples [39]), but at least 100 times. It was not achieved only for a few cases, which are mentioned in Table 7. Moreover, also for these cases, there were ranges on the charts, which showed, with 95% probability, where the mean value was placed (according to the 68-95-99.7 rule [40]). The results are presented in Figures 14-17.

Discussion and Further Work
As one can suppose, the lowest aRMSE values were obtained for the completely linear case Ob403-P (except the highest number of particles, where the Ob403-R case was the best). However, in general, the complexity level of equations (measurement or transition functions) had no impact on the estimation quality, or this influence was insignificant. In most of the cases, one could see the dependence: the greater the effective sample size, the better the estimation quality.
One could see that the results obtained for objects Ob402 and Ob403 in cases with Q-type and S-type measurements were almost the same. This was due to the fact that these systems were very similar (identical transition functions and measurement functions different only in system parameter µ i,j ). In R-type and T-type measurement functions, the line functions f l i,j were used, which were very different between these two plants.
Based on the ESS results in Figures 16 and 17, one can conclude that the number of particles N p impacted the ESS (higher N e f f with higher N p ), but did not impact the order of cases (graphs did not intersect). It was also visible that ESS was relatively high for P-type measurements, which could have an impact on the good estimation quality.
The most interesting was probably the fact that the dependence "the more non-linear transition and measurement functions, the worse the estimation quality" was not always met. This was visible primarily for S-type measurements (see Figure 15): for the Ob401 system, the estimation quality was definitely better than for the two others. Similarly, in cases with T-type and R-type measurements (see Figure 14), a better estimation quality was obtained for the Ob401 system than for the Ob402 plant. Through all the experiments, the internal noise was the same. The following question arises: What was the reason for such results?
For almost all cases, one could see the relationship between the effective sample size and the estimation quality: the higher the N e f f , the lower the aRMSE. The Ob401-Q case was an exception: very high ESS (in comparison to other cases) resulted from the specificity of measurements (products of two or more state variables) and from the fact that for two state variables, "jumps" to the opposite values were possible (see Equation (47)). Consequently, although there was a relatively high number of significant particles in PF, a part of them was incorrect (with opposite signs).
However, the above-mentioned relationship was not satisfied for T-type and R-type measurements. Why was there the smallest ESS number for Ob403 systems, but simultaneously the best estimation quality was obtained (in comparison to the two other objects with the same measurements type)?
In these cases, the influence had a range of measurements, which one can see in Table 8. Variances of the subsequent measurements (from exemplary experiments) are presented in this table. A higher variance meant that, during the simulations, the values of that measurement changed more. One should keep in mind that the measurement noises were the same (had the same distribution) during the experiments for all cases. The higher values the measurement had, the lower the impact the measurement noise had on the estimation process. As a result, the small N e f f values in these cases were associated with more accurate measurements, and not with difficulties in state vector tracking.
This could be presented in a one-dimensional example with measurement function z (k) = ax (k) + n (k) (n (k) ∼ N (0, 1)). For smaller values of a, a higher value of ESS was obtained, whereas for a higher value of a, lower ESS, but in the latter case, particles with significant weights would be definitely closer to the sought value. An example is presented in Figure 18.  N (0, 1)). For a = 1, measured value was assumed to be one, and for a = 20, the measured value was assumed to be 20 (for better readability). (a) Exemplary probability density functions; (b) exemplary particles (with normalized weights) drawn from PDFs. However, regardless of the considered case, one could notice a very low number of significant particles; in most cases from two to 30. How precisely will the four-dimensional PDF be modeled using even 30 points? Using 81 points, one would have accuracy similar to modeling a one-dimensional PDF by three points.
Therefore, the authors considered that further studies should be directed at ESS increasing, and the improvement of estimation quality would be only the result of these actions. They would be started with a comparison of different PF methods and resampling algorithms, in terms of state estimation quality and effective sample size in multivariable systems.
One possible solution may be the usage of an adaptive sample size [41]; however, this would only result in the increase of the particle number, as much as ESS would achieve a set threshold. Another method, which should improve the effective sample size, is the Population Monte Carlo (PMC) algorithm [42,43]. In this approach, it is assumed that particles are divided into a few groups and each one is drawn from different PDFs. Every PDF is calculated based on the drawn, up to this moment, particles and their weights. Yet another solution could be a double drawing: the first time locating an interesting area and the second time drawing mainly from this subspace. The authors are already working on a new PF algorithm that will be based on the latter approach.
The network proposed in this article had of course some disadvantages and limitations: the graphical representation of more complex systems with a higher number of nodes would be illegible (a description in the text would also be required). Despite the fact that many different options were proposed, the one was strongly limited in the creation of transition functions, and measurements were mostly non-parametric (only scale factor α was available). Furthermore, only additive noises were considered (could be also multiplicative and others).
Despite all these disadvantages, the authors believe that the proposed network was sufficient for a variety of particle filter studies, and not only this one. Additionally, for the main application (simulation of power grid behavior), it was totally suitable.
In the future, the authors will prepare and share the source code, which will allow performing experiments, not only for the networks presented in the article, but also for any other one (based on the external file with a system structure description). Acknowledgments: Thanks to Talar Sadalla and Szymon Drgas for their help in the preparation of the conference paper, which was the basis for the creation of the current paper. We also thank Dariusz Horla for his help at the very end.

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

Abbreviations
The following abbreviations are used in this manuscript: · 0.5 ln 0.001 + x