DEVELOPMENT AND EXPERIMENTAL VERIFICATION OF THE MATHEMATICAL MODEL OF THERMAL INERTIA FOR A BRANCHED HEAT SUPPLY SYSTEM

The article describes a new method for making management decisions on heat supply in district heating systems, based on solving a sequence of recurrence relations of fi rst-order differential equations, enabling to synthesize daily schedules of heat supply in such systems. Using fi rst order differential equations, we implement real-time daily heat supply scheduling, predict the time-temperature dependence for heating water in the supply line, and we form a decision on the thermal energy delivery on the basis of this information. The effectiveness of our method is confi rmed by numerical modeling and comparative analysis of daily heat supply scheduling with the help of advanced intelligent decision making tools. For comparative analysis, we considered daily scheduling using a nonlinear regression model, a generalized regression neural network, a radial basis neural network, and a linear neural network. The effectiveness of our method was estimated on the basis of MAPE (mean absolute percentage error) and Accuracy coeffi cients. The model was recognized as most effective for which the MAPE value was maximum, and the Accuracy value tended to one hundred percent. Experimental studies showed that our proposed model has an advantage over the regression model by 1.68 times and over the neural models by more than 10.2 times when modeling for a hundred heating network sections. Thus, the main purpose of our study was to increase the accuracy of the method of making a managerial heat supply decision based on the experimental verifi cation of a mathematical model of thermal inertia of a branched district heating system.


INTRODUCTION
The problem of mathematical support for daily heat supply scheduling in district heating systems (DHS) in the Russian Federation began to take shape in the 1950s.Professor E.Ya.Sokolov [1], a leading scholar in this fi eld, theoretically substantiated the methods for heat supply scheduling.The current state of the DHS is characterized by signifi cant uneven daily consumption.This is caused by a sharp decrease in industrial consumption, increased comfort of housing under construction (an increase in the share of hot water supply), increased share of consumers equipped with sophisticated automatic control systems and the widespread introduction of commercial heat metering (due to the amendments in the regulatory framework).The irregularity of the daily heat consumption in the DHS requires introducing thermal energy regulation during the day -the daily heat control.Currently, the problem of constructing a thermal inertia model for district heating networks has been solved under the condition of averaging the characteristics of branched heating networks and representing the system as a single pipeline.Experimental verifi cation of this model at a number of DHSs of the Trans-Baikal Territory in the Russian Federation showed that its operability only for similar sections in case of insignifi cant temperature perturbations.Since modern DHSs differ in a signifi cant branching, length, and have different thermophysical properties at different sites, the use of averaged characteristics introduces signifi cant errors in the model.Therefore, the problem of daily heat supply scheduling is associated with an insuffi cient amount of research in the fi eld of algorithmization and intellectualization of the managerial decision making process concerning heat supply.In general, the current state of mathematical and algorithmic support for daily heat supply scheduling in the DHS is characterized by the absence of ready-made algorithms that take into account the real thermophysical properties of heating networks and enclosing structures, as well as their operating characteristics [1,2].These algorithms should refl ect the thermal inertia of heating networks with regard to their thermophysical properties, branching and hydraulic modes.In addition, these mathematical models should take into account the heat-accumulating properties of heat consumers [3].In general terms, thermal inertia models of heating networks have been solved only in a simplifi ed form: without regard to losses due to thermal insulation, degree of branching and characteristics of the hydraulic modes at the heating network sections.A whole class of engineering and economic tasks (heat engineering, transport, information, technical and economic optimization, etc.) can be reduced to a sequence of recurrent relations of fi rst order differential equations with a linear dependence in the right-hand side.Such tasks include thermodynamic analysis models [4][5][6][7][8][9][10][11], heat and power system models [12][13][14][15][16], combustion models [17][18][19][20][21][22], models of local heat exchange systems (for example, solar collectors [23][24][25][26][27][28]), and others.Nowadays, systems of homogeneous and inhomogeneous differential equations are one of the main tools for mathematical modeling of physical processes [14].In [15], a system of homogeneous differential equations for controlling adsorption chillers is discussed.The solution is considered by gradient steps and is not represented as a direct analytical dependence.Richard Kicsiny considers modern mathematical models for heating systems with pipes based on the differential equations of Newton's law for cooling a pipe in a number of recent papers [14,16].A fairly common task to improve the effi ciency of urban heating networks is to build their dynamic models.They are based, as a rule, on systems of differential equations, and their solution is reduced either to numerical solutions [29,30], or to a signifi cant simplifi cation of the model and considering it through one or two equations for average parameters, solved by known dependencies [31].Systems of homogeneous and inhomogeneous differential equations underlie the theory describing the combustion processes.For example, when constructing reliable mathematical models for the gasifi cation of solid fuel particles, systems with hundreds of differential equations are compiled [32].As a rule, these systems are reduced by the fi nite difference method to ordinary equations, which leads to a relative loss of accuracy and the absence of direct analytical solutions [19].The worldwide growing use of alternative and renewable energy sources has formed a separate direction in the theory of mathematical modeling concerning solar collectors.As a rule, their modeling is reduced to solving systems and individual differential equations.For example, the multiple linear regression method was improved for solving differential equations simulating a solar collector in [29].A numerical and experimental study of a fl at solar collector is presented in [23], where hydrodynamics and heat transfer in the collector panel are determined by transforming differential equations by the fi nite difference method.A mathematical model described in [33] is used for simulating transient processes that occur in a heated fl at collector fl uid.The model is based on solving a system of coupled differential equations describing the thermal processes of glazing, air gap, fl uid, etc.The system of differential equations is solved iteratively in the MATLAB program.Most of the modeling problems are reduced to the numerical solution of systems of differential equations in this system [34] and similar ones.Such studies can be exemplifi ed by modeling projects (based on numerical solutions) of hybrid solar photovoltaic-thermal systems [35][36][37][38], as well as by the projects on modeling the optimal strategy of heating control planning using green technologies based on wind energy [39].Advanced intelligent tools for modeling daily distribution schedules include neural networks based on various control signal generation algorithms, such as described in [40,41].The main disadvantages of the known approaches in solving systems of differential equations in applied intellectual problems of physical process simulation are the limited accuracy of the solution and the high requirements for computational tools implementing them for multidimensional systems, as well as a limited possibility of using the solutions obtained in more complex models.One of the solutions to this problem is to develop a direct analytical solution on heat supply in district heating systems, which is presented in the second section of this article.

