{'fraction_non_alphanumeric': 0.04350622328150418, 'fraction_numerical': 0.012522226005372073, 'mean_word_length': 4.606362672322375, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Gaussian Mixture Model Based Distributionally Robust Optimal Power Flow With CVaR Constraints

This paper proposes a distributionally robust optimal power flow (OPF) model for transmission grids with wind power generation. The model uses the conditional value-at-risk (CVaR) constraints to control the reserve and branch flow limit violations caused by wind power forecast errors. Meanwhile, the Gaussian mixture model (GMM) is integrated into the CVaR constraints to guard against the non-Gaussian forecast error distributions. Unlike the previous studies considering the GMM with fixed parameters, this paper allows the GMM parameters to be variable within some credible regions and develops a data-driven GMM-based ambiguity set to achieve the distributional robustness. Also, rather than using the traditional sample-based approximation of CVaR with high computational burden, this paper designs a scalable cutting-plane algorithm to handle the distributionally robust CVaR constraints. Case studies on the IEEE 2736-bus system show the effectiveness and scalability of the proposed OPF model.


I. INTRODUCTION
The forecast errors of wind power can pose considerable risk to power system security. To handle the wind power forecast uncertainties, different decision-making approaches have been proposed, such as the stochastic optimization (SO) using discrete wind power scenarios [1] and the robust optimization (RO) leveraging the uncertainty sets of wind power [2].
The distributionally robust optimization (DRO) has been emerging as an intermediate between SO and RO. DRO assumes the true distribution of the forecast error lies in an ambiguity set and makes decisions to hedge against the worstcase distribution. The moment-based ambiguity set has been used for power system dispatch [3], [4]. However, such ambiguity set utilizes only the moment information of the forecast error and may be very conservative. Another type of ambiguity set leverages the statistical distance between distributions, such as the Kullback-Leibler (KL) distance and Wasserstein distance. In [5], a KL-based ambiguity set is used for the unit commitment. In [6], a Wasserstein-based ambiguity set is applied to the optimal power flow (OPF) problem. However, KL-DRO cannot be applied if the uncertainties have heavy tails [7], while the model size of Wasserstein-DRO grows heavily with the data size employed. In [8], [9], a new type of ambiguity set is also constructed based on the cumulative probability distributions of random variables, but the correlations between random variables cannot be captured.
The decision-making under uncertainty also often needs some risk management tools. Value-at-risk (VaR) and conditional value-at-risk (CVaR) are two popular tools and have been successfully applied to power system problems. Based on VaR, the chance constrained OPF and UC models have been widely used to ensure the power system security with specified probability levels [10], [11]. The use of VaR is intuitive since it limits the violation probability of a constraint. However, VaR neglects the aftermath if the violation occurs. Also, the VaR-based chance constraints are often non-linear and non-convex, and thus are computationally intractable.
Unlike VaR, CVaR is a coherent risk measure [12]. It is defined on the mean of the distributional tail exceeding VaR and can thus account for both the probability and magnitude of constraint violation. Besides, the direct optimization of CVaR is a convex problem and is thus more tractable. In power systems, CVaR has been used for demand-side management [13] and unit commitment optimization [14].
Overall, both VaR and CVaR have their rationality and the decision maker may prefer one over the other in different situations. This paper focuses on CVaR. More specific, this paper aims to address two issues in the application of CVaR.
First, the calculation of CVaR may be still not easy since it involves high-dimensional integral. Currently, the samplebased approximation of CVaR is often used [15], but it may cause a huge computational burden under large sample sizes. An analytical reformulation of CVaR is also provided in [15], but it is valid only when the uncertainties follow Gaussian distributions. In reality, the forecast error of wind power often follows non-Gaussian distributions [16]. Recently, Gaussian mixture model (GMM) has been used to fit non-Gaussian forecast errors [17], [18]. However, the CVaR under GMM still lacks an analytical form. It is also noteworthy that though GMM has been used for VaR-based power system scheduling [19], [20], there is no existing work applying GMM to CVaRconstrained power system dispatch problems.
Second, in practice, the precise probability distribution of the forecast error is often unavailable. It could be biased to use the distribution estimated directly from limited historical data or based on expert knowledge. To address this issue, the momentbased ambiguity set with the mean and covariance information has been used to consider the worst-case CVaR [21]. Since the utilization of only moment information may result in very conservative solutions, the unimodality knowledge is further included in [22] to reduce the conservatism. Besides, the Wasserstein-based ambiguity set has also been used for CVaRconstrained OPF [23], but the resulting model is not scalable since the model size is proportional to the number of historical data. Also for the CVaR-constrained OPF, [24] provides an alternative method by combining the sample-based CVaR approximation with an ambiguous discrete probability distribution. However, similar to [23], this method is still not scalable since its size increases heavily with the number of forecast error samples. This paper proposes a new distributionally robust OPF model with CVaR constraints. A data-driven GMM-based ambiguity set is used to hedge against the ambiguous true distribution of wind power forecast error. A scalable solution approach without any sampling procedure is also developed for the proposed OPF model. The contributions of this paper include: 1) This paper extends the application of GMM in power system operational planning. First, previous research often uses GMM in VaR-constrained dispatch problems, while it is the first time for GMM to be applied to CVaR-constrained OPF. Second, the traditional GMM has fixed parameters (e.g., see [19], [20]), while the GMM in this paper is allowed to have uncertain parameters by adopting a data-driven GMM-based ambiguity set. The use of the proposed GMM is more robust to the distributional ambiguity of wind power forecast error.
2) A scalable cutting-plane-based methodology is used to solve the proposed OPF model. Unlike the traditional methods to handle CVaR, the proposed methodology does not rely on any samplings when evaluating the worst-case CVaR and enforcing the distributionally robust (DR)-CVaR constraints and thus avoids the associated high computational burden (see, e.g., [24]). By exploiting the special properties and structure of the proposed GMM-based ambiguity set, a method is designed to quickly evaluate the worst-case CVaR. Then, based on the convexity and sensitivities of CVaR, cutting-planes are used to enforce the DR-CVaR constraints. Finally, an iterative cuttingplane-based solution algorithm is developed for the proposed OPF model. The whole methodology demonstrates high computational efficiency even for very large-scale networks.

