Locating Optimization of an Integrated Energy Supply Centre in a Typical New District Based on the Load Density

Abstract

In view of the obvious distinctions among energy utilization characteristics in various functional districts, a method for optimizing the location of an integrated energy system for a typical new district is proposed. This method integrates the spatial location and energy load density of each functional district. Based on the timing characteristics and seasonal division of electric/thermal/cold load fluctuations, the improved fuzzy C means (FCM) algorithm combined with the density-based spatial clustering of applications with noise (DBSCAN) is applied to build different multi-scenarios in various functional districts, and a more representative ‘typical maximum load value’ is obtained. The architectural characteristics and the different energy utilization habits are combed on the weight of the electric/thermal/cold load for each functional district. By taking the minimum improved integrated load moment of the system and the minimum partition inner interval ratio as the objective functions, the location optimization model of an integrated energy supply centre in a typical new district is established, and the improved adaptive evolutionary immune algorithm is applied to the nested model. Finally, a typical new district in Northern China is taking as an example to verify the correctness and feasibility of the planning method proposed in this paper.

1. Introduction

With the increasing demand for energy and the acceleration of urban construction, a new district acts as a gathering place for energy utilization. Such a district will have gathered many individuals with the following similar energy utilization characteristics: large energy demand, energy utilization concentration, and collected electric/thermal/cold multi-loads [1]. Moreover, in recent years, new districts have developed rapidly because of their suitable density and architectural scale and the full and reasonable utilization of resources. Because of the utility tunnels within new districts, a centralized energy supply mode is mostly used for electric/thermal/cold loads. Thus, the energy loss in the transfer process is reduced, and it is possible to render the electric/ thermal/cold multi-energy flow to be both coordinated and complementary. Typical new district projects include the districts of Hebei Xiongan and Hunan Xiangjiang in China, the Marae-la-vallee district in Brazil, the Milton Keynes district in Britain and the Marne Valley British district in France, and so on [2,3].

Considering the diversity of a new district’s load types and the relative closeness of the traditional electric/thermal/cold industry, the independent planning operation of the electric/ thermal/cold load is not conducive to improving the overall energy efficiency of a system. Multi-energy flow coordination planning is an effective strategy for energy efficient utilization [4,5], and the location of a reasonable integrated energy supply centre is particularly critical [6]. Some research on the optimal location of an integrated energy supply centre has been performed by domestic and foreign scholars [7,8,9,10,11,12,13,14,15]. In ref. [7], taking the key technology involved in multi-energy interconnection and integrated optimization as research objectives, through the mixture-modelling construction of the multi-energy flow, the coordinated planning of a multi- energy system was explored. In ref. [8], based on the distinctions between the seasonal characteristics and architectural types, cold/thermal/electric coordinated planning and energy management methods for different functional demand types were proposed. In ref. [9], the concavity, convexity and directivity of the weighted Voronoi graph algorithm was improved; furthermore, an optimized location method for electric vehicle charging stations based on an improved weighted Voronoi graph and a genetic algorithm was constructed. Considering the geographic information and administrative boundaries, an optimization method for substation supply areas was put forward in ref. [10], where the limitations of geographic information and administrative boundaries were dealt with by a penalty function method. In ref. [11], based on the geographic information system in a distribution network, a new method for active distribution network planning and the improvement of substations with the Voronoi diagram algorithm was proposed. In ref. [12], using the construction of a thermal network and the cold/thermal/electric cogeneration system, a collaborative planning model for multi-district integrated energy systems with heat networks was put forward. In ref. [13], from the perspective of the regional overall energy supply and demand, a two-layer programming model of a regional distributed energy system based on decision-making and the system configuration was constructed. An integrated energy system model for urban commercial residential mixed zones based on thermoelectric coupling was established in ref. [14], and the optimized location of the integrated energy centre of the system was selected. In ref. [15], taking the minimum weighted moment of the electricity/heat/cold as the objective and adopting the Voronoi diagram algorithm and the alternating location allocation algorithm as the criteria, a method that permits one to partition energy local area networks (LANs) and locate energy routers, is proposed.

Considering the volatility and uncertainty of the electric/thermal/cold load output of each functional district within a new district, it is necessary to consider the uncertainty influence of each type of load fluctuation when selecting the maximum load of each functional district [16,17]. Therefore, in ref. [18], the fuzzy membership method was used to describe the uncertainty factors which is made the load fluctuation as the triangular fuzzy variable. In ref. [19], under the environment of electric power market, the uncertainties of the load fluctuation were taken as random variables that obey a normal distribution, and the Monte Carlo simulations were conducted with consideration of the regularity of the load fluctuation in each functional district which is influenced by numerous factors, such as the environment and climate. What’s more, based on the timing characteristics analysis, the typical daily simulation method has been put forward in some researches. In ref. [20,21], the timing types are divided by seasons and weather, and a load multi-scenario construction method that depends on timing characteristics is proposed based on the uncertainty of the load fluctuation. According to the uncertainty of the load fluctuation in an active distribution network, an optimal scheduling method based on a multi-scenario was constructed in ref. [22]. A two-stage robust optimization model was applied in ref. [23], and multi-stage and multi-region uncertainties were established for the load fluctuation in the planning year. Therefore, the obtained multi-scenario set, which is constructed by means of the clustering method based on the timing characteristics, can enhance the volatility of the daily output curves of the load fluctuation when compared with the single typical daily curve simulation.

