Location and Capacity Selection Method for Electric Thermal Storage Heating Equipment Connected to Distribution Network Considering Load Characteristics and Power Quality Management

: High-permeability distributed wind power and photovoltaic systems are connected to the distribution network, which exacerbates the volatility and uncertainty of the distribution network. Furthermore, with the increasing demand of heating in winter and environmental protection, the wide use of electric thermal storage heating equipment (ETSHE) can promote distributed renewable energy utilization. However, an unplanned ETSHE connection to the distribution network may cause serious power quality problems. A new method of equipment location and capacity is proposed, which considered the improvement of power quality and load demand characteristics of the distribution network. First, based on heat load portrait technology, the node’s thermal load classiﬁcation prediction was carried out to provide the data basis for the model solution. Second, the multi-objective optimal location and capacity programming model including harmonic distortion rate, voltage deviation, voltage ﬂuctuation, and ETSHE cost was established. Then, the system nodes were preprocessed based on the sensitivity analysis method to reduce the number of installation nodes to be selected, and a feasible alternative set of installation nodes for the optimal conﬁguration model of ETSHE could be obtained. Finally, the improved multi-objective particle swarm optimization algorithm was used to solve the model, and the data envelope analysis method was used to evaluate the power quality of each access scheme. The analysis of the numerical example shows that it can not only satisfy the user (cid:48) s heat demand, but also e ﬀ ectively improve the power quality by rationally planning the location and capacity of ETSHE, which achieves the safe and e ﬃ cient utilization of energy.


Introduction
In recent years, large-scale use of fossil fuels has led to exacerbating the environmental pollution problems, thus energy reform is imperative. New energy power generation is widely developed worldwide, among them, distributed photovoltaic and wind power generation are widely used because they have the characteristics of abundant available resources and reliable devices [1]. However, from another perspective, distributed photovoltaic and wind power also have the characteristics of fluctuation and stochastic. When the output of a distributed photovoltaic or wind power connected to a node in the power distribution system exceeds the node s self-consumption, the node voltage will exceed the limit. At the same time, a large number of distributed photovoltaics use power electronic equipment, which will bring a lot of harmonics to the distribution network, and can cause serious 1.
First, the minimum intra-class distance is used as the objective function to perform fuzzy clustering on the heat load. Second, on the basis of conventional forecasts, the influence weight of the meteorological factors on the outdoor temperature of the day is calculated to correct the outdoor temperature, improve the linear fitting goodness with heating heat load, and provide initial classified heat load data for the optimal configuration of ETSHE.

2.
The sensitivity coefficient is calculated based on the sensitivity radar chart and sensitivity calculation formula. The sensitivity discrimination of each node in the distribution network can effectively eliminate the nodes that do not meet the sensitivity and prevents the solution from falling into the wrong node, which can improve the efficiency and accuracy of the solution. 3.
Establish a location and capacity planning model for ETSHE and multi-objective functions including harmonic distortion, voltage deviation, voltage fluctuation, and energy storage equipment costs. 4.
In the planning stage, data envelopment analysis (DEA) is used to evaluate the power quality of access schemes with multiple indicators and single indicators, and an effective and unified method for distribution network planning and power quality assessment is proposed.

Load Characteristic Portrait
Traditional load forecasting methods generally only consider the typical daily curve of heat load, but not consider the differences in the heating characteristics of classified heat loads. Under the background of the rapid development of big data technology, portrait technology has been widely used in the fields of e-commerce and the internet [14]. This paper applied portrait technology to the extraction of heating user characteristics for the first time. In this paper, the thermal user load portrait was defined as follows. Different thermal loads were predicted by a high-precision prediction method so that the thermal loads were treated differently according to the different heating characteristics. To comprehensively consider the user heating characteristics and influencing factors, we established a user heating behavior label library. Based on the fuzzy C-means (FCM) algorithm for label clustering, different types of heating users used warming mode portraits, which makes the results of user heating behavior analysis more intuitive [15].

