As the existing cluster partition index is not comprehensive, based on the DN structure and cluster function, a comprehensive cluster partition index is proposed to complete the cluster partition in this paper. The proposed comprehensive cluster partition index includes the modularity index, power balance index, and nodal size index.

To carry out the cluster partition, a hybrid genetic-simulated annealing (HGSA) algorithm is utilized. The HGSA algorithm uses the annealing selection as the individual replacement strategy, while the global information obtained by the genetic algorithm can be completely used, and the premature convergence of the genetic algorithm can be avoided. Then, the global convergence of the algorithm is enhanced. The cluster partition is implemented as follows:

Step 1: Initial optimization parameters are set: population size n, initial temperature T ⁰, termination temperature T ^end, temperature cooling factor r, maximum number of genetic generations M, and objective function for the cluster partitioning index ϕ .

Step 2: The number of temperature updates is set equal to 0. Considering that the cluster partition is carried out based on the original network connectivity, this paper uses the unweighted adjacency matrix   A i j to represent the connectivity of the network, where A i j = 1 indicates that nodes i and j are connected and A i j = 0 indicates that nodes i and j are not connected. The clustering partition is randomly modified when A i j = 1. The result of the modification is 0 or 1, and 0 means that the two nodes are disconnected from each other, and then, the two nodes belong to different clusters. In the iteration of the algorithm, in order to ensure that the result after crossover still satisfies the meaning of the adjacency matrix, only the upper triangle of the adjacency matrix is chosen to carry out the crossover. After the crossover step, the upper triangle of the matrix is symmetrically transferred to the lower triangle to form the newborn individuals, thus generating the initial population P _l (k).

Step 3: The individuals that satisfy the cluster constraints are screened. Considering that the reverse power flow only occurs within the cluster, the net power within each cluster needs to be greater than 0 at each moment among the cluster partition. For the individuals who cannot satisfy the constraint, the population eliminates these individuals.

Step 4: Eq. (12) is chosen as an indicator to calculate the fitness of individuals. Considering that the value of fitness intuitively reflects the superiority or inferiority of cluster partition, the individuals with greater fitness are replicated to the offspring to form the new population P l s k + 1 .

Step 5: The crossover and mutation are carried out for   P l s k + 1 using the traditional genetic algorithm to obtain the population P l v k + 1 . The mutated populations are limited to accepting bad solutions by the simulated annealing algorithm to form new populations P l z k + 1 .

Step 6: P l k = P l z k + 1 is set. If the genetic algebra accumulates to the maximum number M, then step 7 is repeated; otherwise, step 3 is repeated.

Step 7: The temperature is updated, i.e., T _l =rT ⁰ , k = 0, P l + 1 k + 1 = P l k . If the convergence condition T _l < T _end , then the computation to output the optimal solution is terminated. Otherwise, the temperature reduction operation is performed, i.e., T _l+1 = rT _l , and step 3 is repeated.

In the upper layer, the objective function of the line planning model is established as follows: min ⁡ f = ∑ i , j ∈ Ω s C L ∑ j ∈ φ i D i j x i j χ 1 + χ β 1 + χ β − 1 + C e ∑ t = 1 T ∑ k = 1 N k ∑ i = 1 Ν n , k ∑ j ∈ φ i 1 2 x i j ε i j i i j , t , (13) where Ω s is the node set of the branch network. φ i is the child node set of node i. C ^L is the investment cost parameter of lines. D _ij is the line length of line i–j. x _ij is the 0–1 variable, where x _ij = 1 means that the line is installed and x _ij = 0 means that the line is not installed. χ is the bank rate. β is the payback period. C ^e is the electricity price. N _n,k is the node numbers of cluster k.   ε i j is the resistance of the line i–j. i _ij,t is the squared value of the current in the line i–j at time t. The constraints include a power flow constraint and a penetration rate constraint of PVs in the clusters.

(1) Power flow constraint

