Non-Linear Response of Cable-Mass-Spring System in High-Rise Buildings under Stochastic Seismic Excitation

In high-rise buildings earthquake ground motions induce bending deformation of the host structure. Large dynamic displacements at the top of the building can be observed which in turn lead to the excitation of the cables/ropes within lift installations. In this paper, the stochastic dynamics of a cable with a spring-damper and a mass system deployed in a tall cantilever structure under earthquake excitation is considered. The non-linear system is developed to describe lateral displacements of a vertical cable with a concentrated mass attached at its lower end. The system is moving slowly in the vertical direction. The horizontal displacements of the main mass are constrained by a spring-viscous damping element. The earthquake ground motions are modelled as a filtered Gaussian white noise stochastic process. The equivalent linearization technique is then used to replace the original non-linear system with a linear one with the coefficients determined by utilising the minimization of the mean-square error between both systems. Mean values, variances and covariances of particular random state variables have been obtained by using the numerical calculation. The received results were compared with the deterministic response of the system to the harmonic process and were verified against results obtained by Monte Carlo simulation.


Introduction
High-rise buildings are very sensitive to dynamic loads because of their slenderness. The time-dependent types of loading that have a significant influence on this type of the structures are wind and earthquakes. Due to their nature they should be treated as non-deterministic forces. Dynamic wind loads or earthquakes lead to large sway motions of high-rise buildings. If the displacements at the top of the structure are significant, it may lead to damage of structural elements and, in extreme cases, to the building being destroyed. Therefore, in the design process we need to satisfy all requirements to make sure that the structure is safe. Designers are constantly trying to increase the load-bearing capacity of such systems and the comfort of people staying inside through the use of innovative solutions and materials or methods of monitoring. Some examples can be found in [1,2].
In high-rise buildings, an efficient system for transporting people and equipment should also be provided due to the significant height of the object. The sway of the structure caused by dynamic loads leads to the excitation of the cables inside the lift installation. As a result, the transverse and longitudinal vibration of the cable together with the horizontal and vertical displacements of the lift car can be observed. Lift cables are made of highstrength steel wires that are twisted together to form a strand structure. In the resonance region the significant vibrations may disrupt operation of the lift system, result in fatigue of the material or even damage of the rope. In high-rise buildings, special structural elements, such as spring-viscous damping elements or tuned-mass dampers, are used to limit the negative effects in the lift's systems caused by vibrations. However, their correct design depends on a careful examination of the structural system behavior and an estimation of the risk of resonance. If we wish to verify the fatigue resistance of cables or their load-bearing capacity we need to be able to define the response of the lift installation to the dynamic excitation. Very often, external influences are changeable and in the design process the worst-case scenario has to be found. This is why the methods of computation that can quickly give results for structures under stochastic loads are so important.
Various models describing the behavior of the lift installation under deterministic and stochastic dynamic loads can be found for example in [3,4]. Earthquakes can bring similar problems in high-rise building that can cause the damage inside the lift installation [5][6][7]. Even if the source of the long-period ground motion is in great distance from the building it can lead to significant displacements of its base in the resonance region [8]. In [9,10] the problem of the cable-mass system dynamic response to the earthquake excitation idealized as deterministic harmonic process was considered. However due to the nature of the phenomenon the ground motion should be dealt with by stochastic methods. In the presented paper, the proposed non-linear model of cable-mass-spring system is assumed in the form of a cable with a concentrated mass attached to its lower end. The system is moving slowly in the vertical direction. The horizontal displacements of the mass are constrained by the viscous-damping element. The transverse vibrations of the cable are coupled with the longitudinal displacements. The earthquake excitation of the building base is modelled as a filtered Gaussian white noise stochastic process.
Analytical solution of non-stationary and non-linear problems is usually difficult [11]. Therefore, approximate methods and numerical techniques are often applied. One of the techniques that is commonly used to determine the stochastic response of nonlinear systems in different problems is the Monte Carlo direct simulation technique [12][13][14][15]. However, to obtain reliable results, this technique requires usually the generation of a large ensemble of excitation and response sample functions, hence the computational cost may be very high. Therefore, in this paper, the approximate analytical technique is developed, where the nonlinear model is replaced by the approximate linear one by using the equivalent (statistical) linearization technique. The particular coefficients of an equivalent linear system expressed in terms of expectations of non-linear functions of random state variables are obtained from the condition of mean-square minimization of the error between both systems. An implementation of equivalent linearization technique to obtain the response of a simplified model of lift installation under stochastic wind load can be found in [16].
The equivalent statistical linearization technique was first used to consider the nonlinear random vibrations by Caughey [17]. Since that time, the technique has been implemented in many problems of stochastic dynamics, for example, [18][19][20][21]. The statistical linearization technique has been applied to solving various problems in systems under different random excitations, not only Gaussian, for example, [22]. It has also been combined with other advanced methods and techniques to conduct analysis in various domains of stochastic dynamics and more. Some examples can be found, for example, in [23][24][25][26][27][28].
The equivalent linearization technique is relatively easy to apply and allows us to obtain sufficiently accurate estimates of the mean values, variances and covariances of particular random response state variables. The advantages of this approach and its effectiveness are discussed in the present paper on the basis of a comparison with the results of the solution of an original set of nonlinear differential equations governing the behavior of the system under an idealized deterministic harmonic excitation. Also, the obtained approximate results are verified against Monte Carlo simulation for the filtered Gaussian white noise process excitation.

