Planning of High Renewable-Penetrated Distribution Systems Considering Complementarity and Cluster Partitioning

Photovoltaic (PV) and wind power (WT) resources can influence each other in some scenarios, and this influence tends to show that the rise of PV resources may indicate the drop of WT resources, and vice versa. This pattern of PV and WT resources influencing each other is called the complementary characteristics of PV and WT power. The complementary characteristics of the power outputs of different kinds of distributed renewable energy resources (DRERs) and the correlation between DRERs outputs and loads can impact the consumption of DRERs by the loads within the grid, which represents the rate of DRER outputs consumed by loads instead of being reduced. In this regard, this paper investigates a planning strategy for DRERs considering these two factors. An improved co-variance matrix method is applied to generate complementary samples of DRERs and correlated samples of DRERs and loads. The samples generated are used to study the impacts of the degree of correlation between DRERs and loads on the consumption ability of DRERs. The concept of the cluster is introduced as a region including DRERs with complementary characteristics. Based on the cluster partition method and the samples generated, the DRERs planning model is proposed to maximize the profits of different DRER stakeholders. The planning model is transformed into a single objective model through the ideal point method. A Benders decomposition-based method is developed to efficiently solve the proposed model, and an actual network in China is used to illustrate its performance. The results show DRER consumption can be significantly improved by the proposed planning model.


