Pipeline Network Layout Design of Integrated Energy System based on Energy Station Site Selection and Load Complementary Characteristics

In view of the lack of effective energy station site optimization method in the existing integrated energy system (IES) planning, and the failure to consider the load characteristics in the division process of the energy supply area of energy stations, a pipeline network layout method of integrated energy system is proposed based on energy station site selection and load complementary characteristics. Firstly, based on the analysis of the interaction between the station and network planning of IES, an alternated optimal method framework of energy station site, energy supply area division and pipe network layout are proposed. Secondly, based on the energy supply radius of the energy station, this paper proposes a method to determine the number of new energy stations, and based on the topological characteristics of the cellular networks, a method of determine the initial energy station site is proposed. Thirdly, considering the complementary load characteristics, the evaluation index and division method of energy supply area of integrated energy station are proposed. Fourthly, aiming at the problem of pipe network layout on the road, a pipe network optimization method based on Kruskal algorithm is proposed, and optimize the direction of load access to the pipe network by the linear optimization model. Fifthly, the calculation method of equivalent energy distance of load and the optimization method of energy station site based on equivalent energy distance are proposed. Finally, the rationality and effectiveness of this method are verified by an example.

C ES the annual investment cost of the energy station C pipeline the annual investment cost of the energy pipe network N ES the number of new energy station r 0 the discount rate m s the life cycle of the energy station P k the cooling/ heating load demand of the energy station k f (P k ) the investment cost of the energy station k, N v the number of pipes in the planning area c p,i the unit length cost of pipe i l i the length of pipe i c p,o the unit length construction cost of pipe The associate editor coordinating the review of this manuscript and approving it for publication was Qiu Ye Sun . m l the life cycle of pipe c p,loss the operation loss cost of pipe λ the thermal conductivity of pipe insulation material t c the average temperature difference between inside and outside the pipe T p the annual maximum working time of circulating water pump C e the unit electricity price δ i the thickness of pipe i δ ti the thickness of pipe i insulation layer COP the energy efficiency ratio of IES c p,m the maintenance cost of pipe µ the proportion coefficient of maintenance cost S the area of the planning area r min , r max the minimum / maximum energy supply radius of integrated energy station

I. INTRODUCTION A. MOTIVATION
With the rapid development of the global economy, the scale of energy utilization continues to expand, and energy crisis is becoming serious. Optimizing the energy structure has become an important strategy to solve the energy shortage in the world. Through multi-energy complementary mode, IES breaks the mode of independent planning and independent operation of various energy systems, such as cold, heat, electricity, gas, and becomes the main development direction of energy system. Accordingly, improving the flexibility of IES becomes the important objective for coupling among multiple energy carriers in [1]. Therefore, the construction of IES is of great significance to promote the optimization of energy structure in [2]. Energy supply and energy consumption of IES mainly includes energy station, energy distribution network and load in [3] and [4]. The energy station is composed of energy generation, conversion and storage equipment to provide energy supply for users. The purpose of energy station optimal configuration is to determine the construction capacity of various energy generation, conversion and storage equipment based on the energy demand in energy supply area in [5]. The energy supply area coordination of different energy stations can reduce the peak valley difference of the load in the energy station energy supply area, and then improve the utilization efficiency of the equipment, and effectively reduce the equipment investment cost of the energy station. The energy pipe network connects the energy station and load, which has complex topology and high investment cost in [6]. The pipe network layout is closely related to the energy station site and the load distribution, and the investment cost of the pipe network will be affected by the energy station site and the energy supply area division, so the pipe network layout must be considered in energy station siting. Therefore, station and network collaborative planning of IES can save energy system investment and improve energy system economy.

B. RELATED WORKS
At present, there are many references about the pipe network layout optimization of IES. Reference [7] establishes the cost model of energy pipe network, considers the impact of cold and heat load changes on the operation cost of the pipe network, and takes the minimum cost of pipe network construction as the goal, uses genetic algorithm to solve the pipe diameter optimization problem. In [8], the mathematical models of plane pipe network and space pipe network are constructed respectively, and the optimization methods of pipe network layout based on single parent binary genetic algorithm and pipe diameter based on integer genetic algorithm are proposed. In [9], the pipe network layout problem is transformed into the problem of flow distribution in the pipe route to be selected, and a pipe network layout method is proposed based on linear programming model. Combined with the impedance-based approach, the proposed stability assessment approach is utilized to analyze the stability of the microgrid with multiple converters in [10] and [11]. Compared with the existing characteristic equation methods, the proposed approach can be verified that the performance is significantly improved. In the above references, the load node is equivalent to the road node. Only the pipe network VOLUME 8, 2020 layout on the road is considered in the pipe network solution, but the access direction between the load in the block and pipe on the road of the planning area is not considered. In addition, some references proposed energy station access position optimization method based on the existing power and thermal networks. Different scenarios generated by load profile templates are innovatively integrated into the economic planning model to deal with uncertain operational states in [12]. In [13], considering the uncertainty of load growth, an integrated energy station planning model is established to solve the line expansion problem, pipe expansion problem and equipment access position problem. Reference [14] constructs a power flow model, a hydraulic system balance equation and a thermal system energy flow model, optimizes the access position of heat source and the operation strategy of heat source unit based on the existing power network and thermal network, aiming at the minimum loss of electric and thermal energy. These references studied the energy pipe network layout and the access position of energy equipment, but did not consider the interaction between the layout of energy station and the layout of energy network.
At present, there are few researches on station and network collaborative planning. Reference [15] presents a new multiple stages distribution expansion planning model where investments in distribution network assets, RES, ESS, and EV charging stations are jointly considered. Reference [16] proposes a collaborative planning method for single heat source, heat pipe network and access point between heat supply network and power network based on the existing power grid topology. However, it does not analyse the coordination problem among heating supply areas considering multiple heat sources. In [17], by taking load or existing station as control point, the Voronoi diagram is used to divide the planning area with or without existing stations. The Voronoi diagram intersection is taken as the site to be selected for integrated energy station, but the energy supply area division and the energy pipe network layout are not considered. In [3], the load energy distance is defined by the construction cost of pipe between energy station and load. Accordingly, the integrated energy system station and pipe network planning model is established, and the optimization method of energy station site and pipe network layout is proposed. In [18], the calculation method of electric load energy distance and heat load energy distance is proposed, and the method of energy station site and pipe network layout is proposed based on p-median. In [3] and [18], the energy station site selection methods must rely on the preset set of energy station site to be selected, which is difficult to apply to the planning scenario without preset, and the pipe sharing problem is not considered in pipe network layout. In [19], energy station site selection method and energy supply area division method are proposed based on K-means clustering, and pipe network optimization method is proposed based on a minimum spanning tree. However, the linear distance between load and station is taken as the clustering weight in [19], without considering the actual pipe route and the access direction between load and pipe on the road. Without considering the load characteristics, the energy supply area division methods, which are proposed in existing references, may cause the deviation of energy station equipment construction from the demand, which is unable to realize the economic construction of energy system.