Clustering with Warm Behavior Patterns
The goal of users label clustering is to quickly grasp the heating characteristics of different thermal users to achieve different treatments of different heating groups. After combing and screening many times, the following six labels were obtained, as shown in Table 1. In this paper, FCM clustering algorithm was used to calculate each thermal user's label data to be clustered, and four memberships of thermal user clustering centers to realize the clustering division of the sample thermal user's label data. During the iterative process, FCM takes the minimum intra-cluster distance and the maximum distance between clusters as the optimization goals, and constantly adjusts the categories of various thermal user data, so that it optimizes the clustering results continuously. The object to be clustered is the D-dimensional dataset of N samples: . . . . . . . . .
The data collection class is divided into K subsets, where the value of K is 4, the subset is represented as C = {c 1 , · · · , c i , · · · , c 4 }; and c 1 , c 2 , c 3 , and c 4 represent residents, education, business and industry, respectively. Here, we set D = 7, which represents six items of user label data and one category data, respectively. Assuming that the probability of each thermal user j belonging to the ith clustering subset is u ij , then FCM can be defined as an optimization problem with equality constraints [15]: where J is the value function of FCM. m is the membership factor, which is used to adjust the degree of influence of the membership indicator in the optimization target. In most cases, we chose m = 2. The iterative formulas for obtaining membership coefficients and cluster centers are shown in Equations (2) and ( Clustering steps of the warm behavior patterns for thermal users based on the FCM algorithm:

1.
Initialization: Set the number of clustering K and membership factor m, randomly initialize the membership matrix U, U is a matrix made up of u ij , the center of the cluster.

2.
Parameter update: According to Equation (2), calculate the clustering center C of each label of thermal users; according to Equation (3), calculate the membership matrix U of all thermal user labels and each clustering center C, classify thermal users into the clustering center with the highest membership.

3.
Convergence judgment: If the change of the objective function is less than the preset threshold, then output the clustering result, otherwise return to Step 2.

4.
End of clustering: Save and display the clustering results.

Load Portrait
Load portrait refers to the method of high-precision prediction of the heating heat load by considering the main factors affecting heating load such as outdoor temperature, solar radiation, natural wind, etc. Correct outdoor temperature can be performed by using solar radiation and wind Appl. Sci. 2020, 10, 2666 5 of 27 speed to improve the fit between outdoor temperature and thermal load [16]. Wind speed and solar radiation are the main factors affecting the heating heat load. The correction of wind speed for outdoor temperature can be converted into an equivalent cooling temperature [17]. The equivalent cooling temperature formula considering the influence of wind speed is shown in Equation (4): where W wind is the external wind speed value, m/s; ∆T wind is the equivalent temperature of wind speed cooling; • C. α 1 , α 2 , α 3 , and β are the wind speed correction coefficients, where the values are 0.0246, 0.4525, 3.2398, and 7.23, respectively. The correction of outdoor temperature by solar radiation can be converted into an equivalent heating temperature as shown in Equation (5): where S solar is the external light value; W/m 2 , S k is the light conversion coefficient; W/ m 2• C ; according to [17], we set the unit of S k as 100 W/ m 2 • C . T out is the uncorrected outdoor temperature. When considering both the outdoor wind speed and solar radiation of the day, the calculation formula for outdoor temperature T w,s is given as follows: In order to improve the accuracy of forecasting heating load, this study not only considered the influence of wind speed and solar radiation on outdoor temperature of the day, but also considered the influence of wind speed and solar radiation on outdoor temperature of the past day. In this paper, we calculated the influence weight of previous meteorological factors on the outdoor temperature of the day to modify the outdoor temperature and linearly fit it with the heating heat load, so that their fitting goodness was further improved [18][19][20].
According to heat transfer knowledge, the influence of meteorological factors on the outdoor temperature of the day decreases with the increase in days. In order to facilitate the prediction of heating heat load, the influence of meteorological factors of all past days on the outdoor temperature of that day cannot be considered. Statistical results showed that heating heat load was related to meteorological factors of the previous 4-5 days. We used MATLAB programming to find the optimal days and the reasonable weight of past weather factors to correct the outdoor temperature of the day. When the outdoor temperature and the heating heat load were linearly fitted, the goodness of fit value was the maximum. This means that when the number of days considered is n days, the outdoor temperature of the day is given as follows: where , represent the weight value of day i, day i − 1, day i − 2, and day i − n, respectively, and , and the sum of weight values is 1. T w,(i−n) is the outdoor temperature after solar radiation and wind speed correction of day i − n, and T w,i,e is the outdoor temperature of day i after correction on the current day and the previous day. The "trial" method was adopted to determine the index weight of the influence of past meteorological factors on the outdoor temperature of the day to correct the outdoor temperature and obtain the relationship between the heat load and the temperature of the day, then obtained the relationship between the heat load and time. For example, if n = 2 in Equation (7), we can set A i = m 0 /(m 0 + m 1 + m 2 ), A (i−1) = m 1 /(m 0 + m 1 + m 2 ), and A (i−2) = m 2 /(m 0 + m 1 + m 2 ) to simplify the calculation, and m 0 , m 1 , m 2 is an integer between 1 and 10.

Optimal Configuration Model of Electric Thermal Storage Heating Equipment
High-permeability distributed photovoltaic and wind power connected to the distribution network cause a series of power quality problems, which can be improved by installing ETSHE at appropriate nodes of the distribution network. However, the power quality of the power grid involves many aspects such as voltage and current. The location and capacity of the energy storage system needs to consider the effect of improving the power quality, the consumption of new energy, and the economics of the energy storage system. It is a multi-objective optimization plan model.

Total Current Harmonic Distortion
The fundamental current flowing through the power load is equal to the sum of the fundamental current of the power grid and the fundamental current of the inverter. When the power load is unchanged, the current harmonics are unchanged, but the fundamental current of the distribution network is reduced, and the total current harmonic distortion rate at the access point is increased. With the increase of access capacity, the total current harmonic distortion rate of the distribution network may exceed the range of grid access energy. Therefore, the total distortion rate of current harmonics was selected as one of the objective functions of the energy storage system location and capacity. The mathematical expression is [21]: where I 1 is the root mean square value of the fundamental current and I h is the root mean square value of the h-th harmonic current. Here, we set M = 7, where only the 7th and lower harmonics were considered.

Voltage Deviation
Voltage is one of the most important indicators of power quality, and voltage deviation is a main indicator to measure the normal operation of the power supply system [22]. Distributed photovoltaic access to the distribution network can increase the voltage of the nodes, and excessive photovoltaic output will cause some nodes to exceed the limit. In this paper, the node voltage deviation was selected as one of the objective functions of the energy storage system location and capacity determination, its mathematical expression is [13]: where U ij is the voltage of node i at time j; U N is the system standard voltage; and N is the number of system nodes. T is the number of detection moments.

Voltage Fluctuation
As the output power of the device is fluctuating, the voltage fluctuation caused by the connection to the distribution network is ∆d, so the voltage fluctuation was selected as one of the objective functions of this article. The mathematical expression is [23]: where U N is the rated voltage; ∆I is the current change caused by the output power change at the access point; and Z is the equivalent impedance of the two-port network distribution network viewed at the access point. The ETSHE must also consider the cost issue while ensuring the quality of the power. The total cost of the energy storage system was selected as one of the objective functions.

Equipment Cost of Electric Boiler
The equipment cost of an electric boiler is proportional to its rated heating capacity, which is equal to the design daily heat supply rating, thus Equation (11) is: where ρ S is the unit price corresponding to the unit heat supply of the electric boiler. Q e S (Q S ) is the rated heat supply of the electric boiler. Q S = [Q S (1) Q S (2) . . . Q S (24)] T is a vector composed of the hourly heating capacity of the designed electric power boiler.

Equipment Cost of Heat Storage Body
The cost of the heat storage equipment is proportional to the volume of the equipment, then C X (V X ) = ρ X V X , and ρ X , V X are the unit price and volume of the heat storage body, respectively.

Equipment Cost of Heat Storage Pump
The key parameters of heat pump selection are head and flow. The head is only related to the building structure and was selected according to typical values, and the heat pump investment was proportional to the rated flow at a certain head [24,25]. Therefore, the heat pump equipment cost can be expressed by Equation (12): where ρ B is the increased equipment investment for each additional unit of the rated flow of the heat pump. γ is the safety factor when calculating the rated flow, in this paper, γ = 1.05. C 1 is the conversion constant of kilowatt hours and joules, here C 1 = 3.6 × 10 6 . θ = [θ(1) θ(2) . . . θ(24)] T is a vector formed by designing the hourly direct supply ratio of the electric power boiler, and C W ρ W are the specific heat and density of the heat medium. The difference in heat storage temperature is usually ∆τ T = 40 • C [26].

Equipment Cost of Control Device
The equipment cost of the control device is related to the complexity of the electric boiler system. It was found that the equipment cost of the control device is basically proportional to the rated heating capacity of the heating device of the electric boiler, which is shown as follows: where ρ C is the investment cost of the control device for each additional 1 kW of the rated heating capacity of the electric boiler. This article does not consider the operation and maintenance costs of energy storage equipment, only the investment and installation costs of energy storage equipment were considered, so the total cost of ETSHE is given as Equation (14): where k ins is the equipment installation cost coefficient. Usually, the equipment installation cost is about 20% of the total equipment cost, so we set k ins = 1.2.
Therefore, the location and capacity multi-objective optimization model of the energy storage system comprehensively considered the voltage deviation of the distribution network nodes, the total distortion rate of the current harmonics, the voltage fluctuation, and the total cost of the energy storage system. The objective function is shown as Equation (15): where , ϑ, τ, ς is the normalized weight coefficient.

Power flow constraint equation
The initial power flow calculation of the planned distribution network is performed, the voltage value and active power of each node are calculated. The power flow constraint is as follows [27]: where P i and Q i are the active and reactive power injections at node i. U i and U j are the voltage amplitudes of nodes i and j. G ij and B ij are the conductance and susceptance of branch i j. θ ij is the voltage phase angle difference between the i-th and j-th node. 2. The hourly balance of heat supply and heat demand is: where θ = [θ(1) θ(2) . . . θ(24)] T is a vector formed by the design of the hourly direct supply ratio of the electric power boiler. Q S = [Q S (1) Q S (2) . . . Q S (24)] T is the vector for designing the hourly heating capacity of the electric power boiler. Q q (t) is the heat provided by the heat storage body to the user during period t. Q l,max is the maximum design daily heat load. r l (t) is the ratio of the heat load to the maximum value for the t days on the design day. 3. The thermal equilibrium equation of the heat storage body is: where Q N (t) is the amount of heat stored in the heat storage body at the beginning of time t.

Inequality Constraint
1. The maximum heat storage capacity is limited by the volume of the heat storage body: where V X is the volume of the heat storage body, and C W and ρ W are the specific heat and density of the heat medium, respectively. C 2 is the conversion constant of kilowatt-hours and joules, here C 2 = 3.6 × 10 6 . Usually, the difference in heat storage temperature is ∆τ T = 40 • C. υ X is the effective volume ratio of the heat storage body, which refers to the ratio of the volume that can actually contain the heat medium to the total external volume. 2. The upper limit of the heat supply of the heat storage body is: Appl. Sci. 2020, 10, 2666 9 of 27 3. The range of the value of the direct heat supply ratio of electric boilers is constrained as: 4. Constraints of voltage levels before and after the ETSHE is connected to the distribution network: where L DG and L dis are the voltage level indicators before and after the energy storage is connected to the distribution network, respectively; ε 1 is between 0 and 1. 5. Constraints on active loss before and after the ETSHE is connected to the distribution network: where P loss,DG and P loss,dis are the active power loss before and after the ETSHE is connected to the distribution network, respectively. Similarly, ε 2 is between 0 and 1. The objective functions and constraints of the optimization model in this paper were grouped as shown in Figure 1: to the distribution network, respectively; is between 0 and 1.

304
5. Constraints on active loss before and after the ETSHE is connected to the distribution network: where , and , are the active power loss before and after the ETSHE is connected to the 306 distribution network, respectively. Similarly, is between 0 and 1.

307
The objective functions and constraints of the optimization model in this paper were grouped as 308 shown in Figure

312
For power networks with a small number of nodes, the optimal configuration scheme for the 313 location and capacity of ETSHE can be easily calculated by exhaustive methods. When the power 314 network is complicated, there are many installation nodes available for installing energy storage at 315 this time. If the optimal node is still selected according to the exhaustive method, the calculation of 316 the optimal allocation scheme is large, which cannot be solved or the result falls into a local optimum.

317
For this reason, the sensitivity analysis method was proposed to calculate the sensitivity coefficient 318 of electric regenerative heating equipment installed at the node, and the node with higher sensitivity 319 was selected as the alternative set of the installation node of ETHSE, which will greatly improve the 320 optimal allocation efficiency of electrical thermal storage.

321
Sensitivity coefficient refers to the degree of sensitivity of the system, represented by state 322 variables to changes in the control variables and disturbance variables [28]. For the electric storage

Node Sensitivity Analysis and Calculation Method
For power networks with a small number of nodes, the optimal configuration scheme for the location and capacity of ETSHE can be easily calculated by exhaustive methods. When the power network is complicated, there are many installation nodes available for installing energy storage at this time. If the optimal node is still selected according to the exhaustive method, the calculation of the optimal allocation scheme is large, which cannot be solved or the result falls into a local optimum. For this reason, the sensitivity analysis method was proposed to calculate the sensitivity coefficient of electric regenerative heating equipment installed at the node, and the node with higher sensitivity was selected as the alternative set of the installation node of ETHSE, which will greatly improve the optimal allocation efficiency of electrical thermal storage.
Sensitivity coefficient refers to the degree of sensitivity of the system, represented by state variables to changes in the control variables and disturbance variables [28]. For the electric storage thermal heating equipment location and capacity optimization configuration model, it refers to the sensitivity of the objective function of the optimal configuration to changes in energy storage installation capacity. The optimization objective of this paper is multi-objective function and the fundamental purpose of installing energy storage is to effectively control the power quality problem of the distribution system. When a large number of single-phase photovoltaic power sources are connected to the system, it will cause voltage three-phase imbalance problems. However, as single-phase photovoltaic power generation has not been widely used, it has little influence on the three-phase unbalance degree of power grid voltage. Therefore, only the three-phase imbalance is considered when defining sensitivity.
In a three-phase three-wire system, there is no zero-sequence component. For a three-phase system without a zero-sequence component, the simplified calculation method of the three-phase imbalance recommended by the national standard is as follows [29]: where f ε2 represents the three-phase negative sequence imbalance, and a, b, and c represent the effective value of three-phase voltage, respectively. In this paper, the radar chart analysis method was used to visually and comprehensively analyze and evaluate the sensitivity of the nodes. Finally, the nodes with higher sensitivity were selected as candidate nodes. The schematic diagram of the sensitivity radar is shown in Figure 2 where represents the three-phase negative sequence imbalance, and , , and represent the 336 effective value of three-phase voltage, respectively.

337
In this paper, the radar chart analysis method was used to visually and comprehensively analyze 338 and evaluate the sensitivity of the nodes. Finally, the nodes with higher sensitivity were selected as  349 Therefore, in this paper, when the four power quality indicators of any node were defined as 350 "−", it was considered that the sensitivity of this node to power quality reached the standard and 351 could be used as an alternative node. The specific calculation method is: When the sensitivity coefficient = 1, the node can be used as a candidate node. If = 0, 353 the node is eliminated.

354
The sensitivity coefficient reflects the degree of sensitivity of any node to the four aspects 355 of harmonic content, voltage deviation, voltage fluctuation, and three-phase imbalance. If the energy storage system is installed at this node, the higher the value, the better the power quality will be.

357
Therefore, by calculating the sensitivity coefficients of all nodes in the system, it is possible to find in 358 advance a set of alternative nodes that are better at improving the power quality of the system, which Therefore, in this paper, when the four power quality indicators of any node were defined as "−", it was considered that the sensitivity of this node to power quality reached the standard and could be used as an alternative node. The specific calculation method is: When the sensitivity coefficient α sen = 1, the node can be used as a candidate node. If α sen = 0, the node is eliminated.
The sensitivity coefficient α sen reflects the degree of sensitivity of any node to the four aspects of harmonic content, voltage deviation, voltage fluctuation, and three-phase imbalance. If the energy storage system is installed at this node, the higher the value, the better the power quality will be. Therefore, by calculating the sensitivity coefficients of all nodes in the system, it is possible to find in advance a set of alternative nodes that are better at improving the power quality of the system, which can effectively reduce the number of computing nodes, greatly reduce the amount of calculations, and improve the optimization of the economy and efficiency of configuration.
Nodes are ranked according to the sensitivity coefficient, and higher-sensitivity nodes are preferentially selected to form an initial optimized candidate set. Subsequently, the heuristic optimization algorithm can be applied to solve the model quickly. This paper used the adaptive inertial weight IMOPSO algorithm to solve the model.
The main steps of the adaptive inertia weight IMOPSO algorithm are as follows: 1.
The network topology, switch state, node voltage range, distributed photovoltaic, wind parameters, and parameters of ETSHE in the distribution network system are initialized. The sensitivity coefficients of each node are calculated and sorted separately. The top 40% of the sorted nodes are selected to form a candidate set of installation nodes.

2.
The parameters such as the size of the particle swarm, the maximum number of iterations of the algorithm, the acceleration factor, the value range of the adaptive inertia weight, and the speed range of the particles are set in the improved particle swarm optimization. 3.
The particles in the candidate set are evaluated, and the fitness function value of each particle is calculated. After ranking, the individual optimal value and the global optimal value of the particles are determined and updated. 4.
Update the particle s velocity and displacement, calculate the value of the adaptive inertia parameter.

5.
Check if the maximum number of iterations is exceeded. If it is exceeded, stop iteration and output the calculation result, otherwise skip to Step 3 to continue the iterative calculation.
In this algorithm program, the total current harmonic distortion rate, the sum of voltage deviations, the voltage fluctuations, the investment, and installation cost of the energy storage system are optimized objective functions. The number of iterations is the cycle control variable, and the location and capacity of the ETSHE are taken as decision variables. By improving the solution of the particle swarm algorithm, the optimal position and optimal capacity for the energy storage installation can be obtained.
The solution flow chart of IMOPSO algorithm in this paper is shown in Figure 3.
Appl. Sci. 2019, 9, x FOR PEER REVIEW 11 of 29 optimization algorithm can be applied to solve the model quickly. This paper used the adaptive 364 inertial weight IMOPSO algorithm to solve the model.

365
The main steps of the adaptive inertia weight IMOPSO algorithm are as follows:

386
The solution flow chart of IMOPSO algorithm in this paper is shown in Figure 3.

Simulation Analysis
This paper used IEEE-33 and IEEE-118 power distribution systems for verification, respectively, which were used to analyze the impact of high-permeability wind power and photovoltaic access on the power quality of the power distribution network. Through the rational optimization of flexible resources such as ETSHE, the power quality of high-permeability wind power and photovoltaic power distribution networks is improved.

Load Portrait Result
Using FCM clustering algorithm for label clustering, where the value of N was 200 in this paper, we obtained four typical thermal users, as shown in Table 2. Figure 4 shows the distribution diagram of the FCM clustering results. According to the membership matrix in Figure 4, there were some differences in the characteristics of the four types of heating customers, but there were also some overlaps, which will affect the planning of the ETSHE. Equations (27) and (28) were used to calculate the fitting relationship between T w,i,e and Q, and the four types of typical users' heat consumption over time were further calculated, as shown in Figure 5. The results of the thermal load portrait were input into the PSO as the initial data.     It can be seen from the results that when considering the influence of meteorological factors on the outdoor temperature of the day in the past two days, the goodness-of-fit R 2 of the linear fitting with the heating heat load was a maximum value after the correction of the meteorological factors on that day. It shows that the curve at this time can best reflect the relationship between Q and T, and the value of Q can be predicted from the modified T. For this set, the heating heat load can be predicted by the following formula: The predicted values of different heating heat loads were obtained from Equations (27) and (28). Within a day, the temperature changes with time, and the heating heat load changes with air temperature, so the curve of the heating load with time can be obtained equivalently. The energy consumption characteristics of four typical users are shown in Figure 5, where the typical users energy consumption characteristic curves show different trends.

Simulation Scenario Setting
We used two scenarios to analyze the effects of the optimal scheduling of electrical thermal storage heating equipment and load demand characteristics on the improvement in the power quality of the distribution network.
Scenario 1: In the absence of flexible resources, the indicators of power quality include high-permeability wind power and photovoltaic power distribution networks. Scenario 2: In the case of ETSHE and load demand characteristics affecting the optimal dispatching scenario, various indicators of power quality include high-permeability wind power and photovoltaic distribution networks.

Basic data
The distribution network load, photovoltaic, and wind power data are shown in Table 3. The installation locations of wind power and photovoltaic are shown in Figure 6. The distribution network contains twelve distributed photovoltaics with a capacity of 0.5 MW, which are connected to nodes 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, and 24. At the same time, the distribution network also contained four distributed wind power generation with a capacity of 0.5 MW, which were connected to nodes 26, 28, 30, and 32. In order to ensure the maximum consumption of renewable energy, each distributed power supply is operated with actual output. planning settings, we considered that nodes 4 and 13 were residential loads, node 8 was a commercial 451 load, node 25 was an industrial load, and node 31 was an educational load.

452
According to the optimization model constructed in this paper, the IMOPSO algorithm was used  Table 4.

Result analysis (1) Optimal scheduling results of ETSHE considering load characteristics
In principle, the nodes that can be connected to ETSHE are 1-32, but harmonics, voltage deviations, and voltage fluctuations caused by high-permeability distributed wind power and photovoltaic will cause some nodes to cancel each other out. Considering the objective function, the node sensitivity analysis method was applied, and the set of nodes for energy storage access included five nodes: nodes 4, 8, 13, 25, and 31. In order to highlight the impact of load characteristics on optimal planning settings, we considered that nodes 4 and 13 were residential loads, node 8 was a commercial load, node 25 was an industrial load, and node 31 was an educational load.
According to the optimization model constructed in this paper, the IMOPSO algorithm was used to solve the calculation example. The population size of the IMOPSO algorithm was 80, and µ min = 0.5, µ max = 0.8, σ = 0.2 [30,31]. The installation position and capacity of the electric regenerative heating equipment are shown in Table 4. (2) Node sensitivity analysis before and after comparison As can be seen from Table 5, the range of optimized nodes is effectively reduced and the effectiveness of the solution is greatly improved by the use of the node sensitivity analysis method. Table 5. Comparison analysis before and after node sensitivity analysis.

Optimal Planning Model of ETSHE Before Analysis After Analysis Effectiveness Enhancement
Calculation time/s 153.6 121.5 60% (3) Power quality improvement analysis Through the optimization of two flexible resources, the power quality of distributed wind power with high penetration and photovoltaic power distribution networks has been significantly improved. Figures 7-9 are the comparison curves of the total current harmonic distortion rate, voltage deviation, voltage fluctuation of the distribution network with high penetration wind power, and photovoltaic access in two scenarios. Line charts show the comparison of different scenarios of node power quality indicators, and the bar charts show the percentage improvement of the two scenarios relative to the national standard value.
According to the GB/T14549 standard, the total current distortion rate of the electric heat storage heating equipment after being connected to the distribution network should be less than or equal to 2.2% [21,32]. Figure 7a,b are the comparison curves of the total distortion rate of the node current harmonics in the two scenarios, and Figure 7c is the percentage increase of the harmonic distortion rate of the node current relative to the national standard value. As can be seen from Figure 7a,b, in Scenario 1, nodes 4 and 31 have a high light intensity between 8:00 am and 7:00 pm. Due to the introduction of a large number of power electronic equipment by photovoltaics, the current harmonic distortion rate exceeded the standard. To ensure the stability of the system, wind abandonment, light abandonment, and load shedding must occur. In Scenario 2, the current harmonic distortion rate was successfully reduced below the standard value, and the improvement rates were 45% and 47%, respectively. Although nodes 8 and 25 did not exceed the standard, the percentage of improvement through energy storage optimization can reach up to 50% and 26%. As can be seen from Figure 7c, the percentage of improvement of power quality relative to the standard value at each time was a positive value, which means that the power quality had improved. According to the GB/T14549 standard, the total current distortion rate of the electric heat storage heating equipment after being connected to the distribution network should be less than or equal to 2.2% [21,32]. Figures 7a,b are the comparison curves of the total distortion rate of the node current harmonics in the two scenarios, and Figure 7c is the percentage increase of the harmonic distortion rate of the node current relative to the national standard value. As can be seen from Figures 7a,b, in Scenario 1, nodes 4 and 31 have a high light intensity between 8:00 am and 7:00 pm. Due to the introduction of a large number of power electronic equipment by photovoltaics, the current harmonic According to GB/T12325, the absolute value of the voltage deviation should be less than or equal to 7% [33]. As can be seen from Figure 8, taking into account the load characteristics and the reasonable planning of the ETSHE, the voltage deviation of each node was significantly reduced and the voltage deviation of the access point at each moment was within the allowed range of the national standard. As can be seen from Figure 8a,b, at 12:00 am, the light intensity was the highest, but the corresponding voltage deviation was the lowest. The voltage deviation was inversely proportional to both wind speed and light intensity. Furthermore, as can be seen from Figure 8c, the percentage of each node s improvement relative to the standard value was positive, with an average improvement rate of 38%, and the voltage deviation of each node was improved.
Appl. Sci. 2019, 9, x FOR PEER REVIEW 17 of 29 distortion rate exceeded the standard. To ensure the stability of the system, wind abandonment, light abandonment, and load shedding must occur. In Scenario 2, the current harmonic distortion rate was successfully reduced below the standard value, and the improvement rates were 45% and 47%, respectively. Although nodes 8 and 25 did not exceed the standard, the percentage of improvement through energy storage optimization can reach up to 50% and 26%. As can be seen from Figure 7c, the percentage of improvement of power quality relative to the standard value at each time was a positive value, which means that the power quality had improved. From 8 am to 7 pm, ETSHE mainly used the abandoned light resources, and other times mainly used abandoned wind resources. According to the GB/T12326-2008 standard, the voltage fluctuation of a 10 kV distribution network should be less than or equal to 4% [23]. As can be seen from Figure 9a,b, through the optimized configuration of energy storage, the node voltage fluctuation was improved. At 12:00 am every day, the light intensity was the highest, but the voltage fluctuation is the lowest. It can be seen from Figure 9c that at the time when there was more wind power and photovoltaic power generation, the voltage fluctuation was greatly improved relative to the standard value. Among them, node 4 increased by 23% on average, nodes 8 and 25 increased by about 25%, and node 31 increased by 11.2% on average. According to GB/T12325, the absolute value of the voltage deviation should be less than or equal to 7% [33]. As can be seen from Figure 8, taking into account the load characteristics and the reasonable planning of the ETSHE, the voltage deviation of each node was significantly reduced and the voltage deviation of the access point at each moment was within the allowed range of the national standard. As can be seen from Figures 8a,b, at 12:00 am, the light intensity was the highest, but the corresponding voltage deviation was the lowest. The voltage deviation was inversely proportional to both wind speed and light intensity. Furthermore, as can be seen from Figure 8c, the percentage of each node′s improvement relative to the standard value was positive, with an average improvement rate of 38%, and the voltage deviation of each node was improved. According to the GB/T12326-2008 standard, the voltage fluctuation of a 10 kV distribution network should be less than or equal to 4% [23]. As can be seen from Figures 9a,b, through the optimized configuration of energy storage, the node voltage fluctuation was improved. At 12:00 am every day, the light intensity was the highest, but the voltage fluctuation is the lowest. It can be seen from Figure 9c that at the time when there was more wind power and photovoltaic power generation,
According to the optimization model constructed in this paper, the IMOPSO algorithm was used to solve the calculation example. The installation position and capacity of the electric regenerative heating equipment are shown in Table 6. (2) Node sensitivity analysis before and after comparison This system is more complex than IEEE-33, so it takes a longer time to calculate. The comparative analysis before and after the node sensitivity analysis is shown in Table 7. (3) Power quality improvement analysis From the above simulation results, it can be seen that the power quality of residential and educational loads was similarly affected by new energy generation, while commercial and industrial loads were similar. Therefore, the residential and industrial loads were selected to give the results for node 3 and node 44, respectively.
It can be seen from the simulation results that although the IEEE-118 system had a long and complex calculation time, the node sensitivity calculation could obtain a set of candidate nodes in advance, which improved the efficiency and accuracy of the solution. The same results as those of the IEEE-33 nodes' power distribution system can be seen in Figures 10-15. Through the simulation analysis of the above IEEE-33 and IEEE-118 systems, by configuring a certain capacity of ETHSE at specific nodes of the distribution network, it can not only solve the voltage deviations and voltage fluctuations caused by the photovoltaic and wind power, but can also restrain harmonics, which can improve the power quality of distribution network effectively. Finally, it can be seen that the evaluation levels were all above five levels by using the data envelope analysis method, which indicates that the access scheme is feasible and that the indicators of the access point can reach the prescribed standards. For specific evaluation methods and evaluation levels, see Appendix A.

Conclusions
Aiming at the power quality problems caused by high-permeability distributed photovoltaic and wind power connected to the grid, this paper proposed a method for selecting and determining the capacity of the ETSHE connected to the distribution network, which considered the load characteristics and power quality management. A load forecasting method based on load profile technology was proposed, and the node sensitivity analysis method was used to accurately and quickly solve the multi-objective optimization model, which realized the management of power quality problems such as node voltage overrun, voltage fluctuation, and harmonic pollution. Moreover, the method satisfied the high demand for penetration of photovoltaic and wind power into the distribution network. The main innovations of the work are: (a) comprehensively considering the impact of demand-side response and power quality issues on distribution network planning, and establishing an optimization model that includes power quality indicators and economic indicators; (b) based on the node sensitivity analysis method, so that it can be applied to complex power distribution systems; and (c) the heat load portrait technology is defined to provide a data basis for subsequent optimization. It can provide references for the promotion and application of distributed photovoltaic and wind power, and provide solutions for power quality management methods caused by high-power density distributed energy grid-connection. The next work will consider multiple types of energy storage equipment to further optimize the economy and safety.

Conflicts of Interest:
The authors declare no conflicts of interest. The ratio of heat load to maximum value in period t of design day Q N (t)

Abbreviations
Heat stored in the heat storage body at the beginning of time t [kWh] α sen Sensitivity coefficient w j Relative production efficiency u r Weight of the r-th power quality indicator v i Weight of the i-th quality level x rj r-th power quality input indicator in the j-th DMU y ij i-th quality level output in the j-th DMU indicator

Appendix A
• Comprehensive evaluation method of power quality based on data envelope analysis DEA is a comprehensive evaluation method for evaluating the relative efficiency of the same type of multi-index inputs and multi-index outputs [35]. In this method, the weights of the input and output variables are determined from the actual data of each decision unit, and the evaluation is performed from the perspective of the most favorable for the decision unit, which can avoid the influence of subjective factors. Additionally, it is not necessary to determine the display expressions of input and output, which simplifies the internalized relationship between input indicators and output indicators. In the case where the relationship between input indicators and output indicators is complex and difficult to determine, it can reflect its obvious advantages. The data envelopment analysis method was used to evaluate the power quality of each access point, and the access plan was evaluated for power quality at the planning stage, and finally, an ETSHE access plan with a good power quality level was obtained. The basic model and simplified transformation of DEA are as follows: (1) C 2 R Model The C 2 R model is the most basic DEA model. Suppose there are n evaluation objects, each Decision Making Units (DMU) has m types of input indicators and q types of output indicators. Construct the following C 2 R model [35]: where w j is the relative production efficiency; u r is the weight of the r-th power quality indicator (r = 1,2, . . . , m); v i indicates the weight of the i-th quality level (i = 1,2, . . . , q); x r j indicates the r-th power quality input indicator in the j-th DMU; and y i j indicates the i-th quality level output in the j-th DMU indicator (j = 1,2, . . . , n). Choose a reasonable input weight vector v i and output weight vector u r so that w j ≤ 1. When w j = 1, the j-th DMU is considered to be the most effective.
(2) Charnes-Cooper transformation Due to the difficulty of solving the C 2 R model, the fractional programming is transformed into an equivalent linear programming problem using the Charnes-Cooper transformation, and the dual model is obtained using dual theory.
• Power quality evaluation of the ETSHE planning scheme The power quality evaluation is carried out for the planning scheme of the ETSHE to be connected to the distribution network. First, the classification of grades should be carried out. There are many ways of classification, among which domestic literatures are mainly divided into grades 3, 5, 7, 9, and 10. This paper divided the power quality indicators into five grades. The assessment level was Q, which decreased from Q1-Q5 in order: high quality, good, average, qualified, and unqualified.
Considering the influence of the ETSHE connected to the distribution network, the evaluation indicators selected here were the total current harmonic distortion rate, voltage deviation, and voltage fluctuation. X1, X2, and X3 represent the evaluation indicators, respectively. According to the GB/T14549 standard, the total current distortion rate of ETSHE, after it is connected to the distribution network, should be less than or equal to 2.2% [21,32]. According to the standard of GB/T12326-2008, the voltage fluctuation of the 10 kV distribution network should be less than or equal to 4% [23]. According to GB/T12325, the absolute value of the voltage deviation is less than or equal to 7% [33]. According to the national standard, the upper limit is equally divided into grade intervals, and the specific grade standards are shown in Table A1. Table A1. Classification of power quality indicators.

X1 X2 X3
In this paper, the data envelopment analysis method was used to evaluate the power quality, with X1, X2, and X3 as the inputs and the evaluation level Q as the output to establish a DMU. The DMU format was (X1, X2, X3, Q). The median and evaluation level of each level interval in Table 6 were taken as a single decision unit. The data of the standard decision unit of the classification are shown in Table A2. The power quality index grade standard in Table 6 and the DMU of the hierarchical standard in Table A2 were used as the total decision unit to perform the data envelope validity analysis. Calculate the quality indicators of the access node according to the formula and the average of each indicator. According to the calculation formulas of the total current harmonic distortion rate, voltage deviation, voltage fluctuation, the quality indicators of the power quality of the access node on that day were calculated, respectively, and the average values of the indicators were calculated, as shown in Table A3. As the efficiency value of the power quality indicators of each access point relative to each level was different from the decision unit, it was recorded as δ Q , and the level with the smallest absolute value of δ Q was used as the evaluation level. Table A4 shows the power quality evaluation results of each access point calculated by the data envelope analysis method. As can be seen from the result of the power quality evaluation of each access point by the data envelope analysis method, only the power quality evaluation level of access points 4 and 8 was a level 4 deviation. The level of access point 13 was level 2, and the power quality evaluation level of the access points at 25 and 31 was a quality of level 1. The evaluation level without access points was a level 5 failure, indicating that the access scheme is feasible and that all indicators of the access point can meet the requirements of the standard.