Introduction
With the promotion of renewable energy policies in China, more and more distributed renewable energy resources (DRERs) are being connected to power systems. This has led to occurrences of reverse power flows in distribution networks, where energy flows between substations with the same voltage level or from low-voltage substations to an upstream substation. Excessive reverse power flows are usually not allowed in traditional planning strategies of DRERs because they can deteriorate the power quality and result in increased power losses, even breakdown of conventional transformers. Therefore, DRER outputs are reduced or large capacity storage devices are installed to avoid the occurrence of reverse power flows [1]. Recently, some transformers that allow reverse currents have been considered to improve the penetration rate of DRERs [2]. Therefore, it is also possible to limit reverse power flows within a specific region, which includes several substations to increase the overall outputs of DRERs. In this work, the region in the distribution network where reverse power flows are allowed is called a cluster and the process to determine the region is called cluster partitioning.
The consumption of DRERs is affected by two factors: (1) the complementary characteristics of photovoltaic (PV) generation and wind turbine (WT) generation; and (2) the matching characteristics between the power generation and loads. In order to study how these two factors affect the consumption of DRERs, the generation of the data that can represents the complementary characteristics of PV and WT generation and the matching characteristics between power generation and loads using the expected correlation coefficient are essential. The generated two sets of data that are correlated with each other and show the correlation relationship described using correlation coefficients is called correlated samples.
However, there are two obstacles to describe the correlation relationships among PV, WT and load power. One of the difficulties is that the outputs of PVs and WTs are commonly assumed to obey fixed probability distribution functions respectively, which leads to incorrect correlated samples generated by the traditional methods [3][4][5]. The other problem is that it is difficult to choose the suitable parameters to describe the output characteristics of PV and WT generations. The reason for the difficulty of choosing the suitable parameters is that there are no rules to follow of picking the most suitable parameters for the hypothetical distribution functions, and that the hypothetical parameters may not apply to the partial historical data of PV and WT generations due to the randomness of wind speed and solar radiation. Many researchers have discussed the generation of correlated samples. In [3], a traditional correlation coefficient matrix is used to generate the correlated wind speed. This work assumes the wind speed is complied with the normal distribution function, which results in a large forecast error in wind power outputs. A correlated model is established combining the correlation coefficients with the time series in [4], but the assumption remains that the DRERs follow the given probability distribution functions. In [5], a transformation method for the correlated wind speed non-normal distributed based on the Cholesky decomposition is proposed. However, the accuracy of the proposed method is not guaranteed if the distribution functions of WT generation are unknown. Furthermore, the correlated coefficient would be changed during the transformation process. In order to generate the correlated samples satisfying unknown distributed function with expected correlation coefficient, a method of sample generation is proposed in this paper. The proposed method is based on the orthogonal transformation and marginal transformation proposed in [5] as well as the kernel density estimation method in [6], which is an efficient way to directly establish the distribution functions from characteristics of historical data [7].
Many researchers have also addressed the DRER planning problem in recent years. The main goal is to determine the optimal location, capacity, and type of DRERs. The objective function of [8,9] aims to minimize the total network cost. The objective function of [10] aims to minimize distributed generation (DG) installation cost and network losses. In [11], a DG planning scheme is proposed to minimize the network upgrade cost, network loss cost, power interruption cost, and users' power supply cost. A compromised scheme is obtained in the non-inferior solutions. In [12], a planning model considering the profits of both DG owners and distribution companies are established. The investments and operational costs of DGs as well as the power purchase and interruption costs of distribution companies are minimized in the planning model, and the incomes of DG owners are maximized. This model is established from the view of economy without considering the optimal power losses and voltage profile during the power delivery. In [13], several planning models of DG placement are introduced considering the objectives of the minimum costs and the maximum benefits and loadability with constraints of power flow and power quality. It's pointed out in [13] that the active network management is essential for the access of high penetration of DG, which involves the reactive and active power control of DGs. According to the above analysis, the previous planning models are usually proposed from the perspectives of economy and optimal power flows. Some aspects are not considered in the previous studies. Firstly, the complementary characteristics of PV generation and WT generation as well as its impact on the planning results are seldom considered in the previous works. Secondly, the previous planning models are built based on one substation and its low voltage networks, therefore the DG installation capacity in such networks cannot excess its maximum load demands. In fact, reverse currents are allowed in some transformers, therefore, the proper power interaction between different substations and the output characteristics of different kinds of DGs should be considered in the planning model. This factor is considered in this paper, along with the complementarity in power supply among substations. To ensure the maximum consumption of DRERs, a planning model is proposed in this paper not only from the perspectives of economy and optimal power flow which are commonly considered in the previous works, but also with the consideration of active power and reactive power control of DRERs based on the actual output characteristics of DRERs and the power interactions between substations.
The planning strategy for DRERs in distribution power systems is usually a non-convex complex optimization problem due to the power flow equations. Intelligent algorithms are widely applied to solve the distribution system planning problem, such as the improved particle swarm optimization [14], genetic algorithm [15], ant colony algorithm [16], and Tabu search method [17]. The drawbacks of these algorithms include a local optimal solution and long computation times. To make the problem tractable, some researchers reformulated the original non-convex model into a convex model based on second-order cone programming (SOCP) [18]. Reference [19] shows the SOCP model is strictly accurate when the objective function is convex and increases with the branch current. To improve the computational efficiency, the complex DRER planning strategy is modeled based on SOCP and decomposed by the Benders decomposition method [20,21].
The planning model is formulated in this paper to maximize the profits of different stakeholders and the consumption ability of the DRERs. To deal with the programming problem with several objectives, some researchers combine Pareto's multi-objective optimization with modern intelligent algorithms, such as the particle swarm optimization algorithm [22], improved genetic algorithm [23], and evolutionary algorithm [24]. However, this method usually results in a high computational burden to obtain the Pareto solutions. Other researchers established an addable objective by unifying the dimensions of different objectives [25]. The ideal point method, which is applied in this work, can unify the dimensions of different objectives and does not increase the complexity of the original planning problem [26].
A DRER planning model is proposed in this work to improve the consumption ability of the DRERs in a large-scale distribution network, taking the complementarity between the DRERs and correlation characteristics between DRERs and loads into account. First, a sample generation method for correlated WTs, PVs, and loads that follow arbitrary distribution functions is proposed. Then, a planning model for the DRERs based on the clusters is established. The objective function of the planning model is established by through the application of the ideal point method. Benders decomposition is introduced to efficiently solve the proposed model. Finally, an actual distribution system demonstrates the effectiveness of the proposed approach. The contributions of the paper are as follows: (1) A correlated sample generation method is proposed to generate correlated PV, WT, and load samples following any probability distribution functions. The proposed method can guarantee that the generated samples satisfy the expected correlation coefficient. (2) The concept of a cluster is introduced to confine the reverse power flows within a region by utilizing the complementary characteristics of PVs and WTs to improve DRER consumption. The remainder of this paper is organized as follows. The generation of correlated WT, PV, and load samples is presented in Section 2. The DRER planning model based on the cluster partition is shown in Section 3. The Benders decomposition model is presented in Section 4. The simulation results are presented in Section 5. Section 6 concludes the paper.

Generation of Correlated WT, PV, and Load Samples
In order to study how the complementary characteristics of PV and WT generation and the matching characteristics between the DRER generation and loads affect the consumption of DRERs, a method for generating correlated samples is proposed in this section. An improved correlation matrix method is introduced to describe the correlation among PV generations, WT generations, and loads. Combined with the orthogonal transformation, a non-parameter estimation model is proposed to describe the probability distribution of the PV/WT generations and loads.

Method to Generate Correlated Samples
Correlated samples complying with the identical probability distribution function can be generated through a covariance matrix. The covariance matrix is established based on the Pearson correlation coefficient to describe the correlation of random variables that conform to the standard normal distribution functions, which is shown as: where C represents the covariance matrix; γ p,n represents the Pearson correlation coefficient between the pth and nth elements of the data vector. Equation (1) is used to generate a lower triangular matrix L with the Cholesky decomposition method, and the correlated sample H = [H 1 , H 2 , . . . , H n ] T is standard normal distributed function is obtained by (2), where X = [x 1 , x 2 , . . . , x n ] T is an independently and randomly generated vector satisfying the standard normal distribution function. The elements of vector H are normally distributed, so it is necessary to transform the elements of H to the objective sample X = [x 1 , x 2 , . . . , x n ] T [4], which follows the arbitrary probability distribution of the PV/WT generation or load. The expression is given by (3), where f i is the ith cumulative probability distribution function of PV/WT generation or loads; ψ is the cumulative probability function of the standard normal distribution; In this paper, kernel density estimation [5] is utilized to describe the scatterplot of the ith cumulative probability distribution function f i in (3). The cubic spline interpolation method [27] is used to obtain the expression of the cumulative probability functions.