C. CONTRIBUTION
In view of the above problems, the innovation of this paper are as follows. Firstly, based on the analysis of the interaction between energy system station and network planning, an alternated optimal method framework of energy station site, energy supply area division and pipe network layout are put forward to realize station and network collaborative planning. Secondly, based on the analysis of load complementary characteristics on the annual investment cost of the energy station, the evaluation index of the energy supply area division considering load complementary characteristics is constructed, and the energy supply area division method of the energy station is proposed. Through this method, the energy supply area division of the energy station is more reasonable, which can reduce the peak valley difference of the load within the energy supply area and realize the economic construction of the energy station. Thirdly, considering load complementary characteristics, a collaborative optimization method for pipe network and optimal access direction between load and pipe on the road is proposed based on graph theory and linear optimization model, which realized the integrated layout of energy pipe network. Fourthly, an energy station site optimization method is proposed based on the pipe network layout scheme and the load energy distance. This method considers the load distribution and the planning area road information, which is conducive to reducing the annual investment cost of pipe network.

D. ORGANIZATION
In summary, the energy supply area division of energy station will affect the pipe network layout scheme of IES. In the process of energy supply area division, it is necessary to consider the impact of division results on the economy of pipe network planning. Therefore, there is a close coupling relationship between energy supply area division and pipe network planning. At the same time, the energy supply area division result of the energy station is affected by energy station site selection, and the energy station site selection will also affect the energy pipe network layout result. There are mutual influences among the energy station site, the energy supply area division of energy station and the pipe network layout, so the above three need to be coordinated planning. Therefore, the station and network alternated optimal method framework of IES is proposed in Section II.
In order to solve this problem, the number range of new energy stations is determined in Section III firstly, according to the maximum and minimum energy supply radius of energy stations. Secondly, according to the different number scheme of new energy stations, each energy station initial site is determined and the energy supply range of each energy station is divided in Section IV. Thirdly, the pipe network layout of energy station is solved based on energy supply area division results of energy stations in Section V. Fourthly, combined with the pipe network layout scheme, the energy station site selection is optimized and the energy supply area of energy station is divided again in Section VI. The station and pipe network planning scheme of IES can be obtained in different number scheme of new energy stations, by repeatedly dividing the energy supply area of the energy station, optimizing the pipe network layout and the station site selection until the station site, the energy supply area and the pipeling network scheme convergence. Through the above process, this paper decouples the energy station site selection, the energy supply area division and the pipe network layout, and then considers the interaction among them through the alternative optimization of the above three, so as to realize the collaborative planning of energy station site, energy supply area and pipe network layout. This solution process mainly consults the substation planning solution idea, and also follows the solution sequence between the energy station site selection, energy supply area division and pipe network layout. That is the energy supply area division is based on the premise that the energy station site is known, and the pipe network layout is based on the premise that energy station site and energy supply area division are determined. Finally, an example is given to show the effectiveness of the proposed method in Section VII.