Schematic Cable-Mass-Spring Model
Consider a simplified model of a lift system represented by a cable-mass-spring system presented in Figure 1b is considered. The cantilever host structure of the total height Z 0 is exposed to the ground motion s 0 (t). The system is moving downward with the transport speed and acceleration defined as V and a, respectively. The length of the cable L is changing during the motion. The main mass M attached to the bottom end of the cable is constrained in the horizontal direction by a spring-viscous damping element with the coefficients of stiffness and damping denoted as k and c, respectively. The lateral displacements of the main mass M are defined as v M (t) while the longitudinal as u M (t). The ground motion leads to a sway of the structure, as a result the elastic bending deformationsw(z, t) can be observed. The displacements at the top of the structure, where the z = Z 0 , are denoted asw 0 (t) =w(Z 0 , t) and lead to the excitation of the rope inside the lift installation. The lateral vibrations of the cable v(x, t) are coupled with the longitudinal dynamic displacements u(x, t), where x is defined as the Eulerian spatial coordinate. The cable inside the lift installation is made from the high-strength steel wires that are twisted together automatically (see Figure 1c). Its mass per unit length, cross-sectional area and modulus of elasticity are assumed as m, A and E, respectively. The mean quasi-static tension inside the cable is defined by

Non-Linear System
First, let us consider the dynamic response of the system to a base motion s 0 (t). The overall (absolute) displacements w(z, t) of the structure are governed by the equation of motion: where m s (z), C, L are the linear mass density of the structure, damping operator and spatial operator related to the elastic potential energy, respectively. The symbol () t denotes the partial derivative with respect to time.
The overall displacements of the structure presented in Figure 1b can be given by the formula: The first and second partial derivatives with respect to time are obtained as: Using Equations (2) and (3) in Equation (1) brings: Since s 0 (t) is not a function of geometric co-ordinates, Ls 0 (t) = 0. If damping only depends on the relative motion and it does not depend on the absolute motion, then Cṡ 0 (t) = 0. Therefore, Equation (4) can be rewritten to the following form: whose solution can be assumed in the form of modal expansion: W n (z) and p n (t) are the eigenfunctions of the structure and natural (modal) coordinates, respectively. The first and second partial derivatives of Equation (6) with respect to time are given by:w Using Equations (6) and (7) in Equation (5) results in: Multiplying Equation (8) by the W r brings: By integrating the result over the domain 0 ≤ z ≤ Z 0 and using the eigenfunction orthogonality conditions: Equation (9) is transformed to the following form: (11) Dividing all terms of Equation (11) by the modal mass defined by the expression: the modal equations are obtained as: where ω r and ζ r denote the natural frequencies and modal damping ratios, respectively [9]. P r (t) are the modal excitation terms. If the natural frequencies of damped vibrations are determined as ω dr = ω r 1 − ζ 2 r , the steady-state response of the structure is defined by: If we are looking for an approximate solution of the Equation (5), we can assume that: [29] where Ψ n are the approximating functions and p n (t) denote the general coordinates that correspond to the lateral motions of the structure. The first and second partial derivatives of Equation (15) with respect to time are given by: Using Equation (16) in Equation (5) leads to: Multiplying both sides of Equation (17) by Ψ r and integrating the result over the domain Assuming M, C, K and F(t) as the mass, damping and stiffness matrices and the vector of generalized excitation forces, respectively, expressed by the following equations: Equation (18) is rewritten to equation of motion in the matrix form: As is well known, the sway of the structure leads to the excitation of the rope inside the lift installation. The cable axial strain is assumed to be ε = u x + 1 2 v 2 x and the quasi-static tension of the cable is given by the expression T = (M + mL)(g − a). From the Hamilton's Principle for external work of non-conservative forces, kinetic and potential energies of the system, the following partial differential equations of motion are obtained: where ∆ denotes the deformation of the spring expressed by being the bending deformation of the host structure at the level of the main mass. The derivation of Equation (21) in detail can be found in [9]. The total derivatives with respect to time are given by the equations: Because the lateral frequencies of tensioned metallic cables are much lower than the fundamental longitudinal frequencies and in the case of long-period ground motions, the excitation frequencies are much lower than the fundamental longitudinal frequencies, and the longitudinal inertia of the cable in the first equation of (21) can be neglected [30]. Integrating this equation using boundary conditions u(0, t) = 0 and u(L, t) = u M (t) leads to the following expression [9]: where e(t; τ) represents the averaged over the length quasi-static axial strain in the rope expressed by: This results in reducing the dynamic model described by Equation (21) to three equations of motion. The change of L(t) over a period T 0 corresponding to the fundamental frequency of the system f 0 is small in comparison to the total length of the cable [31,32], therefore the slow time scale is defined as τ = t. The small parameter << 1 is expressed by the equation =L(t 0 )/ f 0 L 0 [33], where t 0 denotes a given time instant corresponding to f 0 and L 0 = L(t 0 ).