Process to Obtain the Expected Samples
A nonlinear transformation process is used to obtain the expected correlated samples. The Pearson correlation coefficients in (1) are only suitable for data that follow the identical probability distributions. Therefore, the Spearman rank correlation coefficient [28] shown as Equation (4) is applied to describe the correlation of the objective samples obeying different probability distributions in this work.
where η is the Spearman correlation coefficient; d i represents the rank of data i; n is the total amount of data. Spearman rank correlation coefficients ranging from −1 to 0 are used to describe the complementary characteristics between PV and WT outputs, while those ranging from 0 to 1 are used to describe the degree of correlation between the DRERs and the loads.
Step 1: Input a series of Pearson correlation coefficients γ 1 , γ 2 , . . . , γ n with a range of [−1,1]; Step 2: Each γ i is taken as the value of all the elements of γ p,n within the matrix of Equation (1) to build the covariance matrix C i . Then generate a lower triangular matrix L i using the Cholesky decomposition method; Step 3: Randomly generate the vector X and generate the objective sample vector X i using the method proposed in Section 2.1. Then split X i in half. The first half of X i is X 1,i representing one of the two correlated samples, and the second half is X 2,i representing the other of the two correlated samples.
Step 4: Sort the sequences X 1,i and X 2,i to obtain the rank of all the elements of X 1,i and X 2,i , respectively, then calculate the difference between the two ranks of each element and the Spearman rank correlation coefficient η i .
Step 6: Use the cubic spline interpolation to interpolate the scatters to form the curve of the relationship between the Spearman and Pearson correlation coefficients.
Step 7: Set the expected Spearman correlation coefficient for the expected samples and calculate the corresponding Pearson correlation coefficient based on the curve obtained in Step 6.
Step 8: Add the Pearson correlation coefficient obtained in Step 7 into the C i in Step 2 and redo Steps 2 and 3 to obtain the objective samples satisfying the expected correlation coefficient.
The generated PV and WT outputs with different Spearman correlation coefficient are used as the resource scenarios for the DRER planning model to study the impact of complementary characteristics of DRER on the planning results.

Cluster Partitioning
A cluster in this paper is defined as a set of substations and related lines that form a region within the distribution network. To deal with reverse power flow due to highly penetrated DRERs, this paper divides the distribution network into several clusters, and allows reverse power flow among the substations within a cluster. The clusters are set as the planning units of DRERs within which the complementary outputs of different type of DRERs are used to improve the consumption of hourly outputs. It should be pointed out that the reverse power flows tend to occur simultaneously when all of the substations within a cluster are installed with the same kind of DRERs. In this case, renewable energy consumption in the cluster cannot greatly be improved with the same kind of DRERs within clusters. However, different kinds of DRERs with complementary characteristics installed within a cluster can be utilized to prevent the reverse power of different transformer stations from existing simultaneously. In this way, the matching characteristics between the DRER outputs and the loads within a cluster can be improved.
The above considerations mean that the electrical distance of the nodes included in a cluster should be limited to a relatively suitable range to realize the power interaction. The cluster partition index is used to guarantee the cluster partition results meet the expected standards. Therefore, the cluster partition index proposed in [29], which represents the electrical distance within a cluster, is used in this work. Moreover, the consumption ability of the DRERs installed within a cluster is expected to reach a maximum, which is represented by the index of power balance within each cluster. Therefore, the cluster partition index A is shown as in (5), where A represents the cluster performance index; ρ 1 represents the index represents the electrical distances between all the nodes within the clusters; ρ 2 represents the index represents the average power balance within all the clusters; S represents the set of nodes within the system e ij represents the weight coefficient between node i and j which can be calculated using the (6); L ij represents the electrical distance between node i and j determined by the feeder length; max(L) represents the maximum value of all the L ij ; υ i represents the weight of node i which can be calculated using (7); φ(i), ϕ(i) represent set of upstream nodes and downstream nodes of node i, respectively; M represents the total weight of all nodes; S supplied,c (t) represents the apparent power of the DRERs within cluster c at time t; S demand,c (t) represents the apparent power of the loads within cluster c at time t; Φ(i, j) is 1 if nodes i and j are in the same cluster and 0 otherwise. The power balance in each cluster is the ratio of DRER generation to the total load within the cluster. The Tabu search method proposed in [29] is introduced to solve the cluster partition. The proposed planning model for the DRERs is based on the results of the cluster partition.
In order to deal with the reverse power flow of each substation, this paper divides the distribution networks into clusters, and allows the 35 kV stations within cluster to have reverse power flow which would flow to recent substations with higher load power. However, this reverse power flow will offset the load power of other substations, and the renewable energy access capacity of the cluster still cannot increase greatly. Nevertheless, this problem can be solved if different power plants within the cluster are connected to different kinds of renewable energy with complementary characteristics.
In order to illustrate the necessity and benefit of cluster partition, the schematic diagram of a cluster is shown in Figure 1 During the period that the reverse power due to PV outputs flows into TS 3 and TS 4, the WT power outputs of TS 3 and TS 4 is low and these two transformer stations still are able to consume the reverse flow. Therefore, reverse power can be consumed by the load of this cluster with the complementary outputs of PV and WT and can improve the PV planning capacities of the TS 1 and TS 2. As for excessive WT power, the reverse power exists at night or during the period that PV outputs are low and can also supply power to the load of TS 1-4 without affect the consumption of PVs. In summary, the use of complementary characteristics between different DRERs within the cluster can improve the ability of local renewable Energies 2019, 12, 2090 7 of 22 energy consumption, but the above-mentioned reverse power should not be able to flow in a wide range of networks, otherwise it is prone to appear problems like more network loss and over-voltage, so it is necessary to limit the region of reverse flow according to the result of cluster partition. PVs. In summary, the use of complementary characteristics between different DRERs within the cluster can improve the ability of local renewable energy consumption, but the above-mentioned reverse power should not be able to flow in a wide range of networks, otherwise it is prone to appear problems like more network loss and over-voltage, so it is necessary to limit the region of reverse flow according to the result of cluster partition.

