Reactive power optimization of a distribution network with high-penetration of wind and solar renewable energy and electric vehicles

As high amounts of new energy and electric vehicle (EV) charging stations are connected to the distribution network, the voltage deviations are likely to occur, which will further affect the power quality. It is challenging to manage high quality voltage control of a distribution network only relying on the traditional reactive power control mode. If the reactive power regulation potentials of new energy and EVs can be tapped, it will greatly reduce the reactive power optimization pressure on the network. Keeping this in mind, our reasearch first adds EVs to the traditional distribution network model with new forms of energy, and then a multi-objective optimization model, with achieving the lowest line loss, voltage deviation, and the highest static voltage stability margin as its objectives, is constructed. Meanwihile, the corresponding model parameters are set under different climate and equipment conditions. Ultimately, the optimization model under specific scenarios is obtained. Furthermore, considering the supply and demand relationship of the network, an improved technique for order preference by similarity to an ideal solution decision method is proposed, which aims to judge the adaptability of different algorithms to the optimized model, so as to select a most suitable algorithm for the problem. Finally, a comparison is made between the constructed model and a model without new energy. The results reveal that the constructed model can provide a high quality reactive power regulation strategy.


Introduction
Renewable energy has characteristics of sustainability, cleanliness and, often, inexhaustible supply. Research has shown that renewable/new energy systems can not only meet active load demand of the power grid, but also achieve rapid reactive power regulation using power electronic devices connected to the network [1][2][3].
However, with large-scale renewable energy connecting to the distribution network, the existence of many power electronic devices in the power network can cause severe reactive power impact, which can seriously affect power quality and power system voltage stability [4]. In addition, because of the influence and unpredictability of the natural environment, there will be a great uncertainty in the reactive power regulation of the distribution network for new energy.
To achieve carbon neutrality, the demand for electric vehicles (EVs) increases in association with the need for emission reduction. EVs may participate in regulating the output of the power grid and trading energy with the grid based on vehicle to grid (V2G), to help realize benign interaction between EVs and the power grid. The bi-directional charging technology of EVs has great flexibility. A certain number of EVs connected to the distribution network can not only solve the problem of power storage, but also make a positive impact on the reactive power regulation of the distribution network. Therefore, to accelerate the transformation of energy structure and improve the overall reactive power regulation ability of the network, EVs may be used to optimize the reactive power regulation [5].
Reactive power optimization is a complex multiobjective optimization problem (MOP) with nonlinear, non-convex and discrete optimization variables [6] considering EVs and new energy types. In general, the traditional interior point or Newton method would be employed to optimize such kind of reactive power regulation model. The idea of these methods is to transform the MOP into a single objective optimization problem with a series of linear weights. However, it is impossible to give a relatively optimal control scheme for an MOP [7][8][9][10]. In contrast, it is more convenient to solve a complex multi-objective reactive power optimization problem with a multi-objective evolutionary algorithm (MOEA) based on Pareto front (PF). Also, the results of an MOEA are easier to understand.
Research in different fields has different approaches to reactive power optimization, and there is a variety of studies, especially in modeling. In [11], a three-stage relaxation-weightsum-correction based probabilistic reactive power optimization model is proposed, while [12] employs a mixed integer convex programming model in reactive power optimization. Both are utilized to address the uncertainties of wind power output. The pioneering work in [13] explores a cold-start linear branch flow model named modified DistFlow to obtain results with better accuracy, whereas [14] presents a novel comprehensive planning framework and optimizes from the perspective of technical, economic, and environmental benefits by considering distributed generation (DG) and EVs.
These different approaches have led to different methods being applied in reactive power optimization. In [15], a fully distributed reactive power optimization algorithm is shown to obtain the global optimum solution of nonconvex problems easily, while [16] employs a decentralized voltage regulation method aimed at mitigating real-time nodal voltage variations. An improved particle swarm optimization algorithm is presented to minimize the sum of power consumption and investment cost of new reactive power compensation equipment in [17]. In [18], a multi-objective dimension-based firefly algorithm is applied to optimize power losses, emission and cost in power systems. An efficient optimization method named hybrid multi-objective moth-flame optimization is proposed to solve the established MOP in [19]. To seek the optimal solutions for minimizing voltage deviation and network power losses, a meta-heuristic algorithm is proposed in [20]. However, there is still a lack of optimization on reactive power regulation when EVs are connected to distribution network. This paper presents a reactive power optimization model where there is high penetration of wind, photovoltaic (PV) and EVs. Because of the rapid increase of the utilization rate of new energy and EVs, their impacts on reactive power optimization in a distribution network have reached a level that cannot be ignored. Most researchers are focused on single factor affecting the distribution network rather than multiple factors. This focus may lead to results not representing reality. Therefore, in this paper, the recent Dual-Population-Based Evolutionary Algorithm for Constrained Multi-objective Algorithm (c-DPEA) is applied to tackle the established constrained MOP [21] and compared to other excellent MOEAs. First, on the basis of reactive power optimization considering a single objective, we present a compromise solution with MOEAs on MOPs. Then, an improved decision-making method, namely an improved technique for order preference by similarity to an ideal solution (TOP-SIS), is proposed. It adjusts the importance of each objective according to the environmental impact in practical application. As a result, different objectives can be set with different importance weights to get more realistic solutions.
The rest of this paper is organized as follows. In Sect. 2, an advanced EV model is added to the traditional distribution network with new energy, while Sect. 3 introduces the optimization principle of c-DPEA and other algorithms, and presents the specific reactive power optimization application design process. The various standard node systems are used to simulate and test the distribution network connected with new energy and EVs in Sect. 4. Finally, Sect. 5 summarizes the work.