In summary, the current location optimization of an electric/thermal/cold energy supply centre is mostly contingent on the type of load that is needed and the location optimization of an electric/thermal/cold multi-load integrated energy supply centre is less involved. Moreover, the appropriate maximum load for each type of functional district is generally regarded as the annual maximum load value, and the uncertainty of the load fluctuation is considered less, especially in regard to the uncertainty of the thermal and cold load fluctuations. Therefore, by considering the construction of different multi-scenarios and the weight of the electric/thermal/cold load for each functional district, the integrated load density can be obtained. What’s more, by taking the minimum improved integrated load moment of the system and the minimum partition inner interval ratio as the objective functions, a comprehensive energy system location optimization method for a typical new district is proposed in this paper. This method integrates the spatial location and energy load density of each functional district. By examining a typical new district in Northern China, the correctness and feasibility of the proposed planning model is verified.

2. The Construction of an Electric/Thermal/Cold Load Model in a Typical New District Based on Timing-Phasing and Hybrid Clustering

2.1. Analysis of the Electric/Thermal/Cold Load Characteristics

For a typical new district, it is normal to divide the energy supply system into several functional districts based on the land classification plan, regional energy supply and load density. The main types of district include the residential district, commercial district, industrial district, and municipal district. The functional district refers to one or several users who have similar geographical positions with the same type of energy demand or a block of land with the same energy utilization properties, such as a factory, a residential district, or a commercial district composed of a number of shops. Therefore, the load classification of the new district mainly includes the residential load, the commercial load, the industrial load, and the municipal load [24].

In addition, according to the different energy utilization characteristics, the new district’s load can be divided into the electric load, the thermal load and the cold load, and the differences between the electric/thermal/cold load fluctuations are obvious for the functional districts with different land utilization types. For each functional district, the demand for the electric load is the largest, the thermal load is the next largest, and the cold load is relatively small. Moreover, there is a significant peak in the residential load, commercial load and municipal load, but the peak times and durations are different for the different functional districts. For the industrial load, there are small fluctuations and little difference in the load value throughout the day; that is, there are no obvious peak and valley periods. The typical timing characteristic curves of the electric/thermal/cold loads in various functional subdistricts of a certain new district are shown in Figure 1.

2.2. The Seasonal Characteristics of the Electric/Thermal/Cold Load

For the analysis of the electric/thermal/cold load characteristics in a typical new district, the differences among the electric/thermal/cold load characteristics of each functional district are considered (Section 2.1). For the thermal/cold load, the characteristics of the thermal load in the winter and the cold load in the summer are distinctive because of the existence of summer refrigeration equipment and winter heating equipment. Naturally, the thermal/cold loads are relatively lower in the spring and autumn due to the suitable weather conditions. For the electric load, due to the obvious differences in electric demand in different seasons, there is also a certain difference in its output characteristics. Therefore, there are obvious seasonal differences in the various load types. The typical timing characteristic curves of the electric/thermal/cold loads in different seasons in a typical new district are shown in Figure 2.

2.3. Multi-Scenario Construction of the Electric/Thermal/Cold Load in a Typical New District Based on Timing-Phasing and Hybrid Clustering

The fluctuation of the electric/thermal/cold load is easily influenced by meteorological factors, such as the environment, weather, temperature, humidity, rainfall and so on, which exhibit different patterns with the changing seasons, and the thermal and cold loads are relevant to the winter and summer, respectively. What’s more, considering the difference of the date types, such as weekdays, saturdays and holidays, the tendencies of the load fluctuation on different date types are different. Therefore, in the annual sequence, the four typical timing situations of the load fluctuation can be obtained due to the season difference. Then, the load fluctuation can be further divided into three typical timing situations for each season due to the date type difference. Thus, 12 typical timing situations can be obtained. Here, we will use a typical timing situation as an example. Using the timing-phasing of the load fluctuation for each functional district, the density-based spatial clustering of applications with noise (DBSCAN) [25] is introduced based on the traditional fuzzy C means (FCM) clustering algorithm [26] to complete the clustering of the electric/thermal/ cold load daily output curves and obtain the corresponding clustering centres. The daily output characteristic curves of the clustering centres are the scenarios set up by a certain type of load, and the scenarios’ number constructed by a certain type of load can be regared as the clustering centres’ number. Thus, the rough scenario division can be refined, and more accurate typical timing situations can be obtained. The process of clustering is as follows:

(1)

Initialize the membership matrix’s degree, U, determine the initial clusterings’ number, c, and the initial clustering centres’ number, v_i.

• For the set of clustering data, define the upper limit number of the initial clustering, l_max, and set the neighbourhood radius, Eps, and the density threshold, Minpts.
• Determine whether there are unmarked points in the data set. If no such points exist, the search ends; otherwise, randomly select an unmarked data point, d, to determine whether the data exceed the density threshold Minpts within the neighbourhood radius, Eps. If not, go to step d; otherwise, go to step c.
• Mark point d as a noise point that is no longer processed and return to step b.
• Mark point d as a core point, form the clusters and place the objects within the neighbourhood radius, Eps, into the clusters. Then, find the density reachable point in the cluster and add to the clusters until all of the data clustering is complete or the upper limit is reached. Return to step b.
• The clustering number, c, is the total number in the final cluster, and the initial clustering centre, v_i, is defined in Equation (1).

  $$v_i={\displaystyle \sum _{t=1}^m{P}_t}/m$$

  (1)

  where v_i (i = 1, 2, 3, ..., c) is the vector of the ith clustering centre, m represents the number of the ith class cluster, and P_t is the value of datum t.