Objective Function
The objective functions are proposed from both the perspectives of DRER generators and grid operators. As for DRER generators, the investments are expected to reach the minimum, which is shown as (8). Apart from that, the earnings gained from selling energy to grids and governmental subsidies are expected to be as high as possible. This objective is shown as (9), which minimizes the earning loss due to the power reduction through DRER control. The power reduction is expressed by the calculation of the expected DRER outputs minus the actual outputs of DRERs. As for the grid operators, it is expected the accessed DRER can be ultimately consumed by the loads of clusters to improve the power self-supply ability of clusters. This raises the claim to the power match between DRER outputs and loads, and the less hourly power flowed into clusters the better the power match of clusters, which is shown as the objective function (10). It is also an important goal for grid operator to reduce power losses with the access of DRERs, and this objective is shown as (11).

Objective Function
The objective functions are proposed from both the perspectives of DRER generators and grid operators. As for DRER generators, the investments are expected to reach the minimum, which is shown as (8). Apart from that, the earnings gained from selling energy to grids and governmental subsidies are expected to be as high as possible. This objective is shown as (9), which minimizes the earning loss due to the power reduction through DRER control. The power reduction is expressed by the calculation of the expected DRER outputs minus the actual outputs of DRERs. As for the grid operators, it is expected the accessed DRER can be ultimately consumed by the loads of clusters to improve the power self-supply ability of clusters. This raises the claim to the power match between DRER outputs and loads, and the less hourly power flowed into clusters the better the power match of clusters, which is shown as the objective function (10). It is also an important goal for grid operator to reduce power losses with the access of DRERs, and this objective is shown as (11).
where C PV , C OMPV represent the investment cost and operation cost of photovoltaics (PVs) per MW, respectively; S i PV , S i WT represent the capacity of the PVs and WTs at node i; C WT , C OMWT represent the investment cost and operation cost of wind turbines (WTs) per MW; m represents the number of years in the planning horizon; N S represents the number of candidate substations for the DRERs; R represents the bank rate.
where C sell represents the electricity selling price; C sell represents the governmental subsidies for the DRERs; T represents the planning horizon; P i PV (t), P i WT (t) represent the actual active power outputs of the PVs and WTs at node i at time t, respectively;P i PV (t),P i WT (t) represent the ratio of the actual active power outputs of the PVs and WTs at node i at time t to their rated power outputs, respectively; where N c represents the number of clusters; P c (t) represents the power flowing into cluster c at time t; where r ij represents the resistance and reactance of line l i-j between nodes i and j; i ij (t) represents the square of current flowing from node i to node j at time t; The first objective function F 1 minimizes the annual investments and operation costs of the DRERs and F 2 is the benefit losses due to the reduction of DRER outputs. Both F 1 and F 2 represent revenue brought by the DRERs. The third objective function F 3 minimizes the electricity energy flowing into clusters and the fourth objective function F 4 minimizes power losses within the clusters.
The dimensions of the four objectives are different. Therefore, it is unreasonable to directly add them together to form the single objective function. The ideal point method is utilized to transform these objective function values into the deviation values from their ideal values. The ideal values of the four objectives are assumed as F 0 1 , F 0 2 , F 0 3 , and F 0 4 , respectively. The actual values of four objectives are assumed to be bigger or smaller than the ideal values of objectives. The deviation d χ is introduced into the transformed objective function. Therefore, the four objective functions are transformed into the single objective function as shown in (12) using the d χ to represent the deviation distance between the actual value and idea value of the χ th objective function. Due to the dimension difference, the F 0 χ in (12) is used to normalize the deviation distance of different deviation distances.
where d χ is the positive and negative deviations between the actual value and ideal value in the χth objective function, respectively; F χ ,w χ are the χ th objective function and weight coefficient representing its priority, respectively.