Reactive power regulation model of wind power
The wind turbine model based on a doubly-fed induction generator (DFIG) is modified from [22], and shown in Fig. 1. DFIG converts the mechanical power captured by the wind turbine into electrical power. Meanwhile, it is fed into the grid by the DFIG stator and rotor through a back-to-back converter system, which is comprised of a rotor side converter and a grid side converter. According to the Betz equation, the input mechanical power P m is directly related to the wind speed, as: where ρ is the air density, R means the radius of the wind turbine, V refers to the wind speed, C P is a value related to blade tip speed ratio, λ, and pitch angle, β. λ can be expressed as: where ω m represents the rotation speed of the wind turbine.
In the ideal case, the active power can be considered as the maximum power point for tracking control (MPPT) [23], which is provided as follows.
where P base w and v base w are the rated output and wind speed of the wind turbine, respectively. However, the total reactive power fed into the grid is only related to the stator and the grid side converter rather than the rotor [24]. Thus, the total reactive power is given by: where Q g,max and Q g,min mean the upper and lower bounds of the reactive power injected into the power grid, respectively. Similarly Q s,max and Q s,min refer to the upper and lower bounds of the reactive power regulation range of the wind turbine stator side, while Q c,max and Q c,min are the upper and lower bounds of the reactive power regulation range of the grid side converter.
DFIG capacity limits are obtained by considering the stator and rotor maximum allowable currents, referred to as the rated stator and rotor currents in steady state, whereas the generator total capacity limits are obtained by adding rotor power to the stator power. Thus, the reactive power regulation range of the stator side is affected by the maximum current constraints of the stator and rotor sides, as shown in (5)- (7).
where Q s1,max and Q s1,min refer to the upper and lower bounds of the reactive power regulation range of the wind turbine stator side under the maximum current constraint of the rotor side. Q s2,max and Q s2,min are the upper and lower bounds of the reactive power regulation range of the wind turbine stator side under the maximum current constraint of the stator side. L s and L m mean the stator inductance and excitation inductance, respectively. I r,max and I s,max are the maximum current values of the rotor side and stator side, s is the rotor slip, ω 1 is the angular velocity of synchronous rotation, and U s is the stator voltage.
The reactive power of the grid side converter is limited by its capacity and the transmitted active power: where S c,max is the capacity of the grid side converter.

