Forecasting Congestions in Interconnected Power Systems with High Shares of Renewable Energy: A Probabilistic Approach using Copulas

--The integration of increasing amounts of volatile renewable energy sources (RES) induces a growing demand for reliable forecasts of congestion to support renewable curtailment and redispatch decisions. Therefore, this paper proposes a novel approach to forecast congestions in high-voltage grids with high shares of distributed photovoltaic (PV) infeed. This approach is based on a physical PV model using intra-day numerical weather prediction (NWP) input data. Subsequently, probabilistic forecasts are generated based on Kernel density estimators (KDE) and Copula, describing the multivariate spatial dependencies for the marginal distributions of forecasting and approximation errors. Finally, a probabilistic power flow (PPF) using a linearized AC version is proposed, combining the benefits of high accuracy with high computational performance. To assess and quantify the overall advantages of this approach, a case study is carried out for an existing power system. To benchmark the proposed approach comparisons to standard multivariate normal (MVN) distributions (with and without correlations), as well as to AC and DC PPFs are compared. Non-parametric probabilistic metrics are used for evaluation, comparing accuracy and precision of individual marginal distributions, and taking the spatial dependency between forecasts at different locations into account.


I. INTRODUCTION
VER the past decades the infeed from renewable energy sources into electric power systems all over the world has progressed rapidly. For instance, photovoltaic (PV) became one of the leading technologies in Germany, comprising more than 49 GW of installed capacity by 2019. High shares of RES contribute obviously to reach decarbonization goals. However, most of the RES are decentralized, rather small-scale units, connected to the lower voltage levels and providing highly volatile infeed. Therefore, system operators are presently facing new challenges, especially, in terms of operational tasks such as congestion management. One appropriate support to system operators to handle these challenges are detailed forecasts focusing not only on the infeed patterns but also on the resulting power flows in the grid.
In terms of forecasting PV generation, various methods have been proposed over the past years focusing on different application areas [1]. Generally, those methods can be divided into deterministic and probabilistic forecasts [2]. In contrast to deterministic models, probabilistic forecasts contain additional information on the uncertainty of the prediction. However, when considering probabilistic forecasts of PV generation physical, statistical and hybrid approaches are distinguished, using either parametric or non-parametric techniques [2]. In this context, statistical approaches for PV modelling provide accurate and reliable results, especially when detailed information on individual PV plants is not available but metering data instead [3]. However, in terms of power system operations, e.g., when considering renewable curtailment as part of the system operator's congestion management, detailed information on individual RES units is required [4]. Hence, physical PV models are beneficial for the system operator, predicting the individual probabilistic state of the relevant PV units.
These PV forecasts enable system operators to ensure a stable operation of the corresponding power system, i.e., to anticipate congestions before their occurrence and take countermeasures in advance. In European electricity markets these countermeasures comprise in particular redispatching, countertrading and renewable curtailment of renewable energies [5]. The U.S. wholesale electricity markets instead rely on locational marginal pricing (LMP), i.e., taking congestions into account during market clearing [6]. Regardless of the regulatory framework and market design, it is essential for system operators to predict potential congestions at an early stage, taking the increasing uncertainties into account, in order to derive the necessary countermeasures.
In this context the probabilistic power flow (PPF) offers the possibility to take those uncertainties, e.g., resulting from probabilistic forecasts or network approximations, into account. In recent years, several methods for PPF have been developed and investigated in detail, as discussed by Chen et al. [7]. Generally, present investigations in this context focus on new mathematical approaches [8], network outages [9], power system stability [10] and computational efficiency [11]. In terms of modelling stochastic dependencies of forecasting errors and uncertainties [12], either spatial [13] or temporal dependencies are considered [14]. In addition, it is rather common for power system applications to assume multivariate normal distributions when considering such dependencies [15,16].
Within this paper, a novel approach is proposed for probabilistic congestion forecasting in the presence of large amounts of renewables, in occurrence PV systems. Thereby a hybrid approach is developed to derive multivariate probabilistic RES generation forecasts. In a first step a deterministic PV forecast is carried out using numerical weather prediction (NWP) data, based on a physical PV model with pre-processing of site-specific parameters (cf. Section II). Subsequently, the marginal distributions of the forecasting and approximation errors are simulated (cf. Section III) using kernel density estimators (KDE), after introducing a concept for eliminating the heteroscedasticity resulting from the temporal dependence. Furthermore, a copula approach is used to model the spatial dependence between those errors to cope for the resulting uncertainty. Therefore, different copulas, namely the Gaussian and Student's t, are applied hereinafter and benchmarked to the assumption of multivariate normal distribution (with and without spatial dependencies). Finally, the PPF is specified (cf. Section IV), including empirically validated shift factors for power imbalances to boundary nodes. To meet real-world performance requirements, a fully linearized version is designed that is benchmarked (cf. Section V) to the AC-PPF and the commonly used DC-PPF in terms of computational efficiency and accuracy (cf. Section VI).