METHODOLOGY Obtaining a general solution for a sequence of recurrence relations
The sequence of recurrence relations of the fi rst-order differential equations with a linear dependence in the right part is written as: (1) where A 1 ,...,A N , α 1 ,...,α N , β 1 ,...,β N are constants; -fi rst derivatives, respectively.
A partial solution to the system (1) is usually conditioned by a sequence of initial distributions of parameter y 1 ,...,y N at x=0 : The sequence of recurrence relations is a solution to the system of differential equations (1): (2) It is convenient to represent the general solution as: (4) where γ S P -a system of products of predefi ned constants A 1 ... A N , α 1 ... α N , β 1 ... β N , S1…SN, P1…PN and initial conditions.To determine the general solution (4) γ S P for three components: (5) Components γ S P can be determined as: With regard to (6-8) we will write down the solution to the initial system (1): The values of the components Ψ can be obtained by recurrent substitution with sequential integration according to the methodology of [38].
It should be noted that the values of the components γ S P are a special case at S=0:

Development of a mathematical model of thermal inertia for a branched heating system based on solving a sequence of recurrence relations of the fi rst-order differential equations
The heat balance of the network water fl ow at the end of the heating network section with temperature perturbation at its beginning can be written in the following form: (14) where: c p -heat capacity of the network water, kJ/kg•°С; ρ m -density of the network water, kg/m 3 ; t 1 -water temperature at the end of the heating network (15) section at time dτ, °С; t 1 y -water temperature at the end of the heating network section at time τ→∞ , °С; V 1 -volume of the heating network section, m 3 ; ‫ט‬ 1 -volumetric water discharge at the section, m 3 /s.Heat losses in the considered area can be accounted for by writing the steady state heat balance equation: where: t 0 -temperature perturbation at the beginning of the heating network section, °С; t 1 -y -average temperature of the network water at the heating network section, °С; ‫ט‬ 1 -volumetric water discharge at the section, m 3 /s.μ 1 -local heat loss coeffi cient of the heating network section; K 1 •π -linear heat losses; l 1 -pipeline length of the heating network section, m.Solving equation ( 14) concerning t 1 y with respect to the fact that , we get the following: where is a dimensionless comple that characterizes the ratio of heat loss through the thermal insulation of the heating network section to the water fl ow passing through it.
In view of (15), equation ( 13) can be written as: (16) Or it can be written in the form of a homogeneous fi rst-order differential equation: (17) where: -relative consumption of heating water at the heating network section.The solution to equation ( 17) is written as follows: ( Value of t y nTN can be found from the solution of the system of auxiliary balances: (22) System ( 20) is transformed with regard to the transition to a dimensionless complex: (23) System ( 21) is transformed into the recurrence relation: (24) Solution to relation (22) was obtained as follows: (25) (26) (27) The exact solution to the problem of modeling the unsteady temperature condition of heating networks, taking into account their heat-accumulating properties and the multivariate confi guration (it is necessary to consider the temperature condition of a complex of consecutive sections from the heat supply source to a specifi c consumer) can be obtained by solving a recurrently related system n TN of differential equations: (28) where In addition to system (22), the system of auxiliary heat balances is used.
To simulate an unsteady temperature condition of heating networks, taking into account their heat-accumulating properties and a multivariate confi guration, it is necessary to consider the temperature condition of a complex of consecutive sections from the heat supply source to a specifi c consumer.Let us consider an approximate solution to this problem through the averaged parameters of the heating network.To do this it is required to get a solution in the form of t nTN = f (t 0 , τ) to equation (29) (30) as well as the system of initial conditions: (31) (32) (33) The general solution to system (22) can be written as follows: (34) (35) (36) where: at S = 1, ..., n TN -1 (37) at S = 0 (38) Taking into account the designation: The general solution to system (30) is convenient to represent as: Equation ( 38) is used for daily heat supply scheduling, and in contrast to the well-known methods of solving systems of fi rst-order differential equations, the direct analytical dependence (38) can be used in more complex models.

Description of the experimental research facility
Experimental verifi cation of the mathematical model of thermal inertia for a district heat supply system was carried out in the framework of periodic tests for heat losses in the heating mains of the urban-type settlement of Yasnogorsk, the Trans-Baikal Region.A brief description of the test facility is given in Table 1.

Description
Heating main  Heating water fl ow rate is measured in the test circuit using portable fl ow meters.The heating main characteristics are given in table 2.

Figure 1: Schematic heating circuit,
A brief description of the tested network and the thermal preparation plant is shown in Fig. 1.

Test Preparation and Execution
A technical heating main test program was developed for detecting heat losses from the heat source on the basis The "temperature wave" passage along the test ring was recorded at each control point.Measurement of the heating water fl ow rate simultaneously in the supply and return heat pipelines, as well as feeding the heating network from the deaerator, made it possible to determine and compare the leakage rate from the test ring, which averaged 0.26 t/h; the deviation of the direct measurement of the heating network makeup fl ow and its calculated value determined through the difference in the heating water fl ow rates in the supply and return pipelines was 0.05 t/h.

Test data processing
The test data was processed as follows.For each observation point, taking into account the actual run of water between the observation points (see Table 3), determined by the "temperature wave" method, the water temperature values obtained from 20 consecutive measurements in the main test period were averaged (from 20.00 h. on 16.08.17through 20.00 h. on 17.08.17).
The heating water fl ow rates in the supply and return heat lines of the test ring were averaged by the indications of fl ow meters.
Table 3 also contains the calculation data for actual heat losses Q s.l. and Q r.l., average annual heat losses Q n.s.lav.an and Q n.r.l.av.an and ratios of these heat losses to the normative ones K s , K r along the supply and return lines for each section of the tested ring.The average annual heat losses at the tested network sections and the ratio of these heat losses to the normative ones were determined with an error of 4.3%.Computational error for indicators was calculated based on the measurement results by the following formula: where δ -relative computational error, %; δ Qn.s.lav.anrelativeerror of average annual loss calculation, %; δ Qt relative error of heat loss calculation determined during tests based on the measurement results of heating water fl ow rate and temperature; δG s.l., δG r.l.-relative permissible error of measured heating water fl ow rates.The relative error in water fl ow rate measurement was determined from the fl ow meter characteristics.The relative error in water temperature measurement was determined by formula: where δt -relative permissible error in temperature measurement.* water specifi c gravity is taken according to the average temperature at the beginning and end of the section in the design test mode; ** water specifi c gravity is taken according to the average actual temperature at the beginning and end of the section during the tests.

Comparison of the test data and computational results for the mathematical model of thermal inertia of a branched district heat supply system
Thermal inertia and heating water temperature values were calculated by formulas (30)(31)(32)(33)(34)(35)(36)(37)(38) for a branched district heating network when dividing the heating network into 1, 10, 20 and 100 sections.
Comparison of the test data and computational results (see Fig. 2) shows that the deviation of the computational results of the model (formula 38) from the actual values lies within the measurement error range already when a network 2881 m long was divided into 20 sections.This indicates a high degree reliability of the developed model.To evaluate the managerial decision making method for heat delivery in district heating networks based on solving a sequence of recurrent relations of fi rst order differential equations (formula 38), and compare it with the existing intelligent managerial decision making models, the MAPE and Accuracy coeffi cients were calculated using the formulas: where T'i -design temperature value calculated by formula (38); Ti -actual temperature value (see Fig. 2); n -count rate.
In the course of the analysis, the MAPE and Accuracy coeffi cients were calculated for the following methods: heat supply method proposed in the article, non-linear regression, generalized regression neural network (GRNN), radial basis neural network, and linear neural network.
At fi rst, the data obtained using formula (38) and presented in Figure 2 were divided into 2 days.Then twoday nonlinear regression equations were constructed, which are presented in Figures 3 and 4.
Based on the calculated nonlinear regression equations (see Fig. 3 and 4), the MAPE and Accuracy values were calculated for 2 days.The calculations are summarized in tables 4 and 5.
Step 1. Enter the input, output and control sample.The model of the created neural network is shown in Fig. 5.
The output of the results of the neural network operation is presented in graphical form in Figure 6.After that, the MAPE and Accuracy values were calculated for three neural networks.The calculations are summarizedin table 6.

CONCLUSION
The paper presents a direct analytical solution of the n-dimensional system of recurrently related differential equations taking into account the initial conditions applicable to the mathematical model of thermal inertia in a branched heating system.The obtained system of recurrently related differential equations allows for real-time forecasting of daily heat supply schedules.
Comparison of the test data and computational results shows that the deviation of the model computational results from the actual values lies within the measurement error already when a network of 2,881 m long is divided into 20 sections.Analysis of the data for the case when the network is divided into 100 sections shows that based on the numerical values of MAPE and Accuracy indicators the proposed model has an advantage over a non-linear regression model, a generalized-regression neural model, at least by 10.19 and more times.
Nominal diameter (DN), mm Type of pipeline laying at the section Thermal insulation design Year of pipeline laying Heating schedule, ˚С Temperature chart straightening for the needs of hot water supply, ˚С Pipeline operation time, h Pipeline repair period (month) Heating season length Closed, dead-end line 1 2.881 500.0;400.0; 300.0 above-ground pipelining on low supports and subsurface pipelining in crawlways mineral through AUG 24 From AUG 25 through AUG 14 of the technical documentation.The test program was accepted for implementation by August 10, 2017.Heat loss tests were conducted from September 15 through September 17, in full and in accordance with the technical and operational test programs.At all testing stages parameters were maintained at the heat source as close to the specifi ed values as possible: temperature in the supply heating line: design temperature -87.04 ºС; actual temperature -85.9 ºС; fl ow rate in the supply heating line: design fl ow rate -53.87 t/h; actual fl ow rate -52.91 t/h.At the fi nal testing stage, the water temperature at the SDPP outlet increased within forty minutes by 15.0°C.

Figure 2 :Figure 3 :Figure 4 :
Figure 2: Comparison of the test data and computational results

Figure 5 :
Figure 5:.General layout of the neural network

DISCUSSION
The analysis of the averaged data given in tables 4, 5 and 6, allowed us to draw the following conclusions: First, the proposed method has an advantage in terms of MAPE over the non-linear regression model by 1.63 times, over the radial basis neural network by 2.34 times, over the linear neural network by 1.9 times and has an advantage in terms of Accuracy over the non-linear regression model by 1.02 times, over the radial basis neural network by 1.04 times, over the linear neural network by 1.025 times.Secondly, the proposed method has similar results with the generalized-regression neural network and is inferior to it in terms of MAPE by 0.89 times and in terms of Accuracy by 0.997 times.Thirdly, when analyzing 100 network sections, the proposed method has an advantage over all models; in terms of MAPE over the non-linear regression model by 16.48 times, over the generalized-regression neural network by 10.19 times, over the radial basis neural network by 11.61 times, over linear neural network by 26.33 times.And it has an advantage in terms of Accuracy over the non-linear regression model by 1.03 times, over the generalized-regression neural network by 1.02 times, over the neural network of the radial basis by 1.02 times, and over the linear neural network by 1.06 times.

Table 1 :
A brief description of the test facility

Table 2 :
Heating main characteristics

Table 3 :
Determination of the time of water passage through the tested network sections

Table 4 :
MAPE and Accuracy calculation for 1 day

Table 5 :
MAPE and Accuracy calculation for 2 days

Table 6 :
Calculation of MAPE and Accuracy values for neural models