Reactive power regulation model of PV
The flow of PV generation is shown in Fig. 2. Small-scale PV is related to some environmental factors, such as the current solar irradiation and temperature [25]. The active power output P pv is considered as MPPT in the ideal state. Thus, P pv can be presented as: where P base pv refers to the rated total power of PV connected to the distribution network, α pv is the temperature The electric energy generated by PV is usually direct current (DC), and it is thus necessary to use a PV inverter to convert it into alternating current (AC) with the same frequency as the distribution network. Therefore, the influence of inverter capacity and active power output should be considered when calculating the regulation range of PV reactive power output, shown as: where Q pv,max and Q pv,min refer to the upper and lower bounds of reactive power regulation range of the PV equipment. S pv is the capacity of the PV inverter.

Reactive power regulation model of EVs
EVS can be regarded as small energy storage equipment. In addition, the converters of EV charging piles have good flexibility in regulating active and reactive power. Therefore, the addition of EVs can adapt to the development of reactive power optimization in a distribution network. While establishing a two-way charging relationship with the distribution network, EVs can provide significant reactive power regulation. A simplified working principle of the EVs is shown in Fig. 3. The electric energy in the power grid is three-phase AC. Thus, there are converter circuits in the EV charging piles to convert AC into DC, and then provide DC electric energy to the high-frequency DC-DC power converters. The battery energy management system will further optimize the DC input to the battery. When the EVs are connected to the charging piles, they can act as generators in the distribution network if the remaining energies exceed a certain threshold. Otherwise, the distribution network will transmit energies to charge the EVs [26]. The relationship between active power and converter parameters can be obtained as follows: where P car is the input or output active power of the EV, and V s and V c are the grid and charging piles voltages. δ is the phase difference between V s and V c , ω is the angular frequency of the AC system, and L c is the total interface inductance, including the line inductance and the charger filter inductance between the charger and the grid. The reactive power regulation range of EVs is similar to PV, shown as: where Q car,max and Q car,min represent the upper and lower bounds of reactive power regulation range of EVs. S car is the capacity of the converters in the charging piles.

Objective function
There are three objectives to be optimized in this research. Line loss minimization belongs to the economic index of the power grid, while voltage deviation [27] minimization and static voltage stability margin maximization [28] are the indices affecting power grid performance and security. However, there is a hidden contradiction between the objectives. In order to obtain a PF with excellent convergence and diversity, an MOEA is used to handle the three objectives, and a Pareto optimal set is selected under certain circumstances. Particularly, in order to minimize the three objectives simultaneously, it is necessary to convert the maximized static voltage stability margin first. One of the methods is to convert it into solving the reciprocal of the minimum singular value of the minimum convergent power flow Jacobian matrix. In short, the minimization of function f 3 is equivalent to the maximization of the static voltage stability margin. The equations of the three objective functions are presented as follows. where f 1 , f 2 and f 3 mean the power losses, voltage deviation and singular value reciprocal of the Jacobian matrix of the system, respectively. V i and V j refer to the voltage amplitudes of the ith and jth nodes. θ ij is the phase angle difference between the ith and jth nodes, g ij presents the admittance between the ith and jth nodes, N i and N L mean the total node set and all branch set, V * j is the rated voltage of the jth node, and δ min represents the minimum singular value of the system Jacobian matrix.

Constraint condition
IN the process of reactive power optimization, the scope of each parameter and condition should be taken into consideration [29].
(1) Power flow equality constraints: where P Gi and Q Gi represent the active power and reactive power of the ith generation node, respectively. P Di and Q Di are the active power and reactive power demand of the ith node. b ij refers to the susceptance between the ith and jth nodes, N 0 presents the node set except the balance node, and N PQ is the PQ node set.
(2) Generator constraints: where Q min Gi and Q max Gi are the lower and upper limits of reactive power regulation of the ith generator. Q Gi refers to the reactive power currently injected into the grid by the ith generator and V min Gi and V max Gi are the lower and upper limits of output voltage of the ith generator. V Gi is the current output voltage of the ith generator, and N G is the generator set.
(3) Generation ramp constraint (GRC): In addition to meeting the constraints of reactive power output and voltage, the generation units also need to meet certain constraints on regulation active power output [30,31].
where P rate j means the maximum ramp rate of the jth generation unit, �P out j (c) denotes the regulation power output of the jth generation unit at the cth control interval, and T refers to the time cycle of generation units dispatch.
(4) Start-stop constraint: where O means the set of operating generation units at the cth control interval, and R j (c) is the downward ramp rate of the jth generation unit. �P min j (c) represents the minimum active output of the jth generation unit in the cth control interval, and P D (c + 1) means the total active power demand of the power grid at the (c + 1)th interval.