II. PIPELINE NETWORK PLANNING OF INTEGRATED ENERGY SYSTEM BASED ON ENERGY STATION SITE SELECTION A. PLANNING MATHEMATICAL MODEL
The energy station and network planning problem are to select a number of energy station sites in the planning area, determine the energy supply area of each energy station and the energy pipe network layout between energy station and load, aming at optimize construction and operation benefit for the station and network of IES. Therefore, the optimization objective function for the station and network of IES can be expressed as follow.
In formula, C ES represents the annual investment cost of the energy station. C pipeline represents the annual investment cost of the energy pipeline network. The expressions of each part can be listed as follows based on the life cycle theory.
In formula, N ES is the number of new energy station. r 0 is the discount rate. m s is the life cycle of the energy station. P k is the cooling/ heating load demand of the energy station k, MW. f (P k ) is the investment cost of the energy station k, which can be calculated according to the load within the energy supply area of the energy station, $. N v is the number of pipes in the planning area. c p,i is the unit length cost of pipe i, $. l i is the length of pipe i, m. c p,o is the unit length construction cost of pipe, which is related to pipe diameter d i , and function g(d i ) is used to express the relationship between them, $. m l is the life cycle of pipe, year. c p,loss is the operation loss cost of pipe, $. λ is the thermal conductivity of pipe insulation material, W/(m · • C). t c is the average temperature difference between inside and outside the pipe, • C. T p is the annual maximum working time of circulating water pump, h. C e is the unit electricity price, $/kWh. δ i is the thickness of pipe i, mm. δ ti is the thickness of pipe i insulation layer, mm. COP is the energy efficiency ratio of IES. c p,m is the maintenance cost of pipe, $. µ is the proportion coefficient of maintenance cost.

B. ALTERNATED OPTIMAL METHOD FRAMEWORK
According to collaborative planning mathematical model, the energy station annual investment cost is related to the load demand within energy supply area. The higher the peak load demand is, the larger the capacity of equipment to be built is, and the higher the annual investment cost is. In order to reduce the energy station annual investment cost, the energy supply area division method of the energy station should reduce the peak valley difference of the load supplied by the energy station as much as possible.
At the same time, according to the energy pipe network cost model, the pipe network annual investment cost is related to the pipe construction length and pipe diameter. For a single load, when the load is supplied by the nearest energy station, the construction length of the energy pipe from the energy station to the load reaches the minimum value, so the energy supply area division scheme of the energy station will be affected by the energy station site scheme. In addition, the pipe diameter of each pipeline is determined by the flow in each pipeline. Due to the common phenomenon of pipe sharing in energy piepe network, the pipe flow will depend on the load connected with the pipe. Therefore, the energy supply route between each load within the energy supply area and the energy station should be considered in the pipe network layout. When the energy supply area of each energy station is determined, the energy pipe network annual investment cost will be affected by the energy station site selection scheme.
Through the above analysis, it can be seen that after the energy station site determined, the reasonable energy supply area division of energy station can reduce the energy station annual investment cost and the pipe construction length between load and energy station. Considering the pipe sharing between different loads, the reasonable energy pipe network layout within the energy supply area of each energy station can realize the energy pipe network economic VOLUME 8, 2020 construction. The energy station site is optimized based on the load distribution and pipe routing information within the energy supply area, which can further reduce the cost of pipe network construction. The pipeline network planning method framework of integrated energy system based on energy station site selection is proposed, which is detailed as follows.
The station and network planning optimal method framework of IES mainly includes six subproblems, such as the number determining of new energy station, the new energy station initial site selection, the energy supply area division of energy station, energy pipe network layout, the equipment configuration of energy station and energy station site optimization. At present, many references have studied the equipment optimization configuration of energy station, which is usually solved separately from the other five subproblems. Therefore, this paper adopts the optimization model proposed in [5] to solve the equipment configuration of energy station, aiming at the objective of economic cost optimization. This paper focus on the other five subproblems and puts forward the corresponding solution methods, as described in the following chapters.
The planning model established in this paper is a mixed integer linear model, which is solved by CPLEX, a mature commercial solver in GAMS platform. The solver uses brance-and-cut algorithms to solve the mixed integer linear programming (MILP) problem, which can find the optimal solution in theory. The termination condition of convergence is shown in equation 8.
In formula, BF represents the objective function value of the current integer optimal solution. BP is the objective function value of the most possible integer solution. OptCR represents the convergence precision, which is set as 0.001 in this paper.
For the problem of ''How to ensure convergence''. According to the [20], the mixed integer linear model proposed in this paper is convergent based on the brance-and-cut algorithms in the CPLEX solver. Whether it is convergent for a certain example is related to the boundary conditions in the example. The CPLEX solver will give the solution status of the model in the output report. If it converges, it will give ''Optimal Solution'' or ''integer solution''. If it does not converge, it will give ''Infeasible Solution'', which means that the feasible solution under the convergence precision cannot be found.
At present, the problem is that the solution process of the optimal solution may last for a long time, and the solution efficiency is negatively related to the solution quality. Because the branch-and-bound tree may have as many as 2 n vertices, a problem with only 30 variables can generate a tree with more than 1 billion nodes. If no other termination criteria are set, the solution process may continue indefinitely until the search is completed or the computer memory is exhausted.
As the station and network planning of IES belongs to the planning problem, the off-line calculation is usually carried out before the construction of system, and the requirement for the algorithm solution speed is low. Therefore, the author thinks that the solution of this problem should improve the accuracy of the solution result as much as possible, and the solution time of the algorithm is within permit constraints. Therefore, the author sets time limit and convergence precision before optimization, and the optimal solution can meet the requirements of planning in accuracy and efficiency respectively. This is also the common way to solve large-scale mixed integer linear programming problem at present. In this paper, the time limit is set to 1800s, the convergence precision is set to 0.001, and the solution state of the model is ''Integer Solution'', which represents the optimal solution under the convergence precision.

