EVALUATION OF DYNAMIC PROPERTIES OF GAS PUMPING UNITS ACCORDING TO THE RESULTS OF EXPERIMENTAL RESEARCHES

Experimental studies of the dynamic properties of gas pumping units (GPU) of various types, which allowed us to obtain GPU acceleration curves and determine the parameters of the transfer function through various transmission channels of input effects, are conducted. For this, a method and a software product were developed for implementing the procedure for determining the areas of the k orders through the moments of the auxiliary function. As a result of the experiments performed on operating GPU, the acceleration characteristics of the selected signal transmission channels were obtained. To identify the dynamic properties of the GPU, a software product was developed in the MatLab environment, which in iterative mode allows each acceleration curve to be approximated by the transfer function, the order of which is selected by the user. The iterative mode allows the user to select the order of the polynomials of the numerator and denominator of the transfer function, and also to calculate the numerical values of the parameters of the selected transfer function. The software product was tested on industrial data obtained during normal start-up of theGPU. The obtained results can be used in the development of a new method for monitoringthe reliability and diagnosing the technical state of automatic control systems (ACS) of the gas pumping units and its components. The essence of the method is that a change in the technical state of the ACS or GPU affects their dynamic properties, and this, in turn, causes a change in the parameters of the transfer functions, which can be recorded after a certain period of time. Determining the ranges of values of the transfer function coefficients for various technical states of different types of ACS and GPU, it will be possible to predict the period of their operation in the future

the model and the class of models to which it can be attributed are known. Then the task of identification is to determine the unknown parameters, provided that the specified structure of the model is preset.
Development of methods for diagnosing and determining the technical state of GPU on the basis of creating mathematical models and analyzing changes in their parameters is an important scientific task.

Literature review and problem statement
To ensure reliable operation of the main equipment of the compressor station (CS), the methods for monitoring and identifying the main parameters of centrifugal superchargers (CS) are widely used, which make it possible to detect damtage at the initial stage of their development and to assess the admissibility of further operation of the CS. In this regard, the importance of developing and improving identification methods is increasing, which allows us to construct their mathematical model based on the results of research on the properties of technological objects [3][4][5][6].
A significant number of scientific works are devoted to identification methods [3]. One of the first was the work [4,5], in which the frequency identification methods were initiated according to the results of observations of harmonic signals of different frequencies at the input and output of the object. Such a method is unacceptable for the majority of industrial facilities, since the normal operation of the facility is disrupted, and a considerable amount of time is required to obtain high-quality experimental material. The authors of [7] propose to evaluate the change of the technical condition of the units by the values of the measured technological parameters of the GPU and by the characteristics of the GPU, such as noise, vibration, color of exhaust gases, air leakage along the path and, in particular, by regenerators and etc.
In [6,8], the authors describe a new method for finding damage based on the identification of system complexity. In general, the complexity of a system is defined as the amount of information that needs to be described.
The advantage of this method is the ability to obtain holistic information about the state of the system by processing a significant array of data. The effectiveness of this method is confirmed by the performed studies, however, the error of this method is not indicated, as well as the effect of each measuring signal on the total measured complexity.
The effectiveness of parametric identification methods is determined not only by the quality of functional algorithms, but also by the quality of technical means. The authors of [9] emphasize that to determine the performance losses of centrifugal superchargers, it is necessary to take account of the errors of the measuring equipment, analyze the gas composition and collect data under stable conditions of the compressor operation (when the inlet and outlet temperatures do not change by more than 1 °C for 3-5 minutes).
The authors of the article [10] investigated the relationship between the main parameters of the centrifugal superchargers and its vibration signals by performing a series of experiments on experimental installation. The results of all tests showed that the vibration level increases with increasing flow rate and the presence of a corresponding defect, but to identify damage it is necessary to do a lot of research.
The combination of several methods allows you to identify effectively the parameters of the supercharger. In particular, the use of a genetic algorithm in combination with the wavelet transform and the neural network provides the method with structural stability, global search and high speed [11].
Correlation methods are used [12] in the case when the experimental data are obtained as a result of observations of an object operated in a normal mode. In this case, the input of the object must be given a signal close to white noise, which is not always possible [13].
Regression methods can be applied in the case when the object model is represented in the state space [12,13], which greatly complicates the identification process. In addition, it is assumed that the disturbances that act on the object are additive and in their properties are white noise. Such an assumption too idealizes real objects. Therefore, such methods have more theoretical value than practical.
The identification of an object by acceleration characteristics is widespread among the active methods. For objects whose dynamics are described by differential equations not higher than the second order, the parameters of the model are determined directly from the acceleration characteristic [14]. This method of determining the parameters of the model can be successfully applied in cases where the object is under the influence of minor noise.
Recently, more and more attention has been attracted to the method called the "area method" [14]. In [15], based on the area method, a method was modified in which the degrees of the polynomials of the numerator and denominator of the transfer function are positive fractional numbers. This form of the transfer function makes it difficult to interpret the results of identification and the need arises to introduce such a concept as fractional differentiation [16].
The analysis of literary sources showed that there are a number of scientific problems that need to be further solved: -the possibility of obtaining high-quality experimental data when the desired test signals are applied to an object of study; -the difficulty of obtaining mathematical models by the analytical method and their inconvenience for practical use; -taking into account the effects of interference; -selection and implementation of numerical methods for identifying objects.

