Robust state estimation of electricity-gas-heat integrated energy system based on the bilinear transformations

A comprehensive energy system state estimation is of great signiﬁcance to promote efﬁcient, stable and reliable operation of the system. Based on coupling characteristics of electricity-gas-heat multi-energy ﬂow and bilinear transformation method, this paper proposes a robust state estimation method suitable for integrated energy systems. First, the state estimation models of electricity, gas and heat energy subsystems and coupling com-ponents are constructed. Then, by introducing pseudo state variables and intermediate variables, the nonlinear state estimation process is divided into four stages. Among them, the ﬁrst and third stages are solved by WLAV method and WLS method respectively, and the second and fourth stages are simple nonlinear transformation processes. Finally, the robust state estimation method is veriﬁed on a 40-node simulation example. The bilinear robust state estimation method proposed not only solves the difﬁculty of selecting initial values of heterogeneous networks and the problem of non-convergence of algorithms, but also can search and identify various types of bad data, and has high computational efﬁciency and good robustness.


INTRODUCTION
Since the beginning of the 21st century, with the rapid growth of population and industrial development in all countries, fossil energy is facing exhaustion, environmental pollution is becoming more and more serious, and the institutional contradictions between extensive energy utilization mode and environmental protection have become increasingly apparent [1]. In order to promote the development and consumption of clean energy and improve the utilization efficiency of fossil energy, it is imperative to develop integrated energy system (IES). Since integrated energy systems are composed of a variety of heterogeneous energy subsystems, the network structure and operation modes are highly complex, thus a higher level of automation is required for optimal system dispatch. In order to develop an integrated energy system, it is necessary to accurately, quickly and comprehensively master its real-time operation state, timely find the problems in the system operation, effectively predict the future trend, and put forward accurate solutions to ensure stable and efficient operation of the system. At present, research on the state estimation of integrated energy systems is still unexplored, mainly focusing on the combined analysis of electricity-gas or electricity-heat coupled networks. Few studies have conducted a comprehensive analysis of all three energy networks-electricity, gas and heat. For integrated systems that take into account electricity, heat, gas, and even more types of energy, due to the differences between heterogeneous networks, the coupling system is complex and may result in long calculations or complicated calculation processes. In addition, it is difficult to determine the initial values of heat and gas networks.
The state estimation research of the integrated energy systems mainly focuses on the following two aspects: state estimation modelling and state estimation methods.
1. State estimation modelling: For a single energy system, by combining the predicted power injection with the power balance equation, a new state transition model of the power system is derived in [2], which clearly expresses the relationship between the predicted state vector and the predicted power injection; In [3], the Darcy-Weisbach equation is used as an approximate model for the momentum dissipation of the gas network pipeline to construct a spatially discrete nonlinear partial differential equation. On the basis of summarizing the existing heat network models, the semi-joint static model was adopted as the model of the state estimation of the heat network in [4], and the state estimation was successfully introduced in the heat network. For the two energy coupling networks, the characteristics of the electric-gas coupling elements are considered in [5], and the steady-state estimation model of the electric-gas coupling network for complex gas networks is established. In [6], the mathematical model and corresponding simulation model of combined heat and power plant (CHP) under the constraints of power grid and heat grid are proposed. For multiple energy coupled networks, the energy hub proposed in [7], to a certain extent, can realize the description of multiple energy coupling relations, has strong applicability to the comprehensive energy system of different scales and levels, and has a certain guiding role for the planning, design and operation control of the comprehensive energy system. 2. State estimation methods: In [8], the weighted least squares (WLS) method is used to estimate the state of the electricheat-gas coupled IES multi-energy flow model, but this method is complicated and not applicable to large-scale integrated energy system. In [9], a comprehensive mathematical analysis is presented to obtain accurate distribution network state estimation by using weighted least square estimation to find the most important factors affecting the accuracy of voltage curve. For weighted least absolute value (WLAV), a bilinear robust state estimation method based on WLAV supplementing two energy sources is proposed in [10,11], which provides a good reference for state estimation of multi-energy flow systems. In [12], in order to further improve its applicability, a systematic method was proposed to calculate the optimal rotation angles and scaling factors of the weighted least absolute value state estimation with transformations, so as to derive a robust WLAV estimation method with the optimal transformation. In addition, based on the WLS and WLAV methods, an artificial intelligence (AI) algorithm, the particle swarm optimization (PSO) is proposed in [13,14] to solve the problem of state estimation, but this method is currently only used in power systems.
In light of the conducted literature review, this paper analyses the power flow model of the existing electricity, gas and heat coupling system to construct a corresponding measurement model. In view of the non-linear characteristics of the measurement model, the pseudo-state variable is introduced to linearize the measurement model and convert it into a WLAV model for solution. Then, the estimated value of the state variable could be obtained after two more nonlinear changes and one linear change. In order to test the robustness of the state estimation method, consider using the standard residual search identification method to test. Finally, a case is given to verify the applicability and validity of IES state estimation method.

