An integrated multi-energy flow calculation method for electricity-gas-thermal integrated energy systems

Mengting Zhu,Chengsi Xu,Shufeng Dong*,Kunjie Tang and Chenghong Gu

Abstract The modeling and multi-energy flow calculation of an integrated energy system (IES) are the bases of its operation and planning.This paper establishes the models of various energy sub-systems and the coupling equipment for an electricity-gas-thermal IES, and an integrated multi-energy flow calculation model of the IES is constructed.A simplified calculation method for the compressor model in a natural gas network, one which is not included in a loop and works in constant compression ratio mode, is also proposed based on the concept of model reduction.In addition, a numerical conversion method for dealing with the conflict between nominal value and per unit value in the multi-energy flow calculation of IES is described.A case study is given to verify the correctness and speed of the proposed method, and the electricity-gas-thermal coupling interaction characteristics among sub-systems are studied.

Keywords: Integrated energy system,Integrated multi-energy flow calculation, Newton-Raphson, Compressor, Per unit value and nominal value

1 Introduction

In users or building-level systems, through the coordination of heterogenous generation units, energy storage systems and flexible loads, multiple energies can be simultaneously generated, transmitted, stored and consumed.An integrated energy system (IES) is conducive to the rational planning and optimal operation of various energy sources, so as to give full play to the complementary advantages of various energy sources, improve energy efficiency and promote the consumption of renewable energy [1–4].Multi-energy flow calculation,which can determine the operational state of each energy sub-system in an IES for a given network structure,parameters and boundary conditions, is an important basis for the planning and optimal operation of IES [5].At present, there are complete energy flow calculation models and methods in electricity, gas and thermal energy sub-systems [6].The models and methods of power flow calculation for transmission and distribution networks are introduced in detail in [7, 8], respectively.The flow calculation methods of integrated transmission and distribution networks have also been widely studied[9, 10].In [11], steady-state and dynamic energy flow models and calculation methods of a natural gas network on the basis of thermodynamics and fluid mechanics are established.Reference [12] describes the detailed steady-state hydraulic model and thermal model of a thermal network.

An IES is a physical entity of the energy internet [13].With the relative maturity of each energy system model,the modeling and multi-energy flow calculation of an IES including the electricity, natural gas and thermal networks or just two of them have been widely studied.At present, the Newton-Raphson method [5, 14–16] is mainly used to solve IES energy flow.This can be divided into two solutions, i.e., the integrated solution and the decomposed solution [17].In [18], the two solutions are used to calculate the energy flow of an electricthermal coupling system on Barry Island, and it concludes that the integrated solution has the advantages of fewer iterations and its iterations do not increase with the size of the system.In [19], an improved practical method is proposed to reduce the complexity caused by the traditional method in calculating the natural gas system with compressors, and the integrated calculation model of electricity-gas-heat IES is conducted on this basis.In [20], a traditional method is used to calculate the gas system with compressors, and a hybrid technique is proposed to solve the energy flow problem of an electricity-gas system using a genetic algorithm (GA) to search the initial values of the gas system.In [21], a new alternate iterative calculation method considering the different characteristics of each system and the scalability of IES is proposed, and compared with the integrated solution using an example.In [22], a generalized energy flow (GEF) analysis model is proposed and a hybrid technique combining homotopy and the Newton-Raphson algorithm is used to solve the nonlinear equations of GEF.

Although there has been some progress in the modeling and multi-energy flow calculation of an IES, some challenges remain.First, users or building level IES usually include electricity, gas,and heat /cold energy.While there are many studies on an electricity-thermal coupling or electricity-gas coupling energy system, few studies consider electricity-gas-thermal energy coupling.Second, in the multi-energy flow calculation of an IES,there are many forms of energy.The use of per unit value in a power system has many advantages, while other energy systems usually use nominal values.However, there are no clear steps for data processing to solve this problem.Third, there are few studies on the energy flow calculation of a natural gas pipeline model with a compressor.

When considering users or building level IES, there are also some special requirements.First, the coexistence of electricity, gas, and heat / cold energy means the energy sub-systems are closely coupled.Second, it includes a variety of equipment, and thus the type of models needs to be determined according to specific requirements.For example, when considering the uncertainty of distributed energy output [23, 24] and the volatility of load, it is necessary to establish a dynamic or probabilistic multi-energy flow calculation model.