III. METHODS FOR DETERMINING THE NUMBER OF NEW ENERGY STATION AND INITIAL STATION SITE A. NUMBER OF NEW ENERGY STATION DETERMINING
In the station and network planning of IES, it is necessary to determine the number of new energy stations N ES . For the central cooling / heating supply system, the larger the cooling / heating supply area and the greater the load density in the area, the more conducive to reduce the equipment investment cost in the planning area, and give full play to the scale advantage of central cooling / heating supply. However, with the increase of central cooling / heating supply radius, the pipe length also increases, and the cost of pipe and energy loss is also higher. Considering the construction and operation cost of the equipment and energy pipe network of the energy station, this paper sets the reasonable energy supply radius of the energy station as (500, 1000)m, such as in [21]- [24]. Accordingly, the number range of new energy station in the planning area (n min , n max ) can be determined, as shown in equations 9-10. n max = S/π r 2 min n min = S/π r 2 max (9) n min ≤ N ES ≤ n max (10) In formula, S is the area of the planning area. r min and r max represent the minimum / maximum energy supply radius of integrated energy station.

B. DETERMINATION OF NEW ENERGY STATION INITIAL SITE
According to the topographic characteristics of planning area, this part discuss the initial station site selection from two aspects, such as long and narrow area (LNA) and non long and narrow area (non-LNA) in [25]. Cellular network is formed by tight splicing of polygons and cover the area to be planned. Each station has the same coverage radius in the cellular network topology, and the number of station needed to completely cover the planning area is the least, which can realize the pipe network economic construction. This part proposes a method to determine the energy station initial site in LNA and non-LNA based on the cellular network distribution characteristics.

1) NON-LNA
In the traditional cellular network, each station coverage area is hexagon, which is connected with the surrounding six stations. However, triangle or quadrilateral splicing can also be considered in practical engineering applications in [26]. However, due to the inner angle limitation of the regular polygon, it is impossible to get the non overlapping close splicing scheme considering the restriction of any station number. Therefore, this part only analysis the cellular network site distribution considering the special number of new energy station.
Non-LNA is abstracted as a circle whose center is O and radius is R. Aiming at the restriction of differnt station number, the station distribution is mainly concentrated on the center O of the planning area and several concentric circles with O as the center. According to the cellular network site distribution law, it can be classified into four distribution modes: single ring, center + single ring, double ring and center + double ring.
In the single ring mode, each station is evenly distributed on the circle with O as the center and 1/2 R as the radius. In the center + single ring mode, each station is located at the center O and the circle with O as the center and 2/3 R as the radius. In the double ring mode, the stations are mainly distributed on two circles with O as the center and 1/4 R and 3/4 R as the radius. In center + double ring mode, each station is located at the center O and two circles with O as the center and 2/5 R and 4/5 R as the radius.
From the single ring mode to the double ring + center mode, the number of energy station increases gradually. Each station coverage radius is calculated, when supply area of staions completely cover the total area to be planned. Considering the moderate coverage radius of the station, conclusions can be shown as follows. When N station ∈ (2, 5), the single ring distribution mode is selected. When N station ∈ (6, 10), the center + single ring distribution mode is selected. When N station ∈ (11, 16), the double ring distribution mode is selected. When the number of new station continues to increase, the center + double ring distribution mode will be selected.

2) LNA
When the ratio of the planning area narrowest edge to the average radius is small, the area is defined as a LNA. In LNA, the stations of two neighboring lines adopt the interval distribution, and the crow flies between the adjacent stations is the same as far as possible to reduce the supply energy area overlap degree of the adjacent stations.

IV. ENERGY SUPPLY AREA DIVISION CONSIDERING COMPLEMENTARY LOAD CHARACTERISTICS A. ASSESSMENT INDEX
Energy supply area division scheme of energy station is evaluated from two indicators. One is the average peak valley difference index f pvd based on the load characteristics in the energy supply area of each energy station, which is used to judge whether the energy supply area division results can reduce the peak load and the equipment planned capacity, such as in [27]. The other is the distance index f ed between the load and the energy station, which is used to judge whether the energy network annual investment cost can be reduced within the energy supply area.
In formula, pvdp i,h and pvdp i,c represent the peak valley difference percentage of heating and cooling load in the energy supply area of energy station i. P i,h,max and P i,h,min represent the maximum and minimum heat load in the energy supply area of energy station i. P i,c,max and P i,c,min represent the maximum and minimum cooling load in the energy supply area of energy station i. pvdp i represents the peak valley difference percentage of multi-load in the energy supply area of energy station i. D is the total distance between each energy station and its load in the planning area. N load is the number of loads. d ij is the Euclidean distance between load j and energy station i. x ij is 0 / 1 variable, which indicates that load j is whether or not within the energy supply area of energy station i.
The range of f pvd is [0, 1]. When the peak valley difference of the load within the energy supply area of each energy station is 0, f pvd is 1, which means that the peak valley difference of the load within the supply range of each energy station is the smallest in this scheme. The range of f ed is [0, 1]. When all loads are supplied by the nearest energy station, f ed is 1, which indicates that the total distance between load and energy station is the shortest in this scheme. The comprehensive evaluation index F esa of energy supply area division scheme of energy station can be obtained based on f pvd and f ed , as shown in equation 17.
In formula, α pvd and β ed represent the weight of the two indicators respectively, and the sum of them is 1. The weights of the two indexes can be coordinated according to the geographic information and load information of different planning areas.