Electricity network model
In the power system, ac mode is generally adopted as the electricity network model: where ΔP i and ΔQ i are the unbalance for active power and reactive power; P i and Q i are the active power and reactive power of node i respectively; u i and u j are the voltage amplitude of node i and node j, respectively; i j = i − j ,. i and j are the voltage angles at nodes i and j respectively; G i j and B i j are the conductance and susceptance of branch ij.

Natural gas network model
The steady state model is chosen for the natural gas network model in this paper. The natural gas network models generally include pipeline equation and non-pipeline equation [15,16].

Pipeline equation
The steady state flow equation of a pipeline in the natural gas network is specifically expressed as: where f i j is the flowrate of pipeline i j; p i and p j are the pressure at nodes i and j respectively; K i j is the constant of the pipeline i j; s i j is the flow direction of the pipeline i j, if p i > p j ,s i j = 1, otherwise s i j = −1.
The pipeline constant is usually calculated by the Panhandle 'A' formula, which is specifically expressed as: where T n is the standard temperature; p n is the standard pressure; D i j is the internal diameter of pipeline ij; F is the friction factor of pipeline ij; L i j is the length of pipeline ij ; T a is the average gas temperature; Z a is the average gas compressibility factor; G is the gas specific gravity.

Non-pipeline equation
As the main non-pipeline component in the natural gas network, the driving power of gas turbine compressor (GTC) is provided by the gas turbine consumption of natural gas. The specific consumption of natural gas flow is expressed as follows: where R k is the compression ratio; f com is the flowrate through compressor; q gas is the calorific value of natural gas; a is the specific heat ratio. GTC has four modes: Constant pressure at the inlet node, constant pressure at the outlet node, constant flow through the compressor, and constant compression ratio. The specific compressor model is given in Appendix B. In this paper, only constant pressure at the outlet node is considered: where p o is the pressure of outlet node;p o is the constant value of the pressure at the outlet node.

Balance equation
Analogous to Kirchhoff 's current law and voltage drop relationship, there are also flow balance and pressure balance equations in the natural gas network: where A g is the node-branch incidence matrix; f is the flowrate of pipeline; Π i = p 2 i is the pressure squared vector of node.

Heat network model
Heat network models are usually divided into hydraulic and thermal models. The specific operation mode is given in Appendix B. Variables of the hydraulic model include pressure and mass flow rates. Variables in the thermal model are supply and return temperatures and heat power [17,18].

Hydraulic model
Hydraulic models mainly include flow continuity equation, head loss equation and loop pressure equation [19]: 1. Flow continuity equation: The continuity of flow is shown as the flow that enters into a node is equal to the flow that leaves the node plus the flow consumed at the node: (8) whereṁ in andṁ out are the mass flow rate within a pipe coming into and leaving the node respectively;ṁ qi is the mass flow through each node injected from a source or discharged to a load.
1. Head loss equation: The relationship between the mass flow of each pipeline and the head loss is expressed as: whereṁ is the mass flow within each pipe; is the vector of resistance coefficients of each pipe. The relationship between the mass flow of each pipeline and the head loss is expressed as: where L is the pipe length; f is the friction factor; D is the pipe diameter; is water density; g is gravitational acceleration.
1. For any closed loop in the hydraulic model, the sum of the head loss must be equal to zero:

Thermal model
The thermal model is used to solve the supply temperature and return temperature of the nodes of the heat network. The heat power calculation equation of each node is: where Φ is the vector of heat power consumed or supplied of each node; T s is the supply temperature; T o is the outlet temperature; T r is the return temperature; Cp is the specific heat of water. Similar to the network loss in the power network, the heat network pipeline also has heat loss, and the temperature drop equation of the pipeline is: where T a is the ambient temperature; T end and T start are the temperatures at the end and start node of a pipe; is the heat transfer coefficient of per unit length at each pipe. The calculation equation of mixing temperature at the node is: (14) where T out is the mixture temperature of a node; T in is the temperature of flow at the end of an incoming pipe.

Coupling element model
In order to realize the mutual conversion of electricity energy, natural gas energy and heat energy, this paper uses CHP and power to gas (P2G) as coupling elements.

CHP coupling equation
This paper considers natural gas as the fuel for CHP to achieve the coupling between electricity, heat and gas. The specific expression is: where c m is the conversion efficiency of heat energy and electric energy; P CON is the electric power of CHP unit in full condensing mode; F in is the natural gas consumption; Φ CHP and P CHP are the electrical and thermal power at the CHP connection node respectively; is the conversion ratio of electric energy and gas consumption.

P2G coupling equation
P2G can convert the surplus power in the electricity network into natural gas for storage and transmission to improve the efficiency of energy consumption [20]. The specific expression is: where F s,P2G and P d,P2G are the gas consumption and electric power at the P2G connection node respectively; C LHV is the low The IES network structure coupled by CHP and P2G is shown in Figure 1.

IES measurement model
To conduct state estimation of IES, the vectors of measurements and the vectors of state variables of each subnetwork should be determined first. For electricity network, select x e = [u 2 i , i ] T as the vectors of state variables and z e = [u 2 i , P i , Q i , P i j , Q i j , P ji , Q ji ] T as the vectors of measurements. The electricity network measurement model is specifically expressed as: For natural gas network, select x g = [p i ] T as the vectors of state variables and z g = [p i , f i , f i j ] T as the vectors of measurements. The natural gas network measurement model is specifically expressed as: The constraint model includes pressure constraint of compressor outlet node, CHP and P2G constraint etc.

Solution process of model
For the IES measurement model above, since the relationship between the vectors of measurements and the vectors of state variables is nonlinear, if WLS method is adopted to calculate, although it has high estimation accuracy without bad data, the following problems arise: 1) It is difficult to select the initial values of heterogeneous energy networks such as gas network and heat network, which may lead to non-convergence of calculation; 2) The Jacobian matrix needs to be updated during each iteration, which requires complex calculation and takes a long time, and the calculation efficiency is not high; 3) The WLS method itself does not have robustness, so it is necessary to add a bad data identification link based on the largest normalized residual (LNR) in the future, which further increases the algorithm time. Therefore, it is necessary to find a simpler and effective state estimation method for the integrated energy system. Through the analysis of the advantages and disadvantages of the above two methods, the WLAV method and the WLS method are used to complement each other. Based on the variable substitution method, by introducing intermediate variables, the nonlinear state estimation process of the integrated energy system of electricity, gas and heat is divided into the following four stages for solving respectively. Among them, the first and third stages are linear changes, and the second and fourth stages are simple nonlinear changes, which is the origin of bilinear transformation.