II. PHOTOVOLTAIC POWER FORECASTING
To forecast the net power injection of photovoltaic (PV) units connected to an electric power system, the physical PV model introduced by Schinke-Nendza et al. [17] is taken as basis. This model uses site-specific hourly NWP data as an input and provides hourly forecasts of the power generation of individual PV units. The structure of the physical PV model used by [17] is depicted in Fig. 1. Using site-specific hourly NWP data, the model provides a rather high accuracy on a single plant level compared to other approaches, e.g., [18]. However, empirical analyses indicate that available data regarding some of the model parameters, such as the nominal power and the panel orientation of PV units, may be erroneous. Therefore, a nonlinear least squares (NLS) regression is used here to obtain the best possible fit of the physical model, instead of using a statistical post processing as in [17].

III. MODELING THE SPATIAL INTERDEPENDENCE OF FORECASTING ERRORS
The cumulative net power injection at a certain feed-in location, i.e., a node of the underlying power system, can be obtained by aggregation of multiple PV power forecasts. The relative forecasting error of the PV net power injection at feedin locations ∈ = {1, … , } during timestep is defined as: where ̂, is the deterministic forecast of the physical model using NWP data, 0, is the cumulative rated power, and , is the actual measured power of the PV units at feed-in location .

A. Upscaling the PV forecasts
Considering feed-in location , there are | | PV units ∈ , connected to this node, comprising an actual rated power of 0, . In practice, the available dataset (for measurement and NWP data) yet commonly does not include all | | PV units, but rather a subset ′ of | ′ | PV units (with | ′ | ≤ | | and 0, ′ ≤ 0, ). Therefore, the deterministic forecast ̂ of the physical model as well as the forecasting error ∆ , need to be upscaled to receive reliable forecasts in terms of network congestions, taking all connected PV units into account. In terms of the deterministic forecast, the upscaling is carried out using a linear upscaling coefficient , 0 defined as the ratio of the actual cumulative rated power 0, to the rated power of the subset 0, ′ : Considering the relative forecasting error ∆ , at feed-in location , the question arises to which degree the individual forecasting errors of the single PV units are stochastically dependent. Our empirical analyses indicate the linear upscaling coefficient to perform the best, i.e., assuming a full correlation of the forecasting errors of the PV units located at a certain node. This assumption of fully correlated forecasting errors outperforms both a correlation-based upscaling and the independency assumption.

B. Elimination of Bias and Heteroscedasticity
To reduce the modelling effort in terms of the copulas and the complexity of the marginal distributions, any bias of the forecasts possibly arising from the physical model is quantified and eliminated for each hour separately. Therefore, the forecasting error is hereinafter assumed to be unbiased (notated as ∆̃, ). Yet, a specific heteroscedasticity, i.e., the time dependence of the variance, may still be observed in the relative forecasting errors ∆̃, , leading to higher deviations during the midday hours, but lower deviations during the morning and evening hours [14]. Therefore, a novel approach for reducing the heteroscedasticity of the forecasting error is introduced, using the extraterrestrial irradiance as a basis for normalization, following [19]. The unbiased normalized forecasting error ̃, is then obtained as follows: The spatial dependence for all feed-in locations is hence considered hereinafter based on the marginal distributions of ̃, using a copula approach.