Brief description of the optimization process of c-DPEA, CCMO and WOF
In general, constrained multi-objective optimization problems (CMOP) are more challenging than unconstrained MOPs [32][33][34][35]. In the constrained environment, some search spaces will be directly excluded by the constraint factors, and the algorithm can only search for optimization over a narrow feasible region, or even make all search spaces infeasible in extreme cases, resulting in a final failure to get the optimization results [21,36].
In order to address this kind of CMOP, c-DPEA is applied for reactive power optimization. The main innovation of this algorithm is that two populations are used to co-evolve for the problem. They deal with the infeasible solutions in different ways. Population 1 can find better points in infeasible regions more conveniently and quickly, while the optimization process of population 2 can approach the front faster. This kind of co-evolution obtains a win-win result, as it emphasizes convergence and diversity in two different populations.
A co-evolutionary constrained multi-objective optimization framework (CCMO) also proposes the idea of two populations [37]. This framework evolves one population to solve the original CMOP and evolves another to solve a helper problem derived from the original. It is worth noting that the two populations will share their useful information by their own optimizer separately.
A weighted optimization framework (WOF) is employed to deal with this problem here [38]. A WOF can be seen as a generic method that can be used with any population-based meta-heuristic algorithm. This algorithm is intended to solve multi-objective optimization problems with a large number of decision variables.

Application design of algorithms in reactive power optimization of a distribution network
(1) How to handle variables, especially the discrete ones: Reactive power optimization includes continuous variables and discrete variables. These need to be discussed separately in iteration. Continuous variables can be iterated according to normal optimization, while discrete variables need to be rounded by continuous spatial values. In addition, for the values of the discrete variables, the upper and lower bounds of the continuous spaces correspond to the upper and lower bounds of discrete variables. For example, the decision variables are set to 17 in the 33-bus system, including several new energy outputs, charging and battery change stations, reactive power compensation devices and a group of transformer taps. Then, the upper and lower bounds of reactive power can be calculated from Sect. 2.
(2) How to evaluate the fitness function: It is necessary to consider the objectives and constraints in reactive power optimization in an actual distribution network. In the process of optimization iteration, the individuals must satisfy the condition provided in (13) to calculate the power flow. Consequently, the fitness function can be designed as: where T is the set of objective functions, and t represents the tth objective value in the set of objective functions. f t x i and f fit,t x i are the values of the objective and fitness functions, respectively. η is a penalty coefficient which is often set as a larger constant, and q is the number of objectives that do not meet the constraints.
For example, the reactive power optimization flow chart of c-DPEA for the model is shown in Fig. 4. The power flow to obtain the solution set of the objective function value needs to be calculated first. Then the method of adjusting the fitness function (18) is employed. Finally, the algorithm enters the next iteration.