For an IES including electricity, gas supply and heating networks, this paper makes the following main contributions:

i.Establishes the steady-state models of electricity,heat, natural gas sub-systems and the coupling equipment.The integrated multi-energy flow calculation model and method are presented with the calculation steps given in detail.

ii.Based on the reduction concept, a simplified calculation method is proposed for a compressor which is not included in any loop and whose working mode is constant compression ratio.

iii.Proposes a simple and fast numerical conversion method to deal with the conflicts caused by the use of per unit value in power system and nominal value in thermal and natural gas systems in the process of the integrated multi-energy flow calculation.

iv.Studies the coupling interaction characteristics of the electricity-gas-thermal IES considering the source-load characteristics.

The rest of the paper is organized as follows.Section 2 presents the models of electricity, district heat and natural gas systems, and proposes a simplified calculation method for a natural gas network model with compressor.In Section 3, the integrated multi-energy flow calculation model and method are illustrated.In Section 4, a case study is given to verify the correctness and speed of the proposed method, and the coupling interaction among sub-systems is studied.Finally, Section 5 draws the conclusion.

2 Model of the electricity-gas-thermal IES

2.1 Electricity system model

This paper uses the AC power flow model, in which the node voltage is expressed in the form of polar coordinates.The node power expressions are given as:

where Piand Qirepresent the active power and reactive power of node i, respectively.Viand Vjrepresent the voltage of node i and j, respectively, while Gijand Bijrepresent the conductance and admittance between node i and j, respectively.n is the number of nodes and θijis the phase angle between node i and j.

2.2 District heat system model

Figure 1 shows a closed district heat system, which is mainly composed of heat sources, supply pipelines, return pipelines and heat loads.In the network, the heat medium (usually hot water or steam) transports heat from the heat source to the users through the supply pipelines and then flows back to the heat source through the return pipelines.

In Fig.1, mris the flow of the rth pipe, Tsiand Triare the supply and return temperatures of node i, respectively.T0iis the output temperature of node i.

Fig.1 A closed district heat system schematic diagram

The model of the district heat system is mainly composed of a hydraulic model and a thermal model.

2.2.1 Hydraulic model

The hydraulic model describes the flow law of the heat medium in the network, and also satisfies Kirchhoff’s law.It can be obtained from Kirchhoff’s first law that the node injection flow is equal to the outflow, as shown in the first line of (3).When the heat medium flows in the pipeline, the pressure difference is caused by friction,which is usually referred to as head losses expressed by hf.It can be obtained from Kirchhoff’s second law that in a closed pipe loop, the sum of the head losses of the heat medium flowing in the loop is 0, as shown in the second line of(3).

In (3), Ahis the node-branch incidence matrix of the heating supply network, m is the flow column vector of the heat pipeline, and mqis the output flow column vector of the nodes.Bhis the loop-branch incidence matrix of the heat network, whereas hfis the head losses column vector of each branch.

2.2.2 Thermal model

The thermal model describes the distribution of heat energy in the network and mainly includes the following three equations, in which the included temperatures are the relative ambient temperatures, i.e., T=T-Tawhere Tais the ambient temperature.

Equation (4) expresses the node heat power of the network, in which Φ is the node heat power column vector,Cpis the specific heat of the heat medium,Tsis the node supply temperature column vector, and T0is the user output temperature column vector.

Equation (5) describes the temperature drop law in the heat pipe, in which Tstartand Tendare the temperatures of the start and end of the pipe, respectively.KTis the resistance coefficient, m is the flow of the pipe with unit of kg/s, and L is the pipe length with unit of meter.

Equation (6) describes the temperature rule of the heat medium after mixing at the node, in which minand moutare the input and output flows of the node, respectively.Tinand Toutare the temperatures before and after the flow mixture at the node.

2.3 Natural gas system model and improved gas flow calculation method with compressor

The natural gas network model is mainly composed of gas sources, gas pipelines, compressors and loads.Because of the particularity of the compressor, the pipelines containing compressors are treated separately[19].

2.3.1 Natural gas pipeline model without compressor

The pipelines can be divided into three levels by air pressure: low pressure (0–750 mbar), medium pressure(0.75–7.0 bar) and high pressure (greater than 7.0 bar).The equations of the three levels are different [9], while this paper considers the medium pressure model.

The natural gas pipeline model without compressor can be expressed as:

where piand pjare the pressures (unit: bar) of node i and j, respectively.qris the flow of the rth pipe and Kgis the pipeline constant.sijis the characteristic quantity characterizing the flow direction of natural gas, when pi≥pj, sij=+1, otherwise sij=−1.

According to Kirchhoff’s first law, the flow continuity equation of a natural gas network can be obtained as:

where G is the natural gas load column vector.

2.3.2 Natural gas pipeline model with compressor

According to [19], the natural gas pipeline model with compressors (shown in Fig.2) can be described as:

Fig.2 Compressor model in [19]

where qcomis gas flow through the compressor, and qmaand qbnare the input and output flows of the compressor, respectively.qcpis the natural gas consumption of the compressor, kcpis the compression ratio, and kcp=pb/pa.Kgmaand Kgbnare the pipeline constants of pipes ma and bn, respectively.Tgasis the gas temperature in the pipe, qgasis the calorific value of natural gas, and s is a polytropic index.

As shown in Fig.2, in this paper, m, a, b and n are set as the nodes in the natural gas network, while nodes a and b are the entry and exit nodes of the compressor,respectively.The compressor is located on the pipeline ab.Based on this, the natural gas network model of the compressor is divided into two categories.One is the model in which the compressor is included in a loop,as shown in Fig.3, and the other is the model in which the compressor is not included in a loop, as shown in Fig.4.

Fig.3 Compressor located in a loop

Fig.4 Compressor not located in a loop

For the compressor included in a natural gas pipeline loop, because the nodes on both sides of the compressor are connected to each other, the method proposed in[19] is used to calculate the compressor flow.The specific calculation process is shown in Fig.5.Based on the different control modes of the compressor, it calculates the input and output flows of the compressor, and converts them into the equivalent load of the entry and exit nodes of the compressor, respectively.In the later calculation of gas network energy flow, the pipeline containing compressor is viewed as being cut off.

Fig.5 Flow chart for calculating pipeline flow with compressor used in [19]

For a compressor which is not included in a natural gas pipeline loop, and works in the mode of a constant compression ratio, this paper proposes a method similar to the voltage level reduction method of a power system transformer to simplify the calculation.

As shown in Fig.6, the known compression ratio of the compressor kcp1reduces the pressure of the highpressure side to the low side of the compressor.The relevant node pressure changes are:

Fig.6 Example of compressor not included in loop[19]

According to (10) and from (7) and (8), it can be concluded that the pipe flow and the node gas load also require a corresponding operation as:

As shown in Fig.7, the reduced compressor can be equivalent to node o, and its gas load is qcp’=qcp∕kcp1.However, the model of the compressor is nonlinear and changes need to be made in order to make its flow continuous.Before calculating the deviation, for the element of the node-branch incidence matrix whose row is corresponding to equivalent node o, the column corresponding to the next pipeline on, should multiply kcp1, i.e.,qma=kcp1∙qbn’+Go=qbn+qcp∕kcp1.In addition, the error caused by the reduction of natural gas consumption by the compressor, i.e.(kcp1–1)qcp∕kcp1, can be ignored as it is relatively small.

Fig.7 Equivalent natural gas pipeline model of Fig.6 after reduction

Using the reduced parameters for subsequent calculation, this method needs to multiply the parameters of the corresponding nodes by kcp1after calculation.When there are multiple compressors, it is similar to multi voltage levels of the power system, and the method is unchanged, and is repeated here.

By using the proposed reduction method, the compressor working in the mode of constant compression ratio does not need to go through the iterative calculation proposed in [19].This can simplify the calculation steps and reduce the calculation time.

The reliability of the improved method is proved by an example in Section 4.1.

2.4 Coupling equipment model

In an IES, the energy sub-systems are closely connected by coupling equipment.Common multi-energy coupling equipment includes combined heat and power (CHP)systems, gas turbines, electric boilers, etc.In this paper,the most common CHP system [25, 26] is considered as the coupling equipment.

The CHP system with natural gas as input fuel is considered to operate in the following thermal load (FTL)mode.Using the method proposed by the Public Utilities Regulatory Policies Act of 1978 [27] to calculate the efficiency of the CHP system, the relationship among electric power PCHP, thermal power ΦCHPand natural gas consumption Gincan be obtained as follows:

where cCHPis the power-to-heat ratio of the CHP system and ηCHPis the efficiency of the CHP system.

3 The integrated calculation method and data processing

3.1 Integrated calculation model of IES