C. Spatial Dependence Structure
The interdependence of the forecasting errors at multiple feed-in locations ∈ = {1, … , } can be described using a joint distribution function (̃1, … ,̃). According to [20] this function may also be written as follows: where : [0,1] → [0,1] represents a copula and describes the marginal distribution (cumulative distribution function, CDF) of the forecasting error at feed-in location . In this context, the elements of the input vector can be written as: whereby the elements ~(0,1) are standard uniformly distributed. By means of the inverse probability function −1 , the copula may be expressed as follows: for all feed-in locations ∈ .

1) Gaussian Copula
The spatial dependence structure of the forecasting errors at multiple locations can be modelled, for instance using a Gaussian copula, see [21]. Following [22], the copula : [0,1] → [0,1] is described as the multivariate normal distribution ( , ) in ℝ with the density following to: where ∈ ℝ × is a correlation matrix with diagonal elements equal to one, ∈ ℝ × the identity matrix, = ( −1 ( 1 ), … , −1 ( )) the vector of inverse cumulative standard normal distribution −1 , and the input vector of the Gaussian copula according to (5). The correlation matrix can be estimated approximately, following [22], giving an estimator ̂ close to the maximum likelihood estimation (MLE).

2) Student's t Copula
The spatial interrelation of the forecasting errors at multiple locations can also be modelled using a t copula. Following [23], this type of copula surpasses symmetrical copulas, such as Gaussian copulas, when modelling tail-dependencies of the underlying marginal distributions. Based on the multivariate Student's t distribution for a set of random variables = ( 1 , … , ) ∈ ℝ , written as ~, ( , ), with degrees of freedom, a mean vector equal to and a correlation matrix as introduced above, the density of the t copula is given by: where Γ( )is the gamma function, the input vector of the t copula with elements = −1 ( ) for ∈ and represents the CDF value of the forecasting error at location ∈ according to equation (5), see [23]. Following [22], the correlation matrix estimator ̂ can be obtained using Kendall's tau estimates in a first step. In a second step, the remaining parameter , i.e., the degree of freedoms, can be estimated by fixing ̂ and applying an MLE to maximize the log-likelihood function for the pseudo-sample, which is defined as with ∈ {1, … , } and observations:

D. Marginal Distributions
The marginal distributions (̃) of the forecast errors at the feed-in locations ∈ , can be estimated using a nonparametric approach. Subsequently a kernel density estimation (KDE) is applied to the marginal distributions, as proposed by Chai et al. for the forecasting errors of PV generation [14]. Using the sample of the forecasting error ̃= (̃, 1 , … ,̃, ) at location ∈ with ∈ {1, … , } and observations, the kernel density estimator using a gauss kernel is defined as: where is the density of the standard normal distribution, ̃, corresponds to the mean of the -th kernel function and ℎ is the bandwidth of the KDE [24]. The bandwidth ℎ , representing a free parameter of the kernel density estimation, can be determined analytically, for instance using the normal distribution approximation following [25]. In this case the optimal estimator of the bandwidth can be obtained as: in case that the forecasting errors are unbiased.

IV. PROBABILISTIC POWER FLOW ANALYSIS
A network consisting of nodes (not necessarily equal to the number of feed-in locations ) and lines is investigated, assuming a perfectly balanced power system, hence focussing on a single-phase representation. To determine the power flow along the transmission lines of the network, the equations below can be used: = ∑ ( , sin , − , cos , ) =1 , ∈ where is the net active power injection and is the net reactive power injection, both at node ∈ , is the voltage magnitude, , = − is the difference of voltage angles between node and node , and , = , + , is the admittance of the line from node to node , consisting of the conductance , and the susceptance , , following [7]. In addition, the line parameters are assumed to be time-independent and to not vary with the weather conditions (e.g., temperature, wind speed etc.).