The aim and objectives of the study
The aim of the work is to estimate the dynamic properties of GPU by experimentally taken accelerating characteristics by finding the parameters of the transfer function of the main control loops of the GPU using the area method. Determining the parameters of the transmission functions for different control loops for the GPU with different technical conditions (or operating time) will make it possible to establish the operating conditions for the elements of the GPU ACS and to formulate recommendations for their further operation.
To achieve this goal, it is necessary to perform the following tasks: -to develop a methodology for diagnosing the technical state of GPU automation systems and its structural units based on the use of the method of Sima and the software for its implementation; -to develop a program for conducting experimental studies on the diagnosis of the technical state of GPU automation systems; -to process the obtained experimental data using the proposed method.

Theoretical principles of the area method. Numerical solution method
It is assumed that at a small value of the input effect on the object, the relationship between the input and output of the object is described by a linear differential equation with zero initial conditions, which in the Laplace image area will be represented by the transfer function where k -object transfer coefficient; a i , i n = 1, , b j , j m = 1,constant values -parameters of the transfer function (1).
For real objects, the relation m≤n always holds. By the form of the acceleration curve of the object, the value of m and n could be estimated. Once the order of the object transfer function is selected, the question of determining the parameters a i , i n = 1, and b j , j m = 1, on the experimentally taken acceleration curve of the object arises.
We will consider an object with self-leveling, which, with a single discontinuous input action, will be characterized by a transition characteristic H(t). We assume that at time t 0 , which is the origin of observations, H(t 0 )=0. If at time instant H(t 0 ) is nonzero, then by appropriate transfer of coordinates it is easy to attain the condition We turn from the characteristic H(t) of the object to its normalized characteristic. To do this, expression (1) is presented in the following form: The normalized transfer function w(p)=W(p)/k will correspond to the time characteristic of the object h(t). Let us find its value as t→∞, when a single hopping signal 1(t) is fed to the input. So: The characteristic h(t), for which its ordinate goes to unity as t→∞, is called the normalized characteristic [17].
Thus, it is necessary to determine the parameters of the transfer function w(t): a i , i n = 1, and b j , j m = 1, , according to the normalized transition characteristic of the object h(t).
To obtain the characteristic h(t) it is necessary to divide all the values of H(t) by the steady-state value H init of the initial value of the object.
The idea of the area method is that the inverse transfer function ŵ(p)=1/w(p) is decomposed into a Taylor series of the point p = 0. Moreover, in the function w(p), the degrees of the polynomials of the numerator and the denominator are artificially aligned, assuming that b n =b n-1 =…=b m+1 =0. Then Given the ratio (3), we will have The coefficients of the expansion of S k in the formula (4) h ave the content of areas of the k order [18] and can be calculated by the formula: The values of z(t) and z k (t) correspond to the transfer function: From equation (4), it is easy to determine the relationship between the parameters of the transfer function (3) and the areas S k . The ratio (4) can be written in the following form: ) .
Opening the brackets and equating the coefficients with the same degrees, we will get Equations (6) form a linear algebraic system in which the quantities a i , i n = 1, and b j , j m = 1, are unknown. The size of such a system is m+n.
Thus, if areas S k , k m n = + 1, , are known, then from the system of equations (6) you can determine the parameters a i , i n = 1, and b j , j m = 1, of transfer functions (3). The system of equations (6) is conveniently presented in a matrix-vector form. To do this, we rewrite the system of equations (6) in expanded form: where N=m+n.
Turning to the matrix-vector form of writing the system of equations (7), we will get where In matrix A, you can select a single submatrix of size n×n. Under it, a zero submatrix of size m×n will be placed. It is convenient to fill the rest of the matrix with columns. The first element n+1 column «-1», the second element n+1 column «-S 1 », the third -«-S 2 », fourth -«-S 3 » etc. n+1 column ends with the «-S N-1 » element. The following columns are filled with a zero offset down one element in the previous column of the matrix. To solve the matrix-vector equation (8) it is necessary to determine the area S k , k N = 1, . Formula (5) is unsuitable for practical use, since the calculation of the values S k , k N = 1, requires multiple integration. A more efficient method of S k , k N = 1, calculation is the method of moments [18], which does not require multiple integration.
Form the auxiliary function φ(t) (Fig. 1) If the function Φ(p) is expanded into a Taylor series of the point p=0, then we arrive at the following result: where Differentiating successively the last expression i-1 times by the variable p and equating p to zero, we obtain Equation (9) is reduced to this form: ( 1 2 ) After substitution of the value w(p) and Φ(p) from the formulas (4) and (10)  In order for the left side of the last equation for arbitrary values of p to be zero, the coefficients at the corresponding powers of the polynomial are equated to zero. As a result, we will get Summarizing the result, it can be argued that So, we have a system of linear algebraic equations of dimension N Fig. 1 Since Λ is a lower diagonal matrix of size N×N, it is advisable to form it as follows: place the units on the main diagonal (λ ii =1, i N = 1, ); then under the main diagonal there will be elements λ(i+k, j)=μ(i-1), i = 2, N; ., while the outer cycle will be a cycle at index i. The index k is incremented by one inside the inner loop.
The expansion coefficients μ i , i N = 1, , which are called the moments of the auxiliary function, are calculated by formula (11), replacing the infinite integration interval by the final observation time t s , which introduces a certain error in calculating the values μ i , i N = 1, . Since the function φ(t) is specified in the form of a graphical dependence, in order to find the integrals, which are defined by expression (11), it is necessary to use the methods of numerical integration [19].
Numerical integration can be carried out using the formulas of rectangles, trapezoids and Simpson formula. The smallest integration error will be when using the Simpson formula, based on approximation of the function φ(t) by the Lagrange polynomial. It is assumed that the interpolation points t i are equidistant from each other. Often in practice, the value of t i is calculated with different steps along the time axis t. Therefore, we will perform numerical integration for the case when the value of t i is not equidistant along the time axis t.
Let the value of φ(t) be determined at discrete instants of time t i , i q = 1, , where q is the number of interpolation points. In addition t i+1 -t i =h i , where h i ≠const.
If t i is an interpolation point on the x-axis, and φ(t i ) is the value of the function φ(t) at the interpolation points, then the Lagrange interpolation polynomial is defined by the following formula [20]: where φ k =φ(t k ); N 1 -the number of interpolation points. Construct an interpolation polynomial for three points τ 1 =t i , τ 2 =t i+1 and τ 3 =t i+2 in which the auxiliary function φ(t) takes values y 1 =φ(t i ), y 2 =φ(t i+1 ) and y 3 =φ(t i+2 ): y .