B. SOLUTION METHOD
The larger the load peak valley difference index in the energy supply area of energy station is, the better the effect of load complementary characteristics in energy supply area is. The peak energy consumption of different loads staggers each other, which can improve the equipment utilization efficiency. The larger the distance index of the energy station, the shorter the energy supply distance between the energy station and the load point, which can reduce the energy network annual investment cost. According to the peak valley difference index, a distance weighting factor σ is introduced to balance the above two indexes. The energy supply area division method of the energy station is proposed based on the weighted distance between the load and the energy station, as shown below.
Step 1, calculate the distance d ij from each load point to each energy station, and assign each load to the energy supply area of the nearest energy station.
Step 2, get the overall load characteristics within the energy supply area of each energy station, and calculate the distance weighting factor σ ed,ij between each load and each energy station according to Eqs. 18-20.
In formula, σ ed,ij,c and σ ed,ij,h respectively represent the weighting factor between load j and energy station i only considering the characteristics of cooling load or heating load. When load j does not belong to the energy supply area of energy station i, ω ij,c and ω ij,h respectively represent the overall cooling / heating load peak valley difference within the range after the load j is classified into the energy supply area of energy station i. ω kj,c and ω kj,h respectively represent the overall cooling / heating load peak valley difference within the energy supply area of energy station k. γ is the weighting factor amplification factor.
Step 3, calculate the weighted distance σ ed,ij × d ij between each load point and each energy station. All loads are classified into the energy supply area of the energy station with the minimum weighted distance.
Step 4, repeat step 2 and step 3 until the energy supply area of each energy station no longer changes.

V. PIPE NETWORK PLANNING CONSIDERING LOAD COMPLEMENTARY CHARACTERISTICS
The energy pipe network planning optimize the pipe route to minimize the energy pipe network annual investment cost based on the energy station site scheme and the energy supply area division result. The objective function is shown in equation 21.
min C pipeline s.t. Q T load = q T v m T (21) In formula, Q T load is the load flow vector. q T v is the pipe flow vector. m T is the pipe circuit matrix.
Since the load point is usually located in the buildings surrounded by multiple roads, and the pipe network must be arranged in the road corridor. Therefore, the optimal access direction between load and road should be optimized together with the pipe network in the energy pipe network layout optimization in [3]. There are multiple access directions between each load and road pipe network, which will further increase the pipe network layout complexity. In order to simplify the problem solution difficulty, a pipe network layout method is proposed based on Kruskal algorithm, and then a load access direction optimization method is proposed based on linear model. Finally, particle swarm optimization algorithm is used to optimize the road pipe network and load access direction.

A. PIPE NETWORK LAYOUT METHOD BASED ON KRUSKAL ALGORITHM
According to the planning area road information, the nodes are divided into road node C, load node L, access node T of load and road. In order to meet the pipe laying along the road constraint, the connection scheme between road node is formed only concidering road node C in [28]. There are multiple access directions between the energy station and the energy pipe network. All access nodes between the energy station and the road, both end points of this road is treated as internally connected root node. On this basis, a spanning tree composed of road node C is constructed to meet the radiant energy pipe network topological constraint.
Undirected graph G (V , E, W ) can be constructed based on the road information. Point set V = {v 1 , v 2 , . . . , v n } represents the road nodes in the planning area, edge set E = {e(v i , v j )|v i , v j ∈ V } represents the connection relationship between the road nodes in the planning area, and weight set W = {w(v i , v j )|e(v i , v j ) ∈ E} represents each road weight. Energy station is generally located in the central of the energy supply area. In order to ensure the energy pipe network economic construction, the pipe should be laid on the road within the energy supply area as far as possible, circuitous energy supply from the energy station to the load should be avoided, increase the pipe sharing between different loads to reduce the length of pipe construction. For this purpose, the weight of the road within the energy supply area is set as the road length, the weight of the outside road is twice the road length. The probability of internal road selection can be increased based on the road weight. Considering the irregular energy supply area, outside road may also be constructed.
Kruskal algorithm is used to build the spanning tree based on road information and road weight of the planning area,. The specific steps are as follows.
Step 1, take the road within the energy supply area of the energy station as the edge to be selected, build the edge set E and the weight set W for each road.
Step 2, regard the end points at both ends of the road where all access nodes between the energy station and this road are located as the internal connected root nodes, and delete this road from the edge set E.
Step 3, sort the roads according to the weight, select the line with the least weight in turn, judge it will form a loop whether or not with the selected line and internal connected lines of root nodes. If so, give up the line, otherwise add the line to the selected line set.
Step 4, repeat step 3 until selected lines and internal connected lines of root nodes constitute the spanning tree T 0 .
Step 5, judge each load has access point whether or not in the spanning tree T 0 . if not, select a road adjacent to the load to access the pipe network, and then break the loop formed. Repeat step.5 until there are access points between each load and the pipe network, and record it as the spanning tree T r .