A. Network Congestions
The constraints in a power system for steady-state operation typically depend on the values of the current magnitudes on the transmission lines ∈ ℳ and the voltage magnitudes at the nodes ∈ . Assuming the voltages to be uncritical, congestions in a power system exclusively depend on the transmission line currents. Therefore, the following constraints needs to hold: whereby the maximum permissible current of line may contain some additional security margins. For ease of notation, the time index is dropped here and subsequently.

B. Modelling Boundary Nodes to External Networks
To investigate the impact of the forecasting error on the congestions in an interconnected network, the effects of adjacent networks need to be considered. In case of a forecasting error different from zero, the power imbalance is at least partly shifted to the boundary nodes ∈ connecting the considered network with external networks. Therefore, a multivariate linear regression (MVLR) on available measurement data is proposed in this paper. Building on the adjustment matrix introduced in [26], it models the dependency of power deviations at the boundary nodes from the inner nodes that connect decentralized PV units and loads. The MVLR is specified as follows: where , is the adjustment factor at boundary node given changes at inner node . The adjustment matrix summarizes the shifts in power flows at the boundary nodes as a (linear) function of the changes at the infeed and load nodes: Furthermore, the error of this approximation is denoted as and described by a time-invariant marginal distribution ( ). The overall uncertainty at boundary node given a balanced power flow schedule and some PV forecast errors ̃ can then be described as a superposition of the approximation error and the forecasting errors of the feed-in locations ̃ weighted by the adjustment factors, respectively: To estimate marginal distributions of the approximation error , the changes in power balances between two consecutive timesteps are considered. In the event of forecasting errors, the same adjustment matrix may be used to describe the corresponding changes in the boundary power flows, assuming the initial power flow to be known perfectly for the expected values of the PV power injection. The additional approximation error , then comes on top of the forecasting errors. Since both types of errors have different sources and are derived from different datasets (approximation error: time series of observed power flow changes; forecasting error: mismatches between prediction and measurement for a certain time step) the errors are appropriately deemed to be stochastically independent.

C. Linearized Power Flow
The AC power flow defined by equations (12) and (13) describes a non-linear problem due to the quadratic dependence on the voltage magnitudes and the sinusoidal dependence on the voltage angles. To substantially improve the computation time of the PPF, a linearized AC power flow is used hereinafter.

1) Network Equations
In this context, a certain state of operation in an upcoming interval (index: *) is estimated based on the expected initial state of operation (index: 0), using a Taylor series: where is the Jacobian matrix, consisting of the gradients of the active and reactive power differentiated with respect to the voltage angle and the voltage magnitude , see [8]: Hence the inverse Jacobian matrix can be used to determine the change in the voltage magnitudes ∆ and angles ∆ following a change in the active and reactive powers ∆ and ∆ .

2) Line Current
In addition, also equation (18) can be linearized for the initial state of operation, using the expected values of net power injections, as described above. The current on transmission line ∈ ℳ depends on the voltage magnitudes and angles at the beginning of the line ∈ and the end of the line ∈ (written as ↔ , where ≠ holds). The first-order Taylor approximation for the current on line ∈ ℳ writes Where the elements , of the matrix are defined as

3) Probabilistic Power Flow
Assuming the loads in the power system to be forecasted with a high accuracy, the uncertainties in the net power injection of the PV units ∆ , at feed-in location ∈ are larger than the uncertainties in the load ∆ , : ∆ , ≫ ∆ , .. Furthermore, it is assumed that most of the considered PV units are located in the medium and low voltage levels (without any reactive power management as described in [27]), hence not affecting the reactive power exchange between voltage levels. Therefore, the change in the line currents ∆ = * − depends on the changes of the net power injection of the PV units, i.e., ∆ , = * − ,0 for ∈ , exclusively. Hence, the overall linearized equations for obtaining the transmission line currents based on equations (18) and (21) can be formulated as: with the inverse Jacobian matrix being segmented into the four × matrices −1 to −1 : The problem of solving multiple non-linear equations thus becomes linear and the solution is obtained directly through matrix operations on the change of the active power net injection of the PV units at the different feed-in locations. The PPF then is carried out using the Monte Carlo simulation introduced below.