The first stage
For non-linear vectors of measurement, a pseudo-state vector y can be introduced to construct a linear WLAV model for easy solution. The specific pseudo-state vector of each sub-network is as follows [21,22]: The pseudo-state vector y e is selected in the electricity network as: The pseudo-state vector y g is selected in the natural gas network as: The pseudo-state vector y h is selected in the heat network as: Based on the linear models of the above three subnetworks, the WLAV estimation model at this stage can be constructed as follows: where is the weight coefficient; C is the Jacobian matrix with constant coefficients, and its specific form is given in Appendix C; z is the vector of measurement; m is the number of measurements; n is the number of state variables; l is the total number of constraints; E is the Jacobian matrix corresponding to the compressor constraint, CHP coupling constraint, and P2G coupling constraint, which has been given in Equation (22).
Since the residual = z − Cy, the above formula can be converted into: The WLAV model can be solved by using the linear programming method, but the absolute value of the objective function needs to be removed before this, so relaxation factors s and t are considered to replace the residual in the objective function. By making s=(| | + )∕2 and t =(| | − )∕2, the Equation (27) can be converted to the standard form as follows: where y (s) and y (t ) are auxiliary variables corresponding to s and t , and this kind of transformation is made to facilitate subsequent solutions. n is the number of state variables;l is the total number of constraints.
It is easy to see that Equation (28) is the common form of linear programming, which can be expressed as the following form by matrix form: where f , A eq and b eq are respectively as follows: where I is the identity matrix; O is the zero matrix. Equation (29) can be solved by calling the Linprog function in MATLAB, and the result is: Then the pseudo-state vector y constructed in the first stage can be obtained according to y = y (s) − y (t ) .

3.2.2
The second stage A non-linear change is made to the pseudo-quantity measurement y obtained in the first stage to obtain an estimated value y ′ . The specific transformation is as follows: It can be seen from this formula that although this stage is a nonlinear change process, it only need to do a little transformation of y to get y ′ , and the calculation process is very simple.

The third stage
For the estimated value y ′ obtained in the second stage, linear change is performed again. Set the intermediate variable x ′ at this stage as: Then the relationship between the estimated value y ′ and the intermediate variable x ′ is: y is the calculation error; S = diag(s i j ); A g , A e , and A h are the node-branch correlation matrices in the gas network,

FIGURE 2
The flowchart of the proposed state estimation algorithm electricity network, and heat network, respectively; I g , I e and I h are the unit matrix in the gas network, electricity network and heat network, respectively;Ā e is the node-branch correlation matrix with balanced nodes removed in the electricity network.
At this stage, the WLS method can be used to obtain the value of the intermediate variable x ′ , and then the fourth stage calculation is performed.

The fourth stage
Finally, the non-linear transformation method is used to convert the intermediate variable x ′ obtained in the third stage into state variable x.
In summary, the solution process of the IES state estimation method is shown in Figure 2.

CASE STUDY
This case study considers the electric-gas-heat IES shown in Figure 3. The system has a total of 40 nodes. The natural gas network consists of 12 nodes and 12 branches (1 compressor branch), among which the compressor is a gas turbine with constant outlet pressure. The electricity network consists of 14 nodes and 14 branches. The heat network consists of 14 nodes and 13 branches. The coupling element includes two CHPs and one P2G. Among them, CHP1 works in the mode of heat determined by electricity, and CHP2 works in the mode of electricity determined by heat. The specific structural parameters of the integrated energy system are listed in [23].
The following experiments are carried out using Matlab R2016b environment, with the experimental platform CPU 2.60 GHZ, RAM 4GB and Windows 10, a 64-bit operating system.

Effect analysis for SE
In order to verify the effectiveness of the proposed IES state estimation method considering electricity, gas and heat, the power flow calculation result of the IES obtained from the above system is used as the true value of measurement variables, the specific values are given in Table D1-D3. Consider using several indicators such as measurement error statistics, estimated error statistics and maximum estimated error of state variable to measure accuracy [24].
where z is the measured value of the measurement. h() is the measurement function. x true is the truth value of state variable.x is the is the estimated value of the state variable. is the standard deviation of Gaussian white noise. T is the number of Monte Carlo experiments.