Application for improved TOPSIS in comparing optimization results of different algorithm
In some real situations, the importance of different objectives will change with different factors. Similarly, different algorithms may have different optimization results because of the change of real parameters. We propose an improved adaptive Technique for Order Preference by Similarity to an Ideal Solution (TOPSIS) decision method to compare the optimization results of the different algorithms in different cases. Line loss and voltage deviation can be the minimum index in the decision method. However, the static voltage stability margin should be the maximum index. The maximum index needs to be transformed into the minimum index by the method of finding the reciprocal. If there are n objects to be evaluated and m indices to be evaluated, an evaluation matrix can be obtained, as: The general calculation of the decision method is as follows: where Z + means the maximum value and Z − is the minimum value. D + i refer to the distance between the ith evaluation object and the maximum value, while D − i is the distance between the ith evaluation object and the minimum value. S i is the non-normalized score of the ith evaluation object.
For the objective function of the three very small indicators mentioned above, the value of the objective function needs to be transformed into the maximum index. Then, TOPSIS can be employed to judge the advantages and disadvantages of the different algorithms. It is worth noting that the result range of S i must satisfy the condition of 0 < S i < 1 . The larger S i is, the closer it is to the objective maximum. A threshold will be set in this process. When the S i obtained from TOPSIS is greater than the threshold and less than 1 after the optimization, the scores of this specific algorithm will be increased by 1. The higher the score is, the higher the matching degree between the algorithm and the optimization model is.
The weights of line loss, voltage deviation and static voltage stability margin are set as w 1 , w 2 and w 3 , respectively. There is a PV coefficient k 1 , which changes with the increase or decrease of solar irradiation or temperature difference. Similarly, the wind power coefficient k 2 varies according to the wind speed. k 1 and k 2 affect the weight w 2 of voltage deviation by a feedback link. At the same time, a line loss threshold will be set. If the optimization result is less than this threshold, the line loss weight will be increased according to a certain proportion. Finally, two or three objective weights are imported to calculate the fitness score of the algorithm.

Simulation model
The topology of the 33-bus system is shown in Fig. 5. In this work, the installed capacities of DFIG and PV are set as 300 kW and 21 kW, respectively, while the total capacity of EV charging and battery change station is 80 kW. The active power output of the wind turbine affected by wind speed can be divided into three stages, with v in w , v out w and v base w set as 3 m/s, 25 m/s and 12 m/s, respectively. The initial solar irradiation is 600 W m −2 in the PV system, the temperature conversion coefficient is 0.2 and the temperature difference can be regarded as a fixed value of 2 °C. New types of energy are affected by different environmental factors, so the initial parameter differences are explained in Tables 1 and 2. EV is a great uncontrollable factor, and can be regarded as an uncertain load or a small battery. Here we have 3 and 6 EV charging and battery change stations in the 33-bus system and 118-bus system, respectively. However, the number of EVs in each station is set to 10 and the power of a single EV is 0-10 kW, so the total power of EVs at each station will always be within the range of 0-100 kW. Then, the state of EV batteries will be determined by its own state of charge (SOC). The EVs are in charging state within the SOC of The default population and iteration are 50 and 50, respectively, so the optimization results can be obtained quickly. However, it is necessary to increase the respective population and iteration to 200 and 50 in the three-objective optimization, so that comprehensive results can be obtained and the error can be reduced to some extent. In this work, the three algorithms are independently applied in the same software environment.
The optimization variables include reactive power compensation devices, transformer taps, new energy generators, EV charging and battery change stations, while the variations of wind speed and solar irradiation are only set to several regular stages. More practical problems will be gradually realized in subsequent research.