D. Monte Carlo Simulation
Within this paper a Monte Carlo simulation (MCS) is proposed, following [28]. Hereinafter the MCS is used to stochastically predict network congestions by applying the probabilistic forecasts of the PV generation as input variables to the power flow analysis. Since finding a closed-form analytical solution for the given mathematical problem may either not be feasible or lead to an over-sophisticated solution, the application of the Monte Carlo method is deemed more appropriate, relying on the central limit theorem and the law of large numbers [28]. In the application, the MCS will be applied to all alternatives: The AC power flow, the linearized AC power flow and the DC power flow. Finally, the different approaches are compared using the performance metrics and benchmark model described in the following section, to evaluate the improvements of the overall approach.

A. Probabilistic Metrics
To evaluate the goodness of fit of the spatial dependency and to assess the accuracy of multivariate forecasting models predicting the D-dimensional variable for a certain timestep t, there are two commonly used probabilistic measures: The energy score and the variogram-based score [21]. 1) Energy Score The energy score (ES) represents a proper multivariate generalization of the univariate continuous ranked probability score and is applicable to evaluate the ensembles of forecasts of the variable = ( ,1 , … , , ) : where describes the sample size and ̂( ) is the vector of forecasts for sample . The ES is a negative scoring rule; hence, lower values indicate a higher forecasting accuracy. However, investigations of Scheurer and Hamill [29] have shown that the ES is not adequately sensitive to incorrectly specified correlations of multivariate forecasts. It is worth mentioning that the ES is able to detect mean-biased forecasts correctly, but it is not able to indicate incorrect dependency structures and false variances clearly [29].
2) Variogram-Based Score Especially for evaluating the spatial dependence structure and the spread of multivariate forecasts for the D-dimensional variable , the variogram-based score (VS), according to [29], is considered within this paper as a proper scoring rule: (27) where is the order of the VS (where [29] proposes 0.5 for accurate results) and , is a weighting factor. Just as the ES, this scoring rule is also negative oriented; hence, lower values indicate a higher forecasting accuracy. 3

) Diebold-Mariano Test
In addition, the Diebold-Mariano (DM) test can be carried out on VS and ES, according to [30], for testing the null hypothesis that the performance of two considered models is equal, assuming a standard normal distribution to hold for the test statistic, i.e., ~(0,1). For a given the significance level α, the DM test statistic furthermore indicates the preferred model [30].

B. Benchmark Model
Since the probabilistic performance metrics provide a comparative assessment of different models, a simple benchmark model is introduced hereinafter to assess and quantify the improvement resulting from the overall approach introduced in the sections above.

1) Multivariate Normal Distribution
As a simplified benchmark model, a multivariate normal distribution (MVN) is introduced, following [21]. This allows to identify the benefits of applying the more complex and numerically demanding copula approach for modelling the forecasting errors at multiple feed-in locations and their spatial dependences. A further simplification is to assume the forecasting errors of the multiple feed-in locations to be spatially independent. For both approaches (with and without spatial dependence) the relative forecast errors ∆ can be described as ∆~ℳ ( , ). Hence, for each feed-in location and each hour a set of parameters exists: Whereas the mean relative forecasting error , should be zero for unbiased forecasts, the variance of the forecasting error , 2 and the covariances to the other feed-in locations are summarized in the variance-covariance matrix . When neglecting the spatial dependence between the forecasting errors in the second benchmark, non-diagonal elements become zero.