Analysis for measurement redundancy
The method presented in this paper is called the robust state estimation method (RSE), and the measurement redundancy analysis of RSE and the WLS method are carried out separately. Among them, the number of measurements in the first stage of RSE is the same as the number of measurements in WLS.
The number of state variables in the third stage of RSE is also the same as the number of state variables of WLS. The number of measurements in the third stage of RSE is the same as the number of state variables in the first stage of RSE. The specific measurement redundancy of the two methods is shown in Table 1.

Accuracy test for state estimation
Superposing Gaussian white noise (the mean is 0 and the standard deviation is 0.001) on the basis of the true value as the measurement value of the system is considered to perform 1000 Monte Carlo simulation experiments. The maximum estimated error of state variables of each sub-network is calculated as shown in Table 2.
The state variable error of each sub-network is shown in Figures 4-6.
It can be seen from Table 2 that in terms of the maximum estimated error of the state variables, the state estimation results of the two methods are slightly biased, but overall there is little difference, and RSE is slightly better than WLS. In addition, it can be seen from Figures 4-6 that the overall estimation error of RSE method is small and tends to be flat. However, the overall estimation error of the traditional WLS method is large and the volatility is present. Therefore, the RSE proposed in this paper has better calculation accuracy than the WLS method. The statistical errors of the two state estimation methods are shown in Table 3.
The value of S H ∕S M represents the filtering effect. The smaller the value, the better the state estimation effect. Therefore, it can be seen from Table 3 that, compared with WLS, the RSE method proposed in this paper has better filtering effect for each sub-network, especially for the electricity network and the gas network.

Efficiency test for state estimation
Using the RSE method and WLS method to perform a single state estimation, the statistical results of the calculation time and the number of iterations are as follows: As can be seen from Table 4, the number of iterations in the single state estimation by the RSE method is much lower than the WLS method, and it takes less computational effort. In terms of the time for single state estimation, the calculation time of the two methods is not much different, and both have higher calculation efficiency. Therefore, the RSE method

Robust performance analysis for SE
Robustness refers to the characteristics of a control system that can maintain certain performance under a certain structure and size of parameter perturbation. The parameter perturbation is usually reflected in the bad data in the measured value. The bad data caused by equipment failure, environmental change or human factors in the measurement process has a serious impact on the state estimation. So it is necessary to detect and identify the bad data in state estimation. In this paper, the standard residual r N detection method is used to test the robust performance of RSE. This method greatly reduces the amount of calculation by arranging the measurement according to the standard residual absolute value |r N | from the largest to the smallest, and then conducting detection in turn according to the following criteria [24]. When the probability of misdetection is P e = 0.005, the detection criterion of the standard residual is:

Bad data on single subnet
This section considers the simulation experiment by randomly replacing the true value of the measurement with bad data. The form of bad data includes making the original data equal to 0, 1/2 of the original data, and −1/2 of the original data. For the integrated energy system, the electricity network has the fastest response speed and the measurement system is relatively complete, so the probability of bad data in the vector of measurement is low, while the measurement error of the heat network and natural gas network is relatively large. Therefore, this paper sets the following three types of bad data for the three subnets of electricity, gas and heat: 1. Suppose case 1 is that bad data only exists in the electricity sub-network, and the proportion of bad data in the measurement of the electricity network is 5.1%. In this case, the identification results of the two methods are shown in Table 5: The comparison of state variable estimation result in the electricity network is shown in the Figure 7. 1. Suppose case 2 is that bad data only exists in the natural gas sub-network, and the proportion of bad data in the measurement of the gas network is 8.3%. In this case, the identification results of the two methods are shown in Table 6: When there are bad data in the gas network, WLS+LNR method and RSE method proposed in this paper are respectively used for state estimation, and the comparison of pressure estimation result of nodes in the gas network is shown in Figure 8.  1. Suppose case 3 is that bad data only exists in the heat subnetwork, and the proportion of bad data in the measurement of the heat network is 5.9%. In this case, the identification results of the two methods are shown in Table 7: Results of the comparison of state variable estimation in the power network is shown in the Figure 9.
The maximum estimation error of the state variables of each sub-network in the three cases is shown in Table 8.
It can be seen that in case 1, the estimation error of electricity network using WLS method and RSE method is obviously      different, and the estimation accuracy of WLS method is significantly reduced. In case 2, the estimation error of the gas network using the WLS method has increased compared with before, but the RSE method proposed in this paper has not changed much after adding bad data. In case 3, the RSE method can still maintain good robustness, while the state variable accuracy of the heat network using the WLS method has declined to varying degrees.