Constraints
(1) Constraints on Deviations The deviation variables in (12) should satisfy the equality constraints shown in (14)- (18): Energies 2019, 12, 2090 9 of 22 (2) AC Power Flow Constraints Based on the branch power flow, the AC power flow constraints are relaxed and formulated as per the SOCP model: where P ij (t), Q ij (t) represent the active and reactive power flow from node i to node j at time t, respectively;x ij represents the reactance of line l i-j between nodes i and j; P i L (t),Q i L (t) represent the active and reactive power of loads at node i at time t, respectively; Q i PV (t),Q i WT (t) represent the Reactive power outputs of the PVs and WTs at node i at time t, respectively; v i (t) represents the square of the voltage at node i at time t; v max and v min represent the upper and lower limits of the nodal voltages, respectively; I max represents the maximum value of current; (19) and (20) are the active and reactive power balance constraints; (21) represents the voltage relationships between nodes at two terminals of a connected branch; (22) represents the relaxed branch current expression within the network; (23) represents the minimum and maximum constraints of the nodal voltages; and (24) represents the minimum and maximum constraints of the line currents.
(3) Constraints on the penetration rate of clusters The installed DRER capacities should be limited within a specified threshold. In this paper, the concept of penetration rate is introduced to limit the overall capacities of the DRERs within a cluster. The definition of the penetration rate of a cluster is shown in (25): where S i L represents the apparent power of substation i; ρ c represents the penetration rate of DRERs within a cluster; The penetration rate of any cluster c is constrained to excess a specified value. The constraint is as shown in (26): where ρ c represents the expected value of of the penetration rate of DRERs within a cluster (4) Constraints on the active power outputs of the DRERs Constraints of the reactive power of DRERs that do not participate in the voltage regulation are shown as cos θ = 0.95 Constraints of the reactive power of DRERs that participate in the voltage regulation are (31) and (32).

Establishment and Calculation of the Benders Decomposition Model
The bulleted lists appear as follows: To improve the computational efficiency, the proposed planning model is decomposed into a master problem (MP) and a sub-problem (SP). The master problem is an investment problem, whose decision variables are the installation capacities and the types of DRERs. The sub-problem deals with the optimal power flow problem, and the decision variables include the actual active and reactive power outputs of the DRERs, nodal voltages, currents, and transmission power of each line. The standard form of the proposed planning model is as follows: K T κ y = k κ κ = 1, · · · , N 4 (41) x ≥ 0; where a represents the coefficient matrix of objective function correlated with investment problem; x represents the variables of the installation capacities and the types of DRERs; b represents the coefficient matrix of objective function correlated with optimal power flow problem; y represents the variables including the actual active and reactive power outputs of the DRERs, nodal voltages, currents, and transmission power of each line; U represents the coefficient matrix correlated with variables of penetration rate of clusters of equality constraints; u represents the matrix of expected penetration rate of clusters of equality constraints; E κ represents the coefficient matrix correlated with variables of the installation capacities and the types of DRERs of inequality constraint κ; F κ represents the coefficient matrix correlated with variables of actual active and power outputs of the DRERs of inequality constraint κ; h κ represents the constant of inequality constraint κ; J κ represents the coefficient matrix correlated with variables of actual active and reactive power outputs of the DRERs, nodal voltages, currents, and transmission power of each line of second-order-cone constraint κ; W κ represents the coefficient matrix correlated with variables of nodal voltages and currents of second-order-cone constraint κ; g κ represents the coefficient matrix correlated with variables of the installation capacities of second-order-cone constraint κ; G κ represents the coefficient matrix correlated with variables of actual active and power outputs of the DRERs, nodal voltages, and currents of inequality constraint κ; q κ represents the constraint of nodal voltage, current, and actual active and reactive power output of DRER of inequality constraint κ; K κ represents the coefficient matrix correlated with variables of actual active and power outputs of the DRERs, nodal voltages, positive deviation and negative deviation of objective functions of equality constraint κ; k κ represents the constant including active, reactive load power and ideal value of objective function of equality constant κ; (37) represents (26)-(28); (38) represents (29) (14)- (17) and (19)- (21). The MP is formulated as: The SP is presented as: G T κ y ≤ q κ : δ κ κ = 1, · · · , N 3 (50) The dual form of the SP is expressed as: The steps for solving the Benders decomposition-based model are as follows. Step 1: Solve the MP and obtain the lower boundẐ l of the original problem (36)-(42) and the decision variable x .
where α τ κ represents the relax variable of corresponding with constraint κ of type τ. (b) If the SP is feasible, calculate the upper boundẐ u = a Tx + b Tŷ . Check whether Ẑ u −Ẑ l ≤ ε is satisfied. If so, stop the iteration and output the results. Otherwise, generate the feasibility cut (62), and add it into the MP, then return to Step 1.
Step 3: Solve the RSP, and add the infeasibility cut (63) into the MP, then return to Step 1.
The calculation process of the proposed method is shown in Figure 3.