2) DC Power Flow
The DC power flow is based on the assumption that the voltages at all nodes and are close to their nominal value, the phase angle differences , are small and losses can be neglected, as described in [6]. Therefore, the relation between nodal net power injections and the power flow across the considered transmission lines can be described using the power transfer distribution factor (PTDF) matrix : Based on the previously introduced assumptions = holds for per unit values.

VI. APPLICATION
Within this paper an existing power system on the high voltage level (110 kV) is considered as a case study, see Fig. 2. This network consists of 18 buses and 17 nodes, of which 12 nodes represent PV feed-in locations, accounting for an installed capacity of 527.5 MW. The nodes of the considered power system are serially numbered (1 to 17), whereby the boundary nodes connecting to other networks (nodes 4, 6, 14 and 17) are highlighted with arrows. For the power system, measurement data for each node and each transmission line are available in a quarter-hourly temporal solution for one year (here: 2017), i.e., voltages, active and reactive power. Regarding the net power injection of the PV units, measurement data of the active power is available for a sample of PV units comprising 156.2 MW of installed capacity, see Table I.  The measurement data is available in an hourly temporal solution for three months of three years (May to July, 2015 to 2017). For each PV unit, hourly NWP data (global irradiance and temperature) is available in a forecast horizon of 1-3 hours, i.e., the forecasts are updated in intervals of 3 hours and contain information on the next couple of hours.

A. Evaluation Methodology
In a first step, a reference current is calculated for each line based on the available measurement data. Furthermore, the load, the grid voltages and the transfer power to external networks are assumed to be known exactly, i.e., the power transfer at boundary nodes solely varies due to forecasting and related approximation errors. All models have been trained on the available data of the years 2015 and 2016. In addition, the adjustment matrix and the approximation error for the boundary nodes are determined for January to April 2017 exclusively. Hence, an out-of-sample test of all models is carried out for 2,015 hours of May to July 2017. The MCS hereinafter is based on a sample size of 5,000. Subsequently, the results are evaluated as follows: the KDE-based gaussian copula (GC), the KDE-based student's-t copula (TC), the multivariate normal distribution with correlation (MVN) and without correlation (UMVN) are combined with the different PPF approaches, namely the linearized AC (LAC) PPF as well as AC and DC PPF, as introduced above.

B. Results
In the following, the results of the case study introduced above are presented. Therefore, the sensitivity to the probabilistic PV forecasts, the upscaling method of the forecasting errors, the power flow method and the approximation error are considered. In addition, the computational efficiency of the proposed model is scrutinized.

1) Quality of Probabilistic Forecasts
The multivariate probabilistic PV forecasts based on the UMVN, MVN, GC and TC models are compared in the upper left panel of Fig. 3. Thereby, the average variogram-based score with = 0.5 (VS1: wj,l = 1; VS2: wj,l = ρj,l, where ρj,l is the correlation coefficient of the reference current) and the average energy score (ES) are compared. As depicted in Fig. 3 there are significant differences in VS1, VS2 and ES when considering the PV forecasts, indicating the TC to surpass the UMVN, MVN and GC approaches. When using these forecasts for the AC and LAC PPF, eight test cases are obtained. Thereby, the TC still outperforms the other approaches in VS1 and VS2. However, the results for the ES indicate that the MVN and the TC might be mean-biased. Yet the ES is not indicating the improvements of modelling the dependency structure (and the KDE of TC) correctly as emphasized in [29]. Subsequently, the VS is hence considered exclusively. Since the differences of the scores appear to be rather low for both line current forecasts (AC and   AC-VS1 -13.5 ** -4.01 ** -8.78 ** -0.29 ++ -5.69 ** -5.8 ** AC-VS2 -15.3 ** -2.47 ++ -7.99 ** 1.67 ++ -3.87 ** -5.91 ** LAC-VS1 -9.37 ** -9.43 ** -9.9 ** -5.94 ** -7.24 ** 0.38 ++ LAC-VS2 -9.33 ** -11.21 ** -11.23 ** -7.59 ** -8.45 ** 1.61 ++ In this context, the following notation is used indicating the confidence level on which a null hypothesis H0 of equally performing models is rejected: No asterisk means that the null hypothesis is not rejected, plus sign (+) confidence = 95%, oneasterisk (*) confidence = 99% and two-asterisks (**) confidence = 99.9%. Solely the MVN reveals equal performance to the TC in the AC-VS2 test case

2) Sensitivity to Upscaling of PV Forecasts
To assess the advantages of the proposed upscaling method of the PV forecasts assuming full correlation, a comparative analysis is carried out using the VS. In this context, the results of the upscaling of the PV forecasts assuming full correlation of the forecasting errors is benchmarked to a correlation-based approach, using the average correlation ̅ of the forecasting errors for each node, and to an upscaling assuming independent forecasting errors, i.e., ̅ = 0. For this comparison, the results are given in Table III for the full correlation (FC), the correlation-based upscaling (CB) and the independency assumption (IA). Accordingly, the FC upscaling method (assuming fully correlated forecasting errors at each node) surpasses the correlation-based and independency upscaling methods, for both TC-AC and TC-LAC for VS1 (wij=1) and VS2 (wij=ρij). Although the small differences in scores of the FC and CB approaches seem to indicate a comparable forecasting accuracy, the DM test shows that the results are statistically significant at the 99.9% confidence level.