Bad data on integrated network
The above three cases are the detection and identification of RSE method and WLS method when bad data occurs in a single subnetwork. When bad data occurs in each subnetwork at the same time, the detection and identification of the two methods are shown in Table 9: In the above situation, the maximum estimated error of the state variables of each sub-network is shown in Table 10. In case 4, in addition to analysing the changes in state variables in the IES, it is also necessary to pay attention to the changes in some key measurement variables. Such as the active and reactive power changes of nodes and branches, these variables are shown in Figure 13.
In summary, the following conclusions can be drawn: 1. When there is bad data in the electricity network, the RSE method has a better filtering effect than the WLS method. 2. When there is bad data in the natural gas network, both WLS method and RSE method can correctly find out the location of the bad data, but the RSE method exhibits better robustness. 3. When there is bad data in the heat network system, the estimation error of the RSE method is lower. 4. When there are bad data in the three subnets at the same time, if the WLS method is used for state estimation, the estimation accuracy of each subnet system is significantly reduced. However, the RSE method can overcome the influence of various types of bad data in the measurement data of each subnet system, obtain a state estimation result closer to the true value, and show good robustness and tolerance, which is better than WLS.

CONCLUSION
This paper combines the WLAV state estimation method with the WLS state estimation method and proposes a state estimation method for electricity, gas and heat coupling integrated energy system that provides a reference for the state estimation of an integrated energy system with multiple energy sources. The analysis shows that: 1. IES state estimation method proposed linearizes the nonlinear measurement model by introducing pseudo-state variables, without initial value selection and multiple iterations, and the calculation process is simple. 2. The state estimation method of IES has high calculation efficiency, the estimation accuracy can meet the requirements,  Figure A1.
In Figure A1, where m and n are pipe nodes respectively; i and o are the inlet and outlet nodes of the compressor respectively; p m , p i , p o and p n are the pressures at nodes m, i, o and n respectively; f mi is the flow of the pipeline into the compressor; f com is the flow through the compressor; f cp is the flow consumed to drive the compressor; f on is the flow of the compressor into the pipeline.
In the current study, the compressor usually works in the following four control modes: 1. The pressure at the inlet node of the compressor is constant: 2. The pressure at the outlet node of the compressor is constant: 3. The flow through the compressor is constant: 4. The compression ratio of the compressor is constant wherep i is the constant value of the pressure at the compressor inlet node;p o is the constant value of the pressure at the FIGURE A1 Model of gas turbine compressor compressor outlet node;f com is the constant value of the flow through the compressor;R k is the constant value of the compression ratio.

APPENDIX B: TOPOLOGY STRUCTURE OF HEAT NETWORK
The heat network can be divided into two networks with the same topology structure, which are the supply network and the return network. Through the supply network, hot water flows from the heat source to the heat load point. In this process, due to radiation and other effects, the hot water in the pipeline will dissipate heat outwards, so there will be heat loss from the heat source node 3 to the load nodes 1 and 2, which is similar to the network loss in the electricity network. At the load nodes 1 and 2, hot water provides heat to the user side through the heat exchanger, and then the cooled water is returned to the heat source through the return network for the next heating. Through the supply network and the return network, the heat is transferred and a water cycle is completed. The specific heat network structure is shown in Figure B1.
where I is the identity matrix; O is the zero matrix; * represents the relationship between the quantity measurement and the state variable within the same subnetwork, whose value is constant coefficient (not equal to 0); The remaining six modules represent the coupling relationship between the two subsystems. If there is a coupling relationship between the two sub-networks, the value is a constant coefficient, otherwise it is 0.