The actual output power of the wind turbine is related to the wind speed in the current time period (Qu et al., 2017). The output power P _WT of the fan is related to the wind speed v as follows: P W T = 0 0 ≤ v ≤ v c i or v c o ≤ v P r W T v − v c i v r − v c i v c i < v ≤ v r P r W T v r < v ≤ v c o , (1) where P _rWT is the rated power of the wind turbine (WT). v _ci , v _r , and v _co are the cut-in, rated, and cut-out wind speeds, respectively, at which the wind turbine operates.

The actual output power of PV is related to the solar intensity (Karaki et al., 1999). The relationship between output power P _PV and solar intensity L is shown below: P P V = P r P V L L r L ≤ L r P r P V L > L r , (2) where P _rPV is the rated power of PV. L is the actual value of solar intensity. L _r is the rated value of solar intensity.

In this paper, the load is set to obey the normal distribution (Ding et al., 2021), the active load is denoted by P _L , and P _L obeys a normal distribution. Then, the probability density function of the active load f P L is expressed as shown in Equations 3, 4 below: f P L = 1 2 π σ L exp − P L − μ L 2 2 σ L 2 , (3) Q L = P L ⁡ tan ⁡ φ , (4) where Q _L is the reactive load. φ is the power factor angle of the load. μ L and σ L are the expectation and standard deviation of the active load, respectively.

The charging power and discharging power of the energy storage device are determined by the difference between the output power of wind power and PV and the power of the load (Bei et al., 2022). The difference is denoted as Δ P E S S t . When Δ P E S S t is greater than 0, the energy storage device is charged. When Δ P E S S t is less than 0, the energy storage device is discharged. Its charging and discharging model expression is expressed as shown in Equation 5 below: Δ P E S S t = P W T t + P P V t − P L t , (5) where P W T t , P P V t , and P L t are the turbine, PV output, and load power at time t, respectively.

Since wind, PV, and load have strong randomness and volatility, scenario analysis is often used to reduce the incidence of similar scenarios in order to reduce the computational complexity. Firstly, according to the historical data of wind speed, solar intensity and load in the planning area, in order to determine the probability model parameters of the three. Secondly, LHS (Osawa and Katsura, 2018) is used to obtain M corresponding wind speed, solar intensity, and load data. Finally, Equations 1, 2 are used to transform the generated M wind speeds and solar intensities into M wind and PV output data and combined with M load rates to form the original wind–PV–load output base scenario.

The matrices of the wind–solar–load output generated using LHS are all 500 × 24 matrices. In order to improve the accuracy of planning and reduce the repetitiveness caused by duplicate scenarios, those with high similarity need to be removed. We thus propose the AP-DTW-K-medoids scenario reduction method. The AP-DTW-K-medoids method proceeds in three steps:

1) Initial clustering: affinity propagation (AP) generates k _AP candidate centers (k _AP > k) to avoid K-means’ sensitivity to initial values (Ji et al., 2021).

The AP algorithm is used in the first layer to obtain k _AP (k _AP > k) candidate clustering centers for the initial scenarios. The clustering center at this point is a matrix of k A P * 24 . Each clustering center is a time-series curve. The similarity of each curve is calculated at the second layer by the DTW algorithm to remove similar points. Finally, the obtained clustering center matrix is used as the initial clustering center of K-medoids. The specific steps are as follows:

1) Input the initial scenario matrix X = x 1 , x 2 , ⋯ , x n . Given the number of clusters k, select the AP similarity reference type. The threshold parameter is θ ∈ 0 , 1 .

The clustering evaluation index DBI is used to evaluate the final scenario. The smaller the DBI, the better it represents its clustering effect (Yu et al., 2024). It is expressed as follows: V D B I = 1 k ∑ i = 1 k max j ≠ i M i * ¯ + M j * ¯ x i p − x j p , (12) where M i * ¯ and M j * ¯ are the average distances from the cluster data to the clustering center of the i-th and j-th nest, respectively. x i p and x j p are the clustering centers of the i-th and j-th nest, respectively.