According to[19],the integrated multi-energy flow calculation model of the electricity-gas-thermal IES is given as:

In (14), rows 1–2 are the power balance equations of the electricity system, rows 3–6 are the balance equations of the thermal system, and row 7 is the natural gas flow balance equation.ΔP and ΔQ are the deviations of the active and reactive power of the electrical power system, respectively.ΔΦ, Δhf, ΔTsand ΔTrrepresent the deviations of the nodal heat power, the loop head losses,the supply and return temperatures, respectively.ΔG represents the nodal flow deviation of the natural gas system.PL, QL, ΦLand GLare the given active power,reactive power, heat power and natural gas load, respectively.Ah1is the node-branch incidence matrix of the heating network, which has removed the heat source nodes.Cs, Cr, bsand brare the matrices related to the pipe flow and the output temperature of the thermal network, and their specific calculations can be referred to [12].Agis the node-branch incidence matrix of the natural gas system while Ag1is derived from Agby removing the gas source node and compressor branch.Π represents p2which is the column vector of the square of node pressure in a natural gas network, while -AgTΠ represents the square difference vector of a natural gas network.x=[V, θ, m, Tsload, Trload, Π]Tare the state variables of the IES calculation.

3.2 Data processing of the IES Jacobian matrix

Using the extended Newton-Raphson method to calculate the electricity-gas-thermal IES system, the IES Jacobian matrix can be noted as:

where Jee, Jhhand Jggare the matrices on the diagonal block, which can be derived independently from each energy sub-system.The expressions of Jhhand Jggare given in Section 6.1 and 6.2.

Considering the electricity sub-system in an IES is connected to the bulk power grid, when an internal fluctuation occurs, it will be balanced by the bulk power grid, and thus, Jheand Jgein (15) are zero matrices.Similarly, the natural gas system usually contains gas source nodes so when an internal fluctuation occurs, it will be balanced by the gas source.Thus, Jegand Jhgin (15) are also zero matrices.

Since the CHP system works in FTL mode,a fluctuation in the thermal network will affect the operational state of other systems,and thus,Jehand Jghare non-zero matrices.The model of CHP system can be expressed as:

where Ahsourceis the row related to the heat source in the node-branch incidence matrix of the heating network.

According to(15),it can be obtained that:

In the traditional AC power flow calculation of the electricity system, the per unit value of variables and parameters is adopted, which has unique advantages.For example, it is convenient to observe and compare data,and simplify the reduction of multi voltage level networks.However,the nominal value is generally used in the power flow calculation of thermal and natural gas systems.

Therefore, in the calculation process of the electricitygas-thermal IES, due to the existence of the coupling link, there will be problems when the electricity, thermal and natural gas systems use different value representations.If converting per unit values to nominal values in an electricity system to participate in the multi-energy flow calculation, there exists the problem of tedious calculation of the Jacobian matrix elements.In view of this,this paper proposes a technique that enables the electricity system to use the per unit value while standardizing the coupling elements of the IES Jacobian matrix, as shown below.

Before the deviation is calculated in the Newton-Raphson method, the output thermal power of the CHP system needs to be calculated by the given value of the state variables.According to (16).the output electric power and the gas consumption of the CHP system are then calculated according to (13), and the corresponding active power vector of the electricity system and the nodal gas load of the natural gas system are updated.As the active power of the CHP system calculated by (13) is in nominal value, it is necessary to convert to per unit value, as:

When calculating the non-diagonal elements in the Jacobian matrix, Jehand Jghcalculated from (17) and (18)are also nominal values, and thus are converted into per unit form as:

If the element in Jghgoes through the proposed reduction method process in Section 2.3, the corresponding element should be calculated by:

Using the above conversion method, compared with the power system using nominal values to participate in the calculation, it will lead to a larger difference in value between different matrix blocks.However,standardization is only on the elementary row transformation of the matrix.This will not change the singularity, and thus the results of multi-energy flow calculation remain unchanged.

3.3 Calculation steps

In this paper, the integrated solution method is used to solve the multi-energy flow of the electric-gas-thermal IES.This method has the advantages of fewer solving steps and faster calculation speed.The iterative calculation is carried out based on the Newton-Raphson method,and the specific steps are as follows:

1.Set iteration parameters including the limit on the number of iterations Tmaxand the accuracy condition of iterative convergence εE, εG,εH.