Fig. 2. Auxiliary function with selected interpolation points
where R -the number of intervals of the duration of observations. The moments μ i of the functions were determined similarly.
In practice, the value of N rarely exceeds five. Therefore, knowing the moments of the auxiliary function, from the equation (14) we determine the area ( 1 7 ) Once the areas S i , i N = 1, are determined, from the equaetion (8) determine the parameters a i , i n = 1, , and b j , j m = 1, of the transfer function (1) ( 1 8 ) To visualize the results of the program, it is necessary to move from the transfer function to a system of differential equations in the form of Cauchy.
Let the object transfer function be known In the case when m<n the part of the coefficients of the transfer function, which is in the numerator, is equal to zero: b n , b n-1 ,…, b m+1 .
The transfer function (19) corresponds to such a system of differential equations [20]: where x i , i n = 1, -object state variables. The value of y corresponds to the output of the object with a single transmission coefficient. To obtain the initial value of the object, the dynamic properties of which are characterized by the transfer function (19) it is necessary to increase the value of y by a factor of k. Then Y(t)=ky(t).
Unknown quantities β i , i n = 0, are found from such a system of linear algebraic equations: The vector of coefficients is found from (21). So, By known values β i , i n = 0, and input action u(t), we solve a system of differential equations (20) with zero initial conditions numerically.
To implement the developed methodology, software has been developed that allows determining the parameters of the transfer function using the experimental curve of an object.

