The evolving energy paradigm, driven by the imperative to mitigate climate change and enhance energy security, has ushered in an unprecedented integration of variable renewable energy sources like wind and photovoltaic (PV) generation into modern power grids. Concurrently, the global push towards electrification of transportation is leading to a rapid proliferation of electric vehicle cars. While this transition is crucial for a sustainable future, the inherent variability and uncertainty of these new resources present significant challenges to power system operation, affecting stability, security, and power quality. Conventional deterministic load flow analysis is ill-equipped to quantify these risks, necessitating advanced probabilistic methods. This paper presents a comprehensive framework for dynamic probabilistic load flow (DPLF) analysis that explicitly models the spatial and temporal uncertainties introduced by correlated wind-PV generation and geographically distributed electric vehicle car charging loads.
The core challenge lies in accurately characterizing the input random variables. Wind and solar power outputs are intrinsically intermittent and weather-dependent. Furthermore, geographically proximate wind farms or PV plants exhibit spatial correlation—their outputs do not vary independently but often rise and fall together due to shared weather patterns. Ignoring this correlation can lead to an underestimation of collective power fluctuations and associated grid impacts. Similarly, the charging demand from a fleet of electric vehicle cars is not a simple aggregate; it is a complex function of user behavior, traffic network dynamics, travel patterns, ambient conditions, and the availability of charging infrastructure. A realistic model must capture this spatiotemporal distribution.
To address these challenges, our methodology is built upon three foundational pillars: 1) a hybrid parametric-nonparametric probability model for wind and PV generation, coupled with a Copula function to capture their interdependence; 2) a high-fidelity model for electric vehicle car charging loads that incorporates road network constraints, detailed travel and charging characteristics, ambient temperature effects, and queueing theory for charging station delays; and 3) an efficient DPLF solution algorithm combining the Nataf transformation with Singular Value Decomposition (SVD) for correlation handling, the cumulant method for probabilistic convolution, and the Gram-Charlier series expansion for reconstructing output probability distributions.
1. Probabilistic Modeling of Input Uncertainties
1.1 Correlated Wind and PV Generation Model
Accurate representation of renewable generation uncertainty is paramount. We employ a hybrid approach. For parametric fitting, wind speed is typically modeled by a Weibull distribution, and solar irradiance by a Beta distribution, which are then transformed to power outputs using the respective generator’s power curve. To enhance accuracy, especially when empirical data deviates from standard parametric forms, we augment this with a nonparametric kernel density estimation (KDE). The probability density function (PDF) for a generation output $P$ is estimated as:
$$f(P) = \frac{1}{nh}\sum_{i=1}^{n} K\left(\frac{P – P_i}{h}\right)$$
where $n$ is the number of historical samples, $h$ is the bandwidth, $K(\cdot)$ is the kernel function (e.g., Gaussian), and $P_i$ is the $i$-th historical power output sample. This combined method provides a more robust fit to real-world data.
The spatial correlation between different wind farms or PV plants is modeled using Copula theory. A Copula function $C$ links marginal distribution functions to form a joint distribution. For two generation sources with outputs $x$ and $y$ and marginal Cumulative Distribution Functions (CDFs) $F_X(x)$ and $F_Y(y)$, their joint CDF is given by:
$$H(x, y) = C[F_X(x), F_Y(y); \theta]$$
where $\theta$ represents the Copula parameters capturing the dependence structure. Different Copula families (e.g., Gaussian, t, Gumbel) can model various types of tail dependence. The optimal Copula is selected based on a distance metric, such as comparing the empirical joint distribution from KDE, $F_e(x,y)$, with the candidate Copula $C_k$:
$$d_k = \left\{ \sum_{t=1}^{T} (F_e – C_k)^2 \right\}^{1/2}$$
The Copula with the smallest distance $d_k$ is chosen.
1.2 Electric Vehicle Car Charging Load Model
The aggregation of electric vehicle car charging load is a spatiotemporally distributed process influenced by multiple layers of complexity. Our model integrates the following components:
1.2.1 Dynamic Traffic Network Model: The road network is represented as a graph $G=(V, E, K, W)$, where $V$ is the set of nodes (intersections, zones), $E$ is the set of edges (road segments), $K$ is the set of time intervals, and $W$ is the time-dependent travel impedance (cost) matrix. The impedance $w_{ij}^k(t)$ for a road segment between nodes $i$ and $j$ during interval $k$ combines segment travel time $R_{ij}(t)$ and nodal delay $C_i(t)$ (e.g., at intersections).
$$w_{ij}^k(t) = C_i(t) + R_{ij}(t)$$
1.2.2 Electric Vehicle Car Agent Modeling: Different types of electric vehicle cars (private, taxi, public service) have distinct travel patterns. An Origin-Destination (OD) probability matrix $c_{ij}^{T,T+1}$ defines the likelihood of a trip starting from zone $i$ and ending in zone $j$ during time period $T$ to $T+1$. Vehicle-specific parameters (battery capacity $Cap_i$, initial State of Charge $SOC_{init}$, daily mileage, start/end times) are assigned via Monte Carlo sampling from empirical distributions. For instance, battery capacity may follow a Gamma distribution:
$$f(Cap_i; \alpha_i, \beta_i) = \frac{1}{\beta_i^{\alpha_i}\Gamma(\alpha_i)} Cap_i^{\alpha_i-1} e^{-Cap_i / \beta_i}$$
1.2.3 Energy Consumption and Charging Decision: The instantaneous energy consumption $E_{sum}$ per km is critical and depends on driving speed $v_{ij}$ (from traffic models) and ambient temperature $T_{amb}$ due to climate control usage.
$$E_{sum} = E_T + E_m, \quad \text{where } E_T =
\begin{cases}
\frac{W_c \cdot S}{v_{ij}}, & T_{amb} > T_{heat} \\
\frac{W_h \cdot S}{v_{ij}}, & T_{amb} < T_{cold} \\
0, & \text{otherwise}
\end{cases}$$
Here, $E_m$ is the base mileage consumption, $W_c$ and $W_h$ are cooling and heating powers, $S$ is distance, and $T_{cold}$, $T_{heat}$ are temperature thresholds. The $SOC$ at time $t$ updates as:
$$SOC_t = SOC_{t-1} – \frac{E_{sum} \cdot \Delta l}{Cap_i}$$
where $\Delta l$ is the distance traveled. A charging demand is triggered when $SOC_t$ falls below a driver-defined threshold upon arriving at a destination with charging infrastructure.
1.2.4 Queueing at Charging Stations: To model congestion, we treat charging stations as $M/M/c$ queues. The average waiting time $T_{mmc}$ for an electric vehicle car arriving at a station with $c$ chargers, arrival rate $\lambda$, and service rate $\mu$ per charger is:
$$T_{mmc} = \frac{(c\rho)^{c}\rho}{c!(1-\rho)^2 \lambda} P_0, \quad \text{where } \rho = \frac{\lambda}{c\mu}, \quad P_0 = \left[ \sum_{k=0}^{c-1} \frac{(c\rho)^k}{k!} + \frac{(c\rho)^c}{c!(1-\rho)} \right]^{-1}$$
This waiting time can influence user behavior and the final temporal distribution of the load.
1.3 Conventional Generator and Load Models
Conventional thermal generators are modeled using a two-state (operating/forced outage) model with forced outage rate (FOR) $a_P$. Nodal loads are treated as time-varying Gaussian random variables. The dynamic load $W_l(t)$ at a node is:
$$W_l(t) = P_l(t) + \Delta_l, \quad \Delta_l \sim \mathcal{N}(0, \sigma_l^2)$$
where $P_l(t)$ is the forecasted base load profile and $\Delta_l$ is a Gaussian random fluctuation.
2. Dynamic Probabilistic Load Flow Solution Methodology
The DPLF aims to find the probability distributions of system state variables (bus voltages) and branch flows over a time horizon, given the correlated stochastic inputs. The solution involves three key steps: correlation handling, cumulant calculation, and series expansion.
2.1 Handling Input Correlations via Nataf Transformation and SVD
Samples of correlated wind/PV power ($\mathbf{X}$) are generated using the fitted Copula model. For the cumulant method, inputs must be independent. The Nataf transformation is used to map correlated non-Gaussian variables $\mathbf{X}$ to independent standard normal variables $\mathbf{S}$.
First, each variable is transformed to standard normal via its marginal CDF $F_i$:
$$s_i = \Phi^{-1}[F_i(x_i)]$$
where $\Phi^{-1}$ is the inverse standard normal CDF.
The correlation matrix $\mathbf{\rho_S}$ for $\mathbf{S}$ is related to the original correlation matrix $\mathbf{\rho_X}$ through an integral equation. Often, an empirical approximation is used. Since $\mathbf{\rho_S}$ must be positive definite for the subsequent step, we perform a Singular Value Decomposition (SVD) to ensure robustness:
$$\mathbf{\rho_S} = \mathbf{B} \mathbf{B}^T$$
Finally, the independent samples $\mathbf{S}$ are transformed back to the correlated standard normal space $\mathbf{Y}$:
$$\mathbf{Y} = \mathbf{B} \mathbf{S}$$
These are then transformed back to the original space to obtain independent samples $\mathbf{X’}$ that preserve the correlation structure when aggregated for cumulant calculation.
2.2 Cumulant Method and Gram-Charlier Expansion
The core of the analytical DPLF is the linearization of the power flow equations around the expected operating point:
$$\Delta \mathbf{Z} \approx \mathbf{T} \cdot \Delta \mathbf{S}$$
where $\Delta \mathbf{Z}$ is the vector of output deviations (voltage magnitudes, angles, branch flows), $\Delta \mathbf{S}$ is the vector of input injection deviations (from wind, PV, electric vehicle car load, conventional load), and $\mathbf{T}$ is the sensitivity matrix obtained from the Jacobian.
The cumulant method exploits the properties that the $k$-th order cumulant of a sum of independent random variables is the sum of their $k$-th order cumulants, and that for a linear transformation $\mathbf{Z}=\mathbf{T}\mathbf{S}$, the $k$-th order cumulant of $\mathbf{Z}$ is $\kappa_k(\mathbf{Z}) = \mathbf{T}^k \kappa_k(\mathbf{S})$. The cumulants of the input injections (up to a desired order, e.g., 5th or 6th) are calculated from their respective hybrid probability models. The total cumulants for the output variables are then computed by summing the contributions from all independent input sources.
Finally, the PDF of an output variable $z$ is reconstructed using the Gram-Charlier series expansion, which expresses the PDF in terms of its cumulants and the standard normal distribution $\phi(z)$:
$$f(z) = \phi(z) \left[ 1 + \frac{\gamma_1}{3!} He_3(z) + \frac{\gamma_2}{4!} He_4(z) + \frac{\gamma_3}{5!} He_5(z) + \cdots \right]$$
where $\gamma_i$ are the standardized cumulants and $He_n(z)$ are probabilists’ Hermite polynomials. This process is repeated for each time step in the study horizon, constituting the dynamic analysis.
3. Case Study and Performance Analysis
The proposed framework was tested on a modified IEEE 39-bus system. Wind farms were connected at buses 14, 26, and 27, and PV plants at buses 3, 4, and 6. Historical data from European sites was scaled appropriately. The electric vehicle car fleet comprised 12,000 private commuters, 4,000 taxis, and 4,000 public service vehicles, interacting with a 29-node traffic network mapped onto the grid topology.
3.1 Accuracy of Probabilistic Modeling
The hybrid parametric-nonparametric model for PV generation showed excellent agreement with historical data. Key statistical metrics for a sample PV plant are summarized below:
| Statistical Feature | Actual Data Value | Modeled Value | Relative Error (%) |
|---|---|---|---|
| Mean Power (MW) | 1.215 | 1.213 | 0.13 |
| Standard Deviation (MW) | 0.912 | 0.911 | 0.11 |
| Skewness | -0.156 | -0.155 | 0.58 |
| Kurtosis | 2.84 | 2.82 | 0.54 |
The Average Root Mean Square (ARMS) error between the actual and modeled cumulative distribution function (CDF) was 1.57%, confirming high fidelity.
3.2 Static Probabilistic Load Flow Validation
The accuracy of the overall DPLF algorithm (Method 3: Proposed) was validated against the benchmark Monte Carlo Simulation (MCS) with 20,000 samples. The table below shows the average and maximum relative errors for key output statistics.
| Output Variable | Statistic | Avg. Relative Error | Max Relative Error |
|---|---|---|---|
| Voltage Magnitude | Mean | 0.0252 × 10⁻³ | 0.0969 × 10⁻³ |
| Std. Dev. | 0.0202 × 10⁻³ | 0.0864 × 10⁻³ | |
| Branch Active Power | Mean | 0.0334 × 10⁻³ | 0.1857 × 10⁻³ |
| Std. Dev. | 0.0282 × 10⁻³ | 0.1680 × 10⁻³ |
The ARMS error for CDFs across all buses and branches was minimal (average 0.008%), proving the analytical method’s precision.
3.3 Impact of Wind-PV Correlation
To isolate the effect of correlation, three methods were compared: 1) Traditional cumulant method ignoring correlation, 2) Cumulant method with basic Nataf transformation, and 3) The proposed method (Copula + Nataf-SVD). While the means of voltage and power flows were largely unaffected, their standard deviations—which quantify uncertainty—changed significantly.
| Method | Avg. Increase in Voltage Std. Dev. | Avg. Increase in Branch Flow Std. Dev. |
|---|---|---|
| Method 2 vs Method 1 | 0.011% | 1.15% |
| Method 3 (Proposed) vs Method 1 | 0.025% | 2.55% |
The proposed method, with its more accurate Copula-based correlation model, predicted higher collective variability, indicating that ignoring correlation leads to an underestimation of grid volatility risks.
3.4 Dynamic Analysis of Integration Impacts
PV Integration: The dynamic PDFs of voltage at a bus near a PV plant (Bus 6) revealed a distinct pattern. During peak insolation hours (10:00-15:00), the voltage probability distribution became more dispersed after PV integration, indicating increased variability. The probability mass within the 1.0-1.4 pu range decreased from an average of 46.39% (without PV) to 32.3% (with PV), a 30.37% reduction in concentration within that band.
Wind Integration: Conversely, wind integration at Bus 14 during high-wind periods (11:00-18:00) made the voltage distribution more concentrated. The probability of voltage lying within 1.01-1.04 pu increased from 30.1% to 54.4%, an 80.7% increase, suggesting a stabilizing effect on voltage at that node during those hours.
Electric Vehicle Car Load Impact: The integration of the large-scale electric vehicle car charging fleet had a consistent downward pressure on system voltages. The 24-hour average expected voltage across most network nodes showed a decline. The mean voltage depression across all affected nodes was approximately 1.36%, highlighting the significant impact of uncontrolled charging of a large number of electric vehicle cars on grid voltage profiles.
4. Conclusion
This paper developed a comprehensive and efficient framework for dynamic probabilistic load flow analysis in modern power grids with high penetrations of correlated renewable generation and spatially distributed electric vehicle car charging loads. The hybrid probabilistic modeling of wind/PV outputs provides greater accuracy, and the use of Copula functions effectively captures their critical spatial correlations. The detailed electric vehicle car model, incorporating traffic dynamics, user behavior, ambient temperature, and queueing theory, offers a realistic depiction of charging load distribution. The solution methodology, combining Nataf-SVD transformation, the cumulant method, and Gram-Charlier expansion, proves to be highly accurate compared to Monte Carlo simulation while being computationally efficient for dynamic studies.
The case study results yield critical insights: 1) The spatial correlation of renewable resources significantly affects the uncertainty (standard deviation) of grid states and must be modeled for accurate risk assessment. 2) The integration of PV can increase voltage variability during daytime peaks, while wind may help concentrate voltage distributions during high-output periods. 3) Perhaps most importantly, the aggregate demand from a fleet of electric vehicle cars represents a substantial new load that can cause measurable voltage depression across the grid, underscoring the need for smart charging strategies to mitigate adverse impacts. This integrated DPLF framework serves as a powerful tool for system planners and operators to quantify risks and develop strategies for the secure and stable integration of a decarbonized and electrified future energy system.