Base Excitation and the Nonlinear Response of the Cable-Mass-Spring System
The earthquake leads to the base motion excitation. As a result, the sway of the structure can be observed that causes the bending deformations of the host structure, which is treated as a cantilever beam. In the resonance region, they can be described approximately by the polynomial shape function Ψ(η) = 3η 2 − 2η 3 , where η = z Z 0 . In this consideration, it is assumed that the impact of the cable system vibrations on the structural response can be neglected. The time dependent overall lateral displacements of the cable-mass system can be expressed by the equation: where Ψ L andw 0 (t) are adopted in the form: The relative lateral displacements can be described by the finite series as: with orthogonal trial functions expressed by the formula Φ n [x, L(τ)] = sin[σ n (L(τ))x], n = 1, 2, . . . , N where N denotes the number of considered modes. The q n describes the generalized coordinates that correspond to the lateral motions of the cable system. The eigenvalues σ n (τ) are slow varying and can be defined by the equations: Using Equations (25) and (27) in the reduced system of Equation (21) leads to the differential equations in the form [9]: where the symbolsζ r and ζ M denote the damping ratios for the cable lateral mode and for the mass longitudinal mode, respectively. The particular expressions included in Equation (29) are defined by: To describe the system's behavior in the resonance region, a single-mode approximation for r-th mode with displacements assumed in the form presented by Equation (27) can be considered. Then, the equations of motion are obtained as:

Stochastic Modelling of the Ground Motion Due to an Earthquake
The structure, due to an earthquake, is subjected to a base motion s o (t), which is represented as (see Figure 2): where G(t) is the ground motion relative to the bedrock and R(t) is the absolute motion of the bedrock. Between the bedrock and the ground surface there is a layer of soil. The properties of the soil layer are idealized as a single-degree-of-freedom system with the mass m soil , stiffness coefficient k soil = ω 2 s m soil and the damping coefficient c soil = 2ζ s ω s m soil . In the model presented in Figure 2 subscript r and corresponding modal mass m r , given as: refer to the fundamental mode of the structure, where m s denotes the linear mass density of the structure. Using the approximations given by the equation: and considering the response of the structure in the fundamental resonance, the fundamental modal shape function is expressed by the polynomial shape function W r (z) ≈ Ψ(η).
The modal mode approach leads to the formula (see Equation (13)): therefore, the displacements at the top end of the structure are defined asw(Z 0 , t) = w 0 (t) = p r (t) and consequentlyẇ 0 (t) =ṗ r (t),ẅ 0 (t) =p r (t). The model shown in Figure 2 represents a 2DOF system, where the ground motion is coupled with motion of the building structure and is as follows: Noting that w 0 =w 0 (t) + s 0 = p r + s 0 and s 0 = R + G we obtain: Dividing both sides of the Equation (37) by m soil leads to the following formula: If we assume that the ground motion is not influenced by the building structure we get: Noting that c soil m soil = 2ζ s ω s and k soil m soil = ω 2 s the Equation (40) is transformed to the following form: The negative acceleration −R(t) of the bedrock is assumed as a Gaussian white noise ξ(t), hence: Stochastic differential equations governing the state vector G = [G 1 , G 2 ] T are: where W(t) is a Wiener process and P is the constant spectral density of the white noise process ξ(t). The state vector G must be appended to the state vector of the dynamic system and in the equations of motion the following replacement must be done: Using Equation (36) together with the Equation (44) in Equation (31) leads to a set of equations of motion given by:

Equivalent Linearization Technique Implementation
Converting the second-order differential equations to the first order differential ones is performed by using the following expression: where W(t) denotes the standard Wiener process, c(Y(t), t) is the drift vector and d(t) means the diffusion vector. The augmented state vector Y(t) is defined by: The particular elements of drift vector are expressed as: The augmented state vector transformation to the centralized state vector is required to convert the original nonlinear set of differential equation into the linear one by using the equivalent linearization technique. The centralized state vector is defined as: The differential equations for the mean values are given by: Using the centralized state vector leads to the stochastic equation in the following form: where the centralized drift vector c 0 (Y 0 (t), t) is expressed by: and diffusion vector, independent of the state vector, is defined as: The particular elements of vector c 0 (Y 0 (t), t) are given by the following formulae: In further considerations, the original nonlinear system described by Equation (46) is replaced by the linear one expressed by: The centralized drift terms are defined as a linear function of the state variables: and equivalent coefficients B im , determined from the condition of mean-square minimizations the difference between the Equations (51) and (55), are obtained as: or in matrix form: The centralized state variables Y 0 are jointly Gaussian distributed, therefore in further consideration, the relationship for zero-mean Gaussian random vector X, given by [34], is used: where non-linear function is denoted by f (X) while ∇ is given by the following expression Equation (60) describes the components of matrix B as: The result of applying Equation (61) to the elements of the centralized drift vector defined by the Equation (54) is matrix B obtained as: where the particular terms are given by: To obtain variances and covariances of particular random state variables the following set of differential equations for covariance matrix together with the differential equations for mean values defined by the Equation (50). As a result the set of 44 differential equations is obtained that can be solved numerically.

Main Assumption for Numerical Calculation
In the simulation study and analysis, the case of the lift moving downwards from the building's top to the base level is considered (see Figure 1a). The main mass of the lift car is assumed to be M = 3600 kg and it is attached at the lower end of the cable system comprising n r = 6 wire ropes (compare Figure 1c). For every rope, the mass per unit length is equal to m = 0.872 kg/m, while the longitudinal stiffness is EA = 22.889 MN. The values of the transport speed and acceleration taken into consideration are V = 2.5 m/s and a = 1 m/s 2 , respectively. The total height of the host structure and the initial length of the cable are adopted as Z 0 = 258.66 m and L 0 = 58.66 m. That brings the value of the travel height H = 200 m. The horizontal displacement of main mass is constrained by spring-viscous damping element with the stiffness k = 66.689 kN/m and damping coefficient c = 9.297 kNs/m. The natural frequency of the main mass is estimated as ω M = √ k/M = 4.3043 rad/s, which brings tof M = 0.685 Hz, while the fundamental longitudinal natural frequency is defined as ω M = √ EA/ML. In this case study, the fundamental mode of the cable-mass system has been used in the single-mode approximation. The damping ratio for the cable lateral mode, the structure damping ratio and damping ratio for the mass lateral and longitudinal mode are assumed asζ r = 0.003, ζ r = 0.025 and ζ M = 0.3, respectively.