Results of experimental studies of the dynamics of gas pumping units
As an object of research, GPU type GTK-10-4 and GTN-10 were taken, mounted on the compressor station CS-2 "Dolyna" of the GTM "Prikarpattransgaz". Monitoring of the experiments was carried out using the interfaces of the GPU ACS and the interface of the automated workplace of the shift engineer CS-2. This made it possible to monitor the current values of the technological parameters on the monitors and to monitor the operative messages to the CS operator and the changes on the mnemonic diagrams and the recovery of the retrospective information. The program of experimental research envisaged the regular launch of the GPU type GTK-10. The launch of the GPU took place on April 5, 2018. The change in the current values of the technological parameters when the GPU was launched into operation in accordance with the program of experiments was monitored by the GPU ACS SAT-01 D monitor and recorded in the archive by the standard means of the ACS.
The gas turbine is controlled using software and hardware and GPU start-up schedules: a low-pressure turbine (LPT) speed controller, a high-pressure turbine regulator (HPT) and a temperature controller, and control using the Amot fuel consumption regulator.
According to the results of the experimental study, the accelerating characteristics of GPU were obtained through the channel "fuel gas consumption -rotor speed of the high-pressure turbine". Data are obtained at the start of the GPU. A fragment of the experimental results is given in Table 1. For data processing, a sample of which is presented in Table 1, a program that converts the time format "date. month.year", "hour.minute.second" into the time interval from the moment of observation to its end in seconds (minutes) was developed in MatLab environment. Fig. 3 shows the result of the operation of the software for determining the dynamic properties of the channel "GPU inlet pressure -GPU outlet pressure". The following transfer function was selected: 1 .
The acceleration characteristic was approximated by the following transfer function: To study the dynamic properties of the pressure transmission channel of the natural gas blower "inlet pressureoutlet pressure" in the GPU start-up mode, the pressure change at the GPU outlet was recorded. At the start of the observations, the pressure at the blower outlet was 38.07 bar, and at the end of the observations, the pressure rose to 43.8 bar. So, the pressure increase was 5.71 bar. This change in pressure over time is shown in Fig. 6.
With the help of the developed software, the obtained parameters of the transfer function were obtained where a 1 =129954,5 s, a 2 =4,3465×10 9 s 2 . Thus, as a result of experimental studies, the transfer functions of individual signal transmission channels were obtained, and their parameters were calculated using the developed software.

Discussion of the results of studies of the dynamic properties of GPU
To assess the technical condition of the GPU, it is necessary to form certain diagnostic features. The analysis of the literature and experimental studies have shown that dynamic characteristics change, when the technical state of the GPU changes. Such a change leads to a change in the parameters of the transfer functions of the corresponding signal transmission channels. This means that changing the parameters of the transfer functions of the GPU can be taken as diagnostic signs of the technical state of the GPU.
At the moment there are no analytical models on the basis of which one can obtain the corresponding transfer functions of the GPU. Therefore, a method and appropriate algorithmic support were developed, which allows us to determine the structure of the transfer functions, as well as its parameters, based on the obtained acceleration characteristics.
The determination of the parameters and structure of the transfer functions by the acceleration characteristics of the GPU is based on the area method.
The experimentally obtained acceleration characteristics during the start-up of GPU via signal transmission channels made it possible to determine the parameters of transfer functions with sufficient accuracy for practice.
The results obtained in the work can be used as the basis for a new method for diagnosing the state of an automatic control system for changing the type of transition curve based on direct quality indicators: overshoot, permissible number of oscillations and time of transition.
To do this, with the introduction of the GPU ACS in operation, or after repair work, they remove the acceleration curve of the GPU, which determines the transfer function coefficients, which will correspond to the nominal (operational) state of the ACS. During the operation of the GPU, after certain periods of time, they remove the acceleration curve and compare the obtained values of the transfer function coefficients with their nominal values, which will allow establishing the presence of changes in the technical condition of the ACS or the GPU itself. If there is a change in the technical state of the GPU ACS, the reasons for this change (search for a defect) are determined. The basis of any method of diagnosis lies in the working conditions according to the selected diagnostic features. In order to form the working conditions according to the criteria considered in the article, it is necessary to collect statistical data for each type of control system and with different operating time periods. On the basis of certain working conditions, it will be possible to predict the period of further trouble-free operation of the ACS and the GPU.