B. LOAD ACCESS DIRECTION OPTIMIZATION METHOD
After obtaining the spanning tree T r composed of the the road point C, the corridor information of load point access location can be discussed from three types. There is only one access pipe around the first type of load, and only one branch access pipe around the second type of load. The remaining loads are classified as the third type. Among them, the first type of load will be directly connected to the pipe network. The second type of load can also be directly connected to the branch pipe to avoid redundant pipe. The third type of load access problem is more complex. Due to the phenomenon of pipe sharing, the load access pipe selection will affect the annual investment cost of multiple pipes. It is difficult to directly judge the advantages and disadvantages of different access schemes for the situation that the pipe to be selected contains multiple branches or multiple leaf branches. The pipe flow is determined by all loads supplied by this pipe, so all load access schemes need to be coordinated planning. Therefore, this paper solves the third kind of load access direction problem based on the pipe network linear optimization model proposed in [29], as shown below.
Step 1, increase all access pipes between the third type of load and the pipe network, and transform the pipe network optimization objective function. The annual investment cost of each pipe is expressed as the product of pipe flow and unit flow annual investment cost, as shown in equation 22.
In fomula,c pipe,i represents the pipe i unit flow annual investment cost. q vi represents the pipe i flow.
Step 2, initialize each pipe unit flow annual investment cost c pipe,i , solve the linear model (21) to get each pipe flow q vi based on each node average loadQ load . q vi > 0 means that pipe i will be constructed, and q vi = 0 means that pipe i will not be constructed.
Step 3, according to the pipe planning result, the node correlation matrix A is constructed to calculate each pipe sequential pipe flow.
Step 4, determine the pipe diameter according to the maximum pipe flow, and update each pipe unit flow annual investment costc pipe,i .
Step 5, judge the pipe network annual investment cost converges whether or not. If so, output the pipe network layout scheme. If not, return to step.3.

C. PIPE NETWORK OPTIMIZATION METHOD BASED ON PARTICLE SWARM OPTIMIZATION
Combined with optimization methods of pipe network layout and load access direction, cold and hot pipe network layout is optimized based on particle swarm optimization (PSO).
According to the pipe network layout method, a basic road pipe network layout scheme T 0 can be obtained. Each unselected road needs to be added in the edge set S to be selected. Several unselected lines is selected to join T 0 and form multiple circuits. For each circuit, a pipe with the largest weight except the newly added pipe can be removed to generate a new road network layout scheme T i . If the number of unselected roads is assumed to be N l , the position and velocity of particles is N l -dimensional vector. Each encoding of the particle is a 0 / 1 variable, indicating this pipe is added whether or not in the pipe network layout.
After decoding according to the particle code, it is necessary to check there are access points for each load whether or not in the pipe network layout scheme. If not, it is necessary to be corrected according to step.5 in Part A, Section V. Based on the pipe network layout scheme, the method proposed in Part B, Section V is used to optimize the load access scheme, and each particle fitness is calculated based on the pipe network layout and load access scheme.

VI. ENERGY STATION SITE OPTIMIZATION METHOD BASED ON LOAD ENERGY DISTANCE
Considering the energy supply network annual investment cost from the energy station to each load, the station site optimization objective function is shown in equation 23.
In formula, D S,i represents the investment cost of the energy supply network between the energy station and the load, that is, the energy distance between the load and the energy station, such as in [3] and [18].
Due to the common phenomenon of pipe sharing between loads, the investment cost of pipe between the energy station and a load can not be all included in the energy distance of this load. In this paper, according to the average pipe flow required by each load, the pipe investment cost is apportioned in proportion. According to this, the calculation of each load energy distance is shown in equation 24.
In formula, G(S, i) refers to the pipe set between energy station S and load point i. c k refers to pipe k unit length cost. l k refers to pipe k length.Q i refers to the average pipe flow required by load i.q k refers to pipe k average flow.
Combining the load energy distance and the Euclidean distance between the load and the energy station, the unit energy distance D t S,i between each load and the energy station in pipe network layout scheme can be obtained, as shown in equation 25.
In formula, d iS represents the Euclidean distance between load point i and energy station S.
In the traditional multi-source continuous site selection problem, the principle of minimum load moment is often used to optimize the station site in [29] and [30]. But this method is difficult to reflect the impact of the actual pipe route on the station site optimization. Therefore, an energy station site optimization method is proposed based on the principle of minimum load energy distance, as shown in equation 26. Through the iterative process, the energy station site optimization can be finally obtained.
In formula, x S and y S represent the abscissa value and ordinate value of the energy station site S respectively. x i and y i represent the abscissa value and ordinate value of the load point i respectively.

VII. EXAMPLE ANALYSIS
Tianjin eco city is a strategic cooperation project between the governments of China and Singapore. The construction of eco city provides positive discussion and typical demonstration for the construction of resource-saving and environment-friendly society.
Eco city is located in Binhai New Area of Tianjin, which is 40 km away from the center of Tianjin. In order to meet the needs of China's urbanization development, the eco city covering an area of 30 km 2 will be built into a sustainable and harmonious urban community based on the model of new towns in Singapore and other developed countries.
Two areas are selected from Tianjin eco city as the research object for data collection and planning in the example. The intial data and user information of the example are obtained from three system platforms, i.e. Measurement Automation System, voltage monitoring system and Power Quality Data Center. Other information can be obtained directly from the system platform and working report of state grid Tianjin electric power company.
The original map of the example area and its latitude and longitude coordinates are collected from the internationally recognized open source map system OpenStreetMap, and obtained through the secondary editing of josm map editing software, as shown in Fig.1 and Fig.6. The color nodes in Fig.1 and Fig.6 represent all the nodes in the load area. The symbols L1 and so on indicate the possible pipeline layout direction around the load point, which is determined by the way perpendicular to the surrounding road. At this time, the access pipeline is the shortest. The symbols n1 and so on represent the start node and the end node of the pipeline. The comparison scheme brief description of the example is shown in TABLE 1.
This paper uses matlab programming language to call cplex commercial solver on yalmip platform. Yalmip is a free optimization solution tool, which is an ''integrator''. It not only contains basic linear programming solution algorithms, such as linpro (linear programming), bintprog (binary linear programming), bnb (branch and bound), etc., but also provides a higher-level package for cplex, glpk, lpsolve and other solvers. Yalmip realizes the separation of modeling and algorithm, and provides a unified and simple modeling language, which can be used for all planning problems to reduce the learning cost of learners. The cplex is developed by ILOG company, which is the world's top software package for solving linear programming, integer programming and some nonlinear programming. The cplex is very fast, can solve many large-scale problems with millions of constraints and variables, and has been refreshing the highest performance record of mathematical programming. Because yalmip and cplex belong to the software package with optimization algorithm, the programmers do not need to program the optimization algorithm used. They only need to set reasonable control variables, objective functions, constraints and call appropriate solvers to solve the problem, which greatly simplifies the workload of the programmers. Fig.1 shows the planning area AS 1 used in this example analysis, including 24 communities. Circle represents residential load, diamond represents commercial load and rectangle represents administrative load. The planning area is 3.2km 2 , and an energy station will be built for central cooling / heating.