2.Input the parameters of each sub-system and coupling equipment, and set the iterative initial value of each state quantity.

3.The compressor mode is dealt with according to the compressor category and its working mode.Calculate the output heat power ΦCHPof the CHP system using the initial value, and obtain the output electric power PCHPand natural gas consumption Ginfrom the coupling equipment.Convert the relevant parameters from nominal value to standard unitary value.

4.Calculate the deviation matrix ΔF.

5.Determine whether the number of iterations exceeds the limit Tmax.If not,proceed to the next step.Otherwise, go to step 12.

6.Calculate the Jacobian matrix J.

7.Calculate the correction of the state variables Δx=−JΔF.

8.Update the state variables x=x+Δx.

9.Repeat step 3.

10.Calculate the deviation matrix ΔF.

11.Judge whether the deviation of each sub-system satisfies the iterative convergence accuracy condition,i.e., max|ΔFE|≤εE, max|ΔFG|≤εG, max|ΔFH|≤εH.If all are satisfied, go to the next step.Otherwise,return to step 5.

12.Output results.

4 Results and discussions

The studies carried out in this paper use the R2019b version of MATLAB software, and a desktop computer operating with 64-bit Windows 10, Intel Core i7-6500U CPU, 2.5 GHz frequency and 4 GB memory.Based on the above model and the integrated multi-energy flow calculation method, the park level electricity-gas-thermal IES shown in Fig.8 is simulated.The park IES system consists of an 8-node natural gas system, a 14-node district heating system and an IEEE 39 nodes electric power system.The natural gas network contains a gas turbinedriven compressor operating at constant compression ratio, and there is no loop in the network.The electric system is connected to the bulk power grid, while the CHP system couples the three electricity-gas-heat networks, working in the FTL mode.The balance node of the electric system is the node connected to the bulk power grid, while for the natural gas system it is the node connected to the gas source, and for the thermal system it is node 14 connected to the CHP system.The detailed parameters of each energy network are shown in Section 6.3.

Fig.8 Case of the electricity-gas-thermal IES

4.1 Simulation of the natural gas sub-system

The flow of the natural gas sub-system in Fig.8 is calculated separately using the method proposed in [19] and the reduction method in this paper, and the calculation time is 0.0032 s and 0.0020 s, respectively.The results of energy flow are compared in Table 1.As can be seen,the calculated results of the two methods are very close,and the maximum error is 0.12%.

Table 1 Node pressure of the natural gas sub-system calculated by the two methods

It is common knowledge that in the field of natural gas the pressure measurement error of the instrument is not as accurate as that of the power system, and the error of 0.12% in this case is negligible.It can be concluded that the method proposed in this paper is effective and fast in dealing with the compressor pipe model which is not located in the loop and operated in constant compression ratio mode.

4.2 IES simulation results

The method proposed in this paper is used to calculate the multi-energy flow of the above model, whichconverges after 26 iterations.The steady-state data of the thermal and natural gas sub-systems are shown in Tables 2 and 3, respectively, while the power flow calculation results of the electricity system are shown in Section 6.3.

Table 2 Energy flow of the thermal system

Table 3 Energy flow of the natural gas system

Table 4 Gas load data of the natural gas system

Table 5 Pipeline data of the natural gas system

Table 6 Heat load data of the thermal system

Table 7 Pipeline Data of the thermal system

Table 8 Power flow results of the electricity system

The above results are obtained by standardizing the coupling elements of the IES Jacobian matrix.Compared with the calculation results of converting the per unit system into nominal value in the electricity system, the two results are completely consistent.

4.3 Coupling analysis

To analyze the coupling relationship among the three sub-systems, system behavior when the heat load of node 13 is increased from 1 MW to 2 MW is studied.Because of the coupling between various energy networks, the operating states of the CHP system and each energy system change.

Fig.9 Pipe flow change in the thermal system

Fig.10 Supply temperature change in the thermal system