A. Optimal Power Flow (OPF) Formulation
As the response to wind power forecast errors, we assume the outputs of traditional units are adjusted based on the policy: where the total power imbalance is allocated among the traditional units according to the participation factors . If ∑ = 1, the system-wide power balance is guaranteed.
With the above policy, the OPF dispatch with uncertain wind power can be modelled as (P1): where (2a) minimizes the sum of the fuel cost and reserveprocuring cost; (2b)-(2d) limit the available up-and down-ward reserves; (2e)-(2f) enforce system-wide power balance; (2g)-(2h) bound the deployment of up-and down-ward reserves; (2i) indicates the nominal branch flow while (2j)-(2k) represent the random branch flow under the wind power fluctuations. In the following analysis, the random constraints (2g)-(2h) and (2j)-(2k) are all represented by the general form: where is the decision vector in (P1), ( ) and ℎ ( ) are two affine transformations of , ̃= (̃ , ∀ ∈ ) , is the total number of the random constraints. In the following analysis, ( ) is re-expressed as to simplify notation.

B. Risk Measure: Conditional Value at Risk (CVaR)
To control the violation risk of (3), we resort to CVaR. Given a risk level ∈ (0, 1), the CVaR of ̃ is defined as [12] CVaR 1− ( ̃) = min In (4), CVaR is the conditional expectation of ̃ exceeding the optimal solution of (i.e., the value-at-risk of ̃ at the probability level 1 − ). Then, we use the CVaR constraints in (5) to control the violation risk of (3).
However, the use of (5) has the following two issues: i) Usually, we only have limited historical data of ̃ and it is hard to estimate the precise distribution of ̃. It can be biased to estimate CVaR by assuming any parametric distribution of ̃.
ii) CVaR lacks a general computational efficient form. Though (5) can be handled by a sample-based approximation, the approximation quality could be sensitive to the number of samples and the sampling technique, and the computational burden could be huge when the number of samples is large.