A. SINGLE ENERGY STATION SITE AND PIPE NETWORK LAYOUT COLLABORATIVE OPTIMIZATION 1) EXAMPLE INTRODUCTION
The cooling and heating loads in the planning area are 1.84MW and 1.39MW respectively, the informations of each load are shown in Fig.2.
The planning area mainly includes three types of loads, i.e. administrative, commercial and residential loads. The cooling and heating demand for different load types is shown in TABLE 2.
The temperature difference between the supply and return water of the pipe is 45 • C. The available standard pipe types and other relevant parameters are shown in Table 3, such as in [6] and [31]. The pipe service life is 20 years, the discount rate is 0.08, and the proportion of operation and maintenance cost is 0.25. The electricity price is 0.086$/kWh, the annual  operation hours of water pump is 3360h, the thermal conductivity of the pipe insulation material is 0.06W/(m× • C), the average temperature difference between the inside and outside of pipe is 10 • C, and the energy efficiency ratio of the integrated energy system is 4.2 in [5].

2) ANALYSIS OF RESULTS a: STATION AND NETWORK PLANNING RESULTS
First of all, station and network planning is based on different energy station site optimization schemes without considering load characteristics. Scheme 1, take the minimum load moment as the target to determine the energy station site, and then carry out the pipe network layout. Scheme 2, take the minimum load energy distance as the target to carry out the station and network planning.
The both scheme planning results are shown in Fig.3. The energy station is represented by black-box and marked with ES. Loads of the planning area are supplied by a single energy station in Fig.3, so both energy station investment costs are same, and the energy pipe network planning result of each scheme is only analyzed as follow.
There is a very common phenomenon of pipe sharing between different loads in Fig.3. The pipe actual construction length accounts for only 35% ∼ 40% of the total path length between all loads and the energy station in both schemes, as shown in TABLE 4. Therefore, considering the pipe sharing between loads, it can improve the pipe network utilization rate, reduce the pipe length of the energy pipe network, and realize the economic construction of the energy pipe network. According to the principle of minimum load moment, the site coordinate is (54.93, 109.59), and the energy station is located in block 14. In the station site optimization method based on the load energy distance, due to the twists and turns of the energy supply pipeline at the lower right of the planning  area, the load energy distance in the lower right area is higher. Therefore, in the station site optimization method based on the load energy distance, the station site will move to the load point with the higher load energy distance, and the station site coordinate is (56.81, 108.00), where the energy station is located at block 15 as shown in Fig.3 (b).
Compared with scheme 1 in TABLE 4, the pipe construction length in scheme 2 is reduced by 4.3%, and the energy system pipe network annual investment cost is reduced by 920.63$. The above shows that the load energy distance reflects the actual pipeline construction cost of the unit Euclidean distance between the energy station and each load node, the sum of the load energy distances between each load point and the energy station is reduced, which can reduce the construction cost of the actual pipe network. The station site scheme based on the load energy distance is more conducive to improving the economy of energy pipe network construction.
Normalization result of load energy distance and load demand in each load node is shown in Fig.4.

b: ENERGY NETWORK PLANNING CONSIDERING LOAD COMPLEMENTARY CHARACTERISTICS
In order to analyze the influence of load characteristics on the energy pipe network planning, each pipe diameter will be determined according to the energy flow peak value of the pipe. Scheme 3 of energy station and network planning can be obtained in Fig.5.
Through the comparison of Fig.3 and Fig.5, it can be seen that the selected pipe types of some road are smaller consider-  ing the load characteristics complementation. Scheme 3 can shorten the length of pipe 133 and pipe 159 significantly and increase the length of pipe 76 ∼ pipe 108 correspondingly, as shown in TABLE 5.
Compared with scheme 1 and scheme 2, the overall pipe network annual investment cost in scheme 3 is relatively low and the energy pipe network annual investment cost is reduced by 1987.52$ and 1066.90$, respectively. At the same time, with the change of pipe diameter, the average utilization rate of energy pipe network increased from 31.11% and 31.48% to 33.29%.
As there are dozens of pipelines in the planning area, the time series pipe flow of each pipeline can be solved based on the pipe network layout scheme and the time series demand of each load point. Considering the length limitation of the text, the author did not give the time series pipe flow of each pipeline in the paper, the following five pipelines are selected to show the time series flow in the pipeline, as shown in TABLE 6.
In order to more intuitively show the influence of load characteristics complementation on pipe diameter, TABLE 6 shows the pipe flow in S-T28 at the energy station outlet, the peak sum of all loads supplied by this pipe is 41.25t/h. The maximum pipe flow in pipe S-T28 occurs at 18, accounting for only 71.15% of the peak sum of all loads supplied by this pipe. It can be seen that the pipe diameter selection  Fig.6 shows the planning area AS 2 used in this example analysis, including 68 communities. The planning area is 4.71km 2 . The planning area mainly includes three types of loads, i.e. administrative, commercial and residential loads. The cooling and heating loads in the planning area AS 2 are 3.8MW and 3.2MW respectively. Each load information is shown in Fig.7.
Devices to be selected in the energy station include gas internal combustion engine, absorption chiller, conventional chiller, dual condition chiller, gas boiler, ice storage device, cold water storage device and heat storage device. Parameters of each device is shown in TABLE 7.