In Section 6.3, the power flow results of the electricity system before and after the heat load change are provided.The data shows that before and after the increase of thermal load, the voltages of the load nodes in the electricity system fluctuate from 0.9848–1.0450 p.u.to 0.9848–1.0454 p.u..The active power output of node 30 in the electricity system connected to the CHP system increases from 2.3853 p.u.to 3.7075p .u., while the active power output of node 39 (the bulk power grid) decreases from 1.5520 p.u.to 1.3313 p.u..As shown in Figs.9 and 10, the related pipe flows increase in the thermal network, and the maximum supply temperature increase of the nodes is 1.8°C.From Fig.11, the flows in related natural gas network pipes also increase, while the maximum nodal pressure decrease is 1.0 bar.Figure 12 demonstrates the changes of the operation of the CHP system and the bulk power grid.Because of the CHP system working in the FTL mode, when the heat load increases, the output heat power of the CHP is increased by 53.2% while the output electric power is increased by 55.5%.The natural gas consumption of the CHP system is also increased by 55.4%.At the same time, the electricity system takes less electricity from the bulk power grid.

Fig.11 Pipe flow and node pressure change in the natural gas system

Fig.12 Operation state change of the CHP system and the bulk power grid

4.4 Convergence comparison

The methods proposed in this paper and in [19] are used to separately simulate the compressor in the IES model in Fig.8.The number of iterations of the method proposed in this paper is 26 and the simulation time is 2.440429 s.In comparison, using the method in [19], the number of iterations is 55 and the simulation time is 3.974894 s.Thus, the method presented in this paper can be used to calculate for a large natural gas pipeline model with compressors not included in the loop, as it can effectively reduce the number of iterations and calculation time.

4.5 Discussions and extensions

The model and multi-energy flow calculation method of an IES proposed in this paper can be applied to actual buildings in certain scenarios.In actual buildings, the existing energy is electricity, gas, heat and cold, and its steady-state model is given in this paper.The characteristics of cold energy and heat energy are the same so their models are also the same,as both transfer energy through a heat medium.In view of the diversity of equipment contained in buildings, it is necessary to model flexibly according to the actual situation.If distributed generations,load fluctuation and other conditions need to be considered, the dynamic or probabilistic IES multi-energy flow calculation model should be considered.

5 Conclusions

This paper has established the models of various energy sub-systems and the coupling equipment for park level electricity-gas-thermal IES.On this basis, an integrated multi-energy flow calculation model of IES is constructed.By treating the compressor separately and its pipe flow equivalent to the nodal gas load in accordance with the working mode in a natural gas network,a simplified calculation method of compressor model which is not included in a loop and operates in constant compression ratio mode is proposed using a reduction method.In addition,a numerical conversion method for dealing with the conflict between nominal value and per unit value in the multi-energy flow calculation of IES is described in detail.A case study is given to verify the correctness and speed of the proposed method, and the coupling interaction among the sub-systems is also studied.

General conclusions are as follows:

1.The proposed simplified calculation method for the certain compressor category mentioned above is correct and fast, and this has been confirmed by the simulation results.The system can converge within 30 iterations, which constitutes a reduction of 52.7%, while the calculation time is also reduced by 38.6%, compared to the existing method.Thus, the method can be particularly effective when the scale of the system is large.

2.The technique of standardizing the coupling elements of the IES Jacobian matrix complements some data processing gaps in the calculation of integrated multienergy flow in an IES.This technique enables the electricity system to use per unit value to participate in the calculation and simplifies the calculation process.Results are compared to the method of converting the electricity system from per unit value to nominal value,and show good match,which proves the correctness of the proposed technique.

3.Considering electricity-gas-thermal coupling interaction characteristics,it can be concluded that when the network state of one energy sub-system in an IES changes,the operational states of other energy subsystems will be affected.Through the integrated multi-energy flow calculation,the changes of node and equipment status in each energy sub-system can be known in order to evaluate the rationality of the system,and avoid possible risks such as transformer overload in advance.This is of great significance.

6 Appendix

6.1 Mathematical expression of Jhh

6.2 Mathematical expression of Jgg

6.3 Data of the case study

Abbreviations

IES: Integrated energy system; GA: Genetic algorithm; GEF: Generalized energy flow; CHP: Combined heat and power; FTL: Following the thermal load

Acknowledgments

Not applicable.

Authors’ contributions

M.Zhu, the main author of this study, her contributions included the idea,mathematical and practical design,case study and writing the paper.S.Dong, the corresponding author, he guided the study at all stage and improved the text.C.Xu,K.Tang and C.Gu, the supervisors of the study,their reviewed and improved the text.The author(s) read and approved the final manuscript.

Funding

This work was supported by National Natural Science Foundation of China(52077193).

Availability of data and materials

All data generated or analyzed during this study are included in the published article (and its supplementary information files).

Competing interests

The authors declare that they have no competing interests.