Robust DC optimal power flow with modeling of solar power supply uncertainty via R-vine copulas

We present a robust approximation of joint chance constrained DC optimal power flow in combination with a model-based prediction of uncertain power supply via R-vine copulas. It is applied to optimize the discrete curtailment of solar feed-in in an electrical distribution network and guarantees network stability under fluctuating feed-in. This is modeled by a two-stage mixed-integer stochastic optimization problem proposed by Aigner et al. (Eur J Oper Res (2022) https://doi.org/10.1016/j.ejor.2021.10.051). The solution approach is based on the approximation of chance constraints via robust constraints using suitable uncertainty sets. The resulting robust optimization problem has a known equivalent tractable reformulation. To compute uncertainty sets that lead to an inner approximation of the stochastic problem, an R-vine copula model is fitted to the distribution of the multi-dimensional power forecast error, i.e., the difference between the forecasted solar power and the measured feed-in at several network nodes. The uncertainty sets are determined by encompassing a sufficient number of samples drawn from the R-vine copula model. Furthermore, an enhanced algorithm is proposed to fit R-vine copulas which can be used to draw conditional samples for given solar radiation forecasts. The experimental results obtained for real-world weather and network data demonstrate the effectiveness of the combination of stochastic programming and model-based prediction of uncertainty via copulas. We improve the outcomes of previous work by showing that the resulting uncertainty sets are much smaller and lead to less conservative solutions while maintaining the same probabilistic guarantees.


Introduction
The proportion of renewable energy, such as solar and wind energy, in electrical distribution networks is constantly increasing. Due to these difficult to predict and highly fluctuating energy sources, the operational management of electrical networks becomes very challenging. Transmission system operators (TSO) have to control the feed-in and the power distribution in the network and have to meet safety requirements at the same time. If the network risks a system overload, the feed-in from renewables must be curtailed. However, the curtailed energy has to be minimized for financial and ecological reasons. Therefore, there is a high demand for the combination of advanced forecasting and optimization models. In this work, we show how these models can be applied and combined for the optimal curtailment of solar feed-in in an electrical distribution network. The predominantly used model and optimizing the production and distribution of power in an electrical network is the Optimal Power Flow (OPF) model. In its classic version this is a non-linear non-convex optimization problem which is hard to solve and was originally introduced in Carpentier (1962). For a broad overview of the literature on OPF, we refer to Frank et al. (2012a) and Frank et al. (2012b). Due to the computational difficulty of the OPF problem, there are some approximation approaches in the literature. One of the most frequently used approximations is the DC Optimal Power Flow (DC OPF), see Christie et al. (2000). It results in a power flow model including only linear constraints and can be solved efficiently with standard software. For the optimization of power grids under uncertainty the DC OPF model is also used in this work.
In applications to power grids, it is important to ensure that there is a sufficiently high probability (chosen beforehand) that all safety constraints like transmission limits are satisfied. This can be modeled with a two-stage stochastic optimization model incorporating joint chance constraints that enforce the simultaneous satisfaction of several constraints with a predefined probability. In the first stage, the nominal network operating solution, including generator output, (discrete) curtailment, power flows and voltage angles, is decided before the realization of uncertainty is revealed (here-and-now). After the uncertain parameters manifest themselves, the two-stage variables react to them. In the second stage, the network response to fluctuation ensures that there is a high probability of transmission limits being maintained. From a practical perspective, protection through probabilistic constraints is suitable because short-term overloads in the electrical network are acceptable. In the event of larger or longer lasting overloads, countermeasures will need to be taken, where a TSO will need to (re-) optimize interventions in order to stabilize the network. In our model, curtailment limits the output of renewable power production to a specific percentage proportion of the installed power.
We approximate the probabilistic constraints in the optimization problem using robust constraints within a robust safe approximation, see Nemirovski (2012). By a suitable choice of the uncertainty set we can ensure that all robust feasible solutions are also feasible for the stochastic optimization problem. The constraints of Robust DC optimal power flow with modeling of solar power supply… the robust approximation thus lead to sufficient conditions for the chance constraints being satisfied. In particular, we use a mixed-integer linear reformulation for the approximation introduced in Aigner et al. (2021). Hence, by solving only one mixed integer optimization model to global optimality, a robust solution is computed that is feasible for the chance constrained problem. The respective uncertainty sets are computed with the procedure proposed in Margellos et al. (2014) based on the scenario approach (see Calafiore and Campi (2005)) of stochastic optimization, which uses samples from a suitably chosen probability distribution. The present paper proposes several enhancements of our previous work, which consist in the utilization of R-vine copulas (see e.g. Joe (2015)), a flexible parametric model to construct multivariate probability densities by decomposing them into several bivariate conditional (and univariate) densities to fit distributions to available data. Note that R-vine copulas contain the family of D-vine copulas as special case, which we used so far to model data from meteorology and solar power supply, see Schinke-Nendza et al. (2021);. From the fitted R-vine copula model we draw samples in order to obtain the uncertainty sets with the help of the scenario approach mentioned above. Then, in a second step, we modify the R-vine copula model such that we can draw samples from conditional distributions. This allows us to determine uncertainty sets depending on weather forecasts provided by DWD (German Meteorological Service) which are significantly smaller and lead to a drastic reduction of conservatism and less costly curtailment with same probabilistic guarantees.
There are many research activities regarding OPF under uncertainty. The goal is to determine an optimal network configuration that remains feasible under uncertainty where the approach considered in this paper uses methods and models from stochastic and robust optimization. We refer to Ben-Tal et al. (2009) and Prékopa (1995) for a broad overview of these two paradigms regarding optimization under uncertainty. Note that due to the non-convexity of the nominal AC OPF, only solutions that are approximately protected against uncertainty can be computed as in (Dall'Anese et al. 2017;Roald and Andersson 2018;Zhang and Li 2011) with robust or probabilistic constraints.
Essential for an algorithmically tractable treatment of uncertainty in optimization problems is the possibility to solve the underlying deterministic problems (without uncertainty) efficiently. This is why the linear DC OPF model is suitable and of great interest. Such uncertain optimization problems are usually solved by reformulating them under specific assumptions on the underlying probability distribution or by using approximation techniques from stochastic programming. Most chance constrained OPF problems considered in the literature have separate chance constraints for each engineering limit, including both generation and transmission limits. For example, the authors of Bienstock et al. (2014); Lubin et al. (2016) focus on OPF with individual probabilistic constraints under Gaussian distributions. Uncertainty probabilities for specific classes of probability distributions are considered robustly in Roald et al. (2015); Xie and Ahmed (2018). Furthermore, there is a limited number of papers that deal with joint chance constraints OPF models. They allow much stronger system security guarantees, but are much harder to solve, see Geng and Xie (2019). Most common solution methods are based on the Boolean or Bonferroni approximation (see e.g. Jia et al. (2021)) and on scenario approximations (see e.g. Peña-Ordieres et al. (2021)).
In addition, the curtailment of renewable power is used in practice to reduce the feed-in of renewable energy sources, maintaining network stability and avoiding overloads of transmission lines. The curtailment of uncertain feed-in from renewables has also been considered in several OPF models. Examples can be found in (Aigner et al. 2021;Roald et al. 2016;Qiu and Wang 2014;Wang et al. 2011;Dall'Anese et al. 2017). Note that there are two principal types of curtailment strategies, which are usually modeled by additional discrete or continuous decision variables or fixed parameters. The first and more common type of curtailment uses output capacities, which restrict the maximum possible power input. This limit cannot be exceeded and any potential power production above the limit is cut off. The second type of curtailment reduces the produced energy by a fixed value regardless of how high the feed-in amount is. Chance constraints in combination with curtailment are usually tackled by sampling techniques from stochastic optimization already mentioned above. In the present paper we use discrete curtailment levels as it is common practice in many industrial applications and set by law in Germany.
To construct parametric models for multivariate distributions, vine copulas are a versatile tool which has been used in the literature for similar problems. For example, in Guo et al. (2021); Khuntia et al. (2019); Xiao et al. (2020), copulas are applied for dependency modeling of wind power in conjunction with OPF. Furthermore, in Xu et al. (2021), Gaussian copulas are used to determine uncertainty sets for an OPF problem with chance constraints.
The main contribution of the present paper is an extension of the safe approximation of the joint chance constrained DC OPF model introduced in Aigner et al. (2021), by combining it with a model-based prediction of solar power supply via copulas. Furthermore, additional information gained from weather data can be integrated into the copula approach and thus conditional distributions of solar power supply can be modeled. However, with regard to conditional sampling, vine copulas have some restrictions as described in Cooke et al. (2015), i.e., when drawing conditional samples from a given vine copula model, only some components of the underlying vector data can be taken into account in the conditioning set. To resolve this issue, various algorithms for conditional sampling from D-and C-vine copulas have been considered in the literature, see Bevacqua et al. (2017). In the present paper we propose a modification of the fitting procedure for the more general class of R-vine copulas. This modification allows us to obtain a suitable R-vine copula for any set of components on which we want to condition. To the best of our knowledge, this modification has not yet been considered before. This rest of this paper is structured as follows. Section 2 recalls the joint chance constrained DC OPF model considered in Aigner et al. (2021), together with its robust approximation using box uncertainty sets. Then, in Sect. 3, the modeling of the underlying multivariate probability distribution with the help of R-vine copulas is introduced, where suitable uncertainty sets are constructed via the novel combination of the scenario approach and the fitted R-vine copulas. The numerical results of case studies based on real-world data for the distribution network of N-ENER-GIE GmbH are presented in Sect. 4. They demonstrate the benefit of combining stochastic programming with a model-based prediction of uncertainty via copulas. The computed solutions are robust and lead to relatively small cost increase compared to the nominal optimization model that ignores uncertainty. The consideration of conditional probability distributions further improves the solution quality. Finally, Sect. 5 concludes.

Chance constrained DC optimal power flow model
In this section, we recall the chance constrained DC optimal power flow model with the possibility to curtail feed-in proposed in Aigner et al. (2021), which is based on Bienstock et al. (2014).

Nominal DC optimal power flow with curtailment
We model the electrical distribution network as an undirected graph G = (N, L) where N = {1, … , n} for some integer n > 1 represents the set of vertices and L ⊆ N × N denotes the set of edges. In the context of power system optimization, vertices are also called nodes or buses, and edges are called (transmission) lines. The set of those nodes that are connected with (continuously controllable) slack generators of higher network hierarchies is denoted by N G ⊆ N . Furthermore, for each k ∈ N we denote the set of adjacent nodes with N(k) ⊆ N . For notational ease, we assume that every node is connected to (discretely) controllable solar power generation units. The energy production on a bus without solar feed-in is set equal to zero.
In order to control the solar feed-in, discrete regulation decisions can be made at each node. Curtailment is realized by restricting the maximum feed-in to a certain fraction vector = ( 1 , … , n ) ∈ S = S 1 × … × S n ⊂ [0, 1] n of the installed capacity vector Note that the installed capacity is the intended full-load sustained solar energy production at each node. In practice, sets of curtailment factors with a small number of levels are common. Typical sets of curtailment factors for single nodes are {0, 0.3, 0.6, 1.0} or {0, 0.1, 0.2, … , 1.0}.
Thus, at a node k ∈ N , the power fed into the network cannot exceed k P I k . Any potential feed-in above this value is cut off. We model the curtailed uncertain solar feed-in P in k based on a given solar power production P PV k ≥ 0 via i.e., P in k = min{P PV k , k P I k }. In the following, we briefly recall the DC optimal power flow model with discrete curtailment of solar feed-in proposed in Aigner et al. (2021), where Table 1 summarizes the notation used for decision variables and input parameters.
The equality constraints (1b)-(1d) model the active power flow, which is determined by the power flow equations (1d) and Kirchhoff's first law where we distinguish the two cases with and without generators, see (1b) and (1c) respectively. Note that the power at each node has to be balanced. This means that at each node k ∈ N the active power production P G k + P in k ∈ [0, ∞) from generators and renewables equals the demand (1a) min 1 3 Robust DC optimal power flow with modeling of solar power supply… P D k ≥ 0 plus the active power sent to adjacent nodes ∑ l∈N(k) p kl ∈ R . The active power flow on transmission line (k, l) ∈ L is the product of voltage angle difference k − l ∈ [−2 , 2 ] and susceptance b kl > 0 . At the same time, the transmission limits considered in (1e) must not be exceeded. The vector of generator outputs P G can be continuously controlled within the generator bounds considered in (1f). Furthermore, we assume that there is a bus k 0 ∈ N with a reference angle k 0 = 0.
The optimization task consists in minimizing the objective function given in (1a) which is the sum of power generation costs ( f k ) and curtailment costs ( c k ) subject to the constraints mentioned above. Note that the functions f k for all k ∈ N G and c k for all k ∈ N can be assumed to be linear or convex quadratic in the generator output. Since the minimum expressions in (1b) and (1c) can be linearized by introducing auxiliary variables and additional linear constraints (see e.g. Sherali and Adams (2013)), the optimization problem considered in (1) is a mixed-integer linear or convex quadratic program and can be solved efficiently to global optimality with standard techniques and software using, e.g., the Gurobi optimizer [23].

Uncertainty modeling
In practice, the vector of solar power production P PV = (P PV 1 , … , P PV n ) ∈ [0, ∞) n is not known in advance. In addition, the production of renewable power can be subject to high fluctuations and is therefore an uncertain quantity. Using a network operating strategy that is computed by ignoring such uncertainties, a sudden fluctuation of renewable energy can lead to overloads in the electrical network. In the worst case, this can lead to failure of network elements owing to cascade effects. To prevent this, the optimization model explained in Sect. 2.1 has to be extended in order to take such fluctuations into account, and individual feed-in units may have to be regulated. In particular, we model the vector of produced solar power P PV as the sum of a vector P F = (P F 1 , … , P F n ) ∈ [0, ∞) n of forecasted solar power and a random fluctuation vector X = (X 1 , … , X n ) ∶ Ω → R n defined on some probability space (Ω, F, ℙ) , i.e., However, in a first step, we need to determine a nominal operating solution (P G , , p) together with a curtailment decision that is feasible for the nominal feed-in vector P F (corresponding to X = 0 ), i.e., the decision variables P G , , p, have to fulfill the constraints (1b)-(1e), where P PV is given in (2) with X = 0 . In addition, we require that, with high probability, the network reaction to fluctuating feed-in remains feasible, see the chance constraint given in (6g) below. To model this kind of network reaction, we consider randomized duplicates P G,X ∶ Ω → [0, ∞) |N G | , X ∶ Ω → [− , ] n and p X ∶ Ω → R |L| of the decision variables P G , , p introduced in Sect. 2.1, which depend on the realizations X( ) for ∈ Ω of the random fluctuation vector X. Note that realizations X( ) ≠ 0 of X may lead to a changed distribution of power in the network and, therefore, to an imbalanced network. The generators then change their output to P G,X( ) in order to balance the total active network power. Furthermore, the decision variables X and p X are adjusted correspondingly to ensure feasibility.
(2) P PV k = P F k + X k for all k ∈ N.

3
Thus, in the setting of the two-stage stochastic optimization problem described above (see also Sects. 2.3 and 2.4), the variables P G , , p refer to first-stage (or here-and-now) decisions. They must be decided for the nominal feed-in vector P F (corresponding to X = 0 ), before uncertainty is revealed. Moreover, for fixed first-stage variables P G , , p , any realization X( ) ≠ 0 of X leads to a reaction of the network by choosing optimal second-stage (or wait-and-see) variables P G,X( ) , X( ) , p X( ) , where we assume that the power generation is balanced by the Automatic Generation Control Borkowska (1974). This means that the total power generation mismatch The vector of decision variables X is adjusted in a way that the power balance equations are fulfilled for each ∈ Ω . Furthermore, for each ∈ Ω we put It can be shown, see Aigner et al. (2021), that for each realization X( ) of X the equation system given in (4a)-(4b) has a uniquely determined solution X( ) , i.e., the wait-and-see variables P G,X( ) , X( ) , and p X( ) are uniquely determined by (3), (4a)-(4b), and (5).

Chance constrained DC optimal power flow
By construction, the vectors p X and P G,X of power flows and generator outputs are random variables that depend on the realization X( ) of the random fluctuation vector X and on the values of first-stage decision variables P G , , p, . Thus, we are searching for solutions (P G , , p, ) which satisfy the limits of type (1e) and (1f) for power flows and generators outputs, respectively, with a probability of at least 1 − for some small number ∈ [0, 1]. We model this requirement by a joint chance constraint in order to guarantee network stability. This means that the desired compliance probabilities for all power flows and generator outputs are simultaneously met. Thus, combining all modeling elements considered in the previous sections, we formulate the joint chance constrained DC optimal power flow problem with discrete curtailment as follows: Robust DC optimal power flow with modeling of solar power supply… where the wait-and-see variables (3) and (5), respectively.

Safe approximation of the chance constraints
Chance constrained optimization problems like (6) are in general hard to solve and may not be algorithmically tractable. Therefore, a large number of approximation techniques can be found in the literature, see Prékopa (1995) for a broad overview of the paradigm of stochastic optimization.
Thus, following Nemirovski (2012), we will replace the chance constraint considered in (6g) by a strictly robust protection against a suitably chosen uncertainty set B ∈ B(R n ) that fulfills where B(R n ) denotes the -algebra of Borel sets in the n-dimensional Euclidean space R n .
The robust approximation of (6) is then given by where P G,u k , p u kl are determined as in (3) and (5) replacing X( ) by u. One can show that every feasible solution of the safe approximation (8) is feasible for (6), see Gorissen et al. (2015). To ensure that the safe approximation generates not overly conservative solutions, the uncertainty set B should be chosen as small as possible, but as large as necessary.
Assuming that has been shown in Aigner et al. (2021) that the optimization problem (8) possesses an equivalent mixed-integer linear reformulation which -although being NP-hard in general -can be solved e.g. with the Gurobi optimizer [23] within reasonable time also for huge instances.

Modeling the distribution of the random forecasting error
In order to solve the safe approximation (8) of the stochastic optimization problem (6) described in Sect. 2.3, a suitable uncertainty set B ⊂ R n has to be determined such that (7) holds. For the novel construction of uncertainty sets with the help of copulas, we propose a method for modeling the multivariate probability distribution of the n-dimensional power forecasting error X = P PV − P F introduced in (2). The model for the distribution of X is based on R-vine copulas, which are fitted to empirical data.
To make the paper self-contained, we first give a brief overview of some fundamentals of copula theory in Sect. 3.1. In Sects. 3.2 and 3.3 we explain how R-vine copulas are structured and how they can be fitted to empirical data. Once an R-vine copula is fitted for the distribution of the random fluctuation vector X, in Sect. 3.4 we explain how samples can be drawn from it, in order to determine an uncertainty set B ⊂ R n of the form given in (9) which satisfies a slightly modified version of condition (7), see Sect. 3.6. Furthermore, in Sect. 3.5 we propose a modification of the fitting procedure for R-vine copulas in order to fit the distribution of the (2n)-dimensional random Robust DC optimal power flow with modeling of solar power supply… vector (S, X) to empirical data, where S ∶ Ω → [0, ∞) n models the forecasted solar radiation at the n nodes of the electrical network. This allows for an enhanced modeling of uncertainty sets B s ∈ B(R n ) conditioned on S = s for any given radiation forecast s ∈ [0, ∞) n .

Copulas: definition and sklar's representation formula
the mutual interdependence of the components U 1 and U 2 can be described. For example, the product copula, where models the case that U 1 and U 2 are independent random variables. On the other hand, if C(u 1 , u 2 ) = min{u 1 , u 2 } for all u 1 , u 2 ∈ [0, 1] , then ℙ(U 1 = U 2 ) = 1 , i.e., the components U 1 and U 2 are identical almost surely. Besides these two extreme cases, many further (parametric) families of bivariate copulas C ∶ [0, 1] 2 → [0, 1] can be found in the literature, which model the case that U 1 and U 2 are neither independent nor identical. In particular, for the purposes of the present paper, the following bivariate copula families will be considered: Gaussian, Student t, Clayton, Gumbel, Frank, Joe, BB1, BB6, BB7, BB8 and their rotations, see e.g. Joe (2015); Nelsen (2006) for details.
Note that the notion of a copula is not restricted to the bivariate case. For any integer m ≥ 2 , the function C ∶ [0, 1] m → [0, 1] is called a copula if it is the CDF of an m-dimensional random vector U = (U 1 , … , U m ) ∶ Ω → [0, 1] m such that the (marginal) distributions of U 1 , … , U m are the standard uniform distribution on the unit interval [0, 1]. The importance of copulas results from Sklar's representation formula, see Joe (2015); Nelsen (2006), which states that the CDF of any random vector Y = (Y 1 , … , Y m ) ∶ Ω → R m with arbitrary (not necessarily uniform) marginal distributions can be written as the superposition of the univariate CDFs of Y 1 , … , Y m and a certain copula Vice versa, for any sequence F 1 , … , F m of univariate CDFs and for any copula C, the superposition of F 1 , … , F m and C considered on the right-hand side of (11) is the CDF of an m-dimensional random vector.
for all y 1 , … , y m ∈ R,

R-vine copulas
Note that the representation formula given in (11) can not directly be used in order to fit multivariate probability distributions to data. For this, sufficiently simple and, simultaneously, flexible parametric families of multivariate copulas C ∶ [0, 1] m → [0, 1] are needed. One possible way to construct such parametric copula families is given by so-called R-vine copulas (regular vines), which is a generalization of D-vine copulas recently applied, e.g. in Schinke-Nendza et al. (2021); von Loeper et al. (2021), to model data from meteorology and solar power supply. The structure of R-vine copulas offers the advantage that the probability distribution of the m-dimensional random vector Y = (Y 1 , … , Y m ) to be modelled can be expressed in terms of a number of bivariate copulas. Hereby the structure of an R-vine copula is given by a vector of trees R = (T 1 , … , T m−1 ) with the following properties, see also Fig. 1: only if these edges share one common vertex. Let E(R) denote the set of all edges in R , meaning that E(R) = E 1 ∪ … ∪ E m−1 . Furthermore, we need the following notation. First, for each e = {v 1 , v 2 } ∈ E 1 we define S(e) = � and O(e) = {v 1 , v 2 } . Next, we iterate over i ∈ {2, … , m − 1} and, . We call S(e) the conditioning set and O(e) the conditioned set of edge e. According to Kurowicka and Joe (2010), it holds that |O(e)| = 2 for each e ∈ E(R) and, for each pair of indices be a vector of trees with the properties mentioned above. Then, the following representation formula is true, see Czado (2019); Bedford and Cooke (2001); Joe (2015): For any where Y S(e) denotes the random vector consisting of those components of Y = (Y 1 , … , Y m ) the indices of which belong to the set S(e) ⊂ {1, … , m} , and, analogously, y S(e) is the corresponding subvector of (y 1 , … , y m ) . Furthermore, c o 1 ,o 2 |Y S(e) =y S(e) ∶ R 2 → [0, ∞) denotes the bivariate copula density of the conditional probability distribution of the two-dimensional random vector Note that the right-hand side of (12) is the product of uni-and bivariate functions. Thus, in order to determine the multivariate probability density f 1,…,m , we just have to determine the univariate (marginal) densities f 1 , … , f m , the (conditional) univariate CDFs F o j |Y S(e) =y S(e) , and the (conditional) bivariate copula densities c o 1 ,o 2 |Y S(e) =y S(e) for all e = (o 1 , o 2 ) ∈ E(R) , where the recursion formulas (see Aas et al. (2009)) and are used in order to determine the univariate CDFs F o j |Y S(e) =y S(e) for j = 1, 2. (12)

Fitting R-vine copulas to empirical data
In this section we outline how the representation formula given in (12) can be utilized in order to fit an m-dimensional probability density f 1,…,m to empirical data, i.e., for a given sample of k realizations y (1) = (y (1) 1 , … , y (1) m ), … , y (k) = (y (k) 1 , … , y (k) m ) ∈ R m of the random vector Y = (Y 1 , … , Y m ) , where we use the sequential algorithm proposed in Dissmann et al. (2013). First, for each i ∈ {1, … , m} , we use the sample Then, in the next step, a valid tree T 1 = (V 1 , E 1 ) with V 1 = {1, … , m} is chosen such that the expression is maximized with respect to E 1 , where ̂ denotes an empirical version of Kendall's tau, which is defined for pairs of realizations {(x 1 , y 1 ) … , (x n , y n )} of two random variables X and Y where x = (x 1 , … , x n ) and y = (y 1 , … , y n ).
In other words, the edge set E 1 is chosen such that the sum of pairwise empirical correlations between Y o 1 and Y o 2 is maximized, where the sum extends over all edges e = (o 1 , o 2 ) ∈ E 1 . Subsequently, for each e = (o 1 , o 2 ) ∈ E 1 , a bivariate copula C e is fitted. For this, the independence of Y o 1 and Y o 2 is checked via a statistical test Dissmann et al. (2013). If the null hypothesis (stating that Y o 1 and Y o 2 are independent) is not rejected, then the product copula given in (10) is chosen for C e . Otherwise, an (unconditional) bivariate copula Ĉ e and its parameters are fitted to the data vectors )) with the help of a maximum likelihood method Joe (2015). Now, analogously to (15), a valid tree T 2 = (V 2 , E 2 ) with V 2 = E 1 is selected such that the following expression is maximized: Note that |S(e)| = 1 for all e ∈ E 2 . Thus, using (13) and (14), the conditional CDFs (unconditional) bivariate copula Ĉ o 1 ,o 2 and the (unconditional) CDFs F o 1 and F o 2 , which are determined as described above. Then, for each e ∈ E 2 and o 1 , o 2 ∈ O(e) , a bivariate copula Ĉ o 1 ,o 2 |S(e) and its parameters are fitted to the data vectors (F o j |Y S(e) =y (1) Finally, in the same way as described above, the trees for j = 1, 2 and = 1, … , k , and the bivariate copulas

Sampling from multivariate probability densities
In Sect. 3.3 we showed how the multivariate probability density f ∶ R m → [0, ∞) given by the representation formula for (y 1 , … , y m ) ∈ R m can be fitted to empirical data. We now explain how samples can be drawn from the probability density given in (17).
Recall that the Rosenblatt transform Joe (2015) maps a sample y = (y 1 , … , y m ) of a random vector Y = (Y 1 , … , Y m ) with joint probability density f 1,…,m ∶ R m → (0, ∞) onto a sample u = (u 1 , … , u m ) of a vector of independent and uniformly distributed random variables Assuming that the densities f Y i |Y 1 =y 1 ,…,Y i−1 =y i−1 for i = 1, … , m − 1 are positive, the CDFs F Y i |Y 1 =y 1 ,…,Y i−1 =y i−1 are bijective for i = 1, … , m − 1 and thus, by applying the inverse CDFs to both sides of the above equations, we obtain the inverse Rosenblatt transform: which maps a sample u = (u 1 , … , u m ) of U onto a sample y = (y 1 , … , y m ) of Y. Note that the (inverse) Rosenblatt transform works for any permutation of the indices 1, … , m. Now, consider some sequence of edges e (1) , … , e (m−1) with e (i) ∈ E i for i = 1, … , m − 1 such that e (i) ∈ e (i+1) for i = 1, … , m − 2 . For the given edges, it follows from the third property of the trees T 1 , … , T m−1 introduced in Sect. 3.2 that there is a permutation (o 1 , … , o m ) of (1, … , m) such that o 1 ∈ O(e (1) ) and o i+1 ∈ O(e (i) ) for i = 1, … , m − 1 . Thus, the inverse Rosenblatt transform can be used as follows, in order to draw a sample (y 1 , … , y m ) from the probability density f 1,…,m given in (17): where u = (u 1 , … , u m ) is a sample of a vector of independent and uniformly distributed random variables U = (U 1 , … , U m ) ∶ Ω → [0, 1] m , the (unconditional) CDF F o 1 is given by an integrated kernel density estimator (KDE), and the (conditional) CDFs F o i |Y S(e (i−1) )∪{o i−1 } for i = 2, … , m are determined as described in Sect. 3.3. Later on, in Sect. 4, the algorithms stated in Sects. 3.3 and 3.4 are applied to derive the numerical results presented in this paper, where the implementation provided by the python library pyvinecopulibNagler and Vatter (2021) is used.

Conditional sampling
In the previous section we described a method how to sample from a multivariate distribution with the help of the Rosenblatt transform. This method is used in Sect. 4 below in order to draw samples from the (unconditional) distribution of the forecasting error X = P PV − P F . Furthermore, to model the distribution of the random fluctuation vector X more accurately, we modify the approach considered in Sects. 3.3 and 3.4 such that we can draw samples from the conditional distribution of X for any given radiation forecast S = s . For D-vine copulas, a similar conditional sampling algorithm can be found in Aas et al. (2021) and Bevacqua et al. (2017). Let m, m ′ ≥ 1 be some integers with m ′ < m . We first explain the reasons why the fitting and (unconditional) sampling approach considered in Sects. 3.3 and 3.4 has to be modified such that we can draw samples from arbitrary conditional distributions of a random vector Y = (Y 1 , … , Y m ) , i.e., to draw samples y = (y 1 , … , y m ) from the conditional distribution of Y = (Y 1 , … , Y m ) , given that Y i 1 = y i 1 , … , Y i m � = y i m � for some subset of indices D = {i 1 , … , i m � } ⊂ {1, … , m} and some vector (y i 1 , … , y i m � ) ∈ R m � , Recall that the (direct and inverse) Rosenblatt transform considered in Sect. 3.4 works for arbitrary permutations of the sampling order provided that all conditional CDFs required for this transformation are known. Here, the sampling order refers to the order of the marginal dimensions from which samples are drawn. However, if we want to obtain these CDFs with the help of (13) and (14), the structure of the underlying R-vine copula restricts the choice of possible sampling orders. To understand why this is the case, note that in order to sample in any given order would require the construction of arbitrary (conditional) CDFs, the total number of which is equal to m2 m−1 . However, an R-vine copula of dimension m consists of m(m−1) 2 bivariate copulas. With the help of (13) and (14) two (conditional) CDFs can be obtained from each bivariate copula, i.e., we can obtain m(m − 1) (conditional) CDFs in total from a given R-vine copula, which limits the number of possible sampling orders.
Consider the R-vine copula in Fig. 1 which has (1, 2, 3, 5, 4) as a possible sampling order. To sample in this order with the inverse Rosenblatt transform, we obtain the required inverse CDFs F −1 and F −1 Y 4 |Y 1 =y 1 ,Y 2 =y 2 ,Y 3 =y 3 ,Y 5 =y 5 from the marginal distribution 1 and the copulas 1, 2 , 1, 3 | 2 , 1, 5 | 2, 3 and 1, 4 | 2, 3, 5 respectively. Note that this sampling order is possible because each copula corresponds to an edge connected to the previous copula or marginal distribution, e.g., 1, 3 | 2 corresponds to an edge connected to 1, 2 while 1, 2 corresponds to the edge connected to 1 . This ensures that a suitable copula for the next dimension in the sampling order exists. Now consider the sampling order (1, 2, 3, 4, 5), which is impossible. Analogously to the previous sampling order the inverse CDFs F −1 can be obtained. However, to obtain the 4 th necessary inverse CDF F −1 Y 4 |Y 1 =y 1 ,Y 2 =y 2 ,Y 3 =y 3 for the inverse Rosenblatt transform, the copulas 1, 4 | 2, 3 or 3, 4 | 1, 2 are required which do not exist within the considered R-vine copula.
As shown in Theorem 5.1 in Cooke et al. (2015), an R-vine copula of dimension m has only 2 m−1 possible sampling orders. This is due to the fact that every possible sampling order corresponds to a vector = ( 1 , … , m ) = (v 1 , … , v m−1 , e) with v i ∈ T i and e ∈ E m−1 , i.e., the CDF used in the first equation of the Rosenblatt transform is the marginal CDF F v 1 whereas the conditional CDFs of the equations thereafter are given by the copulas corresponding to 2 , … , m . Recall that for each equation of the Rosenblatt transformation the dimension of the condition of the corresponding conditional CDF grows by one. This restricts the choice of i+1 to copulas for which it holds that i ∈ i+1 for all i ∈ {1, … , m − 2} , i.e., i ∈ T i must be a vertex of the edge i+1 ∈ T i+1 , because only then the copula corresponding to i+1 can be used to construct a conditional CDF with a valid condition for the (i + 1) -th equation of the Rosenblatt transform.
We thus modify the fitting process for vine copulas presented in Sect. 3.3 such that a vector = ( 1 , … , m ) as described above exists for a given set of indices D = {i 1 , … , i m � } ⊂ {1, … , m} . For this, we consider , i.e., T D 1 is a graph with vertex set D and edges e ∈ E 1 which connect two vertices in D. For i > 1 , we recursively define . Note that in general T D i is not a tree but a forest, however, only if all T D i are trees the vector (T D 1 , … , T D m � ) is a valid R-vine copula. This is necessary to construct an inverse Rosenblatt transform for the dimensions in D, or more generally speaking, it is necessary for the construction of an inverse Rosenblatt transform for all dimensions {1, … , m} where the dimensions in D occur at the beginning.
To ensure that there is a sampling order in which all indices in D are in successive order, we choose the graphs T D i in the fitting process of the R-vine copula such that I(E D i ) in (15) is maximized (as in the unmodified fitting process considered in Sect. 3.3), where additionally it must hold that T D i is a tree for all i ∈ {1, … , m � } because only then we can chose a sampling order where i ∈ i+1 holds for all i ∈ {1, … , m − 1}.
Without loss of generality, we now assume that D = {1, … , m � } . Thus, we omit the first m ′ equations of the inverse Rosenblatt transform and sample values y m � +1 , … , y m for the remaining m − m � components via As an example, consider again the R-vine copula in Fig. 1 and the set D = {1, 2, 3} to sample from the conditional distribution of (Y 4 , Y 5 ) | Y 1 =y 1 ,Y 2 =y 2 ,Y 3 =y 3 . The graphs T D 1 , T D 2 and T D 3 with the sets of vertices { 1 , 2 , 3 } , { 1, 2 , 2, 3 } , { 1, 3 | 2 } and the corresponding edges correspond to the lower left part of the diagram. Since the graphs T D 1 , T D 2 and T D 3 are trees and (T D 1 , T D 2 , T D 3 ) is a valid R-vine copula, sampling orders with 1, 2 and 3 at the beginning are possible. Now consider D = (1, 2, 4) for which T D 1 = ({ 1 , 2 , 4 }, {{ 1 , 2 }}) is not a tree and the vector ( is not an R-vine copula. Since 4 is not connected to 1 or 2 in T D 1 there can be neither 1, 4 nor 2, 4 in T D 2 and in turn there can be neither 2, 4 | 1 nor 1, 4 | 2 in T D 3 . Therefore it is not possible to obtain the required inverse CDFs for an inverse Rosenblatt transform for which the sampling order begins with the elements of D.
In the following section, we explain how the construction of uncertainty sets is performed with the scenario approach from stochastic optimization. We then use the copula-based modeling from this section in order to construct high-quality uncertainty sets for given weather situations.

Scenario approach to determine a suitable uncertainty set
In order to determine a suitable uncertainty set of the form given in (9) which satisfies (a slightly modified version of) condition (7), we apply, as in Aigner et al. (2021), an idea described in Margellos et al. (2014) and formulate the estimation of the uncertainty set B = [ 1 , u 1 ] × … × [ n , u n ] ⊂ R n as an auxiliary probabilistic optimization problem. Then, for this problem with chance constraints, we apply the scenario approach proposed in Campi and Garatti (2008), i.e., the chance constraints considered in (7) are replaced by constraints based on a sufficiently large number of samples drawn from the probability distribution of the random forecasting error X = P PV − P F . In this work this distribution is fitted to empirical data, using the algorithm described in Sect. 3.3, and simulation is performed with the technique described in Sect. 3.4. The auxiliary optimization problem in its general form consists of a chance constraint model for the enclosure B ∈ B(R n ) of the probability mass of X = (X 1 , … , X n ) satisfying the condition ℙ({ ∶ X( ) ∈ B}) ≥ 1 − for some ∈ (0, 1) , see (7). At the same time, this problem aims for an uncertainty set B such that its size is as small as possible. Thus, in order to apply the scenario approach proposed in Campi and Garatti (2008) to determine an uncertainty box B = [ 1 , u 1 ] × … × [ n , u n ] ⊂ R n , we consider the probabilistic optimization problem where the minimum in (18a) extends over all = ( 1 , … , n ), u = (u 1 , … , u n ) ∈ R n with k < u k for all k = 1, … , n.
Thus, to control the size of the set B, we minimize the sum of interval lengths u k − k . In contrast, if minimization of the box volume were used instead, this would lead to a non-convex objective. In this case, the scenario approach proposed in Campi and Garatti (2008) is no longer applicable. Although the solution of (18) does not necessarily minimize the box volume, the solution of the following scenario program does. This is why this choice of objective is suitable. We further explain this after introducing our scenario program.
Suppose that N > 0 samples x 1 , … , x N are independently drawn from the probability distribution of X. Instead of (18b), in our scenario approach we want to ensure that the samples x 1 , … , x N are included in the uncertainty set B. The resulting scenario program for computing B = [ , u] is thus given by

3
The solution of this optimization problem can be written explicitly as It is true that set B * = [ * , u * ] also minimizes the volume over all sets [l, u] containing the samples x 1 , … , x N . Although, in general, the solution of problem (18) does not calculate boxes with minimal volume, this is the case for the optimization problem given in (19).
From the results presented in Campi and Garatti (2008), we know that the optimal solution B * = [ * , u * ] of (19) fulfills condition (18b) with a confidence probability of at least 1 − for some small ∈ (0, 1) if N > 0 is chosen such that Note that in the latter inequality, the necessary number of samples N > 0 for a predefined confidence level 1 − ∈ (0, 1) is given implicitly. However, an explicit sufficient condition has been derived in Alamo et al. (2010), which reads as Furthermore, we determine the optimal solution B * s = [ * s , u * s ] of (19) based on samples drawn, as described in Sect. 3.5, from the conditional distribution of X for given radiation forecasts S = s.

Numerical results
In order to derive the results presented in this section we used the library pyvinecopulib Nagler and Vatter (2021). Furthermore, we utilized Gurobi 9.1.2 [23] as solver for mixed-integer linear programs. The computations were carried out by means of a python implementation on a cluster using 4 cores of a machine with

Fig. 2
Sketch of NNG subnetwork, where (slack-)generator nodes are denoted by + , solar feed-in points by ⋆ and load buses by • 1 3 Robust DC optimal power flow with modeling of solar power supply… two Xeon E3-1240 v6 "Kaby Lake" chips (4 cores, HT disabled) running at 3.7 GHz with 32 GB of RAM.

Data description
Data regarding power measurements as well as weather forecasts were provided by the distribution network operator N-ERGIE Netz GmbH (NNG) and the German weather service Deutscher Wetterdienst (DWD). In particular, NNG provided data of solar power supply at more than 150 feed-in points and corresponding active power measurements at 13 network nodes (buses) measured in 15 min intervals. Moreover, NNG provided data regarding the positions of network nodes (buses) and their connections through lines (branches) which include resistance values and transmission limits of each line in the distribution network. A fragment of the NNG distribution network with 34 nodes and 37 lines is visualized in Fig. 2. The solar power forecast P F is provided by a model proposed in Schinke-Nendza et al. (2021). DWD provided hourly forecasts of global horizontal irradiation, which were generated by the ensemble system of the numerical weather prediction model COSMO-DE, called COSMO-DE-EPS, and statistically interpreted based on synoptic observations at weather stations by Ensemble-MOS of DWD, see Hess (2020). The weather forecasts are issued on a 20 km × 20 km grid covering Germany and parts of the neighboring countries at every third hour. The forecasts of global horizontal irradiation were We split the data into a training set and a validation set. The training set is used to fit model parameters and consists of data from the years 2015 and 2016. Based on the validation set from 2017 the accuracy of the predictions generated by the fitted model is evaluated.

Fitting unconditional and conditional distributions of forecasting errors
In this section we discuss the fitting of R-vine copulas, as outlined in Sects. 3.3 and 3.5, in order to determine uncertainty sets B * of the form introduced in Sect. 3.6. First we explain how to model the (unconditional) distribution of the n-dimensional random vector X = P PV − P F of power forecasting errors at the n nodes of the electricity network considered in the present paper, where n = 13 . Besides this, we additionally consider the random vector S = (S 1 , … , S n ) ∶ Ω → [0, ∞) n , which describes the forecasted solar radiation at the n nodes of the electricity network, and we model the conditional distribution of X given that S = s for some s ∈ [0, ∞) n . Moreover, we consider two further types of conditional distributions of X under the condition that S = s and S k = s k , respectively, for some s ≥ 0 , s k ≥ 0 and k ∈ {1, … , n} , where As outlined in Sect. 3, copula theory allows for the modeling of the multivariate distribution of random vectors like the random power forecasting error X ∶ Ω → R n . In order to estimate the univariate (marginal) CDFs F X 1 , … , F X n we use numerically integrated KDEs, with a Gaussian kernel and a bandwidth being equal to the estimated standard deviations k of X k for k = 1, … , n , see the left column of Fig. 3. Once an R-vine copula is fitted to the distribution of X, as descibed in Sect. 3.3, we are able to draw realizations from the fitted distribution of X, with which the uncertainty set B * can be determined as described in Sect. 3.6. This method results in one single uncertainty set B * for all considered hours, since the fitted R-vine copula models the (unconditional) distribution of X, irrespective of other variables, which are possibly correlated with X. Thus, it is sensible to investigate if and to which extent the random vector X of power forecasting errors depends on various other variables, like the random vector S of forecasted solar radiations at the n nodes. For this reason, we also model various conditional distributions of X.
To condition on the forecasted solar radiation vector S, we consider the three cases mentioned above, i.e., S = s , S = s , and S k = s k for some k ∈ {1, … , n} . From a meteorological perspective, the network nodes in N are in close geographical proximity and, therefore, the forecasted solar radiations S 1 , … , S n at the n network nodes are highly correlated. Thus, it might be sufficient to consider either the average solar radiation S or the solar radiation S k for one single node, instead of the random vector S, which reduces the complexity of the copula model without much loss of information.
As can be seen in Fig. 3, the power forecasting errors X k , X k ′ have unimodal distributions which are well approximated by KDEs. For the forecasted solar radiations, S k , S k ′ , however, the values of the densities are significantly larger than zero at the distribution limits. Since the kernel of the KDE would cross the bounds of the distribution for data points close to those bounds, we first transform the components of S, as well as S and S k , using the mapping T ∶ where F N(0,1) is the CDF of the standard normal distribution and F U(a,b) is the CDF of U(a, b), the uniform distribution for the interval [a, b] for some a, b ∈ R with a < b . Thus, T maps the bounded interval [a, b] onto R . Since the endpoints a and b are mapped to −∞ and ∞ , respectively, we choose them to be slightly outside the bounds of the solar radiation distribution such that T does not map any data point to ±∞ . The ranges of values of the transformed random variables T(S), T(S) and T(S k ) are unbounded and we can apply kernel density estimators to their transformed data points T(s), T(s) and T(s k ) , where T(s) = (T(s 1 ), … , T(s n )) . Finally, we transform the density functions f Once the densities of the marginal distributions of X and S, as well as the densities of S and S k are determined, they are numerically integrated to obtain the corresponding CDFs with which an R-vine copula is fitted, as described in Sect. 3.3. Now we can draw samples from the (unconditional and conditional) R-vine copula model with which we construct uncertainty sets B * , as described in Sect. 3.6. Figure 4 shows the histograms of samples drawn from conditional R-vine copula models for different solar radiation forecasts and, in particular, how the conditional error distribution changes for different forecasted solar radiations.
To check how well the R-vine copula model captures the correlations of the dataset of forecasted radiations and power forecasting errors, we compare the values of empirical Kendall's tau (see (16)) for all pairs of components of the vector (S 1 , … , S n , X 1 , … , X n ) . It can be seen in Fig. 5 that the R-vine copula model manages to capture the correlation within the underlying dataset quite well, since the values of empirical Kendall's tau computed from the dataset of forecasted radiations and power forecasting errors (left) and from simulated realizations of the R-vine copula model (right), respectively, show very similar correlation structures.
Note that we consider copulas with up to 26 dimensions while the available dataset contains only 180 data points. This makes it difficult to reliably assess the goodness of fit of the copula model. However, in the following we evaluate the entire model chain with various validation scores in order to assess the additional benefit of the copula model.

Analyzing the size of uncertainty sets
We now analyze the size of uncertainty sets for the robust approximation of chance constraints using the scenario approach described in Sect. 3.6. The resulting sets depend on the samples drawn from the unconditional probability distribution and the three conditional distributions of power forecasting errors, respectively, considered in Sect. 4.2. Note that the minimum number N of samples required for the scenario approach, determined by means of (20), goes  For the numerical results discussed in the present section, we use an average uncertainty set which is obtained from applying the scenario approach 500 times. In this way, our numerical results become reproducible because the average uncertainty set does not change significantly, when the procedure described above is repeated. Figure 6 shows values of the size measure given in (18a), i.e. for the sum of interval lengths, of uncertainty sets computed exemplarily for a usual summer day at noon with an average hourly global horizontal irradiation of 0.63 kWh m 2 , in dependence of different values of the coverage probability 1 − with a confidence of 1 − = 0.99 . Note that smaller confidence levels would lead to smaller uncertainty sets, but the quality of these sets also decreases. In particular, there would no longer be a confidence probability of 0.99 that the computed uncertainty set covers the chosen probability mass of 1 − .
The values displayed in Fig. 6 are normalized by the size of the largest uncertainty set, namely the unconditional uncertainty set for a coverage probability of 0.99. It can be seen that the sizes of the uncertainty sets increase with increasing probabilities 1 − as the confidence regions cover a larger set of realizations of the random vector X of power forecasting errors. In comparison to the uncertainty sets constructed with conditional probability distributions of X, the unconditional distribution of X leads for all coverage probabilities 1 − to larger uncertainty sets. Thus, with knowledge on the forecasted solar radiation, it is possible to adapt the uncertainty sets to the current weather situation, which leads to small sizes. Not surprisingly, the conditional distribution of X with given solar radiation at all n solar feed-in nodes yields the smallest uncertainty sets for all coverage probabilities 1 − . However, the differences between these sizes and those obtained for the other two conditional settings with less complete information on the forecasted solar radiation, i.e. knowledge of average solar radiation (avg), and at one single node (one), are not too large. Furthermore, the size differences between the conditional settings 'avg' and 'one' are negligible. The numerical results presented in the remaining part of this section concern the case 1 − = 0.9 , i.e. the practically most relevant value of the coverage probability 1 − . For this safety margin, we analyze the uncertainty sets obtained for the four (unconditional and conditional) distributions of X described above and for each day in the validation dataset. In particular, we determine the empirical coverage probability by counting how often the realizations drawn from the respective distribution of the random vector X belong to the corresponding uncertainty set. Furthermore, we compute and compare the average size of the uncertainty sets, i.e. the sum of interval lengths, and their average volume, i.e. the product of interval lengths. The results are displayed in Table 2, where it can be seen that the four different settings lead to similar empirical coverage probabilities around the given level of 0.9. On the other hand, the reductions of size and volume of uncertainty sets implied by considering conditional distributions of the power forecasting error X are clearly visible. Again, the case with given solar radiation at all n solar feed-in nodes yields the smallest uncertainty sets, whereas the size differences between the conditional settings 'avg' and 'one' are negligible.
To further analyze the impact of additional knowledge regarding solar radiation forecast on size and location of uncertainty sets, we determined uncertainty sets for a rather sunny day at noon with a high average solar radiation forecast of 0.76 kWh m 2 and a less sunny day at noon with a low average solar radiation forecast of 0.18 kWh m 2 . The results are shown in Fig. 7, where the uncertainty sets are plotted via their confidence intervals (in MW) for each solar feed-in point.
It turned out that the lengths of the confidence intervals significantly shrink by considering conditional distributions of the power forecasting error X, given a high average solar radiation forecast. More precisely, the lower endpoints of the confidence intervals are shifted upwards, i.e., negative power forecasting errors are less likely, whereas the upper endpoints remain almost unchanged, see Fig. 7 (left). On the other hand, for low average radiation forecast, the confidence intervals are shifted downwards by considering conditional distributions of the power forecasting error, but their lengths remain almost unchanged, see Fig. 7 (right).
Finally, we note that also the results of the numerical experiments presented in Aigner et al. (2021) are based on (measured) power feed-in data from NNG and forecasted radiation data from DWD. However, the used database differs from that of the present paper, where, in addition, solar power forecast data are exploited provided by the forecasting model of Schinke-Nendza et al. (2021). In this way, by modeling the multivariate probability distribution of solar power Robust DC optimal power flow with modeling of solar power supply… forecast data via R-vine copulas, it is possible to determine conditional uncertainty sets, which meet the desired coverage probability of 0.9. They have significantly smaller sizes than the corresponding unconditional uncertainty sets from Aigner et al. (2021) which led to an larger empirical coverage of 0.98 although 1 − = 0.9 was required.

Robust curtailment
As important as the size of the computed uncertainty sets is the quality of solutions obtained by solving the robust approximation (8) of the chance constrained optimization problem described in (6). In order to solve (8), we use the network parameters given by the power network operator NNG. The curtailment options for the feedin nodes in the electrical power network of NNG are k ∈ {0, 0.1, 0.2, … , 1.0} . Moreover, the participation factors of the generators are fixed values given by NNG ( 31 = 34 = 0.05 , 32 = 33 = 0.45 ). There are no costs affiliated with the power transfer at the (slack-) generators on the boundary nodes. Hence, there are no generator production costs and the corresponding term in the objective function is given as The minimization of this objective function leads to a minimum curtailment of solar feed-in.
Due to the balanced network situations in the historical data, there is no need to curtail the solar feed-in in the instances from the validation set. There is also no danger of overload and the optimization leads to trivial solutions with a curtailed solar power equal to 0. Thus, in order to generate test cases with critical network situations (and non-trivial solutions), we artificially increased the solar power feed-in, whereas the network topology, transmission line parameters and the power demand remained unchanged. More precisely, based on the data of the validation set, we increased the installed solar power and the feed-in up to the by NNG planned total solar power capacities of the year 2022 and the planned total solar power increase of year 2025. The corresponding scaling of power generation forecast and uncertainty sets creates an oversupply of renewable energy, and therefore it is more likely in these instances that a curtailment will be required. Furthermore, in addition to the up-scaled solar power, we simulated the impact of transmission line failure on the solution of our optimization problem.
Thus, we now discuss further details for the following experimental setups: A: Installed solar power as planned in 2025, B: Installed solar power as planned in 2022 with a failure of lines (6, 19) and (9, 30).
To obtain the results, a mixed-integer optimization problem was solved for each instance and each (unconditional and conditional) uncertainty set. The computing times are very low and, thus, solutions can be generated efficiently. Indeed, the average computing times for the two settings are 2.8s (setting A) and 1.1s (setting B), with a maximal run time of 8.2s (setting A) and 4.0s (setting B).
The robustness of a solution of (8) can be validated by checking if the computed network configuration leads to an overload after the realization of uncertainty. The corresponding entries in Table 3 show that nominal solutions generated without probabilistic constraints (or, in other words, for 1 − = 0 ) lead to overload in a large amount of test instances. In contrast, only up to three robust solutions lead to constraint violation in each setting for the different probabilistic models. The relative frequencies for this is therefore below the given threshold of = 0.1 . This indicates the feasibility of the robust solutions for the chance constraints. This shows that the robust protection against uncertainties is necessary and reasonable, since the number of technical constraint violation could be strongly reduced in the numerical experiments.
To further investigate the quality of solutions of (8), we computed the amount of curtailed solar power of the robust solution in comparison to the solution of the nominal problem (1) without a protection against uncertainty. The increase in curtailed energy of the robust solutions in comparison to the nominal ones can be interpreted as the cost of robust protection. That means how much the curtailment costs increase due to the protection against uncertainties. Figure 8 shows box plots for the increase of relative curtailment costs using the four (unconditional/conditional) types 'no', 'avg', 'one' and 'all' of uncertainty sets. One can see that, again, the Table 3 Number of instances where a nominal solution of (1) without probabilistic constraints leads to overload in the network compared to robust solutions (with security of 1 − = 0.9 ) using the four (unconditional/conditional) types 'no', 'avg', 'one' and 'all' of uncertainty sets addition of further knowledge about the solar radiation improves the performance in both settings. This corresponds to the size reduction of the uncertainty sets recognized in Sect. 4.3. Overall, the relative cost increase in all experiments is relatively small. However, using the samples drawn from the three conditional distributions of power forecasting errors enable us to further reduce the amount of wasted energy under the same solution guarantees, where, again, the conditional settings 'avg' and 'one' have a similar impact. In comparison with the preliminary results obtained in Aigner et al. (2021), the amount of curtailed energy could drastically reduced on average from about 13 to 5% under the same solution quality guarantees. This coincides with the reduction of uncertainty set size discussed at the end of Sect. 3.6. In summary, the obtained results show that the scenario approach for the considered instances in combination with the copula-based stochastic modeling of power forecasting errors leads to high-quality solutions. The addition of further knowledge about the current weather situation allows us to construct more precise uncertainty sets. We are able to produce robust solutions with a relative small increase of curtailment costs, while maintaining the same level of protection.

Conclusion
In this paper, we combine the robust approximation of chance constrained DC Optimal Power Flow with a probabilistic uncertainty model based on R-vine copulas to reduce the curtailment of solar power while keeping the power grid stable. The chance constrained DC Optimal Power Flow determines appropriate levels of curtailment based on a deterministic forecast for the expected solar power feed-in and uncertainty sets, i.e., multidimensional cuboids which contain the forecasting error with a given probability. These uncertainty sets are approximated with the help of the multivariate probability distribution of the forecasting error at all considered power grid nodes. This results in less curtailments and a more stable power grid compared to the results of a model without uncertainty sets. To further improve upon these results, we incorporate knowledge about solar radiation in the solution process by considering the conditional forecasting error distribution for a given solar radiation forecast. This leads to sharper distributions, i.e., the forecasting error can be predicted with higher accuracy, which results in smaller uncertainty sets. Compared to the unconditional case, this leads to even less curtailments and improved stability of the power grid.
Our numerical results demonstrate the applicability of our procedure and the positive effects of incorporating a probabilistic model for the distribution of random solar radiation vectors. Future research can transfer our solution framework to different applications under uncertainty like in energy network optimization.
Future research could add further features and investigate questions arising from the application, for example adding optimal transmission switching under uncertainty or including storage elements and unit commitment constraints over time. From a mathematical point of view, it would be interesting to study different geometries for uncertainty sets to further reduce the conservatism of the robust approximation. The major challenge is to find assumptions where an equivalent reformulation for the resulting problems is possible. In order to improve the copula-based sampling from conditional probability distributions, it might be promising to add more information (e.g. temperature, solar altitude, time) to the model.