III. DISTRIBUTIONALLY ROBUST FORMULATION
To address the aforementioned two issues in Section II-B, we propose a distributionally robust approach combining the CVaR constraints (5) with a GMM-based ambiguity set. This is based on the following considerations: i) GMM is adopted due to its flexibility in modelling various types of multivariate non-Gaussian probability distributions.
ii) The GMM fitted on limited training data may deviate from the unknown true distribution of ̃, so a GMM-based ambiguity set is used to account for a GMM with uncertain parameters and hedge against the distributional ambiguity of ̃.
iii) Though the GMM-based CVaR still lacks an analytical form, the use of GMM makes it possible to apply a cuttingplane method to enforce the CVaR constraints (5) (see Section IV). Thus, the sample-based approximation of CVaR is avoided.

A. Gaussian Mixture Model (GMM)
The probability density function (PDF) of the wind power forecast error ̃ is defined as the following GMM: where ( , ) is the multivariate Gaussian distribution with the mean and covariance matrix , and it is the ℎ Gaussian component of the GMM; ( 1 , …, ) are the mixing weights of the Gaussian components; is the vector form of { 1 , …, , 1 , …, , 1 , …, }. Given the historical data of ̃, the parameters and can be obtained by the maximum likelihood estimation (MLE) [25].

B. Constructing Credible Regions for GMM Parameters
Different to the traditional GMM with fixed parameters, we assign some GMM parameters with certain uncertainty and allow them to be variable within some credible regions. This yields an ambiguous GMM. The credible regions are as follows.

C. Quantifying the Uncertainty of GMM Parameters
Next, we provide two optional methods to quantify the uncertainty of , and and determine the coefficients in the credible regions , and Σ . [26] This method estimates the empirical distribution of a statistic by calculating the statistic on many bootstrap-resamples. This empirical distribution is used to construct the confidence region of the statistic. The bootstrap procedure is as follows [26]:

Method 1: A non-parametric bootstrap method
i) With all historical data of ̃, use MLE to obtain the optimal parameters * and * , where * is chosen according to the Bayesian information criterion; then extract the maximum likelihood posterior group membership probability matrix * .
ii) Construct bootstrap resamples { 1 , 2 , …, } by drawing random observations with replacement from the original historical observations of ̃.
With the obtained parameters * , we can determine the coefficient of the credible regions in Section III-B.
For the credible region , we estimate are the upper 1+ 2 and 1− 2 quantiles of a data series, respectively.
For the credible region , we estimate For the credible region Σ , we estimate

Method 2: A method from a Bayesian viewpoint
This method performs asymptotic inference on the estimator of . More specific, under suitable regularity conditions and for large data sizes of ̃, the posterior distribution of can be approximated by the following normal distribution and it does not depend on the prior distribution of (see Section 4 in [27]): where * is the MLE of , ̂( * ) is the estimated Fisher information matrix evaluated at * and it can be calculated based on the methods in [28], is the data size of ̃. With this normal distribution, it is easy to determine the coefficients in the above credible regions and the details are omitted here.
The Bayesian method does not need the bootstrap resamples and is computationally very cheap. As the data size gets larger, the Bayesian method can be more accurate and is expected to offer a lower level of conservativeness associated with the credible regions. This is desirable since it helps make lessconservative decisions when more data information is at hand.