Case Study
An actual distribution network is applied in this paper to illustrate the efficiency of the proposed model. The network consists of 31 candidate substations. The numerical results are obtained using a laptop computer with a 3.3 GHz Intel Core i5-4590 CPU and 8 GB of memory.

Case Study
An actual distribution network is applied in this paper to illustrate the efficiency of the proposed model. The network consists of 31 candidate substations. The numerical results are obtained using a laptop computer with a 3.3 GHz Intel Core i5-4590 CPU and 8 GB of memory.

Expected Correlated Samples of the PVs and WTs
According to the Sections 2.1 and 2.2, in order to apply the proposed method of generating the expected PV and WT samples into the actual distribution system in China, it is necessary to calculate cumulative probability functions of PV and WT according to the historical data based on the kernel density estimation. The results of cumulative probability functions at different time of the day are shown as Figures 4 and 5.

Case Study
An actual distribution network is applied in this paper to illustrate the efficiency of the proposed model. The network consists of 31 candidate substations. The numerical results are obtained using a laptop computer with a 3.3 GHz Intel Core i5-4590 CPU and 8 GB of memory.

Expected Correlated Samples of the PVs and WTs
According to the Sections 2.1 and 2.2, in order to apply the proposed method of generating the expected PV and WT samples into the  Based on the results shown in Figures 4 and 5, the relationships between the Pearson correlation coefficients used to describe the standard normal distributed variables and the Spearman correlation coefficients used to describe the PV and WT data can be obtained through the steps 1-6. The results of the relationships at two different time periods of the day are shown as Figure 6. It can be seen that the differences of the relationships of correlation coefficients under different time periods of the day are trivial. The expected PV and WT samples can be obtained using the data shown in Figure 6 through Steps 7 and 8 of Section 2.2. Based on the results shown in Figures 4 and 5, the relationships between the Pearson correlation coefficients used to describe the standard normal distributed variables and the Spearman correlation coefficients used to describe the PV and WT data can be obtained through the steps 1-6. The results of the relationships at two different time periods of the day are shown as Figure 6. It can be seen that the differences of the relationships of correlation coefficients under different time periods of the day are trivial. The expected PV and WT samples can be obtained using the data shown in Figure 6 through Steps 7 and 8 of Section 2.2.
correlation coefficients used to describe the standard normal distributed variables and the Spearman correlation coefficients used to describe the PV and WT data can be obtained through the steps 1-6. The results of the relationships at two different time periods of the day are shown as Figure 6. It can be seen that the differences of the relationships of correlation coefficients under different time periods of the day are trivial. The expected PV and WT samples can be obtained using the data shown in Figure 6 through Steps 7 and 8 of Section 2.2. To illustrate the effectiveness of the proposed correlated sample generation method, the comparison of correlation coefficients of generated samples attained by two methods is made. Two cases are defined as follows.
Case 1: Assume that the PV and WT outputs satisfy the Beta and Weibull distribution respectively. The correlated samples of PV and WT outputs are generated using the orthogonal transformation in [5].
Case 2: Correlated samples of PV and WT outputs are generated by the proposed method in this paper.  To illustrate the effectiveness of the proposed correlated sample generation method, the comparison of correlation coefficients of generated samples attained by two methods is made. Two cases are defined as follows.
Case 1: Assume that the PV and WT outputs satisfy the Beta and Weibull distribution respectively. The correlated samples of PV and WT outputs are generated using the orthogonal transformation in [5].
Case 2: Correlated samples of PV and WT outputs are generated by the proposed method in this paper.
The expected PV and WT samples with different complementary characteristics are shown as Figure 7, and the results of these two cases under different correlation coefficients are close to each other. It's seen from Table 1 that the correlation coefficients of samples in Case 2 are closer to the expected correlation coefficients than those in Case 1. Therefore, the effectiveness and the accuracy of the proposed method for generating correlated samples can be verified.  The correlated hourly PV and WT samples in a day are shown in Figure 8, where the increase in absolute value of the negative correlation coefficient between the PV and WT outputs represents the rise in complementary characteristics between the PV and WT outputs.