(2)

Determine the weighting index, w, which can be regared as the final clustering effect’s fuzzy degree.

$$w^{\ast }=\left\{w|\frac{\partial }{\partial w}\left(\frac{\partial J_w\left(U,V\right)}{\partial w}\right)=0\right\}$$

(2)

where “|” is an explanation of the weighting index, w, namely, that there is no practical meaning; J_w(U, V) is the objective function. The distance which is weight the squares sum of each sample to all clustering centres is defined by Equation (3).

$$J_w\left(U,V\right)={\displaystyle \sum _{j=1}^n{\displaystyle \sum _{i=1}^c{u}_{ij}^w}}{\Vert x_j-v_i\Vert }^2$$

(3)

where v_i represents the ith vector of the clustering centre, and u_ij represents the membership degree of the jth sample of the ith clustering centre.

(3)

Update the clustering centre and the membership matrix degree, which are shown in Equation (4).

$$\{\begin{array}{ll}{v}_i=\frac{{\displaystyle \sum _{j=1}^n{u}_{ij}^w{x}_j}}{{\displaystyle \sum _{j=1}^n{u}_{ij}^w}},& 1\le i\le c\\ u_{ij}={\left[{{\displaystyle \sum _{k=1}^c\left(\frac{\Vert x_j-v_i\Vert }{\Vert x_j-v_k\Vert }\right)}}^{\frac{2}{w-1}}\right]}^{-1},& 1\le i\le c,1\le j\le n\end{array}$$

(4)

where x_i (i = 1, 2, 3, ..., c) is the ith load classification set.

(4)

Calculate the objective function through the Equation (3).

(5)

Calculate the iteration error of the two iterations of the objective function, ∆J_w(U, V), then determine whether it is less than the given positive number, $\epsilon$, and if this condition is not met, return to step 3. Otherwise, the clustering process is finished.

Therefore, taking a typical timing situation as an example, assuming that the obtained clustering centres’ number is k and the objects’ number contained in each cluster are n₁, n₂, ..., n_s, ..., n_k, the probability of the sth scenario can be given by Equation (5).

$$p_s=n_s/{\displaystyle \sum _{i=1}^k{n}_k}$$

(5)

If the number of the scenarios for the electric/thermal/cold loads are $N_{L,D}^l$, $N_{L,H}^l$, $N_{L,C}^l$ in a typical timing situation, l, of a certain functional district, the total number of scenarios in this typical timing situation is $N_l=N_{L,D}^l\times N_{L,H}^l\times N_{L,C}^l$. Similarly, the construction of a scenario based on timing-phasing and hybrid clustering can be obtained.

3. The Location Optimization Model of the Integrated Energy Supply Centre in a Typical New District

3.1. The Integrated Load Density of Each Functional District in a Typical New District

3.1.1. The Integrated Load Density

The load density is a quantitative parameter that characterizes the concentrated degree of the load distribution, and its value is closely related to the spatial location, population density and per capita income of the functional districts. For a functional district with multi-energy loads, the integrated load density is introduced to characterize the concentrated degree of the overall distribution of every type of energy load. Due to the detailed stipulation of land utilization types for each functional district in the city planning map, and given the different energy utilization characteristics, the weight ratios of the energy loads in each functional district are determined and analysed. Combined with the load density of each energy load, the integrated load density of the functional district is obtained with Equation (6).

$$D_L^i=w_{L,E}^i{D}_{L,E}^i+w_{L,H}^i{D}_{L,H}^i+w_{L,C}^i{D}_{L,C}^i$$

(6)

where $D_L^i$ is the integrated load density of the ith functional district; $D_{L,E}^i$, $D_{L,E}^i$ and $D_{L,C}^i$ are the respective electric/thermal/cold load densities of the ith functional district; and, $w_{L,E}^i$, $w_{L,H}^i$, and $w_{L,C}^i$ are the respective electric/thermal/cold load weights of the ith functional district.

3.1.2. The Integrated Load Density of Each Functional District in a Typical New District Based on the Construction of a Timing Multi-Scenario

For the calculation of the electric/thermal/cold load density in each functional district, after considering the construction of the timing multi-scenario of the electric/thermal/cold load, the maximum daily load of each scenario is selected. Then, through combination with the corresponding probability, a more representative ‘typical maximum load value’ is obtained. Assuming that the date number of each scenario in a typical timing situation of a certain energy load is L = {L₁, L₂, L₃, ..., L_n}, where n is the number of scenarios, the electric/thermal/ cold load density of the ith functional district can be obtained through combination with the covered area.

$$D_{L,Z}^i=\left({\displaystyle \sum _{i=1}^n\left(L_i/{\displaystyle \sum _{i=1}^n{L}_i}\right)}{\displaystyle \sum _{j=1}^{s_m}{p}_j{P}_{L,Z}^{j\max }}\right)/S_i$$

(7)

where s_m is the number of scenarios for each energy load; p_j represents the corresponding probability of the jth scenario; $P_{L,Z}^{jmax}$ is the maximum daily load curve value of the jth scenario; Z represents the energy utilization type and E, H, and C represent the electric/thermal/cold load types, respectively; finally, S_i is the covered area of the ith functional district.

The weight of the electric/thermal/cold load in each functional district is determined by its ratio relative to the electric/thermal/cold annual total load, which is used to measure the relative importance of the electric/thermal/cold load in each functional district. This weight’s value is related to the energy characteristics, energy use equipment, annual electric/thermal/cold supply days, and electric/thermal/cold transmission efficiencies of the various functional districts [27]. In addition, the weight considers the different architectural structure characteristics of the different functional districts—whereas a residential district mostly consists of middle-rise or low-rise buildings, an industrial district includes mostly low-rise buildings. Because of the limited urban planning area, commercial and municipal districts are mostly composed of high-rise buildings. For the different functional districts, due to the different suitable crowds, industries engaged and quality building materials, there are obvious differences in the buildings’ insulation effects, energy consumptions and enclosure structures in each functional district. Therefore, the architectural structure coefficient [28] is introduced to improve the electric/thermal/cold load weight calculation for each functional district.

$$w_{L,E}^i=\frac{\frac{T_{E,\max }^i{\epsilon }_E^i}{{\eta }^E}}{\frac{T_{E,\max }^i{\epsilon }_E^i}{{\eta }^E}+\frac{T_{H,\max }^i{\epsilon }_H^i}{{\eta }^H}+\frac{T_{C,\max }^i{\epsilon }_C^i}{{\eta }^C}}$$

(8)

where $T_{E,\max }^i$, $T_{H,\max }^i$ and $T_{C,\max }^i$ represent the maximum utilization hours of the electric/thermal/cold loads in functional district i, respectively; ${\epsilon }_E^i$, ${\epsilon }_H^i$ and ${\epsilon }_C^i$ are the building structure coefficients of the electric/thermal/cold loads in the ith functional district, respectively; lastly, η^E, η^H and η^C represent the respective energy transmission efficiencies of the electric/thermal/cold loads.

Similarly, we can obtain the thermal and cold load weights of the ith functional district, or $w_{L,H}^i$ and $w_{L,C}^i$, respectively. On this basis, we standardize each energy load density, and we obtain the integrated load density of the functional district:

$$D_L^i=w_{L,E}^i\frac{D_{L,E}^i}{D_{L,E}^{\max }}+w_{L,H}^i\frac{D_{L,H}^i}{D_{L,H}^{\max }}+w_{L,C}^i\frac{D_{L,C}^i}{D_{L,C}^{\max }}$$

(9)

where $D_{L,E}^{\max }$, $D_{L,H}^{\max }$ and $D_{L,C}^{\max }$ represent the maximum electric/thermal/cold load densities in all of the functional districts.

3.2. The Multi-Objective Location Optimization Model of the Integrated Energy Supply Centre in a Typical New District

3.2.1. The Objective Functions

a. The minimum improved integrated load moment of the energy supply system is

$$\min f_1={\displaystyle \sum _{i=1}^k{\displaystyle \sum _{j=1}^m{D}_L^j}}dis\left(s_i,s_j\right)$$

(10)

In the above equation, k is the number of integrated energy supply centres in a typical new district; m represents the number of functional districts of each energy supply centre partition; $D_L^j$ is the integrated load density of the jth functional district; lastly, dis(s_i, s_j) represents the space fusion distance between the ith and jth functional districts. The Euclidean distance is generally adopted as it is the traditional space fusion distance. Considering that the energy supply route is usually transported along the actual road, we use the Manhattan distance, which is closer to the actual route, to characterize the distance between two functional districts:

$$dis\left(s_i,s_j\right)=\left|x_i-x_j\right|+\left|y_i-y_j\right|$$

(11)

where x_i and y_i are the spatial location attributes of integrated energy supply centre I, and x_j and y_j are the spatial location attributes of integrated energy supply centre j.

b. The minimum inner interval ratio of the energy supply partition is

$$\min f_2=\frac{\frac{1}{k}{\displaystyle \sum _{i=1}^k{d}_i^{in}}}{\frac{1}{k(k-1)}{\displaystyle \sum _{i=1}^k{\displaystyle \sum _{l=1}^k{d}_{il}^{out}(i,l)}}}$$

(12)

where $d_i^{in}$ represents the inner distance for each energy supply centre partition, which characterizes the degree of similarity between each partition. The smaller the value is, the higher the degree of similarity will be. Furthermore, $d_{il}^{out}$ is the interval distance between two energy supply centre partitions, and it represents the degree of dissimilarity between the different partitions. Larger values of $d_{il}^{out}$ imply a greater degree dissimilarity. The calculation formula is Equation (13).

$$\{\begin{array}{l}{d}_i^{in}=\frac{1}{m\left(m-1\right)}{\displaystyle \sum _{p=1}^m{\displaystyle \sum _{q=1}^m\delta \left(z_i^p,z_i^q\right)}}\\ d^{out}\left(i,l\right)=\frac{1}{m_i{m}_l}{\displaystyle \sum _{p=1}^{m_i}{\displaystyle \sum _{q=1}^{m_l}\delta \left(z_i^p,z_l^q\right)}}\end{array}$$

(13)

where $\delta \left(z_i^p,z_i^q\right)$ is the spatial similar fusion distance between functional districts p and q in energy supply centre partition i.

$$\{\begin{array}{l}\delta \left(z_i^p,z_i^q\right)=w_d\times \Vert s_p-s_q\Vert +w_{\rho }\times \Vert {\rho }_p-{\rho }_q\Vert \\ \Vert s_p-s_q\Vert =dis\left(s_p,s_q\right)\\ \Vert {\rho }_p-{\rho }_q\Vert =D_L^p-D_L^q\end{array}$$

(14)

In Equation (14), w_d and w_ρ are the weights of the spatial location property and load density property, respectively. The load density value is 1, $\Vert s_p-s_q\Vert$ is the space fusion distance, which represents the spatial location property of the two functional districts, and $\Vert {\rho }_p-{\rho }_q\Vert$ is the integrated load density difference of the two functional districts, which represents the load density property of the two functional districts.

3.2.2. The Constraint Conditions

For the location optimization of an integrated energy supply centre in a typical new district, the constraint conditions include the maximum energy radius constraint of the energy supply centres, the numerical constraints of these centres, the space position coordinate constraints of these centres, the attribution constraint of each functional district and the spatial location constraints of each functional district.

$$\max \left\{dis\left(s_i,s_j\right)|i=1,2,\dots ,k,j=1,2,\dots ,m\right\}\le r_i^{\max }$$

(15)

$$2\le k\le 1+\lg n/\lg 2$$

(16)

$$\{\begin{array}{l}{x}_{\min }\le x_{i\in k}\le x_{\max }\\ y_{\min }\le y_{i\in k}\le y_{\max }\end{array}$$

(17)

$${\displaystyle \sum _{i=1}^k{\delta }_{ij}=1,j=1,2,\dots ,n}$$

(18)

$${\displaystyle \sum _{i=1}^m{\mu }_{ij}}>1$$

(19)

where $r_i^{\max }$ is the maximum energy radius of energy supply centre i, and n represents the number of functional districts in a typical new district. In addition, x_min, x_max, y_min and y_max are the maximum and minimum values of the coordinate position of energy supply centre i; these values are related to the scale of the planned typical new district. Furthermore, δ_ij is the attribution sign, which judges whether functional district j belongs to partition i—if it does belong, δ_ij = 1, but otherwise, δ_ij = 0, which ensures that only one energy supply centre partition corresponds to the result of each functional district partition. Similarly, µ_ij represents the position attribution sign of functional district i and functional district j—if two functional districts are connected, µ_ij = 1, but otherwise, µ_ij = 0, which ensures that any functional district that belongs to the same energy supply partition has at least one functional district in space connected to it.

3.2.3. The Model Solution

For the location optimization of an integrated energy supply centre in a typical new district in this paper, an improved adaptive evolutionary immune algorithm was used to solve the nested model. Evolutionary immunity is an algorithm based on the principles of the immune response. The main idea is to use the affinity of antibodies and antigens to describe the degree of approximation of the feasible solution and the optimal solution [29]. However, the traditional evolutionary immune algorithm has some shortcomings, such as a lack of diversity and poor local search ability. Therefore, in order to maintain the population diversity and get better global search ability, the adaptive genetic algorithm combined with the sigmoid function which contains a cosine transform is proposed. Then, considering the individual fitness value of the population, the crossover rate, p_c, and the mutation rate, p_m, can be obtained, and the model solution is shown as below:

(1)

Input the original planning data and initialize the population. Then, randomly generate N antibodies and extract m individuals from the memory library to form the initial population.

(2)

Calculate the antibody fitness value and the expected reproductive probability; then arrange the antibodies in descending order.

$$\{\begin{array}{l}P=\alpha \frac{\frac{1}{f_1+C\times f_C}}{{\displaystyle \sum _{i=1}^N\left(\frac{1}{f_1+C\times f_C}\right)}}+\left(1-\alpha \right)\frac{\frac{1}{N}{\displaystyle \sum _{j\in N}{S}_{v,s}}}{{\displaystyle \sum _{i=1}^N\frac{1}{N}{\displaystyle \sum _{j\in N}{S}_{v,s}}}}\\ S_{v,s}=\{\begin{array}{ll}1,& S_{v,s}>T\\ 0,& \mathrm{Other}\end{array}\end{array}$$

(20)

where P is the expected reproductive probability; α is a constant that represents the weight ratio of the affinity between antibody antigens and the concentration of antibodies; N is the number of antibodies; f₁ represents the system integrated load moment; C × f_C is the punishment for solutions that violate the distance constraint, where C is a large positive number; S_v,s represents the affinity between two antibodies; lastly, T is a predetermined threshold.

(3)

Form the parent population, and arrange the initial population in descending order according to the expected reproductive probability, P. Extract the first N individuals to form the parent population. The first m individuals form the memory library.

(4)

Determine whether the number of iterations has reached the upper limit. If not, the next step will be continued; otherwise, the iterative process is finished, and a non-inferior solution set will be output.

(5)

Implement adaptive selection, crossover and mutation through Equations (21) and (22). Then, return to step 3.

$$p_c=\{\begin{array}{l}\frac{p_{c2}+p_{c1}}{2}+\frac{p_{c2}-p_{c1}}{2}\cos \left[\frac{2\pi }{1+\exp \left(\beta \frac{f^{\prime }-f_{avg}}{f_{\max }-f_{avg}}\right)}\right],\hspace{1em}{f}^{\prime }\ge f_{avg}\\ p_{c1},\hspace{1em}{f}^{\prime }\le f_{avg}\end{array}$$

(21)

$$p_m=\{\begin{array}{l}\frac{p_{m2}+p_{m1}}{2}+\frac{p_{m2}-p_{m1}}{2}\cos \left[\frac{2\pi }{1+\exp \left(\beta \frac{f-f_{avg}}{f_{\max }-f_{avg}}\right)}\right],\hspace{1em}f\ge f_{avg}\\ p_{m1},\hspace{1em}f\le f_{avg}\end{array}$$

(22)

where α, β represent the constant parameters; p_m1, p_m2 are the minimum and maximum variance values, respectively; p_c1, p_c2 are the minimum and maximum crossover rate values, respectively; f_max, f_avg are the largest and average population fitness value, respectively; f’ represents the larger fitness value with two individuals to be crossed, and f represents the individual fitness value to be mutated.

(6)

Calculate the inner interval ratio in each partition result according to Equation (13). Then, compare the inner interval ratio’s value, and obtain the optimal partition results of the typical new district.

A flow chart of the integrated energy supply centre’s location optimization based on the improved adaptive evolutionary immune algorithm is shown in Figure 3.

4. Example Analysis

Taking a typical new district in Northern China as an example, the planned size of the new area is approximately 8.6 × 10⁷ m², and it includes eight residential districts, six commercial districts, seven industrial districts, five municipal districts, and four other functional districts (such as parks and greening sites). The functional type, spatial location and covered area of each functional district are given in Appendix A (

Table A1. Planning information for various functional districts in a typical new district.

Label	Functional District Type	Covered Area/×105 m2	Spatial Location		Label	Functional District Type	Covered Area/×105 m2	Spatial Location
			X Coordinate/m	Y Coordinate/m				X Coordinate/m	Y Coordinate/m
1	Residential district-1	22	2550	9100	16	Residential district-6	32.8	5600	5600
2	Municipal district-1	32	3450	9000	17	Residential district-7	41.6	7600	6350
3	Industrial district-1	48.8	4750	8800	18	Industrial district-3	15.6	8700	6800
4	Other district-1	17.2	1100	7850	19	Industrial district-4	33.2	4650	4800
5	Commercial district-1	27.2	2300	7150	20	Municipal district-4	34.8	7350	5100
6	Residential district-2	18.4	3550	7450	21	Commercial district-4	11.6	3800	3600
7	Commercial district-2	28.4	5100	7500	22	Commercial district-5	49.2	6350	3950
8	Residential district-3	30.4	6350	8100	23	Residential district-8	25.6	8700	4400
9	Industrial district-2	17.6	1600	6050	24	Industrial district-5	14.4	8150	3600
10	Municipal district-2	76.8	3050	4900	25	Other district-3	8.8	2600	3300
11	Residential district-4	23.6	3450	6250	26	Commercial district-6	25.6	3300	2700
12	Residential district-5	34.4	4800	6100	27	Municipal district-5	26.4	4250	2800
13	Commercial district-3	28.8	6200	6300	28	Other district-4	28.4	5150	3000
14	Other district-2	20.8	6700	7100	29	Industrial district-6	45.6	4900	1200
15	Municipal district-3	15.2	7850	8000	30	Industrial district-7	29.2	6150	2200

). This district’s planning map is shown in Figure 4.

4.1. Analysis of the Electric/Thermal/Cold Load Characteristics

Because the cold and thermal load mainly apply in summer and winter, respectively, the timing characteristics of the electric/thermal/cold load on weekdays in the summer and winter seasons in the typical functional districts are taken as an example in the analysis. The timing scenarios of the summer electric and cold loads on weekdays for the typical functional districts are shown in Figure 5 and Figure 6, respectively.

The timing scenarios of the winter electric and thermal load on weekdays for typical functional districts are shown in Figure 7 and Figure 8, respectively.

It can be seen from Figure 5, Figure 6, Figure 7 and Figure 8: based on the timing-phasing and hybrid clustering method, the multi-scenarios of the electric/thermal/cold load fluctuatuon can be obtained. For each typical timing characteristic curve of the electric/thermal/cold load fluctuation, its amplitude of variation are different, which represent several discrete scenarios, and its probabilities are also distinct, which indicate the likelihood of each scenario occurring. Therefore, the obtained multi-scenario curves have a larger coverage area and also can enhance its volatility when compared with the single typical daily curve simulation. Thus, the uncertainty of the electric/thermal/cold load fluctuation can be further characterized.

In addition, the timing scenarios graph is obviously different for different types of functional districts. There are two peak periods (early and late period) of the electric/thermal/cold load in the residential district, which are relatively lower in daytime. The peak period of the electric/ thermal/cold load in the commercial district occurs approximately between 6 p.m. and 10 p.m., and the night load value is almost 0. In contrast, the electric load in the industrial area is larger than the thermal and cold loads during the daytime and the peak valley’s value is not very different. However, for the municipal district, the peak period of the electric/thermal/cold load is roughly between 9 a.m. and 5 p.m., and there are two small peaks in the morning and afternoon. Based on the different date types and the seasonal characteristics, the number of scenarios constructed under the typical timing characteristics of each functional district is different. Taking the residential district as an example, the number of scenarios for the summer electric and cold loads on weekdays, weekends and holidays are 5, 4, and 4, and 3, 5, and 4, respectively; the number of scenarios of the winter electric and thermal loads on weekdays, weekends and holidays are 4, 4, and 5 and 4, 5, and 3, respectively. Therefore, for the residential district, 102 scenarios will be constructed. Similarly, 75, 94, and 86 scenarios are constructed for the commercial, industrial and municipal districts, respectively.

4.2. The Location Optimization of the Integrated Energy Supply Centre in a Typical New District

For the location optimization of the integrated energy supply centre in a typical new district, we first consider the different electric/thermal/cold load characteristics in different functional districts. Then, the weight ratios of the electric/thermal/cold loads in the typical functional districts of the new district are explored, as shown in

Table 1. The electric/thermal/cold load weight values in various typical functional districts.

Weight Value	Residential District	Commercial District	Industrial District	Municipal District	Other District
Electric load	0.54	0.41	0.63	0.51	0.44
Thermal load	0.24	0.33	0.21	0.28	0.32
Cold load	0.22	0.26	0.16	0.21	0.24

and Figure 9.

Based on the construction and corresponding probabilities of timing multi-scenarios, the maximum electric/thermal/cold load values for each scenario can be selected. Then, after taking into account the electric/thermal/cold load fluctuation in each typical functional district, the ‘typical maximum load value’ can be obtained. Combined with the weight ratio of the electric/thermal/cold load and the covered area of each typical functional district, the variation tendency of the electric/thermal/cold load density and the integrated load density for each functional district can be obtained. This tendency is shown in Figure 10, and the detailed data are given in Appendix A (

Table A2. Load density values in various functional districts.

Label	Electric Load Density/(w/m2)	Thermal Load Density/(w/m2)	Cold Load Density/(w/m2)	INTEGRATED Load Density/(w/m2)	Label	Electric Load Density/(w/m2)	Thermal Load Density/(w/m2)	Cold Load Density/(w/m2)	Integrated Load Density/(w/m2)
1	122.3	78.8	86.5	103.984	16	142.4	105.3	88.6	121.66
2	204.3	168.6	103.2	173.727	17	204.5	97.2	113.5	158.728
3	407.8	98.6	72.3	289.188	18	253.8	67.8	172.3	201.7
4	78.9	65.4	78.6	74.775	19	396.8	123.5	64.8	286.287
5	186.3	96.2	102.7	134.831	20	212.5	168.2	86.5	174.453
6	142.3	98.7	112.4	125.258	21	196.2	87.2	89.6	132.514
7	172.1	92.4	96.8	126.221	22	147.5	76.5	69.4	103.764
8	123.4	86.2	94.6	108.136	23	90.8	46.3	65.7	74.598
9	507.8	68.7	89.2	348.613	24	262.5	88.4	77.6	196.355
10	192.4	154.2	123.2	167.482	25	66.2	40.7	57.8	56.45
11	176.2	90.8	102.4	139.468	26	104.2	62.3	70.8	81.689
12	184.5	82.2	98.7	141.072	27	192.7	146.8	102.7	161.389
13	164.5	97.8	90.2	123.171	28	96.3	50.8	42.7	69.25
14	96.8	72.4	89.7	87.705	29	396.2	172.5	102.4	302.215
15	162.4	78.6	64.3	118.478	30	276.8	67.5	89.2	202.831

).

In Figure 10, the colour planes represent the load densities of the different energy utilization types for all of the functional districts. The electric load density of each functional district is larger, but the thermal and cold loads are relatively small, and the difference between them is slight. For the functional districts with the same type of energy utilization, the difference in the load density values of each functional district is obvious. For functional districts 3, 9, 18, 19, 24, 29 and 30, the load density values are significantly larger and the functional district type belongs to the industrial category; the load densities of municipal districts 2, 10, 20 and 27, commercial districts 5, 7, and 13 and residential districts 6, 11, 12, and 17 are relatively small, and the variation in the ranges of the load density of the commercial district, residential district and municipal district is very similar. The variation in range of the load density of the municipal districts is slightly higher than the corresponding variation in range for the residential districts, and the residential districts’ variation range is slightly higher than the commercial districts’ variation range. However, the load densities of functional districts 4, 14, 25 and 28 are the smallest, and their functional district types are parks, greening belts, squares and other areas with less obvious energy utilization characteristics. In addition, for functional districts with the same energy utilization type, due to the differences in geographical position, population density, economic level and energy utilization habits, the load densities of different functional districts may be different.

For the solution of the location optimization model of an integrated energy supply centre in a typical new district, the parameters are set as follows: the population size (N) is 70, the memory library capacity (m) is 30, the crossover rate (P_cross) is 0.5, the mutation rate (P_mut) is 0.4 and the diversity evaluation parameter (P_s) is 0.95.

The selection of the integrated energy supply centres’ partition numbers is closely related to the inner interval ratio of the energy supply partition, as shown in Figure 11.

As shown in Figure 11, the optimization number of each integrated energy centre corresponds to the inner interval ratio of the energy supply partition. When the optimization number of the integrated energy centre is 5, the inner interval ratio of the energy supply partition is the smallest. This result implies that the degree of similarity between each partition is higher and the degree of dissimilarity between the different partitions is greater; therefore, the partition result is optimal.

To further verify the correctness of the proposed planning method in this paper, we compared it with the adaptive robust optimization method that is mentioned in ref. [30,31,32]. The partition situation under the planning method based on the optimal energy supply centre number of the typical new district is shown in Figure 12. In addition, the coordinate positions of the integrated energy supply centres and the partition situation are given in

Table 2. The coordinate positions of the integrated energy supply centres and the partition situation.

Method	Partition Marking	X coor dinate/m	Y coor dinate/m	Functional District Marking	Sum of Integrated Load Density/(w/m2)	Iteration Number	Simulink Time/s
Planning method in this paper	Partition one	4550.7	8620.9	1, 2, 3, 7, 8	801.256	1258	160.23
	Partition two	2349.6	6663.4	4, 5, 6, 9, 11	822.945
	Partition three	4096.9	5046.6	10, 12, 16, 19, 21, 25	905.465
	Partition four	7579.9	5862.8	13, 14, 15, 17, 18, 20, 23, 24	1135.188
	Partition five	5177.2	2209.9	22, 26, 27, 28, 29, 30	921.138
Adaptive robust optimization method	Partition one	4476.6	8534.2	1, 2, 3, 6, 7, 8	926.514	1426	1960.86
	Partition two	2157.2	6537.4	4, 5, 9, 11	697.687
	Partition three	4446.7	4977.6	10, 12, 16, 19, 22	820.265
	Partition four	7579.9	5862.8	13, 14, 15, 17, 18, 20, 23, 24	1135.188
	Partition five	4864.9	2003.5	21,25, 26, 27, 28, 29, 30	1006.338

.

As shown in Figure 12 and Table 2, the two planning methods can both achieve location optimization of the typical new district using the optimal energy supply centre number. However, the planning result is slightly different. Specifically, there are differences in the number of functional districts for each partition, the attribution of individual functional districts and the coordinates of the integrated energy supply centre. Indeed, the variation in range of the total integrated load density of each partition for the planning method in this paper (333.932 w/m², which represents the equilibrium of the partition result) is obviously smaller than the result obtained using the adaptive robust optimization method (437.501 w/m²). Therefore, the equilibrium of the partition results for the planning method in this paper is superior to the results derived from the adaptive robust optimization method. In addition, during the process of optimization, considering the planning method proposed in this paper, the iteration number, 1258, and the simulation time, 160.23 s, are both smaller than the iteration number, 1426, and the simulation time, 196.86 s, of the adaptive robust optimization method; they are reduced by 13.5% and 22.86%, respectively. Thus, our method is more suitable for the location optimization of energy supply centres in typical new districts.

Figure 12a and Table 2, regarding the proposed planning method in this paper, we know the following: the number of functional districts and their energy utilization types may be different, and the sum of the integrated load densities is slightly different for functional districts of the same partition. For partitions one and two, the number of functional districts is the smallest, at 5, and partition four contains the largest number of functional districts, namely, 8. This characteristic is due to the obvious differences in the integrated load densities in each functional district. Moreover, the partition with fewer functional districts contains more functional districts with higher integrated load densities, such as functional district 3 in partition one and functional district 9 in partition two. Comparably, the partition with more functional districts contains more functional districts with lower integrated load densities, such as functional districts 14 and 24 in partition four. Thus, we can ensure the rationality and balance of each energy supply partition. We consider the different geographical positions of each functional district, as their distribution densities, energy supply types and energy supply radii are varied. As a result, the sum of the total integrated load density of each partition is slightly different. In addition, for the same partition, the coordinates of the integrated energy supply centre are not at the centre of each partition, but close to a functional district with a larger load density. For example, the coordinates of the integrated energy supply centre in partition one is close to functional districts 2 and 3, which have larger load densities, and the coordinates of the integrated energy supply centre in partition three is close to functional district 19. This pattern occurs because the integrated energy supply centres in different partitions give priority to the functional districts with larger load densities to ensure the optimization of the transmission path of the electric/thermal/cold energy in the process of energy transmission.

5. Conclusions

In this paper, for the location optimization of the integrated energy supply centre in a typical new district, based on the timing characteristics and seasonal division of the electric/thermal/cold load fluctuations, the improved FCM algorithm combined with DBSCAN density clustering was applied to build different multi-scenarios in different functional districts. A location optimization model of the integrated energy supply centre for a typical new district was proposed, which integrates the spatial location and energy load density of each functional district. Taking a typical new district in Northern China as an example, the following conclusions were obtained:

(1)

The timing multi-scenario characteristics of the electric/thermal/cold loads, the load density values and the load weight ratios for different functional districts are distinct.

(2)

The different optimized numbers of integrated energy centres correspond to the different inner interval ratios of energy supply partitions. The smaller the inner interval ratio value is, the better the partition situation will be.

(3)

The location optimization method with the minimum improved integrated load moment of the system and the minimum inner interval ratio of the partition can suitably balance the uncertainty of the electric/thermal/cold load fluctuation and the spatial location in each functional district. Therefore, the location of optimization of the integrated energy supply centre in a typical new district based on the load density has been realized.

In addition, compared with the traditional decentralized energy supply mode, the energy supply mode based on the coupling complementarity of electric/thermal/cold multi-energy flow can reduce the energy loss and achieve the high efficiency of electric/thermal/cold multi-energy complementarity. What is more, considering the coordination and mutual aid of the power grid and heating network, the mismatch problem between peak and valley times of electric/thermal/cold load demands can be effectively alleviated, and the fluctuation and imbalance in the power grid and heating network can be reduced, which may increase its stability. Therefore, the location optimization method proposed in this paper for the typical new district provides the direction for the further planning of energy configuration in each partition based on the electric/thermal/cold multi-energy complementarity and the coordinated operation of power grid and heating network.