The lower layer model is mainly composed of the operation cost and voltage deviation, and the operation cost includes wind and solar abandonment cost, network loss cost, power purchase cost, and operation and maintenance (O&M) cost.

The upper level model objective function represents the total cost, which consists of the operating cost and investment cost, and the expression is shown in Equations 19, 20 below: C T O T A L = C I N V + C O P E , (19) C I N V = r 1 + r n 1 + r n − 1 ∑ j ∈ Ω W T C W T I N V P j W T + ∑ j ∈ Ω P V C P V I N V P j P V + ∑ j ∈ Ω E S S C E S S I N V P j E S S , (20) where Ω W T , Ω P V , and Ω E S S are the set of nodes of the wind turbine, PV, and energy storage to be installed, respectively. C W T I N V , C P V I N V , and C E S S I N V are the investment costs per unit of the wind turbine, PV, and energy storage capacity, respectively. P j W T , P j P V , and P j E S S are the j-th node of the installed capacity of the wind turbine, PV, and energy storage, respectively. n is the service life of the resource. r is the discount rate. C T O T A L is the total cost. C I N V is the investment cost.

The PSO algorithm utilizes a stochastic method to determine the maximum value of a specific region of a multidimensional function. It can be described as follows: there are m particles in D-dimensional space, the position of the i-th particle is a vector xi, the velocity of the i-th particle is a vector V _i , the optimized position searched by the i-th particle is pi, and the optimized position searched by the whole swarm of particles is p _gbest ; the iterative formulas for the particle velocities and positions are expressed as shown in Equations 26, 27 below: V i , d k + 1 = ω V i , d k + c 1 r 1 p i , d − x i , d k + c 2 r 2 p g b e s t − x i , d k , (26) x i , d k + 1 = x i , d k + V i , d k + 1 , (27) where i = 1 , 2 , 3 , ⋯ , m and d = 1 , 2 , 3 , ⋯ , D . ω is the inertia parameter. c ₁ and c ₂ are the learning factors. r ₁ and r ₂ are random numbers.

The optimization performance of the traditional PSO algorithm is highly dependent on the setting of inertia parameters and learning factors. In the iterative process, when the initial particles update their spatial positions and velocities, the differences between different generations of particles are often neglected. To address this problem, this paper imposes constraints on the above parameters and adopts a linearly decreasing parameter setting method to enhance the performance of the traditional PSO algorithm. The specific approach is as follows: first, larger inertia parameters and learning factors are set to start the iteration, then these parameters are adjusted in real time in each iteration, and finally, the iteration is ended with a smaller value. This method enhances the global optimization seeking ability of the initial particle swarm and helps the particle swarm to get rid of the local optimum. The specific update rules of the parameters are expressed as shown in Equation 28 below: ω = ω max − ω max − ω min t max t c 1 = c 2 = c max − c max − c min t max t , (28) where the values of w _max and w _min are 0.9 and 0.7, respectively, and the values of c _max and c _min are 1.5 and 1, respectively.

In the solution process of the PSO algorithm, the historical optimized positions of search particles and global optimized particles are constantly updated to guide other particles to move toward the optimized positions so as to realize the convergence of the algorithm. However, the rapid convergence of the search particles may produce a large number of close-range invalid solutions, and it is easy to make the algorithm fall into the local optimum. We thus introduce an adaptive operator to set the optimization strategy of the particles. When a particle satisfies the position update condition of Equation 30, it can perform position transformation according to Equation 31 to seek a more optimized solution, as shown in Equation 29 below: e = L i , j × r a n d − 0.5 , 0.5 , L i , g b e s t ⩾ e − k / k max L i , j × r a n d − 0.2 , 0.2 , L i , g b e s t < e − k / k max , (29) L i , j L i , g b e a t < Q Q = Q 0 1 − k k max , (30) x i k + 1 = k k max x i k + e 1 − k k max x i k , (31) where e is the adaptive operator. rand (•) is the random number function. e − k / k max is the adaptive threshold. L _i,j is the Euclidean distance between particle i and the nearest particle j in the D-dimensional space formed by the objective function. L _i,gbest is the Euclidean distance between particle i and the optimized particle of the population. Q is the optimization decision threshold.