Analysis of optimization results
The PF obtained by different algorithms for bi-objective and three-objective in the IEEE 33-bus system are shown in Figs. 6 and 7, respectively. It is seen from Fig. 6 that CCMO and c-DPEA can find better PF in the optimization process. By comparison, the optimization degree of WOF is lower than the other two. The PF of WOF is much narrower than those of CCMO and c-DPEA. WOF has the most dominated solutions in several simulation experiments. From Table 3, the HV and spacing also present the outstanding uniformity of c-DPEA and poor uniformity of WOF. For CCMO, it has the characteristic of being mediocre. Its HV and spacing are inferior to c-DPEA although it can obtain good PF. In the threeobjective optimization in Fig. 7, the optimized Pareto surface (PS) of CCMO is concentrated in a small area in the middle while WOF has the advantage of uniform distribution. By contrast, c-DPEA shows the advantages of fast optimization speed, uniform PS and great optimization degree. From Table 4, the indices of c-DPEA cover almost all the optimal values for both the 33-bus and 118-bus systems. In this regard, c-DPEA gives satisfactory results with excellent equilibrium and stability. Its results are moderate while being optimized as much as possible.
Here, 0.7 is set as the threshold and the initial score of each algorithm is 0. Fifty individuals in the optimized population are randomly selected to calculate the fitness degree. If the fitness degree is greater than the threshold, the score is increased by 1. In order to reduce contingency and error, the optimization results of four different systems are added to calculate the fit. All data in the tables are average values. Tables 5 and 6 shows the fit scores of bi-objective and three-objective, respectively. In Table 5, c-DPEA shows the advantages in a small distribution network while WOF shows its advantages in the 118-bus system. Although CCMO has excellent scalability in PF, its score in practical model is not satisfactory. In Table 6, the average fit scores of c-DPEA are slightly higher than those of the other two algorithms. However, CCMO has the best performance in the 39-bus system. These data directly prove that c-DPEA is still more stable and suitable in different situations than the other two algorithms when the number of objectives is different.
In Tables 7 and 8, the work is divided into three different cases. The first case shows the influence of new energy and EVs while the last two consider only one of them. In longitudinal comparison, the indices of c-DPEA  largely become the optimal solutions in the first case for both the 33-bus and 118-bus systems, as C-DPEA has better comprehensive capability and evenness. In other words, this algorithm can exert its excellent optimization ability in the face of these situations. From the    perspective of horizontal comparison, the case with only EVs has good performance in spacing. This is because the addition of EVs may reduce the reactive power optimization pressure of a distribution network. Thus, some losses of the distribution network will be greatly reduced with new energy and Evs. This directly affects the economy and security of the distribution network in practice. In Figs. 8 and 9, the effects of different wind speeds and solar irradiation in the optimization results are shown.    The simulation is carried out when other variables are consistent. In Fig. 8, the three objectives decrease significantly with the increase of wind speed, especially in the range of 5-13 m/s. The change trend of the three objectives is not obvious for the increase of solar irradiation in Fig. 9. However, the three objectives still have an overall downward trend. Their variations occur around the average values, but the deviation rates are not high. Thus, c-DPEA is fully applicable to this distribution network model. Although there are various influencing factors such as climate in real life, this algorithm can get the results stably and quickly. Figures 10 and 11 compare the results with and without the reactive power compensation devices in the distribution network. The initial conditions of the two situations are the same based on c-DPEA. In order to draw a more complete and clearer graph, 5000 individuals and 50 iterations are set in the two different situations. Figure 10 is obtained by ignoring the reactive power compensation   Fig. 11 Optimization results of reactive power compensation device connected to distribution network by c-DPEA device, and its distribution approximates a smooth surface while the solutions are evenly distributed. For the optimization results with a reactive power compensation device shown in Fig. 11, the distribution change of the solution is small and is concentrated in a narrow long band compared with the former. The data show that the optimal values of line loss and voltage deviation after adding a reactive power compensation device are slightly higher than those without the device. In contrast, the optimal results of the latter are obviously better than the former in solving the maximum value of the static voltage stability margin. Consequently, c-DPEA can be used in general situations in the multi-objective optimization task of a real distribution network. It can obtain different optimization priority results in the face of optimization requirements in different situations.

Conclusion
To show the adaptability of MOEAs to a distribution network with new energy and EVs, this paper constructs a large framework, summarized as follows: (1) The distribution network model with new energy and EVs is built to reflect the reactive power regulation potential in various types of distribution network. (2) Several MOEAs are applied to the above model.
Moreover, their operational process and applicability to specific models are described in detail.
Comparing with different algorithms, both CCMO and WOF are slightly inferior to c-DPEA from different indices. c-DPEA has the characteristic of stable adaptation in optimization problems among different environments. The comprehensive capacity of c-DPEA and the evenness of its PF are shown in Sect. 4. However, its advantages in large systems are not as significant as those in small systems, especially in the 39-bus system. It also explains that there are a certain number of dominated solutions in PF for c-DPEA. (3) An adaptive strategy is added to the improved TOPSIS method with the precision of a distribution network model. The strategy can change the priority of an objective according to its actual needs during the optimization, with the more important objectives receiving more attention. (4) This research also attempts to examine and verify whether the proposed model is affected by various environmental factors, e.g., the change of wind speed and solar irradiation, or the addition and deletion of reactive power compensation devices. Taking c-DPEA as an example, it shows that the model has a good ability to adapt to the actual changeable environment.
In summary, the rational use of the proposed method can provide an alternative idea for the multi-objective optimization of a distribution network in the future.