Nonlinear Results under Harmonic Excitation via Equivalent Linearization Technique
In nonlinear computation, the ground motion s 0 (t) is assumed in the form of harmonic motion defined by the equation s 0 (t) = A g sin(Ω g t) where the excitation frequency is adopted as Ω g = 0.68 Hz, which is related to the main mass natural frequency. The acceleration magnitude is equals 0 = 0.1 m/s 2 . Therefore, the amplitude of the ground motion is obtained as A g =s 0 /Ω 2 g = 0.0055 m. Using Equation (13) and assuming that m s (z) ≈ms = const and m r =ms Z 0 0 Ψ 2 (z)dz the building response can be described by the formula: For r = 1 the phase angle ε → π 2 . The amplitude at the top of the building is given by: where R = Ω g ω r . In the resonance region, R → 1, and then Equation (64) is reduced to the following form: That gives a result of maximum displacement at the top of the structure with a value of 0.15 m. The system of differential equations defined by Equation (31) is then solved numerically by using the 4th-5th order Runge-Kutta algorithm. Time-average of the square ofs 0 as the harmonic process is equal.
In the stochastic approach, the ground motion is given by the Equation (44) together with Equation (43). To compare the results obtained from nonlinear solution under harmonic excitation with the expected values obtained by equivalent linearization technique the variance of the stochastic processs 0 (t), that is given by the following equation: should be compared with the time-average s 0 (t) . The form of Equation (67) results from the fact that G 1 and G 2 are uncorrelated. In the steady-state, the variance is written as follows: A comparison of the right hand sides of Equations (66) and (68) leads to the quadratic equation whose discriminant is given by the formula: To satisfy the condition of positive values of quadratic discriminant, the constant spectral density should fulfill the following expression: For selected values of P the damping coefficients are defined by the formula: Therefore, the highest value of P that was taken to the analysis, corresponding to the main data and Equation (70), was 1.87 × 10 −4 m 2 /s 3 , for which ζ s 1 = 0.54 and ζ s 2 = 0.45 were obtained. The equivalent linearization technique was conducted for different values of P and ζ s , that satisfies the Equations (70) and (71), to observe the influence of adopted P and ζ s on the results of expected values ( Figure 3) and variances (Figure 4) of particular random state variables. However, for the clarity of presentation only some of them were presented in the paper. Of course, in one process, the value P and one of the corresponding damping coefficients ζ s 1 or ζ s 2 can be considered. The damping coefficient values greater than 1 were omitted in the analysis in advance. The comparison of the results obtained by non-linear solution under harmonic excitation and by an equivalent linearization technique conducted for Gaussian white noise excitation is presented in Figure 3. The total time that is needed for the lift car to travel from the top of the building to the base floor level at the transport speed specified above is about 82 s. However, for the clarity of the presentation, some diagrams presented in Figures 3-7 show the selected parts of the motion. In Figure 3 it can be seen that the expected values of the generalized coordinates E[q r ] and vertical displacements of the main mass E[u M ] are the same for every case of values P and ζ s that were taken into consideration during the analysis. Additionally, the expected values obtained from the equivalent linearization technique are comparable with the deterministic results received under harmonic excitation. The diagrams of the variances for the particular random state variables presented in Figure 4 show that the lower the value of the damping ζ s is taken into the analysis the higher value of the variance is obtained, which seems to be obvious. Some deviations from this rule can be noticed in the initial phase of the motion.