D. Distributionally Robust CVaR-constrained OPF
After determining the credible regions of the GMM parameters, we can construct the following GMM-based ambiguity set to account for an ambiguous GMM distribution: With the ambiguity set in (13), we can consider the following distributionally robust CVaR (DR-CVaR) constraints: where the worst-case CVaR of ̃ is no higher than ℎ ( ).
IV. SOLUTION METHODOLOGY Next, we design an iterative cutting-plane algorithm to enforce the DR-CVaR constraints (14) and solve problem (P2).

A. A Cutting-plane Method
To handle the DR-CVaR constraints (14), we will utilize two important and very desirable properties of CVaR: convex in (see the Corollary 11 in [12]).
ii) it is easy to verify that ̃ satisfies the Assumptions 1-3 in [29], so the gradients of CVaR 1− ( ̃) with respect to can be calculated as the following conditional expectations (see Theorem 3.1 in [29]): where , denotes the ℎ element of , and VaR 1− ( ̃) is the value-at-risk of ̃ at the probability level 1 − .
In [29], a Monte-Carlo sampling method is used to estimate the VaR, CVaR and its gradients in (16). However, this method is computationally inefficient, and it is difficult to account for the distributional ambiguity of ̃ within such method.
Next, we will show that, thanks to the use of GMM, a precise and scalable method can be applied to calculate the VaR, CVaR and gradients in (16) while considering the distributional ambiguity of ̃ as modelled by (13).
Step 1: identify the worst-case and .
For a given * , the corresponding worst-case CVaR of * ̃ can be reformulated as (17) by combining (4), (13) and (14): In (17) where and Φ are respectively the PDF and cumulative distribution function (CDF) of the normal distribution with mean ̅ = * and standard deviation ̅̅ ̅̅ = √ * * . It can be verified that is monotonically increasing with ̅ and ̅̅ ̅̅ for any . Thus, the optimal solution of (17) corresponds to the maximum values of ̅ and ̅̅ ̅̅ .
The maximum ̅ and ̅̅ ̅̅ can be obtained by identifying the worst-case ∈ and worst-case ∈ Σ . Due to the special structures of the and Σ in (8)-(9), the worstcase and have closed-form solutions: where the superscript means 'the worst-case'.
Step 2: identify the worst-case and the corresponding It is easy to verify that is convex in , so (17) can be rewritten as (20) according to the max-min theorem [30]: To simplify expressions, let ( ) be the function inside the square bracket in (20), and so (20) is equal to Due to the convexity of in , it is easy to proof that ( ) is also convex in . Also, problem (21) involves only one decision variable , and ( ) could be non-smooth. These observations motivate us to design a bisection search method to solve (21): where the structure of allows us to use the sorting method in [31] (i.e., SortM) to identify the ∈ that corresponds to the maximum ∑ ( )| ̅̅ ̅̅ , ̅̅̅̅ * =1 (i.e., ).

B. Iterative Solution Algorithm for Problem (P2)
With the cutting-plane-based constraint (23), we can develop an iterative solution algorithm to gradually enforce the DR-CVaR constraints in (14) and solve the problem (P2) in (15).
We also have one enhancement to accelerate Algorithm 2. The upward-reserve DR-CVaR constraints are rewritten as: We can first calculate the worst-case CVaR of − ̃ through the Steps 1-3 in Section IV-A (the result is denoted as Θ ) and use (25) to replace the original constraints in (24).
Similarly, after calculating the worst-case CVaR of ̃ evaluated at 1 − (the result is Θ ), we can use (26) to replace the downward-reserve DR-CVaR constraints in (P2).
Now, the reserve-related DR-CVaR constraints in (P2) are replaced by linear constraints (25)- (26), and we only need to apply the cutting-planes to branch-flow DR-CVaR constraints.