Cluster Partition Results
The topology and cluster partition results calculated using the method of Section 3.1 are shown in Figure 9, where the light blue blocks represent the different clusters. The distribution network is divided into 8 clusters based on the cluster partition indexes and method in Section 3.1. There are negative power demands at 11 substations, located at nodes 11, 15, 29, 31, 36, 40, 42, 43, 48, 51, and 52. The calculated cluster partition indexes are shown as Table 2. The value of index 2 ρ is calculated based on the power balance of clusters which are also shown in the Table 2. The power balance of the clusters represents the consumption of the DRERs original installed within the clusters. The lower the index of the power balance is, there are more PVs/WTs still can be installed into the cluster, and the power balance of cluster 2, 3, 4, 7, 8 is relatively lower, which represents there are more PV and WT can be installed within these five clusters.  The calculated cluster partition indexes are shown as Table 2. The value of index ρ 2 is calculated based on the power balance of clusters which are also shown in the Table 2. The power balance of the clusters represents the consumption of the DRERs original installed within the clusters. The lower the index of the power balance is, there are more PVs/WTs still can be installed into the cluster, and the power balance of cluster 2, 3, 4, 7, 8 is relatively lower, which represents there are more PV and WT can be installed within these five clusters.  To obtain the planning solution without the cluster partitions, the energy flowing into each cluster in the proposed model is replaced by the energy flowing into the entire power system, and the cluster penetration rate is replaced by the penetration rate of the distribution system. Substations that have negative power demands are not allowed to install DRERs and reverse power flows are not allowed in the system.

Indexes of Cluster Partition
The capacities of PVs and WTs in the planning solution with/without the cluster partitions are shown in Table 3. It's shown in Table 4 that PVs can be installed in the substations with negative power demands in Case 4, while no PVs are installed in the substations with negative power demands in Case 3. Furthermore, the installation of WTs increases in Case 4 compared to the results for Case 3, which could lead to fewer investment costs. It's shown in Table 4 that the investment and operating costs, the power losses of the network, and the energy flow into the system for Case 3 and Case 4. The investment costs in Case 4 are lower by $730,226 compared to Case 3, which indicates that the economic efficiency is enhanced in the planning solution with the cluster partitions. The power losses and energy flow into the system in Case 4 drop by 23.18% and 43.14%, respectively, compared to those in Case 3, which indicates the overall outputs of the DRERs in Case 4 are significantly improved. Therefore, this confirms the proposed planning model taking the cluster partition into account can improve DRER consumption. The voltage characteristics of the DRERs in Case 3 and Case 4 are shown in Table 5. The voltage fluctuation range in the planning solution in Case 4 increases by 2.74% compared to that in Case 3. The power exchanges between substations increase the voltage volatility. The planning solution in Case 1 will therefore be at the expense of voltage fluctuations.