3) Sensitivity to Power Flow Method
Using the DC PPFs as benchmarks for the AC and LAC PPFs, Table IV   As already expected, the results of the PTDF based DC PPF do not come close to the LAC or AC power flow. In opposite, the VS1 and VS2 results of the TC-LAC are nearby the UMVN-AC results, as illustrated by Fig. 3. When applying a DM test to compare both approaches, i.e., TC-LAC and UMVN-AC, the resulting test statistic are -0.1806 for VS1 and 4.3981 for VS2, respectively. Therefore, both approaches provide equal performance in terms of the VS1, but UMVN-AC surpasses the TC-LAC approach for VS2.

4) Sensitivity to Approximation Errors at Boundary Nodes
The approximation error (AE) can be modelled using the same approaches as used for the forecasting errors. Therefore, each of the approaches TC to UMVN has been carried out for AC-and LAC-PPFs, with and without AE., see Table V. As depicted by Table V, the differences in VS1 and VS2 are rather small, correspondingly the DM test indicates equal performance when comparing the models with AE with those models without AE. In the GC test case, the model without AE even surpasses the model with AE, albeit the difference is not statistically significant.

5) Computational Efficiency
The case study has been carried out using MATLAB / MATPOWER on a Windows Server 2012 R2 Standard with 64bit operating system, a two kernel 3.50 GHz processor (Intel ® Xeon ® CPU E5-2637 v2) and 192 GB RAM of memory. Table  VI compares the computational efficiency for the different PPF approaches. Thereby, the MCS of the AC-PPF have been carried out sequentially and parallelized. The results for the computation time of the parallelized MCS AC-PPF opposed the expectations: The additional effort for distributing and allocating large variables (type: 2,015 hours × 19 lines × 5,000 samples, i.e., 1.53 GByte per variable) to the virtual kernels, exceeds the gains of parallel computation. Furthermore, the gains of the LAC-PPF are remarkable, taking less than one minute for computing more than 10 6 scenarios.

VII. CONCLUSION
This paper proposes a novel approach to provide probabilistic congestion forecasts in high voltage power systems utilizing multivariate copula-based RES-infeed simulations. The improvements presented contribute to increase both, the performance of the congestion forecasts (upscaling of forecasting errors, reduction of heteroscedasticity, modelling of boundary nodes) and the computational efficiency (application of linearized power flows). The benefits of the proposed model are demonstrated on a real-world case study against commonly used benchmark models. Thereupon, this approach enables system operators to predict potential congestions at an early stage, taking the increasing number of uncertainties into account, in order to derive the necessary countermeasures.