The power flow containing line variables is constrained by second-order conic relaxation (SOCR) (Shaker et al., 2023) as follows: ∑ i ∈ Ψ j P i j , t − ε i j i i j , t − x i j P j , t L + η P j , t CH − P j , t DS η + P j , t PV + P j , t PVf = ∑ c ∈ φ j P j c , t , (14) ∑ i ∈ Ψ j Q i j , t − δ i j i i j , t − x i j Q j , t L + Q j , t PV + Q j , t PVf = ∑ c ∈ φ j Q j c , t , (15) ∑ j ∈ φ i x i j = 1 , (16) v j , t = v i , t ′ − 2 ε i j P i j , t + δ i j Q i j , t + ε i j 2 + δ i j 2 i i j , t , (17) 0 ≤ ν i , t ′ ≤ W x i j , (18) W 1 − x i j + ν i , t ≤ ν i , t ′ ≤ ν i , t , (19) 2 P i j , t 2 Q i j , t i i j , t − ν i , t 2 ≤ i i j , t + ν i , t , (20) x i j v min ≤ ν i , t ≤ x i j v max , (21) 0 ≤ i i j , t ≤ x i j I max 2 , (22) where ψ _j is the upstream node set of node j. δ i j is the reactance of line i–j. P _ij,t and Q _ij,t are the active and reactive power through line i–j at time t, respectively. η is the ES charging and discharging efficiency. P j , t CH and P j , t DS are the charging and discharging power of ES at node j at time t, respectively. P j , t L and Q j , t L are the active and reactive power of loads at node j at time t, respectively. P j , t PV is the actual active power generated by the PV at node j at time t. P j , t PVf is the additional active power required to be generated by the PV at node j at time t. φ _j is the downstream node set of node j. Q j , t PV is the actual reactive power generated by the PV inverter at node j at time t. Q j , t PVf is the additional reactive power required to be generated by the PV inverter at node j at time t. v _i,t is the square of the voltage at node i at time t. v i , t ′ is the square of the voltage at node i that is constrained by the line variable x _ij . W is a large constant. v ^min and v ^max are the upper and lower voltage limits, respectively. I ^max is the maximum current limit.

(2) Penetration rate constraint of PVs

Many isolated nodes exist in the DN that need to be connected to the system, and the penetration rate of PVs should be constrained when carrying out line planning to access these isolated nodes, which is defined as follows: ∑ j = 1 Ν n , k S j PV ≤ ∑ j = 1 Ν n , k S j L , (23)

where S j PV is the planed capacity of the PV of node j in the cluster and S j L is the load apparent power of node j in the cluster.

In the lower layer, the PV and ES planning model is established to identify the “worst scenario " of the uncertain variables, and based on the worst scenario, the proposed model minimizes the investment and operation costs of source storage within the cluster. Therefore, the objective function for the PV and ES planning model within the cluster is established as min ⁡ F = F 1 − F 2 + F 3 , (24) F 1 = ∑ j = 1 , j ∈ Π k N n , k C PV S j PV , f χ 1 + χ β 1 + χ β − 1 + C OMPV S j PV + S j PV , f , (25) F 4 = C sell + C sub ∑ t = 1 T ∑ j = 1 , j ∈ Π k N n , k ω t PV S j PV + S j PV , f − P j , t L , (26) F 5 = ∑ j = 1 , j ∈ Π k Ν n , k ρ 1 + ρ r 1 + ρ r − 1 c bat P j batt , (27)

where F ₁ is the annual investment and operation costs of PVs in cluster k. F ₂ is the annual revenue of PVs in cluster k. F ₃ is the annual investment costs of ESs in cluster k. C ^PV is the investment cost parameter of PVs. C _OMPV is the annual fixed maintenance cost parameter of PVs. П _k is the node set of cluster k. S j PV is the original installed PV capacity of node j in cluster k. S j PV , f is the planned PV capacity of node j in cluster k. C ^sell and C ^sub are the feed-in tariff and subsidized tariff of the PVs, respectively. ω t PV is the output limit of PVs per megawatt for a given sunshine condition at time t. ρ is the discount rate. r is the discounted number of years. c ^bat is the investment cost parameter of ESs. P _batt,j is the allocated capacity of ES at node j. The constraints include uncertainty constraints for PVs and loads, PV and ES capacity constraints, and power flow constraints.