Impacts of the Complementary Characteristics on Planning Results
In order to illustrate the impact of different complementary characteristics of PV and WT outputs on the access planning of DRERs, there are four cases as followed being compared and analyzed.
Case 5: Penetration rates in each cluster are 0.80 and the Spearman correlation coefficient of the PV and WT outputs is −0.95; Case 6: Penetration rates in each cluster are 0.80, and the Spearman correlation coefficient of the PV and WT outputs is −0.65; Case 7: Penetration rates in each cluster are 0.50, and the Spearman correlation coefficient of the PV and WT outputs is −0.95; Case 8: Penetration rates in each cluster are 0.50, and the Spearman correlation coefficient of the PV and WT outputs is −0.65.
These four cases are accomplished in the IEEE 33 bus system shown in Figure 10 and the actual networks shown in Figure 9.
The planning capacities of DRERs in each cluster of Cases 5-8 for the IEEE 33 bus system is shown in Table 6. Figure 11 shows the planning capacities of PVs/WTs at each candidate node in the eight clusters. The bars in dark blue represent the substations installed with the PVs, and the bars in light blue represent the substations installed with the WTs. Four consecutive bars are regarded as a group, and represent the optimal capacity results in Case 5, Case 6, Case 7, and Case 8, respectively. Table 6 and Figure 11 of two different networks show the capacity ratio of the WTs to all of the DRERs in each cluster increases when the complementary characteristics between the PVs and WTs drop. This is because the installation and operating costs of the WTs are lower than those of the PVs. Moreover, WTs can provide active outputs on all days, which lead to less energy flow into the clusters when the complementary characteristics between the PVs' and WTs' drop.
Case 7: Penetration rates in each cluster are 0.50, and the Spearman correlation coefficient of the PV and WT outputs is −0.95; Case 8: Penetration rates in each cluster are 0.50, and the Spearman correlation coefficient of the PV and WT outputs is −0.65.
These four cases are accomplished in the IEEE 33 bus system shown in Figure 10 and the actual networks shown in Figure 9. The planning capacities of DRERs in each cluster of Cases 5-8 for the IEEE 33 bus system is shown in Table 6. Figure 11 shows the planning capacities of PVs/WTs at each candidate node in the eight clusters. The bars in dark blue represent the substations installed with the PVs, and the bars in light blue represent the substations installed with the WTs. Four consecutive bars are regarded as a group, and represent the optimal capacity results in Case 5, Case 6, Case 7, and Case 8, respectively. Table 6 and Figure 11 of two different networks show the capacity ratio of the WTs to all of the DRERs in each cluster increases when the complementary characteristics between the PVs and WTs drop. This is because the installation and operating costs of the WTs are lower than those of the PVs. Moreover, WTs can provide active outputs on all days, which lead to less energy flow into the clusters when the complementary characteristics between the PVs' and WTs' drop.  Figure 10. Topology of IEEE 33 bus system.  In this work, the reduction rate of the PV and WT outputs shown in Table 7 is defined as the predicted PV and WT outputs divided by their actual outputs. Table 7 shows the reduction rates of the PVs and WTs both decline due to the rise of complementary characteristics between PVs and WTs. It's seen from the results of IEEE 33 bus system that the PV reduction rate and the WT reduction rate in Case 5 decrease by 7.89% and 9.37% compared to those in Case 6, respectively. The PV reduction rate and the WT reduction rate in Case 7 decrease by 0.33% and 1.31% compared to those in Case 8, respectively. According to the results of the actual network, the PV reduction rate and the WT reduction rate for a penetration rate of 0.80 in Case 5 decrease by 5.15% and 5.75% compared to those in Case 6, respectively. The PV reduction rate and the WT reduction rate for a penetration rate of 0.50 in Case 7 decrease by 0.03% and 0.88%, respectively, compared to those in Case 8. Thus, the drop in WT reduction rate exceeds the drop in PV reduction rate when the In this work, the reduction rate of the PV and WT outputs shown in Table 7 is defined as the predicted PV and WT outputs divided by their actual outputs. Table 7 shows the reduction rates of the PVs and WTs both decline due to the rise of complementary characteristics between PVs and WTs. It's seen from the results of IEEE 33 bus system that the PV reduction rate and the WT reduction rate in Case 5 decrease by 7.89% and 9.37% compared to those in Case 6, respectively. The PV reduction rate and the WT reduction rate in Case 7 decrease by 0.33% and 1.31% compared to those in Case 8, respectively. According to the results of the actual network, the PV reduction rate and the WT reduction rate for a penetration rate of 0.80 in Case 5 decrease by 5.15% and 5.75% compared to those in Case 6, respectively. The PV reduction rate and the WT reduction rate for a penetration rate of 0.50 in Case 7 decrease by 0.03% and 0.88%, respectively, compared to those in Case 8. Thus, the drop in WT reduction rate exceeds the drop in PV reduction rate when the complementary characteristics between PVs and WTs rise. The voltage fluctuation range of the IEEE 33 bus system in Case 6 is 3.92% more than that in Case 5. The voltage fluctuation range of the actual network in Case 6 is 2.74% more than that in Case 5. The voltage fluctuation range of the IEEE 33 bus system in Case 8 is 3.37% more than that of Case 7. The voltage fluctuation range of the actual network in Case 8 is 1.50% more than that of Case 7. Therefore, the rise in complementary characteristics between PVs and WTs decreases the voltage fluctuation range. This is because the rise in complementary characteristics between the PVs and WTs reduces the probability that PV and WT outputs simultaneously increase or decrease, and thus results in smaller nodal voltage fluctuation. Table 7 also shows the power losses within clusters rise as the penetration rate increases due to the increase of reverse power flows within clusters. The energy flow into clusters decreases as the penetration rate increases, which leads to the increasing power balance of clusters. If the loads in the clusters are fixed, the decreased energy flowing into clusters exceeds the decreased power losses within clusters when the complementary characteristics between PVs and WTs rise. In this way, the DRER outputs increase with the rise of complementary characteristics between PVs and WTs. Therefore, an increase in complementary characteristics between PVs and WTs within clusters helps improve the consumption abilities of DRERs.

Impacts of the Spearman Correlation Coefficients between the DRERs and Loads on the Planning Results
In order to illustrate the impact of the variation of the correlation relationship between DRERs and loads on DRER access planning results, there are four cases in this subsection are compared and analyzed. In this subsection, the Spearman correlation coefficient between PV and WT outputs is −0.85. The following scenarios are compared and analyzed.
Case 9: Penetration rates of the clusters are 0.85, and the Spearman correlation coefficient between the DRERs and loads is 0.90; Case 10: Penetration rates of the clusters are 0.85, and the Spearman correlation coefficient between the DRERs and loads is 0.30; Case 11: Penetration rates of clusters are 0.45, and the Spearman correlation coefficient between the DRERs and loads is 0.90; Case 12: Penetration rates of clusters are 0.45, and the Spearman correlation coefficient between the DRERs and loads is 0.30. Table 8 shows that the increase in correlation coefficient between the DRERs and loads results in a decrease in the reduction rate of DRER outputs. As the correlation coefficient between the DRERs and loads increases, the power losses and the energy flow into clusters increase. Comparing the power losses within clusters and energy flow into clusters shows the decreased energy flow into clusters