1 Introduction
With the steady growth of the national economy, electricity demand in China has continued to rise, and the need for a low-carbon energy transition has become increasingly urgent. Distributed photovoltaic (DPV) generation has received great attention in this context because of its relatively low investment cost, flexible siting, and obvious environmental benefits. Large numbers of DPV units have been integrated into medium- and low-voltage distribution networks. However, a distribution network is no longer a passive system that only delivers power from the substation to consumers. The inherent randomness and volatility of solar irradiance change the magnitude and direction of power flows in unpredictable ways. When DPV capacity is planned without considering the spatiotemporal characteristics of other resources, the system can suffer from voltage deviation, reverse power flow, increased loss, and even mal-operation of protection.
At the same time, electric vehicles (EVs) have experienced explosive growth in many countries, and the corresponding charging infrastructure has become an important type of load in distribution networks. Electric vehicle charging stations (EVCSs) impose an intensive and time-varying burden on the grid. If charging demand is not coordinated with renewable generation, the coincidence of charging peaks and load peaks could increase the peak-to-valley difference, overload transformers, and reduce power quality. In fact, EV charging load and PV generation exhibit a strong complementary potential. PV production usually reaches its maximum around noon, whereas many EVs remain parked at workplaces and could be charged during this period. Conversely, uncontrolled evening charging coincides with the residual load peak but not with PV production. Therefore, treating DPV and EVCS planning separately often leads to sub-optimal solutions. A joint optimization framework can exploit the temporal complementarity between generation and load, improve renewable energy utilization, defer network upgrades, and lower the total cost of ownership for all stakeholders.
The planning problem is complicated by several modelling and computational challenges. First, both DPV generation and EV charging demand are uncertain. The correlation between them, as well as their seasonal and diurnal patterns, must be preserved, but using one year of hourly data directly in an optimization model makes the problem nearly intractable. Therefore, an effective scenario generation and reduction technique is required. Second, the economic objective function includes investment cost, operating cost, and technical indices such as voltage quality and active power loss. The model is therefore a non-linear, mixed-integer, large-scale multi-objective problem that cannot be reliably solved by conventional mathematical programming methods. Third, heuristic algorithms often suffer from premature convergence and unstable solution quality when applied to this high-dimensional problem. Thus, improvements to the solution algorithm are necessary.
The aim of the present research is to develop a systematic procedure for the joint optimal allocation of DPV plants and EVCSs in a distribution network. The procedure comprises four key elements:
- a data-preprocessing and feature-extraction framework;
- a novel representative scenario extraction method based on an improved density-based spatial clustering algorithm, termed improved DBSCAN;
- a chance-constrained-like multi-objective planning model including investment, operation, and energy-loss terms, together with network constraints;
- a modified grey wolf optimization algorithm, named IGWO, with adaptive position updating and individual re-location, which is used to solve the planning model.
The main contributions of the study are demonstrated by applying the framework to a modified IEEE 37-node feeder and then to a real 62-node industrial park distribution system in north-eastern China. The results show that the proposed method can significantly reduce total cost and network loss, improve terminal voltage profiles, and produce robust planning schemes under different operating conditions.
2 Typical Scenario Extraction Using an Improved DBSCAN Algorithm
2.1 Data Preprocessing and Feature Selection
The data used in this study include hourly DPV generation profiles, conventional load profiles, and EV charging load profiles generated by Monte Carlo simulation. Before conducting cluster analysis, the raw data may contain missing points, outliers, or data recorded in different units. The direct use of low-quality data can seriously affect the quality of scenario extraction. I therefore apply two preliminary steps.
Data cleaning. Outliers, such as negative power values or sudden jumps beyond the physical limit, are detected and replaced by linear interpolation between the two nearest reliable neighbours. If the abnormal value occurs at time \(t\), while the most recent reliable points are \((t_1,y_1)\) and \((t_2,y_2)\), the corrected value \(y\) is
\[
y=y_1+\frac{t-t_1}{t_2-t_1}(y_2-y_1) \qquad\qquad (1)
\]
This interpolation preserves the dynamic trend of the original series and avoids artificial discontinuities.
Normalization. DPV output, base load, and EVCS demand have different physical units and different magnitudes. To avoid the dominance of large numerical values in distance-based clustering, I project every raw value into the range \([0,1]\) using min–max normalization:
\[
x_{\rm norm}=\frac{x-x_{\min}}{x_{\max}-x_{\min}} \qquad\qquad (2)
\]
After preprocessing, I use four statistical features to represent each day in the historical data set. They are the daily average PV output, the daily max-min difference of PV output, the daily average load, and the daily max-min difference of load. The equations are summarized in Table 1. These features encapsulate not only the average operating point but also the volatility of each day. They provide good intra-cluster compactness and inter-cluster separation.
| Feature | Formula | Meaning |
|---|---|---|
| Daily average PV output | \(P_{\mathrm{avg,pv}} = \frac{1}{N}\sum_{t=1}^{N}P_{\mathrm{pv}}(t)\) | Average solar generation level over a day |
| Daily PV range | \(\Delta P_{\mathrm{max,pv}} = \max P_{\mathrm{pv}}(t)-\min P_{\mathrm{pv}}(t)\) | Intra-day fluctuation of PV production |
| Daily average load | \(L_{\mathrm{avg}} = \frac{1}{N}\sum_{t=1}^{N}L(t)\) | Average electricity consumption level |
| Daily load range | \(\Delta L_{\mathrm{max}} = \max L(t)-\min L(t)\) | Peak-valley difference of load |
These four values form a four-dimensional feature vector for each day. The vector is used as an input to the clustering algorithm. Daily records that have similar average values and similar fluctuation levels are mapped to nearby positions in the feature space, while days with unusual variation become isolated points.
2.2 Standard DBSCAN and Its Limitation
DBSCAN is a density-based clustering algorithm that groups points that are closely packed together and identifies points in low-density regions as noise. For a point \(p\), the \(\varepsilon\)-neighbourhood is
\[
N_{\varepsilon}(p)=\{q \in D \mid \operatorname{dist}(p,q) \le \varepsilon\} \qquad\qquad (3)
\]
The point \(p\) is regarded as a core point when the number of points in its neighbourhood satisfies
\[
|N_{\varepsilon}(p)| \ge \mathrm{MinPts} \qquad\qquad (4)
\]
DBSCAN can discover arbitrarily shaped clusters and does not require the number of clusters to be fixed in advance. It is also robust with respect to outliers. However, its final clustering result is sensitive to the two parameters \(\varepsilon\) and \(\mathrm{MinPts}\). When the density of the data varies significantly, as in PV and load patterns, a single global parameter setting cannot capture both dense normal scenarios and sparse but statistically important extreme scenarios. Moreover, DBSCAN does not guarantee stable clustering when the parameter values are not well chosen.
2.3 Improved DBSCAN Combined with K-Means
To overcome this weakness, I combine standard DBSCAN with K-means. The proposed approach proceeds in two stages. In the first stage, an initial density clustering is performed using DBSCAN. The number of clusters and the initial cluster centres are estimated instead of being set randomly. In the second stage, the K-means algorithm is initialized with these centres and then refines the partition. K-means is a fast centroid-based method whose iterative distance minimization is suitable for producing a clean and stable result after DBSCAN has removed the noise. The steps of the improved procedure are as follows:
- Construct the feature matrix from historical PV and load data.
- Apply DBSCAN to the normalized feature matrix to determine a preliminary cluster structure. Noise points are temporarily removed in order to find stable high-density regions.
- Obtain the number of clusters \(K\) and the coordinates of the cluster centres from the DBSCAN result.
- Input these centres as the initial centroids of K-means.
- Iterate K-means by assigning each data point to the nearest centroid and recomputing the centroid until convergence.
- Assign each original sample, including noise points, to its closest final cluster.
I use the Davies–Bouldin index (\(DBI\)) to validate the cluster quality. For two clusters \(i\) and \(j\), the average intra-cluster dispersion is
\[
S_i = \left( \frac{1}{T_i} \sum_{X \in A_i} \Vert X – A_j \Vert^q \right)^{1/q} \qquad\qquad (5)
\]
and the distance between cluster centres is
\[
M_{ij} = \left( \sum_{k=1}^{p} \Vert a_{ki} – a_{kj} \Vert^m \right)^{1/m} \qquad\qquad (6)
\]
The similarity between clusters is defined as
\[
R_{ij} = \frac{S_i + S_j}{M_{ij}} \qquad\qquad (7)
\]
Finally, the \(DBI\) is the average of the maximum similarities:
\[
K_{\mathrm{DBI}} = \frac{1}{K}\sum_{i=1}^{K} \max_{j\neq i} R_{ij} \qquad\qquad (8)
\]
Lower \(DBI\) indicates that the clusters are compact and well separated.
2.4 Numerical Example of Scenario Extraction
I apply the improved DBSCAN method to the PV output and load data of a demonstration park for one complete year. The original data set consists of 365 daily profiles with hourly resolution. After calculating the four features, I test different numbers of clusters. The smallest \(DBI\) value is obtained when \(K=8\), indicating that eight typical days can preserve the main statistical features of the whole year. The clustering result is summarized in Table 2.
| Scenario | PV range (p.u.) | PV mean (p.u.) | Load range (p.u.) | Load mean (p.u.) | Days |
|---|---|---|---|---|---|
| S1 | 0.378 | 0.136 | 0.600 | 0.567 | 17 |
| S2 | 0.822 | 0.257 | 0.643 | 0.593 | 25 |
| S3 | 0.940 | 0.320 | 0.279 | 0.652 | 78 |
| S4 | 0.921 | 0.303 | 0.418 | 0.413 | 37 |
| S5 | 0.457 | 0.187 | 0.294 | 0.307 | 40 |
| S6 | 0.247 | 0.152 | 0.292 | 0.318 | 60 |
| S7 | 0.644 | 0.191 | 0.308 | 0.375 | 49 |
| S8 | 0.692 | 0.219 | 0.285 | 0.272 | 59 |
The corresponding weight \(\eta_i\) of each scenario is obtained by
\[
\eta_i = \frac{n_i}{N} \qquad\qquad (9)
\]
where \(n_i\) is the number of days belonging to scenario \(i\) and \(N=365\). For instance, scenario S6, which represents low PV production and relatively low load, covers 60 days and has a weight of 0.164. Scenario S3 represents high solar irradiance with a relatively heavy but stable load and has a weight of 0.214. This weighted representation is used in all subsequent optimization calculations.
3 Modelling and Mathematical Formulation of the Joint Planning Problem
3.1 PV Output Uncertainty Model
Solar irradiance at a given location can be described by a Beta distribution. If \(S\) is the irradiance and \(S_{\max}\) is its maximum value, the probability density function is
\[
f(S)=\frac{1}{B(\alpha,\beta)}\left(\frac{S}{S_{\max}}\right)^{\alpha-1}
\left(1-\frac{S}{S_{\max}}\right)^{\beta-1} \qquad\qquad (10)
\]
where the Beta function is
\[
B(\alpha,\beta)=\frac{\Gamma(\alpha)\Gamma(\beta)}{\Gamma(\alpha+\beta)} \qquad\qquad (11)
\]
Under the assumption that the PV output is linearly proportional to irradiance, the probability density function of the PV output \(P_{\mathrm{pv}}\) can be transformed into
\[
f(P_{\mathrm{pv}})=\frac{1}{B(\alpha,\beta)}
\left(\frac{P_{\mathrm{pv}}}{P_{\mathrm{pv,max}}}\right)^{\alpha-1}
\left(1-\frac{P_{\mathrm{pv}}}{P_{\mathrm{pv,max}}}\right)^{\beta-1} \qquad\qquad (12)
\]
The parameters \(\alpha\) and \(\beta\) are estimated from historical irradiance records. The Beta distribution is flexible enough to represent different climates by adjusting the two shape parameters. For planning purposes, I first generate many PV output sequences using Monte Carlo sampling and then reduce them to the eight typical scenarios described in Section 2. In load-flow calculation, the DPV installation is modelled as a PQ node with unit power factor, i.e. only active power is injected.
3.2 Load Modelling with Gaussian Mixture Model
The conventional load in the distribution network has multiple peaks and strong stochastic variation. A single Gaussian distribution cannot adequately represent such a bimodal or multimodal probability shape. I therefore construct the probability distribution of the base load at each bus using a Gaussian mixture model:
\[
f_i(P^{\mathrm{base}}_{i,t})=\sum_{k=1}^{K}\pi_k\,
\mathcal{N}\left(P^{\mathrm{base}}_{i,t}\mid \mu_{k,t},\sigma_{k,t}^2\right) \qquad\qquad (13)
\]
with
\[
\mathcal{N}(p \mid \mu,\sigma^2)=\frac{1}{\sqrt{2\pi}\sigma}
\exp\left[-\frac{(p-\mu)^2}{2\sigma^2}\right] \qquad\qquad (14)
\]
The mixture coefficients satisfy \(\sum_{k=1}^{K}\pi_k=1\) and \(0 \le \pi_k \le 1\). I use the expectation maximization algorithm to estimate the means, variances, and weights of the Gaussian components. By adjusting the number of components \(K\), the model can represent bimodal and asymmetric load distributions caused by the difference between working days and weekends.
3.3 Stochastic Charging Model of an EVCS
The charging load of an EVCS is the aggregate demand of many electric vehicles. The charging behavior of each vehicle is described by its arrival time \(t_{\mathrm{arr}}\), initial state of charge \(SOC_{\mathrm{ini}}\), departure time \(t_{\mathrm{dep}}\), and rated charging power. For a single electric vehicle, the instantaneous charging indicator is
\[
S_{\mathrm{EV},k}(t)=
\begin{cases}
1, & t\in[t_{\mathrm{arr},k},\, t_{\mathrm{arr},k}+T_{\mathrm{ch},k}] \\
0, & \text{otherwise}
\end{cases} \qquad\qquad (15)
\]
The required charging duration is
\[
T_{\mathrm{ch},k}=\frac{E_{\mathrm{rated},k}(SOC_{\mathrm{max}}-SOC_{\mathrm{ini},k})}
{P_{\mathrm{rated},k}\,\eta_{\mathrm{ch}}} \qquad\qquad (16)
\]
where \(E_{\mathrm{rated},k}\) is the rated battery capacity, \(P_{\mathrm{rated},k}\) is the rated charging power, and \(\eta_{\mathrm{ch}}\in[0.9,0.95]\) is the charging efficiency. The initial SOC is sampled from a uniform distribution between the allowable lower and upper limits.
For an EVCS \(i\), the total charging power at time \(t\) is the sum of all electric vehicles being served:
\[
P_{\mathrm{EVCS},i}(t)=\sum_{k=1}^{N_{\mathrm{EV}}}
P_{\mathrm{rated},k}\,S_{\mathrm{EV},k}(t) \qquad\qquad (17)
\]
Because a large number of vehicles is considered, Monte Carlo simulation is used to obtain an aggregated charging curve for each EVCS. The resulting charging profile is then added to the base load profile to form the total demand seen by the distribution network. In the optimization model, the EVCS is treated as a new load connected to a candidate bus.
3.4 Power Flow Calculation for a Radial Distribution Network
The distribution network considered in this work is radial. For such a network, the forward-backward sweep method is an efficient and robust power-flow solver. The net injected active and reactive power at bus \(i\), in scenario \(s\) at time \(t\), is written as
\[
P^{\mathrm{net}}_{i,t,s}=P^{\mathrm{pv}}_{i,t,s}-P^{\mathrm{base}}_{i,t,s}-P^{\mathrm{EVCS}}_{i,t,s} \qquad\qquad (18)
\]
\[
Q^{\mathrm{net}}_{i,t,s}=Q^{\mathrm{pv}}_{i,t,s}-Q^{\mathrm{base}}_{i,t,s}-Q^{\mathrm{EVCS}}_{i,t,s} \qquad\qquad (19)
\]
If the DPV inverter operates at unity power factor, the reactive output of the PV system is set to zero. In this formulation, a DPV installation is modelled as a negative load, while an EVCS is modelled as a positive additional load. Therefore, all non-slack buses can be represented as PQ buses with known constant power injections, which simplifies the forward-backward sweep.
The iterative procedure contains a forward sweep and a backward sweep. In the backward sweep, starting from the farthest branch, the current in branch \(i-j\) is calculated from the nodal voltage and the complex power flowing at the receiving end. The voltage at the sending bus is updated as
\[
\dot{V}_i^{(k)}=\dot{V}_j^{(k)}+Z_{ij}\,\frac{\dot{S}_{j}^{\,(k)}}{\dot{V}_{j}^{\,(k)*}} \qquad\qquad (20)
\]
where \(Z_{ij}=R_{ij}+jX_{ij}\) is the branch impedance. In the forward sweep, the branch power loss and the sending power are calculated:
\[
\Delta \dot{S}_{ij}^{(k)}=Z_{ij}\,
\frac{\left(P_j^{\,(k)}\right)^2+\left(Q_j^{\,(k)}\right)^2}
{\left(V_j^{\,(k)}\right)^2} \qquad\qquad (21)
\]
\[
\dot{S}_{ij}^{(k)}=\dot{S}_{j}^{(k)}+\Delta \dot{S}_{ij}^{(k)} \qquad\qquad (22)
\]
The procedure is repeated until the maximum voltage mismatch between two successive iterations is smaller than the convergence tolerance.
3.5 Multi-objective Optimization Model for DPV and EVCS Allocation
The first objective is the minimization of total investment cost. It includes the procurement and installation cost of DPV units, the construction cost of EVCSs, and the cost of additional feeder branches required by the new equipment:
\[
C_{\mathrm{inv}}=C_{\mathrm{PV}}\,P_{\mathrm{PV},i}\,N_{\mathrm{pv}}+
C_{\mathrm{CS}}\,P_{\mathrm{CS},j}\,N_{\mathrm{cs}}+
\sum_{l=1}^{N_{\mathrm{line}}}c_{\mathrm{line}}\,L_l \qquad\qquad (23)
\]
where \(N_{\mathrm{pv}}\) is the number of PV units installed, \(P_{\mathrm{PV},i}\) is the rated capacity of one PV unit, \(C_{\mathrm{PV}}\) is the specific cost per kW, and \(N_{\mathrm{cs}}\) is the number of EVCSs installed. The line cost term depends on the number of new feeder branches, the length of every new branch, and the feeder unit cost.
The second objective is the minimization of the annual operating cost. It contains two parts: yearly operation and maintenance cost of DPV and EVCS equipment, and the net cost of electricity purchased from the upper grid. Since the DPV output reduces the amount of purchased power, the operating cost can be expressed by weighting the eight typical scenarios:
\[
C_{\mathrm{op}}=\sum_{s=1}^{N_s}\eta_s\,
\sum_{t=1}^{24}\left[
c_{\mathrm{buy},t}\,P^{\mathrm{grid}}_{t,s}+
OM_{\mathrm{pv}}P^{\mathrm{pv}}_{t,s}+OM_{\mathrm{cs}}P^{\mathrm{cs}}_{t,s}
\right] \qquad\qquad (24)
\]
where \(P^{\mathrm{grid}}_{t,s}\) is the power drawn from the substation in scenario \(s\), and \(c_{\mathrm{buy},t}\) is the time-of-use electricity price. \(OM_{\mathrm{pv}}\) and \(OM_{\mathrm{cs}}\) are the operation and maintenance costs per kilowatt-hour.
The third objective is the minimization of the annual network energy loss. The active power loss of the network in scenario \(s\) at time \(t\) is calculated by summing over all branches:
\[
P^{\mathrm{loss}}_{t,s}=\sum_{(i,j)\in \Omega}
G_{ij}\left(V_{i,t,s}^2+V_{j,t,s}^2-2V_{i,t,s}V_{j,t,s}\cos\delta_{ij,t,s}\right)
\qquad\qquad (25)
\]
The objective is then
\[
C_{\mathrm{loss}}=365\,\sum_{s=1}^{N_s}\eta_s\,\sum_{t=1}^{24}c_{\mathrm{e}}\,P^{\mathrm{loss}}_{t,s}
\qquad\qquad (26)
\]
where \(c_{\mathrm{e}}\) is the unit cost of energy loss.
The three objectives are converted into a scalar fitness function by assigning the entropy weight coefficients \(\alpha_1\), \(\alpha_2\), \(\alpha_3\). Each normalized sub-criterion is denoted by \(C’_1\), \(C’_2\), and \(C’_3\). The total fitness is
\[
C=\alpha_1\,C’_1+\alpha_2\,C’_2+\alpha_3\,C’_3 \qquad\qquad (27)
\]
In the entropy-weight method, the normalized decision matrix is first constructed. For every objective \(j\), the information entropy is computed as
\[
E_j=-k\sum_{i=1}^{M} \tilde{p}_{ij}\ln(\tilde{p}_{ij}) \qquad\qquad (28)
\]
The deviation degree is \(g_j=1-E_j\), and the weight of the \(j\)-th objective is
\[
\alpha_j=\frac{g_j}{\sum_{j=1}^{3}g_j} \qquad\qquad (29)
\]
The model is subject to several operational constraints:
(1) Power flow constraints: the active and reactive power balance equations must be satisfied at every bus and in every scenario.
(2) Voltage limits: all nodal voltages must remain between the lower and upper bounds
\[
U_{\min}\le U_{i,t,s}\le U_{\max} \qquad\qquad (30)
\]
(3) Branch current limits:
\[
I_{ij,t,s}\le I_{ij,\max} \qquad\qquad (31)
\]
(4) Capacity limits of DPV and EVCS:
\[
P_{i}^{\mathrm{pv}}\le P_{i}^{\mathrm{pv,max}},
\qquad
P_{j}^{\mathrm{cs}}\le P_{j}^{\mathrm{cs,max}} \qquad\qquad (32)
\]
4 Improved Grey Wolf Optimization Algorithm and Its Application
4.1 Standard Grey Wolf Optimizer
The grey wolf optimizer is a swarm-intelligence algorithm inspired by the hierarchical hunting behaviour of grey wolves. In the mathematical model, the best solution is denoted by the alpha wolf \(\alpha\), the second best by the beta wolf \(\beta\), the third best by the delta wolf \(\delta\), and all remaining individuals are regarded as omega wolves \(\omega\). The position update of each individual is based on the following distance and position equations:
\[
\mathbf{D}=|\mathbf{C}\cdot \mathbf{X}_{p}(t)-\mathbf{X}(t)| \qquad\qquad (33)
\]
\[
\mathbf{X}(t+1)=\mathbf{X}_{p}(t)-\mathbf{A}\cdot \mathbf{D} \qquad\qquad (34)
\]
\[
\mathbf{A}=2a\,\mathbf{r}_1-a,\qquad \mathbf{C}=2\mathbf{r}_2 \qquad\qquad (35)
\]
where \(a\) is decreased linearly from 2 to 0 during the iterations. The positions of alpha, beta, and delta wolves are used to guide the whole population. The individual position is updated by
\[
\mathbf{X}(t+1)=\frac{\mathbf{X}_1+\mathbf{X}_2+\mathbf{X}_3}{3} \qquad\qquad (36)
\]
where \(\mathbf{X}_1\), \(\mathbf{X}_2\), and \(\mathbf{X}_3\) are candidate positions calculated from the three leading wolves.
This standard mechanism is simple and effective for many test problems. However, the equal weighting in Eq. (36) does not reflect the different experience levels of the three leaders. Moreover, if alpha, beta, and delta simultaneously fall into a local basin, Eq. (36) misleads the whole population and reduces diversity. Convergence therefore becomes premature when the objective function is strongly non-linear and high-dimensional.
4.2 Improved IGWO Algorithm
I introduce two improvements to overcome the above drawbacks. First, an adaptive weight is introduced in the position update equation. The weights are assigned according to the relative fitness level of the three leading wolves:
\[
\omega_1=\frac{X_{\alpha}}{X_{\alpha}+X_{\beta}+X_{\delta}},
\quad
\omega_2=\frac{X_{\beta}}{X_{\alpha}+X_{\beta}+X_{\delta}},
\quad
\omega_3=\frac{X_{\delta}}{X_{\alpha}+X_{\beta}+X_{\delta}} \qquad\qquad (37)
\]
The new position update is
\[
\mathbf{X}(t+1)=\omega_1\mathbf{X}_1+\omega_2\mathbf{X}_2+\omega_3\mathbf{X}_3 \qquad\qquad (38)
\]
This modification gives the alpha wolf more influence at the beginning of the search and lets the beta and delta wolves provide a corrective effect when the alpha position is not sufficiently accurate.
Second, an individual relocation mechanism is used to improve the performance of low-fitness individuals. Let \(d_1\in[0.5,1]\), \(d_2\in[0.25,0.5]\), and \(d_3=1-d_1-d_2\). The relocation coordinate is
\[
\mathbf{X}_{\mathrm{reloc}}=d_1\mathbf{X}_{\alpha}+d_2\mathbf{X}_{\beta}+d_3\mathbf{X}_{\delta} \qquad\qquad (39)
\]
Because the alpha wolf has the highest leadership level, the relocation point is placed closer to the alpha position. This operation is applied only in the later stage of the iteration so that the population diversity in the early stage is not destroyed. Low-fitness individuals are moved into the spatial region around the elite group and can thereby contribute to a more refined local search.
4.3 Performance Verification on Benchmark Functions
I verify the IGWO algorithm by using four standard benchmark functions. Two are unimodal functions, namely Sphere and Schwefel or Rosenbrock, and two are multimodal functions, namely Griewank and Rastrigin. The expressions and the global minima are provided in Table 3.
| Function | Expression | Minimum |
|---|---|---|
| Sphere | \(f_1(x)=\sum_{i=1}^{D}x_i^2\) | 0 |
| Rosenbrock | \(f_2(x)=\sum_{i=1}^{D-1}\left[100(x_{i+1}-x_i^2)^2+(x_i-1)^2\right]\) | 0 |
| Griewank | \(f_3(x)=\frac{1}{4000}\sum_{i=1}^{D}x_i^2-\prod_{i=1}^{D}\cos{\frac{x_i}{\sqrt{i}}}+1\) | 0 |
| Rastrigin | \(f_4(x)=\sum_{i=1}^{D}\left[x_i^2-10\cos(2\pi x_i)+10\right]\) | 0 |
The population size of each comparative algorithm is set to 93, the maximum number of iterations is 500, and each test runs 26 times. The mean value (MV), standard deviation (SD), and optimal value (OP) obtained by PSO, GA, GWO, and IGWO are listed in Table 4. The results show that IGWO achieves the smallest optimal value and the smallest standard deviation on almost all tested functions. The improvement of the local search ability reduces the possibility of premature convergence and gives the algorithm a greater chance to approach the global optimum. The adaptive leader weights in Eq. (38) allow a smooth transition from exploration to exploitation, while the relocation mechanism in Eq. (39) reinforces the convergence in the later stage.
| Function | Algorithm | MV | SD | OP |
|---|---|---|---|---|
| Sphere | PSO | 1.778 | 0.381 | 0.059 |
| GA | 0.138 | 0.109 | 0.051 | |
| GWO | 0.038 | 0.038 | 0.047 | |
| IGWO | 0.047 | 0.027 | 0.017 | |
| Rosenbrock | PSO | 26.041 | 22.060 | 25.711 |
| GA | 30.180 | 11.388 | 24.407 | |
| GWO | 27.273 | 3.696 | 17.464 | |
| IGWO | 24.403 | 0.812 | 14.300 | |
| Griewank | PSO | 77.694 | 16.958 | 24.368 |
| GA | 124.363 | 67.154 | 21.913 | |
| GWO | 44.173 | 23.151 | 22.131 | |
| IGWO | 31.010 | 9.009 | 19.040 | |
| Rastrigin | PSO | 0.065 | 0.055 | 0.039 |
| GA | 0.046 | 0.361 | 0.041 | |
| GWO | 0.038 | 0.048 | 0.040 | |
| IGWO | 0.003 | 0.016 | 0.025 |
4.4 Entire Solution Procedure of DPV and EVCS Joint Planning
The joint planning problem contains both discrete sizing decisions and continuous placement decisions. The IGWO algorithm is employed as the global optimizer. The detailed calculation steps are:
- Load the eight typical scenarios and their weights obtained from improved DBSCAN.
- Initialize the population of wolves. Each individual encodes the candidate capacities of DPV units and charging stations at the preselected buses.
- Check the capacity limits and run the forward-backward power flow for every scenario and every time interval.
- Calculate the investment cost, operation cost, and network loss for each individual. Normalize all three sub-criteria and calculate the entropy weights.
- Compute the scalar fitness function of each wolf.
- Designate the first three wolves with the best fitness values as alpha, beta, and delta. Update all the other wolves according to Eqs. (7) and (8).
- Apply the individual relocation process when the iteration counter exceeds half of the maximum iteration number.
- Repeat steps 5–7 until the stopping criterion is satisfied.
- Output the position of the alpha wolf as the optimal allocation scheme.
5 Case Study on an IEEE 37-Node Distribution Feeder
Before applying the method to a practical urban park, I first test the joint planning framework on a modified IEEE 37-node distribution feeder. Three extra load buses are incorporated so that the feeder becomes more representative of a real distribution network with different consumer types and dispersed load. The candidate bus data are shown in Table 5. Each bus can receive a DPV unit with a maximum capacity of \(250~\text{kW}\), while each EVCS has a planned capacity range of \([0, 150~\text{kW}]\).
| Start Bus | End Bus | Length (km) |
|---|---|---|
| 38 | 10 | 0.200 |
| 39 | 16 | 0.200 |
| 40 | 35 | 0.200 |
The optimal DPV and EVCS configuration results obtained by GWO, PSO, and IGWO are shown in Table 6. It can be seen that IGWO tends to distribute the DPV capacity among more buses in a balanced way, whereas PSO concentrates a large capacity on one individual bus, and GWO has an intermediate result. A concentrated DPV allocation may be attractive in terms of construction convenience, but it often causes severe voltage rise at the point of connection and larger reverse flows.
| Bus | Device | GWO (kW) | PSO (kW) | IGWO (kW) |
|---|---|---|---|---|
| 18 | DPV | 60.93 | 47.53 | 75.24 |
| 19 | DPV | 101.15 | 107.55 | 62.26 |
| 21 | DPV | 72.54 | 109.23 | 60.85 |
| 28 | DPV | 45.05 | 46.76 | 26.33 |
| 29 | DPV | 150.04 | 123.39 | 79.77 |
| 32 | DPV | 45.53 | 62.04 | 25.02 |
| 33 | DPV | 160.04 | 240.21 | 85.08 |
| 15 | EVCS | 98.45 | 112.83 | 82.14 |
| 25 | EVCS | 75.12 | 63.50 | 48.72 |
| 34 | EVCS | 88.21 | 105.16 | 95.39 |
As shown in Table 7, IGWO obtains a normalized fitness value of only 0.03416, which is roughly 65% smaller than that of GWO and 54% smaller than that of PSO. This substantial improvement arises because the IGWO-based scheme uses a more effective mixture of DPV locations and charging stations, thus reducing both construction cost and energy loss.
| Algorithm | \(C’_1\) | \(C’_2\) | \(C’_3\) | Fitness |
|---|---|---|---|---|
| GWO | 0.6345 | 0.6226 | 0.9145 | 0.09838 |
| PSO | 0.7334 | 0.7236 | 0.9638 | 0.07427 |
| IGWO | 0.4127 | 0.4088 | 0.7749 | 0.03416 |
I also compare the system voltage profiles at the peak-load hour when no DPV or EVCS is connected, when the PSO-based allocation is used, and when the IGWO-based allocation is used. In the original network, some remote buses have voltage magnitudes below \(0.95\;\mathrm{p.u.}\). The PSO-based plan improves the situation but still cannot provide adequate support at the feeder tail because too much DPV is concentrated at a single bus. The IGWO-based plan keeps all voltages comfortably inside the range \(0.95–1.05\;\mathrm{p.u.}\), and the voltage profile is much smoother along the feeder. Therefore, the joint planning framework yields a better trade-off between economic performance and voltage security.
6 Application to a Practical Park Distribution Network
6.1 System Description
For the final application study, I consider a real industrial and commercial park distribution network. The park is supplied by a 10-kV substation and contains a single-line radial distribution system with 62 buses, 61 branches, and a peak demand of \(115\;\text{MW}\). There are seven groups of distributed PV systems, each with a rated capacity of \(35\;\text{MW}\). The feeder is divided into three areas: an industrial area, a residential area, and a commercial area. EVCSs can be built on several candidate nodes that are located near parking lots and workplace destinations.
The specific investment costs of DPV and EVCS are \(7980~\text{yuan/kW}\) and \(3985~\text{yuan/kW}\), respectively. The annual operation and maintenance cost is set to \(1\%\) of the investment cost. The candidate capacities of the DPV and EVCS units are bounded by the feeder thermal limits. The typical charging curve of an EVCS in the park is characterized by two peaks: a modest noon peak, which matches noon-time PV production, and a much higher evening peak. This temporal pattern is suitable for testing the complementarity between distributed solar generation and electric vehicle charging demand.
6.2 Weighting Factors and Algorithm Comparison
In Section 3.5, the entropy weight method is used to determine the objective weights. For this practical case, the calculated entropy weights are \(\alpha_1=0.411\), \(\alpha_2=0.294\), and \(\alpha_3=0.295\). To simplify the engineering calculation and ensure the repeatability of the method, these values are rounded to \(\alpha_1=0.4\), \(\alpha_2=0.3\), and \(\alpha_3=0.3\). The slightly rounded weights remain close to the original entropy solution and give appropriate importance to the investment cost, operating cost, and technical loss index.
Three algorithms, GWO, PSO, and IGWO, are then independently used to solve the planning model. The fitness convergence histories reveal different characteristics. GWO converges after about 63 iterations, PSO after about 112 iterations, while IGWO needs 177 iterations. Although IGWO requires more iterations, it reaches a significantly lower final fitness value, equal to 0.1091, compared with 0.1857 for GWO and 0.1954 for PSO. This indicates that GWO and PSO converge quickly to a local optimum, while IGWO maintains a more effective global search process and finally obtains a better solution.
| Algorithm | \(C’_1\) | \(C’_2\) | \(C’_3\) | Fitness \(C\) |
|---|---|---|---|---|
| GWO | 0.8767 | -0.2577 | 0.9563 | 0.1857 |
| PSO | 0.9305 | -0.2406 | 0.9864 | 0.1954 |
| IGWO | 0.8640 | -0.2256 | 0.9060 | 0.1091 |
The corresponding economic indices are reported in Table 9. IGWO obtains the lowest investment cost of \(939.44\times10^4\) yuan, the lowest operation cost of \(90.91\times10^4\) yuan, and the smallest annual network loss equal to \(3636.27\) yuan. These results confirm that the proposed algorithm improves all the individual objectives simultaneously rather than merely improving an aggregated fitness value.
| Algorithm | Investment cost (\(10^4\) yuan) | Operation cost (\(10^4\) yuan) | Annual energy loss (yuan) |
|---|---|---|---|
| GWO | 942.04 | 91.77 | 4981.56 |
| PSO | 1023.78 | 92.93 | 4659.44 |
| IGWO | 939.44 | 90.91 | 3636.27 |
The allocation scheme produced by IGWO is summarized in Table 10. Three EVCSs are installed at nodes 22, 25, and 53 with charging capacities of around 130.6 kW, 85.3 kW, and 63.5 kW, respectively. Six DPV units are assigned to nodes 28, 31, 33, 46, 47, and 54. The capacities vary from roughly 91.6 kW to 239.5 kW. The less loaded nodes at the middle part of the feeder receive larger PV capacities, while the DPV units close to the feeder terminal have moderate capacities in order not to cause unacceptable voltage rise.
| Bus | Device | IGWO (kW) |
|---|---|---|
| 22 | EVCS | 130.57 |
| 25 | EVCS | 85.35 |
| 53 | EVCS | 63.46 |
| 28 | DPV | 148.56 |
| 31 | DPV | 135.56 |
| 33 | DPV | 213.46 |
| 46 | DPV | 239.46 |
| 47 | DPV | 195.06 |
| 54 | DPV | 91.65 |
6.3 Effect of Weighting Scheme on Planning Results
To investigate how the preference of the decision-maker affects the final solution, I compare three weighting scenarios:
- Scheme 1: balanced preference obtained from the entropy method, \(\alpha_1=0.4\), \(\alpha_2=0.3\), \(\alpha_3=0.3\);
- Scheme 2: equal weights, \(\alpha_1=0.333\), \(\alpha_2=0.333\), \(\alpha_3=0.333\);
- Scheme 3: operational preference, \(\alpha_1=0.2\), \(\alpha_2=0.4\), \(\alpha_3=0.4\).
The resulting configurations are presented in Table 11. Scheme 3, which assigns a lower importance to investment cost, tends to reduce the installed DPV capacity in order to decrease the network loss and the operation cost. Scheme 2 lies in between. Scheme 1, based on entropy weights, gives the highest initial investment but the lowest operation cost and the lowest loss among all schemes, because more DPV capacity is placed in optimal locations to reduce purchased energy and active power loss.
| Scheme | Investment (\(10^4\) yuan) | Operation cost (\(10^4\) yuan) | Network loss (yuan) |
|---|---|---|---|
| Scheme 1 | 939.44 | 90.91 | 3636.27 |
| Scheme 2 | 925.13 | 95.24 | 3580.15 |
| Scheme 3 | 910.76 | 98.76 | 3522.88 |
An important conclusion is that the direct comparison of investment costs alone can be misleading. The entropy-weighted scheme has a higher initial investment than Scheme 3 by about \(28.7\times10^4\) yuan, but its operation cost is \(7.85\times10^4\) yuan lower every year and its network loss is \(113.4\) yuan lower. Over a typical project lifetime of 20–25 years, the annual saving in operation cost will exceed the initial additional investment. Consequently, the entropy-based weighting strategy produces the best long-term economic performance, while still maintaining a safe and stable network operation.
6.4 Robustness and Technical Performance in the Eight Typical Scenarios
A planning solution obtained under one representative scenario may not perform well in other operating conditions. I therefore apply the IGWO-based optimal plan to each of the eight typical scenarios and evaluate its nodal voltage profiles, daily operation cost, and daily network loss. This scenario-level analysis provides a rigorous robustness check.
6.4.1 Nodal Voltage Profiles
The voltage profiles at the daily peak-load hour are compared among the eight scenarios. In scenario S1, which has little PV output and a high evening load, voltage decreases continuously from the substation to the end of the feeder. Nevertheless, the proposed placement of DPV and EVCSs keeps most voltages within the acceptable range, and only the most remote feeder nodes come close to the lower voltage limit. In scenario S3, where PV generation and load are both relatively high, the DPV supply coincides well with daytime demand. As a result, the voltage profile is smooth and nearly flat along the entire feeder. In scenarios S4 and S8, where PV output is high but the load is low, reverse power flow raises the voltage at the feeder end. The proposed plan nevertheless prevents the voltage from exceeding the upper limit because the DPV capacities are not over-concentrated in the tail part of the feeder. In conclusion, the final configuration is sufficiently robust across all the eight operating scenarios and does not require frequent reactive-power adjustment or transformer tap changes under normal operating conditions.
6.4.2 Operating Cost and Energy-Loss Analysis
Table 12 summarizes the daily average operation cost and daily average loss rate for all eight typical scenarios when the IGWO-based plan is implemented. The results reveal substantial variation. Scenario S1 has the highest daily operation cost because the system heavily depends on purchased electricity at the peak-load period. By contrast, scenario S4 achieves a daily net benefit because PV generation can be partly sold back to the upper grid while the base load remains at a moderate level. In addition, scenario S3 has the lowest daily loss rate, namely 1.57%, because the local DPV production matches the local load throughout the daytime. Scenario S6, although having the smallest load, has a higher loss rate than S3 because all the load is supplied over long distances from the substation. This comparison illustrates that the joint allocation of DPV and EVCS is not simply a matter of maximizing installed PV capacity; it requires an excellent match between generation and demand in both time and space.
| Scenario | Daily operation cost (\(10^4\) yuan) | Loss rate (%) |
|---|---|---|
| S1 | 1.85 | 3.10 |
| S2 | 1.35 | 2.60 |
| S3 | 0.63 | 1.57 |
| S4 | -0.21 | 2.42 |
| S5 | 0.95 | 2.81 |
| S6 | 1.05 | 2.86 |
| S7 | 0.86 | 2.48 |
| S8 | 0.78 | 2.74 |
In a further comparison with the GWO-based and PSO-based planning schemes, the IGWO-based scheme yields the smallest network loss in every month of the year. It consistently reduces the monthly loss, especially in the high-insolation months when PV generation has the greatest potential to reduce line flows. This confirms that the solution found by the improved optimization algorithm is stable and superior from both an economic and an operational perspective.
7 Conclusion
In this research, I developed a comprehensive method for the coordinated optimal allocation of distributed photovoltaic generation and electric vehicle charging stations in distribution networks. The proposed method consists of three tightly coupled components: representative scenario extraction, a multi-objective planning model, and an improved grey wolf optimization algorithm.
First, an improved DBSCAN clustering algorithm was proposed by combining density clustering with K-means iteration. The method reduces an entire year of historical data to eight typical scenarios with well-defined weights. Four statistical features, namely the daily average and daily fluctuation range of PV output and load, are used as the clustering feature set. The improved algorithm overcomes the parameter-sensitivity weakness of standard DBSCAN and produces better cluster separability, as indicated by a lower Davies–Bouldin index.
Second, a stochastic model of the PV source, the base load, and the EVCS load was constructed. The PV output is represented by a Beta distribution. The base load is represented by a Gaussian mixture model to account for its multimodal nature, while the electric vehicle charging load is generated by Monte Carlo simulation based on realistic travel behavior. A joint planning model was constructed that minimizes investment cost, annual operating cost, and network loss while ensuring power-flow feasibility, voltage limits, branch-current limits, and capacity restrictions.
Third, the standard grey wolf optimizer was improved by introducing adaptive leader weights and an individual relocation strategy. The resulting IGWO algorithm performs substantially better than GWO and PSO on standard benchmark functions and on the modified IEEE 37-node feeder. When applied to a 62-bus practical park distribution network, IGWO obtains investment, operation, and loss indices of \(939.44\times10^4\) yuan, \(90.91\times10^4\) yuan, and \(3636.27\) yuan, respectively, both of which are the best among the three tested algorithms. The final configuration also maintains acceptable voltage profiles under all eight typical scenarios and gives the lowest network loss in every month of the year.
In conclusion, the joint planning framework developed in this work can provide useful technical references for the planning and operation of distribution networks containing high shares of distributed photovoltaic resources and electric vehicle charging loads. In future research, I intend to expand the framework by considering vehicle-to-grid services, battery storage, and more complex multi-stage planning decisions. It is also worthwhile to validate the proposed method in several other distribution networks with different load characteristics so that the adaptability and limitations of the algorithm can be more comprehensively assessed.