(1) Uncertainty constraints for PVs and loads

The box uncertainty set is utilized to characterize the uncertainty range of active load power, reactive load power, and PV outputs. Meanwhile, an uncertainty parameter is utilized that can be set to adjust the conservativeness of the optimal solution. The larger the uncertainty parameter, the more conservative the solution. The specific formula is as follows: U = { u = P j , t L , Q j , t L , P j , t PV T P j , t L ∈ P j , t L , F − Δ P j , t L , P j , t L , F + Δ P j , t L , ∑ j = 1 N a P j , t L − P j , t L , F Δ P j , t L ≤ Γ PL ; Q j , t L ∈ Q j , t L , F − Δ Q j , t L , Q j , t L , F + Δ Q j , t L , ∑ j = 1 N b Q j , t L − Q j , t L , F Δ Q j , t L ≤ Γ QL ; P j , t PV ∈ P j , t PV , F − Δ P j , t PV , P j , t PV , F + Δ P j , t PV , ∑ e = 1 N e P j , t PV − P j , t PV , F Δ P j , t PV ≤ Γ PV ; , (28) where U represents the boxed uncertainty set. P j , t L , F , Q j , t L , F , and P j , t PV , F are the predicted values of active load power, reactive load power, and PV outputs of node j at time t in cluster k, respectively. Δ P j , t L , Δ Q j , t L , and Δ P j , t PV are the fluctuation deviations of active load power, reactive load power, and PV output of node j at time t in cluster k, respectively. N _a , N _b , and N _e are the total number of active load nodes, the total number of reactive load nodes, and the total number of PV nodes in cluster k, respectively. Г ^PL, Г ^QL, and Г ^PV are the uncertain adjustment parameters for active load power, reactive load power, and PV output S, respectively, which are integers from 0 to N _a , 0 to N _b , and 0 to N _e, respectively. The decision-maker can choose a variety of uncertainty regulation parameters to adjust the scheme flexibly; the larger the value of each uncertainty regulation parameter, the more conservative the resulting planning scheme.

(2) PV and ES constraints

where T k ES is the simulated operating time of the ESs in cluster k. t 1 ¯ is the annual PV generation time. soc _j,t is the state of charge of ES at node j in cluster k at time t. soc ^min and soc ^max are the minimum and maximum charge states of ES, respectively. P ^CH,max and P ^DS,max are the maximum power limits for charging and discharging power of ES, respectively. P j L , ⁡ max is the maximum load demand of node j in the cluster.

(3) Power flow constraints

Eqs (14), (15), and (17) are referred to for the power flow constraints.

In this paper, an iterative solution method is utilized to solve the proposed planning models. For the upper layer, the planning model can be expressed by a specific form as follows: min x c T x + α s . t . P x = W α ≥ k T y l B x + C y l ≤ D R x u l * = O y l S u * y l ≤ V G y l ≤ H T y l ∀ l ≤ k , (36b)

where x is the optimization vectors. c is the coefficient matrices corresponding to the objective functions. P , k , B , C , R _x , S* u , G , and H are the coefficient matrices corresponding to the variables under the constraints. α , D , O , V, and W are the constant column vectors. l is the current number of iterations. k is the maximum iteration. y _l is the variable at the lth iteration. u * l is the value of the uncertain variable u in the “worst” scenario after the lth iteration.

The specific form of the lower-layer model is max u ∈ U min y ∈ Ω x , u d T y s . t . B x * + C y ≤ D R u * u = O y S u y ≤ V G y ≤ H T y , (37)

where y is the optimization vectors. Ω x , u is the feasible domain of y under the given x and u . x ^* is the optimization vectors obtained from the upper layer. d is the coefficient matrices corresponding to the objective functions. S _u and R* u are the coefficient matrices corresponding to the variables under the constraints. Given a set of u , the inner min. problem becomes a second-order cone programming problem, which can be transformed into a “max” form by the duality theory: max u ∈ U , γ , ν , π , μ , μ 1 D − B x * T γ + R x u T ν + u T π s . t . C T γ + O T π + S u T ν + ∑ i G i μ i + H i μ 1 i ≤ d μ i 2 ≤ μ 1 i , ∀ i = 1 , 2 , . . . , j γ , π , ν ≥ 0 , u ∈ U , (38) where[γ,ν,π,μ,μ ₁ ] are the dual variables in the lower-layer model.

After the above transformation, the proposed model can be solved by the iterative solution method as follows:

Step 1: Given a set of u values as the initial worst-case scenario, a lower bound is set on the operating cost LB = ‐ ∞ , an upper bound UB = + ∞ , and the number of iterations l = 1.

Step 2: The upper model is solved based on the worst-case scenario u _l to obtain the optimal solution x * l and α * l. The value of α * l is used as the new lower bound LB = max(LB, α * l).

Step 3: The lower layer is optimized based on the optimization results of the upper layer, and the optimized results f _l ( x * l) and the worst-case scenario x * l are obtained. The upper bound is updated as UB = min(UB, f _l ( x * l)).

Step 4: If UB-LB < ε, where ε is a threshold of convergence, then, the optimal solutions can be obtained, and the iteration is stopped. Otherwise, the variable y ^l+1 and the following constraints are added: α ≥ k T y l + 1 B x + C y l + 1 ≤ D R x u l + 1 * = O y l + 1 S u * l + 1 y l + 1 ≤ V G y l + 1 ≤ H T y l + 1 . (39)

Let l = l + 1, and step 2 is repeated until the algorithm converges.

In order to verify the effectiveness of the proposed method, an actual 35 kV/10 kV DN in China is utilized for analysis. The total loads in this DN are 7.3 MVA, and the total capacity of PVs in this DN is 4 MW. The installed PV capacity is shown in Table 1. The topology of this system is shown in Figure 1, the lines to be planned are shown in Table 2, and the line parameters of the network are shown in Table 3.

The simulated genetic annealing algorithm used in this paper sets the population size n = 40, the initial temperature T ₀ = 100°C, the termination temperature T _end = 1°C, the temperature cooling factor r = 0.8, the maximum number of genetic generations N = 500, the crossover probability p _c = 0.4, and the variation probability p _m = 0.2. This paper uses MATLAB version 2019 to program the algorithm.

In order to illustrate the superiority of the proposed cluster partition index in this paper, the modularity function used by Wang et al. (2021) is selected to compare with the proposed cluster partition index. The weights of the proposed cluster partition index in this paper are taken as ω ¹ = ω ² = ω ³ = ω ⁴ = 0.25. The cluster partition results under the proposed cluster partition index are shown in Figure 2, and the cluster partitioning results under the modularity function are shown in Figure 3.

Figure 2 and Figure 3 show that the size difference among sub-networks attained by the proposed cluster partition index is smaller than that attained by the modularity function, and isolated nodes forming separate clusters (cluster 3 and cluster 6) occur in modularity function. In addition, there are PVs in each cluster attained by the proposed method, while no PVs exist in cluster 1 attained by the modularity function. The differences between the two methods prove that the cluster partition index can result in a more reasonable cluster partition for the subsequent cluster control and operation.

In order to illustrate the superiority of the proposed cluster partitioning algorithm, the traditional genetic algorithm is selected to be compared. The comparison of cluster partition results obtained by the two algorithms is shown in Table 4. Table 4 shows that the optimization results of each cluster partitioning index of the proposed algorithm are greater than those of the traditional genetic algorithm, which indicates that the proposed cluster partitioning algorithm has better optimization performance.

Different combinations of partition index weights are given in Table 5. As shown in Table 5, when increasing the weight of the modularity index, the nodes within the cluster are better connected, but the cluster power complementarity decreases. When increasing the weight of the active and reactive power balance index, the cluster power complementarity characteristics improve, but the nodes within the cluster are significantly less connected. When increasing the weight of the nodal size index, the number of clusters changes. When increasing the weight of the nodal size index, the number of clusters changes, and the cluster partition is more focused on the change in cluster size.

In order to verify the flexibility of the proposed two-layer expansion planning model, three sets of uncertain regulation parameters are selected and compared in terms of network connection, planned capacity of PVs and ESs, and total planning costs. Three sets of uncertain regulation parameters are set as follows:

Case 1:

Г ^PL = Г ^QL = Г ^PV = 0 is set, and in this case, the active load power, reactive load power, and PV output are all equal to the predicted values.

Case 2:

Г ^PL = Г ^QL = 16 and Г ^PV = 4 are set, and in this case, the active and reactive load power of 16 nodes is taken to the minimum value of the prediction interval, while the PV outputs of 4 nodes are taken to the maximum value of the prediction interval.

Case 3:

Г ^PL = Г ^QL = 27 and Г ^PV = 6 are set. In this case, the active and reactive load power of 27 nodes is taken to the minimum value of the prediction interval, while the PV outputs of 4 nodes are taken to the maximum value of the prediction interval, which is the worst scenario.

The comparisons of the three cases in planned PV capacity and ES capacity are shown in Figure 4 and Figure 4, and the planning results and total planning costs are shown in Table 6 and Table 7, respectively.

Figure 4 shows that with the increase in uncertainty regulation parameters, the planned PV capacity decreases. The reason is that the scenario becomes increasingly severe, and to maintain the safe operation of the cluster, the corresponding PVs are reduced.

As shown in Table 6, with the increase in uncertainty regulation parameters, the connection of nodes 21 and 27 is different under the three cases. The reason is that as the operating scenario becomes more severe, and the load power increases, the planning results tend to favor a reduction in the power supply range.

As shown in Figure 5, as the conservatism of the planning increases, the ES installed capacity also increases. The reason is that the active power and reactive power of the load in cases 2 and 3 are smaller than those of case 1, and the PV outputs are larger; hence, more ESs are needed to mitigate the uncertainty of PVs and loads.

In addition, Table 7 shows that as the uncertainty regulation parameter increases, although the annual investment and operation costs of PVs correspondingly decrease, the annual income of PVs decreases, and the annual investment costs of networks, network losses, and the annual investment costs of ESs correspondingly increase. The reason is that the more uncertainty of DNs considered in the planning process, the more conservative the resulting solution becomes, and the corresponding total cost increases.

In order to further illustrate the superiority of the proposed planning method, the centralized planning method is selected to compare with the proposed planning method. The uncertain regulation parameters in this comparison are set as Г ^PL = Г ^QL = 12 and Г ^PV = 3, and in this case, the active and reactive load power of 12 nodes is taken to the minimum value of the prediction interval, while the PV outputs of 3 nodes are taken to the maximum value of the prediction interval. The network connection of the two methods is shown in Figure 6 and Figure 7. The comparison of the calculation time and total planning cost of the two methods is shown in Table 6.

Figure 6 and Figure 7 show that the connection of node 11 and node 27 under the proposed method is different from the centralized planning scheme. The length of the power supply lines for node 11 and node 27 under the proposed method has been reduced by 6.73 km and 6.02 km, respectively; the power supply range is greatly shortened, which can effectively reduce the line investment costs and network loss costs, as well as improve the economic efficiency of the planning scheme.

The calculation method for the difference rate given in Table 8 is to subtract the results of the proposed method from those of the centralized planning method and divide them by the results of the centralized planning method. As shown in Table 8, in the proposed method, the annual PV investment costs and PV operation costs, as well as the annual income of PVs, increased by 5.3% and 14.1%, respectively, compared to the centralized planning method, and the annual investment costs of lines, network losses, and annual investment costs of ESs decreased by 13.5%, 25.7%, and 15.5%, respectively. These different rates indicate that the proposed method can improve the overall economic efficiency of the planning, as well as reduce the costs of planning. In addition, the calculation time of the proposed method is 20.41 s, but that of the centralized planning method is 43.62 s. Owing to the fact that the centralized planning method requires a whole optimization process, the huge variables complex the optimization. However, the proposed method can partition the complex network into smaller sub-networks, and the complex problem is divided into several relatively simple sub-problems that can be solved rapidly and easily by a parallel computing way. Therefore, the computing time can be greatly shortened. The optimization speed of the proposed method is improved by 53.2% compared to centralized planning methods, which indicates that proposed method is more suitable for DN planning with large-scale PVs.