2) ANALYSIS OF RESULTS
Considering the area of the planning area and the energy supply radius of the energy station, 2-6 energy stations will be built in the planning area AS 2 . Station and network planning is based on the following two strategies.  Scheme 4, take the minimum load energy distance as the target to carry out the station and network planning. Weighting factor σ ed,ij is used to considering load complementary characteristics in the energy supply area division of energy station.
Scheme 5, take the minimum load moment as the target to determine the energy station site, and then carry out the pipe network layout, without considering load characteristics complementation in the energy supply area division of energy station.
Based on schemes of each energy station site and energy supply area division, devices combination of each energy station can be obtained aiming at the lowest construction and operation cost of energy system, and then each energy station annual investment cost in each scheme can be calculated. The energy station and network planning result of each scheme is shown in Fig.8, and the annual investment cost of each scheme is shown in TABLE 8 and TABLE 9.
According to the comparison of TABLE 8 and TABLE 9, a part of load is no longer supplied by the nearest energy station considering load complementary characteristics, but belong to the relatively far energy station. Therefore, the annual investment cost of energy pipe network in   4 is increased by 1563.06$ compared with scheme 5. However, considering load complementary characteristics, the overall load characteristics within the energy supply area of each energy station have been improved. On the one hand, the operation loss cost per unit pipe length has been reduced, and the pipe operation cost has slightly decreased compared with that before. On the other hand, the peak valley  difference of cooling and heating load of each station has been reduced from 94.80% and 99.79% to 88.02% and 99.42%. Therefore, the peak value of cooling and heating load of each energy station decreased by 14.17% and 12.13% respectively, and investment cost of devices in each energy station also decreased. Compared with scheme 5, the total station and network annual investment cost of scheme 4 is reduced by 14985.3$.
It can be seen that considering the load complementary characteristics, part of the load may be supplied by a relatively far energy station, which will increase the construction cost and operation cost of the energy pipe network. However, considering the load complementary characteristics can reduce the overall load peak valley difference within the energy supply area of each energy station, reduce the planned capacity of energy station equipment, optimize the pipe diameter type selection, and improve the overall planning scheme economy of IES. In order to consider all the costs comprehensively, the annual construction cost and operation loss cost of the pipe network are taken into account in the calculation of the annual investment cost, and the operation loss cost of the pipe network is shown in equation 6. Therefore, the operation loss has been included in the annual investment cost of pipe network in TABLE 8 and TABLE 9.
In addition, since the construction cost of pipe network and energy station equipment is usually one-off expenditure in the initial period of planning, while the operation cost is annual expenditure. In order to comprehensively consider the construction cost and operation cost, the equal annual value is used in the economic measurement of IES planning scheme, combined with the whole life cycle of pipe network and equipment. It can be seen from TABLE 8 that when the total cooling and heating load is about 7MW, the energy station annual investment cost of IES is 17.92 × 10 4 $. Compared with Scheme 5 in TABLE 9, the energy station annual investment cost is reduced by 8.45% considering the load complementary characteristics. Compared with the overall economy of energy network and energy station equipment of Scheme 5 in TABLE 9, the annual total investment cost of station and network collaborative planning scheme is reduced by 4.7% considering the load complementary characteristics. Therefore, the planning results based on the method proposed in this paper have certain reference value.

VIII. CONCLUSIONS
This paper constructs the mathematical model for station and network planning of IES, and proposes the pipeline network layout method of integrated energy system is proposed based on energy station site selection and load complementary characteristics. The effectiveness of model and method is verified by an example, some conclusions are obtained as follows.
Firstly, based on load energy distance, the station site optimization method can not only consider the load distribution, but also reflect the routing information of pipe network. The station site selection scheme obtained from this method can reduce the pipe construction length and improve the economy of pipe network construction.
Sencondly, the pipe diameter type selection can be optimized in the pipe network layout considering load complementary characteristics, which is conductive to increase the pipe network, utilization efficiency, and improve the economy of pipe network layout scheme.
Finally, the energy supply area division should be realized based on the load complementary characteristics, which is conductive to effectively reduce the peak valley difference of the overall load within the energy supply area of the energy station, reduce the peak load supplied by the energy station, reduce the investment cost of devices in the energy station, and improve the economy of station and network planning scheme of IES.
In the future, the transmission capacity of feeder lines and other factors can also be considered to realize the pipeline network layout method of integrated energy system is proposed based on energy station site selection and load complementary characteristics, and more detailed research will be carried out on methods to deal with the sudden change or slowly varying change in load.