In the early stage of the algorithm, the initial particle is far away from the optimized particle, resulting in a larger L _i,gbest value, and the corresponding adaptive operator e also increases. At this time, the particle position update is mainly dominated by the adaptive operator term, which enhances the diversity of particles. In the later stage of the algorithm, the particles gradually converge to the neighborhood of the global optimized solution, the influence of the adaptive operator term is weakened, and the particles rely more on their own positions to carry out a fine search, which ensures that the algorithm has the accuracy of a small range of solutions.

The multi-objective PSO algorithm solution results in a Pareto optimized solution set. To select the final scheduling plan, this paper uses the fuzzy decision-making method to fuzzy the two objective functions of operating cost and environmental cost. The processing function is expressed as shown in Equation 32 below: T a r i = f i − f min f max − f min , (32) where f _i is the value of the single objective function of the i-th Pareto optimized solution. Tar _i is the value of the objective function after fuzzy processing. f _max and f _min are the maximum and minimum values of the corresponding objective function, respectively.

In order to verify the correctness and effectiveness of the improved PSO algorithm, we set up a specific solution process and optimization calculation method. The process is as follows: the particle swarm is initialized, the speed threshold and the overrun processing mechanism are set according to the equipment model, and the initial calculation results and the historical and global optimized position are recorded. The inertia parameters and the learning factor are updated in each iteration, the adaptive optimization conditions of the particles are judged and calculated, and the Pareto-optimized solution set is recorded. The process is stopped after iterating to a set number of times, and the collection of historical optimized values is output. The specific flow is shown in Figure 1.

The IEEE-33-bus system, as shown in Figure 2, is used for simulation. The base voltage is 12.66 kV, and the specific system parameters are shown in the literature (Baran and Wu, 1989). The DG types are WT and PV. Their specifications are detailed in Table 1. While the IEEE-33-bus system provides a controlled test environment, its radial topology and homogeneous load profile may under-represent complexities of urban distribution networks with multiple feeders and diverse load types.

In comparison to non-dominated sorting genetic algorithm II (NSGA-II), the enhanced MOPSO algorithm demonstrated a 12.4% faster convergence rate and a 15.7% improvement in hypervolume metric, indicating superior exploration–exploitation balance. Against multi-objective evolutionary algorithm based on decomposition (MOEA/D), the proposed algorithm achieved a 9.8% reduction in voltage deviation while maintaining comparable computational efficiency. These results underscore the enhanced MOPSO’s ability to generate high-quality Pareto fronts in complex multi-objective scenarios.

To elucidate relationships between our three objectives, we conducted a comprehensive Pareto front analysis. The “knee” point at 1.05 × baseline cost yields 22% voltage deviation reduction. The Pareto front reveals critical insights: 1) solutions below the knee region (1.05 × baseline cost) show diminishing returns; 2) above the knee, each 1% voltage improvement requires 2.3% additional investment; 3) operating costs increase linearly with voltage deviation reduction.

The normalized Pareto front reveals three operational regions with distinct cost–voltage trade-offs, including a “knee” region at 1.05 × baseline cost that balances 22% voltage improvement with 5% cost increase, beyond which marginal gains diminish rapidly.

The 25.5% improvement in DBI directly translates to enhanced planning decisions through the following: 1) better extreme event representation: it preserves 37% more ramp events (≥5% load change/hour) compared to K-medoids. 2) Improved temporal correlation: it maintains 0.89 Pearson correlation with the original 500-scenario set vs 0.76 for K-means. 3) Reduced over-conservatism: it lowers DG capacity overestimation by 9% through more accurate scenario weighting. These improvements stem from DTW’s ability to capture temporal dynamics missed by Euclidean-based methods (K-means/K-medoids).

Compared to alternative scenario reduction techniques, as shown in Table 6, the proposed method shows the following: 1) A 23% better DBI than the nearest competitor (Fast-SS). 2) A 35% higher ramp event preservation than K-medoids. 3) Maintains 95% computational efficiency of baseline methods. These advantages lead to 7.87% lower voltage deviation and 1.41% cost reduction in planning outcomes.