Verification Using Monte Carlo Simulation
In this section, the results obtained by equivalent linearization technique are compared with the values received from the Monte Carlo Simulation, which was conducted by using 1000 simulations for the different time steps ∆t. The set of differential equations that was solved numerically in Monte Carlo simulation by using the 4th-5th Runge-Kutta algorithm is defined by the Equation (45) together with Equations (42) and (36). The most important part in these computation is the correct generation of the white noise process with the zero mean values. Many programs for generating of random variables with the normal distribution and assumed mean value can be found and used in the numerical analysis.
The question arises what value of the standard deviation of the white noise excitation process should be assumed in the Monte Carlo simulation in order to be able to compare the results with those obtained by equivalent linearization technique. In many papers information can be found that the Monte Carlo Simulation results depend on the value of the time step ∆t taken into the analysis in such a way that, the smaller the step value, the lower the variance. In practice if the value of the standard deviation σ = 1 of the white noise excitation process is assumed in the Monte Carlo Simulation the different results of variances of the same variable for different time steps can be observed. It can be seen then that the shape of the lines is similar but the order of magnitude is completely different. So the incorrect assumption of the standard deviation value can lead to misinterpretation of the final results.
As is well known, the differential increment of the Wiener process W(t) may be represented as: where Z is the zero-mean Gaussian random variable. Hence, the Gaussian white noise process ξ(t), which is the generalized derivative of the Wiener process, is represented as: It means that, to make the results of the Monte Carlo simulations independent of the value of the time step ∆t, the generated value of the variable Z must be divided by √ ∆t.      In the Monte Carlo simulation conducted in this work, the standard deviation of the Gaussian white noise process is assumed according to Equation (73) with the value σ = 1/ √ ∆t. This leads to the variances that are comparable for every time step ∆t. During the analysis, the Monte Carlo simulation was made for different values of P and ζ s and the behavior of the system is similar in each case. Therefore, in this article only results for P = 1.87 × 10 −4 m 2 /s 3 and ζ s = 0.54 that are presented in Figures 5-7 are discussed. Figure 5 shows that the expected values of particular random state variables are almost the same in both methods. The only differences can be seen in case of E[G(t)] where in equivalent linearization technique a solid line with the value 0 is observed while in the case of the Monte Carlo simulation, the results oscillate around zero, because they depend on the Gaussian white noise process, which is generated as a sequence of random variables. From Equations (42) and (43) governing the process G(t) it follows that the exact mean value E[G(t)] = 0. Since these equations are just appended to the dynamic system and are linear, the equivalent linearization technique yields exact mean value.
In Figure 6, good matches between the variances obtained from the simulations with those obtained by the equivalent linearization technique can be seen. These results were obtained for ∆t = 0.05 s and only in the transient part of the motion can some differences between both methods be observed. Therefore, Figure 7 shows the variances obtained for the first 10 s of the motion, to depict that if smaller ∆t is taken for the Monte Carlo simulation then better matches between the diagrams from both methods can be observed. Additionally, adopting ∆t = 0.01 s allows us to avoid disturbances in first second of the motion that for the bigger values of time steps lead to incorrect results of Var[G(t)] which are made by inaccuracies due to numerical analysis. On the other hand, the fivefold reduction of ∆t brings multiple increase of the time that is needed to conduct the Monte Carlo simulation. What is worth mentioning the equivalent linearization technique is not so sensitive for time step variation because obtained results are very similar for different ∆t. Only a significant increase in the time step above quarter of fundamental period of the motion causes disturbances in the obtained values. This remark is important from the point of view of the time that is needed to conduct the computations.
It can be noticed that the values of the variances of particular random state variables depend on the response of the system. As it can be seen in Figure 5  Figures 5-7 present good matches between the diagrams obtained by Monte Carlo simulation and equivalent linearization technique, which show that the proposed linearized model is adequate. It leads to the conclusion that the equivalent linearization technique can be successfully used to obtain the mean values and variances of particular random state variables. The most important advantage of this solution is the much lower time needed to obtain results. For example, for small time steps like ∆t = 0.01 s, the total time to conduct the Monte Carlo simulation is several hundred times longer in comparison to the equivalent linearization technique. The application of the proposed method in the computer codes leads to a procedure that can be easy adapted to objects with different parameters.

Conclusions
Earthquakes lead to motions of the building base that cause bending deformations of the host structure. This results in vibrations of the top of the building, which cause the excitation of the steel wire cables and ropes inside the lift installation. In the resonance region, the amplitude of the vibration increases significantly. This can cause failure of the lift, excessive fatigue of steel wires and damage of the rope. Lifts are pivotal components of the building infrastructure and to prevent damage/ failure to the installation their behavior needs to be thoroughly investigated. Because the nature of earthquakes is nondeterministic it should be considered by using the stochastic method.
The results presented in this paper show that the proposed model of earthquake excitation represented as a filtered Gaussian stochastic process can be used to replace an idealized deterministic harmonic one. The equivalent linearization technique leads to the replacement of the original non-linear system governed by differential equations with an equivalent linear one governed by differential equations whose coefficients are expressed in terms of expectations of the non-linear functions of the response process. It allows us to obtain mean values, variances and covariances of particular random state variables, which gives a better view of an earthquake's influence on the behavior of the entire system. This information is very important for the design process. The presented procedure can be easy implemented in computer codes by using numerical techniques. In this approach, the important factor is that the time needed to obtain the results by using the equivalent linearization technique is much lower compared to that for other statistical methods. Funding: This work was carried out as a part of the project UPB-503-12-088-05/4 (supported by West Pomeranian University of Technology in Szczecin, Poland).

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

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