C. Remarks on the Solution Methodology
1) The whole methodology does not rely on any sampling procedure. First, for the evaluation of the worst-case CVaR, the Steps 1-3 in Section IV-A can be quickly implemented since they do not require any sample-based simulations but only involve the bisection search and some analytical solutions. Second, for the solving of the proposed OPF model, cuttingplanes are used to enforce the DR-CVaR constraints. Compared to the traditional methods with sample-based approximation of CVaR (e.g., [24]), the cutting-plane method is computationally more desirable since it makes the scale of the proposed OPF model irrelevant to the size of historical samples.

V. CASE STUDIES
The proposed OPF model is tested on a modified IEEE 2736bus system, with the original data from [32]. The system has 2736 buses, 289 controllable units, 3504 transmission lines and 10 wind farms. The total capacities of the controllable units and wind farms are 28.88 GW and 3.08 GW, respectively. The flow limits of all transmission lines are reduced by 20% to create transfer congestions. All numerical tests are carried out on a desktop with Intel Core i7-8700 CPU and 16GB RAM.

A. Construction of the Proposed Ambiguity set
In this paper, the underlying marginal distributions of the forecast errors at all wind farms and the associated correlations are extracted from Wind Integration National Dataset Toolkit [33]. Based on the marginal distributions and correlations, the Nataf method in [34] is used to generate some forecast error samples, which forms the historical dataset available to system operators. With the samples as the input, the bootstrap procedure in Section III-C is used to construct the proposed ambiguity set. The number of the bootstrap resamples is 2000. The confidence level in (10)-(12) is 95%. As an example, 4000 forecast error samples are generated as the input. The number of the Gaussian components given by the bootstrap method is six. For the 4 th Gaussian component, the associated bootstrap means of the forecast errors at wind farms 1, 2 and 5 are plotted in Fig. 1. It is seen that there is a correlation between the means for wind farm 1 and 2. This is speculated to be related to the high dependency between the forecast errors themselves (the correlation coefficient of the two forecast errors is around 0.75). As a comparison, the means for wind farms 1 and 5 have no obvious dependency since the corresponding forecast errors are weakly correlated. The potential correlation between the means is another reason why we use ellipsoids as the credible regions of the means (the first reason can be seen in Section IV-C), considering ellipsoids have a good capability of modelling correlated factors [35].

B. Comparison of Different CVaR-constrained OPF Models
The proposed OPF model, denoted by DG-OPF, is compared with three other OPF models with CVaR constraints: i) NA-OPF: the GMM with fixed parameters is used and no distributional ambiguity is considered. Here, the GMM fitted at the first step of the bootstrap procedure is used.
ii) M-OPF: The traditional moment-based ambiguity set is used, which utilizes only the estimated mean ̂ and covariance ∑ of the forecast error (see its structure in [3], [4]). In this case, the worst-case CVaR of ̃ has a simple form [21]: iii) U-OPF: This model applies not only the mean and covariance information but also the unimodality knowledge to construct its ambiguity set. The model can be formulated as a semi-definite programming problem (see more details in [22]).
DG-OPF, NA-OPF and M-OPF are solved by GUROBI while U-OPF is solved by MOSEK.
We consider two sample sizes: 200 and 4000. The risk levels for reserve and branch flow-CVaR constraints are set to 2% and 4%, respectively. For each sample size, we do the following procedure ten times: a) generate samples; b) solve all the four models; c) re-generate 10 6 forecast error scenarios using the Nataf method; d) for each model, keep the optimal operating strategy fixed, and do the out-of-sample test (to check if all CVaR constraints are satisfied under the 10 6 scenarios).
For each OPF model, the average, maximum and minimum cost results are listed in Table I. Meanwhile, for each model, the worst-case violations of various CVaR constraints observed in the out-of-sample tests are reported in Table II. It is seen from Table I and II that NA-OPF has the lowest cost but fails to meet the reliability requirements (the violated CVaR constraints are marked by boldface type). The reason is that NA-OPF makes decisions based on a deterministic GMM distribution, which is estimated from a small dataset with limited distributional information. This GMM could deviate from the underlying true distribution, and the NA-OPF dispatch may become unsafe. In contrast to NA-OPF, each of the other three models accounts for not a single but a family of possible true distributions. This needs additional cost but can hedge against the ambiguity in the true distribution and yield more robust operational strategies.
For DG, M and U-OPF models, they all satisfy the reliability requirements, but this is at the expense of different costs. As shown in Table I, M-OPF needs the highest total cost and the inclusion of the unimodality information in U-OPF decreases the required cost. Compared to M-and U-OPF, DG-OPF needs much lower cost. To better present this point, the additional costs of the three models used to hedge against the distributional ambiguity (i.e., the additional costs in relative to the total cost of NA-OPF) are listed in Table III. Under the data size 200, M and U-OPF need additional 7.97% and 6.74% costs in average, respectively, while such percentage is reduced to 2.21% in DG-OPF. Similar results are observed for the data size 4000. * a positive value indicates the associated CVaR constraint is violated and it is the magnitude of the violation; a negative value indicates the CVaR constraint is not violated and it is the safety gap between the CVaR and operational limit. The difference in the conservativeness of DG, M and U-OPF can be also observed in Table IV, which lists the amounts of the upward reserves procured by these models. In average, the amount for M-OPF is more than 70% higher than that for DG-OPF. Though the amount for M-OPF guarantees no violations of the upward reserve-related CVaR constraints, it could be unnecessarily high. We also compare the worst-case distributions in M, U and DG-OPF models. For the upward reserve CVaR constraints, the associated worst-case distributions of the total wind power forecast error are plotted in Fig. 2. The empirical histogram of the total forecast error formed by the 10 6 scenarios for out-ofsample-tests and the upward reserves procured by the three models are also plotted. The worst-case distribution of M-OPF is a two-supported discrete distribution, and the left support has a total probability mass of 2% and is far from the zero point, which causes the procurement of the large amount of upward reserve; In the U-OPF with the unimodality information, the worst-case distribution consists of two uniform distributions, and the left one has a total probability mass of 2.89%, and the occurrence probability of very large downward forecast errors is reduced; For DG-OPF, its worst-case distribution is still a continuous GMM and it has a shorter stretch at the left-hand tail compared to the cases in M and U-OPF, and thus the occurrence probability of very large downward forecast errors is the lowest.   Table V lists the solving time of DG, M, and U-OPF. As can be seen, the solving time of DG-OPF is longer than that of M-OPF but shorter than that of U-OPF. The reasons include that DG-OPF has more complex reformulations of DR-CVaR PDF constraints compared to M-OPF, and DG-OPF involves only linear programs while U-OPF is a semi-definite model. Still, DG-OPF can be solved within seconds on this 2736-bus system. This proves the usefulness of DG-OPF in short-term dispatch (e.g., 5 minutes ahead) on large networks. The solving time of DG-OPF also changes slightly as the data size increases. This is expected since the scale of DG-OPF is irrelated to the data size. Moreover, DG-OPF uses averagely only 37 ms to evaluate the worst-case CVaRs for totally 7008 branch flow constraints. This shows the scalability of the Steps 1-3 in Section IV-A. In summary, compared to M and U-OPF, DG-OPF uses more economical dispatch to satisfy the same reliability requirements while still maintaining a high computational efficiency.

VI. CONCLUSION
This paper proposes a GMM-based distributionally robust OPF model with CVaR constraints. A new ambiguity set based on uncertain GMM parameters is constructed, and a scalable cutting-plane-based solution methodology is designed. Case studies on a modified 2736-bus system show that the proposed GMM-based ambiguity set is capable of effectively hedging against the distributional ambiguity of wind power forecast error while offering less conservative dispatch than the moment and unimodality-based ambiguity sets. The case studies also reveal the scalability of the proposed OPF model, validating its potential application for short-term dispatch on large networks.